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