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