Merge pull request #628 from EnergySystemsModellingLab/add-constraints

Add commodity balance and demand satisfaction constraints

Author: alexdewar

src/asset.rs

@@ -264,20 +264,15 @@ impl AssetPool {
         self.iter().filter(|asset| asset.region_id == *region_id)
     }
 
-    /// Iterate over only the active assets in a given region that produce or consume a given
-    /// commodity
+    /// Iterate over the active assets in a given region that produce/consume a commodity with the
+    /// associated process flow
     pub fn iter_for_region_and_commodity<'a>(
         &'a self,
         region_id: &'a RegionID,
         commodity_id: &'a CommodityID,
-    ) -> impl Iterator<Item = &'a AssetRef> {
-        self.iter_for_region(region_id).filter(|asset| {
-            asset.process.contains_commodity_flow(
-                commodity_id,
-                &asset.region_id,
-                asset.commission_year,
-            )
-        })
+    ) -> impl Iterator<Item = (&'a AssetRef, &'a ProcessFlow)> {
+        self.iter_for_region(region_id)
+            .filter_map(|asset| Some((asset, asset.get_flow(commodity_id)?)))
     }
 
     /// Replace the active pool with new and/or already commissioned assets

src/input/commodity/demand.rs

@@ -187,7 +187,7 @@ fn compute_demand_maps(
 ) -> HashMap<CommodityID, DemandMap> {
     let mut map = HashMap::new();
     for ((commodity_id, region_id, year), (level, annual_demand)) in demand.iter() {
-        for ts_selection in time_slice_info.iter_selections_for_level(*level) {
+        for ts_selection in time_slice_info.iter_selections_at_level(*level) {
             let slice_key = (
                 commodity_id.clone(),
                 region_id.clone(),

src/input/commodity/demand_slicing.rs

@@ -128,7 +128,7 @@ fn validate_demand_slices(
 ) -> Result<()> {
     for (commodity, region_id) in iproduct!(svd_commodities.values(), region_ids) {
         time_slice_info
-            .iter_selections_for_level(commodity.time_slice_level)
+            .iter_selections_at_level(commodity.time_slice_level)
             .map(|ts_selection| {
                 demand_slices
                     .get(&(

src/input/process.rs

@@ -174,7 +174,7 @@ fn validate_commodities(
             }
             CommodityType::ServiceDemand => {
                 for ts_selection in
-                    time_slice_info.iter_selections_for_level(commodity.time_slice_level)
+                    time_slice_info.iter_selections_at_level(commodity.time_slice_level)
                 {
                     validate_svd_commodity(
                         time_slice_info,

src/output.rs

@@ -130,17 +130,6 @@ struct CommodityBalanceDualsRow {
     value: f64,
 }
 
-/// Represents the fixed asset duals data in a row of the fixed asset duals CSV file
-#[derive(Serialize, Deserialize, Debug, PartialEq)]
-struct FixedAssetDualsRow {
-    pac: CommodityID,
-    pac_flow: f64,
-    commodity_id: CommodityID,
-    commodity_flow: f64,
-    time_slice: TimeSliceID,
-    value: f64,
-}
-
 /// For writing extra debug information about the model
 struct DebugDataWriter {
     commodity_balance_duals_writer: csv::Writer<File>,

src/process.rs

@@ -48,21 +48,6 @@ pub struct Process {
     pub regions: HashSet<RegionID>,
 }
 
-impl Process {
-    /// Whether the process contains a flow for a given commodity
-    pub fn contains_commodity_flow(
-        &self,
-        commodity_id: &CommodityID,
-        region_id: &RegionID,
-        year: u32,
-    ) -> bool {
-        self.flows
-            .get(&(region_id.clone(), year))
-            .unwrap() // all regions and years are covered
-            .contains_key(commodity_id)
-    }
-}
-
 /// Represents a maximum annual commodity coeff for a given process
 #[derive(PartialEq, Debug, Clone)]
 pub struct ProcessFlow {

src/simulation/optimisation/constraints.rs

@@ -1,7 +1,7 @@
 //! Code for adding constraints to the dispatch optimisation problem.
 use super::VariableMap;
 use crate::asset::{AssetPool, AssetRef};
-use crate::commodity::CommodityID;
+use crate::commodity::{CommodityID, CommodityType};
 use crate::model::Model;
 use crate::region::RegionID;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection};
@@ -85,18 +85,57 @@ pub fn add_asset_constraints(
 /// [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,
+    variables: &VariableMap,
+    model: &Model,
+    assets: &AssetPool,
+    year: u32,
 ) -> CommodityBalanceKeys {
     // Row offset in problem. This line **must** come before we add more constraints.
     let offset = problem.num_rows();
 
-    let keys = Vec::new();
-
-    // **TODO:** Add commodity balance constraints:
-    //  https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/577
+    let mut keys = Vec::new();
+    let mut terms = Vec::new();
+    for (commodity_id, commodity) in model.commodities.iter() {
+        if !matches!(
+            commodity.kind,
+            CommodityType::SupplyEqualsDemand | CommodityType::ServiceDemand
+        ) {
+            continue;
+        }
+
+        for region_id in model.iter_regions() {
+            for ts_selection in model
+                .time_slice_info
+                .iter_selections_at_level(commodity.time_slice_level)
+            {
+                for (asset, flow) in assets.iter_for_region_and_commodity(region_id, commodity_id) {
+                    // If the commodity has a time slice level of season/annual, the constraint will
+                    // cover multiple time slices
+                    for (time_slice, _) in ts_selection.iter(&model.time_slice_info) {
+                        let var = variables.get(asset, time_slice);
+                        terms.push((var, flow.coeff));
+                    }
+                }
+
+                // Add constraint. For SED commodities, the RHS is zero and for SVD commodities it
+                // is the exogenous demand supplied by the user.
+                let rhs = if commodity.kind == CommodityType::ServiceDemand {
+                    *commodity
+                        .demand
+                        .get(&(region_id.clone(), year, ts_selection.clone()))
+                        .unwrap()
+                } else {
+                    0.0
+                };
+                problem.add_row(rhs..=rhs, terms.drain(..));
+                keys.push((
+                    commodity_id.clone(),
+                    region_id.clone(),
+                    ts_selection.clone(),
+                ))
+            }
+        }
+    }
 
     CommodityBalanceKeys { offset, keys }
 }

src/time_slice.rs

@@ -292,7 +292,7 @@ impl TimeSliceInfo {
     ///
     /// For example, if [`TimeSliceLevel::Season`] is specified, this function will return an
     /// iterator of [`TimeSliceSelection`]s covering each season.
-    pub fn iter_selections_for_level(
+    pub fn iter_selections_at_level(
         &self,
         level: TimeSliceLevel,
     ) -> Box<dyn Iterator<Item = TimeSliceSelection> + '_> {