Actual source code: ex3opt.c
2: static char help[] = "Finds optimal parameter P_m for the generator system while maintaining generator stability.\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}
13: /*
14: This code demonstrates how to solve a ODE-constrained optimization problem with TAO, TSEvent, TSAdjoint and TS.
15: The problem features discontinuities and a cost function in integral form.
16: The gradient is computed with the discrete adjoint of an implicit theta method, see ex3adj.c for details.
17: */
19: #include <petsctao.h>
20: #include <petscts.h>
21: #include "ex3.h"
23: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
25: PetscErrorCode monitor(Tao tao,AppCtx *ctx)
26: {
27: FILE *fp;
28: PetscInt iterate;
29: PetscReal f,gnorm,cnorm,xdiff;
30: TaoConvergedReason reason;
33: TaoGetSolutionStatus(tao,&iterate,&f,&gnorm,&cnorm,&xdiff,&reason);
35: fp = fopen("ex3opt_conv.out","a");
36: PetscFPrintf(PETSC_COMM_WORLD,fp,"%D %g\n",iterate,(double)gnorm);
37: fclose(fp);
38: return 0;
39: }
41: int main(int argc,char **argv)
42: {
43: Vec p;
44: PetscScalar *x_ptr;
45: PetscErrorCode ierr;
46: PetscMPIInt size;
47: AppCtx ctx;
48: Tao tao;
49: KSP ksp;
50: PC pc;
51: Vec lambda[1],mu[1],lowerb,upperb;
52: PetscBool printtofile;
53: PetscInt direction[2];
54: PetscBool terminate[2];
55: Mat qgrad; /* Forward sesivitiy */
56: Mat sp; /* Forward sensitivity matrix */
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Initialize program
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: PetscInitialize(&argc,&argv,NULL,help);
63: MPI_Comm_size(PETSC_COMM_WORLD,&size);
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Set runtime options
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
70: {
71: ctx.beta = 2;
72: ctx.c = 10000.0;
73: ctx.u_s = 1.0;
74: ctx.omega_s = 1.0;
75: ctx.omega_b = 120.0*PETSC_PI;
76: ctx.H = 5.0;
77: PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
78: ctx.D = 5.0;
79: PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
80: ctx.E = 1.1378;
81: ctx.V = 1.0;
82: ctx.X = 0.545;
83: ctx.Pmax = ctx.E*ctx.V/ctx.X;
84: ctx.Pmax_ini = ctx.Pmax;
85: PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
86: ctx.Pm = 1.06;
87: PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
88: ctx.tf = 0.1;
89: ctx.tcl = 0.2;
90: PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
91: PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
92: printtofile = PETSC_FALSE;
93: PetscOptionsBool("-printtofile","Print convergence results to file","",printtofile,&printtofile,NULL);
94: ctx.sa = SA_ADJ;
95: PetscOptionsEnum("-sa_method","Sensitivity analysis method (adj or tlm)","",SAMethods,(PetscEnum)ctx.sa,(PetscEnum*)&ctx.sa,NULL);
96: }
97: PetscOptionsEnd();
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create necessary matrix and vectors
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: MatCreate(PETSC_COMM_WORLD,&ctx.Jac);
103: MatSetSizes(ctx.Jac,2,2,PETSC_DETERMINE,PETSC_DETERMINE);
104: MatSetType(ctx.Jac,MATDENSE);
105: MatSetFromOptions(ctx.Jac);
106: MatSetUp(ctx.Jac);
107: MatCreate(PETSC_COMM_WORLD,&ctx.Jacp);
108: MatSetSizes(ctx.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
109: MatSetFromOptions(ctx.Jacp);
110: MatSetUp(ctx.Jacp);
111: MatCreateVecs(ctx.Jac,&ctx.U,NULL);
112: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&ctx.DRDP);
113: MatSetUp(ctx.DRDP);
114: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&ctx.DRDU);
115: MatSetUp(ctx.DRDU);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Create timestepping solver context
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: TSCreate(PETSC_COMM_WORLD,&ctx.ts);
121: TSSetProblemType(ctx.ts,TS_NONLINEAR);
122: TSSetType(ctx.ts,TSCN);
123: TSSetRHSFunction(ctx.ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
124: TSSetRHSJacobian(ctx.ts,ctx.Jac,ctx.Jac,(TSRHSJacobian)RHSJacobian,&ctx);
125: TSSetRHSJacobianP(ctx.ts,ctx.Jacp,RHSJacobianP,&ctx);
127: if (ctx.sa == SA_ADJ) {
128: MatCreateVecs(ctx.Jac,&lambda[0],NULL);
129: MatCreateVecs(ctx.Jacp,&mu[0],NULL);
130: TSSetSaveTrajectory(ctx.ts);
131: TSSetCostGradients(ctx.ts,1,lambda,mu);
132: TSCreateQuadratureTS(ctx.ts,PETSC_FALSE,&ctx.quadts);
133: TSSetRHSFunction(ctx.quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);
134: TSSetRHSJacobian(ctx.quadts,ctx.DRDU,ctx.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);
135: TSSetRHSJacobianP(ctx.quadts,ctx.DRDP,DRDPJacobianTranspose,&ctx);
136: }
137: if (ctx.sa == SA_TLM) {
138: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&qgrad);
139: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&sp);
140: TSForwardSetSensitivities(ctx.ts,1,sp);
141: TSCreateQuadratureTS(ctx.ts,PETSC_TRUE,&ctx.quadts);
142: TSForwardSetSensitivities(ctx.quadts,1,qgrad);
143: TSSetRHSFunction(ctx.quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);
144: TSSetRHSJacobian(ctx.quadts,ctx.DRDU,ctx.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);
145: TSSetRHSJacobianP(ctx.quadts,ctx.DRDP,DRDPJacobianTranspose,&ctx);
146: }
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Set solver options
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: TSSetMaxTime(ctx.ts,1.0);
152: TSSetExactFinalTime(ctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
153: TSSetTimeStep(ctx.ts,0.03125);
154: TSSetFromOptions(ctx.ts);
156: direction[0] = direction[1] = 1;
157: terminate[0] = terminate[1] = PETSC_FALSE;
158: TSSetEventHandler(ctx.ts,2,direction,terminate,EventFunction,PostEventFunction,&ctx);
160: /* Create TAO solver and set desired solution method */
161: TaoCreate(PETSC_COMM_WORLD,&tao);
162: TaoSetType(tao,TAOBLMVM);
163: if (printtofile) {
164: TaoSetMonitor(tao,(PetscErrorCode (*)(Tao, void*))monitor,(void *)&ctx,PETSC_NULL);
165: }
166: /*
167: Optimization starts
168: */
169: /* Set initial solution guess */
170: VecCreateSeq(PETSC_COMM_WORLD,1,&p);
171: VecGetArray(p,&x_ptr);
172: x_ptr[0] = ctx.Pm;
173: VecRestoreArray(p,&x_ptr);
175: TaoSetSolution(tao,p);
176: /* Set routine for function and gradient evaluation */
177: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&ctx);
179: /* Set bounds for the optimization */
180: VecDuplicate(p,&lowerb);
181: VecDuplicate(p,&upperb);
182: VecGetArray(lowerb,&x_ptr);
183: x_ptr[0] = 0.;
184: VecRestoreArray(lowerb,&x_ptr);
185: VecGetArray(upperb,&x_ptr);
186: x_ptr[0] = 1.1;
187: VecRestoreArray(upperb,&x_ptr);
188: TaoSetVariableBounds(tao,lowerb,upperb);
190: /* Check for any TAO command line options */
191: TaoSetFromOptions(tao);
192: TaoGetKSP(tao,&ksp);
193: if (ksp) {
194: KSPGetPC(ksp,&pc);
195: PCSetType(pc,PCNONE);
196: }
198: /* SOLVE THE APPLICATION */
199: TaoSolve(tao);
201: VecView(p,PETSC_VIEWER_STDOUT_WORLD);
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Free work space. All PETSc objects should be destroyed when they are no longer needed.
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: MatDestroy(&ctx.Jac);
207: MatDestroy(&ctx.Jacp);
208: MatDestroy(&ctx.DRDU);
209: MatDestroy(&ctx.DRDP);
210: VecDestroy(&ctx.U);
211: if (ctx.sa == SA_ADJ) {
212: VecDestroy(&lambda[0]);
213: VecDestroy(&mu[0]);
214: }
215: if (ctx.sa == SA_TLM) {
216: MatDestroy(&qgrad);
217: MatDestroy(&sp);
218: }
219: TSDestroy(&ctx.ts);
220: VecDestroy(&p);
221: VecDestroy(&lowerb);
222: VecDestroy(&upperb);
223: TaoDestroy(&tao);
224: PetscFinalize();
225: return 0;
226: }
228: /* ------------------------------------------------------------------ */
229: /*
230: FormFunctionGradient - Evaluates the function and corresponding gradient.
232: Input Parameters:
233: tao - the Tao context
234: X - the input vector
235: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
237: Output Parameters:
238: f - the newly evaluated function
239: G - the newly evaluated gradient
240: */
241: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx0)
242: {
243: AppCtx *ctx = (AppCtx*)ctx0;
244: PetscInt nadj;
245: PetscReal ftime;
246: PetscInt steps;
247: PetscScalar *u;
248: PetscScalar *x_ptr,*y_ptr;
249: Vec q;
250: Mat qgrad;
252: VecGetArrayRead(P,(const PetscScalar**)&x_ptr);
253: ctx->Pm = x_ptr[0];
254: VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);
256: /* reinitialize the solution vector */
257: VecGetArray(ctx->U,&u);
258: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
259: u[1] = 1.0;
260: VecRestoreArray(ctx->U,&u);
261: TSSetSolution(ctx->ts,ctx->U);
263: /* reset time */
264: TSSetTime(ctx->ts,0.0);
266: /* reset step counter, this is critical for adjoint solver */
267: TSSetStepNumber(ctx->ts,0);
269: /* reset step size, the step size becomes negative after TSAdjointSolve */
270: TSSetTimeStep(ctx->ts,0.03125);
272: /* reinitialize the integral value */
273: TSGetQuadratureTS(ctx->ts,NULL,&ctx->quadts);
274: TSGetSolution(ctx->quadts,&q);
275: VecSet(q,0.0);
277: if (ctx->sa == SA_TLM) { /* reset the forward sensitivities */
278: TS quadts;
279: Mat sp;
280: PetscScalar val[2];
281: const PetscInt row[]={0,1},col[]={0};
283: TSGetQuadratureTS(ctx->ts,NULL,&quadts);
284: TSForwardGetSensitivities(quadts,NULL,&qgrad);
285: MatZeroEntries(qgrad);
286: TSForwardGetSensitivities(ctx->ts,NULL,&sp);
287: val[0] = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax;
288: val[1] = 0.0;
289: MatSetValues(sp,2,row,1,col,val,INSERT_VALUES);
290: MatAssemblyBegin(sp,MAT_FINAL_ASSEMBLY);
291: MatAssemblyEnd(sp,MAT_FINAL_ASSEMBLY);
292: }
294: /* solve the ODE */
295: TSSolve(ctx->ts,ctx->U);
296: TSGetSolveTime(ctx->ts,&ftime);
297: TSGetStepNumber(ctx->ts,&steps);
299: if (ctx->sa == SA_ADJ) {
300: Vec *lambda,*mu;
301: /* reset the terminal condition for adjoint */
302: TSGetCostGradients(ctx->ts,&nadj,&lambda,&mu);
303: VecGetArray(lambda[0],&y_ptr);
304: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
305: VecRestoreArray(lambda[0],&y_ptr);
306: VecGetArray(mu[0],&x_ptr);
307: x_ptr[0] = -1.0;
308: VecRestoreArray(mu[0],&x_ptr);
310: /* solve the adjont */
311: TSAdjointSolve(ctx->ts);
313: ComputeSensiP(lambda[0],mu[0],ctx);
314: VecCopy(mu[0],G);
315: }
317: if (ctx->sa == SA_TLM) {
318: VecGetArray(G,&x_ptr);
319: MatDenseGetArray(qgrad,&y_ptr);
320: x_ptr[0] = y_ptr[0]-1.;
321: MatDenseRestoreArray(qgrad,&y_ptr);
322: VecRestoreArray(G,&x_ptr);
323: }
325: TSGetSolution(ctx->quadts,&q);
326: VecGetArray(q,&x_ptr);
327: *f = -ctx->Pm + x_ptr[0];
328: VecRestoreArray(q,&x_ptr);
329: return 0;
330: }
332: /*TEST
334: build:
335: requires: !complex !single
337: test:
338: args: -viewer_binary_skip_info -ts_type cn -pc_type lu -tao_monitor
340: test:
341: suffix: 2
342: output_file: output/ex3opt_1.out
343: args: -sa_method tlm -ts_type cn -pc_type lu -tao_monitor
344: TEST*/