Coverage for pySDC/projects/DAE/problems/synchronous_machine.py: 96%
76 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
2import warnings
3from scipy.interpolate import interp1d
5from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
6from pySDC.implementations.datatype_classes.mesh import mesh
9class synchronous_machine_infinite_bus(ptype_dae):
10 r"""
11 Synchronous machine model from Kundur (equiv. circuits fig. 3.18 in [1]_) attached to infinite bus. The machine can be
12 represented as two different circuits at the direct-axis and the quadrature-axis. Detailed information can be found in
13 [1]_. The system of differential-algebraic equations (DAEs) consists of the equations for
15 - the stator voltage equations
17 .. math::
18 \frac{d \Psi_d (t)}{dt} = \omega_b (v_d + R_a i_d (t) + \omega_r \Psi_q (t)),
20 .. math::
21 \frac{d \Psi_q (t)}{dt} = \omega_b (v_q + R_a i_q (t) - \omega_r \Psi_d (t)),
23 .. math::
24 \frac{d \Psi_0 (t)}{dt} = \omega_b (v_0 + R_a i_0 (t)),
26 - the rotor voltage equations
28 .. math::
29 \frac{d \Psi_F (t)}{dt} = \omega_b (v_F - R_F i_F (t)),
31 .. math::
32 \frac{d \Psi_D (t)}{dt} = -\omega_b (R_D i_D (t)),
34 .. math::
35 \frac{d \Psi_{Q1} (t)}{dt} = -\omega_b (R_{Q1} i_{Q1} (t)),
37 .. math::
38 \frac{d \Psi_{Q2} (t)}{dt} = -\omega_b (R_{Q2} i_{Q2} (t)),
40 - the stator flux linkage equations
42 .. math::
43 \Psi_d (t) = L_d i_d (t) + L_{md} i_F (t) + L_{md} i_D (t),
45 .. math::
46 \Psi_q (t) = L_q i_q (t) + L_{mq} i_{Q1} (t) + L_{mq} i_{Q2} (t),
48 .. math::
49 \Psi_0 (t) = L_0 i_0 (t)
51 - the rotor flux linkage equations
53 .. math::
54 \Psi_F = L_F i_F (t) + L_D i_D + L_{md} i_d (t),
56 .. math::
57 \Psi_D = L_F i_F (t) + L_D i_D + L_{md} i_d (t),
59 .. math::
60 \Psi_{Q1} = L_{Q1} i_{Q1} (t) + L_{mq} i_{Q2} + L_{mq} i_q (t),
62 .. math::
63 \Psi_{Q2} = L_{mq} i_{Q1} (t) + L_{Q2} i_{Q2} + L_{mq} i_q (t),
65 - the swing equations
67 .. math::
68 \frac{d \delta (t)}{dt} = \omega_b (\omega_r (t) - 1),
70 .. math::
71 \frac{d \omega_r (t)}{dt} = \frac{1}{2 H}(T_m - T_e - K_D \omega_b (\omega_r (t) - 1)).
73 The voltages :math:`v_d`, :math:`v_q` can be updated via the following procedure. The stator's currents are mapped
74 to the comlex-valued external reference frame current :math:`I` with
76 .. math::
77 \Re(I) = i_d (t) \sin(\delta (t)) + i_q (t) \cos(\delta (t)),
79 .. math::
80 \Im(I) = -i_d (t) \cos(\delta (t)) + i_q (t) \sin(\delta (t)).
82 The voltage V across the stator terminals can then be computed as complex-value via
84 .. math::
85 V_{comp} = E_B + Z_{line} (\Re(I) + i \Im(I))
87 with impedance :math:`Z_{line}\in\mathbb{C}`. Then, :math:`v_d`, :math:`v_q` can be computed via the network equations
89 .. math::
90 v_d = \Re(V_{comp}) \sin(\delta (t)) - \Im(V_{comp}) \cos(\delta (t)),
92 .. math::
93 v_q = \Re(V_{comp}) \cos(\delta (t)) + \Im(V_{comp}) \sin(\delta (t)),
95 which describes the connection between the machine and the infinite bus.
97 Parameters
98 ----------
99 nvars : int
100 Number of unknowns of the system of DAEs.
101 newton_tol : float
102 Tolerance for Newton solver.
104 Attributes
105 ----------
106 L_d: float
107 Inductance of inductor :math:'L_d', see [1]_.
108 L_q: float
109 Inductance of inductor :math:'L_q', see [1]_.
110 L_F: float
111 Inductance of inductor :math:'L_F', see [1]_.
112 L_D: float
113 Inductance of inductor :math:'L_D', see [1]_.
114 L_Q1: float
115 Inductance of inductor :math:'L_{Q1}', see [1]_.
116 L_Q2: float
117 Inductance of inductor :math:'L_{Q2}', see [1]_.
118 L_md: float
119 Inductance of inductor :math:'L_{md}', see [1]_.
120 L_mq: float
121 Inductance of inductor :math:'L_{mq}', see [1]_.
122 R_s: float
123 Resistance of resistor :math:`R_s`, see [1]_.
124 R_F: float
125 Resistance of resistor :math:`R_F`, see [1]_.
126 R_D: float
127 Resistance of resistor :math:`R_D`, see [1]_.
128 R_Q1: float
129 Resistance of resistor :math:`R_{Q1}`, see [1]_.
130 R_Q2: float
131 Resistance of resistor :math:`R_{Q2}`, see [1]_.
132 omega_b: float
133 Base frequency of the rotor in mechanical :math:`rad/s`.
134 H_: float
135 Defines the per unit inertia constant.
136 K_D: float
137 Factor that accounts for damping losses.
138 Z_line: complex
139 Impedance of the transmission line that connects the infinite bus to the generator.
140 E_B: float
141 Voltage of infinite bus.
142 v_F: float
143 Voltage at the field winding.
144 T_m: float
145 Defines the mechanical torque applied to the rotor shaft.
147 References
148 ----------
149 .. [1] P. Kundur, N. J. Balu, M. G. Lauby. Power system stability and control. The EPRI power system series (1994).
150 """
152 def __init__(self, newton_tol):
153 super().__init__(nvars=14, newton_tol=newton_tol)
154 # load reference solution
155 # data file must be generated and stored under misc/data and self.t_end = t[-1]
156 # data = np.load(r'pySDC/projects/DAE/misc/data/synch_gen.npy')
157 # x = data[:, 0]
158 # y = data[:, 1:]
159 # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate')
160 self.t_end = 0.0
162 self.L_d = 1.8099
163 self.L_q = 1.76
164 self.L_F = 1.8247
165 self.L_D = 1.8312
166 self.L_Q1 = 2.3352
167 self.L_Q2 = 1.735
168 self.L_md = 1.6599
169 self.L_mq = 1.61
170 self.R_s = 0.003
171 self.R_F = 0.0006
172 self.R_D = 0.0284
173 self.R_Q1 = 0.0062
174 self.R_Q2 = 0.0237
175 self.omega_b = 376.9911184307752
176 self.H_ = 3.525
177 self.K_D = 0.0
178 # Line impedance
179 self.Z_line = -0.2688022164909709 - 0.15007173591230372j
180 # Infinite bus voltage
181 self.E_B = 0.7
182 # Rotor (field) operating voltages
183 # These are modelled as constants. Intuition: permanent magnet as rotor
184 self.v_F = 8.736809687330562e-4
185 self.T_m = 0.854
187 def eval_f(self, u, du, t):
188 r"""
189 Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`.
191 Parameters
192 ----------
193 u : dtype_u
194 Current values of the numerical solution at time t.
195 du : dtype_u
196 Current values of the derivative of the numerical solution at time t.
197 t : float
198 Current time of the numerical solution.
200 Returns
201 -------
202 f : dtype_f
203 Current value of the right-hand side of f (which includes 14 components).
204 """
206 # simulate torque change at t = 0.05
207 if t >= 0.05:
208 self.T_m = 0.354
210 f = self.dtype_f(self.init)
212 # u = [psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2,
213 # i_d, i_q, i_F, i_D, i_Q1, i_Q2
214 # omega_m,
215 # v_d, v_q,
216 # iz_d, iz_q, il_d, il_q, vl_d, vl_q]
218 # extract variables for readability
219 # algebraic components
220 psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2 = u.diff[0], u.diff[1], u.diff[2], u.diff[3], u.diff[4], u.diff[5]
221 delta_r, omega_m = u.diff[6], u.diff[7]
222 i_d, i_q, i_F, i_D, i_Q1, i_Q2 = u.alg[0], u.alg[1], u.alg[2], u.alg[3], u.alg[4], u.alg[5]
224 # differential components
225 # these result directly from the voltage equations, introduced e.g. pg. 145 Krause
226 dpsi_d, dpsi_q, dpsi_F, dpsi_D, dpsi_Q1, dpsi_Q2 = (
227 du.diff[0],
228 du.diff[1],
229 du.diff[2],
230 du.diff[3],
231 du.diff[4],
232 du.diff[5],
233 )
234 ddelta_r, domega_m = du.diff[6], du.diff[7]
236 # Network current
237 I_Re = i_d * np.sin(delta_r) + i_q * np.cos(delta_r)
238 I_Im = -i_d * np.cos(delta_r) + i_q * np.sin(delta_r)
239 # Machine terminal voltages in network coordinates
240 # Need to transform like this to subtract infinite bus voltage
241 V_comp = self.E_B - self.Z_line * (-1) * (I_Re + 1j * I_Im)
242 # Terminal voltages in dq0 coordinates
243 v_d = np.real(V_comp) * np.sin(delta_r) - np.imag(V_comp) * np.cos(delta_r)
244 v_q = np.real(V_comp) * np.cos(delta_r) + np.imag(V_comp) * np.sin(delta_r)
246 f.diff[:8] = (
247 -dpsi_d + self.omega_b * (v_d - self.R_s * i_d + omega_m * psi_q),
248 -dpsi_q + self.omega_b * (v_q - self.R_s * i_q - omega_m * psi_d),
249 -dpsi_F + self.omega_b * (self.v_F - self.R_F * i_F),
250 -dpsi_D - self.omega_b * self.R_D * i_D,
251 -dpsi_Q1 - self.omega_b * self.R_Q1 * i_Q1,
252 -dpsi_Q2 - self.omega_b * self.R_Q2 * i_Q2,
253 -ddelta_r + self.omega_b * (omega_m - 1),
254 -domega_m
255 + 1 / (2 * self.H_) * (self.T_m - (psi_q * i_d - psi_d * i_q) - self.K_D * self.omega_b * (omega_m - 1)),
256 )
257 f.alg[:6] = (
258 -psi_d + self.L_d * i_d + self.L_md * i_F + self.L_md * i_D,
259 -psi_q + self.L_q * i_q + self.L_mq * i_Q1 + self.L_mq * i_Q2,
260 -psi_F + self.L_md * i_d + self.L_F * i_F + self.L_md * i_D,
261 -psi_D + self.L_md * i_d + self.L_md * i_F + self.L_D * i_D,
262 -psi_Q1 + self.L_mq * i_q + self.L_Q1 * i_Q1 + self.L_mq * i_Q2,
263 -psi_Q2 + self.L_mq * i_q + self.L_mq * i_Q1 + self.L_Q2 * i_Q2,
264 )
265 self.work_counters['rhs']()
266 return f
268 def u_exact(self, t):
269 """
270 Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution.
272 Parameters
273 ----------
274 t : float
275 The time of the reference solution.
277 Returns
278 -------
279 me : dtype_u
280 The reference solution as mesh object. It contains fixed initial conditions at initial time (which includes
281 14 components).
282 """
283 me = self.dtype_u(self.init)
285 if t == 0:
286 psi_d = 0.7770802016688648
287 psi_q = -0.6337183129426077
288 psi_F = 1.152966888216155
289 psi_D = 0.9129958488040036
290 psi_Q1 = -0.5797082294536264
291 psi_Q2 = -0.579708229453273
292 i_d = -0.9061043142342473
293 i_q = -0.36006722326230495
294 i_F = 1.45613494788927
295 i_D = 0.0
296 i_Q1 = 0.0
297 i_Q2 = 0.0
299 delta_r = 39.1 * np.pi / 180
300 omega_0 = 2 * np.pi * 60
301 omega_b = 2 * np.pi * 60
302 omega_m = omega_0 / omega_b # = omega_r since pf = 2 i.e. two pole machine
304 me.diff[:8] = (psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2, delta_r, omega_m)
305 me.alg[:6] = (i_d, i_q, i_F, i_D, i_Q1, i_Q2)
306 elif t < self.t_end:
307 u_ref = self.u_ref(t)
308 me.diff[:8] = u_ref[:8]
309 me.alg[:6] = u_ref[8:]
310 else:
311 self.logger.warning("Requested time exceeds domain of the reference solution. Returning zero.")
312 me.diff[:8] = (0, 0, 0, 0, 0, 0, 0, 0)
313 me.alg[:6] = (0, 0, 0, 0, 0, 0)
315 return me
318# class synchronous_machine_pi_line(ptype_dae):
319# """
320# Synchronous machine model from Kundur (equiv. circuits fig. 3.18)
321# attached to pi line with resistive load
323# This model does not work yet but is included as a starting point for developing similar models and in the hope that somebody will figure out why it does not work
324# """
326# def __init__(self, problem_params, dtype_u=mesh, dtype_f=mesh):
327# super(synchronous_machine_pi_line, self).__init__(problem_params, dtype_u, dtype_f)
328# # load reference solution
329# # data file must be generated and stored under misc/data and self.t_end = t[-1]
330# # data = np.load(r'pySDC/projects/DAE/misc/data/synch_gen.npy')
331# # x = data[:, 0]
332# # y = data[:, 1:]
333# # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate')
334# self.t_end = 0.0
336# self.L_d = 1.8099
337# self.L_q = 1.76
338# self.L_F = 1.8247
339# self.L_D = 1.8312
340# self.L_Q1 = 2.3352
341# self.L_Q2 = 1.735
342# self.L_md = 1.6599
343# self.L_mq = 1.61
344# self.R_s = 0.003
345# self.R_F = 0.0006
346# self.R_D = 0.0284
347# self.R_Q1 = 0.0062
348# self.R_Q2 = 0.0237
349# self.omega_b = 376.9911184307752
350# self.H_ = 3.525
351# self.K_D = 0.0
352# # pi line
353# self.C_pi = 0.000002
354# self.R_pi = 0.02
355# self.L_pi = 0.00003
356# # load
357# self.R_L = 0.75
358# self.v_F = 8.736809687330562e-4
359# self.v_D = 0
360# self.v_Q1 = 0
361# self.v_Q2 = 0
362# self.T_m = 0.854
364# def eval_f(self, u, du, t):
365# """
366# Routine to evaluate the implicit representation of the problem i.e. F(u', u, t)
367# Args:
368# u (dtype_u): the current values. This parameter has been "hijacked" to contain [u', u] in this case to enable evaluation of the implicit representation
369# t (float): current time
370# Returns:
371# Current value of F(), 21 components
372# """
374# # simulate torque change at t = 0.05
375# if t >= 0.05:
376# self.T_m = 0.354
378# f = self.dtype_f(self.init)
380# # u = [psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2,
381# # i_d, i_q, i_F, i_D, i_Q1, i_Q2
382# # omega_m,
383# # v_d, v_q,
384# # iz_d, iz_q, il_d, il_q, vl_d, vl_q]
386# # extract variables for readability
387# # algebraic components
388# psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2 = u[0], u[1], u[2], u[3], u[4], u[5]
389# i_d, i_q, i_F, i_D, i_Q1, i_Q2 = u[6], u[7], u[8], u[9], u[10], u[11]
390# # delta_r = u[12]
391# omega_m = u[12]
392# v_d, v_q = u[13], u[14]
393# iz_d, iz_q, il_d, il_q, vl_d, vl_q = u[15], u[16], u[17], u[18], u[19], u[20]
395# # differential components
396# # these result directly from the voltage equations, introduced e.g. pg. 145 Krause
397# dpsi_d, dpsi_q, dpsi_F, dpsi_D, dpsi_Q1, dpsi_Q2 = du[0], du[1], du[2], du[3], du[4], du[5]
398# # ddelta_r = du[12]
399# domega_m = du[12]
400# dv_d, dv_q = du[13], du[14]
401# diz_d, diz_q, dvl_d, dvl_q = du[15], du[16],du[19], du[20]
403# # algebraic variables are i_d, i_q, i_F, i_D, i_Q1, i_Q2, il_d, il_q
405# f[:] = (
406# # differential generator
407# dpsi_d + self.omega_b * (v_d - self.R_s * i_d + omega_m * psi_q),
408# dpsi_q + self.omega_b * (v_q - self.R_s * i_q - omega_m * psi_d),
409# dpsi_F + self.omega_b * (self.v_F - self.R_F * i_F),
410# dpsi_D + self.omega_b * (self.v_D - self.R_D * i_D),
411# dpsi_Q1 + self.omega_b * (self.v_Q1 - self.R_Q1 * i_Q1),
412# dpsi_Q2 + self.omega_b * (self.v_Q2 - self.R_Q2 * i_Q2),
413# -domega_m + 1 / (2 * self.H_) * (self.T_m - (psi_q * i_d - psi_d * i_q) - self.K_D * self.omega_b * (omega_m-1)),
414# # differential pi line
415# -dv_d + omega_m * v_q + 2/self.C_pi * (i_d - iz_d),
416# -dv_q - omega_m * v_d + 2/self.C_pi * (i_q - iz_q),
417# -dvl_d + omega_m * vl_q + 2/self.C_pi * (iz_d - il_d),
418# -dvl_q - omega_m * vl_d + 2/self.C_pi * (iz_q - il_q),
419# -diz_d - self.R_pi/self.L_pi * iz_d + omega_m * iz_q + (v_d - vl_d) / self.L_pi,
420# -diz_q - self.R_pi/self.L_pi * iz_q - omega_m * iz_d + (v_q - vl_q) / self.L_pi,
421# # algebraic generator
422# psi_d + self.L_d * i_d + self.L_md * i_F + self.L_md * i_D,
423# psi_q + self.L_q * i_q + self.L_mq * i_Q1 + self.L_mq * i_Q2,
424# psi_F + self.L_md * i_d + self.L_F * i_F + self.L_md * i_D,
425# psi_D + self.L_md * i_d + self.L_md * i_F + self.L_D * i_D,
426# psi_Q1 + self.L_mq * i_q + self.L_Q1 * i_Q1 + self.L_mq * i_Q2,
427# psi_Q2 + self.L_mq * i_q + self.L_mq * i_Q1 + self.L_Q2 * i_Q2,
428# # algebraic pi line
429# -il_d + vl_d/self.R_L,
430# -il_q + vl_q/self.R_L,
431# )
432# return f
434# def u_exact(self, t):
435# """
436# Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution.
437# Args:
438# t (float): current time
439# Returns:
440# Mesh containing fixed initial value, 5 components
441# """
442# me = self.dtype_u(self.init)
444# if t == 0:
445# psi_d = 0.3971299
446# psi_q = 0.9219154
447# psi_F = 0.8374232
448# psi_D = 0.5795112
449# psi_Q1 = 0.8433430
450# psi_Q2 = 0.8433430
451# i_d = -1.215876
452# i_q = 0.5238156
453# i_F = 1.565
454# i_D = 0
455# i_Q1 = 0
456# i_Q2 = 0
457# v_d = -0.9362397
458# v_q = 0.4033005
459# omega_m = 1.0
460# # pi line
461# iz_d = -1.215875
462# iz_q = 0.5238151
463# il_d = -1.215875
464# il_q = 0.5238147
465# vl_d = -0.9119063
466# vl_q = 0.3928611
467# me[:] = (psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2,
468# i_d, i_q, i_F, i_D, i_Q1, i_Q2,
469# omega_m,
470# v_d, v_q,
471# iz_d, iz_q, il_d, il_q, vl_d, vl_q)
472# elif t < self.t_end:
473# me[:] = self.u_ref(t)
474# else:
475# self.logger.warning("Requested time exceeds domain of the reference solution. Returning zero.")
476# me[:] = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
478# return me