Parameter Estimation of a 3 node thermal network

Im trying to formulate a parameter estimation problem to estimate the Resistance and Capacitance values of a 3 node thermal network. The thermal network is in the form of DAE with a mass matrix and exogenous inputs. I already have the temperature values of my 3 nodes as an interpolation function and im using it to estimate the loss of the then use Optimization.solve to solve the problem.

using DifferentialEquations, Interpolations,
Optimization, OptimizationPolyalgorithms, SciMLSensitivity,
Zygote, OptimizationOptimisers, OptimizationNLopt
using Plots
import XLSX

function thermal(du, u, p, t, P1, P2)
T1, T2, T3 = u
R_12, R_32, R_35, R_24, T4, T5 = p
du[1] = T2/R_12 - T1/R_12 + P1(t);
du[2] = T2/R_12 + T4/R_24 + T3/R_32 - T2*(1/R_12+1/R_24+1/R_32);
du[3] = T2/R_32 + T5/R_35 - T3*(1/R_32+1/R_35) + P2(t);
end

p=[0.5,0.2,0.1,0.05,30,30,1000,3000,2000];
M = [p[end-2] 0 0
0 p[end-1] 0
0 0 p[end]]

sh = xf[â€śExp1â€ť]
p1 = sh[â€śC5:C1804â€ť]
p2 = sh[â€śD5:D1804â€ť]
t1 = sh[â€śH5:H180004â€ť]
t2 = sh[â€śI5:I180004â€ť]
t3 = sh[â€śJ5:J180004â€ť]
d = 0.01:0.01:1800;
t = 1:1:1800;
tsteps = 0.1:0.1:1800
P1=LinearInterpolation(t,vec(p1â€™));
P2=LinearInterpolation(t,vec(p2â€™));
T1_s=LinearInterpolation(d,vec(t1â€™));
T2_s=LinearInterpolation(d,vec(t2â€™));
T3_s=LinearInterpolation(d,vec(t3â€™));

#experimental power and temperature values

u=[20,20,20];
tspan = (1,1800);
f = ODEFunction(thermal, mass_matrix = M)
prob_mm = ODEProblem(f, u, tspan, p)
sol = solve(prob_mm, Rodas5P(), saveat = 0.01, abstol=1e-8, reltol=1e-8)
plot(sol);

function loss(p)
sol = solve(prob_mm, Rodas5P(), p=p, saveat = 0.01, abstol=1e-8, reltol=1e-8)
T1_hat = [u[1] for u in sol.u]
T2_hat = [u[2] for u in sol.u]
T3_hat = [u[3] for u in sol.u]

``````T1_model = T1_s(sol.t)
T2_model = T2_s(sol.t)
T3_model = T3_s(sol.t)

loss1 = sum(abs2, T1_model .- T1_hat)
loss2 = sum(abs2, T2_model .- T2_hat)
loss3 = sum(abs2, T3_model .- T3_hat)

loss = loss1+loss2+loss3
return loss, sol
``````

end

callback = function (p, l, pred)
display(l)

display(plt)

``````# Tell Optimization.solve to not halt the optimization. If return true, then
# optimization stops.
return false
``````

end

optf = Optimization.OptimizationFunction((x, p) â†’ loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, p)

result_ode = Optimization.solve(optprob, PolyOpt(), callback=callback, maxiters = 100)

I have tried with different AD solvers and different algorithms but this keeps throwing me
â€śERROR: gradient of 180000-element extrapolate(scale(interpolate(::Vector{Any}, BSpline(Linear())), (0.01:0.01:1800.0,)), Throw()) with element type Float64:â€ť
3node_thermal_network.jl (2.2 KB)

Hello and welcome to the community Could you post the data file, such that people can run your code and reproduce the error you are seeing?

Hi Fredrik,
The data file i have in a .xlsx file and i am not able to upload it, is there any other way i could possibly share the file ?

Any standard file sharing service I guess, google drive, dropbox etc.?

Hereâ€™s an example that uses a simple implementation of the prediction-error method. I call the function `nonlinear_pem`, even though this system is actually linear in the state and inputs.

``````function thermal(u, input, p, t)
P1, P2 = input
T1, T2, T3 = u
R_12, R_32, R_35, R_24, T4, T5 = p
SA[(T2/R_12 - T1/R_12 + P1) / p[end-2]
(T2/R_12 + T4/R_24 + T3/R_32 - T2*(1/R_12+1/R_24+1/R_32)) / p[end-1]
(T2/R_32 + T5/R_35 - T3*(1/R_32+1/R_35) + P2) / p[end]]
end

sh = xf["Exp1"]
p1 = Float64.(sh["C5:C1804"][:])
p2 = Float64.(sh["D5:D1804"][:])
t1 = Float64.(sh["H5:H180004"][:])
t2 = Float64.(sh["I5:I180004"][:])
t3 = Float64.(sh["J5:J180004"][:])
d = 0.01:0.01:1800;
t = 1:1:1800;
plot(t, u, layout=(5,1))
plot!(d, y, sp=(3:5)')

Ts = 1
t = 1:Ts:1800

## experimental power and temperature values
u = [p1 p2] # Replace with input  measurement data here
y = [t1 t2 t3] # Replace with output measurement data here

d = iddata(y[1:100:end, :]', u', Ts) # Use every 100th sample of the outputs

# Initial guess
p0 = [0.5,0.2,0.1,0.05,30,30,1000,3000,2000]
u0 = SA[20.0,20,20]

R1 = 1e-4I(3) # Tuning parameters for the internal Kalman filter
R2 = 0.1I(3)
nx, nu, ny = 3, 2, 3 # Number of state variables, inputs and outputs

discrete_dynamics = SeeToDee.Rk4(thermal, Ts) # Discretize the continuous-time dynamics
measurement = (u, input, p, t) -> u # Measurement function is the identity

model = ControlSystemIdentification.nonlinear_pem(
d,
discrete_dynamics,
measurement,
p0,
u0,
R1,
R2,
nu;
optimize_x0 = false,
lower = 0.1*p0, # Constrain the search space to be within 10% of the initial guess
upper = 10*p0,
)

simplot(model, d, layout=(3, 1))
``````

The estimated parameters are available as

``````julia> model.p
9-element Vector{Float64}:
0.5042616495643081
0.22669213158156562
0.10462125464969055
0.447811612771852
300.0
29.79321313071231
1006.2410623845433
30000.0
2004.2052737272784
``````

The model fit is not perfect, in particular, there seem to be some problem with modeling the second output. This can be either due to missing dynamics in the model, due to the gradient-based optimizer getting stuck in a local minimum, or due to me constraining the parameters to be within 0.1-10x of the initial guess. You probably have better insight into what values are reasonable here

If the problem is due to a local minimum, you could try some form of global optimizer instead. However, this is not supported by the function `nonlinear_pem`. @SebastianM-C might have suggestions that can handle local minima better.

1 Like

BTW, a black-box linear model of order 4 fits the data perfectly:

``````bbmodel, x0 = newpem(d, 4, zeroD=true, focus=:simulation)
simplot(bbmodel, dd, x0, layout = (3, 1))
``````

Is the data simulated?

1 Like

This is really great. Actually i have been using Matlabâ€™s global optimization toolbox to estimate the values and been getting decent results but i havenâ€™t tried implementing the same in Julia yet.

Hi
i have kind of modified my code to embed neural networks into my differential equation and it trains well for my training data (which is first 200sec of my simulation) but as soon as i try to extrapolate the prediction to a higher time the results donâ€™t match. Im providing my code below, please take a look at it

using DifferentialEquations, Interpolations,Optimization, OptimizationPolyalgorithms, SciMLSensitivity,
Zygote, OptimizationOptimisers, OptimizationNLopt;
using Plots;
using OptimizationOptimJL;
using LinearAlgebra, Statistics;
using RecursiveArrayTools,CommonSolve;
using Lux, ComponentArrays, Zygote, StableRNGs;
import XLSX

gr()
rng = StableRNG(1111)

function thermal(du, u, p_, t)
T1, T2, T3 = u
R_12, R_32, R_35, R_24, T4, T5 = p_
du[1] = T2/R_12 - T1/R_12 + P1(t);
du[2] = T1/R_12 + T4/R_24 + T3/R_32 - T2*(1/R_12+1/R_24+1/R_32);
du[3] = T2/R_32 + T5/R_35 - T3*(1/R_32+1/R_35) + P2(t);
end
p_=[0.5,0.2,0.1,0.05,30,30,1000,3000,2000];
M = [p_[end-2] 0 0
0 p_[end-1] 0
0 0 p_[end]]

sh = xf[â€śExp1â€ť]
p1 = sh[â€śC5:C1804â€ť]
p2 = sh[â€śD5:D1804â€ť]
t1 = sh[â€śH5:H180004â€ť]
t2 = sh[â€śI5:I180004â€ť]
t3 = sh[â€śJ5:J180004â€ť]
d = 0.01:0.01:1800;
t = 1:1:1800;
tsteps = 0.1:0.1:1800
P1=Interpolations.LinearInterpolation(t,vec(p1â€™));
P2=Interpolations.LinearInterpolation(t,vec(p2â€™));

T1_s=Interpolations.LinearInterpolation(d,vec(t1â€™));
T2_s=Interpolations.LinearInterpolation(d,vec(t2â€™));
T3_s=Interpolations.LinearInterpolation(d,vec(t3â€™));

u=[20,20,20];
tspan = (1,200);
f = ODEFunction(thermal, mass_matrix = M)
prob_mm = ODEProblem(f, u, tspan, p_)
sol = solve(prob_mm, Rodas5P(), saveat = 0.1, abstol=1e-8, reltol=1e-8)
X = Array(sol);
t = sol.t;
plot(t, transpose(X));

rbf(x) = exp.(-(x .^ 2))
const U = Lux.Chain(Lux.Dense(2, 10, rbf), Lux.Dense(10, 10, rbf), Lux.Dense(10, 10, rbf),
Lux.Dense(10, 1))
p, st1 = Lux.setup(rng, U)
const _st1 = st1

const V = Lux.Chain(Lux.Dense(2, 10, rbf), Lux.Dense(10, 10, rbf), Lux.Dense(10, 10, rbf),
Lux.Dense(10, 1))
q, st2 = Lux.setup(rng, V)
const _st2 = st2

p = ComponentArray(p)
q = ComponentArray(q)
r = ComponentArray{Float64}()
r = ComponentArray(r;p)
r = ComponentArray(r;q)

p_=[0.5,0.2,0.1,0.05,30,30,1000,3000,2000];
p_ = ComponentArray([p_[1], p_[2], p_[3], p_[4], p_[5], p_[6], p_[7], p_[8], p_[9]])
M = [p_[end-2] 0 0
0 p_[end-1] 0
0 0 p_[end]]

Define the hybrid model

function ude_dynamics!(du, u, r, t)
uĚ‚ = U(u[1:2], r.p, _st1)[1] # Network prediction
g = V(u[2:3], r.q, _st2)[1]
p_true=[0.5,0.2,0.1,0.05,30,30,1000,3000,2000];
R_12, R_32, R_35, R_24, T4, T5 = p_true
du[1] = uĚ‚[1] + P1(t) ;
du[2] = -uĚ‚[1] + (T4-u[2])/R_24 + g[1];
du[3] = -g[1] + (T5-u[3])/R_35 + P2(t);

end

ude_dynamics = ODEFunction(ude_dynamics!, mass_matrix = M)
#nn_dynamics(du, u, p, t) = ude_dynamics(du, u, p, t)
prob_nn = ODEProblem(ude_dynamics, X[:, 1], tspan, r)

function predict(Î¸, x = X[:, 1], T = t)

``````#_prob = remake(prob_nn, u0 = x, tspan = (T[1], T[end]), p=Î¸)
Array(solve(prob_nn, Rodas5P(), p=Î¸, saveat = 0.1,
abstol = 1e-8, reltol = 1e-8,
``````

end

function loss(Î¸)
XĚ‚ = predict(Î¸)
mean(abs2, X .- XĚ‚)
end

losses = Float64

callback = function (r, l)
push!(losses, l)
if length(losses) % 50 == 0
println(â€śCurrent loss after (length(losses)) iterations: (losses[end])â€ť)
end
return false
end

optf = Optimization.OptimizationFunction((x, r) â†’ loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf,r)

res1 = Optimization.solve(optprob, ADAM(), callback = callback, maxiters = 1000)
println(â€śTraining loss after (length(losses)) iterations: (losses[end])â€ť)

optprob2 = Optimization.OptimizationProblem(optf, res1.u)
res2 = Optimization.solve(optprob2, Optim.LBFGS(), callback = callback, maxiters = 5000)
println(â€śFinal training loss after (length(losses)) iterations: (losses[end])â€ť)

Rename the best candidate

p_trained = res2.u

pl_losses = plot(1:1000, losses[1:1000], yaxis = :log10, xaxis = :log10,
xlabel = â€śIterationsâ€ť, ylabel = â€śLossâ€ť, label = â€śADAMâ€ť, color = :blue)
plot!(1001:length(losses), losses[1001:end], yaxis = :log10, xaxis = :log10,
xlabel = â€śIterationsâ€ť, ylabel = â€śLossâ€ť, label = â€śBFGSâ€ť, color = :red)

#ts = first(sol.t):(mean(diff(sol.t)) / 2):last(sol.t)

XĚ‚ = predict(p_trained, X[:, 1], t)

Trained on noisy data vs real solution

pl_trajectory = plot(t, transpose(XĚ‚), xlabel = â€śtâ€ť, ylabel = â€śT1(t), T2(t),T3(t)â€ť,
label = [â€śUDE Approximationâ€ť nothing])

I have built upon the example that was provided in documentation Automatically Discover Missing Physics by Embedding Machine Learning into Differential Equations Â· Overview of Julia's SciML

i have a different variation of code where i tried to learn one part of each of my 3 equations using one neural network but that too didnâ€™t work.

rbf(x) = exp.(-(x .^ 2))
const U = Lux.Chain(Lux.Dense(3, 5, rbf), Lux.Dense(5, 5, rbf), Lux.Dense(5, 5, rbf),
Lux.Dense(5, 3))
p_, st = Lux.setup(rng, U)
const _st = st

function ude_dynamics!(du, u, p_, t, p_true)
uĚ‚ = U(u, p_, _st)[1] # Network prediction
T1, T2, T3 = u
R_12, R_32, R_35, R_24, T4, T5 = p_true
du[1] = T2/R_12 - uĚ‚[1] + P1(t);
du[2] = T2/R_12 + T4/R_24 + T3/R_32 - uĚ‚[2];
du[3] = T2/R_32 + T5/R_35 - uĚ‚[3] + P2(t);

end

Im not able to understand why my implementation is not working, is it because i have a mass matrix and exogenous inputs involved ? if so how can i embed neural networks into ODE problems with both exogenous inputs and a mass matrix so that i can learn certain parts of my differential equation ??