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;
219: PetscFunctionBegin;
220: PetscCall(VecZeroEntries(F));
221: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
222: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
223: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
224: PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));
226: /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
227: The generator current injection, IG, and load current injection, ID are added later
228: */
229: /* Note that the values in Ybus are stored assuming the imaginary current balance
230: equation is ordered first followed by real current balance equation for each bus.
231: Thus imaginary current contribution goes in location 2*i, and
232: real current contribution in 2*i+1
233: */
234: PetscCall(MatMult(user->Ybus, Xnet, Fnet));
236: PetscCall(VecGetArray(Xgen, &xgen));
237: PetscCall(VecGetArray(Xnet, &xnet));
238: PetscCall(VecGetArray(Fgen, &fgen));
239: PetscCall(VecGetArray(Fnet, &fnet));
241: /* Generator subsystem */
242: for (i = 0; i < ngen; i++) {
243: Eqp = xgen[idx];
244: Edp = xgen[idx + 1];
245: delta = xgen[idx + 2];
246: w = xgen[idx + 3];
247: Id = xgen[idx + 4];
248: Iq = xgen[idx + 5];
249: Efd = xgen[idx + 6];
250: RF = xgen[idx + 7];
251: VR = xgen[idx + 8];
253: /* Generator differential equations */
254: fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
255: fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
256: fgen[idx + 2] = -w + w_s;
257: fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
259: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
260: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
262: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
263: /* Algebraic equations for stator currents */
264: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
266: Zdq_inv[0] = Rs[i] / det;
267: Zdq_inv[1] = Xqp[i] / det;
268: Zdq_inv[2] = -Xdp[i] / det;
269: Zdq_inv[3] = Rs[i] / det;
271: fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
272: fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
274: /* Add generator current injection to network */
275: PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
277: fnet[2 * gbus[i]] -= IGi;
278: fnet[2 * gbus[i] + 1] -= IGr;
280: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
282: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
284: /* Exciter differential equations */
285: fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
286: fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
287: fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
289: idx = idx + 9;
290: }
292: PetscCall(VecGetArray(user->V0, &v0));
293: for (i = 0; i < nload; i++) {
294: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
295: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
296: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
297: Vm2 = Vm * Vm;
298: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
299: PD = QD = 0.0;
300: for (PetscInt k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
301: for (PetscInt k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
303: /* Load currents */
304: IDr = (PD * Vr + QD * Vi) / Vm2;
305: IDi = (-QD * Vr + PD * Vi) / Vm2;
307: fnet[2 * lbus[i]] += IDi;
308: fnet[2 * lbus[i] + 1] += IDr;
309: }
310: PetscCall(VecRestoreArray(user->V0, &v0));
312: PetscCall(VecRestoreArray(Xgen, &xgen));
313: PetscCall(VecRestoreArray(Xnet, &xnet));
314: PetscCall(VecRestoreArray(Fgen, &fgen));
315: PetscCall(VecRestoreArray(Fnet, &fnet));
317: PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
318: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
319: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /* \dot{x} - f(x,y)
324: g(x,y) = 0
325: */
326: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
327: {
328: SNES snes;
329: PetscScalar *f;
330: const PetscScalar *xdot;
331: PetscInt i;
333: PetscFunctionBegin;
334: user->t = t;
336: PetscCall(TSGetSNES(ts, &snes));
337: PetscCall(ResidualFunction(snes, X, F, user));
338: PetscCall(VecGetArray(F, &f));
339: PetscCall(VecGetArrayRead(Xdot, &xdot));
340: for (i = 0; i < ngen; i++) {
341: f[9 * i] += xdot[9 * i];
342: f[9 * i + 1] += xdot[9 * i + 1];
343: f[9 * i + 2] += xdot[9 * i + 2];
344: f[9 * i + 3] += xdot[9 * i + 3];
345: f[9 * i + 6] += xdot[9 * i + 6];
346: f[9 * i + 7] += xdot[9 * i + 7];
347: f[9 * i + 8] += xdot[9 * i + 8];
348: }
349: PetscCall(VecRestoreArray(F, &f));
350: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /* This function is used for solving the algebraic system only during fault on and
355: off times. It computes the entire F and then zeros out the part corresponding to
356: differential equations
357: F = [0;g(y)];
358: */
359: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, PetscCtx ctx)
360: {
361: Userctx *user = (Userctx *)ctx;
362: PetscScalar *f;
363: PetscInt i;
365: PetscFunctionBegin;
366: PetscCall(ResidualFunction(snes, X, F, user));
367: PetscCall(VecGetArray(F, &f));
368: for (i = 0; i < ngen; i++) {
369: f[9 * i] = 0;
370: f[9 * i + 1] = 0;
371: f[9 * i + 2] = 0;
372: f[9 * i + 3] = 0;
373: f[9 * i + 6] = 0;
374: f[9 * i + 7] = 0;
375: f[9 * i + 8] = 0;
376: }
377: PetscCall(VecRestoreArray(F, &f));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
382: {
383: PetscInt *d_nnz;
384: PetscInt i, idx = 0, start = 0;
385: PetscInt ncols;
387: PetscFunctionBegin;
388: PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
389: for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
390: /* Generator subsystem */
391: for (i = 0; i < ngen; i++) {
392: d_nnz[idx] += 3;
393: d_nnz[idx + 1] += 2;
394: d_nnz[idx + 2] += 2;
395: d_nnz[idx + 3] += 5;
396: d_nnz[idx + 4] += 6;
397: d_nnz[idx + 5] += 6;
399: d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
400: d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
402: d_nnz[idx + 6] += 2;
403: d_nnz[idx + 7] += 2;
404: d_nnz[idx + 8] += 5;
406: idx = idx + 9;
407: }
409: start = user->neqs_gen;
411: for (i = 0; i < nbus; i++) {
412: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
413: d_nnz[start + 2 * i] += ncols;
414: d_nnz[start + 2 * i + 1] += ncols;
415: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
416: }
418: PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
420: PetscCall(PetscFree(d_nnz));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*
425: J = [-df_dx, -df_dy
426: dg_dx, dg_dy]
427: */
428: PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, PetscCtx ctx)
429: {
430: Userctx *user = (Userctx *)ctx;
431: Vec Xgen, Xnet;
432: PetscScalar *xgen, *xnet;
433: PetscInt i, idx = 0;
434: PetscScalar Vr, Vi, Vm, Vm2;
435: PetscScalar Eqp, Edp, delta; /* Generator variables */
436: PetscScalar Efd; /* Exciter variables */
437: PetscScalar Id, Iq; /* Generator dq axis currents */
438: PetscScalar Vd, Vq;
439: PetscScalar val[10];
440: PetscInt row[2], col[10];
441: PetscInt net_start = user->neqs_gen;
442: PetscScalar Zdq_inv[4], det;
443: PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
444: PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
445: PetscScalar dSE_dEfd;
446: PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
447: PetscInt ncols;
448: const PetscInt *cols;
449: const PetscScalar *yvals;
450: PetscScalar PD, QD, Vm0, *v0, Vm4;
451: PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
452: PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
454: PetscFunctionBegin;
455: PetscCall(MatZeroEntries(B));
456: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
457: PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
459: PetscCall(VecGetArray(Xgen, &xgen));
460: PetscCall(VecGetArray(Xnet, &xnet));
462: /* Generator subsystem */
463: for (i = 0; i < ngen; i++) {
464: Eqp = xgen[idx];
465: Edp = xgen[idx + 1];
466: delta = xgen[idx + 2];
467: Id = xgen[idx + 4];
468: Iq = xgen[idx + 5];
469: Efd = xgen[idx + 6];
471: /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
472: row[0] = idx;
473: col[0] = idx;
474: col[1] = idx + 4;
475: col[2] = idx + 6;
476: val[0] = 1 / Td0p[i];
477: val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
478: val[2] = -1 / Td0p[i];
480: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
482: /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
483: row[0] = idx + 1;
484: col[0] = idx + 1;
485: col[1] = idx + 5;
486: val[0] = 1 / Tq0p[i];
487: val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
488: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
490: /* fgen[idx+2] = - w + w_s; */
491: row[0] = idx + 2;
492: col[0] = idx + 2;
493: col[1] = idx + 3;
494: val[0] = 0;
495: val[1] = -1;
496: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
498: /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
499: row[0] = idx + 3;
500: col[0] = idx;
501: col[1] = idx + 1;
502: col[2] = idx + 3;
503: col[3] = idx + 4;
504: col[4] = idx + 5;
505: val[0] = Iq / M[i];
506: val[1] = Id / M[i];
507: val[2] = D[i] / M[i];
508: val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
509: val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
510: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
512: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
513: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
514: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
516: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
518: Zdq_inv[0] = Rs[i] / det;
519: Zdq_inv[1] = Xqp[i] / det;
520: Zdq_inv[2] = -Xdp[i] / det;
521: Zdq_inv[3] = Rs[i] / det;
523: dVd_dVr = PetscSinScalar(delta);
524: dVd_dVi = -PetscCosScalar(delta);
525: dVq_dVr = PetscCosScalar(delta);
526: dVq_dVi = PetscSinScalar(delta);
527: dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
528: dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
530: /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
531: row[0] = idx + 4;
532: col[0] = idx;
533: col[1] = idx + 1;
534: col[2] = idx + 2;
535: val[0] = -Zdq_inv[1];
536: val[1] = -Zdq_inv[0];
537: val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
538: col[3] = idx + 4;
539: col[4] = net_start + 2 * gbus[i];
540: col[5] = net_start + 2 * gbus[i] + 1;
541: val[3] = 1;
542: val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
543: val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
544: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
546: /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
547: row[0] = idx + 5;
548: col[0] = idx;
549: col[1] = idx + 1;
550: col[2] = idx + 2;
551: val[0] = -Zdq_inv[3];
552: val[1] = -Zdq_inv[2];
553: val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
554: col[3] = idx + 5;
555: col[4] = net_start + 2 * gbus[i];
556: col[5] = net_start + 2 * gbus[i] + 1;
557: val[3] = 1;
558: val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
559: val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
560: PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
562: dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
563: dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
564: dIGr_dId = PetscSinScalar(delta);
565: dIGr_dIq = PetscCosScalar(delta);
566: dIGi_dId = -PetscCosScalar(delta);
567: dIGi_dIq = PetscSinScalar(delta);
569: /* fnet[2*gbus[i]] -= IGi; */
570: row[0] = net_start + 2 * gbus[i];
571: col[0] = idx + 2;
572: col[1] = idx + 4;
573: col[2] = idx + 5;
574: val[0] = -dIGi_ddelta;
575: val[1] = -dIGi_dId;
576: val[2] = -dIGi_dIq;
577: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
579: /* fnet[2*gbus[i]+1] -= IGr; */
580: row[0] = net_start + 2 * gbus[i] + 1;
581: col[0] = idx + 2;
582: col[1] = idx + 4;
583: col[2] = idx + 5;
584: val[0] = -dIGr_ddelta;
585: val[1] = -dIGr_dId;
586: val[2] = -dIGr_dIq;
587: PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
589: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
591: /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
592: /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
594: dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
596: row[0] = idx + 6;
597: col[0] = idx + 6;
598: col[1] = idx + 8;
599: val[0] = (KE[i] + dSE_dEfd) / TE[i];
600: val[1] = -1 / TE[i];
601: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
603: /* Exciter differential equations */
605: /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
606: row[0] = idx + 7;
607: col[0] = idx + 6;
608: col[1] = idx + 7;
609: val[0] = (-KF[i] / TF[i]) / TF[i];
610: val[1] = 1 / TF[i];
611: PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
613: /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
614: /* Vm = (Vd^2 + Vq^2)^0.5; */
615: dVm_dVd = Vd / Vm;
616: dVm_dVq = Vq / Vm;
617: dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
618: dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
619: row[0] = idx + 8;
620: col[0] = idx + 6;
621: col[1] = idx + 7;
622: col[2] = idx + 8;
623: val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
624: val[1] = -KA[i] / TA[i];
625: val[2] = 1 / TA[i];
626: col[3] = net_start + 2 * gbus[i];
627: col[4] = net_start + 2 * gbus[i] + 1;
628: val[3] = KA[i] * dVm_dVr / TA[i];
629: val[4] = KA[i] * dVm_dVi / TA[i];
630: PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
631: idx = idx + 9;
632: }
634: for (i = 0; i < nbus; i++) {
635: PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
636: row[0] = net_start + 2 * i;
637: for (PetscInt k = 0; k < ncols; k++) {
638: col[k] = net_start + cols[k];
639: val[k] = yvals[k];
640: }
641: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
642: PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
644: PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
645: row[0] = net_start + 2 * i + 1;
646: for (PetscInt k = 0; k < ncols; k++) {
647: col[k] = net_start + cols[k];
648: val[k] = yvals[k];
649: }
650: PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
651: PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
652: }
654: PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
655: PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));
657: PetscCall(VecGetArray(user->V0, &v0));
658: for (i = 0; i < nload; i++) {
659: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
660: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
661: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
662: Vm2 = Vm * Vm;
663: Vm4 = Vm2 * Vm2;
664: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
665: PD = QD = 0.0;
666: dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
667: for (PetscInt k = 0; k < ld_nsegsp[i]; k++) {
668: PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
669: dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vr * PetscPowScalar(Vm, ld_betap[k] - 2);
670: dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vi * PetscPowScalar(Vm, ld_betap[k] - 2);
671: }
672: for (PetscInt k = 0; k < ld_nsegsq[i]; k++) {
673: QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
674: dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vr * PetscPowScalar(Vm, ld_betaq[k] - 2);
675: dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vi * PetscPowScalar(Vm, ld_betaq[k] - 2);
676: }
678: /* IDr = (PD*Vr + QD*Vi)/Vm2; */
679: /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
681: dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
682: dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
684: dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
685: dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
687: /* fnet[2*lbus[i]] += IDi; */
688: row[0] = net_start + 2 * lbus[i];
689: col[0] = net_start + 2 * lbus[i];
690: col[1] = net_start + 2 * lbus[i] + 1;
691: val[0] = dIDi_dVr;
692: val[1] = dIDi_dVi;
693: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
694: /* fnet[2*lbus[i]+1] += IDr; */
695: row[0] = net_start + 2 * lbus[i] + 1;
696: col[0] = net_start + 2 * lbus[i];
697: col[1] = net_start + 2 * lbus[i] + 1;
698: val[0] = dIDr_dVr;
699: val[1] = dIDr_dVi;
700: PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
701: }
702: PetscCall(VecRestoreArray(user->V0, &v0));
704: PetscCall(VecRestoreArray(Xgen, &xgen));
705: PetscCall(VecRestoreArray(Xnet, &xnet));
707: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
709: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
710: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: /*
715: J = [I, 0
716: dg_dx, dg_dy]
717: */
718: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
719: {
720: Userctx *user = (Userctx *)ctx;
722: PetscFunctionBegin;
723: PetscCall(ResidualJacobian(snes, X, A, B, ctx));
724: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
725: PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: /*
730: J = [a*I-df_dx, -df_dy
731: dg_dx, dg_dy]
732: */
734: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
735: {
736: SNES snes;
737: PetscScalar atmp = (PetscScalar)a;
738: PetscInt row;
740: PetscFunctionBegin;
741: user->t = t;
743: PetscCall(TSGetSNES(ts, &snes));
744: PetscCall(ResidualJacobian(snes, X, A, B, user));
745: for (PetscInt i = 0; i < ngen; i++) {
746: row = 9 * i;
747: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
748: row = 9 * i + 1;
749: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
750: row = 9 * i + 2;
751: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
752: row = 9 * i + 3;
753: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
754: row = 9 * i + 6;
755: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
756: row = 9 * i + 7;
757: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
758: row = 9 * i + 8;
759: PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
760: }
761: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
762: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, Userctx *user)
767: {
768: PetscScalar *r;
769: const PetscScalar *u;
770: PetscInt idx;
771: Vec Xgen, Xnet;
772: PetscScalar *xgen;
774: PetscFunctionBegin;
775: PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
776: PetscCall(DMCompositeScatter(user->dmpgrid, U, Xgen, Xnet));
778: PetscCall(VecGetArray(Xgen, &xgen));
780: PetscCall(VecGetArrayRead(U, &u));
781: PetscCall(VecGetArray(R, &r));
782: r[0] = 0.;
784: idx = 0;
785: for (PetscInt i = 0; i < ngen; i++) {
786: 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);
787: idx += 9;
788: }
789: PetscCall(VecRestoreArray(R, &r));
790: PetscCall(VecRestoreArrayRead(U, &u));
791: PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
792: PetscFunctionReturn(PETSC_SUCCESS);
793: }
795: static PetscErrorCode MonitorUpdateQ(TS ts, PetscInt stepnum, PetscReal time, Vec X, PetscCtx ctx0)
796: {
797: Vec C, *Y;
798: PetscInt Nr;
799: PetscReal h, theta;
800: Userctx *ctx = (Userctx *)ctx0;
802: PetscFunctionBegin;
803: theta = 0.5;
804: PetscCall(TSGetStages(ts, &Nr, &Y));
805: PetscCall(TSGetTimeStep(ts, &h));
806: PetscCall(VecDuplicate(ctx->vec_q, &C));
807: /* compute integrals */
808: if (stepnum > 0) {
809: PetscCall(CostIntegrand(ts, time, X, C, ctx));
810: PetscCall(VecAXPY(ctx->vec_q, h * theta, C));
811: PetscCall(CostIntegrand(ts, time + h * theta, Y[0], C, ctx));
812: PetscCall(VecAXPY(ctx->vec_q, h * (1 - theta), C));
813: }
814: PetscCall(VecDestroy(&C));
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: int main(int argc, char **argv)
819: {
820: Userctx user;
821: Vec p;
822: PetscScalar *x_ptr;
823: PetscMPIInt size;
824: KSP ksp;
825: PC pc;
826: PetscInt *idx2;
827: Tao tao;
828: Vec lowerb, upperb;
830: PetscFunctionBeginUser;
831: PetscFunctionBeginUser;
832: PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
833: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
834: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
836: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 1, &user.vec_q));
838: user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
839: user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
840: user.neqs_pgrid = user.neqs_gen + user.neqs_net;
842: /* Create indices for differential and algebraic equations */
843: PetscCall(PetscMalloc1(7 * ngen, &idx2));
844: for (PetscInt i = 0; i < ngen; i++) {
845: idx2[7 * i] = 9 * i;
846: idx2[7 * i + 1] = 9 * i + 1;
847: idx2[7 * i + 2] = 9 * i + 2;
848: idx2[7 * i + 3] = 9 * i + 3;
849: idx2[7 * i + 4] = 9 * i + 6;
850: idx2[7 * i + 5] = 9 * i + 7;
851: idx2[7 * i + 6] = 9 * i + 8;
852: }
853: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
854: PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
855: PetscCall(PetscFree(idx2));
857: /* Set run time options */
858: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
859: {
860: user.tfaulton = 1.0;
861: user.tfaultoff = 1.2;
862: user.Rfault = 0.0001;
863: user.faultbus = 8;
864: PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
865: PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
866: PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
867: user.t0 = 0.0;
868: user.tmax = 1.5;
869: PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
870: PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
871: user.freq_u = 61.0;
872: user.freq_l = 59.0;
873: user.pow = 2;
874: PetscCall(PetscOptionsReal("-frequ", "", "", user.freq_u, &user.freq_u, NULL));
875: PetscCall(PetscOptionsReal("-freql", "", "", user.freq_l, &user.freq_l, NULL));
876: PetscCall(PetscOptionsInt("-pow", "", "", user.pow, &user.pow, NULL));
877: }
878: PetscOptionsEnd();
880: /* Create DMs for generator and network subsystems */
881: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
882: PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
883: PetscCall(DMSetFromOptions(user.dmgen));
884: PetscCall(DMSetUp(user.dmgen));
885: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
886: PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
887: PetscCall(DMSetFromOptions(user.dmnet));
888: PetscCall(DMSetUp(user.dmnet));
889: /* Create a composite DM packer and add the two DMs */
890: PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
891: PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
892: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
893: PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));
895: /* Create TAO solver and set desired solution method */
896: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
897: PetscCall(TaoSetType(tao, TAOBLMVM));
898: /*
899: Optimization starts
900: */
901: /* Set initial solution guess */
902: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 3, &p));
903: PetscCall(VecGetArray(p, &x_ptr));
904: x_ptr[0] = PG[0];
905: x_ptr[1] = PG[1];
906: x_ptr[2] = PG[2];
907: PetscCall(VecRestoreArray(p, &x_ptr));
909: PetscCall(TaoSetSolution(tao, p));
910: /* Set routine for function and gradient evaluation */
911: PetscCall(TaoSetObjective(tao, FormFunction, (void *)&user));
912: PetscCall(TaoSetGradient(tao, NULL, TaoDefaultComputeGradient, (void *)&user));
914: /* Set bounds for the optimization */
915: PetscCall(VecDuplicate(p, &lowerb));
916: PetscCall(VecDuplicate(p, &upperb));
917: PetscCall(VecGetArray(lowerb, &x_ptr));
918: x_ptr[0] = 0.5;
919: x_ptr[1] = 0.5;
920: x_ptr[2] = 0.5;
921: PetscCall(VecRestoreArray(lowerb, &x_ptr));
922: PetscCall(VecGetArray(upperb, &x_ptr));
923: x_ptr[0] = 2.0;
924: x_ptr[1] = 2.0;
925: x_ptr[2] = 2.0;
926: PetscCall(VecRestoreArray(upperb, &x_ptr));
927: PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));
929: /* Check for any TAO command line options */
930: PetscCall(TaoSetFromOptions(tao));
931: PetscCall(TaoGetKSP(tao, &ksp));
932: if (ksp) {
933: PetscCall(KSPGetPC(ksp, &pc));
934: PetscCall(PCSetType(pc, PCNONE));
935: }
937: /* SOLVE THE APPLICATION */
938: PetscCall(TaoSolve(tao));
940: PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
941: /* Free TAO data structures */
942: PetscCall(TaoDestroy(&tao));
943: PetscCall(VecDestroy(&user.vec_q));
944: PetscCall(VecDestroy(&lowerb));
945: PetscCall(VecDestroy(&upperb));
946: PetscCall(VecDestroy(&p));
947: PetscCall(DMDestroy(&user.dmgen));
948: PetscCall(DMDestroy(&user.dmnet));
949: PetscCall(DMDestroy(&user.dmpgrid));
950: PetscCall(ISDestroy(&user.is_diff));
951: PetscCall(ISDestroy(&user.is_alg));
952: PetscCall(PetscFinalize());
953: return 0;
954: }
956: /* ------------------------------------------------------------------ */
957: /*
958: FormFunction - Evaluates the function and corresponding gradient.
960: Input Parameters:
961: tao - the Tao context
962: X - the input vector
963: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
965: Output Parameters:
966: f - the newly evaluated function
967: */
968: PetscErrorCode FormFunction(Tao tao, Vec P, PetscReal *f, PetscCtx ctx0)
969: {
970: TS ts;
971: SNES snes_alg;
972: Userctx *ctx = (Userctx *)ctx0;
973: Vec X;
974: Mat J;
975: /* sensitivity context */
976: PetscScalar *x_ptr;
977: PetscViewer Xview, Ybusview;
978: Vec F_alg;
979: Vec Xdot;
980: PetscInt row_loc, col_loc;
981: PetscScalar val;
983: PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
984: PG[0] = x_ptr[0];
985: PG[1] = x_ptr[1];
986: PG[2] = x_ptr[2];
987: PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
989: ctx->stepnum = 0;
991: PetscCall(VecZeroEntries(ctx->vec_q));
993: /* Read initial voltage vector and Ybus */
994: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
995: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));
997: PetscCall(VecCreate(PETSC_COMM_WORLD, &ctx->V0));
998: PetscCall(VecSetSizes(ctx->V0, PETSC_DECIDE, ctx->neqs_net));
999: PetscCall(VecLoad(ctx->V0, Xview));
1001: PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx->Ybus));
1002: PetscCall(MatSetSizes(ctx->Ybus, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_net, ctx->neqs_net));
1003: PetscCall(MatSetType(ctx->Ybus, MATBAIJ));
1004: /* PetscCall(MatSetBlockSize(ctx->Ybus,2)); */
1005: PetscCall(MatLoad(ctx->Ybus, Ybusview));
1007: PetscCall(PetscViewerDestroy(&Xview));
1008: PetscCall(PetscViewerDestroy(&Ybusview));
1010: PetscCall(DMCreateGlobalVector(ctx->dmpgrid, &X));
1012: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
1013: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx->neqs_pgrid, ctx->neqs_pgrid));
1014: PetscCall(MatSetFromOptions(J));
1015: PetscCall(PreallocateJacobian(J, ctx));
1017: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1018: Create timestepping solver context
1019: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1020: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1021: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1022: PetscCall(TSSetType(ts, TSCN));
1023: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, ctx));
1024: PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobianFn *)IJacobian, ctx));
1025: PetscCall(TSSetApplicationContext(ts, ctx));
1027: PetscCall(TSMonitorSet(ts, MonitorUpdateQ, ctx, NULL));
1028: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1029: Set initial conditions
1030: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1031: PetscCall(SetInitialGuess(X, ctx));
1033: PetscCall(VecDuplicate(X, &F_alg));
1034: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1035: PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, ctx));
1036: PetscCall(MatZeroEntries(J));
1037: PetscCall(SNESSetJacobian(snes_alg, J, J, AlgJacobian, ctx));
1038: PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1039: PetscCall(SNESSetFromOptions(snes_alg));
1040: ctx->alg_flg = PETSC_TRUE;
1041: /* Solve the algebraic equations */
1042: PetscCall(SNESSolve(snes_alg, NULL, X));
1044: /* Just to set up the Jacobian structure */
1045: PetscCall(VecDuplicate(X, &Xdot));
1046: PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, J, J, ctx));
1047: PetscCall(VecDestroy(&Xdot));
1049: ctx->stepnum++;
1051: PetscCall(TSSetTimeStep(ts, 0.01));
1052: PetscCall(TSSetMaxTime(ts, ctx->tmax));
1053: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1054: PetscCall(TSSetFromOptions(ts));
1056: /* Prefault period */
1057: ctx->alg_flg = PETSC_FALSE;
1058: PetscCall(TSSetTime(ts, 0.0));
1059: PetscCall(TSSetMaxTime(ts, ctx->tfaulton));
1060: PetscCall(TSSolve(ts, X));
1062: /* Create the nonlinear solver for solving the algebraic system */
1063: /* Note that although the algebraic system needs to be solved only for
1064: Idq and V, we reuse the entire system including xgen. The xgen
1065: variables are held constant by setting their residuals to 0 and
1066: putting a 1 on the Jacobian diagonal for xgen rows
1067: */
1068: PetscCall(MatZeroEntries(J));
1070: /* Apply disturbance - resistive fault at ctx->faultbus */
1071: /* This is done by adding shunt conductance to the diagonal location
1072: in the Ybus matrix */
1073: row_loc = 2 * ctx->faultbus;
1074: col_loc = 2 * ctx->faultbus + 1; /* Location for G */
1075: val = 1 / ctx->Rfault;
1076: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1077: row_loc = 2 * ctx->faultbus + 1;
1078: col_loc = 2 * ctx->faultbus; /* Location for G */
1079: val = 1 / ctx->Rfault;
1080: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1082: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1083: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1085: ctx->alg_flg = PETSC_TRUE;
1086: /* Solve the algebraic equations */
1087: PetscCall(SNESSolve(snes_alg, NULL, X));
1089: ctx->stepnum++;
1091: /* Disturbance period */
1092: ctx->alg_flg = PETSC_FALSE;
1093: PetscCall(TSSetTime(ts, ctx->tfaulton));
1094: PetscCall(TSSetMaxTime(ts, ctx->tfaultoff));
1095: PetscCall(TSSolve(ts, X));
1097: /* Remove the fault */
1098: row_loc = 2 * ctx->faultbus;
1099: col_loc = 2 * ctx->faultbus + 1;
1100: val = -1 / ctx->Rfault;
1101: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1102: row_loc = 2 * ctx->faultbus + 1;
1103: col_loc = 2 * ctx->faultbus;
1104: val = -1 / ctx->Rfault;
1105: PetscCall(MatSetValues(ctx->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
1107: PetscCall(MatAssemblyBegin(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1108: PetscCall(MatAssemblyEnd(ctx->Ybus, MAT_FINAL_ASSEMBLY));
1110: PetscCall(MatZeroEntries(J));
1112: ctx->alg_flg = PETSC_TRUE;
1114: /* Solve the algebraic equations */
1115: PetscCall(SNESSolve(snes_alg, NULL, X));
1117: ctx->stepnum++;
1119: /* Post-disturbance period */
1120: ctx->alg_flg = PETSC_TRUE;
1121: PetscCall(TSSetTime(ts, ctx->tfaultoff));
1122: PetscCall(TSSetMaxTime(ts, ctx->tmax));
1123: PetscCall(TSSolve(ts, X));
1125: PetscCall(VecGetArray(ctx->vec_q, &x_ptr));
1126: *f = x_ptr[0];
1127: PetscCall(VecRestoreArray(ctx->vec_q, &x_ptr));
1129: PetscCall(MatDestroy(&ctx->Ybus));
1130: PetscCall(VecDestroy(&ctx->V0));
1131: PetscCall(SNESDestroy(&snes_alg));
1132: PetscCall(VecDestroy(&F_alg));
1133: PetscCall(MatDestroy(&J));
1134: PetscCall(VecDestroy(&X));
1135: PetscCall(TSDestroy(&ts));
1137: return 0;
1138: }
1140: /*TEST
1142: build:
1143: requires: double !complex !defined(USE_64BIT_INDICES)
1145: test:
1146: args: -viewer_binary_skip_info -tao_monitor -tao_gttol .2
1147: localrunfiles: petscoptions X.bin Ybus.bin
1149: TEST*/