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