Actual source code: ex9opt.c
2: static char help[] = "Basic equation for generator stability analysis.\n";
\begin{eqnarray}
\frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
\frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
\end{eqnarray}
Ensemble of initial conditions
./ex2 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
Fault at .1 seconds
./ex2 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
Initial conditions same as when fault is ended
./ex2 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
22: /*
23: Include "petscts.h" so that we can use TS solvers. Note that this
24: file automatically includes:
25: petscsys.h - base PETSc routines petscvec.h - vectors
26: petscmat.h - matrices
27: petscis.h - index sets petscksp.h - Krylov subspace methods
28: petscviewer.h - viewers petscpc.h - preconditioners
29: petscksp.h - linear solvers
30: */
32: #include <petsctao.h>
33: #include <petscts.h>
35: typedef struct {
36: TS ts;
37: PetscScalar H,D,omega_b,omega_s,Pmax,Pm,E,V,X,u_s,c;
38: PetscInt beta;
39: PetscReal tf,tcl,dt;
40: } AppCtx;
42: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
43: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
45: /*
46: Defines the ODE passed to the ODE solver
47: */
48: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
49: {
50: PetscScalar *f,Pmax;
51: const PetscScalar *u;
53: /* The next three lines allow us to access the entries of the vectors directly */
54: VecGetArrayRead(U,&u);
55: VecGetArray(F,&f);
56: if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
57: else Pmax = ctx->Pmax;
59: f[0] = ctx->omega_b*(u[1] - ctx->omega_s);
60: f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H);
62: VecRestoreArrayRead(U,&u);
63: VecRestoreArray(F,&f);
64: return 0;
65: }
67: /*
68: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
69: */
70: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx)
71: {
72: PetscInt rowcol[] = {0,1};
73: PetscScalar J[2][2],Pmax;
74: const PetscScalar *u;
76: VecGetArrayRead(U,&u);
77: if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
78: else Pmax = ctx->Pmax;
80: J[0][0] = 0; J[0][1] = ctx->omega_b;
81: J[1][1] = -ctx->D*ctx->omega_s/(2.0*ctx->H); J[1][0] = -Pmax*PetscCosScalar(u[0])*ctx->omega_s/(2.0*ctx->H);
83: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
84: VecRestoreArrayRead(U,&u);
86: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
87: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
88: if (A != B) {
89: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
91: }
92: return 0;
93: }
95: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
96: {
97: PetscInt row[] = {0,1},col[]={0};
98: PetscScalar J[2][1];
99: AppCtx *ctx=(AppCtx*)ctx0;
102: J[0][0] = 0;
103: J[1][0] = ctx->omega_s/(2.0*ctx->H);
104: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
105: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
106: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
107: return 0;
108: }
110: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
111: {
112: PetscScalar *r;
113: const PetscScalar *u;
115: VecGetArrayRead(U,&u);
116: VecGetArray(R,&r);
117: r[0] = ctx->c*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta);
118: VecRestoreArray(R,&r);
119: VecRestoreArrayRead(U,&u);
120: return 0;
121: }
123: static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx *ctx)
124: {
125: PetscScalar ru[1];
126: const PetscScalar *u;
127: PetscInt row[] = {0},col[] = {0};
129: VecGetArrayRead(U,&u);
130: ru[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta-1);
131: VecRestoreArrayRead(U,&u);
132: MatSetValues(DRDU,1,row,1,col,ru,INSERT_VALUES);
133: MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);
134: MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);
135: return 0;
136: }
138: static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx *ctx)
139: {
140: MatZeroEntries(DRDP);
141: MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);
142: MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);
143: return 0;
144: }
146: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
147: {
148: PetscScalar *y,sensip;
149: const PetscScalar *x;
151: VecGetArrayRead(lambda,&x);
152: VecGetArray(mu,&y);
153: sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
154: y[0] = sensip;
155: VecRestoreArray(mu,&y);
156: VecRestoreArrayRead(lambda,&x);
157: return 0;
158: }
160: int main(int argc,char **argv)
161: {
162: Vec p;
163: PetscScalar *x_ptr;
165: PetscMPIInt size;
166: AppCtx ctx;
167: Vec lowerb,upperb;
168: Tao tao;
169: KSP ksp;
170: PC pc;
171: Vec U,lambda[1],mu[1];
172: Mat A; /* Jacobian matrix */
173: Mat Jacp; /* Jacobian matrix */
174: Mat DRDU,DRDP;
175: PetscInt n = 2;
176: TS quadts;
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Initialize program
180: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: PetscInitialize(&argc,&argv,NULL,help);
183: MPI_Comm_size(PETSC_COMM_WORLD,&size);
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Set runtime options
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
190: {
191: ctx.beta = 2;
192: ctx.c = PetscRealConstant(10000.0);
193: ctx.u_s = PetscRealConstant(1.0);
194: ctx.omega_s = PetscRealConstant(1.0);
195: ctx.omega_b = PetscRealConstant(120.0)*PETSC_PI;
196: ctx.H = PetscRealConstant(5.0);
197: PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
198: ctx.D = PetscRealConstant(5.0);
199: PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
200: ctx.E = PetscRealConstant(1.1378);
201: ctx.V = PetscRealConstant(1.0);
202: ctx.X = PetscRealConstant(0.545);
203: ctx.Pmax = ctx.E*ctx.V/ctx.X;
204: PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
205: ctx.Pm = PetscRealConstant(1.0194);
206: PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
207: ctx.tf = PetscRealConstant(0.1);
208: ctx.tcl = PetscRealConstant(0.2);
209: PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
210: PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
212: }
213: PetscOptionsEnd();
215: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: Create necessary matrix and vectors
217: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218: MatCreate(PETSC_COMM_WORLD,&A);
219: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
220: MatSetType(A,MATDENSE);
221: MatSetFromOptions(A);
222: MatSetUp(A);
224: MatCreateVecs(A,&U,NULL);
226: MatCreate(PETSC_COMM_WORLD,&Jacp);
227: MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
228: MatSetFromOptions(Jacp);
229: MatSetUp(Jacp);
230: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&DRDP);
231: MatSetUp(DRDP);
232: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,2,NULL,&DRDU);
233: MatSetUp(DRDU);
235: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: Create timestepping solver context
237: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
238: TSCreate(PETSC_COMM_WORLD,&ctx.ts);
239: TSSetProblemType(ctx.ts,TS_NONLINEAR);
240: TSSetEquationType(ctx.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
241: TSSetType(ctx.ts,TSRK);
242: TSSetRHSFunction(ctx.ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
243: TSSetRHSJacobian(ctx.ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);
244: TSSetExactFinalTime(ctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
246: MatCreateVecs(A,&lambda[0],NULL);
247: MatCreateVecs(Jacp,&mu[0],NULL);
248: TSSetCostGradients(ctx.ts,1,lambda,mu);
249: TSSetRHSJacobianP(ctx.ts,Jacp,RHSJacobianP,&ctx);
251: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
252: Set solver options
253: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254: TSSetMaxTime(ctx.ts,PetscRealConstant(1.0));
255: TSSetTimeStep(ctx.ts,PetscRealConstant(0.01));
256: TSSetFromOptions(ctx.ts);
258: TSGetTimeStep(ctx.ts,&ctx.dt); /* save the stepsize */
260: TSCreateQuadratureTS(ctx.ts,PETSC_TRUE,&quadts);
261: TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);
262: TSSetRHSJacobian(quadts,DRDU,DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);
263: TSSetRHSJacobianP(quadts,DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&ctx);
264: TSSetSolution(ctx.ts,U);
266: /* Create TAO solver and set desired solution method */
267: TaoCreate(PETSC_COMM_WORLD,&tao);
268: TaoSetType(tao,TAOBLMVM);
270: /*
271: Optimization starts
272: */
273: /* Set initial solution guess */
274: VecCreateSeq(PETSC_COMM_WORLD,1,&p);
275: VecGetArray(p,&x_ptr);
276: x_ptr[0] = ctx.Pm;
277: VecRestoreArray(p,&x_ptr);
279: TaoSetSolution(tao,p);
280: /* Set routine for function and gradient evaluation */
281: TaoSetObjective(tao,FormFunction,(void *)&ctx);
282: TaoSetGradient(tao,NULL,FormGradient,(void *)&ctx);
284: /* Set bounds for the optimization */
285: VecDuplicate(p,&lowerb);
286: VecDuplicate(p,&upperb);
287: VecGetArray(lowerb,&x_ptr);
288: x_ptr[0] = 0.;
289: VecRestoreArray(lowerb,&x_ptr);
290: VecGetArray(upperb,&x_ptr);
291: x_ptr[0] = PetscRealConstant(1.1);
292: VecRestoreArray(upperb,&x_ptr);
293: TaoSetVariableBounds(tao,lowerb,upperb);
295: /* Check for any TAO command line options */
296: TaoSetFromOptions(tao);
297: TaoGetKSP(tao,&ksp);
298: if (ksp) {
299: KSPGetPC(ksp,&pc);
300: PCSetType(pc,PCNONE);
301: }
303: /* SOLVE THE APPLICATION */
304: TaoSolve(tao);
306: VecView(p,PETSC_VIEWER_STDOUT_WORLD);
307: /* Free TAO data structures */
308: TaoDestroy(&tao);
309: VecDestroy(&p);
310: VecDestroy(&lowerb);
311: VecDestroy(&upperb);
313: TSDestroy(&ctx.ts);
314: VecDestroy(&U);
315: MatDestroy(&A);
316: MatDestroy(&Jacp);
317: MatDestroy(&DRDU);
318: MatDestroy(&DRDP);
319: VecDestroy(&lambda[0]);
320: VecDestroy(&mu[0]);
321: PetscFinalize();
322: return 0;
323: }
325: /* ------------------------------------------------------------------ */
326: /*
327: FormFunction - Evaluates the function
329: Input Parameters:
330: tao - the Tao context
331: X - the input vector
332: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
334: Output Parameters:
335: f - the newly evaluated function
336: */
337: PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0)
338: {
339: AppCtx *ctx = (AppCtx*)ctx0;
340: TS ts = ctx->ts;
341: Vec U; /* solution will be stored here */
342: PetscScalar *u;
343: PetscScalar *x_ptr;
344: Vec q;
346: VecGetArrayRead(P,(const PetscScalar**)&x_ptr);
347: ctx->Pm = x_ptr[0];
348: VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);
350: /* reset time */
351: TSSetTime(ts,0.0);
352: /* reset step counter, this is critical for adjoint solver */
353: TSSetStepNumber(ts,0);
354: /* reset step size, the step size becomes negative after TSAdjointSolve */
355: TSSetTimeStep(ts,ctx->dt);
356: /* reinitialize the integral value */
357: TSGetCostIntegral(ts,&q);
358: VecSet(q,0.0);
360: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
361: Set initial conditions
362: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
363: TSGetSolution(ts,&U);
364: VecGetArray(U,&u);
365: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
366: u[1] = PetscRealConstant(1.0);
367: VecRestoreArray(U,&u);
369: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
370: Solve nonlinear system
371: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
372: TSSolve(ts,U);
373: TSGetCostIntegral(ts,&q);
374: VecGetArray(q,&x_ptr);
375: *f = -ctx->Pm + x_ptr[0];
376: VecRestoreArray(q,&x_ptr);
377: return 0;
378: }
380: PetscErrorCode FormGradient(Tao tao,Vec P,Vec G,void *ctx0)
381: {
382: AppCtx *ctx = (AppCtx*)ctx0;
383: TS ts = ctx->ts;
384: Vec U; /* solution will be stored here */
385: PetscReal ftime;
386: PetscInt steps;
387: PetscScalar *u;
388: PetscScalar *x_ptr,*y_ptr;
389: Vec *lambda,q,*mu;
391: VecGetArrayRead(P,(const PetscScalar**)&x_ptr);
392: ctx->Pm = x_ptr[0];
393: VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);
395: /* reset time */
396: TSSetTime(ts,0.0);
397: /* reset step counter, this is critical for adjoint solver */
398: TSSetStepNumber(ts,0);
399: /* reset step size, the step size becomes negative after TSAdjointSolve */
400: TSSetTimeStep(ts,ctx->dt);
401: /* reinitialize the integral value */
402: TSGetCostIntegral(ts,&q);
403: VecSet(q,0.0);
405: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
406: Set initial conditions
407: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
408: TSGetSolution(ts,&U);
409: VecGetArray(U,&u);
410: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
411: u[1] = PetscRealConstant(1.0);
412: VecRestoreArray(U,&u);
414: /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
415: TSSetSaveTrajectory(ts);
416: TSSetFromOptions(ts);
418: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
419: Solve nonlinear system
420: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
421: TSSolve(ts,U);
423: TSGetSolveTime(ts,&ftime);
424: TSGetStepNumber(ts,&steps);
426: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
427: Adjoint model starts here
428: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
429: TSGetCostGradients(ts,NULL,&lambda,&mu);
430: /* Set initial conditions for the adjoint integration */
431: VecGetArray(lambda[0],&y_ptr);
432: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
433: VecRestoreArray(lambda[0],&y_ptr);
434: VecGetArray(mu[0],&x_ptr);
435: x_ptr[0] = PetscRealConstant(-1.0);
436: VecRestoreArray(mu[0],&x_ptr);
438: TSAdjointSolve(ts);
439: TSGetCostIntegral(ts,&q);
440: ComputeSensiP(lambda[0],mu[0],ctx);
441: VecCopy(mu[0],G);
442: return 0;
443: }
445: /*TEST
447: build:
448: requires: !complex
450: test:
451: args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason
453: test:
454: suffix: 2
455: args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason -tao_test_gradient
457: TEST*/