Mean field control problems for vaccine distribution.

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 pandemics by the optimal transportation strategy of vaccine distribution. Here, we integrate the vaccine distribution into the mean-field SIR model designed in Lee W, Liu S, Tembine H, Li W, Osher S (2020) Controlling propagation of epidemics via mean-field games. arXiv preprint arXiv:2006.01249. Numerical examples demonstrate that the proposed model provides practical strategies for vaccine distribution in a spatial domain.

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– 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 prioritize the elderly and healthcare workers getting vaccinated. Meanwhile, the shipping of vaccines is 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 to 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:subject toandIn our model, different populations are described using (), representing the susceptible, infectious, and recovered. The term describes the density distribution of the vaccine over the spatial domain at location x and time t. The control variables () create velocity fields over time-space domain that move the corresponding populations. As for vaccines, the control variable 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 and running costs . The terminal costs 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 include the costs of transportation of vaccines and different classes of the populations, etc. We will discuss more details of cost functionals in Sect. 2.3. As for constraints of our optimization problem, the five partial differential equations of , () 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 produced at particular factory locations with a daily maximal production rate . The dynamics of the vaccine density 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, , , 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] which 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 Sect. 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 consideration 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 [12, 13, 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 [32, 35]. 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 [19, 22, 24, 29, 37]. In [22], an optimal vaccine distribution strategy is proposed with a limited total amount of vaccines and maximal daily supply. [29] first uses an 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 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 problem studies the interactions on 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 distribution. Therefore, we obtain a novel mean-field control problem.

The rest of the paper is organized as follows. Section 2 proposes a novel multi-population 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 Sect. 4, we present numerical experiments with different model parameter choices and discuss their implications on mean-field controls.

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.

Classical SIR model

The SIR epidemic model describes an infectious disease epidemic via an ordinary differential equation systemThe 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 homogeneous-mixing, which means that every individual has an equal probability of disease-causing contact. 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 they become infected. This leads to the extension of the SIR model, such as the SEIR model. Our proposed model has a flexible structure and can naturally be generalized to such epidemiological models.

Spatial SIR variational problem with vaccine distribution

In [28], we add the spatial dimension to the S, I, R functions. Let be a bounded domain. Consider the following density functionsHere, , , and represent susceptible, infected, and recovered populations distribution, respectively. We assume for each moves over a spatial domain with a velocity . Here are our controls variables. With change of variables , we define the momentumthat govern the corresponding density flows. In the following, instead of using control variables , we replace them with and regard as the control variables.

We can describe the flows of the densities by the following continuity equations.This system of continuity equations describes the flows of three groups of densities while satisfying the SIR model. The nonnegative constants () are the coefficients for viscosity terms. These terms can also be understood as noise terms generated by the data. is a symmetric positive definite kernel with . In this model, we consider the Gaussian kernelThe kernel convolution describes the spreading rate of 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 , which leads to the following equalityIn this paper, we consider the optimization problem for the distribution of vaccines. We add an extra function which represents the vaccine density in at each time . The vaccine distribution will be described as the following PDE:where is a momentum, represents the utilization rate of vaccines, and represents the production rate of vaccines in at . During , the vaccines are produced with a production rate f and used at a rate . During , the vaccines are delivered to the area where the susceptible population is located, and they are used at a rate of . In summary, the first part of the PDE describes vaccines’ production, and the second part describes the delivery of vaccines. For all time , the susceptible population is 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 and , respectively. The constant represents the vaccine efficiency and represents the vaccinated population at . We denote a set and define a nonlinear operator A as followswhere is a step function that equals 1 on C and 0 otherwise.

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 , satisfies the following criteria: Item (i) can be described byfor wherewhich is convex, lower semi-continuous, and 1-homogeneous with respect to . The parameter characterizes the cost of moving with velocity . Larger means it is more expensive to move . Note that this function comes from the quadratic kinetic energy. To see this, we use the definition and plug into formula (2.5):Item (ii) and (iii) can be described by the terminal costs of the cost functionalwhere functions are convex and lower semi-continuous functions. We also minimize the terminal cost for 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 is far away from 1 for . In this paper, we use a quadratic cost functionwhere is some constant.

minimize the transportation cost for moving each population;

minimize the total number of infected people and the total number of susceptible people by maximizing the usage of the vaccines at time T;

maximize the total number of recovered people at time T;

avoid high concentration of population and vaccines at each time ;

minimize the amount of vaccines produced during ;

minimize the transportation cost for delivering vaccines during .

For Item (iv), the cost functional for the concentration of the total population and vaccines can be represented bywherefor and convex and lower semi-continuous functions . Similar to  (2.6) from Item (ii), we use quadratic functions for and .

Items (v) and (vi) are criteria specific to the vaccine distribution. From PDE (2.2), the vaccines are produced during by a function f. We use the similar functional (2.7) to minimize the amount of vaccines produced by f. Thus, we set the functionalwhere is a convex and lower semi-continuous function.

The vaccines are delivered during . Similar to the Item (i), we setwhere has the same definition as (2.5).

The total cost functional we consider is thenIn 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 , F satisfieswhere is defined asand denotes the convex subdifferential of F. Since , , are convex and lower-semicontinuous, 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).

Constraints for vaccine production

In addition to the constraint from (2.3), we adapt the following constraints to reflect the limited vaccination coverage:where indicates the factory area where vaccines are produced and is a nonnegative constant representing the maximum vaccine production rate. In the third inequality, a nonnegative constant limits the total number of vaccines produced during .Constraints (2.9) can be imposed by having the following functionals for and .where indicates the factory area where vaccines are produced. The functionals and are defined aswhere a, b are constants and is a function. The function is defined asThis function forces if , thus vaccines are produced only in .

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 for the deceased population with a mortality rate .

Properties

From the definition of the cost functional and constraint (2.3), we have the following minimization problem:We first define the inner product of vectors of functions in . Given vectors of functions and with , the inner product of vectors u and v and norm of u are defined bywhere is a inner product such thatWe introduce dual variables 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 is the Lagrangian functional defined asFor brevity, we denoteWe can rewrite the Lagrangian aswhere the nonlinear operator A(u) is defined asAs noted in [28], the dual gap, the difference between the primal solution and dual solution, may not be zero because the nonconvex functions and make the feasible set nonconvex. We circumvent the problem by linearizing the nonlinear operator at a base point In our formulation, the linearized operator can be written as follows.where . We define a linearized Lagrangian asIn the paper [10], the author developed a primal-dual algorithm using the linearized Lagrangian (Algorithm (3.5)) and proves that the sequence from the algorithm converges to the saddle point (in Sect. 3.1, we prove the local convergence to the saddle point given is sufficiently close to the saddle point). By the first-order optimality conditions (also known as Karush–Kuhn–Tucker (KKT) conditions), the saddle point satisfiesIn the next proposition, we present the equations derived from the KKT conditions (2.17).

Proposition 2.2

By KKT conditions, the saddle point of (2.13) satisfies the following equations.The terms , , , , and are the functional derivatives. In other words, given be a smooth functional where is a separable Hilbert space and , we say a map is the functional derivative of F with respect to if it satisfiesfor any arbitrary function .

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.

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 inner product is defined in (2.12) andis the Legendre transform of f. The method solves the saddle point problem by iteratingThe scheme converges if the step sizes and satisfywhere is an operator norm in . 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 algorithm (3.1) with a linear approximation of A at a base point Denote . We have a linearized saddle point problemand the scheme iteratesThe paper [10] proves that the sequence of the algorithm converges to some saddle point that satisfies (2.17). However, the scheme converges if the step sizes satisfySuppose we use an unbounded operator that depends on the grid size, for example, . The discrete approximation of the operator norm of A increases as the grid size increases (Fig. 1 illustrates the relationship between the norm of an unbounded operator and grid sizes). Thus, 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 iterateswhere the norm is defined asBy choosing the proper norms, the step sizes only need to satisfywhich are clearly independent of the grid size.

Fig. 1

The image a shows u on a unit square domain . The plot (b) shows the discrete approximation of with respect to grid sizes. It shows that in the discrete approximation, the norm of an unbounded operator increases as grid size increases

Local convergence of the algorithm

In this section, we show the iterations from algorithm (3.6) locally converges to the saddle point. The local convergence theorem in this paper is mainly based on 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 algorithm (3.6), satisfies the following first-order optimality conditionswhich can be rewritten aswith . Here, the monotone operator is defined asandwhere Id is an identity operator.

Recall that from (2.17), the saddle point has to satisfyThroughout, we assume that

Lemma 3.1

There exists constants and such thatwhere is an operator norm.

Proof

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

Lemma 3.2

Suppose (3.9) holds and let . Then there exist constants such thatwhere

A proof of Lemma 3.2 is provided in the appendix.

With the above Lemmas, we can use Theorem 2.11 from [10] to show the local convergence of the algorithm.

Theorem 3.3

Let be a solution to (2.17) where . Let the step sizes and satisfy for all k. Then there exists such that for any initial point satisfyingthe iterates from (3.6) converges to the saddle point .

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 to satisfy the condition called metric regularity. In our formulation, the constraint makes metrically regular by [11, Section 5.3]. We refer readers to [10, 11, 33] for further details about metric regularity.

Implementation of the algorithm

To implement the algorithm to the minimization problem (2.8), we setWe use (2.15) for the definition of the operator A. Define the Lagrangian functional aswhere is defined in (2.12). We summarize the algorithm as follows.

Here, and norms are defined asfor any . Moreover, the relative error is defined asIn sect. 4, We use quadratic functions for , , , . With definition (2.10), we useThus, we can write the cost functional as followswhere , , , are nonnegative constants. With this cost functional, we find explicit formula for each variable .

Proposition 3.5

The variables (), and from Algorithm 3.1 satisfy the following explicit formulas:where is a positive root of a cubic polynomial and we approximate the as follows

We use FFTW library to compute () and convolution terms by Fast Fourier Transform (FFT), which is operations per iteration with n being the number of points. Thus, the algorithm takes just operations per iteration.

Experiments

In this section, we present several sets of numerical experiments using Algorithm 3.1 with various parameters. We wrote C++ codes to run the numerical experiments. Let be a unit square in and the terminal time . The domain is discretized with the regular Cartesian grid below.where , are the number of discretized points in space and 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 , 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 , the solution will produce the vaccines during and deliver them during . Furthermore, we fix the parameters for the infection rate and recovery rateThe 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 formulation (3.10), we have terminal functionalsThus, 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 (Fig. 2). 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 constraint conditions (2.9).

Fig. 2

Visualization of the flow of three populations. The susceptible transforms to the infected with a rate and the recovered with a rate . The infected transforms to the recovered with a rate

We present two experiments that demonstrate how the various factors in the formulation affect the production and the distribution of vaccines.

Experiment 1

In this experiment, we show that Algorithm 3.1 converges independent of grid sizes when we use the preconditioning operator defined in Proposition 3.5. Consider the initial densities for the () and the factory location aswhere and is a ball of a radius r centered at x. Figure 3 shows the images of initial conditions (4.1).

Fig. 3

Experiment 1: Initial densities of (left) and (right). The green circle indicates

We compute the solution of the SIRV variational problem (2.11) with the above initial conditions using Algorithm 3.1. For simplicity, we assume recovered population density does not move. Thus, we use an arbitrary large number for a parameter to penalize when . The rest of the parameters are identical to the parameters defined in the preceding section. We ran four simulations with same initial conditions and same step sizes (, ) with four different grid sizes:

The result of the experiment is depicted in Fig. 4. The figure shows the convergence plot of the algorithm with respect to the number of iteration for each grid size. The x-axis indicates the iteration number and the y-axis indicates the value of the following Lagrangian functional:Note that this Lagrangian functional is different from  (2.8) and (2.13). The terms with indicator functions , are removed to avoid representing numerically. The absence of the terms may explain that the value increases in the first 500 iterations and then decreases afterward. Figure 5 shows the computed solutions at iteration 3000 from four different spatial grid sizes (, , , ). Each row of the figure shows the evolution of a vaccine density from time to computed from each grid size. These figures clearly show that the algorithm converges to the same saddle point independent of the grid sizes.

Fig. 4

Convergence plot of the algorithm for each grid size () with the same step sizes (, ). The plot shows that the convergence of the algorithm is independent of grid sizes

Fig. 5

Computed solution of a vaccine density variable from Experiment 1. Each row shows the evolution of a vaccine density variable from time to with different grid sizes. Row 1: , Row 2: , Row 3: , Row 4:

Experiment 2

In this experiment, we show how the parameters related to the vaccine density variable () affect the solution. We use the same initial densities for () and f as in Experiment 1. With the initial densities (4.1), we run two simulations with different values for , , and .

Figure 6 shows the comparison between the results from the simulation 1 and the simulation 2. The first three plots (Fig. 6a) show the total mass of (), i.e.,and the last plot (Fig. 6b) shows the total mass of during 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 . Furthermore, the solution from the simulation 1 cannot vaccinate a large number of susceptible due to a small . Thus, there are more susceptible and less recovered at the terminal time in the simulation 1.

Fig. 6

Experiment 2: The comparison between the results from the simulation 1 and the simulation 2. The first three plots a show the total mass of () and the fourth plot b shows the total mass of produced at the factory area during the production time

Experiment 3

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 as obstacles. We use the following functionals in the experiment.The densities () cannot be positive on due to . Thus, the densities transport while avoiding the obstacle . We show two sets of experiments based on this setup.

Single factory

We set the initial densities and as followsand fix the parametersThe initial densities are shown in Fig. 7.

Fig. 7

Experiment 3: The initial densities (left) and (right), and the location of the factory (indicated as a green circle)

Figures 8 and 9 show the evolution of densities with and without obstacles, respectively. In both simulations, the density of vaccines (the fourth row) transports to the areas where the susceptible people are present. In Fig. 9, transports while avoiding the obstacle at the right. Figure 10 shows the comparison between these two solutions and how the presence of the obstacle affects the production and delivery of vaccines quantitatively. Figure 10a shows the total mass of the vaccines in the factory area during the production timeFigure 10b shows the total mass of the vaccines during the delivery time at the left side and the right side of the domainduring . When there is no obstacle, the vaccines are delivered more to the right than to the left (Fig. 10b). The number of susceptible people at the left decreases very fast because there are infected people with a high infection rate. When starts to transport at time , 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, has to bypass the obstacle to reach the susceptible areas. Thus, the kinetic energy cost during the delivery time 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 (Fig. 10a) and delivered to the right during when there is an obstacle (Fig. 10b).

Fig. 8

Experiment 3: The evolution of densities () without the obstacle over time . The first row: the susceptible density . The second row: the infected density . The third row: the recovered density . The fourth row: the vaccine density

Fig. 9

Experiment 3: The evolution of densities () with the obstacle (indicated as a yellow block) over time . The first row: the susceptible density . The second row: the infected density . The third row: the recovered density . The fourth row: the vaccine density

Fig. 10

Experiment 3: The left plot shows the total mass of vaccine density during the production time . The right plot shows the total mass of at the left side of the domain and at the right side of the domain

Multiple factories

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 , and three factory locations in this experiment. We set the initial densities and as followsand fix the parametersThe initial densities are shown in Fig. 11.

Fig. 11

Experiment 3: The initial densities (left) and (right), and the location of the factory (indicated as green circles)

Figures 12 and  13 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 (the fourth row) transports to the areas where the susceptible people are present while avoiding the obstacles.

Fig. 12

Experiment 3: The evolution of densities () without the obstacle over time . The first row: the susceptible density . The second row: the infected density . The third row: the recovered density . The fourth row: the vaccine density

Fig. 13

Experiment 3: The evolution of densities () with the obstacle (colored yellow) over time . The first row: the susceptible density . The second row: the infected density . The third row: the recovered density . The fourth row: the vaccine density

Figure 14a shows the total mass of the vaccines produced during the production time at each factory location. Without the obstacles, the total mass of 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) to reduce the kinetic energy cost during the delivery time . However, the vaccines are produced the most at the middle factory with the obstacles. Since the obstacles block the paths between the top and the bottom factories and the susceptible people, has to bypass them to reach 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, producing more vaccines at the middle factory is more efficient.

Fig. 14

Experiment 3: The top plot shows the total mass of vaccine density at three factory locations during the production time . The bottom plot shows the total mass of at the top left area of the domain , at the bottom left area , at the top right area , and at the bottom right area during the distribution time

Figure 14b shows the total mass of the vaccines during the delivery time at different locations. The lines in the plot represent the following quantities:over . With the obstacles, the kinetic energy cost increases since 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.

Experiment 4

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 are set as followsWe fix the parametersThe initial densities and locations of factories are shown in Fig. 15.

Fig. 15

Experiment 4: The initial densities (left) and (right), the location of the factory (indicated as green circles), and the obstacle (colored yellow)

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 , , . We consider the following PDEs:Furthermore, by taking out the momentum terms from S, I, R groups, the cost functional for this experiment isWith 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 Algorithm 3.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 asFigure 16 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 Fig. 15. The solid lines, the result with the same fixed rates of production, show that all three factories produce identical 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.

Fig. 16

Experiment 4: The plot shows the total mass of vaccine densities during production at each factory location: left, middle, and right. The dotted lines are from the optimal strategy from Algorithm 3.1, and the solid lines are from the fixed production rates

The table below is the quantitative comparison between the two results.

The first row of the table shows that more vaccines are produced with a fixed rate of production. However, the result of 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. The algorithm finds the more efficient strategy with fewer vaccines produced.

N_x_1	N_x_2	N_t
32	32	32
64	64	32
128	128	32
256	256	32

Parameters	Sim 1	Sim 2	Description
θ_1	0.5	0.9	Vaccine efficiency
f_max	0.5	10	Maximum production rate of vaccines
C_factory	0.5	2	Maximum amount of vaccines
			that can be produced at x∈Ω during 0≤t≤12

Quantity	Description	Algorithm 3.1	Fixed rates
∫_Ωρ_V(12,x)dx	The total amount of vaccines produced.	7.997×10^-3	8.411×10^-3
∫_Ωρ_S(1,x)dx	The number of susceptible people at the terminal time.	1.520×10^-2	1.525×10^-2
∫_Ωρ_I(1,x)dx	The number of infected people at the terminal time.	5.133×10^-3	5.134×10^-3
∫121∫_Ω|m_V|^22ρ_Vdxdt	The transportation cost of vaccines.	7.339×10^-3	7.544×10^-3