Ing. Daniel Carbonel
EE616M
CONTROL ÓPTIMO PARA EL PENDULO INVERTIDO El sis sistem temaa de péndul péndulo o invert invertido ido mostr mostrado ado,, la gobier gobiernan nan las ecuaci ecuacione oness mostradas.
Y
m
x
θ
.. .. (M + m) x + m l θ = F .. .. x + l θ = g θ
mg
l F
M
X
Donde: M m
x g l
F
masa del carro (2000 gr) masa de la esfera (500gr) Posición horizontal de la varilla gravedad (10m/sg2) longitud de la varilla (1 mt) fuerza a controlar para mantener la varilla vertical.
Si definimos las variables de estado: x = θ 1 y 1 = θ y las variables de salida como:
. = x 2 θ y 2 = x
x 3 = x
. = x 4 x
El trabajo se divide en dos partes: 1°)
Diseño Diseño de de un siste sistema ma de de contr control ol ópti óptimo, mo, tal tal que que haga haga míni mínimo mo el el índic índice e de comportamiento cuadrático: cuadrático: J =
2°)
∞
∫ (x Q x + u R u ) dt T
T
0
Import Importanc ancia ia de de la elecci elección ón de un pond pondera eramie miento nto adecua adecuado do para para Q y R, R, most mostra rand ndo o en una una simu simula laci ción ón la vari variac ació ión n de las las resp respue uest stas as en función de ellas.
1
Ing. Daniel Carbonel
EE616M
Solución:
Primer Paso:
-
Hallemos las ecuaciones de estado que gobiernan este sistema. Reemplazando los valores dados de M, m y l Entonces tenemos: ..
..
2.5 x + 0.5θ = F ...(1) ..
..
...(2)
x + θ = 10θ
Reemplazando (2) en (1) obtenemos: ... (3)
..
2θ − 25θ + F = 0
... (4)
..
2 x + 5θ − F = 0
La ecuación de estados:
x1 0 x2 25 / 2 = x3 0 − 5/ 2 4 x .
.
.
.
1
0
0
0
0
0
0
0
x1 0 0 x 2 − 1 / 2 + 0 F 1 x3 0 x 4 1 / 2 0
y la ecuación de salida: x1
y1 1 = y 2 0 -
0
0
0 1
0 x3 x 4 0 x 2
Segundo paso:
Debemos diseñar un sistema de control óptimo, tal que haga mínimo el índice de comportamiento cuadrático. ∞
J
= ∫ ( x
T
Qx + uT Ru)dt
0
Para diseñar este sistema de control cuadrático, necesitamos hallar el valor de K, por lo cual usaremos Matlab, para obtener dicho vector. Para usar Matlab, necesitamos las matrices Q y R, que le asignaremos valores adecuados, estos valores son: 2
Ing. Daniel Carbonel
EE616M
Primer Caso Q=I, R=1. A=[0 1 0 0; 25/2 0 0 0; 0 0 0 1; -5/2 0 0 0]; B=[0; -1/2; 0; 1/2]; Q=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]; R=1; [K,P,E] = LQR(A,B,Q,R) donde se obtiene: K = [-65.8202 -19.1914 -1 -3.0272] P=
618.4727 183.6550 19.1914 52.0145
183.6550 54.8753 6.0820 16.4925
E=
-3.7968 -3.3045 -0.4904 + 0.3975i -0.4904 - 0.3975i
19.1914 6.0820 3.0272 4.0820
52.0145 16.4925 4.0820 10.4380
Este vector E son los eigenvalores de (A-BK), donde se ve que realmente el sistema es estable. Analicemos Analicemos la respuesta respuesta temporal, de las variables variables de estado, para ello consideremos la siguiente condición inicial x(0)=[0 0 1 0], y usando Matlab.
A=[0 1 0 0; 25/2 0 0 0;0 0 0 1; -5/2 0 0 0]; B=[0; -1/2; 0; 1/2]; K =[-65.8202 -19.1914 -1 -3.0272]; t=0:0.01:10; AA=A-B*K; sys=ss(AA,eye(4),eye(4),eye(4)); x=initial(sys,[0 0 1 0],t); x1=[1 0 0 0]*x'; x2=[0 1 0 0]*x'; x3=[0 0 1 0]*x'; x4=[0 0 0 1]*x'; subplot(4,1,1); plot(t,x1); grid; title('Respuesta a condición inicial') ylabel('x1') subplot(4,1,2);plot(t,x2);grid; ylabel('x2') subplot(4,1,3);plot(t,x3);grid; ylabel('x3') subplot(4,1,4);plot(t,x4);grid; ylabel('x4') xlabel('t segundos')
3
Ing. Daniel Carbonel
EE616M
Segu Segund ndo o Cas Caso o Q11=100, R=1. A=[0 1 0 0; 25/2 0 0 0;0 0 0 1; -5/2 0 0 0]; B=[0;-1/2;0;1/2]; Q=[100 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]; R=1; [K,P,E] = LQR(A,B,Q,R) donde obtenemos: K =[-67.3904 -19.4360 -1 -3.0787] P= 651.2569 188.3784 19.4360 53.5975
188.3784 56.1362 6.2390 17.2642
19.4360 6.2390 3.0787 4.2390
53.5975 17.2642 4.2390 11.1069
E= -3.6189 + 0.6475i -3.6189 - 0.6475i -0.4704 + 0.3856i -0.4704 - 0.3856i
4
Ing. Daniel Carbonel
EE616M
Este vector E son los eigenvalores de (A-BK), donde se ve que realmente el sistema es estable. Analicemos Analicemos la respuesta respuesta temporal, de las variables variables de estado, para ello consideremos la siguiente condición inicial x(0)=[0 0 1 0], y usando Matlab.
A=[0 1 0 0; 25/2 0 0 0;0 0 0 1; -5/2 0 0 0]; B=[0;-1/2;0;1/2]; K =[-67.3904 -19.4360 -1 -3.0787]; t=0:0.01:10; AA=A-B*K; sys=ss(AA,eye(4),eye(4),eye(4)); x=initial(sys,[0 0 1 0],t); x1=[1 0 0 0]*x'; x2=[0 1 0 0]*x'; x3=[0 0 1 0]*x'; x4=[0 0 0 1]*x'; subplot(2,2,1);plot(t,x1);grid; title('Respuesta a condición inicial para Q11=100 y R=1') ylabel('x1') subplot(2,2,2);plot(t,x2);grid; ylabel('x2') subplot(2,2,3);plot(t,x3);grid; ylabel('x3') subplot(2,2,4);plot(t,x4);grid; ylabel('x4') xlabel('t segundos')
5
Ing. Daniel Carbonel
EE616M
Terc Tercer er Cas Caso o Q11=100, R=0.01. A=[0 1 0 0; 25/2 0 0 0; 0 0 0 1; -5/2 0 0 0]; B=[0;-1/2;0;1/2]; Q=[100 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]; R=0.01; [K,P,E] = LQR(A,B,Q,R) de donde obtenemos: K =[ -193.4239 -51.4573 -10 -20.9009] P= 52.4299 12.7393 5.1457 8.8708
12.7393 4.4528 1.8842 3.4237
5.1457 1.8842 2.0901 1.6842
8.8708 3.4237 1.6842 3.0056
E= -6.8117 + 3.0614i -6.8117 - 3.0614i -0.8274 + 0.4603i -0.8274 - 0.4603i
Este vector E son los eigenvalores de (A-BK), donde se ve que realmente el sistema es estable. Analicemos Analicemos la respuesta respuesta temporal, de las variables variables de estado, para ello consideremos la siguiente condición inicial x(0)=[0,0,1,0], y usando Matlab.
A=[0 1 0 0; 25/2 0 0 0; 0 0 0 1; -5/2 0 0 0]; B=[0;-1/2;0;1/2]; K =[ -193.4239 -51.4573 -10 -20.9009]; t=0:0.01:10; AA=A-B*K; sys=ss(AA,eye(4),eye(4),eye(4)); x=initial(sys,[0 0 1 0],t); x1=[1 0 0 0]*x'; x2=[0 1 0 0]*x'; x3=[0 0 1 0]*x'; x4=[0 0 0 1]*x'; subplot(2,2,1);plot(t,x1);grid; title('Respuesta a condicion inicial para Q11=100 y R=0.01') ylabel('x1') subplot(2,2,2);plot(t,x2);grid; ylabel('x2') subplot(2,2,3);plot(t,x3);grid; ylabel('x3') subplot(2,2,4);plot(t,x4);grid; ylabel('x4') xlabel('t segundos')
6
Ing. Daniel Carbonel
EE616M
Cuart Cua rto o Cas Caso o Q11=100, Q33=100, R=0.01. A=[0 1 0 0; 25/2 0 0 0; 0 0 0 1; -5/2 0 0 0]; B=[0;-1/2;0;1/2]; Q=[100 0 0 0; 0 1 0 0; 0 0 100 0; 0 0 0 1]; R=0.01; [K,P,E] = LQR(A,B,Q,R) de donde obtenemos: K =[-421.0160 -128.0798 -100 -89.5562] P= 262.7280 81.5222 128.0798 73.1019
81.5222 27.0104 41.6016 24.4489
128.0798 41.6016 89.5562 39.6016
73.1019 24.4489 39.6016 22.6577
E= -7.1752 + 4.3526i -7.1752 - 4.3526i -2.4557 + 1.0339i -2.4557 - 1.0339i
7
Ing. Daniel Carbonel
EE616M
Este vector E son los eigenvalores de (A-BK), donde se ve que realmente el sistema es estable. Analicemos Analicemos la respuesta respuesta temporal, de las variables variables de estado, para ello consideremos la siguiente condición inicial x(0)=[0 0 1 0], y usando Matlab.
A=[0 1 0 0; 25/2 0 0 0; 0 0 0 1; -5/2 0 0 0]; B=[0;-1/2;0;1/2]; K=[-421.0160 -128.0798 -100 -89.5562]; t=0:0.01:10; AA=A-B*K; sys=ss(AA,eye(4),eye(4),eye(4)); x=initial(sys,[0 0 1 0],t); x1=[1 0 0 0]*x'; x2=[0 1 0 0]*x'; x3=[0 0 1 0]*x'; x4=[0 0 0 1]*x'; subplot(2,2,1);plot(t,x1);grid; title('Respuesta a condicion inicial para Q11=100, Q33=100 y R=0.01') ylabel('x1') subplot(2,2,2);plot(t,x2);grid; ylabel('x2') subplot(2,2,3);plot(t,x3);grid; ylabel('x3') subplot(2,2,4);plot(t,x4);grid; ylabel('x4') xlabel('t segundos')
8
Ing. Daniel Carbonel
EE616M
OBSERVACIONES: 1.
Se puede ver en el caso 2 (Q 11=100, R=1), que el tiempo de subida y el tiempo de establecimiento establecimiento son largos, y que al poner Q11 un valor grande como lo es 100, no hay mucha diferencia con el primer caso, debido a que solo se puso un valor alto en Q11 con respecto a Q22, Q33, Q44, pero no se le puso un valor valor consi consider derab able lemen mente te alto alto con con resp respect ecto o a R. Dand Dando o un valor valor alto alto a Q11, Q11, se está indicando que la variable de estado x1 (el valor angular) es la más importante para el sistema.
2.
En el tercer caso (Q 11=100, R=0.01), sí se puso un valor alto con respecto a R, pero solo se le puso ese valor alto a la variable x1, sin tener la misma consideración con la otra variable de estado x3 (posición del carro) respecto de la matriz de ponderación Q (Q33), que representa a la segunda variable de salida y2.
3. En este tercer caso también se puede notar que tanto el tempo de establecimiento, como el de subida se han reducido considerablemente con resp respec ecto to a los los caso casos s 1 y 2. Pero Pero un inco inconv nven enie ient nte e es el hech hecho o que que el sobrepico se haya incrementado en factores de 5 a 10 veces.
4.
En el cuarto caso (Q 11=100, Q33=100, =100, R=0.01 R=0.01), ), ocurr ocurre e que que los los tiempo tiempos s de establecimiento establecimiento y de subida, tienen valores muy pequeños debido al alto valor que tienen Q11 y Q33, con respecto a las otros elementos de la matriz de ponderación Q.
5.
Entre el tercer y el cuarto caso, se puede ver un crecimiento brusco de los sobrepicos, en un factor aproximado de 10 veces con respecto al tercer caso, por ejemplo la variable x1 llega a tomar un valor tan alto como –0.3 rad. (que es aproximadamente 17 0) y si aumentase mas ese valor de x1 con un mayor sobrepico, se estaría perdiendo esa linealidad asumida al principio que hizo que sen(θ )=θ (con respecto al modelo del péndulo), por lo que el modelo perdería validez.
9
Ing. Daniel Carbonel
EE616M
CONCLUSI ONCL USIONE ONES: S: 1.
Para obtener una respuesta más rápida de las variables x1 y x3, entonces Q11 y Q33 deben ser lo suficientemente suficientemente grandes comparados con Q22 y Q44 y R. Es decir, la ponderación ponderación debe poner en evidencia evidencia la importancia importancia de una u otra variable de estados, respecto de otras.
2.
No solo basta con aumentar Q11 (como se ve en los casos 1, 2 y 3) se tiene que que hace hacerr que que Q11 Q11 y Q33 Q33 sean sean más más gran grande des s con con resp respec ecto to a las las otra otras s componentes de la matriz Q, estas componentes reflejan la importancia que representan a las variables de estado x1 y x3, quienes a su vez son también las salidas del sistema
3.
Para elegir los valor adecuados de Q y R, se deben tener en cuenta las observaciones anteriores, ya que si queremos tener una respuesta rápida del sistem sistema, a, evid evident enteme ement nte e estar estaremo emos s hacien haciendo do crecer crecer el sobr sobrep epico ico y éste éste aumento del sobrepico tiene que ser un valor tolerable por el sistema, ello para el correcto funcionamiento del mismo.
4.
Notar también que todas las elecciones anteriores son óptimas, que minimizan J, pero la elección correcta de las matrices de ponderación, depende de las cara caract cter erís ísti tica cas s temp tempor oral ales es que que desea eseamo mos s del del sist sistem ema a (tie (tiemp mpo o de establecimiento establecimiento y del sobrepico sobrepico máximo). Evidentemente, Evidentemente, el decir deseamos deseamos es deci decirr lo que que nues nuestr tro o sist sistem ema a es capa capaz z de resp respon onde derr para para un buen buen funcionamiento.
5.
Final Finalmen mente, te, podría podría concl concluir uirse se que que un valor valor adecu adecuado ado de la matric matrices es de ponderación Q y R, tal vez serian Q11= Q33=50 y R = 0.01.
EL PROFESOR
10
