From cfbf3651a336e48886a34ce14c2583c718f9dca6 Mon Sep 17 00:00:00 2001 From: James Tomlinson Date: Mon, 15 Jul 2024 20:14:54 +0100 Subject: [PATCH 1/3] feat: Utilise constant parameters in aggregated node LP creation. Where possible utilise constant parameters to try to get aggregated node factor pairs during LP creation. This makes it possible to support constant factors in LP solvers without an API for updating the constraint matrix. It should also be generally more efficient by only updating non-constant factors. --- pywr-core/src/aggregated_node.rs | 229 +++++++++++++++++- pywr-core/src/metric.rs | 31 +++ pywr-core/src/models/multi.rs | 2 +- pywr-core/src/models/simple.rs | 2 +- pywr-core/src/network.rs | 17 +- pywr-core/src/solvers/builder.rs | 90 ++++--- pywr-core/src/solvers/cbc/mod.rs | 16 +- pywr-core/src/solvers/clp/mod.rs | 11 +- pywr-core/src/solvers/highs/mod.rs | 30 ++- pywr-core/src/solvers/mod.rs | 6 +- pywr-core/src/test_utils.rs | 11 +- .../src/nodes/water_treatment_works.rs | 2 +- 12 files changed, 394 insertions(+), 53 deletions(-) diff --git a/pywr-core/src/aggregated_node.rs b/pywr-core/src/aggregated_node.rs index ddcc3c1e..37e68b3a 100644 --- a/pywr-core/src/aggregated_node.rs +++ b/pywr-core/src/aggregated_node.rs @@ -1,7 +1,7 @@ use crate::metric::MetricF64; use crate::network::Network; use crate::node::{Constraint, FlowConstraints, NodeMeta}; -use crate::state::State; +use crate::state::{ConstParameterValues, State}; use crate::{NodeIndex, PywrError}; use std::ops::{Deref, DerefMut}; @@ -64,6 +64,16 @@ pub enum Factors { Ratio(Vec), } +impl Factors { + /// Returns true if all factors are constant + pub fn is_constant(&self) -> bool { + match self { + Factors::Proportion(factors) => factors.iter().all(|f| f.is_constant()), + Factors::Ratio(factors) => factors.iter().all(|f| f.is_constant()), + } + } +} + #[derive(Debug, PartialEq)] pub struct AggregatedNode { meta: NodeMeta, @@ -102,6 +112,41 @@ impl NodeFactorPair { } } +/// A constant node factor. If the factor is non-constant, the factor value here is `None`. +#[derive(Debug, PartialEq, Copy, Clone)] +pub struct NodeConstFactor { + pub index: NodeIndex, + pub factor: Option, +} + +impl NodeConstFactor { + fn new(node: NodeIndex, factor: Option) -> Self { + Self { index: node, factor } + } +} + +/// A pair of nodes and their factors +#[derive(Debug, PartialEq, Copy, Clone)] +pub struct NodeConstFactorPair { + pub node0: NodeConstFactor, + pub node1: NodeConstFactor, +} + +impl NodeConstFactorPair { + fn new(node0: NodeConstFactor, node1: NodeConstFactor) -> Self { + Self { node0, node1 } + } + + /// Return the ratio of the two factors (node0 / node1). If either factor is `None`, + /// the ratio is also `None`. + pub fn ratio(&self) -> Option { + match (self.node0.factor, self.node1.factor) { + (Some(f0), Some(f1)) => Some(f0 / f1), + _ => None, + } + } +} + impl AggregatedNode { pub fn new( index: &AggregatedNodeIndex, @@ -140,6 +185,15 @@ impl AggregatedNode { self.nodes.to_vec() } + /// Does the aggregated node have factors? + pub fn has_factors(&self) -> bool { + self.factors.is_some() + } + + /// Does the aggregated node have constant factors? + pub fn has_const_factors(&self) -> bool { + self.factors.as_ref().map(|f| f.is_constant()).unwrap_or(false) + } pub fn set_factors(&mut self, factors: Option) { self.factors = factors; } @@ -148,6 +202,7 @@ impl AggregatedNode { self.factors.as_ref() } + /// Return normalised factor pairs pub fn get_factor_node_pairs(&self) -> Option> { if self.factors.is_some() { let n0 = self.nodes[0]; @@ -158,6 +213,21 @@ impl AggregatedNode { } } + /// Return constant normalised factor pairs + pub fn get_const_norm_factor_pairs(&self, values: &ConstParameterValues) -> Option> { + if let Some(factors) = &self.factors { + let pairs = match factors { + Factors::Proportion(prop_factors) => { + get_const_norm_proportional_factor_pairs(prop_factors, &self.nodes, values) + } + Factors::Ratio(ratio_factors) => get_const_norm_ratio_factor_pairs(ratio_factors, &self.nodes, values), + }; + Some(pairs) + } else { + None + } + } + /// Return normalised factor pairs /// pub fn get_norm_factor_pairs(&self, model: &Network, state: &State) -> Option> { @@ -225,7 +295,11 @@ impl AggregatedNode { } } -/// Proportional factors +/// Calculate factor pairs for proportional factors. +/// +/// There should be one less factor than node indices. The factors correspond to each of the node +/// indices after the first. Factor pairs relating the first index to each of the other indices are +/// calculated. This requires the sum of the factors to be greater than 0.0 and less than 1.0. fn get_norm_proportional_factor_pairs( factors: &[MetricF64], nodes: &[NodeIndex], @@ -263,7 +337,71 @@ fn get_norm_proportional_factor_pairs( .collect::>() } -/// Ratio factors +/// Calculate constant factor pairs for proportional factors. +/// +/// There should be one less factor than node indices. The factors correspond to each of the node +/// indices after the first. Factor pairs relating the first index to each of the other indices are +/// calculated. This requires the sum of the factors to be greater than 0.0 and less than 1.0. If +/// any of the factors are not constant, the factor pairs will contain `None` values. +fn get_const_norm_proportional_factor_pairs( + factors: &[MetricF64], + nodes: &[NodeIndex], + values: &ConstParameterValues, +) -> Vec { + if factors.len() != nodes.len() - 1 { + panic!("Found {} proportional factors and {} nodes in aggregated node. The number of proportional factors should equal one less than the number of nodes.", factors.len(), nodes.len()); + } + + // First get the current factor values + let values: Vec> = factors + .iter() + .map(|f| f.try_get_constant_value(values)) + .collect::, PywrError>>() + .expect("Failed to get current factor values."); + + let n0 = nodes[0]; + + // To calculate the factors we require that every factor is available. + if values.iter().any(|v| v.is_none()) { + // At least one factor is not available; therefore we can not calculate "f0" + nodes + .iter() + .skip(1) + .zip(values) + .map(move |(&n1, f1)| { + NodeConstFactorPair::new(NodeConstFactor::new(n0, None), NodeConstFactor::new(n1, f1)) + }) + .collect::>() + } else { + // All factors are available; therefore we can calculate "f0" + // TODO do we need to assert that each individual factor is positive? + let total: f64 = values + .iter() + .map(|v| v.expect("Factor is `None`; this should be impossible.")) + .sum(); + if total < 0.0 { + panic!("Proportional factors are too small or negative."); + } + if total >= 1.0 { + panic!("Proportional factors are too large.") + } + + let f0 = Some(1.0 - total); + + nodes + .iter() + .skip(1) + .zip(values) + .map(move |(&n1, f1)| NodeConstFactorPair::new(NodeConstFactor::new(n0, f0), NodeConstFactor::new(n1, f1))) + .collect::>() + } +} + +/// Calculate factor pairs for ratio factors. +/// +/// The number of node indices and factors should be equal. The factors correspond to each of the +/// node indices. Factor pairs relating the first index to each of the other indices are calculated. +/// This requires that the factors are all non-zero. fn get_norm_ratio_factor_pairs( factors: &[MetricF64], nodes: &[NodeIndex], @@ -290,12 +428,45 @@ fn get_norm_ratio_factor_pairs( .collect::>() } +/// Constant ratio factors using constant values if they are available. If they are not available, +/// the factors are `None`. +fn get_const_norm_ratio_factor_pairs( + factors: &[MetricF64], + nodes: &[NodeIndex], + values: &ConstParameterValues, +) -> Vec { + if factors.len() != nodes.len() { + panic!("Found {} ratio factors and {} nodes in aggregated node. The number of ratio factors should equal the number of nodes.", factors.len(), nodes.len()); + } + + let n0 = nodes[0]; + // Try to convert the factor into a constant + + let f0 = factors[0] + .try_get_constant_value(values) + .unwrap_or_else(|e| panic!("Failed to get constant value for factor: {}", e)); + + nodes + .iter() + .zip(factors) + .skip(1) + .map(move |(&n1, f1)| { + let v1 = f1 + .try_get_constant_value(values) + .unwrap_or_else(|e| panic!("Failed to get constant value for factor: {}", e)); + + NodeConstFactorPair::new(NodeConstFactor::new(n0, f0), NodeConstFactor::new(n1, v1)) + }) + .collect::>() +} + #[cfg(test)] mod tests { use crate::aggregated_node::Factors; use crate::metric::MetricF64; use crate::models::Model; use crate::network::Network; + use crate::parameters::MonthlyProfileParameter; use crate::recorders::AssertionRecorder; use crate::test_utils::{default_time_domain, run_all_solvers}; use ndarray::Array2; @@ -344,6 +515,56 @@ mod tests { let model = Model::new(default_time_domain().into(), network); - run_all_solvers(&model, &["cbc", "highs"]); + run_all_solvers(&model, &[]); + } + + /// Test the factors forcing a simple ratio of flow that varies over time + /// + /// The model has a single input that diverges to two links and respective output nodes. + #[test] + fn test_simple_factor_profile() { + let mut network = Network::default(); + + let input_node = network.add_input_node("input", None).unwrap(); + let link_node0 = network.add_link_node("link", Some("0")).unwrap(); + let output_node0 = network.add_output_node("output", Some("0")).unwrap(); + + network.connect_nodes(input_node, link_node0).unwrap(); + network.connect_nodes(link_node0, output_node0).unwrap(); + + let link_node1 = network.add_link_node("link", Some("1")).unwrap(); + let output_node1 = network.add_output_node("output", Some("1")).unwrap(); + + network.connect_nodes(input_node, link_node1).unwrap(); + network.connect_nodes(link_node1, output_node1).unwrap(); + + let factor_profile = MonthlyProfileParameter::new("factor-profile", [2.0; 12], None); + let factor_profile_idx = network.add_simple_parameter(Box::new(factor_profile)).unwrap(); + + let factors = Some(Factors::Ratio(vec![factor_profile_idx.into(), 1.0.into()])); + + let _agg_node = network.add_aggregated_node("agg-node", None, &[link_node0, link_node1], factors); + + // Setup a demand on output-0 + let output_node = network.get_mut_node_by_name("output", Some("0")).unwrap(); + output_node.set_max_flow_constraint(Some(100.0.into())).unwrap(); + + output_node.set_cost(Some((-10.0).into())); + + // Set-up assertion for "input" node + let idx = network.get_node_by_name("link", Some("0")).unwrap().index(); + let expected = Array2::from_elem((366, 10), 100.0); + let recorder = AssertionRecorder::new("link-0-flow", MetricF64::NodeOutFlow(idx), expected, None, None); + network.add_recorder(Box::new(recorder)).unwrap(); + + // Set-up assertion for "input" node + let idx = network.get_node_by_name("link", Some("1")).unwrap().index(); + let expected = Array2::from_elem((366, 10), 50.0); + let recorder = AssertionRecorder::new("link-0-flow", MetricF64::NodeOutFlow(idx), expected, None, None); + network.add_recorder(Box::new(recorder)).unwrap(); + + let model = Model::new(default_time_domain().into(), network); + + run_all_solvers(&model, &["cbc"]); } } diff --git a/pywr-core/src/metric.rs b/pywr-core/src/metric.rs index 523b6cbc..47449e0c 100644 --- a/pywr-core/src/metric.rs +++ b/pywr-core/src/metric.rs @@ -45,6 +45,22 @@ impl SimpleMetricF64 { SimpleMetricF64::Constant(m) => m.get_value(&values.get_constant_values()), } } + + /// Try to get the constant value of the metric, if it is a constant value. + pub fn try_get_constant_value(&self, values: &ConstParameterValues) -> Result, PywrError> { + match self { + SimpleMetricF64::Constant(c) => c.get_value(values).map(Some), + _ => Ok(None), + } + } + + /// Returns true if the metric is a constant value. + pub fn is_constant(&self) -> bool { + match self { + SimpleMetricF64::Constant(_) => true, + _ => false, + } + } } #[derive(Clone, Debug, PartialEq)] @@ -124,6 +140,21 @@ impl MetricF64 { MetricF64::Simple(s) => s.get_value(&state.get_simple_parameter_values()), } } + + /// Try to get the constant value of the metric, if it is a constant value. + pub fn try_get_constant_value(&self, values: &ConstParameterValues) -> Result, PywrError> { + match self { + MetricF64::Simple(s) => s.try_get_constant_value(values), + _ => Ok(None), + } + } + + pub fn is_constant(&self) -> bool { + match self { + MetricF64::Simple(s) => s.is_constant(), + _ => false, + } + } } impl TryFrom for SimpleMetricF64 { diff --git a/pywr-core/src/models/multi.rs b/pywr-core/src/models/multi.rs index 316c572b..6e1dcad8 100644 --- a/pywr-core/src/models/multi.rs +++ b/pywr-core/src/models/multi.rs @@ -159,7 +159,7 @@ impl MultiNetworkModel { .network .setup_network(timesteps, scenario_indices, entry.parameters.len())?; let recorder_state = entry.network.setup_recorders(&self.domain)?; - let solver = entry.network.setup_solver::(scenario_indices, settings)?; + let solver = entry.network.setup_solver::(scenario_indices, &state, settings)?; states.push(state); recorder_states.push(recorder_state); diff --git a/pywr-core/src/models/simple.rs b/pywr-core/src/models/simple.rs index 5017f9f1..a1b382b8 100644 --- a/pywr-core/src/models/simple.rs +++ b/pywr-core/src/models/simple.rs @@ -78,7 +78,7 @@ impl Model { let state = self.network.setup_network(timesteps, scenario_indices, 0)?; let recorder_state = self.network.setup_recorders(&self.domain)?; - let solvers = self.network.setup_solver::(scenario_indices, settings)?; + let solvers = self.network.setup_solver::(scenario_indices, &state, settings)?; Ok(ModelState { current_time_step_idx: 0, diff --git a/pywr-core/src/network.rs b/pywr-core/src/network.rs index f0f55ca7..f0b3a9d0 100644 --- a/pywr-core/src/network.rs +++ b/pywr-core/src/network.rs @@ -306,6 +306,7 @@ impl Network { pub fn setup_solver( &self, scenario_indices: &[ScenarioIndex], + state: &NetworkState, settings: &S::Settings, ) -> Result>, PywrError> where @@ -317,9 +318,10 @@ impl Network { let mut solvers = Vec::with_capacity(scenario_indices.len()); - for _scenario_index in scenario_indices { + for scenario_index in scenario_indices { // Create a solver for each scenario - let solver = S::setup(self, settings)?; + let const_values = state.state(scenario_index).get_const_parameter_values(); + let solver = S::setup(self, &const_values, settings)?; solvers.push(solver); } @@ -557,10 +559,19 @@ impl Network { } // Aggregated node factors required if any aggregated node has factors defined. - if self.aggregated_nodes.iter().any(|n| n.get_factors().is_some()) { + if self.aggregated_nodes.iter().any(|n| n.has_factors()) { features.insert(SolverFeatures::AggregatedNodeFactors); } + // Aggregated node dynamic factors required if any aggregated node has dynamic factors defined. + if self + .aggregated_nodes + .iter() + .any(|n| n.has_factors() && !n.has_const_factors()) + { + features.insert(SolverFeatures::AggregatedNodeDynamicFactors); + } + // The presence of any virtual storage node requires the VirtualStorage feature. if self.virtual_storage_nodes.len() > 0 { features.insert(SolverFeatures::VirtualStorage); diff --git a/pywr-core/src/solvers/builder.rs b/pywr-core/src/solvers/builder.rs index 65228728..0390f2d6 100644 --- a/pywr-core/src/solvers/builder.rs +++ b/pywr-core/src/solvers/builder.rs @@ -4,7 +4,7 @@ use crate::network::Network; use crate::node::{Node, NodeType}; use crate::solvers::col_edge_map::{ColumnEdgeMap, ColumnEdgeMapBuilder}; use crate::solvers::SolverTimings; -use crate::state::State; +use crate::state::{ConstParameterValues, State}; use crate::timestep::Timestep; use crate::PywrError; use std::collections::BTreeMap; @@ -276,12 +276,18 @@ where } } +struct AggNodeFactorRow { + agg_node_idx: AggregatedNodeIndex, + // Row index for each node-pair. If `None` the row is fixed and does not need updating. + row_indices: Vec>, +} + pub struct BuiltSolver { builder: Lp, col_edge_map: ColumnEdgeMap, node_constraints_row_ids: Vec, agg_node_constraint_row_ids: Vec, - agg_node_factor_constraint_row_ids: Vec<(AggregatedNodeIndex, I)>, + agg_node_factor_constraint_row_ids: Vec>, virtual_storage_constraint_row_ids: Vec, } @@ -413,23 +419,33 @@ where } fn update_aggregated_node_factor_constraints(&mut self, network: &Network, state: &State) -> Result<(), PywrError> { - for (agg_node_id, row_id) in self.agg_node_factor_constraint_row_ids.iter() { - let agg_node = network.get_aggregated_node(agg_node_id)?; + // Update the aggregated node factor constraints which are *not* constant + for agg_node_row in self.agg_node_factor_constraint_row_ids.iter() { + let agg_node = network.get_aggregated_node(&agg_node_row.agg_node_idx)?; // Only create row for nodes that have factors if let Some(node_pairs) = agg_node.get_norm_factor_pairs(network, state) { - for node_pair in node_pairs { - // Modify the constraint matrix coefficients for the nodes - // TODO error handling? - let nodes = network.nodes(); - let node0 = nodes.get(&node_pair.node0.index).expect("Node index not found!"); - let node1 = nodes.get(&node_pair.node1.index).expect("Node index not found!"); - - self.builder - .update_row_coefficients(*row_id, node0, 1.0, &self.col_edge_map); - self.builder - .update_row_coefficients(*row_id, node1, -node_pair.ratio(), &self.col_edge_map); - - self.builder.apply_row_bounds(row_id.to_usize().unwrap(), 0.0, 0.0); + assert_eq!( + agg_node_row.row_indices.len(), + node_pairs.len(), + "Row indices and node pairs do not match!" + ); + + for (node_pair, row_idx) in node_pairs.iter().zip(agg_node_row.row_indices.iter()) { + // Only update pairs with a row index (i.e. not fixed) + if let Some(row_idx) = row_idx { + // Modify the constraint matrix coefficients for the nodes + // TODO error handling? + let nodes = network.nodes(); + let node0 = nodes.get(&node_pair.node0.index).expect("Node index not found!"); + let node1 = nodes.get(&node_pair.node1.index).expect("Node index not found!"); + + self.builder + .update_row_coefficients(*row_idx, node0, 1.0, &self.col_edge_map); + self.builder + .update_row_coefficients(*row_idx, node1, -node_pair.ratio(), &self.col_edge_map); + + self.builder.apply_row_bounds(row_idx.to_usize().unwrap(), 0.0, 0.0); + } } } else { panic!("No factor pairs found for an aggregated node that was setup with factors?!"); @@ -504,7 +520,7 @@ where self.col_edge_map.col_for_edge(edge_index) } - pub fn create(mut self, network: &Network) -> Result, PywrError> { + pub fn create(mut self, network: &Network, values: &ConstParameterValues) -> Result, PywrError> { // Create the columns self.create_columns(network)?; @@ -515,7 +531,7 @@ where // Create the aggregated node constraints let agg_node_constraint_row_ids = self.create_aggregated_node_constraints(network); // Create the aggregated node factor constraints - let agg_node_factor_constraint_row_ids = self.create_aggregated_node_factor_constraints(network); + let agg_node_factor_constraint_row_ids = self.create_aggregated_node_factor_constraints(network, values); // Create virtual storage constraints let virtual_storage_constraint_row_ids = self.create_virtual_storage_constraints(network); @@ -671,31 +687,51 @@ where /// Create aggregated node factor constraints /// /// One constraint is created per node to enforce any factor constraints. - fn create_aggregated_node_factor_constraints(&mut self, network: &Network) -> Vec<(AggregatedNodeIndex, I)> { + fn create_aggregated_node_factor_constraints( + &mut self, + network: &Network, + values: &ConstParameterValues, + ) -> Vec> { let mut row_ids = Vec::new(); for agg_node in network.aggregated_nodes().deref() { // Only create row for nodes that have factors - if let Some(node_pairs) = agg_node.get_factor_node_pairs() { - for (n0, n1) in node_pairs { + if let Some(node_pairs) = agg_node.get_const_norm_factor_pairs(values) { + let mut row_indices_for_agg_node = Vec::with_capacity(node_pairs.len()); + + for node_pair in node_pairs { // Create rows for each node in the aggregated node pair with the first one. let mut row = RowBuilder::default(); // TODO error handling? let nodes = network.nodes(); - let node0 = nodes.get(&n0).expect("Node index not found!"); - let node1 = nodes.get(&n1).expect("Node index not found!"); + let node0 = nodes.get(&node_pair.node0.index).expect("Node index not found!"); + let node1 = nodes.get(&node_pair.node1.index).expect("Node index not found!"); + + let ratio = node_pair.ratio(); self.add_node(node0, 1.0, &mut row); - self.add_node(node1, -1.0, &mut row); + self.add_node(node1, -ratio.unwrap_or(1.0), &mut row); // Make the row fixed at zero RHS row.set_lower(0.0); row.set_upper(0.0); - let row_id = self.builder.add_variable_row(row); - row_ids.push((agg_node.index(), row_id)) + // Row is fixed if we can compute the ratio now + if ratio.is_some() { + self.builder.add_fixed_row(row); + row_indices_for_agg_node.push(None) + } else { + // These rows will be updated with the correct ratio later + let row_idx = self.builder.add_variable_row(row); + row_indices_for_agg_node.push(Some(row_idx)); + } } + + row_ids.push(AggNodeFactorRow { + agg_node_idx: agg_node.index(), + row_indices: row_indices_for_agg_node, + }) } } diff --git a/pywr-core/src/solvers/cbc/mod.rs b/pywr-core/src/solvers/cbc/mod.rs index 802f3d05..d88fad76 100644 --- a/pywr-core/src/solvers/cbc/mod.rs +++ b/pywr-core/src/solvers/cbc/mod.rs @@ -4,7 +4,7 @@ use super::builder::SolverBuilder; use crate::network::Network; use crate::solvers::builder::BuiltSolver; use crate::solvers::{Solver, SolverFeatures, SolverTimings}; -use crate::state::State; +use crate::state::{ConstParameterValues, State}; use crate::timestep::Timestep; use crate::PywrError; use coin_or_sys::cbc::*; @@ -217,12 +217,20 @@ impl Solver for CbcSolver { } fn features() -> &'static [SolverFeatures] { - &[SolverFeatures::AggregatedNode, SolverFeatures::VirtualStorage] + &[ + SolverFeatures::AggregatedNode, + SolverFeatures::VirtualStorage, + SolverFeatures::AggregatedNodeFactors, + ] } - fn setup(model: &Network, _settings: &Self::Settings) -> Result, PywrError> { + fn setup( + model: &Network, + values: &ConstParameterValues, + _settings: &Self::Settings, + ) -> Result, PywrError> { let builder = SolverBuilder::default(); - let built = builder.create(model)?; + let built = builder.create(model, values)?; let solver = CbcSolver::from_builder(built); Ok(Box::new(solver)) diff --git a/pywr-core/src/solvers/clp/mod.rs b/pywr-core/src/solvers/clp/mod.rs index 5f659503..e539b731 100644 --- a/pywr-core/src/solvers/clp/mod.rs +++ b/pywr-core/src/solvers/clp/mod.rs @@ -4,7 +4,7 @@ use super::builder::SolverBuilder; use crate::network::Network; use crate::solvers::builder::BuiltSolver; use crate::solvers::{Solver, SolverFeatures, SolverTimings}; -use crate::state::State; +use crate::state::{ConstParameterValues, State}; use crate::timestep::Timestep; use crate::PywrError; use coin_or_sys::clp::*; @@ -237,13 +237,18 @@ impl Solver for ClpSolver { &[ SolverFeatures::AggregatedNode, SolverFeatures::AggregatedNodeFactors, + SolverFeatures::AggregatedNodeDynamicFactors, SolverFeatures::VirtualStorage, ] } - fn setup(model: &Network, _settings: &Self::Settings) -> Result, PywrError> { + fn setup( + model: &Network, + values: &ConstParameterValues, + _settings: &Self::Settings, + ) -> Result, PywrError> { let builder = SolverBuilder::default(); - let built = builder.create(model)?; + let built = builder.create(model, values)?; let solver = ClpSolver::from_builder(built); Ok(Box::new(solver)) diff --git a/pywr-core/src/solvers/highs/mod.rs b/pywr-core/src/solvers/highs/mod.rs index d52a9b9b..3b18af7e 100644 --- a/pywr-core/src/solvers/highs/mod.rs +++ b/pywr-core/src/solvers/highs/mod.rs @@ -3,11 +3,11 @@ mod settings; use crate::network::Network; use crate::solvers::builder::{BuiltSolver, SolverBuilder}; use crate::solvers::{Solver, SolverFeatures, SolverTimings}; -use crate::state::State; +use crate::state::{ConstParameterValues, State}; use crate::timestep::Timestep; use crate::PywrError; use highs_sys::{ - HighsInt, Highs_addCols, Highs_addRows, Highs_changeColsCostByRange, Highs_changeObjectiveSense, + HighsInt, Highs_addCols, Highs_addRows, Highs_changeCoeff, Highs_changeColsCostByRange, Highs_changeObjectiveSense, Highs_changeRowsBoundsByMask, Highs_create, Highs_getDoubleInfoValue, Highs_getSolution, Highs_run, Highs_setBoolOptionValue, Highs_setStringOptionValue, OBJECTIVE_SENSE_MINIMIZE, STATUS_OK, }; @@ -112,6 +112,13 @@ impl Highs { } } + pub fn change_coeff(&mut self, row: HighsInt, col: HighsInt, value: f64) { + unsafe { + let ret = Highs_changeCoeff(self.ptr, row, col, value); + assert_eq!(ret, STATUS_OK); + } + } + pub fn run(&mut self) { unsafe { let status = Highs_run(self.ptr); @@ -167,12 +174,20 @@ impl Solver for HighsSolver { } fn features() -> &'static [SolverFeatures] { - &[] + &[ + SolverFeatures::AggregatedNode, + SolverFeatures::AggregatedNodeFactors, + SolverFeatures::AggregatedNodeDynamicFactors, + ] } - fn setup(network: &Network, settings: &Self::Settings) -> Result, PywrError> { + fn setup( + network: &Network, + values: &ConstParameterValues, + settings: &Self::Settings, + ) -> Result, PywrError> { let builder: SolverBuilder = SolverBuilder::default(); - let built = builder.create(network)?; + let built = builder.create(network, values)?; let num_cols = built.num_cols(); let num_rows = built.num_rows(); @@ -215,6 +230,11 @@ impl Solver for HighsSolver { self.builder.row_lower(), self.builder.row_upper(), ); + + for (row, column, coefficient) in self.builder.coefficients_to_update() { + self.highs.change_coeff(*row, *column, *coefficient); + } + timings.update_constraints += now.elapsed(); let now = Instant::now(); diff --git a/pywr-core/src/solvers/mod.rs b/pywr-core/src/solvers/mod.rs index 11d6ccaa..c5974674 100644 --- a/pywr-core/src/solvers/mod.rs +++ b/pywr-core/src/solvers/mod.rs @@ -1,5 +1,5 @@ use crate::network::Network; -use crate::state::State; +use crate::state::{ConstParameterValues, State}; use crate::timestep::Timestep; use crate::PywrError; use std::ops::{Add, AddAssign}; @@ -72,6 +72,7 @@ impl AddAssign for SolverTimings { pub enum SolverFeatures { AggregatedNode, AggregatedNodeFactors, + AggregatedNodeDynamicFactors, VirtualStorage, } @@ -87,7 +88,8 @@ pub trait Solver: Send { fn name() -> &'static str; /// An array of features that this solver provides. fn features() -> &'static [SolverFeatures]; - fn setup(model: &Network, settings: &Self::Settings) -> Result, PywrError>; + fn setup(model: &Network, values: &ConstParameterValues, settings: &Self::Settings) + -> Result, PywrError>; fn solve(&mut self, model: &Network, timestep: &Timestep, state: &mut State) -> Result; } diff --git a/pywr-core/src/test_utils.rs b/pywr-core/src/test_utils.rs index 0a9f9666..0fe1de6c 100644 --- a/pywr-core/src/test_utils.rs +++ b/pywr-core/src/test_utils.rs @@ -171,13 +171,20 @@ pub fn run_and_assert_parameter( /// The model will only be run if the solver has the required solver features (and /// is also enabled as a Cargo feature). pub fn run_all_solvers(model: &Model, solvers_without_features: &[&str]) { + println!("Running CLP"); check_features_and_run::(model, !solvers_without_features.contains(&"clp")); #[cfg(feature = "cbc")] - check_features_and_run::(model, !solvers_without_features.contains(&"cbc")); + { + println!("Running CBC"); + check_features_and_run::(model, !solvers_without_features.contains(&"cbc")); + } #[cfg(feature = "highs")] - check_features_and_run::(model, !solvers_without_features.contains(&"highs")); + { + println!("Running Highs"); + check_features_and_run::(model, !solvers_without_features.contains(&"highs")); + } #[cfg(feature = "ipm-simd")] { diff --git a/pywr-schema/src/nodes/water_treatment_works.rs b/pywr-schema/src/nodes/water_treatment_works.rs index c65b42a5..60536cda 100644 --- a/pywr-schema/src/nodes/water_treatment_works.rs +++ b/pywr-schema/src/nodes/water_treatment_works.rs @@ -417,6 +417,6 @@ mod tests { network.add_recorder(Box::new(recorder)).unwrap(); // Test all solvers - run_all_solvers(&model, &["cbc", "highs"]); + run_all_solvers(&model, &[]); } } From db125a869261e19e93efc7ddf646b2c6a9ef4004 Mon Sep 17 00:00:00 2001 From: James Tomlinson Date: Mon, 22 Jul 2024 15:41:23 +0100 Subject: [PATCH 2/3] chore: Address Clippy warnings. --- pywr-core/src/metric.rs | 5 +---- pywr-schema/src/model.rs | 2 +- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/pywr-core/src/metric.rs b/pywr-core/src/metric.rs index 10bda0b6..79aa80b5 100644 --- a/pywr-core/src/metric.rs +++ b/pywr-core/src/metric.rs @@ -56,10 +56,7 @@ impl SimpleMetricF64 { /// Returns true if the metric is a constant value. pub fn is_constant(&self) -> bool { - match self { - SimpleMetricF64::Constant(_) => true, - _ => false, - } + matches!(self, SimpleMetricF64::Constant(_)) } } diff --git a/pywr-schema/src/model.rs b/pywr-schema/src/model.rs index 84fa6b0c..c2a17556 100644 --- a/pywr-schema/src/model.rs +++ b/pywr-schema/src/model.rs @@ -879,7 +879,7 @@ impl PywrMultiNetworkModel { let mut model = pywr_core::models::MultiNetworkModel::new(domain); for (name, network) in networks { - model.add_network(&name, network); + model.add_network(&name, network)?; } for (from_network_idx, from_metric, to_network_idx, initial_value) in inter_network_transfers { From ce4301cadc644be880734d96e8adda7c18f83948 Mon Sep 17 00:00:00 2001 From: James Tomlinson Date: Thu, 25 Jul 2024 16:54:15 +0100 Subject: [PATCH 3/3] fix: Add checks for negative factors in aggregated node. --- pywr-core/src/aggregated_node.rs | 66 +++++++++++++++++++++++++------- pywr-core/src/lib.rs | 2 + 2 files changed, 54 insertions(+), 14 deletions(-) diff --git a/pywr-core/src/aggregated_node.rs b/pywr-core/src/aggregated_node.rs index 37e68b3a..92cec115 100644 --- a/pywr-core/src/aggregated_node.rs +++ b/pywr-core/src/aggregated_node.rs @@ -313,11 +313,17 @@ fn get_norm_proportional_factor_pairs( // First get the current factor values let values: Vec = factors .iter() - .map(|f| f.get_value(model, state)) + .map(|f| { + let v = f.get_value(model, state)?; + if v < 0.0 { + Err(PywrError::NegativeFactor) + } else { + Ok(v) + } + }) .collect::, PywrError>>() - .expect("Failed to get current factor values."); + .expect("Failed to get current factor values. Ensure that all factors are not negative."); - // TODO do we need to assert that each individual factor is positive? let total: f64 = values.iter().sum(); if total < 0.0 { panic!("Proportional factors are too small or negative."); @@ -352,12 +358,23 @@ fn get_const_norm_proportional_factor_pairs( panic!("Found {} proportional factors and {} nodes in aggregated node. The number of proportional factors should equal one less than the number of nodes.", factors.len(), nodes.len()); } - // First get the current factor values + // First get the current factor values, ensuring they are all non-negative let values: Vec> = factors .iter() - .map(|f| f.try_get_constant_value(values)) + .map(|f| { + let v = f.try_get_constant_value(values)?; + if let Some(v) = v { + if v < 0.0 { + Err(PywrError::NegativeFactor) + } else { + Ok(Some(v)) + } + } else { + Ok(None) + } + }) .collect::, PywrError>>() - .expect("Failed to get current factor values."); + .expect("Failed to get current factor values. Ensure that all factors are not negative."); let n0 = nodes[0]; @@ -374,7 +391,6 @@ fn get_const_norm_proportional_factor_pairs( .collect::>() } else { // All factors are available; therefore we can calculate "f0" - // TODO do we need to assert that each individual factor is positive? let total: f64 = values .iter() .map(|v| v.expect("Factor is `None`; this should be impossible.")) @@ -414,18 +430,24 @@ fn get_norm_ratio_factor_pairs( let n0 = nodes[0]; let f0 = factors[0].get_value(model, state).unwrap(); + if f0 < 0.0 { + panic!("Negative factor is not allowed"); + } nodes .iter() .zip(factors) .skip(1) .map(move |(&n1, f1)| { - NodeFactorPair::new( - NodeFactor::new(n0, f0), - NodeFactor::new(n1, f1.get_value(model, state).unwrap()), - ) + let v1 = f1.get_value(model, state)?; + if v1 < 0.0 { + Err(PywrError::NegativeFactor) + } else { + Ok(NodeFactorPair::new(NodeFactor::new(n0, f0), NodeFactor::new(n1, v1))) + } }) - .collect::>() + .collect::, PywrError>>() + .expect("Failed to get current factor values. Ensure that all factors are not negative.") } /// Constant ratio factors using constant values if they are available. If they are not available, @@ -446,6 +468,12 @@ fn get_const_norm_ratio_factor_pairs( .try_get_constant_value(values) .unwrap_or_else(|e| panic!("Failed to get constant value for factor: {}", e)); + if let Some(v0) = f0 { + if v0 < 0.0 { + panic!("Negative factor is not allowed"); + } + } + nodes .iter() .zip(factors) @@ -455,9 +483,19 @@ fn get_const_norm_ratio_factor_pairs( .try_get_constant_value(values) .unwrap_or_else(|e| panic!("Failed to get constant value for factor: {}", e)); - NodeConstFactorPair::new(NodeConstFactor::new(n0, f0), NodeConstFactor::new(n1, v1)) + if let Some(v) = v1 { + if v < 0.0 { + return Err(PywrError::NegativeFactor); + } + } + + Ok(NodeConstFactorPair::new( + NodeConstFactor::new(n0, f0), + NodeConstFactor::new(n1, v1), + )) }) - .collect::>() + .collect::, PywrError>>() + .expect("Failed to get current factor values. Ensure that all factors are not negative.") } #[cfg(test)] diff --git a/pywr-core/src/lib.rs b/pywr-core/src/lib.rs index 7220c577..381786c5 100644 --- a/pywr-core/src/lib.rs +++ b/pywr-core/src/lib.rs @@ -197,6 +197,8 @@ pub enum PywrError { Aggregation(#[from] AggregationError), #[error("cannot simplify metric")] CannotSimplifyMetric, + #[error("Negative factor is not allowed")] + NegativeFactor, } // Python errors