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