Actual source code: ex3opt.c


  2: static char help[] = "Finds optimal parameter P_m for the generator system while maintaining generator stability.\n";

  4: /*F

  6: \begin{eqnarray}
  7:                  \frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
  8:                  \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
  9: \end{eqnarray}

 11: F*/

 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, "%" PetscInt_FMT " %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:   PetscMPIInt  size;
 46:   AppCtx       ctx;
 47:   Tao          tao;
 48:   KSP          ksp;
 49:   PC           pc;
 50:   Vec          lambda[1], mu[1], lowerb, upperb;
 51:   PetscBool    printtofile;
 52:   PetscInt     direction[2];
 53:   PetscBool    terminate[2];
 54:   Mat          qgrad; /* Forward sesivitiy */
 55:   Mat          sp;    /* Forward sensitivity matrix */

 57:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 58:      Initialize program
 59:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 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) TaoSetMonitor(tao, (PetscErrorCode(*)(Tao, void *))monitor, (void *)&ctx, PETSC_NULL);
164:   /*
165:      Optimization starts
166:   */
167:   /* Set initial solution guess */
168:   VecCreateSeq(PETSC_COMM_WORLD, 1, &p);
169:   VecGetArray(p, &x_ptr);
170:   x_ptr[0] = ctx.Pm;
171:   VecRestoreArray(p, &x_ptr);

173:   TaoSetSolution(tao, p);
174:   /* Set routine for function and gradient evaluation */
175:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&ctx);

177:   /* Set bounds for the optimization */
178:   VecDuplicate(p, &lowerb);
179:   VecDuplicate(p, &upperb);
180:   VecGetArray(lowerb, &x_ptr);
181:   x_ptr[0] = 0.;
182:   VecRestoreArray(lowerb, &x_ptr);
183:   VecGetArray(upperb, &x_ptr);
184:   x_ptr[0] = 1.1;
185:   VecRestoreArray(upperb, &x_ptr);
186:   TaoSetVariableBounds(tao, lowerb, upperb);

188:   /* Check for any TAO command line options */
189:   TaoSetFromOptions(tao);
190:   TaoGetKSP(tao, &ksp);
191:   if (ksp) {
192:     KSPGetPC(ksp, &pc);
193:     PCSetType(pc, PCNONE);
194:   }

196:   /* SOLVE THE APPLICATION */
197:   TaoSolve(tao);

199:   VecView(p, PETSC_VIEWER_STDOUT_WORLD);

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
203:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204:   MatDestroy(&ctx.Jac);
205:   MatDestroy(&ctx.Jacp);
206:   MatDestroy(&ctx.DRDU);
207:   MatDestroy(&ctx.DRDP);
208:   VecDestroy(&ctx.U);
209:   if (ctx.sa == SA_ADJ) {
210:     VecDestroy(&lambda[0]);
211:     VecDestroy(&mu[0]);
212:   }
213:   if (ctx.sa == SA_TLM) {
214:     MatDestroy(&qgrad);
215:     MatDestroy(&sp);
216:   }
217:   TSDestroy(&ctx.ts);
218:   VecDestroy(&p);
219:   VecDestroy(&lowerb);
220:   VecDestroy(&upperb);
221:   TaoDestroy(&tao);
222:   PetscFinalize();
223:   return 0;
224: }

226: /* ------------------------------------------------------------------ */
227: /*
228:    FormFunctionGradient - Evaluates the function and corresponding gradient.

230:    Input Parameters:
231:    tao - the Tao context
232:    X   - the input vector
233:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

235:    Output Parameters:
236:    f   - the newly evaluated function
237:    G   - the newly evaluated gradient
238: */
239: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx0)
240: {
241:   AppCtx      *ctx = (AppCtx *)ctx0;
242:   PetscInt     nadj;
243:   PetscReal    ftime;
244:   PetscInt     steps;
245:   PetscScalar *u;
246:   PetscScalar *x_ptr, *y_ptr;
247:   Vec          q;
248:   Mat          qgrad;

250:   VecGetArrayRead(P, (const PetscScalar **)&x_ptr);
251:   ctx->Pm = x_ptr[0];
252:   VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr);

254:   /* reinitialize the solution vector */
255:   VecGetArray(ctx->U, &u);
256:   u[0] = PetscAsinScalar(ctx->Pm / ctx->Pmax);
257:   u[1] = 1.0;
258:   VecRestoreArray(ctx->U, &u);
259:   TSSetSolution(ctx->ts, ctx->U);

261:   /* reset time */
262:   TSSetTime(ctx->ts, 0.0);

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

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

270:   /* reinitialize the integral value */
271:   TSGetQuadratureTS(ctx->ts, NULL, &ctx->quadts);
272:   TSGetSolution(ctx->quadts, &q);
273:   VecSet(q, 0.0);

275:   if (ctx->sa == SA_TLM) { /* reset the forward sensitivities */
276:     TS             quadts;
277:     Mat            sp;
278:     PetscScalar    val[2];
279:     const PetscInt row[] = {0, 1}, col[] = {0};

281:     TSGetQuadratureTS(ctx->ts, NULL, &quadts);
282:     TSForwardGetSensitivities(quadts, NULL, &qgrad);
283:     MatZeroEntries(qgrad);
284:     TSForwardGetSensitivities(ctx->ts, NULL, &sp);
285:     val[0] = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax;
286:     val[1] = 0.0;
287:     MatSetValues(sp, 2, row, 1, col, val, INSERT_VALUES);
288:     MatAssemblyBegin(sp, MAT_FINAL_ASSEMBLY);
289:     MatAssemblyEnd(sp, MAT_FINAL_ASSEMBLY);
290:   }

292:   /* solve the ODE */
293:   TSSolve(ctx->ts, ctx->U);
294:   TSGetSolveTime(ctx->ts, &ftime);
295:   TSGetStepNumber(ctx->ts, &steps);

297:   if (ctx->sa == SA_ADJ) {
298:     Vec *lambda, *mu;
299:     /* reset the terminal condition for adjoint */
300:     TSGetCostGradients(ctx->ts, &nadj, &lambda, &mu);
301:     VecGetArray(lambda[0], &y_ptr);
302:     y_ptr[0] = 0.0;
303:     y_ptr[1] = 0.0;
304:     VecRestoreArray(lambda[0], &y_ptr);
305:     VecGetArray(mu[0], &x_ptr);
306:     x_ptr[0] = -1.0;
307:     VecRestoreArray(mu[0], &x_ptr);

309:     /* solve the adjont */
310:     TSAdjointSolve(ctx->ts);

312:     ComputeSensiP(lambda[0], mu[0], ctx);
313:     VecCopy(mu[0], G);
314:   }

316:   if (ctx->sa == SA_TLM) {
317:     VecGetArray(G, &x_ptr);
318:     MatDenseGetArray(qgrad, &y_ptr);
319:     x_ptr[0] = y_ptr[0] - 1.;
320:     MatDenseRestoreArray(qgrad, &y_ptr);
321:     VecRestoreArray(G, &x_ptr);
322:   }

324:   TSGetSolution(ctx->quadts, &q);
325:   VecGetArray(q, &x_ptr);
326:   *f = -ctx->Pm + x_ptr[0];
327:   VecRestoreArray(q, &x_ptr);
328:   return 0;
329: }

331: /*TEST

333:    build:
334:       requires: !complex !single

336:    test:
337:       args: -viewer_binary_skip_info -ts_type cn -pc_type lu -tao_monitor

339:    test:
340:       suffix: 2
341:       output_file: output/ex3opt_1.out
342:       args: -sa_method tlm -ts_type cn -pc_type lu -tao_monitor
343: TEST*/