Skip to content

Commit a823875

Browse files
phumthepphumthep
committed
Revise reoperation
1 parent 2c10eb9 commit a823875

File tree

7 files changed

+293
-54
lines changed

7 files changed

+293
-54
lines changed

‎src/pownet/builder/hydro.py‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,19 @@ def update_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
157157
hydro_capacity=self.inputs.daily_hydro_capacity,
158158
)
159159

160+
def update_daily_hydropower_capacity(
161+
self, step_k: int, new_capacity: dict[(str, int), float]
162+
) -> None:
163+
self.model.remove(self.c_hydro_limit_daily)
164+
self.c_hydro_limit_daily = nondispatch_constr.add_c_hydro_limit_daily_dict(
165+
model=self.model,
166+
phydro=self.phydro,
167+
step_k=step_k,
168+
sim_horizon=self.inputs.sim_horizon,
169+
hydro_units=self.inputs.daily_hydro_unit_node.keys(),
170+
hydro_capacity_dict=new_capacity,
171+
)
172+
160173
def get_variables(self) -> dict[str, gp.tupledict]:
161174
"""Get hydro unit variables.
162175

‎src/pownet/core/builder.py‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -172,3 +172,15 @@ def update(self, step_k: int, init_conds: dict[str, dict]) -> PowerSystemModel:
172172

173173
self.model.update()
174174
return PowerSystemModel(self.model)
175+
176+
def get_phydro(self) -> gp.tupledict:
177+
"""Get the hydro power variable from the model."""
178+
return self.hydro_builder.phydro
179+
180+
def update_daily_hydropower_capacity(
181+
self, step_k: int, new_capacity: dict[(str, int), float]
182+
) -> None:
183+
"""Update the daily hydro capacity in the model."""
184+
self.hydro_builder.update_daily_hydropower_capacity(step_k, new_capacity)
185+
self.model.update()
186+
return PowerSystemModel(self.model)

‎src/pownet/optim_model/constraints/nondispatch_constr.py‎

Lines changed: 51 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,8 +81,8 @@ def add_c_hydro_limit_daily(
8181
max_day = sim_horizon // 24
8282
for day in range(step_k, step_k + max_day):
8383
for hydro_unit in hydro_units:
84-
cname = f"hydro_limit_daily[{hydro_unit},{day}]"
8584
current_day = day - step_k + 1
85+
cname = f"hydro_limit_daily[{hydro_unit},{current_day}]"
8686
constraints[cname] = model.addConstr(
8787
gp.quicksum(
8888
phydro[hydro_unit, t]
@@ -92,3 +92,53 @@ def add_c_hydro_limit_daily(
9292
name=cname,
9393
)
9494
return constraints
95+
96+
97+
def add_c_hydro_limit_daily_dict(
98+
model: gp.Model,
99+
phydro: gp.tupledict,
100+
step_k: int,
101+
sim_horizon: int,
102+
hydro_units: list,
103+
hydro_capacity_dict: pd.DataFrame,
104+
) -> gp.tupledict:
105+
"""
106+
Add constraints to limit hydropower by the daily amount. The sum of dispatch variables
107+
for each hydro unit over a 24-hour period must be less than or equal to the daily capacity.
108+
109+
Args:
110+
model (gp.Model): The optimization model
111+
phydro (gp.tupledict): The power output of hydro units
112+
step_k (int): The current iteration
113+
sim_horizon (int): The simulation horizon
114+
hydro_units (list): The list of hydro units
115+
hydro_capacity_dict: The daily capacity of the hydro unit
116+
117+
Returns:
118+
gp.tupledict: The constraints for the daily hydro limit
119+
120+
Raises:
121+
ValueError: If the simulation horizon is not divisible by 24
122+
123+
"""
124+
# When formulating with daily hydropower, sim_horizon must be divisible
125+
# by 24 because the hydro_capacity is daily.
126+
if sim_horizon % 24 != 0:
127+
raise ValueError(
128+
"The simulation horizon must be divisible by 24 when using daily hydropower capacity."
129+
)
130+
constraints = gp.tupledict()
131+
max_day = sim_horizon // 24
132+
for day in range(step_k, step_k + max_day):
133+
for hydro_unit in hydro_units:
134+
current_day = day - step_k + 1
135+
cname = f"hydro_limit_daily[{hydro_unit},{current_day}]"
136+
constraints[cname] = model.addConstr(
137+
gp.quicksum(
138+
phydro[hydro_unit, t]
139+
for t in range(1 + (current_day - 1) * 24, current_day * 24 + 1)
140+
)
141+
<= hydro_capacity_dict[hydro_unit, day],
142+
name=cname,
143+
)
144+
return constraints

‎src/pownet/reservoir/coupler.py‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
from pownet.core import ModelBuilder
2+
from .manager import ReservoirManager
3+
4+
import gurobipy as gp
5+
6+
from ..data_utils import get_unit_hour_from_varname
7+
8+
9+
class PowerWaterCoupler:
10+
def __init__(
11+
self,
12+
model_builder: ModelBuilder,
13+
reservoir_manager: ReservoirManager,
14+
solver: str = "gurobi",
15+
mip_gap: float = 0.0001,
16+
timelimit: float = 600,
17+
log_to_console: bool = False,
18+
):
19+
self.model_builder = model_builder
20+
self.reservoir_manager = reservoir_manager
21+
22+
self.solver = solver
23+
self.mipgap = mip_gap
24+
self.timelimit = timelimit
25+
self.log_to_console = log_to_console
26+
27+
self.num_days_in_step = self.model_builder.inputs.sim_horizon // 24
28+
29+
self.reop_iter = []
30+
self.reop_opt_time = 0.0
31+
32+
def get_reop_opt_time(self):
33+
return self.reop_opt_time
34+
35+
def get_reop_iter(self):
36+
return self.reop_iter
37+
38+
def reoperate(
39+
self,
40+
step_k: int,
41+
):
42+
43+
# Assume optimization is rolling horizon of 24 hours
44+
days_in_step = range(step_k, step_k + self.num_days_in_step)
45+
46+
reop_converge = False
47+
reop_k = 0
48+
49+
while not reop_converge:
50+
print(f"\nReservoirs reoperation iteration {reop_k}")
51+
52+
# --- PowNet returns the hydropower dispatch in hourly resolution across the simulation horizon
53+
hydropower_dispatch = {
54+
(unit, day): 0
55+
for unit in self.reservoir_manager.simulation_order
56+
for day in days_in_step
57+
}
58+
for varname, var in self.model_builder.get_phydro().items():
59+
unit = varname[0]
60+
61+
if varname[1] % 24 == 0:
62+
current_day = varname[1] // 24 + step_k - 1
63+
else:
64+
current_day = varname[1] // 24 + step_k
65+
66+
hydropower_dispatch[unit, current_day] += var.X
67+
68+
# --- Reoperate the reservoirs
69+
proposed_capacity = self.reservoir_manager.reoperate(
70+
daily_dispatch=hydropower_dispatch,
71+
days_in_step=days_in_step,
72+
)
73+
74+
# --- Iterate the reoperation process
75+
# Compare the new hydropower capacity with the current dispatch
76+
max_deviation = {
77+
(unit, day): abs(
78+
proposed_capacity[unit, day] - hydropower_dispatch[unit, day]
79+
)
80+
for unit in self.reservoir_manager.simulation_order
81+
for day in days_in_step
82+
}
83+
84+
# Set the tolerance for convergence to 5%
85+
reop_tol = {
86+
(idx): 0.05 * hydropower_dispatch[unit, day]
87+
for idx in max_deviation
88+
for day in days_in_step
89+
}
90+
91+
if all(
92+
max_deviation[unit, day] <= reop_tol[unit, day]
93+
for unit in self.reservoir_manager.simulation_order
94+
for day in days_in_step
95+
):
96+
reop_converge = True
97+
print(
98+
f"PowNet: Day {step_k+1} - Reservoirs converged at iteration {reop_k}"
99+
)
100+
101+
print("Max deviation:", max_deviation)
102+
103+
if reop_k > 50:
104+
raise ValueError(
105+
"Reservoirs reoperation did not converge after 100 iterations"
106+
)
107+
108+
# To reoptimize PowNet with the new hydropower capacity,
109+
# update the builder class
110+
power_system_model = self.model_builder.update_daily_hydropower_capacity(
111+
step_k=step_k, new_capacity=proposed_capacity
112+
)
113+
power_system_model.optimize(
114+
solver=self.solver,
115+
mipgap=self.mipgap,
116+
timelimit=self.timelimit,
117+
log_to_console=self.log_to_console,
118+
)
119+
120+
# Keep track of optimization time oand reoperation iterations
121+
self.reop_opt_time += power_system_model.get_runtime()
122+
reop_k += 1
123+
124+
# Record the number of iterations after convergence
125+
self.reop_iter.append(reop_k)

‎src/pownet/reservoir/manager.py‎

Lines changed: 57 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,22 @@
1212
)
1313

1414

15+
def find_upstream_flow(
16+
reservoir: Reservoir, reservoirs: dict[str, Reservoir]
17+
) -> pd.Series:
18+
unit_name = reservoir.name
19+
total_upstream_flow = pd.Series(
20+
0, index=range(1, reservoir.sim_days + 1), name="upstream_flow"
21+
)
22+
for upstream_unit in reservoir.upstream_units:
23+
upstream_reservoir = reservoirs[upstream_unit]
24+
# Get the release and spill from the upstream reservoir
25+
total_upstream_flow += (
26+
upstream_reservoir.release + upstream_reservoir.spill
27+
) * upstream_reservoir.downstream_flow_fracs[unit_name]
28+
return total_upstream_flow
29+
30+
1531
class ReservoirManager:
1632
def __init__(self):
1733
self.reservoirs: dict[str, Reservoir] = {}
@@ -78,35 +94,21 @@ def load_reservoirs_from_csv(self, input_folder: str) -> None:
7894
#############################################################################
7995
# Process the network topology
8096
##############################################################################
81-
self.simulation_order = find_simulation_order(flow_paths)
82-
# Ensure that all reservoirs are in the simulation order // The two sets must be equal
83-
if set(self.simulation_order) != set(self.reservoirs.keys()):
84-
raise ValueError(
85-
"The simulation order does not include all reservoirs: "
86-
f"{set(self.simulation_order) - set(self.reservoirs.keys())}"
87-
)
97+
self.simulation_order = find_simulation_order(
98+
reservoir_names=self.reservoirs.keys(), flow_paths=flow_paths
99+
)
88100

89101
def simulate(self) -> None:
90102
"""Simulate the reservoir operations to get hydropower time series."""
91103
for unit_name in self.simulation_order:
92104
reservoir = self.reservoirs[unit_name]
93-
# Find the upstream flow
94-
total_upstream_flow = pd.Series(
95-
0, index=range(1, reservoir.sim_days + 1), name="upstream_flow"
96-
)
97-
98-
for upstream_unit in reservoir.upstream_units:
99-
upstream_reservoir = self.reservoirs[upstream_unit]
100-
# Get the release and spill from the upstream reservoir
101-
total_upstream_flow += (
102-
upstream_reservoir.release + upstream_reservoir.spill
103-
) * upstream_reservoir.downstream_flow_fracs[unit_name]
104-
105-
# Simulate
105+
total_upstream_flow = find_upstream_flow(reservoir, self.reservoirs)
106106
reservoir.set_upstream_flow(total_upstream_flow)
107107
reservoir.simulate()
108108

109-
def get_hydropower_ts(self) -> pd.DataFrame:
109+
def get_hydropower_ts(
110+
self, unit_node_mapping: dict[str, str] = None
111+
) -> pd.DataFrame:
110112
"""Get the hydropower time series for all reservoirs."""
111113
df = pd.DataFrame()
112114
for unit_name in self.simulation_order:
@@ -119,10 +121,40 @@ def get_hydropower_ts(self) -> pd.DataFrame:
119121
df = pd.concat([df, temp_df], axis=1)
120122

121123
df.index = range(1, len(df) + 1)
124+
# Create multi-level column index if unit_node_mapping is provided
125+
if unit_node_mapping:
126+
df.columns = pd.MultiIndex.from_tuples(
127+
[(unit, unit_node_mapping[unit]) for unit in df.columns]
128+
)
122129
return df
123130

124-
def write_hydropower_to_csv(self, output_filepath: str) -> None:
131+
def write_hydropower_to_csv(
132+
self, output_filepath: str, unit_node_mapping: dict[str, str] = None
133+
) -> None:
125134
"""Write the hydropower time series to CSV files."""
126-
os.makedirs(output_filepath, exist_ok=True)
127-
for reservoir in self.reservoirs:
128-
reservoir.hydropower.to_csv(output_filepath, index=False)
135+
hydropower_df = self.get_hydropower_ts(unit_node_mapping)
136+
hydropower_df.to_csv(output_filepath, index=False)
137+
138+
def reoperate(
139+
self, daily_dispatch: dict[(str, int), float], days_in_step: range
140+
) -> dict[str, float]:
141+
"""Reoperate the reservoirs based on the daily dispatch of the power system model.
142+
Note that we don't reoperate on the first day of the simulation period.
143+
"""
144+
proposed_capacity = {k: 0 for k in daily_dispatch.keys()}
145+
146+
for unit_name in self.simulation_order:
147+
for day in days_in_step:
148+
reservoir = self.reservoirs[unit_name]
149+
# Find the upstream flow
150+
total_upstream_flow = find_upstream_flow(reservoir, self.reservoirs)
151+
152+
# Simulate
153+
reservoir.set_upstream_flow(total_upstream_flow)
154+
proposed_capacity[unit_name, day] = reservoir.reoperate(
155+
day=day,
156+
daily_dispatch=daily_dispatch[unit_name, day],
157+
upstream_flow_t=total_upstream_flow.loc[day],
158+
)
159+
160+
return proposed_capacity

0 commit comments

Comments
 (0)