Hello everybody,
I am having some issues with the Moving Horizon Estimation (MHE) code I am writing. To me, the task and the code seem quite easy, but I cannot obtain the desired results for some reason.
I attach here the code that you can get also from my GitHub repo.
using MAT
using DataFrames
using Dates
using OrdinaryDiffEq, LinearAlgebra
using Optimization, OptimizationPolyalgorithms, OptimizationOptimJL
using SciMLSensitivity, Zygote, ForwardDiff
using Plots, LaTeXStrings
using Statistics, StatsBase
using OptimizationMOI, Ipopt
Npop = 59240329 # Total Population of Italy
N = 399 # Total data length
N_mhe = 21 # Estimation horizon (3 weeks)
Ts = 1. # Time steps for integration
# Create the "Date" vector
start_date = Date(2020, 8, 31)
end_date = Date(2021, 10, 3)
date = collect(start_date:end_date)
SIDTTHEfile = matopen("SIDTTHE_data_DEF.mat")
SIDTTHE_data = read(SIDTTHEfile, "SIDTTHE_data")
close(SIDTTHEfile)
I_data = SIDTTHE_data[1, 1]["data"] / Npop
D_data = SIDTTHE_data[2, 1]["data"] / Npop
T_data = (SIDTTHE_data[3, 1]["data"] + SIDTTHE_data[4, 1]["data"]) / Npop
H_data= SIDTTHE_data[5, 1]["data"] / Npop
E_data = SIDTTHE_data[6, 1]["data"] / Npop
S_data = ones(1,399)  (I_data + D_data + T_data + H_data + E_data)
# Gather all the data in a single data matrix
measure_mat = vcat(S_data, I_data, D_data, T_data, H_data, E_data)
function SIDTHEfun!(du, u, p, t)
s, i, d, t, h, e = u
Ī±, Ī“, Ī³, Ļ, Ļ = p
Ī» = 0.1 # coefficient Lambda (fixed)
du[1] = s * i * Ī± # dS/dt
du[2] = s * i * Ī±  (Ī³ + Ī») * i # dI/dt
du[3] = i * Ī³  d * (Ī» + Ī“) # dD/dt
du[4] = Ī“ * d  (Ļ + Ļ) * t # dT/dt
du[5] = (i + d) * Ī» + t * Ļ # dH/dt
du[6] = t * Ļ # dE/dt
end
u0 = [S_data[1]; I_data[1]; D_data[1]; T_data[1]; H_data[1]; E_data[1]]
Ī±0 = 0.25
Ī³0 = 0.12
Ī“0 = 0.01
Ļ0 = 0.02
Ļ0 = 0.02
p0 = [Ī±0, Ī³0, Ī“0, Ļ0, Ļ0]
# First simulation with initial condition parameters
tspan = (1.0, N_mhe)
tsteps = collect(range(1, stop = N_mhe, length = N_mhe))
prob = ODEProblem(SIDTHEfun!, u0, tspan, p0)
function simulate(p)
newprob = remake(prob, u0 = u0, p = p)
solsim = solve(newprob, Tsit5(), saveat = tsteps)
return reduce(hcat,solsim.u)
end
function loss(p)
sol = simulate(p)
loss = sum(abs2, ((sol. ymeas) ./ maximum.(eachrow(ymeas)) )) #= + sum(abs2, p  p_tilde) sum(abs2, sol[:,1]. x_tilde) =#
end
p_lb = [0.05, 0.005, 1e4, 1e4, 1e4] # lower bounds on parameters
p_ub = [0.8, 0.6, 0.6, 0.5, 0.5] # upper bound on parametersu0
states_dyn = [ S_data[N_mhe]; I_data[N_mhe]; D_data[N_mhe]; T_data[N_mhe]; H_data[N_mhe]; E_data[N_mhe] ]
params_dyn = [ Ī±0; Ī³0; Ī“0; Ļ0; Ļ0 ]
for k in N_mhe:NN_mhe
global u0
global ymeas = measure_mat[:, kN_mhe+1:k] #available measurement for each time window
if k == N_mhe
# u0 = u0 in the first iteration (for simulate(p) function)
global x_tilde = u0
global p_tilde = p0
else
u0 = reshape(opti_odesol[:,2],6,1) # initial condition for the new simulate(p) function
global x_tilde = reshape(opti_odesol[:,2],6,1)
global p_tilde = reshape(optires.u, 5, 1)
end
adtype = Optimization.AutoForwardDiff()
optf = Optimization.OptimizationFunction((x, p) > loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, p0, lb=p_lb, ub=p_ub) # Initial guess for parameter estimation
optires = Optimization.solve(optprob, Ipopt.Optimizer(), maxiters = 5000) # Ipopot Optimizer
opti_odesol = solve(remake(prob, u0 = u0, p = optires.u), Tsit5(), saveat = tsteps)
opti_odesol = reduce(hcat,opti_odesol.u)
# Save of the optimization values
params_dyn = hcat( params_dyn, reshape(optires.u, 5, 1) )
states_dyn = hcat( states_dyn, reshape(opti_odesol[:,end],6,1) )
# update on the initial conditions
p0 = reshape(optires.u, 5, 1)
end
#  Results Plot 
plot_font = "Computer Modern"
default(fontfamily=plot_font,
linewidth=2, framestyle=:box, label=nothing, grid=true)
scalefontsizes(1)
plotS = plot(date[N_mhe:NN_mhe], measure_mat[1,N_mhe:NN_mhe],label="Real Data")
plotS = plot!(date[N_mhe:NN_mhe], states_dyn[1,2:end],label="Sim Opti Params")
plotS = plot!(ylabel="Population", title="S  Susceptible population")
plotI = plot(date[N_mhe:NN_mhe], measure_mat[2,N_mhe:NN_mhe],label="Real Data")
plotI = plot!(date[N_mhe:NN_mhe], states_dyn[2,2:end],label="Sim Opti Params")
plotI = plot!(ylabel="Population", title="I  Infected population")
plotD = plot(date[N_mhe:NN_mhe], measure_mat[3,N_mhe:NN_mhe],label="Real Data")
plotD = plot!(date[N_mhe:NN_mhe], states_dyn[3,2:end],label="Sim Opti Params")
plotD = plot!(ylabel="Population", title="D  Detected population")
plotT = plot(date[N_mhe:NN_mhe], measure_mat[4,N_mhe:NN_mhe],label="Real Data")
plotT = plot!(date[N_mhe:NN_mhe], states_dyn[4,2:end],label="Sim Opti Params")
plotT = plot!(ylabel="Population", title="T  Threatened population")
plot_Ī± = plot(date[N_mhe:NN_mhe], params_dyn[1,2:end],label="Real Data")
plot_Ī± = plot!(ylabel="Population", title="Ī± parameter")
My goal here is to minimise the LSE between data and model estimation, obtaining a set of parameters p for every window.
The algorithm I am running here for the MHE works like this:

Initial parameter p0 and states u0 guess is made

The first LSE minimisation problem is solved for the past 21day time window

The results from the LSE are collected (in this case the set of parameters obtained) and the ODE system is resimulated in order to obtain state trajectory results.

The second step (t+1) minimisation initial conditions are set as follows: The initial condition for the states u0 is set as the 2nd step of the previously simulated trajectory, since the widows are shifting one step into the future, so the second step of the window @ t = t_k must be close to the initial conditions of the window @ t = t_k+1.
The initial parameter guess is set equal to the previously obtained parameter. 
Steps 345 are repeated shifting the window one day at a time and the FINAL VALUES of the state trajectory in each time window are collected and appended to the
states_dyn
vector at every time step.
The code runs, but it seems that the minimization is working only on one set of data (in this case, all data curves are way off the resulting prediction, but the T_data one), while the others are totally disregarded. I thought this was a normalisation error, but the data fitting process did not work even after normalising the data in the cost function (see sum(abs2, ((sol. ymeas) ./ maximum.(eachrow(ymeas)) ))
)
Moreover, I would like to add some more terms to the cited cost function, such as that initial conditions between consecutive windows must be close to one another and that parameter estimation results from consecutive windows should not vary too much (to avoid overfitting).
I am open to every tip/help on this topic, I am pretty new to the Julia Language and could have made a ton of errors, so if you have any critiques on how the code is structured, please let me know!!
Thanks for the help