key: cord-0187534-bxvrj8fn
authors: Lee, Wonjun; Liu, Siting; Li, Wuchen; Osher, Stanley
title: Mean field control problems for vaccine distribution
date: 2021-04-24
journal: nan
DOI: nan
sha: 340b232a79b53fefdafad52d3280091f3644ed65
doc_id: 187534
cord_uid: bxvrj8fn

With the invention of the COVID-19 vaccine, shipping and distributing are crucial in controlling the pandemic. In this paper, we build a mean-field variational problem in a spatial domain, which controls the propagation of pandemic by the optimal transportation strategy of vaccine distribution. Here we integrate the vaccine distribution into the mean-field SIR model designed in our previous paper arXiv:2006.01249. Numerical examples demonstrate that the proposed model provides practical strategies in vaccine distribution on a spatial domain.

1. Introduction. The COVID-19 pandemic has affected society significantly. Various actions are taken to mitigate the spread of the infections, such as the travel ban, social distancing, and mask-wearing. The recent invention of the vaccine yields breakthroughs in fighting against this infectious disease. According to the recent effectiveness study [17] , vaccines including Pfizer, Moderna, and Janssen (J&J) show approximately 66%-95% efficacy at preventing both mild and severe symptoms of COVID-19. Therefore, the deployment of COVID-19 vaccines is an urgent and timely task. Many countries have implemented phased distribution plans that give the elderly and healthcare workers priority to get vaccinated. Meanwhile, the shipping of vaccines are expensive due to the cold chain transportation [30] . An effective distribution strategy is necessary to eliminate infectious diseases and prevent more death.

In this work, we propose a novel mean-field control model based on [28] . We consider two approaches (controls) to control the pandemic: relocation of populations and distribution of vaccines. The first one has been discussed thoroughly in [28] , where we address the spatial effect in pandemic modeling by introducing a mean-field control problem into the spatial SIR model. By applying spatial velocity into the classical disease model, the model finds the most optimal strategy to relocate the different populations (susceptible, infected, and recovered), controlling the epidemic's propagation. We considered several aspects of the vaccine in our model for vaccine distribution, including manufacturing, delivery, and consumption. Our goal is to find an optimal strategy to move the population and distribute vaccines to minimize the total number of infectious, the amount of movement of the people, and the transportation cost of the vaccine with limited vaccine supply. To tackle this question, we ensemble these two controls and propose the following constrained optimization problem:

In our model, different populations are described using ρ i (i ∈ {S, I, R}), representing the susceptible, infectious, and recovered. The term ρ V (x, t) describes the density distribution of the vaccine over the spatial domain at location x and time t. The control variables v i (i ∈ {S, I, R}) create velocity fields over time-space domain that move the corresponding populations. As for vaccines, the control variable v V represents the vaccine's transportation strategy, and the control variable f (t, x) describes how many vaccines are produced at a specific time and location. The optimization objective function G is the sum of terminal costs E f inal and running costs E running . The terminal costs E f inal represent the goal of our control to achieve at the terminal time, such as minimizing the total number of infectious individuals and maximizing the total number of recovered (immune) persons. The running costs E running include the costs of transportation of vaccines and different classes of the populations, etc. We will discuss more details of cost functionals in Section 2.3. As for constraints of our optimization problem, the five partial differential equations of ρ i , v i (i ∈ {S, I, R, V }) describe the dynamics of the different classes of population and vaccines in terms of densities and velocities. The inequalities of f (t, x) model the limitation of vaccine manufacturing. Vaccines are only produced at particular factory locations Ω f actory with a daily maximal production rate f max . The dynamics of the vaccine density ρ V share some similar aspects to the unnormalized optimal transport [27] . Specifically, they both study mass transportation with a source term that creates masses.

We solve the main problem using the algorithm based on the first-order Primal-Dual Hybrid Gradient (PDHG) method [6, 7] . Due to the multiplicative interaction terms, ρ S K * ρ I , ρ I K * ρ S , ρ V ρ S , the optimization problem is based on nonlinear PDE constraints, whereas the PDHG only considers linear constraints. We use the extension of the PDHG [10] that solves nonsmooth optimization problems with nonlinear operators between function spaces. We extend the method utilizing the preconditioning operator from [21] where it provides a suitable choice of variable norms to achieve a convergence rate independent of the nonlinear operator. As a result, the algorithm converges to the saddle point locally with step length parameters independent of the finite-difference mesh size; see Section 3.1 for details.

Lots of mathematical models have been invented to predict the future of COVID-19 epidemics. Recently proposed models take more real-world situations into considerations and tend to be more effective in quantitative forecasting. Specifically, there have been studies on the impact of actions such as lockdown, social distancing, wearing a mask [13, 12, 16] . Data-driven approach and machine learning techniques are also integrated to estimate the parameters for the epidemic better and boost the prediction of the trend of the pandemic model [35, 32] . Meanwhile, optimal control serves as an important tool in pandemic control. They seek the optimal strategy to minimize the total number of infected people while keeping certain costs at a minimum. There are work focused on mitigating the epidemic with limited medical supply, such as ICU capacity [8] , face masks [31] , and vaccines [37, 19, 24, 29, 22] . In [22] , an optimal vaccine distribution strategy is proposed with a limited total amount of vaccines and maximal daily supply. [29] first uses inverse problem to determine the parameters of the SIR model. Then it formulates two optimal control problems, with mono-and multi-objective, and solves for the optimal strategy of vaccine administration. Other non-pharmaceutical interventions are also considered in the scope of optimal control of epidemics, including social distancing, closing schools, and lockdown [18, 23, 36] . [23] computes the optimal non-pharmaceutical intervention strategy based on an extended SEIR model with the absence of the vaccine. The mean-field control problem can be viewed as a particular type of optimal control, which is applied to an individual in terms of population density.

Mean-field game (control), introduced by [20, 26] , describes the deterministic (stochastic) differential games as the number of players tends to infinity, where a given player interacts through the distribution of all players in the state-space. It is a thriving research direction with applications in economics, crowd motion, industrial engineering, and more [3, 14, 25] . Numerical methods are invented to obtain quantitative information of such mean-field game (control) models, especially when the state-space is in high dimensions [1, 2, 5, 34] . Multi-population mean-field game (control) problems have also drawn lots of attention [4, 9, 15] . This type of problems studies the interactions in two levels: between agents of the same population and between populations. Our model is a multi-population mean-field control problem with population dynamics described using reaction-diffusion equations adopted from the epidemic model and the controls over the vaccine production and distribtuion. Therefore, we obtain a novel mean-field control problem.

The rest of the paper is organized as follows. Section 2 proposes a novel multipopulation mean-field control model and explains how population movement and vaccine distribution are integrated into a constrained optimization problem. Section 3 discusses the challenges in numerically solving this mean-field control model, proposes a first-order primal-dual algorithm to solve it, and shows the local convergence of the algorithm. Lastly, in Section 4, we present numerical experiments with different model parameter choices and discuss their implications on mean-field controls.

2. Models. In this section, we review the classical SIR model. Based on it, we formulate the spatial SIR dynamics with vaccine distribution, namely SIRV dynamics. We then introduce a variational problem to control the SIRV dynamics.

2.1. Classical SIR model. The SIR epidemic model describes an infectious disease epidemic via an ordinary differential equation system

The population is divided into three classes: susceptible, infected and recovered. While assuming a closed population without births or deaths, the model uses S(t), I(t) and R(t) to represent the number in each compartment at time t. The SIR model has two parameters: β is the effective contact rate of the susceptible individual being infected and γ is the recovery rate of the infected individual. The simplicity of this model allows people to predict an infectious disease epidemic by only estimating a few parameters. However, it has limitations by assuming the population is homogeneousmixing, which means that every individual has an equal probability of disease-causing contacts. As a result, the predictions will lack spatial information and may not help the (local) governments make policies or relocate medical resources. Therefore, we are motivated to study the spatial SIR model. On the other hand, the SIR model does not consider the latent period between when a person is exposed to a disease and when he or she becomes infected. This leads to the extension of the SIR model, such as the SEIR model. Our proposed model has a flexible structure and can be generalized to such epidemiological models naturally.

2.2. Spatial SIR variational problem with vaccine distribution. In [28] , we add the spatial dimension to the S, I, R functions. Let Ω ⊂ R d be a bounded domain. Consider the following density functions

Here, ρ S , ρ I , and ρ R represent susceptible, infected, and recovered populations distribution, respectively. We assume ρ i for each i ∈ {S, I, R} moves over a spatial domain Ω with a velocity v i . Here v i are our controls variables. With changes of variables m i = ρ i v i , we define the momentum m S , m I , m R : [0, T ] × Ω → R d that govern the corresponding density flows. In the following, instead of using control variables v i , we replace them with mi ρi and regard m i as the control variables. We can describe the flows of the densities by the following continuity equations.

, ρ I (0, ·), ρ R (0, ·) are given.

This system of continuity equations describes the flows of three groups of densities while satisfying the SIR model. The nonnegative constants η i (i ∈ {S, I, R}) are the coefficients for viscosity terms. These terms can also be understood as noise terms with noise generated by the data. K = K(x, y) is a symmetric positive definite kernel with (K * ρ)(x, t) = Ω K(x, y)ρ(y, t) dy. In this model, we consider the Gaussian kernel

The kernel convolution describes the spreading rate of the infectious disease over the spatial domain. In addition, we assume the Neumann boundary conditions on ∂Ω.

Since we don't consider birth or death in our model, the total population is conserved for all time t ∈ [0, T ], which leads to the following equality

In this paper, we consider the optimization problem for the distribution of vaccines. We add an extra function ρ V : [0, T ] × Ω → [0, ∞) which represents the vaccine density in Ω at each time t ∈ [0, T ]. The vaccine distribution will be described as the following PDE:

where m V : [T , T ) × Ω → R d is a momentum, θ 2 represents the utilization rate of vaccines, and f : (0, T ) × Ω → [0, ∞) represents the production rate of vaccines in x ∈ Ω at 0 < t < T . During 0 < t < T , the vaccines are produced with a production rate f and used at a rate θ 2 ρ V ρ S . During T ≤ t < T , the vaccines are delivered to the area where the susceptible population is located, and they are used at a rate of θ 2 ρ V ρ S . In summary, the first part of the PDE describes vaccines' production, and the second part describes the delivery of vaccines. For all time 0 < t < T , the susceptible population are vaccinated if the vaccines are available in the same area. Now we are ready to introduce the new system of equations for the SIRV model.

In the first and third equations, we add the terms −θ 1 ρ V ρ S and +θ 1 ρ V ρ S , respectively. The constant θ 1 represents the vaccine efficiency and θ 1 ρ V (t, x)ρ S (t, x) represents the vaccinated population at (t, x) ∈ (0, T ) × Ω. We denote a set S := {S, I, R, V } and define a nonlinear operator A as follows

where X C : [0, T ] → R is a step function that equals 1 on C and 0 otherwise.

2.3. The cost functional. The cost functional we propose in this paper is the extension of [28] . We design the cost functional so that the solution (ρ i , m i ), i ∈ S satisfies the following criteria:

(i) minimize the transportation cost for moving each population; (ii) minimize the total number of infected people and the total number of susceptible people by maximizing the usage of the vaccines at time T ; (iii) maximize the total number of recovered people at time T ; (iv) avoid high concentration of population and vaccines at each time t ∈ (0, T ); (v) minimize the amount of vaccines produced during t ∈ (0, T ); (vi) minimize the transportation cost for delivering vaccines during t ∈ (T , T ).

which is convex, lower semi-continuous, and 1-homogeneous with respect to (ρ i , m i ). The parameter α i characterizes the cost of moving ρ i with velocity mi ρi . Larger α i means it is more expensive to move ρ i . Note that this function comes from the quadratic kinetic energy. To see this, we use the definition m i = ρ i v i and plug into the formula (2.5):

where functions e : [0, ∞) → [0, ∞) are convex and lower semi-continuous functions. We also minimize the terminal cost for ρ V because maximizing the usage of vaccines is equivalent to minimizing the number of vaccines left at the terminal time T . The total number of the recovered can be maximized by penalizing the density at the terminal time if the value of ρ R (T, x) is far away from 1 for x ∈ Ω. In this paper, we use a quadratic cost function

where a i is some constant. For Item (iv), the cost functional for the concentration of the total population and vaccines can be represented by

for u : Ω → [0, ∞) and convex and lower semi-continuous functions

Similar to e i (2.6) from Item (ii), we use quadratic functions for g P and g V . Items (v) and (vi) are criteria specific to the vaccine distribution. From the PDE (2.2), the vaccines are produced during 0 < t < T by a function f . We use the similar functional (2.7) to minimize the amount of vaccines produced by f . Thus, we set the functional

is a convex and lower semi-continuous function.

The vaccines are delivered during T < t < T . Similar to the Item (i), we set

where F V has the same definition as (2.5).

The total cost functional we consider is then

In the perspective of a control problem, the first term at the right-hand side in (2.8) is the terminal cost, while the rest of the terms accounts for the running costs. The quadratic terms in the last line is a λ-strongly convex functional. The functional F is λ-strongly convex if for any u = (

and ∂F denotes the convex subdifferential of F . Since E i , F i , G i are convex and lowersemicontinuous, G is λ-strongly convex as the sum of convex and λ-strongly convex functionals. The strong convexity of G is important as the algorithm of the paper requires the objective cost functional to be strongly convex (Theorem 3.3).

In addition to the constraint from (2.3), we adapt the following constraints to reflect the limited vaccination coverage:

where Ω f actory ⊂ Ω indicates the factory area where vaccines are produced and f max is a nonnegative constant representing the maximum vaccine production rate. In the third inequality, a nonnegative constant C f actory limits the total number of vaccines produced during 0 < T < T .

The constraints (2.9) can be imposed by having the following functionals for G V and G 0 .

(2.10)

where Ω f actory ⊂ Ω indicates the factory area where vaccines are produced. The

where a, b are constants and u : Ω → R is a function. The function i Ω f actory (x) is defined as

This function forces f (t, x) = 0 if (t, x) ∈ (0, T ) × (Ω\Ω f actory ), thus vaccines are produced only in Ω f actory .

Remark 2.1. The formulation is not limited to SIR epidemic model. For example, we can describe the SIRD (Susceptible-Infected-Recovered-Deceased) epidemic model by adding an extra population ρ D for the deceased population with a mortality rate µ.

From the definition of the cost functional and the constraint (2.3), we have the following minimization problem:

We first define the inner product of vectors of functions in L 2 . Given vectors of

× Ω → R, the L 2 inner product of vectors u and v and L 2 norm of u are defined by

We introduce dual variables (φ i ) i∈S for each continuity equation from (2.4). Using the dual variables and the definitions of the inner products, we convert the minimization problem into a saddle point problem.

where L is the Lagrangian functional defined as

For brevity, we denote

We can rewrite the Lagrangian as

where the nonlinear operator A(u) is defined as

As noted in [28] , the dual gap, the difference between the primal solution and dual solution, may not be zero because the nonconvex functions (ρ S , ρ I ) → ρ S K * ρ I and (ρ S , ρ V ) → ρ S ρ V make the feasible set nonconvex. We circumvent the problem by linearizing the nonlinear operator at a base pointū

In our formulation, the linearlized operatorĀū(u) can be written as follows.

We define a linearized Lagrangian as

In the paper [10] , the author developed a primal-dual algorithm using the linearized Lagrangian (Algorithm (3.5)) and proves that the sequence (u (k) , p (k) ) ∞ k=1 from the algorithm converges to the saddle point (u * , p * ) (in Section 3.1, we prove the local convergence to the saddle point given (u (0) , p (0) ) is sufficiently close to the saddle point). By the first-order optimality conditions (also known as Karush-Kuhn-Tucker (KKT) conditions), the saddle point satisfies

In the next proposition, we present the equations derived from the KKT conditions (2.16).

Proposition 2.2 (Mean-field control SIRV system). By KKT conditions, the saddle point ((ρ i , m i , φ i ) i∈S , f ) of (2.12) satisfies the following equations.

The terms δG P δρ , δG V δρ , δG P δρ , δG0 δf , and δEi δρ(T,·) are the functional derivatives. In other words, given F : H → R be a smooth functional where H is a separable Hilbert space and ρ ∈ H, we say a map δF δρ is the functional derivative of F with respect to ρ if it satisfies

for any arbitrary function h : Ω → R.

The dynamical system models the optimal vector field strategies for S, I, R populations and the vaccine distribution. It combines both strategies from mean field controls and SIRV models. For this reason, we call it Mean-field control SIRV system. The proof of Proposition 2.2 can be found in the Appendix.

3. Algorithms. In this section, we propose an algorithm to solve the proposed SIRV variational problem. We use the primal-dual hybrid gradient (PDHG) algorithm [6, 7] . The PDHG can solve the following convex optimization problem.

where f and g are convex functions and A is a continuous linear operator. The algorithm solves the problem by converting the problem into a saddle point problem by introducing a dual variable p.

with L 2 inner product is defined in (2.11) and

is the Legendre transform of f . The method solves the saddle point problem by iterating

The scheme converges if the step sizes τ and σ satisfy

where · is an operator norm in L 2 . However, the SIRV variational problem has a nonlinear function A for the constraint. Thus, we use the extension of the algorithm from [10] which solves the nonlinear constrained optimization problem.

where A is a nonlinear function. The scheme iterates the algorithm (3.1) with a linear approximation of A at a base pointū

Denote A u := ∇A(u). We have a linearized saddle point problem

and the scheme iterates

The paper [10] proves that the sequence {u (k) , p (k) } ∞ k=0 of the algorithm converges to some saddle point (u * , p * ) that satisfies (2.16) . However, the scheme converges if the step sizes satisfy σ (k) τ (k) ∇A(u (k) ) 2 L 2 < 1, k = 1, 2, · · · . Suppose we use an unbounded operator that depends on the grid size, for example, A = ∇. Then the scheme can result in a very slow convergence if we use a fine grid resolution. To circumvent the problem, we use the General-proximal Primal-Dual Hybrid Gradient (G-prox PDHG) method from [21] which is another variation of the PDHG algorithm. This variant provides an appropriate choice of norms for the algorithm, and the authors prove that choosing the proper norms allows the algorithm to have larger step sizes than the vanilla PDHG algorithm. The G-prox PDHG iterates (3.6)

where the norm · H (k) is defined as

By choosing the proper norms, the step sizes only need to satisfy

which are clearly independent of the grid size.

In this section, we show the iterations from the algorithm (3.6) locally converges to the saddle point. The local convergence theorem in this paper is mainly based on the Theorem 2.11 from [10] . However, we add a preconditioning operator from the G-prox PDHG method. We show that the method converges locally to the saddle point with the step sizes independent of the nonlinear operator A.

From the algorithm (3.6), (u (k+1) , p (k+1) ) satisfies the following first-order optimality conditions

which can be rewritten as

with q = (u, p). Here, the monotone operator Hū is defined as

where Id is an identity operator. Recall that from (2.16), the saddle point q * = (u * , p * ) has to satisfy 0 ∈ H u * (u * , p * ).

Throughout, we assume that (3.9) ∇A(u * ) > 0.

Lemma 3.1. There exists constants 0 < c < C and R > 0 such that

where A is a nonlinear operator from (2.14) and · is an operator norm.

Proof. This follows immediately from (3.9) and the fact that the derivative ∇A(u) is continuous with respect to u.

Lemma 3.2. Suppose (3.9) holds and let τ (k) σ (k) < 1. Then there exist constants 0 < θ < Θ such that

A proof of Lemma 3.2 is provided in the appendix.

With the above Lemmas, we can use the Theorem 2.11 from [10] to show the local convergence of the algorithm. Theorem 3.3. Let (u * , p * ) ∈ L 2 × H ( * ) be a solution to (2.16) where p 2 H ( * ) = A T u * p 2 L 2 . Let the step sizes τ (k) and σ (k) satisfy τ (k) σ (k) < 1 for all k. Then there exists δ > 0 such that for any initial point (u (0) , p (0) ) ∈ L 2 × H (0) satisfying

the iterates (u (k) , p (k) ) from (3.6) converges to the saddle point (u * , p * ).

Proof. By Lemma 3.1, Lemma 3.2, and strong convexity of the functional G from (2.8), we can use [10, Theorem 2.11], which proves the theorem.

Remark 3.4. [10, Theorem 2.11] requires H u * to satisfy the condition called metric regularity. In our formulation, the constraint A(u) = 0 makes H u * metrically regular by [11, Section 5.3] . We refer readers to [10, 11, 33] for further details about metric regularity.

To implement the algorithm to the minimization problem (2.8), we set

We use (2.14) for the definition of the operator A. Define the Lagrangian functional as

where ·, · L 2 is defined in (2.11). We summarize the algorithm as follows.

Algorithm 1: G-prox PDHG for mean-field control SIRV system

Here, L 2 and H (k) i norms are defined as

for any u : [0, T ] × Ω → [0, ∞). Moreover, the relative error is defined as

In the section 4, We use quadratic functions for E i (i ∈ {S, I, V }), G P , G V , G 0 . With the definitions (2.10), we use

Thus, we can write the cost functional as follows

where a i , d P , d V , d 0 are nonnegative constants. With this cost functional, we find explicit formula for each variable ρ

, and f (k+1) from the Algorithm 1 satisfy the following explicit formulas:

where root + (a, b, c) is a positive root of a cubic polynomial x 3 + ax 2 + bx + c = 0 and we approximate the A i A * i as follows

We use FFTW library to compute (A i A T i ) −1 (i ∈ S) and convolution terms by Fast Fourier Transform (FFT), which is O(n log n) operations per iteration with n being the number of points. Thus, the algorithm takes just O(n log n) operations per iteration.

In this section, we present several sets of numerical experiments using the algorithm 1 with various parameters. We wrote C++ codes to run the numerical experiments. Let Ω = [0, 1] 2 be a unit square in R 2 and the terminal time T = 1. The domain [0, 1]×Ω is discretized with the regular Cartesian grid below.

x kl = ((k + 0.5)∆x 1 , (l + 0.5)∆x 2 ) , k = 0, · · · , N x1 − 1, l = 0, · · · , N x2 − 1 t n = n∆t, n = 0, · · · , N t − 1 where N x1 , N x2 are the number of discretized points in space and N t is the number of discretized points in time. For all the experiments, we use the same set of parameters, By setting a higher value for α I , we penalize the infected population's movement more than other populations. Considering the immobility of the infected individuals, this is a reasonable choice in terms of real-world applications. By setting T = 1/2, the solution will produce the vaccines during 0 ≤ t < 1/2 and deliver them during 1/2 ≤ t ≤ 1. Furthermore, we fix the parameters for the infection rate and recovery rate β = 0.8, γ = 0.1.

The paper [28] describes how the parameters β and γ affect the propagation of the populations. In this paper, we focus on the vaccine productions and distributions.

Recall that from the formulation (3.10), we have terminal functionals

Thus, the solution to the problem has to minimize the total number of susceptible, infected, and vaccines at the terminal time T . The solution reduces the total number of infected by recovering them with a rate γ and decreases the total number of susceptible by transforming the susceptible to the infected with a rate β or to the recovered with a rate θ 1 (Figure 1 ). If the β is large and γ is small, the number of infected will grow since there are more inflows from susceptible than the outflows to the recovered. To minimize the total number of the infected, the solution has to vaccinate the susceptible as much as possible to avoid the susceptible becoming infected. Thus, the vaccines need to be produced and delivered to the susceptible efficiently while satisfying the constraint conditions (2.9). We present two experiments that demonstrate how the various factors in the formulation affect the production and the distribution of vaccines.

In this experiment, we show how the parameters related to ρ V affect the solution. We set the initial densities for the ρ i (i ∈ S) and the factory location Ω f actory as 0.3, 0.3) where (x) + = max(x, 0) and B r x is a ball of a radius r centered at x. With the initial densities (4.1), we run two simulations with different values for θ 1 , θ 2 , and f max .

Sim 1 Sim 2 Description θ 1 0.5 0.9 Vaccine efficiency f max 0.5 10 Maximum production rate of vaccines C f actory 0.5 2 Maximum amount of vaccines that can be produced at x ∈ Ω during 0 ≤ t ≤ 1 2 Figure 3 shows the comparison between the results from the simulation 1 and the simulation 2. The first three plots (Figure 3a) show the total mass of ρ i (i = S, I, R), i.e. 

The total number of vaccines produced from the simulation 1 is smaller than that from the simulation 2 because the solution cannot produce a large amount of vaccines due to the low production rate f max . Furthermore, the solution from the simulation 1 cannot vaccinate a large number of susceptible due to a small θ 1 . Thus, there are more susceptible and less recovered at the terminal time in the simulation 1.

This experiment includes the spatial obstacles and shows how the algorithm effectively finds the solution that utilizes the vaccine production and distribution given spatial barriers. Denote a set Ω obs ⊂ Ω as obstacles. We use the following functionals in the experiment.

The densities ρ i (i ∈ S) cannot be positive on Ω obs due to i Ω obs . Thus, the densities transport while avoiding the obstacle Ω obs . We show two sets of experiments based on this setup.

We set the initial densities and Ω f actory as follows

and fix the parameters

The initial densities are shown in Figure 4 . Figure 6 show the evolution of densities with and without obstacles, respectively. In both simulations, the density of vaccines ρ V (the fourth row) transports to the areas where the susceptible people are present. In Figure 6 , ρ V transports while avoiding the obstacle at the right. Figure 7 shows the comparison between these two solutions and how the presence of the obstacle affects the production and delivery of vaccines quantitatively. Figure 7a shows the total mass of the vaccines in the factory area Ω f actory during the production time ρ V starts to transport at time t = 0.5, the number of susceptible is lower at the left. Thus, the solution distributes fewer vaccines to the left with less susceptible people. When there is an obstacle, ρ V has to bypass the obstacle to reach the susceptible areas. Thus, the kinetic energy cost during the delivery time t ∈ [0.5, 1] increases at the right. The solution cannot deliver the vaccines as much as the case without the obstacle. It results in a fewer number of vaccines produced during t ∈ [0, 0.5) (Figure 7a ) and delivered to the right during t ∈ [0.5, 1] when there is an obstacle (Figure 7b ).

Ω f actory ρ V (t, x) dx, t ∈ [0, 0.5).

Similar to the previous experiment, we show how the obstacles in the spatial domain affect the production and distribution of the vaccines. We use more complex initial densities, an obstacle set Ω obs , and three factory locations in this experiment. We set the initial densities and Ω f actory as follows The initial densities are shown in Figure 8 . Figure 9 and Figure 10 show the evolution of densities with and without obstacles, respectively. The experiment demonstrates that even with the complex initial densities, the algorithm successfully converges to the reasonable solution that coincides with the previous experiments. The density of vaccines ρ V (the fourth row) transports to the areas where the susceptible people are present while avoiding the obstacles. Figure 11a shows the total mass of the vaccines produced during the production time at each factory location. Without the obstacles, the total mass of ρ V at the middle is the lowest at time 0.5 because the factory at the middle is the farthest away from the susceptible people. It is more efficient to produce the vaccines at the factories closer to the susceptible (the top and the bottom factories) to reduce the kinetic energy cost during the delivery time t ∈ [0.5, 1]. However, with the obstacles, the vaccines are produced the most at the middle factory. Since the obstacles are blocking the paths between the top and the bottom factories and the susceptible people, ρ V has to bypass them to reach to the target area. The pathways from the middle factory to the susceptible people are not blocked as much as from the top and the bottom factories. Thus, it is more efficient to produce more vaccines at the middle factory. Figure 11b shows the total mass of the vaccines during the delivery time at dif- ferent locations. The lines in the plot represent the following quantities:

With the obstacles, the kinetic energy cost increases since ρ V has to bypass to reach to the targets when it transports from the top and the bottom factories. As a result, the vaccines are not produced as much as the simulation without the obstacles, and there are less vaccines reached to the targets.

This experiment compares the vaccine production strategy generated by the algorithm and the strategy with the fixed rates of production without using the algorithm. The initial densities and Ω f actory are set as follows We fix the parameters θ 1 = 0.9, f max = 5, C f actory = 1.

The initial densities and locations of factories are shown in Figure 12 .

To fairly compare the effect of the optimal vaccine production strategy, we remove the momentum of S, I, R groups; thus, removing the spatial movements defined by m S , m I , m R . We consider the following PDEs:

Furthermore, by taking out the momentum terms from S, I, R groups, the cost functional for this experiment is

With the PDEs and the cost functionals above, we compare two results. The first result is using the optimal vaccine production and distribution strategy generated by the algorithm 1. The second result is using the fixed vaccine production rate and the algorithm's distribution strategy. In the second result, the factory variable f is fixed as Figure 13 shows the comparison between these two results. The result from the fixed production rate is "without control", and the result from the optimal vaccine production strategy is "with control". The labels "left", "middle", and "right" are the locations of the factories in Figure 12 . The solid lines, the result with the same fixed rates of productions, show that all three factories produce the same amounts of vaccines. The dotted lines show the least amount of vaccines in the middle factory and much more in the left and right factories. When vaccines produce at the middle factory, one needs to pay more transportation costs because they bypass the obstacles. The obstacle does not block the paths from the left and right factories to the susceptible. Thus, it's an optimal choice to utilize the left and right more than the middle to minimize the transportation costs. The table below is the quantitative comparison between the two results. 0 Ω ρV dx dt during production t ∈ [0, T ] at each factory location: left, middle, and right. The dotted lines are from the optimal strategy from the Algorithm 1, and the solid lines are from the fixed production rates.

Description Algorithm 1 Fixed rates

The total amount of vaccines produced. 7.997 × 10 −3 8.411 × 10 −3

The number of susceptible people at the terminal time.

1.520 × 10 −2 1.525 × 10 −2

The number of infected people at the terminal time. 5.133 × 10 −3 5.134 × 10 −3

The transportation cost of vaccines. 7.339 × 10 −3 7.544 × 10 −3

The first row of the table shows that more vaccines are produced with a fixed rate of production. However, the result from the fixed-rate vaccinizes fewer susceptible people; as a result, more infected people at the terminal time. Furthermore, the result from the fixed rate shows higher transportation costs. It shows that the algorithm finds the more efficient strategy with fewer vaccines produced.

Proof of Proposition 2.2. From the saddle point problem (2.12), we can rewrite the problem as inf (ρi,mi) i∈S ,f

where a function Q : (0, T ) × Ω → R is defined as

If ((ρ i , m i , φ i ) i∈S , f ) is the saddle point of the problem, the differential of Lagrangian with respect to ρ i , m i , φ i (i ∈ S), f and ρ i (T, ·) (i ∈ {S, I, V }) equal to zero. Thus, from δL δφi = 0 we have

Using integration by parts, we reformulate the Lagrangian function (5.1) as follows.

From δL δρi = 0 (i ∈ {S, I, R}),

From δL δρi(T,·) = 0 (i ∈ S), δE δρ i (T, ·) (ρ i (T, ·)) = φ i (T, ·).

From δL δf = 0, δG 0 δf (f ) + φ V = 0, (t, x) ∈ (0, T ) × Ω.

From δL δmi = 0 (i ∈ S), α i m i ρ i = −∇φ i (t, x) ∈ (0, T ) × Ω, i ∈ {S, I, R}

By replacing αimi ρi = −∇ρ i in δL δρi = 0 and δL δφi = 0, we derive the result. Proof of Lemma 3.2. Let q = (u, p). By the definition of M (k) , we have

Using Young's inequality and Lemma 3.1,

We are left to show the lower bound. Let > 0 be such that τ (k) σ (k) = (1 − ) 2 . Then using Hölder's inequality,

Again, using Young's inequality and Lemma 3.1,

This proves the claim.

Mean field games: numerical methods

Mean field games and applications: Numerical aspects

Optimal incentives to mitigate epidemics: a stackelberg mean field game approach

Mean field control and mean field game models with several populations

On the implementation of a primal-dual algorithm for second order time-dependent mean field games with local couplings

A first-order primal-dual algorithm for convex problems with applications to imaging

On the ergodic convergence rates of a first-order primaldual algorithm

Covid-19 pandemic control: balancing detection policy and lockdown intervention under icu sustainability

Multi-population mean field games systems with neumann boundary conditions

Primal-dual extragradient methods for nonlinear nonsmooth pde-constrained optimization

Stability of saddle points via explicit coderivatives of pointwise subdifferentials. Set-Valued and Variational Analysis

Impact of lockdown on covid-19 epidemic inîle-de-france and possible exit strategies

Social contacts and the spread of infectious diseases

A mean field game analysis of sir dynamics with vaccination

The derivation of ergodic mean field game equations for several populations of players

Estimating the effects of non-pharmaceutical interventions on covid-19 in europe

Review of covid-19 vaccine subtypes, efficacy and geographical distributions

A control theory approach to optimal pandemic mitigation

Optimal control of epidemics with limited resources

Large population stochastic dynamic games: closed-loop mckean-vlasov systems and the nash certainty equivalence principle

Solving Large-Scale Optimization Problems with a Convergence Rate Independent of Grid Size

Optimal control problem of an sir reactiondiffusion model with inequality constraints

Beyond just "flattening the curve": Optimal control of epidemics with purely non-pharmaceutical interventions

Constrained optimal control applied to vaccination for influenza

Individual vaccination as nash equilibrium in a sir model with application to the 2009-2010 influenza a (h1n1) epidemic in france

Mean field games

Generalized unnormalized optimal transport and its fast algorithms

Controlling propagation of epidemics via mean-field games

Determination of an optimal control strategy for vaccine administration in covid-19 pandemic treatment. Computer methods and programs in biomedicine

Cold chain transportation decision in the vaccine supply chain

Optimal allocation of face masks during the covid-19 pandemic: a case study of the first epidemic wave in the united states

Analysis of the covid-19 pandemic by sir model and machine learning technics for forecasting

A machine learning framework for solving high-dimensional mean field game and mean field control problems

Adjoint-based data assimilation of an epidemiology model for the Covid-19 pandemic in 2020

Estevão Soares Dos Santos, et al. Optimal control of the covid-19 pandemic: controlled sanitary deconfinement in portugal

Stability analysis and optimal vaccination of an sir epidemic model