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]]

xf = XLSX.readxlsx(â€śTempdat.xlsxâ€ť)

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)## plt = plot(pred)

## display(plt)

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

end

#adtype = Optimization.AutoZygote()

#adtype = Optimization.AutoReverseDiff()

adtype = Optimization.AutoForwardDiff()

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)