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}