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, AssetRef};
5use crate::commodity::CommodityID;
6use crate::model::Model;
7use crate::output::DataWriter;
8use crate::region::RegionID;
9use crate::time_slice::{TimeSliceID, TimeSliceInfo};
10use crate::units::{Activity, Flow, Money, MoneyPerActivity, MoneyPerFlow, UnitType};
11use anyhow::{Result, anyhow, ensure};
12use highs::{HighsModelStatus, RowProblem as Problem, Sense};
13use indexmap::IndexMap;
14use itertools::{chain, iproduct};
15use log::debug;
16use std::collections::{HashMap, HashSet};
17use std::ops::Range;
18
19mod constraints;
20use constraints::{ConstraintKeys, add_asset_constraints};
21
22/// A map of commodity flows calculated during the optimisation
23pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), Flow>;
24
25/// A decision variable in the optimisation
26///
27/// Note that this type does **not** include the value of the variable; it just refers to a
28/// particular column of the problem.
29type Variable = highs::Col;
30
31/// A map for easy lookup of variables in the problem.
32///
33/// The entries are ordered (see [`IndexMap`]).
34///
35/// We use this data structure for two things:
36///
37/// 1. In order define constraints for the optimisation
38/// 2. To keep track of the combination of parameters that each variable corresponds to, for when we
39///    are reading the results of the optimisation.
40#[derive(Default)]
41pub struct VariableMap(IndexMap<(AssetRef, TimeSliceID), Variable>);
42
43impl VariableMap {
44    /// Get the [`Variable`] corresponding to the given parameters.
45    fn get(&self, asset: &AssetRef, time_slice: &TimeSliceID) -> Variable {
46        let key = (asset.clone(), time_slice.clone());
47
48        *self
49            .0
50            .get(&key)
51            .expect("No variable found for given params")
52    }
53
54    /// Iterate over the variable map
55    fn iter(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Variable)> {
56        self.0
57            .iter()
58            .map(|((asset, time_slice), var)| (asset, time_slice, *var))
59    }
60}
61
62/// The solution to the dispatch optimisation problem
63pub struct Solution<'a> {
64    solution: highs::Solution,
65    variables: VariableMap,
66    active_asset_var_idx: Range<usize>,
67    candidate_asset_var_idx: Range<usize>,
68    time_slice_info: &'a TimeSliceInfo,
69    constraint_keys: ConstraintKeys,
70    /// The objective value for the solution
71    pub objective_value: Money,
72}
73
74impl Solution<'_> {
75    /// Create a map of commodity flows for each asset's coeffs at every time slice.
76    ///
77    /// Note that this only includes commodity flows which relate to assets, so not every commodity
78    /// in the simulation will necessarily be represented.
79    pub fn create_flow_map(&self) -> FlowMap {
80        // The decision variables represent assets' activity levels, not commodity flows. We
81        // multiply this value by the flow coeffs to get commodity flows.
82        let mut flows = FlowMap::new();
83        for (asset, time_slice, activity) in self.iter_activity_for_active() {
84            for flow in asset.iter_flows() {
85                let flow_key = (asset.clone(), flow.commodity.id.clone(), time_slice.clone());
86                let flow_value = activity * flow.coeff;
87                flows.insert(flow_key, flow_value);
88            }
89        }
90
91        flows
92    }
93
94    /// Activity for each active asset
95    pub fn iter_activity(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
96        self.variables
97            .0
98            .keys()
99            .zip(self.solution.columns())
100            .map(|((asset, time_slice), activity)| (asset, time_slice, Activity(*activity)))
101    }
102
103    /// Activity for each active asset
104    fn iter_activity_for_active(
105        &self,
106    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
107        self.zip_var_keys_with_output(&self.active_asset_var_idx, self.solution.columns())
108    }
109
110    /// Reduced costs for candidate assets
111    pub fn iter_reduced_costs_for_candidates(
112        &self,
113    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
114        self.zip_var_keys_with_output(&self.candidate_asset_var_idx, self.solution.dual_columns())
115    }
116
117    /// Keys and dual values for commodity balance constraints.
118    pub fn iter_commodity_balance_duals(
119        &self,
120    ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, MoneyPerFlow)> {
121        // Each commodity balance constraint applies to a particular time slice
122        // selection (depending on time slice level). Where this covers multiple timeslices,
123        // we return the same dual for each individual timeslice.
124        self.constraint_keys
125            .commodity_balance_keys
126            .zip_duals(self.solution.dual_rows())
127            .flat_map(|((commodity_id, region_id, ts_selection), price)| {
128                ts_selection
129                    .iter(self.time_slice_info)
130                    .map(move |(ts, _)| (commodity_id, region_id, ts, price))
131            })
132    }
133
134    /// Keys and dual values for activity constraints.
135    pub fn iter_activity_duals(
136        &self,
137    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
138        self.constraint_keys
139            .activity_keys
140            .zip_duals(self.solution.dual_rows())
141            .map(|((asset, time_slice), dual)| (asset, time_slice, dual))
142    }
143
144    /// Zip a subset of keys in the variable map with a subset of the given output variable.
145    ///
146    /// # Arguments
147    ///
148    /// * `variable_idx` - The subset of variables to look at
149    /// * `output` - The output variable of interest
150    fn zip_var_keys_with_output<'a, T: UnitType>(
151        &'a self,
152        variable_idx: &Range<usize>,
153        output: &'a [f64],
154    ) -> impl Iterator<Item = (&'a AssetRef, &'a TimeSliceID, T)> + use<'a, T> {
155        let keys = self.variables.0.keys().skip(variable_idx.start);
156        assert!(keys.len() >= variable_idx.len());
157
158        keys.zip(output[variable_idx.clone()].iter())
159            .map(|((asset, time_slice), value)| (asset, time_slice, T::new(*value)))
160    }
161}
162
163/// Try to solve the model, returning an error if the model is incoherent or result is non-optimal
164pub fn solve_optimal(model: highs::Model) -> Result<highs::SolvedModel> {
165    let solved = model
166        .try_solve()
167        .map_err(|err| anyhow!("Incoherent model: {err:?}"))?;
168
169    let status = solved.status();
170    ensure!(
171        status == HighsModelStatus::Optimal,
172        "Could not find optimal result for model: {status:?}"
173    );
174
175    Ok(solved)
176}
177
178/// Sanity check for input prices.
179///
180/// Input prices should only be provided for commodities for which there will be no commodity
181/// balance constraint.
182fn check_input_prices(
183    input_prices: &HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>,
184    commodities: &[CommodityID],
185) {
186    let commodities_set: HashSet<_> = commodities.iter().collect();
187    let has_prices_for_commodity_subset = input_prices
188        .keys()
189        .any(|(commodity_id, _, _)| commodities_set.contains(commodity_id));
190    assert!(
191        !has_prices_for_commodity_subset,
192        "Input prices were included for commodities that are being modelled, which is not allowed."
193    )
194}
195
196/// Provides the interface for running the dispatch optimisation.
197///
198/// For a detailed description, please see the [dispatch optimisation formulation][1].
199///
200/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/model/dispatch_optimisation.html
201pub struct DispatchRun<'model, 'run> {
202    model: &'model Model,
203    existing_assets: &'run [AssetRef],
204    candidate_assets: &'run [AssetRef],
205    commodities: &'run [CommodityID],
206    input_prices: Option<&'run HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>>,
207    year: u32,
208}
209
210impl<'model, 'run> DispatchRun<'model, 'run> {
211    /// Create a new [`DispatchRun`] for the specified model and assets for a given year
212    pub fn new(model: &'model Model, assets: &'run [AssetRef], year: u32) -> Self {
213        Self {
214            model,
215            existing_assets: assets,
216            candidate_assets: &[],
217            commodities: &[],
218            input_prices: None,
219            year,
220        }
221    }
222
223    /// Include the specified candidate assets in the dispatch run
224    pub fn with_candidates(self, candidate_assets: &'run [AssetRef]) -> Self {
225        Self {
226            candidate_assets,
227            ..self
228        }
229    }
230
231    /// Only apply commodity balance constraints to the specified subset of commodities
232    pub fn with_commodity_subset(self, commodities: &'run [CommodityID]) -> Self {
233        assert!(!commodities.is_empty());
234
235        Self {
236            commodities,
237            ..self
238        }
239    }
240
241    /// Explicitly provide prices for certain input commodities
242    pub fn with_input_prices(
243        self,
244        input_prices: &'run HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>,
245    ) -> Self {
246        Self {
247            input_prices: Some(input_prices),
248            ..self
249        }
250    }
251
252    /// Perform the dispatch optimisation.
253    ///
254    /// # Arguments
255    ///
256    /// * `run_description` - Which dispatch run for the current year this is
257    /// * `writer` - For saving output data
258    ///
259    /// # Returns
260    ///
261    /// A solution containing new commodity flows for assets and prices for (some) commodities or an
262    /// error.
263    pub fn run(self, run_description: &str, writer: &mut DataWriter) -> Result<Solution<'model>> {
264        let solution = self.run_no_save()?;
265        writer.write_dispatch_debug_info(self.year, run_description, &solution)?;
266        Ok(solution)
267    }
268
269    /// Run dispatch without saving the results.
270    ///
271    /// This is an internal function as callers always want to save results.
272    fn run_no_save(&self) -> Result<Solution<'model>> {
273        // Set up problem
274        let mut problem = Problem::default();
275        let mut variables = VariableMap::default();
276        let active_asset_var_idx = add_variables(
277            &mut problem,
278            &mut variables,
279            &self.model.time_slice_info,
280            self.input_prices,
281            self.existing_assets,
282            self.year,
283        );
284        let candidate_asset_var_idx = add_variables(
285            &mut problem,
286            &mut variables,
287            &self.model.time_slice_info,
288            self.input_prices,
289            self.candidate_assets,
290            self.year,
291        );
292
293        // If the user provided no commodities, we all use of them
294        let all_commodities: Vec<_>;
295        let commodities = if !self.commodities.is_empty() {
296            self.commodities
297        } else {
298            all_commodities = self.model.commodities.keys().cloned().collect();
299            &all_commodities
300        };
301        if let Some(input_prices) = self.input_prices {
302            check_input_prices(input_prices, commodities);
303        }
304
305        // Add constraints
306        let all_assets = chain(self.existing_assets.iter(), self.candidate_assets.iter());
307        let constraint_keys = add_asset_constraints(
308            &mut problem,
309            &variables,
310            self.model,
311            all_assets,
312            commodities,
313            self.year,
314        );
315
316        // Solve model
317        let solution = solve_optimal(problem.optimise(Sense::Minimise))?;
318
319        let objective_value = Money(solution.objective_value());
320        debug!("Objective value: {objective_value}");
321
322        Ok(Solution {
323            solution: solution.get_solution(),
324            variables,
325            active_asset_var_idx,
326            candidate_asset_var_idx,
327            time_slice_info: &self.model.time_slice_info,
328            constraint_keys,
329            objective_value,
330        })
331    }
332}
333
334/// Add variables to the optimisation problem.
335///
336/// # Arguments
337///
338/// * `problem` - The optimisation problem
339/// * `variables` - The variable map
340/// * `time_slice_info` - Information about assets
341/// * `input_prices` - Optional explicit prices for input commodities
342/// * `assets` - Assets to include
343/// * `year` - Current milestone year
344fn add_variables(
345    problem: &mut Problem,
346    variables: &mut VariableMap,
347    time_slice_info: &TimeSliceInfo,
348    input_prices: Option<&HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>>,
349    assets: &[AssetRef],
350    year: u32,
351) -> Range<usize> {
352    // This line **must** come before we add more variables
353    let start = problem.num_cols();
354
355    for (asset, time_slice) in iproduct!(assets.iter(), time_slice_info.iter_ids()) {
356        let coeff = calculate_cost_coefficient(asset, year, time_slice, input_prices);
357        let var = problem.add_column(coeff.value(), 0.0..);
358        let key = (asset.clone(), time_slice.clone());
359        let existing = variables.0.insert(key, var).is_some();
360        assert!(!existing, "Duplicate entry for var");
361    }
362
363    start..problem.num_cols()
364}
365
366/// Calculate the cost coefficient for a decision variable.
367///
368/// Normally, the cost coefficient is the same as the asset's operating costs for the given year and
369/// time slice. If `input_prices` is provided then those prices are added to the flow costs for the
370/// relevant commodities, if they are input flows for the asset.
371///
372/// # Arguments
373///
374/// * `asset` - The asset to calculate the coefficient for
375/// * `year` - The current milestone year
376/// * `time_slice` - The time slice to which this coefficient applies
377/// * `input_prices` - Optional map of prices to include for input commodities
378///
379/// # Returns
380///
381/// The cost coefficient to be used for the relevant decision variable.
382fn calculate_cost_coefficient(
383    asset: &Asset,
384    year: u32,
385    time_slice: &TimeSliceID,
386    input_prices: Option<&HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>>,
387) -> MoneyPerActivity {
388    let opex = asset.get_operating_cost(year, time_slice);
389    let input_cost = input_prices
390        .map(|prices| asset.get_input_cost_from_prices(prices, time_slice))
391        .unwrap_or_default();
392    opex + input_cost
393}
muse2/simulation/optimisation.rs

muse2/simulation/
optimisation.rs
