math.math-ds math.math-oc math.math-pr math.math-st nlin spi.meca.stru spi.meca.solid spi.meca.geme
Eber Dantas 
  
Americo Cunha 
  
Problemas inversos para identificação de parâmetros em sistemas dinâmicos

 

Problemas Inversos para Identificação de Parâmetros em Sistemas Dinâmicos Eber Dantas de Sá Paiva Americo Barbosa da Cunha Junior eber.paiva@uerj.br americo@ime. • Sistemas dinâmicos → engenharia, medicina, física, finanças, etc.

• Identificação de parâmetros é essencial em qualquer modelagem.

• Problemas inversos: encontrar os parâmetros que fazem o modelo simular os dados corretamente.

ferramenta problema inverso modelo calibrado método numérico otimizar J(x)

dados replicados

Problema direto e inverso.

Dificuldades:

• Eficiência do método numérico.

• Alta não linearidade de J(x) e múltiplas soluções de sua otimização.

• Ruído nos dados. Incertezas de modelo.

Objetivo

• Identificação de parâmetros para dois problemas de engenharia com dados experimentais de bancada.

• Comparação direta do método da Entropia-Cruzada (CE) com Algoritmo Genético (GA), avaliando eficiência e aplicabilidade.

Método da Entropia-Cruzada

Fundamentos e vantagens

Ideia principal

Problema de otimização → Problema de estimação de evento raro

X ∼ f (x, v) , J(X) ≥ γ ≈ max x J(x) é um evento raro Como? Minimizando a Divergência KL entre δ(x -x * ) e f (•, v).

Divergência de Kullback-Leibler

D(g, f ) = R n g(x) ln g(x) dx - R n g(x) ln f (x) dx
• Algoritmo meta-heurístico de busca global com caráter estocástico.

• Amostragem de Importância e técnicas de minimização de variância.

• Não requer diferenciabilidade nem continuidade de J(x).

• Garantia teórica de convergência se apenas um ótimo global existir.

• Fácil implementação. Fórmulas analíticas para família exponencial.

• Poucos e claros parâmetros de controle.

• Popularização dos modelos de ordem reduzida e técnicas de learning.

Algoritmo CE

Algoritmo para maximização de J(X). Estimação em multiníveis.

1 Defina número de amostras N , tamanho N e dos grupos elite t , distribuição de amostragem f (•, v), hiperparâmetro inicial v 0 e t = 0.

2 Atualize t = t + 1.

3 Gere amostras X 1 , • • • , X N ∼ f •, v t-1 (iid) e compute J(X i ).
4 Separe em t as N e amostras que fornecem os maiores J(X i ).

5 Atualize v t = arg max v X k ∈ t ln f (X k , v).
6 Repita 2 -5 até que um dos critérios de parada ocorra.

Vibração em Viga Euler-Bernoulli

Problema direto

ρ A ∂ 2 w ∂t 2 + ∂ 2 ∂x 2 E I ∂ 2 w ∂x 2 = 0
Suportes fletores flexíveis: k t 1 e k t 2 .

Método dos elementos finitos:

[K] φ = ω 2 [M ] φ Cálculo dos auto-pares (ω, φ).
Esquema para viga EB.

Problema inverso

Parâmetros:

x = {k t 1 , k t 2 }
Quantidades de interesse: for the mode shape of the natural frequency of interest. The developed formulation does not require numerical modelling and provides a useful tool to the design phase. The above referred papers deal with slender beams, typically with l/D > 10. As a consequence, the inertial rotation energy of an infinitesimal element dx is negligible compared to its translational energy, justifying the applicability of the Bernoulli-Euler beam theory. Sometimes, when this principle is violated, the problem has to be treated using Timoshenko beam theory. 13,14 In this study, one shall use Rayleigh's quotient, with shape functions that take into account the torsional flexibility of the supports, while assuming infinite rigidity along the vertical direction (Fig. 2). The algorithm will also include the effect of the mass and torsional inertia of the attached discs. Such an approach aims at a more precise model that more accurately represents the real dynamic behaviour of the system, rather than assuming simply supported or perfectly clamped ends. After showing that these ideal cases are far from the experimental values, the updating procedure for identifying the value of the torsional stiffness of both supports will be explained. 

Ω(x) = {ω 1 , ω 2 }

Theoretical Background

To obtain an analytical model that reproduces the physical system (Fig. 1) one considers the Rayleigh's quotient, using the Bernoulli-Euler's beam theory to define the shape functions. On a first attempt, two ideal cases have been tried: simply supported (SS) and clamped (CC) beam at both ends. However, the results for those ideal cases differ significantly from the experimental ones. Therefore, it is proposed to build a model (Fig. 2) considering that the supports of the beam present different torsional stiffnesses k t i [elastic supports (ES)], whose values are iterated until the error between experimental and numerical results satisfies a given tolerance τ . Regarding Rayleigh's method, the energy dissipation is neglected and therefore the Principle of Conservation of Energy holds:

(T + V ) = 0 ⇒ T max = V max (1) 
where T and V are the kinetic and potential energies, respectively. Through the quantification of the work performed by the elastic forces, in relation to the equilibrium position of the beam, the potential energy is

V = 1 2 l 0 M(x, t)dθ (2)
where M(x, t) is the bending moment given by M(x, t) = E(x)I(x)∂ 2 w(x, t)/∂x 2 , l is the length of the beam and θ is the rotation angle given by θ (x, t) = ∂w(x, t)/∂x. E(x) is the Young modulus and I(x) the second moment of area of the beam; w(x, t) represents the lateral displacement of the beam in relation to its equilibrium position.

At each undamped natural frequency, the time variation of w(x, t) can be shown to be harmonic and therefore the displacement response is given by the harmonic variation of a shape function φ(x):

w(x, t) = φ(x) sin(ωt + ϕ) (3) 
As w(x, t) max = φ(x) it follows that

V max = 1 2 l 0 E(x)I(x) φ ′′ (x) 2 dx (4)
On the other hand, the kinetic energy is given by

T = 1 2 l 0 ẇ(x, t) 2 dm ( 5 
)
From Eq. 3, ẇ(x, t) max = ωφ(x). As dm = ρ(x)A(x)dx, where ρ(x) is the density of the material and A(x) the cross sectional area of the beam, one has

T max = ω 2 1 2 l 0 ρ(x)A(x) (φ(x)) 2 dx (6) 2 
Experimental Techniques (2011) © 2011, Society for Experimental Mechanics Bancada de testes para viga EB.

Achar x * que minimize: Dinâmica Torcional de Drillstring Problema direto 

J(x) = || Ω exp -Ω(x)|| 2 + λ || x || 2
J 1 θ1 + d 1 ( θ1 -θ2 ) + k 1 (θ 1 -θ 2 ) + T r1 ( θ1 ) = 0 , J 2 θ2 + d 1 ( θ2 -θ1 ) + d 2 ( θ2 -θ3 ) + k 1 (θ 2 -θ 1 ) + k 2 (θ 2 -θ 3 ) = 0 , d 2 ( θ3 -θ2 ) -k 2 (θ 2 -θ 3 ) = η (K T i -C m η θ3 -T f -J 3 θ3 /η) , L di dt + R i + K E η θ3 = k p (ω ref -θ3 ) + k i t 0 (ω ref -θ3 ) dt , T r1 ( θ1 ) =    N 1 r 1 µ s θ1 /ω s for | θ1 | < ω s , N 1 r 1 µ k + (µ s -µ k ) e -v b | θ1 | sign( θ1 ) for | θ1 | ≥ ω s , Integração numérica → θ1 , θ2 → Filtro ( 2 
ρ i (x) = cov( Θi (x), Θexp i ) σ Θi (x) σ Θexp i
Bancada de testes para drillstring.

Achar x * que maximize: 

J(x) = ρ 1 (x) + ρ 2 (x) -λ|x -x r | Resultados Calibração, periodograma θ1 

Contribuições

• Verificação da viabilidade do método da Entropia-Cruzada -Acurácia suficiente em algoritmo simples.

-Convergência mais rápida do que o GA em milhares.

-Melhor performance que o GA em baixa amostragem.

• Identificação de apoios flexíveis na vibração de uma viga EB.

• Calibração de um modelo para dinâmica torcional de um drillstring.

Trabalhos futuros:

• Avaliar o método CE em problemas mais difíceis, com sofisticações: atualização dinâmica de variância, mistura Gaussiana, etc

• Avaliar o método CE em modelos de ordem reduzida para problemas de engenharia com alta demanda computacional.

Agradecimentos

Agradecimentos aos professores Bruno Cayres (CEFET), Hans I. Weber (PUC-Rio) e Tiago A. N. Silva (Universidade NOVA de Lisboa) pela colaboração, coautoria e fornecimento dos dados.

  uerj.br NUMERICO -Núcleo de Modelagem e Experimentação Computacional Introdução

Figure 1

 1 Figure 1 System under study (MFS, SpectraQuest, Inc.).

Figure 2

 2 Figure 2 Model of a beam with elastic torsional supports and a rigid disc at x = a.

2 (

 2 Hz) ω 1 (Hz) ω 2 (Hz) k t

  Hz), FFT → periodogramas Θ1 , Θ2 . Problema inverso Parâmetros: x = {k 1 , k 2 , d 1 , d 2 , µ s , µ k } Quantidades de interesse: { Θ1 , Θ2 } Coeficiente de Pearson: