Actual source code: ex9busoptfd.c

  1: static char help[] = "Using finite difference for the problem in ex9busopt.c \n\n";

  3: /*
  4:   Use finite difference approximations to solve the same optimization problem as in ex9busopt.c.
  5:  */

  7: #include <petsctao.h>
  8: #include <petscts.h>
  9: #include <petscdm.h>
 10: #include <petscdmda.h>
 11: #include <petscdmcomposite.h>

 13: PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);

 15: #define freq 60
 16: #define w_s  (2 * PETSC_PI * freq)

 18: /* Sizes and indices */
 19: const PetscInt nbus    = 9;         /* Number of network buses */
 20: const PetscInt ngen    = 3;         /* Number of generators */
 21: const PetscInt nload   = 3;         /* Number of loads */
 22: const PetscInt gbus[3] = {0, 1, 2}; /* Buses at which generators are incident */
 23: const PetscInt lbus[3] = {4, 5, 7}; /* Buses at which loads are incident */

 25: /* Generator real and reactive powers (found via loadflow) */
 26: PetscScalar PG[3] = {0.69, 1.59, 0.69};
 27: /* PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};*/
 28: const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
 29: /* Generator constants */
 30: const PetscScalar H[3]    = {23.64, 6.4, 3.01};       /* Inertia constant */
 31: const PetscScalar Rs[3]   = {0.0, 0.0, 0.0};          /* Stator Resistance */
 32: const PetscScalar Xd[3]   = {0.146, 0.8958, 1.3125};  /* d-axis reactance */
 33: const PetscScalar Xdp[3]  = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
 34: 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 */
 35: const PetscScalar Xqp[3]  = {0.0969, 0.1969, 0.25};   /* q-axis transient reactance */
 36: const PetscScalar Td0p[3] = {8.96, 6.0, 5.89};        /* d-axis open circuit time constant */
 37: const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6};       /* q-axis open circuit time constant */
 38: PetscScalar       M[3];                               /* M = 2*H/w_s */
 39: PetscScalar       D[3];                               /* D = 0.1*M */

 41: PetscScalar TM[3]; /* Mechanical Torque */
 42: /* Exciter system constants */
 43: const PetscScalar KA[3] = {20.0, 20.0, 20.0};    /* Voltage regulartor gain constant */
 44: const PetscScalar TA[3] = {0.2, 0.2, 0.2};       /* Voltage regulator time constant */
 45: const PetscScalar KE[3] = {1.0, 1.0, 1.0};       /* Exciter gain constant */
 46: const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
 47: const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
 48: const PetscScalar TF[3] = {0.35, 0.35, 0.35};    /* Feedback stabilizer time constant */
 49: const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
 50: const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */

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

 66:     Note: All loads have the same characteristic currently.
 67: */
 68: const PetscScalar PD0[3]       = {1.25, 0.9, 1.0};
 69: const PetscScalar QD0[3]       = {0.5, 0.3, 0.35};
 70: const PetscInt    ld_nsegsp[3] = {3, 3, 3};
 71: const PetscScalar ld_alphap[3] = {1.0, 0.0, 0.0};
 72: const PetscScalar ld_betap[3]  = {2.0, 1.0, 0.0};
 73: const PetscInt    ld_nsegsq[3] = {3, 3, 3};
 74: const PetscScalar ld_alphaq[3] = {1.0, 0.0, 0.0};
 75: const PetscScalar ld_betaq[3]  = {2.0, 1.0, 0.0};

 77: typedef struct {
 78:   DM          dmgen, dmnet;        /* DMs to manage generator and network subsystem */
 79:   DM          dmpgrid;             /* Composite DM to manage the entire power grid */
 80:   Mat         Ybus;                /* Network admittance matrix */
 81:   Vec         V0;                  /* Initial voltage vector (Power flow solution) */
 82:   PetscReal   tfaulton, tfaultoff; /* Fault on and off times */
 83:   PetscInt    faultbus;            /* Fault bus */
 84:   PetscScalar Rfault;
 85:   PetscReal   t0, tmax;
 86:   PetscInt    neqs_gen, neqs_net, neqs_pgrid;
 87:   Mat         Sol; /* Matrix to save solution at each time step */
 88:   PetscInt    stepnum;
 89:   PetscBool   alg_flg;
 90:   PetscReal   t;
 91:   IS          is_diff;        /* indices for differential equations */
 92:   IS          is_alg;         /* indices for algebraic equations */
 93:   PetscReal   freq_u, freq_l; /* upper and lower frequency limit */
 94:   PetscInt    pow;            /* power coefficient used in the cost function */
 95:   Vec         vec_q;
 96: } Userctx;

 98: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
 99: PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
100: {
101:   *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
102:   *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
103:   return 0;
104: }

106: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
107: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
108: {
109:   *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
110:   *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
111:   return 0;
112: }

114: PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
115: {
116:   Vec          Xgen, Xnet;
117:   PetscScalar *xgen, *xnet;
118:   PetscInt     i, idx = 0;
119:   PetscScalar  Vr, Vi, IGr, IGi, Vm, Vm2;
120:   PetscScalar  Eqp, Edp, delta;
121:   PetscScalar  Efd, RF, VR; /* Exciter variables */
122:   PetscScalar  Id, Iq;      /* Generator dq axis currents */
123:   PetscScalar  theta, Vd, Vq, SE;

125:   M[0] = 2 * H[0] / w_s;
126:   M[1] = 2 * H[1] / w_s;
127:   M[2] = 2 * H[2] / w_s;
128:   D[0] = 0.1 * M[0];
129:   D[1] = 0.1 * M[1];
130:   D[2] = 0.1 * M[2];

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

134:   /* Network subsystem initialization */
135:   VecCopy(user->V0, Xnet);

137:   /* Generator subsystem initialization */
138:   VecGetArray(Xgen, &xgen);
139:   VecGetArray(Xnet, &xnet);

141:   for (i = 0; i < ngen; i++) {
142:     Vr  = xnet[2 * gbus[i]];     /* Real part of generator terminal voltage */
143:     Vi  = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
144:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
145:     Vm2 = Vm * Vm;
146:     IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
147:     IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;

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

151:     theta = PETSC_PI / 2.0 - delta;

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

156:     Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
157:     Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);

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

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

164:     /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
165:     xgen[idx]     = Eqp;
166:     xgen[idx + 1] = Edp;
167:     xgen[idx + 2] = delta;
168:     xgen[idx + 3] = w_s;

170:     idx = idx + 4;

172:     xgen[idx]     = Id;
173:     xgen[idx + 1] = Iq;

175:     idx = idx + 2;

177:     /* Exciter */
178:     Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
179:     SE  = k1[i] * PetscExpScalar(k2[i] * Efd);
180:     VR  = KE[i] * Efd + SE;
181:     RF  = KF[i] * Efd / TF[i];

183:     xgen[idx]     = Efd;
184:     xgen[idx + 1] = RF;
185:     xgen[idx + 2] = VR;

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

189:     idx = idx + 3;
190:   }

192:   VecRestoreArray(Xgen, &xgen);
193:   VecRestoreArray(Xnet, &xnet);

195:   /* VecView(Xgen,0); */
196:   DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet);
197:   DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
198:   return 0;
199: }

201: /* Computes F = [-f(x,y);g(x,y)] */
202: PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
203: {
204:   Vec          Xgen, Xnet, Fgen, Fnet;
205:   PetscScalar *xgen, *xnet, *fgen, *fnet;
206:   PetscInt     i, idx = 0;
207:   PetscScalar  Vr, Vi, Vm, Vm2;
208:   PetscScalar  Eqp, Edp, delta, w; /* Generator variables */
209:   PetscScalar  Efd, RF, VR;        /* Exciter variables */
210:   PetscScalar  Id, Iq;             /* Generator dq axis currents */
211:   PetscScalar  Vd, Vq, SE;
212:   PetscScalar  IGr, IGi, IDr, IDi;
213:   PetscScalar  Zdq_inv[4], det;
214:   PetscScalar  PD, QD, Vm0, *v0;
215:   PetscInt     k;

217:   VecZeroEntries(F);
218:   DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
219:   DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet);
220:   DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);
221:   DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet);

223:   /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
224:      The generator current injection, IG, and load current injection, ID are added later
225:   */
226:   /* Note that the values in Ybus are stored assuming the imaginary current balance
227:      equation is ordered first followed by real current balance equation for each bus.
228:      Thus imaginary current contribution goes in location 2*i, and
229:      real current contribution in 2*i+1
230:   */
231:   MatMult(user->Ybus, Xnet, Fnet);

233:   VecGetArray(Xgen, &xgen);
234:   VecGetArray(Xnet, &xnet);
235:   VecGetArray(Fgen, &fgen);
236:   VecGetArray(Fnet, &fnet);

238:   /* Generator subsystem */
239:   for (i = 0; i < ngen; i++) {
240:     Eqp   = xgen[idx];
241:     Edp   = xgen[idx + 1];
242:     delta = xgen[idx + 2];
243:     w     = xgen[idx + 3];
244:     Id    = xgen[idx + 4];
245:     Iq    = xgen[idx + 5];
246:     Efd   = xgen[idx + 6];
247:     RF    = xgen[idx + 7];
248:     VR    = xgen[idx + 8];

250:     /* Generator differential equations */
251:     fgen[idx]     = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
252:     fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
253:     fgen[idx + 2] = -w + w_s;
254:     fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];

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

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

263:     Zdq_inv[0] = Rs[i] / det;
264:     Zdq_inv[1] = Xqp[i] / det;
265:     Zdq_inv[2] = -Xdp[i] / det;
266:     Zdq_inv[3] = Rs[i] / det;

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

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

274:     fnet[2 * gbus[i]] -= IGi;
275:     fnet[2 * gbus[i] + 1] -= IGr;

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

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

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

286:     idx = idx + 9;
287:   }

289:   VecGetArray(user->V0, &v0);
290:   for (i = 0; i < nload; i++) {
291:     Vr  = xnet[2 * lbus[i]];     /* Real part of load bus voltage */
292:     Vi  = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
293:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
294:     Vm2 = Vm * Vm;
295:     Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
296:     PD = QD = 0.0;
297:     for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
298:     for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);

300:     /* Load currents */
301:     IDr = (PD * Vr + QD * Vi) / Vm2;
302:     IDi = (-QD * Vr + PD * Vi) / Vm2;

304:     fnet[2 * lbus[i]] += IDi;
305:     fnet[2 * lbus[i] + 1] += IDr;
306:   }
307:   VecRestoreArray(user->V0, &v0);

309:   VecRestoreArray(Xgen, &xgen);
310:   VecRestoreArray(Xnet, &xnet);
311:   VecRestoreArray(Fgen, &fgen);
312:   VecRestoreArray(Fnet, &fnet);

314:   DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet);
315:   DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
316:   DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet);
317:   return 0;
318: }

320: /* \dot{x} - f(x,y)
321:      g(x,y) = 0
322:  */
323: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
324: {
325:   SNES               snes;
326:   PetscScalar       *f;
327:   const PetscScalar *xdot;
328:   PetscInt           i;

330:   user->t = t;

332:   TSGetSNES(ts, &snes);
333:   ResidualFunction(snes, X, F, user);
334:   VecGetArray(F, &f);
335:   VecGetArrayRead(Xdot, &xdot);
336:   for (i = 0; i < ngen; i++) {
337:     f[9 * i] += xdot[9 * i];
338:     f[9 * i + 1] += xdot[9 * i + 1];
339:     f[9 * i + 2] += xdot[9 * i + 2];
340:     f[9 * i + 3] += xdot[9 * i + 3];
341:     f[9 * i + 6] += xdot[9 * i + 6];
342:     f[9 * i + 7] += xdot[9 * i + 7];
343:     f[9 * i + 8] += xdot[9 * i + 8];
344:   }
345:   VecRestoreArray(F, &f);
346:   VecRestoreArrayRead(Xdot, &xdot);
347:   return 0;
348: }

350: /* This function is used for solving the algebraic system only during fault on and
351:    off times. It computes the entire F and then zeros out the part corresponding to
352:    differential equations
353:  F = [0;g(y)];
354: */
355: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
356: {
357:   Userctx     *user = (Userctx *)ctx;
358:   PetscScalar *f;
359:   PetscInt     i;

361:   ResidualFunction(snes, X, F, user);
362:   VecGetArray(F, &f);
363:   for (i = 0; i < ngen; i++) {
364:     f[9 * i]     = 0;
365:     f[9 * i + 1] = 0;
366:     f[9 * i + 2] = 0;
367:     f[9 * i + 3] = 0;
368:     f[9 * i + 6] = 0;
369:     f[9 * i + 7] = 0;
370:     f[9 * i + 8] = 0;
371:   }
372:   VecRestoreArray(F, &f);
373:   return 0;
374: }

376: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
377: {
378:   PetscInt *d_nnz;
379:   PetscInt  i, idx = 0, start = 0;
380:   PetscInt  ncols;

382:   PetscMalloc1(user->neqs_pgrid, &d_nnz);
383:   for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
384:   /* Generator subsystem */
385:   for (i = 0; i < ngen; i++) {
386:     d_nnz[idx] += 3;
387:     d_nnz[idx + 1] += 2;
388:     d_nnz[idx + 2] += 2;
389:     d_nnz[idx + 3] += 5;
390:     d_nnz[idx + 4] += 6;
391:     d_nnz[idx + 5] += 6;

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

396:     d_nnz[idx + 6] += 2;
397:     d_nnz[idx + 7] += 2;
398:     d_nnz[idx + 8] += 5;

400:     idx = idx + 9;
401:   }

403:   start = user->neqs_gen;

405:   for (i = 0; i < nbus; i++) {
406:     MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL);
407:     d_nnz[start + 2 * i] += ncols;
408:     d_nnz[start + 2 * i + 1] += ncols;
409:     MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL);
410:   }

412:   MatSeqAIJSetPreallocation(J, 0, d_nnz);

414:   PetscFree(d_nnz);
415:   return 0;
416: }

418: /*
419:    J = [-df_dx, -df_dy
420:         dg_dx, dg_dy]
421: */
422: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, void *ctx)
423: {
424:   Userctx           *user = (Userctx *)ctx;
425:   Vec                Xgen, Xnet;
426:   PetscScalar       *xgen, *xnet;
427:   PetscInt           i, idx = 0;
428:   PetscScalar        Vr, Vi, Vm, Vm2;
429:   PetscScalar        Eqp, Edp, delta; /* Generator variables */
430:   PetscScalar        Efd;             /* Exciter variables */
431:   PetscScalar        Id, Iq;          /* Generator dq axis currents */
432:   PetscScalar        Vd, Vq;
433:   PetscScalar        val[10];
434:   PetscInt           row[2], col[10];
435:   PetscInt           net_start = user->neqs_gen;
436:   PetscScalar        Zdq_inv[4], det;
437:   PetscScalar        dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
438:   PetscScalar        dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
439:   PetscScalar        dSE_dEfd;
440:   PetscScalar        dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
441:   PetscInt           ncols;
442:   const PetscInt    *cols;
443:   const PetscScalar *yvals;
444:   PetscInt           k;
445:   PetscScalar        PD, QD, Vm0, *v0, Vm4;
446:   PetscScalar        dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
447:   PetscScalar        dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;

449:   MatZeroEntries(B);
450:   DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
451:   DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);

453:   VecGetArray(Xgen, &xgen);
454:   VecGetArray(Xnet, &xnet);

456:   /* Generator subsystem */
457:   for (i = 0; i < ngen; i++) {
458:     Eqp   = xgen[idx];
459:     Edp   = xgen[idx + 1];
460:     delta = xgen[idx + 2];
461:     Id    = xgen[idx + 4];
462:     Iq    = xgen[idx + 5];
463:     Efd   = xgen[idx + 6];

465:     /*    fgen[idx]   = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
466:     row[0] = idx;
467:     col[0] = idx;
468:     col[1] = idx + 4;
469:     col[2] = idx + 6;
470:     val[0] = 1 / Td0p[i];
471:     val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
472:     val[2] = -1 / Td0p[i];

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

476:     /*    fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
477:     row[0] = idx + 1;
478:     col[0] = idx + 1;
479:     col[1] = idx + 5;
480:     val[0] = 1 / Tq0p[i];
481:     val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
482:     MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);

484:     /*    fgen[idx+2] = - w + w_s; */
485:     row[0] = idx + 2;
486:     col[0] = idx + 2;
487:     col[1] = idx + 3;
488:     val[0] = 0;
489:     val[1] = -1;
490:     MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);

492:     /*    fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
493:     row[0] = idx + 3;
494:     col[0] = idx;
495:     col[1] = idx + 1;
496:     col[2] = idx + 3;
497:     col[3] = idx + 4;
498:     col[4] = idx + 5;
499:     val[0] = Iq / M[i];
500:     val[1] = Id / M[i];
501:     val[2] = D[i] / M[i];
502:     val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
503:     val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
504:     MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES);

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

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

512:     Zdq_inv[0] = Rs[i] / det;
513:     Zdq_inv[1] = Xqp[i] / det;
514:     Zdq_inv[2] = -Xdp[i] / det;
515:     Zdq_inv[3] = Rs[i] / det;

517:     dVd_dVr    = PetscSinScalar(delta);
518:     dVd_dVi    = -PetscCosScalar(delta);
519:     dVq_dVr    = PetscCosScalar(delta);
520:     dVq_dVi    = PetscSinScalar(delta);
521:     dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
522:     dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);

524:     /*    fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
525:     row[0] = idx + 4;
526:     col[0] = idx;
527:     col[1] = idx + 1;
528:     col[2] = idx + 2;
529:     val[0] = -Zdq_inv[1];
530:     val[1] = -Zdq_inv[0];
531:     val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
532:     col[3] = idx + 4;
533:     col[4] = net_start + 2 * gbus[i];
534:     col[5] = net_start + 2 * gbus[i] + 1;
535:     val[3] = 1;
536:     val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
537:     val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
538:     MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES);

540:     /*  fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
541:     row[0] = idx + 5;
542:     col[0] = idx;
543:     col[1] = idx + 1;
544:     col[2] = idx + 2;
545:     val[0] = -Zdq_inv[3];
546:     val[1] = -Zdq_inv[2];
547:     val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
548:     col[3] = idx + 5;
549:     col[4] = net_start + 2 * gbus[i];
550:     col[5] = net_start + 2 * gbus[i] + 1;
551:     val[3] = 1;
552:     val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
553:     val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
554:     MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES);

556:     dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
557:     dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
558:     dIGr_dId    = PetscSinScalar(delta);
559:     dIGr_dIq    = PetscCosScalar(delta);
560:     dIGi_dId    = -PetscCosScalar(delta);
561:     dIGi_dIq    = PetscSinScalar(delta);

563:     /* fnet[2*gbus[i]]   -= IGi; */
564:     row[0] = net_start + 2 * gbus[i];
565:     col[0] = idx + 2;
566:     col[1] = idx + 4;
567:     col[2] = idx + 5;
568:     val[0] = -dIGi_ddelta;
569:     val[1] = -dIGi_dId;
570:     val[2] = -dIGi_dIq;
571:     MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);

573:     /* fnet[2*gbus[i]+1]   -= IGr; */
574:     row[0] = net_start + 2 * gbus[i] + 1;
575:     col[0] = idx + 2;
576:     col[1] = idx + 4;
577:     col[2] = idx + 5;
578:     val[0] = -dIGr_ddelta;
579:     val[1] = -dIGr_dId;
580:     val[2] = -dIGr_dIq;
581:     MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);

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

585:     /*    fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
586:     /*    SE  = k1[i]*PetscExpScalar(k2[i]*Efd); */

588:     dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);

590:     row[0] = idx + 6;
591:     col[0] = idx + 6;
592:     col[1] = idx + 8;
593:     val[0] = (KE[i] + dSE_dEfd) / TE[i];
594:     val[1] = -1 / TE[i];
595:     MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);

597:     /* Exciter differential equations */

599:     /*    fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
600:     row[0] = idx + 7;
601:     col[0] = idx + 6;
602:     col[1] = idx + 7;
603:     val[0] = (-KF[i] / TF[i]) / TF[i];
604:     val[1] = 1 / TF[i];
605:     MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);

607:     /*    fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
608:     /* Vm = (Vd^2 + Vq^2)^0.5; */
609:     dVm_dVd = Vd / Vm;
610:     dVm_dVq = Vq / Vm;
611:     dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
612:     dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
613:     row[0]  = idx + 8;
614:     col[0]  = idx + 6;
615:     col[1]  = idx + 7;
616:     col[2]  = idx + 8;
617:     val[0]  = (KA[i] * KF[i] / TF[i]) / TA[i];
618:     val[1]  = -KA[i] / TA[i];
619:     val[2]  = 1 / TA[i];
620:     col[3]  = net_start + 2 * gbus[i];
621:     col[4]  = net_start + 2 * gbus[i] + 1;
622:     val[3]  = KA[i] * dVm_dVr / TA[i];
623:     val[4]  = KA[i] * dVm_dVi / TA[i];
624:     MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES);
625:     idx = idx + 9;
626:   }

628:   for (i = 0; i < nbus; i++) {
629:     MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals);
630:     row[0] = net_start + 2 * i;
631:     for (k = 0; k < ncols; k++) {
632:       col[k] = net_start + cols[k];
633:       val[k] = yvals[k];
634:     }
635:     MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES);
636:     MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals);

638:     MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals);
639:     row[0] = net_start + 2 * i + 1;
640:     for (k = 0; k < ncols; k++) {
641:       col[k] = net_start + cols[k];
642:       val[k] = yvals[k];
643:     }
644:     MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES);
645:     MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals);
646:   }

648:   MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY);
649:   MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY);

651:   VecGetArray(user->V0, &v0);
652:   for (i = 0; i < nload; i++) {
653:     Vr  = xnet[2 * lbus[i]];     /* Real part of load bus voltage */
654:     Vi  = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
655:     Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
656:     Vm2 = Vm * Vm;
657:     Vm4 = Vm2 * Vm2;
658:     Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
659:     PD = QD = 0.0;
660:     dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
661:     for (k = 0; k < ld_nsegsp[i]; k++) {
662:       PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
663:       dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vr * PetscPowScalar(Vm, (ld_betap[k] - 2));
664:       dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vi * PetscPowScalar(Vm, (ld_betap[k] - 2));
665:     }
666:     for (k = 0; k < ld_nsegsq[i]; k++) {
667:       QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
668:       dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vr * PetscPowScalar(Vm, (ld_betaq[k] - 2));
669:       dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vi * PetscPowScalar(Vm, (ld_betaq[k] - 2));
670:     }

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

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

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

681:     /*    fnet[2*lbus[i]]   += IDi; */
682:     row[0] = net_start + 2 * lbus[i];
683:     col[0] = net_start + 2 * lbus[i];
684:     col[1] = net_start + 2 * lbus[i] + 1;
685:     val[0] = dIDi_dVr;
686:     val[1] = dIDi_dVi;
687:     MatSetValues(J, 1, row, 2, col, val, ADD_VALUES);
688:     /*    fnet[2*lbus[i]+1] += IDr; */
689:     row[0] = net_start + 2 * lbus[i] + 1;
690:     col[0] = net_start + 2 * lbus[i];
691:     col[1] = net_start + 2 * lbus[i] + 1;
692:     val[0] = dIDr_dVr;
693:     val[1] = dIDr_dVi;
694:     MatSetValues(J, 1, row, 2, col, val, ADD_VALUES);
695:   }
696:   VecRestoreArray(user->V0, &v0);

698:   VecRestoreArray(Xgen, &xgen);
699:   VecRestoreArray(Xnet, &xnet);

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

703:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
704:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
705:   return 0;
706: }

708: /*
709:    J = [I, 0
710:         dg_dx, dg_dy]
711: */
712: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, void *ctx)
713: {
714:   Userctx *user = (Userctx *)ctx;

716:   ResidualJacobian(snes, X, A, B, ctx);
717:   MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
718:   MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL);
719:   return 0;
720: }

722: /*
723:    J = [a*I-df_dx, -df_dy
724:         dg_dx, dg_dy]
725: */

727: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
728: {
729:   SNES        snes;
730:   PetscScalar atmp = (PetscScalar)a;
731:   PetscInt    i, row;

733:   user->t = t;

735:   TSGetSNES(ts, &snes);
736:   ResidualJacobian(snes, X, A, B, user);
737:   for (i = 0; i < ngen; i++) {
738:     row = 9 * i;
739:     MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
740:     row = 9 * i + 1;
741:     MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
742:     row = 9 * i + 2;
743:     MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
744:     row = 9 * i + 3;
745:     MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
746:     row = 9 * i + 6;
747:     MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
748:     row = 9 * i + 7;
749:     MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
750:     row = 9 * i + 8;
751:     MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
752:   }
753:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
754:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
755:   return 0;
756: }

758: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
759: {
760:   PetscScalar       *r;
761:   const PetscScalar *u;
762:   PetscInt           idx;
763:   Vec                Xgen, Xnet;
764:   PetscScalar       *xgen;
765:   PetscInt           i;

767:   DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
768:   DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet);

770:   VecGetArray(Xgen, &xgen);

772:   VecGetArrayRead(U, &u);
773:   VecGetArray(R, &r);
774:   r[0] = 0.;

776:   idx = 0;
777:   for (i = 0; i < ngen; i++) {
778:     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);
779:     idx += 9;
780:   }
781:   VecRestoreArray(R, &r);
782:   VecRestoreArrayRead(U, &u);
783:   DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
784:   return 0;
785: }

787: static PetscErrorCode MonitorUpdateQ(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx0)
788: {
789:   Vec       C, *Y;
790:   PetscInt  Nr;
791:   PetscReal h, theta;
792:   Userctx  *ctx = (Userctx *)ctx0;

794:   theta = 0.5;
795:   TSGetStages(ts, &Nr, &Y);
796:   TSGetTimeStep(ts, &h);
797:   VecDuplicate(ctx->vec_q, &C);
798:   /* compute integrals */
799:   if (stepnum > 0) {
800:     CostIntegrand(ts, time, X, C, ctx);
801:     VecAXPY(ctx->vec_q, h * theta, C);
802:     CostIntegrand(ts, time + h * theta, Y[0], C, ctx);
803:     VecAXPY(ctx->vec_q, h * (1 - theta), C);
804:   }
805:   VecDestroy(&C);
806:   return 0;
807: }

809: int main(int argc, char **argv)
810: {
811:   Userctx      user;
812:   Vec          p;
813:   PetscScalar *x_ptr;
814:   PetscMPIInt  size;
815:   PetscInt     i;
816:   KSP          ksp;
817:   PC           pc;
818:   PetscInt    *idx2;
819:   Tao          tao;
820:   Vec          lowerb, upperb;

824:   PetscInitialize(&argc, &argv, "petscoptions", help);
825:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

828:   VecCreateSeq(PETSC_COMM_WORLD, 1, &user.vec_q);

830:   user.neqs_gen   = 9 * ngen; /* # eqs. for generator subsystem */
831:   user.neqs_net   = 2 * nbus; /* # eqs. for network subsystem   */
832:   user.neqs_pgrid = user.neqs_gen + user.neqs_net;

834:   /* Create indices for differential and algebraic equations */
835:   PetscMalloc1(7 * ngen, &idx2);
836:   for (i = 0; i < ngen; i++) {
837:     idx2[7 * i]     = 9 * i;
838:     idx2[7 * i + 1] = 9 * i + 1;
839:     idx2[7 * i + 2] = 9 * i + 2;
840:     idx2[7 * i + 3] = 9 * i + 3;
841:     idx2[7 * i + 4] = 9 * i + 6;
842:     idx2[7 * i + 5] = 9 * i + 7;
843:     idx2[7 * i + 6] = 9 * i + 8;
844:   }
845:   ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff);
846:   ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg);
847:   PetscFree(idx2);

849:   /* Set run time options */
850:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
851:   {
852:     user.tfaulton  = 1.0;
853:     user.tfaultoff = 1.2;
854:     user.Rfault    = 0.0001;
855:     user.faultbus  = 8;
856:     PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL);
857:     PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL);
858:     PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL);
859:     user.t0   = 0.0;
860:     user.tmax = 1.5;
861:     PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL);
862:     PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL);
863:     user.freq_u = 61.0;
864:     user.freq_l = 59.0;
865:     user.pow    = 2;
866:     PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL);
867:     PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL);
868:     PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL);
869:   }
870:   PetscOptionsEnd();

872:   /* Create DMs for generator and network subsystems */
873:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen);
874:   DMSetOptionsPrefix(user.dmgen, "dmgen_");
875:   DMSetFromOptions(user.dmgen);
876:   DMSetUp(user.dmgen);
877:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet);
878:   DMSetOptionsPrefix(user.dmnet, "dmnet_");
879:   DMSetFromOptions(user.dmnet);
880:   DMSetUp(user.dmnet);
881:   /* Create a composite DM packer and add the two DMs */
882:   DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid);
883:   DMSetOptionsPrefix(user.dmpgrid, "pgrid_");
884:   DMCompositeAddDM(user.dmpgrid, user.dmgen);
885:   DMCompositeAddDM(user.dmpgrid, user.dmnet);

887:   /* Create TAO solver and set desired solution method */
888:   TaoCreate(PETSC_COMM_WORLD, &tao);
889:   TaoSetType(tao, TAOBLMVM);
890:   /*
891:      Optimization starts
892:   */
893:   /* Set initial solution guess */
894:   VecCreateSeq(PETSC_COMM_WORLD, 3, &p);
895:   VecGetArray(p, &x_ptr);
896:   x_ptr[0] = PG[0];
897:   x_ptr[1] = PG[1];
898:   x_ptr[2] = PG[2];
899:   VecRestoreArray(p, &x_ptr);

901:   TaoSetSolution(tao, p);
902:   /* Set routine for function and gradient evaluation */
903:   TaoSetObjective(tao, FormFunction, (void *)&user);
904:   TaoSetGradient(tao, NULL, TaoDefaultComputeGradient, (void *)&user);

906:   /* Set bounds for the optimization */
907:   VecDuplicate(p, &lowerb);
908:   VecDuplicate(p, &upperb);
909:   VecGetArray(lowerb, &x_ptr);
910:   x_ptr[0] = 0.5;
911:   x_ptr[1] = 0.5;
912:   x_ptr[2] = 0.5;
913:   VecRestoreArray(lowerb, &x_ptr);
914:   VecGetArray(upperb, &x_ptr);
915:   x_ptr[0] = 2.0;
916:   x_ptr[1] = 2.0;
917:   x_ptr[2] = 2.0;
918:   VecRestoreArray(upperb, &x_ptr);
919:   TaoSetVariableBounds(tao, lowerb, upperb);

921:   /* Check for any TAO command line options */
922:   TaoSetFromOptions(tao);
923:   TaoGetKSP(tao, &ksp);
924:   if (ksp) {
925:     KSPGetPC(ksp, &pc);
926:     PCSetType(pc, PCNONE);
927:   }

929:   /* SOLVE THE APPLICATION */
930:   TaoSolve(tao);

932:   VecView(p, PETSC_VIEWER_STDOUT_WORLD);
933:   /* Free TAO data structures */
934:   TaoDestroy(&tao);
935:   VecDestroy(&user.vec_q);
936:   VecDestroy(&lowerb);
937:   VecDestroy(&upperb);
938:   VecDestroy(&p);
939:   DMDestroy(&user.dmgen);
940:   DMDestroy(&user.dmnet);
941:   DMDestroy(&user.dmpgrid);
942:   ISDestroy(&user.is_diff);
943:   ISDestroy(&user.is_alg);
944:   PetscFinalize();
945:   return 0;
946: }

948: /* ------------------------------------------------------------------ */
949: /*
950:    FormFunction - Evaluates the function and corresponding gradient.

952:    Input Parameters:
953:    tao - the Tao context
954:    X   - the input vector
955:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

957:    Output Parameters:
958:    f   - the newly evaluated function
959: */
960: PetscErrorCode FormFunction(Tao tao, Vec P, PetscReal *f, void *ctx0)
961: {
962:   TS       ts;
963:   SNES     snes_alg;
964:   Userctx *ctx = (Userctx *)ctx0;
965:   Vec      X;
966:   Mat      J;
967:   /* sensitivity context */
968:   PetscScalar *x_ptr;
969:   PetscViewer  Xview, Ybusview;
970:   Vec          F_alg;
971:   Vec          Xdot;
972:   PetscInt     row_loc, col_loc;
973:   PetscScalar  val;

975:   VecGetArrayRead(P, (const PetscScalar **)&x_ptr);
976:   PG[0] = x_ptr[0];
977:   PG[1] = x_ptr[1];
978:   PG[2] = x_ptr[2];
979:   VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr);

981:   ctx->stepnum = 0;

983:   VecZeroEntries(ctx->vec_q);

985:   /* Read initial voltage vector and Ybus */
986:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview);
987:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview);

989:   VecCreate(PETSC_COMM_WORLD, &ctx->V0);
990:   VecSetSizes(ctx->V0, PETSC_DECIDE, ctx->neqs_net);
991:   VecLoad(ctx->V0, Xview);

993:   MatCreate(PETSC_COMM_WORLD, &ctx->Ybus);
994:   MatSetSizes(ctx->Ybus, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_net, ctx->neqs_net);
995:   MatSetType(ctx->Ybus, MATBAIJ);
996:   /*  MatSetBlockSize(ctx->Ybus,2); */
997:   MatLoad(ctx->Ybus, Ybusview);

999:   PetscViewerDestroy(&Xview);
1000:   PetscViewerDestroy(&Ybusview);

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

1004:   MatCreate(PETSC_COMM_WORLD, &J);
1005:   MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_pgrid, ctx->neqs_pgrid);
1006:   MatSetFromOptions(J);
1007:   PreallocateJacobian(J, ctx);

1009:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1010:      Create timestepping solver context
1011:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1012:   TSCreate(PETSC_COMM_WORLD, &ts);
1013:   TSSetProblemType(ts, TS_NONLINEAR);
1014:   TSSetType(ts, TSCN);
1015:   TSSetIFunction(ts, NULL, (TSIFunction)IFunction, ctx);
1016:   TSSetIJacobian(ts, J, J, (TSIJacobian)IJacobian, ctx);
1017:   TSSetApplicationContext(ts, ctx);

1019:   TSMonitorSet(ts, MonitorUpdateQ, ctx, NULL);
1020:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1021:      Set initial conditions
1022:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1023:   SetInitialGuess(X, ctx);

1025:   VecDuplicate(X, &F_alg);
1026:   SNESCreate(PETSC_COMM_WORLD, &snes_alg);
1027:   SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx);
1028:   MatZeroEntries(J);
1029:   SNESSetJacobian(snes_alg, J, J, AlgJacobian, ctx);
1030:   SNESSetOptionsPrefix(snes_alg, "alg_");
1031:   SNESSetFromOptions(snes_alg);
1032:   ctx->alg_flg = PETSC_TRUE;
1033:   /* Solve the algebraic equations */
1034:   SNESSolve(snes_alg, NULL, X);

1036:   /* Just to set up the Jacobian structure */
1037:   VecDuplicate(X, &Xdot);
1038:   IJacobian(ts, 0.0, X, Xdot, 0.0, J, J, ctx);
1039:   VecDestroy(&Xdot);

1041:   ctx->stepnum++;

1043:   TSSetTimeStep(ts, 0.01);
1044:   TSSetMaxTime(ts, ctx->tmax);
1045:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1046:   TSSetFromOptions(ts);

1048:   /* Prefault period */
1049:   ctx->alg_flg = PETSC_FALSE;
1050:   TSSetTime(ts, 0.0);
1051:   TSSetMaxTime(ts, ctx->tfaulton);
1052:   TSSolve(ts, X);

1054:   /* Create the nonlinear solver for solving the algebraic system */
1055:   /* Note that although the algebraic system needs to be solved only for
1056:      Idq and V, we reuse the entire system including xgen. The xgen
1057:      variables are held constant by setting their residuals to 0 and
1058:      putting a 1 on the Jacobian diagonal for xgen rows
1059:   */
1060:   MatZeroEntries(J);

1062:   /* Apply disturbance - resistive fault at ctx->faultbus */
1063:   /* This is done by adding shunt conductance to the diagonal location
1064:      in the Ybus matrix */
1065:   row_loc = 2 * ctx->faultbus;
1066:   col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1067:   val     = 1 / ctx->Rfault;
1068:   MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1069:   row_loc = 2 * ctx->faultbus + 1;
1070:   col_loc = 2 * ctx->faultbus; /* Location for G */
1071:   val     = 1 / ctx->Rfault;
1072:   MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);

1074:   MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1075:   MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY);

1077:   ctx->alg_flg = PETSC_TRUE;
1078:   /* Solve the algebraic equations */
1079:   SNESSolve(snes_alg, NULL, X);

1081:   ctx->stepnum++;

1083:   /* Disturbance period */
1084:   ctx->alg_flg = PETSC_FALSE;
1085:   TSSetTime(ts, ctx->tfaulton);
1086:   TSSetMaxTime(ts, ctx->tfaultoff);
1087:   TSSolve(ts, X);

1089:   /* Remove the fault */
1090:   row_loc = 2 * ctx->faultbus;
1091:   col_loc = 2 * ctx->faultbus + 1;
1092:   val     = -1 / ctx->Rfault;
1093:   MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
1094:   row_loc = 2 * ctx->faultbus + 1;
1095:   col_loc = 2 * ctx->faultbus;
1096:   val     = -1 / ctx->Rfault;
1097:   MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);

1099:   MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY);
1100:   MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY);

1102:   MatZeroEntries(J);

1104:   ctx->alg_flg = PETSC_TRUE;

1106:   /* Solve the algebraic equations */
1107:   SNESSolve(snes_alg, NULL, X);

1109:   ctx->stepnum++;

1111:   /* Post-disturbance period */
1112:   ctx->alg_flg = PETSC_TRUE;
1113:   TSSetTime(ts, ctx->tfaultoff);
1114:   TSSetMaxTime(ts, ctx->tmax);
1115:   TSSolve(ts, X);

1117:   VecGetArray(ctx->vec_q, &x_ptr);
1118:   *f = x_ptr[0];
1119:   VecRestoreArray(ctx->vec_q, &x_ptr);

1121:   MatDestroy(&ctx->Ybus);
1122:   VecDestroy(&ctx->V0);
1123:   SNESDestroy(&snes_alg);
1124:   VecDestroy(&F_alg);
1125:   MatDestroy(&J);
1126:   VecDestroy(&X);
1127:   TSDestroy(&ts);

1129:   return 0;
1130: }

1132: /*TEST

1134:   build:
1135:       requires: double !complex !defined(USE_64BIT_INDICES)

1137:    test:
1138:       args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1139:       localrunfiles: petscoptions X.bin Ybus.bin

1141: TEST*/