<jupyter_start><jupyter_text><NAME>., & <NAME>. (2012). A model-based approach for small-signal stability assessment of unbalanced power systems. IEEE Transactions on Power Systems, 27(4), 2006–2014. https://doi.org/10.1109/TPWRS.2012.2192943<jupyter_code>psi_1d,psi_1q,psi_2q = sym.symbols('psi_1d,psi_1q,psi_2q ')  
e1d,e1q,e_fd = sym.symbols('e1d,e1q,e_fd')  
X_d,X_q,X1d,X1q,X2d,X2q,X_l,R_s = sym.symbols('X_d,X_q,X1d,X1q,X2d,X2q,X_l,R_s')  
T1d0,T1q0,T2d0,T2q0 = sym.symbols('T1d0,T1q0,T2d0,T2q0')  
i_a_re,i_b_re,i_c_re = sym.symbols('i_a_re,i_b_re,i_c_re') 
i_a_im,i_b_im,i_c_im = sym.symbols('i_a_im,i_b_im,i_c_im') 
i_d,i_q = sym.symbols('i_d,i_q') 
V_a,V_b,V_c,delta,theta_a,theta_b,theta_c = sym.symbols('V_a,V_b,V_c,delta,theta_a,theta_b,theta_c') 

pi = sym.pi
sin = sym.sin
cos = sym.cos 

de1q = 1/T1d0*(e_fd - e1q - (X_d-X1d)*(i_d-(X1d-X2d)/(X1d-X_l)**2)*(psi_1d+(X1d-X_l)*i_d-e1q))
de1d = 1/T1q0*(     - e1d + (X_q-X1q)*(i_q-(X1q-X2q)/(X1q-X_l)**2)*(psi_2q+(X1q-X_l)*i_q+e1q))
dpsi_1d = 1.0/T2d0*(-psi_1d + e1q - (X1d-X_l)*i_d)
dpsi_2q = 1.0/T2q0*(-psi_2q - e1d - (X1q-X_l)*i_q)
    
eq_e1d = e1d - 1/3*(V_a*sin(delta-theta_a)+  
                    V_b*sin(delta-theta_b-2*pi*3)+
                    V_c*sin(delta-theta_c-4*pi*3)) - R_s*i_d + X2q*i_q
eq_e1q = e1q - 1/3*(V_a*cos(delta-theta_a)+  
                    V_b*cos(delta-theta_b-2*pi*3)+
                    V_c*cos(delta-theta_c-4*pi*3)) - R_s*i_q - X2d*i_d
ip_re = 1/3*(i_a_re + cos(2*pi/3)*(i_b_re + i_c_re)-sin(2/3*pi)*(i_b_im - i_c_im)) 
ip_im = 1/3*(i_a_im + cos(2*pi/3)*(i_b_im + i_c_im)+sin(2/3*pi)*(i_b_re - i_c_re)) 

eq_id = i_d - (ip_re*cos(delta-pi/2) + ip_im*sin(delta-pi/2))
eq_iq = i_q - (ip_re*sin(delta-pi/2) - ip_im*cos(delta-pi/2))

#s = sym.solve([eq_ds,eq_qs,eq_dr,eq_qr],[i_ds,i_qs,i_dr,i_qr])

#for item in s:
#    print(item, '=', s[item])



d =dict(
H = 0.7,
X_d = 2.06,
X_q = 2.5,
X1d = 0.398,
X1q = 0.3,
X2d = 0.254,
X2q = 0.254,
X_l = 0.1,
    R_s = 0.01,
T1d0 = 7.8,
T1q0 = 3.0,
T2d0 = 0.066,
T2q0 = 0.075,
V_a = 1.0,
V_b = 1.0,
V_c = 1.0,
theta_a = 0.0,
theta_b = -2*np.pi/3,  
theta_c = -4*np.pi/3,
delta = 0.2,
e_fd = 1.1,
x=np.zeros((4,1)),
f=np.zeros((4,1)),
y=np.zeros((4,1)),
g=np.zeros((4,1))
)



struct = d2np(d)

for item in d:
    print("    {:s} = struct[i]['{:s}']".format(item,item) )<jupyter_output>H = struct[i]['H']
    X_d = struct[i]['X_d']
    X_q = struct[i]['X_q']
    X1d = struct[i]['X1d']
    X1q = struct[i]['X1q']
    X2d = struct[i]['X2d']
    X2q = struct[i]['X2q']
    X_l = struct[i]['X_l']
    R_s = struct[i]['R_s']
    T1d0 = struct[i]['T1d0']
    T1q0 = struct[i]['T1q0']
    T2d0 = struct[i]['T2d0']
    T2q0 = struct[i]['T2q0']
    V_a = struct[i]['V_a']
    V_b = struct[i]['V_b']
    V_c = struct[i]['V_c']
    theta_a = struct[i]['theta_a']
    theta_b = struct[i]['theta_b']
    theta_c = struct[i]['theta_c']
    delta = struct[i]['delta']
    e_fd = struct[i]['e_fd']
    x = struct[i]['x']
    f = struct[i]['f']
    y = struct[i]['y']
    g = struct[i]['g']<jupyter_text>Given voltages V_a,V_b,V_c and<jupyter_code><EMAIL>(nopython=True, cache=True)
def sm_ord6(struct,i=0):
    x_idx = 0
    y_idx = 0
    
    H = struct[i]['H']
    X_d = struct[i]['X_d']
    X_q = struct[i]['X_q']
    X1d = struct[i]['X1d']
    X1q = struct[i]['X1q']
    X2d = struct[i]['X2d']
    X2q = struct[i]['X2q']
    X_l = struct[i]['X_l']
    R_s = struct[i]['R_s']
    T1d0 = struct[i]['T1d0']
    T1q0 = struct[i]['T1q0']
    T2d0 = struct[i]['T2d0']
    T2q0 = struct[i]['T2q0']
    V_a = struct[i]['V_a']
    V_b = struct[i]['V_b']
    V_c = struct[i]['V_c']
    theta_a = struct[i]['theta_a']
    theta_b = struct[i]['theta_b']
    theta_c = struct[i]['theta_c']
    delta = struct[i]['delta']
    e_fd = struct[i]['e_fd']
    x = struct[i]['x']
    f = struct[i]['f']
    y = struct[i]['y']
    g = struct[i]['g']


    e1q =    struct[i]['x'][x_idx+0,0]
    e1d =    struct[i]['x'][x_idx+1,0]
    psi_1d = struct[i]['x'][x_idx+2,0]
    psi_2q = struct[i]['x'][x_idx+3,0]     
   
    i_d =   struct[i]['y'][y_idx+0,0]
    i_q =   struct[i]['y'][y_idx+1,0]
    ip_re = struct[i]['y'][y_idx+2,0]
    ip_im = struct[i]['y'][y_idx+3,0]   
   
    pi = np.pi
    sin = np.sin
    cos = np.cos 

    de1q = 1/T1d0*(e_fd - e1q - (X_d-X1d)*(i_d-(X1d-X2d)/(X1d-X_l)**2)*(psi_1d+(X1d-X_l)*i_d-e1q))
    de1d = 1/T1q0*(     - e1d + (X_q-X1q)*(i_q-(X1q-X2q)/(X1q-X_l)**2)*(psi_2q+(X1q-X_l)*i_q+e1d))
    dpsi_1d = 1.0/T2d0*(-psi_1d + e1q - (X1d-X_l)*i_d)
    dpsi_2q = 1.0/T2q0*(-psi_2q - e1d - (X1q-X_l)*i_q)

    eq_id = e1d - 1/3*(V_a*sin(delta-theta_a)+  
                        V_b*sin(delta-theta_b-2*pi/3)+
                        V_c*sin(delta-theta_c-4*pi/3)) - R_s*i_d + X2q*i_q
    eq_iq = e1q - 1/3*(V_a*cos(delta-theta_a)+  
                        V_b*cos(delta-theta_b-2*pi/3)+
                        V_c*cos(delta-theta_c-4*pi/3)) - R_s*i_q - X2d*i_d
    v_d = 0.0
    v_q = 1.0
    eq_id = e1d - v_d - R_s*i_d + X2q*i_q
    eq_iq = e1q - v_q- R_s*i_q - X2d*i_d
    #ip_re = 1/3*(i_a_re + cos(2*pi/3)*(i_b_re + i_c_re)-sin(2/3*pi)*(i_b_im - i_c_im)) 
    #ip_im = 1/3*(i_a_im + cos(2*pi/3)*(i_b_im + i_c_im)+sin(2/3*pi)*(i_b_re - i_c_re)) 

    eq_ip_re= i_d - (ip_re*cos(delta-pi/2) + ip_im*sin(delta-pi/2))
    eq_ip_im = i_q + (ip_re*sin(delta-pi/2) - ip_im*cos(delta-pi/2))
    
    struct[i]['f'][x_idx+0,0] = float(de1q)
    struct[i]['f'][x_idx+1,0] = float(de1d)
    struct[i]['f'][x_idx+2,0] = float(dpsi_1d)
    struct[i]['f'][x_idx+3,0] = float(dpsi_2q)

    struct[i]['g'][y_idx+0,0] = float(eq_id)
    struct[i]['g'][y_idx+1,0] = float(eq_id)
    struct[i]['g'][y_idx+2,0] = float(eq_ip_re)
    struct[i]['g'][y_idx+3,0] = float(eq_ip_im)

    return 0

class problem(object)
def problem(lam, i=0):
    #print(lam.shape)
    struct[i]['x'] = lam[0:4,:]
    struct[i]['y'] = lam[4:8,:]
    Phi = np.zeros((8,1))
    sm_ord6(struct,i=0)
    Phi[0:4,:] = struct[i]['f']
    Phi[4:8,:] = struct[i]['g']   
    return Phi
    
sm_ord6(struct,i=0)
i=0
lam_0 = np.zeros((8,1))
lam_0[0] = 1.0
struct[i]['delta'] = 0.2
struct[i]['e_fd'] = 1.2
problem(lam)
solution = sopt.broyden2(problem,lam_0,  maxiter=1000)
solution
struct[i]['x']
Omega_b = 2.0*np.pi*50.0
S_b = 1.0e6
U_b = 690.0
Z_b = U_b**2/S_b
#nu_w =np.linspace(0.1,15,N)
H = 0.001

# H = 0.5*J*Omega_t_n**2/S_b
S_b = 2.0e6
Omega_t_n = 1.5
J_t = 2*H*S_b/Omega_t_n**2


#Z_b = 1.0
#Omega_b = 1.0
d =dict(R_r = 0.01*Z_b,
        R_s = 0.01*Z_b,
        L_r = 0.08*Z_b/Omega_b,
        L_s = 0.1*Z_b/Omega_b,
        L_m = 3.0*Z_b/Omega_b,
        psi_ds = 0.0,
        psi_qs = 0.0,
        psi_dr = 0.0,
        psi_qr = 0.0,
        i_ds   = 0.0,
        i_qs   = 0.0,
        i_dr   = 0.0,
        i_qr   = 0.0,
        v_ds   = 0.0,
        v_qs   = 0.0,
        v_dr   = 0.0,
        v_qr   = 0.0,
        omega_r   = Omega_b*0.99,
        omega_s   = Omega_b,
        sigma   = 0.0,
        tau_e   = 0.0,        
        x = np.zeros((6,1)),
        f = np.zeros((6,1)),   
        Dt = 0.0,
        J_t = J_t,  
        omega_t = 0.0,  
        tau_t = 0.0,  
        tau_r = 0.0, 
        N_tr = 20.0, 
        dfim_idx = 0,
        mech_idx = 4,
        )

struct = d2np(d)
dfim(struct,0)
wecs_mech_1(struct,0)


struct = d2np(d)
struct['v_ds'] = 325.0
struct['v_qs'] = 0.0
struct['omega_r'] = Omega_b*1.0
Dt = 1.0e-4
struct[0]['x']= np.zeros((6,1))
struct[0]['x'][4,0] =   Omega_b*1.0
Tau_e = []
T =[]
N_steps = 100000
X = np.zeros((N_steps,6))

def f_eval(struct):
    dfim(struct,0)
    wecs_mech_1(struct,0)
    return struct[0]['f']
    
    
for it in range(N_steps):
    t = Dt*it
    if t>1.0:
        struct[0]['tau_t'] = 1.0e4
    dfim(struct,0)
    f1 = np.copy(f_eval(struct))
    x1 = np.copy(struct[0]['x'])
    struct[0]['x'] = np.copy(x1 + Dt*f1)
    dfim(struct,0)
    struct[0]['x'] = np.copy(x1 + 0.5*Dt*(f1 + f_eval(struct)))
    struct[0]['tau_r'] =  -struct[0]['tau_e']
    Tau_e += [float(struct['tau_e'])]
    X[it,:] = np.copy(struct[0]['x'][:]).T
    T +=[t]    
    
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(8, 5), sharex = True)
axes.plot(T,X[:,0])
axes.plot(T,X[:,1])
axes.plot(T,X[:,2])
axes.plot(T,X[:,3])

fig.savefig('dfim_tau_e.svg', bbox_inches='tight')

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(8, 5), sharex = True)
axes.plot(T,Tau_e)

fig.savefig('dfim_tau_e.svg', bbox_inches='tight')

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(8, 5), sharex = True)
axes.plot(T,X[:,4])
X[:,4]
20*15.68<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22727807172251616, 'fraction_numerical': 0.07078189300411522, 'mean_word_length': 2.857596371882086, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '239421', 'license': 'MIT', 'n_tokens_mistral': 4847, 'n_tokens_neox': 4612, 'n_words': 587, 'path': 'models/sm_ord6.ipynb', 'repo_name': 'pydgrid/pydgrid', 'size': '8519'}