diff --git a/schemas/input/commodities.yaml b/schemas/input/commodities.yaml index f5b72686e..d88ffe993 100644 --- a/schemas/input/commodities.yaml +++ b/schemas/input/commodities.yaml @@ -16,12 +16,29 @@ fields: description: A human-readable label for the commodity - name: type type: string + enum: [svd, sed, oth] description: The type of commodity notes: | Must be one of `svd` (service demand), `sed` (supply equals demand) or `oth` (other) - name: time_slice_level type: string + enum: [annual, season, daynight] description: The time slice level at which constraints for this commodity are applied notes: | Must be one of `annual` (whole year), `season` (whole season) or `daynight` (a particular time of day) + - name: pricing_strategy + type: string + enum: [shadow, marginal, full, scarcity, unpriced] + description: The pricing strategy for this commodity + notes: | + Optional. If specified, must be one of `shadow` (priced using shadow prices from supply-demand + constraints), `marginal` (priced at the marginal cost of production) , `full` (priced at the + full cost of production including capital costs), `scarcity` (priced adjusted for scarcity), + or `unpriced` (not priced at all). + + For `svd` and `sed` commodities, this must be one of `shadow`, `marginal` or `full`. For `oth` + commodities, this must be `unpriced`. + + If unspecified, the commodity will use the default pricing strategy according to the commodity + type: `shadow` for `svd` and `sed` commodities, and `unpriced` for `oth` commodities. diff --git a/schemas/input/model.yaml b/schemas/input/model.yaml index d9c51aa3b..e58784dd4 100644 --- a/schemas/input/model.yaml +++ b/schemas/input/model.yaml @@ -22,17 +22,6 @@ properties: notes: | This is the proportion of the maximum required capacity across time slices (for a given asset/commodity etc. combination). - pricing_strategy: - type: string - enum: [shadow_prices, scarcity_adjusted] - description: Change the algorithm used for calculating commodity prices - default: shadow_prices - notes: | - The `shadow_prices` option just uses the shadow prices for commodity prices. - - The `scarcity_adjusted` option adjusts prices for scarcity. This may cause price instability - for assets with more than one output commodity. Do not use this unless you know what you're - doing! value_of_lost_load: type: number description: The cost applied to unmet demand diff --git a/src/asset.rs b/src/asset.rs index 2240c947d..fe6e75665 100644 --- a/src/asset.rs +++ b/src/asset.rs @@ -1,13 +1,17 @@ //! Assets are instances of a process which are owned and invested in by agents. use crate::agent::AgentID; -use crate::commodity::CommodityID; +use crate::commodity::{CommodityID, CommodityType}; +use crate::finance::annual_capital_cost; use crate::process::{ ActivityLimits, FlowDirection, Process, ProcessFlow, ProcessID, ProcessParameter, }; use crate::region::RegionID; use crate::simulation::CommodityPrices; use crate::time_slice::{TimeSliceID, TimeSliceSelection}; -use crate::units::{Activity, ActivityPerCapacity, Capacity, MoneyPerActivity}; +use crate::units::{ + Activity, ActivityPerCapacity, Capacity, FlowPerActivity, MoneyPerActivity, MoneyPerCapacity, + MoneyPerFlow, +}; use anyhow::{Context, Result, ensure}; use indexmap::IndexMap; use itertools::{Itertools, chain}; @@ -288,6 +292,18 @@ impl Asset { (cap2act * *limits.start())..=(cap2act * *limits.end()) } + /// Get the activity limits for this asset for a given time slice selection + pub fn get_activity_limits_for_selection( + &self, + time_slice_selection: &TimeSliceSelection, + ) -> RangeInclusive { + let activity_per_capacity_limits = self.activity_limits.get_limit(time_slice_selection); + let cap2act = self.process.capacity_to_activity; + let lb = self.capacity * cap2act * *activity_per_capacity_limits.start(); + let ub = self.capacity * cap2act * *activity_per_capacity_limits.end(); + lb..=ub + } + /// Iterate over activity limits for this asset pub fn iter_activity_limits( &self, @@ -318,12 +334,22 @@ impl Asset { }) } + /// Gets the total SED/SVD output per unit of activity for this asset + /// + /// Note: Since we are summing coefficients from different commodities, this ONLY makes sense + /// if these commodities have the same units (e.g., all in PJ). Users are currently not made to + /// give units for commodities, so we cannot possibly enforce this. Something to potentially + /// address in future. + pub fn get_total_output_per_activity(&self) -> FlowPerActivity { + self.iter_output_flows().map(|flow| flow.coeff).sum() + } + /// Get the operating cost for this asset in a given year and time slice pub fn get_operating_cost(&self, year: u32, time_slice: &TimeSliceID) -> MoneyPerActivity { // The cost for all commodity flows (including levies/incentives) - let flows_cost: MoneyPerActivity = self + let flows_cost = self .iter_flows() - .map(|flow| flow.get_total_cost(&self.region_id, year, time_slice)) + .map(|flow| flow.get_total_cost_per_activity(&self.region_id, year, time_slice)) .sum(); self.process_parameter.variable_operating_cost + flows_cost @@ -355,15 +381,16 @@ impl Asset { }) } - /// Get the cost of input flows using the commodity prices in `input_prices`. + /// Get the total cost of purchasing input commodities per unit of activity for this asset. /// /// If a price is missing, there is assumed to be no cost. pub fn get_input_cost_from_prices( &self, - input_prices: &CommodityPrices, + prices: &CommodityPrices, time_slice: &TimeSliceID, ) -> MoneyPerActivity { - -self.get_revenue_from_flows_with_filter(input_prices, time_slice, |x| { + // Revenues of input flows are negative costs, so we negate the result + -self.get_revenue_from_flows_with_filter(prices, time_slice, |x| { x.direction() == FlowDirection::Input }) } @@ -386,12 +413,142 @@ impl Asset { .map(|flow| { flow.coeff * prices - .get(&flow.commodity.id, self.region_id(), time_slice) + .get(&flow.commodity.id, &self.region_id, time_slice) .unwrap_or_default() }) .sum() } + /// Get the generic activity cost per unit of activity for this asset. + /// + /// These are all activity-related costs that are not associated with specific SED/SVD outputs. + /// Includes levies, flow costs, costs of inputs and variable operating costs + fn get_generic_activity_cost( + &self, + prices: &CommodityPrices, + year: u32, + time_slice: &TimeSliceID, + ) -> MoneyPerActivity { + // The cost of purchasing input commodities + let cost_of_inputs = self.get_input_cost_from_prices(prices, time_slice); + + // Flow costs/levies for all flows except SED/SVD outputs + let excludes_sed_svd_output = |flow: &&ProcessFlow| { + !(flow.direction() == FlowDirection::Output + && matches!( + flow.commodity.kind, + CommodityType::SupplyEqualsDemand | CommodityType::ServiceDemand + )) + }; + let flow_costs = self + .iter_flows() + .filter(excludes_sed_svd_output) + .map(|flow| flow.get_total_cost_per_activity(&self.region_id, year, time_slice)) + .sum(); + + cost_of_inputs + flow_costs + self.process_parameter.variable_operating_cost + } + + /// Iterate over marginal costs for a filtered set of SED/SVD output commodities for this asset + /// + /// For each SED/SVD output commodity, the marginal cost is calculated as the sum of: + /// - Generic activity costs (variable operating costs, cost of purchasing inputs, plus all + /// levies and flow costs not associated with specific SED/SVD outputs), which are + /// shared equally over all SED/SVD outputs + /// - Production levies and flow costs for the specific SED/SVD output commodity + pub fn iter_marginal_costs_with_filter<'a>( + &'a self, + prices: &'a CommodityPrices, + year: u32, + time_slice: &'a TimeSliceID, + filter: impl Fn(&CommodityID) -> bool + 'a, + ) -> Box + 'a> { + // Iterator over SED/SVD output flows matching the filter + let mut output_flows_iter = self + .iter_output_flows() + .filter(move |flow| filter(&flow.commodity.id)) + .peekable(); + + // If there are no output flows after filtering, return an empty iterator + if output_flows_iter.peek().is_none() { + return Box::new(std::iter::empty::<(CommodityID, MoneyPerFlow)>()); + } + + // Calculate generic activity costs. + // This is all activity costs not associated with specific SED/SVD outputs, which will get + // shared equally over all SED/SVD outputs. Includes levies, flow costs, costs of inputs and + // variable operating costs + let generic_activity_cost = self.get_generic_activity_cost(prices, year, time_slice); + + // Share generic activity costs equally over all SED/SVD outputs + // We sum the output coefficients of all SED/SVD commodities to get total output, then + // divide costs by this total output to get the generic cost per unit of output. + // Note: only works if all SED/SVD outputs have the same units - not currently checked! + let total_output_per_activity = self.get_total_output_per_activity(); + assert!(total_output_per_activity > FlowPerActivity::EPSILON); // input checks should guarantee this + let generic_cost_per_flow = generic_activity_cost / total_output_per_activity; + + // Iterate over SED/SVD output flows + Box::new(output_flows_iter.map(move |flow| { + // Get the costs for this specific commodity flow + let commodity_specific_costs_per_flow = + flow.get_total_cost_per_flow(&self.region_id, year, time_slice); + + // Add these to the generic costs to get total cost for this commodity + let marginal_cost = generic_cost_per_flow + commodity_specific_costs_per_flow; + (flow.commodity.id.clone(), marginal_cost) + })) + } + + /// Iterate over marginal costs for all SED/SVD output commodities for this asset + /// + /// See `iter_marginal_costs_with_filter` for details. + pub fn iter_marginal_costs<'a>( + &'a self, + prices: &'a CommodityPrices, + year: u32, + time_slice: &'a TimeSliceID, + ) -> Box + 'a> { + self.iter_marginal_costs_with_filter(prices, year, time_slice, move |_| true) + } + + /// Get the annual capital cost per unit of capacity for this asset + pub fn get_annual_capital_cost_per_capacity(&self) -> MoneyPerCapacity { + let capital_cost = self.process_parameter.capital_cost; + let lifetime = self.process_parameter.lifetime; + let discount_rate = self.process_parameter.discount_rate; + annual_capital_cost(capital_cost, lifetime, discount_rate) + } + + /// Get the annual capital cost per unit of activity for this asset + /// + /// Total capital costs (cost per capacity * capacity) are shared equally over the year in + /// accordance with the annual activity. + pub fn get_annual_capital_cost_per_activity( + &self, + annual_activity: Activity, + ) -> MoneyPerActivity { + let annual_capital_cost_per_capacity = self.get_annual_capital_cost_per_capacity(); + let total_annual_capital_cost = annual_capital_cost_per_capacity * self.capacity(); + assert!( + annual_activity > Activity::EPSILON, + "Cannot calculate annual capital cost per activity for an asset with zero annual activity" + ); + total_annual_capital_cost / annual_activity + } + + /// Get the annual capital cost per unit of output flow for this asset + /// + /// Total capital costs (cost per capacity * capacity) are shared equally across all output + /// flows in accordance with the annual activity and total output per unit of activity. + pub fn get_annual_capital_cost_per_flow(&self, annual_activity: Activity) -> MoneyPerFlow { + let annual_capital_cost_per_activity = + self.get_annual_capital_cost_per_activity(annual_activity); + let total_output_per_activity = self.get_total_output_per_activity(); + assert!(total_output_per_activity > FlowPerActivity::EPSILON); // input checks should guarantee this + annual_capital_cost_per_activity / total_output_per_activity + } + /// Maximum activity for this asset pub fn max_activity(&self) -> Activity { self.capacity * self.process.capacity_to_activity @@ -407,6 +564,17 @@ impl Asset { self.flows.values() } + /// Iterate over the asset's output SED/SVD flows + pub fn iter_output_flows(&self) -> impl Iterator { + self.flows.values().filter(|flow| { + flow.direction() == FlowDirection::Output + && matches!( + flow.commodity.kind, + CommodityType::SupplyEqualsDemand | CommodityType::ServiceDemand + ) + }) + } + /// Get the primary output flow (if any) for this asset pub fn primary_output(&self) -> Option<&ProcessFlow> { self.process @@ -792,7 +960,7 @@ impl AssetPool { the start of the simulation", asset.process_id(), asset.commission_year, - asset.process_parameter().lifetime + asset.process_parameter.lifetime ); continue; } diff --git a/src/commodity.rs b/src/commodity.rs index 2e535fd7c..008929ed1 100644 --- a/src/commodity.rs +++ b/src/commodity.rs @@ -24,37 +24,35 @@ pub type DemandMap = HashMap<(RegionID, u32, TimeSliceSelection), Flow>; /// /// Represents a substance (e.g. CO2) or form of energy (e.g. electricity) that can be produced or /// consumed by processes. -#[derive(PartialEq, Debug, Deserialize, Clone)] +#[derive(PartialEq, Debug, Clone)] pub struct Commodity { /// Unique identifier for the commodity (e.g. "ELC") pub id: CommodityID, /// Text description of commodity (e.g. "electricity") pub description: String, /// Commodity balance type - #[serde(rename = "type")] // NB: we can't name a field type as it's a reserved keyword pub kind: CommodityType, /// The time slice level for commodity balance pub time_slice_level: TimeSliceLevel, + /// Defines the strategy used for calculating commodity prices + pub pricing_strategy: PricingStrategy, /// Production levies for this commodity for different combinations of region, year and time slice. /// /// May be empty if there are no production levies for this commodity, otherwise there must be /// entries for every combination of parameters. Note that these values can be negative, /// indicating an incentive. - #[serde(skip)] pub levies_prod: CommodityLevyMap, /// Consumption levies for this commodity for different combinations of region, year and time slice. /// /// May be empty if there are no consumption levies for this commodity, otherwise there must be /// entries for every combination of parameters. Note that these values can be negative, /// indicating an incentive. - #[serde(skip)] pub levies_cons: CommodityLevyMap, /// Demand as defined in input files. Will be empty for non-service-demand commodities. /// /// The [`TimeSliceSelection`] part of the key is always at the same [`TimeSliceLevel`] as the /// `time_slice_level` field. E.g. if the `time_slice_level` is seasonal, then there will be /// keys representing each season (and not e.g. individual time slices). - #[serde(skip)] pub demand: DemandMap, } define_id_getter! {Commodity, CommodityID} @@ -89,6 +87,26 @@ pub enum CommodityType { Other, } +/// The strategy used for calculating commodity prices +#[derive(Debug, PartialEq, Clone, Deserialize, Hash, Eq)] +pub enum PricingStrategy { + /// Take commodity prices directly from the shadow prices + #[serde(rename = "shadow")] + Shadow, + /// Adjust shadow prices for scarcity + #[serde(rename = "scarcity")] + ScarcityAdjusted, + /// Use marginal cost of highest-cost active asset producing the commodity + #[serde(rename = "marginal")] + MarginalCost, + /// Use full cost of highest-cost active asset producing the commodity + #[serde(rename = "full")] + FullCost, + /// Commodities that should not have prices calculated + #[serde(rename = "unpriced")] + Unpriced, +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/fixture.rs b/src/fixture.rs index 9aa67f894..03ae05ef4 100644 --- a/src/fixture.rs +++ b/src/fixture.rs @@ -5,7 +5,9 @@ use crate::agent::{ DecisionRule, }; use crate::asset::{Asset, AssetPool, AssetRef}; -use crate::commodity::{Commodity, CommodityID, CommodityLevyMap, CommodityType, DemandMap}; +use crate::commodity::{ + Commodity, CommodityID, CommodityLevyMap, CommodityType, DemandMap, PricingStrategy, +}; use crate::process::{ ActivityLimits, Process, ProcessActivityLimitsMap, ProcessFlow, ProcessFlowsMap, ProcessInvestmentConstraintsMap, ProcessMap, ProcessParameter, ProcessParameterMap, @@ -70,6 +72,7 @@ pub fn svd_commodity() -> Commodity { description: "".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::DayNight, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -83,6 +86,7 @@ pub fn sed_commodity() -> Commodity { description: "Test SED commodity".into(), kind: CommodityType::SupplyEqualsDemand, time_slice_level: TimeSliceLevel::DayNight, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -96,6 +100,7 @@ pub fn other_commodity() -> Commodity { description: "Test other commodity".into(), kind: CommodityType::Other, time_slice_level: TimeSliceLevel::DayNight, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), diff --git a/src/input/agent/commodity_portion.rs b/src/input/agent/commodity_portion.rs index 8ba0ca950..10dbba414 100644 --- a/src/input/agent/commodity_portion.rs +++ b/src/input/agent/commodity_portion.rs @@ -188,7 +188,9 @@ fn validate_agent_commodity_portions( mod tests { use super::*; use crate::agent::{Agent, AgentObjectiveMap, AgentSearchSpaceMap, DecisionRule}; - use crate::commodity::{Commodity, CommodityID, CommodityLevyMap, CommodityType, DemandMap}; + use crate::commodity::{ + Commodity, CommodityID, CommodityLevyMap, CommodityType, DemandMap, PricingStrategy, + }; use crate::time_slice::TimeSliceLevel; use indexmap::IndexMap; use std::rc::Rc; @@ -216,6 +218,7 @@ mod tests { description: "A commodity".into(), kind: CommodityType::SupplyEqualsDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -260,6 +263,7 @@ mod tests { description: "Another commodity".into(), kind: CommodityType::SupplyEqualsDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), diff --git a/src/input/commodity.rs b/src/input/commodity.rs index 3d95a0ff9..3adf647fe 100644 --- a/src/input/commodity.rs +++ b/src/input/commodity.rs @@ -1,10 +1,17 @@ //! Code for reading in commodity-related data from CSV files. -use super::read_csv_id_file; -use crate::commodity::{BalanceType, Commodity, CommodityID, CommodityMap}; +use super::{input_err_msg, read_csv}; +use crate::ISSUES_URL; +use crate::commodity::{ + BalanceType, Commodity, CommodityID, CommodityLevyMap, CommodityMap, CommodityType, DemandMap, + PricingStrategy, +}; +use crate::model::{ALLOW_BROKEN_OPTION_NAME, broken_model_options_allowed}; use crate::region::RegionID; -use crate::time_slice::TimeSliceInfo; -use anyhow::Result; -use indexmap::IndexSet; +use crate::time_slice::{TimeSliceInfo, TimeSliceLevel}; +use anyhow::{Context, Ok, Result, ensure}; +use indexmap::{IndexMap, IndexSet}; +use log::warn; +use serde::Deserialize; use std::path::Path; mod levy; @@ -15,6 +22,16 @@ mod demand_slicing; const COMMODITY_FILE_NAME: &str = "commodities.csv"; +#[derive(PartialEq, Debug, Deserialize)] +struct CommodityRaw { + pub id: CommodityID, + pub description: String, + #[serde(rename = "type")] // NB: we can't name a field type as it's a reserved keyword + pub kind: CommodityType, + pub time_slice_level: TimeSliceLevel, + pub pricing_strategy: Option, +} + /// Read commodity data from the specified model directory. /// /// # Arguments @@ -33,9 +50,11 @@ pub fn read_commodities( time_slice_info: &TimeSliceInfo, milestone_years: &[u32], ) -> Result { - let commodities = - read_csv_id_file::(&model_dir.join(COMMODITY_FILE_NAME))?; + // Read commodities table + let commodities = read_commodities_file(model_dir)?; let commodity_ids = commodities.keys().cloned().collect(); + + // Read costs table let mut costs = read_commodity_levies( model_dir, &commodity_ids, @@ -44,6 +63,7 @@ pub fn read_commodities( milestone_years, )?; + // Read demand table let mut demand = read_demand( model_dir, &commodities, @@ -72,3 +92,139 @@ pub fn read_commodities( }) .collect()) } + +fn read_commodities_file(model_dir: &Path) -> Result> { + let file_path = model_dir.join(COMMODITY_FILE_NAME); + let commodities_csv = read_csv(&file_path)?; + read_commodities_file_from_iter(commodities_csv).with_context(|| input_err_msg(&file_path)) +} + +fn read_commodities_file_from_iter(iter: I) -> Result> +where + I: Iterator, +{ + let mut commodities = IndexMap::new(); + for commodity_raw in iter { + let pricing_strategy = match commodity_raw.pricing_strategy { + Some(strategy) => strategy, + None => default_pricing_strategy(&commodity_raw.kind), + }; + + let commodity = Commodity { + id: commodity_raw.id.clone(), + description: commodity_raw.description, + kind: commodity_raw.kind, + time_slice_level: commodity_raw.time_slice_level, + pricing_strategy, + levies_prod: CommodityLevyMap::default(), + levies_cons: CommodityLevyMap::default(), + demand: DemandMap::default(), + }; + + validate_commodity(&commodity)?; + + ensure!( + commodities.insert(commodity_raw.id, commodity).is_none(), + "Duplicate commodity ID" + ); + } + + Ok(commodities) +} + +/// Get the default pricing strategy for a given commodity kind. +fn default_pricing_strategy(commodity_kind: &CommodityType) -> PricingStrategy { + match commodity_kind { + CommodityType::Other => PricingStrategy::Unpriced, + CommodityType::SupplyEqualsDemand | CommodityType::ServiceDemand => PricingStrategy::Shadow, + } +} + +fn validate_commodity(commodity: &Commodity) -> Result<()> { + // Check that the pricing strategy is appropriate for the commodity type + match commodity.kind { + CommodityType::Other => { + ensure!( + commodity.pricing_strategy == PricingStrategy::Unpriced, + "Commodity {} of type Other must be unpriced. \ + Update its pricing strategy to 'unpriced' or 'default'.", + commodity.id + ); + } + CommodityType::SupplyEqualsDemand | CommodityType::ServiceDemand => { + ensure!( + commodity.pricing_strategy != PricingStrategy::Unpriced, + "Commodity {} of type {:?} cannot be unpriced. \ + Update its pricing strategy to a valid option.", + commodity.id, + commodity.kind + ); + } + } + + // Gatekeep alternative pricing options + if !matches!( + commodity.pricing_strategy, + PricingStrategy::Shadow | PricingStrategy::Unpriced + ) { + ensure!( + broken_model_options_allowed(), + "Price strategies other than 'shadow' and 'unpriced' are currently experimental. \ + To run anyway, set the {ALLOW_BROKEN_OPTION_NAME} option to true." + ); + } + if commodity.pricing_strategy == PricingStrategy::ScarcityAdjusted { + warn!( + "The pricing strategy for {} is set to 'scarcity'. Commodity prices may be \ + incorrect if assets have more than one output commodity. See: {ISSUES_URL}/677", + commodity.id + ); + } + + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::fixture::assert_error; + use crate::time_slice::TimeSliceLevel; + + fn make_commodity(kind: CommodityType, pricing_strategy: PricingStrategy) -> Commodity { + Commodity { + id: "ELC".into(), + description: "test".into(), + kind, + time_slice_level: TimeSliceLevel::Annual, + pricing_strategy, + levies_prod: CommodityLevyMap::default(), + levies_cons: CommodityLevyMap::default(), + demand: DemandMap::default(), + } + } + + #[test] + fn test_validate_commodity() { + let commodity = make_commodity(CommodityType::SupplyEqualsDemand, PricingStrategy::Shadow); + assert!(validate_commodity(&commodity).is_ok()); + } + + #[test] + fn test_validate_commodity_other_priced() { + let mut commodity = make_commodity(CommodityType::Other, PricingStrategy::MarginalCost); + assert_error!( + validate_commodity(&mut commodity), + "Commodity ELC of type Other must be unpriced. Update its pricing strategy to 'unpriced' or 'default'." + ); + } + + #[test] + fn test_validate_commodity_sed_unpriced() { + let mut commodity = + make_commodity(CommodityType::SupplyEqualsDemand, PricingStrategy::Unpriced); + assert_error!( + validate_commodity(&mut commodity), + "Commodity ELC of type SupplyEqualsDemand cannot be unpriced. Update its pricing strategy to a valid option." + ); + } +} diff --git a/src/model.rs b/src/model.rs index e996ff6e7..14183d2fa 100644 --- a/src/model.rs +++ b/src/model.rs @@ -9,9 +9,7 @@ use std::collections::HashMap; use std::path::PathBuf; pub mod parameters; -pub use parameters::{ - ALLOW_BROKEN_OPTION_NAME, ModelParameters, PricingStrategy, broken_model_options_allowed, -}; +pub use parameters::{ALLOW_BROKEN_OPTION_NAME, ModelParameters, broken_model_options_allowed}; /// Model definition pub struct Model { diff --git a/src/model/parameters.rs b/src/model/parameters.rs index 0311a5c60..f0a5d12d8 100644 --- a/src/model/parameters.rs +++ b/src/model/parameters.rs @@ -1,5 +1,4 @@ //! Defines the `ModelParameters` struct, which represents the contents of `model.toml`. -use crate::ISSUES_URL; use crate::asset::check_capacity_valid_for_asset; use crate::input::{ deserialise_proportion_nonzero, input_err_msg, is_sorted_and_unique, read_toml, @@ -8,7 +7,6 @@ use crate::units::{Capacity, Dimensionless, MoneyPerFlow}; use anyhow::{Context, Result, ensure}; use log::warn; use serde::Deserialize; -use serde_string_enum::DeserializeLabeledStringEnum; use std::path::Path; use std::sync::OnceLock; @@ -67,9 +65,6 @@ pub struct ModelParameters { /// Don't change unless you know what you're doing. #[serde(default = "default_candidate_asset_capacity")] pub candidate_asset_capacity: Capacity, - /// Defines the strategy used for calculating commodity prices - #[serde(default)] - pub pricing_strategy: PricingStrategy, /// Affects the maximum capacity that can be given to a newly created asset. /// /// It is the proportion of maximum capacity that could be required across time slices. @@ -99,18 +94,6 @@ pub struct ModelParameters { pub mothball_years: u32, } -/// The strategy used for calculating commodity prices -#[derive(DeserializeLabeledStringEnum, Debug, PartialEq, Default)] -pub enum PricingStrategy { - /// Take commodity prices directly from the shadow prices - #[default] - #[string = "shadow_prices"] - ShadowPrices, - /// Adjust shadow prices for scarcity - #[string = "scarcity_adjusted"] - ScarcityAdjusted, -} - /// Check that the `milestone_years` parameter is valid fn check_milestone_years(years: &[u32]) -> Result<()> { ensure!(!years.is_empty(), "`milestone_years` is empty"); @@ -200,21 +183,6 @@ impl ModelParameters { // milestone_years check_milestone_years(&self.milestone_years)?; - // pricing_strategy - if self.pricing_strategy == PricingStrategy::ScarcityAdjusted { - ensure!( - self.allow_broken_options, - "The pricing strategy is set to 'scarcity_adjusted', which is known to be broken. \ - If you are sure that you want to enable it anyway, you need to set the \ - {ALLOW_BROKEN_OPTION_NAME} option to true." - ); - - warn!( - "The pricing strategy is set to 'scarcity_adjusted'. Commodity prices may be \ - incorrect if assets have more than one output commodity. See: {ISSUES_URL}/677" - ); - } - // capacity_limit_factor already validated with deserialise_proportion_nonzero // candidate_asset_capacity diff --git a/src/process.rs b/src/process.rs index fa67d1645..981045059 100644 --- a/src/process.rs +++ b/src/process.rs @@ -399,17 +399,28 @@ pub struct ProcessFlow { } impl ProcessFlow { - /// Get the cost for this flow with the given parameters. + /// Get the cost per unit flow for a given region, year, and time slice. + /// + /// Includes flow costs and levies/incentives, if any. + pub fn get_total_cost_per_flow( + &self, + region_id: &RegionID, + year: u32, + time_slice: &TimeSliceID, + ) -> MoneyPerFlow { + self.cost + self.get_levy(region_id, year, time_slice) + } + + /// Get the cost for this flow per unit of activity for a given region, year, and time slice. /// /// This includes cost per unit flow and levies/incentives, if any. - pub fn get_total_cost( + pub fn get_total_cost_per_activity( &self, region_id: &RegionID, year: u32, time_slice: &TimeSliceID, ) -> MoneyPerActivity { - let cost_per_unit = self.cost + self.get_levy(region_id, year, time_slice); - + let cost_per_unit = self.get_total_cost_per_flow(region_id, year, time_slice); self.coeff.abs() * cost_per_unit } @@ -491,7 +502,7 @@ pub struct ProcessInvestmentConstraint { #[cfg(test)] mod tests { use super::*; - use crate::commodity::{CommodityLevyMap, CommodityType, DemandMap}; + use crate::commodity::{CommodityLevyMap, CommodityType, DemandMap, PricingStrategy}; use crate::fixture::{assert_error, region_id, time_slice, time_slice_info2}; use crate::time_slice::TimeSliceLevel; use crate::time_slice::TimeSliceSelection; @@ -555,6 +566,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: levies_prod, levies_cons: levies_cons, demand: DemandMap::new(), @@ -574,6 +586,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: levies, demand: DemandMap::new(), @@ -593,6 +606,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: levies, levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -614,6 +628,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: levies_prod, levies_cons: levies_cons, demand: DemandMap::new(), @@ -627,6 +642,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -641,6 +657,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -662,6 +679,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: levies, levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -683,6 +701,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: levies, levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), @@ -887,7 +906,7 @@ mod tests { time_slice: TimeSliceID, ) { assert_eq!( - flow_with_cost.get_total_cost(®ion_id, 2020, &time_slice), + flow_with_cost.get_total_cost_per_activity(®ion_id, 2020, &time_slice), MoneyPerActivity(5.0) ); } @@ -899,7 +918,7 @@ mod tests { time_slice: TimeSliceID, ) { assert_eq!( - flow_with_cost_and_levy.get_total_cost(®ion_id, 2020, &time_slice), + flow_with_cost_and_levy.get_total_cost_per_activity(®ion_id, 2020, &time_slice), MoneyPerActivity(15.0) ); } @@ -911,7 +930,7 @@ mod tests { time_slice: TimeSliceID, ) { assert_eq!( - flow_with_cost_and_incentive.get_total_cost(®ion_id, 2020, &time_slice), + flow_with_cost_and_incentive.get_total_cost_per_activity(®ion_id, 2020, &time_slice), MoneyPerActivity(2.0) ); } @@ -924,7 +943,7 @@ mod tests { ) { flow_with_cost.coeff = FlowPerActivity(-2.0); assert_eq!( - flow_with_cost.get_total_cost(®ion_id, 2020, &time_slice), + flow_with_cost.get_total_cost_per_activity(®ion_id, 2020, &time_slice), MoneyPerActivity(10.0) ); } @@ -937,7 +956,7 @@ mod tests { ) { flow_with_cost.coeff = FlowPerActivity(0.0); assert_eq!( - flow_with_cost.get_total_cost(®ion_id, 2020, &time_slice), + flow_with_cost.get_total_cost_per_activity(®ion_id, 2020, &time_slice), MoneyPerActivity(0.0) ); } @@ -949,6 +968,7 @@ mod tests { description: "Test commodity".into(), kind: CommodityType::ServiceDemand, time_slice_level: TimeSliceLevel::Annual, + pricing_strategy: PricingStrategy::Shadow, levies_prod: CommodityLevyMap::new(), levies_cons: CommodityLevyMap::new(), demand: DemandMap::new(), diff --git a/src/simulation.rs b/src/simulation.rs index cfd9fe80d..1eb45b49b 100644 --- a/src/simulation.rs +++ b/src/simulation.rs @@ -191,7 +191,8 @@ fn run_dispatch_for_year( // If there were either existing or candidate assets, we can calculate prices. // If not, return empty maps. let prices = solution_for_prices - .map(|solution| calculate_prices(model, &solution)) + .map(|solution| calculate_prices(model, &solution, year)) + .transpose()? .unwrap_or_default(); Ok((flow_map, prices)) diff --git a/src/simulation/investment/appraisal/costs.rs b/src/simulation/investment/appraisal/costs.rs index 354094090..e1761a9a2 100644 --- a/src/simulation/investment/appraisal/costs.rs +++ b/src/simulation/investment/appraisal/costs.rs @@ -1,6 +1,5 @@ //! Costs for the optimisation problem. use crate::asset::{AssetRef, AssetState}; -use crate::finance::annual_capital_cost; use crate::units::{MoneyPerCapacity, Year}; /// Calculates the annual fixed costs per unit of capacity for an asset. @@ -23,16 +22,9 @@ fn annual_fixed_cost_for_existing(asset: &AssetRef) -> MoneyPerCapacity { fixed_operating_cost * Year(1.0) } -fn annual_capital_cost_for_candidate(asset: &AssetRef) -> MoneyPerCapacity { - let capital_cost = asset.process_parameter().capital_cost; - let lifetime = asset.process_parameter().lifetime; - let discount_rate = asset.process_parameter().discount_rate; - annual_capital_cost(capital_cost, lifetime, discount_rate) -} - fn annual_fixed_cost_for_candidate(asset: &AssetRef) -> MoneyPerCapacity { let fixed_operating_cost = asset.process_parameter().fixed_operating_cost; let annual_fixed_operating_cost = fixed_operating_cost * Year(1.0); - let capital_costs = annual_capital_cost_for_candidate(asset); + let capital_costs = asset.get_annual_capital_cost_per_capacity(); annual_fixed_operating_cost + capital_costs } diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 1be56aac8..8e5d13baa 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -54,6 +54,7 @@ type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Var pub struct VariableMap { activity_vars: ActivityVariableMap, existing_asset_var_idx: Range, + candidate_asset_var_idx: Range, capacity_vars: CapacityVariableMap, capacity_var_idx: Range, unmet_demand_vars: UnmetDemandVariableMap, @@ -88,7 +89,7 @@ impl VariableMap { existing_assets, year, ); - add_activity_variables( + let candidate_asset_var_idx = add_activity_variables( problem, &mut activity_vars, &model.time_slice_info, @@ -100,6 +101,7 @@ impl VariableMap { Self { activity_vars, existing_asset_var_idx, + candidate_asset_var_idx, capacity_vars: CapacityVariableMap::new(), capacity_var_idx: Range::default(), unmet_demand_vars: UnmetDemandVariableMap::default(), @@ -190,8 +192,8 @@ pub struct Solution<'a> { impl Solution<'_> { /// Create a map of commodity flows for each asset's coeffs at every time slice. /// - /// Note that this only includes commodity flows which relate to assets, so not every commodity - /// in the simulation will necessarily be represented. + /// Note that this only includes commodity flows which relate to existing assets, so not every + /// commodity in the simulation will necessarily be represented. pub fn create_flow_map(&self) -> FlowMap { // The decision variables represent assets' activity levels, not commodity flows. We // multiply this value by the flow coeffs to get commodity flows. @@ -207,7 +209,7 @@ impl Solution<'_> { flows } - /// Activity for each existing asset + /// Activity for all assets (existing and candidate, if present) pub fn iter_activity(&self) -> impl Iterator { self.variables .activity_var_keys() @@ -216,7 +218,7 @@ impl Solution<'_> { } /// Activity for each existing asset - fn iter_activity_for_existing( + pub fn iter_activity_for_existing( &self, ) -> impl Iterator { let cols = &self.solution.columns()[self.variables.existing_asset_var_idx.clone()]; @@ -227,6 +229,18 @@ impl Solution<'_> { .map(|((asset, time_slice), &value)| (asset, time_slice, Activity(value))) } + /// Activity for each candidate asset + pub fn iter_activity_for_candidates( + &self, + ) -> impl Iterator { + let cols = &self.solution.columns()[self.variables.candidate_asset_var_idx.clone()]; + self.variables + .activity_var_keys() + .skip(self.variables.candidate_asset_var_idx.start) + .zip(cols.iter()) + .map(|((asset, time_slice), &value)| (asset, time_slice, Activity(value))) + } + /// Iterate over unmet demand pub fn iter_unmet_demand( &self, diff --git a/src/simulation/prices.rs b/src/simulation/prices.rs index ebcb47bf9..0874ed368 100644 --- a/src/simulation/prices.rs +++ b/src/simulation/prices.rs @@ -1,71 +1,106 @@ //! Code for updating the simulation state. use crate::asset::AssetRef; -use crate::commodity::CommodityID; -use crate::model::{Model, PricingStrategy}; -use crate::process::FlowDirection; +use crate::commodity::{CommodityID, PricingStrategy}; +use crate::model::Model; use crate::region::RegionID; use crate::simulation::optimisation::Solution; -use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::{Dimensionless, MoneyPerActivity, MoneyPerFlow, Year}; -use std::collections::{BTreeMap, HashMap, btree_map}; +use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; +use crate::units::{Activity, Dimensionless, MoneyPerActivity, MoneyPerFlow, Year}; +use anyhow::Result; +use itertools::iproduct; +use std::collections::{BTreeMap, HashMap, HashSet, btree_map}; + +/// Iterator item type for asset activity iterators +type Item<'a> = (&'a AssetRef, &'a TimeSliceID, Activity); /// Calculate commodity prices. /// -/// Note that the behaviour will be different depending on the [`PricingStrategy`] the user has -/// selected. +/// Prices for each commodity are calculated based on their respective pricing strategies. /// /// # Arguments /// /// * `model` - The model /// * `solution` - Solution to dispatch optimisation -pub fn calculate_prices(model: &Model, solution: &Solution) -> CommodityPrices { +/// * `year` - The year for which prices are being calculated +pub fn calculate_prices(model: &Model, solution: &Solution, year: u32) -> Result { + // Compute shadow prices for all SED/SVD commodities (needed by all strategies) let shadow_prices = CommodityPrices::from_iter(solution.iter_commodity_balance_duals()); - match model.parameters.pricing_strategy { - // Use raw shadow prices - PricingStrategy::ShadowPrices => shadow_prices, - // Adjust prices for scarcity - PricingStrategy::ScarcityAdjusted => shadow_prices - .clone() - .with_scarcity_adjustment(solution.iter_activity_duals()), + + // Partition markets by pricing strategy into a map keyed by `PricingStrategy`. + // For now, commodities use a single strategy for all regions, but this may change in the future. + let mut pricing_sets = HashMap::new(); + for ((commodity_id, commodity), region_id) in + iproduct!(&model.commodities, model.iter_regions()) + { + if commodity.pricing_strategy == PricingStrategy::Unpriced { + continue; + } + pricing_sets + .entry(&commodity.pricing_strategy) + .or_insert_with(HashSet::new) + .insert((commodity_id.clone(), region_id.clone())); } -} -/// A map relating commodity ID + region + time slice to current price (endogenous) -#[derive(Default, Clone)] -pub struct CommodityPrices(BTreeMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>); + // Set up empty prices map + let mut result = CommodityPrices::default(); -impl CommodityPrices { - /// Remove the impact of scarcity on prices. - /// - /// # Arguments - /// - /// * `activity_duals` - Value of activity duals from solution - fn with_scarcity_adjustment<'a, I>(mut self, activity_duals: I) -> Self - where - I: Iterator, - { - let highest_duals = get_highest_activity_duals(activity_duals); - - // Add the highest activity dual for each commodity/region/time slice to each commodity - // balance dual - for (key, highest) in &highest_duals { - if let Some(price) = self.0.get_mut(key) { - // highest is in units of MoneyPerActivity, but this is correct according to Adam - *price += MoneyPerFlow(highest.value()); + // Add prices for shadow-priced commodities + if let Some(shadow_set) = pricing_sets.get(&PricingStrategy::Shadow) { + for (commodity_id, region_id, time_slice) in shadow_prices.keys() { + if shadow_set.contains(&(commodity_id.clone(), region_id.clone())) { + let price = shadow_prices + .get(commodity_id, region_id, time_slice) + .unwrap(); + result.insert(commodity_id, region_id, time_slice, price); } } + } - self + // Add prices for scarcity-adjusted commodities + if let Some(scarcity_set) = pricing_sets.get(&PricingStrategy::ScarcityAdjusted) { + let scarcity_prices = calculate_scarcity_adjusted_prices( + solution.iter_activity_duals(), + &shadow_prices, + scarcity_set, + ); + result.extend(scarcity_prices); } - /// Extend the prices map, possibly overwriting values - pub fn extend(&mut self, iter: T) - where - T: IntoIterator, - { - self.0.extend(iter); + // Add prices for marginal cost commodities + if let Some(marginal_set) = pricing_sets.get(&PricingStrategy::MarginalCost) { + let marginal_cost_prices = calculate_marginal_cost_prices( + solution.iter_activity_for_existing(), + solution.iter_activity_for_candidates(), + &shadow_prices, + year, + marginal_set, + ); + result.extend(marginal_cost_prices); + } + + // Add prices for full cost commodities + if let Some(fullcost_set) = pricing_sets.get(&PricingStrategy::FullCost) { + let annual_activities = calculate_annual_activities(solution.iter_activity_for_existing()); + let full_cost_prices = calculate_full_cost_prices( + solution.iter_activity_for_existing(), + solution.iter_activity_for_candidates(), + &annual_activities, + &shadow_prices, + year, + fullcost_set, + ); + result.extend(full_cost_prices); } + // Return the completed prices map + Ok(result) +} + +/// A map relating commodity ID + region + time slice to current price (endogenous) +#[derive(Default, Clone)] +pub struct CommodityPrices(BTreeMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>); + +impl CommodityPrices { /// Insert a price for the given commodity, region and time slice pub fn insert( &mut self, @@ -78,6 +113,17 @@ impl CommodityPrices { self.0.insert(key, price); } + /// Extend the prices map, panic if any key already exists + pub fn extend(&mut self, iter: T) + where + T: IntoIterator, + { + for (key, price) in iter { + let existing = self.0.insert(key.clone(), price).is_some(); + assert!(!existing, "Key {key:?} already exists in the map"); + } + } + /// Iterate over the map. /// /// # Returns @@ -108,17 +154,6 @@ impl CommodityPrices { self.0.keys() } - /// Remove the specified entry from the map - pub fn remove( - &mut self, - commodity_id: &CommodityID, - region_id: &RegionID, - time_slice: &TimeSliceID, - ) -> Option { - self.0 - .remove(&(commodity_id.clone(), region_id.clone(), time_slice.clone())) - } - /// Calculate time slice-weighted average prices for each commodity-region pair /// /// This method aggregates prices across time slices by weighting each price @@ -209,25 +244,40 @@ impl IntoIterator for CommodityPrices { } } -fn get_highest_activity_duals<'a, I>( +/// Calculate scarcity-adjusted prices for a set of commodities. +/// +/// # Arguments +/// +/// * `activity_duals` - Iterator over activity duals from optimisation solution +/// * `shadow_prices` - Shadow prices for all commodities +/// * `markets_to_price` - Set of markets to calculate scarcity-adjusted prices for +/// +/// # Returns +/// +/// A map of scarcity-adjusted prices for the specified markets in all time slices +fn calculate_scarcity_adjusted_prices<'a, I>( activity_duals: I, -) -> HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerActivity> + shadow_prices: &CommodityPrices, + markets_to_price: &HashSet<(CommodityID, RegionID)>, +) -> HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow> where I: Iterator, { // Calculate highest activity dual for each commodity/region/time slice let mut highest_duals = HashMap::new(); for (asset, time_slice, dual) in activity_duals { - // Iterate over all output flows - for flow in asset - .iter_flows() - .filter(|flow| flow.direction() == FlowDirection::Output) - { + let region_id = asset.region_id(); + + // Iterate over the output flows of this asset + // Only consider flows for commodities we are pricing + for flow in asset.iter_output_flows().filter(|flow| { + markets_to_price.contains(&(flow.commodity.id.clone(), region_id.clone())) + }) { // Update the highest dual for this commodity/time slice highest_duals .entry(( flow.commodity.id.clone(), - asset.region_id().clone(), + region_id.clone(), time_slice.clone(), )) .and_modify(|current_dual| { @@ -239,17 +289,434 @@ where } } - highest_duals + // Add this to the shadow price for each commodity/region/time slice + let mut scarcity_prices = HashMap::new(); + for ((commodity, region, time_slice), highest_dual) in &highest_duals { + // There should always be a shadow price for commodities we are considering here, so it + // should be safe to unwrap + let shadow_price = shadow_prices.get(commodity, region, time_slice).unwrap(); + // highest_dual is in units of MoneyPerActivity, and shadow_price is in MoneyPerFlow, but + // this is correct according to Adam + let scarcity_price = shadow_price + MoneyPerFlow(highest_dual.value()); + scarcity_prices.insert( + (commodity.clone(), region.clone(), time_slice.clone()), + scarcity_price, + ); + } + + scarcity_prices +} + +/// Calculate marginal cost prices for a set of commodities. +/// +/// This pricing strategy aims to incorporate the marginal cost of commodity production into the price. +/// +/// For a given asset, a marginal cost can be calculated for each SED/SVD output, which is the sum of: +/// - Generic activity costs: Activity-related costs not tied to a specific SED/SVD output +/// (variable operating costs, cost of purchasing inputs, plus all levies and flow costs not +/// associated with specific SED/SVD outputs). These are shared over all SED/SVD +/// outputs according to their flow coefficients. +/// - Commodity-specific activity costs: flow costs/levies for the specific SED/SVD output. +/// +/// --- +/// +/// For example, consider an asset A(SED) -> B(SED) + 2C(SED) + D(OTH), with the following costs: +/// - Variable operating cost: 5 per unit activity +/// - Production levy on C: 3 per unit flow +/// - Production levy on D: 4 per unit flow +/// - Shadow price of A: 1 per unit flow +/// +/// Then: +/// - Generic activity cost per activity = (1 + 5 + 4) = 10 +/// - Generic activity cost per SED/SVD output = 10 / (1 + 2) = 3.333 +/// - Marginal cost of B = 3.333 +/// - Marginal cost of C = 3.333 + 3 = 6.333 +/// +/// --- +/// +/// If any existing assets produce a given commodity in a particular region and time slice, the +/// price is taken from the asset with the highest marginal cost among those existing assets. If _no_ +/// existing assets produce the commodity in that region and time slice (in particular, this will +/// occur when there's no demand for the commodity), then candidate assets are considered: we +/// take the price from the candidate asset with the _lowest_ marginal cost, assuming full utilisation +/// (i.e. the single candidate asset that would be most competitive if a small amount of demand was +/// added). +/// +/// Note: this should be similar to the "shadow price" strategy, which is also based on marginal +/// costs of the most expensive producer, but may be more successful in cases where there are +// multiple SED/SVD outputs per asset. +/// +/// # Arguments +/// +/// * `activity_for_existing` - Iterator over activity from optimisation solution for existing +/// assets +/// * `activity_for_candidates` - Iterator over activity from optimisation solution for candidate +/// assets. Note: we only need the keys, since we assume full utilisation for candidates. +/// * `annual_activities` - Map of annual activities for each asset computed by +/// `calculate_annual_activities`. This only needs to include existing assets. +/// * `shadow_prices` - Shadow prices for all commodities +/// * `year` - The year for which prices are being calculated +/// * `markets_to_price` - Set of markets to calculate full cost prices for +/// +/// # Returns +/// +/// A map of marginal cost prices for the specified markets in all time slices +fn calculate_marginal_cost_prices<'a, I, J>( + activity_for_existing: I, + activity_for_candidates: J, + shadow_prices: &CommodityPrices, + year: u32, + markets_to_price: &HashSet<(CommodityID, RegionID)>, +) -> HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow> +where + I: Iterator>, + J: Iterator>, +{ + let mut prices: HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow> = HashMap::new(); + + // Start by looking at existing assets + // Calculate highest marginal cost for each commodity/region/time slice + // Keep track of keys with prices - missing keys will be handled by candidates later + let mut priced_by_existing = HashSet::new(); + for (asset, time_slice, activity) in activity_for_existing { + let region_id = asset.region_id(); + + // Only proceed if the asset has non-zero activity in this time slice + if activity < Activity::EPSILON { + continue; + } + + // Iterate over all the SED/SVD marginal costs for commodities we need prices for + for (commodity_id, marginal_cost) in asset.iter_marginal_costs_with_filter( + shadow_prices, + year, + time_slice, + |commodity_id: &CommodityID| { + markets_to_price.contains(&(commodity_id.clone(), region_id.clone())) + }, + ) { + // Update the highest cost for this commodity/region/time slice + let key = (commodity_id.clone(), region_id.clone(), time_slice.clone()); + prices + .entry(key.clone()) + .and_modify(|c| *c = c.max(marginal_cost)) + .or_insert(marginal_cost); + priced_by_existing.insert(key); + } + } + + // Next, look at candidate assets for any markets not covered by existing assets + // For these, we take the _lowest_ marginal cost + for (asset, time_slice, _activity) in activity_for_candidates { + let region_id = asset.region_id(); + + // Only consider markets not already priced by existing assets + let should_process = |cid: &CommodityID| { + markets_to_price.contains(&(cid.clone(), region_id.clone())) + && !priced_by_existing.contains(&( + cid.clone(), + region_id.clone(), + time_slice.clone(), + )) + }; + + // Iterate over all the SED/SVD marginal costs for markets we need prices for + for (commodity_id, marginal_cost) in asset.iter_marginal_costs_with_filter( + shadow_prices, + year, + time_slice, + |cid: &CommodityID| should_process(cid), + ) { + // Update the _lowest_ cost for this commodity/region/time slice + let key = (commodity_id.clone(), region_id.clone(), time_slice.clone()); + prices + .entry(key.clone()) + .and_modify(|c| *c = c.min(marginal_cost)) + .or_insert(marginal_cost); + } + } + + // Return the calculated marginal prices + prices +} + +/// Calculated annual activities for each asset by summing across all time slices +fn calculate_annual_activities<'a, I>(activities: I) -> HashMap +where + I: IntoIterator>, +{ + activities + .into_iter() + .map(|(asset, _ts, activity)| (asset.clone(), activity)) + .fold(HashMap::new(), |mut acc, (asset, activity)| { + acc.entry(asset) + .and_modify(|e| *e += activity) + .or_insert(activity); + acc + }) +} + +/// Calculate full cost prices for a set of commodities. +/// +/// This pricing strategy aims to incorporate the full cost of commodity production into the price. +/// +/// For a given asset, a full cost can be calculated for each SED/SVD output, which is the sum of: +/// - Annual capital costs/fixed operating costs: Calculated based on the capacity of the asset +/// and the total annual output. If an asset has multiple SED/SVD outputs, then these costs are +/// shared equally over all outputs according to their flow coefficients. +/// - Generic activity costs: Activity-related costs not tied to a specific SED/SVD output +/// (variable operating costs, cost of purchasing inputs, plus all levies and flow costs not +/// associated with specific SED/SVD outputs). As above, these are shared over all SED/SVD +/// outputs according to their flow coefficients. +/// - Commodity-specific activity costs: flow costs/levies for the specific SED/SVD output. +/// +/// --- +/// +/// For example, consider an asset A(SED) -> B(SED) + 2C(SED) + D(OTH), with the following costs: +/// - Annual capital cost + fixed operating cost: 2.5 per unit capacity +/// - Variable operating cost: 5 per unit activity +/// - Production levy on C: 3 per unit flow +/// - Production levy on D: 4 per unit flow +/// - Shadow price of A: 1 per unit flow +/// +/// If capacity is 4 and annual activity is 2: +/// - Annual capital + fixed operating cost per activity = (2.5 * 4) / 2 = 5 +/// - Annual capital + fixed operating cost per SED/SVD output = 5 / (1 + 2) = 1.666 +/// - Generic activity cost per activity = (1 + 5 + 4) = 10 +/// - Generic activity cost per SED/SVD output = 10 / (1 + 2) = 3.333 +/// - Full cost of B = 1.666 + 3.333 = 5.0 +/// - Full cost of C = 1.666 + 3.333 + 3 = 8.0 +/// +/// --- +/// +/// If any existing assets produce a given commodity in a particular region and time slice, the +/// price is taken from the asset with the highest full cost among those existing assets. If _no_ +/// existing assets produce the commodity in that region and time slice (in particular, this will +/// occur when there's no demand for the commodity), then candidate assets are considered: we +/// take the price from the candidate asset with the _lowest_ full cost, assuming maximum +/// possible dispatch (i.e. the single candidate asset that would be most competitive if a small +/// amount of demand was added). +/// +/// # Arguments +/// +/// * `activity_for_existing` - Iterator over activity from optimisation solution for existing +/// assets +/// * `activity_for_candidates` - Iterator over activity from optimisation solution for candidate +/// assets. Note: we only need the keys, since we assume full dispatch for candidates. +/// * `annual_activities` - Map of annual activities for each asset computed by +/// `calculate_annual_activities`. This only needs to include existing assets. +/// * `shadow_prices` - Shadow prices for all commodities +/// * `year` - The year for which prices are being calculated +/// * `markets_to_price` - Set of markets to calculate full cost prices for +/// +/// # Returns +/// +/// A map of full cost prices for the specified markets in all time slices +fn calculate_full_cost_prices<'a, I, J>( + activity_for_existing: I, + activity_for_candidates: J, + annual_activities: &HashMap, + shadow_prices: &CommodityPrices, + year: u32, + markets_to_price: &HashSet<(CommodityID, RegionID)>, +) -> HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow> +where + I: Iterator>, + J: Iterator>, +{ + let mut prices: HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow> = HashMap::new(); + + // Start by looking at existing assets + // Calculate highest full cost for each commodity/region/time slice + // Keep track of keys with prices - missing keys will be handled by candidates later + let mut annual_capital_costs_cache = HashMap::new(); + let mut priced_by_existing = HashSet::new(); + for (asset, time_slice, activity) in activity_for_existing { + let annual_activity = annual_activities[asset]; + let region_id = asset.region_id(); + + // Only proceed if the asset has non-zero activity in this time slice + if activity < Activity::EPSILON { + continue; + } + + // Only proceed if the asset produces at least one commodity we need prices for + if !asset + .iter_output_flows() + .any(|flow| markets_to_price.contains(&(flow.commodity.id.clone(), region_id.clone()))) + { + continue; + } + + // Calculate/cache annual capital cost for this asset + let annual_capital_cost_per_flow = *annual_capital_costs_cache + .entry(asset.clone()) + .or_insert_with(|| asset.get_annual_capital_cost_per_flow(annual_activity)); + + // Iterate over all the SED/SVD marginal costs for commodities we need prices for + for (commodity_id, marginal_cost) in asset.iter_marginal_costs_with_filter( + shadow_prices, + year, + time_slice, + |cid: &CommodityID| markets_to_price.contains(&(cid.clone(), region_id.clone())), + ) { + // Add capital cost per flow to marginal cost to get full cost + let marginal_cost = marginal_cost + annual_capital_cost_per_flow; + + // Update the highest cost for this commodity/region/time slice + let key = (commodity_id.clone(), region_id.clone(), time_slice.clone()); + prices + .entry(key.clone()) + .and_modify(|c| *c = c.max(marginal_cost)) + .or_insert(marginal_cost); + priced_by_existing.insert(key); + } + } + + // Next, look at candidate assets for any markets not covered by existing assets + // For these we assume full utilisation, and take the _lowest_ full cost + for (asset, time_slice, _activity) in activity_for_candidates { + let region_id = asset.region_id(); + + // Only consider markets not already priced by existing assets + let should_process = |cid: &CommodityID| { + markets_to_price.contains(&(cid.clone(), region_id.clone())) + && !priced_by_existing.contains(&( + cid.clone(), + region_id.clone(), + time_slice.clone(), + )) + }; + + // Only proceed if the asset produces at least one commodity we need prices for + if !asset + .iter_output_flows() + .any(|flow| should_process(&flow.commodity.id)) + { + continue; + } + + // Calculate/cache annual capital cost per flow for this asset assuming full dispatch + // (bound by the activity limits of the asset) + let annual_capital_cost_per_flow = *annual_capital_costs_cache + .entry(asset.clone()) + .or_insert_with(|| { + asset.get_annual_capital_cost_per_flow( + *asset + .get_activity_limits_for_selection(&TimeSliceSelection::Annual) + .end(), + ) + }); + + // Iterate over all the SED/SVD marginal costs for markets we need prices for + for (commodity_id, marginal_cost) in asset.iter_marginal_costs_with_filter( + shadow_prices, + year, + time_slice, + |cid: &CommodityID| should_process(cid), + ) { + // Add capital cost per flow to marginal cost to get full cost + let full_cost = marginal_cost + annual_capital_cost_per_flow; + + // Update the _lowest_ cost for this commodity/region/time slice + let key = (commodity_id.clone(), region_id.clone(), time_slice.clone()); + prices + .entry(key) + .and_modify(|c| *c = c.min(full_cost)) + .or_insert(full_cost); + } + } + + // Return the calculated full cost prices + prices } #[cfg(test)] mod tests { use super::*; - use crate::commodity::CommodityID; - use crate::fixture::{commodity_id, region_id, time_slice, time_slice_info}; + use crate::asset::Asset; + use crate::asset::AssetRef; + use crate::commodity::{Commodity, CommodityID}; + use crate::fixture::{ + commodity_id, other_commodity, region_id, sed_commodity, time_slice, time_slice_info, + }; + use crate::process::{ActivityLimits, FlowType, Process, ProcessFlow, ProcessParameter}; use crate::region::RegionID; use crate::time_slice::TimeSliceID; + use crate::units::ActivityPerCapacity; + use crate::units::{ + Activity, Capacity, Dimensionless, FlowPerActivity, MoneyPerActivity, MoneyPerCapacity, + MoneyPerCapacityPerYear, MoneyPerFlow, + }; + use indexmap::{IndexMap, IndexSet}; use rstest::rstest; + use std::collections::{HashMap, HashSet}; + use std::rc::Rc; + + fn build_process_flow(commodity: &Commodity, coeff: f64, cost: MoneyPerFlow) -> ProcessFlow { + ProcessFlow { + commodity: Rc::new(commodity.clone()), + coeff: FlowPerActivity(coeff), + kind: FlowType::Fixed, + cost, + } + } + + fn build_process( + flows: IndexMap, + region_id: &RegionID, + year: u32, + time_slice_info: &TimeSliceInfo, + variable_operating_cost: MoneyPerActivity, + capital_cost: MoneyPerCapacity, + lifetime: u32, + discount_rate: Dimensionless, + ) -> Process { + let mut process_flows_map = HashMap::new(); + process_flows_map.insert((region_id.clone(), year), Rc::new(flows)); + + let mut process_parameter_map = HashMap::new(); + let proc_param = ProcessParameter { + capital_cost, + fixed_operating_cost: MoneyPerCapacityPerYear(0.0), + variable_operating_cost, + lifetime, + discount_rate, + }; + process_parameter_map.insert((region_id.clone(), year), Rc::new(proc_param)); + + let mut activity_limits_map = HashMap::new(); + activity_limits_map.insert( + (region_id.clone(), year), + Rc::new(ActivityLimits::new_with_full_availability(time_slice_info)), + ); + + let regions: IndexSet = IndexSet::from([region_id.clone()]); + + Process { + id: "p1".into(), + description: "test process".into(), + years: 2010..=2020, + activity_limits: activity_limits_map, + flows: process_flows_map, + parameters: process_parameter_map, + regions, + primary_output: None, + capacity_to_activity: ActivityPerCapacity(1.0), + investment_constraints: HashMap::new(), + } + } + + fn assert_price_approx( + prices: &HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow>, + commodity: &CommodityID, + region: &RegionID, + time_slice: &TimeSliceID, + expected: MoneyPerFlow, + ) { + let p = prices[&(commodity.clone(), region.clone(), time_slice.clone())]; + assert!((p - expected).abs() < MoneyPerFlow::EPSILON); + } #[rstest] #[case(MoneyPerFlow(100.0), MoneyPerFlow(100.0), Dimensionless(0.0), true)] // exactly equal @@ -302,4 +769,146 @@ mod tests { // With single time slice (duration=1.0), average should equal the price assert_eq!(averages[&(commodity_id, region_id)], MoneyPerFlow(100.0)); } + + #[rstest] + fn test_marginal_cost_example( + sed_commodity: Commodity, + other_commodity: Commodity, + region_id: RegionID, + time_slice_info: TimeSliceInfo, + time_slice: TimeSliceID, + ) { + // Use the same setup as in the docstring example for marginal cost pricing + let mut a = sed_commodity.clone(); + a.id = "A".into(); + let mut b = sed_commodity.clone(); + b.id = "B".into(); + let mut c = sed_commodity.clone(); + c.id = "C".into(); + let mut d = other_commodity.clone(); + d.id = "D".into(); + + let mut flows = IndexMap::new(); + flows.insert( + a.id.clone(), + build_process_flow(&a, -1.0, MoneyPerFlow(0.0)), + ); + flows.insert(b.id.clone(), build_process_flow(&b, 1.0, MoneyPerFlow(0.0))); + flows.insert(c.id.clone(), build_process_flow(&c, 2.0, MoneyPerFlow(3.0))); + flows.insert(d.id.clone(), build_process_flow(&d, 1.0, MoneyPerFlow(4.0))); + + let process = build_process( + flows, + ®ion_id, + 2015u32, + &time_slice_info, + MoneyPerActivity(5.0), // variable operating cost + MoneyPerCapacity(0.0), // capital cost + 5, // lifetime + Dimensionless(1.0), // discount rate + ); + + let asset = + Asset::new_candidate(Rc::new(process), region_id.clone(), Capacity(1.0), 2015u32) + .unwrap(); + let asset_ref = AssetRef::from(asset); + let shadow_prices = + CommodityPrices::from_iter(vec![(&a.id, ®ion_id, &time_slice, MoneyPerFlow(1.0))]); + let mut markets = HashSet::new(); + markets.insert((b.id.clone(), region_id.clone())); + markets.insert((c.id.clone(), region_id.clone())); + + let existing = vec![(&asset_ref, &time_slice, Activity(1.0))]; + let candidates = Vec::new(); + + let prices = calculate_marginal_cost_prices( + existing.into_iter(), + candidates.into_iter(), + &shadow_prices, + 2015u32, + &markets, + ); + + assert_price_approx( + &prices, + &b.id, + ®ion_id, + &time_slice, + MoneyPerFlow(10.0 / 3.0), + ); + assert_price_approx( + &prices, + &c.id, + ®ion_id, + &time_slice, + MoneyPerFlow(10.0 / 3.0 + 3.0), + ); + } + + #[rstest] + fn test_full_cost_example( + sed_commodity: Commodity, + other_commodity: Commodity, + region_id: RegionID, + time_slice_info: TimeSliceInfo, + time_slice: TimeSliceID, + ) { + // Use the same setup as in the docstring example for full cost pricing + let mut a = sed_commodity.clone(); + a.id = "A".into(); + let mut b = sed_commodity.clone(); + b.id = "B".into(); + let mut c = sed_commodity.clone(); + c.id = "C".into(); + let mut d = other_commodity.clone(); + d.id = "D".into(); + + let mut flows = IndexMap::new(); + flows.insert( + a.id.clone(), + build_process_flow(&a, -1.0, MoneyPerFlow(0.0)), + ); + flows.insert(b.id.clone(), build_process_flow(&b, 1.0, MoneyPerFlow(0.0))); + flows.insert(c.id.clone(), build_process_flow(&c, 2.0, MoneyPerFlow(3.0))); + flows.insert(d.id.clone(), build_process_flow(&d, 1.0, MoneyPerFlow(4.0))); + + let process = build_process( + flows, + ®ion_id, + 2015u32, + &time_slice_info, + MoneyPerActivity(5.0), // variable operating cost + MoneyPerCapacity(2.5), // capital cost per capacity so annualised=2.5 + 1, // lifetime so annualised = capital_cost + Dimensionless(0.0), // discount rate + ); + + let asset = + Asset::new_candidate(Rc::new(process), region_id.clone(), Capacity(4.0), 2015u32) + .unwrap(); + let asset_ref = AssetRef::from(asset); + let shadow_prices = + CommodityPrices::from_iter(vec![(&a.id, ®ion_id, &time_slice, MoneyPerFlow(1.0))]); + let mut markets = HashSet::new(); + markets.insert((b.id.clone(), region_id.clone())); + markets.insert((c.id.clone(), region_id.clone())); + + let existing = vec![(&asset_ref, &time_slice, Activity(2.0))]; + let candidates = Vec::new(); + + let mut annual_activities = HashMap::new(); + annual_activities.insert(asset_ref.clone(), Activity(2.0)); + + let prices = calculate_full_cost_prices( + existing.into_iter(), + candidates.into_iter(), + &annual_activities, + &shadow_prices, + 2015u32, + &markets, + ); + + assert_price_approx(&prices, &b.id, ®ion_id, &time_slice, MoneyPerFlow(5.0)); + assert_price_approx(&prices, &c.id, ®ion_id, &time_slice, MoneyPerFlow(8.0)); + } } diff --git a/src/units.rs b/src/units.rs index c37a7a68b..fe71d59c2 100644 --- a/src/units.rs +++ b/src/units.rs @@ -97,6 +97,9 @@ macro_rules! base_unit_struct { } } impl $name { + /// Small epsilon constant for this unit type. + pub const EPSILON: $name = $name(f64::EPSILON); + /// Create from an f64 value pub fn new(value: f64) -> Self { $name(value)