Modelagem Estocástica e Quantificação de Incertezas Rubens Sampaio
[email protected]
Roberta de Queiroz Lima
[email protected]
Departamento de Engenharia Mecânica
DINCON 2015 PUC-Rio: DEM
Organização do curso 1a aula:
1
Probabilidade básica − Variáveis aleatórias (discretas e contínuas) − Vetores aleatórios − Independência
2
Função de variavéis aleatórias: − Soma de variáveis aleatórias independentes − Soma de variáveis aleatórias dependentes
3
Geração de amostras de variáveis e vetores aleatórios − Método da Transformada Inversa − Método da Rejeição − Condicionamento: regra da cadeia PUC-Rio: DEM
Organização do curso 2a aula:
1
Método de Monte Carlo − Lei dos Grandes Números − Teorema do Limite Central − Aproximação de integrais
2
Processos estocásticos (discretos e contínuos) − Processo de Bernoulli − Processo de Poisson − Processo uniforme − Processo exponencial
3
Modelagem estocástica de sistema dinâmico não linear: stick-slip. PUC-Rio: DEM
Método de Monte Carlo
Simulações estocásticas de um sistema:
1
constrói-se um modelo determinístico para o sistema;
2
constrói-se um modelo estocástico para o sistema; – histogramas – Princípio da Entropia Máxima (PEM)
3
aplica-se o método de Monte Carlo para obter-se inferências estatísticas sobre a resposta do sistema.
PUC-Rio: DEM
Método de Monte Carlo
PUC-Rio: DEM
Método de Monte Carlo
Amostras Variável aleatória X
(1) x =⇒ x(2) .. . (n) x
Processamento das Amostras y(1) =⇒ Estatísticas =⇒ y(2) .. . (n) y
µˆ =
σˆ 2 =
PUC-Rio: DEM
1 n ∑ yi n i=1
1 n ∑ (yi − µˆ )2 n i=1
Método de Monte Carlo Aproximação para o valor π .
Etapas: 1
Gerar m amostras distribuídas uniformemente no cubo. se:
2
3 4
x12 + x22 + x32 ≤ r
amostra está dentro da esfera
Calcular a razão, R, entre o número de amostras na esfera e o número total de amostras, m. Aproximação: πˆ = 6R. Verificar se πˆ está dentro de uma margem de erro prescrita. PUC-Rio: DEM
Método de Monte Carlo Se esse experimento for repetido n vezes, obtém-se uma sequência de aproximações:
πˆ1 , · · · , πˆn
Como saber se elas aproximam de fato π ?
O método de Monte Carlo é fundamentado em dois teoremas:
Lei dos Grandes Números: garante a convergência das aproximações obtidas através do método; Teorema do Limite Central: especifica como é a convergência. PUC-Rio: DEM
Método de Monte Carlo Se esse experimento for repetido n vezes, obtém-se uma sequência de aproximações: 14 12 10 8
πˆ1 , · · · , πˆn – Média do experimento, µ ? – Variância, σ 2 ?
6 4 2 0 3
3.05
3.1
π ˆ
3.15
3.2
Figura: Histograma normalizado
com n = 103 aproximações de π . Cada aproximação é calculada com m = 104 amostras no cubo. PUC-Rio: DEM
3.25
Lei dos grandes números Considere que os experimentos formam uma sequência de variáveis aleatórias Π1 , Π2 , · · · : independentes; identicamente distribuídas, com média µ e variância σ 2 .
Seja Sn = Π1 + Π2 + · · · + Πn . Tem-se que Sn /n converge em média quadrática para µ : v
u "
2 #
u
Sn S n
− µ = tE −→ 0, quando n −→ +∞ . −µ
n n PUC-Rio: DEM
Lei dos grandes números Pela lei dos grandes números, sendo Sn = Π1 + · · · + Πn :
µ Sn n
σ
2
Sn n
Sn 1 =E = (E[Π1 ] + E[Π2 ] + · · · + E[Πn ]) = µ = π , n n σ2 1 Sn 2 . = E ( − µ ) = 2 (var(Π1 ) + · · · + var(Πn )) = n n n
Note que: com o aumento de n, σˆ2 Sn decai com uma razão n proporcional a 1/n.
PUC-Rio: DEM
Lei dos grandes números Aproximações para µ e σ 2 podem ser obtidas através dos resultados πˆ1 , · · · , πˆn :
µˆ =
1 n ∑ πˆi , n i=1
1 n σˆ2 = ∑ (πˆi − µˆ )2 n i=1
A média e variância de Sn /n podem ser aproximadas por:
µˆ Sn = µˆ , n
1 n σˆ2 σˆ2 Sn = = 2 ∑ (πˆi − µˆ )2 . n n n i=1
PUC-Rio: DEM
Lei dos grandes números −3
2
x 10
σ2 Sn = n σˆ2 Sn
σ2 n
n
var(Sn /n)
1.5
1
0.5
0 0
20
40
n
60
80
100
Figura: Variância de Sn /n em função do número de repetições n do experimento.
PUC-Rio: DEM
Lei dos grandes números Experimento: gerar amostras de X1 e X2 com densidades N (0, 1) e calcular a média das amostras geradas. Repete-se o experimento n = 1, 000 vezes. Varia-se o tamanho, m, da amostra: 0.4
m = 100
m = 1,000 0.4 0.2 E[X2 ]
E[X2 ]
0.2 0 −0.2
0 −0.2
−0.4 −0.4 −0.2
0
E[X1 ]
−0.4 −0.4 −0.2
0.2
m = 100,000
0.4
0.4
0.2
0.2 E[X2 ]
E[X2 ]
m = 10,000
0 −0.2 −0.4 −0.4 −0.2
0 0.2 E[X1 ]
0 −0.2
0 0.2 E[X1 ]
−0.4 −0.4 −0.2
0 0.2 E[X1 ]
PUC-Rio: DEM
Teorema do limite central Seja Π1 , Π2 , · · · uma sequência de variáveis aleatórias independentes e identicamente distribuídas, com média µ e variância σ 2 . Seja Sn = Π1 + Π2 + · · · + Πn e seja: Zn =
Sn − µ n √ σ n
Quando n −→ +∞, Zn converge em probabilidade para uma variável aleatória Z com distribuição cumulativa de probabilidade gaussiana de média nula e variância um.
PZn (Zn ≤ x) −→ PZ (Z ≤ x) quando n −→ +∞, ∀x onde PZ é contínua.
PUC-Rio: DEM
Teorema do limite central Verifica-se que com probabilidade 0.954 uma amostra de N (0, 1) está no intervalo (−2, 2): Pr (−2 < Zn < 2) = 0.954 . Sn − µ n √ < 2 = 0.954 . Pr −2 < σ n Assim, com probabilidade 0.954, µ estará no intervalo: √ √ Sn − 2σ n Sn + 2σ n x) = Pr(X1 > x, · · · , Xn > x) = . l l−x n PY1 (x) = Pr(Y1 ≤ x) = 1 − Pr(Y1 > x) = 1 − . l n(l − x)n−1 PY1 (x) (x) = 1[0,l] (x) . pY1 (x) = dx ln Observe que embora as distribuições deX1 · · · Xn sejam uniformes, a distribuição de Y1 não é uniforme! A ordenação feita modificou a distribuição de Y1 com relação a X1 · · · Xn . PUC-Rio: DEM
Processo uniforme
Distribuição de Yk : O evento (Yk ≤ x) significa que pelo menos k dos n pontos escolhidos estão no intervalo [0, x]. Esse evento pode ser realizado de várias formas diferentes dependendo do número de pontos que de fato estão em [0, x]. Depois de algumas contas, determina-se a densidade: pYk (x) = 1[0,l] (x)
xk−1 (l − x)n−k . ln
PUC-Rio: DEM
Processo uniforme Implementação em MATLAB: sorteio e ordenação de n = 4 pontos no intervalo [0, 5]. 1
0.4
Histograma pY1
0.8
Histograma pY2
0.3
0.6 0.2 0.4 0.1
0.2 0
0 0
1
2
y1
3
4
5
0.4
0
1
2
y2
3
4
5
3
4
5
0.8
Histograma pY3
0.3
Histograma pY4
0.6
0.2
0.4
0.1
0.2
0
0 0
1
2
y3
3
4
5
0
1
2
y4
Figura: Histogramas normalizados contruídos com 105 amostras.
PUC-Rio: DEM
Processo uniforme Distribuição conjunta de Yj e Yk : Sendo x < y:
j−1 n x (y − x)k−j−1 (l − y)n−k pYj Yk (x, y) = j − 1, 1, k − j − 1, 1, n − k ln Distribuição conjunta de Y1 · · · Yn :
pY1 ···Yn (x1 , · · · , xn ) =
n! ln
0,
se 0 ≤ x1 < x2 < · · · < xn < l outros casos .
PUC-Rio: DEM
Processo uniforme Implementação em MATLAB: sorteio e ordenação de n = 4 pontos no intervalo [0, 5]. Y 1 , Y2
Y1 , Y3 0.25 0.2
0.4
0.15 0.2
0.1 0
0 0
0.05
1 2 2 4
0
0 0
3 4
2
2
4
y1
6 5
y2
4 6
6
y1
y3
Y 1 , Y4
Y2 , Y3 0.25
0.5 0.4
0.2
0.3
0.15
0.2
0.1
0.1
0.05 0 6
0 6 6
4
4
2
y4
2 0 0
y1
6
4
4
2
y3
2 0 0
y2
Figura: Histogramas normalizadosPUC-Rio: contruídos com 105 amostras. DEM
Espaçamento no processo uniforme Considere as distâncias entre os pontos sucessivos tomados na ordem crescente:
L1 = Y 1 L2 = Y 2 − Y 1
L3 = Y 3 − Y 2 .. . Ln+1 = l − Yn
⇒
⇒ ⇒
Y2 = L1 + L2 Y3 = Y2 + L3 = L1 + L2 + L3
Yn = l − Ln+1 .
PUC-Rio: DEM
Espaçamento do processo uniforme O espaçamento entre Yk e Yk+1 é notado por Lk+1 , e o espaçamento entre Yn e l é notado por Ln+1 , logo: Yk = L1 + · · · + L k . Observe que os espaçamentos não são independentes: L1 + · · · + Ln+1 = l . Porém são equidistribuidos. Sendo L1 = Y1 , então a distribuição de Lk é a mesma que a de Y1 , ∀k, e a densidade é: n pLk (x) = 1[0,l] (x) n (l − x)n−1 . l
PUC-Rio: DEM
Espaçamento do processo uniforme O valor esperado do espaçamento: E [Lk ] =
Z l n 0
x n (l − x)n−1 dx . l
Uma maneira mais elegante da fazer essa conta é usando o fato de que E [L1 ] = · · · = E [Ln+1 ] e L1 + · · · + Ln+1 = l. Assim: l . n+1 Fica agora evidente como calcular: E [Lk ] =
E [Yk ] = E [L1 ] + · · · + E [Lk ] =
PUC-Rio: DEM
kl . n+1
Processos estocásticos contínuos Um processo estocástico, X , é uma função: X :
Ω × T −→ R (w, t) 7−→ X (t, w)
Quando T é um conjunto contínuo, por exemplo R, o processo estocástico pode ser interpretado como um número infinito de variáveis aleatórias indexadas pelo parâmetro t. Fixado um valor de t = tj , a variável aleatória X (tj , w) tem a função distribuição de probabilidade cumulativa: PX (tj ,w) (x) = Pr(X (tj , w) < x) e uma função densidade de probabilidade pX (tj ,w) . PUC-Rio: DEM
Processos estocásticos contínuos A função densidade do processo estocástico é a função densidade conjunta de todas as variáveis aleatórias indexadas pelo parâmetro t.
Especificação completa um processo esto´castico: É necessário conhecer a densidade de probabilidade conjunta de todas as suas variáveis aleatórias. Problema: Como t assume infinitos valores, é impossivel especificar completamente um processo estocástico! Alternativa: Fazer hipóteses adicionais sobre independência das variáveis aleatórias indexadas por t. PUC-Rio: DEM
Propagação de incertezas em um sistema dinâmico
PUC-Rio: DEM
Propagação de incertezas em um sistema dinâmico Sistema real =⇒ Modelo matemático determinístico
Equação da dinâmica: m¨x(t) + kx(t) = f (t), onde f é a fora¸ de atrito (Coulomb) entre o bloco e a esteira f (t) = nµ sgn(v(t) − x˙ (t)).
PUC-Rio: DEM
Propagação de incertezas em um sistema dinâmico
PUC-Rio: DEM
Oscilador harmônico com atrito seco Equação da dinâmica: m¨x(t) + kx(t) = f (t), onde f é a fora¸ de atrito (Coulomb) entre o bloco e a esteira f (t) = nµ sgn(v(t) − x˙ (t)).
Considerando que v e µ são constantes no tempo, e introduzindo y = x˙ e z = x ωn , onde ωn é a frequencia natural do sistema, quando y > v quando y < v
nµ 2 y + z+ =c mωn nµ 2 2 = c. y + z− mωn 2
PUC-Rio: DEM
Oscilador harmônico com atrito seco v positivo → um ponto de equilíbrio em
,
nµ mωn , 0
v negativo → um ponto de equilíbrio em − mnωµn , 0 .
PUC-Rio: DEM
Stick-slip O stick ocorre quando y = v e quando a força de atrito que varia de acordo com f (t) = k x(t), está confinada no intervalo −fmax < f < fmax , onde fmax = µ n. Dessa forma, o segmento horizontal y = v e |x| ≤ dem ao stick.
µn k
correspon-
O stick só pode ocorrer no regime transiente. No regime estacionário, o diagrma de fase é escrito como um círculo e só terá slip. PUC-Rio: DEM
Stick-slip Vamos observar o que acontece quando a esteira tem velocidade descontínua:
Cada mudança do sinal da velocidade pode ser interpretada como uma renicialização do sistema. O parâmetro 2ωb determina a frequência com que o sistema é renicializado. PUC-Rio: DEM
Esteira tem velocidade descontínua
Figura: Diagrama de fase em regime permanente com (a) r = 0.1 e (b) r = 0.4 .
PUC-Rio: DEM
Esteira tem velocidade descontínua
Figura: Diagrama de fase em regime permanente com (a) r = 1.0 e (b) r = 1.3.
PUC-Rio: DEM
Esteira tem velocidade descontínua
Figura: Diagrama de fase em regime permanente com (a) r = 2.0 e (b) r = 2.5.
PUC-Rio: DEM
Modelo estocástico da velocidade da esteira
Fixado um intervalo de tempo [0, T] para analíse da resposta do sistema, o número de descontinuidades na velocidade da esteira é modelada como um processo de Poisson, N (T).
pN (T) (n) =
(λ T)n exp (−λ T) . n!
PUC-Rio: DEM
Modelo estocástico da velocidade da esteira Fixado o intervalo [0, T], e sorteado o número de descontinuidades da velocidade da esteira, n, os instantes das descontinuidades são modelados como um processo uniforme. Sorteio: n pontos são escolhidos aleatoriamente e de forma independente no intervalo [0, T]. Ordenação: a1 , a2 , · · · , an
=⇒
t1 ≤ t2 ≤ · · · ≤ tn
PUC-Rio: DEM
Modelo estocástico da velocidade da esteira Implementação em Matlab: histograma do número de mudanças da esteira. 0.05
λ × T = 5×20
0.04 0.03 0.02 0.01 0 60
80 100 120 numero de mudancas
140
Figura: Histograma do número de mudanças da esteira contruído com 3, 000 amostras.
PUC-Rio: DEM
Uma realização do sistema Simulação em Matlab: uma realização das velocidades da esteita e do bloco. 1
x˙ [m/s]
0.5 0 −0.5 −1 0
1
2 t [s]
3
4
Figura: Realização da velocidade da esteira e da velocidade do bloco. PUC-Rio: DEM
Propagação de incertezas Objetivo: calcular estatísticas da resposta do sistema.
Número de sticks; S Instantes em que começam os sticks: K1 , K2 , · · · , Ks Duração dos sticks: D1 , D2 , · · · , Ds
Instantes em que começam os slips: L1 , L2 , · · · , Ls
Para os cálculos das estatísticas, 3,000 simulações de Monte Carlo foram realizadas. PUC-Rio: DEM
Número de sticks Implementação em Matlab: histograma do número de sticks. 0.12
λ × T = 5×20
0.1 0.08 0.06 0.04 0.02 0 20
30 40 numero de sticks
50
Figura: Histograma do número de sticks contruído com 3, 000 amostras. PUC-Rio: DEM
Instantes em que começam os sticks Implementação em Matlab: histogramas dos instantes em que começam os sticks.
Figura: Histogramas dos instantes em que começam os sticks contruídos com 3, 000 amostras.
PUC-Rio: DEM
Duração dos sticks Implementação em Matlab: histograma das durações dos sticks.
Figura: Histogramas das durações dos sticks contruídos com 3, 000 amostras.
PUC-Rio: DEM
Instantes em que começam os slips Implementação em Matlab: histogramas dos instantes em que começam os slips.
Figura: Histogramas dos instantes em que começam os slips contruídos com 3, 000 amostras. PUC-Rio: DEM
Artigo
PUC-Rio: DEM
Livros publicados / em preparação
2012
2013
2014
PUC-Rio: DEM
2016
Modelagem Estocástica e Quantificação de Incertezas Rubens Sampaio
[email protected]
Roberta de Queiroz Lima
[email protected]
Departamento de Engenharia Mecânica
DINCON 2015 PUC-Rio: DEM
