muse2/simulation/
optimisation.rs

1//! Code for performing dispatch optimisation.
2//!
3//! This is used to calculate commodity flows and prices.
4use crate::asset::{Asset, AssetPool, AssetRef};
5use crate::commodity::CommodityID;
6use crate::model::Model;
7use crate::region::RegionID;
8use crate::time_slice::{TimeSliceID, TimeSliceInfo};
9use crate::units::{Activity, Flow, MoneyPerActivity, MoneyPerFlow, UnitType};
10use anyhow::{anyhow, Result};
11use highs::{HighsModelStatus, RowProblem as Problem, Sense};
12use indexmap::IndexMap;
13use itertools::{chain, iproduct};
14use std::ops::Range;
15
16mod constraints;
17use constraints::{add_asset_constraints, ConstraintKeys};
18
19/// A map of commodity flows calculated during the optimisation
20pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), Flow>;
21
22/// A decision variable in the optimisation
23///
24/// Note that this type does **not** include the value of the variable; it just refers to a
25/// particular column of the problem.
26type Variable = highs::Col;
27
28/// A map for easy lookup of variables in the problem.
29///
30/// The entries are ordered (see [`IndexMap`]).
31///
32/// We use this data structure for two things:
33///
34/// 1. In order define constraints for the optimisation
35/// 2. To keep track of the combination of parameters that each variable corresponds to, for when we
36///    are reading the results of the optimisation.
37#[derive(Default)]
38pub struct VariableMap(IndexMap<(AssetRef, TimeSliceID), Variable>);
39
40impl VariableMap {
41    /// Get the [`Variable`] corresponding to the given parameters.
42    fn get(&self, asset: &AssetRef, time_slice: &TimeSliceID) -> Variable {
43        let key = (asset.clone(), time_slice.clone());
44
45        *self
46            .0
47            .get(&key)
48            .expect("No variable found for given params")
49    }
50
51    /// Iterate over the variable map
52    fn iter(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Variable)> {
53        self.0
54            .iter()
55            .map(|((asset, time_slice), var)| (asset, time_slice, *var))
56    }
57}
58
59/// The solution to the dispatch optimisation problem
60pub struct Solution<'a> {
61    solution: highs::Solution,
62    variables: VariableMap,
63    active_asset_var_idx: Range<usize>,
64    candidate_asset_var_idx: Range<usize>,
65    time_slice_info: &'a TimeSliceInfo,
66    constraint_keys: ConstraintKeys,
67}
68
69impl Solution<'_> {
70    /// Create a map of commodity flows for each asset's coeffs at every time slice.
71    ///
72    /// Note that this only includes commodity flows which relate to assets, so not every commodity
73    /// in the simulation will necessarily be represented.
74    pub fn create_flow_map(&self) -> FlowMap {
75        // The decision variables represent assets' activity levels, not commodity flows. We
76        // multiply this value by the flow coeffs to get commodity flows.
77        let mut flows = FlowMap::new();
78        for (asset, time_slice, activity) in self.iter_activity_for_active() {
79            for flow in asset.iter_flows() {
80                let flow_key = (asset.clone(), flow.commodity.id.clone(), time_slice.clone());
81                let flow_value = activity * flow.coeff;
82                flows.insert(flow_key, flow_value);
83            }
84        }
85
86        flows
87    }
88
89    /// Activity for each active asset
90    fn iter_activity_for_active(
91        &self,
92    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
93        self.zip_var_keys_with_output(&self.active_asset_var_idx, self.solution.columns())
94    }
95
96    /// Reduced costs for candidate assets
97    pub fn iter_reduced_costs_for_candidates(
98        &self,
99    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
100        self.zip_var_keys_with_output(&self.candidate_asset_var_idx, self.solution.dual_columns())
101    }
102
103    /// Keys and dual values for commodity balance constraints.
104    pub fn iter_commodity_balance_duals(
105        &self,
106    ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, MoneyPerFlow)> {
107        // Each commodity balance constraint applies to a particular time slice
108        // selection (depending on time slice level). Where this covers multiple timeslices,
109        // we return the same dual for each individual timeslice.
110        self.constraint_keys
111            .commodity_balance_keys
112            .zip_duals(self.solution.dual_rows())
113            .flat_map(|((commodity_id, region_id, ts_selection), price)| {
114                ts_selection
115                    .iter(self.time_slice_info)
116                    .map(move |(ts, _)| (commodity_id, region_id, ts, price))
117            })
118    }
119
120    /// Keys and dual values for activity constraints.
121    pub fn iter_activity_duals(
122        &self,
123    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
124        self.constraint_keys
125            .activity_keys
126            .zip_duals(self.solution.dual_rows())
127            .map(|((asset, time_slice), dual)| (asset, time_slice, dual))
128    }
129
130    /// Zip a subset of keys in the variable map with a subset of the given output variable.
131    ///
132    /// # Arguments
133    ///
134    /// * `variable_idx` - The subset of variables to look at
135    /// * `output` - The output variable of interest
136    fn zip_var_keys_with_output<'a, T: UnitType>(
137        &'a self,
138        variable_idx: &Range<usize>,
139        output: &'a [f64],
140    ) -> impl Iterator<Item = (&'a AssetRef, &'a TimeSliceID, T)> {
141        assert!(variable_idx.end <= output.len());
142        self.variables
143            .0
144            .keys()
145            .zip(output[variable_idx.clone()].iter())
146            .map(|((asset, time_slice), value)| (asset, time_slice, T::new(*value)))
147    }
148}
149
150/// Perform the dispatch optimisation.
151///
152/// For a detailed description, please see the [dispatch optimisation formulation][1].
153///
154/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html
155///
156/// # Arguments
157///
158/// * `model` - The model
159/// * `asset_pool` - The asset pool
160/// * `candidate_assets` - Candidate assets for inclusion in active pool
161/// * `year` - Current milestone year
162///
163/// # Returns
164///
165/// A solution containing new commodity flows for assets and prices for (some) commodities.
166pub fn perform_dispatch_optimisation<'a>(
167    model: &'a Model,
168    asset_pool: &AssetPool,
169    candidate_assets: &[AssetRef],
170    year: u32,
171) -> Result<Solution<'a>> {
172    // Set up problem
173    let mut problem = Problem::default();
174    let mut variables = VariableMap::default();
175    let active_asset_var_idx = add_variables(
176        &mut problem,
177        &mut variables,
178        &model.time_slice_info,
179        asset_pool.as_slice(),
180        year,
181    );
182    let candidate_asset_var_idx = add_variables(
183        &mut problem,
184        &mut variables,
185        &model.time_slice_info,
186        candidate_assets,
187        year,
188    );
189
190    // Add constraints
191    let all_assets = chain(asset_pool.iter(), candidate_assets.iter());
192    let constraint_keys = add_asset_constraints(&mut problem, &variables, model, all_assets, year);
193
194    // Solve problem
195    let mut highs_model = problem.optimise(Sense::Minimise);
196
197    // **HACK**: Dump output of HiGHS solver to stdout. Among other things, this includes the
198    // objective value for the solution. Sadly it doesn't go via our logger, so this information
199    // will not be included in the log file.
200    //
201    // Should be removed when we write the objective value to output data properly. See:
202    //   https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/428
203    enable_highs_logging(&mut highs_model);
204
205    // Solve model
206    let solution = highs_model.solve();
207    match solution.status() {
208        HighsModelStatus::Optimal => Ok(Solution {
209            solution: solution.get_solution(),
210            variables,
211            active_asset_var_idx,
212            candidate_asset_var_idx,
213            time_slice_info: &model.time_slice_info,
214            constraint_keys,
215        }),
216        status => Err(anyhow!("Could not solve: {status:?}")),
217    }
218}
219
220/// Enable logging for the HiGHS solver
221fn enable_highs_logging(model: &mut highs::Model) {
222    // **HACK**: Skip this step if logging is disabled (e.g. when running tests)
223    if let Ok(log_level) = std::env::var("MUSE2_LOG_LEVEL") {
224        if log_level.eq_ignore_ascii_case("off") {
225            return;
226        }
227    }
228
229    model.set_option("log_to_console", true);
230    model.set_option("output_flag", true);
231}
232
233/// Add variables to the optimisation problem.
234///
235/// # Arguments
236///
237/// * `problem` - The optimisation problem
238/// * `variables` - The variable map
239/// * `time_slice_info` - Information about assets
240/// * `assets` - Assets to include
241/// * `year` - Current milestone year
242fn add_variables(
243    problem: &mut Problem,
244    variables: &mut VariableMap,
245    time_slice_info: &TimeSliceInfo,
246    assets: &[AssetRef],
247    year: u32,
248) -> Range<usize> {
249    // This line **must** come before we add more variables
250    let start = problem.num_cols();
251
252    for (asset, time_slice) in iproduct!(assets.iter(), time_slice_info.iter_ids()) {
253        let coeff = calculate_cost_coefficient(asset, year, time_slice);
254        let var = problem.add_column(coeff.value(), 0.0..);
255        let key = (asset.clone(), time_slice.clone());
256        let existing = variables.0.insert(key, var).is_some();
257        assert!(!existing, "Duplicate entry for var");
258    }
259
260    start..problem.num_cols()
261}
262
263/// Calculate the cost coefficient for a decision variable
264fn calculate_cost_coefficient(
265    asset: &Asset,
266    year: u32,
267    time_slice: &TimeSliceID,
268) -> MoneyPerActivity {
269    // The cost for all commodity flows (including levies/incentives)
270    let flows_cost: MoneyPerActivity = asset
271        .iter_flows()
272        .map(|flow| flow.get_total_cost(&asset.region_id, year, time_slice))
273        .sum();
274
275    asset.process_parameter.variable_operating_cost + flows_cost
276}
muse2/simulation/optimisation.rs

muse2/simulation/
optimisation.rs
