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::input::format_items_with_cap;
7use crate::model::Model;
8use crate::output::DataWriter;
9use crate::region::RegionID;
10use crate::simulation::CommodityPrices;
11use crate::time_slice::{TimeSliceID, TimeSliceInfo};
12use crate::units::{Activity, Flow, Money, MoneyPerActivity, MoneyPerFlow, UnitType};
13use anyhow::{Result, bail, ensure};
14use highs::{HighsModelStatus, HighsStatus, RowProblem as Problem, Sense};
15use indexmap::{IndexMap, IndexSet};
16use itertools::{chain, iproduct};
17use log::debug;
18use std::collections::HashSet;
19use std::error::Error;
20use std::fmt;
21use std::ops::Range;
22
23mod constraints;
24use constraints::{ConstraintKeys, add_asset_constraints};
25
26/// A map of commodity flows calculated during the optimisation
27pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), Flow>;
28
29/// A decision variable in the optimisation
30///
31/// Note that this type does **not** include the value of the variable; it just refers to a
32/// particular column of the problem.
33type Variable = highs::Col;
34
35/// The map of variables related to assets
36type AssetVariableMap = IndexMap<(AssetRef, TimeSliceID), Variable>;
37
38/// Variables representing unmet demand for a given commodity
39type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Variable>;
40
41/// A map for easy lookup of variables in the problem.
42///
43/// The entries are ordered (see [`IndexMap`]).
44///
45/// We use this data structure for two things:
46///
47/// 1. In order define constraints for the optimisation
48/// 2. To keep track of the combination of parameters that each variable corresponds to, for when we
49///    are reading the results of the optimisation.
50pub struct VariableMap {
51    asset_vars: AssetVariableMap,
52    existing_asset_var_idx: Range<usize>,
53    unmet_demand_vars: UnmetDemandVariableMap,
54    unmet_demand_var_idx: Range<usize>,
55}
56
57impl VariableMap {
58    /// Create a new [`VariableMap`] and add variables to the problem
59    ///
60    /// # Arguments
61    ///
62    /// * `problem` - The optimisation problem
63    /// * `model` - The model
64    /// * `input_prices` - Optional explicit prices for input commodities
65    /// * `existing_assets` - The asset pool
66    /// * `candidate_assets` - Candidate assets for inclusion in active pool
67    /// * `year` - Current milestone year
68    fn new_with_asset_vars(
69        problem: &mut Problem,
70        model: &Model,
71        input_prices: Option<&CommodityPrices>,
72        existing_assets: &[AssetRef],
73        candidate_assets: &[AssetRef],
74        year: u32,
75    ) -> Self {
76        let mut asset_vars = AssetVariableMap::new();
77        let existing_asset_var_idx = add_asset_variables(
78            problem,
79            &mut asset_vars,
80            &model.time_slice_info,
81            input_prices,
82            existing_assets,
83            year,
84        );
85        add_asset_variables(
86            problem,
87            &mut asset_vars,
88            &model.time_slice_info,
89            input_prices,
90            candidate_assets,
91            year,
92        );
93
94        Self {
95            asset_vars,
96            existing_asset_var_idx,
97            unmet_demand_vars: UnmetDemandVariableMap::default(),
98            unmet_demand_var_idx: Range::default(),
99        }
100    }
101
102    /// Add unmet demand variables to the map and the problem
103    ///
104    /// # Arguments
105    ///
106    /// * `problem` - The optimisation problem
107    /// * `model` - The model
108    /// * `commodities` - The subset of commodities the problem is being run for
109    fn add_unmet_demand_variables(
110        &mut self,
111        problem: &mut Problem,
112        model: &Model,
113        commodities: &[CommodityID],
114    ) {
115        assert!(!commodities.is_empty());
116
117        // This line **must** come before we add more variables
118        let start = problem.num_cols();
119
120        // Add variables
121        let voll = model.parameters.value_of_lost_load;
122        self.unmet_demand_vars.extend(
123            iproduct!(
124                commodities.iter(),
125                model.iter_regions(),
126                model.time_slice_info.iter_ids()
127            )
128            .map(|(commodity_id, region_id, time_slice)| {
129                let key = (commodity_id.clone(), region_id.clone(), time_slice.clone());
130                let var = problem.add_column(voll.value(), 0.0..);
131
132                (key, var)
133            }),
134        );
135
136        self.unmet_demand_var_idx = start..problem.num_cols();
137    }
138
139    /// Get the asset [`Variable`] corresponding to the given parameters.
140    fn get_asset_var(&self, asset: &AssetRef, time_slice: &TimeSliceID) -> Variable {
141        let key = (asset.clone(), time_slice.clone());
142
143        *self
144            .asset_vars
145            .get(&key)
146            .expect("No asset variable found for given params")
147    }
148
149    /// Get the unmet demand [`Variable`] corresponding to the given parameters.
150    fn get_unmet_demand_var(
151        &self,
152        commodity_id: &CommodityID,
153        region_id: &RegionID,
154        time_slice: &TimeSliceID,
155    ) -> Variable {
156        let key = (commodity_id.clone(), region_id.clone(), time_slice.clone());
157
158        *self
159            .unmet_demand_vars
160            .get(&key)
161            .expect("No unmet demand variable for given params")
162    }
163
164    /// Iterate over the asset variables
165    fn iter_asset_vars(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Variable)> {
166        self.asset_vars
167            .iter()
168            .map(|((asset, time_slice), var)| (asset, time_slice, *var))
169    }
170
171    /// Iterate over the keys for asset variables
172    fn asset_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> {
173        self.asset_vars.keys()
174    }
175}
176
177/// The solution to the dispatch optimisation problem
178#[allow(clippy::struct_field_names)]
179pub struct Solution<'a> {
180    solution: highs::Solution,
181    variables: VariableMap,
182    time_slice_info: &'a TimeSliceInfo,
183    constraint_keys: ConstraintKeys,
184    /// The objective value for the solution
185    pub objective_value: Money,
186}
187
188impl Solution<'_> {
189    /// Create a map of commodity flows for each asset's coeffs at every time slice.
190    ///
191    /// Note that this only includes commodity flows which relate to assets, so not every commodity
192    /// in the simulation will necessarily be represented.
193    pub fn create_flow_map(&self) -> FlowMap {
194        // The decision variables represent assets' activity levels, not commodity flows. We
195        // multiply this value by the flow coeffs to get commodity flows.
196        let mut flows = FlowMap::new();
197        for (asset, time_slice, activity) in self.iter_activity_for_existing() {
198            for flow in asset.iter_flows() {
199                let flow_key = (asset.clone(), flow.commodity.id.clone(), time_slice.clone());
200                let flow_value = activity * flow.coeff;
201                flows.insert(flow_key, flow_value);
202            }
203        }
204
205        flows
206    }
207
208    /// Activity for each existing asset
209    pub fn iter_activity(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
210        self.variables
211            .asset_var_keys()
212            .zip(self.solution.columns())
213            .map(|((asset, time_slice), activity)| (asset, time_slice, Activity(*activity)))
214    }
215
216    /// Activity for each existing asset
217    fn iter_activity_for_existing(
218        &self,
219    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
220        self.zip_var_keys_with_output(
221            &self.variables.existing_asset_var_idx,
222            self.solution.columns(),
223        )
224    }
225
226    /// Iterate over unmet demand
227    pub fn iter_unmet_demand(
228        &self,
229    ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, Flow)> {
230        self.variables
231            .unmet_demand_vars
232            .keys()
233            .zip(self.solution.columns()[self.variables.unmet_demand_var_idx.clone()].iter())
234            .map(|((commodity_id, region_id, time_slice), flow)| {
235                (commodity_id, region_id, time_slice, Flow(*flow))
236            })
237    }
238
239    /// Keys and dual values for commodity balance constraints.
240    pub fn iter_commodity_balance_duals(
241        &self,
242    ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, MoneyPerFlow)> {
243        // Each commodity balance constraint applies to a particular time slice
244        // selection (depending on time slice level). Where this covers multiple time slices,
245        // we return the same dual for each individual time slice.
246        self.constraint_keys
247            .commodity_balance_keys
248            .zip_duals(self.solution.dual_rows())
249            .flat_map(|((commodity_id, region_id, ts_selection), price)| {
250                ts_selection
251                    .iter(self.time_slice_info)
252                    .map(move |(ts, _)| (commodity_id, region_id, ts, price))
253            })
254    }
255
256    /// Keys and dual values for activity constraints.
257    pub fn iter_activity_duals(
258        &self,
259    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
260        self.constraint_keys
261            .activity_keys
262            .zip_duals(self.solution.dual_rows())
263            .map(|((asset, time_slice), dual)| (asset, time_slice, dual))
264    }
265
266    /// Keys and values for column duals.
267    pub fn iter_column_duals(
268        &self,
269    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
270        self.variables
271            .asset_var_keys()
272            .zip(self.solution.dual_columns())
273            .map(|((asset, time_slice), dual)| (asset, time_slice, MoneyPerActivity(*dual)))
274    }
275
276    /// Zip a subset of keys in the variable map with a subset of the given output variable.
277    ///
278    /// # Arguments
279    ///
280    /// * `variable_idx` - The subset of variables to look at
281    /// * `output` - The output variable of interest
282    fn zip_var_keys_with_output<'a, T: UnitType>(
283        &'a self,
284        variable_idx: &Range<usize>,
285        output: &'a [f64],
286    ) -> impl Iterator<Item = (&'a AssetRef, &'a TimeSliceID, T)> + use<'a, T> {
287        let keys = self.variables.asset_var_keys().skip(variable_idx.start);
288        assert!(keys.len() >= variable_idx.len());
289
290        keys.zip(output[variable_idx.clone()].iter())
291            .map(|((asset, time_slice), value)| (asset, time_slice, T::new(*value)))
292    }
293}
294
295/// Defines the possible errors that can occur when running the solver
296#[derive(Debug, Clone)]
297pub enum ModelError {
298    /// The model definition is incoherent.
299    ///
300    /// Users should not be able to trigger this error.
301    Incoherent(HighsStatus),
302    /// An optimal solution could not be found
303    NonOptimal(HighsModelStatus),
304}
305
306impl fmt::Display for ModelError {
307    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
308        match self {
309            ModelError::Incoherent(status) => write!(f, "Incoherent model: {status:?}"),
310            ModelError::NonOptimal(status) => {
311                write!(f, "Could not find optimal result: {status:?}")
312            }
313        }
314    }
315}
316
317impl Error for ModelError {}
318
319/// Try to solve the model, returning an error if the model is incoherent or result is non-optimal
320pub fn solve_optimal(model: highs::Model) -> Result<highs::SolvedModel, ModelError> {
321    let solved = model.try_solve().map_err(ModelError::Incoherent)?;
322
323    match solved.status() {
324        HighsModelStatus::Optimal => Ok(solved),
325        status => Err(ModelError::NonOptimal(status)),
326    }
327}
328
329/// Sanity check for input prices.
330///
331/// Input prices should only be provided for commodities for which there will be no commodity
332/// balance constraint.
333fn check_input_prices(input_prices: &CommodityPrices, commodities: &[CommodityID]) {
334    let commodities_set: HashSet<_> = commodities.iter().collect();
335    let has_prices_for_commodity_subset = input_prices
336        .keys()
337        .any(|(commodity_id, _, _)| commodities_set.contains(commodity_id));
338    assert!(
339        !has_prices_for_commodity_subset,
340        "Input prices were included for commodities that are being modelled, which is not allowed."
341    );
342}
343
344/// Provides the interface for running the dispatch optimisation.
345///
346/// The caller can allow the dispatch run to return without error when demand is not met by calling
347/// the `with_unmet_demand_allowed` method.
348///
349/// For a detailed description, please see the [dispatch optimisation formulation][1].
350///
351/// [1]: https://energysystemsmodellinglab.github.io/MUSE2/model/dispatch_optimisation.html
352pub struct DispatchRun<'model, 'run> {
353    model: &'model Model,
354    existing_assets: &'run [AssetRef],
355    candidate_assets: &'run [AssetRef],
356    commodities: &'run [CommodityID],
357    input_prices: Option<&'run CommodityPrices>,
358    year: u32,
359}
360
361impl<'model, 'run> DispatchRun<'model, 'run> {
362    /// Create a new [`DispatchRun`] for the specified model and assets for a given year
363    pub fn new(model: &'model Model, assets: &'run [AssetRef], year: u32) -> Self {
364        Self {
365            model,
366            existing_assets: assets,
367            candidate_assets: &[],
368            commodities: &[],
369            input_prices: None,
370            year,
371        }
372    }
373
374    /// Include the specified candidate assets in the dispatch run
375    pub fn with_candidates(self, candidate_assets: &'run [AssetRef]) -> Self {
376        Self {
377            candidate_assets,
378            ..self
379        }
380    }
381
382    /// Only apply commodity balance constraints to the specified subset of commodities
383    pub fn with_commodity_subset(self, commodities: &'run [CommodityID]) -> Self {
384        assert!(!commodities.is_empty());
385
386        Self {
387            commodities,
388            ..self
389        }
390    }
391
392    /// Explicitly provide prices for certain input commodities
393    pub fn with_input_prices(self, input_prices: &'run CommodityPrices) -> Self {
394        Self {
395            input_prices: Some(input_prices),
396            ..self
397        }
398    }
399
400    /// Perform the dispatch optimisation.
401    ///
402    /// # Arguments
403    ///
404    /// * `run_description` - Which dispatch run for the current year this is
405    /// * `writer` - For saving output data
406    ///
407    /// # Returns
408    ///
409    /// A solution containing new commodity flows for assets and prices for (some) commodities or an
410    /// error.
411    pub fn run(self, run_description: &str, writer: &mut DataWriter) -> Result<Solution<'model>> {
412        let solution = self.run_no_save()?;
413        writer.write_dispatch_debug_info(self.year, run_description, &solution)?;
414        Ok(solution)
415    }
416
417    /// Run dispatch without saving the results.
418    ///
419    /// This is an internal function as callers always want to save results.
420    fn run_no_save(&self) -> Result<Solution<'model>> {
421        // If the user provided no commodities, we all use of them
422        let all_commodities: Vec<_>;
423        let commodities = if self.commodities.is_empty() {
424            all_commodities = self.model.commodities.keys().cloned().collect();
425            &all_commodities
426        } else {
427            self.commodities
428        };
429        if let Some(input_prices) = self.input_prices {
430            check_input_prices(input_prices, commodities);
431        }
432
433        // Try running dispatch. If it fails because the model is infeasible, it is likely that this
434        // is due to unmet demand, in this case, we rerun dispatch including extra variables to
435        // track the unmet demand so we can report the offending regions/commodities to users
436        match self.run_without_unmet_demand(commodities) {
437            Ok(solution) => Ok(solution),
438            Err(ModelError::NonOptimal(HighsModelStatus::Infeasible)) => {
439                let pairs = self
440                    .get_regions_and_commodities_with_unmet_demand(commodities)
441                    .expect("Failed to run dispatch to calculate unmet demand");
442
443                ensure!(
444                    !pairs.is_empty(),
445                    "Model is infeasible, but there was no unmet demand"
446                );
447
448                bail!(
449                    "The solver has indicated that the problem is infeasible, probably because \
450                    the supplied assets could not meet the required demand. Demand was not met \
451                    for the following region and commodity pairs: {}",
452                    format_items_with_cap(pairs)
453                )
454            }
455            Err(err) => Err(err)?,
456        }
457    }
458
459    /// Run dispatch without unmet demand variables
460    fn run_without_unmet_demand(
461        &self,
462        commodities: &[CommodityID],
463    ) -> Result<Solution<'model>, ModelError> {
464        self.run_internal(commodities, /*allow_unmet_demand=*/ false)
465    }
466
467    /// Run dispatch to diagnose which regions and commodities have unmet demand
468    fn get_regions_and_commodities_with_unmet_demand(
469        &self,
470        commodities: &[CommodityID],
471    ) -> Result<IndexSet<(RegionID, CommodityID)>> {
472        let solution = self.run_internal(commodities, /*allow_unmet_demand=*/ true)?;
473        Ok(solution
474            .iter_unmet_demand()
475            .filter(|(_, _, _, flow)| *flow > Flow(0.0))
476            .map(|(commodity_id, region_id, _, _)| (region_id.clone(), commodity_id.clone()))
477            .collect())
478    }
479
480    /// Run dispatch for specified commodities, optionally including unmet demand variables
481    fn run_internal(
482        &self,
483        commodities: &[CommodityID],
484        allow_unmet_demand: bool,
485    ) -> Result<Solution<'model>, ModelError> {
486        // Set up problem
487        let mut problem = Problem::default();
488        let mut variables = VariableMap::new_with_asset_vars(
489            &mut problem,
490            self.model,
491            self.input_prices,
492            self.existing_assets,
493            self.candidate_assets,
494            self.year,
495        );
496
497        // If unmet demand is enabled for this dispatch run (and is allowed by the model param) then
498        // we add variables representing unmet demand
499        if allow_unmet_demand {
500            variables.add_unmet_demand_variables(&mut problem, self.model, commodities);
501        }
502
503        // Add constraints
504        let all_assets = chain(self.existing_assets.iter(), self.candidate_assets.iter());
505        let constraint_keys = add_asset_constraints(
506            &mut problem,
507            &variables,
508            self.model,
509            &all_assets,
510            commodities,
511            self.year,
512        );
513
514        // Solve model
515        let solution = solve_optimal(problem.optimise(Sense::Minimise))?;
516
517        let objective_value = Money(solution.objective_value());
518        debug!("Objective value: {objective_value}");
519
520        Ok(Solution {
521            solution: solution.get_solution(),
522            variables,
523            time_slice_info: &self.model.time_slice_info,
524            constraint_keys,
525            objective_value,
526        })
527    }
528}
529
530/// Add variables to the optimisation problem.
531///
532/// # Arguments
533///
534/// * `problem` - The optimisation problem
535/// * `variables` - The map of asset variables
536/// * `time_slice_info` - Information about assets
537/// * `input_prices` - Optional explicit prices for input commodities
538/// * `assets` - Assets to include
539/// * `year` - Current milestone year
540fn add_asset_variables(
541    problem: &mut Problem,
542    variables: &mut AssetVariableMap,
543    time_slice_info: &TimeSliceInfo,
544    input_prices: Option<&CommodityPrices>,
545    assets: &[AssetRef],
546    year: u32,
547) -> Range<usize> {
548    // This line **must** come before we add more variables
549    let start = problem.num_cols();
550
551    for (asset, time_slice) in iproduct!(assets.iter(), time_slice_info.iter_ids()) {
552        let coeff = calculate_cost_coefficient(asset, year, time_slice, input_prices);
553        let var = problem.add_column(coeff.value(), 0.0..);
554        let key = (asset.clone(), time_slice.clone());
555        let existing = variables.insert(key, var).is_some();
556        assert!(!existing, "Duplicate entry for var");
557    }
558
559    start..problem.num_cols()
560}
561
562/// Calculate the cost coefficient for a decision variable.
563///
564/// Normally, the cost coefficient is the same as the asset's operating costs for the given year and
565/// time slice. If `input_prices` is provided then those prices are added to the flow costs for the
566/// relevant commodities, if they are input flows for the asset.
567///
568/// # Arguments
569///
570/// * `asset` - The asset to calculate the coefficient for
571/// * `year` - The current milestone year
572/// * `time_slice` - The time slice to which this coefficient applies
573/// * `input_prices` - Optional map of prices to include for input commodities
574///
575/// # Returns
576///
577/// The cost coefficient to be used for the relevant decision variable.
578fn calculate_cost_coefficient(
579    asset: &Asset,
580    year: u32,
581    time_slice: &TimeSliceID,
582    input_prices: Option<&CommodityPrices>,
583) -> MoneyPerActivity {
584    let opex = asset.get_operating_cost(year, time_slice);
585    let input_cost = input_prices
586        .map(|prices| asset.get_input_cost_from_prices(prices, time_slice))
587        .unwrap_or_default();
588    opex + input_cost
589}
muse2/simulation/optimisation.rs

muse2/simulation/
optimisation.rs
