OPF in AMPL
Optimal Power Flow
1########## SETS ##########
2
3set N;
4set G;
5set L;
6
7########## PARAMETERS ##########
8
9param OPF_TYPE symbolic;
10param CONNECTIVITY symbolic;
11param BASEMVA;
12param MAXANGLE;
13param MINANGLE;
14param COST_2 {G};
15param COST_1 {G};
16param COST_0 {G};
17param PG {G};
18param QG {G};
19param PMAX {G};
20param PMIN {G};
21param QMAX {G};
22param QMIN {G};
23param GEN_BUS {G};
24param PD {N};
25param QD {N};
26param GS {N};
27param BS {N};
28param VMAX {N};
29param VMIN {N};
30param AMAX {N};
31param AMIN {N};
32param VOL0 {N};
33param ANG0 {N};
34param VOLR0 {N};
35param VOLI0 {N};
36param CF {L,N};
37param CT {L,N};
38param CG {G,N};
39param GFF {L};
40param BFF {L};
41param GFT {L};
42param BFT {L};
43param GTF {L};
44param BTF {L};
45param GTT {L};
46param BTT {L};
47param BUS_I {N};
48param BUS_TYPE {N};
49param BUS_AREA {N};
50param VM {N};
51param VA {N};
52param BASE_KV {N};
53param ZONE {N};
54param F_BUS {L};
55param T_BUS {L};
56param BR_R {L};
57param BR_X {L};
58param BR_B {L};
59param RATE_A {L};
60param RATE_B {L};
61param RATE_C {L};
62param PFMAX {L};
63param PFMIN {L};
64param QFMAX {L};
65param QFMIN {L};
66param TAP {L};
67param SHIFT {L};
68param BR_STATUS {L};
69param BR_SWITCH {L};
70param ANGMIN {L};
71param ANGMAX {L};
72param VG {G};
73param MBASE {G};
74param GEN_STATUS {G};
75param PC1 {G};
76param PC2 {G};
77param QC1MIN {G};
78param QC1MAX {G};
79param QC2MIN {G};
80param QC2MAX {G};
81param RAMP_AGC {G};
82param RAMP_10 {G};
83param RAMP_30 {G};
84param RAMP_Q {G};
85param APF {G};
86param MODEL {G};
87param STARTUP {G};
88param SHUTDOWN {G};
89param NCOST {G};
90param PG0 {G};
91param QG0 {G};
92param PF0 {L};
93param QF0 {L};
94param PT0 {L};
95param QT0 {L};
96param PFLODC {L};
97param PFUPDC {L};
98param PFLOAC {L};
99param PFUPAC {L};
100param PTLOAC {L};
101param PTUPAC {L};
102param QFLOAC {L};
103param QFUPAC {L};
104param QTLOAC {L};
105param QTUPAC {L};
106param COSFTMAX {L};
107param COSFTMIN {L};
108param SINFTMAX {L};
109param SINFTMIN {L};
110param VRMAX {N};
111param VRMIN {N};
112param VIMAX {N};
113param VIMIN {N};
114
115########## VARIABLES ##########
116var Pg {g in G} >= PMIN[g], <= PMAX[g]:= PG0[g];
117var Qg {g in G} >= QMIN[g], <= QMAX[g]:= QG0[g];
118var Pf {l in L} := PF0[l];
119var Pt {l in L} := PT0[l];
120var Qf {l in L} := QF0[l];
121var Qt {l in L} := QT0[l];
122var Vm {n in N} >= VMIN[n], <= VMAX[n] := VOL0[n];
123var Vr {n in N} >= VRMIN[n], <=VRMAX[n] := VOLR0[n];
124var Vi {n in N} >= VIMIN[n], <=VIMAX[n] := VOLI0[n];
125var V2 {n in N} >= VMIN[n]^2, <= VMAX[n]^2 := VOL0[n]^2;
126var Va {n in N} >= AMIN[n], <= AMAX[n] := ANG0[n];
127var cosft {l in L} >= COSFTMIN[l], <= COSFTMAX[l] := VOL0[F_BUS[l]]*VOL0[T_BUS[l]]*cos(ANG0[F_BUS[l]]-ANG0[T_BUS[l]]);
128var sinft {l in L} >= SINFTMIN[l], <= SINFTMAX[l] := VOL0[F_BUS[l]]*VOL0[T_BUS[l]]*sin(ANG0[F_BUS[l]]-ANG0[T_BUS[l]]);
129var Pfa {l in L} := PF0[l];
130var Pta {l in L} := PT0[l];
131var Qfa {l in L} := QF0[l];
132var Qta {l in L} := QT0[l];
133var u {n in N} >= 1, <= card(N);
134var status {l in L} binary;
135var statusf {l in L} binary;
136var statust {l in L} binary;
137
138########## POWER BALANCE ##########
139
140subject to active_power_balance_1 {n in N:OPF_TYPE='dc'}:
141 sum {g in G} CG[g,n] * Pg[g] - PD[n] = sum {l in L} (CF[l,n] * Pf[l] + CT[l,n] * Pt[l]);
142
143subject to active_power_balance_2 {n in N:OPF_TYPE='acpolar'}:
144 sum {g in G} CG[g,n] * Pg[g] - PD[n] = GS[n]*Vm[n]*Vm[n] + sum {l in L} (CF[l,n] * Pf[l] + CT[l,n] * Pt[l]);
145
146subject to reactive_power_balance_2 {n in N:OPF_TYPE='acpolar'}:
147 sum {g in G} CG[g,n] * Qg[g] - QD[n] = -BS[n]*Vm[n]*Vm[n] + sum {l in L} (CF[l,n] * Qf[l] + CT[l,n] * Qt[l]);
148
149subject to active_power_balance_3_4 {n in N: OPF_TYPE = 'acrect' or OPF_TYPE = 'acjabr'}:
150 sum {g in G} CG[g,n] * Pg[g] - PD[n] = GS[n]*V2[n] + sum {l in L} (CF[l,n] * Pf[l] + CT[l,n] * Pt[l]);
151
152subject to reactive_power_balance_3_4 {n in N: OPF_TYPE = 'acrect' or OPF_TYPE = 'acjabr'}:
153 sum {g in G} CG[g,n] * Qg[g] - QD[n] = -BS[n]*V2[n] + sum {l in L} (CF[l,n] * Qf[l] + CT[l,n] * Qt[l]);
154
155########## POWER FLOW DEFINITIONS (BR_SWITCH = 0) ##########
156
157subject to active_flow_from_0 {l in L:BR_SWITCH[l] == 0}:
158 Pf[l] = 0;
159
160subject to active_flow_to_0 {l in L:BR_SWITCH[l] == 0}:
161 Pt[l] = 0;
162
163subject to reactive_flow_from_0 {l in L:BR_SWITCH[l] == 0}:
164 Qf[l] = 0;
165
166subject to reactive_flow_to_0 {l in L:BR_SWITCH[l] == 0}:
167 Qt[l] = 0;
168
169########## POWER FLOW DEFINITIONS (BR_SWITCH = 1) ##########
170
171subject to active_flow_from_1 {l in L:OPF_TYPE='dc' and BR_SWITCH[l] == 1}:
172 Pf[l] = (1 / BR_X[l]) * (Va[F_BUS[l]] - Va[T_BUS[l]]);
173
174subject to active_flow_to_1 {l in L:OPF_TYPE='dc' and BR_SWITCH[l] == 1}:
175 Pt[l] = (1 / BR_X[l]) * (Va[T_BUS[l]] - Va[F_BUS[l]]);
176
177subject to active_flow_from_2 {l in L:OPF_TYPE='acpolar' and BR_SWITCH[l] == 1}:
178 Pf[l] = GFF[l]*Vm[F_BUS[l]]*Vm[F_BUS[l]] + Vm[F_BUS[l]]*Vm[T_BUS[l]]*(GFT[l]*cos(Va[F_BUS[l]]-Va[T_BUS[l]])+BFT[l]*sin(Va[F_BUS[l]]-Va[T_BUS[l]]));
179
180subject to active_flow_to_2 {l in L:OPF_TYPE='acpolar' and BR_SWITCH[l] == 1}:
181 Pt[l] = GTT[l]*Vm[T_BUS[l]]*Vm[T_BUS[l]] + Vm[F_BUS[l]]*Vm[T_BUS[l]]*(GTF[l]*cos(Va[T_BUS[l]]-Va[F_BUS[l]])+BTF[l]*sin(Va[T_BUS[l]]-Va[F_BUS[l]]));
182
183subject to reactive_flow_from_2 {l in L:OPF_TYPE='acpolar' and BR_SWITCH[l] == 1}:
184 Qf[l] = -BFF[l]*Vm[F_BUS[l]]*Vm[F_BUS[l]] - Vm[F_BUS[l]]*Vm[T_BUS[l]]*(BFT[l]*cos(Va[F_BUS[l]]-Va[T_BUS[l]])-GFT[l]*sin(Va[F_BUS[l]]-Va[T_BUS[l]]));
185
186subject to reactive_flow_to_2 {l in L:OPF_TYPE='acpolar' and BR_SWITCH[l] == 1}:
187 Qt[l] = -BTT[l]*Vm[T_BUS[l]]*Vm[T_BUS[l]] - Vm[F_BUS[l]]*Vm[T_BUS[l]]*(BTF[l]*cos(Va[T_BUS[l]]-Va[F_BUS[l]])-GTF[l]*sin(Va[T_BUS[l]]-Va[F_BUS[l]]));
188
189subject to active_flow_from_3_4 {l in L:(OPF_TYPE='acrect' or OPF_TYPE = 'acjabr') and BR_SWITCH[l] == 1}:
190 Pf[l] = GFF[l]*V2[F_BUS[l]] + GFT[l]*cosft[l] + BFT[l]*sinft[l];
191
192subject to active_flow_to_3_4 {l in L:(OPF_TYPE='acrect' or OPF_TYPE = 'acjabr') and BR_SWITCH[l] == 1}:
193 Pt[l] = GTT[l]*V2[T_BUS[l]] + GTF[l]*cosft[l] - BTF[l]*sinft[l];
194
195subject to reactive_flow_from_3_4 {l in L:(OPF_TYPE='acrect' or OPF_TYPE = 'acjabr') and BR_SWITCH[l] == 1}:
196 Qf[l] = -BFF[l]*V2[F_BUS[l]] - BFT[l]*cosft[l] + GFT[l]*sinft[l];
197
198subject to reactive_flow_to_3_4 {l in L:(OPF_TYPE='acrect' or OPF_TYPE = 'acjabr') and BR_SWITCH[l] == 1}:
199 Qt[l] = -BTT[l]*V2[T_BUS[l]] - BTF[l]*cosft[l] - GTF[l]*sinft[l];
200
201########## POWER FLOW DEFINITIONS (BR_SWITCH = 2) ##########
202
203subject to active_flow_from_1_switch {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 2}:
204 Pf[l] = status[l] * (1 / BR_X[l]) * (Va[F_BUS[l]] - Va[T_BUS[l]]);
205
206subject to active_flow_to_1_switch {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 2}:
207 Pt[l] = status[l] * (1 / BR_X[l]) * (Va[T_BUS[l]] - Va[F_BUS[l]]);
208
209subject to active_flow_from_2_switch {l in L: OPF_TYPE == 'acpolar' and BR_SWITCH[l] == 2}:
210 Pf[l] = status[l] * (GFF[l]*Vm[F_BUS[l]]*Vm[F_BUS[l]] + Vm[F_BUS[l]]*Vm[T_BUS[l]]*(GFT[l]*cos(Va[F_BUS[l]]-Va[T_BUS[l]])+BFT[l]*sin(Va[F_BUS[l]]-Va[T_BUS[l]])));
211
212subject to active_flow_to_2_switch {l in L: OPF_TYPE == 'acpolar' and BR_SWITCH[l] == 2}:
213 Pt[l] = status[l] * (GTT[l]*Vm[T_BUS[l]]*Vm[T_BUS[l]] + Vm[F_BUS[l]]*Vm[T_BUS[l]]*(GTF[l]*cos(Va[T_BUS[l]]-Va[F_BUS[l]])+BTF[l]*sin(Va[T_BUS[l]]-Va[F_BUS[l]])));
214
215subject to reactive_flow_from_2_switch {l in L: OPF_TYPE == 'acpolar' and BR_SWITCH[l] == 2}:
216 Qf[l] = status[l] * (-BFF[l]*Vm[F_BUS[l]]*Vm[F_BUS[l]] - Vm[F_BUS[l]]*Vm[T_BUS[l]]*(BFT[l]*cos(Va[F_BUS[l]]-Va[T_BUS[l]])-GFT[l]*sin(Va[F_BUS[l]]-Va[T_BUS[l]])));
217
218subject to reactive_flow_to_2_switch {l in L: OPF_TYPE == 'acpolar' and BR_SWITCH[l] == 2}:
219 Qt[l] = status[l] * (-BTT[l]*Vm[T_BUS[l]]*Vm[T_BUS[l]] - Vm[F_BUS[l]]*Vm[T_BUS[l]]*(BTF[l]*cos(Va[T_BUS[l]]-Va[F_BUS[l]])-GTF[l]*sin(Va[T_BUS[l]]-Va[F_BUS[l]])));
220
221subject to active_flow_from_3_4_switch {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 2}:
222 Pf[l] = status[l] * (GFF[l]*V2[F_BUS[l]] + GFT[l]*cosft[l] + BFT[l]*sinft[l]);
223
224subject to active_flow_to_3_4_switch {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 2}:
225 Pt[l] = status[l] * (GTT[l]*V2[T_BUS[l]] + GTF[l]*cosft[l] - BTF[l]*sinft[l]);
226
227subject to reactive_flow_from_3_4_switch {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 2}:
228 Qf[l] = status[l] * (-BFF[l]*V2[F_BUS[l]] - BFT[l]*cosft[l] + GFT[l]*sinft[l]);
229
230subject to reactive_flow_to_3_4_switch {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 2}:
231 Qt[l] = status[l] * (-BTT[l]*V2[T_BUS[l]] - BTF[l]*cosft[l] - GTF[l]*sinft[l]);
232
233########## POWER FLOW DEFINITIONS (BR_SWITCH = 3) ##########
234
235subject to active_flow_from_1_bigm {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
236 Pfa[l] = (1 / BR_X[l]) * (Va[F_BUS[l]] - Va[T_BUS[l]]);
237
238subject to active_flow_to_1_bigm {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
239 Pta[l] = (1 / BR_X[l]) * (Va[T_BUS[l]] - Va[F_BUS[l]]);
240
241subject to Pfa_lower_1 {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
242 PFLODC[l] * (1 - status[l]) <= -Pf[l] + Pfa[l];
243
244subject to Pfa_upper_1 {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
245 -Pf[l] + Pfa[l] <= PFUPDC[l] * (1 - status[l]);
246
247subject to Pta_lower_1 {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
248 -PFUPDC[l] * (1 - status[l]) <= -Pt[l] + Pta[l];
249
250subject to Pta_upper_1 {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
251 -Pt[l] + Pta[l] <= -PFLODC[l] * (1 - status[l]);
252
253subject to active_flow_from_2_bigm {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
254 Pfa[l] = GFF[l] * V2[F_BUS[l]] + GFT[l] * cosft[l] + BFT[l] * sinft[l];
255
256subject to active_flow_to_2_bigm {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
257 Pta[l] = GTT[l] * V2[T_BUS[l]] + GTF[l] * cosft[l] - BTF[l] * sinft[l];
258
259subject to reactive_flow_from_2_bigm {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
260 Qfa[l] = -BFF[l] * V2[F_BUS[l]] - BFT[l] * cosft[l] + GFT[l] * sinft[l];
261
262subject to reactive_flow_to_2_bigm {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
263 Qta[l] = -BTT[l] * V2[T_BUS[l]] - BTF[l] * cosft[l] - GTF[l] * sinft[l];
264
265subject to Pfa_lower_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
266 PFLOAC[l] * (1 - status[l]) <= -Pf[l] + Pfa[l];
267
268subject to Pfa_upper_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
269 -Pf[l] + Pfa[l] <= PFUPAC[l] * (1 - status[l]);
270
271subject to Pta_lower_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
272 PTLOAC[l] * (1 - status[l]) <= -Pt[l] + Pta[l];
273
274subject to Pta_upper_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
275 -Pt[l] + Pta[l] <= PTUPAC[l] * (1 - status[l]);
276
277subject to Qfa_lower_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
278 QFLOAC[l] * (1 - status[l]) <= -Qf[l] + Qfa[l];
279
280subject to Qfa_upper_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
281 -Qf[l] + Qfa[l] <= QFUPAC[l] * (1 - status[l]);
282
283subject to Qta_lower_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
284 QTLOAC[l] * (1 - status[l]) <= -Qt[l] + Qta[l];
285
286subject to Qta_upper_2 {l in L: (OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
287 -Qt[l] + Qta[l] <= QTUPAC[l] * (1 - status[l]);
288
289########## POWER FLOW LIMITS ##########
290
291subject to flowf_limits_1 {l in L:OPF_TYPE='dc' and (BR_SWITCH[l] == 1 or BR_SWITCH[l] == 2)}:
292 PFMIN[l] <= Pf[l] <= PFMAX[l];
293
294subject to flowt_limits_1 {l in L:OPF_TYPE='dc' and (BR_SWITCH[l] == 1 or BR_SWITCH[l] == 2)}:
295 -PFMAX[l] <= Pt[l] <= -PFMIN[l];
296
297subject to flowf_limits_dc_lower {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
298 Pf[l] >= PFMIN[l] * status[l];
299
300subject to flowf_limits_dc_upper {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
301 Pf[l] <= PFMAX[l] * status[l];
302
303subject to flowt_limits_dc_lower {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
304 Pt[l] >= -PFMAX[l] * status[l];
305
306subject to flowt_limits_dc_upper {l in L: OPF_TYPE == 'dc' and BR_SWITCH[l] == 3}:
307 Pt[l] <= -PFMIN[l] * status[l];
308
309subject to flow_limits_from_23_4 {l in L:(OPF_TYPE='acpolar' or OPF_TYPE = 'acrect' or OPF_TYPE = 'acjabr') and (BR_SWITCH[l] == 1 or BR_SWITCH[l] == 2)}:
310 Pf[l]^2 + Qf[l]^2 <= RATE_A[l]^2;
311
312subject to flow_limits_to_23_4 {l in L:(OPF_TYPE='acpolar' or OPF_TYPE = 'acrect' or OPF_TYPE = 'acjabr') and (BR_SWITCH[l] == 1 or BR_SWITCH[l] == 2)}:
313 Pt[l]^2 + Qt[l]^2 <= RATE_A[l]^2;
314
315subject to flow_limits_from_acrect {l in L: (OPF_TYPE == 'acpolar' or OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
316 Pf[l]^2 + Qf[l]^2 <= RATE_A[l]^2 * status[l];
317
318subject to flow_limits_to_acrect {l in L: (OPF_TYPE == 'acpolar' or OPF_TYPE == 'acrect' or OPF_TYPE == 'acjabr') and BR_SWITCH[l] == 3}:
319 Pt[l]^2 + Qt[l]^2 <= RATE_A[l]^2 * status[l];
320
321########## RECTANGULAR DEFINITIONS ##########
322
323subject to eq_vol_squared {n in N:OPF_TYPE = 'acrect'}:
324 V2[n] == Vr[n]*Vr[n] + Vi[n]*Vi[n];
325
326subject to eq_cosft {l in L:OPF_TYPE = 'acrect'}:
327 cosft[l] == Vr[F_BUS[l]]*Vr[T_BUS[l]] + Vi[F_BUS[l]]*Vi[T_BUS[l]];
328
329subject to eq_sinft {l in L:OPF_TYPE = 'acrect'}:
330 sinft[l] == Vi[F_BUS[l]]*Vr[T_BUS[l]] - Vr[F_BUS[l]]*Vi[T_BUS[l]];
331
332########## JABR RELAXATION ##########
333
334subject to jabr_relaxation_ft {l in L:OPF_TYPE = 'acrect' or OPF_TYPE = 'acjabr'}:
335 cosft[l]^2 + sinft[l]^2 <= V2[F_BUS[l]]*V2[T_BUS[l]];
336
337# TODO: Add the relaxation by Muñoz where losses are positive?
338
339########## SLACK BUS ##########
340
341subject to eq_slack {n in N: BUS_TYPE[n] == 3}:
342 Va[n] == 0;
343
344subject to eq_slack_imag {n in N: BUS_TYPE[n] == 3}:
345 Vi[n] == 0;
346
347########## CONNECTIVITY CONSTRAINTS ##########
348
349subject to status_split {l in L: CONNECTIVITY = 'on' and (BR_SWITCH[l] == 2 or BR_SWITCH[l] == 3)}:
350 statusf[l] + statust[l] == status[l];
351
352subject to connectivity {n in N: CONNECTIVITY = 'on'}:
353 sum {l in L: F_BUS[l] == n} statusf[l] + sum {l in L: T_BUS[l] == n} statust[l] >= 1;
354
355subject to mtz_connectivity_f {l in L: CONNECTIVITY = 'on' and F_BUS[l] != 0 and T_BUS[l] != 0}:
356 u[F_BUS[l]] - u[T_BUS[l]] + card(N) * statusf[l] <= card(N) - 1;
357
358subject to mtz_connectivity_t {l in L: CONNECTIVITY = 'on' and F_BUS[l] != 0 and T_BUS[l] != 0}:
359 u[T_BUS[l]] - u[F_BUS[l]] + card(N) * statust[l] <= card(N) - 1;
360
361subject to eq_slack_mtz:
362 u[0] == 1;
363
364########## STATUS FIX ##########
365
366subject to fix_status_0 {l in L: BR_SWITCH[l] == 0}:
367 status[l] == 0;
368
369subject to fix_status_1 {l in L: BR_SWITCH[l] == 1}:
370 status[l] == 1;
371
372# TODO: Split AMPL models?