Filtro de Kalman por ensambles

using Random, Distributions
using KalmanFilter
using Plots

Modelo epidemiológico SEII

El modelo permite modelar el desarrollo de una enfermedad regida por la siguiente dinámica:

\[\begin{aligned} S' &= - \alpha(t) S (E + I^m) \\ E' &= \alpha(t) S (E + I^m) - \gamma_E E \\ (I^m)' &= (1-\varphi) \gamma_E E - \gamma_I I^m \\ I' &= \varphi \gamma_E E - \gamma_I I\\ (cI)' &= \varphi \gamma_E E \end{aligned}\]

donde

• $S$ representa a los susceptibles, las personas que podrían contagiarse.
• $E$ corresponde a los que han sido expuestos al virus y lo están incubando.
• $I^m$ denota a los infectados mild, que tienen síntomas leves o ninguno.
• $I$ denota a los infectados sintomáticos.
• $cI$ es el total de infectados sintomáticos acumulados (agregamos este valor porque es un dato con el que se cuenta del sistema real).

Además la dinámica tiene los siguientes parámetros:

• $\alpha(t)$ es la Tasa de contagio en tiempo $t$.
• $\gamma_E$ es Tasa de infección (paso de expuesto a infectado).
• $\gamma_I$ es la Tasa de remoción (dejar de estar infectado, ya sea por recuperación o muerte).
• $\varphi$ es Proporción de personas sintomáticas.

Denotando $\mathbf{x}(t) = (S(t), E, I^m, I, cI)$, el sistema es de la forma

\[\mathbf{x}(t)' = f(\mathbf{x}(t), \alpha(t), p)\]

donde $\alpha(t)$ será un control y $p$ representa a los otros parámetros. Guardamos la información de la dinámica en la función episystem_full(x, α, p).

function episystem_full(x, α, p)
  γᵢ = p.gammai; γₑ = p.gammae; φ = p.phi;
  β = 1e-7
  [-α * β * x[1] * (x[2] + x[3]),
  α * β * x[1] * (x[2] + x[3]) - γₑ * x[2],
  (1 - φ) * γₑ * x[2] - γᵢ * x[3],
  φ * γₑ * x[2] - γᵢ * x[4],
  φ * γₑ * x[2]]
end

episystem_full (generic function with 1 method)

Vemos que este sistema es no lineal, por lo que para poder trabajarlo con el filtro de Kalman necesitaremos linearizarlo. Para eso será necesaria la información

\[D_x f(\mathbf{x}, \alpha(t), p) =\]

function epijacobian_full_x(x, α, p)
  γᵢ = p.gammai; γₑ = p.gammae; φ = p.phi;
  β = 1e-7
  [-α*β*(x[2] + x[3])  -α*β*x[1]     -α*β*x[1]  0.  0.;
  α*β*(x[2] + x[3])  (α*β*x[1]-γₑ)  α*β*x[1]  0.  0.;
  0.                 (1-φ)*γₑ       -γᵢ       0.  0.;
  0.                 φ*γₑ           0.        -γᵢ 0.;
  0.                 φ*γₑ           0.        0.  0.]
end

epijacobian_full_x (generic function with 1 method)

Puesto que consideraremos el input como desconocido, y la dinámica $f$ es no lineal en ese input, también necesitaremos $D_u f$.

function epijacobian_full_u(x, α, p)
  γᵢ = p.gammai; γₑ = p.gammae; φ = p.phi;
  β = 1e-7
  [-β*x[1]*(x[2] + x[3]),
  β*x[1]*(x[2] + x[3]),
  0.,
  0.,
  0.]
end



"""
Permite generar una muestra de condiciones iniciales ``x_{0}^{(1)}, \\dots, x_{0}^{(N)} ``

Genera una muestra a partir de un vector ``x0 = (S_0, E_0, I^m_0, I_0, cI_0)``. Los
elementos de la muestra suman la misma cantidad de personas. La muestra se obtiene
generando valores al azar que distribuyen ``\\mathcal{N}(E_0, \\beta E_0)``.
- `x0`: ``x_0``
- `factor`: ``\\beta``
- `N`: Cantidad de muestras ``N``.
"""
function get_initial_states(x0, factor, N, a0,unknowninput = false)
    total = sum(x0)
    expuestos_base = x0[2]
    expuestos = rand(Normal(expuestos_base, expuestos_base * factor), N)
    susceptibles = total .- expuestos
    if unknowninput
        states = [[susceptibles[i], expuestos[i], 0., 0., 0., rand(Normal(a0, factor * a₀ ))] for i in 1:N]
    else
        states = [[susceptibles[i], expuestos[i], 0., 0., 0.] for i in 1:N]
    end
    states
end

Main.ex-NLFilterENKF.get_initial_states

Definimos un vector de condiciones iniciales $x_0$.

x0 = [7.112808e6, 1046.8508799517147, 0.0, 521.963080420307, 0.0]
F = 50000. * ones(5)

5-element Array{Float64,1}:
 50000.0
 50000.0
 50000.0
 50000.0
 50000.0

Solo consideraremos conocida la cantidad total de infectados $\mathbf{x}_5 = cI$.

H = [0. 0. 0. 0. 1.]

1×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  1.0

Usaremos un error de unas 50.000 personas en las mediciones.

G = [5e6]
dimensions = 5

5

Definimos las matrices $\tilde{H}$ para nuestro sistema auxiliar, que incluirá al input.

tildeH = [H 0.]

tildeP = [F * F' zeros(5); zeros(5)' 1.]

dt = 0.05
T = 50.
N = Int(T/dt)

1000

Definimos los parámetros

γₑ = 0.14
γᵢ = 1/14
φ = 0.3

0.3

Que agrupamos en una tupla p.

p = (gammae = γₑ, gammai = γᵢ, phi = φ)


a₀ = 1.; Fₐ = 1.5

1.5

EL nuevo vector de estados será $\tilde{\mathbf{x}} = \begin{pmatrix} \mathbf{x}_0 \\ 1. \end{pmatrix}$, donde estamos considerando una suposición inicial para el estado $\alpha = 1$.

tildex0 = [x0; a₀]


rkx = KalmanFilter.RK4Dx(episystem_full, epijacobian_full_x, p, dt)
rk = KalmanFilter.RK4Du(rkx, epijacobian_full_u)

nlupdater = NLUpdater(rk,F,x0,a₀)

NLUpdater(KalmanFilter.RK4Du(Main.ex-NLFilterENKF.episystem_full, Main.ex-NLFilterENKF.epijacobian_full_x, Main.ex-NLFilterENKF.epijacobian_full_u, (gammae = 0.14, gammai = 0.07142857142857142, phi = 0.3), 0.05), [50000.0, 50000.0, 50000.0, 50000.0, 50000.0], SimpleLinearUpdater{Float64}([0.9999945820242319 -0.036165553298136734 … 0.0 0.0; 5.399272459857687e-6 1.0290644148908936 … 0.0 0.0; … ; 5.604354065717587e-9 0.0021265599999929084 … 0.9964349413940481 0.0; 5.610992492667885e-9 0.0021303415221729186 … 0.0 1.0], [-38.53691943239371, 38.403886626199196, 0.09301278948255273, 0.03986262406395118, 0.03990984185835405], [50000.0, 50000.0, 50000.0, 50000.0, 50000.0]))

El control para el sistema interno (que supondremos desconocido y es el que intentaremos averiguar) será la función constante por pedazos control_pieces.

function control_pieces(t)
    if t < 5.
        1.
    elseif t >= 5. && t < 12.
        0.8
    elseif t >= 12. && t < 25.
        0.5
    elseif t >= 25. && t < 40.
        1.5
    elseif t >= 40.
        1.
    end
end

control_pieces (generic function with 1 method)

Veamos un gráfico del control usado, para lo que definimos un vector con los tiempos ts.

ts = 0.0:dt:(T-dt)
plot(ts, control_pieces.(ts), label = "Control real")


#nlaugmented = KalmanFilter.NLUpdaterUnknowInput(nlupdater, control_pieces, Fₐ)
observer = KalmanFilter.LinearObserver(H, zeros(1), G)
#observer = KalmanFilter.LinearObserver(tildeH, zeros(1), G, tildex0)

ensemble_size = 50
states = get_initial_states(x0, 1.5, ensemble_size, a₀, false)
system = KalmanFilter.InnerState(x0)
enkf = KalmanFilter.EnKF(states, nlupdater, observer, system, dt)

results, ensamble = KalmanFilter.full_iteration(enkf, dt, N, control_pieces, ensemble_size)
#Juno.@enter KalmanFilter.full_iteration(enkf, dt, N, control_pieces)

Susceptibles $S$

plot(ensamble, ts, 1, ylims = (0, 7.5e6))
plot!(results, ts, 1)

Expuestos $E$

plot(ensamble, ts, 2)
plot!(results, ts, 2, ylims = (0, 2.8e6))

Infectados asintomáticos mild $I^m$

plot(ensamble, ts, 3)
plot!(results, ts, 3)

Infectados $I$

plot(ensamble, ts, 4)
plot!(results, ts, 4)

Infectados acumulados $cI$, y las observaciones que hicimos de él.

plot(ensamble, ts, 5)
plot!(results, ts, 5)

Control

#=
plot(ts, control_pieces.(ts), label = "Control real")
plot!(results, ts, 6)
plot!(ensamble, ts, 6)=#

This page was generated using Literate.jl.