heatingsystem.jl

using OrdinaryDiffEq
using Sundials
using Plots
using BenchmarkTools

const N = 40  # Number of heated units

const Cu = N == 1 ? [2e7] : (ones(N) .+ range(0,1.348,length=N))*1e7 # "Heat capacity of heated units";
const Cd = 2e6*N # "Heat capacity of distribution circuit";
const Gh = 200 # "Thermal conductance of heating elements";
const Gu = 150 # "Thermal conductance of heated units to the atmosphere";
const Qmax = N*3000 # "Maximum power output of heat generation unit";
const Teps = 0.5 # "Threshold of heated unit temperature controllers";
const Td0 = 343.15 # "Set point of distribution circuit temperature";
const Tu0 = 293.15 # "Heated unit temperature set point";
const Kp = Qmax/4 # "Proportional gain of heat generation unit temperature controller";
const a = 50 # "Gain of the histeresis function";
const b = 15 # "Slope of the saturation function at the origin";

function set_up_params(p)

    Cu   = @view p[1:N]
    Cd   = p[N+1]
    Gh   = p[N+2]
    Gu   = p[N+3]
    Qmax = p[N+4]
    Teps = p[N+5]
    Td0  = p[N+6]
    Tu0  = p[N+7]
    Kp   = p[N+8]
    a    = p[N+8]
    b    = p[N+9]

    Cu, Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b
end

function set_up_vars(u)

    Td = u[1]
    Tu = @view u[2:N+1]
    x  = @view u[N+2:2N+1]

    Td, Tu, x
end

function combine_vars!(du, dTd, dTu, dx)

    du[1] = dTd
    du[2:N+1] .= dTu[1:N]
    du[N+2:2N+1] .=  dx[1:N]

    nothing
end

hist(x, p, e = 1) = -(x + 0.5)*(x - 0.5) *
                        x * 1/(0.0474)*e + p

sat(x, xmin, xmax) = tanh(2*(x-xmin)/(xmax-xmin)-1)*
                          (xmax-xmin)/2 +
                          (xmax+xmin)/2


function heating(du, u, (p, Qh, Que), t)

    (Td, Tu, x) = set_up_vars(u)
    (dTd, dTu, dx) = set_up_vars(du)
    (Cu, Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b) = set_up_params(p)

    Text = 278.15 + 8*sin(2π*t/86400)
    Qh .= Gh .* (Td .- Tu) .* (sat.(b .* x, -0.5, 0.5) .+ 0.5)
    Que .= Gu .* (Tu .- Text)

    Qd = sat(Kp*(Td0 - Td),0, Qmax)
    dTd = (Qd - sum(Qh))/Cd
    dTu .= (Qh .- Que) ./ Cu
    dx  .=  a .* hist.(x, Tu0 .- Tu, Teps)

    combine_vars!(du, dTd, dTu, dx)

end
u0 = [Td0, fill(Tu0, N)..., fill(-0.5, N)...]
Qh0 = zeros(N)
Que0 = zeros(N)
p = [Cu..., Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b]
tspan = (0., 50_000.)

prob = ODEProblem(heating, u0, tspan, (p, Qh0, Que0))

GC.gc()

@btime sol = solve(prob, RadauIIA5(autodiff=false), rtol = 1e-7)

length(sol.t)

plot(sol, vars=1)

heatingsystemnew.jl

using OrdinaryDiffEq, LabelledArrays

N = 40  # Number of heated units

Cu = N == 1 ? [2e7] : (ones(N) .+ range(0,1.348,length=N))*1e7 # "Heat capacity of heated units";
Cd = 2e6*N # "Heat capacity of distribution circuit";
Gh = 200 # "Thermal conductance of heating elements";
Gu = 150 # "Thermal conductance of heated units to the atmosphere";
Qmax = N*3000 # "Maximum power output of heat generation unit";
Teps = 0.5 # "Threshold of heated unit temperature controllers";
Td0 = 343.15 # "Set point of distribution circuit temperature";
Tu0 = 293.15 # "Heated unit temperature set point";
Kp = Qmax/4 # "Proportional gain of heat generation unit temperature controller";
a = 50 # "Gain of the histeresis function";
b = 15 # "Slope of the saturation function at the origin";

hist(x, p, e = 1) = -(x + 0.5)*(x - 0.5) *
                        x * 1/(0.0474)*e + p

sat(x, xmin, xmax) = tanh(2*(x-xmin)/(xmax-xmin)-1)*
                          (xmax-xmin)/2 +
                          (xmax+xmin)/2

function heating(du, u, (Cu, Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b, N), t)

    @views Td, Tu, x = u[1], u[2:N+1], u[N+2:2N+1]
    @views dTu, dx = du[2:N+1], du[N+2:2N+1]

    text = 278.15 + 8*sin(2π*t/86400)
    Qh = Gh .* (Td .- Tu) .* (sat.(b .* x, -0.5, 0.5) .+ 0.5)
    Que = Gu .* (Tu .- text)

    Qd = sat(Kp*(Td0 - Td),0, Qmax)
    du[1] = (Qd - sum(Qh))/Cd # dTd
    dTu .= (Qh .- Que) ./ Cu
    dx  .=  a .* hist.(x, Tu0 .- Tu, Teps)
    return nothing
end
tspan = (0., 50_000.)

p = (Cu=Cu, Cd=Cd, Gh=Gh, Gu=Gu, Qmax=Qmax, Teps=Teps, Td0=Td0, Tu0=Tu0, Kp=Kp, a=a, b=b, N=N)
u0 = [Td0, fill(Tu0, N)..., fill(-0.5, N)...]

prob = ODEProblem(heating, u0, tspan, p);
solve(prob, KenCarp4(), rtol=1e-7);
@time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
@time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
@time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
solve(prob, RadauIIA5(), rtol=1e-7);
@time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
@time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
@time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
solve(prob, Rosenbrock23(), rtol=1e-7);
@time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
@time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
@time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);

#=
julia> solve(prob, KenCarp4(), rtol=1e-7);

julia> @time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
  4.648946 seconds (4.05 M allocations: 1.762 GiB, 3.83% gc time)

julia> @time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
  4.703236 seconds (4.05 M allocations: 1.762 GiB, 3.71% gc time)

julia> @time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
  4.552906 seconds (4.05 M allocations: 1.762 GiB, 3.56% gc time)

julia> solve(prob, RadauIIA5(), rtol=1e-7);

julia> @time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
  5.216250 seconds (2.69 M allocations: 878.738 MiB, 5.67% gc time)

julia> @time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
  4.816831 seconds (2.69 M allocations: 878.738 MiB, 1.91% gc time)

julia> @time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
  4.816324 seconds (2.69 M allocations: 878.738 MiB, 1.73% gc time)

julia> solve(prob, Rosenbrock23(), rtol=1e-7);

julia> @time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
  1.896335 seconds (1.58 M allocations: 938.699 MiB, 5.20% gc time)

julia> @time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
  1.895587 seconds (1.58 M allocations: 938.699 MiB, 4.92% gc time)

julia> @time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
  1.889226 seconds (1.58 M allocations: 938.699 MiB, 5.12% gc time)
=#

using Sundials, DiffEqDevTools, Plots
sol = solve(prob, CVODE_BDF(), reltol=1e-12, abstol=1e-12)
test_sol = TestSolution(sol)
abstols = @. 1/10^(6:2:12)
reltols = @. 1/10^(2:2:8)
setups = [
          Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>CVODE_BDF()),
         ]
wp1 = WorkPrecisionSet(prob,abstols,reltols,setups; save_everystep=false,appxsol=test_sol,seconds=1,numruns=1,print_names=true)
wp2 = WorkPrecisionSet(prob,abstols,reltols,setups; save_everystep=false,appxsol=test_sol,seconds=1,numruns=1,print_names=true,error_estimate=:l2)
plot(wp1)
savefig("$(homedir())/Desktop/heating_bench_finial.svg")
plot(wp2)
savefig("$(homedir())/Desktop/heating_bench_l2.svg")