I’m trying to optimize the parameters of this first state of a rc model using NeuralPDE but I get error. what am I doing wrong?
using NeuralPDE, Lux, ModelingToolkit, Optimization, OptimizationOptimJL, OrdinaryDiffEq,
Plots
import ModelingToolkit: Interval, infimum, supremum
using DifferentialEquations, Plots, Flux,Optim, DiffEqFlux, DataInterpolations,Random, ComponentArrays, Lux
using OptimizationOptimisers,OptimizationNLopt
rng = Random.default_rng()
using CSV
using DataFrames
using Statistics: mean
using OptimizationPolyalgorithms
using DiffEqFlux: group_ranges
#exogenous variables
df = CSV.read("/Users/marco/OneDrive/Desktop/Tesi/Julia_files/measured_data_ex_00.csv", DataFrame)
df=repeat(df, outer=10)
registered_gas_flow = df[:, :2]
Registered_Temperature = df[:, 3]
sec_Temp = [mean(Registered_Temperature[i:i+59]) for i in 1:60:length(Registered_Temperature)-59]
sec_gas = [mean(registered_gas_flow[i:i+59]) for i in 1:60:length(registered_gas_flow)-59]
tspan= (0.0f0,3000.0f0)
datasize = 3000
tsteps= range(tspan[1], tspan[2], length = datasize)
Ext_temperature = LinearInterpolation(sec_Temp,tsteps);
power = LinearInterpolation(sec_gas,tsteps);
function ext_Temp(tsteps)
return Ext_temperature(tsteps)
end
function ext_power(tsteps)
return power(tsteps)
end
#defining the problem
@parameters t, Rv, Ci
@variables x(..)
Dt = Differential(t)
eqs = [Dt(x(t)) ~ 1/(Rv*Ci) * (ext_Temp(t) - x(t)) + ext_power(t)/Ci]
bcs = [x(0) ~ 20.0]
domains = [t ∈ Interval(0.0, 3000.0)]
dt = 0.01
input_ = length(domains)
n = 8
chain1 = Lux.Chain(Lux.Dense(input_, n, Lux.σ), Lux.Dense(n, n, Lux.σ), Lux.Dense(n, n, Lux.σ),
Lux.Dense(n, 1))
# calculate the solution of the Lorenz system with OrdinaryDiffEq.jl based on the adaptivity of the ODE solver.
#This is used to introduce non-uniformity to the time series
function RC!(du,u,p,t)
du[1] = 1/(10000*10) * (ext_Temp(t)- u[1]) + ext_power(t)/10
end
u0 = [20.0]
tspan = (0.0, 3000.0)
prob = ODEProblem(RC!, u0, tspan)
sol = solve(prob, Tsit5(), dt = 0.1)
ts = [infimum(d.domain):dt:supremum(d.domain) for d in domains][1]
function getData(sol)
data = []
us = hcat(sol(ts).u...)
ts_ = hcat(sol(ts).t...)
return [us, ts_]
end
data = getData(sol)
(u_, t_) = data
len = length(data[2])
#Loss term based on the data that we have to optimize the parameters
depvars = [:x]
function additional_loss(phi, θ, p)
return sum(sum(abs2, phi[i](t_, θ[depvars[i]]) .- u_[[i], :]) / len for i in 1:1:1)
end
#A discretize algorithm for the ModelingToolkit PDESystem interface, which transforms a PDESystem into an OptimizationProblem using the Physics-Informed Neural Networks (PINN) methodology
discretization = NeuralPDE.PhysicsInformedNN([chain1],
NeuralPDE.GridTraining(dt), param_estim = true,
additional_loss = additional_loss)
@named pde_system = PDESystem(eqs, bcs, domains, [t], [x(t)], [ Rv, Ci ],
defaults = Dict([p .=> 1.0 for p in [ Rv, Ci]]))
prob = NeuralPDE.discretize(pde_system, discretization)
callback = function (p, l)
println("Current loss is: $l")
return false
end
res = Optimization.solve(prob, BFGS(); callback = callback, maxiters = 5000)
p_ = res.u[(end - 1):end]
minimizers = [res.u.depvar[depvars[i]] for i in 1:1]
ts = [infimum(d.domain):(dt / 10):supremum(d.domain) for d in domains][1]
u_predict = [[discretization.phi[i]([t], minimizers[i])[1] for t in ts] for i in 1:1]
plot(sol)
plot!(ts, u_predict, label = ["x(t)"])
println(p_)
say more?