1 static char help[] = "Application of adjoint sensitivity analysis for power grid stability analysis of WECC 9 bus system.\n\
2 This example is based on the 9-bus (node) example given in the book Power\n\
3 Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\
4 The power grid in this example consists of 9 buses (nodes), 3 generators,\n\
5 3 loads, and 9 transmission lines. The network equations are written\n\
6 in current balance form using rectangular coordinates.\n\n";
7
8 /*
9 This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
10 The objectivie is to find optimal parameter PG for each generator to minizie the frequency violations due to faults.
11 The problem features discontinuities and a cost function in integral form.
12 The gradient is computed with the discrete adjoint of an implicit theta method, see ex9busadj.c for details.
13 */
14
15 #include <petsctao.h>
16 #include <petscts.h>
17 #include <petscdm.h>
18 #include <petscdmda.h>
19 #include <petscdmcomposite.h>
20 #include <petsctime.h>
21
22 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
23
24 #define freq 60
25 #define w_s (2 * PETSC_PI * freq)
26
27 /* Sizes and indices */
28 const PetscInt nbus = 9; /* Number of network buses */
29 const PetscInt ngen = 3; /* Number of generators */
30 const PetscInt nload = 3; /* Number of loads */
31 const PetscInt gbus[3] = {0, 1, 2}; /* Buses at which generators are incident */
32 const PetscInt lbus[3] = {4, 5, 7}; /* Buses at which loads are incident */
33
34 /* Generator real and reactive powers (found via loadflow) */
35 PetscScalar PG[3] = {0.69, 1.59, 0.69};
36 /* PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};*/
37
38 const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
39 /* Generator constants */
40 const PetscScalar H[3] = {23.64, 6.4, 3.01}; /* Inertia constant */
41 const PetscScalar Rs[3] = {0.0, 0.0, 0.0}; /* Stator Resistance */
42 const PetscScalar Xd[3] = {0.146, 0.8958, 1.3125}; /* d-axis reactance */
43 const PetscScalar Xdp[3] = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
44 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 */
45 const PetscScalar Xqp[3] = {0.0969, 0.1969, 0.25}; /* q-axis transient reactance */
46 const PetscScalar Td0p[3] = {8.96, 6.0, 5.89}; /* d-axis open circuit time constant */
47 const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6}; /* q-axis open circuit time constant */
48 PetscScalar M[3]; /* M = 2*H/w_s */
49 PetscScalar D[3]; /* D = 0.1*M */
50
51 PetscScalar TM[3]; /* Mechanical Torque */
52 /* Exciter system constants */
53 const PetscScalar KA[3] = {20.0, 20.0, 20.0}; /* Voltage regulartor gain constant */
54 const PetscScalar TA[3] = {0.2, 0.2, 0.2}; /* Voltage regulator time constant */
55 const PetscScalar KE[3] = {1.0, 1.0, 1.0}; /* Exciter gain constant */
56 const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
57 const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
58 const PetscScalar TF[3] = {0.35, 0.35, 0.35}; /* Feedback stabilizer time constant */
59 const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
60 const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */
61
62 PetscScalar Vref[3];
63 /* Load constants
64 We use a composite load model that describes the load and reactive powers at each time instant as follows
65 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
66 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
67 where
68 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
69 ld_alphap,ld_alphap - Percentage contribution (weights) or loads
70 P_D0 - Real power load
71 Q_D0 - Reactive power load
72 V_m(t) - Voltage magnitude at time t
73 V_m0 - Voltage magnitude at t = 0
74 ld_betap, ld_betaq - exponents describing the load model for real and reactive part
75
76 Note: All loads have the same characteristic currently.
77 */
78 const PetscScalar PD0[3] = {1.25, 0.9, 1.0};
79 const PetscScalar QD0[3] = {0.5, 0.3, 0.35};
80 const PetscInt ld_nsegsp[3] = {3, 3, 3};
81 const PetscScalar ld_alphap[3] = {1.0, 0.0, 0.0};
82 const PetscScalar ld_betap[3] = {2.0, 1.0, 0.0};
83 const PetscInt ld_nsegsq[3] = {3, 3, 3};
84 const PetscScalar ld_alphaq[3] = {1.0, 0.0, 0.0};
85 const PetscScalar ld_betaq[3] = {2.0, 1.0, 0.0};
86
87 typedef struct {
88 DM dmgen, dmnet; /* DMs to manage generator and network subsystem */
89 DM dmpgrid; /* Composite DM to manage the entire power grid */
90 Mat Ybus; /* Network admittance matrix */
91 Vec V0; /* Initial voltage vector (Power flow solution) */
92 PetscReal tfaulton, tfaultoff; /* Fault on and off times */
93 PetscInt faultbus; /* Fault bus */
94 PetscScalar Rfault;
95 PetscReal t0, tmax;
96 PetscInt neqs_gen, neqs_net, neqs_pgrid;
97 Mat Sol; /* Matrix to save solution at each time step */
98 PetscInt stepnum;
99 PetscBool alg_flg;
100 PetscReal t;
101 IS is_diff; /* indices for differential equations */
102 IS is_alg; /* indices for algebraic equations */
103 PetscReal freq_u, freq_l; /* upper and lower frequency limit */
104 PetscInt pow; /* power coefficient used in the cost function */
105 PetscBool jacp_flg;
106 Mat J, Jacp;
107 Mat DRDU, DRDP;
108 } Userctx;
109
110 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar * Fr,PetscScalar * Fi)111 PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
112 {
113 PetscFunctionBegin;
114 *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
115 *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
116 PetscFunctionReturn(PETSC_SUCCESS);
117 }
118
119 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar * Fd,PetscScalar * Fq)120 PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
121 {
122 PetscFunctionBegin;
123 *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
124 *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
128 /* Saves the solution at each time to a matrix */
SaveSolution(TS ts)129 PetscErrorCode SaveSolution(TS ts)
130 {
131 Userctx *user;
132 Vec X;
133 PetscScalar *mat;
134 const PetscScalar *x;
135 PetscInt idx;
136 PetscReal t;
137
138 PetscFunctionBegin;
139 PetscCall(TSGetApplicationContext(ts, &user));
140 PetscCall(TSGetTime(ts, &t));
141 PetscCall(TSGetSolution(ts, &X));
142 idx = user->stepnum * (user->neqs_pgrid + 1);
143 PetscCall(MatDenseGetArray(user->Sol, &mat));
144 PetscCall(VecGetArrayRead(X, &x));
145 mat[idx] = t;
146 PetscCall(PetscArraycpy(mat + idx + 1, x, user->neqs_pgrid));
147 PetscCall(MatDenseRestoreArray(user->Sol, &mat));
148 PetscCall(VecRestoreArrayRead(X, &x));
149 user->stepnum++;
150 PetscFunctionReturn(PETSC_SUCCESS);
151 }
152
SetInitialGuess(Vec X,Userctx * user)153 PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
154 {
155 Vec Xgen, Xnet;
156 PetscScalar *xgen, *xnet;
157 PetscInt i, idx = 0;
158 PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
159 PetscScalar Eqp, Edp, delta;
160 PetscScalar Efd, RF, VR; /* Exciter variables */
161 PetscScalar Id, Iq; /* Generator dq axis currents */
162 PetscScalar theta, Vd, Vq, SE;
163
164 PetscFunctionBegin;
165 M[0] = 2 * H[0] / w_s;
166 M[1] = 2 * H[1] / w_s;
167 M[2] = 2 * H[2] / w_s;
168 D[0] = 0.1 * M[0];
169 D[1] = 0.1 * M[1];
170 D[2] = 0.1 * M[2];
171
172 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
173
174 /* Network subsystem initialization */
175 PetscCall(VecCopy(user->V0, Xnet));
176
177 /* Generator subsystem initialization */
178 PetscCall(VecGetArray(Xgen, &xgen));
179 PetscCall(VecGetArray(Xnet, &xnet));
180
181 for (i = 0; i < ngen; i++) {
182 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
183 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
184 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
185 Vm2 = Vm * Vm;
186 IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
187 IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;
188
189 delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
190
191 theta = PETSC_PI / 2.0 - delta;
192
193 Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
194 Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
195
196 Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
197 Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
198
199 Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
200 Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
201
202 TM[i] = PG[i];
203
204 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
205 xgen[idx] = Eqp;
206 xgen[idx + 1] = Edp;
207 xgen[idx + 2] = delta;
208 xgen[idx + 3] = w_s;
209
210 idx = idx + 4;
211
212 xgen[idx] = Id;
213 xgen[idx + 1] = Iq;
214
215 idx = idx + 2;
216
217 /* Exciter */
218 Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
219 SE = k1[i] * PetscExpScalar(k2[i] * Efd);
220 VR = KE[i] * Efd + SE;
221 RF = KF[i] * Efd / TF[i];
222
223 xgen[idx] = Efd;
224 xgen[idx + 1] = RF;
225 xgen[idx + 2] = VR;
226
227 Vref[i] = Vm + (VR / KA[i]);
228
229 idx = idx + 3;
230 }
231
232 PetscCall(VecRestoreArray(Xgen, &xgen));
233 PetscCall(VecRestoreArray(Xnet, &xnet));
234
235 /* PetscCall(VecView(Xgen,0)); */
236 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
237 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
238 PetscFunctionReturn(PETSC_SUCCESS);
239 }
240
InitialGuess(Vec X,Userctx * user,const PetscScalar PGv[])241 PetscErrorCode InitialGuess(Vec X, Userctx *user, const PetscScalar PGv[])
242 {
243 Vec Xgen, Xnet;
244 PetscScalar *xgen, *xnet;
245 PetscInt i, idx = 0;
246 PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
247 PetscScalar Eqp, Edp, delta;
248 PetscScalar Efd, RF, VR; /* Exciter variables */
249 PetscScalar Id, Iq; /* Generator dq axis currents */
250 PetscScalar theta, Vd, Vq, SE;
251
252 PetscFunctionBegin;
253 M[0] = 2 * H[0] / w_s;
254 M[1] = 2 * H[1] / w_s;
255 M[2] = 2 * H[2] / w_s;
256 D[0] = 0.1 * M[0];
257 D[1] = 0.1 * M[1];
258 D[2] = 0.1 * M[2];
259
260 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
261
262 /* Network subsystem initialization */
263 PetscCall(VecCopy(user->V0, Xnet));
264
265 /* Generator subsystem initialization */
266 PetscCall(VecGetArray(Xgen, &xgen));
267 PetscCall(VecGetArray(Xnet, &xnet));
268
269 for (i = 0; i < ngen; i++) {
270 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
271 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
272 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
273 Vm2 = Vm * Vm;
274 IGr = (Vr * PGv[i] + Vi * QG[i]) / Vm2;
275 IGi = (Vi * PGv[i] - Vr * QG[i]) / Vm2;
276
277 delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
278
279 theta = PETSC_PI / 2.0 - delta;
280
281 Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
282 Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
283
284 Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
285 Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
286
287 Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
288 Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
289
290 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
291 xgen[idx] = Eqp;
292 xgen[idx + 1] = Edp;
293 xgen[idx + 2] = delta;
294 xgen[idx + 3] = w_s;
295
296 idx = idx + 4;
297
298 xgen[idx] = Id;
299 xgen[idx + 1] = Iq;
300
301 idx = idx + 2;
302
303 /* Exciter */
304 Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
305 SE = k1[i] * PetscExpScalar(k2[i] * Efd);
306 VR = KE[i] * Efd + SE;
307 RF = KF[i] * Efd / TF[i];
308
309 xgen[idx] = Efd;
310 xgen[idx + 1] = RF;
311 xgen[idx + 2] = VR;
312
313 idx = idx + 3;
314 }
315
316 PetscCall(VecRestoreArray(Xgen, &xgen));
317 PetscCall(VecRestoreArray(Xnet, &xnet));
318
319 /* PetscCall(VecView(Xgen,0)); */
320 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
321 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
322 PetscFunctionReturn(PETSC_SUCCESS);
323 }
324
DICDPFiniteDifference(Vec X,Vec * DICDP,Userctx * user)325 PetscErrorCode DICDPFiniteDifference(Vec X, Vec *DICDP, Userctx *user)
326 {
327 Vec Y;
328 PetscScalar PGv[3], eps;
329 PetscInt i, j;
330
331 PetscFunctionBegin;
332 eps = 1.e-7;
333 PetscCall(VecDuplicate(X, &Y));
334
335 for (i = 0; i < ngen; i++) {
336 for (j = 0; j < 3; j++) PGv[j] = PG[j];
337 PGv[i] = PG[i] + eps;
338 PetscCall(InitialGuess(Y, user, PGv));
339 PetscCall(InitialGuess(X, user, PG));
340
341 PetscCall(VecAXPY(Y, -1.0, X));
342 PetscCall(VecScale(Y, 1. / eps));
343 PetscCall(VecCopy(Y, DICDP[i]));
344 }
345 PetscCall(VecDestroy(&Y));
346 PetscFunctionReturn(PETSC_SUCCESS);
347 }
348
349 /* Computes F = [-f(x,y);g(x,y)] */
ResidualFunction(SNES snes,Vec X,Vec F,Userctx * user)350 PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
351 {
352 Vec Xgen, Xnet, Fgen, Fnet;
353 PetscScalar *xgen, *xnet, *fgen, *fnet;
354 PetscInt i, idx = 0;
355 PetscScalar Vr, Vi, Vm, Vm2;
356 PetscScalar Eqp, Edp, delta, w; /* Generator variables */
357 PetscScalar Efd, RF, VR; /* Exciter variables */
358 PetscScalar Id, Iq; /* Generator dq axis currents */
359 PetscScalar Vd, Vq, SE;
360 PetscScalar IGr, IGi, IDr, IDi;
361 PetscScalar Zdq_inv[4], det;
362 PetscScalar PD, QD, Vm0, *v0;
363 PetscInt k;
364
365 PetscFunctionBegin;
366 PetscCall(VecZeroEntries(F));
367 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
368 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
369 PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
370 PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));
371
372 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
373 The generator current injection, IG, and load current injection, ID are added later
374 */
375 /* Note that the values in Ybus are stored assuming the imaginary current balance
376 equation is ordered first followed by real current balance equation for each bus.
377 Thus imaginary current contribution goes in location 2*i, and
378 real current contribution in 2*i+1
379 */
380 PetscCall(MatMult(user->Ybus, Xnet, Fnet));
381
382 PetscCall(VecGetArray(Xgen, &xgen));
383 PetscCall(VecGetArray(Xnet, &xnet));
384 PetscCall(VecGetArray(Fgen, &fgen));
385 PetscCall(VecGetArray(Fnet, &fnet));
386
387 /* Generator subsystem */
388 for (i = 0; i < ngen; i++) {
389 Eqp = xgen[idx];
390 Edp = xgen[idx + 1];
391 delta = xgen[idx + 2];
392 w = xgen[idx + 3];
393 Id = xgen[idx + 4];
394 Iq = xgen[idx + 5];
395 Efd = xgen[idx + 6];
396 RF = xgen[idx + 7];
397 VR = xgen[idx + 8];
398
399 /* Generator differential equations */
400 fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
401 fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
402 fgen[idx + 2] = -w + w_s;
403 fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
404
405 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
406 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
407
408 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
409 /* Algebraic equations for stator currents */
410 det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
411
412 Zdq_inv[0] = Rs[i] / det;
413 Zdq_inv[1] = Xqp[i] / det;
414 Zdq_inv[2] = -Xdp[i] / det;
415 Zdq_inv[3] = Rs[i] / det;
416
417 fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
418 fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
419
420 /* Add generator current injection to network */
421 PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
422
423 fnet[2 * gbus[i]] -= IGi;
424 fnet[2 * gbus[i] + 1] -= IGr;
425
426 Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
427
428 SE = k1[i] * PetscExpScalar(k2[i] * Efd);
429
430 /* Exciter differential equations */
431 fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
432 fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
433 fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
434
435 idx = idx + 9;
436 }
437
438 PetscCall(VecGetArray(user->V0, &v0));
439 for (i = 0; i < nload; i++) {
440 Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
441 Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
442 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
443 Vm2 = Vm * Vm;
444 Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
445 PD = QD = 0.0;
446 for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
447 for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
448
449 /* Load currents */
450 IDr = (PD * Vr + QD * Vi) / Vm2;
451 IDi = (-QD * Vr + PD * Vi) / Vm2;
452
453 fnet[2 * lbus[i]] += IDi;
454 fnet[2 * lbus[i] + 1] += IDr;
455 }
456 PetscCall(VecRestoreArray(user->V0, &v0));
457
458 PetscCall(VecRestoreArray(Xgen, &xgen));
459 PetscCall(VecRestoreArray(Xnet, &xnet));
460 PetscCall(VecRestoreArray(Fgen, &fgen));
461 PetscCall(VecRestoreArray(Fnet, &fnet));
462
463 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
464 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
465 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
466 PetscFunctionReturn(PETSC_SUCCESS);
467 }
468
469 /* \dot{x} - f(x,y)
470 g(x,y) = 0
471 */
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx * user)472 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
473 {
474 SNES snes;
475 PetscScalar *f;
476 const PetscScalar *xdot;
477 PetscInt i;
478
479 PetscFunctionBegin;
480 user->t = t;
481
482 PetscCall(TSGetSNES(ts, &snes));
483 PetscCall(ResidualFunction(snes, X, F, user));
484 PetscCall(VecGetArray(F, &f));
485 PetscCall(VecGetArrayRead(Xdot, &xdot));
486 for (i = 0; i < ngen; i++) {
487 f[9 * i] += xdot[9 * i];
488 f[9 * i + 1] += xdot[9 * i + 1];
489 f[9 * i + 2] += xdot[9 * i + 2];
490 f[9 * i + 3] += xdot[9 * i + 3];
491 f[9 * i + 6] += xdot[9 * i + 6];
492 f[9 * i + 7] += xdot[9 * i + 7];
493 f[9 * i + 8] += xdot[9 * i + 8];
494 }
495 PetscCall(VecRestoreArray(F, &f));
496 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
497 PetscFunctionReturn(PETSC_SUCCESS);
498 }
499
500 /* This function is used for solving the algebraic system only during fault on and
501 off times. It computes the entire F and then zeros out the part corresponding to
502 differential equations
503 F = [0;g(y)];
504 */
AlgFunction(SNES snes,Vec X,Vec F,PetscCtx ctx)505 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, PetscCtx ctx)
506 {
507 Userctx *user = (Userctx *)ctx;
508 PetscScalar *f;
509 PetscInt i;
510
511 PetscFunctionBegin;
512 PetscCall(ResidualFunction(snes, X, F, user));
513 PetscCall(VecGetArray(F, &f));
514 for (i = 0; i < ngen; i++) {
515 f[9 * i] = 0;
516 f[9 * i + 1] = 0;
517 f[9 * i + 2] = 0;
518 f[9 * i + 3] = 0;
519 f[9 * i + 6] = 0;
520 f[9 * i + 7] = 0;
521 f[9 * i + 8] = 0;
522 }
523 PetscCall(VecRestoreArray(F, &f));
524 PetscFunctionReturn(PETSC_SUCCESS);
525 }
526
PreallocateJacobian(Mat J,Userctx * user)527 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
528 {
529 PetscInt *d_nnz;
530 PetscInt i, idx = 0, start = 0;
531 PetscInt ncols;
532
533 PetscFunctionBegin;
534 PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
535 for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
536 /* Generator subsystem */
537 for (i = 0; i < ngen; i++) {
538 d_nnz[idx] += 3;
539 d_nnz[idx + 1] += 2;
540 d_nnz[idx + 2] += 2;
541 d_nnz[idx + 3] += 5;
542 d_nnz[idx + 4] += 6;
543 d_nnz[idx + 5] += 6;
544
545 d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
546 d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
547
548 d_nnz[idx + 6] += 2;
549 d_nnz[idx + 7] += 2;
550 d_nnz[idx + 8] += 5;
551
552 idx = idx + 9;
553 }
554
555 start = user->neqs_gen;
556 for (i = 0; i < nbus; i++) {
557 PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
558 d_nnz[start + 2 * i] += ncols;
559 d_nnz[start + 2 * i + 1] += ncols;
560 PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
561 }
562
563 PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
564 PetscCall(PetscFree(d_nnz));
565 PetscFunctionReturn(PETSC_SUCCESS);
566 }
567
568 /*
569 J = [-df_dx, -df_dy
570 dg_dx, dg_dy]
571 */
ResidualJacobian(SNES snes,Vec X,Mat J,Mat B,PetscCtx ctx)572 PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, PetscCtx ctx)
573 {
574 Userctx *user = (Userctx *)ctx;
575 Vec Xgen, Xnet;
576 PetscScalar *xgen, *xnet;
577 PetscInt i, idx = 0;
578 PetscScalar Vr, Vi, Vm, Vm2;
579 PetscScalar Eqp, Edp, delta; /* Generator variables */
580 PetscScalar Efd; /* Exciter variables */
581 PetscScalar Id, Iq; /* Generator dq axis currents */
582 PetscScalar Vd, Vq;
583 PetscScalar val[10];
584 PetscInt row[2], col[10];
585 PetscInt net_start = user->neqs_gen;
586 PetscInt ncols;
587 const PetscInt *cols;
588 const PetscScalar *yvals;
589 PetscInt k;
590 PetscScalar Zdq_inv[4], det;
591 PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
592 PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
593 PetscScalar dSE_dEfd;
594 PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
595 PetscScalar PD, QD, Vm0, *v0, Vm4;
596 PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
597 PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
598
599 PetscFunctionBegin;
600 PetscCall(MatZeroEntries(B));
601 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
602 PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
603
604 PetscCall(VecGetArray(Xgen, &xgen));
605 PetscCall(VecGetArray(Xnet, &xnet));
606
607 /* Generator subsystem */
608 for (i = 0; i < ngen; i++) {
609 Eqp = xgen[idx];
610 Edp = xgen[idx + 1];
611 delta = xgen[idx + 2];
612 Id = xgen[idx + 4];
613 Iq = xgen[idx + 5];
614 Efd = xgen[idx + 6];
615
616 /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
617 row[0] = idx;
618 col[0] = idx;
619 col[1] = idx + 4;
620 col[2] = idx + 6;
621 val[0] = 1 / Td0p[i];
622 val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
623 val[2] = -1 / Td0p[i];
624
625 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
626
627 /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
628 row[0] = idx + 1;
629 col[0] = idx + 1;
630 col[1] = idx + 5;
631 val[0] = 1 / Tq0p[i];
632 val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
633 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
634
635 /* fgen[idx+2] = - w + w_s; */
636 row[0] = idx + 2;
637 col[0] = idx + 2;
638 col[1] = idx + 3;
639 val[0] = 0;
640 val[1] = -1;
641 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
642
643 /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
644 row[0] = idx + 3;
645 col[0] = idx;
646 col[1] = idx + 1;
647 col[2] = idx + 3;
648 col[3] = idx + 4;
649 col[4] = idx + 5;
650 val[0] = Iq / M[i];
651 val[1] = Id / M[i];
652 val[2] = D[i] / M[i];
653 val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
654 val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
655 PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
656
657 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
658 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
659 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
660
661 det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
662
663 Zdq_inv[0] = Rs[i] / det;
664 Zdq_inv[1] = Xqp[i] / det;
665 Zdq_inv[2] = -Xdp[i] / det;
666 Zdq_inv[3] = Rs[i] / det;
667
668 dVd_dVr = PetscSinScalar(delta);
669 dVd_dVi = -PetscCosScalar(delta);
670 dVq_dVr = PetscCosScalar(delta);
671 dVq_dVi = PetscSinScalar(delta);
672 dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
673 dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
674
675 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
676 row[0] = idx + 4;
677 col[0] = idx;
678 col[1] = idx + 1;
679 col[2] = idx + 2;
680 val[0] = -Zdq_inv[1];
681 val[1] = -Zdq_inv[0];
682 val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
683 col[3] = idx + 4;
684 col[4] = net_start + 2 * gbus[i];
685 col[5] = net_start + 2 * gbus[i] + 1;
686 val[3] = 1;
687 val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
688 val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
689 PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
690
691 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
692 row[0] = idx + 5;
693 col[0] = idx;
694 col[1] = idx + 1;
695 col[2] = idx + 2;
696 val[0] = -Zdq_inv[3];
697 val[1] = -Zdq_inv[2];
698 val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
699 col[3] = idx + 5;
700 col[4] = net_start + 2 * gbus[i];
701 col[5] = net_start + 2 * gbus[i] + 1;
702 val[3] = 1;
703 val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
704 val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
705 PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
706
707 dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
708 dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
709 dIGr_dId = PetscSinScalar(delta);
710 dIGr_dIq = PetscCosScalar(delta);
711 dIGi_dId = -PetscCosScalar(delta);
712 dIGi_dIq = PetscSinScalar(delta);
713
714 /* fnet[2*gbus[i]] -= IGi; */
715 row[0] = net_start + 2 * gbus[i];
716 col[0] = idx + 2;
717 col[1] = idx + 4;
718 col[2] = idx + 5;
719 val[0] = -dIGi_ddelta;
720 val[1] = -dIGi_dId;
721 val[2] = -dIGi_dIq;
722 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
723
724 /* fnet[2*gbus[i]+1] -= IGr; */
725 row[0] = net_start + 2 * gbus[i] + 1;
726 col[0] = idx + 2;
727 col[1] = idx + 4;
728 col[2] = idx + 5;
729 val[0] = -dIGr_ddelta;
730 val[1] = -dIGr_dId;
731 val[2] = -dIGr_dIq;
732 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
733
734 Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
735
736 /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
737 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
738 dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
739
740 row[0] = idx + 6;
741 col[0] = idx + 6;
742 col[1] = idx + 8;
743 val[0] = (KE[i] + dSE_dEfd) / TE[i];
744 val[1] = -1 / TE[i];
745 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
746
747 /* Exciter differential equations */
748
749 /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
750 row[0] = idx + 7;
751 col[0] = idx + 6;
752 col[1] = idx + 7;
753 val[0] = (-KF[i] / TF[i]) / TF[i];
754 val[1] = 1 / TF[i];
755 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
756
757 /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
758 /* Vm = (Vd^2 + Vq^2)^0.5; */
759 dVm_dVd = Vd / Vm;
760 dVm_dVq = Vq / Vm;
761 dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
762 dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
763 row[0] = idx + 8;
764 col[0] = idx + 6;
765 col[1] = idx + 7;
766 col[2] = idx + 8;
767 val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
768 val[1] = -KA[i] / TA[i];
769 val[2] = 1 / TA[i];
770 col[3] = net_start + 2 * gbus[i];
771 col[4] = net_start + 2 * gbus[i] + 1;
772 val[3] = KA[i] * dVm_dVr / TA[i];
773 val[4] = KA[i] * dVm_dVi / TA[i];
774 PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
775 idx = idx + 9;
776 }
777
778 for (i = 0; i < nbus; i++) {
779 PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
780 row[0] = net_start + 2 * i;
781 for (k = 0; k < ncols; k++) {
782 col[k] = net_start + cols[k];
783 val[k] = yvals[k];
784 }
785 PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
786 PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
787
788 PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
789 row[0] = net_start + 2 * i + 1;
790 for (k = 0; k < ncols; k++) {
791 col[k] = net_start + cols[k];
792 val[k] = yvals[k];
793 }
794 PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
795 PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
796 }
797
798 PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
799 PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));
800
801 PetscCall(VecGetArray(user->V0, &v0));
802 for (i = 0; i < nload; i++) {
803 Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
804 Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
805 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
806 Vm2 = Vm * Vm;
807 Vm4 = Vm2 * Vm2;
808 Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
809 PD = QD = 0.0;
810 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
811 for (k = 0; k < ld_nsegsp[i]; k++) {
812 PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
813 dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vr * PetscPowScalar(Vm, ld_betap[k] - 2);
814 dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vi * PetscPowScalar(Vm, ld_betap[k] - 2);
815 }
816 for (k = 0; k < ld_nsegsq[i]; k++) {
817 QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
818 dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vr * PetscPowScalar(Vm, ld_betaq[k] - 2);
819 dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vi * PetscPowScalar(Vm, ld_betaq[k] - 2);
820 }
821
822 /* IDr = (PD*Vr + QD*Vi)/Vm2; */
823 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
824
825 dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
826 dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
827
828 dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
829 dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
830
831 /* fnet[2*lbus[i]] += IDi; */
832 row[0] = net_start + 2 * lbus[i];
833 col[0] = net_start + 2 * lbus[i];
834 col[1] = net_start + 2 * lbus[i] + 1;
835 val[0] = dIDi_dVr;
836 val[1] = dIDi_dVi;
837 PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
838 /* fnet[2*lbus[i]+1] += IDr; */
839 row[0] = net_start + 2 * lbus[i] + 1;
840 col[0] = net_start + 2 * lbus[i];
841 col[1] = net_start + 2 * lbus[i] + 1;
842 val[0] = dIDr_dVr;
843 val[1] = dIDr_dVi;
844 PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
845 }
846 PetscCall(VecRestoreArray(user->V0, &v0));
847
848 PetscCall(VecRestoreArray(Xgen, &xgen));
849 PetscCall(VecRestoreArray(Xnet, &xnet));
850
851 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
852
853 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
854 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
855 PetscFunctionReturn(PETSC_SUCCESS);
856 }
857
858 /*
859 J = [I, 0
860 dg_dx, dg_dy]
861 */
AlgJacobian(SNES snes,Vec X,Mat A,Mat B,PetscCtx ctx)862 PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
863 {
864 Userctx *user = (Userctx *)ctx;
865
866 PetscFunctionBegin;
867 PetscCall(ResidualJacobian(snes, X, A, B, ctx));
868 PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
869 PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
870 PetscFunctionReturn(PETSC_SUCCESS);
871 }
872
873 /*
874 J = [a*I-df_dx, -df_dy
875 dg_dx, dg_dy]
876 */
877
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx * user)878 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
879 {
880 SNES snes;
881 PetscScalar atmp = (PetscScalar)a;
882 PetscInt i, row;
883
884 PetscFunctionBegin;
885 user->t = t;
886
887 PetscCall(TSGetSNES(ts, &snes));
888 PetscCall(ResidualJacobian(snes, X, A, B, user));
889 for (i = 0; i < ngen; i++) {
890 row = 9 * i;
891 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
892 row = 9 * i + 1;
893 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
894 row = 9 * i + 2;
895 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
896 row = 9 * i + 3;
897 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
898 row = 9 * i + 6;
899 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
900 row = 9 * i + 7;
901 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
902 row = 9 * i + 8;
903 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
904 }
905 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
906 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
907 PetscFunctionReturn(PETSC_SUCCESS);
908 }
909
910 /* Matrix JacobianP is constant so that it only needs to be evaluated once */
RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,PetscCtx ctx0)911 static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, PetscCtx ctx0)
912 {
913 PetscScalar a;
914 PetscInt row, col;
915 Userctx *ctx = (Userctx *)ctx0;
916
917 PetscFunctionBeginUser;
918 if (ctx->jacp_flg) {
919 PetscCall(MatZeroEntries(A));
920
921 for (col = 0; col < 3; col++) {
922 a = 1.0 / M[col];
923 row = 9 * col + 3;
924 PetscCall(MatSetValues(A, 1, &row, 1, &col, &a, INSERT_VALUES));
925 }
926
927 ctx->jacp_flg = PETSC_FALSE;
928
929 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
930 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
931 }
932 PetscFunctionReturn(PETSC_SUCCESS);
933 }
934
CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,Userctx * user)935 static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
936 {
937 const PetscScalar *u;
938 PetscInt idx;
939 Vec Xgen, Xnet;
940 PetscScalar *r, *xgen;
941 PetscInt i;
942
943 PetscFunctionBegin;
944 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
945 PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
946
947 PetscCall(VecGetArray(Xgen, &xgen));
948
949 PetscCall(VecGetArrayRead(U, &u));
950 PetscCall(VecGetArray(R, &r));
951 r[0] = 0.;
952 idx = 0;
953 for (i = 0; i < ngen; i++) {
954 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);
955 idx += 9;
956 }
957 PetscCall(VecRestoreArrayRead(U, &u));
958 PetscCall(VecRestoreArray(R, &r));
959 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
960 PetscFunctionReturn(PETSC_SUCCESS);
961 }
962
DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,Userctx * user)963 static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, Userctx *user)
964 {
965 Vec Xgen, Xnet, Dgen, Dnet;
966 PetscScalar *xgen, *dgen;
967 PetscInt i;
968 PetscInt idx;
969 Vec drdu_col;
970 PetscScalar *xarr;
971
972 PetscFunctionBegin;
973 PetscCall(VecDuplicate(U, &drdu_col));
974 PetscCall(MatDenseGetColumn(DRDU, 0, &xarr));
975 PetscCall(VecPlaceArray(drdu_col, xarr));
976 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
977 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Dgen, &Dnet));
978 PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
979 PetscCall(DMCompositeScatter(user->dmpgrid, drdu_col, Dgen, Dnet));
980
981 PetscCall(VecGetArray(Xgen, &xgen));
982 PetscCall(VecGetArray(Dgen, &dgen));
983
984 idx = 0;
985 for (i = 0; i < ngen; i++) {
986 dgen[idx + 3] = 0.;
987 if (xgen[idx + 3] / (2. * PETSC_PI) > user->freq_u) dgen[idx + 3] = user->pow * PetscPowScalarInt(xgen[idx + 3] / (2. * PETSC_PI) - user->freq_u, user->pow - 1) / (2. * PETSC_PI);
988 if (xgen[idx + 3] / (2. * PETSC_PI) < user->freq_l) dgen[idx + 3] = user->pow * PetscPowScalarInt(user->freq_l - xgen[idx + 3] / (2. * PETSC_PI), user->pow - 1) / (-2. * PETSC_PI);
989 idx += 9;
990 }
991
992 PetscCall(VecRestoreArray(Dgen, &dgen));
993 PetscCall(VecRestoreArray(Xgen, &xgen));
994 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, drdu_col, Dgen, Dnet));
995 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Dgen, &Dnet));
996 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
997 PetscCall(VecResetArray(drdu_col));
998 PetscCall(MatDenseRestoreColumn(DRDU, &xarr));
999 PetscCall(VecDestroy(&drdu_col));
1000 PetscFunctionReturn(PETSC_SUCCESS);
1001 }
1002
DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat drdp,Userctx * user)1003 static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat drdp, Userctx *user)
1004 {
1005 PetscFunctionBegin;
1006 PetscFunctionReturn(PETSC_SUCCESS);
1007 }
1008
ComputeSensiP(Vec lambda,Vec mu,Vec * DICDP,Userctx * user)1009 PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, Vec *DICDP, Userctx *user)
1010 {
1011 PetscScalar *x, *y, sensip;
1012 PetscInt i;
1013
1014 PetscFunctionBegin;
1015 PetscCall(VecGetArray(lambda, &x));
1016 PetscCall(VecGetArray(mu, &y));
1017
1018 for (i = 0; i < 3; i++) {
1019 PetscCall(VecDot(lambda, DICDP[i], &sensip));
1020 sensip = sensip + y[i];
1021 /* PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt %" PetscInt_FMT " th parameter: %g \n",i,(double)sensip)); */
1022 y[i] = sensip;
1023 }
1024 PetscCall(VecRestoreArray(mu, &y));
1025 PetscFunctionReturn(PETSC_SUCCESS);
1026 }
1027
main(int argc,char ** argv)1028 int main(int argc, char **argv)
1029 {
1030 Userctx user;
1031 Vec p;
1032 PetscScalar *x_ptr;
1033 PetscMPIInt size;
1034 PetscInt i;
1035 PetscViewer Xview, Ybusview;
1036 PetscInt *idx2;
1037 Tao tao;
1038 KSP ksp;
1039 PC pc;
1040 Vec lowerb, upperb;
1041
1042 PetscFunctionBeginUser;
1043 PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
1044 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1045 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
1046
1047 user.jacp_flg = PETSC_TRUE;
1048 user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
1049 user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
1050 user.neqs_pgrid = user.neqs_gen + user.neqs_net;
1051
1052 /* Create indices for differential and algebraic equations */
1053 PetscCall(PetscMalloc1(7 * ngen, &idx2));
1054 for (i = 0; i < ngen; i++) {
1055 idx2[7 * i] = 9 * i;
1056 idx2[7 * i + 1] = 9 * i + 1;
1057 idx2[7 * i + 2] = 9 * i + 2;
1058 idx2[7 * i + 3] = 9 * i + 3;
1059 idx2[7 * i + 4] = 9 * i + 6;
1060 idx2[7 * i + 5] = 9 * i + 7;
1061 idx2[7 * i + 6] = 9 * i + 8;
1062 }
1063 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
1064 PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
1065 PetscCall(PetscFree(idx2));
1066
1067 /* Set run time options */
1068 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1069 {
1070 user.tfaulton = 1.0;
1071 user.tfaultoff = 1.2;
1072 user.Rfault = 0.0001;
1073 user.faultbus = 8;
1074 PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
1075 PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
1076 PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
1077 user.t0 = 0.0;
1078 user.tmax = 1.3;
1079 PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
1080 PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
1081 user.freq_u = 61.0;
1082 user.freq_l = 59.0;
1083 user.pow = 2;
1084 PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
1085 PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
1086 PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
1087 }
1088 PetscOptionsEnd();
1089
1090 /* Create DMs for generator and network subsystems */
1091 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
1092 PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
1093 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
1094 PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
1095 PetscCall(DMSetFromOptions(user.dmnet));
1096 PetscCall(DMSetUp(user.dmnet));
1097 /* Create a composite DM packer and add the two DMs */
1098 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
1099 PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
1100 PetscCall(DMSetFromOptions(user.dmgen));
1101 PetscCall(DMSetUp(user.dmgen));
1102 PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
1103 PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));
1104
1105 /* Read initial voltage vector and Ybus */
1106 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
1107 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));
1108
1109 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.V0));
1110 PetscCall(VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net));
1111 PetscCall(VecLoad(user.V0, Xview));
1112
1113 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Ybus));
1114 PetscCall(MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net));
1115 PetscCall(MatSetType(user.Ybus, MATBAIJ));
1116 /* PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
1117 PetscCall(MatLoad(user.Ybus, Ybusview));
1118
1119 PetscCall(PetscViewerDestroy(&Xview));
1120 PetscCall(PetscViewerDestroy(&Ybusview));
1121
1122 /* Allocate space for Jacobians */
1123 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.J));
1124 PetscCall(MatSetSizes(user.J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid));
1125 PetscCall(MatSetFromOptions(user.J));
1126 PetscCall(PreallocateJacobian(user.J, &user));
1127
1128 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
1129 PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 3));
1130 PetscCall(MatSetFromOptions(user.Jacp));
1131 PetscCall(MatSetUp(user.Jacp));
1132 PetscCall(MatZeroEntries(user.Jacp)); /* initialize to zeros */
1133
1134 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, &user.DRDP));
1135 PetscCall(MatSetUp(user.DRDP));
1136 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 1, NULL, &user.DRDU));
1137 PetscCall(MatSetUp(user.DRDU));
1138
1139 /* Create TAO solver and set desired solution method */
1140 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
1141 PetscCall(TaoSetType(tao, TAOBLMVM));
1142 /*
1143 Optimization starts
1144 */
1145 /* Set initial solution guess */
1146 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
1147 PetscCall(VecGetArray(p, &x_ptr));
1148 x_ptr[0] = PG[0];
1149 x_ptr[1] = PG[1];
1150 x_ptr[2] = PG[2];
1151 PetscCall(VecRestoreArray(p, &x_ptr));
1152
1153 PetscCall(TaoSetSolution(tao, p));
1154 /* Set routine for function and gradient evaluation */
1155 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user));
1156
1157 /* Set bounds for the optimization */
1158 PetscCall(VecDuplicate(p, &lowerb));
1159 PetscCall(VecDuplicate(p, &upperb));
1160 PetscCall(VecGetArray(lowerb, &x_ptr));
1161 x_ptr[0] = 0.5;
1162 x_ptr[1] = 0.5;
1163 x_ptr[2] = 0.5;
1164 PetscCall(VecRestoreArray(lowerb, &x_ptr));
1165 PetscCall(VecGetArray(upperb, &x_ptr));
1166 x_ptr[0] = 2.0;
1167 x_ptr[1] = 2.0;
1168 x_ptr[2] = 2.0;
1169 PetscCall(VecRestoreArray(upperb, &x_ptr));
1170 PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));
1171
1172 /* Check for any TAO command line options */
1173 PetscCall(TaoSetFromOptions(tao));
1174 PetscCall(TaoGetKSP(tao, &ksp));
1175 if (ksp) {
1176 PetscCall(KSPGetPC(ksp, &pc));
1177 PetscCall(PCSetType(pc, PCNONE));
1178 }
1179
1180 /* SOLVE THE APPLICATION */
1181 PetscCall(TaoSolve(tao));
1182
1183 PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
1184 /* Free TAO data structures */
1185 PetscCall(TaoDestroy(&tao));
1186
1187 PetscCall(DMDestroy(&user.dmgen));
1188 PetscCall(DMDestroy(&user.dmnet));
1189 PetscCall(DMDestroy(&user.dmpgrid));
1190 PetscCall(ISDestroy(&user.is_diff));
1191 PetscCall(ISDestroy(&user.is_alg));
1192
1193 PetscCall(MatDestroy(&user.J));
1194 PetscCall(MatDestroy(&user.Jacp));
1195 PetscCall(MatDestroy(&user.Ybus));
1196 /* PetscCall(MatDestroy(&user.Sol)); */
1197 PetscCall(VecDestroy(&user.V0));
1198 PetscCall(VecDestroy(&p));
1199 PetscCall(VecDestroy(&lowerb));
1200 PetscCall(VecDestroy(&upperb));
1201 PetscCall(MatDestroy(&user.DRDU));
1202 PetscCall(MatDestroy(&user.DRDP));
1203 PetscCall(PetscFinalize());
1204 return 0;
1205 }
1206
1207 /* ------------------------------------------------------------------ */
1208 /*
1209 FormFunction - Evaluates the function and corresponding gradient.
1210
1211 Input Parameters:
1212 tao - the Tao context
1213 X - the input vector
1214 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
1215
1216 Output Parameters:
1217 f - the newly evaluated function
1218 G - the newly evaluated gradient
1219 */
FormFunctionGradient(Tao tao,Vec P,PetscReal * f,Vec G,PetscCtx ctx0)1220 PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, PetscCtx ctx0)
1221 {
1222 TS ts, quadts;
1223 SNES snes_alg;
1224 Userctx *ctx = (Userctx *)ctx0;
1225 Vec X;
1226 PetscInt i;
1227 /* sensitivity context */
1228 PetscScalar *x_ptr;
1229 Vec lambda[1], q;
1230 Vec mu[1];
1231 PetscInt steps1, steps2, steps3;
1232 Vec DICDP[3];
1233 Vec F_alg;
1234 PetscInt row_loc, col_loc;
1235 PetscScalar val;
1236 Vec Xdot;
1237
1238 PetscFunctionBegin;
1239 PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
1240 PG[0] = x_ptr[0];
1241 PG[1] = x_ptr[1];
1242 PG[2] = x_ptr[2];
1243 PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
1244
1245 ctx->stepnum = 0;
1246
1247 PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));
1248
1249 /* Create matrix to save solutions at each time step */
1250 /* PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ctx->neqs_pgrid+1,1002,NULL,&ctx->Sol)); */
1251 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1252 Create timestepping solver context
1253 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1254 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1255 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1256 PetscCall(TSSetType(ts, TSCN));
1257 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, ctx));
1258 PetscCall(TSSetIJacobian(ts, ctx->J, ctx->J, (TSIJacobianFn *)IJacobian, ctx));
1259 PetscCall(TSSetApplicationContext(ts, ctx));
1260 /* Set RHS JacobianP */
1261 PetscCall(TSSetRHSJacobianP(ts, ctx->Jacp, RHSJacobianP, ctx));
1262
1263 PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
1264 PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, ctx));
1265 PetscCall(TSSetRHSJacobian(quadts, ctx->DRDU, ctx->DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, ctx));
1266 PetscCall(TSSetRHSJacobianP(quadts, ctx->DRDP, (TSRHSJacobianPFn *)DRDPJacobianTranspose, ctx));
1267
1268 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1269 Set initial conditions
1270 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1271 PetscCall(SetInitialGuess(X, ctx));
1272
1273 /* Approximate DICDP with finite difference, we want to zero out network variables */
1274 for (i = 0; i < 3; i++) PetscCall(VecDuplicate(X, &DICDP[i]));
1275 PetscCall(DICDPFiniteDifference(X, DICDP, ctx));
1276
1277 PetscCall(VecDuplicate(X, &F_alg));
1278 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1279 PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
1280 PetscCall(MatZeroEntries(ctx->J));
1281 PetscCall(SNESSetJacobian(snes_alg, ctx->J, ctx->J, AlgJacobian, ctx));
1282 PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1283 PetscCall(SNESSetFromOptions(snes_alg));
1284 ctx->alg_flg = PETSC_TRUE;
1285 /* Solve the algebraic equations */
1286 PetscCall(SNESSolve(snes_alg, NULL, X));
1287
1288 /* Just to set up the Jacobian structure */
1289 PetscCall(VecDuplicate(X, &Xdot));
1290 PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, ctx->J, ctx->J, ctx));
1291 PetscCall(VecDestroy(&Xdot));
1292
1293 ctx->stepnum++;
1294
1295 /*
1296 Save trajectory of solution so that TSAdjointSolve() may be used
1297 */
1298 PetscCall(TSSetSaveTrajectory(ts));
1299
1300 PetscCall(TSSetTimeStep(ts, 0.01));
1301 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1302 PetscCall(TSSetFromOptions(ts));
1303 /* PetscCall(TSSetPostStep(ts,SaveSolution)); */
1304
1305 /* Prefault period */
1306 ctx->alg_flg = PETSC_FALSE;
1307 PetscCall(TSSetTime(ts, 0.0));
1308 PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
1309 PetscCall(TSSolve(ts, X));
1310 PetscCall(TSGetStepNumber(ts, &steps1));
1311
1312 /* Create the nonlinear solver for solving the algebraic system */
1313 /* Note that although the algebraic system needs to be solved only for
1314 Idq and V, we reuse the entire system including xgen. The xgen
1315 variables are held constant by setting their residuals to 0 and
1316 putting a 1 on the Jacobian diagonal for xgen rows
1317 */
1318 PetscCall(MatZeroEntries(ctx->J));
1319
1320 /* Apply disturbance - resistive fault at ctx->faultbus */
1321 /* This is done by adding shunt conductance to the diagonal location
1322 in the Ybus matrix */
1323 row_loc = 2 * ctx->faultbus;
1324 col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1325 val = 1 / ctx->Rfault;
1326 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1327 row_loc = 2 * ctx->faultbus + 1;
1328 col_loc = 2 * ctx->faultbus; /* Location for G */
1329 val = 1 / ctx->Rfault;
1330 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1331
1332 PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1333 PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1334
1335 ctx->alg_flg = PETSC_TRUE;
1336 /* Solve the algebraic equations */
1337 PetscCall(SNESSolve(snes_alg, NULL, X));
1338
1339 ctx->stepnum++;
1340
1341 /* Disturbance period */
1342 ctx->alg_flg = PETSC_FALSE;
1343 PetscCall(TSSetTime(ts, ctx->tfaulton));
1344 PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
1345 PetscCall(TSSolve(ts, X));
1346 PetscCall(TSGetStepNumber(ts, &steps2));
1347 steps2 -= steps1;
1348
1349 /* Remove the fault */
1350 row_loc = 2 * ctx->faultbus;
1351 col_loc = 2 * ctx->faultbus + 1;
1352 val = -1 / ctx->Rfault;
1353 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1354 row_loc = 2 * ctx->faultbus + 1;
1355 col_loc = 2 * ctx->faultbus;
1356 val = -1 / ctx->Rfault;
1357 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1358
1359 PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1360 PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1361
1362 PetscCall(MatZeroEntries(ctx->J));
1363
1364 ctx->alg_flg = PETSC_TRUE;
1365
1366 /* Solve the algebraic equations */
1367 PetscCall(SNESSolve(snes_alg, NULL, X));
1368
1369 ctx->stepnum++;
1370
1371 /* Post-disturbance period */
1372 ctx->alg_flg = PETSC_TRUE;
1373 PetscCall(TSSetTime(ts, ctx->tfaultoff));
1374 PetscCall(TSSetMaxTime(ts, ctx->tmax));
1375 PetscCall(TSSolve(ts, X));
1376 PetscCall(TSGetStepNumber(ts, &steps3));
1377 steps3 -= steps2;
1378 steps3 -= steps1;
1379
1380 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1381 Adjoint model starts here
1382 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1383 PetscCall(TSSetPostStep(ts, NULL));
1384 PetscCall(MatCreateVecs(ctx->J, &lambda[0], NULL));
1385 /* Set initial conditions for the adjoint integration */
1386 PetscCall(VecZeroEntries(lambda[0]));
1387
1388 PetscCall(MatCreateVecs(ctx->Jacp, &mu[0], NULL));
1389 PetscCall(VecZeroEntries(mu[0]));
1390 PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
1391
1392 PetscCall(TSAdjointSetSteps(ts, steps3));
1393 PetscCall(TSAdjointSolve(ts));
1394
1395 PetscCall(MatZeroEntries(ctx->J));
1396 /* Applying disturbance - resistive fault at ctx->faultbus */
1397 /* This is done by deducting shunt conductance to the diagonal location
1398 in the Ybus matrix */
1399 row_loc = 2 * ctx->faultbus;
1400 col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1401 val = 1. / ctx->Rfault;
1402 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1403 row_loc = 2 * ctx->faultbus + 1;
1404 col_loc = 2 * ctx->faultbus; /* Location for G */
1405 val = 1. / ctx->Rfault;
1406 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1407
1408 PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1409 PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1410
1411 /* Set number of steps for the adjoint integration */
1412 PetscCall(TSAdjointSetSteps(ts, steps2));
1413 PetscCall(TSAdjointSolve(ts));
1414
1415 PetscCall(MatZeroEntries(ctx->J));
1416 /* remove the fault */
1417 row_loc = 2 * ctx->faultbus;
1418 col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1419 val = -1. / ctx->Rfault;
1420 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1421 row_loc = 2 * ctx->faultbus + 1;
1422 col_loc = 2 * ctx->faultbus; /* Location for G */
1423 val = -1. / ctx->Rfault;
1424 PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1425
1426 PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1427 PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1428
1429 /* Set number of steps for the adjoint integration */
1430 PetscCall(TSAdjointSetSteps(ts, steps1));
1431 PetscCall(TSAdjointSolve(ts));
1432
1433 PetscCall(ComputeSensiP(lambda[0], mu[0], DICDP, ctx));
1434 PetscCall(VecCopy(mu[0], G));
1435
1436 PetscCall(TSGetQuadratureTS(ts, NULL, &quadts));
1437 PetscCall(TSGetSolution(quadts, &q));
1438 PetscCall(VecGetArray(q, &x_ptr));
1439 *f = x_ptr[0];
1440 x_ptr[0] = 0;
1441 PetscCall(VecRestoreArray(q, &x_ptr));
1442
1443 PetscCall(VecDestroy(&lambda[0]));
1444 PetscCall(VecDestroy(&mu[0]));
1445
1446 PetscCall(SNESDestroy(&snes_alg));
1447 PetscCall(VecDestroy(&F_alg));
1448 PetscCall(VecDestroy(&X));
1449 PetscCall(TSDestroy(&ts));
1450 for (i = 0; i < 3; i++) PetscCall(VecDestroy(&DICDP[i]));
1451 PetscFunctionReturn(PETSC_SUCCESS);
1452 }
1453
1454 /*TEST
1455
1456 build:
1457 requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1458
1459 test:
1460 args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1461 localrunfiles: petscoptions X.bin Ybus.bin
1462
1463 TEST*/
1464