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