rust input validation fix (faf965e2) · Commits · HydroCHiPPs / Powersheds

examples/Cumberland/optimization.ipynb

deleted100644 → 0

+0 −336

Original line number	Diff line number	Diff line
		%% 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
		import multiprocessing
		import copy
		```

		%% 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]]))

		power_sequence = list(7 * (6 * [decision_vars[0]] + 6 * [decision_vars[1]] + 6 * [decision_vars[2]] + 6 * [decision_vars[3]]))

		local_cascade_data = copy.deepcopy(cascade_data) # Make sure to import copy

		# overwrite power with the decision variables
		#cascade_data.reservoirs['WolfCreek'].target_power = 168 * [decision_vars]
		local_cascade_data.reservoirs['WolfCreek'].target_power = power_sequence
		results = powersheds.simulate_cascade(local_cascade_data)

		# compute deviations at Cordell Hull
		x = np.abs(
		np.array(results["CordellHull"]["actual_power"]) - np.array(results["CordellHull"]["target_power"])
		)

		# compute pool elevation deviation from target
		pool_deviation = 230 - results['WolfCreek']['pool_elevation'][167]
		if pool_deviation <= 0:
		penalty = 0
		else:
		penalty = 10000000


		return np.sum(x*x) + penalty

		# Initial guess for variables
		#initial_guess = 0

		# Define bounds as a simple constraint for the optimization (using l-bfgs-b)
		bounds = [(0, 300), (0, 300), (0, 300), (0, 300)]

		# Get the number of available CPU cores
		num_cores = multiprocessing.cpu_count()

		print(num_cores)
		```

		%% Output

		16

		%% Cell type:code id:781a09a1 tags:

		``` python
		result = differential_evolution(
		simulation_function,
		bounds=bounds # Required when using multiple workers
		)
		print(result)
		```

		%% Output

		message: Optimization terminated successfully.
		success: True
		fun: 57234.68577925186
		x: [ 1.819e+02 1.196e+02 1.109e-01 3.294e+01]
		nit: 127
		nfev: 8150
		population: [[ 1.819e+02 1.196e+02 1.109e-01 3.294e+01]
		[ 1.843e+02 1.173e+02 3.724e-01 3.296e+01]
		...
		[ 1.822e+02 1.170e+02 3.962e+00 3.048e+01]
		[ 1.818e+02 1.199e+02 6.714e-01 3.201e+01]]
		population_energies: [ 5.723e+04 5.756e+04 ... 5.861e+04 5.766e+04]

		%% Cell type:code id:abce1fe3 tags:

		``` python
		from joblib import parallel_backend
		```

		%% Cell type:code id:e6819fa1 tags:

		``` python
		with parallel_backend('loky', n_jobs=15): # Try a smaller number like 4
		result = differential_evolution(
		simulation_function,
		bounds=bounds,
		updating='deferred'
		)

		print(result)
		```

		%% Output

		message: Optimization terminated successfully.
		success: True
		fun: 61459.087113666545
		x: [ 2.127e+02 2.946e+01 8.771e+00 8.021e+01]
		nit: 94
		nfev: 5975
		population: [[ 2.127e+02 2.946e+01 8.771e+00 8.021e+01]
		[ 2.059e+02 2.285e+01 3.522e+00 8.870e+01]
		...
		[ 2.003e+02 3.427e+01 1.681e+01 8.632e+01]
		[ 2.065e+02 3.436e+01 1.021e+01 8.449e+01]]
		population_energies: [ 6.146e+04 6.214e+04 ... 6.190e+04 6.258e+04]

		%% Cell type:code id:324f4f8a tags:

		``` python
		result = differential_evolution(
		simulation_function,
		bounds=bounds,
		workers=4,
		updating='immediate'
		)

		print(result)
		```

		%% Cell type:code id:d6bd99c9 tags:

		``` python
		import matplotlib.pyplot as plt
		from mpl_toolkits.mplot3d import Axes3D
		```

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

		``` python
		x = np.arange(0, 100)
		y = np.arange(200, 300)
		xgrid, ygrid = np.meshgrid(x, y)
		xy = np.stack([xgrid, ygrid])
		```

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

		``` python

		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:9e11a63a tags:

		``` python
		# Create the figure and axes object

		def compute_Z(X, Y):

		power_sequence = 7 * (12 * [X] + 12 * [Y])

		# 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"])
		)

		# compute pool elevation deviation from target
		pool_deviation = 230 - results['WolfCreek']['pool_elevation'][167]
		if pool_deviation <= 0:
		penalty = 0
		else:
		penalty = 1000000

		return np.sum(x*x) + penalty
		```

		%% Cell type:code id:42a7a246 tags:

		``` python
		compute_Z(0,0)
		```

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

		``` python
		# Create the meshgrid

		num_points = 50
		x = np.linspace(0, 200, num_points)
		y = np.linspace(0, 200, num_points)
		X, Y = np.meshgrid(x, y)

		# Initialize a Z array for computations
		Z = np.zeros_like(X)

		# Calculate the Z values for the surface
		for i in range(X.shape[0]):
		for j in range(X.shape[1]):
		Z[i, j] = compute_Z(X[i, j], Y[i, j])

		# Plotting
		fig = plt.figure()
		ax = fig.add_subplot(111, projection='3d')
		ax.plot_surface(X, Y, Z, cmap='viridis') # you may choose a different colormap

		```

		%% Cell type:code id:e2f3990e tags:

		``` python
		```

examples/Cumberland/script.qmd

deleted100644 → 0

+0 −207

Original line number	Diff line number	Diff line
		---
		title: "Powersheds Cumberland Example"
		format: html
		---

		```{r setup}
		#\| echo: false
		library(reticulate)
		library(dplyr)
		library(lubridate)
		library(readr)
		library(ggplot2)
		library(patchwork)
		library(purrr)
		#library(httpgd)
		reticulate::use_virtualenv("./.venv", required = TRUE)
		setwd("examples/Cumberland/")
		```


		```{python}
		import pandas as pd
		import powersheds
		import yaml
		from dataclasses import dataclass
		from pathlib import Path
		```

		```{python}

		@dataclass
		class ReservoirData:
		object_type: str
		capacity: float
		initial_pool_elevation: float
		min_power_pool: float
		# storage-headwater elevation tables:
		set_storage: list[float]
		set_elevation: list[float]
		# head-power-flow table
		hpf_h: list[float]
		hpf_p: list[float]
		hpf_q: list[float]
		# static tailwater elevation
		tailwater_elevation: float
		max_release: float
		min_release: float
		catchment_inflow: list[float]
		# target_release: list[float]
		target_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():

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

		# Storage–elevation table lookup
		#se_table = pd.read_csv(f"storage_HWelevation_tables/{name}.csv")

		# Head-Power-Flow table lookup
		#hpf_path = f"head_power_flow_tables/{name}.parquet"
		#hpf_table = pd.read_parquet(hpf_path)

		reservoir_dict[name] = ReservoirData(
		# specifications from yaml config
		**specs,

		# catchment inflow and target power time series
		catchment_inflow=pd.read_csv(f"time_series/{name}.csv")['catchment_inflow'].tolist(),
		target_power=pd.read_csv(f"time_series/{name}.csv")['target_power'].tolist(),

		# Storage-Elevation tables for interpolation...
		# ... (each column of the storage elevation tables is entered into powersheds as a vector of values)
		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-Power-Flow tables for interpolation...
		# ... (each column of the head-power-flow table is entered into powersheds as a vector of values)
		hpf_h = pd.read_parquet(f"head_power_flow_tables/{name}")["H_m"].astype(float).tolist(),
		hpf_p = pd.read_parquet(f"head_power_flow_tables/{name}")["P_MW"].astype(float).tolist(),
		hpf_q = pd.read_parquet(f"head_power_flow_tables/{name}")["Q_cumecs"].astype(float).tolist(),
		)
		elif specs['object_type'] == 'river':
		river_dict[name] = RiverData(
		# specifications for river objects from yaml config
		**specs,

		# legacy flows (i.e., initial conditions in rivers between reservoirs in the system)
		legacy_flows=pd.read_csv(f"time_series/{name}.csv")['legacy_flow'].tolist()
		)
		elif specs['object_type'] == 'confluence':
		confluence_dict[name] = ConfluenceData(
		# specifications for confluence objects from yaml config
		**specs,
		)

		cascade_data = CascadeData(
		reservoirs=reservoir_dict,
		rivers=river_dict,
		confluences=confluence_dict
		)
		```

		```{python}

		# run the powersheds model!

		results = powersheds.simulate_cascade(
		cascade_data
		)

		```

		```{r}
		names(py$reservoir_dict) -> reservoirs_ordered

		reservoirs_ordered \|> map_dfr(function(reservoir){

		py$results[[reservoir]] -> reservoir_results

		tibble(
		reservoir = !!reservoir,
		inflow = reservoir_results$total_inflow,
		target_release = reservoir_results$target_release,
		storage = reservoir_results$storage,
		release_implemented = reservoir_results$release,
		spill = reservoir_results$spill,
		pool_elevation = reservoir_results$pool_elevation,
		head = reservoir_results$head,
		power_actual = reservoir_results$actual_power,
		power_target = reservoir_results$target_power
		) \|> mutate(hr_index = 1:n())

		}) \|> mutate(reservoir = factor(reservoir, levels = reservoirs_ordered)) ->
		ts_data_for_plot

		powersheds_plot_theme <- function(){
		theme_classic() +
		theme(
		axis.text = element_text(size = 18),
		axis.title = element_text(size = 18),
		legend.text = element_text(size = 18),
		legend.position = "top")
		}

		plot_cols <- length(reservoirs_ordered)

		ts_data_for_plot \|>
		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)) +
		geom_line() +
		powersheds_plot_theme() +
		facet_wrap(~reservoir, ncol = plot_cols) -> flows_plot

		ts_data_for_plot \|>
		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, scales = "free") -> elevation_plot

		ts_data_for_plot \|>
		#filter(!grepl("_", reservoir), reservoir != "Celina") \|>
		select(reservoir, hr_index, storage) \|>
		ggplot(aes(hr_index, storage)) +
		geom_line() +
		powersheds_plot_theme() +
		facet_wrap(~reservoir, ncol = plot_cols) -> storage_plot

		#h <- hgd()
		#h$start()

		flows_plot + elevation_plot + plot_layout(ncol = 1)

		```
		No newline at end of file

src/lib.rs

+87 −16

Original line number	Diff line number	Diff line
		@@ -82,6 +82,89 @@ struct CascadeData {
		confluences: HashMap<String, ConfluenceData>
		}

		fn validate_reservoir_data(name: &str, reservoir: &ReservoirData, n: usize) -> PyResult<()> {
		if reservoir.catchment_inflow.len() != n {
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': catchment_inflow must have length {}",
		name, n
		)));
		}

		if reservoir.target_power.len() != n {
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': target_power must have length {}",
		name, n
		)));
		}

		if reservoir.min_power_pool.len() != n
		\|\| reservoir.max_release.len() != n
		\|\| reservoir.min_release.len() != n
		{
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': constraint time series must have length {}",
		name, n
		)));
		}

		if !reservoir.max_operating_elevation.is_empty()
		&& reservoir.max_operating_elevation.len() != n
		{
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': max_operating_elevation must be empty or length {}",
		name, n
		)));
		}

		if reservoir.set_storage.is_empty() \|\| reservoir.set_elevation.is_empty() {
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': set_storage and set_elevation must be non-empty",
		name
		)));
		}

		if reservoir.set_storage.len() != reservoir.set_elevation.len() {
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': set_storage and set_elevation must have the same length",
		name
		)));
		}

		if reservoir.hpf_h.is_empty() \|\| reservoir.hpf_p.is_empty() \|\| reservoir.hpf_q.is_empty() {
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': hpf_h, hpf_p, and hpf_q must be non-empty",
		name
		)));
		}

		if reservoir.hpf_h.len() != reservoir.hpf_p.len()
		\|\| reservoir.hpf_h.len() != reservoir.hpf_q.len()
		{
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': hpf_h, hpf_p, and hpf_q must have the same length",
		name
		)));
		}

		let has_tw_outflow = !reservoir.set_tw_outflow.is_empty();
		let has_tw_elevation = !reservoir.set_tw_elevation.is_empty();
		if has_tw_outflow != has_tw_elevation {
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': set_tw_outflow and set_tw_elevation must both be empty or both be provided",
		name
		)));
		}

		if has_tw_outflow && reservoir.set_tw_outflow.len() != reservoir.set_tw_elevation.len() {
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': set_tw_outflow and set_tw_elevation must have the same length",
		name
		)));
		}

		Ok(())
		}

		// Define structure for reservoir state
		struct ReservoirState {
		storage: f64,
		@@ -409,23 +492,11 @@ fn simulate_cascade(
		.ok_or_else(\|\| PyValueError::new_err("No reservoirs provided"))?;
		let n = first_reservoir.catchment_inflow.len();

		// Validate constraint time series lengths
		// Validate reservoir input shapes before allocating outputs or entering
		// the timestep loop so malformed user inputs raise Python errors instead
		// of panicking inside Rust helpers.
		for (name, reservoir) in &cascade_data.reservoirs {
		if reservoir.min_power_pool.len() != n
		\|\| reservoir.max_release.len() != n
		\|\| reservoir.min_release.len() != n
		{
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': constraint time series must have length {}", name, n
		)));
		}
		if !reservoir.max_operating_elevation.is_empty()
		&& reservoir.max_operating_elevation.len() != n
		{
		return Err(PyValueError::new_err(format!(
		"Reservoir '{}': max_operating_elevation must be empty or length {}", name, n
		)));
		}
		validate_reservoir_data(name, reservoir, n)?;
		}

		let mut results = HashMap::new();

tests/test_cascade.py

+50 −0

File changed.

Preview size limit exceeded, changes collapsed.