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 Solve the same optimization problem as in ex3opt.c. 15 Use finite difference to approximate the gradients. 16 */ 17 #include <petsctao.h> 18 #include <petscts.h> 19 #include "ex3.h" 20 21 PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*); 22 23 PetscErrorCode monitor(Tao tao,AppCtx *ctx) 24 { 25 FILE *fp; 26 PetscInt iterate; 27 PetscReal f,gnorm,cnorm,xdiff; 28 Vec X,G; 29 const PetscScalar *x,*g; 30 TaoConvergedReason reason; 31 PetscErrorCode ierr; 32 33 PetscFunctionBeginUser; 34 ierr = TaoGetSolutionStatus(tao,&iterate,&f,&gnorm,&cnorm,&xdiff,&reason);CHKERRQ(ierr); 35 ierr = TaoGetSolutionVector(tao,&X);CHKERRQ(ierr); 36 ierr = TaoGetGradientVector(tao,&G);CHKERRQ(ierr); 37 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 38 ierr = VecGetArrayRead(G,&g);CHKERRQ(ierr); 39 fp = fopen("ex3opt_fd_conv.out","a"); 40 ierr = PetscFPrintf(PETSC_COMM_WORLD,fp,"%d %g %.12lf %.12lf\n",iterate,gnorm,x[0],g[0]);CHKERRQ(ierr); 41 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 42 ierr = VecRestoreArrayRead(G,&g);CHKERRQ(ierr); 43 fclose(fp); 44 PetscFunctionReturn(0); 45 } 46 47 int main(int argc,char **argv) 48 { 49 Vec p; 50 PetscScalar *x_ptr; 51 PetscErrorCode ierr; 52 PetscMPIInt size; 53 AppCtx ctx; 54 Vec lowerb,upperb; 55 Tao tao; 56 KSP ksp; 57 PC pc; 58 PetscBool printtofile; 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 } 96 ierr = PetscOptionsEnd();CHKERRQ(ierr); 97 98 /* Create TAO solver and set desired solution method */ 99 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 100 ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr); 101 if(printtofile) { 102 ierr = TaoSetMonitor(tao,(PetscErrorCode (*)(Tao, void*))monitor,(void *)&ctx,PETSC_NULL);CHKERRQ(ierr); 103 } 104 ierr = TaoSetMaximumIterations(tao,30);CHKERRQ(ierr); 105 /* 106 Optimization starts 107 */ 108 /* Set initial solution guess */ 109 ierr = VecCreateSeq(PETSC_COMM_WORLD,1,&p);CHKERRQ(ierr); 110 ierr = VecGetArray(p,&x_ptr);CHKERRQ(ierr); 111 x_ptr[0] = ctx.Pm; 112 ierr = VecRestoreArray(p,&x_ptr);CHKERRQ(ierr); 113 114 ierr = TaoSetInitialVector(tao,p);CHKERRQ(ierr); 115 /* Set routine for function and gradient evaluation */ 116 ierr = TaoSetObjectiveRoutine(tao,FormFunction,(void *)&ctx);CHKERRQ(ierr); 117 ierr = TaoSetGradientRoutine(tao,TaoDefaultComputeGradient,(void *)&ctx);CHKERRQ(ierr); 118 119 /* Set bounds for the optimization */ 120 ierr = VecDuplicate(p,&lowerb);CHKERRQ(ierr); 121 ierr = VecDuplicate(p,&upperb);CHKERRQ(ierr); 122 ierr = VecGetArray(lowerb,&x_ptr);CHKERRQ(ierr); 123 x_ptr[0] = 0.; 124 ierr = VecRestoreArray(lowerb,&x_ptr);CHKERRQ(ierr); 125 ierr = VecGetArray(upperb,&x_ptr);CHKERRQ(ierr); 126 x_ptr[0] = 1.1; 127 ierr = VecRestoreArray(upperb,&x_ptr);CHKERRQ(ierr); 128 ierr = TaoSetVariableBounds(tao,lowerb,upperb);CHKERRQ(ierr); 129 130 /* Check for any TAO command line options */ 131 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 132 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 133 if (ksp) { 134 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 135 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 136 } 137 138 /* SOLVE THE APPLICATION */ 139 ierr = TaoSolve(tao); CHKERRQ(ierr); 140 141 ierr = VecView(p,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 142 143 /* Free TAO data structures */ 144 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 145 ierr = VecDestroy(&p);CHKERRQ(ierr); 146 ierr = VecDestroy(&lowerb);CHKERRQ(ierr); 147 ierr = VecDestroy(&upperb);CHKERRQ(ierr); 148 ierr = PetscFinalize(); 149 return ierr; 150 } 151 152 /* ------------------------------------------------------------------ */ 153 /* 154 FormFunction - Evaluates the function and corresponding gradient. 155 156 Input Parameters: 157 tao - the Tao context 158 X - the input vector 159 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine() 160 161 Output Parameters: 162 f - the newly evaluated function 163 */ 164 PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0) 165 { 166 AppCtx *ctx = (AppCtx*)ctx0; 167 TS ts,quadts; 168 Vec U; /* solution will be stored here */ 169 Mat A; /* Jacobian matrix */ 170 PetscErrorCode ierr; 171 PetscInt n = 2; 172 PetscReal ftime; 173 PetscInt steps; 174 PetscScalar *u; 175 const PetscScalar *x_ptr,*qx_ptr; 176 Vec q; 177 PetscInt direction[2]; 178 PetscBool terminate[2]; 179 180 ierr = VecGetArrayRead(P,&x_ptr);CHKERRQ(ierr); 181 ctx->Pm = x_ptr[0]; 182 ierr = VecRestoreArrayRead(P,&x_ptr);CHKERRQ(ierr); 183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 Create necessary matrix and vectors 185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 186 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 187 ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 188 ierr = MatSetType(A,MATDENSE);CHKERRQ(ierr); 189 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 190 ierr = MatSetUp(A);CHKERRQ(ierr); 191 192 ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Create timestepping solver context 196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 198 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 199 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 200 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,ctx);CHKERRQ(ierr); 201 ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,ctx);CHKERRQ(ierr); 202 203 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 204 Set initial conditions 205 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 206 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 207 u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax); 208 u[1] = 1.0; 209 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 210 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 211 212 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 213 Set solver options 214 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 215 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 216 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 217 ierr = TSSetTimeStep(ts,0.03125);CHKERRQ(ierr); 218 ierr = TSCreateQuadratureTS(ts,PETSC_TRUE,&quadts);CHKERRQ(ierr); 219 ierr = TSGetSolution(quadts,&q);CHKERRQ(ierr); 220 ierr = VecSet(q,0.0);CHKERRQ(ierr); 221 ierr = TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,ctx);CHKERRQ(ierr); 222 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 223 224 direction[0] = direction[1] = 1; 225 terminate[0] = terminate[1] = PETSC_FALSE; 226 227 ierr = TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)ctx);CHKERRQ(ierr); 228 229 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 230 Solve nonlinear system 231 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 232 ierr = TSSolve(ts,U);CHKERRQ(ierr); 233 234 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 235 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 236 ierr = VecGetArrayRead(q,&qx_ptr);CHKERRQ(ierr); 237 *f = -ctx->Pm + qx_ptr[0]; 238 ierr = VecRestoreArrayRead(q,&qx_ptr);CHKERRQ(ierr); 239 240 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 241 Free work space. All PETSc objects should be destroyed when they are no longer needed. 242 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 243 ierr = MatDestroy(&A);CHKERRQ(ierr); 244 ierr = VecDestroy(&U);CHKERRQ(ierr); 245 ierr = TSDestroy(&ts);CHKERRQ(ierr); 246 return 0; 247 } 248 249 /*TEST 250 251 build: 252 requires: !complex !single 253 254 test: 255 args: -ts_type cn -pc_type lu -tao_monitor -tao_gatol 1e-3 256 257 TEST*/ 258