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, j;
331: PetscFunctionBegin;
332: eps = 1.e-7;
333: PetscCall(VecDuplicate(X, &Y));
335: for (i = 0; i < ngen; i++) {
336: for (j = 0; j < 3; j++) PGv[j] = PG[j];
337: PGv[i] = PG[i] + eps;
338: PetscCall(InitialGuess(Y, user, PGv));
339: PetscCall(InitialGuess(X, user, PG));
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;
363: PetscInt k;
365: PetscFunctionBegin;
366: PetscCall(VecZeroEntries(F));
367: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
368: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
369: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
370: PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));
372: /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
373: The generator current injection, IG, and load current injection, ID are added later
374: */
375: /* Note that the values in Ybus are stored assuming the imaginary current balance
376: equation is ordered first followed by real current balance equation for each bus.
377: Thus imaginary current contribution goes in location 2*i, and
378: real current contribution in 2*i+1
379: */
380: PetscCall(MatMult(user->Ybus, Xnet, Fnet));
382: PetscCall(VecGetArray(Xgen, &xgen));
383: PetscCall(VecGetArray(Xnet, &xnet));
384: PetscCall(VecGetArray(Fgen, &fgen));
385: PetscCall(VecGetArray(Fnet, &fnet));
387: /* Generator subsystem */
388: for (i = 0; i < ngen; i++) {
389: Eqp = xgen[idx];
390: Edp = xgen[idx + 1];
391: delta = xgen[idx + 2];
392: w = xgen[idx + 3];
393: Id = xgen[idx + 4];
394: Iq = xgen[idx + 5];
395: Efd = xgen[idx + 6];
396: RF = xgen[idx + 7];
397: VR = xgen[idx + 8];
399: /* Generator differential equations */
400: fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
401: fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
402: fgen[idx + 2] = -w + w_s;
403: fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
405: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
406: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
408: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
409: /* Algebraic equations for stator currents */
410: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
412: Zdq_inv[0] = Rs[i] / det;
413: Zdq_inv[1] = Xqp[i] / det;
414: Zdq_inv[2] = -Xdp[i] / det;
415: Zdq_inv[3] = Rs[i] / det;
417: fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
418: fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
420: /* Add generator current injection to network */
421: PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
423: fnet[2 * gbus[i]] -= IGi;
424: fnet[2 * gbus[i] + 1] -= IGr;
426: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
428: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
430: /* Exciter differential equations */
431: fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
432: fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
433: fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
435: idx = idx + 9;
436: }
438: PetscCall(VecGetArray(user->V0, &v0));
439: for (i = 0; i < nload; i++) {
440: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
441: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
442: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
443: Vm2 = Vm * Vm;
444: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
445: PD = QD = 0.0;
446: for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
447: for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
449: /* Load currents */
450: IDr = (PD * Vr + QD * Vi) / Vm2;
451: IDi = (-QD * Vr + PD * Vi) / Vm2;
453: fnet[2 * lbus[i]] += IDi;
454: fnet[2 * lbus[i] + 1] += IDr;
455: }
456: PetscCall(VecRestoreArray(user->V0, &v0));
458: PetscCall(VecRestoreArray(Xgen, &xgen));
459: PetscCall(VecRestoreArray(Xnet, &xnet));
460: PetscCall(VecRestoreArray(Fgen, &fgen));
461: PetscCall(VecRestoreArray(Fnet, &fnet));
463: PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
464: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
465: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /* \dot{x} - f(x,y)
470: g(x,y) = 0
471: */
472: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
473: {
474: SNES snes;
475: PetscScalar *f;
476: const PetscScalar *xdot;
477: PetscInt i;
479: PetscFunctionBegin;
480: user->t = t;
482: PetscCall(TSGetSNES(ts, &snes));
483: PetscCall(ResidualFunction(snes, X, F, user));
484: PetscCall(VecGetArray(F, &f));
485: PetscCall(VecGetArrayRead(Xdot, &xdot));
486: for (i = 0; i < ngen; i++) {
487: f[9 * i] += xdot[9 * i];
488: f[9 * i + 1] += xdot[9 * i + 1];
489: f[9 * i + 2] += xdot[9 * i + 2];
490: f[9 * i + 3] += xdot[9 * i + 3];
491: f[9 * i + 6] += xdot[9 * i + 6];
492: f[9 * i + 7] += xdot[9 * i + 7];
493: f[9 * i + 8] += xdot[9 * i + 8];
494: }
495: PetscCall(VecRestoreArray(F, &f));
496: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /* This function is used for solving the algebraic system only during fault on and
501: off times. It computes the entire F and then zeros out the part corresponding to
502: differential equations
503: F = [0;g(y)];
504: */
505: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
506: {
507: Userctx *user = (Userctx *)ctx;
508: PetscScalar *f;
509: PetscInt i;
511: PetscFunctionBegin;
512: PetscCall(ResidualFunction(snes, X, F, user));
513: PetscCall(VecGetArray(F, &f));
514: for (i = 0; i < ngen; i++) {
515: f[9 * i] = 0;
516: f[9 * i + 1] = 0;
517: f[9 * i + 2] = 0;
518: f[9 * i + 3] = 0;
519: f[9 * i + 6] = 0;
520: f[9 * i + 7] = 0;
521: f[9 * i + 8] = 0;
522: }
523: PetscCall(VecRestoreArray(F, &f));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
528: {
529: PetscInt *d_nnz;
530: PetscInt i, idx = 0, start = 0;
531: PetscInt ncols;
533: PetscFunctionBegin;
534: PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
535: for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
536: /* Generator subsystem */
537: for (i = 0; i < ngen; i++) {
538: d_nnz[idx] += 3;
539: d_nnz[idx + 1] += 2;
540: d_nnz[idx + 2] += 2;
541: d_nnz[idx + 3] += 5;
542: d_nnz[idx + 4] += 6;
543: d_nnz[idx + 5] += 6;
545: d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
546: d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
548: d_nnz[idx + 6] += 2;
549: d_nnz[idx + 7] += 2;
550: d_nnz[idx + 8] += 5;
552: idx = idx + 9;
553: }
555: start = user->neqs_gen;
556: for (i = 0; i < nbus; i++) {
557: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
558: d_nnz[start + 2 * i] += ncols;
559: d_nnz[start + 2 * i + 1] += ncols;
560: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
561: }
563: PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
564: PetscCall(PetscFree(d_nnz));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*
569: J = [-df_dx, -df_dy
570: dg_dx, dg_dy]
571: */
572: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, void *ctx)
573: {
574: Userctx *user = (Userctx *)ctx;
575: Vec Xgen, Xnet;
576: PetscScalar *xgen, *xnet;
577: PetscInt i, idx = 0;
578: PetscScalar Vr, Vi, Vm, Vm2;
579: PetscScalar Eqp, Edp, delta; /* Generator variables */
580: PetscScalar Efd; /* Exciter variables */
581: PetscScalar Id, Iq; /* Generator dq axis currents */
582: PetscScalar Vd, Vq;
583: PetscScalar val[10];
584: PetscInt row[2], col[10];
585: PetscInt net_start = user->neqs_gen;
586: PetscInt ncols;
587: const PetscInt *cols;
588: const PetscScalar *yvals;
589: PetscInt k;
590: PetscScalar Zdq_inv[4], det;
591: PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
592: PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
593: PetscScalar dSE_dEfd;
594: PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
595: PetscScalar PD, QD, Vm0, *v0, Vm4;
596: PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
597: PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
599: PetscFunctionBegin;
600: PetscCall(MatZeroEntries(B));
601: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
602: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
604: PetscCall(VecGetArray(Xgen, &xgen));
605: PetscCall(VecGetArray(Xnet, &xnet));
607: /* Generator subsystem */
608: for (i = 0; i < ngen; i++) {
609: Eqp = xgen[idx];
610: Edp = xgen[idx + 1];
611: delta = xgen[idx + 2];
612: Id = xgen[idx + 4];
613: Iq = xgen[idx + 5];
614: Efd = xgen[idx + 6];
616: /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
617: row[0] = idx;
618: col[0] = idx;
619: col[1] = idx + 4;
620: col[2] = idx + 6;
621: val[0] = 1 / Td0p[i];
622: val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
623: val[2] = -1 / Td0p[i];
625: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
627: /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
628: row[0] = idx + 1;
629: col[0] = idx + 1;
630: col[1] = idx + 5;
631: val[0] = 1 / Tq0p[i];
632: val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
633: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
635: /* fgen[idx+2] = - w + w_s; */
636: row[0] = idx + 2;
637: col[0] = idx + 2;
638: col[1] = idx + 3;
639: val[0] = 0;
640: val[1] = -1;
641: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
643: /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
644: row[0] = idx + 3;
645: col[0] = idx;
646: col[1] = idx + 1;
647: col[2] = idx + 3;
648: col[3] = idx + 4;
649: col[4] = idx + 5;
650: val[0] = Iq / M[i];
651: val[1] = Id / M[i];
652: val[2] = D[i] / M[i];
653: val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
654: val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
655: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
657: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
658: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
659: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
661: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
663: Zdq_inv[0] = Rs[i] / det;
664: Zdq_inv[1] = Xqp[i] / det;
665: Zdq_inv[2] = -Xdp[i] / det;
666: Zdq_inv[3] = Rs[i] / det;
668: dVd_dVr = PetscSinScalar(delta);
669: dVd_dVi = -PetscCosScalar(delta);
670: dVq_dVr = PetscCosScalar(delta);
671: dVq_dVi = PetscSinScalar(delta);
672: dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
673: dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
675: /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
676: row[0] = idx + 4;
677: col[0] = idx;
678: col[1] = idx + 1;
679: col[2] = idx + 2;
680: val[0] = -Zdq_inv[1];
681: val[1] = -Zdq_inv[0];
682: val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
683: col[3] = idx + 4;
684: col[4] = net_start + 2 * gbus[i];
685: col[5] = net_start + 2 * gbus[i] + 1;
686: val[3] = 1;
687: val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
688: val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
689: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
691: /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
692: row[0] = idx + 5;
693: col[0] = idx;
694: col[1] = idx + 1;
695: col[2] = idx + 2;
696: val[0] = -Zdq_inv[3];
697: val[1] = -Zdq_inv[2];
698: val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
699: col[3] = idx + 5;
700: col[4] = net_start + 2 * gbus[i];
701: col[5] = net_start + 2 * gbus[i] + 1;
702: val[3] = 1;
703: val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
704: val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
705: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
707: dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
708: dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
709: dIGr_dId = PetscSinScalar(delta);
710: dIGr_dIq = PetscCosScalar(delta);
711: dIGi_dId = -PetscCosScalar(delta);
712: dIGi_dIq = PetscSinScalar(delta);
714: /* fnet[2*gbus[i]] -= IGi; */
715: row[0] = net_start + 2 * gbus[i];
716: col[0] = idx + 2;
717: col[1] = idx + 4;
718: col[2] = idx + 5;
719: val[0] = -dIGi_ddelta;
720: val[1] = -dIGi_dId;
721: val[2] = -dIGi_dIq;
722: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
724: /* fnet[2*gbus[i]+1] -= IGr; */
725: row[0] = net_start + 2 * gbus[i] + 1;
726: col[0] = idx + 2;
727: col[1] = idx + 4;
728: col[2] = idx + 5;
729: val[0] = -dIGr_ddelta;
730: val[1] = -dIGr_dId;
731: val[2] = -dIGr_dIq;
732: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
734: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
736: /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
737: /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
738: dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
740: row[0] = idx + 6;
741: col[0] = idx + 6;
742: col[1] = idx + 8;
743: val[0] = (KE[i] + dSE_dEfd) / TE[i];
744: val[1] = -1 / TE[i];
745: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
747: /* Exciter differential equations */
749: /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
750: row[0] = idx + 7;
751: col[0] = idx + 6;
752: col[1] = idx + 7;
753: val[0] = (-KF[i] / TF[i]) / TF[i];
754: val[1] = 1 / TF[i];
755: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
757: /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
758: /* Vm = (Vd^2 + Vq^2)^0.5; */
759: dVm_dVd = Vd / Vm;
760: dVm_dVq = Vq / Vm;
761: dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
762: dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
763: row[0] = idx + 8;
764: col[0] = idx + 6;
765: col[1] = idx + 7;
766: col[2] = idx + 8;
767: val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
768: val[1] = -KA[i] / TA[i];
769: val[2] = 1 / TA[i];
770: col[3] = net_start + 2 * gbus[i];
771: col[4] = net_start + 2 * gbus[i] + 1;
772: val[3] = KA[i] * dVm_dVr / TA[i];
773: val[4] = KA[i] * dVm_dVi / TA[i];
774: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
775: idx = idx + 9;
776: }
778: for (i = 0; i < nbus; i++) {
779: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
780: row[0] = net_start + 2 * i;
781: for (k = 0; k < ncols; k++) {
782: col[k] = net_start + cols[k];
783: val[k] = yvals[k];
784: }
785: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
786: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
788: PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
789: row[0] = net_start + 2 * i + 1;
790: for (k = 0; k < ncols; k++) {
791: col[k] = net_start + cols[k];
792: val[k] = yvals[k];
793: }
794: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
795: PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
796: }
798: PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
799: PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));
801: PetscCall(VecGetArray(user->V0, &v0));
802: for (i = 0; i < nload; i++) {
803: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
804: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
805: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
806: Vm2 = Vm * Vm;
807: Vm4 = Vm2 * Vm2;
808: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
809: PD = QD = 0.0;
810: dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
811: for (k = 0; k < ld_nsegsp[i]; k++) {
812: PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
813: dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vr * PetscPowScalar(Vm, (ld_betap[k] - 2));
814: dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vi * PetscPowScalar(Vm, (ld_betap[k] - 2));
815: }
816: for (k = 0; k < ld_nsegsq[i]; k++) {
817: QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
818: dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vr * PetscPowScalar(Vm, (ld_betaq[k] - 2));
819: dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vi * PetscPowScalar(Vm, (ld_betaq[k] - 2));
820: }
822: /* IDr = (PD*Vr + QD*Vi)/Vm2; */
823: /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
825: dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
826: dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
828: dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
829: dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
831: /* fnet[2*lbus[i]] += IDi; */
832: row[0] = net_start + 2 * lbus[i];
833: col[0] = net_start + 2 * lbus[i];
834: col[1] = net_start + 2 * lbus[i] + 1;
835: val[0] = dIDi_dVr;
836: val[1] = dIDi_dVi;
837: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
838: /* fnet[2*lbus[i]+1] += IDr; */
839: row[0] = net_start + 2 * lbus[i] + 1;
840: col[0] = net_start + 2 * lbus[i];
841: col[1] = net_start + 2 * lbus[i] + 1;
842: val[0] = dIDr_dVr;
843: val[1] = dIDr_dVi;
844: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
845: }
846: PetscCall(VecRestoreArray(user->V0, &v0));
848: PetscCall(VecRestoreArray(Xgen, &xgen));
849: PetscCall(VecRestoreArray(Xnet, &xnet));
851: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
853: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
854: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: /*
859: J = [I, 0
860: dg_dx, dg_dy]
861: */
862: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, void *ctx)
863: {
864: Userctx *user = (Userctx *)ctx;
866: PetscFunctionBegin;
867: PetscCall(ResidualJacobian(snes, X, A, B, ctx));
868: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
869: PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*
874: J = [a*I-df_dx, -df_dy
875: dg_dx, dg_dy]
876: */
878: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
879: {
880: SNES snes;
881: PetscScalar atmp = (PetscScalar)a;
882: PetscInt i, row;
884: PetscFunctionBegin;
885: user->t = t;
887: PetscCall(TSGetSNES(ts, &snes));
888: PetscCall(ResidualJacobian(snes, X, A, B, user));
889: for (i = 0; i < ngen; i++) {
890: row = 9 * i;
891: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
892: row = 9 * i + 1;
893: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
894: row = 9 * i + 2;
895: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
896: row = 9 * i + 3;
897: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
898: row = 9 * i + 6;
899: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
900: row = 9 * i + 7;
901: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
902: row = 9 * i + 8;
903: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
904: }
905: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
906: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /* Matrix JacobianP is constant so that it only needs to be evaluated once */
911: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
912: {
913: PetscScalar a;
914: PetscInt row, col;
915: Userctx *ctx = (Userctx *)ctx0;
917: PetscFunctionBeginUser;
919: if (ctx->jacp_flg) {
920: PetscCall(MatZeroEntries(A));
922: for (col = 0; col < 3; col++) {
923: a = 1.0 / M[col];
924: row = 9 * col + 3;
925: PetscCall(MatSetValues(A, 1, &row, 1, &col, &a, INSERT_VALUES));
926: }
928: ctx->jacp_flg = PETSC_FALSE;
930: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
931: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
932: }
933: PetscFunctionReturn(PETSC_SUCCESS);
934: }
936: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
937: {
938: const PetscScalar *u;
939: PetscInt idx;
940: Vec Xgen, Xnet;
941: PetscScalar *r, *xgen;
942: PetscInt i;
944: PetscFunctionBegin;
945: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
946: PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
948: PetscCall(VecGetArray(Xgen, &xgen));
950: PetscCall(VecGetArrayRead(U, &u));
951: PetscCall(VecGetArray(R, &r));
952: r[0] = 0.;
953: idx = 0;
954: for (i = 0; i < ngen; i++) {
955: 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);
956: idx += 9;
957: }
958: PetscCall(VecRestoreArrayRead(U, &u));
959: PetscCall(VecRestoreArray(R, &r));
960: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
961: PetscFunctionReturn(PETSC_SUCCESS);
962: }
964: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, Userctx *user)
965: {
966: Vec Xgen, Xnet, Dgen, Dnet;
967: PetscScalar *xgen, *dgen;
968: PetscInt i;
969: PetscInt idx;
970: Vec drdu_col;
971: PetscScalar *xarr;
973: PetscFunctionBegin;
974: PetscCall(VecDuplicate(U, &drdu_col));
975: PetscCall(MatDenseGetColumn(DRDU, 0, &xarr));
976: PetscCall(VecPlaceArray(drdu_col, xarr));
977: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
978: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Dgen, &Dnet));
979: PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
980: PetscCall(DMCompositeScatter(user->dmpgrid, drdu_col, Dgen, Dnet));
982: PetscCall(VecGetArray(Xgen, &xgen));
983: PetscCall(VecGetArray(Dgen, &dgen));
985: idx = 0;
986: for (i = 0; i < ngen; i++) {
987: dgen[idx + 3] = 0.;
988: 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);
989: 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);
990: idx += 9;
991: }
993: PetscCall(VecRestoreArray(Dgen, &dgen));
994: PetscCall(VecRestoreArray(Xgen, &xgen));
995: PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, drdu_col, Dgen, Dnet));
996: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Dgen, &Dnet));
997: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
998: PetscCall(VecResetArray(drdu_col));
999: PetscCall(MatDenseRestoreColumn(DRDU, &xarr));
1000: PetscCall(VecDestroy(&drdu_col));
1001: PetscFunctionReturn(PETSC_SUCCESS);
1002: }
1004: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat drdp, Userctx *user)
1005: {
1006: PetscFunctionBegin;
1007: PetscFunctionReturn(PETSC_SUCCESS);
1008: }
1010: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, Vec *DICDP, Userctx *user)
1011: {
1012: PetscScalar *x, *y, sensip;
1013: PetscInt i;
1015: PetscFunctionBegin;
1016: PetscCall(VecGetArray(lambda, &x));
1017: PetscCall(VecGetArray(mu, &y));
1019: for (i = 0; i < 3; i++) {
1020: PetscCall(VecDot(lambda, DICDP[i], &sensip));
1021: sensip = sensip + y[i];
1022: /* PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt %" PetscInt_FMT " th parameter: %g \n",i,(double)sensip)); */
1023: y[i] = sensip;
1024: }
1025: PetscCall(VecRestoreArray(mu, &y));
1026: PetscFunctionReturn(PETSC_SUCCESS);
1027: }
1029: int main(int argc, char **argv)
1030: {
1031: Userctx user;
1032: Vec p;
1033: PetscScalar *x_ptr;
1034: PetscMPIInt size;
1035: PetscInt i;
1036: PetscViewer Xview, Ybusview;
1037: PetscInt *idx2;
1038: Tao tao;
1039: KSP ksp;
1040: PC pc;
1041: Vec lowerb, upperb;
1043: PetscFunctionBeginUser;
1044: PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
1045: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1046: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
1048: user.jacp_flg = PETSC_TRUE;
1049: user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
1050: user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
1051: user.neqs_pgrid = user.neqs_gen + user.neqs_net;
1053: /* Create indices for differential and algebraic equations */
1054: PetscCall(PetscMalloc1(7 * ngen, &idx2));
1055: for (i = 0; i < ngen; i++) {
1056: idx2[7 * i] = 9 * i;
1057: idx2[7 * i + 1] = 9 * i + 1;
1058: idx2[7 * i + 2] = 9 * i + 2;
1059: idx2[7 * i + 3] = 9 * i + 3;
1060: idx2[7 * i + 4] = 9 * i + 6;
1061: idx2[7 * i + 5] = 9 * i + 7;
1062: idx2[7 * i + 6] = 9 * i + 8;
1063: }
1064: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
1065: PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
1066: PetscCall(PetscFree(idx2));
1068: /* Set run time options */
1069: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1070: {
1071: user.tfaulton = 1.0;
1072: user.tfaultoff = 1.2;
1073: user.Rfault = 0.0001;
1074: user.faultbus = 8;
1075: PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
1076: PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
1077: PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
1078: user.t0 = 0.0;
1079: user.tmax = 1.3;
1080: PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
1081: PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
1082: user.freq_u = 61.0;
1083: user.freq_l = 59.0;
1084: user.pow = 2;
1085: PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
1086: PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
1087: PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
1088: }
1089: PetscOptionsEnd();
1091: /* Create DMs for generator and network subsystems */
1092: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
1093: PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
1094: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
1095: PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
1096: PetscCall(DMSetFromOptions(user.dmnet));
1097: PetscCall(DMSetUp(user.dmnet));
1098: /* Create a composite DM packer and add the two DMs */
1099: PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
1100: PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
1101: PetscCall(DMSetFromOptions(user.dmgen));
1102: PetscCall(DMSetUp(user.dmgen));
1103: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
1104: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));
1106: /* Read initial voltage vector and Ybus */
1107: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
1108: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));
1110: PetscCall(VecCreate(PETSC_COMM_WORLD, &user.V0));
1111: PetscCall(VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net));
1112: PetscCall(VecLoad(user.V0, Xview));
1114: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Ybus));
1115: PetscCall(MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net));
1116: PetscCall(MatSetType(user.Ybus, MATBAIJ));
1117: /* PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
1118: PetscCall(MatLoad(user.Ybus, Ybusview));
1120: PetscCall(PetscViewerDestroy(&Xview));
1121: PetscCall(PetscViewerDestroy(&Ybusview));
1123: /* Allocate space for Jacobians */
1124: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.J));
1125: PetscCall(MatSetSizes(user.J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid));
1126: PetscCall(MatSetFromOptions(user.J));
1127: PetscCall(PreallocateJacobian(user.J, &user));
1129: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
1130: PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 3));
1131: PetscCall(MatSetFromOptions(user.Jacp));
1132: PetscCall(MatSetUp(user.Jacp));
1133: PetscCall(MatZeroEntries(user.Jacp)); /* initialize to zeros */
1135: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, &user.DRDP));
1136: PetscCall(MatSetUp(user.DRDP));
1137: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 1, NULL, &user.DRDU));
1138: PetscCall(MatSetUp(user.DRDU));
1140: /* Create TAO solver and set desired solution method */
1141: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
1142: PetscCall(TaoSetType(tao, TAOBLMVM));
1143: /*
1144: Optimization starts
1145: */
1146: /* Set initial solution guess */
1147: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
1148: PetscCall(VecGetArray(p, &x_ptr));
1149: x_ptr[0] = PG[0];
1150: x_ptr[1] = PG[1];
1151: x_ptr[2] = PG[2];
1152: PetscCall(VecRestoreArray(p, &x_ptr));
1154: PetscCall(TaoSetSolution(tao, p));
1155: /* Set routine for function and gradient evaluation */
1156: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user));
1158: /* Set bounds for the optimization */
1159: PetscCall(VecDuplicate(p, &lowerb));
1160: PetscCall(VecDuplicate(p, &upperb));
1161: PetscCall(VecGetArray(lowerb, &x_ptr));
1162: x_ptr[0] = 0.5;
1163: x_ptr[1] = 0.5;
1164: x_ptr[2] = 0.5;
1165: PetscCall(VecRestoreArray(lowerb, &x_ptr));
1166: PetscCall(VecGetArray(upperb, &x_ptr));
1167: x_ptr[0] = 2.0;
1168: x_ptr[1] = 2.0;
1169: x_ptr[2] = 2.0;
1170: PetscCall(VecRestoreArray(upperb, &x_ptr));
1171: PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));
1173: /* Check for any TAO command line options */
1174: PetscCall(TaoSetFromOptions(tao));
1175: PetscCall(TaoGetKSP(tao, &ksp));
1176: if (ksp) {
1177: PetscCall(KSPGetPC(ksp, &pc));
1178: PetscCall(PCSetType(pc, PCNONE));
1179: }
1181: /* SOLVE THE APPLICATION */
1182: PetscCall(TaoSolve(tao));
1184: PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
1185: /* Free TAO data structures */
1186: PetscCall(TaoDestroy(&tao));
1188: PetscCall(DMDestroy(&user.dmgen));
1189: PetscCall(DMDestroy(&user.dmnet));
1190: PetscCall(DMDestroy(&user.dmpgrid));
1191: PetscCall(ISDestroy(&user.is_diff));
1192: PetscCall(ISDestroy(&user.is_alg));
1194: PetscCall(MatDestroy(&user.J));
1195: PetscCall(MatDestroy(&user.Jacp));
1196: PetscCall(MatDestroy(&user.Ybus));
1197: /* PetscCall(MatDestroy(&user.Sol)); */
1198: PetscCall(VecDestroy(&user.V0));
1199: PetscCall(VecDestroy(&p));
1200: PetscCall(VecDestroy(&lowerb));
1201: PetscCall(VecDestroy(&upperb));
1202: PetscCall(MatDestroy(&user.DRDU));
1203: PetscCall(MatDestroy(&user.DRDP));
1204: PetscCall(PetscFinalize());
1205: return 0;
1206: }
1208: /* ------------------------------------------------------------------ */
1209: /*
1210: FormFunction - Evaluates the function and corresponding gradient.
1212: Input Parameters:
1213: tao - the Tao context
1214: X - the input vector
1215: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
1217: Output Parameters:
1218: f - the newly evaluated function
1219: G - the newly evaluated gradient
1220: */
1221: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx0)
1222: {
1223: TS ts, quadts;
1224: SNES snes_alg;
1225: Userctx *ctx = (Userctx *)ctx0;
1226: Vec X;
1227: PetscInt i;
1228: /* sensitivity context */
1229: PetscScalar *x_ptr;
1230: Vec lambda[1], q;
1231: Vec mu[1];
1232: PetscInt steps1, steps2, steps3;
1233: Vec DICDP[3];
1234: Vec F_alg;
1235: PetscInt row_loc, col_loc;
1236: PetscScalar val;
1237: Vec Xdot;
1239: PetscFunctionBegin;
1240: PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
1241: PG[0] = x_ptr[0];
1242: PG[1] = x_ptr[1];
1243: PG[2] = x_ptr[2];
1244: PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
1246: ctx->stepnum = 0;
1248: PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));
1250: /* Create matrix to save solutions at each time step */
1251: /* PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ctx->neqs_pgrid+1,1002,NULL,&ctx->Sol)); */
1252: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1253: Create timestepping solver context
1254: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1255: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1256: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1257: PetscCall(TSSetType(ts, TSCN));
1258: PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)IFunction, ctx));
1259: PetscCall(TSSetIJacobian(ts, ctx->J, ctx->J, (TSIJacobian)IJacobian, ctx));
1260: PetscCall(TSSetApplicationContext(ts, ctx));
1261: /* Set RHS JacobianP */
1262: PetscCall(TSSetRHSJacobianP(ts, ctx->Jacp, RHSJacobianP, ctx));
1264: PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
1265: PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunction)CostIntegrand, ctx));
1266: PetscCall(TSSetRHSJacobian(quadts, ctx->DRDU, ctx->DRDU, (TSRHSJacobian)DRDUJacobianTranspose, ctx));
1267: PetscCall(TSSetRHSJacobianP(quadts, ctx->DRDP, (TSRHSJacobianP)DRDPJacobianTranspose, ctx));
1269: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1270: Set initial conditions
1271: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1272: PetscCall(SetInitialGuess(X, ctx));
1274: /* Approximate DICDP with finite difference, we want to zero out network variables */
1275: for (i = 0; i < 3; i++) PetscCall(VecDuplicate(X, &DICDP[i]));
1276: PetscCall(DICDPFiniteDifference(X, DICDP, ctx));
1278: PetscCall(VecDuplicate(X, &F_alg));
1279: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1280: PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
1281: PetscCall(MatZeroEntries(ctx->J));
1282: PetscCall(SNESSetJacobian(snes_alg, ctx->J, ctx->J, AlgJacobian, ctx));
1283: PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1284: PetscCall(SNESSetFromOptions(snes_alg));
1285: ctx->alg_flg = PETSC_TRUE;
1286: /* Solve the algebraic equations */
1287: PetscCall(SNESSolve(snes_alg, NULL, X));
1289: /* Just to set up the Jacobian structure */
1290: PetscCall(VecDuplicate(X, &Xdot));
1291: PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, ctx->J, ctx->J, ctx));
1292: PetscCall(VecDestroy(&Xdot));
1294: ctx->stepnum++;
1296: /*
1297: Save trajectory of solution so that TSAdjointSolve() may be used
1298: */
1299: PetscCall(TSSetSaveTrajectory(ts));
1301: PetscCall(TSSetTimeStep(ts, 0.01));
1302: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1303: PetscCall(TSSetFromOptions(ts));
1304: /* PetscCall(TSSetPostStep(ts,SaveSolution)); */
1306: /* Prefault period */
1307: ctx->alg_flg = PETSC_FALSE;
1308: PetscCall(TSSetTime(ts, 0.0));
1309: PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
1310: PetscCall(TSSolve(ts, X));
1311: PetscCall(TSGetStepNumber(ts, &steps1));
1313: /* Create the nonlinear solver for solving the algebraic system */
1314: /* Note that although the algebraic system needs to be solved only for
1315: Idq and V, we reuse the entire system including xgen. The xgen
1316: variables are held constant by setting their residuals to 0 and
1317: putting a 1 on the Jacobian diagonal for xgen rows
1318: */
1319: PetscCall(MatZeroEntries(ctx->J));
1321: /* Apply disturbance - resistive fault at ctx->faultbus */
1322: /* This is done by adding shunt conductance to the diagonal location
1323: in the Ybus matrix */
1324: row_loc = 2 * ctx->faultbus;
1325: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1326: val = 1 / ctx->Rfault;
1327: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1328: row_loc = 2 * ctx->faultbus + 1;
1329: col_loc = 2 * ctx->faultbus; /* Location for G */
1330: val = 1 / ctx->Rfault;
1331: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1333: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1334: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1336: ctx->alg_flg = PETSC_TRUE;
1337: /* Solve the algebraic equations */
1338: PetscCall(SNESSolve(snes_alg, NULL, X));
1340: ctx->stepnum++;
1342: /* Disturbance period */
1343: ctx->alg_flg = PETSC_FALSE;
1344: PetscCall(TSSetTime(ts, ctx->tfaulton));
1345: PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
1346: PetscCall(TSSolve(ts, X));
1347: PetscCall(TSGetStepNumber(ts, &steps2));
1348: steps2 -= steps1;
1350: /* Remove the fault */
1351: row_loc = 2 * ctx->faultbus;
1352: col_loc = 2 * ctx->faultbus + 1;
1353: val = -1 / ctx->Rfault;
1354: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1355: row_loc = 2 * ctx->faultbus + 1;
1356: col_loc = 2 * ctx->faultbus;
1357: val = -1 / ctx->Rfault;
1358: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1360: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1361: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1363: PetscCall(MatZeroEntries(ctx->J));
1365: ctx->alg_flg = PETSC_TRUE;
1367: /* Solve the algebraic equations */
1368: PetscCall(SNESSolve(snes_alg, NULL, X));
1370: ctx->stepnum++;
1372: /* Post-disturbance period */
1373: ctx->alg_flg = PETSC_TRUE;
1374: PetscCall(TSSetTime(ts, ctx->tfaultoff));
1375: PetscCall(TSSetMaxTime(ts, ctx->tmax));
1376: PetscCall(TSSolve(ts, X));
1377: PetscCall(TSGetStepNumber(ts, &steps3));
1378: steps3 -= steps2;
1379: steps3 -= steps1;
1381: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1382: Adjoint model starts here
1383: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1384: PetscCall(TSSetPostStep(ts, NULL));
1385: PetscCall(MatCreateVecs(ctx->J, &lambda[0], NULL));
1386: /* Set initial conditions for the adjoint integration */
1387: PetscCall(VecZeroEntries(lambda[0]));
1389: PetscCall(MatCreateVecs(ctx->Jacp, &mu[0], NULL));
1390: PetscCall(VecZeroEntries(mu[0]));
1391: PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
1393: PetscCall(TSAdjointSetSteps(ts, steps3));
1394: PetscCall(TSAdjointSolve(ts));
1396: PetscCall(MatZeroEntries(ctx->J));
1397: /* Applying disturbance - resistive fault at ctx->faultbus */
1398: /* This is done by deducting shunt conductance to the diagonal location
1399: in the Ybus matrix */
1400: row_loc = 2 * ctx->faultbus;
1401: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1402: val = 1. / ctx->Rfault;
1403: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1404: row_loc = 2 * ctx->faultbus + 1;
1405: col_loc = 2 * ctx->faultbus; /* Location for G */
1406: val = 1. / ctx->Rfault;
1407: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1409: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1410: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1412: /* Set number of steps for the adjoint integration */
1413: PetscCall(TSAdjointSetSteps(ts, steps2));
1414: PetscCall(TSAdjointSolve(ts));
1416: PetscCall(MatZeroEntries(ctx->J));
1417: /* remove the fault */
1418: row_loc = 2 * ctx->faultbus;
1419: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1420: val = -1. / ctx->Rfault;
1421: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1422: row_loc = 2 * ctx->faultbus + 1;
1423: col_loc = 2 * ctx->faultbus; /* Location for G */
1424: val = -1. / ctx->Rfault;
1425: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1427: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1428: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1430: /* Set number of steps for the adjoint integration */
1431: PetscCall(TSAdjointSetSteps(ts, steps1));
1432: PetscCall(TSAdjointSolve(ts));
1434: PetscCall(ComputeSensiP(lambda[0], mu[0], DICDP, ctx));
1435: PetscCall(VecCopy(mu[0], G));
1437: PetscCall(TSGetQuadratureTS(ts, NULL, &quadts));
1438: PetscCall(TSGetSolution(quadts, &q));
1439: PetscCall(VecGetArray(q, &x_ptr));
1440: *f = x_ptr[0];
1441: x_ptr[0] = 0;
1442: PetscCall(VecRestoreArray(q, &x_ptr));
1444: PetscCall(VecDestroy(&lambda[0]));
1445: PetscCall(VecDestroy(&mu[0]));
1447: PetscCall(SNESDestroy(&snes_alg));
1448: PetscCall(VecDestroy(&F_alg));
1449: PetscCall(VecDestroy(&X));
1450: PetscCall(TSDestroy(&ts));
1451: for (i = 0; i < 3; i++) PetscCall(VecDestroy(&DICDP[i]));
1452: PetscFunctionReturn(PETSC_SUCCESS);
1453: }
1455: /*TEST
1457: build:
1458: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1460: test:
1461: args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1462: localrunfiles: petscoptions X.bin Ybus.bin
1464: TEST*/