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