Actual source code: ex9busopt.c

  1: static char help[] = "Application of adjoint sensitivity analysis for power grid stability analysis of WECC 9 bus system.\n\
  2: This example is based on the 9-bus (node) example given in the book Power\n\
  3: Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\
  4: The power grid in this example consists of 9 buses (nodes), 3 generators,\n\
  5: 3 loads, and 9 transmission lines. The network equations are written\n\
  6: in current balance form using rectangular coordinates.\n\n";

  8: /*
  9:   This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
 10:   The objectivie is to find optimal parameter PG for each generator to minizie the frequency violations due to faults.
 11:   The problem features discontinuities and a cost function in integral form.
 12:   The gradient is computed with the discrete adjoint of an implicit theta method, see ex9busadj.c for details.
 13: */

 15: #include <petsctao.h>
 16: #include <petscts.h>
 17: #include <petscdm.h>
 18: #include <petscdmda.h>
 19: #include <petscdmcomposite.h>
 20: #include <petsctime.h>

 22: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);

 24: #define freq 60
 25: #define w_s  (2 * PETSC_PI * freq)

 27: /* Sizes and indices */
 28: const PetscInt nbus    = 9;         /* Number of network buses */
 29: const PetscInt ngen    = 3;         /* Number of generators */
 30: const PetscInt nload   = 3;         /* Number of loads */
 31: const PetscInt gbus[3] = {0, 1, 2}; /* Buses at which generators are incident */
 32: const PetscInt lbus[3] = {4, 5, 7}; /* Buses at which loads are incident */

 34: /* Generator real and reactive powers (found via loadflow) */
 35: PetscScalar PG[3] = {0.69, 1.59, 0.69};
 36: /* PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};*/

 38: const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
 39: /* Generator constants */
 40: const PetscScalar H[3]    = {23.64, 6.4, 3.01};       /* Inertia constant */
 41: const PetscScalar Rs[3]   = {0.0, 0.0, 0.0};          /* Stator Resistance */
 42: const PetscScalar Xd[3]   = {0.146, 0.8958, 1.3125};  /* d-axis reactance */
 43: const PetscScalar Xdp[3]  = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
 44: const PetscScalar Xq[3]   = {0.4360, 0.8645, 1.2578}; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
 45: const PetscScalar Xqp[3]  = {0.0969, 0.1969, 0.25};   /* q-axis transient reactance */
 46: const PetscScalar Td0p[3] = {8.96, 6.0, 5.89};        /* d-axis open circuit time constant */
 47: const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6};       /* q-axis open circuit time constant */
 48: PetscScalar       M[3];                               /* M = 2*H/w_s */
 49: PetscScalar       D[3];                               /* D = 0.1*M */

 51: PetscScalar TM[3]; /* Mechanical Torque */
 52: /* Exciter system constants */
 53: const PetscScalar KA[3] = {20.0, 20.0, 20.0};    /* Voltage regulartor gain constant */
 54: const PetscScalar TA[3] = {0.2, 0.2, 0.2};       /* Voltage regulator time constant */
 55: const PetscScalar KE[3] = {1.0, 1.0, 1.0};       /* Exciter gain constant */
 56: const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
 57: const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
 58: const PetscScalar TF[3] = {0.35, 0.35, 0.35};    /* Feedback stabilizer time constant */
 59: const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
 60: const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */

 62: PetscScalar Vref[3];
 63: /* Load constants
 64:   We use a composite load model that describes the load and reactive powers at each time instant as follows
 65:   P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
 66:   Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
 67:   where
 68:     ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
 69:     ld_alphap,ld_alphap - Percentage contribution (weights) or loads
 70:     P_D0                - Real power load
 71:     Q_D0                - Reactive power load
 72:     V_m(t)              - Voltage magnitude at time t
 73:     V_m0                - Voltage magnitude at t = 0
 74:     ld_betap, ld_betaq  - exponents describing the load model for real and reactive part

 76:     Note: All loads have the same characteristic currently.
 77: */
 78: const PetscScalar PD0[3]       = {1.25, 0.9, 1.0};
 79: const PetscScalar QD0[3]       = {0.5, 0.3, 0.35};
 80: const PetscInt    ld_nsegsp[3] = {3, 3, 3};
 81: const PetscScalar ld_alphap[3] = {1.0, 0.0, 0.0};
 82: const PetscScalar ld_betap[3]  = {2.0, 1.0, 0.0};
 83: const PetscInt    ld_nsegsq[3] = {3, 3, 3};
 84: const PetscScalar ld_alphaq[3] = {1.0, 0.0, 0.0};
 85: const PetscScalar ld_betaq[3]  = {2.0, 1.0, 0.0};

 87: typedef struct {
 88:   DM          dmgen, dmnet;        /* DMs to manage generator and network subsystem */
 89:   DM          dmpgrid;             /* Composite DM to manage the entire power grid */
 90:   Mat         Ybus;                /* Network admittance matrix */
 91:   Vec         V0;                  /* Initial voltage vector (Power flow solution) */
 92:   PetscReal   tfaulton, tfaultoff; /* Fault on and off times */
 93:   PetscInt    faultbus;            /* Fault bus */
 94:   PetscScalar Rfault;
 95:   PetscReal   t0, tmax;
 96:   PetscInt    neqs_gen, neqs_net, neqs_pgrid;
 97:   Mat         Sol; /* Matrix to save solution at each time step */
 98:   PetscInt    stepnum;
 99:   PetscBool   alg_flg;
100:   PetscReal   t;
101:   IS          is_diff;        /* indices for differential equations */
102:   IS          is_alg;         /* indices for algebraic equations */
103:   PetscReal   freq_u, freq_l; /* upper and lower frequency limit */
104:   PetscInt    pow;            /* power coefficient used in the cost function */
105:   PetscBool   jacp_flg;
106:   Mat         J, Jacp;
107:   Mat         DRDU, DRDP;
108: } Userctx;

110: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
111: PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
112: {
113:   PetscFunctionBegin;
114:   *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
115:   *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
116:   PetscFunctionReturn(PETSC_SUCCESS);
117: }

119: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
120: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
121: {
122:   PetscFunctionBegin;
123:   *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
124:   *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: /* Saves the solution at each time to a matrix */
129: PetscErrorCode SaveSolution(TS ts)
130: {
131:   Userctx           *user;
132:   Vec                X;
133:   PetscScalar       *mat;
134:   const PetscScalar *x;
135:   PetscInt           idx;
136:   PetscReal          t;

138:   PetscFunctionBegin;
139:   PetscCall(TSGetApplicationContext(ts, &user));
140:   PetscCall(TSGetTime(ts, &t));
141:   PetscCall(TSGetSolution(ts, &X));
142:   idx = user->stepnum * (user->neqs_pgrid + 1);
143:   PetscCall(MatDenseGetArray(user->Sol, &mat));
144:   PetscCall(VecGetArrayRead(X, &x));
145:   mat[idx] = t;
146:   PetscCall(PetscArraycpy(mat + idx + 1, x, user->neqs_pgrid));
147:   PetscCall(MatDenseRestoreArray(user->Sol, &mat));
148:   PetscCall(VecRestoreArrayRead(X, &x));
149:   user->stepnum++;
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
154: {
155:   Vec          Xgen, Xnet;
156:   PetscScalar *xgen, *xnet;
157:   PetscInt     i, idx = 0;
158:   PetscScalar  Vr, Vi, IGr, IGi, Vm, Vm2;
159:   PetscScalar  Eqp, Edp, delta;
160:   PetscScalar  Efd, RF, VR; /* Exciter variables */
161:   PetscScalar  Id, Iq;      /* Generator dq axis currents */
162:   PetscScalar  theta, Vd, Vq, SE;

164:   PetscFunctionBegin;
165:   M[0] = 2 * H[0] / w_s;
166:   M[1] = 2 * H[1] / w_s;
167:   M[2] = 2 * H[2] / w_s;
168:   D[0] = 0.1 * M[0];
169:   D[1] = 0.1 * M[1];
170:   D[2] = 0.1 * M[2];

172:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));

174:   /* Network subsystem initialization */
175:   PetscCall(VecCopy(user->V0, Xnet));

177:   /* Generator subsystem initialization */
178:   PetscCall(VecGetArray(Xgen, &xgen));
179:   PetscCall(VecGetArray(Xnet, &xnet));

181:   for (i = 0; i < ngen; i++) {
182:     Vr  = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
183:     Vi  = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
184:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
185:     Vm2 = Vm * Vm;
186:     IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
187:     IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;

189:     delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */

191:     theta = PETSC_PI / 2.0 - delta;

193:     Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
194:     Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */

196:     Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
197:     Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);

199:     Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
200:     Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */

202:     TM[i] = PG[i];

204:     /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
205:     xgen[idx]     = Eqp;
206:     xgen[idx + 1] = Edp;
207:     xgen[idx + 2] = delta;
208:     xgen[idx + 3] = w_s;

210:     idx = idx + 4;

212:     xgen[idx]     = Id;
213:     xgen[idx + 1] = Iq;

215:     idx = idx + 2;

217:     /* Exciter */
218:     Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
219:     SE  = k1[i] * PetscExpScalar(k2[i] * Efd);
220:     VR  = KE[i] * Efd + SE;
221:     RF  = KF[i] * Efd / TF[i];

223:     xgen[idx]     = Efd;
224:     xgen[idx + 1] = RF;
225:     xgen[idx + 2] = VR;

227:     Vref[i] = Vm + (VR / KA[i]);

229:     idx = idx + 3;
230:   }

232:   PetscCall(VecRestoreArray(Xgen, &xgen));
233:   PetscCall(VecRestoreArray(Xnet, &xnet));

235:   /* PetscCall(VecView(Xgen,0)); */
236:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
237:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
238:   PetscFunctionReturn(PETSC_SUCCESS);
239: }

241: PetscErrorCode InitialGuess(Vec X, Userctx *user, const PetscScalar PGv[])
242: {
243:   Vec          Xgen, Xnet;
244:   PetscScalar *xgen, *xnet;
245:   PetscInt     i, idx = 0;
246:   PetscScalar  Vr, Vi, IGr, IGi, Vm, Vm2;
247:   PetscScalar  Eqp, Edp, delta;
248:   PetscScalar  Efd, RF, VR; /* Exciter variables */
249:   PetscScalar  Id, Iq;      /* Generator dq axis currents */
250:   PetscScalar  theta, Vd, Vq, SE;

252:   PetscFunctionBegin;
253:   M[0] = 2 * H[0] / w_s;
254:   M[1] = 2 * H[1] / w_s;
255:   M[2] = 2 * H[2] / w_s;
256:   D[0] = 0.1 * M[0];
257:   D[1] = 0.1 * M[1];
258:   D[2] = 0.1 * M[2];

260:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));

262:   /* Network subsystem initialization */
263:   PetscCall(VecCopy(user->V0, Xnet));

265:   /* Generator subsystem initialization */
266:   PetscCall(VecGetArray(Xgen, &xgen));
267:   PetscCall(VecGetArray(Xnet, &xnet));

269:   for (i = 0; i < ngen; i++) {
270:     Vr  = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
271:     Vi  = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
272:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
273:     Vm2 = Vm * Vm;
274:     IGr = (Vr * PGv[i] + Vi * QG[i]) / Vm2;
275:     IGi = (Vi * PGv[i] - Vr * QG[i]) / Vm2;

277:     delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */

279:     theta = PETSC_PI / 2.0 - delta;

281:     Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
282:     Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */

284:     Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
285:     Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);

287:     Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
288:     Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */

290:     /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
291:     xgen[idx]     = Eqp;
292:     xgen[idx + 1] = Edp;
293:     xgen[idx + 2] = delta;
294:     xgen[idx + 3] = w_s;

296:     idx = idx + 4;

298:     xgen[idx]     = Id;
299:     xgen[idx + 1] = Iq;

301:     idx = idx + 2;

303:     /* Exciter */
304:     Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
305:     SE  = k1[i] * PetscExpScalar(k2[i] * Efd);
306:     VR  = KE[i] * Efd + SE;
307:     RF  = KF[i] * Efd / TF[i];

309:     xgen[idx]     = Efd;
310:     xgen[idx + 1] = RF;
311:     xgen[idx + 2] = VR;

313:     idx = idx + 3;
314:   }

316:   PetscCall(VecRestoreArray(Xgen, &xgen));
317:   PetscCall(VecRestoreArray(Xnet, &xnet));

319:   /* PetscCall(VecView(Xgen,0)); */
320:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
321:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

325: PetscErrorCode DICDPFiniteDifference(Vec X, Vec *DICDP, Userctx *user)
326: {
327:   Vec         Y;
328:   PetscScalar PGv[3], eps;
329:   PetscInt    i, j;

331:   PetscFunctionBegin;
332:   eps = 1.e-7;
333:   PetscCall(VecDuplicate(X, &Y));

335:   for (i = 0; i < ngen; i++) {
336:     for (j = 0; j < 3; j++) PGv[j] = PG[j];
337:     PGv[i] = PG[i] + eps;
338:     PetscCall(InitialGuess(Y, user, PGv));
339:     PetscCall(InitialGuess(X, user, PG));

341:     PetscCall(VecAXPY(Y, -1.0, X));
342:     PetscCall(VecScale(Y, 1. / eps));
343:     PetscCall(VecCopy(Y, DICDP[i]));
344:   }
345:   PetscCall(VecDestroy(&Y));
346:   PetscFunctionReturn(PETSC_SUCCESS);
347: }

349: /* Computes F = [-f(x,y);g(x,y)] */
350: PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
351: {
352:   Vec          Xgen, Xnet, Fgen, Fnet;
353:   PetscScalar *xgen, *xnet, *fgen, *fnet;
354:   PetscInt     i, idx = 0;
355:   PetscScalar  Vr, Vi, Vm, Vm2;
356:   PetscScalar  Eqp, Edp, delta, w; /* Generator variables */
357:   PetscScalar  Efd, RF, VR;        /* Exciter variables */
358:   PetscScalar  Id, Iq;             /* Generator dq axis currents */
359:   PetscScalar  Vd, Vq, SE;
360:   PetscScalar  IGr, IGi, IDr, IDi;
361:   PetscScalar  Zdq_inv[4], det;
362:   PetscScalar  PD, QD, Vm0, *v0;
363:   PetscInt     k;

365:   PetscFunctionBegin;
366:   PetscCall(VecZeroEntries(F));
367:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
368:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
369:   PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
370:   PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));

372:   /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
373:      The generator current injection, IG, and load current injection, ID are added later
374:   */
375:   /* Note that the values in Ybus are stored assuming the imaginary current balance
376:      equation is ordered first followed by real current balance equation for each bus.
377:      Thus imaginary current contribution goes in location 2*i, and
378:      real current contribution in 2*i+1
379:   */
380:   PetscCall(MatMult(user->Ybus, Xnet, Fnet));

382:   PetscCall(VecGetArray(Xgen, &xgen));
383:   PetscCall(VecGetArray(Xnet, &xnet));
384:   PetscCall(VecGetArray(Fgen, &fgen));
385:   PetscCall(VecGetArray(Fnet, &fnet));

387:   /* Generator subsystem */
388:   for (i = 0; i < ngen; i++) {
389:     Eqp   = xgen[idx];
390:     Edp   = xgen[idx + 1];
391:     delta = xgen[idx + 2];
392:     w     = xgen[idx + 3];
393:     Id    = xgen[idx + 4];
394:     Iq    = xgen[idx + 5];
395:     Efd   = xgen[idx + 6];
396:     RF    = xgen[idx + 7];
397:     VR    = xgen[idx + 8];

399:     /* Generator differential equations */
400:     fgen[idx]     = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
401:     fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
402:     fgen[idx + 2] = -w + w_s;
403:     fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];

405:     Vr = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
406:     Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */

408:     PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
409:     /* Algebraic equations for stator currents */
410:     det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];

412:     Zdq_inv[0] = Rs[i] / det;
413:     Zdq_inv[1] = Xqp[i] / det;
414:     Zdq_inv[2] = -Xdp[i] / det;
415:     Zdq_inv[3] = Rs[i] / det;

417:     fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
418:     fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;

420:     /* Add generator current injection to network */
421:     PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));

423:     fnet[2 * gbus[i]] -= IGi;
424:     fnet[2 * gbus[i] + 1] -= IGr;

426:     Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);

428:     SE = k1[i] * PetscExpScalar(k2[i] * Efd);

430:     /* Exciter differential equations */
431:     fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
432:     fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
433:     fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];

435:     idx = idx + 9;
436:   }

438:   PetscCall(VecGetArray(user->V0, &v0));
439:   for (i = 0; i < nload; i++) {
440:     Vr  = xnet[2 * lbus[i]];     /* Real part of load bus voltage */
441:     Vi  = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
442:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
443:     Vm2 = Vm * Vm;
444:     Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
445:     PD = QD = 0.0;
446:     for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
447:     for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);

449:     /* Load currents */
450:     IDr = (PD * Vr + QD * Vi) / Vm2;
451:     IDi = (-QD * Vr + PD * Vi) / Vm2;

453:     fnet[2 * lbus[i]] += IDi;
454:     fnet[2 * lbus[i] + 1] += IDr;
455:   }
456:   PetscCall(VecRestoreArray(user->V0, &v0));

458:   PetscCall(VecRestoreArray(Xgen, &xgen));
459:   PetscCall(VecRestoreArray(Xnet, &xnet));
460:   PetscCall(VecRestoreArray(Fgen, &fgen));
461:   PetscCall(VecRestoreArray(Fnet, &fnet));

463:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
464:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
465:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
466:   PetscFunctionReturn(PETSC_SUCCESS);
467: }

469: /* \dot{x} - f(x,y)
470:      g(x,y) = 0
471:  */
472: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
473: {
474:   SNES               snes;
475:   PetscScalar       *f;
476:   const PetscScalar *xdot;
477:   PetscInt           i;

479:   PetscFunctionBegin;
480:   user->t = t;

482:   PetscCall(TSGetSNES(ts, &snes));
483:   PetscCall(ResidualFunction(snes, X, F, user));
484:   PetscCall(VecGetArray(F, &f));
485:   PetscCall(VecGetArrayRead(Xdot, &xdot));
486:   for (i = 0; i < ngen; i++) {
487:     f[9 * i] += xdot[9 * i];
488:     f[9 * i + 1] += xdot[9 * i + 1];
489:     f[9 * i + 2] += xdot[9 * i + 2];
490:     f[9 * i + 3] += xdot[9 * i + 3];
491:     f[9 * i + 6] += xdot[9 * i + 6];
492:     f[9 * i + 7] += xdot[9 * i + 7];
493:     f[9 * i + 8] += xdot[9 * i + 8];
494:   }
495:   PetscCall(VecRestoreArray(F, &f));
496:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
497:   PetscFunctionReturn(PETSC_SUCCESS);
498: }

500: /* This function is used for solving the algebraic system only during fault on and
501:    off times. It computes the entire F and then zeros out the part corresponding to
502:    differential equations
503:  F = [0;g(y)];
504: */
505: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
506: {
507:   Userctx     *user = (Userctx *)ctx;
508:   PetscScalar *f;
509:   PetscInt     i;

511:   PetscFunctionBegin;
512:   PetscCall(ResidualFunction(snes, X, F, user));
513:   PetscCall(VecGetArray(F, &f));
514:   for (i = 0; i < ngen; i++) {
515:     f[9 * i]     = 0;
516:     f[9 * i + 1] = 0;
517:     f[9 * i + 2] = 0;
518:     f[9 * i + 3] = 0;
519:     f[9 * i + 6] = 0;
520:     f[9 * i + 7] = 0;
521:     f[9 * i + 8] = 0;
522:   }
523:   PetscCall(VecRestoreArray(F, &f));
524:   PetscFunctionReturn(PETSC_SUCCESS);
525: }

527: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
528: {
529:   PetscInt *d_nnz;
530:   PetscInt  i, idx = 0, start = 0;
531:   PetscInt  ncols;

533:   PetscFunctionBegin;
534:   PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
535:   for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
536:   /* Generator subsystem */
537:   for (i = 0; i < ngen; i++) {
538:     d_nnz[idx] += 3;
539:     d_nnz[idx + 1] += 2;
540:     d_nnz[idx + 2] += 2;
541:     d_nnz[idx + 3] += 5;
542:     d_nnz[idx + 4] += 6;
543:     d_nnz[idx + 5] += 6;

545:     d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
546:     d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;

548:     d_nnz[idx + 6] += 2;
549:     d_nnz[idx + 7] += 2;
550:     d_nnz[idx + 8] += 5;

552:     idx = idx + 9;
553:   }

555:   start = user->neqs_gen;
556:   for (i = 0; i < nbus; i++) {
557:     PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
558:     d_nnz[start + 2 * i] += ncols;
559:     d_nnz[start + 2 * i + 1] += ncols;
560:     PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
561:   }

563:   PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
564:   PetscCall(PetscFree(d_nnz));
565:   PetscFunctionReturn(PETSC_SUCCESS);
566: }

568: /*
569:    J = [-df_dx, -df_dy
570:         dg_dx, dg_dy]
571: */
572: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, void *ctx)
573: {
574:   Userctx           *user = (Userctx *)ctx;
575:   Vec                Xgen, Xnet;
576:   PetscScalar       *xgen, *xnet;
577:   PetscInt           i, idx = 0;
578:   PetscScalar        Vr, Vi, Vm, Vm2;
579:   PetscScalar        Eqp, Edp, delta; /* Generator variables */
580:   PetscScalar        Efd;             /* Exciter variables */
581:   PetscScalar        Id, Iq;          /* Generator dq axis currents */
582:   PetscScalar        Vd, Vq;
583:   PetscScalar        val[10];
584:   PetscInt           row[2], col[10];
585:   PetscInt           net_start = user->neqs_gen;
586:   PetscInt           ncols;
587:   const PetscInt    *cols;
588:   const PetscScalar *yvals;
589:   PetscInt           k;
590:   PetscScalar        Zdq_inv[4], det;
591:   PetscScalar        dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
592:   PetscScalar        dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
593:   PetscScalar        dSE_dEfd;
594:   PetscScalar        dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
595:   PetscScalar        PD, QD, Vm0, *v0, Vm4;
596:   PetscScalar        dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
597:   PetscScalar        dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;

599:   PetscFunctionBegin;
600:   PetscCall(MatZeroEntries(B));
601:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
602:   PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));

604:   PetscCall(VecGetArray(Xgen, &xgen));
605:   PetscCall(VecGetArray(Xnet, &xnet));

607:   /* Generator subsystem */
608:   for (i = 0; i < ngen; i++) {
609:     Eqp   = xgen[idx];
610:     Edp   = xgen[idx + 1];
611:     delta = xgen[idx + 2];
612:     Id    = xgen[idx + 4];
613:     Iq    = xgen[idx + 5];
614:     Efd   = xgen[idx + 6];

616:     /*    fgen[idx]   = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
617:     row[0] = idx;
618:     col[0] = idx;
619:     col[1] = idx + 4;
620:     col[2] = idx + 6;
621:     val[0] = 1 / Td0p[i];
622:     val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
623:     val[2] = -1 / Td0p[i];

625:     PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));

627:     /*    fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
628:     row[0] = idx + 1;
629:     col[0] = idx + 1;
630:     col[1] = idx + 5;
631:     val[0] = 1 / Tq0p[i];
632:     val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
633:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

635:     /*    fgen[idx+2] = - w + w_s; */
636:     row[0] = idx + 2;
637:     col[0] = idx + 2;
638:     col[1] = idx + 3;
639:     val[0] = 0;
640:     val[1] = -1;
641:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

643:     /*    fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
644:     row[0] = idx + 3;
645:     col[0] = idx;
646:     col[1] = idx + 1;
647:     col[2] = idx + 3;
648:     col[3] = idx + 4;
649:     col[4] = idx + 5;
650:     val[0] = Iq / M[i];
651:     val[1] = Id / M[i];
652:     val[2] = D[i] / M[i];
653:     val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
654:     val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
655:     PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));

657:     Vr = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
658:     Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
659:     PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));

661:     det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];

663:     Zdq_inv[0] = Rs[i] / det;
664:     Zdq_inv[1] = Xqp[i] / det;
665:     Zdq_inv[2] = -Xdp[i] / det;
666:     Zdq_inv[3] = Rs[i] / det;

668:     dVd_dVr    = PetscSinScalar(delta);
669:     dVd_dVi    = -PetscCosScalar(delta);
670:     dVq_dVr    = PetscCosScalar(delta);
671:     dVq_dVi    = PetscSinScalar(delta);
672:     dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
673:     dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);

675:     /*    fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
676:     row[0] = idx + 4;
677:     col[0] = idx;
678:     col[1] = idx + 1;
679:     col[2] = idx + 2;
680:     val[0] = -Zdq_inv[1];
681:     val[1] = -Zdq_inv[0];
682:     val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
683:     col[3] = idx + 4;
684:     col[4] = net_start + 2 * gbus[i];
685:     col[5] = net_start + 2 * gbus[i] + 1;
686:     val[3] = 1;
687:     val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
688:     val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
689:     PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));

691:     /*  fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
692:     row[0] = idx + 5;
693:     col[0] = idx;
694:     col[1] = idx + 1;
695:     col[2] = idx + 2;
696:     val[0] = -Zdq_inv[3];
697:     val[1] = -Zdq_inv[2];
698:     val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
699:     col[3] = idx + 5;
700:     col[4] = net_start + 2 * gbus[i];
701:     col[5] = net_start + 2 * gbus[i] + 1;
702:     val[3] = 1;
703:     val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
704:     val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
705:     PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));

707:     dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
708:     dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
709:     dIGr_dId    = PetscSinScalar(delta);
710:     dIGr_dIq    = PetscCosScalar(delta);
711:     dIGi_dId    = -PetscCosScalar(delta);
712:     dIGi_dIq    = PetscSinScalar(delta);

714:     /* fnet[2*gbus[i]]   -= IGi; */
715:     row[0] = net_start + 2 * gbus[i];
716:     col[0] = idx + 2;
717:     col[1] = idx + 4;
718:     col[2] = idx + 5;
719:     val[0] = -dIGi_ddelta;
720:     val[1] = -dIGi_dId;
721:     val[2] = -dIGi_dIq;
722:     PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));

724:     /* fnet[2*gbus[i]+1]   -= IGr; */
725:     row[0] = net_start + 2 * gbus[i] + 1;
726:     col[0] = idx + 2;
727:     col[1] = idx + 4;
728:     col[2] = idx + 5;
729:     val[0] = -dIGr_ddelta;
730:     val[1] = -dIGr_dId;
731:     val[2] = -dIGr_dIq;
732:     PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));

734:     Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);

736:     /*    fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
737:     /*    SE  = k1[i]*PetscExpScalar(k2[i]*Efd); */
738:     dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);

740:     row[0] = idx + 6;
741:     col[0] = idx + 6;
742:     col[1] = idx + 8;
743:     val[0] = (KE[i] + dSE_dEfd) / TE[i];
744:     val[1] = -1 / TE[i];
745:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

747:     /* Exciter differential equations */

749:     /*    fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
750:     row[0] = idx + 7;
751:     col[0] = idx + 6;
752:     col[1] = idx + 7;
753:     val[0] = (-KF[i] / TF[i]) / TF[i];
754:     val[1] = 1 / TF[i];
755:     PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));

757:     /*    fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
758:     /* Vm = (Vd^2 + Vq^2)^0.5; */
759:     dVm_dVd = Vd / Vm;
760:     dVm_dVq = Vq / Vm;
761:     dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
762:     dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
763:     row[0]  = idx + 8;
764:     col[0]  = idx + 6;
765:     col[1]  = idx + 7;
766:     col[2]  = idx + 8;
767:     val[0]  = (KA[i] * KF[i] / TF[i]) / TA[i];
768:     val[1]  = -KA[i] / TA[i];
769:     val[2]  = 1 / TA[i];
770:     col[3]  = net_start + 2 * gbus[i];
771:     col[4]  = net_start + 2 * gbus[i] + 1;
772:     val[3]  = KA[i] * dVm_dVr / TA[i];
773:     val[4]  = KA[i] * dVm_dVi / TA[i];
774:     PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
775:     idx = idx + 9;
776:   }

778:   for (i = 0; i < nbus; i++) {
779:     PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
780:     row[0] = net_start + 2 * i;
781:     for (k = 0; k < ncols; k++) {
782:       col[k] = net_start + cols[k];
783:       val[k] = yvals[k];
784:     }
785:     PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
786:     PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));

788:     PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
789:     row[0] = net_start + 2 * i + 1;
790:     for (k = 0; k < ncols; k++) {
791:       col[k] = net_start + cols[k];
792:       val[k] = yvals[k];
793:     }
794:     PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
795:     PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
796:   }

798:   PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
799:   PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));

801:   PetscCall(VecGetArray(user->V0, &v0));
802:   for (i = 0; i < nload; i++) {
803:     Vr  = xnet[2 * lbus[i]];     /* Real part of load bus voltage */
804:     Vi  = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
805:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
806:     Vm2 = Vm * Vm;
807:     Vm4 = Vm2 * Vm2;
808:     Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
809:     PD = QD = 0.0;
810:     dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
811:     for (k = 0; k < ld_nsegsp[i]; k++) {
812:       PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
813:       dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vr * PetscPowScalar(Vm, (ld_betap[k] - 2));
814:       dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vi * PetscPowScalar(Vm, (ld_betap[k] - 2));
815:     }
816:     for (k = 0; k < ld_nsegsq[i]; k++) {
817:       QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
818:       dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vr * PetscPowScalar(Vm, (ld_betaq[k] - 2));
819:       dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vi * PetscPowScalar(Vm, (ld_betaq[k] - 2));
820:     }

822:     /*    IDr = (PD*Vr + QD*Vi)/Vm2; */
823:     /*    IDi = (-QD*Vr + PD*Vi)/Vm2; */

825:     dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
826:     dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;

828:     dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
829:     dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;

831:     /*    fnet[2*lbus[i]]   += IDi; */
832:     row[0] = net_start + 2 * lbus[i];
833:     col[0] = net_start + 2 * lbus[i];
834:     col[1] = net_start + 2 * lbus[i] + 1;
835:     val[0] = dIDi_dVr;
836:     val[1] = dIDi_dVi;
837:     PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
838:     /*    fnet[2*lbus[i]+1] += IDr; */
839:     row[0] = net_start + 2 * lbus[i] + 1;
840:     col[0] = net_start + 2 * lbus[i];
841:     col[1] = net_start + 2 * lbus[i] + 1;
842:     val[0] = dIDr_dVr;
843:     val[1] = dIDr_dVi;
844:     PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
845:   }
846:   PetscCall(VecRestoreArray(user->V0, &v0));

848:   PetscCall(VecRestoreArray(Xgen, &xgen));
849:   PetscCall(VecRestoreArray(Xnet, &xnet));

851:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));

853:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
854:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
855:   PetscFunctionReturn(PETSC_SUCCESS);
856: }

858: /*
859:    J = [I, 0
860:         dg_dx, dg_dy]
861: */
862: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, void *ctx)
863: {
864:   Userctx *user = (Userctx *)ctx;

866:   PetscFunctionBegin;
867:   PetscCall(ResidualJacobian(snes, X, A, B, ctx));
868:   PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
869:   PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
870:   PetscFunctionReturn(PETSC_SUCCESS);
871: }

873: /*
874:    J = [a*I-df_dx, -df_dy
875:         dg_dx, dg_dy]
876: */

878: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
879: {
880:   SNES        snes;
881:   PetscScalar atmp = (PetscScalar)a;
882:   PetscInt    i, row;

884:   PetscFunctionBegin;
885:   user->t = t;

887:   PetscCall(TSGetSNES(ts, &snes));
888:   PetscCall(ResidualJacobian(snes, X, A, B, user));
889:   for (i = 0; i < ngen; i++) {
890:     row = 9 * i;
891:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
892:     row = 9 * i + 1;
893:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
894:     row = 9 * i + 2;
895:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
896:     row = 9 * i + 3;
897:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
898:     row = 9 * i + 6;
899:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
900:     row = 9 * i + 7;
901:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
902:     row = 9 * i + 8;
903:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
904:   }
905:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
906:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
907:   PetscFunctionReturn(PETSC_SUCCESS);
908: }

910: /* Matrix JacobianP is constant so that it only needs to be evaluated once */
911: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
912: {
913:   PetscScalar a;
914:   PetscInt    row, col;
915:   Userctx    *ctx = (Userctx *)ctx0;

917:   PetscFunctionBeginUser;

919:   if (ctx->jacp_flg) {
920:     PetscCall(MatZeroEntries(A));

922:     for (col = 0; col < 3; col++) {
923:       a   = 1.0 / M[col];
924:       row = 9 * col + 3;
925:       PetscCall(MatSetValues(A, 1, &row, 1, &col, &a, INSERT_VALUES));
926:     }

928:     ctx->jacp_flg = PETSC_FALSE;

930:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
931:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
932:   }
933:   PetscFunctionReturn(PETSC_SUCCESS);
934: }

936: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
937: {
938:   const PetscScalar *u;
939:   PetscInt           idx;
940:   Vec                Xgen, Xnet;
941:   PetscScalar       *r, *xgen;
942:   PetscInt           i;

944:   PetscFunctionBegin;
945:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
946:   PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));

948:   PetscCall(VecGetArray(Xgen, &xgen));

950:   PetscCall(VecGetArrayRead(U, &u));
951:   PetscCall(VecGetArray(R, &r));
952:   r[0] = 0.;
953:   idx  = 0;
954:   for (i = 0; i < ngen; i++) {
955:     r[0] += PetscPowScalarInt(PetscMax(0., PetscMax(xgen[idx + 3] / (2. * PETSC_PI) - user->freq_u, user->freq_l - xgen[idx + 3] / (2. * PETSC_PI))), user->pow);
956:     idx += 9;
957:   }
958:   PetscCall(VecRestoreArrayRead(U, &u));
959:   PetscCall(VecRestoreArray(R, &r));
960:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
961:   PetscFunctionReturn(PETSC_SUCCESS);
962: }

964: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, Userctx *user)
965: {
966:   Vec          Xgen, Xnet, Dgen, Dnet;
967:   PetscScalar *xgen, *dgen;
968:   PetscInt     i;
969:   PetscInt     idx;
970:   Vec          drdu_col;
971:   PetscScalar *xarr;

973:   PetscFunctionBegin;
974:   PetscCall(VecDuplicate(U, &drdu_col));
975:   PetscCall(MatDenseGetColumn(DRDU, 0, &xarr));
976:   PetscCall(VecPlaceArray(drdu_col, xarr));
977:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
978:   PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Dgen, &Dnet));
979:   PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
980:   PetscCall(DMCompositeScatter(user->dmpgrid, drdu_col, Dgen, Dnet));

982:   PetscCall(VecGetArray(Xgen, &xgen));
983:   PetscCall(VecGetArray(Dgen, &dgen));

985:   idx = 0;
986:   for (i = 0; i < ngen; i++) {
987:     dgen[idx + 3] = 0.;
988:     if (xgen[idx + 3] / (2. * PETSC_PI) > user->freq_u) dgen[idx + 3] = user->pow * PetscPowScalarInt(xgen[idx + 3] / (2. * PETSC_PI) - user->freq_u, user->pow - 1) / (2. * PETSC_PI);
989:     if (xgen[idx + 3] / (2. * PETSC_PI) < user->freq_l) dgen[idx + 3] = user->pow * PetscPowScalarInt(user->freq_l - xgen[idx + 3] / (2. * PETSC_PI), user->pow - 1) / (-2. * PETSC_PI);
990:     idx += 9;
991:   }

993:   PetscCall(VecRestoreArray(Dgen, &dgen));
994:   PetscCall(VecRestoreArray(Xgen, &xgen));
995:   PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, drdu_col, Dgen, Dnet));
996:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Dgen, &Dnet));
997:   PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
998:   PetscCall(VecResetArray(drdu_col));
999:   PetscCall(MatDenseRestoreColumn(DRDU, &xarr));
1000:   PetscCall(VecDestroy(&drdu_col));
1001:   PetscFunctionReturn(PETSC_SUCCESS);
1002: }

1004: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat drdp, Userctx *user)
1005: {
1006:   PetscFunctionBegin;
1007:   PetscFunctionReturn(PETSC_SUCCESS);
1008: }

1010: PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, Vec *DICDP, Userctx *user)
1011: {
1012:   PetscScalar *x, *y, sensip;
1013:   PetscInt     i;

1015:   PetscFunctionBegin;
1016:   PetscCall(VecGetArray(lambda, &x));
1017:   PetscCall(VecGetArray(mu, &y));

1019:   for (i = 0; i < 3; i++) {
1020:     PetscCall(VecDot(lambda, DICDP[i], &sensip));
1021:     sensip = sensip + y[i];
1022:     /* PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt %" PetscInt_FMT " th parameter: %g \n",i,(double)sensip)); */
1023:     y[i] = sensip;
1024:   }
1025:   PetscCall(VecRestoreArray(mu, &y));
1026:   PetscFunctionReturn(PETSC_SUCCESS);
1027: }

1029: int main(int argc, char **argv)
1030: {
1031:   Userctx      user;
1032:   Vec          p;
1033:   PetscScalar *x_ptr;
1034:   PetscMPIInt  size;
1035:   PetscInt     i;
1036:   PetscViewer  Xview, Ybusview;
1037:   PetscInt    *idx2;
1038:   Tao          tao;
1039:   KSP          ksp;
1040:   PC           pc;
1041:   Vec          lowerb, upperb;

1043:   PetscFunctionBeginUser;
1044:   PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
1045:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1046:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

1048:   user.jacp_flg   = PETSC_TRUE;
1049:   user.neqs_gen   = 9 * ngen; /* # eqs. for generator subsystem */
1050:   user.neqs_net   = 2 * nbus; /* # eqs. for network subsystem   */
1051:   user.neqs_pgrid = user.neqs_gen + user.neqs_net;

1053:   /* Create indices for differential and algebraic equations */
1054:   PetscCall(PetscMalloc1(7 * ngen, &idx2));
1055:   for (i = 0; i < ngen; i++) {
1056:     idx2[7 * i]     = 9 * i;
1057:     idx2[7 * i + 1] = 9 * i + 1;
1058:     idx2[7 * i + 2] = 9 * i + 2;
1059:     idx2[7 * i + 3] = 9 * i + 3;
1060:     idx2[7 * i + 4] = 9 * i + 6;
1061:     idx2[7 * i + 5] = 9 * i + 7;
1062:     idx2[7 * i + 6] = 9 * i + 8;
1063:   }
1064:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
1065:   PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
1066:   PetscCall(PetscFree(idx2));

1068:   /* Set run time options */
1069:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1070:   {
1071:     user.tfaulton  = 1.0;
1072:     user.tfaultoff = 1.2;
1073:     user.Rfault    = 0.0001;
1074:     user.faultbus  = 8;
1075:     PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
1076:     PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
1077:     PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
1078:     user.t0   = 0.0;
1079:     user.tmax = 1.3;
1080:     PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
1081:     PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
1082:     user.freq_u = 61.0;
1083:     user.freq_l = 59.0;
1084:     user.pow    = 2;
1085:     PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
1086:     PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
1087:     PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
1088:   }
1089:   PetscOptionsEnd();

1091:   /* Create DMs for generator and network subsystems */
1092:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
1093:   PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
1094:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
1095:   PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
1096:   PetscCall(DMSetFromOptions(user.dmnet));
1097:   PetscCall(DMSetUp(user.dmnet));
1098:   /* Create a composite DM packer and add the two DMs */
1099:   PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
1100:   PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
1101:   PetscCall(DMSetFromOptions(user.dmgen));
1102:   PetscCall(DMSetUp(user.dmgen));
1103:   PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
1104:   PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));

1106:   /* Read initial voltage vector and Ybus */
1107:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
1108:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));

1110:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user.V0));
1111:   PetscCall(VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net));
1112:   PetscCall(VecLoad(user.V0, Xview));

1114:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Ybus));
1115:   PetscCall(MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net));
1116:   PetscCall(MatSetType(user.Ybus, MATBAIJ));
1117:   /*  PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
1118:   PetscCall(MatLoad(user.Ybus, Ybusview));

1120:   PetscCall(PetscViewerDestroy(&Xview));
1121:   PetscCall(PetscViewerDestroy(&Ybusview));

1123:   /* Allocate space for Jacobians */
1124:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.J));
1125:   PetscCall(MatSetSizes(user.J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid));
1126:   PetscCall(MatSetFromOptions(user.J));
1127:   PetscCall(PreallocateJacobian(user.J, &user));

1129:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
1130:   PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 3));
1131:   PetscCall(MatSetFromOptions(user.Jacp));
1132:   PetscCall(MatSetUp(user.Jacp));
1133:   PetscCall(MatZeroEntries(user.Jacp)); /* initialize to zeros */

1135:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, &user.DRDP));
1136:   PetscCall(MatSetUp(user.DRDP));
1137:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, 1, NULL, &user.DRDU));
1138:   PetscCall(MatSetUp(user.DRDU));

1140:   /* Create TAO solver and set desired solution method */
1141:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
1142:   PetscCall(TaoSetType(tao, TAOBLMVM));
1143:   /*
1144:      Optimization starts
1145:   */
1146:   /* Set initial solution guess */
1147:   PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
1148:   PetscCall(VecGetArray(p, &x_ptr));
1149:   x_ptr[0] = PG[0];
1150:   x_ptr[1] = PG[1];
1151:   x_ptr[2] = PG[2];
1152:   PetscCall(VecRestoreArray(p, &x_ptr));

1154:   PetscCall(TaoSetSolution(tao, p));
1155:   /* Set routine for function and gradient evaluation */
1156:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user));

1158:   /* Set bounds for the optimization */
1159:   PetscCall(VecDuplicate(p, &lowerb));
1160:   PetscCall(VecDuplicate(p, &upperb));
1161:   PetscCall(VecGetArray(lowerb, &x_ptr));
1162:   x_ptr[0] = 0.5;
1163:   x_ptr[1] = 0.5;
1164:   x_ptr[2] = 0.5;
1165:   PetscCall(VecRestoreArray(lowerb, &x_ptr));
1166:   PetscCall(VecGetArray(upperb, &x_ptr));
1167:   x_ptr[0] = 2.0;
1168:   x_ptr[1] = 2.0;
1169:   x_ptr[2] = 2.0;
1170:   PetscCall(VecRestoreArray(upperb, &x_ptr));
1171:   PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));

1173:   /* Check for any TAO command line options */
1174:   PetscCall(TaoSetFromOptions(tao));
1175:   PetscCall(TaoGetKSP(tao, &ksp));
1176:   if (ksp) {
1177:     PetscCall(KSPGetPC(ksp, &pc));
1178:     PetscCall(PCSetType(pc, PCNONE));
1179:   }

1181:   /* SOLVE THE APPLICATION */
1182:   PetscCall(TaoSolve(tao));

1184:   PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
1185:   /* Free TAO data structures */
1186:   PetscCall(TaoDestroy(&tao));

1188:   PetscCall(DMDestroy(&user.dmgen));
1189:   PetscCall(DMDestroy(&user.dmnet));
1190:   PetscCall(DMDestroy(&user.dmpgrid));
1191:   PetscCall(ISDestroy(&user.is_diff));
1192:   PetscCall(ISDestroy(&user.is_alg));

1194:   PetscCall(MatDestroy(&user.J));
1195:   PetscCall(MatDestroy(&user.Jacp));
1196:   PetscCall(MatDestroy(&user.Ybus));
1197:   /* PetscCall(MatDestroy(&user.Sol)); */
1198:   PetscCall(VecDestroy(&user.V0));
1199:   PetscCall(VecDestroy(&p));
1200:   PetscCall(VecDestroy(&lowerb));
1201:   PetscCall(VecDestroy(&upperb));
1202:   PetscCall(MatDestroy(&user.DRDU));
1203:   PetscCall(MatDestroy(&user.DRDP));
1204:   PetscCall(PetscFinalize());
1205:   return 0;
1206: }

1208: /* ------------------------------------------------------------------ */
1209: /*
1210:    FormFunction - Evaluates the function and corresponding gradient.

1212:    Input Parameters:
1213:    tao - the Tao context
1214:    X   - the input vector
1215:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

1217:    Output Parameters:
1218:    f   - the newly evaluated function
1219:    G   - the newly evaluated gradient
1220: */
1221: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx0)
1222: {
1223:   TS       ts, quadts;
1224:   SNES     snes_alg;
1225:   Userctx *ctx = (Userctx *)ctx0;
1226:   Vec      X;
1227:   PetscInt i;
1228:   /* sensitivity context */
1229:   PetscScalar *x_ptr;
1230:   Vec          lambda[1], q;
1231:   Vec          mu[1];
1232:   PetscInt     steps1, steps2, steps3;
1233:   Vec          DICDP[3];
1234:   Vec          F_alg;
1235:   PetscInt     row_loc, col_loc;
1236:   PetscScalar  val;
1237:   Vec          Xdot;

1239:   PetscFunctionBegin;
1240:   PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
1241:   PG[0] = x_ptr[0];
1242:   PG[1] = x_ptr[1];
1243:   PG[2] = x_ptr[2];
1244:   PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));

1246:   ctx->stepnum = 0;

1248:   PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));

1250:   /* Create matrix to save solutions at each time step */
1251:   /* PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ctx->neqs_pgrid+1,1002,NULL,&ctx->Sol)); */
1252:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1253:      Create timestepping solver context
1254:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1255:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1256:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1257:   PetscCall(TSSetType(ts, TSCN));
1258:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)IFunction, ctx));
1259:   PetscCall(TSSetIJacobian(ts, ctx->J, ctx->J, (TSIJacobian)IJacobian, ctx));
1260:   PetscCall(TSSetApplicationContext(ts, ctx));
1261:   /*   Set RHS JacobianP */
1262:   PetscCall(TSSetRHSJacobianP(ts, ctx->Jacp, RHSJacobianP, ctx));

1264:   PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
1265:   PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunction)CostIntegrand, ctx));
1266:   PetscCall(TSSetRHSJacobian(quadts, ctx->DRDU, ctx->DRDU, (TSRHSJacobian)DRDUJacobianTranspose, ctx));
1267:   PetscCall(TSSetRHSJacobianP(quadts, ctx->DRDP, (TSRHSJacobianP)DRDPJacobianTranspose, ctx));

1269:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1270:      Set initial conditions
1271:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1272:   PetscCall(SetInitialGuess(X, ctx));

1274:   /* Approximate DICDP with finite difference, we want to zero out network variables */
1275:   for (i = 0; i < 3; i++) PetscCall(VecDuplicate(X, &DICDP[i]));
1276:   PetscCall(DICDPFiniteDifference(X, DICDP, ctx));

1278:   PetscCall(VecDuplicate(X, &F_alg));
1279:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1280:   PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
1281:   PetscCall(MatZeroEntries(ctx->J));
1282:   PetscCall(SNESSetJacobian(snes_alg, ctx->J, ctx->J, AlgJacobian, ctx));
1283:   PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1284:   PetscCall(SNESSetFromOptions(snes_alg));
1285:   ctx->alg_flg = PETSC_TRUE;
1286:   /* Solve the algebraic equations */
1287:   PetscCall(SNESSolve(snes_alg, NULL, X));

1289:   /* Just to set up the Jacobian structure */
1290:   PetscCall(VecDuplicate(X, &Xdot));
1291:   PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, ctx->J, ctx->J, ctx));
1292:   PetscCall(VecDestroy(&Xdot));

1294:   ctx->stepnum++;

1296:   /*
1297:     Save trajectory of solution so that TSAdjointSolve() may be used
1298:   */
1299:   PetscCall(TSSetSaveTrajectory(ts));

1301:   PetscCall(TSSetTimeStep(ts, 0.01));
1302:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1303:   PetscCall(TSSetFromOptions(ts));
1304:   /* PetscCall(TSSetPostStep(ts,SaveSolution)); */

1306:   /* Prefault period */
1307:   ctx->alg_flg = PETSC_FALSE;
1308:   PetscCall(TSSetTime(ts, 0.0));
1309:   PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
1310:   PetscCall(TSSolve(ts, X));
1311:   PetscCall(TSGetStepNumber(ts, &steps1));

1313:   /* Create the nonlinear solver for solving the algebraic system */
1314:   /* Note that although the algebraic system needs to be solved only for
1315:      Idq and V, we reuse the entire system including xgen. The xgen
1316:      variables are held constant by setting their residuals to 0 and
1317:      putting a 1 on the Jacobian diagonal for xgen rows
1318:   */
1319:   PetscCall(MatZeroEntries(ctx->J));

1321:   /* Apply disturbance - resistive fault at ctx->faultbus */
1322:   /* This is done by adding shunt conductance to the diagonal location
1323:      in the Ybus matrix */
1324:   row_loc = 2 * ctx->faultbus;
1325:   col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1326:   val     = 1 / ctx->Rfault;
1327:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1328:   row_loc = 2 * ctx->faultbus + 1;
1329:   col_loc = 2 * ctx->faultbus; /* Location for G */
1330:   val     = 1 / ctx->Rfault;
1331:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

1333:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1334:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

1336:   ctx->alg_flg = PETSC_TRUE;
1337:   /* Solve the algebraic equations */
1338:   PetscCall(SNESSolve(snes_alg, NULL, X));

1340:   ctx->stepnum++;

1342:   /* Disturbance period */
1343:   ctx->alg_flg = PETSC_FALSE;
1344:   PetscCall(TSSetTime(ts, ctx->tfaulton));
1345:   PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
1346:   PetscCall(TSSolve(ts, X));
1347:   PetscCall(TSGetStepNumber(ts, &steps2));
1348:   steps2 -= steps1;

1350:   /* Remove the fault */
1351:   row_loc = 2 * ctx->faultbus;
1352:   col_loc = 2 * ctx->faultbus + 1;
1353:   val     = -1 / ctx->Rfault;
1354:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1355:   row_loc = 2 * ctx->faultbus + 1;
1356:   col_loc = 2 * ctx->faultbus;
1357:   val     = -1 / ctx->Rfault;
1358:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

1360:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1361:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

1363:   PetscCall(MatZeroEntries(ctx->J));

1365:   ctx->alg_flg = PETSC_TRUE;

1367:   /* Solve the algebraic equations */
1368:   PetscCall(SNESSolve(snes_alg, NULL, X));

1370:   ctx->stepnum++;

1372:   /* Post-disturbance period */
1373:   ctx->alg_flg = PETSC_TRUE;
1374:   PetscCall(TSSetTime(ts, ctx->tfaultoff));
1375:   PetscCall(TSSetMaxTime(ts, ctx->tmax));
1376:   PetscCall(TSSolve(ts, X));
1377:   PetscCall(TSGetStepNumber(ts, &steps3));
1378:   steps3 -= steps2;
1379:   steps3 -= steps1;

1381:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1382:      Adjoint model starts here
1383:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1384:   PetscCall(TSSetPostStep(ts, NULL));
1385:   PetscCall(MatCreateVecs(ctx->J, &lambda[0], NULL));
1386:   /*   Set initial conditions for the adjoint integration */
1387:   PetscCall(VecZeroEntries(lambda[0]));

1389:   PetscCall(MatCreateVecs(ctx->Jacp, &mu[0], NULL));
1390:   PetscCall(VecZeroEntries(mu[0]));
1391:   PetscCall(TSSetCostGradients(ts, 1, lambda, mu));

1393:   PetscCall(TSAdjointSetSteps(ts, steps3));
1394:   PetscCall(TSAdjointSolve(ts));

1396:   PetscCall(MatZeroEntries(ctx->J));
1397:   /* Applying disturbance - resistive fault at ctx->faultbus */
1398:   /* This is done by deducting shunt conductance to the diagonal location
1399:      in the Ybus matrix */
1400:   row_loc = 2 * ctx->faultbus;
1401:   col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1402:   val     = 1. / ctx->Rfault;
1403:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1404:   row_loc = 2 * ctx->faultbus + 1;
1405:   col_loc = 2 * ctx->faultbus; /* Location for G */
1406:   val     = 1. / ctx->Rfault;
1407:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

1409:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1410:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

1412:   /*   Set number of steps for the adjoint integration */
1413:   PetscCall(TSAdjointSetSteps(ts, steps2));
1414:   PetscCall(TSAdjointSolve(ts));

1416:   PetscCall(MatZeroEntries(ctx->J));
1417:   /* remove the fault */
1418:   row_loc = 2 * ctx->faultbus;
1419:   col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1420:   val     = -1. / ctx->Rfault;
1421:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1422:   row_loc = 2 * ctx->faultbus + 1;
1423:   col_loc = 2 * ctx->faultbus; /* Location for G */
1424:   val     = -1. / ctx->Rfault;
1425:   PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));

1427:   PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1428:   PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));

1430:   /*   Set number of steps for the adjoint integration */
1431:   PetscCall(TSAdjointSetSteps(ts, steps1));
1432:   PetscCall(TSAdjointSolve(ts));

1434:   PetscCall(ComputeSensiP(lambda[0], mu[0], DICDP, ctx));
1435:   PetscCall(VecCopy(mu[0], G));

1437:   PetscCall(TSGetQuadratureTS(ts, NULL, &quadts));
1438:   PetscCall(TSGetSolution(quadts, &q));
1439:   PetscCall(VecGetArray(q, &x_ptr));
1440:   *f       = x_ptr[0];
1441:   x_ptr[0] = 0;
1442:   PetscCall(VecRestoreArray(q, &x_ptr));

1444:   PetscCall(VecDestroy(&lambda[0]));
1445:   PetscCall(VecDestroy(&mu[0]));

1447:   PetscCall(SNESDestroy(&snes_alg));
1448:   PetscCall(VecDestroy(&F_alg));
1449:   PetscCall(VecDestroy(&X));
1450:   PetscCall(TSDestroy(&ts));
1451:   for (i = 0; i < 3; i++) PetscCall(VecDestroy(&DICDP[i]));
1452:   PetscFunctionReturn(PETSC_SUCCESS);
1453: }

1455: /*TEST

1457:    build:
1458:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

1460:    test:
1461:       args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1462:       localrunfiles: petscoptions X.bin Ybus.bin

1464: TEST*/