diff --git a/src/asset.rs b/src/asset.rs index 6d40b9e8d..90b5ba570 100644 --- a/src/asset.rs +++ b/src/asset.rs @@ -4,14 +4,15 @@ use crate::commodity::CommodityID; use crate::process::{FlowDirection, Process, ProcessFlow, ProcessID, ProcessParameter}; use crate::region::RegionID; use crate::simulation::CommodityPrices; -use crate::time_slice::TimeSliceID; -use crate::units::{Activity, ActivityPerCapacity, Capacity, Dimensionless, MoneyPerActivity}; +use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; +use crate::units::{ + Activity, ActivityPerCapacityPerYear, ActivityPerYear, Capacity, MoneyPerActivity, PerYear, +}; use anyhow::{Context, Result, ensure}; use indexmap::IndexMap; use itertools::{Itertools, chain}; use log::{debug, warn}; use serde::{Deserialize, Serialize}; -use std::collections::HashMap; use std::hash::{Hash, Hasher}; use std::ops::{Deref, RangeInclusive}; use std::rc::Rc; @@ -84,7 +85,7 @@ pub struct Asset { /// The [`Process`] that this asset corresponds to process: Rc, /// Activity limits for this asset - activity_limits: Rc>>, + activity_limits: Rc>>, /// The commodity flows for this asset flows: Rc>, /// The [`ProcessParameter`] corresponding to the asset's region and commission year @@ -275,8 +276,11 @@ impl Asset { } /// Get the activity limits for this asset in a particular time slice - pub fn get_activity_limits(&self, time_slice: &TimeSliceID) -> RangeInclusive { - let limits = &self.activity_limits[time_slice]; + pub fn get_activity_limits( + &self, + ts_selection: &TimeSliceSelection, + ) -> RangeInclusive { + let limits = &self.activity_limits[ts_selection]; let max_act = self.max_activity(); // limits in real units (which are user defined) @@ -287,12 +291,50 @@ impl Asset { pub fn get_activity_per_capacity_limits( &self, time_slice: &TimeSliceID, - ) -> RangeInclusive { - let limits = &self.activity_limits[time_slice]; + ) -> RangeInclusive { + // Get limits for this time slice specifically + let mut limits = self + .activity_limits + .get(&TimeSliceSelection::Single(time_slice.clone())) + .cloned() + .unwrap_or(PerYear(0.0)..=PerYear(1.0)); + let mut combine_limits = |ts_selection| { + if let Some(other_limits) = self.activity_limits.get(&ts_selection) { + let lower = limits.start().max(*other_limits.start()); + let upper = limits.end().min(*other_limits.end()); + limits = lower..=upper; + } + }; + + // Take into account possibly more restrictive limits at seasonal and annual level + combine_limits(TimeSliceSelection::Season(time_slice.season.clone())); + combine_limits(TimeSliceSelection::Annual); + let cap2act = self.process.capacity_to_activity; (cap2act * *limits.start())..=(cap2act * *limits.end()) } + /// Iterate over all activity limits for this asset + pub fn iter_activity_limits( + &self, + time_slice_info: &TimeSliceInfo, + ) -> impl Iterator)> { + let max_act = self.max_activity(); + self.activity_limits + .iter() + .map(move |(ts_selection, limits)| { + // NB: Unit types are correct here, but we're circumventing checks so we don't have + // to define a bunch of extra `Mul` impls + let mul = max_act.value() * time_slice_info.get_duration(ts_selection).value(); + let to_act = |limit: PerYear| Activity(mul * limit.value()); + + ( + ts_selection, + to_act(*limits.start())..=to_act(*limits.end()), + ) + }) + } + /// Get the operating cost for this asset in a given year and time slice pub fn get_operating_cost(&self, year: u32, time_slice: &TimeSliceID) -> MoneyPerActivity { // The cost for all commodity flows (including levies/incentives) @@ -916,7 +958,6 @@ mod tests { use itertools::{Itertools, assert_equal}; use map_macro::hash_map; use rstest::{fixture, rstest}; - use std::collections::HashMap; use std::iter; use std::rc::Rc; @@ -963,7 +1004,7 @@ mod tests { let flows = hash_map! {(region_id.clone(), 2020) => flow_map.into()}; // Create empty activity limits map - let activity_limits = hash_map! {(region_id.clone(), 2020) => Rc::new(HashMap::new())}; + let activity_limits = hash_map! {(region_id.clone(), 2020) => Rc::new(IndexMap::new())}; let process = Rc::new(Process { id: ProcessID::from("PROC1"), @@ -1054,7 +1095,7 @@ mod tests { .collect(); let activity_limits = years .iter() - .map(|&year| ((region_id.clone(), year), Rc::new(HashMap::new()))) + .map(|&year| ((region_id.clone(), year), Rc::new(IndexMap::new()))) .collect(); let flows = years .iter() @@ -1106,10 +1147,10 @@ mod tests { season: "winter".into(), time_of_day: "day".into(), }; - let fraction_limits = Dimensionless(1.0)..=Dimensionless(2.0); + let fraction_limits = PerYear(0.5)..=PerYear(1.0); let mut flows = ProcessFlowsMap::new(); let mut activity_limits = ProcessActivityLimitsMap::new(); - let limit_map = Rc::new(hash_map! {time_slice => fraction_limits}); + let limit_map = Rc::new(indexmap! {time_slice.into() => fraction_limits}); for year in [2010, 2020] { // empty flows map, but this is fine for our purposes flows.insert((region_id.clone(), year), Rc::new(IndexMap::new())); @@ -1143,8 +1184,8 @@ mod tests { #[rstest] fn test_asset_get_activity_limits(asset_with_activity_limits: Asset, time_slice: TimeSliceID) { assert_eq!( - asset_with_activity_limits.get_activity_limits(&time_slice), - Activity(6.0)..=Activity(12.0) + asset_with_activity_limits.get_activity_limits(&time_slice.into()), + ActivityPerYear(3.0)..=ActivityPerYear(6.0) ); } @@ -1154,8 +1195,8 @@ mod tests { time_slice: TimeSliceID, ) { assert_eq!( - asset_with_activity_limits.get_activity_per_capacity_limits(&time_slice), - ActivityPerCapacity(3.0)..=ActivityPerCapacity(6.0) + asset_with_activity_limits.get_activity_per_capacity_limits(&time_slice.into()), + ActivityPerCapacityPerYear(1.5)..=ActivityPerCapacityPerYear(3.0) ); } diff --git a/src/fixture.rs b/src/fixture.rs index 051957eaa..53d7c1048 100644 --- a/src/fixture.rs +++ b/src/fixture.rs @@ -152,7 +152,7 @@ pub fn process( // Create maps with (empty) entries for every region/year combo let activity_limits = iproduct!(region_ids.iter(), years.iter()) - .map(|(region_id, year)| ((region_id.clone(), *year), Rc::new(HashMap::new()))) + .map(|(region_id, year)| ((region_id.clone(), *year), Rc::new(IndexMap::new()))) .collect(); let flows = iproduct!(region_ids.iter(), years.iter()) .map(|(region_id, year)| ((region_id.clone(), *year), Rc::new(IndexMap::new()))) diff --git a/src/graph/validate.rs b/src/graph/validate.rs index bf5798349..6913d3d37 100644 --- a/src/graph/validate.rs +++ b/src/graph/validate.rs @@ -4,11 +4,54 @@ use crate::commodity::{CommodityMap, CommodityType}; use crate::process::ProcessMap; use crate::region::RegionID; use crate::time_slice::{TimeSliceInfo, TimeSliceLevel, TimeSliceSelection}; -use crate::units::{Dimensionless, Flow}; +use crate::units::{Flow, PerYear}; use anyhow::{Context, Result, ensure}; use indexmap::IndexMap; +use std::ops::RangeInclusive; use strum::IntoEnumIterator; +/// Check whether any availability is permissible for the [`TimeSliceSelection`] and its parents +fn has_availability_with_parents( + limits: &IndexMap>, + ts_selection: &TimeSliceSelection, +) -> bool { + limits + .get(ts_selection) + .is_none_or(|limits| *limits.end() > PerYear(0.0)) + && ts_selection + .get_parent() + .is_none_or(|ts_selection| has_availability_with_parents(limits, &ts_selection)) +} + +/// Check whether any availability is permissible for any children of the specified +/// [`TimeSliceSelection`]. +/// +/// Returns true if there are no children. +fn children_have_availability( + limits: &IndexMap>, + ts_selection: &TimeSliceSelection, + time_slice_info: &TimeSliceInfo, +) -> bool { + ts_selection.iter_children(time_slice_info).all(|child| { + !(limits + .get(ts_selection) + .is_none_or(|limits| *limits.end() > PerYear(0.0)) + && children_have_availability(limits, &child, time_slice_info)) + }) +} + +/// Check whether there is any availability for this [`TimeSliceSelection`]. +/// +/// Searches at both coarser and finer [`TimeSliceLevel`]s to check for limits. +fn has_any_availability( + limits: &IndexMap>, + ts_selection: &TimeSliceSelection, + time_slice_info: &TimeSliceInfo, +) -> bool { + has_availability_with_parents(limits, ts_selection) + && children_have_availability(limits, ts_selection, time_slice_info) +} + /// Prepares a graph for validation with [`validate_commodities_graph`]. /// /// It takes a base graph produced by `create_commodities_graph_for_region_year`, and modifies it to @@ -41,17 +84,10 @@ fn prepare_commodities_graph_for_validation( }; let process = &processes[process_id]; - // Check if the process has availability > 0 in any time slice in the selection - time_slice_selection - .iter(time_slice_info) - .any(|(time_slice, _)| { - let Some(limits_map) = process.activity_limits.get(&key) else { - return false; - }; - limits_map - .get(time_slice) - .is_some_and(|avail| *avail.end() > Dimensionless(0.0)) - }) + // Check whether there is any availability for this time slice selection + process.activity_limits.get(&key).is_some_and(|limits| { + has_any_availability(limits, time_slice_selection, time_slice_info) + }) }); // Add demand edges diff --git a/src/input.rs b/src/input.rs index e82ad6c47..08358facc 100644 --- a/src/input.rs +++ b/src/input.rs @@ -31,6 +31,24 @@ use region::read_regions; mod time_slice; use time_slice::read_time_slice_info; +/// A trait which provides a method to insert a key and value into a map +pub trait Insert { + /// Insert a key and value into the map + fn insert(&mut self, key: K, value: V) -> Option; +} + +impl Insert for HashMap { + fn insert(&mut self, key: K, value: V) -> Option { + HashMap::insert(self, key, value) + } +} + +impl Insert for IndexMap { + fn insert(&mut self, key: K, value: V) -> Option { + IndexMap::insert(self, key, value) + } +} + /// Read a series of type `T`s from a CSV file. /// /// Will raise an error if the file is empty. @@ -163,9 +181,10 @@ where /// Inserts a key-value pair into a `HashMap` if the key does not already exist. /// /// If the key already exists, it returns an error with a message indicating the key's existence. -pub fn try_insert(map: &mut HashMap, key: &K, value: V) -> Result<()> +pub fn try_insert(map: &mut M, key: &K, value: V) -> Result<()> where - K: Eq + Hash + Clone + fmt::Debug, + M: Insert, + K: Eq + Hash + Clone + std::fmt::Debug, { let existing = map.insert(key.clone(), value).is_some(); ensure!(!existing, "Key {key:?} already exists in the map"); diff --git a/src/input/process.rs b/src/input/process.rs index 3352e7d8d..4b552aec4 100644 --- a/src/input/process.rs +++ b/src/input/process.rs @@ -59,8 +59,7 @@ pub fn read_processes( ) -> Result { let base_year = milestone_years[0]; let mut processes = read_processes_file(model_dir, milestone_years, region_ids, commodities)?; - let mut activity_limits = - read_process_availabilities(model_dir, &processes, time_slice_info, base_year)?; + let mut activity_limits = read_process_availabilities(model_dir, &processes, time_slice_info)?; let mut flows = read_process_flows(model_dir, &mut processes, commodities)?; let mut parameters = read_process_parameters(model_dir, &processes, base_year)?; diff --git a/src/input/process/availability.rs b/src/input/process/availability.rs index 61fe44002..8b9fda103 100644 --- a/src/input/process/availability.rs +++ b/src/input/process/availability.rs @@ -1,9 +1,9 @@ //! Code for reading process availabilities CSV file -use super::super::{format_items_with_cap, input_err_msg, read_csv, try_insert}; -use crate::process::{Process, ProcessActivityLimitsMap, ProcessID, ProcessMap}; +use super::super::{input_err_msg, read_csv, try_insert}; +use crate::process::{ProcessActivityLimitsMap, ProcessID, ProcessMap}; use crate::region::parse_region_str; -use crate::time_slice::TimeSliceInfo; -use crate::units::{Dimensionless, Year}; +use crate::time_slice::{TimeSliceInfo, TimeSliceSelection}; +use crate::units::PerYear; use crate::year::parse_year_str; use anyhow::{Context, Result, ensure}; use itertools::iproduct; @@ -24,32 +24,27 @@ struct ProcessAvailabilityRaw { commission_years: String, time_slice: String, limit_type: LimitType, - value: Dimensionless, + value: PerYear, } impl ProcessAvailabilityRaw { fn validate(&self) -> Result<()> { // Check availability value ensure!( - self.value >= Dimensionless(0.0) && self.value <= Dimensionless(1.0), + self.value >= PerYear(0.0) && self.value <= PerYear(1.0), "Value for availability must be between 0 and 1 inclusive" ); Ok(()) } - /// Calculate fraction of annual energy as availability multiplied by time slice length. - /// - /// The resulting limits are max/min energy produced/consumed in each time slice per - /// `capacity_to_activity` units of capacity. - fn to_bounds(&self, ts_length: Year) -> RangeInclusive { + /// Get this limit as a range + fn to_range(&self) -> RangeInclusive { // We know ts_length also represents a fraction of a year, so this is ok. - let ts_frac = ts_length / Year(1.0); - let value = self.value * ts_frac; match self.limit_type { - LimitType::LowerBound => value..=ts_frac, - LimitType::UpperBound => Dimensionless(0.0)..=value, - LimitType::Equality => value..=value, + LimitType::LowerBound => self.value..=PerYear(1.0), + LimitType::UpperBound => PerYear(0.0)..=self.value, + LimitType::Equality => self.value..=self.value, } } } @@ -75,7 +70,6 @@ enum LimitType { /// * `model_dir` - Folder containing model configuration files /// * `processes` - Map of processes /// * `time_slice_info` - Information about seasons and times of day -/// * `base_year` - First milestone year of simulation /// /// # Returns /// @@ -85,17 +79,11 @@ pub fn read_process_availabilities( model_dir: &Path, processes: &ProcessMap, time_slice_info: &TimeSliceInfo, - base_year: u32, ) -> Result> { let file_path = model_dir.join(PROCESS_AVAILABILITIES_FILE_NAME); let process_availabilities_csv = read_csv(&file_path)?; - read_process_availabilities_from_iter( - process_availabilities_csv, - processes, - time_slice_info, - base_year, - ) - .with_context(|| input_err_msg(&file_path)) + read_process_availabilities_from_iter(process_availabilities_csv, processes, time_slice_info) + .with_context(|| input_err_msg(&file_path)) } /// Process raw process availabilities input data into [`ProcessActivityLimitsMap`]s. @@ -105,7 +93,6 @@ pub fn read_process_availabilities( /// * `iter` - Iterator of raw process availability records /// * `processes` - Map of processes /// * `time_slice_info` - Information about seasons and times of day -/// * `base_year` - First milestone year of simulation /// /// # Returns /// @@ -115,7 +102,6 @@ fn read_process_availabilities_from_iter( iter: I, processes: &ProcessMap, time_slice_info: &TimeSliceInfo, - base_year: u32, ) -> Result> where I: Iterator, @@ -147,108 +133,45 @@ where let ts_selection = time_slice_info.get_selection(&record.time_slice)?; // Insert the activity limit into the map - let limits_map = map - .entry(id.clone()) - .or_insert_with(ProcessActivityLimitsMap::new); + let limits_map: &mut ProcessActivityLimitsMap = map.entry(id.clone()).or_default(); for (region_id, year) in iproduct!(&record_regions, &record_years) { - let limits_map_inner = limits_map - .entry((region_id.clone(), *year)) - .or_insert_with(|| Rc::new(HashMap::new())); + let limits_map_inner = limits_map.entry((region_id.clone(), *year)).or_default(); let limits_map_inner = Rc::get_mut(limits_map_inner).unwrap(); - for (time_slice, ts_length) in ts_selection.iter(time_slice_info) { - let bounds = record.to_bounds(ts_length); - try_insert(limits_map_inner, time_slice, bounds.clone())?; - } + try_insert(limits_map_inner, &ts_selection, record.to_range())?; } } - validate_activity_limits_maps(&map, processes, time_slice_info, base_year)?; + add_fallback_limits(&mut map, processes); Ok(map) } -/// Check that the activity limits cover every time slice and all regions/years of the process -fn validate_activity_limits_maps( - all_availabilities: &HashMap, +/// Add limits of 0..=1 for whole year for every region and year, everywhere needed. +/// +/// The reason this is needed is because users are not required to provide limits for every single +/// time slice (in fact, there may be none at all), and we need to add a constraint to the +/// optimisation so that the total activity for the year does not exceed the maximum for a given +/// asset. +fn add_fallback_limits( + all_availabilities: &mut HashMap, processes: &ProcessMap, - time_slice_info: &TimeSliceInfo, - base_year: u32, -) -> Result<()> { +) { for (process_id, process) in processes { // A map of maps: the outer map is keyed by region and year; the inner one by time slice - let map_for_process = all_availabilities - .get(process_id) - .with_context(|| format!("Missing availabilities for process {process_id}"))?; - - check_missing_milestone_years(process, map_for_process, base_year)?; - check_missing_time_slices(process, map_for_process, time_slice_info)?; - } - - Ok(()) -} - -/// Check every milestone year in which the process can be commissioned has availabilities. -/// -/// Entries for non-milestone years in which the process can be commissioned (which are only -/// required for pre-defined assets, if at all) are not required and will be checked lazily when -/// assets requiring them are constructed. -fn check_missing_milestone_years( - process: &Process, - map_for_process: &ProcessActivityLimitsMap, - base_year: u32, -) -> Result<()> { - let process_milestone_years = process - .years - .iter() - .copied() - .filter(|&year| year >= base_year); - let mut missing = Vec::new(); - for (region_id, year) in iproduct!(&process.regions, process_milestone_years) { - if !map_for_process.contains_key(&(region_id.clone(), year)) { - missing.push((region_id, year)); + // selection + let map_for_process = all_availabilities.entry(process_id.clone()).or_default(); + + // Iterate over all regions and years for this process and insert a fallback availability + // limit covering the whole year, in case it's needed + for (region, &year) in iproduct!(&process.regions, &process.years) { + // Should succeed because there should be only one reference + let map_for_region_year = + Rc::get_mut(map_for_process.entry((region.clone(), year)).or_default()).unwrap(); + map_for_region_year + .entry(TimeSliceSelection::Annual) + .or_insert(PerYear(0.0)..=PerYear(1.0)); } } - - ensure!( - missing.is_empty(), - "Process {} is missing availabilities for the following regions and milestone years: {}", - &process.id, - format_items_with_cap(&missing) - ); - - Ok(()) -} - -/// Check that entries for all time slices are provided for any process/region/year combo for which -/// we have any entries at all -fn check_missing_time_slices( - process: &Process, - map_for_process: &ProcessActivityLimitsMap, - time_slice_info: &TimeSliceInfo, -) -> Result<()> { - let mut missing = Vec::new(); - for (region_id, &year) in iproduct!(&process.regions, &process.years) { - if let Some(map_for_region_year) = map_for_process.get(&(region_id.clone(), year)) { - // There are at least some entries for this region/year combo; check if there are - // any time slices not covered - missing.extend( - time_slice_info - .iter_ids() - .filter(|ts| !map_for_region_year.contains_key(ts)) - .map(|ts| (region_id, year, ts)), - ); - } - } - - ensure!( - missing.is_empty(), - "Availabilities supplied for some, but not all time slices, for process {}. The following \ - regions, years and time slices are missing: {}", - &process.id, - format_items_with_cap(&missing) - ); - - Ok(()) } #[cfg(test)] @@ -257,7 +180,7 @@ mod tests { fn create_process_availability_raw( limit_type: LimitType, - value: Dimensionless, + value: PerYear, ) -> ProcessAvailabilityRaw { ProcessAvailabilityRaw { process_id: "process".into(), @@ -272,56 +195,51 @@ mod tests { #[test] fn test_validate() { // Valid - let valid = create_process_availability_raw(LimitType::LowerBound, Dimensionless(0.5)); + let valid = create_process_availability_raw(LimitType::LowerBound, PerYear(0.5)); assert!(valid.validate().is_ok()); - let valid = create_process_availability_raw(LimitType::LowerBound, Dimensionless(0.0)); + let valid = create_process_availability_raw(LimitType::LowerBound, PerYear(0.0)); assert!(valid.validate().is_ok()); - let valid = create_process_availability_raw(LimitType::LowerBound, Dimensionless(1.0)); + let valid = create_process_availability_raw(LimitType::LowerBound, PerYear(1.0)); assert!(valid.validate().is_ok()); // Invalid: negative value - let invalid = create_process_availability_raw(LimitType::LowerBound, Dimensionless(-0.5)); + let invalid = create_process_availability_raw(LimitType::LowerBound, PerYear(-0.5)); assert!(invalid.validate().is_err()); // Invalid: value greater than 1 - let invalid = create_process_availability_raw(LimitType::LowerBound, Dimensionless(1.5)); + let invalid = create_process_availability_raw(LimitType::LowerBound, PerYear(1.5)); assert!(invalid.validate().is_err()); // Invalid: infinity value let invalid = - create_process_availability_raw(LimitType::LowerBound, Dimensionless(f64::INFINITY)); + create_process_availability_raw(LimitType::LowerBound, PerYear(f64::INFINITY)); assert!(invalid.validate().is_err()); // Invalid: negative infinity value - let invalid = create_process_availability_raw( - LimitType::LowerBound, - Dimensionless(f64::NEG_INFINITY), - ); + let invalid = + create_process_availability_raw(LimitType::LowerBound, PerYear(f64::NEG_INFINITY)); assert!(invalid.validate().is_err()); // Invalid: NaN value - let invalid = - create_process_availability_raw(LimitType::LowerBound, Dimensionless(f64::NAN)); + let invalid = create_process_availability_raw(LimitType::LowerBound, PerYear(f64::NAN)); assert!(invalid.validate().is_err()); } #[test] fn test_to_bounds() { - let ts_length = Year(0.1); - // Lower bound - let raw = create_process_availability_raw(LimitType::LowerBound, Dimensionless(0.5)); - let bounds = raw.to_bounds(ts_length); - assert_eq!(bounds, Dimensionless(0.05)..=Dimensionless(0.1)); + let raw = create_process_availability_raw(LimitType::LowerBound, PerYear(0.5)); + let bounds = raw.to_range(); + assert_eq!(bounds, PerYear(0.5)..=PerYear(1.0)); // Upper bound - let raw = create_process_availability_raw(LimitType::UpperBound, Dimensionless(0.5)); - let bounds = raw.to_bounds(ts_length); - assert_eq!(bounds, Dimensionless(0.0)..=Dimensionless(0.05)); + let raw = create_process_availability_raw(LimitType::UpperBound, PerYear(0.5)); + let bounds = raw.to_range(); + assert_eq!(bounds, PerYear(0.0)..=PerYear(0.5)); // Equality - let raw = create_process_availability_raw(LimitType::Equality, Dimensionless(0.5)); - let bounds = raw.to_bounds(ts_length); - assert_eq!(bounds, Dimensionless(0.05)..=Dimensionless(0.05)); + let raw = create_process_availability_raw(LimitType::Equality, PerYear(0.5)); + let bounds = raw.to_range(); + assert_eq!(bounds, PerYear(0.5)..=PerYear(0.5)); } } diff --git a/src/process.rs b/src/process.rs index 97c183eeb..0bec3e62e 100644 --- a/src/process.rs +++ b/src/process.rs @@ -3,10 +3,10 @@ use crate::commodity::{Commodity, CommodityID}; use crate::id::define_id_type; use crate::region::RegionID; -use crate::time_slice::TimeSliceID; +use crate::time_slice::{TimeSliceID, TimeSliceSelection}; use crate::units::{ ActivityPerCapacity, Dimensionless, FlowPerActivity, MoneyPerActivity, MoneyPerCapacity, - MoneyPerCapacityPerYear, MoneyPerFlow, + MoneyPerCapacityPerYear, MoneyPerFlow, PerYear, }; use indexmap::{IndexMap, IndexSet}; use serde_string_enum::DeserializeLabeledStringEnum; @@ -24,7 +24,7 @@ pub type ProcessMap = IndexMap>; /// The value is calculated as availability multiplied by time slice length. The limits are given as /// ranges, depending on the user-specified limit type and value for availability. pub type ProcessActivityLimitsMap = - HashMap<(RegionID, u32), Rc>>>; + HashMap<(RegionID, u32), Rc>>>; /// A map of [`ProcessParameter`]s, keyed by region and year pub type ProcessParameterMap = HashMap<(RegionID, u32), Rc>; diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index b10b0c759..cd64b33a9 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -417,9 +417,9 @@ fn get_demand_limiting_capacity( // Calculate max capacity required for this time slice selection // For commodities with a coarse time slice level, we have to allow the possibility that all // of the demand gets served by production in a single time slice - for (time_slice, _) in time_slice_selection.iter(time_slice_info) { + for (time_slice, duration) in time_slice_selection.iter(time_slice_info) { let max_flow_per_cap = - *asset.get_activity_per_capacity_limits(time_slice).end() * coeff; + *asset.get_activity_per_capacity_limits(time_slice).end() * duration * coeff; if max_flow_per_cap != FlowPerCapacity(0.0) { capacity = capacity.max(demand_for_selection / max_flow_per_cap); } @@ -652,11 +652,10 @@ mod tests { }; use crate::process::{FlowType, ProcessFlow}; use crate::region::RegionID; - use crate::time_slice::{TimeSliceID, TimeSliceInfo}; - use crate::units::{Dimensionless, Flow, FlowPerActivity, MoneyPerFlow}; + use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; + use crate::units::{Flow, FlowPerActivity, MoneyPerFlow, PerYear}; use indexmap::indexmap; use itertools::Itertools; - use map_macro::hash_map; use rstest::rstest; use std::rc::Rc; @@ -695,7 +694,7 @@ mod tests { // Add activity limits process.activity_limits.insert( (region_id.clone(), 2015), - Rc::new(hash_map! {time_slice.clone() => Dimensionless(0.0)..=Dimensionless(1.0)}), + Rc::new(indexmap! {TimeSliceSelection::Single(time_slice.clone()) => PerYear(0.0)..=PerYear(1.0)}), ); // Create asset with the configured process @@ -749,11 +748,11 @@ mod tests { ); // Add activity limits for both time slices with different limits - let limits = hash_map! { + let limits = indexmap! { // Higher limit for day - time_slice1.clone() => Dimensionless(0.0)..=Dimensionless(2.0), + TimeSliceSelection::Single(time_slice1.clone()) => PerYear(0.0)..=PerYear(1.0), // Zero limit for night - should be skipped - time_slice2.clone() => Dimensionless(0.0)..=Dimensionless(0.0) + TimeSliceSelection::Single(time_slice2.clone()) => PerYear(0.0)..=PerYear(0.0) }; process .activity_limits @@ -773,9 +772,6 @@ mod tests { get_demand_limiting_capacity(&time_slice_info2, &asset, &commodity_rc, &demand); // Expected: maximum of the capacity requirements across time slices (excluding zero limit) - // Time slice 1: demand (4.0) / (activity_limit (2.0) * coeff (1.0)) = 2.0 - // Time slice 2: skipped due to zero activity limit - // Maximum = 2.0 - assert_eq!(result, Capacity(2.0)); + assert_eq!(result, Capacity(8.0)); } } diff --git a/src/simulation/investment/appraisal/constraints.rs b/src/simulation/investment/appraisal/constraints.rs index 36473a232..e62004675 100644 --- a/src/simulation/investment/appraisal/constraints.rs +++ b/src/simulation/investment/appraisal/constraints.rs @@ -7,6 +7,7 @@ use crate::time_slice::{TimeSliceID, TimeSliceInfo}; use crate::units::{Capacity, Flow}; use highs::RowProblem as Problem; use indexmap::IndexMap; +use itertools::Itertools; /// Adds a capacity constraint to the problem. /// @@ -51,13 +52,20 @@ pub fn add_activity_constraints( asset: &AssetRef, capacity_var: Variable, activity_vars: &IndexMap, + time_slice_info: &TimeSliceInfo, ) { match asset.state() { AssetState::Commissioned { .. } => { - add_activity_constraints_for_existing(problem, asset, activity_vars); + add_activity_constraints_for_existing(problem, asset, activity_vars, time_slice_info); } AssetState::Candidate => { - add_activity_constraints_for_candidate(problem, asset, capacity_var, activity_vars); + add_activity_constraints_for_candidate( + problem, + asset, + capacity_var, + activity_vars, + time_slice_info, + ); } _ => panic!( "add_activity_constraints should only be called with Commissioned or Candidate assets" @@ -69,11 +77,18 @@ fn add_activity_constraints_for_existing( problem: &mut Problem, asset: &AssetRef, activity_vars: &IndexMap, + time_slice_info: &TimeSliceInfo, ) { - for (time_slice, var) in activity_vars { - let limits = asset.get_activity_limits(time_slice); + for (ts_selection, limits) in asset.iter_activity_limits(time_slice_info) { let limits = limits.start().value()..=limits.end().value(); - problem.add_row(limits, [(*var, 1.0)]); + let terms = ts_selection + .iter(time_slice_info) + .map(|(time_slice, _)| { + let var = *activity_vars.get(time_slice).unwrap(); + (var, 1.0) + }) + .collect_vec(); + problem.add_row(limits, &terms); } } @@ -82,17 +97,25 @@ fn add_activity_constraints_for_candidate( asset: &AssetRef, capacity_var: Variable, activity_vars: &IndexMap, + time_slice_info: &TimeSliceInfo, ) { - for (time_slice, activity_var) in activity_vars { - let limits = asset.get_activity_per_capacity_limits(time_slice); - let lower_limit = limits.start().value(); + for (ts_selection, limits) in asset.iter_activity_limits(time_slice_info) { let upper_limit = limits.end().value(); + let lower_limit = limits.start().value(); + + let mut terms_upper = vec![(capacity_var, -upper_limit)]; + let mut terms_lower = vec![(capacity_var, lower_limit)]; + for (time_slice, _) in ts_selection.iter(time_slice_info) { + let var = *activity_vars.get(time_slice).unwrap(); + terms_upper.push((var, 1.0)); + terms_lower.push((var, -1.0)); + } // Upper bound: activity ≤ capacity * upper_limit - problem.add_row(..=0.0, [(*activity_var, 1.0), (capacity_var, -upper_limit)]); + problem.add_row(..=0.0, &terms_upper); // Lower bound: activity ≥ capacity * lower_limit - problem.add_row(..=0.0, [(*activity_var, -1.0), (capacity_var, lower_limit)]); + problem.add_row(..=0.0, &terms_lower); } } diff --git a/src/simulation/investment/appraisal/optimisation.rs b/src/simulation/investment/appraisal/optimisation.rs index 6216a9808..319ca68d8 100644 --- a/src/simulation/investment/appraisal/optimisation.rs +++ b/src/simulation/investment/appraisal/optimisation.rs @@ -79,6 +79,7 @@ fn add_constraints( asset, variables.capacity_var, &variables.activity_vars, + time_slice_info, ); add_demand_constraints( problem, diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 4e2a3b5c5..8d438d306 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -155,13 +155,6 @@ impl VariableMap { .expect("No unmet demand variable for given params") } - /// Iterate over the asset variables - fn iter_asset_vars(&self) -> impl Iterator { - self.asset_vars - .iter() - .map(|((asset, time_slice), var)| (asset, time_slice, *var)) - } - /// Iterate over the keys for asset variables fn asset_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> { self.asset_vars.keys() @@ -254,7 +247,11 @@ impl Solution<'_> { self.constraint_keys .activity_keys .zip_duals(self.solution.dual_rows()) - .map(|((asset, time_slice), dual)| (asset, time_slice, dual)) + .flat_map(|((asset, ts_selection), dual)| { + ts_selection + .iter(self.time_slice_info) + .map(move |(ts, _)| (asset, ts, dual)) + }) } /// Keys and values for column duals. diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs index d0a229009..f16e931c5 100644 --- a/src/simulation/optimisation/constraints.rs +++ b/src/simulation/optimisation/constraints.rs @@ -4,9 +4,10 @@ use crate::asset::{AssetIterator, AssetRef}; use crate::commodity::{CommodityID, CommodityType}; use crate::model::Model; use crate::region::RegionID; -use crate::time_slice::{TimeSliceID, TimeSliceSelection}; +use crate::time_slice::{TimeSliceInfo, TimeSliceSelection}; use crate::units::UnitType; use highs::RowProblem as Problem; +use itertools::Itertools; /// Corresponding variables for a constraint along with the row offset in the solution pub struct KeysWithOffset { @@ -35,7 +36,7 @@ impl KeysWithOffset { pub type CommodityBalanceKeys = KeysWithOffset<(CommodityID, RegionID, TimeSliceSelection)>; /// Indicates the asset ID and time slice covered by each activity constraint -pub type ActivityKeys = KeysWithOffset<(AssetRef, TimeSliceID)>; +pub type ActivityKeys = KeysWithOffset<(AssetRef, TimeSliceSelection)>; /// The keys for different constraints pub struct ConstraintKeys { @@ -76,7 +77,8 @@ where let commodity_balance_keys = add_commodity_balance_constraints(problem, variables, model, assets, markets, year); - let activity_keys = add_activity_constraints(problem, variables); + let activity_keys = + add_activity_constraints(problem, variables, &model.time_slice_info, assets.clone()); // Return constraint keys ConstraintKeys { @@ -180,17 +182,29 @@ where /// See description in [the dispatch optimisation documentation][1]. /// /// [1]: https://energysystemsmodellinglab.github.io/MUSE2/model/dispatch_optimisation.html#a4-constraints-capacity--availability-for-standard-assets--a-in-mathbfastd- -fn add_activity_constraints(problem: &mut Problem, variables: &VariableMap) -> ActivityKeys { +fn add_activity_constraints<'a, I>( + problem: &mut Problem, + variables: &VariableMap, + time_slice_info: &TimeSliceInfo, + assets: I, +) -> ActivityKeys +where + I: Iterator + 'a, +{ // Row offset in problem. This line **must** come before we add more constraints. let offset = problem.num_rows(); let mut keys = Vec::new(); - for (asset, time_slice, var) in variables.iter_asset_vars() { - let limits = asset.get_activity_limits(time_slice); - let limits = limits.start().value()..=limits.end().value(); - - problem.add_row(limits, [(var, 1.0)]); - keys.push((asset.clone(), time_slice.clone())); + for asset in assets { + for (ts_selection, limits) in asset.iter_activity_limits(time_slice_info) { + let limits = limits.start().value()..=limits.end().value(); + let vars = ts_selection + .iter(time_slice_info) + .map(|(time_slice, _)| (variables.get_asset_var(asset, time_slice), 1.0)) + .collect_vec(); + problem.add_row(limits, &vars); + keys.push((asset.clone(), ts_selection.clone())); + } } ActivityKeys { offset, keys } diff --git a/src/time_slice.rs b/src/time_slice.rs index ee84c785b..1efedd92a 100644 --- a/src/time_slice.rs +++ b/src/time_slice.rs @@ -11,11 +11,24 @@ use serde::de::Error; use serde::{Deserialize, Serialize}; use serde_string_enum::DeserializeLabeledStringEnum; use std::fmt::Display; -use std::iter; +use std::{cmp, iter}; define_id_type! {Season} define_id_type! {TimeOfDay} +impl Season { + /// Compare this [`Season`] with another for sorting. + /// + /// Compares based on index in which it is stored in `time_slice_info` (ultimately: the order in + /// the input file). + pub fn compare(&self, other: &Season, time_slice_info: &TimeSliceInfo) -> cmp::Ordering { + // Get index of each - they will exist because all season are stored there + let idx1 = time_slice_info.seasons.get_index_of(self).unwrap(); + let idx2 = time_slice_info.seasons.get_index_of(other).unwrap(); + idx1.cmp(&idx2) + } +} + /// An ID describing season and time of day #[derive(Hash, Eq, PartialEq, Ord, PartialOrd, Clone, Debug)] pub struct TimeSliceID { @@ -25,6 +38,19 @@ pub struct TimeSliceID { pub time_of_day: TimeOfDay, } +impl TimeSliceID { + /// Compare this [`TimeSliceID`] with another for sorting. + /// + /// Compares based on index in which it is stored in `time_slice_info` (ultimately: the order in + /// the input file). + pub fn compare(&self, other: &TimeSliceID, time_slice_info: &TimeSliceInfo) -> cmp::Ordering { + // Get index of each - it will exist because all time slices are stored there + let idx1 = time_slice_info.time_slices.get_index_of(self).unwrap(); + let idx2 = time_slice_info.time_slices.get_index_of(other).unwrap(); + idx1.cmp(&idx2) + } +} + /// Only implement for tests as this is a bit of a footgun #[cfg(test)] impl From<&str> for TimeSliceID { @@ -94,6 +120,57 @@ impl TimeSliceSelection { } } + /// Compare this [`TimeSliceSelection`] with another for sorting. + /// + /// If the [`TimeSliceSelection`]s are both individual time slices or seasons, then they are + /// compared based on their order in `time_slice_info`, else they are compared based on + /// [`TimeSliceLevel`]. + pub fn compare(&self, other: &Self, time_slice_info: &TimeSliceInfo) -> cmp::Ordering { + match (self, other) { + (Self::Single(ts1), Self::Single(ts2)) => ts1.compare(ts2, time_slice_info), + (Self::Season(season1), Self::Season(season2)) => { + season1.compare(season2, time_slice_info) + } + (a, b) => a.level().cmp(&b.level()), + } + } + + /// Get the parent [`TimeSliceSelection`], if any. + /// + /// This means the [`TimeSliceSelection`] corresponding to a coarser [`TimeSliceLevel`] than + /// this one. + pub fn get_parent(&self) -> Option { + match self { + Self::Annual => None, + Self::Season(_) => Some(Self::Annual), + Self::Single(time_slice) => Some(Self::Season(time_slice.season.clone())), + } + } + + /// Iterate over this [`TimeSliceSelection`]'s children, if any. + /// + /// This means the [`TimeSliceSelection`]s corresponding to a fine [`TimeSliceLevel`] than + /// this one, if any. + pub fn iter_children<'a>( + &'a self, + time_slice_info: &'a TimeSliceInfo, + ) -> Box + 'a> { + let ts_info = time_slice_info; + match self { + TimeSliceSelection::Annual => { + Box::new(time_slice_info.seasons.keys().cloned().map(Self::Season)) + } + TimeSliceSelection::Season(season) => Box::new( + ts_info + .iter_ids() + .filter(move |ts| &ts.season == season) + .cloned() + .map(Self::Single), + ), + TimeSliceSelection::Single(_) => Box::new(iter::empty()), + } + } + /// Iterate over the subset of time slices in this selection pub fn iter<'a>( &'a self, @@ -192,7 +269,15 @@ impl Display for TimeSliceSelection { /// The time granularity for a particular operation #[derive( - PartialEq, PartialOrd, Copy, Clone, Debug, DeserializeLabeledStringEnum, strum::EnumIter, + PartialEq, + PartialOrd, + Eq, + Ord, + Copy, + Clone, + Debug, + DeserializeLabeledStringEnum, + strum::EnumIter, )] pub enum TimeSliceLevel { /// Treat individual time slices separately @@ -259,6 +344,15 @@ impl TimeSliceInfo { }) } + /// Get the duration covered by a `TimeSliceSelection` + pub fn get_duration(&self, ts_selection: &TimeSliceSelection) -> Year { + match ts_selection { + TimeSliceSelection::Annual => Year(1.0), + TimeSliceSelection::Season(season) => self.seasons[season], + TimeSliceSelection::Single(time_slice) => self.time_slices[time_slice], + } + } + /// Get a `TimeSliceSelection` from the specified string. /// /// If the string is empty, the default value is `TimeSliceSelection::Annual`. diff --git a/src/units.rs b/src/units.rs index c37a7a68b..cc0a4505d 100644 --- a/src/units.rs +++ b/src/units.rs @@ -244,12 +244,15 @@ unit_struct!(Capacity); unit_struct!(Year); // Derived quantities +unit_struct!(PerYear); unit_struct!(MoneyPerYear); unit_struct!(MoneyPerFlow); unit_struct!(MoneyPerCapacity); unit_struct!(MoneyPerCapacityPerYear); unit_struct!(MoneyPerActivity); +unit_struct!(ActivityPerYear); unit_struct!(ActivityPerCapacity); +unit_struct!(ActivityPerCapacityPerYear); unit_struct!(FlowPerActivity); unit_struct!(FlowPerCapacity); @@ -302,7 +305,41 @@ impl_div!(Money, Flow, MoneyPerFlow); impl_div!(Money, Capacity, MoneyPerCapacity); impl_div!(Money, Activity, MoneyPerActivity); impl_div!(Activity, Capacity, ActivityPerCapacity); +impl_div!(Activity, Year, ActivityPerYear); +impl_div!(ActivityPerCapacity, Year, ActivityPerCapacityPerYear); impl_div!(MoneyPerYear, Capacity, MoneyPerCapacityPerYear); impl_div!(MoneyPerActivity, FlowPerActivity, MoneyPerFlow); impl_div!(MoneyPerCapacity, Year, MoneyPerCapacityPerYear); impl_div!(FlowPerCapacity, ActivityPerCapacity, FlowPerActivity); + +impl std::ops::Mul for Year { + type Output = Dimensionless; + + fn mul(self, by: PerYear) -> Self::Output { + Dimensionless(self.0 * by.0) + } +} + +impl std::ops::Mul for Activity { + type Output = ActivityPerYear; + + fn mul(self, by: PerYear) -> Self::Output { + ActivityPerYear(self.0 * by.0) + } +} + +impl std::ops::Mul for ActivityPerCapacity { + type Output = ActivityPerCapacityPerYear; + + fn mul(self, by: PerYear) -> Self::Output { + ActivityPerCapacityPerYear(self.0 * by.0) + } +} + +impl std::ops::Div for Dimensionless { + type Output = PerYear; + + fn div(self, rhs: Year) -> Self::Output { + PerYear(self.0 / rhs.0) + } +}