HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
FPPowerSystem3BusFunctions.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <basic.h>
9 #include <math.h>
10 #include <arrayfunctions.h>
11 #include <physicalmodels/fppowersystem3bus.h>
12 
14 static void ComputeElectricalPower(
15  double theta1,
16  double theta2,
17  void *p,
18  double *Pe1,
19  double *Pe2,
20  double *Pe3
21  )
22 {
23  FPPowerSystem3Bus *params = (FPPowerSystem3Bus*) p;
24 
25  double E1 = params->E1;
26  double E2 = params->E2;
27  double E3 = params->Eref;
28 
29  double *G = params->G;
30  double *B = params->B;
31 
32  double Eph[3][2];
33  Eph[0][0] = E1*cos(theta1); Eph[0][1] = E1*sin(theta1);
34  Eph[1][0] = E2*cos(theta2); Eph[1][1] = E2*sin(theta2);
35  Eph[2][0] = E3; Eph[2][1] = 0.0;
36 
37  double Y[3][3][2];
38  int i,j;
39  for (i=0; i<3; i++) {
40  for (j=0; j<3; j++) {
41  Y[i][j][0] = G[i*3+j];
42  Y[i][j][1] = B[i*3+j];
43  }
44  }
45 
46  double YEph[3][2];
47  YEph[0][0] = Y[0][0][0]*Eph[0][0] - Y[0][0][1]*Eph[0][1]
48  + Y[0][1][0]*Eph[1][0] - Y[0][1][1]*Eph[1][1]
49  + Y[0][2][0]*Eph[2][0] - Y[0][2][1]*Eph[2][1];
50  YEph[0][1] = Y[0][0][0]*Eph[0][1] + Y[0][0][1]*Eph[0][0]
51  + Y[0][1][0]*Eph[1][1] + Y[0][1][1]*Eph[1][0]
52  + Y[0][2][0]*Eph[2][1] + Y[0][2][1]*Eph[2][0];
53  YEph[1][0] = Y[1][0][0]*Eph[0][0] - Y[1][0][1]*Eph[0][1]
54  + Y[1][1][0]*Eph[1][0] - Y[1][1][1]*Eph[1][1]
55  + Y[1][2][0]*Eph[2][0] - Y[1][2][1]*Eph[2][1];
56  YEph[1][1] = Y[1][0][0]*Eph[0][1] + Y[1][0][1]*Eph[0][0]
57  + Y[1][1][0]*Eph[1][1] + Y[1][1][1]*Eph[1][0]
58  + Y[1][2][0]*Eph[2][1] + Y[1][2][1]*Eph[2][0];
59  YEph[2][0] = Y[2][0][0]*Eph[0][0] - Y[2][0][1]*Eph[0][1]
60  + Y[2][1][0]*Eph[1][0] - Y[2][1][1]*Eph[1][1]
61  + Y[2][2][0]*Eph[2][0] - Y[2][2][1]*Eph[2][1];
62  YEph[2][1] = Y[2][0][0]*Eph[0][1] + Y[2][0][1]*Eph[0][0]
63  + Y[2][1][0]*Eph[1][1] + Y[2][1][1]*Eph[1][0]
64  + Y[2][2][0]*Eph[2][1] + Y[2][2][1]*Eph[2][0];
65 
66  YEph[0][1] = - YEph[0][1];
67  YEph[1][1] = - YEph[1][1];
68  YEph[2][1] = - YEph[2][1];
69 
70  *Pe1 = Eph[0][0]*YEph[0][0] - Eph[0][1]*YEph[0][1];
71  *Pe2 = Eph[1][0]*YEph[1][0] - Eph[1][1]*YEph[1][1];
72  *Pe3 = Eph[2][0]*YEph[2][0] - Eph[2][1]*YEph[2][1];
73 }
74 
76 int FPPowerSystem3BusDriftFunction(
77  int dir,
78  void *p,
79  double *x,
80  double t,
81  double *drift
82  )
83 {
84  FPPowerSystem3Bus *params = (FPPowerSystem3Bus*) p;
85 
86  double theta1 = x[0];
87  double theta2 = x[1];
88  double Omega1 = x[2];
89  double Omega2 = x[3];
90 
91  double omegaB = params->omegaB;
92  double Pm1_avg = params->Pm1_avg;
93  double Pm2_avg = params->Pm2_avg;
94  double Pmref_avg = params->Pmref_avg;
95  double H1 = params->H1;
96  double H2 = params->H2;
97  double Href = params->Href;
98  double gamma = params->gamma;
99 
100  double Pe1, Pe2, Peref;
101  ComputeElectricalPower(theta1,theta2,params,&Pe1,&Pe2,&Peref);
102 
103  double F1 = Pm1_avg / (2*H1) - Pmref_avg / (2*Href);
104  double F2 = Pm2_avg / (2*H2) - Pmref_avg / (2*Href);
105  double S1 = Pe1 / (2*H1) - Peref / (2*Href);
106  double S2 = Pe2 / (2*H2) - Peref / (2*Href);
107 
108  drift[0] = omegaB * Omega1;
109  drift[1] = omegaB * Omega2;
110  drift[2] = F1 - gamma*Omega1 - S1;
111  drift[3] = F2 - gamma*Omega2 - S2;
112 
113  return(0);
114 }
115 
117 int FPPowerSystem3BusDissipationFunction(
118  int dir1,
119  int dir2,
120  void *p,
121  double t,
122  double *dissp
124  )
125 {
126  FPPowerSystem3Bus *params = (FPPowerSystem3Bus*) p;
127  _ArraySetValue_(dissp,_MODEL_NDIMS_*_MODEL_NDIMS_,0.0);
128 
129  double sigma11 = params->sigma[0][0];
130  double sigma12 = params->sigma[0][1];
131  double sigma21 = params->sigma[1][0];
132  double sigma22 = params->sigma[1][1];
133 
134  double lambda11 = params->lambda[0][0];
135  double lambda12 = params->lambda[0][1];
136  double lambda21 = params->lambda[1][0];
137  double lambda22 = params->lambda[1][1];
138 
139  double gamma = params->gamma;
140  double omegaB = params->omegaB;
141 
142 #if 0
143  /* steady state coefficients */
144  dissp[2*_MODEL_NDIMS_+0] = sigma11*sigma11*lambda11*lambda11*omegaB;
145  dissp[2*_MODEL_NDIMS_+1] = sigma12*sigma12*lambda12*lambda12*omegaB;
146  dissp[3*_MODEL_NDIMS_+0] = sigma21*sigma21*lambda21*lambda21*omegaB;
147  dissp[3*_MODEL_NDIMS_+1] = sigma22*sigma22*lambda22*lambda22*omegaB;
148 
149  dissp[2*_MODEL_NDIMS_+2] = sigma11*sigma11*lambda11*(1.0-gamma*lambda11);
150  dissp[2*_MODEL_NDIMS_+3] = sigma12*sigma12*lambda12*(1.0-gamma*lambda12);
151  dissp[3*_MODEL_NDIMS_+2] = sigma21*sigma21*lambda21*(1.0-gamma*lambda21);
152  dissp[3*_MODEL_NDIMS_+3] = sigma22*sigma22*lambda22*(1.0-gamma*lambda22);
153 #endif
154 
155  /* time-dependent coefficients */
156  dissp[2*_MODEL_NDIMS_+0] = sigma11*sigma11*lambda11*omegaB*(lambda11*(1-exp(-t/lambda11))-t*exp(-t/lambda11));
157  dissp[2*_MODEL_NDIMS_+1] = sigma12*sigma12*lambda12*omegaB*(lambda12*(1-exp(-t/lambda12))-t*exp(-t/lambda12));
158  dissp[3*_MODEL_NDIMS_+0] = sigma21*sigma21*lambda21*omegaB*(lambda21*(1-exp(-t/lambda21))-t*exp(-t/lambda21));
159  dissp[3*_MODEL_NDIMS_+1] = sigma22*sigma22*lambda22*omegaB*(lambda22*(1-exp(-t/lambda22))-t*exp(-t/lambda22));
160 
161  dissp[2*_MODEL_NDIMS_+2] = sigma11*sigma11*(lambda11*(1-exp(-t/lambda11))+gamma*lambda11*(t*exp(-t/lambda11)-lambda11*(1-exp(-t/lambda11))));
162  dissp[2*_MODEL_NDIMS_+3] = sigma12*sigma12*(lambda12*(1-exp(-t/lambda12))+gamma*lambda12*(t*exp(-t/lambda12)-lambda12*(1-exp(-t/lambda12))));
163  dissp[3*_MODEL_NDIMS_+2] = sigma21*sigma21*(lambda21*(1-exp(-t/lambda21))+gamma*lambda21*(t*exp(-t/lambda21)-lambda21*(1-exp(-t/lambda21))));
164  dissp[3*_MODEL_NDIMS_+3] = sigma22*sigma22*(lambda22*(1-exp(-t/lambda22))+gamma*lambda22*(t*exp(-t/lambda22)-lambda22*(1-exp(-t/lambda22))));
165 
166  return(0);
167 }
ComputeElectricalPower
static void ComputeElectricalPower(double theta1, double theta2, void *p, double *Pe1, double *Pe2, double *Pe3)
Definition: FPPowerSystem3BusFunctions.c:14
basic.h
Some basic definitions and macros.
drift
#define drift(x)
Definition: fpdoublewell.h:31
FPPowerSystem3BusDissipationFunction
int FPPowerSystem3BusDissipationFunction(int dir1, int dir2, void *p, double t, double *dissp)
Definition: FPPowerSystem3BusFunctions.c:117
fppowersystem3bus.h
3-Bus Power System model
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
FPPowerSystem3Bus
Structure containing variable and parameters specific to the 3-bus power system model. This structure contains the physical parameters and variables for the Fokker-Planck model for a 3-bus power system.
Definition: fppowersystem3bus.h:49
arrayfunctions.h
Contains macros and function definitions for common array operations.
FPPowerSystem3BusDriftFunction
int FPPowerSystem3BusDriftFunction(int dir, void *p, double *x, double t, double *drift)
Definition: FPPowerSystem3BusFunctions.c:76