From 46e7a401dd59ba75588cdd0b882e6e821af00602 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 12:12:56 +0100 Subject: [PATCH 01/13] MILP modification to ensure consistency in second pass --- src/graph/investment.rs | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index 5f30ac1ec..44ba78b19 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -272,7 +272,8 @@ fn order_sccs( // Record whether any edge inside the original SCC goes from market i to market j; these become penalties. let mut penalties = vec![vec![0.0f64; n]; n]; - let mut has_external_outgoing = vec![false; n]; + let mut has_external_outgoing: Vec = vec![false; n]; + let mut has_external_incoming: Vec = vec![false; n]; for (i, &idx) in original_indices.iter().enumerate() { // Loop over the edges going out of this node for edge in original_graph.edges_directed(idx, Direction::Outgoing) { @@ -285,6 +286,13 @@ fn order_sccs( has_external_outgoing[i] = true; } } + + // Check whether this node has any incoming edges from outside the SCC + for edge in original_graph.edges_directed(idx, Direction::Incoming) { + if !index_position.contains_key(&edge.source()) { + has_external_incoming[i] = true; + } + } } // Bias: if market j has outgoing edges to nodes outside this SCC, we prefer to place it earlier. @@ -345,6 +353,30 @@ fn order_sccs( } } + // Every SCC node must retain at least one correctly-ordered incoming edge. + for j in 0..n { + let mut incoming_terms = Vec::new(); + for i in 0..n { + if i == j { + continue; + } + + // We need to know whether the original graph contains i -> j. + // If so, x[i][j] represents that edge being retained. + if original_graph + .find_edge(original_indices[i], original_indices[j]) + .is_some() + { + incoming_terms.push((vars[i][j].unwrap(), 1.0)); + } + } + + // If the node has an incoming edge from outside the SCC, then it doesn't need a + // correctly-ordered internal incoming edge. Otherwise it does. + let required = if has_external_incoming[j] { 0.0 } else { 1.0 }; + problem.add_row(required.., incoming_terms); + } + let model = problem.optimise(Sense::Minimise); let solved = match model.try_solve() { Ok(solved) => solved, From 683c8661ff71fd1f092dea0a85e4db10ea6bdbb7 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Wed, 9 Sep 2026 10:04:25 +0100 Subject: [PATCH 02/13] Remove flexible capacity balancing for circularities --- schemas/input/model.yaml | 4 - src/model/parameters.rs | 10 -- src/simulation/market.rs | 78 +----------- src/simulation/optimisation.rs | 138 +-------------------- src/simulation/optimisation/constraints.rs | 81 +++--------- 5 files changed, 25 insertions(+), 286 deletions(-) diff --git a/schemas/input/model.yaml b/schemas/input/model.yaml index 5db20013a..fc8fb347c 100644 --- a/schemas/input/model.yaml +++ b/schemas/input/model.yaml @@ -95,10 +95,6 @@ properties: type: integer description: Number of iterations to perform when calculating prices for cyclically-dependent markets default: 1 - capacity_margin: - type: number - description: Slack proportion for assets selected during cycle balancing to absorb small demand shifts - default: 0.2 mothball_years: type: integer default: 0 diff --git a/src/model/parameters.rs b/src/model/parameters.rs index 6a50ec1e8..b425fff09 100644 --- a/src/model/parameters.rs +++ b/src/model/parameters.rs @@ -115,13 +115,6 @@ pub struct ModelParameters { pub price_tolerance: Dimensionless, /// Number of iterations to perform when calculating prices for cyclically-dependent markets. pub price_cycle_iterations: u32, - /// Slack applied during cycle balancing, allowing newly selected assets to flex their capacity - /// by this proportion. - /// - /// Existing assets remain fixed; this gives newly selected assets the wiggle-room to absorb - /// small demand changes before we would otherwise need to break for re-investment. - #[serde(deserialize_with = "deserialise_finite_non_negative")] - pub capacity_margin: Dimensionless, /// Number of years an asset can remain unused before being decommissioned pub mothball_years: u32, /// Absolute tolerance when checking if remaining demand is close enough to zero @@ -155,7 +148,6 @@ impl Default for ModelParameters { max_ironing_out_iterations: 1, price_tolerance: Dimensionless(1e-6), price_cycle_iterations: 1, - capacity_margin: Dimensionless(0.2), mothball_years: 0, remaining_demand_absolute_tolerance: DEFAULT_REMAINING_DEMAND_ABSOLUTE_TOLERANCE, highs: HighsOptions::default(), @@ -368,8 +360,6 @@ impl ModelParameters { // price_cycle_iterations check_price_cycle_iterations(self.price_cycle_iterations)?; - // capacity_margin already validated with deserialise_finite_non_negative - // remaining_demand_absolute_tolerance already validated with // deserialise_finite_non_negative; check remaining constraints here check_remaining_demand_absolute_tolerance( diff --git a/src/simulation/market.rs b/src/simulation/market.rs index 717a2cf77..46852e23b 100644 --- a/src/simulation/market.rs +++ b/src/simulation/market.rs @@ -1,7 +1,7 @@ //! Code for creating sets of markets. use super::optimisation::DispatchRun; use crate::agent::Agent; -use crate::asset::{Asset, AssetCapacity, AssetIterator, AssetRef, AssetState}; +use crate::asset::{Asset, AssetCapacity, AssetIterator, AssetRef}; use crate::commodity::{Commodity, CommodityID}; use crate::model::Model; use crate::output::DataWriter; @@ -236,17 +236,9 @@ pub fn select_assets_for_single_market( /// Iterates through the a pre-ordered set of markets forming a cycle, selecting assets for each /// market in turn. /// -/// Dispatch optimisation is performed after each market is visited to rebalance demand. -/// While dispatching, newly selected (`Ready`) assets are given flexible capacity (bounded by -/// `capacity_margin`) so small demand shifts caused by later markets can be absorbed. After all -/// markets have been visited once, the final set of assets is returned, applying any capacity -/// adjustments from the final full-system dispatch optimisation. +/// Dispatch optimisation is performed after each market is visited. /// -/// Dispatch may fail at any point if new demands are encountered for previously visited markets, -/// and the `capacity_margin` is not sufficient to absorb the demand shift. At this point, the -/// simulation is terminated with an error prompting the user to increase the `capacity_margin`. -/// A longer-term solution (TODO) may be to trigger re-investment for the affected markets. Other -/// yet-to-implement features may also help to stabilise the cycle, such as capacity growth limits. +/// Dispatch may fail at any point if new demands are encountered for previously visited markets. #[allow(clippy::too_many_arguments)] pub fn select_assets_for_cycle( model: &Model, @@ -265,7 +257,6 @@ pub fn select_assets_for_cycle( // Iterate over the markets to select assets let mut current_demand = demand.clone(); let mut assets_for_cycle = IndexMap::new(); - let mut last_solution = None; for (idx, (commodity_id, region_id)) in markets.iter().enumerate() { // Select assets for this market let assets = select_assets_for_single_market( @@ -292,53 +283,15 @@ pub fn select_assets_for_cycle( let mut markets_to_balance = seen_markets.to_vec(); markets_to_balance.extend_from_slice(&markets[0..=idx]); - // We allow all `Ready` state assets to have flexible capacity - let flexible_capacity_assets: Vec<_> = assets_for_cycle_flat - .iter() - .filter(|asset| matches!(asset.state(), AssetState::Ready { .. })) - .cloned() - .collect(); - - // Retrieve installable capacity limits for flexible capacity assets. - let mut agent_share_cache = HashMap::new(); - let capacity_limits = flexible_capacity_assets - .iter() - .filter_map(|asset| { - let agent_id = asset.agent_id().unwrap(); - let commodity_id = asset.primary_output_commodity().unwrap(); - let agent_share = *agent_share_cache - .entry((agent_id, commodity_id)) - .or_insert_with(|| { - model.agents[agent_id].commodity_portions[&(commodity_id.clone(), year)] - }); - asset - .process() - .agent_addition_limit(asset.region_id(), asset.commission_year(), agent_share) - .map(|max_capacity| (asset.clone(), max_capacity)) - }) - .collect::>(); - // Run dispatch let solution = DispatchRun::new(model, &all_assets, year) .without_commodity_constraints() .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, - ) .run( &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment"), writer, ) - .with_context(|| { - format!( - "Cycle balancing failed for cycle ({markets_str}), capacity_margin: {}. \ - Try increasing the capacity_margin.", - model.parameters.capacity_margin - ) - })?; + .with_context(|| format!("Dispatch failed for cycle ({markets_str})"))?; // Calculate new net demand map with all assets selected so far current_demand.clone_from(demand); @@ -347,29 +300,10 @@ pub fn select_assets_for_cycle( &solution.create_flow_map(), &assets_for_cycle_flat, ); - last_solution = Some(solution); } - // Finally, update flexible capacity assets based on the final solution - let mut all_cycle_assets: Vec<_> = assets_for_cycle.into_values().flatten().collect(); - if let Some(solution) = last_solution { - let new_capacities: HashMap<_, _> = solution.iter_capacity().collect(); - for asset in &mut all_cycle_assets { - if let Some(new_capacity) = new_capacities.get(asset) { - debug!( - "Capacity of asset '{}' modified during cycle balancing ({} to {})", - asset.process_id(), - asset.total_capacity(), - new_capacity.total_capacity() - ); - asset.make_mut().set_capacity(*new_capacity); - } - } - } - - // Drop any assets who's capacities were dropped to zero - all_cycle_assets.retain(|asset| asset.num_tranches() > 0); - + // Collect assets + let all_cycle_assets: Vec<_> = assets_for_cycle.into_values().flatten().collect(); Ok(all_cycle_assets) } diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 3b42e591e..7c464b085 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -1,25 +1,20 @@ //! Code for performing dispatch optimisation. //! //! This is used to calculate commodity flows and prices. -use crate::asset::{Asset, AssetCapacity, AssetRef, AssetState}; +use crate::asset::{Asset, AssetRef}; use crate::commodity::CommodityID; -use crate::finance::annual_capital_cost; use crate::input::format_items_with_cap; use crate::model::Model; use crate::output::DataWriter; use crate::region::RegionID; use crate::simulation::PriceMap; use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; -use crate::units::{ - Activity, Capacity, Dimensionless, Flow, Money, MoneyPerActivity, MoneyPerCapacity, - MoneyPerFlow, Year, -}; +use crate::units::{Activity, Flow, Money, MoneyPerActivity, MoneyPerFlow}; use anyhow::{Context, Result, anyhow, bail}; use highs::{HighsModelStatus, RowProblem as Problem, Sense}; use indexmap::{IndexMap, IndexSet}; use itertools::{chain, iproduct}; use log::warn; -use std::collections::HashMap; use std::error::Error; use std::ops::Range; @@ -38,9 +33,6 @@ type Variable = highs::Col; /// The map of activity variables for assets type ActivityVariableMap = IndexMap<(AssetRef, TimeSliceID), Variable>; -/// A map of capacity variables for assets -type CapacityVariableMap = IndexMap; - /// Variables representing unmet demand for a given market type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Variable>; @@ -57,8 +49,6 @@ pub struct VariableMap { activity_vars: ActivityVariableMap, existing_asset_var_idx: Range, candidate_asset_var_idx: Range, - capacity_vars: CapacityVariableMap, - capacity_var_idx: Range, unmet_demand_vars: UnmetDemandVariableMap, unmet_demand_var_idx: Range, } @@ -104,8 +94,6 @@ impl VariableMap { activity_vars, existing_asset_var_idx, candidate_asset_var_idx, - capacity_vars: CapacityVariableMap::new(), - capacity_var_idx: Range::default(), unmet_demand_vars: UnmetDemandVariableMap::default(), unmet_demand_var_idx: Range::default(), } @@ -173,11 +161,6 @@ impl VariableMap { fn activity_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> { self.activity_vars.keys() } - - /// Iterate over capacity variables - fn iter_capacity_vars(&self) -> impl Iterator { - self.capacity_vars.iter().map(|(asset, var)| (asset, *var)) - } } /// The solution to the dispatch optimisation problem @@ -261,22 +244,6 @@ impl Solution<'_> { }) } - /// Iterate over capacity values - pub fn iter_capacity(&self) -> impl Iterator { - self.variables - .capacity_vars - .keys() - .zip(self.solution.columns()[self.variables.capacity_var_idx.clone()].iter()) - .map(|(asset, capacity_var)| { - // The capacity variable represents number of tranches - let tranche_size = asset.capacity().tranche_size(); - #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] - let asset_capacity = AssetCapacity::new(capacity_var.round() as u32, tranche_size); - - (asset, asset_capacity) - }) - } - /// Keys and dual values for commodity balance constraints. pub fn iter_commodity_balance_duals( &self, @@ -296,11 +263,7 @@ impl Solution<'_> { /// Keys and dual values for activity constraints. /// - /// Note: if there are any flexible capacity assets, these will have two duals with identical - /// keys, and there will be no way to distinguish between them in the resulting iterator. - /// Recommended for now only to use this function when there are no flexible capacity assets. - /// - /// Also note: this excludes seasonal and annual constraints. Recommended for now not to use + /// Note: this excludes seasonal and annual constraints. Recommended for now not to use /// this for models that include seasonal or annual availability constraints. pub fn iter_activity_duals( &self, @@ -436,14 +399,11 @@ fn filter_input_prices( 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 PriceMap>, include_commodity_constraints: bool, year: u32, - capacity_margin: Dimensionless, } impl<'model, 'run> DispatchRun<'model, 'run> { @@ -452,29 +412,11 @@ impl<'model, 'run> DispatchRun<'model, 'run> { Self { model, existing_assets: assets, - flexible_capacity_assets: &[], - capacity_limits: None, candidate_assets: &[], markets_to_balance: &[], input_prices: None, include_commodity_constraints: true, year, - capacity_margin: Dimensionless(0.0), - } - } - - /// Include the specified flexible capacity assets in the dispatch run - pub fn with_flexible_capacity_assets( - self, - flexible_capacity_assets: &'run [AssetRef], - capacity_limits: Option<&'run HashMap>, - capacity_margin: Dimensionless, - ) -> Self { - Self { - flexible_capacity_assets, - capacity_limits, - capacity_margin, - ..self } } @@ -741,25 +683,6 @@ impl<'model, 'run> DispatchRun<'model, 'run> { variables.add_unmet_demand_variables(&mut problem, self.model, markets_to_balance); } - // Check flexible capacity assets is a subset of existing assets - for asset in self.flexible_capacity_assets { - assert!( - self.existing_assets.contains(asset), - "Flexible capacity assets must be a subset of existing assets. Offending asset: {asset:?}" - ); - } - - // Add capacity variables for flexible capacity assets - if !self.flexible_capacity_assets.is_empty() { - variables.capacity_var_idx = add_capacity_variables( - &mut problem, - &mut variables.capacity_vars, - self.flexible_capacity_assets, - self.capacity_limits, - self.capacity_margin, - ); - } - // Add constraints let all_assets = chain(self.existing_assets.iter(), self.candidate_assets.iter()); let constraint_keys = add_model_constraints( @@ -821,51 +744,6 @@ fn add_activity_variables( start..problem.num_cols() } -fn add_capacity_variables( - problem: &mut Problem, - variables: &mut CapacityVariableMap, - assets: &[AssetRef], - capacity_limits: Option<&HashMap>, - capacity_margin: Dimensionless, -) -> Range { - let capacity_margin = capacity_margin.value(); - - // This line **must** come before we add more variables - let start = problem.num_cols(); - - for asset in assets { - // Can only have flexible capacity for `Ready` assets - assert!( - matches!(asset.state(), AssetState::Ready { .. }), - "Flexible capacity can only be assigned to `Ready` type assets. Offending asset: {asset:?}" - ); - - // Coefficient: cost per capacity - let coeff = calculate_capacity_coefficient(asset); - - // 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. - // Since capacity variables are numbers of tranches, we apply constraints to the tranche count - let tranche_size = asset.capacity().tranche_size(); - let current_tranches = asset.capacity().num_tranches(); - - let lower = (current_tranches as f64 * (1.0 - capacity_margin)).max(0.0); - - let mut upper = current_tranches as f64 * (1.0 + capacity_margin); - if let Some(limit) = capacity_limits.and_then(|limits| limits.get(asset)) { - upper = upper.min((*limit / tranche_size).value()); - } - - let var = problem.add_integer_column((coeff * tranche_size).value(), lower..=upper); - - let existing = variables.insert(asset.clone(), var).is_some(); - assert!(!existing, "Duplicate entry for var"); - } - - start..problem.num_cols() -} - /// Calculate the cost coefficient for an activity variable. /// /// Normally, the cost coefficient is the same as the asset's operating costs for the given year and @@ -895,13 +773,3 @@ fn calculate_activity_coefficient( opex } } - -/// Calculate the cost coefficient for a capacity variable (for flexible capacity assets only). -/// -/// This includes both the annual fixed operating cost and the annual capital cost. -fn calculate_capacity_coefficient(asset: &AssetRef) -> MoneyPerCapacity { - let param = asset.process_parameter(); - let annual_fixed_operating_cost = param.fixed_operating_cost * Year(1.0); - annual_fixed_operating_cost - + annual_capital_cost(param.capital_cost, param.lifetime, param.discount_rate) -} diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs index fd233408a..6b5172bd2 100644 --- a/src/simulation/optimisation/constraints.rs +++ b/src/simulation/optimisation/constraints.rs @@ -9,7 +9,7 @@ use crate::time_slice::{Season, TimeSliceInfo, TimeSliceSelection}; use crate::units::{Flow, MoneyPerCapacityPerYear, UnitType, Year}; use highs::RowProblem as Problem; use indexmap::IndexMap; -use std::collections::{HashMap, HashSet}; +use std::collections::HashMap; /// Corresponding variables for a constraint along with the row offset in the solution pub struct KeysWithOffset { @@ -456,10 +456,7 @@ fn candidate_balance_epsilon( /// /// Returns an `ActivityKeys` where `offset` is the row index of the first /// activity constraint added and `keys` enumerates the `(asset, time_selection)` -/// entries in the same row order. Note that for flexible-capacity assets two rows -/// (upper and lower bounds) are added per selection; in that case the same key is -/// stored twice to match the solver ordering. -/// +/// entries in the same row order. #[doc = concat!("[1]: ", crate::docs_url!("model/dispatch_optimisation.html#asset-activity-limits"))] fn add_activity_constraints<'a, I>( problem: &mut Problem, @@ -476,61 +473,23 @@ where let offset = problem.num_rows(); let mut keys = Vec::new(); - let capacity_vars: IndexMap<&AssetRef, highs::Col> = variables.iter_capacity_vars().collect(); // Create constraints for each asset for asset in assets { - if let Some(&capacity_var) = capacity_vars.get(asset) { - // Asset with flexible capacity - 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(); - - // The capacity variable represents number of tranches, so we need to multiply the - // per-capacity limits by the tranche size. - let tranche_size = asset.capacity().tranche_size(); - upper_limit *= tranche_size.value(); - lower_limit *= tranche_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 = variables.get_activity_var(asset, time_slice); - 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); + for (ts_selection, limits) in asset.iter_activity_limits() { + let limits = limits.start().value()..=limits.end().value(); - // Store keys for retrieving duals later. - // TODO: a bit of a hack pushing identical keys twice. Safe for now so long as we don't - // use the activity duals for anything important when using flexible capacity assets. - keys.push((asset.clone(), ts_selection.clone())); - keys.push((asset.clone(), ts_selection.clone())); - } - } else { - // Fixed-capacity asset: simple absolute activity limits. - for (ts_selection, limits) in asset.iter_activity_limits() { - let limits = limits.start().value()..=limits.end().value(); - - // Collect activity terms for the time slices in this selection - let terms = ts_selection - .iter(time_slice_info) - .map(|(time_slice, _)| (variables.get_activity_var(asset, time_slice), 1.0)) - .collect::>(); + // Collect activity terms for the time slices in this selection + let terms = ts_selection + .iter(time_slice_info) + .map(|(time_slice, _)| (variables.get_activity_var(asset, time_slice), 1.0)) + .collect::>(); - // Constraint: sum of activities in selection within limits - problem.add_row(limits, &terms); + // Constraint: sum of activities in selection within limits + problem.add_row(limits, &terms); - // Store keys for retrieving duals later. - keys.push((asset.clone(), ts_selection.clone())); - } + // Store keys for retrieving duals later. + keys.push((asset.clone(), ts_selection.clone())); } } @@ -544,7 +503,7 @@ where /// which is the authoritative check for equivalence. This also handles hash collisions correctly. /// /// The caller must ensure that `assets` contains only assets eligible for equal-utilisation -/// constraints (i.e. flexible-capacity assets have already been filtered out). +/// constraints. fn group_dispatch_equivalent_assets<'a, I>(assets: I) -> Vec> where I: Iterator, @@ -581,8 +540,7 @@ where /// Add constraints requiring dispatch-equivalent assets to have equal utilisation in each time /// slice. /// -/// Flexible-capacity assets are excluded because their maximum activity depends on a decision -/// variable. The constraints added here are not included in [`ConstraintKeys`], as their duals +/// The constraints added here are not included in [`ConstraintKeys`], as their duals /// are not currently used. fn add_equal_utilisation_constraints<'a, I>( problem: &mut Problem, @@ -592,14 +550,7 @@ fn add_equal_utilisation_constraints<'a, I>( ) where I: Iterator + 'a, { - // Identify flexible-capacity assets so we can exclude them from the constraints - let flexible_assets: HashSet<_> = variables - .iter_capacity_vars() - .map(|(asset, _)| asset) - .collect(); - - let asset_groups = - group_dispatch_equivalent_assets(assets.filter(|asset| !flexible_assets.contains(asset))); + let asset_groups = group_dispatch_equivalent_assets(assets); // For each group of assets, add constraints to force equal utilisation in each time slice // This is done by anchoring each asset to the first asset in the group (-> (n-1) constraints From 8931ef3182f5816ab6d911ea76f1468b1c562e30 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Wed, 9 Sep 2026 10:26:18 +0100 Subject: [PATCH 03/13] Fix doc build --- src/simulation/optimisation/constraints.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs index 6b5172bd2..17d89ed96 100644 --- a/src/simulation/optimisation/constraints.rs +++ b/src/simulation/optimisation/constraints.rs @@ -457,6 +457,7 @@ fn candidate_balance_epsilon( /// Returns an `ActivityKeys` where `offset` is the row index of the first /// activity constraint added and `keys` enumerates the `(asset, time_selection)` /// entries in the same row order. +/// #[doc = concat!("[1]: ", crate::docs_url!("model/dispatch_optimisation.html#asset-activity-limits"))] fn add_activity_constraints<'a, I>( problem: &mut Problem, From 87033cb02ebc15d88dd3eb76fd2d700135fc14bd Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 14:35:26 +0100 Subject: [PATCH 04/13] Two-pass algorithm for solving circularities --- benches/assets.rs | 12 +- src/simulation.rs | 10 +- src/simulation/investment.rs | 11 +- src/simulation/market.rs | 193 +++++++++++++++++---- src/simulation/optimisation.rs | 11 +- src/simulation/optimisation/constraints.rs | 24 ++- 6 files changed, 203 insertions(+), 58 deletions(-) diff --git a/benches/assets.rs b/benches/assets.rs index 9a2b0ba02..ba6414777 100644 --- a/benches/assets.rs +++ b/benches/assets.rs @@ -95,13 +95,15 @@ fn calculate_seed_prices( candidates: &[AssetRef], writer: &mut DataWriter, ) -> Prices { - let solution_existing = DispatchRun::new(model, base_year_assets, BASE_YEAR) + let market_demands = collect_preset_demands_for_year(&model.commodities, BASE_YEAR); + let solution_existing = DispatchRun::new(model, base_year_assets, BASE_YEAR, &market_demands) .run("bench setup: without candidates", writer) .expect("Dispatch without candidates failed"); - let solution_with_candidates = DispatchRun::new(model, base_year_assets, BASE_YEAR) - .with_candidates(candidates) - .run("bench setup: with candidates", writer) - .expect("Dispatch with candidates failed"); + let solution_with_candidates = + DispatchRun::new(model, base_year_assets, BASE_YEAR, &market_demands) + .with_candidates(candidates) + .run("bench setup: with candidates", writer) + .expect("Dispatch with candidates failed"); calculate_prices( model, diff --git a/src/simulation.rs b/src/simulation.rs index 914a8a328..cc65a2221 100644 --- a/src/simulation.rs +++ b/src/simulation.rs @@ -15,7 +15,7 @@ use std::sync::Arc; pub mod optimisation; use optimisation::{DispatchRun, FlowMap}; pub mod investment; -use investment::perform_agent_investment; +use investment::{collect_preset_demands_for_year, perform_agent_investment}; pub mod market; pub mod prices; pub use prices::PriceMap; @@ -178,13 +178,15 @@ fn run_dispatch_for_year( debug_assert!(assets.iter().all(|asset| !asset.is_candidate())); debug_assert!(candidates.iter().all(|asset| asset.is_candidate())); + let market_demands = collect_preset_demands_for_year(&model.commodities, year); + // Run dispatch optimisation with existing assets only, if there are any. If not, then assume no // flows (i.e. all are zero) let (solution_existing, flow_map) = if assets.is_empty() { (None, FlowMap::default()) } else { - let solution = - DispatchRun::new(model, assets, year).run("final without candidates", writer)?; + let solution = DispatchRun::new(model, assets, year, &market_demands) + .run("final without candidates", writer)?; let flow_map = solution.create_flow_map(); (Some(solution), flow_map) }; @@ -195,7 +197,7 @@ fn run_dispatch_for_year( None } else { Some( - DispatchRun::new(model, assets, year) + DispatchRun::new(model, assets, year, &market_demands) .with_candidates(candidates) .run("final with candidates", writer)?, ) diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index 1340c402a..f8175ffe1 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -56,7 +56,8 @@ pub fn perform_agent_investment( writer: &mut DataWriter, ) -> Result> { // Initialise net demand map - let mut net_demand = collect_preset_demands_for_year(&model.commodities, year); + let preset_demands = collect_preset_demands_for_year(&model.commodities, year); + let mut net_demand = preset_demands.clone(); // Keep a list of all the assets selected // This includes Commissioned assets that are selected for retention, and new Ready assets @@ -110,7 +111,7 @@ pub fn perform_agent_investment( // As upstream markets by definition will not yet have producers, we explicitly set // their prices using external values so that they don't appear free - let solution = DispatchRun::new(model, &all_selected_assets, year) + let solution = DispatchRun::new(model, &all_selected_assets, year, &preset_demands) .without_commodity_constraints() .with_market_balance_subset(&seen_markets) .with_input_prices(&prices.shadow) @@ -179,10 +180,14 @@ pub fn update_net_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: 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 + let value = demand .entry(key) .and_modify(|value| *value -= *flow) .or_insert(-*flow); + + if *value < Flow(0.0) { + *value = Flow(0.0); + } } } } diff --git a/src/simulation/market.rs b/src/simulation/market.rs index 46852e23b..7f5c019b1 100644 --- a/src/simulation/market.rs +++ b/src/simulation/market.rs @@ -1,11 +1,11 @@ //! Code for creating sets of markets. -use super::optimisation::DispatchRun; +use super::optimisation::{DispatchRun, FlowMap}; use crate::agent::Agent; use crate::asset::{Asset, AssetCapacity, AssetIterator, AssetRef}; use crate::commodity::{Commodity, CommodityID}; use crate::model::Model; use crate::output::DataWriter; -use crate::process::{Process, ProcessID}; +use crate::process::{FlowDirection, Process, ProcessID}; use crate::region::RegionID; use crate::simulation::investment::{ AllDemandMap, DemandMap, calculate_candidate_asset_capacity_scale, select_best_assets, @@ -15,7 +15,6 @@ use crate::simulation::prices::Prices; use crate::time_slice::TimeSliceInfo; use crate::units::{Capacity, Dimensionless, Flow}; use anyhow::{Context, Result}; -use indexmap::IndexMap; use itertools::{Itertools, chain}; use log::debug; use std::collections::HashMap; @@ -80,6 +79,7 @@ impl MarketSet { demand, existing_assets, prices, + &[], writer, ), MarketSet::Cycle(markets) => { @@ -159,6 +159,7 @@ pub fn select_assets_for_single_market( demand: &AllDemandMap, existing_assets: &[AssetRef], prices: &Prices, + excluded_processes: &[ProcessID], writer: &mut DataWriter, ) -> Result> { let commodity = &model.commodities[commodity_id]; @@ -182,7 +183,7 @@ pub fn select_assets_for_single_market( ); // Existing and candidate assets from which to choose - let opt_assets = get_asset_options( + let mut opt_assets = get_asset_options( existing_assets, &demand_portion_for_market, agent, @@ -193,6 +194,9 @@ pub fn select_assets_for_single_market( ) .collect::>(); + // Exclude certain processes from the options + opt_assets.retain(|asset| !excluded_processes.contains(asset.process_id())); + // Calculate the agent's share of addition limits for candidate processes let agent_addition_limits = collect_agent_limits( agent, @@ -233,13 +237,13 @@ pub fn select_assets_for_single_market( Ok(selected_assets) } -/// Iterates through the a pre-ordered set of markets forming a cycle, selecting assets for each +/// Iterates through a pre-ordered set of markets forming a cycle, selecting assets for each /// market in turn. /// /// Dispatch optimisation is performed after each market is visited. /// /// Dispatch may fail at any point if new demands are encountered for previously visited markets. -#[allow(clippy::too_many_arguments)] +#[allow(clippy::too_many_arguments, clippy::too_many_lines)] pub fn select_assets_for_cycle( model: &Model, markets: &[(CommodityID, RegionID)], @@ -247,64 +251,181 @@ pub fn select_assets_for_cycle( demand: &AllDemandMap, existing_assets: &[AssetRef], prices: &Prices, - seen_markets: &[(CommodityID, RegionID)], - previously_selected_assets: &[AssetRef], + _seen_markets: &[(CommodityID, RegionID)], + _previously_selected_assets: &[AssetRef], writer: &mut DataWriter, ) -> Result> { // Precompute a joined string for logging let markets_str = markets.iter().map(|(c, r)| format!("{c}|{r}")).join(", "); - // Iterate over the markets to select assets - let mut current_demand = demand.clone(); - let mut assets_for_cycle = IndexMap::new(); - for (idx, (commodity_id, region_id)) in markets.iter().enumerate() { + // Collect feedback process/region pairs. These are process-region combinations that are the + // wrong way around in the investment order (i.e. any incoming commodity appears before the + // primary output commodity in that region) + let market_order: HashMap<_, _> = markets.iter().enumerate().map(|(i, m)| (m, i)).collect(); + let mut feedback_processes: Vec<(ProcessID, RegionID)> = Vec::new(); + for (process_id, process) in &model.processes { + let Some(primary_output) = &process.primary_output else { + continue; + }; + + // Iterate over the regions the process operates in and assess feedback status. + // For now, we assume that all incoming/outgoing flows occur in the region that the + // process operates in - this will need to be revisited when we implement trade + for region in &process.regions { + let primary_market = (primary_output.clone(), region.clone()); + let Some(&primary_order) = market_order.get(&primary_market) else { + continue; + }; + + // For any (input, primary output) pair, mark (process, region) as feedback if both + // input and primary output are part of the SCC (in the market_order), and input + // comes before primary output in the market order + let flows = &process.flows[&(region.clone(), year)]; + let is_feedback = flows.iter().any(|(commodity_id, flow)| { + flow.direction() == FlowDirection::Input + && market_order + .get(&(commodity_id.clone(), region.clone())) + .is_some_and(|&input_order| input_order < primary_order) + }); + if is_feedback { + feedback_processes.push((process_id.clone(), region.clone())); + } + } + } + + // STEP 1 + // Iterate over the markets in order1, considering all processes + let mut net_demand = demand.clone(); + let mut assets_for_first_pass = Vec::new(); + let mut retained_first_pass_assets = Vec::new(); + let mut excluded_flows = FlowMap::new(); + for market in markets { + let (commodity_id, region_id) = market.clone(); + // Select assets for this market - let assets = select_assets_for_single_market( + debug!("Running {commodity_id}|{region_id} selection pass 1"); + let selected_assets = select_assets_for_single_market( model, - commodity_id, - region_id, + &commodity_id, + ®ion_id, year, - ¤t_demand, + &net_demand, existing_assets, prices, + &[], writer, )?; - assets_for_cycle.insert((commodity_id.clone(), region_id.clone()), assets); + debug!("Completed {commodity_id}|{region_id} selection pass 1"); + + // If no assets have been selected, skip to the next market + if selected_assets.is_empty() { + debug!("No assets selected for '{commodity_id}|{region_id}'"); + continue; + } + assets_for_first_pass.extend(selected_assets.iter().cloned()); + retained_first_pass_assets.extend( + selected_assets + .iter() + .filter(|asset| { + feedback_processes + .contains(&(asset.process_id().clone(), asset.region_id().clone())) + }) + .cloned(), + ); + + // Run dispatch + debug!("Running cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 1"); + let solution = DispatchRun::new(model, &selected_assets, year, &net_demand) + .without_commodity_constraints() + .with_market_balance_subset(std::slice::from_ref(market)) + .run( + &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 1"), + writer, + ) + .with_context(|| format!("Dispatch failed for cycle ({markets_str})"))?; + debug!("Completed cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 1"); + + // Update demand map with flows from newly selected assets + let flows = solution.create_flow_map(); + for ((asset, commodity_id, time_slice), flow) in &flows { + if feedback_processes.contains(&(asset.process_id().clone(), asset.region_id().clone())) + { + excluded_flows + .entry((asset.clone(), commodity_id.clone(), time_slice.clone())) + .and_modify(|stored_flow| *stored_flow += *flow) + .or_insert(*flow); + } + } + update_net_demand_map(&mut net_demand, &flows, &selected_assets); + } - // Assemble full list of assets for dispatch (previously selected + all chosen so far) - let mut all_assets = previously_selected_assets.to_vec(); - let assets_for_cycle_flat: Vec<_> = assets_for_cycle - .values() - .flat_map(|v| v.iter().cloned()) + // STEP 2 + // Iterate over the markets in order2, excluding feedback processes + let mut net_demand = demand.clone(); + update_net_demand_map(&mut net_demand, &excluded_flows, &assets_for_first_pass); + let mut assets_for_second_pass = Vec::new(); + for market in markets { + let (commodity_id, region_id) = market.clone(); + let feedback_processes: Vec = feedback_processes + .iter() + .filter_map(|(p, r)| (*r == region_id).then_some(p.clone())) .collect(); - all_assets.extend_from_slice(&assets_for_cycle_flat); - // We balance all previously seen markets plus all cycle markets up to and including this one - let mut markets_to_balance = seen_markets.to_vec(); - markets_to_balance.extend_from_slice(&markets[0..=idx]); + // Select assets for this market + debug!("Running {commodity_id}|{region_id} selection pass 2"); + let selected_assets = select_assets_for_single_market( + model, + &commodity_id, + ®ion_id, + year, + &net_demand, + existing_assets, + prices, + &feedback_processes, + writer, + )?; + debug!("Completed {commodity_id}|{region_id} selection pass 2"); + + // If no assets have been selected, skip to the next market + if selected_assets.is_empty() { + debug!("No assets selected for '{commodity_id}|{region_id}'"); + continue; + } + assets_for_second_pass.extend(selected_assets.iter().cloned()); + + // // Assemble full list of assets for dispatch (previously selected + all chosen so far) + // let mut all_assets = previously_selected_assets.to_vec(); + // all_assets.extend(assets_for_first_pass.iter().cloned()); + // all_assets.extend(assets_for_second_pass.iter().cloned()); + + // // We balance all previously seen markets plus all cycle markets up to and including this one + // let mut markets_to_balance = seen_markets.to_vec(); + // markets_to_balance.extend_from_slice(&investment_order[0..=idx]); // Run dispatch - let solution = DispatchRun::new(model, &all_assets, year) + debug!("Running cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 2"); + let solution = DispatchRun::new(model, &selected_assets, year, &net_demand) .without_commodity_constraints() - .with_market_balance_subset(&markets_to_balance) + .with_market_balance_subset(std::slice::from_ref(market)) .run( - &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment"), + &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 2"), writer, ) .with_context(|| format!("Dispatch failed for cycle ({markets_str})"))?; + debug!("Completed cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 2"); - // Calculate new net demand map with all assets selected so far - current_demand.clone_from(demand); + // Update demand map with flows from newly selected assets update_net_demand_map( - &mut current_demand, + &mut net_demand, &solution.create_flow_map(), - &assets_for_cycle_flat, + &selected_assets, ); } - // Collect assets - let all_cycle_assets: Vec<_> = assets_for_cycle.into_values().flatten().collect(); - Ok(all_cycle_assets) + Ok(retained_first_pass_assets + .into_iter() + .chain(assets_for_second_pass) + .collect()) } /// Get a portion of the demand profile for this market diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 7c464b085..5e63af895 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -8,6 +8,7 @@ use crate::model::Model; use crate::output::DataWriter; use crate::region::RegionID; use crate::simulation::PriceMap; +use crate::simulation::investment::AllDemandMap; use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; use crate::units::{Activity, Flow, Money, MoneyPerActivity, MoneyPerFlow}; use anyhow::{Context, Result, anyhow, bail}; @@ -401,6 +402,7 @@ pub struct DispatchRun<'model, 'run> { existing_assets: &'run [AssetRef], candidate_assets: &'run [AssetRef], markets_to_balance: &'run [(CommodityID, RegionID)], + market_demands: &'run AllDemandMap, input_prices: Option<&'run PriceMap>, include_commodity_constraints: bool, year: u32, @@ -408,12 +410,18 @@ pub struct DispatchRun<'model, 'run> { impl<'model, 'run> DispatchRun<'model, 'run> { /// Create a new [`DispatchRun`] for the specified model and assets for a given year - pub fn new(model: &'model Model, assets: &'run [AssetRef], year: u32) -> Self { + pub fn new( + model: &'model Model, + assets: &'run [AssetRef], + year: u32, + market_demands: &'run AllDemandMap, + ) -> Self { Self { model, existing_assets: assets, candidate_assets: &[], markets_to_balance: &[], + market_demands, input_prices: None, include_commodity_constraints: true, year, @@ -691,6 +699,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { self.model, &all_assets, markets_to_balance, + self.market_demands, self.year, self.candidate_assets, include_commodity_constraints, diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs index 17d89ed96..07620ec36 100644 --- a/src/simulation/optimisation/constraints.rs +++ b/src/simulation/optimisation/constraints.rs @@ -5,6 +5,7 @@ use crate::commodity::{BalanceType, CommodityID, CommodityType}; use crate::model::Model; use crate::process::FlowDirection; use crate::region::RegionID; +use crate::simulation::investment::AllDemandMap; use crate::time_slice::{Season, TimeSliceInfo, TimeSliceSelection}; use crate::units::{Flow, MoneyPerCapacityPerYear, UnitType, Year}; use highs::RowProblem as Problem; @@ -88,6 +89,7 @@ pub fn add_model_constraints<'a, I>( model: &'a Model, assets: &I, markets_to_balance: &'a [(CommodityID, RegionID)], + market_demands: &AllDemandMap, year: u32, candidate_assets: &'a [AssetRef], include_commodity_constraints: bool, @@ -101,7 +103,7 @@ where model, assets, markets_to_balance, - year, + market_demands, candidate_assets, ); @@ -330,7 +332,7 @@ fn add_commodity_balance_constraints<'a, I>( model: &'a Model, assets: &I, markets_to_balance: &'a [(CommodityID, RegionID)], - year: u32, + market_demands: &AllDemandMap, candidate_assets: &'a [AssetRef], ) -> CommodityBalanceKeys where @@ -394,14 +396,18 @@ where model.parameters.commodity_balance_epsilon, ); - // For SVD commodities, the lower bound is the exogenous demand (or epsilon if larger). - // For SED commodities, the lower bound is just epsilon. - let min = match commodity.kind { - CommodityType::ServiceDemand => { - commodity.demand[&(region_id.clone(), year, ts_selection.clone())].max(epsilon) - } - _ => epsilon, + // For SVD commodities, the demand must be present in the map; for SED commodities, + // missing entries mean no demand for this selection. + let key = ( + commodity_id.clone(), + region_id.clone(), + ts_selection.clone(), + ); + let demand_for_selection = match commodity.kind { + CommodityType::ServiceDemand => market_demands[&key], + _ => market_demands.get(&key).copied().unwrap_or(Flow(0.0)), }; + let min = demand_for_selection.max(epsilon); // Consume collected terms into a row. `terms.drain(..)` ensures the vector is // emptied for the next selection. From 4e4ea0b568c3e0c3a3942e8a59541ebea6a51855 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 14:35:43 +0100 Subject: [PATCH 05/13] Fix MILP error --- src/graph/investment.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index 44ba78b19..f987a7abd 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -354,26 +354,26 @@ fn order_sccs( } // Every SCC node must retain at least one correctly-ordered incoming edge. - for j in 0..n { + for i in 0..n { let mut incoming_terms = Vec::new(); - for i in 0..n { + for j in 0..n { if i == j { continue; } - // We need to know whether the original graph contains i -> j. - // If so, x[i][j] represents that edge being retained. + // We need to know whether the original graph contains j -> i. if original_graph - .find_edge(original_indices[i], original_indices[j]) + .find_edge(original_indices[j], original_indices[i]) .is_some() { + // Variable saying whether i comes before j incoming_terms.push((vars[i][j].unwrap(), 1.0)); } } // If the node has an incoming edge from outside the SCC, then it doesn't need a // correctly-ordered internal incoming edge. Otherwise it does. - let required = if has_external_incoming[j] { 0.0 } else { 1.0 }; + let required = if has_external_incoming[i] { 0.0 } else { 1.0 }; problem.add_row(required.., incoming_terms); } From 19ef0a08bfbe1b5311c599fb077597af8f8ee88b Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 14:41:01 +0100 Subject: [PATCH 06/13] Fix for select_assets_for_cycle --- src/simulation/market.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/simulation/market.rs b/src/simulation/market.rs index 7f5c019b1..a911bd300 100644 --- a/src/simulation/market.rs +++ b/src/simulation/market.rs @@ -264,6 +264,9 @@ pub fn select_assets_for_cycle( let market_order: HashMap<_, _> = markets.iter().enumerate().map(|(i, m)| (m, i)).collect(); let mut feedback_processes: Vec<(ProcessID, RegionID)> = Vec::new(); for (process_id, process) in &model.processes { + if !process.active_for_year(year) { + continue; + } let Some(primary_output) = &process.primary_output else { continue; }; From a2158c4d46cf828b41693726e7b12f7406dbe88b Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 14:48:50 +0100 Subject: [PATCH 07/13] Comments to order_scc --- src/graph/investment.rs | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index f987a7abd..f1b9bd30a 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -162,6 +162,8 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { /// `i` comes before `j`, then `j` cannot be before `i`). /// * Transitivity constraints prevent 3-cycles, ensuring the resulting relation is acyclic (i.e. if /// `i` comes before `j` and `j` comes before `k`, then `k` cannot come before `i`). +/// * Incoming-edge constraints require every market without an external incoming edge to retain at +/// least one correctly ordered incoming edge within the SCC. /// * The objective minimises the number of “forward” edges (edges that would point from an earlier /// market to a later one), counted within the original SCC and treated as unit penalties. A small /// bias (<1) is added to nudge exporters earlier without outweighing the main objective (a bias @@ -361,18 +363,20 @@ fn order_sccs( continue; } - // We need to know whether the original graph contains j -> i. + // Check if the original graph contains j -> i. if original_graph .find_edge(original_indices[j], original_indices[i]) .is_some() { - // Variable saying whether i comes before j + // Get the variable saying whether i comes before j (i.e. "correct" ordering + // for a j -> i edge) and add it to the terms incoming_terms.push((vars[i][j].unwrap(), 1.0)); } } // If the node has an incoming edge from outside the SCC, then it doesn't need a - // correctly-ordered internal incoming edge. Otherwise it does. + // correctly-ordered internal incoming edge. Otherwise it does (i.e. a least one + // variable in `incoming_terms` must be 1.0) let required = if has_external_incoming[i] { 0.0 } else { 1.0 }; problem.add_row(required.., incoming_terms); } From a769b1162999d941ebaa67a68fa0d1870843085e Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 15:08:33 +0100 Subject: [PATCH 08/13] order_sccs returns Result --- src/graph/investment.rs | 63 ++++++++++++++++++++++++++--------------- src/input.rs | 3 +- 2 files changed, 42 insertions(+), 24 deletions(-) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index f1b9bd30a..6efdc0def 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -3,9 +3,9 @@ use super::{CommoditiesGraph, GraphEdge, GraphNode}; use crate::commodity::{CommodityMap, CommodityType}; use crate::region::RegionID; use crate::simulation::market::MarketSet; +use anyhow::{Result, bail}; use highs::{Col, HighsModelStatus, RowProblem, Sense}; use indexmap::IndexMap; -use log::warn; use petgraph::algo::{condensation, toposort}; use petgraph::graph::Graph; use petgraph::prelude::NodeIndex; @@ -41,12 +41,12 @@ fn solve_investment_order_for_year( graphs: &IndexMap<(RegionID, u32), CommoditiesGraph>, commodities: &CommodityMap, year: u32, -) -> Vec { +) -> Result> { // Initialise InvestmentGraph for this year from the set of original `CommodityGraph`s let mut investment_graph = init_investment_graph_for_year(graphs, year, commodities); // Condense strongly connected components - investment_graph = compress_cycles(&investment_graph); + investment_graph = compress_cycles(&investment_graph)?; // Perform a topological sort on the condensed graph // We can safely unwrap because `toposort` will only return an error in case of cycles, which @@ -54,7 +54,7 @@ fn solve_investment_order_for_year( let order = toposort(&investment_graph, None).unwrap(); // Compute layers for investment - compute_layers(&investment_graph, &order) + Ok(compute_layers(&investment_graph, &order)) } /// Initialise an `InvestmentGraph` for the given year from a set of `CommodityGraph`s @@ -115,15 +115,15 @@ fn init_investment_graph_for_year( } /// Compresses cycles into `MarketSet::Cycle` nodes -fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { +fn compress_cycles(graph: &InvestmentGraph) -> Result { // Detect strongly connected components let mut condensed_graph = condensation(graph.clone(), true); // Order nodes within each strongly connected component - order_sccs(&mut condensed_graph, graph); + order_sccs(&mut condensed_graph, graph)?; // Map to a new InvestmentGraph - condensed_graph.map( + Ok(condensed_graph.map( // Map nodes to MarketSet // If only one member, keep as-is; if multiple members, create Cycle |_, node_weight| match node_weight.len() { @@ -139,7 +139,7 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { }, // Keep edges the same |_, edge_weight| edge_weight.clone(), - ) + )) } /// Order the members of each strongly connected component using a mixed-integer linear program. @@ -162,8 +162,8 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { /// `i` comes before `j`, then `j` cannot be before `i`). /// * Transitivity constraints prevent 3-cycles, ensuring the resulting relation is acyclic (i.e. if /// `i` comes before `j` and `j` comes before `k`, then `k` cannot come before `i`). -/// * Incoming-edge constraints require every market without an external incoming edge to retain at -/// least one correctly ordered incoming edge within the SCC. +/// * Incoming-edge constraints require every SCC node to retain at least one correctly-ordered +/// incoming edge within the SCC unless it has an external incoming edge. /// * The objective minimises the number of “forward” edges (edges that would point from an earlier /// market to a later one), counted within the original SCC and treated as unit penalties. A small /// bias (<1) is added to nudge exporters earlier without outweighing the main objective (a bias @@ -237,7 +237,7 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { fn order_sccs( condensed_graph: &mut Graph, GraphEdge>, original_graph: &InvestmentGraph, -) { +) -> Result<()> { const EXTERNAL_BIAS: f64 = 0.1; // Map each market set back to the node index in the original graph so we can inspect edges. @@ -297,6 +297,10 @@ fn order_sccs( } } + // Make sure the SCC has at least one external incoming edge - panic if not as this should + // not happen + assert!(has_external_incoming.contains(&true), "SCC has no inputs"); + // Bias: if market j has outgoing edges to nodes outside this SCC, we prefer to place it earlier. for (j, has_external) in has_external_outgoing.iter().enumerate() { if *has_external { @@ -381,19 +385,30 @@ fn order_sccs( problem.add_row(required.., incoming_terms); } + // String representation of SCC for error messages + let scc_display = format!( + "({})", + scc.iter() + .map(ToString::to_string) + .collect::>() + .join(", ") + ); + + // Solve model let model = problem.optimise(Sense::Minimise); let solved = match model.try_solve() { Ok(solved) => solved, Err(status) => { - warn!("HiGHS failed while ordering an SCC: {status:?}"); - continue; + bail!("HiGHS failed while ordering SCC {scc_display}: {status:?}"); } }; + // Check status if solved.status() != HighsModelStatus::Optimal { let status = solved.status(); - warn!("HiGHS returned a non-optimal status while ordering an SCC: {status:?}"); - continue; + bail!( + "HiGHS returned a non-optimal status while ordering SCC {scc_display}: {status:?}" + ); } let solution = solved.get_solution(); @@ -421,6 +436,8 @@ fn order_sccs( .map(|idx| original_order[idx].clone()) .collect(); } + + Ok(()) } /// Compute layers of market sets from the topological order @@ -524,13 +541,13 @@ pub fn solve_investment_order_for_model( commodity_graphs: &IndexMap<(RegionID, u32), CommoditiesGraph>, commodities: &CommodityMap, years: &[u32], -) -> HashMap> { +) -> Result>> { let mut investment_orders = HashMap::new(); for year in years { - let order = solve_investment_order_for_year(commodity_graphs, commodities, *year); + let order = solve_investment_order_for_year(commodity_graphs, commodities, *year)?; investment_orders.insert(*year, order); } - investment_orders + Ok(investment_orders) } #[cfg(test)] @@ -571,7 +588,7 @@ mod tests { let mut condensed: Graph, GraphEdge> = Graph::new(); let component = condensed.add_node(markets.to_vec()); - order_sccs(&mut condensed, &original); + order_sccs(&mut condensed, &original).unwrap(); // Expected order corresponds to the example in the doc comment. // Note that C should be first, as it has an outgoing edge to the external market. @@ -602,7 +619,7 @@ mod tests { commodities.insert("C".into(), Arc::new(svd_commodity)); let graphs = IndexMap::from([(("GBR".into(), 2020), graph)]); - let result = solve_investment_order_for_year(&graphs, &commodities, 2020); + let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap(); // Expected order: C, B, A (leaf nodes first) // No cycles or layers, so all market sets should be `Single` @@ -630,7 +647,7 @@ mod tests { commodities.insert("B".into(), Arc::new(sed_commodity)); let graphs = IndexMap::from([(("GBR".into(), 2020), graph)]); - let result = solve_investment_order_for_year(&graphs, &commodities, 2020); + let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap(); // Should be a single `Cycle` market set containing both commodities assert_eq!(result.len(), 1); @@ -669,7 +686,7 @@ mod tests { commodities.insert("D".into(), Arc::new(svd_commodity)); let graphs = IndexMap::from([(("GBR".into(), 2020), graph)]); - let result = solve_investment_order_for_year(&graphs, &commodities, 2020); + let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap(); // Expected order: D, Layer(B, C), A assert_eq!(result.len(), 3); @@ -708,7 +725,7 @@ mod tests { (("GBR".into(), 2020), graph.clone()), (("FRA".into(), 2020), graph), ]); - let result = solve_investment_order_for_year(&graphs, &commodities, 2020); + let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap(); // Expected order: Should have three layers, each with two commodities (one per region) assert_eq!(result.len(), 3); diff --git a/src/input.rs b/src/input.rs index 7e90e5b08..b322d2424 100644 --- a/src/input.rs +++ b/src/input.rs @@ -296,7 +296,8 @@ pub fn load_model>(model_dir: P) -> Result { )?; // Solve investment order for each region/year - let investment_order = solve_investment_order_for_model(&commodity_graphs, &commodities, years); + let investment_order = + solve_investment_order_for_model(&commodity_graphs, &commodities, years)?; let model_path = model_dir .as_ref() From 541dc2aa765523601bfb2961e4c86f1dd809b79c Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 15:32:28 +0100 Subject: [PATCH 09/13] Fix tests --- src/graph/investment.rs | 56 +++++++++++++++++++++++++++-------------- 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index 6efdc0def..ee840590b 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -183,7 +183,7 @@ fn compress_cycles(graph: &InvestmentGraph) -> Result { /// A ← B ← C ← A /// ``` /// -/// Additionally, C has an outgoing edge to a node outside the cycle. +/// Additionally, C has an outgoing edge to a node outside the cycle, and B has an incoming edge. /// /// The costs matrix in the MILP is set up to penalise any edge that points “forward” in the final /// order: if there's an edge from X to Y we prefer to place Y before X so the edge points backwards: @@ -254,6 +254,15 @@ fn order_sccs( continue; } + // String representation of SCC for error messages + let scc_display = format!( + "({})", + scc.iter() + .map(ToString::to_string) + .collect::>() + .join(", ") + ); + // Capture current order and resolve each market set back to its original graph index. let original_order = scc.clone(); let original_indices = original_order @@ -299,7 +308,10 @@ fn order_sccs( // Make sure the SCC has at least one external incoming edge - panic if not as this should // not happen - assert!(has_external_incoming.contains(&true), "SCC has no inputs"); + assert!( + has_external_incoming.contains(&true), + "SCC {scc_display} has no inputs" + ); // Bias: if market j has outgoing edges to nodes outside this SCC, we prefer to place it earlier. for (j, has_external) in has_external_outgoing.iter().enumerate() { @@ -385,15 +397,6 @@ fn order_sccs( problem.add_row(required.., incoming_terms); } - // String representation of SCC for error messages - let scc_display = format!( - "({})", - scc.iter() - .map(ToString::to_string) - .collect::>() - .join(", ") - ); - // Solve model let model = problem.optimise(Sense::Minimise); let solved = match model.try_solve() { @@ -576,14 +579,24 @@ mod tests { GraphEdge::Primary("process1".into()), ); } - // External market receiving exports from C; encourages C to appear early. - let external = original.add_node(MarketSet::Single(("X".into(), "GBR".into()))); + // Downstream market receives exports from C; encourages C to appear early. + let downstream = original.add_node(MarketSet::Single(("X".into(), "GBR".into()))); original.add_edge( node_indices[2], - external, + downstream, GraphEdge::Primary("process2".into()), ); + // B receives imports from an upstream market, which permits internal B-producers to come + // before B in the investment order. Note: without at least one upstream input SCC ordering + // will fail + let upstream = original.add_node(MarketSet::Single(("Y".into(), "GBR".into()))); + original.add_edge( + upstream, + node_indices[1], + GraphEdge::Primary("process3".into()), + ); + // Single SCC containing all markets. let mut condensed: Graph, GraphEdge> = Graph::new(); let component = condensed.add_node(markets.to_vec()); @@ -592,6 +605,7 @@ mod tests { // Expected order corresponds to the example in the doc comment. // Note that C should be first, as it has an outgoing edge to the external market. + // B is last because it receives inputs from upstream let expected = ["C", "A", "B"] .map(|id| MarketSet::Single((id.into(), "GBR".into()))) .to_vec(); @@ -631,30 +645,34 @@ mod tests { #[rstest] fn solve_investment_order_cyclic_graph(sed_commodity: Commodity) { - // Create a simple cyclic graph: A -> B -> A + // Create a simple cyclic graph: X -> A -> B -> A let mut graph = Graph::new(); + let node_x = graph.add_node(GraphNode::Commodity("X".into())); let node_a = graph.add_node(GraphNode::Commodity("A".into())); let node_b = graph.add_node(GraphNode::Commodity("B".into())); - // Add edges creating a cycle: A -> B -> A + // Add edges creating a cycle: X -> A -> B -> A + graph.add_edge(node_x, node_a, GraphEdge::Primary("process0".into())); graph.add_edge(node_a, node_b, GraphEdge::Primary("process1".into())); graph.add_edge(node_b, node_a, GraphEdge::Primary("process2".into())); // Create commodities map using fixtures let mut commodities = CommodityMap::new(); + commodities.insert("X".into(), Arc::new(sed_commodity.clone())); commodities.insert("A".into(), Arc::new(sed_commodity.clone())); commodities.insert("B".into(), Arc::new(sed_commodity)); let graphs = IndexMap::from([(("GBR".into(), 2020), graph)]); let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap(); - // Should be a single `Cycle` market set containing both commodities - assert_eq!(result.len(), 1); + // Should be a `Cycle` market set containing both commodities, followed by a `Single` + assert_eq!(result.len(), 2); assert_eq!( result[0], - MarketSet::Cycle(vec![("A".into(), "GBR".into()), ("B".into(), "GBR".into())]) + MarketSet::Cycle(vec![("B".into(), "GBR".into()), ("A".into(), "GBR".into())]) ); + assert_eq!(result[1], MarketSet::Single(("X".into(), "GBR".into()))); } #[rstest] From 0038735eec57a59a59fff8cf4def4f4fe73cbebd Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 22 Sep 2026 16:12:19 +0100 Subject: [PATCH 10/13] Improve docstring --- src/graph/investment.rs | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index ee840590b..d65fe35da 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -183,7 +183,8 @@ fn compress_cycles(graph: &InvestmentGraph) -> Result { /// A ← B ← C ← A /// ``` /// -/// Additionally, C has an outgoing edge to a node outside the cycle, and B has an incoming edge. +/// Additionally, C has an outgoing edge to a node downstream of the cycle, and B has an incoming +/// edge from upstream. /// /// The costs matrix in the MILP is set up to penalise any edge that points “forward” in the final /// order: if there's an edge from X to Y we prefer to place Y before X so the edge points backwards: @@ -206,8 +207,19 @@ fn compress_cycles(graph: &InvestmentGraph) -> Result { /// | C | 0 | 1 | 0 | /// ``` /// -/// Solving this problem with binary decision variables for each `x[i][j]`, and constraints to enforce -/// antisymmetry and transitivity, yields optimal decision variables of: +/// Additionally, each node must retain at least one incoming edge from a node that appears after it +/// in the final investment order. This includes edges within the SCC and any external incoming +/// edges. In this example, since B has an incoming edge from outside the SCC, this constraint +/// implies that only A and C must retain a correctly-ordered incoming edge from within the SCC: +/// +/// ```text +/// x[A][B] >= 1 +/// x[B][C] >= 0 +/// x[C][A] >= 1 +/// ``` +/// +/// Solving this problem with binary decision variables for each `x[i][j]`, and additional +/// constraints to enforce antisymmetry and transitivity, yields optimal decision variables of: /// /// ```text /// x[A][B] = 1 (A before B) @@ -230,6 +242,7 @@ fn compress_cycles(graph: &InvestmentGraph) -> Result { /// * The preference towards having exporter markets early in the order keeps C at the front. /// * As with any SCC, at least one pairwise violation is guaranteed. In this ordering, the only /// pairwise violation is between B and C, as C is solved before B, but B may consume C. +/// * This pairwise violation is permitted because B has an additional producer upstream of the SCC. /// /// The resulting order replaces the original `MarketSet::Cycle` entry inside the condensed /// graph, providing a deterministic processing sequence for downstream logic. From 7857219822861ad2ca34fef07a2f9f0b686f17b6 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Wed, 23 Sep 2026 15:21:20 +0100 Subject: [PATCH 11/13] Use shadow prices in intermediate dispatch --- src/simulation/market.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/simulation/market.rs b/src/simulation/market.rs index a911bd300..1ab49861b 100644 --- a/src/simulation/market.rs +++ b/src/simulation/market.rs @@ -341,6 +341,7 @@ pub fn select_assets_for_cycle( let solution = DispatchRun::new(model, &selected_assets, year, &net_demand) .without_commodity_constraints() .with_market_balance_subset(std::slice::from_ref(market)) + .with_input_prices(&prices.shadow) .run( &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 1"), writer, @@ -410,6 +411,7 @@ pub fn select_assets_for_cycle( let solution = DispatchRun::new(model, &selected_assets, year, &net_demand) .without_commodity_constraints() .with_market_balance_subset(std::slice::from_ref(market)) + .with_input_prices(&prices.shadow) .run( &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 2"), writer, From 18bd64004d097cf82ebe800e86e23d3923d3a1e8 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Wed, 23 Sep 2026 15:58:53 +0100 Subject: [PATCH 12/13] Add TODO --- src/graph/investment.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index d65fe35da..e62b2a824 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -312,6 +312,7 @@ fn order_sccs( } // Check whether this node has any incoming edges from outside the SCC + // TODO: this does not currently consider SOURCE nodes - it should! for edge in original_graph.edges_directed(idx, Direction::Incoming) { if !index_position.contains_key(&edge.source()) { has_external_incoming[i] = true; From b5b0b292ebd380c24476911b8e5bccb8e580ed43 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Thu, 24 Sep 2026 09:57:33 +0100 Subject: [PATCH 13/13] TEMPORARY HACK: allow order_sccs to consider source nodes --- src/graph/investment.rs | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/src/graph/investment.rs b/src/graph/investment.rs index e62b2a824..ec4937e9b 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -46,7 +46,7 @@ fn solve_investment_order_for_year( let mut investment_graph = init_investment_graph_for_year(graphs, year, commodities); // Condense strongly connected components - investment_graph = compress_cycles(&investment_graph)?; + investment_graph = compress_cycles(&investment_graph, graphs, year)?; // Perform a topological sort on the condensed graph // We can safely unwrap because `toposort` will only return an error in case of cycles, which @@ -115,12 +115,16 @@ fn init_investment_graph_for_year( } /// Compresses cycles into `MarketSet::Cycle` nodes -fn compress_cycles(graph: &InvestmentGraph) -> Result { +fn compress_cycles( + graph: &InvestmentGraph, + graphs: &IndexMap<(RegionID, u32), CommoditiesGraph>, + year: u32, +) -> Result { // Detect strongly connected components let mut condensed_graph = condensation(graph.clone(), true); // Order nodes within each strongly connected component - order_sccs(&mut condensed_graph, graph)?; + order_sccs(&mut condensed_graph, graph, graphs, year)?; // Map to a new InvestmentGraph Ok(condensed_graph.map( @@ -250,6 +254,8 @@ fn compress_cycles(graph: &InvestmentGraph) -> Result { fn order_sccs( condensed_graph: &mut Graph, GraphEdge>, original_graph: &InvestmentGraph, + commodity_graphs: &IndexMap<(RegionID, u32), CommoditiesGraph>, + year: u32, ) -> Result<()> { const EXTERNAL_BIAS: f64 = 0.1; @@ -299,6 +305,20 @@ fn order_sccs( let mut has_external_outgoing: Vec = vec![false; n]; let mut has_external_incoming: Vec = vec![false; n]; for (i, &idx) in original_indices.iter().enumerate() { + let (commodity_id, region_id) = original_graph + .node_weight(idx) + .unwrap() + .iter_markets() + .next() + .unwrap(); + let commodity_graph = commodity_graphs.get(&(region_id.clone(), year)); + let commodity_node = commodity_graph.and_then(|commodity_graph| { + commodity_graph.node_indices().find(|&node| { + commodity_graph.node_weight(node) + == Some(&GraphNode::Commodity(commodity_id.clone())) + }) + }); + // Loop over the edges going out of this node for edge in original_graph.edges_directed(idx, Direction::Outgoing) { // If the target j is inside this SCC, record a penalty for putting i before j @@ -311,8 +331,15 @@ fn order_sccs( } } - // Check whether this node has any incoming edges from outside the SCC - // TODO: this does not currently consider SOURCE nodes - it should! + // Check whether this node has any incoming edges from outside the SCC or from SOURCE. + has_external_incoming[i] = + commodity_graph + .zip(commodity_node) + .is_some_and(|(commodity_graph, node)| { + commodity_graph + .edges_directed(node, Direction::Incoming) + .any(|edge| commodity_graph[edge.source()] == GraphNode::Source) + }); for edge in original_graph.edges_directed(idx, Direction::Incoming) { if !index_position.contains_key(&edge.source()) { has_external_incoming[i] = true; @@ -327,7 +354,7 @@ fn order_sccs( "SCC {scc_display} has no inputs" ); - // Bias: if market j has outgoing edges to nodes outside this SCC, we prefer to place it earlier. + // Bias: if market j has outgoing edges to nodes outside the SCC, we prefer to place it earlier. for (j, has_external) in has_external_outgoing.iter().enumerate() { if *has_external { for (row_idx, row) in penalties.iter_mut().enumerate() { @@ -615,7 +642,7 @@ mod tests { let mut condensed: Graph, GraphEdge> = Graph::new(); let component = condensed.add_node(markets.to_vec()); - order_sccs(&mut condensed, &original).unwrap(); + order_sccs(&mut condensed, &original, &IndexMap::new(), 2020).unwrap(); // Expected order corresponds to the example in the doc comment. // Note that C should be first, as it has an outgoing edge to the external market.