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

cumberland eg update

parent 790ad13f
Loading
Loading
Loading
Loading
+28 −10
Original line number Diff line number Diff line
---
title: "Untitled"
title: "Powersheds Cumberland Example"
format: html
---

@@ -44,17 +44,33 @@ df = pd.read_csv("inflow_release_ex1.csv")
inflow = df['inflow'].tolist()
release = df['release'].tolist()

storage, actual_release, spill = powersheds.simulate_cascade(inflow, release, cordell_hull)
wc_storage, wc_release, wc_spill, ch_storage, ch_release, ch_spill = powersheds.simulate_cascade(
    inflow, 
    release, 
    wolf_creek, 
    cordell_hull
)
```


```{r}
tibble(
    reservoir = "Wolf Creek",
    inflow = py$inflow,
    release_target = py$release,
    storage = py$storage,
    release_implemented = py$actual_release,
    spill = py$spill
    storage = py$wc_storage,
    release_implemented = py$wc_release,
    spill = py$wc_spill
    ) |>
    bind_rows(
        tibble(
            reservoir = "Cordell Hull",
            inflow = py$wc_release,
            release_target = py$release,
            storage = py$ch_storage,
            release_implemented = py$ch_release,
            spill = py$ch_spill
        )
    ) |>
    mutate(hr_index = 1:n()) ->
    ts_data_for_plot
@@ -69,17 +85,19 @@ powersheds_plot_theme <- function(){
}

ts_data_for_plot |>
    select(hr_index, inflow, release_target, release_implemented, spill) |>
    tidyr::pivot_longer(-hr_index, names_to = "variable") |>
    select(reservoir, hr_index, inflow, release_target, release_implemented, spill) |>
    tidyr::pivot_longer(-c(hr_index, reservoir), names_to = "variable") |>
    ggplot(aes(hr_index, value, col = variable)) +
    geom_line() +
    powersheds_plot_theme() -> flows_plot
    powersheds_plot_theme() +
    facet_wrap(~reservoir, ncol = 2) -> flows_plot

ts_data_for_plot |>
    select(hr_index, storage) |>
    select(reservoir, hr_index, storage) |>
    ggplot(aes(hr_index, storage)) +
    geom_line() +
    powersheds_plot_theme() -> storage_plot
    powersheds_plot_theme() +
    facet_wrap(~reservoir, ncol = 2) -> storage_plot

flows_plot + storage_plot + plot_layout(ncol = 1)

+70 −32
Original line number Diff line number Diff line
@@ -8,53 +8,91 @@ struct ReservoirSpecs {
    max_release: f64,
    min_release: f64,
}
struct ReservoirState {
    storage: f64,
    release: f64,
    spill: f64,
}

fn simulate_timestep(
    specs: &ReservoirSpecs,
    inflow: f64,
    release_target: f64,
    previous_storage: f64,
) -> ReservoirState {
    let available_water = previous_storage + inflow;

    let release = release_target
        .max(specs.min_release)
        .min(specs.max_release)
        .min(available_water);

    let potential_storage = available_water - release;
    
    let (storage, spill) = if potential_storage > specs.capacity {
        (specs.capacity, potential_storage - specs.capacity)
    } else {
        (potential_storage, 0.0)
    };

    ReservoirState {
        storage,
        release,
        spill,
    }
}

/// Simulates a cascade of reservoirs
#[pyfunction]
fn simulate_cascade(inflow: Vec<f64>,
                    release: Vec<f64>,
                    specs: ReservoirSpecs) -> PyResult<(Vec<f64>, Vec<f64>, Vec<f64>)> {
                    release_target: Vec<f64>,
                    wolf_creek: ReservoirSpecs,
                    cordell_hull: ReservoirSpecs) -> PyResult<(Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>)> {


    if inflow.len() != release.len() {
    if inflow.len() != release_target.len() {
        return Err(PyValueError::new_err("Input vectors must have the same length"));
    }

    let n = inflow.len();
    let mut storage = vec![0.0; n];  // Storage at end of timestep
    let mut actual_release = vec![0.0; n];  // Actual release for timestep
    let mut spill = vec![0.0; n];  // Spill for timestep
    let mut wc_storage = vec![0.0; n];
    let mut wc_release = vec![0.0; n];
    let mut wc_spill = vec![0.0; n];
    let mut ch_storage = vec![0.0; n];
    let mut ch_release = vec![0.0; n];
    let mut ch_spill = vec![0.0; n];

    // Initialize first timestep storage
    storage[0] = specs.initial_storage;
    wc_storage[0] = wolf_creek.initial_storage;
    ch_storage[0] = cordell_hull.initial_storage;

    for t in 0..n {
        let available_water = if t == 0 {
            specs.initial_storage + inflow[t]
        } else {
            storage[t-1] + inflow[t]
        };
        // Wolf Creek simulation
        let wc_state = simulate_timestep(
            &wolf_creek,
            inflow[t],
            release_target[t],
            if t == 0 { wolf_creek.initial_storage } else { wc_storage[t-1] },
        );
        
        // Calculate actual release (constrained by available water and max/min release)
        actual_release[t] = release[t]
            .max(specs.min_release)  // Ensure we meet minimum release
            .min(specs.max_release)  // Cannot exceed maximum release
            .min(available_water);  // Cannot release more than available
        wc_storage[t] = wc_state.storage;
        wc_release[t] = wc_state.release;
        wc_spill[t] = wc_state.spill;

        // Calculate new storage
        let potential_storage = available_water - actual_release[t];
        // Cordell Hull simulation
        let ch_inflow = wc_state.release + wc_state.spill;
        let ch_state = simulate_timestep(
            &cordell_hull,
            ch_inflow,
            release_target[t],  // Note: might want separate release targets later
            if t == 0 { cordell_hull.initial_storage } else { ch_storage[t-1] },
        );

        // Calculate spill if storage would exceed capacity
        if potential_storage > specs.capacity {
            storage[t] = specs.capacity;
            spill[t] = potential_storage - specs.capacity;
        } else {
            storage[t] = potential_storage;
            spill[t] = 0.0;
        }
        ch_storage[t] = ch_state.storage;
        ch_release[t] = ch_state.release;
        ch_spill[t] = ch_state.spill;
    }

    Ok((storage, actual_release, spill))
    Ok((wc_storage, wc_release, wc_spill, ch_storage, ch_release, ch_spill))

}