Apuntes Control Optimo Cap 5-6

Apuntes de Control Óptimo Dr.
Juan Ángel Rodríguez Liñán
CAPÍTULO 5 - Control óptimo de seguimiento en sistemas

lineales utilizando criterios cuadráticos
5.1 Solución numérica de la ecuación diferencial de Riccati
Cuando no se cumplen las condiciones establecidas por Kalman, no es viable aplicar la ecuación
algebraica de Riccati, por lo que sería necesario resolver la ecuación diferencial de Riccati
K  KA  AT K  Q  KBR1BT K . Existen algoritmos que permiten calcular numéricamente la
solución de ecuaciones diferenciales a partir de conocer su derivada y sus condiciones iniciales. Como en
este caso se conocen las condiciones finales, entonces puede calcularse la solución numéricamente, pero
en tiempo inverso, partiendo de las condiciones finales como si fueran iniciales.
K (t )
K
K
K (t0 )  ?
t
K (t f )  H
Figura 4.0. La solución K(t) tiene derivada K (t ) en tiempo t, pero su derivada es  K (t ) en tiempo
inverso; K(t0) es desconocida, pero se conoce K(tf)=H. Ejemplo: K (0.5)  4 ,  K (0.5)  4 y K(tf)=0.2.
1 1
El proceso consiste en que la integración numérica de K(t) inicie en K(tf)=H, procediendo en “tiempo
inverso” ti = tf-t+t0 desde t=tf hasta t=t0, usando  K (t ) .
Puesto que  K (t )  K (ti ) , entonces la ecuación por integrar en tiempo inverso requiere cambio de
signo, es decir
K (ti )   K (ti ) A(ti )  AT (ti ) K (ti )  Q(ti )  K (ti ) B(ti ) R 1 (ti ) BT (ti ) K (ti )
36
Departamento de Electrónica y Automatización, FIME, UANL.
Apuntes de Control Óptimo Dr. Juan Ángel Rodríguez Liñán
Los valores de K(ti) son almacenados durante el proceso de solución y posteriormente reordenados
para tener K(t) en tiempo “directo”. La matriz de ganancia de retroalimentación es determinada por
K (t )  R1 (t ) BT (t )K (t ) para u *(t )  K (t ) x(t ) .
Ejemplo 1: Las ecuaciones diferenciales


k11  2k11  4k12  1  k12 2

k12   k11  4k12  2k 22  k12 k22

k22  2k12  6k 22  1  k22 2
se pueden resolver computacionalmente en tiempo inverso para K(tf)=[0 0;0 1], con t0=0, tf=10 y
un paso de integración de 0.1, de la siguiente forma:
clear;clc;
//Datos de simulacion
t0=0;
tf=10;
paso=0.1;
ti=(t0:paso:tf)'; //tiempo para simulacion
k0=[0;0;1];
//Ec de Riccati
function dk=ecriccati(ti, k)
k11=k(1); k12=k(2); k22=k(3);
dk(1) = -2*k11-4*k12+1-k12^2
dk(2) = k11-4*k12-2*k22-k12*k22
dk(3) = 2*k12-6*k22+1-k22^2
endfunction
//Simulacion de modelo
k=(ode('rk',k0,t0,ti,ecriccati))';//solucion de la ec de estado
t=tf-ti+t0; //tiempo real
//Gráficas
figure(1)
plot(t,k(:,1),t,k(:,2),t,k(:,3));
xgrid
xtitle('Solución a la Ec. de Riccati','Tiempo t (s)','k(t)')
legend('$K_{11}$','$K_{12}$','$K_{22}$')
La solución numérica de las ecuaciones de Riccati se ilustra en la Figura 2.
37
Figura 5.2. Solución numérica de las ecuaciones de Riccati.
5.2 Descripción del problema de seguimiento
Ahora se extienden los resultados obtenidos del problema del regulador lineal al problema de
seguimiento. La ecuación de estado es
x (t )  A(t ) x(t )  B(t )u (t )
y el criterio de desempeño a minimizar es
tf
1 1
J  eT (t f ) He(t f )   eT (t )Q(t )e(t )  u T (t ) R(t )u (t ) dt
2 2 t0
1
  H11e12 (t f )  H 22 e22 (t f )  ...  H nn en2 (t f )  (5.1)
2
tf
1
 
  Q11e12 (t )  Q22 e22 (t )  ...  Qnn en2 (t )    R11u12 (t )  R22u22 (t )  ...  Qmmum2 (t )  dt
2 t0
donde e(t )  x (t )  r (t ) es el error de seguimiento, r(t) es el valor deseado o referencia del vector de
estados, H(t) y Q(t) son matrices simétricas reales semidefinidas positivas, y R(t) es una matriz
simétrica real definida positiva. La minimización del criterio de desempeño se puede interpretar como
mantener el vector de estado cerca a la trayectoria de referencia sin un gasto excesivo del esfuerzo de
control.
38
5.3 Resultados de tiempo finito
Supóngase que el tiempo final es fijo. El Hamiltoniano es

1 T 1
H ( x(t ), u (t ), J *x , t )  e (t )Q (t )e(t )  u T (t ) R (t )u (t )  pT (t ) A(t ) x(t )  pT (t ) B(t )u (t )
2 2
Recordando que las condiciones necesarias para minimizar el Hamiltoniano (optimalidad) son13,14,15:
H H
0   R (t )u *(t )  BT (t ) p * (t )  0
u u
 2H
 R (t )  0
u 2
y a partir de la ecuación (1.10),
H
p * (t )    p *(t )  Q(t ) x *(t )  AT (t ) p *(t )  Q(t )r(t )
x
Resolviendo para el control u:
u *(t )  R1 (t ) BT (t ) p *(t )
y sustituyéndola en x *(t )  A(t ) x *(t )  B(t )u *(t ) se tiene
x *(t )  A(t ) x *(t )  B(t )R1 (t )BT (t ) p *(t ) .
Entonces, se pueden rescribir como un conjunto de ecuaciones diferenciales lineales variantes

en el tiempo no homogénea
 x * (t )   A(t )  B (t ) R 1 (t ) BT (t )   x * (t )   0 
 p * (t )       (0.25)
   Q (t )  AT (t )   p *(t )  Q (t )r (t ) 
La solución de estas ecuaciones es de la forma

 x * (t f )   x * (t )  t f  0 
 p * (t )    (t f , t )      (t f ,  )  d
 f   p *( t )  t0
 Q ( ) r ( ) 
 f1 (t ) 
donde la matriz de transición puede ser particionada y la integral remplazada por  .
 f 2 (t ) 
Entonces,
x * (t f )  11 (t f , t ) x *(t )  12 (t f , t ) p *(t )  f1 (t )
(0.26)
p * (t f )   21 (t f , t ) x *(t )   22 (t f , t ) p *(t )  f 2 (t )
De la condición de frontera de la ecuación HJB en t=tf para el funcional de costo (5.1),
encontramos que
p * ( t f )  Hx * (t f )  Hr (t f ) (0.27)
13 [ yT Mz] z  M T y
14 [ yT My] y  My  M T y
15 [ yT Mz] y  Mz
39
Sustituyéndola en (0.26), eliminando x*(tf), y resolviendo para p*(t), obtenemos

p *(t )  [22 (t f , t )  H 12 (t f , t )]1[ H 11 (t f , t )   21 (t f , t )]x *(t )
[22 (t f , t )  H 12 (t f , t )]1[ Hf1 (t )  Hr (t f )  f 2 (t )]
 K (t ) x *(t )  s(t )
Entonces, la ley de control es
u * (t )   R 1 (t ) BT (t ) K (t ) x (t )  R 1 (t ) BT (t ) s(t )
(0.28)
 K (t ) x(t )  s (t )
donde K (t ) es la matriz de ganancias de retroalimentación y s (t ) es la señal de control
auxiliar. Note que s (t ) depende de parámetros del sistema y de la señal de referencia. De
hecho, s (t ) depende de valores futuros de la señal de referencia. Entonces, el control óptimo es
anticipatorio en el sentido de que debemos determinar nuestra estrategia presente basados en la
posición actual y hacia donde se desea llevar al sistema.
Nuevamente, para evitar la necesidad de determinar la matriz de transición, se elegirá un
camino computacional. Diferenciando
p *(t )  K (t ) x *(t )  s(t ) (0.29)
sustituyendo de (0.25) para p *(t ) y x *(t ) , y utilizando (0.29) para eliminar p *(t ) , obtenemos
[ K (t )  K (t ) A(t )  AT (t ) K (t )  Q (t )  K (t ) B (t ) R 1 (t ) BT (t ) K (t )]x * (t )
[ s(t )  AT (t ) s (t )  Q (t ) r (t )  K (t ) B (t ) R 1 (t ) B T (t ) s (t )]  0
Puesto que esto debe satisfacerse para todo x(t) y r(t), se concluye que,
K (t )  K (t ) A(t )  AT (t ) K (t )  Q(t )  K (t )B(t ) R1 (t ) BT (t ) K (t ) (0.30)
s(t )  [ AT (t )  K (t )B(t ) R1 (t )BT (t )]s(t )  Q(t )r (t ) (1.31)
Puesto que (0.27) y (0.29) son

p * (t f )  Hx * (t f )  Hr (t f )
,
p * (t f )  K (t f ) x * (t f )  s (t f )
entonces las condiciones de frontera son
K (t f )  H
.
s (t f )   Hr (t f )
Para obtener las matrices de ganancias, se integran en tiempo “inverso” (0.30) y (1.31), tomando las
condiciones de frontera, y se almacenan los valores para K(t) y s(t).
Ejemplo 2:
x1  x2 (t )
0.5 x (t )  1  0.0025u (t)dt
tf
J (u )  
2 2
x2  2 x1 (t )  x2 (t )  u (t ) t0
1
Solución: Para este problema

0 1  0 0 1 0 0  1 1
A ; H   ;Q    ; B    ; R  0.0050  ; 𝑟(𝑡) =
 2 1 0 0 0 0 1  200 𝑟
40
Partiendo de la ec. de Riccati:

𝐾̇ (𝑡) = −𝐾(𝑡)𝐴(𝑡) − 𝐴 (𝑡)𝐾(𝑡) − 𝑄(𝑡) + 𝐾(𝑡)𝐵(𝑡)𝑅 (𝑡)𝐵 (𝑡)𝐾(𝑡)
𝑘̇ 𝑘̇ 𝑘 𝑘 0 1 0 2 𝑘 𝑘 1 0
=− − −
𝑘̇ 𝑘̇ 𝑘 𝑘 2 −1 1 −1 𝑘 𝑘 0 0
𝑘 𝑘 0 [200][ 𝑘 𝑘
+ 0 1]
𝑘 𝑘 1 𝑘 𝑘
𝑘̇ 𝑘̇ 2𝑘 𝑘 −𝑘 2𝑘 2𝑘 1 0 𝑘 𝑘 𝑘 𝑘
=− − − + 200
𝑘̇ 𝑘̇ 2𝑘 𝑘 −𝑘 𝑘 −𝑘 𝑘 −𝑘 0 0 𝑘 𝑘 𝑘 𝑘
Se obtienen las ecs. diferenciales:

𝑘 ̇ = −4𝑘 − 2 + 200𝑘
𝑘 ̇ = −𝑘 + 𝑘 − 2𝑘 + 200𝑘 𝑘
𝑘 ̇ = −2𝑘 + 2𝑘 + 200𝑘
Para la ecuación de la referencia:

𝑠̇ (𝑡) = −𝐴 (𝑡)𝑠(𝑡) + 𝐾(𝑡)𝐵(𝑡)𝑅 (𝑡)𝐵 (𝑡)𝑠(𝑡) + 𝑄(𝑡)𝑟(𝑡)
𝑠̇ 0 2 𝑘 𝑘 0 [200][ 𝑠 1 0 1
= − + 0 1] 𝑠 +
𝑠̇ 1 −1 𝑘 𝑘 1 0 0 𝑟2
𝑠̇ 0 −2 𝑘 𝑠 1
= + [200][0 1] 𝑠 +
𝑠̇ −1 1 𝑘 0
𝑠̇ 0 −2 0 𝑘 𝑠 1
= + 200 𝑠 + 0
𝑠̇ −1 1 0 𝑘
𝑠̇ 0 −2 𝑠 0 200𝑘 𝑠 1
= 𝑠 + 𝑠 +
𝑠̇ −1 1 0 200𝑘 0
𝑠̇ −2𝑠 200𝑘 𝑠 1
= + +
𝑠̇ −𝑠 + 𝑠 200𝑘 𝑠 0
Se obtienen las ecs diferenciales:
𝑠 ̇ = −2𝑠 + 200𝑘 𝑠 + 1
𝑠 ̇ = −𝑠 + 𝑠 + 200𝑘 𝑠
La ley de control óptimo es:

𝑢∗ (𝑡) = −𝑅 (𝑡)𝐵 (𝑡)𝐾(𝑡)𝑥(𝑡) − 𝑅 (𝑡)𝐵 (𝑡)𝑠(𝑡)
𝑘 𝑘 𝑥 𝑠
𝑢∗ (𝑡) = −[200][0 1] 𝑥 − [200][ 0 1 ] 𝑠
𝑘 𝑘
𝑘 𝑥 + 𝑘 𝑥
𝑢∗ (𝑡) = −[200][0 1] − [200][𝑠 ]
𝑘 𝑥 +𝑘 𝑥
𝑢∗ (𝑡) = −200𝑘 𝑥 − 200𝑘 𝑥 − 200𝑠
Puesto que el sistema es completamente controlable, las matrices son constantes, H=0 y la
referencia es constante, se calculan las ganancias óptimas para el caso de tf en horizonte infinito:
𝑘 = 0.4226
𝑘 = 0.0814
𝑘 = 0.0239
𝑠 = −1.2245
𝑠 = −0.2118
sustituyendo los valores calculados de las ganancias, la ley de control resulta en
𝑢(𝑡) = −200(0.0814)𝑥 − 200(0.0239)𝑥 − 200(−0.2118)
41
𝑢(𝑡) = −16.28𝑥 − 4.78𝑥 + 42.36
Esto se simula computacionalmente de la siguiente forma:
clear;clc;
//Datos de simulación
t0=0;
tf=10;
paso=0.1;
ti=(t0:paso:tf)'; //tiempo para simulación
k0=[0;0;0;0;0];
x0=[10;-10];
//Ec de Riccati
function dk=ecriccati(ti, k)
k11=k(1); k12=k(2); k22=k(3); s1=k(4); s2=k(5);
dk(1) = 4*k12+1-200*k12^2
dk(2) = k11-k12+2*k22-200*k12*k22
dk(3) = 2*k12-2*k22-200*k22^2
dk(4) = 2*s2-200*k12*s2-1
dk(5) = s1-s2-200*k22*s2
endfunction
k=(ode('rk',k0,t0,ti,ecriccati))';//solucion de la ec de estado
t=tf-ti+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);
k12=k(1,2); k22=k(1,3); s2=k(1,5); //Ganancias calculadas:
u=-200*k12*x1-200*k22*x2-200*s2;
dx(1) = x2
dx(2) = 2*x1-x2+u
endfunction
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));
42
xgrid
legend('$K_{11}$','$K_{12}$','$K_{22}$','$S_{1}$','$S_{2}$')
figure(2)
N=size(t,1);
plot(t,x(:,1),t,x(:,2),t,ones(N,1),'--')
xgrid
xtitle('Estados del sistema','Tiempo t (s)','x(t)')
legend('$x_1$','$x_2$','$r_1(t)$')
figure(3)
plot(t,-200*k(1,2)*x(:,1)-200*k(1,3)*x(:,2)-200*k(1,5));
xgrid
xtitle('Entrada del sistema','Tiempo t (s)','u(t)')
La solución numérica de las ecuaciones de Riccati se ilustra en la Figura 5.1.
La respuesta de los estados x1 y x2 se ilustra en la Figura 5.2 y la señal de control en la Figura 5.3 .
43
Figura 5.2. La respuesta de los estados x1 y x2.
44
Figura 5.3. Señal de control óptima u(t).

Ejemplo 3:
𝑥̇ = 𝑥 (𝑡)
𝑥̇ = 2𝑥 (𝑡) − 𝑥 (𝑡) + 𝑢(𝑡)
𝐽 = 0.5[𝑥 − 0.2𝑡 ] + ([𝑥 − 0.2𝑡] + 0.25𝑢 ) 𝑑𝑡
Solución: Para este problema

0 1 1 0
𝐴= ;𝐻=
2 −1 0 0
0 2 0 0.2𝑡
𝐵= ; 𝑅 = 0.50 = ;𝑄 = ; 𝑟(𝑡) =
1 0 0 𝑟
Partiendo de la ec. de Riccati, se obtienen las ecs. diferenciales:

𝑘 ̇ = −4𝑘 − 1 + 2𝑘
𝑘 ̇ = −𝑘 + 𝑘 − 2𝑘 + 2𝑘 𝑘
𝑘 ̇ = −2𝑘 + 2𝑘 + 2𝑘
𝑠 ̇ = −2𝑠 + 2𝑘 𝑠 + 0.4𝑡
𝑠 ̇ = −𝑠 + 𝑠 + 2𝑘 𝑠
La ley de control óptimo es:

𝑢∗ (𝑡) = −𝑅 (𝑡)𝐵 (𝑡)𝐾(𝑡)𝑥(𝑡) − 𝑅 (𝑡)𝐵 (𝑡)𝑠(𝑡)
𝑢∗ (𝑡) = −2𝑘 (𝑡)𝑥 − 2𝑘 (𝑡)𝑥 − 2𝑠 (𝑡)
45
Esto se simula computacionalmente de la siguiente forma:
clear; clc;
t0=0;
tf=10;
paso=0.1;
ts=(t0:paso:tf)'; //tiempo para simulación
k0=[1;0;0;-0.2*tf;0];
x0=[10;-10];
//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) = 4*k12+1-2*k12^2
dk(2) = k11-k12+2*k22-2*k12*k22
dk(3) = 2*k12-2*k22-2*k22^2
dk(4) = 2*s2-2*k12*s2-0.4*(tf-ts+t0)
dk(5) = s1-s2-2*k22*s2
endfunction
k=(ode('rk',k0,t0,ts,ecriccati))';//solucion de la ec de estado
t=tf-ts+t0; //tiempo real

//Modelo de estado
function dx=planta(t, x)//Ecs. de estado
x1=x(1); x2=x(2);
k12=interp1(ts,k(:,2),t,'linear'); //Interpolate the data (ts,k) at time t

s2=interp1(ts,k(:,5),t,'linear'); //Interpolate the data (ts,k) at time t
u=-2*k12*x1-2*k22*x2-2*s2;
dx(1) = x2
dx(2) = 2*x1-x2+u
endfunction
46
//Gráficas
figure(1)
plot(t,k(:,1),t,k(:,2),t,k(:,3),t,k(:,4),t,k(:,5));
xgrid
legend('$K_{11}$','$K_{12}$','$K_{22}$','$S_{1}$','$S_{2}$')
figure(2)
plot(t,x(:,1),t,x(:,2),t,0.2*t,'--')
xgrid
xtitle('Estados del sistema','Tiempo t (s)','x(t)')
legend('$x_1$','$x_2$','$r_1(t)$')
figure(3)
plot(t,-2*k(:,2).*x(:,1)-2*k(:,3).*x(:,2)-2*k(:,5));
xgrid
La solución numérica de las ecuaciones de Riccati se ilustra en la Figura 5.4.
La respuesta de los estados x1 y x2 se ilustra en la Figura 5.5 y la señal de control en la Figura 5.6 .
47

48
5.4 Control óptimo para seguimiento de casos particulares
Las ecuaciones (1.30) y (1.31) deben resolverse fuera de línea, lo que implica que la referencia r(t) debe
estar disponible para resolver la ecuación (1.31) en tiempo inverso, y posteriormente usar la solución en el
controlador (1.28). Sin embargo, en algunos casos es posible simplificar el control de seguimiento como
un problema de regulación en el origen, tal que el controlador (1.28) no requiera disponer de la referencia
r(t) previamente, sino sólo retroalimentarla en tiempo real. Lo cual se detalla a continuación:
Considere la ecuación de estado dada por

x (t )  A(t ) x (t )  B (t )u (t )
con el mismo criterio de desempeño dado por la ecuación (5.1):
tf
1 1
J  e(t f )T He(t f )   e(t )T Q(t )e(t )  uT (t ) R(t )u(t ) dt
2 2 t0
donde e(t )  x (t )  r (t ) es el vector de error del estado x(t ) respecto a la referencia deseada r (t ) .
Entonces con x (t )  e(t )  r (t ) , la dinámica del error puede escribirse en la forma
e(t )  x (t )  r(t )
 A(t ) x (t )  B (t )u (t )  r(t )
 A(t )[e(t )  r (t )]  B(t )u (t )  r(t )
e(t )  A(t )e(t )  B(t )u(t )  A(t )r (t )  r(t ) (1.32)
Puesto que el seguimiento de x(t ) a r (t ) es equivalente a que e(t )  0 , un control óptimo de

seguimiento puede resolverse como un problema de regulación en el origen de (1.32) en los siguientes
casos:
CASO 1)
El control óptimo es de la forma

u *(t )  R1 (t ) BT (t ) K (t )[ x(t )  r (t )]  B1 (t )[r(t )  A(t )r (t )] (1.33)
si B(t ) es una matriz invertible, donde K (t ) es solución de la ecuación de Riccati (1.30).
Ejemplo 4:
Considere un modelo de la forma

a a12   b1 0 
x   11 x  u , para todo b1  0 y b2  0
 a21 a22   0 b2 
49
CASO 2)
El control óptimo es de la forma

u *(t )  R1 (t ) BT (t ) K (t )[ x(t )  r(t )] (1.34)
si se cumple que A(t )r (t )  0 con r (t ) constante, o si se satisface r(t )  A(t )r (t ) ; donde K (t ) es
solución de la ecuación de Riccati (1.30).
Ejemplo 5:
Considere un modelo lineal del motor de corriente directa en vacío descrito en la forma
0 1   0 
x    x   Ks  u
 0  1   
 Ts   Ts 
donde Ts y Ks son parámetros de fabricante dependientes de la resistencia de armadura,
coeficiente de fuerza contraelectromotriz, momento de inercia y fricción viscosa. Suponga que el
criterio de costo es de la forma: t
J (u )  
t0
f
25 x (t )  r (t )  0.5 x (t )  500u (t )dt
1 1
2
2
2 2
donde r1 (t ) es una función constante a tramos dada por una señal cuadrada con amplitud entre
0 y 40°, con frecuencia de 0.5 Hz. Calcule el control óptimo.
Solución:
0 1   0 
A=  , B=   , 𝑅 = 1000, 𝑄 = 50 0 , 𝐻 = 0 0 , r   r1  .
 0  1   Ks  0 1 0 0  
0
 Ts   Ts 
Ya que este ejemplo corresponde al caso 2, entonces el control óptimo es de la forma

u *(t )  R1 (t ) BT (t ) K (t )[ x(t )  r(t )] ,
ya que A(t )r (t )  0 y r (t )  0 , y donde K (t ) debe calcularse como la solución de la ecuación
de Riccati.
Para calcular numéricamente las ganancias óptimas de este controlador, se puede utilizar
nuevamente la función lqr(·) de Scilab, para obtener automáticamente la solución a la ecuación
algebraica de Riccati. Para los parámetros del motor tome los valores KS= 186 [°/V-s] y TS= 1.04 s.
clear;clc
//Valores predeterminados
Ks=186;
Ts=1.04;
//Matrices del sistema

A=[0 1;0 -1/Ts];
B=[0;Ks/Ts];
50
//Matrices del funcional de costo

Q=diag([50,1]);
R=1000;
//Cálculo de ganancias del LQR

C=eye(2,2); //Matriz de salidas medibles
P = syslin('c',A,B,C); //Sistema continuo en espacio de estado
[Kb,K]=lqr(P,Q,R); //Cálculo de ganancias óptimas del LQR
k11=K(1,1); k12=K(1,2); k22=K(2,2); //Valor de cada ganancia
t0=0;
tf=10;
paso=0.1;
t=(t0:paso:tf)';
x0=[10;0]; // condiciones iniciales
//Señal de referencia
N=size(t,1);
for i=1:N;
if t(i)<2 then
r(:,i)=[0;0];
elseif t(i)<4 then
r(:,i)=[40;0];
elseif t(i)<6 then
r(:,i)=[0;0];
elseif t(i)<8 then
r(:,i)=[40;0];
else r(:,i)=[0;0];
end
end
r=r'; ts=t;
//Modelo de estado
function dx=planta(t, x)
r1=interp1(ts,r(:,1),t);
r2=interp1(ts,r(:,2),t);
u=Kb*[x-[r1;r2]];
dx=A*x+B*u
endfunction
x=(ode('rk',x0,t0,t,planta))';//Solución de la ec de estado
51
//Gráficas
figure(1)
plot(t,x(:,1),t,x(:,2),t,r(:,1),'--');
xgrid
xtitle('Estados del sistema continuo','Tiempo t (s)','x(t)')
legend('x1','x2','r1(t)')
figure(2)
plot(t,(Kb*(x-r)')');
xgrid
La respuesta de los estados x1 y x2 se ilustra en la Figura 5.7 y la señal de control en la Figura 5.8.
52
Ejemplo 6:
Considere un modelo de la forma

 0 1 0  sin(t ) 
x    x    u y referencia r   
 1 0  1  cos(t ) 
Proponga un criterio de costo J y calcule el control óptimo correspondiente.
Solución:
 0 1  0
A=   B=   ,
 1 0  1
Se propone un criterio de costo
5 x (t )  r (t )  0.5 x (t)  r (t)  0.5u (t )dt

tf
J (u )  
2 2 2
1 1 2 2
t0
Ya que este ejemplo corresponde al caso 2, entonces el control óptimo es de la forma

u *(t )  R1 (t ) BT (t ) K (t )[ x(t )  r(t )] ,
ya que A(t )r (t )  r(t ) , y donde K (t ) debe calcularse como la solución de la ecuación de
Riccati, para la simulación numérica de este ejemplo, se calcula la solución algebraica, como se
ilustra a continuación:
clear;clc
//Matrices del sistema

A=[0 1;-1 0];
53
B=[0;1];
//Matrices del funcional de costo

Q=diag([10,1]);
R=1;
//Cálculo de ganancias del LQR

k11=K(1,1); k12=K(1,2); k22=K(2,2); //Valor de cada ganancia
t0=0;
tf=20;
paso=0.1;
t=(t0:paso:tf)';
x0=[10;-10]; // condiciones iniciales
//Modelo de estado
function dx=plantac(t, x)//Ecs. de estado
x1=x(1);
x2=x(2);
r=[sin(t);cos(t)];
u=Kb*[x-r];
dx=A*x+B*u
endfunction
x=(ode('rk',x0,t0,t,plantac))';//Solucion de la ec de estado
//Gráficas
figure(1);clf,
plot(t,x(:,1),t,x(:,2),t,sin(t),'--',t,cos(t),':');
xgrid
xtitle('Estados del sistema continuo','Tiempo t (s)','x(t)')
legend('x1','x2','r1(t)','r2(t)')
figure(2);clf,
plot(t,(Kb*(x-[sin(t),cos(t)])')');
xgrid
La respuesta de los estados x1 y x2 se ilustra en la Figura 5.9 y la señal de control en la Figura

5.10.
54
55
El caso de control “u” no restringido no puede ser satisfecha para un verdadero sistema físico.
Debe recordarse que, si el control u se satura, entonces no satisface el problema de LQ . Para
volver al problema de LQ , la amplitud de la señal u debe disminuirse. Para ello, debe
considerarse algún criterio de diseño (matrices Q y R) acerca de los valores máximos aceptables
de las señales, como en el criterio establecido en el capítulo anterior. Para realizar eso las
herramientas de simulación son recomendadas.
56
CAPÍTULO 6 - Ejemplos de aplicación
En este capítulo se muestra la aplicación del control óptimo en algunos sistemas físicos,
ilustrando el proceso de cálculo e implementación.
Caso de circuito eléctrico
Considere el siguiente circuito RC en cascada de la Figura 7.1.
Figura 6.1. Circuito RC en cascada.

Cuyas ecuaciones son
ei  R1i1  vC1
,
vC1  R2i2  eO ,
dvC 1
C1  i1  i2
dt ,
dvC 2
C2  i2
dt .
x1  v c1 x 2  vc2 u  ei
Obtenga un modelo de estado si las variables de estado son , , es la
eO  vC 2
entrada y es la salida.
Calcule e implemente el control óptimo u* tal que minimice al criterio de costo

1000[ x  3  sin(0.4 t )]2  u 2 dt
tf
J  2
t0
Solución
Un modelo en espacio de estado es obtenido como sigue:
57
ei  R1i1  vC1
u  R1i1  x1
1
i1  (u  x1 )
R1
vC1  R2i2  vC 2
x1  R2i2  x2
1
i2  ( x1  x2 )
R2
dvC1 1 11 1 
 (i1  i2 )   (u  x1 )  ( x1  x2 ) 
dt C1 C1  R1 R2 
1  1 1  1
x1  u    x1  x2
R1C1  R1C1 R2C1  R2C1
dvC 2 1 1 1
 i2  ( x1  x2 )
dt C2 C2 R2
1 1
x2  x1  x2
R2C2 R2C2
  1 1  1 
      1 
 x1    R1C1 R2C1  R2C1   x1  
 x       R1C1  u
 2 1 1   x2   
    0 
 R2C2 R2C2 
Para ilustrar el control, se implementa un circuito con R₁= R₂= 1 kΩ y C₁= C₂= 470 μF. Entonces
las matrices para el regulador lineal cuadrático con seguimiento son:
 4.2553191 2.1276596  2.1276596  0 0  0 0   r1 

A  ,B    ,Q    , R   2 , H    ,r   
 2.1276596 2.1276596   0  0 2000  0 0  3  sin(0.4 t ) 
Aplicando la metodología de control óptimo para seguimiento, el controlador a implementar es

u(t )   R1BT  Kx  s   0.5 2.1276596 0  Kx  s 
u(t )  1.0638298  k11 x1  k12 x2  s1   k1 x1  k2 x2  k0
con k11, k12 y s1 como soluciones de la ecuación de Riccati.
El cual se puede simular computacionalmente para k (t f )   0 0; 0 0 , x(t0 )   0;0 y t f  30 s .

La solución de la ecuación de Riccati se ilustra en la Figura 6.2. La respuesta de los estados en el
tiempo en la Figura 6.3 y la señal de control en la Figura 6.4.
58
Figura 6.2. Solución de las ecuaciones de Riccati.
59
Para la implementación experimental de este controlador se propone utilizar una tarjeta

Arduino UNO. En la figura 6.5, se muestra un diagrama esquemático de conexiones.
Figura 6.5. Diagrama esquemático de conexiones del circuito con la tarjeta Arduino.
60
Por simplicidad, se considera implementar electrónicamente el controlador

u(t )  k1 x1  k2 x2  k0 con las ganancias de la ecuación algebraica de Riccati (para t f   ), que
fueron calculadas como k1  5.3831479 , k2  25.255435 y aproximadamente
k0  94.820931  31.55086532sin(0.4 t  0.1570162946) en estado estacionario. Con ello se
programa el siguiente algoritmo de control (que luego se compila y descarga en Arduino):
#include <math.h>
unsigned long t0;

float t;
int entrada_analogica1= A0; // declara numeros de pines de tipo entero
int entrada_analogica2= A1;
int salida = 6; // declara numeros de pines de salida PWM tipo entero
int valor_entrada_analogica1=0; // Declara conversion de valores entradas y salida a entero
int valor_entrada_analogica2=0;
float V1;
float V2;
float u=0;
int U=0;
float k1=-5.3831479;
float k2=-25.255435;
float k0=0;
float pi=3.1416;
void setup() {
pinMode(salida, OUTPUT); // Pines de salida y entrada

pinMode(entrada_analogica1, INPUT);
Serial.begin(9600); // Inicializar la comunicación serial a 9600 bits por segundo
}
void loop() { //Entra a un ciclo repetitivo
t0= millis(); //tiempo actual tomado del timer en milisegundos

t=t0*0.001;
k0=94.820931+31.55086532*sin(0.4*pi*t+0.1570162946);
valor_entrada_analogica1= analogRead(entrada_analogica1); //Lee el valor de la entrada analógica A0

valor_entrada_analogica2 = analogRead(entrada_analogica2); //Lee el valor de la entrada analógica A1
V1= (valor_entrada_analogica1*5*0.0009775171); // Convierte s un valor de entre 0 y 1023,
V2= (valor_entrada_analogica2*5*0.0009775171); //a otro de entre 0 y 5 para poder hacer el cálculo de la u.
u= k1*V1 +k2*V2 +k0; //Ecuación control óptimo de 0 a 5 volts
u=u*255*0.2; // Convierte el valor de U de 0 a 5V a un rango de entre 0 y 255,se verá en la salida el valor en pwm de U
U=round(u);
if(U<0) { // Si la U es menor a 0 el valor en pwm mínimo,que podrá obtener a la salida es de 0 volts.

U=0;
}
else if(U>255) { // Si la U es mayor a 5, el valor en pwm máximo que podrá obtener a la salida es de 5 volts.
U=255;
}
else { // Si el valor de U está entre 0 y 5 (incluyendolos),
}
analogWrite(salida1,U);
//Serial.println(U); //señal de control

//Serial.println(V1); //Estado1 voltaje Cap1
Serial.println(V2); //Estado2 voltaje Cap2
//delay(50);
}
61
Como puede verse en la Figura 6.5, se alimenta al circuito RC mediante el pin 6, que es la salida
digital PWM elegida como la señal de control u y aterrizándola mediante el pin Ground de la
placa Arduino. Se miden los voltajes x1 y x2 mediante el pin A0 y A1 de las entradas analógicas.
Analizando los diagramas de las figuras 6.1 y 6.5, se puede calcular la impedancia mínima del
circuito como 1 k𝛺, que para un voltaje máximo de 5V que brinda la salida PWM de Arduino,
entonces la corriente máxima que se le demanda a Arduino es 5 mA. En la figura 6.6 se muestra
el circuito implementado físicamente.
Figura 6.6. Circuito con la tarjeta Arduino.
Los resultados de la implementación experimental electrónica pueden verse en las Figuras 6.7 a
6.9.
Figura 6.7. Respuesta del estado x1.
62
Figura 6.8. Respuesta del estado x2.
Figura 6.9. Señal de control óptima u(t) de 0 a 255 para salida PWM de Arduino.
63
Caso de motor de CD
Considere un motor de CD cuyos diagramas esquemáticos eléctrico y mecánico están en la

figura 6.10.
Figura 6.10. Diagrama esquemático de un motor de CD.
La ecuación de la parte eléctrica es

dia
K1ev  Ra ia  La E
dt
donde ev es el voltaje de alimentación al motor, 𝑖 es la corriente en la armadura, E representa la
fuerza contraelectromotriz; K1, 𝑅 y 𝐿 son parámetros eléctricos constantes del motor. Las
ecuaciones de la parte mecánica son
d m 
Jm  T  bmm  L
dt c
m  c
d

dt
donde 𝜔 es la velocidad angular de la flecha del motor antes del engranaje reductor, 𝜔 es la
velocidad angular de la flecha a la salida del engranaje reductor, 𝜃 es el desplazamiento angular
a la salida del engranaje reductor, 𝑇 es el par de torsión de la flecha del motor antes del
engranaje reductor, 𝜏 el par causado por la carga (a la salida del engranaje reductor); 𝐽 , 𝑏 y 𝑐
son parámetros constantes del motor. Las ecuaciones del acoplamiento electromecánico son
E  K m
T  K ia
donde 𝑘 y 𝑘 son parámetros constantes del motor. Supóngase que la carga  L es modelada
d
como  L  J L  bL , donde 𝐽 y 𝑏 son coeficientes de inercia y fricción constantes en la
dt
carga.
64
Considérese que el circuito electrónico para control del motor CD fuese mediante un transistor
BJT NPN como el de la figura 6.11, donde RB es la resistencia de base, iC es la corriente de
colector, Vbe es la tensión de caída constante entre base y emisor, VCC es un voltaje constante de
alimentación y PWM es el voltaje de control mediante modulación por ancho de pulso
(expresada de 0 a 255 en Arduino equivalente a voltaje ‘analógico’ de salida entre 0V y V Amax).
Figura 6.11. Diagrama esquemático del circuito de control.
Obtenga un modelo de estados de todo el sistema si las variables de estado son x1=θ y x2=ω, el
PWM es la entrada y y=ω es la salida. Luego, calcule e implemente el control óptimo u* tal que
la salida gire a las velocidades angulares de la tabla 6.1.
Tabla 6.1. Los valores deseados para la velocidad del motor

Tiempo (segundos) Desplazamiento angular deseado r1 Velocidad angular
(rev) deseada r2 (rpm)
0<t<10 30t/60 30
10<t<20 60(t-10)/60+30(10)/60 60
20<t<30 90(t-20)/60+(60+30)(10)/60 90
30<t<40 120(t-30)/60+(90+60+30)(10)/60 120
40<t<50 150(t-40)/60+(120+90+60+30)(10)/60 150
Solución
La corriente de base depende del voltaje de PWM y la caída de voltaje entre base y emisor por la
expresión
𝑉Amax
⎧ 255 𝑃𝑊𝑀 − 𝑉 𝑉Amax
⎪ , si 𝑃𝑊𝑀 > 𝑉
𝑖 = 𝑅 255
⎨ 𝑉Amax
⎪ 0, si 𝑃𝑊𝑀 ≤ 𝑉
⎩ 255
65
La corriente de colector toma un valor amplificado de la corriente de base con ganancia 𝛽, esto
es
𝑖 = 𝛽𝑖 .
Por lo que, 𝑖 = Amax

𝑃𝑊𝑀 − 𝑉 si Amax
𝑃𝑊𝑀 > 𝑉 . Un modelo en espacio de estado del
motor de CD se obtiene bajo la teoría de perturbaciones singulares (ver Kelly et al. (2005) y
Khalil (1996)), ya que partiendo de que La  0 la dinámica del subsistema eléctrico es mucho
más rápida (aparentemente instantánea) que el resto de la dinámica del sistema. Por lo que la
corriente de armadura se aproxima por una relación algebraica
K1ev  Ra ia  E
K1ev  Ra ia  Km
1
ia   K1ev  Km 
Ra
Sin embargo, en este circuito se puede ver que la corriente de armadura es controlada
aproximadamente como la corriente de colector; esto es,
𝑖 =𝑖
Por lo que se deduce:

d m 
Jm  K ia  bmm  L
dt c
d J d  bL
cJ m  K iC  cbm  L  
dt c dt c
d d
c2 J m  cK iC  c 2bm  J L  bL
dt dt
d
c J
2
m  JL 
dt
  c 2bm  bL    cK iC
 c 2 J m  J L  ddt   c 2bm  bL    cKR   V255

A max 
PWM  Vbe  , si 𝑃𝑊𝑀 >
 Amax
𝑉
B
d 
 2
c 2bm  bL 

cK   VA max 
PWM  Vbe 

dt  c J m  J L  RB  c J m  J L   255
2

Entonces, el modelo en espacio de estados es

0 1   0 
 x1     x1   
 x    1    Ks u
 2  0    x2   
 Ts   Ts 
1 c 2bm  bL Ks cK VA max 𝑃𝑊𝑀 − 𝑉 , si 𝑃𝑊𝑀 > 𝑉
donde  2 ,  y 𝑢 = Amax Amax
.
Ts c J m  J L Ts 255RB  c 2 J m  J L  0, si 𝑃𝑊𝑀 ≤ 𝑉
Amax
66
Los parámetros KS y TS se pueden estimar a partir de un experimento en el que se aplica una

entrada escalón constante 𝑃𝑊𝑀 > 𝑉 , midiendo la velocidad x2=ω(t) y desplazamiento
Amax
angular x1=θ(t) de la flecha hasta un tiempo t en que el motor está en régimen estacionario. Esto
es,
x2 (t ) x 2 (t )
K S  lim  lim
t  u t  255Vbe
PWM 
VA max
 255Vbe 
t K S ·t· PWM    x1 (t )  x1 (0)
K S ·u·t   x2 ( )d   VA max 
TS  lim 0
 lim
t  x 2 (t ) t  x2 ( t )
Como ejemplo, se realiza un experimento con un circuito y motor de CD como en la figura 6.15,
aplicando un PWM=255 se obtiene una curva de respuesta como la figura 6.12. En el Apéndice
se explica cómo configurar Proteus y Arduino para poder adquirir las señales (tiempo,
desplazamiento y velocidad angular) en una hoja de cálculo como MS Excel.
Figura 6.12. Curva de velocidad angular de la flecha de un motor de CD ante entrada escalón.
Tabla 6.2. Parámetros clave de la curva de respuesta del motor de CD de la figura 6.12
Algún tiempo tf en que está Velocidad angular Desplazamiento angular de la
en régimen estacionario de la flecha en tf flecha en tf
50 s 155 rpm 114.8 rev
Usando los valores de la tabla 6.2 y conociendo que 𝑉Amax = 4.97𝑉 y 𝑉 = 0.91𝑉 en el circuito
de la figura 6.15, se calcula que KS=0.744084 [rpm/bits] y TS=0.092688 min. En la figura 6.12 se
ilustra una validación del modelo propuesto para los valores calculados de KS y TS.
67
Por otro lado, para lograr el objetivo de optimización de seguimiento de velocidad angular, se
propone minimizar el índice de costo
0.5( x  r )  1000( x2  r2 ) 2  0.5u 2 dt
tf
J  1 1
2
t0
Entonces las matrices son

0 1   0 
A  1  , B   0.744084  , Q  1 0 
, R  1, H 
0 0 
0      0 2000  0 0 
   
 0.092688   0.092688 
Aplicando la metodología de control óptimo para seguimiento, similar al caso 2 de la sección

5.4, pero sumando el término 𝑟 /𝐾 al control para que la dinámica de error (1.32) sea
asintóticamente estable en el origen, el controlador a implementar es
0.744084 r2 r2
u
0.092688
 k12  x1  r1   k 22  x2  r2   
0.744084
 k1  x1  r1   k 2  x2  r2  
0.744084
con k12 y k22 son soluciones de la ecuación algebraica de Riccati, para este ejemplo resulta
k1  1, k 2  43.400398 . Luego, el valor del PWM se calcula por 𝑃𝑊𝑀 = 𝑢 + 𝑉 .
Amax
El comportamiento del motor se puede simular numéricamente en Scilab para k (t f )   0 0; 0 0

x(t0 )   0.03;0 y t f  50 s. La respuesta de los estados en el tiempo está en la figura 6.13 y la
señal de control en la figura 6.14, donde se ve que por momentos el control se satura (no puede
ser mayor) y por ello tarda más en lograr el seguimiento.
68
Figura 6.15. Implementación en Proteus.
69
Para ilustrar los resultados con el motor de CD, se implementa un circuito de control como se
muestra en la figura 6.15, donde se requiere conocer el 𝑉Amax y 𝑉 . En este ejemplo fueron
𝑉Amax = 4.97𝑉, 𝑉 = 0.91𝑉 para el transistor utilizado, resolución del encoder de 100pasos/rev
en cualquiera de las terminales del encoder nombradas como Q1 o Q2, la terminal IDX es para
transmisión de datos. Para el caso de simulación en Proteus se configuraron los siguientes:
Voltaje nominal del motor de 12 V, resistencia de armadura de 100𝛺, inductancia de armadura
de 100mH, 𝑅B = 2.4 kΩ, ganancia 𝛽 = hfe = 30. Finalmente, se programa el control óptimo en
Arduino, mediante un código como el siguiente:
#include <FreqCount.h>
int pin2=9;
float RPM;
float Des=0;
float T;
float time;
float u1=0;
float u2=0;
float u3=0;
float r1=0;
float r2=0;
float k1=-1; //ganancia óptima Kb(1)
float k2=-43.400398; //ganancia óptima Kb(2)
float Ks=0.7918; //valor de Ks del modelo
float Vbe=0.91; //Voltaje de base-emisor en el transistor
float Vamax=4.97; //voltaje de pin en Arduino para PWM=255
float R=30; //incremento de velocidad deseada de cada intervalo
float A=10; //tiempo en segundos de cada velocidad deseada
void setup() {
Serial.begin(9600);
FreqCount.begin(100);
pinMode(pin2, OUTPUT);
}
void loop() {
time = millis();//tiempo instantaneo en milisegundos
T= time/1000; //tiempo en segundos
if(T<A){
r1=R*(T-0)/60+0;
r2=R;
} else if (T<2*A){
r1=2*R*(T-A)/60+R*A/60;
r2=2*R;
} else if (T<3*A){
r1=3*R*(T-2*A)/60+(1+2)*R*A/60;
r2=3*R;
} else if (T<4*A){
r1=4*R*(T-3*A)/60+(1+2+3)*R*A/60;
r2=4*R;
} else if (T>4*A){
r1=5*R*(T-4*A)/60+(1+2+3+4)*R*A/60;
r2=5*R;
}
if (FreqCount.available()) {
float count = FreqCount.read();//count es el incremento en rev cada 10 seg
RPM= count*6;
Des = Des + count*0.1*0.1;//el muestreo es cada 0.1 segundo
u1=k1*(Des-r1)+k2*(RPM-r2)+r2/Ks; //controlador óptimo
u2=u1+255*Vbe/Vamax;
u3=round(u2);
if(u3<0){
u3=0;
70
}
else if(u3>255){
u3=255;
}
else{
}
analogWrite(pin2,u3);
//Serial.print(" ");
Serial.print(T);
//Serial.print(" Des = ");
Serial.print(",");
Serial.print(Des,3);
Serial.print(",");
Serial.print(RPM);
//Serial.print(" u3 = "); //señal de control
Serial.print(",");
Serial.print(u3,3);
Serial.println();
}
}
En el Apéndice se explica cómo configurar Proteus y Arduino para poder adquirir las señales
(tiempo, desplazamiento y velocidad angular) en una hoja de cálculo como MS Excel. La
respuesta del motor de Proteus en lazo cerrado con el control óptimo implementado en Arduino
se muestra en la figura 6.16, donde se aprecia que la velocidad angular (rpm) vs el tiempo (s)
cumple el objetivo de la tabla 1.
velocidad (rpm) vs tiempo (s)

160
140
120
100
80
60
40
20
0
0 5 10 15 20 25 30 35 40 45 50
Figura 6.16. Velocidad controlada en el motor de CD.
71
Caso de control de nivel en tanques
Considere el sistema de nivel de tanques de la Figura 6.17.
Figura 6.17. Sistema de 3 tanques de nivel.
Cuyas ecuaciones son

dh1
A1  Q1  Q13
dt
dh
A2 2  Q2  Q32  Q20
dt
dh
A3 3  Q13  Q32
dt
Q13  az1S n sign( h1  h3 ) 2 g h1  h3
Q32  az 3 Sn sign( h3  h2 ) 2 g h3  h2
Q20  az 2 S n 2 gh2
donde A1, A2 y A3 son las áreas de los tanques, h1, h2 y h3 son los niveles medibles de líquido en
cada tanque, Q1 y Q2 son los caudales de líquido de entrada en los tanques 1 y 2, Sn es el área de
apertura nominal de conexión entre tanques, az1, az2 y az3 son factores de corrección de apertura
entre 0 y 1, g es la constante de aceleración gravitacional. Si los modelos de caudales de
interconexión se linealizan alrededor de un punto de operación, entonces
72
h1  h3
Q13 
R13
h3  h2
Q32 
R32
h2
Q20 
R20
donde R13, R32 y R20 son constantes de oposición al flujo de corriente, calculadas al linealizar.
Calcule e implemente un control óptimo u* si se desea que el primer tanque mantenga un nivel
constante de 300 mm de altura, el tanque 2 un nivel constante de 100 mm de altura y el tanque 3
un nivel constante de 200 mm de altura.
Solución
El criterio de costo se propone como
10 10 [h  0.30]  10  106 [ h2  0.10]2  10  106 [ h3  0.20]2  u12  u2 2 dt
tf
J  6
1
2
t0
Si las variables de estado son x1=h1, x2=h2 y x3=h3, las entradas u1=Q1/KS y u2=Q2/KS, donde KS es la
relación de voltaje de entrada a las bombas con el caudal de bombeo de entrada. Entonces un
modelo en espacio de estado en el caso linealizado es obtenido como
dh1 h h h h
A1  Q1  1 3  Q1  1  3
dt R13 R13 R13
dh2 h h h h h h
A2  Q2  3 2  2  Q2  3  2  2
dt R32 R20 R32 R32 R20
dh3 h1  h3 h3  h2 h1 h h h
A3     3  3  2
dt R13 R32 R13 R13 R32 R32
 1 1 
 0   KS 
A 0 
 R13 A1 R13 A1

 x1     x1   1 
 1 1  K S   u1 
 x    0   
1   x2    0
 2 
 R32 A2 R20 A2  R32 A2    A2  u2 
 x3     x3   
 1 1  1 1   0 0 
 R A     
R32 A3  R13 A3 R32 A3    
 13 3
Supóngase que los parámetros de la planta son Sn=5x10-5 s/m2, A₁=A₂=A3=0.0154 m2, R13=2850
s/m2, R32=2020 s/m2, R20=3500 s/m2, KS=100/(5x106) m3/(s V), entonces las matrices son
73
 0.0227842 0 0.0227842  0.0012987 0 


A 0  
0.0506989 0.0321461  , B   0 0.0012987  ,
 0.0227842 0.0321461 0.0549303  0 0 
 20 106 0 0  0 0 0   0.30 
  2 0
Q 0 20 10 6
0 ,R    , H  0 0 0  , r   0.10
 
 0 0 2
 0 20 106  0 0 0   0.20 

u1 (t )  0.00064935  k11 x1  k12 x2  k13 x3  s1   k1 x1  k2 x2  k3 x3  k7
u2 (t )  0.00064935  k12 x1  k22 x2  k23 x3  s2   k2 x1  k4 x2  k5 x3  k8
con k11, k12, k13, k22, k23, s1 y s2 como soluciones de la ecuación de Riccati.
Este controlador se puede simular computacionalmente para conocer su desempeño en la

planta no lineal, con k (t f )   0 0 0;0 0 0; 0 0 0 , x(t0 )   0; 0;0 y tf =240 s=4 min. La solución
de la ecuación de Riccati se ilustra en la Figura 6.18. La respuesta de los estados en el tiempo de
la planta no lineal en la Figura 6.19 y la señal de control en la Figura 6.20.
74
Figura 6.19. La respuesta de los estados x1, x2 y x3.
Figura 6.20. Señales de control óptimas u1(t) y u2(t).
75
Para la implementación experimental de este ejemplo se cuenta con el equipo experimental de 3

tanques, que se muestra en la figura 6.21. El cual tiene un módulo de potencia (que se muestra
en la figura 6.22) que permite adquirir las señales de sensores de nivel para cada tanque y
proporcionarle voltajes para control de las bombas hidráulicas.
Figura 6.21. Equipo experimental de 3 tanques.
Figura 6.22. Módulo de potencia del equipo experimental.
Se propone programar el controlador para este equipo mediante la tarjeta Arduino UNO,
usando algunos circuitos para acondicionamiento de las señales eléctricas entre el módulo de
potencia y la tarjeta Arduino, como en la figura 6.23. Los voltajes de salida de los sensores de
nivel y de entrada a las bombas hidráulicas del módulo de potencia están acotados entre +10 V
y -10 V, mientras que las entradas analógicas y las salidas PWM en Arduino están acotadas
entre 0 V y 5 V.
76
Figura 6.23. Esquema de control del equipo experimental.
Para el acondicionamiento de las señales de cada sensor se utiliza el circuito de la figura 6.24.
Figura 6.24 Diagrama esquemático del circuito de adecuación de señal para cada sensor de
nivel.
A partir de la figura 6.24, puede notarse que el voltaje de salida está dado por
𝑅 𝑅 𝑅
𝑉 = 𝑉𝑖 − 𝑉
𝑅 𝑅 𝑅
77
donde VS es el voltaje del sensor, Vi es un voltaje constante (puede leerse de un pin PWM de
Arduino) y el Vsa es el voltaje hacia un pin de entrada analógica en Arduino. Se elige Vi=2.5V y
( )
las resistencias tales que =1y = , con (𝑉 −𝑉 ) como el intervalo de
( )
voltajes de mínimo a máximo del sensor y (ℎ − ℎ ) como el intervalo de nivel de agua en
el tanque asociado al intervalo de voltajes de mínimo a máximo del sensor, para que 0mm de
nivel de agua equivalga a 0V en Arduino y 620mm de nivel del agua a 5V en Arduino. Este
mismo circuito se implementa para el acondicionamiento de señal para cada uno de los 3
sensores, como se ilustra en la figura 6.25.
Figura 6.25 Circuito electrónico de adecuación de señal para los 3 sensores de nivel.
Para el acondicionamiento de las señales de control de las bombas hidráulicas se utiliza el

circuito de la figura 6.26.
78
Figura 6.26 Diagrama del circuito de acondicionamiento de señal para cada bomba hidráulica.
A partir de la figura 6.26, puede notarse que el voltaje de salida está dado por
𝑅 𝑅 𝑅
𝑉 = 𝑉 − 𝑉
𝑅 𝑅 𝑅
donde VB es el voltaje de entrada a la bomba hidráulica, Vi es un voltaje constante (puede leerse

de un pin PWM de Arduino) y el Voa es el voltaje desde un pin de salida PWM de Arduino. Se
( )
elige Vi=2.5V y las resistencias tales que = = , con (𝑉 −𝑉 ) como el
intervalo de voltajes de mínimo a máximo de entrada a las bombas hidráulicas, para que 0ml/s
de caudal de agua bombeada equivalga a 0V en Arduino y 100ml/s de caudal a 5V en Arduino.
Este mismo circuito se implementa para el acondicionamiento de señal para cada una de las 2
bombas, como se ilustra en la figura 6.27.
Figura 6.27 Circuito electrónico de acondicionamiento de señal para cada bomba hidráulica.
79
En Arduino se programan los controladores u1 (t )  k1 x1  k2 x2  k3 x3  k7 y

u2 (t )  k2 x1  k4 x2  k5 x3  k8 , con las ganancias de la ecuación algebraica de Riccati (para
estado estacionario), que fueron calculadas como k1=3148.0924, k2=4.6692496, k3=597.56158,
k4=3130.0672, k5=842.8796, k7=1074.4324, k8=493.63218; mediante la siguiente codificación del
algoritmo de control:
#include <math.h>
unsigned long t0;

float t;
int entrada_analogica1= A0; // declara numeros de pines de entrada de tipo entero
int valor_entrada_analogica1=0; // Declara valores entradas
int salida1 = 5; // declara numeros de pines de salida PWM tipo entero
int salida2 = 6; // declara numeros de pines de salida PWM tipo entero
float h1;
float h2;
float h3;
float u1=0;
float u2=0;
int U1=0;
int U2=0;
float k1=3148.0924;
float k2=4.6692496;
float k3=597.56158;
float k4=3130.0672;
float k5=842.8796;
float k7=1074.4324;
float k8=493.63218;
void setup() {
pinMode(entrada_analogica1, INPUT); // Pines de salida y entrada

pinMode(salida1, OUTPUT);
pinMode(salida2, OUTPUT);
Serial.begin(9600); // Inicializar la comunicación serial a 9600 bits por segundo
}
void loop() { //Entra a un ciclo repetitivo
valor_entrada_analogica1= analogRead(entrada_analogica1); //Lee el valor de la entrada analógica A0

h1= (valor_entrada_analogica1*0.0009775171*0.62); // Convierte un valor de entre 0 y 1023,
h2= (valor_entrada_analogica2*0.0009775171*0.62); //a otro de entre 0 y 0.62m para poder hacer cálculos.
h3= (valor_entrada_analogica3*0.0009775171*0.62);
u1= -k1*h1 -k2*h2 -k3*h3 +k7; //Ecuación control óptimo en Volts

u2= -k2*h1 -k4*h2 -k5*h3 +k8; //Ecuación control óptimo en Volts
u1=u1*255*0.2; // Convierte el valor de U de 0 a 5V a un rango de entre 0 y 255,se verá en la salida el valor en pwm de U
u2=u2*255*0.2; // Convierte el valor de U de 0 a 5V a un rango de entre 0 y 255,se verá en la salida el valor en pwm de U
U1=round(u1);
U2=round(u2);
if(U1<0) { // Si la U es menor a 0 el valor en pwm mínimo,que podrá obtener a la salida es de 0 volts.

U1=0;
}
else if(U1>255) { // Si la U es mayor a 5, el valor en pwm máximo que podrá obtener a la salida es de 5 volts.
U1=255;
}
80

}
if(U2<0) { // Si la U es menor a 0 el valor en pwm mínimo,que podrá obtener a la salida es de 0 volts.

U2=0;
}
else if(U2>255) { // Si la U es mayor a 5, el valor en pwm máximo que podrá obtener a la salida es de 5 volts.
U2=255;
}
}
analogWrite(salida1,U1);
analogWrite(salida2,U2);
Serial.println(U1); //señal de control bomba 1

//Serial.println(U2); //señal de control bomba 2
//delay(50);
}
Los resultados de la implementación experimental pueden verse en las Figuras 6.28.
¿FIGURA?
Figura 6.28. La respuesta de los estados x1, x2 y x3.
81
Caso de robot de 2 grados de libertad
Figura. 6.29. Robot de 2 grados de libertad.
Considérese el robot planar de 2 eslabones rotacionales con estructura mecánica como se ilustra
en la figura 6.29, descrito en la página 15 del Cap. 10 del libro Mecatrónica: control y
automatización de J.F. Reyes Cortés et al., donde se indican los valores de sus parámetros físicos,
así como los elementos y matrices M ( x), C  x  , G  x  , F  x  de su modelo dinámico, el cual
está dado por
 x1   x3 
 x    x 
 2  4
,
 x3 
 x   M ( x )   H  x   G  x   F  x 
1
 4
donde x1 y x2 son las posiciones articulares y x3 y x4 son las velocidades articulares,
respectivamente; τ =[τ1 τ2]T es el vector de pares de entrada en cada articulación.
Puede ilustrarse el desempeño de dicho robot en lazo cerrado con el controlador óptimo, para
cumplir el objetivo de regulación de posición articular en r1  120º, r2  30º , con gasto mínimo
de energía, mediante el siguiente procedimiento:
Puesto que el modelo no es lineal, primero se usa una ley de control de par calculado dada por
𝜏 = 𝐻(𝑥) + 𝐺(𝑥) + 𝐹(𝑥) + 𝑀(𝑥)𝑢
donde u=[u1 u2]T es un vector de controles auxiliares óptimos para el sistema linealizado de lazo
cerrado. Así, la ecuación de lazo cerrado es lineal dada por
82
𝑥̇ 0 0 1 0 𝑥 0 0
𝑥̇ 𝑥 0 𝑢
= 0 0 0 1
𝑥 + 0
𝑥̇ 0 0 0 0 1 0 𝑢
𝑥̇ 0 0 0 0 𝑥 0 1
El funcional de costo se propone como

tf
 1 1 
J    x1 (t )  r1    x2 (t )  r2   u12 (t )  u22 (t )  dt
2 2
t0 
2 2 
Por lo que las matrices para el cálculo de ganancias óptimas son

0 0 1 0 0 0 0 0 0 0 2 0 0 0
𝐴= 0 0 0 1 ,𝐵 = 0 0 , 𝐻 = 0 0 0 0 ,𝑄 = 0 2 0 0 ,𝑅 = 1 0 ,
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 1 0 0 0 0 0 0 0 0
120°
𝑟(𝑡) = 30°
0
0
Ya que este ejemplo corresponde al caso 2 del Capítulo 5 (donde A(t )r (t )  0 y r (t )  0 ),

entonces el control auxiliar óptimo se puede calcular como
u *(t )   R 1 (t ) BT (t ) K (t )[ x(t )  r (t )]
 k11 k12 k13 k14    x1  120º  
  
1 0  0 0 1 0  k21 k22 k23 k24    x2   30º  
    
  k34    x3   0  
0 1  0 0 0 1   k31 k32 k33
    
 k41 k42 k43 k44    x4   0  
  x1  120º  
 
u 
*
 k31 k32 k33 k34    x2   30º  
   
1
  
*
u   k41 k42 k43 k44    x3   0  
2
     
  x4   0  
u1*  k31  x1  120º   k32  x2  30º   k33 x3  k34 x4
u2*  k41  x1  120º   k42  x2  30º   k43 x3  k44 x4
Puesto que se cumplen las condiciones de ganancias óptimas en horizonte infinito, k31, k32, k33,
k34, k41, k42, k43, k44 pueden calcularse como soluciones de la ecuación algebraica de Riccati:
0 = −𝐾(𝑡)𝐴(𝑡) − 𝐴 (𝑡)𝐾(𝑡) − 𝑄(𝑡) + 𝐾(𝑡)𝐵(𝑡)𝑅 (𝑡)𝐵 (𝑡)𝐾(𝑡)
83
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
𝑘 𝑘 𝑘 𝑘 0 0 1 0
𝑘 𝑘 𝑘 𝑘 0 0 0 1
=−
𝑘 𝑘 𝑘 𝑘 0 0 0 0
𝑘 𝑘 𝑘 𝑘 0 0 0 0
0 0 0 0 𝑘 𝑘 𝑘 𝑘 2 0 0 0
𝑘 𝑘 𝑘 𝑘
− 0 0 0 0 − 0 2 0 0
1 0 0 0 𝑘 𝑘 𝑘 𝑘 0 0 0 0
0 1 0 0 𝑘 𝑘 𝑘 𝑘 0 0 0 0
𝑘 𝑘 𝑘 𝑘 0 0 𝑘 𝑘 𝑘 𝑘
𝑘 𝑘 𝑘 𝑘 0 0 0 0 1 0 𝑘 𝑘 𝑘 𝑘
+
𝑘 𝑘 𝑘 𝑘 1 0 0 0 0 1 𝑘 𝑘 𝑘 𝑘
𝑘 𝑘 𝑘 𝑘 0 1 𝑘 𝑘 𝑘 𝑘
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 𝑘 𝑘 0 0 0 0 2 0 0 0
0 0 𝑘 𝑘 0 0 0 0 0 2 0 0
=− − 𝑘 𝑘 𝑘 𝑘 −
0 0 𝑘 𝑘 0 0 0 0
0 0 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘 0 0 0 0
𝑘 𝑘 +𝑘 𝑘 𝑘 𝑘 +𝑘 𝑘 𝑘 𝑘 +𝑘 𝑘 𝑘 𝑘 +𝑘 𝑘
+
Se obtiene
0 = −2 + 𝑘 +𝑘
0=𝑘 𝑘 +𝑘 𝑘
0 = −𝑘 + 𝑘 𝑘 + 𝑘 𝑘
0 = −𝑘 + 𝑘 𝑘 + 𝑘 𝑘
0 = −2 + 𝑘 +𝑘
0 = −𝑘 + 𝑘 𝑘 + 𝑘 𝑘
0 = −𝑘 + 𝑘 𝑘 + 𝑘 𝑘
0 = −𝑘 − 𝑘 + 𝑘 𝑘 + 𝑘 𝑘
0 = −𝑘 − 𝑘 + 𝑘 𝑘 + 𝑘 𝑘
0 = −𝑘 − 𝑘 + 𝑘 +𝑘
Las cuales se puede resolver computacionalmente y simular la respuesta del robot con
x(t0 )  0º ;0º ; 0 sº ; 0 sº  y tf =10 s, utilizando código en Scilab como el siguiente:
84
clear;clc //limpia memoria y pantalla

//----ESPECIFICACIONES DEL OBJETIVO DE CONTROL ------
Qd1 = 120; //posición articular deseada 1 [en º]
Qd2 = 30; //posición articular deseada 2 [en º]
//----PARÁMETROS DE SIMULACIÓN ------

t0=0;//tiempo inicial [s]
tf=10; //tiempo final de simulación [s]
paso=0.1;//paso de integración [s]
t=(t0:paso:tf)';//vector de tiempo
x0=[0;0;0;0;0]*%pi/180;//condiciones iniciales=[pos art1,pos art2,vel art1,vel art2,J]
//----PARÁMETROS FÍSICOS DEL ROBOT DE 2 GDL ------

g=9.81;//gravedad m/s^2
m1=23.902;//masa kg
m2=3.88;//masa kg
l1 = 0.45+%eps;//longitud m
l2 = 0.45;//longitud m
lc1 = 0.091;//centro de masa m
lc2 = 0.048;//centro de masa m
I1=1.266;//inercia N-m·s^/rad
I2=0.093;//inercia N-m·s^/rad
b1=2.288; //coef fricción N-m·s/rad
b2=0.175; //coef fricción N-m·s/rad
//----ELEMENTOS PARA CONTROL ÓPTIMO------

A=[0,0,1,0;0,0,0,1;0,0,0,0;0,0,0,0]; //Matriz de estados A
B=[0,0;0,0;1,0;0,1]; //Matriz de entrada B
Q=[2,0,0,0;0,2,0,0;0,0,0,0;0,0,0,0]; //Matriz de seguimiento del funcional de costo
R=[1,0;0,1]; //Matriz de entrada del funcional de costo
//K1 = -[1,0;0,1]; K2 = -[2,0;0,2]; Kb=[K1,K2]; //Ganancias arbitrarias
//----SIMULACIÓN DE MODELO Y CONTROL EN LAZO CERRADO-----

function dx=robot2gdl(t,x)
q_d1=Qd1*%pi/180; //posición deseada en rad
q_d2=Qd2*%pi/180; //posición deseada en rad
qd = [q_d1;q_d2]; //vector de posiciones articulares deseados
q1=x(1); q2=x(2);//posiciones articulares reales
qp1=x(3); qp2=x(4);//velocidades articulares reales
q = [q1;q2]; //vector de posiciones articular q˙=[q1 q2]T
qp = [qp1;qp2]; //vector de velocidad articular q˙=[q˙1 q˙2]T
qt = q-qd; //Vector de error de posición
//MODELO DINÁMICO DEL ROBOT:
M=[m1*lc1^2+m2*l1^2+m2*lc2^2+2*m2*l1*lc2^2*cos(q2)+I1+I2,
m2*lc2^2+m2*l1*lc2*cos(q2)+I2; //matriz de inercias M(q)
m2*lc2^2+m2*l1*lc2*cos(q2)+I2, m2*lc2^2+I2];
C=[-2*m2*l1*lc2*sin(q2)*qp2, -m2*l1*lc2*sin(q2)*qp2; //matriz de fuerzas centrífugas y
de Coriolis C(q,q˙)
m2*l1*lc2*sin(q2)*qp1, 0];
G=g*[(lc1*m1+m2*l1)*sin(q1)+m2*lc2*sin(q1+q2); //vector de efectos gravitacionales
g(q)
m2*lc2*sin(q1+q2)];
Fr=[b1,0;0,b2]*[qp1;qp2]; //Vector de pares de fricción viscosa
//ELEMENTOS PARA EL CONTROLADOR:

u=Kb*[qt;qp]; //=-R\B'*K*[qt;qp]; Controlador auxiliar para el sistema linealizado
tau = M*(u) + G + C*qp +Fr; //Ley de control por par calculado
qpp = M\(tau-C*qp-G-Fr); //Ec. dif ord en lazo cerrado

Jp=[qt;qp]'*Q*[qt;qp] +u'*R*u; //integrando del funcional de costo
dx = [qp; qpp; Jp];//Ec. en espacio de estado
85
endfunction
x=(ode('rk',x0,t0,t,robot2gdl))';//comando para solución de Ec de estado

q1=x(:,1);q2=x(:,2);qp1=x(:,3);qp2=x(:,4);J=x(:,5);
//-----COORDENADAS DE CINEMÁTICA PARA DIBUJO DE ROBOT EN CARTESIANAS

l1=l1*1000;l2=l2*1000;//para expresar las cartesianas en mm
NT=size(t,1);
//coordenadas de la base:
bx = zeros(NT,1);
by = zeros(NT,1);
bz = zeros(NT,1);
//coordenadas del codo:
cx = l1*sin(q1);
cy = zeros(NT,1);
cz = -l1*cos(q1);
//coordenadas del extremo:
ex = l1*sin(q1)+l2*sin(q1+q2);
ey = zeros(NT,1);
ez = -l1*cos(q1)-l2*cos(q1+q2);
//Puntos de unión de cada eslabón
rx=[bx,cx,ex];
ry=[by,cy,ey];
rz=[bz,cz,ez];
//coordenadas deseadas en cartesianas del extremo:
Px = l1*sind(Qd1)+l2*sind(Qd1+Qd2);
Pz = -l1*cosd(Qd1)-l2*cosd(Qd1+Qd2);
//-----MUESTRA RESULTADOS DE OBJETIVO Y CONTROL DE POSICIÓN------

mprintf('\n\nLas posiciones articulares deseadas (en grados sexagesimales)
son:\nq_d1=%f \nq_d2=%f',Qd1,Qd2)
mprintf('\n\nLas posiciones articulares reales alcanzadas (en grados sexagesimales)
son:\nq1=%f \nq2=%f',q1($)*180/%pi,q2($)*180/%pi)
mprintf('\n\nEl valor final del funcional de costo es:\nJ(t_f)=%f',J($))
//-----GRÁFICAS DEL DESEMPEÑO EN LAZO CERRADO------

eboxC=[-l1-l2,l1+l2,-l1-l2,l1+l2,0,0];//área de gráfica en cartesianas
figure
for j=1:NT//cantidad de muestras de tiempo
drawlater()
clf;//limpia figura
subplot(2,2,1);
plot(t(1:j),q1(1:j)*180/%pi,t(1:j),q2(1:j)*180/%pi,t,Qd1,'m--',t,Qd2,'c--
',t(1:j),qp1(1:j)*180/%pi,t(1:j),qp2(1:j)*180/%pi); //Gráfica de posiciones
articulares vs tiempo
xgrid//dibuja malla
p = get("hdl");p.children.thickness = 2;//Grosor de línea
xtitle('Movimiento de articulaciones','tiempo t [s]','Coordenadas articulares
[grados]') //título
legend('q1(t)','q2(t)','qd1(t)','qd2(t)','$\dot{q}_1$','$\dot{q}_2$',opt=-5)
//etiquetas
subplot(2,2,2);
plot(rx(j,:),rz(j,:),'b',ex(1:j),ez(1:j),'go',Px,Pz,'rx')//posiciones del robot en
coordenadas cartesianas
xgrid//dibuja malla
P=gca();P.data_bounds=eboxC;P = get("hdl");P.children.thickness = 3;//grosor de línea
xtitle('Movimiento del robot en espacio cartesiano, en t='+msprintf('%3.1f',t(j))+'
s','eje X [m]','eje Z [m]') //título
legend('Robot','Posiciones reales de TCP','Posiciones deseadas de TCP',opt=2)
//etiquetas
subplot(2,2,3);
86
plot(t(1:j),Kb(1,1)*(q1(1:j)-Qd1*%pi/180)+Kb(1,2)*(q2(1:j)-
Qd2*%pi/180)+Kb(1,3)*qp1(1:j)+Kb(1,4)*qp2(1:j),t(1:j),Kb(2,1)*(q1(1:j)-
Qd1*%pi/180)+Kb(2,2)*(q2(1:j)-Qd2*%pi/180)+Kb(2,3)*qp1(1:j)+Kb(2,4)*qp2(1:j));
//gráfica de entradas de control
xgrid//dibuja malla
p = get("hdl");p.children.thickness = 2; //Grosor de línea
xtitle('Entradas auxiliares optimas','tiempo t [s]','u(t) [N-m]')//títulos
legend('u1','u2') //etiquetas
subplot(2,2,4);
plot(t(1:j),J(1:j)); //Gráfica del funcional de costo
xgrid//dibuja malla
p = get("hdl");p.children.thickness = 2; //Grosor de línea
xtitle('Funcional de costo','tiempo t [s]','J(t)')//títulos
drawnow()
end
En la Figura 6.30 se ilustra la respuesta en el tiempo de los estados, posiciones cartesianas, las
señales de control y funcional de costo.
Figura 6.30. La respuesta de los estados x1 a x4, posiciones cartesianas, señales de control óptima
u(t) y funcional de costo J(t).
87
Caso de intercepción de aeronave
Considérese una aeronave que inicialmente está en el punto { xA  d , yA  0 , z A  200 } y se

desea interceptar por un misil que inicialmente está en { xM  0 , yM  0 , zM  0 } con gasto
mínimo de energía. Se conoce que la aeronave se desplaza en una trayectoria { x A (t )  d  0.1t 3 ,
yA (t )  0 , zA (t )  200 }. Despreciando fuerzas gravitacionales y aerodinámicas, supóngase que el
modelo de movimiento del misil está dado por  xM  u , donde u es la aceleración por propulsión
 x (t )
 M 
y con zM  z A 1  e 4420  .
 
Para el modelo de estado considere que las variables de estado son el desplazamiento x1  xM y
velocidad x2  xM , la aceleración neta u es la entrada y xA (t ) es la referencia para seguir.
Resuelva numéricamente las ecuaciones diferenciales para obtener las ganancias K(t) y S(t)
mediante un programa computacional. Luego, implemente el control óptimo que logre el
objetivo.
Solución
Un modelo en espacio de estado del misil es
 x1   0 1   x1   0 
 x    0 0   x   1  u
 2   2  
Para lograr el objetivo de seguimiento con ahorro de energía se propone el índice de costo
[ x (t )  x A (t )]2  1000u 2 dt

tf
J ( x, u , t )  100[ xM (t f )  x A (t f )]2   M
t0
entonces las matrices son:
0 1  0  2 0  200 0   d  0.1t 3 
A  , B  1  , Q   0 0  , R   2000 , H   0 0  , r (t )   
0 0         r2 

1
 k12 x1  k 22 x2  s2 
u (t )  
2000
con k12, k22 y s2 como soluciones de la ecuación de Riccati. El cual se puede simular
computacionalmente para k (t f )   200 0;0 0  , x(t0 )   0;0 , d=1,000 m y t f  60 s, mediante
un código como el siguiente:
clear;clc;
//Datos de simulacion
t0=0;
tf=60;
paso=1;
ts=(t0:paso:tf)';
x0=[0;0;0];//condiciones iniciales de estado e indice de costo
d=1e3;
88
// Matrices de la funcion de costo J (Q y R):

H=[200 0; 0 0]; Q=[2 0; 0 0]; R=2000;
// Matrices de la planta
A=[0 1;0 0]; B=[0;1];
//Ec de Riccati
k0=[H(1,1);H(1,2);H(2,2)]; //K(tf)=H, condicion final de ganancias K
s0=-H*[d+0.1*tf^3;0]; //S(tf)=-H*r(tf), condicion final de ganancias S
K0=[k0;s0]; //ambas Condiciones inciales
a11=A(1,1);a12=A(1,2);a21=A(2,1);a22=A(2,2); b1=B(1);b2=B(2);
q11=Q(1,1);q12=Q(1,2);q21=Q(2,1);q22=Q(2,2);
function dk=ecriccati(ts,k)//Ecs. de estado
r1=d+0.1*(tf-ts+t0)^3; //SENAL DE REFERENCIA DESEADA PARA x_1
r2=0; //SENAL DE REFERENCIA DESEADA PARA x_2 con t=tf-ts+t0
k11=k(1);k12=k(2);k22=k(3); s1=k(4); s2=k(5);
dk(1) = q11 + 2*a11*k11 + 2*a21*k12 - (b1*k11*(b1*k11 + b2*k12))/R - (b2*k12*(b1*k11
+ b2*k12))/R
dk(2) = q12 + a11*k12 + a12*k11 + a22*k12 + a21*k22 - (b1*k12*(b1*k11 + b2*k12))/R -
(b2*k22*(b1*k11 + b2*k12))/R
dk(3) = q22 + 2*a12*k12 + 2*a22*k22 - (b1*k12*(b1*k12 + b2*k22))/R - (b2*k22*(b1*k12
+ b2*k22))/R
dk(4) = s1*(a11 - (b1*(b1*k11 + b2*k12))/R) - q12*r2 - q11*r1 + s2*(a21 -
(b2*(b1*k11 + b2*k12))/R)
dk(5) = s1*(a12 - (b1*(b1*k12 + b2*k22))/R) - q22*r2 - q21*r1 + s2*(a22 -
(b2*(b1*k12 + b2*k22))/R)
endfunction
//Solucion de la ec Riccati
k=(ode('rk',K0,t0,ts,ecriccati))';//solucion de la ec de estado
t=tf-ts+t0; //tiempo real

//Modelo de estado
function dx=planta(t,x)//Ecs. de estado
x1=x(1);x2=x(2);

r1=d+0.1*(t)^3; //SENAL DE REFERENCIA DESEADA PARA x_1
r2=0; //SENAL DE REFERENCIA DESEADA PARA x_2 con t=tf-ts+t0
e1=x1-r1;e2=x2-r2;//errores de seguimiento
u=-inv(R)*B.'*[k11,k12;k12,k22]*[x1;x2]-inv(R)*B.'*[s1;s2];
dx=[A*[x1;x2]+B*u;[e1;e2]'*Q*[e1;e2]+u'*R*u];
endfunction
//---Simulacion de modelo y generacion de senales --------

NT=size(t,1);//cantidad de muestras de tiempo

xA=d*ones(NT,1)+0.1*t.^3; //Coordenada X de aeronave
yA=zeros(NT,1); //Coordenada Y de aeronave
zA=200*ones(NT,1); //Coordenada Z de aeronave
xM=x(:,1); //Coordenada X de misil
yM=yA; //Coordenada Y de misil
zM=zA.*(1-exp(-xM/(0.2*0.7*max(xM)))); //Coordenada Z de misil
r1f=d+0.1*(tf)^3; //SENAL DE REFERENCIA DESEADA PARA x_1 en tf
r2f=0; //SENAL DE REFERENCIA DESEADA PARA x_2 en tf
89
e1f=x($,1)-r1f;e2f=x($,2)-r2f;//errores de seguimiento
J=x(:,3);J($)=J($)+[e1f;e2f]'*H*[e1f;e2f];//Valores del indice de costo
for j=1:NT//muestras de tiempo
u(j)=-inv(R)*B(2)*k(j,2).*x(j,1)-inv(R)*B(2)*k(j,3).*x(j,2)-
inv(R)*B(2)*k(j,5);//reconstruccion de la entrada u
end
function plot_suelo(ebox)//funcion para el plano del suelo en z=0

x=[ebox(1);ebox(2)]
y=[ebox(3);ebox(4)]
z=zeros(2,2)
surf(x,y,z)
endfunction
//-----GRAFICAS DEL DESEMPEÃ‘O EN LAZO CERRADO------

figure(1);clf;//limpia figura
plot(t,k(:,1),t,k(:,2),t,k(:,3),t,k(:,4),t,k(:,5));//grafica de las ganancias
xgrid//malla
p = get("hdl");p.children.thickness = 2; //Grosor de linea
xtitle('Solucion a la Ec. de Riccati','Tiempo t [s]','k(t)')//titulos
legend('$K_{11}$','$K_{12}$','$K_{22}$','$S_{1}$','$S_{2}$',opt=3) //etiquetas
eboxC=[min(xM),max(xA),-10,10,0,1.1*max(zA)];//Ã¡rea de grafica en cartesianas

figure(2)
for j=1:NT//cantidad de muestras de tiempo
drawlater()//necesario para animacion
clf;//limpia figura
subplot(2,2,1);
plot(t(1:j),xM(1:j),t(1:j),x(1:j,2),t(1:j),xA(1:j),'--');//grafica de estados
xgrid//malla
P=gca();P.data_bounds=[0,tf,min(xM),max(xA)];p = get("hdl");p.children.thickness = 2;
//Grosor de linea
xtitle('Estados del sistema vs tiempo','Tiempo t [s]','x(t)')//titulos
legend('Desplazamiento en X de misil','Velocidad en X de misil','Desplazamiento en X
de aeronave',opt=2) //etiquetas
subplot(2,2,2);
plot(t(1:j),u(1:j));//grafica del control
xgrid//malla
P=gca();P.data_bounds=[0,tf,min(u),max(u)];p = get("hdl");p.children.thickness = 2;
//Grosor de linea
xtitle('Entrada del sistema','Tiempo t [s]','u(t) [m/s^2]') //titulos
subplot(2,2,3);
plot_suelo(eboxC) //Plano del suelo
param3d1(xM(j),yM(j),zM(j)) //misil
curve = gce();curve.mark_mode = "on";curve.mark_style = 12;curve.mark_foreground =
color("blue")
param3d1(xM(1:j),yM(1:j),zM(1:j)) //Trayectoria del misil
param3d1(xA(j),yA(j),zA(j)) //Aeronave
curve = gce();curve.mark_mode = "on";curve.mark_style = 12;curve.mark_foreground =
color("red")
param3d1(xA(1:j),yA(1:j),zA(1:j)) //Trayectoria del aeronave
P=gca();P.data_bounds=eboxC;P.rotation_angles=[80 -110]; //cotas y rotacion de vista
xtitle('Trayectorias en espacio cartesiano, en t='+msprintf('%3.1f',t(j))+' s','eje X
[m]','eje Y [m]','eje Z [m]') //tÃ-tulo
legend('Misil','Trayectoria del misil','Aeronave','Trayectoria de aeronave',opt=4)
//etiquetas
subplot(2,2,4);
plot(t(1:j),J(1:j)); //grafica de funcion de costo
xgrid//dibuja malla
P=gca();P.data_bounds=[0,tf,min(J),max(J)];p = get("hdl");p.children.thickness = 2;
//Grosor de linea
xtitle('Indice de costo','tiempo t [s]','J(t)')//titulos
drawnow()//necesario para animacion
end
90
La solución de la ecuación de Riccati se ilustra en la figura 6.31. La respuesta de los estados en el

tiempo, la señal de control, las trayectorias en el espacio cartesiano y acumulación de costo
están en la figura 6.32.
Figura 6.32. La respuesta de los estados x1 y x2, señal de control óptima u(t), trayectorias en
espacio Cartesiano y acumulación de índice de costo.
Un video que ilustra la representación de modelos en espacio de estados está en:

https://youtu.be/hpeKrMG-WP0=
91
ACTIVIDAD FUNDAMENTAL 4.– Regulación y seguimiento óptimos
Suba a NEXUS un ÚNICO archivo PDF en que muestre:

a) El procedimiento para obtener modelos en espacio de estado y condiciones de las ganancias
óptimas mediante las ecuaciones de Riccati,
b) Simulaciones numéricas y gráficas de las ganancias K(t) y S(t),
c) Simulaciones numéricas y gráficas de la ley de control óptima u* y las variables de estado x
del sistema,
d) y simulaciones numéricas y gráficas del valor del funcional de costo J, de los siguientes
problemas de seguimiento de referencias:
1. El sistema es
 x1  x2

 x2  x1  x2  u
y el índice de desempeño es dado por
 x  sin(t )  0.25u dt

tf
J ( x, u , t )  
2 2
1 .
0
2. El sistema corresponde al circuito:
dv c di
Las variable de estado son x1  vc1 , x 2  i L , y x 3  v c2 (Nota: ic  c
y v L  L L ).
dt dt
Proponga un funcional de costo para que el voltaje vc1 siga a una señal de referencia
r (t )  1  cos(2 t ) . Proponga valores a las resistencias, capacitancias e inductancia.
92
3. Un sistema de nivel de tanques:
Cuyas ecuaciones en el caso linealizado son

h1  h2
 q1
R1
dh1
C1  q  q1
dt
h2
 q2
R2
dh2
C2  q1  q2
dt
Las variables de estado son x1=h1 y x2=h2, u=q es la entrada y y=q2 es la salida.
Proponga valores de los parámetros de la planta R₁ [s/m2], R₂ [s/m2], C₁ [m2], C₂ [m2].
Un funcional de costo, para que el primer tanque mantenga un nivel constante de 3 m y el
 10[h  3]  10[ h2  2]2  u 2 dt .

tf
segundo tanque un nivel constante de 2 m, es J  1
2
t0
93
PRODUCTO INTEGRADOR DE APRENDIZAJE: CASO DE ESTUDIO
Analice un caso de estudio donde requiere optimizar el control de un sistema (NO pueden ser los
mismos ejemplos de los apuntes, particularmente NO los del capítulo 6), proponiendo utilizar los
métodos de control (regulación o seguimiento) y programación del controlador vistas en el curso.
Plasme el planteamiento, soluciones propuestas y resultados de este caso de estudio en un reporte que
debe subir a NEXUS, como un archivo de MS WORD.
Este archivo debe tener portada con nombre de la unidad de aprendizaje, leyenda que diga ‘Examen
ordinario’, nombre y número de matrícula de estudiante, hora de clase, nombre de profesor, fecha de
entrega. Debe estar escrito en fuente 'Palatino Linotype' de 11 puntos a espacio sencillo. Las
expresiones matemáticas deben estar escritas claramente usando MathType, Microsoft Equation o el
Editor de Ecuaciones (NO insertadas como imágenes).
Dicho reporte debe cumplir los siguientes puntos:

1. Planteamiento del problema, sistema o proceso a controlar (10pts)
2. Procedimiento para obtención de un modelo en espacio de estado (10pts)
3. Análisis de propiedades del sistema y revisión de especificaciones del objetivo de control (10pts)
4. Propuesta de un funcional de costo, en base a las especificaciones (10pts)
5. Desarrollo de cálculos de la técnica seleccionada para control óptimo (10pts)
6. Primera aproximación basada en resultados numéricos y simulaciones, y comparación de las
capacidades vs especificaciones del sistema (10pts)
7. Ajustes a la propuesta y nuevos resultados numéricos y simulaciones (10pts)
8. Implementación final en simulaciones o prototipos con la solución definitiva, mostrando gráficas
de la respuesta de los estados del sistema y del funcional de costo (30pts)
9. Conclusiones y recomendaciones (obligatoria)
10. Originalidad >70% (obligatoria)
Además, se debe adjuntar un archivo ZIP debe contenga todos los demás archivos de código de
programación utilizados, así como un video corto que ilustre el funcionamiento en simulaciones o
experimentos reportados.
*Se acepta SÓLO por NEXUS sin excepción, por lo que debe subirlo antes de la fecha y hora de
vencimiento para evitar contratiempos o fallas de la plataforma. La aclaración de DUDAS será SÓLO
durante las clases síncronas.
Finalmente, los alumnos asisten a exponer su proyecto en la fecha de ordinario, y se revisa que las
simulaciones funcionen correctamente y cumplan con las especificaciones propuestas.
*Como ayuda del inciso 8, algunos videos que ilustran el enlace de software de simulación de circuitos
(como Proteus) con tarjetas de control (como Arduino) son:
https://youtu.be/z3okYUvfqmg
https://youtu.be/IDLiubs_rc4
94
APÉNDICE
Lectura de datos del puerto COM virtual a MS Excel
Para adquisición de datos de Arduino virtual desde Proteus (en su versión 8.11) y
leerlos en MS Excel, es necesario descargar e instalar librerías que no se encuentran
incluidas en el software. Dichas librerías pueden encontrarse en
www.TheEngineeringProjects.com
Para la transmisión de datos de Proteus a Excel con el “Data Streamer” se instala dicho
complemento en la siguiente forma:
 Es necesario usar Excel para el vertido de datos requiere activar el complemento

Excel Data Streamer ya que esta herramienta no requiere más configuraciones de
adicionales para su instalación.
 Para activar el plug in en las opciones de Excel seleccionar Complementos COM
(Archivo > Opciones > Complementos > Administrar: Complementos COM > Ir)
y en la ventana emergente seleccionar la casilla Microsoft Data Streamer for
Excel y aceptar.
Fig. A.1. Activación de plug in para opciones avanzadas
 Windows 10 no detecta los puertos seriales virtuales por defecto por lo que se
tuvo que conseguir un programa capaz de hacer que el sistema lo reconozca, en
este caso el programa resulto ser un controlador para puertos virtuales. El
programa seleccionado fue el Virtual Serial Port Driver 6.9.
100
Fig. A.2. Interfaz del Virtual Serial Port Driver
 Este programa crea dos puertos virtuales cuya utilidad se ve reflejada en la

implementación realizada, pues permite crear los puertos que sean necesarios
para adquirir información.
Fig. A.3. Visualización de los puertos virtuales creados
 Es importante verificar que la instalación haya concluido de forma correcta ya

que de lo contrario no se puede hacer comunicación con Proteus. La forma de
verificarlo es que aparezcan los puertos creados en el administrador de
dispositivos mencionado anteriormente.
101
 De la parte del Proteus con el COMPIM hay que configurar el Physical port COM
1 como puerto de salida de datos y tanto en Physical como en Virtual Baud Rate
9600 para que pueda comunicarse de forma adecuada con el Excel.
Fig. A.4. Configuración del COMPIM
 Ya corriendo la simulación de Proteus se puede constatar que en el puerto

configurado hay salida de información por parte del puerto en Sent del Virtual
Serial Port Driver.
102
Fig. A.5. Visualización de datos enviados por medio de puerto virtual
 Regresando a la pestaña Transmisor de datos en Excel, se selecciona el puerto

COM2 debido a que el puerto COM1 ya está siendo utilizado por Proteus.
Fig. A.6. Visualización de Puertos virtuales en Data Stream o “Transmisor de Datos” de

Excel
 Antes de iniciar a simular en Proteus, para poder recibir información se debe ir a

la pestaña hoja de cálculo creada llamada para modificar los valores por defecto
para poder recibir la cantidad de datos que se requiera.
103
Fig. A.7. Configuración del libro
 Con los pasos seguidos anteriormente, se puede realizar la transmisión de datos

de Proteus a Excel, por medio del programa, cabe añadir que en la parte superior,
se tiene que dato clic la opción de Iniciar Datos, para que empiece a recibir
información del programa de Proteus mientras ya está simulando
104
Fig. A.8. Datos transmitidos de Proteus a Excel
 Por último, así es como debe de visualizarse la transmisión de ambos puertos a

través de la interfaz del Virtual Serial Port Driver, uno enviando datos y el otro
recibiendo datos.
Fig. A.9. Interfaz de Datos Transmitidos
105

Apuntes Control Optimo Cap 5-6

Cargado por

Copyright:

Formatos disponibles

Apuntes Control Optimo Cap 5-6

Cargado por

Información del documento

Descripción original:

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Apuntes Control Optimo Cap 5-6

Cargado por

Copyright:

Formatos disponibles

Apuntes de Control Óptimo Dr.

Juan Ángel Rodríguez Liñán

CAPÍTULO 5 - Control óptimo de seguimiento en sistemas

Ejemplo 1: Las ecuaciones diferenciales

La solución numérica de las ecuaciones de Riccati se ilustra en la Figura 2.

Figura 5.2. Solución numérica de las ecuaciones de Riccati.

5.2 Descripción del problema de seguimiento

5.3 Resultados de tiempo finito

Supóngase que el tiempo final es fijo. El Hamiltoniano es

Entonces, se pueden rescribir como un conjunto de ecuaciones diferenciales lineales variantes

La solución de estas ecuaciones es de la forma

Sustituyéndola en (0.26), eliminando x*(tf), y resolviendo para p*(t), obtenemos

Puesto que (0.27) y (0.29) son

Solución: Para este problema

Partiendo de la ec. de Riccati:

Se obtienen las ecs. diferenciales:

Para la ecuación de la referencia:

La ley de control óptimo es:

𝑢(𝑡) = −16.28𝑥 − 4.78𝑥 + 42.36

Esto se simula computacionalmente de la siguiente forma:

k=k($:-1:1,:);// Reordena ganancias

La solución numérica de las ecuaciones de Riccati se ilustra en la Figura 5.1.

Figura 5.1. Solución numérica de las ecuaciones de Riccati.

Figura 5.2. La respuesta de los estados x1 y x2.

Figura 5.3. Señal de control óptima u(t).

𝐽 = 0.5[𝑥 − 0.2𝑡 ] + ([𝑥 − 0.2𝑡] + 0.25𝑢 ) 𝑑𝑡

Solución: Para este problema

Partiendo de la ec. de Riccati, se obtienen las ecs. diferenciales:

La ley de control óptimo es:

Esto se simula computacionalmente de la siguiente forma:

k=k($:-1:1,:);// Reordena ganancias

k12=interp1(ts,k(:,2),t,'linear'); //Interpolate the data (ts,k) at time t

La solución numérica de las ecuaciones de Riccati se ilustra en la Figura 5.4.

Figura 5.4. Solución numérica de las ecuaciones de Riccati.

Figura 5.5. La respuesta de los estados x1 y x2.

Figura 5.6. Señal de control óptima u(t).

5.4 Control óptimo para seguimiento de casos particulares

Considere la ecuación de estado dada por

Puesto que el seguimiento de x(t ) a r (t ) es equivalente a que e(t )  0 , un control óptimo de

El control óptimo es de la forma

Considere un modelo de la forma

El control óptimo es de la forma

Ya que este ejemplo corresponde al caso 2, entonces el control óptimo es de la forma

//Matrices del sistema

//Matrices del funcional de costo

//Cálculo de ganancias del LQR

Figura 5.7. La respuesta de los estados x1 y x2.

Figura 5.8. Señal de control óptima u(t).

Considere un modelo de la forma

Proponga un criterio de costo J y calcule el control óptimo correspondiente.

5 x (t )  r (t )  0.5 x (t)  r (t)  0.5u (t )dt

Ya que este ejemplo corresponde al caso 2, entonces el control óptimo es de la forma

//Matrices del sistema

//Matrices del funcional de costo

//Cálculo de ganancias del LQR

La respuesta de los estados x1 y x2 se ilustra en la Figura 5.9 y la señal de control en la Figura

Figura 5.9. La respuesta de los estados x1 y x2.

Figura 5.10. Señal de control óptima u(t).

CAPÍTULO 6 - Ejemplos de aplicación

Sustituyéndola en (0.26), eliminando x(tf), y resolviendo para p(t), obtenemos

u= k1V1 +k2V2 +k0; //Ecuación control óptimo de 0 a 5 volts

u1= -k1h1 -k2h2 -k3*h3 +k7; //Ecuación control óptimo en Volts