Universidad Autónoma de Nuevo León
Facultad de Ingeniería Mecánica y Eléctrica
Control Optimo
Producto Integrador de Aprendizaje
Integrantes del equipo
Nombre Matricula Carrera
Gonzalo Alejandro Sánchez 1723034 IEA
Soto
Ramiro de Jesús Carbajal 1619867 IEA
Martínez
Ángel Daniel Pérez Delabra 1747673 IEA
Hora: N1-N3 Grupo: 001
Frecuencia: martes
Nombre del Catedrático: Juan Ángel Rodríguez Liñán
Planteamiento del problema, sistema o proceso a controlar
En el siguiente sistema de masa-resorte-amortiguador se va a calcular la fuerza
mínima para que la masa se mueva de una posición inicial cero a una posición de
referencia de 5 metros (considerando condiciones ideales), la entrada será la
fuerza aplicada y la salida la posición.
N
K=1.1
m
m=1 Kg
Ns
c=1.25
m
Obtención de un modelo en espacio de estado
Fuerza de la masa Fuerza del amortiguador Fuerza del resorte
2
d x dx
m c kx
dt dt
m ẍ+ c ẋ +kx=F Donde
x 1=x −c k 1
ẋ 2= x 2− x1 + F
m m m
x 2= ẋ
x 1=x ( posición)
ẋ 1=x 2
x 2= ẋ (velocidad)
ẋ 2= ẍ
u= F(fuerza)
[ ][ ] [ ] [ ][ ] [ ]
0 1 0 0 1 0
[ ]
x˙1
x˙2
=−k
m
−c
m
x1
x2
+ 1 u=¿
m
x˙1
x˙2 [ ]
= − (1.1) −1.25
1 1
x1
x2
+ 1 u
1
[ xx˙˙12]=[−1.1
0 1
][ ] [ ]
x1 + 0
−1.25 x 2 1
Como se puede apreciar en la ecuación anterior, se trata de un sistema de segundo
orden, por lo que el sistema tendrá dos variables de estado. Asignando una
variable de estado a la posición (x1) y otra variable a la velocidad (x2).
Análisis de propiedades del sistema y revisión de especificaciones del
objetivo de control
Las propiedades del sistema son las siguientes:
N
K=1.1 Constante del resorte
m
m=1 Kg Masa del sistema
Ns
c=1.25 Constante de amortiguación
m
El objetivo para este sistema es lograr un control de posición mediante la fuerza
aplicada, obteniendo una respuesta optima sin tener muchos sobrepasos tanto en
la entrada (fuerza) como en las condiciones iniciales (posición y velocidad) y
también lograr una rápida estabilización del valor de referencia.
Para iniciar el análisis del sistema, primero que nada, hay que verificar que este
sea un sistema controlable entonces realizamos la siguiente operación matricial.
det ⌈ CE ⌉ ≠ 0 CE=[B AB]
B=
[ 01] AB=
[−1.1
0 1 0
][ ] [
=
1
−1.25 1 −1.25 ]
det ⌈ CE ⌉=
[ 01−1.25
1
]= ( 0) (−1.25 )−( 1) ( 1)=−1
Como el determinante de CE es ≠ 0 entonces el sistema es completamente
controlable
Propuesta de una función de costo, en base a las especificaciones
Para poder estabilizar la posición de la masa en 5 metros mediante la entrada
(fuerza) se propuso la siguiente función de costo:
tf
J ( x , u , t ) =∫ ( 0.5 [ x 1−5 ] +0.5 u ) dt
2 2
t0
Cálculos de la técnica seleccionada para control óptimo
A=
[−1.1
0 1
−1.25 ] B=[ 01] H=
[ 00 00] Q=[ 10 00] R=[1] r (t )= [r52]
Ecuaciones de Riccati
T −1 T
0=−KA− A K−Q+ KB R B K
0=−
[ K 11 k 12
][
0 1
k 21 k 22 −1.1 −1.25
− ][
0 −1.1 K 11 k 12
1 −1.25 k 21 k 22
−
0 0
+ ][
1 0 K 11 k 12 0
k 21 k 22 1 ] [ ][
K k
[1][0 1] 11 12
k 21 k 22 ][ ] [ ]
0=¿−
[
−1.1 k 12 k 11 −1.25 k 12
−1.1 k 22 k 21−1.25 k 22
−
][
−1.1 k 21 −1.1 k 22
−
1 0 k 12
+
k 11 −1.25 k 21 k 12−1.25 k 22 0 0 k 22 ][ ][ ]
[1 k 21 1 k 22 ]
0=¿ [ 1.1 k 12 −k 11 +1.25 k 12
+
][
1.1 k 21 1.1 k 22 k k k k
][ ][
− 1 0 + 12 21 22 221
1.1 k 22 −k 21 +1.25 k 22 −k 11 +1.25 k 21 −k 12+1.25 k 22 0 0 k 22 k 21 k 22 ]
0=
[ 1.1 k 12 +1.1 k 21+ k 12 k 21−1 −k 11 +1.25 k 12+ 1.1k 22+ k 22 k 21
1.1 k 22−k 11 +1.25 k 21 +k 22 k 21 k 222−k 21+1.25 k 22−k 12 +1.25 k 22 ]
k 12=k 21
[ ]
2
k 12 + 2.2 k 12−1 −k 11 +1.25 k 12+1.1 k 22 +k 22 k 12
0= 2
1.1 k 22−k 11 +1.25 k 12 +2 k 22 k 12 k 22 −2k 12+2.5 k 22
k˙11=k 12 +2.2 k 12−1
2
k˙12= k 21=¿ 1.1 k 22−k 11˙+1.25 k 12 +k 22 k 12 ¿
k˙22=k 22 −2 k 12+ 2.5 k 22
2
Ecuaciones de referencia
T −1 T
ṡ=−A s+ KB R B s+Qr (t )
ṡ=− [ 01 −1.25 s2
+ ][ ] [
−1.1 s1 K 11 k 12 0
k 21 k 22 1 ][ ]s
[ 1 ][ 0 1 ] 1 + 1 0 5
s2 0 0 r2 [ ] [ ][ ]
ṡ=
[ 1.1 s2 k
][ ] []
+ 12 [ s 2 ]+
1.25 s 2−s1 k 22
5
0
ṡ=
[ 1.1 s 2+ k 12 s2 +5
1.25 s 2−s1 +k 22 s 2 ]
ṡ1=1.1 s2 +k 12 s 2+ ¿
ṡ2=1.25 s2−s 1+ k 22 s2
Ley de Control Optimo
−1 T −1 T
u ( t )=−R B Kx−R B s
u ( t )=−[ 1 ][ 0 1 ]
[ ][ ]
K 11 k 12 x 1
k 21 k 22 x 2 []
s
−[ 1 ][ 0 1 ] 1
s2
u ( t )=[−k 21 x1 −k 22 x 2 ] −s2
u ( t )=−k 21 x1 −k 22 x 2−s2
Primera aproximación basada en resultados numéricos y simulaciones y
comparación de las capacidades vs especificaciones del sistema
Programación
clear;
//Datos de simulación
t0=0;
tf=10;
paso=0.1;
ts=(t0:paso:tf)'; //tiempo para simulación
k0=[0;0;0;0;0];
x0=[0;0];
//Ec de Riccati
function dk=ecriccati(ts, k)//Ecs. de estado
k11=k(1); k12=k(2); k22=k(3); s1=k(4); s2=k(5);
dk(1) = (-1)*(k12^2+(11/5)*k12-1)
dk(2) = (-1)*((5/4)*k12-k11+(11/10)*k22+k12*k22)
dk(3) = (-1)*(k22^2+(5/2)*k22-2*k12)
dk(4) = (-1)*((11/10)*s2+k12*s2+5)
dk(5) = (-1)*((5/4)*s2-s1+k22*s2)
endfunction
//Simulacion de modelo
k=(ode('rk',k0,t0,ts,ecriccati))';//solucion de la ec de estado
t=tf-ts+t0; //tiempo real
k=k($:-1:1,:);// Reordena ganancias
t=t($:-1:1,:);// Reordena tiempo
//Modelo de estado
function dx=planta(t, x)//Ecs. de estado
x1=x(1); x2=x(2);
k11=k(1,1); k12=k(1,2); k22=k(1,3); s1=k(1,4); s2=k(1,5);
u=-k12*x1-k22*x2-s2;
dx(1) =x2
dx(2) =-1.1*x1-1.25*x2+u
mprintf('Las soluciones de la Ec. algebraica de Riccati son:\nk_11=%f\nk_12=%f\
nk_22=%f',k11,k12,k22)
mprintf('Las soluciones de la Ec. algebraica de Riccati son:\ns_1=%f\ns_2=
%f',s1,s2)
endfunction
//Simulacion de modelo
x=(ode('rk',x0,t0,t,planta))';//solución de la ec de estado
//Gráficas
figure(1)
plot(t,k(:,1),t,k(:,2),t,k(:,3),t,k(:,4),t,k(:,5));
xgrid
xtitle('Solución a la Ec. de Riccati','Tiempo t (s)','k(t)')
legend('$K_{11}$','$K_{12}$','$K_{22}$','$S_{1}$','$S_{2}$')
figure(2)
plot(t,x(:,1),t,x(:,2),t,5,'--');
xgrid
xtitle('Estados del sistema','Tiempo t (s)','x(t)')
legend('$x_1$','$x_2$','$r_1(t)$','$r_2(t)$')
figure(3)
plot(t,-k(1,2)*x(:,1)-k(1,3)*x(:,2)-k(1,5));
xgrid
xtitle('Entrada del sistema','Tiempo t (s)','u(t)')
Las soluciones de la Ec. algebraica de Riccati son:
K 11 =0.896987
k 12=0.386607
k 22=0.278304
s1=−5.142284
s2=−3.364946
Gráficas
u ( t )=−0.386607 x 1−0.278304 x 2 +3.364946
La fuerza requerida fue mínima pero la masa no logro llegar al valor de referencia
deseado por lo que hay que ir cambiando los valores de las matrices Q y R para así
lograr que la masa llegue al punto deseado con la fuerza mínima posible y con un
sobrepaso mínimo.
Ajustes a la propuesta y nuevos resultados numéricos y simulaciones
Se ajusto la función de costo ya que esta nos daba algunos problemas de sobrepaso
y sobre todo la función de costo planteada con anterioridad estaba mal planteada
Nueva función de costo:
tf
( 2 1
)
J ( x , u , t ) =∫ 12.5 [ x 1−5 ] + u2 dt
0 t
4
A=
[−1.1
0 1
−1.25 ] B=[ 01] H=
[ 00 00] Q=[ 250 00] R=[.5] r (t )= [r52]
Ecuación de Riccati
T −1 T
0=−KA− A K−Q+ KB R B K
0=−
[ K 11 k 12
k 21 k 22 ][
0
−1.1
1
−1.25
−
1 ][
−1.25 k 21 k 22
−
0 0
+ ][
0 −1.1 K 11 k 12 25 0 K 11 k 12 0
k 21 k 22 ][ ][
1
K k
[2][0 1] 11 12
k 21 k 22 ][ ] [ ]
0=¿−
[
−1.1 k 12 k 11 −1.25 k 12
−1.1 k 22 k 21−1.25 k 22
−
][ −1.1 k 21 −1.1 k 22
k 11 −1.25 k 21 k 12−1.25 k 22
− +
][ ][ ]
25 0 k 12
0 0 k 22
[2 k 21 2 k 22 ]
0=¿ [ 1.1 k 12 −k 11 +1.25 k 12
+
][ 1.1 k 21 1.1 k 22
1.1 k 22 −k 21 +1.25 k 22 −k 11 +1.25 k 21 −k 12+1.25 k 22
− 25 0 +
][ ][
2 k 12 k 21 2 k 22 k 21
0 0 2 k 22 k 21 2 k 222 ]
0=
[ 1.1 k 12+ 1.1 k 21+2 k 12 k 21−25 −k 11 +1.25 k 12+1.1 k 22 +2 k 22 k 21
1.1 k 22−k 11 +1.25 k 21 +2 k 22 k 21 2 k 222−k 21 +1.25 k 22−k 12+1.25 k 22 ]
k 12=k 21
[ ]
2
2 k 12 +2.2 k 12−25 −k 11 + 1.25 k 12+1.1 k 22 +2 k 22 k 12
0= 2
1.1 k 22−k 11 +1.25 k 12 +2 k 22 k 21 2 k 22 −2 k 12 +2.5 k 22
k˙11=2 k 12 + 2.2 k 12−25
2
k˙12= k 21=¿ 1.1 k 22−k 11 +1.25
˙ k 21 +2 k 22 k 21 ¿
k˙22=2k 22 −2 k 12+2.5 k 22
2
Ecuación de referencia
T −1 T
ṡ=−A s+ KB R B s+Qr (t )
ṡ=− [ 1 −1.25 s2
+ ][ ] [
0 −1.1 s1 K 11 k 12 0
k 21 k 22 1 ][ ]
[ 2 ][ 0 1 ] + 25 0 5
0 0 r2 [ ][ ]
ṡ=
[ 1.1 s2 k
][ ] [ ]
+ 12 [ 2 s2 ] +
1.25 s 2−s1 k 22
125
0
ṡ=
[ 1.1 s 2+2 k 12 s 2 +125
1.25 s 2−s 1 +2 k 22 s2 ]
ṡ1=1.1 s2 +2 k 12 s2 +125
ṡ2=1.25 s2−s 1+ 2k 22 s 2
Ley de control optimo
−1 T −1 T
u ( t )=−R B Kx−R B s
u ( t )=−[ 2 ][ 0 1 ]
[ K 11 k 12 x1
k 21 k 22 x2 ][ ] s
−[ 2 ][ 0 1 ] 1
s2 []
u ( t )=[−2 k 21 x 1−2 k 22 x2 ]−2 s2
u ( t )=−2 k 21 x 1−2 k 22 x2 −2 s 2
Simulaciones
clear;
//Datos de simulación
t0=0;
tf=10;
paso=0.1;
ts=(t0:paso:tf)'; //tiempo para simulación
k0=[0;0;0;0;0];
x0=[0;0];
//Ec de Riccati
function dk=ecriccati(ts, k)//Ecs. de estado
k11=k(1); k12=k(2); k22=k(3); s1=k(4); s2=k(5);
dk(1) = (-1)*(2*k12^2+(11/5)*k12-(25/1))
dk(2) = (-1)*((5/4)*k12-k11+(11/10)*k22+2*k12*k22)
dk(3) = (-1)*(2*k22^2+(5/2)*k22-2*k12)
dk(4) = (-1)*((11/10)*s2+2*k12*s2+(125/1))
dk(5) = (-1)*((5/4)*s2-s1+2*k22*s2)
endfunction
//Simulacion de modelo
k=(ode('rk',k0,t0,ts,ecriccati))';//solucion de la ec de estado
t=tf-ts+t0; //tiempo real
k=k($:-1:1,:);// Reordena ganancias
t=t($:-1:1,:);// Reordena tiempo
//Modelo de estado
function dx=planta(t, x)//Ecs. de estado
x1=x(1); x2=x(2);
k11=k(1,1); k12=k(1,2); k22=k(1,3); s1=k(1,4); s2=k(1,5);
u=-2*k12*x1-2*k22*x2-2*s2;
dx(1) =x2
dx(2) =-x1-x2+u
mprintf('Las soluciones de la Ec. algebraica de Riccati son:\nk_11=%f\nk_12=%f\
nk_22=%f',k11,k12,k22)
mprintf('Las soluciones de la Ec. algebraica de Riccati son:\ns_1=%f\ns_2=
%f',s1,s2)
endfunction
//Simulacion de modelo
x=(ode('rk',x0,t0,t,planta))';//solucion de la ec de estado
//Gráficas
figure(1)
plot(t,k(:,1),t,k(:,2),t,k(:,3),t,k(:,4),t,k(:,5));
xgrid
xtitle('Solución a la Ec. de Riccati','Tiempo t (s)','k(t)')
legend('$K_{11}$','$K_{12}$','$K_{22}$','$S_{1}$','$S_{2}$')
figure(2)
plot(t,x(:,1),t,x(:,2),t,5,'--');
xgrid
xtitle('Estados del sistema','Tiempo t (s)','x(t)')
legend('$x_1$','$x_2$','$r_1(t)$','$r_2(t)$')
figure(3)
plot(t,-2*k(1,2)*x(:,1)-2*k(1,3)*x(:,2)-2*k(1,5));
xgrid
xtitle('Entrada del sistema','Tiempo t (s)','u(t)')
Las soluciones de la Ec. algebraica de Riccati son:
k 11=12.543931
k 12=3.028058
k 22=1.223968
s1=−64.593980
s2=−17.467575
Graficas
u ( t )=−6.056116 x 1−2.447936 x 2−34.93515
Con los nuevos ajustes en las matrices Q y R se logró estabilizar la masa en la
posición de referencia con una mínima entrada (fuerza), y se logró tener una mejor
respuesta en la velocidad de la masa.
Implementación final en el sistema con la solución definitiva
Para poder realizar la implementación del sistema masa-resorte-amortiguador, se
utilizó el software de working model. A continuación, se mostrará el sistema:
Se obtuvieron los siguientes resultados de la simulación del sistema:
En las propiedades de la fuerza se introdujo el sistema de control de tal manera
que:
Body[1].p.x es el estado x1 es decir la posición
Body[1].v.x es el estado x2 es decir la velocidad
Conclusiones y recomendaciones
Conclusión y recomendaciones de Ramiro:
Se eligieron los valores de Q y R que se implementaron en la solución final por
que, en este sistema, si se le da preferencia a los estados con la matriz Q el
sistema logra llegar casi a la posición de referencia, pero el error es considerable
respecto a esta, por eso se decidió disminuir la R para tener una mayor fuerza de
entrada y así lograr llegar a la posición deseada, aunque incrementando un poco
más la velocidad.
Sen recomienda establecer primero una Q grande para observar el
comportamiento del sistema y si la entrada es satisfactoria dejarla en la unidad,
por el contrario, si aún no es suficiente se debería disminuir para compensar un
poco más la entrada.
Conclusión y recomendaciones de Angel:
Gracias a los resultados obtenidos por el control optimo, se puede observar que el
sistema requiere de una entrada (fuerza) variable para lograr un mejor tiempo de
respuesta sin tener tantas oscilaciones en las condiciones iniciales. A comparación
de una entrada constante, el sistema tardaría más tiempo en estabilizarse ya que
tendría una menor velocidad y por ende tardaría más en llegar al valor de
referencia deseado.
Conclusión y recomendaciones de Gonzalo:
Conforme estuvimos avanzando en el desarrollo del proyecto, pudimos observar
que para que el sistema se lograra mantener en una posición fija y sobre todo para
evitar los sobresaltos en este mismo se tenía que aplicar una fuerza variable.
Una recomendación que es que al principio se plante bien la función de costo ya
que esta si afecta a nuestro sistema y si se ponen valores al azar esta nos podrá
generar problemas desde el principio lo correcto es que se use la fórmula de
función de costo y saber si se le dará más preferencia a los estados o a la entrada,
en el tema de la simulación se presentó un problema el cual nos llevó días
solucionar este problema consistía en que no se podía vincular el matlab ni el
sacilab con workig model ya se necesitaba tener una cierta versión de matlab y
original para poder hacer esto así que para meter una fuerza variable lo preferente
es meter la función de entrada al working model.