1 static char help[] = "Using finite difference for the problem in ex9busopt.c \n\n";
2
3 /*
4 Use finite difference approximations to solve the same optimization problem as in ex9busopt.c.
5 */
6
7 #include <petsctao.h>
8 #include <petscts.h>
9 #include <petscdm.h>
10 #include <petscdmda.h>
11 #include <petscdmcomposite.h>
12
13 PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
14
15 #define freq 60
16 #define w_s (2 * PETSC_PI * freq)
17
18 /* Sizes and indices */
19 const PetscInt nbus = 9; /* Number of network buses */
20 const PetscInt ngen = 3; /* Number of generators */
21 const PetscInt nload = 3; /* Number of loads */
22 const PetscInt gbus[3] = {0, 1, 2}; /* Buses at which generators are incident */
23 const PetscInt lbus[3] = {4, 5, 7}; /* Buses at which loads are incident */
24
25 /* Generator real and reactive powers (found via loadflow) */
26 PetscScalar PG[3] = {0.69, 1.59, 0.69};
27 /* PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};*/
28 const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
29 /* Generator constants */
30 const PetscScalar H[3] = {23.64, 6.4, 3.01}; /* Inertia constant */
31 const PetscScalar Rs[3] = {0.0, 0.0, 0.0}; /* Stator Resistance */
32 const PetscScalar Xd[3] = {0.146, 0.8958, 1.3125}; /* d-axis reactance */
33 const PetscScalar Xdp[3] = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
34 const PetscScalar Xq[3] = {0.4360, 0.8645, 1.2578}; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
35 const PetscScalar Xqp[3] = {0.0969, 0.1969, 0.25}; /* q-axis transient reactance */
36 const PetscScalar Td0p[3] = {8.96, 6.0, 5.89}; /* d-axis open circuit time constant */
37 const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6}; /* q-axis open circuit time constant */
38 PetscScalar M[3]; /* M = 2*H/w_s */
39 PetscScalar D[3]; /* D = 0.1*M */
40
41 PetscScalar TM[3]; /* Mechanical Torque */
42 /* Exciter system constants */
43 const PetscScalar KA[3] = {20.0, 20.0, 20.0}; /* Voltage regulartor gain constant */
44 const PetscScalar TA[3] = {0.2, 0.2, 0.2}; /* Voltage regulator time constant */
45 const PetscScalar KE[3] = {1.0, 1.0, 1.0}; /* Exciter gain constant */
46 const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
47 const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
48 const PetscScalar TF[3] = {0.35, 0.35, 0.35}; /* Feedback stabilizer time constant */
49 const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
50 const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */
51
52 PetscScalar Vref[3];
53 /* Load constants
54 We use a composite load model that describes the load and reactive powers at each time instant as follows
55 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
56 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
57 where
58 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
59 ld_alphap,ld_alphap - Percentage contribution (weights) or loads
60 P_D0 - Real power load
61 Q_D0 - Reactive power load
62 V_m(t) - Voltage magnitude at time t
63 V_m0 - Voltage magnitude at t = 0
64 ld_betap, ld_betaq - exponents describing the load model for real and reactive part
65
66 Note: All loads have the same characteristic currently.
67 */
68 const PetscScalar PD0[3] = {1.25, 0.9, 1.0};
69 const PetscScalar QD0[3] = {0.5, 0.3, 0.35};
70 const PetscInt ld_nsegsp[3] = {3, 3, 3};
71 const PetscScalar ld_alphap[3] = {1.0, 0.0, 0.0};
72 const PetscScalar ld_betap[3] = {2.0, 1.0, 0.0};
73 const PetscInt ld_nsegsq[3] = {3, 3, 3};
74 const PetscScalar ld_alphaq[3] = {1.0, 0.0, 0.0};
75 const PetscScalar ld_betaq[3] = {2.0, 1.0, 0.0};
76
77 typedef struct {
78 DM dmgen, dmnet; /* DMs to manage generator and network subsystem */
79 DM dmpgrid; /* Composite DM to manage the entire power grid */
80 Mat Ybus; /* Network admittance matrix */
81 Vec V0; /* Initial voltage vector (Power flow solution) */
82 PetscReal tfaulton, tfaultoff; /* Fault on and off times */
83 PetscInt faultbus; /* Fault bus */
84 PetscScalar Rfault;
85 PetscReal t0, tmax;
86 PetscInt neqs_gen, neqs_net, neqs_pgrid;
87 Mat Sol; /* Matrix to save solution at each time step */
88 PetscInt stepnum;
89 PetscBool alg_flg;
90 PetscReal t;
91 IS is_diff; /* indices for differential equations */
92 IS is_alg; /* indices for algebraic equations */
93 PetscReal freq_u, freq_l; /* upper and lower frequency limit */
94 PetscInt pow; /* power coefficient used in the cost function */
95 Vec vec_q;
96 } Userctx;
97
98 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar * Fr,PetscScalar * Fi)99 PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
100 {
101 PetscFunctionBegin;
102 *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
103 *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106
107 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar * Fd,PetscScalar * Fq)108 PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
109 {
110 PetscFunctionBegin;
111 *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
112 *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
SetInitialGuess(Vec X,Userctx * user)116 PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
117 {
118 Vec Xgen, Xnet;
119 PetscScalar *xgen, *xnet;
120 PetscInt i, idx = 0;
121 PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
122 PetscScalar Eqp, Edp, delta;
123 PetscScalar Efd, RF, VR; /* Exciter variables */
124 PetscScalar Id, Iq; /* Generator dq axis currents */
125 PetscScalar theta, Vd, Vq, SE;
126
127 PetscFunctionBegin;
128 M[0] = 2 * H[0] / w_s;
129 M[1] = 2 * H[1] / w_s;
130 M[2] = 2 * H[2] / w_s;
131 D[0] = 0.1 * M[0];
132 D[1] = 0.1 * M[1];
133 D[2] = 0.1 * M[2];
134
135 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
136
137 /* Network subsystem initialization */
138 PetscCall(VecCopy(user->V0, Xnet));
139
140 /* Generator subsystem initialization */
141 PetscCall(VecGetArray(Xgen, &xgen));
142 PetscCall(VecGetArray(Xnet, &xnet));
143
144 for (i = 0; i < ngen; i++) {
145 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
146 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
147 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
148 Vm2 = Vm * Vm;
149 IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
150 IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;
151
152 delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
153
154 theta = PETSC_PI / 2.0 - delta;
155
156 Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
157 Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
158
159 Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
160 Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
161
162 Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
163 Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
164
165 TM[i] = PG[i];
166
167 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
168 xgen[idx] = Eqp;
169 xgen[idx + 1] = Edp;
170 xgen[idx + 2] = delta;
171 xgen[idx + 3] = w_s;
172
173 idx = idx + 4;
174
175 xgen[idx] = Id;
176 xgen[idx + 1] = Iq;
177
178 idx = idx + 2;
179
180 /* Exciter */
181 Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
182 SE = k1[i] * PetscExpScalar(k2[i] * Efd);
183 VR = KE[i] * Efd + SE;
184 RF = KF[i] * Efd / TF[i];
185
186 xgen[idx] = Efd;
187 xgen[idx + 1] = RF;
188 xgen[idx + 2] = VR;
189
190 Vref[i] = Vm + (VR / KA[i]);
191
192 idx = idx + 3;
193 }
194
195 PetscCall(VecRestoreArray(Xgen, &xgen));
196 PetscCall(VecRestoreArray(Xnet, &xnet));
197
198 /* PetscCall(VecView(Xgen,0)); */
199 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
200 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
201 PetscFunctionReturn(PETSC_SUCCESS);
202 }
203
204 /* Computes F = [-f(x,y);g(x,y)] */
ResidualFunction(SNES snes,Vec X,Vec F,Userctx * user)205 PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
206 {
207 Vec Xgen, Xnet, Fgen, Fnet;
208 PetscScalar *xgen, *xnet, *fgen, *fnet;
209 PetscInt i, idx = 0;
210 PetscScalar Vr, Vi, Vm, Vm2;
211 PetscScalar Eqp, Edp, delta, w; /* Generator variables */
212 PetscScalar Efd, RF, VR; /* Exciter variables */
213 PetscScalar Id, Iq; /* Generator dq axis currents */
214 PetscScalar Vd, Vq, SE;
215 PetscScalar IGr, IGi, IDr, IDi;
216 PetscScalar Zdq_inv[4], det;
217 PetscScalar PD, QD, Vm0, *v0;
218 PetscInt k;
219
220 PetscFunctionBegin;
221 PetscCall(VecZeroEntries(F));
222 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
223 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
224 PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
225 PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));
226
227 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
228 The generator current injection, IG, and load current injection, ID are added later
229 */
230 /* Note that the values in Ybus are stored assuming the imaginary current balance
231 equation is ordered first followed by real current balance equation for each bus.
232 Thus imaginary current contribution goes in location 2*i, and
233 real current contribution in 2*i+1
234 */
235 PetscCall(MatMult(user->Ybus, Xnet, Fnet));
236
237 PetscCall(VecGetArray(Xgen, &xgen));
238 PetscCall(VecGetArray(Xnet, &xnet));
239 PetscCall(VecGetArray(Fgen, &fgen));
240 PetscCall(VecGetArray(Fnet, &fnet));
241
242 /* Generator subsystem */
243 for (i = 0; i < ngen; i++) {
244 Eqp = xgen[idx];
245 Edp = xgen[idx + 1];
246 delta = xgen[idx + 2];
247 w = xgen[idx + 3];
248 Id = xgen[idx + 4];
249 Iq = xgen[idx + 5];
250 Efd = xgen[idx + 6];
251 RF = xgen[idx + 7];
252 VR = xgen[idx + 8];
253
254 /* Generator differential equations */
255 fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
256 fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
257 fgen[idx + 2] = -w + w_s;
258 fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
259
260 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
261 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
262
263 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
264 /* Algebraic equations for stator currents */
265 det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
266
267 Zdq_inv[0] = Rs[i] / det;
268 Zdq_inv[1] = Xqp[i] / det;
269 Zdq_inv[2] = -Xdp[i] / det;
270 Zdq_inv[3] = Rs[i] / det;
271
272 fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
273 fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
274
275 /* Add generator current injection to network */
276 PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
277
278 fnet[2 * gbus[i]] -= IGi;
279 fnet[2 * gbus[i] + 1] -= IGr;
280
281 Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
282
283 SE = k1[i] * PetscExpScalar(k2[i] * Efd);
284
285 /* Exciter differential equations */
286 fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
287 fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
288 fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
289
290 idx = idx + 9;
291 }
292
293 PetscCall(VecGetArray(user->V0, &v0));
294 for (i = 0; i < nload; i++) {
295 Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
296 Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
297 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
298 Vm2 = Vm * Vm;
299 Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
300 PD = QD = 0.0;
301 for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
302 for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
303
304 /* Load currents */
305 IDr = (PD * Vr + QD * Vi) / Vm2;
306 IDi = (-QD * Vr + PD * Vi) / Vm2;
307
308 fnet[2 * lbus[i]] += IDi;
309 fnet[2 * lbus[i] + 1] += IDr;
310 }
311 PetscCall(VecRestoreArray(user->V0, &v0));
312
313 PetscCall(VecRestoreArray(Xgen, &xgen));
314 PetscCall(VecRestoreArray(Xnet, &xnet));
315 PetscCall(VecRestoreArray(Fgen, &fgen));
316 PetscCall(VecRestoreArray(Fnet, &fnet));
317
318 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
319 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
320 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
321 PetscFunctionReturn(PETSC_SUCCESS);
322 }
323
324 /* \dot{x} - f(x,y)
325 g(x,y) = 0
326 */
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx * user)327 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
328 {
329 SNES snes;
330 PetscScalar *f;
331 const PetscScalar *xdot;
332 PetscInt i;
333
334 PetscFunctionBegin;
335 user->t = t;
336
337 PetscCall(TSGetSNES(ts, &snes));
338 PetscCall(ResidualFunction(snes, X, F, user));
339 PetscCall(VecGetArray(F, &f));
340 PetscCall(VecGetArrayRead(Xdot, &xdot));
341 for (i = 0; i < ngen; i++) {
342 f[9 * i] += xdot[9 * i];
343 f[9 * i + 1] += xdot[9 * i + 1];
344 f[9 * i + 2] += xdot[9 * i + 2];
345 f[9 * i + 3] += xdot[9 * i + 3];
346 f[9 * i + 6] += xdot[9 * i + 6];
347 f[9 * i + 7] += xdot[9 * i + 7];
348 f[9 * i + 8] += xdot[9 * i + 8];
349 }
350 PetscCall(VecRestoreArray(F, &f));
351 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
352 PetscFunctionReturn(PETSC_SUCCESS);
353 }
354
355 /* This function is used for solving the algebraic system only during fault on and
356 off times. It computes the entire F and then zeros out the part corresponding to
357 differential equations
358 F = [0;g(y)];
359 */
AlgFunction(SNES snes,Vec X,Vec F,PetscCtx ctx)360 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, PetscCtx ctx)
361 {
362 Userctx *user = (Userctx *)ctx;
363 PetscScalar *f;
364 PetscInt i;
365
366 PetscFunctionBegin;
367 PetscCall(ResidualFunction(snes, X, F, user));
368 PetscCall(VecGetArray(F, &f));
369 for (i = 0; i < ngen; i++) {
370 f[9 * i] = 0;
371 f[9 * i + 1] = 0;
372 f[9 * i + 2] = 0;
373 f[9 * i + 3] = 0;
374 f[9 * i + 6] = 0;
375 f[9 * i + 7] = 0;
376 f[9 * i + 8] = 0;
377 }
378 PetscCall(VecRestoreArray(F, &f));
379 PetscFunctionReturn(PETSC_SUCCESS);
380 }
381
PreallocateJacobian(Mat J,Userctx * user)382 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
383 {
384 PetscInt *d_nnz;
385 PetscInt i, idx = 0, start = 0;
386 PetscInt ncols;
387
388 PetscFunctionBegin;
389 PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
390 for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
391 /* Generator subsystem */
392 for (i = 0; i < ngen; i++) {
393 d_nnz[idx] += 3;
394 d_nnz[idx + 1] += 2;
395 d_nnz[idx + 2] += 2;
396 d_nnz[idx + 3] += 5;
397 d_nnz[idx + 4] += 6;
398 d_nnz[idx + 5] += 6;
399
400 d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
401 d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
402
403 d_nnz[idx + 6] += 2;
404 d_nnz[idx + 7] += 2;
405 d_nnz[idx + 8] += 5;
406
407 idx = idx + 9;
408 }
409
410 start = user->neqs_gen;
411
412 for (i = 0; i < nbus; i++) {
413 PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
414 d_nnz[start + 2 * i] += ncols;
415 d_nnz[start + 2 * i + 1] += ncols;
416 PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
417 }
418
419 PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
420
421 PetscCall(PetscFree(d_nnz));
422 PetscFunctionReturn(PETSC_SUCCESS);
423 }
424
425 /*
426 J = [-df_dx, -df_dy
427 dg_dx, dg_dy]
428 */
ResidualJacobian(SNES snes,Vec X,Mat J,Mat B,PetscCtx ctx)429 PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, PetscCtx ctx)
430 {
431 Userctx *user = (Userctx *)ctx;
432 Vec Xgen, Xnet;
433 PetscScalar *xgen, *xnet;
434 PetscInt i, idx = 0;
435 PetscScalar Vr, Vi, Vm, Vm2;
436 PetscScalar Eqp, Edp, delta; /* Generator variables */
437 PetscScalar Efd; /* Exciter variables */
438 PetscScalar Id, Iq; /* Generator dq axis currents */
439 PetscScalar Vd, Vq;
440 PetscScalar val[10];
441 PetscInt row[2], col[10];
442 PetscInt net_start = user->neqs_gen;
443 PetscScalar Zdq_inv[4], det;
444 PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
445 PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
446 PetscScalar dSE_dEfd;
447 PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
448 PetscInt ncols;
449 const PetscInt *cols;
450 const PetscScalar *yvals;
451 PetscInt k;
452 PetscScalar PD, QD, Vm0, *v0, Vm4;
453 PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
454 PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
455
456 PetscFunctionBegin;
457 PetscCall(MatZeroEntries(B));
458 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
459 PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
460
461 PetscCall(VecGetArray(Xgen, &xgen));
462 PetscCall(VecGetArray(Xnet, &xnet));
463
464 /* Generator subsystem */
465 for (i = 0; i < ngen; i++) {
466 Eqp = xgen[idx];
467 Edp = xgen[idx + 1];
468 delta = xgen[idx + 2];
469 Id = xgen[idx + 4];
470 Iq = xgen[idx + 5];
471 Efd = xgen[idx + 6];
472
473 /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
474 row[0] = idx;
475 col[0] = idx;
476 col[1] = idx + 4;
477 col[2] = idx + 6;
478 val[0] = 1 / Td0p[i];
479 val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
480 val[2] = -1 / Td0p[i];
481
482 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
483
484 /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
485 row[0] = idx + 1;
486 col[0] = idx + 1;
487 col[1] = idx + 5;
488 val[0] = 1 / Tq0p[i];
489 val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
490 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
491
492 /* fgen[idx+2] = - w + w_s; */
493 row[0] = idx + 2;
494 col[0] = idx + 2;
495 col[1] = idx + 3;
496 val[0] = 0;
497 val[1] = -1;
498 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
499
500 /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
501 row[0] = idx + 3;
502 col[0] = idx;
503 col[1] = idx + 1;
504 col[2] = idx + 3;
505 col[3] = idx + 4;
506 col[4] = idx + 5;
507 val[0] = Iq / M[i];
508 val[1] = Id / M[i];
509 val[2] = D[i] / M[i];
510 val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
511 val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
512 PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
513
514 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
515 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
516 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
517
518 det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
519
520 Zdq_inv[0] = Rs[i] / det;
521 Zdq_inv[1] = Xqp[i] / det;
522 Zdq_inv[2] = -Xdp[i] / det;
523 Zdq_inv[3] = Rs[i] / det;
524
525 dVd_dVr = PetscSinScalar(delta);
526 dVd_dVi = -PetscCosScalar(delta);
527 dVq_dVr = PetscCosScalar(delta);
528 dVq_dVi = PetscSinScalar(delta);
529 dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
530 dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
531
532 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
533 row[0] = idx + 4;
534 col[0] = idx;
535 col[1] = idx + 1;
536 col[2] = idx + 2;
537 val[0] = -Zdq_inv[1];
538 val[1] = -Zdq_inv[0];
539 val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
540 col[3] = idx + 4;
541 col[4] = net_start + 2 * gbus[i];
542 col[5] = net_start + 2 * gbus[i] + 1;
543 val[3] = 1;
544 val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
545 val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
546 PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
547
548 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
549 row[0] = idx + 5;
550 col[0] = idx;
551 col[1] = idx + 1;
552 col[2] = idx + 2;
553 val[0] = -Zdq_inv[3];
554 val[1] = -Zdq_inv[2];
555 val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
556 col[3] = idx + 5;
557 col[4] = net_start + 2 * gbus[i];
558 col[5] = net_start + 2 * gbus[i] + 1;
559 val[3] = 1;
560 val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
561 val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
562 PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
563
564 dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
565 dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
566 dIGr_dId = PetscSinScalar(delta);
567 dIGr_dIq = PetscCosScalar(delta);
568 dIGi_dId = -PetscCosScalar(delta);
569 dIGi_dIq = PetscSinScalar(delta);
570
571 /* fnet[2*gbus[i]] -= IGi; */
572 row[0] = net_start + 2 * gbus[i];
573 col[0] = idx + 2;
574 col[1] = idx + 4;
575 col[2] = idx + 5;
576 val[0] = -dIGi_ddelta;
577 val[1] = -dIGi_dId;
578 val[2] = -dIGi_dIq;
579 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
580
581 /* fnet[2*gbus[i]+1] -= IGr; */
582 row[0] = net_start + 2 * gbus[i] + 1;
583 col[0] = idx + 2;
584 col[1] = idx + 4;
585 col[2] = idx + 5;
586 val[0] = -dIGr_ddelta;
587 val[1] = -dIGr_dId;
588 val[2] = -dIGr_dIq;
589 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
590
591 Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
592
593 /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
594 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
595
596 dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
597
598 row[0] = idx + 6;
599 col[0] = idx + 6;
600 col[1] = idx + 8;
601 val[0] = (KE[i] + dSE_dEfd) / TE[i];
602 val[1] = -1 / TE[i];
603 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
604
605 /* Exciter differential equations */
606
607 /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
608 row[0] = idx + 7;
609 col[0] = idx + 6;
610 col[1] = idx + 7;
611 val[0] = (-KF[i] / TF[i]) / TF[i];
612 val[1] = 1 / TF[i];
613 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
614
615 /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
616 /* Vm = (Vd^2 + Vq^2)^0.5; */
617 dVm_dVd = Vd / Vm;
618 dVm_dVq = Vq / Vm;
619 dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
620 dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
621 row[0] = idx + 8;
622 col[0] = idx + 6;
623 col[1] = idx + 7;
624 col[2] = idx + 8;
625 val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
626 val[1] = -KA[i] / TA[i];
627 val[2] = 1 / TA[i];
628 col[3] = net_start + 2 * gbus[i];
629 col[4] = net_start + 2 * gbus[i] + 1;
630 val[3] = KA[i] * dVm_dVr / TA[i];
631 val[4] = KA[i] * dVm_dVi / TA[i];
632 PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
633 idx = idx + 9;
634 }
635
636 for (i = 0; i < nbus; i++) {
637 PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
638 row[0] = net_start + 2 * i;
639 for (k = 0; k < ncols; k++) {
640 col[k] = net_start + cols[k];
641 val[k] = yvals[k];
642 }
643 PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
644 PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
645
646 PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
647 row[0] = net_start + 2 * i + 1;
648 for (k = 0; k < ncols; k++) {
649 col[k] = net_start + cols[k];
650 val[k] = yvals[k];
651 }
652 PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
653 PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
654 }
655
656 PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
657 PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));
658
659 PetscCall(VecGetArray(user->V0, &v0));
660 for (i = 0; i < nload; i++) {
661 Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
662 Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
663 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
664 Vm2 = Vm * Vm;
665 Vm4 = Vm2 * Vm2;
666 Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
667 PD = QD = 0.0;
668 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
669 for (k = 0; k < ld_nsegsp[i]; k++) {
670 PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
671 dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vr * PetscPowScalar(Vm, ld_betap[k] - 2);
672 dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vi * PetscPowScalar(Vm, ld_betap[k] - 2);
673 }
674 for (k = 0; k < ld_nsegsq[i]; k++) {
675 QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
676 dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vr * PetscPowScalar(Vm, ld_betaq[k] - 2);
677 dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vi * PetscPowScalar(Vm, ld_betaq[k] - 2);
678 }
679
680 /* IDr = (PD*Vr + QD*Vi)/Vm2; */
681 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
682
683 dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
684 dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
685
686 dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
687 dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
688
689 /* fnet[2*lbus[i]] += IDi; */
690 row[0] = net_start + 2 * lbus[i];
691 col[0] = net_start + 2 * lbus[i];
692 col[1] = net_start + 2 * lbus[i] + 1;
693 val[0] = dIDi_dVr;
694 val[1] = dIDi_dVi;
695 PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
696 /* fnet[2*lbus[i]+1] += IDr; */
697 row[0] = net_start + 2 * lbus[i] + 1;
698 col[0] = net_start + 2 * lbus[i];
699 col[1] = net_start + 2 * lbus[i] + 1;
700 val[0] = dIDr_dVr;
701 val[1] = dIDr_dVi;
702 PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
703 }
704 PetscCall(VecRestoreArray(user->V0, &v0));
705
706 PetscCall(VecRestoreArray(Xgen, &xgen));
707 PetscCall(VecRestoreArray(Xnet, &xnet));
708
709 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
710
711 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
712 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
713 PetscFunctionReturn(PETSC_SUCCESS);
714 }
715
716 /*
717 J = [I, 0
718 dg_dx, dg_dy]
719 */
AlgJacobian(SNES snes,Vec X,Mat A,Mat B,PetscCtx ctx)720 PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
721 {
722 Userctx *user = (Userctx *)ctx;
723
724 PetscFunctionBegin;
725 PetscCall(ResidualJacobian(snes, X, A, B, ctx));
726 PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
727 PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
728 PetscFunctionReturn(PETSC_SUCCESS);
729 }
730
731 /*
732 J = [a*I-df_dx, -df_dy
733 dg_dx, dg_dy]
734 */
735
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx * user)736 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
737 {
738 SNES snes;
739 PetscScalar atmp = (PetscScalar)a;
740 PetscInt i, row;
741
742 PetscFunctionBegin;
743 user->t = t;
744
745 PetscCall(TSGetSNES(ts, &snes));
746 PetscCall(ResidualJacobian(snes, X, A, B, user));
747 for (i = 0; i < ngen; i++) {
748 row = 9 * i;
749 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
750 row = 9 * i + 1;
751 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
752 row = 9 * i + 2;
753 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
754 row = 9 * i + 3;
755 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
756 row = 9 * i + 6;
757 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
758 row = 9 * i + 7;
759 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
760 row = 9 * i + 8;
761 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
762 }
763 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
764 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
765 PetscFunctionReturn(PETSC_SUCCESS);
766 }
767
CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,Userctx * user)768 static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
769 {
770 PetscScalar *r;
771 const PetscScalar *u;
772 PetscInt idx;
773 Vec Xgen, Xnet;
774 PetscScalar *xgen;
775 PetscInt i;
776
777 PetscFunctionBegin;
778 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
779 PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
780
781 PetscCall(VecGetArray(Xgen, &xgen));
782
783 PetscCall(VecGetArrayRead(U, &u));
784 PetscCall(VecGetArray(R, &r));
785 r[0] = 0.;
786
787 idx = 0;
788 for (i = 0; i < ngen; i++) {
789 r[0] += PetscPowScalarInt(PetscMax(0., PetscMax(xgen[idx + 3] / (2. * PETSC_PI) - user->freq_u, user->freq_l - xgen[idx + 3] / (2. * PETSC_PI))), user->pow);
790 idx += 9;
791 }
792 PetscCall(VecRestoreArray(R, &r));
793 PetscCall(VecRestoreArrayRead(U, &u));
794 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
795 PetscFunctionReturn(PETSC_SUCCESS);
796 }
797
MonitorUpdateQ(TS ts,PetscInt stepnum,PetscReal time,Vec X,PetscCtx ctx0)798 static PetscErrorCode MonitorUpdateQ(TS ts, PetscInt stepnum, PetscReal time, Vec X, PetscCtx ctx0)
799 {
800 Vec C, *Y;
801 PetscInt Nr;
802 PetscReal h, theta;
803 Userctx *ctx = (Userctx *)ctx0;
804
805 PetscFunctionBegin;
806 theta = 0.5;
807 PetscCall(TSGetStages(ts, &Nr, &Y));
808 PetscCall(TSGetTimeStep(ts, &h));
809 PetscCall(VecDuplicate(ctx->vec_q, &C));
810 /* compute integrals */
811 if (stepnum > 0) {
812 PetscCall(CostIntegrand(ts, time, X, C, ctx));
813 PetscCall(VecAXPY(ctx->vec_q, h * theta, C));
814 PetscCall(CostIntegrand(ts, time + h * theta, Y[0], C, ctx));
815 PetscCall(VecAXPY(ctx->vec_q, h * (1 - theta), C));
816 }
817 PetscCall(VecDestroy(&C));
818 PetscFunctionReturn(PETSC_SUCCESS);
819 }
820
main(int argc,char ** argv)821 int main(int argc, char **argv)
822 {
823 Userctx user;
824 Vec p;
825 PetscScalar *x_ptr;
826 PetscMPIInt size;
827 PetscInt i;
828 KSP ksp;
829 PC pc;
830 PetscInt *idx2;
831 Tao tao;
832 Vec lowerb, upperb;
833
834 PetscFunctionBeginUser;
835 PetscFunctionBeginUser;
836 PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
837 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
838 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
839
840 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 1, &user.vec_q));
841
842 user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
843 user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
844 user.neqs_pgrid = user.neqs_gen + user.neqs_net;
845
846 /* Create indices for differential and algebraic equations */
847 PetscCall(PetscMalloc1(7 * ngen, &idx2));
848 for (i = 0; i < ngen; i++) {
849 idx2[7 * i] = 9 * i;
850 idx2[7 * i + 1] = 9 * i + 1;
851 idx2[7 * i + 2] = 9 * i + 2;
852 idx2[7 * i + 3] = 9 * i + 3;
853 idx2[7 * i + 4] = 9 * i + 6;
854 idx2[7 * i + 5] = 9 * i + 7;
855 idx2[7 * i + 6] = 9 * i + 8;
856 }
857 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
858 PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
859 PetscCall(PetscFree(idx2));
860
861 /* Set run time options */
862 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
863 {
864 user.tfaulton = 1.0;
865 user.tfaultoff = 1.2;
866 user.Rfault = 0.0001;
867 user.faultbus = 8;
868 PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
869 PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
870 PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
871 user.t0 = 0.0;
872 user.tmax = 1.5;
873 PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
874 PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
875 user.freq_u = 61.0;
876 user.freq_l = 59.0;
877 user.pow = 2;
878 PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
879 PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
880 PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
881 }
882 PetscOptionsEnd();
883
884 /* Create DMs for generator and network subsystems */
885 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
886 PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
887 PetscCall(DMSetFromOptions(user.dmgen));
888 PetscCall(DMSetUp(user.dmgen));
889 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
890 PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
891 PetscCall(DMSetFromOptions(user.dmnet));
892 PetscCall(DMSetUp(user.dmnet));
893 /* Create a composite DM packer and add the two DMs */
894 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
895 PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
896 PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
897 PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));
898
899 /* Create TAO solver and set desired solution method */
900 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
901 PetscCall(TaoSetType(tao, TAOBLMVM));
902 /*
903 Optimization starts
904 */
905 /* Set initial solution guess */
906 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
907 PetscCall(VecGetArray(p, &x_ptr));
908 x_ptr[0] = PG[0];
909 x_ptr[1] = PG[1];
910 x_ptr[2] = PG[2];
911 PetscCall(VecRestoreArray(p, &x_ptr));
912
913 PetscCall(TaoSetSolution(tao, p));
914 /* Set routine for function and gradient evaluation */
915 PetscCall(TaoSetObjective(tao, FormFunction, (void *)&user));
916 PetscCall(TaoSetGradient(tao, NULL, TaoDefaultComputeGradient, (void *)&user));
917
918 /* Set bounds for the optimization */
919 PetscCall(VecDuplicate(p, &lowerb));
920 PetscCall(VecDuplicate(p, &upperb));
921 PetscCall(VecGetArray(lowerb, &x_ptr));
922 x_ptr[0] = 0.5;
923 x_ptr[1] = 0.5;
924 x_ptr[2] = 0.5;
925 PetscCall(VecRestoreArray(lowerb, &x_ptr));
926 PetscCall(VecGetArray(upperb, &x_ptr));
927 x_ptr[0] = 2.0;
928 x_ptr[1] = 2.0;
929 x_ptr[2] = 2.0;
930 PetscCall(VecRestoreArray(upperb, &x_ptr));
931 PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));
932
933 /* Check for any TAO command line options */
934 PetscCall(TaoSetFromOptions(tao));
935 PetscCall(TaoGetKSP(tao, &ksp));
936 if (ksp) {
937 PetscCall(KSPGetPC(ksp, &pc));
938 PetscCall(PCSetType(pc, PCNONE));
939 }
940
941 /* SOLVE THE APPLICATION */
942 PetscCall(TaoSolve(tao));
943
944 PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
945 /* Free TAO data structures */
946 PetscCall(TaoDestroy(&tao));
947 PetscCall(VecDestroy(&user.vec_q));
948 PetscCall(VecDestroy(&lowerb));
949 PetscCall(VecDestroy(&upperb));
950 PetscCall(VecDestroy(&p));
951 PetscCall(DMDestroy(&user.dmgen));
952 PetscCall(DMDestroy(&user.dmnet));
953 PetscCall(DMDestroy(&user.dmpgrid));
954 PetscCall(ISDestroy(&user.is_diff));
955 PetscCall(ISDestroy(&user.is_alg));
956 PetscCall(PetscFinalize());
957 return 0;
958 }
959
960 /* ------------------------------------------------------------------ */
961 /*
962 FormFunction - Evaluates the function and corresponding gradient.
963
964 Input Parameters:
965 tao - the Tao context
966 X - the input vector
967 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
968
969 Output Parameters:
970 f - the newly evaluated function
971 */
FormFunction(Tao tao,Vec P,PetscReal * f,PetscCtx ctx0)972 PetscErrorCode FormFunction(Tao tao, Vec P, PetscReal *f, PetscCtx ctx0)
973 {
974 TS ts;
975 SNES snes_alg;
976 Userctx *ctx = (Userctx *)ctx0;
977 Vec X;
978 Mat J;
979 /* sensitivity context */
980 PetscScalar *x_ptr;
981 PetscViewer Xview, Ybusview;
982 Vec F_alg;
983 Vec Xdot;
984 PetscInt row_loc, col_loc;
985 PetscScalar val;
986
987 PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
988 PG[0] = x_ptr[0];
989 PG[1] = x_ptr[1];
990 PG[2] = x_ptr[2];
991 PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
992
993 ctx->stepnum = 0;
994
995 PetscCall(VecZeroEntries(ctx->vec_q));
996
997 /* Read initial voltage vector and Ybus */
998 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
999 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));
1000
1001 PetscCall(VecCreate(PETSC_COMM_WORLD, &ctx->V0));
1002 PetscCall(VecSetSizes(ctx->V0, PETSC_DECIDE, ctx->neqs_net));
1003 PetscCall(VecLoad(ctx->V0, Xview));
1004
1005 PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx->Ybus));
1006 PetscCall(MatSetSizes(ctx->Ybus, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_net, ctx->neqs_net));
1007 PetscCall(MatSetType(ctx->Ybus, MATBAIJ));
1008 /* PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
1009 PetscCall(MatLoad(ctx->Ybus, Ybusview));
1010
1011 PetscCall(PetscViewerDestroy(&Xview));
1012 PetscCall(PetscViewerDestroy(&Ybusview));
1013
1014 PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));
1015
1016 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
1017 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_pgrid, ctx->neqs_pgrid));
1018 PetscCall(MatSetFromOptions(J));
1019 PetscCall(PreallocateJacobian(J, ctx));
1020
1021 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1022 Create timestepping solver context
1023 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1024 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1025 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1026 PetscCall(TSSetType(ts, TSCN));
1027 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, ctx));
1028 PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobianFn *)IJacobian, ctx));
1029 PetscCall(TSSetApplicationContext(ts, ctx));
1030
1031 PetscCall(TSMonitorSet(ts, MonitorUpdateQ, ctx, NULL));
1032 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1033 Set initial conditions
1034 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1035 PetscCall(SetInitialGuess(X, ctx));
1036
1037 PetscCall(VecDuplicate(X, &F_alg));
1038 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1039 PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
1040 PetscCall(MatZeroEntries(J));
1041 PetscCall(SNESSetJacobian(snes_alg, J, J, AlgJacobian, ctx));
1042 PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1043 PetscCall(SNESSetFromOptions(snes_alg));
1044 ctx->alg_flg = PETSC_TRUE;
1045 /* Solve the algebraic equations */
1046 PetscCall(SNESSolve(snes_alg, NULL, X));
1047
1048 /* Just to set up the Jacobian structure */
1049 PetscCall(VecDuplicate(X, &Xdot));
1050 PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, J, J, ctx));
1051 PetscCall(VecDestroy(&Xdot));
1052
1053 ctx->stepnum++;
1054
1055 PetscCall(TSSetTimeStep(ts, 0.01));
1056 PetscCall(TSSetMaxTime(ts, ctx->tmax));
1057 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1058 PetscCall(TSSetFromOptions(ts));
1059
1060 /* Prefault period */
1061 ctx->alg_flg = PETSC_FALSE;
1062 PetscCall(TSSetTime(ts, 0.0));
1063 PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
1064 PetscCall(TSSolve(ts, X));
1065
1066 /* Create the nonlinear solver for solving the algebraic system */
1067 /* Note that although the algebraic system needs to be solved only for
1068 Idq and V, we reuse the entire system including xgen. The xgen
1069 variables are held constant by setting their residuals to 0 and
1070 putting a 1 on the Jacobian diagonal for xgen rows
1071 */
1072 PetscCall(MatZeroEntries(J));
1073
1074 /* Apply disturbance - resistive fault at ctx->faultbus */
1075 /* This is done by adding shunt conductance to the diagonal location
1076 in the Ybus matrix */
1077 row_loc = 2 * ctx->faultbus;
1078 col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1079 val = 1 / ctx->Rfault;
1080 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1081 row_loc = 2 * ctx->faultbus + 1;
1082 col_loc = 2 * ctx->faultbus; /* Location for G */
1083 val = 1 / ctx->Rfault;
1084 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1085
1086 PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1087 PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1088
1089 ctx->alg_flg = PETSC_TRUE;
1090 /* Solve the algebraic equations */
1091 PetscCall(SNESSolve(snes_alg, NULL, X));
1092
1093 ctx->stepnum++;
1094
1095 /* Disturbance period */
1096 ctx->alg_flg = PETSC_FALSE;
1097 PetscCall(TSSetTime(ts, ctx->tfaulton));
1098 PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
1099 PetscCall(TSSolve(ts, X));
1100
1101 /* Remove the fault */
1102 row_loc = 2 * ctx->faultbus;
1103 col_loc = 2 * ctx->faultbus + 1;
1104 val = -1 / ctx->Rfault;
1105 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1106 row_loc = 2 * ctx->faultbus + 1;
1107 col_loc = 2 * ctx->faultbus;
1108 val = -1 / ctx->Rfault;
1109 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1110
1111 PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1112 PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1113
1114 PetscCall(MatZeroEntries(J));
1115
1116 ctx->alg_flg = PETSC_TRUE;
1117
1118 /* Solve the algebraic equations */
1119 PetscCall(SNESSolve(snes_alg, NULL, X));
1120
1121 ctx->stepnum++;
1122
1123 /* Post-disturbance period */
1124 ctx->alg_flg = PETSC_TRUE;
1125 PetscCall(TSSetTime(ts, ctx->tfaultoff));
1126 PetscCall(TSSetMaxTime(ts, ctx->tmax));
1127 PetscCall(TSSolve(ts, X));
1128
1129 PetscCall(VecGetArray(ctx->vec_q, &x_ptr));
1130 *f = x_ptr[0];
1131 PetscCall(VecRestoreArray(ctx->vec_q, &x_ptr));
1132
1133 PetscCall(MatDestroy(&ctx->Ybus));
1134 PetscCall(VecDestroy(&ctx->V0));
1135 PetscCall(SNESDestroy(&snes_alg));
1136 PetscCall(VecDestroy(&F_alg));
1137 PetscCall(MatDestroy(&J));
1138 PetscCall(VecDestroy(&X));
1139 PetscCall(TSDestroy(&ts));
1140
1141 return 0;
1142 }
1143
1144 /*TEST
1145
1146 build:
1147 requires: double !complex !defined(USE_64BIT_INDICES)
1148
1149 test:
1150 args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1151 localrunfiles: petscoptions X.bin Ybus.bin
1152
1153 TEST*/
1154