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]:
G=zpk([],[-20 -10 10],-40)

A=[-20 100 2000
    1 0 0
    0 1 0];
B=[1;0;0];
C=[0 0 -40];
D=0;
sys=ss(A,B,C,D)


G =
 
          -40
  --------------------
  (s+20) (s+10) (s-10)
 
Continuous-time zero/pole/gain model.


sys =
 
  a = 
         x1    x2    x3
   x1   -20   100  2000
   x2     1     0     0
   x3     0     1     0
 
  b = 
       u1
   x1   1
   x2   0
   x3   0
 
  c = 
        x1   x2   x3
   y1    0    0  -40
 
  d = 
       u1
   y1   0
 
Continuous-time state-space model.

Verificamos controlabilidad y observabilidad.


In [2]:
MatCtr=[B A*B A^2*B]
MatObs=[C
    C*A
    C*A^2]
rank(MatCtr)
rank(MatObs)


MatCtr =

     1   -20   500
     0     1   -20
     0     0     1


MatObs =

     0     0   -40
     0   -40     0
   -40     0     0


ans =

     3


ans =

     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 [4]:
pol=[-10+20j,-10-20j,-40]
K1=place(A,B,[-10+20j,-10-20j,-40])
pol=[-10+20j,-10-20j,-20]
K2=acker(A,B,pol)
K2=place(A,B,pol)


pol =

 -10.0000 +20.0000i -10.0000 -20.0000i -40.0000          


K1 =

   1.0e+04 *

    0.0040    0.1400    2.2000


pol =

 -10.0000 +20.0000i -10.0000 -20.0000i -20.0000          


K2 =

          20        1000       12000


K2 =

   1.0e+04 *

    0.0020    0.1000    1.2000

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 roden 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=tf(G);
Gtfms=zpk([],[20 10 -10], 40);
%rlocus(Gtf*Gtfms)
p=rlocus(Gtf*Gtfms,38700);
pol3=p(1:3)
K3=place(A,B,pol3)


pol3 =

 -24.9963          
 -12.4981 +12.9850i
 -12.4981 -12.9850i


K3 =

   1.0e+04 *

    0.0030    0.1050    1.0119

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:

  1. 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
  2. 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:

  1. Ubicación de los polos muy rápidos comparados con la dinámica de la planta deseada.
  2. Ubicación de los polos algo más rápida que la dinímica de la planta deseada.
  3. 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=place(A',C',polObs1)'


polObs1 =

   -14   -15   -16


L1 =

  -59.5000
   -6.8500
   -0.6250

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


In [7]:
polObs2=[-27,-30,-33]
L2=place(A',C',polObs2)'


polObs2 =

   -27   -30   -33


L2 =

 -197.7500
  -34.7750
   -1.7500

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


In [13]:
B1=[10 0 0]';
sysGe=ss(A,B1,C,D);
Ge=tf(sysGe)
Gesim=tf(-400,[-1 20 100 -2000])
%rlocus(Ge*Gesim)
polObs3=rlocus(Ge*Gesim,10);
polObs3=polObs3(1:3)
L3=place(A',C',polObs3)'


Ge =
 
             -400
  ---------------------------
  s^3 + 20 s^2 - 100 s - 2000
 
Continuous-time transfer function.


Gesim =
 
              400
  ---------------------------
  s^3 - 20 s^2 - 100 s + 2000
 
Continuous-time transfer function.


polObs3 =

 -20.3966          
 -10.1983 + 3.4664i
 -10.1983 - 3.4664i


L3 =

  -53.0546
   -5.4045
   -0.5198

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 [21]:
T=[0 0 -40;1 0 0;0 1 0]
An=T*A*inv(T)
Bn=T*B
Cn=C*inv(T)
Dn=D


T =

     0     0   -40
     1     0     0
     0     1     0


An =

     0     0   -40
   -50   -20   100
     0     1     0


Bn =

     0
     1
     0


Cn =

     1     0     0


Dn =

     0

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 [16]:
Abb=An([2:end],[2:end]);
Aba=An([2:end],1);
Aab=An(1,[2:end]);
Aaa=An(1,1);

Ba=Bn(1);
Bb=Bn(2:end);

Lred1=place(Abb',Aab',polObs3(2:3))';

aux=[An Bn
    Cn Dn];

Nn=inv(aux)*[zeros(3,1);1];
Nxn=Nn(1:3);
Nun=Nn(4);

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


In [17]:
K3n=place(An,Bn,pol3)


K3n =

   1.0e+03 *

   -0.2530    0.0300    1.0496

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 [22]:
aux=[A B
    C D];
N=inv(aux)*[zeros(3,1);1];
Nx=N(1:3)
Nu=N(4)
Nvib=Nu+K3*Nx

Acomp=A-B*K3-L3*C;
Bcomp=[B*Nvib L3];
Ccomp=-K3;
Dcomp=[Nvib 0];
comp=ss(Acomp,Bcomp,Ccomp,Dcomp)


Nx =

         0
         0
   -0.0250


Nu =

   50.0000


Nvib =

 -202.9778


comp =
 
  a = 
               x1          x2          x3
   x1      -49.99      -949.6  -1.024e+04
   x2           1           0      -216.2
   x3           0           1      -20.79
 
  b = 
            u1       u2
   x1     -203   -53.05
   x2        0   -5.404
   x3        0  -0.5198
 
  c = 
               x1          x2          x3
   y1      -29.99       -1050  -1.012e+04
 
  d = 
         u1    u2
   y1  -203     0
 
Continuous-time state-space model.

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 [23]:
tf(comp)


ans =
 
  From input 1 to output:
  -203 s^3 - 8280 s^2 - 1.08e05 s - 4.803e05
  ------------------------------------------
     s^3 + 70.79 s^2 + 2205 s + 4.079e04
 
  From input 2 to output:
  1.252e04 s^2 + 3.764e05 s + 2.52e06
  -----------------------------------
  s^3 + 70.79 s^2 + 2205 s + 4.079e04
 
Continuous-time transfer function.

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 [24]:
Aaug=[zeros(1,1) C
    zeros(3,1) A]
Baug=[zeros(1,1)
    B]
Caug=[0 C]
pol4=[-5 pol3']

K4=place(Aaug,Baug,pol4)


Aaug =

           0           0           0         -40
           0         -20         100        2000
           0           1           0           0
           0           0           1           0


Baug =

     0
     1
     0
     0


Caug =

     0     0     0   -40


pol4 =

  -5.0000           -24.9963           -12.4981 -12.9850i -12.4981 +12.9850i


K4 =

   1.0e+04 *

   -0.1015    0.0035    0.1300    1.4867