Add L1 regularisation to dispatch

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

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)?;
 

src/simulation/optimisation.rs

@@ -659,7 +659,7 @@
             self.model,
             input_prices,
             &existing_assets_with_parents,
             self.candidate_assets,
             self.year,
         );
 
@@ -705,6 +705,18 @@
             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 @@
     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<Item = &'a AssetRef>,
+{
+    // 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 @@
     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<f64>,
+    }
+
+    // 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<String, Vec<f64>> = IndexMap::new();
+        let mut reader =
+            csv::Reader::from_path(output_path.join("debug_dispatch_assets.csv")).unwrap();
+        for result in reader.deserialize::<DispatchRow>() {
+            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}"
+                );
+            }
+        }
+    }
+}