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..610a7a169 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: 1e-6, 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..096ac97e4 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -705,6 +705,18 @@ impl<'model, 'run> DispatchRun<'model, 'run> { self.candidate_assets, ); + // 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_regularisation( + &mut problem, + &variables, + existing_assets_with_parents.iter(), + 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 +773,90 @@ fn add_activity_variables( start..problem.num_cols() } +/// Add a regularisation term to the dispatch objective to discourage "corner" solutions. +/// +/// 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| +/// ``` +/// +/// 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 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 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, + 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 { + // 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; + + 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 + 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, @@ -868,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}" + ); + } + } + } +}