# Data-driven model for hydraulic fracturing design optimization. Part II:
Inverse problem

V.M. Duplyakov  A.D. Morozov  D.O. Popkov  E.V. Shel  A.L. Vainshtein  E.V.
Burnaev  A.A. Osiptsov<EMAIL_ADDRESS>G.V. Paderin  Skolkovo
Institute of Science and Technology (Skoltech), 3 Nobel Street, 143026,
Moscow, Russian Federation  Gazpromneft Science & Technology Center, 75-79
liter D Moika River emb., St Petersburg, 190000, Russian Federation

###### Abstract

We describe a stacked model for predicting the cumulative fluid production for
an oil well with a multistage-fracture completion based on a combination of
Ridge Regression and CatBoost algorithms. The model is developed based on an
extended digital field data base of reservoir, well and fracturing design
parameters. The database now includes more than 5000 wells from 23 oilfields
of Western Siberia (Russia), with 6687 fracturing operations in total.
Starting with 387 parameters characterizing each well, including construction,
reservoir properties, fracturing design features and production, we end up
with 38 key parameters used as input features for each well in the model
training process. The model demonstrates physically explainable dependencies
plots of the target on the design parameters (number of stages, proppant mass,
average and final proppant concentrations and fluid rate). We developed a set
of methods including those based on the use of Euclidean distance and
clustering techniques to perform similar (offset) wells search, which is
useful for a field engineer to analyze earlier fracturing treatments on
similar wells. These approaches are also adapted for obtaining the
optimization parameters boundaries for the particular pilot well, as part of
the field testing campaign of the methodology. An inverse problem (selecting
an optimum set of fracturing design parameters to maximize production) is
formulated as optimizing a high dimensional black box approximation function
constrained by boundaries and solved with four different optimization methods:
surrogate-based optimization, sequential least squares programming, particle
swarm optimization and differential evolution. A recommendation system
containing all the above methods is designed to advise a production
stimulation engineer on an optimized fracturing design.

###### keywords:

hydraulic fracturing , machine learning , predictive modelling , design
optimization , gradient-free optimization , probability of improvement ,
surrogate optimization , multistage fracturing

††journal: Journal of Petroleum Science & Engineering

##  1 Introduction and problem formulation

In this paper, we continue the work on the development of a workflow on
hydraulic fracturing (HF) design optimization with machine learning on field
data, which was initiated in [1]. In Part I of this project, we presented the
detailed methodology for development of a digital database on reservoir and
well parameters, design parameters of multistage fracturing treatments and
production data. We also discussed in detail the forward problem of production
forecast based on reservoir, well and fracturing design data. Here we continue
the effort and now move on to the development of the optimization workflow
(inverse problem).

Before we move to our results, let us review the recent developments in the
area of machine learning-assisted optimization of hydraulic fracturing
designs. In [1], we made a thorough review of the relevant literature
published on this subject, so the reader is referred to [1] for the state of
the art, whereas here we will only mention new papers appeared after Part I of
our study has been published.

The paper [2] uses candidate selection for refracturing operations based on
data analytics. Gradient boosting is used for production forecast, and the
feature importance analysis indicates the current production (before
refracturing) to be the most important parameter. We would like to emphasize
that the problem of candidate selection for refracturing treatments and the
problem of primary fracturing design optimization are the two very different
problems by formulation: the prior production is known when it comes to
refracturing, which makes the task of predicting subsequent production much
easier.

In [3], the multiobjective random forest method is proposed to predict the
dynamic production data for a shale gas well, with geological and hydraulic
fracturing properties used as input features. The performance of multi-
objective random forest (MORF) and multi-output regression chain (MORC)
methods are compared. A dynamic regression model is constructed in [4] using a
machine learning method referred to as the sliding window regression. Features
of the model are divided into four groups: integral, local, pressure and
autoregressive. The idea of sliding window is used to obtain a model with
stable coefficient dynamics and ability for long-term forecast.

In [5], variable-length particle-swarm optimization (Modified Variable-length
PSO, MVPSO) was proposed to automatically select the optimal fracturing
parameters: the number of fractures as well as the corresponding fracture
properties. Then, MVPSO was verified and compared with VPSO by several
benchmarks. In addition, a gas/water two-phase model considering gas-
adsorption and Knudsen-diffusion effects was used to describe the shale-gas
flow in matrix and fracture domains.

About this work: the paper is organized as follows. The status of the
development of a digital database is summarized in Sec. 2. We describe the
model for production forecast (forward problem) in Sec. 3. A study on offset
wells selection for finding optimization intervals for fracturing design
parameters is presented in Sec. 4. Inverse problem of finding the optimized
fracturing design parameters is discussed in Sec. 5. Overall discussion and
analysis of the results are given in Sec. 6. The paper ends up with summary
and conclusions in Sec. 7.

##  2 Data

###  2.1 Database overview

In this work, we use the version 7.9 of the database. Compared to the database
version 5.5 used in our previous publication [1], the overall architecture and
the sources of data used to compile the database for production forecast have
not been changed. The database still contains information on more than six
thousands miscellaneous fracturing operations (single-stage and multi-stage,
primary and re-fracturing) or about 17 000 injections on 5 425 wells from 23
different fields in Western Siberia for the period 2013-2019, including well
profile, fracturing design, geological and production reports. Also, the
centralized approach to data organization is remained unchanged.

Our work on improvements of the database have been continued and the quality
of the preprocessed data was increased so far through consultations with
engineers ingrained in the raw data collection process.

Current version of the database 7.9 contains the following changes:

  * 1.

Target preprocessing is updated: monthly data consolidated in 3-, 6- and
12-months slices are summed also for all active reservoirs on the well and not
only the stimulated ones. The shape of the target distribution, predictably,
has become smoother by moving the peak from almost 1,000 to near 2,000 cubic
metres of fluid, as shown in Fig. 1;

  * 2.

There was an issue with the lack of standardization in reservoir’s and well’s
namings among the data sources. Similarity between labels is now achieved
through morphological analysis;

  * 3.

Added new features: proppant properties (grain size and density of prevailing
type by proppant manufacturer);

  * 4.

Added facies for the dominant field. The feature is described by five zones
from lithological maps: slope zone, proximal fan, shallow area, alluvial fan
and depocenters.

###  2.2 Target variable

Selection of the right target variable is crucial for the success of the
entire optimization workflow, so we think it deserves a separate discussion in
this subsection. We analyzed the data made available to us and also discussed
in detail the current strategies of fracturing design development among
production stimulation engineers and come to a conclusion that fracturing
operations are optimized in the field based on the metrics of maximizing
reservoir contact: the larger the fracture the better. Hence, bigger fractures
yield larger cumulative volume of recovered total fluid (both oil and
reservoir water), and there is a strong correlation. Rarely this consideration
takes into account the fact the excessively large fractures may breakthrough
into upper or lower layers, thereby causing additional production of water.
Thus, a model targeted to predict the cumulative volume of total produced
fluid was trained on existing data with slightly higher accuracy (as there is
a direct correlation between the input parameters characterizing the
fracturing design and the total cumulative produced fluid), as compared to a
model targeted on cumulative oil only, where the prediction accuracy was
lower, as the oil production is not directly correlated with the design
parameters governing the fracture dimensions. The present realization of the
model does not take into account (yet) the presence of water-bearing layers,
so the model better predicts the total volume of produced fluid. During the
pilot works, we deal with the particular oilfield where there are no bottom
waters nearby target formation, so the targets of cumulative total fluid and
cumulative oil are practically equivalent.

At the same time, we realize that generalization of the present workflow to
oilfields with the presence of aqueous layers will require modifications: we
will need to add some features characterizing the presence of upper/lower
water bearing formations to be able to predict the production of oil and water
separately, and also the target variable should be composed of two components:
maximum total fluid and maximum pure oil (or minimum water cut, which is
equivalent).

We use the 3-month cumulative production data to utilize all the data
available, including the most recently fractured wells, where the production
history is short.

An interpolation of the monthly data was used in order to come up with 90-day
production period. This approach is suitable, for instance, if a well worked
for 20 days in the first month after fracturing treatment and 30 days in the
following months, the desired value is in 80-110 days or 3-4 months range and
can be found through the interpolation.

Figure 1: Target probability distributions for database v. 5.5 used in [1]
(red) & the current database v. 7.9 (blue).

###  2.3 Feature selection

In the first part of this study [1], we made efforts to utilize the entire
database with primary fracturing operations, as well as refracturing to solve
the forward problem of predicting the production rate without any preliminary
separation of the dataset.

Currently, the database is divided into two parts: primary and repeated
stimulations. We fitted models for prediction of the fluid production for both
these types of operations. Further, during the field tests, we focused only on
primary stimulation treatments optimization.

To construct the models and increase interpretability of the problem we use
feature importance analysis. It can be done with or without the involvement of
an approximation model.

One method is a sensitivity analysis via Sobol indices [6], which can be
calculated without constructing any approximation model and decompose the
variance of the target variable into parts attributed to input features. The
most important features can be seen in Fig. 2.

On the contrary, SHAP method [7] is based on an approximation model, utilizes
the concept of Shapley values [8] and measures features importances in terms
of predictive power of each feature. The use of this method more accurately
shows the true importance of the features [9]. The SHAP values can be
calculated for tree-based models (which we use).

Indeed, the difference between the two models for primary stimulation vs
refracturing can be seen through the SHAP feature importance analysis (Figs.
3, 4). (In case of categorical features we use one-hot encoding. The
corresponding features are prefixed with “cat.”.) Particularly, for
refracturing operations the key feature having major impact on the target is
the level of production before the refracturing treatment (which was not
captured in previous analysis), which is absent in case of primary operations.
Having data on production prior to refracturing makes the production forecast
problem easier to solve, compared to the case of production forecast after
primary fracturing operations. This has also been noted in other studies [10].

Based on the analysis, we may conclude that the features describing proppant
properties (introduced in Sec. 2.1) are not important for predicting the fluid
production.

Another way of analyzing parameters is the feature elimination procedure which
implies a reduction of feature space which in its turn may improve the
performance of the approximation model. Feature elimination can be achieved by
applying a pair-wise correlation of parameters via the Spearman correlation
and by the Recursive Feature Elimination (RFE) method which involves the use
of the approximation model. Speaking of the former, the reason of choosing the
Spearman correlation instead of the Pearson correlation is because most of the
features correlations are non-linear. So, by estimating correlations between
features we can remove perfectly correlated features and features with zero
variance. Regarding the second method, RFE is a procedure for backward
selection of features which works as follows. First, a model is built on an
the entire set of input features, and features’ importance is calculated.
Then, the least important features are removed. After that, the model is
rebuilt on the reduced set of features and we repeat the process.

Eventually, we use the following procedure for feature elimination:

  1. 1.

Select 3-month slices; (6- and 12-months production & geological and technical
data are removed due to the lack of data from the latest treatments);

  2. 2.

Remove parameters which are not relevant for the considered target variable,
or for which more than 80% of the observations are missing;

  3. 3.

Estimate Spearman correlations to identify perfectly correlated features;
remove features with almost zero variance;

  4. 4.

Apply RFE to select features, which are the most important for the target
prediction.

As the result, we reduced our initial set of features from 387 to 38 features.

Figure 2: Sobol sensitivity for the entire database. Figure 3: SHAP feature
importance: refracturing operations. Figure 4: Feature importance: primary
fracturing on new wells.

##  3 Forward problem of oil production forecast after fracturing.

###  3.1 Problem formulation

In this section, we describe the methods developed for prediction of the
cumulative fluid production. Previously [1], we created prediction model based
on gradient boosting algorithm which was based on various features including
reservoir parameters, well construction and HF design data (which includes
many features like pumping pressure, fluid efficiency, fracture geometry,
etc). Those features were useful for predicting the efficiency of HF. However,
in this paper, we focus on choosing the proper HF design only for the new
wells. We need to know the optimal design before the treatment, therefore we
cannot use all the design features, as many of them are known only after the
fracturing. Chosen features for the optimization are:

  * 1.

Number of stages,

  * 2.

Pad share,

  * 3.

Fracturing fluid volume (or average proppant concentration),

  * 4.

Proppant mass,

  * 5.

Fluid rate,

  * 6.

Final proppant concentration (proppant concentration at the last stage of
pumping).

These parameters are essential for planning the HF treatment. Through those
features one can estimate a pumping schedule for a particular well.

One may notice the absence of one of the most important parameters, which
determines the fracturing fluid viscosity — polymer concentration (in water-
based polymer gels). This is due to the very low variability of this parameter
in the initial database.

###  3.2 Stacking

The most successful algorithms in terms of error metrics were decision-tree
based ones for the problems with heterogeneous set of features. For our
prediction task, the CatBoost outperformed all others.

Besides, the model predictions should be physically explainable and robust. To
check these properties we plotted dependencies of the target on the input
design parameters. We expect the plots to demonstrate smooth and robust
behaviour. One can see in Figure 17 (green line) that CatBoost-based plots are
not smooth, as there is not so much initial data and decision-tree based
algorithms create piece-wise constant approximations.

To correct this behaviour of the model we used a stacking approach. In the
first step, we trained a ridge regression model, then we subtracted
predictions of the model from the true target values and then fitted a
CatBoost regression model on the residuals. Hyperparameters tuning was
performed using 5-folds cross validation. Overall model uses both ridge and
CatBoost regression algorithms, summing the predictions. One can see the
improved dependencies plots in Figure 17 (red line). We believe that these
plots are more physically explainable. To characterize quality of the model we
used scores presented in Table 1. These metrics were calculated on the hold-
out set that is 30% of the initial data. See also the scatter plot in Fig. 5.

Quality score | Value
---|---
R2 | 0.64
MAE | 1131
MedianAE | 723.7
MAPE | 36.08
Weighted MAPE [11] | 29.07
MSE | 2802631
RMSE | 1674.11
Table 1: Quality scores for the stacked model Figure 5: Scatter plot for the
stacked model

##  4 Selection of optimization intervals

After constructing the model, the inverse problem is formulated as finding a
set of optimal fracturing design parameters to maximize our target. Since the
model is multi-dimensional, a valid selection of the design parameter values
is only possible if they change within the relevant intervals during the
optimization procedure. These intervals are caused by various constraints,
arising in the field.

The additional restrictions on design parameter values can be divided into
geological and technological constraints. Geological constraints include those
related to the geological structure of the formation. For example, proximity
of gas or water bearing formations leads to the necessity to limit fracture
height growth, which automatically leads to limitation of the maximum volume
of injected proppant and requires change of the perforation strategy.
Geological constraints can also be related to waterflooding cases. For
example, when an injection well, operating at bottomhole pressures higher than
the formation breakdown pressure, is located relatively close to the
production well and is in the direction of fracture propagation, it is
necessary to limit the maximum fracture half-length on the production well.

Technological limitations include those related to the technical capabilities
of the equipment and chemicals used. For instance, the maximum fracturing
pressure (first of all, it is connected with the capabilities of the wellhead
equipment and pumping units) can lead to the restriction on the maximum
fracture width. Those limitations at a first approximation could be obtained
by offset (similar) well search methods.

In the field, a fracturing job is pumped with variable proppant concentration,
so another limitation we add is a parameter characterising the rate of
increase in the proppant concentration in the fracturing fluid from initial
($c_{\\_}{start}$) to final ($c_{\\_}{fin}$) concentration $[kg/m^{3}]$ with
respect to increase to average $c_{\\_}{avg}$. The parameter is denoted as
$\epsilon$ and is defined as:

$\epsilon=\frac{c_{\\_}{fin}-c_{\\_}{start}}{c_{\\_}{avg}-c_{\\_}{start}}-1.$
(1)

Boundaries for this parameter are the same for all pilot wells: from 0.5 to
1.5. This is not a parameter that needs to be optimized, we add these bounds
as a constraint for optimization algorithms.

###  4.1 Clustering for offset wells selection

To find optimization intervals for fracturing design parameters we define a
pilot cluster as a set of wells which are similar to the pilot one (the same
field, layer, face and direction). We estimate design parameters limits as
$5^{th}$ and $95^{th}$ percentiles of the values of the parameters for wells,
belonging to the constructed cluster. Thus, we believe, the results of the
optimization would be more robust, as we do not use extrapolated model’s
predictions.

Another approach to obtain optimization boundaries for the specified pilot
well is to utilize unsupervised ML methods like clustering and dimensionality
reduction. The procedure is as follows (See also figs. 6, 7):

Figure 6: Cluster procedure: steps 1 – 5, 9 Figure 7: Cluster procedure: steps
6 – 8

  1. 1.

Remove all features from the database except the environment parameters: PVT,
geomechanics, well log interpretation, well id;

  2. 2.

Filter the obtained subset by the number of stages, layer id and face of the
pilot well. This step is needed in order to collect wells technically similar
to the pilot well;

  3. 3.

Apply a clustering algorithm to the environment parameters of the filtered
data subset. Hyperparameters of the clustering method are optimized via a
gradient-free optimization algorithm by maximizing the mean silhouette
coefficient, calculated for a particular cluster as follows:

$S_{\\_}i=\frac{b_{\\_}{i}-a_{\\_}{i}}{max(a_{\\_}i,b_{\\_}i)},$ (2)

where $a_{\\_}i$ is a mean intra-cluster distance, $b_{\\_}i$ is a mean
nearest-cluster distance;

  4. 4.

Find a cluster to which the pilot well belongs to. Assert this cluster as a
pilot well cluster;

  5. 5.

Using id-s of the objects from the pilot cluster, add the design parameters
values to the objects’ descriptions;

  6. 6.

Analyze the pilot cluster statistics: minimum, maximum and mean values of the
design parameters. The minimum values and the maximum values serve as the
optimization boundaries while the mean values of the design parameters are
used as an initialization for the optimization problem.

  7. 7.

Perform a visual analysis via the t-SNE algorithm based on the environment
parameters. t-SNE can only be applied to data without missing values. Hence,
the missing values are imputed by the matrix factorization algorithm [1].
After that we obtain the t-SNE embedding for the selected subset of data. Here
t-SNE embedding is a mapping of the multidimensional input features to the
two-dimensional $x$ and $y$ coordinates for each observation of the data
subset. These two-dimensional coordinates can be used for visualization of the
considered data subset.

  8. 8.

After that, using the considered subset of the data and the obtained t-SNE
coordinates build an ML regression model to predict t-SNE $x$ and $y$
coordinates for any new object. The corresponding ML model takes as input
environment parameters and returns as output t-SNE coordinates $x$ and $y$.
The reason of predicting the t-SNE coordinates is because any change in the
data like appending a new pilot well to the data set will require running the
t-SNE algorithm from scratch;

  9. 9.

Predict the t-SNE coordinates $x$ and $y$ for the pilot well using the
constructed regression model. This allows to visualize the position of a new
well with respect to the selected data subset;

  10. 10.

Create a t-SNE scatter plot from the $x$ and $y$ t-SNE coordinates. Label the
corresponding clusters of each object by some colors and label the pilot well
by a star symbol. This procedure visually verifies how the data is clustered
and how the pilot well is located with respect to the remaining data set.

The example of this method implementation is shown in Fig. 8.

Figure 8: t-SNE scatter plot for cluster visualization

###  4.2 Missing pilot parameters imputation

In practice some of the parameters values of a well can be missing. Imputation
is needed for wells parameters values because gradient-free optimization
algorithms are based on constructing regression models and so they cannot work
with missing values. For better results, we need to impute these values not
just by filling them with corresponding mean values, but we should do it in a
smarter way. We propose several strategies to do it. The matrix factorization
method was described in the first part of our study [1]. The second method is
to impute missing values by averaging parameters values of the top-_N_ similar
wells. The method to find similar wells is described in Sec. 4.3. Another way
to make the imputation is to use the mean pilot cluster parameters values as
described in step 6 of the algorithm in the previous subsection.

###  4.3 Offset wells selection by Euclidean distance

Offset wells are the wells, similar to the pilot one in terms of their
geological surroundings. One may look for these wells both in terms of their
geological and geographical (closest wells within certain radius from the
pilot one) similarities. These analogue wells search is very useful for a
petroleum engineer as it allows to analyse fracturing operations, conducted
previously, and the design parameters values, check whether an operation was
successful or not, etc. We can also extract additional features from the
neighbouring wells, which increase predictive power of the models [10]. For
example, in our work, we used features, such as average fluid production
divided by distance from the pilot well, using wells within 1 km from the
pilot one.

As a similarity metric we can use the Euclidean distance between the pilot
well and the other wells. To calculate it, we firstly need to normalize values
in the database. We perform a linear min-max normalization, where min and max
values are the $1^{st}$ percentile and the $99^{th}$ percentile respectively.
We use percentiles to discard possible outliers in the initial data set. Then,
the Euclidean distance between two vectors of parameters, characterizing
wells, is calculated as

$d(p,q)={\sqrt{(p_{\\_}{1}-q_{\\_}{1})^{2}+\cdots+(p_{\\_}{n}-q_{\\_}{n})^{2}}},$
(3)

where $p_{\\_}i$ and $q_{\\_}i$ are the $i$-th parameter’s values of the
corresponding wells. The distances are calculated between the pilot well and
all other wells, belonging to the same cluster (within the same field, layer
and face). The results of such similar wells search are considered robust and
sustainable by field geologists. One can see an example of a result of such
search in Fig. 18.

In this work we combine the clustering method for selecting offset wells with
the euclidean distance search. The clustering is used for obtaining the set of
similar wells, then we reduce the size of this set, leaving only top-N wells
by the euclidean distance. This method allows us to look for optimal values of
the design parameters in a certain vicinity where the prediction model works
well. Also, in some cases the top-N similar wells can be used for imputing
missing parameters for the pilot well.

##  5 Inverse problem of finding frac design parameters to maximize production

###  5.1 Problem formulation

An inverse problem can be formulated as optimizing a high dimensional black
box (BB) with respect to inputs constrained by boundaries. In our case BB is a
function with unknown expression or internal structure that, given a list of
inputs, returns corresponding outputs. The high dimensionality of the input
presents an exponential difficulty for problem modeling and optimization (so-
called “curse-of-dimensionality”) [12]. To optimize a high-dimensional
computationally expensive black box (HEB) function, it is required to
iteratively evaluate an objective function, which can be costly and so becomes
unacceptable. In our case, the HEB function is represented by the constructed
ML regression. To optimize HEB, we used the following optimization methods:
surrogate-based optimization (Sec. 5.3), sequential least squares programming
[13], particle swarm optimization [14] and differential evolution [15]
algorithms. The advantages of these methods are that they make no assumptions
about the problem being optimized and can perform searches in very large
spaces of candidate solutions.

For our methodology, the goal is to maximize the cumulative fluid production
by finding a set of optimal design parameters constrained by boundaries (see
Sec. 4.1) for the specified parameters of the pilot well environment.

###  5.2 Review of known optimization procedures

There are studies published in open literature on optimization of HF
parameters based on some relatively moderate datasets with hundreds of wells,
e.g., see a study on 570 wells in Clinton sand [16]. The first stage is a
neural network which takes as input frac completions data and production
history and predicts post-frac deliverability. The model serves as a screening
tool for excluding “dog wells”, which cannot be considerably enhanced after a
frac job. The second stage is based on a neural network, where the branch of
the network, responsible for the fracture design parameters, is connected to
the optimization algorithm to obtain the optimized frac design for each well
and the expected post frac deliverability (see Fig. 9).

Figure 9: A schematic diagram of the neuro-genetic approach

The work [17] also utilizes a neuro-genetic algorithm in the optimization
pipeline instead of a computationally complex implicit numerical hydrodynamic
model. In particular, the article presents a case study of neuro-genetic
optimization using the example of a cross-flow microturbine. In another paper
[18] an ML regression analysis on more than 3500 wells was performed, the
authors optimized the design by visual analysis with taking a pair of the most
important parameters and selecting a region of the most optimal design.

Speaking about other industries, the article [19] uses neural networks and a
hybrid multisubgradient descent method with adaptive learning rate to solve
multicriteria optimization of multiple tasks (generation of too large droplets
and too low droplet speed) in the field of bioprinting. The proposed method
can improve both printing accuracy and stability, and is useful for realizing
precise cell arrays and complex biological functions.

Another application of machine learning and optimization is the design
optimization of thin multilayer solar cells to maximize external quantum
efficiency [20]. The authors utilize the concept of transfer learning which
implies the use of experience gained in solving one problem to solve another,
similar problem. The reason why the transfer learning is applied to the
problem is because the problem involves the change of design specifications.
Therefore, the transfer learning model acts as a function of surrogate
optimization, which refits the surrogate more efficiently. In particular, the
procedure improved the results by 2-3 times using only half of the training
samples compared to the usual model.

###  5.3 Surrogate-based optimization (SBO)

Surrogate models (approximation models) [21] can be used with multi-
dimensional input design spaces. As we know, the more parameters the surrogate
model takes as input the more computational resources (training time) it is
required to construct the model.

A particular example of Surrogate-based Optimization (SBO) or sequential
model-based global optimization is Bayesian optimization. The algorithm
utilizes a probabilistic surrogate function which approximates the expensive
objective function and an acquisition function which, in its turn, allows to
select a new candidate input design point within the domain. At the beginning,
a surrogate function is built on a small number of samples from the original
objective function, which provides a target value. Then, the algorithm
maximizes the acquisition function and the surrogate is updated with a new
input point and its actual output value. After repeating the process, an
optimum can be achieved by taking the information about the samples from the
past iterations. A typical choice of the surrogate is Gaussian processes and
random forest, while typical acquisition functions are the Expected
Improvement (EI) and the Probability of Improvement (PI).

We used the SBO algorithm implemented in the industrial software [22]. The SBO
methodology is based on Gaussian processes modeling technique [23, 24, 25,
26]. The particular numerical realization roots in the scientific works
published in [23].

##  6 Results and Discussion

During the research, we trained the model on wells from the field of interest
only, thereby reducing the database from 6687 to 3308 fracturing operations.

Field tests have been carried out on 21 wells. 9 wells were horizontal, 7 —
vertical multilateral, operating on several layers simultaneously, the rest 5
wells were regular vertical ones. We were testing the accuracy of our
prediction models, as well as HF design optimization overall pipeline. This
pipeline includes:

  1. 1.

Obtaining design parameters optimization boundaries by similar wells search
with DBSCAN clustering (Sec. 4.1);

  2. 2.

Imputing missing parameters with mean corresponding parameters values,
calculated using top-10 similar wells, obtained within the pilot cluster via
the euclidean distance search (Sec. 4.3);

  3. 3.

Performing the design optimization, using the predictive model (Sec. 3) and
optimization algorithms (Sec. 6.2)

###  6.1 Forward model prediction accuracy test

The production prediction errors are shown in Fig. 19. Here we see a
relatively low percentage error for horizontal wells, which can be explained
by the higher production rates of such wells. What is important here is the
high error for multilateral wells, which is most likely caused by data
distortion for this type of wells. Currently, the data point for a
multilateral well is represented as a multistage fracture treatment with the
number of stages equal to the number of laterals with fractures. The
disadvantages of this method are obvious when hydraulic fracturing is
performed at different points in time. In addition, different reservoir
parameters must be considered for each operational production formation.

Generally, the accuracy of our model ($MAPE=37.28\%$, $wMAPE=27.46\%$) in
field tests is close to the hold-out set accuracy check ($MAPE=36.08\%$,
$wMAPE=29.06\%$). The distribution of well types (vertical, horizontal) for
field tests is close to that presented in the hold-out set.

###  6.2 Design optimization test

Overall, we formulated four approaches to the problem of design optimization:

  1. 1.

SBO: Surrogate-Based Optimization (Sec. 5.3),

  2. 2.

SLSQ: Sequential Least Squares Programming,

  3. 3.

PSO: Particle Swamp Optimization,

  4. 4.

DE: Differential Evolution.

We compared these methods with each other in terms of maximum produced fluid
and physically-grounded recommendations for the design parameters for 3 well
types: horizontal, vertical and vertical multilateral.

Comparison of the efficiency of optimization algorithms in Figure 10 shows the
average cumulative fluid production across all pilot wells. A larger value
means higher efficiency of an algorithm. Figure 21 shows the extended results
of the optimization for each well individually.

The results of the optimization for design parameters are shown in Figures
11-16. Here we have indicated the recommended (optimum) values in percentages.
0% is a lower and 100% is an upper bound for each well, obtained by the offset
wells search. Color bars represent the average value of the analyzed parameter
for all wells in each of the three groups.

Figure 10:  Average optimized production for all wells by different
optimization algorithms Figure 11: Average recommended fluid rate for the
wells Figure 12: Average recommended mean proppant concentration for the wells
Figure 13: Average recommended proppant masses (per stage) for the wells
Figure 14: Average recommended final proppant concentration for the wells
Figure 15: Average recommended pad share for the wells Figure 16: Average
recommended calculated epsilon

All the methods were operating under similar conditions: boundaries,
constraints, maximum allowable number of function evaluations (200). Under
these conditions, SBO proved to be the most efficient algorithm in terms of
maximizing the objective function (Figures 10 and 21). It is also possible to
draw conclusions about general optimization trends. First of all, the most of
the approaches for all well types maximized the target by increasing the
proppant mass and reducing the average proppant concentration. It makes sense
to increase the amount of pumping fluid and proppant to get maximum production
while we have neither WOC nor GOC.

We can clearly see the recommendation trends difference between horizontal and
both types of vertical wells even though vertical multilateral wells are
represented in our database as multistage fracturing operations (like
horizontal wells). Pad share here is the most interesting parameter, which
does not follow any particular trend. The final proppant concentration is
lower for horizontal wells, which can be explained by the high probability of
ineffective treatments with high concentration values in these types of wells,
as such fluids are difficult to pump through well sections with high wellbore
curvature. The value of the fluid flow rate varies widely and is probably not
important for the purpose of maximizing production. This can be proved by the
low feature importance of this parameter, and may also be due to the fact that
the values of this parameter are chosen properly in most of HF treatments
presented in the database.

Once we have the results of the optimization, we can carry out a kind of
retrospective analysis to find out how we can change usual HF designs to
improve its results: if the value is below 50% – the optimum value is less
than usual one used during HF treatments and vice versa.

We tested the SBO method to find the optimal set of design parameters for each
well from the pilot tests. In addition, we limited the proppant mass to a
value from the actual fracture design (the actual values of the design
parameters have been chosen by the contractor). We then predicted fluid
production for the optimal parameters and for the actual ones. Optimal
parameters showed a higher production rate, theoretically increased by 38%.
(Fig. 20).

##  7 Summary and Conclusions

Since the first part of this work was published in [1], both the digital
database and the production forecast model have undergone several changes: the
database was updated (now v. 7.9, compare with v. 5.5 in [1]) and new models
were developed for testing. The models for predicting cumulative fluid
production now use a stacked approach with Ridge Regression and CatBoost
trained on the residuals to get smooth dependencies of target values relative
to design parameters.

Field test were performed with 21 wells. The accuracy of our model in pilot
tests ($MAPE=37.28\%$, $wMAPE=27.46\%$) is close to the hold-out set accuracy
check while model training ($MAPE=36.08\%$, $wMAPE=29.06\%$).

A method for similar wells search was developed. The method uses a clustering
technique and the euclidean distance as a metric of similarity. The results
help engineers to look back at previously conducted treatments. This method is
also used to estimate the interval limits for the design parameters to be
optimised. Lower and upper limits are $5^{th}$ and $95^{th}$ percentiles of
the cluster parameter values respectively. One can reduce the size of this set
of wells, leaving only top-N wells by the euclidean distance. This allows us
to look for optimal values of the design parameters in a certain vicinity
where the prediction model works well. Also, the top-N similar wells can be
used for imputing missing parameters for the pilot well.

A number of optimization techniques were tested during the pilot testing:
differential evolution, sequential least squares programming, particle swarm
optimization and surrogate-based optimization.

The SBO approach on average maximized the 3-month fluid production better than
other methods. General optimization trend of 21 pilot wells is to increase the
amount of fracturing fluid and proppant mass. Trends for final, average
proppant concentration and fluid rate are different for each type of well
(horizontal, vertical and vertical multilateral). The calculated production
from the fracturing with the optimal set of design parameters, compared to the
treatments with an actual HF design, gives a theoretical target improvement of
38%.

In future work, we plan the following extensions of the presented workflow:
(i) extend the list of features to include presence of upper/lower water-
bearing layers to be able to predict separately the production of pure oil and
water, and extend the target to maximum total fluid and maximum pure oil (or
minimum water cut); (ii) implement economics model to be able to work under
the metrics in terms of Q/CAPEX; (iii) consider the influence of injection
wells on production rates; (iv) extend the target further to be able to treat
the problem as a multi-criteria optimization where the goal is to maximize the
production and to simultaneously minimize the total proppant load.
Furthermore, it can be very useful to combine synthetic data from fracture
design simulations (i.e, using commercial simulators). With primary fracture
design and reservoir properties data, we could simulate fracture propagation
and obtain fracture length, width and height, which directly affect
production. This should improve the predictive power of our models.

Finally, we plan to perform optimization for each well in the database to
generalize the overall optimization trend.

## Acknowledgements

The authors are grateful to the management of LLC “Gazpromneft-STC” for
organizational and financial support of this work. The authors are
particularly grateful to M.M. Khasanov, A.A. Pustovskikh, I.G. Fayzullin and
A.S. Margarit for organizational support of this initiative. The help from
P.K. Kabanova and A.R. Mukhametov in data gathering is greatefully
appreciated.

We are grateful to the management of the DATADVANCE company, namely, to S.M.
Morozov, and to the engineers of the company, A. Saratov and Yu. Bogdanova,
for the scientific support with Surrogate Based Optimization methods from the
pSeven Core Suite, a product of DATADVANCE.

We would like to state it explicitly that the models presented in this work
are solely based on the field data provided by JSC Gazpromneft and we are
grateful for the permission to publish.

## References

  * [1] A. D. Morozov, D. O. Popkov, V. M. Duplyakov, R. F. Mutalova, A. A. Osiptsov, A. L. Vainshtein, E. V. Burnaev, E. V. Shel, G. V. Paderin, Data-driven model for hydraulic fracturing design optimization: Focus on building digital database and production forecast, Journal of Petroleum Science and Engineering 194 (2020) 107504.
  * [2] A. Azbukhanov, I. Kostrigin, K. Bondarenko, M. Semenova, I. Sereda, D. Yulmukhametov, et al., Selection of wells for hydraulic fracturing based on mathematical modeling using machine learning methods (russian), Oil Industry Journal 2019 (11) (2019) 38–42.
  * [3] L. Xue, Y. Liu, Y. Xiong, Y. Liu, X. Cui, G. Lei, A data-driven shale gas production forecasting method based on the multi-objective random forest regression, Journal of Petroleum Science and Engineering 196 107801.
  * [4] A. Davtyan, A. Rodin, I. Muchnik, A. Romashkin, Oil production forecast models based on sliding window regression, Journal of Petroleum Science and Engineering 195 (2020) 107916.
  * [5] J. Yao, Z. Li, L. Liu, W. Fan, M. Zhang, K. Zhang, Optimization of fracturing parameters by modified variable-length particle-swarm optimization in shale-gas reservoir, SPE Journal (2021) 1–18.
  * [6] I. Sobol, Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Mathematics and Computers in Simulation doi:10.1016/S0378-4754(00)00270-6.
  * [7] Lundberg, Scott M., and Su-In Lee, A unified approach to interpreting model predictions, Advances in Neural Information Processing Systems.
  * [8] L. S. Shapley, A value for n-person games, Contributions to the Theory of Games 2 (28) (1953) 307–317.
  * [9] E. Song, B. L. Nelson, J. Staum, Shapley effects for global sensitivity analysis: Theory and computation, SIAM/ASA Journal on Uncertainty Quantification 4 (1) (2016) 1060–1083.
  * [10] A. Erofeev, D. Orlov, D. Perets, D. Koroteev, AI-Based Estimation of Hydraulic Fracturing Effect, SPE Journal doi:10.2118/205479-PA.
  * [11] Stephan Kolassa and Wolfgang Schütz, Advantages of the MAD/Mean Ratio over the MAPE, Foresight: The International Journal of Applied Forecasting (6) (2007) 40–43.
  * [12] R. E. Bellman, S. E. Dreyfus, Applied dynamic programming, Vol. 2050, Princeton university press, 2015.
  * [13] Z. Fu, G. Liu, L. Guo, Sequential Quadratic Programming Method for Nonlinear Least Squares Estimation and Its Application, Mathematical Problems in Engineering 2019 (2019) 1–8.  doi:10.1155/2019/3087949.
  * [14] M. R. Bonyadi, Z. Michalewicz, Particle swarm optimization for single objective continuous space problems: A review, Evolutionary Computation 25 (1) (2017) 1–54.  doi:10.1162/EVCO_r_00180.
  * [15] R. Storn, K. Price, Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces, Journal of Global Optimization 11 (1997) 341–359.  doi:10.1023/A:1008202821328.
  * [16] A. S. Mohaghegh S, Balan B, A hybrid, neuro-genetic approach to hydraulic fracture treatment design and optimization.
  * [17] E. T. Woldemariam, H. G. Lemu, A machine learning based framework for model approximation followed by design optimization for expensive numerical simulation-based optimization problems, in: Proceedings of the Twenty-ninth International Ocean and Polar Engineering Conference www.isope.org Honolulu, Hawaii, USA, June 16-21, International Society of Offshore and Polar Engineers (ISOPE), 2019.
  * [18] S. Wang, S. Chen, Insights to fracture stimulation design in unconventional reservoirs based on machine learning modeling, Journal of Petroleum Science and Engineering.
  * [19] J. Shi, J. Song, B. Song, W. F. Lu, Multi-objective optimization design through machine learning for drop-on-demand bioprinting, Engineering 5 (3) (2019) 586 – 593.  doi:https://doi.org/10.1016/j.eng.2018.12.009.
  * [20] M. Kaya, S. Hajimirza, Using a novel transfer learning method for designing thin film solar cells with enhanced quantum efficiencies, Nature.
  * [21] M. Belyaev, E. Burnaev, E. Kapushev, M. Panov, P. Prikhodko, D. Vetrov, D. Yarotsky, Gtapprox: Surrogate modeling for industrial design, Advances in Engineering Software 102 (2016) 29 – 39.
  * [22] DATADVANCE© pSeven, Release 6.20, User Manual, 2021.
URL https://www.datadvance.net/product/pseven/

  * [23] E. Burnaev, M. Panov, Adaptive design of experiments based on gaussian processes, in: A. Gammerman, V. Vovk, H. Papadopoulos (Eds.), Statistical Learning and Data Sciences, Springer International Publishing, Cham, 2015, pp. 116–125.
  * [24] E. Burnaev, M. Panov, A. Zaytsev, Regression on the basis of nonstationary gaussian processes with bayesian regularization, Journal of Communications Technology and Electronics 61 (6) (2016) 661–671.
  * [25] A. Zaytsev, E. Burnaev, Large scale variable fidelity surrogate modeling, Annals of Mathematics and Artificial Intelligence 81 (1) (2017) 167–186.
  * [26] E. Burnaev, A. Zaytsev, Surrogate modeling of multifidelity data for large samples, Journal of Communications Technology and Electronics 60 (12) (2015) 1348–1355.

## Appendix

Figure 17: CatBoost (green) vs Stacked (Ridge+CatBoost) dependences Figure 18:
Example of offset wells selection Figure 19: Real vs Predicted production on
the real design Figure 20: Production increase with optimal parameters set
Figure 21: Optimized production comparison between methods


# Application of Deep Q Learning with Simulation Results for Elevator
Optimization

Zheng Cao1, Raymond Guo2, Caesar M. Tuguinay3,
Mark Pock4, Jiayi Gao5, Ziyu Wang6 University of Washington, Seattle, USA
Department of Mathematics<EMAIL_ADDRESS>Department of Computer Science & Math
<EMAIL_ADDRESS>Department of Mathematics<EMAIL_ADDRESS>Department of
Computer Science<EMAIL_ADDRESS>Academy for Young Scholars
<EMAIL_ADDRESS>Department of Economics<EMAIL_ADDRESS>
###### Abstract

This paper presents a methodology for combining programming and mathematics to
optimize elevator wait times. Based on simulated user data generated according
to the canonical three-peak model of elevator traffic, we first develop a
naïve model from an intuitive understanding of the logic behind elevators. We
take into consideration a general array of features including capacity,
acceleration, and maximum wait time thresholds to adequately model realistic
circumstances. Using the same evaluation framework, we proceed to develop a
Deep Q Learning model in an attempt to match the hard-coded naïve approach for
elevator control. Throughout the majority of the paper, we work under a Markov
Decision Process (MDP) schema, but later explore how the assumption fails to
characterize the highly stochastic overall Elevator Group Control System
(EGCS).

Keywords: Deep Q Learning, Optimization, Simulation, Markov Decision Process,
Temporal Difference, Elevator Group Control Systems

###### Contents

  1. 1 Introduction
    1. 1.1 Framing and Literature
    2. 1.2 Strategy
  2. 2 Theoretical Background
    1. 2.1 Reinforcement Learning, Markov Decision Processes
    2. 2.2 Q-Learning
  3. 3 Simulation
    1. 3.1 Individual Person Values
    2. 3.2 Office Building Values
    3. 3.3 Simulated Table
  4. 4 Modeling
    1. 4.1 Elevator
    2. 4.2 Environment
    3. 4.3 Model Interactions
    4. 4.4 Design
    5. 4.5 Naïve Approach to Elevator Control:
  5. 5 Modeling via Deep Q Network
    1. 5.1 Motivations
    2. 5.2 EGCS Data Encoding and Decoding
    3. 5.3 DQN Action Model
  6. 6 Results
    1. 6.1 Model Results
    2. 6.2 Model Failure and Why The Model Failed
    3. 6.3 Possible Solutions
  7. 7 Future Work
  8. 8 Conclusion
  9. 9 Acknowledgement

##  1 Introduction

Elevators figure strongly into the daily life of the ordinary urbanite. The
role of elevator wait times may seem negligible, but this comes as a result of
decades of optimizations and improvements in EGCS facilitated by a wide array
of fields. The minimization of elevator wait times becomes especially crucial
during the down-peak and up-peak, when many elevators are crowded. Poor
algorithmic design results in frustrated and tired workers cramming around
elevator doors, wasting valuable time.

This research group has approached elevator optimization via two angles,
addressed by two separate teams – explicit mathematical modelling and machine
learning for approximate optimization. This paper focuses on the latter
approach. Our team’s source code is contained in the “Elevator Project”
GitHub. [4] We first generate data according to the canonical three-peak
model, and use it to build a base-case model to analyze the performance of
traditional elevator design. We subsequently turn towards Deep Q Learning to
attempt optimization over the naïve base-case.

###  1.1 Framing and Literature

We frame our discussion of elevator optimization through the hierarchical
paradigm of Elevator Group Control Systems (EGCS), the central mechanisms in
multi-elevator buildings which control and monitor elevator motion. Where
elevators stay by default, which elevators will be dispatched to various hall
calls, etc. are managed by EGCS. Its importance to internal transportation has
led to an array of research where innovations from across engineering
disciplines have been combined and synthesized to produce the modern elevator
– a far cry from the early days of elevators which ran on schedules.
Throughout the paper, we interpret the idea of EGCS in a way abstracted from
its material implementations – as the overall state, algorithm for responding
to that state, and state transition functions in a particular building.

Several authors before us have attempted to apply machine learning to EGCS –
particle swarm optimizations [1], Convolutional Neural Networks (CNNs), and
neuro-fuzzy systems [5], amongst other approaches. Combining machine learning
approaches with more rigorous mathematical approaches has been a particularly
fruitful approach to other problems – for example, Zheng Cao’s previous paper
“Application of Convolutional Neural Networks with Quasi-Reversibility Method
Results for Option Forecasting.” [3]

###  1.2 Strategy

Our approach to optimization is straightforward – an application of Deep Q
Learning to what we characterize as a classification problem for an arbitrary
decision algorithm taking a simplified version of a building’s current state
and outputting commands to the elevator(s) therein.

##  2 Theoretical Background

###  2.1 Reinforcement Learning, Markov Decision Processes

The generic Reinforcement Learning (RL) problem is framed as an interaction
between an agent and an environment. At each time-step, the agent selects an
action out of a set of possibilities. The environment responds by shifting to
a different state and presenting that state to the agent alongside a scalar
reward. This interaction continues until the environment reaches a terminal
state (there are also RL problems involving environments without a terminal
state, but the problem discussed here is not one of them). A complete sequence
of actions from the agent and responses from the environment, from start to
terminal state, is known as an “episode”. We denote the $n$th state by
$S_{n}$, the $n$th action by $A_{n}$, and the $n$th reward by $R_{n}$ (where
$S_{0}$ is the initial state, $A_{0}$ is the first action, and $R_{0}$ is the
reward given in response to that action).

In this sense, RL problems can be thought of as a series of classification
problems where the agent, at every time step $n$, is tasked with choosing the
action that maximizes some function of the rewards following the $n$th action.
This function is usually (and is here),

$\sum_{i=n}^{t}\lambda^{i-n}S_{i}$

known as the discounted return, where $t$ is the time step after which the
terminal state is reached. $0\leq\lambda<1$ is the “discount factor” and
chosen as a hyperparameter in training, where lower values will make the model
prioritize increasing immediate rewards, and higher values will make the model
have a more “long term” value.

A finite Markov Decision Process (MDP) is a special case of these RL problems
where the number of states and possible rewards is finite, the number of
actions that can be chosen in response to each state is finite, and most
importantly, the probability of any state-reward pair given in response to any
action and previous state is dependent only on that action and previous state
(and not any of the actions and states that preceded them). This is known as
the “Markov Property” and can be expressed symbolically by

$\text{Pr}\\{R_{t+1}=r,S_{t+1}=s|S_{0},A_{0},R_{1},...,S_{t-1},A_{t-1},R_{t},S_{t},A_{t}\\}=\text{Pr}\\{R_{t+1}=r,S_{t+1}=s|S_{t},A_{t}\\}$

for any $r,s$ that lie in the set of possible rewards and states respectively.
A more formal definition is given in the appendix. MDPs are important because
most proofs that provide guarantees that RL methods converge are based on the
MDP case, although empirically great success has been achieved in applying
these methods to non-MDP RL problems. Its importance will be expounded upon in
a later section. [8]

###  2.2 Q-Learning

Finite MDPs are solved by finding a good policy $\pi(a|s)$, which is a
probabilistic function that defines the actions of the agent. In particular,
$\pi(a|s)$ denotes the probability of the agent choosing action $a$ next, if
the last state were to be $s$. With respect to such a policy $\pi$, we can
define a value function $v_{\pi}(s)$, which defines the expected return of an
agent starting at state $s$ and following policy $\pi$. Similarly, we can
define an action-value function $q_{\pi}(s,a)$ which defines the expected
return of an agent starting at state $s$, taking action $a$, and then
following policy $\pi$ thereafter. In finite MDPs, there is at least one
optimal policy [8] $\pi_{*}$, whose value function $v_{*}$ has the property
that for any other value function $\pi$, $v_{*}(s)\geq v_{\pi}(s)$ for all
states $s$. Its action-value function is denoted $q_{*}$. The optimal action-
value function must satisfy the Bellman equations [8], shown below:

$\forall\text{ states
}s,q_{*}(s,a)\mathbb{E}\left[R_{0}+\gamma\max_{a^{\prime}}q_{*}(S_{1},a^{\prime})|S_{0}=s,A_{0}=a\right]$

In fact, the optimal action-value function is the only action-value function
capable of satisfying this. [8]

It essentially states that the action-value function must be the expected
value of the first reward provided added to the action-value function applied
to the next state and best action given that state.

Many RL methods try to approximate solutions to this equation through
alternating policy evaluation (computing the action-value function of the
policy you currently have) and policy improvement (actually making the policy
better).

Early versions of Monte-Carlo methods performed policy evaluation by computing
a collection of episodes and approximating the action-value function at state
$s$ and action $a$ by the average of the return from all episodes passing
through $s$ using $a$, and performed policy improvement by changing the policy
to pick $\pi(s)$ to be the argmax of $q(s,a)$ over all possible actions $a$,
where $q$ is the approximate action-value function found in policy evaluation.

These methods have the problem of only updating predictions after (possibly
extremely long) episodes complete, so they are improved upon by Q-Learning
[8], which instead, in policy evaluation, uses the update rule

$Q(S_{t},A_{t})\longleftarrow
Q(S_{t},A_{t})+\alpha[R_{t+1}+\lambda\max_{a}Q(S_{t},a)-Q(S_{t},A_{t})]$

where $R_{t+1}$ is the actual reward given by the environment, $Q$ is our
current estimate for the action-value function, and $\alpha$ is a
hyperparameter which decides how heavily new data should be weighted. It has
been proven that Q-learning on an MDP converges almost surely to the optimal
value function. [7] This method is still impractical because they require
storing estimates for $Q(s,a)$ for every state-action pair $s$, $a$, or at
least the ones that are often seen in practice. In our case, and many others,
there are simply far too many state-action pairs for this to be done.

As a result, we instead rely on a method first introduced in a paper using RL
to play Atari [6] which utilizes a Neural Network to evaluate $Q$. The state,
actions, reward, and next state from each time-step are stored in a queue, and
every time a finite number of time-steps have occurred, the Neural Network is
trained on data in the queue with loss function

$Q(S_{n},A_{n})-(r_{n}+\max_{a^{\prime}}Q(S_{n+1},a^{\prime}))$

where $Q$ is the current action-value function estimate, $S_{n}$ is the
current state, $A_{n}$ is the current action, $r_{n}$ is the reward given for
that action, and $S_{n+1}$ is the next state. We note that
$(r_{n}+\max_{a^{\prime}}Q(S_{n+1},a^{\prime}))$ is the update value from
Q-learning with $\alpha$ set to 1. The policy we follow is dictated by a
second neural network, and policy improvement is performed by updating the
second neural network’s weights with a weighted average of weights from both
networks, essentially “nudging” the policy in the right direction. More
details are shown below.

##  3 Simulation

Before tackling EGCS, we seek to develop a system for generation of user
traffic data to feed into the system. We modularize the generation into
individual people who select times to make hall calls from a truncated Poisson
distribution.

We chose to simulate data rather than directly observe it in the real world in
order to have the flexibility to work with an arbitrary amount of people while
still accounting for the three-peak model. Moreover, we are able to adjust the
characteristics of the building with which we are working rather than being
bound to observed characteristics. We also save a considerable amount of time
and energy.

###  3.1 Individual Person Values

To have the widest application of our model, we generated raw data based on
the standard 9-to-5 work schedule with a 30 minute lunch break in between. We
assumed that all the workers would arrive and leave the building with a mean
time of 32400 seconds (9AM) and 61200 seconds (5PM) with a standard deviation
of 1800 seconds (30 minutes). To simulate real life, we also randomly
generated each person’s weight based on a normal distribution based off of an
individual’s sex.

###  3.2 Office Building Values

Similarly to individual person values, we simulated office building values for
simulation; the following describe our office building values and are based
off of aggregate personal observation: There are 8 floors in the building.
There are approximately 200 workers per elevator in each building. It takes
approximately 15 seconds for an elevator to open, load people, and close. It
takes approximately 5 seconds for an elevator to ascend or descend a single
floor.

###  3.3 Simulated Table

Table 1. Sample Simulated User Data

time (second) | start floor | destination floor) | weight
---|---|---|---
… | … | … | …
27000 | 1 | 2 | 77.9…
… | … | … | …
27042.9… | 1 | 7 | 78.3…
… | … | … | …
60754.1… | 8 | 1 | 101.6…
… | … | … | …

##  4 Modeling

###  4.1 Elevator

We formulate a basic representation of a single-elevator control system which
will serve as the later basis for our expansions to more complex EGCS.

Using one elevator is both the easiest to mathematically model and the easiest
to attempt to bring out an optimal method. Therefore, we decided to use a
single elevator to evaluate a base-line value.

When there are more elevators, say $n$ elevators, a very simple threshold to
compare our results with would simply be to see if our models can improve the
time such that it is better than base-line / n seconds, which would be the
result if we were to distribute the people evenly across all elevators.

###  4.2 Environment

The environment consists of an object called a Time List containing all the
events, an object called a State which describes the certain aspects of the
simulation.

In particular, the State consists of the upward and downward calls of
passengers looking to enter the elevator and go a specific direction, the time
of the simulation, the elevator’s speed, the wait time that it takes to open
and close the door to let passengers in, and the elevator.

The environment details will be expounded upon in the following naïve models
we will detail.

###  4.3 Model Interactions

Attached to the environment is the model which decides an action based off of
each State. Each action will update the State and Time List accordingly;
however, the State may be updated prior to each action based off of events
within the Time List.

At the start of the simulation, we initialize the state so that the elevator
starts at the first floor and the state’s time starts at $27000$ seconds (7:30
AM). During the simulation, we iterate through the Time List until there are
no more events within the Time List and when there are no more passengers
within the State that wish to be moved. When this occurs, the simulation ends
and we calculate the total time all it took for all the passengers to reach
their destinations.

The model takes actions based off of the State; however, it only gets passed 3
objects from the State:

> 1\. The state of the up / down button of each floor – Whether it is pressed
> or not.

> 2\. The state of the buttons in the elevator – whether it is pressed or not.

> 3\. The current floor of the elevator.

###  4.4 Design

Our system centers around the management of a global state variable by a main
event loop which processes events in a Time List, an ordered procedure for all
the relevant events in the period over which the simulation occurs – for us a
day – generated beforehand.

To be specific, we determine before the start of a simulation when people will
enter and exit the building and which floors they will want to go to using the
classical three-peak model and truncated Poisson distributions around each
peak. The various hall calls thereby generated are fed into our Time List to
be processed by the event loop as the simulation progresses.

Treating each Time List event as a stimulus, we manage the current State and
request an action from the model, passing in our State as a parameter, to
precipitate the succeeding state. Any model we use for this single-elevator
system has a simple choice between three commands – Idle, Move, and Open
Doors. This classification problem forms the premise of the optimization we
attempt to achieve. Our naïve model deals with State in a hard coded and
’common sense’ way – one which is of course deeply inapplicable to more
complex EGCS, but allows us to verify the functionality of the system.

See Naïve Model with Environment 1.0.

###  4.5 Naïve Approach to Elevator Control:

We created a basic algorithm for elevator control that approximates how many
elevators practically function. The algorithm alternates between two phases.
In phase 1, there are no passengers in the elevator, and the elevator travels
to the closest floor with a hall call (if there are no hall calls at all, the
elevator will instead idle). Once the elevator travels to a floor with a hall
call, it will randomly choose one of the (at most two) directions in which
hall calls have been made on that floor, open its doors, and signal its intent
to travel in that direction.

This begins phase 2, where the elevator will continuously travel in that
direction, letting people off at their intended destination and picking people
up when it reaches a floor with hall calls in the same direction it is
travelling in. This phase continues until there are no remaining passengers on
the elevator, in which case phase 1 will start again. We tested the results of
this model on both the 1.0 (see image Naïve Model with Environment 1.0) and
2.0 (see image Naïve Model with Environment 2.0) versions of the naïve model.

##  5 Modeling via Deep Q Network

To further optimize total run time, we set a Deep Q Networks (DQN) to be the
EGCS’s inner decision making model.

###  5.1 Motivations

A DQN model was chosen to be the machine learning model to optimize the
environment because there exist numerous examples of successful DQN decision
making models trained and tested in similar environments to our EGCS.

The original DQN Atari Paper [4] has shown that DQNs can be trained on
environments with large continuous State Spaces and discrete Action Spaces.
The EGCS State Space is mildly large and continuous, and the Action Space is
discrete.

###  5.2 EGCS Data Encoding and Decoding

The following describe how data flows in and out of the inner decision making
DQN model within the EGCS.

Algorithm 1 *

Input Encoding

Get capacity, current weight, current position, up buttons, down buttons,
buttons pressed from State.

tensor = $[\text{cacapacity)},\text{float(current weight)},\text{float(current
position)}]$

for entry in buttons pressed (8 entries) do

append float(entry) to tensor

end for

for entry in up buttons (8 entries) do

append float(entry) to tensor

end for

for entry in down buttons (8 entries) do

append float(entry) to tensor

end for

(Note: float(boolean) $==0$ for false and float(boolean) $==1$ for true)

Algorithm 2 *

Output Decoding

Get Action Value: Natural Number between 0 and 4.

if Action Value == 0 then

Action = Idle

else if Action Value == 1 then

Action = Open Close Doors for Up Calls

else if Action Value == 2 then

Action = Open Close Doors for Down Calls

else if Action Value == 3 then

Action = Move Up

else if Action Value == 4 then

Action = Move Down

end if

###  5.3 DQN Action Model

The following describes how the EGCS DQN Action Model is trained.

Algorithm 3 *

Training the EGCS Action Taking DQN

Set $\epsilon$ and $\lambda$ values: Number between 0 and 1.

Initialize Sample Hall Calls (SHC).

Initialize Updating $Q^{*}$ Steps as $C$.

Initialize Replay Memory as D.

Initialize Action Value Model $Q$ with weights $W$.

Initialize Target Action Value Model $Q^{*}$ with weights $W^{*}=W$.

Set epoch range as RANGE.

for Epoch in RANGE do

Initialize Time List and Initialize State $S_{t}$ with SHC.

Total Time = 0, Number Steps = 0.

while no more events in Time List do

if action can be taken then

if $\epsilon>$ random number between 0 and 1 then

$A_{t}$ = Choose Random Action.

else

$A_{t}=\text{ArgMax}_{A}Q(S_{t},A,W)$.

end if

Take action $A_{t}$ and observe $Q_{\text{Val}}$ = $Q(S_{t},A,W)$, and
$S_{t+1}$.

end if

Update Total Time, Number Steps, and $S_{t}$.

$R_{t}$ based off of environment: $R_{t}=-1*$ number of people waiting * sum
added time.

Store transition $(S_{t},A_{t},R_{t},S_{t+1})$ in D.

Sample Random Minibatch of Transitions $(S_{j},A_{j},R_{j},S_{j+1})$ from D as
Mini.

Minibatch Count $=0$.

for $(S_{j},A_{j},R_{j},S_{j+1})$ in Mini do

$M_{j}=\lambda*\text{max}_{A}Q^{*}(S_{j+1},A,W^{*})$ or $0$ if $S_{j+1}$ is
terminal. Then set $Y_{j}=R_{j}-M_{j}$.

Loss = $(Y_{j}-Q(S_{j},A_{j},W))^{2}$. Perform a Gradient Descent Step with
Loss Value, with respect to $Q(W)$. Minibatch Count $+=1$.

if Minibatch Count % $C==0$ then

$Q^{*}(W^{*})=Q(W)$.

end if

end for

end while

end for

The following describes how the EGCS DQN Action Model is used for inference.

Algorithm 4 *

Inference with the EGCS Action Taking DQN

Initialize Sample Hall Calls (SHC).

Initialize Action Value Model $Q^{*}$ with trained weights $W^{*}$.

Initialize Time List and Initialize State $S_{t}$ with SHC.

Total Time = 0, Number Steps = 0.

while no more events in Time List do

if action can be taken then

$A_{t}=\text{ArgMax}_{A}Q(S_{t},A,W)$.

Take action $A_{t}$ and observe $Q_{\text{Val}}$ = $Q(S_{t},A,W)$, and
$S_{t+1}$.

end if

Update Total Time, Number Steps, and $S_{t}$.

end while

##  6 Results

###  6.1 Model Results

Figure 1: *

Naïve Model with Environment 1.0

Figure 2: *

Naïve Model with Environment 2.0

Figure 3: *

Output logs displaying the training of the DQN Model.

###  6.2 Model Failure and Why The Model Failed

As seen by the above output logs, there is no improvement in Number of Events,
Number of People Moved, Mean Total Time, Median Total Time, Max Total Time,
Sum Total Time, when examining the model over ten epochs/full simulation
cycles; almost all improvements were abrupt and did not follow a trend of
improvement.

The EGCS we are working with is not an MDP. Observe the following: Take any
state $s$ such that the immediate state after, $s^{{}^{\prime}}$, contains a
Hall Call not previously introduced into the EGCS. The Hall Call is placed
within the EGCS and the state $s$ is updated to $s^{{}^{\prime}}$ at random
with no effect from action $a$. Then, there does not exist a function $P_{a}$
for $s,s^{{}^{\prime}},a$ as described by Definition 9.3. Therefore, the EGCS
is not an MDP.

We know that such states $s,s^{{}^{\prime}}$ exists by examining the
simulation at the beginning of time: Time $27000$ seconds (7:30 AM). Progress
through the simulation by taking a set of sequential actions
$A=\\{a_{1},...,a_{n}\\}$ with each action being taken on a set of sequential
states $S=\\{s_{1},...,s_{n}\\}$ until we reach a single Hall Call $h$ and no
$s\in S$ contains $h$. The Hall Call will then be placed in the next state:
$s^{{}^{\prime}}$.

An MDP is necessary for any analytic assurance of convergence within
Q-learning. [7] DQN is an implementation of Q-learning where we approximate
the reward function with a neural network. Extrapolating from [7], there is no
assurance for convergence within the DQN model proposed in this paper, from
the perspective of the paper from [7]. Additionally, the authors of the
original DQN Atari Paper [6] said that their method lacks any convergence
guarantees.

###  6.3 Possible Solutions

To fix the problem of having no analytic assurance of convergence, a modified
EGCS where the inner decision making model has knowledge on all the future
Hall Calls can be created.

The complication with this is that the EGCS no longer models real-life
scenarios: Elevators should not know when and how many people will make Hall
Calls. Therefore, this model will not be useful outside of simulations.

A possible fix to this modified EGCS inner decision making model previously
described is to create a model which will predict future Hall Calls. This will
attach to the inner decision making model and will feed it predicted future
Hall Calls. It could back-propagate and update its weights based off of the
Hall Calls the EGCS actually receives during each run.

##  7 Future Work

The largest piece of future work we would like to see is a more successful RL-
based model that offers improvements over the Naïve algorithm we created. We
would also like to see the environment expanded to handle multi-elevator
systems and create a more generalized Naïve algorithm that works on these
multi-elevator systems.

##  8 Conclusion

In this project, we created a simulated environment to test the effectiveness
of a given elevator control algorithm and demonstrated its correctness with a
naïve control algorithm that achieved expected results. In the future, this
environment can be used to test various other Reinforcement Learning
approaches to our elevator optimization problem. We also attempted to employ a
Deep Q-Learning approach to create a ML model to improve upon the naïve
algorithm discussed above. This goal will be addressed more completely in
future work, and we hope that the environment we created can be used as
groundwork for future attempts based on Reinforcement Learning to solve the
same problem.

##  9 Acknowledgement

The original inspiration for elevator optimization came from Zheng Cao and his
experiences with the antiquated elevators on the University of Washington
campus. From this catalyst, we have attempted to approach elevator
optimization using reinforcement learning and explicit mathematical modelling.
Our second approach, facilitated by our math research team, we are currently
outlaying in a separate paper, “Application of Spatial Process and Gibbs
Random Field Approaches for Dumbwaiter Modeling for Elevator Optimization,”
which applies methods from pure math such as Spatial Process and Gibbs Random
Field approaches to optimize elevator wait times. [2]

Both projects produce independent results, and potential combinations of both
approaches may help achieve a better optimization in future research.
Nevertheless, we acknowledge the help the main authors of the additional paper
have serviced in creating this reinforcement learning approach: Wanchaloem
Wunkaew, Xiyah Chang, and Benjamin Davis.

## References

  * [1] B Bolat, O Altun, P Cortés, _A particle swarm optimization algorithm for optimal car-call allocation in elevator group control systems_. Applied Soft Computing, 2013.
  * [2] Z Cao, X Chang, B Davis, W Wunkaew, _Elevator Optimization: Application of Spatial Process and Gibbs Random Field Approaches for Dumbwaiter Modeling and Multi-Dumbwaiter Systems_ , https://arxiv.org/abs/2209.12401
  * [3] Z Cao, W Du, K V. Golubnichiy, Application of Convolutional Neural Networks with Quasi-Reversibility Method Results for Option Forecasting, https://arxiv.org/pdf/2208.14385.pdf.
  * [4] C M. Tuguinay, R Guo, M Pock, J Gao, Z Wang, Z Cao, _Elevator Project_ , https://github.com/ctuguinay/ElevatorProject. Source code, 2022.
  * [5] Y L. Dong, X M. Wang, _The Intelligent Scheduling Method of Elevator Group Control Based on Fuzzy Neural Network_. Applied Mechanics and Materials, 2012.
  * [6] V Mnih, K Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, Martin Riedmiller, _Playing Atari with Deep Reinforcement Learning,_ https://arxiv.org/pdf/1312.5602.pdf.
  * [7] M T. Regehr, A Ayoub, An Elementary Proof that Q-learning Converges Almost Surely, https://arxiv.org/pdf/2108.02827.pdf
  * [8] R Sutton, A Barto, _Reinforcement Learning: An Introduction,_ https://inst.eecs.berkeley.edu/ cs188/sp20/assets/files/SuttonBartoIPRLBook2ndEd.pdf.

## Appendix

### Markov Decision Processes

A Markov Decision Process (MDP) is defined as a four-tuple $(S,A,P_{a},R_{a})$
with the following properties:

######  Definition 9.1.

S is a set of states, which is called a State Space.

######  Definition 9.2.

A is a set of actions, called the Action Space. The Action Space is (usually)
derived from S.

######  Definition 9.3.

$\forall a\in A$ and $\forall s,s^{{}^{\prime}}\in S$,
$P_{a}(s,s^{{}^{\prime}})=Pr(s_{t+1}=s^{{}^{\prime}}|s_{t}=s,a_{t}=a)$; the
probability of taking action $a$ at state $s$ and transitioning to state
$s^{{}^{\prime}}$.

######  Definition 9.4.

$\forall a\in A$ and $\forall s,s^{{}^{\prime}}\in
S,R_{a}(s,s^{{}^{\prime}},a)$ is the immediate reward received after taking
action $a$ and transitioning from state $s$ to state $s^{{}^{\prime}}$.


possible mitigation strategies (e.g., gated release of models, providing
defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms
to monitor how a system learns from feedback over time, improving the
efficiency and accessibility of ML).

  51. 11.

Safeguards

  52. Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?

  53. Answer: [N/A]

  54. Justification: no models or data are in this paper.

  55. Guidelines:

     * •

The answer NA means that the paper poses no such risks.

     * •

Released models that have a high risk for misuse or dual-use should be
released with necessary safeguards to allow for controlled use of the model,
for example by requiring that users adhere to usage guidelines or restrictions
to access the model or implementing safety filters.

     * •

Datasets that have been scraped from the Internet could pose safety risks. The
authors should describe how they avoided releasing unsafe images.

     * •

We recognize that providing effective safeguards is challenging, and many
papers do not require this, but we encourage authors to take this into account
and make a best faith effort.

  56. 12.

Licenses for existing assets

  57. Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?

  58. Answer: [N/A]

  59. Justification: no code or models are used.

  60. Guidelines:

     * •

The answer NA means that the paper does not use existing assets.

     * •

The authors should cite the original paper that produced the code package or
dataset.

     * •

The authors should state which version of the asset is used and, if possible,
include a URL.

     * •

The name of the license (e.g., CC-BY 4.0) should be included for each asset.

     * •

For scraped data from a particular source (e.g., website), the copyright and
terms of service of that source should be provided.

     * •

If assets are released, the license, copyright information, and terms of use
in the package should be provided. For popular datasets,
paperswithcode.com/datasets has curated licenses for some datasets. Their
licensing guide can help determine the license of a dataset.

     * •

For existing datasets that are re-packaged, both the original license and the
license of the derived asset (if it has changed) should be provided.

     * •

If this information is not available online, the authors are encouraged to
reach out to the asset’s creators.

  61. 13.

New Assets

  62. Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?

  63. Answer: [N/A]

  64. Justification: No released assets

  65. Guidelines:

     * •

The answer NA means that the paper does not release new assets.

     * •

Researchers should communicate the details of the dataset/code/model as part
of their submissions via structured templates. This includes details about
training, license, limitations, etc.

     * •

The paper should discuss whether and how consent was obtained from people
whose asset is used.

     * •

At submission time, remember to anonymize your assets (if applicable). You can
either create an anonymized URL or include an anonymized zip file.

  66. 14.

Crowdsourcing and Research with Human Subjects

  67. Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?

  68. Answer: [N/A]

  69. Justification: no human subjects

  70. Guidelines:

     * •

The answer NA means that the paper does not involve crowdsourcing nor research
with human subjects.

     * •

Including this information in the supplemental material is fine, but if the
main contribution of the paper involves human subjects, then as much detail as
possible should be included in the main paper.

     * •

According to the NeurIPS Code of Ethics, workers involved in data collection,
curation, or other labor should be paid at least the minimum wage in the
country of the data collector.

  71. 15.

Institutional Review Board (IRB) Approvals or Equivalent for Research with
Human Subjects

  72. Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?

  73. Answer: [N/A]

  74. Justification: No human subjects

  75. Guidelines:

     * •

The answer NA means that the paper does not involve crowdsourcing nor research
with human subjects.

     * •

Depending on the country in which research is conducted, IRB approval (or
equivalent) may be required for any human subjects research. If you obtained
IRB approval, you should clearly state this in the paper.

     * •

We recognize that the procedures for this may vary significantly between
institutions and locations, and we expect authors to adhere to the NeurIPS
Code of Ethics and the guidelines for their institution.

     * •

For initial submissions, do not include any information that would break
anonymity (if applicable), such as the institution conducting the review.


# An alternating proof of sharp inequalities related with Burnside’s formula

necdet batir  department of mathematics, nevşehir hbv university, nevşehir,
50300 turkey<EMAIL_ADDRESS>
###### Abstract.

We provide an alternating proof of the following sharp inequalities related
with Burnside’s formula for $n!$ proved by the author in [A double inequality
related with Burnside’s, formula, Proyecciones J. Math., Vol. 37, No 1, pp.
55-59, 2018.]

$\sqrt{2\pi}\left(\frac{n+a_{*}}{e}\right)^{n+a_{*}}<n!<\sqrt{2\pi}\left(\frac{n+a^{*}}{e}\right)^{n+a^{*}}\,(n\in\mathbb{N}),$

where the constants $a_{*}=0.428844044...$ and $a^{*}=0.5$ are the best
possible.

###### Key words and phrases:

Burnside’s formula, Stirling’s formula, factorial function, inequalities

###### 1991 Mathematics Subject Classification:

33B15, 26D07

##  1\. introduction

It is of interest to study how the factorial function $n!$ behaves when $n$ is
sufficiently large. The most well known approximation formula to approximate
factorial function is the Stirling’s formula given by

$n!\sim n^{n}e^{-n}\sqrt{2\pi n}.$ (1.1)

The origin of this formula is based on a study of the French mathematician
Abraham de Moivre. In 1733 de Moivre developped the formula

$n!\sim C\cdot\sqrt{n}n^{n}e^{-n},$

where $C$ is a constant. He was unable, however, to evaluate exact numerical
value of this constant; this task befell a Scot mathematician James Stirling
(1692-1770), who found $C=\sqrt{2\pi}$. Formula (1.1) is known as Stirling’s
formula today. This formula has many applications in probability theory,
statistical physics. number theory, and special functions. The most important
known simple approximation formula for factorial function after Stirling’s
formula is Burnside’s formula [3], which is given by

$n!\sim\sqrt{2\pi}\left(\frac{n+\frac{1}{2}}{e}\right)^{n+\frac{1}{2}}.$

It is well known that this formula is more accurate than the Stirling’s
formula. Please refer to [4, 5, 6, 7] and the references therein for related
inequalities. In very new paper the author [2] proved that

$\sqrt{2\pi}\left(\frac{n+a_{*}}{e}\right)^{n+a_{*}}<n!<\sqrt{2\pi}\left(\frac{n+a^{*}}{e}\right)^{n+a^{*}}\,(n\in\mathbb{N}),$

where the constants $a_{*}=0.428844044...$ and $a^{*}=0.5$ are the best
possible. Our aim in this note is to supply a completely different proof of
these inequalities.

##  2\. main result

Our main result is a new proof of the following theorem.

######  Theorem 2.1.

For all $n\in\mathbb{N}$ we have

$\sqrt{2\pi}\left(\frac{n+a_{*}}{e}\right)^{n+a_{*}}<n!<\sqrt{2\pi}\left(\frac{n+a^{*}}{e}\right)^{n+a^{*}},$
(2.1)

where the constants $a_{*}=0.428844044...$ and $a^{*}=0.5$ are the best
possible.

###### Proof.

Let $n$ be a fixed positive integer and $x$ be a non-negative real number. We
define

$f(x):=\sqrt{2\pi}\left(\frac{x+n}{e}\right)^{x+n}.$

Differentiation gives $f^{\prime}(x)=f(x)\log(x+n)>0$. So, $f$ is strictly
increasing for all $x>0$. On the other hand, we have

$f(0)=n^{n}e^{-n}\sqrt{2\pi}<n!$

and

$f(1/2)=\sqrt{2\pi}\left(\frac{n+\frac{1}{2}}{e}\right)^{n+\frac{1}{2}}>n!,$

see [1], that is,

$f(0)<n!<f(1/2).$

Since $f$ is continuous on $(0,1/2)$, the intermediate value theorem yields
that there exist a sequence $(a_{n})$ with $a_{n}\subset(0,1/2)$ such that

$f(a_{n})=n!.$ (2.2)

Since $f$ is continuous and monotonic increasing on $(0,1/2)$, it has an
inverse $f^{-1}$ and $f^{-1}$ is also monotonic increasing on the interval

$(\sqrt{2\pi}e^{-1},\infty).$

We therefore get from (2.2) that

$a_{n+1}=f^{-1}((n+1)!)>f^{-1}(n!)=a_{n},$

which implies that $(a_{n})$ is strictly increasing. Since $(a_{n})$ is
bounded and strictly increasing, it has a limit and we have for all
$n\in\mathbb{N}$

$a_{1}<a_{n}<\lim\limits_{n\to\infty}a_{n}.$ (2.3)

If we take $n=1$ in (2.2) we obtain

$1=\sqrt{2\pi}\left(\frac{a_{1}+1}{e}\right)^{a_{1}+1}$

or taking the logarithm of both side

$a_{1}+1-(a_{1}+1)\log(a_{1}+1)-\log\sqrt{2\pi}=0.$

By the help of computer program Mathematica, we can solve this equation and we
find that $a_{1}=0.428844044...$. Now we want to show that
$\lim\limits_{n\to\infty}a_{n}=1/2$. It is difficult to evaluate it in the
usual way. Instead, we shall follow a different way. From

$n!=\sqrt{2\pi}\left(\frac{n+a_{n}}{e}\right)^{n+a_{n}},$

we get

$\lim\limits_{n\to\infty}\frac{\sqrt{2\pi}\left(\frac{n+a_{n}}{e}\right)^{n+a_{n}}}{n!}=1.$
(2.4)

Using Stirling’s formula this becomes

$\lim\limits_{n\to\infty}\frac{(n+a_{n})^{a_{n}}e^{-a_{n}}\left(1+\frac{a_{n}}{n}\right)^{n}}{\sqrt{n}}=1.$

Since $a_{n}$ is convergent it is easy to see that
$\lim\limits_{n\to\infty}e^{-a_{n}}\left(1+\frac{a_{n}}{n}\right)^{n}=1$, so
that

$\lim\limits_{n\to\infty}\frac{(n+a_{n})^{a_{n}}}{\sqrt{n}}=\lim\limits_{n\to\infty}\left(1+\frac{a_{n}}{n}\right)^{a_{n}}\lim\limits_{n\to\infty}n^{a_{n}-\frac{1}{2}}=1.$

The first limit on the right hand side goes to 1 as $n$ approaches to
infinity, so we have $\lim\limits_{n\to\infty}n^{a_{n}-\frac{1}{2}}=1$, which
is possible only if $\lim\limits_{n\to\infty}a_{n}=\frac{1}{2}$. Let us set
$a_{*}=a_{1}=0.428844044...$ and $a^{*}=0.5$. Then by (2.3) we have
$a_{*}<a_{n}<a^{*}$ for all $n\in\mathbb{N}$. Applying $f$ to each side leads
to

$f(a_{*})<f(a_{n})=n!<f(a^{*}),$

which is equivalent to (2.1)with the possible constants $a_{*}=0.428844044...$
and $a^{*}=0.5$. This completes the proof. ∎

## References

  * [1] N. Batir, Inequalities for the gamma function, Archiv Math. (Basel), 91(2008), 554-563.
  * [2] N. Batir, A double inequality related with Burnside’s, formula, Proyecciones J. Math., Vol. 37, No 1, pp. 55-59, 2018.
  * [3] W. Burnside, A rapidly convergent series for $\log N!$, Messenger Math., 46(1917), 157-159.
  * [4] C-P Chen, On the asymptotic expansions of the gamma function related to the Nemes, Gosper and Burnside formulas, Appl. Math. Comp., 276 (2016) 417–431.
  * [5] D. Lu, X. Wang, A new asymptotic expansion and some inequalities for the gamma function, J. Number Theory 140 (2014) 314–323.
  * [6] Z-H. Yang, J-F Tian, On Burnside type approximation for the gamma function, RACSAM, https://doi.org/10.1007/s13398-019-00651-2, 2019.
  * [7] X. Shi, F. Liu and H. Qu, The Burnside Approximation of the gamma function, Analysis and Applications, V. 08, No. 03, pp. 315-322, 2010.