From f26d12148118ac0fd4dfc7819d0c9b0393748e91 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Thu, 30 Jul 2026 16:03:13 +0100 Subject: [PATCH 1/6] Add L1 regularisation to dispatch --- schemas/input/model.yaml | 8 ++++ src/model/parameters.rs | 13 +++++++ src/simulation/optimisation.rs | 68 ++++++++++++++++++++++++++++++++-- 3 files changed, 86 insertions(+), 3 deletions(-) diff --git a/schemas/input/model.yaml b/schemas/input/model.yaml index 5120a4a9a..31d80d437 100644 --- a/schemas/input/model.yaml +++ b/schemas/input/model.yaml @@ -21,6 +21,14 @@ properties: default: 1e-6 notes: | The default value should work. Do not change this value unless you know what you're doing! + dispatch_activity_equalisation_epsilon: + type: number + description: Small penalty on L1 utilisation deviation to encourage even dispatch across assets + default: 1e-6 + notes: | + Discourages corner solutions where all activity is concentrated on a single asset when + equally good alternatives exist. Must be much smaller than the smallest meaningful cost + difference between assets. Set to zero to disable. capacity_limit_factor: type: number description: Scales the capacity assigned to candidate assets diff --git a/src/model/parameters.rs b/src/model/parameters.rs index ef542d3cb..b328d3515 100644 --- a/src/model/parameters.rs +++ b/src/model/parameters.rs @@ -82,6 +82,11 @@ pub struct ModelParameters { /// /// Don't change unless you know what you're doing. pub commodity_balance_epsilon: Flow, + /// Small penalty on L1 utilisation deviation to encourage even dispatch across assets. + /// + /// Set to zero to disable. Should be much smaller than the smallest meaningful cost + /// difference between assets. + pub dispatch_activity_equalisation_epsilon: f64, /// Affects the maximum capacity that can be given to a newly created asset. /// /// It is the proportion of maximum capacity that could be required across time slices. @@ -130,6 +135,7 @@ impl Default for ModelParameters { allow_dangerous_options: false, candidate_asset_capacity: Capacity(1e-4), commodity_balance_epsilon: Flow(1e-6), + dispatch_activity_equalisation_epsilon: 1e6, capacity_limit_factor: Dimensionless(0.05), fallback_pricing_strategy: PricingStrategy::FullCostAverage, value_of_lost_load: MoneyPerFlow(1e9), @@ -356,6 +362,13 @@ impl ModelParameters { "commodity_balance_epsilon must be a finite number greater than or equal to zero" ); + // dispatch_activity_equalisation_epsilon + ensure!( + self.dispatch_activity_equalisation_epsilon.is_finite() + && self.dispatch_activity_equalisation_epsilon >= 0.0, + "dispatch_activity_equalisation_epsilon must be a finite number greater than or equal to zero" + ); + // value_of_lost_load check_value_of_lost_load(self.value_of_lost_load)?; diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index e58b59abf..568767587 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -691,20 +691,33 @@ impl<'model, 'run> DispatchRun<'model, 'run> { } // Add constraints - let all_assets = chain( + let all_assets: Vec<_> = chain( existing_assets_with_parents.iter(), self.candidate_assets.iter(), - ); + ) + .collect(); let constraint_keys = add_model_constraints( &mut problem, &variables, self.model, - &all_assets, + &all_assets.iter().copied(), markets_to_balance, self.year, self.candidate_assets, ); + // Add secondary objective to encourage even dispatch across assets, if enabled + let epsilon = self.model.parameters.dispatch_activity_equalisation_epsilon; + if epsilon > 0.0 { + add_activity_equalisation( + &mut problem, + &variables, + all_assets.iter().copied(), + self.model, + epsilon, + ); + } + // Create model and apply any user-supplied HiGHS options to it let mut model = problem.optimise(Sense::Minimise); apply_highs_options_from_toml(&mut model, &self.model.parameters.highs.dispatch_options) @@ -761,6 +774,55 @@ fn add_activity_variables( start..problem.num_cols() } +/// Add secondary objective variables and constraints to encourage even dispatch across assets. +/// +/// The secondary objective is `epsilon * mean(d)`, where each `d[a,t]` is the L1 deviation of +/// asset `a`'s utilisation fraction from a single free centre variable `m`: +/// +/// ```text +/// d[a,t] = |act[a,t] / cap[a] - m| +/// ``` +/// +/// `|x|` is linearised via two constraints per `d`: `d >= x` and `d >= -x` (both with `d >= 0`). +/// Both are needed because we minimise `d` — without the second constraint HiGHS would set `d = 0` +/// whenever `x < 0`, missing half the deviation. +/// +/// `m` is a free variable (bounded below at zero since all utilisation fractions are +/// non-negative). The LP sets it to the L1-optimal centre (the median of utilisations), so the +/// penalty measures spread rather than deviation from any fixed target. This makes `epsilon` +/// independent of the absolute utilisation level, which is determined by the primary objective. +fn add_activity_equalisation<'a, I>( + problem: &mut Problem, + variables: &VariableMap, + assets: I, + model: &Model, + epsilon: f64, +) where + I: Iterator, +{ + // m is the free L1-centre variable in utilisation-fraction space; bounded below at zero + // since all utilisation fractions are non-negative + let m = problem.add_column(0.0, 0.0..); + + for asset in assets { + // Flexible capacity assets use their current capacity as an approximation + let cap = asset.total_capacity().value(); + if cap <= 1e-9 { + continue; + } + let inv_cap = 1.0 / cap; + + for time_slice in model.time_slice_info.iter_ids() { + let act = variables.get_activity_var(asset, time_slice); + let d = problem.add_column(epsilon, 0.0..); + // d >= act/cap - m + problem.add_row(0.0.., [(d, 1.0), (act, -inv_cap), (m, 1.0)]); + // d >= m - act/cap + problem.add_row(0.0.., [(d, 1.0), (act, inv_cap), (m, -1.0)]); + } + } +} + fn add_capacity_variables( problem: &mut Problem, variables: &mut CapacityVariableMap, From f481607d02adf919990f179379f1df8af41a5086 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Thu, 30 Jul 2026 16:13:03 +0100 Subject: [PATCH 2/6] Fix typo, correct docstring, and exclude candidates from regularisation --- src/model/parameters.rs | 2 +- src/simulation/optimisation.rs | 16 +++++++++------- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/model/parameters.rs b/src/model/parameters.rs index b328d3515..610a7a169 100644 --- a/src/model/parameters.rs +++ b/src/model/parameters.rs @@ -135,7 +135,7 @@ impl Default for ModelParameters { allow_dangerous_options: false, candidate_asset_capacity: Capacity(1e-4), commodity_balance_epsilon: Flow(1e-6), - dispatch_activity_equalisation_epsilon: 1e6, + dispatch_activity_equalisation_epsilon: 1e-6, capacity_limit_factor: Dimensionless(0.05), fallback_pricing_strategy: PricingStrategy::FullCostAverage, value_of_lost_load: MoneyPerFlow(1e9), diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 568767587..668f9ad2a 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -691,28 +691,27 @@ impl<'model, 'run> DispatchRun<'model, 'run> { } // Add constraints - let all_assets: Vec<_> = chain( + let all_assets = chain( existing_assets_with_parents.iter(), self.candidate_assets.iter(), - ) - .collect(); + ); let constraint_keys = add_model_constraints( &mut problem, &variables, self.model, - &all_assets.iter().copied(), + &all_assets, markets_to_balance, self.year, self.candidate_assets, ); - // Add secondary objective to encourage even dispatch across assets, if enabled + // Add secondary objective to encourage even dispatch across existing assets, if enabled let epsilon = self.model.parameters.dispatch_activity_equalisation_epsilon; if epsilon > 0.0 { add_activity_equalisation( &mut problem, &variables, - all_assets.iter().copied(), + existing_assets_with_parents.iter(), self.model, epsilon, ); @@ -776,7 +775,7 @@ fn add_activity_variables( /// Add secondary objective variables and constraints to encourage even dispatch across assets. /// -/// The secondary objective is `epsilon * mean(d)`, where each `d[a,t]` is the L1 deviation of +/// The secondary objective is `epsilon * sum(d)`, where each `d[a,t]` is the L1 deviation of /// asset `a`'s utilisation fraction from a single free centre variable `m`: /// /// ```text @@ -791,6 +790,9 @@ fn add_activity_variables( /// non-negative). The LP sets it to the L1-optimal centre (the median of utilisations), so the /// penalty measures spread rather than deviation from any fixed target. This makes `epsilon` /// independent of the absolute utilisation level, which is determined by the primary objective. +/// +/// Candidate assets are excluded: they exist only to receive shadow prices and have artificially +/// small capacity, so including them would distort `m` without meaningfully affecting dispatch. fn add_activity_equalisation<'a, I>( problem: &mut Problem, variables: &VariableMap, From c71cf114cd287e2636c45b18929c338130da35ed Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Thu, 30 Jul 2026 16:45:37 +0100 Subject: [PATCH 3/6] Improve docstring --- src/simulation/optimisation.rs | 46 +++++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 15 deletions(-) diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 668f9ad2a..b51f6ad0a 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -708,7 +708,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { // Add secondary objective to encourage even dispatch across existing assets, if enabled let epsilon = self.model.parameters.dispatch_activity_equalisation_epsilon; if epsilon > 0.0 { - add_activity_equalisation( + add_activity_regularisation( &mut problem, &variables, existing_assets_with_parents.iter(), @@ -773,27 +773,35 @@ fn add_activity_variables( start..problem.num_cols() } -/// Add secondary objective variables and constraints to encourage even dispatch across assets. +/// Add a regularisation term to the dispatch objective to discourage "corner" solutions. /// -/// The secondary objective is `epsilon * sum(d)`, where each `d[a,t]` is the L1 deviation of -/// asset `a`'s utilisation fraction from a single free centre variable `m`: +/// Without regularisation, when multiple assets have identical costs the LP may concentrate all +/// activity on a single asset rather than spreading it evenly — both solutions are equally optimal +/// from the primary objective's perspective, but the former is less physically realistic. +/// +/// The regularisation adds `epsilon * sum(d[a,t])` to the objective, where `d[a,t]` is a +/// variable representing the L1 deviation of asset `a`'s utilisation fraction from a free centre +/// variable `m` (shared across all assets and time slices): /// /// ```text /// d[a,t] = |act[a,t] / cap[a] - m| /// ``` /// -/// `|x|` is linearised via two constraints per `d`: `d >= x` and `d >= -x` (both with `d >= 0`). -/// Both are needed because we minimise `d` — without the second constraint HiGHS would set `d = 0` -/// whenever `x < 0`, missing half the deviation. +/// The absolute value is linearised as two inequality constraints per `(a, t)`: +/// +/// ```text +/// d[a,t] >= act[a,t] / cap[a] - m +/// d[a,t] >= m - act[a,t] / cap[a] +/// ``` /// -/// `m` is a free variable (bounded below at zero since all utilisation fractions are -/// non-negative). The LP sets it to the L1-optimal centre (the median of utilisations), so the -/// penalty measures spread rather than deviation from any fixed target. This makes `epsilon` -/// independent of the absolute utilisation level, which is determined by the primary objective. +/// `m` is optimised by the LP and settles at the median utilisation across all assets and time +/// slices, so the regularisation penalises *spread* rather than deviation from any fixed target. +/// This means it acts as a pure tiebreaker: `epsilon` should be chosen small enough that no +/// change to the primary objective cost is acceptable as a trade for reduced spread. /// -/// Candidate assets are excluded: they exist only to receive shadow prices and have artificially -/// small capacity, so including them would distort `m` without meaningfully affecting dispatch. -fn add_activity_equalisation<'a, I>( +/// Candidate assets are excluded because they have artificially small capacity and exist only to +/// receive shadow prices, not to make real dispatch decisions. +fn add_activity_regularisation<'a, I>( problem: &mut Problem, variables: &VariableMap, assets: I, @@ -807,9 +815,17 @@ fn add_activity_equalisation<'a, I>( let m = problem.add_column(0.0, 0.0..); for asset in assets { - // Flexible capacity assets use their current capacity as an approximation + // Get capacity for the asset, since we're regularising utilisation (activity/capacity). + // + // Note: for flexible capacity assets this is their _initial_ capacity, which may end up + // different from their final capacity in the solution - therefore, this is an + // approximation, but shouldn't be too far off since final capacities are bound within a + // tolerance of the initial capacity. We can't use the actual capacity variable here because + // that would make the problem non-linear. let cap = asset.total_capacity().value(); if cap <= 1e-9 { + // Skip very small capacities to avoid numerical issues. Unlikely to ever have assets + // this small in practice. continue; } let inv_cap = 1.0 / cap; From 0fc0754a62694ddce3d330bf4ac52a97299392e5 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 31 Jul 2026 15:28:51 +0100 Subject: [PATCH 4/6] Exclude pinned activity variables from regularisation --- src/simulation/optimisation.rs | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index b51f6ad0a..bea611aaf 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -831,6 +831,12 @@ fn add_activity_regularisation<'a, I>( let inv_cap = 1.0 / cap; for time_slice in model.time_slice_info.iter_ids() { + let limits = asset.get_activity_per_capacity_limits(time_slice); + // Skip variables with no freedom — they distort m without contributing to spreading + if limits.start() == limits.end() { + continue; + } + let act = variables.get_activity_var(asset, time_slice); let d = problem.add_column(epsilon, 0.0..); // d >= act/cap - m From 0ea8f0270dbc57d96bf20d7999e20b5a0824b2f2 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 31 Jul 2026 15:33:29 +0100 Subject: [PATCH 5/6] Add limitations comment --- src/simulation/optimisation.rs | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index bea611aaf..2ece37cf9 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -801,6 +801,16 @@ fn add_activity_variables( /// /// Candidate assets are excluded because they have artificially small capacity and exist only to /// receive shadow prices, not to make real dispatch decisions. +/// +/// # Limitations +/// +/// `m` settles at the median utilisation, so the regularisation is only effective when fewer than +/// half of the included `(asset, timeslice)` pairs are at zero. If the majority are zero (e.g. +/// many off-peak timeslices where most assets are idle), the median collapses to zero and the +/// tiebreaking effect for the active assets is lost. Assets with hard-constrained activity limits +/// (`min == max`) are excluded from the regularisation to avoid distorting `m`; cost-driven zeros +/// (assets at zero due to high cost or absent demand) are not excluded and contribute to this +/// limitation. fn add_activity_regularisation<'a, I>( problem: &mut Problem, variables: &VariableMap, From 5d0a172c5fb9ae69ea9ef539a0ce975a562d7170 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 31 Jul 2026 18:01:49 +0100 Subject: [PATCH 6/6] Add test for regularisation (currently failing) --- src/simulation/optimisation.rs | 69 ++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 2ece37cf9..096ac97e4 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -964,3 +964,72 @@ fn calculate_capacity_coefficient(asset: &AssetRef) -> MoneyPerCapacity { annual_fixed_operating_cost + annual_capital_cost(param.capital_cost, param.lifetime, param.discount_rate) } + +#[cfg(test)] +mod tests { + use crate::patch::{FilePatch, ModelPatch}; + use indexmap::IndexMap; + use serde::Deserialize; + use std::fs; + + #[derive(Debug, Deserialize)] + struct DispatchRow { + milestone_year: u32, + run_description: String, + process_id: String, + time_slice: String, + activity: Option, + } + + // Verifies that two identical assets receive equal dispatch under L1 regularisation. + // With epsilon=0 this test is expected to fail (LP corner solutions). + #[test] + fn two_identical_assets_dispatch_evenly() { + let tmp = ModelPatch::from_example("simple") + .with_file_patches([ + FilePatch::new("assets.csv").with_addition("GASCGT,GBR,A0_ELC,2.430,2020") + ]) + .build_to_tempdir() + .unwrap(); + + let output_path = tmp.path().join("output"); + fs::create_dir_all(&output_path).unwrap(); + + let model = crate::input::load_model(tmp.path()).unwrap(); + crate::simulation::run(&model, &output_path, /*debug_model=*/ true).unwrap(); + + // Read per-asset activity for GASCGT, grouped by time slice + let mut activities: IndexMap> = IndexMap::new(); + let mut reader = + csv::Reader::from_path(output_path.join("debug_dispatch_assets.csv")).unwrap(); + for result in reader.deserialize::() { + let row = result.unwrap(); + if row.milestone_year == 2020 + && row.run_description == "final without candidates" + && row.process_id == "GASCGT" + && let Some(act) = row.activity + { + activities.entry(row.time_slice).or_default().push(act); + } + } + + // Every time slice where at least one GASCGT is active should have equal activities + for (time_slice, acts) in &activities { + assert_eq!( + acts.len(), + 2, + "Expected 2 GASCGT entries for time slice {time_slice}" + ); + let [a, b] = acts.as_slice() else { + unreachable!() + }; + if *a > 0.0 || *b > 0.0 { + let diff = (a - b).abs() / a.max(*b); + assert!( + diff < 1e-6, + "GASCGT activities not equal in time slice {time_slice}: {a} vs {b}" + ); + } + } + } +}