Introdução a PINNs

3 Redes informadas pela física 3.2 Problemas transientes 3.4 Saiba mais

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

3.3 Problemas inversos

Talvez os problemas mais interessantes que podemos resolver com PINNs sejam os problemas inversos, em que queremos determinar parâmetros desconhecidos de um modelo matemático a partir de dados experimentais. Ao contrário de redes baseadas em dados, as PINNs podem incorporar o conhecimento do modelo matemático, o que permite estimar parâmetros desconhecidos com menos dados e, em alguns casos, com maior precisão.

Vamos definir o seguinte problema inverso: dada uma equação diferencial

\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% L(u;\lambda)=f,\leavevmode\nobreak\ \bm{x}\in D\subset\mathbb{R}^{n},

(3.54)

queremos estimar o parâmetro $\lambda$ a partir de um conjunto de amostras

\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \mathcal{D}:=\left\{\left(\bm{x}^{(s)},u^{(s)}\right)\right\}_{s=1}^{n_{s}},

(3.55)

da solução $u=u(\bm{x})$ da equação (3.54). Aqui, $L$ é um operador diferencial, $\lambda\in\mathbb{R}$ é o parâmetro a determinar e $f$ é uma função fonte conhecida. Assumimos conhecidas as condições iniciais ou de contorno, conforme o caso seja estacionário ou transiente.

Denotemos a PINN com parâmetro a determinar por

\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \tilde{u}=\mathcal{N}(\bm{x};(\bm{\theta},\lambda)),

(3.56)

em que $\tilde{u}$ é a solução estimada da equação diferencial (3.54), $\bm{\theta}$ são os parâmetros da rede neural (pesos e biases) e $\lambda$ é o parâmetro a determinar. Ou seja, uma rede neural em que $\lambda$ é apenas mais um parâmetro a ser treinado juntamente com os pesos e biases da rede.

O treinamento da PINN é feito minimizando a função de perda que, agora, além de considerar o resíduo da equação e as condições inicial/contorno, também considera o erro entre a solução estimada $\tilde{u}$ e as amostras $\mathcal{D}$ da solução $u$ . A função de perda é dada por

\displaystyle\varepsilon_{\lambda}:=\underbrace{\frac{1}{n_{s,\text{in}}}\sum_% {s=1}^{n_{s,\text{in}}}\left|f-L\left(u_{\text{in}}^{(s)};\lambda\right)\right% |^{2}}_{\text{res\'{\i}duo}}+\text{c.i./c.c.}+\text{amostras},

(3.57)

em que os termos c.i. e c.c. são definidos conforme as condições iniciais e de contorno. O termo das amostras é dado pelo erro médio quadrático entre a solução estimada $\tilde{u}$ e as amostras $\mathcal{D}$ da solução $u$ .

\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \frac{1}{n_{s}}\sum_{s=1}^{n_{s}}\left|\tilde{u}^{(s)}-u^{(s)}\right|^{2}.

(3.58)

Observemos que a dependência do parâmetro $\lambda$ na função de perda se dá pelo resíduo da equação diferencial, que é calculado em pontos de amostragem internos no domínio $D$ .

Observação 3.3.1.(Outros tipos de problemas inversos)

A aplicação de PINNs não é restrita a problemas inversos de determinação de parâmetros. Podemos também estimar funções desconhecidas, como por exemplo, a função fonte $f$ ou o operador $L$ . Neste caso, podemos representar a função desconhecida por uma rede neural e treinar a rede juntamente com a rede que representa a solução $u$ . Consultemos, por exemplo, o E.3.3.3.

Exemplo 3.3.1.(Taxa populacional)

O modelo de Fisher⁶⁶6Ronald Aylmer Fisher, 1890-1962, biólogo inglês. Fonte: Wikipédia: Ronald Fisher. descreve a evolução de uma população com taxa de crescimento $\lambda$ e difusão espacial. Trata-se de um modelo de reação-difusão, em que a população cresce logisticamente e se difunde no espaço. Aqui, vamos assumir a equação de Fisher com condições de contorno de Dirichlet não homogêneas:

	$\displaystyle u_{t}=u_{xx}+\lambda u(1-u),\leavevmode\nobreak\ (t,x)\in(0,t_{f% })\times(0,1),$		(3.59)
	$\displaystyle u(0,x)=1/\left(1+e^{\sqrt{\frac{\lambda}{6}}x}\right)^{2},% \leavevmode\nobreak\ x\in[0,1],$		(3.60)
	$\displaystyle u(t,0)=1/\left(1+e^{-\frac{5}{6}\lambda t}\right)^{2},$		(3.61)
	$\displaystyle u(t,1)=1/\left(1+e^{\sqrt{\frac{\lambda}{6}}-\frac{5}{6}\lambda t% }\right)^{2}.$		(3.62)

Queremos determinar a taxa de crescimento populacional $\lambda$ assumindo o seguinte conjunto de amostras da solução $u$

\displaystyle\mathcal{D}=\left\{\left(\left(t^{(s)},x^{(s)}\right),u^{(s)}% \right)\right\}_{s=1}^{n_{s}},

(3.63)

com $\left(t^{(s)},x^{(s)}\right)\in\{0.1,0.2,0.3\}\times\{0.25,0.5,0.75\}$ e $u^{(s)}=u_{a}\left(t^{(s)},x^{(s)}\right)$ , sendo

\displaystyle u_{a}(t,x)=1/\left(1+e^{\sqrt{\frac{\lambda}{6}}x-\frac{5}{6}% \lambda t}\right)^{2}

(3.64)

a solução analítica do problema [2]. Consulte a Figura 3.4 para uma ilustração da solução e dos dados conhecidos para o valor esperado de $\lambda=6$ .

Refer to caption — Figura 3.4: Solução PINN versus analítica para $\lambda=6$ . Os pontos correspondem às amostras da solução conhecidas.

O Código 20 é uma implementação de uma PINN para resolver este problema inverso.

Código 20: ex_pinn_fisher.py

⬇

1import torch

3# dispositivo GPU/CPU ?

4device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

5print(f'device: {device}')

7# modelo

9n_h = 4 # num de camadas escondidas

10n_n = 50 # num de neurônios por camada escondida

11act_fun = torch.nn.Tanh() # função de ativação

13model = torch.nn.Sequential()

14# camada de entrada/

15model.add_module('layer_1', torch.nn.Linear(2,n_n))

16model.add_module('fun_1', act_fun)

17# camadas escondidas

18for i in range(n_h-1):

19 model.add_module(f'layer_{i+2}', torch.nn.Linear(n_n,n_n))

20 model.add_module(f'fun_{i+2}', act_fun)

21# camada de saída

22model.add_module(f'layer_{n_h+1}', torch.nn.Linear(n_n,1))

24# adiciona param lambda

25model.lmbda = torch.nn.Parameter(12.0 * torch.rand(1, device=device))

27model.to(device) # envia modelo para o dispositivo

29# treinamento

31## optimizador

32optim = torch.optim.Adam(model.parameters(), lr=0.001)

34## num de amostras por época

35ns_t = 20 # num de amostras no tempo

36ns_x = 40 # num de amostras no espaço

38## num max épocas

39nepochs = 50_000

41## critério de parada (lambda)

42window_tol = 1e-3

43window_size = 100

44lmbda_window = [model.lmbda.item()]

46# dados

47ts_data = torch.tensor([0.1, 0.2, 0.3], device=device)

48xs_data = torch.tensor([0.25, 0.5, 0.75], device=device)

49Ts_data, Xs_data = torch.meshgrid(ts_data, xs_data, indexing='ij')

50S_data = torch.hstack((Ts_data.reshape(-1,1), Xs_data.reshape(-1,1)))

51lmbda_exp = torch.tensor(6.0, device=device)

52Us_data = 1.0/(1.0 + torch.exp(torch.sqrt(lmbda_exp/6.0)*S_data[:,1:2] \

53 - 5.0/6.0*lmbda_exp*S_data[:,0:1]))**2

56# laço de treinamento

57for epoch in range(nepochs):

59 # amostras de pontos internos

60 t = torch.rand(ns_t, 1, device=device)

61 x = torch.rand(ns_x, 1, device=device)

62 T, X = torch.meshgrid(t.flatten(), x.flatten(), indexing='ij')

63 S = torch.hstack((T.reshape(-1,1), X.reshape(-1,1)))

64 S.requires_grad = True

66 # propagação nos pontos internos

67 U = model(S)

69 # derivadas de primeira ordem

70 DU = torch.autograd.grad(U, S,

71 grad_outputs=torch.ones_like(U),

72 create_graph=True)[0]

73 DUDT = DU[:,0:1] # du/dt

74 DUDX = DU[:,1:2] # du/dx

76 # derivadas de segunda ordem

77 # d2u/dx^2

78 D2UDX = torch.autograd.grad(DUDX, S,

79 grad_outputs=torch.ones_like(DUDX),

80 create_graph=True)[0][:,1:2]

82 # param

83 lmbda = model.lmbda

85 # resíduo

86 R = DUDT - D2UDX - lmbda*U*(1.0 - U)

88 # função de perda: resíduo

89 loss_res = torch.mean(R**2)

91 # condição inicial

92 Sci = torch.hstack((torch.zeros_like(x), x))

93 Uci = model(Sci)

94 Uexp_ci = 1.0/(1.0 + torch.exp(torch.sqrt(lmbda/6.0)*x))**2

95 loss_ci = torch.mean((Uci - Uexp_ci)**2)

97 # condições de contorno

99 # contorno x = 0

100 Scc1 = torch.hstack((t, torch.zeros_like(t)))

101 Ucc1 = model(Scc1)

102 Ucc1_exp = 1.0/(1.0 + torch.exp(-5.0/6.0*lmbda*t))**2

103 loss_cc1 = torch.mean((Ucc1 - Ucc1_exp)**2)

104

105 # contorno x = 1

106 Scc2 = torch.hstack((t, torch.ones_like(t)))

107 Ucc2 = model(Scc2)

108 Ucc2_exp = 1.0/(1.0 + torch.exp(torch.sqrt(lmbda/6.0) \

109 - 5.0/6.0*lmbda*t))**2

110 loss_cc2 = torch.mean((Ucc2 - Ucc2_exp)**2)

111

112 # dados

113 Udata_est = model(S_data)

114 loss_data = torch.mean((Udata_est - Us_data)**2)

115

116 # função de perda total

117 loss = loss_res + loss_ci + loss_cc1 + loss_cc2 + loss_data

118

119

120 if (epoch % 100 == 0):

121 print(f'{epoch}: lambda: {model.lmbda.item():.2e} loss_res: {loss.item():.2e} ci: {loss_ci.item():.1e} cc1: {loss_cc1.item():.1e} cc2: {loss_cc2.item():.1e}')

122 torch.save(model, f'model.pt')

123

124 # critério de parada

125 if (epoch < window_size):

126 lmbda_window.append(model.lmbda.item())

127 else:

128 lmbda_window.append(model.lmbda.item())

129 lmbda_window.pop(0)

130 lmbda_var = max(lmbda_window) - min(lmbda_window)

131

132 if lmbda_var < window_tol and epoch > 5000:

133 print(f'Critério de parada atingido na época {epoch}.')

134 print(f'\tlambda: {model.lmbda.item():.2e}')

135 break

136

137 # backward

138 optim.zero_grad()

139 loss.backward()

140 optim.step()

3.3.1 Teste de sensibilidade

Uma característica importante de todo o método para problemas inversos é sua sensibilidade em relação aos dados. Mais precisamente, queremos saber como a estimativa do parâmetro $\lambda$ muda quando adicionamos ruídos aos dados.

Exemplo 3.3.2.

Vamos considerar o mesmo problema do Exemplo 3.3.1, mas agora vamos adicionar ruídos aos dados. Para isso, vamos gerar amostras da solução $u$ com ruído gaussiano aditivo

\displaystyle u^{(s)}=u_{a}\left(t^{(s)},x^{(s)}\right)+\sigma\epsilon^{(s)},% \leavevmode\nobreak\ s=1,\dotsc,n_{s},

(3.65)

em que $\epsilon^{(s)}\sim\mathcal{N}(0,1)$ é uma variável aleatória com distribuição normal padrão e $\sigma>0$ é o desvio padrão do ruído. Vamos considerar diferentes valores de $\sigma$ e verificar como a estimativa do parâmetro $\lambda$ muda com o aumento do ruído.

Para isso, vamos fazer a seguinte alteração no Código 20:

Código 21: ex_pinn_fisher_ruido.py

⬇

1# dados com ruído

2sigma = 0.1 # desvio padrão do ruído

3Us_data = 1.0/(1.0 + torch.exp(torch.sqrt(lmbda_exp/6.0)*S_data[:,1:2] \

4 - 5.0/6.0*lmbda_exp*S_data[:,0:1]))**2

5Us_data += sigma * torch.randn_like(Us_data)

A Tabela 3.1 mostra o menor valor e o maior valor estimado do parâmetro $\lambda$ em 5 tentativas para cada valor de $\sigma$ . O percentual de erro é calculado como a diferença entre o maior e o menor valor estimado de $\lambda$ em relação ao valor esperado de $\lambda=6$ . Aqui, fizemos apenas 5 tentativas para cada valor de $\sigma$ e descartamos o menor e o maior valor das estimativas. Uma análise mais completa necessitaria de mais tentativas para obter uma estimativa mais fidedigna do erro. O erro relativo foi calculado como

\displaystyle\text{Erro}=\frac{\tilde{\lambda}_{\max}-\tilde{\lambda}_{\min}}{% \lambda_{\text{esperado}}}\times 100\%.

(3.66)

Tabela 3.1: Estimativa do parâmetro

\lambda

(média de 5 tentativas) para diferentes níveis de ruído

\sigma

$\sigma$	$\tilde{\lambda}$	Erro
$0.00$	$[6.0,6.0]$	$0.0$ %
$0.01$	$[5.9,6.1]$	$1.7$ %
$0.05$	$[5.3,6.4]$	$12$ %
$0.1$	$[6.0,6.7]$	$12$ %

3.3.2 Exercícios

E. 3.3.1.

Crie uma PINN para estimar o parâmetro $\lambda>0$ do seguinte problema de Poisson

	$\displaystyle-\lambda u^{\prime\prime}=\pi^{2}\operatorname{sen}(\pi x),% \leavevmode\nobreak\ 0<x<1,$		(3.67)
	$\displaystyle u(0)=u(1)=0,$		(3.68)

dado o conjunto de amostras $D=\{(x^{(s)},u^{(s)})\}_{s=1}^{3}$

D=\left\{\left(\frac{1}{6},2\right),\left(\frac{1}{4},\sqrt{2}\right),\left(% \frac{1}{3},\sqrt{3}\right)\right\}.

(3.69)

Faça uma análise de sensibilidade adicionando ruídos aos dados e verificando como a estimativa do parâmetro $\lambda$ varia com o aumento do ruído.

E. 3.3.2.

Crie uma PINN para resolver o seguinte problema inverso de Poisson

	$\displaystyle-\Delta u=\lambda(2-x^{2}-y^{2}),\leavevmode\nobreak\ (x,y)\in D=% (-1,1)^{2},$		(3.70)
	$\displaystyle u=0,(x,y)\in\partial D,$		(3.71)

para estimar o parâmetro $\lambda>0$ dado que $u(0,0)=25$ . Faça uma análise de sensibilidade adicionando ruído ao dado. O aumento do conjunto de dados melhora a estimativa do parâmetro $\lambda$ ? Justifique sua resposta.

E. 3.3.3.

Crie uma PINN para resolver o problema inverso de estimar a fonte de calor para o modelo físico

	$\displaystyle u_{t}=u_{xx}+e^{-t}f(x),\leavevmode\nobreak\ (t,x)\in(0,1)^{2},$		(3.72)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ x\in[0,1],$		(3.73)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ t\in[0,1],$		(3.74)

em que $f(x)$ é a função a ser estimada dado que a solução no tempo final $t=1$ é

\displaystyle u(1,x)=e^{-1}\operatorname{sen}(\pi x),\leavevmode\nobreak\ x\in% [0,1].

(3.75)

E. 3.3.4.

Crie uma PINN para estimar os parâmetros escalares $\lambda_{1},\lambda_{2}>0$ do seguinte sistema de EDPs

	$\displaystyle-\Delta u+\lambda_{1}v=f_{1},\leavevmode\nobreak\ (x,y)\in D=(-1,% 1)^{2},$		(3.76)
	$\displaystyle-\Delta v+\lambda_{2}u=f_{2},\leavevmode\nobreak\ (x,y)\in D=(-1,% 1)^{2},$		(3.77)
	$\displaystyle u=v=0,\leavevmode\nobreak\ (x,y)\in\partial D,$		(3.78)

com fontes

	$\displaystyle f_{1}(x,y)=4-2x^{2}-2y^{2}+\operatorname{sen}(\pi x)% \operatorname{sen}(\pi y),$		(3.79)
	$\displaystyle f_{2}(x,y)=2\pi^{2}\operatorname{sen}(\pi x)\operatorname{sen}(% \pi y)+2(1-x^{2})(1-y^{2}).$		(3.80)

tendo, ainda, o conjunto de amostras da solução

\displaystyle\mathcal{D}=\left\{\left((0,0),1,0\right),\left((0.5,0.5),0.5625,% 1\right),\left((-0.5,0.5),0.5625,-1\right)\right\},

(3.81)

em que cada elemento é da forma $\left((x^{(s)},y^{(s)}),u^{(s)},v^{(s)}\right)$ . Faça uma análise de sensibilidade adicionando ruído aos dados e verifique como as estimativas de $\lambda_{1}$ e $\lambda_{2}$ variam com o aumento do ruído.

E. 3.3.5.

Considere o escoamento incompressível de Navier-Stokes no domínio $D=(0,2\pi)^{2}$

	$\displaystyle u_{t}+uu_{x}+vu_{y}=-p_{x}+\nu\left(u_{xx}+u_{yy}\right),% \leavevmode\nobreak\ (t,x,y)\in(0,1]\times D,$		(3.82)
	$\displaystyle v_{t}+uv_{x}+vv_{y}=-p_{y}+\nu\left(v_{xx}+v_{yy}\right),% \leavevmode\nobreak\ (t,x,y)\in(0,1]\times D,$		(3.83)
	$\displaystyle u_{x}+v_{y}=0,\leavevmode\nobreak\ (t,x,y)\in(0,1]\times D,$		(3.84)

em que $u=u(t,x,y)$ e $v=v(t,x,y)$ são as componentes da velocidade, $p=p(t,x,y)$ é a pressão e $\nu>0$ é a viscosidade cinemática do fluido, sendo este o parâmetro que queremos estimar. Condições inicial e de contorno são dadas pelo vórtice de Taylor⁷⁷7Brook Taylor, 1685 - 1731, matemático britânico. Fonte: Wikipédia: Brook Taylor.-Green⁸⁸8George Green, 1793 - 1841, matemático britânico. Fonte: Wikipédia: George Green., uma solução analítica das equações de Navier-Stokes

	$\displaystyle u_{a}(t,x,y)=-\cos(x)\operatorname{sen}(y)\,e^{-2\nu t},$		(3.85)
	$\displaystyle v_{a}(t,x,y)=\operatorname{sen}(x)\cos(y)\,e^{-2\nu t},$		(3.86)
	$\displaystyle p_{a}(t,x,y)=-\frac{1}{4}\left(\cos(2x)+\cos(2y)\right)e^{-4\nu t}.$		(3.87)

Crie uma PINN

\displaystyle\left(\tilde{u},\tilde{v},\tilde{p}\right)=\mathcal{N}(t,x,y;(\bm% {\theta},\nu)),

(3.88)

para estimar $\nu$ . Estude diferentes conjuntos de dados a partir da solução analítica e verifique como a estimativa de $\nu$ muda com o aumento do ruído nos dados. Como se comporta a sensibilidade da estimativa de $\nu$ para valores esperados de $\nu$ muito pequenos? Justifique sua resposta.

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Política de uso de dados