Selección de Polos del Controlador

Para selecccionar del sistema a lazo cerrado debemos tener en cuenta los requerimientos de control.

Para cumplir tenemos dos opciones:

Si los requeriemientos son solo de performance dinámica a lazo cerrado podemos pensar al sistema a lazo cerrado como dominante de segundo orden, cumplir lo que se nos pide de requerimientos y ubicar dos polos dominantes que coumplan con lo requerido y un tercero suficientemente alejado.
Realizar un root locus y implicitamente ubicar los tres polos de forma de que el sistema cumpla de la mejor manera posible con los requerimientos pero además teniendo en cuenta el esfuerzo de control

De los requerimientos del problema surge que la posición de los polos deseada para un segundo orden puro es: $10\pm j20$.

Entonces primero definimos el sistema:



In [1]:

    
import numpy as np
import control as ctrl
import matplotlib.pylab 
%matplotlib inline
A=np.matrix("-20 100 2000;1 0 0;0 1 0")
B=np.matrix("1;0;0")
C=np.matrix("0 0 -40")
D=np.matrix("0");
sys=ctrl.ss(A,B,C,D)

Verificamos controlabilidad y observabilidad.



In [2]:

    
MatCtr=np.concatenate((B,A*B,np.linalg.matrix_power(A,2)*B),axis=1)
MatObs=np.concatenate((C,C*A,C*np.linalg.matrix_power(A,2)),axis=0)
print MatCtr
print np.linalg.matrix_rank(MatCtr)
print MatObs
print np.linalg.matrix_rank(MatObs)









    



[[  1 -20 500]
 [  0   1 -20]
 [  0   0   1]]
3
[[  0   0 -40]
 [  0 -40   0]
 [-40   0   0]]
3

Podemos ver que el sismtea es controlable y observable. Es era de esperar ya que no modelamos ninguna cancelación de polos y ceros en la transferencia a partir de la cual obtuvimos el esapcio de estado.



In [3]:

    
pol1=[-40.,-10-20j,-10+20j]
K1=ctrl.acker(A,B,pol1)
print "Ubicando los polos en", np.linalg.eigvals(A-B*K1),"se tiene una ganancia ", K1
pol2=[-20.,-10-20j,-10+20j]
K2=ctrl.acker(A,B,pol2)
print "Ubicando los polos en", np.linalg.eigvals(A-B*K2),"se tiene una ganancia ", K2









    



Ubicando los polos en [-40. +0.j -10.+20.j -10.-20.j] se tiene una ganancia  [[    40.   1400.  22000.]]
Ubicando los polos en [-10.+20.j -10.-20.j -20. +0.j] se tiene una ganancia  [[    20.   1000.  12000.]]

Vemos que las ganancias son más grandes cuando tratamos de alejar mucho el tercer polo. La ventaja que se tiene es que haciendo esto el sistema se comporta de forma más parecida a un segundo orden puro. El precio que se paga es un mayor esfuerzo de control.

Si hacemos una optimización usando el lugar de las raíces simétrico se puede obtener una ubicación de polo que cumple de forma aproximada el requerimiento de control y bajan bastante el esfuerzo de control.



In [4]:

    
Gtf=ctrl.tf(sys);
Gtfms=ctrl.tf(-40,[-1,20,100,-2000])
print Gtf, Gtfms
rlist,klist=ctrl.rlocus(Gtf*Gtfms,np.logspace(-3,5,1000),Plot=True,PrintGain=False)









    



            -40
---------------------------
s^3 + 20 s^2 - 100 s - 2000
 
             -40
----------------------------
-s^3 + 20 s^2 + 100 s - 2000



In [5]:

    
rlist,klist=ctrl.rlocus(Gtf*Gtfms,klist=[38700],Plot=False)
pol3 = rlist.reshape(6)[0:3]
print pol3.size
print pol3
K3=ctrl.acker(A,B,pol3)
print np.linalg.eigvals(A-B*K3)









    



3
[-24.99626172 +0.j         -12.49813086-12.98498459j
 -12.49813086+12.98498459j]
[-24.99626172 +0.j         -12.49813086+12.98498459j
 -12.49813086-12.98498459j]

Finalemte utilizamos esta última selección de polos para implementación ya que además aseguran estabilidad y buenos margenes de ganancia y fase .

Seleccion polos observador

Para esto también existen dos criterios:

Ubicar los polos del observador usfientemente más rápido que los polos dominantes del controlador para que sean estos últimos los que gobiernen la dinámica
Tienendo en cuenta un criterio de optimización en el cual se optimize la relación entre el ruido en los sensores a partir de los cuales se hace el estimar y el rechazo a perturbaciones.

Para esta clase hicimos observadores con 3 ubicaciones diferentres:

Ubicación de los polos muy rápidos comparados con la dinámica de la planta deseada.
Ubicación de los polos algo más rápida que la dinímica de la planta deseada.
Mediante la optimización usando el lugar de las raices simétrico.

Para esto se definio una perturbación al sistema.

Definimos los polos del observador el la ubicación algo más rápido.



In [6]:

    
polObs1=[-14,-15,-16]
L1=ctrl.acker(A.T,C.T,polObs1).T
print L1









    



[[-59.5  ]
 [ -6.85 ]
 [ -0.625]]

Definimos los polos del observador el la ubicación mucho más rápido.



In [7]:

    
polObs2=[-27,-30,-33]
L2=ctrl.acker(A.T,C.T,polObs2).T

Definimos la ubicación de los polo meidante optimización:



In [8]:

    
B1="10;0;0"
sysGe=ctrl.ss(A,B1,C,D);
Ge=ctrl.tf(sysGe)
print Ge









    



            -400
---------------------------
s^3 + 20 s^2 - 100 s - 2000



In [9]:

    
Gesim=ctrl.tf(-400,[-1,20, 100, -2000])
ctrl.rlocus(Ge*Gesim,Plot=True)
rlist,klist=ctrl.rlocus(Ge*Gesim,[10],Plot=False)
polObs3=rlist.reshape(6)[0:3]



In [10]:

    
L3=ctrl.acker(A.T,C.T,polObs3).T
print L3









    



[[-53.05463891]
 [ -5.40445585]
 [ -0.51982958]]

Observador Reducido

Como la salida del sismtema no es la primera variable de estado, sino una propocionalidad de la tercera, proponemos una transformación que ubique la variable de estado tercera en la primera y su valor sea igual a la salida.

Esta transformación se obtiene haciendo:



In [11]:

    
T=np.matrix("0 0 -40;1 0 0;0 1 0")
An=T*A*np.linalg.inv(T)
Bn=T*B
Cn=C*np.linalg.inv(T)
Dn=D

Estas matrices las que usaremos para calcular el estimador. La planta no se tocará para nada ya que no es parte del diseño.



In [12]:

    
Abb=An[1:3,1:3];
Aba=An[1:3,0];
Aab=An[0,1:3];
Aaa=An[0,0];

Ba=Bn[0];
Bb=Bn[1:3];

polRed1=polObs3[1:3]

Lred1=ctrl.acker(Abb.T,Aab.T,polRed1).T;

aux1=np.concatenate((np.concatenate((An,Bn),axis=1),np.concatenate((Cn,Dn),axis=1)),axis=0);
aux2=np.concatenate((np.zeros((3,1)),np.matrix([1])),axis=0)
Nn=np.linalg.inv(aux1)*aux2;
Nxn=Nn[0:3];
Nun=Nn[3];

Recalculo la ganancia del controlador para el obwervador reducido para que la salida tenga el valor de ganancia requerido (escaleado).



In [13]:

    
K3n=ctrl.acker(An,Bn,pol3)
print K3n









    



[[ -252.9778313     29.99252343  1049.62619963]]

Notar que las ganacias son las mismas que K3, salvo que la tercera ganacia aparece ahora en primer lugar cambiada de signoy dividida por 40 que conicide con el escaleo y reodenamiento de eatados que hicimos para poder usar la teor'ia y hacer el estimador reducido.

Implementación del controlador obtenido por Ecuanciones de estado

Definimos el controlador en forma de ecuaciones de estado



In [14]:

    
aux3=np.concatenate((np.concatenate((A,B),axis=1),np.concatenate((C,D),axis=1)),axis=0);
N= np.linalg.inv(aux3)*aux2;
Nx=N[0:3]
Nu=N[3]
Nvib=Nu+K3*Nx

Acomp=A-B*K3-L3*C;
Bcomp=np.concatenate((B*Nvib,L3),axis=1);
Ccomp=-K3;
Dcomp=np.concatenate((Nvib,np.zeros((1,1))),axis=1);
comp=ctrl.ss(Acomp,Bcomp,Ccomp,Dcomp)
print Nx
print Nu









    



[[ -3.55271368e-16]
 [  0.00000000e+00]
 [ -2.50000000e-02]]
[[ 50.]]

Podemos ver que este compensador no es exatamente el compensador que habíamos obtenido en clase. Este compensado, así definido, tiene 2 entradas:

la entrada de la referencia (la primera)
la salida del modelo de la planta (la segunda)

En clase generamos una función transferencia que solo tenía como entrada la salida de la planta. Ya que el compensador que usamos no está pensado para seguimiento a referencias, sino solo para la ubicar su dinámica a lazo cerrado.

La forma en que se agrega la referencia al controlador es teniendola en cuenta. Para esto se debe tomar la ecuación (7.191a) y reemplazarla en (7.191b) del Franfklyn 6ta edición. Podemos ver que se desarrollar ese termino queda parecido a lo que hicimos en clase pero con un término extra debido a la entrada de la referencia. Este término es $N_xBr$. Para la salida del controlador, debemos ver la figura 7.49b, donde se ve claramente que la referencia entra directamente a la planta múltiplicada por $\overline N$.Esto sería la parte de la matriz $D$ para la ecuación de estado del controlador. Además se tiene que sumar $K\hat x$ que sería la matriz C de compensador.

Podemos ver que se hacemos $r=0$ las ecuaciones de estado nuevas quedan igual que las que hicimos en clase.

Por otro lado, el análisis de los diagramas de Bode realizados en clases están bien hechos, ya que la dinamica de la planta a lazo cerrado (o de lazo), es la misma si uno considera una entrada o las dos.

Podemos ver claramente esto transformando este sistema en espacio de estado en funciones transferencia (en este caso dos funciones transferencias una para la transferencia de $r$ hacia $u$ y la otra para la transferencia de $y$ hacia $u$. Está útima es la que habiamso obtenido.



In [15]:

    
ctrl.tf(comp)









    Out[15]:





Input 1 to output 1:
-203 s^3 - 8280 s^2 - 1.08e+05 s - 4.803e+05
--------------------------------------------
    s^3 + 70.79 s^2 + 2205 s + 4.079e+04

Input 2 to output 1:
1.252e+04 s^2 + 3.764e+05 s + 2.52e+06
--------------------------------------
 s^3 + 70.79 s^2 + 2205 s + 4.079e+04

Control Integral Para Seguimiento Robusto

Para agregar el contro integral podemos usar cualquier observador de estados y aumentar la planta para que este prensente el estado del integrador y darle una posici'on al polo ese.



In [16]:

    
Aaug=np.concatenate((np.concatenate((np.zeros((1,1)),C),axis=1),
                     np.concatenate((np.zeros((3,1)),A),axis=1)),axis=0)
Baug=np.concatenate((np.zeros((1,1)),B),axis=0)
Caug=np.concatenate((np.zeros((1,1)),C),axis=1)
pol4=np.concatenate((np.array([-5.]), pol3),axis=0)
K4=ctrl.acker(Aaug,Baug,pol4)
print K4









    



[[ -1014.88915651     34.99252343   1299.5888168   14867.24425024]]