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