Actual source code: ex3opt.c
petsc-3.11.3 2019-06-26
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*/