Actual source code: ex5.c


  2: static char help[] = "Basic equation for an induction generator driven by a wind turbine.\n";

  4: /*F
  5: \begin{eqnarray}
  6:           T_w\frac{dv_w}{dt} & = & v_w - v_we \\
  7:           2(H_t+H_m)\frac{ds}{dt} & = & P_w - P_e
  8: \end{eqnarray}
  9: F*/
 10: /*
 11:  - Pw is the power extracted from the wind turbine given by
 12:            Pw = 0.5*\rho*cp*Ar*vw^3

 14:  - The wind speed time series is modeled using a Weibull distribution and then
 15:    passed through a low pass filter (with time constant T_w).
 16:  - v_we is the wind speed data calculated using Weibull distribution while v_w is
 17:    the output of the filter.
 18:  - P_e is assumed as constant electrical torque

 20:  - This example does not work with adaptive time stepping!

 22: Reference:
 23: Power System Modeling and Scripting - F. Milano
 24: */

 26: #include <petscts.h>

 28: #define freq    50
 29: #define ws      (2 * PETSC_PI * freq)
 30: #define MVAbase 100

 32: typedef struct {
 33:   /* Parameters for wind speed model */
 34:   PetscInt  nsamples;  /* Number of wind samples */
 35:   PetscReal cw;        /* Scale factor for Weibull distribution */
 36:   PetscReal kw;        /* Shape factor for Weibull distribution */
 37:   Vec       wind_data; /* Vector to hold wind speeds */
 38:   Vec       t_wind;    /* Vector to hold wind speed times */
 39:   PetscReal Tw;        /* Filter time constant */

 41:   /* Wind turbine parameters */
 42:   PetscScalar Rt;  /* Rotor radius */
 43:   PetscScalar Ar;  /* Area swept by rotor (pi*R*R) */
 44:   PetscReal   nGB; /* Gear box ratio */
 45:   PetscReal   Ht;  /* Turbine inertia constant */
 46:   PetscReal   rho; /* Atmospheric pressure */

 48:   /* Induction generator parameters */
 49:   PetscInt    np; /* Number of poles */
 50:   PetscReal   Xm; /* Magnetizing reactance */
 51:   PetscReal   Xs; /* Stator Reactance */
 52:   PetscReal   Xr; /* Rotor reactance */
 53:   PetscReal   Rs; /* Stator resistance */
 54:   PetscReal   Rr; /* Rotor resistance */
 55:   PetscReal   Hm; /* Motor inertia constant */
 56:   PetscReal   Xp; /* Xs + Xm*Xr/(Xm + Xr) */
 57:   PetscScalar Te; /* Electrical Torque */

 59:   Mat      Sol;     /* Solution matrix */
 60:   PetscInt stepnum; /* Column number of solution matrix */
 61: } AppCtx;

 63: /* Initial values computed by Power flow and initialization */
 64: PetscScalar s = -0.00011577790353;
 65: /*Pw = 0.011064344110238; %Te*wm */
 66: PetscScalar vwa  = 22.317142184449754;
 67: PetscReal   tmax = 20.0;

 69: /* Saves the solution at each time to a matrix */
 70: PetscErrorCode SaveSolution(TS ts)
 71: {
 72:   AppCtx            *user;
 73:   Vec                X;
 74:   PetscScalar       *mat;
 75:   const PetscScalar *x;
 76:   PetscInt           idx;
 77:   PetscReal          t;

 79:   TSGetApplicationContext(ts, &user);
 80:   TSGetTime(ts, &t);
 81:   TSGetSolution(ts, &X);
 82:   idx = 3 * user->stepnum;
 83:   MatDenseGetArray(user->Sol, &mat);
 84:   VecGetArrayRead(X, &x);
 85:   mat[idx] = t;
 86:   PetscArraycpy(mat + idx + 1, x, 2);
 87:   MatDenseRestoreArray(user->Sol, &mat);
 88:   VecRestoreArrayRead(X, &x);
 89:   user->stepnum++;
 90:   return 0;
 91: }

 93: /* Computes the wind speed using Weibull distribution */
 94: PetscErrorCode WindSpeeds(AppCtx *user)
 95: {
 96:   PetscScalar *x, *t, avg_dev, sum;
 97:   PetscInt     i;

 99:   user->cw       = 5;
100:   user->kw       = 2; /* Rayleigh distribution */
101:   user->nsamples = 2000;
102:   user->Tw       = 0.2;
103:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Wind Speed Options", "");
104:   {
105:     PetscOptionsReal("-cw", "", "", user->cw, &user->cw, NULL);
106:     PetscOptionsReal("-kw", "", "", user->kw, &user->kw, NULL);
107:     PetscOptionsInt("-nsamples", "", "", user->nsamples, &user->nsamples, NULL);
108:     PetscOptionsReal("-Tw", "", "", user->Tw, &user->Tw, NULL);
109:   }
110:   PetscOptionsEnd();
111:   VecCreate(PETSC_COMM_WORLD, &user->wind_data);
112:   VecSetSizes(user->wind_data, PETSC_DECIDE, user->nsamples);
113:   VecSetFromOptions(user->wind_data);
114:   VecDuplicate(user->wind_data, &user->t_wind);

116:   VecGetArray(user->t_wind, &t);
117:   for (i = 0; i < user->nsamples; i++) t[i] = (i + 1) * tmax / user->nsamples;
118:   VecRestoreArray(user->t_wind, &t);

120:   /* Wind speed deviation = (-log(rand)/cw)^(1/kw) */
121:   VecSetRandom(user->wind_data, NULL);
122:   VecLog(user->wind_data);
123:   VecScale(user->wind_data, -1 / user->cw);
124:   VecGetArray(user->wind_data, &x);
125:   for (i = 0; i < user->nsamples; i++) x[i] = PetscPowScalar(x[i], (1 / user->kw));
126:   VecRestoreArray(user->wind_data, &x);
127:   VecSum(user->wind_data, &sum);
128:   avg_dev = sum / user->nsamples;
129:   /* Wind speed (t) = (1 + wind speed deviation(t) - avg_dev)*average wind speed */
130:   VecShift(user->wind_data, (1 - avg_dev));
131:   VecScale(user->wind_data, vwa);
132:   return 0;
133: }

135: /* Sets the parameters for wind turbine */
136: PetscErrorCode SetWindTurbineParams(AppCtx *user)
137: {
138:   user->Rt  = 35;
139:   user->Ar  = PETSC_PI * user->Rt * user->Rt;
140:   user->nGB = 1.0 / 89.0;
141:   user->rho = 1.225;
142:   user->Ht  = 1.5;
143:   return 0;
144: }

146: /* Sets the parameters for induction generator */
147: PetscErrorCode SetInductionGeneratorParams(AppCtx *user)
148: {
149:   user->np = 4;
150:   user->Xm = 3.0;
151:   user->Xs = 0.1;
152:   user->Xr = 0.08;
153:   user->Rs = 0.01;
154:   user->Rr = 0.01;
155:   user->Xp = user->Xs + user->Xm * user->Xr / (user->Xm + user->Xr);
156:   user->Hm = 1.0;
157:   user->Te = 0.011063063063251968;
158:   return 0;
159: }

161: /* Computes the power extracted from wind */
162: PetscErrorCode GetWindPower(PetscScalar wm, PetscScalar vw, PetscScalar *Pw, AppCtx *user)
163: {
164:   PetscScalar temp, lambda, lambda_i, cp;

166:   temp     = user->nGB * 2 * user->Rt * ws / user->np;
167:   lambda   = temp * wm / vw;
168:   lambda_i = 1 / (1 / lambda + 0.002);
169:   cp       = 0.44 * (125 / lambda_i - 6.94) * PetscExpScalar(-16.5 / lambda_i);
170:   *Pw      = 0.5 * user->rho * cp * user->Ar * vw * vw * vw / (MVAbase * 1e6);
171:   return 0;
172: }

174: /*
175:      Defines the ODE passed to the ODE solver
176: */
177: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *user)
178: {
179:   PetscScalar       *f, wm, Pw, *wd;
180:   const PetscScalar *u, *udot;
181:   PetscInt           stepnum;

183:   TSGetStepNumber(ts, &stepnum);
184:   /*  The next three lines allow us to access the entries of the vectors directly */
185:   VecGetArrayRead(U, &u);
186:   VecGetArrayRead(Udot, &udot);
187:   VecGetArray(F, &f);
188:   VecGetArray(user->wind_data, &wd);

190:   f[0] = user->Tw * udot[0] - wd[stepnum] + u[0];
191:   wm   = 1 - u[1];
192:   GetWindPower(wm, u[0], &Pw, user);
193:   f[1] = 2.0 * (user->Ht + user->Hm) * udot[1] - Pw / wm + user->Te;

195:   VecRestoreArray(user->wind_data, &wd);
196:   VecRestoreArrayRead(U, &u);
197:   VecRestoreArrayRead(Udot, &udot);
198:   VecRestoreArray(F, &f);
199:   return 0;
200: }

202: int main(int argc, char **argv)
203: {
204:   TS                 ts; /* ODE integrator */
205:   Vec                U;  /* solution will be stored here */
206:   Mat                A;  /* Jacobian matrix */
207:   PetscMPIInt        size;
208:   PetscInt           n = 2, idx;
209:   AppCtx             user;
210:   PetscScalar       *u;
211:   SNES               snes;
212:   PetscScalar       *mat;
213:   const PetscScalar *x, *rmat;
214:   Mat                B;
215:   PetscScalar       *amat;
216:   PetscViewer        viewer;

218:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219:      Initialize program
220:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222:   PetscInitialize(&argc, &argv, (char *)0, help);
223:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

226:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227:     Create necessary matrix and vectors
228:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
229:   MatCreate(PETSC_COMM_WORLD, &A);
230:   MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
231:   MatSetFromOptions(A);
232:   MatSetUp(A);

234:   MatCreateVecs(A, &U, NULL);

236:   /* Create wind speed data using Weibull distribution */
237:   WindSpeeds(&user);
238:   /* Set parameters for wind turbine and induction generator */
239:   SetWindTurbineParams(&user);
240:   SetInductionGeneratorParams(&user);

242:   VecGetArray(U, &u);
243:   u[0] = vwa;
244:   u[1] = s;
245:   VecRestoreArray(U, &u);

247:   /* Create matrix to save solutions at each time step */
248:   user.stepnum = 0;

250:   MatCreateSeqDense(PETSC_COMM_SELF, 3, 2010, NULL, &user.Sol);

252:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253:      Create timestepping solver context
254:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255:   TSCreate(PETSC_COMM_WORLD, &ts);
256:   TSSetProblemType(ts, TS_NONLINEAR);
257:   TSSetType(ts, TSBEULER);
258:   TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &user);

260:   TSGetSNES(ts, &snes);
261:   SNESSetJacobian(snes, A, A, SNESComputeJacobianDefault, NULL);
262:   /*  TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&user); */
263:   TSSetApplicationContext(ts, &user);

265:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266:      Set initial conditions
267:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
268:   TSSetSolution(ts, U);

270:   /* Save initial solution */
271:   idx = 3 * user.stepnum;

273:   MatDenseGetArray(user.Sol, &mat);
274:   VecGetArrayRead(U, &x);

276:   mat[idx] = 0.0;

278:   PetscArraycpy(mat + idx + 1, x, 2);
279:   MatDenseRestoreArray(user.Sol, &mat);
280:   VecRestoreArrayRead(U, &x);
281:   user.stepnum++;

283:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
284:      Set solver options
285:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
286:   TSSetMaxTime(ts, 20.0);
287:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
288:   TSSetTimeStep(ts, .01);
289:   TSSetFromOptions(ts);
290:   TSSetPostStep(ts, SaveSolution);
291:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
292:      Solve nonlinear system
293:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
294:   TSSolve(ts, U);

296:   MatCreateSeqDense(PETSC_COMM_SELF, 3, user.stepnum, NULL, &B);
297:   MatDenseGetArrayRead(user.Sol, &rmat);
298:   MatDenseGetArray(B, &amat);
299:   PetscArraycpy(amat, rmat, user.stepnum * 3);
300:   MatDenseRestoreArray(B, &amat);
301:   MatDenseRestoreArrayRead(user.Sol, &rmat);

303:   PetscViewerBinaryOpen(PETSC_COMM_SELF, "out.bin", FILE_MODE_WRITE, &viewer);
304:   MatView(B, viewer);
305:   PetscViewerDestroy(&viewer);
306:   MatDestroy(&user.Sol);
307:   MatDestroy(&B);
308:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
309:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
310:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
311:   VecDestroy(&user.wind_data);
312:   VecDestroy(&user.t_wind);
313:   MatDestroy(&A);
314:   VecDestroy(&U);
315:   TSDestroy(&ts);

317:   PetscFinalize();
318:   return 0;
319: }

321: /*TEST

323:    build:
324:       requires: !complex

326:    test:

328: TEST*/