Hi All!
I am trying to solve a simple ODE with the PINN technology.
There might be better ways to do that, but I am using it to test and understand few practical things.
This code worked perfectly as long as the problem was a nice toy problem, well scaled between 0 and 1 in therms of input and output.
I then changed the scale to see how to generalize the framework to different more realistic problems, and… the PINN simply does not work at all.
Any suggestion 
?
I am planning to use this knowledge for solving other 3D PDE problems, that might not be naturally nicely scaled.
I am trying to solve the differential equation for the equilibrium of a simply supported beam with uniform distributed load (could be self weight), Wikipedia reference.
I compare the results of the PINN with the analytic and finite difference results:
# Packages
using NeuralPDE, Lux, ModelingToolkit, Optimization, OptimizationOptimisers
import ModelingToolkit: Interval, infimum, supremum
using Plots
using MethodOfLines, NonlinearSolve
# Functions
analytic_sol_func(x, p, EI, L) = p * x * ( x^3 - 2*L*x^2 + L^3)/24/EI
# function runPINNbeam()
# Parameters
D = 100.0 # mm, pipe outer diameter
t = 20.0 # mm, pipe wall thickness
R1 = D/2-t # mm, internal radius
R2 = D/2 # mm, external radius
Ix = 1/4*π*(R2^4-R1^4) # mm^4, moment of inertia
A = π*(R2^2-R1^2) # mm^2, cross sectional area of the pipe
E = 210.0 # kN/mm2, elastcity modulus
pLoad = -0.01 # kN/mm
Length = 3000.0 # mm
EI = E*Ix
scaleCoeff = 1.0 # to bring up the NN to the right order of magnitude form a 0-1 range
@parameters x
@variables u(..)
Dxxxx = Differential(x)^4
Dxx = Differential(x)^2
# PDE
eq = EI * Dxxxx(u(x))*scaleCoeff ~ pLoad
# Boundary conditions
bcs = [u(0) ~ 0.0, u(Length) ~ 0.0,
Dxx(u(0)) ~ 0.0, Dxx(u(Length)) ~ 0.0]
# Space and time domains
domains = [x ∈ Interval(0.0, Length)]
# MethodOfLines, for FD solution
dx = 10.0
discretization = MOLFiniteDifference([x => dx])
@named pde_system = PDESystem(eq, bcs, domains, [x], [u(x)])
prob = discretize(pde_system, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson())
xsFD = sol[x]
u_MOL = sol[u(x)]
# Neural network
dim = 1 # number of dimensions
neurons = 10
chain = Lux.Chain(Lux.LayerNorm((dim,),identity), Dense(dim, neurons, Lux.tanh_fast), [Dense(neurons, neurons, Lux.tanh_fast) for _ =1:2], Dense(neurons, 1))
#strategy = QuadratureTraining(; batch = 200, abstol = 1e-6, reltol = 1e-6)
dx = 10.0
strategy = GridTraining(dx)
discretization = PhysicsInformedNN(chain, strategy)
@named pde_system = PDESystem(eq, bcs, domains, [x], [u(x)])
prob = discretize(pde_system, discretization)
lossHist = []
iter = 0
callback = function (p, l)
global iter += 1
println("Iter: $iter Current loss is: $l")
push!(lossHist, l)
return false
end
res = Optimization.solve(prob, ADAM(1e-1); callback = callback, maxiters = 100)
prob = remake(prob, u0 = res.minimizer)
res = Optimization.solve(prob, ADAM(1e-2); callback = callback, maxiters = 300)
prob = remake(prob, u0 = res.minimizer)
res = Optimization.solve(prob, ADAM(1e-5); callback = callback, maxiters = 500)
phi = discretization.phi
dx = 10
xs = [infimum(d.domain):(dx / 1):supremum(d.domain) for d in domains][1]
u_predict = [first(phi(x*ones(1,1), res.u)) for x in xs]
u_real = [analytic_sol_func(x, pLoad, EI, Length) for x in xs]
# Plots
x_plot = collect(xs)
plot(x_plot, u_real, label = "analytic", w=2, linestyle=:dash)
plot!(xsFD, scaleCoeff*u_MOL, label = "finite diff.", w=2)
plot!(x_plot, scaleCoeff*u_predict, label = "PINN", w=2)
xlabel!("x [mm]")
ylabel!("u [mm]")
savefig("beamPINN.png")
plot(lossHist, yscale=:log10, w=2, label="loss")
savefig("beamPINN_loss.png")
# return xs, u_real, u_predict
# end
# Main
# xs, u_real, u_predict = runPINNbeam()