-
Notifications
You must be signed in to change notification settings - Fork 0
/
pwl_opt.py
215 lines (177 loc) · 9.64 KB
/
pwl_opt.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
"""
An optimization module based on piecewise linear formulation
Not in use in the current study and may contain errors
"""
import numpy as np
import rsome
from rsome import ro
from rsome import grb_solver as grb
import graphs
import utils
from pds import PDS
from wds import WDSPWL
class Opt:
def __init__(self, pds_data: str, wds_data: str, t: int, n: int):
self.pds_data = pds_data
self.wds_data = wds_data
self.t = t
self.n = n # number of breakpoints for piecewise linearization
self.pds, self.wds = self.init_distribution_systems()
self.model = ro.Model()
self.x = self.declare_vars()
self.pl_flow_mat = self.build_piecewise_matrices('flow')
self.pl_head_mat = self.build_piecewise_matrices('head')
self.pl_power_mat = self.build_piecewise_matrices('power')
self.build_opt_problem()
def init_distribution_systems(self):
pds = PDS(self.pds_data)
wds = WDSPWL(self.wds_data, self.n)
return pds, wds
def declare_vars(self):
gen_p = self.model.dvar((self.pds.n_bus, self.t)) # active power generation
gen_q = self.model.dvar((self.pds.n_bus, self.t)) # reactive power generation
psh_y = self.model.dvar((self.pds.n_psh, self.t)) # psh_y = pumped storage hydropower injection
psh_h = self.model.dvar((self.pds.n_psh, self.t)) # psh_y = pumped storage hydropower consumption
v = self.model.dvar((self.pds.n_bus, self.t)) # buses squared voltage
I = self.model.dvar((self.pds.n_lines, self.t)) # lines squared current flow
p = self.model.dvar((self.pds.n_lines, self.t)) # lines active power flow
q = self.model.dvar((self.pds.n_lines, self.t)) # lines reactive power flow
#
self.model.st(gen_p >= 0)
self.model.st(gen_q >= 0)
self.model.st(psh_y >= 0)
self.model.st(psh_h >= 0)
self.model.st(I >= 0)
alpha = self.model.dvar((self.wds.n_pipes, self.t, self.n)) # pipe segments relative flows
beta = self.model.dvar((self.wds.n_pipes, self.t, self.n - 1), vtype='B')
h = self.model.dvar((self.wds.n_nodes, self.t)) # nodes head
pf = self.model.dvar((self.wds.n_pumps, self.t)) # pump flows
vol = self.model.dvar((self.wds.n_tanks, self.t)) # tanks volume
self.model.st(alpha <= 1)
self.model.st(alpha >= 0)
self.model.st(alpha.sum(axis=-1) == 1)
self.model.st(beta.sum(axis=-1) == 1)
# piecewise linear constraints
mat = np.zeros((self.n - 1, self.n))
np.fill_diagonal(mat, 1)
np.fill_diagonal(mat[:, 1:], 1)
self.model.st(alpha - beta @ mat <= 0)
return {'gen_p': gen_p, 'gen_q': gen_q, 'psh_y': psh_y, 'psh_h': psh_h, 'v': v, 'I': I, 'p': p, 'q': q,
'vol': vol, 'alpha': alpha, 'h': h, 'pf': pf, 'beta': beta}
def build_opt_problem(self):
self.objective_func()
self.bus_balance()
self.energy_conservation()
self.voltage_bounds()
self.power_flow_constraint()
self.water_mass_balance()
self.no_pumps_backflow()
self.head_boundaries()
self.head_conservation()
def build_piecewise_matrices(self, param):
""" Build matrices of the x, and y values at the piecewise linear breakpoints """
mat = np.zeros((self.wds.n_pipes, self.t, self.n))
for p in self.wds.pipes.index:
for point in range(self.n):
mat[p, :, point] = self.wds.pipes_pl[p][point][param]
return mat
def objective_func(self):
pds_cost = (self.pds.gen_mat @ (self.pds.pu_to_kw * self.x['gen_p']) @ self.pds.tariffs.values).sum()
pumps = self.pumps_power().sum(axis=0)
if self.wds.n_turbines > 0:
turbine = self.turbine_power().sum(axis=0)
else:
turbine = 0
wds_cost = (pumps - 0.9 * turbine) @ self.wds.tariffs.sum(axis=1).values
# psh_cost = (self.pds.psh['fill_tariff'].values @ self.x['psh_y']).sum()
self.model.min(pds_cost + wds_cost)
def bus_balance(self):
r = utils.connectivity_mat(self.pds.lines, from_col='from_bus', to_col='to_bus', direction='in',
param='r_pu')
x = utils.connectivity_mat(self.pds.lines, from_col='from_bus', to_col='to_bus', direction='in',
param='x_pu')
a = utils.connectivity_mat(self.pds.lines, from_col='from_bus', to_col='to_bus')
self.model.st(self.pds.gen_mat @ self.x['gen_p'] + a @ self.x['p']
- r @ self.x['I']
- self.pds.dem_active.values
+ self.pds.bus.loc[:, 'G'].values @ self.x['v']
== 0)
self.model.st(self.pds.gen_mat @ self.x['gen_q'] + a @ self.x['q']
- x @ self.x['I']
- self.pds.dem_reactive.values
+ self.pds.bus.loc[:, 'B'].values @ self.x['v']
== 0)
def energy_conservation(self):
r = self.pds.lines['r_pu'].values.reshape(1, -1)
x = self.pds.lines['x_pu'].values.reshape(1, -1)
a = utils.connectivity_mat(self.pds.lines, from_col='from_bus', to_col='to_bus')
self.model.st(a.T @ self.x['v']
+ 2 * ((self.x['p'].T * r).T + (self.x['q'].T * x).T)
- (self.x['I'].T * (r ** 2 + x ** 2)).T
== 0)
def voltage_bounds(self):
self.model.st(self.x['v'] - self.pds.bus['Vmax_pu'].values.reshape(-1, 1) <= 0)
self.model.st(self.pds.bus['Vmin_pu'].values.reshape(-1, 1) - self.x['v'] <= 0)
def power_flow_constraint(self):
for t in range(self.t):
for line in range(self.pds.n_lines):
b_id = self.pds.lines.loc[line, 'to_bus']
self.model.st(rsome.rsocone(self.x['p'][line, t] + self.x['q'][line, t],
self.x['v'][b_id, t],
self.x['I'][line, t]))
def water_mass_balance(self):
not_source = utils.get_mat_for_type(self.wds.nodes, 'reservoir', inverse=True) # all but reservoirs
a = utils.connectivity_mat(self.wds.pipes, from_col='from_node', to_col='to_node')
tanks_mat = utils.get_mat_for_type(self.wds.nodes, 'tank')
# rows - all nodes. columns - only tanks
# tanks_mat @ x_volume results in a matrix with the desired shape
tanks_mat = tanks_mat[:, ~np.all(tanks_mat == 0, axis=0)]
dt = utils.get_dt_mat(self.t)
init_vol = np.zeros((self.wds.n_tanks, self.t))
init_vol[:, 0] = self.wds.tanks['init_vol'].values
init_vol = tanks_mat @ init_vol
self.model.st(not_source @ a @ ((self.pl_flow_mat * self.x['alpha']).sum(axis=-1))
- ((tanks_mat @ self.x['vol']) @ dt) + init_vol
- self.wds.demands.values == 0)
self.model.st(self.x['vol'] <= self.wds.tanks['max_vol'].values.reshape(-1, 1))
self.model.st(self.x['vol'] >= self.wds.tanks['min_vol'].values.reshape(-1, 1))
def no_pumps_backflow(self):
pumps_idx = self.wds.pipes.loc[self.wds.pipes['type'] == 'pump'].index
self.model.st(self.x['alpha'][pumps_idx, :, :] >= 0)
def head_boundaries(self):
self.model.st(self.x['h'] >= self.wds.nodes['elevation'].values.reshape(-1, 1))
tanks_idx = self.wds.nodes.loc[self.wds.nodes['type'] == 'tank'].index.values
tanks_diam = self.wds.tanks.loc[tanks_idx, 'diameter'].values.reshape(-1, 1)
self.model.st((self.x['h'])[tanks_idx, :] - (self.x['vol']) * (4 / (np.pi * tanks_diam ** 2)) >= 0)
res_idx = self.wds.nodes.loc[self.wds.nodes['type'] == 'reservoir'].index
self.model.st(self.x['h'][res_idx, :] == self.wds.nodes.loc[res_idx, 'elevation'].values.reshape(-1, 1))
def head_conservation(self):
"""
numpy allows for "broadcasting" matrix multiplication, to avoid loops
This will work when the n columns in the first matrix is equal to the n rows in the second matrix.
numpy has a function called tensordot that can perform the dot product over specified axes for arrays of
any dimensionality.
"""
a = utils.connectivity_mat(self.wds.pipes, from_col='from_node', to_col='to_node')
# exclude turbines from headloss constraint
b = utils.get_mat_for_type(data=self.wds.pipes, element_type='turbine', inverse=True)
bb = np.tensordot(b, self.pl_head_mat, axes=([1], [0]))
self.model.st((a @ b).T @ self.x['h'] - ((bb * self.x['alpha']).sum(axis=-1)) == 0)
def pumps_power(self):
idx = self.wds.pipes.loc[self.wds.pipes['type'] == 'pump'].index
pl_power = np.take(self.pl_power_mat, idx, axis=0)
return (self.pl_power_mat * self.x['alpha']).sum(axis=-1)
def turbine_power(self):
idx = self.wds.pipes.loc[self.wds.pipes['type'] == 'turbine'].index
return (self.pl_power_mat[idx, :, :] * self.x['alpha'][idx, :, :]).sum(axis=-1)
def solve(self):
self.model.solve(grb, display=True, params={'TimeLimit': 45, 'MIPGap': 0.01})
obj, status = self.model.solution.objval, self.model.solution.status
print(obj, status)
def extract_res(self, x, elem_type, series_type='time', elem_idx=None, t_idx=None, dec=1, factor=1):
n = {'nodes': self.pds.n_bus, 'lines': self.pds.n_lines}
if series_type == 'time':
values = {t: self.x[x].get()[elem_idx, t] * factor for t in range(self.t)}
if series_type == 'elements':
values = {i: self.x[x].get()[i, t_idx] * factor for i in range(n[elem_type])}
return {t: round(val, dec) for t, val in values.items()}