Add fallback code to calculate metrics for non-dispatched options

docs/release_notes/upcoming.md

@@ -30,10 +30,13 @@ ready to be released, carry out the following steps:
 - Fix parsing and validation of agent search space file ([#1293])
 - Use shadow prices rather than market prices for appraisal optimisations and dispatch runs during
   investment ([#1349])
+- If all investment options fail to dispatch, reappraise assuming maximum activity rather than
+  erroring out ([#1363])
 
 [highs-opts-docs]: ../developer_guide/custom_highs_options.md
 [#1259]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1259
 [#1276]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1276
 [#1281]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1281
 [#1293]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1293
 [#1349]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1349
+[#1363]: https://github.com/EnergySystemsModellingLab/MUSE2/pull/1363

schemas/input/model.yaml

@@ -59,6 +59,11 @@ properties:
       investment cycle. Changing the value of this parameter is potentially dangerous,
       so it requires setting `please_give_me_broken_results` to true.
     default: 1e-12
+  allow_investment_in_non_dispatched_options:
+    type: boolean
+    description: |
+      Whether to allow investment in non-dispatched options.
+    default: true
   highs:
     type: object
     description: |

src/finance.rs

@@ -57,6 +57,11 @@ pub fn profitability_index(
     activity: &IndexMap<TimeSliceID, Activity>,
     activity_surpluses: &IndexMap<TimeSliceID, MoneyPerActivity>,
 ) -> ProfitabilityIndex {
+    assert!(
+        annual_fixed_cost >= MoneyPerCapacity(0.0),
+        "The current NPV calculation does not support negative annual fixed costs"
+    );
+
     // Calculate the annualised fixed costs
     let annualised_fixed_cost = annual_fixed_cost * capacity;
 

src/model/parameters.rs

@@ -100,6 +100,8 @@ pub struct ModelParameters {
     pub mothball_years: u32,
     /// Absolute tolerance when checking if remaining demand is close enough to zero
     pub remaining_demand_absolute_tolerance: Flow,
+    /// Whether to allow investment in non-dispatched options.
+    pub allow_investment_in_non_dispatched_options: bool,
     /// Options for the HiGHS solver.
     ///
     /// For a full list of options, see [the HiGHS documentation].
@@ -125,6 +127,7 @@ impl Default for ModelParameters {
             capacity_margin: Dimensionless(0.2),
             mothball_years: 0,
             remaining_demand_absolute_tolerance: DEFAULT_REMAINING_DEMAND_ABSOLUTE_TOLERANCE,
+            allow_investment_in_non_dispatched_options: true,
             highs: HighsOptions::default(),
         }
     }

src/simulation/investment.rs

@@ -1,6 +1,6 @@
 //! Code for performing agent investment.
 use super::optimisation::{DispatchRun, FlowMap};
-use crate::agent::{Agent, AgentID};
+use crate::agent::{Agent, AgentID, ObjectiveType};
 use crate::asset::{Asset, AssetCapacity, AssetIterator, AssetRef, AssetState};
 use crate::commodity::{Commodity, CommodityID, CommodityMap};
 use crate::model::Model;
@@ -12,16 +12,17 @@ use crate::units::{ActivityPerCapacity, Capacity, Dimensionless, Flow, FlowPerCa
 use anyhow::{Context, Result, bail, ensure};
 use indexmap::IndexMap;
 use itertools::{Itertools, chain};
-use log::debug;
+use log::{debug, warn};
 use std::collections::{HashMap, HashSet};
 use std::fmt::Display;
+use std::rc::Rc;
 use strum::IntoEnumIterator;
 
 pub mod appraisal;
-use appraisal::coefficients::calculate_coefficients_for_assets;
+use appraisal::coefficients::{ObjectiveCoefficients, calculate_coefficients_for_assets};
 use appraisal::{
-    AppraisalOutput, appraise_investment, count_equal_and_best_appraisal_outputs,
-    sort_and_filter_appraisal_outputs,
+    AppraisalOutput, appraise_investment, appraise_investment_assuming_max_activity,
+    count_equal_and_best_appraisal_outputs, sort_and_filter_appraisal_outputs,
 };
 
 /// A map of demand across time slices for a specific market
@@ -772,6 +773,7 @@ fn calculate_investment_limits_for_candidates(
 
 /// Get the best assets for meeting demand for the given commodity
 #[allow(clippy::too_many_arguments)]
+#[allow(clippy::too_many_lines)]
 fn select_best_assets(
     model: &Model,
     mut opt_assets: Vec<AssetRef>,
@@ -851,15 +853,36 @@ fn select_best_assets(
         }
 
         // Save appraisal results
-        writer.write_appraisal_debug_info(
-            year,
-            &format!("{} {} round {}", &commodity.id, &agent.id, round),
-            &outputs_for_opts,
-            &demand,
-        )?;
+        let run_description = format!("{} {} round {}", &commodity.id, &agent.id, round);
+        writer.write_appraisal_debug_info(year, &run_description, &outputs_for_opts, &demand)?;
+
+        if model.parameters.allow_investment_in_non_dispatched_options
+            && !outputs_for_opts
+                .iter()
+                .any(|output| output.metric.is_some())
+        {
+            warn!(
+                "None of the investment options dispatched. Reappraising assuming maximum activity."
+            );
+            outputs_for_opts = reappraise_outputs_assuming_max_activity(
+                &outputs_for_opts,
+                &model.time_slice_info,
+                commodity,
+                *objective_type,
+                &coefficients,
+                &demand,
+            );
+
+            writer.write_appraisal_debug_info(
+                year,
+                &format!("{run_description}; max activity"),
+                &outputs_for_opts,
+                &demand,
+            )?;
+        }
 
         // 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 +895,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")
             );
         }
@@ -926,6 +951,45 @@ fn is_any_remaining_demand(demand: &DemandMap, absolute_tolerance: Flow) -> bool
     demand.values().any(|flow| *flow > absolute_tolerance)
 }
 
+/// Reappraise the specified outputs assuming max activity in each time slice
+fn reappraise_outputs_assuming_max_activity(
+    outputs: &[AppraisalOutput],
+    time_slice_info: &TimeSliceInfo,
+    commodity: &Commodity,
+    objective_type: ObjectiveType,
+    coefficients: &HashMap<AssetRef, Rc<ObjectiveCoefficients>>,
+    demand: &DemandMap,
+) -> Vec<AppraisalOutput> {
+    // Since all assets with the same `group_id` are identical, we only need to appraise one
+    // from each group.
+    let mut seen_groups = HashSet::new();
+
+    let mut new_outputs = Vec::new();
+    for output in outputs {
+        let asset = &output.asset;
+
+        // Skip any assets from groups we've already seen
+        if let Some(group_id) = asset.group_id()
+            && !seen_groups.insert(group_id)
+        {
+            continue;
+        }
+
+        let new_output = appraise_investment_assuming_max_activity(
+            time_slice_info,
+            asset,
+            output.capacity,
+            commodity,
+            objective_type,
+            &coefficients[asset],
+            demand,
+        );
+        new_outputs.push(new_output);
+    }
+
+    new_outputs
+}
+
 /// Update capacity of chosen asset, if needed, and update both asset options and chosen assets
 fn update_assets(
     mut best_asset: AssetRef,

src/simulation/investment/appraisal.rs

@@ -5,8 +5,8 @@ use crate::asset::{Asset, AssetCapacity, AssetRef};
 use crate::commodity::Commodity;
 use crate::finance::{ProfitabilityIndex, lcox, profitability_index};
 use crate::model::Model;
-use crate::time_slice::TimeSliceID;
-use crate::units::{Activity, Capacity, Money, MoneyPerActivity, MoneyPerCapacity};
+use crate::time_slice::{TimeSliceID, TimeSliceInfo};
+use crate::units::{Activity, Capacity, Dimensionless, Flow, Money, MoneyPerActivity};
 use anyhow::Result;
 use costs::annual_fixed_cost;
 use erased_serde::Serialize as ErasedSerialize;
@@ -307,15 +307,9 @@ fn calculate_npv(
         highs::Sense::Maximise,
     )?;
 
-    let annual_fixed_cost = annual_fixed_cost(asset);
-    assert!(
-        annual_fixed_cost >= MoneyPerCapacity(0.0),
-        "The current NPV calculation does not support negative annual fixed costs"
-    );
-
     let profitability_index = profitability_index(
         max_capacity.total_capacity(),
-        annual_fixed_cost,
+        annual_fixed_cost(asset),
         &results.activity,
         &coefficients.market_costs,
     );
@@ -352,6 +346,110 @@ pub fn appraise_investment(
     appraisal_method(model, asset, max_capacity, commodity, coefficients, demand)
 }
 
+/// Computes remaining unmet demand per time slice.
+///
+/// For each time-slice selection at the commodity's balance level, the selection-level residual
+/// (`demand_total - supply_total`, clamped to zero) is divided equally across the time slices in
+/// the selection.
+///
+/// The exact per-time-slice distribution is arbitrary: all downstream consumers sum values back up
+/// to the selection level before using them (e.g. the next round's demand constraint, and
+/// `is_any_remaining_demand`), so only the selection-level total needs to be correct.
+fn compute_unmet_demand(
+    demand: &DemandMap,
+    activity: &IndexMap<TimeSliceID, Activity>,
+    commodity: &Commodity,
+    asset: &Asset,
+    time_slice_info: &TimeSliceInfo,
+) -> DemandMap {
+    let mut unmet_demand = DemandMap::new();
+    let flow_coeff = asset.get_flow(&commodity.id).unwrap().coeff;
+    for ts_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level) {
+        let time_slices: Vec<_> = ts_selection.iter(time_slice_info).collect();
+        let demand_for_selection: Flow = time_slices.iter().map(|(ts, _)| demand[*ts]).sum();
+        let supply_for_selection: Flow = time_slices
+            .iter()
+            .map(|(ts, _)| activity[*ts] * flow_coeff)
+            .sum();
+
+        #[allow(clippy::cast_precision_loss)]
+        let unmet_per_slice = (demand_for_selection - supply_for_selection).max(Flow(0.0))
+            / Dimensionless(time_slices.len() as f64);
+        for (time_slice, _) in &time_slices {
+            unmet_demand.insert((*time_slice).clone(), unmet_per_slice);
+        }
+    }
+    unmet_demand
+}
+
+/// Get the maximum allowed activity per time slice for an asset
+fn iter_max_activity_per_time_slice(
+    asset: &Asset,
+    capacity: Capacity,
+    time_slice_info: &TimeSliceInfo,
+) -> impl Iterator<Item = (TimeSliceID, Activity)> {
+    time_slice_info.iter_ids().cloned().map(move |time_slice| {
+        let max_act = *asset.get_activity_per_capacity_limits(&time_slice).end();
+        (time_slice, max_act * capacity)
+    })
+}
+
+/// Appraise an investment assuming maximum activity in every time slice
+pub fn appraise_investment_assuming_max_activity(
+    time_slice_info: &TimeSliceInfo,
+    asset: &AssetRef,
+    capacity: AssetCapacity,
+    commodity: &Commodity,
+    objective_type: ObjectiveType,
+    coefficients: &Rc<ObjectiveCoefficients>,
+    demand: &DemandMap,
+) -> AppraisalOutput {
+    let activity =
+        iter_max_activity_per_time_slice(asset, capacity.total_capacity(), time_slice_info)
+            .collect();
+    let unmet_demand = compute_unmet_demand(demand, &activity, commodity, asset, time_slice_info);
+
+    let results = ResultsMap {
+        capacity,
+        activity,
+        unmet_demand,
+    };
+    match objective_type {
+        ObjectiveType::LevelisedCostOfX => {
+            let cost_index = lcox(
+                capacity.total_capacity(),
+                coefficients.capacity_coefficient,
+                &results.activity,
+                &coefficients.market_costs,
+            );
+
+            AppraisalOutput::new(
+                asset.clone(),
+                capacity,
+                results,
+                cost_index.map(LCOXMetric::new),
+                coefficients.clone(),
+            )
+        }
+        ObjectiveType::NetPresentValue => {
+            let profitability_index = profitability_index(
+                capacity.total_capacity(),
+                annual_fixed_cost(asset),
+                &results.activity,
+                &coefficients.market_costs,
+            );
+
+            AppraisalOutput::new(
+                asset.clone(),
+                capacity,
+                results,
+                Some(NPVMetric::new(profitability_index)),
+                coefficients.clone(),
+            )
+        }
+    }
+}
+
 /// Compare assets as a fallback if metrics are equal.
 ///
 /// Commissioned assets are ordered before uncommissioned and newer before older.
@@ -373,13 +471,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 +512,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, MoneyPerCapacity, MoneyPerFlow};
     use float_cmp::assert_approx_eq;
     use rstest::rstest;
     use std::rc::Rc;