Facultad de Ingenieria UNC
Sistemas de Control II
Analisis y simulacion de
Sistemas de Control
No Lineales
con Matlab - Simulink
Autor: Ing Sergio O. Laboret
Electronica Industrial
Sistemas de Control
Mail: slaboret@[Link]
Objetivo: Presentar y utilizar las herramientas de Matlab 6 para resolver
analiticamente (Symbolic toolbox) y simular numericamente (Simulink) sistemas
dinamicos de control No Lineales
Requisitos: Se presupone un conocimiento basico de Matlab, Sistemas de
control no lineales, Plano de fase, Calculo y Ecuaciones Diferenciales
Temario
Introducion a Simulink y Modelado de Sistemas dinamicos no lineales
Introduccion a Symbolic toolbox y ejemplos aplicables al control
Controladores No lineales
Control todo-nada
Histeresis
Zona muerta
Control proporcional con saturacion
Alinealidades en la planta
Friccion No Lineal o de coulomb
Comparacion con friccion viscosa
Backlash (juego)
A la Salida
A al Entrada
Retardo de Tiempo puro
1
Acerca de Simulink
Simulink es una herramienta de Matlab para modelar y simular Sistemas Dinamicos es
decir basados en ecuaciones diferenciales ordinarias (ODEs) donde el tiempo es la
variable independiente que a su vez pueden ser lineales o no y considerados en tiempo
continuo o discreto (muestreado)
Proporciona una excelente interfaz grafica (GUI) para construir modelos en forma de
diagramas de bloques utilizando operaciones de tipo Drag and Drop, y no requiere
escribir codigo (excepto en algunos casos muy especificos)
Luego de definir un modelo se puede simular usando uno de los multiples sistemas
numericos de resolucion de ecuaciones diferenciales de que dispone
Creacion de los Modelos: Los modelos constan de: Bloques (Funciones) que se
interconectan mediante lineas de union (Seales) y comentarios aclaratorios
Librera: Simulink incluye una amplia biblioteca de bloques de fuentes, salidas,
componentes matematicos lineales y no lineales de tiempo continuo y discreto
Se abre con View->Show Library Browser , desde la barra de herramientas de
matlab o simulink o ejecutando simulink en la linea de ordenes de matlab y desde
aqu se navega en forma jerarquica sobre las categorias y los bloques
Corrida de la Simulacion: Se inicia comienza con el boton Start ( ) de la barra de
herramientas, se pausa con (ll) y se detiene con el boton stop ( ) o cuando finaliza,
tambien se puede utilizar el menu simulation->start, pause y stop
2
Modelado de un sistema de orden N que pueda ser escrito en FORMA NORMAL
dNx d N 1 x dx d N 1 x dx
Sea
N
= F ( N 1
,...... , x , t , u ), y = G ( N 1
,...... , x, t , u ) Con
d m x dt dt dt dt dt
= X m +1 para m=1,2,.N-1 (estados en forma canonica de variables de fase)
dt m
Ecuaciones de estado Integrando Obtenemos
t
dX 1
= X2 X 1 = X 1(0) + X 2(t )dt
dt .. 0t
dXN XN = XN (0) + F ( X 1, X 2,.. XN , U , t ) dt
= F ( X 1, X 2... XN , U , t )
dt 0
Ecuacion de Salida Y = G ( X 1, X 2... XN , U , t )
Las variables de estado se obtienen conectando integradores en cascada y la entrada
al primer integrador es funcion de dichas variables y la entrada al sistema
Libreria simbolica de matlab (Symbolic toolbox)
Permite efectuar calculos con variables y parametros simbolicos
Para definir cadenas de caracteres como objetos simbolicos se debe usar
variable=sym(expresion) o bien
syms lista de variables separadas por espacios
Calculo y ecuaciones diferenciales
El toolbox intenta encontrar una solucion explicita y si no puede una implicita
La variable independiente por defecto es t
Solucion de ecuaciones diferenciales ordinarias
solucion=dsolve(ecuacion,CCII,variable ind) Las derivadas se expresan en
ecuacion como D (Mayuscula) seguido del orden y el nombre de la variable a derivar
Integracion
integral indefinida: solucion=Int(funcion,variable indep)
integral definida entre a y b: solucion=Int(funcion,variable indep,a,b)
Derivacion
solucion=diff(funcion,orden,variable ind) Si se omite orden es 1
Nota: a veces estas funciones no toman todas la soluciones posibles
Por ejemplo dsolve(Dy=1/x,x) da log(x)+C y la solucion es log|x|+C
Ecuaciones algebraicas y funciones
solucion=solve(ecuacion,variable ind) Solucion de ecuaciones algebraicas
solucion=finverse(funcion,variable ind) Funcion inversa
solucion=simplify(expresion) simplifica expresiones (si es posible)
Aclaracion: para matlab log es logaritmo natural , el logaritmo base 10 es log10
3
Ejemplos Declaracion >>x=sym(x) declara a x variable simbolica
de variables >>syms a x y declara como simbolicas a,x,y
Ecuacion d 2Y dY
Diferencial + +Y = 0 >>y=dsolve('D2y+Dy+y=0','y(0)=1,Dy(0)=0)
Con CCII dt 2 dt
dY y=1/3*3^(1/2)*exp(-1/2*t)*sin(1/2*3^(1/2)*t)
Y (0) = 1, ( 0) = 0 +exp(-1/2*t)*cos(1/2*3^(1/2)*t)
dt
Integral Indefinida Y = AX 2 dX >>y=int(A*x^2,x)
1/3*A*x^3+C
Integral definida
con otra t
variable Y = f (t ) = AX 2 dX >>y=int('A*x^2','x'
,0,t)
o 1/3*A*t^3
d 2(X 3) >>y=diff('x^3',2,'x')
Diferenciacion y= 6*x
dX 2
Raices de >>solve('x^2+3*x+4=0')
ecuaciones X 2 + 3X + 4 = 0
[ -3/2+1/2*i*7^(1/2)] [ -3/2-1/2*i*7^(1/2)]
Funcion inversa >>syms a x; g=finverse(a*x^2,x)
Warning: finverse(a*x^2) is not unique
f ( X ) = AX 2 g = f 1 ( X ) g=1/a*(a*x)^(1/2)
Legibilidad
pretty(expresion) muestra la expresion en una forma mas legible en pantalla
Hacer >>g=1/a*(a*x)^(1/2);pretty(g) y comparar con el anterior resultado
Graficacion
ezplot(expresion,[desde,hasta]) grafica expresiones explicitas e implicitas en un
intervalo dado, si no se especifica es entre[ 2*pi,2*pi]
No admite parametros simbolicos, solo las variables a graficar: si hay un valor simbolico
en expresion es una grafica explicita si hay dos 2 implicita
Graficacion explicita con intervalo Graficacion implicita
>>y=log(x);ezplot(y,[0,100]) >>ezplot('x^2+y^2=10')
4
Ejercicio: Modelado matematico Modelo: masa_res_amort.mdl
Se desea modelar como diagrama de bloques
un sistema masa-resorte-amortiguador
sometido a la gravedad que parte de una K
posicion inicial en reposo Xo=15 X(t)
La ecuacion es
d2X dX
M + C + KX = Mg
dt 2 dt M
donde G=aceleracion gravedad=10,
M=masa, C= coef friccion, K=cte resorte C
2
dX d X dV Coeficiente
Definiendo: V= a= =
t
dt dt 2 dt P=Mg Rozamiento
t
V = Vo + a (t ) dt X = X 0 + V (t ) dt Son salidas de integradores en cascada
o o
Con esto construiremos el diagrama de bloques expresando la entrada al primer
integrador como funcion de las salidas de los mismos:
C K
Ma + CV + KX = Mg a=g V X
Posicion Final:
M M
t a = 0, V = 0, KX () = gM X () = gM = 1.11
K
Enfoque simbolico
Se resuelve la ODE
mediante dsolve y se
obtiene X
mediante derivacion con
diff se obtiene V y luego se
grafican ambas con con
ezplot escalando los ejes
con axis y dibujando la
grilla
M File: m_r_am.m
5
Comparativa de los dos enfoques: libreria simbolica y simulacion
Control No lineal en dos niveles
r(t) e(t) M U(t) Y(t)
Planta
+ -M
-
La Ecuacion de la Planta es
d 2Y dY A, M > 0
+A = U (t ) Y + AY = U (t )
dt dt
Asignando variables de estado La entrada a la planta es
X1 = Y X1 = X 2 + M , ( e > 0)
U (t ) =
X2 =Y X 2 = AX 2 + U (t ) M , ( e < 0)
6
Tenemos 2 sistemas lineales segun el signo del error
a )e < 0 U (t ) = M b)e > 0 U (t ) = + M
X1 = X 2 X1 = X 2
X 2 = ( AX 2 + M ) X 2 = M AX 2
Plano de fase
X2 =Y = X Y = Y Y+ = Y
haciendo
e<0 e>0
X 1 dY dt dY
Dividiendo las ecuaciones = =
y considerando
X 2 dt dX dX
dY X X
a) = Y = dX + C = F( X ) + C
dX M + AX M + AX
dY X X
b) + = Y+ = dX + C
dX M AX M AX
( X )
Y+ = d ( X ) + C = F ( X ) + C
M + A( X )
Esto significa que ambas Trayectorias de fase son simetricas respecto del origen
Entonces basta con analizar uno de los casos
lo haremos para el a) e<0 (*) utilizando la integracion simbolica de matlab escribiendo
>>y=int(-x/(m+a*x),x);pretty(y)
X M log(M + AX) M
Y = + +C La cual vale para X >
A A2 A
Debemos considerar la otra solucion:
X M log (M + AX ) M
Y = + +C La cual vale para X <
A A2 A
Asignando Valores: M=10, A=2
Resolviendo con matlab la primera ecuacion para que pase por el origen Y(0)=0
>>c1=solve('10/4*log(10)+c=0') da C1=-10/4*log(10),
Haciendo lo mismo para la segunda para que pase por Y(-20)=0
>>c2=solve('10+10/4*log(50)+c=0') da C2=-10-5/2*log(50)
(*) Nota: se hace para el caso de error negativo porque el otro caso conduce a una
expresion con valores imaginarios al aplicar y(0)=0, probar hacer
>>y=dsolve(Dy=x/10-2*x),y(0)=0,x);ezplot(y)
7
Trayectorias de fase con M=10, A=2
Para X>-5 con y(0)=0 Para X<-5 con y(-20)=0
>>y='-x/2+10/4*log(10+2*x)- >>y='-x/2+10/4*log(-10-2*x)-10-
10/4*log(10);ezplot(y,[-5,20]);grid 5/2*log(30);ezplot(y,[-20,-5]);grid
Nota: Aqui lo que obtenemos es la funcion inversa de la trayectoria de fase es decir
y en funcion de dy/dt , se hace de esta manera porque si hacemos dy/dt en funcion
de y da una relacion no expresable en terminos de funciones elementales, de hecho
no es una funcion univaluada, probar hacer:
>> y='-x/2+10/4*log(10+2*x)-10/4*log(10); x=finverse(y)
para obtener el plano de fase a) hay que girar las figuras 90 grados y reflejar sobre Y
y para el caso b) e>0 reflejar (a) sobre el origen
8
El plano de fase generico con M>A quedaria
X2=dy/dt
M/A
e<0
X1=Y
e>0
-M/A
Modelado de la planta con simulink
2
d 2Y dY dY
a=
d Y dV
=
+A = U (t ) Definiendo: V =
dt dt dt dt dt
t t
Son salidas de integradores
V = Vo + a(t )dt Y = Y 0 + V (t ) dt en cascada
o o
Con esto construiremos el diagrama de bloques expresando la entrada al primer
integrador como funcion de las salidas de los mismos:
a + AV = U (t ) a = U (t ) AV
9
Simulacion del modelo y Graficacion del plano de fase
Entrada=+10 Vo=-10, Yo=-2.5 Entrada=+10 Vo=10, Yo=-2.5
Entrada=-10 Vo=-10, Yo=2.5 Entrada=-10 Vo=-10, Yo=-2.5
10
Simulacion control de nivel discreto
Resultados de la simulacion
Vemos si bien el error tiende a 0 el
punto e=0 es un punto de equilibrio
inestable, ya que cualquier perturbacion
por minima que sea lo saca de ese
estado quedando un sistema que esta
permanentemente conmutando con una
frecuencia impredecible
Para ello se hace necesario hacer algunas modificaciones
11
Control discreto con histeresis
Con la Histeresis le agregamos memoria al controlador mediante un
elemento relay, el ancho de histeresis es 2h
U(t)
Elemento relay con M=10, h=0.5
-h h
M
e(t)
-M
e(t ) < h
U (t ) = M
h < e(t ) < h, U (t ) = M
e(t ) > h
U (t ) = + M
h < e(t ) < h, U (t ) = M
Las trayectorias de fase son como antes, lo que cambia son los valores de e(t) para los
cuales se producen las conmutaciones, es como si se solaparan los planos de fase
sobre el eje dy/dt
El plano de fase generico seria
X2=dy/dt
M/A
X1=Y
-M/A
12
Simulacion del control con histeresis
Resultados de la simulacion
Con h=0.5
Vemos que en este caso hay un
ciclo limite
El error es fluctuante y esta acotado
por la histeresis y la dinamica
La frecuencia de conmutacion del
controlador depende de dicha
histeresis y de la dinamica de la
planta pero ahora es previsible y/o Si bajamos la histeresis baja el error pero
calculable aumenta la frecuencia de conmutacion
13
Control con Zona Muerta
Ahora el controlador tiene tres niveles a)e(t ) < a U (t ) = M
U(t)
b)e(t ) > a U (t ) = + M
M Estos casos han sido ya analizados
-a Analizaremos el que falta:
a
e(t) c) e(t ) < a U (t ) = 0
-M
X1 = X 2
si X2 =Y = X X 2 = AX 2
X 1 dY dt dY
Dividiendo las ecuaciones = =
y considerando
X 2 dt dX dX
dY0 X Y
= = A Y0 = AdX + C = AX + C Y = +C
dX AX A
Son rectas de pendiente -1/A, dependiendo de las CCII el error puede no llegar a 0
Ahora el plano de fase seria con 3 zonas
X2=dy/dt
M/A
e<0
-a a
X1=Y
e>0
-M/A
14
Simulacion controlador con zona muerta
Nota: el orden en que van la ZM y el signo no puede ser alterado
Resultados de la simulacion
Con a=0.5
Vemos que en este caso no hay un
ciclo limite
El error es constante, finito y acotado
en a<e(inf)<a
Una vez que el sistema cae en la
zona muerta con condiciones
iniciales que no le permitan escapar
de ella el controlador ya no responde
U(t)=0 y el sistema queda librado a
si mismo
15
Control Proporcional con saturacion
Tambien tiene tres niveles K
si K es la ganancia proporcional
a)e(t ) < U (t ) = M
M
U(t) K
b)e(t ) > U (t ) = + M
M M
Estos casos han sido ya analizados
-M/K
Analizaremos el que falta por matlab y
e(t) luego por simulacion:
M/K
-M K
c) e(t ) < U (t ) = Ke(t )
M
Aqui consideraremos el sistema sin entrada e(t)=-y(t) y sujeto a distintas condiciones
iniciales resolviendo directamente el sistema y luego derivando con A=2, K=5
d 2Y dY d 2Y dY
+2 = 5Y +2 + 5Y = 0
dt dt dt dt
Plano de Fase
Utilizaremos el metodo de resolver directamente la ecuacion diferencial del
sistema y luego obtener la derivada
Lo haremos mediante la orden dsolve con condiciones iniciales y luego
mediante diff obtendremos la derivada y graficaremos y versus diff(y)
Programa matlab para trayectorias
16
Trayectorias superpuestas con CCII en los 4 cuadrantes
Simulacion
Planteamos el sistema equivalente con condiciones iniciales en los
cuatro cuadrantes
17
Resultado de la simulacion con K=5
Yo=-10, Vo=-10 Yo=10, Vo=10
Yo=-10, Vo=10 Yo=10, Vo=-10
Simulacion Control Proporcional con saturacion
18
Resultados con k=5 m=10 r(t)=10
La respuesta se lentifica durante el
periodo de saturacion ya que la
seal U(t) queda acotada en +M pero
este fenomeno es inevitable pues los
actuadores reales disponen de una
potencia limitada ya sea por razones
inherentes a su estructura fisica o
por el consumo de energia que
demandarian de la fuente de poder
Alinealidades en la planta
Huelgo (juego) o backlash
Es inherente a los sistemas de transmision mecanicos, correderas, engranajes,
tornillos, etc pude ser lineal o angular
Friccion No lineal
Freccion de coulomb: Es una fuerza contraria a la de actuacion que depende del
signo de la velocidad y no de su valor como la friccion viscosa, Fc=k*sign(V) con
Fc<=F
Friccion estatica: se opone al movimiento inicial entre dos superficies por ejemplo
la rueda de un auto contra el asfalto (que permite el movimiento) o un objeto en un
plano inclinado que no se cae hasta que opera en el una cierta fuerza
Generalmente se presentan combinados estos efectos con la friccion viscosa siendo
fenomenos complejos y muy difilies de modelar, normalmente se hacen ciertas
simplificaciones
Aparte dependen de varios factores: Lubricacion, temperatura, estado superficies,
fuerza normal
Retardo de tiempo
Tiempos muertos en la respuesta o en la medicion, por ejemplo un tren de
laminacion con control de espesor el cual se mide una cierta distancia D despues de
los rolos lo que da un tiempo de retardo de D*V donde V es la velocidad del material
Todas estas alinealidades son indeseables y se las trata de minimizar pues pueden
provocar mal funcionamiento, errores e incluso inestabilidad
19
Friccion de coulomb
H
r(t) e(t)
-
m(t) U(t) C(t)
K 1/Js 1/s
+ +
-
Quedan en dos semiplanos segun Si la entrada es un escalon
signo de velocidad
a )C > 0, Ke H = JC r (t ) = Ru (t ), e(t ) = R C (t )
b)C < 0, Ke + H = JC e(0) = R, e(0) = C (0) = 0
a )semiplano inferior : e < 0, Je + Ke H = 0
b)semiplano superior : e > 0, Je + Ke + H = 0
e
a) +ea = 0
=K/J
2
2
a = H / J 2 , a = H / K e
b) +e+a = 0
2
Haciendo en matlab
ea=dsolve('D2e/w^2+e-a=0','De(0)=0') a+C2*cos(w*t)
dea=diff(ea,'t') -C2*sin(w*t)*w
eb=dsolve('D2e/w^2+e+a=0','De(0)=0') -a+C2*cos(w*t)
deb=diff(eb,'t') -C2*sin(w*t)*w
se cumple que en el plano normalizado
2
e
+ (e a) = r 2
2
a) circulo radio r centrado en a donde r
depende de las CCII
2
e
+ (e + a) = r2
2
b) idem centrado en -a
20
Plano de fase del e /
error
(Normalizado)
y>0
-a
a e
y<0
Simulacion Con J=1, K=2 H=3 R=10u(t)
21
Resultados de la simulacion, Parametros: a=1.5 w=1.41
Las trayectorias son elipses en el
plano de fase no normalizado
Una vez que el modulo de la fuerza
ejercida es menor que la friccion y la
velocidad se anula el movimiento
cesa y hay un error de regimen, finito
y acotado en a<erp<a
Pico maximo Calculado del PF
Cmax=2(R-a)=17
tp=pi/w=2.22
Comparacion entre friccion de coulomb y viscosa
La envolvente de los picos o atenuacion en el caso coulomb es una recta y en el
caso viscoso es una exponencial, Esto puede servir para detectar el tipo
predominante de friccion presente en un sistema
22
Backlash: huelgo o juego Explicacion grafica
(poca Inercia)
Respuesta
temporal ante
exitacion senoidal
Caracteristica
entrada salida
Simulacion backlash a la Salida
Inercia del actuador alta comparada con la carga o juego en el elemento
de medicion
23
Resultados de la simulacion, Juego=0.1
la salida del elemento controlado
( medicion) se desconecta
momentaneamente de la salida
del actuador (salida) durante la
inversion de sentido del
movimiento
Resultados de la simulacion, Juego=0.5
24
Simulacion backlash a la Entrada
Inercia dela carga alta comparada con el actuador
Resultados de la simulacion, Juego=0.1
25
Resultados de la simulacion, Juego=0.5
Retardo de tiempo puro (Transport Delay) y Saturacion del elemento de
ganancia
26
Resultados de la simulacion, Retardo=0.1
Un retardo pequeo no
produce inestabilidad
Resultados de la simulacion, Retardo=0.5
Un retardo excesivo produce
inestabilidad y la saturacion
hace que se llegue a un ciclo
limite, si no estuviera habria
oscilaciones de amplitud
creciente
27
