Small penalties to prevent volatile asset utilisation

docs/model/dispatch_optimisation.md

@@ -112,6 +112,63 @@ satisfying an additional unit of demand for that commodity in region \\( r \\) d
 
 ---
 
+## Seasonal/Annual Utilisation Penalties
+
+MUSE2 optionally applies small penalties to the peak capacity required by each asset within a season
+and across the whole year. These penalties encourage activity to be distributed across time slices
+within a season and across seasons, respectively. They are particularly useful when commodities are
+balanced at the seasonal/annual levels and the balance constraint otherwise leaves the
+intra-seasonal or inter-seasonal production profile undetermined. In real-world terms, this
+represents a preference to avoid concentrating an asset's operation into short periods of high
+utilisation, which may reduce cycling, wear, start-up requirements, or the need to maintain capacity
+for seasonal peaks.
+
+For asset \\(a\\) and time slice \\(t\\), the capacity required to support its activity is
+
+\\[
+  \mathrm{RequiredCapacity}\_{a,t} =
+  \frac{\mathrm{Activity}\_{a,t}}
+  {\mathrm{cap2act}\_a \cdot \Delta\_t}
+\\]
+
+For each asset and season, MUSE2 introduces an auxiliary variable \\(U_{a,s}\\), representing the
+peak capacity required by the asset during that season. It is constrained by
+
+\\[
+  U_{a,s} \geq \mathrm{RequiredCapacity}_{a,t}
+\\]
+
+for every time slice \\(t\\) in season \\(s\\). In addition, MUSE2 introduces an auxiliary variable
+\\(U_{a,\\mathrm{annual}}\\), representing the greatest seasonal peak capacity required by the asset
+during the year. It is constrained by
+
+\\[
+  U_{a,\\mathrm{annual}} \\geq U_{a,s}
+\\]
+
+for every season \\(s\\).
+
+When enabled, the penalties add the following term to the optimisation objective:
+
+\\[
+  \\lambda_{\\mathrm{seasonal}}
+  \\sum_{a \\in \\mathbf{A}} \\sum_{s \\in \\mathbf{S}}
+  \\Delta_s U_{a,s}
+  +
+  \\lambda_{\\mathrm{annual}}
+  \\sum_{a \\in \\mathbf{A}} U_{a,\\mathrm{annual}}
+\\]
+
+Here, \\(\\lambda_{\\mathrm{seasonal}}\\) and \\(\\lambda_{\\mathrm{annual}}\\) are set by the
+`seasonal_utilisation_penalty` and `annual_utilisation_penalty` model parameters, respectively.
+Setting either parameter to zero disables its corresponding penalty. The seasonal parameter controls
+how strongly activity is spread within seasons, while the annual parameter controls how strongly it
+is spread across seasons. Both are weighted objective terms, so their values should be small enough
+that smoothing dispatch does not outweigh meaningful differences in operating cost (default for
+both = `1e-6`).
+
+---
+
 ## Candidate Dispatch Run
 
 After the primary dispatch run, MUSE2 performs a second dispatch run that includes

schemas/input/model.yaml

@@ -56,6 +56,22 @@ properties:
     default: 1e9
     notes: |
       Currently this only applies to the LCOX appraisal.
+  seasonal_utilisation_penalty:
+    type: number
+    description: Penalty per unit of capacity used within a season
+    default: 1e-6
+    notes: |
+      A small additive penalty to discourage assets from using a high capacity-equivalent peak
+      within a season. This should be very small, as it is only intended to break ties between
+      otherwise equivalent solutions.
+  annual_utilisation_penalty:
+    type: number
+    description: Penalty per unit of capacity used across the whole year
+    default: 1e-6
+    notes: |
+      A small additive penalty to discourage assets from using a high capacity-equivalent peak
+      across the year. This should be very small, as it is only intended to break ties between
+      otherwise equivalent solutions.
   max_ironing_out_iterations:
     type: integer
     description: The maximum number of iterations to run the "ironing out" step of agent investment for

src/model/parameters.rs

@@ -9,7 +9,7 @@ use crate::input::{
     deserialise_finite_non_negative, deserialise_proportion_nonzero, input_err_msg,
     is_sorted_and_unique, read_toml,
 };
-use crate::units::{Capacity, Dimensionless, Flow, MoneyPerFlow};
+use crate::units::{Capacity, Dimensionless, Flow, MoneyPerCapacityPerYear, MoneyPerFlow};
 use anyhow::{Context, Result, ensure};
 use itertools::Itertools;
 use log::warn;
@@ -99,6 +99,12 @@ pub struct ModelParameters {
     ///
     /// Currently this only applies to the LCOX appraisal.
     pub value_of_lost_load: MoneyPerFlow,
+    /// Additive penalty per unit of capacity used within a season.
+    #[serde(deserialize_with = "deserialise_finite_non_negative")]
+    pub seasonal_utilisation_penalty: MoneyPerCapacityPerYear,
+    /// Additive penalty per unit of capacity used across the whole year.
+    #[serde(deserialize_with = "deserialise_finite_non_negative")]
+    pub annual_utilisation_penalty: MoneyPerCapacityPerYear,
     /// The maximum number of iterations to run the "ironing out" step of agent investment for
     pub max_ironing_out_iterations: u32,
     /// The relative tolerance for price convergence in the ironing out loop
@@ -138,6 +144,8 @@ impl Default for ModelParameters {
             capacity_limit_factor: Dimensionless(0.05),
             fallback_pricing_strategy: PricingStrategy::FullCostAverage,
             value_of_lost_load: MoneyPerFlow(1e9),
+            seasonal_utilisation_penalty: MoneyPerCapacityPerYear(1e-6),
+            annual_utilisation_penalty: MoneyPerCapacityPerYear(1e-6),
             max_ironing_out_iterations: 1,
             price_tolerance: Dimensionless(1e-6),
             capacity_margin: Dimensionless(0.2),

src/simulation/optimisation/constraints.rs

@@ -4,8 +4,8 @@ use crate::asset::{AssetCapacity, AssetIterator, AssetRef};
 use crate::commodity::{CommodityID, CommodityType};
 use crate::model::Model;
 use crate::region::RegionID;
-use crate::time_slice::{TimeSliceInfo, TimeSliceSelection};
-use crate::units::{Flow, UnitType};
+use crate::time_slice::{Season, TimeSliceInfo, TimeSliceSelection};
+use crate::units::{Flow, MoneyPerCapacityPerYear, UnitType, Year};
 use highs::RowProblem as Problem;
 use indexmap::IndexMap;
 
@@ -48,6 +48,12 @@ pub type CommodityBalanceKeys = KeysWithOffset<(CommodityID, RegionID, TimeSlice
 /// Indicates the asset ID and time slice covered by each activity constraint
 pub type ActivityKeys = KeysWithOffset<(AssetRef, TimeSliceSelection)>;
 
+/// Map containing the seasonal peak variables for each (asset, season) pair
+type SeasonalPeakVariableMap = IndexMap<(AssetRef, Season), highs::Col>;
+
+/// Map containing the annual peak variables for each asset
+type AnnualPeakVariableMap = IndexMap<AssetRef, highs::Col>;
+
 /// The keys for different constraints
 pub struct ConstraintKeys {
     /// Keys for commodity balance constraints
@@ -98,13 +104,147 @@ where
     let activity_keys =
         add_activity_constraints(problem, variables, &model.time_slice_info, assets.clone());
 
+    add_utilisation_peak_constraints(problem, model, assets.clone(), variables);
+
     // Return constraint keys
     ConstraintKeys {
         commodity_balance_keys,
         activity_keys,
     }
 }
 
+/// Add seasonal and annual utilisation peak constraints to the problem.
+fn add_utilisation_peak_constraints<'a, I>(
+    problem: &mut Problem,
+    model: &Model,
+    assets: I,
+    variables: &VariableMap,
+) where
+    I: Iterator<Item = &'a AssetRef> + Clone,
+{
+    let has_seasonal_penalty =
+        model.parameters.seasonal_utilisation_penalty > MoneyPerCapacityPerYear(0.0);
+    let has_annual_penalty =
+        model.parameters.annual_utilisation_penalty > MoneyPerCapacityPerYear(0.0);
+
+    // If neither penalties are applied, we don't need to add any variables and constraints
+    if !has_seasonal_penalty && !has_annual_penalty {
+        return;
+    }
+
+    // So long as either penalty is applied, we need to add seasonal peak variables and constraints
+    let seasonal_peak_vars = add_seasonal_peak_variables(problem, model, assets.clone());
+    add_seasonal_peak_constraints(
+        problem,
+        variables,
+        &model.time_slice_info,
+        &seasonal_peak_vars,
+    );
+
+    // If the annual penalty is applied, we also need to add annual peak variables and constraints
+    if has_annual_penalty {
+        let annual_peak_vars = add_annual_peak_variables(problem, model, assets);
+        add_annual_peak_constraints(
+            problem,
+            &model.time_slice_info,
+            &annual_peak_vars,
+            &seasonal_peak_vars,
+        );
+    }
+}
+
+/// Add seasonal peak variables to the problem for each (asset, season) pair.
+fn add_seasonal_peak_variables<'a, I>(
+    problem: &mut Problem,
+    model: &Model,
+    assets: I,
+) -> SeasonalPeakVariableMap
+where
+    I: Iterator<Item = &'a AssetRef>,
+{
+    let mut seasonal_peak_vars = SeasonalPeakVariableMap::new();
+    for asset in assets {
+        for (season, duration) in &model.time_slice_info.seasons {
+            // Scale penalty by season duration
+            let col_factor = model.parameters.seasonal_utilisation_penalty * *duration;
+            let variable = problem.add_column(col_factor.value(), 0.0..);
+            seasonal_peak_vars.insert((asset.clone(), season.clone()), variable);
+        }
+    }
+    seasonal_peak_vars
+}
+
+/// Add annual peak variables to the problem for each asset.
+fn add_annual_peak_variables<'a, I>(
+    problem: &mut Problem,
+    model: &Model,
+    assets: I,
+) -> AnnualPeakVariableMap
+where
+    I: Iterator<Item = &'a AssetRef>,
+{
+    // Penalty is applied over the whole year, so scale by 1 year
+    let col_factor = model.parameters.annual_utilisation_penalty * Year(1.0);
+    assets
+        .map(|asset| {
+            let variable = problem.add_column(col_factor.value(), 0.0..);
+            (asset.clone(), variable)
+        })
+        .collect()
+}
+
+/// Add constraints linking seasonal peak variables to activity variables for each (asset, season) pair.
+fn add_seasonal_peak_constraints(
+    problem: &mut Problem,
+    variables: &VariableMap,
+    time_slice_info: &TimeSliceInfo,
+    seasonal_peak_vars: &SeasonalPeakVariableMap,
+) {
+    for ((asset, season), &peak_variable) in seasonal_peak_vars {
+        let activity_per_capacity = asset.process().capacity_to_activity;
+        let season_selection = TimeSliceSelection::Season(season.clone());
+        for (time_slice, ts_length) in season_selection.iter(time_slice_info) {
+            let time_slice_fraction = ts_length / Year(1.0);
+            let activity_per_capacity_in_time_slice = activity_per_capacity * time_slice_fraction;
+            let capacity_required_per_activity = 1.0 / activity_per_capacity_in_time_slice.value();
+
+            // One unit of capacity supports `activity_per_capacity_in_time_slice` activity in
+            // this time slice. The peak variable therefore measures the capacity required by
+            // the activity in the time slice.
+            problem.add_row(
+                0.0..,
+                [
+                    (peak_variable, 1.0),
+                    (
+                        variables.get_activity_var(asset, time_slice),
+                        -capacity_required_per_activity,
+                    ),
+                ],
+            );
+        }
+    }
+}
+
+/// Add constraints linking seasonal peak variables to annual peak variables for each asset.
+fn add_annual_peak_constraints(
+    problem: &mut Problem,
+    time_slice_info: &TimeSliceInfo,
+    annual_peak_vars: &AnnualPeakVariableMap,
+    seasonal_peak_vars: &SeasonalPeakVariableMap,
+) {
+    for (asset, &annual_peak_variable) in annual_peak_vars {
+        for season in time_slice_info.seasons.keys() {
+            let seasonal_peak_variable = seasonal_peak_vars
+                .get(&(asset.clone(), season.clone()))
+                .expect("Missing seasonal peak variable for annual peak constraint");
+            problem.add_row(
+                0.0..,
+                [(annual_peak_variable, 1.0), (*seasonal_peak_variable, -1.0)],
+            );
+        }
+    }
+}
+
 /// Add asset-level input-output commodity balances.
 ///
 /// These constraints fix the supply-demand balance for the whole system.

tests/data/circularity/asset_capacities.csv

@@ -24,23 +24,23 @@ milestone_year,asset_id,capacity,num_units
 2030,15,777.3007440837426,
 2030,16,1849.2839907535792,
 2030,17,3680.183145405078,
-2030,18,113.37276042798105,
-2030,19,2944.146516324062,
-2030,20,119.46854244531636,
+2030,18,298.6475820926018,
+2030,19,2944.1465163240623,
+2030,20,163.3092156793555,
 2040,1,1738.05,
 2040,5,3.964844,
 2040,6,2.999,
 2040,16,1849.2839907535792,
-2040,18,113.37276042798105,
-2040,20,119.46854244531636,
+2040,18,298.6475820926018,
+2040,20,163.3092156793555,
 2040,21,912.8939641298446,
 2040,22,2162.373384722901,
 2040,23,33.51213708228713,
-2040,24,4.285571618385969,
-2040,25,9.849164928608296,
-2040,26,115.45981236997261,
-2040,27,462.124901316292,
+2040,24,5.471034970878489,
+2040,25,8.663701576115777,
+2040,26,151.10373493751322,
+2040,27,481.3573712526796,
 2040,28,972.4856516656134,
-2040,29,2492.1031454342215,
-2040,30,485.2311463821062,
-2040,31,2627.922766860386,
+2040,29,2349.697251409409,
+2040,30,505.4252398153141,
+2040,31,2592.9322599984084,