Skip to content

Move constraints code to new submodule #470

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tsmbland merged 1 commit into from
Apr 4, 2025
Merged

Move constraints code to new submodule #470

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
228 changes: 5 additions & 223 deletions src/simulation/optimisation.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -64,12 +67,6 @@ impl VariableMapKey {
}
}

/// Indicates the commodity ID and time slice selection covered by each commodity balance constraint
type CommodityBalanceConstraintKeys = Vec<(Rc<str>, 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,
Expand Down Expand Up @@ -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::*;
Expand Down
Loading