-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathopf_lag.mod
294 lines (254 loc) · 10.4 KB
/
opf_lag.mod
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
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
reset;
param baseMVA;
set bus_i;
param type {bus_i};
param Pd {bus_i};
param Qd {bus_i};
param Gs {bus_i};
param Bs {bus_i};
param area {bus_i};
param Vm {bus_i};
param Va {bus_i};
param baseKV{bus_i};
param zone {bus_i};
param Vmax {bus_i};
param Vmin {bus_i};
set gen_i within bus_i;
param Pg {gen_i};
param Qg {gen_i};
param Qmax {gen_i};
param Qmin {gen_i};
param Vg {gen_i};
param mBase {gen_i};
param gen_status{gen_i};
param Pmax {gen_i};
param Pmin {gen_i};
set branch dimen 2;
param r {branch};
param x {branch};
param b {branch};
param rateA {branch};
param rateB {branch};
param rateC {branch};
param ratio {branch};
param angle {branch};
param branch_status{branch};
param startup {gen_i};
param shutdown{gen_i};
param n {gen_i};
param c0 {gen_i};
param c1 {gen_i};
param c2 {gen_i};
###
# other parameter
###
param line_G{(l,m) in branch} := r[l,m]/(r[l,m]^2+x[l,m]^2);
param line_B{(l,m) in branch} := -x[l,m]/(r[l,m]^2+x[l,m]^2);
param pst_ratio{(l,m) in branch} := 1;
param pst_angle{(l,m) in branch} := 0;
param pst_cos{(l,m) in branch} := cos(pst_angle[l,m]);
param pst_sin{(l,m) in branch} := sin(pst_angle[l,m]);
param y11_re{(l,m) in branch}:= line_G[l,m]/(pst_ratio[l,m]^2);
param y11_im{(l,m) in branch}:= (line_B[l,m]+b[l,m]/2)/(pst_ratio[l,m]^2);
param y22_re{(l,m) in branch}:= line_G[l,m];
param y22_im{(l,m) in branch}:= line_B[l,m]+b[l,m]/2;
# -ys/(rho exp +jTheta)
param y12_re{(l,m) in branch}:= (-1/pst_ratio[l,m])*( line_G[l,m]*pst_cos[l,m]-line_B[l,m]*pst_sin[l,m]);
param y12_im{(l,m) in branch}:= (-1/pst_ratio[l,m])*( line_B[l,m]*pst_cos[l,m]+line_G[l,m]*pst_sin[l,m]);
# -ys/(rho exp -jTheta)
param y21_re{(l,m) in branch}:= (-1/pst_ratio[l,m])*( line_G[l,m]*pst_cos[l,m]+line_B[l,m]*pst_sin[l,m]);
param y21_im{(l,m) in branch}:= (-1/pst_ratio[l,m])*( line_B[l,m]*pst_cos[l,m]-line_G[l,m]*pst_sin[l,m]);
###
###
# OPF formulation
###
param dual_p{bus in bus_i} default 0;
param dual_q{bus in bus_i} default 0;
#param dual_vmin{bus in bus_i} default 0;
#param dual_vmax{bus in bus_i} default 0;
param dual_convexity default 0;
param nIte default 0;
param V_re{ite in 1..nIte, bus in bus_i};
param V_im{ite in 1..nIte, bus in bus_i};
param P_or{ite in 1..nIte, (bus1,bus2) in branch} = 100*(
+y11_re[bus1, bus2]*V_re[ite,bus1]*V_re[ite,bus1]
+y11_re[bus1, bus2]*V_im[ite,bus1]*V_im[ite,bus1]
+y12_re[bus1, bus2]*V_re[ite,bus1]*V_re[ite,bus2]
+y12_re[bus1, bus2]*V_im[ite,bus1]*V_im[ite,bus2]
+y12_im[bus1, bus2]*V_im[ite,bus1]*V_re[ite,bus2]
-y12_im[bus1, bus2]*V_re[ite,bus1]*V_im[ite,bus2]);
param P_ex{ite in 1..nIte, (bus1,bus2) in branch} = 100*(
+y22_re[bus1,bus2]*V_re[ite,bus2]*V_re[ite,bus2]
+y22_re[bus1,bus2]*V_im[ite,bus2]*V_im[ite,bus2]
+y21_re[bus1,bus2]*V_re[ite,bus1]*V_re[ite,bus2]
+y21_re[bus1,bus2]*V_im[ite,bus1]*V_im[ite,bus2]
+y21_im[bus1,bus2]*V_re[ite,bus1]*V_im[ite,bus2]
-y21_im[bus1,bus2]*V_im[ite,bus1]*V_re[ite,bus2]);
param Q_or{ite in 1..nIte, (bus1,bus2) in branch} = 100*(
-y11_im[bus1,bus2]*V_re[ite,bus1]*V_re[ite,bus1]
-y11_im[bus1,bus2]*V_im[ite,bus1]*V_im[ite,bus1]
+y12_re[bus1,bus2]*V_im[ite,bus1]*V_re[ite,bus2]
-y12_re[bus1,bus2]*V_re[ite,bus1]*V_im[ite,bus2]
-y12_im[bus1,bus2]*V_re[ite,bus1]*V_re[ite,bus2]
-y12_im[bus1,bus2]*V_im[ite,bus1]*V_im[ite,bus2]);
param Q_ex{ite in 1..nIte, (bus1,bus2) in branch} = 100*(
-y22_im[bus1,bus2]*V_re[ite,bus2]*V_re[ite,bus2]
-y22_im[bus1,bus2]*V_im[ite,bus2]*V_im[ite,bus2]
+y21_re[bus1,bus2]*V_re[ite,bus1]*V_im[ite,bus2]
-y21_re[bus1,bus2]*V_im[ite,bus1]*V_re[ite,bus2]
-y21_im[bus1,bus2]*V_re[ite,bus1]*V_re[ite,bus2]
-y21_im[bus1,bus2]*V_im[ite,bus1]*V_im[ite,bus2]);
set master_basis default {};
###
# master
###
problem master;
param FEASIBILITY default 1;
var p_gen{gen in gen_i};
subject to p_gen_min{gen in gen_i}: p_gen[gen] >= Pmin[gen]*100;
subject to p_gen_max{gen in gen_i}: p_gen[gen] <= Pmax[gen]*100;
var q_gen{gen in gen_i};
subject to q_gen_min{gen in gen_i}: q_gen[gen] >= Qmin[gen]*100;
subject to q_gen_max{gen in gen_i}: q_gen[gen] <= Qmax[gen]*100;
var lambda{ite in 1..nIte} >= 0;
var fake_lambda >= 0;
var p_pos{bus in bus_i} >= 0;
var p_neg{bus in bus_i} >= 0;
var q_pos{bus in bus_i} >= 0;
var q_neg{bus in bus_i} >= 0;
#var v_min_slack{bus in bus_i} >= 0;
#var v_max_slack{bus in bus_i} >= 0;
subject to convexity:
+(if FEASIBILITY == 1 then fake_lambda)+sum{ite in 1..nIte} lambda[ite]=1;
subject to p_balance{bus in bus_i}:
+(if FEASIBILITY == 1 then +p_pos[bus]-p_neg[bus])
-(if bus in gen_i then p_gen[bus])/100
+sum{ite in 1..nIte, (bus1, bus2) in branch:bus==bus1} P_or[ite,bus1,bus2]*lambda[ite]
+sum{ite in 1..nIte, (bus1, bus2) in branch:bus==bus2} P_ex[ite,bus1,bus2]*lambda[ite]
=
-Pd[bus]
;
subject to q_balance{bus in bus_i}:
+(if FEASIBILITY == 1 then +q_pos[bus]-q_neg[bus])
-(if bus in gen_i then q_gen[bus])/100
+sum{ite in 1..nIte, (bus1, bus2) in branch:bus==bus1} Q_or[ite,bus1,bus2]*lambda[ite]
+sum{ite in 1..nIte, (bus1, bus2) in branch:bus==bus2} Q_ex[ite,bus1,bus2]*lambda[ite]
=
-Qd[bus]
;
minimize master_obj:
+if FEASIBILITY == 1 then(
+10e5*fake_lambda
+10e5*sum{bus in bus_i} (
+p_pos[bus]+p_neg[bus]
+q_pos[bus]+q_neg[bus]
#+v_min_slack[bus]+v_max_slack[bus]
)
)else
+sum{gen in gen_i}(c0[gen]+c1[gen]*p_gen[gen]+c2[gen]*p_gen[gen]^2)
;
#subject to ctr_v_min{bus in bus_i}: sum{ite in 1..nIte}(V_re[ite,bus]^2+V_im[ite,bus]^2)*lambda[ite] + v_min_slack[bus] >= Vmin[bus]^2 ;
#subject to ctr_v_max{bus in bus_i}: sum{ite in 1..nIte}(V_re[ite,bus]^2+V_im[ite,bus]^2)*lambda[ite] - v_max_slack[bus] <= Vmax[bus]^2;
###
# sub
###
problem subproblem:
p_gen, p_gen_min, p_gen_max,
q_gen, q_gen_min, q_gen_max,
lambda, p_pos, p_neg, q_pos, q_neg;
var v_re{bus in bus_i}, >= -Vmax[bus], <= Vmax[bus];
var v_im{bus in bus_i}, >= -Vmax[bus], <= Vmax[bus];
var p_or{(bus1, bus2) in branch};
var p_ex{(bus1, bus2) in branch};
var q_or{(bus1, bus2) in branch};
var q_ex{(bus1, bus2) in branch};
subject to p_or_eq{(bus1, bus2) in branch}: p_or[bus1, bus2] =
+y11_re[bus1, bus2]*v_re[bus1]*v_re[bus1]
+y11_re[bus1, bus2]*v_im[bus1]*v_im[bus1]
+y12_re[bus1, bus2]*v_re[bus1]*v_re[bus2]
+y12_re[bus1, bus2]*v_im[bus1]*v_im[bus2]
+y12_im[bus1, bus2]*v_im[bus1]*v_re[bus2]
-y12_im[bus1, bus2]*v_re[bus1]*v_im[bus2];
subject to p_ex_eq{(bus1, bus2) in branch}: p_ex[bus1, bus2] =
+y22_re[bus1,bus2]*v_re[bus2]*v_re[bus2]
+y22_re[bus1,bus2]*v_im[bus2]*v_im[bus2]
+y21_re[bus1,bus2]*v_re[bus1]*v_re[bus2]
+y21_re[bus1,bus2]*v_im[bus1]*v_im[bus2]
+y21_im[bus1,bus2]*v_re[bus1]*v_im[bus2]
-y21_im[bus1,bus2]*v_im[bus1]*v_re[bus2];
subject to q_or_eq{(bus1, bus2) in branch}: q_or[bus1, bus2] =
-y11_im[bus1,bus2]*v_re[bus1]*v_re[bus1]
-y11_im[bus1,bus2]*v_im[bus1]*v_im[bus1]
+y12_re[bus1,bus2]*v_im[bus1]*v_re[bus2]
-y12_re[bus1,bus2]*v_re[bus1]*v_im[bus2]
-y12_im[bus1,bus2]*v_re[bus1]*v_re[bus2]
-y12_im[bus1,bus2]*v_im[bus1]*v_im[bus2];
subject to q_ex_eq{(bus1, bus2) in branch}: q_ex[bus1, bus2] =
-y22_im[bus1,bus2]*v_re[bus2]*v_re[bus2]
-y22_im[bus1,bus2]*v_im[bus2]*v_im[bus2]
+y21_re[bus1,bus2]*v_re[bus1]*v_im[bus2]
-y21_re[bus1,bus2]*v_im[bus1]*v_re[bus2]
-y21_im[bus1,bus2]*v_re[bus1]*v_re[bus2]
-y21_im[bus1,bus2]*v_im[bus1]*v_im[bus2];
# convex
var ips_lambda >= 0;
subject to ref_bus_re{bus in bus_i:type[bus]==3}:v_re[bus]>= 0;
subject to ref_bus_im{bus in bus_i:type[bus]==3}:v_im[bus] = 0;
subject to ctr_v_min{bus in bus_i}: v_re[bus]*v_re[bus]+v_im[bus]*v_im[bus] >= Vmin[bus]^2;
subject to ctr_v_max{bus in bus_i}: v_re[bus]*v_re[bus]+v_im[bus]*v_im[bus] <= Vmax[bus]^2;
minimize sub_obj2:
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus1} dual_p[bus]*100*(
+y11_re[bus1, bus2]*v_re[bus1]*v_re[bus1]
+y11_re[bus1, bus2]*v_im[bus1]*v_im[bus1]
+y12_re[bus1, bus2]*v_re[bus1]*v_re[bus2]
+y12_re[bus1, bus2]*v_im[bus1]*v_im[bus2]
+y12_im[bus1, bus2]*v_im[bus1]*v_re[bus2]
-y12_im[bus1, bus2]*v_re[bus1]*v_im[bus2])
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus2} dual_p[bus]*100*(
+y22_re[bus1,bus2]*v_re[bus2]*v_re[bus2]
+y22_re[bus1,bus2]*v_im[bus2]*v_im[bus2]
+y21_re[bus1,bus2]*v_re[bus1]*v_re[bus2]
+y21_re[bus1,bus2]*v_im[bus1]*v_im[bus2]
+y21_im[bus1,bus2]*v_re[bus1]*v_im[bus2]
-y21_im[bus1,bus2]*v_im[bus1]*v_re[bus2])
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus1} dual_q[bus]*100*(
-y11_im[bus1,bus2]*v_re[bus1]*v_re[bus1]
-y11_im[bus1,bus2]*v_im[bus1]*v_im[bus1]
+y12_re[bus1,bus2]*v_im[bus1]*v_re[bus2]
-y12_re[bus1,bus2]*v_re[bus1]*v_im[bus2]
-y12_im[bus1,bus2]*v_re[bus1]*v_re[bus2]
-y12_im[bus1,bus2]*v_im[bus1]*v_im[bus2])
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus2} dual_q[bus]*100*(
-y22_im[bus1,bus2]*v_re[bus2]*v_re[bus2]
-y22_im[bus1,bus2]*v_im[bus2]*v_im[bus2]
+y21_re[bus1,bus2]*v_re[bus1]*v_im[bus2]
-y21_re[bus1,bus2]*v_im[bus1]*v_re[bus2]
-y21_im[bus1,bus2]*v_re[bus1]*v_re[bus2]
-y21_im[bus1,bus2]*v_im[bus1]*v_im[bus2])
-dual_convexity;
minimize sub_obj:
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus1} dual_p[bus]*100*p_or[bus1, bus2]
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus2} dual_p[bus]*100*p_ex[bus1, bus2]
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus1} dual_q[bus]*100*q_or[bus1, bus2]
-sum{bus in bus_i, (bus1, bus2) in branch:bus==bus2} dual_q[bus]*100*q_ex[bus1, bus2]
-dual_convexity;
subject to ips_convexity{tmp in 1..1:FEASIBILITY==0 and card(master_basis)>0}:
ips_lambda+sum{ite in master_basis} lambda[ite]=1;
subject to ips_p_balance{bus in bus_i: FEASIBILITY==0 and card(master_basis)>0}:
-(if bus in gen_i then p_gen[bus])/100
+sum{ite in master_basis, (bus1, bus2) in branch:bus==bus1} P_or[ite,bus1,bus2]*lambda[ite]
+sum{ite in master_basis, (bus1, bus2) in branch:bus==bus2} P_ex[ite,bus1,bus2]*lambda[ite]
+sum{(bus1, bus2) in branch:bus==bus1}100*p_or[bus1, bus2]
+sum{(bus1, bus2) in branch:bus==bus2}100*p_ex[bus1, bus2]
=
-Pd[bus]
;
subject to ips_q_balance{bus in bus_i: FEASIBILITY==0 and card(master_basis)>0}:
-(if bus in gen_i then q_gen[bus])/100
+sum{ite in master_basis, (bus1, bus2) in branch:bus==bus1} Q_or[ite,bus1,bus2]*lambda[ite]
+sum{ite in master_basis, (bus1, bus2) in branch:bus==bus2} Q_ex[ite,bus1,bus2]*lambda[ite]
+sum{(bus1, bus2) in branch:bus==bus1}100*q_or[bus1, bus2]
+sum{(bus1, bus2) in branch:bus==bus2}100*q_ex[bus1, bus2]
=
-Qd[bus]
;