diff --git a/Cargo.toml b/Cargo.toml index ee6ad13cf..6f7303465 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -78,7 +78,6 @@ missing_panics_doc = "allow" implicit_hasher = "allow" cast_lossless = "allow" return_self_not_must_use = "allow" -from_iter_instead_of_collect = "allow" [[bench]] name = "examples" diff --git a/benches/assets.rs b/benches/assets.rs index 9a2b0ba02..ba6414777 100644 --- a/benches/assets.rs +++ b/benches/assets.rs @@ -95,13 +95,15 @@ fn calculate_seed_prices( candidates: &[AssetRef], writer: &mut DataWriter, ) -> Prices { - let solution_existing = DispatchRun::new(model, base_year_assets, BASE_YEAR) + let market_demands = collect_preset_demands_for_year(&model.commodities, BASE_YEAR); + let solution_existing = DispatchRun::new(model, base_year_assets, BASE_YEAR, &market_demands) .run("bench setup: without candidates", writer) .expect("Dispatch without candidates failed"); - let solution_with_candidates = DispatchRun::new(model, base_year_assets, BASE_YEAR) - .with_candidates(candidates) - .run("bench setup: with candidates", writer) - .expect("Dispatch with candidates failed"); + let solution_with_candidates = + DispatchRun::new(model, base_year_assets, BASE_YEAR, &market_demands) + .with_candidates(candidates) + .run("bench setup: with candidates", writer) + .expect("Dispatch with candidates failed"); calculate_prices( model, diff --git a/examples/circularity/processes.csv b/examples/circularity/processes.csv index 1afc63624..314fbc7f4 100644 --- a/examples/circularity/processes.csv +++ b/examples/circularity/processes.csv @@ -1,19 +1,19 @@ -id,description,regions,primary_output,start_year,end_year,capacity_to_activity -GASDRV,Dry gas extraction,GBR,GASPRD,2020,,1.0 -OAGRSV,Oil and associated gas extraction,GBR,OILCRD,2020,,1.0 -GASPRC,Gas processing,GBR,GASNAT,2020,,1.0 -OILREF,Oil refinery,GBR,DIESEL,2020,,1.0 -OILRF2,O2G refinery,GBR,GASOLI,2020,,1.0 -WNDFRM,Wind farm,GBR,ELCTRI,2020,,31.54 -GASCGT,Gas combined cycle turbine,GBR,ELCTRI,2020,,31.54 -H2YGEN,H2 turbine,GBR,ELCTRI,2020,,31.54 -H2YPRO,H2 electrolyer,GBR,H2YPRD,2020,,1.0 -TPETCR,Petrol car,GBR,TPASKM,2020,,1.0 -TDIECR,Diesel car,GBR,TPASKM,2020,,1.0 -TELCCR,Electric car,GBR,TPASKM,2020,,1.0 -THYBCR,Plug-in hybrid car,GBR,TPASKM,2020,,1.0 -RGASBR,Gas boiler,GBR,RSHEAT,2020,,1.0 -RELCHP,Heat pump,GBR,RSHEAT,2020,,1.0 -BIOPRO,Biomass production,all,BIOPRD,2020,,1.0 -BIOPLL,Biomass pelletiser,all,BIOPEL,2020,,1.0 -RBIOBL,Biomass boiler,all,RSHEAT,2020,,1.0 +id,description,regions,primary_output,start_year,end_year,capacity_to_activity,feedback_process +GASDRV,Dry gas extraction,GBR,GASPRD,2020,,1.0, +OAGRSV,Oil and associated gas extraction,GBR,OILCRD,2020,,1.0, +GASPRC,Gas processing,GBR,GASNAT,2020,,1.0, +OILREF,Oil refinery,GBR,DIESEL,2020,,1.0, +OILRF2,O2G refinery,GBR,GASOLI,2020,,1.0, +WNDFRM,Wind farm,GBR,ELCTRI,2020,,31.54, +GASCGT,Gas combined cycle turbine,GBR,ELCTRI,2020,,31.54, +H2YGEN,H2 turbine,GBR,ELCTRI,2020,,31.54,true +H2YPRO,H2 electrolyer,GBR,H2YPRD,2020,,1.0,true +TPETCR,Petrol car,GBR,TPASKM,2020,,1.0, +TDIECR,Diesel car,GBR,TPASKM,2020,,1.0, +TELCCR,Electric car,GBR,TPASKM,2020,,1.0, +THYBCR,Plug-in hybrid car,GBR,TPASKM,2020,,1.0, +RGASBR,Gas boiler,GBR,RSHEAT,2020,,1.0, +RELCHP,Heat pump,GBR,RSHEAT,2020,,1.0, +BIOPRO,Biomass production,all,BIOPRD,2020,,1.0, +BIOPLL,Biomass pelletiser,all,BIOPEL,2020,,1.0, +RBIOBL,Biomass boiler,all,RSHEAT,2020,,1.0, diff --git a/schemas/input/model.yaml b/schemas/input/model.yaml index 5db20013a..fc8fb347c 100644 --- a/schemas/input/model.yaml +++ b/schemas/input/model.yaml @@ -95,10 +95,6 @@ properties: type: integer description: Number of iterations to perform when calculating prices for cyclically-dependent markets default: 1 - capacity_margin: - type: number - description: Slack proportion for assets selected during cycle balancing to absorb small demand shifts - default: 0.2 mothball_years: type: integer default: 0 diff --git a/src/fixture.rs b/src/fixture.rs index 6331a557c..1dfd4ed01 100644 --- a/src/fixture.rs +++ b/src/fixture.rs @@ -325,6 +325,7 @@ pub fn process( capacity_to_activity: ActivityPerCapacity(1.0), investment_constraints: process_investment_constraints, tranche_size: None, + feedback_process: false, } } diff --git a/src/graph.rs b/src/graph.rs index def0ddb06..d5d6389e8 100644 --- a/src/graph.rs +++ b/src/graph.rs @@ -14,6 +14,7 @@ use std::io::Write as IoWrite; use std::path::Path; use std::sync::Arc; +pub mod feedback_suggest; pub mod investment; pub mod validate; @@ -41,11 +42,21 @@ pub enum GraphNode { /// An edge in the commodity graph pub enum GraphEdge { /// An edge representing a primary flow of a process - #[display("{_0}")] - Primary(ProcessID), + #[display("{process_id}")] + Primary { + /// The process responsible for the flow. + process_id: ProcessID, + /// Whether the process participates in a feedback conversion. + feedback: bool, + }, /// An edge representing a secondary (non-primary) flow of a process - #[display("{_0}")] - Secondary(ProcessID), + #[display("{process_id}")] + Secondary { + /// The process responsible for the flow. + process_id: ProcessID, + /// Whether the process participates in a feedback conversion. + feedback: bool, + }, /// An edge representing a service demand #[display("DEMAND")] Demand, @@ -150,9 +161,15 @@ fn create_commodities_graph_for_region_year( source_node_index, target_node_index, if is_primary { - GraphEdge::Primary(process.id.clone()) + GraphEdge::Primary { + process_id: process.id.clone(), + feedback: process.feedback_process, + } } else { - GraphEdge::Secondary(process.id.clone()) + GraphEdge::Secondary { + process_id: process.id.clone(), + feedback: process.feedback_process, + } }, ); } @@ -181,7 +198,7 @@ pub fn build_commodity_graphs_for_model( fn get_edge_attributes(_: &CommoditiesGraph, edge_ref: EdgeReference) -> String { match edge_ref.weight() { // Use dashed lines for secondary flows - GraphEdge::Secondary(_) => "style=dashed".to_string(), + GraphEdge::Secondary { .. } => "style=dashed".to_string(), // Other edges use default attributes _ => String::new(), } diff --git a/src/graph/feedback_suggest.rs b/src/graph/feedback_suggest.rs new file mode 100644 index 000000000..871a8c073 --- /dev/null +++ b/src/graph/feedback_suggest.rs @@ -0,0 +1,560 @@ +//! Suggests the correct `feedback_process` configuration when a model's current settings leave the +//! commodity network unresolvable. +//! +//! When [`validate_non_feedback_commodity_graphs_for_model`](super::validate) fails, the network +//! cannot be resolved with the current feedback flags \u2014 either because a loop is not fully marked, +//! or because processes that are not part of any loop have been marked and removed needed edges. +//! This module searches for the smallest complete feedback sets that leave a network which is both +//! acyclic and structurally valid (which may be the empty set, meaning nothing should be marked). +//! +//! Validity is stricter than mere acyclicity: after the feedback edges are removed, a SED commodity +//! must be either fully connected (produced and consumed) or fully disconnected (neither). Breaking +//! a cycle by removing a single edge can leave a SED commodity half-connected (e.g. produced but no +//! longer consumed), which is invalid; resolving such a loop therefore requires marking every +//! process that keeps it half-connected. For the electricity/hydrogen loop this means both the +//! electrolyser and the turbine must be marked. +//! +//! # Algorithm +//! +//! 1. Enumerate every simple cycle in the commodity conversion graph with Johnson's algorithm and +//! reduce each to the *set of processes* on its edges. Cycles are detected over the full +//! topology, ignoring existing feedback flags, so a partially-marked loop is still found. +//! 2. Search for the smallest sets of processes to mark. A HiGHS MILP proposes a minimum-size set +//! that touches ("hits") every cycle; each proposal is verified by actually marking those +//! processes as feedback and re-running +//! `validate_non_feedback_commodity_graphs_for_region_year`. Reusing the real validator means +//! the suggestions can never drift from the validation rules. A proposal that fails validation +//! (for instance because it breaks a cycle but leaves a commodity half-connected) is excluded +//! with a no-good cut and the search continues. +//! +//! The search reports every minimal valid set (so the user sees genuine alternatives), or an empty +//! result if none exists — for example when the sole producer of an externally-consumed commodity +//! is the only way to break a cycle, which commonly arises with multiple-output processes. Each set +//! is the *complete* list of processes to mark, independent of what is already marked, so a process +//! that is currently marked but missing from a suggestion should be unmarked. +//! +//! Suggestions are generated for the single region/year graph whose validation failed. A suggestion +//! that fixes that graph may break another region/year, which the user resolves iteratively. +use super::validate::validate_non_feedback_commodity_graphs_for_region_year; +use super::{CommoditiesGraph, GraphEdge, GraphNode}; +use crate::commodity::CommodityMap; +use crate::process::{ProcessID, ProcessMap}; +use crate::region::RegionID; +use crate::time_slice::TimeSliceInfo; +use highs::{HighsModelStatus, RowProblem, Sense}; +use petgraph::visit::EdgeRef; +use std::collections::{BTreeSet, HashMap, HashSet}; + +/// The maximum number of alternative suggestions to report. +const MAX_ALTERNATIVES: usize = 5; + +/// A bound on solver iterations, to guard against pathological search spaces. +const MAX_ITERATIONS: usize = 1000; + +/// A candidate set of processes to mark as `feedback_process` to make the graph valid. +pub type FeedbackSuggestion = BTreeSet; + +/// Suggests the complete `feedback_process` sets that would make the given graph valid. +/// +/// `base_graph` is the failing region/year graph. Existing feedback flags on it are ignored when +/// evaluating candidates, so each suggestion is a complete feedback set. The remaining arguments +/// are what the validator needs to re-check a candidate. +/// +/// Returns every minimal valid set. This may be the empty set (the network has no loops, so nothing +/// should be marked), one or more non-empty sets (the processes that should be marked), or an empty +/// `Vec` if no valid configuration exists at all. +pub fn suggest_feedback_processes( + base_graph: &CommoditiesGraph, + processes: &ProcessMap, + commodities: &CommodityMap, + time_slice_info: &TimeSliceInfo, + region_id: &RegionID, + year: u32, +) -> Vec { + // Find every cycle, expressed as the set of processes that could break it. If there are none, + // the search still runs and reports the empty set (nothing should be marked) when that is valid. + let mut cycle_sets: HashSet> = HashSet::new(); + collect_cycle_sets(base_graph, &mut cycle_sets); + + // Search for minimal sets that break every cycle *and* pass full validation. + search_valid_sets(&cycle_sets, |candidate| { + candidate_passes_validation( + base_graph, + processes, + commodities, + time_slice_info, + region_id, + year, + candidate, + ) + }) +} + +/// Tests one candidate set as the *complete* feedback set, reusing the real validator. +/// +/// Existing feedback flags on `base_graph` are ignored: every edge is marked as feedback exactly +/// when its process is in the candidate, and cleared otherwise. A suggestion is therefore the full +/// set of processes that should be marked, so a process that is currently marked but absent from a +/// suggestion is one the user should unmark. Delegating to +/// `validate_non_feedback_commodity_graphs_for_region_year` keeps this check identical to real +/// validation. +fn candidate_passes_validation( + base_graph: &CommoditiesGraph, + processes: &ProcessMap, + commodities: &CommodityMap, + time_slice_info: &TimeSliceInfo, + region_id: &RegionID, + year: u32, + candidate: &BTreeSet, +) -> bool { + let mut graph = base_graph.clone(); + for weight in graph.edge_weights_mut() { + // Demand edges carry no process, so they can never be feedback edges; skip them. + let (GraphEdge::Primary { + process_id, + feedback, + } + | GraphEdge::Secondary { + process_id, + feedback, + }) = weight + else { + continue; + }; + // Overwrite existing marks: the candidate is the complete feedback set under test. + *feedback = candidate.contains(process_id); + } + validate_non_feedback_commodity_graphs_for_region_year( + &graph, + region_id, + year, + processes, + commodities, + time_slice_info, + ) + .is_ok() +} + +/// Builds the top-level guidance for a `feedback_process` configuration failure. +/// +/// The wording distinguishes the three outcomes of [`suggest_feedback_processes`]: mark a specific +/// set (or one of several alternatives), mark nothing at all (an acyclic network), or no valid +/// configuration exists. +pub fn feedback_error_message(suggestions: &[FeedbackSuggestion]) -> String { + let base = "The commodity network cannot be resolved with the current `feedback_process` \ + settings."; + match suggestions { + // A cyclic network with no valid marking (e.g. a sole producer must break the loop). + [] => format!("{base} No valid `feedback_process` configuration could be found."), + // An acyclic network: the only valid set is empty, so nothing should be marked. + [only] if only.is_empty() => { + format!("{base} No processes in this model should be marked as `feedback_process`.") + } + _ => format!( + "{base} Recommended to set `feedback_process = true` for: {}.", + format_feedback_suggestions(suggestions) + ), + } +} + +/// Formats the non-empty suggestion sets for the error message. +/// +/// A single option is rendered as a bare list (`A, B`); multiple options are each wrapped in +/// parentheses and joined with `" or "` (`(A, B) or (C)`) so the alternatives read unambiguously. +fn format_feedback_suggestions(suggestions: &[FeedbackSuggestion]) -> String { + let multiple = suggestions.len() > 1; + suggestions + .iter() + .map(|processes| { + let ids = processes + .iter() + .map(ToString::to_string) + .collect::>() + .join(", "); + if multiple { format!("({ids})") } else { ids } + }) + .collect::>() + .join(" or ") +} + +/// Collects, for a single graph, the set of processes on every simple cycle in the non-feedback +/// network. A valid feedback set must remove at least one process from each. +fn collect_cycle_sets(graph: &CommoditiesGraph, cycle_sets: &mut HashSet>) { + // Processes on each commodity-to-commodity, non-feedback edge. + let mut edge_procs: HashMap<(usize, usize), BTreeSet> = HashMap::new(); + for edge in graph.edge_references() { + let Some((process_id, _)) = edge_label(edge.weight()) else { + continue; + }; + let is_commodity = |idx| matches!(graph.node_weight(idx), Some(GraphNode::Commodity(_))); + if is_commodity(edge.source()) && is_commodity(edge.target()) { + edge_procs + .entry((edge.source().index(), edge.target().index())) + .or_default() + .insert(process_id.clone()); + } + } + + let mut adj = vec![Vec::new(); graph.node_count()]; + for &(a, b) in edge_procs.keys() { + adj[a].push(b); + } + // Johnson requires deterministic successor order for reproducible output. + for row in &mut adj { + row.sort_unstable(); + } + + // Reduce each cycle to the union of processes on its edges. Removing any one of them breaks the + // cycle, so the cover constraint over this set is "mark at least one". + for cycle in Johnson::new(adj).elementary_circuits() { + let procs = (0..cycle.len()) + .flat_map(|i| &edge_procs[&(cycle[i], cycle[(i + 1) % cycle.len()])]) + .cloned() + .collect(); + cycle_sets.insert(procs); + } +} + +/// Searches for the minimal sets of processes that break every cycle and pass `is_valid`. +/// +/// This is a minimum hitting-set problem solved by repeated MILP calls ("verify and cut"): +/// [`solve_once`] returns a smallest process set that hits every cycle, we verify it with +/// `is_valid`, and — valid or not — forbid that exact assignment before solving again. Because the +/// MILP always returns a minimum-size solution, results come out in non-decreasing size, so once we +/// have a valid set of size `n` we can stop as soon as the solver is forced to a larger size. This +/// yields all minimal valid sets (the alternatives) without exploring larger ones. +/// +/// Indices are used throughout because the MILP works over integer columns: `candidates[i]` is the +/// process for column `i`, and `cover_constraints`/`forbidden` are expressed in those indices. +fn search_valid_sets( + cycle_sets: &HashSet>, + is_valid: impl Fn(&BTreeSet) -> bool, +) -> Vec { + // The processes on any cycle are the only ones worth marking; number them for the MILP. + let candidates: Vec = cycle_sets + .iter() + .flatten() + .cloned() + .collect::>() + .into_iter() + .collect(); + + // With no cycles the only possible feedback set is empty: nothing should be marked. Report it + // if valid so the caller can tell the user to unmark everything. + if candidates.is_empty() { + return is_valid(&BTreeSet::new()) + .then(BTreeSet::new) + .into_iter() + .collect(); + } + + let index: HashMap<&ProcessID, usize> = + candidates.iter().enumerate().map(|(i, p)| (p, i)).collect(); + + // One cover constraint per cycle, as candidate indices: at least one of them must be marked. + let cover_constraints: Vec> = cycle_sets + .iter() + .map(|set| set.iter().map(|p| index[p]).collect()) + .collect(); + + let mut valid: Vec> = Vec::new(); + let mut forbidden: Vec> = Vec::new(); + let mut target_size: Option = None; + + // MAX_ITERATIONS bounds validator calls; MAX_ALTERNATIVES bounds how many options we report. + for _ in 0..MAX_ITERATIONS { + if valid.len() >= MAX_ALTERNATIVES { + break; + } + let Some(marked) = solve_once(candidates.len(), &cover_constraints, &forbidden) else { + break; // No further cycle-breaking set exists. + }; + // Solutions only grow in size, so once we exceed the first valid size we have them all. + if target_size.is_some_and(|size| marked.len() > size) { + break; + } + // Forbid this exact set next time: if valid we want a *different* alternative, if invalid we + // must not see it again. + forbidden.push(marked.clone()); + + let candidate: BTreeSet = + marked.iter().map(|&i| candidates[i].clone()).collect(); + if is_valid(&candidate) { + target_size.get_or_insert(marked.len()); + valid.push(candidate); + } + } + + valid +} + +/// Solves one instance of the hitting-set MILP: pick the fewest candidate processes to mark. +/// +/// Each candidate is a binary column (1 = mark as feedback) with objective coefficient 1, so +/// minimising the objective minimises the number of processes marked. Two families of rows +/// constrain the solution: +/// +/// * **Cover** — for every cycle, `sum(x_i) >= 1`, forcing at least one of its processes to be +/// marked. +/// * **No-good cuts** — each previously returned assignment `S` is excluded *exactly* (supersets +/// stay reachable, which matters when a small set breaks a cycle but fails validation and a larger +/// one is required). Requiring at least one variable to differ from `S` gives +/// `sum_{i in S}(1 - x_i) + sum_{i not in S}(x_i) >= 1`, which rearranges to a single linear row +/// with coefficient -1 for columns in `S`, +1 for the rest, and lower bound `1 - |S|`. +/// +/// Returns the marked candidate indices, or `None` if no assignment satisfies the constraints. +fn solve_once( + num_candidates: usize, + cover_constraints: &[Vec], + forbidden: &[Vec], +) -> Option> { + let mut problem = RowProblem::default(); + // One binary column per candidate; cost 1 each so the objective counts marked processes. + let cols: Vec<_> = (0..num_candidates) + .map(|_| problem.add_integer_column(1.0, 0..=1)) + .collect(); + + // Cover: every cycle must lose at least one process. + for cycle in cover_constraints { + let terms: Vec<_> = cycle.iter().map(|&i| (cols[i], 1.0)).collect(); + problem.add_row(1.0.., terms); + } + + // No-good cut per forbidden assignment: -1 for its members, +1 for the rest, bound 1 - |S|. + for assignment in forbidden { + let mut in_assignment = vec![false; num_candidates]; + for &i in assignment { + in_assignment[i] = true; + } + let terms: Vec<_> = (0..num_candidates) + .map(|i| (cols[i], if in_assignment[i] { -1.0 } else { 1.0 })) + .collect(); + problem.add_row((1.0 - count_as_f64(assignment.len())).., terms); + } + + let solved = problem.optimise(Sense::Minimise).try_solve().ok()?; + if solved.status() != HighsModelStatus::Optimal { + return None; + } + let solution = solved.get_solution(); + // Columns above 0.5 are the marked processes (values are binary, so effectively == 1). + Some( + (0..num_candidates) + .filter(|&i| solution[cols[i]] > 0.5) + .collect(), + ) +} + +/// Converts a small count to an `f64` without triggering precision-loss lints. +fn count_as_f64(count: usize) -> f64 { + f64::from(u32::try_from(count).expect("count fits in u32")) +} + +/// Returns the process and feedback flag for a conversion edge, or `None` for demand edges. +fn edge_label(edge: &GraphEdge) -> Option<(&ProcessID, bool)> { + match edge { + GraphEdge::Primary { + process_id, + feedback, + } + | GraphEdge::Secondary { + process_id, + feedback, + } => Some((process_id, *feedback)), + GraphEdge::Demand => None, + } +} + +/// Johnson's algorithm for enumerating the elementary circuits (simple cycles) of a directed graph. +/// +/// Reference: Donald B. Johnson, "Finding all the elementary circuits of a directed graph" (1975). +/// Each circuit is reported exactly once, canonicalised by its least-indexed vertex: the search is +/// restricted to the subgraph on vertices `>= s` and rooted at `s`, so a circuit is only emitted +/// when `s` is its smallest vertex. The `blocked`/`b` bookkeeping is what keeps the algorithm from +/// being exponential in dead-end exploration — a vertex stays blocked until unblocking it could +/// actually yield a new circuit. +struct Johnson { + /// Successor lists, indexed by vertex (a commodity node's index). + adj: Vec>, + /// Whether each vertex is currently blocked from being revisited on the active search. + blocked: Vec, + /// Deferred-unblock lists: `b[w]` holds the vertices to unblock once `w` is unblocked. + b: Vec>, + /// The vertices on the path currently being explored. + stack: Vec, + /// The circuits found so far, each as a list of vertices. + circuits: Vec>, +} + +impl Johnson { + fn new(adj: Vec>) -> Self { + let n = adj.len(); + Self { + adj, + blocked: vec![false; n], + b: vec![BTreeSet::new(); n], + stack: Vec::new(), + circuits: Vec::new(), + } + } + + /// Runs the search once from each start vertex `s` and returns every circuit found. + fn elementary_circuits(mut self) -> Vec> { + for s in 0..self.adj.len() { + // Reset per-root state; only vertices in the subgraph induced by `>= s` are visited. + for v in s..self.adj.len() { + self.blocked[v] = false; + self.b[v].clear(); + } + self.stack.clear(); + self.circuit(s, s); + } + self.circuits + } + + /// Depth-first search for circuits rooted at `s`, extending the current path with `v`. + /// + /// Returns whether a circuit back to `s` was found through `v`. This drives whether `v` is + /// unblocked immediately (a circuit closed, so `v` may yield more) or its unblocking is deferred + /// via `b` (retried only once one of its successors is unblocked). + fn circuit(&mut self, v: usize, s: usize) -> bool { + let mut found = false; + self.stack.push(v); + self.blocked[v] = true; + + // Only consider the subgraph induced by vertices with index >= s. + for w in self.adj[v].clone() { + if w < s { + continue; + } + if w == s { + // Reaching the root closes an elementary circuit: the current path. + self.circuits.push(self.stack.clone()); + found = true; + } else if !self.blocked[w] && self.circuit(w, s) { + found = true; + } + } + + if found { + self.unblock(v); + } else { + // No circuit through `v` yet; defer its unblocking to when a successor is unblocked. + for w in self.adj[v].clone() { + if w >= s { + self.b[w].insert(v); + } + } + } + + self.stack.pop(); + found + } + + /// Marks `u` explorable again and recursively unblocks the vertices that were waiting on it. + fn unblock(&mut self, u: usize) { + self.blocked[u] = false; + for w in std::mem::take(&mut self.b[u]) { + if self.blocked[w] { + self.unblock(w); + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use petgraph::graph::Graph; + + /// Builds a non-feedback primary edge for the given process. + fn primary(process_id: &str) -> GraphEdge { + GraphEdge::Primary { + process_id: process_id.into(), + feedback: false, + } + } + + #[test] + fn collect_cycle_sets_finds_loop_processes() { + // A<->B loop, plus a Source edge that must be ignored (not part of any cycle). + let mut graph: CommoditiesGraph = Graph::new(); + let src = graph.add_node(GraphNode::Source); + let a = graph.add_node(GraphNode::Commodity("A".into())); + let b = graph.add_node(GraphNode::Commodity("B".into())); + graph.add_edge(src, a, primary("ext")); + graph.add_edge(a, b, primary("p_ab")); + graph.add_edge(b, a, primary("p_ba")); + + let mut cycle_sets = HashSet::new(); + collect_cycle_sets(&graph, &mut cycle_sets); + + assert_eq!( + cycle_sets, + HashSet::from([BTreeSet::from(["p_ab".into(), "p_ba".into()])]) + ); + } + + #[test] + fn search_reports_all_minimal_alternatives() { + // Single 3-cycle where any one process is a valid fix. + let cycle_sets = HashSet::from([BTreeSet::from(["p1".into(), "p2".into(), "p3".into()])]); + + let options: HashSet> = + search_valid_sets(&cycle_sets, |c| c.len() == 1) + .into_iter() + .collect(); + + assert_eq!( + options, + HashSet::from([ + BTreeSet::from(["p1".into()]), + BTreeSet::from(["p2".into()]), + BTreeSet::from(["p3".into()]), + ]) + ); + } + + #[test] + fn search_finds_superset_when_singletons_invalid() { + // Only the whole loop is valid, so singletons are cut and the pair is found. + let cycle_sets = HashSet::from([BTreeSet::from(["p_ab".into(), "p_ba".into()])]); + let both: BTreeSet = BTreeSet::from(["p_ab".into(), "p_ba".into()]); + + let suggestions = search_valid_sets(&cycle_sets, |c| *c == both); + + assert_eq!(suggestions.len(), 1); + assert_eq!(suggestions[0], both); + } + + #[test] + fn search_returns_empty_when_nothing_valid() { + let cycle_sets = HashSet::from([BTreeSet::from(["p1".into(), "p2".into()])]); + assert!(search_valid_sets(&cycle_sets, |_| false).is_empty()); + } + + #[test] + fn search_reports_empty_set_when_acyclic() { + // No cycles: the only candidate feedback set is empty, and it is valid. + let suggestions = search_valid_sets(&HashSet::new(), |_| true); + assert_eq!(suggestions, vec![BTreeSet::new()]); + } + + #[test] + fn message_recommends_unmarking_for_empty_set() { + let message = feedback_error_message(&[BTreeSet::new()]); + assert!(message.contains("No processes in this model should be marked")); + } + + #[test] + fn message_lists_processes_to_mark() { + let message = feedback_error_message(&[BTreeSet::from(["A".into(), "B".into()])]); + assert!(message.contains("Set `feedback_process = true` for exactly: A, B")); + } + + #[test] + fn message_reports_no_valid_configuration() { + let message = feedback_error_message(&[]); + assert!(message.contains("No valid `feedback_process` configuration")); + } +} diff --git a/src/graph/investment.rs b/src/graph/investment.rs index 5f30ac1ec..61096083c 100644 --- a/src/graph/investment.rs +++ b/src/graph/investment.rs @@ -162,10 +162,12 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { /// `i` comes before `j`, then `j` cannot be before `i`). /// * Transitivity constraints prevent 3-cycles, ensuring the resulting relation is acyclic (i.e. if /// `i` comes before `j` and `j` comes before `k`, then `k` cannot come before `i`). -/// * The objective minimises the number of “forward” edges (edges that would point from an earlier -/// market to a later one), counted within the original SCC and treated as unit penalties. A small -/// bias (<1) is added to nudge exporters earlier without outweighing the main objective (a bias -/// >1 would instead prioritise exporters even if it created extra conflicts in the final order). +/// * The objective minimises unit penalties for violating conversion ordering. For a normal edge +/// `i -> j`, the penalty favours `j` before `i`; for a feedback edge, it favours `i` before `j`. +/// Feedback labels are checked for consistency across each conversion before this function is +/// called. A small bias (<1) is added to nudge exporters earlier without outweighing the main +/// objective (a bias >1 would instead prioritise exporters even if it created extra conflicts in +/// the final order). /// /// Once the MILP is solved, markets are scored by the number of pairwise “wins” (how many other /// markets they precede). Sorting by this score — using the original index as a tiebreaker to keep @@ -183,8 +185,9 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { /// /// Additionally, C has an outgoing edge to a node outside the cycle. /// -/// The costs matrix in the MILP is set up to penalise any edge that points “forward” in the final -/// order: if there's an edge from X to Y we prefer to place Y before X so the edge points backwards: +/// The costs matrix in the MILP is set up to penalise any non-feedback edge that points “forward” +/// in the final order: if there's an edge from X to Y we prefer to place Y before X so the edge +/// points backwards. Feedback edges reverse this preference, placing X before Y: /// /// ```text /// | | A | B | C | @@ -193,6 +196,10 @@ fn compress_cycles(graph: &InvestmentGraph) -> InvestmentGraph { /// | C | 0 | 1 | 0 | /// ``` /// +/// For feedback edges, the corresponding penalty is placed in the opposite matrix entry. If all +/// edges in a cycle are feedback edges, the directional penalties can therefore balance, leaving +/// the external-outgoing bias to determine the order. +/// /// On top of this, we give a small preference to markets that export outside the SCC, so nodes with /// outgoing edges beyond the cycle are pushed earlier. This is done via an `EXTERNAL_BIAS` /// parameter (B) applied to the cost matrix: @@ -276,9 +283,23 @@ fn order_sccs( for (i, &idx) in original_indices.iter().enumerate() { // Loop over the edges going out of this node for edge in original_graph.edges_directed(idx, Direction::Outgoing) { - // If the target j is inside this SCC, record a penalty for putting i before j + // If the target j is inside this SCC, record a penalty for the preferred order. if let Some(&j) = index_position.get(&edge.target()) { - penalties[i][j] = 1.0; + let feedback = match edge.weight() { + GraphEdge::Primary { feedback, .. } + | GraphEdge::Secondary { feedback, .. } => *feedback, + GraphEdge::Demand => unreachable!( + "Demand edges should not be present in the investment graph" + ), + }; + if feedback { + // Feedback conversions prefer the source before the target, so penalise + // the opposite ordering variable from a normal conversion. + penalties[j][i] = 1.0; + } else { + // Normal conversions retain the downstream-first preference. + penalties[i][j] = 1.0; + } // Otherwise, mark that i has an outgoing edge to outside the SCC } else { @@ -520,7 +541,10 @@ mod tests { original.add_edge( node_indices[src], node_indices[dst], - GraphEdge::Primary("process1".into()), + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, ); } // External market receiving exports from C; encourages C to appear early. @@ -528,7 +552,10 @@ mod tests { original.add_edge( node_indices[2], external, - GraphEdge::Primary("process2".into()), + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, ); // Single SCC containing all markets. @@ -556,8 +583,22 @@ mod tests { let node_c = graph.add_node(GraphNode::Commodity("C".into())); // Add edges: A -> B -> C - graph.add_edge(node_a, node_b, GraphEdge::Primary("process1".into())); - graph.add_edge(node_b, node_c, GraphEdge::Primary("process2".into())); + graph.add_edge( + node_a, + node_b, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); + graph.add_edge( + node_b, + node_c, + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, + ); // Create commodities map using fixtures let mut commodities = CommodityMap::new(); @@ -585,8 +626,22 @@ mod tests { let node_b = graph.add_node(GraphNode::Commodity("B".into())); // Add edges creating a cycle: A -> B -> A - graph.add_edge(node_a, node_b, GraphEdge::Primary("process1".into())); - graph.add_edge(node_b, node_a, GraphEdge::Primary("process2".into())); + graph.add_edge( + node_a, + node_b, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); + graph.add_edge( + node_b, + node_a, + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, + ); // Create commodities map using fixtures let mut commodities = CommodityMap::new(); @@ -620,10 +675,38 @@ mod tests { let node_d = graph.add_node(GraphNode::Commodity("D".into())); // Add edges - graph.add_edge(node_a, node_b, GraphEdge::Primary("process1".into())); - graph.add_edge(node_a, node_c, GraphEdge::Primary("process2".into())); - graph.add_edge(node_b, node_d, GraphEdge::Primary("process3".into())); - graph.add_edge(node_c, node_d, GraphEdge::Primary("process4".into())); + graph.add_edge( + node_a, + node_b, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); + graph.add_edge( + node_a, + node_c, + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, + ); + graph.add_edge( + node_b, + node_d, + GraphEdge::Primary { + process_id: "process3".into(), + feedback: false, + }, + ); + graph.add_edge( + node_c, + node_d, + GraphEdge::Primary { + process_id: "process4".into(), + feedback: false, + }, + ); // Create commodities map using fixtures let mut commodities = CommodityMap::new(); @@ -658,8 +741,22 @@ mod tests { let node_c = graph.add_node(GraphNode::Commodity("C".into())); // Add edges: A -> B -> C - graph.add_edge(node_a, node_b, GraphEdge::Primary("process1".into())); - graph.add_edge(node_b, node_c, GraphEdge::Primary("process2".into())); + graph.add_edge( + node_a, + node_b, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); + graph.add_edge( + node_b, + node_c, + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, + ); // Create commodities map using fixtures let mut commodities = CommodityMap::new(); diff --git a/src/graph/validate.rs b/src/graph/validate.rs index b750ae9db..0864580fd 100644 --- a/src/graph/validate.rs +++ b/src/graph/validate.rs @@ -1,14 +1,89 @@ //! Module for validating commodity graphs +use super::feedback_suggest::{feedback_error_message, suggest_feedback_processes}; use super::{CommoditiesGraph, GraphEdge, GraphNode}; -use crate::commodity::{CommodityMap, CommodityType}; -use crate::process::{Process, ProcessMap}; +use crate::commodity::{CommodityID, CommodityMap, CommodityType}; +use crate::process::{Process, ProcessID, ProcessMap}; use crate::region::RegionID; use crate::time_slice::{TimeSliceInfo, TimeSliceLevel, TimeSliceSelection}; use crate::units::{Dimensionless, Flow}; -use anyhow::{Context, Result, ensure}; +use anyhow::{Context, Result, bail, ensure}; use indexmap::IndexMap; +use petgraph::algo::toposort; +use petgraph::visit::EdgeRef; +use std::collections::HashMap; use strum::IntoEnumIterator; +/// Checks that all processes representing a SED-to-SED conversion agree on whether they are +/// feedback processes. +fn validate_feedback_processes(graph: &CommoditiesGraph, commodities: &CommodityMap) -> Result<()> { + // Parallel process edges represent the same commodity conversion. Grouping them lets us + // compare every process contributing to that conversion rather than choosing one edge. + let mut conversions: HashMap<(CommodityID, CommodityID), Vec<(ProcessID, bool)>> = + HashMap::new(); + + for edge in graph.edge_references() { + // Source, sink, and demand edges do not describe a commodity-to-commodity conversion. + let (Some(GraphNode::Commodity(source)), Some(GraphNode::Commodity(target))) = ( + graph.node_weight(edge.source()), + graph.node_weight(edge.target()), + ) else { + continue; + }; + + // Only SED markets can be both consumed and produced, and therefore only SED-to-SED + // conversions can participate in an SCC and influence the investment ordering. + if commodities[source].kind != CommodityType::SupplyEqualsDemand + || commodities[target].kind != CommodityType::SupplyEqualsDemand + { + continue; + } + + let process = match edge.weight() { + GraphEdge::Primary { + process_id, + feedback, + } + | GraphEdge::Secondary { + process_id, + feedback, + } => (process_id.clone(), *feedback), + GraphEdge::Demand => continue, + }; + + conversions + .entry((source.clone(), target.clone())) + .or_default() + .push(process); + } + + // Check each conversion independently. A conversion is valid when every contributing process + // agrees, whether that shared label is feedback or non-feedback. + for ((source, target), processes) in conversions { + let Some((_, first_feedback)) = processes.first() else { + continue; + }; + if processes + .iter() + .all(|(_, feedback)| feedback == first_feedback) + { + continue; + } + + // Include every contributing process in the error so the conflicting input rows are easy + // to locate and correct. + let process_details = processes + .iter() + .map(|(process_id, feedback)| format!("{process_id}={feedback}")) + .collect::>() + .join(", "); + bail!( + "Conversion {source} -> {target} has conflicting feedback_process values: {process_details}" + ); + } + + Ok(()) +} + /// Prepares a graph for validation with [`validate_commodities_graph`]. /// /// It takes a base graph produced by `create_commodities_graph_for_region_year`, and modifies it to @@ -35,7 +110,9 @@ fn prepare_commodities_graph_for_validation( filtered_graph.retain_edges(|graph, edge_idx| { // Get the process for the edge let process_id = match graph.edge_weight(edge_idx).unwrap() { - GraphEdge::Primary(process_id) | GraphEdge::Secondary(process_id) => process_id, + GraphEdge::Primary { process_id, .. } | GraphEdge::Secondary { process_id, .. } => { + process_id + } GraphEdge::Demand => panic!("Demand edges should not be present in the base graph"), }; let process = &processes[process_id]; @@ -208,6 +285,7 @@ pub fn validate_commodity_graphs_for_model( ) -> Result<()> { // Validate graphs at all time slice levels (taking into account process availability and demand) for ((region_id, year), base_graph) in commodity_graphs { + // Validate the original graph for ts_level in TimeSliceLevel::iter() { for ts_selection in time_slice_info.iter_selections_at_level(ts_level) { let graph = prepare_commodities_graph_for_validation( @@ -230,6 +308,88 @@ pub fn validate_commodity_graphs_for_model( Ok(()) } +/// Validate graphs without feedback processes +pub fn validate_non_feedback_commodity_graphs_for_model( + commodity_graphs: &IndexMap<(RegionID, u32), CommoditiesGraph>, + processes: &ProcessMap, + commodities: &CommodityMap, + time_slice_info: &TimeSliceInfo, +) -> Result<()> { + // Validate graphs at all time slice levels (taking into account process availability and demand) + for ((region_id, year), base_graph) in commodity_graphs { + validate_non_feedback_commodity_graphs_for_region_year( + base_graph, + region_id, + *year, + processes, + commodities, + time_slice_info, + ) + .with_context(|| format!("Error for {region_id} in {year}.")) + .with_context(|| { + let suggestions = suggest_feedback_processes( + base_graph, + processes, + commodities, + time_slice_info, + region_id, + *year, + ); + feedback_error_message(&suggestions) + })?; + } + Ok(()) +} + +pub(super) fn validate_non_feedback_commodity_graphs_for_region_year( + base_graph: &CommoditiesGraph, + region_id: &RegionID, + year: u32, + processes: &ProcessMap, + commodities: &CommodityMap, + time_slice_info: &TimeSliceInfo, +) -> Result<()> { + // Check for consistency amongst feedback processes + validate_feedback_processes(base_graph, commodities)?; + + // Create a filtered graph without feedback processes + let mut graph_without_feedback = base_graph.clone(); + graph_without_feedback.retain_edges(|graph, edge_idx| match graph.edge_weight(edge_idx) { + Some(GraphEdge::Primary { feedback, .. } | GraphEdge::Secondary { feedback, .. }) => { + !feedback + } + Some(GraphEdge::Demand) => true, + None => unreachable!("Retained graph edge must have a weight"), + }); + + // Check that this graph is acyclic + if let Err(cycle) = toposort(&graph_without_feedback, None) { + let cycle_node = graph_without_feedback + .node_weight(cycle.node_id()) + .expect("cycle node must exist"); + bail!("Graph contains a cycle involving {cycle_node}"); + } + + // Validate the network that remains after feedback processes are removed. + // Do not need to quote the time slice selection in the error message, as `feedback_process` + // is a whole year property, so any breakages here are irrelevant of the ts selection + for ts_level in TimeSliceLevel::iter() { + for ts_selection in time_slice_info.iter_selections_at_level(ts_level) { + let graph = prepare_commodities_graph_for_validation( + &graph_without_feedback, + processes, + commodities, + region_id, + year, + &ts_selection, + ); + validate_commodities_graph(&graph, commodities, ts_level)?; + } + } + + Ok(()) +} + #[cfg(test)] mod tests { use super::*; @@ -258,8 +418,22 @@ mod tests { let node_b = graph.add_node(GraphNode::Commodity("B".into())); let node_c = graph.add_node(GraphNode::Commodity("C".into())); let node_d = graph.add_node(GraphNode::Demand); - graph.add_edge(node_a, node_b, GraphEdge::Primary("process1".into())); - graph.add_edge(node_b, node_c, GraphEdge::Primary("process2".into())); + graph.add_edge( + node_a, + node_b, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); + graph.add_edge( + node_b, + node_c, + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, + ); graph.add_edge(node_c, node_d, GraphEdge::Demand); // Validate the graph at DayNight level @@ -284,8 +458,22 @@ mod tests { let node_c = graph.add_node(GraphNode::Commodity("C".into())); let node_a = graph.add_node(GraphNode::Commodity("A".into())); let node_b = graph.add_node(GraphNode::Commodity("B".into())); - graph.add_edge(node_c, node_a, GraphEdge::Primary("process1".into())); - graph.add_edge(node_a, node_b, GraphEdge::Primary("process2".into())); + graph.add_edge( + node_c, + node_a, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); + graph.add_edge( + node_a, + node_b, + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, + ); // Validate the graph at DayNight level assert_error!( @@ -326,7 +514,14 @@ mod tests { // Build invalid graph: B(SED) -> A(SED) let node_a = graph.add_node(GraphNode::Commodity("A".into())); let node_b = graph.add_node(GraphNode::Commodity("B".into())); - graph.add_edge(node_b, node_a, GraphEdge::Primary("process1".into())); + graph.add_edge( + node_b, + node_a, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); // Validate the graph at DayNight level assert_error!( @@ -352,8 +547,22 @@ mod tests { let node_a = graph.add_node(GraphNode::Commodity("A".into())); let node_b = graph.add_node(GraphNode::Commodity("B".into())); let node_c = graph.add_node(GraphNode::Commodity("C".into())); - graph.add_edge(node_b, node_a, GraphEdge::Primary("process1".into())); - graph.add_edge(node_a, node_c, GraphEdge::Primary("process2".into())); + graph.add_edge( + node_b, + node_a, + GraphEdge::Primary { + process_id: "process1".into(), + feedback: false, + }, + ); + graph.add_edge( + node_a, + node_c, + GraphEdge::Primary { + process_id: "process2".into(), + feedback: false, + }, + ); // Validate the graph at DayNight level assert_error!( diff --git a/src/input.rs b/src/input.rs index 7e90e5b08..246518310 100644 --- a/src/input.rs +++ b/src/input.rs @@ -1,6 +1,8 @@ //! Common routines for handling input data. use crate::graph::investment::solve_investment_order_for_model; -use crate::graph::validate::validate_commodity_graphs_for_model; +use crate::graph::validate::{ + validate_commodity_graphs_for_model, validate_non_feedback_commodity_graphs_for_model, +}; use crate::graph::{CommoditiesGraph, build_commodity_graphs_for_model}; use crate::id::{HasID, ID}; use crate::model::{Model, ModelParameters}; @@ -294,6 +296,12 @@ pub fn load_model>(model_dir: P) -> Result { &commodities, &time_slice_info, )?; + validate_non_feedback_commodity_graphs_for_model( + &commodity_graphs, + &processes, + &commodities, + &time_slice_info, + )?; // Solve investment order for each region/year let investment_order = solve_investment_order_for_model(&commodity_graphs, &commodities, years); diff --git a/src/input/process.rs b/src/input/process.rs index 0a4f21894..a99636623 100644 --- a/src/input/process.rs +++ b/src/input/process.rs @@ -38,6 +38,7 @@ struct ProcessRaw { end_year: Option, capacity_to_activity: Option, tranche_size: Option, + feedback_process: Option, } define_id_getter! {ProcessRaw, ProcessID} @@ -174,6 +175,7 @@ where capacity_to_activity, investment_constraints: ProcessInvestmentConstraintsMap::new(), tranche_size: process_raw.tranche_size, + feedback_process: process_raw.feedback_process.unwrap_or(false), }; ensure!( diff --git a/src/model/parameters.rs b/src/model/parameters.rs index 6a50ec1e8..b425fff09 100644 --- a/src/model/parameters.rs +++ b/src/model/parameters.rs @@ -115,13 +115,6 @@ pub struct ModelParameters { pub price_tolerance: Dimensionless, /// Number of iterations to perform when calculating prices for cyclically-dependent markets. pub price_cycle_iterations: u32, - /// Slack applied during cycle balancing, allowing newly selected assets to flex their capacity - /// by this proportion. - /// - /// Existing assets remain fixed; this gives newly selected assets the wiggle-room to absorb - /// small demand changes before we would otherwise need to break for re-investment. - #[serde(deserialize_with = "deserialise_finite_non_negative")] - pub capacity_margin: Dimensionless, /// Number of years an asset can remain unused before being decommissioned pub mothball_years: u32, /// Absolute tolerance when checking if remaining demand is close enough to zero @@ -155,7 +148,6 @@ impl Default for ModelParameters { max_ironing_out_iterations: 1, price_tolerance: Dimensionless(1e-6), price_cycle_iterations: 1, - capacity_margin: Dimensionless(0.2), mothball_years: 0, remaining_demand_absolute_tolerance: DEFAULT_REMAINING_DEMAND_ABSOLUTE_TOLERANCE, highs: HighsOptions::default(), @@ -368,8 +360,6 @@ impl ModelParameters { // price_cycle_iterations check_price_cycle_iterations(self.price_cycle_iterations)?; - // capacity_margin already validated with deserialise_finite_non_negative - // remaining_demand_absolute_tolerance already validated with // deserialise_finite_non_negative; check remaining constraints here check_remaining_demand_absolute_tolerance( diff --git a/src/process.rs b/src/process.rs index d592ad95d..e979fb0b2 100644 --- a/src/process.rs +++ b/src/process.rs @@ -70,6 +70,8 @@ pub struct Process { /// how their total capacity is represented as equal-capacity tranches; the number of tranches /// is rounded up when necessary. pub tranche_size: Option, + /// Whether this process participates in feedback loops between commodities. + pub feedback_process: bool, } impl Process { diff --git a/src/simulation.rs b/src/simulation.rs index 914a8a328..cc65a2221 100644 --- a/src/simulation.rs +++ b/src/simulation.rs @@ -15,7 +15,7 @@ use std::sync::Arc; pub mod optimisation; use optimisation::{DispatchRun, FlowMap}; pub mod investment; -use investment::perform_agent_investment; +use investment::{collect_preset_demands_for_year, perform_agent_investment}; pub mod market; pub mod prices; pub use prices::PriceMap; @@ -178,13 +178,15 @@ fn run_dispatch_for_year( debug_assert!(assets.iter().all(|asset| !asset.is_candidate())); debug_assert!(candidates.iter().all(|asset| asset.is_candidate())); + let market_demands = collect_preset_demands_for_year(&model.commodities, year); + // Run dispatch optimisation with existing assets only, if there are any. If not, then assume no // flows (i.e. all are zero) let (solution_existing, flow_map) = if assets.is_empty() { (None, FlowMap::default()) } else { - let solution = - DispatchRun::new(model, assets, year).run("final without candidates", writer)?; + let solution = DispatchRun::new(model, assets, year, &market_demands) + .run("final without candidates", writer)?; let flow_map = solution.create_flow_map(); (Some(solution), flow_map) }; @@ -195,7 +197,7 @@ fn run_dispatch_for_year( None } else { Some( - DispatchRun::new(model, assets, year) + DispatchRun::new(model, assets, year, &market_demands) .with_candidates(candidates) .run("final with candidates", writer)?, ) diff --git a/src/simulation/investment.rs b/src/simulation/investment.rs index 1340c402a..f8175ffe1 100644 --- a/src/simulation/investment.rs +++ b/src/simulation/investment.rs @@ -56,7 +56,8 @@ pub fn perform_agent_investment( writer: &mut DataWriter, ) -> Result> { // Initialise net demand map - let mut net_demand = collect_preset_demands_for_year(&model.commodities, year); + let preset_demands = collect_preset_demands_for_year(&model.commodities, year); + let mut net_demand = preset_demands.clone(); // Keep a list of all the assets selected // This includes Commissioned assets that are selected for retention, and new Ready assets @@ -110,7 +111,7 @@ pub fn perform_agent_investment( // As upstream markets by definition will not yet have producers, we explicitly set // their prices using external values so that they don't appear free - let solution = DispatchRun::new(model, &all_selected_assets, year) + let solution = DispatchRun::new(model, &all_selected_assets, year, &preset_demands) .without_commodity_constraints() .with_market_balance_subset(&seen_markets) .with_input_prices(&prices.shadow) @@ -179,10 +180,14 @@ pub fn update_net_demand_map(demand: &mut AllDemandMap, flows: &FlowMap, assets: let selection = level.containing_selection(time_slice); let key = (commodity_id.clone(), asset.region_id().clone(), selection); // Note: we use the negative of the flow as input flows are negative in the flow map. - demand + let value = demand .entry(key) .and_modify(|value| *value -= *flow) .or_insert(-*flow); + + if *value < Flow(0.0) { + *value = Flow(0.0); + } } } } diff --git a/src/simulation/market.rs b/src/simulation/market.rs index 717a2cf77..909a153a3 100644 --- a/src/simulation/market.rs +++ b/src/simulation/market.rs @@ -1,7 +1,7 @@ //! Code for creating sets of markets. -use super::optimisation::DispatchRun; +use super::optimisation::{DispatchRun, FlowMap}; use crate::agent::Agent; -use crate::asset::{Asset, AssetCapacity, AssetIterator, AssetRef, AssetState}; +use crate::asset::{Asset, AssetCapacity, AssetIterator, AssetRef}; use crate::commodity::{Commodity, CommodityID}; use crate::model::Model; use crate::output::DataWriter; @@ -15,7 +15,6 @@ use crate::simulation::prices::Prices; use crate::time_slice::TimeSliceInfo; use crate::units::{Capacity, Dimensionless, Flow}; use anyhow::{Context, Result}; -use indexmap::IndexMap; use itertools::{Itertools, chain}; use log::debug; use std::collections::HashMap; @@ -80,6 +79,7 @@ impl MarketSet { demand, existing_assets, prices, + &[], writer, ), MarketSet::Cycle(markets) => { @@ -159,6 +159,7 @@ pub fn select_assets_for_single_market( demand: &AllDemandMap, existing_assets: &[AssetRef], prices: &Prices, + excluded_processes: &[ProcessID], writer: &mut DataWriter, ) -> Result> { let commodity = &model.commodities[commodity_id]; @@ -182,7 +183,7 @@ pub fn select_assets_for_single_market( ); // Existing and candidate assets from which to choose - let opt_assets = get_asset_options( + let mut opt_assets = get_asset_options( existing_assets, &demand_portion_for_market, agent, @@ -193,6 +194,9 @@ pub fn select_assets_for_single_market( ) .collect::>(); + // Exclude certain processes from the options + opt_assets.retain(|asset| !excluded_processes.contains(asset.process_id())); + // Calculate the agent's share of addition limits for candidate processes let agent_addition_limits = collect_agent_limits( agent, @@ -233,21 +237,13 @@ pub fn select_assets_for_single_market( Ok(selected_assets) } -/// Iterates through the a pre-ordered set of markets forming a cycle, selecting assets for each +/// Iterates through a pre-ordered set of markets forming a cycle, selecting assets for each /// market in turn. /// -/// Dispatch optimisation is performed after each market is visited to rebalance demand. -/// While dispatching, newly selected (`Ready`) assets are given flexible capacity (bounded by -/// `capacity_margin`) so small demand shifts caused by later markets can be absorbed. After all -/// markets have been visited once, the final set of assets is returned, applying any capacity -/// adjustments from the final full-system dispatch optimisation. +/// Dispatch optimisation is performed after each market is visited. /// -/// Dispatch may fail at any point if new demands are encountered for previously visited markets, -/// and the `capacity_margin` is not sufficient to absorb the demand shift. At this point, the -/// simulation is terminated with an error prompting the user to increase the `capacity_margin`. -/// A longer-term solution (TODO) may be to trigger re-investment for the affected markets. Other -/// yet-to-implement features may also help to stabilise the cycle, such as capacity growth limits. -#[allow(clippy::too_many_arguments)] +/// Dispatch may fail at any point if new demands are encountered for previously visited markets. +#[allow(clippy::too_many_arguments, clippy::too_many_lines)] pub fn select_assets_for_cycle( model: &Model, markets: &[(CommodityID, RegionID)], @@ -255,122 +251,149 @@ pub fn select_assets_for_cycle( demand: &AllDemandMap, existing_assets: &[AssetRef], prices: &Prices, - seen_markets: &[(CommodityID, RegionID)], - previously_selected_assets: &[AssetRef], + _seen_markets: &[(CommodityID, RegionID)], + _previously_selected_assets: &[AssetRef], writer: &mut DataWriter, ) -> Result> { // Precompute a joined string for logging let markets_str = markets.iter().map(|(c, r)| format!("{c}|{r}")).join(", "); - // Iterate over the markets to select assets - let mut current_demand = demand.clone(); - let mut assets_for_cycle = IndexMap::new(); - let mut last_solution = None; - for (idx, (commodity_id, region_id)) in markets.iter().enumerate() { + // Feedback processes are excluded from the second pass while their first-pass flows are + // retained as part of the demand carried into that pass. + let feedback_processes = model + .processes + .iter() + .filter(|(_, process)| process.feedback_process) + .map(|(process_id, _)| process_id.clone()) + .collect::>(); + + // STEP 1 + // Iterate over the markets in order1, considering all processes + let mut net_demand = demand.clone(); + let mut assets_for_first_pass = Vec::new(); + let mut retained_first_pass_assets = Vec::new(); + let mut excluded_flows = FlowMap::new(); + for market in markets { + let (commodity_id, region_id) = market.clone(); + // Select assets for this market - let assets = select_assets_for_single_market( + debug!("Running {commodity_id}|{region_id} selection pass 1"); + let selected_assets = select_assets_for_single_market( model, - commodity_id, - region_id, + &commodity_id, + ®ion_id, year, - ¤t_demand, + &net_demand, existing_assets, prices, + &[], writer, )?; - assets_for_cycle.insert((commodity_id.clone(), region_id.clone()), assets); - - // Assemble full list of assets for dispatch (previously selected + all chosen so far) - let mut all_assets = previously_selected_assets.to_vec(); - let assets_for_cycle_flat: Vec<_> = assets_for_cycle - .values() - .flat_map(|v| v.iter().cloned()) - .collect(); - all_assets.extend_from_slice(&assets_for_cycle_flat); - - // We balance all previously seen markets plus all cycle markets up to and including this one - let mut markets_to_balance = seen_markets.to_vec(); - markets_to_balance.extend_from_slice(&markets[0..=idx]); - - // We allow all `Ready` state assets to have flexible capacity - let flexible_capacity_assets: Vec<_> = assets_for_cycle_flat - .iter() - .filter(|asset| matches!(asset.state(), AssetState::Ready { .. })) - .cloned() - .collect(); - - // Retrieve installable capacity limits for flexible capacity assets. - let mut agent_share_cache = HashMap::new(); - let capacity_limits = flexible_capacity_assets - .iter() - .filter_map(|asset| { - let agent_id = asset.agent_id().unwrap(); - let commodity_id = asset.primary_output_commodity().unwrap(); - let agent_share = *agent_share_cache - .entry((agent_id, commodity_id)) - .or_insert_with(|| { - model.agents[agent_id].commodity_portions[&(commodity_id.clone(), year)] - }); - asset - .process() - .agent_addition_limit(asset.region_id(), asset.commission_year(), agent_share) - .map(|max_capacity| (asset.clone(), max_capacity)) - }) - .collect::>(); + debug!("Completed {commodity_id}|{region_id} selection pass 1"); + + // If no assets have been selected, skip to the next market + if selected_assets.is_empty() { + debug!("No assets selected for '{commodity_id}|{region_id}'"); + continue; + } + assets_for_first_pass.extend(selected_assets.iter().cloned()); + retained_first_pass_assets.extend( + selected_assets + .iter() + .filter(|asset| feedback_processes.contains(asset.process_id())) + .cloned(), + ); // Run dispatch - let solution = DispatchRun::new(model, &all_assets, year) + debug!("Running cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 1"); + let solution = DispatchRun::new(model, &selected_assets, year, &net_demand) .without_commodity_constraints() - .with_market_balance_subset(&markets_to_balance) - .with_flexible_capacity_assets( - &flexible_capacity_assets, - Some(&capacity_limits), - // Gives newly selected cycle assets limited capacity wiggle-room; existing assets stay fixed. - model.parameters.capacity_margin, + .with_market_balance_subset(std::slice::from_ref(market)) + .with_input_prices(&prices.shadow) + .run( + &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 1"), + writer, ) + .with_context(|| format!("Dispatch failed for cycle ({markets_str})"))?; + debug!("Completed cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 1"); + + // Update demand map with flows from newly selected assets + let flows = solution.create_flow_map(); + for ((asset, commodity_id, time_slice), flow) in &flows { + if feedback_processes.contains(asset.process_id()) { + excluded_flows + .entry((asset.clone(), commodity_id.clone(), time_slice.clone())) + .and_modify(|stored_flow| *stored_flow += *flow) + .or_insert(*flow); + } + } + update_net_demand_map(&mut net_demand, &flows, &selected_assets); + } + + // STEP 2 + // Iterate over the markets in order2, excluding feedback processes + let mut net_demand = demand.clone(); + update_net_demand_map(&mut net_demand, &excluded_flows, &assets_for_first_pass); + let mut assets_for_second_pass = Vec::new(); + for market in markets { + let (commodity_id, region_id) = market.clone(); + + // Select assets for this market + debug!("Running {commodity_id}|{region_id} selection pass 2"); + let selected_assets = select_assets_for_single_market( + model, + &commodity_id, + ®ion_id, + year, + &net_demand, + existing_assets, + prices, + &feedback_processes, + writer, + )?; + debug!("Completed {commodity_id}|{region_id} selection pass 2"); + + // If no assets have been selected, skip to the next market + if selected_assets.is_empty() { + debug!("No assets selected for '{commodity_id}|{region_id}'"); + continue; + } + assets_for_second_pass.extend(selected_assets.iter().cloned()); + + // // Assemble full list of assets for dispatch (previously selected + all chosen so far) + // let mut all_assets = previously_selected_assets.to_vec(); + // all_assets.extend(assets_for_first_pass.iter().cloned()); + // all_assets.extend(assets_for_second_pass.iter().cloned()); + + // // We balance all previously seen markets plus all cycle markets up to and including this one + // let mut markets_to_balance = seen_markets.to_vec(); + // markets_to_balance.extend_from_slice(&investment_order[0..=idx]); + + // Run dispatch + debug!("Running cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 2"); + let solution = DispatchRun::new(model, &selected_assets, year, &net_demand) + .without_commodity_constraints() + .with_market_balance_subset(std::slice::from_ref(market)) + .with_input_prices(&prices.shadow) .run( - &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment"), + &format!("cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 2"), writer, ) - .with_context(|| { - format!( - "Cycle balancing failed for cycle ({markets_str}), capacity_margin: {}. \ - Try increasing the capacity_margin.", - model.parameters.capacity_margin - ) - })?; + .with_context(|| format!("Dispatch failed for cycle ({markets_str})"))?; + debug!("Completed cycle ({markets_str}) post {commodity_id}|{region_id} investment pass 2"); - // Calculate new net demand map with all assets selected so far - current_demand.clone_from(demand); + // Update demand map with flows from newly selected assets update_net_demand_map( - &mut current_demand, + &mut net_demand, &solution.create_flow_map(), - &assets_for_cycle_flat, + &selected_assets, ); - last_solution = Some(solution); } - // Finally, update flexible capacity assets based on the final solution - let mut all_cycle_assets: Vec<_> = assets_for_cycle.into_values().flatten().collect(); - if let Some(solution) = last_solution { - let new_capacities: HashMap<_, _> = solution.iter_capacity().collect(); - for asset in &mut all_cycle_assets { - if let Some(new_capacity) = new_capacities.get(asset) { - debug!( - "Capacity of asset '{}' modified during cycle balancing ({} to {})", - asset.process_id(), - asset.total_capacity(), - new_capacity.total_capacity() - ); - asset.make_mut().set_capacity(*new_capacity); - } - } - } - - // Drop any assets who's capacities were dropped to zero - all_cycle_assets.retain(|asset| asset.num_tranches() > 0); - - Ok(all_cycle_assets) + Ok(retained_first_pass_assets + .into_iter() + .chain(assets_for_second_pass) + .collect()) } /// Get a portion of the demand profile for this market diff --git a/src/simulation/optimisation.rs b/src/simulation/optimisation.rs index 3b42e591e..5e63af895 100644 --- a/src/simulation/optimisation.rs +++ b/src/simulation/optimisation.rs @@ -1,25 +1,21 @@ //! Code for performing dispatch optimisation. //! //! This is used to calculate commodity flows and prices. -use crate::asset::{Asset, AssetCapacity, AssetRef, AssetState}; +use crate::asset::{Asset, AssetRef}; use crate::commodity::CommodityID; -use crate::finance::annual_capital_cost; use crate::input::format_items_with_cap; use crate::model::Model; use crate::output::DataWriter; use crate::region::RegionID; use crate::simulation::PriceMap; +use crate::simulation::investment::AllDemandMap; use crate::time_slice::{TimeSliceID, TimeSliceInfo, TimeSliceSelection}; -use crate::units::{ - Activity, Capacity, Dimensionless, Flow, Money, MoneyPerActivity, MoneyPerCapacity, - MoneyPerFlow, Year, -}; +use crate::units::{Activity, Flow, Money, MoneyPerActivity, MoneyPerFlow}; use anyhow::{Context, Result, anyhow, bail}; use highs::{HighsModelStatus, RowProblem as Problem, Sense}; use indexmap::{IndexMap, IndexSet}; use itertools::{chain, iproduct}; use log::warn; -use std::collections::HashMap; use std::error::Error; use std::ops::Range; @@ -38,9 +34,6 @@ type Variable = highs::Col; /// The map of activity variables for assets type ActivityVariableMap = IndexMap<(AssetRef, TimeSliceID), Variable>; -/// A map of capacity variables for assets -type CapacityVariableMap = IndexMap; - /// Variables representing unmet demand for a given market type UnmetDemandVariableMap = IndexMap<(CommodityID, RegionID, TimeSliceID), Variable>; @@ -57,8 +50,6 @@ pub struct VariableMap { activity_vars: ActivityVariableMap, existing_asset_var_idx: Range, candidate_asset_var_idx: Range, - capacity_vars: CapacityVariableMap, - capacity_var_idx: Range, unmet_demand_vars: UnmetDemandVariableMap, unmet_demand_var_idx: Range, } @@ -104,8 +95,6 @@ impl VariableMap { activity_vars, existing_asset_var_idx, candidate_asset_var_idx, - capacity_vars: CapacityVariableMap::new(), - capacity_var_idx: Range::default(), unmet_demand_vars: UnmetDemandVariableMap::default(), unmet_demand_var_idx: Range::default(), } @@ -173,11 +162,6 @@ impl VariableMap { fn activity_var_keys(&self) -> indexmap::map::Keys<'_, (AssetRef, TimeSliceID), Variable> { self.activity_vars.keys() } - - /// Iterate over capacity variables - fn iter_capacity_vars(&self) -> impl Iterator { - self.capacity_vars.iter().map(|(asset, var)| (asset, *var)) - } } /// The solution to the dispatch optimisation problem @@ -261,22 +245,6 @@ impl Solution<'_> { }) } - /// Iterate over capacity values - pub fn iter_capacity(&self) -> impl Iterator { - self.variables - .capacity_vars - .keys() - .zip(self.solution.columns()[self.variables.capacity_var_idx.clone()].iter()) - .map(|(asset, capacity_var)| { - // The capacity variable represents number of tranches - let tranche_size = asset.capacity().tranche_size(); - #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] - let asset_capacity = AssetCapacity::new(capacity_var.round() as u32, tranche_size); - - (asset, asset_capacity) - }) - } - /// Keys and dual values for commodity balance constraints. pub fn iter_commodity_balance_duals( &self, @@ -296,11 +264,7 @@ impl Solution<'_> { /// Keys and dual values for activity constraints. /// - /// Note: if there are any flexible capacity assets, these will have two duals with identical - /// keys, and there will be no way to distinguish between them in the resulting iterator. - /// Recommended for now only to use this function when there are no flexible capacity assets. - /// - /// Also note: this excludes seasonal and annual constraints. Recommended for now not to use + /// Note: this excludes seasonal and annual constraints. Recommended for now not to use /// this for models that include seasonal or annual availability constraints. pub fn iter_activity_duals( &self, @@ -436,45 +400,31 @@ fn filter_input_prices( pub struct DispatchRun<'model, 'run> { model: &'model Model, existing_assets: &'run [AssetRef], - flexible_capacity_assets: &'run [AssetRef], - capacity_limits: Option<&'run HashMap>, candidate_assets: &'run [AssetRef], markets_to_balance: &'run [(CommodityID, RegionID)], + market_demands: &'run AllDemandMap, input_prices: Option<&'run PriceMap>, include_commodity_constraints: bool, year: u32, - capacity_margin: Dimensionless, } impl<'model, 'run> DispatchRun<'model, 'run> { /// Create a new [`DispatchRun`] for the specified model and assets for a given year - pub fn new(model: &'model Model, assets: &'run [AssetRef], year: u32) -> Self { + pub fn new( + model: &'model Model, + assets: &'run [AssetRef], + year: u32, + market_demands: &'run AllDemandMap, + ) -> Self { Self { model, existing_assets: assets, - flexible_capacity_assets: &[], - capacity_limits: None, candidate_assets: &[], markets_to_balance: &[], + market_demands, input_prices: None, include_commodity_constraints: true, year, - capacity_margin: Dimensionless(0.0), - } - } - - /// Include the specified flexible capacity assets in the dispatch run - pub fn with_flexible_capacity_assets( - self, - flexible_capacity_assets: &'run [AssetRef], - capacity_limits: Option<&'run HashMap>, - capacity_margin: Dimensionless, - ) -> Self { - Self { - flexible_capacity_assets, - capacity_limits, - capacity_margin, - ..self } } @@ -741,25 +691,6 @@ impl<'model, 'run> DispatchRun<'model, 'run> { variables.add_unmet_demand_variables(&mut problem, self.model, markets_to_balance); } - // Check flexible capacity assets is a subset of existing assets - for asset in self.flexible_capacity_assets { - assert!( - self.existing_assets.contains(asset), - "Flexible capacity assets must be a subset of existing assets. Offending asset: {asset:?}" - ); - } - - // Add capacity variables for flexible capacity assets - if !self.flexible_capacity_assets.is_empty() { - variables.capacity_var_idx = add_capacity_variables( - &mut problem, - &mut variables.capacity_vars, - self.flexible_capacity_assets, - self.capacity_limits, - self.capacity_margin, - ); - } - // Add constraints let all_assets = chain(self.existing_assets.iter(), self.candidate_assets.iter()); let constraint_keys = add_model_constraints( @@ -768,6 +699,7 @@ impl<'model, 'run> DispatchRun<'model, 'run> { self.model, &all_assets, markets_to_balance, + self.market_demands, self.year, self.candidate_assets, include_commodity_constraints, @@ -821,51 +753,6 @@ fn add_activity_variables( start..problem.num_cols() } -fn add_capacity_variables( - problem: &mut Problem, - variables: &mut CapacityVariableMap, - assets: &[AssetRef], - capacity_limits: Option<&HashMap>, - capacity_margin: Dimensionless, -) -> Range { - let capacity_margin = capacity_margin.value(); - - // This line **must** come before we add more variables - let start = problem.num_cols(); - - for asset in assets { - // Can only have flexible capacity for `Ready` assets - assert!( - matches!(asset.state(), AssetState::Ready { .. }), - "Flexible capacity can only be assigned to `Ready` type assets. Offending asset: {asset:?}" - ); - - // Coefficient: cost per capacity - let coeff = calculate_capacity_coefficient(asset); - - // Add a capacity variable for each asset - // Bounds are calculated based on current capacity with wiggle-room defined by - // `capacity_margin`, and limited by `capacity_limit` if provided. - // Since capacity variables are numbers of tranches, we apply constraints to the tranche count - let tranche_size = asset.capacity().tranche_size(); - let current_tranches = asset.capacity().num_tranches(); - - let lower = (current_tranches as f64 * (1.0 - capacity_margin)).max(0.0); - - let mut upper = current_tranches as f64 * (1.0 + capacity_margin); - if let Some(limit) = capacity_limits.and_then(|limits| limits.get(asset)) { - upper = upper.min((*limit / tranche_size).value()); - } - - let var = problem.add_integer_column((coeff * tranche_size).value(), lower..=upper); - - let existing = variables.insert(asset.clone(), var).is_some(); - assert!(!existing, "Duplicate entry for var"); - } - - start..problem.num_cols() -} - /// Calculate the cost coefficient for an activity variable. /// /// Normally, the cost coefficient is the same as the asset's operating costs for the given year and @@ -895,13 +782,3 @@ fn calculate_activity_coefficient( opex } } - -/// Calculate the cost coefficient for a capacity variable (for flexible capacity assets only). -/// -/// This includes both the annual fixed operating cost and the annual capital cost. -fn calculate_capacity_coefficient(asset: &AssetRef) -> MoneyPerCapacity { - let param = asset.process_parameter(); - let annual_fixed_operating_cost = param.fixed_operating_cost * Year(1.0); - annual_fixed_operating_cost - + annual_capital_cost(param.capital_cost, param.lifetime, param.discount_rate) -} diff --git a/src/simulation/optimisation/constraints.rs b/src/simulation/optimisation/constraints.rs index fd233408a..07620ec36 100644 --- a/src/simulation/optimisation/constraints.rs +++ b/src/simulation/optimisation/constraints.rs @@ -5,11 +5,12 @@ use crate::commodity::{BalanceType, CommodityID, CommodityType}; use crate::model::Model; use crate::process::FlowDirection; use crate::region::RegionID; +use crate::simulation::investment::AllDemandMap; use crate::time_slice::{Season, TimeSliceInfo, TimeSliceSelection}; use crate::units::{Flow, MoneyPerCapacityPerYear, UnitType, Year}; use highs::RowProblem as Problem; use indexmap::IndexMap; -use std::collections::{HashMap, HashSet}; +use std::collections::HashMap; /// Corresponding variables for a constraint along with the row offset in the solution pub struct KeysWithOffset { @@ -88,6 +89,7 @@ pub fn add_model_constraints<'a, I>( model: &'a Model, assets: &I, markets_to_balance: &'a [(CommodityID, RegionID)], + market_demands: &AllDemandMap, year: u32, candidate_assets: &'a [AssetRef], include_commodity_constraints: bool, @@ -101,7 +103,7 @@ where model, assets, markets_to_balance, - year, + market_demands, candidate_assets, ); @@ -330,7 +332,7 @@ fn add_commodity_balance_constraints<'a, I>( model: &'a Model, assets: &I, markets_to_balance: &'a [(CommodityID, RegionID)], - year: u32, + market_demands: &AllDemandMap, candidate_assets: &'a [AssetRef], ) -> CommodityBalanceKeys where @@ -394,14 +396,18 @@ where model.parameters.commodity_balance_epsilon, ); - // For SVD commodities, the lower bound is the exogenous demand (or epsilon if larger). - // For SED commodities, the lower bound is just epsilon. - let min = match commodity.kind { - CommodityType::ServiceDemand => { - commodity.demand[&(region_id.clone(), year, ts_selection.clone())].max(epsilon) - } - _ => epsilon, + // For SVD commodities, the demand must be present in the map; for SED commodities, + // missing entries mean no demand for this selection. + let key = ( + commodity_id.clone(), + region_id.clone(), + ts_selection.clone(), + ); + let demand_for_selection = match commodity.kind { + CommodityType::ServiceDemand => market_demands[&key], + _ => market_demands.get(&key).copied().unwrap_or(Flow(0.0)), }; + let min = demand_for_selection.max(epsilon); // Consume collected terms into a row. `terms.drain(..)` ensures the vector is // emptied for the next selection. @@ -456,9 +462,7 @@ fn candidate_balance_epsilon( /// /// Returns an `ActivityKeys` where `offset` is the row index of the first /// activity constraint added and `keys` enumerates the `(asset, time_selection)` -/// entries in the same row order. Note that for flexible-capacity assets two rows -/// (upper and lower bounds) are added per selection; in that case the same key is -/// stored twice to match the solver ordering. +/// entries in the same row order. /// #[doc = concat!("[1]: ", crate::docs_url!("model/dispatch_optimisation.html#asset-activity-limits"))] fn add_activity_constraints<'a, I>( @@ -476,61 +480,23 @@ where let offset = problem.num_rows(); let mut keys = Vec::new(); - let capacity_vars: IndexMap<&AssetRef, highs::Col> = variables.iter_capacity_vars().collect(); // Create constraints for each asset for asset in assets { - if let Some(&capacity_var) = capacity_vars.get(asset) { - // Asset with flexible capacity - for (ts_selection, limits) in asset.iter_activity_per_capacity_limits() { - let mut upper_limit = limits.end().value(); - let mut lower_limit = limits.start().value(); - - // The capacity variable represents number of tranches, so we need to multiply the - // per-capacity limits by the tranche size. - let tranche_size = asset.capacity().tranche_size(); - upper_limit *= tranche_size.value(); - lower_limit *= tranche_size.value(); - - // Collect capacity and activity terms - // We have a single capacity term, and activity terms for all time slices in the selection - let mut terms_upper = vec![(capacity_var, -upper_limit)]; - let mut terms_lower = vec![(capacity_var, -lower_limit)]; - for (time_slice, _) in ts_selection.iter(time_slice_info) { - let var = variables.get_activity_var(asset, time_slice); - terms_upper.push((var, 1.0)); - terms_lower.push((var, 1.0)); - } + for (ts_selection, limits) in asset.iter_activity_limits() { + let limits = limits.start().value()..=limits.end().value(); - // Upper bound: sum(activity) - (capacity * upper_limit_per_capacity) ≤ 0 - problem.add_row(..=0.0, &terms_upper); - - // Lower bound: sum(activity) - (capacity * lower_limit_per_capacity) ≥ 0 - problem.add_row(0.0.., &terms_lower); - - // Store keys for retrieving duals later. - // TODO: a bit of a hack pushing identical keys twice. Safe for now so long as we don't - // use the activity duals for anything important when using flexible capacity assets. - keys.push((asset.clone(), ts_selection.clone())); - keys.push((asset.clone(), ts_selection.clone())); - } - } else { - // Fixed-capacity asset: simple absolute activity limits. - for (ts_selection, limits) in asset.iter_activity_limits() { - let limits = limits.start().value()..=limits.end().value(); - - // Collect activity terms for the time slices in this selection - let terms = ts_selection - .iter(time_slice_info) - .map(|(time_slice, _)| (variables.get_activity_var(asset, time_slice), 1.0)) - .collect::>(); + // Collect activity terms for the time slices in this selection + let terms = ts_selection + .iter(time_slice_info) + .map(|(time_slice, _)| (variables.get_activity_var(asset, time_slice), 1.0)) + .collect::>(); - // Constraint: sum of activities in selection within limits - problem.add_row(limits, &terms); + // Constraint: sum of activities in selection within limits + problem.add_row(limits, &terms); - // Store keys for retrieving duals later. - keys.push((asset.clone(), ts_selection.clone())); - } + // Store keys for retrieving duals later. + keys.push((asset.clone(), ts_selection.clone())); } } @@ -544,7 +510,7 @@ where /// which is the authoritative check for equivalence. This also handles hash collisions correctly. /// /// The caller must ensure that `assets` contains only assets eligible for equal-utilisation -/// constraints (i.e. flexible-capacity assets have already been filtered out). +/// constraints. fn group_dispatch_equivalent_assets<'a, I>(assets: I) -> Vec> where I: Iterator, @@ -581,8 +547,7 @@ where /// Add constraints requiring dispatch-equivalent assets to have equal utilisation in each time /// slice. /// -/// Flexible-capacity assets are excluded because their maximum activity depends on a decision -/// variable. The constraints added here are not included in [`ConstraintKeys`], as their duals +/// The constraints added here are not included in [`ConstraintKeys`], as their duals /// are not currently used. fn add_equal_utilisation_constraints<'a, I>( problem: &mut Problem, @@ -592,14 +557,7 @@ fn add_equal_utilisation_constraints<'a, I>( ) where I: Iterator + 'a, { - // Identify flexible-capacity assets so we can exclude them from the constraints - let flexible_assets: HashSet<_> = variables - .iter_capacity_vars() - .map(|(asset, _)| asset) - .collect(); - - let asset_groups = - group_dispatch_equivalent_assets(assets.filter(|asset| !flexible_assets.contains(asset))); + let asset_groups = group_dispatch_equivalent_assets(assets); // For each group of assets, add constraints to force equal utilisation in each time slice // This is done by anchoring each asset to the first asset in the group (-> (n-1) constraints diff --git a/src/simulation/prices.rs b/src/simulation/prices.rs index 87be3f0e6..b2abe76a3 100644 --- a/src/simulation/prices.rs +++ b/src/simulation/prices.rs @@ -1408,6 +1408,7 @@ mod tests { capacity_to_activity: ActivityPerCapacity(1.0), investment_constraints: HashMap::new(), tranche_size: None, + feedback_process: false, } }