A. Rojas M - 2
A. Rojas M - 2
MULTIVARIABLE
APLICACIONES EN TIEMPO REAL
Modelado de Procesos
Implementación en Tiempo Real
Control Óptimo Cuadrático
Control Adaptativo con Modelo Referencial
Linealización de la Realimentación
Control con Modos Deslizantes
Programas Fuente en MATLAB y LabVIEW
II
Copyright
c 2007 Arturo Rojas-Moreno. Todos los derechos reservados.
ISBN
Queda rigurosamente prohibida la reproducción total o parcial de esta obra por cualquier
medio o procedimiento, sin la autorización escrita del propietario del “Copyright”.
A la Memoria de mis Padres
Índice general
III
Prefacio VII
2. Control Óptimo 67
2.1. Configuración del Sistema de Control Óptimo . . . . . . . . . . . . . . 67
2.2. El Sistema Dinámico No Lineal Multivariable . . . . . . . . . . . . . . 68
2.3. El Controlador PI de Realimentación de Estados . . . . . . . . . . . . 69
2.4. El Observador No Lineal Multivariable . . . . . . . . . . . . . . . . . . 72
2.5. Procedimiento de Diseño . . . . . . . . . . . . . . . . . . . . . . . . . . 73
. Bibliografı́a 139
Modelado de Sistemas No
Lineales
posición incorporado, el cual se emplea para medir la posición angular del brazo del
manipulador en cada instante de tiempo. El servomotor posee una caja de engranajes
para reducir la velocidad en su eje de salida; de de esta manera, se facilita el control
de posición del manipulador.
El subsistema mecánico consiste de un brazo accionado gracias al torque rota-
cional generado en el actuador (el servomotor D.C.). En el extremo del brazo robótico
se puede acoplar un efector final, el cual puede ser una pinza para asir objetos, una
herramienta para soldar, una herramienta para pintar, etc. En nuestro caso usamos
una pinza con dos grados de libertad: un grado para rotar la pinza y otro para abrirla
y cerrarla. Con propósitos de modelado, vamos a suponer que el efector final y su
carga se pueden modelar mediante una masa mh variable. La tabla 1.1 describe las
variables y los valores de los parámetros del manipulador mostrado en la figura 1.1.
El manipulador de 1GDL es del tipo SISO (“Single Input Single Output”)
ya que sólo posee una entrada: el voltaje de control u aplicado a la armadura del
servomotor, y una salida: la posición angular θ del brazo.
La ecuación (1.2) es cierta debido a que el espacio angular recorrido por el engranaje
de menor radio debe de ser n veces mayor que el espacio recorrido por el el engranaje
de radio mayor. Por otra parte, el principio de la conservación de la energı́a establece
que el trabajo realizado por el engranaje de la izquierda debe ser igual al trabajo
realizado por el engranaje de la derecha, es decir:
Tg2 = Jg θ̈ + Bg θ̇ + TL (1.4)
TL = JL θ̈ + BL θ̇ + τL (1.5)
L
τL = mb g senθ + mh g (L + rh )senθ = Q sen θ (1.6)
2
L
Q = mb g + mh g (L + rh )
2
donde JL y BL representan el momento de inercia y la constante de fricción viscosa de
la carga no lineal (brazo más efector final), g es la constante gravitacional, mb y mh
denotan las masas del brazo y del efector final (esta masa también incluye la masa de
la carga en el efector) respectivamente, y rh denota la distancia desde el extremo del
brazo al centro de masa de mh . Para fines prácticos asumiremos que rh ≈ 0. Notar
en (1.5) que el torque τL se debe a las fuerzas ejercidas por los pesos del brazo y de
la esfera. Ası́, el torque mb g L senθ
2 se debe al peso mb g del brazo, mientras que el
torque mh g (L + rh )senθ ≈ mh g Lsenθ es causado por el peso mh g del efector, Las
distancias L senθ
2 y (L + rh )senθ ≈ Lsenθ son los correspondientes brazos de palanca
del brazo y del efector respectivamente.
El momento de inercia JL de la carga es la suma del momento de inercia del
brazo Jb más el momento de inercia del efector Jh . Asumiendo que la masa mh
está concentrada en su C.M. (centro de masa), entonces:
Jh = mh (L + rh )2 ≈ mh L2 (1.7)
Del mismo modo, asumiendo que la masa mb del brazo se concentra en su C.M.,
entonces: 2
L
Jb = mb (1.8)
2
Sin embargo, si se considera que la masa del brazo está distribuida a lo largo de su
longitud, entonces debemos aplicar el teorema de los ejes paralelos, el cual establece
que el momento de inercia de una masa m alrededor de un eje que no pasa por su
centro de masa está dado por:
J = Jo + ma2 (1.9)
1.1 Modelado Empleando las Leyes de la Fı́sica 5
donde a es la distancia entre el eje que pasa por el centro de masa de m y el eje
paralelo, y Jo es el momento de inercia alrededor del eje que pasa por su centro de
masa. Aplicando el teorema de los ejes paralelos al brazo, el momento de inercia Jb
con respecto al punto de articulación se formula como:
2
1 L 1
Jb = mb L2 + mb = mb L2 (1.10)
12 2 3
Por simplicidad, sin embargo, emplearemos las fórmulas dadas en (1.7) y (1.8).
donde:
Jeq = n2 Jm + Jg + JL Beq = n2 Bm + Bg + BL
Las expresiones de Q, Jh y Jb (tener en cuenta que JL = Jh + Jb ) se dan en (1.6),
(1.7) y (1.8) respectivamente.
dia
Va = ia Ra + La + Vb (1.12)
dt
donde ia , Ra y La son la corriente, la resistencia y la inductancia en la armadu-
ra del servomotor respectivamente, y Vb es el voltaje de fuerza contraelectromotriz
gobernado por la relación:
Vb = Kb ωm = Kb nω = Kb nθ̇ (1.13)
6 Modelado de Sistemas No Lineales
Va = uKA (1.14)
Tm = Km ia (1.15)
dia KA Kb n Ra
= u− ω− ia (1.16)
dt La La La
dω Q Beq nKm
=− senθ − ω+ ia (1.17)
dt Jeq Jeq Jeq
Ecuación de Estado
Las ecuaciones (1.16) y (1.17) describen el modelo no lineal del proceso. Eligien-
do como variables de estado: x1 = θ (posición angular), x2 = ω (velocidad angular)
y x3 = ia (corriente de armadura), se obtiene:
ẋ1 = x2
Q Beq nKm
ẋ2 = − senx1 − x2 + x3
Jeq Jeq Jeq
nKb Ra KA
ẋ3 = − x2 − x3 + u (1.18)
La La La
donde la salida es la posición x1 y la señal de control es u.
En la Tabla 1.1 podemos observar que la inductancia de armadura La del ser-
vomotor es bastante pequeña, de modo tal que puede despreciarse sin que se pierda
exactitud en los resultados. Considerando el producto ẋ3 La = 0 en la tercera ecuación
de (1.18) y despejando la corriente de armadura x3 obtenemos:
KA nKb
x3 = u− x2 (1.19)
Ra Ra
ẋ1 = x2
P Beq Ra + n2 Km Kb nKm KA
ẋ2 = − senx1 − x2 + u (1.20)
Jeq Jeq Ra Jeq Ra
1.1 Modelado Empleando las Leyes de la Fı́sica 7
ẋ = Ax + Bu (1.21)
ẋ1 0 1 x1 0
= Beq Ra +n2 Km Kb + nKm KA u
ẋ2 − JQeq − Jeq Ra
x2 Jeq Ra
Ejemplo 1.1
Graficar la respuesta (la posición angular θ) del modelo no lineal de segundo orden y
de los modelos lineales continuo y discreto derivados de dicho modelo. Para fines de
análisis asumir que la señal de entrada u es de 1.4 voltios de magnitud, que BL = Bg
y que el efector final es una esfera de masa mh y de radio rh . En esta situación:
2
Jh = mh rh2 + mh (L + rh )2
5
Solución: El programa C3rpta.m determina tales respuestas (ver figura 1.3) ası́ como
también las funciones de transferencia, la estabilidad, la controlabilidad y la obser-
vabilidad de los modelos lineales continuo y discreto. ♣
0.8
0.4
0.2
0
0 50 100 150
6
y(t) LINEAL [rad]
0
0 50 100 150
6
y(k) LINEAL [rad]
0
0 50 100 150
Tiempo en segundos
Figura 1.3: Respuestas de los modelos no lineal, lineal continuo y lineal discreto del
brazo robótico a una entrada u escalón de magnitud 1.4 voltios.
M2
Jm mH
Ra JL
+ + + Tm Tg1 L
θm BL
u Va Vb n Bg K ω m
- - -
KA La θ m Bm
Gearbox
Jg
n Tg2 θ
Ia
Jeq = n2 Jm + Jg Beq = n2 Bm + Bg
La ecuación que gobierna la dinámica del brazo del manipulador se puede expresar
como:
θm 1
Kw − θ = JL θ̈ + BL θ̇ + mb gLsen θ + mh gLsin θ (1.24)
n 2
donde JL y BL representan el momento de inercia y la constante de fricción viscosa
de la carga no lineal (brazo más efector), g es la constante gravitacional, mb y mh
(esta masa también incluye la masa de la carga) denotan la masa del del brazo y
del efector respectivamente y, 12 mb gLsen θ y mh gLsin θ son los torques debido a los
pesos del brazo y del efector respectivamente. El momento de inercia de JL = Jh +Jb ,
asumiendo que las masas mh y mb se concentran en sus respectivos C.M. se da en las
ecuaciones (1.7) y (1.8):
2
2 L
JL = mh L + mb
2
Para completar el modelado de la parte eléctrica del proceso MRAE, podemos
aseverar que:
dia
Ra ia + La + Vb = KA u (1.25)
dt
donde KA es la ganancia del amplificador y Vb es el voltaje de la fuerza contra
electromotrı́z y responde a la relación:
Vb = Kb θ̇m (1.26)
Tm = Km ia (1.27)
donde:
Jeq Beq
Tm 0 θ̈m
0 θ̇m Kw θm
−θ
n n n n n n
= + +
0 0 JL θ̈ 0 BL θ̇ d2
1 θm
d2 = mb + mh Lg sin θ − Kw −θ
2 n
Despreciando la inductancia de armadura La en (1.25), y despejando ia se obtiene:
Kb KA
ia = − θ̇m + u (1.29)
Ra Ra
Sustituyendo ia en Tm = Km ia de (1.28) y despejando u, el modelo de Lagrange
toma una nueva forma:
θ̈ θ̇
u m11 0 m p 11 0 m d1
= n n
+ + (1.30)
0 0 JL θ̈ 0 BL θ̇ d2
Ra Jeq nKb Ra Beq Ra Kw θm
m11 = p11 = + d1 = −θ
nKA Km KA nKA Km nKA Km n
f1 (x, u) = x2
Kw BL Kw Lg mb
f2 (x, u) = − x1 − x2 + x3 − + mh sen x1
JL JL JL JL 2
f3 (x, u) = x4
Kw Kw Beq nKm
f4 (x, u) = x1 − x3 − x4 + x5
Jeq Jeq Jeq Jeq
nKb Ra KA
f5 (x, u) = − x4 − x5 + u
La La La
donde hemos usado el hecho de que ẋ1 = x2 y ẋ3 = x4 . Si la salida del sistema es
y = θ, entonces la ecuación de salida del MRAE resulta:
y = h(x) = Cx = [1 0 0 0 0] x (1.32)
1.1 Modelado Empleando las Leyes de la Fı́sica 11
f1 (x) = x2
Kw BL Kw Lg m
f2 (x) = − x1 − x2 + x3 − + mH sin x1
JL JL JL JL 2
f3 (x) = x4
2
Kw Kw n Km Kb Beq nKm KA
f4 (x) = x1 − x3 − + x4 + u
Jeq Jeq Jeq Ra Jeq Jeq Ra
La salida del sistema en este caso se expresa como:
y = h(x) = Cx = [1 0 0 0] x (1.34)
+ Fuerza de
u control y’
-
Servomotor
D.C. y
θ Pendulo
z F
Carro
11
00
00
11
00
11
l v /2
me g
θ
mvg lv
le
z’
0
z
P
F
lv
ze = z + le senθ zv = z + senθ (1.35)
2
Para simplificar el procedimiento de modelado del péndulo, aplicaremos las leyes
de Newton para los movimientos lineal y rotacional, considerando al péndulo como
un todo. Para el movimiento lineal, dicha ley establece que para un sistema de N
partı́culas:
N
d2 M
mi 2 ri = Fj (1.36)
dt
i=1 j=1
d2 d2 d2
mc z + m z
e 2 e + mv 2 zv = F (1.37)
dt2 dt dt
Sustituyendo ze y zv (ecuación (1.35) en (1.37)) resulta:
d2 d2 d2 lv
mc 2
z + me 2
(z + l e sen θ) + mv 2
(z + sen θ) = F (1.38)
dt dt dt 2
14 Modelado de Sistemas No Lineales
Vb = Kb θ̇m (1.44)
El torque Tg2 a la salida del tren de engranajes debe vencer a la inercia Jg de los
engranajes, a la inercia Jp de la polea y al torque de polea F rp (rp es el radio de la
polea):
Tg2 = Jg θ̈ + Bg θ̇ + Jp θ̈ + Bp θ̇ + F rp = nTg1 (1.48)
Sustituyendo (1.48) en (1.46) se obtiene:
donde:
Jeq = Jm + n2 (Jg + Jp ) Beq = Bm + n2 (Bg + Bp )
Para transformar el desplazamiento angular θ del servomotor en el desplazamiento
horizontal z del carro, empleamos:
z = rp θ (1.50)
Ra θm
+
u Va Vb Jo rp
-
La Jm F
Bm Bo
Las ecuaciones (1.42) y (1.53) representan el modelo dinámico no lineal del sistema
péndulo invertido controlado por la corriente de armadura. Tales ecuaciones pueden
ser escritas en forma compacta como:
donde:
lv
M1 = mc + me + mv ; M2 = me le + mv ; J1 = Je + Jv
2
y puesto que en nuestro sistema las salidas disponibles son el desplazamiento angular
x1 de la varilla y el desplazamiento x3 del carro, la ecuación de salida toma la forma:
y = h(x) = Cx (1.56)
donde:
1 0 0 0
C=
0 0 1 0
1.1 Modelado Empleando las Leyes de la Fı́sica 17
Para la operación del proceso alrededor del estado de equilibrio xo = [0, 0, 0, 0]T = 0
y uo = 0, las matrices A, B, C y D pueden ser determinadas evaluando las siguientes
matrices jacobianas:
∂f1 ∂f1 ∂f1 ∂f1
∂X1 · · · ∂Xn ∂U1 · · · ∂Um
.. ..
A = ... ..
. . B = ... ..
. .
∂fn ∂fn ∂fn ∂fn
∂X1 ··· ∂Xn (xo , uo ) ∂U1 ··· ∂Um (xo , uo )
∂h1
∂X1 ··· ∂h1
∂Xn
∂h1
∂U1 ··· ∂h1
∂Um
.. .. .. .. .. ..
C = . . . D= . . . (1.58)
∂hr
∂X1 ··· ∂hr
∂Xn (xo , uo )
∂hr
∂U1 ··· ∂hr
∂Um (xo , uo )
ẋ = Ax + Bu (1.59)
donde:
0 1 0 0 0
(M1 +J2 )M2 g Bx M2 −Kx M2 KA
(M1 +J2 )J1 −M22
0 0 (M1 +J2 )J1 −M22 (M1 +J2 )J1 −M22
A=
B=
0 0 0 1 0
−M22 g −J1 Bx J1 Kx KA
(M1 +J2 )J1 −M22
0 0 (M1 +J2 )J1 −M22 (M1 +J2 )J1 −M22
Sistema Grúa–Puente
El modelado del sistema grúa-puente es similar al modelado del péndulo inver-
tido. En este caso caso el péndulo debe apuntar hacia abajo, tal como se muestra en
la figura 1.8. Al igual que en el caso del péndulo invertido, para mayor facilidad, el
sistema grúa puente se puede subdividir en dos subsistemas: carro–varilla y motor–
polea. El subsistema carro-varilla está representado en la figura 1.9. De dicha figura
podemos observar que los centros de gravedad zv de la varilla y ze de la esfera son:
lv
ze = z − le sen θ zv = z − sen θ
2
+ Fuerza de
u control
-
Servomotor
D.C. y
z F
Carro
θ
’ puente
Grua
y’
y y’
z
P z’
0
F z
lv
θ
2
111
000 le
000
111
000
111
000 m v g
111 lv
me g
donde:
lv
M1 = mc + me + mv ; M2 = me le + mv ; J1 = Je + Jv
2
Como era de esperarse, las relaciones y (1.54) (1.60) sólo se diferencian en los signos.
dh
S = S ḣ = qi − qo (1.63)
dt
Considerando un flujo laminar de salida:
h
qo = (1.64)
Rh
H
Rh = (1.65)
Q
20 Modelado de Sistemas No Lineales
θo Temperatura perturbacional: θo = Θo − Θo oC
Determinación Experimental de Rh
La válvula de control empleada para regular la entrada de agua al tanque es del
tipo VXN015F250, con diámetro nominal DN15, conexión G1B y actuador motórico.
La abertura máxima se obtiene alimentando con 10 V al actuador, la cual corresponde
a un flujo de 0.4 m3 /h, de acuerdo al manual del fabricante. La mı́nima abertura,
con 0 V, corresponde a un flujo de 0 m3 /h. Esto significa que para una abertura de
1 V el flujo que pasa por la válvula es de 1/90000 m3 /s.
Asumiendo una variación lineal entre el flujo Qi que pasa por la válvula y
la altura H del tanque, se realizó el siguiente experimento. Con una abertura de
válvula para 4 V (0.16 m3 /h), se abrió convenientemente la válvula de descarga
hasta lograr una altura estable de 0.12 m. Luego, con una abertura de válvula para
6 V (0.24 m3 /h), se siguió abriendo la válvula de descarga hasta lograr una altura
de 0.18 m. Empleando la relación (1.65), la resistencia hidráulica para cada punto
resultó aproximadamente Rt = 2700 s/m2 . Se asume también que los valores estables
de H y Qi son H = 0.12 m y Q = 0.16 m3 /h.
22 Modelado de Sistemas No Lineales
dΘo
La segunda ecuación de estado se obtiene despejando dt = Θ̇o de (1.77):
a Θo 1 Θo Θa 1 Θi Qi 1 Φi
Θ̇o = − √ − + + + (1.82)
S H SρCp Rt H SρCp Rt H A h SρCp H
Para linealizar el proceso tanque emplearemos la técnica del Jacobiano. Para ello se
definen los siguientes vectores de estado y de control:
X1 H Ẋ1 f1 U1 Qi
X= = ; Ẋ = = ; U= = (1.83)
X2 Θo Ẋ2 f2 U2 Φi
∂f2 a 1 1 1
= − −
∂X2 S X 1 SρCp Rt X 1
∂f2 θi 1
=
∂U1 S X1
∂f2 1 1
= (1.86)
∂U2 SρCp X 1
ẋ = Ax + Bu (1.87)
donde:
∂f1 ∂f1 ∂f1 ∂f1
A= ∂X1 ∂X2 B= ∂U1 ∂U2
∂f2 ∂f2 ∂f2 ∂f2
∂X1 ∂X2 ∂U1 ∂U2
24 Modelado de Sistemas No Lineales
De la ecuación (1.88), es claro que: L = L(q1 , . . . , qr , q̇1 , . . . , q̇r ). Por otra parte, de
acuerdo al principio de de la mı́nima acción de Hamilton para sistemas conservativos,
la integral I definida por:
t2
I= L(q1 , . . . , qr , q̇1 , . . . , q̇r )
t1
θ1 L1cos θ 1
M 1 Pulley Gripper
(lateral view)
M2 Link
Cart x Rp
F F1 Pulley
r L1 sinθ 1
Tabla 1.4: Parámetros valorados del sistema MRT. La abreviatura C.M. significa
Centro de Masa.
Sı́mbolo Descripción Valor Unidades
u 1 , u2 Voltaje de entrada al sistema V
KA Ganancia del amplificador 15
Va1 , Va2 , Voltaje de armadura V
Ra Resistencia de armadura 2 Ω
La Inductancia de armadura 0.002 H
ia 1, ia2 Corriente de armadura A
Km Constante del torque motor 31.071×10−3 N-m/A
Tm1 , Tm2 Torque motor N-m
Tg1 Torque de entrada a los engranajes N-m
Tg2 Torque de salida de los engranajes N-m
Jm Momento de inercia del motor 1.9062×10−6 kg-m2
Jg Momento de inercia de los engranajes ∼
=0 kg-m2
Jp Momento de inercia de la polea kg-m2
Jh Momento de inercia del efector kg-m2
JL Momento de inercia en la carga kg-m2
Jb Momento de inercia del brazo kg-m2
Jeq1 , Jeq2 Momento de inercia equivalente 4.2×10−6 kg-m2
Beq1 , Beq2 Constante de fricción equivalente 3.36×10−6 N-m-s/rad
Bm Constante de fricción del motor 1.8338×10−6 N-m/rad/s
Bg Constante de fricción en engranajes N-m/rad/s
BL Constante de fricción en la carga N-m-s/rad
Bp Constante de fricción (punto pivote) 1.92×10−3 kg-m2 /s
Bj Constante de fricción de la articulación 1.92×10−3 kg-m2 /s
F Fuerza aplicada al carro N
Ff Fuerza de rozamiento N
Fc Constante de fricción del carro 2.81 kg/s
mh Masa del efector kg
mb Masa del brazo 0.103 kg
mc masa del carro 0.9574 kg
mp masa de la polea 0.2 kg
L Longitud del brazo 0.225 m
rh Distancia: C.M. del efector al brazo ∼
=0 m
rp Radio de la polea 0.05 m
r Posición del carro m
L Longitud del brazo 0.225 m
xh , xb Posiciones horizontales m
yh , yb Posiciones verticales m
Vb1 , Vb2 Voltaje contraelectromotriz V
Kb Constante contraelectromotriz 25×10−3 V-s/rad
g Aceleración de la gravedad 9.81 m/s2
N 1 , N2 N o de dientes de los engranajes N2 > N 1
n Relación de engranajes (n = N2 /N1 ) 12.5
θm Posición angular del motor rad
θ Posición angular del brazo rad
ω Velocidad angular de la carga rad/s
ωm Velocidad angular del motor rad/s
1.2 Método de Las Ecuaciones de Lagrange 27
L = V − U = (V1 + V2 ) − (U1 + U2 )
Ra La θ m1
+ + 111
000
Tm1 Tg1
J
000 B gg Pulley
111
u1 Vb 1
- KA - 111
000
B
000
111
n 000
111
Ia 1 Jm m F1
Rp
Tg2
θ̈m1 θ̇m1
Tg2 = nTg1 = (Jg + Jp ) + (Bg + Bp ) + F Rp (1.101)
n n
donde n > 1 es la relación de dientes de los engranajes del mecanismo de reducción,
Jm , Jg y Jp son los momentos of inercia de la armadura, del mecanismo de reducción
y de la polea respectivamente, mientras que Bm , Bg y Bp son las constantes de
fricción de la armadura, del mecanismo de reducción y de la polea respectivamente.
1.2 Método de Las Ecuaciones de Lagrange 29
La relación Tg2 = nTg1 se obtiene asumiendo que los engranajes son ideales. En esta
situación, el principio de conservación de energı́a requiere que:
θm1
Tg1 θm1 = Tg2
n
El torque servomotor Tm1 es proporcional a ia1 :
θm1
r= rp (1.103)
n
Usando las ecuaciones (1.98), (1.99), (1.100), (1.101), (1.102) y (1.103) podemos
obtener:
Km KA Jeq1 Beq1 nKm Kb
F = u1 − r̈ − + ṙ (1.104)
rp Ra nrp2 nrp2 nrp2
donde:
Jeq1 = n2 Jm + Jg + Jp Beq1 = n2 Bm + Bg + Bp
donde:
Vb2 = Kb nθ̇ (1.106)
La ecuación del torque motor Tm2 es (ver figura 1.13):
Ra La
nθ1
Tm2 Tg1
Jg
000
111
+ + +
Bg
u Va2 Vb2
-
2
KA - - 111
000
Bm
00 θ 1
11
Ia2
Jm n 11T
00
Tg2 L
Empleando las ecuaciones (1.105), (1.106), (1.107), (1.108) y (1.109) se puede de-
mostrar que:
n2 Km Kb nKm Kb
TL = −Jeq2 θ̈1 − Beq2 + θ̇1 + u2 (1.110)
Ra Ra
donde:
Jeq2 = n2 Jm + Jg Beq2 = n2 Bm + Bg
donde:
m11 m12 cos θ p11 p12 θ̇sin θ 0
M= P= d(q) =
m21 cos θ m22 0 p22 d2 sin θ
Ra r p Jeq1 Ra r p L mb
m11 = mc + mh + mb + 2 2 ; m12 = mh +
Km KA n rp Km KA 2
Ra L mb Ra mb 2
m21 = mh + ; m22 = Jeq2 + mh L2 + L
nKA Km 2 nKA Km 4
Ra r p Beq1 nKm Kb Ra rp L mb
p11 = Fc + 2 2 + ; p 12 = − m h +
KA Km n rp Ra rp2 KA Km 2
Ra n2 Km Kb Ra Lg mb
p22 = Beq2 + Bp + ; d2 = − mh +
nKa Km Ra nKa Km 2
Ejemplo 1.2
1.2 Método de Las Ecuaciones de Lagrange 31
Para el sistema péndulo mostrado en las figuras 1.5 y 1.6 asuma que lv = le = l y
mv = 0. Determinar el modelo en el espacio de estado del sistema para variaciones
pequeñas de la posición θ, empleando primero las ecuaciones de Lagrange y luego las
leyes de la fı́sica.
ze = z + l senθ ye = l cosθ
Ue = me gye
sen θ ∼
=θ cos θ ∼
=1 θ̇2 θ ∼
=0
ẋ = Ax + Bu
0 0 1 0 0
0 0 0 1 0
A=
0 −me g
B=
1
mc 0 0 mc
0 g(mc +me )
M cl 0 0 − m1c l
32 Modelado de Sistemas No Lineales
Aplicando las leyes fı́sicas, para los movimientos horizontal y rotacional, de acuerdo
a las relaciones (1.37) y (1.41) se obtiene:
d2 d2
mc z + me 2 ze = F
dt2 dt
Je θ̈ = me gl sen θ − me z̈l cos θ
Sustituyendo las relaciones ze = z + l senθ, ye = l cosθ en la primera ecuación y
Je = me l2 en la segunda, se obtienen las ecuaciones dadas en (1.114). ♣
1. T(z, θi ): rotación de un ángulo θi alrededor del eje zi−1 para alinear el eje xi−1
con el eje xi .
2. T(z, di ): traslación a lo largo del eje zi−1 de una distancia di para hacer coincidir
los ejes xi−1 y xi .
3. T(x, ai ): traslación a lo largo del eje xi de una distancia ai para hacer coincidir
los orı́genes (xi−1 , yi−1 , zi−1 ) y (xi , yi , zi ) ası́ como también el eje xi .
4. T(x, αi ): rotación de un ángulo αi alrededor del eje xi para hacer coincidir los
dos sistemas de coordenadas.
Ai−1
i = T(z, θi )T(z, di )T(x, ai )T (x, αi )
Cθi −Sθi 0 0 1 0 0 0 1 0 0 ai 1 0 0 0
Sθi Cθi 0 0 0 0 0 −Sαi 0
= 0 1 0 1 0 0 Cαi
0 0 1 0 0 0 1 di 0 0 1 0 0 Sαi Cαi 0
0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1
Cθi −Cαi Sθi Sαi Sθi ai Cθi
Sθi Cαi Cθi −Sαi Cθi ai Sθi
=
0
(1.115)
Sαi Cαi di
0 0 0 1
donde C y S representan las funciones cos y sen respectivamente. La siguiente matriz
de transformación homogénea:
i
A0i = Ajj−1 = A01 A12 · · · Ai−1
i i = 1, 2, . . . , n (1.116)
j=1
Ejemplo 1.3
Solución: Notar que la figura 1.14 ilustra los sistemas de coordenadas D–H del
manipulador esférico. La tabla 1.5 muestra los correspondientes parámetros D–H
extraı́dos de dicha figura. Por consiguiente, las matrices de transformación Ai−1
i (i
= 1, 2) del manipulador esférico toman la forma:
Cθ1 0 Sθ1 0 Cθ2 −Sθ2 0 L2 Cθ2
Sθ1 0 −Cθ1 0 Sθ2 Cθ2 0 L2 Sθ2
A01 =
0
A12 =
1 0 b+h 0 0 1 0
0 0 0 1 0 0 0 1
Cθ1 Cθ2 −Cθ1 Sθ2 Sθ1 L2 Cθ1 Cθ2
Sθ1 Cθ2 −Sθ1 Sθ2 −Cθ1 L2 Sθ1 Cθ2
A02 = A01 A12 =
Sθ2
♣
Cθ2 0 L2 Sθ2 + b + h
0 0 0 1
con A00 = I (la matriz identidad) y donde A0j−1 relaciona el (j −1)-ésimo sistema
de coordenadas con el sistema de coordenadas base, y:
0 −1 0 0
1 0 0 0
Qi =
0 0 0 0
si la articulación i es de rotación
0 0 0 0
0 0 0 0
0 0 0 0
Qi =
0
si la articulación i es prismática
0 0 1
0 0 0 0
36 Modelado de Sistemas No Lineales
donde (x̄i , ȳi , z̄i ) es el centro de gravedad (CG) del eslabón i con respecto al
i-ésimo sistema de coordenadas, Ixi , Iyi , Izi son los momentos de inercia con
respecto al sistema de coordenadas (xi , yi , zi ), mi es la masa del cuerpo i, x̄i ,
y¯i and z¯i son las distancias desde el CG del cuerpo i hacia el sistema de coor-
denadas (xi , yi , zi ), y Ixi yi , Iyi zi , Ixi zi , Ixi xi , Iyi yi , Izi zi son los correspondientes
momentos de inercia producto.
6. Obtener la matriz simétrica inercial de aceleración H de orden n × n, cuyos
elementos son:
Hik = T r Ujk Jj UTji i, k = 1, 2, . . . , n (1.120)
j=máx(i,k)
c = Cq̇ (1.125)
Tabla 1.6: Variables y parámetros del sistema MRE. Las abreviaciones C.G. y c.r.
significan centro de gravedad y con respecto, respectivamente.
Sı́mbolo Descripción Valor Unidad
md Masa del disco 0.4 kg
mp Masa de la barra o prisma 1.0 kg
mb Masa del brazo (eslabón) 0.1 kg
mh Masa del efector final kg
d Espesor del disco 0.01 m
rd Radio del disco 0.07 m
b Longitud de la barra 0.21 m
a Lado de la sección de la barra 0.044 m
Lb Longitud del brazo (eslabón) 0.3 m
Ly1 Longitud del C.G. c.r. a (x1 , y1 , z1 ) m
Lx2 Longitud del C.G. c.r. a (x2 , y2 , z2 ) m
Jh Momento de inercia del efector final kg-m2
Jm Momento de inercia de M1 y M2 1.8×10−4 kg-m2
Bm Constante de fricción de M1 y M2 1.83×10−6 N-m-s/rad
Jg1 , Jg2 Momento de inercia de los engranajes 5.63×10−5 kg-m2
Bg1 , Bg2 Constante de fricción en engranajes 7.05×10−5 N-m-s/rad
n Relación de engranajes de M1 y M2 18.5
Ra Resistencia de armadura de M1 y M2 2.8 Ω
La Inductancia de armadura de M1 y M2 1.5×10−3 H
Vb1 , Vb2 Voltaje contraelectromotriz V
Va1 , Va2 Voltaje de armadura 24 V
ia1 , ia2 Corriente de armadura 0.7 A
KA Ganancia del amplificador 14.9
Km Constante del motor de M1 y M2 0.9 N-m/A
Kb Constante contraelectromotriz 0.9 V-s/rad
u 1 , u2 Voltaje de entrada 1.4 V
g Constante gravitacional 9.81 m/s2
1.3 El Método de Lagrange–Euler 39
y
y x2 2
2 M3 x2
z2 Gripper M3
y
1111
0000
1 Link
L2 z2
q
2
M y
2 Lx2 1
z1 z1 x1 M2
x
1
Ly1 a a
Metallic bar b z
z0 0
of lenght b
: Axis is pointing
Disk away from the plane
h x0 y x y0
0 0
q q M1
1 Servomotor M 1 1
1 h
Iy1d y1d = md rd2 + md ( + b)2
2 2
1 b
Iy1p y1p = mp (a2 + a2 ) + mp ( )2
12 2
donde 12 md rd2 es el M.I. del disco girando en su eje de rotación, md ( h2 + b)2 es el M.I.
1
del disco con respecto al sistema (x1 , y1 , z1 ), 12 mp (a2 + a2 ) es el M.I. de la barra
1.3 El Método de Lagrange–Euler 41
1 Lb
Iz2 z2 = mb L2b + mb ( )2 + Jh = 2I2
12 2
1
donde 12 mb L2b es el M.I. del brazo girando en un eje perpendicular a la barra y que
pasa por su C.M. y mb ( L2b )2 es le M.I. del brazo referido al sistema de coordenadas
(x2 , y2 , z2 ) (ver figura 1.14) y Jh es el M.I. de la mano del manipulador de masa mh ,
la cual incluye también la masa de la carga aplicada a la mano. Luego, la matriz J2
para i = 2 viene a ser:
I2 0 0 −m2 Lx2
0 I2 0 0
J2 =
(1.129)
0 0 −I2 0
−m2 Lx2 0 0 m2
T
h112 = T r(U212 J2 U21 ) = (m2 L22 − 2m2 L2 lc2 − I2 )senq2 cosq2
T
h121 = T r(U221 J2 U21 ) = (m2 L22 − 2m2 L2 lc2 − I2 )senq2 cosq2
T
h122 = T r(U222 J2 U21 )=0
T
h211 = T r(U211 J2 U22 ) = (2m2 L2 lc2 − m2 L22 + I2 )senq2 cosq2
T T T
h212 = T r(U212 J2 U22 ) = h221 = T r(U221 J2 U22 ) = h222 = T r(U222 J2 U22 )=0
c1 = h111 q̇12 + h112 q̇1 q̇2 + h121 q̇2 q̇1 + h122 q̇22
= 2(m2 L22 − 2m2 L2 Lx2 − I2 )q̇1 q̇2 sin q2 cos q2 = −2H22 q̇1 q̇2 sin q2 cos q2
c2 = h211 q̇12 + h212 q̇1 q̇2 + h221 q̇2 q̇1 + h222 q̇22
= (2m2 L2 Lx2 − m2 L22 + I2 )q̇12 sin q2 cosq2 = H22 q̇12 sin q2 cosq2
Ra1 La1
+ + + 111
000
Tm1 Tg1
000 B g1g1
111
J
u Va1 Vb1
-
1
KA1 - - 111
000
B m1
00 θ 1
11
I a1
J m1 n
1
T
11T
00
g2 1
Va1 = KA1 u1 = La1 I˙a1 + Ra1 Ia1 + Vb1 = La1 I˙a1 + Ra1 Ia1 + Kb1 n1 q̇1 (1.133)
mientras que la ecuación del torque Tg2 requerido para accionar el primer eslabón se
expresa como:
Tg2 = nTg1 = Jg1 q̈1 + Bg1 q̇1 + T1 (1.135)
donde T1 es el torque de carga. Substituyendo (1.135) en (1.134) se obtiene:
donde:
Jeq1 = n21 Jm1 + Jg1 Beq1 = n21 Bm1 + Bg1
Reemplazando (1.131) en (1.136) nos conduce a:
donde:
La1 Ra1 n1 Kb1
LT 1 = RT 1 = NT 1 =
n1 KA1 Km1 n1 KA1 Km1 Km1
44 Modelado de Sistemas No Lineales
donde:
Jeq2 = n22 Jm2 + Jg2 Beq2 = n22 Bm2 + Bg2
Ra2 La2
+ + + 111
000
Tm2 Tg1
000
111
J g2
B g2
u Va2 Vb2
-
2
KA1 - - 111
000
B m2
00 θ 2
11
Ia2
J m2 n2
Tg2
11T
00
2
donde:
RT 1 (Jeq1 + 2I1 + H22 cos2 q2 ) 0
M(q) =
0 RT 2 (Jeq2 + H22 )
RT 1 (Beq1 − H22 q̇2 sin q2 cos q2 ) + NT 1 −RT 1 H22 q̇1 sin q2 cos q2
P(q, q̇) =
RT 2 H22 q̇1 sin q2 cos q2 RT 2 Beq2 + NT 2
0
d(q) =
RT 2 m2 g Lx2 cos q2
1.3 El Método de Lagrange–Euler 45
P1 = LT 1 (Ḣ11 q̈1 + Beq1 q̈1 + Ċ1 ) + RT 1 (Jeq1 q̈1 + Beq1 q̇1 + H11 q̈ + C1 ) + NT 1 q̇1
P2 = LT 2 (Beq2 q̈2 + Ċ2 + d˙2 ) + RT 2 (Jeq2 q̈2 + Beq2 q̇2 + H22 q̈ + C2 + d2 ) + NT 2 q̇2
La ecuación de estado para las salidas y1 = x1 and y2 = x2 resulta:
y1 h1 (x) 1 0 0 0 0 0
y= = = x (1.147)
y2 h2 (x) 0 1 0 0 0 0
46 Modelado de Sistemas No Lineales
1.4. PROBLEMAS
Problema 1.1
Problema 1.3
Determinar la descripción no lineal en el espacio de estado del manipulador robótico
traslacional descrito en la subsección 1.2.2, para los casos siguientes:
(a) Despreciando la inductancia de armadura La de los servomotores.
(b) Sin despreciar la inductancia de armadura La de los servomotores.
Problema 1.4
Para el manipulador robótico esférico descrito en la subsección 1.2.2:
(a) Verificar su modelo L–E (ecuación (1.130)) aplicando las leyes fı́sicas.
(b) Verificar su modelo L–E (ecuación (1.130)) usando las ecuaciones de Lagrange.
Problema 1.5
La figura 1.17 ilustra un Manipulador Robótico Traslacional Subactuado (MRTS)
de 3 GDL. M1 es un servomotor DC que posee un mecanismo de reducción por
engranajes y un codificador óptico articulado a una polea de radio Rp . Esta polea
usa un cable para transmitir la fuerza F1 para accionar el movimiento de traslación
de un carro de masa Mc montado sobre un par de rieles a lo largo de un eje x. M2 es
también un servomotor con codificador óptimo y mecanismo de reducción, empleado
para accionar el movimiento rotatorio del primer eslabón del MRTS alrededor de
un pivote ubicado en el CG del carro. Este primer eslabón es también articulado a
un segundo eslabón. En esta segunda articulación se monta un codificador óptico de
masa despreciables para detectar la posición angular del segundo eslabón.
En la figura 1.17 θ1 es la posición angular del primer eslabón longitud L1 y
masa m1 , θ2 es la posición angular del segundo eslabón longitud L2 y masa m2 , r es
la posición longitudinal del carro y F es la fuerza de fricción opuesta al movimiento
del carro. El sistema MRTS a ser controlado a lazo cerrado, representa un sistema no
cuadrado porque posee dos entradas de control: los voltajes KA1 u1 y KA2 u2 aplicados
los terminales de las armaduras de M1 and M2 respectivamente, y tres salidas por
controlar: r, θ1 y θ2 . Los parámetros KA1 y KA2 son las ganancias de los amplifi-
cadores. La tabla 1.4 muestra los valores de los parámetros del sistema. El MRTS
1.4 PROBLEMAS 47
L2 Gripper M3
y
θ2
y M 3 Gripper of mass m H
u1 1 (lateral view)
M1 θ1 Link 2
Pulley Pulley
M2 L1
u2 Link 1
x Rp
F F1
r
x1
Problema 1.6
x3
y x
3 M3 3 y M3
Gripper 3
Link L 3
y
q
3
z
3 z3 1111
0000
y
2 Lx3 2
x2
y z2 x2
1 Link
L2 Encoder z2 Encoder
q
2
M y
2 Lx2 1
z1 z1 x1 M2
x
1
Problema 1.7
dci
Qi = A = ke Gei − kP Pi
dt
donde A es la sección del pistón y G es la ganancia del amplificador. Para el sistema
descrito:
(a) Determine su modelo de Lagrange.
(b) Determine su modelo en el espacio de estado (ecuaciones de estado y de salida) .
1.4 PROBLEMAS 49
0110 e1
+ 1010 +
u1
0110 1010 0110
P1 0 0110 y
e1
r1 10 or Electric
10
1 1 e2 actuator
0110- R 1
-
Y1
u2
Return
K1 B1
1010+ +
0110 Pressure
P1 0 0110 y source
r2
0110 2
10
2
Return
0110- R2
1010 -
Y2
P 1 or P2
1010 K2 B2
10
y 1 or y 2
e2
111111111111111
000000000000000
Figura 1.19: Plataformas acopladas por resorte (izquierda) y su actuador hidráulico
(derecha).
Problema 1.8
Considere ahora que los dos actuadores hidráulicos del problema 1.7 se reemplazan
por cilindros neumáticos. Las ecuaciones que describen la dinámica de cada cilindro
neumático son:
ẏi = vi
−kyi − bvi + Api
v̇i =
m
RT Rv psupply ui − RT Rv pi − pi Avi
ṗi = i = 1, 2
Ayi
Problema 1.9
La figura 1.20 ilustra una columna de destilación para fraccionar petróleo pesado. El
modelo con tiempos muertos del fraccionador (ecuación (1.148)), en el cual p es el
operador de Laplace, se describe en [28] y [27]. Tal modelo considera tres variables
que deben de ser controladas: la composición de los productos de la parte superior
50 Modelado de Sistemas No Lineales
PC
T
LC
FC
Upper Top
Reflux Draw
T
Side
Intermediate
Reflux Stripper
T
T
T LC FC
Side
A Draw
LC
Bottoms
Reflux
F T
Feed Bottoms
Problema 1.10
El modelo dinámico del reactor quı́mico con chaqueta de enfriamiento ha sido tomado
de [29] y [27]. La figura 1.21 muestra este reactor en el cual ingresa un fluido liquido
que contiene el producto A. Este liquido se revuelve dentro del tanque mediante
un agitador para formar una mezcla perfecta. El producto A en esta condición va a
experimentar una reacción irreversible exotérmica. Debido a que este tipo de reacción
libera calor, es entonces necesario que la temperatura en el interior del tanque sea
controlada por medio del agua de refrigeración que circula en la chaqueta que rodea
1.4 PROBLEMAS 51
al reactor. La reacción quı́mica del producto o componente A dentro del tanque para
formar el componente B se formula como:
A → B
Product A
Ca 0 Fl Tl 0
11
00 11
00
Reactor
00
11 00
11
00
11 00
11
00
11 00
11
Tc Fc
00
11 Tl
00
11
00
11 00
11
00
11 00
11 Jacket
00
11 00
11
Coolant 00
11 A
Ca
B
Cb 00
11
Fc Tc 0 Products A and B
Ca Cb Fl Tl
describen la dinámica del sistema. Para ello se aplican las leyes de conservación de la
masa y de la energı́a. Para hacer esto, se supone que no existe liquido acumulado en
el reactor, que las concentraciones y temperaturas son homogéneas y que las pérdidas
de energı́a hacia el exterior son despreciables. Las ecuaciones de balance de masa son:
d(V Ca )
= F Ca0 − V kCa − F Ca (1.149)
dt
d(V Cb )
= V kCa − F Cb (1.150)
dt
donde: F Ca0 es el flujo del componente A en kmol/h que ingresa al sistema, F Ca es
el flujo de A que sale del sistema, V es el volumen del liquido en el sistema, −V kCa
es la velocidad de formación de A (el signo negativo indica que el componente A
Ca )
se está consumiendo), d(Vdt es el flujo de A (en kmol/h) en transición, F Cb es el
flujo de B que sale del sistema, +V kCb es la velocidad de formación de B a partir
Cb )
de A (por esta razón posee signo positivo) y d(Vdt es el flujo de B (en kmol/h) en
transición. Las ecuaciones de balance de energı́a se formulan como:
d(V ρ Cp T )
= F ρ Cp T0 − F ρ Cp T − Q + V kCa H (1.151)
dt
d(Vc ρc Cp c Tc )
= Fc ρc Cp c (Tc0 − Tc ) + Q (1.152)
dt
donde: F ρ Cp T0 es el flujo calorı́fico entregado al sistema, F ρ Cp T es el flujo
calorı́fico que sale del sistema, Q es el flujo calorı́fico absorbido por el agua de re-
frigeración, V kCa H es flujo calorı́fico producido en la reacción debido a la entalpı́a
H de la reacción, d(V ρ Cp T )/dt es el flujo calorı́fico en transición (acumulado) en
el interior del tanque, Fc ρc Cp c (Tc0 − Tc ) es el flujo calorı́fico absorbido en el sistema
(refrigeración) y d(Vc ρc Cp c Tc )/dt es el flujo calorı́fico en transición (acumulado) en
la chaqueta del tanque. Todos los flujos calorı́ficos están en kJ/h, mientras que la
entalpı́a H posee unidades de kJ/kmol.
El balance de energı́a en el sistema asume que la temperatura Tc es uniforme en
toda la chaqueta. La transferencia de calor entre el sistema de reacción (que se realiza
a la temperatura T ) y el agua de refrigeración (que se realiza a la temperatura Tc ),
se describe mediante la relación:
Q = U S(T − Tc ) (1.153)
Tabla 1.7: Parámetros nominales del reactor quı́mico con chaqueta de enfriamiento.
Problema 1.11
Los ángulos ψ, θ y φ que caracterizan la posición del avión con respecto a los
54 Modelado de Sistemas No Lineales
donde w∗ es el vector velocidad angular expresado con respecto a los ejes del viento
y M es una matriz de la forma:
0 sin φ sec θ cos φ sec θ
M(ψ, θ, φ) = 0 cos φ − sin φ
1 sin φ tan θ cos φ tan θ
J es la matriz de inercia:
Ix 0 −Ixz
J= 0 Iy 0
−Ixz 0 Iz
V̇ = −(D/m) − g sin θ
α̇ = q − q ∗ sec β − (p cos α + r sin α) tan β
β̇ = r∗ + p sin α − r cos α
1.4 PROBLEMAS 55
en donde aij ’s y los bij ’s son los parámetros aerodinámicos fijos que dependen de la
geometrı́a del avión, de la densidad del aire, etc., y δa , δe y δr denotan las deflexiones
del alerón, del elevador y del timón. El vector columna con elementos D, L y S se
puede expresar como:
D c11 V 2 + c12 V 2 cos α − cos α cos β
L = c21 V 2 + c22 V 2 sin 2α + P sin α δP
S 2
c31 V sin 2β cos α cos β
donde los cij ’s son parámetros de ganancia fija, P indica el máximo empuje y δp la
posición del acelerador. Las ecuaciones presentadas describen un sistema cuyo estado
está definido en una cierta vecindad abierta U de R9 , sujeto a la acción del vector de
entrada u con elementos δa , δe , δr y δP .
Problema 1.12
56 Modelado de Sistemas No Lineales
I Ls Φ
+ Rr Lr
Vr
Ω Rs
- I
- Vs +
dIs
Ls + Rs If = Vs
dt
donde Is es la corriente del estator, Vs es el voltaje del estator, Rs y Ls son la re-
sistencia e inductancia del devanado del estator respectivamente. El balance eléctrico
en el devanado del rotor resulta:
dIr
Lr + Rr Ir = Vr − E
dt
donde Ir es la corriente del rotor, Vr es el voltaje del rotor, Rr y Lr son la resistencia
e inductancia del devanado del rotor respectivamente, y E es el voltaje de fuerza
contra electromotriz. Considerando sólo fricción viscosa, el balance mecánico en la
carga está dado por:
dΩ
J + FΩ = T
dt
donde Ω representsa la velocidad angular del eje del motor, J denota la inercia de la
carga, F es la constante de fricción viscosa, y T es el torque de carga. Si Φ denota el
flijo asociado con el devanado del estator, entonces:
E = Ke ΦΩ T = Km ΦIr Φ = Ls Is
Problema 1.13
Problema 1.14
Problema 1.15
Problema 1.16
La figura 1.26 muestra un servomotor con carga no lineal de dos grados de libertad
(2GDL). Notar que la carga posee dos grados de libertad. La unión de la carga
con el eje del servomotor no es flexible. Asumir que la inductancia de armadura no
es despreciable. Formular el modelo dinámico no lineal del sistema empleando las
ecuaciones de Lagrange y las ecuaciones de la fı́sica.
Mo Ro
L1 τ
m1
’
Union
R θm no
+ Ν1 flexible
Lo θ
u Va eb
m
-
Jm
L Bm θ
Ν2 BL JL
La figura 1.27 muestra el proceso péndulo doble no lineal. El acoplamiento entre los
dos péndulos no es flexible. Las varillas poseen longitud L1 y L2 . Despreciando la in-
ductancia de armadura, formular el modelo dinámico no lineal del sistema empleando
las ecuaciones de Lagrange y las ecuaciones de la fı́sica.
+ Fuerza de
u control
- τ
Servomotor
D.C. y
θ
z F
Carro
La figura 1.28 muestra el proceso doble grúa puente no lineal. El acoplamiento entre
las dos secciones de la varilla no es flexible. Las varillas poseen longitud L1 y L2 .
Despreciando la inductancia de armadura, formular el modelo dinámico no lineal del
sistema empleando las ecuaciones de Lagrange y las ecuaciones de la fı́sica.
+ Fuerza de
u control
-
Servomotor
D.C. y
z F
Carro
θ
z
Doble
y’ ’ puente
τ grua
Figura 1.30: Sistema de control en tiempo real para el manipulador robótico de 1GDL.
en tiempo real para controlar la posición del brazo robótico. Tales programas se
encuentran en el CD adjunto a este trabajo.
El algoritmo de control grabado en la PC procesa tal medición y genera una
señal de control. Esta señal resulta débil para poder alimentar directamente la ar-
madura el servomotor D.C. Por ello es que se procesa empleando la técnica de mod-
ulación por ancho de pulsos (en inglés PWM: Pulse Width Modulation) y luego pasa
a un amplificador con configuración en H para darle el nivel de potencia adecuado
para hacer girar al servomotor en la dirección adecuada.
H, obtienen el nivel adecuado de potencia para conmutar el giro del motor tal como
se explica en la siguiente sección. La subsección ?? muestra el código en lenguaje
ensamblador del programa para generar las señales PWM, mientras que la figura
1.33 muestra el diagrama circuital para realizar la placa del generador de pulsos
PWM empleando ORCAD.
Figura 1.33: Diseño ORCAD para implementar la placa del generador de pulsos
PWM.
en H es que sólo necesita una fuente. Los dispositivos conmutadores son conmutados
en pares generando una tensión bipolar a la salida. Para disparar los MOSFET se
necesita circuiterı́a adicional para generar la tensión de disparo en cada conmutador.
En sı́ntesis, cuando el sistema de disparo cierra el conmutador (A1 A4 ) y abre
(A2 A3 ), el sentido de la corriente es la lı́nea punteada, induciendo de esta forma una
tensión +Vcc en el servomotor. Luego, si el sistema de disparo abre al conmutador A1,4
y cierra A2,3 , el sentido de la corriente es la lı́nea continua, induciendo ası́ una tensión
−Vcc en el servomotor. Por lo tanto, el servomotor ve en sus terminales una onda de
voltaje cuadrada variando entre ±Vcc y la corriente que pueda absorber dependerá de
la capacidad en corriente de los dispositivos usados en los conmutadores y de la fuente.
la figura 1.34.
La figura 1.35 muestra los circuitos de disparo de los MOSFET Q14, Q5, Q6
y Q7. Tales circuitos emplean compuertas lógicas de tecnologı́a TTL (Transistor
Transistor Logic) conformado por los integrados 7404, 7486, 7408 y dos resistencias de
500 Ω. El circuito de disparo de Q14 comprende los transistores Q1 y Q3 formando un
amplificador que trabaja en clase B. Este circuito de disparo trabaja simultáneamente
con el circuito de disparo del MOSFET Q6 (transistores Q8 y Q9). Del mismo modo,
los circuitos de disparo de los MOSFET Q5 (transistores Q2 y Q4) y Q7 (transistores
Q10 y Q11) trabajan en conjunto.
Por ejemplo para el MOSFET Q14, cuando en la salida de la compuerta 7417
existe una señal ON, el transistor Q1 se corta mientras que Q3 se satura. De este
modo a la base del MOSFET llega una corriente que lo habilita. En caso contrario,
cuando en la salida de la compuerta 7417 existe una señal OFF, el transistor Q1 se
satura mientras que Q3 se corta. De este modo a la base del MOSFET no llega una
corriente y lo inhabilita. El mismo análisis se puede hacer para todos los circuitos de
disparo.
Control Óptimo
Este capı́tulo trata el problema del control óptimo cuadrático gaussiano, denomi-
nado ası́ porque el ı́ndice de rendimiento o función de costo que emplea es una fun-
ción cuadrática de los estados y de las señales de control. La solución del problema de
control planteado consiste en determinar un extremo de la función de costo mediante
minimización con el propósito de generar la ley de control óptima.
La configuración del sistema de control óptimo no lineal desarrollado en este capı́tu-
lo comprende el modelo no lineal multivariable del sistema a controlar, un observador no
lineal para estimar los estados del sistema y un controlador de realimentación de esta-
dos del tipo proporcional–integral. El diseño del sistema de control óptimo consiste en
producir una fuerza de control u que sea capaz de hacer que el vector de salida y del
sistema (la salida controlada) siga al vector de referencias deseadas r cumpliendo ciertas
especificaciones de diseño, no obstante la presencia de incertidumbres en los parámetros
y de disturbios estocásticos gaussianos actuando sobre el sistema en operación.
Nonlinear state
x^ observer
r -
z y
+
ẋ = f (x, u) + v x(0− ) = 0
y = h(x) + w (2.1)
x1 u1 v1 y1 w1
.. ..
x = ... u = ... v= . y = ... w= .
xn um vn yp wp
f1 (x, u) h1 (x)
.. ..
f = . h= .
fn (x, u) hp (x)
donde x es el vector de estado, u es el vector de control (la entrada al sistema) e y
es el vector de salida. Asumiremos que las funciones vectoriales de variable vectorial
f (x, u) y h(x) son operadores no lineales diferenciables que representan la dinámica
del sistema y la dinámica de la salida respectivamente. Los vectores v de orden n × 1
y w de orden p × 1 son los disturbios en los estados del sistema (ruidos en los estados)
y en sus salidas (ruidos de medición) respectivamente.
Sean x0 , u0 e y0 los vectores de referencia (o nominales) correspondientes a x,
u e y respectivamente. Si la entrada del sistema se selecciona exactamente igual a
u0 , su respuesta será x0 . Entonces x0 satisface:
x = x0 + δx u = u0 + δu y = y0 + δy (2.3)
2.3 El Controlador PI de Realimentación de Estados 69
donde los vectores residuales δx, δu y δy representan desviaciones con respecto a los
correspondientes vectores de estado, de control y de salida respectivamente. Reem-
plazando (2.3) en (2.1) produce:
Como hemos asumido que las desviaciones actuales son pequeñas, entonces el sistema
(2.4) admite ser linealizado alrededor de un vector de trayectorias nominales o de
referencia x0 . La expansión de (2.4) en series de Taylor alrededor de x0 y u0 resulta:
donde:
∂f1 ∂f1 ∂f1 ∂f1
∂x1 ··· ∂xn ∂u1 ··· ∂um
.. .. .. .. .. ..
A = . . . B= . . .
∂fn ∂fn ∂fn ∂fn
∂x1 ··· ∂xn (x0 ,u0 ) ∂u1 ··· ∂um (x0 ,u0 )
∂h1
∂x1 ··· ∂h1
∂xn
.. .. ..
C = . . . (2.6)
∂hp ∂hp
∂x1 ··· ∂xn (x0 )
Observar que las matrices jacobianas Ann , Bnm , Cpn y Dpm necesitan ser evaluadas
alrededor de los vectores de referencia x0 y u0 . Para resolver el problema del control
óptimo en consideración, se requiere diseñar una apropiada ley de control de reali-
mentación de estados de la forma u = −Gx, donde G de dimensión m × n, es la
la matriz de ganancia de realimentación de estados, la cual se diseña empleando el
siguiente sistema linealizado:
u x y
g h
+
+
f(.)
f(x)
u
-D
G = R−1 BT K (2.12)
δ ż = δy = Cδx + w (2.14)
δ ẋa = Aa δxa + Ba δu + va
δy = Ca δxa + w (2.15)
2.3 El Controlador PI de Realimentación de Estados 71
x
u x y z
g h z
+
+
f(.)
f(x) x
u a z
-D
cumpla la condición de controlabilidad dada en (2.11), que para este caso es:
rango Ma = rango Ba Aa Ba · · · (Aa )n−1 Ba = n + p (2.18)
α ≥ 0 es una constante de peso exponencial. Observar que las matrices para ponderar
el rendimiento del sistema (las matrices de sintonización) son R, Q y Qa .
De acuerdo al teorema de la separación [12], los controladores de realimentación
de estado se pueden implementar empleando los estimados de los estados en lugar de
los estados actuales del sistema. Por consiguiente, las leyes de control (2.8) y (2.21)
se pueden implementar también como:
δu = −Gδ
x
δ
x
xa = −
δu = −Ga δ Gx Gz (2.23)
δ
z
v + w
+
u + x + y z
g h
+
+ -
f(.) r=0
f(x)
u ^
x
-D a z^ = z
z^
^
x ^
y y
u g C
+ + - +
+ +
f(.)
f(x)
E
Figura 2.4: Estructura del sistema de control óptimo para estimación de estados y
control.
donde y(t) es la salida real del sistema que debe de seguir a la trayectoria de referencia
r(t) cumpliendo ciertas especificaciones de diseño previamente establecidas. Es claro
que cuando se cumpla el objetivo de control del sistema, es decir, cuando δy(t) ∼ = 0,
entonces: y(t) ∼
= r(t).
El procedimiento de diseño del sistema de control óptimo desarrollado sigue los
pasos siguientes:
Ejemplo 2.1 Design and Simulation of an Optimal Control System for the
EJRM.
IN PROGRESS
Ejemplo 2.2 Design and Simulation of an Optimal Control System for the
TRM.
IN PROGRESS
0 0
0 0 1 0 0 0
B= 1 C=
RT 1 (J1 +a1 +a2 ) 0 0 1 0 0
1
0 RT 2 (J2 +a2 )
Por consiguiente:
a A 0 a B
A = B = Ca = C 0
C 0 0
Las matrices de sintonización se pueden elegir como: Q = 4I, R = 10I, Qa = 0,01I,
V = 0,1I, W = 0,1I. La ganancia Da del controlador y la ganancia H del observador
fueron determinadas usando los comandos de MATLAB lqr y lqe, respectivamente,
usando el hecho de que el sistema es completamente controlable y completamente
observable. El parámetro α del ı́ndice de rendimiento se fijó en 5. La simulación se
llevó a cabo ejecutando el archivo srm4opt.m (que se encuentra en el CD adjunto),
para las siguientes trayectorias deseadas: xd1 (t) = sin t + 0,1t y xd2 (t) = cos t. La
t
relación z = 0 ydτ se aproximó en el dominio discreto como z(k + 1) = z(k) + T x(k),
donde T = 0,0014 s es el tiempo de muestreo y k = t/T es el tiempo discreto. Los
resultados de la simulación se muestran en las figuras 2.5 y 2.6.
2.5
2
Base trayectory (rad)
1.5
0.5
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s)
2
Control signal (volt)
−1
−2
−3
−4
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s)
PROBLEMAS
Problema 2.1
Problema 2.2
Diseñe un sistema de control óptimo cuya salida sea capaz de seguir a una trayectoria
determinada a las plataformas acopladas con actuadores neumáticos descrito en el
problema 1.8 y mostrado en la figura 1.19. Asumir valores mı́nimo y máximo de las
masas de las plataformas m1 and m2 .
Problema 2.3
Diseñe un sistema de control óptimo cuya salida sea capaz de seguir a una trayectoria
determinada, para la plataforma inercial descrita en el problema ?? e ilustrada en la
figura ??. Asumir valores mı́nimo y máximo de ... IN PROGRESS
Problema 2.4
Design a robust nonlinear optimal control system for the heavy oil fractionator de-
scribed in the problem 1.9 and depicted in figure 1.20. Assume minimum and maxi-
mum values of the ... IN PROGRESS
Problema 2.5
Design a robust nonlinear optimal control system for the chemical jacket reactor
described in the problem 1.10 and depicted in figure 1.21. Assume minimum and
maximum values of the ... IN PROGRESS
Problema 2.6
Design a robust nonlinear optimal control system for the high purity distillation
column described in the problem ?? and depicted in figure ??. Assume minimum
and maximum values of the ... IN PROGRESS
1.5
Arm trajectory (rad)
0.5
−0.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s)
15
Control signal (volt)
10
−5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time (s)
Problema 2.7
Design a robust nonlinear optimal control system for the aircraft described in the
problem 1.11 and depicted in figure 1.22. Assume minimum and maximum values of
the ... IN PROGRESS
Problema 2.8
Design a robust nonlinear optimal control system for the spacecraft described in the
problem ?? and depicted in figure ??. Assume minimum and maximum values of the
... IN PROGRESS
Problema 2.9
Design a robust nonlinear optimal control system for the nonlinear DC motor de-
scribed in the problem 1.12 and depicted in figure 1.23. Assume minimum and max-
imum values of the ... IN PROGRESS
Capı́tulo 3
Reference
model
r -
Adaptive +
Process
controller
Adaptation
mechanism
ẋ = f (x, t) (3.1)
se dice que es no autónomo si f depende del tiempo, por ejemplo, si f posee parámetros
variantes con el tiempo. por consiguiente, un sistema autónomo puede ser descrito
por: ẋ = f (x). Las trayectorias de estado para procesos autónomos son independientes
del tiempo inicial, mientras que para los no autónomos generalmente no lo son.
Un estado o punto de equilibrio xe (realmente un vector constante) de un
sistema autónomo se puede determinar de:
0 = f (xe ) (3.2)
||x|| = R en sı́. La región esf0rica anular cerrada r ≤ ||x|| ≤ R será denotada co-
mo BrR . Asumiremos que en una cierta región esférica Ω : ||x|| < B(R), todas las
derivadas parciales ∂xi /∂xj existen y son continuas en Ω. Entonces diremos que el
origen es:
estable si alguna trayectoria que empieza en B(r) en un punto arbitrario xo
nunca logra alcanzar la esfera frontera S(R) de B(R);
asintóticamente estable cuando es estable y en adición cada trayectoria de es-
tado que empieza en B(R) en un punto arbitrario xo , tiende hacia le origen
conforme el tiempo se incremente indefinidamente;
inestable cuando para algún R y r, ya sea grande o pequeño, alguna trayectoria
de estado que empieza en B(R) en un punto arbitrario x0 logra alcanzar la
esfera frontera S(R). Observar en la figura 3.2 que para la trayectoria T el
origen es inestable en el sentido de Lyapunov, a pesar de que tal trayectoria de
estado muestre convergencia.
Asymptotically
stable
Stable
B(r) S(R)
0 r
x0 R
Unstable
H(A)
B(R) T
S(A)
(b) V (0) = 0
(c) Fuera del origen, pero siempre en Ω, V (x) es positiva. Por consiguiente, el origen
es un mı́nimo aislado de V (x).
Ejemplo 3.1
Ejemplo 3.2
3.2 Estabilidad vı́a el Método Directo de Lyapunov 83
V(x1 , x 2) x2
V(x1 , x 2)
x2 x1
x1
R i
+
+
u C e RN
- -
Solución: Sumando las corrientes que salen del nodo superior derecho (ver la figura
3.4) produce:
e−u
C ė + + N e3 = 0
R
La energı́a almacenada en el capacitor está dada por V (e) = 12 Ce2 . Para el sistema
no actuado (u = 0) tenemos:
2 1 2
V̇ (e) = Ceė = −e + Ne ≤ 0
R
Por consiguiente, la función V (e) es una función de Lyapunov.
Ejemplo 3.3
K0 B
K1 b(t)
111111111111
000000000000
Figura 3.5: Sistema masa–amortiguador–resorte.
dV ∂V ∂V ∂V ∂V
V̇ = = + ẋ = + f (x, t)
dt ∂t ∂x ∂t ∂x
Se dice que una función V (x, t) es definida positiva si satisface las condiciones
(a)–(c). También se dice que V (x, t) es una función definida positiva si esta función
domina a otra función definida positiva W (x), por ejemplo, cuando V (x, t) ≥ W (x).
También, V (x, t) es definida negativa si −V (x, t) es definida positiva; V (x, t) es
semidefinida positiva si esta función domina a otra función semidefinida positiva
W (x); V (x, t) es semidefinida negativa si −V (x, t) es semidefinida positiva.
Ejemplo 3.4
3.2 Estabilidad vı́a el Método Directo de Lyapunov 85
V=0
. . 0
Ω
x0
V = constant
Ω1
VII. Teorema de Inestabilidad. Sea Ω una vecindad del origen. Sea Ω una
región en Ω. Dada una funci2n V (x, t) en Ω, entonces el equilibrio del origen en el
tiempo t0 es inestable si:
Ejemplo 3.5
Determine la estabilidad del sistema autónomo descrito en los ejemplos 3.2 y 3.3
aplicando el método directo de Lyapunov.
Solución: Del ejemplo 3.2 podemos establecer que la función de Lyapunov V (e) =
1 2
2 Ce → ∞ cuando ||e|| → ∞. Por consiguiente, el sistema autónomo no lineal
descrito en dicho ejemplo es completamente (globalmente) asintóticamente estable.
Del ejemplo 3.3 podemos establecer que la función de Lyapunov:
1 1 1
V (x) = M ẋ2 + K0 x2 + K1 x4 → ∞
2 2 4
cuando ||x|| → ∞. Por consiguiente, el sistema autónomo no lineal descrito en dicho
ejemplo es completamente (globalmente) asintóticamente estable.
Ejemplo 3.6
conforme ||x|| → ∞, siempre que α > 0, b(t) > M α, y ḃ(t) < 2K0 . Por consigu-
iente, el sistema no autónomo descrito en tal ejemplo es completamente (globalmente)
asintóticamente estable.
V̇ (x) ≤ 0 (3.5)
para todo x en Ω . Sea R el conjunto de todos los puntos dentro de Ω donde V̇ (x) = 0,
y sea M el más grande conjunto invariante en R. Entonces, cada solución x(t) en Ω
tiende a M conforme t → ∞.
V
V=l
Ωl
x2
x0
x1
Ejemplo 3.7
Observe que V̇ (x) < 0 dentro de un cı́rculo de radio 2. Por consiguiente, usando el
teorema de estabilidad II de Lyapunov, podemos inferir que el origen es asintótica-
mente estable. Para = 4, la región Ω definida por V (x) = x21 +x22 < 4 es acotada y el
conjunto R es el origen, el cual es un conjunto invariante. Por consiguiente, cualquier
trayectoria que se inicia dentro del cı́rculo de radio 2 converge hacia el origen y esta
región constituye el dominio de atracción.
m
dn y
(n)
y + αi∗ fi (x, t) = bu y (n) (3.7)
dtn
i=1
m
1 αi∗
h y (n) + αi fi (x, t) = u h= αi = (3.8)
b b
i=1
90 Control Adaptativo con Modelo Referencial
e = y − yd
donde:
∆(p) = pn−1 + λn−2 pn−2 + · · · + λ0 = (p − p1 ) . . . (p − pn )
(n−1)
yr(n−1) = yd − λn−2 e(n−2) − · · · − λ0 e
en el cual p es el operador de Laplace y ∆(p) es un polinomio estable Hurwitz, lo
cual significa que todas las raı́ces complejas pi = σi + jωi , i = 1, . . . , n de ∆(p) = 0
verifican la condición σi < 0. Asumamos la siguiente ley de control:
m
u = h yr(n) − k s + αi fi (x, t) (3.10)
i=1
(n)
donde las constantes k y h poseen el mismo signo. La variable yr , denotada como
la variable de referencia de y (n) , se determina de:
(n−1)
(n) dyr
yr(n) = yd − λn−2 e(n−1)
− ··· − λ0 ė yr(n) =
dt
Sustituyendo (3.10) en (3.8) nos conduce al error de seguimiento dinámico:
h ṡ + k s = 0 (3.11)
˙
h = −γ sgn(h)s yr(n)
α̂˙ = −γ sgn(h)s fi (3.15)
Usando el hecho de que h sign(h) = |h| y k sign(k) = |k| y dado que k y h poseen el
mismo signo, se puede demostrar fácilmente que:
donde las ganancias de adaptación γh y γk son diferentes para cada parámetro des-
conocido. Seleccionando la siguiente candidata para función de Lyapunov:
m
V = |h|s2 + γh−1%
h2 + γk−1 %i2
α (3.19)
i=1
se puede demostrar fácilmente que V̇ nos conduce a (3.17), lo cual garantiza conver-
gencia de seguimiento global del SCAMR.
por:
&
(n)
˙
−γh sgn(h)s yr |s| > ∆y
h =
0 |s| < ∆y
−γα sgn(h)s fi |s| > ∆f
˙ i =
α (3.20)
0 |s| < ∆f
Ejemplo 3.8
donde:
Jeq La Beq La Ra Jeq
h= α1 = +
KA KA nKA Km
nKb Ra Beq meq gLLa meq gLRa
α2 = + α3 = α4 =
KA nKA Km KA nKA Km
1
meq = m + mH f1 = ÿ f2 = ẏ f3 = ẏ cos y f4 = sen y
2
3.4 SCAMR para Sistemas No Lineales Multivariables 93
donde:
Jeq La Beq La Ra Jeq
h= α1 = +
KA KA nKA Km
nKb Ra Beq meq gLLa meq gLRa
α2 = + α3 = α4 =
KA nKA Km KA nKA Km
1
meq = m + mH f1 = ÿ f2 = ẏ f3 = ẏ cos y f4 = sen y
2
donde las matrices M, P y d representan la inercia del sistema a controlar, los torques
centrı́petos y de Coriolis, y los torques gravitacionales respectivamente. También, q
es el vector de coordenadas generalizadas y u es el vector de control.
El objetivo de control del SCAMR es diseñar una ley de control u capaz de
hacer que la salida del sistema q(t) siga a la taryectoria deseada qd (t) con velocidad
suficiente a pesar de la presencia de incertidumbres en los parámetros. Asumamos
que todos los términos en (3.21) dependen linealmente de un vector de parámetros a
con elementos conocidos, a saber:
Ya = u (3.22)
donde Y es una matriz conocida. Considere la siguiente ley de control:
a − KD s
u = Y (3.23)
d
si = ( + λi )n−1 q%i = (p + λi )n−1 q%i (3.24)
dt
donde λi > 0 es una constante (el ancho de banda), p es el operador de Laplace y:
q%i = qi − qdi
%̇ + Λ%
s=q q = q̇ − q̇r q̃ = q − qd q̇r = q̇d − Λ%
q (3.25)
donde:
1 d T
[q̇ Mq̇] = q̇T (u − d) (3.27)
2 dt
donde q̇T Mq̇ es la energı́a cinética del sistema y q̇T (u − d) es la potencia de entrada
generada por el actuador. Diferenciando el miembro izquierdo de (3.27):
1
q̇T Mq̈ + q̇T Ṁq̇ = q̇T (u − d) (3.28)
2
De (3.21) obtenemos: Mq̈ = u − d − Pq̇. Sustituyendo este término en (3.28):
Ha sido establecido en [24], [3] que (Ṁ − 2P) = J es antisimétrica. Por consiguiente:
Ṁ = 2P + J (3.30)
1 T
˙ Γ−1 a
V̇ (t) = sT Mṡ + sT Ṁ + a %
2
3.4 SCAMR para Sistemas No Lineales Multivariables 95
T
˙ Γ−1 a
a − KD s − Mq̈r − Pq̇r − d) + a
V̇ (t) = sT (Y % (3.32)
˙ = −ΓYT s
a (3.33)
Dado que s = q %̇ + Λ%q (ver (3.25)), entonces (3.35) garantiza que los errores de
seguimiento de posición q% y de velocidad q%̇ tiendan a 0 conforme t → ∞. En otras
palabras, los errores de seguimiento convergen en la superficie s = 0.
El resultado establecido en (3.35) es también válido para el modelo L–E dado
en (1.130). Para este caso, las relaciones (3.22), (3.23) y (3.34) toman las formas
siguientes:
Ya = T (3.36)
a − KD s
u = Y (3.37)
0.5
Cart position (m)
−0.5
−1
0 1 2 3 4 5 6 7 8 9 10
Time (s)
3
Control voltage u1
−1
−2
0 1 2 3 4 5 6 7 8 9 10
Time (s)
0.5
−0.5
−1
−1.5
0 1 2 3 4 5 6 7 8 9 10
Time (s)
6
Control voltage u2
−2
−4
−6
0 1 2 3 4 5 6 7 8 9 10
Time (s)
0.04 0.05
0
0.02
ae1
ae2
−0.05
0
−0.1
−0.02 −0.15
0 2 4 6 8 10 0 2 4 6 8 10
0.2 0.05
0.1
ae3
ae4
0 0
−0.1
−0.2 −0.05
0 2 4 6 8 10 0 2 4 6 8 10
0.6 0.02
0.4
0.01
ae5
ae6
0.2
0
0
−0.2 −0.01
0 2 4 6 8 10 0 2 4 6 8 10
Time (s) Time (s)
Definiendo:
donde:
Y11 = q̈r1 Y12 = q̈r1 cos2 q2 − q̇r1 q̇2 sin q2 cos q2 − q̇r2 q̇1 sin q2 cos q2
Y13 = 0 Y21 = 0 Y22 = q̈r2 + q̇r1 q̇1 sin q2 cos q2 Y23 = cos q2
La ley de control está dada por (3.36):
a − KD s
T = Y
La simulación del SCAMR diseñado (archivo srmmrac.m) asume las siguientes condi-
ciones iniciales: q1 (0) = 0 rad (posición de la base) y q2 (0) = −π/2 rad (posición del
brazo). Para un tiempo de muestreo de T = 0,01 s, los parámetros de la ley de control
se fijaron en KD = diag[1,5]. Las trayectorias deseadas están dadas por:
PROBLEMAS
IN PROGRESS
3.4 SCAMR para Sistemas No Lineales Multivariables 99
1.5
−0.5
−1
−1.5
0 1 2 3 4 5 6 7 8 9 10
6
Voltage control u1
−2
−4
0 1 2 3 4 5 6 7 8 9 10
Time (s)
1.5
1
Link position (rad)
0.5
−0.5
−1
−1.5
−2
0 1 2 3 4 5 6 7 8 9 10
10
8
Voltagr control u2
−2
−4
0 1 2 3 4 5 6 7 8 9 10
Time (s)
0.15
0.1
ae1
0.05
−0.05
0 1 2 3 4 5 6 7 8 9 10
0.6
0.4
ae2
0.2
−0.2
0 1 2 3 4 5 6 7 8 9 10
1.3
ae3
1.25
0 1 2 3 4 5 6 7 8 9 10
Time (s)
Linealización por
Realimentación de Estados
∂h ' (
∇h(x) = = ∂h
∂x1 ··· ∂h
∂xn (4.2)
∂x
La derivada de Lie de una función escalar h(x) con respecto al campo vectorial
f (x) es una nueva función escalar Lf h definida como:
∂h ∂h
Lf h = ∇h f = f1 + · · · + fn (4.4)
∂x1 ∂xn
Observar que la derivada de Lie es el producto interno entre ∇h(x) y f (x). Derivadas
de Lie repetidas se pueden formular en forma recursiva como:
L0f h = h
f h) = ∇(Lf h) f
Lif h = Lf (Li−1 i−1
for i = 1, 2, . . . (4.5)
Lg Lf = ∇(Lf ) g (4.6)
Corchetes de Lie
adf0 g = g
adfi g = [f , adfi−1 g] para i = 1, 2, . . . (4.8)
1) Bilinealidad:
2) Anti-conmutatividad:
[f , g] = −[g, f ] (4.10)
3) Identidad de Jacoby:
Ladf g h = Lf Lg h − Lg Lf h (4.11)
4.1 Herramientas de la Geometrı́a Diferencial 103
lo cual significa que estos campos se pueden representar como una combinación
lineal.
2) Un conjunto que contiene un único campo vectorial f es involutivo porque:
[f , f ] = 0.
3) El conjunto de campos vectoriales f1 , f2 , . . ., fm es involutivo si ∀ x y ∀ i, j
' ( ' (
rank f1 (x) . . . fm (x) = rank f1 (x) . . . fm (x) [fi (x), fj ](x)
'(
donde la notación rank . representa el rango de una matriz de campos vecto-
riales columna.
1
u= (−Lf h(x) + v)
Lg h(x)
d iy
y (i) = = Lif h(x) + Lg Li−1
f h(x)u
dt i
Lg Lr−1
f h(x) = 0
1
u= r−1 (−Lrf h(x) + v)
Lg Lf h(x)
en:
y (r) = Lrf h(x) + Lg Lr−1
f h(x)u
donde:
donde v es una nueva entrada a ser determinada, a(z) y b(z) están dadas en (4.23),
y z es el nuevo estado para linealización. Substituyendo (4.24) en (4.22) produce la
siguiente representación canónica invariante con el tiempo:
ż = Mz + Nv y = Cz (4.25)
donde:
0 1 0 ··· 0 0
0 0 1 ··· 0 0
.. .. .. .. .. ..
M= . . . . . N= . C 1 0 ··· 0 0
0 0 0 0 1 0
0 0 0 0 0 1
d ny
v= = y (n) = ρ(n) + K1 ȳ (n−1) + · · · + Kn−1 ȳ˙ + Kn ȳ (4.27)
dt n
m
ẋ = f (x) + G(x)u = f (x) + gj (x)uj (x) y = h(x) (4.29)
j=1
x1 f1 (x) u1
x = ... f (x) = ..
. u = ...
xn fn (x) um
y1 h1 (x)
y = ... h(x) = ..
.
ym hm (x)
G11 (x) . . . G1m (x)
..
G(x) = . ... = g1 (x) · · · gm (x)
Gm1 (x) . . . Gmm (x)
Observando las filas de A, podemos postular que cada entero ri está relacionado con
la i-ésima salida del sistema. También, notar que ri es el número de diferenciaciones
ejecutadas en la salida yi requeridas para que aparezca al menos uno de los compo-
nentes del vector de entrada u. La no singularidad de A(x) en x = x0 es la versión
multivariable de la condición impuesta por la ecuación (4.20).
Ejemplo 4.1
La descripción en el espacio de estado del sistema MRE (manipulador robótico esféri-
co) de dos entradas y dos salidas se da en (1.146) y (1.147). Determinar el grado
relativo total del MRE.
r = r1 + r2 + · · · + rm ≤ n
Un sistema MIMO descrito por (4.29) posee linealización exacta, si su grado relativo
total r = r1 + . . . + rm es igual al orden n del sistema, es decir, la dimensión n del
vector de estados. Para linealización exacta, la transformación de estados:
φ1 (x) z1
.. ..
. .
φr1 (x) z
r 1
ψ1 (x) zr1 +1
.. ..
. .
z = Φ(x) = ψr (x)
=
(4.31)
2 zr1 +r2
.. ..
. .
ξ1 (x) zr +r +···+1
1 2
.. ..
. .
ξrm (x) zr1 +r2 +···+rm
donde
(r )
m
y1 1 = φ̇r1 (x) = Lrf 1 h1 (Φ−1 (z)) + Lgj Lfr1 −1 h1 (Φ−1 (z))uj
j=1
(r )
m
y2 2 = ψ̇r2 (x) = Lrf 2 h2 (Φ−1 (z)) + Lgj Lfr2 −1 h2 (Φ−1 (z))uj
j=1
..
.
110 Linealización por Realimentación de Estados
Solución: Del ejemplo 4.1, la condición para linealización exacta se cumple para el
sistema MRE porque r = 6 = n (orden del sistema). Aplicación de las relaciones
(4.31) y (4.32) en el sistema MRE nos conduce a:
z1 h1 x1 y1 z1
z2 Lf h1 (x) x3 ẏ1 z4
z3 L2f h1 (x) x5 ÿ1 z2
z=
= Φ(x) =
=
=
x = Φ (z) =
−1
z4 h2 (x) x2 y2 z5
z5 Lf h2 (x) x4 ẏ2 z3
z6 L2f h2 (x) x6 ÿ2 z6
(4.33)
ż1 z2
ż2 z3
ż3 v1 y1 z1
ż =
=
= Mz + Nv
y= = = Cz (4.34)
ż4 z5 y2 z4
ż5 z6
ż6 v2
donde:
0 1 0 0 0 0
0 0 1 0 0 0
v1 0 0 0 0 0 0
v= M=
v2 0 0 0 0 1 0
0 0 0 0 0 1
0 0 0 0 0 0
0 0
0 0
1 0 1 0 0 0 0 0
N=
C=
0 0 0 0 0 1 0 0
0 0
0 1
4.3 Linealización por Realimentación: Caso MIMO 111
m
yiri = Lrf i hi (Φ−1 (z)) + Lgj Lfri −1 hi (Φ−1 (z))uj
j=1
donde v = [v1 . . . vm ]T es una nueva entrada por determinar, y A(z) (ver (4.30)) y
B(z) son matrices con elementos:
r −1
aij = Lgi Lf j hj (Φ−1 (z)) i, j = 1, . . . , m
(r )
lo cual significa que yi i = d ri yi /dt ri = vi , i = i, . . . , m, donde cada entrada vi se
diseña para ubicar los polos del correspondiente subsistema lineal equivalente. Tal
subsistema se puede seleccionar como:
d r i yi (r ) (r ) (r −1)
vi = r
= yi i = ρi i +Ki,1 ȳi i +· · ·+Ki,ri −1 ȳ˙ i +Ki,ri ȳi i = 1, . . . , m (4.37)
dt i
donde ρi (t) es la i-ésima trayectoria deseada y, ȳi (t) = ρi (t) − yi (t) es la i-ésima señal
de error de seguimiento. La i-ésima ecuación caracterı́stica del subsistema resulta:
(ri ) (ri −1)
ȳi + Ki,1 ȳi + · · · + Ki,ri −1 ȳ˙ i + Ki,ri ȳi = 0 i = 1, . . . , m (4.38)
Ejemplo 4.3
112 Linealización por Realimentación de Estados
(3) (3)
ż3 = z1 = y13 = ẋ5 = q1 = v1
(3) (3)
ż6 = z2 = y23 = ẋ6 = q2 = v2 (4.39)
(3)
v1 = ρ1 + K13 ρ̈12 + K12 ρ̇1 + K11 ρ1 − K13 ÿ12 − K2 ẏ1 − K11 y1
(3)
v2 = ρ2 + K23 ρ̈2 + K22 ρ̇2 + K21 ρ2 − K23 ÿ2 − K22 ẏ2 − K21 y2 (4.40)
u1 = v1 LT1 (H11 + J1 ) + P1
u2 = v2 LT2 (H22 + J2 ) + P2 (4.41)
El diseño de controladores por ubicación de polos para sistemas lineales [13], [14]
asume que todos los estados del sistema están disponibles en la realimentación. En
este caso, es posible forzar para que el sistema posea polos predeterminados a lazo
cerrado, es decir, polos prescritos (o eigenvalores) en localizaciones deseadas. En la
práctica, sólo parte de los estados son disponibles para procesamiento, y además, no
debemos de dejarnos tentar en la diferenciación de variables de estado con el fin de
generar una nueva. Es bien conocido que la simple diferenciación de una señal pude
hacer decrecer por varias veces la relación señal a ruido.
Los observadores de estado para sistemas lineales estiman estados no medibles
sin un sistema de diferenciación. El diseño de un observador de estado lineal es muy
similar al diseño de un controlador por ubicación de polos. En otras palabras, la
expresión de la matriz de ganancia del observador de estados es la expresión dual de
la matriz de ganancia del controlador de realimentación de estados.
Vimos en la subsección 4.3.3 que linealización exacta nos permite diseñar una
ley de control desacoplada (subsección 4.3.4) usando eigenvalores prescritos para
asegurar la estabilización del sistema descrito por (4.29). Apelando al concepto de
dualidad como en el caso lineal discutido lı́neas arriba, nosotros enfrentaremos el
problema de sı́ntesis de observadores no lineales con eigenvalores prescritos para la
estimación de estados. estimation.
4.4 Observadores No Lineales con Polos Prescritos 113
∂Φ
ż = ẋ = O(x)(f (x) + g(x)u) = Mz + Nv
∂x
v = Lnf h(Φ−1 (z)) + Lg Lfn−1 h(Φ−1 (z))u
y = Cz (4.45)
El observador no lineal para el sistema descrito por las ecuaciones 4.22) y (4.25) posee
la forma:
ė = (M − PC)e (4.55)
Como en el caso SISO, la matriz de ganancia P requiere ser seleccionado de modo tal
que los eigenvalores de la ecuación caracterı́stica del observador no lineal: eigenvalues
of the characteristic equation of the nonlinear observer
donde µ1 , . . . , µn son los polos deseados a lazo cerrado del observador no lineal.
1 0 0 0 0 0 P11 P12
0 0 0 1 0 0 P21 P22
0 1 0 0 0 0 P31 P32 x1 − x̂1
0 0 0 0 1 0 P41 P42 x2 − x̂2
0 0 1 0 0 0 P51 P52
0 0 0 0 0 1 P61 P62
Ejemplo 4.5
Solución: La figura 4.1 ilustra el diagrama de bloques usado para la simulación del
sistema de control de trayectoria no lineal MIMO para el MRE. Dicha simulación
emplea los resultados obtenidos en los ejemplos 4.1 a 4.4. Los polos del sistema
a lazo cerrado controlado se ubican en −2,2 ± 0,85 and −6500. Las trayectorias
deseadas seleccionadas son: ρ1 = sin(0,75kT ) y ρ2 = cos(0,75kT ). Los valores de los
parámetros se dan en la tabla 1.6. En adición, considere una masa de la muñeca de la
mano del manipulador de 0.2 kg y una masa de la pinza de of 0.4 kg. La pinza puede
llevar consigo una carga de 0.2 kg of mass. Las condiciones iniciales para la base es
de o rad, mientras que para el brazo (el eslabón) es de π/2 rad.
Las figures 4.1 y 4.2 muestran los resultados de la simulación obtenidos con
el programa srm4fl.m. El tiempo de estabilización y la fuerza de control máxima en
la salida y1 del sistema son 2 s y 2.2 V, mientras que para la segunda salida y2 son
4 s y 26 V. Una segunda simulación sin considerar la carga revela que se obtienen
aproximadamente los mismos resultados..
Linear controller
ρ u x
v = v ( ρ , z) u = u( ^
x, v) x=fx+Gu
Nonlinear process
Pole-placement
loop Nonlinear
observer
^ ^
x
z
z = ψ (x)
^ ^
Linearization
Transformation loop
Figura 4.1: Diagrama de bloques del sistema de control no lineal por realimentación
de estados.
1.5
−0.5
−1
−1.5
0 1 2 3 4 5 6 7 8 9 10
TIME (s)
2.5
CONTROL VOLTAGE u1 (vol)
1.5
0.5
−0.5
−1
0 1 2 3 4 5 6 7 8 9 10
TIME (s)
1.5
LINK POSITION (rad)
0.5
−0.5
−1
−1.5
0 1 2 3 4 5 6 7 8 9 10
TIME (s)
30
CONTROL VOLTAGE u1 (vol)
20
10
−10
0 1 2 3 4 5 6 7 8 9 10
TIME (s)
PROBLEMAS
IN PROGRESS
Capı́tulo 5
Control Deslizante
dxn
x(n) = = f (x) + b(x)u (5.1)
dtn
donde el escalar x es ;la salida de interés del sistema, el escalar u es la entrada de
control, x = [x ẋ x(2) . . . x(n−1) ] es el vector de estado, y las funciones f (x)
y g(x) no son exactamente conocidas pero son acotadas por funciones conocidas
dependientes de x. En adición, b(x) es de signo conocido. Dada una trayectoria de
(2) (n−1)
estado variante con el tiempo xd = [xd ẋd xd . . . xd ], con:
d
s(x, t) = ( + λ)n−1 x̃ = (p + λ)n−1 x̃ = x(n−1) − xr(n−1) (5.4)
dt
donde λ es una constante positiva cuya selección se discutirá luego, p es el operador
(n−1)
de L:aplace y xr es una función que puede ser computada de x y xd . Por ejemplo,
para un sistema de orden n = 3, s toma la forma:
S(t)
(a)
dx
dt s=0
. xd
x
0
(b)
dx
s=0
dt
. xd
x
0
(c)
Figura 5.1: (a) La superficie de deslizamiento S(t). (b) Convergencia exponencial. (c)
Fenómeno del “chattering”.
nos conduce a thit ≤ s(t = 0)/η. Se puede obtener el mismo resultado si se arranca
con s(t = 0) < 0. Por consiguiente:
La figure 5.1(b) ilustra el caso de una trayectoria de estado que evoluciona con una
condición inicial arbitraria, y luego golpea a la superficie s(t) = 0 en un tiempo finito
thit ≤ |s(t = 0)|/η. En el modo de deslizamiento, tal trayectoria “se desliza.a lo largo
de S(t) con el objeto de alcanzar exponencialmente a xd con una constante de tiempo
igual a 1/λ. Por cierto, la expresión de (5.4) para n = 2:
1/λ
s = (p + λ)x̃ → x̃ = s
(1/λ)p + 1
Ejemplo 5.1
sgn(x) = +1 if x > 0
sgn(x) = −1 if x < 0
5.1 Control Deslizante para Sistemas de una Entrada 123
La figura 5.2 muestra las trayectorias del sistema x(t) = −t y x(t) = t originadas
por la ley de control u+ y u− respectivamente. La superficie s(t) = x(t) = 0 es
una superficie de deslizamiento porque satisface la condición (5.6) para η ≤ 1. Tal
superficie será golpeada por alguna trayectoria de estado del sistema en estudio en
un tiempo (5.7).
x
s=0
5.1.2. Control deslizante para Sistemas de la Forma x(n) = fi + u
Considere el sistema de orden n:
dxn
x(n) = = fi + u (5.9)
dtn
i
s = (p + λ)(n−1) x
% = x(n−1) − xr(n−1) (5.11)
Entonces:
ṡ = x(n) − x(n)
r = fi + u − x(n)
r (5.12)
i
− k(x)sgn(s)
u=u (5.14)
Finalmente, sustituyendo:
k= Fi + η (5.16)
i
Ejemplo 5.2
donde a1 = 2,1, a2 (t) y a3 (t) son desconocidas pero verifican que 2,5 < a2 (t) < 5 y
2 < a3 (t) < 4 respectivamente. Por consiguiente,
POSITION x(t)
1
−1
−2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
20
CONTROL INPUT u(t)
10
−10
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0.2
s(t)−TRAJECTORIES
0.1
−0.1
−0.2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
TIME (s)
− k(x)sat(s/Φ)
u=u (5.17)
d
s≥Φ → (s − Φ) ≤ −η
dt
d
s≤Φ → (s − (−Φ)) ≤ η
dt
126 Control Deslizante
^u
−φ φ
s
B(t)
Therefore
1d 2
|s| ≥ Φ → s ≤ (Φ̇ − η)|s| (5.20)
2 dt
Para verificar (5.20), se requiere adicionar la cantidad −Φ̇ a la ganancia de control
discontinua k(x). La ganancia de control modificada k̄(x) es entonces:
− k̄(x)sat(s/Φ)
u=u (5.22)
k̄(xd )
λ= (5.25)
Φ
Reemplazando (5.25) en (5.24) se obtiene:
1
s= (−T (∆fi (xd ), x̃, xd )) (5.26)
p+λ
Sabemos que la variable s es una medida de la distancia algebraica a la superficie de
deslizamiento S(t) y, de acuerdo a (5.26), constituye la salida de un filtro de primer
5.1 Control Deslizante para Sistemas de una Entrada 127
orden filter, cuya dinámica (p + λ) depende de xd (ver (5.25)), y cuyas entradas son
perturbaciones e incertidumbres. Por consiguiente, tal filtro puede librar al sistema
del fenómeno chattering.
La estructura del error dinámico a lazo cerrado cuando Φ = Φ(t) se muestra
en la figura 5.5, donde el primer filtro se diseña de acuerdo a (5.26) para eliminar
perturbaciones e incertidumbres, de modo tal que el segundo filtro pasa bajo, de
acuerdo a la definición (5.4) pueda proporcionar el error de seguimiento x̃. De (5.25)
podemos deducir que el espesor Φ de la capa de frontera requiere ser sintonizada tal
que (5.26) pueda representar un filtro de primer orden con ancho de banda λ. Ahora,
~
T(∆f (x d ),x~ , x d ) 1 s 1 x
(p + λ) (p + λ)
n-1
Figura 5.5: Estructura del error dinámico a lazo cerrado para Φ = Φ(t).
Φ̇ + λΦ = k(xd ) (5.27)
Ejemplo 5.3
POSITION x(t)
1
−1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
20
CONTROL INPUT u(t)
10
−10
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
1
s(t)−TRAJECTORIES
a
0.5
0
b
−0.5
c
−1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
TIME (s)
Control Integral
t
Let 0 x̃(q)dq instead of x̃ be the variable of interest in the surface (5.4). Therefore,
the process (5.9) turns into a nth -order process with respect to the variable of interest,
namely
t
d
s(x, t) = ( + λ) n
x̃(q)dq = (p + λ)n x̃ = x(n−1) − xr(n−1) (5.29)
dt 0
Since equations (5.29) and (5.11) possess the same form, then relations (5.14) and
(n−1)
(5.16) remain unchanged. Note that the term xr in (5.29) contains the variable
of interest.
IN PROGRESS
5.1.3. Control Deslizante para Sistemas de la Forma x(n) = fi + bu
Considere ahora el sistema de orden n:
dxn
x(n) = = fi + bu (5.30)
dtn
i
Recordar que |fi − fi | es una incertidumbre tipo aditivo. Dado que la ganancia de
control b multiplica la entrada de control u, parece entonces razonable definir un
estimado de tipo multiplicativo b:
b = (bmin bmax )1/2 (5.32)
Entonces:
ṡ = x(n) − x(n)
r = fi + bu − x(n)
r (5.35)
i
b el control equivalente que puede lograr que ṡ = 0; por consiguiente:
Sea u
b = b−1 ( −
u fi + x(n) −1
r )= b u (5.36)
i
(n)
=−
donde u fi + xr (ver ecuación (5.13)). Con u definida como:
u = b−1 (
u − k sgn(s)) (5.37)
el lı́mite:
k≥β Fi + η + |β − 1||
u| β = b b−1 (5.38)
i
satisface la condición de deslizamiento (5.6). Para demostrar este hecho, reemplazar
(5.37) en (5.35):
ṡ = fi − b b−1 fi − (1 − b b−1 )x(n) −1
r − b b k sgn(s)
i i
= (fi − fi ) + fi − b b−1 fi − (1 − b b−1 )x(n) −1
r − b b k sgn(s)
i i i
= (fi − fi ) − (1 − β −1
u − β −1 k sign(s)
) (5.39)
i
Ejemplo 5.4
El modelo dinámico del sistema MRAF mostrado en la figura ?? y dado por (??) se
puede reformular con x = θ como:
4
x(3) = fi + bu = −a1 ẍ − a2 ẋ − a3 cos x ẋ − a4 sin x + bu (5.40)
i=1
donde:
Beq Ra nKb Ra Beq meq gL
a1 = + a2 = + a3 =
Jeq nKm La Jeq La nKm La Jeq Jeq
meq gL Ka
a4 = b=
nKm Jeq La Jeq La
IN PROGRESS
u = b−1 (
u − k(x)sat(s/Φ)) (5.41)
donde la función de saturación sat(.) fue definida en (5.18). Con el objeto de verificar
(5.20) para el caso β = 0 en la ganancia de control, la ganancia de control modificada
k̄(x) dada en (5.42) necesita ser reformulada como (ver problema 5.1):
λΦ
Φ̇ > 0 → = k̄(x) − Φ̇/βd
βd
λΦ
Φ̇ < 0 → = k̄(x) − βd Φ̇ (5.45)
βd
5.2 Control Deslizante para Sistemas Multivariables 131
λΦ
k̄(xd ) ≥ → Φ̇ + λΦ = βd (xd )
βd
λΦ λΦ
k̄(xd ) ≤ → Φ̇ + = βd (xd ) (5.46)
βd βd2
λΦ
k̄(x) = (k̄(x) − k̄(xd )) + k̄(xd ) = k(x) − k(xd ) + (5.47)
βd
Ejemplo 5.5
IN PROGRESS
λn ε ≈ βd k(xd )
bandwidthn × tracking precision ≈ parametric uncertainty along xd
cuales se suponen ser funciones del tiempo continuamente diferenciables. Los vectores
de error se definen como:
%(t) = q − qdi
q %̇ = q̇ − q̇di
q(t)
q%i (t) = qi (t) − qdi (t) q%̇i (t) = q̇i (t) − q̇di (t)
o en forma matricial:
%̇ =
s(x, t) = s(q, q̇, t) = Lq̃ + q (5.52)
s1 (x, t) s1 (q, q̇, t) 11 · · · 0 q%1 q%̇1
.. .. .. . . .. .. + ..
. = . = . . . . .
sm (x, t) sm (q, q̇, t) 0 · · · mm q%m q%̇m
donde las constantes positivas ii son los elementos de una matriz diagonal L de
orden m × m. Asumiendo que una fuerza de control diseñada es capaz de confinar
todas las trayectorias que se originan en la intersección de tales superficies y hacerlas
permanecer allı́, entonces en tal situación se cumple, de acuerdo a 5.51, que:
Esta última relación nos indica que q̃i (t) y q̃˙ i (t) deben converger exponencialmente a
cero, esto es:
q̃i (t) = qi (t) − qdi (t) = 0 q̃˙ i (t) = q̇i (t) − q̇di (t) = 0
Por consiguiente: qi (t) = qdi (t) y q̇i (t) = q̇di (t), con lo cual se logra el objetivo de
control.
u = u0 − Usgn(s) (5.53)
5.2 Control Deslizante para Sistemas Multivariables 133
+ +
u1 u1 − u−
1 u1 + u− 1 ··· 0 sign(s1 )
.. 1 .. 1 .. .. .. ..
. = . − . . . .
2 + −
2 + −
um um − um 0 · · · um − um sign(sm )
o en función de sus elementos:
1 + 1 +
u0i = ui + u−
i Ui = ui − u−
i
2 2
La derivada de s (ecuación 5.52) produce:
%̇ + q
ṡ = Lq %̈ = Lq
%̇ + (q̈ − q̈d )
Seleccionando:
Ui ≤ |[u0 + Ps − ueq ]j | + ε ε>0 (5.57)
y reemplazando esta última expresión en (5.56) se obtiene:
m
V̇ ≤ −ε |sj | ε>0
j=1
u− +
i + ε ≤ [ueq − Ps]i ≤ ui − ε (5.59)
ū− +
i + ε ≤ [ūeq − P̄ s]i ≤ ūi − ε (5.61)
eq
ūeq = ueq − u
P̄ = P − P
lo cual verifique (5.59). Si seleccionamos ū− +
i = Ki y ūi = − Ki , entonces de (5.60)
obtenemos:
1 i
u0i = (u+ + u−
i ) = [ueq − Ps] Ui = Ki
2 i
o lo que es lo mismo, en forma matricial:
+
u1 + u− 1
1 ..
u0 = . =u eq − Ps
2
u+m + um
−
+
u1 − u− 1 ··· 0 K1 · · · 0
1 .. .. .. .. . . ..
U = . . . = . . . (5.62)
2 + −
0 · · · um − um 0 · · · Km
d=u
eq − Ps (5.63)
Ki ≥ [ueq − Ps − d]i + ε
Ahora, dado que ueq − Ps = −MLq̃˙ + Pq̇ + d + Mq̇d − Ps, entonces las ganancias
Ki se pueden seleccionar siempre que:
Ejemplo 5.6 Design and Simulation of a Sliding Control for the Process
TRM
The elements of the matrices M and P given by (1.113) may be upper-bounded as
|M11 | ≤ m11 M 11 |M12 | ≤ m12 M 12 |M21 | ≤ m21 M 21
|M22 | ≤ m22 M 22 |P11 | ≤ p11 P 11 |P12 | ≤ |p12 θ̇1 | |P 12 |
|P21 | = 0 P 21 |P22 | ≤ m22 P 22
Using (5.65), the variable Ki is computed from
2
Ki = (M ij |q̈dj − jj q̃ j |
˙ + P ij |q̇j − sj |) + ε
j=1
where q1 is the car position r and q2 is the angular position θ1 of the arm. The
switching control force given by (5.64) is found to be
u1 0 K1 sgn(s1 )
= −
u2 d2 sin θ1 K2 sgn(s2 )
in which d2 was found in (1.113).
The simulation of the sliding control system was performed for the following
initial conditions: r(0) = 0 m, θ1 (0) = −π/2 rad. The desired trajectories are set
to: rd (t) = 1 m and θd1 (t) = π4 cos(2πt) rad. The parameter ε takes on the value of
0.1 while the gripper carries out a payload of 0.2 kg of mass. Figure 5.7 depicts the
simulation results. Such results can be obtained executing the m-file trmsc.m.
Ejemplo 5.7 Design and Simulation of a Sliding Control for the Process
SRM
The elements of matrices M and P given by (1.143) can be upper-bounded as
follows
|M11 | ≤ RT 1 (Jeq1 + 2I1 + H22 ) M 11 |M12 | = 0 M̄12
|M21 | = 0 M̄21 |M22 | ≤ RT 2 (Jeq2 + H22 ) M 22
|P11 | ≤ |RT 1 (Beq1 − H22 q̇2 ) + NT 1 | P 11 |P12 | ≤ |RT 1 H22 q̇1 | P 12
|p21 | ≤ |RT 2 H22 q̇1 | P̄21 |P22 | ≤ RT 2 Beq2 + NT 2 P 22
Using (5.65), the variable Ki can be computed from
2
Ki = (M ij |q̈dj − %̇j |
jj q + P ij |q̇j − sj |) + ε
j=1
The switching control force given by (5.53) (see also (1.143))is found to be
u1 0 K1 sat(s1 /φ)
= −
u2 RT 2 m2 g Lx2 cos q2 ) K2 sat(s2 /phi)
where φ = 0.5. The simulation of the sliding control system was performed for the
following initial conditions: q1 (0) = q2 (0) = 0,1 rad. The desired trajectories are:
qd1 (t) = 0.1745 rad and qd2 (t) = (0,1745/2)t2 rad. The parameter ε takes on the
value of 0.2 while the gripper carries out a payload of 0.2 kg of mass.
Figure 5.8 shows the simulation results. Such results can be obtained executing
the m-file srmsc.m.
136 Control Deslizante
1.5 5
1 0
0.5 −5
0 −10
0 2 4 6 8 10 0 2 4 6 8 10
Entradas
CONTROL VOLTAGE u1
CONTROL VOLTAGE u2
5 4
0 0
−2
−5 −4
0 2 4 6 8 10 0 2 4 6 8 10
5 5
SURFACE s1
SURFACE s2
0 0
−5 −5
−10 −10
0 2 4 6 8 10 0 2 4 6 8 10
TIME (seconds) TIME (seconds)
PROBLEMS
Problema 5.1
1.5 0.18
POSITION q1 (rad)
POSITION q2 (rad)
0.16
1
0.14
0.5
0.12
0 0.1
0 1 2 3 4 0 1 2 3 4
0.02 0.04
Control u1 (V)
Control u2 (V)
0 0.03
−0.02 0.02
−0.04 0.01
−0.06 0
0 1 2 3 4 0 1 2 3 4
0.3 0.05
SURFACE s1
SURFACE s2
0.2 0
0.1 −0.05
0 −0.1
−0.1 −0.15
0 1 2 3 4 0 1 2 3 4
TIME (seconds) TIME (seconds)
[1] H. Goldstain, Mecánica Clásica, segunda edición. Editorial Reverté, S.A., 1996,
ISBN: 84-291-4306-8
[3] Jean-Jacques E. Slotine and Weiping Li, Applied Nonlinear Control, Prentice
Hall, Englewood Cliffs, New Jersey 07632, 1991.
[7] Mohsen Shahinpoor, A Robot Engineering Textbook, Harper & Row, Publishers,
New York, Cambridge, Philadelphia, and others, 1987.
[8] V. I. Utkin, “Variable structure systems with sliding modes,” IEEE Trans. Au-
tomat. Contr., vol. 22, pp. 2-22, Feb. 1993.
[13] Katsuhiko Ogata, Designing Linear Control Systems with MATLAB, Prentice
Hall, Englewood Cliffs, New Jersey 07632, 1944.
[16] Joseph La Salle and Solomon Lefschetz, Stability by Liapunov’s Direct Method
With Applications, Academic Press new York, London, 1961.
[21] Olle I. Elgerd, Control Systems Theory, McGraw-Hill Kogakusha, Ltd., Tokyo,
Auckland, Dusseldorf, and others, 1967.
[24] D. E. Koditschek, Proc. 23rd I.E.E.E. Conf. on Decision and Control, Las Vegas,
p. 733, 1984.
[25] E. Bailey and A. Arapostathis, “Simple sliding mode control scheme applied to
robot manipulators,” Int. J. Control, 1987, vol. 45, No. 4, p. 1197–1209.
[30] Manfred Morari y Evanghelos Zafiriou, Robust Process Control. PTR Prentice
Hall, Englewood Cliffs, New Jersey 07632, 1989.
[32] MathWorks, Inc., MATLAB Reference Guide, Prentice Hall, Englewood Cliffs,
New Jersey, primera edición, 1998.
[33] MathWorks, Inc., SIMULINK User’s Guide, Prentice Hall, Englewood Cliffs,
New Jersey, primera edición, 1996.