###### Abstract

The study considers the usage of a probabilistic optimization method called
Cross-Entropy (CE). This is the version of the Monte Carlo method created by
Reuven Rubinstein (1997). It was developed in the context of determining rare
events. Here we will present the way in which the CE method can be used for
problems of optimization of epidemiological models, and more specifically the
optimization of the SIRC (Susceptible - Infectious - Recovered - Cross-immune)
model based on the functions supervising the care of specific groups in the
model. With the help of weighted sampling, an attempt was made to find the
fastest and most accurate version of the algorithm.

###### keywords:

optimal stopping; counting process; cross-entropy method; epidemiological
models; SIR and SIRC models; Cross-immunity and boosting.

xx 1 5 Academic Editor: Georgy Sofronov Received: 22nd of September 2020;
Accepted: date; Published: date Cross-entropy method in application to SIRC
model Maria Katarzyna Stachowiak 1 and Krzysztof Józef Szajowski 2Maria
Katarzyna Stachowiak and Krzysztof Józef Szajowski Correspondence:
mkwstachowiak@gmail.com(MKS<EMAIL_ADDRESS>Tel.: +48-71-3203185(KJS)
Primary: 90C59; Secondary: 49N90; 92D30; 68Q60

##  1 Introduction.

In this study our aim is to develop the possibilities of numerically solving
variational problems that appear in epidemic dynamics models. The algorithm
that we will use for this purpose will be a modified Monte Carlo method. In a
few sentences, we recall those stages of the Monte Carlo method development
and its improvement, which constitute the essence of the applied approach (cf.
Martino et al. (Martino et al., 2018, in sec. 1.1)).

Before constructing the computer, the numerical analysis of deterministic
models used approximations that aimed, among other things, at reducing and
simplifying the calculations. The ability to perform a significantly larger
number of operations in a short time led Ulam (cf. Metropolis Metropolis
(1987)) to the idea of transforming a deterministic task to an appropriate
stochastic task and applying the numerical analysis using the simulation of
random quantities of the equivalent model (cf. Metropolis and Ulam Metropolis
(1987),). The problem Ulam was working on was part of a project headed by von
Neumann, who accepted the idea, and Metropolis gave it the name of the Monte
Carlo Method (Martino et al., 2018, in sec. 1.1).

The resulting idea allowed us to free oneself from difficult calculations,
replacing them with a large number of easier calculations. The problem was the
errors that could not be solved by increasing the number of iterations. Work
on improving the method consisted in reducing the variance of simulated
samples, which led to the development of the Importance Sampling (IS) method
in 1956 (cf. Marshall Marshall (1956) ).

Computational algorithms in this direction are a rapidly growing field.
Currently, the computing power of computers is already large enough, to
provide an opportunity to solve problems that can’t be solved analytically.
However, some problems and methods still require a lot of power. A lot of
research is devoted to finding or improving such methods. The development of
this field is currently very important. There are many models describing
current problems. Not all offer the possibility of an analytical solution.
That is why various studies are appearing to find best methods for accurate
results. the Cross-Entropy method helps to realize IS and it may, among other
things, be an alternative to current methods in problems in epidemiological
models.

###  1.1 Sequential Monte Carlo methods.

The CE method was proposed as an algorithm for estimation of probabilities of
rare events. Modifications of cross entropy to minimize variance were used for
this by Rubinstein in his seminal paper Rubinstein (1997) in 1997. This is
done by translating the "deterministic" optimization problem in the related
"stochastic" optimization and then simulating a rare event. However, to
determine this low probability well, the system would need a large sample and
a long time to simulate.

###  1.2 Importance sampling.

Therefore, it was decided to apply the importance sampling (IS), which is a
technique used for reducing the variance. The system is simulated using other
parameter sets (more precisely, a different probability distribution) that
helps increase the likelihood of this occurrence. Optimal parameters, however,
are very often difficult to obtain. This is where the important advantage of
the CE method comes in handy, which is a simple procedure for estimating these
optimal parameters using an advanced simulation theory. Sani and Kroese Sani
and Kroese (2007) have used the Cross-Entropy method as the support of the
importance sampling to solve the problem of the optimal epidemic intervention
of HIV spread. The idea of his paper will be adopted to treat the SIRC models,
which has a wide application to modeling the division of population to four
groups of society members based on the resistance to infections (v. modeling
of bacterial Meningitis transmission dynamics analyzed by Asamoah et al.
Asamoah et al. (2018), Vereen Vereen (2008)). The result will be compared with
an analytical solution of Casagrandi et al. Casagrandi et al. (2006).

###  1.3 Cross-Entropy method.

This part will contain information on the methods used in the paper. The
Cross-Entropy method developed in Rubinstein (1997) is one of the versions of
the Monte Carlo method developed for problems requiring the estimation of
events with low probabilities. It is an alternative approach to combinatorial
and continuous, multi-step, optimization tasks. In the field of rare event
simulation, this method is used in conjunction with weighted sampling, a known
technique of variance reduction, in which the system is simulated as part of
another set of parameters, called reference parameters, to increase the
likelihood of a rare event. The advantage of the relative entropy method is
that it provides a simple and fast procedure for estimating optimal reference
parameters in IS.

Cross entropy is the term from the information theory (v. Parry (1969), (Cover
and Thomas, 2006, Chapter 2)). It is a consequence of measuring the distance
of the random variables based on the information provided by their
observation. The relative entropy, the close idea to the cross entropy, was
first defined by Kullback and Leibler Kullback and Leibler (1951) as a measure
of the distance $\mathcal{D}(p||q)=\textbf{E}_{p}\log\frac{p(X)}{q(X)}$ (v.
(Cover and Thomas, 2006, (1.6) on p. 9; Definition on Sec. 2.3 p. 19), Inglot
(2014)) for two distributions $p$, $q$ and is known also as Kullback-Leibler
distance. This idea was studied in detail by Csiszár Csiszár (1967) and Amari
Amari (1985). The application of this distance measure in the Monte Carlo
refinements is related to a realization of the importance sampling technique.

The CE algorithm can be seen as a self-learning code covering the following
two iterative phases:

  1. 1.

Generating random data samples (trajectories, vectors, etc.) according to a
specific random mechanism.

  2. 2.

Updating parameters of the random mechanism based on data to obtain a “better”
sample in the next iteration.

Now the main general algorithm behind the Cross-Entropy method will be
presented on the examples.

###  1.4 Application to optimization problems.

This section is intended to show how extensive the CE method is. The examples
will show the transformation of the optimization problems to the Monte Carlo
task of estimation for some expected values. The equivalent MC task uses the
CE method. Other very good examples can be found in Rubinstein (1999).

The first variational problem presented in this context is the specific
multiple stopping model. Consider the so-called secretary problem (v. Ferguson
Ferguson, T.S. (1989) or the second author’s paper Szajowski, K. (1982)) for
multiple choice, i.e. the issue of selecting the best proposals from a finite
set with at most $k$ attempts. There is a set with $N$ objects numbered from
$1$ to $N$. By convention, the object with the number 1 is classified as the
best and the object with the number $N$ as the worst. Objects arrive one at a
time in a random order. We can accept this object or reject it. However, if we
reject an object, we cannot return. The task is to find such an optimal
stopping rule that the sum of the ranks of all selected objects is the lowest
(their value is then the highest). So we strive for a designation

$\inf_{\tau}\mathbf{E}X_{\tau}=\mathbf{E}X_{\tau^{*}}.$

Let $a_{i}$ be the rank of the selected object. The goal is to find a routine
for which the value of $\mathbf{E}(a_{\tau_{1}}+\dots+a_{\tau_{k}})$ for
$k\geqslant 2$, is the smallest. More details are moved to A.1.

Next, the problem was transformed to the minimization of the mean of sum of
ranks. The details are presented in the section A.2. This problem and its
solution by CE was described by Polushina Polushina (2010). Its correctness
was checked and programmed in a different environment by Stachowiak in
Stachowiak (2019).

The second example is presented in details in Section A.3. It is a formulation
of the vehicle routing problem (v. Dror (2002), Chepuri and Homem-de Mello
(2005)). This is an example of stochastic optimization where the Cross-Entropy
is used.

###  1.5 Goal and organization of this paper.

In the following parts of the paper, we will focus on models of the dynamics
of the spread of infection over time when the population consists of
individuals susceptible to infection, sick, immune to vaccination or past
infection, and partially susceptible i.e. the population is assigned to four
compartments of individuals. Actions taken and their impact on the dynamics of
the population are important. It is precisely the analysis of the impact of
the preventive diagnosis and treatment that is particularly interesting in the
model. How the mathematical model covers such actions is presented in Section
2.1. The model created in this way is then adapted for Monte Carlo analysis in
2.2 and the results obtained on this basis are found in 2.3. The analysis of
computational experiments concludes this discussion in 3.

##  2  Optimization of control for SIRC model.

Let us focus the attention on the introduction of SIRC model and presentation
of the logic behind using the CE method to determine functions that optimize
the spread of epidemics. The ideology behind creation of the SIRC model and
its interpretation will be presented in order to match the right parameters
for calculating the cost functions.

The CE method here is used to solve the variational deterministic problem.
Solving such problems has long been undertaken with the help of numerical
methods with a positive result. In the position Ekeland and Temam (1999)
proposals of the route for variational problems for optimization problems in
time are presented. The main focus was on the non-convex problems, which often
occur with optimal control problems. In the book by Glowinski Glowinski (2008)
a review of the methods for solving variational problems was made. And then
Mumford and Shah (1989) presented a complicated method of approximating
functions for the problem optimization.

###  2.1 SIRC model.

The subject of the work is to solve the problem optimizing the spread of the
disease, which can be modeled with the SIRC model. This model was proposed by
Casagrandi et al.Casagrandi et al. (2006). Its creation was intended to create
a better model that would describe the course of influenza type A. Contrary to
appearances, this is an important topic because the disease, although it is
widely regarded as weak, is a huge problem for healthcare. In the US alone,
the cost associated with the influenza epidemic in the season is estimated to
exceed USD 10 billion (v. Klimov et al. (1999)) and the number of deaths is
over 21,000 per season (v. Simonsen et al. (1997)).

Various articles have previously proposed many different mathematical models
to describe a pandemic of influenza type A. An overview of such models was
made by Earn Earn et al. (2002). In general, it came down to combining SIR
models with the use of cross-resistance parameters (cf. details in the papers
by Andreasen et al. Andreasen et al. (1997) and Lin et al. Lin et al. (1999)).
The authors of the article write that the main disadvantage of this approach
is the difficulty of the analysis and calculations with a large number of
strains. The results of their research showed that the classic SIR model
cannot be used to model and study influenza epidemics. Instead, they proposed
the extension of the SIR model to a new class C that will simulate a state
between a susceptible state and a fully protected state. This extension aims
to cover situations where vulnerable carriers are exposed to similar stresses
as they had before. As a result, their immune system is stronger when fighting
this disease.

The SIRC model divides the community into 4 groups:

  * $\sqbullet$

S - persons susceptible to infection, who have not previously had contact with
this disease and have no immune defense against this strain

  * $\sqbullet$

I - people infected with the current disease

  * $\sqbullet$

R - people who have had this disease and are completely immune to this strain

  * $\sqbullet$

C - people partially resistant to the current strain (e.g. vaccinated or those
who have had a different strain)

Figure 1 contains a general scheme of this model. The four rhombus represent
the four compartments of individuals, the movement between the compartments is
indicated by the continuous arrows. The main advantage of this model over
other SIR users is the fact that after recovering from a given strain, in
addition to being completely resistant to this strain, they are also partly
immune to a new virus that will appear later. This allows you to model the
resistance and response of people in the group to different types of disease.

Figure 1: Schematic representation of the SIRC model.

Parameters $\alpha,\delta$ i $\gamma$ can be interpreted as the reciprocal of
the average time spent by a person in order of ranges $I,R,C$. Parameters
$\mu$ with indices $S,I,R,C$ represent the natural mortality in each group,
respectively. In some versions of the model, additional mortality rates are
considered for the group of infected people. Here we assume that this factor
is not affected by the disease and we will denote it as $\mu$ in later
formulas ($\mu_{S}=\mu_{I}=\mu_{R}=\mu_{C}=\mu$). The next one $\sigma$ is the
likelihood of reinfection of a person who has cross-resistance while the
parameter $\beta$ describes the contact indicator.

The SIRC model is represented as a set of four ordinary differential
equations. Let $S(t)$, $I(t)$ be the number of people in adequate
compartments. We have (v. Casagrandi et al. (2006))

$\displaystyle\frac{dS}{dt}$ $\displaystyle=\mu(1-S(t))-\beta S(t)I(t)+\gamma
C(t),$ (1) $\displaystyle\frac{dI}{dt}$ $\displaystyle=\beta
S(t)I(t)+\sigma\beta C(t)I(t)-(\mu+\alpha)I(t),$ $\displaystyle\frac{dR}{dt}$
$\displaystyle=(1-\sigma)\beta C(t)I(t)+\alpha I(t)-(\mu+\delta)R(t),$
$\displaystyle\frac{dC}{dt}$ $\displaystyle=\delta R(t)-\beta
C(t)I(t)-(\mu+\gamma)C(t),$

with initial conditions

$\begin{split}S(0)=S_{0},\,I(0)=I_{0},\,R(0)=R_{0},\,C(0)=C_{0}\end{split}$
(2)

###  2.2 Derivation of optimization functions

An approach to optimize this model was proposed by Iacoviello and Stasio
Iacoviello and Stasio (2013). In this article one can find the suggestion
concerning performing the calculations as outlined by Zaman et al. in Zaman et
al. (2008). Two functions were proposed in the approach as the parameter or
controls. One relates to the _susceptible people_ and the other describes the
number of _sick people_. The method described by Kamien and Schwartz Kamien
and Schwartz (1991) was used to determine the optimal solution. In order to
apply it, the set of equations (1) should be updated accordingly. After these
modifications, it has the following form:

$\displaystyle\frac{dS}{dt}$ $\displaystyle=\mu(1-S)-\beta SI+\gamma
C-g(S(t),u(t)),$ (3) $\displaystyle\frac{dI}{dt}$ $\displaystyle=\beta
S(t)I(t)+\sigma\beta C(t)I(t)-(\mu+\alpha)I(t)-h(I(t),v(t)),$
$\displaystyle\frac{dR}{dt}$ $\displaystyle=(1-\sigma)\beta C(t)I(t)+\alpha
I(t)-(\mu+\delta)R(t)+g(S(t),u(t))+h(I(t),v(t)),$ $\displaystyle\frac{dC}{dt}$
$\displaystyle=\delta R(t)-\beta C(t)I(t)-(\mu+\gamma)C(t),$

where

$\displaystyle g(S(t),u(t))$ $\displaystyle=\rho_{1}S(t)u(t),$ (4)
$\displaystyle h(I(t),v(t))$ $\displaystyle=\rho_{2}I(t)v(t).$

$u(t)$ represent the percentage of those susceptible who have been taken care
of thanks to using control on population and $v(t)$ represent the percentage
of infected people with the same description, respectively. $\rho_{1}$ and
$\rho_{2}$ are weights that optimize the proportion of given control options.
Both functions $g(S(t),u(t))$ $h(I(t),v(t))$ and in order they can be
interpreted as actions performed on people susceptible to the disease and
infected. In addition, two new conditions related to optimization functions
are added to the initial conditions.

$\begin{split}u_{min}\leqslant u(t)\leqslant u_{max},\,v_{min}\leqslant
v(t)\leqslant v_{max}\end{split}$ (5)

They describe the limits on what part of the population care can be given to.
The smallest values $u_{min}$ and $v_{min}$ is zero. The existence of the
solution is shown by Iacoviello et al. Iacoviello and Stasio (2013) and they
are functions that minimize the following cost index:

$J(u,v)=\int_{t_{1}}^{t_{2}}\varphi(S(t),I(t),u(t),v(t))dt=\int_{t_{1}}^{t_{2}}\left[\alpha_{1}S(t)+\alpha_{2}I(t)+\frac{1}{2}\tau_{1}u^{2}(t)+\frac{1}{2}\tau_{2}v^{2}(t)\right]dt$
(6)

$\alpha_{1},\alpha_{2}$ are used to maintain the balance in the susceptible
and infected group. $\tau_{1},\tau_{2}$ is interpreted as weighting in the
cost index. $t_{1},t_{2}$ determine the time interval. The square next to the
functions $u(t)$ and $v(t)$ indicates increasing intensity of functions Zaman
et al. (2008). Thus the function $J(u,v)$ reflects the human value of
susceptible and infected persons, taking into account the growing value of
funds used over a specified period of time.

The objective function can be changed as long as there is a minimum. The CE
method does not impose any restrictions here. However, to compare the results
with Iacoviello and Stasio (2013), the function proposed by them is used also
here. There, the objective function must be square in order for the quadratic
programming techniques to be used. Usually quadratic function is good in
mathematics, but not in practice. Changing the objective function and impact
of it can be considered further.

###  2.3 Optimization of the epistemological model by CE method.

Sani and Kroese Sani and Kroese (2007) present the way in which the Cross-
Entropy method can be used to solve the optimization of the epistemological
model consisting of ordinary differential equations. The main task was to
minimize the objective function $J(u)$ depending on one optimization function
$u(t)$ over a certain set $U$ consisting of continuous functions $u$. The
minimum can be saved as:

$\gamma^{*}=J(u^{*})=\min_{u\in U}J(u)$ (7)

Parameterizing the minimum problem looks as follows:

$\min_{c\in C}J(u_{c})$ (8)

where $u_{c}$ is a function from $U$, which is parameterized by a certain
control vector $c\in\mathbb{R}^{m}$. Collection $C$ is a set of vectors $c$.
It should be noted that the selected set $C$ should be big enough to get the
enough precise solution. Then in such a set there are such $c^{*}$– which will
be the optimal control vector and $\gamma^{*}$–optimal value, respectively.

One way to parameterize a problem is to divide the time interval
$[t_{1},t_{2}]$ at small intervals and using these intervals together with
control vector points $c$ to define a function $u_{c}$. This function can be
created by interpolating between points. Such interpolation can be done using
e.g. the finite element method, finite difference method, finite volume method
or cubic B-Spline . Fang et al. in Fang et al. (2002) and Caglar et al. in
Caglar et al. (2006) show that all of this methods can be used to the two-
point boundary value problems. Here the FEM method was used (v. Papadopoulos
and Giovanis (2018)). The results for the other methods have not been checked.
The choice of the approximation method is a difficult and interesting issue,
but this is not what is consider here. With this function $u_{c}$, it takes
form

$u_{c}(t)=\sum_{t=0}^{n}c_{i}k_{i}(t)$ (9)

where

$k_{i}(t)=\begin{cases}\frac{t-t_{i-1}}{t_{i}-t_{i-1}}&\text{for
$t\in[t_{i-1},t_{i}]$}\\\ \frac{t_{i+1}-t}{t_{i+1}-t_{i}}&\text{for
$t\in[t_{i},t_{i+1}]$}\\\ 0&\text{otherwise.}\end{cases}$ (10)

Now when $u_{c}$ is represented as a parameterized function, the value can be
counted$J(u_{c})$ by solving the system of ordinary differential equations (3)
e.g. by the Runge-Kutta method. The idea of using the CE method in this
problem is to use a multi-level algorithm with generated checkpoints $c_{i}$,
here, from a normal distribution. The CE method does not impose the
distribution. Usually it is chosen from the family of density which is
expected. Here, it’s expected that there will be one high peak at the peak of
the epidemic. Hence the normal distribution. Then, update the distribution
parameters using the CE procedure using the target indicator as a condition
for selecting the best samples. The calculations are carried out until the
empirical variance of optimal control function is smaller than the given
$\epsilon$. This means that the function values are close to the optimal
expected value The idea of applying the CE method to the SIRC model is
similar. The only difficulty is the fact that when optimizing the SIRC model
there are two functions instead of one as in the case of the algorithm
described here. In the form of an algorithm, it looks as follows:
Algorithm ( _Modification of theSani and Kroese’s algorthm (v. Sani and Kroese
(2007)) to two optimal functions case_):

  1. 1.

Initialize: Choose $\mu^{(u)}_{0}=\\{\mu_{i0},i=0,\dots,n\\}$,
$\sigma^{(u)}_{0}=\\{\sigma_{i0},i=0,\dots,n\\}$,
$\mu^{(v)}_{0}=\\{\mu_{i0},i=0,\dots,n\\}$ and
$\sigma^{(v)}_{0}=\\{\sigma_{i0},i=0,\dots,n\\}$, set $k=1$

  2. 2.

Draw: Generate a random samples $C_{1},\dots,C_{N}\sim
N(\mu^{(u)}_{k-1},{\sigma^{(u)}_{k-1}}^{2})$ and $D_{1},\dots,D_{N}\sim
N(\mu^{(v)}_{k-1},{\sigma^{(v)}_{k-1}}^{2})$ with
$C_{m}=\\{C_{mi},i=0,\dots,n\\}$ and $D_{l}=\\{D_{li},i=0,\dots,n\\}$.

  3. 3.

Evaluate: For each control vector $C_{m}$ and $D_{m}$ evaluate the objective
function $J(u_{C_{m}},u_{D_{l}})$ by solving the ODE system

  4. 4.

Find the $p$ best performing samples, based on the values
$\\{J(u_{C_{m},D_{l}})\\}$. Let $I$ be the corresponding set of indices.

  5. 5.

Update: for all $i=0,\dots,n$ let

$\displaystyle\hat{\mu_{ki}}^{(u)}=\frac{1}{p}\sum_{m\in
I}C_{mi};\qquad(\hat{\sigma_{ki}}^{(u)})^{2}=\frac{1}{p}\sum_{m\in
I}(C_{mi}-\mu^{(u)}_{ki}),$
$\displaystyle\hat{\mu_{ki}}^{(v)}=\frac{1}{p}\sum_{l\in
I}D_{li};\qquad(\hat{\sigma_{ki}}^{(v)})^{2}=\frac{1}{p}\sum_{l\in
I}(D_{li}-\mu^{(v)}_{ki}).$

  6. 6.

Smooth: For a fixed smoothing parameter $0<\alpha\leqslant 1$ let

$\displaystyle\hat{\mu_{k}}^{(u)}$
$\displaystyle=\alpha\mu_{k}^{(u)}+(1-\alpha)\hat{\mu_{k-1}}^{(u)};\qquad\hat{\sigma_{k}}^{(u)}$
$\displaystyle=\alpha\sigma_{k}^{(u)}+(1-\alpha)\hat{\sigma_{k-1}}^{(u)}$
$\displaystyle\hat{\mu_{k}}^{(v)}$
$\displaystyle=\alpha\mu_{k}^{(v)}+(1-\alpha)\hat{\mu_{k-1}}^{(v)};\qquad\hat{\sigma_{k}}^{(v)}$
$\displaystyle=\alpha\sigma_{k}^{(v)}+(1-\alpha)\hat{\sigma_{k-1}}^{(v)}$

  7. 7.

Repeat 2–6 until $max_{i}\sigma_{ki}^{(u)}<\epsilon$ and
$max_{i}\sigma_{ki}^{(v)}<\epsilon$ with $\epsilon=10^{-5}$ . Let $L$ be the
final iteration number. Return $\mu_{L}$ as an estimate of the optimal control
parameter $c^{*}$

The value of the criterion $J(\cdot,\cdot)$ is calculated in steps 3 and 4.
The selection of the function $\varphi$ allows to fit the cost function to the
modeled case. Here, the criterion given by (6) is adopted.

##  3 Description of the numerical results.

This section will present the results obtained using the Cross-Entropy method.
For the correct comparison the same parameters for the SIRC model were used.
These values are:

$\displaystyle\mu$ $\displaystyle=\frac{1}{75}[\text{year}^{-1}],\,$
$\displaystyle\gamma$ $\displaystyle=\frac{1}{2}[\text{year}^{-1}],\,$
$\displaystyle\alpha$ $\displaystyle=\frac{365}{5}[\text{year}^{-1}],$
$\displaystyle\delta$ $\displaystyle=1[\text{year}^{-1}],\,$
$\displaystyle\beta$ $\displaystyle\approx 146[\text{year}^{-1}],\,$
$\displaystyle\sigma$ $\displaystyle\approx 0.078.$

For optimization functions $u(t)$ and $v(t)$ the following restrictions have
been applied:

$0\leqslant u(t)\leqslant 0.9,\,0\leqslant v(t)\leqslant 0.9$

Restrictions were proposed by Lenhart and Workman Lenhart and Workman (2007)).
Their value was explained by the fact that the whole group cannot be
controlled. In addition, the weights used for functions have values
$\rho_{1}=2$ and $\rho_{2}=2$
Now, all that’s left is to propose variable values for the objective function:

$\alpha_{1}=10^{-3},\,\alpha_{2}=0.997,\,\tau_{1}=10^{-3},\,\tau_{2}=10^{-3}.$

The parameter values are adjusted using historical data like parameter
mortality rate $\mu$, which is counted as average lifetime of a host or
parameters $\alpha,\delta,\gamma$ mean the inverted time of belonging to the
group $I,R$ and $C$. A detailed description of determining parameters along
with their limitations is presented in the article by Casagrandi et al.
Casagrandi et al. (2006).

###  3.1 A remark about adjusting the parameters of the control determination
procedure.

Here these parameters are given and describes the influenza A epidemics well,
because of long historical data sets. However, this may not be the case, and
then these parameters are subject to adjustment in the control determination
procedure. An overview of applications and accuracy of calibration methods,
which can be used for it is presented in the article by Hazelbag et al.
Hazelbag et al. (2020), which provides an overview of the model calibration
methods. All parameters in the model are constant and independent of time, so
methods which try to optimizes a goodness-of-fit (GOF) can be use to this
example. GOF is a measure that assesses consistency between model the output
and goals. As a result, it gives the best combination of parameters. Examples
of such methods are Grid search, Nelder-Mead method, the Iterative, descent-
guided optimisation algorithm, Sampling from a tolerable range. After finding
the appropriate parameters, other algorithms like the profile likelihood
method or Fisher information can be used to calculate the confidence intervals
for these coefficients. If the epidemic described by the model consists of
transition probabilities that cannot be estimated from currently available
data, calibrations can be performed to many end points. Then GOF is measured
as the mean percentage deviation of the values obtained at the endpoints (v.
Taylor et al. (2010)).

###  3.2 Proposed optimization methods in the model analysis.

The results will be presented for two moments in the model: at the beginning
and in the middle of the epidemic. This is initialized with other initial
parameters. For the beginning of the epidemic, they are as follows:

$S_{0}=1-I_{0},\,I_{0}=10^{-6},\,R_{0}=0,\,C_{0}=0.$

And for the widespread epidemic:

$S_{0}=0.99,\,I_{0}=5*10^{-3},\,R_{0}=3*10^{-3},\,C_{0}=2*10^{-3}.$

The values in the model have been normalized for the entire population, and
they add up to 1. The results for both models without the use of controls are
shown in 2. They were obtained using the Runge-Kutty method. The cost index in
this case was in order $0.00799$ and $0.00789$ for the model of the beginning
and the development of an epidemic.

In Sections 3.3 and 3.4 two ways of calculating optimal functions $u(t)$ and
$v(t)$ for the SIRC model will be considered. In the first version, the
functions will be calculated separately. In the next both $u(t)$ and $v(t)$
will be counted at the same time. All versions will be presented at two
moments of the epidemic: when it began and when it has already spread. The
results will be compared with those obtained in the article by Iacoviello et
al. Iacoviello and Stasio (2013), where the problem was solved using the
sequential quadratic programming method using a tool from Matlab.
Unfortunately, the article does not specify the cost index value for solving
the obtained sequential quadratic programming method. Therefore, in the paper,
there were attempts to recreate the form of the function $u(t)$ and $v(t)$ and
the cost index obtained for them ($0.003308$ for the first situation and
$0.006489$ for the second one) and the model solution.

Figure 2: Solutions for SIRC models without optimal control (in order the
epidemic started and developed)

As can be seen in Figure 3 in the first situation, the control functions
helped reduce peak infection, which was previously without control (Figure 2)
occurred between 2 and 4 months. Previously, the largest peak was around
$20\%$ occurred between 2 and 4 months. Previously, the largest peak was
around $10\%$. In the case of the second situation, the peak also decreased,
but not very clearly (it is also visible in the value of the cost index, which
for the first situation is $0.003308$, and for the other one $0.006489$). It
can be seen that timely control is also very important to reduce the harmful
effects of an epidemic. For more conclusions on how to apply control properly
in this model, see Iacoviello and Stasio (2013).

Figure 3: Values of control functions and solutions obtained for them. Solved
by the sequential quadratic programming method

###  3.3 CE method version 1.

In the first method, the algorithm shown at the end of the section 2.3 has
been applied twice. Once for the calculation of the optimal function $u(t)$
with a fixed $v(t)$, and then $v(t)$ with a fixed $u(t)$. In the next step,
the results were combined and the result for the SIRC model was calculated
using these two functions. This approach is possible due to the lack of a
combined restriction on $u(t)$ and $v(t)$. Features are not dependent on each
other with resources.

The cost index was $0.003086$ for the start of epidemic and $0.006443$ a
developed epidemic. The figure 4 show the graphs of functions $u(t)$ and
$v(t)$ along with a graph of solutions for both moments during the epidemic.

Figure 4: Values of control functions and solutions obtained for them. Solved
by the cross-entropy method - version 1

For functions $u(t)$ the results are similar, only the decrease here is more
linear. However, the solution of the SIRC model gives the same results.
Function $v(t)$ is much more interesting. The CE method showed that it should
have very low values, even close to $0$ throughout the entire period in both
cases. However, this does not affect the cost index (it is even smaller than
in the previous case) and the SIRC model. More on this subject will be
included in the summary.

###  3.4 CE method version 2.

In the second case, it was necessary to change slightly the algorithm from
Section 2.3 to be able to use it for two control functions. In point 1 of the
algorithm a variable is generated $C_{m}\sim N(\mu_{k},\sigma_{k}^{2})$, which
serves to designate $u_{C_{m}}$, where $\mu_{k},\sigma_{k}$ are requested by
the user. Here 2 $\mu$ will need to be initialized
$\mu_{k}^{(u)},\mu_{k}^{(v)}$ and 2 $\sigma$:
$\sigma_{k}^{(u)},\sigma_{k}^{(v)}$. Then designate
$C_{1}^{(u)},\dots,C_{N}^{(u)}\sim N(\mu_{k}^{(u)},\sigma_{k}^{(u)2})$ and
similarly $C_{1}^{(v)},\dots,C_{N}^{(v)}\sim
N(\mu_{k}^{(v)},\sigma_{k}^{(v)2})$. Now $J(u,v)$ can be count without
problems. Further the algorithm goes the same way, and in point 6 all 4
parameters $\mu_{k}^{(u)},\mu_{k}^{(v)},\sigma_{k}^{(u)},\sigma_{k}^{(v)}$ are
updated. The graphs obtained for these functions are below.

Figure 5: Values of control functions and solutions obtained for them. Solved
by the cross-entropy method - version 2

In this case, we have very similar results as in the previous section. This
can be seen in both the charts and the cost index values ($0.003044$ and
$0.006451$). The only problem here is the need for a large sample to bring the
results closer.

###  3.5 Comparison of results.

This section gathers the results from each version and compares them with each
other. At first, the behavior of each group of the SIRC model was looked at
and compared with the solution of the uncontrolled model. This is shown in
Figures 6, 7, 8, 9. The results are presented only for the first situation,
where the epidemic begins to spread, because there you can see the results of
using control much better.

Figure 6: Comparison of results obtained for group S among various methods
Figure 7: Comparison of results obtained for group I among various methods

The most interesting is Figure 7 comparing infected people. You can see here
very well how functions $u(t)$ and $v(t)$ help in the fight against the
epidemic. Let two versions of the CE method overlap. The functions obtained by
the sequential quadratic programming method give slightly better results.
However, they have a worse cost index (the comparison is presented in Table
1). This is probably due to the high value of the function $v(t)$, which
increases the cost index. The question remains why there is such a discrepancy
between the obtained function values $v(t)$ with both methods. In this article
Iacoviello and Stasio (2013) weight and impact of functions $u(t)$ and $v(t)$
were compared. The chart in the article compares the results obtained for the
application of both functions, one at a time and none. It turned out that the
main influence on the number of infected has control over susceptible persons,
i.e. the function $u(t)$. Choosing the optimal values for $v(t)$ is not very
important for the solution of the model. Therefore, such discrepancies are
possible when choosing optimal values.

Figure 8: Comparison of results obtained for group R among various methods
Figure 9: Comparison of results obtained for group C among various methods

The table 1 contains a comparison of the cost index values and the duration of
the algorithms for each method. The first value is when the epidemic started,
and the second is when it has already developed. The article does not specify
the algorithm time of the method used.

Method | Cost index | Time
---|---|---
Results without control functions | 0.00799 | -
0.00789
Sequential quadratic programming method | 0.003308 | -
0.006489
Cross-entropy method version 1 | 0.003086 | 30s
0.006443
Cross-entropy method version 2 | 0.003044 | 4m54s
0.006451
Table 1: Summary of cost indexes and times for all situations

It can be seen that for each version the cost index values came out very
similar. The differences can be seen in time. This is because in the first
version of the CE method a smaller sample was simulated because both functions
were considered separately. Fortunately, the functions were not closely
related and it was possible to consider them separately.

##  4 Summary.

The study examined the possibility of using the CE method to determine the
optimal control in the SIRC model. Similar models are being developed to model
the propagation of malicious software in computer networks (v. Taynitskiy et
al. Taynitskiy et al. (2017)). However, the number of works with this approach
is still relatively small, although its universality should encourage checking
its effectiveness. There are two control functions in the model, which when
properly selected optimize the cost function. The CE method was used in two
versions: considering the two functions separately and together. The results
were then compared with those obtained by Iacoviello and Stasio Iacoviello and
Stasio (2013) using routing of Matlab111Matlab is a paid environment with a
number of functions. Unlike other environments and programming languages,
additional features are created by developers in a closed environment and can
only be obtained through purchase.. The Cross-Entropy method had little
problems when considering two functions at the same time. Due to the large
number of possibilities, it was necessary to simulate a large sample, which
significantly extended the algorithm’s time. The second way came out much
better. However, it only worked because the functions were not dependent on
each other. In the opposite situation there is also a possibility to get
result via CE results however, it will be necessary to change the algorithm
and the dependence of how the functions depend on each other and add this to
the algorithm when initiating optimization functions in the code. The same
situation applies to changing the objective function and the density used for
approximation. Changing the objective function gives the possibility to use
other methods. If it stays quadratic and the constraints are linear, quadratic
programming techniques like the sequential quadratic programming method can
still be used. If the objective function and the constraints stays convex, we
use general methods from convex optimization. The CE method is based on well-
known and simply classical rules and is therefore quite problem-independent.
There are no rigid formulas and therefore it requires consideration for each
problem separately. So it is very possible to reconsider it in another way,
which can be interesting for further work. The CE method may also facilitate
the analysis of modifications of epidemiological models related to virus
mutation or delayed expression (v. Gubar et al. Gubar et al. (2018),
Kochańczyk et al. Kochańczyk et al. (2020)).

KSz gives an idea of application CE to SIRC model; MKS realizes implementation
and simulation.

KSz are thankful to Faculty of Pure and Applied Mathematics for research
funding.

The authors declare no conflict of interest.

The following abbreviations are used in this article:
CE Cross Entropy Method (v. page 1.2) FEM the finite element methods (v. page
2.3) GOF a goodness-of-fit (v. page 3.1) IS the importance sampling (v. page
1.2) SIR The SIR model is one of the simplest _compartmental_ models. The
letters means the number of Susceptible-Infected-Removed (and immune) or
deceased individuals. (v. Kermack and McKendrick (1991), Harko et al. (2014)).
SIRC The SIR model with the additional group of partially resistant to the
current strain people: Susceptible - Infectious - Recovered - Cross-immune (v.
page 2.1).

yes

##  Appendix A Optimization by the method of cross-entropy.

###  A.1 Multiple selection to minimize the sum of ranks.

Let $(a_{1},\dots,a_{N})$ be a permutation of integers $(1,2,\dots,N)$, where
all permutations are equally probable. For every $i=1,2,\dots,N$ $Y_{i}$ will
be the number of values $a_{1},\dots,a_{i}$, which are $\leqslant a_{i}$.
$Y_{i}=\text{card}\\{1\leqslant j\leqslant i:a_{j}\leqslant i\\}$ is the
relative rank of the $i$-th object. The decision-maker observe the relative
ranks and the realite ranks define filtration:
$\mathcal{F}_{i}=\sigma\\{Y_{1},Y_{2},\ldots,Y_{i}\\}$. Let $\mathfrak{S}$ be
the set of Markov time with respect to $\\{\mathcal{F}_{i}\\}_{i=1}^{N}$.
Define

$\displaystyle v$
$\displaystyle=\inf_{\tau}\mathbf{E}(a_{\tau_{1}}+\dots+a_{\tau_{k}}),\qquad\text{where
$\tau=(\tau_{1},\dots,\tau_{k})$, $\tau_{i}\in\mathfrak{S}$.}$

We want to find the optimal procedure
$\tau^{*}=(\tau_{1}^{*},\dots,\tau_{k}^{*})$ and the value of the optimization
problem $v$. Let $\mathcal{F}_{(m)_{i}}$ be $\sigma$-algebra generated by
$(Y_{1},\dots,Y_{m_{i}})$. If we accept

$\displaystyle Z_{(m)_{k}}$
$\displaystyle=\mathbf{E}(a_{\tau_{1}}+\dots+a_{\tau_{k}}|F_{(m)_{k}}),$ then
$\displaystyle v$
$\displaystyle=\inf_{\tau}\mathbf{E}Z_{\tau},\qquad\tau=(\tau_{1},\dots,\tau_{k}).$

The problem is reduced to the problem of stopping sequences multiple times
$Z_{(m)_{k}}$. As shown in Nikolaev, M. L. (1998) and Nikolaev (1998) the
solution to the problem is the following strategy:
If there are integer vectors

$\displaystyle\delta^{(k)}=(\delta^{(k)}_{1},\dots,\delta^{(k)}_{N-k+1}),$
$\displaystyle
0\leqslant\delta^{(k)}_{1}\leqslant\dots\leqslant\delta^{(k)}_{N-k}<\delta^{(k)}_{N-k+1}=N,$
(11)
$\displaystyle\delta^{(2)}=(\delta^{(2)}_{k-1},\dots,\delta^{(2)}_{N-1}),$
$\displaystyle
0\leqslant\delta^{(2)}_{k-1}\leqslant\dots\leqslant\delta^{(2)}_{N-2}<\delta^{(2)}_{N-1}=N,$
$\displaystyle\delta^{(1)}=(\delta^{(1)}_{k},\dots,\delta^{(1)}_{N}),$
$\displaystyle
0\leqslant\delta^{(1)}_{k}\leqslant\dots\leqslant\delta^{(1)}_{N-1}<\delta^{(1)}_{N}=N,$
$\displaystyle\delta^{i_{1}}_{j}\leqslant\delta^{i_{2}}_{j},\qquad 1\leqslant
i_{1}\leqslant i_{2}\leqslant k,$ $\displaystyle k-i_{1}+1\leqslant j\leqslant
N-i_{2}+1$

then

$\begin{split}&\tau_{1}^{*}=min\\{m_{1}:y_{m_{1}}\leqslant\delta^{(k)}_{m_{1}}\\},\\\
&\tau_{i}^{*}=min\\{m_{i}>m_{i-1}:y_{m_{i}}\leqslant\delta^{(k-i+1)}_{m_{i}}\\},\end{split}$

on the set
$F_{i-1}=\\{\omega:\tau_{1}^{*}=m_{1},\dots,\tau_{i-1}^{*}=m_{i-1}\\}$.
$i=2,\dots,k$, $F_{0}=\Omega$. For small $N$ there is a possibility to get
accurate values by the analytical method. The CE method solves this problem by
changing the estimation problem to the optimization problem. Then by
randomizing this problem using a defined family of probability density
functions. With this the CE method solves this efficiently by making adaptive
changes to this pdf and go in the direction of the theoretically optimal
density.

###  A.2 Minimization of mean sum of ranks.

Let’s consider the following problem of minimizing the mean sum of ranks:

$\min_{x\in\chi}\mathbf{E}S(x,R),$ (12)

where $\chi=\\{x=(x^{(1)},\dots,x^{(k)}):$ conditions from (11) are
preserved$\\}$ is some defined set. $R=(R_{1},\dots,R_{N})$ is a random
permutation of numbers $1,\dots,N$. $\hat{S}$ is an unbiased estimator
$\mathbf{E}S(x,R)$ with the following formula:

$\hat{S}(x)=\frac{1}{N_{1}}\sum_{n=1}^{N_{1}}(R_{n\tau_{1}}+\dots+R_{n\tau_{k}}),$

where $(R_{n1},\dots,R_{nN})$ is the nth copy of a random permutation $R$.
Now a cross-entropy algorithm can be applied (see 1.3, Rubinstein (1997)).
Let’s define the indicator collections $\\{I_{\\{S(x)\leqslant\gamma\\}}\\}$
on $\chi$ for different levels $\gamma\in R$. Let $\\{f((\cdotp;u)\\}$ be the
density family on $\chi$ parameterized by the actual parameter value $u$. For
a specific $u$ we can combine (12) with the problem of estimation

$l(\gamma)=\mathbf{P}_{u}(S(X)\leqslant\gamma)=\sum_{x}I_{\\{S(x)\leqslant\gamma\\}}f(x,u)=\mathbf{E}_{u}I_{\\{S(x)\leqslant\gamma\\}},$

where $\mathbf{P}_{u}$ is a measure of the probability in which the random
state $X$ has a density $\\{f((\cdotp;u)\\}$ and $\gamma$ is a known or
unknown parameter. $l$ we estimate using the Kullback-Leibler distance. The
Kullback-Leibler distance is defined as follows:

$D(g,h)=\int g(x)\ln g(x)dx-\int g(x)h(x)dx$

Usually $g(x)$ is chosen from the family of density $f(\cdotp;v)$. So here $D$
for $g$ and $f(x;v)$ it comes down to selecting such a reference parameter $v$
for which $-\int g(x)\ln f(x;v)dx$ is the smallest, that is, it comes down to
maximization:

$\max_{v}\int g(x)f(x;v)dx.$

After some transformations:

$\max_{v}D(v)=\max_{v}\frac{1}{N}\sum_{i=1}^{N}\mathbb{I}_{\\{S(X_{i})\geqslant\gamma\\}}W(X_{i};u;w)\ln
f(X_{i};v)$

So firstly, let’s generate a pair $\\{(\gamma_{t},u_{t})\\}$, which we then
update until the stopping criterion is met and the optimal pair is obtained
$\\{(\gamma^{*},u^{*})\\}$. More precisely, arrange $u_{0}$ and choose not too
small $\varrho$ and we proceed as follows:

1) Updating $\gamma_{t}$

Generate a sample $X_{1},\dots,X_{N_{2}}$ from $\\{f((\cdotp;u_{t-1})\\}$.
Calculate $\hat{S}(X_{1}),\dots,\hat{S}(X_{N_{2}})$ and sort in ascending
order. For $\hat{\gamma}_{t}$ choose

$\hat{\gamma}_{t}=\hat{S}_{(\lceil\varrho N\rceil)}$

2) Updating $u_{t}$

$\hat{u}_{t}$ obtain from the Kullback-Leibler distance, that is, from
maximization

$\displaystyle\max_{u}D(u)$
$\displaystyle=\max_{u}\mathbf{E}_{u_{t-1}}I_{\\{\hat{S}(x)\leqslant\gamma_{t}\\}}lnf(X;u),$
(13) so $\displaystyle\max_{u}\hat{D}(u)$
$\displaystyle=\max_{u}\frac{1}{N_{2}}\sum_{n=1}^{N_{2}}I_{\\{\hat{S}(x)\leqslant\hat{\gamma_{t}}\\}}lnf(X_{n};u).$
(14)

As in Safronov et al. (2006), here a 3-dimensional matrix of parameters
$u=\\{u_{ijl}\\}$ is consider.

$\displaystyle u_{ijl}$ $\displaystyle=\Pr\\{X_{j}^{(i)}=l\\},\text{
$i=1,\dots,k;$}\quad\text{ $j=k-i+1,\dots,N-i;l=0,\dots,N-1$.}$ It seems that
$\displaystyle f(x^{(i)}_{j};u)$
$\displaystyle=\sum_{l=0}^{N-1}u_{ijl}I_{\\{x_{j}^{(i)}=l\\}}.$ and then after
some transformations $\displaystyle\hat{u}^{(t)}_{ijl}$
$\displaystyle=\dfrac{\sum_{n=1}^{N_{2}}I_{\\{\hat{S}(X_{n})\leqslant\hat{\gamma}_{t}\\}}W^{(t-1)}_{nij}I_{\\{X_{nij}=l\\}}}{\sum_{n=1}^{N_{2}}I_{\\{\hat{S}(X_{n})\leqslant\hat{\gamma}_{t}\\}}},$
$\displaystyle W_{nij}^{(t-1)}$
$\displaystyle=\frac{\hat{u}_{ijX_{nij}}^{(0)}}{\hat{u}_{ijX_{nij}}^{(t-1)}},$

where $X_{n}=\\{X_{nij}\\}$, $X_{nij}$ is a random variable from
$f(x_{j}^{(i)};\hat{u}_{t-1})$, corresponding to the formula (14). Instead of
updating a parameter use the following smoothed version

$\hat{u}_{t}=\alpha\hat{u}_{t}+(1-\alpha)\hat{u}_{t-1}.$

3) Stopping Criterion

The criterion is from Rubinstein (1999), which stop the algorithm when
$\hat{\gamma_{T}}$ ($T$ is last step) has reached stationarity. To identify
the stopping point of $T$, consider the following moving average process

$B_{t}(K)=\frac{1}{K}\sum_{s=t-K+1}^{t}\hat{\gamma}_{t},\,t=s,s+1,\dots,\,s\geqslant
K,$

where $K$ is fixed.

$C_{t}(K)=\dfrac{\dfrac{1}{K-1}\\{\sum_{s=t-K+1}^{t}(\hat{\gamma}_{t}-B_{t}(K))^{2}\\}}{B_{t}(K)^{2}}.$

Then let’s define

$C_{t}^{-}(K,R)=\min_{j=1,\dots,R}C_{t+j}(K)$

and

$C_{t}^{+}(K,R)=\max_{j=1,\dots,R}C_{t+j}(K),$

where $R$ is fixed.
Then the stopping criterion is defined as follows

$T=\min\\{t:\dfrac{C_{t}^{+}(K,R)-C_{t}^{-}(K,R)}{C_{t}^{+}(K,R)}\leqslant\epsilon\\}$
(15)

where $K$ and $R$ are fixed and $\epsilon$ is not too small.

###  A.3  The vehicle routing problem.

The classical vehicle routing problem (VRP) is defined on a graph
$\mathbb{G}=(\mathbf{V},\mathcal{A})$, where
$\mathbf{V}=\\{v_{0},v_{1},\ldots,v_{n}\\}$ is a set of vertices and
$\mathcal{A}=\\{(v_{i},v_{j}):i,j\in\\{0,\ldots,n\\}$,
$v_{i},v_{j}\in\mathbf{V}\\}$ is the arc set. A matrix $\mathbf{L}=(L_{i,j})$
can be defined on $\mathcal{A}$, where the coefficient $L_{i,j}$ defines the
distance between the nodes $v_{i}$ and $v_{j}$ and is proportional to the cost
of travelling by the corresponding arc. There can be one or more vehicles
starting off from the depot $v_{0}$ with a given capacity, visiting all or a
subset of the vertices, and returning to the depot after having satisfied the
demands at the vertices. The _Stochastic Vehicle Routing Problem_ (SVRP)
arises when elements of the vehicle routing problem are stochastic – the set
of customers visited, the demands at the vertices, or the travel times.

Let us consider that a certain type of a product is distributed from a plant
to $N$ customers, using a _single_ vehicle, having a fixed capacity $Q$. The
vehicle strives to visit all the customers periodically to supply the product
and replenish their inventories. On a given periodical trip through the
network, on visiting a customer, an amount equal to the demand of that
customer is downloaded from the vehicle, which then moves to the next site.
The demands of a given customer during each period are modelled as independent
and identically distributed random variables with known distribution. A
reasonable assumption is that all the customers’ demands belong to a certain
distribution (say normal) with varying parameters for different customers.

The predominant approach for solving the SVRP class of problems is to use a
“here-and-now” optimization technique, where the sequence of customers to be
visited is decided in advance. On the given route, if the vehicle fails to
meet the demands of a customer, there is a _recourse action_ taken. The
recourse action could be in the form of _going back to the depot_ to replenish
and fulfil the customers’ demand, and continue with the remaining sequence of
the customers to be visited or any other meaningful formulation. The problem
then reduces to a stochastic optimization problem where the sequence with the
minimum expected distance of travel (or equivalently, the minimum expected
cost) has to be arrived at. The alternative is to use a re-optimization
strategy where upon failure at a node, the optimum route for the remaining
nodes is recalculated. The degree of re-optimization varies. At one extreme is
the use of a dynamic approach where one can re-optimize at any point, using
the newly obtained data about customer demands, or to re-optimize after
failure. Neuro Dynamic Programming has been used to implement techniques based
on re-optimization (cf., e.g., Secomandi (2000)).

The vehicle departs at constant capacity Q and has no knowledge of the
requirements it will encounter on the route, except for the probability
distributions of individual orders. Hence, there is a positive probability
that the vehicle runs out of the product along the route, in which case the
remaining demands of that customer and the remaining customers further along
the route are not satisfied (a failure route). Such failures are discouraged
with penalties, which are functions of the recourse actions taken. Each
customer in the set can have a unique penalty cost for not satisfying the
demand. The cost function for a particular route travelled by the vehicle
during a period is calculated as the sum of all the arcs visited and the
penalties (if any) imposed. If the vehicle satisfies all the demands on that
route, the cost of that route will simply be the sum of the arcs visited
including the arc from the plant to the first customer visited, and the arc
from the last customer visited back to the plant.

Alternatively, if the vehicle fails to meet the demands of a particular
customer, the vehicle heads back to the plant at that point, terminating the
remaining route. The cost function is then the sum of all the arcs visited
(including the arc from the customer where the failure occurred back to the
plant) and the penalty for that customer. In addition, the penalties for the
remaining customers who were not visited will also be imposed. Thus, a given
route can have a range of cost function values associated with it. The
objective is to find the route for which the expected value of the cost
function is minimum compared to all other routes.

Let us adopt the “here-and-now” optimization approach. It means that the
operator decides the sequence of vertices to be visited in advance and
independent of the demands encountered. Such approach leads to a _discrete
stochastic optimization problem_ with respect to the discrete set of finite
routes that can be taken. Let $G(r):=\mathbf{E}[H(r,D_{r})]$, where
$H(r,D_{r})$ is the deterministic cost function of the route
$r=(r_{1},\ldots,r_{N})$, $R$ – is the discrete and finite (or countable
finite) feasible set of values that r can take, and demands $D_{r_{i}}$ are a
random variable that may or may not depend on the parameters $r_{i}$. The
class of discrete stochastic optimization problems we are considering here are
those of the form

$\displaystyle\min_{r\in R}\\{G(r)\\}.$

Thus, $H(r,D_{r})$ is a random variable whose expected value, $G(r)$ is
usually estimated by Monte Carlo simulation. The global optimum solution set
can then be denoted by

$\displaystyle R^{*}=\\{r*\in R:G(r^{*})\leq G(r),\forall_{r\in R}\\}.$

The determination of $R^{*}$ can be done by CE (v. Rubinstein Rubinstein
(1997, 1999, 2002)). In this case also the basic idea is to connect the
underlying optimization problem to a problem of the estimating rare-event
probabilities (v. de Boer et al. (2005) in the context of the deterministic
optimization). Here the CE method is used in the context of the _discrete
stochastic optimization_ and the Monte Carlo techniques are needed to estimate
the objective function. In solving the problem by CE methods some practical
issues should be solved, as to when to draw new samples and how many samples
to use.

References

## References

  * Martino et al. (2018) Martino, L.; Luengo, D.; Míguez, J., Introduction.  In Independent Random Sampling Methods; Springer International Publishing: Cham, 2018; pp. 1–26.  doi:black10.1007/978-3-319-72634-2_1.
  * Metropolis (1987) Metropolis, N.  The beginning of the Monte Carlo method.  Los Alamos Sci. 1987, 15, 125–130.  Special Issue: Stanislaw Ulam 1909–1984.
  * Metropolis and Ulam (1949) Metropolis, N.; Ulam, S.  The Monte Carlo method.  J. Am. Stat. Assoc. 1949, 44, 335–341.  Zbl 0033.28807, doi:black10.2307/2280232.
  * Marshall (1956) Marshall, A.  The use of multistage sampling schemes in Monte Carlo computations.  Symposium on Monte Carlo; Wiley: New York, 1956; p. 123–140.
  * Rubinstein (1997) Rubinstein, R.Y.  Optimization of computer simulation models with rare events.  European Journal of Operational Research 1997, 99, 89–112.  doi:black10.1016/S0377-2217(96)00385-2.
  * Sani and Kroese (2007) Sani, A.; Kroese, D.  Optimal Epidemic Intervention of HIV Spread Using Cross-Entropy Method.  Proceedings of the International Congress on Modelling and Simulation (MODSIM); Oxley, L.; Kulasiri, D., Eds.; Modelling and Simulation Society of Australia and New Zeeland: Christchurch, New Zeeland, 2007; pp. 448–454.
  * Asamoah et al. (2018) Asamoah, J.; Nyabadza, F.; Seidu, B.; Chand, M.; Dutta, H.  Mathematical Modelling of Bacterial Meningitis Transmission Dynamics with Control Measures.  Computational and Mathematical Methods in Medicine 2018, Article ID 2657461, 21 page.  PMID: 29780431 , PMCID: PMC5892307 , doi:black10.1155/2018/2657461.
  * Vereen (2008) Vereen, K.  An SCIR Model of Meningococcal Meningitis.  Master’s thesis, Virginia Commonwealth University, Richmond, Virginia 23284, 2008.  Electronic theses and dissertations (ETDs) of VCU.
  * Casagrandi et al. (2006) Casagrandi, R.; Bolzoni, L.; Levin, S.A.; Andreasen, V.  The SIRC model and influenza A.  Math. Biosci. 2006, 200, 152–169.  MR 2225509, PMID: 16504214 , doi:black10.1016/j.mbs.2005.12.029.
  * Parry (1969) Parry, W.  Entropy and generators in ergodic theory; Mathematics Lecture Note Series, W. A. Benjamin, Inc.: New York-Amsterdam, 1969; pp. xii, 124 p.  Zbl 0175.34001.
  * Cover and Thomas (2006) Cover, T.M.; Thomas, J.A.  Elements of information theory, second ed.; Wiley-Interscience [John Wiley & Sons], Hoboken, NJ, 2006; pp. xxiv+748.  doi:black10.1002/047174882X.
  * Kullback and Leibler (1951) Kullback, S.; Leibler, R.A.  On information and sufficiency.  Ann. Math. Stat. 1951, 22, 79–86.  Zbl 0042.38403, doi:black10.1214/aoms/1177729694.
  * Inglot (2014) Inglot, T.  Teoria informacji a statystyka matematyczna.  Mathematica Applicanda 2014, 42, 115–174.  doi:black10.14708/ma.v42i1.521.
  * Csiszár (1967) Csiszár, I.  Information-type measures of difference of probability distributions and indirect observations.  Stud. Sci. Math. Hung. 1967, 2, 299–318.  Zbl 0157.25802.
  * Amari (1985) Amari, S.i.  Differential-geometrical methods in statistics; Vol. 28, Lecture Notes in Statistics, Springer-Verlag, New York, 1985; pp. v+290.  MR 788689, doi:black10.1007/978-1-4612-5056-2.
  * Rubinstein (1999) Rubinstein, R.  The cross-entropy method for combinatorial and continuous optimization.  Methodol. Comput. Appl. Probab. 1999, 1, 127–190.  Zbl 0941.65061, doi:black10.1023/A:1010091220143.
  * Ferguson, T.S. (1989) Ferguson, T.S..  Who solved the secretary problem?  Statistical Science 1989, 4, 282–296.
  * Szajowski, K. (1982) Szajowski, K..  Optimal choice problem of a-th object.  Matem. Stos. 1982, 19, 51–65.  in Polish, doi:black10.14708/ma.v10i19.1533.
  * Polushina (2010) Polushina, T.V.  Estimating Optimal Stopping Rules in the Multiple Best Choice Problem with Minimal Summarized Rank via the Cross-Entropy Method. In Exploitation of Linkage Learning in Evolutionary Algorithms; Chen, Ying-Ping., Ed.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2010; pp. 227–241.  doi:black10.1007/978-3-642-12834-9_11.
  * Stachowiak (2019) Stachowiak, M.  The cross-entropy method and its applications.  Technical report, Faculty of Pure and Applied Mathematics, Wrocław University of Science and Technology, Wroclaw, 2019.  33p. Engineering Thesis.
  * Dror (2002) Dror, M.  Vehicle Routing with Stochastic Demands: Models & Computational Methods. In Modeling Uncertainty; Dror, M. and L’Ecuyer, P. and Szidarovszky, F.., Ed.; Springer: New York, NY, 2002; Vol. 46, International Series in Operations Research & Management Science, pp. 625–649.  doi:black10.1007/0-306-48102-2_25.
  * Chepuri and Homem-de Mello (2005) Chepuri, K.; Homem-de Mello, T.  Solving the vehicle routing problem with stochastic demands using the cross-entropy method.  Ann. Oper. Res. 2005, 134, 153–181.  doi:black10.1007/s10479-005-5729-7.
  * Ekeland and Temam (1999) Ekeland, I.; Temam, R.  Convex analysis and variational problems; Vol. 28, SIAM, 1999.
  * Glowinski (2008) Glowinski, R.  Lectures on numerical methods for non-linear variational problems; Springer Science and Business Media, 2008.
  * Mumford and Shah (1989) Mumford, D.; Shah, J.  Optimal approximations by piecewise smooth functions and associated variational problems.  Communications on Pure and Applied Mathematics 1989, 42, 577–685.
  * Klimov et al. (1999) Klimov, A.; Simonsen, L.; Fukuda, K.; Cox, N.  Surveillance and impact of influenza in the United States.  Vaccine 1999, 17, S42–S46.  Suplement 1. PMID: 10471179 , doi:black10.1016/S0264-410X(99)00104-8.
  * Simonsen et al. (1997) Simonsen, L.; Clarke, M.J.; Williamson, G.D.; Stroup, D.F.; Arden, N.H.; Schonberger, L.B.  The impact of influenza epidemics on mortality: introducing a severity index.  American Journal of Public Health 1997, 87, 1944–1950.  PMID: 9431281 , doi:black10.2105/AJPH.87.12.1944.
  * Earn et al. (2002) Earn, D.J.; Dushoff, J.; Levin, S.A.  Ecology and evolution of the flu.  Trends in Ecology & Evolution 2002, 17, 334–340.
  * Andreasen et al. (1997) Andreasen, V.; Lin, J.; Levin, S.A.  The dynamics of cocirculating influenza strains conferring partial cross-immunity.  J. Math. Biol. 1997, 35, 825–842.  MR 1479341, doi:black10.1007/s002850050079.
  * Lin et al. (1999) Lin, J.; Andreasen, V.; Levin, S.A.  Dynamics of influenza A drift: the linear three-strain model.  Mathematical Biosciences 1999, 162, 33–51.
  * Iacoviello and Stasio (2013) Iacoviello, D.; Stasio, N.  Optimal control for SIRC epidemic outbreak.  Computer Methods and Programs in Biomedicine 2013, 110, 333–342.  PMID: 23399104 , doi:black10.1016/j.cmpb.2013.01.006.
  * Zaman et al. (2008) Zaman, G.; Kang, Y.H.; Jung, I.H.  Stability analysis and optimal vaccination of an SIR epidemic model.  BioSystems 2008, 93, 240–249.
  * Kamien and Schwartz (1991) Kamien, M.I.; Schwartz, N.L.  Dynamic Optimization: The Calculus of Variations and Optimal Control in Economics and Management, second ed.; Vol. 31, Advanced Textbooks in Economics, North-Holland Publishing Co., Amsterdam, 1991; pp. xviii+377.  MR 1159711.
  * Fang et al. (2002) Fang, Q.; Tsuchiya, T.; Yamamoto, T.  Finite difference, finite element and finite volume methods applied to two-point boundary value problems.  J. Comput. Appl. Math. 2002, 139, 9–19.  MR 1876870, doi:black10.1016/S0377-0427(01)00392-2.
  * Caglar et al. (2006) Caglar, H.; Caglar, N.; Elfaituri, K.  B-spline interpolation compared with finite difference, finite element and finite volume methods which applied to two-point boundary value problems.  Applied Mathematics and Computation 2006, 175, 72–79.  doi:black10.1016/j.amc.2005.07.019.
  * Papadopoulos and Giovanis (2018) Papadopoulos, V.; Giovanis, D.G.  Stochastic finite element methods; Mathematical Engineering, Springer, Cham, 2018; pp. xxi+138.  An introduction, doi:black10.1007/978-3-319-64528-5.
  * Lenhart and Workman (2007) Lenhart, S.; Workman, J.T.  Optimal control applied to biological models; CRC Press, 2007.
  * Hazelbag et al. (2020) Hazelbag, C.M.; Dushoff, J.; Dominic, E.M.; Mthombothi, Z.E.; Delva, W.  Calibration of individual-based models to epidemiological data: A systematic review.  PLoS Computational Biology 2020, 16, e1007893.
  * Taylor et al. (2010) Taylor, D.C.; Pawar, V.; Kruzikas, D.; Gilmore, K.E.; Pandya, A.; Iskandar, R.; Weinstein, M.C.  Methods of model calibration.  Pharmacoeconomics 2010, 28, 995–1000.
  * Taynitskiy et al. (2017) Taynitskiy, V.; Gubar, E.; Zhu, Q.  Optimal Impulsive Control of Epidemic Spreading of Heterogeneous Malware.  IFAC-PapersOnLine 2017, 50, 15038 – 15043.  20th IFAC World Congress, doi:black10.1016/j.ifacol.2017.08.2515.
  * Gubar et al. (2018) Gubar, E.; Taynitskiy, V.; Zhu, Q.  Optimal Control of Heterogeneous Mutating Viruses.  Games 2018, 9.  doi:black10.3390/g9040103.
  * Kochańczyk et al. (2020) Kochańczyk, M.; Grabowski, F.; Lipniacki, T.  Dynamics of COVID-19 pandemic at constant and time-dependent contact rates.  Math. Model. Nat. Phenom. 2020, 15, 28.  doi:black10.1051/mmnp/2020011.
  * Kermack and McKendrick (1991) Kermack, W.; McKendrick, A.  Contributions to the mathematical theory of epidemics–I. 1927.  Bull Math Biol. 1991, 53, 35–55.  Reprint of the Proc. of the Royal Society of London. Series A, 1927, Vol. 115, No. 772, pp. 700-721. 10.1098/rspa.1927.0118. The same journal published the parts II and III of the article in 1932 and 1933, respectively. PMID: 2059741 ; JFM 53.0517.01, doi:black10.1007/BF02464423.
  * Harko et al. (2014) Harko, T.; Lobo, F.S.N.; Mak, M.K.  Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates.  Applied Mathematics and Computation 2014, 236, 184 – 194.  MR 3197716, doi:black10.1016/j.amc.2014.03.030.
  * Nikolaev, M. L. (1998) Nikolaev, M. L. .  On optimal multiple stopping of Markov sequences.  Theory Probab. Appl. 1998, 43, 298–306.  Translation from Teor. Veroyatn. Primen. 43, No. 2, 374-382 (1998). Zbl 0971.60046, doi:black10.1137/S0040585X9797691X.
  * Nikolaev (1998) Nikolaev, M.L.  Optimal multi-stopping rules.  Obozr. Prikl. Prom. Mat. 1998, 5, 309–348.  Zbl 1075.91511.
  * Safronov et al. (2006) Safronov, G.Y.; Kroese, D.P.; Keith, J.M.; Nikolaev, M.L.  Simulations of thresholds in multiple best choice problem.  Obozr. Prikl. Prom. Mat. 2006, 13, 975–982.  Zbl 1161.60333.
  * Secomandi (2000) Secomandi, N.  Comparing neuro-dynamic programming algorithms for the vehicle routing problem with stochastic demands.  Computers & Operations Research 2000, 27, 1201 – 1225\.  doi:black10.1016/S0305-0548(99)00146-X.
  * Rubinstein (2002) Rubinstein, R.Y.  Cross-entropy and rare events for maximal cut and partition problems.  ACM Trans. Model. Comput. Simul. 2002, 12, 27–53.  Zbl 1390.90482, doi:black10.1145/511442.511444.
  * de Boer et al. (2005) de Boer, P.T.; Kroese, D.P.; Mannor, S.; Rubinstein, R.Y.  A tutorial on the cross-entropy method.  Ann. Oper. Res. 2005, 134, 19–67.  MR 2136658, doi:black10.1007/s10479-005-5724-z.

The implementation codes for the algorithms used in the examples are available
from the authors.