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