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: *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
114: *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
115: return 0;
116: }
118: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
119: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
120: {
121: *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
122: *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
123: return 0;
124: }
126: /* Saves the solution at each time to a matrix */
127: PetscErrorCode SaveSolution(TS ts)
128: {
129: Userctx *user;
130: Vec X;
131: PetscScalar *mat;
132: const PetscScalar *x;
133: PetscInt idx;
134: PetscReal t;
136: TSGetApplicationContext(ts, &user);
137: TSGetTime(ts, &t);
138: TSGetSolution(ts, &X);
139: idx = user->stepnum * (user->neqs_pgrid + 1);
140: MatDenseGetArray(user->Sol, &mat);
141: VecGetArrayRead(X, &x);
142: mat[idx] = t;
143: PetscArraycpy(mat + idx + 1, x, user->neqs_pgrid);
144: MatDenseRestoreArray(user->Sol, &mat);
145: VecRestoreArrayRead(X, &x);
146: user->stepnum++;
147: return 0;
148: }
150: PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
151: {
152: Vec Xgen, Xnet;
153: PetscScalar *xgen, *xnet;
154: PetscInt i, idx = 0;
155: PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
156: PetscScalar Eqp, Edp, delta;
157: PetscScalar Efd, RF, VR; /* Exciter variables */
158: PetscScalar Id, Iq; /* Generator dq axis currents */
159: PetscScalar theta, Vd, Vq, SE;
161: M[0] = 2 * H[0] / w_s;
162: M[1] = 2 * H[1] / w_s;
163: M[2] = 2 * H[2] / w_s;
164: D[0] = 0.1 * M[0];
165: D[1] = 0.1 * M[1];
166: D[2] = 0.1 * M[2];
168: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
170: /* Network subsystem initialization */
171: VecCopy(user->V0, Xnet);
173: /* Generator subsystem initialization */
174: VecGetArray(Xgen, &xgen);
175: VecGetArray(Xnet, &xnet);
177: for (i = 0; i < ngen; i++) {
178: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
179: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
180: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
181: Vm2 = Vm * Vm;
182: IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
183: IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;
185: delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
187: theta = PETSC_PI / 2.0 - delta;
189: Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
190: Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
192: Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
193: Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
195: Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
196: Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
198: TM[i] = PG[i];
200: /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
201: xgen[idx] = Eqp;
202: xgen[idx + 1] = Edp;
203: xgen[idx + 2] = delta;
204: xgen[idx + 3] = w_s;
206: idx = idx + 4;
208: xgen[idx] = Id;
209: xgen[idx + 1] = Iq;
211: idx = idx + 2;
213: /* Exciter */
214: Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
215: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
216: VR = KE[i] * Efd + SE;
217: RF = KF[i] * Efd / TF[i];
219: xgen[idx] = Efd;
220: xgen[idx + 1] = RF;
221: xgen[idx + 2] = VR;
223: Vref[i] = Vm + (VR / KA[i]);
225: idx = idx + 3;
226: }
228: VecRestoreArray(Xgen, &xgen);
229: VecRestoreArray(Xnet, &xnet);
231: /* VecView(Xgen,0); */
232: DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet);
233: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
234: return 0;
235: }
237: PetscErrorCode InitialGuess(Vec X, Userctx *user, const PetscScalar PGv[])
238: {
239: Vec Xgen, Xnet;
240: PetscScalar *xgen, *xnet;
241: PetscInt i, idx = 0;
242: PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
243: PetscScalar Eqp, Edp, delta;
244: PetscScalar Efd, RF, VR; /* Exciter variables */
245: PetscScalar Id, Iq; /* Generator dq axis currents */
246: PetscScalar theta, Vd, Vq, SE;
248: M[0] = 2 * H[0] / w_s;
249: M[1] = 2 * H[1] / w_s;
250: M[2] = 2 * H[2] / w_s;
251: D[0] = 0.1 * M[0];
252: D[1] = 0.1 * M[1];
253: D[2] = 0.1 * M[2];
255: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
257: /* Network subsystem initialization */
258: VecCopy(user->V0, Xnet);
260: /* Generator subsystem initialization */
261: VecGetArray(Xgen, &xgen);
262: VecGetArray(Xnet, &xnet);
264: for (i = 0; i < ngen; i++) {
265: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
266: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
267: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
268: Vm2 = Vm * Vm;
269: IGr = (Vr * PGv[i] + Vi * QG[i]) / Vm2;
270: IGi = (Vi * PGv[i] - Vr * QG[i]) / Vm2;
272: delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
274: theta = PETSC_PI / 2.0 - delta;
276: Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
277: Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
279: Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
280: Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
282: Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
283: Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
285: /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
286: xgen[idx] = Eqp;
287: xgen[idx + 1] = Edp;
288: xgen[idx + 2] = delta;
289: xgen[idx + 3] = w_s;
291: idx = idx + 4;
293: xgen[idx] = Id;
294: xgen[idx + 1] = Iq;
296: idx = idx + 2;
298: /* Exciter */
299: Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
300: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
301: VR = KE[i] * Efd + SE;
302: RF = KF[i] * Efd / TF[i];
304: xgen[idx] = Efd;
305: xgen[idx + 1] = RF;
306: xgen[idx + 2] = VR;
308: idx = idx + 3;
309: }
311: VecRestoreArray(Xgen, &xgen);
312: VecRestoreArray(Xnet, &xnet);
314: /* VecView(Xgen,0); */
315: DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet);
316: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
317: return 0;
318: }
320: PetscErrorCode DICDPFiniteDifference(Vec X, Vec *DICDP, Userctx *user)
321: {
322: Vec Y;
323: PetscScalar PGv[3], eps;
324: PetscInt i, j;
326: eps = 1.e-7;
327: VecDuplicate(X, &Y);
329: for (i = 0; i < ngen; i++) {
330: for (j = 0; j < 3; j++) PGv[j] = PG[j];
331: PGv[i] = PG[i] + eps;
332: InitialGuess(Y, user, PGv);
333: InitialGuess(X, user, PG);
335: VecAXPY(Y, -1.0, X);
336: VecScale(Y, 1. / eps);
337: VecCopy(Y, DICDP[i]);
338: }
339: VecDestroy(&Y);
340: return 0;
341: }
343: /* Computes F = [-f(x,y);g(x,y)] */
344: PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
345: {
346: Vec Xgen, Xnet, Fgen, Fnet;
347: PetscScalar *xgen, *xnet, *fgen, *fnet;
348: PetscInt i, idx = 0;
349: PetscScalar Vr, Vi, Vm, Vm2;
350: PetscScalar Eqp, Edp, delta, w; /* Generator variables */
351: PetscScalar Efd, RF, VR; /* Exciter variables */
352: PetscScalar Id, Iq; /* Generator dq axis currents */
353: PetscScalar Vd, Vq, SE;
354: PetscScalar IGr, IGi, IDr, IDi;
355: PetscScalar Zdq_inv[4], det;
356: PetscScalar PD, QD, Vm0, *v0;
357: PetscInt k;
359: VecZeroEntries(F);
360: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
361: DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet);
362: DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);
363: DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet);
365: /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
366: The generator current injection, IG, and load current injection, ID are added later
367: */
368: /* Note that the values in Ybus are stored assuming the imaginary current balance
369: equation is ordered first followed by real current balance equation for each bus.
370: Thus imaginary current contribution goes in location 2*i, and
371: real current contribution in 2*i+1
372: */
373: MatMult(user->Ybus, Xnet, Fnet);
375: VecGetArray(Xgen, &xgen);
376: VecGetArray(Xnet, &xnet);
377: VecGetArray(Fgen, &fgen);
378: VecGetArray(Fnet, &fnet);
380: /* Generator subsystem */
381: for (i = 0; i < ngen; i++) {
382: Eqp = xgen[idx];
383: Edp = xgen[idx + 1];
384: delta = xgen[idx + 2];
385: w = xgen[idx + 3];
386: Id = xgen[idx + 4];
387: Iq = xgen[idx + 5];
388: Efd = xgen[idx + 6];
389: RF = xgen[idx + 7];
390: VR = xgen[idx + 8];
392: /* Generator differential equations */
393: fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
394: fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
395: fgen[idx + 2] = -w + w_s;
396: fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
398: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
399: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
401: ri2dq(Vr, Vi, delta, &Vd, &Vq);
402: /* Algebraic equations for stator currents */
403: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
405: Zdq_inv[0] = Rs[i] / det;
406: Zdq_inv[1] = Xqp[i] / det;
407: Zdq_inv[2] = -Xdp[i] / det;
408: Zdq_inv[3] = Rs[i] / det;
410: fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
411: fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
413: /* Add generator current injection to network */
414: dq2ri(Id, Iq, delta, &IGr, &IGi);
416: fnet[2 * gbus[i]] -= IGi;
417: fnet[2 * gbus[i] + 1] -= IGr;
419: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
421: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
423: /* Exciter differential equations */
424: fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
425: fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
426: fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
428: idx = idx + 9;
429: }
431: VecGetArray(user->V0, &v0);
432: for (i = 0; i < nload; i++) {
433: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
434: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
435: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
436: Vm2 = Vm * Vm;
437: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
438: PD = QD = 0.0;
439: for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
440: for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
442: /* Load currents */
443: IDr = (PD * Vr + QD * Vi) / Vm2;
444: IDi = (-QD * Vr + PD * Vi) / Vm2;
446: fnet[2 * lbus[i]] += IDi;
447: fnet[2 * lbus[i] + 1] += IDr;
448: }
449: VecRestoreArray(user->V0, &v0);
451: VecRestoreArray(Xgen, &xgen);
452: VecRestoreArray(Xnet, &xnet);
453: VecRestoreArray(Fgen, &fgen);
454: VecRestoreArray(Fnet, &fnet);
456: DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet);
457: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
458: DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet);
459: return 0;
460: }
462: /* \dot{x} - f(x,y)
463: g(x,y) = 0
464: */
465: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
466: {
467: SNES snes;
468: PetscScalar *f;
469: const PetscScalar *xdot;
470: PetscInt i;
472: user->t = t;
474: TSGetSNES(ts, &snes);
475: ResidualFunction(snes, X, F, user);
476: VecGetArray(F, &f);
477: VecGetArrayRead(Xdot, &xdot);
478: for (i = 0; i < ngen; i++) {
479: f[9 * i] += xdot[9 * i];
480: f[9 * i + 1] += xdot[9 * i + 1];
481: f[9 * i + 2] += xdot[9 * i + 2];
482: f[9 * i + 3] += xdot[9 * i + 3];
483: f[9 * i + 6] += xdot[9 * i + 6];
484: f[9 * i + 7] += xdot[9 * i + 7];
485: f[9 * i + 8] += xdot[9 * i + 8];
486: }
487: VecRestoreArray(F, &f);
488: VecRestoreArrayRead(Xdot, &xdot);
489: return 0;
490: }
492: /* This function is used for solving the algebraic system only during fault on and
493: off times. It computes the entire F and then zeros out the part corresponding to
494: differential equations
495: F = [0;g(y)];
496: */
497: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
498: {
499: Userctx *user = (Userctx *)ctx;
500: PetscScalar *f;
501: PetscInt i;
503: ResidualFunction(snes, X, F, user);
504: VecGetArray(F, &f);
505: for (i = 0; i < ngen; i++) {
506: f[9 * i] = 0;
507: f[9 * i + 1] = 0;
508: f[9 * i + 2] = 0;
509: f[9 * i + 3] = 0;
510: f[9 * i + 6] = 0;
511: f[9 * i + 7] = 0;
512: f[9 * i + 8] = 0;
513: }
514: VecRestoreArray(F, &f);
515: return 0;
516: }
518: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
519: {
520: PetscInt *d_nnz;
521: PetscInt i, idx = 0, start = 0;
522: PetscInt ncols;
524: PetscMalloc1(user->neqs_pgrid, &d_nnz);
525: for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
526: /* Generator subsystem */
527: for (i = 0; i < ngen; i++) {
528: d_nnz[idx] += 3;
529: d_nnz[idx + 1] += 2;
530: d_nnz[idx + 2] += 2;
531: d_nnz[idx + 3] += 5;
532: d_nnz[idx + 4] += 6;
533: d_nnz[idx + 5] += 6;
535: d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
536: d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
538: d_nnz[idx + 6] += 2;
539: d_nnz[idx + 7] += 2;
540: d_nnz[idx + 8] += 5;
542: idx = idx + 9;
543: }
545: start = user->neqs_gen;
546: for (i = 0; i < nbus; i++) {
547: MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL);
548: d_nnz[start + 2 * i] += ncols;
549: d_nnz[start + 2 * i + 1] += ncols;
550: MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL);
551: }
553: MatSeqAIJSetPreallocation(J, 0, d_nnz);
554: PetscFree(d_nnz);
555: return 0;
556: }
558: /*
559: J = [-df_dx, -df_dy
560: dg_dx, dg_dy]
561: */
562: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, void *ctx)
563: {
564: Userctx *user = (Userctx *)ctx;
565: Vec Xgen, Xnet;
566: PetscScalar *xgen, *xnet;
567: PetscInt i, idx = 0;
568: PetscScalar Vr, Vi, Vm, Vm2;
569: PetscScalar Eqp, Edp, delta; /* Generator variables */
570: PetscScalar Efd; /* Exciter variables */
571: PetscScalar Id, Iq; /* Generator dq axis currents */
572: PetscScalar Vd, Vq;
573: PetscScalar val[10];
574: PetscInt row[2], col[10];
575: PetscInt net_start = user->neqs_gen;
576: PetscInt ncols;
577: const PetscInt *cols;
578: const PetscScalar *yvals;
579: PetscInt k;
580: PetscScalar Zdq_inv[4], det;
581: PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
582: PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
583: PetscScalar dSE_dEfd;
584: PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
585: PetscScalar PD, QD, Vm0, *v0, Vm4;
586: PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
587: PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
589: MatZeroEntries(B);
590: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
591: DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);
593: VecGetArray(Xgen, &xgen);
594: VecGetArray(Xnet, &xnet);
596: /* Generator subsystem */
597: for (i = 0; i < ngen; i++) {
598: Eqp = xgen[idx];
599: Edp = xgen[idx + 1];
600: delta = xgen[idx + 2];
601: Id = xgen[idx + 4];
602: Iq = xgen[idx + 5];
603: Efd = xgen[idx + 6];
605: /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
606: row[0] = idx;
607: col[0] = idx;
608: col[1] = idx + 4;
609: col[2] = idx + 6;
610: val[0] = 1 / Td0p[i];
611: val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
612: val[2] = -1 / Td0p[i];
614: MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);
616: /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
617: row[0] = idx + 1;
618: col[0] = idx + 1;
619: col[1] = idx + 5;
620: val[0] = 1 / Tq0p[i];
621: val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
622: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
624: /* fgen[idx+2] = - w + w_s; */
625: row[0] = idx + 2;
626: col[0] = idx + 2;
627: col[1] = idx + 3;
628: val[0] = 0;
629: val[1] = -1;
630: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
632: /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
633: row[0] = idx + 3;
634: col[0] = idx;
635: col[1] = idx + 1;
636: col[2] = idx + 3;
637: col[3] = idx + 4;
638: col[4] = idx + 5;
639: val[0] = Iq / M[i];
640: val[1] = Id / M[i];
641: val[2] = D[i] / M[i];
642: val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
643: val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
644: MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES);
646: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
647: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
648: ri2dq(Vr, Vi, delta, &Vd, &Vq);
650: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
652: Zdq_inv[0] = Rs[i] / det;
653: Zdq_inv[1] = Xqp[i] / det;
654: Zdq_inv[2] = -Xdp[i] / det;
655: Zdq_inv[3] = Rs[i] / det;
657: dVd_dVr = PetscSinScalar(delta);
658: dVd_dVi = -PetscCosScalar(delta);
659: dVq_dVr = PetscCosScalar(delta);
660: dVq_dVi = PetscSinScalar(delta);
661: dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
662: dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
664: /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
665: row[0] = idx + 4;
666: col[0] = idx;
667: col[1] = idx + 1;
668: col[2] = idx + 2;
669: val[0] = -Zdq_inv[1];
670: val[1] = -Zdq_inv[0];
671: val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
672: col[3] = idx + 4;
673: col[4] = net_start + 2 * gbus[i];
674: col[5] = net_start + 2 * gbus[i] + 1;
675: val[3] = 1;
676: val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
677: val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
678: MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES);
680: /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
681: row[0] = idx + 5;
682: col[0] = idx;
683: col[1] = idx + 1;
684: col[2] = idx + 2;
685: val[0] = -Zdq_inv[3];
686: val[1] = -Zdq_inv[2];
687: val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
688: col[3] = idx + 5;
689: col[4] = net_start + 2 * gbus[i];
690: col[5] = net_start + 2 * gbus[i] + 1;
691: val[3] = 1;
692: val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
693: val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
694: MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES);
696: dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
697: dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
698: dIGr_dId = PetscSinScalar(delta);
699: dIGr_dIq = PetscCosScalar(delta);
700: dIGi_dId = -PetscCosScalar(delta);
701: dIGi_dIq = PetscSinScalar(delta);
703: /* fnet[2*gbus[i]] -= IGi; */
704: row[0] = net_start + 2 * gbus[i];
705: col[0] = idx + 2;
706: col[1] = idx + 4;
707: col[2] = idx + 5;
708: val[0] = -dIGi_ddelta;
709: val[1] = -dIGi_dId;
710: val[2] = -dIGi_dIq;
711: MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);
713: /* fnet[2*gbus[i]+1] -= IGr; */
714: row[0] = net_start + 2 * gbus[i] + 1;
715: col[0] = idx + 2;
716: col[1] = idx + 4;
717: col[2] = idx + 5;
718: val[0] = -dIGr_ddelta;
719: val[1] = -dIGr_dId;
720: val[2] = -dIGr_dIq;
721: MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);
723: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
725: /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
726: /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
727: dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
729: row[0] = idx + 6;
730: col[0] = idx + 6;
731: col[1] = idx + 8;
732: val[0] = (KE[i] + dSE_dEfd) / TE[i];
733: val[1] = -1 / TE[i];
734: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
736: /* Exciter differential equations */
738: /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
739: row[0] = idx + 7;
740: col[0] = idx + 6;
741: col[1] = idx + 7;
742: val[0] = (-KF[i] / TF[i]) / TF[i];
743: val[1] = 1 / TF[i];
744: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
746: /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
747: /* Vm = (Vd^2 + Vq^2)^0.5; */
748: dVm_dVd = Vd / Vm;
749: dVm_dVq = Vq / Vm;
750: dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
751: dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
752: row[0] = idx + 8;
753: col[0] = idx + 6;
754: col[1] = idx + 7;
755: col[2] = idx + 8;
756: val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
757: val[1] = -KA[i] / TA[i];
758: val[2] = 1 / TA[i];
759: col[3] = net_start + 2 * gbus[i];
760: col[4] = net_start + 2 * gbus[i] + 1;
761: val[3] = KA[i] * dVm_dVr / TA[i];
762: val[4] = KA[i] * dVm_dVi / TA[i];
763: MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES);
764: idx = idx + 9;
765: }
767: for (i = 0; i < nbus; i++) {
768: MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals);
769: row[0] = net_start + 2 * i;
770: for (k = 0; k < ncols; k++) {
771: col[k] = net_start + cols[k];
772: val[k] = yvals[k];
773: }
774: MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES);
775: MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals);
777: MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals);
778: row[0] = net_start + 2 * i + 1;
779: for (k = 0; k < ncols; k++) {
780: col[k] = net_start + cols[k];
781: val[k] = yvals[k];
782: }
783: MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES);
784: MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals);
785: }
787: MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY);
788: MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY);
790: VecGetArray(user->V0, &v0);
791: for (i = 0; i < nload; i++) {
792: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
793: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
794: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
795: Vm2 = Vm * Vm;
796: Vm4 = Vm2 * Vm2;
797: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
798: PD = QD = 0.0;
799: dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
800: for (k = 0; k < ld_nsegsp[i]; k++) {
801: PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
802: dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vr * PetscPowScalar(Vm, (ld_betap[k] - 2));
803: dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vi * PetscPowScalar(Vm, (ld_betap[k] - 2));
804: }
805: for (k = 0; k < ld_nsegsq[i]; k++) {
806: QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
807: dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vr * PetscPowScalar(Vm, (ld_betaq[k] - 2));
808: dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vi * PetscPowScalar(Vm, (ld_betaq[k] - 2));
809: }
811: /* IDr = (PD*Vr + QD*Vi)/Vm2; */
812: /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
814: dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
815: dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
817: dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
818: dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
820: /* fnet[2*lbus[i]] += IDi; */
821: row[0] = net_start + 2 * lbus[i];
822: col[0] = net_start + 2 * lbus[i];
823: col[1] = net_start + 2 * lbus[i] + 1;
824: val[0] = dIDi_dVr;
825: val[1] = dIDi_dVi;
826: MatSetValues(J, 1, row, 2, col, val, ADD_VALUES);
827: /* fnet[2*lbus[i]+1] += IDr; */
828: row[0] = net_start + 2 * lbus[i] + 1;
829: col[0] = net_start + 2 * lbus[i];
830: col[1] = net_start + 2 * lbus[i] + 1;
831: val[0] = dIDr_dVr;
832: val[1] = dIDr_dVi;
833: MatSetValues(J, 1, row, 2, col, val, ADD_VALUES);
834: }
835: VecRestoreArray(user->V0, &v0);
837: VecRestoreArray(Xgen, &xgen);
838: VecRestoreArray(Xnet, &xnet);
840: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
842: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
843: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
844: return 0;
845: }
847: /*
848: J = [I, 0
849: dg_dx, dg_dy]
850: */
851: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, void *ctx)
852: {
853: Userctx *user = (Userctx *)ctx;
855: ResidualJacobian(snes, X, A, B, ctx);
856: MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
857: MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL);
858: return 0;
859: }
861: /*
862: J = [a*I-df_dx, -df_dy
863: dg_dx, dg_dy]
864: */
866: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
867: {
868: SNES snes;
869: PetscScalar atmp = (PetscScalar)a;
870: PetscInt i, row;
872: user->t = t;
874: TSGetSNES(ts, &snes);
875: ResidualJacobian(snes, X, A, B, user);
876: for (i = 0; i < ngen; i++) {
877: row = 9 * i;
878: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
879: row = 9 * i + 1;
880: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
881: row = 9 * i + 2;
882: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
883: row = 9 * i + 3;
884: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
885: row = 9 * i + 6;
886: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
887: row = 9 * i + 7;
888: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
889: row = 9 * i + 8;
890: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
891: }
892: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
893: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
894: return 0;
895: }
897: /* Matrix JacobianP is constant so that it only needs to be evaluated once */
898: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
899: {
900: PetscScalar a;
901: PetscInt row, col;
902: Userctx *ctx = (Userctx *)ctx0;
906: if (ctx->jacp_flg) {
907: MatZeroEntries(A);
909: for (col = 0; col < 3; col++) {
910: a = 1.0 / M[col];
911: row = 9 * col + 3;
912: MatSetValues(A, 1, &row, 1, &col, &a, INSERT_VALUES);
913: }
915: ctx->jacp_flg = PETSC_FALSE;
917: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
918: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
919: }
920: return 0;
921: }
923: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
924: {
925: const PetscScalar *u;
926: PetscInt idx;
927: Vec Xgen, Xnet;
928: PetscScalar *r, *xgen;
929: PetscInt i;
931: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
932: DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet);
934: VecGetArray(Xgen, &xgen);
936: VecGetArrayRead(U, &u);
937: VecGetArray(R, &r);
938: r[0] = 0.;
939: idx = 0;
940: for (i = 0; i < ngen; i++) {
941: 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);
942: idx += 9;
943: }
944: VecRestoreArrayRead(U, &u);
945: VecRestoreArray(R, &r);
946: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
947: return 0;
948: }
950: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, Userctx *user)
951: {
952: Vec Xgen, Xnet, Dgen, Dnet;
953: PetscScalar *xgen, *dgen;
954: PetscInt i;
955: PetscInt idx;
956: Vec drdu_col;
957: PetscScalar *xarr;
959: VecDuplicate(U, &drdu_col);
960: MatDenseGetColumn(DRDU, 0, &xarr);
961: VecPlaceArray(drdu_col, xarr);
962: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
963: DMCompositeGetLocalVectors(user->dmpgrid, &Dgen, &Dnet);
964: DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet);
965: DMCompositeScatter(user->dmpgrid, drdu_col, Dgen, Dnet);
967: VecGetArray(Xgen, &xgen);
968: VecGetArray(Dgen, &dgen);
970: idx = 0;
971: for (i = 0; i < ngen; i++) {
972: dgen[idx + 3] = 0.;
973: 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);
974: 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);
975: idx += 9;
976: }
978: VecRestoreArray(Dgen, &dgen);
979: VecRestoreArray(Xgen, &xgen);
980: DMCompositeGather(user->dmpgrid, INSERT_VALUES, drdu_col, Dgen, Dnet);
981: DMCompositeRestoreLocalVectors(user->dmpgrid, &Dgen, &Dnet);
982: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
983: VecResetArray(drdu_col);
984: MatDenseRestoreColumn(DRDU, &xarr);
985: VecDestroy(&drdu_col);
986: return 0;
987: }
989: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat drdp, Userctx *user)
990: {
991: return 0;
992: }
994: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, Vec *DICDP, Userctx *user)
995: {
996: PetscScalar *x, *y, sensip;
997: PetscInt i;
999: VecGetArray(lambda, &x);
1000: VecGetArray(mu, &y);
1002: for (i = 0; i < 3; i++) {
1003: VecDot(lambda, DICDP[i], &sensip);
1004: sensip = sensip + y[i];
1005: /* PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt %" PetscInt_FMT " th parameter: %g \n",i,(double)sensip); */
1006: y[i] = sensip;
1007: }
1008: VecRestoreArray(mu, &y);
1009: return 0;
1010: }
1012: int main(int argc, char **argv)
1013: {
1014: Userctx user;
1015: Vec p;
1016: PetscScalar *x_ptr;
1017: PetscMPIInt size;
1018: PetscInt i;
1019: PetscViewer Xview, Ybusview;
1020: PetscInt *idx2;
1021: Tao tao;
1022: KSP ksp;
1023: PC pc;
1024: Vec lowerb, upperb;
1027: PetscInitialize(&argc, &argv, "petscoptions", help);
1028: MPI_Comm_size(PETSC_COMM_WORLD, &size);
1031: user.jacp_flg = PETSC_TRUE;
1032: user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
1033: user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
1034: user.neqs_pgrid = user.neqs_gen + user.neqs_net;
1036: /* Create indices for differential and algebraic equations */
1037: PetscMalloc1(7 * ngen, &idx2);
1038: for (i = 0; i < ngen; i++) {
1039: idx2[7 * i] = 9 * i;
1040: idx2[7 * i + 1] = 9 * i + 1;
1041: idx2[7 * i + 2] = 9 * i + 2;
1042: idx2[7 * i + 3] = 9 * i + 3;
1043: idx2[7 * i + 4] = 9 * i + 6;
1044: idx2[7 * i + 5] = 9 * i + 7;
1045: idx2[7 * i + 6] = 9 * i + 8;
1046: }
1047: ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff);
1048: ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg);
1049: PetscFree(idx2);
1051: /* Set run time options */
1052: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1053: {
1054: user.tfaulton = 1.0;
1055: user.tfaultoff = 1.2;
1056: user.Rfault = 0.0001;
1057: user.faultbus = 8;
1058: PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL);
1059: PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL);
1060: PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL);
1061: user.t0 = 0.0;
1062: user.tmax = 1.3;
1063: PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL);
1064: PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL);
1065: user.freq_u = 61.0;
1066: user.freq_l = 59.0;
1067: user.pow = 2;
1068: PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL);
1069: PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL);
1070: PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL);
1071: }
1072: PetscOptionsEnd();
1074: /* Create DMs for generator and network subsystems */
1075: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen);
1076: DMSetOptionsPrefix(user.dmgen, "dmgen_");
1077: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet);
1078: DMSetOptionsPrefix(user.dmnet, "dmnet_");
1079: DMSetFromOptions(user.dmnet);
1080: DMSetUp(user.dmnet);
1081: /* Create a composite DM packer and add the two DMs */
1082: DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid);
1083: DMSetOptionsPrefix(user.dmpgrid, "pgrid_");
1084: DMSetFromOptions(user.dmgen);
1085: DMSetUp(user.dmgen);
1086: DMCompositeAddDM(user.dmpgrid, user.dmgen);
1087: DMCompositeAddDM(user.dmpgrid, user.dmnet);
1089: /* Read initial voltage vector and Ybus */
1090: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview);
1091: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview);
1093: VecCreate(PETSC_COMM_WORLD, &user.V0);
1094: VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net);
1095: VecLoad(user.V0, Xview);
1097: MatCreate(PETSC_COMM_WORLD, &user.Ybus);
1098: MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net);
1099: MatSetType(user.Ybus, MATBAIJ);
1100: /* MatSetBlockSize(ctx->Ybus,2); */
1101: MatLoad(user.Ybus, Ybusview);
1103: PetscViewerDestroy(&Xview);
1104: PetscViewerDestroy(&Ybusview);
1106: /* Allocate space for Jacobians */
1107: MatCreate(PETSC_COMM_WORLD, &user.J);
1108: MatSetSizes(user.J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid);
1109: MatSetFromOptions(user.J);
1110: PreallocateJacobian(user.J, &user);
1112: MatCreate(PETSC_COMM_WORLD, &user.Jacp);
1113: MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 3);
1114: MatSetFromOptions(user.Jacp);
1115: MatSetUp(user.Jacp);
1116: MatZeroEntries(user.Jacp); /* initialize to zeros */
1118: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, &user.DRDP);
1119: MatSetUp(user.DRDP);
1120: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 1, NULL, &user.DRDU);
1121: MatSetUp(user.DRDU);
1123: /* Create TAO solver and set desired solution method */
1124: TaoCreate(PETSC_COMM_WORLD, &tao);
1125: TaoSetType(tao, TAOBLMVM);
1126: /*
1127: Optimization starts
1128: */
1129: /* Set initial solution guess */
1130: VecCreateSeq(PETSC_COMM_WORLD, 3, &p);
1131: VecGetArray(p, &x_ptr);
1132: x_ptr[0] = PG[0];
1133: x_ptr[1] = PG[1];
1134: x_ptr[2] = PG[2];
1135: VecRestoreArray(p, &x_ptr);
1137: TaoSetSolution(tao, p);
1138: /* Set routine for function and gradient evaluation */
1139: TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user);
1141: /* Set bounds for the optimization */
1142: VecDuplicate(p, &lowerb);
1143: VecDuplicate(p, &upperb);
1144: VecGetArray(lowerb, &x_ptr);
1145: x_ptr[0] = 0.5;
1146: x_ptr[1] = 0.5;
1147: x_ptr[2] = 0.5;
1148: VecRestoreArray(lowerb, &x_ptr);
1149: VecGetArray(upperb, &x_ptr);
1150: x_ptr[0] = 2.0;
1151: x_ptr[1] = 2.0;
1152: x_ptr[2] = 2.0;
1153: VecRestoreArray(upperb, &x_ptr);
1154: TaoSetVariableBounds(tao, lowerb, upperb);
1156: /* Check for any TAO command line options */
1157: TaoSetFromOptions(tao);
1158: TaoGetKSP(tao, &ksp);
1159: if (ksp) {
1160: KSPGetPC(ksp, &pc);
1161: PCSetType(pc, PCNONE);
1162: }
1164: /* SOLVE THE APPLICATION */
1165: TaoSolve(tao);
1167: VecView(p, PETSC_VIEWER_STDOUT_WORLD);
1168: /* Free TAO data structures */
1169: TaoDestroy(&tao);
1171: DMDestroy(&user.dmgen);
1172: DMDestroy(&user.dmnet);
1173: DMDestroy(&user.dmpgrid);
1174: ISDestroy(&user.is_diff);
1175: ISDestroy(&user.is_alg);
1177: MatDestroy(&user.J);
1178: MatDestroy(&user.Jacp);
1179: MatDestroy(&user.Ybus);
1180: /* MatDestroy(&user.Sol); */
1181: VecDestroy(&user.V0);
1182: VecDestroy(&p);
1183: VecDestroy(&lowerb);
1184: VecDestroy(&upperb);
1185: MatDestroy(&user.DRDU);
1186: MatDestroy(&user.DRDP);
1187: PetscFinalize();
1188: return 0;
1189: }
1191: /* ------------------------------------------------------------------ */
1192: /*
1193: FormFunction - Evaluates the function and corresponding gradient.
1195: Input Parameters:
1196: tao - the Tao context
1197: X - the input vector
1198: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
1200: Output Parameters:
1201: f - the newly evaluated function
1202: G - the newly evaluated gradient
1203: */
1204: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx0)
1205: {
1206: TS ts, quadts;
1207: SNES snes_alg;
1208: Userctx *ctx = (Userctx *)ctx0;
1209: Vec X;
1210: PetscInt i;
1211: /* sensitivity context */
1212: PetscScalar *x_ptr;
1213: Vec lambda[1], q;
1214: Vec mu[1];
1215: PetscInt steps1, steps2, steps3;
1216: Vec DICDP[3];
1217: Vec F_alg;
1218: PetscInt row_loc, col_loc;
1219: PetscScalar val;
1220: Vec Xdot;
1222: VecGetArrayRead(P, (const PetscScalar **)&x_ptr);
1223: PG[0] = x_ptr[0];
1224: PG[1] = x_ptr[1];
1225: PG[2] = x_ptr[2];
1226: VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr);
1228: ctx->stepnum = 0;
1230: DMCreateGlobalVector(ctx->dmpgrid, &X);
1232: /* Create matrix to save solutions at each time step */
1233: /* MatCreateSeqDense(PETSC_COMM_SELF,ctx->neqs_pgrid+1,1002,NULL,&ctx->Sol); */
1234: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1235: Create timestepping solver context
1236: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1237: TSCreate(PETSC_COMM_WORLD, &ts);
1238: TSSetProblemType(ts, TS_NONLINEAR);
1239: TSSetType(ts, TSCN);
1240: TSSetIFunction(ts, NULL, (TSIFunction)IFunction, ctx);
1241: TSSetIJacobian(ts, ctx->J, ctx->J, (TSIJacobian)IJacobian, ctx);
1242: TSSetApplicationContext(ts, ctx);
1243: /* Set RHS JacobianP */
1244: TSSetRHSJacobianP(ts, ctx->Jacp, RHSJacobianP, ctx);
1246: TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts);
1247: TSSetRHSFunction(quadts, NULL, (TSRHSFunction)CostIntegrand, ctx);
1248: TSSetRHSJacobian(quadts, ctx->DRDU, ctx->DRDU, (TSRHSJacobian)DRDUJacobianTranspose, ctx);
1249: TSSetRHSJacobianP(quadts, ctx->DRDP, (TSRHSJacobianP)DRDPJacobianTranspose, ctx);
1251: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1252: Set initial conditions
1253: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1254: SetInitialGuess(X, ctx);
1256: /* Approximate DICDP with finite difference, we want to zero out network variables */
1257: for (i = 0; i < 3; i++) VecDuplicate(X, &DICDP[i]);
1258: DICDPFiniteDifference(X, DICDP, ctx);
1260: VecDuplicate(X, &F_alg);
1261: SNESCreate(PETSC_COMM_WORLD, &snes_alg);
1262: SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx);
1263: MatZeroEntries(ctx->J);
1264: SNESSetJacobian(snes_alg, ctx->J, ctx->J, AlgJacobian, ctx);
1265: SNESSetOptionsPrefix(snes_alg, "alg_");
1266: SNESSetFromOptions(snes_alg);
1267: ctx->alg_flg = PETSC_TRUE;
1268: /* Solve the algebraic equations */
1269: SNESSolve(snes_alg, NULL, X);
1271: /* Just to set up the Jacobian structure */
1272: VecDuplicate(X, &Xdot);
1273: IJacobian(ts, 0.0, X, Xdot, 0.0, ctx->J, ctx->J, ctx);
1274: VecDestroy(&Xdot);
1276: ctx->stepnum++;
1278: /*
1279: Save trajectory of solution so that TSAdjointSolve() may be used
1280: */
1281: TSSetSaveTrajectory(ts);
1283: TSSetTimeStep(ts, 0.01);
1284: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1285: TSSetFromOptions(ts);
1286: /* TSSetPostStep(ts,SaveSolution); */
1288: /* Prefault period */
1289: ctx->alg_flg = PETSC_FALSE;
1290: TSSetTime(ts, 0.0);
1291: TSSetMaxTime(ts, ctx->tfaulton);
1292: TSSolve(ts, X);
1293: TSGetStepNumber(ts, &steps1);
1295: /* Create the nonlinear solver for solving the algebraic system */
1296: /* Note that although the algebraic system needs to be solved only for
1297: Idq and V, we reuse the entire system including xgen. The xgen
1298: variables are held constant by setting their residuals to 0 and
1299: putting a 1 on the Jacobian diagonal for xgen rows
1300: */
1301: MatZeroEntries(ctx->J);
1303: /* Apply disturbance - resistive fault at ctx->faultbus */
1304: /* This is done by adding shunt conductance to the diagonal location
1305: in the Ybus matrix */
1306: row_loc = 2 * ctx->faultbus;
1307: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1308: val = 1 / ctx->Rfault;
1309: MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1310: row_loc = 2 * ctx->faultbus + 1;
1311: col_loc = 2 * ctx->faultbus; /* Location for G */
1312: val = 1 / ctx->Rfault;
1313: MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1315: MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1316: MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1318: ctx->alg_flg = PETSC_TRUE;
1319: /* Solve the algebraic equations */
1320: SNESSolve(snes_alg, NULL, X);
1322: ctx->stepnum++;
1324: /* Disturbance period */
1325: ctx->alg_flg = PETSC_FALSE;
1326: TSSetTime(ts, ctx->tfaulton);
1327: TSSetMaxTime(ts, ctx->tfaultoff);
1328: TSSolve(ts, X);
1329: TSGetStepNumber(ts, &steps2);
1330: steps2 -= steps1;
1332: /* Remove the fault */
1333: row_loc = 2 * ctx->faultbus;
1334: col_loc = 2 * ctx->faultbus + 1;
1335: val = -1 / ctx->Rfault;
1336: MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1337: row_loc = 2 * ctx->faultbus + 1;
1338: col_loc = 2 * ctx->faultbus;
1339: val = -1 / ctx->Rfault;
1340: MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1342: MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1343: MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1345: MatZeroEntries(ctx->J);
1347: ctx->alg_flg = PETSC_TRUE;
1349: /* Solve the algebraic equations */
1350: SNESSolve(snes_alg, NULL, X);
1352: ctx->stepnum++;
1354: /* Post-disturbance period */
1355: ctx->alg_flg = PETSC_TRUE;
1356: TSSetTime(ts, ctx->tfaultoff);
1357: TSSetMaxTime(ts, ctx->tmax);
1358: TSSolve(ts, X);
1359: TSGetStepNumber(ts, &steps3);
1360: steps3 -= steps2;
1361: steps3 -= steps1;
1363: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1364: Adjoint model starts here
1365: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1366: TSSetPostStep(ts, NULL);
1367: MatCreateVecs(ctx->J, &lambda[0], NULL);
1368: /* Set initial conditions for the adjoint integration */
1369: VecZeroEntries(lambda[0]);
1371: MatCreateVecs(ctx->Jacp, &mu[0], NULL);
1372: VecZeroEntries(mu[0]);
1373: TSSetCostGradients(ts, 1, lambda, mu);
1375: TSAdjointSetSteps(ts, steps3);
1376: TSAdjointSolve(ts);
1378: MatZeroEntries(ctx->J);
1379: /* Applying disturbance - resistive fault at ctx->faultbus */
1380: /* This is done by deducting shunt conductance to the diagonal location
1381: in the Ybus matrix */
1382: row_loc = 2 * ctx->faultbus;
1383: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1384: val = 1. / ctx->Rfault;
1385: MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1386: row_loc = 2 * ctx->faultbus + 1;
1387: col_loc = 2 * ctx->faultbus; /* Location for G */
1388: val = 1. / ctx->Rfault;
1389: MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1391: MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1392: MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1394: /* Set number of steps for the adjoint integration */
1395: TSAdjointSetSteps(ts, steps2);
1396: TSAdjointSolve(ts);
1398: MatZeroEntries(ctx->J);
1399: /* remove the fault */
1400: row_loc = 2 * ctx->faultbus;
1401: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1402: val = -1. / ctx->Rfault;
1403: 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: MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1409: MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1410: MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1412: /* Set number of steps for the adjoint integration */
1413: TSAdjointSetSteps(ts, steps1);
1414: TSAdjointSolve(ts);
1416: ComputeSensiP(lambda[0], mu[0], DICDP, ctx);
1417: VecCopy(mu[0], G);
1419: TSGetQuadratureTS(ts, NULL, &quadts);
1420: TSGetSolution(quadts, &q);
1421: VecGetArray(q, &x_ptr);
1422: *f = x_ptr[0];
1423: x_ptr[0] = 0;
1424: VecRestoreArray(q, &x_ptr);
1426: VecDestroy(&lambda[0]);
1427: VecDestroy(&mu[0]);
1429: SNESDestroy(&snes_alg);
1430: VecDestroy(&F_alg);
1431: VecDestroy(&X);
1432: TSDestroy(&ts);
1433: for (i = 0; i < 3; i++) VecDestroy(&DICDP[i]);
1434: return 0;
1435: }
1437: /*TEST
1439: build:
1440: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1442: test:
1443: args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1444: localrunfiles: petscoptions X.bin Ybus.bin
1446: TEST*/