1 static char help[] = "Basic equation for an induction generator driven by a wind turbine.\n";
2
3 /*F
4 \begin{eqnarray}
5 T_w\frac{dv_w}{dt} & = & v_w - v_we \\
6 2(H_t+H_m)\frac{ds}{dt} & = & P_w - P_e
7 \end{eqnarray}
8 F*/
9 /*
10 - Pw is the power extracted from the wind turbine given by
11 Pw = 0.5*\rho*cp*Ar*vw^3
12
13 - The wind speed time series is modeled using a Weibull distribution and then
14 passed through a low pass filter (with time constant T_w).
15 - v_we is the wind speed data calculated using Weibull distribution while v_w is
16 the output of the filter.
17 - P_e is assumed as constant electrical torque
18
19 - This example does not work with adaptive time stepping!
20
21 Reference:
22 Power System Modeling and Scripting - F. Milano
23 */
24
25 #include <petscts.h>
26
27 #define freq 50
28 #define ws (2 * PETSC_PI * freq)
29 #define MVAbase 100
30
31 typedef struct {
32 /* Parameters for wind speed model */
33 PetscInt nsamples; /* Number of wind samples */
34 PetscReal cw; /* Scale factor for Weibull distribution */
35 PetscReal kw; /* Shape factor for Weibull distribution */
36 Vec wind_data; /* Vector to hold wind speeds */
37 Vec t_wind; /* Vector to hold wind speed times */
38 PetscReal Tw; /* Filter time constant */
39
40 /* Wind turbine parameters */
41 PetscScalar Rt; /* Rotor radius */
42 PetscScalar Ar; /* Area swept by rotor (pi*R*R) */
43 PetscReal nGB; /* Gear box ratio */
44 PetscReal Ht; /* Turbine inertia constant */
45 PetscReal rho; /* Atmospheric pressure */
46
47 /* Induction generator parameters */
48 PetscInt np; /* Number of poles */
49 PetscReal Xm; /* Magnetizing reactance */
50 PetscReal Xs; /* Stator Reactance */
51 PetscReal Xr; /* Rotor reactance */
52 PetscReal Rs; /* Stator resistance */
53 PetscReal Rr; /* Rotor resistance */
54 PetscReal Hm; /* Motor inertia constant */
55 PetscReal Xp; /* Xs + Xm*Xr/(Xm + Xr) */
56 PetscScalar Te; /* Electrical Torque */
57
58 Mat Sol; /* Solution matrix */
59 PetscInt stepnum; /* Column number of solution matrix */
60 } AppCtx;
61
62 /* Initial values computed by Power flow and initialization */
63 PetscScalar s = -0.00011577790353;
64 /*Pw = 0.011064344110238; %Te*wm */
65 PetscScalar vwa = 22.317142184449754;
66 PetscReal tmax = 20.0;
67
68 /* Saves the solution at each time to a matrix */
SaveSolution(TS ts)69 PetscErrorCode SaveSolution(TS ts)
70 {
71 AppCtx *user;
72 Vec X;
73 PetscScalar *mat;
74 const PetscScalar *x;
75 PetscInt idx;
76 PetscReal t;
77
78 PetscFunctionBegin;
79 PetscCall(TSGetApplicationContext(ts, &user));
80 PetscCall(TSGetTime(ts, &t));
81 PetscCall(TSGetSolution(ts, &X));
82 idx = 3 * user->stepnum;
83 PetscCall(MatDenseGetArray(user->Sol, &mat));
84 PetscCall(VecGetArrayRead(X, &x));
85 mat[idx] = t;
86 PetscCall(PetscArraycpy(mat + idx + 1, x, 2));
87 PetscCall(MatDenseRestoreArray(user->Sol, &mat));
88 PetscCall(VecRestoreArrayRead(X, &x));
89 user->stepnum++;
90 PetscFunctionReturn(PETSC_SUCCESS);
91 }
92
93 /* Computes the wind speed using Weibull distribution */
WindSpeeds(AppCtx * user)94 PetscErrorCode WindSpeeds(AppCtx *user)
95 {
96 PetscScalar *x, *t, avg_dev, sum;
97 PetscInt i;
98
99 PetscFunctionBegin;
100 user->cw = 5;
101 user->kw = 2; /* Rayleigh distribution */
102 user->nsamples = 2000;
103 user->Tw = 0.2;
104 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Wind Speed Options", "");
105 {
106 PetscCall(PetscOptionsReal("-cw", "", "", user->cw, &user->cw, NULL));
107 PetscCall(PetscOptionsReal("-kw", "", "", user->kw, &user->kw, NULL));
108 PetscCall(PetscOptionsInt("-nsamples", "", "", user->nsamples, &user->nsamples, NULL));
109 PetscCall(PetscOptionsReal("-Tw", "", "", user->Tw, &user->Tw, NULL));
110 }
111 PetscOptionsEnd();
112 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->wind_data));
113 PetscCall(VecSetSizes(user->wind_data, PETSC_DECIDE, user->nsamples));
114 PetscCall(VecSetFromOptions(user->wind_data));
115 PetscCall(VecDuplicate(user->wind_data, &user->t_wind));
116
117 PetscCall(VecGetArray(user->t_wind, &t));
118 for (i = 0; i < user->nsamples; i++) t[i] = (i + 1) * tmax / user->nsamples;
119 PetscCall(VecRestoreArray(user->t_wind, &t));
120
121 /* Wind speed deviation = (-log(rand)/cw)^(1/kw) */
122 PetscCall(VecSetRandom(user->wind_data, NULL));
123 PetscCall(VecLog(user->wind_data));
124 PetscCall(VecScale(user->wind_data, -1 / user->cw));
125 PetscCall(VecGetArray(user->wind_data, &x));
126 for (i = 0; i < user->nsamples; i++) x[i] = PetscPowScalar(x[i], 1 / user->kw);
127 PetscCall(VecRestoreArray(user->wind_data, &x));
128 PetscCall(VecSum(user->wind_data, &sum));
129 avg_dev = sum / user->nsamples;
130 /* Wind speed (t) = (1 + wind speed deviation(t) - avg_dev)*average wind speed */
131 PetscCall(VecShift(user->wind_data, 1 - avg_dev));
132 PetscCall(VecScale(user->wind_data, vwa));
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
136 /* Sets the parameters for wind turbine */
SetWindTurbineParams(AppCtx * user)137 PetscErrorCode SetWindTurbineParams(AppCtx *user)
138 {
139 PetscFunctionBegin;
140 user->Rt = 35;
141 user->Ar = PETSC_PI * user->Rt * user->Rt;
142 user->nGB = 1.0 / 89.0;
143 user->rho = 1.225;
144 user->Ht = 1.5;
145 PetscFunctionReturn(PETSC_SUCCESS);
146 }
147
148 /* Sets the parameters for induction generator */
SetInductionGeneratorParams(AppCtx * user)149 PetscErrorCode SetInductionGeneratorParams(AppCtx *user)
150 {
151 PetscFunctionBegin;
152 user->np = 4;
153 user->Xm = 3.0;
154 user->Xs = 0.1;
155 user->Xr = 0.08;
156 user->Rs = 0.01;
157 user->Rr = 0.01;
158 user->Xp = user->Xs + user->Xm * user->Xr / (user->Xm + user->Xr);
159 user->Hm = 1.0;
160 user->Te = 0.011063063063251968;
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
164 /* Computes the power extracted from wind */
GetWindPower(PetscScalar wm,PetscScalar vw,PetscScalar * Pw,AppCtx * user)165 PetscErrorCode GetWindPower(PetscScalar wm, PetscScalar vw, PetscScalar *Pw, AppCtx *user)
166 {
167 PetscScalar temp, lambda, lambda_i, cp;
168
169 PetscFunctionBegin;
170 temp = user->nGB * 2 * user->Rt * ws / user->np;
171 lambda = temp * wm / vw;
172 lambda_i = 1 / (1 / lambda + 0.002);
173 cp = 0.44 * (125 / lambda_i - 6.94) * PetscExpScalar(-16.5 / lambda_i);
174 *Pw = 0.5 * user->rho * cp * user->Ar * vw * vw * vw / (MVAbase * 1e6);
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
178 /*
179 Defines the ODE passed to the ODE solver
180 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx * user)181 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *user)
182 {
183 PetscScalar *f, wm, Pw, *wd;
184 const PetscScalar *u, *udot;
185 PetscInt stepnum;
186
187 PetscFunctionBegin;
188 PetscCall(TSGetStepNumber(ts, &stepnum));
189 /* The next three lines allow us to access the entries of the vectors directly */
190 PetscCall(VecGetArrayRead(U, &u));
191 PetscCall(VecGetArrayRead(Udot, &udot));
192 PetscCall(VecGetArray(F, &f));
193 PetscCall(VecGetArray(user->wind_data, &wd));
194
195 f[0] = user->Tw * udot[0] - wd[stepnum] + u[0];
196 wm = 1 - u[1];
197 PetscCall(GetWindPower(wm, u[0], &Pw, user));
198 f[1] = 2.0 * (user->Ht + user->Hm) * udot[1] - Pw / wm + user->Te;
199
200 PetscCall(VecRestoreArray(user->wind_data, &wd));
201 PetscCall(VecRestoreArrayRead(U, &u));
202 PetscCall(VecRestoreArrayRead(Udot, &udot));
203 PetscCall(VecRestoreArray(F, &f));
204 PetscFunctionReturn(PETSC_SUCCESS);
205 }
206
main(int argc,char ** argv)207 int main(int argc, char **argv)
208 {
209 TS ts; /* ODE integrator */
210 Vec U; /* solution will be stored here */
211 Mat A; /* Jacobian matrix */
212 PetscMPIInt size;
213 PetscInt n = 2, idx;
214 AppCtx user;
215 PetscScalar *u;
216 SNES snes;
217 PetscScalar *mat;
218 const PetscScalar *x, *rmat;
219 Mat B;
220 PetscScalar *amat;
221 PetscViewer viewer;
222
223 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
224 Initialize program
225 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
226 PetscFunctionBeginUser;
227 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
228 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
229 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
230
231 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232 Create necessary matrix and vectors
233 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
235 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
236 PetscCall(MatSetFromOptions(A));
237 PetscCall(MatSetUp(A));
238
239 PetscCall(MatCreateVecs(A, &U, NULL));
240
241 /* Create wind speed data using Weibull distribution */
242 PetscCall(WindSpeeds(&user));
243 /* Set parameters for wind turbine and induction generator */
244 PetscCall(SetWindTurbineParams(&user));
245 PetscCall(SetInductionGeneratorParams(&user));
246
247 PetscCall(VecGetArray(U, &u));
248 u[0] = vwa;
249 u[1] = s;
250 PetscCall(VecRestoreArray(U, &u));
251
252 /* Create matrix to save solutions at each time step */
253 user.stepnum = 0;
254
255 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 3, 2010, NULL, &user.Sol));
256
257 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258 Create timestepping solver context
259 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
261 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
262 PetscCall(TSSetType(ts, TSBEULER));
263 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &user));
264
265 PetscCall(TSGetSNES(ts, &snes));
266 PetscCall(SNESSetJacobian(snes, A, A, SNESComputeJacobianDefault, NULL));
267 /* PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobianFn *)IJacobian,&user)); */
268 PetscCall(TSSetApplicationContext(ts, &user));
269
270 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271 Set initial conditions
272 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
273 PetscCall(TSSetSolution(ts, U));
274
275 /* Save initial solution */
276 idx = 3 * user.stepnum;
277
278 PetscCall(MatDenseGetArray(user.Sol, &mat));
279 PetscCall(VecGetArrayRead(U, &x));
280
281 mat[idx] = 0.0;
282
283 PetscCall(PetscArraycpy(mat + idx + 1, x, 2));
284 PetscCall(MatDenseRestoreArray(user.Sol, &mat));
285 PetscCall(VecRestoreArrayRead(U, &x));
286 user.stepnum++;
287
288 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289 Set solver options
290 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291 PetscCall(TSSetMaxTime(ts, 20.0));
292 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
293 PetscCall(TSSetTimeStep(ts, .01));
294 PetscCall(TSSetFromOptions(ts));
295 PetscCall(TSSetPostStep(ts, SaveSolution));
296 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
297 Solve nonlinear system
298 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
299 PetscCall(TSSolve(ts, U));
300
301 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 3, user.stepnum, NULL, &B));
302 PetscCall(MatDenseGetArrayRead(user.Sol, &rmat));
303 PetscCall(MatDenseGetArray(B, &amat));
304 PetscCall(PetscArraycpy(amat, rmat, user.stepnum * 3));
305 PetscCall(MatDenseRestoreArray(B, &amat));
306 PetscCall(MatDenseRestoreArrayRead(user.Sol, &rmat));
307
308 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "out.bin", FILE_MODE_WRITE, &viewer));
309 PetscCall(MatView(B, viewer));
310 PetscCall(PetscViewerDestroy(&viewer));
311 PetscCall(MatDestroy(&user.Sol));
312 PetscCall(MatDestroy(&B));
313 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
314 Free work space. All PETSc objects should be destroyed when they are no longer needed.
315 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
316 PetscCall(VecDestroy(&user.wind_data));
317 PetscCall(VecDestroy(&user.t_wind));
318 PetscCall(MatDestroy(&A));
319 PetscCall(VecDestroy(&U));
320 PetscCall(TSDestroy(&ts));
321
322 PetscCall(PetscFinalize());
323 return 0;
324 }
325
326 /*TEST
327
328 build:
329 requires: !complex
330
331 test:
332 output_file: output/empty.out
333
334 TEST*/
335