375
IDENTIFICAÇÃO EM MALHA FECHADA DE PARÂMETROS DE MODELOS DE ATRITO
EM UM PÊNDULO INVERTIDO
Hugo Tanzarella Teixeira1 , Victor Semedo de Mattos Siqueira2 , Celso José Munaro3
1
2
Universidade Federal do Espírito Santo, Vitória, Brasil, [email protected]
Universidade Federal do Espírito Santo, Vitória, Brasil, [email protected]
3
Universidade Federal do Espírito Santo, Vitória, Brasil, [email protected]
Resumo: Neste artigo, os atritos presentes em um carrinho
que atua sobre um pêndulo invertido são estimados através de
duas abordagens, usando dados obtidos em malha fechada.
O primeiro método utiliza os ciclos limite originados pelo
atrito para obter os parâmetros do modelo de Coulomb mais
atrito estático, através do ajuste de elipses aos dados medidos. O segundo método estima os parâmetros do modelo de
atrito de Karnopp via mínimos quadrados, usando para isto
dados medidos e parâmetros do pêndulo. Os dois métodos
são comparados através de sua aplicação para um modelo de
simulação e para dados medidos no pêndulo real.
Palavras-chave: Aplicações de sistemas dinâmicos, pêndulo
invertido, estimação de atrito.
1. INTRODUÇÃO
Todo sistema mecânico está sujeito ao atrito e sua presença produz ciclos limite que podem impedir sua estabilização caso compensações adequadas não sejam utilizadas [1].
Por se tratar de um sistema mecânico subatuado, inerentemente instável em malha aberta e de dinâmica não linear, o
pêndulo invertido é um problema clássico e já bastante estudado para aplicações de técnicas de controle e por isso foi
escolhido neste trabalho como plataforma de testes para o
presente estudo.
Modelos de atrito têm sido utilizados para gerar estimativas do atrito a ser compensado. Em [2] o modelo de
Karnopp foi utilizado para compensar o atrito em um servomecanismo. A compensação de atrito também pode ser realizada através da adição de um sinal pulsado ao sinal de controle usando o modelo de Coulomb. Em [3] esta abordagem
foi utilizada para compensar o atrito em uma válvula de controle. Em ambos os casos é necessário obter os parâmetros
do modelo de atrito para serem usados pelo compensador.
A proposta deste trabalho é apresentar métodos para identificação em malha fechada dos parâmetros dos modelos de
atrito. Assim é possível quantificar e compensar o atrito em
processos de maneira automática, evitando que o processo
seja parado para que métodos invasivos sejam realizados para
identificação dos parâmetros dos modelos de atrito. Como
o atrito surge com o uso, esta abordagem também permite
detectar o momento em que o atrito se torna significativo e
precisa ser compensado.
2. AMBIENTE DE TESTES
O sistema experimental utilizado, produzido pela
QuanserTM , consiste em um carrinho com uma haste móvel
acoplada através de uma articulação. O carrinho move-se
horizontalmente em um trilho e a haste móvel por sua vez
pode realizar movimento rotacional de 360◦ no plano vertical. O carrinho consiste de uma base de alumínio impulsionada por um motor CC de 400 W, trifásico sem escova. O
sistema físico está interligado a um ambiente de controle em
tempo real (QuaRC) integrado ao MatlabTM /SimulinkTM
que permite seu controle e monitoramento. Uma biblioteca
de funções permite a interface entre as variáveis do ambiente
Simulink e o sistema real.
2.1. Modelo do pêndulo invertido
Pelo método de Newton-Euler obtemos as equações que
descrevem a dinâmica do sistema do pêndulo invertido
(M + m)ẍ = Fe − mL sen(θ)θ2 + mL cos(θ)θ̈ − Fat
(J + mL2 )θ̈ = mLg sen(θ) + mL cos(θ)ẍ − fat
(1)
onde M é a massa do carrinho, m a massa da haste, L a
distância do pivô da haste até seu centro de gravidade, J o
momento de inércia da haste, Fe a força aplicada ao carrinho,
Fat a força de atrito no carrinho, fat o atrito resistente ao
momento de rotação da haste, x é a posição do carrinho e θ a
posição angular da haste. Em nosso sistema a força aplicada
ao carrinho é dada por
Fe = βi
(2)
376
onde β é a constante de conversão entre a corrente aplicada
(i) e a força que atua no motor CC utilizado. Neste trabalho o
atrito (fat ) na haste será desconsiderado. Os valores, fornecidos pelo fabricante, dos parâmetros destas equações estão listados na Tabela 1(a).
3. MODELOS DE ATRITO
3.1. Modelo de Coulomb mais atrito estático
Uma maneira simples de modelar o atrito é através do
modelo de Coulomb mais atrito estático

 Fc sgn(ẋ) se ẋ 6= 0
Fe
se ẋ = 0 |Fe | < Fs
Fat =
(3)

Fs sgn(Fe ) caso contrário
onde Fc é o coeficiente de atrito de Coulomb e Fs o coeficiente de atrito estático. Embora seja um modelo muito simples, captura o comportamento essencial do atrito e é útil em
muitas aplicações.
3.2. Modelo de Karnopp
O modelo de Karnopp [4] é adequado para a detecção da
velocidade nula em simuladores. A criação de uma faixa de
valores, dentro da qual a velocidade do movimento é considerada nula, define a velocidade zero no intervalo |ẋ| < DV .
Devido à assimetria na força de atrito o modelo de Karnopp
será descrito através das equações
Fat

−
Fatrito
(ẋ)



−

 Fs sgn(Fe )
Fe
=


 Fs+ sgn(Fe )


+
Fatrito
(ẋ)
4. MÉTODOS DE IDENTIFICAÇÃO DE PARÂMETROS DOS MODELOS DE ATRITO
4.1. Parâmetros do modelo de atrito de Coulomb
A presença de atrito em processos sob certas condições
produz ciclos limite [1] que tem sido explorados na literatura para realizar uma estimativa dos atritos em válvulas
de controle [5]. Estes métodos são agora aplicados ao pêndulo, tendo em vista o comportamento oscilatório em malha
fechada (Figura 1). De modo particular, o gráfico i versus x
obtido, é típico de sistemas sob influência de atrito estático
e que apresentam oscilações e evidencia o comportamento
da posição do carrinho de acordo com a corrente aplicada.
As regiões em aproximadamente -50 mm e 15 mm representam o carrinho parado e estão diretamente relacionadas com
o atrito do sistema.
Algoritmo proposto para estimar o atrito, baseado em [5].
1. Aplicar ao conjunto de dados (i, x) a um filtro passa
baixa a fim de eliminar componentes de alta frequência
e gerar os sinais filtrados (if , xf )
2. Escolher um segmento de dados (if , xf ) com base na
regularidade de oscilações desses sinais. Tal regularidade é verificada a partir do critério apresentado em [6].
3. Ajustar uma elipse ao gráfico if versus xf .
4. Traçar um segmento de reta paralelo ao eixo i passando
pelo centro da elipse.
ẋ ≤ −DV
5. Calcular a distância entre os interceptos do segmento de
|ẋ| < DV e |Fe | ≤ Fs−
−
+
reta traçado com a reta i = 0 e com a elipse.
|ẋ| < DV e Fs < Fe < Fs
+
|ẋ| < DV e |Fe | ≥ Fs
Este algoritmo pode ser executado periodicamente, de
ẋ ≥ DV
forma
automática, usando dados coletados que indicaram a
(4)
presença
de oscilações [6]
−
Fatrito
(ẋ) = Fc− sgn(ẋ) + Fv− ẋ
(5)
se
se
se
se
se
+
Fatrito
(ẋ) = Fc+ sgn(ẋ) + Fv+ ẋ
(6)
onde Fv é o coeficiente de atrito viscoso e os coeficientes
positivo e negativo denotam a direção do movimento.
Para realização dos testes em simulação, serão usados o
modelo do pêndulo invertido mais os modelos de atrito, com
parâmetros listados na Tabela 1(b).
4.1.1
Identificação dos parâmetros
Para a simulação utilizou-se o modelo de atrito de Coulomb
mais atrito estático com coeficientes de atrito listados na
Tabela 1(b), o resultado obtido pode ser observado na Figura
2(a). É importante evidenciar que o uso do algoritmo que
Tabela 1 – Parâmetros
(b) Modelos de atrito
(a) Pêndulo
Parâmetro
Valor
Parâmetro
Valor
M
3,220 kg
Fs+
6,0 N
m
0,230 kg
Fs−
10,0 N
L
0,641 m
Fc+
5,0 N
J
0,008 kg.m2
Fc−
9,0 N
g
2
Fv+
30,0 N.s/m
32,396 N/A
Fv−
25,0 N.s/m
DV
0,002 m/s
β
9,810 m/s
Figura 1 – Comportamento oscilatório no pêndulo invertido
377
descrito a seguir é aplicado apenas quando o carrinho está
em movimento. Portanto, o coeficiente de atrito estático Fs
não é estimado, ao final desta seção discutiremos alternativas
para sua estimação.
O primeiro passo é definir o vetor de parâmetros a serem
estimados
λ = Fc Fv
(7)
Quando a velocidade do carrinho tem módulo maior que
DV , a força de atrito é dada por Fatrito , substituindo θ̈ em ẍ
no sistema (1), obtemos a seguinte equação
m2 L2
˙
cos(θ) − (M + m) ẍ − mL sen(θ)(θ)
Fe +
J + mL2
m2 L2
+
g sen(θ) cos(θ) = Fc sgn(ẋ) + Fv ẋ (8)
J + mL2
(a) Simulação
Note que, sendo os estados medidos no processo e os parâmetros do pêndulo invertido previamente conhecidos, (8)
é linear em relação ao vetor de parâmetros λ. Em seguida
define-se o vetor de regressores como
ϕ = sgn(ẋ(t)) ẋ(t)
(9)
Chamando
Figura 2 – Elipse ajustada ao gráfico i - x
ajusta a elipse por mínimos quadrados aos dados do gráfico
i - x claramente estima o atrito de Coulumb. Para estimar o
atrito estático, a elipse deve envolver todos os dados. A estratégia aqui proposta e expandir os raios da elipse através de
uma busca unidimensional até que todos os pontos pertençam
ao interior da mesma. No caso mostrado na Figura 2(a) os
raios da elipse externa foram aumentados 33% em relação à
elipse interna. Na Tabela 2(a) temos a comparação entre os
coeficientes de atrito estimados e os utilizados para simulação. Na Tabela 2(b) temos a média e o desvio padrão (σ) dos
coeficientes de atrito de Coulomb estimados com dados do
pêndulo real. O desvio padrão foi foi calculado após a realização de sete testes consecutivos, e seu valor pequeno indica
que a estimativa é consistente. Neste caso os rais da elipse
foram aumentados 46% para estimar Fs+ e Fc− .
Tabela 2 – Coeficiente de atrito de Coulomb estimado
(b) Testes na planta
(a) Simulado
Fs−
F̂s−
Fc−
F̂c−
Média
σ
Fs−
14,12
0,68
F̂c+
Fs+
10,28
0,46
5,00 5,01
Fc−
10,86
0,30
Fc+
7,14
0,62
10,00 11,44 9,00 9,02
Fs−
6,00
F̂s−
7,47
Fc+
m2 L2
cos(θ(t)) − (M + m) ẍ(t)
∆(t) = Fe +
J + mL2
m2 L2
˙
− mL sen(θ(t))(θ)(t)
+
g sen(θ(t)) cos(θ(t))
J + mL2
(10)
(b) Planta
o vetor de parâmetros pode ser obtido minimizando o erro
quadrático
X
2
λ̂ = argλ min
∆(t) − ϕ(t)λT
(11)
t
Entretanto, os períodos em que (8) é válida são desconhecidos, pois a velocidade limite DV é uma incógnita do
modelo. Para lidar com este problema, define-se a variável
δv(s) [4], tal que
s
(12)
δv(s) = max (|ẋ(t)|)
Z
sendo Z >> 1 e S < Z, para s = 1, 2, ..., S, ou seja,
δv(s) assume S valores distintos, obrigatoriamente menores
do que a maior velocidade verificada nos dados experimentais. Para cada valor de δv(s) o vetor de regressores ϕ(t)
e ∆(t) são escolhidos dos dados experimentais segundo a
condição |ẋ| > δv(s) e em seguida o vetor de parâmetros é
estimado resolvendo (11).
Como Fs é necessário para o modelo, as alternativas para
estimá-lo são: a) usar estimativas obtidas offline, b) adotar
Fs = αFc , com α > 1, c) usar o procedimento discutido na
seção 4.1.1.
4.2.1
4.2. Parâmetros do modelo de Karnopp
O método descrito a seguir foi inicialmente proposto para
aplicação em válvulas de controle em [4], sendo aqui adaptado para o pêndulo invertido. Entretanto o método que será
Identificação dos parâmetros
A fim de validar o método descrito na seção anterior, o
mesmo foi aplicado a dados simulados e dados reais do pêndulo invertido. Para isso utilizou-se o modelo de atrito de
Karnopp com coeficientes de atrito listados na Tabela 1(b).
378
O comportamento da estimação do vetor de parâmetros para diferentes valores de δv(s), usando o modelo é
mostrado na Figura 3(a). Nota-se que para δv(s) << DV o
valor dos parâmetros estimados varia significantemente para
diferentes valores de δv(s), pois os dados correspondem a
instantes em que (8) não se aplica. Por outro lado, à medida
que s aumenta e δv(s) ≈ DV , as estimativas aproximam-se
dos valores reais dos parâmetros e se mantém quase constantes, mesmo quando δv(s) ultrapassa o valor ideal de DV .
Para os dados reais, devido aos ruídos de medição, os
estados devem ser filtrados antes de serem submetidos ao
método de identificação. Para isso utilizou-se um filtro FIR,
de ordem 1200 e frequência de corte 1 Hz e tempo de
amostragem 0,5 ms, por apresentar um atraso linear para
todos os sinais. Na Figura 3(b) podemos ver o comportamento da estimação do vetor de parâmetro para testes realizados na planta real. Nota-se que para δv(s) > 0, 02 m/s
não há variação significativa no valor dos parâmetros estimados. Na Tabela 3 temos o valor médio dos parâmetros estimados e seu desvio padrão (σ) para os sete testes realizados,
onde observa-se a similaridade (para Fc ) com os resultados
da Tabela 2(b).
(a) Simulação
Tabela 3 – Parâmetros do modelo de Karnopp estimados
(b) Testes na planta
(a) Simulado
Média
σ
9,0
Fc−
9,49
0,30
25,0
Fv−
6,841
1,54
5,0
Fc+
5,91
0,33
30,0
Fv+
10,40
2,68
Modelo Estimado
Fc−
9,0
Fv−
25,0
Fc+
Fv+
5,0
30,0
(b) Planta
Figura 3 – Comportamento do vetor de parâmetros estimado
5. CONCLUSÃO
Dois métodos para estimação de parâmetros de modelos
de atrito foram apresentados e analisados via simulação e
com dados reais. O método que estima os parâmetros do
modelo de Coulomb e atrito estático faz uso das oscilações
geradas pelo atrito para através delas estimar elipses a gráficos de corrente versus posição. O método se destaca pela
simplicidade e por necessitar apenas de dados de operação.
A obtenção dos parâmetros do modelo de Karnopp requer os
parâmetros do modelo do pêndulo, além de não fornecer uma
estimativa para o atrito estático. Sua utilização torna-se vantajosa caso o atrito não produza ciclos limite, além de conter
também um parâmetro relacionados à velocidade.
AGRADECIMENTOS
Os autores agradecem ao CNPq e a FAPES pelo apoio aos
programas de pós-graduação e iniciação cientifica.
REFERÊNCIAS
friction compensation based on karnopp model,” Proceedings of the European Control conference, pp. 2558
– 2563, 2001.
[3]DOI L. Z. X. Ivan and S. Lakshminarayanan, “A new unified
approach to valve stiction quantification and compensation,” Industrial & Engineering Chemistry Research,
vol. 48, pp. 3474 – 3483, 2009.
[4]PUB R. A. Romano and C. Garcia, “Karnopp friction model
identification for a real control valve,” Proceedings of
the 17th World Congress the International Federation of
Automatic Control, 2008.
[5]DOI M. A. A. S. Choudhury, S. L. Shah, N. F. Thornhill, and
D. S. Shook, “Automatic detection and quantification of
stiction in control valves,” Control Engineering Practice,
vol. 14, pp. 1395 – 1412, 2006.
[6]DOI N. F. Thornhill, B. Huang, and H. Zhang, “Detection of
multiple oscillations in control loops,” Journal of Pro[1]DOI H. Olsson and K. J. Aström,“Friction generated limit
cess Control, vol. 13, pp. 91 – 100, 2003.
cycles,” IEEE Trans. Contr. Syst. Technol, vol. 9, no. 4,
pp. 629 – 636, 2001.
[2]PUB L. Ravanbod-Shirazi and A.Besançon-Voda,“Robust