Introdução a PINNs

3 Redes informadas pela física 3 Redes informadas pela física 3.2 Problemas transientes

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

3.1 Problemas estacionários

Como um problema fundamental em equações diferenciais parciais estacionárias, vamos estudar como podemos criar uma PINN para resolver o problema de Poisson¹¹1Siméon Denis Poisson, 1781 - 1840, matemático francês. Fonte: Wikipédia: Siméon Denis Poisson.

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% -\Delta u=f,\leavevmode\nobreak\ \bm{x}\in\mathcal{D},$		(3.6)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u=g,\leavevmode\nobreak\ \bm{x}\in\partial D,$		(3.7)

no domínio $\mathcal{D}\in\mathbb{R}^{2}$ , com dadas fonte $f$ e condições de contorno de Dirichlet²²2Johann Peter Gustav Lejeune Dirichlet, 1805 - 1859, matemático alemão. Fonte: Wikipédia: Johann Peter Gustav Lejeune Dirichlet..

Vamos considerar uma PINN dada por uma MLP

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

(3.8)

tendo como entrada as coordenadas $(x_{1},x_{2})$ e como saída a estimativa da solução do problema (3.6)-(3.7).

Como uma PINN, o treinamento da rede é feito minimizando uma função de perda que incorpora a equação diferencial e a condição de contorno. A função de perda é dada por

\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \varepsilon:=\underbrace{\frac{1}{n_{s,in}}\sum_{s=1}^{n_{s,in}}\left|f^{(s)}+% \Delta\tilde{u}^{(s)}\right|^{2}}_{\text{res\'{\i}duo}}+\underbrace{\frac{1}{n% _{s,cc}}\sum_{s=1}^{n_{s,cc}}\left|\tilde{u}^{s}-g^{(s)}\right|^{2}}_{\text{c.% c.}},

(3.9)

onde $f^{(s)}$ , $\tilde{u}^{(s)}$ e $g^{(s)}$ são os valores da fonte, da solução estimada e da condição de contorno no $s$ -ésimo ponto de amostragem $\bm{x}^{(s)}=(x_{1}^{(s)},x_{2}^{(s)})\in\overline{\mathcal{D}}$ . A cada época de treinamento, conjuntos randômicos de pontos de amostragem internos $\left\{\bm{x}^{(s)}\right\}_{s=1}^{n_{s,in}}\subset\mathcal{D}$ e de pontos de contorno $\left\{\bm{x}^{(s)}\right\}_{s=1}^{n_{s,cc}}\subset\partial\mathcal{D}$ são gerados com distribuição uniforme.

Exemplo 3.1.1.(Equação de Poisson)

Vamos criar uma PINN para resolver o problema de Poisson (3.6)-(3.7) com fonte

\displaystyle f(x_{1},x_{2})=2\pi^{2}\operatorname{sen}(\pi x_{1})% \operatorname{sen}(\pi x_{2}).

(3.10)

Para fins de comparação, a solução analítica do problema é

\displaystyle u(x_{1},x_{2})=\operatorname{sen}(\pi x_{1})\operatorname{sen}(% \pi x_{2}).

(3.11)

Observemos que a solução é a mesma função que estudamos no Exemplo LABEL:cap_mlp_sec_eqpoisson:ex:mlp_apfun_2d. Desta forma, aqui, adaptamos o código do Exemplo LABEL:cap_mlp_sec_eqpoisson:ex:mlp_apfun_2d para criar uma PINN para resolver o problema de Poisson. A diferença é que, agora, a função de perda da rede é baseada no resíduo da equação de Poisson e na condição de contorno, em vez de ser baseada na diferença entre a solução estimada e a solução analítica.

Refer to caption — Figura 3.1: Comparação entre a solução PINN (iso-cores) e a solução analítica (iso-linhas). Exemplo de pontos de amostragem para treinamento da PINN (pontos vermelhos).

O Código 17 é uma implementação de uma PINN para este problema de Poisson. Usa uma MLP de arquitetura $2-2\times 30-1$ (duas camadas escondidas, cada uma com 30 neurônios) e função de ativação tangente hiperbólica. O método de otimização Adam foi escolhido para o treinamento. O critério de parada do treinamento é baseado na função de perda total, que deve ser menor que uma tolerância especificada por um número consecutivo de épocas.

A Figura 3.1 mostra uma comparação entra a solução PINN (iso-cores) obtida pelo código e a solução analítica (iso-linhas). A cada época de treinamento, os pontos de amostragem são gerados de forma randômica por uma distribuição uniforme. A figura mostra um exemplo de pontos de amostragem (pontos vermelhos) para treinamento da PINN.

Código 17: py_pinn_poisson.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 = 2 # num de camadas escondidas

10n_n = 30 # 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))

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

25# treinamento

27## fonte

28def f(X):

29 return 2.0 * torch.pi**2 * torch.sin(torch.pi*X[:,0:1]) * \

30 torch.sin(torch.pi*X[:,1:2])

32## optimizador

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

35## num de amostras por época

36ns_x0 = 40

37ns_x1 = 40

39## num max épocas

40nepochs = 50_000

42## tolerância

43tol = 1e-3

44persist = 5

46# contador de persistência

47n_conv = 0

49# laço de treinamento

50for epoch in range(nepochs):

52 model.train()

54 # amostras de pontos internos

55 x0 = 2.0 * torch.rand(ns_x0, 1, device=device) - 1.0

56 x1 = 2.0 * torch.rand(ns_x1, 1, device=device) - 1.0

57 X0, X1 = torch.meshgrid(x0.flatten(), x1.flatten(), indexing='ij')

58 S = torch.hstack((X0.reshape(-1,1), X1.reshape(-1,1)))

59 S.requires_grad = True

61 # propagação nos pontos internos

62 U = model(S)

64 # derivadas de primeira ordem

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

66 grad_outputs=torch.ones_like(U),

67 create_graph=True)[0]

68 DUDX0 = DU[:,0:1] # du/dx0

69 DUDX1 = DU[:,1:2] # du/dx1

71 # derivadas de segunda ordem

72 # d2u/dx0^2

73 D2UDX0 = torch.autograd.grad(DUDX0, S,

74 grad_outputs=torch.ones_like(DUDX0),

75 create_graph=True)[0][:,0:1]

76 # d2u/dx1^2

77 D2UDX1 = torch.autograd.grad(DUDX1, S,

78 grad_outputs=torch.ones_like(DUDX1),

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

81 # laplaciano

82 L = D2UDX0 + D2UDX1

84 # fonte

85 F = f(S)

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

88 loss_res = torch.mean((F + L)**2)

90 # amostras de pontos de contorno

92 # contorno x1 = -1

93 Scc1 = torch.hstack((x0, -torch.ones_like(x0)))

94 Ucc1 = model(Scc1)

95 loss_cc1 = torch.mean((Ucc1 - 0.0)**2)

97 # contorno x0 = 1

98 Scc2 = torch.hstack((torch.ones_like(x1), x1))

99 Ucc2 = model(Scc2)

100 loss_cc2 = torch.mean((Ucc2 - 0.0)**2)

101

102 # contorno x1 = 1

103 Scc3 = torch.hstack((x0, torch.ones_like(x0)))

104 Ucc3 = model(Scc3)

105 loss_cc3 = torch.mean((Ucc3 - 0.0)**2)

106

107 # contorno x0 = -1

108 Scc4 = torch.hstack((-torch.ones_like(x1), x1))

109 Ucc4 = model(Scc4)

110 loss_cc4 = torch.mean((Ucc4 - 0.0)**2)

111

112 # função de perda total

113 loss = loss_res + loss_cc1 + loss_cc2 + loss_cc3 + loss_cc4

114

115 if (epoch % 100 == 0):

116 print(f'{epoch}: loss_res: {loss.item():.2e}, cc1: {loss_cc1.item():.2e}, cc2: {loss_cc2.item():.2e}, cc3: {loss_cc3.item():.2e}, cc4: {loss_cc4.item():.2e}')

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

118

119 if (loss.item() < tol):

120 n_conv += 1

121 print(f'\t{epoch}: loss: {loss.item():.2e}, {n_conv}/{persist}')

122 if (n_conv >= persist):

123 print(f'\t{epoch}: critério de parada atingido.')

124 break

125 else:

126 n_conv = 0

127

128 # backward

129 optim.zero_grad()

130 loss.backward()

131 optim.step()

Observemos que o treinamento de uma PINN naturalmente envolve uma função de perda composta por múltiplos termos, cada um representando uma parte do problema físico. Isto pode levar a dificuldades de convergência, especialmente quando os termos da função de perda têm escalas diferentes. Uma abordagem comum para lidar com isso é a utilização de pesos de penalidade para equilibrar a contribuição de cada termo na função de perda. Consultemos o E.LABEL:exer:mlp_poisson_penalidade.

Outro aspecto importante para se observar, é que, em geral, não é possível determinar a relação entre o erro de uma solução aproximada e o valor de seu resíduo em uma equação diferencial. É claro que um valor de resíduo igual a zero implica que a solução aproximada é uma solução exata da equação diferencial. No entanto, um valor de resíduo pequeno não garante que a solução aproximada seja próxima da solução exata.

Exemplo 3.1.2.

Vamos verificar o comportamento do erro da solução estimada pela PINN em relação ao valor da função de perda da rede treinada no Exemplo 3.1.1. Para isso, podemos fazer uma pequena moficação no Código 17 para calcular e armazenar, a cada época, os valores da função de perda e do erro $L^{2}$ relativo da solução estimada em relação à solução analítica.

A Figura 3.2 mostra a média móvel (com 100 épocas) da evolução da função de perda da PINN e do erro $L^{2}$ relativo da solução estimada em relação à solução analítica. Do critério de parada do treinamento, temos que a função de perda da PINN atingiu o valor de $10^{-3}$ , enquanto o erro é da ordem de $10^{-2}$ , uma ordem de grandeza maior. Verifique!

3.1.1 Exercícios

E. 3.1.1.

Considere o problema de Poisson no quadrado $D=(0,2)^{2}$

\displaystyle-\Delta u=-4,\leavevmode\nobreak\ (x,y)\in D,

(3.12)

com condições de contorno de Dirichlet não homogêneas nos lados $x=0$ e $x=2$ ,

\displaystyle u(0,y)=y^{2},\leavevmode\nobreak\ u(2,y)=4+y^{2},\leavevmode% \nobreak\ 0\leq y\leq 2,

(3.13)

e condições de contorno de Neumann³³3Carl Gottfried Neumann, 1832 - 1925, matemático alemão. Fonte: Wikipédia: Carl Neumann. nos lados $y=0$ e $y=2$ ,

\displaystyle u_{y}(x,0)=0,\leavevmode\nobreak\ u_{y}(x,2)=4,\leavevmode% \nobreak\ 0\leq x\leq 2.

(3.14)

Crie uma PINN $\tilde{u}=\mathcal{N}(x,y;\bm{\theta})$ para resolver este problema. Estude diferentes arquiteturas de MLP e funções de ativação. Compare a solução estimada com a solução analítica $u(x,y)=x^{2}+y^{2}$ .

E. 3.1.2.

Modifique o Código 17 para incluir pesos de penalidade na função de perda da PINN. Experimente diferentes valores de pesos e observe como eles afetam a convergência do treinamento e a precisão da solução estimada. Compare os resultados com a versão sem pesos de penalidade.

E. 3.1.3.

Considere o seguinte problema de Poisson

\displaystyle-\Delta u=f,\bm{x}\in[-1,1]^{2},

(3.15)

com fonte não suave

\displaystyle f(x_{1},x_{2})=\begin{cases}1,&\text{se }x_{1}^{2}+x_{2}^{2}<0.2% 5,\\ 0,&\text{caso contr\'{a}rio},\end{cases}

(3.16)

e condições de contorno de Dirichlet homogêneas. Crie uma PINN para resolver este problema e compare a solução estimada com a solução obtida por métodos numéricos tradicionais, como o método de diferenças finitas. Discuta as dificuldades encontradas na aproximação da solução devido à não suavidade da fonte.

E. 3.1.4.

Considere o problema estacionário de difusão-convecção-reação no domínio retangular $D=(0,1)\times(-2,1)$

\displaystyle-\Delta u+u_{x}+u_{y}+u=-e^{x-y},\leavevmode\nobreak\ (x,y)\in D,

(3.17)

e condições de contorno de Robin⁴⁴4Victor Gustave Robin, 1855 - 1897, matemático francês. Fonte: Wikipedia: Victor Gustave Robin.

\displaystyle u+\frac{\partial u}{\partial n}=g,\leavevmode\nobreak\ (x,y)\in% \partial D,

(3.18)

onde $\partial u/\partial n$ é a derivada normal exterior a $\partial D$ e

	$\displaystyle g(0,y)=0,\leavevmode\nobreak\ -2\leq y\leq 1,$		(3.19)
	$\displaystyle g(1,y)=2e^{1-y},\leavevmode\nobreak\ -2\leq y\leq 1,$		(3.20)
	$\displaystyle g(x,-2)=2e^{x+2},\leavevmode\nobreak\ 0\leq x\leq 1,$		(3.21)
	$\displaystyle g(x,1)=0,\leavevmode\nobreak\ 0\leq x\leq 1.$		(3.22)

Crie uma PINN $\tilde{u}=\mathcal{N}(x,y;\bm{\theta})$ para resolver este problema. Compare seus resultados com a solução analítica $u(x,y)=e^{x-y}$ .

E. 3.1.5.

Considere o problema de Poisson na coroa circular

\displaystyle D=\left\{(x,y)\in\mathbb{R}^{2}:1<x^{2}+y^{2}<4\right\},

(3.23)

dado por

\displaystyle-\Delta u=0,\leavevmode\nobreak\ (x,y)\in D,

(3.24)

com condições de contorno de Dirichlet não homogêneas

	$\displaystyle u(x,y)=0,\leavevmode\nobreak\ x^{2}+y^{2}=1,$		(3.25)
	$\displaystyle u(x,y)=\ln 2,\leavevmode\nobreak\ x^{2}+y^{2}=4.$		(3.26)

Crie uma PINN $\tilde{u}=\mathcal{N}(x,y;\bm{\theta})$ para resolver este problema, adaptando a estratégia de amostragem de pontos internos (por exemplo, por coordenadas polares ou rejeição) e de pontos de contorno para a geometria da coroa circular. Compare a solução estimada com a solução analítica

\displaystyle u(x,y)=\frac{1}{2}\ln\left(x^{2}+y^{2}\right).

(3.27)

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