Actual source code: ex3opt.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

  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>

 22: typedef enum {SA_ADJ, SA_TLM} SAMethod;
 23: static const char *const SAMethods[] = {"ADJ","TLM","SAMethod","SA_",0};

 25: typedef struct {
 26:   PetscScalar H,D,omega_b,omega_s,Pmax,Pmax_ini,Pm,E,V,X,u_s,c;
 27:   PetscInt    beta;
 28:   PetscReal   tf,tcl;
 29:   /* Solver context */
 30:   TS          ts;
 31:   Vec         U;    /* solution will be stored here */
 32:   Mat         Jac;  /* Jacobian matrix */
 33:   Mat         Jacp; /* Jacobianp matrix */
 34:   SAMethod    sa;
 35: } AppCtx;

 37: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);

 39: /* Event check */
 40: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec X,PetscScalar *fvalue,void *ctx)
 41: {
 42:   AppCtx *user=(AppCtx*)ctx;

 45:   /* Event for fault-on time */
 46:   fvalue[0] = t - user->tf;
 47:   /* Event for fault-off time */
 48:   fvalue[1] = t - user->tcl;

 50:   return(0);
 51: }

 53: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec X,PetscBool forwardsolve,void* ctx)
 54: {
 55:   AppCtx *user=(AppCtx*)ctx;

 58:   if (event_list[0] == 0) {
 59:     if (forwardsolve) user->Pmax = 0.0; /* Apply disturbance - this is done by setting Pmax = 0 */
 60:     else user->Pmax = user->Pmax_ini; /* Going backward, reversal of event */
 61:   } else if(event_list[0] == 1) {
 62:     if (forwardsolve) user->Pmax = user->Pmax_ini; /* Remove the fault  - this is done by setting Pmax = Pmax_ini */
 63:     else user->Pmax = 0.0; /* Going backward, reversal of event */
 64:   }
 65:   return(0);
 66: }

 68: /*
 69:      Defines the ODE passed to the ODE solver
 70: */
 71: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx)
 72: {
 73:   PetscErrorCode    ierr;
 74:   PetscScalar       *f,Pmax;
 75:   const PetscScalar *u,*udot;

 78:   /*  The next three lines allow us to access the entries of the vectors directly */
 79:   VecGetArrayRead(U,&u);
 80:   VecGetArrayRead(Udot,&udot);
 81:   VecGetArray(F,&f);
 82:   Pmax = ctx->Pmax;
 83:   f[0] = udot[0] - ctx->omega_b*(u[1] - ctx->omega_s);
 84:   f[1] = 2.0*ctx->H/ctx->omega_s*udot[1] +  Pmax*PetscSinScalar(u[0]) + ctx->D*(u[1] - ctx->omega_s)- ctx->Pm;

 86:   VecRestoreArrayRead(U,&u);
 87:   VecRestoreArrayRead(Udot,&udot);
 88:   VecRestoreArray(F,&f);
 89:   return(0);
 90: }

 92: /*
 93:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
 94: */
 95: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx)
 96: {
 97:   PetscErrorCode    ierr;
 98:   PetscInt          rowcol[] = {0,1};
 99:   PetscScalar       J[2][2],Pmax;
100:   const PetscScalar *u,*udot;

103:   VecGetArrayRead(U,&u);
104:   VecGetArrayRead(Udot,&udot);
105:   Pmax = ctx->Pmax;
106:   J[0][0] = a;                       J[0][1] = -ctx->omega_b;
107:   J[1][1] = 2.0*ctx->H/ctx->omega_s*a + ctx->D;   J[1][0] = Pmax*PetscCosScalar(u[0]);

109:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
110:   VecRestoreArrayRead(U,&u);
111:   VecRestoreArrayRead(Udot,&udot);

113:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
114:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
115:   if (A != B) {
116:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
117:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
118:   }
119:   return(0);
120: }

122: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
123: {
125:   PetscInt       row[] = {0,1},col[]={0};
126:   PetscScalar    J[2][1];

129:   J[0][0] = 0;
130:   J[1][0] = 1.;
131:   MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
132:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
133:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
134:   return(0);
135: }

137: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
138: {
139:   PetscErrorCode    ierr;
140:   PetscScalar       *r;
141:   const PetscScalar *u;

144:   VecGetArrayRead(U,&u);
145:   VecGetArray(R,&r);
146:   r[0] = ctx->c*PetscPowScalarInt(PetscMax(0.,u[0]-ctx->u_s),ctx->beta);
147:   VecRestoreArray(R,&r);
148:   VecRestoreArrayRead(U,&u);
149:   return(0);
150: }

152: static PetscErrorCode DRDYFunction(TS ts,PetscReal t,Vec U,Vec *drdy,AppCtx *ctx)
153: {
154:   PetscErrorCode     ierr;
155:   PetscScalar        *ry;
156:   const PetscScalar  *u;

159:   VecGetArrayRead(U,&u);
160:   VecGetArray(drdy[0],&ry);
161:   ry[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0.,u[0]-ctx->u_s),ctx->beta-1);
162:   VecRestoreArray(drdy[0],&ry);
163:   VecRestoreArrayRead(U,&u);
164:   return(0);
165: }

167: static PetscErrorCode DRDPFunction(TS ts,PetscReal t,Vec U,Vec *drdp,AppCtx *ctx)
168: {
170:   PetscScalar    *rp;

173:   VecGetArray(drdp[0],&rp);
174:   rp[0] = 0.;
175:   VecRestoreArray(drdp[0],&rp);
176:   return(0);
177: }

179: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
180: {
181:   PetscErrorCode    ierr;
182:   PetscScalar       *y,sensip;
183:   const PetscScalar *x;

186:   VecGetArrayRead(lambda,&x);
187:   VecGetArray(mu,&y);
188:   sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
189:   y[0] = sensip;
190:   VecRestoreArray(mu,&y);
191:   VecRestoreArrayRead(lambda,&x);
192:   return(0);
193: }

195: PetscErrorCode monitor(Tao tao,AppCtx *ctx)
196: {
197:   FILE               *fp;
198:   PetscInt           iterate;
199:   PetscReal          f,gnorm,cnorm,xdiff;
200:   TaoConvergedReason reason;
201:   PetscErrorCode     ierr;

204:   TaoGetSolutionStatus(tao,&iterate,&f,&gnorm,&cnorm,&xdiff,&reason);

206:   fp = fopen("ex3opt_conv.out","a");
207:   PetscFPrintf(PETSC_COMM_WORLD,fp,"%D %g\n",iterate,(double)gnorm);
208:   fclose(fp);
209:   return(0);
210: }

212: int main(int argc,char **argv)
213: {
214:   Vec                p;
215:   PetscScalar        *x_ptr;
216:   PetscErrorCode     ierr;
217:   PetscMPIInt        size;
218:   AppCtx             ctx;
219:   Tao                tao;
220:   KSP                ksp;
221:   PC                 pc;
222:   Vec                lambda[1],mu[1],lowerb,upperb;
223:   PetscBool          printtofile;
224:   PetscInt           direction[2];
225:   PetscBool          terminate[2];
226:   Vec                qgrad[1];      /* Forward sesivitiy */
227:   Mat                sp;            /* Forward sensitivity matrix */

229:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230:      Initialize program
231:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
232:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
234:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
235:   if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");

237:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238:     Set runtime options
239:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
240:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
241:   {
242:     ctx.beta    = 2;
243:     ctx.c       = 10000.0;
244:     ctx.u_s     = 1.0;
245:     ctx.omega_s = 1.0;
246:     ctx.omega_b = 120.0*PETSC_PI;
247:     ctx.H       = 5.0;
248:     PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
249:     ctx.D       = 5.0;
250:     PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
251:     ctx.E       = 1.1378;
252:     ctx.V       = 1.0;
253:     ctx.X       = 0.545;
254:     ctx.Pmax    = ctx.E*ctx.V/ctx.X;;
255:     ctx.Pmax_ini = ctx.Pmax;
256:     PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
257:     ctx.Pm      = 1.06;
258:     PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
259:     ctx.tf      = 0.1;
260:     ctx.tcl     = 0.2;
261:     PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
262:     PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
263:     printtofile = PETSC_FALSE;
264:     PetscOptionsBool("-printtofile","Print convergence results to file","",printtofile,&printtofile,NULL);
265:     ctx.sa      = SA_ADJ;
266:     PetscOptionsEnum("-sa_method","Sensitivity analysis method (adj or tlm)","",SAMethods,(PetscEnum)ctx.sa,(PetscEnum*)&ctx.sa,NULL);
267:   }
268:   PetscOptionsEnd();

270:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271:     Create necessary matrix and vectors
272:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
273:   MatCreate(PETSC_COMM_WORLD,&ctx.Jac);
274:   MatSetSizes(ctx.Jac,2,2,PETSC_DETERMINE,PETSC_DETERMINE);
275:   MatSetType(ctx.Jac,MATDENSE);
276:   MatSetFromOptions(ctx.Jac);
277:   MatSetUp(ctx.Jac);
278:   MatCreate(PETSC_COMM_WORLD,&ctx.Jacp);
279:   MatSetSizes(ctx.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
280:   MatSetFromOptions(ctx.Jacp);
281:   MatSetUp(ctx.Jacp);
282:   MatCreateVecs(ctx.Jac,&ctx.U,NULL);

284:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
285:      Create timestepping solver context
286:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
287:   TSCreate(PETSC_COMM_WORLD,&ctx.ts);
288:   TSSetProblemType(ctx.ts,TS_NONLINEAR);
289:   TSSetType(ctx.ts,TSCN);
290:   TSSetIFunction(ctx.ts,NULL,(TSIFunction) IFunction,&ctx);
291:   TSSetIJacobian(ctx.ts,ctx.Jac,ctx.Jac,(TSIJacobian)IJacobian,&ctx);
292:   TSSetRHSJacobianP(ctx.ts,ctx.Jacp,RHSJacobianP,&ctx);

294:   if (ctx.sa == SA_ADJ) {
295:     MatCreateVecs(ctx.Jac,&lambda[0],NULL);
296:     MatCreateVecs(ctx.Jacp,&mu[0],NULL);
297:     TSSetSaveTrajectory(ctx.ts);
298:     TSSetCostGradients(ctx.ts,1,lambda,mu);
299:     TSSetCostIntegrand(ctx.ts,1,NULL,(PetscErrorCode (*)(TS,PetscReal,Vec,Vec,void*))CostIntegrand,
300:                               (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDYFunction,
301:                               (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDPFunction,PETSC_TRUE,&ctx);
302:   }
303:   if (ctx.sa == SA_TLM) {
304:     VecCreate(PETSC_COMM_WORLD,&qgrad[0]);
305:     VecSetSizes(qgrad[0],PETSC_DECIDE,1);
306:     VecSetFromOptions(qgrad[0]);
307:     MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&sp);
308:     TSForwardSetSensitivities(ctx.ts,1,sp);
309:     TSForwardSetIntegralGradients(ctx.ts,1,qgrad);
310:     TSSetCostIntegrand(ctx.ts,1,NULL,(PetscErrorCode (*)(TS,PetscReal,Vec,Vec,void*))CostIntegrand,
311:                               (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDYFunction,
312:                               (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDPFunction,PETSC_TRUE,&ctx);
313:   }

315:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
316:      Set solver options
317:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
318:   TSSetMaxTime(ctx.ts,1.0);
319:   TSSetExactFinalTime(ctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
320:   TSSetTimeStep(ctx.ts,0.03125);
321:   TSSetFromOptions(ctx.ts);

323:   direction[0] = direction[1] = 1;
324:   terminate[0] = terminate[1] = PETSC_FALSE;
325:   TSSetEventHandler(ctx.ts,2,direction,terminate,EventFunction,PostEventFunction,&ctx);

327:   /* Create TAO solver and set desired solution method */
328:   TaoCreate(PETSC_COMM_WORLD,&tao);
329:   TaoSetType(tao,TAOBLMVM);
330:   if(printtofile) {
331:     TaoSetMonitor(tao,(PetscErrorCode (*)(Tao, void*))monitor,(void *)&ctx,PETSC_NULL);
332:   }
333:   /*
334:      Optimization starts
335:   */
336:   /* Set initial solution guess */
337:   VecCreateSeq(PETSC_COMM_WORLD,1,&p);
338:   VecGetArray(p,&x_ptr);
339:   x_ptr[0] = ctx.Pm;
340:   VecRestoreArray(p,&x_ptr);

342:   TaoSetInitialVector(tao,p);
343:   /* Set routine for function and gradient evaluation */
344:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&ctx);

346:   /* Set bounds for the optimization */
347:   VecDuplicate(p,&lowerb);
348:   VecDuplicate(p,&upperb);
349:   VecGetArray(lowerb,&x_ptr);
350:   x_ptr[0] = 0.;
351:   VecRestoreArray(lowerb,&x_ptr);
352:   VecGetArray(upperb,&x_ptr);
353:   x_ptr[0] = 1.1;
354:   VecRestoreArray(upperb,&x_ptr);
355:   TaoSetVariableBounds(tao,lowerb,upperb);

357:   /* Check for any TAO command line options */
358:   TaoSetFromOptions(tao);
359:   TaoGetKSP(tao,&ksp);
360:   if (ksp) {
361:     KSPGetPC(ksp,&pc);
362:     PCSetType(pc,PCNONE);
363:   }

365:   /* SOLVE THE APPLICATION */
366:   TaoSolve(tao);

368:   VecView(p,PETSC_VIEWER_STDOUT_WORLD);

370:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
372:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
373:   MatDestroy(&ctx.Jac);
374:   MatDestroy(&ctx.Jacp);
375:   VecDestroy(&ctx.U);
376:   if (ctx.sa == SA_ADJ) {
377:     VecDestroy(&lambda[0]);
378:     VecDestroy(&mu[0]);
379:   }
380:   if (ctx.sa == SA_TLM) {
381:     VecDestroy(&qgrad[0]);
382:     MatDestroy(&sp);
383:   }
384:   TSDestroy(&ctx.ts);
385:   VecDestroy(&p);
386:   VecDestroy(&lowerb);
387:   VecDestroy(&upperb);
388:   TaoDestroy(&tao);
389:   PetscFinalize();
390:   return ierr;
391: }

393: /* ------------------------------------------------------------------ */
394: /*
395:    FormFunctionGradient - Evaluates the function and corresponding gradient.

397:    Input Parameters:
398:    tao - the Tao context
399:    X   - the input vector
400:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

402:    Output Parameters:
403:    f   - the newly evaluated function
404:    G   - the newly evaluated gradient
405: */
406: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx0)
407: {
408:   AppCtx         *ctx = (AppCtx*)ctx0;
409:   PetscInt       nadj;
410:   PetscReal      ftime;
411:   PetscInt       steps;
412:   PetscScalar    *u;
413:   PetscScalar    *x_ptr,*y_ptr;
414:   Vec            q,*qgrad;

417:   VecGetArrayRead(P,(const PetscScalar**)&x_ptr);
418:   ctx->Pm = x_ptr[0];
419:   VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);

421:   /* reinitialize the solution vector */
422:   VecGetArray(ctx->U,&u);
423:   u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
424:   u[1] = 1.0;
425:   VecRestoreArray(ctx->U,&u);
426:   TSSetSolution(ctx->ts,ctx->U);

428:   /* reset time */
429:   TSSetTime(ctx->ts,0.0);

431:   /* reset step counter, this is critical for adjoint solver */
432:   TSSetStepNumber(ctx->ts,0);

434:   /* reset step size, the step size becomes negative after TSAdjointSolve */
435:   TSSetTimeStep(ctx->ts,0.03125);

437:   /* reinitialize the integral value */
438:   TSGetCostIntegral(ctx->ts,&q);
439:   VecSet(q,0.0);

441:   if (ctx->sa == SA_TLM) { /* reset the forward sensitivities */
442:     Mat            sp;
443:     PetscScalar    val[2];
444:     const PetscInt row[]={0,1},col[]={0};

446:     TSForwardGetIntegralGradients(ctx->ts,NULL,&qgrad);
447:     VecSet(qgrad[0],0.0);
448:     TSForwardGetSensitivities(ctx->ts,NULL,&sp);
449:     val[0] = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax;
450:     val[1] = 0.0;
451:     MatSetValues(sp,2,row,1,col,val,INSERT_VALUES);
452:     MatAssemblyBegin(sp,MAT_FINAL_ASSEMBLY);
453:     MatAssemblyEnd(sp,MAT_FINAL_ASSEMBLY);
454:   }

456:   /* reset the trajectory */
457:   TSResetTrajectory(ctx->ts);

459:   /* solve the ODE */
460:   TSSolve(ctx->ts,ctx->U);
461:   TSGetSolveTime(ctx->ts,&ftime);
462:   TSGetStepNumber(ctx->ts,&steps);

464:   if (ctx->sa == SA_ADJ) {
465:     Vec *lambda,*mu;
466:     /* reset the terminal condition for adjoint */
467:     TSGetCostGradients(ctx->ts,&nadj,&lambda,&mu);
468:     VecGetArray(lambda[0],&y_ptr);
469:     y_ptr[0] = 0.0; y_ptr[1] = 0.0;
470:     VecRestoreArray(lambda[0],&y_ptr);
471:     VecGetArray(mu[0],&x_ptr);
472:     x_ptr[0] = -1.0;
473:     VecRestoreArray(mu[0],&x_ptr);

475:     /* solve the adjont */
476:     TSAdjointSolve(ctx->ts);

478:     ComputeSensiP(lambda[0],mu[0],ctx);
479:     VecCopy(mu[0],G);
480:   }

482:   if (ctx->sa == SA_TLM) {
483:     VecGetArray(G,&x_ptr);
484:     VecGetArray(qgrad[0],&y_ptr);
485:     x_ptr[0] = y_ptr[0]-1.;
486:     VecRestoreArray(qgrad[0],&y_ptr);
487:     VecRestoreArray(G,&x_ptr);
488:   }

490:   TSGetCostIntegral(ctx->ts,&q);
491:   VecGetArray(q,&x_ptr);
492:   *f   = -ctx->Pm + x_ptr[0];
493:   VecRestoreArray(q,&x_ptr);
494:   return 0;
495: }

497: /*TEST

499:    build:
500:       requires: !complex !single

502:    test:
503:       args: -viewer_binary_skip_info -ts_type cn -pc_type lu -tao_monitor

505:    test:
506:       suffix: 2
507:       output_file: output/ex3opt_1.out
508:       args: -sa_method tlm -ts_type cn -pc_type lu -tao_monitor
509: TEST*/