Commit fa5b923c authored by William Gurecky's avatar William Gurecky
Browse files

adds lin solve for the i_qd armeture currents

parent 3f3cee65

dthydro_cfd/generator_models.py

+97 −84
Original line number Diff line number Diff line
@@ -2,6 +2,8 @@ 
Implements sychronous generator models from:



Machowski. Power System Dynamics: Stability and Control. 2008.



Author: William Gurecky, wll@ornl.gov

"""

import numpy as np

from copy import deepcopy
@@ -180,8 +182,8 @@ class SyncGen5(object): 
    @property

    def power_e(self):

        """ electrical power W """

        I_d = self.I_d

        I_q = self.I_q

        I_d = self.I_d()

        I_q = self.I_q()

        X_d_pp = self.machine_params["Xd2"]

        X_q_pp = self.machine_params["Xq2"]

        P_e = (self.x[4]*I_d + self.x[3]*I_q) + (X_d_pp - X_q_pp)*I_d*I_q
@@ -206,7 +208,7 @@ class SyncGen5(object): 
    @property

    def delta_angle(self):

        """ delta phasor angle """

        return self.x[1] % (np.pi / 2)

        return self.x[1] # % (np.pi / 2)



    @property

    def E_q_p(self):
@@ -240,72 +242,32 @@ class SyncGen5(object): 
            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 """

    def I_qd(self, R=0.):

        """ Solve linear system for I_q, I_d """

        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)

        X_q_pp = self.machine_params["Xq2"]

        X_d_pp = self.machine_params["Xd2"]

        M_tr = np.array([

            [0., 1.],

            [1.,  0.]])

        Arm_tr = np.array([

            [R, X_q_pp],

            [-X_d_pp, R],

            ])

        V_dq = self.V_dq_from_vs()

        K_dq = np.array([

            [-1*V_dq[0] - self.E_d_pp],

            [V_dq[1] - self.E_q_pp]

            ])

        B_tr = Xs - M_tr @ Arm_tr

        I_dq = np.linalg.inv(B_tr) @ M_tr @ K_dq

        return np.array([I_dq[1], -I_dq[0]])



    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, R=0.):

        return self.I_qd(R)[0]



    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 I_d(self, R=0.):

        return self.I_qd(R)[1]



    def arm_v(self, R=0.0):

        """ armature_voltage """
@@ -315,21 +277,75 @@ class SyncGen5(object): 
        X_q = self.machine_params["Xq"]

        X_q_p = self.machine_params["Xq1"]

        X_q_pp = self.machine_params["Xq2"]

        I_d = float(self.I_d(R))

        I_q = float(self.I_q(R))

        return armature_voltage(

            R, X_q_pp, X_d_pp, self.E_d_pp,

            self.E_q_pp, self.I_d, self.I_q)

            self.E_q_pp, I_d, I_q)



    def V_gdq(self, R=0.):

        """ alias to arm_v """

        return self.arm_v(R)





    # 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 test_gen5():

    """ Test of the 5th order sychronous generator model """

    from scipy.integrate import solve_ivp

    # Driving voltage, generator main field voltage

    e_f = 10.0

    delta_w = 2.0 * np.pi * 1e-4  # rad / s

    e_f = 220.0

    delta_w = 2.0 * np.pi * 1e-14  # rad / s

    w_0 = W_S + delta_w

    # System voltage (13.5 kV)

    V_s_0 = 13.5e3

    # Mechanical power input

    P_m = 25.0e6  # W

    P_m = 10.0e6  # W

    # Driving voltage from excitation system

    E_f = e_f



    # Unit test parks transform

    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.

    I_q_0 = 9.5e3 * 1.
@@ -344,28 +360,24 @@ 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))



    # Init generator model

    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, E_f)



    # test time integration of the system

    dt = 1e-1

    dt = 0.5e-1

    x_new = x_0

    p_e_list = []

    t_list = []

    w_list = []

    delta_list = []

    I_q_list, I_d_list = [], []

    for i in range(300):

    for i in range(120):

        sol = solve_ivp(gen5, [t, t+dt], x_new, args=(P_m, E_f))

        x_new = sol.y[:, -1]

        # reset rpm
@@ -374,12 +386,13 @@ def test_gen5(): 
        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)

        I_d = float(gen5.I_d())

        I_q = float(gen5.I_q())



        I_q_list.append(I_q)

        I_d_list.append(I_d)

        t_list.append(t)



        aV = gen5.arm_v(R=0.0)

        V_0, V_d, V_q = 0., aV[0], aV[1]