diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index c18097c5..29b86645 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -2,15 +2,18 @@ //! //! This is used to calculate commodity flows and prices. use crate::asset::{Asset, AssetID, AssetPool}; -use crate::commodity::{BalanceType, CommodityType}; +use crate::commodity::BalanceType; use crate::model::Model; use crate::process::ProcessFlow; -use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; +use crate::time_slice::{TimeSliceID, TimeSliceInfo}; use anyhow::{anyhow, Result}; use highs::{HighsModelStatus, RowProblem as Problem, Sense}; use indexmap::IndexMap; use std::rc::Rc; +mod constraints; +use constraints::{add_asset_constraints, CapacityConstraintKeys, CommodityBalanceConstraintKeys}; + /// A decision variable in the optimisation /// /// Note that this type does **not** include the value of the variable; it just refers to a @@ -64,12 +67,6 @@ impl VariableMapKey { } } -/// Indicates the commodity ID and time slice selection covered by each commodity balance constraint -type CommodityBalanceConstraintKeys = Vec<(Rc, TimeSliceSelection)>; - -/// Indicates the asset ID and time slice covered by each capacity constraint -type CapacityConstraintKeys = Vec<(AssetID, TimeSliceID)>; - /// The solution to the dispatch optimisation problem pub struct Solution<'a> { solution: highs::Solution, @@ -278,221 +275,6 @@ fn calculate_cost_coefficient( } } -/// Add asset-level constraints -/// -/// Note: the ordering of constraints is important, as the dual values of the constraints must later -/// be retrieved to calculate commodity prices. -/// -/// # Arguments: -/// -/// * `problem` - The optimisation problem -/// * `variables` - The variables in the problem -/// * `model` - The model -/// * `assets` - The asset pool -/// * `year` - Current milestone year -/// -/// # Returns: -/// -/// * A vector of keys for commodity balance constraints -/// * A vector of keys for capacity constraints -fn add_asset_constraints( - problem: &mut Problem, - variables: &VariableMap, - model: &Model, - assets: &AssetPool, - year: u32, -) -> (CommodityBalanceConstraintKeys, CapacityConstraintKeys) { - let commodity_balance_constraint_keys = - add_commodity_balance_constraints(problem, variables, model, assets, year); - - let capacity_constraint_keys = add_asset_capacity_constraints( - problem, - variables, - assets, - &model.time_slice_info, - &commodity_balance_constraint_keys, - ); - - // **TODO**: Currently it's safe to assume all process flows are non-flexible, as we enforce - // this when reading data in. Once we've added support for flexible process flows, we will - // need to add different constraints for assets with flexible and non-flexible flows. - // - // See: https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/360 - add_fixed_asset_constraints(problem, variables, assets, &model.time_slice_info); - - // Return constraint keys which are required to calculate commodity prices - (commodity_balance_constraint_keys, capacity_constraint_keys) -} - -/// Add asset-level input-output commodity balances. -/// -/// These constraints fix the supply-demand balance for the whole system. -/// -/// See description in [the dispatch optimisation documentation][1]. -/// -/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html#commodity-balance-constraints -fn add_commodity_balance_constraints( - problem: &mut Problem, - variables: &VariableMap, - model: &Model, - assets: &AssetPool, - year: u32, -) -> CommodityBalanceConstraintKeys { - // Sanity check: we rely on the first n values of the dual row values corresponding to the - // commodity constraints, so these must be the first rows - assert!( - problem.num_rows() == 0, - "Commodity balance constraints must be added before other constraints" - ); - - let mut terms = Vec::new(); - let mut keys = CommodityBalanceConstraintKeys::new(); - for commodity in model.commodities.values() { - if commodity.kind != CommodityType::SupplyEqualsDemand - && commodity.kind != CommodityType::ServiceDemand - { - continue; - } - - for region_id in model.iter_regions() { - for ts_selection in model - .time_slice_info - .iter_selections_for_level(commodity.time_slice_level) - { - // Note about performance: this loop **may** prove to be a bottleneck as - // `time_slice_info.iter_selection` returns a `Box` and so requires a heap - // allocation each time. For commodities with a `TimeSliceLevel` of `TimeSlice` (the - // worst case), this means the number of additional heap allocations will equal the - // number of time slices, which for this function could be in the - // hundreds/thousands. - for (time_slice, _) in model.time_slice_info.iter_selection(&ts_selection) { - // Add terms for this asset + commodity at this time slice. The coefficient for - // each variable is one. - terms.extend( - assets - .iter_for_region_and_commodity(region_id, commodity) - .map(|asset| (variables.get(asset.id, &commodity.id, time_slice), 1.0)), - ); - } - - // Get the RHS of the equation for a commodity balance constraint. For SED - // commodities, the RHS will be zero and for SVD commodities it will be equal to the - // demand for the given time slice selection. - let rhs = match commodity.kind { - CommodityType::SupplyEqualsDemand => 0.0, - CommodityType::ServiceDemand => { - match ts_selection { - TimeSliceSelection::Single(ref ts) => { - commodity.demand.get(region_id, year, ts) - } - // We currently only support specifying demand at the time slice level: - // https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/391 - _ => panic!( - "Currently SVD commodities must have a time slice level of time slice" - ), - } - } - _ => unreachable!(), - }; - - // Add constraint (sum of terms must equal rhs) - problem.add_row(rhs..=rhs, terms.drain(0..)); - - // Keep track of the order in which constraints were added - keys.push((Rc::clone(&commodity.id), ts_selection)); - } - } - } - - keys -} - -/// Add constraints for non-flexible assets. -/// -/// Non-flexible assets are those which have a fixed ratio between inputs and outputs. -/// -/// See description in [the dispatch optimisation documentation][1]. -/// -/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html#non-flexible-assets -fn add_fixed_asset_constraints( - problem: &mut Problem, - variables: &VariableMap, - assets: &AssetPool, - time_slice_info: &TimeSliceInfo, -) { - for asset in assets.iter() { - // Get first PAC. unwrap is safe because all processes have at least one PAC. - let pac1 = asset.process.iter_pacs().next().unwrap(); - - for time_slice in time_slice_info.iter_ids() { - let pac_var = variables.get(asset.id, &pac1.commodity.id, time_slice); - let pac_term = (pac_var, -1.0 / pac1.flow); - - for flow in asset.process.flows.iter() { - // Don't add a constraint for the PAC itself - if Rc::ptr_eq(&flow.commodity, &pac1.commodity) { - continue; - } - - // We are enforcing that (var / flow) - (pac_var / pac_flow) = 0 - let var = variables.get(asset.id, &flow.commodity.id, time_slice); - problem.add_row(0.0..=0.0, [(var, 1.0 / flow.flow), pac_term]); - } - } - } -} - -/// Add asset-level capacity and availability constraints. -/// -/// For every asset at every time slice, the sum of the commodity flows for PACs must not exceed the -/// capacity limits, which are a product of the annual capacity, time slice length and process -/// availability. -/// -/// See description in [the dispatch optimisation documentation][1]. -/// -/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html#asset-level-capacity-and-availability-constraints -fn add_asset_capacity_constraints( - problem: &mut Problem, - variables: &VariableMap, - assets: &AssetPool, - time_slice_info: &TimeSliceInfo, - commodity_balance_constraint_keys: &CommodityBalanceConstraintKeys, -) -> CapacityConstraintKeys { - // Sanity check: we rely on the dual rows corresponding to the capacity constraints being - // immediately after the commodity balance constraints in `problem`. - assert!( - problem.num_rows() == commodity_balance_constraint_keys.len(), - "Capacity constraints must be added immediately after commodity constraints." - ); - - let mut terms = Vec::new(); - let mut keys = CapacityConstraintKeys::new(); - for asset in assets.iter() { - for time_slice in time_slice_info.iter_ids() { - let mut is_input = false; // NB: there will be at least one PAC - for flow in asset.process.iter_pacs() { - is_input = flow.flow < 0.0; // NB: PACs will be all inputs or all outputs - - let var = variables.get(asset.id, &flow.commodity.id, time_slice); - terms.push((var, 1.0)); - } - - let mut limits = asset.get_activity_limits(time_slice); - - // If it's an input flow, the q's will be negative, so we need to invert the limits - if is_input { - limits = -limits.end()..=-limits.start(); - } - - problem.add_row(limits, terms.drain(0..)); - - // Keep track of the order in which constraints were added - keys.push((asset.id, time_slice.clone())); - } - } - keys -} - #[cfg(test)] mod tests { use super::*; diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs new file mode 100644 index 00000000..f77ea45c --- /dev/null +++ b/src/simulation/optimisation/constraints.rs @@ -0,0 +1,230 @@ +//! Code for adding constraints to the dispatch optimisation problem. +use crate::asset::{AssetID, AssetPool}; +use crate::commodity::CommodityType; +use crate::model::Model; +use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; +use highs::RowProblem as Problem; +use std::rc::Rc; + +use super::VariableMap; + +/// Indicates the commodity ID and time slice selection covered by each commodity balance constraint +pub type CommodityBalanceConstraintKeys = Vec<(Rc, TimeSliceSelection)>; + +/// Indicates the asset ID and time slice covered by each capacity constraint +pub type CapacityConstraintKeys = Vec<(AssetID, TimeSliceID)>; + +/// Add asset-level constraints +/// +/// Note: the ordering of constraints is important, as the dual values of the constraints must later +/// be retrieved to calculate commodity prices. +/// +/// # Arguments: +/// +/// * `problem` - The optimisation problem +/// * `variables` - The variables in the problem +/// * `model` - The model +/// * `assets` - The asset pool +/// * `year` - Current milestone year +/// +/// # Returns: +/// +/// * A vector of keys for commodity balance constraints +/// * A vector of keys for capacity constraints +pub fn add_asset_constraints( + problem: &mut Problem, + variables: &VariableMap, + model: &Model, + assets: &AssetPool, + year: u32, +) -> (CommodityBalanceConstraintKeys, CapacityConstraintKeys) { + let commodity_balance_constraint_keys = + add_commodity_balance_constraints(problem, variables, model, assets, year); + + let capacity_constraint_keys = add_asset_capacity_constraints( + problem, + variables, + assets, + &model.time_slice_info, + &commodity_balance_constraint_keys, + ); + + // **TODO**: Currently it's safe to assume all process flows are non-flexible, as we enforce + // this when reading data in. Once we've added support for flexible process flows, we will + // need to add different constraints for assets with flexible and non-flexible flows. + // + // See: https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/360 + add_fixed_asset_constraints(problem, variables, assets, &model.time_slice_info); + + // Return constraint keys which are required to calculate commodity prices + (commodity_balance_constraint_keys, capacity_constraint_keys) +} + +/// Add asset-level input-output commodity balances. +/// +/// These constraints fix the supply-demand balance for the whole system. +/// +/// See description in [the dispatch optimisation documentation][1]. +/// +/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html#commodity-balance-constraints +fn add_commodity_balance_constraints( + problem: &mut Problem, + variables: &VariableMap, + model: &Model, + assets: &AssetPool, + year: u32, +) -> CommodityBalanceConstraintKeys { + // Sanity check: we rely on the first n values of the dual row values corresponding to the + // commodity constraints, so these must be the first rows + assert!( + problem.num_rows() == 0, + "Commodity balance constraints must be added before other constraints" + ); + + let mut terms = Vec::new(); + let mut keys = CommodityBalanceConstraintKeys::new(); + for commodity in model.commodities.values() { + if commodity.kind != CommodityType::SupplyEqualsDemand + && commodity.kind != CommodityType::ServiceDemand + { + continue; + } + + for region_id in model.iter_regions() { + for ts_selection in model + .time_slice_info + .iter_selections_for_level(commodity.time_slice_level) + { + // Note about performance: this loop **may** prove to be a bottleneck as + // `time_slice_info.iter_selection` returns a `Box` and so requires a heap + // allocation each time. For commodities with a `TimeSliceLevel` of `TimeSlice` (the + // worst case), this means the number of additional heap allocations will equal the + // number of time slices, which for this function could be in the + // hundreds/thousands. + for (time_slice, _) in model.time_slice_info.iter_selection(&ts_selection) { + // Add terms for this asset + commodity at this time slice. The coefficient for + // each variable is one. + terms.extend( + assets + .iter_for_region_and_commodity(region_id, commodity) + .map(|asset| (variables.get(asset.id, &commodity.id, time_slice), 1.0)), + ); + } + + // Get the RHS of the equation for a commodity balance constraint. For SED + // commodities, the RHS will be zero and for SVD commodities it will be equal to the + // demand for the given time slice selection. + let rhs = match commodity.kind { + CommodityType::SupplyEqualsDemand => 0.0, + CommodityType::ServiceDemand => { + match ts_selection { + TimeSliceSelection::Single(ref ts) => { + commodity.demand.get(region_id, year, ts) + } + // We currently only support specifying demand at the time slice level: + // https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/391 + _ => panic!( + "Currently SVD commodities must have a time slice level of time slice" + ), + } + } + _ => unreachable!(), + }; + + // Add constraint (sum of terms must equal rhs) + problem.add_row(rhs..=rhs, terms.drain(0..)); + + // Keep track of the order in which constraints were added + keys.push((Rc::clone(&commodity.id), ts_selection)); + } + } + } + + keys +} + +/// Add constraints for non-flexible assets. +/// +/// Non-flexible assets are those which have a fixed ratio between inputs and outputs. +/// +/// See description in [the dispatch optimisation documentation][1]. +/// +/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html#non-flexible-assets +fn add_fixed_asset_constraints( + problem: &mut Problem, + variables: &VariableMap, + assets: &AssetPool, + time_slice_info: &TimeSliceInfo, +) { + for asset in assets.iter() { + // Get first PAC. unwrap is safe because all processes have at least one PAC. + let pac1 = asset.process.iter_pacs().next().unwrap(); + + for time_slice in time_slice_info.iter_ids() { + let pac_var = variables.get(asset.id, &pac1.commodity.id, time_slice); + let pac_term = (pac_var, -1.0 / pac1.flow); + + for flow in asset.process.flows.iter() { + // Don't add a constraint for the PAC itself + if Rc::ptr_eq(&flow.commodity, &pac1.commodity) { + continue; + } + + // We are enforcing that (var / flow) - (pac_var / pac_flow) = 0 + let var = variables.get(asset.id, &flow.commodity.id, time_slice); + problem.add_row(0.0..=0.0, [(var, 1.0 / flow.flow), pac_term]); + } + } + } +} + +/// Add asset-level capacity and availability constraints. +/// +/// For every asset at every time slice, the sum of the commodity flows for PACs must not exceed the +/// capacity limits, which are a product of the annual capacity, time slice length and process +/// availability. +/// +/// See description in [the dispatch optimisation documentation][1]. +/// +/// [1]: https://energysystemsmodellinglab.github.io/MUSE_2.0/dispatch_optimisation.html#asset-level-capacity-and-availability-constraints +fn add_asset_capacity_constraints( + problem: &mut Problem, + variables: &VariableMap, + assets: &AssetPool, + time_slice_info: &TimeSliceInfo, + commodity_balance_constraint_keys: &CommodityBalanceConstraintKeys, +) -> CapacityConstraintKeys { + // Sanity check: we rely on the dual rows corresponding to the capacity constraints being + // immediately after the commodity balance constraints in `problem`. + assert!( + problem.num_rows() == commodity_balance_constraint_keys.len(), + "Capacity constraints must be added immediately after commodity constraints." + ); + + let mut terms = Vec::new(); + let mut keys = CapacityConstraintKeys::new(); + for asset in assets.iter() { + for time_slice in time_slice_info.iter_ids() { + let mut is_input = false; // NB: there will be at least one PAC + for flow in asset.process.iter_pacs() { + is_input = flow.flow < 0.0; // NB: PACs will be all inputs or all outputs + + let var = variables.get(asset.id, &flow.commodity.id, time_slice); + terms.push((var, 1.0)); + } + + let mut limits = asset.get_activity_limits(time_slice); + + // If it's an input flow, the q's will be negative, so we need to invert the limits + if is_input { + limits = -limits.end()..=-limits.start(); + } + + problem.add_row(limits, terms.drain(0..)); + + // Keep track of the order in which constraints were added + keys.push((asset.id, time_slice.clone())); + } + } + keys +}