# Coverage for pySDC/projects/DAE/problems/synchronousMachine.py: 96%

## 75 statements

, created at 2024-09-09 14:59 +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)

154 # data file must be generated and stored under misc/data and self.t_end = t[-1]

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)

328# # data file must be generated and stored under misc/data and self.t_end = t[-1]

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

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