diff --git a/docs/model/dispatch_optimisation.md b/docs/model/dispatch_optimisation.md index 6609d9dad..0030b3ed5 100644 --- a/docs/model/dispatch_optimisation.md +++ b/docs/model/dispatch_optimisation.md @@ -125,6 +125,28 @@ where: - For **Service Demand** (`SVD`): \\( \mathrm{Demand}\_{c, r, s} \\) - For **Supply-Equals-Demand** (`SED`): \\( 0 \\) +### Commodity Consumption/Production Constraints + +Commodity constraints impose lower and upper limits on the total production or consumption of a +commodity in a region over a specified time slice selection: + +\\[ + L\_{c,r,s} \leq + \sum\_{a \in \mathbf{A}\_r^d} |f\_{\mathrm{coeff},a,c}| \cdot + \sum\_{t \in s} \mathrm{Activity}\_{a,t} + \leq U\_{c,r,s} +\\] + +where: + +- \\( d \\) is the balance type: production (`prod`) or consumption (`cons`). +- \\( \mathbf{A}\_r^d \\) contains assets in region \\( r \\) with flows in direction \\( d \\). +- \\( L\_{c,r,s} \\) and \\( U\_{c,r,s} \\) are the lower and upper limits. + +These constraints are defined in the optional `commodity_constraints.csv` file. They can apply to +`SED` and `OTH` commodities, but not `SVD` commodities. The feature is experimental and requires +`please_give_me_broken_results = true` in `model.toml`. + ## Shadow Prices The dual values (shadow prices) of the commodity balance constraints represent the marginal cost of @@ -230,24 +252,23 @@ candidate dispatch run are then used to seed and guide investment appraisal in s ## Diagnosing Infeasible Models -In practice, a dispatch optimisation run can fail if the problem is **infeasible** — typically -because the installed asset capacity in the region is insufficient to meet the required exogenous or -intermediate commodity demands. +In practice, a dispatch optimisation run may be **infeasible** for several reasons, such as +insufficient installed asset capacity to meet demand or incompatible commodity production/consumption +constraints. When this occurs, MUSE2 performs additional dispatch runs with modified optimisation +problems to help identify the cause. The resulting diagnostic information is included in error +messages and saved in the dispatch debug files. -To help debug and pinpoint the exact source of failure, MUSE2 employs a diagnostic mechanism using -**unmet demand variables**: +### Unmet Demand Diagnostic -1. **First-Pass Run:** MUSE2 first attempts to solve the dispatch model in its standard form -(without unmet demand variables). -2. **Diagnostic Re-Run:** If the solver reports that the problem is infeasible, MUSE2 automatically -spawns a second, diagnostic dispatch run. In this run, a set of slack variables representing unmet -demand, \\( \mathrm{UnmetD}\_{c, r, t} \ge 0 \\), is added to the commodity balance constraints: +1. **Diagnostic Re-Run:** MUSE2 reruns the dispatch optimisation with a set of slack variables +representing unmet demand, \\( \\mathrm{UnmetD}\_{c, r, t} \\ge 0 \\), added to the commodity balance +constraints: \\[ \sum_{a \in \mathbf{A}\_r} f\_{\mathrm{coeff},a,c} \cdot \sum\_{t \in s} \mathrm{Activity}\_{a, t} + \sum\_{t \in s} \mathrm{UnmetD}\_{c, r, t} \ge \mathrm{Bound}\_{c, r, s} \\] -3. **Objective Penalty:** To ensure the solver only leaves demand unmet if it is mathematically +1. **Objective Penalty:** To ensure the solver only leaves demand unmet if it is mathematically impossible to satisfy it, these variables are heavily penalised in the diagnostic objective function using the `value_of_lost_load` parameter (\\( \mathrm{VoLL} \\)): \\[ @@ -255,9 +276,19 @@ using the `value_of_lost_load` parameter (\\( \mathrm{VoLL} \\)): \mathrm{Cost}\_{\mathrm{Activity},a,t} + \mathrm{VoLL} \cdot \sum\_{c, r, t} \mathrm{UnmetD}\_{c, r, t} \\] -4. **Isolating Shortfalls:** The addition of \\( \mathrm{UnmetD}\_{c, r, t} \\) guarantees that the +1. **Isolating Shortfalls:** The addition of \\( \mathrm{UnmetD}\_{c, r, t} \\) guarantees that the LP remains mathematically feasible. When solved, any time slice, region, or commodity with a shortfall will have \\( \mathrm{UnmetD}_{c, r, t} > 0 \\). -5. **Error Reporting:** MUSE2 scans the solution, identifies all balanced markets \\( (c, r) \\) +1. **Error Reporting:** MUSE2 scans the solution, identifies all balanced markets \\( (c, r) \\) where unmet demand occurred, outputs detailed diagnostic CSV files, and aborts the simulation with an error identifying the exact out-of-balance markets. + +### Commodity Constraints Diagnostic + +If the dispatch optimisation remains infeasible, MUSE2 reruns it with the commodity +consumption/production constraints disabled. If this rerun succeeds, the infeasibility is likely +caused by one or more constraints defined in `commodity_constraints.csv`. + +If the rerun remains infeasible, commodity constraints are not identified as the cause. Since +commodity constraints are an experimental feature, this diagnosis should be treated as indicative +rather than definitive. diff --git a/schemas/input/commodity_constraints.yaml b/schemas/input/commodity_constraints.yaml index 6a65446b1..ed3293765 100644 --- a/schemas/input/commodity_constraints.yaml +++ b/schemas/input/commodity_constraints.yaml @@ -3,6 +3,10 @@ description: | Specifies the limits on the total amount of consumption/production of a given commodity in a given region, year and time slice(s). +notes: + - Commodity constraints are currently experimental. To use this file, you must enable the + `please_give_me_broken_results` option in `model.toml`. + fields: - name: commodity_id type: string diff --git a/src/commodity.rs b/src/commodity.rs index 585cf5fdc..24bbbc4c5 100644 --- a/src/commodity.rs +++ b/src/commodity.rs @@ -17,8 +17,8 @@ pub type CommodityMap = IndexMap>; /// A map of [`MoneyPerFlow`]s, keyed by region ID, year and time slice ID for a specific levy pub type CommodityLevyMap = HashMap<(RegionID, u32, TimeSliceID), MoneyPerFlow>; -/// A map of vectors of [`CommodityConstraint`]s, keyed by region ID and year -pub type CommodityConstraintsMap = HashMap<(RegionID, u32), Vec>; +/// A map of vectors of [`CommodityConstraint`]s, keyed by year +pub type CommodityConstraintsMap = HashMap>; /// A map of demand values, keyed by region ID, year and time slice selection pub type DemandMap = HashMap<(RegionID, u32, TimeSliceSelection), Flow>; @@ -127,6 +127,8 @@ pub enum PricingStrategy { /// A constraint imposed on commodity values #[derive(PartialEq, Debug, Clone)] pub struct CommodityConstraint { + /// Region to which the commodity constraint applies + pub region_id: RegionID, /// The balance type for the commodity constraint pub balance_type: BalanceType, /// The time slice selection for the commodity constraint diff --git a/src/input/commodity/constraints.rs b/src/input/commodity/constraints.rs index 0c0401f5e..db40d173c 100644 --- a/src/input/commodity/constraints.rs +++ b/src/input/commodity/constraints.rs @@ -1,4 +1,8 @@ //! Code for reading commodity constraints from a CSV file. +//! +//! The `commodity_constraints.csv` file is optional. If it is provided, the +//! `please_give_me_broken_results` option in `model.toml` must be set to `true` because commodity +//! constraints are experimental. use super::super::{input_err_msg, read_csv_optional}; use crate::commodity::{ BalanceType, Commodity, CommodityConstraint, CommodityConstraintsMap, CommodityID, @@ -6,6 +10,7 @@ use crate::commodity::{ }; use crate::id::{GetIDValue, IDCollection}; use crate::input::{parse_range, parse_year_str}; +use crate::model::{ALLOW_DANGEROUS_OPTION_NAME, dangerous_model_options_enabled}; use crate::region::RegionID; use crate::time_slice::TimeSliceInfo; use crate::units::Flow; @@ -69,14 +74,22 @@ pub fn read_commodity_constraints( ) -> Result> { let file_path = model_dir.join(COMMODITY_CONSTRAINTS_FILE_NAME); let commodity_constraints_csv = read_csv_optional(&file_path)?; - read_commodity_constraints_from_iter( + let commodity_constraints = read_commodity_constraints_from_iter( commodity_constraints_csv, commodities, region_ids, time_slice_info, milestone_years, ) - .with_context(|| input_err_msg(&file_path)) + .with_context(|| input_err_msg(&file_path))?; + + ensure!( + commodity_constraints.is_empty() || dangerous_model_options_enabled(), + "Commodity constraints are currently experimental. To use them, set the \ + {ALLOW_DANGEROUS_OPTION_NAME} option to true." + ); + + Ok(commodity_constraints) } /// Process raw commodity-constraint records into a constraints map. @@ -84,7 +97,7 @@ pub fn read_commodity_constraints( /// # Arguments /// /// * `iter` - Iterator over `CommodityConstraintRaw` records -/// * `commodities` - The commodoties in the model +/// * `commodities` - The commodities in the model /// * `region_ids` - All possible region IDs /// * `time_slice_info` - Information about time slices /// * `milestone_years` - All milestone years @@ -124,12 +137,13 @@ where let commodity_map = map.entry(commodity_id.clone()).or_default(); for year in &years { let constraint = CommodityConstraint { + region_id: region_id.clone(), balance_type: record.balance_type.clone(), ts_selection: ts_selection.clone(), limits: limits.clone(), }; commodity_map - .entry((region_id.clone(), *year)) + .entry(*year) .and_modify(|constraints| constraints.push(constraint.clone())) .or_insert(vec![constraint]); } @@ -162,8 +176,7 @@ mod tests { #[test] fn validate_constraints_valid() { - let valid = validate_raw_constraint(BalanceType::Production); - valid.unwrap(); + validate_raw_constraint(BalanceType::Production).unwrap(); } #[test] @@ -177,128 +190,136 @@ mod tests { } #[test] - fn read_commodity_constraints_success() -> Result<()> { - // Create a model dir and write simple CSV files - let dir = tempdir()?; + #[allow(clippy::too_many_lines)] + fn read_commodity_constraints_from_iter_success() { + // Create a model dir and write a simple commodities CSV file + let dir = tempdir().unwrap(); let model_dir = dir.path(); - - // Create simple commodity constraints CSV file - let constraints_csv = concat!( - "commodity_id,region_id,balance_type,years,time_slice,limits\n", - "ELCTRI,GBR,cons,2030,summer,12.34..56.78\n", - "CO2EMT,GBR,cons,2030,winter,..9.99\n", - "CO2EMT,GBR,prod,2030,summer,9.99..\n", - ); - fs::write( - model_dir.join(COMMODITY_CONSTRAINTS_FILE_NAME), - constraints_csv, - )?; - - // Create simple commodities CSV to simplify creating `Commodity`s let commodities_csv = concat!( "id,description,type,time_slice_level,units\n", "ELCTRI,Electricity,sed,season,PJ\n", "CO2EMT,CO2 emitted,oth,season,ktCO2\n", ); - fs::write(model_dir.join(COMMODITY_FILE_NAME), commodities_csv)?; + fs::write(model_dir.join(COMMODITY_FILE_NAME), commodities_csv).unwrap(); // Create basic model inputs let commodities = read_commodities_file(model_dir).unwrap(); - let mut region_ids: IndexSet = IndexSet::new(); - region_ids.insert(RegionID::from("GBR")); - - let time_slice1 = TimeSliceID { - season: "summer".into(), - time_of_day: "day".into(), - }; - let time_slice2 = TimeSliceID { - season: "summer".into(), - time_of_day: "night".into(), - }; - let time_slice3 = TimeSliceID { - season: "winter".into(), - time_of_day: "day".into(), - }; - let time_slice4 = TimeSliceID { - season: "winter".into(), - time_of_day: "night".into(), - }; + let region_ids: IndexSet = ["GBR".into()].into_iter().collect(); let time_slice_info = TimeSliceInfo { seasons: [("summer".into(), Year(0.5)), ("winter".into(), Year(0.5))].into(), times_of_day: ["day".into(), "night".into()].into(), time_slices: [ - (time_slice1.clone(), Year(0.25)), - (time_slice2.clone(), Year(0.25)), - (time_slice3.clone(), Year(0.25)), - (time_slice4.clone(), Year(0.25)), + ( + TimeSliceID { + season: "summer".into(), + time_of_day: "day".into(), + }, + Year(0.25), + ), + ( + TimeSliceID { + season: "summer".into(), + time_of_day: "night".into(), + }, + Year(0.25), + ), + ( + TimeSliceID { + season: "winter".into(), + time_of_day: "day".into(), + }, + Year(0.25), + ), + ( + TimeSliceID { + season: "winter".into(), + time_of_day: "night".into(), + }, + Year(0.25), + ), ] .into(), }; let milestone_years = vec![2030]; + let constraints = [ + CommodityConstraintRaw { + commodity_id: "ELCTRI".into(), + region_id: "GBR".into(), + balance_type: BalanceType::Consumption, + years: "2030".into(), + time_slice: "summer".into(), + limits: "12.34..56.78".into(), + }, + CommodityConstraintRaw { + commodity_id: "CO2EMT".into(), + region_id: "GBR".into(), + balance_type: BalanceType::Consumption, + years: "2030".into(), + time_slice: "winter".into(), + limits: "..9.99".into(), + }, + CommodityConstraintRaw { + commodity_id: "CO2EMT".into(), + region_id: "GBR".into(), + balance_type: BalanceType::Production, + years: "2030".into(), + time_slice: "summer".into(), + limits: "9.99..".into(), + }, + ]; + // Create the constraints map - let constraints_map = read_commodity_constraints( - model_dir, + let constraints_map = read_commodity_constraints_from_iter( + constraints.into_iter(), &commodities, ®ion_ids, &time_slice_info, &milestone_years, - )?; - - // Check the constraints map contains the expected constraint, keyed by the expected - // commodity id - assert!(constraints_map.contains_key(&CommodityID::from("ELCTRI"))); - assert!(constraints_map.contains_key(&CommodityID::from("CO2EMT"))); + ) + .unwrap(); // ELCTRI constraint - let elctri_constraint = &constraints_map[&CommodityID::from("ELCTRI")]; - let elctri_gbr_2030 = elctri_constraint - .get(&(RegionID::from("GBR"), 2030)) - .unwrap(); - assert_eq!(elctri_gbr_2030[0].balance_type, BalanceType::Consumption); + let elctri = &constraints_map[&CommodityID::from("ELCTRI")][&2030][0]; + assert_eq!(elctri.balance_type, BalanceType::Consumption); assert_eq!( - elctri_gbr_2030[0].ts_selection, - TimeSliceSelection::Season("summer".into()), + elctri.ts_selection, + TimeSliceSelection::Season("summer".into()) ); - assert_approx_eq!(f64, elctri_gbr_2030[0].limits.start().value(), 12.34); - assert_approx_eq!(f64, elctri_gbr_2030[0].limits.end().value(), 56.78); + assert_approx_eq!(f64, elctri.limits.start().value(), 12.34); + assert_approx_eq!(f64, elctri.limits.end().value(), 56.78); // CO2EMT constraints - let co2emt_constraint = &constraints_map[&CommodityID::from("CO2EMT")]; - let co2emt_gbr_2030 = co2emt_constraint - .get(&(RegionID::from("GBR"), 2030)) - .unwrap(); - assert_eq!(co2emt_gbr_2030[0].balance_type, BalanceType::Consumption); + let co2emt = &constraints_map[&CommodityID::from("CO2EMT")][&2030]; + assert_eq!(co2emt[0].balance_type, BalanceType::Consumption); assert_eq!( - co2emt_gbr_2030[0].ts_selection, - TimeSliceSelection::Season("winter".into()), + co2emt[0].ts_selection, + TimeSliceSelection::Season("winter".into()) ); - assert_approx_eq!(f64, co2emt_gbr_2030[0].limits.start().value(), 0.0); - assert_approx_eq!(f64, co2emt_gbr_2030[0].limits.end().value(), 9.99); + assert_approx_eq!(f64, co2emt[0].limits.start().value(), 0.0); + assert_approx_eq!(f64, co2emt[0].limits.end().value(), 9.99); - assert_eq!(co2emt_gbr_2030[1].balance_type, BalanceType::Production); + assert_eq!(co2emt[1].balance_type, BalanceType::Production); assert_eq!( - co2emt_gbr_2030[1].ts_selection, - TimeSliceSelection::Season("summer".into()), + co2emt[1].ts_selection, + TimeSliceSelection::Season("summer".into()) ); - assert_approx_eq!(f64, co2emt_gbr_2030[1].limits.start().value(), 9.99); - assert_approx_eq!(f64, co2emt_gbr_2030[1].limits.end().value(), f64::INFINITY); - - Ok(()) + assert_approx_eq!(f64, co2emt[1].limits.start().value(), 9.99); + assert_approx_eq!(f64, co2emt[1].limits.end().value(), f64::INFINITY); } #[test] - fn read_commodity_constraints_fails_with_invalid_csv() -> Result<()> { + fn read_commodity_constraints_fails_with_invalid_csv() { // Create a model dir and write invalid CSV content to force // read_commodity_constraints_from_iter failure - let dir = tempdir()?; + let dir = tempdir().unwrap(); let model_dir = dir.path(); // Create invalid commodity constraints CSV content let file_path = model_dir.join(COMMODITY_CONSTRAINTS_FILE_NAME); - fs::write(&file_path, "invalid,commodity,constraints\nbad,row\n")?; + fs::write(&file_path, "invalid,commodity,constraints\nbad,row\n").unwrap(); // Create empty model inputs let commodities: IndexMap = IndexMap::new(); @@ -322,7 +343,5 @@ mod tests { err_text.contains(COMMODITY_CONSTRAINTS_FILE_NAME), "error message should include file name context, got: {err_text}" ); - - Ok(()) } } diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index 06418126f..771f9dac9 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -112,6 +112,7 @@ pub fn perform_agent_investment( // As upstream markets by definition will not yet have producers, we explicitly set // their prices using external values so that they don't appear free let solution = DispatchRun::new(model, &all_selected_assets, year) + .without_commodity_constraints() .with_market_balance_subset(&seen_markets) .with_input_prices(&prices.shadow) .run(&format!("post {market_set} investment"), writer)?; diff --git a/src/simulation/market.rs b/src/simulation/market.rs index 9ed4b6b64..5c997e7d3 100644 --- a/src/simulation/market.rs +++ b/src/simulation/market.rs @@ -303,6 +303,7 @@ pub fn select_assets_for_cycle( // Run dispatch let solution = DispatchRun::new(model, &all_assets, year) + .without_commodity_constraints() .with_market_balance_subset(&markets_to_balance) .with_flexible_capacity_assets( &flexible_capacity_assets, diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index bf7bdef48..879745574 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -14,10 +14,11 @@ use crate::units::{ Activity, Capacity, Dimensionless, Flow, Money, MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow, Year, }; -use anyhow::{Context, Result, anyhow, bail, ensure}; +use anyhow::{Context, Result, anyhow, bail}; use highs::{HighsModelStatus, RowProblem as Problem, Sense}; use indexmap::{IndexMap, IndexSet}; use itertools::{chain, iproduct}; +use log::warn; use std::collections::HashMap; use std::error::Error; use std::ops::Range; @@ -440,6 +441,7 @@ pub struct DispatchRun<'model, 'run> { candidate_assets: &'run [AssetRef], markets_to_balance: &'run [(CommodityID, RegionID)], input_prices: Option<&'run PriceMap>, + include_commodity_constraints: bool, year: u32, capacity_margin: Dimensionless, } @@ -455,6 +457,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { candidate_assets: &[], markets_to_balance: &[], input_prices: None, + include_commodity_constraints: true, year, capacity_margin: Dimensionless(0.0), } @@ -483,6 +486,14 @@ impl<'model, 'run> DispatchRun<'model, 'run> { } } + /// Exclude explicit production and consumption constraints from the dispatch run. + pub fn without_commodity_constraints(self) -> Self { + Self { + include_commodity_constraints: false, + ..self + } + } + /// Only apply commodity balance constraints to the specified subset of markets pub fn with_market_balance_subset( self, @@ -531,30 +542,139 @@ impl<'model, 'run> DispatchRun<'model, 'run> { .map(|prices| filter_input_prices(prices, markets_to_balance)); let input_prices = input_prices_owned.as_ref(); - // Try running dispatch. If it fails because the model is infeasible, it is likely that this - // is due to unmet demand, in this case, we rerun dispatch including extra variables to - // track the unmet demand so we can report the offending markets to users + // First solve the configured dispatch problem. If it is infeasible, run diagnostic solves + // below to distinguish unmet demand from infeasibility caused by explicit constraints. match self.run_without_unmet_demand_variables(markets_to_balance, input_prices) { + // If the run is successful, we write debug info and return the solution Ok(solution) => { - // Normal successful run: write debug info and return writer.write_dispatch_debug_info(self.year, run_description, &solution)?; Ok(solution) } + + // If the problem is infeasible, we run diagnostics to identify the cause and provide a + // more helpful error message. Err(ModelError::NonOptimal(HighsModelStatus::Infeasible)) => { - // Re-run including unmet demand variables so we can record detailed unmet-demand - // debug output before returning an error to the caller. - let solution = self - .run_internal( + // Generic message for infeasibility, to be augmented with more specific diagnostics + // below + let mut diagnoses = vec![ + "The solver has indicated that the dispatch problem is infeasible".to_string(), + ]; + + // Get diagnostic information for unmet demand + if let Some(diagnosis) = self.run_unmet_demand_diagnostic( + markets_to_balance, + input_prices, + run_description, + writer, + )? { + diagnoses.push(diagnosis); + } + + // Get diagnostic information for commodity constraints, if any apply this year. + if self.has_commodity_constraints() + && let Some(diagnosis) = self.run_commodity_constraints_diagnosis( markets_to_balance, - /*allow_unmet_demand=*/ true, input_prices, - ) - .expect("Failed to run dispatch to calculate unmet demand"); + run_description, + writer, + )? + { + diagnoses.push(diagnosis); + } + + // Assemble and return the final error message, which may include multiple diagnoses + bail!("{}.", diagnoses.join(". ")); + } - // Write debug CSVs to help diagnosis - writer.write_dispatch_debug_info(self.year, run_description, &solution)?; + // Other errors are propagated up to the caller + Err(err) => Err(err.into_anyhow()), + } + } + + /// Check whether any explicit commodity constraints apply in the current year. + fn has_commodity_constraints(&self) -> bool { + self.model.commodities.values().any(|commodity| { + commodity + .constraints + .get(&self.year) + .is_some_and(|constraints| !constraints.is_empty()) + }) + } + + /// Diagnose whether explicit commodity constraints cause infeasibility. + fn run_commodity_constraints_diagnosis( + &self, + markets_to_balance: &[(CommodityID, RegionID)], + input_prices: Option<&PriceMap>, + run_description: &str, + writer: &mut DataWriter, + ) -> Result> { + if !self.include_commodity_constraints { + return Ok(None); + } + + warn!("Dispatch optimisation was infeasible; running commodity constraints diagnostic"); - // Collect markets with unmet demand from the solution + match self.run_internal( + markets_to_balance, + /*include_commodity_constraints=*/ false, + /*allow_unmet_demand=*/ false, + input_prices, + ) { + Ok(solution) => { + let diagnostic_run_description = + format!("{run_description} COMMODITY_CONSTRAINTS_DIAGNOSTIC"); + writer.write_dispatch_debug_info( + self.year, + &diagnostic_run_description, + &solution, + )?; + + Ok(Some( + "The infeasibility is likely caused by one or more constraints defined in \ + `commodity_constraints.csv`. Please note that commodity constraints are \ + currently an experimental feature, so this is not necessarily unexpected" + .to_string(), + )) + } + + // The problem remains infeasible without explicit commodity constraints, so they are + // not identified as the cause. Don't return a diagnostic message in this case. + Err(ModelError::NonOptimal(HighsModelStatus::Infeasible)) => Ok(None), + + // Other errors are propagated up to the caller + Err(error) => Err(error.into_anyhow()), + } + } + + /// Re-run the configured problem with unmet-demand variables to identify unmet demand. + fn run_unmet_demand_diagnostic( + &self, + markets_to_balance: &[(CommodityID, RegionID)], + input_prices: Option<&PriceMap>, + run_description: &str, + writer: &mut DataWriter, + ) -> Result> { + warn!("Dispatch optimisation was infeasible; running unmet demand diagnostic"); + + match self.run_internal( + markets_to_balance, + self.include_commodity_constraints, + /*allow_unmet_demand=*/ true, + input_prices, + ) { + Ok(solution) => { + // The diagnostic solution is written only to provide debugging information; it is + // never returned as the result of the original dispatch run. + let diagnostic_run_description = + format!("{run_description} UNMET_DEMAND_DIAGNOSTIC"); + writer.write_dispatch_debug_info( + self.year, + &diagnostic_run_description, + &solution, + )?; + + // Collect markets where the diagnostic solution uses positive unmet demand. let markets: IndexSet<_> = solution .iter_unmet_demand() .filter(|(_, _, _, flow)| *flow > Flow(0.0)) @@ -563,19 +683,22 @@ impl<'model, 'run> DispatchRun<'model, 'run> { }) .collect(); - ensure!( - !markets.is_empty(), - "Model is infeasible, but there was no unmet demand" - ); - - bail!( - "The solver has indicated that the problem is infeasible, probably because \ - the supplied assets could not meet the required demand. Demand was not met \ - for the following markets: {}", - format_items_with_cap(markets) - ); + Ok(Some(if markets.is_empty() { + "No unmet demand was identified".to_string() + } else { + format!( + "Demand was not met for the following markets: {}", + format_items_with_cap(markets) + ) + })) } - Err(err) => Err(err.into_anyhow()), + + // The problem remains infeasible even with unmet demand variables, so unmet demand is + // not identified as the cause. Don't return a diagnostic message in this case. + Err(ModelError::NonOptimal(HighsModelStatus::Infeasible)) => Ok(None), + + // Other errors are propagated up to the caller + Err(error) => Err(error.into_anyhow()), } } @@ -587,6 +710,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { ) -> Result, ModelError> { self.run_internal( markets_to_balance, + self.include_commodity_constraints, /*allow_unmet_demand=*/ false, input_prices, ) @@ -596,6 +720,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { fn run_internal( &self, markets_to_balance: &[(CommodityID, RegionID)], + include_commodity_constraints: bool, allow_unmet_demand: bool, input_prices: Option<&PriceMap>, ) -> Result, ModelError> { @@ -645,6 +770,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { markets_to_balance, self.year, self.candidate_assets, + include_commodity_constraints, ); // Create model and apply any user-supplied HiGHS options to it diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs index c718af91c..a20c7983e 100644 --- a/src/simulation/optimisation/constraints.rs +++ b/src/simulation/optimisation/constraints.rs @@ -1,8 +1,9 @@ //! Code for adding constraints to the dispatch optimisation problem. use super::VariableMap; use crate::asset::{AssetIterator, AssetRef}; -use crate::commodity::{CommodityID, CommodityType}; +use crate::commodity::{BalanceType, CommodityID, CommodityType}; use crate::model::Model; +use crate::process::FlowDirection; use crate::region::RegionID; use crate::time_slice::{Season, TimeSliceInfo, TimeSliceSelection}; use crate::units::{Flow, MoneyPerCapacityPerYear, UnitType, Year}; @@ -80,6 +81,7 @@ pub struct ConstraintKeys { /// # Returns /// /// Keys for the different constraints. +#[allow(clippy::too_many_arguments)] pub fn add_model_constraints<'a, I>( problem: &mut Problem, variables: &VariableMap, @@ -88,6 +90,7 @@ pub fn add_model_constraints<'a, I>( markets_to_balance: &'a [(CommodityID, RegionID)], year: u32, candidate_assets: &'a [AssetRef], + include_commodity_constraints: bool, ) -> ConstraintKeys where I: Iterator + Clone + 'a, @@ -102,6 +105,10 @@ where candidate_assets, ); + if include_commodity_constraints { + add_commodity_constraints(problem, variables, model, assets, year); + } + let activity_keys = add_activity_constraints(problem, variables, &model.time_slice_info, assets.clone()); @@ -116,6 +123,59 @@ where } } +/// Add explicit production and consumption constraints for commodities. +fn add_commodity_constraints<'a, I>( + problem: &mut Problem, + variables: &VariableMap, + model: &'a Model, + assets: &I, + year: u32, +) where + I: Iterator + Clone + 'a, +{ + // Commodity constraints are indexed by milestone year, so commodities without a constraint + // for this year do not contribute any rows. + for commodity in model.commodities.values() { + let Some(constraints) = commodity.constraints.get(&year) else { + continue; + }; + + for constraint in constraints { + // Select the flow direction represented by the constraint and normalise the + // coefficient sign used in the solver row. + let (flow_direction, coefficient_sign) = match constraint.balance_type { + BalanceType::Production => (FlowDirection::Output, 1.0), + // Input flow coefficients are negative, but the constraint limits describe + // consumption as positive. + BalanceType::Consumption => (FlowDirection::Input, -1.0), + BalanceType::Net => unreachable!("Net commodity constraints are invalid"), + }; + + // Build one term for every matching asset and every time slice in the selection. + let terms = assets + .clone() + .filter_region(&constraint.region_id) + .flows_for_commodity(&commodity.id) + .filter(|(_, flow)| flow.direction() == flow_direction) + .flat_map(|(asset, flow)| { + let coefficient = coefficient_sign * flow.coeff.value(); + constraint.ts_selection.iter(&model.time_slice_info).map( + move |(time_slice, _)| { + (variables.get_activity_var(asset, time_slice), coefficient) + }, + ) + }) + .collect::>(); + + // Apply the configured inclusive lower and upper limits to the sum of the terms. + problem.add_row( + constraint.limits.start().value()..=constraint.limits.end().value(), + terms, + ); + } + } +} + /// Add seasonal and annual utilisation peak constraints to the problem. fn add_utilisation_peak_constraints<'a, I>( problem: &mut Problem, diff --git a/tests/model.rs b/tests/model.rs new file mode 100644 index 000000000..facf0cc3b --- /dev/null +++ b/tests/model.rs @@ -0,0 +1,31 @@ +//! Integration tests for model loading and simulation. +use muse2::input::load_model; +use muse2::patch::{FilePatch, ModelPatch}; +use muse2::simulation; +use tempfile::tempdir; + +#[test] +fn commodity_constraints_infeasibility_is_reported() { + // The `missing_commodity` model has no BIOPRD-producing assets in the base year, so + // enforcing positive production of BIOPRD should make the model infeasible. + let model_dir = ModelPatch::from_example("missing_commodity") + .with_toml_patch("please_give_me_broken_results = true") + .with_file_patch( + FilePatch::new("commodity_constraints.csv").with_replacement(&[ + "commodity_id,region_id,balance_type,years,time_slice,limits", + "BIOPRD,GBR,prod,2020,annual,0.0001..", + ]), + ) + .build_to_tempdir() + .unwrap(); + let model = load_model(model_dir.path()).unwrap(); + let output_dir = tempdir().unwrap(); + + let error = simulation::run(&model, output_dir.path(), true).unwrap_err(); + let message = format!("{error:#}"); + + assert!( + message.contains("The infeasibility is likely caused by one or more constraints defined in `commodity_constraints.csv`"), + "{message}" + ); +}