From 5794ee1b01140dc1832e8ff2dac171b03aa2df30 Mon Sep 17 00:00:00 2001 From: Alex Dewar Date: Mon, 15 Jun 2026 18:04:35 +0100 Subject: [PATCH 1/6] fixture.rs: Unwrap `Result`s rather than asserting OK in tests This way, you get a more informative error message if a test fails. --- src/fixture.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fixture.rs b/src/fixture.rs index 8ac2e2bdf..23a1649a2 100644 --- a/src/fixture.rs +++ b/src/fixture.rs @@ -75,7 +75,7 @@ pub(crate) use patch_and_validate_simple; /// Check whether validation succeeds for simple example with patches macro_rules! assert_validate_ok_simple { ($file_patches:expr) => { - assert!(crate::fixture::patch_and_validate_simple!($file_patches).is_ok()) + crate::fixture::patch_and_validate_simple!($file_patches).unwrap(); }; } pub(crate) use assert_validate_ok_simple; @@ -116,7 +116,7 @@ pub(crate) use patch_and_run_simple; /// Check whether the simple example runs successfully after applying file patches macro_rules! assert_patched_runs_ok_simple { ($file_patches:expr) => { - assert!(crate::fixture::patch_and_run_simple!($file_patches).is_ok()) + crate::fixture::patch_and_run_simple!($file_patches).unwrap(); }; } pub(crate) use assert_patched_runs_ok_simple; From c7fa47a62d517ebc8f6fdf43145714e3ccae8771 Mon Sep 17 00:00:00 2001 From: Alex Dewar Date: Mon, 15 Jun 2026 18:23:51 +0100 Subject: [PATCH 2/6] regenerate_test_data.sh: Don't bail early if a model fails --- tests/regenerate_test_data.sh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/regenerate_test_data.sh b/tests/regenerate_test_data.sh index 052a33063..7fade8103 100755 --- a/tests/regenerate_test_data.sh +++ b/tests/regenerate_test_data.sh @@ -25,7 +25,8 @@ run_example() { env MUSE2_LOG_LEVEL=error MUSE2_USE_DEFAULT_SETTINGS=1 \ cargo -q run example run -o "data/$example" "$example" \ - --overwrite --no-copy-input-files $@ + --overwrite --no-copy-input-files $@ \ + || echo ERROR: Failed to run example "$example" > /dev/stderr } for example in $examples; do From 6f804755e23c8929f351387c24ced5afc068e5a9 Mon Sep 17 00:00:00 2001 From: Alex Dewar Date: Tue, 16 Jun 2026 09:08:47 +0100 Subject: [PATCH 3/6] Print number of non-feasible assets ignored when investment fails --- src/simulation/investment.rs | 6 ++++-- src/simulation/investment/appraisal.rs | 15 ++++++++++++--- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index 025ae59cd..3f4a021f3 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -807,7 +807,7 @@ fn select_best_assets( )?; // Sort by investment priority and discard non-feasible options - sort_and_filter_appraisal_outputs(&mut outputs_for_opts); + let num_nonfeasible = sort_and_filter_appraisal_outputs(&mut outputs_for_opts); // Check if there are any remaining options. If not, we cannot meet demand, so have to bail // out. @@ -820,12 +820,14 @@ fn select_best_assets( bail!( "No feasible investment options left for \ - commodity '{}', region '{}', year '{}', agent '{}' after appraisal.\n\ + commodity '{}', region '{}', year '{}', agent '{}' after appraisal. \ + {} non-feasible options were not considered.\n\ Remaining unmet demand (time_slice : flow):\n{}", &commodity.id, region_id, year, agent.id, + num_nonfeasible, remaining_demands.join("\n") ); } diff --git a/src/simulation/investment/appraisal.rs b/src/simulation/investment/appraisal.rs index b509b1a7b..d81fe14f1 100644 --- a/src/simulation/investment/appraisal.rs +++ b/src/simulation/investment/appraisal.rs @@ -373,13 +373,22 @@ fn compare_asset_fallback(asset1: &Asset, asset2: &Asset) -> Ordering { /// with invalid metrics (e.g. `None`) as well as zero capacity. This avoids meaningless or `NaN` /// appraisal metrics that could cause the program to panic, so the length of the returned vector /// may be less than the input. -pub fn sort_and_filter_appraisal_outputs(outputs_for_opts: &mut Vec) { - outputs_for_opts.retain(AppraisalOutput::is_valid); - outputs_for_opts.sort_by(|output1, output2| match output1.compare_metric(output2) { +/// +/// # Returns +/// +/// Returns the number of non-feasible assets which were removed. +pub fn sort_and_filter_appraisal_outputs(outputs: &mut Vec) -> usize { + let old_len = outputs.len(); + outputs.retain(AppraisalOutput::is_valid); + let num_nonfeasible = old_len - outputs.len(); + + outputs.sort_by(|output1, output2| match output1.compare_metric(output2) { // If equal, we fall back on comparing asset properties Ordering::Equal => compare_asset_fallback(&output1.asset, &output2.asset), cmp => cmp, }); + + num_nonfeasible } /// Counts the number of top appraisal outputs in a sorted slice that are indistinguishable From 919a10bf0a5efbd13803a74a48f0238895fbc8d6 Mon Sep 17 00:00:00 2001 From: Alex Dewar Date: Thu, 28 May 2026 13:08:08 +0100 Subject: [PATCH 4/6] Change LCOX to use same optimisation as NPV Closes #1199. Co-authored-by: Tom Bland --- src/fixture.rs | 4 +- src/output.rs | 7 +- src/simulation/investment/appraisal.rs | 32 ++--- .../investment/appraisal/coefficients.rs | 133 ++++++------------ .../investment/appraisal/constraints.rs | 116 ++------------- .../investment/appraisal/optimisation.rs | 131 ++++++++--------- 6 files changed, 120 insertions(+), 303 deletions(-) diff --git a/src/fixture.rs b/src/fixture.rs index 23a1649a2..e212a0c12 100644 --- a/src/fixture.rs +++ b/src/fixture.rs @@ -21,7 +21,7 @@ use crate::simulation::investment::appraisal::{ use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceLevel}; use crate::units::{ Activity, ActivityPerCapacity, Capacity, Dimensionless, Flow, MoneyPerActivity, - MoneyPerCapacity, MoneyPerCapacityPerYear, MoneyPerFlow, Year, + MoneyPerCapacity, MoneyPerCapacityPerYear, Year, }; use anyhow::Result; use indexmap::indexmap; @@ -408,10 +408,8 @@ pub fn appraisal_output(asset: Asset, time_slice: TimeSliceID) -> AppraisalOutpu asset: AssetRef::from(asset), capacity: AssetCapacity::Continuous(Capacity(42.0)), coefficients: Rc::new(ObjectiveCoefficients { - capacity_coefficient: MoneyPerCapacity(2.14), activity_coefficients, market_costs, - unmet_demand_coefficient: MoneyPerFlow(10000.0), }), activity, unmet_demand, diff --git a/src/output.rs b/src/output.rs index bb323c0ed..29d00263d 100644 --- a/src/output.rs +++ b/src/output.rs @@ -8,9 +8,7 @@ use crate::simulation::investment::appraisal::AppraisalOutput; use crate::simulation::optimisation::{FlowMap, Solution}; use crate::simulation::prices::PriceMap; use crate::time_slice::TimeSliceID; -use crate::units::{ - Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow, -}; +use crate::units::{Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerFlow}; use anyhow::{Context, Result, ensure}; use csv; use indexmap::IndexMap; @@ -261,7 +259,6 @@ struct AppraisalResultsRow { process_id: ProcessID, region_id: RegionID, capacity: Capacity, - capacity_coefficient: MoneyPerCapacity, metric: Option, } @@ -495,7 +492,6 @@ impl DebugDataWriter { process_id: result.asset.process_id().clone(), region_id: result.asset.region_id().clone(), capacity: result.capacity.total_capacity(), - capacity_coefficient: result.coefficients.capacity_coefficient, metric: result.metric.as_ref().map(|m| m.value()), }; self.appraisal_results_writer.serialize(row)?; @@ -1195,7 +1191,6 @@ mod tests { process_id: asset.process_id().clone(), region_id: asset.region_id().clone(), capacity: Capacity(42.0), - capacity_coefficient: MoneyPerCapacity(2.14), metric: Some(4.14), }; let records: Vec = diff --git a/src/simulation/investment/appraisal.rs b/src/simulation/investment/appraisal.rs index d81fe14f1..30f8cd760 100644 --- a/src/simulation/investment/appraisal.rs +++ b/src/simulation/investment/appraisal.rs @@ -258,26 +258,19 @@ fn calculate_lcox( coefficients: &Rc, demand: &DemandMap, ) -> Result { - let results = perform_optimisation( - model, - asset, - max_capacity, - commodity, - coefficients, - demand, - highs::Sense::Minimise, - )?; + let results = + perform_optimisation(model, asset, max_capacity, commodity, coefficients, demand)?; let cost_index = lcox( - results.capacity.total_capacity(), - coefficients.capacity_coefficient, + max_capacity.total_capacity(), + annual_fixed_cost(asset), &results.activity, &coefficients.market_costs, ); Ok(AppraisalOutput::new( asset.clone(), - results.capacity, + max_capacity, results, cost_index.map(LCOXMetric::new), coefficients.clone(), @@ -297,15 +290,8 @@ fn calculate_npv( coefficients: &Rc, demand: &DemandMap, ) -> Result { - let results = perform_optimisation( - model, - asset, - max_capacity, - commodity, - coefficients, - demand, - highs::Sense::Maximise, - )?; + let results = + perform_optimisation(model, asset, max_capacity, commodity, coefficients, demand)?; let annual_fixed_cost = annual_fixed_cost(asset); assert!( @@ -414,7 +400,7 @@ mod tests { use crate::fixture::{agent_id, asset, process, region_id}; use crate::process::Process; use crate::region::RegionID; - use crate::units::{Money, MoneyPerActivity, MoneyPerFlow}; + use crate::units::{Money, MoneyPerActivity}; use float_cmp::assert_approx_eq; use rstest::rstest; use std::rc::Rc; @@ -583,10 +569,8 @@ mod tests { fn objective_coeffs() -> Rc { Rc::new(ObjectiveCoefficients { - capacity_coefficient: MoneyPerCapacity(0.0), activity_coefficients: IndexMap::new(), market_costs: IndexMap::new(), - unmet_demand_coefficient: MoneyPerFlow(0.0), }) } diff --git a/src/simulation/investment/appraisal/coefficients.rs b/src/simulation/investment/appraisal/coefficients.rs index fe59c3de7..b85d330c3 100644 --- a/src/simulation/investment/appraisal/coefficients.rs +++ b/src/simulation/investment/appraisal/coefficients.rs @@ -1,12 +1,11 @@ //! Calculation of cost coefficients for investment tools. -use super::costs::annual_fixed_cost; use crate::agent::ObjectiveType; use crate::asset::AssetRef; use crate::model::Model; use crate::simulation::PriceMap; use crate::simulation::prices::Prices; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::{MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow}; +use crate::units::MoneyPerActivity; use indexmap::IndexMap; use std::collections::HashMap; use std::rc::Rc; @@ -18,17 +17,22 @@ use std::rc::Rc; /// coefficients used in the appraisal optimisation, together with the unmet-demand penalty. #[derive(Clone)] pub struct ObjectiveCoefficients { - /// Cost per unit of capacity - pub capacity_coefficient: MoneyPerCapacity, /// Cost per unit of activity in each time slice pub activity_coefficients: IndexMap, /// Market costs associated with asset for each time slice pub market_costs: IndexMap, - /// Unmet demand coefficient - pub unmet_demand_coefficient: MoneyPerFlow, } /// Calculates cost coefficients for a set of assets for a given objective type. +/// +/// Activity coefficients are revenue (including primary output) minus operating cost; a small +/// positive epsilon is added to activity coefficients so that assets with near-zero net value still +/// appear in dispatch. Capacity costs and unmet-demand penalties are set to zero for the NPV +/// objective. These activity coefficients are calculated using shadow prices. +/// +/// For NPV, "market costs" are calculated in the same way, except thew prices used are market +/// prices. For LCOX, a slightly different calculation is performed: the sign is inverted (as it +/// represents a cost) and the primary output (commodity of interest) is excluded. pub fn calculate_coefficients_for_assets( model: &Model, objective_type: &ObjectiveType, @@ -39,73 +43,22 @@ pub fn calculate_coefficients_for_assets( assets .iter() .map(|asset| { - let coefficient = match objective_type { - ObjectiveType::LevelisedCostOfX => calculate_coefficients_for_lcox( - asset, - &model.time_slice_info, - prices, - model.parameters.value_of_lost_load, - year, - ), - ObjectiveType::NetPresentValue => { - calculate_coefficients_for_npv(asset, &model.time_slice_info, prices, year) - } - }; + let coefficient = calculate_coefficients_for_asset( + asset, + objective_type, + &model.time_slice_info, + prices, + year, + ); (asset.clone(), Rc::new(coefficient)) }) .collect() } -/// Calculates the cost coefficients for LCOX. -/// -/// For LCOX the activity coefficient is calculated as operating cost minus revenue from -/// non-primary flows. The unmet demand coefficient is set from the model parameter -/// `value_of_lost_load`. -pub fn calculate_coefficients_for_lcox( - asset: &AssetRef, - time_slice_info: &TimeSliceInfo, - prices: &Prices, - value_of_lost_load: MoneyPerFlow, - year: u32, -) -> ObjectiveCoefficients { - // Capacity coefficient - let capacity_coefficient = annual_fixed_cost(asset); - - // Activity coefficients - let mut activity_coefficients = IndexMap::new(); - let mut market_costs = IndexMap::new(); - for time_slice in time_slice_info.iter_ids() { - // Get the operating cost of the asset. This includes the variable operating cost, levies and - // flow costs, but excludes costs/revenues from commodity consumption/production. - let operating_cost = asset.get_operating_cost(year, time_slice); - - let coefficient = - calculate_asset_costs_for_lcox(asset, operating_cost, time_slice, &prices.shadow); - activity_coefficients.insert(time_slice.clone(), coefficient); - let market_cost = - calculate_asset_costs_for_lcox(asset, operating_cost, time_slice, &prices.market); - market_costs.insert(time_slice.clone(), market_cost); - } - - // Unmet demand coefficient - let unmet_demand_coefficient = value_of_lost_load; - - ObjectiveCoefficients { - capacity_coefficient, - activity_coefficients, - market_costs, - unmet_demand_coefficient, - } -} - -/// Calculates the cost coefficients for NPV. -/// -/// For NPV the activity coefficient is revenue (including primary output) minus operating -/// cost; a small positive epsilon is added to activity coefficients so that assets with -/// near-zero net value still appear in dispatch. Capacity costs and unmet-demand penalties -/// are set to zero for the NPV objective. -pub fn calculate_coefficients_for_npv( +/// Calculates cost coefficients for a single asset +pub fn calculate_coefficients_for_asset( asset: &AssetRef, + objective_type: &ObjectiveType, time_slice_info: &TimeSliceInfo, prices: &Prices, year: u32, @@ -123,51 +76,55 @@ pub fn calculate_coefficients_for_npv( let operating_cost = asset.get_operating_cost(year, time_slice); let coefficient = - calculate_asset_costs_for_npv(asset, operating_cost, time_slice, &prices.shadow); + calculate_asset_revenues(asset, operating_cost, time_slice, &prices.shadow); activity_coefficients.insert( time_slice.clone(), coefficient + EPSILON_ACTIVITY_COEFFICIENT, ); - let market_cost = - calculate_asset_costs_for_npv(asset, operating_cost, time_slice, &prices.market); + + let market_cost = match objective_type { + ObjectiveType::LevelisedCostOfX => { + calculate_asset_costs_for_lcox(asset, operating_cost, time_slice, &prices.market) + } + ObjectiveType::NetPresentValue => { + calculate_asset_revenues(asset, operating_cost, time_slice, &prices.market) + } + }; market_costs.insert(time_slice.clone(), market_cost); } - // Unmet demand coefficient (we don't apply a cost to unmet demand, so we set this to zero) - let unmet_demand_coefficient = MoneyPerFlow(0.0); - ObjectiveCoefficients { - capacity_coefficient: MoneyPerCapacity(0.0), - market_costs, activity_coefficients, - unmet_demand_coefficient, + market_costs, } } -/// Calculate a single activity coefficient for the LCOX objective for a given time slice. -fn calculate_asset_costs_for_lcox( +/// Calculate the revenue from all flows minus operating cost +fn calculate_asset_revenues( asset: &AssetRef, operating_cost: MoneyPerActivity, time_slice: &TimeSliceID, prices: &PriceMap, ) -> MoneyPerActivity { - // Revenue from flows excluding the primary output - let revenue_from_flows = asset.get_revenue_from_flows_excluding_primary(prices, time_slice); + // Revenue from flows including the primary output + let revenue_from_flows = asset.get_revenue_from_flows(prices, time_slice); - // The activity coefficient is the operating cost minus the revenue from non-primary flows - operating_cost - revenue_from_flows + // The activity coefficient is the revenue from flows minus the operating cost (net revenue) + revenue_from_flows - operating_cost } -/// Calculate a single activity coefficient for the NPV objective for a given time slice. -fn calculate_asset_costs_for_npv( +/// Calculate asset costs for LCOX objective. +/// +/// Excludes revenues from the primary output (commodity of interest). +fn calculate_asset_costs_for_lcox( asset: &AssetRef, operating_cost: MoneyPerActivity, time_slice: &TimeSliceID, prices: &PriceMap, ) -> MoneyPerActivity { - // Revenue from flows including the primary output - let revenue_from_flows = asset.get_revenue_from_flows(prices, time_slice); + // Revenue from flows excluding the primary output + let revenue_from_flows = asset.get_revenue_from_flows_excluding_primary(prices, time_slice); - // The activity coefficient is the revenue from flows minus the operating cost (net revenue) - revenue_from_flows - operating_cost + // The activity coefficient is the operating cost minus the revenue from non-primary flows + operating_cost - revenue_from_flows } diff --git a/src/simulation/investment/appraisal/constraints.rs b/src/simulation/investment/appraisal/constraints.rs index 750cb85a3..01270ace0 100644 --- a/src/simulation/investment/appraisal/constraints.rs +++ b/src/simulation/investment/appraisal/constraints.rs @@ -1,45 +1,13 @@ //! Constraints for the optimisation problem. use super::DemandMap; use super::optimisation::Variable; -use crate::asset::{AssetCapacity, AssetRef, AssetState}; +use crate::asset::{AssetCapacity, AssetRef}; use crate::commodity::Commodity; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; use crate::units::Flow; use highs::RowProblem as Problem; use indexmap::IndexMap; -/// Adds a capacity constraint to the problem. -/// -/// The behaviour depends on whether the asset is commissioned or a candidate: -/// - For a commissioned asset, the capacity is fixed. -/// - For a candidate asset, the capacity is variable between zero and an upper bound. -pub fn add_capacity_constraint( - problem: &mut Problem, - asset: &AssetRef, - max_capacity: AssetCapacity, - capacity_var: Variable, -) { - let capacity_limit = match max_capacity { - AssetCapacity::Continuous(cap) => cap.value(), - AssetCapacity::Discrete(units, _) => units as f64, - }; - - let bounds = match asset.state() { - AssetState::Commissioned { .. } => { - // Fixed capacity for commissioned assets - capacity_limit..=capacity_limit - } - AssetState::Candidate => { - // Variable capacity between 0 and max for candidate assets - 0.0..=capacity_limit - } - _ => panic!( - "add_capacity_constraint should only be called with Commissioned or Candidate assets" - ), - }; - problem.add_row(bounds, [(capacity_var, 1.0)]); -} - /// Adds activity constraints to the problem. /// /// Constrains the activity variables to be within the asset's activity limits. @@ -53,37 +21,13 @@ pub fn add_capacity_constraint( pub fn add_activity_constraints( problem: &mut Problem, 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, time_slice_info); - } - AssetState::Candidate => { - 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" - ), - } -} - -fn add_activity_constraints_for_existing( - problem: &mut Problem, - asset: &AssetRef, + max_capacity: AssetCapacity, activity_vars: &IndexMap, time_slice_info: &TimeSliceInfo, ) { - for (ts_selection, limits) in asset.iter_activity_limits() { - let limits = limits.start().value()..=limits.end().value(); + let capacity = max_capacity.total_capacity(); + for (ts_selection, limits) in asset.iter_activity_per_capacity_limits() { + let limits = (capacity * *limits.start()).value()..=(capacity * *limits.end()).value(); // Collect activity terms for the time slices in this selection let terms = ts_selection @@ -96,48 +40,11 @@ fn add_activity_constraints_for_existing( } } -fn add_activity_constraints_for_candidate( - problem: &mut Problem, - asset: &AssetRef, - capacity_var: Variable, - activity_vars: &IndexMap, - time_slice_info: &TimeSliceInfo, -) { - for (ts_selection, limits) in asset.iter_activity_per_capacity_limits() { - let mut upper_limit = limits.end().value(); - let mut lower_limit = limits.start().value(); - - // If the asset capacity is discrete, the capacity variable represents number of - // units, so we need to multiply the per-capacity limits by the unit size. - if let AssetCapacity::Discrete(_, unit_size) = asset.capacity() { - upper_limit *= unit_size.value(); - lower_limit *= unit_size.value(); - } - - // Collect capacity and activity terms - // We have a single capacity term, and activity terms for all time slices in the selection - 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: sum(activity) - (capacity * upper_limit_per_capacity) ≤ 0 - problem.add_row(..=0.0, &terms_upper); - - // Lower bound: sum(activity) - (capacity * lower_limit_per_capacity) ≥ 0 - problem.add_row(0.0.., &terms_lower); - } -} - /// Adds demand constraints to the problem. /// -/// Constrains supply to be less than or equal to demand. This is implemented as an equality -/// across each time-slice selection: supply (activity terms, scaled by flow coefficients) plus -/// the `unmet_demand` variables equals the total demand for that selection, so non-negative -/// `unmet_demand` enforces supply ≤ demand. The selections follow the commodity's balance level. +/// Constrains supply to be less than or equal to demand. One inequality constraint is added per +/// time-slice selection at the commodity's balance level, capping the sum of activity (scaled by +/// flow coefficients) to the total demand for that selection. pub fn add_demand_constraints( problem: &mut Problem, asset: &AssetRef, @@ -145,7 +52,6 @@ pub fn add_demand_constraints( time_slice_info: &TimeSliceInfo, demand: &DemandMap, activity_vars: &IndexMap, - unmet_demand_vars: &IndexMap, ) { for ts_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) { let mut demand_for_ts_selection = Flow(0.0); @@ -154,11 +60,7 @@ pub fn add_demand_constraints( demand_for_ts_selection += demand[time_slice]; let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff; terms.push((activity_vars[time_slice], flow_coeff.value())); - terms.push((unmet_demand_vars[time_slice], 1.0)); } - problem.add_row( - demand_for_ts_selection.value()..=demand_for_ts_selection.value(), - terms, - ); + problem.add_row(0.0..=demand_for_ts_selection.value(), terms); } } diff --git a/src/simulation/investment/appraisal/optimisation.rs b/src/simulation/investment/appraisal/optimisation.rs index 0cec3f760..d5cfe9e42 100644 --- a/src/simulation/investment/appraisal/optimisation.rs +++ b/src/simulation/investment/appraisal/optimisation.rs @@ -1,17 +1,16 @@ //! Optimisation problem for investment tools. use super::DemandMap; use super::ObjectiveCoefficients; -use super::constraints::{ - add_activity_constraints, add_capacity_constraint, add_demand_constraints, -}; -use crate::asset::{AssetCapacity, AssetRef}; +use super::constraints::{add_activity_constraints, add_demand_constraints}; +use crate::asset::AssetCapacity; +use crate::asset::AssetRef; use crate::commodity::Commodity; use crate::model::Model; use crate::simulation::optimisation::ModelError; use crate::simulation::optimisation::apply_highs_options_from_toml; use crate::simulation::optimisation::solve_optimal; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::{Activity, Capacity, Flow}; +use crate::units::{Activity, Dimensionless, Flow}; use anyhow::{Context, Result}; use highs::{RowProblem as Problem, Sense}; use indexmap::IndexMap; @@ -24,15 +23,8 @@ pub type Variable = highs::Col; /// Map storing variables for the optimisation problem struct VariableMap { - /// Capacity variable. - /// - /// This represents absolute capacity for indivisible assets and number of units for - /// divisible assets. - capacity_var: Variable, /// Activity variables in each time slice activity_vars: IndexMap, - /// Unmet demand variables - unmet_demand_vars: IndexMap, } impl VariableMap { @@ -44,24 +36,7 @@ impl VariableMap { /// /// # Returns /// A new `VariableMap` containing all created decision variables - fn add_to_problem( - problem: &mut Problem, - cost_coefficients: &ObjectiveCoefficients, - capacity_unit_size: Option, - ) -> Self { - // Create capacity variable with its associated cost - let capacity_coefficient = cost_coefficients.capacity_coefficient.value(); - let capacity_var = match capacity_unit_size { - Some(unit_size) => { - // Divisible asset: capacity variable represents number of units - problem.add_integer_column(capacity_coefficient * unit_size.value(), 0.0..) - } - None => { - // Indivisible asset: capacity variable represents total capacity - problem.add_column(capacity_coefficient, 0.0..) - } - }; - + fn add_to_problem(problem: &mut Problem, cost_coefficients: &ObjectiveCoefficients) -> Self { // Create activity variables for each time slice let mut activity_vars = IndexMap::new(); for (time_slice, cost) in &cost_coefficients.activity_coefficients { @@ -69,28 +44,15 @@ impl VariableMap { activity_vars.insert(time_slice.clone(), var); } - // Create unmet demand variables for each time slice - let mut unmet_demand_vars = IndexMap::new(); - for time_slice in cost_coefficients.activity_coefficients.keys() { - let var = problem.add_column(cost_coefficients.unmet_demand_coefficient.value(), 0.0..); - unmet_demand_vars.insert(time_slice.clone(), var); - } - - Self { - capacity_var, - activity_vars, - unmet_demand_vars, - } + Self { activity_vars } } } /// Map containing optimisation results and coefficients pub struct ResultsMap { - /// Capacity variable - pub capacity: AssetCapacity, /// Activity variables in each time slice pub activity: IndexMap, - /// Unmet demand variables + /// Remaining unmet demand per time slice, computed post-solve pub unmet_demand: DemandMap, } @@ -104,11 +66,10 @@ fn add_constraints( demand: &DemandMap, time_slice_info: &TimeSliceInfo, ) { - add_capacity_constraint(problem, asset, max_capacity, variables.capacity_var); add_activity_constraints( problem, asset, - variables.capacity_var, + max_capacity, &variables.activity_vars, time_slice_info, ); @@ -119,10 +80,45 @@ fn add_constraints( time_slice_info, demand, &variables.activity_vars, - &variables.unmet_demand_vars, ); } +/// Computes remaining unmet demand per time slice after a solve. +/// +/// For each time-slice selection at the commodity's balance level, the selection-level residual +/// (`demand_total - supply_total`, clamped to zero) is divided equally across the time slices in +/// the selection. +/// +/// The exact per-time-slice distribution is arbitrary: all downstream consumers sum values back up +/// to the selection level before using them (e.g. the next round's demand constraint, and +/// `is_any_remaining_demand`), so only the selection-level total needs to be correct. +fn compute_unmet_demand( + demand: &DemandMap, + activity: &IndexMap, + commodity: &Commodity, + asset: &AssetRef, + time_slice_info: &TimeSliceInfo, +) -> DemandMap { + let mut unmet_demand = DemandMap::new(); + let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff; + for ts_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) { + let time_slices: Vec<_> = ts_selection.iter(time_slice_info).collect(); + let demand_for_selection: Flow = time_slices.iter().map(|(ts, _)| demand[*ts]).sum(); + let supply_for_selection: Flow = time_slices + .iter() + .map(|(ts, _)| activity[*ts] * flow_coeff) + .sum(); + + #[allow(clippy::cast_precision_loss)] + let unmet_per_slice = (demand_for_selection - supply_for_selection).max(Flow(0.0)) + / Dimensionless(time_slices.len() as f64); + for (time_slice, _) in &time_slices { + unmet_demand.insert((*time_slice).clone(), unmet_per_slice); + } + } + unmet_demand +} + /// Performs optimisation for an asset, given the coefficients and demand. /// /// Will either maximise or minimise the objective function, depending on the `sense` parameter. @@ -135,11 +131,10 @@ pub fn perform_optimisation( commodity: &Commodity, coefficients: &ObjectiveCoefficients, demand: &DemandMap, - sense: Sense, ) -> Result { // Create problem and add variables let mut problem = Problem::default(); - let variables = VariableMap::add_to_problem(&mut problem, coefficients, asset.unit_size()); + let variables = VariableMap::add_to_problem(&mut problem, coefficients); // Add constraints add_constraints( @@ -153,37 +148,23 @@ pub fn perform_optimisation( ); // Solve model - let mut highs_model = problem.optimise(sense); + let mut highs_model = problem.optimise(Sense::Maximise); apply_highs_options_from_toml(&mut highs_model, &model.parameters.highs.appraisal_options) .context("Failed to apply custom HiGHS options to appraisal optimisation")?; let solution = solve_optimal(highs_model) .map_err(ModelError::into_anyhow)? .get_solution(); let solution_values = solution.columns(); + let activity = variables + .activity_vars + .keys() + .zip(solution_values.iter()) + .map(|(time_slice, &value)| (time_slice.clone(), Activity::new(value))) + .collect(); + let unmet_demand = + compute_unmet_demand(demand, &activity, commodity, asset, &model.time_slice_info); Ok(ResultsMap { - // If the asset has a defined unit size, the capacity variable represents number of units, - // otherwise it represents absolute capacity - #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] - capacity: match asset.unit_size() { - Some(unit_size) => { - AssetCapacity::Discrete(solution_values[0].round() as u32, unit_size) - } - None => AssetCapacity::Continuous(Capacity::new(solution_values[0])), - }, - // The mapping below assumes the column ordering documented on `VariableMap::add_to_problem`: - // index 0 = capacity, next `n` entries = activities (in the same key order as - // `cost_coefficients.activity_coefficients`), remaining entries = unmet demand. - activity: variables - .activity_vars - .keys() - .zip(solution_values[1..].iter()) - .map(|(time_slice, &value)| (time_slice.clone(), Activity::new(value))) - .collect(), - unmet_demand: variables - .unmet_demand_vars - .keys() - .zip(solution_values[variables.activity_vars.len() + 1..].iter()) - .map(|(time_slice, &value)| (time_slice.clone(), Flow::new(value))) - .collect(), + activity, + unmet_demand, }) } From 8fd063471206413e5a355f0d072d154de4fbb339 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 29 May 2026 15:09:41 +0100 Subject: [PATCH 5/6] Remove VariableMap --- .../investment/appraisal/optimisation.rs | 63 +++++++------------ 1 file changed, 23 insertions(+), 40 deletions(-) diff --git a/src/simulation/investment/appraisal/optimisation.rs b/src/simulation/investment/appraisal/optimisation.rs index d5cfe9e42..87b94afd3 100644 --- a/src/simulation/investment/appraisal/optimisation.rs +++ b/src/simulation/investment/appraisal/optimisation.rs @@ -21,33 +21,6 @@ use indexmap::IndexMap; /// in which columns are added to the problem when extracting solution values. pub type Variable = highs::Col; -/// Map storing variables for the optimisation problem -struct VariableMap { - /// Activity variables in each time slice - activity_vars: IndexMap, -} - -impl VariableMap { - /// Creates a new variable map by adding variables to the optimisation problem. - /// - /// # Arguments - /// * `problem` - The optimisation problem to add variables to - /// * `cost_coefficients` - Objective function coefficients for each variable - /// - /// # Returns - /// A new `VariableMap` containing all created decision variables - fn add_to_problem(problem: &mut Problem, cost_coefficients: &ObjectiveCoefficients) -> Self { - // Create activity variables for each time slice - let mut activity_vars = IndexMap::new(); - for (time_slice, cost) in &cost_coefficients.activity_coefficients { - let var = problem.add_column(cost.value(), 0.0..); - activity_vars.insert(time_slice.clone(), var); - } - - Self { activity_vars } - } -} - /// Map containing optimisation results and coefficients pub struct ResultsMap { /// Activity variables in each time slice @@ -56,30 +29,41 @@ pub struct ResultsMap { pub unmet_demand: DemandMap, } +/// Adds activity variables to the problem, one per time slice. +/// +/// Returns a map from time slice to the corresponding decision variable. +fn add_activity_vars( + problem: &mut Problem, + cost_coefficients: &ObjectiveCoefficients, +) -> IndexMap { + cost_coefficients + .activity_coefficients + .iter() + .map(|(time_slice, cost)| { + let var = problem.add_column(cost.value(), 0.0..); + (time_slice.clone(), var) + }) + .collect() +} + /// Adds constraints to the problem. fn add_constraints( problem: &mut Problem, asset: &AssetRef, max_capacity: AssetCapacity, commodity: &Commodity, - variables: &VariableMap, + activity_vars: &IndexMap, demand: &DemandMap, time_slice_info: &TimeSliceInfo, ) { - add_activity_constraints( - problem, - asset, - max_capacity, - &variables.activity_vars, - time_slice_info, - ); + add_activity_constraints(problem, asset, max_capacity, activity_vars, time_slice_info); add_demand_constraints( problem, asset, commodity, time_slice_info, demand, - &variables.activity_vars, + activity_vars, ); } @@ -134,7 +118,7 @@ pub fn perform_optimisation( ) -> Result { // Create problem and add variables let mut problem = Problem::default(); - let variables = VariableMap::add_to_problem(&mut problem, coefficients); + let activity_vars = add_activity_vars(&mut problem, coefficients); // Add constraints add_constraints( @@ -142,7 +126,7 @@ pub fn perform_optimisation( asset, max_capacity, commodity, - &variables, + &activity_vars, demand, &model.time_slice_info, ); @@ -155,8 +139,7 @@ pub fn perform_optimisation( .map_err(ModelError::into_anyhow)? .get_solution(); let solution_values = solution.columns(); - let activity = variables - .activity_vars + let activity = activity_vars .keys() .zip(solution_values.iter()) .map(|(time_slice, &value)| (time_slice.clone(), Activity::new(value))) From 7d32c2ed546292151e747096d4f06b44071fcc76 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 29 May 2026 15:57:59 +0100 Subject: [PATCH 6/6] Use `TimeSliceSelection` for demand maps --- src/fixture.rs | 4 +- src/output.rs | 40 ++++--- src/simulation/investment.rs | 101 ++++++++---------- .../investment/appraisal/constraints.rs | 16 +-- .../investment/appraisal/optimisation.rs | 42 +++----- src/time_slice.rs | 14 +++ 6 files changed, 112 insertions(+), 105 deletions(-) diff --git a/src/fixture.rs b/src/fixture.rs index e212a0c12..8eea77637 100644 --- a/src/fixture.rs +++ b/src/fixture.rs @@ -18,7 +18,7 @@ use crate::simulation::investment::appraisal::LCOXMetric; use crate::simulation::investment::appraisal::{ AppraisalOutput, coefficients::ObjectiveCoefficients, }; -use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceLevel}; +use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceLevel, TimeSliceSelection}; use crate::units::{ Activity, ActivityPerCapacity, Capacity, Dimensionless, Flow, MoneyPerActivity, MoneyPerCapacity, MoneyPerCapacityPerYear, Year, @@ -403,7 +403,7 @@ pub fn appraisal_output(asset: Asset, time_slice: TimeSliceID) -> AppraisalOutpu let activity_coefficients = indexmap! { time_slice.clone() => MoneyPerActivity(0.5) }; let market_costs = indexmap! { time_slice.clone() => MoneyPerActivity(0.4) }; let activity = indexmap! { time_slice.clone() => Activity(10.0) }; - let unmet_demand = indexmap! { time_slice.clone() => Flow(5.0) }; + let unmet_demand = indexmap! { TimeSliceSelection::Single(time_slice.clone()) => Flow(5.0) }; AppraisalOutput { asset: AssetRef::from(asset), capacity: AssetCapacity::Continuous(Capacity(42.0)), diff --git a/src/output.rs b/src/output.rs index 29d00263d..d7a1b9ad1 100644 --- a/src/output.rs +++ b/src/output.rs @@ -7,7 +7,7 @@ use crate::region::RegionID; use crate::simulation::investment::appraisal::AppraisalOutput; use crate::simulation::optimisation::{FlowMap, Solution}; use crate::simulation::prices::PriceMap; -use crate::time_slice::TimeSliceID; +use crate::time_slice::{TimeSliceID, TimeSliceLevel, TimeSliceSelection}; use crate::units::{Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerFlow}; use anyhow::{Context, Result, ensure}; use csv; @@ -271,10 +271,11 @@ struct AppraisalResultsTimeSliceRow { process_id: ProcessID, region_id: RegionID, time_slice: TimeSliceID, + time_slice_level: TimeSliceLevel, activity: Activity, activity_coefficient: MoneyPerActivity, - demand: Flow, - unmet_demand: Flow, + demand_for_selection: Flow, + unmet_demand_for_selection: Flow, } /// For writing extra debug information about the model @@ -506,13 +507,21 @@ impl DebugDataWriter { milestone_year: u32, run_description: &str, appraisal_results: &[AppraisalOutput], - demand: &IndexMap, + demand: &IndexMap, + time_slice_level: TimeSliceLevel, ) -> Result<()> { for result in appraisal_results { for (time_slice, activity) in &result.activity { let activity_coefficient = result.coefficients.activity_coefficients[time_slice]; - let demand = demand[time_slice]; - let unmet_demand = result.unmet_demand[time_slice]; + // Map the individual time slice back to its containing selection so we can look + // up selection-level demand and unmet demand. + let selection = match time_slice_level { + TimeSliceLevel::Annual => TimeSliceSelection::Annual, + TimeSliceLevel::Season => TimeSliceSelection::Season(time_slice.season.clone()), + TimeSliceLevel::DayNight => TimeSliceSelection::Single(time_slice.clone()), + }; + let demand = demand[&selection]; + let unmet_demand = result.unmet_demand[&selection]; let row = AppraisalResultsTimeSliceRow { milestone_year, run_description: self.with_context(run_description), @@ -520,10 +529,11 @@ impl DebugDataWriter { process_id: result.asset.process_id().clone(), region_id: result.asset.region_id().clone(), time_slice: time_slice.clone(), + time_slice_level, activity: *activity, activity_coefficient, - demand, - unmet_demand, + demand_for_selection: demand, + unmet_demand_for_selection: unmet_demand, }; self.appraisal_results_time_slice_writer.serialize(row)?; } @@ -606,7 +616,8 @@ impl DataWriter { milestone_year: u32, run_description: &str, appraisal_results: &[AppraisalOutput], - demand: &IndexMap, + demand: &IndexMap, + time_slice_level: TimeSliceLevel, ) -> Result<()> { if let Some(wtr) = &mut self.debug { wtr.write_appraisal_results(milestone_year, run_description, appraisal_results)?; @@ -615,6 +626,7 @@ impl DataWriter { run_description, appraisal_results, demand, + time_slice_level, )?; } @@ -756,7 +768,7 @@ mod tests { appraisal_output, asset, asset_divisible, assets, commodity_id, region_id, time_slice, }; use crate::simulation::investment::appraisal::AppraisalOutput; - use crate::time_slice::TimeSliceID; + use crate::time_slice::{TimeSliceID, TimeSliceLevel, TimeSliceSelection}; use indexmap::indexmap; use itertools::{Itertools, assert_equal}; use rstest::rstest; @@ -1211,7 +1223,7 @@ mod tests { let milestone_year = 2020; let run_description = "test_run".to_string(); let dir = tempdir().unwrap(); - let demand = indexmap! {time_slice.clone() => Flow(100.0) }; + let demand = indexmap! {TimeSliceSelection::Single(time_slice.clone()) => Flow(100.0) }; // Write appraisal time slice results { @@ -1222,6 +1234,7 @@ mod tests { &run_description, &[appraisal_output], &demand, + TimeSliceLevel::DayNight, ) .unwrap(); writer.flush().unwrap(); @@ -1235,10 +1248,11 @@ mod tests { process_id: asset.process_id().clone(), region_id: asset.region_id().clone(), time_slice: time_slice.clone(), + time_slice_level: TimeSliceLevel::DayNight, activity: Activity(10.0), activity_coefficient: MoneyPerActivity(0.5), - demand: Flow(100.0), - unmet_demand: Flow(5.0), + demand_for_selection: Flow(100.0), + unmet_demand_for_selection: Flow(5.0), }; let records: Vec = csv::Reader::from_path(dir.path().join(APPRAISAL_RESULTS_TIME_SLICE_FILE_NAME)) diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index 3f4a021f3..1fb2c3fd8 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -7,7 +7,7 @@ use crate::model::Model; use crate::output::DataWriter; use crate::region::RegionID; use crate::simulation::prices::Prices; -use crate::time_slice::{TimeSliceID, TimeSliceInfo}; +use crate::time_slice::{TimeSliceInfo, TimeSliceSelection}; use crate::units::{Capacity, Dimensionless, Flow, FlowPerCapacity}; use anyhow::{Context, Result, bail, ensure}; use indexmap::IndexMap; @@ -23,11 +23,11 @@ use appraisal::{ sort_and_filter_appraisal_outputs, }; -/// A map of demand across time slices for a specific market -type DemandMap = IndexMap; +/// A map of demand across time-slice selections for a specific market +type DemandMap = IndexMap; -/// Demand for a given combination of commodity, region and time slice -type AllDemandMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Flow>; +/// Demand for a given combination of commodity, region and time-slice selection +type AllDemandMap = IndexMap<(CommodityID, RegionID, TimeSliceSelection), Flow>; /// Represents a set of markets which are invested in together. #[derive(PartialEq, Debug, Clone, Eq, Hash)] @@ -177,8 +177,7 @@ pub fn perform_agent_investment( writer: &mut DataWriter, ) -> Result> { // Initialise net demand map - let mut net_demand = - flatten_preset_demands_for_year(&model.commodities, &model.time_slice_info, year); + let mut net_demand = flatten_preset_demands_for_year(&model.commodities, year); // Keep a list of all the assets selected // This includes Commissioned assets that are selected for retention, and new Ready assets @@ -277,7 +276,7 @@ fn select_assets_for_single_market( let demand_portion_for_market = get_demand_portion_for_market( &model.time_slice_info, demand, - commodity_id, + commodity, region_id, commodity_portion, ); @@ -452,39 +451,27 @@ fn select_assets_for_cycle( Ok(all_cycle_assets) } -/// Flatten the preset commodity demands for a given year into a map of commodity, region and -/// time slice to demand. +/// Build the initial net demand map for a given year. /// -/// Since demands for some commodities may be specified at a coarser time slice level, we need to -/// distribute these demands over all time slices. Note: the way that we do this distribution is -/// irrelevant, as demands will only be balanced to the appropriate level, but we still need to do -/// this for the solver to work. +/// Demand for each commodity is stored at its natural time-slice selection level, matching the +/// balance level at which the investment appraisal operates. /// /// **TODO**: these assumptions may need to be revisited, e.g. when we come to storage technologies -fn flatten_preset_demands_for_year( - commodities: &CommodityMap, - time_slice_info: &TimeSliceInfo, - year: u32, -) -> AllDemandMap { +fn flatten_preset_demands_for_year(commodities: &CommodityMap, year: u32) -> AllDemandMap { let mut demand_map = AllDemandMap::new(); for (commodity_id, commodity) in commodities { for ((region_id, data_year, time_slice_selection), demand) in &commodity.demand { if *data_year != year { continue; } - - // We split the demand equally over all time slices in the selection - // NOTE: since demands will only be balanced to the time slice level of the commodity - // it doesn't matter how we do this distribution, only the total matters. - #[allow(clippy::cast_precision_loss)] - let n_time_slices = time_slice_selection.iter(time_slice_info).count() as f64; - let demand_per_slice = *demand / Dimensionless(n_time_slices); - for (time_slice, _) in time_slice_selection.iter(time_slice_info) { - demand_map.insert( - (commodity_id.clone(), region_id.clone(), time_slice.clone()), - demand_per_slice, - ); - } + demand_map.insert( + ( + commodity_id.clone(), + region_id.clone(), + time_slice_selection.clone(), + ), + *demand, + ); } } demand_map @@ -502,12 +489,6 @@ fn flatten_preset_demands_for_year( fn update_net_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: &[AssetRef]) { for ((asset, commodity_id, time_slice), flow) in flows { if assets.contains(asset) { - let key = ( - commodity_id.clone(), - asset.region_id().clone(), - time_slice.clone(), - ); - // Only consider input flows and output flows from the primary output commodity // (excluding secondary outputs) if (flow < &Flow(0.0)) @@ -515,6 +496,13 @@ fn update_net_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: &[A .primary_output() .is_some_and(|p| &p.commodity.id == commodity_id) { + let level = asset + .get_flow(commodity_id) + .unwrap() + .commodity + .time_slice_level; + let selection = level.containing_selection(time_slice); + let key = (commodity_id.clone(), asset.region_id().clone(), selection); // Note: we use the negative of the flow as input flows are negative in the flow map. demand .entry(key) @@ -529,20 +517,21 @@ fn update_net_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: &[A fn get_demand_portion_for_market( time_slice_info: &TimeSliceInfo, demand: &AllDemandMap, - commodity_id: &CommodityID, + commodity: &Commodity, region_id: &RegionID, commodity_portion: Dimensionless, ) -> DemandMap { time_slice_info - .iter_ids() - .map(|time_slice| { - ( - time_slice.clone(), - commodity_portion - * *demand - .get(&(commodity_id.clone(), region_id.clone(), time_slice.clone())) - .unwrap_or(&Flow(0.0)), - ) + .iter_selections_at_level(commodity.time_slice_level) + .map(|ts_selection| { + let demand_for_selection = *demand + .get(&( + commodity.id.clone(), + region_id.clone(), + ts_selection.clone(), + )) + .unwrap_or(&Flow(0.0)); + (ts_selection, commodity_portion * demand_for_selection) }) .collect() } @@ -582,10 +571,7 @@ fn get_demand_limiting_capacity( for time_slice_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) { - let demand_for_selection: Flow = time_slice_selection - .iter(time_slice_info) - .map(|(time_slice, _)| demand[time_slice]) - .sum(); + let demand_for_selection = demand[&time_slice_selection]; // 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 @@ -804,6 +790,7 @@ fn select_best_assets( &format!("{} {} round {}", &commodity.id, &agent.id, round), &outputs_for_opts, &demand, + commodity.time_slice_level, )?; // Sort by investment priority and discard non-feasible options @@ -815,7 +802,7 @@ fn select_best_assets( let remaining_demands: Vec<_> = demand .iter() .filter(|(_, flow)| **flow > Flow(0.0)) - .map(|(time_slice, flow)| format!("{} : {:e}", time_slice, flow.value())) + .map(|(ts_selection, flow)| format!("{} : {:e}", ts_selection, flow.value())) .collect(); bail!( @@ -933,7 +920,7 @@ mod tests { ProcessInvestmentConstraint, ProcessInvestmentConstraintsMap, ProcessParameterMap, }; use crate::region::RegionID; - use crate::time_slice::{TimeSliceID, TimeSliceInfo}; + use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; use crate::units::Dimensionless; use crate::units::{ActivityPerCapacity, Capacity, Flow, FlowPerActivity, MoneyPerFlow}; use indexmap::{IndexSet, indexmap}; @@ -963,7 +950,7 @@ mod tests { let asset = asset(process); // Create demand map - demand of 10.0 for our time slice - let demand = indexmap! { time_slice.clone() => Flow(10.0)}; + let demand = indexmap! { TimeSliceSelection::Single(time_slice.clone()) => Flow(10.0)}; // Call the function let result = get_demand_limiting_capacity(&time_slice_info, &asset, &commodity_rc, &demand); @@ -1007,8 +994,8 @@ mod tests { // Create demand map with different demands for each time slice let demand = indexmap! { - time_slice1.clone() => Flow(4.0), // Requires capacity of 4.0/0.2 = 20.0 - time_slice2.clone() => Flow(3.0), // Would require infinite capacity, but should be skipped + TimeSliceSelection::Single(time_slice1.clone()) => Flow(4.0), // Requires capacity of 4.0/0.2 = 20.0 + TimeSliceSelection::Single(time_slice2.clone()) => Flow(3.0), // Would require infinite capacity, but should be skipped }; // Call the function diff --git a/src/simulation/investment/appraisal/constraints.rs b/src/simulation/investment/appraisal/constraints.rs index 01270ace0..74beeaeb7 100644 --- a/src/simulation/investment/appraisal/constraints.rs +++ b/src/simulation/investment/appraisal/constraints.rs @@ -4,7 +4,6 @@ use super::optimisation::Variable; use crate::asset::{AssetCapacity, AssetRef}; use crate::commodity::Commodity; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::Flow; use highs::RowProblem as Problem; use indexmap::IndexMap; @@ -54,13 +53,14 @@ pub fn add_demand_constraints( activity_vars: &IndexMap, ) { for ts_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) { - let mut demand_for_ts_selection = Flow(0.0); - let mut terms = Vec::new(); - for (time_slice, _) in ts_selection.iter(time_slice_info) { - demand_for_ts_selection += demand[time_slice]; - let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff; - terms.push((activity_vars[time_slice], flow_coeff.value())); - } + let demand_for_ts_selection = demand[&ts_selection]; + let terms: Vec<_> = ts_selection + .iter(time_slice_info) + .map(|(time_slice, _)| { + let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff; + (activity_vars[time_slice], flow_coeff.value()) + }) + .collect(); problem.add_row(0.0..=demand_for_ts_selection.value(), terms); } } diff --git a/src/simulation/investment/appraisal/optimisation.rs b/src/simulation/investment/appraisal/optimisation.rs index 87b94afd3..150ea8d8b 100644 --- a/src/simulation/investment/appraisal/optimisation.rs +++ b/src/simulation/investment/appraisal/optimisation.rs @@ -10,7 +10,7 @@ use crate::simulation::optimisation::ModelError; use crate::simulation::optimisation::apply_highs_options_from_toml; use crate::simulation::optimisation::solve_optimal; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::{Activity, Dimensionless, Flow}; +use crate::units::{Activity, Flow}; use anyhow::{Context, Result}; use highs::{RowProblem as Problem, Sense}; use indexmap::IndexMap; @@ -67,15 +67,12 @@ fn add_constraints( ); } -/// Computes remaining unmet demand per time slice after a solve. +/// Computes remaining unmet demand per time-slice selection after a solve. /// /// For each time-slice selection at the commodity's balance level, the selection-level residual -/// (`demand_total - supply_total`, clamped to zero) is divided equally across the time slices in -/// the selection. -/// -/// The exact per-time-slice distribution is arbitrary: all downstream consumers sum values back up -/// to the selection level before using them (e.g. the next round's demand constraint, and -/// `is_any_remaining_demand`), so only the selection-level total needs to be correct. +/// (`demand - supply`, clamped to zero) is stored directly. Downstream consumers operate at the +/// selection level (e.g. the next round's demand constraint and `is_any_remaining_demand`), so +/// there is no need to distribute values across individual time slices. fn compute_unmet_demand( demand: &DemandMap, activity: &IndexMap, @@ -83,24 +80,19 @@ fn compute_unmet_demand( asset: &AssetRef, time_slice_info: &TimeSliceInfo, ) -> DemandMap { - let mut unmet_demand = DemandMap::new(); let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff; - for ts_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) { - let time_slices: Vec<_> = ts_selection.iter(time_slice_info).collect(); - let demand_for_selection: Flow = time_slices.iter().map(|(ts, _)| demand[*ts]).sum(); - let supply_for_selection: Flow = time_slices - .iter() - .map(|(ts, _)| activity[*ts] * flow_coeff) - .sum(); - - #[allow(clippy::cast_precision_loss)] - let unmet_per_slice = (demand_for_selection - supply_for_selection).max(Flow(0.0)) - / Dimensionless(time_slices.len() as f64); - for (time_slice, _) in &time_slices { - unmet_demand.insert((*time_slice).clone(), unmet_per_slice); - } - } - unmet_demand + time_slice_info + .iter_selections_at_level(commodity.time_slice_level) + .map(|ts_selection| { + let demand_for_selection = demand[&ts_selection]; + let supply_for_selection: Flow = ts_selection + .iter(time_slice_info) + .map(|(ts, _)| activity[ts] * flow_coeff) + .sum(); + let unmet = (demand_for_selection - supply_for_selection).max(Flow(0.0)); + (ts_selection, unmet) + }) + .collect() } /// Performs optimisation for an asset, given the coefficients and demand. diff --git a/src/time_slice.rs b/src/time_slice.rs index aaefabe7d..aca8df075 100644 --- a/src/time_slice.rs +++ b/src/time_slice.rs @@ -217,6 +217,20 @@ impl TimeSliceLevel { } } +impl Serialize for TimeSliceLevel { + fn serialize(&self, serializer: S) -> std::result::Result + where + S: serde::Serializer, + { + let s = match self { + Self::DayNight => "daynight", + Self::Season => "season", + Self::Annual => "annual", + }; + serializer.serialize_str(s) + } +} + /// Information about the time slices in the simulation, including names and durations #[derive(PartialEq, Debug)] pub struct TimeSliceInfo {