diff --git a/src/fixture.rs b/src/fixture.rs index 0cefa572b..aacce2140 100644 --- a/src/fixture.rs +++ b/src/fixture.rs @@ -21,7 +21,7 @@ use crate::simulation::investment::appraisal::{ use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceLevel}; use crate::units::{ Activity, ActivityPerCapacity, Capacity, Dimensionless, Flow, MoneyPerActivity, - MoneyPerCapacity, MoneyPerCapacityPerYear, MoneyPerFlow, Year, + MoneyPerCapacity, MoneyPerCapacityPerYear, Year, }; use anyhow::Result; use indexmap::indexmap; @@ -407,9 +407,8 @@ pub fn appraisal_output(asset: Asset, time_slice: TimeSliceID) -> AppraisalOutpu asset: AssetRef::from(asset), capacity: AssetCapacity::Continuous(Capacity(42.0)), coefficients: Rc::new(ObjectiveCoefficients { - capacity_coefficient: MoneyPerCapacity(2.14), activity_coefficients, - unmet_demand_coefficient: MoneyPerFlow(10000.0), + lcox_costs: IndexMap::new(), }), activity, unmet_demand, diff --git a/src/output.rs b/src/output.rs index 8df0843db..ce9d94fe7 100644 --- a/src/output.rs +++ b/src/output.rs @@ -8,9 +8,7 @@ use crate::simulation::CommodityPrices; use crate::simulation::investment::appraisal::AppraisalOutput; use crate::simulation::optimisation::{FlowMap, Solution}; use crate::time_slice::TimeSliceID; -use crate::units::{ - Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow, -}; +use crate::units::{Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerFlow}; use anyhow::{Context, Result, ensure}; use csv; use indexmap::IndexMap; @@ -261,7 +259,6 @@ struct AppraisalResultsRow { process_id: ProcessID, region_id: RegionID, capacity: Capacity, - capacity_coefficient: MoneyPerCapacity, metric: Option, } @@ -495,7 +492,6 @@ impl DebugDataWriter { process_id: result.asset.process_id().clone(), region_id: result.asset.region_id().clone(), capacity: result.capacity.total_capacity(), - capacity_coefficient: result.coefficients.capacity_coefficient, metric: result.metric.as_ref().map(|m| m.value()), }; self.appraisal_results_writer.serialize(row)?; @@ -1195,7 +1191,6 @@ mod tests { process_id: asset.process_id().clone(), region_id: asset.region_id().clone(), capacity: Capacity(42.0), - capacity_coefficient: MoneyPerCapacity(2.14), metric: Some(4.14), }; let records: Vec = diff --git a/src/simulation/investment/appraisal.rs b/src/simulation/investment/appraisal.rs index 47662a60c..3462da246 100644 --- a/src/simulation/investment/appraisal.rs +++ b/src/simulation/investment/appraisal.rs @@ -258,26 +258,19 @@ fn calculate_lcox( coefficients: &Rc, demand: &DemandMap, ) -> Result { - let results = perform_optimisation( - model, - asset, - max_capacity, - commodity, - coefficients, - demand, - highs::Sense::Minimise, - )?; + let results = + perform_optimisation(model, asset, max_capacity, commodity, coefficients, demand)?; let cost_index = lcox( - results.capacity.total_capacity(), - coefficients.capacity_coefficient, + max_capacity.total_capacity(), + annual_fixed_cost(asset), &results.activity, - &coefficients.activity_coefficients, + &coefficients.lcox_costs, ); Ok(AppraisalOutput::new( asset.clone(), - results.capacity, + max_capacity, results, cost_index.map(LCOXMetric::new), coefficients.clone(), @@ -297,15 +290,8 @@ fn calculate_npv( coefficients: &Rc, demand: &DemandMap, ) -> Result { - let results = perform_optimisation( - model, - asset, - max_capacity, - commodity, - coefficients, - demand, - highs::Sense::Maximise, - )?; + let results = + perform_optimisation(model, asset, max_capacity, commodity, coefficients, demand)?; let annual_fixed_cost = annual_fixed_cost(asset); assert!( @@ -314,7 +300,7 @@ fn calculate_npv( ); let profitability_index = profitability_index( - results.capacity.total_capacity(), + max_capacity.total_capacity(), annual_fixed_cost, &results.activity, &coefficients.activity_coefficients, @@ -405,7 +391,7 @@ mod tests { use crate::fixture::{agent_id, asset, process, region_id}; use crate::process::Process; use crate::region::RegionID; - use crate::units::{Money, MoneyPerActivity, MoneyPerFlow}; + use crate::units::{Money, MoneyPerActivity}; use float_cmp::assert_approx_eq; use rstest::rstest; use std::rc::Rc; @@ -574,9 +560,8 @@ mod tests { fn objective_coeffs() -> Rc { Rc::new(ObjectiveCoefficients { - capacity_coefficient: MoneyPerCapacity(0.0), activity_coefficients: IndexMap::new(), - unmet_demand_coefficient: MoneyPerFlow(0.0), + lcox_costs: IndexMap::new(), }) } diff --git a/src/simulation/investment/appraisal/coefficients.rs b/src/simulation/investment/appraisal/coefficients.rs index 286b93535..90bebc03f 100644 --- a/src/simulation/investment/appraisal/coefficients.rs +++ b/src/simulation/investment/appraisal/coefficients.rs @@ -1,11 +1,11 @@ //! Calculation of cost coefficients for investment tools. -use super::costs::annual_fixed_cost; use crate::agent::ObjectiveType; use crate::asset::AssetRef; use crate::model::Model; + use crate::simulation::CommodityPrices; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::{MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow}; +use crate::units::MoneyPerActivity; use indexmap::IndexMap; use std::collections::HashMap; use std::rc::Rc; @@ -17,12 +17,10 @@ use std::rc::Rc; /// coefficients used in the appraisal optimisation, together with the unmet-demand penalty. #[derive(Clone)] pub struct ObjectiveCoefficients { - /// Cost per unit of capacity - pub capacity_coefficient: MoneyPerCapacity, /// Cost per unit of activity in each time slice pub activity_coefficients: IndexMap, - /// Unmet demand coefficient - pub unmet_demand_coefficient: MoneyPerFlow, + /// Costs for LCOX + pub lcox_costs: IndexMap, } /// Calculates cost coefficients for a set of assets for a given objective type. @@ -36,93 +34,64 @@ pub fn calculate_coefficients_for_assets( assets .iter() .map(|asset| { - let coefficient = match objective_type { - ObjectiveType::LevelisedCostOfX => calculate_coefficients_for_lcox( - asset, - &model.time_slice_info, - prices, - model.parameters.value_of_lost_load, - year, - ), - ObjectiveType::NetPresentValue => { - calculate_coefficients_for_npv(asset, &model.time_slice_info, prices, year) - } - }; - (asset.clone(), Rc::new(coefficient)) + let mut coeff = + calculate_coefficients_for_asset(asset, &model.time_slice_info, prices, year); + + // For LCOX, we also store per-time slice costs, which are used to calculate the cost + // index + if objective_type == &ObjectiveType::LevelisedCostOfX { + coeff.lcox_costs = + calculate_lcox_costs(asset, &model.time_slice_info, prices, year); + } + + (asset.clone(), Rc::new(coeff)) }) .collect() } -/// Calculates the cost coefficients for LCOX. +/// Calculates the cost coefficients for a single asset. /// -/// For LCOX the activity coefficient is calculated as operating cost minus revenue from -/// non-primary flows. The unmet demand coefficient is set from the model parameter -/// `value_of_lost_load`. -pub fn calculate_coefficients_for_lcox( +/// The activity coefficient is revenue minus operating cost; a small positive epsilon is added to +/// activity coefficients so that assets with near-zero net value still appear in dispatch. Capacity +/// costs and unmet-demand penalties are set to zero. +pub fn calculate_coefficients_for_asset( asset: &AssetRef, time_slice_info: &TimeSliceInfo, prices: &CommodityPrices, - value_of_lost_load: MoneyPerFlow, year: u32, ) -> ObjectiveCoefficients { - // Capacity coefficient - let capacity_coefficient = annual_fixed_cost(asset); - // Activity coefficients let mut activity_coefficients = IndexMap::new(); for time_slice in time_slice_info.iter_ids() { - let coefficient = calculate_activity_coefficient_for_lcox(asset, time_slice, prices, year); + let coefficient = calculate_activity_coefficient(asset, time_slice, prices, year); activity_coefficients.insert(time_slice.clone(), coefficient); } - // Unmet demand coefficient - let unmet_demand_coefficient = value_of_lost_load; - ObjectiveCoefficients { - capacity_coefficient, activity_coefficients, - unmet_demand_coefficient, + lcox_costs: IndexMap::new(), } } -/// Calculates the cost coefficients for NPV. -/// -/// For NPV the activity coefficient is revenue (including primary output) minus operating -/// cost; a small positive epsilon is added to activity coefficients so that assets with -/// near-zero net value still appear in dispatch. Capacity costs and unmet-demand penalties -/// are set to zero for the NPV objective. -pub fn calculate_coefficients_for_npv( +/// Calculate costs for LCOX for all time slices +fn calculate_lcox_costs( asset: &AssetRef, time_slice_info: &TimeSliceInfo, prices: &CommodityPrices, year: u32, -) -> ObjectiveCoefficients { - // Small constant added to each activity coefficient to ensure break-even/slightly negative - // assets are still dispatched - const EPSILON_ACTIVITY_COEFFICIENT: MoneyPerActivity = MoneyPerActivity(f64::EPSILON * 100.0); - - // Activity coefficients - let mut activity_coefficients = IndexMap::new(); - for time_slice in time_slice_info.iter_ids() { - let coefficient = calculate_activity_coefficient_for_npv(asset, time_slice, prices, year); - activity_coefficients.insert( - time_slice.clone(), - coefficient + EPSILON_ACTIVITY_COEFFICIENT, - ); - } - - // Unmet demand coefficient (we don't apply a cost to unmet demand, so we set this to zero) - let unmet_demand_coefficient = MoneyPerFlow(0.0); - - ObjectiveCoefficients { - capacity_coefficient: MoneyPerCapacity(0.0), - activity_coefficients, - unmet_demand_coefficient, - } +) -> IndexMap { + time_slice_info + .iter_ids() + .cloned() + .map(|time_slice| { + let coefficient = calculate_lcox_cost_for_time_slice(asset, &time_slice, prices, year); + (time_slice, coefficient) + }) + .collect() } -/// Calculate a single activity coefficient for the LCOX objective for a given time slice. -fn calculate_activity_coefficient_for_lcox( +/// Calculate cost per unit activity excluding the primary commodity for a given time slice. +fn calculate_lcox_cost_for_time_slice( asset: &AssetRef, time_slice: &TimeSliceID, prices: &CommodityPrices, @@ -139,13 +108,17 @@ fn calculate_activity_coefficient_for_lcox( operating_cost - revenue_from_flows } -/// Calculate a single activity coefficient for the NPV objective for a given time slice. -fn calculate_activity_coefficient_for_npv( +/// Calculate a single activity coefficient for a given time slice +fn calculate_activity_coefficient( asset: &AssetRef, time_slice: &TimeSliceID, prices: &CommodityPrices, year: u32, ) -> MoneyPerActivity { + // Small constant added to each activity coefficient to ensure break-even/slightly negative + // assets are still dispatched + const EPSILON_ACTIVITY_COEFFICIENT: MoneyPerActivity = MoneyPerActivity(f64::EPSILON * 100.0); + // Get the operating cost of the asset. This includes the variable operating cost, levies and // flow costs, but excludes costs/revenues from commodity consumption/production. let operating_cost = asset.get_operating_cost(year, time_slice); @@ -154,5 +127,5 @@ fn calculate_activity_coefficient_for_npv( let revenue_from_flows = asset.get_revenue_from_flows(prices, time_slice); // The activity coefficient is the revenue from flows minus the operating cost (net revenue) - revenue_from_flows - operating_cost + revenue_from_flows - operating_cost + EPSILON_ACTIVITY_COEFFICIENT } diff --git a/src/simulation/investment/appraisal/constraints.rs b/src/simulation/investment/appraisal/constraints.rs index 750cb85a3..e61089239 100644 --- a/src/simulation/investment/appraisal/constraints.rs +++ b/src/simulation/investment/appraisal/constraints.rs @@ -1,45 +1,13 @@ //! Constraints for the optimisation problem. use super::DemandMap; use super::optimisation::Variable; -use crate::asset::{AssetCapacity, AssetRef, AssetState}; +use crate::asset::{AssetCapacity, AssetRef}; use crate::commodity::Commodity; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; use crate::units::Flow; use highs::RowProblem as Problem; use indexmap::IndexMap; -/// Adds a capacity constraint to the problem. -/// -/// The behaviour depends on whether the asset is commissioned or a candidate: -/// - For a commissioned asset, the capacity is fixed. -/// - For a candidate asset, the capacity is variable between zero and an upper bound. -pub fn add_capacity_constraint( - problem: &mut Problem, - asset: &AssetRef, - max_capacity: AssetCapacity, - capacity_var: Variable, -) { - let capacity_limit = match max_capacity { - AssetCapacity::Continuous(cap) => cap.value(), - AssetCapacity::Discrete(units, _) => units as f64, - }; - - let bounds = match asset.state() { - AssetState::Commissioned { .. } => { - // Fixed capacity for commissioned assets - capacity_limit..=capacity_limit - } - AssetState::Candidate => { - // Variable capacity between 0 and max for candidate assets - 0.0..=capacity_limit - } - _ => panic!( - "add_capacity_constraint should only be called with Commissioned or Candidate assets" - ), - }; - problem.add_row(bounds, [(capacity_var, 1.0)]); -} - /// Adds activity constraints to the problem. /// /// Constrains the activity variables to be within the asset's activity limits. @@ -53,37 +21,13 @@ pub fn add_capacity_constraint( pub fn add_activity_constraints( problem: &mut Problem, asset: &AssetRef, - capacity_var: Variable, - activity_vars: &IndexMap, - time_slice_info: &TimeSliceInfo, -) { - match asset.state() { - AssetState::Commissioned { .. } => { - add_activity_constraints_for_existing(problem, asset, activity_vars, time_slice_info); - } - AssetState::Candidate => { - add_activity_constraints_for_candidate( - problem, - asset, - capacity_var, - activity_vars, - time_slice_info, - ); - } - _ => panic!( - "add_activity_constraints should only be called with Commissioned or Candidate assets" - ), - } -} - -fn add_activity_constraints_for_existing( - problem: &mut Problem, - asset: &AssetRef, + max_capacity: AssetCapacity, activity_vars: &IndexMap, time_slice_info: &TimeSliceInfo, ) { - for (ts_selection, limits) in asset.iter_activity_limits() { - let limits = limits.start().value()..=limits.end().value(); + let capacity = max_capacity.total_capacity(); + for (ts_selection, limits) in asset.iter_activity_per_capacity_limits() { + let limits = (capacity * *limits.start()).value()..=(capacity * *limits.end()).value(); // Collect activity terms for the time slices in this selection let terms = ts_selection @@ -96,42 +40,6 @@ fn add_activity_constraints_for_existing( } } -fn add_activity_constraints_for_candidate( - problem: &mut Problem, - asset: &AssetRef, - capacity_var: Variable, - activity_vars: &IndexMap, - time_slice_info: &TimeSliceInfo, -) { - for (ts_selection, limits) in asset.iter_activity_per_capacity_limits() { - let mut upper_limit = limits.end().value(); - let mut lower_limit = limits.start().value(); - - // If the asset capacity is discrete, the capacity variable represents number of - // units, so we need to multiply the per-capacity limits by the unit size. - if let AssetCapacity::Discrete(_, unit_size) = asset.capacity() { - upper_limit *= unit_size.value(); - lower_limit *= unit_size.value(); - } - - // Collect capacity and activity terms - // We have a single capacity term, and activity terms for all time slices in the selection - let mut terms_upper = vec![(capacity_var, -upper_limit)]; - let mut terms_lower = vec![(capacity_var, -lower_limit)]; - for (time_slice, _) in ts_selection.iter(time_slice_info) { - let var = *activity_vars.get(time_slice).unwrap(); - terms_upper.push((var, 1.0)); - terms_lower.push((var, 1.0)); - } - - // Upper bound: sum(activity) - (capacity * upper_limit_per_capacity) ≤ 0 - problem.add_row(..=0.0, &terms_upper); - - // Lower bound: sum(activity) - (capacity * lower_limit_per_capacity) ≥ 0 - problem.add_row(0.0.., &terms_lower); - } -} - /// Adds demand constraints to the problem. /// /// Constrains supply to be less than or equal to demand. This is implemented as an equality diff --git a/src/simulation/investment/appraisal/optimisation.rs b/src/simulation/investment/appraisal/optimisation.rs index 0cec3f760..075374dc0 100644 --- a/src/simulation/investment/appraisal/optimisation.rs +++ b/src/simulation/investment/appraisal/optimisation.rs @@ -1,17 +1,16 @@ //! Optimisation problem for investment tools. use super::DemandMap; use super::ObjectiveCoefficients; -use super::constraints::{ - add_activity_constraints, add_capacity_constraint, add_demand_constraints, -}; -use crate::asset::{AssetCapacity, AssetRef}; +use super::constraints::{add_activity_constraints, add_demand_constraints}; +use crate::asset::AssetCapacity; +use crate::asset::AssetRef; use crate::commodity::Commodity; use crate::model::Model; use crate::simulation::optimisation::ModelError; use crate::simulation::optimisation::apply_highs_options_from_toml; use crate::simulation::optimisation::solve_optimal; use crate::time_slice::{TimeSliceID, TimeSliceInfo}; -use crate::units::{Activity, Capacity, Flow}; +use crate::units::{Activity, Flow}; use anyhow::{Context, Result}; use highs::{RowProblem as Problem, Sense}; use indexmap::IndexMap; @@ -24,11 +23,6 @@ pub type Variable = highs::Col; /// Map storing variables for the optimisation problem struct VariableMap { - /// Capacity variable. - /// - /// This represents absolute capacity for indivisible assets and number of units for - /// divisible assets. - capacity_var: Variable, /// Activity variables in each time slice activity_vars: IndexMap, /// Unmet demand variables @@ -44,24 +38,7 @@ impl VariableMap { /// /// # Returns /// A new `VariableMap` containing all created decision variables - fn add_to_problem( - problem: &mut Problem, - cost_coefficients: &ObjectiveCoefficients, - capacity_unit_size: Option, - ) -> Self { - // Create capacity variable with its associated cost - let capacity_coefficient = cost_coefficients.capacity_coefficient.value(); - let capacity_var = match capacity_unit_size { - Some(unit_size) => { - // Divisible asset: capacity variable represents number of units - problem.add_integer_column(capacity_coefficient * unit_size.value(), 0.0..) - } - None => { - // Indivisible asset: capacity variable represents total capacity - problem.add_column(capacity_coefficient, 0.0..) - } - }; - + fn add_to_problem(problem: &mut Problem, cost_coefficients: &ObjectiveCoefficients) -> Self { // Create activity variables for each time slice let mut activity_vars = IndexMap::new(); for (time_slice, cost) in &cost_coefficients.activity_coefficients { @@ -72,12 +49,11 @@ impl VariableMap { // Create unmet demand variables for each time slice let mut unmet_demand_vars = IndexMap::new(); for time_slice in cost_coefficients.activity_coefficients.keys() { - let var = problem.add_column(cost_coefficients.unmet_demand_coefficient.value(), 0.0..); + let var = problem.add_column(0.0, 0.0..); unmet_demand_vars.insert(time_slice.clone(), var); } Self { - capacity_var, activity_vars, unmet_demand_vars, } @@ -86,8 +62,6 @@ impl VariableMap { /// Map containing optimisation results and coefficients pub struct ResultsMap { - /// Capacity variable - pub capacity: AssetCapacity, /// Activity variables in each time slice pub activity: IndexMap, /// Unmet demand variables @@ -104,11 +78,10 @@ fn add_constraints( demand: &DemandMap, time_slice_info: &TimeSliceInfo, ) { - add_capacity_constraint(problem, asset, max_capacity, variables.capacity_var); add_activity_constraints( problem, asset, - variables.capacity_var, + max_capacity, &variables.activity_vars, time_slice_info, ); @@ -135,11 +108,10 @@ pub fn perform_optimisation( commodity: &Commodity, coefficients: &ObjectiveCoefficients, demand: &DemandMap, - sense: Sense, ) -> Result { // Create problem and add variables let mut problem = Problem::default(); - let variables = VariableMap::add_to_problem(&mut problem, coefficients, asset.unit_size()); + let variables = VariableMap::add_to_problem(&mut problem, coefficients); // Add constraints add_constraints( @@ -153,7 +125,7 @@ pub fn perform_optimisation( ); // Solve model - let mut highs_model = problem.optimise(sense); + let mut highs_model = problem.optimise(Sense::Maximise); apply_highs_options_from_toml(&mut highs_model, &model.parameters.highs.appraisal_options) .context("Failed to apply custom HiGHS options to appraisal optimisation")?; let solution = solve_optimal(highs_model) @@ -161,28 +133,19 @@ pub fn perform_optimisation( .get_solution(); let solution_values = solution.columns(); Ok(ResultsMap { - // If the asset has a defined unit size, the capacity variable represents number of units, - // otherwise it represents absolute capacity - #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] - capacity: match asset.unit_size() { - Some(unit_size) => { - AssetCapacity::Discrete(solution_values[0].round() as u32, unit_size) - } - None => AssetCapacity::Continuous(Capacity::new(solution_values[0])), - }, // The mapping below assumes the column ordering documented on `VariableMap::add_to_problem`: - // index 0 = capacity, next `n` entries = activities (in the same key order as + // first `n` entries = activities (in the same key order as // `cost_coefficients.activity_coefficients`), remaining entries = unmet demand. activity: variables .activity_vars .keys() - .zip(solution_values[1..].iter()) + .zip(solution_values.iter()) .map(|(time_slice, &value)| (time_slice.clone(), Activity::new(value))) .collect(), unmet_demand: variables .unmet_demand_vars .keys() - .zip(solution_values[variables.activity_vars.len() + 1..].iter()) + .zip(solution_values[variables.activity_vars.len()..].iter()) .map(|(time_slice, &value)| (time_slice.clone(), Flow::new(value))) .collect(), })