Pandemic Control
A stochastic optimal control problem were we seek to control the spread of a contagion with uncertain parameters.
Problem Statement and Model
We wish to determine optimal social distancing policy to control spread of a contagion and minimize the economic impact of doing so. We'll model the spread of the virus through a given population using the SEIR Model which considers 4 population subsets that follow:
\[\text{Susceptible} \rightarrow \text{Exposed} \rightarrow \text{Infectious} \rightarrow \text{Recovered}\]
This model is formalized as:
\[\begin{gathered} \frac{ds(t)}{dt} = (u(t) - 1)\beta s(t)i(t) \\ \frac{de(t)}{dt} = (1 - u(t))\beta s(t)i(t) - \xi e(t) \\ \frac{di(t)}{dt} = \xi e(t) - \gamma i(t)\\ \frac{dr(t)}{dt} = \gamma i(t) \\ \end{gathered}\]
where $s(t)$ is the susceptible population, $e(t)$ is the exposed population, $i(t)$ is the infectious population, $r(t)$ is the recovered population, and $u(t) \in [0, 1]$ is the enforced population isolation (social distancing). The other values denote model parameters that will be specific to the contagion in question.
For our case study, we'll consider the incubation constant $\xi$ to be an uncertain parameter $\xi \sim \mathcal{U}(\underline{\xi}, \overline{\xi})$. This introduces to infinite dependencies into our model: time $t$ and uncertain incubation $\xi$. Moreover we'll seek to minimize the isolation measures $u(t)$ that are implemented while limiting the amount of infected individuals $i(t)$ to be below a threshold $i_{max}$. Thus, the optimization problem becomes:
\[\begin{aligned} &&\min_{} &&& \int_{t \in \mathcal{D}_{t}} u(t) dt \\ && \text{s.t.} &&& \frac{\partial s(t, \xi)}{\partial t} = (u(t) - 1)\beta s(t, \xi)i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&&&& \frac{\partial e(t, \xi)}{\partial t} = (1 - u(t))\beta s(t, \xi)i(t, \xi) - \xi e(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&&&& \frac{\partial i(t, \xi)}{\partial t} = \xi e(t, \xi) - \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&&&& \frac{\partial r(t, \xi)}{\partial t} = \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&&&& s(0, \xi) = s_0, e(0, \xi) = e_0, i(0, \xi) = i_0, r(0, \xi) = r_0, && \forall \xi \in \mathcal{D}_{\xi} \\ &&&&& i(t, \xi) \leq i_{max}, && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&&&& u(t) \in [0, 0.8] \\ &&&&& \xi \sim \mathcal{U}(\underline{\xi}, \overline{\xi}) \end{aligned}\]
Notice that the SEIR model variables now all depend on both $t$ and $\xi$, except $u(t)$ which means we need to decide our policy before knowing the true value of $\xi$.
Modeling in InfiniteOpt
Traditional modeling frameworks like JuMP cannot be used to solve this model directly because it is infinite, contains partial difference equations, and contains a time expectation. We would first have to transform it into a transcripted (discretized) variant by applying all the necessary techniques (e.g., orthogonal collocation over finite elements, trapezoid rule, etc.) which is very combersome and nontrivial in this case.
However, we can directly model the above form in InfiniteOpt and it will take care of the rest! Let's get started by importing the needed packages and specifying the model parameters that we'll need.
using InfiniteOpt, Ipopt, Distributions, Plots
# Set the SEIR parameters
γ = 0.303
β = 0.727
N = 1e5
ξ_min = 0.1
ξ_max = 0.6
# Set the domain information
i_max = 0.02
t0 = 0
tf = 200
extra_ts = [0.001, 0.002, 0.004, 0.008, 0.02, 0.04, 0.08, 0.2, 0.4, 0.8]
num_samples = 5
# Set the intial condition values
e0 = 1 / N
i0 = 0
r0 = 0
s0 = 1 - 1 / N;Model Initialization
Now let's setup our infinite model and select Ipopt as our optimizer that will be used to solve it. This is accomplished making an InfiniteModel:
model = InfiniteModel(Ipopt.Optimizer);Infinite Parameter Definition
We now need to define the infinite parameters $t \in [t_0, t_f]$ and $\xi \sim \mathcal{U}(\underline{\xi}, \overline{\xi})$. This is accomplished with @infinite_parameter. We'll also include the following specifications:
- use 51 equidistant time points
- specify that orgothonal collocation using 2 nodes should be used for time derivatives
- specify that the number of random scenarios should equal
num_samples - add
extra_tsas extra time points
@infinite_parameter(model, t ∈ [t0, tf], num_supports = 51,
derivative_method = OrthogonalCollocation(3))
@infinite_parameter(model, ξ ~ Uniform(ξ_min, ξ_max), num_supports = num_samples)
add_supports(t, extra_ts)Infinite Variable Definition
With our infinite parameters defined, we can now define our infinite variables:
- $s(t, \xi) \geq 0$
- $e(t, \xi) \geq 0$
- $i(t, \xi) \geq 0$
- $r(t, \xi) \geq 0$
- $0 \leq u(t) \leq 0.8$
@variable(model, s ≥ 0, Infinite(t, ξ))
@variable(model, e ≥ 0, Infinite(t, ξ))
@variable(model, i ≥ 0, Infinite(t, ξ))
@variable(model, r ≥ 0, Infinite(t, ξ))
@variable(model, 0 ≤ u ≤ 0.8, Infinite(t), start = 0.2)u(t)Objective Definition
Now its time to add the objective $\min \ \int_{t \in \mathcal{D}_{t}} u(t) dt$ using @objective:
@objective(model, Min, ∫(u, t))∫{t ∈ [0, 200]}[u(t)]Constraint Definition
The last step now to defining our model is to define the constraints using @constraint. This will involve defining the initial conditions:
\[s(0, \xi) = s_0, e(0, \xi) = e_0, i(0, \xi) = i_0, r(0, \xi) = r_0, \ \forall \xi \in \mathcal{D}_{\xi}\]
the model equations:
\[\begin{aligned} &&& \frac{\partial s(t, \xi)}{\partial t} = (u(t) - 1)\beta si(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&& \frac{\partial e(t, \xi)}{\partial t} = (1 - u(t))\beta si(t, \xi) - \xi e(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&& \frac{\partial i(t, \xi)}{\partial t} = \xi e(t, \xi) - \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ &&& \frac{\partial r(t, \xi)}{\partial t} = \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\ \end{aligned}\]
and the infection limit constraint:
\[i(t, \xi) \leq i_{max}, \ \forall t \in \mathcal{D}_t, \xi \in \mathcal{D}_{\xi}.\]
# Define the initial conditions
@constraint(model, s(0, ξ) == s0)
@constraint(model, e(0, ξ) == e0)
@constraint(model, i(0, ξ) == i0)
@constraint(model, r(0, ξ) == r0)
# Define the SEIR equations
@constraint(model, s_constr, ∂(s, t) == -(1 - u) * β * s * i)
@constraint(model, e_constr, ∂(e, t) == (1 - u) * β * s * i - ξ * e)
@constraint(model, i_constr, ∂(i, t) == ξ * e - γ * i)
@constraint(model, r_constr, ∂(r, t) == γ * i)
# Define the infection rate limit
@constraint(model, imax_constr, i ≤ i_max)imax_constr : i(t, ξ) ≤ 0.02, ∀ t ∈ [0, 200], ξ ~ UniformAdjust Degrees of Freedom
Since we are using OrthogonalCollocation, this introduces additional collocation points for t. These in turn, will artificially increase the degrees of freedom for u. Thus, we will treat u as a piecewise constant function that is held constant over the internal collocation nodes via constant_over_collocation:
constant_over_collocation(u, t)Display the Infinite Model
Let's display model now that it is fully defined:
print(model)Min ∫{t ∈ [0, 200]}[u(t)]
Subject to
s(t, ξ) ≥ 0, ∀ t ∈ [0, 200], ξ ~ Uniform
e(t, ξ) ≥ 0, ∀ t ∈ [0, 200], ξ ~ Uniform
i(t, ξ) ≥ 0, ∀ t ∈ [0, 200], ξ ~ Uniform
r(t, ξ) ≥ 0, ∀ t ∈ [0, 200], ξ ~ Uniform
u(t) ≥ 0, ∀ t ∈ [0, 200]
u(t) ≤ 0.8, ∀ t ∈ [0, 200]
s(0, ξ) = 0.99999, ∀ ξ ~ Uniform
e(0, ξ) = 1.0e-5, ∀ ξ ~ Uniform
i(0, ξ) = 0, ∀ ξ ~ Uniform
r(0, ξ) = 0, ∀ ξ ~ Uniform
s_constr : ∂/∂t[s(t, ξ)] - ((0.727 u(t)*s(t, ξ) - 0.727 s(t, ξ)) * i(t, ξ)) = 0, ∀ t ∈ [0, 200], ξ ~ Uniform
e_constr : ∂/∂t[e(t, ξ)] - (((-0.727 u(t)*s(t, ξ) + 0.727 s(t, ξ)) * i(t, ξ)) - (ξ*e(t, ξ))) = 0, ∀ t ∈ [0, 200], ξ ~ Uniform
i_constr : -ξ*e(t, ξ) + ∂/∂t[i(t, ξ)] + 0.303 i(t, ξ) = 0, ∀ t ∈ [0, 200], ξ ~ Uniform
r_constr : ∂/∂t[r(t, ξ)] - 0.303 i(t, ξ) = 0, ∀ t ∈ [0, 200], ξ ~ Uniform
imax_constr : i(t, ξ) ≤ 0.02, ∀ t ∈ [0, 200], ξ ~ Uniform
Optimize the Model
Let's solve our model and find the optimal policy. All we have to do is invoke optimize! and the model will automatically be transformed solved behind the scenes:
optimize!(model)This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.
Number of nonzeros in equality constraint Jacobian...: 17110
Number of nonzeros in inequality constraint Jacobian.: 605
Number of nonzeros in Lagrangian Hessian.............: 7260
Total number of variables............................: 4961
variables with only lower bounds: 2420
variables with lower and upper bounds: 121
variables with only upper bounds: 0
Total number of equality constraints.................: 4900
Total number of inequality constraints...............: 605
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 605
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0000000e+01 9.90e-01 2.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.7019358e+01 5.20e-01 1.08e+03 -1.0 9.90e-01 - 1.10e-02 4.75e-01f 1
2 1.4837666e+01 5.11e-01 1.06e+03 -1.0 2.31e+00 - 1.02e-01 1.72e-02h 1
3 1.4304342e+01 4.70e-01 9.75e+02 -1.0 1.22e+00 - 2.50e-01 7.98e-02h 1
4 1.4308729e+01 4.70e-01 9.75e+02 -1.0 7.75e-01 - 2.27e-02 7.90e-04h 1
5 1.4513243e+01 4.65e-01 9.64e+02 -1.0 3.35e+00 - 5.03e-04 1.09e-02f 1
6 1.4705078e+01 4.60e-01 1.57e+03 -1.0 3.31e+00 - 5.08e-03 9.99e-03h 1
7 1.5265087e+01 4.49e-01 1.34e+03 -1.0 2.29e+00 - 5.16e-03 2.46e-02f 1
8 1.5465044e+01 4.45e-01 1.34e+03 -1.0 2.67e+00 - 1.68e-02 8.06e-03h 1
9 1.5869343e+01 4.39e-01 1.32e+03 -1.0 2.29e+00 - 1.99e-02 1.49e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.6980915e+01 4.23e-01 1.24e+03 -1.0 2.27e+00 - 1.53e-02 3.59e-02h 1
11 1.8364740e+01 4.06e-01 1.18e+03 -1.0 1.97e+00 - 3.13e-02 3.93e-02h 1
12 2.2494526e+01 3.66e-01 1.02e+03 -1.0 1.41e+00 - 3.81e-02 9.97e-02h 1
13 3.4254284e+01 2.67e-01 6.79e+02 -1.0 1.01e+00 - 1.04e-01 2.70e-01h 1
14 6.8304296e+01 2.67e-03 1.07e+03 -1.0 4.77e-01 - 2.75e-01 9.90e-01h 1
15 6.5792608e+01 2.20e-03 8.43e+02 -1.0 6.05e-01 - 3.62e-01 1.77e-01h 1
16 6.1991025e+01 1.42e-03 5.39e+02 -1.0 8.55e-01 - 3.77e-01 3.53e-01h 1
17 5.9821080e+01 1.24e-03 8.14e+02 -1.0 8.93e-01 - 3.53e-01 2.91e-01h 1
18 5.9304981e+01 1.13e-03 4.86e+02 -1.0 1.05e+00 - 6.21e-02 8.82e-02h 1
19 5.9177595e+01 1.08e-03 5.16e+02 -1.0 1.01e+00 - 4.87e-02 4.39e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 5.7945313e+01 2.20e-03 3.40e+03 -1.0 7.75e-01 - 2.63e-01 5.85e-01h 1
21 5.8023146e+01 2.00e-03 1.22e+03 -1.0 9.10e-01 - 2.18e-01 1.57e-01h 1
22 5.8120016e+01 1.84e-03 2.74e+03 -1.0 9.51e-01 - 3.15e-02 8.23e-02h 1
23 5.9813444e+01 1.36e-03 6.44e+03 -1.0 8.73e-01 - 3.06e-01 4.47e-01h 1
24 6.0596980e+01 1.57e-03 1.40e+04 -1.0 8.49e-01 0.0 1.96e-01 3.53e-01h 1
25 6.0602918e+01 1.55e-03 1.23e+04 -1.0 3.39e-01 0.4 2.90e-02 1.38e-02h 1
26 6.0793877e+01 1.49e-03 8.79e+03 -1.0 3.94e+00 -0.1 6.87e-02 3.72e-02h 1
27 6.0792843e+01 1.47e-03 7.12e+03 -1.0 9.69e-01 0.4 2.90e-02 1.37e-02h 1
28 6.3564723e+01 1.58e-03 4.46e+04 -1.0 6.83e-01 -0.1 7.98e-02 4.77e-01h 1
29 6.3462869e+01 1.52e-03 8.08e+03 -1.0 8.10e-01 1.2 2.41e-01 3.72e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 6.3513417e+01 1.40e-03 4.36e+03 -1.0 1.17e+00 0.7 2.19e-01 1.86e-01h 1
31 6.3716045e+01 1.31e-03 4.93e+04 -1.0 1.01e+00 0.3 2.97e-01 1.11e-01h 1
32 6.3888462e+01 1.25e-03 5.21e+04 -1.0 1.44e+00 -0.2 6.76e-02 4.85e-02h 1
33 6.3872033e+01 1.21e-03 4.49e+04 -1.0 3.22e-01 0.2 1.65e-02 3.53e-02h 1
34 6.6215237e+01 8.08e-04 1.67e+04 -1.0 7.26e-01 -0.3 2.28e-01 3.79e-01h 1
35 6.7028304e+01 1.03e-03 4.23e+04 -1.0 7.94e-01 0.2 2.49e-01 3.15e-01h 1
36 6.7281056e+01 9.75e-04 2.32e+05 -1.0 7.80e-01 -0.3 5.41e-01 1.32e-01h 1
37 6.7352457e+01 9.35e-04 2.19e+05 -1.0 4.84e-01 -0.8 3.67e-02 4.04e-02h 1
38 6.7757245e+01 8.52e-04 2.02e+05 -1.0 1.95e+00 -0.4 9.30e-02 8.98e-02h 1
39 6.8134941e+01 7.47e-04 9.95e+04 -1.0 5.47e-01 -0.8 3.33e-02 1.29e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 6.9091265e+01 5.41e-04 3.23e+05 -1.0 6.77e-01 -1.3 5.55e-01 2.79e-01h 2
41 6.9620574e+01 6.74e-04 7.75e+04 -1.0 6.22e-01 - 3.93e-01 4.74e-01H 1
42 6.9533232e+01 6.01e-04 4.26e+04 -1.0 5.30e-01 - 8.69e-02 1.09e-01h 1
43 6.9559798e+01 5.52e-04 4.52e+05 -1.0 1.27e+00 - 2.97e-01 8.14e-02h 1
44 6.9706079e+01 4.08e-04 4.37e+04 -1.0 6.70e-01 - 7.25e-02 2.65e-01h 2
45 6.9349717e+01 4.17e-04 1.40e+05 -1.0 4.79e-01 - 6.85e-01 6.19e-01H 1
46 6.9127281e+01 3.64e-04 1.74e+04 -1.0 2.99e-01 - 3.50e-01 3.83e-01h 1
47 6.8931737e+01 2.80e-04 1.17e+05 -1.0 2.12e-01 - 1.94e-01 2.28e-01h 1
48 6.8903465e+01 2.03e-04 1.22e+06 -1.0 5.35e-01 - 7.23e-01 2.67e-01h 1
49 6.8738895e+01 1.35e-04 5.64e+05 -1.0 7.06e-02 - 2.74e-01 3.94e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 6.8698308e+01 1.37e-04 7.87e+02 -1.0 1.12e-01 - 1.00e+00 1.00e+00h 1
51 6.7722046e+01 8.79e-05 8.46e+05 -2.5 8.29e-02 - 8.62e-01 9.51e-01f 1
52 6.6701561e+01 1.90e-04 7.99e+04 -2.5 1.93e-01 0.0 7.87e-01 5.71e-01f 1
53 6.5746846e+01 1.34e-04 1.75e+05 -2.5 1.37e-01 0.4 2.15e-01 7.80e-01f 1
54 6.5164073e+01 2.05e-04 1.12e+05 -2.5 4.88e-01 -0.0 3.35e-01 2.31e-01f 1
55 6.4855737e+01 1.51e-04 2.48e+04 -2.5 1.18e-01 0.4 6.69e-01 3.08e-01f 1
56 6.4517658e+01 1.00e-04 2.14e+04 -2.5 1.05e-01 -0.1 2.44e-01 3.57e-01f 1
57 6.4365123e+01 4.91e-05 2.50e+04 -2.5 7.79e-02 0.3 1.13e-01 5.17e-01f 1
58 6.4194886e+01 7.93e-06 1.36e+03 -2.5 6.58e-02 -0.1 9.22e-01 8.40e-01f 1
59 6.4198166e+01 5.22e-07 5.74e+01 -2.5 2.24e-02 0.3 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 6.4179181e+01 6.20e-06 2.68e+01 -2.5 4.50e-02 -0.2 1.00e+00 1.00e+00f 1
61 6.4180167e+01 9.22e-06 2.66e+01 -2.5 6.48e+00 -0.7 1.60e-02 9.01e-03f 3
62 6.4167786e+01 8.50e-06 2.61e+00 -2.5 6.61e-02 -0.2 1.00e+00 1.00e+00f 1
63 6.3981334e+01 3.34e-06 1.49e+04 -3.8 3.38e-02 0.2 7.84e-01 7.28e-01f 1
64 6.3947507e+01 1.86e-05 5.89e+03 -3.8 3.20e-01 -0.3 1.62e-01 2.53e-01f 1
65 6.3904838e+01 5.22e-06 1.47e+04 -3.8 2.07e-02 0.1 1.00e+00 7.30e-01f 1
66 6.3888883e+01 1.41e-07 2.54e+00 -3.8 7.86e-03 0.6 1.00e+00 1.00e+00f 1
67 6.3882735e+01 1.81e-06 8.82e+00 -3.8 3.56e-02 0.1 1.00e+00 5.87e-01f 1
68 6.3880927e+01 2.24e-07 2.49e+00 -3.8 1.15e-02 0.5 1.00e+00 1.00e+00f 1
69 6.3871069e+01 5.31e-06 2.44e+00 -3.8 6.29e-02 0.0 1.00e+00 8.69e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 6.3868516e+01 3.67e-06 3.22e+01 -3.8 7.50e-02 -0.4 7.53e-01 3.25e-01f 1
71 6.3858974e+01 2.43e-05 1.43e+00 -3.8 9.15e-02 -0.0 1.00e+00 1.00e+00f 1
72 6.3844608e+01 2.96e-05 1.69e+01 -3.8 6.45e-01 -0.5 1.42e-01 8.01e-02f 1
73 6.3837642e+01 2.77e-05 1.50e+01 -3.8 1.94e-01 -0.1 1.00e+00 9.12e-02f 1
74 6.3827020e+01 2.59e-05 1.41e+01 -3.8 3.81e+00 - 1.47e-01 6.38e-02f 1
75 6.3821466e+01 2.06e-05 1.17e+01 -3.8 9.50e-01 - 9.01e-01 2.06e-01f 1
76 6.3812666e+01 6.11e-09 2.10e+00 -3.8 1.71e-01 - 1.00e+00 1.00e+00f 1
77 6.3812870e+01 3.18e-12 4.78e-04 -3.8 3.95e-02 - 1.00e+00 1.00e+00h 1
78 6.3802369e+01 2.58e-08 1.04e+03 -5.7 1.18e-01 - 8.44e-01 7.74e-01f 1
79 6.3799788e+01 2.62e-09 3.98e+02 -5.7 1.96e-01 - 7.82e-01 9.68e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 6.3799715e+01 6.11e-11 3.86e-04 -5.7 5.44e-03 - 1.00e+00 1.00e+00f 1
81 6.3799564e+01 3.38e-12 4.42e+00 -8.6 7.23e-03 - 9.92e-01 9.68e-01f 1
82 6.3799559e+01 3.14e-15 1.90e-07 -8.6 2.97e-04 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 82
(scaled) (unscaled)
Objective...............: 6.3799559298047420e+01 6.3799559298047420e+01
Dual infeasibility......: 1.9008555687816120e-07 1.9008555687816120e-07
Constraint violation....: 3.1378979276075469e-15 3.1378979276075469e-15
Variable bound violation: 8.7651495543848049e-09 8.7651495543848049e-09
Complementarity.........: 2.5172702981179954e-09 2.5172702981179954e-09
Overall NLP error.......: 2.5172702981179954e-09 1.9008555687816120e-07
Number of objective function evaluations = 98
Number of objective gradient evaluations = 83
Number of equality constraint evaluations = 98
Number of inequality constraint evaluations = 98
Number of equality constraint Jacobian evaluations = 83
Number of inequality constraint Jacobian evaluations = 83
Number of Lagrangian Hessian evaluations = 82
Total seconds in IPOPT = 4.972
EXIT: Optimal Solution Found.
Retrieve and Plot the Results
Now we can retrieve the optimal results and plot them to visualize the optimal policy. Note that the values of infinite variables will be returned as arrays corresponding to how the supports were used to discretize our model. We can retrieve our values using value.
Get the results:
r_opt = value(r, ndarray = true) * 100 # make the population fractions into percentages
s_opt = value(s, ndarray = true) * 100
i_opt = value(i, ndarray = true) * 100
e_opt = value(e, ndarray = true) * 100
u_opt = value(u)
obj_opt = objective_value(model)
ts = value(t)
ξs = value(ξ);Plot the values of $r(t, \xi)$ and $s(t, \xi)$ over time with confidence bounds:
r_mean = mean(r_opt, dims = 2)
r_std = std(r_opt, dims = 2)
plot(ts, r_mean, label = "r(t, ξ)", linecolor = :red, background_color = :transparent)
plot!(ts, r_mean + r_std, linecolor = :red, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, r_mean - r_std, linecolor = :red, linestyle = :dash, linealpha = 0.4, label = "")
s_mean = mean(s_opt, dims = 2)
s_std = std(s_opt, dims = 2)
plot!(ts, s_mean, label = "s(t, ξ)", linecolor = :blue)
plot!(ts, s_mean + s_std, linecolor = :blue, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, s_mean - s_std, linecolor = :blue, linestyle = :dash, linealpha = 0.4, label = "")
ylabel!("Pop. (%)")
xlabel!("Time (Days)")Plot the values of $i(t, \xi)$ and $e(t, \xi)$ over time with confidence bounds:
i_mean = mean(i_opt, dims = 2)
i_std = std(i_opt, dims = 2)
plot(ts, i_mean, label = "i(t, ξ)", linecolor = :green, background_color = :transparent)
plot!(ts, i_mean + i_std, linecolor = :green, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, i_mean - i_std, linecolor = :green, linestyle = :dash, linealpha = 0.4, label = "")
e_mean = mean(e_opt, dims = 2)
e_std = std(e_opt, dims = 2)
plot!(ts, e_mean, label = "e(t, ξ)", linecolor = :purple)
plot!(ts, e_mean + e_std, linecolor = :purple, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, e_mean - e_std, linecolor = :purple, linestyle = :dash, linealpha = 0.4, label = "")
ylabel!("Pop. (%)")
xlabel!("Time (Days)")Plot the values of $u(t)$ over time:
plot(ts, u_opt, linecolor = :orange, label = "u(t)", ylims = (-0.02, 1.02),
background_color = :transparent)
xlabel!("Time (Days)")
ylabel!("Distancing Ratio")Maintenance Tests
These are here to ensure this example stays up to date.
using Test
@test termination_status(model) == MOI.LOCALLY_SOLVED
@test has_values(model)
@test u_opt isa Vector{<:Real}Test PassedThis page was generated using Literate.jl.