Embed Size (px)
Citation preview
CALIBRAÇÃO DE MODELO CONSTITUTIVO DE INTERFACE
POR MEIO DE FILTRO DE KALMAN ESTENDIDO
Gabriel Lucas Sousa da Silva
Projeto de Graduação apresentado ao
Curso de Engenharia Mecânica da Escola
Politécnica, Universidade Federal do Rio de
Janeiro, como parte dos requisitos
necessários à obtenção do título de
Engenheiro.
Orientador: Daniel Alves Castello
Rio de Janeiro
Março/2016
CALIBRAÇÃO DE MODELO CONSTITUTIVO DE INTERFACE POR MEIO DE FILTRO DE KALMAN ESTENDIDO
Gabriel Lucas Sousa da Silva
PROJETO DE GRADUAÇÃO SUBMETIDO AO CORPO DOCENTE DO CURSO DE
ENGENHARIA DE MECÂNICA DA ESCOLA POLITÉCNICA DA UNIVERSIDADE
FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS NECESSÁRIOS
PARA A OBTENÇÃO DO GRAU DE ENGENHEIRO MECÂNICO.
Examinada por:
________________________________________
Prof.º Daniel Alves Castello, DSc.
________________________________________
Prof.ª Lavinia Maria Sanabio Alves Borges, DSc.
________________________________________
Prof.º Hélcio Rangel Barreto Orlande, PhD.
RIO DE JANEIRO – RJ, BRASIL
Março/2016
iii
S. da Silva, Gabriel Lucas
Calibração De Modelo Constitutivo De Interface Por
Meio De Filtro De Kalman Estendido/ Gabriel Lucas Sousa
da Silva - Rio de Janeiro: UFRJ/ Escola Politécnica, 2016.
IX, 44 p.: il.; 29, 7cm.
Orientador: Daniel Alves Castello
Projeto de Graduação – UFRJ/ Escola Politécnica/
Curso de Engenharia Mecânica, 2016.
Referências Bibliográficas: p. 43-44.
1.Adesão 2.Modelo de Interface 3.Problemas
Inversos 4.Filtro de Kalman 5.Calibração I. Alves Castello,
Daniel. II. Universidade Federal do Rio de Janeiro, Escola
Politécnica, Curso de Engenharia Mecânica. III.
Calibração De Modelo Constitutivo De Interface Por Meio
De Filtro De Kalman Estendido
iv
AGRADECIMENTOS
De início, eu gostaria de agradecer meus amigos Beatriz, Filipe, Camilla, e
Nathana, os quais em alguns (muitos) momentos me apoiaram (e suportaram) ao longo
de cinco anos de faculdade, em especial nesses últimos meses. Aos amigos da
mecânica, Jorge, Daniel e Felipe. E é claro a Fabiana, pessoa que conheci muito em
tão pouco tempo, que conviveu com dificuldades similares para concluir este projeto.
Gostaria de agradecer aos meus professores orientadores, Daniel e Lavínia, pela
oportunidade de realizar este trabalho. E nessa mesma linha, ao novo amigo William, e
aos técnicos do Laboratório de Acústica e Vibrações, que ajudaram na confecção de
um aparato experimental que infelizmente não pode ser concluído a tempo.
Reconhecimento claro deve-se também a todos os meus professores do
Departamento de Engenharia Mecânica da UFRJ, sem os quais, obviamente, isso não
seria possível.
E por fim, mas não menos importante, à minha família. Em especial minha mãe,
meu pai e minha irmã, que esperaram ansiosos pelo dia em que eu finalmente receberia
o diploma.
v
Resumo Do Projeto De Graduação Apresentado À Escola Politécnica/ UFRJ Como
Parte Dos Requisitos Necessários Para A Obtenção Do Grau De Engenheiro Mecânico.
CALIBRAÇÃO DE MODELO CONSTITUTIVO DE INTERFACE POR MEIO DE FILTRO DE KALMAN ESTENDIDO
Gabriel Lucas Sousa da Silva
Março/2016
Orientador: Daniel Alves Castello
Curso: Engenharia Mecânica
Este trabalho apresenta um modelo constitutivo de interface simplificado, o qual é capaz de fornecer estimativas quanto a integridade da interface adesiva.
O objetivo deste trabalho é apresentar uma estratégia que pode ser utilizada para a calibração deste modelo; aqui utilizando-se um algoritmo recursivo de um Filtro de Kalman Estendido.
É proposto um sistema físico, o qual pode ser utilizado como aparato experimental viável para a identificação dos parâmetros constitutivos do adesivo. Para este sistema é possível obter relações matemáticas explícitas entre os parâmetros do modelo e as variáveis de estado do sistema.
É realizada uma análise global de sensibilidade de modo a se avaliar a influência da incerteza nos parâmetros constitutivos em uma dessas variáveis de estado, no caso, a força de reação. Observou-se que para o sistema proposto, as variáveis às quais o sistema era mais sensível são rigidez e energia de decoesão, sendo a viscosidade apenas influente em casos onde seu valor é alto.
Considerando-se o problema de calibração do modelo, uma série de dados sintéticos é gerada e utilizada como base para testes da viabilidade do procedimento. São avaliadas a influência da condição inicial escolhida para o processo de estimação, e a robustez dos resultados em uma condição específica de ruído e condição inicial.
Em regime de viscosidade baixa, obteve-se resultados satisfatórios para os parâmetros estimados, dada a condição inicial considerada.
Palavras-chave: Adesão, Modelo de Interface, Problemas Inversos, Filtro de Kalman, Calibração
vi
Abstract of Undergraduate Project presented to POLI/UFRJ as a partial fulfillment of
the requirements for the degree of Mechanical Engineer
CALIBRATION OF CONSTITUTIVE ADHESION MODEL BASED ON EXTENDED KALMAN FILTERING
Gabriel Lucas Sousa da Silva
March/2016
Advisor: Daniel Alves Castello
Course: Mechanical Engineering
In this work, a simplified adhesion model is presented. The model considered is capable of generating estimations about the integrity of the adhesive.
This works objective is to present a strategy to be used for the calibration of this adhesion model, by means of an Extended Kalman Filter algorithm.
A physical system is proposed as a viable experimental setup, which is able to produce explicit relations between the model parameters and the system’s state variables. The effect of the uncertainties in the constitutive parameters is analyzed for one of these state variables.
This analysis showed that the stiffness and the decohesion energy parameters are the ones most influent, while the viscosity parameter only seemed to be influential when its value was high enough.
Considering the calibration of the model, a series of numerical experiments synthetic data is presented. In these experiments, the influence of the initial conditions on the recursive EKF algorithm is analyzed for the parameter estimation problem, and a short study on the variability of the code due to different realizations of the noisy data is presented, given predefined initial conditions.
Given these initial conditions, the results are unsatisfactory only for a high viscosity regime.
keywords: Adhesion, Interface Model, Inverse Problems, Kalman Filter, Calibration
Lista de Figuras
1 Região de adesão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Comportamentos Previstos - 1 . . . . . . . . . . . . . . . . . . . . . . 83 Comportamentos Previstos - 2 . . . . . . . . . . . . . . . . . . . . . . 94 Viga engastada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Histórico de carregamento . . . . . . . . . . . . . . . . . . . . . . . . 136 Deslocamento segundo análise global . . . . . . . . . . . . . . . . . . 147 Influência da rigidez c . . . . . . . . . . . . . . . . . . . . . . . . . . 158 Influência da energia w . . . . . . . . . . . . . . . . . . . . . . . . . . 169 Influência da viscosidade b . . . . . . . . . . . . . . . . . . . . . . . . 1710 Influência de todos os parâmetros . . . . . . . . . . . . . . . . . . . . 1811 PDF da reação máxima rmáx . . . . . . . . . . . . . . . . . . . . . . . 1912 PDF do Tempo de Falha T0.2 . . . . . . . . . . . . . . . . . . . . . . 2013 Influência das Condições Iniciais - 1 . . . . . . . . . . . . . . . . . . . 3614 Influência das Condições Iniciais - 2 . . . . . . . . . . . . . . . . . . . 3715 Influência das Condições Iniciais - 3 . . . . . . . . . . . . . . . . . . . 3816 PDF dos parâmetros em ruído . . . . . . . . . . . . . . . . . . . . . . 3917 Incerteza propagada . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
vii
Lista de Tabelas
1 Valores de referência para análise de velocidade . . . . . . . . . . . . 72 Variáveis adimensionais de referência. . . . . . . . . . . . . . . . . . . 133 Comportamento da reação máxima rmáx . . . . . . . . . . . . . . . . 194 Comportamento do tempo de falha T0.2 . . . . . . . . . . . . . . . . . 21
viii
Sumário
1 Introdução 11.1 Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Revisão Bibliográfica . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.4 Organização do Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Modelo Constitutivo 52.1 Formulação Matemática . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Comportamentos previstos . . . . . . . . . . . . . . . . . . . . . . . . 7
3 Análises Numéricas 103.1 Sistema Físico Proposto . . . . . . . . . . . . . . . . . . . . . . . . . 103.2 Formulação Adimensional . . . . . . . . . . . . . . . . . . . . . . . . 123.3 Análise Global de Sensibilidade . . . . . . . . . . . . . . . . . . . . . 13
4 Filtro de Kalman 234.1 O Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 234.2 Filtro de Kalman Discreto . . . . . . . . . . . . . . . . . . . . . . . . 254.3 Filtro de Kalman Estendido Discreto . . . . . . . . . . . . . . . . . . 274.4 EKF Contínuo-Discreto . . . . . . . . . . . . . . . . . . . . . . . . . . 284.5 Estimação Conjunta de Parâmetros . . . . . . . . . . . . . . . . . . . 31
5 Estimação de Parâmetros 325.1 Formulação em Espaço de Estados . . . . . . . . . . . . . . . . . . . 325.2 Influência das Condições Iniciais . . . . . . . . . . . . . . . . . . . . . 345.3 Influência do Ruído de Medição . . . . . . . . . . . . . . . . . . . . . 39
6 Conclusão 42
Referências Bibliográficas 43
ix
1 Introdução
1.1 Motivação
Existem situações em engenharia complexas demais para se analisar de outraforma que não através de simulação. Análises essas que se baseiam em modelosmecânicos de comportamento do sistema analisado. Num contexto de adesivosindustriais, esses modelos envolvem relações constitutivas que devem ser capazes deprever o comportamento do adesivo em diversas condições de aplicação.
Sauer [1] muito recentemente publicou um survey de métodos computacionaispara contato adesivo no qual distingue os modelos em três abordagens:
1. Modelos locais de material, os quais descrevem o comportamento do materialdo adesivo. Tais modelos buscam relacionar o tensor de tensão σ (x, t) com otensor de deformação ε (x, t);
2. Modelos locais de interface, os quais descrevem o comportamento adesivo nainterface entre os substratos. Estes modelos relacionam a força de tração t(x, t)com o deslocamento relativo entre os substratos;
3. Modelos globais, os quais descrevem o comportamento efetivo previsto pelosoutros dois modelos em seus domínios de atuação. Estes podem ser utilizadospara construir modelos hierárquicos multi-escala.
Entretanto, antes da utilização de quaisquer desde modelos em situações reais,há a necessidade de calibrá-los, de modo a que o comportamento previsto por elesseja condizente com a realidade. Considerando-se o modelo constitutivo de interfaceproposto por Raous et al[2], que acopla contato unilateral, adesão, e a integridade dainterface adesiva, existe a necessidade de se calibrar parâmetros correspondentes àviscosidade do adesivo, rigidez e energia máxima absorvida antes do início de dano.
Do ponto de vista de controle de uma junta adesiva, informações sobre aintegridade do adesivo são atrativas. A partir disso, foi escolhido um algoritmo defiltragem recursivo que pode facilmente ser utilizado para aplicações de monitoraçãoem tempo real, um Filtro de Kalman [3].
O filtro de Kalman é um algoritmo recursivo em duas etapas que utiliza umasérie de observações ao longo do tempo, contendo ruído e outras incertezas, junto aum modelo dinâmico de um sistema para produzir estimativas de seus estados quetendem a ser mais “precisas” que aquelas obtidas puramente pela medição. Nessesentido, o filtro se presta tanto a estimativas sequenciais (ou online), quando as
1
medições são tomadas e processadas conforme chegam, ou em batelada (ou offline),quando um conjunto de medições é recolhido e processado por inteiro posteriormente.
Aplicações comuns são em navegação e controle de veículos terrestres e aéreos,tracking de objetos ou (controle de) variáveis de processo, fusão de informação vindade diferentes fontes, e processamento de imagem, som e sinais em geral, dentre outros.
No contexto de problemas inversos não-estacionários, as medições queconstituem o conjunto de dados para o problema inverso são tomadas num ambientenão-estacionário. Os parâmetros físicos que queremos estimar podem serdependentes do tempo, e as medições dependem dos valores destes parâmetros emdiferentes instantes de tempo. Filtros Bayesianos, como o Filtro de Kalman, seprestam a esta tarefa por possibilitar a incorporação de informação recolhida emdiferentes instantes [4].
1.2 Objetivos
O objetivo deste texto é tratar de uma versão simplificada do modelo propostopor Raous et al[2]. São abordados os comportamentos previstos por este modelosimplificado e é proposto um método de calibração do modelo. A viabilidade dométodo será analisada em diferentes de regimes de sensibilidade do sistema comrelação a um dos parâmetros a serem determinados.
Com relação à estimação de parâmetros, esta deve ser realizada através deanálise inversa [4], conjugando medições de um sistema físico com conhecimento apriori dos parâmetros. Para tanto é utilizada como ferramenta deste trabalho umFiltro de Kalman Estendido.
1.3 Revisão Bibliográfica
Diversos autores no contexto de calibração de modelos de interface, dentre elesestão Wang et al. [5], que apresentou um processo de identificação para um modelode interface viscoelástico através da utilização de um algoritmo genético.
Valoroso e Fedele [6] apresentaram o caso de calibração dos parâmetrosconstitutivos para um modelo coesivo de fratura de Modo 1 em juntas adesivas,onde a análise inversa foi realizada utilizando-se um algoritmo Gradiente, baseadoem região de confiança, e fazendo-se a otimização de uma função custo.
Albuquerque et al. [7] utilizaram uma abordagem Bayesiana, utilizando umalgoritmo MCMC modificado, para a solução do problema de identificação dosparâmetros constitutivos para um modelo de Raous simplificado, tal como é o destetrabalho.
2
Dentre os autores que utilizaram Filtro de Kalman como abordagem paraestimação de parâmetros podemos citar Corigliano e Mariani [8], que discutiram oproblema de identificação de parâmetros de modelos de interface em compósitosutilizando Filtro de Kalman Estendido. Posteriormente, investigaram duasformulações distintas do Filtro de Kalman Estendido e sua influência naidentificação de parâmetros [9].
Outros trabalhos que utilizaram Filtro de Kalman são os de Bolzon et al.[10], que apresenta um procedimento utilizando-se iterações globais do filtro para aidentificação dos parâmetros de um modelo coesivo bilinear para materiais similaresao concreto; e de Xu et al. [11] para um modelo coesivo exponencial para propagaçãode trinca, para o qual os parâmetros foram estimados utilizando-se tanto Filtro deKalman, algoritmo genético e um procedimento misto.
Quanto a filtros de Kalman especificamente, diversos trabalhos podem sercitados. Entretanto, valem citar inicialmente os trabalhos de Odelson et al [12] eSolonen et al [13]. Em princípio, para se utilizar o filtro corretamente é necessário umcerto “tuning” de alguns parâmetros, como as matrizes de covariância do processoobservado. Odelson et al [12] apresentam um método de estimar essas matrizesatravés de dados operacionais comuns, assim como garantir que as matrizes sejampositivas semidefinidas. Enquanto que Solonen et al [13] apresentam um método quepermite, através da parametrização das componentes, complexidade variada para aestimação dessas matrizes.
Almagbile et al [14] apresentam o efeito de se utilizar métodos adaptativospara as matrizes de covariância presentes no filtro. Os métodos foram implementadosem um sistema GPS/INS, e testados utilizando-se dados reais.
Miller e Howard [15] apresentam um Filtro de Kalman escrito na forma“contínua-discreta” (como é de acordo com a nomenclatura adotada neste trabalho)para a identificação dos coeficientes de um modelo equivalente de segunda ordempara um mancal, incluindo amortecimento e rigidez.
Ding et al [16], também usando uma formulação mista “contínua-discreta”,utilizam o filtro de Kalman para estimar simultaneamente a excitação, decompostaem uma base de funções ortogonais, quanto os parâmetros de modelos linear enão-linear de uma estrutura.
Em sua tese, Foun [17] apresenta um estudo utilizando-se filtros de Kalmancom características adaptativas para monitorar a integridade estrutural baseado emdados de vibração em diversas intensidades diferentes de ruído de medição. Foramconsiderados modelos tanto lineares quanto não-lineares.
3
1.4 Organização do Trabalho
Este trabalho é organizado como segue. Na Seção 2, é apresentada umaformulação matemática simplificada do modelo constitutivo de interface propostopor Raous et al. [2], seguida de um breve estudo quanto aos comportamentosreproduzíveis pelo modelo num cenário de tração uniaxial.
Na Seção 3, um sistema físico é proposto como meio de se estudar a influênciados parâmetros constitutivos do modelo de interface em um cenário mais elaborado.Uma análise de sensibilidade global é realizada com o intuito de se quantificar ainfluência da incerteza nos parâmetros constitutivos em variáveis de interesse comoforça de reação da interface adesiva e tempo de falha.
Na Seção 4, é apresentado o algoritmo padrão do Filtro de Kalman, egradualmente as modificações feitas ao mesmo até se chegar no algoritmo utilizadoneste trabalho.
Na Seção 5, o problema de calibração do modelo constitutivo é abordadoutilizando o Filtro de Kalman apresentado na Seção 4.4, utilizando-se dadosexperimentais sintéticos. A formulação matemática em Espaço de Estadosconsiderando a evolução das variáveis constitutivas do adesivo é apresentada, e sãoanalisadas a influência de um chute inicial para o processo de calibração sob diversascondições, e sob uma dessas condições a estabilidade dos resultados quanto a ruídonas observações.
Na Seção 6, as conclusões são apresentadas.
4
2 Modelo Constitutivo
Um modelo de adesão deve ser capaz de descrever:
1. A relação entre o deslocamento relativo das superfícies em contato e odesenvolvimento das forças de reação
2. A transição entre adesão completa e separação total
3. O efeito da perda de adesão na capacidade de carga
Esta seção apresenta uma versão simplificada do modelo mecânico de adesãoproposto por Raous et al.[2], o qual é um modelo bastante geral que leva emconsideração contato unilateral na fronteira, o modelo de fricção de Coulomb (“atritoseco”) e/ou adesão. Considera-se no entanto que esta versão simplificada é suficientepara descrever a interface analisada, onde fricção e esforço tangencial não sãoaplicados, e portanto desconsiderados.
2.1 Formulação Matemática
Para apresentar o modelo simplificado, primeiro considere a fig.1.
Figura 1: Região de adesão
Seja A um corpo em contato com uma superfície rígida B através da interfaceS, com espessura desprezível. Sob a hipótese de pequenas deformações podemosdefinir um vetor unitário n, normal ao corpo A em cada ponto de contato com S, demesma direção e sentido contrário à normal da superfície B.
5
O movimento vertical de A a partir da superfície rígida é descrita pelodeslocamento relativo un = u · n, onde u é o deslocamento de A. Comoconsequência, temos que a condição de não penetração é escrita como un ≤ 0. Demodo similar, rn = R · n.
As variáveis de estado escolhidas devem ser tais que o processo definido sejatermodinamicamente admissível. Para isso é introduzida a variável interna β, quecaracteriza intensidade de adesão, uma variável de dano. β é livre para assumirquaisquer valores em [0, 1], onde β = 1 significa adesão completa, 0 ≤ β ≤ 1 significaadesão parcial, e β = 0 significa perda completa de adesão. Neste modelo, dano éconsiderado um fenômeno irreversível, assim β deve apenas diminuir.
Como estamos desconsiderando atrito e forças tangenciais, o modelo ficareduzido a contato unilateral com adesão:
rn + β2Cnun ≤ 0, un ≤ 0, (rn + β2Cnun)un = 0 (2.1a)
β = −[1b
⟨βCnu
2n − w
⟩]1/p
(2.1b)
Onde Cn é um parâmetro positivo representando a rigidez original da interfaceS, isto é, com β = 1; sendo a rigidez em qualquer momento definida por Cnβ(t)2. b érelacionada com a viscosidade da interface. w é o limite de decoesão, energia quedetermina o início da falha do adesivo. E 〈g〉 é notação de Macaulay para a funçãorampa:
〈g〉 =
g se g > 0
0 se g ≤ 0
6
2.2 Comportamentos previstos
Tabela 1: Valores de referência para análise de velocidade
wT (N ·mm) bT (N ·mm · s) Cn(N/mm) p1.025 3.25 4.9× 10−2 1
Ao avaliar a eq.(2.1a), observa-se que em compressão a condição de nãopenetração é estritamente verificada, un = 0. Porém quando a interface está sobtração, un < 0, a eq.(2.1a) impõe que a resistência elástica seja rn = −β2Cnun > 0.Se eventualmente un for tal que a energia elástica seja maior que a energia decoesão, βCnu2
n > w, a interface começa a se degradar, isto é, β começa a diminuirsegundo eq.(2.1b).
O modelo descrito pelas equações (2.1a)-(2.1b) é capaz de prever uma sériede comportamentos observáveis em testes padrão. E, ao se adotar um conjuntode parâmetros, ainda é possível obter comportamentos distintos dependendo davelocidade de carregamento.
Considere a força de reação entre as superfícies A e B, rn(t), quando o corpo Aé submetido a um deslocamento monotonicamente crescente à velocidade constante:un(t) = −v t, v > 0. Considere também os efeitos sobre a integridade da interfaceadesiva β. Assumimos como condições iniciais adesão completa (β = 1), e adotamosos valores para os parâmetros constitutivos presentes na Tabela 1. Observando-se oscomportamentos apresentados na fig. 2, podemos chegar as seguintes conclusões:
1. Um aumento na velocidade v causa um aumento tanto na força (rn) quanto navelocidade de degradação (β).
2. Um aumento do parâmetro de viscosidade b leva a uma redução na velocidadede degradação (β) e a um aumento na força (rn).
Em consideração ao valor máximo de rn(t), ao se comparar as figuras 2a e2e pode-se observar a influência clara do parâmetro de viscosidade no modelo: osistema mais viscoso demora mais para falhar e apresenta uma resistência maiorà grande velocidade deformação. A figura 2a descreve um comportamento quaseelástico durante a fase saudável da interface, onde um aumento de v tem pouquíssimainfluência sobre o valor máximo da reação. Por outro lado, da fig. 2e, para altaviscosidade sob as mesmas variações de v, observações elevada variação na força rn.
Adicionalmente, é possível alterar o expoente p para se obter diferentescomportamentos da interface com base na velocidade de degradação.
7
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
t(s)
rn(N
)
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(a) rn(t) (b = 0.1bT ).
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
t(s)
β
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(b) β(t) (b = 0.1bT ).
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
t(s)
rn(N
)
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(c) rn(t) (b = bT ).
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
t(s)
β
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(d) β(t) (b = bT ).
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
t(s)
rn(N
)
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(e) rn(t) (b = 10bT ).
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
t(s)
β
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(f) β(t) (b = 10bT )
Figura 2: Comportamentos previstos para força rn(t) (à esquerda) e intensidade de adesão β (àdireita), contra o tempo t, para diversas velocidades de deslocamento prescrito.
8
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
un(mm)
rn(N
)
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(a) rn(t) (b = 0.1bT ).
0 5 10 15 20 25−0.07
−0.06
−0.05
−0.04
−0.03
−0.02
−0.01
0
0.01
un(mm)
β
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(b) β(t) (b = 0.1bT ).
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
un(mm)
rn(N
)
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(c) rn(t) (b = bT ).
0 5 10 15 20 25−0.035
−0.03
−0.025
−0.02
−0.015
−0.01
−0.005
0
0.005
un(mm)
β
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(d) β(t) (b = bT ).
0 5 10 15 20 250
0.1
0.2
0.3
0.4
0.5
0.6
un(mm)
rn(N
)
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(e) rn(t) (b = 10bT ).
0 5 10 15 20 25−16
−14
−12
−10
−8
−6
−4
−2
0
2x 10
−3
un(mm)
β
V = 1 mm/sV = 2 mm/sV = 5 mm/sV = 10 mm/s
(f) β(t) (b = 10bT )
Figura 3: Comportamentos previstos para força rn(t) (à esquerda) e taxa de degradação de adesãoβ(t) (à direita), contra o deslocamento un, para diversas velocidades de deslocamento prescrito.
A fig. 3 apresenta o comportamento da reação do adesivo, em função dodeslocamento. É apresentada também a taxa de degradação β. Observa-se queindependentemente da velocidade, o início do dano se dá em un ≈ 4.57mm.
9
3 Análises Numéricas
Esta seção apresenta sistema físico proposto como base de um experimentode identificação dos parâmetros do modelo constitutivo apresentado na Seção 2.Uma formulação adimensional é apresentada, e um estudo de sensibilidade globaldo sistema físico proposto é realizado. Tal análise é feita através de um método deMonte Carlo, com Nmc = 1000 realizações do problema direto.
3.1 Sistema Físico Proposto
Considere o seguinte sistema mecânico, proposto em [7, 18], para a avaliaçãodas propriedades do adesivo:
Figura 4: Viga fina engastada em uma ponta e com extremidade “livre” conectada a uma superfícierígida por um adesivo. Acima: diagrama de corpo livre de uma seção da viga.
Este sistema é utilizado para gerar dados sintéticos a serem usados na calibraçãodo modelo.
A viga esbelta de comprimento L, com material homogeneamente distribuído,com seção transversal uniforme ao longo do comprimento e simétrica com relação aoplano xy. A seção transversal possui momento de inércia I. A viga é supostamentefeita de um material elástico, com coeficiente de elasticidade E. A viga é engastadaem x = 0 numa parede rígida, enquanto que a extremidade x = L é colada por umacamada fina de adesivo em um suporte rígido.
10
A fig. 4 apresenta um esquema desta construção e o diagrama de corpo livre deuma seção espacial da viga em [x, L] em um determinado instante de tempo. Ondeu(x) corresponde ao campo de deslocamentos verticais de cada ponto da viga, V (x)e M(x) são respectivamente a força cortante e o momento fletor atuantes na seçãotransversal. rn corresponde à reação entre a ponta da viga e o suporte. F é umcarregamento imposto à viga, aplicado em x = L/2.
Fazemos as seguintes hipóteses:
1. O suporte é rígido de modo que o deslocamento relativo entre este a ponta daviga coincide com o deslocamento da ponta, isto é, un(L, t) = −u(L, t).
2. As dimensões da interface são desprezíveis em comparação com o comprimentoda viga. O que nos permite um modelo local de adesão, isto é, apenas a pontada viga (x = L) está sujeita à restrição imposta pelas equações (2.1a)-(2.1b).
3. O histórico de carregamento F (t) é quasi estático.
4. A viga mantém seu comportamento elástico linear ao longo de todocarregamento.
Nestas condições, a relação entre a curvatura da viga e o momento fletor M(x),definida pela hipótese de Euler-Bernoulli, é:
EId2uydx2 = −rn(L− x) + F 〈L/2− x〉 , 0 ≤ x ≤ L (3.1a)
uy(0) = 0, duydx
(0) = 0 (3.1b)
Integrando (3.1a) sujeita as condições de contorno (3.1b), temos:
u(x, t) = L3
48EI
{−8rn(t)
(3x
2
L2 −x3
L3
)+ F (t)
[⟨1− 2x
L
⟩3+ 6x
L− 1
]}
Para x = L:u(L, t) = L3
3EI
[−rn(t) + 5
16F (t)]
(3.2)
Impomos a condição de adesão, eq. (2.1a), em x = L:
rn(t)− β2Cnu(L, t) ≤ 0, u(L, t) ≥ 0,(β2Cnu(L, t)− rn(t)
)u(L, t) = 0 (3.3)
A partir de (3.2) e (3.3) é possível chegar uma relação fechada entre a reaçãode contato rn(t), deslocamento da ponta u(L, t) e o carregamento aplicado F (t):
11
rn(t) = 5F16
1 se F (t) ≤ 0
β2CnL3
3EI(1 + β2CnL3
3EI
) se F (t) > 0(3.4)
u(L, t) = 5F16
0 se F (t) ≤ 0
L3
3EI(1 + β2CnL3
3EI
) se F (t) > 0(3.5)
Portanto, para o sistema da fig. 4, o problema direto consistem em determinaro deslocamento u(t) dado: um carregamento prescrito F (t), as características daviga [E, I, L], e um conjunto de parâmetros do modelo de adesão, [w,Cn, b, p]. Parao problema inverso, nos interessa determinar [w,Cn, b, p] dado um conjunto deobservações y e as características da viga.
3.2 Formulação Adimensional
Considere que a rigidez EI/L3 da viga é conhecida e perfeitamente definida,de modo que é possível eliminar sua influência na calibração do modelo definindouma variável de escala. E por simplicidade, considere que o carregamento F (t) podeser escrito como F (t) = F0f(t), onde F0 é um parâmetro de escala. Sendo assim,definimos as variáveis adimensionais:
u(t) = 3u(L, t)EIF0L3 , r(t) = rn(t)
F0(3.6)
As quais fornecem as propriedades adimensionais do adesivo:
c = CnL3
3EI , w = 3wEIF0L3 , b = 3bEI
F0L3 (3.7)
Portanto é possível reescrever eq. (2.1b) em termos das propriedadesadimensionais:
β = −[1b
⟨βcu2 − w
⟩]1/p
com 0 ≤ β ≤ 1 (3.8)
Para o processo de calibração, basta considerar um carregamentomonotonicamente crescente f(t), uma vez que o dano age apenas sob tração. Nessascondições reescrevemos (3.4) e (3.5) como:
[r , u] = 5f(t)16 (1 + β2c)
[β2c , 1
]com f(t) > 0 (3.9)
12
3.3 Análise Global de Sensibilidade
Para nossa análise, consideramos os valores de referência apresentados naTabela 2. Estes parâmetros foram escolhidos de modo a reproduzir uma fase inicialelástica seguida de perda de adesão. Quanto à força aplicada sobre a viga na fig. 4,utilizamos dois históricos de carregamento, apresentados na fig. 5.
Tabela 2: Variáveis adimensionais de referência para análise global e estimação de parâmetros.
wT bT cT pT1.025 3.25 4.9× 10−2 1
0 20 40 60 80 1000
20
40
60
80
100
120
t(s)
f(t)
F
1
F2
Figura 5: Histórico de carregamento
A análise de sensibilidade é realizada modelando-se os parâmetros constitutivosdo adesivo como variáveis aleatórias cuja função densidade de probabilidade π(θ)é conhecida. São extraídas Nmc = 1000 realizações dos parâmetros constitutivospara cada um dos cenários apresentados, θi ∼ π(θ), e simulações de resposta dosistema aos carregamentos da fig. 5 são geradas para cada uma dessas realizações,{r1(θ1, t), ..., rN (θNmc , t)}. Os cenários construídos com base nos valores de referênciada Tabela 2 são:
1. Cenário S1: c ∼ Nc(cT , (0.1 cT )2), mantendo os outros parâmetros fixos eiguais aos de referência.
2. Cenário S2: w ∼ Nw(wT , (0.1 wT )2), mantendo os outros parâmetros fixos eiguais aos de referência.
13
3. Cenário S3: b∗ ∼ Nb∗(b∗, (0.1 b∗)2), mantendo os outros parâmetros fixos eiguais aos de referência.
4. Cenário S4: Todas as variáveis acima, assumindo que b∗, w e c são variáveisaleatórias independentes de distribuição gaussiana.
Os quatro cenários são avaliados para b∗ = 0.1bT , b∗ = bT e b∗ = 10bT , paraambos os carregamentos apresentados na fig. 5. O interesse de se utilizar doiscarregamentos está em podermos avaliar a influência na incerteza dos parâmetrosem diferentes velocidades de deslocamento.
Os resultados para a reação r segundo sensibilidade aos parâmetros sãoapresentados nas figuras de 7 a 10, páginas 15 a 18: a resposta média do sistema e oenvelope de resposta definido pelos percentis de 1% e 99%. Por simplicidade, na fig.6 apresenta-se a influência que a variação simultânea em todos parâmetros produzno deslocamento, para o caso de alta viscosidade; os resultados para os outros casossão semelhantes, e apresentam envelope de resposta mais estreito. Portanto, épossível concluir-se que o deslocamento u é muito pouco sensível a variações nosparâmetros para o modelo considerado.
Figura 6: Sensibilidade do deslocamento para variações na rigidez, com b = 10 bT e p = 1.
14
(a) F1 - b∗ = 0.1 bT (b) F2 - b∗ = 0.1 bT
(c) F1 - b∗ = bT (d) F2 - b∗ = bT
(e) F1 - b∗ = 10 bT (f) F2 - b∗ = 10 bT
Figura 7: Cenário S1: Influência da rigidez.
15
(a) F1 - b∗ = 0.1 bT (b) F2 - b∗ = 0.1 bT
(c) F1 - b∗ = bT (d) F2 - b∗ = bT
(e) F1 - b∗ = 10 bT (f) F2 - b∗ = 10 bT
Figura 8: Cenário S2: Influência da energia.
16
(a) F1 - b∗ = 0.1 bT (b) F2 - b∗ = 0.1 bT
(c) F1 - b∗ = bT (d) F2 - b∗ = bT
(e) F1 - b∗ = 10 bT (f) F2 - b∗ = 10 bT
Figura 9: Cenário S3: Influência da viscosidade.
17
(a) F1 - b∗ = 0.1 bT (b) F2 - b∗ = 0.1 bT
(c) F1 - b∗ = bT (d) F2 - b∗ = bT
(e) F1 - b∗ = 10 bT (f) F2 - b∗ = 10 bT
Figura 10: Cenário S4: Influência de todos os parâmetros.
18
0.2 0.25 0.3 0.35 0.40
5
10
15
20
25
30
35
40
rmax
F1 − 0.1 bT
F2 − 0.1 bT
F1 − 10 bT
F2 − 10 bT
(a) π(rmax) para Cenário S1
0.18 0.2 0.22 0.24 0.26 0.28 0.30
20
40
60
80
100
120
140
160
180
200
rmax
F1 − 0.1 bT
F2 − 0.1 bT
F1 − 10 bT
F2 − 10 bT
(b) π(rmax) para Cenário S2
0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30
200
400
600
800
1000
1200
1400
rmax
F1 − 0.1 bT
F2 − 0.1 bT
F1 − 10 bT
F2 − 10 bT
(c) π(rmax) para Cenário S3
0.2 0.25 0.3 0.35 0.40
5
10
15
20
25
rmax
F1 − 0.1 bT
F2 − 0.1 bT
F1 − 10 bT
F2 − 10 bT
(d) π(rmax) para Cenário S4
Figura 11: Função Densidade de Probabilidade da reação máxima rmáx, para carregamentos F1 eF2, para os diversos cenários.
Tabela 3: Média µr e desvio padrão σr da reação máxima. rmáx.
Cenário 0.1 bT bT 10 bTµr σr µr σr µr σr
F1
S1 0.2244 0.0116 0.2298 0.0112 0.2662 0.0169S2 0.2246 0.0113 0.2297 0.0110 0.2669 0.0074S3 0.2246 0.0003 0.2310 0.0012 0.2666 0.0027S4 0.2242 0.0158 0.2304 0.0160 0.2653 0.0194
F2
S1 0.2265 0.0113 0.2448 0.0128 0.2804 0.0235S2 0.2265 0.0110 0.2454 0.0095 0.2794 0.0043S3 0.2269 0.0005 0.2446 0.0018 0.2832 0.0046S4 0.2258 0.0157 0.2454 0.0172 0.2831 0.0252
19
80 82 84 86 88 90 920
0.5
1
1.5
2
2.5
3
3.5
4
tβ≈0.2
cwbc, w, b
(a) π(T0.2) para F1 - b∗ = 0.1 bT
40.6 40.8 41 41.2 41.4 41.6 41.8 42 42.2 42.40
5
10
15
tβ≈0.2
cwbc, w, b
(b) π(T0.2) para F2 - b∗ = 0.1 bT
83 84 85 86 87 88 89 90 91 920
0.5
1
1.5
2
2.5
3
3.5
4
4.5
tβ≈0.2
cwbc, w, b
(c) π(T0.2) para F1 - b∗ = bT
40.5 41 41.5 42 42.5 43 43.50
0.5
1
1.5
2
2.5
3
3.5
4
tβ≈0.2
cwbc, w, b
(d) π(T0.2) para F2 - b∗ = bT
86 88 90 92 94 96 98 1000
0.2
0.4
0.6
0.8
1
1.2
1.4
tβ≈0.2
cwbc, w, b
(e) π(T0.2) para F1 - b∗ = 10 bT
43 43.5 44 44.5 45 45.5 46 46.5 470
0.5
1
1.5
2
2.5
3
3.5
tβ≈0.2
cwbc, w, b
(f) π(T0.2) para F2 - b∗ = 10 bT
Figura 12: Função Densidade de Probabilidade do tempo de falha T0.2, para carregamentos F1 eF2, para os diversos cenários.
20
Tabela 4: Média µT e desvio padrão σT do tempo T0.2.
Cenário 0.1 bT bT 10 bTµT σT µT σT µT σT
F1
S1 85.521 0.867 87.121 0.858 92.713 1.141S2 86.445 0.838 87.184 0.687 92.545 0.639S3 86.443 0.058 87.018 0.117 92.528 0.470S4 86.505 1.164 87.144 1.050 92.727 1.500
F2
S1 41.318 0.165 41.986 0.196 44.877 0.387S2 41.320 0.160 41.968 0.138 44.858 0.190S3 41.302 0.015 41.964 0.065 44.823 0.204S4 41.324 0.215 41.966 0.250 44.845 0.500
Os resultados das análises apresentados nas figuras 7-10 estão condensadosnas Tabelas 3 e 4. As principais diferenças entre os carregamentos F1 e F2 são omáximo da força de reação r, rmáx = máx{r(t)}, e o tempo de falha “total” doadesivo; sendo o segundo essencialmente diferente pelas diferentes velocidades dosdois carregamentos. Definimos o tempo T0.2 tal que β(T0.2) = 0.2, equivalente a umdano acumulado de 80%, como o tempo de falha “parcial” do adesivo.
Nas figuras 11 e 12, são apresentadas as distribuições resultantes nestas duasvariáveis para os cenários analisados. Na fig. 11, são omitidas as PDFs para o casode viscosidade mediana (b∗ = bT ) em vista da similaridade com o caso de baixaviscosidade (b∗ = 0.1 bT ), de modo a que a figura ficasse mais clara.
Com base nos valores apresentados na Tabela 3, com relação a rmax, observamos(que):
• Aumento monótono da variância da reação máxima rmáx com o aumento daviscosidade, para todos os cenários e carregamentos considerados.
• Um aumento na taxa de aplicação da força, equivalente à imposição de umdeslocamento com maior velocidade, gera uma redução na variância de rmáx,quando variações na energia w são consideradas.
• Para ambos os carregamentos, variações na viscosidade para os casos de baixa(b∗ = 0.1 bT ) e média (b∗ = bT ) viscosidade apresentam desvios padrão similares;a diferença apenas sendo acentuada para o carregamento F2.
21
Vale notar também comparando-se os casos de baixa (b∗ = 0.1 bT ) e alta(b∗ = 10 bT ) viscosidade, se observa o aparecimento de um período de relaxamento noqual o adesivo se acomoda à força nele aplicada pela viga. Assim como o surgimentode um segundo pico de resistência (em t ≈ 85s para F1, e t ≈ 42s para F2) tão maiorquanto maior for a taxa de aplicação da força.
Com base nos valores apresentados na Tabela 4, com relação a T0.2, observamos(que):
• Embora os tempos médios µT sejam muito similares entre si para os diversoscenários, consistentemente os cenários S1 e S2 apresentam maior desvio padrão.
• Um aumento na taxa de aplicação da força reduz a variância em todos oscenários considerados.
A análise global sugere que para modelos com p = 1 os parâmetros maisinfluentes em todos os regimes de viscosidade considerados são a rigidez c e energiade decoesão w — uma pequena variação em seus valores muda drasticamente aresposta do sistema. Em especial há uma alteração do valor pico de resistência, oqual geralmente é uma variável de projeto importante.
22
4 Filtro de Kalman
4.1 O Problema
Sejam {xk}∞k=0 e {yk}∞k=1 dois processos estocásticos. Os vetores colunaaleatórios xk ∈ Rm×1 e yk ∈ Rn×1 são respetivamente os vetores de estados eobservação no instante k. Dadas as distribuições de probabilidade absolutamentecontínuas, postulamos o seguinte:
1. O processo {xk}∞k=0 é um processo de Markov, isto é,
p(xk+1|x0, ...,xk) = p(xk+1|xk)
2. O processo {yk}∞k=1 é um processo de Markov com respeito a {xk}, isto é,
p(yk|x0, ...,xk) = p(yk|xk)
3. O processo {xk}∞k=0 é dependente das observações passadas apenas pelo seupróprio histórico, isto é,
p(xk+1|xk,y1, ...,yk) = p(xk+1|xk)
Se as condições acima são satisfeitas, dizemos que {xk}∞k=0 e {yk}∞k=1 formamum modelo de evolução e observação. Evidentemente, para o modelo ser completo, énecessário especificar:
1. A função de densidade de probabilidade do estado inicial x0.
2. As funções de transição de Markov p(xk+1|xk), k = 0, 1, 2, ...
3. As funções de verossimilhança p(yk|xk), k = 1, 2, ...
As funções de transição p(xk+1|xk) podem variar com o tempo, assim como asfunções p(yk|xk). Assuma agora que temos um modelo de Markov descrevendo aevolução dos estados xk e um modelo de observação para yk, para k = 0, 1, 2, ...
xk+1 = f(xk,wk) (4.1a)
yk+1 = h(xk+1,vk+1) (4.1b)
Onde f : Rm×1 → Rm×1 é uma função conhecida, possivelmente não-linear, doestado xk e do ruído de processo wk, é chamada de função de evolução de estado.
23
Enquanto que h : Rm×1 → Rn×1 é chamada de função de observação, uma função dexk+1 e de um ruído de medição vk+1, também possivelmente não-linear. Fazemosagora as seguintes hipóteses com relação aos ruídos:
1. Para k 6= l, os vetores de ruído wk e wl, assim como vk e vl, são mutuamenteindependentes da condição inicial x0.
2. Os vetores wk e vl são mutuamente independentes para todos k e l.
Procuramos estimativas do estado xk baseado na sequência de todas asobservações Dk = {yi}ki=1 feitas até o momento tk. Do ponto de vista Bayesiano,queremos encontrar a Função Densidade de Probabilidade (PDF -ProbabilityDensity Function) p(xk|yk) recursivamente em duas etapas:
Previsão:Dado p(xk|Dk), encontrar p(xk+1|Dk) baseado na função de transição de Markovp(xk+1|xk). Temos:
p(xk+1|Dk) =∫
Xk
p(xk+1|xk)p(xk|Dk)dxk (4.2)
Correção:Dado p(xk+1|Dk), encontrar p(xk+1|Dk+1) baseado na nova observação yk+1 ena função de verossimilhança p(yk+1|xk+1). Temos:
p(xk+1|Dk+1) = p(yk+1|xk+1)p(xk+1|Dk)p(yk+1|Dk)
(4.3)
A prova dessas relações não faz parte deste trabalho, mas pode ser encontradaem detalhes na referência [19]. A Eq. (4.2) é a equação de Chapman-Kolmogorov,representando a distribuição de probabilidade marginal de xk+1, enquanto que (4.3)é simplesmente a fórmula de Bayes. Estas duas equações formam a base de umasolução ótima para a determinação da densidade de probabilidade recursivamente.
24
4.2 Filtro de Kalman Discreto
No caso linear e gaussiano, equações (4.2) e (4.3) formam as equações do Filtrode Kalman. Assumindo que as PDF são sempre gaussianas, estas são caracterizadaspor apenas dois parâmetros, a média µ e a matriz de covariância Σ. Dizer x ∼ N (µ,Σ)é dizer que a função densidade de probabilidade de x é:
PDF (x) = ‖2πΣ‖− 12 exp
(−1
2(x− µ)T Σ−1(x− µ))
É possível mostrar que se p(xk−1|Dk−1) é gaussiana, então p(xk|Dk) também oé. Para isso, as seguintes condições devem ser atendidas:
1. wk e vk+1 são sequências mutuamente independentes de ruído branco gaussianode média zero, com E[wkwT
k ] = Qk, E[vkvTk ] = Rk e E[viwTj ] = 0, ∀i, j.
2. f(xk,wk) é uma função linear.
3. h(xk+1,vk+1) é uma função linear.
Isto é,xk+1 = Fk xk + wk (4.4a)
yk+1 = Hk+1xk+1 + vk+1 (4.4b)
onde wk é o ruído de processo com covariância Qk ∈ Rm×m, e vk+1 é o ruído demedição com covariância Rk+1 ∈ Rn×n. Note que as matrizes Fk, Hk+1, Qk e Rk+1
estão livres para variar com o tempo, isto é, Fk(tk), Hk+1(tk+1), etc.É comum em controle haver um termo adicional Bk uk, composto de um vetor
de controle uk ∈ R p×1 e uma mapeamento Bk ∈ Rm×p. Neste caso, reescrevemosequação (4.4a) como:
xk+1 = Fk xk +Bkuk (4.5)
E modificamos (4.6) de acordo.Denotando agora xk|l = E(xk|Dl) e Pk|l = cov(xk|Dl), e definindo p(x0) =
p(x0|D0), continuamos com o desenvolvimento:
Previsão: Assumindo conhecida a distribuição gaussiana
p(xk|Dk) ∼ N (xk|k, Pk|k)
Entãop(xk+1|Dk) ∼ N (xk+1|k, Pk+1|k)
25
Ondexk+1|k = Fk xk|k (4.6)
Pk+1|k = FkPk|kFTk +Qk (4.7)
Correção: Assumindo conhecida a distribuição gaussiana
p(xk+1|Dk) ∼ N (xk+1|k, Pk+1|k)
Entãop(xk+1|Dk+1) ∼ N (xk+1|k+1, Pk+1|k+1)
OndeK = Pk+1|kHk+1(Hk+1Pk+1|kH
Tk+1 +Rk+1)−1 (4.8)
xk+1|k+1 = xk+1|k +K (yk+1 −Hk+1xk+1|k) (4.9)
Pk+1|k+1 = (1−KHk+1)Pk+1|k (4.10)
As equações (4.6)-(4.10) formam o algoritmo padrão do filtro de Kalman.Iniciamos o processo iterativo fornecendo x0|0, a condição inicial, e uma medida dequão confiantes estamos neste valor, sua matriz de covariância P0|0. Os aspectospráticos da filtragem têm suas sutilezas, tais como a forma de fatoração matricialadequada, valores fornecidos na inicialização do filtro, dentre outras coisas.
A determinação da matrizes de covariância Qk e Rk é um problema à parte.Para a maioria das aplicações, as matrizes de covariância são determinadas em umafase laboratorial a partir de ensaios ou simulações [12, 13], mas em outros casosconhecimento a priori destas é difícil de se obter, ou simplesmente o ajuste (tuning)ótimo do filtro pode variar com o tempo. Nestas situações, um esquema adaptativopode ser utilizado para estimar Qk e/ou Rk conforme o filtro opera [14, 20].
26
4.3 Filtro de Kalman Estendido Discreto
O desenvolvimento do filtro apresentado na seção anterior assumia um sistemalinear. Em realidade esse quase nunca é o caso, e somos forçados a utilizar algumasaproximações. Ao longo dos anos diversas aproximações foram propostas, baseadasem métodos analíticos (EKF, EKF de 2a ordem), numéricos ou de amostragem (UKF,Filtros de Partículas).
A consequência prática mais relevante de se considerar o sistema linear é que aPDF após a propagação do vetor de estado pelo modelo dinâmico (4.4a). No casonão-linear isso não é realista, assim como não é realista imaginar que (4.7) descrevasua covariância. Além disso outros momentos seriam necessários para descrever aPDF, uma vez que esta não é mais gaussiana.
O Filtro de Kalman Estendido (EKF - Extended Kalman Filter) é derivadopara sistemas não-lineares com ruído aditivo, isto é:
xk+1 = f(xk) + wk (4.11a)
yk+1 = h(xk+1) + vk+1 (4.11b)
vk+1 e wk são sequências aleatórias mutuamente independentes de ruídosgaussianos, brancos de média nula e covariâncias Rk+1 e Qk, respectivamente. Asfunções (4.11a) e (4.11b) são aproximadas fazendo-se uma expansão em série deTaylor de 1a ordem. Assumimos então que esta linearização local é suficiente paradescrever as não-linearidades envolvidas.
O EKF é o seguinte algoritmo formado por (4.12)-(4.18):
Previsão:xk+1|k = f(xk|k) (4.12)
Fk = ∂f
∂x
∣∣∣∣∣xk|k
(4.13)
Pk+1|k = Fk Pk|k FTk +Qk (4.14)
Correção:
Hk+1 = ∂h
∂x
∣∣∣∣∣xk+1|k
(4.15)
K = Pk+1|kHk+1 (Hk+1 Pk+1|kHTk+1 +Rk+1)−1 (4.16)
xk+1|k+1 = xk+1|k +K (yk+1 − h(xk+1|k)) (4.17)
Pk+1|k+1 = (1−KHk+1)Pk+1|k (4.18)
27
Fk e Hk+1 são as matrizes jacobianas de f e h, avaliadas em xk|k e xk+1|k,respectivamente. A matriz Jacobiana de uma função vetorial f : Rm → Rn é a matrizJn×m formada pelas suas derivadas parciais.
Ji,j = ∂fi∂xj
=
∂f1
∂x1· · · ∂f1
∂xm... . . . ...∂fn∂x1
· · · ∂fn∂xm
(4.19)
Para algumas funções, avaliar o jacobiano é computacionalmente custoso sefeito numericamente, ou simplesmente difícil (ou impossível) de ser calculadoanaliticamente. Isto e o fato conhecido de haver instabilidade, até divergência, emcondições de não-linearidade elevada e modelagem incerta do sistema [21] são asprincipais desvantagens desta formulação, e levaram ao desenvolvimento diversosesquemas para superar estas dificuldades [14, 22, 23].
4.4 EKF Contínuo-Discreto
Os valores de Rk e Qk utilizados no filtro dependem da forma do sistema e dadiscretização. Se estivéssemos lidando com um sistema linear invariante no tempo,por exemplo:
x = Ax(t) +Bu(t) + w(t), w(t) ∼ N (0, Q(t))
y(t) = Cx(t) +Du(t) + v(t), v(t) ∼ N (0, R(t))
Poderia ser facilmente discretizado por Ordem Zero (zero-order hold) para:
xk+1 = Akxk +Bkuk + wk,wk ∼ N (0, Qk)
yk = Cxk +Buk + vk,vk ∼ N (0, Rk)
Onde ∆T é o intervalo de tempo de discretização e
Ak = exp {A∆t}
Bk = A−1(Ak − 1)B, se det(A) 6= 0
Qk =∆t∫0
eAτQeτ AT
dτ
Rk = R
∆t
28
A partir disso é comum se fazer as seguintes aproximações Rk ≈ R/∆t eQk ≈ Q∆t, assim como aproximações para Ad, mesmo para sistemas que não sejamlineares invariantes, pela dificuldade que é estimar estas matrizes.
Muitas das vezes os sistema de engenharia são formulados no domínio dotempo, de modo contínuo, e posteriormente passam por uma discretização, enquantoque medições são tomadas em intervalos definidos de tempo. O Espaço de Estadoscorrespondente é apresentado em (4.20a)-(4.20b):
x = f(t,x) + w(t), w(t) ∼ N (0, Q(t)) (4.20a)
yk = h(tk,xk) + vk, vk ∼ N (0, Rk) (4.20b)
Modificamos o algoritmo do EKF para melhor se encaixar nesta formulação,também é conhecida como “híbrida”.. Note que a etapa de previsão é feitaresolvendo-se duas equações diferenciais (4.21) e (4.23), cujos métodos de soluçãoagora estão desvinculados do algoritmo do filtro, possibilitando a utilizando deesquemas mais robustos de integração. Equação (4.23) é a equação matricial deRiccati. Note também que não há necessidade de corrigir Q pela discretização, umavez que este está dentro da integração.
Em busca de uma notação mais simplificada, omitimos o vetor de controle ucomo argumento da função f(t,x,u). Mas vale ressaltar que isso em nada interfereno cálculo do jacobiano em (4.22) se lembrarmos de avaliá-lo em [t,x(t),u(t)].
O algoritmo utilizado neste trabalho é formado por (4.21)-(4.30). Algoritmossimilares que levam em consideração a formulção em EE da forma (4.20a)-(4.20b)podem ser encontrados nas referências [15, 16, 17].
Previsão:xk+1|k = xk|k +
∫ tk+1
tk
f(t,x) dt (4.21)
F (t) = ∂f
∂x
∣∣∣∣∣t,x(t)
(4.22)
Pk+1|k = Pk|k +∫ tk+1
tk
(FP + PF T +Q) dt (4.23)
Correção:ek+1 = yk+1 − h(tk+1,xk+1|k) (4.24)
Hk+1 = ∂h
∂x
∣∣∣∣∣tk+1,xk+1|k
(4.25)
K = Pk+1|kHk+1 (Hk+1 Pk+1|kHTk+1 +Rk+1)−1 (4.26)
29
xk+1|k+1 = xk+1|k +K ek+1 (4.27)
Pk+1|k+1 = (1−KHk+1)Pk+1|k (4.28)
Adaptação:
Cν = 1n
k+1∑j=(k+1)−n
ejeTj (4.29)
Q = KCνKT (4.30)
Embora, segundo a eq. (4.20a), Q(t) esteja livre para variar, assumimos queQ(t) = Qc,k, um valor constante ao longo da integração na eq. (4.23) que é corrigidoa cada iteração, eq. (4.30).
O processo adaptativo aqui apresentado segue um esquema baseado eminovações — simplesmente a diferença entre a medição real e a esperada, eq. (4.24)— (IAE - Innovation-based Adaptative Estimation) [14, 20]. Ao se fazer IAE énecessário fornecer um valor para Q para primeiras n iterações, enquanto a série deestatísticas, eq. (4.29), está sendo calculada. Um n pequeno pode gerar umaestimação instável, devido ao ruído na medição, entretanto um n elevado causa umatraso da inicialização do processo de adaptação e também uma redução dasensibilidade a mudanças de comportamento do sistema, em outras palavrascomprometendo a capacidade de adaptação. A escolha equivocada de um Q inicialou um n pode não produzir resultados melhores que o algoritmo sem adaptação.
Escolheu-se fazer a correção apenas da matriz de covariância do processo Q,pois valores muito baixos de Rk+1 podem levar à instabilidade numérica em umambiente muito ruidoso.
A adaptação de Rk+1 pode ser feita antes da eq. (4.26) por:
Rk+1 = Cν −Hk+1 Pk+1|kHTk+1
Entretanto, isto não garante que Rk+1 é positiva definida devido à subtração,o que pode facilmente causar divergência do filtro. Foi proposto [24] então umaadaptação com base em resíduos e — a diferença entre a medição real e a mediçãodo valor estimado do estado, ao invés do valor previsto:
e = yk+1 − h(tk+1,xk+1|k+1)
Cν = 1n
k+1∑j=(k+1)−n
ej eTj
Rk+1 = Cν +Hk+1 Pk+1|kHTk+1
30
4.5 Estimação Conjunta de Parâmetros
Com o objetivo de determinar os parâmetros θ de uma função f(t,x,θ), épossível fazer uso de estados aumentados. Para este fim, supomos que a variação deθ é nula, ou muito mais lenta que a variação dos estados [15, 16, 22]. Essa suposiçãoelimina a necessidade de um modelo para propagação dos parâmetros (como é o casoda variável interna β).
Portanto a propagação dos parâmetros assume simplesmente que θ ≈ 0 (oualternativamente em uma formulação discreta, θk+1 ≈ θk). Nesse sentido, é possívelconsiderar ou não um random walk do parâmetro na implementação [15, 16, 22]:θ = a (ou θk+1 = θk + ak), onde a (ou ak) é uma variável aleatória de distribuiçãogaussiana, de média zero e covariância baixa.
Ao se trabalhar com estados aumentados, é necessário reescrever a função deevolução f(t,x,θ) como f(t, z), e a função de observação h(t,x,θ) como h(t, z).Onde z = [x,θ]T é o vetor de estados aumentados. Deste modo a identificação dosparâmetros é obtida aplicando-se o EKF ao sistema descrito pela variável z. O filtroentão produz uma estimativa dos estados e parâmetros de forma recursiva, tratandoas medições de forma sequencial e “melhorando” a estimativa em cada instante.
Esta formulação não é a melhor. Desconsiderando-se os problemas tradicionaisrelacionados ao EKF, esta formulação é conhecida por fornecer resultados inferiores,sendo mais sensível a não linearidades, ruído e a variação das condições iniciais queuma formulação “desacoplada” (DKF - Dual Extended Kalman Filter) [9], razãopela qual a formulação de estados aumentados também costuma ser chamar de“joint”. Digressões a parte, este é um método de fácil implementação, o qual nãonecessita modificar o algoritmo do filtro.
31
5 Estimação de Parâmetros
Nesta seção, dados ideais yT são gerados com os parâmetros adimensionais daTable 2 (wT = 1.025, cT = 4.9× 10−2, pT = 1). Em seguida, estes são poluídos comruído branco gaussiano ε, de variância σ2
ε = 0.025 max(yT ), para formar um vetorde dados corrompidos y = yT + ε, com SNR ≈ 20. Considera-se que medições de re u são tomadas ambas com uma frequência de 20 Hz. São analisados os casos paraos quais b∗ = 0.1 b T , b∗ = b T e b∗ = 10 b T , sob o carregamento F2 apresentado nafig. 5.
5.1 Formulação em Espaço de Estados
Porém, antes de atacar o problema da estimação dos parâmetros, é precisocolocar as equações (3.9) numa forma em que possam ser utilizadas no algoritmorecursivo. Para tal, estas equações são derivadas em função do tempo, para quea posterior integração na eq. (4.21) introduza uma condição inicial, que fornece aparte recursiva do processo.
Seja f = f(t) o carregamento, e f = df/dt sua derivada. Então:
du
dt= 5f(t)
16cβ2 + 16 −160cβf(t)β
(16cβ2 + 16)2 (5.1)
dr
dt= 5cβ2f(t)
16cβ2 + 16 −160c2β3f(t)β(16cβ2 + 16)2 + 10cβf(t)β
16cβ2 + 16 (5.2)
Portanto, nossa função de evolução de estados é dada pelas equações (5.1),(5.2)e (3.8), e a função de observação é simplesmente da forma yk = Cxk. Em espaço deestados:
x =
˙u˙rβ
=
5f(t)16cβ2 + 16 −
160cβf(t)β(16cβ2 + 16)2
5cβ2f(t) + 10cβf(t)β16cβ2 + 16 − 160c2β3f(t)β
(16cβ2 + 16)2
−[1b
⟨βcu2 − w
⟩]1/p
yk =1 0 0
0 1 0
x(tk)
32
Nesta seção, foi escolhido estimar o vetor de parâmetros θ = [g, w] T . Ondeé feito g = b−1, para evitar a introdução de não-linearidades desnecessárias para ofiltro. Fazendo uso da formulação de estados aumentados, o espaço de estados queentra no algoritmo (4.21)-(4.30) é dado por:
z =
˙u˙rβ
˙g˙w
=
5f(t)16cβ2 + 16 −
160cβf(t)β(16cβ2 + 16)2
5cβ2f(t) + 10cβf(t)β16cβ2 + 16 − 160c2β3f(t)β
(16cβ2 + 16)2
−[g⟨βcu2 − w
⟩]1/p
νg
νw
(5.3a)
yk =1 0 0 0 0
0 1 0 0 0
z(tk) (5.3b)
Onde νg ∼ N (0, (0.05 g(t))2) e νw ∼ N (0, (0.05 w(t))2).
33
5.2 Influência das Condições Iniciais
Aqui os parâmetro p e c são utilizados como parâmetros de controle para acalibração do modelo, embora por motivos distintos: a estimação conjunta de p e θ
se mostrou problemática; enquanto que c é um parâmetro que espera-se pode serestimado com facilidade durante um simples ensaio de tração feito numa etapa inicial,e assim foi deixado “de lado” com o intuito de facilitar a apresentação visual dosresultados. Lembrando que é sempre possível incluir a estimação de c, e fornecer aomesmo tempo um bom chute inicial para este parâmetro segundo a breve discussãodeste parágrafo.
Escolheu-se utilizar o carregamento F2, apresentado na fig. 5, para o processode estimação. Para a inicialização do filtro, fazemos:
1. A condição inicial z0|0 = [x0|0,θ0|0] T , onde x0|0 = [0, 0, 1] T . Recorde dasequações (3.4) e (3.5) que para F (t) ≤ 0 temos rn(t) = 5F/16 e u(L, t) = 0.Em nosso carregamento F (0) = 0, e assumimos que o adesivo está em perfeitascondições, β = 1.
2. A matriz de covariância de processo Q = diag( [1, 1, 10−2, 25, 25] T ). Estevalor é esperado mudar segundo eq. (4.30).
3. A matriz de covariância diagonal P0|0 = Q. Uma P0|0 diagonal significaassumir que as condições iniciais não são relacionadas. As entradas da matrizcorrespondentes ao parâmetros a estimar possuem variância elevada de modoa que o filtro rapidamente abandone as condições iniciais e efetue as correçõesnecessárias. De modo análogo, o valor inicial de β é assumido conhecido,portanto associado a uma variância bem baixa.
4. A matriz de covariância de medição R = 1.5 diag(σ2ε). O ruído introduzido
para gerar as medições de r é relativamente menor, portanto associado a umvalor menor em R. Estes valores são 50% maiores que as variâncias dos ruídosintroduzidos para gerar as medições sintéticas para estabilizar o filtro contrapossíveis valores elevados de ruído.
5. No que diz respeito à estimação de parâmetros, considera-se a influência dacondição inicial θ0|0 fazendo w e g variarem no domínio D = {(w, g) ∈ R2|0.5 ≤w0/wT ≤ 1.5, 0.5 ≤ g0/gT ≤ 1.5}.
6. Consideram-se os casos de baixa (b = 0.1 bT ), média (b = bT ) e altaviscosidade (b = 10 bT ). Fazendo a estimação para c conhecido e igual a:0.9 c T , c T , 1.1 c T .
34
Os resultados da estimação de parâmetros para as condições iniciais em D sãoapresentados nas figuras 13, 14 e 15. Nos eixos estão a razão entre a condição inicialdos parâmetros e seus valores true — g 0/g T × w 0/w T —, e é plotada a relaçãoentre a média das últimas i = 200 iterações do filtro e seus valores true — w∗/w T (àesquerda) e g∗/g T (à direita).
A fig. 13 mostra que, para o regime de baixa viscosidade (b = 0.1 bT ), é possívelestimar consistentemente a energia w dado que o chute inicial seja menor ou igualao valor real do parâmetro. A região onde é possível se estimar w se estende para adireita (maiores valores para w0) conforme a rigidez aumenta. Entretanto, a regiãoonde se é possível estimar g permanece estreita; nesta região a tendência é subestimaro valor de g, ou seja, sobrestima b.
A fig. 14 mostra que, para o regime de viscosidade mediana (b = bT ), é possívelestimar a energia w para uma faixa ampla de condições iniciais somente quandoc = 1.1 cT , similar ao caso apresentado na fig. 13e. Nesta mesma região, os resultadospara g se mostraram melhores que para o caso de baixa viscosidade. Porém a regiãoem D onde se é possível obter estimação satisfatória de ambos os parâmetros ficareduzida no caso em que c = cT paraD∗ = {(w, g) ∈ R2|0.8 wT ≤ w0 ≤ 1.3 wT , gT ≤g0 ≤ 1.5 gT , }, como mostram as figuras 14c e 14d.
A fig. 15 mostra que, para o regime de alta viscosidade (b = 10 bT ), não épossível identificar simultaneamente w e g. A figura não possui indicação das curvasde nível para não ficar muito poluída visualmente.
35
0
0
0
0
0.2
0.2
0.2
0.4
0.4
0.4
0.6
0.6
0.6
0.8
0.8
0.8
0.20.21
1
1
1 0.40.20.40.6
0
0
0.8
1.2
1
1
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(a) w∗/wT (c = 0.9 cT )
0
0.2
0.2
0.2
0.4
0.4
0.4
0.4
0.4
0.6
0.6
0.6
0.6
0.8
0.8
0.8
0.8
0.2
0.80.8
0.81
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(b) g∗/gT (c = 0.9 cT )0
0
0
0
0.2
0.2
0.2
0.2
0.4
0.4
0.4
0.6
0.6
0.6
0.8
0.8
0.8
1
1
1
1
1
1
0.4
0.4
0.2
1
0.6
1.2
1
1
0.8
0.4
1 0.6
1
1
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(c) w∗/wT (c = cT )
0
0.2
0.2
0.2
0.2
0.4
0.40.4
0.4
0.6
0.6
0.6
0.6
0.6
0.8
0
0.2
0.6
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(d) g∗/gT (c = cT )
000
0
0.2
0.2
0.2
0.4
0.4
0.4
0.6
0.6
0.6
0.8
0.8
0.8
1
1
1
1
1
0.2
0.4
0.60.8
1
11.
2
1.2
1.2
1.4
1.6
1.82
0.41.4
1.6
1.2
1.4
1.8
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(e) w∗/wT (c = 1.1 cT )
0
00.
2
0.2
0.2
0.2
0.4
0.4
0.4
0.4
0.4
0.6
0.6
0.6
0.6
0.6
0.6
0.8
0.2
0
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(f) g∗/gT (c = 1.1 cT )
Figura 13: Influência das Condições Iniciais para b = 0.1 bT .
36
00
0
0.2
0.2
0.2
0.4
0.4
0.4
0.6
0.6
0.6
0.8
0.8
0.8
0.8
0.8
0.80.8
1
1
1
1
1
1.21.41.6
1.2
1.2
1.822
1.4
1.4
1.6
1.6
1.8
1.8
2
2
0.8
0.8
12
11
1
1
1
1
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(a) w∗/wT (c = 0.9 cT )
0
0.2
0.2
0.4
0.4
0.40.4
0.4
0.4
0.4
0.4
0.6
0.6
0.6
0.6
0.60.6
0.6
0.6
0.6
0.6
0.6
0.8
0.8
0.8 0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.2
0.2
0
0
0.2 0.20
0
1
1
1
1
1
1
0.80.8
11
11
1
1
0.6
0.6
1.2
0.4
1.2
0.2
0.4
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(b) g∗/gT (c = 0.9 cT )
0
0
0
0.2
0.2
0.2
0.4
0.4
0.4
0.6
0.6
0.60.8
0.8
0.8
0.8
1
1
1
1
1
1
1
1
1.21.2
1.2
1.2
1.41.4
1.4
1.4
1.61.6
1.6
1.6
1.81.8
1.8
1.8 2
2
22
1.2
1.2
1.4
1.4
1.6
1.6
1.8
1.8 2
2
21
10
0.20.40.60.81
1.21.41.61.82
1
0.8
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(c) w∗/wT (c = cT )
00
00
0.2
0.2
0.2
0.2
0.2
0.4
0.4
0.4
0.4
0.4
0.4
0.6
0.6
0.60.60.6
0.6
0.8
0.8
0.8
0.8
0
00
0.2
0.2
0.2
0.6 0.60.4 0.40.2 0.20 0
0.40.
40.
4
0.8
0.8
11
0.6
0.6
0.8
0.80.6
1
0.4
1
1.20.8
0.6
1.4
0.6
1.61
1.2 1
1
2
1
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(d) g∗/gT (c = cT )
00
00.
20.
20.
20.
40.
40.
60.
6
0.80.8
11
1
1
0.8
0.8
0.8
0.8
0.8
1
1
1
1
1
1.21.2
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(e) w∗/wT (c = 1.1 cT )
0.2
0.2
0.2
0.4
0.4
0.4
0.4
0.4
0.4
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.6
0.8
0.8
0.8
0.8
0.8
0.6
0.8
0.80.8
0.8
0.811.2
0.4
1.4
1.6
1.820.6
0.6
11
0.8
0.4
0.4
1
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(f) g∗/gT (c = 1.1 cT )
Figura 14: Influência das Condições Iniciais para b = bT .
37
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(a) w∗/wT (c = 0.9 cT )
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(b) g∗/gT (c = 0.9 cT )
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(c) w∗/wT (c = cT )
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(d) g∗/gT (c = cT )
w∗/wT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(e) w∗/wT (c = 1.1 cT )
g∗/gT
g 0/g T
w0/wT
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(f) g∗/gT (c = 1.1 cT )
Figura 15: Influência das Condições Iniciais para b = 10 bT .
38
5.3 Influência do Ruído de Medição
Com o intuito de avaliar localmente a variabilidade do código, foi escolhida apartir das figuras 13, 14 e 15 uma condição inicial comparável para os três casos deviscosidade, a saber: c = c T , com w0 = 0.9 w T e g0 = 1.2 g T .
São geradas N = 300 realizações de ruído gaussiano ε i, as quais se somam a yT
para produzir N sequências de medição {y i}, i = 1, 2, ..., N , com SNR ≈ 20. Paracada y i, tem-se uma rodada do algoritmo EKF (4.21)-(4.30), da qual se obtem umaestimativa θ∗i = [w∗i , g∗i ]T dos parâmetros. A partir desses θ∗i parâmetros é geradauma distribuição de probabilidade π(θ∗) válida para cada caso, isto é, são descartadosvalores w∗i ≤ 0 e g∗i ≤ 0. Em seguida, é realizada uma análise de propagação deincertezas similar à realizada na Seção 3.3, considerando-se o conjunto de Nmc = 1000amostras dos parâmetros considerando-se θ∗j ∼ π(θ∗).
As PDF encontradas para os três casos de viscosidade são apresentadas nafig.16.
0 0.5 1 1.5 2 2.5 3 3.5 40
5
10
15
20
25
30
w∗/wT
b = 0.1 bTb = bTb = 10 bT
(a) π(w∗/wT )
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50
5
10
15
20
25
30
w∗/wT
b = 0.1 bTb = bTb = 10 bT
(b) π(w∗/wT ) (zoom)
0 0.5 1 1.5 2 2.5 3 3.5 40
1
2
3
4
5
6
7
8
9
10
g∗/gT
b = 0.1 bTb = bTb = 10 bT
(c) π(g∗/gT )
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.50
1
2
3
4
5
6
7
8
9
10
g∗/gT
b = 0.1 bTb = bTb = 10 bT
(d) π(g∗/gT ) (zoom)
Figura 16: Função Densidade de Probabilidade dos parâmetros estimados ao se considerar diversascondições de ruído. À esquerda, as PDFs. À direita, um zoom nas PDFs de menor suporte.
39
(a) r(t) (b = 0.1 bT ) (b) β(t) (b = 0.1 bT )
(c) r(t) (b = bT ) (d) β(t) (b = bT )
(e) r(t) (b = 10 bT ) (f) β(t) (b = 10 bT )
Figura 17: Comparação entre a resposta true do sistema e a resposta com os parâmetros estimadosao se considerar diversas condições de ruído. Em cada figura: à esquerda, em cinza, deslocamentoprescrito de v = 10/s; à direita, em preto, v = 1/s.
40
É possível notar claramente pela fig.16a que o ponto escolhido(θ0 = [0.9 w T , 1.2 g T ]T ) é bastante estável no que se diz respeito à estimação de w,enquanto que para g este ponto se encontra numa fronteira das curvas de nível eportanto é esperado que mude bastante com o ruído. De modo geral, os resultadospara os casos de baixa e média viscosidade foram satisfatórios dada a condiçãoinicial considerada, mesmo que sobrestimando a viscosidade do real do sistema(subestimando g).
Na fig.17 são comparadas a resposta true do sistema com a resposta obtida comos parâmetros estimados. É considerada uma condição de deslocamento prescrito,como feito na Seção 2.2. Conforme indicado na fig.17a, as curvas claras à esquerdasão para um deslocamento prescrito a v = 10/s, enquanto que as curvas escuras àdireita são para v = 1/s. Comparando-se as figuras 16 e 17, nota-se que emborao suporte das PDFs para os casos de baixa e média viscosidade sejam similares, oenvelope de resposta para o caso de média viscosidade é mais espesso, como mostramas figuras 17c e 17d. Isso é coerente com a análise de sensibilidade realizada na Seção3.3, ver fig. 10d em comparação a fig. 10b.
Entretanto, mesmo para esta condição mais “favorável”, estimação no caso dealta viscosidade se mostrou insatisfatória, como observado pelo amplo suporte daPDF. Fica claro que a ampla incerteza nos valores estimados de g, e inclusive w,resulta no envelope de resposta apresentado nas figuras 17e e 17f.
41
6 Conclusão
Este trabalho apresentou um breve estudo do comportamento previsto para umacamada de adesivo pelo modelo simplificado de interface; e a influência que incertezasnos parâmetros constitutivos desde modelo causam em variáveis de interesse, atravésde uma análise global de sensibilidade.
É possível fazer uma relação entre o comportamento observado na análise desensibilidade e os resultados obtidos para as simulações de calibração. Isto é, emregime de baixa sensibilidade do parâmetro viscoso b, embora os valores estimadosapresentassem um viés inferior ao valor real utilizado para criação dos dados sintéticos,ainda é possível se observar uma resposta prevista satisfatória. Em regime de altasensibilidade (alta viscosidade), entretanto, observaram enormes dificuldades naestimação do parâmetro b.
Para trabalhos futuros é de interesse utilizar métodos mais apropriados para selidar com não-linearidades, com o intuito de melhorar os resultados da estimação doparâmetro viscoso em regime de alta sensibilidade. É também de interesse, dado umadesivo conhecido, estimar o expoente p na lei de evolução do dano que melhor seadequa à resposta do mesmo.
42
Referências
[1] SAUER, R.A. A survey of computational models for adhesion. Journal ofAdhesion, 92(2):81–120, 2016.
[2] RAOUS, M., CANGEMI, L., COCU, M. A consistent model coupling adhesion,friction, and unilateral contact. vol. 177(3-4):383–399, Julho 1999.
[3] KALMAN, R.E. A new approach to linear filtering and prediction problems.Transactions of the ASME–Journal of Basic Engineering, 82(Series D):35–45,1960.
[4] KAIPIO, J.P., SOMERSALO, E. Statistical and Computational InverseProblems. Ed. Springer, 2005.
[5] WANG J., QIN Q.H., KANG Y.L., et al. Viscoelastic adhesive interfacial modeland experimental characterization for interfacial parameters. Mechanics ofMaterials, 42(5):537–547, 2010.
[6] VALOROSO N., FEDELE R. Characterization of a cohesive-zone modeldescribing damage and de-cohesion at bonded interfaces. sensitivity analysis andmode-i parameter identification. International Journal of Solids and Structures,47(13):1666–1677, 2010.
[7] ALBUQUERQUE, E.B., CASTELLO, D.A., BORGES, L.A. Calibration ofadhesion models based on bayesian inference. A ser publicado em InverseProblems in Sciece and Engineering, 2016, ,
[8] CORIGLIANO A., MARIANI, S. Simulation of damage in composites bymeans of interface models: parameter identification. Composites Science andTechnology, 61(15):2299 – 2315, 2001.
[9] MARIANI, S., CORIGLIANO, A. Impact induced composite delamination:state and parameter identification via joint and dual extended kalman filters.vol. 194(50-52):5242–5272, Dezembro 2005.
[10] BOLZON G., FEDELE R., MAIER G. Parameter identification of a cohesivecrack model by kalman filter. Computer Methods in Applied Mechanics andEngineering, 191(25–26):2847–2871, 2002.
[11] WANG, X. LIANG, L. XU, Y., LI, X. Inverse parameter identification of cohesivezone model for simulating mixed-mode crack propagation. 51(13):2400–2410,Junho 2014.
43
[12] ODELSON, B.J., RAJAMANI, M.R., RAWLINGS, J.B. A new autocovarianceleast-squares method for estimating noise covariances. vol. 42(2):303–308,Fevereiro 2006.
[13] SOLONEN, A., HAKKARAINEN, J., ILIN , A. et al. Estimating model errorcovariance matrix parameters in extended kalman filtering. vol. 21:919–927,2014.
[14] ALMAGBILE, A., WANG, J., DING, W. Evaluating the performances ofadaptive kalman filter methods in gps/ins integration. vol. 9(1), 2010.
[15] MILLER, B.A., HOWARD, S.A. Identifying bearing rotor-dynamic coefficientsusing an extended kalman filter. vol. 52(5):671–679, Setembro 2009.
[16] DING, Y., ZHAO, B.Y., WU, B. Structural system identification with extendedkalman filter and orthogonal decomposition of excitation. vol. 2014, Fevereiro2014.
[17] FOUN, K. Identification of civil structural parameters using extended kalmanfilter. Master’s thesis, Massachusetts Institute of Technology, Cambridge,Massachusetts, EUA, 2010.
[18] ALBUQUERQUE, E.B. Abordagens Bayesianas na Calibração de Modelos deAdesão. PhD thesis, COPPE/UFRJ, Rio de Janeiro, RJ. Brasil, 2013.
[19] THRUN, S., BURGARD, W., FOX, D. Probabilistic Robotics. The MIT Press,2005.
[20] GROVES, P.D. Principles of GNSS, Inertial, and Multisensor IntegratedNavigation Systems, 2nd Edition. Artech House, 2008.
[21] REIF, K., GUNTHER, S., YAZ, E., et al. Stochastic stability of thediscrete-time extended kalman filter. IEEE Transactions on Automatic Control,vol. 44(4):714–728, Abril 1999.
[22] BIÇER, C., KOKSAL, E., OSBEK, L. Stability of the adaptive fading extendedkalman filter with the matrix forgetting factor. vol. 20(5):303–308, Janeiro 2012.
[23] XIONG, K., ZHANG, H., LIU, L. Adaptive robust extended kalman filter,kalman filter. In MORENO, V., PIGAZO, A., editor, Kalman Filter: RecentAdvances and Applications. InTech, 2009.
[24] WANG, J. Stochastic modeling for real-time kinematic gps/glonass positioning.vol. 46:297–305, 1999.
44
