Merge pull request #1137 from EnergySystemsModellingLab/add-unmet-demand-threshold

add remaining_demand_absolute_tolerance param

Author: Aurashk

schemas/input/model.yaml

@@ -50,5 +50,12 @@ properties:
     type: boolean
     description: |
       Allows other options that are known to be broken to be used. Please don't ever enable this.
+  remaining_demand_absolute_tolerance:
+    type: number
+    description: |
+      Absolute tolerance when checking if remaining demand is close enough to zero in the
+      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
 
 required: [milestone_years]

src/model/parameters.rs

@@ -7,7 +7,7 @@ use crate::asset::check_capacity_valid_for_asset;
 use crate::input::{
     deserialise_proportion_nonzero, input_err_msg, is_sorted_and_unique, read_toml,
 };
-use crate::units::{Capacity, Dimensionless, MoneyPerFlow};
+use crate::units::{Capacity, Dimensionless, Flow, MoneyPerFlow};
 use anyhow::{Context, Result, ensure};
 use log::warn;
 use serde::Deserialize;
@@ -76,6 +76,7 @@ define_unit_param_default!(default_candidate_asset_capacity, Capacity, 0.0001);
 define_unit_param_default!(default_capacity_limit_factor, Dimensionless, 0.1);
 define_unit_param_default!(default_value_of_lost_load, MoneyPerFlow, 1e9);
 define_unit_param_default!(default_price_tolerance, Dimensionless, 1e-6);
+define_unit_param_default!(default_remaining_demand_absolute_tolerance, Flow, 1e-12);
 define_param_default!(default_max_ironing_out_iterations, u32, 10);
 define_param_default!(default_capacity_margin, f64, 0.2);
 define_param_default!(default_mothball_years, u32, 0);
@@ -123,6 +124,9 @@ pub struct ModelParameters {
     /// Number of years an asset can remain unused before being decommissioned
     #[serde(default = "default_mothball_years")]
     pub mothball_years: u32,
+    /// Absolute tolerance when checking if remaining demand is close enough to zero
+    #[serde(default = "default_remaining_demand_absolute_tolerance")]
+    pub remaining_demand_absolute_tolerance: Flow,
 }
 
 /// Check that the `milestone_years` parameter is valid
@@ -164,6 +168,29 @@ fn check_price_tolerance(value: Dimensionless) -> Result<()> {
     Ok(())
 }
 
+fn check_remaining_demand_absolute_tolerance(
+    allow_broken_options: bool,
+    value: Flow,
+) -> Result<()> {
+    ensure!(
+        value.is_finite() && value >= Flow(0.0),
+        "remaining_demand_absolute_tolerance must be a finite number greater than or equal to zero"
+    );
+
+    let default_value = default_remaining_demand_absolute_tolerance();
+    if !allow_broken_options {
+        ensure!(
+            value == default_value,
+            "Setting a remaining_demand_absolute_tolerance different from the default value of {:e} \
+             is potentially dangerous, set please_give_me_broken_results to true \
+             if you want to allow this.",
+            default_value.0
+        );
+    }
+
+    Ok(())
+}
+
 /// Check that the `capacity_margin` parameter is valid
 fn check_capacity_margin(value: f64) -> Result<()> {
     ensure!(
@@ -229,6 +256,12 @@ impl ModelParameters {
         // capacity_margin
         check_capacity_margin(self.capacity_margin)?;
 
+        // remaining_demand_absolute_tolerance
+        check_remaining_demand_absolute_tolerance(
+            self.allow_broken_options,
+            self.remaining_demand_absolute_tolerance,
+        )?;
+
         Ok(())
     }
 }
@@ -356,6 +389,57 @@ mod tests {
         );
     }
 
+    #[rstest]
+    #[case(true, 0.0, true)] // Valid minimum value broken options allowed
+    #[case(true, 1e-10, true)] // Valid value with broken options allowed
+    #[case(true, 1e-15, true)] // Valid value with broken options allowed
+    #[case(false, 1e-12, true)] // Valid value same as default, no broken options needed
+    #[case(true, 1.0, true)] // Valid larger value with broken options allowed
+    #[case(true, f64::MAX, true)] // Valid maximum finite value with broken options allowed
+    #[case(true, -1e-10, false)] // Invalid: negative value
+    #[case(true, f64::INFINITY, false)] // Invalid: positive infinity
+    #[case(true, f64::NEG_INFINITY, false)] // Invalid: negative infinity
+    #[case(true, f64::NAN, false)] // Invalid: NaN
+    #[case(false, -1e-10, false)] // Invalid: negative value
+    #[case(false, f64::INFINITY, false)] // Invalid: positive infinity
+    #[case(false, f64::NEG_INFINITY, false)] // Invalid: negative infinity
+    #[case(false, f64::NAN, false)] // Invalid: NaN
+    fn check_remaining_demand_absolute_tolerance_works(
+        #[case] allow_broken_options: bool,
+        #[case] value: f64,
+        #[case] expected_valid: bool,
+    ) {
+        let flow = Flow::new(value);
+        let result = check_remaining_demand_absolute_tolerance(allow_broken_options, flow);
+
+        assert_validation_result(
+            result,
+            expected_valid,
+            value,
+            "remaining_demand_absolute_tolerance must be a finite number greater than or equal to zero",
+        );
+    }
+
+    #[rstest]
+    #[case(0.0)] // smaller than default
+    #[case(1e-10)] // Larger than default (1e-12)
+    #[case(1.0)] // Well above default
+    #[case(f64::MAX)] // Maximum finite value
+    fn check_remaining_demand_absolute_tolerance_requires_broken_options_if_non_default(
+        #[case] value: f64,
+    ) {
+        let flow = Flow::new(value);
+        let result = check_remaining_demand_absolute_tolerance(false, flow);
+        assert_validation_result(
+            result,
+            false,
+            value,
+            "Setting a remaining_demand_absolute_tolerance different from the default value \
+             of 1e-12 is potentially dangerous, set \
+             please_give_me_broken_results to true if you want to allow this.",
+        );
+    }
+
     #[rstest]
     #[case(0.0, true)] // Valid minimum value
     #[case(0.2, true)] // Valid default value

src/simulation/investment.rs

@@ -9,7 +9,7 @@ use crate::region::RegionID;
 use crate::simulation::CommodityPrices;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
 use crate::units::{Capacity, Dimensionless, Flow, FlowPerCapacity};
-use anyhow::{Context, Result, ensure};
+use anyhow::{Context, Result, bail, ensure};
 use indexmap::IndexMap;
 use itertools::{Itertools, chain};
 use log::debug;
@@ -743,7 +743,10 @@ fn select_best_assets(
     // Iteratively select the best asset until demand is met
     let mut round = 0;
     let mut best_assets: Vec<AssetRef> = Vec::new();
-    while is_any_remaining_demand(&demand) {
+    while is_any_remaining_demand(
+        &demand,
+        model.parameters.remaining_demand_absolute_tolerance,
+    ) {
         ensure!(
             !opt_assets.is_empty(),
             "Failed to meet demand for commodity '{}' in region '{}' with provided investment \
@@ -805,11 +808,24 @@ fn select_best_assets(
         // demand.
         // - known issue with the NPV objective
         // (see https://github.com/EnergySystemsModellingLab/MUSE2/issues/716).
-        ensure!(
-            !outputs_for_opts.is_empty(),
-            "No feasible investment options for commodity '{}' after appraisal",
-            &commodity.id
-        );
+        if outputs_for_opts.is_empty() {
+            let remaining_demands: Vec<_> = demand
+                .iter()
+                .filter(|(_, flow)| **flow > Flow(0.0))
+                .map(|(time_slice, flow)| format!("{} : {:e}", time_slice, flow.value()))
+                .collect();
+
+            bail!(
+                "No feasible investment options left for \
+                commodity '{}', region '{}', year '{}', agent '{}' after appraisal.\n\
+                Remaining unmet demand (time_slice : flow):\n{}",
+                &commodity.id,
+                region_id,
+                year,
+                agent.id,
+                remaining_demands.join("\n")
+            );
+        }
 
         // Warn if there are multiple equally good assets
         warn_on_equal_appraisal_outputs(&outputs_for_opts, &agent.id, &commodity.id, region_id);
@@ -851,8 +867,8 @@ fn select_best_assets(
 }
 
 /// Check whether there is any remaining demand that is unmet in any time slice
-fn is_any_remaining_demand(demand: &DemandMap) -> bool {
-    demand.values().any(|flow| *flow > Flow(0.0))
+fn is_any_remaining_demand(demand: &DemandMap, absolute_tolerance: Flow) -> bool {
+    demand.values().any(|flow| *flow > absolute_tolerance)
 }
 
 /// Update capacity of chosen asset, if needed, and update both asset options and chosen assets