Epidemiological modeling

[Adrilihan]

Hello everyone,

I was trying to use the OptimalControl.jl package to solve some epidemiological optimal control problems (OCP) and have encountered some issues. I considered the following OCP:

$$ \begin{cases} J(u(\cdot)) = \int_0^T (s(t) - i(t) - u(t)) \ dt \rightarrow \text{max}, \ \\ s'(t) = b(1 - s(t)) - (1 - u(t)) \beta (i(t) + \eta_C c(t) + \eta_A a(t))s(t) + d a(t) s(t), \\ i'(t) = (1 - u(t)) \beta (i(t) + \eta_C c(t) + \eta_A a(t)) s(t) - (\rho + \phi + b)i(t) + \alpha a(t) + \omega c(t) + d a(t) i(t), \\ c'(t) = \phi i(t) - (\omega + b)c(t) + d a(t) c(t), \\ a'(t) = \rho i(t) - (\alpha + b + d) a(t) + d a^2(t), \ \\ u \in [0, u_{\text{max}}], \ t \in [t_0 , t_f] \ \text{a.e.},\\ s(t_0) = s_0, \ i(t_0) = i_0, \ c(t_0) = c_0, \ a(t_0) = a_0. \end{cases} $$

Problem 1: Strange Behavior with Concatenation/Plotting

I implemented both the direct method and the shooting method. The shooting method consists of two parts: the first one approximates the switching time, and the second one obtains the initial adjoint state values. Despite obtaining a small shooting error after the second shooting, the plots of the adjoint variables do not agree with the error provided by the shooting function. Specifically, the plots of the adjoint at the final time are not close to zero (indicating that the transversality condition is not satisfied), even though the shooting variable norm and values are close to zero.

Interestingly, if I plot the solution without using the concatenation operator *, and instead plot the upper bang part and then modify the plot (using plot!) with the lower bang part, the plot appears as expected. Is this a bug?

Problem 2: Flow Function Convergence Error

I considered the same OCP with a quadratic cost function. This resulted in a bang-singular control solution. The direct method was implemented successfully, but the indirect method was not, because the resulting Hamiltonian system could not be solved by the Flow function. More specifically the Flow function was unable to compute the “singular flow” (flow with singular control), even when different solvers were used. Does anyone have suggestions on how to manage this issue?

Thank you in advance for your help!

[ocots]

For the first problem, it seems that using the MINPACK package without giving the Jacobian of the shooting function is not a good idea. Please consider to use: NonLinearSolve instead of MINPACK.