Actual source code: ex3opt_fd.c
petsc-3.12.2 2019-11-22
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>
19: #include "ex3.h"
21: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
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;
34: TaoGetSolutionStatus(tao,&iterate,&f,&gnorm,&cnorm,&xdiff,&reason);
35: TaoGetSolutionVector(tao,&X);
36: TaoGetGradientVector(tao,&G);
37: VecGetArrayRead(X,&x);
38: VecGetArrayRead(G,&g);
39: fp = fopen("ex3opt_fd_conv.out","a");
40: PetscFPrintf(PETSC_COMM_WORLD,fp,"%d %g %.12lf %.12lf\n",iterate,gnorm,x[0],g[0]);
41: VecRestoreArrayRead(X,&x);
42: VecRestoreArrayRead(G,&g);
43: fclose(fp);
44: return(0);
45: }
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: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
64: MPI_Comm_size(PETSC_COMM_WORLD,&size);
65: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Set runtime options
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
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: PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
79: ctx.D = 5.0;
80: PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
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: PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
87: ctx.Pm = 1.06;
88: PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
89: ctx.tf = 0.1;
90: ctx.tcl = 0.2;
91: PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
92: PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
93: printtofile = PETSC_FALSE;
94: PetscOptionsBool("-printtofile","Print convergence results to file","",printtofile,&printtofile,NULL);
95: }
96: PetscOptionsEnd();
98: /* Create TAO solver and set desired solution method */
99: TaoCreate(PETSC_COMM_WORLD,&tao);
100: TaoSetType(tao,TAOBLMVM);
101: if(printtofile) {
102: TaoSetMonitor(tao,(PetscErrorCode (*)(Tao, void*))monitor,(void *)&ctx,PETSC_NULL);
103: }
104: TaoSetMaximumIterations(tao,30);
105: /*
106: Optimization starts
107: */
108: /* Set initial solution guess */
109: VecCreateSeq(PETSC_COMM_WORLD,1,&p);
110: VecGetArray(p,&x_ptr);
111: x_ptr[0] = ctx.Pm;
112: VecRestoreArray(p,&x_ptr);
114: TaoSetInitialVector(tao,p);
115: /* Set routine for function and gradient evaluation */
116: TaoSetObjectiveRoutine(tao,FormFunction,(void *)&ctx);
117: TaoSetGradientRoutine(tao,TaoDefaultComputeGradient,(void *)&ctx);
119: /* Set bounds for the optimization */
120: VecDuplicate(p,&lowerb);
121: VecDuplicate(p,&upperb);
122: VecGetArray(lowerb,&x_ptr);
123: x_ptr[0] = 0.;
124: VecRestoreArray(lowerb,&x_ptr);
125: VecGetArray(upperb,&x_ptr);
126: x_ptr[0] = 1.1;
127: VecRestoreArray(upperb,&x_ptr);
128: TaoSetVariableBounds(tao,lowerb,upperb);
130: /* Check for any TAO command line options */
131: TaoSetFromOptions(tao);
132: TaoGetKSP(tao,&ksp);
133: if (ksp) {
134: KSPGetPC(ksp,&pc);
135: PCSetType(pc,PCNONE);
136: }
138: /* SOLVE THE APPLICATION */
139: TaoSolve(tao);
141: VecView(p,PETSC_VIEWER_STDOUT_WORLD);
143: /* Free TAO data structures */
144: TaoDestroy(&tao);
145: VecDestroy(&p);
146: VecDestroy(&lowerb);
147: VecDestroy(&upperb);
148: PetscFinalize();
149: return ierr;
150: }
152: /* ------------------------------------------------------------------ */
153: /*
154: FormFunction - Evaluates the function and corresponding gradient.
156: Input Parameters:
157: tao - the Tao context
158: X - the input vector
159: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
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];
180: VecGetArrayRead(P,&x_ptr);
181: ctx->Pm = x_ptr[0];
182: VecRestoreArrayRead(P,&x_ptr);
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Create necessary matrix and vectors
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: MatCreate(PETSC_COMM_WORLD,&A);
187: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
188: MatSetType(A,MATDENSE);
189: MatSetFromOptions(A);
190: MatSetUp(A);
192: MatCreateVecs(A,&U,NULL);
194: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: Create timestepping solver context
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197: TSCreate(PETSC_COMM_WORLD,&ts);
198: TSSetProblemType(ts,TS_NONLINEAR);
199: TSSetType(ts,TSCN);
200: TSSetIFunction(ts,NULL,(TSIFunction) IFunction,ctx);
201: TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,ctx);
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Set initial conditions
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: VecGetArray(U,&u);
207: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
208: u[1] = 1.0;
209: VecRestoreArray(U,&u);
210: TSSetSolution(ts,U);
212: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213: Set solver options
214: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
215: TSSetMaxTime(ts,1.0);
216: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
217: TSSetTimeStep(ts,0.03125);
218: TSCreateQuadratureTS(ts,PETSC_TRUE,&quadts);
219: TSGetSolution(quadts,&q);
220: VecSet(q,0.0);
221: TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,ctx);
222: TSSetFromOptions(ts);
224: direction[0] = direction[1] = 1;
225: terminate[0] = terminate[1] = PETSC_FALSE;
227: TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)ctx);
229: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230: Solve nonlinear system
231: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
232: TSSolve(ts,U);
234: TSGetSolveTime(ts,&ftime);
235: TSGetStepNumber(ts,&steps);
236: VecGetArrayRead(q,&qx_ptr);
237: *f = -ctx->Pm + qx_ptr[0];
238: VecRestoreArrayRead(q,&qx_ptr);
240: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241: Free work space. All PETSc objects should be destroyed when they are no longer needed.
242: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
243: MatDestroy(&A);
244: VecDestroy(&U);
245: TSDestroy(&ts);
246: return 0;
247: }
249: /*TEST
251: build:
252: requires: !complex !single
254: test:
255: args: -ts_type cn -pc_type lu -tao_monitor -tao_gatol 1e-3
257: TEST*/