diff --git a/docs/release_notes/upcoming.md b/docs/release_notes/upcoming.md index 99432f3f0..c77ad6cb1 100644 --- a/docs/release_notes/upcoming.md +++ b/docs/release_notes/upcoming.md @@ -30,6 +30,8 @@ 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 @@ -37,3 +39,4 @@ ready to be released, carry out the following steps: [#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 diff --git a/schemas/input/model.yaml b/schemas/input/model.yaml index 1380ed978..04d90d97a 100644 --- a/schemas/input/model.yaml +++ b/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: | diff --git a/src/finance.rs b/src/finance.rs index bdb296576..5b3de0a59 100644 --- a/src/finance.rs +++ b/src/finance.rs @@ -57,6 +57,11 @@ pub fn profitability_index( activity: &IndexMap, activity_surpluses: &IndexMap, ) -> 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; diff --git a/src/model/parameters.rs b/src/model/parameters.rs index bd0dd3c9b..e2b9da69c 100644 --- a/src/model/parameters.rs +++ b/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(), } } diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index 959a22f75..489bfe3a3 100644 --- a/src/simulation/investment.rs +++ b/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, @@ -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>, + demand: &DemandMap, +) -> Vec { + // 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, diff --git a/src/simulation/investment/appraisal.rs b/src/simulation/investment/appraisal.rs index b509b1a7b..9c082a89e 100644 --- a/src/simulation/investment/appraisal.rs +++ b/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, + 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 { + 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, + 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) { - 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) -> 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;