Idea for simple implementation (not working...)

Author: tsmbland

src/model/parameters.rs

@@ -47,7 +47,7 @@ define_unit_param_default!(default_candidate_asset_capacity, Capacity, 0.0001);
 define_unit_param_default!(default_capacity_limit_factor, Dimensionless, 0.1);
 define_unit_param_default!(default_value_of_lost_load, MoneyPerFlow, 1e9);
 define_unit_param_default!(default_price_tolerance, Dimensionless, 1e-6);
-define_param_default!(default_max_ironing_out_iterations, u32, 10);
+define_param_default!(default_max_ironing_out_iterations, u32, 1);
 
 /// Model parameters as defined in the `model.toml` file.
 ///

src/simulation/investment.rs

@@ -9,11 +9,11 @@ use crate::region::RegionID;
 use crate::simulation::CommodityPrices;
 use crate::time_slice::{TimeSliceID, TimeSliceInfo};
 use crate::units::{Capacity, Dimensionless, Flow, FlowPerCapacity};
-use anyhow::{Result, bail, ensure};
+use anyhow::{Result, ensure};
 use indexmap::IndexMap;
 use itertools::{Itertools, chain, iproduct};
-use log::debug;
-use std::collections::HashMap;
+use log::{debug, warn};
+use std::collections::{HashMap, VecDeque};
 use std::fmt::Display;
 
 pub mod appraisal;
@@ -55,6 +55,8 @@ impl InvestmentSet {
         demand: &AllDemandMap,
         existing_assets: &[AssetRef],
         prices: &CommodityPrices,
+        seen_markets: &[MarketID],
+        previously_selected_assets: &[AssetRef],
         writer: &mut DataWriter,
     ) -> Result<Vec<AssetRef>> {
         match self {
@@ -67,15 +69,20 @@ impl InvestmentSet {
                 prices,
                 writer,
             ),
-            InvestmentSet::Cycle(markets) => select_assets_for_cycle(
-                model,
-                markets,
-                year,
-                demand,
-                existing_assets,
-                prices,
-                writer,
-            ),
+            InvestmentSet::Cycle(markets) => {
+                debug!("Starting investment for cycle '{self}'");
+                select_assets_for_cycle(
+                    model,
+                    markets,
+                    year,
+                    demand,
+                    existing_assets,
+                    prices,
+                    seen_markets,
+                    previously_selected_assets,
+                    writer,
+                )
+            }
             InvestmentSet::Layer(investment_sets) => {
                 debug!("Starting investment for layer '{self}'");
                 let mut all_assets = Vec::new();
@@ -86,6 +93,8 @@ impl InvestmentSet {
                         demand,
                         existing_assets,
                         prices,
+                        seen_markets,
+                        previously_selected_assets,
                         writer,
                     )?;
                     all_assets.extend(assets);
@@ -157,7 +166,9 @@ pub fn perform_agent_investment(
             year,
             &net_demand,
             existing_assets,
-            prices,
+            &external_prices,
+            &seen_markets,
+            &all_selected_assets,
             writer,
         )?;
 
@@ -277,16 +288,114 @@ fn select_assets_for_market(
     Ok(selected_assets)
 }
 
+#[allow(clippy::too_many_arguments)]
 fn select_assets_for_cycle(
-    _model: &Model,
+    model: &Model,
     markets: &[MarketID],
-    _year: u32,
-    _demand: &AllDemandMap,
-    _existing_assets: &[AssetRef],
-    _prices: &CommodityPrices,
-    _writer: &mut DataWriter,
+    year: u32,
+    demand: &AllDemandMap,
+    existing_assets: &[AssetRef],
+    prices: &CommodityPrices,
+    seen_markets: &[MarketID],
+    previously_selected_assets: &[AssetRef],
+    writer: &mut DataWriter,
 ) -> Result<Vec<AssetRef>> {
-    bail!("Investment cycles are not yet supported. Found cycle for commodities: {markets:?}");
+    // Get markets to balance: all seen so far plus all in loop
+    let mut markets_to_balance = seen_markets.to_vec();
+    markets_to_balance.extend_from_slice(markets);
+
+    // Initialise list of unbalanced markets: start with all markets in the cycle
+    let mut unbalanced_markets = VecDeque::from(markets.to_vec());
+
+    // Iterate while there are markets in the unbalanced list
+    let mut current_demand = demand.clone();
+    let mut selected_assets = HashMap::new();
+    let mut loop_iter = 0;
+    loop {
+        // If max iterations reached, break with a warning
+        if loop_iter == 10 {
+            warn!(
+                "Max cycle investment iterations reached. Unable to find stable solution for the following markets: {}",
+                unbalanced_markets.iter().join(", ")
+            );
+            break;
+        }
+
+        // Pop off the first market
+        // If there are no markets with unmet demand, we're done
+        let Some(current_market) = unbalanced_markets.pop_front() else {
+            debug!(
+                "Cycle investment for '{}' converged after {} iterations",
+                markets.iter().join(", "),
+                loop_iter
+            );
+            break;
+        };
+
+        // Select assets for this market
+        // This will replace any previously selected assets for this market
+        let assets = select_assets_for_market(
+            model,
+            &current_market,
+            year,
+            &current_demand,
+            existing_assets,
+            prices,
+            writer,
+        )?;
+        selected_assets.insert(current_market.clone(), assets);
+
+        // Assemble full list of assets for dispatch
+        let mut all_assets = previously_selected_assets.to_vec();
+        let flat_selected: Vec<AssetRef> = selected_assets
+            .values()
+            .flat_map(|v| v.iter().cloned())
+            .collect();
+        all_assets.extend_from_slice(&flat_selected);
+
+        // Run dispatch
+        // We allow unmet demand for all markets except the one we're currently balancing
+        let markets_with_unmet_demand: Vec<MarketID> = markets
+            .iter()
+            .filter(|m| *m != &current_market)
+            .cloned()
+            .collect();
+        let solution = DispatchRun::new(model, &all_assets, year)
+            .with_market_balance_subset(&markets_to_balance)
+            .with_unmet_demand_vars(&markets_with_unmet_demand)
+            .run(
+                &format!(
+                    "cycle ({}) iteration {}",
+                    markets.iter().join(", "),
+                    loop_iter
+                ),
+                writer,
+            )?;
+
+        // Find markets with unmet demand and add these to the unbalanced list if not already present
+        let coms = solution.get_markets_with_unmet_demand();
+        for market in coms {
+            if !unbalanced_markets.contains(&market) {
+                unbalanced_markets.push_back(market);
+            }
+        }
+        println!(
+            "Unbalanced markets: {}",
+            unbalanced_markets.iter().join(", ")
+        );
+
+        // Calculate a new demand map for the next iteration
+        current_demand.clone_from(demand);
+        update_demand_map_with_inputs(
+            &mut current_demand,
+            &solution.create_flow_map(),
+            &flat_selected,
+        );
+
+        loop_iter += 1;
+    }
+
+    Ok(selected_assets.into_values().flatten().collect())
 }
 
 /// Flatten the preset commodity demands for a given year into a map of commodity, region and
@@ -362,6 +471,28 @@ fn update_net_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: &[A
     }
 }
 
+/// Update demand map with input flows from a set of assets
+fn update_demand_map_with_inputs(demand: &mut AllDemandMap, flows: &FlowMap, assets: &[AssetRef]) {
+    for ((asset, commodity_id, time_slice), flow) in flows {
+        if assets.contains(asset) {
+            let key = (
+                commodity_id.clone(),
+                asset.region_id().clone(),
+                time_slice.clone(),
+            );
+
+            // Only consider input flows
+            if flow < &Flow(0.0) {
+                // Note: we use the negative of the flow as input flows are negative in the flow map.
+                demand
+                    .entry(key)
+                    .and_modify(|value| *value -= *flow)
+                    .or_insert(-*flow);
+            }
+        }
+    }
+}
+
 /// Get a portion of the demand profile for this commodity and region
 fn get_demand_portion_for_commodity(
     time_slice_info: &TimeSliceInfo,

src/simulation/optimisation.rs

@@ -280,6 +280,14 @@ impl Solution<'_> {
         keys.zip(output[variable_idx.clone()].iter())
             .map(|((asset, time_slice), value)| (asset, time_slice, T::new(*value)))
     }
+
+    /// Get the markets for which unment demand is positive in any time slice.
+    pub fn get_markets_with_unmet_demand(&self) -> IndexSet<MarketID> {
+        self.iter_unmet_demand()
+            .filter(|(_, _, flow)| *flow > Flow(0.0))
+            .map(|(market, _, _)| market.clone())
+            .collect()
+    }
 }
 
 /// Defines the possible errors that can occur when running the solver
@@ -461,19 +469,15 @@ impl<'model, 'run> DispatchRun<'model, 'run> {
     }
 
     /// Run dispatch to diagnose which markets have unmet demand
-    fn get_markets_with_unmet_demand(
+    pub fn get_markets_with_unmet_demand(
         &self,
         markets_to_balance: &[MarketID],
     ) -> Result<IndexSet<MarketID>> {
         // Run dispatch including unmet demand variables for all markets being balanced
         let solution = self.run_internal(markets_to_balance, markets_to_balance)?;
 
         // Collect markets with unmet demand
-        Ok(solution
-            .iter_unmet_demand()
-            .filter(|(_, _, flow)| *flow > Flow(0.0))
-            .map(|(market, _, _)| market.clone())
-            .collect())
+        Ok(solution.get_markets_with_unmet_demand())
     }
 
     /// Run dispatch to balance the specified markets, allowing unmet demand for a subset of these