From b4771fe5a464047c75fc3de95e3c3584655a4bea Mon Sep 17 00:00:00 2001 From: James Tomlinson Date: Thu, 9 Nov 2023 21:30:50 +0000 Subject: [PATCH] feat: Initial commit of DerivedMetric (#68) A derived metric is updated after solve and is for values that could be dependent on state AND paramter values. These need updating and storing for use between time-steps. Adds a test for piecewise storage model that shows how the volume is updated at the end of time-step, and the derived proportional volume is calculated and retained for use in the next time-step despite the max volumes being updated. Fixes #63. --- pywr-core/src/derived_metric.rs | 96 +++++++++++++++ pywr-core/src/lib.rs | 8 ++ pywr-core/src/metric.rs | 89 +++----------- pywr-core/src/model.rs | 110 +++++++++++++++++- pywr-core/src/node.rs | 2 +- .../src/parameters/control_curves/mod.rs | 2 + .../control_curves/volume_between.rs | 58 +++++++++ pywr-core/src/parameters/mod.rs | 2 +- pywr-core/src/parameters/py.rs | 4 +- pywr-core/src/parameters/rhai.rs | 2 +- pywr-core/src/recorders/csv.rs | 28 +---- pywr-core/src/recorders/hdf.rs | 19 +-- pywr-core/src/recorders/mod.rs | 58 ++++++++- pywr-core/src/state.rs | 21 ++++ pywr-schema/src/nodes/piecewise_storage.rs | 55 +++++++-- pywr-schema/src/parameters/mod.rs | 4 +- pywr-schema/src/parameters/polynomial.rs | 9 +- .../src/test_models/piecewise_storage2.json | 18 +++ 18 files changed, 442 insertions(+), 143 deletions(-) create mode 100644 pywr-core/src/derived_metric.rs create mode 100644 pywr-core/src/parameters/control_curves/volume_between.rs diff --git a/pywr-core/src/derived_metric.rs b/pywr-core/src/derived_metric.rs new file mode 100644 index 00000000..13bd1763 --- /dev/null +++ b/pywr-core/src/derived_metric.rs @@ -0,0 +1,96 @@ +use crate::aggregated_storage_node::AggregatedStorageNodeIndex; +use crate::model::Model; +use crate::node::NodeIndex; +use crate::state::State; +use crate::timestep::Timestep; +use crate::virtual_storage::VirtualStorageIndex; +use crate::PywrError; +use std::fmt; +use std::fmt::{Display, Formatter}; +use std::ops::Deref; + +#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Debug)] +pub struct DerivedMetricIndex(usize); + +impl Deref for DerivedMetricIndex { + type Target = usize; + + fn deref(&self) -> &Self::Target { + &self.0 + } +} + +impl DerivedMetricIndex { + pub fn new(idx: usize) -> Self { + Self(idx) + } +} + +impl Display for DerivedMetricIndex { + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { + write!(f, "{}", self.0) + } +} + +/// Derived metrics are updated after the model is solved. +/// +/// These metrics are "derived" from node states (e.g. volume, flow) and must be updated +/// after those states have been updated. This should happen after the model is solved. The values +/// are then available in this state for the next time-step. +#[derive(Clone, Debug, PartialEq)] +pub enum DerivedMetric { + NodeInFlowDeficit(NodeIndex), + NodeProportionalVolume(NodeIndex), + AggregatedNodeProportionalVolume(AggregatedStorageNodeIndex), + VirtualStorageProportionalVolume(VirtualStorageIndex), +} + +impl DerivedMetric { + pub fn before(&self, timestep: &Timestep, model: &Model, state: &State) -> Result, PywrError> { + // On the first time-step set the initial value + if timestep.is_first() { + self.compute(model, state).map(|v| Some(v)) + } else { + Ok(None) + } + } + + pub fn compute(&self, model: &Model, state: &State) -> Result { + match self { + Self::NodeProportionalVolume(idx) => { + let max_volume = model.get_node(idx)?.get_current_max_volume(model, state)?; + Ok(state + .get_network_state() + .get_node_proportional_volume(idx, max_volume)?) + } + Self::VirtualStorageProportionalVolume(idx) => { + let max_volume = model.get_virtual_storage_node(idx)?.get_max_volume(model, state)?; + Ok(state + .get_network_state() + .get_virtual_storage_proportional_volume(idx, max_volume)?) + } + Self::AggregatedNodeProportionalVolume(idx) => { + let node = model.get_aggregated_storage_node(idx)?; + let volume: f64 = node + .nodes + .iter() + .map(|idx| state.get_network_state().get_node_volume(idx)) + .sum::>()?; + + let max_volume: f64 = node + .nodes + .iter() + .map(|idx| model.get_node(idx)?.get_current_max_volume(model, state)) + .sum::>()?; + // TODO handle divide by zero + Ok(volume / max_volume) + } + Self::NodeInFlowDeficit(idx) => { + let node = model.get_node(idx)?; + let flow = state.get_network_state().get_node_in_flow(idx)?; + let max_flow = node.get_current_max_flow(model, state)?; + Ok(max_flow - flow) + } + } + } +} diff --git a/pywr-core/src/lib.rs b/pywr-core/src/lib.rs index 9c82d52b..1eafe82b 100644 --- a/pywr-core/src/lib.rs +++ b/pywr-core/src/lib.rs @@ -2,6 +2,7 @@ extern crate core; +use crate::derived_metric::DerivedMetricIndex; use crate::node::NodeIndex; use crate::parameters::{IndexParameterIndex, MultiValueParameterIndex, ParameterIndex}; use crate::recorders::RecorderIndex; @@ -11,6 +12,7 @@ use thiserror::Error; pub mod aggregated_node; mod aggregated_storage_node; +pub mod derived_metric; pub mod edge; pub mod metric; pub mod model; @@ -55,6 +57,10 @@ pub enum PywrError { RecorderIndexNotFound, #[error("recorder not found")] RecorderNotFound, + #[error("derived metric not found")] + DerivedMetricNotFound, + #[error("derived metric index {0} not found")] + DerivedMetricIndexNotFound(DerivedMetricIndex), #[error("node name `{0}` already exists")] NodeNameAlreadyExists(String), #[error("parameter name `{0}` already exists at index {1}")] @@ -127,6 +133,8 @@ pub enum PywrError { ParameterVariableValuesIncorrectLength, #[error("missing solver features")] MissingSolverFeatures, + #[error("parameters do not provide an initial value")] + ParameterNoInitialValue, } // Python errors diff --git a/pywr-core/src/metric.rs b/pywr-core/src/metric.rs index 21452aac..6431bff3 100644 --- a/pywr-core/src/metric.rs +++ b/pywr-core/src/metric.rs @@ -1,47 +1,25 @@ use crate::aggregated_node::AggregatedNodeIndex; use crate::aggregated_storage_node::AggregatedStorageNodeIndex; +use crate::derived_metric::DerivedMetricIndex; use crate::edge::EdgeIndex; use crate::model::Model; use crate::node::NodeIndex; -use crate::parameters::{MultiValueParameterIndex, ParameterIndex}; +use crate::parameters::{IndexParameterIndex, MultiValueParameterIndex, ParameterIndex}; use crate::state::State; use crate::virtual_storage::VirtualStorageIndex; use crate::PywrError; - -#[derive(Clone, Debug, PartialEq)] -pub struct VolumeBetweenControlCurves { - max_volume: Box, - upper: Option>, - lower: Option>, -} - -impl VolumeBetweenControlCurves { - pub fn new(max_volume: Metric, upper: Option, lower: Option) -> Self { - Self { - max_volume: Box::new(max_volume), - upper: upper.map(Box::new), - lower: lower.map(Box::new), - } - } -} - #[derive(Clone, Debug, PartialEq)] pub enum Metric { NodeInFlow(NodeIndex), NodeOutFlow(NodeIndex), NodeVolume(NodeIndex), - NodeInFlowDeficit(NodeIndex), - NodeProportionalVolume(NodeIndex), AggregatedNodeInFlow(AggregatedNodeIndex), AggregatedNodeOutFlow(AggregatedNodeIndex), AggregatedNodeVolume(AggregatedStorageNodeIndex), - AggregatedNodeProportionalVolume(AggregatedStorageNodeIndex), EdgeFlow(EdgeIndex), ParameterValue(ParameterIndex), MultiParameterValue((MultiValueParameterIndex, String)), VirtualStorageVolume(VirtualStorageIndex), - VirtualStorageProportionalVolume(VirtualStorageIndex), - VolumeBetweenControlCurves(VolumeBetweenControlCurves), MultiNodeInFlow { indices: Vec, name: String, @@ -49,6 +27,7 @@ pub enum Metric { }, // TODO implement other MultiNodeXXX variants Constant(f64), + DerivedMetric(DerivedMetricIndex), } impl Metric { @@ -73,22 +52,12 @@ impl Metric { .map(|idx| state.get_network_state().get_node_out_flow(idx)) .sum::>() } - Metric::NodeProportionalVolume(idx) => { - let max_volume = model.get_node(idx)?.get_current_max_volume(model, state)?; - Ok(state - .get_network_state() - .get_node_proportional_volume(idx, max_volume)?) - } + Metric::EdgeFlow(idx) => Ok(state.get_network_state().get_edge_flow(idx)?), Metric::ParameterValue(idx) => Ok(state.get_parameter_value(*idx)?), Metric::MultiParameterValue((idx, key)) => Ok(state.get_multi_parameter_value(*idx, key)?), Metric::VirtualStorageVolume(idx) => Ok(state.get_network_state().get_virtual_storage_volume(idx)?), - Metric::VirtualStorageProportionalVolume(idx) => { - let max_volume = model.get_virtual_storage_node(idx)?.get_max_volume(model, state)?; - Ok(state - .get_network_state() - .get_virtual_storage_proportional_volume(idx, max_volume)?) - } + Metric::DerivedMetric(idx) => state.get_derived_metric_value(*idx), Metric::Constant(v) => Ok(*v), Metric::AggregatedNodeVolume(idx) => { let node = model.get_aggregated_storage_node(idx)?; @@ -97,22 +66,7 @@ impl Metric { .map(|idx| state.get_network_state().get_node_volume(idx)) .sum::>() } - Metric::AggregatedNodeProportionalVolume(idx) => { - let node = model.get_aggregated_storage_node(idx)?; - let volume: f64 = node - .nodes - .iter() - .map(|idx| state.get_network_state().get_node_volume(idx)) - .sum::>()?; - let max_volume: f64 = node - .nodes - .iter() - .map(|idx| model.get_node(idx)?.get_current_max_volume(model, state)) - .sum::>()?; - // TODO handle divide by zero - Ok(volume / max_volume) - } Metric::MultiNodeInFlow { indices, .. } => { let flow = indices .iter() @@ -120,26 +74,21 @@ impl Metric { .sum::>()?; Ok(flow) } - Metric::NodeInFlowDeficit(idx) => { - let node = model.get_node(idx)?; - let flow = state.get_network_state().get_node_in_flow(idx)?; - let max_flow = node.get_current_max_flow(model, state)?; - Ok(max_flow - flow) - } - Metric::VolumeBetweenControlCurves(vol) => { - let max_volume = vol.max_volume.get_value(model, state)?; - let lower = vol - .lower - .as_ref() - .map_or(Ok(0.0), |metric| metric.get_value(model, state))?; - let upper = vol - .upper - .as_ref() - .map_or(Ok(1.0), |metric| metric.get_value(model, state))?; + } + } +} - // TODO handle invalid bounds - Ok(max_volume * (upper - lower)) - } +#[derive(Clone, Debug, PartialEq)] +pub enum IndexMetric { + IndexParameterValue(IndexParameterIndex), + Constant(usize), +} + +impl IndexMetric { + pub fn get_value(&self, model: &Model, state: &State) -> Result { + match self { + Self::IndexParameterValue(idx) => state.get_parameter_index(*idx), + Self::Constant(i) => Ok(*i), } } } diff --git a/pywr-core/src/model.rs b/pywr-core/src/model.rs index abf363d1..da236409 100644 --- a/pywr-core/src/model.rs +++ b/pywr-core/src/model.rs @@ -1,5 +1,6 @@ use crate::aggregated_node::{AggregatedNode, AggregatedNodeIndex, AggregatedNodeVec, Factors}; use crate::aggregated_storage_node::{AggregatedStorageNode, AggregatedStorageNodeIndex, AggregatedStorageNodeVec}; +use crate::derived_metric::{DerivedMetric, DerivedMetricIndex}; use crate::edge::{EdgeIndex, EdgeVec}; use crate::metric::Metric; use crate::node::{ConstraintValue, Node, NodeVec, StorageInitialVolume}; @@ -138,6 +139,7 @@ enum ComponentType { Node(NodeIndex), VirtualStorageNode(VirtualStorageIndex), Parameter(ParameterType), + DerivedMetric(DerivedMetricIndex), } #[derive(Default)] @@ -151,6 +153,7 @@ pub struct Model { parameters: Vec>, index_parameters: Vec>, multi_parameters: Vec>, + derived_metrics: Vec, metric_sets: Vec, resolve_order: Vec, recorders: Vec>, @@ -206,6 +209,7 @@ impl Model { initial_values_states.len(), initial_indices_states.len(), initial_multi_param_states.len(), + self.derived_metrics.len(), ); states.push(state); @@ -748,6 +752,16 @@ impl Model { } } } + ComponentType::DerivedMetric(idx) => { + // Compute derived metrics in before + let m = self + .derived_metrics + .get(*idx.deref()) + .ok_or(PywrError::DerivedMetricIndexNotFound(*idx))?; + if let Some(value) = m.before(timestep, self, state)? { + state.set_derived_metric_value(*idx, value)?; + } + } } } @@ -821,6 +835,15 @@ impl Model { } } } + ComponentType::DerivedMetric(idx) => { + // Compute derived metrics in "after" + let m = self + .derived_metrics + .get(*idx.deref()) + .ok_or(PywrError::DerivedMetricIndexNotFound(*idx))?; + let value = m.compute(self, state)?; + state.set_derived_metric_value(*idx, value)?; + } } } @@ -1089,7 +1112,7 @@ impl Model { } pub fn get_storage_node_metric( - &self, + &mut self, name: &str, sub_name: Option<&str>, proportional: bool, @@ -1097,19 +1120,25 @@ impl Model { if let Ok(idx) = self.get_node_index_by_name(name, sub_name) { // A regular node if proportional { - Ok(Metric::NodeProportionalVolume(idx)) + // Proportional is a derived metric + let dm_idx = self.add_derived_metric(DerivedMetric::NodeProportionalVolume(idx)); + Ok(Metric::DerivedMetric(dm_idx)) } else { Ok(Metric::NodeVolume(idx)) } } else if let Ok(idx) = self.get_aggregated_storage_node_index_by_name(name, sub_name) { if proportional { - Ok(Metric::AggregatedNodeProportionalVolume(idx)) + // Proportional is a derived metric + let dm_idx = self.add_derived_metric(DerivedMetric::AggregatedNodeProportionalVolume(idx)); + Ok(Metric::DerivedMetric(dm_idx)) } else { Ok(Metric::AggregatedNodeVolume(idx)) } } else if let Ok(node) = self.get_virtual_storage_node_by_name(name, sub_name) { if proportional { - Ok(Metric::VirtualStorageProportionalVolume(node.index())) + // Proportional is a derived metric + let dm_idx = self.add_derived_metric(DerivedMetric::VirtualStorageProportionalVolume(node.index())); + Ok(Metric::DerivedMetric(dm_idx)) } else { Ok(Metric::VirtualStorageVolume(node.index())) } @@ -1118,6 +1147,36 @@ impl Model { } } + /// Get a [`DerivedMetricIndex`] for the given derived metric + pub fn get_derived_metric_index(&self, derived_metric: &DerivedMetric) -> Result { + let idx = self + .derived_metrics + .iter() + .position(|dm| dm == derived_metric) + .ok_or(PywrError::DerivedMetricNotFound)?; + + Ok(DerivedMetricIndex::new(idx)) + } + + /// Get a [`DerivedMetricIndex`] for the given derived metric + pub fn get_derived_metric(&self, index: &DerivedMetricIndex) -> Result<&DerivedMetric, PywrError> { + self.derived_metrics + .get(*index.deref()) + .ok_or(PywrError::DerivedMetricNotFound) + } + + pub fn add_derived_metric(&mut self, derived_metric: DerivedMetric) -> DerivedMetricIndex { + match self.get_derived_metric_index(&derived_metric) { + Ok(idx) => idx, + Err(_) => { + self.derived_metrics.push(derived_metric); + let idx = DerivedMetricIndex::new(self.derived_metrics.len() - 1); + self.resolve_order.push(ComponentType::DerivedMetric(idx)); + idx + } + } + } + /// Get a `Parameter` from a parameter's name pub fn get_parameter(&self, index: &ParameterIndex) -> Result<&dyn parameters::Parameter, PywrError> { match self.parameters.get(*index.deref()) { @@ -1476,7 +1535,7 @@ mod tests { use crate::metric::Metric; use crate::model::Model; use crate::node::{Constraint, ConstraintValue}; - use crate::parameters::{ActivationFunction, Parameter, VariableParameter}; + use crate::parameters::{ActivationFunction, InterpolatedParameter, Parameter, VariableParameter}; use crate::recorders::AssertionRecorder; use crate::scenario::{ScenarioGroupCollection, ScenarioIndex}; #[cfg(feature = "clipm")] @@ -1674,6 +1733,47 @@ mod tests { run_all_solvers(&model, ×tepper); } + /// Test proportional storage derived metric. + /// + /// Proportional storage is a derived metric that is updated after each solve. However, a + /// parameter may required a value for the initial time-step based on the initial volume. + #[test] + fn test_storage_proportional_volume() { + let mut model = simple_storage_model(); + let timestepper = default_timestepper(); + + let idx = model.get_node_by_name("reservoir", None).unwrap().index(); + let dm_idx = model.add_derived_metric(DerivedMetric::NodeProportionalVolume(idx)); + + // These are the expected values for the proportional volume at the end of the time-step + let expected = Array2::from_shape_fn((15, 10), |(i, _j)| (90.0 - 10.0 * i as f64).max(0.0) / 100.0); + let recorder = AssertionRecorder::new( + "reservoir-proportion-volume", + Metric::DerivedMetric(dm_idx), + expected, + None, + None, + ); + model.add_recorder(Box::new(recorder)).unwrap(); + + // Set-up a control curve that uses the proportional volume + // This should be use the initial proportion (100%) on the first time-step, and then the previous day's end value + let cc = InterpolatedParameter::new( + "interp", + Metric::DerivedMetric(dm_idx), + vec![], + vec![Metric::Constant(100.0), Metric::Constant(0.0)], + ); + let p_idx = model.add_parameter(Box::new(cc)).unwrap(); + let expected = Array2::from_shape_fn((15, 10), |(i, _j)| (100.0 - 10.0 * i as f64).max(0.0)); + + let recorder = AssertionRecorder::new("reservoir-cc", Metric::ParameterValue(p_idx), expected, None, None); + model.add_recorder(Box::new(recorder)).unwrap(); + + // Test all solvers + run_all_solvers(&model, ×tepper); + } + #[test] /// Test `ScenarioGroupCollection` iteration fn test_scenario_iteration() { diff --git a/pywr-core/src/node.rs b/pywr-core/src/node.rs index ad9f3d92..17dbeac1 100644 --- a/pywr-core/src/node.rs +++ b/pywr-core/src/node.rs @@ -616,7 +616,7 @@ impl StorageConstraints { ConstraintValue::Metric(m) => m.get_value(model, state), } } - /// Return the current maximum volume from the parameter state + /// Return the current maximum volume from the metric state /// /// Defaults to f64::MAX if no parameter is defined. pub fn get_max_volume(&self, model: &Model, state: &State) -> Result { diff --git a/pywr-core/src/parameters/control_curves/mod.rs b/pywr-core/src/parameters/control_curves/mod.rs index 3f566f22..af4778b2 100644 --- a/pywr-core/src/parameters/control_curves/mod.rs +++ b/pywr-core/src/parameters/control_curves/mod.rs @@ -3,12 +3,14 @@ mod index; mod interpolated; mod piecewise; mod simple; +mod volume_between; pub use apportion::ApportionParameter; pub use index::ControlCurveIndexParameter; pub use interpolated::InterpolatedParameter; pub use piecewise::PiecewiseInterpolatedParameter; pub use simple::ControlCurveParameter; +pub use volume_between::VolumeBetweenControlCurvesParameter; /// Interpolate fn interpolate(value: f64, lower_bound: f64, upper_bound: f64, lower_value: f64, upper_value: f64) -> f64 { diff --git a/pywr-core/src/parameters/control_curves/volume_between.rs b/pywr-core/src/parameters/control_curves/volume_between.rs new file mode 100644 index 00000000..79907eec --- /dev/null +++ b/pywr-core/src/parameters/control_curves/volume_between.rs @@ -0,0 +1,58 @@ +use crate::metric::Metric; +use crate::model::Model; +use crate::parameters::{Parameter, ParameterMeta}; +use crate::scenario::ScenarioIndex; +use crate::state::State; +use crate::timestep::Timestep; +use crate::PywrError; +use std::any::Any; + +/// A parameter that returns the volume that is the proportion between two control curves +pub struct VolumeBetweenControlCurvesParameter { + meta: ParameterMeta, + total: Metric, + upper: Option, + lower: Option, +} + +impl VolumeBetweenControlCurvesParameter { + pub fn new(name: &str, total: Metric, upper: Option, lower: Option) -> Self { + Self { + meta: ParameterMeta::new(name), + total, + upper, + lower, + } + } +} + +impl Parameter for VolumeBetweenControlCurvesParameter { + fn as_any_mut(&mut self) -> &mut dyn Any { + self + } + + fn meta(&self) -> &ParameterMeta { + &self.meta + } + fn compute( + &self, + _timestep: &Timestep, + _scenario_index: &ScenarioIndex, + model: &Model, + state: &State, + _internal_state: &mut Option>, + ) -> Result { + let total = self.total.get_value(model, state)?; + + let lower = self + .lower + .as_ref() + .map_or(Ok(0.0), |metric| metric.get_value(model, state))?; + let upper = self + .upper + .as_ref() + .map_or(Ok(1.0), |metric| metric.get_value(model, state))?; + + Ok(total * (upper - lower)) + } +} diff --git a/pywr-core/src/parameters/mod.rs b/pywr-core/src/parameters/mod.rs index 8c0a2ef5..76d1986c 100644 --- a/pywr-core/src/parameters/mod.rs +++ b/pywr-core/src/parameters/mod.rs @@ -36,7 +36,7 @@ pub use asymmetric::AsymmetricSwitchIndexParameter; pub use constant::ConstantParameter; pub use control_curves::{ ApportionParameter, ControlCurveIndexParameter, ControlCurveParameter, InterpolatedParameter, - PiecewiseInterpolatedParameter, + PiecewiseInterpolatedParameter, VolumeBetweenControlCurvesParameter, }; pub use delay::DelayParameter; pub use division::DivisionParameter; diff --git a/pywr-core/src/parameters/py.rs b/pywr-core/src/parameters/py.rs index 2ade76e8..3b1ddba3 100644 --- a/pywr-core/src/parameters/py.rs +++ b/pywr-core/src/parameters/py.rs @@ -355,7 +355,7 @@ class MyParameter: }, ]; - let state = State::new(vec![], 0, vec![], 1, 0, 0); + let state = State::new(vec![], 0, vec![], 1, 0, 0, 0); let mut internal_p_states: Vec<_> = scenario_indices .iter() @@ -423,7 +423,7 @@ class MyParameter: }, ]; - let state = State::new(vec![], 0, vec![], 1, 0, 0); + let state = State::new(vec![], 0, vec![], 1, 0, 0, 0); let mut internal_p_states: Vec<_> = scenario_indices .iter() diff --git a/pywr-core/src/parameters/rhai.rs b/pywr-core/src/parameters/rhai.rs index 5010f11c..e88e6adc 100644 --- a/pywr-core/src/parameters/rhai.rs +++ b/pywr-core/src/parameters/rhai.rs @@ -161,7 +161,7 @@ mod tests { }, ]; - let state = State::new(vec![], 0, vec![], 1, 0, 0); + let state = State::new(vec![], 0, vec![], 1, 0, 0, 0); let mut internal_p_states: Vec<_> = scenario_indices .iter() diff --git a/pywr-core/src/recorders/csv.rs b/pywr-core/src/recorders/csv.rs index 17a9dabc..ea3f6223 100644 --- a/pywr-core/src/recorders/csv.rs +++ b/pywr-core/src/recorders/csv.rs @@ -77,12 +77,8 @@ impl Recorder for CSVRecorder { (name.to_string(), sub_name, "volume".to_string()) } - Metric::NodeProportionalVolume(idx) => { - let node = model.get_node(idx)?; - let (name, sub_name) = node.full_name(); - let sub_name = sub_name.map_or("".to_string(), |sn| sn.to_string()); - - (name.to_string(), sub_name, "proportional-volume".to_string()) + Metric::DerivedMetric(idx) => { + todo!("Derived metrics are not yet supported in CSV recorders"); } Metric::AggregatedNodeVolume(idx) => { let node = model.get_aggregated_storage_node(idx)?; @@ -91,13 +87,6 @@ impl Recorder for CSVRecorder { (name.to_string(), sub_name, "volume".to_string()) } - Metric::AggregatedNodeProportionalVolume(idx) => { - let node = model.get_aggregated_storage_node(idx)?; - let (name, sub_name) = node.full_name(); - let sub_name = sub_name.map_or("".to_string(), |sn| sn.to_string()); - - (name.to_string(), sub_name, "proportional-volume".to_string()) - } Metric::EdgeFlow(_) => { continue; // TODO } @@ -109,9 +98,6 @@ impl Recorder for CSVRecorder { Metric::VirtualStorageVolume(_) => { continue; // TODO } - Metric::VirtualStorageProportionalVolume(_) => { - continue; // TODO - } Metric::Constant(_) => { continue; // TODO } @@ -132,16 +118,6 @@ impl Recorder for CSVRecorder { (name.to_string(), sub_name, "outflow".to_string()) } - Metric::NodeInFlowDeficit(idx) => { - let node = model.get_node(idx)?; - let (name, sub_name) = node.full_name(); - let sub_name = sub_name.map_or("".to_string(), |sn| sn.to_string()); - - (name.to_string(), sub_name, "inflow-deficit".to_string()) - } - Metric::VolumeBetweenControlCurves(_) => { - todo!("Recording VolumeBetweenControlCurves not implemented.") - } Metric::MultiNodeInFlow { name, sub_name, .. } => ( name.to_string(), sub_name.clone().unwrap_or("".to_string()), diff --git a/pywr-core/src/recorders/hdf.rs b/pywr-core/src/recorders/hdf.rs index 8f859eb6..9bd0f148 100644 --- a/pywr-core/src/recorders/hdf.rs +++ b/pywr-core/src/recorders/hdf.rs @@ -95,18 +95,13 @@ impl Recorder for HDF5Recorder { let node = model.get_node(idx)?; require_node_dataset(root_grp, shape, node.name(), node.sub_name(), "volume")? } - Metric::NodeProportionalVolume(idx) => { - let node = model.get_node(idx)?; - require_node_dataset(root_grp, shape, node.name(), node.sub_name(), "proportional-volume")? + Metric::DerivedMetric(idx) => { + todo!("Derived metrics are not yet supported in HDF recorders"); } Metric::AggregatedNodeVolume(idx) => { let node = model.get_aggregated_storage_node(idx)?; require_node_dataset(root_grp, shape, node.name(), node.sub_name(), "volume")? } - Metric::AggregatedNodeProportionalVolume(idx) => { - let node = model.get_aggregated_storage_node(idx)?; - require_node_dataset(root_grp, shape, node.name(), node.sub_name(), "proportional-volume")? - } Metric::EdgeFlow(_) => { continue; // TODO } @@ -118,9 +113,6 @@ impl Recorder for HDF5Recorder { Metric::VirtualStorageVolume(_) => { continue; // TODO } - Metric::VirtualStorageProportionalVolume(_) => { - continue; // TODO - } Metric::Constant(_) => { continue; // TODO } @@ -135,13 +127,6 @@ impl Recorder for HDF5Recorder { let node = model.get_aggregated_node(idx)?; require_node_dataset(root_grp, shape, node.name(), node.sub_name(), "outflow")? } - Metric::NodeInFlowDeficit(idx) => { - let node = model.get_node(idx)?; - require_node_dataset(root_grp, shape, node.name(), node.sub_name(), "flow-deficit")? - } - Metric::VolumeBetweenControlCurves(_) => { - continue; // TODO - } Metric::MultiNodeInFlow { name, sub_name, .. } => { require_node_dataset(root_grp, shape, name, sub_name.as_deref(), "inflow")? } diff --git a/pywr-core/src/recorders/mod.rs b/pywr-core/src/recorders/mod.rs index 45399d04..105c71c1 100644 --- a/pywr-core/src/recorders/mod.rs +++ b/pywr-core/src/recorders/mod.rs @@ -5,7 +5,7 @@ mod metric_set; mod py; pub use self::csv::CSVRecorder; -use crate::metric::Metric; +use crate::metric::{IndexMetric, Metric}; use crate::model::Model; use crate::scenario::ScenarioIndex; use crate::state::State; @@ -283,6 +283,62 @@ where } } +pub struct IndexAssertionRecorder { + meta: RecorderMeta, + expected_values: Array2, + metric: IndexMetric, +} + +impl IndexAssertionRecorder { + pub fn new(name: &str, metric: IndexMetric, expected_values: Array2) -> Self { + Self { + meta: RecorderMeta::new(name), + expected_values, + metric, + } + } +} + +impl Recorder for IndexAssertionRecorder { + fn meta(&self) -> &RecorderMeta { + &self.meta + } + + fn save( + &self, + timestep: &Timestep, + scenario_indices: &[ScenarioIndex], + model: &Model, + state: &[State], + _internal_state: &mut Option>, + ) -> Result<(), PywrError> { + // This panics if out-of-bounds + + for scenario_index in scenario_indices { + let expected_value = match self.expected_values.get([timestep.index, scenario_index.index]) { + Some(v) => *v, + None => panic!("Simulation produced results out of range."), + }; + + let actual_value = self.metric.get_value(model, &state[scenario_index.index])?; + + if actual_value != expected_value { + panic!( + r#"assertion failed: (actual eq expected) +recorder: `{}` +timestep: `{:?}` ({}) +scenario: `{:?}` +actual: `{:?}` +expected: `{:?}`"#, + self.meta.name, timestep.date, timestep.index, scenario_index.index, actual_value, expected_value + ) + } + } + + Ok(()) + } +} + pub enum RecorderAggregation { Min, Max, diff --git a/pywr-core/src/state.rs b/pywr-core/src/state.rs index 6aff528d..54dae7cc 100644 --- a/pywr-core/src/state.rs +++ b/pywr-core/src/state.rs @@ -1,3 +1,4 @@ +use crate::derived_metric::DerivedMetricIndex; use crate::edge::{Edge, EdgeIndex}; use crate::model::Model; use crate::node::{Node, NodeIndex}; @@ -503,6 +504,7 @@ impl NetworkState { pub struct State { network: NetworkState, parameters: ParameterValues, + derived_metrics: Vec, } impl State { @@ -513,10 +515,12 @@ impl State { num_parameter_values: usize, num_parameter_indices: usize, num_multi_parameters: usize, + num_derived_metrics: usize, ) -> Self { Self { network: NetworkState::new(initial_node_states, num_edges, initial_virtual_storage_states), parameters: ParameterValues::new(num_parameter_values, num_parameter_indices, num_multi_parameters), + derived_metrics: vec![0.0; num_derived_metrics], } } @@ -572,4 +576,21 @@ impl State { ) -> Result<(), PywrError> { self.network.reset_virtual_storage_volume(idx, volume, timestep) } + + pub fn get_derived_metric_value(&self, idx: DerivedMetricIndex) -> Result { + match self.derived_metrics.get(*idx.deref()) { + Some(s) => Ok(*s), + None => Err(PywrError::DerivedMetricIndexNotFound(idx)), + } + } + + pub fn set_derived_metric_value(&mut self, idx: DerivedMetricIndex, value: f64) -> Result<(), PywrError> { + match self.derived_metrics.get_mut(*idx.deref()) { + Some(s) => { + *s = value; + Ok(()) + } + None => Err(PywrError::DerivedMetricIndexNotFound(idx)), + } + } } diff --git a/pywr-schema/src/nodes/piecewise_storage.rs b/pywr-schema/src/nodes/piecewise_storage.rs index 1489e456..d9554fd0 100644 --- a/pywr-schema/src/nodes/piecewise_storage.rs +++ b/pywr-schema/src/nodes/piecewise_storage.rs @@ -2,8 +2,9 @@ use crate::data_tables::LoadedTableCollection; use crate::error::SchemaError; use crate::nodes::NodeMeta; use crate::parameters::DynamicFloatValue; -use pywr_core::metric::{Metric, VolumeBetweenControlCurves}; +use pywr_core::metric::Metric; use pywr_core::node::{ConstraintValue, StorageInitialVolume}; +use pywr_core::parameters::VolumeBetweenControlCurvesParameter; use std::path::Path; #[derive(serde::Deserialize, serde::Serialize, Clone)] @@ -79,9 +80,14 @@ impl PiecewiseStorageNode { let upper = step.control_curve.load(model, tables, data_path)?; - let max_volume = ConstraintValue::Metric(Metric::VolumeBetweenControlCurves( - VolumeBetweenControlCurves::new(max_volume.clone(), Some(upper), lower), - )); + let max_volume_parameter = VolumeBetweenControlCurvesParameter::new( + format!("{}-{}-max-volume", self.meta.name, Self::step_sub_name(i).unwrap()).as_str(), + max_volume.clone(), + Some(upper), + lower, + ); + let max_volume_parameter_idx = model.add_parameter(Box::new(max_volume_parameter))?; + let max_volume = ConstraintValue::Metric(Metric::ParameterValue(max_volume_parameter_idx)); // Each store has min volume of zero let min_volume = ConstraintValue::Scalar(0.0); @@ -113,11 +119,19 @@ impl PiecewiseStorageNode { let upper = None; - let max_volume = ConstraintValue::Metric(Metric::VolumeBetweenControlCurves(VolumeBetweenControlCurves::new( + let max_volume_parameter = VolumeBetweenControlCurvesParameter::new( + format!( + "{}-{}-max-volume", + self.meta.name, + Self::step_sub_name(self.steps.len()).unwrap() + ) + .as_str(), max_volume.clone(), upper, lower, - ))); + ); + let max_volume_parameter_idx = model.add_parameter(Box::new(max_volume_parameter))?; + let max_volume = ConstraintValue::Metric(Metric::ParameterValue(max_volume_parameter_idx)); // Each store has min volume of zero let min_volume = ConstraintValue::Scalar(0.0); @@ -181,9 +195,9 @@ impl PiecewiseStorageNode { #[cfg(test)] mod tests { use crate::model::PywrModel; - use ndarray::Array2; - use pywr_core::metric::Metric; - use pywr_core::recorders::AssertionRecorder; + use ndarray::{concatenate, Array, Array2, Axis}; + use pywr_core::metric::{IndexMetric, Metric}; + use pywr_core::recorders::{AssertionRecorder, IndexAssertionRecorder}; use pywr_core::test_utils::run_all_solvers; use pywr_core::timestep::Timestepper; @@ -251,7 +265,7 @@ mod tests { .get_aggregated_storage_node_index_by_name("storage1", None) .unwrap(); - let expected = Array2::from_shape_fn((366, 1), |(i, _)| { + let expected_volume = Array2::from_shape_fn((366, 1), |(i, _)| { if i < 49 { // Draw-down top store at 5 until it is emptied (control curve starts at 75%) 995.0 - i as f64 * 5.0 @@ -278,15 +292,34 @@ mod tests { } }); + // The drought index should register the time-step after storage is below 500 + let expected_drought_index = expected_volume.mapv(|v| if v < 500.0 { 1 } else { 0 }); + // The initial time-step has zero drought index because the initial volume is above the control curve + let initial_drought_index: Array = Array2::zeros((1, 1)); + + let expected_drought_index = + concatenate(Axis(0), &[initial_drought_index.view(), expected_drought_index.view()]).unwrap(); + let recorder = AssertionRecorder::new( "storage1-volume", Metric::AggregatedNodeVolume(idx), - expected, + expected_volume, None, None, ); model.add_recorder(Box::new(recorder)).unwrap(); + let idx = model + .get_index_parameter_index_by_name("storage1-drought-index") + .unwrap(); + + let recorder = IndexAssertionRecorder::new( + "storage1-drought-index", + IndexMetric::IndexParameterValue(idx), + expected_drought_index, + ); + model.add_recorder(Box::new(recorder)).unwrap(); + // Test all solvers run_all_solvers(&model, ×tepper); } diff --git a/pywr-schema/src/parameters/mod.rs b/pywr-schema/src/parameters/mod.rs index ad2228ed..ff18c7a7 100644 --- a/pywr-schema/src/parameters/mod.rs +++ b/pywr-schema/src/parameters/mod.rs @@ -44,6 +44,7 @@ use crate::error::{ConversionError, SchemaError}; use crate::parameters::core::DivisionParameter; pub use crate::parameters::data_frame::DataFrameParameter; pub use offset::OffsetParameter; +use pywr_core::derived_metric::DerivedMetric; use pywr_core::metric::Metric; use pywr_core::node::NodeIndex; use pywr_core::parameters::{IndexParameterIndex, IndexValue, ParameterType}; @@ -478,7 +479,8 @@ impl MetricFloatValue { Self::NodeOutFlow(node_ref) => Ok(Metric::NodeOutFlow(node_ref.get_node_index(model)?)), Self::NodeVolume(node_ref) => Ok(Metric::NodeVolume(node_ref.get_node_index(model)?)), Self::NodeProportionalVolume(node_ref) => { - Ok(Metric::NodeProportionalVolume(node_ref.get_node_index(model)?)) + let dm = DerivedMetric::NodeProportionalVolume(node_ref.get_node_index(model)?); + Ok(Metric::DerivedMetric(model.add_derived_metric(dm))) } Self::Parameter { name, key } => { match key { diff --git a/pywr-schema/src/parameters/polynomial.rs b/pywr-schema/src/parameters/polynomial.rs index 54a3dfcc..05a71839 100644 --- a/pywr-schema/src/parameters/polynomial.rs +++ b/pywr-schema/src/parameters/polynomial.rs @@ -1,6 +1,5 @@ use crate::error::{ConversionError, SchemaError}; use crate::parameters::{DynamicFloatValueType, IntoV2Parameter, ParameterMeta, TryFromV1Parameter}; -use pywr_core::metric::Metric; use pywr_core::parameters::ParameterIndex; use pywr_v1_schema::parameters::Polynomial1DParameter as Polynomial1DParameterV1; use std::collections::HashMap; @@ -25,12 +24,8 @@ impl Polynomial1DParameter { } pub fn add_to_model(&self, model: &mut pywr_core::model::Model) -> Result { - let node_idx = model.get_node_index_by_name(&self.storage_node, None)?; - let metric = if self.use_proportional_volume.unwrap_or(true) { - Metric::NodeProportionalVolume(node_idx) - } else { - Metric::NodeVolume(node_idx) - }; + let metric = + model.get_storage_node_metric(&self.storage_node, None, self.use_proportional_volume.unwrap_or(true))?; let p = pywr_core::parameters::Polynomial1DParameter::new( &self.meta.name, diff --git a/pywr-schema/src/test_models/piecewise_storage2.json b/pywr-schema/src/test_models/piecewise_storage2.json index 04a9fe49..a081341a 100644 --- a/pywr-schema/src/test_models/piecewise_storage2.json +++ b/pywr-schema/src/test_models/piecewise_storage2.json @@ -54,5 +54,23 @@ "from_node": "storage1", "to_node": "demand1" } + ], + "parameters": [ + { + "name": "storage1-drought-curve", + "type": "Constant", + "value": 0.5 + }, + { + "name": "storage1-drought-index", + "type": "ControlCurveIndex", + "storage_node": "storage1", + "control_curves": [ + { + "type": "Parameter", + "name": "storage1-drought-curve" + } + ] + } ] }