1 2 static char help[] = "Finds optimal parameter P_m for the generator system while maintaining generator stability.\n"; 3 4 /*F 5 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} 10 11 F*/ 12 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 */ 18 19 #include <petsctao.h> 20 #include <petscts.h> 21 #include "ex3.h" 22 23 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 24 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; 32 33 PetscFunctionBeginUser; 34 ierr = TaoGetSolutionStatus(tao,&iterate,&f,&gnorm,&cnorm,&xdiff,&reason);CHKERRQ(ierr); 35 36 fp = fopen("ex3opt_conv.out","a"); 37 ierr = PetscFPrintf(PETSC_COMM_WORLD,fp,"%D %g\n",iterate,(double)gnorm);CHKERRQ(ierr); 38 fclose(fp); 39 PetscFunctionReturn(0); 40 } 41 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 */ 58 59 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 60 Initialize program 61 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 62 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 63 PetscFunctionBeginUser; 64 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 65 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 66 67 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 68 Set runtime options 69 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 70 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");CHKERRQ(ierr); 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 ierr = PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);CHKERRQ(ierr); 79 ctx.D = 5.0; 80 ierr = PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);CHKERRQ(ierr); 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 ierr = PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);CHKERRQ(ierr); 87 ctx.Pm = 1.06; 88 ierr = PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);CHKERRQ(ierr); 89 ctx.tf = 0.1; 90 ctx.tcl = 0.2; 91 ierr = PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);CHKERRQ(ierr); 92 ierr = PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);CHKERRQ(ierr); 93 printtofile = PETSC_FALSE; 94 ierr = PetscOptionsBool("-printtofile","Print convergence results to file","",printtofile,&printtofile,NULL);CHKERRQ(ierr); 95 ctx.sa = SA_ADJ; 96 ierr = PetscOptionsEnum("-sa_method","Sensitivity analysis method (adj or tlm)","",SAMethods,(PetscEnum)ctx.sa,(PetscEnum*)&ctx.sa,NULL);CHKERRQ(ierr); 97 } 98 ierr = PetscOptionsEnd();CHKERRQ(ierr); 99 100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 Create necessary matrix and vectors 102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 103 ierr = MatCreate(PETSC_COMM_WORLD,&ctx.Jac);CHKERRQ(ierr); 104 ierr = MatSetSizes(ctx.Jac,2,2,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 105 ierr = MatSetType(ctx.Jac,MATDENSE);CHKERRQ(ierr); 106 ierr = MatSetFromOptions(ctx.Jac);CHKERRQ(ierr); 107 ierr = MatSetUp(ctx.Jac);CHKERRQ(ierr); 108 ierr = MatCreate(PETSC_COMM_WORLD,&ctx.Jacp);CHKERRQ(ierr); 109 ierr = MatSetSizes(ctx.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);CHKERRQ(ierr); 110 ierr = MatSetFromOptions(ctx.Jacp);CHKERRQ(ierr); 111 ierr = MatSetUp(ctx.Jacp);CHKERRQ(ierr); 112 ierr = MatCreateVecs(ctx.Jac,&ctx.U,NULL);CHKERRQ(ierr); 113 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&ctx.DRDP);CHKERRQ(ierr); 114 ierr = MatSetUp(ctx.DRDP);CHKERRQ(ierr); 115 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&ctx.DRDU);CHKERRQ(ierr); 116 ierr = MatSetUp(ctx.DRDU);CHKERRQ(ierr); 117 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Create timestepping solver context 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 ierr = TSCreate(PETSC_COMM_WORLD,&ctx.ts);CHKERRQ(ierr); 122 ierr = TSSetProblemType(ctx.ts,TS_NONLINEAR);CHKERRQ(ierr); 123 ierr = TSSetType(ctx.ts,TSCN);CHKERRQ(ierr); 124 ierr = TSSetRHSFunction(ctx.ts,NULL,(TSRHSFunction)RHSFunction,&ctx);CHKERRQ(ierr); 125 ierr = TSSetRHSJacobian(ctx.ts,ctx.Jac,ctx.Jac,(TSRHSJacobian)RHSJacobian,&ctx);CHKERRQ(ierr); 126 ierr = TSSetRHSJacobianP(ctx.ts,ctx.Jacp,RHSJacobianP,&ctx);CHKERRQ(ierr); 127 128 if (ctx.sa == SA_ADJ) { 129 ierr = MatCreateVecs(ctx.Jac,&lambda[0],NULL);CHKERRQ(ierr); 130 ierr = MatCreateVecs(ctx.Jacp,&mu[0],NULL);CHKERRQ(ierr); 131 ierr = TSSetSaveTrajectory(ctx.ts);CHKERRQ(ierr); 132 ierr = TSSetCostGradients(ctx.ts,1,lambda,mu);CHKERRQ(ierr); 133 ierr = TSCreateQuadratureTS(ctx.ts,PETSC_FALSE,&ctx.quadts);CHKERRQ(ierr); 134 ierr = TSSetRHSFunction(ctx.quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);CHKERRQ(ierr); 135 ierr = TSSetRHSJacobian(ctx.quadts,ctx.DRDU,ctx.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);CHKERRQ(ierr); 136 ierr = TSSetRHSJacobianP(ctx.quadts,ctx.DRDP,DRDPJacobianTranspose,&ctx);CHKERRQ(ierr); 137 } 138 if (ctx.sa == SA_TLM) { 139 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&qgrad);CHKERRQ(ierr); 140 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&sp);CHKERRQ(ierr); 141 ierr = TSForwardSetSensitivities(ctx.ts,1,sp);CHKERRQ(ierr); 142 ierr = TSCreateQuadratureTS(ctx.ts,PETSC_TRUE,&ctx.quadts);CHKERRQ(ierr); 143 ierr = TSForwardSetSensitivities(ctx.quadts,1,qgrad);CHKERRQ(ierr); 144 ierr = TSSetRHSFunction(ctx.quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);CHKERRQ(ierr); 145 ierr = TSSetRHSJacobian(ctx.quadts,ctx.DRDU,ctx.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);CHKERRQ(ierr); 146 ierr = TSSetRHSJacobianP(ctx.quadts,ctx.DRDP,DRDPJacobianTranspose,&ctx);CHKERRQ(ierr); 147 } 148 149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 Set solver options 151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152 ierr = TSSetMaxTime(ctx.ts,1.0);CHKERRQ(ierr); 153 ierr = TSSetExactFinalTime(ctx.ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 154 ierr = TSSetTimeStep(ctx.ts,0.03125);CHKERRQ(ierr); 155 ierr = TSSetFromOptions(ctx.ts);CHKERRQ(ierr); 156 157 direction[0] = direction[1] = 1; 158 terminate[0] = terminate[1] = PETSC_FALSE; 159 ierr = TSSetEventHandler(ctx.ts,2,direction,terminate,EventFunction,PostEventFunction,&ctx);CHKERRQ(ierr); 160 161 /* Create TAO solver and set desired solution method */ 162 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 163 ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr); 164 if (printtofile) { 165 ierr = TaoSetMonitor(tao,(PetscErrorCode (*)(Tao, void*))monitor,(void *)&ctx,PETSC_NULL);CHKERRQ(ierr); 166 } 167 /* 168 Optimization starts 169 */ 170 /* Set initial solution guess */ 171 ierr = VecCreateSeq(PETSC_COMM_WORLD,1,&p);CHKERRQ(ierr); 172 ierr = VecGetArray(p,&x_ptr);CHKERRQ(ierr); 173 x_ptr[0] = ctx.Pm; 174 ierr = VecRestoreArray(p,&x_ptr);CHKERRQ(ierr); 175 176 ierr = TaoSetInitialVector(tao,p);CHKERRQ(ierr); 177 /* Set routine for function and gradient evaluation */ 178 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&ctx);CHKERRQ(ierr); 179 180 /* Set bounds for the optimization */ 181 ierr = VecDuplicate(p,&lowerb);CHKERRQ(ierr); 182 ierr = VecDuplicate(p,&upperb);CHKERRQ(ierr); 183 ierr = VecGetArray(lowerb,&x_ptr);CHKERRQ(ierr); 184 x_ptr[0] = 0.; 185 ierr = VecRestoreArray(lowerb,&x_ptr);CHKERRQ(ierr); 186 ierr = VecGetArray(upperb,&x_ptr);CHKERRQ(ierr); 187 x_ptr[0] = 1.1; 188 ierr = VecRestoreArray(upperb,&x_ptr);CHKERRQ(ierr); 189 ierr = TaoSetVariableBounds(tao,lowerb,upperb);CHKERRQ(ierr); 190 191 /* Check for any TAO command line options */ 192 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 193 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 194 if (ksp) { 195 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 196 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 197 } 198 199 /* SOLVE THE APPLICATION */ 200 ierr = TaoSolve(tao);CHKERRQ(ierr); 201 202 ierr = VecView(p,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 203 204 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 205 Free work space. All PETSc objects should be destroyed when they are no longer needed. 206 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 207 ierr = MatDestroy(&ctx.Jac);CHKERRQ(ierr); 208 ierr = MatDestroy(&ctx.Jacp);CHKERRQ(ierr); 209 ierr = MatDestroy(&ctx.DRDU);CHKERRQ(ierr); 210 ierr = MatDestroy(&ctx.DRDP);CHKERRQ(ierr); 211 ierr = VecDestroy(&ctx.U);CHKERRQ(ierr); 212 if (ctx.sa == SA_ADJ) { 213 ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr); 214 ierr = VecDestroy(&mu[0]);CHKERRQ(ierr); 215 } 216 if (ctx.sa == SA_TLM) { 217 ierr = MatDestroy(&qgrad);CHKERRQ(ierr); 218 ierr = MatDestroy(&sp);CHKERRQ(ierr); 219 } 220 ierr = TSDestroy(&ctx.ts);CHKERRQ(ierr); 221 ierr = VecDestroy(&p);CHKERRQ(ierr); 222 ierr = VecDestroy(&lowerb);CHKERRQ(ierr); 223 ierr = VecDestroy(&upperb);CHKERRQ(ierr); 224 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 225 ierr = PetscFinalize(); 226 return ierr; 227 } 228 229 /* ------------------------------------------------------------------ */ 230 /* 231 FormFunctionGradient - Evaluates the function and corresponding gradient. 232 233 Input Parameters: 234 tao - the Tao context 235 X - the input vector 236 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine() 237 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; 252 PetscErrorCode ierr; 253 254 ierr = VecGetArrayRead(P,(const PetscScalar**)&x_ptr);CHKERRQ(ierr); 255 ctx->Pm = x_ptr[0]; 256 ierr = VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);CHKERRQ(ierr); 257 258 /* reinitialize the solution vector */ 259 ierr = VecGetArray(ctx->U,&u);CHKERRQ(ierr); 260 u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax); 261 u[1] = 1.0; 262 ierr = VecRestoreArray(ctx->U,&u);CHKERRQ(ierr); 263 ierr = TSSetSolution(ctx->ts,ctx->U);CHKERRQ(ierr); 264 265 /* reset time */ 266 ierr = TSSetTime(ctx->ts,0.0);CHKERRQ(ierr); 267 268 /* reset step counter, this is critical for adjoint solver */ 269 ierr = TSSetStepNumber(ctx->ts,0);CHKERRQ(ierr); 270 271 /* reset step size, the step size becomes negative after TSAdjointSolve */ 272 ierr = TSSetTimeStep(ctx->ts,0.03125);CHKERRQ(ierr); 273 274 /* reinitialize the integral value */ 275 ierr = TSGetQuadratureTS(ctx->ts,NULL,&ctx->quadts);CHKERRQ(ierr); 276 ierr = TSGetSolution(ctx->quadts,&q);CHKERRQ(ierr); 277 ierr = VecSet(q,0.0);CHKERRQ(ierr); 278 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}; 284 285 ierr = TSGetQuadratureTS(ctx->ts,NULL,&quadts);CHKERRQ(ierr); 286 ierr = TSForwardGetSensitivities(quadts,NULL,&qgrad);CHKERRQ(ierr); 287 ierr = MatZeroEntries(qgrad);CHKERRQ(ierr); 288 ierr = TSForwardGetSensitivities(ctx->ts,NULL,&sp);CHKERRQ(ierr); 289 val[0] = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax; 290 val[1] = 0.0; 291 ierr = MatSetValues(sp,2,row,1,col,val,INSERT_VALUES);CHKERRQ(ierr); 292 ierr = MatAssemblyBegin(sp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 293 ierr = MatAssemblyEnd(sp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 294 } 295 296 /* solve the ODE */ 297 ierr = TSSolve(ctx->ts,ctx->U);CHKERRQ(ierr); 298 ierr = TSGetSolveTime(ctx->ts,&ftime);CHKERRQ(ierr); 299 ierr = TSGetStepNumber(ctx->ts,&steps);CHKERRQ(ierr); 300 301 if (ctx->sa == SA_ADJ) { 302 Vec *lambda,*mu; 303 /* reset the terminal condition for adjoint */ 304 ierr = TSGetCostGradients(ctx->ts,&nadj,&lambda,&mu);CHKERRQ(ierr); 305 ierr = VecGetArray(lambda[0],&y_ptr);CHKERRQ(ierr); 306 y_ptr[0] = 0.0; y_ptr[1] = 0.0; 307 ierr = VecRestoreArray(lambda[0],&y_ptr);CHKERRQ(ierr); 308 ierr = VecGetArray(mu[0],&x_ptr);CHKERRQ(ierr); 309 x_ptr[0] = -1.0; 310 ierr = VecRestoreArray(mu[0],&x_ptr);CHKERRQ(ierr); 311 312 /* solve the adjont */ 313 ierr = TSAdjointSolve(ctx->ts);CHKERRQ(ierr); 314 315 ierr = ComputeSensiP(lambda[0],mu[0],ctx);CHKERRQ(ierr); 316 ierr = VecCopy(mu[0],G);CHKERRQ(ierr); 317 } 318 319 if (ctx->sa == SA_TLM) { 320 ierr = VecGetArray(G,&x_ptr);CHKERRQ(ierr); 321 ierr = MatDenseGetArray(qgrad,&y_ptr);CHKERRQ(ierr); 322 x_ptr[0] = y_ptr[0]-1.; 323 ierr = MatDenseRestoreArray(qgrad,&y_ptr);CHKERRQ(ierr); 324 ierr = VecRestoreArray(G,&x_ptr);CHKERRQ(ierr); 325 } 326 327 ierr = TSGetSolution(ctx->quadts,&q);CHKERRQ(ierr); 328 ierr = VecGetArray(q,&x_ptr);CHKERRQ(ierr); 329 *f = -ctx->Pm + x_ptr[0]; 330 ierr = VecRestoreArray(q,&x_ptr);CHKERRQ(ierr); 331 return 0; 332 } 333 334 /*TEST 335 336 build: 337 requires: !complex !single 338 339 test: 340 args: -viewer_binary_skip_info -ts_type cn -pc_type lu -tao_monitor 341 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*/ 347