Diseño del sistema de control del pédulo invertido#

Seguimos con el ejemplo de teoría, pero con un pequeño cambio:

w0=2
A=[[0,1],[-w0**2,0]]
B=[[0],[1]]
C=[1,0]
D=0

pendulo=ctrl.ss(A,B,C,D)
pendulo
\[\begin{split} \left(\begin{array}{rllrll|rll} 0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ -4\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ \hline 1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ \end{array}\right) \end{split}\]

Obtención de la ley de control#

K=ctrl.acker(pendulo.A,pendulo.B,[-2*w0,-2*w0])
K

array([[12.,  8.]])

Diseño del estimador completo#

L=ctrl.acker(pendulo.A.T,pendulo.C.T,[-4*w0,-4*w0]).T

L

array([[16.],
       [60.]])

np.linalg.eigvals(pendulo.A-pendulo.B@K)

array([-4., -4.])

np.linalg.eigvals(pendulo.A-L@pendulo.C)

array([-8., -8.])

Definimos el sistema compensador ( entrada \(y\), salida \(u\))#

s=pendulo
Acomp=s.A-s.B@K-L@s.C
Bcomp=L
Ccomp=-K
comp=ctrl.ss(Acomp,Bcomp,Ccomp,0)
comp
\[\begin{split} \left(\begin{array}{rllrll|rll} -16\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&16\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ -76\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-8\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&60\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ \hline -12\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-8\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ \end{array}\right) \end{split}\]
np.linalg.eigvals(comp.A)

array([-12.+7.74596669j, -12.-7.74596669j])

Vemos que los polos del compensador son estables!

 comp.zero()

array([-0.30952381+0.j])

Verificamos que los polos a lazo cerrado del sistema compensado son los adecuados:

sys_comp = ctrl.feedback(-s*comp)
sys_comp.pole()

array([-8.00000055+0.00000000e+00j, -7.99999945+0.00000000e+00j,
       -4.        +2.82767841e-07j, -4.        -2.82767841e-07j])

Ahora verificamos que el sistema a lazo cerrado tiene efectivamente los polos donde se pretendía. Igualmente será necesario rediseñar para lograr un controlador que sea correcto.

Sistema completo.#

sys_completo =  ctrl.connect(ctrl.append(s,comp),[[1,2],[2,1]],[1],[1,2])
t,y = ctrl.initial_response(sys_completo, X0=[1,0,0,0], T=np.linspace(0,3,1001))
Hide code cell source 
fig, ax = plt.subplots(2,1,figsize=(12,8))
ax[0].plot(t,y[0,:],label="sys completo")
ax[0].set_title('Respuesta al escalón en la referencia')
ax[0].set_xlabel("Tiempo")
ax[0].set_ylabel("Salida")
ax[0].grid()
ax[1].plot(t,y[1,:],label="sys completo")
ax[1].set_title('Respuesta al escalón en la referencia')
ax[1].set_xlabel("Tiempo")
ax[1].set_ylabel("Salida")
ax[1].grid()
fig.tight_layout()
../../../_images/85d113fa20115a792068b73888b22352e2fbb1cd6c31b039268e35a0b60de0e5.png 
np.linalg.eigvals(sys_completo.A)

array([-8.+6.17544035e-07j, -8.-6.17544035e-07j, -4.+1.20399686e-07j,
       -4.-1.20399686e-07j])

Ubicación óptima de los autovalores del estimador#

B1=B

def conjugate_tf(G):
    num = ctrl.tf(G).num[0][0]
    den = ctrl.tf(G).den[0][0]
    nume = [num[i]*((-1)**(len(num)%2+1-i)) for i in range(0, len(num))]
    dene = [den[i]*((-1)**(len(den)%2+1-i)) for i in range(0, len(den))]
    return ctrl.tf(nume,dene)

aux1 = ctrl.ss(A,B1,C,D)
aux2 = conjugate_tf(aux1)
G_srle = ctrl.tf(aux1)*aux2
ctrl.rlocus(G_srle, xlim=[-8,8], ylim=[-10,10]); #aumento los límites para ver la zona donde quiero ubicar el polo
plt.gcf().set_size_inches(8,6)
../../../_images/1f4473a05d516a5e1a480c3c422a57b6a2410672c56e8d9d1e310919444c56a0.png 
k=12000 # lo busco haciendo clcik sobre las lineas del root locus simétrico

r,k = ctrl.rlocus(G_srle, kvect=[8000],Plot=False)
r

/usr/share/miniconda/lib/python3.11/site-packages/control/rlocus.py:132: FutureWarning: 'Plot' keyword is deprecated in root_locus; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in root_locus; "

array([[-6.53957633-6.83857138j, -6.53957633+6.83857138j,
         6.53957633-6.83857138j,  6.53957633+6.83857138j]])

rsel = r[np.real(r)<0]
rsel

array([-6.53957633-6.83857138j, -6.53957633+6.83857138j])

L = ctrl.place(pendulo.A.T, pendulo.C.T, rsel).T
L

array([[13.07915266],
       [85.53211714]])

comp_srl = ctrl.ss(s.A-s.B@K-L@s.C, L, -K, 0)
comp_srl.pole()

array([-10.53957633+9.75103426j, -10.53957633-9.75103426j])

sys_completo_srl =  ctrl.connect(ctrl.append(s,comp_srl),[[1,2],[2,1]],[1],[1,2])
t_srl, y_srl = ctrl.initial_response(sys_completo_srl, X0=[1,0,0,0], T=np.linspace(0,3,1001))

L

array([[13.07915266],
       [85.53211714]])

K

array([[12.,  8.]])
Hide code cell source 
fig, ax = plt.subplots(2,1,figsize=(12,8))
ax[0].plot(t,y[0,:],label="sys completo 2do orden dominante")
ax[0].plot(t_srl,y_srl[0,:],label="sys completo SRL")
ax[0].set_title('Respuesta al escalón en la referencia')
ax[0].set_xlabel("Tiempo")
ax[0].set_ylabel("Salida")
ax[0].grid()
ax[1].plot(t,y[1,:],label="sys completo 2do orden dominante")
ax[1].plot(t_srl,y_srl[1,:],label="sys completo SRL")
ax[1].set_title('Respuesta al escalón en la referencia')
ax[1].set_xlabel("Tiempo")
ax[1].set_ylabel("Salida")
ax[1].grid()
fig.tight_layout()
../../../_images/1d3a38bf0a4cbb62f4742de23a0b38c48d2da8de33c5dd54a7c40d43fc05915f.png