muse2/simulation/optimisation/
constraints.rs

1//! Code for adding constraints to the dispatch optimisation problem.
2use super::VariableMap;
3use crate::asset::{AssetIterator, AssetRef};
4use crate::commodity::{CommodityID, CommodityType};
5use crate::model::Model;
6use crate::region::RegionID;
7use crate::time_slice::{TimeSliceID, TimeSliceSelection};
8use crate::units::UnitType;
9use highs::RowProblem as Problem;
10
11/// Corresponding variables for a constraint along with the row offset in the solution
12pub struct KeysWithOffset<T> {
13    offset: usize,
14    keys: Vec<T>,
15}
16
17impl<T> KeysWithOffset<T> {
18    /// Zip the keys with the corresponding dual values in the solution, accounting for the offset
19    pub fn zip_duals<'a, U>(&'a self, duals: &'a [f64]) -> impl Iterator<Item = (&'a T, U)>
20    where
21        U: UnitType,
22    {
23        assert!(
24            self.offset + self.keys.len() <= duals.len(),
25            "Bad constraint keys: dual rows out of range"
26        );
27
28        self.keys
29            .iter()
30            .zip(duals[self.offset..].iter().copied().map(U::new))
31    }
32}
33
34/// Indicates the commodity ID and time slice selection covered by each commodity balance constraint
35pub type CommodityBalanceKeys = KeysWithOffset<(CommodityID, RegionID, TimeSliceSelection)>;
36
37/// Indicates the asset ID and time slice covered by each activity constraint
38pub type ActivityKeys = KeysWithOffset<(AssetRef, TimeSliceID)>;
39
40/// The keys for different constraints
41pub struct ConstraintKeys {
42    /// Keys for commodity balance constraints
43    pub commodity_balance_keys: CommodityBalanceKeys,
44    /// Keys for activity constraints
45    pub activity_keys: ActivityKeys,
46}
47
48/// Add asset-level constraints
49///
50/// Note: the ordering of constraints is important, as the dual values of the constraints must later
51/// be retrieved to calculate commodity prices.
52///
53/// # Arguments
54///
55/// * `problem` - The optimisation problem
56/// * `variables` - The variables in the problem
57/// * `model` - The model
58/// * `assets` - The asset pool
59/// * `commodities` - The subset of commodities to apply constraints to
60/// * `year` - Current milestone year
61///
62/// # Returns
63///
64/// Keys for the different constraints.
65pub fn add_asset_constraints<'a, I>(
66    problem: &mut Problem,
67    variables: &VariableMap,
68    model: &'a Model,
69    assets: &I,
70    commodities: &'a [CommodityID],
71    year: u32,
72) -> ConstraintKeys
73where
74    I: Iterator<Item = &'a AssetRef> + Clone + 'a,
75{
76    let commodity_balance_keys =
77        add_commodity_balance_constraints(problem, variables, model, assets, commodities, year);
78
79    let activity_keys = add_activity_constraints(problem, variables);
80
81    // Return constraint keys
82    ConstraintKeys {
83        commodity_balance_keys,
84        activity_keys,
85    }
86}
87
88/// Add asset-level input-output commodity balances.
89///
90/// These constraints fix the supply-demand balance for the whole system.
91///
92/// See description in [the dispatch optimisation documentation][1].
93///
94/// [1]: https://energysystemsmodellinglab.github.io/MUSE2/model/dispatch_optimisation.html#commodity-balance-for--cin-mathbfcmathrmsed-
95fn add_commodity_balance_constraints<'a, I>(
96    problem: &mut Problem,
97    variables: &VariableMap,
98    model: &'a Model,
99    assets: &I,
100    commodities: &'a [CommodityID],
101    year: u32,
102) -> CommodityBalanceKeys
103where
104    I: Iterator<Item = &'a AssetRef> + Clone + 'a,
105{
106    // Row offset in problem. This line **must** come before we add more constraints.
107    let offset = problem.num_rows();
108
109    let mut keys = Vec::new();
110    let mut terms = Vec::new();
111    for commodity_id in commodities {
112        let commodity = &model.commodities[commodity_id];
113        if !matches!(
114            commodity.kind,
115            CommodityType::SupplyEqualsDemand | CommodityType::ServiceDemand
116        ) {
117            continue;
118        }
119
120        for region_id in model.iter_regions() {
121            for ts_selection in model
122                .time_slice_info
123                .iter_selections_at_level(commodity.time_slice_level)
124            {
125                for (asset, flow) in assets
126                    .clone()
127                    .filter_region(region_id)
128                    .flows_for_commodity(commodity_id)
129                {
130                    // If the commodity has a time slice level of season/annual, the constraint will
131                    // cover multiple time slices
132                    for (time_slice, _) in ts_selection.iter(&model.time_slice_info) {
133                        let var = variables.get_asset_var(asset, time_slice);
134                        terms.push((var, flow.coeff.value()));
135                    }
136                }
137
138                // It is possible that a commodity may not be produced or consumed by anything in a
139                // given milestone year, in which case it doesn't make sense to add a commodity
140                // balance constraint
141                if terms.is_empty() {
142                    continue;
143                }
144
145                // Also include unmet demand variables if required
146                if !variables.unmet_demand_var_idx.is_empty() {
147                    for (time_slice, _) in ts_selection.iter(&model.time_slice_info) {
148                        let var =
149                            variables.get_unmet_demand_var(commodity_id, region_id, time_slice);
150                        terms.push((var, 1.0));
151                    }
152                }
153
154                // For SED commodities, the LHS must be >=0 and for SVD commodities, it must be >=
155                // the exogenous demand supplied by the user
156                let min = if commodity.kind == CommodityType::ServiceDemand {
157                    commodity
158                        .demand
159                        .get(&(region_id.clone(), year, ts_selection.clone()))
160                        .unwrap()
161                        .value()
162                } else {
163                    0.0
164                };
165                problem.add_row(min.., terms.drain(..));
166                keys.push((
167                    commodity_id.clone(),
168                    region_id.clone(),
169                    ts_selection.clone(),
170                ));
171            }
172        }
173    }
174
175    CommodityBalanceKeys { offset, keys }
176}
177
178/// Add constraints on the activity of different assets.
179///
180/// This ensures that assets do not exceed their specified capacity and availability for each time
181/// slice.
182///
183/// See description in [the dispatch optimisation documentation][1].
184///
185/// [1]: https://energysystemsmodellinglab.github.io/MUSE2/model/dispatch_optimisation.html#a4-constraints-capacity--availability-for-standard-assets--a-in-mathbfastd-
186fn add_activity_constraints(problem: &mut Problem, variables: &VariableMap) -> ActivityKeys {
187    // Row offset in problem. This line **must** come before we add more constraints.
188    let offset = problem.num_rows();
189
190    let mut keys = Vec::new();
191    for (asset, time_slice, var) in variables.iter_asset_vars() {
192        let limits = asset.get_activity_limits(time_slice);
193        let limits = limits.start().value()..=limits.end().value();
194
195        problem.add_row(limits, [(var, 1.0)]);
196        keys.push((asset.clone(), time_slice.clone()));
197    }
198
199    ActivityKeys { offset, keys }
200}
muse2/simulation/optimisation/constraints.rs

muse2/simulation/optimisation/
constraints.rs
