diff --git a/examples/circularity/agent_commodity_portions.csv b/examples/circularity/agent_commodity_portions.csv new file mode 100644 index 000000000..7b34ad1d2 --- /dev/null +++ b/examples/circularity/agent_commodity_portions.csv @@ -0,0 +1,12 @@ +agent_id,commodity_id,years,commodity_portion +A0_OAG,OILCRD,all,1 +A0_OAG,GASPRD,all,1 +A0_OAG,GASNAT,all,1 +A0_REF,GASOLI,all,1 +A0_REF,DIESEL,all,1 +A0_BPD,BIOPRD,all,1 +A0_BPL,BIOPEL,all,1 +A0_ELC,ELCTRI,all,1 +A0_ELC,H2YPRD,all,1 +A0_RES,RSHEAT,all,1 +A0_TRA,TPASKM,all,1 diff --git a/examples/circularity/agent_cost_limits.csv b/examples/circularity/agent_cost_limits.csv new file mode 100644 index 000000000..9926ae85e --- /dev/null +++ b/examples/circularity/agent_cost_limits.csv @@ -0,0 +1,8 @@ +agent_id,years,capex_limit,annual_cost_limit +A0_OAG,all,, +A0_ELC,all,, +A0_RES,all,, +A0_TRA,all,, +A0_BPD,all,, +A0_BPL,all,, +A0_REF,all,, diff --git a/examples/circularity/agent_objectives.csv b/examples/circularity/agent_objectives.csv new file mode 100644 index 000000000..c1482b68b --- /dev/null +++ b/examples/circularity/agent_objectives.csv @@ -0,0 +1,8 @@ +agent_id,years,objective_type,decision_weight,decision_lexico_order +A0_OAG,all,lcox,, +A0_ELC,all,lcox,, +A0_RES,all,lcox,, +A0_TRA,all,lcox,, +A0_BPD,all,lcox,, +A0_BPL,all,lcox,, +A0_REF,all,lcox,, diff --git a/examples/circularity/agent_search_space.csv b/examples/circularity/agent_search_space.csv new file mode 100644 index 000000000..13f6cf903 --- /dev/null +++ b/examples/circularity/agent_search_space.csv @@ -0,0 +1 @@ +agent_id,commodity_id,years,search_space diff --git a/examples/circularity/agents.csv b/examples/circularity/agents.csv new file mode 100644 index 000000000..ceebf51cd --- /dev/null +++ b/examples/circularity/agents.csv @@ -0,0 +1,8 @@ +id,description,regions,decision_rule,decision_lexico_tolerance +A0_OAG,Oil and gas IOC,all,single, +A0_ELC,Electricity generators,all,single, +A0_RES,Residential consumer,all,single, +A0_TRA,Consumer who drives,all,single, +A0_BPD,Biomass producer,all,single, +A0_BPL,Biomass pelletiser,all,single, +A0_REF,Refiners,all,single, diff --git a/examples/circularity/assets.csv b/examples/circularity/assets.csv new file mode 100644 index 000000000..5ff134914 --- /dev/null +++ b/examples/circularity/assets.csv @@ -0,0 +1,16 @@ +process_id,region_id,agent_id,capacity,commission_year +GASDRV,GBR,A0_OAG,4010.2,2020 +OAGRSV,GBR,A0_OAG,1738.05,2020 +GASPRC,GBR,A0_OAG,3789.64,2020 +OILREF,GBR,A0_OAG,716.5,2020 +OILRF2,GBR,A0_OAG,573.9,2020 +WNDFRM,GBR,A0_ELC,3.964844,2020 +GASCGT,GBR,A0_ELC,2.999,2020 +H2YGEN,GBR,A0_ELC,1.0,2020 +H2YPRO,GBR,A0_ELC,31.54,2020 +TPETCR,GBR,A0_TRA,451.91,2020 +TDIECR,GBR,A0_TRA,176.58,2020 +TELCCR,GBR,A0_TRA,20.99,2020 +THYBCR,GBR,A0_TRA,15.0,2020 +RGASBR,GBR,A0_RES,2900,2020 +RELCHP,GBR,A0_RES,399.98,2020 diff --git a/examples/circularity/commodities.csv b/examples/circularity/commodities.csv new file mode 100644 index 000000000..092ca5baa --- /dev/null +++ b/examples/circularity/commodities.csv @@ -0,0 +1,13 @@ +id,description,type,time_slice_level +OILCRD,Crude oil,sed,season +GASPRD,Gas produced,sed,season +GASOLI,Gasoline,sed,season +DIESEL,Diesel,sed,season +GASNAT,Natural gas,sed,season +ELCTRI,Electricity,sed,daynight +H2YPRD,Hydrogen,sed,season +TPASKM,Passenger kms,svd,daynight +RSHEAT,Residential heating,svd,daynight +CO2EMT,CO2 emitted,oth,annual +BIOPRD,Biomass produced,sed,season +BIOPEL,Biomass pellets,sed,season diff --git a/examples/circularity/commodity_levies.csv b/examples/circularity/commodity_levies.csv new file mode 100644 index 000000000..00e64fa89 --- /dev/null +++ b/examples/circularity/commodity_levies.csv @@ -0,0 +1,2 @@ +commodity_id,regions,years,time_slice,balance_type,value +CO2EMT,all,all,annual,net,0.04 diff --git a/examples/circularity/demand.csv b/examples/circularity/demand.csv new file mode 100644 index 000000000..4d31a7be9 --- /dev/null +++ b/examples/circularity/demand.csv @@ -0,0 +1,7 @@ +commodity_id,region_id,year,demand +RSHEAT,GBR,2020,927.38 +RSHEAT,GBR,2030,1027.38 +RSHEAT,GBR,2040,1127.38 +TPASKM,GBR,2020,473.2593 +TPASKM,GBR,2030,573.2593 +TPASKM,GBR,2040,673.2593 diff --git a/examples/circularity/demand_slicing.csv b/examples/circularity/demand_slicing.csv new file mode 100644 index 000000000..9957c6097 --- /dev/null +++ b/examples/circularity/demand_slicing.csv @@ -0,0 +1,33 @@ +commodity_id,region_id,time_slice,fraction +RSHEAT,GBR,winter.night,0.065498087 +RSHEAT,GBR,winter.day,0.226148097 +RSHEAT,GBR,winter.peak,0.111199678 +RSHEAT,GBR,winter.evening,0.076688397 +RSHEAT,GBR,peak.night,0.045882959 +RSHEAT,GBR,peak.day,0.121987853 +RSHEAT,GBR,peak.peak,0.068135438 +RSHEAT,GBR,peak.evening,0.041334529 +RSHEAT,GBR,summer.night,0.005698366 +RSHEAT,GBR,summer.day,0.009494171 +RSHEAT,GBR,summer.peak,0.003126411 +RSHEAT,GBR,summer.evening,0.002327863 +RSHEAT,GBR,autumn.night,0.037269058 +RSHEAT,GBR,autumn.day,0.098086985 +RSHEAT,GBR,autumn.peak,0.054500569 +RSHEAT,GBR,autumn.evening,0.032621538 +TPASKM,GBR,winter.peak,0.042372881 +TPASKM,GBR,winter.night,0.009887006 +TPASKM,GBR,winter.evening,0.056497175 +TPASKM,GBR,winter.day,0.141242938 +TPASKM,GBR,summer.peak,0.042372881 +TPASKM,GBR,summer.night,0.009887006 +TPASKM,GBR,summer.evening,0.056497175 +TPASKM,GBR,summer.day,0.141242938 +TPASKM,GBR,peak.peak,0.042372881 +TPASKM,GBR,peak.night,0.009887006 +TPASKM,GBR,peak.evening,0.056497175 +TPASKM,GBR,peak.day,0.141242938 +TPASKM,GBR,autumn.peak,0.042372881 +TPASKM,GBR,autumn.night,0.009887006 +TPASKM,GBR,autumn.evening,0.056497175 +TPASKM,GBR,autumn.day,0.141242938 diff --git a/examples/circularity/model.toml b/examples/circularity/model.toml new file mode 100644 index 000000000..98be065bf --- /dev/null +++ b/examples/circularity/model.toml @@ -0,0 +1 @@ +milestone_years = [2020, 2030, 2040] diff --git a/examples/circularity/process_availabilities.csv b/examples/circularity/process_availabilities.csv new file mode 100644 index 000000000..652cff992 --- /dev/null +++ b/examples/circularity/process_availabilities.csv @@ -0,0 +1,34 @@ +process_id,regions,commission_years,time_slice,limit_type,value +GASDRV,all,all,annual,up,0.9 +OAGRSV,all,all,annual,up,0.9 +GASPRC,all,all,annual,up,0.9 +OILREF,all,all,annual,up,0.9 +OILRF2,all,all,annual,up,0.9 +GASCGT,all,all,annual,up,0.9 +TPETCR,all,all,annual,up,1 +TDIECR,all,all,annual,up,1 +TELCCR,all,all,annual,up,1 +THYBCR,all,all,annual,up,1 +RGASBR,all,all,annual,up,1 +RELCHP,all,all,annual,up,1 +WNDFRM,all,all,winter.night,up,0.486418015 +WNDFRM,all,all,winter.day,up,0.543166784 +WNDFRM,all,all,winter.peak,up,0.504433498 +WNDFRM,all,all,winter.evening,up,0.493173821 +WNDFRM,all,all,peak.night,up,0.312697296 +WNDFRM,all,all,peak.day,up,0.489120338 +WNDFRM,all,all,peak.peak,up,0.454890922 +WNDFRM,all,all,peak.evening,up,0.331034483 +WNDFRM,all,all,summer.night,up,0.17951141 +WNDFRM,all,all,summer.day,up,0.349950739 +WNDFRM,all,all,summer.peak,up,0.342294159 +WNDFRM,all,all,summer.evening,up,0.202674173 +WNDFRM,all,all,autumn.night,up,0.3513019 +WNDFRM,all,all,autumn.day,up,0.460745954 +WNDFRM,all,all,autumn.peak,up,0.396340605 +WNDFRM,all,all,autumn.evening,up,0.364813512 +H2YGEN,all,all,annual,up,1 +H2YPRO,all,all,annual,up,1 +BIOPRO,all,all,annual,up,1.0 +BIOPLL,all,all,annual,up,0.95 +RBIOBL,all,all,annual,up,1.0 diff --git a/examples/circularity/process_flows.csv b/examples/circularity/process_flows.csv new file mode 100644 index 000000000..51b586741 --- /dev/null +++ b/examples/circularity/process_flows.csv @@ -0,0 +1,45 @@ +process_id,commodity_id,regions,commission_years,coeff,type,cost +GASDRV,GASPRD,all,all,1,fixed, +OAGRSV,OILCRD,all,all,1,fixed, +OAGRSV,GASPRD,all,all,0.1,fixed, +GASPRC,GASPRD,all,all,-1.05,fixed, +GASPRC,GASNAT,all,all,1,fixed, +OILREF,OILCRD,all,all,-1.3,fixed, +OILREF,GASOLI,all,all,0.5,fixed, +OILREF,DIESEL,all,all,0.5,fixed, +OILRF2,OILCRD,all,all,-1.31,fixed, +OILRF2,GASOLI,all,all,1,fixed, +WNDFRM,ELCTRI,all,all,1,fixed, +GASCGT,GASNAT,all,all,-1.5,fixed, +GASCGT,ELCTRI,all,all,1,fixed, +H2YGEN,H2YPRD,all,all,-1.5,fixed, +H2YGEN,ELCTRI,all,all,1,fixed, +H2YPRO,ELCTRI,all,all,-1.333,fixed, +H2YPRO,H2YPRD,all,all,1,fixed, +TPETCR,GASOLI,all,all,-2.702702703,fixed, +TPETCR,TPASKM,all,all,1,fixed, +TDIECR,DIESEL,all,all,-2.222222222,fixed, +TDIECR,TPASKM,all,all,1,fixed, +TELCCR,ELCTRI,all,all,-1.754385965,fixed, +TELCCR,TPASKM,all,all,1,fixed, +THYBCR,ELCTRI,all,all,-0.854385965,fixed, +THYBCR,GASOLI,all,all,-1.0,fixed, +THYBCR,TPASKM,all,all,1,fixed, +RGASBR,GASNAT,all,all,-1.15,fixed, +RGASBR,RSHEAT,all,all,1,fixed, +RELCHP,ELCTRI,all,all,-0.33,fixed, +RELCHP,RSHEAT,all,all,1,fixed, +GASDRV,CO2EMT,all,all,5.113,fixed, +OAGRSV,CO2EMT,all,all,7.333,fixed, +GASPRC,CO2EMT,all,all,2.5565,fixed, +OILREF,CO2EMT,all,all,21.999,fixed, +OILRF2,CO2EMT,all,all,22,fixed, +GASCGT,CO2EMT,all,all,76.695,fixed, +TPETCR,CO2EMT,all,all,179.7297297,fixed, +TDIECR,CO2EMT,all,all,154.0222222,fixed, +RGASBR,CO2EMT,all,all,58.7995,fixed, +BIOPRO,BIOPRD,all,all,1.0,fixed, +BIOPLL,BIOPRD,all,all,-1.05,fixed, +BIOPLL,BIOPEL,all,all,1.0,fixed, +RBIOBL,BIOPEL,all,all,-1.2,fixed, +RBIOBL,RSHEAT,all,all,1,fixed, diff --git a/examples/circularity/process_parameters.csv b/examples/circularity/process_parameters.csv new file mode 100644 index 000000000..8af68f420 --- /dev/null +++ b/examples/circularity/process_parameters.csv @@ -0,0 +1,19 @@ +process_id,regions,commission_years,capital_cost,fixed_operating_cost,variable_operating_cost,lifetime,discount_rate +GASDRV,GBR,all,10,0.3,2,25,0.1 +OAGRSV,GBR,all,15,0.45,3,25,0.1 +GASPRC,GBR,all,7,0.21,0.5,25,0.1 +OILREF,GBR,all,9,0.27,0.66,25,0.1 +OILRF2,GBR,all,9.1,0.3,0.68,25,0.1 +WNDFRM,GBR,all,1000,30,0.4,25,0.1 +GASCGT,GBR,all,700,21,0.55,30,0.1 +H2YGEN,GBR,all,1200,1,0.6,30,0.1 +H2YPRO,GBR,all,15,0.9,0.6,20,0.1 +TPETCR,GBR,all,2,0.06,0.2,10,0.1 +TDIECR,GBR,all,2.1,0.063,0.18,10,0.1 +TELCCR,GBR,all,6,0.18,0.15,10,0.1 +THYBCR,GBR,all,5,0.17,0.17,10,0.1 +RGASBR,GBR,all,55.56,1.6668,0.16,15,0.1 +RELCHP,GBR,all,138.9,4.167,0.17,15,0.1 +BIOPRO,all,all,1,0.2,0.25,20,0.09 +BIOPLL,all,all,2,0.22,0.26,20,0.1 +RBIOBL,all,all,60,1.05,0.2,20,0.1 diff --git a/examples/circularity/processes.csv b/examples/circularity/processes.csv new file mode 100644 index 000000000..1afc63624 --- /dev/null +++ b/examples/circularity/processes.csv @@ -0,0 +1,19 @@ +id,description,regions,primary_output,start_year,end_year,capacity_to_activity +GASDRV,Dry gas extraction,GBR,GASPRD,2020,,1.0 +OAGRSV,Oil and associated gas extraction,GBR,OILCRD,2020,,1.0 +GASPRC,Gas processing,GBR,GASNAT,2020,,1.0 +OILREF,Oil refinery,GBR,DIESEL,2020,,1.0 +OILRF2,O2G refinery,GBR,GASOLI,2020,,1.0 +WNDFRM,Wind farm,GBR,ELCTRI,2020,,31.54 +GASCGT,Gas combined cycle turbine,GBR,ELCTRI,2020,,31.54 +H2YGEN,H2 turbine,GBR,ELCTRI,2020,,31.54 +H2YPRO,H2 electrolyer,GBR,H2YPRD,2020,,1.0 +TPETCR,Petrol car,GBR,TPASKM,2020,,1.0 +TDIECR,Diesel car,GBR,TPASKM,2020,,1.0 +TELCCR,Electric car,GBR,TPASKM,2020,,1.0 +THYBCR,Plug-in hybrid car,GBR,TPASKM,2020,,1.0 +RGASBR,Gas boiler,GBR,RSHEAT,2020,,1.0 +RELCHP,Heat pump,GBR,RSHEAT,2020,,1.0 +BIOPRO,Biomass production,all,BIOPRD,2020,,1.0 +BIOPLL,Biomass pelletiser,all,BIOPEL,2020,,1.0 +RBIOBL,Biomass boiler,all,RSHEAT,2020,,1.0 diff --git a/examples/circularity/regions.csv b/examples/circularity/regions.csv new file mode 100644 index 000000000..a1d2d74c5 --- /dev/null +++ b/examples/circularity/regions.csv @@ -0,0 +1,2 @@ +id,description +GBR,United Kingdom diff --git a/examples/circularity/time_slices.csv b/examples/circularity/time_slices.csv new file mode 100644 index 000000000..76583f12c --- /dev/null +++ b/examples/circularity/time_slices.csv @@ -0,0 +1,17 @@ +season,time_of_day,fraction +winter,night,0.072916667 +winter,day,0.104166667 +winter,peak,0.03125 +winter,evening,0.041666667 +peak,night,0.072916667 +peak,day,0.104166667 +peak,peak,0.03125 +peak,evening,0.041666667 +summer,night,0.072916667 +summer,day,0.104166667 +summer,peak,0.03125 +summer,evening,0.041666667 +autumn,night,0.072916667 +autumn,day,0.104166667 +autumn,peak,0.03125 +autumn,evening,0.041666667 diff --git a/src/asset.rs b/src/asset.rs index 6d40b9e8d..af1df7bd8 100644 --- a/src/asset.rs +++ b/src/asset.rs @@ -451,11 +451,14 @@ impl Asset { self.capacity } - /// Set the capacity for this asset (only for Candidate assets) + /// Set the capacity for this asset (only for Candidate or Selected assets) pub fn set_capacity(&mut self, capacity: Capacity) { assert!( - self.state == AssetState::Candidate, - "set_capacity can only be called on Candidate assets" + matches!( + self.state, + AssetState::Candidate | AssetState::Selected { .. } + ), + "set_capacity can only be called on Candidate or Selected assets" ); assert!(capacity >= Capacity(0.0), "Capacity must be >= 0"); self.capacity = capacity; diff --git a/src/log.rs b/src/log.rs index 75baca563..d659a8b3e 100644 --- a/src/log.rs +++ b/src/log.rs @@ -100,6 +100,7 @@ pub fn init(log_level_from_settings: &str, log_file_path: Option<&Path>) -> Resu // Configure the logger let mut dispatch = Dispatch::new() + .level_for("highs", LevelFilter::Off) // turn off logging for highs .chain( // Write non-error messages to stdout Dispatch::new() diff --git a/src/output.rs b/src/output.rs index 2c1b9e4af..5faa9ae33 100644 --- a/src/output.rs +++ b/src/output.rs @@ -41,6 +41,9 @@ const ACTIVITY_ASSET_DISPATCH: &str = "debug_dispatch_assets.csv"; /// The output file name for commodity balance duals const COMMODITY_BALANCE_DUALS_FILE_NAME: &str = "debug_commodity_balance_duals.csv"; +/// The output file name for unmet demand values +const UNMET_DEMAND_FILE_NAME: &str = "debug_unmet_demand.csv"; + /// The output file name for extra solver output values const SOLVER_VALUES_FILE_NAME: &str = "debug_solver.csv"; @@ -192,6 +195,17 @@ struct CommodityBalanceDualsRow { value: MoneyPerFlow, } +/// Represents the unmet demand data in a row of the unmet demand CSV file +#[derive(Serialize, Deserialize, Debug, PartialEq)] +struct UnmetDemandRow { + milestone_year: u32, + run_description: String, + commodity_id: CommodityID, + region_id: RegionID, + time_slice: TimeSliceID, + value: Flow, +} + /// Represents solver output values #[derive(Serialize, Deserialize, Debug, PartialEq)] struct SolverValuesRow { @@ -232,6 +246,7 @@ struct AppraisalResultsTimeSliceRow { struct DebugDataWriter { context: Option, commodity_balance_duals_writer: csv::Writer, + unmet_demand_writer: csv::Writer, solver_values_writer: csv::Writer, appraisal_results_writer: csv::Writer, appraisal_results_time_slice_writer: csv::Writer, @@ -253,6 +268,7 @@ impl DebugDataWriter { Ok(Self { context: None, commodity_balance_duals_writer: new_writer(COMMODITY_BALANCE_DUALS_FILE_NAME)?, + unmet_demand_writer: new_writer(UNMET_DEMAND_FILE_NAME)?, solver_values_writer: new_writer(SOLVER_VALUES_FILE_NAME)?, appraisal_results_writer: new_writer(APPRAISAL_RESULTS_FILE_NAME)?, appraisal_results_time_slice_writer: new_writer( @@ -290,6 +306,11 @@ impl DebugDataWriter { run_description, solution.iter_commodity_balance_duals(), )?; + self.write_unmet_demand( + milestone_year, + run_description, + solution.iter_unmet_demand(), + )?; self.write_solver_values(milestone_year, run_description, solution.objective_value)?; Ok(()) } @@ -377,6 +398,31 @@ impl DebugDataWriter { Ok(()) } + /// Write unmet demand values to file + fn write_unmet_demand<'a, I>( + &mut self, + milestone_year: u32, + run_description: &str, + iter: I, + ) -> Result<()> + where + I: Iterator, + { + for (commodity_id, region_id, time_slice, value) in iter { + let row = UnmetDemandRow { + milestone_year, + run_description: self.with_context(run_description), + commodity_id: commodity_id.clone(), + region_id: region_id.clone(), + time_slice: time_slice.clone(), + value, + }; + self.unmet_demand_writer.serialize(row)?; + } + + Ok(()) + } + /// Write additional solver output values to file fn write_solver_values( &mut self, @@ -453,6 +499,7 @@ impl DebugDataWriter { /// Flush the underlying streams fn flush(&mut self) -> Result<()> { self.commodity_balance_duals_writer.flush()?; + self.unmet_demand_writer.flush()?; self.solver_values_writer.flush()?; self.appraisal_results_writer.flush()?; self.appraisal_results_time_slice_writer.flush()?; @@ -760,6 +807,48 @@ mod tests { assert_equal(records, iter::once(expected)); } + #[rstest] + fn test_write_unmet_demand( + commodity_id: CommodityID, + region_id: RegionID, + time_slice: TimeSliceID, + ) { + let milestone_year = 2020; + let run_description = "test_run".to_string(); + let value = Flow(0.5); + let dir = tempdir().unwrap(); + + // Write unmet demand + { + let mut writer = DebugDataWriter::create(dir.path()).unwrap(); + writer + .write_unmet_demand( + milestone_year, + &run_description, + iter::once((&commodity_id, ®ion_id, &time_slice, value)), + ) + .unwrap(); + writer.flush().unwrap(); + } + + // Read back and compare + let expected = UnmetDemandRow { + milestone_year, + run_description, + commodity_id: commodity_id, + region_id: region_id, + time_slice, + value, + }; + let records: Vec = + csv::Reader::from_path(dir.path().join(UNMET_DEMAND_FILE_NAME)) + .unwrap() + .into_deserialize() + .try_collect() + .unwrap(); + assert_equal(records, iter::once(expected)); + } + #[rstest] fn test_write_activity(assets: AssetPool, time_slice: TimeSliceID) { let milestone_year = 2020; diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index b10b0c759..a3791cb95 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -11,9 +11,9 @@ use crate::time_slice::{TimeSliceID, TimeSliceInfo}; use crate::units::{Capacity, Dimensionless, Flow, FlowPerCapacity}; use anyhow::{Result, bail, ensure}; use indexmap::IndexMap; -use itertools::{Itertools, chain, iproduct}; -use log::debug; -use std::collections::HashMap; +use itertools::{Itertools, chain}; +use log::{debug, warn}; +use std::collections::{HashMap, VecDeque}; use std::fmt::Display; pub mod appraisal; @@ -57,6 +57,8 @@ impl InvestmentSet { demand: &AllDemandMap, existing_assets: &[AssetRef], prices: &CommodityPrices, + seen_markets: &[(CommodityID, RegionID)], + previously_selected_assets: &[AssetRef], writer: &mut DataWriter, ) -> Result> { match self { @@ -70,9 +72,20 @@ impl InvestmentSet { prices, writer, ), - InvestmentSet::Cycle(_) => bail!( - "Investment cycles are not yet supported. Found cycle for commodities: {self}" - ), + InvestmentSet::Cycle(markets) => { + debug!("Starting investment for cycle '{self}'"); + select_assets_for_cycle( + model, + markets, + year, + demand, + existing_assets, + prices, + seen_markets, + previously_selected_assets, + writer, + ) + } InvestmentSet::Layer(investment_sets) => { debug!("Starting asset selection for layer '{self}'"); let mut all_assets = Vec::new(); @@ -83,6 +96,8 @@ impl InvestmentSet { demand, existing_assets, prices, + seen_markets, + previously_selected_assets, writer, )?; all_assets.extend(assets); @@ -103,7 +118,7 @@ impl Display for InvestmentSet { InvestmentSet::Cycle(markets) => { write!( f, - "{}", + "({})", markets.iter().map(|(c, r)| format!("{c}|{r}")).join(", ") ) } @@ -144,11 +159,6 @@ pub fn perform_agent_investment( investment_order.iter().join(" -> ") ); - // External prices to be used in dispatch optimisation - // Once investments are performed for a market, the dispatch system will be able to produce - // endogenous prices for that market, so we'll gradually remove these external prices. - let mut external_prices = prices.clone(); - // Keep track of the markets that have been seen so far. This will be used to apply // balance constraints in the dispatch optimisation - we only apply balance constraints for // markets that have been seen so far. @@ -163,6 +173,8 @@ pub fn perform_agent_investment( &net_demand, existing_assets, prices, + &seen_markets, + &all_selected_assets, writer, )?; @@ -171,22 +183,11 @@ pub fn perform_agent_investment( seen_markets.push(market.clone()); } - // Remove prices for already-seen markets. Markets which are produced by at - // least one asset in the dispatch run will have prices produced endogenously (via the - // commodity balance constraints), but markets for which investment has not yet been - // performed will, by definition, not have any producers. For these, we provide prices - // from the previous dispatch run otherwise they will appear to be free to the model. - for ((commodity_id, region_id), time_slice) in iproduct!( - investment_set.iter_markets(), - model.time_slice_info.iter_ids() - ) { - external_prices.remove(commodity_id, region_id, time_slice); - } - // If no assets have been selected, skip dispatch optimisation // **TODO**: this probably means there's no demand for the market, which we could // presumably preempt if selected_assets.is_empty() { + debug!("No assets selected for '{investment_set}'"); continue; } @@ -201,8 +202,8 @@ 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) - .with_market_subset(&seen_markets) - .with_input_prices(&external_prices) + .with_market_balance_subset(&seen_markets) + .with_input_prices(prices) .run(&format!("post {investment_set} investment"), writer)?; // Update demand map with flows from newly added assets @@ -279,6 +280,137 @@ fn select_assets_for_single_market( Ok(selected_assets) } +#[allow(clippy::too_many_arguments)] +fn select_assets_for_cycle( + model: &Model, + markets: &[(CommodityID, RegionID)], + year: u32, + demand: &AllDemandMap, + existing_assets: &[AssetRef], + prices: &CommodityPrices, + seen_markets: &[(CommodityID, RegionID)], + previously_selected_assets: &[AssetRef], + writer: &mut DataWriter, +) -> Result> { + const MAX_ITERATIONS: u32 = 5; + const CAPACITY_MARGIN: f64 = 0.1; + + // Precompute a joined string for logging + let markets_str = markets.iter().map(|(c, r)| format!("{c}|{r}")).join(", "); + + // Get markets to balance: all seen so far plus all in the cycle + let mut markets_to_balance = seen_markets.to_vec(); + markets_to_balance.extend_from_slice(markets); + + // Initialise queue of unbalanced markets: start with all markets in the cycle + // TODO: currently the order is unspecified, but we may be able to speed up convergence by + // pre-sorting this in some way + let mut unbalanced_markets = VecDeque::from(markets.to_vec()); + + // Keep track of the number of times we've visited each market in the cycle + let mut markets_visit_count: HashMap<(CommodityID, RegionID), u32> = + markets.iter().cloned().map(|m| (m, 0u32)).collect(); + + // Iterate while there are still markets in the queue, or until max iterations reached + let mut current_demand = demand.clone(); + let mut assets_for_cycle = HashMap::new(); + let mut last_solution = None; + while !unbalanced_markets.is_empty() { + // Pop the first market off the queue + let (commodity_id, region_id) = unbalanced_markets + .pop_front() + .expect("Already checked that the queue is non-empty"); + let count = markets_visit_count + .get_mut(&(commodity_id.clone(), region_id.clone())) + .expect("Market must be in visit count map"); + + // Break if max iterations reached for this market + if *count == MAX_ITERATIONS { + warn!("Max cycle investment iterations reached. Unable to find stable solution"); + break; + } + + // Select/retain assets for this market + // This will replace any previously selected/retained assets for this market + let assets = select_assets_for_single_market( + model, + &commodity_id, + ®ion_id, + year, + ¤t_demand, + existing_assets, + prices, + writer, + )?; + assets_for_cycle.insert((commodity_id.clone(), region_id.clone()), assets.clone()); + + // If no assets were selected/retained, no need to run dispatch + if assets.is_empty() { + continue; + } + + // Assemble full list of assets for dispatch + // This is all previously selected assets plus assets selected/retained as part of the cycle + let mut all_assets = previously_selected_assets.to_vec(); + let assets_for_cycle_flat: Vec<_> = assets_for_cycle + .values() + .flat_map(|v| v.iter().cloned()) + .collect(); + all_assets.extend_from_slice(&assets_for_cycle_flat); + + // We allow all `Selected` state assets to have flexible capacity + let flexible_capacity_assets: Vec<_> = assets_for_cycle_flat + .iter() + .filter(|asset| matches!(asset.state(), AssetState::Selected { .. })) + .cloned() + .collect(); + + // We allow unmet demand for all markets except the one we're currently balancing + let markets_to_allow_unmet_demand: Vec<_> = markets + .iter() + .filter(|m| *m != &(commodity_id.clone(), region_id.clone())) + .cloned() + .collect(); + + // Run dispatch + let solution = DispatchRun::new(model, &all_assets, year) + .with_market_balance_subset(&markets_to_balance) + .with_unmet_demand_vars(&markets_to_allow_unmet_demand) + .with_flexible_capacity_assets(&flexible_capacity_assets, CAPACITY_MARGIN) + .run( + &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment (round {count})"), + writer, + )?; + + // Find markets with unmet demand and add these to the queue if not already present + let coms = solution.get_markets_with_unmet_demand(); + for market in coms { + if !unbalanced_markets.contains(&market) { + unbalanced_markets.push_back(market); + } + } + + // Calculate new demand map including unmet demand from this solution + current_demand.clone_from(demand); + update_demand_map_with_unmet_demand(&mut current_demand, solution.iter_unmet_demand()); + last_solution = Some(solution); + *count += 1; + } + + // Finally, update flexible capacity assets based on the final solution + let mut all_cycle_assets: Vec<_> = assets_for_cycle.into_values().flatten().collect(); + if let Some(solution) = last_solution { + let new_capacities: HashMap<_, _> = solution.iter_capacity().collect(); + for asset in &mut all_cycle_assets { + if let Some(new_capacity) = new_capacities.get(asset) { + asset.make_mut().set_capacity(*new_capacity); + } + } + } + + Ok(all_cycle_assets) +} + /// Flatten the preset commodity demands for a given year into a map of commodity, region and /// time slice to demand. /// @@ -352,6 +484,25 @@ fn update_net_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: &[A } } +/// Update demand map with unmmet demand from dispatch solution +fn update_demand_map_with_unmet_demand<'a, I>(demand: &mut AllDemandMap, unmet_demand: I) +where + I: Iterator, +{ + for (c_id, r_id, time_slice, flow) in unmet_demand { + // Only consider positive unmet demand + if flow <= Flow(0.0) { + continue; + } + + let key = (c_id.clone(), r_id.clone(), time_slice.clone()); + demand + .entry(key) + .and_modify(|value| *value += flow) + .or_insert(flow); + } +} + /// Get a portion of the demand profile for this market fn get_demand_portion_for_market( time_slice_info: &TimeSliceInfo, diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 4e2a3b5c5..b32531b08 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -1,7 +1,7 @@ //! Code for performing dispatch optimisation. //! //! This is used to calculate commodity flows and prices. -use crate::asset::{Asset, AssetRef}; +use crate::asset::{Asset, AssetRef, AssetState}; use crate::commodity::CommodityID; use crate::input::format_items_with_cap; use crate::model::Model; @@ -9,19 +9,19 @@ use crate::output::DataWriter; use crate::region::RegionID; use crate::simulation::CommodityPrices; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::{Activity, Flow, Money, MoneyPerActivity, MoneyPerFlow, UnitType}; -use anyhow::{Result, bail, ensure}; +use crate::units::{Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerFlow}; +use anyhow::{Result, bail}; use highs::{HighsModelStatus, HighsStatus, RowProblem as Problem, Sense}; use indexmap::{IndexMap, IndexSet}; use itertools::{chain, iproduct}; -use log::debug; +use log::warn; use std::collections::HashSet; use std::error::Error; use std::fmt; use std::ops::Range; mod constraints; -use constraints::{ConstraintKeys, add_asset_constraints}; +use constraints::{ConstraintKeys, add_model_constraints}; /// A map of commodity flows calculated during the optimisation pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), Flow>; @@ -32,8 +32,11 @@ pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), Flow>; /// particular column of the problem. type Variable = highs::Col; -/// The map of variables related to assets -type AssetVariableMap = IndexMap<(AssetRef, TimeSliceID), Variable>; +/// The map of activity variables for assets +type ActivityVariableMap = IndexMap<(AssetRef, TimeSliceID), Variable>; + +/// A map of capacity variables for assets +type CapacityVariableMap = IndexMap; /// Variables representing unmet demand for a given market type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Variable>; @@ -48,14 +51,16 @@ type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Var /// 2. To keep track of the combination of parameters that each variable corresponds to, for when we /// are reading the results of the optimisation. pub struct VariableMap { - asset_vars: AssetVariableMap, + activity_vars: ActivityVariableMap, existing_asset_var_idx: Range, + capacity_vars: CapacityVariableMap, + capacity_var_idx: Range, unmet_demand_vars: UnmetDemandVariableMap, unmet_demand_var_idx: Range, } impl VariableMap { - /// Create a new [`VariableMap`] and add variables to the problem + /// Create a new [`VariableMap`] and add activity variables to the problem /// /// # Arguments /// @@ -65,7 +70,7 @@ impl VariableMap { /// * `existing_assets` - The asset pool /// * `candidate_assets` - Candidate assets for inclusion in active pool /// * `year` - Current milestone year - fn new_with_asset_vars( + fn new_with_activity_vars( problem: &mut Problem, model: &Model, input_prices: Option<&CommodityPrices>, @@ -73,18 +78,18 @@ impl VariableMap { candidate_assets: &[AssetRef], year: u32, ) -> Self { - let mut asset_vars = AssetVariableMap::new(); - let existing_asset_var_idx = add_asset_variables( + let mut activity_vars = ActivityVariableMap::new(); + let existing_asset_var_idx = add_activity_variables( problem, - &mut asset_vars, + &mut activity_vars, &model.time_slice_info, input_prices, existing_assets, year, ); - add_asset_variables( + add_activity_variables( problem, - &mut asset_vars, + &mut activity_vars, &model.time_slice_info, input_prices, candidate_assets, @@ -92,8 +97,10 @@ impl VariableMap { ); Self { - asset_vars, + activity_vars, existing_asset_var_idx, + capacity_vars: CapacityVariableMap::new(), + capacity_var_idx: Range::default(), unmet_demand_vars: UnmetDemandVariableMap::default(), unmet_demand_var_idx: Range::default(), } @@ -105,14 +112,14 @@ impl VariableMap { /// /// * `problem` - The optimisation problem /// * `model` - The model - /// * `markets` - The subset of markets the problem is being run for + /// * `markets_to_allow_unmet_demand` - The subset of markets to assign unmet demand variables to fn add_unmet_demand_variables( &mut self, problem: &mut Problem, model: &Model, - markets: &[(CommodityID, RegionID)], + markets_to_allow_unmet_demand: &[(CommodityID, RegionID)], ) { - assert!(!markets.is_empty()); + assert!(!markets_to_allow_unmet_demand.is_empty()); // This line **must** come before we add more variables let start = problem.num_cols(); @@ -120,24 +127,26 @@ impl VariableMap { // Add variables let voll = model.parameters.value_of_lost_load; self.unmet_demand_vars.extend( - iproduct!(markets.iter(), model.time_slice_info.iter_ids()).map( - |((commodity_id, region_id), time_slice)| { - let key = (commodity_id.clone(), region_id.clone(), time_slice.clone()); - let var = problem.add_column(voll.value(), 0.0..); - (key, var) - }, - ), + iproduct!( + markets_to_allow_unmet_demand.iter(), + model.time_slice_info.iter_ids() + ) + .map(|((commodity_id, region_id), time_slice)| { + let key = (commodity_id.clone(), region_id.clone(), time_slice.clone()); + let var = problem.add_column(voll.value(), 0.0..); + (key, var) + }), ); self.unmet_demand_var_idx = start..problem.num_cols(); } - /// Get the asset [`Variable`] corresponding to the given parameters. - fn get_asset_var(&self, asset: &AssetRef, time_slice: &TimeSliceID) -> Variable { + /// Get the activity [`Variable`] corresponding to the given parameters. + fn get_activity_var(&self, asset: &AssetRef, time_slice: &TimeSliceID) -> Variable { let key = (asset.clone(), time_slice.clone()); *self - .asset_vars + .activity_vars .get(&key) .expect("No asset variable found for given params") } @@ -155,16 +164,21 @@ impl VariableMap { .expect("No unmet demand variable for given params") } - /// Iterate over the asset variables - fn iter_asset_vars(&self) -> impl Iterator { - self.asset_vars + /// Iterate over the activity variables + fn iter_activity_vars(&self) -> impl Iterator { + self.activity_vars .iter() .map(|((asset, time_slice), var)| (asset, time_slice, *var)) } - /// Iterate over the keys for asset variables - fn asset_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> { - self.asset_vars.keys() + /// Iterate over the keys for activity variables + fn activity_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> { + self.activity_vars.keys() + } + + /// Iterate over capacity variables + fn iter_capacity_vars(&self) -> impl Iterator { + self.capacity_vars.iter().map(|(asset, var)| (asset, *var)) } } @@ -202,7 +216,7 @@ impl Solution<'_> { /// Activity for each existing asset pub fn iter_activity(&self) -> impl Iterator { self.variables - .asset_var_keys() + .activity_var_keys() .zip(self.solution.columns()) .map(|((asset, time_slice), activity)| (asset, time_slice, Activity(*activity))) } @@ -211,10 +225,12 @@ impl Solution<'_> { fn iter_activity_for_existing( &self, ) -> impl Iterator { - self.zip_var_keys_with_output( - &self.variables.existing_asset_var_idx, - self.solution.columns(), - ) + let cols = &self.solution.columns()[self.variables.existing_asset_var_idx.clone()]; + self.variables + .activity_var_keys() + .skip(self.variables.existing_asset_var_idx.start) + .zip(cols.iter()) + .map(|((asset, time_slice), &value)| (asset, time_slice, Activity(value))) } /// Iterate over unmet demand @@ -230,6 +246,15 @@ impl Solution<'_> { }) } + /// Iterate over capacity values + pub fn iter_capacity(&self) -> impl Iterator { + self.variables + .capacity_vars + .keys() + .zip(self.solution.columns()[self.variables.capacity_var_idx.clone()].iter()) + .map(|(asset, capacity)| (asset, Capacity(*capacity))) + } + /// Keys and dual values for commodity balance constraints. pub fn iter_commodity_balance_duals( &self, @@ -248,9 +273,19 @@ impl Solution<'_> { } /// Keys and dual values for activity constraints. + /// + /// Note: if there are any flexible capacity assets, these will have two duals with identical + /// keys, and there will be no way to distinguish between them in the resulting iterator. + /// Recommended for now only to use this function when there are no flexible capacity assets. pub fn iter_activity_duals( &self, ) -> impl Iterator { + if self.iter_capacity().next().is_some() { + warn!( + "Warning: iter_activity_duals called on solution of model with flexible capacity assets. + The resulting duals may be ambiguous." + ); + } self.constraint_keys .activity_keys .zip_duals(self.solution.dual_rows()) @@ -262,27 +297,17 @@ impl Solution<'_> { &self, ) -> impl Iterator { self.variables - .asset_var_keys() + .activity_var_keys() .zip(self.solution.dual_columns()) .map(|((asset, time_slice), dual)| (asset, time_slice, MoneyPerActivity(*dual))) } - /// Zip a subset of keys in the variable map with a subset of the given output variable. - /// - /// # Arguments - /// - /// * `variable_idx` - The subset of variables to look at - /// * `output` - The output variable of interest - fn zip_var_keys_with_output<'a, T: UnitType>( - &'a self, - variable_idx: &Range, - output: &'a [f64], - ) -> impl Iterator + use<'a, T> { - let keys = self.variables.asset_var_keys().skip(variable_idx.start); - assert!(keys.len() >= variable_idx.len()); - - keys.zip(output[variable_idx.clone()].iter()) - .map(|((asset, time_slice), value)| (asset, time_slice, T::new(*value))) + /// Get the markets for which unmet demand is positive in any time slice. + pub fn get_markets_with_unmet_demand(&self) -> IndexSet<(CommodityID, RegionID)> { + self.iter_unmet_demand() + .filter(|(_, _, _, flow)| *flow > Flow(0.0)) + .map(|(commodity_id, region_id, _, _)| (commodity_id.clone(), region_id.clone())) + .collect() } } @@ -320,19 +345,19 @@ pub fn solve_optimal(model: highs::Model) -> Result = markets.iter().collect(); - let has_prices_for_market_subset = input_prices.keys().any(|(commodity_id, region_id, _)| { - markets_set.contains(&(commodity_id.clone(), region_id.clone())) - }); - assert!( - !has_prices_for_market_subset, - "Input prices were included for commodities that are being modelled, which is not allowed." - ); +/// Select prices for commodities not being balanced +fn select_input_prices( + input_prices: &CommodityPrices, + markets_to_balance: &[(CommodityID, RegionID)], +) -> CommodityPrices { + let commodity_regions_set: HashSet<(&CommodityID, &RegionID)> = + markets_to_balance.iter().map(|m| (&m.0, &m.1)).collect(); + input_prices + .iter() + .filter(|(commodity_id, region_id, _, _)| { + !commodity_regions_set.contains(&(commodity_id, region_id)) + }) + .collect() } /// Provides the interface for running the dispatch optimisation. @@ -346,10 +371,13 @@ fn check_input_prices(input_prices: &CommodityPrices, markets: &[(CommodityID, R pub struct DispatchRun<'model, 'run> { model: &'model Model, existing_assets: &'run [AssetRef], + flexible_capacity_assets: &'run [AssetRef], candidate_assets: &'run [AssetRef], markets_to_balance: &'run [(CommodityID, RegionID)], + markets_to_allow_unmet_demand: &'run [(CommodityID, RegionID)], input_prices: Option<&'run CommodityPrices>, year: u32, + capacity_margin: f64, } impl<'model, 'run> DispatchRun<'model, 'run> { @@ -358,10 +386,26 @@ impl<'model, 'run> DispatchRun<'model, 'run> { Self { model, existing_assets: assets, + flexible_capacity_assets: &[], candidate_assets: &[], markets_to_balance: &[], + markets_to_allow_unmet_demand: &[], input_prices: None, year, + capacity_margin: 0.0, + } + } + + /// Include the specified flexible capacity assets in the dispatch run + pub fn with_flexible_capacity_assets( + self, + flexible_capacity_assets: &'run [AssetRef], + capacity_margin: f64, + ) -> Self { + Self { + flexible_capacity_assets, + capacity_margin, + ..self } } @@ -374,7 +418,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { } /// Only apply commodity balance constraints to the specified subset of markets - pub fn with_market_subset(self, markets: &'run [(CommodityID, RegionID)]) -> Self { + pub fn with_market_balance_subset(self, markets: &'run [(CommodityID, RegionID)]) -> Self { assert!(!markets.is_empty()); Self { @@ -391,6 +435,14 @@ impl<'model, 'run> DispatchRun<'model, 'run> { } } + /// Allow unmet demand variables for the specified subset of markets + pub fn with_unmet_demand_vars(self, markets: &'run [(CommodityID, RegionID)]) -> Self { + Self { + markets_to_allow_unmet_demand: markets, + ..self + } + } + /// Perform the dispatch optimisation. /// /// # Arguments @@ -403,15 +455,6 @@ impl<'model, 'run> DispatchRun<'model, 'run> { /// A solution containing new commodity flows for assets and prices for (some) commodities or an /// error. pub fn run(self, run_description: &str, writer: &mut DataWriter) -> Result> { - let solution = self.run_no_save()?; - writer.write_dispatch_debug_info(self.year, run_description, &solution)?; - Ok(solution) - } - - /// Run dispatch without saving the results. - /// - /// This is an internal function as callers always want to save results. - fn run_no_save(&self) -> Result> { // If the user provided no markets to balance, we all use of them let all_markets: Vec<_>; let markets_to_balance = if self.markets_to_balance.is_empty() { @@ -420,103 +463,120 @@ impl<'model, 'run> DispatchRun<'model, 'run> { } else { self.markets_to_balance }; - if let Some(input_prices) = self.input_prices { - check_input_prices(input_prices, markets_to_balance); - } + + // Select prices for commodities not being balanced + let input_prices_owned = self + .input_prices + .map(|prices| select_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 regions/commodities to users - match self.run_without_unmet_demand(markets_to_balance) { - Ok(solution) => Ok(solution), + // is due to unmet demand, in this case, we rerun dispatch including with unmet demand + // variables for all markets so we can report the offending markets to users + match self.run_internal( + markets_to_balance, + self.markets_to_allow_unmet_demand, + input_prices, + ) { + Ok(solution) => { + writer.write_dispatch_debug_info(self.year, run_description, &solution)?; + Ok(solution) + } Err(ModelError::NonOptimal(HighsModelStatus::Infeasible)) => { - let markets = self - .get_markets_with_unmet_demand(markets_to_balance) - .expect("Failed to run dispatch to calculate unmet demand"); - - ensure!( - !markets.is_empty(), - "Model is infeasible, but there was no unmet demand" - ); - + let solution = + self.run_internal(markets_to_balance, markets_to_balance, input_prices)?; + writer.write_dispatch_debug_info(self.year, run_description, &solution)?; + let markets_with_unmet_demand: IndexSet = solution + .iter_unmet_demand() + .filter(|(_, _, _, flow)| *flow > Flow(0.0)) + .map(|(commodity_id, region_id, _, _)| format!("{commodity_id}|{region_id}")) + .collect(); 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 region and commodity pairs: {}", - format_items_with_cap(markets) + for the following markets: {}", + format_items_with_cap(markets_with_unmet_demand) ) } Err(err) => Err(err)?, } } - /// Run dispatch without unmet demand variables - fn run_without_unmet_demand( - &self, - markets: &[(CommodityID, RegionID)], - ) -> Result, ModelError> { - self.run_internal(markets, /*allow_unmet_demand=*/ false) - } - - /// Run dispatch to diagnose which markets have unmet demand - fn get_markets_with_unmet_demand( - &self, - markets: &[(CommodityID, RegionID)], - ) -> Result> { - let solution = self.run_internal(markets, /*allow_unmet_demand=*/ true)?; - Ok(solution - .iter_unmet_demand() - .filter(|(_, _, _, flow)| *flow > Flow(0.0)) - .map(|(commodity_id, region_id, _, _)| (commodity_id.clone(), region_id.clone())) - .collect()) - } - - /// Run dispatch for specified commodities, optionally including unmet demand variables + /// Run dispatch to balance the specified markets, allowing unmet demand for a subset of these fn run_internal( &self, - markets: &[(CommodityID, RegionID)], - allow_unmet_demand: bool, + markets_to_balance: &[(CommodityID, RegionID)], + markets_to_allow_unmet_demand: &[(CommodityID, RegionID)], + input_prices: Option<&CommodityPrices>, ) -> Result, ModelError> { // Set up problem let mut problem = Problem::default(); - let mut variables = VariableMap::new_with_asset_vars( + let mut variables = VariableMap::new_with_activity_vars( &mut problem, self.model, - self.input_prices, + input_prices, self.existing_assets, self.candidate_assets, self.year, ); - // If unmet demand is enabled for this dispatch run (and is allowed by the model param) then - // we add variables representing unmet demand - if allow_unmet_demand { - variables.add_unmet_demand_variables(&mut problem, self.model, markets); + // Check flexible capacity assets is a subset of existing assets + for asset in self.flexible_capacity_assets { + assert!( + self.existing_assets.contains(asset), + "Flexible capacity assets must be a subset of existing assets. Offending asset: {asset:?}" + ); + } + + // Add capacity variables for flexible capacity assets + if !self.flexible_capacity_assets.is_empty() { + variables.capacity_var_idx = add_capacity_variables( + &mut problem, + &mut variables.capacity_vars, + self.flexible_capacity_assets, + self.capacity_margin, + ); + } + + // Check that markets_to_allow_unmet_demand is a subset of markets_to_balance + for market in markets_to_allow_unmet_demand { + assert!( + markets_to_balance.contains(market), + "markets_to_allow_unmet_demand must be a subset of markets_to_balance. \ + Offending market: {market:?}" + ); + } + + // Add variables representing unmet demand + if !markets_to_allow_unmet_demand.is_empty() { + variables.add_unmet_demand_variables( + &mut problem, + self.model, + markets_to_allow_unmet_demand, + ); } // Add constraints let all_assets = chain(self.existing_assets.iter(), self.candidate_assets.iter()); - let constraint_keys = add_asset_constraints( + let constraint_keys = add_model_constraints( &mut problem, &variables, self.model, &all_assets, - markets, + markets_to_balance, + markets_to_allow_unmet_demand, self.year, ); // Solve model let solution = solve_optimal(problem.optimise(Sense::Minimise))?; - let objective_value = Money(solution.objective_value()); - debug!("Objective value: {objective_value}"); - Ok(Solution { solution: solution.get_solution(), variables, time_slice_info: &self.model.time_slice_info, constraint_keys, - objective_value, + objective_value: Money(solution.objective_value()), }) } } @@ -531,9 +591,9 @@ impl<'model, 'run> DispatchRun<'model, 'run> { /// * `input_prices` - Optional explicit prices for input commodities /// * `assets` - Assets to include /// * `year` - Current milestone year -fn add_asset_variables( +fn add_activity_variables( problem: &mut Problem, - variables: &mut AssetVariableMap, + variables: &mut ActivityVariableMap, time_slice_info: &TimeSliceInfo, input_prices: Option<&CommodityPrices>, assets: &[AssetRef], @@ -553,6 +613,37 @@ fn add_asset_variables( start..problem.num_cols() } +fn add_capacity_variables( + problem: &mut Problem, + variables: &mut CapacityVariableMap, + assets: &[AssetRef], + capacity_margin: f64, +) -> Range { + // This line **must** come before we add more variables + let start = problem.num_cols(); + + for asset in assets { + // Can only have flexible capacity for Candidate or Selected assets + if !matches!( + asset.state(), + AssetState::Candidate | AssetState::Selected { .. } + ) { + panic!( + "Flexible capacity can only be assigned to Candidate or Selected assets. Offending asset: {asset:?}" + ); + } + + let current_capacity = asset.capacity().value(); + let bounds = + (1.0 - capacity_margin) * current_capacity..=(1.0 + capacity_margin) * current_capacity; + let var = problem.add_column(0.0, bounds); + let existing = variables.insert(asset.clone(), var).is_some(); + assert!(!existing, "Duplicate entry for var"); + } + + start..problem.num_cols() +} + /// Calculate the cost coefficient for a decision variable. /// /// Normally, the cost coefficient is the same as the asset's operating costs for the given year and diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs index d0a229009..f5484c2e0 100644 --- a/src/simulation/optimisation/constraints.rs +++ b/src/simulation/optimisation/constraints.rs @@ -7,6 +7,7 @@ use crate::region::RegionID; use crate::time_slice::{TimeSliceID, TimeSliceSelection}; use crate::units::UnitType; use highs::RowProblem as Problem; +use indexmap::IndexMap; /// Corresponding variables for a constraint along with the row offset in the solution pub struct KeysWithOffset { @@ -45,7 +46,7 @@ pub struct ConstraintKeys { pub activity_keys: ActivityKeys, } -/// Add asset-level constraints +/// Add constraints for the dispatch model. /// /// Note: the ordering of constraints is important, as the dual values of the constraints must later /// be retrieved to calculate commodity prices. @@ -56,25 +57,34 @@ pub struct ConstraintKeys { /// * `variables` - The variables in the problem /// * `model` - The model /// * `assets` - The asset pool -/// * `markets` - The subset of markets to apply constraints to +/// * `markets_to_balance` - The subset of markets to apply balance constraints to +/// * `markets_to_allow_unmet_demand` - The subset of markets to assign unmet demand variables to /// * `year` - Current milestone year /// /// # Returns /// /// Keys for the different constraints. -pub fn add_asset_constraints<'a, I>( +pub fn add_model_constraints<'a, I>( problem: &mut Problem, variables: &VariableMap, model: &'a Model, assets: &I, - markets: &'a [(CommodityID, RegionID)], + markets_to_balance: &'a [(CommodityID, RegionID)], + markets_to_allow_unmet_demand: &'a [(CommodityID, RegionID)], year: u32, ) -> ConstraintKeys where I: Iterator + Clone + 'a, { - let commodity_balance_keys = - add_commodity_balance_constraints(problem, variables, model, assets, markets, year); + let commodity_balance_keys = add_commodity_balance_constraints( + problem, + variables, + model, + assets, + markets_to_balance, + markets_to_allow_unmet_demand, + year, + ); let activity_keys = add_activity_constraints(problem, variables); @@ -97,7 +107,8 @@ fn add_commodity_balance_constraints<'a, I>( variables: &VariableMap, model: &'a Model, assets: &I, - markets: &'a [(CommodityID, RegionID)], + markets_to_balance: &'a [(CommodityID, RegionID)], + markets_to_allow_unmet_demand: &'a [(CommodityID, RegionID)], year: u32, ) -> CommodityBalanceKeys where @@ -108,7 +119,7 @@ where let mut keys = Vec::new(); let mut terms = Vec::new(); - for (commodity_id, region_id) in markets { + for (commodity_id, region_id) in markets_to_balance { let commodity = &model.commodities[commodity_id]; if !matches!( commodity.kind, @@ -129,7 +140,7 @@ where // If the commodity has a time slice level of season/annual, the constraint will // cover multiple time slices for (time_slice, _) in ts_selection.iter(&model.time_slice_info) { - let var = variables.get_asset_var(asset, time_slice); + let var = variables.get_activity_var(asset, time_slice); terms.push((var, flow.coeff.value())); } } @@ -142,7 +153,7 @@ where } // Also include unmet demand variables if required - if !variables.unmet_demand_var_idx.is_empty() { + if markets_to_allow_unmet_demand.contains(&(commodity_id.clone(), region_id.clone())) { for (time_slice, _) in ts_selection.iter(&model.time_slice_info) { let var = variables.get_unmet_demand_var(commodity_id, region_id, time_slice); terms.push((var, 1.0)); @@ -185,12 +196,34 @@ fn add_activity_constraints(problem: &mut Problem, variables: &VariableMap) -> A let offset = problem.num_rows(); let mut keys = Vec::new(); - for (asset, time_slice, var) in variables.iter_asset_vars() { - let limits = asset.get_activity_limits(time_slice); - let limits = limits.start().value()..=limits.end().value(); - - problem.add_row(limits, [(var, 1.0)]); - keys.push((asset.clone(), time_slice.clone())); + let capacity_vars: IndexMap<&AssetRef, highs::Col> = variables.iter_capacity_vars().collect(); + for (asset, time_slice, activity_var) in variables.iter_activity_vars() { + if let Some(&capacity_var) = capacity_vars.get(asset) { + // Asset has flexible capacity. + let per_capacity_limits = asset.get_activity_per_capacity_limits(time_slice); + let lower_limit = per_capacity_limits.start().value(); + let upper_limit = per_capacity_limits.end().value(); + + // Upper bound: activity ≤ capacity * upper_limit + problem.add_row(..=0.0, [(activity_var, 1.0), (capacity_var, -upper_limit)]); + + // Lower bound: activity ≥ capacity * lower_limit + problem.add_row(..=0.0, [(activity_var, -1.0), (capacity_var, lower_limit)]); + + // Store keys for retrieving duals later. + // TODO: a bit of a hack pushing identical keys twice. Safe for now so long as we don't + // use the activity duals for anything important when using flexible capacity assets. + keys.push((asset.clone(), time_slice.clone())); + keys.push((asset.clone(), time_slice.clone())); + } else { + // Fixed-capacity asset: simple absolute activity limits. + let activity_limits = asset.get_activity_limits(time_slice); + let range = activity_limits.start().value()..=activity_limits.end().value(); + problem.add_row(range, [(activity_var, 1.0)]); + + // Store keys for retrieving duals later. + keys.push((asset.clone(), time_slice.clone())); + } } ActivityKeys { offset, keys } diff --git a/tests/data/simple/debug_unmet_demand.csv b/tests/data/simple/debug_unmet_demand.csv new file mode 100644 index 000000000..e69de29bb