Linealización y control SISO

Table of Contents
Proceso de linealización Control Función de transferencia SISO La forma básica de un sistema de control SISO La función de transferencia de un sistema de control Aplicación del método de Routh-Hurwitz Explicación del método Caso de un cero en la fila Caso de una fila de ceros Aplicación del método al caso del péndulo Control de un sistema discreto Equivalencia de Z y S Transformación Bilineal

Proceso de linealización

Primero generamos una tabla para controlar nuestro sistema de simulink. Este sistema debe estar en la misma forma que utilizamos para realizar el análisis de sensibilidad
clear;clc
model='PenduloC';% nombre del modelo
ParIn={'g','l','m','f','x10','x20', 'u'}; %nombres de los parámetros (tal y como aparecen en la máscara del sistema)
Ranges=[9.5 10.5; 2 4; 3 7; 0 1; -1 1;-1 1; -20 20]; % Rangos para los parámetros
out_names={'theta','dtheta'};%nombres de las salidas del modelo
T = gsua_dataprep(model,Ranges,ParIn,'out_names',out_names)%Creación de la tabla que controla el sistema
Setting environment to work with simulink Setting nominal values All done!
T = 7×2 table
 RangeNominal
1 g9.500010.500010
2 l243
3 m375
4 f010.5000
5 x10-110
6 x20-110
7 u-20200
Creamos una parametrización para el modelo que usaremos como valores nominales
T.clase=[9.8 3 5 0.5 0 0 10]'
T = 7×3 table
 RangeNominalclase
1 g9.500010.5000109.8000
2 l2433
3 m3755
4 f010.50000.5000
5 x10-1100
6 x20-1100
7 u-2020010
Evaluamos el comportamiento del modelo para dicho conjunto nominal de parámetros
gsua_eval(T.clase,T,0:0.1:150);%comprobamos la respuesta del sistema
Extraemos los valores de parámetros que estamos usando como nominales para crear una matriz de diseño experimental
nom=T.clase;
u_nominal=10;%entrada nominal
Creamos la matriz de diseño repitiendo los valores de la nominal y variando los valores de la entrada
Mdis=T.clase.*ones(1,101)+u_nominal*(-1:0.02:1).*[0 0 0 0 0 0 1]'
Mdis = 7×101
9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 9.8000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 2.0000 2.2000 2.4000 2.6000 2.8000 3.0000 3.2000 3.4000 3.6000 3.8000 4.0000 4.2000 4.4000 4.6000 4.8000 5.0000 5.2000 5.4000 5.6000 5.8000 6.0000 6.2000 6.4000 6.6000 6.8000 7.0000 7.2000 7.4000 7.6000 7.8000 8.0000 8.2000 8.4000 8.6000 8.8000 9.0000 9.2000 9.4000 9.6000 9.8000
Obtenemos la respuesta del sistema a cada entrada
y=gsua_eval(Mdis,T,0:0.1:100);
Progress: 100% Estimated processing time (h:m:s): 0:0:15 Remaining time (h:m:s): 0:0:0 Elapsed time (h:m:s): 0:0:15 Estimated stop time (h:m:s): 19:27:6 Number of simulations: 101
El tiempo que elegimos debe ser suficiente para que todas las respuestas del sistema se estabilicen (si alguna no se estabiliza es necesario cambiar el rango en el que se varía la entrada). Por tal razón, decidimos aumentar el tiempo de simulación haciendo
T.Properties.CustomProperties.Domain=[0 300];
y=gsua_eval(Mdis,T,0:0.1:300);
Progress: 100% Estimated processing time (h:m:s): 0:0:15 Remaining time (h:m:s): 0:0:0 Elapsed time (h:m:s): 0:0:15 Estimated stop time (h:m:s): 19:27:24 Number of simulations: 101
Ahora extraemos los valores donde se estabilizó el sistema
yss1=y(:,end,1);
yss2=y(:,end,2);
Y procedemos a trazar la curva de linealidad
uin=u_nominal*(-1:0.02:1)+u_nominal;
figure(1)
clf
plot(uin,yss1)
hold on
plot(uin,yss2)
legend({'Estado 1','Estado 2'})
xlabel('U')
ylabel('Y_{ss}')
Volvemos a escribir las ecuaciones del modelo para encontrar el punto de equilibrio
syms x_1 x_2 m f g l u(t) %symbolic state variables and parameters
%Defining the system of differential equations
eq1 = x_2;
eq2 = -f*x_2/m - g*sin(x_1)/l + u(t)/(m*l);
%Array with the system
eqs=[eq1; eq2]
eqs = 
Especificamos quiénes son las variables e igualamos las ecuaciones a 0 para encontrar los puntos de equilibrio
vars = [x_1 x_2];
sol=solve(eqs==zeros(2,1),vars)
sol = struct with fields:
x_1: [2×1 sym] x_2: [2×1 sym]
Exploramos las soluciones
sol.x_1
ans = 
sol.x_2
ans = 
Sustituimos en las expresiones anteriores los valores de los parámetros requeridos
sol1=double(subs(sol.x_1,{'g' 'm','u'},{9.8 5,10}));
%valores puntos de equilibrio
x1_0=sol1(2)
x1_0 = 0.2055
x2_0=0
x2_0 = 0
Miremos si estos valores coinciden con los valores de la simulación
ynom=gsua_eval(T.clase,T,0:300);
ynom(end,end,:)
ans =
ans(:,:,1) = 0.2055 ans(:,:,2) = 9.4666e-08
Para obtener las matrices del sistema lineal podemos proceder con linmod o calcular el jacobiano directamente
ja1=jacobian(eqs,vars)
ja1 = 
ja1=subs(ja1,{'g','l','m','f','x_1'},{9.8,3,5,0.5,x1_0})
ja1 = 
ja2=jacobian(eqs,u)
ja2(t) = 
ja2=subs(ja2,{'l','m'},{3,5})
ja2(t) = 
Valores numéricos del jacobiano
Esta sería la matriz A:
double(ja1)
ans = 2×2
0 1.0000 -3.1979 -0.1000
Esta sería la matriz B
double(ja2)
ans = 2×1
0 0.0667
Las matrices C y D son fácilmente deducibles de la estructura del modelo no lineal
Obtenemos todas las matrices necesarias para la linealización, nótese cómo las matrices coinciden con las obtenidas anteriormente
%linmod
[A,B,C,D]=linmod('PenduloC',[x1_0 x2_0],u_nominal)
A = 2×2
0 1.0000 -3.2632 -0.1000
B = 2×1
0 0.0667
C = 2×2
1 0 0 1
D = 2×1
0 0
Revisemos también la estabilidad de este punto crítico mediante los autovalores de la matriz A
eig(A)
ans = 2×1 complex
-0.0500 + 1.8057i -0.0500 - 1.8057i
Se observa que todos los autovalores tienen parte real negativa y por ende el punto de equilibrio es estable
Ahora, debemos tener un modelo para realizar la comparación tal y como se señala en la guía (PenduloLineal.slx). También es posible comparar con un modelo de espacio de estado (ss).

Control

Función de transferencia SISO

Para hacer el control SISO debemos elegir una de las salidas y una de las entradas del modelo lineal. Es recomdable elegir como salida uno de los estados que esté siendo afectado por la entrada. Si miramos la matrix B que se obtuvo en la linealización, nos damos cuenta de que la entrada actúa sobre el segundo estado y por ende lo elegiremos para controlarlo definiendo . Ahora, recordemos de las notas de clase que la función de transferencia de un sistema lineal en la forma de espacio de estado viene dado por
De forma tal, hacemos
syms s t
G=[0 1]*inv(s*eye(2)-A)*B
G = 
Esta sería la función de transferencia para nuestro péndulo si solo consideramos una entrada y tomamos como salida la velocidad angular (estado 2).

La forma básica de un sistema de control SISO

La función de transferencia de un sistema de control

Para obtener la función de transferencia del sistema de control debemos encontrar la relación , para ello nos valemos de los principios vistos en clase y notamos que
Ahora, recordemos que , donde y son polinomios en s, siendo el polonimio del numerador de grado menor o igual al polinomio del denominador. De forma tal

Aplicación del método de Routh-Hurwitz

El objetivo de un sistema de control es determinar una función de transferencia para el controlador que otorgue ciertas propiedades a la función de transferencia resultante para el sistema de control. En nuestro caso, la función de transferencia del sistema de control es una constante desconocida lo cual constituye el ejemplo más simple de control. Ahora bien, lo que podemos pedirle a este controlador es que mantenga el sistema estable. Resulta que no todos los valores de K necesariamente resultarán en un sistema estable. Para analizar qué valores de K son apropiados podemos aplicar el método de Routh-Hurwitz.

Explicación del método

sea un polinomio en S, el arreglo de Routh-Hurwitz viene dado por
Para el caso de n par, y por
Para el caso de n impar.
Ahora, según el teorema demostrado por Routh y Hurwitz, el número de cambios de signo en la primera columna de coeficientes equivale al número de raíces en el semiplano derecho. Adicionalmente, cada fila representa un polinomio auxiliar del grado que se señala en la primera columna y con los coeficientes dados en la fila. El grado de cada polinomio auxiliar se reduce de dos en dos para cada término siguiente. Existen dos casos especiales que pueden ocurrir en el arreglo de R-H que trataremos a continuación

Caso de un cero en la fila

Partamos del polinomio , el arreglo de R-H sería
Nótese que el 0 se toma como y se reemplaza por ϵ. Finalmente, para identificar el comportamiento de las raíces del polinomio, se hace
Por tanto, se concluye que este polinomio tiene 2 raíces en el semiplano derecho (si fuera un polinomio característico el sistema sería inestable).

Caso de una fila de ceros

Partamos del polinomio
Para este caso lo que hacemos es derivar el polinomio auxiliar que se encuentra directamente sobre la fila de ceros y reemplazamos dicha fila de ceros con los coeficientes del polinomio derivado. Tenemos entonces que el polinomio auxiliar sobre la fila de ceros es y así, . De forma tal
Se concluye entonces que el polinomio tiene 2 raíces en el semiplano derecho. El caso de la fila de ceros se presenta cuando el polinomio tiene raíces que son simétricas con respecto al origen. El arreglo de R-H nos da una condición necesaria para que todas las raíces del polinomio se encuentren en el semiplano izquierdo: que todos los coeficientes del polinomio original sean tengan el mismo signo.

Aplicación del método al caso del péndulo

Tenemos que la función de transferencia SISO de nuestro péndulo viene dado por
G
G = 
De forma tal, si implementamos un sistema simple de control tendremos que dicho sistema tendrá la siguiente función de transferencia
cuya estabilidad depende del polinomio resultante en el numerador . Para obtener este polinomio en MatLab hacemos
[Pn,Pd]=numden(G)%extraemos el numerador y el denominador
Pn = 
Pd = 
syms k
Pch=Pd+k*Pn
Pch = 
Polinomio característico =
Por tanto, para que no tengamos cambios de signo en la primera columna de coeficientes tendremos que cumplir que
Con esta información podemos implementar el sistema de control para el péndulo (ver PenduloControl en simulink). Nótese que la condición que encontramos para k corresponde al sistema lineal, sin embargo, el péndulo es un sistema no lineal, por tanto, el rango de linealidad podría ser mayor o menor. La siguiente imagen es un ejemplo de aplicación de un control basado en una linealización a un modelo no lineal.

Control de un sistema discreto

Para un sistema discreto tendríamos algo de la forma
El análisis de estabilidad de un controlador discreto es un poco diferente al control de un sistema continuo en lo que respecta a la determinación de la región de estabilidad del controlador. De hecho, es considerablemente más complejo controlar un sistema en tiempo discreto dado que el controlador solo nota los cambios en el sistema tras cada instante de muestreo. Si bien existe un método análogo a R-H para determinar la estabilidad del polinomio característico en tiempo discreto (llamado método de Jury), éste tiene condiciones considerablemente engorrosas de comprobar, por tal razón, en este curso optaremos por algunas transformaciones que nos lleven del plano Z al plano S para aplicar nuevamente el método de R-H.

Equivalencia de Z y S

Recordando que podemos hacer, mediante una aproximación por series de Taylor
Sin embargo, esta aproximación tiene algún grado de error con el que tendremos que lidiar al aplicar el control.

Transformación Bilineal

Otra alternativa consiste en llevar z a una nueva variable (w) cuyo plano cumple con algunas características importantes. Sea
Resulta que el interior del círculo unitario coincide con el semiplano izquierdo en w. De esta forma, aunque el plano w no es equivalente al plano s sí se cumple que la región estable de ambos es equivalente y por ende se puede apicar el método de R-H al polinomio en w para determinar el rango de estabilidad. La desventaja de este método es que no permite ubicar los polos de la función de transferencia en un sitio específico para que el sistema tenga ciertas propiedades más allá de la estabilidad.