f
8º CONGRESSO IBEROAMERICANO DE ENGENHARIA
MECANICA
Cusco, 23 a 25 de Outubro de 2007
CONTROLE ROBUSTO APLICADO AO MODELO MATEMÁTICO DE
SISTEMA MÚSCULO ESQUELÉTICO INFERIOR HUMANO
Pergher R.1 , Bottega V.1 e Nabinger E.2
1: Departamento de Matemática e Estatística - DEME
Centro de Ciências Exatas e Tecnologia - CCET
Universidade de Caxias do Sul - UCS
Rua Francisco Getúlio Vargas, 1130 – Caxias do Sul – RS - Brasil
e-mail: [email protected] , [email protected]
2:Instituto Brasileiro de Tecnologia do Couro, Calçado e Artefatos
Rua Araxá, 750 - Bairro Ideal - Novo Hamburgo – RS - Brasil
e-mail: [email protected]
RESUMO
O modelo biomecânico de um sistema músculo-esquelético humano e a simulação do comportamento do
sistema em movimento pode ser aplicado nas mais diversas áreas, com destaque na medicina, no tratamento
de doenças motoras e no esporte, para um melhor desempenho de atletas.
Este trabalho tem como objetivo obter um modelo dinâmico representando o sistema músculo-esquelético de
uma perna. A descrição da cinemática e da dinâmica dos movimentos dos membros é baseada na formulação
de Newton-Euler e Euler-Lagrange. Os movimentos e forças resultantes são produzidos por conjuntos de
atuadores músculo-tendinídeos. Estes modelos dinâmicos são não-lineares, com múltiplas entradas e múltiplas
saídas, muitos graus de liberdade e com fortes acoplamentos. Isto resulta em grandes dificuldades na
determinação de parâmetros como forças e momentos. Por isso, o modelo deve ser genérico o suficiente para
aceitar diversos modelos de músculo e técnicas de controle. Neste trabalho, foi utilizado um controle baseado
na teoria VSS (variable structure system) que tem como característica principal, robustez com relação a
variações nos parâmetros.
Um modelo geométrico para simulações é obtido e implementado com o software Matlab/Simulink. Este
modelo permite simulações de diversos movimentos como pedalar ou andar.
PALAVRAS CHAVE: modelagem, biomecânica, controle, sistema músculo-esquelético.
INTRODUÇÃO
Modelos matemáticos de sistemas biológicos constituem um dos campos de grande crescimento na área
científica. Apesar de todas as técnicas de modelagem matemática e simulação disponíveis, sua aplicação
generalizada em sistemas músculo-esquelético exige um aporte de esforços e recursos proporcional à sua
complexidade. Simulações computacionais podem ser vistas como uma extensão do experimento de análise
de movimento em dois aspectos: primeiro, modelagem e simulação podem dar informações que não são
diretamente acessíveis por experimentos em humanos e, segundo, os dados de simulação do modelo podem
ser muito úteis para explicar os resultados obtidos do experimento de análise de movimento. Dados obtidos de
ambos experimentos, in vitro e in vivo, são totalmente desenvolvidos e usados em modelos computacionais.
O modelo músculo-esquelético é uma representação matemática idealizada do corpo humano, compreendendo
ossos, músculos, juntas e estruturas passivas de vários graus de complexidade. O modelo de membros
inferiores, inclui uma base formada pela pelve, duas pernas (cada uma contendo fêmur, tíbia, patela, tálus,
calcâneo e dedos). Este modelo compreende nove graus de liberdade incluindo as rotações e translações do
quadril, joelho e tornozelo [5]. O sistema músculo-esquelético transforma sinais elétricos do sistema nervoso
periférico em grandezas mecânicas, como força e deslocamento. Existem dois modelos tradicionais para o
estudo de movimentos humanos, dinâmica direta e dinâmica inversa, ambas podem ser usadas para
determinar variáveis cinemáticas como momentos e forças nas articulações durante o movimento.
MODELO DINÂMICO
Posição e Orientação de um Corpo Rígido
Uma perna pode ser esquematizada, do ponto de vista mecânico, como uma cadeia cinemática aberta,
formada por corpos rígidos conexos em cascata por meio de juntas rotacionais (articulações). Um extremo da
cadeia é vinculado a uma base (pelve) e, o outro, ao elemento terminal (pé). O movimento da estrutura (perna)
é realizado mediante a composição dos movimentos elementares de cada membro, com respeito ao
precedente. Para isto, faz-se necessária uma derivação das equações cinemáticas da perna, descrevendo a
posição e a orientação do membro terminal (pé), em função das variáveis de junta (articulações), em relação a
um sistema de coordenadas de referência. Estas equações podem ser obtidas através da convenção de
Denavit-Hartenberg [1].
Figura 1: Representação de um ponto P em sistemas coordenados diferentes.
A matriz ortogonal de rotação do sistema O1, com respeito ao sistema O0, é dada por
⎡ x1T x0
⎢
R10 = ⎢ x1T y 0
⎢ x1T z 0
⎣
(1)
y1T x0
y1T y 0
y1T z 0
z1T x0 ⎤
⎥
z1T y 0 ⎥
z1T z 0 ⎥⎦
A posição do ponto P, com relação ao sistema de referência O0, é expressa como
p 0 = o10 + R10 p 1
(2)
onde p1 é o vetor das coordenadas de P.
Jacobiano Geométrico
p& e a velocidade angular w do elemento
&
terminal, em função da velocidade das variáveis de junta q , mediante relações do tipo
Quer-se expressar, como vetores livres, a velocidade linear
p& = J p (q )q& ,
w = J 0 (q )q&
(3)
ou, na forma compacta,
⎡J ⎤
v = ⎢ p ⎥ q& = J (q )q&
⎣J0 ⎦
(4)
onde a matriz de transformação J6xn é denominada jacobiano geométrico.
Particionando o jacobiano geométrico na equação (4) em vetores coluna (3x1), tem-se
⎡J p
J =⎢ 1
⎣ J O1
L J pn ⎤
,
L J On ⎥⎦
(5)
onde o termo q& i J p representa a contribuição da junta i à velocidade linear do membro terminal, enquanto
i
que o termo q& i J O i representa a contribuição desta junta à velocidade angular do membro terminal.
Formulação de Lagrange
A dedução do modelo dinâmico do sistema é de grande importância para a simulação do movimento, para
a análise de estruturas de manipulação e para a elaboração de algoritmos de controle. Este modelo descreve as
relações existentes entre o torque aplicado às juntas e o movimento da estrutura. Na formulação de Lagrange,
que é conceitualmente simples e sistemática, a derivação das equações do movimento é independente do
sistema de coordenadas de referência. Um conjunto de variáveis λ i , i=1,...,n , denominadas coordenadas
generalizadas, descrevem a posição dos elementos mecânicos que constituem uma estrutura com n graus de
liberdade. Neste caso, λ i = q i , onde qi representa o ângulo das juntas do sistema. Define-se como
lagrangeano do sistema mecânico a função, dependente das coordenadas generalizadas,
L = Τ −U
(6)
onde T é a energia cinética e U é a energia potencial total do sistema.
A partir da função de Lagrange, obtém-se a equação de Lagrange, dada por
d ∂L ∂L
−
= ξ fi
dt ∂q& i ∂qi
,
i=1,...,n
(7)
onde
ξf
i
são as forças generalizadas associadas às coordenadas generalizadas qi , as quais provêm das forças
não conservativas, devido ao torque gerado pelos atuadores e torque induzido por situações de interação com
o ambiente.
A equação (7) define as relações existentes entre as forças generalizadas aplicadas no sistema e a velocidade e
a aceleração de suas juntas.
Energia Cinética
A contribuição da energia cinética do elemento i é expressa em função das coordenadas generalizadas do
sistema, utilizando o jacobiano geométrico,
T=
1 n n
1
bij (q)q& i q& j = q& T B(q)q& ,
∑∑
2 i =1 j =1
2
(8)
onde
n
(
B (q) = ∑ mli J lpiT J lpi + J 0( li )T I li J 0( li )
)
i =1
(9)
é denominada a matriz de inércia (nxn) simétrica, positiva definida e dependente da configuração do sistema.
Energia Potencial
A energia potencial total é dada por
(
U = −∑ mli g 0T pli
)
(10)
onde g0 é o vetor aceleração da gravidade, com relação ao sistema de coordenadas base.
Equação do Movimento
Considerando as expressões (8) e (10), a derivação do lagrangeano do sistema, equação (6), resultam as
equações do movimento do sistema [2]
n
∑b
j =1
ij
n
n
(q )q&& j + ∑∑ hijk (q )q& k q& j + g i (q ) = ξ f i
i=1,...,n
j =1 k =1
(11)
onde
hijk =
∂bij
∂q k
−
∂b jk
∂qi
(12)
Uma interpretação física da equação do movimento evidencia que o coeficiente bii representa o momento de
inércia na junta i quando as outras juntas estão bloqueadas, enquanto que o coeficiente bij leva em conta o
efeito da aceleração da junta j sobre a junta i. Os termos
hijk (q)q& k q& j representam o efeito centrífugo e de
coriolis na junta i, provocado pela velocidade da junta j e k. Os termos gi representam o torque gerado no eixo
da junta i devido ao efeito da gravidade.
Se o sistema está em contato com um ambiente, parte dos torques de atuação devem compensar os torques
induzidos nas juntas pelas forças de contato. Tais torques são dados por JT(q)ha, onde ha denota o vetor de
forças exercidas pelo elemento terminal do sistema sobre o ambiente.
Enfim, as equações do movimento (11) podem ser reescritas na forma matricial compacta, representando o
modelo dinâmico no espaço das juntas
B (q)q&& + C (q, q& )q& + g (q ) = τ − J T (q )ha
(13)
onde Cnxn satisfaz a relação
n
n
n
j =1
j =1 k =1
∑ cij q& j = ∑∑ hijk (q)q& k q& j .
(14)
MODELO DINÂMICO DO MEMBRO INFERIOR HUMANO
Para efeito de simplicidade, considera-se o sistema em duas dimensões formado por três articulações no
quadril, joelho e tornozelo, livre de forças externas de contato.
Figura 2: Sistema músculo-esquelético
Considera-se o sistema da figura 21, para o qual o vetor das variáveis generalizadas resulta em
q = [q1 , q 2 , q3 ]T que
representam os ângulos das articulações do quadril, joelho e tornozelo,
respectivamente. Sejam m1 , m 2 e m3 as massas dos elementos, l1 , l 2 e l 3 as medidas dos baricentros da
coxa, perna e pé, respectivamente, e a1, a2 e a3 os comprimentos da coxa, perna e pé, respectivamente.
Escolhidos os sistemas de coordenadas, o cálculo dos jacobianos de velocidade linear (5) dos baricentros
fornece
J
( l1 )
p
⎡− l1 s1
= ⎢⎢ l1c1
⎢⎣ 0
J
( l3 )
p
0 0⎤
0 0⎥⎥ ,
0 0⎥⎦
J
( l2 )
p
⎡− l1 s1 − l 2 s12
= ⎢⎢ l1c1 + l 2 c12
⎢⎣
0
⎡− l 3 s123 − l 2 s12 − l1 s1
= ⎢⎢ l 3 c123 + l 2 c12 + l1c1
⎢⎣
0
− l 2 s12
l 2 c12
− l 3 s123 − l 2 s12
l3 c123 + l 2 c12
0
0
0⎤
0⎥⎥ ,
0⎥⎦
− l 3 s123 ⎤
l 3 c123 ⎥⎥
0 ⎥⎦
(15)
onde c(12...n) e s(12...n) indicam, respectivamente, cosseno e seno de ( θ 1 + θ 2 + ... + θ n ).
Os jacobianos da velocidade angular dos baricentros fornecem
J
( l1 )
0
⎡0 0 0 ⎤
= ⎢⎢0 0 0⎥⎥ ,
⎢⎣1 0 0⎥⎦
J
( l2 )
0
⎡0 0 0 ⎤
= ⎢⎢0 0 0⎥⎥ ,
⎢⎣1 1 0⎥⎦
J
( l3 )
0
⎡0 0 0 ⎤
= ⎢⎢0 0 0⎥⎥
⎢⎣1 1 1⎥⎦
(16)
CONTROLE
O sistema músculo-esquelético em humanos é controlado por uma complexa rede neuronal biológica,
responsáveis pela formulação de algumas estratégias e pela geração de sinais de controle motor. Controlar a
distribuição de forças entre diversos músculos que cruzam uma articulação e que produzem uma determinada
ação é, do ponto de vista matemático, um problema de otimização. Neste sentido, o controle neuronal, pode
ser substituído por um modelo de controle, que visa determinar as excitações neuro-musculares que levam o
sistema de um estado para outro, partindo apenas das equações de movimento e de vínculos ao mesmo tempo
que é minimizado um funcional, geralmente integral, que depende dos estados dinâmicos e das variáveis de
controle ao longo de todo o movimento [3].
Controladores robustos têm sido intensamente estudados por possuirem numerosas vantagens comparadas
com controladores ótimos convencionais. Sua principal característica é a robustez contra erros de modelagem,
como distúrbios externos e variações nos parâmetros do sistema físico. Estas técnicas de controle fornecem
ações naturais contra estes distúrbios. Utiliza-se uma técnica de controle robusto (sliding mode control),
baseada na teoria VSS (variable structure system).
Linearização
A idéia básica de linearização do modelo matemático do sistema é construir uma lei de controle não linear,
que, no caso ideal, linearize exatamente o modelo matemático em questão. Assim, podemos utilizar um
segundo estágio no modelo para o controle do sistema [4].
Uma das técnicas de linearização mais conhecida é a de controle de dinâmica inversa que lineariza o sistema,
resultando num sistema linear composto de subsistemas, onde cada subsistema é afetado por apenas uma
entrada de controle, tornando o sistema linear e desacoplado.
Considera-se o problema de seguir uma trajetória específica qd no espaço dos ângulos das articulações. Por
simplificação, escreve-se o modelo dinâmico, equação (13), na forma
B (q)q&& + n(q, q& ) = u
(17)
onde
n(q, q& ) = C (q, q& )q& + Fv q& + g (q ) .
(18)
Assume-se u dependente do modelo músculo-esquelético e
u = B(q ) y + n(q, q& ).
(19)
Usando a propriedade da matriz B ser inversível , o sistema, na equação (17), com a lei de controle, equação
(19), denominada controle de dinâmica inversa, resulta num sistema linear desacoplado descrito por
q&& = y .
(20)
A equação (19) reduz o problema de controle à determinação de uma lei y que torne o sistema estável.
Controle Robusto
A idéia dos controladores robusto VSS é obter características robustas, alterando sua lei de controle, baseada
em algum critério, restringindo a trajetória para uma superfície específica no espaço estado. Como resultado,
a resposta do sistema fica insensível a variações de parâmetros e distúrbios externos. A superfície é definida
de tal maneira que a sua trajetória possua algumas características, tal como, estabilidade assintótica. No
entanto, requer a estimativa de um limite superior nas incertezas para especificar o ganho da retroalimentação
u = Bˆ (q) y + nˆ (q, q& ) u
.
y = q&&d + K D q~& + K P q~ + Ω
.
(21)
Escolhe-se
(22)
⎛ ρr ⎞
⎟
⎜ z ⎟ z r representa a contribuição de robustez contrastando os distúrbios não modelados e
⎝ r ⎠
O termo Ω = ⎜
incertezas nos parâmetros da equação do movimento do sistema músculo-esquelético.
SIMULAÇÕES E RESULTADOS
Modelo e Dados
Para testar as leis de controle obtidas anteriormente, considera-se um modelo simplificado tendo o torque
como variável de entrada de controle
B(θ )q&& + C (θ , q& )q& + g = u
.
(23)
Trajetória Desejada
Para verificar o desempenho dos controladores apresentados, utilizou-se uma trajetória circular para a
extremidade do pé com centro em x=0.4 m e y=0.2m e raio de r=0.2m simulando um movimento de pedalada
livre, obtida através da expressão parametrizada para 0 ≤ t ≤ 2π
x = r cos t + 0.4
y = r sin t + 0.2
(24)
A partir da trajetória da extremidade do pé dada acima obtém-se a trajetória do ângulo das articulações
resolvendo um problema de cinemática inversa mostrada na figura 3(a).
tragetoria percorrida
2.5
angulo quadril
angulo joelho
2
1.5
[rad]
1
0.5
0
-0.5
-1
-1.5
0
1
2
3
4
5
6
7
[t]
(a)
( b)
Figura 3: (a) Trajetória desejada para ângulo das articulações do sistema músculo-esquelético.
(b) Trajetória percorrida pelo ângulo das articulações do sistema músculo-esquelético.
Simulações
O movimento do sistema músculo esquelético foi simulado em PC, utilizando MatLab/Simulink com uma
implementação em tempo discreto, com período de tempo ∆t = 1 ms, com o método numérico para solução
de equações diferenciais Runge Kutta de quarta ordem, por um período de 6 segundos. Na figura 3(b), podese verificar a trajetória percorrida pelo modelo simulado. Na figura 4, tem-se o erro entre a trajetória desejada
para o ângulo das articulações e a trajetória traçada pelo modelo matemático. Pode-se verificar que o erro de
trajetória permanece próximo de zero com uma pequena oscilação provocada pelas forças gravitacionais
agindo sobre o transcorrer da trajetória.
Figura 4: Erro de trajetória do sistema músculo-esquelético.
Esta simulação mostra o bom desempenho do sistema de controle robusto apresentado tanto no estado
transiente como no estado estacionário.
CONCLUSÕES
Neste trabalho, obteve-se o modelo cinemático, através da convenção de Denavit-Hatemberg, que descreve a
posição e a orientação do membro terminal (pé) em função das variáveis de junta (articulações), em relação a
um sistema de coordenadas de referência, situado na base (pelve). As relações entre a velocidade das
articulações e as velocidades linear e angular dos membros, foram descritas através do jacobiano geométrico e
a derivação das equações do movimento do sistema, numa forma fechada, baseado na formulação de
Lagrange. Posteriormente, obteve-se uma lei de controle robusto baseada na teoria VSS (variable structure
system). Esta técnica de controle fornece ações naturais contra variações nos parâmetros do sistema. Para
validar o modelo matemático do sistema músculo-esquelético e a lei de controle robusto, simulações foram
com o software Matlab/Simulink. Optou-se por um modelo simplificado onde a variável de controle é o
próprio torque das articulações. Uma trajetória circular foi usada para simular um movimento de pedalada
livre, desconsiderando forças de contato. As simulações se mostraram satisfatórias com o sistema
percorrendo a trajetória desejada com um erro próximo de zero. No entanto, apesar de apresentar bons
resultados, um modelo que tenha o torque como variável de controle esta distante do sistema fisiológico. Em
virtude disto, a este modelo descrito, deve-se, em trabalho futuro, considerar a ação de forças externas e
acrescentar um modelo de sistema atuador músculo-tendíneo, que relacione ativação, velocidade e
comprimento do atuador com a força na extremidade do tendão. Também se faz necessário a comparação dos
resultados com dados experimentais com os quais, a resposta do modelo pode ser verificada, como, por
exemplo, dados de vídeo, freqüentemente usados para verificar cálculos de movimento de juntas e
eletromiografia de músculos, importante para validar cálculos de força do músculo.
REFERÊNCIAS
1.
V. Bottega, Controle e Otimização Estrutural de Manipuladores Robóticos com Elementos Flexíveis
usando Atuadores Piezelétricos, Tese de Doutorado, UFRGS/PROMEC, Porto Alegre, 2005.
2.
S. C. Cannon and G.I. Zahalak, The Mechanical Behavior of Active Human Skeletal Muscles in Small
Oscillations, Journal of Biomechanics, vol.15, pp. 111-121, 1982.
3.
A. Kuo, An Optimal Control Model for Analyzing Human Postural Balance, IEEE Transactions on
Biomedical Engineering, Vol. 42, No. 1, (1995).
4.
M. W. Spong, M. Vidysagar, Robot Dynamics and Control, John Wiley and Sons, USA, (1989).
5.
D.G. Thelen and F.C. Anderson, Using computed muscle control to generate forward dynamic
simulations of human walking from experimental data, Journal of Biomechanics, vol. 39, pp. 1107-1115,
2006.
NOMENCLATURA
B(q) matriz de inércia (adimensional)
g0 vetor aceleração da gravidade (m/s2)
J
jacobiano geométrico (adimensional)
L
lagrangeano (adimensional)
p1 vetor das coordenadas do ponto (m)
p& velocidade linear (m/s)
q&
velocidade das variáveis de junta (adimensional)
coordenadas generalizadas (adimensional)
qi
R10 matriz de rotação ortogonal (adimensional)
T
U
w
ξf
i
energia cinética (adimensional)
energia potencial (adimensional)
velocidade angular (rd/s)
forças generalizadas (adimensional)