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*/