Merge pull request #611 from EnergySystemsModellingLab/commodity-flows

Add code for getting commodity flows

Author: alexdewar

src/output.rs

@@ -4,7 +4,7 @@ use crate::asset::{Asset, AssetRef};
 use crate::commodity::CommodityID;
 use crate::process::ProcessID;
 use crate::region::RegionID;
-use crate::simulation::optimisation::Solution;
+use crate::simulation::optimisation::{FlowMap, Solution};
 use crate::simulation::CommodityPrices;
 use crate::time_slice::TimeSliceID;
 use anyhow::{Context, Result};
@@ -269,16 +269,13 @@ impl DataWriter {
     }
 
     /// Write commodity flows to a CSV file
-    pub fn write_flows<'a, I>(&mut self, milestone_year: u32, flows: I) -> Result<()>
-    where
-        I: Iterator<Item = (&'a AssetRef, &'a CommodityID, &'a TimeSliceID, f64)>,
-    {
-        for (asset, commodity_id, time_slice, flow) in flows {
+    pub fn write_flows(&mut self, milestone_year: u32, flow_map: &FlowMap) -> Result<()> {
+        for ((asset, commodity_id, time_slice), flow) in flow_map {
             let asset_row = AssetRow::new(milestone_year, asset);
             let flow_row = CommodityFlowRow {
                 commodity_id: commodity_id.clone(),
                 time_slice: time_slice.clone(),
-                flow,
+                flow: *flow,
             };
             self.flows_writer.serialize((asset_row, flow_row))?;
         }
@@ -330,6 +327,7 @@ mod tests {
     use crate::asset::AssetPool;
     use crate::fixture::{assets, commodity_id, region_id, time_slice};
     use crate::time_slice::TimeSliceID;
+    use indexmap::indexmap;
     use itertools::{assert_equal, Itertools};
     use rstest::rstest;
     use std::iter;
@@ -362,15 +360,15 @@ mod tests {
     fn test_write_flows(assets: AssetPool, commodity_id: CommodityID, time_slice: TimeSliceID) {
         let milestone_year = 2020;
         let asset = assets.iter().next().unwrap();
-        let flow_item = (asset, &commodity_id, &time_slice, 42.0);
+        let flow_map = indexmap! {
+            (asset.clone(), commodity_id.clone(), time_slice.clone()) => 42.0
+        };
 
         // Write a flow
         let dir = tempdir().unwrap();
         {
             let mut writer = DataWriter::create(dir.path(), false).unwrap();
-            writer
-                .write_flows(milestone_year, iter::once(flow_item))
-                .unwrap();
+            writer.write_flows(milestone_year, &flow_map).unwrap();
             writer.flush().unwrap();
         }
 

src/simulation.rs

@@ -29,22 +29,17 @@ pub fn run(
 ) -> Result<()> {
     let mut writer = DataWriter::create(output_path, debug_model)?;
 
-    let mut opt_results = None; // all results of dispatch optimisation
-    for year in model.iter_years() {
-        info!("Milestone year: {year}");
-
-        // Assets that have been decommissioned cannot be selected by agents
-        assets.decommission_old(year);
+    // Iterate over milestone years
+    let mut year_iter = model.iter_years();
+    let mut year = year_iter.next().unwrap(); // NB: There will be at least one year
 
-        // NB: Agent investment is not carried out in first milestone year
-        if let Some((solution, prices)) = opt_results {
-            perform_agent_investment(&model, &solution, &prices, &mut assets);
+    // There shouldn't be assets already commissioned, but let's do this just in case
+    assets.decommission_old(year);
 
-            // **TODO:** Remove this when we implement at least some of the agent investment code
-            //   See: https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/304
-            error!("Agent investment is not yet implemented. Exiting...");
-            return Ok(());
-        }
+    // **TODO:** Remove annotation when the loop actually loops
+    #[allow(clippy::never_loop)]
+    loop {
+        info!("Milestone year: {year}");
 
         // Newly commissioned assets will be included in optimisation for at least one milestone
         // year before agents have the option of decommissioning them
@@ -56,14 +51,31 @@ pub fn run(
 
         // Dispatch optimisation
         let solution = perform_dispatch_optimisation(&model, &assets, year)?;
+        let flow_map = solution.create_flow_map();
         let prices = CommodityPrices::from_model_and_solution(&model, &solution, &assets);
 
         // Write result of dispatch optimisation to file
         writer.write_debug_info(year, &solution)?;
-        writer.write_flows(year, solution.iter_commodity_flows_for_assets())?;
+        writer.write_flows(year, &flow_map)?;
         writer.write_prices(year, &prices)?;
 
-        opt_results = Some((solution, prices));
+        if let Some(next_year) = year_iter.next() {
+            year = next_year;
+
+            // NB: Agent investment is not carried out in first milestone year
+            perform_agent_investment(&model, &flow_map, &prices, &mut assets);
+
+            // Decommission assets whose lifetime has passed
+            assets.decommission_old(year);
+        } else {
+            // No more milestone years. Simulation is finished.
+            break;
+        }
+
+        // **TODO:** Remove this when we implement at least some of the agent investment code
+        //   See: https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/304
+        error!("Agent investment is not yet implemented. Exiting...");
+        break;
     }
 
     writer.flush()?;

src/simulation/investment.rs

@@ -1,5 +1,5 @@
 //! Code for performing agent investment.
-use super::optimisation::Solution;
+use super::optimisation::FlowMap;
 use super::CommodityPrices;
 use crate::asset::AssetPool;
 use crate::model::Model;
@@ -10,12 +10,12 @@ use log::info;
 /// # Arguments
 ///
 /// * `model` - The model
-/// * `solution` - The solution to the dispatch optimisation
+/// * `flow_map` - Map of commodity flows
 /// * `prices` - Commodity prices
 /// * `assets` - The asset pool
 pub fn perform_agent_investment(
     _model: &Model,
-    _solution: &Solution,
+    _flow_map: &FlowMap,
     _prices: &CommodityPrices,
     assets: &mut AssetPool,
 ) {

src/simulation/optimisation.rs

@@ -13,6 +13,9 @@ use indexmap::IndexMap;
 mod constraints;
 use constraints::{add_asset_constraints, ConstraintKeys};
 
+/// A map of commodity flows calculated during the optimisation
+pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), f64>;
+
 /// A decision variable in the optimisation
 ///
 /// Note that this type does **not** include the value of the variable; it just refers to a
@@ -48,26 +51,30 @@ impl VariableMap {
 /// The solution to the dispatch optimisation problem
 pub struct Solution<'a> {
     solution: highs::Solution,
-    _variables: VariableMap,
+    variables: VariableMap,
     time_slice_info: &'a TimeSliceInfo,
     constraint_keys: ConstraintKeys,
 }
 
 impl Solution<'_> {
-    /// Iterate over the newly calculated commodity flows for assets.
+    /// Create a map of commodity flows for each asset's coeffs at every time slice.
     ///
     /// Note that this only includes commodity flows which relate to assets, so not every commodity
     /// in the simulation will necessarily be represented.
-    ///
-    /// # Returns
-    ///
-    /// An iterator of tuples containing an asset ID, commodity, time slice and flow.
-    pub fn iter_commodity_flows_for_assets(
-        &self,
-    ) -> impl Iterator<Item = (&AssetRef, &CommodityID, &TimeSliceID, f64)> {
-        // **TODO:** Need to calculate flows by multiplying coeffs by asset activity:
-        //  https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/593
-        std::iter::empty()
+    pub fn create_flow_map(&self) -> FlowMap {
+        // The decision variables represent assets' activity levels, not commodity flows. We
+        // multiply this value by the flow coeffs to get commodity flows.
+        let mut flows = FlowMap::new();
+        for ((asset, time_slice), activity) in self.variables.0.keys().zip(self.solution.columns())
+        {
+            for flow in asset.iter_flows() {
+                let flow_key = (asset.clone(), flow.commodity.id.clone(), time_slice.clone());
+                let flow_value = activity * flow.coeff;
+                flows.insert(flow_key, flow_value);
+            }
+        }
+
+        flows
     }
 
     /// Keys and dual values for commodity balance constraints.
@@ -139,7 +146,7 @@ pub fn perform_dispatch_optimisation<'a>(
     match solution.status() {
         HighsModelStatus::Optimal => Ok(Solution {
             solution: solution.get_solution(),
-            _variables: variables,
+            variables,
             time_slice_info: &model.time_slice_info,
             constraint_keys,
         }),