Aspectos Pr
´
acticos del Problema de Identificaci
´
on
de la Din
´
amica de Robots Industriales
Practical Aspects of the Dynamic Model Identification of Industrial Robots
Pablo S. Gonzalez
Facultad de Ingenier
´
ıa, Universidad de Buenos Aires
Paseo Col
´
on 850 (1063) Buenos Aires, Argentina
pgonzal@fi.uba.ar
Abstract—The practical aspects of the dynamical parameter
identification method of an ABB IRB140 robot (off-line
calibration procedure) are presented, such as identifiability,
reduction to base parameters, inclusion of the measurements
noise model and their influence on estimation uncertainty,
scaling, and prior knowledge incorporation. A remarkable
result is the calculation of the trajectories that provide a
set of persistently excitatory signals of the dynamics in
order to improve the estimation. The signals preprocessing is
another relevant aspect since, to reduce the effect of noise,
the excitation is limited in bandwidth and the derivatives of
the joint positions are calculated in the transformed space.
Keywords: robot; parametric identification; robot dynamics;
trajectory optimization.
Resumen— Se presentan los aspectos pr
´
acticos del m
´
etodo
de identificaci
´
on de los par
´
ametros din
´
amicos de un robot
ABB IRB140 (procedimiento de calibraci
´
on off-line), tales
como la identificabilidad, la reducci
´
on a par
´
ametros base,
la inclusi
´
on del modelo del ruido en las mediciones y
su influencia en la incertidumbre de la estimaci
´
on, el
escalamiento, y la incorporaci
´
on del conocimiento previo.
Un resultado destacable es el c
´
alculo de las trayectorias que
brindan un conjunto de se
˜
nales persistentemente excitatorias
de las din
´
amicas para mejorar la estimaci
´
on, pasando la
planificaci
´
on del ensayo a ser una variable m
´
as del proceso
de identificaci
´
on. El tratamiento de las se
˜
nales es otro aspecto
relevante pues para reducir el efecto del ruido se limita en
frecuencia a la excitaci
´
on y las derivadas de las posiciones
articulares se calculan en el espacio transformado.
Palabras clave: robot; identificaci
´
on param
´
etrica; din
´
amica;
optimizaci
´
on de trayectorias.
I. INTRODUCCI
´
ON.
La din
´
amica de los robots articulados, es decir aquellos
que tienen juntas de revoluci
´
on que articulan los eslabones,
presenta un comportamiento fuertemente no lineal. Los
controladores m
´
as simples (por ejemplo PID) permiten
estabilizar la planta y dar una respuesta satisfactoria en
cuanto las prestaciones requeridas sean exigentes.
Para mejorar la performance es necesario considerar el
comportamiento no lineal y el acoplamiento en la din
´
amica
de los ejes. Por lo tanto una estrategia basada en N
controladores SISO, tantos como ejes se tengan, no es
suficiente. La t
´
ecnica de control por torque computado [1],
permite desacoplar y linealizar en forma exacta el modelo de
la planta de manera que puede ser expresado por la siguiente
ecuaci
´
on,
v =
¨
q
donde la entrada es v y
¨
q es el vector de aceleraciones
articulares. Luego las salidas del sistema son las posiciones
de los ejes q y las velocidades articulares
˙
q.
Entonces el problema de control se reduce a reproducir
una trayectoria de referencia en una planta compuesta de N
doble integradores, para lo cual las t
´
ecnicas de control lineal
dan sobradas alternativas de elevado desempe
˜
no.
Este escenario se consigue siempre que la cancelaci
´
on de
las no-linealidades sea exacta, lo que implica contar con un
modelo de la din
´
amica que describa de forma exacta las
relaciones entre las variables.
El planteo del modelo de la din
´
amica expresa estas
relaciones en funci
´
on de una serie de par
´
ametros, tales
como las masas de los eslabones o los momentos y
productos de inercia entre otros. La estructura del modelo
obedece a los efectos considerados; por ejemplo, un modelo
que contemple las elasticidades de las juntas tendr
´
a una
estructura diferente a aquel que no las incluye. Por lo tanto
para tener un modelo exacto
1
se necesita no solo incluir
todos los efectos sino tambi
´
en establecer los valores de los
par
´
ametros para tal fin.
Es claro que no es posible acceder a un conocimiento
absoluto y lo mejor que se puede hacer es acercarse en la
menor medida posible, hecho que podr
´
ıamos traducir en el
conocimiento con menor incertidumbre de los par
´
ametros
dada una relaci
´
on estructural que incluya los efectos
preponderantes.
El control robusto y el control adaptativo proponen dos
caminos para lidiar con este problema. En el primero de
los casos se plantea un control que estabilice la familia de
modelos descriptos por el conjunto de valores que pueden
tomar los par
´
ametros, dada la incertidumbre. En el segundo
de los casos se busca identificarlos durante el movimiento
adecuando el comportamiento del controlador.
La incertidumbre en los par
´
ametros no solamente se
explica en los t
´
erminos tratados. Tambi
´
en ocurre en
situaciones tales como cuando el robot manipula un objeto
desconocido; como el robot es una m
´
aquina que puede
realizar tareas que no est
´
en completamente determinadas
a-priori, es natural encontrarse en esta situaci
´
on. En este
caso es importante partir de un valor inicial de los
par
´
ametros din
´
amicos para ajustar la estimaci
´
on durante
la ejecuci
´
on de la tarea. As
´
ı se contempla una estrategia
1
El sentido de esta expresi
´
on es referirnos a una construcci
´
on abstracta
que intenta describir el comportamiento de un sistema de la manera m
´
as
fiel posible. Tener presente la diferencia entre modelo y realidad
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
83
Recibido: 01/11/21; Aceptado: 01/12/21
Creative Commons License - Attribution-NonCommercial-
NoDerivatives 4.0 International (CC BY-NC-ND 4.0)
https://doi.org/10.37537/rev.elektron.5.2.138.2021
Original Article
off-line de caracterizaci
´
on de par
´
ametros propios del brazo
complementado con un algoritmo on-line que identifique la
carga durante el movimiento.
El marco te
´
orico en temas de identificaci
´
on est
´
a nutrido
por el trabajo de muchos grupos de investigaci
´
on. En
particular, en el Laboratorio de Rob
´
otica de FIUBA
se cuentan como antecedentes dos tesis de grado
correspondientes a Diego Zeida [2] y Alejandro Oubi
˜
na [3],
la tesis de doctorado de Mauricio Anigstein [4], y varios
art
´
ıculos con referato. Los primeros trabajos publicados en
el tema datan de mediados de la d
´
ecada de 1980. Hasta el
d
´
ıa de la fecha se encuentra como un
´
area activa, donde el
inter
´
es reciente se enfoca en la mejora de los m
´
etodos de
estimaci
´
on y el tratamiento de las se
˜
nales medidas.
A. Organizaci
´
on del art
´
ıculo.
En la secci
´
on II se repasa el m
´
etodo de modelizaci
´
on
de la din
´
amica, que establece las relaciones entre variables
y par
´
ametros. Luego en III se presenta el algoritmo de
identificaci
´
on, el c
´
alculo del sesgo y la incertidumbre de
la estimaci
´
on. Tambi
´
en se eval
´
uan distintas t
´
ecnicas para
mejorar los resultados consistentes en el escalamiento, la
reducci
´
on del modelo y la incorporaci
´
on de conocimiento
previo. En la secci
´
on IV-B se aborda el c
´
alculo de la
trayectoria
´
optima. Los resultados se presentan en VI
seguidamente de las conclusiones.
II. MODELO DE LA DIN
´
AMICA DEL ROBOT
Uno de los m
´
etodos para estudiar la din
´
amica de los
manipuladores es mediante la resoluci
´
on de la ecuaci
´
on de
Lagrange (ec. 1), donde el Lagrangiano L se calcula como
la resta entre la energ
´
ıa cin
´
etica T (q,
˙
q) y la potencial
U(q), siendo las variables τ
e
s
y q
s
las fuerza y posici
´
on
generalizada sobre el eje s.
τ
e
s
=
d
dt
∂L
∂ ˙q
s
−
∂L
∂q
s
(1)
El m
´
etodo de Lagrange permite obtener anal
´
ıticamente un
modelo en forma cerrada donde el principal desaf
´
ıo consiste
en calcular T y U.
La energ
´
ıa cin
´
etica puede expresarse entonces como una
forma cuadr
´
atica seg
´
un la ec. 2.
T =
1
2
˙
q
t
M
˙
q (2)
La matriz de inercia o de energ
´
ıa cin
´
etica M se calcula
seg
´
un la ec. 3, donde A
i
0
es la matriz homog
´
enea de
rototraslaci
´
on que relaciona la terna asociada al eslab
´
on i
con la terna inercial
´
o 0, y J
i
es la matriz de pseudoinercia
del eslab
´
on i que incluye como par
´
ametros a los momentos
de inercia (I
oxxi
, I
oyyi
, I
ozzi
), los productos de inercia
(I
oxyi
, I
oxzi
, I
oyzi
), la posici
´
on del centro de masas (r
i
g
=
x
Gi
y
Gi
z
Gi
t
) y su masa (m
i
), resultando en un total
de 10.
m
sj
=
N
X
i=max(s,j)
Tr
"
∂A
i
0
∂q
s
J
i
∂A
i
0
∂q
j
t
#
(3)
La energ
´
ıa potencial, desestimando el efecto el
´
astico de
articulaciones y considerando a los eslabones como cuerpos
r
´
ıgidos, se expresa seg
´
un la ec. 4, donde g es la aceleraci
´
on
de la gravedad y r
i
g
es la posici
´
on del centro de masas del
eslab
´
on i respecto de su terna.
U = m
i
g
t
A
i
0
r
i
g
(4)
Reemplazando las ecs. 2 y 4 en 1 se obtiene en forma
general el IDIM, o modelo de la din
´
amica inversa
2
del
manipulador, seg
´
un la ec. 5, donde se indica la dependencia
de cada uno de los elementos M , C y G que resultan
mayormente en funciones no lineales de las posiciones y
velocidades generalizadas.
τ
e
= M(q)
¨
q + C(q,
˙
q)
˙
q + G(q) (5)
La matriz C se obtiene derivando los componentes de M
seg
´
un se expresa en la ec. 6 y los elementos del vector de
torques producidos por peso propio G se calculan con la
ec. 7.
c
sj
=
N
X
k=1
1
2
∂m
sj
∂q
k
+
∂m
sk
∂q
j
−
∂m
jk
∂q
s
˙q
k
(6)
g
s
= −
N
X
i=s
m
i
g
T
∂A
i
0
∂q
s
r
i
g
(7)
Luego al considerar los motores, las transmisiones y los
reductores, se suman los momentos de inercia de los rotores,
de los engranajes, y sus masas. Adem
´
as por el hecho de
que la unidad mec
´
anica debe transportar los motores, la
masa y la ubicaci
´
on de los estatores y de los reductores
incide sobre los par
´
ametros din
´
amicos de los eslabones.
Finalmente la p
´
erdida producida en los rodamientos y en
los componentes del reductor pueden modelizarse como
rozamiento seco y viscoso. Todos estos efectos se deber
´
ıan
incorporar al planteo, previo a la resoluci
´
on de la ec. 1. Sin
embargo, considerando que los ejes est
´
an actuados a trav
´
es
de reducciones de orden N
r
∼ 100 se puede desestimar
el d
´
ebil acoplamiento din
´
amico entre el rotor del motor y
el eslab
´
on
3
, simplificando el planteo al resolver la ec. de
Lagrange para el mecanismo del brazo sumando luego el
efecto de los motores y las reducciones [5] . Se considera
adem
´
as que cada motor act
´
ua solo sobre un eje, hecho que
produce que la matriz N
r
resulte diagonal.
El modelo din
´
amico del actuador puede expresarse
entonces seg
´
un la ec. 8, donde τ es el vector de torques
motor y los nuevos par
´
ametros que se agregan son el
momento de inercia equivalente (J
m
), el coeficiente de
rozamiento viscoso (B
m
) y el de rozamiento seco o de
Coulomb (F
c
).
N
r
τ = J
m
N
2
r
¨
q + B
m
N
2
r
˙
q + F
c
tanh(
˙
q) + τ
e
(8)
As
´
ı las ecs. 5 y 8 definen el modelo de la din
´
amica del
brazo con las simplificaciones propuestas. Ambas pueden
2
El modelo directo permite implementar el simulador, y se obtiene
f
´
acilmente pues M siempre tiene inversa
3
Esta simplificaci
´
on es equivalente a considerar que la velocidad angular
del rotor del motor se debe
´
unicamente a su propio giro, pues el efecto
producido sobre
´
el por el movimiento del mecanismo es al menos 2 ordenes
de magnitud menor
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
84
http://elektron.fi.uba.ar
resumirse en la ec. 9, donde el conjunto de par
´
ametros tiene
un cardinal de 13 por eje.
τ = M(q)
¨
q + h(q,
˙
q) (9)
La ec. 9 vincula los torques activos en los ejes con
las variables de estado x
t
=
q
t
˙
q
t
y sus derivadas.
Observando a continuaci
´
on las ecs. 3, 6, 7 y 8, se concluye
que las operaciones sobre los par
´
ametros din
´
amicos son
todas lineales si se considera las posiciones de los centros de
masa multiplicadas por la masa (m
i
r
i
g
). Entonces el IDIM
se puede reescribir seg
´
un la ec. 10, donde la matriz φ es
conocida como el regresor de la din
´
amica y p almacena los
N × 13 par
´
ametros enunciados.
τ = φ(q,
˙
q,
¨
q)p (10)
Como el m
´
etodo de c
´
alculo presentado es general y
sistematizado, se implement
´
o un algoritmo que genera en
forma autom
´
atica el modelo din
´
amico para simular en
Matlab y la funci
´
on que implementa el regresor, a partir
de los par
´
ametros cinem
´
aticos y din
´
amicos expresados en
forma simb
´
olica.
III. IDENTIFICACI
´
ON PARAM
´
ETRICA POR M
´
INIMOS
CUADRADOS
La identificaci
´
on del IDIM puede tener como objetivo
hallar un modelo estructural obteniendo as
´
ı los par
´
ametros
f
´
ısicos con la menor incertidumbre posible, o bien conseguir
un modelo para la predicci
´
on donde el foco est
´
a puesto en el
error de la salida. La identificaci
´
on se plantea en este trabajo
como un problema de ajuste o estimaci
´
on de par
´
ametros
donde la meta es la minimizaci
´
on del error de predicci
´
on.
Una t
´
ecnica habitual consiste en la aplicaci
´
on de
cuadrados m
´
ınimos. El planteo en este caso lleva a
minimizar el error de predicci
´
on del IDIM seg
´
un un ensayo
donde se adquieren K muestras de los torques aplicados
τ
k
y sus posiciones y derivadas (q
k
,
˙
q
k
,
¨
q
k
). El regresor φ
k
se calcula para cada una de las muestras y se plantea un
problema lineal con KN ecuaciones y R inc
´
ognitas seg
´
un
la ec. 11, que puede expresarse en forma resumida en la ec.
12, donde Φ es una matriz de KN × R elementos con los
resultados apilados del regresor para las K muestras, e es
el vector de error y T es el vector de torques medidos.
τ
1
.
.
.
τ
K
=
φ
1
.
.
.
φ
K
p (11)
T = Φp + e (12)
La soluci
´
on se obtiene, siempre que los R valores
singulares de Φ sean mayores que 0, por medio de la
pseudoinversa seg
´
un la ec. 13.
ˆ
p = (Φ
t
Φ)
−1
Φ
t
T (13)
El vector
ˆ
p obtenido de esta manera minimiza kek
2
.
En este punto surge investigar una serie de cuestiones
sobre la posibilidad de realizar el c
´
alculo de la ec. 13, y
qu
´
e tan buenos son los resultados en caso de poder hallarse.
Entonces se estudia la identificabilidad, la reducci
´
on a
par
´
ametros base, la incorporaci
´
on de conocimiento previo,
el efecto del ruido en la medici
´
on, y la adecuaci
´
on de las
muestras para mejorar el c
´
alculo num
´
erico.
A. Identificabilidad y reducci
´
on del modelo
En forma general no todos los par
´
ametros pueden ser
identificados debido a cuestiones referentes a la estructura
de los manipuladores. Por ejemplo, en el robot IRB140 se
observa que la masa del primer eslab
´
on (m
1
) no puede ser
hallada (aunque se puede obtener su momento de inercia
sobre el eje de giro). Sin embargo si el manipulador se
montara en un carro se podr
´
ıa excitar para observar m
1
aplicando fuerzas de desplazamiento en alguna direcci
´
on.
Vale decir que a medida que uno avanza en la cadena
cinem
´
atica m
´
as son los par
´
ametros identificables pues las
posibilidades de producir un movimiento m
´
as general son
mayores cu
´
anto m
´
as cerca se est
´
a de la herramienta. Esta
limitaci
´
on inicial se pone de manifiesto inclusive al observar
la formulaci
´
on del IDIM (ec. 9), pues de los 13 × N
par
´
ametros que podr
´
ıan estar presentes solo se observan R
de ellos.
Por otro lado, una deficiente selecci
´
on de las muestras
puede impedir la identificaci
´
on de alg
´
un par
´
ametro
identificable. El regresor φ
k
est
´
a compuesto de funciones
que no forman una base ortogonal y por lo tanto se puede
tener dependencia lineal entre ecuaciones. Sin llegar al
caso extremo, si no se cuida que la selecci
´
on de muestras
produzcan que todas las funciones base sean expuestas
se tiene un problema de p
´
erdida de generalidad en la
identificaci
´
on que se manifiesta con menor performance del
brazo al seguir otras trayectorias que las utilizadas para la
identificaci
´
on.
El n
´
umero de condici
´
on de Φ indica cuan identificable es
el sistema para un conjunto de muestras dadas. Habiendo
sido lo suficientemente cuidadosos en la elecci
´
on de ellas,
la descomposici
´
on en valores singulares (ec. 14) muestra los
par
´
ametros que no pueden identificarse (aquellos asociados
a valores singulares s
j
= 0) o los que de hacerlo se incurre
en grandes errores (s
j
< ).
Φ = UΣV
t
(14)
Sin embargo, se debe apuntar que en este contexto
los par
´
ametros est
´
an en un espacio transformado por V
t
(ec. 15). Por lo que al detectar un valor singular muy
peque
˜
no no se descartar
´
ıa un par
´
ametro del modelo sino
una combinaci
´
on de ellos.
T =
r
X
j=1
s
j
v
t
j
p
u
j
+ e
j
(15)
La reducci
´
on del modelo consiste en todas las
acciones tendientes a eliminar los par
´
ametros que no
pueden identificarse o son dif
´
ıcilmente identificables, y la
discriminaci
´
on procede en la observaci
´
on de la matriz Σ.
En [6] se propone un algoritmo donde para aquellos valores
singulares s
j
tal que s
1
/s
j
> 100 se inspecciona el vector
columna j de V , es decir v
j
, en busca de una componente
dominante, que de hallarse indicar
´
ıa un par
´
ametro candidato
a eliminar. Si se encontraran r
j
componentes dominantes,
repiti
´
endose el patr
´
on en r
j
ecuaciones entonces se est
´
a
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
85
http://elektron.fi.uba.ar
frente a una situaci
´
on en la que los correspondientes
par
´
ametros pueden hallarse solo como combinaci
´
on lineal,
pudiendo plantearse una reducci
´
on a par
´
ametros base
mediante el m
´
etodo basado en la descomposici
´
on QR de
la matriz de observaci
´
on. Es claro que el nuevo conjunto
reducido de par
´
ametros p
b
pierde el sentido f
´
ısico. Los
par
´
ametros base tambi
´
en pueden obtenerse mediante la
descomposici
´
on de Φ en valores singulares [7], [6] . En
[3] y [2] se presenta un m
´
etodo iterativo para hallar las
expresiones simb
´
olicas de agrupamiento.
Una cr
´
ıtica al algoritmo de reducci
´
on del orden surge del
hecho de trabajar con distintas escalas en los par
´
ametros,
siendo conveniente realizar un escalamiento previo (ver
III-C) para no tomar decisiones err
´
oneas. En todo caso
siempre es mejor desestimar un par
´
ametro cuando su
influencia es peque
˜
na que cuando es dif
´
ıcil de identificar.
La reducci
´
on de par
´
ametros lleva a un modelo para la
predicci
´
on, a menos que la informaci
´
on suprimida pueda
suplirse de alguna forma, como por ejemplo incorporando
el conocimiento previo conseguido por otros m
´
etodos.
Si se tienen R −r par
´
ametros que no pueden identificarse
pero pueden ser determinados por otros m
´
etodos, partiendo
de la ec. 10 se expresa la ec. 16, donde se han agrupado
los t
´
erminos conocidos con los supra
´
ındices kn y los
desconocidos con uk, se replantea el problema de estimaci
´
on
seg
´
un la ec. 18.
τ
k
= φ
uk
k
p
uk
+ φ
kn
k
p
kn
+ e (16)
τ
k
− φ
kn
k
p
kn
= φ
uk
k
p
uk
+ e (17)
τ
0
k
= φ
uk
k
p
uk
+ e (18)
B. Modelo del ruido
Durante el ensayo se relevan las posiciones q y sus
derivadas
˙
q y
¨
q. Si se considera que la variable articular
es tomada por un encoder de CPR pulsos por vuelta (∆q =
2π/CPR) y que el ruido no se debe a ning
´
un otro factor m
´
as
que a la discretizaci
´
on entonces se propone para el error
de posici
´
on un modelo Gaussiano con media nula y desv
´
ıo
σ
p
= ∆q/(2
√
3). Luego al aplicar un filtro derivador por
diferencias sucesivas se tiene para la velocidad que el desv
´
ıo
es σ
v
=
√
2σ
p
/T
s
, donde T
s
es el tiempo de muestreo
y para la aceleraci
´
on σ
a
= 2σ
p
/T
2
s
. Siendo T
s
1 es
evidente la preponderancia del t
´
ermino 1/T
s
por lo que
es indispensable acotar las componentes en frecuencia de
las se
˜
nales derivadas. Este primer grupo de incertidumbres
impactan en el c
´
alculo de Φ.
En un segundo grupo se citan otras fuentes de error
como las debidas a incertidumbres en las mediciones de los
torques aplicados (
˜
e) caracterizadas por tener media nula y
matriz de covarianza diagonal, y a din
´
amicas no modelizadas
(e
um
), como rozamientos o comportamientos el
´
asticos de
ejes y eslabones, caracterizadas por variables aleatorias de
media no nula y matriz de covarianza no diagonal. El error
e en la ec. 12 es suma de ambas.
Si se considera solo el efecto del segundo grupo (es decir
Φ determin
´
ıstico), la esperanza del error de estimaci
´
on
˜
p se
calcula [8] seg
´
un la ec. 19.
E{
˜
p} = −
Φ
t
Φ
−1
Φ
t
E{e} (19)
Por otro lado, considerando solo las del primer grupo y
ruido de medici
´
on del torque, la covarianza del error de
estimaci
´
on resulta en la ec. 20, siendo su norma seg
´
un la
ec. 21 [9], donde se ha reemplazado el desv
´
ıo est
´
andar del
ruido de medici
´
on.
cov{
˜
p} = E
Φ
t
Φ
−1
E{ee
t
} (20)
kcov{
˜
p}k =
s
σ
2
s
min
(21)
Si bien los resultados de las ecs. 19 y 20 son bien
conocidos, son dif
´
ıcilmente aplicables en el caso que se
presenta, ya que los efectos en ambos grupos se dan
simult
´
aneamente y el regresor es una funci
´
on no lineal
de las entradas por lo que se debilita la hip
´
otesis de
estacionariedad.
Sin embargo es esperable que la estimaci
´
on tenga un
sesgo (ec. 19), como as
´
ı tambi
´
en que se pueda definir un
intervalo de confianza a partir de la covarianza (ec. 20). En
particular, esto
´
ultimo es especialmente
´
util porque permite
acotar las estimaciones.
Partiendo de un conjunto de muestras e modelizada por
un proceso de ruido blanco Gaussiano, la varianza muestral
s
2
tiene una distribuci
´
on χ
2
con ν = KN − r grados de
libertad y se puede definir un intervalo de confianza α (ec.
22) para la varianza.
σ
2
min
=
(KN − 1)s
2
χ
2
α
2
< σ
2
<
(KN − 1)s
2
χ
2
1−
α
2
= σ
2
max
(22)
Por otro lado la estimaci
´
on de la media tiene una
distribuci
´
on T-Student con ν grados de libertad, por
lo que las cotas en las estimaciones de par
´
ametros
puede aproximarse para una determinada confianza α
par
combinando los resultados de la ec. 22 para σ
max
en 20.
Finalmente, un punto destacable es el resultado [10] que
muestra la forma en que se propagan los errores de medici
´
on
(ec. 23) y de modelizaci
´
on hacia las estimaciones (ec. 24),
dejando de manifiesto la importancia de realizar el ensayo
con Φ lo mejor condicionada posible (cond(Φ) es el n
´
umero
de condici
´
on de la matriz de observaci
´
on).
k
˜
pk
k
ˆ
pk
≤ cond (Φ)
k
˜
T k
kT k
(23)
k
˜
pk
kpk
≤ cond (Φ)
kδΦk
kΦk
(24)
C. Escalamiento
Los par
´
ametros representan distintas magnitudes
que cubren un rango amplio de valores, aumentando
consecuentemente el n
´
umero de condici
´
on de Φ sin
significar esto una pobre selecci
´
on de las muestras. Es decir
que la resoluci
´
on del problema num
´
erico puede mejorarse
definiendo un nuevo conjunto de valores que cubran la
misma escala, aunque el efecto de propagaci
´
on de los
errores notado en las ecs. 23 y 24 quede inalterado por este
artificio.
Adicionalmente se observa que si las fuerzas
generalizadas tienen un rango muy dis
´
ımil en cada
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
86
http://elektron.fi.uba.ar
eje tambi
´
en es conveniente realizar una transformaci
´
on
antes de la identificaci
´
on. Esto lleva a una modificaci
´
on de
los par
´
ametros identificados, hecho que no es significativo
en un modelo para la predicci
´
on.
La propuesta habitual sobre el escalamiento de par
´
ametros
consiste en hallar el factor de escala como la inversa de la
norma del correspondiente vector columna de Φ. La ec. 12
se puede reescribir en la ec. 25, donde se ha multiplicado
y dividido por la norma del vector columna. Agrupando los
t
´
erminos y definiendo h
j
= 1/kΦ
j
k se tiene la ec. 26, donde
H = diag(h).
T =
r
X
j=1
Φ
j
1
kΦ
j
k
kΦ
j
kp
j
+ e (25)
T = ΦHH
−1
p + e (26)
La nueva matriz de observaci
´
on resulta Φ
e
= ΦH y tiene
un menor n
´
umero de condici
´
on.
D. Incorporaci
´
on de estimaci
´
on previa
Otro aspecto interesante es la posibilidad de incluir una
estimaci
´
on inicial de los par
´
ametros, que puede ser obtenida
por medici
´
on directa, extra
´
ıdas de modelos CAD, o bien
de otras experiencias previas. Esto es especialmente
´
util
cuando se trabaja con algoritmos iterativos ya que aumenta
la velocidad de convergencia. En cambio en el m
´
etodo off-
line podr
´
ıa servir para balancear los resultados pesando en
mayor o menor medida la informaci
´
on previa con un m
´
etodo
denominado como m
´
ınimos cuadrados amortiguados.
Contando con un valor inicial p
0
se replantea la ec. 12
agregando el conocimiento inicial a continuaci
´
on seg
´
un la
ec. 27, donde I es la matriz identidad de dimensi
´
on dada
por el espacio de par
´
ametros.
T
p
0
=
Φ
I
p +
e
0
(27)
Definiendo
˘
p = p − p
0
,
˘
T = T − Φp
0
, e incorporando
un par
´
ametro λ para regular el peso de p
0
en la estimaci
´
on
final se tiene la ec. 28, cuya soluci
´
on por m
´
ınimos cuadrados
resulta en la ec. 29 mostrando que el t
´
ermino λ
2
I ayuda
al c
´
alculo de la inversa. Expres
´
andolo de otro modo, el
aporte de la informaci
´
on puede suplir los problemas que
se presentan cuando alg
´
un valor singular de Φ es muy bajo
[6].
˘
T
0
=
Φ
λI
˘
p +
e
0
(28)
˜
˘
p =
Φ
t
Φ + λ
2
I
−1
Φ
t
˘
T (29)
IV. PLANIFICACI
´
ON DEL ENSAYO
En III-B se explor
´
o la importancia de contar con
mediciones con la mayor relaci
´
on se
˜
nal a ruido posible, y
la dificultad que presenta la determinaci
´
on de las derivadas
de la posici
´
on de los ejes en cuanto a esto.
En un esquema de identificaci
´
on on-line la
´
unica
alternativa es aplicar un filtro pasa bajos digital. La principal
dificultad siempre presente es la decisi
´
on en cuanto a la
frecuencia de corte, ya que como la trayectoria no puede
fijarse de antemano entonces no se puede asumir nada
respecto al ancho de banda de las se
˜
nales. El criterio
adoptado consiste en tomar la frecuencia de corte cinco
veces por arriba de la frecuencia el
´
astica [11] de manera de
capturar en la mejor medida el comportamiento din
´
amico
y atenuar el ruido de medici
´
on de frecuencias superiores.
Una segunda limitante es que los filtros aplicados deben ser
causales y por lo tanto inducen una distorsi
´
on de fase [12].
Por otro lado, en la identificaci
´
on off-line se cuenta
con la realizaci
´
on completa, permitiendo que las se
˜
nales
reciban otro tratamiento. En principio el ancho de
banda de la excitaci
´
on puede limitarse a una frecuencia
conocida facilitando la derivaci
´
on y el filtrado. Adem
´
as
las trayectorias de prueba se pueden repetir innumerables
veces, y si el control del robot es lo suficientemente preciso
el promedio de las realizaciones arroja una medici
´
on de
la posici
´
on con sesgo nulo y varianza tan baja como sea
posible.
A. Trayectorias peri
´
odicas
Considerando los beneficios de contar con una excitaci
´
on
peri
´
odica [13] [14], la trayectoria deseada se puede plantear
(ec. 30) como suma de senoidales de una frecuencia base
ω
per
= 2π/T
per
, donde T
per
es el per
´
ıodo.
q(t) = q
0
+
N
h
X
h=1
A
h
sin(hω
per
t) + B
h
cos(hω
per
t)
(30)
El contenido arm
´
onico, y consecuentemente el ancho
de banda, se limita con la variable N
h
. Las derivadas
temporales de la posici
´
on tambi
´
en resultan peri
´
odicas de
ancho de banda acotado por N
h
. Por otro lado la se
˜
nal
de torque tambi
´
en es peri
´
odica y su espectro limitado en
frecuencia, pues el regresor φ incluye t
´
erminos no lineales
como multiplicaci
´
on de velocidades y funciones senoidales
de senoidales que a lo sumo duplican el contenido arm
´
onico.
Los par
´
ametros de dise
˜
no son los valores medios q
0
y
los vectores de amplitud A
h
y B
h
para cada arm
´
onica. La
pulsaci
´
on ω
per
es preferentemente com
´
un a los ejes. En total
se tienen N(2N
h
+1)+1 grados de libertad para especificar
la trayectoria.
Para implementar esta t
´
ecnica o bien se tiene acceso al
controlador [14] y se llega con las referencias calculadas
en la ec. 30 para t = kT
s
, o se genera una trayectoria
segmentada decimando el vector de muestras para usarlo
como referencias en las intrucciones de movimientos
joint cuidando que los errores de seguimiento sobre
el trayecto deseado no sean significativos. Este aspecto
es sumamente importante porque cualquier apartamiento
implica un incremento del ancho de banda de las se
˜
nales.
B. Dise
˜
no
´
optimo de la trayectoria
Una de las ventajas en cuanto a la se
˜
nal de excitaci
´
on en
el esquema de identificaci
´
on off-line es que forma parte del
dise
˜
no del ensayo, por lo que puede ajustarse la trayectoria
para excitar adecuadamente todas las caracter
´
ısticas de la
din
´
amica.
Un enfoque cuando se plantea un ensayo de identificaci
´
on
consiste en realizar movimientos simples en forma
secuencial a la vez que se calculan las din
´
amicas que
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
87
http://elektron.fi.uba.ar
ellos excitan y se incorporan en el procedimiento como
valores conocidos [15]. Si bien el m
´
etodo es sencillo y
permite dirigir el an
´
alisis a unos pocos par
´
ametros a la vez,
se produce una acumulaci
´
on de error. As
´
ı, resulta mejor
realizar un movimiento sobre una trayectoria, tomando una
cantidad de muestras sobre la misma que permita una
minimizaci
´
on global del error [16]. Combinando ambas
ideas, en [17] se propone un procedimiento basado en una
serie de trayectorias que minimizan el n
´
umero de condici
´
on
de Φ y que excitan en forma separada los componentes
de inercia, los acoplamientos inerciales, los t
´
erminos de
Coriolis, fuerza centr
´
ıpeta y los de peso propio, donde la
estimaci
´
on de los par
´
ametros se realiza al final considerando
todos los datos colectados.
En la secci
´
on III-A se discuti
´
o la importancia de
seleccionar adecuadamente las muestras del ensayo para
poder excitar las din
´
amicas a identificar. Luego en III-B
se mostr
´
o que la influencia de las perturbaciones sobre las
estimaciones guardan relaci
´
on a las muestras utilizadas para
conformar el ensayo. Por lo tanto es conveniente preguntarse
si se puede elegir el mejor conjunto de entrenamiento que
maximice la identificabilidad y minimice la varianza de la
estimaci
´
on.
Se plantea as
´
ı un problema de optimizaci
´
on, donde la
funci
´
on de costo se relaciona con los objetivos subrayados
y las variables de ajuste son los par
´
ametros que definen
la trayectoria (ec. 30), que se encuentran restringidos a un
dominio dado para que sean realizables por el robot. El
problema de buscar la mejor trayectoria se traslada entonces
a plantear la funci
´
on de costo por un lado y expresar
convenientemente las restricciones por otro.
As
´
ı se define primero un funcional, o funci
´
on de costo
O, cuyo extremo resuelve el problema. Una posibilidad es
buscar el m
´
ınimo del n
´
umero de condici
´
on de Φ [9] pues
ambos objetivos se relacionan con
´
el. De esta manera se
busca que la excentricidad s
max
/s
min
sea m
´
ınima y que la
hiperelipsoide definida se asemeje a una hiperesfera. Otra
opci
´
on equivalente es maximizar el menor de los valores
singulares, resultando en la mejora del peor caso. En [10] y
[16] se propone adem
´
as, sumar en los funcionales anteriores
t
´
erminos relacionados a los m
´
aximos y m
´
ınimos de Φ para
regular su escala. En [6] se reinterpreta el costo a trav
´
es
del alfabeto de optimalidades en relaci
´
on con los criterios
utilizados en rob
´
otica.
C. Procesamiento de se
˜
nales adquiridas
Dada la naturaleza peri
´
odica de la se
˜
nal de posici
´
on y
torque, se puede considerar cada ciclo como la realizaci
´
on
de un proceso estoc
´
astico. En cada instante de tiempo se
tiene una variable aleatoria cuyos par
´
ametros estad
´
ısticos se
estiman como,
¯q
i,k
=
1
N
c
N
c
X
j=1
q
i,k+jN
s
¯τ
i,k
=
1
N
c
N
c
X
j=1
τ
i,k+jN
s
ˆσ
2
{q
i,k
} =
1
N
c
− 1
N
c
X
j=1
(q
i,k+jN
s
− ¯q
i,k+jN
s
)
2
(31)
ˆσ
2
{τ
i,k
} =
1
N
c
− 1
N
c
X
j=1
(τ
i,k+jN
s
− ¯τ
i,k+jN
s
)
2
(32)
donde el
´
ındice i indica el eje, k es el n
´
umero de muestra
dentro del ciclo y va desde 1 hasta N
s
y N
c
es la cantidad
de repeticiones.
A medida que la cantidad de ciclos aumenta es esperable
que la varianza de la estimaci
´
on de las se
˜
nales promedio
disminuya (σ
2
¯q
i,k
, σ
2
¯τ
i,k
), por lo que se puede considerar a
las curvas de posici
´
on y torque promedios como punto de
partida para la identificaci
´
on.
Si se asume que las se
˜
nales son suma de una funci
´
on
determin
´
ıstica, dadas por las curvas promedio, y un procesos
estoc
´
astico estacionario de media nula, la varianza es
constante y su estimaci
´
on se encuentra promediando las
halladas en las ecs. 31 y 32.
ˆσ
2
q
i
=
1
N
s
N
c
− 1
N
s
X
k=1
N
c
X
j=1
(q
i,k+jN
s
− ¯q
i,k+jN
s
)
2
ˆσ
2
τ
i
=
1
N
s
N
c
− 1
N
s
X
k=1
N
c
X
j=1
(τ
i,k+jN
s
− ¯τ
i,k+jN
s
)
2
El intervalo de confianza de las varianzas σ
2
q
i
y σ
2
τ
i
se
puede construir a partir de sus estimaciones seg
´
un la ec. 22.
Asumiendo que el controlador del robot es de muy alta
performance, se suma como hip
´
otesis que la dispersi
´
on en
q se debe a ruido blanco Gaussiano aditivo, tal como se
plantea en la secci
´
on III-B, y se puede limitar su influencia al
m
´
ınimo con un filtro pasa bajos que aproveche el particular
contenido en frecuencia propuesto para la trayectoria de
excitaci
´
on. Se calcula entonces el espectro utilizando la
transformada discreta de Fourier y se elimina todo el
contenido en frecuencias superiores a la pulsaci
´
on N
h
ω
per
.
Las derivadas de la posici
´
on articular se pueden calcular
multiplicando su espectro filtrado Q(jω) por la variable
compleja jω y luego antitransformando para tener
˙
q, o
por −ω
2
y antitransformando para tener
¨
q. As
´
ı se evita
la aproximaci
´
on de la derivada por restas sucesivas y la
limitaci
´
on que puede presentarse al no poder muestrear a
altas frecuencias.
En cuanto a la se
˜
nal de torque, el contenido arm
´
onico
es mayor resultando en teor
´
ıa como m
´
aximo el doble que
el de la se
˜
nal de posici
´
on. Se sigue entonces el mismo
procedimiento duplicando la frecuencia de corte.
V. VALIDACI
´
ON
Si el objetivo de la identificaci
´
on es obtener un
modelo para la predicci
´
on, la validaci
´
on puede realizarse
comparando los resultados que arroja el modelo con las
mediciones obtenidas para una tarea tipo. Por el contrario si
el objetivo es determinar una estimaci
´
on de los valores de los
par
´
ametros din
´
amicos la validaci
´
on se realiza comparando
con los valores convencionales verdaderos si son accesibles
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
88
http://elektron.fi.uba.ar
Fig. 1: Robot IRB140 del Laboratorio de Rob
´
otica de
FIUBA
de alg
´
un modo
4
o estimando las incertidumbres en la
identificaci
´
on y procurando que sean m
´
ınimas.
En [16] se recopilan otras alternativas frecuentemente
adoptadas que se basan en comparar los resultados obtenidos
para el mismo robot variando en un caso el conjunto de
trayectorias de estimaci
´
on, o la formulaci
´
on del modelo en
otro, como por ejemplo el modelo de potencia. Por otro
lado se subraya que si el objetivo buscado es la mejora
del control, una buena forma de validar la estimaci
´
on es
calificando el comportamiento del sistema en lazo cerrado.
En el presente trabajo la validaci
´
on se realiza comparando
medici
´
on contra predicci
´
on del modelo. En este punto
es importante replicar el comportamiento de la planta
considerando las tareas para las que fue pensada, aunque
puede resultar enga
˜
noso ya que un robot en principio
no est
´
a dise
˜
nado para resolver un problema en particular
pues por definici
´
on es una m
´
aquina multiprop
´
osito. Por lo
tanto, y para simplificar la evaluaci
´
on, se consideran tareas
b
´
asicas en las que tenga sentido definir la performance
del movimiento [18], como por ejemplo desplazarse en
l
´
ınea recta entre dos puntos o describir una U invertida
tal como se realiza en una aplicaci
´
on tipo pick and place.
Es claro que cualquiera sea la tarea tipo propuesta para
validar, ninguna de ellas tendr
´
a las caracter
´
ısticas impuestas
a aquellas utilizadas durante la identificaci
´
on.
La exigencia sobre la periodicidad puede ser mantenida
repitiendo N
c
veces la tarea tipo y promediando las
realizaciones para mejorar la relaci
´
on se
˜
nal a ruido de las
mediciones de posici
´
on de los ejes y torques motor. Aunque
no es posible tener un control fino sobre el ancho de banda
por lo que se dificulta luego el filtrado.
VI. RESULTADOS
Los temas te
´
oricos y los aspectos pr
´
acticos presentados
hasta aqu
´
ı se aplican en la identificaci
´
on del IDIM de un
robot ABB-IRB140 (Fig. 1).
Para evaluar los aspectos de implementaci
´
on se limita el
estudio al brazo inferior, compuesto por los tres primeros
4
Hecho sumamente improbable pues en tal caso no tendr
´
ıa demasiado
sentido la tarea de identificaci
´
on
Fig. 2: C
´
alculo del volumen del eslab
´
on 2 del robot IRB140
mediante FreeCAD
ejes, agrupando as
´
ı los componentes articulados en la
mu
˜
neca con el eslab
´
on 3.
Al plantear las ecuaciones que definen el IDIM, se
encuentra solo la influencia de 32 de los 13 × 3 = 39
par
´
ametros iniciales. Luego se avanza en la reducci
´
on del
conjunto de par
´
ametros a estimar tomando del controlador
del robot los momentos de inercia de los rotores de los
motores (J
m
) y estimando los coeficientes de rozamiento
viscoso (B
m
) especificados en motores de similares
caracter
´
ısticas.
En repetidas experiencias se observaron mejores
resultados al sumar hip
´
otesis sobre la masas de los
eslabones. Se obtiene entonces el volumen de los eslabones
0, 1 (Fig. 2), 2, y agrupados el 3, 4, 5 y 6. Como en la
hoja de datos figura que el peso de la unidad mec
´
anica es
98 Kg, se calcula la densidad promedio y luego la masa de
los grupos de eslabones en forma proporcional al volumen.
El regresor se separa en una parte conocida φ
kn
(matriz
de 9 columnas) y una desconocida φ
un
(de 32 − 9 = 23
columnas).
Los par
´
ametros base se formulan armando la matriz de
observaci
´
on Φ a partir del regresor φ
un
para 100 valores
aleatorios de posici
´
on, velocidad y aceleraci
´
on, aplicando
luego el m
´
etodo discutido en la secci
´
on III-A. As
´
ı el
conjunto de par
´
ametros a identificar se reducen a 18 (ec.
33).
p
b
=
I
oyy1
+ I
oyy2
+ I
oyy3
+ 0.14mx
G1
F
c1
I
oxx2
− I
oyy2
I
ozz2
I
oxy2
I
oxz2
− 0.36mz
G2
− 0.36mz
G3
I
oyz2
mx
G2
my
G2
F
c2
I
oxx3
− I
oyy3
I
ozz3
I
oxy3
I
oxz3
− 0.38mz
G3
I
oyz3
mx
G3
my
G3
F
c3
(33)
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
89
http://elektron.fi.uba.ar
A. Trayectoria
´
optima
Se propone una trayectoria de per
´
ıodo T
per
= 5 s
compuesta por hasta 4 arm
´
onicos. El ancho de banda se
limita entonces a 0.8 Hz, permitiendo obtener en las se
˜
nales
de referencia del movimiento la suficiente riqueza como para
excitar las din
´
amicas, siendo a la vez aceptable en cuanto a
la frecuencia de muestreo de 20 Hz de un controlador de la
serie 5.
Los par
´
ametros de amplitud y de valor medio de las
se
˜
nales (ec. 30) se calculan de manera de conseguir que
la matriz de observaci
´
on tenga un m
´
ınimo en el n
´
umero de
condici
´
on, acotando la excursi
´
on de la posici
´
on en un 80%
de su rango y la de la velocidad y aceleraci
´
on en un 70%.
En la Fig. 3 se muestra la evoluci
´
on de cond(Φ) a la
izquierda, y las se
˜
nales que definen la trayectoria del ensayo
a la derecha para cuando el proceso de optimizaci
´
on alcanz
´
o
la cota m
´
ınima del vector de incrementos.
0 20 40 60 80 100 120 140 160 180
Iteration
0
50
100
150
200
250
300
350
400
Function value
Current Function Value: 36.3743
Stop
Pause
10 20 30 40 50 60 70 80 90 100
-2
0
2
Pos
10 20 30 40 50 60 70 80 90 100
-2
0
2
Vel
10 20 30 40 50 60 70 80 90 100
-10
0
10
Acel
Fig. 3: Izq: evoluci
´
on de cond(Φ) a medida que se ajusta
la trayectoria de prueba. Der: se
˜
nales obtenidas
La trayectoria
´
optima hallada se puede reproducir con
movimientos joint entre 20 destinos dados en coordenadas
articulares mostrados en la tabla I y distanciados 250 ms.
TABLE I: Coordenadas joint de la Trayectoria Optima
Tiempo q
1
q
2
q
3
[s] deg deg deg
0.25 24.098 84.657 −181.197
0.50 17.286 60.213 −159.677
0.75 −13.739 30.426 −137.475
1.00 −37.406 13.867 −140.872
1.25 −25.228 4.204 −166.747
1.50 7.936 −17.576 −184.000
1.75 19.932 −51.261 −168.327
2.00 −4.650 −71.724 −126.987
2.25 −34.243 −59.435 −85.765
2.50 −32.008 −26.323 −59.771
2.75 −4.401 −1.225 −46.267
3.00 8.146 6.087 −42.475
3.25 −14.428 14.048 −56.444
3.50 −45.251 39.973 −92.940
3.75 −46.080 73.231 −136.702
4.00 −15.988 88.000 −162.384
4.25 11.171 77.988 −162.589
4.50 13.406 65.175 −156.354
4.75 5.658 70.291 −163.972
5.00 11.268 85.232 −180.023
B. Identificaci
´
on sobre el simulador de ABB
El entorno de programaci
´
on off-line que complementa los
sistemas rob
´
oticos de ABB incluye un completo simulador
cinem
´
atico y din
´
amico para reproducir fielmente las tareas
de las aplicaciones rob
´
oticas antes de que ellas sean
instaladas en planta. Por tratarse de software cerrado no se
tiene acceso a los modelos ni los par
´
ametros utilizados. La
experiencia de trabajar con el software es similar en muchos
aspectos a hacerlo con el robot real.
Al igual que los robotistas ensayan sus aplicaciones en el
entorno de software, en el presente trabajo la identificaci
´
on
del IDIM se realiza en el ambiente virtual. Habiendo
probado la reproducci
´
on de las trayectorias y teniendo un
modelo identificado, se puede replicar la experiencia con
el robot real utilizando los mismos programas. En nuestro
caso surge adicionalmente una serie de dificultades [19] por
tratarse de una versi
´
on de controlador anterior:
• Per
´
ıodo de muestreo. el modelo del Laboratorio de
Rob
´
otica limita las interrupciones c
´
ıclicas a 250 ms,
aunque permite realizar interrupciones no c
´
ıclicas a 50
ms. Esto impacta en la tasa de captura de datos.
• Acceso a la se
˜
nal de torque. La versi
´
on SC4-M2000
del controlador permite leer el torque de referencia a
trav
´
es de la funci
´
on TestSignRead, aunque limita su
aplicaci
´
on a ejes externos. Se necesita entonces medir
externamente la corriente en los motores sincronizando
la adquisici
´
on con el robot.
La trayectoria del ensayo se ejecuta en el entorno virtual
(Fig. 4), y los datos adquiridos de posici
´
on y torque de
los 3 primeros ejes se salvan en un archivo. El primer y
´
ultimo ciclo se descartan porque durante parte de ellos el
robot acelera desde el reposo en un caso y en el otro frena
hasta detenerse. En la figura se observa que los trazos de
las distintas realizaciones son coincidentes, salvo el del ciclo
inicial.
Fig. 4: Ejecuci
´
on del ensayo en RobotStudio
1) Procesamiento de los datos adquiridos: El tiempo de
la interrupci
´
on c
´
ıclica destinada a tomar las muestras tiene
una variabilidad que llega a 2 ms. Se puede modelizar como
una dispersi
´
on de 0.5 ms al rededor de 51.5 ms. As
´
ı los datos
adquiridos son tomados en intervalos de muestreo variables,
por lo que el primer paso consiste en regularizarlo. Como
adem
´
as se deben identificar con el menor error los instantes
de inicio de cada ciclo para no tener realizaciones desfasadas
que influyan negativamente en la estimaci
´
on de los datos de
posici
´
on y torque, se submuestrean los datos con T
s
= 10
ms. En la Fig. 5 se muestran los datos de posici
´
on y torque
interpolados.
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
90
http://elektron.fi.uba.ar
5 10 15 20 25 30 35
Tiempo [s]
-4
-2
0
2
Posición [rad]
Datos adquiridos
5 10 15 20 25 30 35
Tiempo [s]
-2
-1
0
1
2
Torque Motor [Nm]
Fig. 5: Datos submuestreados con un per
´
ıodo T
s
= 10 ms
El inicio de cada ciclo se obtiene buscando los m
´
aximos
de la autocorrelaci
´
on de la se
˜
nal de posici
´
on en un entorno
correspondiente al tiempo de ciclo propuesto para la se
˜
nal
del ensayo. En el caso presentado bas
´
andose en la trayectoria
propuesta cada ciclo debiera tener 500 muestras espaciadas
10 ms, pero luego de analizar los datos se encuentra que son
en realidad 504 muestras por ciclo.
El espectro de las curvas de posici
´
on y torque promedios
se muestran en la Fig. 6. La posici
´
on presenta un ancho de
banda de 0.8 Hz, aunque el espectro del torque resulta m
´
as
disperso. Se proponen las frecuencias de corte f
cq
= 2.4 Hz
para la se
˜
nal de posici
´
on y f
cu
= 4.8 Hz para la del torque.
0
0.2
0.4
0.6
0.8
1
1.2
|q(f)|
Composición en frecuencia de la posición
0 1 2 3 4 5
f (Hz)
0
0.2
0.4
0.6
0.8
1
1.2
|u(f)|
Composición en frecuencia del torque
0 1 2 3 4 5
f (Hz)
Fig. 6: Espectros de las se
˜
nales promediadas de posici
´
on y
torque
En la Fig. 7 se observa a derecha la estimaci
´
on de la
posici
´
on superpuesta a los datos de los ciclos muestreados y
a la izquierda la curva de velocidad calculada en el dominio
transformado.
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-1
-0.5
0
0.5
q
1
[rad]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-3
-2
-1
0
q
2
[rad]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo [s]
0
1
2
q
3
[rad]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-2
0
2
qp
1
[rad/s]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-2
0
2
qp
2
[rad/s]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo [s]
-2
0
2
qp
3
[rad/s]
Fig. 7: Se
˜
nales de posici
´
on y velocidad estimadas
En la Fig. 8 se presentan la aceleraci
´
on y el torque
motor estimados. Se superpone en la gr
´
afica del torque las
mediciones tomadas en todos los ciclos. La estimaci
´
on del
torque se grafica junto con las bandas de m
´
aximo y m
´
ınimo
valor caracterizadas a partir de la varianza muestral de las
mediciones para incluir el 95% de las muestras.
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-20
0
20
q2p
1
[rad/s
2
]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-10
0
10
20
q2p
2
[rad/s
2
]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo [s]
-10
0
10
q2p
3
[rad/s
2
]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-1
0
1
2
u
1
[Nm]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-2
0
2
u
2
[Nm]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo [s]
-0.5
0
0.5
1
u
3
[Nm]
Fig. 8: Se
˜
nales de aceleraci
´
on y torque estimadas
2) Resultado de la identificaci
´
on: Con los datos
preprocesados se corre la identificaci
´
on. La estimaci
´
on lineal
sirve de condici
´
on inicial de un proceso de identificaci
´
on
no lineal e iterativo basado en la estimaci
´
on previa de la
covarianza del torque. As
´
ı se refinan los par
´
ametros dando
mayor peso a las mediciones m
´
as confiables. En la tabla II
se muestran los valores estimados de los par
´
ametros base
junto con su intervalo de confianza del 95%.
TABLE II: Par
´
ametros Base Estimados
Param Ident Min Max
I
oyy1
+ I
oyy2
+ I
oyy3
+ 0.14mx
G1
6.2672 6.2651 6.2693
F
c1
0.45214 0.45203 0.45225
I
oxx2
− I
oyy2
-2.1315 -2.1343 -2.1287
I
ozz2
1.9237 1.9218 1.9256
I
oxy2
0.10048 0.098863 0.10209
I
oxz2
− 0.36mz
G2
− 0.36mz
G3
0.34712 0.3459 0.34834
I
oyz2
0.18018 0.17905 0.18131
mx
G2
-8.1397 -8.1408 -8.1386
my
G2
-0.40546 -0.40657 -0.40434
F
c2
0.40672 0.40663 0.40681
I
oxx3
− I
oyy3
-3.3542 -3.3552 -3.3532
I
ozz3
4.0412 4.0409 4.0414
I
oxy3
-0.62208 -0.62266 -0.62151
I
oxz3
− 0.38mz
G3
-0.29632 -0.29656 -0.29609
I
oyz3
0.17383 0.17367 0.17399
mx
G3
-9.6518 -9.6519 -9.6516
my
G3
-2.2081 -2.2083 -2.2078
F
c3
0.28581 0.2858 0.28582
Para evaluar la bondad del ajuste, se corre la validaci
´
on
sobre los datos de identificaci
´
on (Fig. 9).
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-1
0
1
2
Torque [Nm]
Validación del modelo identificado
Medido
Prediccion
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-2
0
2
Torque [Nm]
Medido
Prediccion
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo [s]
-0.5
0
0.5
Torque [Nm]
Medido
Prediccion
Fig. 9: Comparaci
´
on del torque estimado para la tarea de
identificaci
´
on y su predicci
´
on por el modelo identificado
para los 3 ejes
Se realiza una nueva identificaci
´
on para evaluar la
incidencia de la frecuencia de corte de los filtros aplicados al
preprocesamiento, limitando el espectro de los datos hasta
las frecuencias f
cq
= 1.2 Hz y f
cu
= 2.4 Hz. Con esta
nueva identificaci
´
on se busca que el modelo refleje mejor
las tendencias. En la Fig. 10 se observa que el ajuste resulta
m
´
as grosero.
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
91
http://elektron.fi.uba.ar
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-1
0
1
Torque [Nm]
Validación del modelo identificado
Medido
Prediccion
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
-2
0
2
Torque [Nm]
Medido
Prediccion
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Tiempo [s]
-0.5
0
0.5
Torque [Nm]
Medido
Prediccion
Fig. 10: Comparaci
´
on del torque estimado para la tarea de
identificaci
´
on y su predicci
´
on considerando f
cq
= 1.2 Hz
C. Validaci
´
on
El modelo identificado se valida evaluando su capacidad
de predicci
´
on del torque necesario para realizar la tarea
t
´
ıpica que se muestra en la Fig. 11.
Fig. 11: Ejecuci
´
on de la trayectoria de validaci
´
on en
RobotStudio
Las curvas capturadas en este caso est
´
an compuestas
por menor cantidad de muestras. El procedimiento para
procesarlas es el mismo al detallado para la trayectoria
de identificaci
´
on. En la Fig. 12 se muestran las se
˜
nales
estimadas y su intervalo de confianza del 95% comparadas
con los datos preprocesados.
0 0.2 0.4 0.6 0.8 1 1.2
-0.5
0
0.5
q
1
[rad]
0 0.2 0.4 0.6 0.8 1 1.2
-1.4
-1.3
-1.2
-1.1
q
2
[rad]
0 0.2 0.4 0.6 0.8 1 1.2
Tiempo [s]
3.2
3.25
3.3
q
3
[rad]
0 0.2 0.4 0.6 0.8 1 1.2
-2
0
2
u
1
[Nm]
0 0.2 0.4 0.6 0.8 1 1.2
-2
-1
0
u
2
[Nm]
0 0.2 0.4 0.6 0.8 1 1.2
Tiempo [s]
-1
-0.5
0
0.5
u
3
[Nm]
Fig. 12: Se
˜
nales de posici
´
on y torque estimadas comparadas
con los datos adquiridos
En la Fig. 13 se muestra la curva de torque estimada
comparada con la predicci
´
on. En todos los ejes la predicci
´
on
resulta dentro del intervalo de confianza de la estimaci
´
on del
torque. Se observa que la predicci
´
on aparece como una curva
promedio de espectro acotado en relaci
´
on a la estimada a
partir de los datos.
0 0.2 0.4 0.6 0.8 1 1.2
-1
0
1
Torque [Nm]
Validación del modelo identificado
Medido
Prediccion
0 0.2 0.4 0.6 0.8 1 1.2
-2
-1
0
Torque [Nm]
Medido
Prediccion
0 0.2 0.4 0.6 0.8 1 1.2
Tiempo [s]
-0.5
0
0.5
Torque [Nm]
Medido
Prediccion
Fig. 13: Comparaci
´
on del torque estimado y su predicci
´
on
Luego se eval
´
ua la capacidad de predicci
´
on del modelo
identificado usando filtros con frecuencias de corte f
cq
=
1.2 Hz y f
cu
= 2.4 Hz. En la Fig 14 se observa que si
bien este modelo no ajusta adecuadamente las muestras de
identificaci
´
on, produce mejores resultados en cuanto a la
predicci
´
on.
0 0.2 0.4 0.6 0.8 1 1.2
-1
0
1
Torque [Nm]
Validación del modelo identificado
Medido
Prediccion
0 0.2 0.4 0.6 0.8 1 1.2
-2
-1
0
1
Torque [Nm]
Medido
Prediccion
0 0.2 0.4 0.6 0.8 1 1.2
Tiempo [s]
-0.5
0
0.5
Torque [Nm]
Medido
Prediccion
Fig. 14: Comparaci
´
on del torque estimado y su predicci
´
on
por el modelo identificado para los 3 ejes, considerando
f
cq
= 1.2 Hz en la identificaci
´
on
VII. CONCLUSIONES
Se present
´
o el marco te
´
orico que permite derivar el
modelo param
´
etrico de la din
´
amica de un manipulador
serie en forma sistem
´
atica. Atendiendo a la imposibilidad
de realizar mediciones directas de ciertos par
´
ametros ya
sea porque implique desensamblar un robot o bien porque
ellos se manifiestan solo durante el movimiento, se plante
´
o
el problema de la identificaci
´
on. Un primer enfoque se
bas
´
o en el ajuste por m
´
ınimos cuadrados aprovechando la
posibilidad de reformular el modelo din
´
amico inverso en una
versi
´
on lineal en los par
´
ametros. Para refinar la estimaci
´
on
se incluyeron luego algoritmos no-lineales iterativos. Como
el objetivo del ajuste est
´
a puesto en minimizar el error
cuadr
´
atico medio entre la estimaci
´
on y el conjunto de datos
medidos se produce un modelo con una buena capacidad de
predicci
´
on, sustentando de esta forma el objetivo del trabajo
que consiste en generar un modelo
´
util para el estudio del
control, dejando en un segundo plano la idea de tener una
buena medici
´
on de ciertos par
´
ametros din
´
amicos.
Se detect
´
o la influencia en el modelo de un conjunto de
par
´
ametros y se concluy
´
o que correspond
´
ıan al conjunto
m
´
ınimo identificable o par
´
ametros base.
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
92
http://elektron.fi.uba.ar
Se estudi
´
o el modelo del ruido y la incidencia en las
estimaciones de los par
´
ametros. Se present
´
o un marco
estad
´
ıstico que justifica la validez de las estimaciones y
permite evaluar la calidad de las mismas.
Para mejorar el tratamiento de las se
˜
nales, atendiendo a
que el muestreo est
´
a limitado a una muy baja frecuencia,
y que solo se cuenta con se
˜
nales de posici
´
on teniendo
que estimar las velocidades y aceleraciones de los ejes,
se incorpor
´
o un esquema de trayectorias de identificaci
´
on
peri
´
odicas limitadas en ancho de banda. As
´
ı se redujo el
ruido de medici
´
on de las se
˜
nales de torque y posici
´
on a
la vez que se implement
´
o el c
´
alculo de las velocidades y
aceleraciones en el espacio transformado.
Se incorpor
´
o el dise
˜
no del ensayo al procedimiento de
identificaci
´
on, generando la trayectoria para maximizar la
excitaci
´
on de las din
´
amicas asociada a los par
´
ametros por
un lado, y adecu
´
andolas a las posibilidades de reproducci
´
on
del robot. Se comprob
´
o que la capacidad de predicci
´
on del
modelo identificado mejora en virtud de la selecci
´
on de la
trayectoria de identificaci
´
on.
Finalmente se trabaj
´
o sobre el modelo IRB140 colectando
los datos mediante el simulador del RobotStudio, e
identificando luego su IDIM que puede ser aplicado en la
investigaci
´
on de t
´
ecnicas de control entre otros.
Como trabajos futuros se destacan, la ejecuci
´
on de las
trayectorias de prueba en el robot IRB140 del laboratorio.
En este caso aparece un desaf
´
ıo adicional debido a una
limitaci
´
on en relaci
´
on a la adquisici
´
on de datos.
REFERENCES
[1] M. W. Spong and M. Vidyasagar, Robot Dynamics and Control. NY:
Wiley & Sons, 1989.
[2] D. Zeida, “Desarrollo de un paquete de simulaci
´
on para robots en
un ambiente con manejo de ecuaciones simb
´
olicas (mathem
´
atica).”
Tesina de Graduaci
´
on en Ingenier
´
ıa Electr
´
onica, FIUBA, 1996.
[3] A. Oubi
˜
na, “Identificaci
´
on de par
´
ametros din
´
amicos de robots para
calibraci
´
on y control adaptativo.” Tesis de Graduaci
´
on en Ingenier
´
ıa
Electr
´
onica, FIUBA, 1996.
[4] M. Anigstein, “Din
´
amica de manipuladores rob
´
oticos.” Tesis
Doctoral, FIUBA, 2002.
[5] A. Luca and W. Book, Robots with Flexible Elements, pp. 287–319.
01 2008.
[6] J. M. Hollerbach, W. Khalil, and M. Gautier, “Model identification,”
in Springer Handbook of Robotics, pp. 113–138, 2016.
[7] C. Atkeson, C. An, and J. Hollerbach, “Estimation of inertial
parameters of manipulator loads and links,” Int J Robot Res, vol. 5,
no. 3, pp. 101–119, 1986.
[8] M. Gautier, “Identification of robots dynamics,” IFAC Proceedings
Volumes, vol. 19, no. 14, pp. 125–130, 1986.
[9] B. Armstrong, “On finding exciting trajectories for identification
experiments involving systems with nonlinear dynamics,” The
International Journal of Robotics Research, vol. 8, no. 6, pp. 28–
48, 1989.
[10] M. Gautier and W. Khalil, “Exciting trajectories for the identification
of base inertial parameters of robots,” The International Journal of
Robotics Research, vol. 11, no. 4, pp. 362–375, 1992.
[11] M. Brunot, A. Janot, P. Young, and F. Carrillo, “An improved
instrumental variable method for industrial robot model
identification,” Control Engineering Practice, vol. 74, 03 2018.
[12] W. Khalil and E. Dombre, Modeling, Identification and Control of
Robots. USA: Taylor & Francis, Inc., 2002.
[13] J. Swevers, C. Ganseman, D. B. Tukel, J. de Schutter, and
H. Van Brussel, “Optimal robot excitation and identification,” IEEE
Transactions on Robotics and Automation, vol. 13, pp. 730–740, Oct
1997.
[14] J. Swevers, W. Verdonck, and J. Schutter, “Dynamic model
identification for industrial robots,” Control Systems, IEEE, vol. 27,
pp. 58 – 71, 11 2007.
[15] K. Yoshida, N. Ikeda, and H. Mayeda, “Experimental study of
the identification methods for an industrial robot manipulator,” in
Proceedings of the IEEE/RSJ International Conference on Intelligent
Robots and Systems, vol. 1, pp. 263–270, July 1992.
[16] J. Wu, J. Wang, and Z. You, “An overview of dynamic
parameter identification of robots,” Robotics and Computer-Integrated
Manufacturing, vol. 26, no. 5, pp. 414–419, 2010.
[17] P. O. Vandanjon, M. Gautier, and P. Desbats, “Identification of robots
inertial parameters by means of spectrum analysis,” in Proceedings
of 1995 IEEE International Conference on Robotics and Automation,
vol. 3, pp. 3033–3038 vol.3, May 1995.
[18] ISO, “Manipulating industrial robots – performance criteria and
related test methods,” standard, International Organization for
Standardization, Geneva, CH, Apr. 1998.
[19] ABB, System Data Types and Routines. ABB Robotics AB, S-721
68 Vasteras, Sweden, 4.0 ed., 2000.
Revista elektron, Vol. 5, No. 2, pp. 83-93 (2021)
ISSN 2525-0159
93
http://elektron.fi.uba.ar
