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