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
monitor(Tao tao,AppCtx * ctx)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
main(int argc,char ** argv)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 */
FormFunction(Tao tao,Vec P,PetscReal * f,void * ctx0)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