MECATRONICA
MECATRONICA
MECATRONICA
Editores:
I
Prefacio
III
IV Prefacio
Esta obra consta de ocho capı́tulos invitados que se presentan con una estruc-
tura que facilita el acercamiento y la comprensión mediante el análisis matemático,
métodos de modelado, ejemplos prácticos y casos de estudio de las áreas antes men-
cionadas.
El libro inicia con el área de Robótica que se presenta en los tres primeros capı́tu-
los. En el capı́tulo 1, Garcı́a-Rodrı́guez propone un control descentralizado continuo
y libre de modelo para el caso articular y cartesiano para un sistema de robots mani-
puladores cooperativos en la tarea de manipulación de un objeto rı́gido. La prueba
de estabilidad asintótica de los errores de seguimiento de fuerza y posición se rea-
liza aplicando la teorı́a de Lyapunov. Los resultados de la simulación demuestran
su desempeño. En el capı́tulo 2, Ordaz-Oliver, et al., hacen una aportación en la
técnica de diseño de estrategias de control de movimiento libre basadas en pasivi-
dad aplicadas a sistemas Euler-Lagrange, particularmente en robots manipuladores.
Lo anterior ilustrado con casos de estudio en donde se incluyen pruebas de estabili-
dad en el sentido Lyapunov, ası́ como la presentación de resultados experimentales
donde se muestra el desempeño de los controladores diseñados. En el capı́tulo 3,
Silva-Ortigoza, et al., se centran en los aspectos del análisis, diseño, construcción
e instrumentación de un robot móvil de tipo diferencial. Para el diseño emplea el
potencial de herramientas de software de desarrollo. En la etapa de construcción
describe paso a paso la evolución de la estructura mecánica con explicaciones e
ilustraciones prácticas. La instrumentación la lleva a cabo mediante el uso de una
tarjeta de adquisición de datos para el control del robot.
A continuación se presentan dos aplicaciones de la mecatrónica, en el capı́tu-
lo 4, Muñoz-Hernández, et al., proponen el modelado de una planta hidroeléctrica
con control predictivo basado en el modelo hı́brido. Para lograr este propósito de-
fine las reglas heurı́sticas correspondientes y los resultados obtenidos son evaluados
haciendo una comparación contra el desempeño del control clásico, demostrando
que el control predictivo es más eficiente. Y en el capı́tulo 5, Hernández-Guzmán,
hace un análisis metodológico que justifica el uso de lazos de corriente en sistemas
electromecánicos, cuando su control es de par o de corriente. Presenta además un
estudio para establecer las condiciones que garanticen la estabilidad asintótica del
sistema, partiendo de los problemas de control de posición de un motor de CD con
escobillas y de robots manipuladores equipados con motores de CD sin escobillas.
Los resultados presentados en simulación demuestran el desempeño en ambos ca-
sos.
La instrumentación es parte intrı́nseca de la mecatrónica, en el capı́tulo 6, Juárez-
Guerra et al., presentan la importancia que tiene la instrumentación virtual aplicada
a procesos industriales o eventos desarrollados con fines didácticos en los labo-
ratorios de enseñanza. Se describe la instrumentación tradicional y la virtual para
exponer en forma clara el concepto de instrumentación virtual, ası́ como los compo-
nentes que integran la base de la misma. Una descripción de las diferentes platafor-
mas de desarrollo de la instrumentación virtual se muestra, para finalmente exponer
el lugar que está adquiriendo la misma en la ingenierı́a. Y en el capı́tulo 7, Saldaña-
González presenta los principios básicos de la visión por computadora, ilustrando
con ejemplos prácticos su relación con los sistemas mecatrónicos, especı́ficamente
Prefacio V
VII
VIII Índice general
Barrientos-Sotelo V. R.
CIDETEC-IPN, Departamento de Posgrado. Área de Mecatrónica. Av. Juan de
Dios Bátiz s/n. Esq. Miguel Othón de Mendizábal, “Unidad Profesional Adolfo
López Mateos”. C.P. 07700, México, D.F., MÉXICO
Cano-Corona A.
Universidad Politécnica de Tlaxcala, Carretera Federal Tlaxcala-Zacatelco,
Km. 9.5, San Lorenzo Axocomanitla, Tlaxcala, México. CP. e-mail: ari-
co 96@yahoo.com
Cetina-Domı́nguez O.
Laboratorio Nacional de Informática Avanzada (LANIA A.C.), Rébsamen 80,
Centro, Xalapa, Veracruz, 91000 MEXICO, e-mail: omarcetina@hotmail.com
Domı́nguez-Ramı́rez Omar A
Universidad Autónoma del Estado de Hidalgo, Instituto de Ciencias Básicas e
Ingenierı́a, Centro deInvestigación en Tecnologı́as de Información y Sistemas,Carr.
Pachuca-Tulancingo, Km. 4.5, México, (Tel: +52 771-7172000 ext. 6734), e-mail:
omar@uaeh.edu.mx
Garcı́a-Rodrı́guez R.
Universidad de Chile, Departamento de Ingenierı́a Eléctrica, Av. Tupper 2007,
Santiago, CHILE, e-mail: rgarciar@ing.uchile.cl
Garcı́a-Sánchez J. R.
CIDETEC-IPN, Departamento de Posgrado. Área de Mecatrónica. Av. Juan de
Dios Bátiz s/n. Esq. Miguel Othón de Mendizábal, “Unidad Profesional Adolfo
López Mateos”. C.P. 07700, México, D.F., MÉXICO
Hernández-Guzmán V. M.
Universidad Autónoma de Querétaro. Facultad de Ingenierı́a. A.P. 3-24, Querétaro,
Qro. 76150, México, e-mail: vmhg@uaq.mx
XIII
XIV Lista de Autores
Hernández-Ocaña B.
Universidad Juárez Autónoma de Tabasco, División Académica de In-
formática y Sistemas, Km. 1 Carr. Cunduacán-Jalpa de Méndez, e-mail:
betania h o@hotmail.com
Jarillo-Silva A.
Universidad Autónoma del Estado de Hidalgo, Instituto de Ciencias Básicas e
Ingenierı́a, Centro deInvestigación en Tecnologı́as de Información y Sistemas, Carr.
Pachuca-Tulancingo, Km. 4.5, México, (Tel: +52 771-7172000 ext. 6734), e-mail:
jarillo@uaeh.edu.mx
Jones D.I.
Centre for Advanced Software Technology, Technium Wales network Ffordd Pen-
lan Parc Menai Bangor Paı́s de Gales, LL574HJ U.K., e-mail: dewi@gwefr.co.uk
Juárez-Guerra E.
Universidad Autónoma de Tlaxcala, Facultad de Ciencias Básicas. Ingenierı́a
y Tecnologı́a, Calzada Apizaquito s/n Km. 1.5, Apizaco, Tlaxcala, México.
CP.90300, e-mail: ejuarez@ingenieria.uatx.mx
Mansoor S. P.
University of Bangor School of Informatics, Dean Street, Bangor Paı́s de Gales,
LL57 1UT U.K., e-mail: s.mansoor@bangor.ac.uk
Marciano-Melchor M.
CIDETEC-IPN, Departamento de Posgrado. Área de Mecatrónica, Av. Juan de
Dios Bátiz s/n. Esq. Miguel Othón de Mendizábal, “Unidad Profesional Adolfo
López Mateos”. C.P. 07700, México, D.F., MÉXICO, e-mail: mmarciano@ipn.mx
Mezura-Montes E.
Laboratorio Nacional de Informática Avanzada (LANIA A.C.), Rébsamen 80,
Centro, Xalapa, Veracruz, 91000 MEXICO, e-mail: emezura@lania.mx
Molina-Vilchis M. A.
CIDETEC-IPN, Departamento de Posgrado. Área de Telemática. Av, Juan de Dios
Bátiz s/n. Esq. Miguel Othón de Mendizábal, “Unidad Profesional Adolfo López
Mateos”. C.P. 07700, México D.F., MÉXICO, e-mail: mamolinav@ipn.mx
Muñoz Hernández G. A.
Instituto Tecnológico de Puebla, Av. Tecnológico 420, Col. Maravillas, Puebla,
Puebla, 72220 MEXICO, e-mail: gmunoz@ieee.org
Ordaz-Oliver J. Patricio,
Universidad Politécnica de Pachuca, Ex Hacienda Sta. Bárbara, Municipio de
Zempoala, Hgo., México, (Tel: +52 771-5477510), e-mail: patricio@upp.edu.mx
Lista de Autores XV
Parra-Vega Vicente
Centro de Investigación y Estudios Avanzados, Grupo de Robótica y Manufactura
Avanzada - Cinvestav Saltillo - Carretera Saltillo-Monterrey, Km 13.5, CP
25000, Ramos Arizpe, Coahuila, México (Tel y Fax +52 844-4390305), e-mail:
vparra@cinvestav.edu.mx.
Portilla-Flores E. A.
CIDETEC-IPN, Departamento de Posgrado. Área de Mecatrónica Av. Juan de Dios
Bátiz s/n. Esq. Miguel Othón de Mendizábal, “Unidad Profesional Adolfo López
Mateos”. C.P. 07700, México, D.F., MÉXICO, e-mail: aportilla@ipn.mx
Saldaña González G.
Universidad Tecnológica de Puebla, Antiguo Camino a la Resurrección
1002-A, Zona Industrial, Puebla, Puebla, 72300 MEXICO, e-mail: grisel-
da.saldana@utpuebla.edu.mx
Sánchez-Olavarrı́a C.
Universidad del Altiplano, Eucalipto No. 1, Col. El Sabinal. Tlaxcala, Tlaxcala.
90102, México, e-mail: cesar ari@hotmail.com
Silva-Ortigoza G.
Benemérita Universidad Autónoma de Puebla, Facultad de Ciencias Fı́sico
Matemáticas. Ap. Postal 1152. C.P. 72001, Puebla, Pue., MÉXICO, e-mail:
gsilva@fcfm.buap.mx
Silva-Ortigoza R.
CIDETEC-IPN, Departamento de Posgrado. Área de Mecatrónica. Av. Juan de
Dios Bátiz s/n. Esq. Miguel Othón de Mendizábal, “Unidad Profesional Adolfo
López Mateos”. C.P. 07700, México D.F., e-mail: rsilvao@ipn.mx
Capı́tulo 1
Esquemas de Control Descentralizados para
Manipuladores Robóticos Cooperativos
R. Garcı́a-Rodrı́guez
1.1. Introducción
Desde su aparición los robots manipuladores han sido utilizados para ejecutar di-
versas tareas como: ensamblaje, inspección, soldadura, transferencia de materiales
entre otras. Cada robot manipulador en la mayorı́a de los casos tiene incorporado
un elemento que le permite agarrar/sostener diversos materiales, accesorios y/o he-
rramentales1 para realizar una tarea especı́fica, a través del cual mantiene contacto
R. Garcı́a-Rodrı́guez
Universidad de Chile, Departamento de Ingenierı́a Eléctrica,
Av. Tupper 2007, Santiago, CHILE
e-mail: rgarciar@ing.uchile.cl
1 two-jaw gripper en inglés
1
2 R. Garcı́a-Rodrı́guez
2 Si los parámetros varı́an con el tiempo se requerirá una mayor demanda de computo
3 Ecuaciones Diferenciales Algebraicas
4 R. Garcı́a-Rodrı́guez
nocimiento de las caracterı́sticas del objeto, de forma que los errores dinámicos
residuales son cancelados por un modo deslizante de segundo orden para garantizar
la convergencia de los errores de seguimiento de fuerza y posición. Por otro lado, el
modo deslizante de segundo orden nos permite obtener un control suave, continuo
y libre de castañeo.
Es importante mencionar que a lo largo de este capı́tulo el objeto será consi-
derado como rı́gido e indeformable. Adicionalmente, supondremos que no existe
deslizamiento ni rotación entre cada robot y el objeto, es decir, no hay un movimien-
to relativo entre ambos.
Este capı́tulo está organizado de la siguiente forma: en la Sección 1.2 se pre-
senta el modelo dinámico del sistema cooperativo de manipuladores robóticos, la
dinámica del error en lazo abierto se presenta en la Sección 1.3. La definición de las
referencias nominales de fuerza y posición -cartesianas y articulares se presentan en
la Sección 1.4, mientras en la Sección 1.5 se presenta el diseño del controlador. Fi-
nalmente en la Sección 1.6 y 1.7 se presentan los resultados de simulación obtenidos
y las conclusiones, respectivamente.
Comenzaremos esta sección con una serie de definiciones que serán de utilidad a
lo largo de este capı́tulo. Posteriormente se presentará el modelo dinámico para un
sistema cooperativo de manipuladores robóticos.
Definición 1. Se llama restricción holonómica aquella ecuación algebraica que
aparece en un sistema a nivel cinemático y es descrita como
ϕ (x,t) = 0 (1.1)
ϕ (q1 , q2 , . . . , ql ) ≡ Σi=1
l
ϕi (qi ) = 0 (1.2)
donde ϕi (qi ) es un vector de funciones que depende solo de las variables articulares
del robot manipulador-i.
Se asumirá en este capı́tulo que el número de restricciones del sistema coopera-
tivo es igual al número de robots manipuladores que lo conforman, es decir, se tiene
una restricción de movimiento para cada robot manipulador del sistema. Además de
que cada robot manipulador sujeta al objeto en un solo punto de contacto.
Hi (qi )q̈i +Ci (qi , q̇i )q̇i + gi (qi ) = τi + JϕTi + (qi )λi , i = 1, . . . , l (1.3)
ϕi (qi ) = 0 (1.4)
donde qi , q̇i , q̈i ∈ ℜn son los vectores de variables de posición, velocidad y acele-
ración articular respectivamente del robot i, Hi (qi ) ∈ ℜnxn representa la matriz de
inercia, Ci (qi , q̇i ) ∈ ℜnxn es la matriz Coriolis la cual modela las fuerzas centrı́petas
y centrı́fugas, gi (qi ) ∈ ℜn es un vector que representa los pares gravitacionales, τi ∈
ℜn representa los torques de control, JϕTi + (qi )λi ∈ ℜn es un vector que representa las
fuerzas de reacción o de contacto ejercidas en el robot manipulador-i por el sistema
cooperativo y λi ∈ ℜ1 es el multiplicador de Lagrange.
Note que las fuerzas de contacto entre los robots manipuladores y el objeto son
mapeadas a la dinámica de cada manipulador vı́a el jacobiano de la restricción
definido como
∂ ϕi (qi ) ∂ ϕi (x(qi ))
Jϕi (qi ) = = ∈ ℜ1xn (1.5)
∂ qi ∂ qi
de tal forma para un sistema cooperativo formado por l manipuladores robóticos
tenemos que
ϕ̇i (qi ) = Σi=1
l
Jϕi (qi )q̇i = 0 i = 1, . . . , l. (1.6)
Con la finalidad de que los multiplicadores de Lagrange puedan ser considerados
como la magnitud de la fuerza de contacto, usaremos el jacobiano de la restricción
JϕTi + (qi ) ∈ ℜ1xn normalizado, es decir
∂ ϕi (q)
JϕTi
JϕTi + (qi ) = ° ∂ ϕ∂ q(qi ° = (1.7)
° i )° Jϕi JϕTi
° ∂ qi °
Dado que la dinámica de cada robot manipulador del sistema cooperativo está su-
jeta a una restricción de la forma (1.4), la dinámica queda restringida a una variedad
de la forma © ª
M = (qi , q̇i ) : ϕi (qi ) = 0, JϕTi (qi )q̇i = 0 (1.8)
es decir, cuando la dinámica de cada robot manipulador evoluciona en su corres-
pondiente variedad (1.8), podemos asegurar que dicho robot manipulador cumple
con su restricción asociada ası́ como con la trayectoria de posición deseada en todo
instante de tiempo, manteniendo el contacto con el objeto manipulado.
Por otro lado, las restricciones definidas en (1.2) modelarán al objeto causando
una reducción de grados de libertad del sistema cooperativo, esto se debe a que los
elementos terminales de los robots manipuladores deben mantener el contacto con
el objeto, por lo que no se pueden mover en todas las direcciones. De tal forma
el sistema de manipuladores robóticos pierde tantos grados de libertad como res-
tricciones independientes de movimiento sean impuestas sobre ellos. Los grados de
libertad perdidos se convierten en grados de libertad de fuerza dado que en las di-
6 R. Garcı́a-Rodrı́guez
Hi (qi )q̈i + (Ci (qi , q̇i ) + Boi )q̇i + gi (qi ) = Yi (qi , q̇i , q̈i )Θi
donde Yi (qi , q̇i , q̈i ) ∈ ℜnxp es una matriz de funciones conocidas, Θi ∈ ℜ p es el vector
de parámetros fı́sicos tales como masas de los eslabones, momentos de inercia y
productos de estos.
Propiedad 4. Las fuerzas de reacción JϕTi + (qi )λi ejercidas sobre el manipulador
i por el sistema cooperativo se encuentran en un plano normal al plano tangente
del punto de contacto i, por lo que el span(JϕTi (qi )) es exactamente esta imagen. El
plano tangente al punto de contacto i es llamado plano de movimiento, y en el se
encuentran las variables de velocidad articular q̇i , por lo que
d ϕi
= Jϕi (qi )q̇i = 0
dt
Propiedad 5. Cada robot manipulador involucrado en un sistema cooperativo
establece un mapeo estrictamente pasivo entre los torques τi y las velocidades de ar-
ticulación q̇i [19], [24]. La propiedad de pasividad es muy importante ya que implica
que el sistema no puede generar energı́a, y en caso de que el sistema sea estricta-
mente pasivo [18], se tiene que el sistema además de no generar energı́a, la disipa y
por lo tanto es asintóticamente estable en lazo abierto.
Dado que los robots manipuladores que forman un sistema cooperativo están su-
jetos a restricciones holonómicas, existen ciertos grados de libertad de posición y
velocidad, en los cuales los robots manipuladores no pueden evolucionar. En es-
1. Control Descentralizado de Manipuladores Robóticos Cooperativos 7
tas direcciones restringidas para cada robot manipulador aparecen los conceptos de
velocidad y posición restringida.
Para lograr un trabajo cooperativo eficiente de cada uno de los robots manipula-
dores, se han definido algunas variables que nos indicarán los efectos de las fuerzas
que cada robot manipulador ejerce sobre los otros a través del objeto, es decir, son
variables que contienen información de los efectos de las fuerzas de contacto en
términos de lo que se ha denominado posición y velocidad restringida [19].
Definición 2. Dada una restricción holonómica se tiene que la variable de posi-
ción y velocidad restringida para el manipulador-i es definida como
Z t
pi = Jϕi (qi )q̇i dt, ṗi = Jϕi (qi )q̇i (1.9)
0
donde ∆ pi = pi − pdi , ∆ ṗi = ṗi − ṗdi . De tal forma, tanto los errores de posición y
velocidad restringida deben ser tomados en cuenta en el diseño del controlador.
donde el regresor YriΘi = Yri (qi , q̇i , q̇ri , q̈ri ) está compuesto de funciones no lineales
conocidas y Θi ∈ ℜ p representa parámetros constantes desconocidos, con (q̇ri , q̈ri )
a ser definidas posteriormente. Si sumamos y restamos (1.12) a (1.3), la ecuación
del sistema en lazo abierto está dada como
JϕTi (qi )q̇i = 0, Qi (qi )q̇i = q̇i Qi (qi )JϕTi (qi ) = 0 (1.15)
De esta manera q̇i puede ser escrita como la suma de un componente tangente y
otro normal, de la forma
q̇i = q̇i + (Jϕi + (qi )Jϕi (qi )q̇i − Jϕi + (qi )Jϕi (qi )q̇i )
= (In − Jϕi + (qi )Jϕi (qi ))q̇i +Jϕi + (qi ) Jϕi (qi )q̇i
| {z } | {z }
Qi ṗi
donde Jϕi + (qi ) es la pseudoinversa de Penrose y Qi (qi ) = In − Jϕi + (qi )Jϕi (qi ) es la
matriz de proyección. Esto nos indica que Qi (qi ) y Jϕi (qi ) son dos espacios ortogo-
nales tal que ℜni puede ser escrito como la suma dada entre rank(im(Qi )) = mi (=
ni − ri ) y rank(im(Jϕi )) = ri tal que mi + ri = ni . Es decir, surge una descomposición
natural de q̇i en el punto de contacto.
1. Control Descentralizado de Manipuladores Robóticos Cooperativos 9
con
donde Svxpi y Svx fi son las variedades ortogonales extendidas de posición y fuerza,
respectivamente.
A partir de la forma general para la señal de referencia nominal qri para el robot
manipulador-i, a continuación se presentarán dos señales de referencia nominal de
posición tanto para el caso cartesiano como para el caso articular. Estas señales de
referencia involucrán un cambio de coordendas, el cual permite que el sistema en
lazo cerrado mantenga la propiedad de pasividad propia del sisietma en lazo abierto,
5Se asume que cada robot manipulador cuenta con un sensor de posición en cada eslabón ası́ como
un sensor de fuerza en su elemento terminal
10 R. Garcı́a-Rodrı́guez
= Q(qi )Svxpi
Xi = fi (qi ) (1.25)
donde Ji (qi ) ∈ ℜnxm es la matriz jacobiana del robot manipulador-i. Dado que el
rank(Ji (qi )) = n tenemos que
donde
12 R. Garcı́a-Rodrı́guez
Sxi = S pi − ρ Sd pi , (1.31)
−k p (t−t0 )
Sd pi = S pi (t0 )exp . (1.32)
con µi > 0 ηi > 0, Kdi = Kdi T ∈ ℜnxn es una matriz diagonal simétrica definida posi-
Hi (qi )Ṡri = −[Kdi +Ci (qi , q̇i )]Sri −YriΘi + JϕTi + [Ṡv fi + ηi Sv fi ] (1.34)
1. Control Descentralizado de Manipuladores Robóticos Cooperativos 13
1 l ¡ T ¢
V= ∑
2 i=1
Sri Hi (qi )Sri + βi SvTfi Sv fi ,
Si Kdi , βi , ηi son suficientemente grandes y los errores iniciales son los suficiente-
mente pequeños, Sri es estable y convergerá en una vecindad del equilibrio trivial de
tal forma
Sri , Sv fi ∈ L∞ (1.37)
Puesto que las trayectorias deseadas son C2 y las ganancias de retroalimentación
son acotadas, tenemos que (q̇ri , q̈r ) ∈ L∞ , lo cual implica que Yri ∈ L∞ . De esta
forma, kSri kδ1 ≤ ε1 con ε1 > 0 es acotado. Dado que la inversa de la matriz de
inercia, la matriz de Coriolis son acotadas y el vector de pares gravitacionales es
acotado tenemos que existe un ε2 > 0 tal que
14 R. Garcı́a-Rodrı́guez
De esta forma concluimos el acotamiento de todas las señales del error en lazo
cerrado y una estabilidad local de Sri y Ṡri .
Parte 2. Convergencia exponencial para los errores de seguimiento de posición
cartesianos
Multiplicando (1.19) por QT y sustituyendo (1.30) tenemos
Z t
Sxi = (QTi (qi )Qi (qi ))−1 Ji QTi (qi )Sri − γ pi sgn(Sxi (ζ ))d ζ (1.39)
t0
donde QTi (qi )Qi (qi ) ∈ R(n−m)×(n−m) tiene columna de rango completo rank(n − m).
Si multiplicamos la derivadas de (1.39) por SxTi obtenemos
d³ T ´
SxTi Ṡxi = −γ pi |Sxi | + SxTi (Qi (qi )Qi (qi ))−1 Ji QTi (qi )Sri
dt¯
¯d³ T ´¯¯
¯
≤ −γ pi |Sxi | + |Sxi | ¯ (Qi (qi )Qi (qi )) Ji Qi (qi )Sri ¯¯
−1 T
dt
≤ −γ pi |Sxi | + |Sxi |ε3
≤ −µ1 |Sxi | (1.40)
donde para una k p elegida suficientemente grande, tenemos que Sd pi ≈ 0 para algún
tiempo pequeño 0 < td ¿ 1 tal que S pi = 0 ∀ t≥ td > 0.
Es importante notar que a pesar de la incertidumbre en los parámetros, la con-
vergencia de los errores de posición no depende de algún parámetro de la dinámica
del robot.
Por consiguiente los errores de seguimiento cartesianos tienen una solución ex-
ponencial hacia la trayectoria deseada, es decir, se establece la convergencia
xi → xdi (1.42)
ẋi → ẋdi . (1.43)
10 −5
[Nm]
[Nm]
5 −10
0 −15
0 2 4 6 0 2 4 6
t [s] t [s]
Control R12 Control R22
10 −2
8 −4
[Nm]
[Nm]
6 −6
4 −8
2 −10
0 2 4 6 0 2 4 6
t [s] t [s]
d
SvTf i Ṡv f i = −γ fi |Sv f i | + SvTf i (Jϕ Sr )
dt¯ i i ¯
¯d ¯
≤ −γ fi |Sv f i | + |Sv f i | ¯¯ (Jϕi Sri )¯¯ (1.45)
dt
≤ −γ fi |Sv f i | + |Sv f i |ε4 (1.46)
≤ −µ2 |Sv f i | (1.47)
14 14
12 12
[Nm]
[Nm]
10 10
8 8
6 6
4 4
0 1 2 3 4 0 1 2 3 4
t [s] t [s]
0 0
[Nm]
−1 −1
0 1 2 3 4 0 1 2 3 4
t [s] t [s]
Figura 1.3 Fuerza real y fuerza deseada (arriba), Error de fuerza (abajo) .
−3 −3
x 10 Error Cartesiano R11 x 10 Error Cartesiano R21
2 6
0 4
[m]
[m]
−2 2
−4 0
−6 −2
0 2 4 6 0 2 4 6
t [s] t [s]
−3 −3
x 10 Error Cartesiano R12 x 10 Error Cartesiano R22
8 2
6 0
[m]
[m]
4 −2
2 −4
0 −6
0 2 4 6 0 2 4 6
t [s] t [s]
Xd VS X−R11 X Vs X−R21
d
0.51 0.505
0.505 0.5
[m]
[m]
0.5 0.495
0.495 0.49
0 2 4 6 0 2 4 6
t [s] t [s]
Yd VS Y−R12 Yd VS Y−R22
0.53 0.53
0.52
0.52
[m]
[m]
0.51
0.51
0.5
0.49 0.5
0 2 4 6 0 2 4 6
t [s] t [s]
1.7. Conclusiones
Referencias
1. Alford C. O., Belyeau S. M., Coordinated Control of two Robot Arms. IEEE Int. Conf. on
Robotics, 468-473, 1984.
2. Arimoto S., Miyazaki F., Kawamura S., Cooperative Motion Control of Multiple Robot Arms
or Fingers. IEEE Int. Conf, on Robotics and Automation, 4, 1407-1412, 1987.
3. Arimoto S., Theory of Nonlinear Mechanical Systems- A Passivity Based Approach and Cir-
cuit Theoretic Approach. Oxford Clarendon Press, 1996.
4. Bonitz R.G., Hsia T.C., Internal force-based impedance control for cooperating manipulators,
IEEE Trans. on Robotics and Automat. 12(1), 78-89, 1996.
5. Buss M., Hashimoto H. , Moore, J.B., Dextrous Hand Grasping Force Optimization. IEEE
Trans. on Robotics and Automat, 12(3), 406-418, 1996.
6. Chiu C.-S., Lian K.-Y., Wu, T.-C., Robust Adaptive Motion/Force Tracking Control Design
for Uncertain Constrained Robot Manipulators. Automatica, 40(12), 2111-21, 2004
7. Damaren, C. J., An Adaptive Controller for Two Cooperating Flexible Manipulators. Journal
of Robotic Systems, 20(1), 15-21, 2003.
8. Deshpande A., Luntz J., Decentralized Control for a Team of Physically Cooperating Robots.
Int. Conf. on Intelligent Robots and Systems, 2, 1757-1762, 2003.
20 R. Garcı́a-Rodrı́guez
9. Fujii, S., Kurono, S. (1975). Coordinated Computer Control of a Pair of Manipulators. Pro-
ceedings of the Fourth World Congress of the International Federation for Theory of Machines
and Mechanisms, 411-417, 1975.
10. Huang, H. P., Chen, R. S., Modeling and Adaptive Coordination Control of Two-robot System.
Journal of Robotic Systems, 9(1), 65-92, 1992.
11. Hu,Y. R., Goldenberg A. A., An Adaptive Approach to Motion and Force Control of Multiple
Coordinated Robot Arms. IEEE Int. Conf. on Robotics and Automation, 2, 1091-1096, 1989.
12. Kawasaki H., T. Shimizu, S. Ito, Adaptive Coordinated Control of Multiple Robot Arms. Proc.
of the 6th IFAC Symposium on Robot Control, 261-266, 2000.
13. Kawasaki H., S. Ito, Rizaudddin B. R., Adaptive Coordinated Control of Multiple Robot Arms.
Proc. of the 6th IFAC Symposium on Robot Control, 261-266, 2000.
14. Kosuge K., Koga M., Furuta K., Nosaki K., Coordinated Motion Control of Robot Arm Based
on Virtual Internal Model. Proc IEEE Int. Conf. on Robotics and Automation, 1097-1102,
1989.
15. Kwon W, Lee BH. New Optimal Force Distribution Scheme of Multiple Cooperating Robots
using Dual Method. Journal of Intelligent Robotic Systems: Theory Appl., 21(4), 301-326,
1998.
16. Lasky, T.A., Hsia T.C., On Force-Tracking Impedance Control of Robot Manipulators. IEEE
Int. Conf. on Robotics and Automation, 274-280, 1991.
17. Morel G., Bidaud P., A Reactive External Force Loop Approach to Control Manipulators in
the Presence of Environmental Disturbances. IEEE Int. Conf. on Robotics and Automation,
1229-1234, 1996.
18. Lewis F.L., Abdallah C.T., Dawson D. M., Control of Robot Manipulators. New York Macmil-
lan, 1993.
19. Liu Y.H., Arimoto S., Parra-Vega V., Kitagaki S., Decentralized Adaptive Control Of Multiple
Manipulators in Cooperations. International Journal of Control, 67(5), 649-673, 1997.
20. Liu Y.H., Arimoto S., Distributively Controlling two Robots Handling an Object in the Task
Space without any Communication. IEEE Trans. on Automatic Control, 41(8), 1193-1198,
1996.
21. Liu Y.H., Arimoto S., Kitagaki S., Parra-Vega V., Model-based Adaptive and Distributed Con-
troller for Holonomic Cooperations of Multiple Manipulators. Int. Journal of Control, 67(5),
649-673, 1997.
22. Liu J-F, Abdel-Malek K., Robust Control of Planar Dual-arm Cooperative Manipulators.
Robotics and Computer Integrated Manufacturing, 16, 109-119, 2000.
23. Parra-Vega V., Rodrı́guez A., Aromoto S., Hirzinger G., High Precision Constrained Grasping
with Cooperative Adaptive Handcontrol. Journal of Intelligent and Robotic Systems 32(39),
235-254, 2001.
24. Rodrı́guez A., Control de Multi-Robots Manipuladores. Tesis de Maestrı́a, Cinvestav, 1997.
25. Woon L.C., Ge S.S., Chen X.Q., Zhang C. Adaptive Neural Network Control of Coordinated
Manipulators. Journal of Robotics Systems, 16(4), 195-211, 1999.
26. Ge S.S., Huang L.,Lee T.H., Model-based and Neural-network-based Adaptive Control of Two
Robotic Arms Manipulating an Object with Relative Motion. International Journal of Systems
Science, 32(1), 9-23, 2001.
27. Scewczyk J., Plumet F., Bidauf P., Planing and Controlling Cooperating Robots through Dis-
tributed Impedance. Journal of Robotics Systems, 19(6), 283-297, 2002.
28. Sun, D., Mills, J. K., Adaptive Synchronized Sontrol for Coordination of Multirobot Assembly
Tasks. IEEE Trans. on Robotics and Automation, 18(4), 498-510, 2002.
29. Tarn T.J., Bejczy, A.K., Yun, X., Design of Dynamic Control of Two Cooperating Robot Arms:
Closed Chain Formulation. IEEE Int. Conf. on Robotics and Automation pp. 7-13, 1987.
30. Yao B., M. Tomizuka, Adaptive Coordinated Control of Multiple Manipulators Handling a
Constrained Object. IEEE Int. Conf. on Robotics and Automation, 624-629, 1993
31. Zuo B.R., Qian W.H., On the Equivalence of Internal and Interaction Forces in Multifingered
Grasping. IEEE Trans. on Robotics and Automat, 15(5), 934-940, 1999.
Capı́tulo 2
Diseño de Estrategias de Control Basadas en
Pasividad para Sistemas Euler-Lagrange
Aplicado en Robots Manipuladores
J. P. Ordaz-Oliver
Universidad Politécnica de Pachuca
Ex Hacienda Sta. Bárbara,
Municipio de Zempoala, Hgo., México, (Tel: +52 771-5477510)
e-mail: patricio@upp.edu.mx
O. A. Domı́nguez-Ramı́rez, A. Jarillo-Silva
Universidad Autónoma del Estado de Hidalgo. Instituto de Ciencias Básicas e Ingenierı́a.
Centro de Investigación en Tecnologı́as de Información y Sistemas
Carr. Pachuca-Tulancingo, Km. 4.5, México, (Tel: +52 771-7172000 ext. 6734)
e-mail: omar@uaeh.edu.mx, jarillo@uaeh.edu.mx
V. Parra-Vega
Centro de Investigación y Estudios Avanzados.
Grupo de Robótica y Manufactura Avanzada-Cinvestav Saltillo.
Carretera Saltillo-Monterrey Km 13.5, C.P. 25000, Ramos Arizpe, Coahuila, México, (Tel y Fax
+52 844-4390305)
e-mail: vparra@cinvestav.edu.mx
21
22 J. P. Ordaz-Oliver et al.
Lista de acrónimos
descritos de la forma:
Considere que el sistema (2.1) donde x son las variables de estado, f es un campo
vectorial que representa las dinámicas del sistema, en la cual están involucradas las
fuerzas centrı́petas y de Coriolis, las fuerzas tribológicas, y fuerzas de gravedad, los
estados del sistema involucran al vector de posiciones y al vector de velocidades,
g(x) es una matriz de fuerzas de control, y u es el vector de entradas de control o
par de torsión (τ ).
h (q,t) = 0 (2.4)
Definición 2.4. Por otra parte, cuando no es posible la reducción de las ecuaciones
de restricción en la forma (2.4), son llamadas no holónomas. En un sistema no
holónomo, las coordenadas generalizadas no dependen una de otra.
30 J. P. Ordaz-Oliver et al.
Por ejemplo, el péndulo sobre el carro (ver Figura 2.4), el cual tiene tres coor-
denadas generalizadas, x, x1 , y1 , sin embargo solo tiene una sola restricción que es
el movimiento del péndulo (restringido a una circunferencia o llamada restricción
holónoma). El péndulo sobre el carro no solo es un ejemplo de restricciones, sino
también es un ejemplo de grados de libertad nGDL ya que contiene 3 coordenadas
generalizadas y una restricción, que por Definición 2.2 el sistema es de 2GDL.
arbitrario BR , un valor r(R) puede ser encontrado tal que se inicie el estado desde el
interior de la bola Br en un tiempo cero, se garantiza que el estado permanecerá den-
tro de la bola BR .
Definición 2.7. El punto de equilibrio x = 0 de (2.1) es llamado atractor si existe un
número η > 0 que tiene la propiedad de que
lı́m x(t) = 0
t→∞
Definición 2.11. Una función escalar continua V (x) es llamada definida positiva lo-
calmente, si V (0) = 0 y sobre una bola BR0
x 6= 0; ⇒ V (x) > 0 ∀x
x 6= 0 ⇒ V (x) < 0
Si se denota con x el estado del sistema (2.1), una función escalar V (x) representa
actualmente una función implı́cita del tiempo t. Asumiendo que V (x) es diferencia-
ble, y se usa la regla de la cadena [20].
dV (x) ∂ v ∂v
V̇ (x) = = ẋ = f (x) (2.6)
dx ∂x ∂x
2. Control Basado en Pasividad 33
Definición 2.13. Si en una bola BR0 , la función V (x) es definida positiva y tiene
derivadas parciales continuas, y si su derivada a lo largo de la trayectoria del sistema
(2.1) es semi definida negativa, entonces V (x) es llamada función de Lyapunov para
el sistema (2.1).
Una función de Lyapunov puede ser dada en una simple interpretación geométrica.
En la Figura 2.6, el valor del punto V (x1 , x2 ) es visto en un punto hacia abajo del
cono [20]. El término Local en un punto de equilibrio es de particular interés ya que
Br = x ∈ Rn : kxk < r ⊂ D,
donde r ∈ (0, ε ].
Sea β ∈ (0, α ] donde
α = máx V (x),
kxk=r
y defina
Ωβ = x ∈ Br : V (x) ≤ β ⊂ Br .
Elija x(0) ∈ Ωβ y sea x(t) cualquier solución de ẋ = f (x). Por hipótesis [D f V ](x(t))
≤ 0 lo cual indica que
V (x(t)) ≤ V (x(0)) ≤ β ,
para todo t ≥ 0. Esto indica que el conjunto Ωβ es invariante bajo f .
Ahora entonces, existe un δ > 0 tal que
Teorema 2.2. Si en una bola BR0 existe una función escalar V (x) con la primer
derivada parcial continual tal que:
Una condición adicional para V es que tiene que satisfacer el no tener un radio
acotado, porque indica que V (x) → ∞ cuando kxk → ∞ (en otras palabras x tiende a
infinito en cualquier dirección). De aquı́ se obtiene el siguiente resultado.
Teorema 2.3. Suponga que existe una función escalar V (x) del estado x, con la
primer derivada parcial continua tal que
V1 = ρ V α
donde ρ es una constante estrictamente positiva y α > 1 es cualquier escalar, no
necesariamente es entero, entonces V1 también es una función de Lyapunov.
Teorema 2.4. Bajo la hipótesis del Teorema 2.1, V̇ (x) < 0 en D − {0} implica que
x = 0 es asintóticamente estable [7].
V (x(t)) → c ≥ 0,
cuando t → ∞.
Rt
V (x(t)) = V (x(0)) + V̇ (x(τ )) d τ
0
= V (x(0)) − γ (t)
para un t grande, esto debe ser negativo, lo cual contradice el hecho de que V ≥ 0.
Por lo tanto c puede ser cero.
La idea básica es encontrada por la prueba original: retomando que r > 0 tal que
Br ⊂ D. Para esto se elige encontrar
β = mı́n V (x)
kxk=r
Ω = {x ∈ Br |kV (x)k ≤ β }
Definición 2.14. Sea V : Rn → R una función tal que V (x) → ∞ cuando kxk → ∞,
entonces V es llamada como radialmente no acotada.
entonces x es inestable.3
{x|kxk = r},
y entonces existe un ε > 0 tal que para algún x(0) arbitrariamente cerrado en el
punto de equilibrio, no se puede permanecer dentro de ε . Por lo tanto el punto de
equilibrio es inestable.
Teorema 2.8. Sea L(x, τ ) un escalar con una función de variables de estado x y τ .
∂ 2 L(x,τ ∗ )
∂ u2
existe y es acotada y continua. También se asume que
∂ L(x, τ ∗ )
=0 (2.10)
∂u
y
∂ 2 L(x, τ ∗ )
> 0. (2.11)
∂ τ2
Entonces u∗ es un mı́nimo local [18].
2.3. Sistemas EL
Modelo dinámico
donde E es la suma total de las energı́as,6 las energı́as cinéticas K (q, q̇) y las ener-
£ ¤T
gı́as potenciales U (q), y q = q1 . . . qn , es una función explı́cita del tiempo que
corresponde a un vector de coordenadas generalizadas.
El Lagrangiano L (q, q̇) de un robot manipulador es la diferencia entre su energı́a
cinética K (q, q̇) y su energı́a potencial U (q) de tal forma que:
1 n 1
K (q, q̇) = ∑ mi v2i = 2 q̇T D(q)q̇
2 i=1
(2.14)
1 ∂ £ T ¤
Ċ(q, q̇)q̇ = Ḋ(q)q̇ − q̇ D(q)q̇ , (2.19)
2 ∂q
∂ U (q)
G(q) = (2.20)
∂q
de hecho (2.26) es una consecuencia de la definición de C(q, q̇). También las de-
sigualdades (2.24) y (2.25) son una definición de los sı́mbolos de Christoffel.
Zt
V (x) −V (x0 ) ≤ yT (s)u(s)ds, (2.29)
0
Una vez que se ha comprendido la teorı́a anterior, ahora se puede realizar una
conexión con la teorı́a de control para el desarrollo de un algoritmo de control como
se presenta en la siguiente sección.
En esté capı́tulo solo se aborda a sistemas completamente actuados, ya que los sis-
temas subactuados son un caso de estudio que particularmente involucra un análisis
más complejo.
Figura 2.7 Dispositivo háptico PHANToM Premium 1.0 empleado como mecanismo robótico
experimental [17].
Del modelado por EL se obtiene que las ecuaciones que describen el compor-
tamiento del sistema son definidas por (2.21), donde la matriz de inercias se define
como D(q) ∈ Rn×n , la matriz de fuerzas centrı́petas y de Coriolis como C(q, q̇)
∈ Rn×n y finalmente el vector de fuerzas de gravedad de la forma G(q) ∈ Rn×n y
sus elementos están determinados de la siguiente forma:
m11 0 0 θ̈1 c11 c12 c13 θ̇1 0 τ1
0 m22 m23 θ̈2 + c21 0 c23 θ̇2 + g2 = τ2 (2.34)
0 m32 m33 θ̈3 c31 c32 0 θ̇3 g3 τ3
g1 = 0
g2 = g{L1(ma + 0.5mc) + L5mbe } cos(θ2 )
g3 = g{0.5L2ma + L3mc − L5md f }sen(θ3 )
Los parámetros aproximados del modelo dinámico que representan masas y lon-
gitudes de los eslabones y elementos del dispositivo experimental son:
Parámetro Valor Unidad
ma 17.5x10−3 Kg
mc 10.4x10−3 Kg
mbe 0.2214 Kg
md f 0.1106 Kg
L1 13.97 cm
L2 13.97 cm
L3 0.325 cm
L4 0.368 cm
L5 0.527 cm
g 9.81 m/s2
Para mayor detalle consultar [2].
0.5
Z [m]
−0.5
−1
0.03
0.025 0.04
0.02 0.035
0.03
0.015 0.025
Y [m] 0.01 0.02
X [m]
donde x representa los estados del sistema (2.1), V (x) ∈ R la función candidata a ser
Lyapunov, por lo tanto su derivada temporal tendrá la siguiente estructura:
dx
V̇ (x) = x = xẋ
dt
entonces se puede cerrar el lazo de la función derivada de Lyapunov con (2.1) y
obtener lo siguiente:
V̇ (x) = x [ f (x) + g(x)u]
con esta última ecuación y propiedades explı́citas del sistema no lineal, como son
sus propiedades dinámicas (en el caso de estudio, robots manipuladores, son las
propiedades dinámicas de un sistema EL), entonces se puede generar una función de
control de la función de Lyapunov con el uso de la entrada u que induzca estabilidad,
esto genera un gran potencial para el diseño de técnicas de control.
El caso de estudio de este capı́tulo prácticamente es basado en robots manipu-
ladores de nGDL, sin embargo la teorı́a es aplicada de forma general a sistemas
dinámicos EL.
8Se propone que las dinámicas no incluyen efectos tribológicos, ya que estos inducen disipación
de energı́a y en consecuencia estabilidad.
50 J. P. Ordaz-Oliver et al.
en el caso de control de posición se tiene que q̃˙ = q̇− q̇d donde q̇d = 0 por ende q̃˙ = q̇
y por propiedades de matrices definidas positivas (2.38) se expresa de la siguiente
forma:
V̇ (x) = q̇T τ − q̇T G(q) + q̇T Kp q̃ = q̇T [τ − G(q) + Kp q̃] (2.39)
ya que (2.39) no tiene caracterı́sticas de ser una función definida positiva, entonces
se propone que la derivada V̇ (x) tenga la estructura de una función de Lyapunov, i.e.
que V̇ (x) < 0, entonces se propone lo siguiente:
Como se puede ver el la Figura 14(a) el error tiende a ser cero, sin embargo no
alcanza el equilibrio, esto debido a que el control es insuficiente para atenuar las no-
52 J. P. Ordaz-Oliver et al.
linealidades que se presentan en el desempeño del mismo, esto puede ser visto com
más claridad en las Figuras 15(b) y 14(b) que no sigue la trayectoria eficientemente.
Por la estructura fı́sica del controlador, se puede concluir que es un controlador
PD con compensación de energı́a potencial no es apto para el seguimiento de trayec-
torias, su interpretación puede ser descrita por la Figura 2.12.
Figura 2.11 Interpretación fı́sica de un sistema de control, robot manipulador con tres grados de
libertad.
Figura 2.12 Interpretación fı́sica de un sistema de control por compensación de energı́a potencial.
V̇ (x) = 12 q̈T D(q)q̃˙ + 12 q̃˙T Ḋ(q)q̃˙ + 12 q̃˙T D(q)q̈ + 21 q̃˙T K p q̃ + 12 q̃T Kp˙˜q
1 ˙T
= q̃˙T D(q) ˙ T
£ q̈ + 2 q̃ Ḋ(q)q̃ + ¤q̃ K1p q̃T
˙
= q̃˙ τ −C (q, q̇) q̃˙ − G(q) + 2 q̃˙ Ḋ(q)q̃˙ + q̃T K p q̃˙
T
Pos. Deseada
Error[rad]
0.25
0
0.2
−10 0.15
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Señal de error (Segundo eslabón) Posición articular (Segundo eslabón)
Pos. Phantom
Posición[rad]
2
Pos. Deseada
Error[rad]
0.065
0.06
0
0.055
0.05
−2
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Pos. Phantom
Señal de error (Tercer eslabón) Posición articular (Tercer eslabón)
10 Pos. Deseada
Posición[rad]
Error[rad]
0.2
0
0.15
−10 0.1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Del control (2.41) se puede determinar que (2.37) y (2.40) son las funciones de
Lyapunov del sistema (2.18), por simplicidad y dado que V (x) > 0 (excepto en el
equilibrio q = qd ) y V̇ (x) < 0, entonces se cumple con el teorema (2.5), por lo tanto
se concluye que la ley de control (2.41) induce estabilidad asintótica global en el
equilibrio deseado qd para el sistema EL (2.18). Por lo tanto mediante esté análisis se
resuelve de una forma muy sencilla el problema de control de posición. Sin embargo
esté desarrollo da pauta para encontrar una infinidad de controladores.
2. Control Basado en Pasividad 55
Señal de control (Primer eslabón)
Trayectoria en el espacio de trabajo (X vs Y)
1
Par[Nm]
0.5 px−py−pz Deseada
X−Y−Z Real
0
0.02
−0.5
0 2 4 6 8 10
Tiempo [seg] 0.015
Señal de control (Segundo eslabón)
Y−PY[m]
0.1
Par[Nm]
0 0.01
−0.1
−0.2 0.005
0 2 4 6 8 10
Tiempo [seg]
Señal de control (Tercer eslabón) 0
0.4 0
0.05
Par[Nm]
Como se puede observar el control (2.41) y (2.44) son controladores del tipo PD
pero en esté caso tienen compensación de gravedad. Sin embargo es bien sabido
por lo que se propone una función candidata a ser Lyapunov la cual involucre direc-
tamente las referencias anteriores . Entonces la función de Lyapunov, que hasta el
momento se ha manejado (2.42), tendrá la siguiente estructura:
1 1
V (x) = q̃˙T D(q)q̃˙ + q̃T K p q̃
2 2
56 J. P. Ordaz-Oliver et al.
sin embargo se sabe que C (q, q̇) q̃˙ = C (q, q̇) q̇ − C (q, q̇) q̇d , lo cual implica que la
matriz de fuerzas centrı́petas y de Coriolis puede ser expresada de la siguiente for-
ma:
C (q, q̇) q̇ = C (q, q̇) q̃˙ +C (q, q̇) q̇d
por lo que (2.45) puede ser reescrita de la siguiente forma:
£ ¤ 1
V̇ (x) = q̃˙T τ −C (q, q̇) q̃˙ −C (q, q̇) q̇d − G(q) − q̃˙T D(q)q̈d + q̃˙T Ḋ(q)q̃˙ + q̃T K p q̃˙
2
está ultima ecuación da pauta para un diseño más especifico, como se presenta a
continuación:
V̇ (x) = q̃˙T τ − q̃˙T C (q, q̇) q̃˙ − q̃˙T C (q, q̇) q̇d − q̃˙T G(q) − q̃˙T D(q)q̈d
+ 21 q̃˙T Ḋ(q)q̃˙ + q̃T K p q̃˙
= q̃˙T τ· − q̃˙T C (q, q̇) q̇d −¸q̃˙T G(q) − q̃˙T D(q)q̈d
1 (2.46)
+ q̃˙T Ḋ(q) −C (q, q̇) q̃˙ +q̃˙T K p q̃
2
| {z }
propiedad de anti−simetria
= q̃˙T [τ − D(q)q̈d −C (q, q̇) q̇d − G(q) + K p q̃]
como en los controladores anteriores, está función candidata a ser Lyapunov no tiene
una estructura especı́fica, por lo que la aproximamos a la función deseada (2.40), y
obviamente se reduce a lo siguiente:
V̇ (x) = −q̃˙T Kd q̃˙ = q̃˙T [τ − D(q)q̈d −C (q, q̇) q̇d − G(q) + K p q̃]
la entrada de control o el par que hace que se cumpla la ecuación anterior es:
donde esté control es mejor conocido como controlador par calculado PD (ver Figu-
ra 2.16). La ley de control (2.47) representada en la Figura 2.16 al aplicarse sobre el
2. Control Basado en Pasividad 57
sistema cuyas dinámicas son dadas por (2.34), con K p = 0.5, y Kd = 0.4, genera los
siguientes resultados experimentales.
Error[rad]
0.25 1
0.2 0
0.15 −1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Posición articular (Segundo eslabón) Señal de error (Segundo eslabón)
Pos. Phantom
Posición[rad]
Error[rad]
0.06 0
−0.5
0.05 −1
−1.5
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Pos. Phantom
Posición articular (Tercer eslabón) Señal de error (Tercer eslabón)
0.25 Pos. Deseada
Posición[rad]
Error[rad]
0.2 1
0.15 0
0.1 −1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
px−py−pz Deseada
−0.2 X−Y−Z Real
−0.4 0.02
0 2 4 6 8 10
Tiempo [seg] 0.015
Señal de control (Segundo eslabón)
Y−PY[m]
0.3
Par[Nm]
0.2 0.01
0.1
0
−0.1
0.005
0 2 4 6 8 10
Tiempo [seg]
Señal de control (Tercer eslabón) 0
0
0.2 0.05
Par[Nm]
−0.01 0.04
0 0.03
−0.02 0.02
−0.2
0.01
Z−PZ[m] −0.03 0
0 2 4 6 8 10 X−PX[m]
Tiempo [seg]
Comentario 2.2. Como se puede ver en la configuración de estos controles, con una
simple modificación de estos se puede llegar al controlador propuesto por Slotine-
Li,9 es decir, introducir un error dinámico S ∈ Rn×1 de la forma S = q̇ − q̇r donde
q̇r = q̇d − α q̃, q̃ = q − qd , entonces:
S = q̇ − q̇r ,
q̇r = q̇d − α q̃,
q̃ = q − qd , ∴ S = q̇ − q̇d + α q̃
Por lo tanto los controladores (2.44) y (2.47) puede reescribirse de la siguiente for-
ma:
τ = G(q) − Kd S
(2.48)
= D(q)q̈d +C (q, q̇) q̇d + G(q) − Kd S
Ya que la segunda ecuación de (2.48) presenta cierta estructura, entonces se puede
hacer una factorización especial en la que los parámetros estén desacoplados de las
dinámicas (como lo sugiere la propiedad 2.3), por lo que la siguiente factorización
es viable:
Yd (q, q̇, q̇d , q̈d ) Φ = D(q)q̈d +C (q, q̇) q̇d + G(q). (2.49)
dejando el control de la forma:
Pos. Deseada
Error[rad]
0.2 0
0 −0.1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Posición articular (Segundo eslabón) Señal de error (Segundo eslabón)
Pos. Phantom
Posición[rad]
0.2 0.1
Pos. Deseada
Error[rad]
0.1 0
0 −0.1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Pos. Phantom
Posición articular (Tercer eslabón) Señal de error (Trecer eslabón)
0.4 Pos. Deseada 0.2
Posición[rad]
Error[rad]
0.2 0
0 −0.2
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
9I.e. simplemente con cambiar la matriz de ganancias α ∈ Rn×n con el producto matricial Kd−1 K p ∈
Rn×n sobre los controladores (2.44) y (2.47)
2. Control Basado en Pasividad 59
Señal de control (Primer eslabón)
Trayectoria en el espacio de trabajo (X vs Y)
10
Par[Nm]
X−Y−Z Real
0 px−py−pz Deseada
0.03
−10
0 2 4 6 8 10
Tiempo [seg] 0.02
Señal de control (Segundo eslabón)
10
Y−PY[m]
Par[Nm]
0.01
0
−10 0
0 2 4 6 8 10
Tiempo [seg]
Señal de control (Tercer eslabón) −0.01
4
0.02
0.06
Par[Nm]
2 0
0.04
0 −0.02
0.02
−0.04
−2 Z−PZ[m] 0
0 2 4 6 8 10 X−PX[m]
Tiempo [seg]
para dar una explicación más apropiada a estas funciones, las mismas están repre-
sentadas por la Figura 2.21.
f(x)=sinh(x) f(x)=cosh(x)
100 80
50 60
40
f(x)
f(x)
0
20
−50
0
−100
−5 0 5 −5 0 5
x x
f(x)=tanh(x) f(x)=sech(x)
2 1.5
1 1
f(x)
f(x)
0 0.5
−1 0
−2 −0.5
−5 0 5 −5 0 5
x x
vechadas sobre el diseño, en este capı́tulo solo se hará mención de algunas de las
propiedades más importantes:
tanh(x) = 0, ⇔ x = 0.
entonces los controladores que hacen que se cumplan las dos ecuaciones anteriores
serán los siguientes:
Por el tipo de dinámica que se presenta dentro de la ley de control (2.54) y (2.55),
se dice que la aproximación de estos es relativa a la de controladores de estructura
variable o mejor conocidos como Modos Deslizantes de primer orden.
Sin embargo los controladores desarrollados tienen la ventaja que no tienen in-
determinación en el origen, dando ası́ un desarrollo más suave sobre el efecto de
chattering que se presenta dentro de este tipo de controladores, ya que en este caso
β = 1 a diferencia de la aproximación de β → ∞ que se hace comúnmente para
suplir la función signo.
Para la realización experimental, se propone aplicar las formas del controlador
tipo deslizante, i.e. el controlador que contiene la dinámica completa más la estruc-
tura variable (2.55), y el controlador deslizante simple (2.54), proponiendo ganan-
cias KL = 1.1diag(3) y Kd = diag(3), se obtiene lo siguiente.
En la Figura 23(a) se puede ver como el error del controlador simple de estructura
variable se aproxima a cero, i.e. las dinámicas del controlador por tipo deslizante
siguen la trayectoria propuesta de forma eficiente 23(b), sin embargo se produce
el efecto de chattering o cacheteo que es propio de los controladores de estructura
variable, haciendo que el control oscile en cierto intervalo, este efecto puede ser
visto en las Figuras 22(b), 22(a), 23(a).
Como se puede ver en la Figura 25(a) se observa el efecto del controlador de
estructura variable con compensación de la dinámica completa, ası́ como en los
resultados obtenidos del controlador deslizante simple, en la Figura 24(b), se puede
2. Control Basado en Pasividad 63
Posición articular (Primer eslabón) Señal de error (Primer eslabón)
Pos. Phantom 2
Posición[rad]
0.25 Pos. Deseada
Error[rad]
0
0.2
0.15 −2
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Posición articular (Segundo eslabón) Señal de error (Segundo eslabón)
Posición[rad] Pos. Phantom
0.07
Pos. Deseada
Error[rad]
2
0.06 1
0
0.05 −1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Pos. Phantom
Posición articular (Tercer eslabón) Señal de error (Tercer eslabón)
Pos. Deseada
2
Posición[rad]
0.2
Error[rad]
0.15 0
0.1 −2
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Figura 2.22 Estados del robot ante control (deslizante primer orden simple).
px−py−pz Deseada
0.2 X−Y−Z Real
0 0.02
0 2 4 6 8 10
Tiempo [seg] 0.015
Señal de control (Segundo eslabón)
Y−PY[m]
0.2
Par[Nm]
0 0.01
−0.2
−0.4
−0.6
0.005
0 2 4 6 8 10
Tiempo [seg]
Señal de control (Tercer eslabón) 0
0
0.5 0.05
Par[Nm]
−0.01 0.04
0 0.03
−0.02 0.02
−0.5 0.01
Z−PZ[m] −0.03 0
0 2 4 6 8 10 X−PX[m]
Tiempo [seg]
ver el efecto de chattering o cacheteo que es propio de este tipo de controles, y las
Figuras 25(b), y 24(a) se muestra como el la trayectoria descrita por el manipulador
se aproxima a la trayectoria deseada, sin embargo el comportamiento es mucho
mejor que su antecesor.
Posición[rad]
0.25 Pos. Deseada 5
Error[rad]
0.2 0
0.15 −5
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Posición articular (Segundo eslabón) Señal de error (Segundo eslabón)
Posición[rad] Pos. Phantom
Pos. Deseada 1
Error[rad]
0.06 0.5
0
−0.5
0.04
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Pos. Phantom
Posición articular (Tercer eslabón) Señal de error (Tercer eslabón)
Pos. Deseada
0.25 4
Posición[rad]
Error[rad]
0.2 2
0.15 0
−2
0.1
−4
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Figura 2.24 Estados del robot ante control (deslizante primer orden).
px−py−pz Deseada
0
X−Y−Z Real
−1
0.02
0 2 4 6 8 10
Tiempo [seg] 0.015
Señal de control (Segundo eslabón)
0.2
Y−PY[m]
Par[Nm]
0 0.01
−0.2
−0.4 0.005
0 2 4 6 8 10
Tiempo [seg]
Señal de control (Tercer eslabón) 0
0
0.5 0.05
Par[Nm]
−0.01 0.04
0 0.03
−0.02 0.02
−0.5
0.01
Z−PZ[m] −0.03 0
0 2 4 6 8 10 X−PX[m]
Tiempo [seg]
do:
1 1 1 1
V̇ (x) = ṠT D(q)S + ST D(q)Ṡ + ST Ḋ(q)S = ST D(q)Ṡ + ST Ḋ(q)S (2.57)
2 2 2 2
de una extensión para el seguimiento de trayectorias, Slotine-Li maneja la siguiente
propiedad en pasividad de la dinámica del error:
donde Kd (q, q̇) = Kd (q, q̇)T > 0 es una matriz de inyección de amortiguamiento.
Basados en está propiedad y la propiedad de anti-simetrı́a se tiene el siguiente lema:
donde D(q) y KD son matrices definidas positivas, C(q, q̇) satisface (2.58) y define
una salida estrictamente pasiva ∑d : Ψ 7→ S. Consecuentemente, si Ψ ≡ 0 se tiene
S ∈ L2 .
Con el lema anterior y las propiedades de la dinámica del error se puede obtener un
control basado en pasividad y estabilidad en el sentido de Lyapunov, entonces:
donde q̇r y q̈r son referencias nominales. Ası́, introduciendo (2.61), (2.60) en (2.57)
se obtiene el siguiente resultado:
1
V̇ (x) = ST [τ −Yr (q, q̇, q̇r , q̈r )Φ −C(q, q̇)S + Kd S] − ST Ḋ(q)S (2.62)
2
1
V̇ (x) = ST [τ −Yr (q, q̇, q̇r , q̈r )Φ − Kd S] −ST C(q, q̇)S + ST Ḋ(q)S (2.63)
| {z 2 }
propiedad de anti−simetria
donde la matriz de ganancias de control KD por definición tendrá que ser mayor que
la matriz Kd .
Ahora bien, si se tiene una configuración más completa del cambio de variable que
hace Slotine-Li, entonces se puede aproximar un controlador tipo modos deslizastes
de la siguiente forma. Por lo tanto si tomamos la función candidata a ser Lyapunov,
y la aproximamos a una función de la siguiente forma:
Posición[rad]
Pos. Deseada 1.5
0.25
Error[rad]
1
0.5
0.2 0
−0.5
0.15
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Posición articular (Segundo eslabón) Señal de error extendido (Segundo eslabón)
Posición[rad] Pos. Phantom
1
Pos. Deseada
Error[rad]
0.06 0
−1
0.04
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Pos. Phantom
Posición articular (Tercer eslabón) Señal de error extendido (Tercer eslabón)
Pos. Deseada 1
Posición[rad]
Error[rad]
0.2 0.5
0.15 0
−0.5
0.1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Figura 2.26 Estados del robot ante control (deslizante primer orden).
0 2 4 6 8 10
Tiempo [seg] 0.015
Señal de control (Segundo eslabón)
Y−PY[m]
0.4
Par[Nm]
0.01
0.2
0
−0.2
0.005
0 2 4 6 8 10
Tiempo [seg]
Señal de control (Tercer eslabón) 0
0
0.1
0.05
Par[Nm]
−0.01 0.04
0
0.03
−0.02 0.02
−0.1
0.01
Z−PZ[m] −0.03 0
0 2 4 6 8 10 X−PX[m]
Tiempo [seg]
Posición[rad]
Pos. Deseada 2
Error[rad]
0.2 0
−2
0.15
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Posición articular (Segundo eslabón) Señal de error extendido (Segundo eslabón)
Posición[rad] Pos. Phantom
Pos. Deseada
Error[rad]
0.5
0.06 0
−0.5
0.04 −1
−1.5
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Pos. Phantom
Posición articular (Tercer eslabón) Señal de error extendido (Tercer eslabón)
0.25 Pos. Deseada
2
Posición[rad]
Error[rad]
0.2 0
0.15 −2
0.1 −4
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo[seg] Tiempo[seg]
Figura 2.28 Estados del robot ante control (deslizante primer orden).
px−py−pz Deseada
1
X−Y−Z Real
0
0.02
0 2 4 6 8 10
Tiempo [seg] 0.015
Señal de control (Segundo eslabón)
Y−PY[m]
0.6
Par[Nm]
0.4 0.01
0.2
0
0.005
0 2 4 6 8 10
Tiempo [seg]
Señal de control (Tercer eslabón) 0
0
1 0.05
Par[Nm]
−0.01 0.04
0.5 0.03
−0.02 0.02
0 0.01
Z−PZ[m] −0.03 0
0 2 4 6 8 10 X−PX[m]
Tiempo [seg]
como se puede ver en (2.74), para seguir con el algoritmo de diseño que se ha
manejado hasta el momento es necesario que en la parte derecha de la ecuación (el
término que involucra Φ̂˙ (t)Θ −1 Φ̂ (t)) se encuentre el error extendido S, entonces
se requiere que la variación del parámetro de adaptación Φ̂˙ (t) tenga una estructura
especifica la cual involucre el error extendido y como se desea que la estructura de
la función de Lyapunov se mantenga de la forma (2.65), entonces se debe cumplir
lo siguiente:
Φ̂˙ (t) = −Θ YrT S (2.75)
por lo tanto, la función (2.74) se reduce a la siguiente expresión:
£ ¤
V̇ (x) = ST τ −Yr (q, q̇, q̇r , q̈r )Φ̂ − Kd S − Θ YrT SΘ −1 Φ̂ (t)
por lo tanto, para que se cumpla con la condición de que la derivada de la función de
Lyapunov sea semi-definida negativa, es necesario que la ley de control en lazo ce-
rrado, elimine las dinámicas no deseadas, y para obtener el mismo resultado que en
la sección 2.6.3, i.e que V̇ (x) = −ST KD S, la ley de control se escribe de la siguiente
forma:
τ = Yr (q, q̇, q̇r , q̈r )Φ̂ (t) − KL S (2.76)
como se puede observar, la ley de control (2.76) en lazo cerrado indica estabilidad
asintótica.
que al igual que la estructura del controlador (2.76), estos también presentan la
propiedad de estabilidad asintótica.
V̇ (x) → 0 ⇒ S → 0
Ya que V (x) es definida positiva, entonces el lema de Barbalat [7] indica que V̈ (x)
tiende a cero si esta es uniformemente continua, en particular esto se cumple ya
que:
dV̇ (x) T T
dt = −S KD Ṡ − Ṡ KD S
T
= −2S Kd Ṡ
lo cual involucra que se cumple con lo siguiente:
entonces S → 0
como V̇ (x) y V (x) son acotadas, esto implica que S es acotada. Después, esto
implica que q, q̇ son todas acotadas. Ası́ S → 0 y t → ∞, implican que q̃, y q̇r → 0
cuando t → ∞, por lo tanto
2. Control Basado en Pasividad 71
S = q̃˙ + α q̃
Entonces:
q̃˙
q̃˙ + α q̃ = 0 ⇒ q̃˙ = −α q̃, = −α
q̃
Al resolver está ecuación diferencial se obtiene:
Z Z
d q̃
= −α dt, ln q̃ = −α t + c ⇒ q̃ = ce−α t
q̃
entonces, cuando t = 0, q̃ = C y t → ∞, q̃ → 0 ⇒ s → 0.
2. Si V (x) = 0, entonces se tiene que:
1
V (x) = ST D(q)S = 0. (2.78)
2
para todo ε > 0, kzk < r, β ∈ (0, ζ ], donde r ∈ (0, ε ], ζ = maxV (x)kzk=r . Entonces
el segundo método
£ de
¤ Lyapunov [7], para sistemas estables, y por hipótesis se
considera que D f V (z(t)) ≤ 0 (ver (2.56)), esto implica
V (z(t)) ≤ V (z(0)) ≤ β ,
Con esto se puede concluir que la ley de control descrita por (2.66) es asintótica-
mente estable ya que la la función de error sobre la función de Lyapunov indica que
lı́m q̃(t) = lı́m [q(t) − qd (t)] = 0
t→∞ t→∞
como se puede ver en (2.81) y si se aplica el lema 2.1, se obtiene el siguiente resul-
tado:
D(q)Ṡ = Ψ − [C(q, q̇) + Kd ] S (2.82)
donde Ψ son dinámicas del sistema:
Ψ = τ − (D(q)q̈r +C(q, q̇)q̇r + g(q))
(2.83)
= τ −Yr (q, q̇, q̇r , q̈r )Φ = τ −Yr Φ
donde q̇r y q̈r son referencias nominales. Entonces al aplicar el principio de progra-
mación dinámica (2.8), se requiere derivar (2.79), y se obtiene el siguiente resultado:
· ¸ · ¸ · ¸
dV (x̃) 1 T D(q) k 1 D(q) k 1 Ḋ(q) 0
= x̃ x̃˙ + x̃˙T x̃ + x̃T x̃
dt 2 k D(q) 2 k D(q) 2 0 Ḋ(q)
es decir: · ¸ · ¸
dV (x̃) D(q) k 1 Ḋ(q) 0
= x̃T x̃˙ + x̃T x̃
dt k D(q) 2 0 Ḋ(q)
introduciendo (2.80) y (2.81) en la ecuación anterior, se obtiene:
· T ¸T · ¸· ¸ · ¸ · ¸· ¸
dV (x̃) S D(q) k Ṡ 1 ST T Ḋ(q) 0 S
= ˙T +
dt q̃ k D(q) q̃¨ 2 q̃˙T 0 Ḋ(q) q̃˙
es decir:
dV (x̃) 1 1
= ST D(q)Ṡ + q̃˙T kṠ + ST kq̃¨ + q̃˙T D(q)q̃¨ + ST Ḋ(q)S + q̃˙T Ḋ(q)q̃˙ (2.84)
dt 2 2
y por propiedades de sistemas EL y de pasividad en la dinámica de error S, se obtiene
lo siguiente:
2. Control Basado en Pasividad 73
es decir:
£ ¡1 ¢ ¤
V̇ (x̃) = ST£ τ −Y ¡ r Φ + 2 Ḋ(q) −C(q,
¢ q̇) S ¤− Kd S
+q̃˙T τ + 12 Ḋ(q) −C(q, q̇) q̃˙ − G(q)
+q̃˙T kD−1 (q) (τ −Yr Φ −C(q, q̇)S − Kd S)
(2.85)
+ST kD−1 (q) (τ −C(q, q̇)q̇ − G(q))
dV (x̃)
dt = q̃ [τ − G(q)] + q̃ kD (q) (τ −Yr Φ −C(q, q̇)s − Kd S)
˙T ˙T −1
+S [τ −Yr Φ − Kd S] + ST kD−1 (q) (τ −C(q, q̇)q̇ − G(q))
T
donde:
Q1 = R−1 + R−1 kD−1 (q) + kD−1 (q)R−1 kD−1 (q)
Q12 = Q1
Q2 = Q1 + Kd
como det(Q(x)) = Kd , entonces la condición suficiente para que se cumpla Q(x) > 0
es Kd > 0, la norma de la función Yr Φ en (2.88) es acotada superiormente, y con
algún escalar positivo de la forma βi (i = 0, 1..., 5) tal que [6]:
Yr Φ ≤ kD(q)k
° °kq̈r k + k{Kd +C(q, q̇)}¡ q̇°r k°+ kG(q)k ¢
≤ β1¡α °q̃˙° + λ¢M (Kd ) + β2 kq̇k α °q̃˙° + β4 + γ kσ k + β̄3
≤ η q̃, q̃˙, σ , βi
¡ ¢
donde β̄3 = β1 β5 + β3 , y η q̃, q̃˙, σ , βi es un escalar. Acorde con la ecuación ante-
rior y (2.89) se obtiene:
10 ∂2
£ dV ¤
i.e. Tomando ∂ 2τ dt + f0 esto es cero.
2. Control Basado en Pasividad 75
¡ ¢ ° ° ° °
ϕ (x) ≤ ksk η q̃, q̃˙, σ , βi +°°q̃°˙° β3 + °q̃˙¡° ksk λM kβ¢0 ° °
+kqk ksk¡λM kβ0 β2 + °q̃˙°¢λM°kβ°0 η q̃, ° ˙°
¡ q̃, σ , βi +ksk¢λM kβ0 β3 + q̃ ksk λM kβ2
˙
˙ ° ˙ °
≤ ksk η2 q̃, q̃, σ , βi , λM k + q̃ η3 q̃, q̃, σ , βi , λM k
˙
se sabe que:
¡ ¢ ° ° ¡ ¢
ksk η2 q̃, q̃˙, σ , βi , λM k > °q̃˙° η3 q̃, q̃˙, σ , βi , λM k
entonces ¡ ¢
ϕ (x) ≤ ksk η2 q̃, q̃˙, σ , βi , λM k (2.90)
y se tiene que βi > 0 y λM Kd > 0, esto implica:
¡ ¢
V̇ (x) = − ksk2 λM Kd + ksk η2 q̃, q̃˙, σ , βi , λM k (2.91)
¡ ¢
entonces Kd es mayor que η2 q̃, q̃˙, σ , βi , λM k , y esto es una condición suficiente
para concluir que el sistema está en un conjunto compacto el cual involucra los
argumentos de Lyapunov, i.e. involucra V (x) y a su derivada V̇ (x), λM k < λm D(q)
tal que x̃ converge dentro de una vecindad ε > 0 con radio r > 0.
2.8. Conclusiones
Error[rad]
Par[Nm]
0
0 −2
−4
−1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo [seg] Tiempo[seg]
Señal de control (Segundo eslabón) Señal de error extendido (Segundo eslabón)
0.1
Error[rad]
0.5
Par[Nm]
0
−0.1 0
−0.2
−0.3 −0.5
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo [seg] Tiempo[seg]
Señal de control (Tercer eslabón) Señal de error extendido (Tercer eslabón)
2
0.2
Error[rad]
Par[Nm]
1
0
−0.2 0
−0.4 −1
0 2 4 6 8 10 0 2 4 6 8 10
Tiempo [seg] Tiempo[seg]
Posición[rad]
Pos. Deseada
px−py−pz Deseada 0.25
X−Y−Z Real
0.2
0.02
0.15
0 2 4 6 8 10
0.015 Tiempo[seg]
Posición articular (Segundo eslabón)
Pos. Phantom
Posición[rad]
Y−PY[m]
0.05 0.2
−0.01 0.04
0.03 0.15
−0.02 0.02
0.01
−0.03 0 0.1
Z−PZ[m] X−PX[m] 0 2 4 6 8 10
Tiempo[seg]
Referencias
1. Bayard D., Wen J.T., “A New Class of Stabilizing Control Laws for Robotic Manipulators,
Part II: Adaptive Case”, International Journal of Control, 47(5), 1988, pp.1387-1406.
2. Domı́nguez-Ramı́rez Omar A., and Parra-Vega Vicente, “Active Haptic Interface for Remote
Training Purposes”, 11T H International Conference on Advanced Robotics, Coimbra, Portu-
gal, June 30 - July 3, 2003.
3. Hans P. Geering, Optimal Control with Engineering Applications, Springer Berlin Heidelberg
New York, 2007.
4. Jean-Jacques E. Slotine and Weiping Li, Applied Nonlinear Control, Prentice-Hall, Engle-
wood Cliffs, New Jersey, 1991.
5. Kelly R., “A tuning procedure for stable PID control of robot manipulators”, Robotica, vol.
13, pp. 141 - 148, 1995.
6. Kelly Rafael, Santibáñez Victor, Control de Movimiento de Robots Manipuladores, Prentice
Hall, 2003.
7. Khalil, H., Nonlinear Systems, Segunda edición, Prentice Hall, NJ, 1996.
8. Ordaz-Oliver J. Patricio, Domı́nguez-Ramı́rez Omar A., and Muñoz-Palacios Filiberto, “Sub-
Optimal Control Based on Passivity for Euler-Lagrange Sysems”, International Conference
on Control 2008 (UKACC), Manchester UK, 4-6 September, 2008.
9. Ordaz-Oliver J. Patricio, Domı́nguez-Ramı́rez Omar A., and Espinoza-Quesada Eduardo S.,
“Trajectory Tracking Control for Robotics Manipulators Based on Passivity”, Electronics,
Robotics, and Automotive Mechanics Conference, September 25-28, 2008.
10. Ordaz-Oliver J. Patricio, Domı́nguez-Ramı́rez Omar A. “Global Regulation for Robotics Ma-
nipulators Based on Kinetic Energy”, Electronics, Robotics, and Automotive Mechanics Con-
ference, Cuernavaca, Morelos, México, September 25-28, 2007.
11. Ordaz-Oliver J. Patricio, Domı́nguez-Ramı́rez Omar A., “Novel Sliding-mode Control for
Euler-Lagrange Systems: Nonadaptive and Adaptive Case.”, 4th IEEE Latinoamerican
Robotics Symposium, Monterey, Nuevo León, México, Nobemver 9-14, 2007.
12. Ortega Romeo, Lorı́a Antonio, Nicklasson Per Johan, and Sira-Ramı́rez Hebertt, Passivity-
based Control of Euler-Lagrange Systems, Mechanical, Electrical and Electromechanical Ap-
plications, Spriger, 1998.
13. Parra-Vega V., Arimoto S., Yun-Hui Liu, Hirzinger G., Akella P., “Dynamic sliding PID con-
trol for tracking of robot manipulators: theory and experiments”, Robotics and Automation,
IEEE Transactions, Volume 19, Issue 6, pp. 967 - 976, 2003.
14. Parra V., and Arimoto S., “Nonlinear PID control with Sliding modes for Tracking of Robots
Manipulators”, IEEE Conference on Control Aplications, Mexico 2001.
15. Parra-Vega V., Arimoto S., “Adaptive control for robot manipulators with sliding mode error
coordinate system: free and constrained motions”, Robotics and Automation, IEEE Interna-
tional Conference on Volume 1, May 21-27, 1995.
16. Parra-Vega V., and Arimoto S., “Adaptive control for robot manipulators with sliding mode
error coordinate system: Free and constrained motions”, IEEE Int. Conf. Robotics and Au-
tomation, Japan, 1995.
17. http://www.sensable.com/
18. Simon Dan, Optimal State Estimation, A JOHN WILEY and SONS, INC. 2006
19. Spong Mark W., Vidyasagar M., Robot Dynamics and Control, John Wiley and Sons, 1989.
20. Vidyasagar M., Nonlinear Systems Analysis, Prentice-Hall, Englewood Cliffs, 1993.
21. Wen J.T., Bayard D., “A New Class of Stabilizing Control Laws for Robotic Manipulators,
Part I: Nonadaptive Case”, International Journal of Control, vol. 47, pp.1361-1385, 1988.
Capı́tulo 3
Diseño, Construcción e Instrumentación de un
Robot Móvil Diferencial
79
80 R. Silva-Ortigoza et al.
3.1. Introducción
robot humanoide para fines de investigación sobre el movimiento del cuerpo hu-
mano. En 1999 en el CMU, Zeglin propuso un nuevo diseño de robot con una pata
llamado Bow Leg Hopper [15], un diseño que por su naturaleza permite almacenar
la energı́a potencial de la pata, ya que es una articulación de forma curva que al dar
cada paso comprime un sistema de resorte y de esta forma el consumo de energı́a
es mı́nimo, permitiendo un ahorro significativo de la misma. En 2006, Hollis et al.
desarrollaron el robot Ballbot [16], un sistema holónomo cuyo movimiento es pro-
porcionado por una esfera ubicada en la parte inferior de la estructura; el control
de la esfera, es similar al de un mouse de computadora, empleando encoders para
determinar cada posición de la misma. Sin embargo, el estudio de este tipo de robots
con una esfera, fue iniciado por Koshiyama y Yamafuji [17] en 1991 (diversas in-
vestigaciones han sido efectuadas sobre este tipo de robots en los últimos años [18]).
Actualmente los robots teleoperados Spirit rover y Opportunity rover, (ver [19]), se
encuentran explorando la superficie del planeta marte en busca de mantos acuı́feros.
Los robots aquı́ mencionados, son únicamente una porción de los tantos que se han
diseñado, sin embargo, es posible notar que las aplicaciones de éstos son vastas y
que las mismas son ilimitadas debido al desarrollo cada vez más vertiginoso de la
tecnologı́a.
Robot Móvil
Protección
3.3.1.1. Motores
Puesto que se emplearán dos motores, cada uno tiene una fuerza asociada, Fmder
para el motor derecho y Fmizq para el motor izquierdo, en donde cada uno de estos
proporcionará la mitad de la fuerza necesaria, es decir:
F 46.745 N
Fmder = Fmizq = = = 23.373 N,
2 2
y, sabiendo que el diámetro de las ruedas es de 0.15 m (0.075 m de radio), se calcula
el par requerido para cada motor:
υmáx = ωrueda r,
de donde:
υmáx 1 m/s
ωrueda = = = 13.3333 rad/seg = 127.3267 rpm.
r 0.075 m
Esta velocidad angular no es la asociada a cada motor, ya que al acoplarse a la
caja de engranes, ésta es la que determinará dicha velocidad. Teniendo los datos
anteriores, se procede a calcular la potencia requerida Prueda para cada rueda:
Esta potencia, es la que debe ser transmitida a cada rueda del robot móvil, por lo
que la potencia del mismo debe ser superior debido a la pérdida en el reductor (caja
de engranes). La eficiencia promedio de los reductores, es aproximadamente del
70 %, lo que permite calcular una estimación de la potencia deseada de cada motor
Pdeseada :
Prueda = 0.7Pdeseada ,
de donde:
Prueda 23.373 W
Pdeseada = = = 33.39 W.
0.7 0.7
En base a la documentación técnica proporcionada por el fabricante [23], la primera
regla de selección de un motor establece que este debe de tener entre 1.5 y 2 veces
la potencia deseada Pdeseada , por lo tanto:
3. Diseño, Construcción e Instrumentación de un Robot Móvil Diferencial 87
se tiene:
Pdeseada 33.39 W
Mmotor = = = 0.125 Nm.
ωmotor 266.666 rad/seg
Para lograr el movimiento del robot móvil, es necesario acoplar una caja reduc-
tora a la salida de cada uno de los actuadores (motores) para que, de esta forma, sea
posible aumentar el par de salida (torque) y reducir la velocidad de cada actuador.
Una caja reductora, esta compuesta internamente por engranes que pueden variar en
tamaño, forma y número de dientes. Actualmente, el uso de las cajas reductoras de
engranes, es casi en su totalidad para amplificar el torque, mejorar la inercia permi-
sible y la reducción de vibración interna en un motor, respectivamente. En la Figura
3.3, se presenta la relación entre dos engranes acoplados entre sı́ para apreciar la
relación de reducción de velocidad que existe entre ambos.
La determinación del factor de reducción de velocidad, se calcula mediante la
razón entre la salida respecto a la entrada, es decir:
es decir, que el factor de reducción es de 2:1, lo que significa que por 2 revoluciones
del piñón, a la salida se tendrá una revolución.
Para el caso en que se tengan más de dos engranes, la determinación del factor
de reducción de velocidad, se realiza multiplicando la razón de las relaciones entre
cada engrane (ver Figura 3.4) mediante (3.1):
µ ¶ µ ¶ µ ¶
48 60 96
reducción = · · = 8.57 ≈ 9
14 32 72
es decir, que el factor de reducción es de 9:1, lo que significa que por 9 revoluciones
del piñón, a la salida se tendrá una revolución.
96 dientes
48 dientes
32 dientes
Eje de salida
Eje de entrada
14 dientes
72 dientes
60 dientes
(a) Vista frontal (b) Vista isométrica
3.3.1.3. Sensores
Las aplicaciones de los sensores son muy amplias. Los parámetros que se deseen
monitorear dentro de un sistema de cualquier tipo, son precisamente aquellos que
permitirán una adecuada selección del dispositivo. Los parámetros más comunes
que se toman en cuenta en la selección adecuada de un sensor, se muestran a con-
tinuación (ver [25]). En [26], se detalla cada una de estas caracterı́sticas.
Rango.- diferencia entre el valor máximo y el valor mı́nimo del parámetro a
detectar.
Resolución.- la variación más pequeña del sistema que el sensor puede detectar.
Precisión.- diferencia entre el valor medido y el valor real.
Sensibilidad.- razón de cambio a la salida por cada unidad de cambio a la entrada.
Cero offset.- valor nulo a la salida para una valor nulo a la entrada.
Linealidad.- porcentaje de desviación respecto a la curva de linealidad ideal del
dispositivo.
Tiempo de respuesta.- el retardo temporal entre la entrada y la salida.
Ancho de banda.- frecuencia en la que la magnitud de la salida desciende por
debajo de 3 DB.
Resonancia.- frecuencia a la que se presenta la magnitud de la salida más alta.
Temperatura de operación.- rango en el que el sensor opera de forma segura.
Banda muerta.- rango de valores a la entrada para los que no existe salida.
Sensibilidad al ruido.- razón entre la magnitud de la señal y el ruido de salida.
3. Diseño, Construcción e Instrumentación de un Robot Móvil Diferencial 91
Los sensores a emplear para el prototipo y que irán acoplados directamente a los
actuadores, son encoders ópticos. Los encoders, son actualmente los dispositivos
más empleados en robótica móvil para determinar la odometrı́a, ya que éstos per-
miten medir la velocidad angular y la posición del eje de un motor. Básicamente, un
encoder está formado por un arreglo de fototransistor y fotodiodo, dispuestos en am-
bos extremos de un disco que tiene sectores opacos y libres, que al ser acoplado a un
actuador rotatorio (motor), interrumpe el haz de luz que recibe el fototransistor del
fotodiodo y, de esta manera, genera un tren de pulsos (señal digital) que son envia-
dos a un controlador. El número de pulsos por revolución, determina la resolución
del encoder. En la Figura 3.6, se muestra a grandes rasgos como está constituido
internamente un encoder óptico.
Disco perforado
Circuito de
procesamiento Señal digital
digital de
señales
Emisor de luz Detector de luz
3.3.1.4. Baterı́a
Es preciso puntualizar que el diseño presentado en la Figura 3.9, sufrió una evolu-
ción paulatina conforme se fue adentrando más en las caracterı́sticas mecánicas
propiamente dichas. Esto permitió de forma inherente optimizar completamente el
diseño y lograr un aprovechamiento total del espacio. En los siguientes párrafos se
presentan detalladamente los pasos requeridos para efectuar el maquinado de las
piezas; se presentan los diseños en SolidWorksr y su contraparte maquinada.
Inicialmente, se adquirió el siguiente material:
Hoja de aluminio para la base del robot.
Barra de aluminio para los las ruedas motrices.
Barra de latón para buje que une la flecha de cada motor con cada rueda motriz.
Solera de aluminio para soportes de las ruedas de bola.
3. Diseño, Construcción e Instrumentación de un Robot Móvil Diferencial 95
Barras de aluminio para unir los soportes de las ruedas de bola a la base del robot.
Solera de aluminio para sujetadores de rodamientos y sujetadores de encoders.
En la Figura 3.10, se muestra el material utilizado a lo largo del proceso de diseño
para conseguir manufacturar la estructura cinemática, las ruedas motrices y los com-
ponentes mecánicos extra (soportes, sujetadores, bujes, etc.).
de manera significativa a la complejidad del maquinado. Por otro lado, ambas ruedas
constan de una pequeña caja en sus respectivos centros en donde irá acoplado el buje
(tapón de latón) que conectará ambas ruedas con las flechas de los motores. En la
Figura 3.12, se muestran los detalles antes expuestos.
Se observa en la Figura 3.12(b), las dos ruedas antes y después del vulcanizado.
También es posible apreciar en la rueda vulcanizada el detalle mencionado en el
párrafo anterior, el referido a la siglas del Instituto.
La pieza mecánica que permite el giro de las ruedas motrices, es el buje de unión.
La intención de éste, es acoplar las ruedas motrices y la flecha de los motores.
3. Diseño, Construcción e Instrumentación de un Robot Móvil Diferencial 97
(c) Diseño final de los bujes (d) Rueda con buje acoplado
Los rodamientos permiten la rotación sin traslación de las ruedas motrices al mo-
mento del giro de la flecha de cada motor. Los sujetadores tienen una doble función,
la de empotramiento de los rodamientos y la de unión con la base del robot móvil.
Para su maquinado, se empleó una solera de aluminio que fue cortada a las dimen-
siones requeridas (ver Figuras 3.16(b), 3.16(c)); asimismo, se empleó una máquina
fresadora para delimitar la caja en la que serán empotrados los rodamientos.
Para concluir esta sección se presenta la Figura 3.18, en la que se observan tanto
la evolución en el diseño propuesto en SolidWorksr como el diseño final de la
estructura mecánica del prototipo de robot móvil.
Hasta este punto se ha explicado a grandes rasgos el diseño mecánico del robot
móvil, lo que permite abordar en este apartado el diseño electrónico del sistema.
El diseño electrónico es imprescindible en esta etapa de la construcción, debido
a que de éste depende el funcionamiento eficaz del prototipo del robot móvil. Se
describirán los aspectos más importantes relacionados a la interconexión entre la
etapa de potencia y las etapas tanto de adquisición de datos ası́ como de actuación
(i. e. subsistemas) respectivamente. Cabe mencionar que los cálculos matemáticos
necesarios para la selección de algunos componentes (i. e. resistencias, capacitores)
se basaron en la ley de Ohm y en las leyes de Kirchoff (ver [29]) y que no se detallan
en este capı́tulo.
Se retoma en este apartado la Figura 3.2 en la que se muestra un diagrama a
bloques general del sistema. Es de nuestro interés en este punto la etapa central del
robot móvil, es decir, la etapa de potencia. A grandes rasgos, el diseño de la etapa
de potencia tiene dos vertientes, las cuales se muestran en los diagramas a bloques
de las Figuras 3.19 y 3.20 respectivamente.
Cabe señalar que para este prototipo, un diagrama a bloques más preciso es el
que se presenta en la Figura 3.20. La diferencia entre ambos diagramas, radica en
3. Diseño, Construcción e Instrumentación de un Robot Móvil Diferencial 101
Etapa de potencia
Subetapa 1
Circuito
Fuente
Etapa de potencia
Subetapa 1
Circuito
Fuente
JP7
F1 3
2
1
F2
JP1 24V
4
JP8
3
2 1
1 2
3
Bateria 24V
GNA24
SW1
JP3
1
2P-1T-2P 2
3
JP2 D1 VR1 Volt Reg 12V
4
4 Vin Vout
12V
3 GND
1N4007 JP4
2
1 C1 C3 1
0.33uF 0.1uF
2
Bateria 13V
3
4
GNA12
VCC5
R3
C1
330
.1uF U3A
R1
1 2
330
VCC5 MC74HC04AN
U1
6n137 14
DSpace VDD
7 PWM - Dir
GND
1
VCC5 1
2
2
3
4 Puente H
C2 R4
PWM - Dir
.1uF 330
U3D
R2
9 8
330
VCC5
5V MC74HC04AN
U2 VCC5
1 6n137
2 14
VDD
3 7
GND
4
++--
C7
Cap
.1uF
VCC VCC
B0 B1
T1 T2
MOTOR
B1 B0
T3 T4
B0 B1 B0 B1
T1 Sentido T2 T1 Sentido T2
de de
giro giro
MOTOR MOTOR
B1 B0 B1 B0
T3 T4 T3 T4
B0 B1 B0 B1
T1 Sentido T2 T1 T2
de
giro
Sin giro
MOTOR MOTOR
B1 B0 B1 B0
T3 T4 T3 T4
Para llevar a cabo la tarea de monitoreo respecto a la posición del móvil dentro
del espacio de trabajo, se hace uso de encoders (ver Figura 3.7). Dentro del diseño
mecánico, es importante la distribución de todos los componentes, y en el caso par-
ticular de los encoders (puesto que se acoplan directamente a la flecha del motor a
través de una pieza de acoplamiento), la disposición final se muestra en la Figura
3.26 resaltados en un rectángulo.
El circuito que gobierna el funcionamiento de estos encoders, se muestra en la
Figura 3.27. Se observa el uso de cuatro conectores etiquetados como 5V, Señal En-
coder, 12V, DSP. El conector etiquetado como 5V, es el voltaje que se requiere para
activar la etapa de optoacoplamiento y es proporcionado por el DSP. Cabe men-
cionar que en este circuito de encoders, es igualmente empleado el optoacoplador
108 R. Silva-Ortigoza et al.
VCC12
R5
VCC24
330
U4
D1
6
LMD18200
MUR640
1 C3
1
2
Bootstrap1 10nF M1 M2
PWM - Dir 5
PWM 2
Out1 1
Ventilador 1
2
-+ 2 10
Out2
3 Motor
VCC12 Puente H Dir C4
11 10nF D2
12V Bootstrap2
9 MUR640
1 4 Flag
Brake
2
8
Sense R8
++
10K VCC5
R7
7
VCC24
24V R6 2K7
VCC5 10K
4
3
2
VCC24
1
+ + - -
C5 C6
1uF 1200uF
6N137, el cual entrega una señal invertida a la salida respecto a la entrada y por
ello se hace uso de una compuerta lógica NOT. El conector etiquetado por Señal
Encoder, corresponde a las señales que envı́a el sensor hacia el DSP para lograr el
monitoreo del robot móvil. El conector 12V, es el voltaje proveniente de la fuente de
poder y permite energizar al encoder. Finalmente el conector DSP, es el que recibe
las señales provenientes de la etapa de optoacoplamiento y son enviadas directa-
mente a la tarjeta DSP. Se reitera el uso de dos circuitos de encoder, pues al tener
dos actuadores, es indispensable monitorear ambos.
3. Diseño, Construcción e Instrumentación de un Robot Móvil Diferencial 109
VCC VDD
14
VDD
VDD
R1 R7 7
GND
VDD 4K7 4K7
5V U1A
1 R4
1 2
2
3 1K C1
4 MC54HC04AJ
5V (VDD-DSP) U2 Cap
6n137 0.1uF
VCC VDD
Señal Encoder VDD
1 R2 R8 DSP
2 4K7 4K7
3 6
4 U1D 5
R5 4
5 9 8
6 3
1K C2 2
Encoder VCC MC54HC04AJ 1
U3 Cap Señales a DSP
6n137 0.1uF
VCC VCC VDD
12V
1 VDD
R3 R9
2
4K7 4K7
3
4 U1F
12V (VCC) R6
13 12
1K C3
MC54HC04AJ
U4 Cap
6n137 0.1uF
Una vez que se tienen todos y cada uno de los circuitos que permiten el fun-
cionamiento del prototipo de robot móvil, es necesaria una interfaz que permita la
comunicación eficaz entre la tarjeta DSP y el sistema. Para ello se ha diseñado un
3. Diseño, Construcción e Instrumentación de un Robot Móvil Diferencial 111
Pte. H Der.
PWM D. OutBF
1
Dir D. OutBF
2
3
4
Pt.H Der.
Se observa que la interfaz se lleva a cabo mediante un conector tipo DB25 cuyas
lı́neas se encuentran directamente acopladas a la tarjeta DSP. Se emplean igualmente
conectores que permiten direccionar las señales hacia los actuadores, sensores e in-
dicadores. Los conectores etiquetados como Fte. Enc. Izq., Enc. Izq., Pte. H Izq.,
son los correspondientes a las señales del voltaje requerido por el encoder izquier-
do, las señales enviadas desde el encoder hacia el DSP y finalmente las señales de
activación del puente H izquierdo respectivamente. De igual forma, los conectores
etiquetados como Fte. Enc. Der., Enc. Der., Pte. H Der., son los correspondientes a
las señales del voltaje requerido por el encoder izquierdo, las señales enviadas desde
el encoder hacia el DSP y finalmente las señales de activación del puente H dere-
cho respectivamente. Se hace uso de un buffer 74HC541 para la activación de los
distintos indicadores que permiten visualizar la dirección de cada motor, ası́ como
el promedio de voltaje enviado a cada motor (PWM) y a los puentes H izquier-
do y derecho, respectivamente. La configuración del conector DB25 respecto a las
señales que requiere el sistema para su funcionamiento, se presenta en la tabla de la
Figura 3.30.
Cada una de las señales provenientes del DSP es interceptada por un panel de
pruebas (rack), el cual permite no solo visualizar los led de las señales que son
112 R. Silva-Ortigoza et al.
empleadas, sino también facilita la conexión hacia la tarjeta interfaz del prototipo.
En la Figura 3.31 se observa el panel de pruebas conectado al robot móvil.
3.7. Conclusiones
Agradecimientos
Referencias
18. A. Meghdari, M. Rajaei, A. Tadayon and Y. Radparvar, “Geometric and control of a spheri-
cal mobile robot”, Proceedings of the 2nd. WSEAS International Conference on Applied and
Theoritical Mechanics, Venice, Italy, November 20-22, 2006.
19. Z. H. Aronson, T. Lechler, R. R. Reilly and A. J. Shenhar, “Project spirit, a strategic concept,
management of engineering and technology”, PICMET 2001, Portland International Confer-
ence, Vol. 1, July 29 - Aug 2, 2001.
20. P. F. Muir, “Modeling and control of wheeled mobile robots”, Ph.D. dissertation, The Robotics
Institute Carnegie Mellon University, Pittsburgh, August 1988.
21. D. Todd,Walking machines: An introduction to legged robotics, Kogan-Page, London, 1985.
22. T. Iwamoto, H. Yamamoto and K. Honma, “Transformable Crawler mechanism with adap-
tability to terrain variations”, International Conference on Advanced Robotics, Tokyo Japan,
September 1983.
23. MicroMo Electronics, Faulhaber Group, http://www.micromo.com.
24. R. Siegwart and I. R. Nourbakhsh, Introduction to mobile robots, The MIT Press, Cambridge,
MA, 2004.
25. R. H. Bishop, Mechatronics an Introduction, CRC Press, 2005.
26. J. S. Wilson, Sensor Technology Handbook, Elsevier, 2005.
27. J. R. Garcı́a Sánchez, “Diseño y construcción de un robot móvil aplicando el método de cam-
pos potenciales en la evasión de obstáculos”, Tesis de Maestrı́a, Sección Mecatrónica, CIDE-
TEC, México, 2008.
28. V. R. Barrientos Sotelo, “Análisis, diseño, construcción y control en tiempo real de un robot
móvil tipo Shakey en el seguimiento de trayectoria”, Tesis de Maestrı́a, Sección Mecatrónica,
CIDETEC, México, 2008.
29. R. L. Boylestad, L. Nashelsky, Electrónica: Teorı́a de circuitos, Pearson Educación, 1997.
30. J. Lovine, Robot, Androids and Animatrons: 12 Incredible Projects You Can Build, McGraw-
Hill, 2002.
Capı́tulo 4
Control Predictivo de una Planta Hidroeléctrica
4.1. Introducción
G. A. Muñoz-Hernández
Instituto Tecnológico de Puebla.
Av. Tecnológico 420, Col. Maravillas, Puebla, Puebla, 72220 MEXICO
e-mail: gmunoz@ieee.org
D. I. Jones
GWEFR, Centre for Advanced Software Technology. Technium Wales network
Ffordd Penlan Parc Menai Bangor Paı́s de Gales, LL574HJ U.K.
e-mail: dewi@gwefr.co.uk
S. P. Mansoor
University of Bangor. School of Informatics.
Dean Street, Bangor Paı́s de Gales, LL57 1UT U.K.
e-mail: s.mansoor@bangor.ac.uk
117
118 G. A. Muñoz-Hernandez et al.
X1 ∨ X2 ⇒ δ1 + δ2 ≥ 1
X1 ∧ X2 ⇒ δ1 = 1, δ2 = 1 (4.1)
X1 ⊕ X2 ⇒ δ1 + δ2 = 1
donde X1 y X2 son variables lógicas (verdadero o falso) y δ1 y δ2 son variables
enteras (0 ó 1). Un conjunto muy completo de estas relaciones ha sido establecido
por Cavalier et al., [9].
4. Control Predictivo de una Planta Hidroeléctrica 121
Para ilustrar las caracterı́sticas del MLD, permı́tanos considerar un modelo li-
neal simple del subsistema hidráulico de Dinorwig, el cual consistirá de sólo dos
unidades, con funciones de transferencia idénticas G j (s) con salida (∆ Pmech ), con-
trolada por la apertura en la válvula (∆ G). Aunque la dinámica de acoplamiento
cruzado no es modelada, dicho acoplamiento hidráulico aparece como resultado del
cambio en las condiciones de operación y por tanto en el tiempo de arranque del
túnel principal (Twti ). Todas las variables están en el sistema por unidad (p.u.) nor-
malizado a 300 MW y 50 Hz.
δ1 = 1 si p1 = 0 ∧ p2 = 0
δ2 = 1 si p1 = 0 ∧ p2 = 1
(4.7)
δ3 = 1 si p1 = 1 ∧ p2 = 0
δ4 = 1 si p1 = 1 ∧ p2 = 1
las cuales siguen la condición: ⊕4i=1 [δi = 1]. Definiendo:
entonces el sistema (4.2) puede ser descrito como un sistema MLD en tiempo dis-
creto con T = 0.25, ver ecuación (4.8).
4
y j (k) = ∑ zi (k) (4.8)
i=1
z1 (k) = 0
z2 (k) = 0
z3 (k) = 0.45y(k − 1) − 2.36u(k) + 3.01u(k − 1)
z4 (k) = 0
entonces:
y1,2 (k) = 0.45y(k − 1) − 2.36u(k) + 3.01u(k − 1)
donde 1 ≥ y1,2 (k) ≥ 0.
Por otro lado, si h = 0.91 y Uo = 1, entonces: δ1 = 1, y
entonces:
y1 (k) = 0.45y(k − 1) − 2.10u(k) + 2.68u(k − 1)
y2 (k) = 0
donde 1 ≥ y1 (k) ≥ 0.
El sistema es por tanto descrito de forma diferente dependiendo del valor de δ i .
Por ejemplo si δ 1 = 1 sólo habrı́a una unidad activa y la cabeza hidráulica serı́a baja,
causando una ganancia de lazo abierto baja y por lo tanto un valor final también
bajo. En cambio si δ 3 = 1 la cabeza hidráulica serı́a más alta y ambas unidades
estarı́an activas, incrementándose el tiempo de arranque, dando como resultado una
constante de tiempo más larga, originando una respuesta más lenta. En uno de los
modelos de la siguiente sección se incluirán en el modelo de la planta, la válvula
4. Control Predictivo de una Planta Hidroeléctrica 123
en MVA, del Generador. Es conveniente recordar que los modelos son expresados
en el sistema por unidad (p.u.), normalizado a 300 MW y 50 Hz.
# $
! ! ! "! %
& ' ! ( )
* + , % %
%
- %
$ %! !
#
" "
&
$ '
∆ Pes (s) 1
= (4.10)
∆ Pe (s) s+1
Este modelo fue validado a través de todo el rango operacional con una simu-
lación no lineal independiente de Dinorwig, la cual ha sido validada usando res-
puestas medidas experimentalmente [1, 12].
En un sistema real algunos de sus parámetros varı́an aun cerca de un punto de
operación, G0 y q0 , por lo tanto los resultados de simulación pueden ser precisos
solamente cerca del punto de operación seleccionado durante el diseño del proceso.
Para cambios grandes de velocidad y potencia, un modelo no lineal, que asume un
fluido no comprimible, un conducto rı́gido y una cabeza hidráulica sin restricciones,
fue considerada [3]. La potencia mecánica de salida de la turbina ∆ Pm mueve al
subsistema eléctrico, cuya función de transferencia (4.11) es representada por la
bien conocida ecuación de swing [2]:
∆ Pe (s) Ks ω0 /2H
= (4.11)
∆ Pm (s) s2 + KD s/2H + Ks ω0 /2H
Los modelos no lineales son usados en este estudio para representar a la planta,
para el control predictivo se utilizará como base de cálculo un modelo lineal, esto
permitirá mantener al mı́nimo el tiempo de cómputo. En la planta de potencia con
múltiples unidades las unidades que están en lı́nea reaccionan, mediante sus regu-
ladores, a perturbaciones en la presión y flujo (y por lo tanto potencia) causadas por
otras unidades. Las unidades que no están en lı́nea tienen sus válvulas de control ce-
rradas y por tanto no interactúan. La estructura de la planta, sin embargo, varı́a con
el tiempo dependiendo del número de unidades activas. Para este modelo multiva-
riable (MIMO) linealizado de Dinorwig, una columna de agua no elástica y pérdidas
insignificantes son consideradas, además todas las seis unidades son tratadas como
idénticas (aunque es sabido que menores diferencias ocurren en la práctica debido
a las tolerancias en la manufactura). La razón de cambio en el flujo en el tubo de
conexión (penstock) puede ser determinado como:
dq gA
= (h0 − h − hi ) (4.12)
dt l
donde (p.u. es la notación por unidad):
A - sección de área a través del túnel.
ho - cabeza fija de la columna de agua p.u.
126 G. A. Muñoz-Hernandez et al.
La suma de las corrientes en las conexiones individuales debe ser igual al total
de la corriente en el túnel común, entonces usando (4.12) la ecuación de momento
para el agua en el túnel común es:
µ ¶
l dq1 dq2 dqn
h0 − h = + +...+ (4.14)
Aga dt dt dt
En la ecuación (4.17) los dos tubos de conexión son considerados iguales, esto es
Twt = Twt1 = Twt2 . Un modelo MIMO, Figura 4.6, el cual incluye las seis máquinas,
será usado más tarde en simulaciones. En esta Figura ∆ G es la apertura de la válvu-
la de control, ∆ Pmech es la potencia mecánica de las turbinas y ∆ Pe es la potencia
eléctrica. El regulador, que está actualmente en operación, incluye un control indi-
vidual clásico PI en cada turbina. El bloque de hidrodinámica del modelo extendido
4. Control Predictivo de una Planta Hidroeléctrica 127
Unidad n está en lı́nea cuando Pdi > 0 y Unidad n está fuera de lı́nea cuando
[(∆ Pdi = 0) y (∆ Pmi = 0)].
U0 Gi (s) Xi (s)
−2.358s+3.395
1 0.076s3 +0.8204s2 +2.788s+3.031
0
−2.358s2 −5.454s+14.96 −8.559s
2 0.076s +1.26s +7.213s +16.69s+13.35 0.076s +1.26s +7.213s2 +16.69s+13.35
4 3 2 4 3
−2.358s2 −1.986s+11.01 −6.301s
3 0.076s4 +1.221s3 +6.643s2 +14.1s+9.83 0.076s4 +1.221s3 +6.643s2 +14.1s+9.83
−2.358s2 +0.03428s+8.711 −4.985s
4 0.076s4 +1.198s3 +6.311s2 +12.59s+7.778 0.076s4 +1.198s3 +6.311s2 +12.59s+7.778
−2.358s2 +1.357s+7.207 −4.124s
5 0.076s4 +1.183s3 +6.093s2 +11.6s+6.435 0.076s4 +1.183s3 +6.093s2 +11.6s+6.435
2
−2.358s +2.289s+6.145 −3.517s
6 0.076s4 +1.173s3 +5.94s2 +10.9s+5.487 0.076s4 +1.173s3 +5.94s2 +10.9s+5.487
el punto de operación fijo a 0.8 p.u. La Figura 4.8a muestra las respuestas con una
unidad activa, mientras que la Figura 4.8b muestra las respuestas con seis unidades
en operación. Observe que los modelos tienen un comportamiento diferente depen-
diendo del número de unidades activas. La caracterı́stica de fase no mı́nima es más
marcada en ambos casos lineales; esta y otras diferencias entre los modelos son prin-
cipalmente porque el modelo lineal no está tomando en cuenta otros factores, tales
como las pérdidas por fricción. Sin embargo, el modelo MLD tiene un desempeño
mejor o igual al modelo lineal fijo y este sigue al modelo no lineal con casi la mis-
ma precisión a pesar del número de unidades activas. Esta caracterı́stica adaptiva es
extremadamente importante cuando un control es diseñado para mantener el mismo
desempeño bajo diferentes condiciones de operación. El modelo MLD muestra una
respuesta con errores máximos de 8 % (6 unidades) y 6 % (1 unidad) con respecto
a la respuesta del modelo no lineal, mientras que aunque el modelo fijo presenta el
mismo error que el modelo MLD para el caso de seis unidades, que fue donde se
fijó, este modelo fijo presenta un error máximo de 13 % cuando una unidad está en
operación. Lo cual es más del doble de lo que presenta el modelo MLD en estas
mismas condiciones.
Figura 4.8 Respuesta al escalón de lazo abierto de los modelos MLD y no lineal con el punto de
operación en 0.8. (a) una unidad en operación y (b) seis unidades en operación.
4. Control Predictivo de una Planta Hidroeléctrica 131
Figura 4.9 Respuesta directa y de acoplamiento-cruzado. (a) Unidad uno y (b) unidades 2-6 en
sincronı́a.
N2 Nu
J(N1 , N2 , Nu ) = ∑ [ŷ(t + j|t) − w(t + j)]2 Q + ∑ [∆ u(t + j − 1)]2 R (4.20)
j=N1 j=1
e(t)
A(q−1 )y(t) = q−d B(q−1 )u(t − 1) +C(q−1 ) (4.21)
1 − q−1
134 G. A. Muñoz-Hernandez et al.
donde A(q−1 ) y C(q−1 ) son matrices de polinomios mónicos de orden nxn, mientras
que B(q−1 ) es una matriz polinomial de orden nxm; e(t) es el error. Las matrices A,
B y C fueron obtenidas usando un modelo MLD. Este modelo MLD produce una
representación CARIMA de intervalos lineales, cuya estructura y valores exactos
dependen del número de unidades activas de la estación hidroeléctrica bajo estu-
dio. El modelo MLD también depende del punto de operación. Es común formular
las leyes de control con la suposición de que los estados de las plantas y los con-
troles poseen un rango ilimitado. Para resolver este problema, GPC suministra una
solución computacionalmente eficiente en el control óptimo de la señal:
El modelo MLD representa de mejor manera a la planta, esto nos permite usar
un horizonte de predicción más corto, decrementando con esto el tiempo de cálculo,
además los factores de peso del control son reducidos (λ ), acelerando la respuesta
del sistema. El modelo MLD de Dinorwig, discutido en las secciones 4.3 y 4.4, es
empleado en el control predictivo MLD-GPC para el cálculo de las predicciones.
El MLD-GPC fue sintonizado con N = 10 (horizonte de predicción) y Nu = 10
4. Control Predictivo de una Planta Hidroeléctrica 135
Figura 4.10 Respuesta a entradas rampas grandes con diferentes cabezas hidráulicas del lago (H),
para una (a) y seis (b) unidades en operación.
p.u. para la Figura 4.13b. Con este controlador, la respuesta primaria permanece
prácticamente constante en todos los casos.
4.7. Conclusiones
Los resultados de simulación han mostrado como el MLD-GPC puede ser apli-
cado para mejorar las caracterı́sticas de respuesta de Dinorwig. El tomar en cuenta
la naturaleza multivariable de la planta mejora directamente la respuesta a tran-
sientes, con respecto al desempeño de la planta bajo el controlador Scheduled PI.
La selección de diferentes valores de λ dependiendo del número de unidades activas
permite al MLD-GPC producir una respuesta rápida en cada modo de operación. La
inclusión de restricciones de saturación en la señal de salida y en el crecimiento de
la señal de control, lleva a una rápida respuesta sin sobretiros en el caso de una so-
la unidad en operación, sin comprometer la estabilidad cuando múltiples unidades
están en lı́nea. El MLD-GPC permite también un buen desempeño cuando diferentes
cabezas hidráulicas del lago son probadas en la planta de poder. En general estos re-
sultados son muy prometedores y nos permiten visualizar trabajos futuros emplean-
4. Control Predictivo de una Planta Hidroeléctrica 137
Figura 4.11 Respuesta a escalones pequeños con diferentes cabezas hidráulicas del lago (H), para
una (a) y seis (b) unidades en operación.
Referencias
Figura 4.12 Respuesta a escalones pequeños bajo el controlador Scheduled PI con diferentes
cabezas hidráulicas del lago (H), para una (a) y seis (b) unidades en operación.
Figura 4.13 Respuesta a escalones pequeños bajo el controlador MLD-GPC con diferentes
cabezas hidráulicas del lago (H), para una (a) y seis (b) unidades en operación.
16. J. M. Goncalves, A. Megretski and M. A. Dahleh, “Global Analysis of Piecewise Linear Sys-
tems using impact maps and surface Lyapunov functions”, IEEE Transactions on Automatic
Control, 2003, Vol. 48, No. 12, pp. 2089-2106.
17. M. Johannson and A. Rantzer, “Computation of piecewise quadratic Lyapunov functions for
hybrid systems”, IEEE Transaction on Automatic Control, 1998, Vol. 43, pp. 555-559.
18. S. Yokokawa, Y. Ueki, H. Tanaka, H. Doi, K. Ueda and N. Taniguchi, “Multivariable Adaptive
control for a thermal generator”, IEEE Transactions on Energy Conversion, July 1988, Vol. 3,
No. 3, pp. 479-486.
19. P. Idowu and A. Ghandakly, “A coordinated adaptive stabilizer for multimachine power sys-
tems”, Proceedings of the Twenty-First Annual North-American Power Symposium, 1989, pp.
175-184.
20. J. E. Lansberry and L. Wozniak, “Adaptive Hydrogenerator Governor tuning with a Genetic
Algorithm”, IEEE Transaction on Energy Conversion, 1994, Vol. 9, No. 1, pp. 179-186.
21. S. Bonnet and L. Wozniak, “Adaptive speed control of Hydrogenerators by recursive least
squares identification algorithm ”, IEEE Transaction on Energy Conversion, 1995, Vol. 10,
No. 1, pp. 162-168.
22. O. P. Malik and Y. Zeng, “Design a robust adaptive controller for a water turbine governing
system”, IEEE Transactions on Energy Conversion, 1995, Vol. 10, No. 2, pp. 354-359.
23. Q. Lu, “SDM hybrid control approach for nonlinear systems and application to power sys-
tems”, International Conference on Power Control, 2002, China, pp.19-23.
140 G. A. Muñoz-Hernandez et al.
24. J. Thomas, D. Dumur and D. Buisson, “Moving horizon state estimation of hybrid systems.
Application to fault detection of sensors of a steam generator”, IEEE Int. Conf. on Control
Applications, 2003, pp. 1375-1380.
25. G. Ramond, D. Dumur, A. Libaux and P. Boucher, “Direct Adaptive Predictive Control of
an hydro-electric plant”, 2001, IEEE Int. Conf. on Control Applications, Mexico, D. F., pp.
606-611.
26. G. Ferrari-Trecate, E. Gallestey, P. Leticia, M. Spedicato, M. Morari and M. Antoine, “Mod-
eling and Control of Co-Generation Power Plants: A Hybrid System Approach”, IEEE Trans-
actions on Control Systems Technology, 2004, Vol. 12, No. 5, pp. 694-705.
27. E. F. Camacho and C. Bordons, “Springer-Verlag“Model Predictive Control ”, Springer-
Verlag, 1999.
28. D. W. Clarke and C. Mohtadi P. S. Tuffs, “Generalized Predictive Control (Part I The basic
algorithm) ”, Automatica, 1987, Vol. 23, pp. 137-148.
29. G. A. Munoz-Hernandez, D. I. Jones and S. I. Fuentes-Goiz, “Modelling and Simulation of a
hydroelectric power station using MLD”, CONIELECOMP, 2005, pp. 83-88.
30. D. I. Jones, S. P. Mansoor, F. C. Aris, G. R. Jones, D. A. Bradley and D. J. King, “A standard
method for specifying the response of hydroelectric plant in frequency-control mode”, Electric
Power Systems Research, 2004, Vol. 68, pp. 19-32.
Capı́tulo 5
Lazos de Corriente en Sistemas
Electromecánicos
V. M. Hernández-Guzmán
5.1. Introducción
V. M. Hernández-Guzmán
Universidad Autónoma de Querétaro.
Facultad de Ingenierı́a.
A.P. 3-24, Querétaro, Qro. 76150, MEXICO
e-mail: vmhg@uaq.mx
141
142 V. M. Hernández-Guzmán
[5], [6]. Esta técnica de control utiliza un lazo interno de corriente con el que se
trata de conseguir que el par generado por los motores alcance rápidamente el par
deseado. Esta consigna de par está especificada por un lazo externo que actúa sobre
la variable a controlar, es decir, por un controlador PID o PD de posición. Aunque
en la literatura académica se argumenta que la inductancia es pequeña, sin embargo
en la práctica industrial esto no es necesario y el control de par o de corriente es
el método que permite suponer que el par es la señal de control y que la dinámica
eléctrica de los motores puede ser despreciada.
Sin embargo, aún en el presente, el diseño usando control de par o de corriente se
realiza utilizando ideas intuitivas porque no se han propuesto métodos analı́ticos que
presenten justificaciones formales para dicha estrategia de control. En este capı́tu-
lo se presenta un estudio formal que el autor ha propuesto en la literatura [7], [8]
el cual prueba la estabilidad en lazo cerrado de sistemas electromecánicos cuando
son controlados usando control de par o de corriente. Este estudio de estabilidad es
desarrollado tomando en consideración la dinámica eléctrica de los motores utiliza-
dos. Esto significa que el análisis aquı́ presentado evita la posibilidad de los malos
desempeños y de las inestabilidades que pueden aparecer cuando no se considera la
dinámica eléctrica de los motores durante el diseño.
El presente capı́tulo se divide en dos partes. En la sección 5.2 se presenta el
problema de controlar motores de CD con escobillas mientras que en la sección 5.3
se presenta el problema de controlar robots manipuladores equipados con motores
de CD sin escobillas. En ambos casos se presentan resultados en simulación que
permiten apreciar el desempeño obtenido.
Ra L
+ i + T Ê T c n1
u ea Jm bL
Tp
à à bm JL
n2 T L ò
El modelo matemático del conjunto motor-carga está dado por [9], pp. 139:
di
L = u − Ra i − n ke θ̇ (5.1)
dt
J θ̈ = −b θ̇ + n km i − Tp (5.2)
donde:
n2
Θ = n θ, n= (5.3)
n1
J = n2 Jm + JL , b = n2 bm + bL
Por otro lado, nótese que la parte mecánica (5.2) se puede escribir como:
µ Z t ¶
J θ̈ + b θ̇ = n km i − n km R−1 k p θ̃ + ki θ̃ (σ )d σ
0
µ Z t ¶
+n km R−1 k p θ̃ + ki θ̃ (σ )d σ − Tp
0
Z t
= n km ρ + k p θ̃ + ki θ̃ (σ )d σ − Tp (5.8)
0
Nótese que θ̃˙ = −θ̇ porque θd es una constante y, por tanto, sθ̃ (s) = −sθ (s). Apli-
cando la transformada de Laplace a (5.7) y (5.10) y usando la observación anterior:
+ nk m òe(s)
0 Js 3+bs 2+k p s+k i
s
à
+ 1 òe(s)
à ú(s) nk m Js 2+bs
à
k ps+k i
s
Esto significa que los polos de la función de transferencia en (5.12) pueden ser
determinados usando la técnica del lugar de las raı́ces en el sistema en lazo cerrado
de la figura 5.3. En la figura 5.4 se muestra el lugar de las raı́ces correspondiente, el
cual se obtiene de la siguiente manera. La función de transferencia de lazo abierto
en la figura 5.3 está dada como:
s − (−ki /k p ) 1/J
kp (5.14)
s s(s − (−b/J))
146 V. M. Hernández-Guzmán
Im (s)
kp ! 1
!n
kp ! 1
+ 90 î
à b=J ûa
ï Re (s)
0
à k i=k p à 90 î
àa
!n
ï
kp ! 1
Figura 5.4 Lugar de las raı́ces del sistema en lazo cerrado de la figura 5.3.
Nótese que k p es el valor que el método del lugar de las raı́ces varı́a desde 0 hasta
+∞ y que mientras esto sucede también cambia el valor de ki de manera que el
cociente ki /k p permanece constante. Usando fórmulas propias del método del lugar
de las raı́ces, [9] Caps. 6 y 7, se obtiene:
Por otro lado, aplicando el criterio de Routh, [9] Cap. 5, al polinomio Js3 + b s2 +
k p s + ki se obtiene la tabla 5.1. De acuerdo a (5.15), (5.17), la tabla 5.1 y la figura
5.4, se concluye lo siguiente para los polos de la función de transferencia en (5.13)
y, por tanto, de la función de transferencia en (5.12):
Los tres polos tienen parte real negativa de modo que dos serán complejos con-
jugados y uno real si:
5. Lazos de Corriente 147
s3 J kp
s2 b ki
bk p −Jki
s1 b 0
s0 ki
Tabla 5.1 Aplicación del criterio de Routh al polinomio Js3 + b s2 + k p s + ki .
b ki
> , ki > 0 (5.18)
J kp
Dos polos complejos conjugados tienen parte real positiva mientras que el tercero
es real y negativo si:
b ki
< , ki > 0 (5.19)
J kp
Hay dos polos con parte real cero sobre el eje imaginario y el tercero tiene parte
negativa si:
b ki
= , ki > 0 (5.20)
J kp
Es conveniente resaltar que las condiciones bJ < kkpi y ki > 0 son las que normal-
mente se cumplen debido a que el coeficiente de fricción viscosa b comúnmente
es pequeño mientras que la ganancia integral ki comúnmente se selecciona grande
para obtener respuestas transitorias rápidas. Ası́ que se puede afirmar que la fun-
ción de transferencia en (5.12) tiene dos polos complejos conjugados con parte real
positiva y otro real y negativo. Suponga que los polos complejos conjugados están
a una distancia ωn del origen y que el polo real está ubicado en s = −a donde
b/J < a < ki /k p . En el lugar de las raı́ces de la figura 5.4 se puede apreciar que para
valores suficientemente grandes de k p se tiene que ωn > a.
Con estos antecedentes, ahora se estudiará la estabilidad del sistema en lazo ce-
rrado (5.11), (5.12), es decir, el mostrado en la figura 5.2. Para esto se hará uso del
criterio de estabilidad de Nyquist [9] Cap. 8. De acuerdo a la discusión anterior, la
función de transferencia en (5.12) se puede escribir como sigue si bJ < kkpi y ki > 0:
n km s n km a ωn2
= s (5.21)
Js3 + b s2 + k p s + ki Jaωn2 (s + a) (s2 − 2ζ ωn s + ωn2 )
En la figura 5.5 se muestran las gráficas de Bode, [9] Cap.8, de cada una de las
funciones de transferencia que componen a la función de transferencia en (5.22)
ası́ como la gráfica de Bode de G(s)H(s). A partir de esto se obtiene la gráfica
polar, [9] Cap. 8, que se muestra en la figura 5.6. Nótese que el término de segundo
orden introduce un adelanto de fase, lo cual es debido a que contiene dos polos
inestables. Como G(s)H(s) tiene dos ceros sobre el eje imaginario el criterio de
estabilidad de Nyquist requiere que se dibuje la gráfica polar de G(s)H(s) obtenida
al recorrer la trayectoria de Nyquist modificada, [9] Cap. 8, que en este caso incluye
un semicı́rculo al rededor de s = 0. No es difı́cil comprobar que la gráfica polar
correspondiente es idéntica a la gráfica mostrada en la figura 5.6 ya que los dos
ceros que tiene G(s)H(s) sobre el eje imaginario están en s = 0.
dB
ii)
vi)
i)
a !n R=L
0 !
1
iii) iv) v)
fase vi)
[grados]
ii)
+ 180
iv)
a !n R=L
0 !
iii) v)
à 90
Im (G(j!)H(j!))
ï
(à 1; j0) 0 Re (G(j!)H(j!))
es necesario que N = −2, es decir, que el punto (−1, j0) del plano mostrado en
la figura 5.6 sea rodeado dos veces en sentido antihorario. De acuerdo a la gráfica
polar en la figura 5.6 y a las gráficas de Bode en la figura 5.5, esto es posible si
ωn > a, es decir si k p es suficientemente grande, porque esto introduce el atraso de
fase necesario para conseguir N = −2.
Por otro lado, es común que en la práctica se utilice un valor grande de ra , lo cual
resulta en un valor grande de R. Bajo dichas condiciones se tiene que kd R À n2 ke km
y, por tanto, la función de transferencia en (5.22) se reduce a:
kd R/L a ωn2
G(s)H(s) = s2 (5.23)
Jaωn s + R/L (s + a) (s − 2ζ ωn s + ωn2 )
2 2
Nótese que para conseguir dos rodeos antihorarios al punto (−1, j0) en la figura
5.6 es necesario un valor suficientemente grande del factor Jakωd 2 . Por otro lado, tal
n
como se ha mencionado, el valor de ωn2 crece al aumentar k p , lo cual es necesario
para conseguir la condición ωn > a. Por tanto, es necesario utilizar también un valor
suficiente grande de kd .
De acuerdo a lo anterior, se ha conseguido probar el siguiente resultado.
Proposición 5.1. El polinomio caracterı́stico del sistema en lazo cerrado (5.1),
(5.2), (5.4) tiene todas sus raı́ces con parte real negativa si las ganancias ra , k p ,
kd , ki se seleccionan de manera que la gráfica polar de la siguiente función de
transferencia:
150 V. M. Hernández-Guzmán
rodea dos veces en sentido antihorario al punto (−1, j0). Una manera de conseguir
esto es la siguiente:
Defı́nase:
b ki
Elı́jase k p y ki de manera que J < kp y ki > 0.
Elı́jase k p suficientemente grande de modo que ωn > a, donde Js3 + b s2 + k p s +
ki = J(s + a)(s2 − 2ζ ωn s + ωn2 ) para algún 0 < ζ < 1.
Elı́jase kd suficientemente grande de manera que haya dos rodeos en sentido
antihorario al rededor del punto (−1, j0).
u = γ (i∗ − i) (5.26)
donde γ es una constante positiva mientras que i∗ representa el valor que la corriente
eléctrica i debe alcanzar para generar el par deseado τ ∗ , es decir:
1 ∗
i∗ = τ (5.27)
nkm
En el caso en el que se utiliza un controlador PID para establecer el par deseado se
tiene que:
Z t
τ ∗ = κ p θ̃ − κd θ̇ + κi θ̃ (σ )d σ (5.28)
0
Nótese que se puede obtener el controlador (5.4) a partir de (5.26), (5.27) y (5.28)
si se establece lo siguiente:
γ γ γ
γ = ra , kp = κp, kd = κd , ki = κi (5.29)
nkm nkm nkm
Es importante señalar que en la práctica se seleccionan valores muy grandes de ra ,
[10] pp.76, de manera que ra À Ra y, por tanto, R ≈ ra = γ . Por otro lado, usando
las expresiones en (5.25) se obtiene:
5. Lazos de Corriente 151
R L R L R
kp = kp + ki , kd = kd + kp, ki = ki (5.30)
nkm nkm nkm nkm nkm
De acuerdo al párrafo anterior R ≈ γ es muy grande y, por tanto, se pueden des-
preciar los segundos términos que están en el lado derecho de estas expresiones.
Entonces, comparando (5.30) y (5.29) se concluye que κ p ≈ k p , κd ≈ kd , κi ≈ ki .
Ahora se presenta un ejemplo en simulación que, sin embargo, utiliza los valo-
res numéricos de un motor real. Este motor es el modelo MT-4060 ALYBE [11]:
L = 9.6 × 10−3 [H], km = 0.573[Nm/A], ke = 0.573[V/(rad/s)], Ra = 2.3[Ohm],
Jm = 12.43 × 10−4 [kgm2 ], bm = 0.001[Nm/(rad/s)]. Además se supone que exis-
te una carga cuyos valores son los siguientes JL = 120.43 × 10−4 [kgm2 ], bL =
0.001[Nm/(rad/seg)], n = 10. Se utilizan las siguientes ganancias: ra = 50, k p =
182.6318, kd = 18.2883, ki = 456.37. Aunque estas ganancias pudieran parecer muy
grandes, sin embargo con ellas se obtienen los valores kd = 2, k p = 20, ki = 50 cuan-
do se utilizan las relaciones en (5.25). Recuérdese que, de acuerdo al comentario 5.1,
las ganancias kd , k p , ki son las que obtienen cuando el diseño se hace suponiendo
que el par es la señal de entrada. En las figuras 5.7 y 5.8 se muestra la gráfica polar
y las gráficas de Bode de la función de transferencia en (5.24) cuando se usan los
valores numéricos aquı́ propuestos. En la figura 5.8 se aprecia que se consiguen dos
rodeos antihorarios al rededor del punto (−1, j0), por lo que el sistema de control
es estable. Más aún, en la figura 5.7 se puede observar que el margen de fase es del
orden de 50 grados por lo que es muy aceptable. Es conveniente aclarar lo siguien-
te. Aunque la gráfica 5.8 parece ser diferente a la bosquejada en la figura 5.6, sin
embargo todo es debido a la visualización de la figura: si se hace un acercamiento
al rededor del origen se podrá observar que la curva es similar a la de la figura 5.6.
Prueba de ello es que las gráficas de Bode de las figuras 5.7 y 5.5 son muy parecidas.
En la figura 5.9 se presentan resultados en simulación cuando se usa θd = 0.8[rad] y
todas las condiciones iniciales iguales a cero. Es conveniente mencionar que existe
un pico de voltaje u muy grande (al rededor de 150[V]) al principio de la simulación
pero tiene una duración muy corta. Este pico se debe a las grandes ganancias que
se utilizan. Sin embargo, se reduce rápidamente porque la diferencia i∗ − i disminu-
ye rápidamente (véase el comentario 5.1) con lo que se obtiene de nuevo un valor
(pequeño) razonable de u. Debido a la corta duración que tiene este pico de volta-
je los resultados no tienen ningún cambio cuando se satura el voltaje a valores de
±30[V]. En estas simulaciones se aplica una perturbación constante Tp = 2[Nm] a
152 V. M. Hernández-Guzmán
Bode Diagrams
20
−20
−40
Phase (deg); Magnitude (dB)
−60
−80
−100
250
200
150
−2 0 2 4 5
10 10 10 10 10
Frequency (rad/sec)
Figura 5.7 Gráficas de Bode de la función de transferencia en (5.24) usando las ganancias uti-
lizadas por el controlador diseñado.
Nyquist Diagrams
2
Imaginary Axis
−2
−4
−6
−10 −8 −6 −4 −2 0
Real Axis
Figura 5.8 Gráfica polar de la función de transferencia en (5.24) usando las ganancias utilizadas
por el controlador diseñado.
5. Lazos de Corriente 153
òe
[rad]
u
[V]
t [s]
define como kAk = λmax (AT A). El sı́mbolo p = (d/dt) se utiliza para representar
el operador derivada.
donde:
r
3
KT 1 = N p (Lb − La ), KT 2 = Np KB (5.34)
2
Va = [va1 , va2 , . . . , van ]T ∈ Rn
Vb = [vb1 , vb2 , . . . , vbn ]T ∈ Rn
Ia = [ia1 , ia2 , . . . , ian ]T ∈ Rn
Ib = [ib1 , ib2 , . . . , ibn ]T ∈ Rn
IA = diag{ia1 , ia2 , . . . , ian } ∈ Rn×n
IB = diag{ib1 , ib2 , . . . , ibn } ∈ Rn×n
5. Lazos de Corriente 155
El vector q ∈ R n representa las posiciones en cada una de las uniones del robot,
M(q) es una matriz de n × n simétrica y definida positiva conocida como la matriz
de inercia, C(q, q̇)q̇ es el término de Coriolis y de efectos centrı́fugos, g(q) = ∂ U(q) ∂q
es el término de efectos gravitacionales, donde U(q) es una función escalar que
representa la energı́a potencial del robot, y F es una matriz de n × n constante,
diagonal y definida positiva que representa los coeficientes de fricción viscosa en
cada una de las uniones. A lo largo de esta sección se usa el vector q̃ = qd − q para
representar el error de posición en las uniones donde qd ∈ R n es un vector constante
que representa las posiciones que se desean en cada unión. También se supone que
el robot que se estudia está equipado solamente con uniones rotativas.
El modelo (5.31)-(5.34) se obtiene después de que se ha aplicado una transforma-
ción DQ (también conocida como transformación de Park) sobre el modelo original
de cada uno de los motores los cuales se supone tienen tres fases conectadas en es-
trella [17], [14], [18]. Por tanto, Va y Vb representan, respectivamente, los voltajes
de las fases DQ asociadas con cada uno de los motores. Ia y Ib son las corrientes
eléctricas correspondientes. La , Lb , Ra , Rb , N p , KB son matrices constantes, diago-
nales, definidas positivas. Se recomienda leer las referencias [12], [17] para conocer
una descripción completa de estas matrices. Finalmente, KT 1 y KT 2 son matrices
diagonales (KT 2 es definida positiva) que contienen las constantes de par de los mo-
tores, mientras que τ = [KT 1 IB + KT 2 ]Ia es el par aplicado por los motores en las
uniones del robot.
Sean V j = [v j1 , v j2 , v j3 ]T ∈ R 3 e I j = [i j1 , i j2 , i j3 ]T ∈ R 3 , respectivamente, los
voltajes y las corrientes de fase del motor de CD sin escobillas colocado en la unión
j, el cual es trifásico y está conectado en estrella. La aplicación de la transformación
DQ (o de Park) implica que (véase la referencia [17], pp. 373):
· ¸ r · ¸
ζa j 2 ζ j1 cos(n p j q j ) + ζ j2 cos(n p j q j − 23π ) + ζ j3 cos(n p j q j + 23π )
=
ζb j 3 ζ j1 sin(n p j q j ) + ζ j2 sin(n p j q j − 23π ) + ζ j3 sin(n p j q j + 23π )
r
v j1 va j cos(n p j q j ) + vb j sin(n p j q j )
v j2 = 2 va j cos(n p j q j − 2π ) + vb j sin(n p j q j − 2π )
3 3 3
v j3 va j cos(n p j q j + 23π ) + vb j sin(n p j q j + 23π )
∂ g(q)
Propiedad 5.2. [21] pp. 120, [22] pp. 170. La matriz ∂q es simétrica porque
∂ U(q)
g(q) = ∂q es el gradiente de una función escalar U(q) que es dos veces continua-
mente diferenciable. Por otro lado, todos los eigenvalores de una matriz simétrica
son reales, es decir, todos los eigenvalores de ∂ g(q) n
∂ q son reales para toda q ∈ R .
Propiedad 5.3. [19], [23], [20] pp. 102. Existen constantes positivas kg y kc tales
que para toda w, y, z, q ∈ R n , se tiene que:
A continuación se listan algunos resultados que son útiles para el desarrollo pre-
sentado en esta sección.
Teorema 5.1. (Weyl)[22] pp. 184. Sean A y B dos matrices simétricas de n × n.
Suponga que los eigenvalores de A de B y de A + B se ordenan de manera creciente,
es decir de menor a mayor: λmin = λ1 ≤ λ2 ≤ · · · ≤ λn−1 ≤ λn = λmax . Entonces,
para cada par de enteros j, k tales que 1 ≤ j ≤ n, 1 ≤ k ≤ n y j + k ≤ n + 1 se tiene
que:
Lema 5.2. (Barbalat)[24] pp. 21. Si g(t) es una función vectorial del tiempo de
n componentes tal que g(t) ∈ L∞ , ġ(t) ∈ L∞ y g(t) ∈ L p para alguna p ∈ [1, ∞),
entonces g(t) → 0 conforme t → ∞. Se dice que g(t) ∈ L p para p ∈ [1, ∞) si el valor
de la siguiente expresión existe:
5. Lazos de Corriente 157
µZ ∞
¶1/p
kg(t)k p dt (5.42)
0
Teorema 5.2. [22] pp. 402. Una matriz es definida positiva si y sólo si todos sus
eigenvalores son positivos.
Va = −ra Ia + KP q̃ − KD ϑ + RKT−1
2 g(qd ) (5.47)
b b
V = −Q̇∆a θ1 + ε Q̃IA θ2 (5.48)
" #b · ¸
d θb1 IB Q̇ δa
=Γ (5.49)
dt θb2 −ε IB Q̃ Ia
donde:
158 V. M. Hernández-Guzmán
ε0
ε= , R = Ra + ra
1 + kq̃k
½ ¾
bi p
ϑ = diag q (5.50)
p + ai
Q̃ = diag{q̃1 , q̃2 , . . . , q̃n } ∈ R n×n
Q̇ = diag{q̇1 , q̇2 , . . . , q̇n } ∈ R n×n
∆a = diag{δa1 , δa2 , . . . , δan } ∈ R n×n (5.51)
δa = KP q̃ − KD ϑ + RKT−1 2 g(qd ) = [δa1 , δa2 , . . . , δan ] ∈ R
T n
Comentario 5.2. Es importante subrayar que el filtro en (5.50) puede ser construido
en la práctica sin necesidad de medir la velocidad utilizando la siguiente realización:
ẋ = −Ax − ABq
ϑ = x + Bq
Nótese que ϑ̇ = −Aϑ + Bq̇, es una realización del filtro definido en (5.50). Usando
estas expresiones y sustituyendo (5.47) en (5.32) se obtiene:
La ρ̇ = −Rρ − N p Lb IB q̇ − KT 2 q̇
+La R (KP + KD B)q̇ − La R−1 KD Aϑ
−1
(5.54)
donde θ̃1 = θb1 − θ1∗ , θ̃2 = θb2 − θ2∗ . Nótese que si se define:
KP = KT 2 R−1 KP (5.56)
KD = KT 2 R−1 KD (5.57)
Por tanto, la dinámica en lazo cerrado está dada por (5.61), (5.54), (5.59), (5.62)
junto con:
1 T 2 1 2
q̃ [ KP − ε 2 M(q)]q̃ ≥ ε02 q̃T [ 2 KP − M(q)]q̃
2 ε2 2 ε2 ε0
porque esto asegura que la matriz [ ε 2ε 2 KP −M(q)] tiene todos sus eigenvalores posi-
2 0
tivos y, por tanto, es una matriz definida positiva (véase el teorema 5.2). Para obtener
5. Lazos de Corriente 161
(5.64) se ha usado el hecho de que siendo M(q) una matriz definida positiva entonces
−M(q) es definida negativa y su eigenvalor menor es igual a −λmax (M(q)). Por otro
lado, como KP es una matriz diagonal, entonces (5.64) se puede escribir como:
2 2
0< λmin (KP ) − λmax (M(q)) ≤ λmin ( 2 KP − M(q)) (5.65)
ε2 ε02 ε2 ε0
2 ∂ g(q)
KP +
ε1 ∂q
Usando de nuevo el teorema 5.1 y la propiedad 5.2 se encuentra que esto es cierto
si:
½ ¾
2 ∂ g(q)
0 < λmin (KP ) − λmax
ε1 ∂q
Por tanto, usando (5.38) se concluye que f (q̃) es positiva definida de manera global
si:
2λmin (KP )
> ε1 (5.67)
kg
asegura que ε11 + ε12 = 12 , lo cual es importante para el cálculo que sigue a con-
tinuación.
Usando (5.34) y el hecho de que todas las matrices involucradas son diagonales
se encuentra que q̇T KT 1 IB ρ − ρ T Np Lb IB q̇ + IbT N p La Q̇ρ = 0 y q̇T KT 1 IB KT−1
2 δa +
∗
−1 ∗ −1 ∗
Ib N p La Q̇KT 2 δa = Ib N p Lb Q̇KT 2 δa . Estos resultados, ası́ como las propiedades en
T T
largo de las trayectorias del sistema en lazo cerrado (5.61), (5.54), (5.59), (5.62),
(5.63) como:
+θ̃ T Γ −1 θ̃˙ + G,
G = q̇T N p Lb IB KT−1 ∗
2 δa − ε q̃ KT 1 IB ρ
T
−ε q̃T KT 1 IB KT−1 ∗
2 δa
Utilizando (5.37), (5.39) y (5.43)-(5.46) se encuentra que (véase también [20] pp.
193):
donde:
5. Lazos de Corriente 163
donde z = [kϑ k, kq̃k, kρ k]T , λmin (Ā) > 0 y ε > δ > 0 para alguna constante finita δ
porque q̃ está acotada. Integrando directamente se encuentra:
Z t
ν (0) − ν (t) ν (0)
kq̃(r)k2 dr ≤ ≤ < ∞, ∀t ≥ 0 (5.73)
0 δ λmin (Ā) δ λmin (Ā)
164 V. M. Hernández-Guzmán
Comentario 5.3. En la práctica es común suponer que el par generado por los mo-
tores de CD sin escobillas es proporcional a la corriente eléctrica. Por esta razón los
manejadores de estos motores incluyen algún controlador de corriente que asegure
que se genera el par deseado [6]. Suponiendo que La = Lb entonces el par aplicado
por los motores en las uniones del robot está dado como τ = KT 2 Ia . Entonces, los
manejadores de los motores entregan el siguiente voltaje:
Va = γ (Ia∗ − Ia ) (5.74)
donde γ es una matriz diagonal definida positiva e Ia∗ representa el valor de la corri-
ente eléctrica Ia que se necesita para general el par deseado τ ∗ , es decir:
Ia∗ = KT−1
2τ
∗
(5.75)
τ ∗ = κ p q̃ − κd ϑ + g(qd ) (5.76)
considerar una situación más cercana a una aplicación real. Es interesante resaltar
que la introducción de dicha saturación no afecta en nada a los resultados obtenidos.
0,4
[rad]
0,3
qe1
0,2
0,1
qe2
t [s]
Figura 5.10 Resultados en simulación. Errores de posición.
Referencias
1. A. Ailon, R. Lozano, and M.I. Gil’, “Iterative regulation of an electrically driven flexible-joint
robot with model uncertainty”, IEEE Trans. Robot. Autom., vol. 16, no. 6, pp. 863-870, Dec
2000.
2. M.C. Good, L.M. Sweet, and K.L. Strobel, “Dynamic models for control system design of
integrated robot and drive systems”, Journal of Dynamic Systems, Measurement, and Control,
Vol. 107, pp. 53-59, 1985.
3. T.-J. Tarn, A.K. Bejczy, X. Yun, and Z. Li, “Effect of motor dynamics on nonlinear feedback
robot arm control”, IEEE Trans. Robot. Autom., vol. 7, no. 1, pp. 114-122, Feb. 1991.
4. S.D. Eppinger and W.P. Seering, “Introduction to dynamic models for robot force control”,
IEEE Control Systems Magazine, pp. 48-52, 1987.
5. R. Campa, E. Torres, V. Santibáñez and R. Vargas, “Electromechanical dynamics characteriza-
tion of brushless direct-drive servomotor”, Proc. VII Mexican Congress on Robotics, COMRob
2005, México, D.F., October 27-28, 2005.
5. Lazos de Corriente 167
20
[V]
18
16
14
12
10
6 V a1
4
2
V a2
0
t [s]
Figura 5.11 Resultados en simulación. Voltajes en la fase a.
V b2
0
-1
-2
-3
-4
-5
-6
V b1
-7
-8
0 0,2 0,4 0,6 0,8 1 1,2 1,4
t [s]
Figura 5.12 Resultados en simulación. Voltajes en la fase b.
15. N. Hemati, J. Thorp and M. C. Leu, Robust Nonlinear control of brushless DC motors for
direct-drive robotic applications, IEEE Transactions on Industrial Electronics, Vol. 37, No. 6,
pp. 460-468, December 1990.
16. R. Ortega, A Lorı́a, P. Nicklasson and H. Sira-Ramı́rez, Passivity-based control of Euler-
Lagrange Systems, Springer, London, 1998.
17. D. M. Dawson, J. Hu and T. C. Burg, Nonlinear control of electric machinery, Marcel Dekker,
New York, 1998.
18. P.C. Krause, O. Wasynczuk, S.D. Sudhoff, Analysis of electic machinery and drive systems,
IEEE Press, 2002.
19. R. Kelly, A tuning procedure for stable PID control of robot manipulators, Robotica, Vol. 13,
pp. 141-148, 1995.
20. R. Kelly, V. Santibáñez, and A. Lorı́a. Control of robot manipulators in joint space. London:
Springer, 2005.
21. H. K. Khalil, Nonlinear Systems, 3rd. Edition, Prentice Hall, Upper Saddle River, 2002.
22. R.A. Horn and C.R. Johnson, Matrix analysis, Cambridge University Press, New York, 1993.
23. R. Kelly, Global positioning for robot manipulators via PD control plus a class of nonlinear
integral actions, IEEE Transactions on Automatic Control, Vol. 43, No. 7, pp. 934-938, July
1998.
24. S. Sastry and M. Bodson, “Adaptive control: stability, convergence and robustness”, Prentice
Hall, Englewood Cliffs, 1989.
25. R. Kelly, Comments on: Adaptive PD controller for robot manipulators, IEEE Transcaction
on Robotics and Automation, Vol. 9, No. 1, pp. 117-119, February, 1993.
Capı́tulo 6
Instrumentación Virtual
E. Juárez-Guerra
Universidad Autónoma de Tlaxcala.
Facultad de Ciencias Básicas. Ingenierı́a y Tecnologı́a.
Calzada Apizaquito s/n Km. 1.5, Apizaco, Tlaxcala, México. CP.90300
e-mail: ejuarez@ingenieria.uatx.mx
E. A. Portilla-Flores
CIDETEC-IPN. Departamento de Posgrado. Área de Mecatrónica
Av. Juan de Dios Bátiz s/n. Esq. Miguel Othón de Mendizábal,
“Unidad Profesional Adolfo López Mateos”. C.P. 07700, México, D.F., MÉXICO
e-mail: aportilla@ipn.mx
A. Cano-Corona
Universidad Politécnica de Tlaxcala.
Carretera Federal Tlaxcala-Zacatelco, Km. 9.5, San Lorenzo Axocomanitla, Tlaxcala, México. CP.
e-mail: arico 96@yahoo.com
C. Sánchez-Olavarrı́a
Universidad del Altiplano.
Eucalipto No. 1, Col. El Sabinal. Tlaxcala, Tlaxcala. 90102, México.
e-mail: cesar ari@hotmail.com
169
170 E. Juárez-Guerra et al.
6.1. Introducción
Detectan y corrigen fallas que ocurren con el tiempo en los elementos de proce-
samiento analógico, como son los capacitores, OP-AMPS, etc.
Codifican y transmiten en forma digital la información, además de detectar y
corregir errores que surgen en la transmisión de información.
La figura 6.4 muestra las imágenes de algunos instrumentos inteligentes.
La tendencia de hace unos 10 años fue que las computadoras iniciaran un proce-
so de servir como mecanismo de instrumentación. Actualmente se refuerza esta ten-
dencia y se puede observar en cualquier equipo electrónico, el cual ahora cuenta con
6. Instrumentación Virtual 175
todos los elementos necesarios para que el instrumento funcione de acuerdo a las
necesidades del usuario.
Bajos costo y reusabilidad, ya que el software minimiza el desarrollo y los costos
de mantenimiento.
La figura 6.7 muestra las imágenes de algunos instrumentos virtuales.
tienen los mismos resultados que con un instrumento real tradicional, siendo ası́ que
el usuario está trabajando con un instrumento virtual.
La ventaja de sustituir elementos de hardware por software es la de permitir que
el instrumento o sistema virtual desarrollado sea flexible en su funcionalidad, pudi-
endo ser modificado en el momento en que se desee. Para realizar un instrumento
virtual es necesario hardware especı́fico y un software con el cual se podrá progra-
mar la funcionalidad del instrumento y configurar al hardware para leer los datos
requeridos para realizar una medición. En la mayorı́a de las ocasiones el usuario
sólo tiene acceso a una representación gráfica o panel de control en pantalla de la
computadora, donde se cuenta con controles e indicadores para manipular el instru-
mento virtual.
El concepto de Instrumentación virtual básicamente involucra la adquisición de
señales, el análisis, procesamiento, almacenamiento, distribución de los datos y pre-
sentación de los resultados o información al usuario, todo esto relacionado con la
medición de una o varias señales o variables fı́sicas. Además se le puede agregar
el diseño de la interfaz hombre-máquina, la visualización, monitoreo y supervisión
remota del proceso, la comunicación con otros dispositivos o equipos, etc.
A partir de esto, se puede decir que la Instrumentación virtual está enfocada al
desarrollo de instrumentos o sistemas de medición encargados de medir señales,
registrar datos y tomar decisiones (acciones de control), basados en un sistema de
cómputo, mediante hardware especializado y software especı́fico para la progra-
mación de la funcionalidad del sistema o instrumento.
Es evidente que para la instrumentación virtual se requiere al menos de dos etapas
de hardware esenciales, una para la entrada y otra para la salida. En la etapa de
entrada se requieren generalmente transductores de entrada, mejor conocidos como
sensores, que son los encargados de transformar la variable fı́sica a medir a una señal
eléctrica. En la etapa de salida se requieren transductores de salida, mejor conocidos
como actuadores, que son los encargados de transformar la señal eléctrica en otra
señal que tal vez no necesariamente sea eléctrica, esto dependerá de la señal de
salida que se desee. Tanto en la etapa de entrada como en la de salida se pueden
requerir circuitos de acondicionamiento de señal.
La parte fundamental de la Instrumentación virtual es el software. Una de las
empresas que inició con el concepto de la instrumentación virtual es National Ins-
truments. En 1983 Truchard y Kodosky de National Instruments, decidieron crear
un software que permitiera utilizar la computadora personal como un instrumento
para realizar mediciones. Tres años más tarde crearon la primera versión de soft-
ware que permitió, de una manera gráfica y sencilla, diseñar un instrumento de
medición basado en la computadora [4]. La figura 6.8 muestra el concepto de la
instrumentación virtual, donde se puede observar que el usuario interactúa con un
sistema de cómputo y a través de él se realizan mediciones y en su caso se toman
decisiones en un proceso en el cuál no se está en forma presencial.
El software que desarrolló National Instruments para la programación de los
instrumentos virtuales se conoce como Laboratory Virtual Instrument Engineering
Workbench, más comúnmente conocido por las siglas LabVIEW. También existen
178 E. Juárez-Guerra et al.
otros software dedicados a la Instrumentación virtual tales como Cyber tools, Dasy-
Lab, Matlab Simulink, DAP Measurement Studio, Agilent VEE y otros.
Un instrumento tradicional se caracteriza por realizar una o varias funciones es-
pecı́ficas que no pueden ser modificadas, mientras que un Instrumento virtual es
una combinación de elementos de hardware y software usados en un sistema de
cómputo, que cumple las mismas funciones que el instrumento tradicional. A di-
ferencia de un instrumento convencional o tradicional, un instrumento virtual es
altamente flexible y puede ser diseñado por el usuario de acuerdo a sus necesidades
y sus funciones pueden ser cambiadas a voluntad del mismo usuario modificando
el programa. Estas caracterı́sticas de los instrumentos virtuales los convierten en
una herramienta didáctica muy importante para aplicarse en el aprendizaje de los
estudiantes de las ciencias y de la ingenierı́a [4].
6.3.2.9. Hardware
6.4.2. LabVIEW
especı́ficas sobre éstas áreas, por lo que logra obtener la aplicación sin sacrificar
rendimiento. Estas bibliotecas a la vez se dividen en grupos para la adquisición,
procesamiento, presentación y almacenamiento de datos. Otra caracterı́stica impor-
tante de LabVIEW es que se pueden crear interfaces interactivas de usuario, de tal
manera que la apariencia sea amigable y de fácil comprensión para el usuario y
generalmente muy similar al sistema o instrumento tradicional de medición.
Clark, Cockrum, Ibrahim y Smith (1994), señalan que LabVIEW es un lenguaje
de programación gráfica, que se ejecuta a velocidades comparables con programas
compilados en C; igualmente mencionan que un instrumento virtual es un módulo
de software, realizado gráficamente para que parezca un instrumento fı́sico; tiene
un panel frontal que sirve como interface interactiva para entradas y salidas, un
diagrama de bloques que determina la funcionalidad del instrumento virtual [4].
Estos autores resaltan como caracterı́stica muy importante del LabVIEW que, por
ser conceptualmente simple, los estudiantes se pueden centrar en el contenido básico
del experimento, no perdiendo tiempo en actividades menos importantes, como la
recolección de datos.
A continuación se muestra un breve resumen de la evolución de LabVIEW:
1986 LabVIEW 1.0 para Macintosh.
1990 LabVIEW 2.0 para Macintosh.
1992 LabVIEW para Sun y para Windows.
1993 LabVIEW 3.0. Versión para múltiples plataformas.
1994 LabVIEW 3.1. Versiones para HP-UX y PowerMac.
1997 LabVIEW 4.1. DAQ Wizards sobre Windows.
1998 LabVIEW 5.0. ActiveX, multihilos.
1999 LabVIEW 5.1. Gráficas 3D, herramientas Web, Linux.
2000 LabVIEW 6i. Mediciones inteligentes por internet.
2002 LabVIEW 6.1. Paneles remotos, manejo de eventos.
2003 LabVIEW 7.0. VI Express, Modulo de tiempo real.
2004. LabVIEW 7.1. VI Express adicionales.
2005. LabVIEW 8.0. VI Express e Inteligencia distribuida.
2006. LabVIEW 8.2. Programación orientada a objetos.
2007. LabVIEW 8.5. Herramientas de control de memoria y funciones Math-
script.
2008. LabVIEW 8.6. Programación paralela, Wi-Fi, Multinúcleo y Tecnologı́a
FPGA.
186 E. Juárez-Guerra et al.
LabVIEW puede comunicarse de forma sencilla con todo tipo de hardware, des-
de un Asistente digital personal (PDA, por sus siglas en inglés Personal Digital
Assistant) hasta una computadora o bien hasta con controladores lógicos programa-
bles. Ası́ mismo, también es compatible con otras aplicaciones, pudiendo compartir
datos a través del internet, AcviveX, DLLs, librerı́as compartidas, SQL, TCP/IP,
XML, OPC y otros. De la misma manera que LabVIEW es un software compilado
para optimizar el desempeño del sistema, dado que la velocidad de ejecución resulta
vital en cualquier proceso.
6. Instrumentación Virtual 189
Los programas en VEE son creados seleccionando y conectando entre sı́ objetos
desde la barra de menú. El resultado de esto es un diagrama de flujo de datos, el
cual es mucho más fácil de entender que un programa tradicional basado en texto.
VEE cuenta con un panel principal, el cual muestra sólo los objetos que el usuario
necesita para ejecutar el programa y para ver los resultados.
Este software está diseñado para ser una herramienta que pueda automatizar in-
herentemente tareas complejas. VEE es compatible con los programas de Microsoft,
tal es el caso de Word, Excel, Outlook y otros, de la misma manera reconoce lengua-
jes de programación como Visual Basic, Visual C++, LabVIEW y Matlab. VEE
utiliza LAN/GPIB y también RS-232 como puertos de enlace permitiendo comuni-
carse con otros dispositivos en cuanto a tecnologı́as de redes. Este software puede
aumentar la productividad reduciendo dramáticamente el tiempo de desarrollo en
aplicaciones.
Finalmente otra de las herramientas importantes de VEE es que permite simpli-
ficar la comunicación con instrumentos y otros componentes. Mediante el ”Instru-
ment Manager”del VEE se puede hacer una búsqueda en el bus de equipos conec-
tados a la computadora, sin la necesidad de crear enlaces de direcciones hacia los
instrumentos que están conectados al bus o leer controladores [5].
190 E. Juárez-Guerra et al.
6.4.4. CyberTools
Las licencias son totalmente modulares y de bajo costo, incluyendo librerı́as para
almacenamiento de datos, control y monitoreo de hardware y controladores indus-
triales, interconexión de aplicaciones Cyber Tool por medio de internet o intranet.
6.4.5. DASYLab
el tiempo para señales de control complejas. Además, los datos pueden guardarse
para ser analizados posteriormente por cualquier tipo de aplicación externa.
El programa está disponible en cuatro ediciones diferentes: DASYLab Ligth,
contiene las funciones básicas para la adquisición y representación de datos vı́a
computadora; DASYLab Basic, incluye funciones matemático-estadı́sticas, ası́ co-
mo un conjunto de módulos de control básicos; DASYLab Full, proporciona un nue-
vo conjunto de bloques adicionales que permiten automatizar las medidas, además
de tareas de análisis; y DASYLab Professional que además contiene funcionalidades
de red, análisis de frecuencias, amplitudes y el módulo generador de setpoints [8].
DASYLab soporta una amplia variedad de diferentes dispositivos de adquisición
de datos utilizando cualquier tipo de interfaz que disponga la computadora. Entre
estos se pueden mencionar: PCI, PXI/Compact-PCI, USB, PC-Card, CAN, Ether-
net, RS-232, IEEE, SPS Simatic S7, etc. Además cabe mencionar que dispone de
módulos de extensión, tales como: Analysis Toolkit, Sequence generator, Acústica,
entre otros.
El nombre Matlab proviene de Matrix Laboratory. Este software que fue ini-
cialmente desarrollado para realizar operaciones con matrices muy fácilmente y ha
evolucionado hasta convertirse en una herramienta muy popular en diversos campos
de la ingenierı́a y la Ciencia. Matlab es un lenguaje de alto rendimiento para cálculo
técnico. El mismo integra cálculo, visualización y programación en un entorno de
fácil utilización en donde los problemas y las soluciones son expresadas en una no-
tación matemática familiar. Los usos más tı́picos incluyen: Cálculos matemáticos,
Desarrollo de algoritmos, Modelado, simulación y prototipos y Gráficas cientı́ficas
y de ingenierı́a.
Simulink es un ambiente interactivo para modelar, simular y analizar una amplia
variedad de sistemas dinámicos, pudiendo ser estos lineales, no lineales, discretos,
de tiempo continuo y sistemas mixtos. Permite realizar diagramas de bloques con
operaciones click-and-drag, cambiar parámetros del modelo y visualizar resultados
durante una simulación. Es también un sistema abierto, que permite al usuario es-
coger, adaptar y crear componentes o subistemas. Simulink se apoya en el ambiente
Técnico Computacional de Matlab. Este software y su grupo de Toolboxes ofre-
cen un conjunto amplio de herramientas de ingenierı́a y matemática para definir
algoritmos, analizando datos y visualizando resultados. Simulink y Matlab juntos
proveen un entorno integrado para construir modelos versátiles y simular dinámi-
cos, diseñando y probando ideas nuevas. La figura 6.20 muestra en ambiente de
desarrollo de Simulink [9].
Simulink es un paquete desarrollado por The MathWorks Inc., y es soportado
en ambientes como Unix, Macintosh y Windows. En este software los sistemas son
dibujados en pantalla como diagramas de bloques, tales como funciones de trans-
6. Instrumentación Virtual 193
ferencia, sumadores, uniones, etc., ası́ como entradas y salidas virtuales de equipos
tales como generadores, voltı́metros, osciloscopios, entre otros.
Los instrumentos virtuales brindan significativas ventajas en cada etapa del pro-
ceso de ingenierı́a, desde la investigación y el diseño hasta el ensayo de manufactura.
6.5.3. Manufactura
enseñanza que de ellos se derivan. Ertugrul afirma que, de acuerdo a las experiencias
de enseñanza-aprendizaje que se tienen actualmente en tecnologı́a utilizando com-
putadoras, éstas se clasifican en cuatro grupos: entrenamiento basado en computa-
doras, aprendizaje asistido por computadoras, instrucción asistida por computadoras
y experimentación asistida por computadora. Schär y Krueger definen el aprendiza-
je asistido por computadora como: “diferentes formas de métodos de enseñanza por
computadora en los cuales el estudiante tiene a la computadora como un profesor
virtual”[4].
Definitivamente, las mejores posibilidades para aprovechar las ventajas que
ofrece la instrumentación virtual se encuentra en la implementación de laboratorios.
Como se ha señalado, ella permite la realización de sistemas de medición basados
en sistemas de cómputo, que hacen posible a los ingenieros, los profesores, investi-
gadores y estudiantes resolver problemas de ingenierı́a o de las ciencias. Ası́ como
una hoja de cálculo le permite a un administrador solucionar problemas de adminis-
tración, la instrumentación virtual es también una solución a los problemas de cos-
tos y obsolescencia de los equipos en los laboratorios. Reemplazar los instrumentos
tradicionales por instrumentos virtuales que se ejecutan en computadoras, permite
que las funciones de los mismos vayan a la par del desarrollo de las nuevas tec-
nologı́as de las computadoras cuyos costos siguen una tendencia decreciente.
Los laboratorios son un elemento clave en la formación integral y actualizada de
un ingeniero. No se puede concebir un ingeniero que no haya realizado prácticas
de laboratorio en su trayectoria de formación inicial. Los avances tecnológicos de
los últimos años han abierto posibilidades para cambiar la estructura rı́gida de los
laboratorios tradicionales, por una estructura flexible que se apoya en las computa-
doras, circuitos de acondicionamiento, hardware de adquisición de datos y software.
Se puede afirmar que cada año, aumenta el número de universidades que se acogen
a esta propuesta de laboratorios virtuales. Ertugrul reporta más de 50 aplicaciones
de laboratorios basados en la instrumentación virtual en ingenierı́a eléctrica, elec-
trónica, mecánica, biomédica, control e instrumentación, quı́mica, medio ambiente,
instrumentación y control a través de internet, lo cuál garantiza que las universidades
formen profesionales con competencias para los nuevos desafı́os [4].
6.6. Conclusiones
Agradecimientos
Referencias
G. Saldaña-González
Resumen Este capı́tulo esta dirigido a lectores con poca experiencia en el área de
visión por computadora, para los más experimentados, parte del contenido expuesto
resultará familiar, sin embargo se presentará material relacionado a la utilización
de dispositivos programables para la implementación de los sistemas de visión, que
resulta ser un tópico novedoso que permite apreciar el potencial de esta tecnologı́a
para la implementación de aplicaciones en tiempo real. En este capı́tulo se introduce
al lector a los principios básicos de la visión por computadora, se presentan además
algunos ejemplos de la importancia de este tipo de sistemas en relación con otras
áreas del conocimiento, tales como la robótica y la mecatrónica. Se aborda la forma
en la cual se representan las imágenes, pasando por algunas caracterı́sticas de los
sensores más utilizados para la obtención y formación de la imagen. Finalmente
se presentan algunos algoritmos para extraer información a partir de una imagen
haciendo énfasis en las caracterı́sticas que hacen que dichos algoritmos sean ideales
para ser implementados en dispositivos digitales programables en campo FPGAs
(Field Programmable Gate Arrays).
7.1. Introducción
G. Saldaña-González
Universidad Tecnológica de Puebla
Antiguo Camino a la Resurrección 1002-A, Zona Industrial, Puebla, Puebla, 72300 MEXICO
e-mail: griselda.saldana@utpuebla.edu.mx
201
202 G. Saldaña-González
La visión por computadora es una rama de las ciencias computacionales que in-
volucra la extracción automática, manipulación, análisis y clasificación de imágenes
o secuencias de imágenes, generalmente dentro de un sistema de cómputo especı́fico
o de propósito general. El objetivo de estos sistemas consiste en obtener información
útil acerca del mundo utilizando imágenes para realizar alguna tarea.
A menudo el término visión por computadora se utiliza como un sinónimo
de máquina de visión, sin embargo existen elementos que los diferencian. Las
máquinas de visión generalmente se usan en ambientes industriales y pretenden dar
solución a problemas especı́ficos. Por otro lado un sistema de visión por computa-
dora, generalmente trata de duplicar el efecto de la visión humana mediante la per-
cepción electrónica y la comprensión de una imagen. La investigación en esta área
a menudo tiene que ver con la percepción psicológica y el análisis de los humanos
y otros sistemas de visión biológicos para generar un comportamiento cercano al
sistema de visión humano.
Una máquina de visión analiza una imagen y proporciona una descripción de su
contenido [7].
Esta descripción debe capturar los aspectos de los objetos en la imagen que son
útiles para realizar alguna tarea. Entonces se considera a una máquina de visión
como parte de una entidad más grande que interactúa con su medioambiente. El
sistema de visión se puede considerar como un elemento de retroalimentación que
se relaciona con el sensado mientras otros elementos están dedicados a la toma de
decisiones y a la implementación de esas decisiones. En la Figura 7.1 se presenta un
diagrama que muestra la relación existente entre una máquina de visión y un sistema
de visión.
204 G. Saldaña-González
Figura 7.1 Descripción de un sistema de una máquina de visión en relación a un sistema de visión.
mente es necesario procesar los datos para asegurar que satisfacen ciertas su-
posiciones implicadas en el método. Ejemplos del tipo de procesamiento que
puede realizarse son:
• Re-muestreo para asegurar que el sistema de coordenadas de la imagen es
correcto.
• Reducción de ruido para asegurar que el ruido del sensor no introduce infor-
mación falsa.
• Mejoramiento del contraste para asegurar que la información relevante se
puede detectar.
• Representación espacial de escala para mejorar la estructura de las imágenes
en escalas locales apropiadas.
Extracción de caracterı́sticas. Las caracterı́sticas de la imagen a varios niveles
de complejidad son extraı́das a partir de la imagen. Algunos ejemplos tı́picos
son:
• Lı́neas, bordes, esquinas.
• Localización de puntos de interés tales como esquinas o regiones.
• Caracterı́sticas más complejas tienen que ver con la textura, la forma o el
movimiento.
• Detección/Segmentación.
Detección/Segmentación. En algún punto en el proceso se debe tomar una de-
cisión sobre cuales puntos o regiones de la imagen son relevantes para realizar
más procesamiento. Algunos ejemplos son:
• Selección de un conjunto especı́fico de puntos de interés.
• Segmentación de una o múltiples regiones de la imagen que contienen un
objeto especı́fico de interés.
Procesamiento de alto nivel. En este paso la entrada es tı́picamente un pequeño
conjunto de datos, por ejemplo un conjunto de puntos o una región de la imagen
que se asume contiene un objeto especı́fico. El procesamiento restante trata, por
ejemplo con:
• Verificación que los datos satisfacen un modelo propuesto y suposiciones es-
pecı́ficas de la aplicación.
• Estimación de los parámetros especı́ficos de la aplicación, tales como la posi-
ción del objeto o su tamaño.
• Clasificación de un objeto detectado en diferentes categorı́as.
técnicas para estimar los detalles en las imágenes, relacionar las mediciones de de-
talles a la geometrı́a de objetos en espacio e interpretar la información geométrica.
Figura 7.4 Principales áreas del conocimiento relacionadas con la Visión por Computadora.
Figura 7.6 Cámara CCD donde las celdas convierten la energı́a luminosa en cargas eléctricas que
representan entradas para la PC.
Figura 7.7 Rol del Frame Grabber en el sistema de procesamiento de una imagen.
Entrada de Video.
Buffer de la imagen (Frame buffer).
Procesador digital de señales DSP.
Unidad de salida de video.
La estructura conformada por estos elementos se muestra en la Figura 7.8. Sin
embargo en el mercado existen disponibles una gran variedad de tipos de frame
grabbers por lo que es difı́cil describir una estructura tı́pica.
Figura 7.9 Digitalización de una imagen donde un pı́xel representa el valor de brillo en ese punto.
A menudo es posible encontrar casos en los que M=N=2K donde {K=8, 9, 10}.
Esto pude ser motivado por un circuito digital o por el uso de ciertos algoritmos
tales como la trasformada de Fourier y la transformada rápida de Fourier.
El número de niveles de gris usualmente es una potencia de 2, esto es, L=2B
donde B es el número de bits de representación binaria de la brillantes de los niveles.
Cuando B>1 se habla de una imagen en niveles de gris (o en escala de grises);
cuando B=1 se habla de una imagen binaria. En una imagen binaria sólo hay dos
niveles de gris lo cuales se pueden referenciar, por ejemplo como “negro” y “blanco”
o “0” y “1”.
Una imagen en niveles de gris (monocromática) se puede representar mediante
un arreglo de valores. Sean i y j dos números enteros donde 1 ≤i ≤M y 1 ≤j ≤N.
Adicionalmente f (i, j) denota una función entera tal que 0 ≤f (i, j) ≤W, donde W
denota el nivel de blanco de una imagen en niveles de gris. Un arreglo F se denomina
imagen digital.
f (1, 1) f (1, 2) ... f (1, n)
f (2, 1) f (2, 2) ... f (2, n)
F = . . ... . (7.1)
. . ... .
f (m, 1) f (m, 2) ... f (m, n)
Una dirección (i,j) define una posición en F, llamada pı́xel, pel o elemento de
la imagen. Los elementos de F denotan la intensidad dentro de algunas regiones
rectangulares pequeñas dentro de una imagen como se observa en la Figura 7.10.
Estrictamente hablando, f (i,j) mide la intensidad en un solo punto, pero si la re-
gión rectangular correspondiente es lo suficientemente pequeña, la aproximación
será lo suficientemente precisa para la mayorı́a de los propósitos. El arreglo F con-
tiene un total de M×N elementos y este producto se denomina la resolución espacial
de F. En la Figura 7.11 se muestra una asignación simple de 5 niveles de gris para
una imagen; es posible asignar arbitrariamente intensidades de acuerdo al siguiente
esquema:
Considérese cuantos datos se requieren para representar una imagen en escala de
grises en esta forma. Cada pı́xel requiere el almacenamiento de log2 (1+W) bits. Se
asume que (1+W) es un entero potencia de 2. De no ser ası́, entonces log2 (1+W)
debe redondearse al siguiente entero. Esto se puede representar usando la fun-
214 G. Saldaña-González
Figura 7.10 Imagen digital consistente en un arreglo de M×N, donde el renglón i-ésimo y la
columna j-ésima representan la intensidad igual a f (i,j).
Figura 7.11 Esquema para la asignación de los niveles de gris de una imagen.
ción ceiling (techo) <...>. Entonces una imagen en escala de grises requiere al-
macenamiento de <log2 (1+W)>bits. Dado que hay M×N pı́xeles, el total de los
datos almacenados para la imagen F entera es igual al M×N <log2 (1+W)>bits. Si
M=N≥128 y W≥64, se puede obtener una buena imagen de una cara humana. Mu-
chos de los sistemas de procesamiento de imágenes industriales en uso actualmente
manipulan imágenes en las cuales M=N=512 y W=255. Esto conduce a requerimi-
entos de almacenamiento de 256 Kbytes/imagen. Una imagen binaria requiere un
almacenamiento de M×N bits/imagen.
La señal de una cámara de color se puede representar utilizando tres compo-
nentes:R={r(i,j)}; G={g(i,j)}; B={b(i,j)}, donde R, G y B se definen de forma simi-
lar a F. El vector {r(i,j), g(i,j), b(i,j)} define la intensidad y el color en el punto (i,j)
en la imagen de color.
Las imágenes multiespectrales también se pueden representar usando varias
imágenes monocromáticas. La cantidad total de datos requeridos para codificar una
imagen de color con r componentes es igual a M×N×r <log2 (1+W)>bits, donde
W es simplemente el máximo nivel de la señal en cada canal.
7. Visión por Computadora 215
Estos pasos se pueden dividir en tres categorı́as, basado en la naturaleza del mo-
delo requerido para realizar el análisis: el análisis de la escena en bajo nivel se basa
en las propiedades locales de la imagen, el análisis de la escena a nivel intermedio
utiliza modelos geométricos genéricos y modelos fotométricos y el análisis en al-
to nivel de la imagen se basa en relaciones y en modelos semánticos orientados a
metas.
Este tipo de análisis trata con la integración de las caracterı́sticas locales o pun-
tuales en construcciones globales, por ejemplo: convertir los puntos de un borde en
contornos continuos, dividir la imagen en regiones coherentes, asignar etiquetas a
las entidades detectadas en la escena, evaluación de los parámetros del sensor (de-
terminación de la localización espacial del sensor), derivación de un modelo de las
fuentes de iluminación, etc. Son necesarias diferentes representaciones y técnicas
para el análisis de la escena a nivel intermedio ya que no es posible extender las
técnicas utilizadas en el análisis a bajo nivel a fenómenos globales más complejos
que requieren comprender e interpretar las imágenes.
El análisis a nivel intermedio tiene que ver con la selección de modelos y la
asignación de valores a estos modelos (“instanciación”). Algunas de las tareas a
nivel intermedio son:
División de la imagen. Consiste en separar la imagen en unidades significativas
coherentes.
Enlazado de bordes y dibujo de bosquejos. Organización de los pı́xeles in-
dividuales en segmentos continuos; transformación de las lı́neas primitivas en
representaciones más abstractas.
Recuperación 3D de la geometrı́a de la escena. Utilización de las coacciones
entre bordes y superficies para deducir la geometrı́a tridimensional de la escena
a partir de trazos bidimensionales.
7. Visión por Computadora 217
Los tipos de operaciones que se pueden aplicar a las imágenes digitales para
transformar una imagen de entrada f [m,n] a una imagen de salida b[m,n], se pueden
clasificar en tres categorı́as como se observa en la Tabla 7.2.
Operación Caracteristica
Punto El valor de salida en una coordenada especı́fica
depende sólo del valor de entrada en la misma coordenada
Local (Ventanas) El valor de salida en una coordenada especı́fica
depende de los valores de una vecindad de la
misma coordenada
Global El valor de salida en una coordenada especı́fica
depende de todos los valores en la imagen de entrada
7.7.2. Notación
Figura 7.13 Operadores de imágenes: (a) operadores de punto, (b) operadores locales o de ven-
tanas, (c) operadores globales.
A B C
D E F
G H I
Figura 7.14 Operadores punto a punto. El pı́xel (i,j) en la imagen de entrada tiene intensidad a(i,j)
y se usa para calcular la intensidad del pı́xel c(i,j) en la imagen de salida.
Corrimiento de la intensidad
0 a(i,j)+k < 0
c(i, j) ⇐= a(i, j) + k 0 ≤ a(i,j)+k≤W (7.3)
W W< a(i,j)+k
K es una constante, cuyo valor depende del usuario. Nótese que esta definición
fue cuidadosamente diseñada para mantener c(i,j) dentro del mismo rango que la en-
trada, esto es [0,W]. Esto es un ejemplo de un proceso denominado normalización
de la intensidad. La Normalización es importante porque permite el procesamiento
iterativo por este y otros operadores en una máquina que tiene precisión limitada
para la aritmética (por ejemplo 8 bits).
7. Visión por Computadora 221
Multiplicación de la intensidad
0 a(i,j)·k < 0
c(i, j) ⇐= a(i, j) · k 0 ≤ a(i,j)·k≤W (7.4)
W W< a(i,j)·k
Logaritmo
(
0 a(i,j)=0
c(i, j) ⇐= log(a(i,j)) (7.5)
W[ log(W) ] de otro modo
Esta definición arbitrariamente reemplaza el valor infinito de log(0) por cero, y
de ese modo evita una dificultad de escalamiento del problema.
Antilogaritmo (exponencial)
Resaltado
(
k3 k1≤a(i,j)≤k2
c(i, j) ⇐= (7.9)
a(i, j) de otro modo
Elevación al cuadrado
Figura 7.15 Operador binario punto a punto. Las intensidades de los pı́xeles a(i,j) y b(i,j) en dos
imágenes de entrada se combinan para calcular la intensidad c(i,j) en la imagen de salida.
Es importante notar que c(i,j) depende sólo de a(i,j) y b(i,j). Algunos ejemplos
de operadores binarios son:
Suma
Figura 7.16 Operador local, las intensidades de 9 pı́xeles en una ventana de 3×3 se combinan
para obtener un pı́xel de salida.
Nótese que las intensidades de varios pı́xeles se combinan juntas, para calcular la
intensidad de solo un pı́xel. Entre los operadores locales más simples está aquellos
en los que se utiliza un conjunto de 9 pı́xeles acomodados en un cuadrado de 3×3.
Estos tienen una ecuación caracterı́stica de la forma:
c(i, j) ⇐= g(a(i − 1, j − 1), a(i − 1, j), a(i − 1, j + 1), a(i, j − 1), a(i, j),
E ⇐= k1 · (A · W1 + B · W2 + C · W3 + D · W4 + E · W5+
F · W6 + G · W7 + H · W8 + I · W9) + k2 (7.18)
Donde W1, W2,...W3 son pesos, los cuales pueden ser positivos, negativos o
cero. Los valores para los constantes normalizadas k1 y k2 se darán más tarde. La
matriz ilustrada a continuación se denomina matriz de pesos y es importante porque
determina las propiedades del operador local lineal.
W1 W2 W3
W4 W5 W6
W7 W8 W9
0 V1 0
Q= 0 V2 0
0 V3 0
Se dice que el operador es de tipo Separable. La importancia de esto es que es
posible aplicar dos operadores más simples sucesivamente, con matrices de pesos
P y Q, para obtener el mismo efecto que produce el operador separable.
6. La aplicación sucesiva de los operadores lineales que usan ventanas de 3×3 pı́xe-
les producen el mismo resultado que los operadores locales lineales con ventanas
más grandes. Por ejemplo, aplicando el operador que usa la siguiente matriz de
pesos:
1 1 1
1 1 1
1 1 1
Dos veces sucesivamente resulta en una imagen similar a la obtenida con un ope-
rador de 5×5 con la siguiente matriz de pesos (por simplicidad, se ha ignorado
la normalización):
1 2 3 2 1
2 4 6 4 2
3 6 9 6 3
2 4 6 4 2
1 2 3 2 1
1 3 6 7 6 3 1
3 9 18 21 18 9 3
6 18 36 42 36 18 6
7 21 42 49 42 21 7
6 18 36 42 36 18 6
3 9 18 21 18 9 3
1 3 6 7 6 3 1
226 G. Saldaña-González
Nótese que todos los operadores también son separables. Por lo tanto serı́a posi-
ble reemplazar el último operador de 7×7 mencionado con cuatro operadores
más simples: 3×1, 3×1, 1×3 y 1×3, aplicados en cualquier orden. No siempre
es posible reemplazar un operador de ventana grande con una sucesión de 3×3
operadores. Esto se vuelve obvio cuando se considera, por ejemplo, que un ope-
rador de 7×7 que usa 49 pesos y que tres operadores de 3×3 proporcionan sólo
27 grados de libertad. Sin embargo la separación a menudo es posible, cuando
un operador más grande tiene una matriz de pesos con alguna redundancia, por
ejemplo cuando es simétrica.
7. Para realizar la normalización, se usan los siguientes valores para k1 y k2:
¯ ¯
k1 ← 1/ ∑ ¯Wp,q ¯ (7.19)
p,q
¯ ¯
k2 ← [1 − ∑ Wp,q / ∑ ¯Wp,q ¯]W /2 (7.20)
p,q p,q
8. Un filtro que usa la siguiente matriz de pesos realiza la función de promedio local
sobre una ventana de 11×11:
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1
Detector de esquinas
Filtro de la media
Figura 7.17 Detección de contornos. (a) Imagen original, (b) Detector Roberts, (c) Detector Sobel.
Nótese que mientras el operador Roberts produce contornos más delgados, estos
contornos tienden a romperse en las regiones de alta curvatura. La principal desven-
taja del operador Roberts es su alta sensibilidad al ruido, ya que se utilizan pocos
pı́xeles para calcular el gradiente del contorno. También hay un ligero corrimiento
en la imagen, cuando se usa el detector de contornos Roberts. El detector de con-
tornos Sobel no produce tal corrimiento.
Detector Prewitt
El detector de contornos Prewitt es similar al operador Sobel, pero es más sen-
sible al ruido ya que no posee el mismo suavizado inherente. Este operador usa las
dos ventanas de 3×3 que se observan para determinar el gradiente del contorno:
7. Visión por Computadora 229
-1 -1 -1 -1 0 1
0 0 0 -1 0 1
1 1 1 -1 0 1
P1 P2
c(i, j) ⇐= k1 · (A 0 · W1 + B 0 · W2 + C 0 · W3 + D 0 · W4 + E 0 · W5+
F 0 · W6 + G 0 · W7 + H 0 · W8 + I 0 · W9) + k2 (7.29)
Donde:
A ’=LARGEST(A,B,C,D,E,F,G,H,I)
B ’=SECOND LARGEST(A,B,C,D,E,F,G,H,I)
C ’=THIRD LARGEST(A,B,C,D,E,F,G,H,I)
...
I ’=NINTH LARGEST(A,B,C,D,E,F,G,H,I)
Códigos de dirección
Esta función puede utilizarse para detectar la dirección de la intensidad del gra-
diente. La función DIR CODE se define como:
1 si A≥MAX(B,C,D,E,F,G,H,I)
2 si B≥MAX(A,C,D,E,F,G,H,I)
3 si C≥MAX(A,B,D,E,F,G,H,I)
4 si D≥MAX(A,B,C,E,F,G,H,I)
DIR CODE(A, B,C, D, E, F, G, H, I) ⇐=
5 si E≥MAX(A,B,C,D,F,G,H,I)
6 si F≥MAX(A,B,C,D,E,G,H,I)
7 si G≥MAX(A,B,C,D,E,F,H,I)
8 si H≥MAX(A,B,C,D,E,F,G,H)
(7.30)
Usando esta definición el operador se define como:
Los operadores N-tuple están muy relacionados con los operadores locales y
tienen un gran número de variaciones lineales y no lineales. Los operadores N-
tuple se pueden ver como una versión generalizada de los operadores locales, sus
caracterı́sticas se observan en la Figura 7.18. Para entender los operadores N-tuple,
primero considérese un operador lineal el cual usa una ventana de procesamiento
grande, (r×s pı́xeles) con la mayorı́a de sus pesos iguales a cero. Sólo N de sus
pesos no son cero, donde N<<r×s, este es un filtro N-tuple.
Figura 7.18 Operador N-tuple, donde se observa que a diferencia de los operadores locales, los
pı́xeles que se combinan no forman parte de un conjunto compacto.
7. Visión por Computadora 231
Todos los operadores locales y los filtros N-tuple son susceptibles de producir
efectos peculiares alrededor de los bordes de una imagen. La razón es simplemente
que, para calcular la intensidad de un punto cerca del borde de una imagen, se re-
quiere información acerca de algunos pı́xeles fuera de la imagen, los cuales por
supuesto no están presentes. Para calcular los valores en los pı́xeles de los bordes,
es necesario hacer algunas asunciones, por ejemplo que todos los puntos fuera de la
imagen son negros o tienen el mismo valor que los pı́xeles del borde. Esta estrate-
gia, o alguna otra que se adopte, es perfectamente arbitraria y habrá ocasiones en
que los efectos de borde sean tan pronunciados que no hay nada que se pueda hacer
para removerlos por enmascaramiento. Los efectos de borde son importantes ya que
requieren tomar provisiones cuando se intenta remendar (pegar) varias imágenes de
baja resolución.
Entre las operaciones más comunes para las imágenes binarias se encuentran:
Inversa
OR Exclusiva
Encuentra las diferencias entre las regiones blancas.
E ⇐= A + B +C + D + E + F + G + H + I (7.42)
c(i, j) ⇐= A · B ·C · D · E · F · G · H · I (7.44)
234 G. Saldaña-González
Detector de bordes
Detector de conectividad
Considérese el siguiente patrón:
1 0 1
1 X 1
1 0 1
Si X=1, entonces todos los 1’s tienen 8 conexiones entre si. Alternativamente,
si X=0, entonces no están conectados. En este sentido, el punto marcado con X es
crı́tico para la conectividad. Este es el mismo caso para los siguientes ejemplos:
1 0 0 1 1 0 0 0 1
0 X 1 0 X 0 1 X 0
0 0 0 0 0 1 1 0 1
Sin embargo, estos puntos marcados con X en la parte inferior no son crı́ticos
para la conectividad, ya que puestos a X=0 en lugar de 1 no tiene ningún efecto en
la conectividad de los 1s.
1 1 1 0 1 1 0 1 1
1 X 1 1 X 0 1 X 0
0 0 1 1 1 1 0 1 1
7. Visión por Computadora 235
Número de Euler
El número de Euler se define como el número de componentes conectados
(Blobs) menos el número de huecos en la imagen binaria. El número de Euler repre-
senta un método simple para contar los blobs en una imagen binaria, con tal que no
contengan huecos en ellos. Alternativamente, se puede usar para contar los huecos
en un objeto dado, con tal que no tengan “islas” en ellos. La razón por la cual se uti-
liza esta aproximación para contar blobs, es la simpleza con que puede calcularse el
número de Euler. También es un medio útil para clasificar formas en una imagen. El
número de Euler se puede calcular usando tres operadores locales. Sean los números
N1, N2 y N3 donde N α indica el número de veces que un patrón en el conjunto de
patrones α (α =1, 2 ó 3) ocurre en la imagen de entrada:
0 0 0 0 1 0 0 1
0 1 1 0 0 0 0 0
Conjunto de Patrones 1 (N1)
0 1 1 0
1 0 0 0
Conjunto de Patrones 2 (N2)
1 1 1 1 0 1 1 0
1 0 0 1 1 1 1 1
Conjunto de Patrones 3 (N3)
Llenado de huecos
Considérese una figura de forma semejante a un blob que contiene un hueco
(lago), contra un fondo negro. La aplicación del operador de llenado de huecos
causará que se llenen todos los huecos al poner todos los pı́xeles del hueco en blan-
co. Este operador no altera los bordes externos de la figura.
Etiquetado de regiones
Considérese una imagen que contiene un número de figuras separadas semejantes
a un blob. Un operador de etiquetado de regiones sombreará la imagen de salida de
tal forma que cada blob obtiene un valor de intensidad separado. Se podrı́an som-
brear los blobs de acuerdo al orden en el que son encontrados, durante un escaneo
236 G. Saldaña-González
Figura 7.19 Sombreado de los blobs en una imagen binaria: (a) de acuerdo a las áreas, (b) de
acuerdo al orden en que son encontrados (izquierda a derecha; arriba hacia abajo).
Los blobs pequeños también pueden ser eliminados de una imagen usando este
operador. El etiquetado de regiones puede usarse para contar el número de blobs
diferentes en una imagen. A diferencia del número de Euler, el conteo basado en el
etiquetado de regiones no se afecta por la presencia de huecos.
Figura 7.20 Elemento estructurante. B cabe dentro de la imagen, A no cabe en la imagen dada.
Figura 7.21 Dilatación de una imagen por un elemento estructurante en forma de cruz.
\
E(A, B) = A ª B = Ab (7.49)
b∈B
T
Donde representa la intersección de un conjunto de puntos. La Erosión de la
imagen A por B es el conjunto de todos los puntos para los cuales B trasladado a un
punto x esta contenido en A. Esto consiste en deslizar el elemento estructurante B a
través de A, y donde B esta completamente contenido en A (al colocar el origen del
elemento estructurante en el punto x) después x pertenece a la imagen erosionada
AªB. Por ejemplo en la Figura 7.22 la imagen en la rejilla es erosionada por un
elemento estructurante con forma de cruz contenido en una rejilla de pı́xeles de
3×3.
Figura 7.22 Erosión de una imagen por un elemento estructurante en forma de cruz.
Existe una relación dual entre ciertos operadores morfológicos, tales como la
erosión y la dilatación. Esto significa que el equivalente de tal operación se puede
realizar por su dual en la imagen complemento (negativo) y tomando el comple-
mento del resultado. Aunque duales las operaciones erosión y dilatación no son
inversas la una de la otra. Más bien están relacionadas de acuerdo con las siguientes
relaciones:
Apertura
Es una combinación de las operaciones erosión y dilatación que tiene el efecto
de remover puntos aislados en la imagen que son más pequeños que el elemento es-
tructurante B y aquellas secciones de la imagen A que son más angostas que B. Esto
también puede verse como una operación geométrica de redondeo como se observa
en la Figura 7.23. La apertura de una imagen A por un elemento estructurante B, se
denota por A◦B, y se define como:
(A ª B) ⊕ B (7.52)
Cerradura
Es la operación morfológica dual a la apertura. Esta transformación tiene el efec-
to de llenado de los huecos y bloqueo de valles angostos en la imagen A, donde un
elemento estructurante B (de tamaño similar a los hueco y los valles) es aplicado. El
resultado de esta operación puede verse en la Figura 7.23.
La cerradura de la imagen A por un elemento estructurante B se denota por A•B,
y se define como (A ⊕ B) ª B.
Figura 7.23 Aplicación de un elemento estructurante de 3×3 a una imagen binaria. (a) Imagen
original, (b) Resultado de la operación de Apertura, (c) Resultado de la operación de Cierre.
Esta señal de entrada se representa por la curva delgada y la salida por la curva
gruesa. En este simple ejemplo el elemento estructurante tiene una forma aproxima-
damente parabólica. Para calcular un valor en la señal de salida, el elemento estruc-
turante se empuja hacia arriba, desde abajo en la curva de entrada. La altura de la
cima del elemento estructurante se nota. Este proceso entonces se repite, al deslizar
el elemento estructurante lateralmente. Nótese como este operador particular atenúa
el pico de intensidad pero sigue la señal de entrada de forma precisa en los demás
puntos. Restando la señal de salida de la de entrada se produce un resultado en el
cual el pico de intensidad se enfatiza y todas las demás variaciones se reducen.
El efecto de los operadores morfológicos básicos en imágenes de dos dimen-
siones en escala de grises también se puede explicar en estos términos. Imagı́nese
que la imagen en escala de grises es un paisaje, en el cual cada pı́xel puede verse en
3-D. La dimensión extra de la altura representa el valor en la escala de grises de un
pı́xel. Se generan nuevas imágenes al pasar el elemento estructurante por encima o
por debajo de este paisaje. Figura 7.24.
7. Visión por Computadora 241
Esta operación se define como la erosión de nivel de gris de la imagen seguida por
la dilatación en nivel de gris de la imagen erosionada. El efecto producido cortará los
picos en la topografı́a de niveles de gris al nivel más alto para el cual el elemento
cabe bajo la superficie.
Las funciones que debe realizar un sistema de visión imponen muchas demandas
a los elementos de procesamiento, los cuales el sistema debe manejar satisfactoria-
mente. Entre las demandas principales están requerimientos substanciales para el
poder de procesamiento, gran ancho de banda de E/S, y una latencia bien definida
en el procesamiento.
Un ejemplo tı́pico de estas necesidades computacionales es la convolución, uti-
lizando una máscara programable para aplicaciones tales como filtrado, detección y
reconocimiento de caracterı́sticas. Trabajando con una imagen de 4K×4K a 30 fps,
se requiere que el sistema maneje velocidades de datos por arriba de 490M pı́xe-
les/segundo. Una máscara de convolución tı́pica de 5×5 tiene 25 coeficientes que
deben multiplicarse con pı́xeles almacenados para cada nuevo pı́xel que llega. La
demanda de procesamiento resultante, por tanto, es al menos 12 GFLOPS.
La necesidad de alto desempeño de E/S puede no ser tan obvia, pero puede ser
incluso más crı́tica. Un ejemplo es el histograma, donde para cada pı́xel de la ima-
gen el sistema incrementa un contador correspondiente con el valor del pı́xel para
construir la descripción estadı́stica de la distribución de intensidad de la imagen.
Aunque esta tarea no requiere ningún cálculo, necesita un rápido acceso aleatorio a
los contadores en la memoria para construir un histograma en tiempo real.
7. Visión por Computadora 243
Una latencia bien definida y determinismo son necesarios para obtener resultados
apropiados. Un sistema de visión que inspecciona las partes que se mueven a través
de una banda transportadora debe decidir rápidamente si aceptar o rechazar una
parte. Entre más rápido pueda operar, más rápido operará la lı́nea de manufactura.
Pero el sistema también necesita una baja latencia y el tiempo consistente permitido
por determinismo en el proceso para tomar una decisión. Entre mas tiempo requiera
el sistema para tomar una decisión, la pieza viajará más lejos a lo largo de la banda
entre el punto de inspección y el punto de rechazo de la pieza ası́ que se requiere
una baja latencia para mantener la distancia de viaje de la pieza en lı́mites prácticos.
Los tiempos consistentes permiten al sistema alcanzar sus decisiones en tiempo para
rechazar la parte cuando es necesario.
La estructura de un CLB puede ser tan simple con un transistor o tan compleja
como un procesador. Los CLBs se pueden ordenar en un renglón o más comúnmente
en una matriz. El número de CLBs disponible varı́a de compañı́a en compañı́a. La
red de interconexión sirve como una fábrica implı́cita para proporcionar interco-
nexiones flexibles entre los CLBs para la sı́ntesis de la lógica digital. De los tres
bloques de un FPGA, la red de interconexión tı́picamente ocupa el máximo espacio.
7. Visión por Computadora 245
Figura 7.27 Resultado de los operadores locales aplicados a diversas imágenes. (a) Filtrado, (b)
Operadores morfológicos, Erosión y Dilatación, (c) Pirámide Gaussiana, (e) Multiplicación de
matrices.
que mejoran el funcionamiento de las aplicaciones los FPGA ofrecen la mejor mez-
cla de alto desempeño y flexibilidad en el diseño.
Agradecimientos
Referencias
1. K. R. Castleman, Digital Image Processing, 2nd. Edition, Prentice Hall, New Jersey 1996.
2. K. Compton and S. Hauck, “Reconfigurable Computing: A Survey of Systems and Software”,
ACM Computing Surveys, Vol. 34, no. 2, pp. 171-210, June 2002.
3. E. R. Davies, Machine Vision : Theory, Algorithms, Practicalities, 3er. Edition, Morgan Kauf-
mann, USA 2004.
248 G. Saldaña-González
4. A. DeHon and J. Wawrzynek, “Reconfigurable Computing: What, Why, and Implications for
Design Automation”, Proceedings of the 1999 IEEE Design Automation Conference, pp. 610-
615, June 1999.
5. J. G. Eldredge and B. L. Hutchings, “Run-Time Reconfiguration: A Method for Enhancing
the Functional Density of SRAM-Based FPGAs”, Journal of VLSI Signal Processing, Vol.
12, no. 1, pp. 67-86, January 1996.
6. R. C. González and R. E. Woods, Tratamiento Digital de Imágenes, 1a. Edition, Addison-
Wesley, Delaware 1996.
7. Berthold K. P. Horn, Robot Vision, The MIT Press, 1986.
8. J. R. Parker, Algorithms for Image Processing and Computer Vision, Wiley, USA 1996.
9. G. X. Ritter and J. N. Wilson, Handbook of Computer Vision Algorithms in Image Algebra,
CRC Press, USA 2000.
10. G. Saldaña-González and M. Arias-Estrada, “Compact FPGA-based Systolic Array Architec-
ture Suitable for Vision Systems”, International Journal of High Performance Systems Archi-
tecture, Vol 1, no. 2, pp. 124-132, October 2007.
11. L. G. Shapiro and G. Stockman, Computer Vision, Prentice-Hall, USA 2001.
12. S. E. Umbaugh, Computer Imaging: Digital Image Analysis and Processing, CRC Press, USA
2005.
13. Xilinx, url: http://www.xilinx.com, 2008.
Capı́tulo 8
Nuevas Heurı́sticas Inspiradas en la Naturaleza
para Optimización Numérica
8.1. Introducción
E. Mezura-Montes, O. Cetina-Domı́nguez
Laboratorio Nacional de Informática Avanzada (LANIA A.C.)
Rébsamen 80, Centro, Xalapa, Veracruz, 91000 MEXICO
e-mail: emezura@lania.mx, omarcetina@hotmail.com
B. Hernández-Ocaña
Universidad Juárez Autónoma de Tabasco
División Académica de Informática y Sistemas
Km. 1 Carr. Cunduacán-Jalpa de Méndez
e-mail: betania h o@hotmail.com
249
250 E. Mezura-Montes et al.
que satisfaga de mejor manera una serie de condiciones planteadas. Los problemas
de optimización pueden clasificarse en dos tipos:
1. Problemas de optimización numérica: Se busca un conjunto de valores para
las variables del problema de manera que al sustituirse en la función objetivo
se maximice o minimice el valor de esta función. Un ejemplo de este problema
puede ser el diseño de una pieza mecánica donde se buscan valores óptimos de
sus dimensiones para minimizar su costo de fabricación.
2. Problemas de optimización combinatoria: Se busca encontrar el ordenamiento
de un conjunto de elementos de manera que se maximice o minimice el valor
de la función objetivo. Un ejemplo de este tipo de problemas es el del agente
viajero, que debe recorrer un conjunto de ciudades, pasando por ellas sólo una
vez, de manera que regrese a su punto de salida y se minimice el costo del viaje.
Aquı́ se desea encontrar el orden óptimo de recorrido de las ciudades.
Este capı́tulo se centra en el primer tipo de problemas, conocido como el pro-
blema general de programación no lineal, que se define, sin pérdida de generalidad,
como: Encontrar x que minimiza
f (x) (8.1)
sujeta a:
gi (x) ≤ 0, i = 1, . . . , m (8.2)
h j (x) = 0, j = 1, . . . , p (8.3)
n
]T
donde x ∈ IR es un vector de n variables de decisión x = [x1 , x2 , . . . , xn y cada
xi , i = 1, ..., n está acotado por lı́mites inferiores y superiores Li ≤ xi ≤ Ui , los
cuales definen el espacio de búsqueda S , F es el conjunto de todas las soluciones
que satisfacen las restricciones del problema y se le llama zona factible, siendo claro
que F ∈ S , m es el número de restricciones de desigualdad y p es el número de
restricciones de igualdad. Las restricciones y la función objetivo pueden ser lineales
y no lineales.
En este problema se pretende encontrar una solución que minimice una cierta
medida de calidad, conocida como función objetivo (8.1) que además debe satis-
facer un conjunto de restricciones asociadas al problema (8.2)-(8.3). Detallando el
ejemplo propuesto anteriormente sobre el diseño de la pieza mecánica, se quiere
encontrar su costo mı́nimo de construcción, el cual se modela mediante una fun-
ción (no lineal tal vez) y que depende de variables como el largo y el ancho de la
pieza, ası́ como el grosor del material. Además, pueden existir restricciones de es-
pacio, de manera que la pieza no puede ser ni muy grande ni muy chica, pues debe
incorporarse a un mecanismo más grande.
8. Nuevas Heurı́sticas para Optimización Numérica 251
Seleccion
Reemplazo Soluciones
seleccionadas
Nuevas Operadores
Soluciones de variacion
se basará en el modelo propuesto por Karaboga [13] que resuelve problemas de op-
timización numérica mediante dos comportamientos: El reclutamiento de abejas en
una fuente de alimento y el abandono de una fuente.
El modelo biológico de recolección de alimento en abejas melı́feras consta de los
siguientes elementos mı́nimos [13]:
Fuentes de alimento: El valor de una fuente de alimento depende de muchos
factores, como lo son la cercanı́a a la colmena, la concentración de alimento y la
facilidad para extraerlo. Por simplicidad es posible representar la rentabilidad de
una fuente en un solo valor numérico.
Recolectoras empleadas: Estas abejas están asociadas con una fuente particular
de alimento la cual están explotando, es decir, están empleadas para con ella. Las
abejas empleadas comparten la información sobre su fuente de alimento, como
la ubicación y rentabilidad hacia las abejas recolectoras.
Recolectoras desempleadas: Están constantemente buscando una fuente de ali-
mento a la cual explotar. Hay dos tipos de recolectora desempleada: la explo-
radora, que busca en las cercanı́as de la colmena por nuevas fuentes de alimento
y la observadora, que espera en la colmena y elige una fuente de alimento con
base en la información compartida por una recolectora empleada.
La información sobre las fuentes de alimento es compartida por las recolectoras
empleadas por medio de una danza cuya duración indica la rentabilidad de la fuente,
el ángulo con respecto al sol indica la dirección de la fuente y el número de movi-
mientos en zig-zag durante la danza indica la distancia. Dado que las danzas de las
fuentes más rentables son de mayor duración, existe una mayor probabilidad de ser
observada por recolectoras desempleadas, aumentando la probabilidad de que una
recolectora elija dicha fuente de alimento.
Cuando una fuente de alimento se agota, la abeja o abejas empleadas en ellas se
convierten en abejas desempleadas y tendrán que decidir entre convertirse en abejas
exploradoras o regresar a la colmena como abejas observadoras.
En el algoritmo ABC la colonia de abejas artificiales se compone también de 3
grupos de abejas: abejas empleadas, abejas observadoras y abejas exploradoras. El
número de abejas empleadas es usualmente igual al número de fuentes de alimen-
to y se asignará una abeja empleada a cada una de las fuentes. Al llegar a dicha
fuente, la abeja calculará una nueva solución (volará hacia otra fuente de alimento
cercana) a partir de ésta y conservará la mejor solución. El número de abejas obser-
vadoras es también usualmente igual al de abejas empleadas y se asignarán a una
fuente de alimento con base en la aptitud de éstas, al igual que las abejas empleadas,
calculará una nueva solución a partir de su fuente de alimento. Cuando una fuente
no mejora después de un cierto número de iteraciones se abandona, siendo reem-
plazado por aquella encontrada por una abeja exploradora, lo cual implica calcular
una nueva solución aleatoriamente. Una de las ventajas de este algoritmo es el bajo
número de parámetros que requiere como puede verse en la Tabla 8.1.
256 E. Mezura-Montes et al.
Cada una de las variables que forma parte de la solución está asociada a un rango
(Li ≤ xi ≤ Ui ), el cual se debe de tomar en cuenta al momento de generar de ma-
nera aleatoria, con una distribución uniforme, las soluciones (fuentes de alimento)
iniciales.
En el ABC, las abejas son vistas como operadores de variación, pues cuando una
de ellas llega a una fuente de alimento, calcula una nueva solución candidata vi,g
utilizando la fórmula dada en (8.5). En donde xi,g representa la solución en la que
la abeja se encuentra en ese momento, xk,g es una fuente de alimento aleatoria (y
distinta de xi,g ), g es el número de ciclo actual del programa y φ es un número real
aleatorio en el intervalo [-1,1].
¡ ¢
vi,g = xi,g + φ · xi,g − xk,g (8.5)
En la Figura 8.3 puede verse un ejemplo gráfico de este cálculo. En primer lugar
se observa en la Figura 8.3(a) el vector que es generado mediante la diferencia entre
xi y xk . Posteriormente se genera la solución candidata utilizando (8.5) en la Figura
8.3(b). Se nota que v y v0 se generaron utilizando el mismo valor de φ pero con signo
opuesto. Ésto es posible porque el intervalo de φ lo permite ([-1,1]).
Y Y
xk xk
xi− xk v
xi
xi
v’
X X
(a) (b)
Figura 8.3 (a) El vector generado por la diferencia entre xi y xk . (b) Dos posibles soluciones
candidatas generadas utilizando (8.5) y el mismo valor φ con signo opuesto.
1 Begin
2 Inicializar la población de soluciones xi,0 , i = 1, . . . , SN
3 Evaluar la población
4 g=1
5 Repeat
6 Producir nuevas soluciones vi,g para las abejas empleadas
utilizando (8.5) y evaluarlas
7 Conservar la mejor solución entre la actual y la candidata
8 Seleccionar las soluciones que serán visitadas por una abeja
observadora según su aptitud
9 Producir nuevas soluciones vi,g para las abejas observadoras
utilizando (8.5) y evaluarlas
10 Conservar la mejor solución entre la actual y la candidata
11 Determinar si existe una fuente abandonada y reemplazarla
utilizando una abeja exploradora
12 Memorizar la mejor solución encontrada hasta este momento
13 g = g+1
14 Until g = MCN
15 End
con una nueva encontrada por una abeja exploradora, que en este caso es una nueva
solución calculada aleatoriamente dentro del espacio de búsqueda.
Finalmente la mejor solución encontrada durante el ciclo se compara con la mejor
solución en memoria y si tiene una mejor aptitud la reemplazará. Se repite el ciclo
MCN veces.
−5 ≤ x1 ≤ 5
−5 ≤ x2 ≤ 5 (8.7)
x1 , x2 ∈ R
Los valores de los parámetros del ABC para este ejemplo se muestran en la Tabla
8.2.
Parámetro Valor
SN 5
MCN 2
limit 10
Fuente x1 x2 f (xi )
1 (x1,0 ) 0.5 0.3 0.34
2 (x2,0 ) -0.2 0.1 0.05
3 (x3,0 ) -0.4 -0.7 0.65
4 (x4,0 ) 0.2 0.05 0.0425
5 (x5,0 ) 0.1 0.9 0.82
Tabla 8.3 Población inicial, el valor de la función para cada individuo se encuentra en la cuarta
columna.
260 E. Mezura-Montes et al.
Empleada k φ1 φ2
1 3 -0.1 0.2
2 1 0.2 0.5
3 5 0.7 -0.3
4 2 -0.5 0.7
5 1 0.3 -0.1
Fuente x1 x2 f (x)
1 (x1,1a ) 0.5 0.3 0.34
2 (x2,1a ) -0.2 0.1 0.05
3 (x3,1a ) -0.75 -0.22 0.6109
4 (x4,1a ) 0 0.015 0.000225
5 (x5,1a ) -0.02 0.84 0.706
El siguiente paso es enviar las abejas observadoras a las mejores fuentes de ali-
mento (con base en su aptitud), utilizando algún mecanismo de selección propor-
8. Nuevas Heurı́sticas para Optimización Numérica 261
cional. Por motivos de simplicidad, para este ejemplo se determinó enviar tres abe-
jas observadoras a la mejor solución (fuente de alimento 4) y las dos restantes a la
segunda mejor (fuente 2). Los valores aleatorios para calcular las soluciones candi-
datas pueden verse en la Tabla 8.6.
Fuente x1 x2 f (x)
1 (x1,1 ) 0.5 0.3 0.34
2 (x2,1 ) -0.088 0.06 0.011344
3 (x3,1 ) -0.75 -0.22 0.6109
4 (x4,1 ) 0 0.015 0.000225
5 (x5,1 ) -0.02 0.84 0.706
Empleada k φ1 φ2
1 2 0.1 -0.9
2 4 -0.3 0.25
3 5 0.4 -0.27
4 2 -0.05 0.2
5 3 0.3 -0.1
Fuente x1 x2 f (x)
1 (x1,2a ) 0.5588 0.084 0.31931344
2 (x2,2a ) -0.0616 0.07125 0.008871123
3 (x3,2a ) -0.75 -0.22 0.6109
4 (x4,2a ) -0.0044 0.006 5.536E-05
5 (x5,2a ) 0.199 0.734 0.578357
Para las abejas observadoras nuevamente se utilizan tres abejas en la mejor solu-
ción (fuente 4) y dos abejas en la segunda mejor (fuente 2) con la intención de re-
presentar una selección proporcional con base en la aptitud de la fuente de alimento.
Los valores aleatorios utilizados se muestran en la Tabla 8.10.
Fuente x1 x2 f (x)
1 (x1,2 ) 0.5588 0.084 0.31931344
2 (x2,2 ) -0.032428 0.0843 0.008158065
3 (x3,2 ) -0.75 -0.22 0.6109
4 (x4,2 ) -0.0044 0.006 5.536E-05
5 (x5,2 ) 0.199 0.734 0.578357
El BFOA original, fue propuesto por K.M. Passino [14], el cual simula el pro-
ceso completo de forrajeo de la bacteria E. Coli: Chemotaxis o paso quimiotáctico
(movimientos de giro y nado), reproducción y eliminación-dispersión. El modelado
de la búsqueda de alimento (concentración de nutrientes óptima) que realizan las
bacterias contempla las siguientes etapas:
las bacterias y finalmente el recorrido que tienen que hacer para encontrar el punto
óptimo de concentración de nutrientes, entre otros.
Debido a que el algoritmo de Passino (BFOA) requiere de muchos parámetros y
su costo computacional es alto, la versión modificada del BFOA, propuesta en [15]
y llamada MBFOA es la utilizada en este capı́tulo.
La representación de las soluciones, al igual que para el caso del ABC, se lleva
a cabo mediante vectores de números reales llamados bacterias y representados en
(8.8) por la letra θ , donde θ i ( j, G) representa la bacteria i, en su ciclo quimiotácti-
co j, en la generación G y Sb es el número total de bacterias en la población. La
representación gráfica es similar a la de las soluciones del ABC en la Figura 8.2.
θ i ( j, G), ∀i = 1, . . . , Sb (8.8)
Para el caso de las bacterias, el proceso de selección se realiza dentro del paso
quimiotáctico, donde cada bacteria se moverá de manera constante en una misma
dirección, generada de manera aleatoria, hacia una nueva posición siempre y cuando
sea de mejor calidad que la previa. A este proceso dentro del ciclo quimiotáctico se
le conoce como nado. Si al pasar a una nueva posición y ésta es de menor calidad,
la bacteria realizará un giro para buscar una nueva dirección de búsqueda siempre
buscando mejorar la calidad de su posición.
∆ (i)
φ (i) = p (8.9)
∆ (i)T ∆ (i)
donde ∆ (i) es un vector del mismo número de dimensiones que las variables del
problema a resolver y generado aleatoriamente con sus elementos dentro del inter-
valo: [−1, 1]. A este paso se le conoce como giro. Después de ello, cada bacteria
θ i ( j, G) modifica su posición mediante el nado como se indica en (8.10).
8. Nuevas Heurı́sticas para Optimización Numérica 265
Begin
Inicializar la población de soluciones (bacterias) θ i (0, 0) ∀i, i = 1, . . . , Sb
Evaluar la población de bacterias f (θ i (0, 0)) ∀i, i = 1, . . . , Sb
For G=1 to GMAX Do
For i=1 to Sb Do
For j=1 to Nc Do
Realizar el paso quimiotáctico (giro-nado y swarming)
para la bacteria θ i ( j, G) usando (8.10) y (8.12)
End For
End For
Realizar la reproducción eliminando las Sr peores
bacterias y duplicando cada una de las (Sb − Sr ) mejores bacterias
Eliminar la peor bacteria θ w ( j, G) en la población actual
End For
End
Parámetro Valor
Sb 3
Nc 3
GMAX 4
Sr 2
R 2.1E-2
β 0.44
Tabla 8.14 Población inicial, el valor de la función para cada bacteria se encuentra en la cuarta
columna.
A continuación, se tiene que evaluar cada una de las bacterias en la función ob-
jetivo, estos valores se muestran en la cuarta columna de la Tabla 8.14.
El siguiente paso consiste en inicializar el contador del número de generaciones
(G), el contador del número de bacterias Sb y el contador del número de iteraciones
del ciclo quimiotáctico Nc con el valor de 1. En el primer ciclo quimiotáctico, para
la bacteria 1 θ 1 ( j, G) se aplica el operador de nado-giro, (8.9) y (8.10), recordando
que a la mitad del ciclo quimiotáctico y al final del mismo se aplicará el swarming
(8.12). Supóngase que los valores que se emplean para la dirección aleatoria son
(8.9):
φ1 (1) = −0.289841
φ2 (1) = −0.029217
Se calcula la nueva posición de la bacteria 1 usando el operador de nado (8.9) de
la siguiente manera:
268 E. Mezura-Montes et al.
Nc = 1
θ11 (1, 1) = −1.728921 + (0.1485 ∗ −0.289841) = −1.77190
θ21 (1, 1) = 0.564320 + (0.1485 ∗ −0.029217) = 0.55998
A continuación se evalúa la calidad de la nueva posición de la bacteria 1, quedan-
do de la siguiente manera:
Nc = 1
θ12 (1, 1) = (1.595993) + (0.1485 ∗ 0.825212) = 1.718537
θ22 (1, 1) = (1.107934) + (0.1485 ∗ −0.564822) = 1.024058
f (θ 2 (1, 1)) = (1.718537)2 + (1.024058)2 = 4.002064
Nc = 2
θ12 (2, 1) = (1.595993) + (0.44 ∗ (1.082393 − 1.595993)) = 1.370009
θ22 (2, 1) = (1.107934) + (0.44 ∗ (1.011439 − (1.107934))) = 1.065476
f (θ 2 (2, 1)) = (1.370009)2 + (1.065476)2 = 3.012164
Nc = 3
θ12 (3, 1) = (1.370009) + (0.44 ∗ (1.082393 − 1.370009)) = 1.243458
θ22 (3, 1) = (1.065476) + (0.44 ∗ (1.011439 − 1.065476)) = 1.041700
f (θ 2 (3, 1)) = (1.243458)2 + (1.041700)2 = 2.631326
Los cálculos para la tercera bacteria se presentan a continuación:
Nc = 1
θ13 (1, 1) = (1.082393) + (0.1485 ∗ −0.519574) = 1.005236
θ23 (1, 1) = (1.011439) + (0.1485 ∗ 0.854426) = 1.138321
f (θ 3 (1, 1)) = (1.005236)2 + (1.138321)2 = 2.306275
Nc = 2
θ13 (2, 1) = (1.082393) + (0.44 ∗ (1.082393 − 1.082393)) = 1.082393
θ23 (2, 1) = (1.011439) + (0.44 ∗ (1.011439 − 1.011439)) = 1.011439
f (θ 3 (2, 1)) = (1.082393)2 + (1.011439)2 = 2.194583
Nc = 3
θ13 (3, 1) = (1.082393) + (0.44 ∗ (1.082393 − 1.082393)) = 1.082393
θ23 (3, 1) = (1.011439) + (0.44 ∗ (1.011439 − 1.011439)) = 1.011439
f (θ 3 (3, 1)) = (1.082393)2 + (1.011439)2 = 2.194583
Los valores para cada una de las bacterias después de haber pasado el proceso
quimiotáctico es presentado a continuación en la Tabla 8.15
La solución a reportar será la bacteria 1 y 2 que tienen los mismo valores, con
x1 = −0.042711, x2 = −0.040718 y f (x) = 0.003482.
8. Nuevas Heurı́sticas para Optimización Numérica 271
8.5. Conclusiones
Agradecimientos
Referencias
3. Z. Michalewicz and D.B. Fogel. How to Solve It: Modern Heuristics. Springer-Verlag, second
edition, 2004.
4. E. Mezura-Montes, C.A. Coello Coello, and R. Landa-Becerra. Engineering Optimization
Using a Simple Evolutionary Algorithm. In Proceedings of the Fiftheenth International Con-
ference on Tools with Artificial Intelligence (ICTAI’2003), pages 149–156, Los Alamitos, CA,
November 2003. Sacramento, California, IEEE Computer Society.
5. A. Angantyr and J.O. Aidanpää. A Pareto-Based Genetic Algorithm Search Approach to
Handle Damped Natural Frequency Constraints in Turbo Generator Rotor System Design.
Journal of Engineering for Gas Turbines and Power, 126(3):619–625, July 2004.
6. J. Alvarez-Gallegos, C.A. Cruz-Villar, and E.A. Portilla-Flores. Evolutionary Dynamic Op-
timization of a Continuously Variable Transmission for Mechanical Efficiency Maximiza-
tion. In Alexander Gelbukh, Álvaro de Albornoz, and Hugo Terashima-Marı́n, editors, MICAI
2005: Advances in Artificial Intelligence, pages 1093–1102, Monterrey, México, November
2005. Springer. Lecture Notes in Artificial Intelligence Vol. 3789,.
7. S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi. Optimization by Simulated Annealing.
Science, 220(4598):671–680, May 1983.
8. F.W. Glover and M. Laguna. Tabu Search. Kluwer Academic Publishers, London, UK, 1997.
9. A.P. Engelbrecht. Fundamentals of Computational Swarm Ingelligence. John Wiley & Sons,
Sussex, England, 2005.
10. A. E. Eiben and J. E. Smith. Introduction to Evolutionary Computing. Springer, first edition,
2003.
11. D.B. Fogel. An introduction to simulated evolutionary optimization. Neural Networks, IEEE
Transactions on, 5:3–14, 1994.
12. A. Baykasoglu, L. Ozbakir, and P. Tapkan. Artificial bee colony algorithm and its application
to generalized assignment problem. In Felix T.S. Chan and Manoj Kumar Tiwari, editors,
Swarm Intelligence: Focus on Ant and Particle Swarm Optimization, pages 113–144. Itech
Education and Pub., Vienna, Austria, 2007. ISBN 978-3-902613-09-7.
13. D. Karaboga and B. Basturk. A powerful and efficient algorithm for numerical function opti-
mization: Artificial bee colony (ABC) algorithm. Journal of Global Optimization, 39(3):459–
471, 2007.
14. K.M. Passino. Biomimicry of bacterial foraging for distributed optimization andcontrol. IEEE
Control Systems Magazine, 22(3):52–67, 2002.
15. E. Mezura-Montes and B. Hernández-Ocaña. Bacterial foraging for engineering design prob-
lems: Preliminary results. In Arturo Hernández-Aguirre et al., editor, Fourth Mexican Con-
ference on Evolutionary Computation (COMCEV’2008), pages 33–38, Guanajuato, Mexico,
October 2008. CIMAT.