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: PetscFunctionBegin;
102: *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
103: *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
108: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
109: {
110: PetscFunctionBegin;
111: *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
112: *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
117: {
118: Vec Xgen, Xnet;
119: PetscScalar *xgen, *xnet;
120: PetscInt i, idx = 0;
121: PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
122: PetscScalar Eqp, Edp, delta;
123: PetscScalar Efd, RF, VR; /* Exciter variables */
124: PetscScalar Id, Iq; /* Generator dq axis currents */
125: PetscScalar theta, Vd, Vq, SE;
127: PetscFunctionBegin;
128: M[0] = 2 * H[0] / w_s;
129: M[1] = 2 * H[1] / w_s;
130: M[2] = 2 * H[2] / w_s;
131: D[0] = 0.1 * M[0];
132: D[1] = 0.1 * M[1];
133: D[2] = 0.1 * M[2];
135: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
137: /* Network subsystem initialization */
138: PetscCall(VecCopy(user->V0, Xnet));
140: /* Generator subsystem initialization */
141: PetscCall(VecGetArray(Xgen, &xgen));
142: PetscCall(VecGetArray(Xnet, &xnet));
144: for (i = 0; i < ngen; i++) {
145: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
146: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
147: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
148: Vm2 = Vm * Vm;
149: IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
150: IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;
152: delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
154: theta = PETSC_PI / 2.0 - delta;
156: Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
157: Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
159: Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
160: Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
162: Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
163: Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
165: TM[i] = PG[i];
167: /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
168: xgen[idx] = Eqp;
169: xgen[idx + 1] = Edp;
170: xgen[idx + 2] = delta;
171: xgen[idx + 3] = w_s;
173: idx = idx + 4;
175: xgen[idx] = Id;
176: xgen[idx + 1] = Iq;
178: idx = idx + 2;
180: /* Exciter */
181: Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
182: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
183: VR = KE[i] * Efd + SE;
184: RF = KF[i] * Efd / TF[i];
186: xgen[idx] = Efd;
187: xgen[idx + 1] = RF;
188: xgen[idx + 2] = VR;
190: Vref[i] = Vm + (VR / KA[i]);
192: idx = idx + 3;
193: }
195: PetscCall(VecRestoreArray(Xgen, &xgen));
196: PetscCall(VecRestoreArray(Xnet, &xnet));
198: /* PetscCall(VecView(Xgen,0)); */
199: PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
200: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /* Computes F = [-f(x,y);g(x,y)] */
205: PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
206: {
207: Vec Xgen, Xnet, Fgen, Fnet;
208: PetscScalar *xgen, *xnet, *fgen, *fnet;
209: PetscInt i, idx = 0;
210: PetscScalar Vr, Vi, Vm, Vm2;
211: PetscScalar Eqp, Edp, delta, w; /* Generator variables */
212: PetscScalar Efd, RF, VR; /* Exciter variables */
213: PetscScalar Id, Iq; /* Generator dq axis currents */
214: PetscScalar Vd, Vq, SE;
215: PetscScalar IGr, IGi, IDr, IDi;
216: PetscScalar Zdq_inv[4], det;
217: PetscScalar PD, QD, Vm0, *v0;
218: PetscInt k;
220: PetscFunctionBegin;
221: PetscCall(VecZeroEntries(F));
222: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
223: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
224: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
225: PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));
227: /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
228: The generator current injection, IG, and load current injection, ID are added later
229: */
230: /* Note that the values in Ybus are stored assuming the imaginary current balance
231: equation is ordered first followed by real current balance equation for each bus.
232: Thus imaginary current contribution goes in location 2*i, and
233: real current contribution in 2*i+1
234: */
235: PetscCall(MatMult(user->Ybus, Xnet, Fnet));
237: PetscCall(VecGetArray(Xgen, &xgen));
238: PetscCall(VecGetArray(Xnet, &xnet));
239: PetscCall(VecGetArray(Fgen, &fgen));
240: PetscCall(VecGetArray(Fnet, &fnet));
242: /* Generator subsystem */
243: for (i = 0; i < ngen; i++) {
244: Eqp = xgen[idx];
245: Edp = xgen[idx + 1];
246: delta = xgen[idx + 2];
247: w = xgen[idx + 3];
248: Id = xgen[idx + 4];
249: Iq = xgen[idx + 5];
250: Efd = xgen[idx + 6];
251: RF = xgen[idx + 7];
252: VR = xgen[idx + 8];
254: /* Generator differential equations */
255: fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
256: fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
257: fgen[idx + 2] = -w + w_s;
258: fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
260: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
261: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
263: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
264: /* Algebraic equations for stator currents */
265: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
267: Zdq_inv[0] = Rs[i] / det;
268: Zdq_inv[1] = Xqp[i] / det;
269: Zdq_inv[2] = -Xdp[i] / det;
270: Zdq_inv[3] = Rs[i] / det;
272: fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
273: fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
275: /* Add generator current injection to network */
276: PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
278: fnet[2 * gbus[i]] -= IGi;
279: fnet[2 * gbus[i] + 1] -= IGr;
281: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
283: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
285: /* Exciter differential equations */
286: fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
287: fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
288: fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
290: idx = idx + 9;
291: }
293: PetscCall(VecGetArray(user->V0, &v0));
294: for (i = 0; i < nload; i++) {
295: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
296: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
297: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
298: Vm2 = Vm * Vm;
299: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
300: PD = QD = 0.0;
301: for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
302: for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
304: /* Load currents */
305: IDr = (PD * Vr + QD * Vi) / Vm2;
306: IDi = (-QD * Vr + PD * Vi) / Vm2;
308: fnet[2 * lbus[i]] += IDi;
309: fnet[2 * lbus[i] + 1] += IDr;
310: }
311: PetscCall(VecRestoreArray(user->V0, &v0));
313: PetscCall(VecRestoreArray(Xgen, &xgen));
314: PetscCall(VecRestoreArray(Xnet, &xnet));
315: PetscCall(VecRestoreArray(Fgen, &fgen));
316: PetscCall(VecRestoreArray(Fnet, &fnet));
318: PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
319: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
320: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /* \dot{x} - f(x,y)
325: g(x,y) = 0
326: */
327: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
328: {
329: SNES snes;
330: PetscScalar *f;
331: const PetscScalar *xdot;
332: PetscInt i;
334: PetscFunctionBegin;
335: user->t = t;
337: PetscCall(TSGetSNES(ts, &snes));
338: PetscCall(ResidualFunction(snes, X, F, user));
339: PetscCall(VecGetArray(F, &f));
340: PetscCall(VecGetArrayRead(Xdot, &xdot));
341: for (i = 0; i < ngen; i++) {
342: f[9 * i] += xdot[9 * i];
343: f[9 * i + 1] += xdot[9 * i + 1];
344: f[9 * i + 2] += xdot[9 * i + 2];
345: f[9 * i + 3] += xdot[9 * i + 3];
346: f[9 * i + 6] += xdot[9 * i + 6];
347: f[9 * i + 7] += xdot[9 * i + 7];
348: f[9 * i + 8] += xdot[9 * i + 8];
349: }
350: PetscCall(VecRestoreArray(F, &f));
351: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /* This function is used for solving the algebraic system only during fault on and
356: off times. It computes the entire F and then zeros out the part corresponding to
357: differential equations
358: F = [0;g(y)];
359: */
360: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
361: {
362: Userctx *user = (Userctx *)ctx;
363: PetscScalar *f;
364: PetscInt i;
366: PetscFunctionBegin;
367: PetscCall(ResidualFunction(snes, X, F, user));
368: PetscCall(VecGetArray(F, &f));
369: for (i = 0; i < ngen; i++) {
370: f[9 * i] = 0;
371: f[9 * i + 1] = 0;
372: f[9 * i + 2] = 0;
373: f[9 * i + 3] = 0;
374: f[9 * i + 6] = 0;
375: f[9 * i + 7] = 0;
376: f[9 * i + 8] = 0;
377: }
378: PetscCall(VecRestoreArray(F, &f));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
383: {
384: PetscInt *d_nnz;
385: PetscInt i, idx = 0, start = 0;
386: PetscInt ncols;
388: PetscFunctionBegin;
389: PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
390: for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
391: /* Generator subsystem */
392: for (i = 0; i < ngen; i++) {
393: d_nnz[idx] += 3;
394: d_nnz[idx + 1] += 2;
395: d_nnz[idx + 2] += 2;
396: d_nnz[idx + 3] += 5;
397: d_nnz[idx + 4] += 6;
398: d_nnz[idx + 5] += 6;
400: d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
401: d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
403: d_nnz[idx + 6] += 2;
404: d_nnz[idx + 7] += 2;
405: d_nnz[idx + 8] += 5;
407: idx = idx + 9;
408: }
410: start = user->neqs_gen;
412: for (i = 0; i < nbus; i++) {
413: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
414: d_nnz[start + 2 * i] += ncols;
415: d_nnz[start + 2 * i + 1] += ncols;
416: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
417: }
419: PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
421: PetscCall(PetscFree(d_nnz));
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*
426: J = [-df_dx, -df_dy
427: dg_dx, dg_dy]
428: */
429: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, void *ctx)
430: {
431: Userctx *user = (Userctx *)ctx;
432: Vec Xgen, Xnet;
433: PetscScalar *xgen, *xnet;
434: PetscInt i, idx = 0;
435: PetscScalar Vr, Vi, Vm, Vm2;
436: PetscScalar Eqp, Edp, delta; /* Generator variables */
437: PetscScalar Efd; /* Exciter variables */
438: PetscScalar Id, Iq; /* Generator dq axis currents */
439: PetscScalar Vd, Vq;
440: PetscScalar val[10];
441: PetscInt row[2], col[10];
442: PetscInt net_start = user->neqs_gen;
443: PetscScalar Zdq_inv[4], det;
444: PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
445: PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
446: PetscScalar dSE_dEfd;
447: PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
448: PetscInt ncols;
449: const PetscInt *cols;
450: const PetscScalar *yvals;
451: PetscInt k;
452: PetscScalar PD, QD, Vm0, *v0, Vm4;
453: PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
454: PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
456: PetscFunctionBegin;
457: PetscCall(MatZeroEntries(B));
458: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
459: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
461: PetscCall(VecGetArray(Xgen, &xgen));
462: PetscCall(VecGetArray(Xnet, &xnet));
464: /* Generator subsystem */
465: for (i = 0; i < ngen; i++) {
466: Eqp = xgen[idx];
467: Edp = xgen[idx + 1];
468: delta = xgen[idx + 2];
469: Id = xgen[idx + 4];
470: Iq = xgen[idx + 5];
471: Efd = xgen[idx + 6];
473: /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
474: row[0] = idx;
475: col[0] = idx;
476: col[1] = idx + 4;
477: col[2] = idx + 6;
478: val[0] = 1 / Td0p[i];
479: val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
480: val[2] = -1 / Td0p[i];
482: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
484: /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
485: row[0] = idx + 1;
486: col[0] = idx + 1;
487: col[1] = idx + 5;
488: val[0] = 1 / Tq0p[i];
489: val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
490: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
492: /* fgen[idx+2] = - w + w_s; */
493: row[0] = idx + 2;
494: col[0] = idx + 2;
495: col[1] = idx + 3;
496: val[0] = 0;
497: val[1] = -1;
498: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
500: /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
501: row[0] = idx + 3;
502: col[0] = idx;
503: col[1] = idx + 1;
504: col[2] = idx + 3;
505: col[3] = idx + 4;
506: col[4] = idx + 5;
507: val[0] = Iq / M[i];
508: val[1] = Id / M[i];
509: val[2] = D[i] / M[i];
510: val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
511: val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
512: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
514: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
515: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
516: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
518: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
520: Zdq_inv[0] = Rs[i] / det;
521: Zdq_inv[1] = Xqp[i] / det;
522: Zdq_inv[2] = -Xdp[i] / det;
523: Zdq_inv[3] = Rs[i] / det;
525: dVd_dVr = PetscSinScalar(delta);
526: dVd_dVi = -PetscCosScalar(delta);
527: dVq_dVr = PetscCosScalar(delta);
528: dVq_dVi = PetscSinScalar(delta);
529: dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
530: dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
532: /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
533: row[0] = idx + 4;
534: col[0] = idx;
535: col[1] = idx + 1;
536: col[2] = idx + 2;
537: val[0] = -Zdq_inv[1];
538: val[1] = -Zdq_inv[0];
539: val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
540: col[3] = idx + 4;
541: col[4] = net_start + 2 * gbus[i];
542: col[5] = net_start + 2 * gbus[i] + 1;
543: val[3] = 1;
544: val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
545: val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
546: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
548: /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
549: row[0] = idx + 5;
550: col[0] = idx;
551: col[1] = idx + 1;
552: col[2] = idx + 2;
553: val[0] = -Zdq_inv[3];
554: val[1] = -Zdq_inv[2];
555: val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
556: col[3] = idx + 5;
557: col[4] = net_start + 2 * gbus[i];
558: col[5] = net_start + 2 * gbus[i] + 1;
559: val[3] = 1;
560: val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
561: val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
562: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
564: dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
565: dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
566: dIGr_dId = PetscSinScalar(delta);
567: dIGr_dIq = PetscCosScalar(delta);
568: dIGi_dId = -PetscCosScalar(delta);
569: dIGi_dIq = PetscSinScalar(delta);
571: /* fnet[2*gbus[i]] -= IGi; */
572: row[0] = net_start + 2 * gbus[i];
573: col[0] = idx + 2;
574: col[1] = idx + 4;
575: col[2] = idx + 5;
576: val[0] = -dIGi_ddelta;
577: val[1] = -dIGi_dId;
578: val[2] = -dIGi_dIq;
579: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
581: /* fnet[2*gbus[i]+1] -= IGr; */
582: row[0] = net_start + 2 * gbus[i] + 1;
583: col[0] = idx + 2;
584: col[1] = idx + 4;
585: col[2] = idx + 5;
586: val[0] = -dIGr_ddelta;
587: val[1] = -dIGr_dId;
588: val[2] = -dIGr_dIq;
589: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
591: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
593: /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
594: /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
596: dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
598: row[0] = idx + 6;
599: col[0] = idx + 6;
600: col[1] = idx + 8;
601: val[0] = (KE[i] + dSE_dEfd) / TE[i];
602: val[1] = -1 / TE[i];
603: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
605: /* Exciter differential equations */
607: /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
608: row[0] = idx + 7;
609: col[0] = idx + 6;
610: col[1] = idx + 7;
611: val[0] = (-KF[i] / TF[i]) / TF[i];
612: val[1] = 1 / TF[i];
613: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
615: /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
616: /* Vm = (Vd^2 + Vq^2)^0.5; */
617: dVm_dVd = Vd / Vm;
618: dVm_dVq = Vq / Vm;
619: dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
620: dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
621: row[0] = idx + 8;
622: col[0] = idx + 6;
623: col[1] = idx + 7;
624: col[2] = idx + 8;
625: val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
626: val[1] = -KA[i] / TA[i];
627: val[2] = 1 / TA[i];
628: col[3] = net_start + 2 * gbus[i];
629: col[4] = net_start + 2 * gbus[i] + 1;
630: val[3] = KA[i] * dVm_dVr / TA[i];
631: val[4] = KA[i] * dVm_dVi / TA[i];
632: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
633: idx = idx + 9;
634: }
636: for (i = 0; i < nbus; i++) {
637: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
638: row[0] = net_start + 2 * i;
639: for (k = 0; k < ncols; k++) {
640: col[k] = net_start + cols[k];
641: val[k] = yvals[k];
642: }
643: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
644: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
646: PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
647: row[0] = net_start + 2 * i + 1;
648: for (k = 0; k < ncols; k++) {
649: col[k] = net_start + cols[k];
650: val[k] = yvals[k];
651: }
652: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
653: PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
654: }
656: PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
657: PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));
659: PetscCall(VecGetArray(user->V0, &v0));
660: for (i = 0; i < nload; i++) {
661: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
662: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
663: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
664: Vm2 = Vm * Vm;
665: Vm4 = Vm2 * Vm2;
666: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
667: PD = QD = 0.0;
668: dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
669: for (k = 0; k < ld_nsegsp[i]; k++) {
670: PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
671: dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vr * PetscPowScalar(Vm, (ld_betap[k] - 2));
672: dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vi * PetscPowScalar(Vm, (ld_betap[k] - 2));
673: }
674: for (k = 0; k < ld_nsegsq[i]; k++) {
675: QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
676: dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vr * PetscPowScalar(Vm, (ld_betaq[k] - 2));
677: dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vi * PetscPowScalar(Vm, (ld_betaq[k] - 2));
678: }
680: /* IDr = (PD*Vr + QD*Vi)/Vm2; */
681: /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
683: dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
684: dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
686: dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
687: dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
689: /* fnet[2*lbus[i]] += IDi; */
690: row[0] = net_start + 2 * lbus[i];
691: col[0] = net_start + 2 * lbus[i];
692: col[1] = net_start + 2 * lbus[i] + 1;
693: val[0] = dIDi_dVr;
694: val[1] = dIDi_dVi;
695: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
696: /* fnet[2*lbus[i]+1] += IDr; */
697: row[0] = net_start + 2 * lbus[i] + 1;
698: col[0] = net_start + 2 * lbus[i];
699: col[1] = net_start + 2 * lbus[i] + 1;
700: val[0] = dIDr_dVr;
701: val[1] = dIDr_dVi;
702: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
703: }
704: PetscCall(VecRestoreArray(user->V0, &v0));
706: PetscCall(VecRestoreArray(Xgen, &xgen));
707: PetscCall(VecRestoreArray(Xnet, &xnet));
709: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
711: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
712: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: /*
717: J = [I, 0
718: dg_dx, dg_dy]
719: */
720: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, void *ctx)
721: {
722: Userctx *user = (Userctx *)ctx;
724: PetscFunctionBegin;
725: PetscCall(ResidualJacobian(snes, X, A, B, ctx));
726: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
727: PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: /*
732: J = [a*I-df_dx, -df_dy
733: dg_dx, dg_dy]
734: */
736: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
737: {
738: SNES snes;
739: PetscScalar atmp = (PetscScalar)a;
740: PetscInt i, row;
742: PetscFunctionBegin;
743: user->t = t;
745: PetscCall(TSGetSNES(ts, &snes));
746: PetscCall(ResidualJacobian(snes, X, A, B, user));
747: for (i = 0; i < ngen; i++) {
748: row = 9 * i;
749: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
750: row = 9 * i + 1;
751: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
752: row = 9 * i + 2;
753: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
754: row = 9 * i + 3;
755: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
756: row = 9 * i + 6;
757: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
758: row = 9 * i + 7;
759: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
760: row = 9 * i + 8;
761: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
762: }
763: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
764: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
769: {
770: PetscScalar *r;
771: const PetscScalar *u;
772: PetscInt idx;
773: Vec Xgen, Xnet;
774: PetscScalar *xgen;
775: PetscInt i;
777: PetscFunctionBegin;
778: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
779: PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
781: PetscCall(VecGetArray(Xgen, &xgen));
783: PetscCall(VecGetArrayRead(U, &u));
784: PetscCall(VecGetArray(R, &r));
785: r[0] = 0.;
787: idx = 0;
788: for (i = 0; i < ngen; i++) {
789: 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);
790: idx += 9;
791: }
792: PetscCall(VecRestoreArray(R, &r));
793: PetscCall(VecRestoreArrayRead(U, &u));
794: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: static PetscErrorCode MonitorUpdateQ(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx0)
799: {
800: Vec C, *Y;
801: PetscInt Nr;
802: PetscReal h, theta;
803: Userctx *ctx = (Userctx *)ctx0;
805: PetscFunctionBegin;
806: theta = 0.5;
807: PetscCall(TSGetStages(ts, &Nr, &Y));
808: PetscCall(TSGetTimeStep(ts, &h));
809: PetscCall(VecDuplicate(ctx->vec_q, &C));
810: /* compute integrals */
811: if (stepnum > 0) {
812: PetscCall(CostIntegrand(ts, time, X, C, ctx));
813: PetscCall(VecAXPY(ctx->vec_q, h * theta, C));
814: PetscCall(CostIntegrand(ts, time + h * theta, Y[0], C, ctx));
815: PetscCall(VecAXPY(ctx->vec_q, h * (1 - theta), C));
816: }
817: PetscCall(VecDestroy(&C));
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: int main(int argc, char **argv)
822: {
823: Userctx user;
824: Vec p;
825: PetscScalar *x_ptr;
826: PetscMPIInt size;
827: PetscInt i;
828: KSP ksp;
829: PC pc;
830: PetscInt *idx2;
831: Tao tao;
832: Vec lowerb, upperb;
834: PetscFunctionBeginUser;
835: PetscFunctionBeginUser;
836: PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
837: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
838: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
840: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 1, &user.vec_q));
842: user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
843: user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
844: user.neqs_pgrid = user.neqs_gen + user.neqs_net;
846: /* Create indices for differential and algebraic equations */
847: PetscCall(PetscMalloc1(7 * ngen, &idx2));
848: for (i = 0; i < ngen; i++) {
849: idx2[7 * i] = 9 * i;
850: idx2[7 * i + 1] = 9 * i + 1;
851: idx2[7 * i + 2] = 9 * i + 2;
852: idx2[7 * i + 3] = 9 * i + 3;
853: idx2[7 * i + 4] = 9 * i + 6;
854: idx2[7 * i + 5] = 9 * i + 7;
855: idx2[7 * i + 6] = 9 * i + 8;
856: }
857: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
858: PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
859: PetscCall(PetscFree(idx2));
861: /* Set run time options */
862: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
863: {
864: user.tfaulton = 1.0;
865: user.tfaultoff = 1.2;
866: user.Rfault = 0.0001;
867: user.faultbus = 8;
868: PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
869: PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
870: PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
871: user.t0 = 0.0;
872: user.tmax = 1.5;
873: PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
874: PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
875: user.freq_u = 61.0;
876: user.freq_l = 59.0;
877: user.pow = 2;
878: PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
879: PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
880: PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
881: }
882: PetscOptionsEnd();
884: /* Create DMs for generator and network subsystems */
885: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
886: PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
887: PetscCall(DMSetFromOptions(user.dmgen));
888: PetscCall(DMSetUp(user.dmgen));
889: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
890: PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
891: PetscCall(DMSetFromOptions(user.dmnet));
892: PetscCall(DMSetUp(user.dmnet));
893: /* Create a composite DM packer and add the two DMs */
894: PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
895: PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
896: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
897: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));
899: /* Create TAO solver and set desired solution method */
900: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
901: PetscCall(TaoSetType(tao, TAOBLMVM));
902: /*
903: Optimization starts
904: */
905: /* Set initial solution guess */
906: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
907: PetscCall(VecGetArray(p, &x_ptr));
908: x_ptr[0] = PG[0];
909: x_ptr[1] = PG[1];
910: x_ptr[2] = PG[2];
911: PetscCall(VecRestoreArray(p, &x_ptr));
913: PetscCall(TaoSetSolution(tao, p));
914: /* Set routine for function and gradient evaluation */
915: PetscCall(TaoSetObjective(tao, FormFunction, (void *)&user));
916: PetscCall(TaoSetGradient(tao, NULL, TaoDefaultComputeGradient, (void *)&user));
918: /* Set bounds for the optimization */
919: PetscCall(VecDuplicate(p, &lowerb));
920: PetscCall(VecDuplicate(p, &upperb));
921: PetscCall(VecGetArray(lowerb, &x_ptr));
922: x_ptr[0] = 0.5;
923: x_ptr[1] = 0.5;
924: x_ptr[2] = 0.5;
925: PetscCall(VecRestoreArray(lowerb, &x_ptr));
926: PetscCall(VecGetArray(upperb, &x_ptr));
927: x_ptr[0] = 2.0;
928: x_ptr[1] = 2.0;
929: x_ptr[2] = 2.0;
930: PetscCall(VecRestoreArray(upperb, &x_ptr));
931: PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));
933: /* Check for any TAO command line options */
934: PetscCall(TaoSetFromOptions(tao));
935: PetscCall(TaoGetKSP(tao, &ksp));
936: if (ksp) {
937: PetscCall(KSPGetPC(ksp, &pc));
938: PetscCall(PCSetType(pc, PCNONE));
939: }
941: /* SOLVE THE APPLICATION */
942: PetscCall(TaoSolve(tao));
944: PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
945: /* Free TAO data structures */
946: PetscCall(TaoDestroy(&tao));
947: PetscCall(VecDestroy(&user.vec_q));
948: PetscCall(VecDestroy(&lowerb));
949: PetscCall(VecDestroy(&upperb));
950: PetscCall(VecDestroy(&p));
951: PetscCall(DMDestroy(&user.dmgen));
952: PetscCall(DMDestroy(&user.dmnet));
953: PetscCall(DMDestroy(&user.dmpgrid));
954: PetscCall(ISDestroy(&user.is_diff));
955: PetscCall(ISDestroy(&user.is_alg));
956: PetscCall(PetscFinalize());
957: return 0;
958: }
960: /* ------------------------------------------------------------------ */
961: /*
962: FormFunction - Evaluates the function and corresponding gradient.
964: Input Parameters:
965: tao - the Tao context
966: X - the input vector
967: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
969: Output Parameters:
970: f - the newly evaluated function
971: */
972: PetscErrorCode FormFunction(Tao tao, Vec P, PetscReal *f, void *ctx0)
973: {
974: TS ts;
975: SNES snes_alg;
976: Userctx *ctx = (Userctx *)ctx0;
977: Vec X;
978: Mat J;
979: /* sensitivity context */
980: PetscScalar *x_ptr;
981: PetscViewer Xview, Ybusview;
982: Vec F_alg;
983: Vec Xdot;
984: PetscInt row_loc, col_loc;
985: PetscScalar val;
987: PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
988: PG[0] = x_ptr[0];
989: PG[1] = x_ptr[1];
990: PG[2] = x_ptr[2];
991: PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
993: ctx->stepnum = 0;
995: PetscCall(VecZeroEntries(ctx->vec_q));
997: /* Read initial voltage vector and Ybus */
998: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
999: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));
1001: PetscCall(VecCreate(PETSC_COMM_WORLD, &ctx->V0));
1002: PetscCall(VecSetSizes(ctx->V0, PETSC_DECIDE, ctx->neqs_net));
1003: PetscCall(VecLoad(ctx->V0, Xview));
1005: PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx->Ybus));
1006: PetscCall(MatSetSizes(ctx->Ybus, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_net, ctx->neqs_net));
1007: PetscCall(MatSetType(ctx->Ybus, MATBAIJ));
1008: /* PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
1009: PetscCall(MatLoad(ctx->Ybus, Ybusview));
1011: PetscCall(PetscViewerDestroy(&Xview));
1012: PetscCall(PetscViewerDestroy(&Ybusview));
1014: PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));
1016: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
1017: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_pgrid, ctx->neqs_pgrid));
1018: PetscCall(MatSetFromOptions(J));
1019: PetscCall(PreallocateJacobian(J, ctx));
1021: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1022: Create timestepping solver context
1023: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1024: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1025: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1026: PetscCall(TSSetType(ts, TSCN));
1027: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, ctx));
1028: PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobianFn *)IJacobian, ctx));
1029: PetscCall(TSSetApplicationContext(ts, ctx));
1031: PetscCall(TSMonitorSet(ts, MonitorUpdateQ, ctx, NULL));
1032: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1033: Set initial conditions
1034: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1035: PetscCall(SetInitialGuess(X, ctx));
1037: PetscCall(VecDuplicate(X, &F_alg));
1038: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1039: PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
1040: PetscCall(MatZeroEntries(J));
1041: PetscCall(SNESSetJacobian(snes_alg, J, J, AlgJacobian, ctx));
1042: PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1043: PetscCall(SNESSetFromOptions(snes_alg));
1044: ctx->alg_flg = PETSC_TRUE;
1045: /* Solve the algebraic equations */
1046: PetscCall(SNESSolve(snes_alg, NULL, X));
1048: /* Just to set up the Jacobian structure */
1049: PetscCall(VecDuplicate(X, &Xdot));
1050: PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, J, J, ctx));
1051: PetscCall(VecDestroy(&Xdot));
1053: ctx->stepnum++;
1055: PetscCall(TSSetTimeStep(ts, 0.01));
1056: PetscCall(TSSetMaxTime(ts, ctx->tmax));
1057: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1058: PetscCall(TSSetFromOptions(ts));
1060: /* Prefault period */
1061: ctx->alg_flg = PETSC_FALSE;
1062: PetscCall(TSSetTime(ts, 0.0));
1063: PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
1064: PetscCall(TSSolve(ts, X));
1066: /* Create the nonlinear solver for solving the algebraic system */
1067: /* Note that although the algebraic system needs to be solved only for
1068: Idq and V, we reuse the entire system including xgen. The xgen
1069: variables are held constant by setting their residuals to 0 and
1070: putting a 1 on the Jacobian diagonal for xgen rows
1071: */
1072: PetscCall(MatZeroEntries(J));
1074: /* Apply disturbance - resistive fault at ctx->faultbus */
1075: /* This is done by adding shunt conductance to the diagonal location
1076: in the Ybus matrix */
1077: row_loc = 2 * ctx->faultbus;
1078: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1079: val = 1 / ctx->Rfault;
1080: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1081: row_loc = 2 * ctx->faultbus + 1;
1082: col_loc = 2 * ctx->faultbus; /* Location for G */
1083: val = 1 / ctx->Rfault;
1084: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1086: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1087: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1089: ctx->alg_flg = PETSC_TRUE;
1090: /* Solve the algebraic equations */
1091: PetscCall(SNESSolve(snes_alg, NULL, X));
1093: ctx->stepnum++;
1095: /* Disturbance period */
1096: ctx->alg_flg = PETSC_FALSE;
1097: PetscCall(TSSetTime(ts, ctx->tfaulton));
1098: PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
1099: PetscCall(TSSolve(ts, X));
1101: /* Remove the fault */
1102: row_loc = 2 * ctx->faultbus;
1103: col_loc = 2 * ctx->faultbus + 1;
1104: val = -1 / ctx->Rfault;
1105: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1106: row_loc = 2 * ctx->faultbus + 1;
1107: col_loc = 2 * ctx->faultbus;
1108: val = -1 / ctx->Rfault;
1109: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1111: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1112: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1114: PetscCall(MatZeroEntries(J));
1116: ctx->alg_flg = PETSC_TRUE;
1118: /* Solve the algebraic equations */
1119: PetscCall(SNESSolve(snes_alg, NULL, X));
1121: ctx->stepnum++;
1123: /* Post-disturbance period */
1124: ctx->alg_flg = PETSC_TRUE;
1125: PetscCall(TSSetTime(ts, ctx->tfaultoff));
1126: PetscCall(TSSetMaxTime(ts, ctx->tmax));
1127: PetscCall(TSSolve(ts, X));
1129: PetscCall(VecGetArray(ctx->vec_q, &x_ptr));
1130: *f = x_ptr[0];
1131: PetscCall(VecRestoreArray(ctx->vec_q, &x_ptr));
1133: PetscCall(MatDestroy(&ctx->Ybus));
1134: PetscCall(VecDestroy(&ctx->V0));
1135: PetscCall(SNESDestroy(&snes_alg));
1136: PetscCall(VecDestroy(&F_alg));
1137: PetscCall(MatDestroy(&J));
1138: PetscCall(VecDestroy(&X));
1139: PetscCall(TSDestroy(&ts));
1141: return 0;
1142: }
1144: /*TEST
1146: build:
1147: requires: double !complex !defined(USE_64BIT_INDICES)
1149: test:
1150: args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1151: localrunfiles: petscoptions X.bin Ybus.bin
1153: TEST*/