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 = add_commodity_balance_constraints(
77 problem,
78 variables,
79 model,
80 assets.clone(),
81 commodities,
82 year,
83 );
84
85 let activity_keys = add_activity_constraints(problem, variables);
86
87 // Return constraint keys
88 ConstraintKeys {
89 commodity_balance_keys,
90 activity_keys,
91 }
92}
93
94/// Add asset-level input-output commodity balances.
95///
96/// These constraints fix the supply-demand balance for the whole system.
97///
98/// See description in [the dispatch optimisation documentation][1].
99///
100/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/model/dispatch_optimisation.html#commodity-balance-for--cin-mathbfcmathrmsed-
101fn add_commodity_balance_constraints<'a, I>(
102 problem: &mut Problem,
103 variables: &VariableMap,
104 model: &'a Model,
105 assets: I,
106 commodities: &'a [CommodityID],
107 year: u32,
108) -> CommodityBalanceKeys
109where
110 I: Iterator<Item = &'a AssetRef> + Clone + 'a,
111{
112 // Row offset in problem. This line **must** come before we add more constraints.
113 let offset = problem.num_rows();
114
115 let mut keys = Vec::new();
116 let mut terms = Vec::new();
117 for commodity_id in commodities {
118 let commodity = &model.commodities[commodity_id];
119 if !matches!(
120 commodity.kind,
121 CommodityType::SupplyEqualsDemand | CommodityType::ServiceDemand
122 ) {
123 continue;
124 }
125
126 for region_id in model.iter_regions() {
127 for ts_selection in model
128 .time_slice_info
129 .iter_selections_at_level(commodity.time_slice_level)
130 {
131 for (asset, flow) in assets
132 .clone()
133 .filter_region(region_id)
134 .flows_for_commodity(commodity_id)
135 {
136 // If the commodity has a time slice level of season/annual, the constraint will
137 // cover multiple time slices
138 for (time_slice, _) in ts_selection.iter(&model.time_slice_info) {
139 let var = variables.get(asset, time_slice);
140 terms.push((var, flow.coeff.value()));
141 }
142 }
143
144 // It is possible that a commodity may not be produced or consumed by anything in a
145 // given milestone year, in which case it doesn't make sense to add a commodity
146 // balance constraint
147 if terms.is_empty() {
148 continue;
149 }
150
151 // Add constraint. For SED commodities, the RHS is zero and for SVD commodities it
152 // is the exogenous demand supplied by the user.
153 let rhs = if commodity.kind == CommodityType::ServiceDemand {
154 commodity
155 .demand
156 .get(&(region_id.clone(), year, ts_selection.clone()))
157 .unwrap()
158 .value()
159 } else {
160 0.0
161 };
162 problem.add_row(rhs..=rhs, terms.drain(..));
163 keys.push((
164 commodity_id.clone(),
165 region_id.clone(),
166 ts_selection.clone(),
167 ))
168 }
169 }
170 }
171
172 CommodityBalanceKeys { offset, keys }
173}
174
175/// Add constraints on the activity of different assets.
176///
177/// This ensures that assets do not exceed their specified capacity and availability for each time
178/// slice.
179///
180/// See description in [the dispatch optimisation documentation][1].
181///
182/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/model/dispatch_optimisation.html#a4-constraints-capacity--availability-for-standard-assets--a-in-mathbfastd-
183fn add_activity_constraints(problem: &mut Problem, variables: &VariableMap) -> ActivityKeys {
184 // Row offset in problem. This line **must** come before we add more constraints.
185 let offset = problem.num_rows();
186
187 let mut keys = Vec::new();
188 for (asset, time_slice, var) in variables.iter() {
189 let limits = asset.get_activity_limits(time_slice);
190 let limits = limits.start().value()..=limits.end().value();
191
192 problem.add_row(limits, [(var, 1.0)]);
193 keys.push((asset.clone(), time_slice.clone()))
194 }
195
196 ActivityKeys { offset, keys }
197}