Change LCOX to use same optimisation as NPV

src/fixture.rs

@@ -21,7 +21,7 @@ use crate::simulation::investment::appraisal::{
 use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceLevel};
 use crate::units::{
     Activity, ActivityPerCapacity, Capacity, Dimensionless, Flow, MoneyPerActivity,
-    MoneyPerCapacity, MoneyPerCapacityPerYear, MoneyPerFlow, Year,
+    MoneyPerCapacity, MoneyPerCapacityPerYear, Year,
 };
 use anyhow::Result;
 use indexmap::indexmap;
@@ -407,9 +407,8 @@ pub fn appraisal_output(asset: Asset, time_slice: TimeSliceID) -> AppraisalOutpu
         asset: AssetRef::from(asset),
         capacity: AssetCapacity::Continuous(Capacity(42.0)),
         coefficients: Rc::new(ObjectiveCoefficients {
-            capacity_coefficient: MoneyPerCapacity(2.14),
             activity_coefficients,
-            unmet_demand_coefficient: MoneyPerFlow(10000.0),
+            lcox_costs: IndexMap::new(),
         }),
         activity,
         unmet_demand,

src/output.rs

@@ -8,9 +8,7 @@ use crate::simulation::CommodityPrices;
 use crate::simulation::investment::appraisal::AppraisalOutput;
 use crate::simulation::optimisation::{FlowMap, Solution};
 use crate::time_slice::TimeSliceID;
-use crate::units::{
-    Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow,
-};
+use crate::units::{Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerFlow};
 use anyhow::{Context, Result, ensure};
 use csv;
 use indexmap::IndexMap;
@@ -261,7 +259,6 @@ struct AppraisalResultsRow {
     process_id: ProcessID,
     region_id: RegionID,
     capacity: Capacity,
-    capacity_coefficient: MoneyPerCapacity,
     metric: Option<f64>,
 }
 
@@ -495,7 +492,6 @@ impl DebugDataWriter {
                 process_id: result.asset.process_id().clone(),
                 region_id: result.asset.region_id().clone(),
                 capacity: result.capacity.total_capacity(),
-                capacity_coefficient: result.coefficients.capacity_coefficient,
                 metric: result.metric.as_ref().map(|m| m.value()),
             };
             self.appraisal_results_writer.serialize(row)?;
@@ -1195,7 +1191,6 @@ mod tests {
             process_id: asset.process_id().clone(),
             region_id: asset.region_id().clone(),
             capacity: Capacity(42.0),
-            capacity_coefficient: MoneyPerCapacity(2.14),
             metric: Some(4.14),
         };
         let records: Vec<AppraisalResultsRow> =

src/simulation/investment/appraisal.rs

@@ -258,26 +258,19 @@ fn calculate_lcox(
     coefficients: &Rc<ObjectiveCoefficients>,
     demand: &DemandMap,
 ) -> Result<AppraisalOutput> {
-    let results = perform_optimisation(
-        model,
-        asset,
-        max_capacity,
-        commodity,
-        coefficients,
-        demand,
-        highs::Sense::Minimise,
-    )?;
+    let results =
+        perform_optimisation(model, asset, max_capacity, commodity, coefficients, demand)?;
 
     let cost_index = lcox(
-        results.capacity.total_capacity(),
-        coefficients.capacity_coefficient,
+        max_capacity.total_capacity(),
+        annual_fixed_cost(asset),
         &results.activity,
-        &coefficients.activity_coefficients,
+        &coefficients.lcox_costs,
     );
 
     Ok(AppraisalOutput::new(
         asset.clone(),
-        results.capacity,
+        max_capacity,
         results,
         cost_index.map(LCOXMetric::new),
         coefficients.clone(),
@@ -297,15 +290,8 @@ fn calculate_npv(
     coefficients: &Rc<ObjectiveCoefficients>,
     demand: &DemandMap,
 ) -> Result<AppraisalOutput> {
-    let results = perform_optimisation(
-        model,
-        asset,
-        max_capacity,
-        commodity,
-        coefficients,
-        demand,
-        highs::Sense::Maximise,
-    )?;
+    let results =
+        perform_optimisation(model, asset, max_capacity, commodity, coefficients, demand)?;
 
     let annual_fixed_cost = annual_fixed_cost(asset);
     assert!(
@@ -314,7 +300,7 @@ fn calculate_npv(
     );
 
     let profitability_index = profitability_index(
-        results.capacity.total_capacity(),
+        max_capacity.total_capacity(),
         annual_fixed_cost,
         &results.activity,
         &coefficients.activity_coefficients,
@@ -405,7 +391,7 @@ mod tests {
     use crate::fixture::{agent_id, asset, process, region_id};
     use crate::process::Process;
     use crate::region::RegionID;
-    use crate::units::{Money, MoneyPerActivity, MoneyPerFlow};
+    use crate::units::{Money, MoneyPerActivity};
     use float_cmp::assert_approx_eq;
     use rstest::rstest;
     use std::rc::Rc;
@@ -574,9 +560,8 @@ mod tests {
 
     fn objective_coeffs() -> Rc<ObjectiveCoefficients> {
         Rc::new(ObjectiveCoefficients {
-            capacity_coefficient: MoneyPerCapacity(0.0),
             activity_coefficients: IndexMap::new(),
-            unmet_demand_coefficient: MoneyPerFlow(0.0),
+            lcox_costs: IndexMap::new(),
         })
     }
 

src/simulation/investment/appraisal/coefficients.rs

@@ -1,11 +1,11 @@
 //! Calculation of cost coefficients for investment tools.
-use super::costs::annual_fixed_cost;
 use crate::agent::ObjectiveType;
 use crate::asset::AssetRef;
 use crate::model::Model;
+
 use crate::simulation::CommodityPrices;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
-use crate::units::{MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow};
+use crate::units::MoneyPerActivity;
 use indexmap::IndexMap;
 use std::collections::HashMap;
 use std::rc::Rc;
@@ -17,12 +17,10 @@ use std::rc::Rc;
 /// coefficients used in the appraisal optimisation, together with the unmet-demand penalty.
 #[derive(Clone)]
 pub struct ObjectiveCoefficients {
-    /// Cost per unit of capacity
-    pub capacity_coefficient: MoneyPerCapacity,
     /// Cost per unit of activity in each time slice
     pub activity_coefficients: IndexMap<TimeSliceID, MoneyPerActivity>,
-    /// Unmet demand coefficient
-    pub unmet_demand_coefficient: MoneyPerFlow,
+    /// Costs for LCOX
+    pub lcox_costs: IndexMap<TimeSliceID, MoneyPerActivity>,
 }
 
 /// Calculates cost coefficients for a set of assets for a given objective type.
@@ -36,93 +34,64 @@ pub fn calculate_coefficients_for_assets(
     assets
         .iter()
         .map(|asset| {
-            let coefficient = match objective_type {
-                ObjectiveType::LevelisedCostOfX => calculate_coefficients_for_lcox(
-                    asset,
-                    &model.time_slice_info,
-                    prices,
-                    model.parameters.value_of_lost_load,
-                    year,
-                ),
-                ObjectiveType::NetPresentValue => {
-                    calculate_coefficients_for_npv(asset, &model.time_slice_info, prices, year)
-                }
-            };
-            (asset.clone(), Rc::new(coefficient))
+            let mut coeff =
+                calculate_coefficients_for_asset(asset, &model.time_slice_info, prices, year);
+
+            // For LCOX, we also store per-time slice costs, which are used to calculate the cost
+            // index
+            if objective_type == &ObjectiveType::LevelisedCostOfX {
+                coeff.lcox_costs =
+                    calculate_lcox_costs(asset, &model.time_slice_info, prices, year);
+            }
+
+            (asset.clone(), Rc::new(coeff))
         })
         .collect()
 }
 
-/// Calculates the cost coefficients for LCOX.
+/// Calculates the cost coefficients for a single asset.
 ///
-/// For LCOX the activity coefficient is calculated as operating cost minus revenue from
-/// non-primary flows. The unmet demand coefficient is set from the model parameter
-/// `value_of_lost_load`.
-pub fn calculate_coefficients_for_lcox(
+/// The activity coefficient is revenue minus operating cost; a small positive epsilon is added to
+/// activity coefficients so that assets with near-zero net value still appear in dispatch. Capacity
+/// costs and unmet-demand penalties are set to zero.
+pub fn calculate_coefficients_for_asset(
     asset: &AssetRef,
     time_slice_info: &TimeSliceInfo,
     prices: &CommodityPrices,
-    value_of_lost_load: MoneyPerFlow,
     year: u32,
 ) -> ObjectiveCoefficients {
-    // Capacity coefficient
-    let capacity_coefficient = annual_fixed_cost(asset);
-
     // Activity coefficients
     let mut activity_coefficients = IndexMap::new();
     for time_slice in time_slice_info.iter_ids() {
-        let coefficient = calculate_activity_coefficient_for_lcox(asset, time_slice, prices, year);
+        let coefficient = calculate_activity_coefficient(asset, time_slice, prices, year);
         activity_coefficients.insert(time_slice.clone(), coefficient);
     }
 
-    // Unmet demand coefficient
-    let unmet_demand_coefficient = value_of_lost_load;
-
     ObjectiveCoefficients {
-        capacity_coefficient,
         activity_coefficients,
-        unmet_demand_coefficient,
+        lcox_costs: IndexMap::new(),
     }
 }
 
-/// Calculates the cost coefficients for NPV.
-///
-/// For NPV the activity coefficient is revenue (including primary output) minus operating
-/// cost; a small positive epsilon is added to activity coefficients so that assets with
-/// near-zero net value still appear in dispatch. Capacity costs and unmet-demand penalties
-/// are set to zero for the NPV objective.
-pub fn calculate_coefficients_for_npv(
+/// Calculate costs for LCOX for all time slices
+fn calculate_lcox_costs(
     asset: &AssetRef,
     time_slice_info: &TimeSliceInfo,
     prices: &CommodityPrices,
     year: u32,
-) -> ObjectiveCoefficients {
-    // Small constant added to each activity coefficient to ensure break-even/slightly negative
-    // assets are still dispatched
-    const EPSILON_ACTIVITY_COEFFICIENT: MoneyPerActivity = MoneyPerActivity(f64::EPSILON * 100.0);
-
-    // Activity coefficients
-    let mut activity_coefficients = IndexMap::new();
-    for time_slice in time_slice_info.iter_ids() {
-        let coefficient = calculate_activity_coefficient_for_npv(asset, time_slice, prices, year);
-        activity_coefficients.insert(
-            time_slice.clone(),
-            coefficient + EPSILON_ACTIVITY_COEFFICIENT,
-        );
-    }
-
-    // Unmet demand coefficient (we don't apply a cost to unmet demand, so we set this to zero)
-    let unmet_demand_coefficient = MoneyPerFlow(0.0);
-
-    ObjectiveCoefficients {
-        capacity_coefficient: MoneyPerCapacity(0.0),
-        activity_coefficients,
-        unmet_demand_coefficient,
-    }
+) -> IndexMap<TimeSliceID, MoneyPerActivity> {
+    time_slice_info
+        .iter_ids()
+        .cloned()
+        .map(|time_slice| {
+            let coefficient = calculate_lcox_cost_for_time_slice(asset, &time_slice, prices, year);
+            (time_slice, coefficient)
+        })
+        .collect()
 }
 
-/// Calculate a single activity coefficient for the LCOX objective for a given time slice.
-fn calculate_activity_coefficient_for_lcox(
+/// Calculate cost per unit activity excluding the primary commodity for a given time slice.
+fn calculate_lcox_cost_for_time_slice(
     asset: &AssetRef,
     time_slice: &TimeSliceID,
     prices: &CommodityPrices,
@@ -139,13 +108,17 @@ fn calculate_activity_coefficient_for_lcox(
     operating_cost - revenue_from_flows
 }
 
-/// Calculate a single activity coefficient for the NPV objective for a given time slice.
-fn calculate_activity_coefficient_for_npv(
+/// Calculate a single activity coefficient for a given time slice
+fn calculate_activity_coefficient(
     asset: &AssetRef,
     time_slice: &TimeSliceID,
     prices: &CommodityPrices,
     year: u32,
 ) -> MoneyPerActivity {
+    // Small constant added to each activity coefficient to ensure break-even/slightly negative
+    // assets are still dispatched
+    const EPSILON_ACTIVITY_COEFFICIENT: MoneyPerActivity = MoneyPerActivity(f64::EPSILON * 100.0);
+
     // Get the operating cost of the asset. This includes the variable operating cost, levies and
     // flow costs, but excludes costs/revenues from commodity consumption/production.
     let operating_cost = asset.get_operating_cost(year, time_slice);
@@ -154,5 +127,5 @@ fn calculate_activity_coefficient_for_npv(
     let revenue_from_flows = asset.get_revenue_from_flows(prices, time_slice);
 
     // The activity coefficient is the revenue from flows minus the operating cost (net revenue)
-    revenue_from_flows - operating_cost
+    revenue_from_flows - operating_cost + EPSILON_ACTIVITY_COEFFICIENT
 }