Skip to content

Commit

Permalink
Fix two basin states issue
Browse files Browse the repository at this point in the history
  • Loading branch information
luca-tomasini committed Dec 12, 2024
1 parent d76f39a commit 672ffc2
Show file tree
Hide file tree
Showing 15 changed files with 41,936 additions and 4,621 deletions.
6 changes: 3 additions & 3 deletions docs/first_stage_baseline.rst
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ Sets
----

.. csv-table::
:file: /tables/sets_docs.csv
:file: /tables/first_stage/sets_docs.csv
:header-rows: 1

.. automodule:: pyomo_models.baseline.first_stage.sets
Expand All @@ -18,7 +18,7 @@ Variables
:no-index:

.. csv-table::
:file: /tables/variables_docs.csv
:file: /tables/first_stage/variables_docs.csv
:header-rows: 1

Parameters
Expand All @@ -28,7 +28,7 @@ Parameters
:no-index:

.. csv-table::
:file: /tables/parameters_docs.csv
:file: /tables/first_stage/parameters_docs.csv
:header-rows: 1

Objective
Expand Down
5 changes: 3 additions & 2 deletions docs/tables/baseline/second_stage/sets_docs.csv
Original file line number Diff line number Diff line change
Expand Up @@ -4,5 +4,6 @@
:math:`B`,1,,Water basins
:math:`S_B`,2,:math:`B`,State of each basins
:math:`S_H`,2,:math:`H`,State of each hydro powerplants (related to upstream basin state)
:math:`F`,2,:math:`H`,State of each hydro powerplant water flow
:math:`S_{BH}`,4,:math:`\{ H; ~B; ~S_H \;~S_B \}`,"Index to make the correspondence between the states, basins and hydro powerplants"
:math:`S_Q`,3,:math:`S_H`,State of each hydro powerplant water flow
:math:`B_{H}`,2,:math:`H`,Index to make the correspondence between hydro powerplants ant their upper stream bassin
:math:`SB_{H}`,3,:math:`S_H`,Index to make the correspondence between basins' and hydro powerplants' states
20 changes: 10 additions & 10 deletions docs/tables/baseline/second_stage/variables_docs.csv
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
Name,Description,Type,Units
":math:`V_\text{BAS}^{t,~b}`",Basins actual water volume,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`V_\text{SPIL}^{t,~b}`",Spilled water volumes when basins are full,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`V_\text{BAS, END}^{b}`",,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`dV_\text{+}^{t,~b}`",,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`dV_\text{-}^{t,~b}`",,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`P^{t,~h}`",Electrical power produced (or consumed) by a hydropower plant,:math:`\mathbb{R}^{+}`,:math:`\mathrm{MW}`
":math:`Q^{t,~h}`",Amount of water that flows through a hydropower plant,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
":math:`V_\text{BAS, END}^{b}`",Basins end volume,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`dV_\text{+}^{t,~h}`",Positive difference bertween powered water volume sepoints and calculated,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`dV_\text{-}^{t,~h}`",Positive difference bertween powered water volume sepoints and calculated,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
":math:`P^{t,~h}`",Electrical power produced (or consumed) by a hydropower plant,:math:`\mathbb{R}`,:math:`\mathrm{MW}`
":math:`Q^{t,~h}`",Amount of water that flows through a hydropower plant,:math:`\mathbb{R}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
":math:`S_\text{BAS}^{t,~b,~s_b}`","Binary variable that indicates the basin's volume state derived from its volume: :math:`S_\text{BAS}^{t, b, s_b} = \begin{cases} 1 \text{ if } V_\text{BAS}^{t, b} \in \left[ V_\text{BAS, MIN}^{b,~s_b} ;V_\text{BAS, MAX}^{b,~s_b} \right] \\ 0 \text{ otherwise} \end{cases}`",:math:`\mathbb{B}`,
":math:`S_\text{FLOW}^{t,~h,~s_q}`",Binary variable that indicates the hydo power plant's flow state (for turbine and pump),:math:`\mathbb{B}`,
":math:`P_\text{CAL}^{t,~h,~s_q}`",Calculated electrical power produced (or consumed) by a hydropower plant for each flow state,:math:`\mathbb{R}^{+}`,:math:`\mathrm{MW}`
":math:`Q_\text{CAL}^{t,~h,~s_q}`",Calculated water that flows through a hydropower plant for each flow state,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
":math:`P_\text{S}^{t,~h,~s_q}`","Real electrical power produced (or consumed) by a hydropower plant for each flow state: :math:`P_\text{S}^{t,~h,~s_q} = \begin{cases} P_\text{CAL}^{t,~h,~s_q} \text{ if } S_\text{FLOW}^{t,~h,~s_q} = 1 \\ 0 \text{ otherwise} \end{cases}`",:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
":math:`Q_\text{S}^{t,~h,~s_q}`","Real amount of water that flows through a pump for each hydro powerplant state: :math:`Q_\text{S}^{t,~h,~s_q} = \begin{cases} Q_\text{CAL}^{t,~h,~s_q} \text{ if } S_\text{FLOW}^{t,~h,~s_q} = 1 \\ 0",:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
":math:`S_\text{FLOW}^{t,~h,~s_h,~s_q}`",Binary variable that indicates the hydo power plant's flow state (for turbine and pump),:math:`\mathbb{B}`,
":math:`P_\text{CAL}^{t,~h,~s_h,~s_q}`",Calculated electrical power produced (or consumed) by a hydropower plant for each flow state,:math:`\mathbb{R}`,:math:`\mathrm{MW}`
":math:`Q_\text{CAL}^{t,~h,~s_h,~s_q}`",Calculated water that flows through a hydropower plant for each flow state,:math:`\mathbb{R}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
":math:`P_\text{S}^{t,~h,~s_h,~s_q}`","Real electrical power produced (or consumed) by a hydropower plant for each flow state: :math:`P_\text{S}^{t,~h,~s_h,~s_q} = \begin{cases} P_\text{CAL}^{t,~h,~s_h,~s_q} \text{ if } S_\text{FLOW}^{t,~h,~s_h,~s_q} = 1 \\ 0 \text{ otherwise} \end{cases}`",:math:`\mathbb{R}`,:math:`\mathrm{MW}`
":math:`Q_\text{S}^{t,~h,~s_h,~s_q}`","Real amount of water that flows through a pump for each hydro powerplant state: :math:`Q_\text{S}^{t,~h,~s_h,~s_q} = \begin{cases} Q_\text{CAL}^{t,~h,~s_h,~s_q} \text{ if } S_\text{FLOW}^{t,~h,~s_h,~s_q} = 1 0 \text{ otherwise} \end{cases}`",:math:`\mathbb{R}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
120 changes: 72 additions & 48 deletions scripts/baseline_optimizazion.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
import json
from datetime import timedelta
import polars as pl
from typing import Union
from multiprocessing import Pool

from polars import col as c
import shutil
Expand All @@ -19,90 +21,112 @@
YEARS = [2020, 2021, 2022, 2023]

TURBINE_FACTORS = {0.7, 0.8, 0.9}
SIMULATION_SETTING = {
"1": {"quantile": 0, "buffer": 0.2, "powered_volume_enabled": True, "with_penalty": True},
"2": {"quantile": 0.15, "buffer": 0.3, "powered_volume_enabled": True, "with_penalty": True},
"3": {"quantile": 0.25, "buffer": 0.3, "powered_volume_enabled": False, "with_penalty": True},
# "4": {"quantile": 0, "buffer": 0.2, "powered_volume_enabled": True, "with_penalty": False},
}
SIMULATION_SETTING = [
{"quantile": 0, "buffer": 0.2, "powered_volume_enabled": True, "with_penalty": True},
{"quantile": 0.15, "buffer": 0.3, "powered_volume_enabled": True, "with_penalty": True},
{"quantile": 0.25, "buffer": 0.3, "powered_volume_enabled": False, "with_penalty": True},
]
REAL_TIMESTEP = timedelta(hours=1)
FIRST_STAGE_TIMESTEP = timedelta(days=1)
SECOND_STAGE_TIME_SIM = timedelta(days=4)
TIME_LIMIT = 20 # in seconds
VOLUME_FACTOR = 1e-6

output_file_names: dict[str, str] = json.load(open(settings.OUTPUT_FILE_NAMES))

log = generate_log(name=__name__)

def solve_second_stage_model(
second_stage: BaselineSecondStage, plot_name: str
) -> tuple[BaselineSecondStage, float]:

second_stage.solve_model()
income: float = round(second_stage.simulation_results["income"].sum()/1e6, 3)
fig = plot_second_stage_result(
simulation_results=second_stage.simulation_results, time_divider=7*24
)

fig.write_html(plot_name)
return second_stage, income

if __name__=="__main__":
baseline_folder = output_file_names["baseline"]
if os.path.exists(baseline_folder):
shutil.rmtree(baseline_folder)
for year in YEARS:
year_folder = f"{baseline_folder}/year_{year}"
plot_folder = f"{baseline_folder}/plots_year_{year}"

income_result_list: list[dict] = []
log_book_final: pl.DataFrame = pl.DataFrame()
for turbine_factor in TURBINE_FACTORS:
test_folder = f"{year_folder}/turbine_factor_{turbine_factor}"
plot_folder = f"{test_folder}/plots"
log_book_final: pl.DataFrame = pl.DataFrame()

income_result: dict = {}
income_result["turbine_factor"] = turbine_factor

build_non_existing_dirs(plot_folder)

sim_results: dict[str, pl.DataFrame] = dict()
sim_summary: list = []

baseline_input: BaseLineInput = BaseLineInput(
input_schema_file_name=output_file_names["duckdb_input"],
real_timestep=timedelta(hours=1),
real_timestep=REAL_TIMESTEP,
year=year,
hydro_power_mask = c("name").is_in(["Aegina hydro"]),
volume_factor=1e-6
volume_factor=VOLUME_FACTOR
)

baseline_first_stage: BaselineFirstStage = BaselineFirstStage(
baseline_input, timestep=timedelta(days=1), turbine_factor=turbine_factor
first_stage: BaselineFirstStage = BaselineFirstStage(
input_instance=baseline_input,
timestep=timedelta(days=1),
turbine_factor=turbine_factor
)

baseline_first_stage.solve_model()

sim_results["first_stage"] = baseline_first_stage.simulation_results

sim_summary.append(["first_stage", baseline_first_stage.simulation_results["income"].sum()])
value = round(baseline_first_stage.simulation_results["income"].sum()/1e6, 3)
print(f"{year} year, {turbine_factor} turbine_factor and first_stage : {value}")
first_stage.solve_model()

sim_results["first_stage"] = first_stage.simulation_results
income_result["first_stage"] = round(first_stage.simulation_results["income"].sum()/1e6, 3)

fig = plot_first_stage_result(
simulation_results=baseline_first_stage.simulation_results, time_divider=7
simulation_results=first_stage.simulation_results, time_divider=7
)

fig.write_html(f"{plot_folder}/first_stage.html")

fig.write_html(f"{plot_folder}/{turbine_factor}_turbine_factor_first_stage.html")

optimization_inputs: list[list[Union[BaselineSecondStage, str]]] = []
income_results: list[dict] = []

for name, sim_setting in SIMULATION_SETTING.items():
baseline_second_stage: BaselineSecondStage = BaselineSecondStage(
for model_nb, sim_setting in enumerate(SIMULATION_SETTING):
second_stage: BaselineSecondStage = BaselineSecondStage(
input_instance=baseline_input,
first_stage=baseline_first_stage,
first_stage=first_stage,
timestep=timedelta(days=4),
time_limit=20,
time_limit=TIME_LIMIT,
model_nb=model_nb,
**sim_setting
)
baseline_second_stage.solve_model()

sim_results[f"second_stage_{name}"] = baseline_second_stage.simulation_results
sim_summary.append([name, baseline_second_stage.simulation_results["income"].sum()])
fig_path = f"{plot_folder}/{turbine_factor}_turbine_factor_model_{model_nb}.html"
optimization_inputs.append([second_stage, fig_path])

value = round(baseline_second_stage.simulation_results["income"].sum()/1e6, 3)
print(f"Results for {year} year, {turbine_factor} turbine_factor and second stage {name}\n")
print(f"Total income: {value}")
log_book = baseline_second_stage.log_book
with Pool(processes=len(optimization_inputs)) as pool:
results = pool.starmap(solve_second_stage_model, optimization_inputs)
for result in results:
second_stage, income = result
name = f"second_stage_{second_stage.model_nb}"
sim_results[f"{turbine_factor}_turbine_factor_model_{name}"] = second_stage.simulation_results
income_result[f"model_{name}"] = income
log_book = second_stage.log_book
if not log_book.is_empty():
print(f"Non optimal solutions:\n{log_book.to_pandas().to_string()}")

log_book_final = pl.concat([
log_book_final,
log_book.with_columns(pl.lit(name).alias("sim_name"))
], how="diagonal_relaxed")

fig = plot_second_stage_result(
simulation_results=baseline_second_stage.simulation_results, time_divider=7*24
)

fig.write_html(f"{plot_folder}/second_stage_{name}.html")

sim_results["summary"] = pl.DataFrame(sim_summary, schema=["name", "income"], orient="row")
sim_results["log_book"] = log_book_final
income_result_list.append(income_result)

log.info(f"Income results for {year} year:\n{sim_results["income_result"]}")


dict_to_duckdb(data=sim_results, file_path= f"{test_folder}/optimization_results.duckdb")
sim_results["income_result"] = pl.from_dicts(income_result_list)
sim_results["log_book"] = log_book_final

dict_to_duckdb(data=sim_results, file_path= f"{baseline_folder}/{year}_year_results.duckdb")
36 changes: 20 additions & 16 deletions src/pyomo_models/baseline/second_stage/constraints/basin_volume.py
Original file line number Diff line number Diff line change
Expand Up @@ -120,36 +120,40 @@ def basin_state_constraint(model, t, b):
## Basin volume state constraints used to determine the state of each basin #######################################
###################################################################################################################

@model.Constraint(model.T, model.BS) # type: ignore
def basin_state_max_active_constraint(model, t, b, s_b):
@model.Constraint(model.T, model.HQS) # type: ignore
def basin_state_max_active_constraint(model, t, h, s_h, s_q):
b = model.B_H[h].first()
return (
model.basin_volume_by_state[t, b, s_b] <=
model.max_basin_volume[b, model.S_B[b].last()] * model.basin_state[t, b, s_b]
model.basin_volume_by_state[t, h, s_h, s_q] <=
model.max_basin_volume[b, model.S_B[b].last()] * model.flow_state[t, h, s_h, s_q]
)

@model.Constraint(model.T, model.BS) # type: ignore
def basin_state_max_inactive_constraint(model, t, b, s_b):
@model.Constraint(model.T, model.HQS) # type: ignore
def basin_state_max_inactive_constraint(model, t, h, s_h, s_q):
b = model.B_H[h].first()
return (
model.basin_volume_by_state[t, b, s_b] >=
model.basin_volume_by_state[t, h, s_h, s_q] >=
model.basin_volume[t, b] -
model.max_basin_volume[b, model.S_B[b].last()] * (1 - model.basin_state[t, b, s_b])
model.max_basin_volume[b, model.S_B[b].last()] * (1 - model.flow_state[t, h, s_h, s_q])
)

@model.Constraint(model.T, model.BS) # type: ignore
def basin_state_min_active_constraint(model, t, b, s_b):
@model.Constraint(model.T, model.HQS) # type: ignore
def basin_state_min_active_constraint(model, t, h, s_h, s_q):
b = model.B_H[h].first()
return (
model.basin_volume_by_state[t, b, s_b] >=
model.min_basin_volume[b, model.S_B[b].first()] * model.basin_state[t, b, s_b]
model.basin_volume_by_state[t, h, s_h, s_q] >=
model.min_basin_volume[b, model.S_B[b].first()] * model.flow_state[t, h, s_h, s_q]
)



@model.Constraint(model.T, model.BS) # type: ignore
def basin_state_min_inactive_constraint(model, t, b, s_b):
@model.Constraint(model.T, model.HQS) # type: ignore
def basin_state_min_inactive_constraint(model, t, h, s_h, s_q):
b = model.B_H[h].first()
return (
model.basin_volume_by_state[t, b, s_b] <=
model.basin_volume_by_state[t, h, s_h, s_q] <=
model.basin_volume[t, b] -
model.min_basin_volume[b, model.S_B[b].first()] * (1 - model.basin_state[t, b, s_b])
model.min_basin_volume[b, model.S_B[b].first()] * (1 - model.flow_state[t, h, s_h, s_q])
)

return model
Loading

0 comments on commit 672ffc2

Please sign in to comment.