optimization tests

  - pandas>=2.2,<2.3

  - numpy>=2.2,<2.3

  - pip>=24.3,<24.4

  - matplotlib>=3.10

  - pip:

    - maturin>=1.8.1,<1.9

    - scipy>=1.15

%% Cell type:markdown id:514e7b0a tags:

## Example optimization

This script illustrates how to optimize the power generating schedule at Wolf Creek dam in order to minimize downstream violations in meeting the power target and without dropping below a specified target storage at Wolf Creek. There are 168 decision variables (release from Wolf Creek) and the objective function is

%% Cell type:code id:1ed54f38 tags:

``` python

# import required libraries

import pandas as pd

import powersheds

import yaml

from dataclasses import dataclass

from pathlib import Path

```

%% Cell type:code id:7a67ff6f tags:

``` python

# set up the model input

@dataclass

class ReservoirData:

    object_type: str

    capacity: float

    initial_storage: float

    set_storage: list[float]

    set_elevation: list[float]

    tailwater_elevation: float

    max_release: float

    min_release: float

    catchment_inflow: list[float]

    target_release: list[float]

    target_power: list[float]

    hrpt_head: list[float]

    hrpt_release: list[float]

    hrpt_power: list[float]

    simulation_order: int

    downstream_object: str

@dataclass

class RiverData:

    object_type: str

    simulation_order: int

    downstream_object: str

    lag: int

    legacy_flows: list[float]

@dataclass

class ConfluenceData:

    object_type: str

    simulation_order: int

    downstream_object: str

@dataclass

class CascadeData:

    reservoirs: dict[str, ReservoirData]

    rivers: dict[str, RiverData]

    confluences: dict[str, ConfluenceData]

# Read the reservoir configuration

with open('cascade_config.yaml', 'r') as file:

    config_dict = yaml.safe_load(file)

# Initialize dictionaries for reservoirs and rivers

reservoir_dict = {}

river_dict = {}

confluence_dict = {}

for name, specs in config_dict.items():

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

    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(),

            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(

            **specs,

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

        )

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

        confluence_dict[name] = ConfluenceData(

            **specs,

        )

cascade_data = CascadeData(

    reservoirs=reservoir_dict,

    rivers=river_dict,

    confluences=confluence_dict

)

```

%% Cell type:code id:083d3c27 tags:

``` python

# check that the model runs

results = powersheds.simulate_cascade(cascade_data)

```

%% Cell type:code id:ace3af70 tags:

``` python

# import tools for optimization

import numpy as np

from scipy.optimize import basinhopping, differential_evolution, dual_annealing

```

%% Cell type:code id:1e4bde66 tags:

``` python

decision_vars = 0

168 *  [decision_vars]

```

%% Cell type:code id:ac8675b5 tags:

``` python

# set up the function to optimize.

# Define the same simulation function as above

def simulation_function(decision_vars):

    power_sequence = list(7 * (12 * [decision_vars[0]] + 12 * [decision_vars[1]]))

    # overwrite power with the decision variables

    #cascade_data.reservoirs['WolfCreek'].target_power = 168 * [decision_vars]

    cascade_data.reservoirs['WolfCreek'].target_power = power_sequence

    results = powersheds.simulate_cascade(cascade_data)

    # compute deviations at Old Hickory

    x = np.abs(

    np.array(results["CordellHull"]["actual_power"]) - np.array(results["CordellHull"]["target_power"])

    )

    return np.sum(x*x)

# Initial guess for variables

#initial_guess = 0

# Define bounds as a simple constraint for the optimization (using l-bfgs-b)

bounds = [(0, 50), (100, 300)]

# Run the basinhopping optimization

#result = basinhopping(simulation_function, initial_guess, niter=100, T=1.0, stepsize=0.5, minimizer_kwargs={"bounds": bounds})

result = differential_evolution(simulation_function, bounds=bounds)

print(result)

```

%% Cell type:code id:37a6e404 tags:

``` python

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

x = np.arange(0, 100)

y = np.arange(200, 300)

xgrid, ygrid = np.meshgrid(x, y)

xy = np.stack([xgrid, ygrid])

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

ax.view_init(45, -45)

ax.plot_surface(xgrid, ygrid, simulation_function(xy), cmap='terrain')

ax.set_xlabel('x')

ax.set_ylabel('y')

ax.set_zlabel('eggholder(x, y)')

plt.show()

```

%% Cell type:code id:4bbed8d3 tags:

``` python

cascade_data.reservoirs['WolfCreek'].target_power = 168 * [decision_vars]

cascade_data.reservoirs['WolfCreek'].target_power

```