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