I am Trying to model some temperature curves with a physically informed NN (PINN).
I started solve the FOPDT first, to be able to apply it later.
function FOPDT(z)
τ, Κ, θ = z
return [T1[i] .- (1 ./ τ) .* (-T1[i] .+ Κ .* 100 .*(i .- θ)) for i in [1, 75, 130]]
end
initial_x = [1.0, 1.0, 1.0]
sol = nlsolve(FOPDT, initial_x, xtol=1e-2, ftol=1e-2, iterations=1_000_000)
sol
The temperatures at the time steps are:
0 s → 21.5 °C
75 s → 46.0 °C
130 s → 60.3 °C
But I was not able to find a static solution. Hence I used Scipy.optimize.fsolve. I tried to use
this solution as starting parameters but this also did not converge. Are there some other hints, how to solve the equation? Did I do any other mistake, e.g. in function?
T1 isn’t defined in your example.
using NonlinearSolve
function FOPDT(dz,z,p)
τ, Κ, θ = z
return dz .= [T1[i] .- (1 ./ τ) .* (-T1[i] .+ Κ .* 100 .*(i .- θ)) for i in [1, 75, 130]]
end
initial_x = [1.0, 1.0, 1.0]
prob = NonlinearProblem(FOPDT, initial_x)
sol = solve(prob, NonlinearSolve.NewtonRaphson(), abstol=1e-2)
Can you please give enough information to run this?
Hi,
I added the temperatures from the text.
function FOPDT(z)
τ, Κ, θ = z
return [T - (1 / τ) * (-T + Κ * 100 *(i - θ)) for (T, i) in zip([21.543, 46.035,60.279], [1, 75, 130])]
end
initial_x = [145.7, 66.3, 73.4] # from scipy.optimize.fsolve
sol = nlsolve(FOPDT, initial_x, xtol=1e-1, ftol=1e-2)
The result is:
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [145.7, 66.3, 73.4]
* Zero: [-7.290966107926511e11, -2.2018393026492186e9, -72.45793039799061]
* Inf-norm of residuals: 1.503362
* Iterations: 183
* Convergence: true
* |x - x'| < 1.0e-01: true
* |f(x)| < 1.0e-02: false
* Function Calls (f): 184
* Jacobian Calls (df/dx): 100
That point very much evaluates to non-zero. But can you share the SciPy code? It would take the comparison in order to see what’s different in the function.
Here is the scipy variant:
from scipy.optimize import fsolve
def f(z):
tau, kappa, theta = z
return [T-1/tau*(-T + kappa*100*(i - theta)) for T, i in zip([21.543, 46.035, 60.279], [0, 75, 130])]
tau, kappa, theta = fsolve(f, [1, 1, 1], xtol=1e-01, maxfev=100_000)
A solution was only found with adding the xtol parameter and setting it to: 1e-1. I don’t know how much to trust the result.
Well that’s easy to check:
>>> def f(z):
... tau, kappa, theta = z
... return [T-1/tau*(-T + kappa*100*(i - theta)) for T, i in zip([21.543, 46.035, 60.279], [0, 75, 130])]
...
>>> tau, kappa, theta = fsolve(f, [1, 1, 1], xtol=1e-01, maxfev=100_000)
>>> f([tau, kappa, theta])
[3360.9606761659174, -27.294401084985765, -2515.757097939305]
It’s also not a zero in the Python code.
Strange. But I assume I did a mistake in the formula.
In addition, I solved the solution for cooling phase and than for the heating phase.
Cooling:
def f_off(z):
_tau = z[0]
def k(c: float):
return c + 273.15
return [dT + k(T)/_tau for dT, T in zip([-0.2425,], [76.006,] )]
r = fsolve(f_off, [1.0,], xtol=1e-6, maxfev=1000_000)
Heating:
def f_on(z):
_tau = 1439.818557
_kappa, _theta = z
def k(c: float):
return c + 273.15
return [dT + 1/_tau*(T + _kappa*100*(_theta -i)) for dT, T, i in zip([0.279, 0.3065], [46.035, 60.279], [75, 130])]
kappa, theta = fsolve(f_on, [66, 73], xtol=5e-00, maxfev=1000_000)
I got:
tau = 1439.818557,
kapp = 0.0097885
theta = 72.93
Again no solution found in the Julia. Should I try another library? Which?
I did a bug fix and compared the python and the julia solution.
Both are stable and results are the same with:
function FOPDT_off(z)
τ = z[1]
k(c) = c + 273.15
return [dT + k(T) / τ for (dT, T) in zip([-0.2425,], [76.006,])]
end
function FOPDT_on(z)
τ = 1439.8
Κ, θ = z
k(c) = c + 273.15
return [dT + (1 / τ) * (k(T) + Κ * 100 *(θ -i)) for (dT, T, i) in zip(
[0.279, 0.3065],
[46.035,60.279],
[75, 130])]
end
The solutions are: τ=1439.8, Κ=0.009788, θ=-661.4 .
1 Like