muse2/simulation/
optimisation.rs

1//! Code for performing dispatch optimisation.
2//!
3//! This is used to calculate commodity flows and prices.
4use 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
27/// A map of commodity flows calculated during the optimisation
28pub type FlowMap = IndexMap<(AssetRef, CommodityID, TimeSliceID), Flow>;
29
30/// A decision variable in the optimisation
31///
32/// Note that this type does **not** include the value of the variable; it just refers to a
33/// particular column of the problem.
34type Variable = highs::Col;
35
36/// The map of activity variables for assets
37type ActivityVariableMap = IndexMap<(AssetRef, TimeSliceID), Variable>;
38
39/// A map of capacity variables for assets
40type CapacityVariableMap = IndexMap<AssetRef, Variable>;
41
42/// Variables representing unmet demand for a given market
43type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Variable>;
44
45/// A map for easy lookup of variables in the problem.
46///
47/// The entries are ordered (see [`IndexMap`]).
48///
49/// We use this data structure for two things:
50///
51/// 1. In order define constraints for the optimisation
52/// 2. To keep track of the combination of parameters that each variable corresponds to, for when we
53///    are reading the results of the optimisation.
54pub 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    /// Create a new [`VariableMap`] and add activity variables to the problem
66    ///
67    /// # Arguments
68    ///
69    /// * `problem` - The optimisation problem
70    /// * `model` - The model
71    /// * `input_prices` - Optional explicit prices for input commodities
72    /// * `existing_assets` - The asset pool
73    /// * `candidate_assets` - Candidate assets for inclusion in active pool
74    /// * `year` - Current milestone year
75    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    /// Add unmet demand variables to the map and the problem
113    ///
114    /// # Arguments
115    ///
116    /// * `problem` - The optimisation problem
117    /// * `model` - The model
118    /// * `markets_to_allow_unmet_demand` - The subset of markets to add unmet demand variables for
119    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        // This line **must** come before we add more variables
128        let start = problem.num_cols();
129
130        // Add variables
131        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    /// Get the activity [`Variable`] corresponding to the given parameters.
148    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    /// Get the unmet demand [`Variable`] corresponding to the given parameters.
158    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    /// Iterate over the keys for activity variables
171    fn activity_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> {
172        self.activity_vars.keys()
173    }
174
175    /// Iterate over capacity variables
176    fn iter_capacity_vars(&self) -> impl Iterator<Item = (&AssetRef, Variable)> {
177        self.capacity_vars.iter().map(|(asset, var)| (asset, *var))
178    }
179}
180
181/// The solution to the dispatch optimisation problem
182#[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    /// The objective value for the solution
189    pub objective_value: Money,
190}
191
192impl Solution<'_> {
193    /// Create a map of commodity flows for each asset's coeffs at every time slice.
194    ///
195    /// Note that this only includes commodity flows which relate to existing assets, so not every
196    /// commodity in the simulation will necessarily be represented.
197    pub fn create_flow_map(&self) -> FlowMap {
198        // The decision variables represent assets' activity levels, not commodity flows. We
199        // multiply this value by the flow coeffs to get commodity flows.
200        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    /// Activity for all assets (existing and candidate, if present)
213    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    /// Activity for each existing asset
221    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    /// Activity for each candidate asset
233    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    /// Iterate over unmet demand
245    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    /// Iterate over capacity values
258    ///
259    /// Will return `AssetCapacity::Continuous` or `AssetCapacity::Discrete` depending on whether
260    /// the asset has a defined unit size.
261    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                // If the asset has a defined unit size, the capacity variable represents number of
268                // units, otherwise it represents absolute capacity
269                #[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    /// Keys and dual values for commodity balance constraints.
280    pub fn iter_commodity_balance_duals(
281        &self,
282    ) -> impl Iterator<Item = (&CommodityID, &RegionID, &TimeSliceID, MoneyPerFlow)> {
283        // Each commodity balance constraint applies to a particular time slice
284        // selection (depending on time slice level). Where this covers multiple time slices,
285        // we return the same dual for each individual time slice.
286        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    /// Keys and dual values for activity constraints.
297    ///
298    /// Note: if there are any flexible capacity assets, these will have two duals with identical
299    /// keys, and there will be no way to distinguish between them in the resulting iterator.
300    /// Recommended for now only to use this function when there are no flexible capacity assets.
301    ///
302    /// Also note: this excludes seasonal and annual constraints. Recommended for now not to use
303    /// this for models that include seasonal or annual availability constraints.
304    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                // `unwrap` is safe here because we just matched Single(_)
315                let (time_slice, _) = ts_selection.iter(self.time_slice_info).next().unwrap();
316                (asset, time_slice, dual)
317            })
318    }
319
320    /// Keys and values for column duals.
321    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/// Defines the possible errors that can occur when running the solver
332#[derive(Debug, Clone)]
333pub enum ModelError {
334    /// The model definition is incoherent.
335    ///
336    /// Users should not be able to trigger this error.
337    Incoherent(HighsStatus),
338    /// An optimal solution could not be found
339    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
355/// Try to solve the model, returning an error if the model is incoherent or result is non-optimal
356pub 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
365/// Filter prices data to only include prices for markets not being balanced
366///
367/// Markets being balanced (i.e. with commodity balance constraints) will have prices calculated
368/// internally by the solver, so we need to remove them to prevent double-counting.
369fn 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
383/// Provides the interface for running the dispatch optimisation.
384///
385/// The run will attempt to meet unmet demand: if the solver reports infeasibility
386/// the implementation will rerun including unmet-demand variables to identify offending
387/// markets and provide a clearer error message.
388///
389/// For a detailed description, please see the [dispatch optimisation formulation][1].
390///
391/// [1]: https://energysystemsmodellinglab.github.io/MUSE2/model/dispatch_optimisation.html
392pub 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    /// Create a new [`DispatchRun`] for the specified model and assets for a given year
405    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    /// Include the specified flexible capacity assets in the dispatch run
419    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    /// Include the specified candidate assets in the dispatch run
432    pub fn with_candidates(self, candidate_assets: &'run [AssetRef]) -> Self {
433        Self {
434            candidate_assets,
435            ..self
436        }
437    }
438
439    /// Only apply commodity balance constraints to the specified subset of markets
440    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    /// Explicitly provide prices for certain input commodities
453    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    /// Perform the dispatch optimisation.
461    ///
462    /// # Arguments
463    ///
464    /// * `run_description` - Which dispatch run for the current year this is
465    /// * `writer` - For saving output data
466    ///
467    /// # Returns
468    ///
469    /// A solution containing new commodity flows for assets and prices for (some) commodities or an
470    /// error.
471    pub fn run(&self, run_description: &str, writer: &mut DataWriter) -> Result<Solution<'model>> {
472        // If the user provided no markets to balance, we use all of them
473        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        // Select prices for markets not being balanced
482        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        // Try running dispatch. If it fails because the model is infeasible, it is likely that this
488        // is due to unmet demand, in this case, we rerun dispatch including extra variables to
489        // track the unmet demand so we can report the offending markets to users
490        match self.run_without_unmet_demand_variables(markets_to_balance, input_prices) {
491            Ok(solution) => {
492                // Normal successful run: write debug info and return
493                writer.write_dispatch_debug_info(self.year, run_description, &solution)?;
494                Ok(solution)
495            }
496            Err(ModelError::NonOptimal(HighsModelStatus::Infeasible)) => {
497                // Re-run including unmet demand variables so we can record detailed unmet-demand
498                // debug output before returning an error to the caller.
499                let solution = self
500                    .run_internal(
501                        markets_to_balance,
502                        /*allow_unmet_demand=*/ true,
503                        input_prices,
504                    )
505                    .expect("Failed to run dispatch to calculate unmet demand");
506
507                // Write debug CSVs to help diagnosis
508                writer.write_dispatch_debug_info(self.year, run_description, &solution)?;
509
510                // Collect markets with unmet demand from the solution
511                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    /// Run dispatch without unmet demand variables
536    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            /*allow_unmet_demand=*/ false,
544            input_prices,
545        )
546    }
547
548    /// Run dispatch to balance the specified markets, optionally including unmet demand variables
549    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        // Set up problem
556        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 unmet demand is enabled for this dispatch run (and is allowed by the model param) then
567        // we add variables representing unmet demand for all markets being balanced
568        if allow_unmet_demand {
569            variables.add_unmet_demand_variables(&mut problem, self.model, markets_to_balance);
570        }
571
572        // Check flexible capacity assets is a subset of existing assets
573        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        // Add capacity variables for flexible capacity assets
581        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        // Add constraints
591        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        // Solve model
602        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
614/// Add variables to the optimisation problem.
615///
616/// # Arguments
617///
618/// * `problem` - The optimisation problem
619/// * `variables` - The map of asset variables
620/// * `time_slice_info` - Information about assets
621/// * `input_prices` - Optional explicit prices for input commodities
622/// * `assets` - Assets to include
623/// * `year` - Current milestone year
624fn 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    // This line **must** come before we add more variables
633    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    // This line **must** come before we add more variables
653    let start = problem.num_cols();
654
655    for asset in assets {
656        // Can only have flexible capacity for `Selected` assets
657        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                // Continuous capacity: capacity variable represents total capacity
668                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                // Discrete capacity: capacity variable represents number of units
674                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
687/// Calculate the cost coefficient for an activity variable.
688///
689/// Normally, the cost coefficient is the same as the asset's operating costs for the given year and
690/// time slice. If `input_prices` is provided then those prices are added to the flow costs for the
691/// relevant commodities, if they are input flows for the asset.
692///
693/// # Arguments
694///
695/// * `asset` - The asset to calculate the coefficient for
696/// * `year` - The current milestone year
697/// * `time_slice` - The time slice to which this coefficient applies
698/// * `input_prices` - Optional map of prices to include for input commodities
699///
700/// # Returns
701///
702/// The cost coefficient to be used for the relevant decision variable.
703fn 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
716/// Calculate the cost coefficient for a capacity variable (for flexible capacity assets only).
717///
718/// This includes both the annual fixed operating cost and the annual capital cost.
719fn 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}