muse2/simulation/
investment.rs

1//! Code for performing agent investment.
2use super::optimisation::{DispatchRun, FlowMap};
3use super::prices::ReducedCosts;
4use crate::agent::Agent;
5use crate::asset::{Asset, AssetIterator, AssetRef, AssetState};
6use crate::commodity::{Commodity, CommodityID, CommodityMap};
7use crate::model::Model;
8use crate::output::DataWriter;
9use crate::region::RegionID;
10use crate::simulation::CommodityPrices;
11use crate::time_slice::{TimeSliceID, TimeSliceInfo};
12use crate::units::{Capacity, Dimensionless, Flow, FlowPerCapacity, MoneyPerFlow};
13use anyhow::{Result, ensure};
14use indexmap::IndexMap;
15use itertools::{chain, iproduct};
16use log::debug;
17use std::collections::HashMap;
18
19pub mod appraisal;
20use appraisal::appraise_investment;
21
22/// A map of demand across time slices for a specific commodity and region
23type DemandMap = IndexMap<TimeSliceID, Flow>;
24
25/// Demand for a given combination of commodity, region and time slice
26type AllDemandMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Flow>;
27
28/// Perform agent investment to determine capacity investment of new assets for next milestone year.
29///
30/// # Arguments
31///
32/// * `model` - The model
33/// * `year` - Current milestone year
34/// * `assets` - The asset pool
35/// * `prices` - Commodity prices
36/// * `reduced_costs` - Reduced costs for assets
37/// * `writer` - Data writer
38pub fn perform_agent_investment(
39    model: &Model,
40    year: u32,
41    existing_assets: &[AssetRef],
42    prices: &CommodityPrices,
43    reduced_costs: &ReducedCosts,
44    writer: &mut DataWriter,
45) -> Result<Vec<AssetRef>> {
46    // Initialise demand map
47    let mut demand =
48        flatten_preset_demands_for_year(&model.commodities, &model.time_slice_info, year);
49
50    // Keep a list of all the assets selected
51    // This includes Commissioned assets that are selected for retention, and new Selected assets
52    let mut all_selected_assets = Vec::new();
53
54    for region_id in model.iter_regions() {
55        let cur_commodities = &model.commodity_order[&(region_id.clone(), year)];
56
57        // Prices to be used for input flows for commodities not produced in dispatch run
58        let mut external_prices =
59            get_prices_for_commodities(prices, &model.time_slice_info, region_id, cur_commodities);
60        let mut seen_commodities = Vec::new();
61        for commodity_id in cur_commodities.iter() {
62            seen_commodities.push(commodity_id.clone());
63            let commodity = &model.commodities[commodity_id];
64
65            // Remove prices for already-seen commodities. Commodities which are produced by at
66            // least one asset in the dispatch run will have prices produced endogenously (via the
67            // commodity balance constraints), but commodities for which investment has not yet been
68            // performed will, by definition, not have any producers. For these, we provide prices
69            // from the previous dispatch run otherwise they will appear to be free to the model.
70            for time_slice in model.time_slice_info.iter_ids() {
71                external_prices.remove(&(
72                    commodity_id.clone(),
73                    region_id.clone(),
74                    time_slice.clone(),
75                ));
76            }
77
78            // List of assets selected/retained for this region/commodity
79            let mut selected_assets = Vec::new();
80
81            for (agent, commodity_portion) in
82                get_responsible_agents(model.agents.values(), commodity_id, region_id, year)
83            {
84                debug!(
85                    "Running investment for agent '{}' with commodity '{}' in region '{}'",
86                    &agent.id, commodity_id, region_id
87                );
88
89                // Get demand portion for this commodity for this agent in this region/year
90                let demand_portion_for_commodity = get_demand_portion_for_commodity(
91                    &model.time_slice_info,
92                    &demand,
93                    commodity_id,
94                    region_id,
95                    commodity_portion,
96                );
97
98                // Existing and candidate assets from which to choose
99                let opt_assets = get_asset_options(
100                    &model.time_slice_info,
101                    existing_assets,
102                    &demand_portion_for_commodity,
103                    agent,
104                    commodity,
105                    region_id,
106                    year,
107                )
108                .collect();
109
110                // Choose assets from among existing pool and candidates
111                let best_assets = select_best_assets(
112                    model,
113                    opt_assets,
114                    commodity,
115                    agent,
116                    reduced_costs,
117                    demand_portion_for_commodity,
118                    year,
119                    writer,
120                )?;
121                selected_assets.extend(best_assets);
122            }
123
124            // If no assets have been selected for this region/commodity, skip dispatch optimisation
125            // **TODO**: this probably means there's no demand for the commodity, which we could
126            // presumably preempt
127            if selected_assets.is_empty() {
128                continue;
129            }
130
131            // Add the selected assets to the list of all selected assets
132            all_selected_assets.extend(selected_assets.clone());
133
134            // Perform dispatch optimisation with assets that have been selected so far
135            // **TODO**: presumably we only need to do this for selected_assets, as assets added in
136            // previous iterations should not change
137            debug!(
138                "Running post-investment dispatch for commodity '{commodity_id}' in region '{region_id}'"
139            );
140
141            // As upstream commodities by definition will not yet have producers, we explicitly set
142            // their prices using previous values so that they don't appear free
143            let solution = DispatchRun::new(model, &all_selected_assets, year)
144                .with_commodity_subset(&seen_commodities)
145                .with_input_prices(&external_prices)
146                .run(
147                    &format!("post {commodity_id}/{region_id} investment"),
148                    writer,
149                )?;
150
151            // Update demand map with flows from newly added assets
152            update_demand_map(&mut demand, &solution.create_flow_map(), &selected_assets);
153        }
154    }
155
156    Ok(all_selected_assets)
157}
158
159/// Flatten the preset commodity demands for a given year into a map of commodity, region and
160/// time slice to demand.
161///
162/// Since demands for some commodities may be specified at a coarser timeslice level, we need to
163/// distribute these demands over all timeslices. Note: the way that we do this distribution is
164/// irrelevant, as demands will only be balanced to the appropriate level, but we still need to do
165/// this for the solver to work.
166///
167/// **TODO**: these assumptions may need to be revisited, e.g. when we come to storage technologies
168fn flatten_preset_demands_for_year(
169    commodities: &CommodityMap,
170    time_slice_info: &TimeSliceInfo,
171    year: u32,
172) -> AllDemandMap {
173    let mut demand_map = AllDemandMap::new();
174    for (commodity_id, commodity) in commodities.iter() {
175        for ((region_id, data_year, time_slice_selection), demand) in commodity.demand.iter() {
176            if *data_year != year {
177                continue;
178            }
179
180            // We split the demand equally over all timeslices in the selection
181            // NOTE: since demands will only be balanced to the timeslice level of the commodity
182            // it doesn't matter how we do this distribution, only the total matters.
183            let n_timeslices = time_slice_selection.iter(time_slice_info).count() as f64;
184            let demand_per_slice = *demand / Dimensionless(n_timeslices);
185            for (time_slice, _) in time_slice_selection.iter(time_slice_info) {
186                demand_map.insert(
187                    (commodity_id.clone(), region_id.clone(), time_slice.clone()),
188                    demand_per_slice,
189                );
190            }
191        }
192    }
193    demand_map
194}
195
196/// Update demand map with flows from a set of assets
197fn update_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: &[AssetRef]) {
198    for ((asset, commodity_id, time_slice), flow) in flows.iter() {
199        if assets.contains(asset) {
200            let key = (
201                commodity_id.clone(),
202                asset.region_id().clone(),
203                time_slice.clone(),
204            );
205
206            // Note: we use the negative of the flow as input flows are negative in the flow map.
207            demand
208                .entry(key)
209                .and_modify(|value| *value -= *flow)
210                .or_insert(-*flow);
211        }
212    }
213}
214
215/// Get a portion of the demand profile for this commodity and region
216fn get_demand_portion_for_commodity(
217    time_slice_info: &TimeSliceInfo,
218    demand: &AllDemandMap,
219    commodity_id: &CommodityID,
220    region_id: &RegionID,
221    commodity_portion: Dimensionless,
222) -> DemandMap {
223    time_slice_info
224        .iter_ids()
225        .map(|time_slice| {
226            (
227                time_slice.clone(),
228                commodity_portion
229                    * *demand
230                        .get(&(commodity_id.clone(), region_id.clone(), time_slice.clone()))
231                        .unwrap_or(&Flow(0.0)),
232            )
233        })
234        .collect()
235}
236
237/// Get the agents responsible for a given commodity in a given year along with the commodity
238/// portion for which they are responsible
239fn get_responsible_agents<'a, I>(
240    agents: I,
241    commodity_id: &'a CommodityID,
242    region_id: &'a RegionID,
243    year: u32,
244) -> impl Iterator<Item = (&'a Agent, Dimensionless)>
245where
246    I: Iterator<Item = &'a Agent>,
247{
248    agents.filter_map(move |agent| {
249        if !agent.regions.contains(region_id) {
250            return None;
251        }
252        let portion = agent
253            .commodity_portions
254            .get(&(commodity_id.clone(), year))?;
255
256        Some((agent, *portion))
257    })
258}
259
260/// Get the maximum required capacity across time slices
261fn get_demand_limiting_capacity(
262    time_slice_info: &TimeSliceInfo,
263    asset: &Asset,
264    commodity: &Commodity,
265    demand: &DemandMap,
266) -> Capacity {
267    let coeff = asset.get_flow(&commodity.id).unwrap().coeff;
268    let mut capacity = Capacity(0.0);
269
270    for time_slice_selection in time_slice_info.iter_selections_at_level(commodity.time_slice_level)
271    {
272        let demand_for_selection: Flow = time_slice_selection
273            .iter(time_slice_info)
274            .map(|(time_slice, _)| demand[time_slice])
275            .sum();
276
277        // Calculate max capacity required for this time slice selection
278        // For commodities with a coarse timeslice level, we have to allow the possibility that all
279        // of the demand gets served by production in a single timeslice
280        for (time_slice, _) in time_slice_selection.iter(time_slice_info) {
281            let max_flow_per_cap =
282                *asset.get_activity_per_capacity_limits(time_slice).end() * coeff;
283            if max_flow_per_cap != FlowPerCapacity(0.0) {
284                capacity = capacity.max(demand_for_selection / max_flow_per_cap);
285            }
286        }
287    }
288
289    capacity
290}
291
292/// Get options from existing and potential assets for the given parameters
293fn get_asset_options<'a>(
294    time_slice_info: &'a TimeSliceInfo,
295    all_existing_assets: &'a [AssetRef],
296    demand: &'a DemandMap,
297    agent: &'a Agent,
298    commodity: &'a Commodity,
299    region_id: &'a RegionID,
300    year: u32,
301) -> impl Iterator<Item = AssetRef> + 'a {
302    // Get existing assets which produce the commodity of interest
303    let existing_assets = all_existing_assets
304        .iter()
305        .filter_agent(&agent.id)
306        .filter_region(region_id)
307        .filter_primary_producers_of(&commodity.id)
308        .cloned();
309
310    // Get candidates assets which produce the commodity of interest
311    let candidate_assets =
312        get_candidate_assets(time_slice_info, demand, agent, region_id, commodity, year);
313
314    chain(existing_assets, candidate_assets)
315}
316
317/// Get candidate assets which produce a particular commodity for a given agent
318fn get_candidate_assets<'a>(
319    time_slice_info: &'a TimeSliceInfo,
320    demand: &'a DemandMap,
321    agent: &'a Agent,
322    region_id: &'a RegionID,
323    commodity: &'a Commodity,
324    year: u32,
325) -> impl Iterator<Item = AssetRef> + 'a {
326    agent
327        .iter_possible_producers_of(region_id, &commodity.id, year)
328        .map(move |process| {
329            let mut asset =
330                Asset::new_candidate(process.clone(), region_id.clone(), Capacity(0.0), year)
331                    .unwrap();
332            asset.set_capacity(get_demand_limiting_capacity(
333                time_slice_info,
334                &asset,
335                commodity,
336                demand,
337            ));
338
339            asset.into()
340        })
341}
342
343/// Get a map of prices for a subset of commodities
344fn get_prices_for_commodities(
345    prices: &CommodityPrices,
346    time_slice_info: &TimeSliceInfo,
347    region_id: &RegionID,
348    commodities: &[CommodityID],
349) -> HashMap<(CommodityID, RegionID, TimeSliceID), MoneyPerFlow> {
350    iproduct!(commodities.iter(), time_slice_info.iter_ids())
351        .map(|(commodity_id, time_slice)| {
352            let price = prices.get(commodity_id, region_id, time_slice).unwrap();
353            (
354                (commodity_id.clone(), region_id.clone(), time_slice.clone()),
355                price,
356            )
357        })
358        .collect()
359}
360
361/// Get the best assets for meeting demand for the given commodity
362#[allow(clippy::too_many_arguments)]
363fn select_best_assets(
364    model: &Model,
365    mut opt_assets: Vec<AssetRef>,
366    commodity: &Commodity,
367    agent: &Agent,
368    reduced_costs: &ReducedCosts,
369    mut demand: DemandMap,
370    year: u32,
371    writer: &mut DataWriter,
372) -> Result<Vec<AssetRef>> {
373    let mut best_assets: Vec<AssetRef> = Vec::new();
374
375    let mut remaining_candidate_capacity = HashMap::from_iter(
376        opt_assets
377            .iter()
378            .filter(|asset| !asset.is_commissioned())
379            .map(|asset| (asset.clone(), asset.capacity())),
380    );
381
382    let mut round = 0;
383    while is_any_remaining_demand(&demand) {
384        ensure!(
385            !opt_assets.is_empty(),
386            "Failed to meet demand for commodity '{}' with provided assets",
387            &commodity.id
388        );
389
390        // Appraise all options
391        let mut outputs_for_opts = Vec::new();
392        for asset in opt_assets.iter() {
393            let max_capacity = (!asset.is_commissioned()).then(|| {
394                let max_capacity = model.parameters.capacity_limit_factor * asset.capacity();
395                let remaining_capacity = remaining_candidate_capacity[asset];
396                max_capacity.min(remaining_capacity)
397            });
398
399            let output = appraise_investment(
400                model,
401                asset,
402                max_capacity,
403                commodity,
404                &agent.objectives[&year],
405                reduced_costs,
406                &demand,
407            )?;
408
409            // Store the appraisal results if the capacity is positive. If capacity is zero,
410            // this means the asset is infeasible for investment. This can happen if the asset has
411            // zero activity limits for all time slices with demand. This can also happen due to a
412            // known issue with the NPV objective, for which we do not currently have a solution
413            // (see https://github.com/EnergySystemsModellingLab/MUSE_2.0/issues/716).
414            if output.capacity > Capacity(0.0) {
415                outputs_for_opts.push(output);
416            } else {
417                debug!(
418                    "Skipping candidate '{}' with zero capacity",
419                    asset.process_id()
420                );
421            }
422        }
423
424        // Make sure there are some options to consider
425        ensure!(
426            !outputs_for_opts.is_empty(),
427            "No feasible investment options for commodity '{}'",
428            &commodity.id
429        );
430
431        // Save appraisal results
432        writer.write_appraisal_debug_info(
433            year,
434            &format!("{} {} round {}", &commodity.id, &agent.id, round),
435            &outputs_for_opts,
436        )?;
437
438        // Select the best investment option
439        let best_output = outputs_for_opts
440            .into_iter()
441            .min_by(|a, b| a.metric.partial_cmp(&b.metric).unwrap())
442            .unwrap();
443
444        // Log the selected asset
445        let commissioned_txt = match best_output.asset.state() {
446            AssetState::Commissioned { .. } => "existing",
447            AssetState::Candidate => "candidate",
448            _ => panic!("Selected asset should be either Commissioned or Candidate"),
449        };
450        debug!(
451            "Selected {} asset '{}' (capacity: {})",
452            commissioned_txt,
453            &best_output.asset.process_id(),
454            best_output.capacity
455        );
456
457        // Update the assets
458        update_assets(
459            best_output.asset,
460            best_output.capacity,
461            &mut opt_assets,
462            &mut remaining_candidate_capacity,
463            &mut best_assets,
464        );
465
466        demand = best_output.unmet_demand;
467        round += 1;
468    }
469
470    // Convert Candidate assets to Selected
471    // At this point we also assign the agent ID to the asset
472    for asset in best_assets.iter_mut() {
473        if let AssetState::Candidate = asset.state() {
474            asset
475                .make_mut()
476                .select_candidate_for_investment(agent.id.clone());
477        }
478    }
479
480    Ok(best_assets)
481}
482
483/// Check whether there is any remaining demand that is unmet in any time slice
484fn is_any_remaining_demand(demand: &DemandMap) -> bool {
485    demand.values().any(|flow| *flow > Flow(0.0))
486}
487
488/// Update capacity of chosen asset, if needed, and update both asset options and chosen assets
489fn update_assets(
490    mut best_asset: AssetRef,
491    capacity: Capacity,
492    opt_assets: &mut Vec<AssetRef>,
493    remaining_candidate_capacity: &mut HashMap<AssetRef, Capacity>,
494    best_assets: &mut Vec<AssetRef>,
495) {
496    match best_asset.state() {
497        AssetState::Commissioned { .. } => {
498            // Remove this asset from the options
499            opt_assets.retain(|asset| *asset != best_asset);
500            best_assets.push(best_asset);
501        }
502        AssetState::Candidate => {
503            // Remove this capacity from the available remaining capacity for this asset
504            let remaining_capacity = remaining_candidate_capacity.get_mut(&best_asset).unwrap();
505            *remaining_capacity -= capacity;
506
507            // If there's no capacity remaining, remove the asset from the options
508            if *remaining_capacity <= Capacity(0.0) {
509                let old_idx = opt_assets
510                    .iter()
511                    .position(|asset| *asset == best_asset)
512                    .unwrap();
513                opt_assets.swap_remove(old_idx);
514                remaining_candidate_capacity.remove(&best_asset);
515            }
516
517            if let Some(existing_asset) = best_assets.iter_mut().find(|asset| **asset == best_asset)
518            {
519                // If the asset is already in the list of best assets, add the additional required capacity
520                existing_asset.make_mut().increase_capacity(capacity);
521            } else {
522                // Otherwise, update the capacity of the chosen asset and add it to the list of best assets
523                best_asset.make_mut().set_capacity(capacity);
524                best_assets.push(best_asset);
525            };
526        }
527        _ => panic!("update_assets should only be called with Commissioned or Candidate assets"),
528    }
529}
530
531#[cfg(test)]
532mod tests {
533    use super::*;
534    use crate::commodity::Commodity;
535    use crate::fixture::{
536        asset, process, process_parameter_map, region_id, svd_commodity, time_slice,
537        time_slice_info, time_slice_info2,
538    };
539    use crate::process::{FlowType, ProcessFlow, ProcessParameter};
540    use crate::region::RegionID;
541    use crate::time_slice::{TimeSliceID, TimeSliceInfo};
542    use crate::units::{
543        ActivityPerCapacity, Dimensionless, Flow, FlowPerActivity, MoneyPerActivity,
544        MoneyPerCapacity, MoneyPerCapacityPerYear, MoneyPerFlow,
545    };
546    use indexmap::indexmap;
547    use itertools::Itertools;
548    use rstest::rstest;
549    use std::rc::Rc;
550
551    /// Custom fixture for process parameters with non-zero capacity_to_activity
552    fn process_parameter_with_capacity_to_activity() -> Rc<ProcessParameter> {
553        Rc::new(ProcessParameter {
554            capital_cost: MoneyPerCapacity(0.0),
555            fixed_operating_cost: MoneyPerCapacityPerYear(0.0),
556            variable_operating_cost: MoneyPerActivity(0.0),
557            lifetime: 1,
558            discount_rate: Dimensionless(1.0),
559            capacity_to_activity: ActivityPerCapacity(1.0), // Non-zero value
560        })
561    }
562
563    #[rstest]
564    fn test_get_demand_limiting_capacity(
565        time_slice: TimeSliceID,
566        region_id: RegionID,
567        time_slice_info: TimeSliceInfo,
568        svd_commodity: Commodity,
569    ) {
570        // Create a process flow using the existing commodity fixture
571        let commodity_rc = Rc::new(svd_commodity);
572        let process_flow = ProcessFlow {
573            commodity: Rc::clone(&commodity_rc),
574            coeff: FlowPerActivity(2.0), // 2 units of flow per unit of activity
575            kind: FlowType::Fixed,
576            cost: MoneyPerFlow(0.0),
577        };
578
579        // Create a process with the flows and activity limits
580        let mut process = process(
581            [region_id.clone()].into_iter().collect(),
582            process_parameter_map([region_id.clone()].into_iter().collect()),
583        );
584
585        // Add the flow to the process
586        process.flows.insert(
587            (region_id.clone(), 2015), // Using default commission year from fixture
588            [(commodity_rc.id.clone(), process_flow)]
589                .into_iter()
590                .collect(),
591        );
592
593        // Add activity limits
594        process.activity_limits.insert(
595            (region_id.clone(), 2015, time_slice.clone()),
596            Dimensionless(0.0)..=Dimensionless(1.0),
597        );
598
599        // Update process parameters to have non-zero capacity_to_activity
600        let updated_parameter = process_parameter_with_capacity_to_activity();
601        process
602            .parameters
603            .insert((region_id.clone(), 2015), updated_parameter);
604
605        // Create asset with the configured process
606        let asset = asset(process);
607
608        // Create demand map - demand of 10.0 for our time slice
609        let demand = indexmap! { time_slice.clone() => Flow(10.0)};
610
611        // Call the function
612        let result = get_demand_limiting_capacity(&time_slice_info, &asset, &commodity_rc, &demand);
613
614        // Expected calculation:
615        // max_flow_per_cap = activity_per_capacity_limit (1.0) * coeff (2.0) = 2.0
616        // required_capacity = demand (10.0) / max_flow_per_cap (2.0) = 5.0
617        assert_eq!(result, Capacity(5.0));
618    }
619
620    #[rstest]
621    fn test_get_demand_limiting_capacity_multiple_time_slices(
622        time_slice_info2: TimeSliceInfo,
623        svd_commodity: Commodity,
624        region_id: RegionID,
625    ) {
626        // Create time slices from the fixture (day and night)
627        let (time_slice1, time_slice2) =
628            time_slice_info2.time_slices.keys().collect_tuple().unwrap();
629
630        // Create a process flow using the existing commodity fixture
631        let commodity_rc = Rc::new(svd_commodity);
632        let process_flow = ProcessFlow {
633            commodity: Rc::clone(&commodity_rc),
634            coeff: FlowPerActivity(1.0), // 1 unit of flow per unit of activity
635            kind: FlowType::Fixed,
636            cost: MoneyPerFlow(0.0),
637        };
638
639        // Create a process with the flows and activity limits
640        let mut process = process(
641            [region_id.clone()].into_iter().collect(),
642            process_parameter_map([region_id.clone()].into_iter().collect()),
643        );
644
645        // Add the flow to the process
646        process.flows.insert(
647            (region_id.clone(), 2015), // Using default commission year from fixture
648            [(commodity_rc.id.clone(), process_flow)]
649                .into_iter()
650                .collect(),
651        );
652
653        // Add activity limits for both time slices with different limits
654        process.activity_limits.insert(
655            (region_id.clone(), 2015, time_slice1.clone()),
656            Dimensionless(0.0)..=Dimensionless(2.0), // Higher limit for day
657        );
658        process.activity_limits.insert(
659            (region_id.clone(), 2015, time_slice2.clone()),
660            Dimensionless(0.0)..=Dimensionless(0.0), // Zero limit for night - should be skipped
661        );
662
663        // Update process parameters to have non-zero capacity_to_activity
664        let updated_parameter = process_parameter_with_capacity_to_activity();
665        process
666            .parameters
667            .insert((region_id.clone(), 2015), updated_parameter);
668
669        // Create asset with the configured process
670        let asset = asset(process);
671
672        // Create demand map with different demands for each time slice
673        let demand = indexmap! {
674            time_slice1.clone() => Flow(4.0), // Requires capacity of 4.0/2.0 = 2.0
675            time_slice2.clone() => Flow(3.0), // Would require infinite capacity, but should be skipped
676        };
677
678        // Call the function
679        let result =
680            get_demand_limiting_capacity(&time_slice_info2, &asset, &commodity_rc, &demand);
681
682        // Expected: maximum of the capacity requirements across time slices (excluding zero limit)
683        // Time slice 1: demand (4.0) / (activity_limit (2.0) * coeff (1.0)) = 2.0
684        // Time slice 2: skipped due to zero activity limit
685        // Maximum = 2.0
686        assert_eq!(result, Capacity(2.0));
687    }
688}