Adaptive reweighting or modify neural network to solve complex system.

[HuynhTran0301]

I tried to use NeuralPDE to solve the ODE system. The results are not good when I used the code as a tutorial then I tried the adaptive_loss function and additional_loss function and got better results. My code is follow:

using NeuralPDE, Lux, ModelingToolkit, Optimization, OptimizationOptimJL, OptimizationOptimisers
import ModelingToolkit: Interval
using CSV
using DataFrames
data = CSV.File("3gens.csv");
Y1 = CSV.read("Y_init_change3.CSV", DataFrame, types=Complex{Float64});

#Input of the system.
E1 = 1.054;
E2 = 1.050;
E3 = 1.017;

#Define equation of the system
@variables delta1(..) delta2(..) delta3(..) omega1(..) omega2(..) omega3(..) TM1(..) TM2(..) TM3(..) Psv1(..) Psv2(..) Psv3(..)
@parameters t
D = Differential(t)
omega_s = round(120*pi,digits = 4)

eq1 = [
    D(delta1(t)) ~ omega1(t),
    D(delta2(t)) ~ omega2(t),
    D(delta3(t)) ~ omega3(t),
    D(omega1(t)) ~ (TM1(t)-(E1^2*Y1[1,1].re + E1*E2*sin(delta1(t)-delta2(t))*Y1[1,2].im + E1*E2*cos(delta1(t)-delta2(t))*Y1[1,2].re + E1*E3*sin(delta1(t)-delta3(t))*Y1[1,3].im + 
        E1*E3*cos(delta1(t)-delta3(t))*Y1[1,3].re))/(2*data["H"][1])*omega_s,
    #
    D(omega2(t)) ~ (TM2(t)-(E2^2*Y1[2,2].re + E1*E2*sin(delta2(t)-delta1(t))*Y1[2,1].im + E1*E2*cos(delta2(t)-delta1(t))*Y1[2,1].re + E2*E3*sin(delta2(t)-delta3(t))*Y1[2,3].im + 
        E2*E3*cos(delta2(t)-delta3(t))*Y1[2,3].re))/(2*data["H"][2])*omega_s,
    #
    D(omega3(t)) ~ (TM3(t)-(E3^2*Y1[3,3].re + E3*E1*sin(delta3(t)-delta1(t))*Y1[3,1].im + E3*E1*cos(delta3(t)-delta1(t))*Y1[3,1].re + E3*E2*sin(delta3(t)-delta2(t))*Y1[3,2].im + 
        E3*E2*cos(delta3(t)-delta2(t))*Y1[3,2].re))/(2*data["H"][3])*omega_s];
eq2 = [D(TM1(t)) ~ (-TM1(t)+Psv1(t))/data["TCH"][1],
    D(TM2(t)) ~ (-TM2(t)+Psv2(t))/data["TCH"][2],
    D(TM3(t)) ~ (-TM3(t)+Psv3(t))/data["TCH"][3]]

eq3 = [D(Psv1(t)) ~ (-Psv1(t) + 0.70945 + 0.335*(-0.0833) - (omega1(t)/omega_s)/data["RD"][1])/data["TSV"][1],
    D(Psv2(t)) ~ (-Psv2(t) + 1.62342 + 0.33*(-0.0833) - (omega2(t)/omega_s)/data["RD"][2])/data["TSV"][2],
    D(Psv3(t)) ~ (-Psv3(t) + 0.84843 + 0.335*(-0.0833) - (omega3(t)/omega_s)/data["RD"][3])/data["TSV"][3]];
eqs = [eq1;eq2;eq3];

bcs = [delta1(0.0) ~ 0.03957, delta2(0.0) ~ 0.3447, delta3(0.0) ~ 0.23038,
    omega1(0.0) ~ 0.0, omega2(0.0) ~ 0.0, omega3(0.0) ~ 0.0,
    TM1(0.0) ~ 0.70945, TM2(0.0) ~ 1.62342, TM3(0.0) ~ 0.848433,
    Psv1(0.0) ~ 0.70945, Psv2(0.0) ~ 1.62342, Psv3(0.0)~ 0.848433]

domains = [t ∈ Interval(0.0,25.0)]

chain =[Lux.Chain(Dense(1,10,Lux.tanh),Dense(10,20,Lux.tanh),Dense(20,10,Lux.tanh),Dense(10,1)) for _ in 1:12]

dvs = [delta1(t),delta2(t),delta3(t),omega1(t),omega2(t),omega3(t),TM1(t),TM2(t),TM3(t),Psv1(t),Psv2(t),Psv3(t)]

@named pde_system = PDESystem(eqs,bcs,domains,[t],dvs)

discretization = PhysicsInformedNN(chain, strategy, additional_loss = NeuralPDE.GradientScaleAdaptiveLoss(20), adaptive_loss  = NeuralPDE.GradientScaleAdaptiveLoss(20))

sym_prob = NeuralPDE.symbolic_discretize(pde_system, discretization)

pde_loss_functions = sym_prob.loss_functions.pde_loss_functions
bc_loss_functions = sym_prob.loss_functions.bc_loss_functions

callback = function (p, l)
    println("loss: ", l)
    return false
end
loss_functions =  [pde_loss_functions;bc_loss_functions]

function loss_function(θ,p)
    sum(map(l->l(θ) ,loss_functions))
end

f_ = OptimizationFunction(loss_function, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f_, sym_prob.flat_init_params);
phi = sym_prob.phi;

res = Optimization.solve(prob,OptimizationOptimJL.BFGS(); callback = callback, maxiters = 1500)

ts = 0.0:0.01:25.0
minimizers_ = [res.u.depvar[sym_prob.depvars[i]] for i in 1:12]
u_predict  = [[phi[i]([t],minimizers_[i])[1] for t in ts] for i in 1:12];

I got the results for the phase angle below:

However, when I solved the system by ODEPlroblem I got the results as follows:

We can see that the results from the two methods are approx in the first 5 seconds but after that, the neural network can not capture the nonlinear of the phase angle.

Please help me on how could I modify the neural network to capture the nonlinear functions for whole domains.
Thank you all.

[sdesai1287]

You could modify the strategy you use to train the neural network to better fit the problem. This can be done like this:

using NeuralPDE, OrdinaryDiffEq, OptimizationPolyalgorithms, Lux, OptimizationOptimJL, Test, Statistics, Plots, Optimisers

function f(u, p, t)
    [p[1] * u[1] - p[2] * u[1] * u[2], -p[3] * u[2] + p[4] * u[1] * u[2]]
end

p = [1.5, 1.0, 3.0, 1.0]
u0 = [1.0, 1.0]
prob_oop = ODEProblem{false}(f, u0, (0.0, 3.0), p)
true_sol = solve(prob_oop, Tsit5(), saveat = 0.01)
func = Lux.σ
N = 12
chain = Lux.Chain(Lux.Dense(1, N, func), Lux.Dense(N, N, func), Lux.Dense(N, N, func),
                    Lux.Dense(N, N, func), Lux.Dense(N, length(u0)))

opt = Optimisers.Adam(0.01)
weights = [0.7, 0.2, 0.1]
samples = 200
alg = NeuralPDE.NNODE(chain, opt, autodiff = false, strategy = NeuralPDE.WeightedIntervalTraining(weights, samples))
sol = solve(prob_oop, alg, verbose=true, maxiters = 100000, saveat = 0.01)

test(abs(mean(true_sol .- sol)))

using Plots

plot(sol)
plot!(true_sol)
ylims!(0,8)

This method allows you to select a weighting and number of samples for your problem. I found this to work for this specific problem, but you may have to experiment with some other weightings or samples for your problem to fit well. Please let me know if you have any questions about this and I will do my best to help you out

[HuynhTran0301]

I have tried with your code but I got the message UndefVarError: `WeightedIntervalTraining` not defined. Then, I modify your code as follow:

using NeuralPDE, OrdinaryDiffEq, OptimizationPolyalgorithms, Lux, OptimizationOptimJL, Test, Statistics, Plots, Optimisers, WeightedIntervalTraining

function f(u, p, t)
    [p[1] * u[1] - p[2] * u[1] * u[2], -p[3] * u[2] + p[4] * u[1] * u[2]]
end

p = [1.5, 1.0, 3.0, 1.0]
u0 = [1.0, 1.0]
prob_oop = ODEProblem{false}(f, u0, (0.0, 3.0), p)
true_sol = solve(prob_oop, Tsit5(), saveat = 0.01)
func = Lux.σ
N = 12
chain = Lux.Chain(Lux.Dense(1, N, func), Lux.Dense(N, N, func), Lux.Dense(N, N, func),
                    Lux.Dense(N, N, func), Lux.Dense(N, length(u0)))

opt = Optimisers.Adam(0.01)
weights = [0.7, 0.2, 0.1]
samples = 200
alg = NeuralPDE.NNODE(chain, opt, autodiff = false, strategy = NeuralPDE.WeightedIntervalTraining(weights, samples))
sol = solve(prob_oop, alg, verbose=true, maxiters = 100000, saveat = 0.01)

test(abs(mean(true_sol .- sol)))

using Plots

plot(sol)
plot!(true_sol)
ylims!(0,8)

And got another message:

ArgumentError: Package WeightedIntervalTraining not found in current path.
Run `import Pkg; Pkg.add("WeightedIntervalTraining")` to install the WeightedIntervalTraining package.

After using Pkg to add WeightedIntervalTraining package and got The following package names could not be resolved: WeightedIntervalTraining (not found in project, manifest or registry)
Please explant how could I use the WeightedIntervalTraining, because in the tutorial of NeuralPDE, I do not see the WeightedIntervalTraining strategy. So maybe it had been deleted?

[sdesai1287]

WeightedIntervalTraining is a new method of training that I developed. Can you verify that you have the latest version of NeuralPDE? There should be mentions of WeightedIntervalTraining in NeuralPDE.jl, ode_solve.jl and training_strategies.jl

[HuynhTran0301]

Oh, I forgot about that. Right now, I just use version 1.9.0. I will update the version and update the packages.

[HuynhTran0301]

I have updated the latest version, and your code is running well. For my code, I got a new error as follows:

 Can't differentiate foreigncall expression $(Expr(:foreigncall, :(Core.tuple("cos", NaNMath.libm)), Float64, svec(Float64), 0, :(:ccall), %4, %3)).
You might want to check the Zygote limitations documentation.

My equation is:

eq1 = [D(delta[1]) ~ omega[1],
    #
    D(delta[2]) ~ omega[2],
    #
    D(delta[3]) ~ omega[3],
    #
    D(omega[1]) ~ (TM[1]-Pe[1])/(2*data["H"][1])*omega_s,
    #
    D(omega[2]) ~ (TM[2]-Pe[2])/(2*data["H"][2])*omega_s,
    #
    D(omega[3]) ~ (TM[3]-Pe[3])/(2*data["H"][3])*omega_s,
    #
    0 ~ Pe[1] - (E1^2*G11 + E1*E2* sin(delta[1]-delta[2])*B12 + E1*E2* cos(delta[1]-delta[2])*G12 + E1*E3* sin(delta[1]-delta[3])*B13 + 
        E1*E3* cos(delta[1]-delta[3])*G13),
    #
    0 ~ Pe[2] - (E2^2*G22 + E1*E2* sin(delta[2]-delta[1])*B21 + E1*E2* cos(delta[2]-delta[1])*G21 + E2*E3* sin(delta[2]-delta[3])*B23 + 
        E2*E3* cos(delta[2]-delta[3])*G23),
    #
    0 ~ Pe[3] - (E3^2*G33 + E3*E1* sin(delta[3]-delta[1])*B31 + E3*E1* cos(delta[3]-delta[1])*G31 + E3*E2* sin(delta[3]-delta[2])*B32 + 
        E3*E2* cos(delta[3]-delta[2])*G32)
    ];

I used ODESystem to convert this equation to an ODEs system and then also used ODEProblem{false}() to create an out-of-place ODE problem.
If you need a full code, I will update it.

Besides, I would like to ask that if my initial conditions are a large number (200 or 300), the network can work with it. because the sigmoid function has a limited value of the input (like it will work well in intervals from 0 to 1). If not, is there any way to normalize the input data?

[xtalax]

Can you post a full error, with the code that causes it?

[HuynhTran0301]

The error had been solved with the latest version of Julia and packages.