diff --git a/examples/muse1_default/process_investment_constraints.csv b/examples/muse1_default/process_investment_constraints.csv new file mode 100644 index 000000000..c8f9aeca3 --- /dev/null +++ b/examples/muse1_default/process_investment_constraints.csv @@ -0,0 +1,6 @@ +process_id,regions,commission_years,addition_limit +gassupply1,R1,all,10 +gasCCGT,R1,all,10 +windturbine,R1,all,10 +gasboiler,R1,all,10 +heatpump,R1,all,10 diff --git a/examples/two_regions/process_investment_constraints.csv b/examples/two_regions/process_investment_constraints.csv new file mode 100644 index 000000000..09b9404eb --- /dev/null +++ b/examples/two_regions/process_investment_constraints.csv @@ -0,0 +1,6 @@ +process_id,regions,commission_years,addition_limit +gassupply1,R1;R2,all,10 +gasCCGT,R1;R2,all,10 +windturbine,R1;R2,all,10 +gasboiler,R1;R2,all,10 +heatpump,R1;R2,all,10 diff --git a/src/asset.rs b/src/asset.rs index f7fe42873..11b963377 100644 --- a/src/asset.rs +++ b/src/asset.rs @@ -194,6 +194,22 @@ impl AssetCapacity { } } + /// Create an `AssetCapacity` from a total capacity and optional unit size + /// + /// If a unit size is provided, the capacity is represented as a discrete number of units, + /// calculated as the floor of (capacity / `unit_size`). If no unit size is provided, the + /// capacity is represented as continuous. + pub fn from_capacity_floor(capacity: Capacity, unit_size: Option) -> Self { + #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] + match unit_size { + Some(size) => { + let num_units = (capacity / size).value().floor() as u32; + AssetCapacity::Discrete(num_units, size) + } + None => AssetCapacity::Continuous(capacity), + } + } + /// Returns the total capacity represented by this `AssetCapacity`. pub fn total_capacity(&self) -> Capacity { match self { @@ -1046,6 +1062,30 @@ impl Asset { let child_asset = AssetRef::from(Rc::new(child_asset)); std::iter::repeat_n(child_asset, n_units as usize).collect() } + + /// For uncommissioned assets, get the maximum capacity permitted to be installed based on the + /// investment constraints for the asset's process. + /// + /// The limit is taken from the process's investment constraints for the asset's region and + /// commission year, and the portion of the commodity demand being considered. + /// + /// For divisible assets, the returned capacity will be rounded down to the nearest multiple of + /// the asset's unit size. + pub fn max_installable_capacity( + &self, + commodity_portion: Dimensionless, + ) -> Option { + assert!( + !self.is_commissioned(), + "max_installable_capacity can only be called on uncommissioned assets" + ); + + self.process + .investment_constraints + .get(&(self.region_id.clone(), self.commission_year)) + .and_then(|c| c.get_addition_limit().map(|l| l * commodity_portion)) + .map(|limit| AssetCapacity::from_capacity_floor(limit, self.unit_size())) + } } #[allow(clippy::missing_fields_in_debug)] @@ -1524,6 +1564,27 @@ mod tests { assert_eq!(got.total_capacity(), expected_total); } + #[rstest] + #[case::exact_multiple(Capacity(12.0), Some(Capacity(4.0)), Some(3), Capacity(12.0))] + #[case::rounded_down(Capacity(11.0), Some(Capacity(4.0)), Some(2), Capacity(8.0))] + #[case::unit_size_greater_than_capacity( + Capacity(3.0), + Some(Capacity(4.0)), + Some(0), + Capacity(0.0) + )] + #[case::continuous(Capacity(5.5), None, None, Capacity(5.5))] + fn from_capacity_floor( + #[case] capacity: Capacity, + #[case] unit_size: Option, + #[case] expected_n: Option, + #[case] expected_total: Capacity, + ) { + let got = AssetCapacity::from_capacity_floor(capacity, unit_size); + assert_eq!(got.n_units(), expected_n); + assert_eq!(got.total_capacity(), expected_total); + } + #[rstest] #[case::round_up(3u32, Capacity(4.0), Dimensionless(0.5), 2u32)] #[case::exact(3u32, Capacity(4.0), Dimensionless(0.33), 1u32)] @@ -2303,4 +2364,25 @@ mod tests { "Agent A0_GEX has asset with commission year 2060, not within process GASDRV commission years: 2020..=2050" ); } + + #[rstest] + fn max_installable_capacity(mut process: Process, region_id: RegionID) { + // Set an addition limit of 3 for (region, year 2015) + process.investment_constraints.insert( + (region_id.clone(), 2015), + Rc::new(crate::process::ProcessInvestmentConstraint { + addition_limit: Some(Capacity(3.0)), + }), + ); + let process_rc = Rc::new(process); + + // Create a candidate asset with commission year 2015 + let asset = + Asset::new_candidate(process_rc.clone(), region_id.clone(), Capacity(1.0), 2015) + .unwrap(); + + // commodity_portion = 0.5 -> limit = 3 * 0.5 = 1.5 + let result = asset.max_installable_capacity(Dimensionless(0.5)); + assert_eq!(result, Some(AssetCapacity::Continuous(Capacity(1.5)))); + } } diff --git a/src/input/process/investment_constraints.rs b/src/input/process/investment_constraints.rs index aaacf07bf..ca5490b38 100644 --- a/src/input/process/investment_constraints.rs +++ b/src/input/process/investment_constraints.rs @@ -5,6 +5,7 @@ use crate::process::{ ProcessID, ProcessInvestmentConstraint, ProcessInvestmentConstraintsMap, ProcessMap, }; use crate::region::parse_region_str; +use crate::units::{Capacity, Dimensionless}; use crate::year::parse_year_str; use anyhow::{Context, Result, ensure}; use itertools::iproduct; @@ -21,7 +22,7 @@ struct ProcessInvestmentConstraintRaw { process_id: String, regions: String, commission_years: String, - addition_limit: f64, + addition_limit: Capacity, } impl ProcessInvestmentConstraintRaw { @@ -29,7 +30,7 @@ impl ProcessInvestmentConstraintRaw { fn validate(&self) -> Result<()> { // Validate that value is finite ensure!( - self.addition_limit.is_finite() && self.addition_limit >= 0.0, + self.addition_limit.is_finite() && self.addition_limit >= Capacity(0.0), "Invalid value for addition constraint: '{}'; must be non-negative and finite.", self.addition_limit ); @@ -118,11 +119,29 @@ where })?; // Create constraints for each region and year combination - let constraint = Rc::new(ProcessInvestmentConstraint { - addition_limit: Some(record.addition_limit), - }); + // For a given milestone year, the addition limit should be multiplied + // by the number of years since the previous milestone year. Any + // addition limits specified for the first milestone year are ignored. let process_map = map.entry(process_id.clone()).or_default(); for (region, &year) in iproduct!(&record_regions, &constraint_years) { + // Calculate years since previous milestone year + // We can ignore constraints in the first milestone year as no investments are performed then + let idx = milestone_years.iter().position(|y| *y == year).expect( + "Year should be in milestone_years since it was validated by parse_year_str", + ); + if idx == 0 { + continue; + } + let prev_year = milestone_years[idx - 1]; + let years_since_prev = year - prev_year; + + // Multiply the addition limit by the number of years since previous milestone. + let scaled_limit = record.addition_limit * Dimensionless(years_since_prev as f64); + + let constraint = Rc::new(ProcessInvestmentConstraint { + addition_limit: Some(scaled_limit), + }); + try_insert(process_map, &(region.clone(), year), constraint.clone())?; } } @@ -136,7 +155,7 @@ mod tests { use crate::region::RegionID; use rstest::rstest; - fn validate_raw_constraint(addition_limit: f64) -> Result<()> { + fn validate_raw_constraint(addition_limit: Capacity) -> Result<()> { let constraint = ProcessInvestmentConstraintRaw { process_id: "test_process".into(), regions: "ALL".into(), @@ -155,7 +174,7 @@ mod tests { process_id: "process1".into(), regions: "GBR".into(), commission_years: "ALL".into(), // Should apply to milestone years [2012, 2016] - addition_limit: 100.0, + addition_limit: Capacity(100.0), }]; let result = read_process_investment_constraints_from_iter( @@ -201,19 +220,19 @@ mod tests { process_id: "process1".into(), regions: "GBR".into(), commission_years: "2010".into(), - addition_limit: 100.0, + addition_limit: Capacity(100.0), }, ProcessInvestmentConstraintRaw { process_id: "process1".into(), regions: "ALL".into(), commission_years: "2015".into(), - addition_limit: 200.0, + addition_limit: Capacity(200.0), }, ProcessInvestmentConstraintRaw { process_id: "process1".into(), regions: "USA".into(), commission_years: "2020".into(), - addition_limit: 50.0, + addition_limit: Capacity(50.0), }, ]; @@ -234,32 +253,32 @@ mod tests { let gbr_region: RegionID = "GBR".into(); let usa_region: RegionID = "USA".into(); - // Check GBR 2010 constraint - let gbr_2010 = process_constraints - .get(&(gbr_region.clone(), 2010)) - .expect("GBR 2010 constraint should exist"); - assert_eq!(gbr_2010.addition_limit, Some(100.0)); + // GBR 2010 constraint is for the first milestone year and should be ignored + assert!( + !process_constraints.contains_key(&(gbr_region.clone(), 2010)), + "GBR 2010 constraint should not exist" + ); - // Check GBR 2015 constraint (from ALL regions) + // Check GBR 2015 constraint (from ALL regions), scaled by years since previous milestone (5 years) let gbr_2015 = process_constraints .get(&(gbr_region, 2015)) .expect("GBR 2015 constraint should exist"); - assert_eq!(gbr_2015.addition_limit, Some(200.0)); + assert_eq!(gbr_2015.addition_limit, Some(Capacity(200.0 * 5.0))); - // Check USA 2015 constraint (from ALL regions) + // Check USA 2015 constraint (from ALL regions), scaled by 5 years let usa_2015 = process_constraints .get(&(usa_region.clone(), 2015)) .expect("USA 2015 constraint should exist"); - assert_eq!(usa_2015.addition_limit, Some(200.0)); + assert_eq!(usa_2015.addition_limit, Some(Capacity(200.0 * 5.0))); - // Check USA 2020 constraint + // Check USA 2020 constraint, scaled by years since previous milestone (5 years) let usa_2020 = process_constraints .get(&(usa_region, 2020)) .expect("USA 2020 constraint should exist"); - assert_eq!(usa_2020.addition_limit, Some(50.0)); + assert_eq!(usa_2020.addition_limit, Some(Capacity(50.0 * 5.0))); - // Verify total number of constraints (2 GBR + 2 USA = 4) - assert_eq!(process_constraints.len(), 4); + // Verify total number of constraints (GBR 2015, USA 2015, USA 2020 = 3) + assert_eq!(process_constraints.len(), 3); } #[rstest] @@ -272,7 +291,7 @@ mod tests { process_id: "process1".into(), regions: "ALL".into(), commission_years: "ALL".into(), - addition_limit: 75.0, + addition_limit: Capacity(75.0), }]; // Read constraints into the map @@ -292,21 +311,22 @@ mod tests { let gbr_region: RegionID = "GBR".into(); let usa_region: RegionID = "USA".into(); - // Verify constraint exists for all region-year combinations - for &year in &milestone_years { + // Verify constraint exists for all region-year combinations except the first milestone year + for &year in &milestone_years[1..] { let gbr_constraint = process_constraints .get(&(gbr_region.clone(), year)) .unwrap_or_else(|| panic!("GBR {year} constraint should exist")); - assert_eq!(gbr_constraint.addition_limit, Some(75.0)); + // scaled by years since previous milestone (5 years) + assert_eq!(gbr_constraint.addition_limit, Some(Capacity(75.0 * 5.0))); let usa_constraint = process_constraints .get(&(usa_region.clone(), year)) .unwrap_or_else(|| panic!("USA {year} constraint should exist")); - assert_eq!(usa_constraint.addition_limit, Some(75.0)); + assert_eq!(usa_constraint.addition_limit, Some(Capacity(75.0 * 5.0))); } - // Verify total number of constraints (2 regions × 3 years = 6) - assert_eq!(process_constraints.len(), 6); + // Verify total number of constraints (2 regions × 2 years = 4) + assert_eq!(process_constraints.len(), 4); } #[rstest] @@ -319,7 +339,7 @@ mod tests { process_id: "process1".into(), regions: "GBR".into(), commission_years: "2025".into(), // Outside milestone years (2010-2020) - addition_limit: 100.0, + addition_limit: Capacity(100.0), }]; // Should fail with milestone year validation error @@ -337,15 +357,15 @@ mod tests { #[test] fn validate_addition_with_finite_value() { // Valid: addition constraint with positive value - let valid = validate_raw_constraint(10.0); + let valid = validate_raw_constraint(Capacity(10.0)); valid.unwrap(); // Valid: addition constraint with zero value - let valid = validate_raw_constraint(0.0); + let valid = validate_raw_constraint(Capacity(0.0)); valid.unwrap(); // Not valid: addition constraint with negative value - let invalid = validate_raw_constraint(-10.0); + let invalid = validate_raw_constraint(Capacity(-10.0)); assert_error!( invalid, "Invalid value for addition constraint: '-10'; must be non-negative and finite." @@ -355,14 +375,14 @@ mod tests { #[test] fn validate_addition_rejects_infinite() { // Invalid: infinite value - let invalid = validate_raw_constraint(f64::INFINITY); + let invalid = validate_raw_constraint(Capacity(f64::INFINITY)); assert_error!( invalid, "Invalid value for addition constraint: 'inf'; must be non-negative and finite." ); // Invalid: NaN value - let invalid = validate_raw_constraint(f64::NAN); + let invalid = validate_raw_constraint(Capacity(f64::NAN)); assert_error!( invalid, "Invalid value for addition constraint: 'NaN'; must be non-negative and finite." diff --git a/src/process.rs b/src/process.rs index f12986e16..6cea357f1 100644 --- a/src/process.rs +++ b/src/process.rs @@ -497,12 +497,25 @@ pub struct ProcessParameter { } /// A constraint imposed on investments in the process +/// +/// Constraints apply to a specific milestone year, and have been pre-scaled from annual limits +/// defined in the input data to account for the number of years since the previous milestone year. #[derive(PartialEq, Debug, Clone)] pub struct ProcessInvestmentConstraint { - /// Addition constraint: Yearly limit an agent can invest - /// in the process, shared according to the agent's - /// proportion of the processes primary commodity demand - pub addition_limit: Option, + /// Addition constraint: Limit an agent can invest in the process, shared according to the + /// agent's proportion of the process's primary commodity demand + pub addition_limit: Option, +} + +impl ProcessInvestmentConstraint { + /// Calculate the effective addition limit + /// + /// For now, this just returns `addition_limit`, but in the future when we add growth + /// limits and total capacity limits, this will have more complex logic which will depend on the + /// current total capacity. + pub fn get_addition_limit(&self) -> Option { + self.addition_limit + } } #[cfg(test)] diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index d7a927923..e5d9e9088 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -290,12 +290,17 @@ fn select_assets_for_single_market( region_id, year, ) - .collect(); + .collect::>(); + + // Calculate investment limits for candidate assets + let investment_limits = + calculate_investment_limits_for_candidates(&opt_assets, commodity_portion); // Choose assets from among existing pool and candidates let best_assets = select_best_assets( model, opt_assets, + investment_limits, commodity, agent, region_id, @@ -376,11 +381,28 @@ fn select_assets_for_cycle( .cloned() .collect(); + // Retrieve installable capacity limits for flexible capacity assets. + let key = (commodity_id.clone(), year); + let mut agent_share_cache = HashMap::new(); + let capacity_limits = flexible_capacity_assets + .iter() + .filter_map(|asset| { + let agent_id = asset.agent_id().unwrap().clone(); + let agent_share = *agent_share_cache + .entry(agent_id.clone()) + .or_insert_with(|| model.agents[&agent_id].commodity_portions[&key]); + asset + .max_installable_capacity(agent_share) + .map(|max_capacity| (asset.clone(), max_capacity)) + }) + .collect::>(); + // Run dispatch let solution = DispatchRun::new(model, &all_assets, year) .with_market_balance_subset(&markets_to_balance) .with_flexible_capacity_assets( &flexible_capacity_assets, + Some(&capacity_limits), // Gives newly selected cycle assets limited capacity wiggle-room; existing assets stay fixed. model.parameters.capacity_margin, ) @@ -668,11 +690,40 @@ fn warn_on_equal_appraisal_outputs( } } +/// Calculate investment limits for an agent's candidate assets in a given year +/// +/// Investment limits are based on demand for the commodity (capacity cannot exceed that needed to +/// meet demand), and any annual addition limits specified by the process (scaled according to the +/// agent's portion of the commodity demand and the number of years elapsed since the previous +/// milestone year). +fn calculate_investment_limits_for_candidates( + opt_assets: &[AssetRef], + commodity_portion: Dimensionless, +) -> HashMap { + // Calculate limits for each candidate asset + opt_assets + .iter() + .filter(|asset| !asset.is_commissioned()) + .map(|asset| { + // Start off with the demand-limiting capacity (pre-calculated when creating candidate) + let mut cap = asset.capacity(); + + // Cap by the addition limits of the process, if specified + if let Some(limit_capacity) = asset.max_installable_capacity(commodity_portion) { + cap = cap.min(limit_capacity); + } + + (asset.clone(), cap) + }) + .collect() +} + /// Get the best assets for meeting demand for the given commodity #[allow(clippy::too_many_arguments)] fn select_best_assets( model: &Model, mut opt_assets: Vec, + investment_limits: HashMap, commodity: &Commodity, agent: &Agent, region_id: &RegionID, @@ -682,25 +733,22 @@ fn select_best_assets( writer: &mut DataWriter, ) -> Result> { let objective_type = &agent.objectives[&year]; + let mut remaining_candidate_capacity = investment_limits; // Calculate coefficients for all asset options according to the agent's objective let coefficients = calculate_coefficients_for_assets(model, objective_type, &opt_assets, prices, year); - let mut remaining_candidate_capacity = HashMap::from_iter( - opt_assets - .iter() - .filter(|asset| !asset.is_commissioned()) - .map(|asset| (asset.clone(), asset.capacity())), - ); - + // Iteratively select the best asset until demand is met let mut round = 0; let mut best_assets: Vec = Vec::new(); while is_any_remaining_demand(&demand) { ensure!( !opt_assets.is_empty(), - "Failed to meet demand for commodity '{}' with provided assets", - &commodity.id + "Failed to meet demand for commodity '{}' in region '{}' with provided investment \ + options. This may be due to overly restrictive process investment constraints.", + &commodity.id, + region_id ); // Since all assets with the same `group_id` are identical, we only need to appraise one @@ -710,13 +758,14 @@ fn select_best_assets( // Appraise all options let mut outputs_for_opts = Vec::new(); for asset in &opt_assets { + // For candidates, determine the maximum capacity that can be invested in this round, + // according to the tranche size and remaining capacity limits. let max_capacity = (!asset.is_commissioned()).then(|| { - let max_capacity = asset + let tranche_capacity = asset .capacity() .apply_limit_factor(model.parameters.capacity_limit_factor); - let remaining_capacity = remaining_candidate_capacity[asset]; - max_capacity.min(remaining_capacity) + tranche_capacity.min(remaining_capacity) }); // Skip any assets from groups we've already seen @@ -788,7 +837,7 @@ fn select_best_assets( best_output.capacity.total_capacity() ); - // Update the assets + // Update the assets and remaining candidate capacity update_assets( best_output.asset, best_output.capacity, @@ -885,6 +934,7 @@ mod tests { use crate::process::{ActivityLimits, FlowType, Process, ProcessFlow}; use crate::region::RegionID; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; + use crate::units::Dimensionless; use crate::units::{Capacity, Flow, FlowPerActivity, MoneyPerFlow}; use indexmap::indexmap; use rstest::rstest; diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 0797e4161..7b082c873 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -17,6 +17,7 @@ use anyhow::{Result, bail, ensure}; use highs::{HighsModelStatus, HighsStatus, RowProblem as Problem, Sense}; use indexmap::{IndexMap, IndexSet}; use itertools::{chain, iproduct}; +use std::collections::HashMap; use std::error::Error; use std::fmt; use std::ops::Range; @@ -393,6 +394,7 @@ pub struct DispatchRun<'model, 'run> { model: &'model Model, existing_assets: &'run [AssetRef], flexible_capacity_assets: &'run [AssetRef], + capacity_limits: Option<&'run HashMap>, candidate_assets: &'run [AssetRef], markets_to_balance: &'run [(CommodityID, RegionID)], input_prices: Option<&'run CommodityPrices>, @@ -407,6 +409,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { model, existing_assets: assets, flexible_capacity_assets: &[], + capacity_limits: None, candidate_assets: &[], markets_to_balance: &[], input_prices: None, @@ -419,10 +422,12 @@ impl<'model, 'run> DispatchRun<'model, 'run> { pub fn with_flexible_capacity_assets( self, flexible_capacity_assets: &'run [AssetRef], + capacity_limits: Option<&'run HashMap>, capacity_margin: f64, ) -> Self { Self { flexible_capacity_assets, + capacity_limits, capacity_margin, ..self } @@ -583,6 +588,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { &mut problem, &mut variables.capacity_vars, self.flexible_capacity_assets, + self.capacity_limits, self.capacity_margin, ); } @@ -647,6 +653,7 @@ fn add_capacity_variables( problem: &mut Problem, variables: &mut CapacityVariableMap, assets: &[AssetRef], + capacity_limits: Option<&HashMap>, capacity_margin: f64, ) -> Range { // This line **must** come before we add more variables @@ -662,17 +669,41 @@ fn add_capacity_variables( let current_capacity = asset.capacity(); let coeff = calculate_capacity_coefficient(asset); + // Retrieve capacity limit if provided + let capacity_limit = capacity_limits.and_then(|limits| limits.get(asset)); + + // Sanity check: make sure capacity_limit is compatible with current_capacity + if let Some(limit) = capacity_limit { + assert!( + matches!( + (current_capacity, limit), + (AssetCapacity::Continuous(_), AssetCapacity::Continuous(_)) + | (AssetCapacity::Discrete(_, _), AssetCapacity::Discrete(_, _)) + ), + "Incompatible capacity types for asset capacity limit" + ); + } + + // Add a capacity variable for each asset + // Bounds are calculated based on current capacity with wiggle-room defined by + // `capacity_margin`, and limited by `capacity_limit` if provided. let var = match current_capacity { AssetCapacity::Continuous(cap) => { // Continuous capacity: capacity variable represents total capacity let lower = ((1.0 - capacity_margin) * cap.value()).max(0.0); - let upper = (1.0 + capacity_margin) * cap.value(); + let mut upper = (1.0 + capacity_margin) * cap.value(); + if let Some(limit) = capacity_limit { + upper = upper.min(limit.total_capacity().value()); + } problem.add_column(coeff.value(), lower..=upper) } AssetCapacity::Discrete(units, unit_size) => { // Discrete capacity: capacity variable represents number of units let lower = ((1.0 - capacity_margin) * units as f64).max(0.0); - let upper = (1.0 + capacity_margin) * units as f64; + let mut upper = (1.0 + capacity_margin) * units as f64; + if let Some(limit) = capacity_limit { + upper = upper.min(limit.n_units().unwrap() as f64); + } problem.add_integer_column((coeff * unit_size).value(), lower..=upper) } };