<jupyter_start><jupyter_text>SMIB system as in Milano's book example 8.1<jupyter_code>%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as sopt
import ipywidgets
from pydae import ssa
import json<jupyter_output><empty_output><jupyter_text>Import system module<jupyter_code>from smib_milano_ex8p1_4ord_avr_pss import smib_milano_ex8p1_4ord_avr_pss_class<jupyter_output><empty_output><jupyter_text>Instantiate system<jupyter_code>smib = smib_milano_ex8p1_4ord_avr_pss_class()

xy_0_dict = {
'omega':1,'v_ref':1,'v_c':1
}<jupyter_output><empty_output><jupyter_text>Initialize the system (backward and foreward)<jupyter_code>events=[{'p_m':0.8, 'v_t':1.0, 'K_a':100, 'T_e':0.1}]
smib.initialize(events,xy_0_dict)

smib.save_0()
smib.ss()
smib.report_u()
smib.report_x()
smib.report_y()

ssa.eval_A(smib)
ssa.damp_report(smib)
smib = smib_milano_ex8p1_4ord_avr_pss_class()
smib.load_0('xy_0.json')
smib.ss()
smib.eval_jacobians()
smib.eval_A()
ssa.damp_report(smib)
smib.report_params()
smib.load_0('xy_0.json')
def obj(x):
    T_1 = x[0]
    #K_stab = x[1]
    smib.set_value('T_1',T_1)
    freq = 1.2
    T_2 = 0.1
    cplx = (1j*2*np.pi*freq*T_1 + 1)/(1j*2*np.pi*freq*T_2 + 1)

    smib.set_value('K_stab',1)
    smib.set_value('K_a',100)
    smib.set_value('H',6)
    smib.ss()
    smib.eval_jacobians()
    ssa.eval_A(smib)
    eig_values,eig_vectors = np.linalg.eig(smib.A)
    zetas = -eig_values.real/np.abs(eig_values)
    return -100*np.min(zetas)
sopt.differential_evolution(obj,bounds=[(0.1,5)])

events=[{'t_end':1.0},
        {'t_end':5.0, 'v_ref':v_ref_0*1.05}
       ]

smib.simulate(events,xy0='prev');

plt.close('all')
fig, axes = plt.subplots(nrows=2,ncols=2, figsize=(10, 5), frameon=False, dpi=70)

axes[0,0].plot(smib.T, smib.get_values('omega'), label=f'$\omega$')
axes[0,1].plot(smib.T, smib.get_values('v_t'),   label=f'$v_t$')
axes[1,0].plot(smib.T, smib.get_values('p_t'),   label=f'$p_t$')
axes[1,1].plot(smib.T, smib.get_values('q_t'),   label=f'$q_t$')

for ax in axes.flatten():
    ax.legend()<jupyter_output><empty_output><jupyter_text>Simulation<jupyter_code>smib = smib_milano_ex8p1_4ord_avr_pss_class()


events=[{'p_t':0.8, 'v_t':1.0, 'K_a':200, 'T_e':0.2, 'H':6, 'K_stab':0, 'T_1':0.1}]
smib.initialize(events,xy0=1)
 
v_ref_0 = smib.get_value('v_ref')
events=[{'t_end':1.0},
        {'t_end':15.0, 'v_ref':v_ref_0*1.05}
       ]

smib.simulate(events,xy0='prev');
plt.close('all')
fig, axes = plt.subplots(nrows=2,ncols=2, figsize=(10, 5), frameon=False, dpi=70)

axes[0,0].plot(smib.T, smib.get_values('omega'), label=f'$\omega$')
axes[0,1].plot(smib.T, smib.get_values('v_t'),   label=f'$v_t$')
axes[1,0].plot(smib.T, smib.get_values('p_t'),   label=f'$p_t$')
axes[1,1].plot(smib.T, smib.get_values('q_t'),   label=f'$q_t$')

for ax in axes.flatten():
    ax.legend()
smib = smib_milano_ex8p1_4ord_avr_pss_class()


events=[{'p_t':0.8, 'v_t':1.0, 'K_a':200, 'T_e':0.2, 'H':6, 'K_stab':0, 'T_1':0.1}]
smib.initialize(events,xy_0_dict)
 
ssa.eval_A(smib)
ssa.damp_report(smib)<jupyter_output><empty_output><jupyter_text>Run in two time intervals<jupyter_code>events=[{'t_end':1.0}]
syst.run(events)
events=[{'t_end':2.0}]
syst.run(events)
syst.get_value('omega')
events=[{'p_t':0.8, 'v_t':1.0, 'K_a':100, 'T_e':0.2, 'H':6, 'K_stab':0, 'T_1':0.1}]
smib.initialize(events,xy0=1)
ssa.eval_A(smib)
ssa.damp_report(smib)
ssa.participation(smib).abs().round(2)
smib.report_params()
Ts_control = 0.010
times = np.arange(0.0,10,Ts_control)

# Calculate second references
events=[{'P_t':0.9, 'Q_t':0.0}]
syst.initialize(events,xy0=1.0)
x_ref = np.copy(syst.struct[0].x)
v_f_ref = syst.struct[0]['v_f']
p_m_ref = syst.struct[0]['p_m']

# Calculate initial references
events=[{'P_t':0.0, 'Q_t':0.0}]
syst.initialize(events,xy0=1.0)
x_0 = np.copy(syst.struct[0].x)
v_f_0 = syst.get_value('v_f')
p_m_0 = syst.get_value('p_m')

# Control design
ssa.eval_ss(syst)
Q = np.eye(syst.N_x)*100
R = np.eye(syst.N_u)

K = ctrl.place(syst.A,syst.B,[-2.0+1j*6,-2.0-1j*6,-100,-101])
K,S,E = ctrl.lqr(syst.A,syst.B,Q,R)
Ad,Bd = ssa.discretise_time(syst.A,syst.B,Ts_control)
Kd,S,E = ssa.dlqr(Ad,Bd,Q,R)

for t in times:
    
    x = np.copy(syst.struct[0].x)
    v_f = v_f_0 
    p_m = p_m_0 

    if t>1.0: 
        u_ctrl = K*(x_ref - x)
        p_m = p_m_ref + u_ctrl[0]
        v_f = v_f_ref + u_ctrl[1]
        
        
    events=[{'t_end':t,'v_f':v_f,'p_m':p_m}]
    syst.run(events)

syst.post();
    
plt.close('all')
fig, axes = plt.subplots(nrows=2,ncols=2, figsize=(10, 5), frameon=False, dpi=50)

axes[0,0].plot(syst.T, syst.get_values('omega'), label=f'$\omega$')
axes[0,1].plot(syst.T, syst.get_values('v_1'),   label=f'$v_1$')
axes[1,0].plot(syst.T, syst.get_values('P_t'),   label=f'$P_t$')
axes[1,1].plot(syst.T, syst.get_values('Q_t'),   label=f'$Q_t$')
ssa.eval_ss(syst)

from scipy.signal import ss2tf,lti,bode
num,den =ss2tf(syst.A,syst.B,syst.C,syst.D,input=0)
G = lti(num[1],den)
w, mag, phase = G.bode()
plt.figure()
plt.semilogx(w, mag)    # Bode magnitude plot
plt.figure()
plt.semilogx(w, phase)  # Bode phase plot
plt.show()
events=[{'t_end':1.0,'P_t':0.8, 'Q_t':0.5},
        {'t_end':10.0, 'p_m':0.9}
       ]

syst.simulate(events,xy0=1.0);
syst.inputs_run_list

0.01/6
syst.B
syst.struct[0]['Fu']<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2101299190359631, 'fraction_numerical': 0.0495198644323103, 'mean_word_length': 5.493887530562347, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 0, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '612301', 'license': 'MIT', 'n_tokens_mistral': 2784, 'n_tokens_neox': 2680, 'n_words': 315, 'path': 'examples/grids/smib_milano_ex8p1/smib_milano_ex8p1_4ord_avr/smib_milano_ex8p1_4ord_avr_pss.ipynb', 'repo_name': 'pydae/pydae', 'size': '5311'}