```python
import numpy as np
import sympy as sym
import numba
import pydae.build as db
```


```python

```

\begin{eqnarray}
\dot \delta &=& \Omega_b \left(\omega - \omega_s\right)\\
\dot \omega &=& 1/(2 H) \left(p_m - p_e - D (\omega - \omega_s) \right)\\
\end{eqnarray}

$$ \sf
\Omega_{b} \left(\sf \omega - \omega_{s}\right)
$$

## System definition 


```python
params = {'X_d':1.81,'X1d':0.3,'T1d0':8.0,   # synnchronous machine d-axis parameters
          'X_q':1.76,'X1q':0.65,'T1q0':1.0,  # synnchronous machine q-axis parameters
          'R_a':0.003,'X_l': 0.05, 
          'H':6.0,'D':1.0,'S_n':1.0,
          'Omega_b':2*np.pi*50,'omega_s':1.0,
          'K_a':100, 'T_r':0.1, 'v_pss':0.0,
          'Droop':0.05,'T_m':5.0,'K_sec':0.001,
          'K_delta':0.01}


u_ini_dict = {'P_l':0.8, 'Q_l':0.2,'P_cig':0.1,'Q_cig':0.0,'v_ref':1.0}  # for the initialization problem
u_run_dict = {'P_l':0.8, 'Q_l':0.2,'P_cig':0.1,'Q_cig':0.0,'v_ref':1.0}  # for the running problem (here initialization and running problem are the same)


x_list = ['delta','omega','e1q','e1d','v_c','p_m','xi_m']    # [inductor current, PI integrator]
y_ini_list = ['i_d','i_q','v_1','theta_1','P_t','Q_t','v_f','p_m_ref','v_0','theta_0'] # for the initialization problem
y_run_list = ['i_d','i_q','v_1','theta_1','P_t','Q_t','v_f','p_m_ref','v_0','theta_0'] # for the running problem (here initialization and running problem are the same)

sys_vars = {'params':params,
            'u_list':u_run_dict,
            'x_list':x_list,
            'y_list':y_run_list}

exec(db.sym_gen_str())  # exec to generate the required symbolic varables and constants
```


```python
v_d = v_1*sin(delta - theta_1) 
v_q = v_1*cos(delta - theta_1) 

p_e = i_d*(v_d + R_a*i_d) + i_q*(v_q + R_a*i_q) 


ddelta = Omega_b*(omega - omega_s) - K_delta*delta
domega = 1/(2*H)*(p_m - p_e - D*(omega - omega_s))
de1q = 1/T1d0*(-e1q - (X_d - X1d)*i_d + v_f)
de1d = 1/T1q0*(-e1d + (X_q - X1q)*i_q)
dv_c =   (v_1 - v_c)/T_r
dp_m =   (p_m_ref - p_m)/T_m
dxi_m =   omega - 1

g_1  = v_q + R_a*i_q + X1d*i_d - e1q
g_2  = v_d + R_a*i_d - X1q*i_q - e1d
g_3  = P_t*S_n + P_cig - (v_1*v_0*sin(theta_1 - theta_0))/X_l
g_4  = Q_t*S_n + Q_cig + (v_1*v_0*cos(theta_1 - theta_0))/X_l - v_1**2/X_l
g_5  = i_d*v_d + i_q*v_q - P_t
g_6  = i_d*v_q - i_q*v_d - Q_t
g_7  = K_a*(v_ref - v_c + v_pss) - v_f 
g_8  = -p_m_ref - K_sec*xi_m - 1/Droop*(omega - 1)
g_9  = -P_l - (v_0*v_1*sin(theta_0 - theta_1))/X_l
g_10 = -Q_l + (v_0*v_1*cos(theta_0 - theta_1))/X_l - v_0**2/X_l

sys = {'name':'iso_milano_ex8p1_4ord_avr_gov',
       'params_dict':params,
       'f_list':[ddelta,domega,de1q,de1d,dv_c,dp_m,dxi_m],
       'g_list':[g_1,g_2,g_3,g_4,g_5,g_6,g_7,g_8,g_9,g_10],
       'x_list':x_list,
       'y_ini_list':y_ini_list,
       'y_run_list':y_run_list,
       'u_run_dict':u_run_dict,
       'u_ini_dict':u_ini_dict,
       'h_dict':{'p_m':p_m,'P_l':P_l,'Q_l':Q_l,'P_cig':P_cig,'Q_cig':Q_cig}}

sys = db.system(sys)
db.sys2num(sys)
                
```


```python

```


```python

```


```python

```
