wip: 5th order sync gen model updates (3f3cee65) · Commits · dt-hydro / dt-hydro-cfd

dthydro_cfd/generator_models.py

+154 −21

Original line number	Diff line number	Diff line
		@@ -62,7 +62,7 @@ def armature_voltage(R, X_q_pp, X_d_pp, E_d_pp, E_q_pp, I_d, I_q):
		return v_arr


		class sync_gen5(object):
		class SyncGen5(object):
		"""
		5th order sychronous generator model.
		Machowski. Power System Dynamics: Stability and Control. 2008.
		@@ -71,17 +71,23 @@ class sync_gen5(object):
		E_d_p := transient emf on d
		E_q_pp := subtransient emf on q
		E_d_pp := subtransient emf on d

		Args:
		w_0: Initial angular velocity [rad/s]
		w_s: Synchronous angular velocity [rad/s]
		"""
		def __init__(self, w_0, w_s=W_S, **kwargs):
		def __init__(self, w_0: float, w_s: float=W_S, **kwargs):
		# Setpoint rotation rate
		self.w_s = w_s
		# System state vector
		self.x = np.zeros(5)
		# rotation rate deviation
		self.x[0] = w_0 - w_s
		# System voltage (V), this is V_g in the PSD book
		self.V_s = kwargs.get("V_s", 13.5e3)
		# Initial currents
		self.I_d = kwargs.get("I_d", 0.)
		self.I_q = kwargs.get("I_q", 0.)
		#self.I_d = kwargs.get("I_d", 0.)
		#self.I_q = kwargs.get("I_q", 0.)
		# generator constants
		Xd = 0.768
		Xd1 = 0.249
		@@ -96,6 +102,7 @@ class sync_gen5(object):
		#Tj = 7.
		Xs = 0.3
		self.machine_params = {
		"Xs": kwargs.get("Xs", Xs),
		"Xd": kwargs.get("Xd", Xd),
		"Xd1": kwargs.get("Xd", Xd1),
		"Xd2": kwargs.get("Xd", Xd2),
		@@ -117,13 +124,15 @@ class sync_gen5(object):
		self.x = x
		# Mechanical power input
		P_m = args[0]
		# Driving voltage (field voltage)
		E_f = args[1]
		# Current inputs
		I_d = args[1]
		I_q = args[2]
		self.I_d = I_d
		self.I_q = I_q
		# Driving voltage
		E_f = args[3]
		# I_d = args[1]
		# I_q = args[2]
		# self.I_d = I_d
		# self.I_q = I_q
		I_d = self.I_d()
		I_q = self.I_q()
		# Machine constants
		T_do_p = self.machine_params["Tdo1"]
		T_do_pp = self.machine_params["Tdo2"]
		@@ -160,6 +169,14 @@ class sync_gen5(object):
		def __call__(self, t, x, args, *kwargs):
		return self.d_sys_dt(t, x, args, *kwargs)

		#@property
		#def I_d(self):
		# return self.i_d(0, 0)

		#@property
		#def I_q(self):
		# return self.i_q(0, 0)

		@property
		def power_e(self):
		""" electrical power W """
		@@ -189,13 +206,18 @@ class sync_gen5(object):
		@property
		def delta_angle(self):
		""" delta phasor angle """
		return self.x[1]
		return self.x[1] % (np.pi / 2)

		@property
		def E_q_p(self):
		""" emf """
		return self.x[2]

		@property
		def E_d_p(self):
		""" emf """
		return 0.0

		@property
		def E_q_pp(self):
		return self.x[3]
		@@ -204,6 +226,87 @@ class sync_gen5(object):
		def E_d_pp(self):
		return self.x[4]

		@property
		def angle(self):
		""" alias to delta_angle mod pi/2"""
		return self.delta_angle

		def V_dq_from_vs(self, V_s=None):
		""" compute [v_sd, v_sq] from system voltage """
		if V_s is not None:
		return np.array([-V_s * np.sin(self.angle),
		V_s * np.cos(self.angle)])
		else:
		return np.array([-self.V_s * np.sin(self.angle),
		self.V_s * np.cos(self.angle)])

		def V_dq_from_idq(self, R=0.):
		return self.arm_v(R) + self.machine_params["Xs"] * \
		np.array([self.I_q(), -self.I_d()])

		def solve_idq_vdq_sys(self):
		pass

		def V_gdq(self, R=0.):
		""" alias to arm_v """
		return self.arm_v(R)

		def I_qd(self):
		""" equation 3.24 from book """
		Xs = self.machine_params["Xs"]
		I_dq = (self.V_dq_from_vs() - self.V_gdq()) / Xs
		I_q = I_dq[0]
		I_d = I_dq[1]
		return np.array([I_q, -I_d])

		def I_q(self):
		return self.I_qd()[0]

		def I_d(self):
		return self.I_qd()[1]

		# Aux fns from DynPSSymPy
		def d(self, x, v):
		return np.exp(1j * (self.angle - np.pi / 2))

		def q(self, x, v):
		return np.exp(1j * self.angle)

		def v_t(self, x, v):
		""" terminal voltage """
		# return v[self.bus_idx_red['terminal']]
		return self.machine_params.get("v_t", 1.0)

		def e_q_st(self, x, v):
		return self.x[3]

		def e_d_st(self, x, v):
		return self.x[4]

		def e_q_t(self, x, v):
		return self.x[2]

		def e_d_t(self, x, v):
		return self.E_d_p

		def e_st(self, x, v):
		return self.e_q_st(x, v)self.q(x, v) + self.e_d_st(x, v)self.d(x, v)

		def e_t(self, x, v):
		return self.e_q_t(x, v)self.q(x, v) + self.e_d_t(x, v)self.d(x, v)

		def i(self, x, v):
		X_d_st = self.machine_params["Xd2"]
		return (self.e_st(x, v) - self.v_t(x, v)) / (1j * X_d_st)

		def i_d(self, x, v):
		i_dq = self.i(x, v)np.exp(1j(np.pi/2 - self.angle))
		return i_dq.real

		def i_q(self, x, v):
		i_dq = self.i(x, v)np.exp(1j(np.pi/2 - self.angle))
		return i_dq.imag

		def arm_v(self, R=0.0):
		""" armature_voltage """
		X_d = self.machine_params["Xd"]
		@@ -213,20 +316,19 @@ class sync_gen5(object):
		X_q_p = self.machine_params["Xq1"]
		X_q_pp = self.machine_params["Xq2"]
		return armature_voltage(
		R, X_q_pp, X_d_pp, self.E_d_pp, self.E_q_pp, self.I_d, self.I_q)
		R, X_q_pp, X_d_pp, self.E_d_pp,
		self.E_q_pp, self.I_d, self.I_q)


		def test_gen5():
		""" Test of the 5th order sychronous generator model """
		from scipy.integrate import solve_ivp
		# Driving voltage, generator main field voltage
		e_f = 860.0
		e_f = 10.0
		delta_w = 2.0 * np.pi * 1e-4 # rad / s
		w_0 = W_S + delta_w
		# Line Currents in system frame (amps)
		#i_a = 1.0
		#i_b = 1.0
		#i_c = 1.0
		# System voltage (13.5 kV)
		V_s_0 = 13.5e3

		I_0_0 = 0. # Must be 0 bc, (i_a + i_b + i_c) / sqrt(3) = 0 = i_0
		I_d_0 = 9.5e3 * 0.
		@@ -242,17 +344,18 @@ def test_gen5():
		i_0dq = np.array([I_0, I_d, I_q]) * np.sqrt(3.)
		assert np.allclose((I_0), (0.0))
		assert np.allclose((I_0, I_d, I_q), (I_0_0, I_d_0, I_q_0))
		gen5 = sync_gen5(w_0, W_S, I_d=I_d, I_q=I_q)
		gen5 = SyncGen5(w_0, W_S, V_s=V_s_0)
		x_0 = deepcopy(gen5.x)
		t = 0.

		# Mechanical power input
		P_m = 25.0e6 # W
		P_m = 30.0 # W
		# Driving voltage from excitation system
		E_f = e_f

		# test call to system derivative
		test_dx_dt = gen5(t, x_0, P_m, I_d, I_q, E_f)
		test_dx_dt = gen5(t, x_0, P_m, E_f)

		# test time integration of the system
		dt = 1e-1
		@@ -260,16 +363,24 @@ def test_gen5():
		p_e_list = []
		t_list = []
		w_list = []
		for i in range(8000):
		sol = solve_ivp(gen5, [t, t+dt], x_new, args=(P_m, I_d, I_q, E_f))
		delta_list = []
		I_q_list, I_d_list = [], []
		for i in range(300):
		sol = solve_ivp(gen5, [t, t+dt], x_new, args=(P_m, E_f))
		x_new = sol.y[:, -1]
		# reset rpm
		# gen5.x[0] = 0.
		t = t + dt
		p_e_list.append(gen5.power_e)
		w_list.append(gen5.w)
		delta_list.append(gen5.delta_angle)
		I_q_list.append(gen5.I_q)
		I_d_list.append(gen5.I_d)
		t_list.append(t)

		I_d = gen5.i_d(0, 0)
		I_q = gen5.i_q(0, 0)

		aV = gen5.arm_v(R=0.0)
		V_0, V_d, V_q = 0., aV[0], aV[1]
		# Eq. 11.80
		@@ -297,6 +408,9 @@ def test_gen5():
		t_list = np.asarray(t_list)
		p_e_list = np.asarray(p_e_list)
		w_list = np.asarray(w_list)
		delta_list = np.asarray(delta_list)
		I_q_list = np.asarray(I_q_list)
		I_d_list = np.asarray(I_d_list)
		print(sol.t)
		print(sol.y)
		sol.y[:, -1][-1]
		@@ -323,6 +437,25 @@ def test_gen5():
		plt.savefig("sync_gen5_w.png")
		plt.close()

		plt.figure()
		plt.plot(t_list, delta_list, label="Delta angle")
		plt.xlabel("Time [s]")
		plt.ylabel("Radians")
		plt.grid(ls="--", alpha=0.5)
		plt.legend()
		plt.savefig("sync_gen5_delta.png")
		plt.close()

		plt.figure()
		plt.plot(t_list, I_d_list, label="I_d")
		plt.plot(t_list, I_q_list, label="I_q")
		plt.xlabel("Time [s]")
		plt.ylabel("Current")
		plt.grid(ls="--", alpha=0.5)
		plt.legend()
		plt.savefig("sync_gen5_Iqd.png")
		plt.close()


		if __name__ == "__main__":
		test_gen5()

Admin message