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