2: static char help[] = "Power grid stability analysis of a large power system.\n\
  3: The base case is the WECC 9 bus system based on the example given in the book Power\n\
  4: Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\
  5: The base power grid in this example consists of 9 buses (nodes), 3 generators,\n\
  6: 3 loads, and 9 transmission lines. The network equations are written\n\
  7: in current balance form using rectangular coordiantes. \n\
  8: DMNetwork is used to manage the variables and equations of the entire system.\n\
  9: Input parameters include:\n\
 10:   -nc : number of copies of the base case\n\n";

 12: /* T
 13:    Concepts: DMNetwork
 14:    Concepts: PETSc TS solver

 16:    This example was contributed by Bikiran Guha and Jianqiao Huang, Illinois Institute of Technology, 2017.
 17: */

 19: #include <petscts.h>
 20: #include <petscdmnetwork.h>

 22: #define FREQ 60
 23: #define W_S (2*PETSC_PI*FREQ)
 24: #define NGEN    3   /* No. of generators in the 9 bus system */
 25: #define NLOAD   3   /* No. of loads in the 9 bus system */
 26: #define NBUS    9   /* No. of buses in the 9 bus system */
 27: #define NBRANCH 9   /* No. of branches in the 9 bus system */

 29: typedef struct {
 30:   PetscInt    id;    /* Bus Number or extended bus name*/
 31:   PetscScalar mbase; /* MVA base of the machine */
 32:   PetscScalar PG;    /* Generator active power output */
 33:   PetscScalar QG;    /* Generator reactive power output */

 35:   /* Generator constants */
 36:   PetscScalar H;    /* Inertia constant */
 37:   PetscScalar Rs;   /* Stator Resistance */
 38:   PetscScalar Xd;   /* d-axis reactance */
 39:   PetscScalar Xdp;  /* d-axis transient reactance */
 40:   PetscScalar Xq;   /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
 41:   PetscScalar Xqp;  /* q-axis transient reactance */
 42:   PetscScalar Td0p; /* d-axis open circuit time constant */
 43:   PetscScalar Tq0p; /* q-axis open circuit time constant */
 44:   PetscScalar M;    /* M = 2*H/W_S */
 45:   PetscScalar D;    /* D = 0.1*M */
 46:   PetscScalar TM;   /* Mechanical Torque */

 48:   /* Exciter system constants */
 49:   PetscScalar KA ;   /* Voltage regulator gain constant */
 50:   PetscScalar TA;    /* Voltage regulator time constant */
 51:   PetscScalar KE;    /* Exciter gain constant */
 52:   PetscScalar TE;    /* Exciter time constant */
 53:   PetscScalar KF;    /* Feedback stabilizer gain constant */
 54:   PetscScalar TF;    /* Feedback stabilizer time constant */
 55:   PetscScalar k1,k2; /* calculating the saturation function SE = k1*exp(k2*Efd) */
 56:   PetscScalar Vref;  /* Voltage regulator voltage setpoint */
 57: } Gen;

 59: typedef struct {
 60:   PetscInt     id;      /* node id */
 61:   PetscInt     nofgen;  /* Number of generators at the bus*/
 62:   PetscInt     nofload; /*  Number of load at the bus*/
 63:   PetscScalar  yff[2]; /* yff[0]= imaginary part of admittance, yff[1]=real part of admittance*/
 64:   PetscScalar  vr;     /* Real component of bus voltage */
 65:   PetscScalar  vi;     /* Imaginary component of bus voltage */
 66: } Bus;

 68:   /* Load constants
 69:   We use a composite load model that describes the load and reactive powers at each time instant as follows
 70:   P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
 71:   Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
 72:   where
 73:     id                  - index of the load
 74:     ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
 75:     ld_alphap,ld_alphap - Percentage contribution (weights) or loads
 76:     P_D0                - Real power load
 77:     Q_D0                - Reactive power load
 78:     Vm(t)              - Voltage magnitude at time t
 79:     Vm0                - Voltage magnitude at t = 0
 80:     ld_betap, ld_betaq  - exponents describing the load model for real and reactive part

 82:     Note: All loads have the same characteristic currently.
 83:   */
 84: typedef struct {
 85:   PetscInt    id;           /* bus id */
 86:   PetscInt    ld_nsegsp,ld_nsegsq;
 87:   PetscScalar PD0, QD0;
 88:   PetscScalar ld_alphap[3]; /* ld_alphap=[1,0,0], an array, not a value, so use *ld_alphap; */
 89:   PetscScalar ld_betap[3],ld_alphaq[3],ld_betaq[3];
 90: } Load;

 92: typedef struct {
 93:   PetscInt    id;     /* node id */
 94:   PetscScalar yft[2]; /* yft[0]= imaginary part of admittance, yft[1]=real part of admittance*/
 95: } Branch;

 97: typedef struct {
 98:   PetscReal   tfaulton,tfaultoff; /* Fault on and off times */
 99:   PetscReal   t;
100:   PetscReal   t0,tmax;            /* initial time and final time */
101:   PetscInt    faultbus;           /* Fault bus */
102:   PetscScalar Rfault;             /* Fault resistance (pu) */
103:   PetscScalar *ybusfault;
104:   PetscBool   alg_flg;
105: } Userctx;

107: /* Used to read data into the DMNetwork components */
108: PetscErrorCode read_data(PetscInt nc, Gen **pgen,Load **pload,Bus **pbus, Branch **pbranch, PetscInt **pedgelist)
109: {
110:   PetscErrorCode    ierr;
111:   PetscInt          i,j,row[1],col[2];
112:   PetscInt          *edgelist;
113:   PetscInt          nofgen[9] = {1,1,1,0,0,0,0,0,0}; /* Buses at which generators are incident */
114:   PetscInt          nofload[9] = {0,0,0,0,1,1,0,1,0}; /* Buses at which loads are incident */
115:   PetscScalar       *varr,M[3],D[3];
116:   Bus               *bus;
117:   Branch            *branch;
118:   Gen               *gen;
119:   Load              *load;
120:   Mat               Ybus;
121:   Vec               V0;

123:   /*10 parameters*/
124:   /* Generator real and reactive powers (found via loadflow) */
125:   static const PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};
126:   static const PetscScalar QG[3] = {0.270702180178785,0.066120127797275,-0.108402221791588};

128:   /* Generator constants */
129:   static const PetscScalar H[3]    = {23.64,6.4,3.01};   /* Inertia constant */
130:   static const PetscScalar Rs[3]   = {0.0,0.0,0.0}; /* Stator Resistance */
131:   static const PetscScalar Xd[3]   = {0.146,0.8958,1.3125};  /* d-axis reactance */
132:   static const PetscScalar Xdp[3]  = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */
133:   static 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 */
134:   static const PetscScalar Xqp[3]  = {0.0969,0.1969,0.25}; /* q-axis transient reactance */
135:   static const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */
136:   static const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */

138:   /* Exciter system constants (8 parameters)*/
139:   static const PetscScalar KA[3] = {20.0,20.0,20.0};  /* Voltage regulartor gain constant */
140:   static const PetscScalar TA[3] = {0.2,0.2,0.2};     /* Voltage regulator time constant */
141:   static const PetscScalar KE[3] = {1.0,1.0,1.0};     /* Exciter gain constant */
142:   static const PetscScalar TE[3] = {0.314,0.314,0.314}; /* Exciter time constant */
143:   static const PetscScalar KF[3] = {0.063,0.063,0.063};  /* Feedback stabilizer gain constant */
144:   static const PetscScalar TF[3] = {0.35,0.35,0.35};    /* Feedback stabilizer time constant */
145:   static const PetscScalar k1[3] = {0.0039,0.0039,0.0039};
146:   static const PetscScalar k2[3] = {1.555,1.555,1.555};  /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */

148:   /* Load constants */
149:    static const PetscScalar       PD0[3]       = {1.25,0.9,1.0};
150:    static const PetscScalar       QD0[3]       = {0.5,0.3,0.35};
151:    static const PetscScalar       ld_alphaq[3] = {1,0,0};
152:    static const PetscScalar       ld_betaq[3]  = {2,1,0};
153:    static const PetscScalar       ld_betap[3]  = {2,1,0};
154:    static const PetscScalar       ld_alphap[3] = {1,0,0};
155:    PetscInt                       ld_nsegsp[3] = {3,3,3};
156:    PetscInt                       ld_nsegsq[3] = {3,3,3};
157:    PetscViewer                    Xview,Ybusview;
158:    PetscInt                       neqs_net,m,n;

161:    /* Read V0 and Ybus from files */
162:    PetscViewerBinaryOpen(PETSC_COMM_SELF,"X.bin",FILE_MODE_READ,&Xview);
163:    PetscViewerBinaryOpen(PETSC_COMM_SELF,"Ybus.bin",FILE_MODE_READ,&Ybusview);
164:    VecCreate(PETSC_COMM_SELF,&V0);
165:    VecLoad(V0,Xview);

167:    MatCreate(PETSC_COMM_SELF,&Ybus);
168:    MatSetType(Ybus,MATBAIJ);
169:    MatLoad(Ybus,Ybusview);

171:    /* Destroy unnecessary stuff */
172:    PetscViewerDestroy(&Xview);
173:    PetscViewerDestroy(&Ybusview);

175:    MatGetLocalSize(Ybus,&m,&n);
176:    neqs_net = 2*NBUS; /* # eqs. for network subsystem   */
177:    if (m != neqs_net || n != neqs_net) SETERRQ(PETSC_COMM_SELF,0,"matrix Ybus is in wrong sizes");

179:    M[0] = 2*H[0]/W_S;
180:    M[1] = 2*H[1]/W_S;
181:    M[2] = 2*H[2]/W_S;
182:    D[0] = 0.1*M[0];
183:    D[1] = 0.1*M[1];
184:    D[2] = 0.1*M[2];

186:    /* Alocate memory for total number of buses, generators, loads and branches */
187:    PetscCalloc4(NBUS*nc,&bus,NGEN*nc,&gen,NLOAD*nc,&load,NBRANCH*nc+(nc-1),&branch);

189:    VecGetArray(V0,&varr);

191:    /* read bus data */
192:    for (i = 0; i < nc; i++) {
193:      for (j = 0; j < NBUS; j++) {
194:        bus[i*9+j].id      = i*9+j;
195:        bus[i*9+j].nofgen  = nofgen[j];
196:        bus[i*9+j].nofload = nofload[j];
197:        bus[i*9+j].vr      = varr[2*j];
198:        bus[i*9+j].vi      = varr[2*j+1];
199:        row[0]             = 2*j;
200:        col[0]             = 2*j;
201:        col[1]             = 2*j+1;
202:        /* real and imaginary part of admittance from Ybus into yff */
203:        MatGetValues(Ybus,1,row,2,col,bus[i*9+j].yff);
204:      }
205:    }

207:    /* read generator data */
208:    for (i = 0; i<nc; i++){
209:      for (j = 0; j < NGEN; j++) {
210:        gen[i*3+j].id   = i*3+j;
211:        gen[i*3+j].PG   = PG[j];
212:        gen[i*3+j].QG   = QG[j];
213:        gen[i*3+j].H    = H[j];
214:        gen[i*3+j].Rs   = Rs[j];
215:        gen[i*3+j].Xd   = Xd[j];
216:        gen[i*3+j].Xdp  = Xdp[j];
217:        gen[i*3+j].Xq   = Xq[j];
218:        gen[i*3+j].Xqp  = Xqp[j];
219:        gen[i*3+j].Td0p = Td0p[j];
220:        gen[i*3+j].Tq0p = Tq0p[j];
221:        gen[i*3+j].M    = M[j];
222:        gen[i*3+j].D    = D[j];

224:        /* exciter system */
225:        gen[i*3+j].KA = KA[j];
226:        gen[i*3+j].TA = TA[j];
227:        gen[i*3+j].KE = KE[j];
228:        gen[i*3+j].TE = TE[j];
229:        gen[i*3+j].KF = KF[j];
230:        gen[i*3+j].TF = TF[j];
231:        gen[i*3+j].k1 = k1[j];
232:        gen[i*3+j].k2 = k2[j];
233:      }
234:    }

236:    /* read load data */
237:    for (i = 0; i<nc; i++){
238:      for (j = 0; j < NLOAD; j++) {
239:        load[i*3+j].id        = i*3+j;
240:        load[i*3+j].PD0       = PD0[j];
241:        load[i*3+j].QD0       = QD0[j];
242:        load[i*3+j].ld_nsegsp = ld_nsegsp[j];

244:        load[i*3+j].ld_alphap[0] = ld_alphap[0];
245:        load[i*3+j].ld_alphap[1] = ld_alphap[1];
246:        load[i*3+j].ld_alphap[2] = ld_alphap[2];

248:        load[i*3+j].ld_alphaq[0] = ld_alphaq[0];
249:        load[i*3+j].ld_alphaq[1] = ld_alphaq[1];
250:        load[i*3+j].ld_alphaq[2] = ld_alphaq[2];

252:        load[i*3+j].ld_betap[0] = ld_betap[0];
253:        load[i*3+j].ld_betap[1] = ld_betap[1];
254:        load[i*3+j].ld_betap[2] = ld_betap[2];
255:        load[i*3+j].ld_nsegsq   = ld_nsegsq[j];

257:        load[i*3+j].ld_betaq[0] = ld_betaq[0];
258:        load[i*3+j].ld_betaq[1] = ld_betaq[1];
259:        load[i*3+j].ld_betaq[2] = ld_betaq[2];
260:      }
261:    }
262:    PetscCalloc1(2*NBRANCH*nc+2*(nc-1),&edgelist);

264:    /* read edgelist */
265:    for (i = 0; i<nc; i++){
266:      for (j = 0; j < NBRANCH; j++) {
267:        switch (j) {
268:        case 0:
269:          edgelist[i*18+2*j]    = 0+9*i;
270:          edgelist[i*18+2*j+1]  = 3+9*i;
271:          break;
272:        case 1:
273:          edgelist[i*18+2*j]    = 1+9*i;
274:          edgelist[i*18+2*j+1]  = 6+9*i;
275:          break;
276:        case 2:
277:          edgelist[i*18+2*j]    = 2+9*i;
278:          edgelist[i*18+2*j+1]  = 8+9*i;
279:          break;
280:        case 3:
281:          edgelist[i*18+2*j]    = 3+9*i;
282:          edgelist[i*18+2*j+1]  = 4+9*i;
283:          break;
284:        case 4:
285:          edgelist[i*18+2*j]    = 3+9*i;
286:          edgelist[i*18+2*j+1]  = 5+9*i;
287:          break;
288:        case 5:
289:          edgelist[i*18+2*j]    = 4+9*i;
290:          edgelist[i*18+2*j+1]  = 6+9*i;
291:          break;
292:        case 6:
293:          edgelist[i*18+2*j]    = 5+9*i;
294:          edgelist[i*18+2*j+1]  = 8+9*i;
295:          break;
296:        case 7:
297:          edgelist[i*18+2*j]     = 6+9*i;
298:          edgelist[i*18+2*j+1]   = 7+9*i;
299:          break;
300:        case 8:
301:          edgelist[i*18+2*j]     = 7+9*i;
302:          edgelist[i*18+2*j+1]   = 8+9*i;
303:          break;
304:        default:
305:          break;
306:        }
307:      }
308:    }

310:    /* for connecting last bus of previous network(9*i-1) to first bus of next network(9*i), the branch admittance=-0.0301407+j17.3611 */
311:     for (i = 1; i<nc; i++){
312:         edgelist[18*nc+2*(i-1)]   = 8+(i-1)*9;
313:         edgelist[18*nc+2*(i-1)+1] = 9*i;

315:         /* adding admittances to the off-diagonal elements */
316:         branch[9*nc+(i-1)].id     = 9*nc+(i-1);
317:         branch[9*nc+(i-1)].yft[0] = 17.3611;
318:         branch[9*nc+(i-1)].yft[1] = -0.0301407;

320:         /* subtracting admittances from the diagonal elements */
321:         bus[9*i-1].yff[0] -= 17.3611;
322:         bus[9*i-1].yff[1] -= -0.0301407;

324:         bus[9*i].yff[0]   -= 17.3611;
325:         bus[9*i].yff[1]   -= -0.0301407;
326:     }

328:     /* read branch data */
329:     for (i = 0; i<nc; i++){
330:       for (j = 0; j < NBRANCH; j++) {
331:         branch[i*9+j].id  = i*9+j;

333:         row[0] = edgelist[2*j]*2;
334:         col[0] = edgelist[2*j+1]*2;
335:         col[1] = edgelist[2*j+1]*2+1;
336:         MatGetValues(Ybus,1,row,2,col,branch[i*9+j].yft);/*imaginary part of admittance*/
337:       }
338:     }

340:    *pgen      = gen;
341:    *pload     = load;
342:    *pbus      = bus;
343:    *pbranch   = branch;
344:    *pedgelist = edgelist;

346:    VecRestoreArray(V0,&varr);

348:    /* Destroy unnecessary stuff */
349:    MatDestroy(&Ybus);
350:    VecDestroy(&V0);
351:    return(0);
352: }

354: PetscErrorCode SetInitialGuess(DM networkdm, Vec X)
355: {
357:   Bus            *bus;
358:   Gen            *gen;
359:   PetscInt       v,vStart,vEnd,offset;
360:   PetscInt       key,numComps,j;
361:   PetscInt       idx=0;
362:   PetscBool      ghostvtex;
363:   Vec            localX;
364:   PetscScalar    *xarr;
365:   PetscScalar    Vr=0,Vi=0,Vm,Vm2;  /* Terminal voltage variables */
366:   PetscScalar    IGr, IGi;          /* Generator real and imaginary current */
367:   PetscScalar    Eqp,Edp,delta;     /* Generator variables */
368:   PetscScalar    Efd,RF,VR;         /* Exciter variables */
369:   PetscScalar    Vd,Vq;             /* Generator dq axis voltages */
370:   PetscScalar    Id,Iq;             /* Generator dq axis currents */
371:   PetscScalar    theta;             /* Generator phase angle */
372:   PetscScalar    SE;
373:   void*          component;

376:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
377:   DMGetLocalVector(networkdm,&localX);

379:   VecSet(X,0.0);
380:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
381:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

383:   VecGetArray(localX,&xarr);

385:   for (v = vStart; v < vEnd; v++) {
386:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
387:     if (ghostvtex) continue;

389:     DMNetworkGetVariableOffset(networkdm,v,&offset);
390:     DMNetworkGetNumComponents(networkdm,v,&numComps);
391:     for (j=0; j < numComps; j++) {
392:       DMNetworkGetComponent(networkdm,v,j,&key,&component);
393:       if (key == 1) {
394:         bus = (Bus*)(component);

396:         xarr[offset]   = bus->vr;
397:         xarr[offset+1] = bus->vi;

399:         Vr = bus->vr;
400:         Vi = bus->vi;
401:       } else if(key == 2) {
402:         gen = (Gen*)(component);

404:         Vm  = PetscSqrtScalar(Vr*Vr + Vi*Vi);
405:         Vm2 = Vm*Vm;
406:         /* Real part of gen current */
407:         IGr = (Vr*gen->PG + Vi*gen->QG)/Vm2;
408:         /* Imaginary part of gen current */
409:         IGi = (Vi*gen->PG - Vr*gen->QG)/Vm2;

411:         /* Machine angle */
412:         delta = atan2(Vi+gen->Xq*IGr,Vr-gen->Xq*IGi);
413:         theta = PETSC_PI/2.0 - delta;

415:         /* d-axis stator current */
416:         Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta);

418:         /* q-axis stator current */
419:         Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta);

421:         Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta);
422:         Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta);

424:         /* d-axis transient EMF */
425:         Edp = Vd + gen->Rs*Id - gen->Xqp*Iq;

427:         /* q-axis transient EMF */
428:         Eqp = Vq + gen->Rs*Iq + gen->Xdp*Id;

430:         gen->TM = gen->PG;
431:         idx     = offset+2;

433:         xarr[idx]   = Eqp;
434:         xarr[idx+1] = Edp;
435:         xarr[idx+2] = delta;
436:         xarr[idx+3] = W_S;
437:         xarr[idx+4] = Id;
438:         xarr[idx+5] = Iq;

440:         /* Exciter */
441:         Efd = Eqp + (gen->Xd - gen->Xdp)*Id;
442:         SE  = gen->k1*PetscExpScalar(gen->k2*Efd);
443:         VR  = gen->KE*Efd + SE;
444:         RF  = gen->KF*Efd/gen->TF;

446:         xarr[idx+6] = Efd;
447:         xarr[idx+7] = RF;
448:         xarr[idx+8] = VR;

450:         gen->Vref = Vm + (VR/gen->KA);
451:       }
452:     }
453:   }
454:   VecRestoreArray(localX,&xarr);
455:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
456:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
457:   DMRestoreLocalVector(networkdm,&localX);
458:   return(0);
459:  }

461:  /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
462: PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr,PetscScalar *Fi)
463: {
465:   *Fr =  Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta);
466:   *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta);
467:   return(0);
468: }

470: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
471: PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd,PetscScalar *Fq)
472: {
474:   *Fd =  Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta);
475:   *Fq =  Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta);
476:   return(0);
477: }

479: /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
480: PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx *user)
481: {
482:   PetscErrorCode                     ierr;
483:   DM                                 networkdm;
484:   Vec                                localX,localXdot,localF;
485:   PetscInt                           vfrom,vto,offsetfrom,offsetto;
486:   PetscInt                           v,vStart,vEnd,e;
487:   PetscScalar                        *farr;
488:   PetscScalar                        Vd,Vq,SE;
489:   const PetscScalar                  *xarr,*xdotarr;
490:   void*                              component;

493:   VecSet(F,0.0);

495:   TSGetDM(ts,&networkdm);
496:   DMGetLocalVector(networkdm,&localF);
497:   DMGetLocalVector(networkdm,&localX);
498:   DMGetLocalVector(networkdm,&localXdot);
499:   VecSet(localF,0.0);

501:   /* update ghost values of localX and localXdot */
502:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
503:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

505:   DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
506:   DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);

508:   VecGetArrayRead(localX,&xarr);
509:   VecGetArrayRead(localXdot,&xdotarr);
510:   VecGetArray(localF,&farr);

512:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

514:   for (v=vStart; v < vEnd; v++) {
515:     PetscInt     i,j,offset,key;
516:     PetscScalar  Vr, Vi;
517:     Bus          *bus;
518:     Gen          *gen;
519:     Load         *load;
520:     PetscBool    ghostvtex;
521:     PetscInt     numComps;
522:     PetscScalar  Yffr,Yffi; /* Real and imaginary fault admittances */
523:     PetscScalar  Vm,Vm2,Vm0;
524:     PetscScalar  Vr0=0,Vi0=0;
525:     PetscScalar  PD,QD;

527:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
528:     DMNetworkGetNumComponents(networkdm,v,&numComps);
529:     DMNetworkGetVariableOffset(networkdm,v,&offset);

531:     for (j = 0; j < numComps; j++) {
532:       DMNetworkGetComponent(networkdm,v,j,&key,&component);
533:       if (key == 1) {
534:         PetscInt       nconnedges;
535:         const PetscInt *connedges;

537:         bus = (Bus*)(component);
538:         if (!ghostvtex) {
539:           Vr   = xarr[offset];
540:           Vi   = xarr[offset+1];

542:           Yffr = bus->yff[1];
543:           Yffi = bus->yff[0];

545:           if (user->alg_flg){
546:             Yffr += user->ybusfault[bus->id*2+1];
547:             Yffi += user->ybusfault[bus->id*2];
548:           }
549:           Vr0 = bus->vr;
550:           Vi0 = bus->vi;

552:           /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
553:            The generator current injection, IG, and load current injection, ID are added later
554:            */
555:           farr[offset] +=  Yffi*Vr + Yffr*Vi; /* imaginary current due to diagonal elements */
556:           farr[offset+1] += Yffr*Vr - Yffi*Vi; /* real current due to diagonal elements */
557:         }

559:         DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);

561:         for (i=0; i < nconnedges; i++) {
562:           Branch         *branch;
563:           PetscInt       keye;
564:           PetscScalar    Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
565:           const PetscInt *cone;

567:           e = connedges[i];
568:           DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch);

570:           Yfti = branch->yft[0];
571:           Yftr = branch->yft[1];

573:           DMNetworkGetConnectedVertices(networkdm,e,&cone);

575:           vfrom = cone[0];
576:           vto   = cone[1];

578:           DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
579:           DMNetworkGetVariableOffset(networkdm,vto,&offsetto);

581:           /* From bus and to bus real and imaginary voltages */
582:           Vfr     = xarr[offsetfrom];
583:           Vfi     = xarr[offsetfrom+1];
584:           Vtr            = xarr[offsetto];
585:           Vti     = xarr[offsetto+1];

587:           if (vfrom == v) {
588:             farr[offsetfrom]   += Yftr*Vti + Yfti*Vtr;
589:             farr[offsetfrom+1] += Yftr*Vtr - Yfti*Vti;
590:           } else {
591:             farr[offsetto]   += Yftr*Vfi + Yfti*Vfr;
592:             farr[offsetto+1] += Yftr*Vfr - Yfti*Vfi;
593:           }
594:         }
595:       } else if (key == 2){
596:         if (!ghostvtex) {
597:           PetscScalar    Eqp,Edp,delta,w; /* Generator variables */
598:           PetscScalar    Efd,RF,VR; /* Exciter variables */
599:           PetscScalar    Id,Iq;  /* Generator dq axis currents */
600:           PetscScalar    IGr,IGi,Zdq_inv[4],det;
601:           PetscScalar    Xd,Xdp,Td0p,Xq,Xqp,Tq0p,TM,D,M,Rs; /* Generator parameters */
602:           PetscScalar    k1,k2,KE,TE,TF,KA,KF,Vref,TA; /* Generator parameters */
603:           PetscInt       idx;

605:           gen = (Gen*)(component);
606:           idx = offset + 2;

608:           /* Generator state variables */
609:           Eqp   = xarr[idx];
610:           Edp   = xarr[idx+1];
611:           delta = xarr[idx+2];
612:           w     = xarr[idx+3];
613:           Id    = xarr[idx+4];
614:           Iq    = xarr[idx+5];
615:           Efd   = xarr[idx+6];
616:           RF    = xarr[idx+7];
617:           VR    = xarr[idx+8];

619:           /* Generator parameters */
620:           Xd   = gen->Xd;
621:           Xdp  = gen->Xdp;
622:           Td0p = gen->Td0p;
623:           Xq   = gen->Xq;
624:           Xqp  = gen->Xqp;
625:           Tq0p = gen->Tq0p;
626:           TM   = gen->TM;
627:           D    = gen->D;
628:           M    = gen->M;
629:           Rs   = gen->Rs;
630:           k1   = gen->k1;
631:           k2   = gen->k2;
632:           KE   = gen->KE;
633:           TE   = gen->TE;
634:           TF   = gen->TF;
635:           KA   = gen->KA;
636:           KF   = gen->KF;
637:           Vref = gen->Vref;
638:           TA   = gen->TA;

640:           /* Generator differential equations */
641:           farr[idx]   = (Eqp + (Xd - Xdp)*Id - Efd)/Td0p + xdotarr[idx];
642:           farr[idx+1] = (Edp - (Xq - Xqp)*Iq)/Tq0p  + xdotarr[idx+1];
643:           farr[idx+2] = -w + W_S + xdotarr[idx+2];
644:           farr[idx+3] = (-TM + Edp*Id + Eqp*Iq + (Xqp - Xdp)*Id*Iq + D*(w - W_S))/M  + xdotarr[idx+3];

646:           Vr = xarr[offset]; /* Real part of generator terminal voltage */
647:           Vi = xarr[offset+1]; /* Imaginary part of the generator terminal voltage */

649:           ri2dq(Vr,Vi,delta,&Vd,&Vq);

651:           /* Algebraic equations for stator currents */
652:           det = Rs*Rs + Xdp*Xqp;

654:           Zdq_inv[0] = Rs/det;
655:           Zdq_inv[1] = Xqp/det;
656:           Zdq_inv[2] = -Xdp/det;
657:           Zdq_inv[3] = Rs/det;

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

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

665:           farr[offset]   -= IGi;
666:           farr[offset+1] -= IGr;

668:           Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);
669:           SE = k1*PetscExpScalar(k2*Efd);

671:           /* Exciter differential equations */
672:           farr[idx+6] = (KE*Efd + SE - VR)/TE + xdotarr[idx+6];
673:           farr[idx+7] = (RF - KF*Efd/TF)/TF + xdotarr[idx+7];
674:           farr[idx+8] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[idx+8];
675:         }
676:       } else if (key ==3){
677:         if (!ghostvtex) {
678:           PetscInt    k;
679:           PetscInt    ld_nsegsp;
680:           PetscInt    ld_nsegsq;
681:           PetscScalar *ld_alphap;
682:           PetscScalar *ld_betap,*ld_alphaq,*ld_betaq,PD0, QD0, IDr,IDi;

684:           load = (Load*)(component);

686:           /* Load Parameters */
687:           ld_nsegsp = load->ld_nsegsp;
688:           ld_alphap = load->ld_alphap;
689:           ld_betap  = load->ld_betap;
690:           ld_nsegsq = load->ld_nsegsq;
691:           ld_alphaq = load->ld_alphaq;
692:           ld_betaq  = load->ld_betaq;
693:           PD0       = load->PD0;
694:           QD0       = load->QD0;

696:           Vr  = xarr[offset]; /* Real part of generator terminal voltage */
697:           Vi  = xarr[offset+1]; /* Imaginary part of the generator terminal voltage */
698:           Vm  = PetscSqrtScalar(Vr*Vr + Vi*Vi);
699:           Vm2 = Vm*Vm;
700:           Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0);
701:           PD  = QD = 0.0;
702:           for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]);
703:           for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]);

705:           /* Load currents */
706:           IDr = (PD*Vr + QD*Vi)/Vm2;
707:           IDi = (-QD*Vr + PD*Vi)/Vm2;

709:           farr[offset]   += IDi;
710:           farr[offset+1] += IDr;
711:         }
712:       }
713:     }
714:   }

716:   VecRestoreArrayRead(localX,&xarr);
717:   VecRestoreArrayRead(localXdot,&xdotarr);
718:   VecRestoreArray(localF,&farr);
719:   DMRestoreLocalVector(networkdm,&localX);
720:   DMRestoreLocalVector(networkdm,&localXdot);

722:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
723:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
724:   DMRestoreLocalVector(networkdm,&localF);
725:   return(0);
726: }

728: /* This function is used for solving the algebraic system only during fault on and
729:    off times. It computes the entire F and then zeros out the part corresponding to
730:    differential equations
731:  F = [0;g(y)];
732: */
733: PetscErrorCode AlgFunction (SNES snes, Vec X, Vec F, void *ctx)
734: {
736:   DM             networkdm;
737:   Vec            localX,localF;
738:   PetscInt       vfrom,vto,offsetfrom,offsetto;
739:   PetscInt       v,vStart,vEnd,e;
740:   PetscScalar    *farr;
741:   Userctx        *user=(Userctx*)ctx;
742:   const PetscScalar *xarr;
743:   void*          component;

746:   VecSet(F,0.0);
747:   SNESGetDM(snes,&networkdm);
748:   DMGetLocalVector(networkdm,&localF);
749:   DMGetLocalVector(networkdm,&localX);
750:   VecSet(localF,0.0);

752:   /* update ghost values of locaX and locaXdot */
753:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
754:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

756:   VecGetArrayRead(localX,&xarr);
757:   VecGetArray(localF,&farr);

759:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

761:   for (v=vStart; v < vEnd; v++) {
762:     PetscInt      i,j,offset,key,numComps;
763:     PetscScalar   Vr, Vi, Yffr, Yffi, Vm, Vm2, Vm0, Vr0=0, Vi0=0, PD, QD;
764:     Bus           *bus;
765:     Gen           *gen;
766:     Load          *load;
767:     PetscBool     ghostvtex;

769:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
770:     DMNetworkGetNumComponents(networkdm,v,&numComps);
771:     DMNetworkGetVariableOffset(networkdm,v,&offset);

773:     for (j = 0; j < numComps; j++) {
774:       DMNetworkGetComponent(networkdm,v,j,&key,&component);
775:       if (key == 1) {
776:         PetscInt       nconnedges;
777:         const PetscInt *connedges;

779:         bus = (Bus*)(component);
780:         if (!ghostvtex) {
781:           Vr = xarr[offset];
782:           Vi = xarr[offset+1];

784:           Yffr = bus->yff[1];
785:           Yffi = bus->yff[0];
786:           if (user->alg_flg){
787:             Yffr += user->ybusfault[bus->id*2+1];
788:             Yffi += user->ybusfault[bus->id*2];
789:           }
790:           Vr0 = bus->vr;
791:           Vi0 = bus->vi;

793:           farr[offset]   += Yffi*Vr + Yffr*Vi;
794:           farr[offset+1] += Yffr*Vr - Yffi*Vi;
795:         }
796:         DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);

798:         for (i=0; i < nconnedges; i++) {
799:           Branch         *branch;
800:           PetscInt       keye;
801:           PetscScalar    Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
802:           const PetscInt *cone;

804:           e = connedges[i];
805:           DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch);

807:           Yfti = branch->yft[0];
808:           Yftr = branch->yft[1];

810:           DMNetworkGetConnectedVertices(networkdm,e,&cone);
811:           vfrom = cone[0];
812:           vto   = cone[1];

814:           DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
815:           DMNetworkGetVariableOffset(networkdm,vto,&offsetto);