1use crate::asset::{Asset, AssetCapacity, AssetRef, AssetState};
5use crate::commodity::CommodityID;
6use crate::finance::annual_capital_cost;
7use crate::input::format_items_with_cap;
8use crate::model::Model;
9use crate::output::DataWriter;
10use crate::region::RegionID;
11use crate::simulation::CommodityPrices;
12use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection};
13use crate::units::{
14 Activity, Capacity, Flow, Money, MoneyPerActivity, MoneyPerCapacity, MoneyPerFlow, Year,
15};
16use anyhow::{Result, bail, ensure};
17use highs::{HighsModelStatus, HighsStatus, RowProblem as Problem, Sense};
18use indexmap::{IndexMap, IndexSet};
19use itertools::{chain, iproduct};
20use std::error::Error;
21use std::fmt;
22use std::ops::Range;
23
24mod constraints;
25use constraints::{ConstraintKeys, add_model_constraints};
26
27pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), Flow>;
29
30type Variable = highs::Col;
35
36type ActivityVariableMap = IndexMap<(AssetRef, TimeSliceID), Variable>;
38
39type CapacityVariableMap = IndexMap<AssetRef, Variable>;
41
42type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Variable>;
44
45pub struct VariableMap {
55 activity_vars: ActivityVariableMap,
56 existing_asset_var_idx: Range<usize>,
57 candidate_asset_var_idx: Range<usize>,
58 capacity_vars: CapacityVariableMap,
59 capacity_var_idx: Range<usize>,
60 unmet_demand_vars: UnmetDemandVariableMap,
61 unmet_demand_var_idx: Range<usize>,
62}
63
64impl VariableMap {
65 fn new_with_activity_vars(
76 problem: &mut Problem,
77 model: &Model,
78 input_prices: Option<&CommodityPrices>,
79 existing_assets: &[AssetRef],
80 candidate_assets: &[AssetRef],
81 year: u32,
82 ) -> Self {
83 let mut activity_vars = ActivityVariableMap::new();
84 let existing_asset_var_idx = add_activity_variables(
85 problem,
86 &mut activity_vars,
87 &model.time_slice_info,
88 input_prices,
89 existing_assets,
90 year,
91 );
92 let candidate_asset_var_idx = add_activity_variables(
93 problem,
94 &mut activity_vars,
95 &model.time_slice_info,
96 input_prices,
97 candidate_assets,
98 year,
99 );
100
101 Self {
102 activity_vars,
103 existing_asset_var_idx,
104 candidate_asset_var_idx,
105 capacity_vars: CapacityVariableMap::new(),
106 capacity_var_idx: Range::default(),
107 unmet_demand_vars: UnmetDemandVariableMap::default(),
108 unmet_demand_var_idx: Range::default(),
109 }
110 }
111
112 fn add_unmet_demand_variables(
120 &mut self,
121 problem: &mut Problem,
122 model: &Model,
123 markets_to_allow_unmet_demand: &[(CommodityID, RegionID)],
124 ) {
125 assert!(!markets_to_allow_unmet_demand.is_empty());
126
127 let start = problem.num_cols();
129
130 let voll = model.parameters.value_of_lost_load;
132 self.unmet_demand_vars.extend(
133 iproduct!(
134 markets_to_allow_unmet_demand.iter(),
135 model.time_slice_info.iter_ids()
136 )
137 .map(|((commodity_id, region_id), time_slice)| {
138 let key = (commodity_id.clone(), region_id.clone(), time_slice.clone());
139 let var = problem.add_column(voll.value(), 0.0..);
140 (key, var)
141 }),
142 );
143
144 self.unmet_demand_var_idx = start..problem.num_cols();
145 }
146
147 fn get_activity_var(&self, asset: &AssetRef, time_slice: &TimeSliceID) -> Variable {
149 let key = (asset.clone(), time_slice.clone());
150
151 *self
152 .activity_vars
153 .get(&key)
154 .expect("No asset variable found for given params")
155 }
156
157 fn get_unmet_demand_var(
159 &self,
160 commodity_id: &CommodityID,
161 region_id: &RegionID,
162 time_slice: &TimeSliceID,
163 ) -> Variable {
164 *self
165 .unmet_demand_vars
166 .get(&(commodity_id.clone(), region_id.clone(), time_slice.clone()))
167 .expect("No unmet demand variable for given params")
168 }
169
170 fn activity_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> {
172 self.activity_vars.keys()
173 }
174
175 fn iter_capacity_vars(&self) -> impl Iterator<Item = (&AssetRef, Variable)> {
177 self.capacity_vars.iter().map(|(asset, var)| (asset, *var))
178 }
179}
180
181#[allow(clippy::struct_field_names)]
183pub struct Solution<'a> {
184 solution: highs::Solution,
185 variables: VariableMap,
186 time_slice_info: &'a TimeSliceInfo,
187 constraint_keys: ConstraintKeys,
188 pub objective_value: Money,
190}
191
192impl Solution<'_> {
193 pub fn create_flow_map(&self) -> FlowMap {
198 let mut flows = FlowMap::new();
201 for (asset, time_slice, activity) in self.iter_activity_for_existing() {
202 for flow in asset.iter_flows() {
203 let flow_key = (asset.clone(), flow.commodity.id.clone(), time_slice.clone());
204 let flow_value = activity * flow.coeff;
205 flows.insert(flow_key, flow_value);
206 }
207 }
208
209 flows
210 }
211
212 pub fn iter_activity(&self) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
214 self.variables
215 .activity_var_keys()
216 .zip(self.solution.columns())
217 .map(|((asset, time_slice), activity)| (asset, time_slice, Activity(*activity)))
218 }
219
220 pub fn iter_activity_for_existing(
222 &self,
223 ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
224 let cols = &self.solution.columns()[self.variables.existing_asset_var_idx.clone()];
225 self.variables
226 .activity_var_keys()
227 .skip(self.variables.existing_asset_var_idx.start)
228 .zip(cols.iter())
229 .map(|((asset, time_slice), &value)| (asset, time_slice, Activity(value)))
230 }
231
232 pub fn iter_activity_for_candidates(
234 &self,
235 ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, Activity)> {
236 let cols = &self.solution.columns()[self.variables.candidate_asset_var_idx.clone()];
237 self.variables
238 .activity_var_keys()
239 .skip(self.variables.candidate_asset_var_idx.start)
240 .zip(cols.iter())
241 .map(|((asset, time_slice), &value)| (asset, time_slice, Activity(value)))
242 }
243
244 pub fn iter_unmet_demand(
246 &self,
247 ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, Flow)> {
248 self.variables
249 .unmet_demand_vars
250 .keys()
251 .zip(self.solution.columns()[self.variables.unmet_demand_var_idx.clone()].iter())
252 .map(|((commodity_id, region_id, time_slice), flow)| {
253 (commodity_id, region_id, time_slice, Flow(*flow))
254 })
255 }
256
257 pub fn iter_capacity(&self) -> impl Iterator<Item = (&AssetRef, AssetCapacity)> {
262 self.variables
263 .capacity_vars
264 .keys()
265 .zip(self.solution.columns()[self.variables.capacity_var_idx.clone()].iter())
266 .map(|(asset, capacity_var)| {
267 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
270 let asset_capacity = if let Some(unit_size) = asset.unit_size() {
271 AssetCapacity::Discrete(capacity_var.round() as u32, unit_size)
272 } else {
273 AssetCapacity::Continuous(Capacity(*capacity_var))
274 };
275 (asset, asset_capacity)
276 })
277 }
278
279 pub fn iter_commodity_balance_duals(
281 &self,
282 ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, MoneyPerFlow)> {
283 self.constraint_keys
287 .commodity_balance_keys
288 .zip_duals(self.solution.dual_rows())
289 .flat_map(|((commodity_id, region_id, ts_selection), price)| {
290 ts_selection
291 .iter(self.time_slice_info)
292 .map(move |(ts, _)| (commodity_id, region_id, ts, price))
293 })
294 }
295
296 pub fn iter_activity_duals(
305 &self,
306 ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
307 self.constraint_keys
308 .activity_keys
309 .zip_duals(self.solution.dual_rows())
310 .filter(|&((_asset, ts_selection), _dual)| {
311 matches!(ts_selection, TimeSliceSelection::Single(_))
312 })
313 .map(|((asset, ts_selection), dual)| {
314 let (time_slice, _) = ts_selection.iter(self.time_slice_info).next().unwrap();
316 (asset, time_slice, dual)
317 })
318 }
319
320 pub fn iter_column_duals(
322 &self,
323 ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID, MoneyPerActivity)> {
324 self.variables
325 .activity_var_keys()
326 .zip(self.solution.dual_columns())
327 .map(|((asset, time_slice), dual)| (asset, time_slice, MoneyPerActivity(*dual)))
328 }
329}
330
331#[derive(Debug, Clone)]
333pub enum ModelError {
334 Incoherent(HighsStatus),
338 NonOptimal(HighsModelStatus),
340}
341
342impl fmt::Display for ModelError {
343 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
344 match self {
345 ModelError::Incoherent(status) => write!(f, "Incoherent model: {status:?}"),
346 ModelError::NonOptimal(status) => {
347 write!(f, "Could not find optimal result: {status:?}")
348 }
349 }
350 }
351}
352
353impl Error for ModelError {}
354
355pub fn solve_optimal(model: highs::Model) -> Result<highs::SolvedModel, ModelError> {
357 let solved = model.try_solve().map_err(ModelError::Incoherent)?;
358
359 match solved.status() {
360 HighsModelStatus::Optimal => Ok(solved),
361 status => Err(ModelError::NonOptimal(status)),
362 }
363}
364
365fn filter_input_prices(
370 input_prices: &CommodityPrices,
371 markets_to_balance: &[(CommodityID, RegionID)],
372) -> CommodityPrices {
373 input_prices
374 .iter()
375 .filter(|(commodity_id, region_id, _, _)| {
376 !markets_to_balance
377 .iter()
378 .any(|(c, r)| c == *commodity_id && r == *region_id)
379 })
380 .collect()
381}
382
383pub struct DispatchRun<'model, 'run> {
393 model: &'model Model,
394 existing_assets: &'run [AssetRef],
395 flexible_capacity_assets: &'run [AssetRef],
396 candidate_assets: &'run [AssetRef],
397 markets_to_balance: &'run [(CommodityID, RegionID)],
398 input_prices: Option<&'run CommodityPrices>,
399 year: u32,
400 capacity_margin: f64,
401}
402
403impl<'model, 'run> DispatchRun<'model, 'run> {
404 pub fn new(model: &'model Model, assets: &'run [AssetRef], year: u32) -> Self {
406 Self {
407 model,
408 existing_assets: assets,
409 flexible_capacity_assets: &[],
410 candidate_assets: &[],
411 markets_to_balance: &[],
412 input_prices: None,
413 year,
414 capacity_margin: 0.0,
415 }
416 }
417
418 pub fn with_flexible_capacity_assets(
420 self,
421 flexible_capacity_assets: &'run [AssetRef],
422 capacity_margin: f64,
423 ) -> Self {
424 Self {
425 flexible_capacity_assets,
426 capacity_margin,
427 ..self
428 }
429 }
430
431 pub fn with_candidates(self, candidate_assets: &'run [AssetRef]) -> Self {
433 Self {
434 candidate_assets,
435 ..self
436 }
437 }
438
439 pub fn with_market_balance_subset(
441 self,
442 markets_to_balance: &'run [(CommodityID, RegionID)],
443 ) -> Self {
444 assert!(!markets_to_balance.is_empty());
445
446 Self {
447 markets_to_balance,
448 ..self
449 }
450 }
451
452 pub fn with_input_prices(self, input_prices: &'run CommodityPrices) -> Self {
454 Self {
455 input_prices: Some(input_prices),
456 ..self
457 }
458 }
459
460 pub fn run(&self, run_description: &str, writer: &mut DataWriter) -> Result<Solution<'model>> {
472 let all_markets: Vec<_>;
474 let markets_to_balance = if self.markets_to_balance.is_empty() {
475 all_markets = self.model.iter_markets().collect();
476 &all_markets
477 } else {
478 self.markets_to_balance
479 };
480
481 let input_prices_owned = self
483 .input_prices
484 .map(|prices| filter_input_prices(prices, markets_to_balance));
485 let input_prices = input_prices_owned.as_ref();
486
487 match self.run_without_unmet_demand_variables(markets_to_balance, input_prices) {
491 Ok(solution) => {
492 writer.write_dispatch_debug_info(self.year, run_description, &solution)?;
494 Ok(solution)
495 }
496 Err(ModelError::NonOptimal(HighsModelStatus::Infeasible)) => {
497 let solution = self
500 .run_internal(
501 markets_to_balance,
502 true,
503 input_prices,
504 )
505 .expect("Failed to run dispatch to calculate unmet demand");
506
507 writer.write_dispatch_debug_info(self.year, run_description, &solution)?;
509
510 let markets: IndexSet<_> = solution
512 .iter_unmet_demand()
513 .filter(|(_, _, _, flow)| *flow > Flow(0.0))
514 .map(|(commodity_id, region_id, _, _)| {
515 (commodity_id.clone(), region_id.clone())
516 })
517 .collect();
518
519 ensure!(
520 !markets.is_empty(),
521 "Model is infeasible, but there was no unmet demand"
522 );
523
524 bail!(
525 "The solver has indicated that the problem is infeasible, probably because \
526 the supplied assets could not meet the required demand. Demand was not met \
527 for the following markets: {}",
528 format_items_with_cap(markets)
529 )
530 }
531 Err(err) => Err(err)?,
532 }
533 }
534
535 fn run_without_unmet_demand_variables(
537 &self,
538 markets_to_balance: &[(CommodityID, RegionID)],
539 input_prices: Option<&CommodityPrices>,
540 ) -> Result<Solution<'model>, ModelError> {
541 self.run_internal(
542 markets_to_balance,
543 false,
544 input_prices,
545 )
546 }
547
548 fn run_internal(
550 &self,
551 markets_to_balance: &[(CommodityID, RegionID)],
552 allow_unmet_demand: bool,
553 input_prices: Option<&CommodityPrices>,
554 ) -> Result<Solution<'model>, ModelError> {
555 let mut problem = Problem::default();
557 let mut variables = VariableMap::new_with_activity_vars(
558 &mut problem,
559 self.model,
560 input_prices,
561 self.existing_assets,
562 self.candidate_assets,
563 self.year,
564 );
565
566 if allow_unmet_demand {
569 variables.add_unmet_demand_variables(&mut problem, self.model, markets_to_balance);
570 }
571
572 for asset in self.flexible_capacity_assets {
574 assert!(
575 self.existing_assets.contains(asset),
576 "Flexible capacity assets must be a subset of existing assets. Offending asset: {asset:?}"
577 );
578 }
579
580 if !self.flexible_capacity_assets.is_empty() {
582 variables.capacity_var_idx = add_capacity_variables(
583 &mut problem,
584 &mut variables.capacity_vars,
585 self.flexible_capacity_assets,
586 self.capacity_margin,
587 );
588 }
589
590 let all_assets = chain(self.existing_assets.iter(), self.candidate_assets.iter());
592 let constraint_keys = add_model_constraints(
593 &mut problem,
594 &variables,
595 self.model,
596 &all_assets,
597 markets_to_balance,
598 self.year,
599 );
600
601 let solution = solve_optimal(problem.optimise(Sense::Minimise))?;
603
604 Ok(Solution {
605 solution: solution.get_solution(),
606 variables,
607 time_slice_info: &self.model.time_slice_info,
608 constraint_keys,
609 objective_value: Money(solution.objective_value()),
610 })
611 }
612}
613
614fn add_activity_variables(
625 problem: &mut Problem,
626 variables: &mut ActivityVariableMap,
627 time_slice_info: &TimeSliceInfo,
628 input_prices: Option<&CommodityPrices>,
629 assets: &[AssetRef],
630 year: u32,
631) -> Range<usize> {
632 let start = problem.num_cols();
634
635 for (asset, time_slice) in iproduct!(assets.iter(), time_slice_info.iter_ids()) {
636 let coeff = calculate_activity_coefficient(asset, year, time_slice, input_prices);
637 let var = problem.add_column(coeff.value(), 0.0..);
638 let key = (asset.clone(), time_slice.clone());
639 let existing = variables.insert(key, var).is_some();
640 assert!(!existing, "Duplicate entry for var");
641 }
642
643 start..problem.num_cols()
644}
645
646fn add_capacity_variables(
647 problem: &mut Problem,
648 variables: &mut CapacityVariableMap,
649 assets: &[AssetRef],
650 capacity_margin: f64,
651) -> Range<usize> {
652 let start = problem.num_cols();
654
655 for asset in assets {
656 assert!(
658 matches!(asset.state(), AssetState::Selected { .. }),
659 "Flexible capacity can only be assigned to `Selected` type assets. Offending asset: {asset:?}"
660 );
661
662 let current_capacity = asset.capacity();
663 let coeff = calculate_capacity_coefficient(asset);
664
665 let var = match current_capacity {
666 AssetCapacity::Continuous(cap) => {
667 let lower = ((1.0 - capacity_margin) * cap.value()).max(0.0);
669 let upper = (1.0 + capacity_margin) * cap.value();
670 problem.add_column(coeff.value(), lower..=upper)
671 }
672 AssetCapacity::Discrete(units, unit_size) => {
673 let lower = ((1.0 - capacity_margin) * units as f64).max(0.0);
675 let upper = (1.0 + capacity_margin) * units as f64;
676 problem.add_integer_column((coeff * unit_size).value(), lower..=upper)
677 }
678 };
679
680 let existing = variables.insert(asset.clone(), var).is_some();
681 assert!(!existing, "Duplicate entry for var");
682 }
683
684 start..problem.num_cols()
685}
686
687fn calculate_activity_coefficient(
704 asset: &Asset,
705 year: u32,
706 time_slice: &TimeSliceID,
707 input_prices: Option<&CommodityPrices>,
708) -> MoneyPerActivity {
709 let opex = asset.get_operating_cost(year, time_slice);
710 let input_cost = input_prices
711 .map(|prices| asset.get_input_cost_from_prices(prices, time_slice))
712 .unwrap_or_default();
713 opex + input_cost
714}
715
716fn calculate_capacity_coefficient(asset: &AssetRef) -> MoneyPerCapacity {
720 let param = asset.process_parameter();
721 let annual_fixed_operating_cost = param.fixed_operating_cost * Year(1.0);
722 annual_fixed_operating_cost
723 + annual_capital_cost(param.capital_cost, param.lifetime, param.discount_rate)
724}