Fix prices calculation

src/asset.rs

@@ -10,7 +10,7 @@ use crate::simulation::CommodityPrices;
 use crate::time_slice::{TimeSliceID, TimeSliceSelection};
 use crate::units::{
     Activity, ActivityPerCapacity, Capacity, Dimensionless, FlowPerActivity, MoneyPerActivity,
-    MoneyPerCapacity, MoneyPerFlow,
+    MoneyPerCapacity, MoneyPerFlow, Year,
 };
 use anyhow::{Context, Result, ensure};
 use indexmap::IndexMap;
@@ -756,33 +756,35 @@ impl Asset {
         annual_capital_cost(capital_cost, lifetime, discount_rate)
     }
 
-    /// Get the annual capital cost per unit of activity for this asset
+    /// Get the annual fixed costs (AFC) per unit of activity for this asset
     ///
-    /// Total capital costs (cost per capacity * capacity) are shared equally over the year in
-    /// accordance with the annual activity.
-    pub fn get_annual_capital_cost_per_activity(
+    /// Total capital costs and fixed opex are shared equally over the year in accordance with the
+    /// annual activity.
+    pub fn get_annual_fixed_costs_per_activity(
         &self,
         annual_activity: Activity,
     ) -> MoneyPerActivity {
         let annual_capital_cost_per_capacity = self.get_annual_capital_cost_per_capacity();
-        let total_annual_capital_cost = annual_capital_cost_per_capacity * self.total_capacity();
+        let annual_fixed_opex = self.process_parameter.fixed_operating_cost * Year(1.0);
+        let total_annual_fixed_costs =
+            (annual_capital_cost_per_capacity + annual_fixed_opex) * self.total_capacity();
         assert!(
             annual_activity > Activity::EPSILON,
-            "Cannot calculate annual capital cost per activity for an asset with zero annual activity"
+            "Cannot calculate annual fixed costs per activity for an asset with zero annual activity"
         );
-        total_annual_capital_cost / annual_activity
+        total_annual_fixed_costs / annual_activity
     }
 
-    /// Get the annual capital cost per unit of output flow for this asset
+    /// Get the annual fixed costs (AFC) per unit of output flow for this asset
     ///
-    /// Total capital costs (cost per capacity * capacity) are shared equally across all output
-    /// flows in accordance with the annual activity and total output per unit of activity.
-    pub fn get_annual_capital_cost_per_flow(&self, annual_activity: Activity) -> MoneyPerFlow {
-        let annual_capital_cost_per_activity =
-            self.get_annual_capital_cost_per_activity(annual_activity);
+    /// Total capital costs and fixed opex are shared equally across all output flows in accordance
+    /// with the annual activity and total output per unit of activity.
+    pub fn get_annual_fixed_costs_per_flow(&self, annual_activity: Activity) -> MoneyPerFlow {
+        let annual_fixed_costs_per_activity =
+            self.get_annual_fixed_costs_per_activity(annual_activity);
         let total_output_per_activity = self.get_total_output_per_activity();
         assert!(total_output_per_activity > FlowPerActivity::EPSILON); // input checks should guarantee this
-        annual_capital_cost_per_activity / total_output_per_activity
+        annual_fixed_costs_per_activity / total_output_per_activity
     }
 
     /// Maximum activity for this asset

src/graph/investment.rs

@@ -1,8 +1,9 @@
 //! Module for solving the investment order of commodities
 use super::{CommoditiesGraph, GraphEdge, GraphNode};
-use crate::commodity::{CommodityMap, CommodityType};
+use crate::commodity::{CommodityMap, CommodityType, PricingStrategy};
 use crate::region::RegionID;
 use crate::simulation::investment::InvestmentSet;
+use anyhow::{Result, ensure};
 use highs::{Col, HighsModelStatus, RowProblem, Sense};
 use indexmap::IndexMap;
 use log::warn;
@@ -41,22 +42,22 @@ fn solve_investment_order_for_year(
     graphs: &IndexMap<(RegionID, u32), CommoditiesGraph>,
     commodities: &CommodityMap,
     year: u32,
-) -> Vec<InvestmentSet> {
+) -> Result<Vec<InvestmentSet>> {
     // Initialise InvestmentGraph for this year from the set of original `CommodityGraph`s
     let mut investment_graph = init_investment_graph_for_year(graphs, year, commodities);
 
     // TODO: condense sibling commodities (commodities that share at least one producer)
 
     // Condense strongly connected components
-    investment_graph = compress_cycles(&investment_graph);
+    investment_graph = compress_cycles(&investment_graph, commodities)?;
 
     // Perform a topological sort on the condensed graph
     // We can safely unwrap because `toposort` will only return an error in case of cycles, which
     // should have been detected and compressed with `compress_cycles`
     let order = toposort(&investment_graph, None).unwrap();
 
     // Compute layers for investment
-    compute_layers(&investment_graph, &order)
+    Ok(compute_layers(&investment_graph, &order))
 }
 
 /// Initialise an `InvestmentGraph` for the given year from a set of `CommodityGraph`s
@@ -117,15 +118,39 @@ fn init_investment_graph_for_year(
 }
 
 /// Compresses cycles into `InvestmentSet::Cycle` nodes
-fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph {
+fn compress_cycles(graph: &InvestmentGraph, commodities: &CommodityMap) -> Result<InvestmentGraph> {
     // Detect strongly connected components
     let mut condensed_graph = condensation(graph.clone(), true);
 
     // Order nodes within each strongly connected component
     order_sccs(&mut condensed_graph, graph);
 
+    // Pre-scan SCCs for offending pricing strategies (FullCost / MarginalCost).
+    for node_weight in condensed_graph.node_weights() {
+        if node_weight.len() <= 1 {
+            continue;
+        }
+        let offenders: Vec<_> = node_weight
+            .iter()
+            .flat_map(|s| s.iter_markets())
+            .filter(|(cid, _)| {
+                matches!(
+                    &commodities[cid].pricing_strategy,
+                    PricingStrategy::MarginalCost | PricingStrategy::FullCost
+                )
+            })
+            .map(|(cid, _)| cid.clone())
+            .collect();
+
+        ensure!(
+            offenders.is_empty(),
+            "Cannot use FullCost/MarginalCost pricing strategies for commodities with circular \
+            dependencies. Offending commodities: {offenders:?}"
+        );
+    }
+
     // Map to a new InvestmentGraph
-    condensed_graph.map(
+    let mapped = condensed_graph.map(
         // Map nodes to InvestmentSet
         // If only one member, keep as-is; if multiple members, create Cycle
         |_, node_weight| match node_weight.len() {
@@ -141,7 +166,9 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph {
         },
         // Keep edges the same
         |_, edge_weight| edge_weight.clone(),
-    )
+    );
+
+    Ok(mapped)
 }
 
 /// Order the members of each strongly connected component using a mixed-integer linear program.
@@ -490,13 +517,13 @@ pub fn solve_investment_order_for_model(
     commodity_graphs: &IndexMap<(RegionID, u32), CommoditiesGraph>,
     commodities: &CommodityMap,
     years: &[u32],
-) -> HashMap<u32, Vec<InvestmentSet>> {
+) -> Result<HashMap<u32, Vec<InvestmentSet>>> {
     let mut investment_orders = HashMap::new();
     for year in years {
-        let order = solve_investment_order_for_year(commodity_graphs, commodities, *year);
+        let order = solve_investment_order_for_year(commodity_graphs, commodities, *year)?;
         investment_orders.insert(*year, order);
     }
-    investment_orders
+    Ok(investment_orders)
 }
 
 #[cfg(test)]
@@ -568,7 +595,7 @@ mod tests {
         commodities.insert("C".into(), Rc::new(svd_commodity));
 
         let graphs = IndexMap::from([(("GBR".into(), 2020), graph)]);
-        let result = solve_investment_order_for_year(&graphs, &commodities, 2020);
+        let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap();
 
         // Expected order: C, B, A (leaf nodes first)
         // No cycles or layers, so all investment sets should be `Single`
@@ -596,7 +623,7 @@ mod tests {
         commodities.insert("B".into(), Rc::new(sed_commodity));
 
         let graphs = IndexMap::from([(("GBR".into(), 2020), graph)]);
-        let result = solve_investment_order_for_year(&graphs, &commodities, 2020);
+        let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap();
 
         // Should be a single `Cycle` investment set containing both commodities
         assert_eq!(result.len(), 1);
@@ -635,7 +662,7 @@ mod tests {
         commodities.insert("D".into(), Rc::new(svd_commodity));
 
         let graphs = IndexMap::from([(("GBR".into(), 2020), graph)]);
-        let result = solve_investment_order_for_year(&graphs, &commodities, 2020);
+        let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap();
 
         // Expected order: D, Layer(B, C), A
         assert_eq!(result.len(), 3);
@@ -674,7 +701,7 @@ mod tests {
             (("GBR".into(), 2020), graph.clone()),
             (("FRA".into(), 2020), graph),
         ]);
-        let result = solve_investment_order_for_year(&graphs, &commodities, 2020);
+        let result = solve_investment_order_for_year(&graphs, &commodities, 2020).unwrap();
 
         // Expected order: Should have three layers, each with two commodities (one per region)
         assert_eq!(result.len(), 3);

src/input.rs

@@ -263,7 +263,8 @@ pub fn load_model<P: AsRef<Path>>(model_dir: P) -> Result<Model> {
     )?;
 
     // Solve investment order for each region/year
-    let investment_order = solve_investment_order_for_model(&commodity_graphs, &commodities, years);
+    let investment_order =
+        solve_investment_order_for_model(&commodity_graphs, &commodities, years)?;
 
     let model_path = model_dir
         .as_ref()

src/simulation/optimisation.rs

@@ -263,6 +263,14 @@ impl Solution<'_> {
             .map(|((asset, time_slice), &value)| (asset, time_slice, Activity(value)))
     }
 
+    /// Iterate over the keys for activity for each existing asset
+    pub fn iter_activity_keys_for_existing(
+        &self,
+    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID)> {
+        self.iter_activity_for_existing()
+            .map(|(asset, time_slice, _activity)| (asset, time_slice))
+    }
+
     /// Activity for each candidate asset
     pub fn iter_activity_for_candidates(
         &self,
@@ -275,6 +283,14 @@ impl Solution<'_> {
             .map(|((asset, time_slice), &value)| (asset, time_slice, Activity(value)))
     }
 
+    /// Iterate over the keys for activity for each candidate asset
+    pub fn iter_activity_keys_for_candidates(
+        &self,
+    ) -> impl Iterator<Item = (&AssetRef, &TimeSliceID)> {
+        self.iter_activity_for_candidates()
+            .map(|(asset, time_slice, _activity)| (asset, time_slice))
+    }
+
     /// Iterate over unmet demand
     pub fn iter_unmet_demand(
         &self,