feat: Add Turbine node and parameter (#144)

* Added HydropowerTargetParameter
* Added PowerFromNodeFlow derived metric
* Added TurbineNode in pywr-schema

pywr-core/src/derived_metric.rs

@@ -1,8 +1,10 @@
 use crate::aggregated_storage_node::AggregatedStorageNodeIndex;
+use crate::metric::MetricF64;
 use crate::network::Network;
 use crate::node::NodeIndex;
 use crate::state::State;
 use crate::timestep::Timestep;
+use crate::utils::hydropower_calculation;
 use crate::virtual_storage::VirtualStorageIndex;
 use crate::PywrError;
 use std::fmt;
@@ -32,6 +34,23 @@ impl Display for DerivedMetricIndex {
     }
 }
 
+// Turbine data to calculate the power in the `PowerFromNodeFlow` metric.
+#[derive(Clone, Debug, PartialEq)]
+pub struct TurbineData {
+    // The turbine elevation
+    pub elevation: f64,
+    // The turbine relative efficiency (0-1)
+    pub efficiency: f64,
+    // The water elevation above the turbine
+    pub water_elevation: Option<MetricF64>,
+    // The water density
+    pub water_density: f64,
+    /// A factor used to transform the units of flow to be compatible with the hydropower equation
+    pub flow_unit_conversion: f64,
+    /// A factor used to transform the units of total energy
+    pub energy_unit_conversion: f64,
+}
+
 /// Derived metrics are updated after the model is solved.
 ///
 /// These metrics are "derived" from node states (e.g. volume, flow) and must be updated
@@ -43,6 +62,7 @@ pub enum DerivedMetric {
     NodeProportionalVolume(NodeIndex),
     AggregatedNodeProportionalVolume(AggregatedStorageNodeIndex),
     VirtualStorageProportionalVolume(VirtualStorageIndex),
+    PowerFromNodeFlow(NodeIndex, TurbineData),
 }
 
 impl DerivedMetric {
@@ -91,20 +111,42 @@ impl DerivedMetric {
                 let max_flow = node.get_current_max_flow(network, state)?;
                 Ok(max_flow - flow)
             }
+            Self::PowerFromNodeFlow(idx, turbine_data) => {
+                let flow = state.get_network_state().get_node_in_flow(idx)?;
+
+                // Calculate the head (the head may be negative)
+                let head = if let Some(water_elevation) = &turbine_data.water_elevation {
+                    water_elevation.get_value(network, state)? - turbine_data.elevation
+                } else {
+                    turbine_data.elevation
+                }
+                .max(0.0);
+
+                Ok(hydropower_calculation(
+                    flow,
+                    head,
+                    turbine_data.efficiency,
+                    turbine_data.flow_unit_conversion,
+                    turbine_data.energy_unit_conversion,
+                    turbine_data.water_density,
+                ))
+            }
         }
     }
 
     pub fn name<'a>(&self, network: &'a Network) -> Result<&'a str, PywrError> {
         match self {
-            Self::NodeInFlowDeficit(idx) | Self::NodeProportionalVolume(idx) => network.get_node(idx).map(|n| n.name()),
+            Self::NodeInFlowDeficit(idx) | Self::NodeProportionalVolume(idx) | Self::PowerFromNodeFlow(idx, _) => {
+                network.get_node(idx).map(|n| n.name())
+            }
             Self::AggregatedNodeProportionalVolume(idx) => network.get_aggregated_storage_node(idx).map(|n| n.name()),
             Self::VirtualStorageProportionalVolume(idx) => network.get_virtual_storage_node(idx).map(|v| v.name()),
         }
     }
 
     pub fn sub_name<'a>(&self, network: &'a Network) -> Result<Option<&'a str>, PywrError> {
         match self {
-            Self::NodeInFlowDeficit(idx) | Self::NodeProportionalVolume(idx) => {
+            Self::NodeInFlowDeficit(idx) | Self::NodeProportionalVolume(idx) | Self::PowerFromNodeFlow(idx, _) => {
                 network.get_node(idx).map(|n| n.sub_name())
             }
             Self::AggregatedNodeProportionalVolume(idx) => {
@@ -120,6 +162,7 @@ impl DerivedMetric {
             Self::NodeProportionalVolume(_) => "proportional_volume",
             Self::AggregatedNodeProportionalVolume(_) => "proportional_volume",
             Self::VirtualStorageProportionalVolume(_) => "proportional_volume",
+            Self::PowerFromNodeFlow(_, _) => "power_from_flow",
         }
     }
 }

pywr-core/src/lib.rs

@@ -28,6 +28,7 @@ pub mod solvers;
 pub mod state;
 pub mod test_utils;
 pub mod timestep;
+pub mod utils;
 pub mod virtual_storage;
 
 #[derive(Error, Debug, PartialEq, Eq)]

pywr-core/src/parameters/hydropower.rs

@@ -0,0 +1,109 @@
+use crate::metric::MetricF64;
+use crate::network::Network;
+use crate::parameters::{Parameter, ParameterMeta};
+use crate::scenario::ScenarioIndex;
+use crate::state::{ParameterState, State};
+use crate::timestep::Timestep;
+use crate::utils::inverse_hydropower_calculation;
+use crate::PywrError;
+
+pub struct HydropowerTargetData {
+    pub target: MetricF64,
+    pub elevation: Option<f64>,
+    pub min_head: Option<f64>,
+    pub max_flow: Option<MetricF64>,
+    pub min_flow: Option<MetricF64>,
+    pub efficiency: Option<f64>,
+    pub water_elevation: Option<MetricF64>,
+    pub water_density: Option<f64>,
+    pub flow_unit_conversion: Option<f64>,
+    pub energy_unit_conversion: Option<f64>,
+}
+
+pub struct HydropowerTargetParameter {
+    pub meta: ParameterMeta,
+    pub target: MetricF64,
+    pub max_flow: Option<MetricF64>,
+    pub min_flow: Option<MetricF64>,
+    pub turbine_min_head: f64,
+    pub turbine_elevation: f64,
+    pub turbine_efficiency: f64,
+    pub water_elevation: Option<MetricF64>,
+    pub water_density: f64,
+    pub flow_unit_conversion: f64,
+    pub energy_unit_conversion: f64,
+}
+
+impl HydropowerTargetParameter {
+    pub fn new(name: &str, turbine_data: HydropowerTargetData) -> Self {
+        Self {
+            meta: ParameterMeta::new(name),
+            target: turbine_data.target,
+            water_elevation: turbine_data.water_elevation,
+            turbine_elevation: turbine_data.elevation.unwrap_or(0.0),
+            turbine_min_head: turbine_data.min_head.unwrap_or(0.0),
+            turbine_efficiency: turbine_data.efficiency.unwrap_or(1.0),
+            max_flow: turbine_data.max_flow,
+            min_flow: turbine_data.min_flow,
+            water_density: turbine_data.water_density.unwrap_or(1000.0),
+            flow_unit_conversion: turbine_data.flow_unit_conversion.unwrap_or(1.0),
+            energy_unit_conversion: turbine_data.energy_unit_conversion.unwrap_or(1e-6),
+        }
+    }
+}
+
+impl Parameter<f64> for HydropowerTargetParameter {
+    fn meta(&self) -> &ParameterMeta {
+        &self.meta
+    }
+    fn compute(
+        &self,
+        _timestep: &Timestep,
+        _scenario_index: &ScenarioIndex,
+        model: &Network,
+        state: &State,
+        _internal_state: &mut Option<Box<dyn ParameterState>>,
+    ) -> Result<f64, PywrError> {
+        // Calculate the head
+        let mut head = if let Some(water_elevation) = &self.water_elevation {
+            water_elevation.get_value(model, state)? - self.turbine_elevation
+        } else {
+            self.turbine_elevation
+        };
+
+        // the head may be negative
+        head = head.max(0.0);
+
+        // apply the minimum head threshold
+        if head <= self.turbine_min_head {
+            return Ok(0.0);
+        }
+
+        // Get the flow from the current node
+        let power = self.target.get_value(model, state)?;
+        let mut q = inverse_hydropower_calculation(
+            power,
+            head,
+            self.turbine_efficiency,
+            self.flow_unit_conversion,
+            self.energy_unit_conversion,
+            self.water_density,
+        );
+
+        // Bound the flow if required
+        if let Some(max_flow) = &self.max_flow {
+            q = q.min(max_flow.get_value(model, state)?);
+        }
+        if let Some(min_flow) = &self.min_flow {
+            q = q.max(min_flow.get_value(model, state)?);
+        }
+
+        if q < 0.0 {
+            return Err(PywrError::InternalParameterError(format!(
+                "The calculated flow in the hydro power parameter named {} is negative",
+                self.name()
+            )));
+        }
+        Ok(q)
+    }
+}

pywr-core/src/parameters/mod.rs

@@ -8,6 +8,7 @@ mod control_curves;
 mod delay;
 mod discount_factor;
 mod division;
+mod hydropower;
 mod indexed_array;
 mod interpolate;
 mod interpolated;
@@ -45,6 +46,7 @@ pub use control_curves::{
 pub use delay::DelayParameter;
 pub use discount_factor::DiscountFactorParameter;
 pub use division::DivisionParameter;
+pub use hydropower::{HydropowerTargetData, HydropowerTargetParameter};
 pub use indexed_array::IndexedArrayParameter;
 pub use interpolate::{interpolate, linear_interpolation, InterpolationError};
 pub use interpolated::InterpolatedParameter;

pywr-core/src/utils.rs

@@ -0,0 +1,23 @@
+/// Calculate the flow required to produce power using the hydropower equation
+pub fn inverse_hydropower_calculation(
+    power: f64,
+    head: f64,
+    efficiency: f64,
+    flow_unit_conversion: f64,
+    energy_unit_conversion: f64,
+    density: f64,
+) -> f64 {
+    power / (energy_unit_conversion * density * 9.81 * head * efficiency * flow_unit_conversion)
+}
+
+/// Calculate the produced power from the flow using the hydropower equation
+pub fn hydropower_calculation(
+    flow: f64,
+    head: f64,
+    efficiency: f64,
+    flow_unit_conversion: f64,
+    energy_unit_conversion: f64,
+    density: f64,
+) -> f64 {
+    flow * (energy_unit_conversion * density * 9.81 * head * efficiency * flow_unit_conversion)
+}

pywr-schema/src/metric.rs

@@ -268,7 +268,7 @@ impl NodeReference {
             .get_node_by_name(&self.name)
             .ok_or_else(|| SchemaError::NodeNotFound(self.name.clone()))?;
 
-        node.create_metric(network, self.attribute)
+        node.create_metric(network, self.attribute, args)
     }
 
     /// Return the attribute of the node. If the attribute is not specified then the default