Filtro de Kalman Extendido con input desconocido

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
epijacobian_full_u (generic function with 1 method)

Definimos un vector de condiciones iniciales $x_0$.

x0 = [7.112808e6, 1046.8508799517147, 0.0, 521.963080420307, 0.0]
F = 100. * ones(5)
5-element Array{Float64,1}:
 100.0
 100.0
 100.0
 100.0
 100.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 = [500.]
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.]
6×6 Array{Float64,2}:
 10000.0  10000.0  10000.0  10000.0  10000.0  0.0
 10000.0  10000.0  10000.0  10000.0  10000.0  0.0
 10000.0  10000.0  10000.0  10000.0  10000.0  0.0
 10000.0  10000.0  10000.0  10000.0  10000.0  0.0
 10000.0  10000.0  10000.0  10000.0  10000.0  0.0
     0.0      0.0      0.0      0.0      0.0  1.0

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

dt = 0.05
T = 30.
N = Int(T/dt)
600

Definimos los parámetros

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

Que agrupamos en una tupla p.

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

Usaremos un Discretizer RungeKutta de orden 4, diferenciable tanto en el estado $x$ como en el control $u$. Lo construimos a partir de un discretizador solo diferenciable en $x$ por comodidad.

rkx = KalmanFilter.RK4Dx(episystem_full, epijacobian_full_x, p, dt)
rk = KalmanFilter.RK4Du(rkx, epijacobian_full_u)
KalmanFilter.RK4Du(Main.ex-NLFilterUnknownInput.episystem_full, Main.ex-NLFilterUnknownInput.epijacobian_full_x, Main.ex-NLFilterUnknownInput.epijacobian_full_u, (gammae = 0.14, gammai = 0.07142857142857142, phi = 0.3), 0.05)

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")

Definimos las estructuras necesarias para crear un LinearKalmanIterator.

nlupdater = NLUpdater(rk,F,x0,0.4)
nlaugmented = KalmanFilter.NLUpdaterUnknowInput(nlupdater, control_pieces, 1.)
observer = KalmanFilter.LinearObserver(tildeH, zeros(1), G)
system = KalmanFilter.InnerState(tildex0)
iterator = KalmanFilter.LinearKalmanIterator(tildex0, tildeP, nlaugmented, observer, system, dt)
KalmanFilter.LinearKalmanIterator{Float64}(0, 1.0, KalmanFilter.InnerState([7.112808e6, 1046.8508799517147, 0.0, 521.963080420307, 0.0, 0.4]), KalmanFilter.ObservedState{Float64}([7.112808e6, 1046.8508799517147, 0.0, 521.963080420307, 0.0, 0.4], [10000.0 10000.0 … 10000.0 0.0; 10000.0 10000.0 … 10000.0 0.0; … ; 10000.0 10000.0 … 10000.0 0.0; 0.0 0.0 … 0.0 1.0]), KalmanFilter.ObservedState{Float64}([7.112808e6, 1046.8508799517147, 0.0, 521.963080420307, 0.0, 0.4], [10000.0 10000.0 … 10000.0 0.0; 10000.0 10000.0 … 10000.0 0.0; … ; 10000.0 10000.0 … 10000.0 0.0; 0.0 0.0 … 0.0 1.0]), KalmanFilter.NLUpdaterUnknowInput(NLUpdater(KalmanFilter.RK4Du(Main.ex-NLFilterUnknownInput.episystem_full, Main.ex-NLFilterUnknownInput.epijacobian_full_x, Main.ex-NLFilterUnknownInput.epijacobian_full_u, (gammae = 0.14, gammai = 0.07142857142857142, phi = 0.3), 0.05), [100.0, 100.0, 100.0, 100.0, 100.0], SimpleLinearUpdater{Float64}([0.9999978785448516 -0.014312165382640435 … 0.0 0.0; 2.114079715813942e-6 1.0072867332314537 … 0.0 0.0; … ; 2.210002323216013e-9 0.002103875043620381 … 0.9964349413940481 0.0; 2.2126297734661627e-9 0.0021076296453560557 … 0.0 1.0], [-37.72371824043088, 37.59256824275232, 0.09169598189541002, 0.03929827795517572, 0.039344999303566136], [100.0, 100.0, 100.0, 100.0, 100.0])), SimpleLinearUpdater{Float64}([0.9999978785448516 -0.014312165382640435 … 0.0 -37.72371824043088; 2.114079715813942e-6 1.0072867332314537 … 0.0 37.59256824275232; … ; 2.2126297734661627e-9 0.0021076296453560557 … 1.0 0.039344999303566136; 0.0 0.0 … 0.0 1.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [100.0, 100.0, 100.0, 100.0, 100.0, 1.0]), Main.ex-NLFilterUnknownInput.control_pieces, 0, 1.0), KalmanFilter.LinearObserver([0.0 0.0 … 1.0 0.0], [0.0], [500.0]), Normal{Float64}(μ=0.0, σ=0.0025000000000000005))

Y realizamos un total de N iteraciones, guardando los estamos intermedios en las variables que aparecen abajo.

results, ensamble = KalmanFilter.full_iteration(iterator, dt, N, t -> 0., 1)
(KalmanFilter.InnerStateSeries(600, 6, [3.828311111530479, 3.6694028193120056, 5.70602634512438, 9.100864342172564, 14.15076380671475, 13.628987049438187, 15.84537226233002, 19.653886357796473, 21.893973544021595, 26.684616690751543  …  1.5958371092777655e6, 1.599493104313237e6, 1.6031233162743628e6, 1.6067300228502874e6, 1.6103171293456592e6, 1.6138797196200336e6, 1.6174211692825188e6, 1.620942747341184e6, 1.6244379425942793e6, 1.6279121297325417e6], [7.112808e6 1046.8508799517147 … 0.0 0.4; 7.112770139114543e6 1077.2762477680553 … 2.2293077595865007 1.0; … ; 40289.05036544819 1.6704714945316857e6 … 1.6209408027960467e6 1.5; 39149.095383671454 1.6599544716559795e6 … 1.6244374895933757e6 1.5], [7.112808e6 1046.8508799517147 … 0.0 0.4; 7.112730940369932e6 1108.7214082512146 … 4.472800770438519 0.9997421027066561; … ; 40465.97091537528 1.670530082328433e6 … 1.6209411109312712e6 1.5034015696812968; 39320.824750573265 1.660018306440141e6 … 1.6244381217480376e6 1.5029905119023066], [NaN 0.0 … 0.0 0.0; 7.11279301728122e6 1054.4790188652134 … 2.2063739598856293 0.4; … ; 40485.87165569458 1.6705100395038363e6 … 1.620940552090514e6 1.4996440717859791; 39318.711191811184 1.6600204371379013e6 … 1.6244381824055e6 1.5034015696812968], [10000.0 10000.0 … 10000.0 1.0; 19416.98892362457 20331.857772943982 … 18554.948263728536 1.9629396461323803; … ; 1.9164828037539008e8 1.95748830596869e8 … 63642.02102919365 6.9201003050939605; 1.8412108287613156e8 1.8818035919905135e8 … 63234.59568147501 6.936881155939154], [NaN 8.4e-323 … 0.0 0.0; 20858.9481890159 21848.65486140269 … 20042.49834305342 2.0; … ; 2.1920974257152197e8 2.2370525599417603e8 … 85376.03458229633 7.902665283892239; 2.1011509032953084e8 2.145976504471678e8 … 84644.4178355826 7.9201003050939605], [0.0, NaN, 9.53333592e-315, 4.2439915824e-314, 5.141399935e-315, 0.25, NaN, 0.35, 0.4, 0.45  …  29.5, NaN, 29.6, 29.65, 29.7, 29.75, NaN, 29.85, 29.9, 29.95]), KalmanFilter.EnsamblesStoring([0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0]

...

[0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0]))

Resultados

Graficaremos los estados internos considerados, y los resultados obtenidos.

Ahora podemos graficar los diferentes estados de nuestro sistema, así como las aproximaciones obtenidas con filtro de Kalman.

Susceptibles $S$

plot(results, ts, 1)

Expuestos $E$

plot(results, ts, 2)

Infectados asintomáticos mild $I^m$

plot(results, ts, 3)

Infectados $I$

plot(results, ts, 4)

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

plot(results, ts, 5)

Finalmente, veamos el control usado y el aproximado

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

Notamos que tras una cierta cantidad de tiempo es posible averiguar el control Aunque hay bastante incerteza de los resultados.

This page was generated using Literate.jl.