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>
20: #include <petsctime.h>
22: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
24: #define freq 60
25: #define w_s (2 * PETSC_PI * freq)
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 */
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};*/
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 */
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) */
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
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};
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;
110: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
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: }
119: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
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: }
128: /* Saves the solution at each time to a matrix */
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;
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: }
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;
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];
172: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
174: /* Network subsystem initialization */
175: PetscCall(VecCopy(user->V0, Xnet));
177: /* Generator subsystem initialization */
178: PetscCall(VecGetArray(Xgen, &xgen));
179: PetscCall(VecGetArray(Xnet, &xnet));
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;
189: delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
191: theta = PETSC_PI / 2.0 - delta;
193: Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
194: Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
196: Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
197: Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
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 */
202: TM[i] = PG[i];
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;
210: idx = idx + 4;
212: xgen[idx] = Id;
213: xgen[idx + 1] = Iq;
215: idx = idx + 2;
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];
223: xgen[idx] = Efd;
224: xgen[idx + 1] = RF;
225: xgen[idx + 2] = VR;
227: Vref[i] = Vm + (VR / KA[i]);
229: idx = idx + 3;
230: }
232: PetscCall(VecRestoreArray(Xgen, &xgen));
233: PetscCall(VecRestoreArray(Xnet, &xnet));
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: }
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;
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];
260: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
262: /* Network subsystem initialization */
263: PetscCall(VecCopy(user->V0, Xnet));
265: /* Generator subsystem initialization */
266: PetscCall(VecGetArray(Xgen, &xgen));
267: PetscCall(VecGetArray(Xnet, &xnet));
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;
277: delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
279: theta = PETSC_PI / 2.0 - delta;
281: Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
282: Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
284: Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
285: Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
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 */
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;
296: idx = idx + 4;
298: xgen[idx] = Id;
299: xgen[idx + 1] = Iq;
301: idx = idx + 2;
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];
309: xgen[idx] = Efd;
310: xgen[idx + 1] = RF;
311: xgen[idx + 2] = VR;
313: idx = idx + 3;
314: }
316: PetscCall(VecRestoreArray(Xgen, &xgen));
317: PetscCall(VecRestoreArray(Xnet, &xnet));
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: }
325: PetscErrorCode DICDPFiniteDifference(Vec X, Vec *DICDP, Userctx *user)
326: {
327: Vec Y;
328: PetscScalar PGv[3], eps;
329: PetscInt i;
331: PetscFunctionBegin;
332: eps = 1.e-7;
333: PetscCall(VecDuplicate(X, &Y));
335: for (i = 0; i < ngen; i++) {
336: for (PetscInt 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));
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: }
349: /* Computes F = [-f(x,y);g(x,y)] */
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;
364: PetscFunctionBegin;
365: PetscCall(VecZeroEntries(F));
366: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
367: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
368: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
369: PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));
371: /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
372: The generator current injection, IG, and load current injection, ID are added later
373: */
374: /* Note that the values in Ybus are stored assuming the imaginary current balance
375: equation is ordered first followed by real current balance equation for each bus.
376: Thus imaginary current contribution goes in location 2*i, and
377: real current contribution in 2*i+1
378: */
379: PetscCall(MatMult(user->Ybus, Xnet, Fnet));
381: PetscCall(VecGetArray(Xgen, &xgen));
382: PetscCall(VecGetArray(Xnet, &xnet));
383: PetscCall(VecGetArray(Fgen, &fgen));
384: PetscCall(VecGetArray(Fnet, &fnet));
386: /* Generator subsystem */
387: for (i = 0; i < ngen; i++) {
388: Eqp = xgen[idx];
389: Edp = xgen[idx + 1];
390: delta = xgen[idx + 2];
391: w = xgen[idx + 3];
392: Id = xgen[idx + 4];
393: Iq = xgen[idx + 5];
394: Efd = xgen[idx + 6];
395: RF = xgen[idx + 7];
396: VR = xgen[idx + 8];
398: /* Generator differential equations */
399: fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
400: fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
401: fgen[idx + 2] = -w + w_s;
402: fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
404: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
405: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
407: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
408: /* Algebraic equations for stator currents */
409: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
411: Zdq_inv[0] = Rs[i] / det;
412: Zdq_inv[1] = Xqp[i] / det;
413: Zdq_inv[2] = -Xdp[i] / det;
414: Zdq_inv[3] = Rs[i] / det;
416: fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
417: fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
419: /* Add generator current injection to network */
420: PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
422: fnet[2 * gbus[i]] -= IGi;
423: fnet[2 * gbus[i] + 1] -= IGr;
425: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
427: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
429: /* Exciter differential equations */
430: fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
431: fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
432: fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
434: idx = idx + 9;
435: }
437: PetscCall(VecGetArray(user->V0, &v0));
438: for (i = 0; i < nload; i++) {
439: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
440: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
441: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
442: Vm2 = Vm * Vm;
443: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
444: PD = QD = 0.0;
445: for (PetscInt k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
446: for (PetscInt k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
448: /* Load currents */
449: IDr = (PD * Vr + QD * Vi) / Vm2;
450: IDi = (-QD * Vr + PD * Vi) / Vm2;
452: fnet[2 * lbus[i]] += IDi;
453: fnet[2 * lbus[i] + 1] += IDr;
454: }
455: PetscCall(VecRestoreArray(user->V0, &v0));
457: PetscCall(VecRestoreArray(Xgen, &xgen));
458: PetscCall(VecRestoreArray(Xnet, &xnet));
459: PetscCall(VecRestoreArray(Fgen, &fgen));
460: PetscCall(VecRestoreArray(Fnet, &fnet));
462: PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
463: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
464: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /* \dot{x} - f(x,y)
469: g(x,y) = 0
470: */
471: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
472: {
473: SNES snes;
474: PetscScalar *f;
475: const PetscScalar *xdot;
476: PetscInt i;
478: PetscFunctionBegin;
479: user->t = t;
481: PetscCall(TSGetSNES(ts, &snes));
482: PetscCall(ResidualFunction(snes, X, F, user));
483: PetscCall(VecGetArray(F, &f));
484: PetscCall(VecGetArrayRead(Xdot, &xdot));
485: for (i = 0; i < ngen; i++) {
486: f[9 * i] += xdot[9 * i];
487: f[9 * i + 1] += xdot[9 * i + 1];
488: f[9 * i + 2] += xdot[9 * i + 2];
489: f[9 * i + 3] += xdot[9 * i + 3];
490: f[9 * i + 6] += xdot[9 * i + 6];
491: f[9 * i + 7] += xdot[9 * i + 7];
492: f[9 * i + 8] += xdot[9 * i + 8];
493: }
494: PetscCall(VecRestoreArray(F, &f));
495: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /* This function is used for solving the algebraic system only during fault on and
500: off times. It computes the entire F and then zeros out the part corresponding to
501: differential equations
502: F = [0;g(y)];
503: */
504: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, PetscCtx ctx)
505: {
506: Userctx *user = (Userctx *)ctx;
507: PetscScalar *f;
508: PetscInt i;
510: PetscFunctionBegin;
511: PetscCall(ResidualFunction(snes, X, F, user));
512: PetscCall(VecGetArray(F, &f));
513: for (i = 0; i < ngen; i++) {
514: f[9 * i] = 0;
515: f[9 * i + 1] = 0;
516: f[9 * i + 2] = 0;
517: f[9 * i + 3] = 0;
518: f[9 * i + 6] = 0;
519: f[9 * i + 7] = 0;
520: f[9 * i + 8] = 0;
521: }
522: PetscCall(VecRestoreArray(F, &f));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
527: {
528: PetscInt *d_nnz;
529: PetscInt i, idx = 0, start = 0;
530: PetscInt ncols;
532: PetscFunctionBegin;
533: PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
534: for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
535: /* Generator subsystem */
536: for (i = 0; i < ngen; i++) {
537: d_nnz[idx] += 3;
538: d_nnz[idx + 1] += 2;
539: d_nnz[idx + 2] += 2;
540: d_nnz[idx + 3] += 5;
541: d_nnz[idx + 4] += 6;
542: d_nnz[idx + 5] += 6;
544: d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
545: d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
547: d_nnz[idx + 6] += 2;
548: d_nnz[idx + 7] += 2;
549: d_nnz[idx + 8] += 5;
551: idx = idx + 9;
552: }
554: start = user->neqs_gen;
555: for (i = 0; i < nbus; i++) {
556: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
557: d_nnz[start + 2 * i] += ncols;
558: d_nnz[start + 2 * i + 1] += ncols;
559: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
560: }
562: PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
563: PetscCall(PetscFree(d_nnz));
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /*
568: J = [-df_dx, -df_dy
569: dg_dx, dg_dy]
570: */
571: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, PetscCtx ctx)
572: {
573: Userctx *user = (Userctx *)ctx;
574: Vec Xgen, Xnet;
575: PetscScalar *xgen, *xnet;
576: PetscInt i, idx = 0;
577: PetscScalar Vr, Vi, Vm, Vm2;
578: PetscScalar Eqp, Edp, delta; /* Generator variables */
579: PetscScalar Efd; /* Exciter variables */
580: PetscScalar Id, Iq; /* Generator dq axis currents */
581: PetscScalar Vd, Vq;
582: PetscScalar val[10];
583: PetscInt row[2], col[10];
584: PetscInt net_start = user->neqs_gen;
585: PetscInt ncols;
586: const PetscInt *cols;
587: const PetscScalar *yvals;
588: PetscScalar Zdq_inv[4], det;
589: PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
590: PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
591: PetscScalar dSE_dEfd;
592: PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
593: PetscScalar PD, QD, Vm0, *v0, Vm4;
594: PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
595: PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
597: PetscFunctionBegin;
598: PetscCall(MatZeroEntries(B));
599: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
600: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
602: PetscCall(VecGetArray(Xgen, &xgen));
603: PetscCall(VecGetArray(Xnet, &xnet));
605: /* Generator subsystem */
606: for (i = 0; i < ngen; i++) {
607: Eqp = xgen[idx];
608: Edp = xgen[idx + 1];
609: delta = xgen[idx + 2];
610: Id = xgen[idx + 4];
611: Iq = xgen[idx + 5];
612: Efd = xgen[idx + 6];
614: /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
615: row[0] = idx;
616: col[0] = idx;
617: col[1] = idx + 4;
618: col[2] = idx + 6;
619: val[0] = 1 / Td0p[i];
620: val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
621: val[2] = -1 / Td0p[i];
623: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
625: /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
626: row[0] = idx + 1;
627: col[0] = idx + 1;
628: col[1] = idx + 5;
629: val[0] = 1 / Tq0p[i];
630: val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
631: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
633: /* fgen[idx+2] = - w + w_s; */
634: row[0] = idx + 2;
635: col[0] = idx + 2;
636: col[1] = idx + 3;
637: val[0] = 0;
638: val[1] = -1;
639: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
641: /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
642: row[0] = idx + 3;
643: col[0] = idx;
644: col[1] = idx + 1;
645: col[2] = idx + 3;
646: col[3] = idx + 4;
647: col[4] = idx + 5;
648: val[0] = Iq / M[i];
649: val[1] = Id / M[i];
650: val[2] = D[i] / M[i];
651: val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
652: val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
653: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
655: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
656: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
657: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
659: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
661: Zdq_inv[0] = Rs[i] / det;
662: Zdq_inv[1] = Xqp[i] / det;
663: Zdq_inv[2] = -Xdp[i] / det;
664: Zdq_inv[3] = Rs[i] / det;
666: dVd_dVr = PetscSinScalar(delta);
667: dVd_dVi = -PetscCosScalar(delta);
668: dVq_dVr = PetscCosScalar(delta);
669: dVq_dVi = PetscSinScalar(delta);
670: dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
671: dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
673: /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
674: row[0] = idx + 4;
675: col[0] = idx;
676: col[1] = idx + 1;
677: col[2] = idx + 2;
678: val[0] = -Zdq_inv[1];
679: val[1] = -Zdq_inv[0];
680: val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
681: col[3] = idx + 4;
682: col[4] = net_start + 2 * gbus[i];
683: col[5] = net_start + 2 * gbus[i] + 1;
684: val[3] = 1;
685: val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
686: val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
687: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
689: /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
690: row[0] = idx + 5;
691: col[0] = idx;
692: col[1] = idx + 1;
693: col[2] = idx + 2;
694: val[0] = -Zdq_inv[3];
695: val[1] = -Zdq_inv[2];
696: val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
697: col[3] = idx + 5;
698: col[4] = net_start + 2 * gbus[i];
699: col[5] = net_start + 2 * gbus[i] + 1;
700: val[3] = 1;
701: val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
702: val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
703: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
705: dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
706: dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
707: dIGr_dId = PetscSinScalar(delta);
708: dIGr_dIq = PetscCosScalar(delta);
709: dIGi_dId = -PetscCosScalar(delta);
710: dIGi_dIq = PetscSinScalar(delta);
712: /* fnet[2*gbus[i]] -= IGi; */
713: row[0] = net_start + 2 * gbus[i];
714: col[0] = idx + 2;
715: col[1] = idx + 4;
716: col[2] = idx + 5;
717: val[0] = -dIGi_ddelta;
718: val[1] = -dIGi_dId;
719: val[2] = -dIGi_dIq;
720: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
722: /* fnet[2*gbus[i]+1] -= IGr; */
723: row[0] = net_start + 2 * gbus[i] + 1;
724: col[0] = idx + 2;
725: col[1] = idx + 4;
726: col[2] = idx + 5;
727: val[0] = -dIGr_ddelta;
728: val[1] = -dIGr_dId;
729: val[2] = -dIGr_dIq;
730: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
732: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
734: /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
735: /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
736: dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
738: row[0] = idx + 6;
739: col[0] = idx + 6;
740: col[1] = idx + 8;
741: val[0] = (KE[i] + dSE_dEfd) / TE[i];
742: val[1] = -1 / TE[i];
743: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
745: /* Exciter differential equations */
747: /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
748: row[0] = idx + 7;
749: col[0] = idx + 6;
750: col[1] = idx + 7;
751: val[0] = (-KF[i] / TF[i]) / TF[i];
752: val[1] = 1 / TF[i];
753: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
755: /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
756: /* Vm = (Vd^2 + Vq^2)^0.5; */
757: dVm_dVd = Vd / Vm;
758: dVm_dVq = Vq / Vm;
759: dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
760: dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
761: row[0] = idx + 8;
762: col[0] = idx + 6;
763: col[1] = idx + 7;
764: col[2] = idx + 8;
765: val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
766: val[1] = -KA[i] / TA[i];
767: val[2] = 1 / TA[i];
768: col[3] = net_start + 2 * gbus[i];
769: col[4] = net_start + 2 * gbus[i] + 1;
770: val[3] = KA[i] * dVm_dVr / TA[i];
771: val[4] = KA[i] * dVm_dVi / TA[i];
772: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
773: idx = idx + 9;
774: }
776: for (i = 0; i < nbus; i++) {
777: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
778: row[0] = net_start + 2 * i;
779: for (PetscInt k = 0; k < ncols; k++) {
780: col[k] = net_start + cols[k];
781: val[k] = yvals[k];
782: }
783: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
784: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
786: PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
787: row[0] = net_start + 2 * i + 1;
788: for (PetscInt k = 0; k < ncols; k++) {
789: col[k] = net_start + cols[k];
790: val[k] = yvals[k];
791: }
792: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
793: PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
794: }
796: PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
797: PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));
799: PetscCall(VecGetArray(user->V0, &v0));
800: for (i = 0; i < nload; i++) {
801: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
802: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
803: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
804: Vm2 = Vm * Vm;
805: Vm4 = Vm2 * Vm2;
806: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
807: PD = QD = 0.0;
808: dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
809: for (PetscInt k = 0; k < ld_nsegsp[i]; k++) {
810: PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
811: dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vr * PetscPowScalar(Vm, ld_betap[k] - 2);
812: dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vi * PetscPowScalar(Vm, ld_betap[k] - 2);
813: }
814: for (PetscInt k = 0; k < ld_nsegsq[i]; k++) {
815: QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
816: dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vr * PetscPowScalar(Vm, ld_betaq[k] - 2);
817: dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vi * PetscPowScalar(Vm, ld_betaq[k] - 2);
818: }
820: /* IDr = (PD*Vr + QD*Vi)/Vm2; */
821: /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
823: dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
824: dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
826: dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
827: dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
829: /* fnet[2*lbus[i]] += IDi; */
830: row[0] = net_start + 2 * lbus[i];
831: col[0] = net_start + 2 * lbus[i];
832: col[1] = net_start + 2 * lbus[i] + 1;
833: val[0] = dIDi_dVr;
834: val[1] = dIDi_dVi;
835: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
836: /* fnet[2*lbus[i]+1] += IDr; */
837: row[0] = net_start + 2 * lbus[i] + 1;
838: col[0] = net_start + 2 * lbus[i];
839: col[1] = net_start + 2 * lbus[i] + 1;
840: val[0] = dIDr_dVr;
841: val[1] = dIDr_dVi;
842: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
843: }
844: PetscCall(VecRestoreArray(user->V0, &v0));
846: PetscCall(VecRestoreArray(Xgen, &xgen));
847: PetscCall(VecRestoreArray(Xnet, &xnet));
849: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
851: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
852: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
853: PetscFunctionReturn(PETSC_SUCCESS);
854: }
856: /*
857: J = [I, 0
858: dg_dx, dg_dy]
859: */
860: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
861: {
862: Userctx *user = (Userctx *)ctx;
864: PetscFunctionBegin;
865: PetscCall(ResidualJacobian(snes, X, A, B, ctx));
866: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
867: PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: /*
872: J = [a*I-df_dx, -df_dy
873: dg_dx, dg_dy]
874: */
876: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
877: {
878: SNES snes;
879: PetscScalar atmp = (PetscScalar)a;
880: PetscInt i, row;
882: PetscFunctionBegin;
883: user->t = t;
885: PetscCall(TSGetSNES(ts, &snes));
886: PetscCall(ResidualJacobian(snes, X, A, B, user));
887: for (i = 0; i < ngen; i++) {
888: row = 9 * i;
889: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
890: row = 9 * i + 1;
891: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
892: row = 9 * i + 2;
893: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
894: row = 9 * i + 3;
895: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
896: row = 9 * i + 6;
897: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
898: row = 9 * i + 7;
899: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
900: row = 9 * i + 8;
901: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
902: }
903: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
904: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: /* Matrix JacobianP is constant so that it only needs to be evaluated once */
909: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, PetscCtx ctx0)
910: {
911: PetscScalar a;
912: PetscInt row;
913: Userctx *ctx = (Userctx *)ctx0;
915: PetscFunctionBeginUser;
916: if (ctx->jacp_flg) {
917: PetscCall(MatZeroEntries(A));
919: for (PetscInt col = 0; col < 3; col++) {
920: a = 1.0 / M[col];
921: row = 9 * col + 3;
922: PetscCall(MatSetValues(A, 1, &row, 1, &col, &a, INSERT_VALUES));
923: }
925: ctx->jacp_flg = PETSC_FALSE;
927: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
928: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
929: }
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
934: {
935: const PetscScalar *u;
936: PetscInt idx;
937: Vec Xgen, Xnet;
938: PetscScalar *r, *xgen;
939: PetscInt i;
941: PetscFunctionBegin;
942: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
943: PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
945: PetscCall(VecGetArray(Xgen, &xgen));
947: PetscCall(VecGetArrayRead(U, &u));
948: PetscCall(VecGetArray(R, &r));
949: r[0] = 0.;
950: idx = 0;
951: for (i = 0; i < ngen; i++) {
952: 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);
953: idx += 9;
954: }
955: PetscCall(VecRestoreArrayRead(U, &u));
956: PetscCall(VecRestoreArray(R, &r));
957: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, Userctx *user)
962: {
963: Vec Xgen, Xnet, Dgen, Dnet;
964: PetscScalar *xgen, *dgen;
965: PetscInt idx;
966: Vec drdu_col;
967: PetscScalar *xarr;
969: PetscFunctionBegin;
970: PetscCall(VecDuplicate(U, &drdu_col));
971: PetscCall(MatDenseGetColumn(DRDU, 0, &xarr));
972: PetscCall(VecPlaceArray(drdu_col, xarr));
973: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
974: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Dgen, &Dnet));
975: PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
976: PetscCall(DMCompositeScatter(user->dmpgrid, drdu_col, Dgen, Dnet));
978: PetscCall(VecGetArray(Xgen, &xgen));
979: PetscCall(VecGetArray(Dgen, &dgen));
981: idx = 0;
982: for (PetscInt i = 0; i < ngen; i++) {
983: dgen[idx + 3] = 0.;
984: 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);
985: 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);
986: idx += 9;
987: }
989: PetscCall(VecRestoreArray(Dgen, &dgen));
990: PetscCall(VecRestoreArray(Xgen, &xgen));
991: PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, drdu_col, Dgen, Dnet));
992: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Dgen, &Dnet));
993: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
994: PetscCall(VecResetArray(drdu_col));
995: PetscCall(MatDenseRestoreColumn(DRDU, &xarr));
996: PetscCall(VecDestroy(&drdu_col));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat drdp, Userctx *user)
1001: {
1002: PetscFunctionBegin;
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, Vec *DICDP, Userctx *user)
1007: {
1008: PetscScalar *x, *y, sensip;
1010: PetscFunctionBegin;
1011: PetscCall(VecGetArray(lambda, &x));
1012: PetscCall(VecGetArray(mu, &y));
1014: for (PetscInt i = 0; i < 3; i++) {
1015: PetscCall(VecDot(lambda, DICDP[i], &sensip));
1016: sensip = sensip + y[i];
1017: /* PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt %" PetscInt_FMT " th parameter: %g \n",i,(double)sensip)); */
1018: y[i] = sensip;
1019: }
1020: PetscCall(VecRestoreArray(mu, &y));
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: int main(int argc, char **argv)
1025: {
1026: Userctx user;
1027: Vec p;
1028: PetscScalar *x_ptr;
1029: PetscMPIInt size;
1030: PetscViewer Xview, Ybusview;
1031: PetscInt *idx2;
1032: Tao tao;
1033: KSP ksp;
1034: PC pc;
1035: Vec lowerb, upperb;
1037: PetscFunctionBeginUser;
1038: PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
1039: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1040: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
1042: user.jacp_flg = PETSC_TRUE;
1043: user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
1044: user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
1045: user.neqs_pgrid = user.neqs_gen + user.neqs_net;
1047: /* Create indices for differential and algebraic equations */
1048: PetscCall(PetscMalloc1(7 * ngen, &idx2));
1049: for (PetscInt i = 0; i < ngen; i++) {
1050: idx2[7 * i] = 9 * i;
1051: idx2[7 * i + 1] = 9 * i + 1;
1052: idx2[7 * i + 2] = 9 * i + 2;
1053: idx2[7 * i + 3] = 9 * i + 3;
1054: idx2[7 * i + 4] = 9 * i + 6;
1055: idx2[7 * i + 5] = 9 * i + 7;
1056: idx2[7 * i + 6] = 9 * i + 8;
1057: }
1058: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
1059: PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
1060: PetscCall(PetscFree(idx2));
1062: /* Set run time options */
1063: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1064: {
1065: user.tfaulton = 1.0;
1066: user.tfaultoff = 1.2;
1067: user.Rfault = 0.0001;
1068: user.faultbus = 8;
1069: PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
1070: PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
1071: PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
1072: user.t0 = 0.0;
1073: user.tmax = 1.3;
1074: PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
1075: PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
1076: user.freq_u = 61.0;
1077: user.freq_l = 59.0;
1078: user.pow = 2;
1079: PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
1080: PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
1081: PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
1082: }
1083: PetscOptionsEnd();
1085: /* Create DMs for generator and network subsystems */
1086: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
1087: PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
1088: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
1089: PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
1090: PetscCall(DMSetFromOptions(user.dmnet));
1091: PetscCall(DMSetUp(user.dmnet));
1092: /* Create a composite DM packer and add the two DMs */
1093: PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
1094: PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
1095: PetscCall(DMSetFromOptions(user.dmgen));
1096: PetscCall(DMSetUp(user.dmgen));
1097: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
1098: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));
1100: /* Read initial voltage vector and Ybus */
1101: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
1102: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));
1104: PetscCall(VecCreate(PETSC_COMM_WORLD, &user.V0));
1105: PetscCall(VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net));
1106: PetscCall(VecLoad(user.V0, Xview));
1108: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Ybus));
1109: PetscCall(MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net));
1110: PetscCall(MatSetType(user.Ybus, MATBAIJ));
1111: /* PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
1112: PetscCall(MatLoad(user.Ybus, Ybusview));
1114: PetscCall(PetscViewerDestroy(&Xview));
1115: PetscCall(PetscViewerDestroy(&Ybusview));
1117: /* Allocate space for Jacobians */
1118: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.J));
1119: PetscCall(MatSetSizes(user.J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid));
1120: PetscCall(MatSetFromOptions(user.J));
1121: PetscCall(PreallocateJacobian(user.J, &user));
1123: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
1124: PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 3));
1125: PetscCall(MatSetFromOptions(user.Jacp));
1126: PetscCall(MatSetUp(user.Jacp));
1127: PetscCall(MatZeroEntries(user.Jacp)); /* initialize to zeros */
1129: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, &user.DRDP));
1130: PetscCall(MatSetUp(user.DRDP));
1131: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 1, NULL, &user.DRDU));
1132: PetscCall(MatSetUp(user.DRDU));
1134: /* Create TAO solver and set desired solution method */
1135: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
1136: PetscCall(TaoSetType(tao, TAOBLMVM));
1137: /*
1138: Optimization starts
1139: */
1140: /* Set initial solution guess */
1141: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
1142: PetscCall(VecGetArray(p, &x_ptr));
1143: x_ptr[0] = PG[0];
1144: x_ptr[1] = PG[1];
1145: x_ptr[2] = PG[2];
1146: PetscCall(VecRestoreArray(p, &x_ptr));
1148: PetscCall(TaoSetSolution(tao, p));
1149: /* Set routine for function and gradient evaluation */
1150: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user));
1152: /* Set bounds for the optimization */
1153: PetscCall(VecDuplicate(p, &lowerb));
1154: PetscCall(VecDuplicate(p, &upperb));
1155: PetscCall(VecGetArray(lowerb, &x_ptr));
1156: x_ptr[0] = 0.5;
1157: x_ptr[1] = 0.5;
1158: x_ptr[2] = 0.5;
1159: PetscCall(VecRestoreArray(lowerb, &x_ptr));
1160: PetscCall(VecGetArray(upperb, &x_ptr));
1161: x_ptr[0] = 2.0;
1162: x_ptr[1] = 2.0;
1163: x_ptr[2] = 2.0;
1164: PetscCall(VecRestoreArray(upperb, &x_ptr));
1165: PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));
1167: /* Check for any TAO command line options */
1168: PetscCall(TaoSetFromOptions(tao));
1169: PetscCall(TaoGetKSP(tao, &ksp));
1170: if (ksp) {
1171: PetscCall(KSPGetPC(ksp, &pc));
1172: PetscCall(PCSetType(pc, PCNONE));
1173: }
1175: /* SOLVE THE APPLICATION */
1176: PetscCall(TaoSolve(tao));
1178: PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
1179: /* Free TAO data structures */
1180: PetscCall(TaoDestroy(&tao));
1182: PetscCall(DMDestroy(&user.dmgen));
1183: PetscCall(DMDestroy(&user.dmnet));
1184: PetscCall(DMDestroy(&user.dmpgrid));
1185: PetscCall(ISDestroy(&user.is_diff));
1186: PetscCall(ISDestroy(&user.is_alg));
1188: PetscCall(MatDestroy(&user.J));
1189: PetscCall(MatDestroy(&user.Jacp));
1190: PetscCall(MatDestroy(&user.Ybus));
1191: /* PetscCall(MatDestroy(&user.Sol)); */
1192: PetscCall(VecDestroy(&user.V0));
1193: PetscCall(VecDestroy(&p));
1194: PetscCall(VecDestroy(&lowerb));
1195: PetscCall(VecDestroy(&upperb));
1196: PetscCall(MatDestroy(&user.DRDU));
1197: PetscCall(MatDestroy(&user.DRDP));
1198: PetscCall(PetscFinalize());
1199: return 0;
1200: }
1202: /* ------------------------------------------------------------------ */
1203: /*
1204: FormFunction - Evaluates the function and corresponding gradient.
1206: Input Parameters:
1207: tao - the Tao context
1208: X - the input vector
1209: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
1211: Output Parameters:
1212: f - the newly evaluated function
1213: G - the newly evaluated gradient
1214: */
1215: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, PetscCtx ctx0)
1216: {
1217: TS ts, quadts;
1218: SNES snes_alg;
1219: Userctx *ctx = (Userctx *)ctx0;
1220: Vec X;
1221: /* sensitivity context */
1222: PetscScalar *x_ptr;
1223: Vec lambda[1], q;
1224: Vec mu[1];
1225: PetscInt steps1, steps2, steps3;
1226: Vec DICDP[3];
1227: Vec F_alg;
1228: PetscInt row_loc, col_loc;
1229: PetscScalar val;
1230: Vec Xdot;
1232: PetscFunctionBegin;
1233: PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
1234: PG[0] = x_ptr[0];
1235: PG[1] = x_ptr[1];
1236: PG[2] = x_ptr[2];
1237: PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
1239: ctx->stepnum = 0;
1241: PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));
1243: /* Create matrix to save solutions at each time step */
1244: /* PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ctx->neqs_pgrid+1,1002,NULL,&ctx->Sol)); */
1245: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1246: Create timestepping solver context
1247: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1248: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1249: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1250: PetscCall(TSSetType(ts, TSCN));
1251: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, ctx));
1252: PetscCall(TSSetIJacobian(ts, ctx->J, ctx->J, (TSIJacobianFn *)IJacobian, ctx));
1253: PetscCall(TSSetApplicationContext(ts, ctx));
1254: /* Set RHS JacobianP */
1255: PetscCall(TSSetRHSJacobianP(ts, ctx->Jacp, RHSJacobianP, ctx));
1257: PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
1258: PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, ctx));
1259: PetscCall(TSSetRHSJacobian(quadts, ctx->DRDU, ctx->DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, ctx));
1260: PetscCall(TSSetRHSJacobianP(quadts, ctx->DRDP, (TSRHSJacobianPFn *)DRDPJacobianTranspose, ctx));
1262: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1263: Set initial conditions
1264: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1265: PetscCall(SetInitialGuess(X, ctx));
1267: /* Approximate DICDP with finite difference, we want to zero out network variables */
1268: for (PetscInt i = 0; i < 3; i++) PetscCall(VecDuplicate(X, &DICDP[i]));
1269: PetscCall(DICDPFiniteDifference(X, DICDP, ctx));
1271: PetscCall(VecDuplicate(X, &F_alg));
1272: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1273: PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
1274: PetscCall(MatZeroEntries(ctx->J));
1275: PetscCall(SNESSetJacobian(snes_alg, ctx->J, ctx->J, AlgJacobian, ctx));
1276: PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1277: PetscCall(SNESSetFromOptions(snes_alg));
1278: ctx->alg_flg = PETSC_TRUE;
1279: /* Solve the algebraic equations */
1280: PetscCall(SNESSolve(snes_alg, NULL, X));
1282: /* Just to set up the Jacobian structure */
1283: PetscCall(VecDuplicate(X, &Xdot));
1284: PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, ctx->J, ctx->J, ctx));
1285: PetscCall(VecDestroy(&Xdot));
1287: ctx->stepnum++;
1289: /*
1290: Save trajectory of solution so that TSAdjointSolve() may be used
1291: */
1292: PetscCall(TSSetSaveTrajectory(ts));
1294: PetscCall(TSSetTimeStep(ts, 0.01));
1295: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1296: PetscCall(TSSetFromOptions(ts));
1297: /* PetscCall(TSSetPostStep(ts,SaveSolution)); */
1299: /* Prefault period */
1300: ctx->alg_flg = PETSC_FALSE;
1301: PetscCall(TSSetTime(ts, 0.0));
1302: PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
1303: PetscCall(TSSolve(ts, X));
1304: PetscCall(TSGetStepNumber(ts, &steps1));
1306: /* Create the nonlinear solver for solving the algebraic system */
1307: /* Note that although the algebraic system needs to be solved only for
1308: Idq and V, we reuse the entire system including xgen. The xgen
1309: variables are held constant by setting their residuals to 0 and
1310: putting a 1 on the Jacobian diagonal for xgen rows
1311: */
1312: PetscCall(MatZeroEntries(ctx->J));
1314: /* Apply disturbance - resistive fault at ctx->faultbus */
1315: /* This is done by adding shunt conductance to the diagonal location
1316: in the Ybus matrix */
1317: row_loc = 2 * ctx->faultbus;
1318: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1319: val = 1 / ctx->Rfault;
1320: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1321: row_loc = 2 * ctx->faultbus + 1;
1322: col_loc = 2 * ctx->faultbus; /* Location for G */
1323: val = 1 / ctx->Rfault;
1324: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1326: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1327: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1329: ctx->alg_flg = PETSC_TRUE;
1330: /* Solve the algebraic equations */
1331: PetscCall(SNESSolve(snes_alg, NULL, X));
1333: ctx->stepnum++;
1335: /* Disturbance period */
1336: ctx->alg_flg = PETSC_FALSE;
1337: PetscCall(TSSetTime(ts, ctx->tfaulton));
1338: PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
1339: PetscCall(TSSolve(ts, X));
1340: PetscCall(TSGetStepNumber(ts, &steps2));
1341: steps2 -= steps1;
1343: /* Remove the fault */
1344: row_loc = 2 * ctx->faultbus;
1345: col_loc = 2 * ctx->faultbus + 1;
1346: val = -1 / ctx->Rfault;
1347: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1348: row_loc = 2 * ctx->faultbus + 1;
1349: col_loc = 2 * ctx->faultbus;
1350: val = -1 / ctx->Rfault;
1351: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1353: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1354: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1356: PetscCall(MatZeroEntries(ctx->J));
1358: ctx->alg_flg = PETSC_TRUE;
1360: /* Solve the algebraic equations */
1361: PetscCall(SNESSolve(snes_alg, NULL, X));
1363: ctx->stepnum++;
1365: /* Post-disturbance period */
1366: ctx->alg_flg = PETSC_TRUE;
1367: PetscCall(TSSetTime(ts, ctx->tfaultoff));
1368: PetscCall(TSSetMaxTime(ts, ctx->tmax));
1369: PetscCall(TSSolve(ts, X));
1370: PetscCall(TSGetStepNumber(ts, &steps3));
1371: steps3 -= steps2;
1372: steps3 -= steps1;
1374: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1375: Adjoint model starts here
1376: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1377: PetscCall(TSSetPostStep(ts, NULL));
1378: PetscCall(MatCreateVecs(ctx->J, &lambda[0], NULL));
1379: /* Set initial conditions for the adjoint integration */
1381: PetscCall(MatCreateVecs(ctx->Jacp, &mu[0], NULL));
1382: PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
1384: PetscCall(TSAdjointSetSteps(ts, steps3));
1385: PetscCall(TSAdjointSolve(ts));
1387: PetscCall(MatZeroEntries(ctx->J));
1388: /* Applying disturbance - resistive fault at ctx->faultbus */
1389: /* This is done by deducting shunt conductance to the diagonal location
1390: in the Ybus matrix */
1391: row_loc = 2 * ctx->faultbus;
1392: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1393: val = 1. / ctx->Rfault;
1394: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1395: row_loc = 2 * ctx->faultbus + 1;
1396: col_loc = 2 * ctx->faultbus; /* Location for G */
1397: val = 1. / ctx->Rfault;
1398: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1400: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1401: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1403: /* Set number of steps for the adjoint integration */
1404: PetscCall(TSAdjointSetSteps(ts, steps2));
1405: PetscCall(TSAdjointSolve(ts));
1407: PetscCall(MatZeroEntries(ctx->J));
1408: /* remove the fault */
1409: row_loc = 2 * ctx->faultbus;
1410: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1411: val = -1. / ctx->Rfault;
1412: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1413: row_loc = 2 * ctx->faultbus + 1;
1414: col_loc = 2 * ctx->faultbus; /* Location for G */
1415: val = -1. / ctx->Rfault;
1416: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1418: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1419: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1421: /* Set number of steps for the adjoint integration */
1422: PetscCall(TSAdjointSetSteps(ts, steps1));
1423: PetscCall(TSAdjointSolve(ts));
1425: PetscCall(ComputeSensiP(lambda[0], mu[0], DICDP, ctx));
1426: PetscCall(VecCopy(mu[0], G));
1428: PetscCall(TSGetQuadratureTS(ts, NULL, &quadts));
1429: PetscCall(TSGetSolution(quadts, &q));
1430: PetscCall(VecGetArray(q, &x_ptr));
1431: *f = x_ptr[0];
1432: x_ptr[0] = 0;
1433: PetscCall(VecRestoreArray(q, &x_ptr));
1435: PetscCall(VecDestroy(&lambda[0]));
1436: PetscCall(VecDestroy(&mu[0]));
1438: PetscCall(SNESDestroy(&snes_alg));
1439: PetscCall(VecDestroy(&F_alg));
1440: PetscCall(VecDestroy(&X));
1441: PetscCall(TSDestroy(&ts));
1442: for (PetscInt i = 0; i < 3; i++) PetscCall(VecDestroy(&DICDP[i]));
1443: PetscFunctionReturn(PETSC_SUCCESS);
1444: }
1446: /*TEST
1448: build:
1449: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1451: test:
1452: args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1453: localrunfiles: petscoptions X.bin Ybus.bin
1455: TEST*/