Merge pull request #882 from EnergySystemsModellingLab/exclude-coi-price-for-lcox

Fix: Exclude commodity of interest from reduced costs calculation for LCOX assets

Author: alexdewar

Cargo.lock

@@ -39,7 +39,6 @@ dirs = "6.0.0"
 edit = "0.1.5"
 
 [dev-dependencies]
-current_dir = "0.1.2"
 map-macro = "0.3.0"
 rstest = {version = "0.26.1", default-features = false, features = ["crate-name"]}
 

Cargo.toml

@@ -100,3 +100,14 @@ pub enum ObjectiveType {
     #[string = "npv"]
     NetPresentValue,
 }
+
+impl ObjectiveType {
+    /// Whether to exclude the price of the primary output commodity from reduced cost calculation
+    pub fn exclude_primary_output_price_from_reduced_costs(&self) -> bool {
+        // Deliberately written as a `match` block, in case we add more objective types in future
+        match self {
+            ObjectiveType::LevelisedCostOfX => true,
+            ObjectiveType::NetPresentValue => false,
+        }
+    }
+}

src/agent.rs

@@ -1,12 +1,11 @@
 //! Assets are instances of a process which are owned and invested in by agents.
-use crate::agent::AgentID;
+use crate::agent::{AgentID, AgentMap};
 use crate::commodity::CommodityID;
 use crate::process::{Process, ProcessFlow, ProcessID, ProcessParameter};
 use crate::region::RegionID;
+use crate::simulation::CommodityPrices;
 use crate::time_slice::TimeSliceID;
-use crate::units::{
-    Activity, ActivityPerCapacity, Capacity, Dimensionless, MoneyPerActivity, MoneyPerFlow,
-};
+use crate::units::{Activity, ActivityPerCapacity, Capacity, Dimensionless, MoneyPerActivity};
 use anyhow::{Context, Result, ensure};
 use indexmap::IndexMap;
 use itertools::{Itertools, chain};
@@ -272,23 +271,67 @@ impl Asset {
         self.process_parameter.variable_operating_cost + flows_cost
     }
 
-    /// Get the cost of input flows using the commodity prices in `input_prices`
+    /// Get the total revenue from all flows for this asset, accounting for the parent agent's
+    /// objective.
+    ///
+    /// We need to account for the agent's objective when calculating reduced costs, because if it
+    /// is LCOX then we should exclude the primary output from the calculation.
+    ///
+    /// If a price is missing from `prices`, then it is assumed to be zero.
+    ///
+    /// # Panics
+    ///
+    /// Panics if this asset has no parent agent (i.e. it's a candidate).
+    pub fn get_revenue_from_flows_for_objective(
+        &self,
+        agents: &AgentMap,
+        prices: &CommodityPrices,
+        year: u32,
+        time_slice: &TimeSliceID,
+    ) -> MoneyPerActivity {
+        let exclude_commodity = self.primary_output().and_then(|flow| {
+            let agent = &agents[self.agent_id().unwrap()];
+            let exclude_coi =
+                agent.objectives[&year].exclude_primary_output_price_from_reduced_costs();
+            exclude_coi.then_some(&flow.commodity.id)
+        });
+
+        self.get_revenue_from_flows_with_filter(prices, time_slice, |flow| {
+            exclude_commodity.is_none_or(|commodity_id| commodity_id != &flow.commodity.id)
+        })
+    }
+
+    /// Get the cost of input flows using the commodity prices in `input_prices`.
+    ///
+    /// If a price is missing, there is assumed to be no cost.
     pub fn get_input_cost_from_prices(
         &self,
-        input_prices: &HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>,
+        input_prices: &CommodityPrices,
         time_slice: &TimeSliceID,
     ) -> MoneyPerActivity {
+        -self.get_revenue_from_flows_with_filter(input_prices, time_slice, ProcessFlow::is_input)
+    }
+
+    /// Get the total revenue from a subset of flows.
+    ///
+    /// Takes a function as an argument to filter the flows. If a price is missing, it is assumed to
+    /// be zero.
+    fn get_revenue_from_flows_with_filter<F>(
+        &self,
+        prices: &CommodityPrices,
+        time_slice: &TimeSliceID,
+        mut filter_for_flows: F,
+    ) -> MoneyPerActivity
+    where
+        F: FnMut(&ProcessFlow) -> bool,
+    {
         self.iter_flows()
-            .filter_map(|flow| {
-                if !flow.is_input() {
-                    return None;
-                }
-                let price = *input_prices.get(&(
-                    flow.commodity.id.clone(),
-                    self.region_id.clone(),
-                    time_slice.clone(),
-                ))?;
-                Some(-flow.coeff * price)
+            .filter(|flow| filter_for_flows(flow))
+            .map(|flow| {
+                flow.coeff
+                    * prices
+                        .get(&flow.commodity.id, self.region_id(), time_slice)
+                        .unwrap_or_default()
             })
             .sum()
     }
@@ -905,11 +948,8 @@ mod tests {
         let asset = Asset::new_candidate(process, region_id.clone(), Capacity(1.0), 2020).unwrap();
 
         // Set input prices
-        let mut input_prices = HashMap::new();
-        input_prices.insert(
-            (commodity_id.clone(), region_id.clone(), time_slice.clone()),
-            MoneyPerFlow(3.0),
-        );
+        let mut input_prices = CommodityPrices::default();
+        input_prices.insert(&commodity_id, &region_id, &time_slice, MoneyPerFlow(3.0));
 
         // Call function
         let cost = asset.get_input_cost_from_prices(&input_prices, &time_slice);

src/asset.rs

@@ -9,7 +9,7 @@ use crate::output::DataWriter;
 use crate::region::RegionID;
 use crate::simulation::CommodityPrices;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
-use crate::units::{Capacity, Dimensionless, Flow, FlowPerCapacity, MoneyPerFlow};
+use crate::units::{Capacity, Dimensionless, Flow, FlowPerCapacity};
 use anyhow::{Result, ensure};
 use indexmap::IndexMap;
 use itertools::{chain, iproduct};
@@ -68,11 +68,7 @@ pub fn perform_agent_investment(
             // performed will, by definition, not have any producers. For these, we provide prices
             // from the previous dispatch run otherwise they will appear to be free to the model.
             for time_slice in model.time_slice_info.iter_ids() {
-                external_prices.remove(&(
-                    commodity_id.clone(),
-                    region_id.clone(),
-                    time_slice.clone(),
-                ));
+                external_prices.remove(commodity_id, region_id, time_slice);
             }
 
             // List of assets selected/retained for this region/commodity
@@ -347,14 +343,11 @@ fn get_prices_for_commodities(
     time_slice_info: &TimeSliceInfo,
     region_id: &RegionID,
     commodities: &[CommodityID],
-) -> HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow> {
+) -> CommodityPrices {
     iproduct!(commodities.iter(), time_slice_info.iter_ids())
         .map(|(commodity_id, time_slice)| {
             let price = prices.get(commodity_id, region_id, time_slice).unwrap();
-            (
-                (commodity_id.clone(), region_id.clone(), time_slice.clone()),
-                price,
-            )
+            (commodity_id, region_id, time_slice, price)
         })
         .collect()
 }

src/simulation/investment.rs

@@ -6,14 +6,15 @@ use crate::commodity::CommodityID;
 use crate::model::Model;
 use crate::output::DataWriter;
 use crate::region::RegionID;
+use crate::simulation::CommodityPrices;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
 use crate::units::{Activity, Flow, Money, MoneyPerActivity, MoneyPerFlow, UnitType};
 use anyhow::{Result, anyhow, ensure};
 use highs::{HighsModelStatus, RowProblem as Problem, Sense};
 use indexmap::IndexMap;
 use itertools::{chain, iproduct};
 use log::debug;
-use std::collections::{HashMap, HashSet};
+use std::collections::HashSet;
 use std::ops::Range;
 
 mod constraints;
@@ -180,10 +181,7 @@ pub fn solve_optimal(model: highs::Model) -> Result<highs::SolvedModel> {
 ///
 /// Input prices should only be provided for commodities for which there will be no commodity
 /// balance constraint.
-fn check_input_prices(
-    input_prices: &HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>,
-    commodities: &[CommodityID],
-) {
+fn check_input_prices(input_prices: &CommodityPrices, commodities: &[CommodityID]) {
     let commodities_set: HashSet<_> = commodities.iter().collect();
     let has_prices_for_commodity_subset = input_prices
         .keys()
@@ -204,7 +202,7 @@ pub struct DispatchRun<'model, 'run> {
     existing_assets: &'run [AssetRef],
     candidate_assets: &'run [AssetRef],
     commodities: &'run [CommodityID],
-    input_prices: Option<&'run HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>>,
+    input_prices: Option<&'run CommodityPrices>,
     year: u32,
 }
 
@@ -240,10 +238,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> {
     }
 
     /// Explicitly provide prices for certain input commodities
-    pub fn with_input_prices(
-        self,
-        input_prices: &'run HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>,
-    ) -> Self {
+    pub fn with_input_prices(self, input_prices: &'run CommodityPrices) -> Self {
         Self {
             input_prices: Some(input_prices),
             ..self
@@ -346,7 +341,7 @@ fn add_variables(
     problem: &mut Problem,
     variables: &mut VariableMap,
     time_slice_info: &TimeSliceInfo,
-    input_prices: Option<&HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>>,
+    input_prices: Option<&CommodityPrices>,
     assets: &[AssetRef],
     year: u32,
 ) -> Range<usize> {
@@ -384,7 +379,7 @@ fn calculate_cost_coefficient(
     asset: &Asset,
     year: u32,
     time_slice: &TimeSliceID,
-    input_prices: Option<&HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>>,
+    input_prices: Option<&CommodityPrices>,
 ) -> MoneyPerActivity {
     let opex = asset.get_operating_cost(year, time_slice);
     let input_cost = input_prices

src/simulation/optimisation.rs

@@ -9,7 +9,7 @@ use crate::time_slice::{TimeSliceID, TimeSliceInfo};
 use crate::units::{Dimensionless, MoneyPerActivity, MoneyPerFlow, Year};
 use indexmap::IndexMap;
 use itertools::iproduct;
-use std::collections::{BTreeMap, HashMap};
+use std::collections::{BTreeMap, HashMap, btree_map};
 
 /// A map of reduced costs for different assets in different time slices
 ///
@@ -146,7 +146,7 @@ pub fn calculate_prices_and_reduced_costs(
     // Add new reduced costs, using old values if not provided
     reduced_costs.extend(reduced_costs_for_candidates);
     reduced_costs.extend(reduced_costs_for_existing(
-        &model.time_slice_info,
+        model,
         existing_assets,
         &prices,
         year,
@@ -256,6 +256,22 @@ impl CommodityPrices {
             .copied()
     }
 
+    /// Iterate over the price map's keys
+    pub fn keys(&self) -> btree_map::Keys<'_, (CommodityID, RegionID, TimeSliceID), MoneyPerFlow> {
+        self.0.keys()
+    }
+
+    /// Remove the specified entry from the map
+    pub fn remove(
+        &mut self,
+        commodity_id: &CommodityID,
+        region_id: &RegionID,
+        time_slice: &TimeSliceID,
+    ) -> Option<MoneyPerFlow> {
+        self.0
+            .remove(&(commodity_id.clone(), region_id.clone(), time_slice.clone()))
+    }
+
     /// Calculate time slice-weighted average prices for each commodity-region pair
     ///
     /// This method aggregates prices across time slices by weighting each price
@@ -419,22 +435,15 @@ fn get_scarcity_adjustment(
 
 /// Calculate reduced costs for existing assets
 fn reduced_costs_for_existing<'a>(
-    time_slice_info: &'a TimeSliceInfo,
+    model: &'a Model,
     assets: &'a [AssetRef],
     prices: &'a CommodityPrices,
     year: u32,
 ) -> impl Iterator<Item = ((AssetRef, TimeSliceID), MoneyPerActivity)> + 'a {
-    iproduct!(assets, time_slice_info.iter_ids()).map(move |(asset, time_slice)| {
+    iproduct!(assets, model.time_slice_info.iter_ids()).map(move |(asset, time_slice)| {
         let operating_cost = asset.get_operating_cost(year, time_slice);
-        let revenue_from_flows = asset
-            .iter_flows()
-            .map(|flow| {
-                flow.coeff
-                    * prices
-                        .get(&flow.commodity.id, asset.region_id(), time_slice)
-                        .unwrap()
-            })
-            .sum();
+        let revenue_from_flows =
+            asset.get_revenue_from_flows_for_objective(&model.agents, prices, year, time_slice);
         let reduced_cost = operating_cost - revenue_from_flows;
 
         ((asset.clone(), time_slice.clone()), reduced_cost)