#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""This module contains utility functions related to the head iteration.
"""
import datetime
import logging
from typing import Tuple, List, Dict, Any
import numpy as np
import pandas as pd
import pyoptinterface as poi
from prepshot.utils import interpolate_z_by_q_or_s
from prepshot.utils import cartesian_product
[docs]def initialize_waterhead(
stations : List[str], year : List[int],
month : List[int], hour : List[int],
params : Dict[str, Any]
) -> Tuple[pd.DataFrame, pd.DataFrame]:
"""Initialize water head.
Parameters
----------
stations : List[str]
List of stations.
year : List[int]
List of years.
month : List[int]
List of months.
hour : List[int]
List of hours.
params : Dict[str, Any]
Dictionary of parameters for the model.
Returns
-------
Tuple[pd.DataFrame, pd.DataFrame]
A tuple of two pandas.DataFrame objects, the first one is
the old water head, the second one is the new water head.
"""
old_waterhead = pd.DataFrame(
index=stations,
columns=pd.MultiIndex.from_product(
[year, month, hour], names=['year', 'month', 'hour']
)
)
new_waterhead = old_waterhead.copy(deep=True)
for s in stations:
old_waterhead.loc[s, :] = [
params['reservoir_characteristics']['head', s]
] * (len(hour) * len(month) * len(year))
return old_waterhead, new_waterhead
[docs]def compute_error(
old_waterhead : pd.DataFrame, new_waterhead : pd.DataFrame
) -> float:
"""Calculate the error of the water head.
Parameters
----------
old_waterhead : pandas.DataFrame
The water head before the solution.
new_waterhead : pandas.DataFrame
The water head after the solution.
Returns
-------
float
The error of the water head.
"""
new_waterhead[new_waterhead <= 0] = 1
error = (
abs(new_waterhead - old_waterhead)
/ new_waterhead
).mean(axis='columns').mean()
return error
[docs]def process_model_solution(
model : object, stations : List[str], year : List[int], month : List[int],
hour : List[int], params : Dict[str, Any],
old_waterhead : pd.DataFrame, new_waterhead : pd.DataFrame
) -> bool:
"""Process the solution of the model, updating the water head data.
Parameters
----------
model : object
Model to be solved.
stations : list
List f hydropower stations.
year : list
List of years.
month : list
List of months.
hour : list
List of hours.
params : dict
Dictionary of parameters for the model.
old_waterhead : pandas.DataFrame
The water head before the solution.
new_waterhead : pandas.DataFrame
The water head after the solution.
Returns
-------
bool
True if the model is solved, False otherwise.
"""
idx = pd.IndexSlice
for s, h, m, y in cartesian_product(stations, hour, month, year):
efficiency = params['reservoir_characteristics']['coeff', s]
model.set_normalized_coefficient(
model.output_calc_cons[s, h, m, y],
model.genflow[s, h, m, y],
- efficiency * 1e-3 * old_waterhead.loc[s, idx[y, m, h]]
)
# Solve the model and check the solution status.
model.set_model_attribute(poi.ModelAttribute.Silent, False)
model.optimize() # add log into log file
status = model.get_model_attribute(poi.ModelAttribute.TerminationStatus)
if status != poi.TerminationStatusCode.OPTIMAL:
return False
if params['iteration_number'] <= 1:
# If fixed head is True, do not update water head.
new_waterhead = old_waterhead
return True
# Iterate over each station to update water head data.
for stcd in stations:
outflow = np.array([[
[model.get_value(model.outflow[int(stcd), h, m, y]) for h in hour]
for m in month] for y in year]
)
storage = np.array([[
[model.get_value(model.storage_reservoir[int(stcd), h, m, y])
for h in model.hour_p] for m in month] for y in year]
)
tail = interpolate_z_by_q_or_s(
str(stcd), outflow,
params['reservoir_tailrace_level_discharge_function']
)
storage = interpolate_z_by_q_or_s(
str(stcd), storage, params['reservoir_forebay_level_volume_function']
)
# Calculate the new water head.
fore = (storage[:, :, :hour[-1]] + storage[:, :, 1:]) / 2
h = np.maximum(fore - tail, 0)
new_waterhead.loc[int(stcd), :] = h.ravel()
return True
[docs]def run_model_iteration(
model : object, params : Dict[str, Any],
error_threshold=0.001, max_iterations=5
) -> bool:
"""Run the model iteratively.
Parameters
----------
model : object
Model to be solved.
params : dict
Dictionary of parameters for the model.
error_threshold : float, optional
The error threshold, by default 0.001
max_iterations : int, optional
The maximum number of iterations, by default 5
Returns
-------
bool
True if the model is solved, False otherwise.
"""
logging.info(
'Starting iteration recorded at %s.',
datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
)
# Initialize water head.
stations, years, months, hours = \
params['stcd'], params['year'], params['month'], params['hour']
old_waterhead, new_waterhead = initialize_waterhead(
stations, years, months, hours, params
)
# Variables for iteration.
error = 1
errors = []
# Perform water head iteration.
for iteration in range(1, max_iterations+1):
alpha = 1 / iteration
success = process_model_solution(
model, stations, years, months, hours, params,
old_waterhead, new_waterhead
)
if not success:
return False
# Calculate error.
if max_iterations <= 1:
logging.warning(
"Maximum iteration is set to 1 and "
+ "the model will be solved with fixed head."
)
error = 0
else:
error = compute_error(old_waterhead, new_waterhead)
errors.append(error)
logging.info('Water head error: %.2f%%', error)
if error < error_threshold:
return True
# Update old water head for next iteration.
old_waterhead += alpha * (new_waterhead - old_waterhead)
logging.warning(
"Ending iteration recorded at %s."
+ "Failed to converge. Maximum iteration exceeded.",
datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
)
return True