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

muse2/simulation/
optimisation.rs
