Actual source code: ex9opt.c
2: static char help[] = "Basic equation for generator stability analysis.\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}
Ensemble of initial conditions
./ex2 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
Fault at .1 seconds
./ex2 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
Initial conditions same as when fault is ended
./ex2 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
25: /*
26: Include "petscts.h" so that we can use TS solvers. Note that this
27: file automatically includes:
28: petscsys.h - base PETSc routines petscvec.h - vectors
29: petscmat.h - matrices
30: petscis.h - index sets petscksp.h - Krylov subspace methods
31: petscviewer.h - viewers petscpc.h - preconditioners
32: petscksp.h - linear solvers
33: */
35: #include <petsctao.h>
36: #include <petscts.h>
38: typedef struct {
39: TS ts;
40: PetscScalar H,D,omega_b,omega_s,Pmax,Pm,E,V,X,u_s,c;
41: PetscInt beta;
42: PetscReal tf,tcl,dt;
43: } AppCtx;
45: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
46: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
48: /*
49: Defines the ODE passed to the ODE solver
50: */
51: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
52: {
53: PetscErrorCode ierr;
54: PetscScalar *f,Pmax;
55: const PetscScalar *u;
58: /* The next three lines allow us to access the entries of the vectors directly */
59: VecGetArrayRead(U,&u);
60: VecGetArray(F,&f);
61: if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
62: else Pmax = ctx->Pmax;
64: f[0] = ctx->omega_b*(u[1] - ctx->omega_s);
65: f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H);
67: VecRestoreArrayRead(U,&u);
68: VecRestoreArray(F,&f);
69: return(0);
70: }
72: /*
73: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
74: */
75: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx)
76: {
77: PetscErrorCode ierr;
78: PetscInt rowcol[] = {0,1};
79: PetscScalar J[2][2],Pmax;
80: const PetscScalar *u;
83: VecGetArrayRead(U,&u);
84: if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
85: else Pmax = ctx->Pmax;
87: J[0][0] = 0; J[0][1] = ctx->omega_b;
88: J[1][1] = -ctx->D*ctx->omega_s/(2.0*ctx->H); J[1][0] = -Pmax*PetscCosScalar(u[0])*ctx->omega_s/(2.0*ctx->H);
90: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
91: VecRestoreArrayRead(U,&u);
93: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
95: if (A != B) {
96: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
98: }
99: return(0);
100: }
102: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
103: {
105: PetscInt row[] = {0,1},col[]={0};
106: PetscScalar J[2][1];
107: AppCtx *ctx=(AppCtx*)ctx0;
110: J[0][0] = 0;
111: J[1][0] = ctx->omega_s/(2.0*ctx->H);
112: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
113: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
114: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
115: return(0);
116: }
118: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
119: {
120: PetscErrorCode ierr;
121: PetscScalar *r;
122: const PetscScalar *u;
125: VecGetArrayRead(U,&u);
126: VecGetArray(R,&r);
127: r[0] = ctx->c*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta);
128: VecRestoreArray(R,&r);
129: VecRestoreArrayRead(U,&u);
130: return(0);
131: }
133: static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx *ctx)
134: {
135: PetscErrorCode ierr;
136: PetscScalar ru[1];
137: const PetscScalar *u;
138: PetscInt row[] = {0},col[] = {0};
141: VecGetArrayRead(U,&u);
142: ru[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta-1);
143: VecRestoreArrayRead(U,&u);
144: MatSetValues(DRDU,1,row,1,col,ru,INSERT_VALUES);
145: MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);
146: MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);
147: return(0);
148: }
150: static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx *ctx)
151: {
155: MatZeroEntries(DRDP);
156: MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);
157: MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);
158: return(0);
159: }
161: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
162: {
163: PetscErrorCode ierr;
164: PetscScalar *y,sensip;
165: const PetscScalar *x;
168: VecGetArrayRead(lambda,&x);
169: VecGetArray(mu,&y);
170: sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
171: y[0] = sensip;
172: VecRestoreArray(mu,&y);
173: VecRestoreArrayRead(lambda,&x);
174: return(0);
175: }
177: int main(int argc,char **argv)
178: {
179: Vec p;
180: PetscScalar *x_ptr;
182: PetscMPIInt size;
183: AppCtx ctx;
184: Vec lowerb,upperb;
185: Tao tao;
186: KSP ksp;
187: PC pc;
188: Vec U,lambda[1],mu[1];
189: Mat A; /* Jacobian matrix */
190: Mat Jacp; /* Jacobian matrix */
191: Mat DRDU,DRDP;
192: PetscInt n = 2;
193: TS quadts;
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Initialize program
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
200: MPI_Comm_size(PETSC_COMM_WORLD,&size);
201: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Set runtime options
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
207: {
208: ctx.beta = 2;
209: ctx.c = PetscRealConstant(10000.0);
210: ctx.u_s = PetscRealConstant(1.0);
211: ctx.omega_s = PetscRealConstant(1.0);
212: ctx.omega_b = PetscRealConstant(120.0)*PETSC_PI;
213: ctx.H = PetscRealConstant(5.0);
214: PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
215: ctx.D = PetscRealConstant(5.0);
216: PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
217: ctx.E = PetscRealConstant(1.1378);
218: ctx.V = PetscRealConstant(1.0);
219: ctx.X = PetscRealConstant(0.545);
220: ctx.Pmax = ctx.E*ctx.V/ctx.X;
221: PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
222: ctx.Pm = PetscRealConstant(1.0194);
223: PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
224: ctx.tf = PetscRealConstant(0.1);
225: ctx.tcl = PetscRealConstant(0.2);
226: PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
227: PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
229: }
230: PetscOptionsEnd();
232: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233: Create necessary matrix and vectors
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: MatCreate(PETSC_COMM_WORLD,&A);
236: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
237: MatSetType(A,MATDENSE);
238: MatSetFromOptions(A);
239: MatSetUp(A);
241: MatCreateVecs(A,&U,NULL);
243: MatCreate(PETSC_COMM_WORLD,&Jacp);
244: MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
245: MatSetFromOptions(Jacp);
246: MatSetUp(Jacp);
247: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&DRDP);
248: MatSetUp(DRDP);
249: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,2,NULL,&DRDU);
250: MatSetUp(DRDU);
252: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253: Create timestepping solver context
254: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255: TSCreate(PETSC_COMM_WORLD,&ctx.ts);
256: TSSetProblemType(ctx.ts,TS_NONLINEAR);
257: TSSetEquationType(ctx.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
258: TSSetType(ctx.ts,TSRK);
259: TSSetRHSFunction(ctx.ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
260: TSSetRHSJacobian(ctx.ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);
261: TSSetExactFinalTime(ctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
263: MatCreateVecs(A,&lambda[0],NULL);
264: MatCreateVecs(Jacp,&mu[0],NULL);
265: TSSetCostGradients(ctx.ts,1,lambda,mu);
266: TSSetRHSJacobianP(ctx.ts,Jacp,RHSJacobianP,&ctx);
268: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: Set solver options
270: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: TSSetMaxTime(ctx.ts,PetscRealConstant(1.0));
272: TSSetTimeStep(ctx.ts,PetscRealConstant(0.01));
273: TSSetFromOptions(ctx.ts);
275: TSGetTimeStep(ctx.ts,&ctx.dt); /* save the stepsize */
277: TSCreateQuadratureTS(ctx.ts,PETSC_TRUE,&quadts);
278: TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);
279: TSSetRHSJacobian(quadts,DRDU,DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);
280: TSSetRHSJacobianP(quadts,DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&ctx);
281: TSSetSolution(ctx.ts,U);
283: /* Create TAO solver and set desired solution method */
284: TaoCreate(PETSC_COMM_WORLD,&tao);
285: TaoSetType(tao,TAOBLMVM);
287: /*
288: Optimization starts
289: */
290: /* Set initial solution guess */
291: VecCreateSeq(PETSC_COMM_WORLD,1,&p);
292: VecGetArray(p,&x_ptr);
293: x_ptr[0] = ctx.Pm;
294: VecRestoreArray(p,&x_ptr);
296: TaoSetInitialVector(tao,p);
297: /* Set routine for function and gradient evaluation */
298: TaoSetObjectiveRoutine(tao,FormFunction,(void *)&ctx);
299: TaoSetGradientRoutine(tao,FormGradient,(void *)&ctx);
301: /* Set bounds for the optimization */
302: VecDuplicate(p,&lowerb);
303: VecDuplicate(p,&upperb);
304: VecGetArray(lowerb,&x_ptr);
305: x_ptr[0] = 0.;
306: VecRestoreArray(lowerb,&x_ptr);
307: VecGetArray(upperb,&x_ptr);
308: x_ptr[0] = PetscRealConstant(1.1);
309: VecRestoreArray(upperb,&x_ptr);
310: TaoSetVariableBounds(tao,lowerb,upperb);
312: /* Check for any TAO command line options */
313: TaoSetFromOptions(tao);
314: TaoGetKSP(tao,&ksp);
315: if (ksp) {
316: KSPGetPC(ksp,&pc);
317: PCSetType(pc,PCNONE);
318: }
320: /* SOLVE THE APPLICATION */
321: TaoSolve(tao);
323: VecView(p,PETSC_VIEWER_STDOUT_WORLD);
324: /* Free TAO data structures */
325: TaoDestroy(&tao);
326: VecDestroy(&p);
327: VecDestroy(&lowerb);
328: VecDestroy(&upperb);
330: TSDestroy(&ctx.ts);
331: VecDestroy(&U);
332: MatDestroy(&A);
333: MatDestroy(&Jacp);
334: MatDestroy(&DRDU);
335: MatDestroy(&DRDP);
336: VecDestroy(&lambda[0]);
337: VecDestroy(&mu[0]);
338: PetscFinalize();
339: return ierr;
340: }
342: /* ------------------------------------------------------------------ */
343: /*
344: FormFunction - Evaluates the function
346: Input Parameters:
347: tao - the Tao context
348: X - the input vector
349: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
351: Output Parameters:
352: f - the newly evaluated function
353: */
354: PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0)
355: {
356: AppCtx *ctx = (AppCtx*)ctx0;
357: TS ts = ctx->ts;
358: Vec U; /* solution will be stored here */
360: PetscScalar *u;
361: PetscScalar *x_ptr;
362: Vec q;
364: VecGetArrayRead(P,(const PetscScalar**)&x_ptr);
365: ctx->Pm = x_ptr[0];
366: VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);
368: /* reset time */
369: TSSetTime(ts,0.0);
370: /* reset step counter, this is critical for adjoint solver */
371: TSSetStepNumber(ts,0);
372: /* reset step size, the step size becomes negative after TSAdjointSolve */
373: TSSetTimeStep(ts,ctx->dt);
374: /* reinitialize the integral value */
375: TSGetCostIntegral(ts,&q);
376: VecSet(q,0.0);
378: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
379: Set initial conditions
380: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
381: TSGetSolution(ts,&U);
382: VecGetArray(U,&u);
383: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
384: u[1] = PetscRealConstant(1.0);
385: VecRestoreArray(U,&u);
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: Solve nonlinear system
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390: TSSolve(ts,U);
391: TSGetCostIntegral(ts,&q);
392: VecGetArray(q,&x_ptr);
393: *f = -ctx->Pm + x_ptr[0];
394: VecRestoreArray(q,&x_ptr);
395: return 0;
396: }
398: PetscErrorCode FormGradient(Tao tao,Vec P,Vec G,void *ctx0)
399: {
400: AppCtx *ctx = (AppCtx*)ctx0;
401: TS ts = ctx->ts;
402: Vec U; /* solution will be stored here */
404: PetscReal ftime;
405: PetscInt steps;
406: PetscScalar *u;
407: PetscScalar *x_ptr,*y_ptr;
408: Vec *lambda,q,*mu;
410: VecGetArrayRead(P,(const PetscScalar**)&x_ptr);
411: ctx->Pm = x_ptr[0];
412: VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr);
414: /* reset time */
415: TSSetTime(ts,0.0);
416: /* reset step counter, this is critical for adjoint solver */
417: TSSetStepNumber(ts,0);
418: /* reset step size, the step size becomes negative after TSAdjointSolve */
419: TSSetTimeStep(ts,ctx->dt);
420: /* reinitialize the integral value */
421: TSGetCostIntegral(ts,&q);
422: VecSet(q,0.0);
424: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
425: Set initial conditions
426: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
427: TSGetSolution(ts,&U);
428: VecGetArray(U,&u);
429: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
430: u[1] = PetscRealConstant(1.0);
431: VecRestoreArray(U,&u);
433: /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
434: TSSetSaveTrajectory(ts);
435: TSSetFromOptions(ts);
437: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
438: Solve nonlinear system
439: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
440: TSSolve(ts,U);
442: TSGetSolveTime(ts,&ftime);
443: TSGetStepNumber(ts,&steps);
445: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
446: Adjoint model starts here
447: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
448: TSGetCostGradients(ts,NULL,&lambda,&mu);
449: /* Set initial conditions for the adjoint integration */
450: VecGetArray(lambda[0],&y_ptr);
451: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
452: VecRestoreArray(lambda[0],&y_ptr);
453: VecGetArray(mu[0],&x_ptr);
454: x_ptr[0] = PetscRealConstant(-1.0);
455: VecRestoreArray(mu[0],&x_ptr);
457: TSAdjointSolve(ts);
458: TSGetCostIntegral(ts,&q);
459: ComputeSensiP(lambda[0],mu[0],ctx);
460: VecCopy(mu[0],G);
461: return 0;
462: }
465: /*TEST
467: build:
468: requires: !complex
470: test:
471: args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason
473: test:
474: suffix: 2
475: args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason -tao_test_gradient
477: TEST*/