Actual source code: ex5.c
1: static char help[] = "Basic equation for an induction generator driven by a wind turbine.\n";
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
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
19: - This example does not work with adaptive time stepping!
21: Reference:
22: Power System Modeling and Scripting - F. Milano
23: */
25: #include <petscts.h>
27: #define freq 50
28: #define ws (2 * PETSC_PI * freq)
29: #define MVAbase 100
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 */
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 */
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 */
58: Mat Sol; /* Solution matrix */
59: PetscInt stepnum; /* Column number of solution matrix */
60: } AppCtx;
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;
68: /* Saves the solution at each time to a matrix */
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;
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: }
93: /* Computes the wind speed using Weibull distribution */
94: PetscErrorCode WindSpeeds(AppCtx *user)
95: {
96: PetscScalar *x, *t, avg_dev, sum;
98: PetscFunctionBegin;
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: PetscCall(PetscOptionsReal("-cw", "", "", user->cw, &user->cw, NULL));
106: PetscCall(PetscOptionsReal("-kw", "", "", user->kw, &user->kw, NULL));
107: PetscCall(PetscOptionsInt("-nsamples", "", "", user->nsamples, &user->nsamples, NULL));
108: PetscCall(PetscOptionsReal("-Tw", "", "", user->Tw, &user->Tw, NULL));
109: }
110: PetscOptionsEnd();
111: PetscCall(VecCreate(PETSC_COMM_WORLD, &user->wind_data));
112: PetscCall(VecSetSizes(user->wind_data, PETSC_DECIDE, user->nsamples));
113: PetscCall(VecSetFromOptions(user->wind_data));
114: PetscCall(VecDuplicate(user->wind_data, &user->t_wind));
116: PetscCall(VecGetArray(user->t_wind, &t));
117: for (PetscInt i = 0; i < user->nsamples; i++) t[i] = (i + 1) * tmax / user->nsamples;
118: PetscCall(VecRestoreArray(user->t_wind, &t));
120: /* Wind speed deviation = (-log(rand)/cw)^(1/kw) */
121: PetscCall(VecSetRandom(user->wind_data, NULL));
122: PetscCall(VecLog(user->wind_data));
123: PetscCall(VecScale(user->wind_data, -1 / user->cw));
124: PetscCall(VecGetArray(user->wind_data, &x));
125: for (PetscInt i = 0; i < user->nsamples; i++) x[i] = PetscPowScalar(x[i], 1 / user->kw);
126: PetscCall(VecRestoreArray(user->wind_data, &x));
127: PetscCall(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: PetscCall(VecShift(user->wind_data, 1 - avg_dev));
131: PetscCall(VecScale(user->wind_data, vwa));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: /* Sets the parameters for wind turbine */
136: PetscErrorCode SetWindTurbineParams(AppCtx *user)
137: {
138: PetscFunctionBegin;
139: user->Rt = 35;
140: user->Ar = PETSC_PI * user->Rt * user->Rt;
141: user->nGB = 1.0 / 89.0;
142: user->rho = 1.225;
143: user->Ht = 1.5;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /* Sets the parameters for induction generator */
148: PetscErrorCode SetInductionGeneratorParams(AppCtx *user)
149: {
150: PetscFunctionBegin;
151: user->np = 4;
152: user->Xm = 3.0;
153: user->Xs = 0.1;
154: user->Xr = 0.08;
155: user->Rs = 0.01;
156: user->Rr = 0.01;
157: user->Xp = user->Xs + user->Xm * user->Xr / (user->Xm + user->Xr);
158: user->Hm = 1.0;
159: user->Te = 0.011063063063251968;
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /* Computes the power extracted from wind */
164: PetscErrorCode GetWindPower(PetscScalar wm, PetscScalar vw, PetscScalar *Pw, AppCtx *user)
165: {
166: PetscScalar temp, lambda, lambda_i, cp;
168: PetscFunctionBegin;
169: temp = user->nGB * 2 * user->Rt * ws / user->np;
170: lambda = temp * wm / vw;
171: lambda_i = 1 / (1 / lambda + 0.002);
172: cp = 0.44 * (125 / lambda_i - 6.94) * PetscExpScalar(-16.5 / lambda_i);
173: *Pw = 0.5 * user->rho * cp * user->Ar * vw * vw * vw / (MVAbase * 1e6);
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*
178: Defines the ODE passed to the ODE solver
179: */
180: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *user)
181: {
182: PetscScalar *f, wm, Pw, *wd;
183: const PetscScalar *u, *udot;
184: PetscInt stepnum;
186: PetscFunctionBegin;
187: PetscCall(TSGetStepNumber(ts, &stepnum));
188: /* The next three lines allow us to access the entries of the vectors directly */
189: PetscCall(VecGetArrayRead(U, &u));
190: PetscCall(VecGetArrayRead(Udot, &udot));
191: PetscCall(VecGetArray(F, &f));
192: PetscCall(VecGetArray(user->wind_data, &wd));
194: f[0] = user->Tw * udot[0] - wd[stepnum] + u[0];
195: wm = 1 - u[1];
196: PetscCall(GetWindPower(wm, u[0], &Pw, user));
197: f[1] = 2.0 * (user->Ht + user->Hm) * udot[1] - Pw / wm + user->Te;
199: PetscCall(VecRestoreArray(user->wind_data, &wd));
200: PetscCall(VecRestoreArrayRead(U, &u));
201: PetscCall(VecRestoreArrayRead(Udot, &udot));
202: PetscCall(VecRestoreArray(F, &f));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: int main(int argc, char **argv)
207: {
208: TS ts; /* ODE integrator */
209: Vec U; /* solution will be stored here */
210: Mat A; /* Jacobian matrix */
211: PetscMPIInt size;
212: PetscInt n = 2, idx;
213: AppCtx user;
214: PetscScalar *u;
215: SNES snes;
216: PetscScalar *mat;
217: const PetscScalar *x, *rmat;
218: Mat B;
219: PetscScalar *amat;
220: PetscViewer viewer;
222: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
223: Initialize program
224: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
225: PetscFunctionBeginUser;
226: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
227: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
228: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
230: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: Create necessary matrix and vectors
232: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
234: PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
235: PetscCall(MatSetFromOptions(A));
236: PetscCall(MatSetUp(A));
238: PetscCall(MatCreateVecs(A, &U, NULL));
240: /* Create wind speed data using Weibull distribution */
241: PetscCall(WindSpeeds(&user));
242: /* Set parameters for wind turbine and induction generator */
243: PetscCall(SetWindTurbineParams(&user));
244: PetscCall(SetInductionGeneratorParams(&user));
246: PetscCall(VecGetArray(U, &u));
247: u[0] = vwa;
248: u[1] = s;
249: PetscCall(VecRestoreArray(U, &u));
251: /* Create matrix to save solutions at each time step */
252: user.stepnum = 0;
254: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 3, 2010, NULL, &user.Sol));
256: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257: Create timestepping solver context
258: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
259: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
260: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
261: PetscCall(TSSetType(ts, TSBEULER));
262: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &user));
264: PetscCall(TSGetSNES(ts, &snes));
265: PetscCall(SNESSetJacobian(snes, A, A, SNESComputeJacobianDefault, NULL));
266: /* PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobianFn *)IJacobian,&user)); */
267: PetscCall(TSSetApplicationContext(ts, &user));
269: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: Set initial conditions
271: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
272: PetscCall(TSSetSolution(ts, U));
274: /* Save initial solution */
275: idx = 3 * user.stepnum;
277: PetscCall(MatDenseGetArray(user.Sol, &mat));
278: PetscCall(VecGetArrayRead(U, &x));
280: mat[idx] = 0.0;
282: PetscCall(PetscArraycpy(mat + idx + 1, x, 2));
283: PetscCall(MatDenseRestoreArray(user.Sol, &mat));
284: PetscCall(VecRestoreArrayRead(U, &x));
285: user.stepnum++;
287: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288: Set solver options
289: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
290: PetscCall(TSSetMaxTime(ts, 20.0));
291: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
292: PetscCall(TSSetTimeStep(ts, .01));
293: PetscCall(TSSetFromOptions(ts));
294: PetscCall(TSSetPostStep(ts, SaveSolution));
295: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
296: Solve nonlinear system
297: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
298: PetscCall(TSSolve(ts, U));
300: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 3, user.stepnum, NULL, &B));
301: PetscCall(MatDenseGetArrayRead(user.Sol, &rmat));
302: PetscCall(MatDenseGetArray(B, &amat));
303: PetscCall(PetscArraycpy(amat, rmat, user.stepnum * 3));
304: PetscCall(MatDenseRestoreArray(B, &amat));
305: PetscCall(MatDenseRestoreArrayRead(user.Sol, &rmat));
307: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "out.bin", FILE_MODE_WRITE, &viewer));
308: PetscCall(MatView(B, viewer));
309: PetscCall(PetscViewerDestroy(&viewer));
310: PetscCall(MatDestroy(&user.Sol));
311: PetscCall(MatDestroy(&B));
312: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313: Free work space. All PETSc objects should be destroyed when they are no longer needed.
314: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
315: PetscCall(VecDestroy(&user.wind_data));
316: PetscCall(VecDestroy(&user.t_wind));
317: PetscCall(MatDestroy(&A));
318: PetscCall(VecDestroy(&U));
319: PetscCall(TSDestroy(&ts));
321: PetscCall(PetscFinalize());
322: return 0;
323: }
325: /*TEST
327: build:
328: requires: !complex
330: test:
331: output_file: output/empty.out
333: TEST*/