Actual source code: ex5.c

      1: static char help[] = "Basic equation for an induction generator driven by a wind turbine.\n";

      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

     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

     19:  - This example does not work with adaptive time stepping!

     21: Reference:
     22: Power System Modeling and Scripting - F. Milano
     23: */

     25: #include <petscts.h>

     27: #define freq    50
     28: #define ws      (2 * PETSC_PI * freq)
     29: #define MVAbase 100

     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 */

     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 */

     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 */

     58:   Mat      Sol;     /* Solution matrix */
     59:   PetscInt stepnum; /* Column number of solution matrix */
     60: } AppCtx;

     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;

     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;

     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: }

     93: /* Computes the wind speed using Weibull distribution */
     94: PetscErrorCode WindSpeeds(AppCtx *user)
     95: {
     96:   PetscScalar *x, *t, avg_dev, sum;
     97:   PetscInt     i;

     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));

    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));

    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: }

    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: }

    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: }

    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;

    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: }

    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;

    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));

    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;

    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: }

    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;

    223:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    224:      Initialize program
    225:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    226:   PetscFunctionBeginUser;
    227:   PetscCall(PetscInitialize(&argc, &argv, NULL, 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");

    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));

    239:   PetscCall(MatCreateVecs(A, &U, NULL));

    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));

    247:   PetscCall(VecGetArray(U, &u));
    248:   u[0] = vwa;
    249:   u[1] = s;
    250:   PetscCall(VecRestoreArray(U, &u));

    252:   /* Create matrix to save solutions at each time step */
    253:   user.stepnum = 0;

    255:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, 3, 2010, NULL, &user.Sol));

    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));

    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));

    270:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    271:      Set initial conditions
    272:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    273:   PetscCall(TSSetSolution(ts, U));

    275:   /* Save initial solution */
    276:   idx = 3 * user.stepnum;

    278:   PetscCall(MatDenseGetArray(user.Sol, &mat));
    279:   PetscCall(VecGetArrayRead(U, &x));

    281:   mat[idx] = 0.0;

    283:   PetscCall(PetscArraycpy(mat + idx + 1, x, 2));
    284:   PetscCall(MatDenseRestoreArray(user.Sol, &mat));
    285:   PetscCall(VecRestoreArrayRead(U, &x));
    286:   user.stepnum++;

    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));

    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));

    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));

    322:   PetscCall(PetscFinalize());
    323:   return 0;
    324: }

    326: /*TEST

    328:    build:
    329:       requires: !complex

    331:    test:
    332:      output_file: output/empty.out

    334: TEST*/