Introdução a PINNs

3 Redes informadas pela física 3.1 Problemas estacionários 3.3 Problemas inversos

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

3.2 Problemas transientes

Como um problema fundamental em equações diferenciais parciais transientes, vamos estudar como podemos criar uma PINN para resolver o problema de calor

	$\displaystyle u_{t}=u_{xx}+f,(t,x)\in(0,1]\times(0,1),$		(3.28)
	$\displaystyle u(0,x)=g(x),x\in[0,1],$		(3.29)
	$\displaystyle u(t,-1)=u(t,1)=0,t\in(t_{0},tf].$		(3.30)

Como um exemplo, vamos assumir a fonte $f(t,x)=(\pi^{2}-1)e^{-t}\operatorname{sen}(\pi x)$ e a condição inicial $g(x)=\operatorname{sen}(\pi x)$ . Este problema foi manufaturado com a solução analítica

\displaystyle u(t,x)=e^{-t}\operatorname{sen}(\pi x).

(3.31)

3.2.1 PINN sem discretização temporal

A primeira abordagem que iremos adotar para resolver o problema de calor é criar uma PINN

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

(3.32)

com uma rede do tipo MLP, em que a entrada da rede são as coordenadas $(t,x)$ e a saída é a estimativa da solução do problema de calor (3.28)-(3.30).

Refer to caption — Figura 3.3: 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 18 é uma implementação da PINNs para esse resolver o problema de calor (verifique!). A função de perda da rede é baseada no resíduo da equação diferencial, na condição inicial e nas condições de contorno. O treinamento da rede é feito minimizando a função de perda, que é composta por

	$\displaystyle\varepsilon:=\underbrace{\frac{1}{n_{s,in}}\sum_{s=1}^{n_{s,in}}% \left\|f^{(s)}-\tilde{u}_{t}^{(s)}+\tilde{u}_{xx}^{(s)}\right\|^{2}}_{\text{res% \'{\i}duo}}$
	$\displaystyle\text{}\qquad+\underbrace{\frac{1}{n_{s,ci}}\sum_{s=1}^{n_{s,ci}}% \left\|\tilde{u}^{(s)}-g^{(s)}\right\|^{2}}_{\text{c.i.}}+\underbrace{\frac{1}{n% _{s,cc}}\sum_{s=1}^{n_{s,cc}}\left\|\tilde{u}^{s}-0\right\|^{2}}_{\text{c.c.}}.$		(3.33)

Aqui, o primeiro termo representa o resíduo da equação diferencial e é calculado em pontos de amostragem internos no domínio $(0,1]\times(0,1)$ . O segundo, representa a condição inicial e é calculado em pontos de amostragem no domínio com $t=0$ . O terceiro termo representa as condições de contorno e é calculado em pontos de amostragem no domínio com $x=0$ e $x=1$ . A cada época de treinamento, conjuntos randômicos de pontos de amostragem internos, de pontos de contorno e de pontos da condição inicial são gerados com distribuição uniforme.

A Figura 3.3 mostra uma comparação entre 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 18: pinn_calor.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))

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

25# treinamento

27## fonte

28def f(t, x):

29 return (torch.pi**2 - 1) * torch.exp(-t) * torch.sin(torch.pi*x)

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## tolerância

42tol = 1e-4

43persist = 5

45# contador de persistência

46n_conv = 0

48# validação

49t_v = torch.linspace(0, 1, 25, device=device)

50x_v = torch.linspace(0, 1, 50, device=device)

51T_v, X_v = torch.meshgrid(t_v, x_v, indexing='ij')

52S_v = torch.hstack((T_v.reshape(-1,1), X_v.reshape(-1,1)))

53S_v.requires_grad = True

55# laço de treinamento

56for epoch in range(nepochs):

58 model.train()

60 # amostras de pontos internos

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

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

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

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

65 S.requires_grad = True

67 # propagação nos pontos internos

68 U = model(S)

70 # derivadas de primeira ordem

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

72 grad_outputs=torch.ones_like(U),

73 create_graph=True)[0]

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

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

77 # derivadas de segunda ordem

78 # d2u/dx^2

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

80 grad_outputs=torch.ones_like(DUDX),

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

83 # rhs

84 L = DUDT - D2UDX

86 # fonte

87 F = f(S[:,0:1], S[:,1:2])

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

90 loss_res = torch.mean((F - L)**2)

92 # condição inicial

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

94 Uci = model(Sci)

95 Uexp_ci = torch.sin(torch.pi*x)

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

98 # condições de contorno

100 # contorno x = 0

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

102 Ucc1 = model(Scc1)

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

104

105 # contorno x = 1

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

107 Ucc2 = model(Scc2)

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

109

110 # função de perda total

111 loss = loss_res + loss_ci + loss_cc1 + loss_cc2

112

113 if (epoch % 100 == 0):

114 print(f'{epoch}: loss_res: {loss.item():.2e}, ci: {loss_ci.item():.2e}, cc1: {loss_cc1.item():.2e}, cc2: {loss_cc2.item():.2e}')

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

116 torch.save(track, f'track.pt')

117

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

119 n_conv += 1

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

121 if (n_conv >= persist):

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

123 break

124 else:

125 n_conv = 0

126

127 # backward

128 optim.zero_grad()

129 loss.backward()

130 optim.step()

3.2.2 PINN com discretização temporal

Uma abordagem alternativa para resolver o problema de calor é criar uma PINN que incorpora a discretização temporal. Neste caso, a rede neural tem apenas a variável espacial como entrada, e as saídas da rede são as estimativas da solução em diferentes instantes de tempo. A função de perda da rede é baseada no resíduo da equação diferencial, na condição inicial e nas condições de contorno, mas agora a equação diferencial é avaliada usando uma discretização temporal.

Como um exemplo, vamos considerar a discretização temporal pelo método de Euler implícito para o problema de calor (3.28)-(3.30). A equação diferencial é discretizada no tempo como

\displaystyle u^{n+1}-h_{t}u_{xx}^{(n+1)}=u^{(n)}+h_{t}f^{(n)},

(3.34)

onde $u^{(n)}$ é a solução no instante de tempo $t_{n}=nh_{t}$ , $f^{(n)}=f(t_{n},x)$ , $n=0,1,\dotsc,n_{t}-1$ , $n_{t}$ é o número de passos no tempo de tamanho $h_{t}=t_{f}/n_{t}$ .

O Código 19 é uma implementação da PINN com esta discretização temporal. A função de perda da rede é

	$\displaystyle\varepsilon:=\underbrace{\frac{1}{n_{s,in}(n_{t}-1)}\sum_{s=1}^{n% _{s,in}}\sum_{n=1}^{n_{t}-1}\left\|\tilde{u}^{(n),(s)}+h_{t}f^{(n),(s)}-\tilde{% u}^{(n+1),(s)}+h_{t}\tilde{u}_{xx}^{(n),(s)}\right\|^{2}}_{\text{res\'{\i}duo}}$
	$\displaystyle\text{}\qquad+\underbrace{\frac{1}{n_{s,ci}}\sum_{s=1}^{n_{s,ci}}% \left\|\tilde{u}^{(0),(s)}-g^{(s)}\right\|^{2}}_{\text{c.i.}}+\underbrace{\frac{% 1}{n_{s,cc}(nt-1)}\sum_{s=1}^{n_{s,cc}}\sum_{n=1}^{n_{t}-1}\left\|\tilde{u}^{(n% ),(s)}-0\right\|^{2}}_{\text{c.c.}}.$		(3.35)

Ao rodarmos o código, é notável que a convergência do treinamento é mais rápida do que na abordagem sem discretização temporal pelo código anterior (Código 18). Ainda, observamos que esta abordagem nos permite usar uma rede neural com menos camadas escondidas e menos neurônios por camada escondida, o que reduz o custo computacional. Verifique!

Código 19: pinn_calor_discretizado.py

⬇

1import torch

3# dispositivo GPU/CPU ?

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

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

7# discretização temporal

8t0 = 0.0 # tempo inicial

9tf = 1.0 # tempo final

10nt = 10 # num de passos no tempo

11ts = torch.linspace(t0, tf, nt+1, device=device) # instantes de tempo

12ht = (tf - t0)/nt # tamanho do passo no tempo

14# modelo

15n_h = 2 # num de camadas escondidas

16n_n = 25 # num de neurônios por camada escondida

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

19model = torch.nn.Sequential()

20# camada de entrada

21model.add_module('layer_1', torch.nn.Linear(1,n_n))

22model.add_module('fun_1', act_fun)

23# camadas escondidas

24for i in range(n_h-1):

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

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

27# camada de saída

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

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

31# treinamento

33## fonte

34def f(t, x):

35 return (torch.pi**2 - 1) * torch.exp(-t) * torch.sin(torch.pi*x)

37## optimizador

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

40## num de amostras por época

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

43## num max épocas

44nepochs = 50_000

46## tolerância

47tol = 1e-4

48persist = 5

50# contador de persistência

51n_conv = 0

53# laço de treinamento

54for epoch in range(nepochs):

56 model.train()

58 # amostras de pontos internos

59 Xs = torch.rand(ns_x, 1, device=device)

60 Xs.requires_grad = True

62 # propagação nos pontos internos

63 U = model(Xs)

65 loss_res = 0.0

66 for n in range(0, nt):

68 # derivadas de primeira ordem

69 DUDX = torch.autograd.grad(U[:,n+1:n+2], Xs,

70 grad_outputs=torch.ones_like(U[:,n+1:n+2]),

71 create_graph=True)[0]

73 # derivadas de segunda ordem

74 D2UDX = torch.autograd.grad(DUDX, Xs,

75 grad_outputs=torch.ones_like(DUDX),

76 create_graph=True)[0]

78 # resíduo da equação discretizada (Euler implícito)

79 R = U[:,n:n+1] + ht * f(ts[n+1], Xs) - U[:,n+1:n+2] + ht * D2UDX

81 loss_res += torch.sum(R**2)

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

84 loss_res = loss_res / (ns_x * nt)

86 # condição inicial

87 Uexp_ci = torch.sin(torch.pi*Xs)

88 loss_ci = torch.mean((U[:,0:1] - Uexp_ci)**2)

90 # condições de contorno

92 # contorno x = 0

93 Ucc1 = model(torch.zeros((1,1), device=device))

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

96 # contorno x = 1

97 Ucc2 = model(torch.ones((1,1), device=device))

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

100 # função de perda total

101 loss = loss_res + loss_ci + loss_cc1 + loss_cc2

102

103

104 if (epoch % 100 == 0):

105 print(f'{epoch}: loss_res: {loss.item():.2e}, ci: {loss_ci.item():.2e}, cc1: {loss_cc1.item():.2e}, cc2: {loss_cc2.item():.2e}')

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

107

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

109 n_conv += 1

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

111 if (n_conv >= persist):

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

113 break

114 else:

115 n_conv = 0

116

117 # backward

118 optim.zero_grad()

119 loss.backward()

120 optim.step()

3.2.3 Exercícios

E. 3.2.1.

Considere o problema de calor com condições de contorno de Robin⁵⁵5Victor Gustave Robin, 1855 - 1897, matemático francês. Fonte: Wikipedia: Victor Gustave Robin. não homogêneas

	$\displaystyle u_{t}=u_{xx}+f(t,x),\leavevmode\nobreak\ (t,x)\in(0,1)\times(0,1),$		(3.36)
	$\displaystyle u(0,x)=2x+1,\leavevmode\nobreak\ x\in[0,1],$		(3.37)
	$\displaystyle u+\frac{\partial u}{\partial n}=g,\leavevmode\nobreak\ x\in\{0,1% \},\leavevmode\nobreak\ t\in(0,1],$		(3.38)

onde $\partial u/\partial n$ é a derivada normal exterior ao domínio $(0,1)$ em $x=0$ e $x=1$ , com fonte

\displaystyle f(t,x)=-(2x+1)e^{-t},

(3.39)

e condições de contorno com

\displaystyle g(t,0)=-e^{-t},\leavevmode\nobreak\ g(t,1)=5e^{-t},\leavevmode% \nobreak\ t\in(0,1].

(3.40)

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

E. 3.2.2.

Refaça o E.3.2.1 para criar PINNs da forma

\displaystyle\tilde{\bm{u}}=\mathcal{N}(x;\bm{\theta}),

(3.41)

tendo como entrada apenas a coordenada espacial $x$ e como saída as estimativas da solução em diferentes instantes de tempo. Faça uma PINN para cada uma das seguintes discretizações temporais:

a)

Euler implícito.
b)

Euler explícito.
c)

Crank-Nicolson.
d)

Runge-Kutta de quarta ordem.

Quais dos esquemas de discretização temporal você considera mais eficiente para criar a PINN para este problema? Justifique sua resposta.

E. 3.2.3.

Considere o problema de onda

	$\displaystyle u_{tt}-u_{xx}=0,\leavevmode\nobreak\ (t,x)\in(0,2)\times(0,1),$		(3.42)
	$\displaystyle u(0,x)=0,\leavevmode\nobreak\ u_{t}(0,x)=\pi\operatorname{sen}(% \pi x),\leavevmode\nobreak\ x\in[0,1],$		(3.43)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ t\in(0,2].$		(3.44)

Este problema tem solução analítica $u(t,x)=\operatorname{sen}(\pi x)\sin(\pi t)$ .

a)

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

Considerando uma discretização temporal por diferenças finitas centrais de ordem 2, crie uma PINN $\tilde{u}=\mathcal{N}(x;\bm{\theta})$ tendo como entrada apenas a coordenada espacial $x$ e como saída as estimativas da solução em diferentes instantes de tempo. Compare seus resultados com a solução analítica.
c)

Qual as duas abordagens apresentadas para resolver o problema de onda você considera mais eficiente? Justifique sua resposta.

E. 3.2.4.

Crie uma PINN para resolver o seguinte problema de Burgers

	$\displaystyle u_{t}+uu_{x}-\nu u_{xx}=0,\leavevmode\nobreak\ (t,x)\in(0,1]% \times(0,1),$		(3.45)
	$\displaystyle u(0,x)=2\nu\pi\frac{\operatorname{sen}(\pi x)}{2+\cos(\pi x)},% \leavevmode\nobreak\ x\in[0,1],$		(3.46)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ t\in(0,1],$		(3.47)

com viscosidade $\nu=0.01/\pi$ . Este problema tem a seguinte solução analítica [19]

\displaystyle u(t,x)=2\nu\pi\frac{e^{-\pi^{2}\nu t}\operatorname{sen}(\pi x)}{% 2+e^{-\pi^{2}\nu t}\cos(\pi x)}.

(3.48)

E. 3.2.5.