ESTIMAÇÃO DE PARÂMETROS DE UM MODELO CONSTITUTIVO PARA
MATERIAIS VISCOELÁSTICOS
Hudson Viegas Alves Fernandes de Souza
Dissertação
de
mestrado
apresentada
ao
Programa de Pós-graduação em Engenharia
Mecânica, COPPE, da Universidade Federal do
Rio de Janeiro, como parte dos requisitos
necessários à obtenção do título de Mestre em
Engenharia Mecânica.
Orientador: Daniel Alves Castello
Rio de Janeiro
Outubro de 2011
ESTIMAÇÃO DE PARÂMETROS DE UM MODELO CONSTITUTIVO PARA
MATERIAIS VISCOELÁSTICOS
Hudson Viegas Alves Fernandes de Souza
DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO
LUIS COIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA
(COPPE) DA UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE
DOS REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE
EM CIÊNCIAS EM ENGENHARIA MECÂNICA.
Examinada por:
_____________________________________________
Prof. Daniel Alves Castello, D.Sc.
_____________________________________________
Prof. Fernando Alves Rochinha, D.Sc.
_____________________________________________
Prof. Ricardo Leiderman, D.sc.
RIO DE JANEIRO, RJ - BRASIL
OUTUBRO DE 2011
Souza, Hudson Viegas Alves Fernandes de
Estimação de Parâmetros de um Modelo Constitutivo
para Materiais Viscoelásticos / Hudson Viegas Alves
Fernandes de Souza. - Rio de Janeiro: UFRJ/COPPE,
2011.
VII, 82 p.: il.; 29,7 cm
Orientador: Daniel Alves Castello
Dissertação (mestrado) - UFRJ/COPPE, Programa de
Engenharia Mecânica, 2011.
Referências Bibliográficas: p. 78-82
1. Viscoelasticidade. 2. Estimação de Parâmetros. 3.
Inferência Bayesiana. I. Castello, Daniel Alves.
II.
Universidade Federal do Rio de Janeiro, COPPE,
Programa de engenharia Mecânica. III. Título.
iii
AGRADECIMENTOS
Em primeiro lugar agradeço a Deus por tudo que tem feito na minha vida. Ao
meu pai Ary, minha mãe Gete e aos meus irmãos Igor e Natalle, por sempre estarem
comigo me dando todo o tipo de apoio possível, pois sem eles eu não chegaria até aqui.
Ao meu professor e orientador Daniel Alves Castello, por toda a ajuda ao longo
deste período. Após esse tempo sendo orientado por ele, passei a admirá-lo, não
somente como um excelente profissional, mas também pela pessoa que é. Ao meu
“irmão” Vitor, que teve uma grande contribuição para esse projeto.
Ao Pesquisador Carlos Frederico Trotta Matt, do Centro de Pesquisas de Energia
Elétrica (CEPEL), pelas conversas construtivas que tivemos durante congressos e
simpósios, nos quais participamos. Aos professores Ney Roitman e Carlos Magluta, do
programa de engenharia civil da COPPE e a professora Lavínia Borges, por sempre me
atenderem forma tão atenciosa. A Flávia Cavalcante pelas diversas contribuições e
ajudas durante a realização deste trabalho.
Aos meus amigos de laboratório de Acústica e Vibrações: Jefferson, Wallace,
Vanderson, Manoela, Guilherme, Bianca, Aline, Sergio e Flávio. Aos amigos de
graduação: Cláudio, Eduardo Veloso, Christovam e Maycon, por estarem comigo desde
o início dessa jornada. Aos amigos de infância Vinícius Alves Portela e Daniel Neves
pelos momentos de descontração e pelas palavras de incentivo nos momentos difíceis.
Ao engenheiro Anderson (LAVI), por toda a ajuda. Aos professores do
Departamento de Engenharia Mecânica, por contribuírem para a minha formação, em
especial a Hélcio Orlande, Fernando Castro Pinto, Marcelo Amorim Savi, Fernando
Alves Rochinha e Nestor A. Z. Pereira, por transmitirem conhecimentos aplicados nesse
trabalho.
iv
Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos
necessários para a obtenção do grau de Mestre em Ciências (M.Sc.)
ESTIMAÇÃO DE PARÂMETROS DE UM MODELO CONSTITUTIVO PARA
MATERIAIS VISCOELÁSTICOS
Hudson Viegas Alves Fernandes de Souza
Outubro/2011
Orientadores: Daniel Alves Castello
Programa: Engenharia Mecânica
O presente trabalho tem como objetivo a estimação de parâmetros constitutivos
de um modelo para materiais viscoelásticos. O modelo constitutivo escolhido para
representar a dinâmica destes materiais é baseado no conceito de variáveis internas e o
problema inverso é formulado a partir de pseudo-experimentos tanto no domínio do
tempo como no domínio da frequência. Um procedimento de inversão estatística,
baseado na técnica de inferência bayesiana, é utilizado para estimar os parâmetros na
forma de distribuições de probabilidade marginal a posteriori. O método de Monte
Carlo via cadeia de Markov (MCMC) é aplicado através do algoritmo Metropólis
Hastings visando extrair amostras oriundas destas distribuições. Um sistema massamola-amortecedor de um grau de liberdade anexado a um componente viscoelástico é
estudado, e são investigadas dificuldades relacionadas tanto com o processo de
estimação de parâmetros, assim como com o cálculo de quantidades a posteriori.
Resultados obtidos a partir de um modelo de estrutura do tipo viga sanduíche também
são apresentados, e alguns conceitos básicos de projeto ótimo de experimentos são
utilizados objetivando-se maximizar o nível de informação sobre os parâmetros
estimados.
v
Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the
requirements for the degree of Master of Science (M.Sc.)
PARAMETER ESTIMATION OF A CONSTITUTIVE MODEL FOR
VISCOELASTIC MATERIALS
Hudson Viegas Alves Fernandes de Souza
October/2011
Advisor: Daniel Alves Castello
Department: Mechanical Engineering
This work is aimed at estimating constitutive parameters of a model for
viscoelastic materials. The constitutive model that was chosen to represent the dynamics
of these materials is based on the concept of internal variables and the associated
inverse problem is built on data from either time or frequency domain. A statistical
inversion procedure, based on bayesian inference technique, is used to estimate the
parameters as marginal posterior probability distributions. Markov chain Monte Carlo
sampling methods (MCMC), with the Metropolis-Hastings algorithm, are employed in
order to draw samples from the posterior distributions. A mass-spring-damper system
with one degree of freedom, attached to a viscoelastic core, is studied, and problems
related to the parameter estimation process as well as the computation of posterior
quantities are investigated. A sandwich structure is also used for the analyses, whose
model equations are solved by the finite element method. Finally, optimum experiment
design procedures are applied to this system in order to maximize the information level
about the estimated parameters
vi
Sumário
1. Introdução -------------------------------------------------------------------------------
1
2. Viscoelasticidade ------------------------------------------------------------------------
5
2.1 - Fluência ------------------------------------------------------------------------
5
2.2 - Relaxação ----------------------------------------------------------------------
7
2.4 - Módulo Complexo ------------------------------------------------------------
7
2.6 - Modelo de Variáveis Internas -----------------------------------------------
8
3. Calibração de Modelos -----------------------------------------------------------------
10
3.1 - Problema Inverso -------------------------------------------------------------
10
3.2 - Fontes de Incerteza -----------------------------------------------------------
11
3.3 - Inferência Bayesiana ---------------------------------------------------------
12
3.4 - Problema Inverso Determinístico ------------------------------------------
17
3.4.1 - Enxame de Partículas ----------------------------------------------
19
3.5 - Projeto Ótimo de Experimento ---------------------------------------------
20
3.6 - Propagação de Incertezas via Método de Monte Carlo -----------------
22
3.7 - Validação de Modelos --------------------------------------------------------
23
4. Resultados Preliminares -----------------------------------------------------------------
26
4.1 - Sistema Investigado -----------------------------------------------------------
26
4.2 – Definição do Problema Inverso ---------------------------------------------
30
4.3 – Resultados e Discussões -----------------------------------------------------
32
4.3.1 – Influência do Nível de Ruído e do --------------------------------
32
Número de Parâmetros
4.3.2 – Propagação de Incertezas via Método de Monte Carlo --------
42
4.3.3 – Cadeias de Markov e Análises de Convergência ----------------
48
4.3.4 – Estratégias Preliminares de Validação ----------------------------
55
4.3.5 – Modelagem Incorreta dos Erros Experimentais -----------------
57
5 – Viga Sanduíche -------------------------------------------------------------------------
63
5.1 – Estrutura Investigada ---------------------------------------------------------
63
5.2 – Experimentos Simulados ----------------------------------------------------
65
5.3 – Projeto Ótimo de Experimento Aplicado a Viga Sanduíche -----------
66
5.4 – Aplicação da Inferência Bayesiana ----------------------------------------
73
6 – Conclusões -----------------------------------------------------------------------------Referências Bibliográficas ---------------------------------------------------------vii
77
78
1 - Introdução
Em aplicações de engenharia, muitas estruturas mecânicas operam sob
influência de carregamentos dinâmicos, os quais, em determinadas condições podem
causar níveis excessivos de vibração podendo, em último caso, conduzir a falhas
estruturais. A fim de reduzir estes níveis de vibração, é uma prática comum anexar a
estas estruturas componentes feitos de materiais dissipativos tais como os materiais
viscoelásticos. Estes materiais apresentam uma alta capacidade de dissipação de
energia, e por este motivo são amplamente utilizados como atenuadores de vibração nas
indústrias civil, automobilística, aeronáutica e petrolífera [1-3].
Neste contexto, destaca-se o fato de que a capacidade de projetar sistemas que
permitam explorar ao máximo as características dissipativas destes materiais é
condicionada ao nível de compreensão que cientistas e engenheiros possuem sobre a
resposta mecânica dos mesmos. Adicionalmente, a exigência pela aceleração das etapas
envolvidas na especificação de projetos é uma tendência devido, principalmente, a
motivações de origem econômica. Devido a estes fatores, simulações computacionais
vem sendo cada vez mais empregadas dentro de engenharia, substituindo, em muitos
casos, o papel desempenhado pelos protótipos experimentais [4].
A partir das simulações computacionais, sistemas cada vez mais sofisticados e
realistas vêm sendo abordados em virtude da disponibilidade de computadores mais
potentes. Isso aumentou as expectativas de resolução de problemas complexos, mas
gerou uma demanda para o tratamento e compreensão de erros e incertezas inerentes a
modelos computacionais. Esta demanda está diretamente relacionada à busca pela
determinação do nível de confiabilidade das predições do modelo escolhido, a fim de
determinar se os resultados previstos nas simulações de fato podem ser levados em
conta durante o processo de tomada de decisão [5].
Dentro deste contexto surge o conceito de validação de modelos computacionais,
o qual estabelece procedimentos a serem seguidos a fim de determinar o grau com o
qual um modelo é uma representação acurada do “mundo real” a partir de uma
perspectiva de aplicação. Esta avaliação deve ser feita a partir de uma comparação
criteriosa entre experimentos de validação e predições do modelo, e sendo o resultado
1
de tais comparações considerado satisfatório, de acordo com um conjunto de métricas, o
modelo é declarado validado dentro de um domínio de aplicação [6].
No que diz respeito a construção de modelos computacionais, é importante
mencionar que estes devem ser desenvolvidos considerando as seguintes etapas:
modelagem conceitual, modelagem mecânica, projeto de experimento, estimação de
parâmetros e processos de validação. O processo de estimação de parâmetros consiste
em ajustar os parâmetros do modelo computacional a fim de melhorar a concordância
entre os dados experimentais e as predições computacionais [6]. Este procedimento
também é denominado na literatura como problema inverso e, em geral pode ser
abordado fazendo-se uso de abordagens determinísticas [7,8] ou estatísticas [9]. A
abordagem determinística é formulada a partir de um problema otimização, onde o
mínimo de uma função objetivo, a qual leva em conta uma dada métrica para avaliar a
discrepância entre os dados experimentais e simulações do modelo, é procurado. Devese destacar que esta abordagem vem sendo amplamente empregada na caracterização de
materiais viscoelásticos [10-16].
Retomando a discussão sobre processos de validação, é importante mencionar
que as comparações entre simulações do modelo e experimentos de validação devem
levar em conta as incertezas envolvidas [4]. Diversas fontes de incertezas estão
presentes na modelagem computacional de sistemas físicos, as quais podem ser
representadas através de distribuições de probabilidade ou, de acordo com as
argumentações de um trabalho recente, fazendo-se uso de intervalos de confiança sem
função densidade de probabilidade associada [17].
Portanto uma comparação justa entre predições do modelo e experimentos de
validação somente é possível a partir de uma abordagem probabilística [6]. Neste
sentido a solução do processo de estimação de parâmetros deve levar em conta as
incertezas que surgem tanto dos processos de medição assim como da modelagem
imperfeita do sistema de interesse e, sendo possível, esta solução deve fornecer
informação acerca das funções de densidade probabilidade dos parâmetros para que,
posteriormente, técnicas de quantificação e propagação de incertezas possam ser
aplicadas a fim de determinar o nível de confiança associado às predições obtidas [6].
Deste modo, sendo possível, o processo de validação deve se basear na comparação de
regiões de probabilidade de ocorrência das quantidades físicas de interesse e não na
comparação simplista de eventos únicos das predições do modelo e medições
experimentais, os quais são oriundos de processos estocásticos [4].
2
Para tornar possível a obtenção dos cenários contendo as predições possíveis,
dado um certo regime de operações, a etapa de estimação de parâmetros deve ser
realizada a partir de uma teoria de inversão estatística. Ressalta-se aqui o fato de que em
geral, as abordagens determinísticas apresentam como solução um único ponto no
espaço dos parâmetros e, desta forma, em geral não fornecem informação adequada para
ser incorporada aos processos de quantificação de incertezas. Em problemas inversos
estocásticos todas as quantidades envolvidas são modeladas como variáveis aleatórias e
suas respectivas incertezas são quantificadas a partir de distribuições de probabilidade.
A teoria de inversão estatística é baseada na inferência bayesiana, onde informações a
priori são combinadas com medições experimentais a fim de extrair o máximo de
informação sobre os parâmetros de modelagem na forma de uma distribuição de
probabilidade a posteriori [9,18]. Adicionalmente, esta ferramenta é útil tanto para
abordar processos que envolvem competição de modelos [20], assim como para
incorporar incertezas de modelagem no resultado do problema inverso, tais como nas
abordagens hierárquicas [21].
Considerando a discussão acima sobre construção e validação de modelos
computacionais, se verifica que as medições experimentais desempenham um papel
crucial tanto no desenvolvimento assim como no refinamento de modelos para sistemas
físicos. Por exemplo, os dados podem ser utilizados para atualizar o estado de
conhecimento sobre os parâmetros de um modelo, ou mesmo, para processos de
competição de modelos. Particularmente, no que tange os experimentos de calibração,
ou seja, aqueles utilizados durante a etapa de estimação de parâmetros, vale destacar
que estes devem ser realizados de forma a obter modelos computacionais que forneçam
predições tão confiáveis quanto possível, a fim de facilitar o processo de tomada de
decisão de engenheiros e analistas responsáveis por projetos mecânicos baseados em
simulações computacionais.
Dentro deste contexto, surge o conceito de projeto ótimo experimento [22], que
consiste na especificação de condições de contorno, condições iniciais e localização dos
sensores, na forma de um protocolo experimental, a fim de maximizar o estado de
conhecimento sobre os parâmetros estimados [22-24].
Uma possível abordagem para a determinação do projeto ótimo de experimento
é baseada na matriz de informação de Fisher, a qual assume as hipóteses de modelo
linear em relação aos parâmetros [23]. É importante mencionar que esta análise é local
e, sendo assim, não leva em conta a variabilidade dos parâmetros estimados. Entretanto,
3
existem abordagens formuladas dentro de uma filosofia bayesiana onde são levadas em
consideração, tanto informação priori assim como a distribuição de probabilidade dos
parâmetros [25,26].
O objetivo do presente trabalho consiste em realizar a estimação de parâmetros
de um modelo constitutivo para materiais viscoelásticos. Portanto, a primeira etapa
consiste em descrever a resposta mecânica dos materiais viscoelásticos, além de
apresentar o modelo constitutivo de variáveis internas [11], o qual é utilizado para
representar a dinâmica destes materiais. Na segunda etapa do trabalho a inferência
bayesiana é aplicada a fim de investigar dificuldades inerentes ao problema inverso.
Estas dificuldades são discutidas a partir dos resultados para distribuições a posteriori
dos parâmetros além de outras quantidades físicas a posteriori. Processos preliminares
de validação são apresentados e o efeito da modelagem incorreta dos erros
experimentais é avaliado.
A terceira e última etapa do trabalho consiste na aplicação de estratégias
probabilísticas simplificadas que visam definir um projeto ótimo de experimento. A
métrica D-optimally é usada como critério de tomada de decisão e o sistema investigado
consiste de uma estrutura do tipo viga sanduíche. Estas análises são realizadas
simulando
algumas
dificuldades
comuns
encontradas
em
procedimentos
de
caracterização tais como, ausência de informação a priori sobre os parâmetros
constitutivos e conhecimento limitado sobre o comportamento estatístico dos erros
experimentais.
A apresentação deste trabalho é dividida em 5 capítulos além da introdução. O
segundo capítulo apresenta detalhes sobre o comportamento mecânico dos materiais
viscoelásticos, além do modelo constitutivo de variáveis internas. O terceiro capítulo
apresenta: o conceito de problema inverso, a técnica de inferência bayesiana, a
abordagem determinística para problemas inversos, a abordagem de projeto ótimo de
experimento baseada na matriz de informação de Fisher, a técnica de propagação de
incerteza via monte Carlo e conceitos preliminares sobre validação de modelos. O
quarto capítulo apresenta resultados preliminares oriundos da aplicação da técnica de
inferência bayesiana em um sistema massa-mola amortecedor de um grau de liberdade
anexado a um componente viscoelástico. O quinto capítulo aplica os conceitos de
projeto ótimo de experimento em uma estrutura do tipo viga sanduíche, cujas equações
do modelo são resolvidas através do método de elementos finitos. No sexto capítulo são
apresentadas as conclusões do trabalho.
4
2- Viscoelasticidade
Viscoelasticidade é um tipo de comportamento mecânico observado em uma
variedade de materiais, compreendendo desde tecidos humanos até polímeros
industrializados. O comportamento destes materiais engloba tanto características de
sólidos elásticos assim como de fluidos viscosos. Pode-se, também, defini-los como
materiais com memória, ou seja, o estado de tensões, em um dado instante, não depende
apenas do estado de deformações naquele instante de tempo, mas também de toda
história de deformação a qual o material esteve sujeito. Outra característica exibida por
estes materiais é o comportamento dissipativo, o qual pode ser verificado analisando-se
uma curva tensão por deformação, quando são sujeitos a um carregamento do tipo
cíclico [27-29]. A frequência de excitação, ou a taxa de carregamento, também possui
influencia sobre a resposta mecânica dos materiais viscoelásticos, podendo amplificar
ou reduzir as propriedades de amortecimento [1,3]. Existem testes padronizados, tais
como os testes de fluência, relaxação e módulo complexo, cujo objetivo é de verificar se
o comportamento de um material é ou não viscoelástico.
2.1- Fluência
A fluência pode ser compreendida como a deformação contínua do material,
quando o mesmo é submetido a uma tensão constante. A resposta de fluência é
fundamental para a caracterização da dependência dos materiais viscoelásticos em
relação ao tempo. Deste modo, define-se a função de fluência, em cada instante de
tempo, como sendo a razão entre a deformação e a tensão em um ponto específico do
material, quando o mesmo é submetido a uma tensão constante σ 0 , como é apresentado
a seguir.
Q F (t ) =
ε (t )
σ0
(2.1)
5
onde QF representa a função de fluência do material e ε (t ) a deformação. Nos
materiais viscoelásticos, a resposta de fluência é caracterizada por um salto instantâneo
de deformação no momento da aplicação da carga, o que configura uma característica
de um sólido elástico. Posteriormente o material se deforma continuamente
apresentando uma combinação de efeitos elásticos e viscosos [27]. A partir da figura
(2.1) observa-se, que nos materiais viscoelásticos, a resposta de fluência apresenta três
estágios: primário, secundário e terciário. No estágio primário a taxa de deformação é
relativamente grande, decrescendo ao longo do tempo até atingir um estado de
equilíbrio, a partir deste ponto se dá início o estágio secundário, no qual a taxa de
deformação se torna constante. Por último, ocorre o estágio terciário, onde a taxa de
deformação tende a aumentar com o decorrer do tempo, e onde ocorre a ruptura do
material por fluência [30].
Figura 2.1- Os três estágios da fluência [30].
6
2.2- Relaxação
A relaxação pode ser compreendida como a evolução temporal da tensão em um
ponto específico do material, quando o mesmo é submetido a uma deformação
constante ε 0 . Portanto, defini-se a função de relaxação, em cada instante de tempo,
como segue.
Q R (t ) =
σ (t )
ε0
(2.2)
onde QR representa a função de fluência do material e σ (t ) a tensão. No momento em
que é imposta a deformação, a resposta de relaxação de um material viscoelástico é
caracterizada por um salto instantâneo no valor da tensão, fato que configura uma
característica de um sólido elástico [27]. Posteriormente a tensão necessária para manter
a deformação constante cai gradualmente com o tempo, conforme pode ser observado
na figura (2.2).
Figura 2.2- Tensão de relaxação.
2.3- Módulo Complexo
É uma abordagem comum representar as equações constitutivas para materiais
viscoelásticos no domínio da frequência. Desta forma, o módulo complexo é definido
como na equação abaixo [27].
7
G (ω ) =
σ (ω )
(2.3)
ε (ω )
onde σ (ω ) e ε (ω ) representam, respectivamente, os valores de tensão e deformação,
no domínio da frequência. É importante mencionar a maior facilidade para medir o
módulo complexo quando comparado com os experimentos de fluência e relaxação de
um material viscoelástico específico. Por outro lado uma vantagem de se trabalhar no
domínio do tempo é que, em geral é possível, após a devida transformação,
trabalharmos também no domínio da frequência. Ressalta-se também o fato de que, ao
longo dos processos de calibração de modelos que, dependendo da abordagem
escolhida, existe a possibilidade do modelo calibrado, no domínio da freqüência ser não
causal.
2.4- Modelo de Variáveis Internas
Neste item são apresentadas as equações constitutivas do modelo de variáveis
internas, as quais satisfazem a desigualdade de Clausius-Duhem, ou seja, a sua
formulação satisfaz todas as restrições impostas pela termodinâmica de processos
irreverssíveis [11,31]. Deste modo, as equações (2.4) e (2.5) representam as equações
constitutivas para o modelo de variáveis internas.
I
•
σ = Eε + ∑ E r (ε − ξ r ) + η ε
(2.4)
r =1
•
ξ r = br (ε − ξ r ),
r = 1,..., I
(2.5)
onde ξ r corresponde a r-ésima variável interna, E e E r correspondem aos parâmetros
de elasticidade, η e br representam os parâmetros viscosos e I o número total de
variáveis internas do modelo. Deve ser enfatizado que a inclusão de um conjunto de
variáveis internas nas equações constitutivas é baseado na termodinâmica de processos
irreversíveis, e baseado nestes princípios pode ser provado que a inclusão destas
variáveis sempre dissipa energia. Portanto, como a tensão na equação (2.4) depende das
8
variáveis internas, as quais representam o fenômeno dissipativo no interior do material,
estas naturalmente contribuem para a dinâmica do sistema [11].
Com o intuito de obter uma representação analítica para a propriedade do
módulo complexo é adotada a hipótese de que todos os parâmetros que caracterizam o
material são constantes em todo o domínio do corpo. Portanto, efetuando a
transformada de Laplace nas equações (2.4) e (2.5) e admitindo-se condições iniciais
nulas, chega-se ao módulo complexo do material viscoelástico, como apresentado na
equação abaixo.
∧
I
I
E r br
r =1 br + s
G (s) = E + ∑ E r − ∑
r =1
(2.6)
É importante destacar que, na literatura, o modelo de variáveis internas vem
sendo utilizado com o objetivo de construir modelos computacionais para os materiais
viscoelásticos, de modo que são mencionados alguns trabalhos mais recentes nesta linha
de pesquisa. CASTELLO et al. [11], utilizaram este modelo constitutivo na
caracterização experimental de uma fita viscoelástica, a qual estava inserida no interior
de uma estrutura do tipo viga sanduíche. SOUZA et al. [15], utilizaram experimentos de
um sistema específico para gerar situações de operação de cisalhamento do material
viscoelástico a partir da equação (2.6), onde também foram calculadas regiões de
confiança aproximadas para os parâmetros constitutivos utilizando-se o método de
amostragem hit and run. BORGES et al. [32] utilizaram este modelo constitutivo para
construir modelos computacionais de estrutura do tipo riser em catenária onde os
resultados obtidos apresentaram um bom nível de concordância com os dados medidos.
BORGES et al. [33], aplicaram estratégias de validação a este modelo constitutivo em
condições de operação de viga sanduíche para condições de contorno variadas.
Portanto, neste capítulo foram apresentadas, de forma bastante sucinta, algumas
informações básicas referentes a viscoelasticidade, de modo que informações mais
detalhadas podem ser obtidas em algumas referências clássicas [27-29].
9
3 - Calibração de Modelos
Este capítulo apresenta a teoria envolvida nos problemas de estimação
abordados neste trabalho. As incertezas envolvidas na computação científica são
discutidas e apresentadas de acordo com as classificações mais recentes encontradas na
literatura [17]. A inferência bayesiana é apresentada, dentro de uma filosofia de
problemas inversos estocásticos, como uma ferramenta capaz de incorporar estas
incertezas no modelo computacional através da definição de funções de densidade
probabilidade para os parâmetros estimados. O problema de estimação determinístico é
apresentado baseando-se no estimador de máxima verossimilhança, onde o método de
otimização do enxame de partículas é proposto para encontrar o ponto ótimo. A
abordagem de Fisher é apresentada com o objetivo de aplicar estratégias que visam
determinar o projeto ótimo de experimento. Por último, são apresentados conceitos
preliminares tanto sobre a técnica de propagação de incertezas via Monte Carlo, assim
como sobre procedimentos de validação de modelos computacionais.
3.1 - Problema Inverso
Este item tem como foco a apresentação da formulação inversa relativa aos
sistemas que serão analisados no presente trabalho.
O problema inverso pode ser definido uma vez que o problema direto tenha sido
apresentado. Para o problema direto, o fenômeno físico e os parâmetros que
caracterizam as equações do movimento de um sistema são conhecidos, e se deseja
determinar a resposta deste sistema y ∈ R
Ny
quando este é submetido à determinada
excitação. Para o problema inverso, embora a resposta do sistema seja parcialmente
conhecida, alguns parâmetros podem ser desconhecidos. Portanto, a resposta medida do
sistema y obs ∈ R
Ny
é utilizada para determinar os parâmetros desconhecidos. Deste
ponto em diante, por questão de simplicidade, é definido um vetor θ o qual contém
informação acerca dos parâmetros desconhecidos do sistema como segue.
{
θ = θ 1 , θ 2 , K , θ Nθ
}
T
(3.1)
10
3.2 - Fontes de Incerteza
A maioria dos sistemas de interesse prático em engenharia são representados por
equações diferenciais parciais, cujas soluções são dadas por valores exatos das
quantidades de interesse. Entretanto diversas fontes de incertezas estão envolvidas na
modelagem destes sistemas e, por isto, as predições computacionais podem ser
consideradas e tratadas como não determinísticas. Visando inferir sobre a variabilidade
e o nível de confiança das predições do modelo deve-se solucionar, através de técnicas
de quantificação de incertezas, equações diferenciais parciais estocásticas [5]. Nestas
equações as incertezas envolvidas são propagadas através do modelo e a solução é dada
na forma de funções densidade de probabilidade [17].
Portanto o primeiro passo na aplicação desta abordagem consiste em identificar
e caracterizar todas as incertezas envolvidas na modelagem computacional. Neste
sentido é adotada uma classificação de incertezas de acordo com a sua essência [17],
dividindo as incertezas em dois grupos: aleatórias e epistêmicas. A incerteza aleatória é
definida como a variação inerente de uma quantidade que, dado um número suficiente
de amostras do processo estocástico, pode ser caracterizada a partir de uma função
densidade de probabilidade. A incerteza epistêmica resulta de informações incompletas
acerca do modelo. No que se refere às incertezas epistêmicas, os autores ROY e
OBERKAMPF [17], argumentam que as mesmas poderiam ser caracterizadas fazendose uso de intervalos de confiança sem distribuição de probabilidade associada. Deve-se
mencionar que estas incertezas podem surgir de diversas fontes, tais como as entradas
do modelo, as aproximações numéricas e dos erros devido à forma do modelo [17].
As entradas do modelo podem ser divididas em dois conjuntos: parâmetros de
modelagem computacional e dados que descrevem a interação entre sistema e
vizinhança. Os parâmetros de modelagem computacional constituem fonte de incerteza
porque são estimados através de dados experimentais os quais são incertos devido a
limitações inerentes dos processos de medição [17]. Os dados que descrevem a
interação entre sistema e vizinha, tais como as condições de contorno e a excitação
sobre o sistema, são incertos devido impossibilidade de medir seus valores de forma
exata. As entradas do modelo, usualmente, são modeladas como incertezas aleatórias, e
11
a abordagem bayesiana é utilizada para inferir sobre as suas distribuições de
probabilidade.
Os erros de aproximação numérica ocorrem devido a necessidade de se obter
soluções aproximadas para as equações do modelo, pois estas, para problemas de
engenharia, são normalmente dadas por sistemas de equações diferenciais parciais que
raramente apresentam soluções analíticas. Deve-se destacar que existem diversas
estratégias empregadas no sentido de caracterizar este tipo de incerteza, as quais
constituem procedimentos de verificação [4].
Por último, os erros devido à forma do modelo têm origem nas hipóteses
conceituais, abstrações e formulações matemáticas sobre as quais o modelo é
construído. A caracterização da incerteza devido a forma do modelo é normalmente
estimada utilizando procedimentos de validação [4,17].
3.3 - Inferência Bayesiana
Portanto, considerando a discussão apresentada no item (3.2), se conclui que o
resultado do problema inverso deve incorporar as incertezas envolvidas. Na teoria de
inversão estatística todas as quantidades envolvidas são modeladas como variáveis
aleatórias e suas respectivas incertezas são quantificadas a partir de distribuições de
probabilidade. Portanto, definindo θ, Y e E como variáveis aleatórias que caracterizam
o nível de conhecimento sobre os parâmetros do modelo, saída do sistema e erros
experimentais, respectivamente, tem-se a seguinte relação.
Y = f (θ, E)
(3.2)
onde f : R n × R k é um operador que depende do modelo matemático. A teoria de
inversão estatística é baseada na inferência bayesiana, onde informações a priori sobre
os parâmetros são combinadas com medições experimentais a fim de extrair o máximo
de informação sobre θ na forma de uma distribuição de probabilidade a posteriori
π (θ / y obs ) [9]. Portanto, assumindo que as variáveis aleatórias do problema são
contínuas podemos escrever uma extensão do teorema de bayes para densidades.
12
π (θ/y obs ) =
π (y obs /θ) × π pr (θ)
π (y obs )
(3.3)
onde y obs e θ representam, respectivamente, eventos das variáveis aleatórias Y e θ .
Na equação (3.3) o termo π pr (θ) é denominado de distribuição a priori, pois quantifica
na forma de uma distribuição de probabilidade, toda informação disponível sobre os
parâmetros do modelo antes que qualquer medição sobre o sistema tenha sido
incorporada no processo de inferência. O termo π ( y obs / θ) é conhecido por
verossimilhança e representa a probabilidade dos dados observados condicionada a um
vetor de entrada θ . Uma situação bastante explorada em problemas inversos é a
ocorrência de erros aditivos, a qual é caracterizada pela independência da variável
aleatória E em relação as variáveis aleatórias θ e Y [9]. Para erros aditivos a equação
(3.2) toma a seguinte forma.
Y = f (θ) + E
(3.4)
Considerando a equação acima e após algumas manipulações matemáticas é possível
provar que a verossimilhança pode ser avaliada como segue.
π (y obs / θ) = π E (y obs − f (θ))
(3.5)
onde π E representa a distribuição de probabilidade dos erros experimentais. O termo
π ( y obs ) corresponde à probabilidade de ocorrência dos dados observados e na prática
não pode ser diretamente avaliado, seu valor independe do vetor de parâmetros θ e,
portanto pode ser considerado como um fator de normalização na equação (3.3). Como
consequência o valor da distribuição marginal a posteriori π (θ / y obs ) é proporcional ao
produto entre as distribuições prior e verossimilhança como indicado na equação (3.6).
π (θ / y obs ) ∝ π ( y obs / θ) × π pr (θ)
(3.6)
Entretanto, uma vez calculada a distribuição a posteriori é possível calcular o fator de
normalização π ( y obs ) a partir da equação abaixo.
13
π (y obs ) = ∫ π (y obs / θ)π (θ)dθ
Sendo π ( ~
y / y obs )
(3.7)
a distribuição de probabilidade das predições do modelo
condicionada a observações experimentais, a sua determinação é de grande importância,
seja como parte de processos de validação ou com o objetivo gerar predições
probabilísticas sobre o sistema físico de interesse [18]. Seu calculo é dado através da
equação abaixo.
π (~y / y obs ) = ∫ π (~
y / θ)π (θ / y obs )dθ
(3.8)
Neste ponto é importante mencionar que o principal problema prático na
aplicação da abordagem bayesiana está relacionado com questões computacionais
referentes ao calculo de integrais que envolvem a distribuição π (θ / y obs ) , pois a
determinação de estimadores a posteriori envolve a solução de integrais tais como as
definidas na equação (3.9).
J = ∫ G (θ)π (θ / y obs )dθ
(3.9)
•
Se G (θ) = θ , J representa o valor esperado a posteriori dos parâmetros do modelo.
•
Se G(θ) = (θ − µ)(θ − µ) T , sendo µ o valor esperado a posteriori, J representa a
•
matriz de covariância dos parâmetros do modelo.
y / θ, y obs ) , sendo ~
Se G(θ) = π (~
y uma observação futura, J representa o valor
esperado a posteriori distribuição de probabilidade da predição do modelo.
Portanto, a integração desempenha papel essencial na análise bayesiana,
substituindo o papel desempenhado pela otimização na abordagem clássica de
problemas inversos. As integrais apresentadas na equação (3.9) são multidimensionais e
muitas vezes de difícil solução, principalmente quando o número de parâmetros
envolvidos é elevado. Entretanto nas últimas duas décadas as dificuldades de integração
têm sido bastante reduzidas devido ao desenvolvimento de novas abordagens baseadas
14
em simulação, que visam gerar amostras aleatórias e utilizar estas amostras para
resolver problemas de inferência através da equação (3.9). Métodos de Monte Carlo via
cadeia de Markov (MCMC) podem ser utilizados para a simulação da distribuição de
probabilidade alvo. Estes métodos são baseados na seleção de amostras dependentes que
admitem como distribuição estacionária a função densidade alvo. Em particular o tipo
de dependência considerado é uma cadeia de Markov, onde um novo estado da amostra
depende somente do estado atual da cadeia. Portanto os métodos MCMC permitem
calcular inferências a posteriori do tipo dado na equação (3.9) utilizando a seguinte
aproximação [18].
J≈
1
k − k0
k
∑ G (θ
k
)
(3.10)
k = k 0 +1
onde θ k são amostras dependentes geradas por uma cadeia de Markov cuja distribuição
estacionária é a função densidade procurada [18], k 0 representa o final do período
transiente, denominado como burn-in e k representa o número de amostras em que
ocorre a convergência da cadeia. Na aplicação dos métodos MCMC existe uma
preocupação em se utilizar técnicas estatísticas para determinar a iteração na qual a
convergência da cadeia possa ser considerada alcançada. Também se deve mencionar
que nos métodos MCMC a cadeia de Markov deve satisfazer a propriedade de
ergodicidade, entretanto esta propriedade somente é alcançada quando a cadeia é
irredutível e aperiódica. Irredutibilidade significa que qualquer estado da cadeia pode
ser alcançado a partir de outro estado em um número finito de movimentos. Esta
propriedade é exigida para permitir a cadeia esquecer o ponto inicial [19]. Portanto, um
resumo dos procedimentos realizados para calcular estimadores a posteriori utilizando
métodos MCMC é apresentado abaixo [9,18,34].
•
Passo 1: Iniciar a cadeia a partir de valores arbitrários de θ
•
Passo 2: Gerar k amostras a partir de uma cadeia de markov cujo kernel de
transição é dado por T (θ k +1 θ k ) e admite como distribuição estacionária a
densidade procurada.
15
•
Remover k 0 amostras a partir do período de burn-in da cadeia e calcular os
estimadores requeridos utilizando a equação (3.9).
Portanto o núcleo de qualquer abordagem MCMC é o kernel de transição
T (θ k +1 θ k ) que representa a probabilidade do movimento de θ k da cadeia para o
próximo estado θ k +1 , portanto definindo a forma na qual o novo estado da cadeia é
obtido a partir do estado atual. Todos os algoritmos de MCMC utilizam a mesma
estrutura base apresentada acima, diferindo somente na forma com que o kernel de
transição é definido. Neste trabalho o kernel de transição é definido, baseado na
condição de reversibilidade, a partir do algoritmo Metropolis hastings.
Uma cadeia de Markov satisfaz a condição de reversibilidade se existe uma
distribuição π (.) tal que.
T (θ k +1 θ k )π (θ k +1 ) = T (θ k θ k +1 )π (θ k )
(3.11)
Quando esta condição é satisfeita a cadeia é reversível e a distribuição π (.) é
estacionária. Entretanto é possível propor uma pequena alteração, utilizando um kernel
de transição definido por α . T ∗ (θ k +1 θ k ) , onde a probabilidade de que um movimento
seja aceito não é mais determinada por T ∗ (θ k +1 θ k ) mas unicamente por um fator
α (θ k +1 θ k ) que é calculado como segue.
T ∗ (θ k θ k +1 )π (θ k +1 )

, 1
∗
k +1
k
k
 T (θ θ )π (θ )

α (θ k +1 , θ k ) = min 
(3.11)
No presente trabalho é utilizado um kernel de transição do tipo passo aleatório, de
forma que as amostras da cadeia são geradas da seguinte forma.
θ k +1 = θ k + p.w
(3.12)
onde p representa o passo de busca, enquanto w é um vetor que contém números
aleatórios com distribuição uniforme, cujos limites inferior e superior são dados,
16
respectivamente, por -1 e 1. Portanto, considerando o tipo de kernel empregado, tem-se
que a equação (3.11) toma a seguinte forma.
π (θ k +1 )

, 1
k
 π (θ )

α (θ k +1 , θ k ) = min 
(3.13)
3.4 - Problema Inverso Determinístico
Este item tem como foco a apresentação da formulação do problema inverso
determinístico, o qual será empregado, posteriormente, para auxiliar a aplicação tanto da
inferência bayesiana, assim como do projeto ótimo de experimento.
Considerando que a formulação geral do problema inverso já foi apresentada,
pode-se formular o problema inverso determinístico como aquele que fornece como
solução um vetor de parâmetros θ que maximiza a função de verossimilhança, que para
condição de erros aditivos e gaussianos toma a seguinte forma.
π (y obs / θ) =
1
2π ∑ e
1/ 2
 1

T
−1
exp− [y obs − y (θ)] ∑ e [y obs − y (θ)]
 2

(3.14)
onde Σ e representa a matriz de covariância dos erros experimentais. Tomando a
derivada da equação (3.14) se verifica que o ponto de máximo da função de
verossimilhança corresponde ao ponto de mínimo de uma função objetivo definida
abaixo.
S (θ) = [y obs − y(θ)] Σ e−1 [y obs − y(θ)]
T
(3.15)
Assumindo que os erros experimentais são não correlacionados e possuem variância
constante, a função objetivo pode ser simplificada como segue.
S (θ) = [y obs − y(θ)] [y obs − y(θ)]
T
(3.16)
17
Portanto, uma vez definida a função objetivo, deve-se realizar o procedimento de
minimização, o qual consiste na busca de um vetor de parâmetros θ para os quais a
função objetivo alcança o mínimo global. A princípio, a minimização da função
objetivo pode ser feita por qualquer método de otimização, porém fatos como a
frequente não linearidade dos modelos e a presença de mínimos locais tornam, muitas
vezes, a minimização da função objetivo uma tarefa difícil. Assim, a escolha adequada
do método pode ser determinante para o sucesso do procedimento de estimação de
parâmetros.
Basicamente, existem dois tipos de métodos de minimização: determinísticos e
heurísticos. Os métodos determinísticos são baseados no uso de técnicas clássicas do
cálculo diferencial para buscar os pontos onde o gradiente da função objetivo é nulo.
Entretanto para sistemas não lineares a convergência para o mínimo global depende de
uma estimativa inicial do vetor θ que seja próxima do mínimo global [35].
Por outro lado, os métodos heurísticos tendem a imitar processos de otimização
observados na natureza para minimizar a função objetivo, de forma que seus
procedimentos iterativos não fazem uso do calculo de gradiente [36]. Estes métodos
também são capazes de resolver problemas não lineares devido à habilidade em driblar
os mínimos locais.
Também deve-se destacar, que os métodos heurísticos apresentam um elevado
custo computacional quando comparado com os métodos determinísticos, e isto se deve
ao caráter aleatório da busca pelo mínimo global, que se reflete em um grande número
de avaliações da função objetivo. Dentre os diversos métodos heurísticos destacam-se:
Algortimo genético, Evolução diferencial, Enxame de partículas, Recozimento
Simulado entre outros.
No presente trabalho decidiu-se utilizar o método do Enxame de Partículas nas
análises onde a abordagem determinística de problemas inversos é empregada. Tal
escolha foi baseada na facilidade de implementação do algoritmo, além da capacidade
deste método em apresentar bons resultados em diferentes aplicações. [37].
18
3.4.1- Enxame de Partículas
O método do Enxame de Partículas consiste em um algoritmo iterativo de busca
do mínimo global o qual se baseia no comportamento gregário de animais (peixes,
pássaros, etc.). Este método foi proposto [38] de forma a representar a troca de
informações entre os elementos de um grupo. Neste método, o movimento de cada
partícula, em cada iteração, corresponde à soma de três termos distintos: o primeiro
termo é relativo à inércia das partículas e traduz o modo com que as partícula vem se
movendo ao longo do espaço de busca. O segundo termo está relacionado com a atração
da partícula pelo melhor ponto que já encontrou. O terceiro termo está relacionado com
a atração da partícula para o melhor ponto que todo o grupo já encontrou. Desta forma,
as seguintes equações foram propostas [38].
θ ik +1 = θ ik,d + v ik,+d1
(3.17)
v ik,+d1 = v ik,d + c1 .r1 .(p ik,d − θ ik,d ) + c2 .r2 .(p kg ,d − θ ik,d )
(3.18)
onde os índices k, i e d denotam, respectivamente, a iteração, a partícula e a direção de
busca, v é a velocidade e θ representa a posição das partículas no espaço de busca, c1 e
c2 são duas constantes positivas, chamadas respectivamente de parâmetro cognitivo e
social, r1 e r2 são vetores, gerados a cada iteração, formados por números aleatórios
com distribuição uniforme no intervalo [0, 1]; pi é o melhor ponto encontrado pela
partícula i, e pg responde pelo melhor valor encontrado por todo enxame.
O Método do Enxame de Partículas apresenta um comportamento característico
ao longo das iterações. Nas iterações iniciais, o caráter aleatório é alto e as partículas
realizam uma busca global sobre toda região de busca. Com a evolução das iterações,
dá-se início à etapa local da busca, onde as partículas se concentram em torno das
regiões mais promissoras encontradas durante a etapa global [35].
Com o objetivo de balancear o caráter global e local da busca pelo ponto ótimo,
foi proposta [39] uma alteração na equação do método, que consistiu na introdução de
um novo parâmetro, chamado de peso de inércia ou fator de inércia, o qual pondera o
termo relativo à velocidade prévia da partícula, de acordo com a seguinte equação:
19
v ik,+d1 = w.v ik,d + c1 .r1 .(p ik,d − θ ik,d ) + c2 .r2 .(p kg ,d − θ ik,d )
(3.19)
Este peso de inércia pode ser uma constante positiva ou mesmo uma função do número
de iterações, podendo ser positivo linear ou não linear. Foi verificado [35], que valores
altos de w aumentam o caráter global da busca, enquanto que valores pequenos de w
aumentam o caráter local da busca, fazendo com que as partículas convirjam
rapidamente, sem que haja uma boa exploração do espaço de busca. Desta forma, uma
alternativa interessante é iniciar a busca com um valor alto de w , e diminuir este valor
ao longo das iterações [35].
A configuração do enxame de partículas utilizada neste trabalho considera a
inserção do peso de inércia, o qual varia linearmente de 1 e 0 ao longo das iterações.
Vale mencionar que visando evitar a explosão das partículas durante as iterações foi
utilizada uma barreira, a qual reposiciona as partículas, no interior de uma região prédefinida, de forma aleatória, quando as mesmas ultrapassam uma determinada região
[40].
3.5 - Projeto Ótimo de Experimento
Uma vez que a técnica de estimação que será utilizada no presente trabalho foi
adequadamente apresentada, deseja-se focar no projeto ótimo de experimento. O que se
deseja fazer é conduzir um experimento no qual a resposta do sistema seja medida de tal
maneira que permita deduzir o máximo de informação possível sobre o vetor de
parâmetros. Deve-se destacar que neste contexto, a definição do experimento envolve a
especificação de condições de contorno, condições iniciais, tempos de amostragem,
localização dos sensores e outros, na forma de um protocolo experimental [23].
Portanto, o objetivo geral do projeto ótimo de experimento é definir um protocolo
experimental que permita estimar o vetor de parâmetros com a máxima confiança
estatística
Deve-se destacar que a formulação aqui apresentada assume a hipótese de que
∧
em uma certa região ao redor do ponto de ótimo θ , o modelo possui uma aproximação
linear. Portanto, uma forma interessante de extrair informação acerca do problema de
20
otimização associado ao projeto ótimo de experimento pode ser obtida através de
∧
análises da função objetivo em torno do vetor de parâmetros estimados S (θ) [24].
Portanto, consideremos a expansão de segunda ordem da série de Taylor da função
objetivo.
∧
∧
S (θ) = S (θ) + (θ − θ)∇S ∧ +
θ
∧
∧
1
(θ − θ ) T H ∧ (θ − θ )
θ
2
(3.20)
onde H e ∇ representam, respectivamente a matriz hessiana e o vetor gradiente. Devese destacar que a matriz hessiana pode ser calculada, de forma aproximada, como na
equação abaixo.
 ∂y 
 ∂y 
H ∧ = 2  Σ e−1  
θ
 ∂θ 
 ∂θ 
T
(3.21)
Considerando também, a relação entre a matriz de covariância dos parâmetros
estimados e a matriz hessiana [22].
H ∧ = 2 Σ θ−1
(3.22)
θ
defini-se a matriz de informação de Fisher, como o inverso da matriz de covariância dos
parâmetros, a qual toma a forma apresentada na equação a seguir.
 ∂y 
 ∂y 
M =   Σ e−1  
 ∂θ 
 ∂θ 
T
(3.23)
Na abordagem de Fisher, a qual é utilizada no presente trabalho para definição
do projeto ótimo de experimento, o protocolo experimental é escolhido de forma a
minimizar os componentes da matriz de covariância dos parâmetros através da
maximização dos componentes da matriz de informação.
21
Entretanto, como M é uma matriz, o projeto ótimo de experimento é
determinado pela exigência de que alguma medida desta matriz, a qual usualmente é
dada por um escalar, seja maximizada. Dentre as diversas métricas propostas decidiu-se
trabalhar com a D-optmally, que é apresentada na equação abaixo.
Dopt = (det (M ))
(3.24)
A estratégia utilizada neste trabalho para determinar o ponto de máximo, é
baseada na enumeração de possíveis protocolos experimentais e posterior cálculo da
métrica D-optmally para cada configuração proposta, de modo que o projeto ótimo é
definido a partir de uma comparação dos valores obtidos.
É importante mencionar que, como apresentado na equação (3.23), o cálculo da
matriz de Fisher requer somente informação acerca do comportamento estatístico dos
erros, e não necessita de informação acerca dos dados medidos. Outra observação
importante diz respeito ao fato de que, a formulação aqui apresentada é baseada em
análises locais, as quais não levam em consideração a variabilidade dos parâmetros.
Entretanto no presente trabalho, é proposta e aplicada uma abordagem probabilística
simplificada para o projeto ótimo de experimento.
3.6 - Método de Monte Carlo
As diversas incertezas presentes na modelagem computacional impossibilitam
predições determinísticas, e sendo assim, sendo possível, devem ser utilizadas
estratégias que visem incorporar estas incertezas dentro do modelo computacional. Um
tipo de estratégia utilizada para abordar este problema consiste em modelar as incertezas
como variáveis aleatórias e reformular os sistemas de equações diferenciais dos
modelos determinísticos originais como sistemas de equações diferenciais estocásticas
[5]. Dentre as possíveis formas de solucionar estas equações, escolheu-se trabalhar com
um método baseado em amostragem, o qual é denominado método de Monte Carlo.
Neste método é gerado um número k de amostras independentes, as quais são oriundas
de funções densidade de probabilidade, as quais representam o nível de informação
22
acerca dos parâmetros de modelagem. A partir destas amostras, são geradas k predições
determinísticas, as quais podem ser usadas para o cálculo de estimadores de valor
esperado e de variância por exemplo. O método de Monte Carlo é de simples aplicação,
pois requer apenas soluções repetidas para as equações determinísticas do sistema. Por
outro lado, deve-se destacar que as soluções obtidas são aproximações da resposta
desejada, devido ao fato de se utilizar um número finito de simulações, e que o nível de
precisão desta aproximação depende do número de amostras consideradas. Portanto
devem ser empregados critérios de convergência a fim de avaliar o número mínimo de
amostras, para as quais pode ser considerado que a quantidade calculada convergiu [5].
3.7 -Validação de Modelos
Este item tem como foco a apresentação dos conceitos básicos da filosofia de
validação de modelos. Neste ponto é importante destacar que embora a discussão aqui
apresentada englobe diversos procedimentos que devem ser seguidos a fim de validar
um modelo, neste trabalho serão abordados somente procedimentos preliminares de
validação. Informações mais detalhadas acerca de estratégias de Verificação e
Validação (V&V) de modelos computacionais podem ser obtidas em diversas
referencias [4,6].
A fim de introduzir adequadamente o conceito de validação, deve-se apresentar
as etapas envolvidas na construção de um modelo computacional. Dentro deste
contexto, primeiro deve ser considerada a modelagem conceitual do problema, onde
uma série de descrições e hipóteses sobre os processos físicos envolvidos na resposta
mecânica do sistema são adotadas. Posteriormente tem-se a modelagem matemática,
onde são definidas equações, valores de contorno e condições iniciais que descrevem o
modelo conceitual. Deste modo, a modelagem física representa uma interpretação do
mundo real baseada em observações experimentais e, no melhor dos casos, corresponde
apenas a uma aproximação. E, por último, defini-se o modelo computacional o qual
responde pela etapa de implementação numérica do modelo matemático, usualmente na
forma de uma discretização numérica, com algoritmo de solução e critério de
convergência [4,6].
Uma vez apresentados estes conceitos, pode-se definir os processos de validação
como aqueles que visam determinar o grau com o qual um modelo é uma representação
23
acurada do “mundo real” a partir de uma perspectiva de utilização e, deste modo, visam
verificar a adequação dos modelos conceitual e matemático para a realidade de interesse
[6]. Por outro lado, um dos principais motivos pelo qual se constrói um modelo
computacional é para fornecer predições para aplicações onde dados experimentais não
estão disponíveis. Portanto, neste sentido a aplicação de procedimentos de validação é
fundamental, pois estes permitem inferir sobre o nível de confiança, de um modelo
proposto para representar o comportamento real de um sistema para situações diferentes
das quais o modelo foi calibrado [6].
Figura 3.1 – Validação, Calibração e predição [41].
A figura (3.1) apresenta aspectos importantes da validação, assim como de
questões relativas a predições e calibração de modelo. Deste ponto em diante, por
24
questão de simplicidade, define-se (RFS) como uma resposta física do sistema.
Retomando a discussão sobre a figura acima, tem-se que, num procedimento de
validação, a mesma RFS deve ser obtida, tanto a partir do modelo computacional assim
como do experimento. Deve-se destacar que algumas vezes a RFS experimental é obtida
a partir de medições indiretas as quais podem demandar um significante processamento.
Nestes casos, é importante que os dados sejam processados da mesma maneira, tanto no
modelo computacional, assim como no experimento [41]. É importante mencionar que
durante o planejamento e projeto do experimento, deve ocorrer uma freqüente
comunicação entre modeladores e experimentalistas, e ainda após a realização do
experimento, os experimentalistas devem fornecer informações aos modeladores sobre
todas as quantidades importantes de entrada necessárias para conduzir as simulações
computacionais [41].
Analisando a parte central da figura (3.1) se observa que as RFSs experimental e
computacional são entradas de um operador que retorna como saída uma métrica de
interesse ou, em algumas situações, um conjunto de métricas. O segundo passo na
validação lida com comparações entre as métricas resultantes e valores préestabelecidos a partir de uma determinada exigência de acerácea. Se estas comparações
forem consideradas satisfatórias conclui-se que não existem fortes evidências para
invalidar o modelo. Em caso contrário, algumas medidas devem ser tomadas tais como :
(i) revisar os protótipos experimentais, condições de operação e sistemas de mediação,
(ii) reiniciar o processo de calibração de modelos, (iii) possivelmente até mesmo
revisitar a estrutura do modelo adotado e etc.
A definição da exigência de acerácea deve levar em consideração diversos
fatores e, de acordo com ROY e OBERKAMPF [4], tais fatores englobam a
complexidade do modelo, o aumento da incerteza devido a extrapolação do modelo,
tolerância de risco nas decisões envolvidas e consequências da falha ou subdesempenho do sistema de interesse.
25
4 - Resultados Preliminares
A teoria envolvida na inferência bayesiana foi apresentada no capítulo anterior,
fornecendo o embasamento necessário para a aplicação desta ferramenta na solução de
problemas inversos de interesse e, deste modo, este ferramental será utilizado na
estimação dos parâmetros constitutivos do modelo de variáveis internas na forma de
uma distribuição de probabilidade marginal a posteriori. Entretanto, é importante
mencionar, que sendo possível, a inferência bayesiana deve ser utilizada de forma a
maximizar o estado de conhecimento, tanto dos parâmetros estimados assim como das
quantidades físicas calculadas a partir destes parâmetros através do modelo.
Neste sentido, este capítulo é de grande importância, pois apresenta análises
sobre a influência de diversos fatores, tais como o nível de poluição experimental,
modelagem dos erros experimentais e número de parâmetros estimados, sobre os
resultados, tanto do problema de inferência estatística, assim como de processos de
propagação de incerteza. Também são abordadas estratégias preliminares que visam
avaliar a capacidade de predição do modelo dentro de uma filosofia de validação de
modelos.
Por outro lado, métodos que visam reduzir o custo computacional das
simulações envolvidas não serão aplicados, principalmente devido a complexidade
envolvida nestas abordagens que constituem, por si só, um tema de pesquisa [5,42,43].
Por este motivo decidiu-se trabalhar, neste capítulo, com um problema computacional
simples, onde a resposta do sistema pode ser simulada milhares de vezes em períodos de
tempo inferiores a dez minutos, isto considerando que o desktop utilizado nas
simulações é um Core I7 de 3.7GHz e 16Gb de memória.
4.1 - Sistema Investigado
O problema investigado compreende um sistema massa-mola-amortecedor de
um grau de liberdade anexado a um componente viscoelástico, onde o comportamento
mecânico do componente dissipativo é representado pela representação discreta da
equação constitutiva apresentada nas equações (2.4) e (2.5).
26
Figura 4.1 – Sistema dinâmico investigado.
Os parâmetros do sistema estão ilustrados na figura (4.1) e seus valores apresentados na
tabela (4.1), onde M , k , c , representam, respectivamente, massa, rigidez e
amortecimento, enquanto que
E , E1 e b1 , são parâmetros do modelo de variáveis
internas, o qual apresenta apenas uma variável interna. As quantidades y (t ) e f(t)
representam, respectivamente, a saída do modelo dada pela solução do problema direto
e a excitação sobre o sistema.
M
1 [ Kg]
k
10 [N/m]
c
0.04 [ Nsm-1]
E
1 [N/m]
E1
6 [N/m]
b1
5 [Nsm2]
Tabela 4.1 – Parâmetros do Sistema investigado.
Portanto, o problema direto para o sistema acima compreende uma equação diferencial
ordinária com coeficientes constantes, que é apresentada abaixo.
••
•
M z (t ) + c z (t ) + k z (t ) = F (t ) − FR (t )
(4.1)
27
•
onde
••
representam, respectivamente, deslocamento, velocidade e
z (t ), z (t ), z (t )
aceleração do sistema, enquanto que FR (t ) responde pela força de restituição devido ao
componente dissipativo, a qual pode ser calculada como segue.
FR (t ) = E z (t ) + E1 ( z (t ) − ξ (t ) )
(4.2)
onde ξ (t ) representa a variável interna do modelo para viscoelasticidade. Observe que a
força de restituição depende tanto do deslocamento do sistema assim como da variável
interna. Entretanto a variável interna evolui no tempo de acordo equação (2.5), de modo
o problema direto passa a compreender um sistema acoplado de equações diferenciais, o
qual pode ser escrito de forma matricial a partir da definição das quantidades
•
z1 (t ) = z (t ) e z 2 (t ) = z (t ) como segue.
•

0
 z• 1 (t )  

 
 z 2 (t ) = − (k + E + E1 ) / M
•  
b1
 ξ (t )  
1
−c/M
0


0   z 1 (t )   0 


E1 / M   z 2 (t ) +  F (t ) / M 
− b1   ξ (t )   0 




(4.3)
Esta equação matricial é resolvida através da formulação de espaço de estados,
utilizando a plataforma computacional Matlab®. As figuras (4.2) e (4.3) apresentam,
respectivamente, a amplitude da função de resposta em frequência e a resposta do
sistema a um estímulo do tipo impulso, onde, em ambos os casos, a saída do problema
direto é dada pelo deslocamento. Nestas figuras são apresentados os resultados tanto
para o sistema com componente dissipativo, assim como sem componente dissipativo.
A partir destas representações do comportamento do sistema em domínios
diferentes se conclui que o nível de dissipação do sistema aumentou de forma
significativa com a inclusão do componente viscoelástico. Na figura (4.2) é possível
observar este aumento devido ao maior espalhamento de função de transferência na
região próxima a frequência natural. Na figura (4.3) se observa que o sistema com
componente viscoelástico para de oscilar em aproximadamente dez segundos, enquanto
que o outro sistema mantém altas amplitudes de vibração em tempos superiores a trinta
28
segundos, indicando uma alta dissipação de energia por parte do componente
dissipativo.
Figura 4.2 – Função de transferência, saída do modelo: Deslocamento.
Figura 4.3 – Resposta a um estímulo tipo impulso, saída do sistema: Deslocamento.
29
4.2 - Definição do problema inverso
Uma vez definido o modelo constitutivo para o núcleo viscoelástico e o
problema direto para o sistema investigado, visando aplicar a inferência bayesiana,
deve-se especificar o experimento de calibração, a função de verossimilhança e a
distribuição a priori para os problemas analisados.
O experimento de calibração deve fornecer tanta informação quanto possível
sobre o sistema, de forma a permitir o calculo de distribuições a posteriori mais
informativas, ou seja, com menores dispersões em torno dos valores médios. Dentre as
diversas quantidades que podem ser medidas como deformação, aceleração,
deslocamento, frequencias naturais e funções de resposta em frequencia, decidiu-se
trabalhar com a evolução temporal da velocidade do sistema, quando o mesmo é
submetido a condições iniciais nulas e a um determinado tipo de forçamento. Estes
dados experimentais são organizados dentro de um vetor de dados observados
•
•

y obs =  z (t1 ) L z (t N y ) , onde N y representa o número de pontos gravados em um


arquivo digital. Visando interpretações conclusivas a partir dos resultados do problema
de inversão estatística, escolheu-se trabalhar com experimentos simulados, onde
resultados oriundos do problema direto são corrompidos por um termo de perturbação
de acordo com a equação (4.4). Desta forma, com esta escolha de pseudo-experimento,
não estamos considerando erros de modelagem, pois a estrutura de modelo utilizada
para a formulação do problema inversa é a mesma utilizada ao longo do processo de
calibração e, sendo assim, não fosse pela perturbação aletória dos dados pseudoexperimentais, teríamos o tão comentado crime inverso na sua definição pura [9]
y obs (t ) = y (t ) + v(t )
(4.4)
Uma vantagem de se trabalhar com experimentos simulados, é o fato de esta abordagem
permite assegurar, no mínimo, uma solução para o problema inverso por construção
[20], além de possibilitar a comparação entre os resultados obtidos e os valores exatos
dos parâmetros, assim permitindo avaliar a eficiência do método de inversão adotado.
O vetor de dados observados é gravado em 100 pontos (N y = 100 ) e com
30
frequência de amostragem de 10hz. O sistema é submetido a um forçamento do tipo
varredura de senos como segue.
f (t ) = 100 cos(2πΩ(t ))
(4.5)
onde a frequência de excitação depende do tempo e cresce linearmente de
Ω min = 0 hertz até Ω max = 1 hertz . Observando a equação (4.4) se verifica que os erros
são aditivos. Portanto, supondo que os erros experimentais seguem uma distribuição
gaussiana, tem-se que a função de verossimilhança pode ser avaliada pela equação
(3.14).
Uma vez que a função de verossimilhança tenha sido adequadamente definida,
deve-se escolher a distribuição a priori π pr (θ ) que melhor representa o estado de
conhecimento sobre os parâmetros antes da execução do processo de estimação. Esta
definição representa um desafio do ponto de vista filosófico, pois se deve transformar
informação objetiva a respeito do sistema analisado em informação não objetiva através
de distribuições de probabilidade. A situação abordada neste capítulo é aquela em que o
nível de informação sobre os parâmetros é pequeno. Neste caso, uma possível estratégia
consiste em utilizar distribuições prior uniformes. Levando isto em consideração
decidiu-se utilizar como prior distribuições uniformes com parâmetros independentes,
de modo que a função densidade de probabilidade toma a seguinte forma.
Nθ
π pr (θ) = ∏
i =1
1
L − LiI
(4.6)
i
S
onde LiS e LiI representam, respectivamente, os limites superior e inferior da
distribuição a priori do parâmetro θ i , dado que o vetor de parâmetros é definido como
{
}
θ = θ 1 ,θ 2 ,θ 3 , onde θ 1 = E , θ 2 = E1 e θ 3 = b1 . Os limites inferior e superior são
dados por LI = {0,0,0} e LS = {10,10,10}. Vale destacar que os parâmetros M , k e c são
assumidos como conhecidos e que embora a distribuição prior definida acima se baseie
na estimação de três parâmetros constitutivos, em algumas análises um número menor
de parâmetros será estimado. Quando isto ocorrer, o parâmetro constitutivo não
estimado será assumido como conhecido.
31
O método de Monte Carlo via Cadeia de Markov (MCMC) é utilizado para
amostrar vetores de parâmetros oriundos da distribuição a posteriori de interesse. Para
tal foi aplicado o algoritimo Metropólis Hastings. O critério utilizado para determinar o
número mínimo de amostras para o qual pode ser considerado que a cadeia de Markov
atingiu o regime permanente, se baseou em análises gráficas das amostras cadeia. A
determinação da convergência da cadeia é obtida a partir de comparações gráficas das
funções densidade de probabilidade dos parâmetros para diferentes números de
amostras da cadeia de Markov. Neste ponto, é importante mencionar que existem
métodos mais pragmáticos, os quais levam em consideração métricas para determinar a
convergência da cadeia.
4.3 - Resultados e discussões
Este item visa abordar algumas questões comuns que surgem durante a aplicação
de problemas inversos, tais como: Como a precisão de um processo de medição pode
influenciar nos valores da distribuição a posteriori? Quais as complicações numéricas
de se estimar diversos parâmetros simultaneamente? Qual o efeito da modelagem
incorreta dos erros experimentais sobre os resultados do problema inverso?
Estes questionamentos serão abordados através da aplicação da inferência
bayesiana em conjunto com análises de sensibilidade, análises de correlação e utilização
de técnicas de propagação de incerteza via o método de Monte Carlo.
4.3.1 - Influência do nível de ruído e do
número de parâmetros estimados
Neste subitem são investigados o efeito do nível de poluição sobre os dados
experimentais, assim como do número de parâmetros estimados sobre os resultados do
problema inverso. Para tal, o vetor de dados observados y obs é gerado com erros
experimentais aditivos, não correlacionados, gaussianos com média zero e com
32
variância constante dada por σ 2 . Portanto a função de verossimilhança é dada pela
equação (3.14) e a matriz de covariância dos erros é dada por:
∑ e = σ 2 II
(4.7)
onde II representa uma matriz identidade N y × N y e o valor da variância dos erros
experimentais é calculada para uma razão sinal ruído (SNR) definida como segue.
σ 
SNR = 20 log 10  R 
 σ 
(4.8)
onde σ R representa o desvio-padrão do sinal de referência, o qual é dado pelo resultado
•
do problema direto y (θ) = z (θ, t ) . Repare que menores valores de SNR implicam em
maiores variâncias dos erros experimentais e vice-versa. A título de ilustração a figura
(4.4) apresenta tanto a excitação f (t ) sobre o sistema investigado, assim como o vetor
de dados observados y obs (t ) , para uma razão sinal ruído SNR=20db. São gerados dois
vetores de dados observados para serem utilizados no processo de estimação de
2
parâmetros, os quais são designados por y 1obs e y obs
. Em y 1obs o valor do desvio padrão é
2
o valor do desvio-padrão é
definido para uma SNR=10db, enquanto que em y obs
definido para SNR=20db.
Figura 4.4 – Força de excitação e velocidade observada (SNR=20 dB).
33
Dados observados
Vetor de parâmetros estimados
Caso A
2
y obs
( SNR = 20 dB )
θ A = {E}
Caso B
2
y obs
( SNR = 20 dB )
θ B = {E1 }
Caso C
2
y obs
( SNR = 20 dB )
θ C = {b1 }
Caso D
2
y obs
( SNR = 20 dB )
θ D = {E, E1 , b1 }
Caso E
y 1obs ( SNR = 10 dB )
θ E = {E}
Caso F
y 1obs ( SNR = 10 dB )
θ F = {E1 }
Caso G
y 1obs ( SNR = 10 dB )
θ G = {b1 }
Caso H
y 1obs ( SNR = 10 dB )
θ H = {E , E1 , b1 }
Tabela 4.2 – Características dos casos A até H.
Visando avaliar tanto o efeito do nível de ruído, assim como o efeito do número
de parâmetros estimados sobre os resultados do problema inverso, foram rodados oito
casos, os quais foram designados por A, B, C, D, E, F, G, H. As características de cada
caso são apresentadas na tabela (4.2) e os resultados obtidos para as distribuições a
posteriori, a partir da aplicação da inferência bayesiana nestes casos, são apresentados
na tabela (4.3).
E [N/m]
E1 [N/m]
b1 [Nsm2]
1
6
5
Valor exato
Caso A
µ
σ
0.9538
0.0452
Caso B
µ
σ
5.9316
0.0694
Caso C
Caso D
0.9226
0.3127
Caso E
1.0470
0.1237
Caso F
6.0213
0.3246
5.8728
0.2057
Caso G
Caso H
1.2987
0.5494
6.2969
0.9209
µ
σ
4.9922
0.0738
5.0107
0.7810
5.0824
0.0738
5.9657
1.8368
Tabela 4.3 – Resultados para as distribuições a posteriori: Casos A até H.
34
(µ(E)-E)/E)
Caso A
(µ(E1)-E1)/E1)
(µ(b1)- b1)/ b1)
4.62%
Caso B
1.14%
Caso C
0.1560%
Caso D
7.74%
Caso E
4.7%
Caso F
0.3550%
0.2140%
2.12%
Caso G
Caso H
1.6480%
29.87%
4.9483%
19.3140%
Tabela 4.4 – Desvio da média a posteriori em relação aos valores exatos: casos A até H.
A partir dos resultados apresentados nas tabelas (4.3) e (4.4), pode-se concluir
que, para os casos apresentados, em geral, processos de inferência que fazem uso de
medições experimentais cujo nível de poluição é maior, apresentam distribuições a
posteriori com maiores desvios-padrão, ou seja, o estado de conhecimento sobre os
parâmetros estimados é menor. Isto provavelmente implica no calculo de quantidades
físicas a posteriori, a partir de técnicas de propagação de incertezas, com maior
variabilidade e menor nível de confiança. Ainda mais, para o problema estudado,
quando o nível de ruído aumenta, a média a posteriori dos parâmetros estimados
apresenta maiores desvios em relação aos valores exatos. Se observa que este efeito é
amplificado na medida em que o número de parâmetros aumenta. De certa forma, estes
resultados eram esperados, pois analisando o teorema de bayes para densidades se
verifica que para uma mesma distribuição prior, funções de verossimilhança que
apresentam maiores desvios-padrão implicam em distribuições a posteriori mais
espalhadas. Portanto, se conclui que quanto mais preciso for o processo de medição
mais acuras serão os resultados obtidos a partir do procedimento de inversão estatística.
Também se verifica, a partir dos resultados das tabelas (4.3) e (4.4), os quais são melhor
ilustrados pelas figuras (4.5),(4.6) e (4,7), que um aumento no número de parâmetros
estimados implica em distribuições a posteriori com maiores desvios-padrão.
Este aumento ocorre devido a um maior do nível de dependência linear entre os
parâmetros, o que, por sua vez, implica em um maior número de combinações de
parâmetros que conduzem a altos valores da função de verossimilhança. Esta
característica pode ser entendida da seguinte forma, enquanto determinado parâmetro
35
caminha numa direção que tende a diminuir o valor da função de verossimilhança,
outros parâmetros podem caminhar em direções que anulem este efeito, tornando maior
a região no espaço dos parâmetros que apresenta valores razoáveis da função de
verossimilhança.
Figura 4.5 - Distribuição a posteriori para o parâmetro E: Caso A x Caso D.
Figura 4.6 - Distribuição a posteriori para o parâmetro E1: Caso B x Caso D.
36
Um ponto de destaque é o fato de que, dada uma resposta física do sistema, a
existência de dependência linear entre alguns grupos de parâmetros provavelmente tem
implicações no nível de correlação entre tais grupos na distribuição a posteriori
Portanto, a fim de realizar análises de correlação entre os parâmetros foram rodados
mais seis casos, os quais são definidos na tabela (4.5). Os resultados das análises estão
apresentados na tabela (4.6), (4.7) e (4.8), onde valores negativos de correlação podem
significar que o aumento no valor de determinado parâmetro tende a causar uma
diminuição no valor do outro. Vale destacar que apesar da distribuição a priori definida
no item anterior se basear na hipótese de parâmetros independentes, o processo de
inversão estatística conferiu uma estrutura correlacionada aos mesmos.
Figura 4.7 - Distribuição a posteriori para o parâmetro E1: Caso C versus Caso D.
Caso I
θ I = {E, E1 }
2
y obs
( SNR = 20 dB )
Caso J
θ J = {E,b1 }
2
y obs
( SNR = 20 dB )
Caso L
θ L = {E1 ,b1 }
2
y obs
( SNR = 20 dB )
Caso M
θ M = {E, E1 }
y 1obs ( SNR = 10 dB )
Caso N
θ N = {E,b1 }
y 1obs ( SNR = 10 dB )
Caso O
θ O = {E1 ,b1 }
y 1obs ( SNR = 10 dB )
Tabela 4.5 – Casos rodados para análise de correlação.
37
A partir destes resultados verificou-se que, para o problema em questão, o
aumento tanto no nível de ruído assim como no número de parâmetros estimados
implica em maiores valores de correlação. Como pode ser observado a partir das figuras
(4.8), (4.9) e (4.10), e tabela (4.7) as distribuições a posteriori apresentaram grande
alteração, dependendo dos parâmetros contidos dentro do vetor θ . Vale destacar que
nestes casos tanto o número de parâmetros estimados assim como o nível de
perturbação sobre o processo de medição foram os mesmos. Portanto, este aumento na
variabilidade das distribuições ocorreu devido a estrutura matemática do modelo, a qual
confere maior nível de dependência linear para algumas configurações de vetores θ do
que para outras, o que inclusive pode ser observado a partir das análises de correlação
apresentadas na tabela (4.6).
Corr (E,E1)
Caso I
Corr (E,b1)
Corr (E1,b1)
-0.0027
Caso J
0.0262
Caso L
0.0130
Caso D
0.0819
Caso M
-0.0233
0.2288
Caso N
0.2244
0.2236
Caso O
0.0865
Caso H
0.3639
0.8806
1.4430
Tabela 4.6 – Resultados da análise de correlação.
E [N/m]
b1 [Nsm2]
E1 [N/m]
Caso I
µ
0.9644
σ
0.0551
µ
5.9718
σ
0.0894
µ
σ
Caso M
1.1291
0.1514
5.7577
0.2617
Caso J
0.9659
0.1240
5.0288
0.2191
Caso N
1.2429
0.3678
5.4651
0.7067
Caso L
6.0058
0.1237
5.0904
0.1295
Caso O
5.6939
0.3217
4.7486
0.3429
Tabela 4.7 – Resultados das distribuições a posteriori: Casos I até O
38
Tabela 4.7 – Resultados das distribuições a posteriori: Casos I até O
Figura 4.8 – Distribuição a posteriori para parâmetro E: casos J x Caso I.
Figura 4.9 - Distribuição a posteriori para o parâmetro E1: Caso I x Caso L.
39
Figura 4.10 - Distribuição a posteriori para o parâmetro b1: Caso I x Caso L.
(µ(E)-E)/E)
(µ(E1)-E1)/E1)
(µ(b1)- b1)/ b1)
Caso I
3.56 %
0.47 %
Caso M
12.91 %
4.04 %
Caso J
3.41 %
0.58 %
Caso N
24.29 %
9.30 %
Caso L
0.10 %
0.0181 %
Caso O
5.10 %
5.03 %
Tabela 4.8 – Desvio da média a posteriori em relação aos valores exatos: casos I até O.
Outra forma de analisar a dificuldade para se obter informação sobre os
parâmetros de interesse a partir do procedimento de calibração do modelo é baseada no
cálculo de derivadas parciais da saída do problema direto. Este calculo possibilita tanto
realizar análises de sensibilidade, assim como avaliar a dependência linear entre os
parâmetros. A fim de comparar a sensibilidade dos parâmetros em relação ao modelo,
decidiu-se multiplicar o valor das derivadas parciais pelos respectivos parâmetros
40
constitutivos, como apresentado na equação (4.9), a qual define a matriz de
sensibilidade ou matriz jacobiano.
J =
 ∂y (t1 )
 ∂E × E

 ∂y (t 2 ) × E
 ∂E

M
 ∂y (t )
Ny

×E
 ∂E
∂y (t1 )
× E1
∂E1
∂y (t 2 )
E1
∂E1
M
∂y (t N y )
E1
∂E1
∂y (t1 )

× b1 
∂b1

∂y (t 2 )
× b1 

∂b1

M

∂y (t N y )
× b1 

∂b1
(4.9)
Figura 4.11 - Colunas Matriz de sensibilidade.
Na figura (4.11), que apresenta as colunas da matriz de sensibilidade relativas
aos parâmetros E, E1 e b1, se verifica que a saída do problema direto é mais sensível ao
parâmetro E1 do que para os outros parâmetros constitutivos. Considerando estes
resultados se conclui que, uma vez anulado o efeito da dependência linear, as
distribuições a posteriori relativas ao parâmetro E1 devem ser mais informativas do que
as distribuições relativas aos outros parâmetros. Para verificar esta conclusão, são
apresentados abaixo alguns resultados da tabela (4.3), onde se estimou somente um
parâmetro.
41
σ (E) / E
Caso A
σ (E1) / E1
σ (b1) / b1
4.52 %
Caso B
1.16 %
Caso C
1.48 %
Caso E
54.94 %
Caso F
15.35 %
Caso G
36.74 %
Tabela 4.9 – Efeito da sensibilidade dos parâmetros
4.3.2 - Propagação de incertezas via
Método de Monte Carlo
A necessidade de avaliar a confiabilidade e variabilidade das quantidades físicas
calculadas através de modelos computacionais justifica o emprego de problemas
inversos estocásticos no processo de estimação de parâmetros. Portanto, este item tem
por objetivo a aplicação do método de Monte Carlo, a fim de obter informação acerca
da variabilidade das predições numéricas fornecidas por modelos cujos parâmetros
foram estimados no item anterior. A resposta física do sistema (RFS), sobre a qual serão
efetuadas todas as análises, é dada pela velocidade do sistema, quando o mesmo é
submetido a condições iniciais nulas e uma história de forçamento que é dada pela
equação (4.5). Portanto se deseja calcular a função densidade de probabilidade a
•
posteriori desta quantidade π (z(t ) / θ, y obs ) . O critério de convergência aplicado baseiase na comparação gráfica das funções de densidade obtidas para diferentes números de
amostras. Entretanto deve-se mencionar que análises mais rigorosas devem se basear no
cálculo de métricas que levem em consideração quantidades estatísticas das
42
distribuições tais como média e desvio-padrão. Novamente, neste item se procurou
avaliar o efeito de algumas variáveis do problema inverso sobre a variabilidade das
predições. Estas variáveis são dadas pelo nível de poluição, número de parâmetros
estimados e nível de correlação entre os parâmetros. A fim de representar a
variabilidade das simulações computacionais se decidiu utilizar envelopes, os quais
englobam todos os valores de velocidade que se encontram dentro de uma região
definida entre dois desvios-padrões acima e abaixo da média a posteriori da velocidade.
Neste ponto, é importante destacar que os resultados aqui apresentados não associam
nenhum valor de probabilidade com os envelopes obtidos. Na tabela abaixo são
apresentadas as distribuições de probabilidade utilizadas como entrada para os
processos de propagação de incertezas nos três casos rodados.
Distribuições Consideradas
π (θ c / y obs )
Caso 1
π (θ D / y obs )
π (θ J / y obs )
Caso 2
π (θ L / y obs )
π (θ D / y obs )
Caso 3
π (θ H / y obs )
Tabela 4.10 – Distribuições a posteriori consideradas nos casos rodados.
Analisando a figura (4.12) se verifica que, para o problema investigado, a
variabilidade da resposta fornecida pelo modelo computacional aumenta quando o
número de parâmetros estimados é maior. Este fato ocorre devido a maior dependência
linear entre os parâmetros. A partir destes resultados, poder-se-ia pensar que uma
estratégia vantajosa seria sempre considerar a estimação probabilística de apenas um
parâmetro, o que tornaria possível o calculo de quantidades a posteriori com menores
variações. Entretanto tal estratégia estaria incorreta, simplesmente porque um modelo
computacional, em grande parte das vezes, é construído para fornecer predições além
dos dados experimentais com os quais foi calibrado [6], de modo que a definição de
quais parâmetros devem ser estimados, necessariamente envolve uma análise criteriosa
43
sobre o domínio de aplicação do modelo e dos fenômenos físicos que se destacam
dentro deste domínio, além da identificação dos parâmetros que representam de forma
mais significativa estes fenômenos através de análises de sensibilidade.
Figura 4.12 - Resultado do Caso 1.
Figura 4.13 - Resultado do Caso 2.
44
Figura 4.14 - Resultado do Caso 3 .
A partir da figura (4.13) se observa que, para o problema investigado, uma vez
constante, o número de parâmetros estimados e o nível de poluição nos dados
experimentais, a velocidade a posteriori apresenta maior variabilidade quando a
distribuição considerada para propagar as incertezas apresenta maior nível de
correlação.
Na figura (4.14) se observa que, para o problema investigado, o aumento no
nível de poluição do processo de medição implica em aumento na variabilidade da
velocidade a posteriori, entretanto tal aumento é menor que o observado nos casos 1 e
2. Este resultado mostra que apesar de os parâmetros serem estimados de forma
probabilística com o objetivo de incorporar as incertezas experimentais dentro do
modelo computacional, se observa que a maior fonte de variabilidade nas simulações
computacionais surge devido ao mau condicionamento numérico do problema inverso,
o qual ocorre devido a dependência linear entre os parâmetros.
A figura (4.15) além de apresentar o envelope de variabilidade para o caso D,
também apresenta o valor observado y obs da velocidade e a média a posteriori desta
quantidade. Portanto, se observa que o valor esperado para velocidade está em total
discordância com os valores medidos. Isto ocorre porque a função densidade de
45
probabilidade (F.d.p) é não simétrica para a maioria dos instantes de tempo. Para ilustrar
esta situação, nas figuras (4.16) e (4.17) são apresentadas comparações, para dois
instantes de tempo distintos, entre a F.d.p da velocidade, os valores observados e a
média a posteriori.
Figura 4.15 – Comparação entre resultados da propagação e dados observados: Caso D.
Na partir da figura (4.16) se verifica que, no instante de tempo de 5,78 segundos,
a distribuição de probabilidade da velocidade é bimodal e que o valor observado é
bastante diferente dos valores mais prováveis da distribuição. Também se observa que a
variabilidade da velocidade neste instante de tempo é elevada.
Na figura (4.17) se verifica que, no instante de tempo de 1,68 segundos, a
variabilidade da velocidade é menor, e que a distribuição de probabilidade da
velocidade é unimodal. Também se deve destacar que, neste instante de tempo, ocorreu
um bom nível de concordância entre o valor médio da velocidade e do valor observado.
Análises de convergência para estes processos de propagação de incertezas serão
apresentadas e discutidas no próximo item.
46
Figura 4.16 – F.d.p x Experimento, Caso D: tempo = 5.78s.
Figura 5.14 – Análise probabilística das predições: Caso D.
Figura 4.17 – F.d.p x Experimento, Caso D: tempo = 1.68s.
47
4.3.3 - Cadeias de Markov e Análises de
Convergência
Este item tem como foco a apresentação de cadeias de Markov relativas aos
processos de estimação apresentados no item (4.3.1), além da apresentação de análises
de convergência para distribuições a posteriori dos parâmetros e para velocidade a
posteriori que foi calculada no item anterior.
Portanto, as figuras (4.18), (4.19) e (4.20) apresentam cadeias de Markov
referentes, respectivamente, aos parâmetros E, E1 e b1, onde se observa uma maior
variabilidade das amostras tanto quando o número de parâmetros estimados aumenta,
assim como quando o nível de perturbação do processo de medição aumenta. Também
se observa que o aumento no número de parâmetros estimados causa maiores distorções
na cadeia. Nos casos analisados, com a abordagem escolhida para o presente trabalho, indicam
que o aumento do numero de parâmetros estimados tende a tornar o processo de
convergência da cadeia mais demorado. Por último deve-se destacar que o número de
amostras apresentadas nos gráficos abaixo não é necessariamente o considerado para
convergência da cadeia.
Figura 4.18 – Cadeias de Markov para parâmetro E.
48
Figura 4.19 – Cadeias de Markov para parâmetro E1.
Figura 4.20 – Cadeias de Markov para parâmetro b1.
49
Nas figuras (4.21), (4.22) e (4.23) são apresentadas cadeias de Markov
bidimensionais oriundas de processos de inferência onde se considerou valores de SNR
iguais a 10db e 20db. A partir destas figuras é possível observar que o ponto inicial das
cadeias está distante da região de maior probabilidade, entretanto esta é capaz de
esquecer posição inicial e caminhar para região de convergência da distribuição a
posteriori. Também se observa que o conjunto inicial de amostras não é representativo
da distribuição a posteriori de interesse e que, portanto não deve ser considerada para o
cálculo de quantidades estatísticas como média e variância. Outra constatação
importante é que regiões de confiança dos parâmetros estimados podem ser obtidas a
partir de subconjuntos de amostras da cadeia da Markov, no que configura uma
vantagem inerente do problema inverso estatístico quando comparado com a abordagem
clássica, onde a determinação de regiões de confiança se mostra uma tarefa mais
trabalhosa.
Figura 4.21 – Cadeias de Markov bidimensionais: Caso M x Caso I.
50
Figura 4.22 – Cadeias de Markov bidimensionais: Caso N x Caso J.
Figura 4.23 – Cadeias de Markov bidimensionais: Caso O x Caso L.
Nas figuras (4.24), (4.25) e (4.26) são apresentadas, para os parâmetros E, E1 e
b1, funções densidade de probabilidade, as quais ilustram a evolução desta quantidade
51
em função do número de amostras considerado. Observando estas figuras, se verifica
que não houve grandes mudanças nas F.d.ps quando o número de amostras aumentou
drasticamente. No presente trabalho este tipo de análise gráfica foi utilizado para
verificar a convergência da cadeia. Embora esta não seja a estratégia mais adequada, se
mostrou bastante útil para os problemas aqui abordados
Figura 4.24 – Convergência parâmetro E: Caso D.
Figura 4.25 – Convergência parâmetro E1: Caso D.
52
Figura 4.26 – Convergência parâmetro b1: Caso D
Figura 4.27 – Convergência velocidade (tempo=5.78): Caso D.
53
O mesmo tipo de análise gráfica foi utilizada para verificar a convergência da
velocidade a posteriori calculada no item anterior através do método de Monte Carlo.
As figuras (4.27) e (4.28) apresentam a evolução da F.d.p para os instantes de tempo de
5.78 segundos e 1.68 segundos, respectivamente. Novamente se observa que a forma
das F.d.ps não se alterara muito para um incremento de 90.000 amostras.
Figura 4.28 – Convergência velocidade (tempo=5.78): Caso D.
O número de amostras para o qual se considerou que as cadeias de Markov
convergiram variou bastante em função do número de parâmetros estimados. A título de
ilustração a ordem de grandeza destes valores é apresentada na tabela (4.11). A taxa de
aceitação das amostras geradas pelo algoritmo Metropolis Hastings é uma quantidade
importante e, portanto seus valores são apresentados na tabela (4.12)
1 parâmetro estimado
10.000 amostras
2 parâmetros estimados
50.000 amostras
3 parâmetros estimados
100.000 amostras
Tabela 4.11 – Ordem de grandeza do número de amostras para a convergência.
54
Caso A
43.71%
Caso H
69.74
Caso B
60.42%
Caso I
32.86%
Caso C
62.36%
Caso J
38.60%
Caso D
31.10%
Caso L
49.55%
Caso E
76.83%
Caso M
71.94%
Caso F
85.44%
Caso N
74.27%
Caso G
86.39%
Caso O
79.23%
Tabela 4.12 – Taxa de Aceitação do Método Metropolis Hastings.
4.3.4 - Estratégias Preliminares de
Validação
Até o presente momento foram discutidas as diversas incertezas envolvidas no
problema inverso. Entretanto as análises estiveram sempre relacionadas com aspectos
numéricos do processo de calibração, que por sua vez, também apresentaram impacto
no cálculo de quantidades a posteriori. Por outro lado, agora se deseja quantificar as
incertezas relativas a forma do modelo, e para isto se deve realizar processos de
validação. Este tipo de abordagem é essencial para creditação de modelos
computacionais. Neste ponto, deve-se destacar que não existe validação de modelo sem
dados genuinamente experimentais, e como as análises deste capítulo são baseadas em
pseudo-experimentos, não faz sentido validar o modelo obtido. Entretanto, visando
ilustrar estratégias preliminares de validação são apresentadas duas análises nas figuras
(4.29) e (4.30). Na primeira análise se compara a função densidade de probabilidade da
frequência natural do sistema com o valor exato da mesma. Na segunda análise se
compara a função densidade de probabilidade do fator de amortecimento do sistema,
com o valor exato do mesmo. Repare na figura (4.33) que o fator de amortecimento não
ocorre na moda da distribuição de probabilidade, mostrando que o valor mais provável
não é necessariamente o valor real da grandeza analisada. Os resultados encontrados
demonstram um nível razoável de correlação entre as distribuições de probabilidade e os
valores exatos das grandezas analisadas.
55
Figura 4.29 – Pseudo-validação com freqüência natural: Caso D.
Figura 4.30 – Pseudo-validação com fator de amortecimento: Caso D.
56
4.3.5 - Modelagem incorreta dos erros
experimentais
Neste item se deseja investigar o efeito da modelagem incorreta dos erros
experimentais no processo de estimação de parâmetros, e para tal são gerados dois
conjuntos de dados observados. Em y 3obs os erros são aditivos, gaussianos e
correlacionados de modo que termo de poluição da equação (4.4) é dado por.
v (t ) = αε (t ) + βε (t − 1)
(4.10)
Onde ε ~ N (0, σ corr ) segue uma distribuição normal com média zero e variância
2
σ corr = 1 , e os parâmetros α = 0.0285 e β = 0.0095 determinam o grau de correlação
do processo de poluição. Portanto a função de verossimilhança é dada pela equação
(3.14), onde a matriz de covariância dos erros experimentais é uma matriz tridiagonal
N y × N y como segue.
α 2 + β 2

 αβ
2
∑ e = σ corr  0

M

 0

αβ


O O O
M

O O O
0 

O O O
αβ 
L 0 αβ α 2 + β 2 
0
L
0
(4.11)
4
Por outro lado, o vetor de dados observados y obs
é gerado a partir de erros
aditivos, gaussianos, não correlacionados que apresentam média zero e variância
constante. Portanto a função de verossimilhança é dada pela equação (3.14) onde a
matriz de covariância dos erros é dada pela equação (4.7). Com o intuito avaliar
somente a influência da modelagem incorreta da estrutura dos erros, deve-se escolher o
4
valor de σ de modo que tanto y 3obs quanto y obs
a possuam a mesma magnitude de
4
perturbação. Portanto o valor do desvio-padrão do vetor de dados observados y obs
é
dado em cada instante de tempo por:
57
σ=
(α
2
)
2
+ β 2 σ corr
= 0.030
(4.12)
Os dois casos rodados nesta análise estimam além dos três parâmetros
constitutivos E, E1 e b1
o valor do desvio padrão σ dos erros, logo o vetor de
parâmetros é dado por θ = {E , E1 , b1 , Desvio − padrão}. A distribuição a priori escohida
para o desvio padrão é uniforme, onde os limites inferior e superior são dados,
respectivamente, por L4I = 0 e L4S = 0.1 . A função de verossimilhança utilizada nos dois
casos para o procedimento de inversão estatística se baseia na hipótese de erros não
correlacionados. No caso P o processo de inferência é realizado a partir do vetor de
dados observados y 3obs , constitutindo uma modelagem errada dos erros exsperimentais.
Enquanto que no caso Q o processo de inferência é realizado a partir do vetor de dados
4
observados y obs
, constituindo uma modelagem correta dos erros. As tabelas (4.13) e
(4.14) apresentam os resultados oriundos da aplicação da inferência bayesiana nos casos
P e Q. Se observa que, para o problema investigado, a modelagem incorreta dos erros
implica em maiores desvios das médias a posteriori em relação aos valores exatos dos
parâmetros. Esta indicação é bastante importante, pois muitas vezes a grandeza física
observada é definida no domínio da freqüência, e nestes casos o processamento de
sinais necessário para obter-la confere uma estrutura correlacionada aos erros
experimentais. Portanto o esforço envolvido na caracterização da estrutura dos erros
experimentais tem como recompensa resultados mais acurados.
Caso Q
Caso P
(Modelagem correta)
(Modelagem errada)
Val. exatos
µ
σ
µ
σ
E [N/m]
1
0.9820
0.3379
1.3027
0.2610
E1 [N/m]
6
6.2130
0.3927
6.5617
0.4739
b1 [Nsm2]
5
5.1915
0.8733
6.0966
0.8694
Desvio-Padrão
0.03
0.0310
0.0022
0.0293
0.0021
Tabela 4.13 – Resultados para distribuição a posteriori casos P e Q.
58
Caso Q
Caso P
(Modelagem correta) (Modelagem errada)
(µ(E)-E)/E)
1.8%
30.27%
(µ(E1)-E1)/E1)
3.7167%
9.3617%
(µ(b1)- b1)/ b1)
3.83%
21.9320%
Tabela 4.14 – Resultados para os casos P e Q.
As figuras (4.31), (4.32), (4.33) e (4.34) apresentam as distribuições a posteriori
dos parâmetros estimados para os casos P e Q. As figuras (4.35), (4.36) ilustram
cadeias de Markov bidimensionais obtidas a partir dos resultados dos casos P e Q. As
regiões definidas pelas cadeias de Markov mostram que a modelagem incorreta dos
erros experimentais não implica em maior variabilidade dos parâmetros estimados.
Figura 4.31 – Distribuição a posteriori, parâmetro E: caso P versus caso Q.
59
Figura 4.32 – Distribuição a posteriori, parâmetro E1: caso P versus caso Q.
Figura 4.33 – Distribuição a posteriori, parâmetro b1: caso P x caso Q.
60
Figura 4.34 – Distribuição a posteriori, Desvio-padrão: caso P x caso Q
Figura 4.35 – Cadeias de Markov Bidimensionais: E x E1.
61
Figura 4.36 – Cadeias de Markov Bidimensionais: E x b1.
62
5 - Viga Sanduíche
Este capítulo tem como foco a aplicação de estratégias probabilísticas
simplificadas que visam determinar um projeto ótimo de experimento para um modelo
de estrutura do tipo viga sanduíche, cujo núcleo é composto de um tipo de material
viscoelástico. O problema inverso é formulado no domínio da frequência e a abordagem
determinística é empregada através do método híbrido. São propostas quatro
configurações de sensores para compor o protocolo experimental. A abordagem de
Fisher é aplicada, e o critério de decisão é baseado na função de densidade
probabilidade da métrica D-optmally. A inferência bayesiana é aplicada e os resultados
são comparados para as possíveis configurações de sensores.
5.1 - Estrutura Investigada
Figura 5.1 – Modelo do Protótipo experimental.
A figura acima apresenta o modelo do protótipo experimental para a viga
sanduíche. As camadas de base e de restrição são de alumínio e as dimensões de tais
camadas se encontram na tabela (5.1). A camada intermediária corresponde a uma fita
viscoelástica cujas dimensões também se encontram na tabela (5.1). Quatro
acelerômetros foram fixados na parte superior da camada de restrição. A estrutura está
submetida, em suas extremidades, a condições de contorno do tipo pinada-pinada e um
shaker aplica um forçamento prescrito na posição L/4, onde L representa o
63
comprimento da viga. É importante destacar que como as camadas de base e de
restrição são totalmente elásticas, a fita viscoelástica é responsável por toda dissipação
de energia que ocorre no interior da viga sanduíche. Outro ponto a ser destacado, e que
os parâmetros físicos e geométricos aqui adotados foram retirados de um paper, onde
este sistema foi utilizado para efetuar a caracterização experimental de um material
viscoelástico específico [11].
Tabela 5.1 – Parâmetros físicos e geométricos da viga sanduíche.
A respeito da descrição matemática da estrutura sanduíche, foi utilizado um
modelo baseado no método de elementos finitos, o qual faz uso de um elemento
equivalente unidimensional a fim de evitar o aumento excessivo do custo computacional
[11,31]. O modelo constitutivo de variáveis internas é utilizado para representar o
comportamento mecânico do núcleo viscoelástico e, para fins de simplicidade, o modelo
utilizado considerou uma única variável interna, entretanto é importante destacar que
em procedimentos de caracterização, o número de variáveis internas também é um
parâmetro a ser estimado, o qual pode ser determinado através de processos de
competição de modelos.
64
5.2 – Experimento Simulado
Visando estimar os parâmetros constitutivos do modelo de variáveis internas se
deve realizar o experimento de calibração. Neste pseudo-experimento foi medida a
função de resposta em frequência H (ωj ) dos acelerômetros AC1, AC2, AC3 e AC4, onde
as respectivas posições destes acelerômetros são L/4, L/3, L/2 e 3L/4. A utilização de
pseudo-experimentos ou experimentos simulados é motivada pelo fato de permitir a
comparação dos resultados do processo de estimação com os valores exatos dos
parâmetros e, deste modo, avaliar a eficiência do método de inversão aplicado para o
problema analisado. Como mencionado no capítulo anterior, neste tipo de abordagem o
vetor de dados observados y obs e obtido a partir do resultado do problema direto
y = H (θ ) corrompido por um termo de perturbação como na equação abaixo.
y obs = y (θ) + v
(5.1)
onde o vetor de frequências ω é definido de [0 100] Hz, com N y = 801 pontos e taxa
de amostragem constante. Os componentes do vetor v , que tem dimensão N y × 1 , são
números oriundos de um processo aleatório v ~ N (0, σ 2 ) , onde o desvio-padrão é
avaliado de acordo com a equação (4.8), considerando uma razão sinal ruído (SNR=30
db), sendo que o sinal de referência é dado pela saída do problema direto. Neste ponto, é
importante ressaltar que, quando se trabalha com dados experimentais definidos no
domínio da frequência, a forma mais pragmática de efetuar a poluição deve levar em
consideração todo o processamento de sinais que ocorre em um analisador de espectros.
Para fins de estimação serão considerados cinco vetores de dados observados, os
quais são apresentados na tabela (5.2), onde n em H n (ωj ) representa o número do
acelerômetro considerado. Considerando a forma de poluição dos dados pseudoexperimentais, tem-se que a matriz de covariância dos erros para y 1obs é dada pela
equação (4.7), enquanto que para y 2 obs até y 5 obs a matriz de covariância dos erros é
dada pela equação (5.2).
65
y 1obs = {H 3 (ωj)}
T
y 2 obs = {H 3 (ωj ) H 4 (ωj )}
T
y 3 obs = {H 3 (ωj ) H1 (ωj )}
T
y 4 obs = {H1 (ωj ) H 2 (ωj )}
T
y 5 obs = {H 1 (ωj ) H 4 (ωj )}
T
Tabela 5.2 – Dados pseudo-experimentais utilizados nos processos de estimação.
σ n2 0 L L L 0 


M 
0 O O
 M O σ n2 O
M 
Σe = 

O σ m2 O M 
 M
 M
O O 0 


2
 0 L L L 0 σ m 
(5.2)
onde n e m em σ n2 e σ m2 representam os números dos dois acelerômetros considerados,
respectivamente.
5.3 - Projeto Ótimo de Experimento para
Viga Sanduíche
Este item tem como foco a aplicação de estratégias probabilísticas simplificadas que
visam definir um projeto ótimo de experimento para estrutura do tipo viga sanduíche.
Portanto o objetivo é determinar um vetor dentre os vetores y 2 obs até y 5 obs ,o qual
permite avaliar com maior confiança estatística os parâmetros constitutivos do modelo
de variáveis internas.
Portanto são analisadas quatro configurações de localização dos sensores, as
quais são denominadas por protocolos experimentais. Para fins de clareza estes
protocolos são apresentados na tabela abaixo, entretanto repare que estas posições são
as mesmas para as quais foram gerados os dados observados da tabela (5.2)
66
Protocolo 1
AC3
AC4
Protocolo 2
AC3
AC1
Protocolo 3
AC1
AC2
Protocolo 4
AC1
AC4
Tabela 5.2 – Configurações de acelerômetros avaliadas no Projeto ótimo de experimento.
O critério de escolha do protocolo experimental é baseado na abordagem de
Fisher, onde se procura determinar um protocolo experimental que maximize a métrica
D-optmally. Entretanto, como esta abordagem é baseada na hipótese de modelo linear,
deve-se determinar uma região no espaço dos parâmetros onde esta hipótese de modelo
∧
constitua uma boa aproximação. Esta região se encontra em torno do ponto ótimo θ e
pode ser encontrada a partir de uma abordagem determinística de problemas inversos.
Portanto é realizado um processo de estimação determinístico, o qual é
denominado de caso A, onde o vetor y 1obs é utilizado como experimento de calibração.
Considerando a forma de poluição dos dados pseudo-experimentais tem-se que a função
objetivo toma a forma da equação (3.16). A fim de encontrar o ponto de mínimo global
da função objetivo é utilizado um método de otimização híbrido. Este método tenta
combinar a capacidade de explorar o espaço dos parâmetros, que é característica dos
métodos heurísticos, com a velocidade de convergência observada nos métodos de
gradiente, a fim de encontrar o ponto de ótimo com um número menor de iterações.
A figura (5.2) ilustra as etapas envolvidas na abordagem híbrida,onde se observa
que primeiramente se deve definir uma região de busca para os parâmetros através dos
vetores XI e XS, os quais representam respectivamente, os limites inferiores e superiores
dos parâmetros. A segunda etapa consiste utilizar esta região previamente definida para
aplicar o método do Enxame de Partículas e obter como resultado um vetor θ EP . Por
último, o vetor θ EP é utilizado como chute inicial para o método do gradiente que é
aplicado, através da rotina “Lsqnonlin”, da plataforma computacional Matlab,
fornecendo, portanto o resultado do método de otimização híbrido θ hib . As
características do método do Enxame de Partículas e os resultados do problema inverso
determinístico são apresentados nas tabelas (5.4) e (5.5)
67
Figura 5.3 – Esquema para aplicação do método híbrido.
Número de iterações
500
Número de Partículas
50
Peso de inércia ( w )
Varia linearmente de 1 a 0
Parâmetro cognitivo ( c1 )
2
Parâmetro social ( c2 )
2
Espaço de busca
Parâmetro b1 → [0,100]
Parâmetros E , E1 → [0,50]
Tabela 5.4 – Características do Método do Enxame de partículas [15].
θ EP
θ hib
Exatos
E [N/m2]
1.2941
1.4028
1.403
E1 [N/m2]
15.2674
10.5526
10.558
b1 [Nsm2]
707.3791
536.95
537.42
Tabela 5.5 – Resultados do Problema Inverso Determinístico: Caso A.
Uma vez estimado o vetor de parâmetros para o caso A, pode-se utilizar o valor
da função objetivo avaliada em θ hib para calcular um estimador para variância dos erros
experimentais como segue.
σ2 =
S (θ hib )
N y − Nθ
(5.3)
68
onde N θ representa o número de parâmetros estimados. Dado que a matriz de covariância
dos erros para o vetor y 1obs é dada pela equação (4.7), também se pode calcular um
estimador para a matriz de covariância dos parâmetros como na equação (5.4) [35]. Os
procedimentos seguidos para calcular a matriz de covariância dos parâmetros são
ilustrados de forma esquemática na figura (5.3).
  ∂y  T −1  ∂y  
Σθ =   Σe  
  ∂θ 
 ∂θ  

−1
(5.4)
Figura 5.3 – Esquema para caçulo da matriz de covariância dos parâmetros.
Uma vez que o vetor de parâmetros estimados para o caso A está bastante
próximo do ponto de ótimo, poder-se-ia proceder com a aplicação da abordagem de
Fisher a fim de determinar o projeto ótimo de experimento. Entretanto, como
mencionado no capítulo 3, esta abordagem é local e, sendo assim, não leva em conta a
variabilidade dos parâmetros. Portanto, o que se deseja fazer é propor um procedimento
simples que permita incorporar as incertezas dos parâmetros no critério de decisão. A
figura (5.5) ilustra estes procedimentos de forma esquemática.
Figura 5.4 – Estratégia probabilística para projeto ótimo de experimento.
69
Portanto, primeiramente deve-se gerar k amostras oriundas de um processo
aleatório dado por θ ~ N (θ hib , Σ θ ) , onde Σ θ representa a matriz de covariância dos
parâmetros. Neste ponto deve-se destacar que foi necessário utilizar uma decomposição
de Cholesky para amostrar estes vetores devido ao fato da matriz Σ θ possuir termos
fora da diagonal.
Posteriormente, deve-se estimar a matriz de covariância dos erros experimentais
relativos aos vetores y 2 obs até
y 5 obs . Entretanto, como a princípio, não se tem
informação acerca destas quantidades se decidiu estimar os valores dos desvios-padrão,
em cada função de resposta em frequência de aceleração, como sendo igual a 3% da
magnitude máxima da saída do problema direto e utilizar a equação (5.2) para calcular
esta matriz.
A terceira etapa consiste em utilizar as k amostras geradas de θ para calcular k
valores determinísticos da métrica D-optmally e utilizar estes valores para estimar a
função densidade de probabilidade da métrica D-optmally. Este procedimento de ser
realizado para os quatro protocolos experimentais propostos.
Por último as distribuições de probabilidades obtidas devem ser comparadas e a
determinação do protocolo experimental que permite estimar os parâmetros com maior
confiança estatística é feita baseando-se no maior valor encontrado para a métrica Doptmally.
As figuras (5.5)-(5.8) apresentam os resultados das análises, a partir das quais se
define que o projeto ótimo de experimento como aquele proposto pelo protocolo
experimental 1.
Figura 5.5 – Métrica D-optmally: Protocolo 1.
70
Figura 5.6 – Métrica D-optmally: Protocolo 2.
Figura 5.7 – Métrica D-optmally: Protocolo 3.
71
Figura 5.8 – Métrica D-optmally: Protocolo 4.
Nos resultados apresentados acima, foram utilizadas 5000 amostras para calcular
as funções densidade de probabilidade da métrica D-optmally. A figura (5.9) mostra os
valores para esta distribuição considerando 1000 e 5000 amostras. Os resultados
indicam que houve uma convergência.
Figura 5.9 – Métrica D-optmally: Análise de convergência.
72
5.4 – Aplicação da Inferência Bayesiana
Este item tem como foco a aplicação da inferência bayesiana considerando os
quatro protocolos experimentais enumerados no item anterior, a fim de verificar as
indicações obtidas a respeito do projeto ótimo. Portanto serão rodados quatro casos, os
quais são apresentados na tabela abaixo.
Pseudo-experimento de calibração
Caso1
y 2 obs
Caso 2
y 3 obs
Caso 3
y 4 obs
Caso 4
y 5 obs
Tabela 5.6 – Casos rodados para Inferência bayesiana.
O vetor de parâmetros estimados é dado por θ = {θ1 , θ 2 , θ 3 , θ 4 , θ 5 }, onde θ 1 = E ,
θ 2 = E1 , θ 3 = b1 , θ 4 = β e θ 4 = β . Os parâmetros β 1 e β 2 são definidos pelas
equações abaixo.
β1 =
β2 =
σn
(5.6)
max(abs( H n (ωj )))
σm
(5.7)
max(abs( H m (ωj )))
onde n e m representam os números dos acelerômetros considerados para o calculo da
função de resposta em freqüência. Portanto o desvio-padrão dos erros também é
estimado de forma indireta. Considerando a forma como os dados pseudo-experimentais
foram gerados no item (5.2) se verifica que a função de verossimilhança toma a forma
da equação (3.14), onde a matriz de covariância toma a forma da equação (5.2).
Visando definir uma distribuição prior mais informativa é aplicada uma
estratégia que consiste em utilizar o conjunto de amostras θ k ,que foram geradas no item
73
anterior durante os procedimentos referentes ao projeto ótimo de experimento, para
calcular quantidades estatísticas tais como média e variância, e utilizar estas quantidades
para se obter um melhor entendimento das regiões de maior probabilidade dos
parâmetros estimados. Portanto a distribuição a priori é dada por uma distribuição
uniforme com parâmetros independentes, cujos limites inferiores e superiores são dados
,respectivamente, por LI e L S . Portanto, as três primeiras componentes deste vetor são
determinadas de acordo com as equações abaixo.
LiI = µ (θ i ) − 20σ (θ i )
i = 1: 3
(5.8)
LiS = µ (θ i ) + 20σ (θ i )
i = 1: 3
(5.9)
k
k
k
k
onde LiS e LiI representam, respectivamente, os limites superior e inferior da
distribuição a priori do parâmetro θ i . Os resultados obtidos são apresentados na tabela
abaixo em conjunto com os valores utilizados como prior para os parâmetros β 1 e β 2 .
L1I
1.3743
L1S
1.4313
L2I
10.4236
L2S
10.6812
L3I
534.9413
L3S
538.9604
L4I
0
L4S
0.1
L5I
0
L5S
0.1
Tabela 5.7 – Limites da distribuição prior uniforme.
Uma vez definida a distribuição prior e a função de verossimilhança, pode-se
aplicar a inferência bayesiana considerando experimentos simulados oriundos de quatro
protocolos experimentais como apresentado na tabela (5.6). Os resultados obtidos são
apresentados nas figuras (5.10) até (5.12).
74
Figura 5.10 – Comparações entre protocolos: E [N/m2].
Figura 5.11 – Comparações entre protocolos: E1 [N/m2].
75
Figura 5.13 – Comparações entre protocolos: b [Nsm2].
A partir dos resultados obtidos se verifica que quando os parâmetros são
estimados fazendo-se uso dos protocolos experimentais 3 e 4, os quais apresentaram os
menores valores paras distribuições da métrica D-optmally, as distribuições a posteriori
apresentaram maiores desvios-padrão. Por outro lado, quando os parâmetros são
estimados fazendo-se uso dos protocolos experimentais 1 e 2, as distribuições a
posteriori apresentam menores variâncias. Portanto estes resultados confirmam as
indicações obtidas nas análises de projeto ótimo de experimento.
76
6 – Conclusão
Este trabalho considerou o problema de estimação de parâmetros constitutivos
de um modelo para materiais viscoelásticos. O modelo de variáveis internas foi
utilizado para representar o comportamento mecânico destes materiais e o problema
inverso foi formulado utilizando dados pseudo-experimentais tanto no domínio do
tempo, assim como no domínio da frequência.
As diversas fontes de incerteza na computação científica foram discutidas e a
teoria de inversão estatística foi apresentada como uma alternativa para incorporar estas
incertezas no modelo computacional. A técnica de inferência bayesiana foi aplicada
com o objetivo de estimar os parâmetros na forma de distribuições de probabilidade
marginal a posteriori, e processos de propagação de incerteza via método de Monte
Carlo foram realizados a fim de inferir sobre a variabilidade das predições
computacionais dos modelos obtidos. Foram apresentadas estratégias preliminares que
visam avaliar a capacidade de predição dos modelos computacionais dentro de uma
filosofia de validação.
Um sistema massa-mola-amortecedor de um grau de liberdade anexado a um
componente viscoelástico foi investigado, e foram avaliadas as influências de diversas
variáveis do problema inverso, tais como, nível de poluição do processo de medição,
número de parâmetros estimados e modelagem estatística dos erros experimentais, sobre
os valores das distribuições a posteriori.
Por último, foi estudado o modelo de uma estrutura do tipo viga sanduíche. As
equações do modelo foram resolvidas através do método de elementos finitos e a saída
observada foi definida pela função de resposta em frequência do sistema. O projeto
ótimo de experimento foi definido a partir da abordagem de Fisher, onde a função
densidade de probabilidade da métrica D-optmally foi empregada como critério de
tomada de decisão. Foram adotadas estratégias no sentido de definir distribuições prior
mais informativas. Tais estratégias fizeram uso da abordagem determinística de
problemas inversos em conjunto com a hipótese de modelo linear em relação aos
parâmetros. A inferência bayesiana foi aplicada e os resultados obtidos confirmaram as
indicações obtidas a partir das análises de projeto ótimo de experimento.
77
Referências Bibliográficas
[1]- SANTOS, E. F., 2003, Atenuadores Viscoelásticos para Redução de Oscilações
Aeroelásticas de Edifícios Altos. Tese* de M.Sc., COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil.
[2] – RAO, M. D., “Recent Applications of Viscoelastic Damping for Noise Control in
Automobiles and Commercial Airplanes”, v.262, Journal of Sound and
eeeeeeeVibration, Academic Press, pp. 457-474, 2001.
[3] - PITELA, B. D. A., 2006, Investigação da Eficiência de Materiais Viscoelásticos
para Redução de Vibrações em Risers. Tese* de M.Sc., COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil.
[4] - OBERKAMPF W, L., ROY C, J., Verification and Validation in Scientific
Computing. 1 ed. New York, Cambridge, 2010.
[5] - GUERRA, G. M., 2011, Quantificação de Incertezas em Problemas de Interação
Fluido Estrutura via Método de Colocação Estocástica. Tese* de D.Sc.,
eeeeeeeCOPPE/UFRJ, Rio de Janeiro,RJ, Brasil.
[6] - American Society of Mechanical Engineers, ASME, 2006, “ Guide for Verification
and Validation in Computational Solid Mechanics”.
[7] - OZISIC, M. N., ORLANDE, H, R. B., Inverse Heat Transfer: Fundamentals and
Applications, Taylor and Francis, London, 2000
[8] - BECK, J.V., ARNOLD K. J., Parameter Estimation in Engineering and Science,
John Wiley & Sons, New York, 1977.
[9] - KAIPIO, J., SOMERSALO, R., Statistical and Computacional Inverse
Problems. Springer, 2004.
78
[10] - AGIRRE, M. M., ELEJABARRIETA, M, J., "Characterization and Modelling of
Viscoelastically Damped Sandwich Structures", v.52, International Journal of
Mechanical Sciences, Elsevier, pp. 1225-1233, 2010.
[11] - CASTELLO, D. A., ROCHINHA, F. A., ROITMAN, N., et al " Constitutive
Parameter Estimation of a Viscoelastic model with internal variables", v.22,
Mechanical Systems and Signal Processing, Elsevier, pp. 1840-1857, 2008.
[12] - WANG, J., QIN, Q. H., KANG, Y. L., et al., "Viscoelastic Adhesive Interfacial
Model and Experimental Characterization for Interfacial Parameters", v.42,
Mechanics of Materials, Elsevier, pp. 537-547, 2010.
[13]- KAUER, M., VUSKOVIC, V., DUAL, J., et al., "Inverse Finite Element
ddddddCharacterization of soft tissues", v.6, Medical Image Analysis, Elsevier, pp.
275-287, 2002.
[14] - Oius, D., "Combination of Standard Viscoelastic Model and Fractional Derivate
Calculos
to
the
Characterization
of
Polymers",
v.7,
Materials
Research Innovations, Springer-Verlag, pp. 42-46, 2003.
[15] - SOUZA, H. V. A. F., CASTELLO, D, A., ROITMAN, N., et al., " Parameter
Estimation and Aproximate Confidence Regions applied on Viscoelastic
Material Characterization". In: Proceedings of Inverse Problems, Design and
Optimization Symposium, João Pessoa, Agosto. 2010.
[16] - SOUZA, H. V. A. F., CASTELLO, D, A., ROITMAN, N., et al., " Constitutive
Parameter Estimation of Viscoelastic Materials". In: Proceedings of the XXX
Iberan-Latin-American Congress on Computational Methods in Engineering,
Armação de Búzios, Novembro. 2009.
[17] - CHRISTOPHER J. R; WILLIAM, L. O., 2011, “A comprehensive framework for
verification, validation, and uncertainty quantification in scientific computing”,
v.200, Comput. Methods Appl. Mech. Engrg, Elsevier, pp. 2131–2144, 2011.
79
[18] - COLOSIMO, B. M; CASTILLO, E., Bayesian Process Monitoring Control and
Optimization, CRC, 2007.
[20] - SHANBHAG, S., " Analytical Rheology of Blends of Linear and Star Polymers
using a Bayesian Formulation", v.49, Rheologica Acta, Springer-Verlag, pp.
411-422, 2010.
[21] - RUDOY, D., YUEN, S. G., HOWE, R. D., et al., " Bayesian Change Point
Analysis for Atomic Force Microscopy and Soft Material Indentation", v.59,
Applied Statistics, Cambridge, pp. 573-593, 2010.
[22] - BARD, Y., “Nonlinear Parameter Estimation”. New York, Academic
Press,1974.
[23] - EMERY, A. F., NENAROKOMOV, A. V., " Optimal Experiment Design", v.9,
Measurement Science and Technology, UK, pp. 864-876, 1998.
[24] - CASTELLO, D. A., ROCHINHA, F. A., " On The Optimum Experiment Design
Applied for Mechanical System Identification". In: Proceedings of the 2nd LNCC
Meeting on Computational Modelling, Rio de janeiro, Agosto. 2006.
[25]- HUAN, X., MARZOUK, M, Y., "Simulation-Based Optimal Bayesian
tttttttttttttExperimental Design for Nonlinear Systems", Submitted to Journal of
eeeeeeeComputational Physics, Agosto. 2011.
[26] – KHODJA, M, R., PRANGE, M, D., “Guided Bayesian Optimal Experimental
Design” v.26 , Inverse problems, 2010.
[27] - WINEMAN, A.S; RAJAGOPAL K. R., Mechanical Response of
Polymers.
Cambribge, 2000.
[28] - FLÜGGE, W., Viscoelasticity. 1ed. London, Blaisdell Publishing Company, 1967
80
[29] - CHRISTENSEN, R, M., Theory of Viscoelasticity: An Introduction, 2º ed.,
Academic Press, New York, 1994.
[30] - COSTI, F., 2006, Metodologia Numérica Aplicada a Viscoelasticidade em
Polímeros. Tese* de M.Sc., PUC-PR, Brasil.
[31] – CASTELLO, D. A., 2004, Modelagem e Identificação de Materiais
Viscoelásticos no Domínio do Tempo. Tese* de D.Sc., COPPE/UFRJ, Rio de
Janeiro,RJ,
Brasil.
[32] – BORGES, F., ROTMAN, N., MAGLUTA, C., et al. “ Vibration Reduction in
Steel Catenary Risers by the Use of Viscoelastic Materials”. In: Proceedings
of the 30th
International
Conference on
Ocean,
Offshore and
Artic
Engineering, Rotterdan, Junho. 2011.
[33] – BORGES, F., ROITMAN, N., MAGLUTA, C., “ A Model Validation Strategy
Applied for Viscoelastic Constitutive Model Calibration”. In: Proceedings of the
18º International Congress on Sound & Vibration
[34] – FOX, C., NICHOLLS, G. K., TAN, S, M., An Introduction to Inverse Problems.
2010.
[35]- SCHWAAB, M; BISCAIA, E. C; MONTEIRO, J. L; PINTO, J. C., 2007,
"Nonlinear parameter Estimation through particle swarm optimization”.
000000Chemical Engineering Science 63 (2008) 1542 – 1552.
[36] - COLAÇO, M. J; ORLANDE. H, R, B; DULIKRAVICH. G, S., 2006, “Inverse
and Optimization Problems in Heat Transfer”. J. of the Braz. Soc. of Mech. Sci.
& Eng. Vol. XXVIII.
[37] - PARSOPOULOS, K. E; VRAHATIS, M.N., 2002, “Recent Approaches to Global
Optimization Problems through Particle Swarm Optimization”, Natural
Computing, v.1, pp. 235-306.
81
[38]- KENNEDY, J; EBERHART, J., 1995, “Particle Swarm Optimization”. In:
gggyyyProceedings of the IEEE International Conference on Neural Networks, Perth,
Australia, vol.4, pp. 1942-1948.
[39] - SHI, Y; EBERHART, R., 1998, “A Modified Particle Swarm Optimizer”. InProc.
ddddddConference on Evolutionary Computation, Anchorage, Alaska, pp. 69-73.
[40] - SCHWAAB, M., 2005, “Avaliação de Algoritmos Heurísticos de Otimização em
Problemas de Estimação de Parâmetros”. Tese de M.Sc. COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil.
[41] - OBERKAMPF, W. L., BARONE, M. F., " Measures of Agreement between
Computational and Experiment: Validation Metrics", v.217, Journal of
Computational Physics, Elsevier, pp 5-36, 2006.
[42] – GHANEM, R., SPANOS, P., Stochastic Finite Elements: A Spectral Approach,
Springer-Verlag. 1991.
[43] – NAJM, H., “Uncertainty quantification and polynomial chaos techniques in
computational fluid dynamics”, v.41, Ann. Rev. Fluid Mech, pp. 35-52, 2009
82