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; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s;
126:   D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2];

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

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

133:   /* Generator subsystem initialization */
134:   VecGetArray(Xgen,&xgen);
135:   VecGetArray(Xnet,&xnet);

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

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

146:     theta = PETSC_PI/2.0 - delta;

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

151:     Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta);
152:     Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta);

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

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

159:     /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
160:     xgen[idx]   = Eqp;
161:     xgen[idx+1] = Edp;
162:     xgen[idx+2] = delta;
163:     xgen[idx+3] = w_s;

165:     idx = idx + 4;

167:     xgen[idx]   = Id;
168:     xgen[idx+1] = Iq;

170:     idx = idx + 2;

172:     /* Exciter */
173:     Efd = Eqp + (Xd[i] - Xdp[i])*Id;
174:     SE  = k1[i]*PetscExpScalar(k2[i]*Efd);
175:     VR  =  KE[i]*Efd + SE;
176:     RF  =  KF[i]*Efd/TF[i];

178:     xgen[idx]   = Efd;
179:     xgen[idx+1] = RF;
180:     xgen[idx+2] = VR;

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

184:     idx = idx + 3;
185:   }

187:   VecRestoreArray(Xgen,&xgen);
188:   VecRestoreArray(Xnet,&xnet);

190:   /* VecView(Xgen,0); */
191:   DMCompositeGather(user->dmpgrid,INSERT_VALUES,X,Xgen,Xnet);
192:   DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);
193:   return 0;
194: }

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

212:   VecZeroEntries(F);
213:   DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);
214:   DMCompositeGetLocalVectors(user->dmpgrid,&Fgen,&Fnet);
215:   DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);
216:   DMCompositeScatter(user->dmpgrid,F,Fgen,Fnet);

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

228:   VecGetArray(Xgen,&xgen);
229:   VecGetArray(Xnet,&xnet);
230:   VecGetArray(Fgen,&fgen);
231:   VecGetArray(Fnet,&fnet);

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

245:     /* Generator differential equations */
246:     fgen[idx]   = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i];
247:     fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i];
248:     fgen[idx+2] = -w + w_s;
249:     fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i];

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

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

258:     Zdq_inv[0] = Rs[i]/det;
259:     Zdq_inv[1] = Xqp[i]/det;
260:     Zdq_inv[2] = -Xdp[i]/det;
261:     Zdq_inv[3] = Rs[i]/det;

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

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

269:     fnet[2*gbus[i]]   -= IGi;
270:     fnet[2*gbus[i]+1] -= IGr;

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

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

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

281:     idx = idx + 9;
282:   }

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

294:     /* Load currents */
295:     IDr = (PD*Vr + QD*Vi)/Vm2;
296:     IDi = (-QD*Vr + PD*Vi)/Vm2;

298:     fnet[2*lbus[i]]   += IDi;
299:     fnet[2*lbus[i]+1] += IDr;
300:   }
301:   VecRestoreArray(user->V0,&v0);

303:   VecRestoreArray(Xgen,&xgen);
304:   VecRestoreArray(Xnet,&xnet);
305:   VecRestoreArray(Fgen,&fgen);
306:   VecRestoreArray(Fnet,&fnet);

308:   DMCompositeGather(user->dmpgrid,INSERT_VALUES,F,Fgen,Fnet);
309:   DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);
310:   DMCompositeRestoreLocalVectors(user->dmpgrid,&Fgen,&Fnet);
311:   return 0;
312: }

314: /* \dot{x} - f(x,y)
315:      g(x,y) = 0
316:  */
317: PetscErrorCode IFunction(TS ts,PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
318: {
319:   SNES              snes;
320:   PetscScalar       *f;
321:   const PetscScalar *xdot;
322:   PetscInt          i;

324:   user->t = t;

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

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

355:   ResidualFunction(snes,X,F,user);
356:   VecGetArray(F,&f);
357:   for (i=0; i < ngen; i++) {
358:     f[9*i]   = 0;
359:     f[9*i+1] = 0;
360:     f[9*i+2] = 0;
361:     f[9*i+3] = 0;
362:     f[9*i+6] = 0;
363:     f[9*i+7] = 0;
364:     f[9*i+8] = 0;
365:   }
366:   VecRestoreArray(F,&f);
367:   return 0;
368: }

370: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
371: {
372:   PetscInt       *d_nnz;
373:   PetscInt       i,idx=0,start=0;
374:   PetscInt       ncols;

376:   PetscMalloc1(user->neqs_pgrid,&d_nnz);
377:   for (i=0; i<user->neqs_pgrid; i++) d_nnz[i] = 0;
378:   /* Generator subsystem */
379:   for (i=0; i < ngen; i++) {

381:     d_nnz[idx]   += 3;
382:     d_nnz[idx+1] += 2;
383:     d_nnz[idx+2] += 2;
384:     d_nnz[idx+3] += 5;
385:     d_nnz[idx+4] += 6;
386:     d_nnz[idx+5] += 6;

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

391:     d_nnz[idx+6] += 2;
392:     d_nnz[idx+7] += 2;
393:     d_nnz[idx+8] += 5;

395:     idx = idx + 9;
396:   }

398:   start = user->neqs_gen;

400:   for (i=0; i < nbus; i++) {
401:     MatGetRow(user->Ybus,2*i,&ncols,NULL,NULL);
402:     d_nnz[start+2*i]   += ncols;
403:     d_nnz[start+2*i+1] += ncols;
404:     MatRestoreRow(user->Ybus,2*i,&ncols,NULL,NULL);
405:   }

407:   MatSeqAIJSetPreallocation(J,0,d_nnz);

409:   PetscFree(d_nnz);
410:   return 0;
411: }

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

444:   MatZeroEntries(B);
445:   DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);
446:   DMCompositeScatter(user->dmpgrid,X,Xgen,Xnet);

448:   VecGetArray(Xgen,&xgen);
449:   VecGetArray(Xnet,&xnet);

451:   /* Generator subsystem */
452:   for (i=0; i < ngen; i++) {
453:     Eqp   = xgen[idx];
454:     Edp   = xgen[idx+1];
455:     delta = xgen[idx+2];
456:     Id    = xgen[idx+4];
457:     Iq    = xgen[idx+5];
458:     Efd   = xgen[idx+6];

460:     /*    fgen[idx]   = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
461:     row[0] = idx;
462:     col[0] = idx;           col[1] = idx+4;          col[2] = idx+6;
463:     val[0] = 1/ Td0p[i]; val[1] = (Xd[i] - Xdp[i])/ Td0p[i]; val[2] = -1/Td0p[i];

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

467:     /*    fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
468:     row[0] = idx + 1;
469:     col[0] = idx + 1;       col[1] = idx+5;
470:     val[0] = 1/Tq0p[i]; val[1] = -(Xq[i] - Xqp[i])/Tq0p[i];
471:     MatSetValues(J,1,row,2,col,val,INSERT_VALUES);

473:     /*    fgen[idx+2] = - w + w_s; */
474:     row[0] = idx + 2;
475:     col[0] = idx + 2; col[1] = idx + 3;
476:     val[0] = 0;       val[1] = -1;
477:     MatSetValues(J,1,row,2,col,val,INSERT_VALUES);

479:     /*    fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
480:     row[0] = idx + 3;
481:     col[0] = idx; col[1] = idx + 1; col[2] = idx + 3;       col[3] = idx + 4;                  col[4] = idx + 5;
482:     val[0] = Iq/M[i];  val[1] = Id/M[i];      val[2] = D[i]/M[i]; val[3] = (Edp + (Xqp[i]-Xdp[i])*Iq)/M[i]; val[4] = (Eqp + (Xqp[i] - Xdp[i])*Id)/M[i];
483:     MatSetValues(J,1,row,5,col,val,INSERT_VALUES);

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

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

491:     Zdq_inv[0] = Rs[i]/det;
492:     Zdq_inv[1] = Xqp[i]/det;
493:     Zdq_inv[2] = -Xdp[i]/det;
494:     Zdq_inv[3] = Rs[i]/det;

496:     dVd_dVr    = PetscSinScalar(delta); dVd_dVi = -PetscCosScalar(delta);
497:     dVq_dVr    = PetscCosScalar(delta); dVq_dVi = PetscSinScalar(delta);
498:     dVd_ddelta = Vr*PetscCosScalar(delta) + Vi*PetscSinScalar(delta);
499:     dVq_ddelta = -Vr*PetscSinScalar(delta) + Vi*PetscCosScalar(delta);

501:     /*    fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
502:     row[0] = idx+4;
503:     col[0] = idx;         col[1] = idx+1;        col[2] = idx + 2;
504:     val[0] = -Zdq_inv[1]; val[1] = -Zdq_inv[0];  val[2] = Zdq_inv[0]*dVd_ddelta + Zdq_inv[1]*dVq_ddelta;
505:     col[3] = idx + 4; col[4] = net_start+2*gbus[i];                     col[5] = net_start + 2*gbus[i]+1;
506:     val[3] = 1;       val[4] = Zdq_inv[0]*dVd_dVr + Zdq_inv[1]*dVq_dVr; val[5] = Zdq_inv[0]*dVd_dVi + Zdq_inv[1]*dVq_dVi;
507:     MatSetValues(J,1,row,6,col,val,INSERT_VALUES);

509:     /*  fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
510:     row[0] = idx+5;
511:     col[0] = idx;         col[1] = idx+1;        col[2] = idx + 2;
512:     val[0] = -Zdq_inv[3]; val[1] = -Zdq_inv[2];  val[2] = Zdq_inv[2]*dVd_ddelta + Zdq_inv[3]*dVq_ddelta;
513:     col[3] = idx + 5; col[4] = net_start+2*gbus[i];                     col[5] = net_start + 2*gbus[i]+1;
514:     val[3] = 1;       val[4] = Zdq_inv[2]*dVd_dVr + Zdq_inv[3]*dVq_dVr; val[5] = Zdq_inv[2]*dVd_dVi + Zdq_inv[3]*dVq_dVi;
515:     MatSetValues(J,1,row,6,col,val,INSERT_VALUES);

517:     dIGr_ddelta = Id*PetscCosScalar(delta) - Iq*PetscSinScalar(delta);
518:     dIGi_ddelta = Id*PetscSinScalar(delta) + Iq*PetscCosScalar(delta);
519:     dIGr_dId    = PetscSinScalar(delta);  dIGr_dIq = PetscCosScalar(delta);
520:     dIGi_dId    = -PetscCosScalar(delta); dIGi_dIq = PetscSinScalar(delta);

522:     /* fnet[2*gbus[i]]   -= IGi; */
523:     row[0] = net_start + 2*gbus[i];
524:     col[0] = idx+2;        col[1] = idx + 4;   col[2] = idx + 5;
525:     val[0] = -dIGi_ddelta; val[1] = -dIGi_dId; val[2] = -dIGi_dIq;
526:     MatSetValues(J,1,row,3,col,val,INSERT_VALUES);

528:     /* fnet[2*gbus[i]+1]   -= IGr; */
529:     row[0] = net_start + 2*gbus[i]+1;
530:     col[0] = idx+2;        col[1] = idx + 4;   col[2] = idx + 5;
531:     val[0] = -dIGr_ddelta; val[1] = -dIGr_dId; val[2] = -dIGr_dIq;
532:     MatSetValues(J,1,row,3,col,val,INSERT_VALUES);

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

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

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

541:     row[0] = idx + 6;
542:     col[0] = idx + 6;                     col[1] = idx + 8;
543:     val[0] = (KE[i] + dSE_dEfd)/TE[i];  val[1] = -1/TE[i];
544:     MatSetValues(J,1,row,2,col,val,INSERT_VALUES);

546:     /* Exciter differential equations */

548:     /*    fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
549:     row[0] = idx + 7;
550:     col[0] = idx + 6;       col[1] = idx + 7;
551:     val[0] = (-KF[i]/TF[i])/TF[i];  val[1] = 1/TF[i];
552:     MatSetValues(J,1,row,2,col,val,INSERT_VALUES);

554:     /*    fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
555:     /* Vm = (Vd^2 + Vq^2)^0.5; */
556:     dVm_dVd    = Vd/Vm; dVm_dVq = Vq/Vm;
557:     dVm_dVr    = dVm_dVd*dVd_dVr + dVm_dVq*dVq_dVr;
558:     dVm_dVi    = dVm_dVd*dVd_dVi + dVm_dVq*dVq_dVi;
559:     row[0]     = idx + 8;
560:     col[0]     = idx + 6;           col[1] = idx + 7; col[2] = idx + 8;
561:     val[0]     = (KA[i]*KF[i]/TF[i])/TA[i]; val[1] = -KA[i]/TA[i];  val[2] = 1/TA[i];
562:     col[3]     = net_start + 2*gbus[i]; col[4] = net_start + 2*gbus[i]+1;
563:     val[3]     = KA[i]*dVm_dVr/TA[i];         val[4] = KA[i]*dVm_dVi/TA[i];
564:     MatSetValues(J,1,row,5,col,val,INSERT_VALUES);
565:     idx        = idx + 9;
566:   }

568:   for (i=0; i<nbus; i++) {
569:     MatGetRow(user->Ybus,2*i,&ncols,&cols,&yvals);
570:     row[0] = net_start + 2*i;
571:     for (k=0; k<ncols; k++) {
572:       col[k] = net_start + cols[k];
573:       val[k] = yvals[k];
574:     }
575:     MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);
576:     MatRestoreRow(user->Ybus,2*i,&ncols,&cols,&yvals);

578:     MatGetRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);
579:     row[0] = net_start + 2*i+1;
580:     for (k=0; k<ncols; k++) {
581:       col[k] = net_start + cols[k];
582:       val[k] = yvals[k];
583:     }
584:     MatSetValues(J,1,row,ncols,col,val,INSERT_VALUES);
585:     MatRestoreRow(user->Ybus,2*i+1,&ncols,&cols,&yvals);
586:   }

588:   MatAssemblyBegin(J,MAT_FLUSH_ASSEMBLY);
589:   MatAssemblyEnd(J,MAT_FLUSH_ASSEMBLY);

591:   VecGetArray(user->V0,&v0);
592:   for (i=0; i < nload; i++) {
593:     Vr      = xnet[2*lbus[i]]; /* Real part of load bus voltage */
594:     Vi      = xnet[2*lbus[i]+1]; /* Imaginary part of the load bus voltage */
595:     Vm      = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm; Vm4 = Vm2*Vm2;
596:     Vm0     = PetscSqrtScalar(v0[2*lbus[i]]*v0[2*lbus[i]] + v0[2*lbus[i]+1]*v0[2*lbus[i]+1]);
597:     PD      = QD = 0.0;
598:     dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
599:     for (k=0; k < ld_nsegsp[i]; k++) {
600:       PD      += ld_alphap[k]*PD0[i]*PetscPowScalar((Vm/Vm0),ld_betap[k]);
601:       dPD_dVr += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vr*PetscPowScalar(Vm,(ld_betap[k]-2));
602:       dPD_dVi += ld_alphap[k]*ld_betap[k]*PD0[i]*PetscPowScalar((1/Vm0),ld_betap[k])*Vi*PetscPowScalar(Vm,(ld_betap[k]-2));
603:     }
604:     for (k=0; k < ld_nsegsq[i]; k++) {
605:       QD      += ld_alphaq[k]*QD0[i]*PetscPowScalar((Vm/Vm0),ld_betaq[k]);
606:       dQD_dVr += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vr*PetscPowScalar(Vm,(ld_betaq[k]-2));
607:       dQD_dVi += ld_alphaq[k]*ld_betaq[k]*QD0[i]*PetscPowScalar((1/Vm0),ld_betaq[k])*Vi*PetscPowScalar(Vm,(ld_betaq[k]-2));
608:     }

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

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

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

619:     /*    fnet[2*lbus[i]]   += IDi; */
620:     row[0] = net_start + 2*lbus[i];
621:     col[0] = net_start + 2*lbus[i];  col[1] = net_start + 2*lbus[i]+1;
622:     val[0] = dIDi_dVr;               val[1] = dIDi_dVi;
623:     MatSetValues(J,1,row,2,col,val,ADD_VALUES);
624:     /*    fnet[2*lbus[i]+1] += IDr; */
625:     row[0] = net_start + 2*lbus[i]+1;
626:     col[0] = net_start + 2*lbus[i];  col[1] = net_start + 2*lbus[i]+1;
627:     val[0] = dIDr_dVr;               val[1] = dIDr_dVi;
628:     MatSetValues(J,1,row,2,col,val,ADD_VALUES);
629:   }
630:   VecRestoreArray(user->V0,&v0);

632:   VecRestoreArray(Xgen,&xgen);
633:   VecRestoreArray(Xnet,&xnet);

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

637:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
638:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
639:   return 0;
640: }

642: /*
643:    J = [I, 0
644:         dg_dx, dg_dy]
645: */
646: PetscErrorCode AlgJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
647: {
648:   Userctx        *user=(Userctx*)ctx;

650:   ResidualJacobian(snes,X,A,B,ctx);
651:   MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
652:   MatZeroRowsIS(A,user->is_diff,1.0,NULL,NULL);
653:   return 0;
654: }

656: /*
657:    J = [a*I-df_dx, -df_dy
658:         dg_dx, dg_dy]
659: */

661: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx *user)
662: {
663:   SNES           snes;
664:   PetscScalar    atmp = (PetscScalar) a;
665:   PetscInt       i,row;

667:   user->t = t;

669:   TSGetSNES(ts,&snes);
670:   ResidualJacobian(snes,X,A,B,user);
671:   for (i=0;i < ngen;i++) {
672:     row = 9*i;
673:     MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);
674:     row  = 9*i+1;
675:     MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);
676:     row  = 9*i+2;
677:     MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);
678:     row  = 9*i+3;
679:     MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);
680:     row  = 9*i+6;
681:     MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);
682:     row  = 9*i+7;
683:     MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);
684:     row  = 9*i+8;
685:     MatSetValues(A,1,&row,1,&row,&atmp,ADD_VALUES);
686:   }
687:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
688:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
689:   return 0;
690: }

692: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,Userctx *user)
693: {
694:   PetscScalar       *r;
695:   const PetscScalar *u;
696:   PetscInt          idx;
697:   Vec               Xgen,Xnet;
698:   PetscScalar       *xgen;
699:   PetscInt          i;

701:   DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);
702:   DMCompositeScatter(user->dmpgrid,U,Xgen,Xnet);

704:   VecGetArray(Xgen,&xgen);

706:   VecGetArrayRead(U,&u);
707:   VecGetArray(R,&r);
708:   r[0] = 0.;

710:   idx = 0;
711:   for (i=0;i<ngen;i++) {
712:     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);
713:     idx  += 9;
714:   }
715:   VecRestoreArray(R,&r);
716:   VecRestoreArrayRead(U,&u);
717:   DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);
718:   return 0;
719: }

721: static PetscErrorCode MonitorUpdateQ(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx0)
722: {
723:   Vec            C,*Y;
724:   PetscInt       Nr;
725:   PetscReal      h,theta;
726:   Userctx        *ctx=(Userctx*)ctx0;

728:   theta = 0.5;
729:   TSGetStages(ts,&Nr,&Y);
730:   TSGetTimeStep(ts,&h);
731:   VecDuplicate(ctx->vec_q,&C);
732:   /* compute integrals */
733:   if (stepnum>0) {
734:     CostIntegrand(ts,time,X,C,ctx);
735:     VecAXPY(ctx->vec_q,h*theta,C);
736:     CostIntegrand(ts,time+h*theta,Y[0],C,ctx);
737:     VecAXPY(ctx->vec_q,h*(1-theta),C);
738:   }
739:   VecDestroy(&C);
740:   return 0;
741: }

743: int main(int argc,char **argv)
744: {
745:   Userctx            user;
746:   Vec                p;
747:   PetscScalar        *x_ptr;
748:   PetscErrorCode     ierr;
749:   PetscMPIInt        size;
750:   PetscInt           i;
751:   KSP                ksp;
752:   PC                 pc;
753:   PetscInt           *idx2;
754:   Tao                tao;
755:   Vec                lowerb,upperb;

758:   PetscInitialize(&argc,&argv,"petscoptions",help);
759:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

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

764:   user.neqs_gen   = 9*ngen; /* # eqs. for generator subsystem */
765:   user.neqs_net   = 2*nbus; /* # eqs. for network subsystem   */
766:   user.neqs_pgrid = user.neqs_gen + user.neqs_net;

768:   /* Create indices for differential and algebraic equations */
769:   PetscMalloc1(7*ngen,&idx2);
770:   for (i=0; i<ngen; i++) {
771:     idx2[7*i]   = 9*i;   idx2[7*i+1] = 9*i+1; idx2[7*i+2] = 9*i+2; idx2[7*i+3] = 9*i+3;
772:     idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8;
773:   }
774:   ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff);
775:   ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg);
776:   PetscFree(idx2);

778:   /* Set run time options */
779:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");
780:   {
781:     user.tfaulton  = 1.0;
782:     user.tfaultoff = 1.2;
783:     user.Rfault    = 0.0001;
784:     user.faultbus  = 8;
785:     PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL);
786:     PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL);
787:     PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL);
788:     user.t0        = 0.0;
789:     user.tmax      = 1.5;
790:     PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL);
791:     PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL);
792:     user.freq_u    = 61.0;
793:     user.freq_l    = 59.0;
794:     user.pow       = 2;
795:     PetscOptionsReal("-frequ","","",user.freq_u,&user.freq_u,NULL);
796:     PetscOptionsReal("-freql","","",user.freq_l,&user.freq_l,NULL);
797:     PetscOptionsInt("-pow","","",user.pow,&user.pow,NULL);

799:   }
800:   PetscOptionsEnd();

802:   /* Create DMs for generator and network subsystems */
803:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);
804:   DMSetOptionsPrefix(user.dmgen,"dmgen_");
805:   DMSetFromOptions(user.dmgen);
806:   DMSetUp(user.dmgen);
807:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);
808:   DMSetOptionsPrefix(user.dmnet,"dmnet_");
809:   DMSetFromOptions(user.dmnet);
810:   DMSetUp(user.dmnet);
811:   /* Create a composite DM packer and add the two DMs */
812:   DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid);
813:   DMSetOptionsPrefix(user.dmpgrid,"pgrid_");
814:   DMCompositeAddDM(user.dmpgrid,user.dmgen);
815:   DMCompositeAddDM(user.dmpgrid,user.dmnet);

817:   /* Create TAO solver and set desired solution method */
818:   TaoCreate(PETSC_COMM_WORLD,&tao);
819:   TaoSetType(tao,TAOBLMVM);
820:   /*
821:      Optimization starts
822:   */
823:   /* Set initial solution guess */
824:   VecCreateSeq(PETSC_COMM_WORLD,3,&p);
825:   VecGetArray(p,&x_ptr);
826:   x_ptr[0] = PG[0]; x_ptr[1] = PG[1]; x_ptr[2] = PG[2];
827:   VecRestoreArray(p,&x_ptr);

829:   TaoSetSolution(tao,p);
830:   /* Set routine for function and gradient evaluation */
831:   TaoSetObjective(tao,FormFunction,(void *)&user);
832:   TaoSetGradient(tao,NULL,TaoDefaultComputeGradient,(void *)&user);

834:   /* Set bounds for the optimization */
835:   VecDuplicate(p,&lowerb);
836:   VecDuplicate(p,&upperb);
837:   VecGetArray(lowerb,&x_ptr);
838:   x_ptr[0] = 0.5; x_ptr[1] = 0.5; x_ptr[2] = 0.5;
839:   VecRestoreArray(lowerb,&x_ptr);
840:   VecGetArray(upperb,&x_ptr);
841:   x_ptr[0] = 2.0; x_ptr[1] = 2.0; x_ptr[2] = 2.0;
842:   VecRestoreArray(upperb,&x_ptr);
843:   TaoSetVariableBounds(tao,lowerb,upperb);

845:   /* Check for any TAO command line options */
846:   TaoSetFromOptions(tao);
847:   TaoGetKSP(tao,&ksp);
848:   if (ksp) {
849:     KSPGetPC(ksp,&pc);
850:     PCSetType(pc,PCNONE);
851:   }

853:   /* SOLVE THE APPLICATION */
854:   TaoSolve(tao);

856:   VecView(p,PETSC_VIEWER_STDOUT_WORLD);
857:   /* Free TAO data structures */
858:   TaoDestroy(&tao);
859:   VecDestroy(&user.vec_q);
860:   VecDestroy(&lowerb);
861:   VecDestroy(&upperb);
862:   VecDestroy(&p);
863:   DMDestroy(&user.dmgen);
864:   DMDestroy(&user.dmnet);
865:   DMDestroy(&user.dmpgrid);
866:   ISDestroy(&user.is_diff);
867:   ISDestroy(&user.is_alg);
868:   PetscFinalize();
869:   return 0;
870: }

872: /* ------------------------------------------------------------------ */
873: /*
874:    FormFunction - Evaluates the function and corresponding gradient.

876:    Input Parameters:
877:    tao - the Tao context
878:    X   - the input vector
879:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

881:    Output Parameters:
882:    f   - the newly evaluated function
883: */
884: PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0)
885: {
886:   TS             ts;
887:   SNES           snes_alg;
888:   Userctx        *ctx = (Userctx*)ctx0;
889:   Vec            X;
890:   Mat            J;
891:   /* sensitivity context */
892:   PetscScalar    *x_ptr;
893:   PetscViewer    Xview,Ybusview;
894:   Vec            F_alg;
895:   Vec            Xdot;
896:   PetscInt       row_loc,col_loc;
897:   PetscScalar    val;

899:   VecGetArrayRead(P,(const PetscScalar**)&x_ptr);
900:   PG[0] = x_ptr[0];
901:   PG[1] = x_ptr[1];
902:   PG[2] = x_ptr[2];
903:   VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);

905:   ctx->stepnum = 0;

907:   VecZeroEntries(ctx->vec_q);

909:   /* Read initial voltage vector and Ybus */
910:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"X.bin",FILE_MODE_READ,&Xview);
911:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Ybus.bin",FILE_MODE_READ,&Ybusview);

913:   VecCreate(PETSC_COMM_WORLD,&ctx->V0);
914:   VecSetSizes(ctx->V0,PETSC_DECIDE,ctx->neqs_net);
915:   VecLoad(ctx->V0,Xview);

917:   MatCreate(PETSC_COMM_WORLD,&ctx->Ybus);
918:   MatSetSizes(ctx->Ybus,PETSC_DECIDE,PETSC_DECIDE,ctx->neqs_net,ctx->neqs_net);
919:   MatSetType(ctx->Ybus,MATBAIJ);
920:   /*  MatSetBlockSize(ctx->Ybus,2); */
921:   MatLoad(ctx->Ybus,Ybusview);

923:   PetscViewerDestroy(&Xview);
924:   PetscViewerDestroy(&Ybusview);

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

928:   MatCreate(PETSC_COMM_WORLD,&J);
929:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,ctx->neqs_pgrid,ctx->neqs_pgrid);
930:   MatSetFromOptions(J);
931:   PreallocateJacobian(J,ctx);

933:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
934:      Create timestepping solver context
935:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
936:   TSCreate(PETSC_COMM_WORLD,&ts);
937:   TSSetProblemType(ts,TS_NONLINEAR);
938:   TSSetType(ts,TSCN);
939:   TSSetIFunction(ts,NULL,(TSIFunction) IFunction,ctx);
940:   TSSetIJacobian(ts,J,J,(TSIJacobian)IJacobian,ctx);
941:   TSSetApplicationContext(ts,ctx);

943:   TSMonitorSet(ts,MonitorUpdateQ,ctx,NULL);
944:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
945:      Set initial conditions
946:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
947:   SetInitialGuess(X,ctx);

949:   VecDuplicate(X,&F_alg);
950:   SNESCreate(PETSC_COMM_WORLD,&snes_alg);
951:   SNESSetFunction(snes_alg,F_alg,AlgFunction,ctx);
952:   MatZeroEntries(J);
953:   SNESSetJacobian(snes_alg,J,J,AlgJacobian,ctx);
954:   SNESSetOptionsPrefix(snes_alg,"alg_");
955:   SNESSetFromOptions(snes_alg);
956:   ctx->alg_flg = PETSC_TRUE;
957:   /* Solve the algebraic equations */
958:   SNESSolve(snes_alg,NULL,X);

960:   /* Just to set up the Jacobian structure */
961:   VecDuplicate(X,&Xdot);
962:   IJacobian(ts,0.0,X,Xdot,0.0,J,J,ctx);
963:   VecDestroy(&Xdot);

965:   ctx->stepnum++;

967:   TSSetTimeStep(ts,0.01);
968:   TSSetMaxTime(ts,ctx->tmax);
969:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
970:   TSSetFromOptions(ts);

972:   /* Prefault period */
973:   ctx->alg_flg = PETSC_FALSE;
974:   TSSetTime(ts,0.0);
975:   TSSetMaxTime(ts,ctx->tfaulton);
976:   TSSolve(ts,X);

978:   /* Create the nonlinear solver for solving the algebraic system */
979:   /* Note that although the algebraic system needs to be solved only for
980:      Idq and V, we reuse the entire system including xgen. The xgen
981:      variables are held constant by setting their residuals to 0 and
982:      putting a 1 on the Jacobian diagonal for xgen rows
983:   */
984:   MatZeroEntries(J);

986:   /* Apply disturbance - resistive fault at ctx->faultbus */
987:   /* This is done by adding shunt conductance to the diagonal location
988:      in the Ybus matrix */
989:   row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1; /* Location for G */
990:   val     = 1/ctx->Rfault;
991:   MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);
992:   row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus; /* Location for G */
993:   val     = 1/ctx->Rfault;
994:   MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);

996:   MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY);
997:   MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY);

999:   ctx->alg_flg = PETSC_TRUE;
1000:   /* Solve the algebraic equations */
1001:   SNESSolve(snes_alg,NULL,X);

1003:   ctx->stepnum++;

1005:   /* Disturbance period */
1006:   ctx->alg_flg = PETSC_FALSE;
1007:   TSSetTime(ts,ctx->tfaulton);
1008:   TSSetMaxTime(ts,ctx->tfaultoff);
1009:   TSSolve(ts,X);

1011:   /* Remove the fault */
1012:   row_loc = 2*ctx->faultbus; col_loc = 2*ctx->faultbus+1;
1013:   val     = -1/ctx->Rfault;
1014:   MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);
1015:   row_loc = 2*ctx->faultbus+1; col_loc = 2*ctx->faultbus;
1016:   val     = -1/ctx->Rfault;
1017:   MatSetValues(ctx->Ybus,1,&row_loc,1,&col_loc,&val,ADD_VALUES);

1019:   MatAssemblyBegin(ctx->Ybus,MAT_FINAL_ASSEMBLY);
1020:   MatAssemblyEnd(ctx->Ybus,MAT_FINAL_ASSEMBLY);

1022:   MatZeroEntries(J);

1024:   ctx->alg_flg = PETSC_TRUE;

1026:   /* Solve the algebraic equations */
1027:   SNESSolve(snes_alg,NULL,X);

1029:   ctx->stepnum++;

1031:   /* Post-disturbance period */
1032:   ctx->alg_flg = PETSC_TRUE;
1033:   TSSetTime(ts,ctx->tfaultoff);
1034:   TSSetMaxTime(ts,ctx->tmax);
1035:   TSSolve(ts,X);

1037:   VecGetArray(ctx->vec_q,&x_ptr);
1038:   *f   = x_ptr[0];
1039:   VecRestoreArray(ctx->vec_q,&x_ptr);

1041:   MatDestroy(&ctx->Ybus);
1042:   VecDestroy(&ctx->V0);
1043:   SNESDestroy(&snes_alg);
1044:   VecDestroy(&F_alg);
1045:   MatDestroy(&J);
1046:   VecDestroy(&X);
1047:   TSDestroy(&ts);

1049:   return 0;
1050: }

1052: /*TEST

1054:   build:
1055:       requires: double !complex !defined(USE_64BIT_INDICES)

1057:    test:
1058:       args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1059:       localrunfiles: petscoptions X.bin Ybus.bin

1061: TEST*/