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};
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(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(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, f64)> {
91        self.zip_var_keys_with_output(&self.active_asset_var_idx, self.solution.columns())
92    }
93
94    /// Reduced costs for candidate assets
95    pub fn iter_reduced_costs_for_candidates(
96        &self,
97    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, f64)> {
98        self.zip_var_keys_with_output(&self.candidate_asset_var_idx, self.solution.dual_columns())
99    }
100
101    /// Keys and dual values for commodity balance constraints.
102    pub fn iter_commodity_balance_duals(
103        &self,
104    ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, f64)> {
105        // Each commodity balance constraint applies to a particular time slice
106        // selection (depending on time slice level). Where this covers multiple timeslices,
107        // we return the same dual for each individual timeslice.
108        self.constraint_keys
109            .commodity_balance_keys
110            .zip_duals(self.solution.dual_rows())
111            .flat_map(|((commodity_id, region_id, ts_selection), price)| {
112                ts_selection
113                    .iter(self.time_slice_info)
114                    .map(move |(ts, _)| (commodity_id, region_id, ts, price))
115            })
116    }
117
118    /// Keys and dual values for activity constraints.
119    pub fn iter_activity_duals(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, f64)> {
120        self.constraint_keys
121            .activity_keys
122            .zip_duals(self.solution.dual_rows())
123            .map(|((asset, time_slice), dual)| (asset, time_slice, dual))
124    }
125
126    /// Zip a subset of keys in the variable map with a subset of the given output variable.
127    ///
128    /// # Arguments
129    ///
130    /// * `variable_idx` - The subset of variables to look at
131    /// * `output` - The output variable of interest
132    fn zip_var_keys_with_output<'a>(
133        &'a self,
134        variable_idx: &Range<usize>,
135        output: &'a [f64],
136    ) -> impl Iterator<Item = (&'a AssetRef, &'a TimeSliceID, f64)> {
137        assert!(variable_idx.end <= output.len());
138        self.variables
139            .0
140            .keys()
141            .zip(output[variable_idx.clone()].iter())
142            .map(|((asset, time_slice), value)| (asset, time_slice, *value))
143    }
144}
145
146/// Perform the dispatch optimisation.
147///
148/// For a detailed description, please see the [dispatch optimisation formulation][1].
149///
150/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html
151///
152/// # Arguments
153///
154/// * `model` - The model
155/// * `asset_pool` - The asset pool
156/// * `candidate_assets` - Candidate assets for inclusion in active pool
157/// * `year` - Current milestone year
158///
159/// # Returns
160///
161/// A solution containing new commodity flows for assets and prices for (some) commodities.
162pub fn perform_dispatch_optimisation<'a>(
163    model: &'a Model,
164    asset_pool: &AssetPool,
165    candidate_assets: &[AssetRef],
166    year: u32,
167) -> Result<Solution<'a>> {
168    // Set up problem
169    let mut problem = Problem::default();
170    let mut variables = VariableMap::default();
171    let active_asset_var_idx = add_variables(
172        &mut problem,
173        &mut variables,
174        &model.time_slice_info,
175        asset_pool.as_slice(),
176        year,
177    );
178    let candidate_asset_var_idx = add_variables(
179        &mut problem,
180        &mut variables,
181        &model.time_slice_info,
182        candidate_assets,
183        year,
184    );
185
186    // Add constraints
187    let all_assets = chain(asset_pool.iter(), candidate_assets.iter());
188    let constraint_keys = add_asset_constraints(&mut problem, &variables, model, all_assets, year);
189
190    // Solve problem
191    let mut highs_model = problem.optimise(Sense::Minimise);
192
193    // **HACK**: Dump output of HiGHS solver to stdout. Among other things, this includes the
194    // objective value for the solution. Sadly it doesn't go via our logger, so this information
195    // will not be included in the log file.
196    //
197    // Should be removed when we write the objective value to output data properly. See:
198    //   https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/428
199    enable_highs_logging(&mut highs_model);
200
201    // Solve model
202    let solution = highs_model.solve();
203    match solution.status() {
204        HighsModelStatus::Optimal => Ok(Solution {
205            solution: solution.get_solution(),
206            variables,
207            active_asset_var_idx,
208            candidate_asset_var_idx,
209            time_slice_info: &model.time_slice_info,
210            constraint_keys,
211        }),
212        status => Err(anyhow!("Could not solve: {status:?}")),
213    }
214}
215
216/// Enable logging for the HiGHS solver
217fn enable_highs_logging(model: &mut highs::Model) {
218    // **HACK**: Skip this step if logging is disabled (e.g. when running tests)
219    if let Ok(log_level) = std::env::var("MUSE2_LOG_LEVEL") {
220        if log_level.eq_ignore_ascii_case("off") {
221            return;
222        }
223    }
224
225    model.set_option("log_to_console", true);
226    model.set_option("output_flag", true);
227}
228
229/// Add variables to the optimisation problem.
230///
231/// # Arguments
232///
233/// * `problem` - The optimisation problem
234/// * `variables` - The variable map
235/// * `time_slice_info` - Information about assets
236/// * `assets` - Assets to include
237/// * `year` - Current milestone year
238fn add_variables(
239    problem: &mut Problem,
240    variables: &mut VariableMap,
241    time_slice_info: &TimeSliceInfo,
242    assets: &[AssetRef],
243    year: u32,
244) -> Range<usize> {
245    // This line **must** come before we add more variables
246    let start = problem.num_cols();
247
248    for (asset, time_slice) in iproduct!(assets.iter(), time_slice_info.iter_ids()) {
249        let coeff = calculate_cost_coefficient(asset, year, time_slice);
250        let var = problem.add_column(coeff.value(), 0.0..);
251        let key = (asset.clone(), time_slice.clone());
252        let existing = variables.0.insert(key, var).is_some();
253        assert!(!existing, "Duplicate entry for var");
254    }
255
256    start..problem.num_cols()
257}
258
259/// Calculate the cost coefficient for a decision variable
260fn calculate_cost_coefficient(
261    asset: &Asset,
262    year: u32,
263    time_slice: &TimeSliceID,
264) -> MoneyPerActivity {
265    // The cost for all commodity flows (including levies/incentives)
266    let flows_cost: MoneyPerActivity = asset
267        .iter_flows()
268        .map(|flow| flow.get_total_cost(&asset.region_id, year, time_slice))
269        .sum();
270
271    asset.process_parameter.variable_operating_cost + flows_cost
272}
muse2/simulation/optimisation.rs

muse2/simulation/
optimisation.rs
