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