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*/