muse2/simulation/investment/appraisal/
optimisation.rs

1//! Optimisation problem for investment tools.
2use super::DemandMap;
3use super::ObjectiveCoefficients;
4use super::constraints::{add_activity_constraints, add_demand_constraints};
5use crate::asset::AssetRef;
6use crate::commodity::Commodity;
7use crate::model::Model;
8use crate::simulation::optimisation::ModelError;
9use crate::simulation::optimisation::apply_highs_options_from_toml;
10use crate::simulation::optimisation::solve_optimal;
11use crate::time_slice::{TimeSliceID, TimeSliceInfo};
12use crate::units::{Activity, Dimensionless, Flow};
13use anyhow::{Context, Result};
14use highs::{RowProblem as Problem, Sense};
15use indexmap::IndexMap;
16
17/// A decision variable in the optimisation
18///
19/// This alias represents a column created in the `highs` solver. Callers rely on the order
20/// in which columns are added to the problem when extracting solution values.
21pub type Variable = highs::Col;
22
23/// Map containing optimisation results and coefficients
24pub struct ResultsMap {
25    /// Activity variables in each time slice
26    pub activity: IndexMap<TimeSliceID, Activity>,
27    /// Remaining unmet demand per time slice, computed post-solve
28    pub unmet_demand: DemandMap,
29}
30
31/// Adds activity variables to the problem, one per time slice.
32///
33/// Returns a map from time slice to the corresponding decision variable.
34fn add_activity_vars(
35    problem: &mut Problem,
36    cost_coefficients: &ObjectiveCoefficients,
37) -> IndexMap<TimeSliceID, Variable> {
38    cost_coefficients
39        .activity_coefficients
40        .iter()
41        .map(|(time_slice, cost)| {
42            let var = problem.add_column(cost.value(), 0.0..);
43            (time_slice.clone(), var)
44        })
45        .collect()
46}
47
48/// Adds constraints to the problem.
49fn add_constraints(
50    problem: &mut Problem,
51    asset: &AssetRef,
52    commodity: &Commodity,
53    activity_vars: &IndexMap<TimeSliceID, Variable>,
54    demand: &DemandMap,
55    time_slice_info: &TimeSliceInfo,
56) {
57    add_activity_constraints(problem, asset, activity_vars, time_slice_info);
58    add_demand_constraints(
59        problem,
60        asset,
61        commodity,
62        time_slice_info,
63        demand,
64        activity_vars,
65    );
66}
67
68/// Computes remaining unmet demand per time slice after a solve.
69///
70/// For each time-slice selection at the commodity's balance level, the selection-level residual
71/// (`demand_total - supply_total`, clamped to zero) is divided equally across the time slices in
72/// the selection.
73///
74/// The exact per-time-slice distribution is arbitrary: all downstream consumers sum values back up
75/// to the selection level before using them (e.g. the next round's demand constraint, and
76/// `is_any_remaining_demand`), so only the selection-level total needs to be correct.
77fn compute_unmet_demand(
78    demand: &DemandMap,
79    activity: &IndexMap<TimeSliceID, Activity>,
80    commodity: &Commodity,
81    asset: &AssetRef,
82    time_slice_info: &TimeSliceInfo,
83) -> DemandMap {
84    let mut unmet_demand = DemandMap::new();
85    let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff;
86    for ts_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) {
87        let time_slices: Vec<_> = ts_selection.iter(time_slice_info).collect();
88        let demand_for_selection: Flow = time_slices.iter().map(|(ts, _)| demand[*ts]).sum();
89        let supply_for_selection: Flow = time_slices
90            .iter()
91            .map(|(ts, _)| activity[*ts] * flow_coeff)
92            .sum();
93
94        #[allow(clippy::cast_precision_loss)]
95        let unmet_per_slice = (demand_for_selection - supply_for_selection).max(Flow(0.0))
96            / Dimensionless(time_slices.len() as f64);
97        for (time_slice, _) in &time_slices {
98            unmet_demand.insert((*time_slice).clone(), unmet_per_slice);
99        }
100    }
101    unmet_demand
102}
103
104/// Performs optimisation for an asset, given the coefficients and demand.
105pub fn perform_optimisation(
106    model: &Model,
107    asset: &AssetRef,
108    commodity: &Commodity,
109    coefficients: &ObjectiveCoefficients,
110    demand: &DemandMap,
111) -> Result<ResultsMap> {
112    // Create problem and add variables
113    let mut problem = Problem::default();
114    let activity_vars = add_activity_vars(&mut problem, coefficients);
115
116    // Add constraints
117    add_constraints(
118        &mut problem,
119        asset,
120        commodity,
121        &activity_vars,
122        demand,
123        &model.time_slice_info,
124    );
125
126    // Solve model
127    let mut highs_model = problem.optimise(Sense::Maximise);
128    apply_highs_options_from_toml(&mut highs_model, &model.parameters.highs.appraisal_options)
129        .context("Failed to apply custom HiGHS options to appraisal optimisation")?;
130    let solution = solve_optimal(highs_model)
131        .map_err(ModelError::into_anyhow)?
132        .get_solution();
133    let solution_values = solution.columns();
134    let activity = activity_vars
135        .keys()
136        .zip(solution_values.iter())
137        .map(|(time_slice, &value)| (time_slice.clone(), Activity::new(value)))
138        .collect();
139    let unmet_demand =
140        compute_unmet_demand(demand, &activity, commodity, asset, &model.time_slice_info);
141    Ok(ResultsMap {
142        activity,
143        unmet_demand,
144    })
145}
muse2/simulation/investment/appraisal/optimisation.rs

muse2/simulation/investment/appraisal/
optimisation.rs
