Add seasonal/annual utilisation penalties to the appraisal mini-dispatch

docs/model/investment.md

@@ -256,6 +256,11 @@ constraints:
   \max \sum_t \alpha_t \cdot \mathrm{Activity}_t
 \\]
 
+In addition, seasonal/annual utilisation penalties are included to encourage even spread of
+utilisation within and across seasons. These are equivalent to the
+[Seasonal/Annual Utilisation Penalties][dispatch-utilisation-penalties] used in the full-system
+dispatch optimisation.
+
 ## Metric Calculation
 
 After the dispatch LP is solved, an investment metric is calculated from the resulting activity
@@ -518,3 +523,4 @@ assets shaping the demand seen by those upstream.
 [model-toml]: ../file_formats/input_files.md#model-parameters-modeltoml
 [process-investment-constraints-csv]: ../file_formats/input_files.md#process_investment_constraintscsv
 [dispatch-optimisation]: ./dispatch_optimisation.md
+[dispatch-utilisation-penalties]: ./dispatch_optimisation.md#seasonalannual-utilisation-penalties

docs/release_notes/upcoming.md

@@ -2,6 +2,10 @@
 
 ## New features
 
+- Seasonal/annual utilisation penalties are now included in the investment mini-dispatch ([#1548])
+
 ## Breaking changes
 
 ## Bug fixes
+
+[#1548]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1548

src/simulation/investment/appraisal/constraints.rs

@@ -3,8 +3,9 @@ use super::DemandMap;
 use super::optimisation::Variable;
 use crate::asset::AssetRef;
 use crate::commodity::Commodity;
-use crate::time_slice::{TimeSliceID, TimeSliceInfo};
-use crate::units::Flow;
+use crate::model::Model;
+use crate::time_slice::{Season, TimeSliceID, TimeSliceInfo, TimeSliceSelection};
+use crate::units::{Flow, MoneyPerCapacityPerYear, Year};
 use highs::RowProblem as Problem;
 use indexmap::IndexMap;
 
@@ -27,7 +28,7 @@ pub fn add_activity_constraints(
         // Collect activity terms for the time slices in this selection
         let terms = ts_selection
             .iter(time_slice_info)
-            .map(|(time_slice, _)| (*activity_vars.get(time_slice).unwrap(), 1.0))
+            .map(|(time_slice, _)| (activity_vars[time_slice], 1.0))
             .collect::<Vec<_>>();
 
         // Constraint: sum of activities in selection within limits
@@ -48,14 +49,127 @@ pub fn add_demand_constraints(
     demand: &DemandMap,
     activity_vars: &IndexMap<TimeSliceID, Variable>,
 ) {
+    let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff;
     for ts_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) {
         let mut demand_for_ts_selection = Flow(0.0);
         let mut terms = Vec::new();
         for (time_slice, _) in ts_selection.iter(time_slice_info) {
             demand_for_ts_selection += demand[time_slice];
-            let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff;
             terms.push((activity_vars[time_slice], flow_coeff.value()));
         }
         problem.add_row(0.0..=demand_for_ts_selection.value(), terms);
     }
 }
+
+/// Add seasonal and annual utilisation peak constraints to the problem.
+///
+/// This is almost identical to equivalent constraints in the full-system dispatch, but here we
+/// make the penalties negative as this is a maximise optimisation (whereas dispatch is a minimise
+/// optimisation).
+pub fn add_utilisation_peak_constraints(
+    problem: &mut Problem,
+    model: &Model,
+    asset: &AssetRef,
+    activity_vars: &IndexMap<TimeSliceID, Variable>,
+) {
+    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);
+    add_seasonal_peak_constraints(
+        problem,
+        asset,
+        activity_vars,
+        &model.time_slice_info,
+        &seasonal_peak_vars,
+    );
+
+    // If the annual penalty is applied, we also need to add an annual peak variable and constraints
+    if has_annual_penalty {
+        let annual_peak_var = add_annual_peak_variable(problem, model);
+        add_annual_peak_constraints(
+            problem,
+            &model.time_slice_info,
+            annual_peak_var,
+            &seasonal_peak_vars,
+        );
+    }
+}
+
+/// Add seasonal peak variables to the problem.
+fn add_seasonal_peak_variables(
+    problem: &mut Problem,
+    model: &Model,
+) -> IndexMap<Season, highs::Col> {
+    let mut seasonal_peak_vars = IndexMap::new();
+    for (season, duration) in &model.time_slice_info.seasons {
+        // Scale penalty by season duration
+        // We make the penalty negative as we're doing a max optimisation
+        let col_factor = -(model.parameters.seasonal_utilisation_penalty * *duration);
+        let variable = problem.add_column(col_factor.value(), 0.0..);
+        seasonal_peak_vars.insert(season.clone(), variable);
+    }
+    seasonal_peak_vars
+}
+
+/// Add annual peak variable to the problem.
+fn add_annual_peak_variable(problem: &mut Problem, model: &Model) -> highs::Col {
+    // Penalty is applied over the whole year, so scale by 1 year
+    // We make the penalty negative as we're doing a max optimisation
+    let col_factor = -(model.parameters.annual_utilisation_penalty * Year(1.0));
+    problem.add_column(col_factor.value(), 0.0..)
+}
+
+/// Add constraints linking seasonal peak variables to activity variables for each (asset, season) pair.
+fn add_seasonal_peak_constraints(
+    problem: &mut Problem,
+    asset: &AssetRef,
+    activity_vars: &IndexMap<TimeSliceID, Variable>,
+    time_slice_info: &TimeSliceInfo,
+    seasonal_peak_vars: &IndexMap<Season, highs::Col>,
+) {
+    let activity_per_capacity = asset.process().capacity_to_activity;
+    for (season, &peak_variable) in seasonal_peak_vars {
+        for (time_slice, ts_length) in
+            TimeSliceSelection::Season(season.clone()).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),
+                    (activity_vars[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_var: highs::Col,
+    seasonal_peak_vars: &IndexMap<Season, highs::Col>,
+) {
+    for season in time_slice_info.seasons.keys() {
+        problem.add_row(
+            0.0..,
+            [(annual_peak_var, 1.0), (seasonal_peak_vars[season], -1.0)],
+        );
+    }
+}

src/simulation/investment/appraisal/optimisation.rs

@@ -1,7 +1,9 @@
 //! Optimisation problem for investment tools.
 use super::DemandMap;
 use super::ObjectiveCoefficients;
-use super::constraints::{add_activity_constraints, add_demand_constraints};
+use super::constraints::{
+    add_activity_constraints, add_demand_constraints, add_utilisation_peak_constraints,
+};
 use crate::asset::AssetRef;
 use crate::commodity::Commodity;
 use crate::model::Model;
@@ -48,6 +50,7 @@ fn add_activity_vars(
 /// Adds constraints to the problem.
 fn add_constraints(
     problem: &mut Problem,
+    model: &Model,
     asset: &AssetRef,
     commodity: &Commodity,
     activity_vars: &IndexMap<TimeSliceID, Variable>,
@@ -63,6 +66,7 @@ fn add_constraints(
         demand,
         activity_vars,
     );
+    add_utilisation_peak_constraints(problem, model, asset, activity_vars);
 }
 
 /// Computes remaining unmet demand per time slice after a solve.
@@ -116,6 +120,7 @@ pub fn perform_optimisation(
     // Add constraints
     add_constraints(
         &mut problem,
+        model,
         asset,
         commodity,
         &activity_vars,

tests/data/circularity/asset_capacities.csv

@@ -23,15 +23,15 @@ milestone_year,asset_id,capacity,num_tranches,mothballed_capacity,mothballed_tra
 2030,14,399.98,1,0.0,0
 2030,15,802.563016789748,28,0.0,0
 2030,16,1849.28399075358,36,0.0,0
-2030,17,3680.183145405077,26,0.0,0
+2030,17,2830.9101118500594,20,0.0,0
 2030,18,298.6475820926018,57,0.0,0
 2030,19,2944.1465163240614,16,0.0,0
 2030,20,165.3505808753475,27,0.0,0
 2040,5,3.964844,1,0.0,0
 2040,6,2.999,1,0.0,0
 2040,16,1849.28399075358,36,0.0,0
 2040,18,298.6475820926018,57,0.0,0
-2040,19,1288.064100891777,7,0.0,0
+2040,19,1104.054943621523,6,0.0,0
 2040,20,165.3505808753475,27,0.0,0
 2040,21,942.563016229748,28,0.0,0
 2040,22,2198.390989008045,39,0.0,0