Add commodity prices

Closes #589.

Author: alexdewar

src/simulation.rs

@@ -52,7 +52,7 @@ 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);
+        let prices = CommodityPrices::from_model_and_solution(&model, &solution);
 
         // Write result of dispatch optimisation to file
         writer.write_debug_info(year, &solution)?;

src/simulation/prices.rs

@@ -1,12 +1,11 @@
 //! Code for updating the simulation state.
 use super::optimisation::Solution;
-use crate::asset::AssetPool;
 use crate::commodity::CommodityID;
 use crate::model::Model;
 use crate::region::RegionID;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
 use log::warn;
-use std::collections::{BTreeMap, HashSet};
+use std::collections::{BTreeMap, HashMap, HashSet};
 
 /// A map relating commodity ID + region + time slice to current price (endogenous)
 #[derive(Default)]
@@ -16,9 +15,9 @@ impl CommodityPrices {
     /// Calculate commodity prices based on the result of the dispatch optimisation.
     ///
     /// Missing prices will be calculated directly from the input data
-    pub fn from_model_and_solution(model: &Model, solution: &Solution, assets: &AssetPool) -> Self {
+    pub fn from_model_and_solution(model: &Model, solution: &Solution) -> Self {
         let mut prices = CommodityPrices::default();
-        let commodity_regions_updated = prices.add_from_solution(solution, assets);
+        let commodity_regions_updated = prices.add_from_solution(solution);
 
         // Find commodity/region combinations not updated in last step
         let mut remaining_commodity_regions = HashSet::new();
@@ -44,20 +43,43 @@ impl CommodityPrices {
     /// # Arguments
     ///
     /// * `solution` - The solution to the dispatch optimisation
-    /// * `assets` - The asset pool
     ///
     /// # Returns
     ///
     /// The set of commodities for which prices were added.
-    fn add_from_solution(
-        &mut self,
-        _solution: &Solution,
-        _assets: &AssetPool,
-    ) -> HashSet<(CommodityID, RegionID)> {
-        // **TODO:** Calculate commodity prices here:
-        //  https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/589
+    fn add_from_solution(&mut self, solution: &Solution) -> HashSet<(CommodityID, RegionID)> {
+        let mut commodity_regions_updated = HashSet::new();
+
+        // Calculate highest activity dual for each commodity/region/timeslice
+        let mut highest_duals = HashMap::new();
+        for (asset, time_slice, dual) in solution.iter_activity_duals() {
+            // Iterate over all output flows
+            for flow in asset.iter_flows().filter(|flow| flow.coeff > 0.0) {
+                // Update the highest dual for this commodity/timeslice
+                highest_duals
+                    .entry((
+                        flow.commodity.id.clone(),
+                        asset.region_id.clone(),
+                        time_slice.clone(),
+                    ))
+                    .and_modify(|current_dual| {
+                        if dual > *current_dual {
+                            *current_dual = dual;
+                        }
+                    })
+                    .or_insert(dual);
+            }
+        }
+
+        // Add the highest capacity dual for each commodity/timeslice to each commodity balance dual
+        for (commodity_id, region_id, time_slice, dual) in solution.iter_commodity_balance_duals() {
+            let key = (commodity_id.clone(), region_id.clone(), time_slice.clone());
+            let price = dual + highest_duals.get(&key).unwrap_or(&0.0);
+            self.insert(commodity_id, region_id, time_slice, price);
+            commodity_regions_updated.insert((commodity_id.clone(), region_id.clone()));
+        }
 
-        HashSet::new()
+        commodity_regions_updated
     }
 
     /// Add prices for any commodity not updated by the dispatch step.