Skip to main content

muse2/simulation/investment/appraisal/
optimisation.rs

1//! Optimisation problem for investment tools.
2use super::DemandMap;
3use super::ObjectiveCoefficients;
4use super::constraints::{
5    add_activity_constraints, add_demand_constraints, add_utilisation_peak_constraints,
6};
7use crate::asset::AssetRef;
8use crate::commodity::Commodity;
9use crate::model::Model;
10use crate::simulation::optimisation::ModelError;
11use crate::simulation::optimisation::apply_highs_options_from_toml;
12use crate::simulation::optimisation::solve_optimal;
13use crate::time_slice::{TimeSliceID, TimeSliceInfo};
14use crate::units::{Activity, Flow};
15use anyhow::{Context, Result};
16use highs::{RowProblem as Problem, Sense};
17use indexmap::IndexMap;
18
19/// A decision variable in the optimisation
20///
21/// This alias represents a column created in the `highs` solver. Callers rely on the order
22/// in which columns are added to the problem when extracting solution values.
23pub type Variable = highs::Col;
24
25/// Map containing optimisation results and coefficients
26pub struct ResultsMap {
27    /// Activity variables in each time slice
28    pub activity: IndexMap<TimeSliceID, Activity>,
29    /// Remaining unmet demand per time slice, computed post-solve
30    pub unmet_demand: DemandMap,
31}
32
33/// Adds activity variables to the problem, one per time slice.
34///
35/// Returns a map from time slice to the corresponding decision variable.
36fn add_activity_vars(
37    problem: &mut Problem,
38    cost_coefficients: &ObjectiveCoefficients,
39) -> IndexMap<TimeSliceID, Variable> {
40    cost_coefficients
41        .activity_coefficients
42        .iter()
43        .map(|(time_slice, cost)| {
44            let var = problem.add_column(cost.value(), 0.0..);
45            (time_slice.clone(), var)
46        })
47        .collect()
48}
49
50/// Adds constraints to the problem.
51fn add_constraints(
52    problem: &mut Problem,
53    model: &Model,
54    asset: &AssetRef,
55    commodity: &Commodity,
56    activity_vars: &IndexMap<TimeSliceID, Variable>,
57    demand: &DemandMap,
58    time_slice_info: &TimeSliceInfo,
59) {
60    add_activity_constraints(problem, asset, activity_vars, time_slice_info);
61    add_demand_constraints(
62        problem,
63        asset,
64        commodity,
65        time_slice_info,
66        demand,
67        activity_vars,
68    );
69    add_utilisation_peak_constraints(problem, model, asset, activity_vars);
70}
71
72/// Computes remaining unmet demand per time-slice selection after a solve.
73///
74/// For each time-slice selection at the commodity's balance level, the selection-level residual
75/// (`demand - supply`, clamped to zero) is stored directly. Downstream consumers operate at the
76/// selection level (e.g. the next round's demand constraint and `is_any_remaining_demand`), so
77/// there is no need to distribute values across individual time slices.
78fn compute_unmet_demand(
79    demand: &DemandMap,
80    activity: &IndexMap<TimeSliceID, Activity>,
81    commodity: &Commodity,
82    asset: &AssetRef,
83    time_slice_info: &TimeSliceInfo,
84) -> DemandMap {
85    let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff;
86    time_slice_info
87        .iter_selections_at_level(commodity.time_slice_level)
88        .map(|ts_selection| {
89            let demand_for_selection = demand[&ts_selection];
90            let supply_for_selection: Flow = ts_selection
91                .iter(time_slice_info)
92                .map(|(ts, _)| activity[ts] * flow_coeff)
93                .sum();
94            let unmet = (demand_for_selection - supply_for_selection).max(Flow(0.0));
95            (ts_selection, unmet)
96        })
97        .collect()
98}
99
100/// Performs optimisation for an asset, given the coefficients and demand.
101pub fn perform_optimisation(
102    model: &Model,
103    asset: &AssetRef,
104    commodity: &Commodity,
105    coefficients: &ObjectiveCoefficients,
106    demand: &DemandMap,
107) -> Result<ResultsMap> {
108    // Create problem and add variables
109    let mut problem = Problem::default();
110    let activity_vars = add_activity_vars(&mut problem, coefficients);
111
112    // Add constraints
113    add_constraints(
114        &mut problem,
115        model,
116        asset,
117        commodity,
118        &activity_vars,
119        demand,
120        &model.time_slice_info,
121    );
122
123    // Solve model
124    let mut highs_model = problem.optimise(Sense::Maximise);
125    apply_highs_options_from_toml(&mut highs_model, &model.parameters.highs.appraisal_options)
126        .context("Failed to apply custom HiGHS options to appraisal optimisation")?;
127    // Enforce single-threaded solving: when appraisals run in parallel via Rayon each HiGHS
128    // instance must not spawn additional worker threads. Setting `parallel="off"` disables
129    // HiGHS's concurrent simplex strategy without touching the global thread-pool scheduler
130    // (setting `threads=N` would fail if the scheduler was already initialised on this thread
131    // by a previous solve with a different count).
132    highs_model.set_option("parallel", "off");
133    let solution = solve_optimal(highs_model)
134        .map_err(ModelError::into_anyhow)?
135        .get_solution();
136    let solution_values = solution.columns();
137    let activity = activity_vars
138        .keys()
139        .zip(solution_values.iter())
140        .map(|(time_slice, &value)| (time_slice.clone(), Activity::new(value)))
141        .collect();
142    let unmet_demand =
143        compute_unmet_demand(demand, &activity, commodity, asset, &model.time_slice_info);
144    Ok(ResultsMap {
145        activity,
146        unmet_demand,
147    })
148}