Commit 35897009 authored by Turner, Sean's avatar Turner, Sean
Browse files

initial and minimum pool elevation

parent b13fc6ae

examples/Cumberland/cascade_config.yaml

+20 −12
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@ 


# *Reservoir* objects have the following additional, unique attributes:

#   capacity: The storage capacity of the reservoir in Million Cubic Meters (MCM)

#   initial_storage: The initial storage in MCM

#   initial_pool_elevation: The initial elevation of the headwater (m)

#   max_release: the maximum release volume over one hour (MCM)

#   min_release: the minimum release volume over one hour (MCM)

#   power_params: parameters for multi-level model relating power (Watts)...
@@ -36,10 +36,11 @@ WolfCreek: 
  object_type: reservoir

  simulation_order: 1

  capacity: 8307

  initial_storage: 8000

  initial_pool_elevation: 220

  min_power_pool: 199

  tailwater_elevation: 168

  max_release: 0.94

  min_release: 0

  max_release: 2.5

  min_release: 0.1

  power_params:

    a: 0.893

    b: 0
@@ -50,7 +51,8 @@ DaleHollow: 
  object_type: reservoir

  simulation_order: 2

  capacity: 2181

  initial_storage: 2000

  initial_pool_elevation: 200

  min_power_pool: 192

  tailwater_elevation: 154

  max_release: 100

  min_release: 0
@@ -64,7 +66,8 @@ CordellHull: 
  object_type: reservoir

  simulation_order: 7

  capacity: 541

  initial_storage: 400

  initial_pool_elevation: 155

  min_power_pool: 152

  tailwater_elevation: 137

  max_release: 3

  min_release: 0
@@ -78,7 +81,8 @@ CenterHill: 
  object_type: reservoir

  simulation_order: 9

  capacity: 3017

  initial_storage: 3000

  initial_pool_elevation: 210

  min_power_pool: 188

  tailwater_elevation: 147

  max_release: 0.93

  min_release: 0
@@ -92,7 +96,8 @@ OldHickory: 
  object_type: reservoir

  simulation_order: 11

  capacity: 857

  initial_storage: 800

  initial_pool_elevation: 147

  min_power_pool: 146

  tailwater_elevation: 121

  max_release: 3

  min_release: 0
@@ -106,7 +111,8 @@ JPercyPriest: 
  object_type: reservoir

  simulation_order: 13

  capacity: 1414

  initial_storage: 1400

  initial_pool_elevation: 150

  min_power_pool: 129

  tailwater_elevation: 121

  max_release: 3

  min_release: 0
@@ -120,7 +126,8 @@ Cheatham: 
  object_type: reservoir

  simulation_order: 15

  capacity: 1223

  initial_storage: 1200

  initial_pool_elevation: 117

  min_power_pool: 116

  tailwater_elevation: 111

  max_release: 3

  min_release: 0
@@ -134,7 +141,8 @@ Barkley: 
  object_type: reservoir

  simulation_order: 16

  capacity: 2460

  initial_storage: 2000

  initial_pool_elevation: 109

  min_power_pool: 107

  tailwater_elevation: 96

  max_release: 3

  min_release: 0

examples/Cumberland/script.qmd

+20 −22
Original line number Diff line number Diff line
@@ -31,25 +31,22 @@ from pathlib import Path 
class ReservoirData:

    object_type: str

    capacity: float

    initial_storage: float

    initial_pool_elevation: float

    min_power_pool: float

    # storage-headwater elevation tables:

    set_storage: list[float]

    set_elevation: list[float]

    # outflow-tailwater elevation tables:

    #   otet_outflow: list[float]

    #   otet_elevation: list[float]

    # currently we use static tailwater elevation

    tailwater_elevation: float

    max_release: float

    min_release: float

    catchment_inflow: list[float]

    target_release: list[float]

    # target_release: list[float]

    target_power: list[float]

    # head-release-power tables:

    hrpt_head: list[float]

    hrpt_release: list[float]

    hrpt_power: list[float]

    simulation_order: int

    downstream_object: str

    power_params_a: float

    power_params_b: float



@dataclass

class RiverData:
@@ -80,22 +77,26 @@ reservoir_dict = {} 
river_dict = {}

confluence_dict = {}



# Process power_params for all reservoirs

for name, specs in config_dict.items():

    if specs['object_type'] == 'reservoir' and 'power_params' in specs:

        power_params = specs['power_params']

        specs['power_params_a'] = float(power_params['a'])

        specs['power_params_b'] = float(power_params['b'])

        del specs['power_params']



for name, specs in config_dict.items():

    #specs['simulation_order'] = simulation_orders[name]  # Add simulation order to specs

    if specs['object_type'] == 'reservoir':



    if specs['object_type'] == 'reservoir':

        reservoir_dict[name] = ReservoirData(

            **specs,

            catchment_inflow=pd.read_csv(f"time_series/{name}.csv")['catchment_inflow'].tolist(),

            target_release=pd.read_csv(f"time_series/{name}.csv")['target_release'].tolist(),

            catchment_inflow = [0] * 168,

            #catchment_inflow=pd.read_csv(f"time_series/{name}.csv")['catchment_inflow'].tolist(),

            # target_release=pd.read_csv(f"time_series/{name}.csv")['target_release'].tolist(),

            target_power=pd.read_csv(f"time_series/{name}.csv")['target_power'].tolist(),

            # storage elevation tables

            set_storage = pd.read_csv(f"storage_HWelevation_tables/{name}.csv")['storage_Mm3'].tolist(),

            set_elevation = pd.read_csv(f"storage_HWelevation_tables/{name}.csv")['elevation_m'].tolist(),

            # head-release-power tables

            hrpt_head = pd.read_csv(f"head_release_power_tables/{name}.csv")['head_m'].tolist(),

            hrpt_release = pd.read_csv(f"head_release_power_tables/{name}.csv")['release_Mm3'].tolist(),

            hrpt_power = pd.read_csv(f"head_release_power_tables/{name}.csv")['power_MW'].tolist(),

        )

    elif specs['object_type'] == 'river':

        river_dict[name] = RiverData(
@@ -133,8 +134,6 @@ reservoirs_ordered |> map_dfr(function(reservoir){ 


    tibble(

        reservoir = !!reservoir,

        #flow = reservoir_results$inflow,

        #outflow = reservoir_results$outflow,

        inflow = reservoir_results$total_inflow,

        target_release = reservoir_results$target_release,

        storage = reservoir_results$storage,
@@ -159,7 +158,6 @@ powersheds_plot_theme <- function(){ 
plot_cols <- length(reservoirs_ordered)



ts_data_for_plot |>

    #filter(!grepl("_", reservoir), reservoir != "Celina") |>

    select(reservoir, hr_index, inflow, target_release, release_implemented, spill) |>

    tidyr::pivot_longer(-c(hr_index, reservoir), names_to = "variable") |>

    ggplot(aes(hr_index, value, col = variable)) +
@@ -168,12 +166,12 @@ ts_data_for_plot |> 
    facet_wrap(~reservoir, ncol = plot_cols) -> flows_plot



ts_data_for_plot |>

    select(reservoir, hr_index, head, pool_elevation) |>

    select(reservoir, hr_index, pool_elevation) |>

    tidyr::pivot_longer(-c(hr_index, reservoir), names_to = "variable") |>

    ggplot(aes(hr_index, value, col = variable)) +

    geom_line() +

    powersheds_plot_theme() +

    facet_wrap(~reservoir, ncol = plot_cols) -> elevation_plot

    facet_wrap(~reservoir, ncol = plot_cols, scales = "free") -> elevation_plot



ts_data_for_plot |>

    #filter(!grepl("_", reservoir), reservoir != "Celina") |>

examples/Cumberland/time_series/WolfCreek.csv

+72 −72
Original line number Diff line number Diff line 
datetime,catchment_inflow,target_release,target_power

2007-07-01T00:00:00Z,8.42,0,0

2007-07-01T01:00:00Z,8.42,0,0

2007-07-01T02:00:00Z,8.42,0,0

2007-07-01T03:00:00Z,8.42,0,0

2007-07-01T04:00:00Z,8.42,0,0

2007-07-01T05:00:00Z,8.42,0,0

2007-07-01T06:00:00Z,8.42,0,0

2007-07-01T07:00:00Z,8.42,0,0

2007-07-01T08:00:00Z,8.42,0,0

2007-07-01T09:00:00Z,8.42,0,0

2007-07-01T10:00:00Z,8.42,0,0

2007-07-01T11:00:00Z,8.42,0,0

2007-07-01T12:00:00Z,8.42,56.27,276.67

2007-07-01T13:00:00Z,8.42,56.27,276.67

2007-07-01T14:00:00Z,8.42,56.27,276.67

2007-07-01T15:00:00Z,8.42,56.27,276.67

2007-07-01T16:00:00Z,8.42,56.27,276.67

2007-07-01T17:00:00Z,8.42,56.27,276.67

2007-07-01T18:00:00Z,8.42,13.9,73

2007-07-01T19:00:00Z,8.42,13.9,73

2007-07-01T20:00:00Z,8.42,13.9,73

2007-07-01T21:00:00Z,8.42,13.9,73

2007-07-01T22:00:00Z,8.42,13.9,73

2007-07-01T23:00:00Z,8.42,13.9,73

2007-07-02T00:00:00Z,11.4,0,0

2007-07-02T01:00:00Z,11.4,0,0

2007-07-02T02:00:00Z,11.4,0,0

2007-07-02T03:00:00Z,11.4,0,0

2007-07-02T04:00:00Z,11.4,0,0

2007-07-02T05:00:00Z,11.4,0,0

2007-07-02T06:00:00Z,11.4,0,0

2007-07-02T07:00:00Z,11.4,0,0

2007-07-02T08:00:00Z,11.4,0,0

2007-07-02T09:00:00Z,11.4,0,0

2007-07-02T10:00:00Z,11.4,0,0

2007-07-02T11:00:00Z,11.4,0,0

2007-07-02T12:00:00Z,11.4,54.78,270

2007-07-02T13:00:00Z,11.4,54.78,270

2007-07-02T14:00:00Z,11.4,54.78,270

2007-07-02T15:00:00Z,11.4,54.78,270

2007-07-02T16:00:00Z,11.4,54.78,270

2007-07-02T17:00:00Z,11.4,54.78,270

2007-07-02T18:00:00Z,11.4,27.3,140.17

2007-07-02T19:00:00Z,11.4,27.3,140.17

2007-07-02T20:00:00Z,11.4,27.3,140.17

2007-07-02T21:00:00Z,11.4,27.3,140.17

2007-07-02T22:00:00Z,11.4,27.3,140.17

2007-07-02T23:00:00Z,11.4,27.3,140.17

2007-07-03T00:00:00Z,11.55,0,0

2007-07-03T01:00:00Z,11.55,0,0

2007-07-03T02:00:00Z,11.55,0,0

2007-07-03T03:00:00Z,11.55,0,0

2007-07-03T04:00:00Z,11.55,0,0

2007-07-03T05:00:00Z,11.55,0,0

2007-07-03T06:00:00Z,11.55,0,0

2007-07-03T07:00:00Z,11.55,0,0

2007-07-03T08:00:00Z,11.55,0,0

2007-07-03T09:00:00Z,11.55,0,0

2007-07-03T10:00:00Z,11.55,0,0

2007-07-03T11:00:00Z,11.55,0,0

2007-07-03T12:00:00Z,11.55,54.85,270

2007-07-03T13:00:00Z,11.55,54.85,270

2007-07-03T14:00:00Z,11.55,54.85,270

2007-07-03T15:00:00Z,11.55,54.85,270

2007-07-03T16:00:00Z,11.55,54.85,270

2007-07-03T17:00:00Z,11.55,54.85,270

2007-07-03T18:00:00Z,11.55,27.84,142.67

2007-07-03T19:00:00Z,11.55,27.84,142.67

2007-07-03T20:00:00Z,11.55,27.84,142.67

2007-07-03T21:00:00Z,11.55,27.84,142.67

2007-07-03T22:00:00Z,11.55,27.84,142.67

2007-07-03T23:00:00Z,11.55,27.84,142.67

2007-07-01T00:00:00Z,0,0,0

2007-07-01T01:00:00Z,0,0,0

2007-07-01T02:00:00Z,0,0,0

2007-07-01T03:00:00Z,0,0,0

2007-07-01T04:00:00Z,0,0,0

2007-07-01T05:00:00Z,0,0,0

2007-07-01T06:00:00Z,0,0,0

2007-07-01T07:00:00Z,0,0,0

2007-07-01T08:00:00Z,0,0,0

2007-07-01T09:00:00Z,0,0,0

2007-07-01T10:00:00Z,0,0,0

2007-07-01T11:00:00Z,0,0,0

2007-07-01T12:00:00Z,0,56.27,276.67

2007-07-01T13:00:00Z,0,56.27,276.67

2007-07-01T14:00:00Z,0,56.27,276.67

2007-07-01T15:00:00Z,0,56.27,276.67

2007-07-01T16:00:00Z,0,56.27,276.67

2007-07-01T17:00:00Z,0,56.27,276.67

2007-07-01T18:00:00Z,0,13.9,73

2007-07-01T19:00:00Z,0,13.9,73

2007-07-01T20:00:00Z,0,13.9,73

2007-07-01T21:00:00Z,0,13.9,73

2007-07-01T22:00:00Z,0,13.9,73

2007-07-01T23:00:00Z,0,13.9,73

2007-07-02T00:00:00Z,0,0,0

2007-07-02T01:00:00Z,0,0,0

2007-07-02T02:00:00Z,0,0,0

2007-07-02T03:00:00Z,0,0,0

2007-07-02T04:00:00Z,0,0,0

2007-07-02T05:00:00Z,0,0,0

2007-07-02T06:00:00Z,0,0,0

2007-07-02T07:00:00Z,0,0,0

2007-07-02T08:00:00Z,0,0,0

2007-07-02T09:00:00Z,0,0,0

2007-07-02T10:00:00Z,0,0,0

2007-07-02T11:00:00Z,0,0,0

2007-07-02T12:00:00Z,0,54.78,270

2007-07-02T13:00:00Z,0,54.78,270

2007-07-02T14:00:00Z,0,54.78,270

2007-07-02T15:00:00Z,0,54.78,270

2007-07-02T16:00:00Z,0,54.78,270

2007-07-02T17:00:00Z,0,54.78,270

2007-07-02T18:00:00Z,0,27.3,140.17

2007-07-02T19:00:00Z,0,27.3,140.17

2007-07-02T20:00:00Z,0,27.3,140.17

2007-07-02T21:00:00Z,0,27.3,140.17

2007-07-02T22:00:00Z,0,27.3,140.17

2007-07-02T23:00:00Z,0,27.3,140.17

2007-07-03T00:00:00Z,0,0,0

2007-07-03T01:00:00Z,0,0,0

2007-07-03T02:00:00Z,0,0,0

2007-07-03T03:00:00Z,0,0,0

2007-07-03T04:00:00Z,0,0,0

2007-07-03T05:00:00Z,0,0,0

2007-07-03T06:00:00Z,0,0,0

2007-07-03T07:00:00Z,0,0,0

2007-07-03T08:00:00Z,0,0,0

2007-07-03T09:00:00Z,0,0,0

2007-07-03T10:00:00Z,0,0,0

2007-07-03T11:00:00Z,0,0,0

2007-07-03T12:00:00Z,0,54.85,270

2007-07-03T13:00:00Z,0,54.85,270

2007-07-03T14:00:00Z,0,54.85,270

2007-07-03T15:00:00Z,0,54.85,270

2007-07-03T16:00:00Z,0,54.85,270

2007-07-03T17:00:00Z,0,54.85,270

2007-07-03T18:00:00Z,0,27.84,142.67

2007-07-03T19:00:00Z,0,27.84,142.67

2007-07-03T20:00:00Z,0,27.84,142.67

2007-07-03T21:00:00Z,0,27.84,142.67

2007-07-03T22:00:00Z,0,27.84,142.67

2007-07-03T23:00:00Z,0,27.84,142.67

2007-07-04T00:00:00Z,7.05,0,0

2007-07-04T01:00:00Z,7.05,0,0

2007-07-04T02:00:00Z,7.05,0,0

src/helpers.rs

+39 −1
Original line number Diff line number Diff line
@@ -12,7 +12,10 @@ const MEGAWATT_TO_WATT: f64 = 1_000_000.0; 
/// Elevation computation based on storage.

/// Performs linear interpolation to find a y-value (elevation)

/// corresponding to an x-value (storage)

pub fn interpolate_elevation(storage: f64, storage_values: &[f64], elevation_values: &[f64]) -> f64 {

pub fn interpolate_elevation(

    storage: f64,

    storage_values: &[f64],

    elevation_values: &[f64]) -> f64 {

    // Handle edge cases

    if storage_values.is_empty() || elevation_values.is_empty() {

        return 0.0; // Default value if tables are empty
@@ -43,6 +46,41 @@ pub fn interpolate_elevation(storage: f64, storage_values: &[f64], elevation_val 
    elevation_values[0] // Fallback

}



pub fn interpolate_storage(

    elevation: f64,

    storage_values: &[f64],

    elevation_values: &[f64]) -> f64 {

    // Handle edge cases

    if storage_values.is_empty() || elevation_values.is_empty() {

        return 0.0; // Default value if tables are empty

    }

    

    if elevation <= elevation_values[0] {

        return storage_values[0]; // Below minimum elevation

    }

    

    let last_idx = elevation_values.len() - 1;

    if elevation >= elevation_values[last_idx] {

        return storage_values[last_idx]; // Above maximum elevation

    }

    

    // Find the bounding indices

    for i in 0..last_idx {

        if elevation_values[i] <= elevation && elevation <= elevation_values[i + 1] {

            // Linear interpolation formula: y = y1 + (x - x1) * (y2 - y1) / (x2 - x1)

            let x1 = elevation_values[i];

            let x2 = elevation_values[i + 1];

            let y1 = storage_values[i];

            let y2 = storage_values[i + 1];

            

            return y1 + (elevation - x1) * (y2 - y1) / (x2 - x1);

        }

    }

    

    storage_values[0] // Fallback

}







/// Computes volumetric release (Mm3) required to meet target power

/// ...

src/lib.rs

+30 −21
Original line number Diff line number Diff line
@@ -6,6 +6,7 @@ use std::str::FromStr; 


mod helpers;

use helpers::interpolate_elevation;

use helpers::interpolate_storage;

use helpers::compute_release_to_meet_target_power;

use helpers::compute_power;
@@ -32,14 +33,14 @@ impl FromStr for ObjectType { 
struct ReservoirData {

    object_type: String,

    capacity: f64,

    initial_storage: f64,

    initial_pool_elevation: f64,

    min_power_pool: f64,

    set_storage: Vec<f64>,

    set_elevation: Vec<f64>,

    tailwater_elevation: f64,

    max_release: f64,

    min_release: f64,

    catchment_inflow: Vec<f64>,

    ///target_release: Vec<f64>,

    target_power: Vec<f64>,

    power_params_a: f64,

    power_params_b: f64,
@@ -121,21 +122,20 @@ fn simulate_timestep( 
    previous_storage: f64,

) -> ReservoirState {



    // THIS WILL NEED TO BE IMPROVED WITH AN ITERATIVE PROCEDURE



    // STEP 1. Use previous storage to estimate pool elevation, then head

    // Compute pool elevation by interpolating storage elevation table for this reservoir

    // Use previous storage to estimate pool elevation, then head

    let previous_pool_elevation = interpolate_elevation(

        previous_storage, 

        &reservoir.set_storage, 

        &reservoir.set_elevation

    );



    // Determine whether the pool has dropped below the power pool minimum

    let below_power_pool = previous_pool_elevation < reservoir.min_power_pool;



    // Calculate head as the difference between pool elevation and tailwater elevation

    let head = previous_pool_elevation - reservoir.tailwater_elevation;



    // STEP 2. COMPUTE THE TARGET RELEASE GIVEN THE POWER TARGETED



    // Compute the target release, given the power generation desired

    let target_release = compute_release_to_meet_target_power(

        target_power,

        head,
@@ -143,14 +143,18 @@ fn simulate_timestep( 
        reservoir.power_params_b,

    );

    

    // STEP 3. COMPUTE AVAILABLE WATER GIVEN INFLOW, THEN DETERMINE ACTUAL RELEASE AND STORAGE

    // Compute available water given inflow and previous storage...

    // ... then determine actual release and resulting storage

    let available_water = previous_storage + inflow;



    let release = target_release

    let release = if below_power_pool {

        0.0

    } else {

        target_release

            .max(reservoir.min_release)

            .min(reservoir.max_release)

        .min(available_water);

    

            .min(available_water)

    };

    

    let potential_storage = available_water - release;
@@ -162,13 +166,14 @@ fn simulate_timestep( 
        reservoir.power_params_b,

    );



    // Determine spill--driven solely by storage capacity limit

    let (storage, spill) = if potential_storage > reservoir.capacity {

        (reservoir.capacity, potential_storage - reservoir.capacity)

    } else {

        (potential_storage, 0.0)

    };



    // STEP 4. COMPUTE RESULTING POOL ELEVATION

    // Compute final pool elevation



    let pool_elevation = interpolate_elevation(

        storage,
@@ -222,10 +227,10 @@ fn simulate_cascade( 
        match obj_type {

            ObjectType::Reservoir => {

                let reservoir = &cascade_data.reservoirs[*name];

                let mut storage = vec![0.0; n];

                storage[0] = reservoir.initial_storage;

                //let mut storage = vec![0.0; n];

                //storage[0] = reservoir.initial_storage;

                results.insert((*name).to_string(), CascadeResults::Reservoir(ReservoirResults {

                    storage,

                    storage: vec![0.0; n],

                    pool_elevation: vec![0.0; n],

                    head: vec![0.0; n],

                    release: vec![0.0; n],
@@ -300,14 +305,18 @@ fn simulate_cascade( 
            match obj_type {

                ObjectType::Reservoir => {

                    let reservoir = &cascade_data.reservoirs[*name];

                    let initial_storage = interpolate_storage(

                        reservoir.initial_pool_elevation,

                        &reservoir.set_storage,

                        &reservoir.set_elevation,                                

                        );

                    if let CascadeResults::Reservoir(current_results) = results.get_mut(*name).unwrap() {

                        let current_inflow = reservoir.catchment_inflow[t] + upstream_inflow;

                        let previous_storage = if t == 0 {

                            reservoir.initial_storage 

                            initial_storage

                        } else { 

                            current_results.storage[t-1] 

                        };



                        let state = simulate_timestep(

                            reservoir,

                            current_inflow,