Actual source code: ex3opt_fd.c
petsc-3.11.4 2019-09-28
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: 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>
20: typedef struct {
21: PetscScalar H,D,omega_b,omega_s,Pmax,Pmax_ini,Pm,E,V,X,u_s,c;
22: PetscInt beta;
23: PetscReal tf,tcl;
24: } AppCtx;
26: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
28: /* Event check */
29: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec X,PetscScalar *fvalue,void *ctx)
30: {
31: AppCtx *user=(AppCtx*)ctx;
34: /* Event for fault-on time */
35: fvalue[0] = t - user->tf;
36: /* Event for fault-off time */
37: fvalue[1] = t - user->tcl;
39: return(0);
40: }
42: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec X,PetscBool forwardsolve,void* ctx)
43: {
44: AppCtx *user=(AppCtx*)ctx;
48: if (event_list[0] == 0) {
49: if (forwardsolve) user->Pmax = 0.0; /* Apply disturbance - this is done by setting Pmax = 0 */
50: else user->Pmax = user->Pmax_ini; /* Going backward, reversal of event */
51: } else if(event_list[0] == 1) {
52: if (forwardsolve) user->Pmax = user->Pmax_ini; /* Remove the fault - this is done by setting Pmax = Pmax_ini */
53: else user->Pmax = 0.0; /* Going backward, reversal of event */
54: }
55: return(0);
56: }
58: /*
59: Defines the ODE passed to the ODE solver
60: */
61: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx)
62: {
63: PetscErrorCode ierr;
64: PetscScalar *f,Pmax;
65: const PetscScalar *u,*udot;
68: /* The next three lines allow us to access the entries of the vectors directly */
69: VecGetArrayRead(U,&u);
70: VecGetArrayRead(Udot,&udot);
71: VecGetArray(F,&f);
72: Pmax = ctx->Pmax;
74: f[0] = udot[0] - ctx->omega_b*(u[1] - ctx->omega_s);
75: f[1] = 2.0*ctx->H/ctx->omega_s*udot[1] + Pmax*PetscSinScalar(u[0]) + ctx->D*(u[1] - ctx->omega_s)- ctx->Pm;
77: VecRestoreArrayRead(U,&u);
78: VecRestoreArrayRead(Udot,&udot);
79: VecRestoreArray(F,&f);
80: return(0);
81: }
83: /*
84: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
85: */
86: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx)
87: {
88: PetscErrorCode ierr;
89: PetscInt rowcol[] = {0,1};
90: PetscScalar J[2][2],Pmax;
91: const PetscScalar *u,*udot;
94: VecGetArrayRead(U,&u);
95: VecGetArrayRead(Udot,&udot);
96: Pmax = ctx->Pmax;
98: J[0][0] = a; J[0][1] = -ctx->omega_b;
99: J[1][1] = 2.0*ctx->H/ctx->omega_s*a + ctx->D; J[1][0] = Pmax*PetscCosScalar(u[0]);
101: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
102: VecRestoreArrayRead(U,&u);
103: VecRestoreArrayRead(Udot,&udot);
105: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
106: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
107: if (A != B) {
108: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
109: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
110: }
111: return(0);
112: }
114: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
115: {
116: PetscErrorCode ierr;
117: PetscScalar *r;
118: const PetscScalar *u;
121: VecGetArrayRead(U,&u);
122: VecGetArray(R,&r);
123: r[0] = ctx->c*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta);
124: VecRestoreArray(R,&r);
125: VecRestoreArrayRead(U,&u);
126: return(0);
127: }
129: PetscErrorCode monitor(Tao tao,AppCtx *ctx)
130: {
131: FILE *fp;
132: PetscInt iterate;
133: PetscReal f,gnorm,cnorm,xdiff;
134: Vec X,G;
135: const PetscScalar *x,*g;
136: TaoConvergedReason reason;
137: PetscErrorCode ierr;
140: TaoGetSolutionStatus(tao,&iterate,&f,&gnorm,&cnorm,&xdiff,&reason);
141: TaoGetSolutionVector(tao,&X);
142: TaoGetGradientVector(tao,&G);
143: VecGetArrayRead(X,&x);
144: VecGetArrayRead(G,&g);
145: fp = fopen("ex3opt_fd_conv.out","a");
146: PetscFPrintf(PETSC_COMM_WORLD,fp,"%d %g %.12lf %.12lf\n",iterate,gnorm,x[0],g[0]);
147: VecRestoreArrayRead(X,&x);
148: VecRestoreArrayRead(G,&g);
149: fclose(fp);
150: return(0);
151: }
153: int main(int argc,char **argv)
154: {
155: Vec p;
156: PetscScalar *x_ptr;
157: PetscErrorCode ierr;
158: PetscMPIInt size;
159: AppCtx ctx;
160: Vec lowerb,upperb;
161: Tao tao;
162: KSP ksp;
163: PC pc;
164: PetscBool printtofile;
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Initialize program
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
170: MPI_Comm_size(PETSC_COMM_WORLD,&size);
171: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Set runtime options
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
177: {
178: ctx.beta = 2;
179: ctx.c = 10000.0;
180: ctx.u_s = 1.0;
181: ctx.omega_s = 1.0;
182: ctx.omega_b = 120.0*PETSC_PI;
183: ctx.H = 5.0;
184: PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
185: ctx.D = 5.0;
186: PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
187: ctx.E = 1.1378;
188: ctx.V = 1.0;
189: ctx.X = 0.545;
190: ctx.Pmax = ctx.E*ctx.V/ctx.X;
191: ctx.Pmax_ini = ctx.Pmax;
192: PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
193: ctx.Pm = 1.06;
194: PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
195: ctx.tf = 0.1;
196: ctx.tcl = 0.2;
197: PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
198: PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
199: printtofile = PETSC_FALSE;
200: PetscOptionsBool("-printtofile","Print convergence results to file","",printtofile,&printtofile,NULL);
201: }
202: PetscOptionsEnd();
204: /* Create TAO solver and set desired solution method */
205: TaoCreate(PETSC_COMM_WORLD,&tao);
206: TaoSetType(tao,TAOBLMVM);
207: if(printtofile) {
208: TaoSetMonitor(tao,(PetscErrorCode (*)(Tao, void*))monitor,(void *)&ctx,PETSC_NULL);
209: }
210: TaoSetMaximumIterations(tao,30);
211: /*
212: Optimization starts
213: */
214: /* Set initial solution guess */
215: VecCreateSeq(PETSC_COMM_WORLD,1,&p);
216: VecGetArray(p,&x_ptr);
217: x_ptr[0] = ctx.Pm;
218: VecRestoreArray(p,&x_ptr);
220: TaoSetInitialVector(tao,p);
221: /* Set routine for function and gradient evaluation */
222: TaoSetObjectiveRoutine(tao,FormFunction,(void *)&ctx);
223: TaoSetGradientRoutine(tao,TaoDefaultComputeGradient,(void *)&ctx);
225: /* Set bounds for the optimization */
226: VecDuplicate(p,&lowerb);
227: VecDuplicate(p,&upperb);
228: VecGetArray(lowerb,&x_ptr);
229: x_ptr[0] = 0.;
230: VecRestoreArray(lowerb,&x_ptr);
231: VecGetArray(upperb,&x_ptr);
232: x_ptr[0] = 1.1;
233: VecRestoreArray(upperb,&x_ptr);
234: TaoSetVariableBounds(tao,lowerb,upperb);
236: /* Check for any TAO command line options */
237: TaoSetFromOptions(tao);
238: TaoGetKSP(tao,&ksp);
239: if (ksp) {
240: KSPGetPC(ksp,&pc);
241: PCSetType(pc,PCNONE);
242: }
244: /* SOLVE THE APPLICATION */
245: TaoSolve(tao);
247: VecView(p,PETSC_VIEWER_STDOUT_WORLD);
249: /* Free TAO data structures */
250: TaoDestroy(&tao);
251: VecDestroy(&p);
252: VecDestroy(&lowerb);
253: VecDestroy(&upperb);
254: PetscFinalize();
255: return ierr;
256: }
258: /* ------------------------------------------------------------------ */
259: /*
260: FormFunction - Evaluates the function and corresponding gradient.
262: Input Parameters:
263: tao - the Tao context
264: X - the input vector
265: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
267: Output Parameters:
268: f - the newly evaluated function
269: */
270: PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0)
271: {
272: AppCtx *ctx = (AppCtx*)ctx0;
273: TS ts;
275: Vec U; /* solution will be stored here */
276: Mat A; /* Jacobian matrix */
277: PetscErrorCode ierr;
278: PetscInt n = 2;
279: PetscReal ftime;
280: PetscInt steps;
281: PetscScalar *u;
282: const PetscScalar *x_ptr,*qx_ptr;
283: Vec q;
284: PetscInt direction[2];
285: PetscBool terminate[2];
287: VecGetArrayRead(P,&x_ptr);
288: ctx->Pm = x_ptr[0];
289: VecRestoreArrayRead(P,&x_ptr);
290: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291: Create necessary matrix and vectors
292: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
293: MatCreate(PETSC_COMM_WORLD,&A);
294: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
295: MatSetType(A,MATDENSE);
296: MatSetFromOptions(A);
297: MatSetUp(A);
299: MatCreateVecs(A,&U,NULL);
301: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
302: Create timestepping solver context
303: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
304: TSCreate(PETSC_COMM_WORLD,&ts);
305: TSSetProblemType(ts,TS_NONLINEAR);
306: TSSetType(ts,TSCN);
307: TSSetIFunction(ts,NULL,(TSIFunction) IFunction,ctx);
308: TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,ctx);
310: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
311: Set initial conditions
312: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
313: VecGetArray(U,&u);
314: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
315: u[1] = 1.0;
316: VecRestoreArray(U,&u);
317: TSSetSolution(ts,U);
319: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
320: Set solver options
321: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
322: TSSetMaxTime(ts,1.0);
323: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
324: TSSetTimeStep(ts,0.03125);
325: TSSetCostIntegrand(ts,1,NULL,(PetscErrorCode (*)(TS,PetscReal,Vec,Vec,void*))CostIntegrand,PETSC_NULL,PETSC_NULL,PETSC_TRUE,ctx);
326: TSSetFromOptions(ts);
328: direction[0] = direction[1] = 1;
329: terminate[0] = terminate[1] = PETSC_FALSE;
331: TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)ctx);
333: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
334: Solve nonlinear system
335: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
336: TSSolve(ts,U);
338: TSGetSolveTime(ts,&ftime);
339: TSGetStepNumber(ts,&steps);
341: TSGetCostIntegral(ts,&q);
343: VecGetArrayRead(q,&qx_ptr);
344: *f = -ctx->Pm + qx_ptr[0];
345: VecRestoreArrayRead(q,&qx_ptr);
347: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: Free work space. All PETSc objects should be destroyed when they are no longer needed.
349: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350: MatDestroy(&A);
351: VecDestroy(&U);
352: TSDestroy(&ts);
354: return 0;
355: }
357: /*TEST
359: build:
360: requires: !complex !single
362: test:
363: args: -ts_type cn -pc_type lu -tao_monitor -tao_gatol 1e-3
365: TEST*/