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;
@@ -75,7 +75,7 @@ pub(crate) use patch_and_validate_simple;
 /// Check whether validation succeeds for simple example with patches
 macro_rules! assert_validate_ok_simple {
     ($file_patches:expr) => {
-        assert!(crate::fixture::patch_and_validate_simple!($file_patches).is_ok())
+        crate::fixture::patch_and_validate_simple!($file_patches).unwrap();
     };
 }
 pub(crate) use assert_validate_ok_simple;
@@ -116,7 +116,7 @@ pub(crate) use patch_and_run_simple;
 /// Check whether the simple example runs successfully after applying file patches
 macro_rules! assert_patched_runs_ok_simple {
     ($file_patches:expr) => {
-        assert!(crate::fixture::patch_and_run_simple!($file_patches).is_ok())
+        crate::fixture::patch_and_run_simple!($file_patches).unwrap();
     };
 }
 pub(crate) use assert_patched_runs_ok_simple;
@@ -408,10 +408,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,
             market_costs,
-            unmet_demand_coefficient: MoneyPerFlow(10000.0),
         }),
         activity,
         unmet_demand,

src/output.rs

@@ -8,9 +8,7 @@ use crate::simulation::investment::appraisal::AppraisalOutput;
 use crate::simulation::optimisation::{FlowMap, Solution};
 use crate::simulation::prices::PriceMap;
 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>,
 }
 
@@ -490,7 +487,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)?;
@@ -1188,7 +1184,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.rs

@@ -859,7 +859,7 @@ fn select_best_assets(
         )?;
 
         // Sort by investment priority and discard non-feasible options
-        sort_and_filter_appraisal_outputs(&mut outputs_for_opts);
+        let num_nonfeasible = sort_and_filter_appraisal_outputs(&mut outputs_for_opts);
 
         // Check if there are any remaining options. If not, we cannot meet demand, so have to bail
         // out.
@@ -872,12 +872,14 @@ fn select_best_assets(
 
             bail!(
                 "No feasible investment options left for \
-                commodity '{}', region '{}', year '{}', agent '{}' after appraisal.\n\
+                commodity '{}', region '{}', year '{}', agent '{}' after appraisal. \
+                {} non-feasible options were not considered.\n\
                 Remaining unmet demand (time_slice : flow):\n{}",
                 &commodity.id,
                 region_id,
                 year,
                 agent.id,
+                num_nonfeasible,
                 remaining_demands.join("\n")
             );
         }

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.market_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!(
@@ -373,13 +359,22 @@ fn compare_asset_fallback(asset1: &Asset, asset2: &Asset) -> Ordering {
 /// with invalid metrics (e.g. `None`) as well as zero capacity. This avoids meaningless or `NaN`
 /// appraisal metrics that could cause the program to panic, so the length of the returned vector
 /// may be less than the input.
-pub fn sort_and_filter_appraisal_outputs(outputs_for_opts: &mut Vec<AppraisalOutput>) {
-    outputs_for_opts.retain(AppraisalOutput::is_valid);
-    outputs_for_opts.sort_by(|output1, output2| match output1.compare_metric(output2) {
+///
+/// # Returns
+///
+/// Returns the number of non-feasible assets which were removed.
+pub fn sort_and_filter_appraisal_outputs(outputs: &mut Vec<AppraisalOutput>) -> usize {
+    let old_len = outputs.len();
+    outputs.retain(AppraisalOutput::is_valid);
+    let num_nonfeasible = old_len - outputs.len();
+
+    outputs.sort_by(|output1, output2| match output1.compare_metric(output2) {
         // If equal, we fall back on comparing asset properties
         Ordering::Equal => compare_asset_fallback(&output1.asset, &output2.asset),
         cmp => cmp,
     });
+
+    num_nonfeasible
 }
 
 /// Counts the number of top appraisal outputs in a sorted slice that are indistinguishable
@@ -405,7 +400,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,10 +569,8 @@ mod tests {
 
     fn objective_coeffs() -> Rc<ObjectiveCoefficients> {
         Rc::new(ObjectiveCoefficients {
-            capacity_coefficient: MoneyPerCapacity(0.0),
             activity_coefficients: IndexMap::new(),
             market_costs: IndexMap::new(),
-            unmet_demand_coefficient: MoneyPerFlow(0.0),
         })
     }
 

src/simulation/investment/appraisal/coefficients.rs

@@ -1,12 +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::PriceMap;
 use crate::simulation::prices::Prices;
 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;
@@ -18,17 +17,22 @@ 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>,
     /// Market costs associated with asset for each time slice
     pub market_costs: IndexMap<TimeSliceID, MoneyPerActivity>,
-    /// Unmet demand coefficient
-    pub unmet_demand_coefficient: MoneyPerFlow,
 }
 
 /// Calculates cost coefficients for a set of assets for a given objective type.
+///
+/// Activity coefficients are 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. These activity coefficients are calculated using shadow prices.
+///
+/// For NPV, "market costs" are calculated in the same way, except thew prices used are market
+/// prices. For LCOX, a slightly different calculation is performed: the sign is inverted (as it
+/// represents a cost) and the primary output (commodity of interest) is excluded.
 pub fn calculate_coefficients_for_assets(
     model: &Model,
     objective_type: &ObjectiveType,
@@ -39,73 +43,22 @@ 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)
-                }
-            };
+            let coefficient = calculate_coefficients_for_asset(
+                asset,
+                objective_type,
+                &model.time_slice_info,
+                prices,
+                year,
+            );
             (asset.clone(), Rc::new(coefficient))
         })
         .collect()
 }
 
-/// Calculates the cost coefficients for LCOX.
-///
-/// 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(
-    asset: &AssetRef,
-    time_slice_info: &TimeSliceInfo,
-    prices: &Prices,
-    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();
-    let mut market_costs = IndexMap::new();
-    for time_slice in time_slice_info.iter_ids() {
-        // 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);
-
-        let coefficient =
-            calculate_asset_costs_for_lcox(asset, operating_cost, time_slice, &prices.shadow);
-        activity_coefficients.insert(time_slice.clone(), coefficient);
-        let market_cost =
-            calculate_asset_costs_for_lcox(asset, operating_cost, time_slice, &prices.market);
-        market_costs.insert(time_slice.clone(), market_cost);
-    }
-
-    // Unmet demand coefficient
-    let unmet_demand_coefficient = value_of_lost_load;
-
-    ObjectiveCoefficients {
-        capacity_coefficient,
-        activity_coefficients,
-        market_costs,
-        unmet_demand_coefficient,
-    }
-}
-
-/// 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(
+/// Calculates cost coefficients for a single asset
+pub fn calculate_coefficients_for_asset(
     asset: &AssetRef,
+    objective_type: &ObjectiveType,
     time_slice_info: &TimeSliceInfo,
     prices: &Prices,
     year: u32,
@@ -123,51 +76,55 @@ pub fn calculate_coefficients_for_npv(
         let operating_cost = asset.get_operating_cost(year, time_slice);
 
         let coefficient =
-            calculate_asset_costs_for_npv(asset, operating_cost, time_slice, &prices.shadow);
+            calculate_asset_revenues(asset, operating_cost, time_slice, &prices.shadow);
         activity_coefficients.insert(
             time_slice.clone(),
             coefficient + EPSILON_ACTIVITY_COEFFICIENT,
         );
-        let market_cost =
-            calculate_asset_costs_for_npv(asset, operating_cost, time_slice, &prices.market);
+
+        let market_cost = match objective_type {
+            ObjectiveType::LevelisedCostOfX => {
+                calculate_asset_costs_for_lcox(asset, operating_cost, time_slice, &prices.market)
+            }
+            ObjectiveType::NetPresentValue => {
+                calculate_asset_revenues(asset, operating_cost, time_slice, &prices.market)
+            }
+        };
         market_costs.insert(time_slice.clone(), market_cost);
     }
 
-    // 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),
-        market_costs,
         activity_coefficients,
-        unmet_demand_coefficient,
+        market_costs,
     }
 }
 
-/// Calculate a single activity coefficient for the LCOX objective for a given time slice.
-fn calculate_asset_costs_for_lcox(
+/// Calculate the revenue from all flows minus operating cost
+fn calculate_asset_revenues(
     asset: &AssetRef,
     operating_cost: MoneyPerActivity,
     time_slice: &TimeSliceID,
     prices: &PriceMap,
 ) -> MoneyPerActivity {
-    // Revenue from flows excluding the primary output
-    let revenue_from_flows = asset.get_revenue_from_flows_excluding_primary(prices, time_slice);
+    // Revenue from flows including the primary output
+    let revenue_from_flows = asset.get_revenue_from_flows(prices, time_slice);
 
-    // The activity coefficient is the operating cost minus the revenue from non-primary flows
-    operating_cost - revenue_from_flows
+    // The activity coefficient is the revenue from flows minus the operating cost (net revenue)
+    revenue_from_flows - operating_cost
 }
 
-/// Calculate a single activity coefficient for the NPV objective for a given time slice.
-fn calculate_asset_costs_for_npv(
+/// Calculate asset costs for LCOX objective.
+///
+/// Excludes revenues from the primary output (commodity of interest).
+fn calculate_asset_costs_for_lcox(
     asset: &AssetRef,
     operating_cost: MoneyPerActivity,
     time_slice: &TimeSliceID,
     prices: &PriceMap,
 ) -> MoneyPerActivity {
-    // Revenue from flows including the primary output
-    let revenue_from_flows = asset.get_revenue_from_flows(prices, time_slice);
+    // Revenue from flows excluding the primary output
+    let revenue_from_flows = asset.get_revenue_from_flows_excluding_primary(prices, time_slice);
 
-    // The activity coefficient is the revenue from flows minus the operating cost (net revenue)
-    revenue_from_flows - operating_cost
+    // The activity coefficient is the operating cost minus the revenue from non-primary flows
+    operating_cost - revenue_from_flows
 }