# Clustering Superconductors Using Unsupervised Machine Learning

B. Roter∗ Department of Physics, The University of Akron, Akron, OH 44325, USA
∗Present address: Program of Applied Physics, Northwestern University,
Evanston, IL 60208, USA  N. Ninkovic  Department of Physics, The University of
Akron, Akron, OH 44325, USA  S.V. Dordevic  Department of Physics, The
University of Akron, Akron, OH 44325, USA<EMAIL_ADDRESS>
###### Abstract

In this work we used unsupervised machine learning methods in order to find
possible clustering structures in superconducting materials data sets. We used
the SuperCon database, as well as our own data sets complied from literature,
in order to explore how machine learning algorithms groups superconductors.
Both conventional clustering methods like k-means, hierarchical or Gaussian
mixtures, as well as clustering methods based on artificial neural networks
like self-organizing maps, were used. For dimensionality reduction and
visualization t-SNE was found to be the best choice. Our results indicate that
machine learning techniques can achieve, and in some cases exceed, human level
performance. Calculations suggest that the clustering of superconducting
materials works best when machine learning techniques are used in concert with
human knowledge of superconductors. We also show that in order to resolve fine
subcluster structure in the data, clustering of superconducting materials
should be done in stages.

###### keywords:

clustering , high temperature superconductors , unsupervised machine learning

††journal: Physica C

##  1 Introduction

In recent years a number of artificial intelligence (AI) studies of
superconductors have been performed [1, 2, 3, 4, 5, 6, 7], using either
machine learning (ML) or deep learning techniques. In some of these studies,
the superconducting critical temperature of known superconductors was
predicted, i.e. calculated. Most of these studies achieved impressive values
of statistical parameters, such as R2 or RMSE. However, none of them have
(yet) led to the discovery of new superconducting materials, which should be
the ultimate goal of AI studies.

It was recently pointed out [8] that AI models achieved high values of
statistical parameters because they used traditional k-fold cross-validation
procedures, in particular leave-one-out cross validation which is its special
case. These cross validation techniques have been shown to result in
unrealistically high values of statistical parameters when applied to highly
clustered data sets [8]. Moreover, they have been known to have problems when
extrapolating to new classes of materials [1]. As an alternative, an approach
called leave-one-cluster-out cross validation was proposed [8], which was
supposed to alleviate the problem. The idea is that instead of searching for
new superconducting materials, AI is used to search for new classes of
superconducting materials. However, the problem is that there is currently no
universally accepted classification of superconductors. There have been
numerous attempts to classify superconductors based on their crystal
structure, physical properties, pairing mechanics, etc. [9, 10] but none of
them are universally accepted.

Over the past century, superconductivity has been found in a variety of
different systems, such as metals, oxides, alloys, organic materials, heavy
fermions, etc. These systems can have vastly different physical, chemical and
all other properties, which makes it difficult for humans to classify. In a
recent study Hirsch, Maple and Marcillio [10] made the most comprehensive
attempt yet to separate superconductors into different classes, based on the
physical mechanism of superconductivity. They identified 32 classes, separated
into three big groups: conventional (with 12 classes), potentially
unconventional (with 9 classes) and unconventional (with 11 classes). However,
some superconductors in this study can be found in several different classes.

In this work we used a ML technique called clustering in order to separate
superconductors into classes (or clusters). Our hope is that ML might be able
to separate superconductors into physically meaningful groups, which would in
turn lead to the discovery of entirely new groups of superconductors.
Additionally, ML techniques have been known to both provide additional
insights and reveal hidden patterns in other groups of materials not
previously known to humans [11].

Our results reveal that data-driven ML techniques can indeed separate
superconductors into groups, achieving and in some cases exceeding human level
performance. However, there are also limitations, which will need to be
addressed in the future. In particular, constructing a more reliable database,
which will undoubtedly result in better ML calculations. Clustering results
obtained on the test data sets we compiled from literature, confirm this
claim. Our calculations indicate that, at least for superconducting materials
data, clustering in stages is the best approach.

##  2 Clustering

Clustering is one of the most common tasks of unsupervised machine learning
[12, 13]. The main goal of unsupervised learning algorithms is to find
patterns and learn meaningful relationships in data in order to describe its
underlying structure. However, in unsupervised machine learning one does not
know the correct output; thus the results of calculations cannot be assessed
in terms of accuracy or correctness. Instead, clustering is an exploratory
technique. One must not assume there is a single ”true” clustering. One should
rather be interested in exploring different clusterings of the same data set
to learn more about it. Using commercially available clustering programs as a
”black box” should not be done, as it may lead to erroneous results and
conclusions.

Before we attempt to cluster superconductors, there are several issues that
need to be addressed, such as the choice of a clustering technique, the number
of clusters, the choice of the distance function, dimensionality reduction and
visualizations. All these issues are discussed below in separate subsections
before clustering results are presented.

###  2.1 Conventional Clustering Methods

We used several conventional clustering methods [12, 13]: k-means, k-medoids,
hierarchical, Gaussian mixtures, DBSCAN, etc. For most of them the number of
clusters is an input parameter, which is selected based on the data set.
DBSCAN is one of the few that does not require the number of clusters to be
specified. Another important factor that affects clustering is the distance
function used to calculate similarly/dissimilarity between data points. We
have tried several standard distance functions such as Euclidean, cityblock,
cosine, Chebychev, etc. In our calculations with superconductors, the cosine
distance function has usually produced the best results.

A procedure for selecting the best clustering method was recently proposed in
Ref. [14]. It is referred to as Iterative Label Spreading, and it allows one
to visually assess the quality of clusters, before any clustering is done. It
was shown to be particularly suitable for noisy, high dimensional data
typically encountered in materials science.

###  2.2 Neural Network–Based Clustering Methods

In recent years artificial neural networks, especially deep neural networks,
have become widely used in a number of disciplines such as image and speech
processing, chemistry, biology, material science, etc [15]. Although neural
networks are used primarily for supervised learning tasks [15] such as
regression and classification, they can also be used for unsupervised
learning, i.e. clustering.

In this work we used a neural network based clustering method known as the
Self-Organizing Map (SOM) algorithm. SOMs perform nonlinear transformations of
a multidimensional data set into a low-dimensional set that retains the
intrinsic topological properties of the original data set. The network is
trained to transform a so-called ”input space”, whose dimension equals the
number of predictors, into a low-dimensional output called ”map space”. The
most important parameter of these calculations is the number of nodes, which
is equivalent to the number of clusters. In addition, one also needs to
carefully choose the topology of the input layer. Two most commonly used
topologies are rectangular and hexagonal. In our calculations with
superconducting materials, SOMs have consistently performed better compared to
conventional clustering algorithms.

###  2.3 The Choice of Predictors

One of the most important issues in machine learning calculations is the
choice of predictors. Similar to what we did in our previous work [6], we used
chemical composition as the only predictor for clustering calculations. We
previously showed that the use of a large number of predictors does not
necessarily lead to any improvements in ML models. We suggested that the
predictors must be carefully selected based on the machine learning task. In
calculations with superconductors, one should use those predictors that are
known (or at least believed) to be closely related to superconductivity, such
as crystal symmetry, electronic band structure, number of valence electrons,
etc. The problem is that these predictors are not systematically reported in
the existing databases. They should be included in future databases, as they
might result in more reliable machine learning calculations.

###  2.4 The Number of Clusters

As mentioned above, clustering is an exploratory technique, and the ”correct”
number of clusters does not exist. Most clustering techniques treat the number
of clusters as a parameter of calculations, which can be adjusted based on the
data set. There are several algorithms one can use in order to estimate the
optimal number of clusters. The most commonly used are: Calinski-Harabasz
criterion, Davies-Bouldin index, silhouette and dendrogram plots [12]. A
dendrogram plot in particular allows one to explore not only cluster, but also
sub-cluster structure in the data. Iterative Label Spreading can also be used
to estimate the number of clusters [14]. However, in our experience, the best
way to estimate the optimal number of clusters is by combining machine
learning calculations with human knowledge of superconducting materials. This
approach can produce the most physically meaningful clustering.

Another important issue that one must be aware of is outliers in the data set
used. In most AI studies of superconductors the SuperCon database [16] was
used, which contains a number of wrong entries, either in terms of chemical
composition or in terms of superconducting critical temperature Tc. We
previously estimated that they might account for up to 20$\%$ of the entire
database, in particular for the cuprate entries [6]. These wrong entries can
have strong effect on AI calculations, as they can create false clusters.
Future calculations would benefit from an improved database, as we show here
by using our own data sets compiled from literature.

###  2.5 Dimensionality Reduction and Data Visualization

Similar to what we did in our previous publication [6], we used element-
vectors to mathematically represent superconductors based on their chemical
composition. In this approach, each superconductor is represented with a
1$\times$96 vector. As is sometimes done in machine learning, this dimension
can be reduced to allow more efficient calculations. A variety of different
techniques have been developed over the years: PCA (or SVD), random
projections, MDS, etc. [17]. However, our calculation with superconductors
indicate that the best results are achieved when using a nonlinear
dimensionality reduction technique called t-distributed Stochastic Neighbor
Embedding (t-SNE) [18].

t-SNE has been used extensively for machine learning calculations in materials
science, and our calculations also indicate that it is the best choice for
superconductors. It is a non-linear transformation that reduces the number of
dimensions down to two or three, none of which have any physical meaning. This
not only allows more efficient calculations, but also visualization of high
dimensional data in a low dimensional space. Another useful feature of t-SNE
is that it allows easy manual backward validation, which we used extensively
to investigate interesting data points, as well as outliers.

One must be careful when using t-SNE, however. Due to its stochastic nature,
its results can be difficult to interpret [19]. The most important parameter
for t-SNE calculations is the so-called perplexity. Perplexity value must be
carefully chosen to avoid misinterpretations [19]. In our calculations we
performed a number of runs with different values of perplexity in order to
eliminate spurious features. Depending on the size of analyzed data set, the
values of perplexity we used were between 5 and 500. Fig.1 shows schematically
the work flow chart for the whole project.

Figure 1: (Color online). The work flow chart for the calculations done in
this project.

##  3 Clustering Test Data Sets

In order to test the feasibility of clustering methods, we first complied
several test sets of superconductors from literature [9, 10], each consisting
of 25–40 materials. Such a small number of entries allows us to look into each
superconductor individually and verify its cluster assignment. In the first
set, we included a variety of different superconductors from literature,
including elements, cuprates, pnictides, heavy fermions, transition metal
dichalcogenides, Chevrel compounds, etc. In Fig. 2 we display the results of
clustering analysis and plot them using t-SNE. The two axes (Y1 and Y2) are
the reduced dimensions of the original data set; they do not have any physical
meaning. One can readily identify three major clusters, which correspond to
cuprates (red circles), pnictides (green circles) and what we previously [6]
called ”others” (blue circles). The three clusters are well separated from
each other, and no data points can be found between them. It is also possible
to identify some sub-cluster structures in Fig. 2. However in our experience,
it is better to do clustering in steps: once we identify the three main
clusters, we cluster each of them separately.

It is interesting to notice that in Fig. 2 oxide superconductors Sr2RuO4,
B0.6K0.4BiO3, LiTi2O4 and the recently discovered [20] nickelate
Nd0.8Sr0.2NiO2 are clustered with the cuprates. This should not be surprising,
as they are chemically much more similar to cuprates (copper-oxides) than to
any other group of superconductors shown in Fig. 2.

Figure 2: (Color online). t-SNE plot of the test set consisting of 37
superconductors complied from literature. The two axes do not have any
physical meaning. Three clusters can be identified and they correspond to
cuprates (red circles), pnictides (green circles) and other (blue circles)
superconductors.

The second test set includes only cuprates [9, 10, 21] and has 39 materials.
Fig. 3 shows the results of clustering, where one can identify eight major
cuprate families: BSCCO (red), LSCO (black), YBCO (purple), NCCO (orange), Hg-
based (dark yellow), Tl-based (blue), Pb-based (green) and infinite layer
(magenta). We notice that the hole-doped LSCO family is separated from the
electron-doped NCCO family, even though they are chemically and
crystallographically very similar and are often assumed to belong to the same
family [9]. Similarly, Tl-based and Hg-based cuprates are considered to be in
the same family, but the algorithm was able to separate. Not surprisingly,
they are located close to each other in Fig. 3.

Figure 3: (Color online). t-SNE plot of the cuprate test set consisting of 39
cuprate superconductors. Machine learning clustering identified eight major
cuprate families: BSCCO (red), LSCO (black), YBCO (purple), NCCO (orange), Hg-
based (dark yellow), Tl-based (blue), Pb-based (green) and infinite layer
(magenta).

The last test set contains 25 pnictide superconductors [22, 23]. The results
of clustering are shown in Fig. 4. One can identify five major pnictide
families: 11/111 (green), 122 (magenta), 1111 (red), 245 (blue), and 42622
(black). All five families are well separated, with no data points in between
clusters. We notice that there might be some sub-cluster structure within the
11/111 family, which might be resolvable with a larger data set.

Figure 4: (Color online). t-SNE plot of the pnictide test set consisting of 25
pnictide superconductors. Machine Learning clustering was able to identify
five major families: 11/111 (green), 122 (magenta), 1111 (red), 245 (blue),
and 42622 (black).

##  4 Clustering SuperCon Database

The results presented in the previous section (Figs. 2, 3 and 4) indicated
that clustering methods were capable of distinguishing between different types
of superconductors. We then proceeded to cluster the SuperCon database [16].
After removing incorrect, incomplete and multiple entries, we were left with
about 16,000 superconducting materials. In order to improve the learning
process, we removed entries with eight or more elements. We also note that the
number of pnictide entries in the database is much smaller compared with
cuprates and others. Our calculations indicate that clustering works better
when the number of constituents in each cluster is similar. Therefore, for our
calculations we randomly selected between 1,400 and 1,500 superconductors of
each type, for a total of approximately 4,500 entries.

In Fig. 5 we show the results of clustering applied to this data set. Similar
to Fig. 2, one can identify three major clusters, which correspond to
cuprates, pnictides and others. However, unlike clustering achieved for the
test set in Fig. 2, the separation of clusters is not as clear, and there are
points between the clusters. This is not surprising considering the size and
quality of SuperCon database.

Figure 5: (Color online). t-SNE plot of 4,500 randomly selected
superconductors from the SuperCon database. Three clusters are clearly seen,
and they correspond to cuprate (red), pnictide (green) and other (blue)
superconductors.

After that, we extracted cuprates from the SuperCon database as
superconductors containing both copper and oxygen. The results of clustering
are shown in Fig. 6. We can identify seven clusters, which correspond to YBCO
(red), LSCO (dark yellow), BSCCO (orange), Hg-based (blue), Tl-based
(magenta), NCCO (green) and a cluster (black) that contains a variety of
different cuprates, which the program was unable to cluster with any other
family. Manual backward validation indicates that this cluster contains mostly
cuprates with rare-earth elements like lutetium, dysprosium, erbium, holmium,
etc. We also notice that, unlike our test set in Fig. 3, we are not able to
identify the Pb-based and infinite layer families in Fig. 6, which we
speculate is because of only a small number of their entries existing in the
database.

Figure 6: (Color online). t-SNE plot of cuprates. Seven clusters are clearly
identified, and they correspond to the following families: YBCO (red), LSCO
(dark yellow), BSCCO (orange), Hg-based (blue), Tl-based (magenta), NCCO
(green) and a cluster (black) that contains a variety of different cuprates,
which the program was unable to cluster with any other family.

Next, we extracted pnictides from the SuperCon database as entries that have
iron and at least one of the following: arsenic, phosphorous, selenium or
tellurium. This resulted in a set with only about 1,400 entries. In Fig. 7 we
show the result of clustering. The t-SNE plot indicates that there are four
well separated clusters: 11/111 (red), 122 (black), 1111 (green) and 245
(blue). We notice that the 11/111 and 245 clusters are much more compact
compared to 122 and 1111 clusters. This might be due to the fact that there
are significantly fewer entries of 11/111 and 245 superconductors in the
SuperCon database. It is also clear that there is some sub-cluster structure
in the 122 and 1111 clusters, which might be resolvable with a larger data
set. ML clustering was not able to resolve the 42622 family, as there were
very few entries of them in the database.

Figure 7: (Color online). t-SNE plot of 1,400 pnictides from the SuperCon
database. Four clusters are clearly identified, and they correspond to the
following families: 11/111 (red), 122 (black), 1111 (green) and 245 (blue).

After removing cuprates and pnictes, we were left with approximately 7,000
superconductors, which we refer to as others. This is a very diverse set of
materials, and clustering it is not an easy task. As mentioned above, for most
clustering algorithms the number of desired clusters is an input parameter.
However, for other superconductors we do not a priori know the number of
clusters. As a result, we used DBSCAN algorithm which does not require the
number of clusters to be pre-specified [24].

The results of calculations on other superconductors are shown in Fig. 8. The
algorithm separated this data set into 24 clusters. We were able to identify
them and they are listed in Table 1. For each cluster several representative
examples are shown. We notice that in each cluster there is predominant
element or a group of elements, which is the common thread for that particular
cluster. This is not surprising, as the chemical composition was used as the
only predictor in these calculations.

DBSCAN also identify a number of outliers among other superconductors, shown
with open circles in Fig. 8. These points are usually separated from the main
clusters. Some outliers are: BaHg, LiHg3, Ni0.3Th0.7, Th0.9988U0.0012,
Cr0.75Ru0.25, Co0.96Mo0.04U6, Co0.94Rh0.06U6, Gd0.001Th0.999, Ag0.7Zn0.3,
Ag0.625Al0.375, etc. In our opinion, these materials, and in particular their
derivatives, should be explored further as seeds for potentially new classes
of superconducting materials.

Figure 8: (Color online). t-SNE plot of other supercondutors from the SuperCon database. DBSCAN was used for clustering, as it does not require the number of clusters to be specified. The program separated the data set into 24 clusters. We were able to identify them, and they are listed in Table. 1. DBSCAN was also able to identify several outliers, which are shown with open circles.  Table 1: Clusters identified in Fig. 8 using DBSCAN. Several examples of superconductors from each cluster are shown in the last column. Cluster number | Common thread | Examples
---|---|---
1 | Zr- and Ti-based | Zr2Co0.9Cu0.1, Zr0.8Mo0.2, Ir0.1Ti0.9 , Ti0.576V0.384Ru0.04
2 | V-based | PtV3, Ga0.3V0.7, V3Si0.6Sn0.4
3 | Si-based | SrGa0.5Si1.5, CaAl0.8Ga0.2Si, Lu2Ir3Si5
4 | S-based | Cu1.5Mo4.5S6, NbPbS3, Ni0.02TaS2
5 | In-based | CeIn3, ThIn1.5Sn1.5, Bi0.343In0.657
6 | Pd-based | Li2Pd3B, As0.5Ni0.06Pd0.44, ZrPd2Al
7 | Te-based | IrTe3, Pb0.985Tl0.015Te, Ir0.95Pt0.05Te2
8 | La-based | Bi3La4, La3Rh2Ge2, AlGd0.003La2
9 | Se-based | TiSe2, Sr0.2WSe2, La2.85Pr0.15Se4
10 | C-based | C1.35La, Cs2RbC60, La0.2Th0.8NiC2
11 | Pt-based | SrPt3P, Pt5Th, PrPt2B2C
12 | Ni-based | GdSmNi2, HoNi2B2C, MgC1.25Ni3
13 | Nb-based | Nb3Sb, AlGeNb3, Nb4FeSi
14 | Ge-based | ThPt2Ge2, Pr0.5Eu0.5Pt4Ge12, CeNiGe3
15 | B-based | MgB2, DyRh2Ir2B4, Ca3Rh8B2
16 | Be-based | LuBe13, Be0.9Ni0.1, MgReBe12
17 | Oxides | Ba0.19K0.81BiO3, Ag5Pb2O6, Bi4S2.91Se0.09O4
18 | Bi-based | Bi0.6Tl0.4, Ba2Bi3, Bi0.525Pb0.32Sn0.155
19 | Au-based | Au0.84In0.16, AlScAu2, Au0.85Pd0.15Ga2
20 | Ir-based | CaIr2, SrIr2As2, Ir0.8Pt0.2
21 | Sn-based | NbSn3, Ca2.25Sr0.75Rh4Sn13, RhSn2
22 | Ta-based | Hf0.4Ta0.6, Be2Ta3, Nd1.72Ta3.28S2
23 | Ru-based | Ce0.85La0.15Ru2, Ir0.43Ru0.57, Cr0.5Ru0.5
24 | NbN-based | NbC0.5N0.47, Nb2BN, NbN0.87O0.13

##  5 Summary

We have made the first attempt at clustering superconductors using ML methods.
It was shown that ML clustering can achieve, and in some cases (cuprates, for
example) exceeded human level performance. Clustering was done based on
chemical composition only. We speculate that when additional physical
parameters are included as predictors, even better results might be obtained.
In particular, one should consider empirical predictors that are believed [28]
to have strong relations with superconductivity, such as crystal symmetry,
layered structure, density of states at the Fermi level, the number of valence
electrons, etc. We showed that the best results are achieved when
superconducting materials are clustered in steps. A similar approach might be
useful in other data-driven materials calculations.

We also discussed the limiting factors, i.e. wrong entries in the database.
One can hope that with a more reliable database, which will also include more
predictors, better clustering results can be achieved. In that regard, the
recently initiated SuperMat project [29], which uses text mining and other
natural language processing techniques to construct a superconductor database,
might turn out to be fruitful.

##  6 References

## References

  * [1] V. Stanev, C. Oses, A.G. Kusne, E. Rodriguez, J. Paglione, S. Curtarolo and I. Takeuchi, NPJ Computational Materials, 4, 29 (2018).
  * [2] K. Hamidieh, Computational Materials Science, 154, 346 (2018).
  * [3] T. Konno, H. Kurokawa, F. Nabeshima, Y. Sakishita, R. Ogawa, I. Hosako, A. Maeda, Physical Review B, 103, 014509 (2021).
  * [4] K. Matsumoto and T. Horide, Applied Physics Express, 12, 073003 (2019).
  * [5] S. Zeng, Y. Zhao, G Li, R. Wang, X. Wang and J. Ni, NPJ Computational Materials, 5, 84 (2019).
  * [6] B. Roter and S.V. Dordevic, Physica C, 575, 1353689 (2020).
  * [7] Z.-L. Liu, P. Kang, Y. Zhu, L. Liu and H. Guo, APL Materials, 8, 061104 (2020).
  * [8] B. Meredig, E. Antono, C. Church, M. Hutchinson, J. Ling, S. Paradiso, B. Blaiszik, I. Foster, B. Gibbons, J. Hattrick-Simpers, A. Mehtaf and L. Ward, Molecular Systems Design and Engineering, 3, 819 (2018).
  * [9] C.P. Poole Jr., H.A. Farach, R.J. Creswick, R. Prozorov, Superconductivity, Elsevier, 3rd edition, (2014).
  * [10] J.E. Hirsch, M.B. Maple and F. Marsiglio (editors), Superconducting Materials: Conventional, Unconventional and Undetermined, Physica C, Vol.514, (2015), and references therein.
  * [11] J. Hill, G. Mulholland, K. Persson, R. Seshadri, C. Wolverton and B. Meredig, MRS Bulletin, 41, 399 (2016).
  * [12] R. Xu and D. Wunsch, Clustering, Wiley-IEEE Press (2008).
  * [13] J. Kogan, C. Nicholas and M. Teboulle (eds.), Grouping multidimensional data: recent advances in clustering, Berlin; New York, Springer (2006).
  * [14] A.J. Parker and A.S. Barnard, Advanced Theory and Simulations, 2, 1900145 (2019).
  * [15] Y. LeCun, Y. Bengio and G. Hinton, Nature, 521, 436 (2015), and references therein.
  * [16] SuperCon database: https://supercon.nims.go.jp
  * [17] A. Matasov and V. Krasavina, SN Applied Sciences, 2, 1463 (2020).
  * [18] L.v.D. Maaten and G. Hinton, Journal of Machine Learning Research, 9, 2579 (2008).
  * [19] M. Wattenberg, F. Viegas and I. Johnson, How to Use t-SNE Effectively, Distill, (2016). http://doi.org/10.23915/distill.00002
  * [20] D. Li, K. Lee, B.Y. Wang, M. Osada, S. Crossley, H.R. Lee, Y. Cui, Y. Hikita and H.Y. Hwang, Nature, 572, 624 (2019).
  * [21] C.W.Chu, L.Z.Deng and B. Lv, Physica C, 514, 290 (2015).
  * [22] H. Hosono and K. Kuroki, Physica C, 514, 399 (2015).
  * [23] K. Ido, Y. Motoyama, K. Yoshimi and T. Misawa, arXiv:2109.09121
  * [24] DBSCAN requires two input parameters to be specified: $\epsilon$ and minPts. As has been pointed out before, specifying these two parameters is not trivial [12]. Their values determine both the number and content of clusters. In recent years, there have been numerous attempts to algorithmically determine the values of $\epsilon$ and minPts [25, 26, 27]. However, as stated above, in our experience the best results when clustering superconductors are obtained when machine learning techniques are used in concert with human knowledge of superconductors. Therefore, we have manually adjusted the values of $\epsilon$ and minPts until the most physically meaningful clusters were obtained. The values used for clustering in Fig.8 were $\epsilon$ = 3 and minPts = 30 and the distance function was cityblock.
  * [25] E. Schubert, J. Sander, M. Ester, H.-P. Kriegel and X. Xu, ACM Transactions on Database Systems, 42, 19 (2017).
  * [26] N. Rahmah and I.S. Sitanggang, IOP Conference Series: Earth and Environmental Science, 31, 012012 (2016).
  * [27] A. Starczewski, P. Goetzen and M.J. Er, Journal of Artificial Intelligence and Soft Computing Research, 10, 209 (2020).
  * [28] K. Conder, Superconductor Science and Technology, 29 080502 (2016).
  * [29] L.Foppiano, S. Dieb, A. Suzuki, P.B. Castro, S, Iwasaki, A. Uzuki, M.G.E. Echevarria, Y. Meng, K, Terashima, L, Romary, Y, Takano and M. Ishii, Science and Technology of Advanced Materials: Methods, 1, 34 (2021).


# Maximum principle for discrete time mean-field stochastic optimal control
problems

Arzu Ahmadova1 Nazim I. Mahmudov2

1Faculty of Mathematics, University of Duisburg-Essen, 45127, Essen, Germany,
e-mail<EMAIL_ADDRESS>
2Department of Mathematics, Eastern Mediterranean University, 99628, T.R.
North Cyprus,
e-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, we study the optimal control of a discrete-time stochastic
differential equation (SDE) of mean-field type, where the coefficients can
depend on both a function of the law and the state of the process. We
establish a new version of the maximum principle for discrete-time stochastic
optimal control problems. Moreover, the cost functional is also of the mean-
field type. This maximum principle differs from the classical principle since
we introduce new discrete-time backward (matrix) stochastic equations. Based
on the discrete-time backward stochastic equations where the adjoint equations
turn out to be discrete backward SDEs with mean field, we obtain necessary
first-order and sufficient optimality conditions for the stochastic discrete
optimal control problem. To verify, we apply the result to production and
consumption choice optimization problem.

Keywords: Discrete time stochastic maximum principle, backward stochastic
difference equations, mean-field theory, optimal control problem, necessary
and sufficient conditions

###### Contents

  1. 1 Introduction
  2. 2 Mathematical description
  3. 3 Statement of main results
  4. 4 Backward Stochastic Difference Equations
  5. 5 Proof of Theorem 2
    1. 5.1 Duality analysis
  6. 6 Sufficient conditions for optimality
  7. 7 Application
  8. 8 Conclusions

##  1 Introduction

A large number of problems, interesting from a theoretical point of view and
important from a practical one, has attracted the attention of many
mathematicians and engineers. It is not surprising that there is no field in
which extremal problems do not arise, and in which it is not essential to the
development of these fields that such problems should be solved. The
development of the necessary conditions for an extremum was the elaboration of
convex programming theory. A central place in this theory is occupied by the
Kuhn-Tucker theorem. The embedding of the theory of optimal control in a
general theory of necessary conditions was first carried out by Milyutin and
Dubovitskii [1]. The great importance of their work lies in the fact that they
succeeded in formulating in a refined form necessary conditions for an
extremum which can be applied to a wide class of problems.

The maximum principle for discrete-time systems has become a subject of great
interest. We begin this section by summarizing some seminal articles in this
field. In [3] it was shown that the convexity requirement is not applicable to
many practical systems. Holtzman and Halkin [4] extend the applicability to
much broader classes of practical systems under the condition of directional
convexity being weaker than convexity. Moreover, Gamkrelidze [12] proved a
maximum principle for systems with phase constraints under a number of
assumptions. A number of original ideas related to proving the maximum
principle can be found in the works of Rozenoer [6].

Jordon and Polak [7] have also considered the problem for optimal discrete
systems and derived a stationary principle. They applied similar arguments to
those used in deriving the Pontryagin maximum principle for continuous-time
systems [5]. Butkovski [8] first showed that, in contrast to the continuous
case, a direct extension of Pontryagin’s maximum principle to discrete systems
is in general impossible. Of course, such a property of these systems is of
theoretical interest to researchers. He clearly demonstrated some errors in
the existing works. The intrinsic reason for the errors is that the
significance of convexity has been ignored. Generally speaking, the discrete-
time maximum principle fails unless a certain convexity precondition is
imposed on the control system. However, in this connection, some researchers
have established additional conditions, such as convexity of the set of
admissible velocities of the system, directional convexity and $z$-directional
convexity, etc., and found that under these conditions the maximum principle
is valid for discrete control systems.

Many results have been done on this topic for different kinds of continuous-
time stochastic optimal control problems, for example [10, 14, 17, 25, 9, 27,
28, 29], and discrete-time stochastic optimal control problems, see [16, 11,
muller, 18, 19, 20, 21, 22, 13, 26] and the references therein). The main
difficulty of the stochastic maximum principle for an optimal control problem
governed by continuous-time stochastic Itô equations is that the stochastic
Itô integral is only of order $\varepsilon$ (”hidden convexity” fails).
Therefore, the usual method of first-order needle variation fails. To overcome
this difficulty, one has to study both the first and second order terms in the
Taylor expansion of the needle variation, and establish a stochastic maximum
principle consisting of two backward stochastic differential equations and a
maximum condition with an additional quadratic term in the diffusion
coefficients, see [14, 15]. It should be noted that Lin and Zhang [21] used
spike variations to show that the necessary condition for discrete-time
stochastic optimal problems is associated with the solutions of a pair of
discrete-time backward stochastic equations. On this basis, they obtained the
maximum principle for the discrete-time stochastic optimal control problem.

As for the discrete maximum principle for mean-field stochastic optimal
control problems framework, there are a few papers dealing with discrete-time
mean-field stochastic optimal control. Unlike the classical stochastic control
problem, mean-field terms appear in the system dynamics and cost function,
connecting mean-field theory to stochastic control problems. The stochastic
mean-field control problem has been an important research topic since the
1950s. The system state is described by a controlled mean-field stochastic
differential equation (MF-SDE), which was first proposed in [5], and the first
study on MF-SDEs was published in [2]. Since then, many researchers have made
numerous contributions to the study of MF-SDEs and related topics, see, e.g.
[3, 9, 11, 15, 23, 24] and the references cited therein.

Among the many scientific articles on discrete stochastic maximum principle,
we will mention only a few with comparison and relation that motivate this
work:

  * •

Song and Liu considered in [10] the optimal control problem for fully coupled
forward–backward stochastic difference equations of mean-field type under weak
convexity assumption. Note that the form of (3.6) as an adjoint equation which
was introduced in [10] is one kind of backward stochastic difference equation.
This adjoint equation is quite different from our adjoint equation (3) studied
in this paper. One the one hand, they have different forms, on the other hand,
the adjoint equation (3.6) is $\mathcal{F}_{t+1}$-measurable;

  * •

In [26], Wu and Zhang studied recently discrete-time stochastic optimal
control problem with convex control domains, for which necessary condition in
the form of Pontryagin’s maximum principle and sufficient condition of
optimality are derived. They also pointed out that how to overcome of
integrability problem of the solution to the adjoint equation which was not
taken into consideration in [10].

  * •

Recently, Mahmudov [13] derived the first-order and second-order necessary
optimality conditions for discrete-time stochastic optimal control problems by
virtue of new discrete-time backward stochastic equation and backward
stochastic matrix equation under assumption of the set
$(f,\sigma_{1},\sigma_{2},\ldots,\sigma_{d},l)(t,\bar{x}(t),U(t))$ being
convex. Unlike [13] and [26] we study mean-field type discrete-time stochastic
maximum principle.

Based on the above considerations, the main purpose of this paper is to
construct a rigorous mathematical framework for a mean-field type of discrete-
time stochastic optimal control problems and to obtain a rigorous maximum
principle in an understandable way. We study the maximum principle for the
optimal control of discrete-time systems described by mean-field stochastic
difference equations. As far as we know, there are few results on such
stochastic control problems. In fact, discrete-time control systems are of
great value in practice. For example, digital control can be formulated as a
discrete-time control problem in which the sampled data are obtained at
discrete times. In a discrete-time system, the Riccati difference equation
plays an important role in synthesizing the optimal control. As pointed out in
[26], the integrability of the solution to the adjoint equation in discrete-
time stochastic optimal control problem is completely different from that in
the continuous-time case. However, we also prove that the solution of the
adjoint equation has no problem of integrability.

The main perspectives of our work are systematized as below:

  * •

First, to study discrete stochastic optimal control problems, we use the
finite approximation method applied in [20]. We extend this method to study
the discrete-time stochastic backward equation and introduce the discrete-time
stochastic backward matrix equation;

  * •

Second, we prove that the solution to the adjoint equation has no problem of
integrability;

  * •

Next, a constructive method is that when the necessary optimality condition
are also sufficient under certain assumptions;

  * •

Finally, as an application, we adapt the practical application based on
Theorem 2 and consider the discrete-time system with some risk in the
investment process.

The structure of the paper is as follows. In Section 2, we formulate the main
results and give an example to show the applicability of our results. Section
3 is devoted to stating main results of this paper. In Section 4, we introduce
the discrete-time backward stochastic equation and the discrete-time backward
stochastic matrix equation and present the solutions in terms of a fundamental
stochastic matrix. In Section 5, we prove the discrete-time stochastic maximum
principle: a first-order necessary condition for optimality. Section 6 is
devoted to the sufficient condition for optimality. Section 7 is devoted to
the application of production and consumption choice optimization problems.

##  2 Mathematical description

In Section 2 we present in Setting 1 the mathematical framework which we use
to study the discrete-time stochastic optimal control problems of mean-field
type.

######  Setting 1.

Let $\left\|\cdot\right\|$ be a norm, $\langle\cdot,\cdot\rangle$ be an inner
product, let $n_{1},n_{2}\in\mathbb{N}$ and denote the space of ($n_{1}\times
n_{2}$)-matrices by $\mathbb{R}^{n_{1}\times n_{2}}$, and let
$\mathbb{R}^{n_{1}}\coloneqq\mathbb{R}^{n_{1}\times 1}$, that is, each element
of $\mathbb{R}^{n_{1}}$ is understood as a column vector, let $I$ be the unit
matrix with appropriate dimension. For each matrix $A$, $A^{\intercal}$
denotes the transpose of $A$. Moreover, the forward difference operator
$\Delta$ is defined for all $h>0$ as $\Delta f(t)=f(t+h)-f(t)$. For a vector
$x\in\mathbb{R}^{n}$ denote by $x^{\intercal}$ its transpose. For a symmetric
matrix $A$ and vectors $y,y_{1},y_{2}$ of matching dimensions, we denote
$A[y]^{2}\coloneqq y^{\intercal}Ay$, $A[y_{1},y_{2}]\coloneqq
y_{1}^{\intercal}Ay_{2}$, $\widehat{f}[t]\coloneqq
f(t,\widehat{x}(t),\textbf{E}\widehat{x}(t),\widehat{u}(t))$.

Let $(\Omega,\mathfrak{F},\mathbb{P})$ be a complete probability space and $N$
be a positive integer. $\mathbb{T}:=\left\\{t_{k}=t_{0}+kh,\
h>0\right\\}_{k=0}^{N}$, let $\left\\{w(t_{k}):k=1,\ldots,N+1\right\\}$ be a
sequence of $\mathfrak{F}_{k}$-measurable $\mathbb{R}^{d}$-valued random
variables, and let $\mathfrak{F}_{k}\subseteq\mathfrak{F}$ be the
$\sigma$-field generated by $w(t_{1}),\ldots,w(t_{k})$, i.e.,
$\mathfrak{F}_{k}=\sigma\left\\{w(t_{1}),\ldots,w(t_{k})\right\\}$,
$k=1,\ldots,N+1$, and $\mathfrak{F}_{0}=\left\\{\emptyset,\Omega\right\\}$.
Let the expectation operator E be denoted by
$\textbf{E}x(t)=\int_{\Omega}x(t)\mathbb{P}(\mathrm{d}w)$ for each
$w\in\Omega$. For each $t\geq 0$,
$\textbf{E}\left\\{\cdot\mid\mathfrak{F}_{t}\right\\}$ is the conditional
expectation given by $\mathfrak{F}_{t}$. Assume for all $k\in\mathbb{N}$ that
$w_{h}(t_{k})\coloneqq w(t_{k+1})-w(t_{k})$ satisfies the following
conditions:

  1. (wi)

For every
$w_{h}\left(t_{k}\right)=\left(w_{h}^{1}\left(t_{k}\right),...,w_{h}^{d}\left(t_{k}\right)\right),\
w_{h}^{1}\left(t_{k}\right),...,w_{h}^{d}\left(t_{k}\right)$ are independent
$\mathbb{R}$-valued random variables.

  2. (wii)

$\textbf{E}\left\\{w_{h}\left(t_{k}\right)\mid\mathfrak{F}_{k}\right\\}=0$,
$\textbf{E}\left\\{\left(w_{h}^{j}\left(t_{k}\right)\right)^{2}\mid\mathfrak{F}_{k}\right\\}=h,\\\
\textbf{E}\left(w_{h}^{j}\left(t_{k}\right)\right)^{4}<\infty,$
$\textbf{E}\left(w_{h}^{m}\left(t_{k}\right)w_{h}^{l}\left(t_{k}\right)\right)=\left(t_{k+1}-t_{k}\right)\delta_{ml}I$.

Moreover, let
$\widehat{\mathfrak{F}}_{k}=\sigma\left\\{w_{h}\left(t_{k+1}\right),\ldots,w_{h}\left(t_{N}\right)\right\\}$.
Note that $\mathfrak{F}_{k}$ and $\widehat{\mathfrak{F}}_{k}$ are independent.
Let for all $v\in\mathbb{R}^{r}$ $\Delta f(t,v)\coloneqq
f(t,\widehat{x}(t),\textbf{E}\widehat{x}(t),v)-f(t,\widehat{x}(t),\textbf{E}\widehat{x}(t),\widehat{u}(t))$,
and let $\mathbb{F}=\left\\{\mathfrak{F}_{k}:k=0,1,...,N\right\\}$ be the set.
A random variable $z=\left\\{z_{k}:k=0,1,...,N\right\\}$ is called
$\mathbb{F}$-predictable if the random variable $z_{k}$ is
$\mathfrak{F}_{k}$-measurable for every $k=0,1,...,N,$. Let
$L^{2}(\Omega,\mathfrak{F}_{t_{k}},\mathbb{R}^{n})$ be the set of all
$\mathbb{R}^{n}$-valued $\mathfrak{F}_{t_{k}}$-measurable random variables
$x(t_{k})$ with $\textbf{E}\|x(t_{k})\|^{2}<\infty$.

##  3 Statement of main results

In Section 3, we establish a class of discrete-time stochastic nonlinear
optimal control problems of mean-field type. The system equation is the
following nonlinear stochastic difference equation:

$\begin{cases}x\left(t+h\right)=x(t)+hf\left(t,x(t),\mathbf{E}x(t),u(t)\right)+\\\
\sum\limits_{j=1}^{d}\sigma^{j}\left(t,x(t),\mathbf{E}x(t),u(t)\right)w_{h}^{j}\left(t\right),\\\
\textbf{E}x\left(t+h\right)=\textbf{E}x\left(t\right)+h\textbf{E}f\left(t,x\left(t\right),\mathbf{E}x\left(t\right),u\left(t\right)\right),\\\
x\left(t_{0}\right)=\xi\in\mathbb{R}^{n},\
u\left(t\right)\in\mathcal{U}(t)\subset\mathbb{R}^{r},\quad
t\in\mathbb{T}.\end{cases}$ (1)

Note that the initial value $\xi$ and $\left\\{w(t_{k}),k\in 0,1,\ldots
N\right\\}$ are assumed to be independent of each other. Let
$\varphi:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R},\
\left(l,f,\sigma^{j}\right):\mathbb{T}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{r}\rightarrow\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n}.$
$\left\\{w\left(t_{k}\right):k=1,\dots,N+1\right\\}$ be a sequence of
$\mathfrak{F}_{k}$-measurable $\mathbb{R}^{d}$-valued random variables. Then
the optimal control problem minimizes the following expected cost functional
defined by:

$\mathfrak{J}(u)=\textbf{E}\varphi\left(x\left(t_{N+1}\right),\mathbf{E}x\left(t_{N+1}\right)\right)+\textbf{E}\sum\limits_{t=t_{0}}^{t_{N}}l\left(t,x\left(t\right),\mathbf{E}x\left(t\right),u\left(t\right)\right)\longrightarrow\min.\\\
$ (2)

The cost functional (2) is also of mean-field type since both the running and
terminal cost functions $l$ and $\varphi$ depend on the state process through
their expected values.
Observe from (1) and (2) that
$\mathbf{x}=\left\\{x\left(t_{k}\right)\right\\}_{k=0}^{N+1}$ and
$\mathbf{u}=\left\\{u\left(t_{k}\right)\right\\}_{k=0}^{N}$ are the state
process and control process, respectively. Let
$\left\\{\mathcal{U}(t_{k})\right\\}_{k=0}^{N}$ be a sequence of nonempty
convex subset of $\mathbb{R}^{r}$. We introduce the following admissible
control set

$\mathfrak{U}_{ad}=\left\\{u=\left\\{u\left(t_{k}\right)\right\\}_{k=0}^{N}:u\left(t_{k}\right)\in
L^{2}\left(\Omega,\mathfrak{F}_{t_{k}},\mathbb{R}^{r}\right)\quad\text{and}\quad
u(t_{k})\in\mathcal{U}(t_{k})\right\\}.$

The pair $\left(\mathbf{x},\mathbf{u}\right)$ satisfying the constraints (1)
is called an admissible pair, and the pair
$\left(\widehat{\mathbf{x}},\widehat{\mathbf{u}}\right)$ which is a solution
of the problem (1)-(2), is called an optimal pair. Our optimal control problem
can be stated as follows:

Problem ($\mathcal{DOPC}$). Minimize (2) over $\mathcal{U}(t)$. Any
$\widehat{u}(\cdot)\in\mathcal{U}(t)$ satisfying

$\mathfrak{J}(\widehat{u}(\cdot))=\inf_{v\in\mathcal{U}(t)}\mathfrak{J}(v(\cdot))$

is called an optimal control. The corresponding
$\widehat{x}(\cdot)=\widehat{x}(\cdot,\widehat{u})$ and
$(\widehat{x}(\cdot),\widehat{u}(\cdot))$ are called optimal state process and
optimal pair, respectively.

Throughout the paper we use the following assumptions.

(A1)

Let $\psi\coloneqq l,f,\sigma^{j},\varphi$. There exists a constant $L>0$ such
that

$\displaystyle\left\|\psi\left(t,x_{1},y_{1},u_{1}\right)-\psi\left(t,x_{2},y_{2},u_{2}\right)\right\|$
$\displaystyle\leq
L\left(\left\|x_{1}-x_{2}\right\|+\left\|y_{1}-y_{2}\right\|+\left\|u_{1}-u_{2}\right\|\right),\
$ $\displaystyle\left\|\psi\left(t,0,0,0\right)\right\|$ $\displaystyle\leq
L,\quad t\in\mathbb{T},\quad x_{1},x_{2},y_{1},y_{2}\in\mathbb{R}^{n},\
u_{1},u_{2}\in\mathbb{R}^{r}.$

(A2)

Let $\psi\coloneqq l,f,\sigma^{j},\varphi$ be continuously differentiable with
respect to $x$. Moreover, there exists a constant $L_{1}>0$ such that

$\displaystyle\left\|\psi_{x}\left(t,x_{1},y_{1},u_{1}\right)-\psi_{x}\left(t,x_{2},y_{2},u_{2}\right)\right\|$
$\displaystyle\leq
L_{1}\left(\left\|x_{1}-x_{2}\right\|+\left\|y_{1}-y_{2}\right\|+\left\|u_{1}-u_{2}\right\|\right),\
$ $\displaystyle\left\|\psi_{x}\left(t,0,0,0\right)\right\|$
$\displaystyle\leq L_{1},\quad t\in\mathbb{T},\quad
x_{1},x_{2},y_{1},y_{2}\in\mathbb{R}^{n},\ u_{1},u_{2}\in\mathbb{R}^{r};$

(A3)

The set $\left\\{\mathcal{U}(t):t\in\mathbb{T}\right\\}$ is convex;

(A4)

The set
$(f,\sigma_{1},\sigma_{2},\ldots,\sigma_{d},l)(t,\bar{x}(t),\mathcal{U}(t))$
is convex.

Under the assumptions (A1) and (A2), $J$ is well-defined on $\mathcal{U}$.
Now we state the first main result of the paper: first-order necessary
conditions for the problem (1)-(2).

######  Theorem 2.

Assume Setting 1, let $\left(\widehat{\mathbf{x}},\widehat{\mathbf{u}}\right)$
be an optimal pair in problem (1), and assume that assumptions (A1)-(A4) hold.
Moreover, assume that the function $\varphi$ satisfy the Lipschitz condition
in a neighborhood of the point $\widehat{x}\left(t_{N+1}\right)$ and
differentiable at that point. Then there exists a solution
$\left(p,q^{1},...,q^{d}\right):\mathbb{T}\times\Omega\rightarrow\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
d}$ of the discrete-time backward stochastic equation

$\displaystyle\begin{cases}p\left(t\right)=\left(I+h\widehat{f}_{x}^{\intercal}\left[t\right]\right)\mathbf{E}\left\\{p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\}+\mathbf{E}\left\\{\widehat{f}_{y}^{\intercal}\left[t\right]p\left(t+h\right)\right\\}\\\
+\sum_{j=1}^{d}\left(\widehat{\sigma}_{x}^{j}\left[t\right]\right)^{\intercal}q^{j}(t)+\sum_{j=1}^{d}\mathbf{E}\left\\{\left(\widehat{\sigma}_{y}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)\right\\}-\widehat{l}_{x}[t]-\textbf{E}\widehat{l}_{y}[t],\\\
p\left(t_{N+1}\right)=-\varphi_{x}\left(x\left(t_{N+1}\right),\mathbf{E}x\left(t_{N+1}\right)\right)-\textbf{E}\varphi_{y}\left(x\left(t_{N+1}\right),\mathbf{E}x\left(t_{N+1}\right)\right),\\\
q^{j}\left(t\right)=\mathbf{E}\left\\{p\left(t+h\right)w_{h}^{j}\left(t\right)\mid\mathfrak{F}_{t}\right\\},\end{cases}$
(3)

Note that for any $t\in\mathbb{T}$ and for any $v\in\mathcal{U}\left(t\right)$
the following Hamiltonian function $H$ as follows:

$\left\langle
H_{u}\left(t,\widehat{u}\left(t\right)\right),v-\widehat{u}\left(t\right)\right\rangle\leq
0,\ a.s.,$ (4)

where

$\displaystyle H\left(t,v\right)$ $\displaystyle\coloneqq
H\left(t,p\left(t+h\right),q\left(t\right),\widehat{x}\left(t\right),v\right)$
$\displaystyle=\left\langle\mathbf{E}\left\\{p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\},hf\left(t,\widehat{x}\left(t\right),\mathbf{E}\widehat{x}\left(t\right),v\right)\right\rangle+\sum_{J=1}^{d}\left\langle
q^{j}\left(t\right),\sigma^{j}\left(t,\widehat{x}\left(t\right),\mathbf{E}\widehat{x}\left(t\right),v\right)\right\rangle$
$\displaystyle-l\left(t,\widehat{x}\left(t\right),\mathbf{E}\widehat{x}\left(t\right),v\right).$

######  Remark 3.

It is well-known that the adjoint equation in the continuous-time case admits
a square-integrable solution under the classical assumptions. However, the
solution $\left\\{(p(t_{k}),q(t_{k}))\right\\},k\in\mathbb{T}$ of (3) has a
problem of integrability. Nevertheless, we overcome this problem based on the
following discussion. Since $w(t_{k})$ are square-integrable, then
$x(t_{k}),k\in\mathbb{T}$ are square integrable.
By Assumption $(A2)$ and (3), we obtain that

$\displaystyle\textbf{E}\left|p\left(t_{N+1}\right)\right|^{2}<\infty.$ (5)

The fact that property (wii) and (5) imply that

$\displaystyle\textbf{E}\left|q^{j}\left(t_{N}\right)\right|^{2}<\infty.$ (6)

Substituting (5) and (6) into (3) and assumption $(A2)$ ensure for all
positive integer $N$ that

$\displaystyle\textbf{E}\left|p(t_{N})\right|^{2}<\infty.$ (7)

Therefore, for all $k\in\mathbb{T}$, $\left\\{(p(t_{k}),q(t_{k}))\right\\}$ is
square integrable, so $\mathbf{E}|p(t_{k})||x(t_{k})|$ and
$\mathbf{E}|q(t_{k})||x(t_{k})|$ are well-defined, which implies that the
sufficient condition (Theorem 25 works). On the other hand,
$\left\\{u(t_{k}),k\in\mathbb{T}\right\\}$ are square-integrable, then in
general we can still get the square-integrability of
$\left\\{x(t_{k}),k\in\mathbb{T}\right\\}$, and
$\left\\{(p(t_{k}),q(t_{k}))\right\\},k\in\mathbb{T}$ are square-integrable
which also means that the proof of the sufficient condition works. In
addition, $\mathbf{E}|p(t_{k})||u(t_{k})|$ and
$\mathbf{E}|q(t_{k})||u(t_{k})|$ are well-defined for all admissible control
$u=\left\\{u(t_{k}),k\in\mathbb{T}\right\\}$ which implies that the proof of
necessary condition (Theorem 2) works. Therefore, all the expectations
involving $\left\\{(p(t_{k}),q(t_{k}))\right\\},k\in\mathbb{T}$ in this paper
are well-defined.

##  4 Backward Stochastic Difference Equations

In this section, we first define the discrete-time mean-field type backward
stochastic equations.
Let $\phi=\left\\{\phi\left(t_{k}\right)\right\\}_{k=0}^{N},$ $\
\psi^{j}=\left\\{\psi^{j}\left(t_{k}\right)\right\\}_{k=0}^{N},$
$\phi\left(t\right)\in
L^{2}\left(\Omega,\mathfrak{F}_{t},\mathbb{R}^{n}\right),$ and
$\psi^{j}\left(t\right)\in
L^{2}\left(\Omega,\mathfrak{F}_{t},\mathbb{R}^{n}\right)$, let
$A(t),A_{1}(t)\in L^{2}\left(\Omega,\mathfrak{F}_{t},\mathbb{R}^{n\times
n}\right)$, and for all $j=1,...,d$ $B^{j}(t),B_{1}^{j}(t)\in
L^{2}\left(\Omega,\mathfrak{F}_{t},\mathbb{R}^{n\times n}\right)$ be uniformly
bounded $\mathbb{R}^{n\times n}$-valued random matrices. Assume that
$z\left(t\right)$ satisfies the following discrete time difference equation

$\displaystyle\begin{cases}z\left(t+h\right)=\left(I+A\left(t\right)\right)z\left(t\right)+A_{1}\left(t\right)\mathbf{E}z\left(t\right)+\phi\left(t\right)\\\
+\sum_{j=1}^{d}\left(B^{j}\left(t\right)z\left(t\right)+B_{1}^{j}\left(t\right)\mathbf{E}z\left(t\right)+\psi^{j}\left(t\right)\right)w_{h}^{j}\left(t\right)-l(t),\\\
z\left(t_{0}\right)=\xi\in\mathbb{R}^{n},\ t\in\mathbb{T}.\end{cases}$ (8)

For arbitrary indices $t_{l},t_{k}\in\mathbb{T}\cup\left\\{t_{N+1}\right\\}$,
we introduce the $n\times n$ matrix $\Phi\left(t_{l},t_{k}\right)$,

$\displaystyle\Phi\left(t_{l},t_{k}\right)$
$\displaystyle=\left\\{\begin{tabular}[c]{lll}$0$&for&$l<k,$\\\
$I$&for&$l=k,$\\\
$\Theta\left(t_{l-1}\right)\Theta\left(t_{l-2}\right)...\Theta\left(t_{k}\right)$&for&$l>k,$\end{tabular}\
\ \ \ \ \ \right.$ (12) $\displaystyle\Theta\left(t_{k}\right)$
$\displaystyle=\left(I+A\left(t\right)\right)+A_{1}\left(t\right)\mathbf{E}\left\\{\cdot\right\\}+\sum_{j=1}^{d}\left(B^{j}\left(t\right)+B_{1}^{j}\left(t\right)\mathbf{E}\left\\{\cdot\right\\}\right)w_{h}^{j}\left(t\right).$
(13)

Observe from (12) and (13) that

$\Phi\left(t_{l},t_{0}\right)=\Phi\left(t_{l},t_{k}\right)\Phi\left(t_{k},t_{0}\right),\\\
l\geq k\geq 0.$ (14)

Then the solution of difference equation (8) can be written as

$z\left(t\right)=\Phi\left(t,t_{0}\right)\xi+\sum_{\tau=t_{0}}^{t-h}\Phi\left(t,\tau+h\right)\left(\phi\left(\tau\right)+\sum_{j=1}^{d}\psi^{j}\left(\tau\right)w_{h}^{j}\left(\tau\right)\right),\
t\in\mathbb{T}.$ (15)

Let $h_{N+1}\in
L^{2}\left(\Omega,\mathfrak{F}_{t_{N+1}},\mathbb{R}^{n}\right)$ and
$l\left(t\right)\in L^{2}\left(\Omega,\mathfrak{F}_{t},\mathbb{R}^{n}\right)$
for all $t\in\mathbb{T}$. Construct a pair of discrete-time backward
stochastic equations
$p=\left(p\left(t_{0}\right),p\left(t_{1}\right),\ldots,p\left(t_{N+1}\right)\right)$,
$q^{j}=\left(q^{j}\left(t_{0}\right),q^{j}\left(t_{1}\right),\ldots,q^{j}\left(t_{N}\right)\right)$
corresponding to (8) as follows:

$\displaystyle\begin{cases}p\left(t\right)=\left(I+A^{\intercal}\left(t\right)\right)\mathbf{E}\left\\{p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\}+\mathbf{E}\left\\{A_{1}^{\intercal}\left(t\right)p\left(t+h\right)\right\\}\\\
+\sum\limits_{j=1}^{d}\left(\left(\sigma_{x}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)+\mathbf{E}\left\\{\left(\sigma_{y}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)\right\\}\right)-l(t),\\\
p(t_{N+1})=-h_{N+1},\quad
q^{j}(t)=\textbf{E}\left\\{p(t+h)w^{j}_{h}(t)\mid\mathfrak{F}_{t}\right\\},\quad
t\in\mathbb{T}.\end{cases}$ (16)

######  Lemma 4 ([22]).

Discrete-time backward stochastic equation (16) has a unique solution
$\left(p,q\right)$ such that $p\left(t\right)\in
L^{2}\left(\Omega,\mathfrak{F}_{t},\mathbb{R}^{n}\right),q\left(t\right)\in
L^{2}\left(\Omega,\mathfrak{F}_{t},\mathbb{R}^{n\times d}\right)$, and has the
following representation

$\displaystyle p\left(t\right)$
$\displaystyle=-\mathbf{E}\left\\{\Phi^{\intercal}\left(t_{N+1},t\right)h_{N+1}+\sum_{s=t}^{t_{N}}\Phi^{\intercal}\left(s,t\right)l\left(s\right)\mid\mathfrak{F}_{t}\right\\},\
\ $ $\displaystyle\ p\left(t_{N+1}\right)$ $\displaystyle=-h_{N+1},\
q^{j}\left(t\right)=\mathbf{E}\left\\{p\left(t+h\right)w_{h}^{j}\left(t\right)\mid\mathfrak{F}_{t}\right\\},\
$ (17)

where $t\in\mathbb{T},\ j=1,\ldots,d.$

######  Remark 5.

Assume that $\Phi^{\intercal}\left(t_{k},0\right),$ $t_{k}\in\mathbb{T},$ is
invertible we have

$\displaystyle p\left(t_{k}\right)$
$\displaystyle=-\mathbf{E}\left\\{\Phi^{\intercal}\left(t_{N+1},t\right)h_{N+1}+\sum_{s=t}^{t_{N}}\Phi^{\intercal}\left(s,t\right)l\left(s\right)\mid\mathfrak{F}_{t_{k}}\right\\}$
$\displaystyle=-\left(\Phi^{\intercal}\left(t,0\right)\right)^{-1}\mathbf{E}\left\\{\Phi^{\intercal}\left(t_{N+1},0\right)h_{N+1}+\sum_{s=t}^{t_{N}}\Phi^{\intercal}\left(s,0\right)l\left(s\right)\mid\mathfrak{F}_{t_{k}}\right\\}.$

##  5 Proof of Theorem 2

Assume that
$\widehat{u}=\left\\{\widehat{u}\left(t\right)\right\\}_{t=t_{0}}^{t_{N}}$ is
the optimal control of the problem (1) and
$\widehat{x}=\left\\{\widehat{x}\left(t\right)\right\\}_{t=t_{0}}^{t_{N+1}}$
is the corresponding optimal trajectory. We fix a time $t_{0}\leq\theta\leq
t_{N}$ and choose $v\in
L^{2}\left(\Omega,\mathfrak{F}_{\theta},\mathbb{R}^{r}\right)$ such that
$u\left(\theta\right)+\Delta v\in\mathcal{U}\left(\theta\right)$. For any
$\varepsilon>0$, we define the perturbed admissible control

$u^{\varepsilon}\left(t\right)=\widehat{u}\left(t\right)+\delta_{t\theta}\varepsilon\Delta
v,\ \ t\in\mathbb{T},$

where $\delta_{t\theta}=1$ for $t=\theta,$ $\delta_{t\theta}=0$ for
$t\neq\theta$. Convexity of $\mathcal{U}\left(\theta\right)$ implies that the
control
$u^{\varepsilon}=\left\\{u^{\varepsilon}\left(t\right)\right\\}_{t=t_{0}}^{t_{N}}$
is admissible. Let $x^{\varepsilon}$ be a solution of (1) corresponding to the
control $u^{\varepsilon}$.

We introduce the following short-hand notations for $f,\ \sigma^{j}$, and $l$:

$\displaystyle\widehat{\psi}\left[t\right]$
$\displaystyle\coloneqq\psi\left(t,\widehat{x}\left(t\right),\mathbf{E}\widehat{x}\left(t\right),\widehat{u}\left(t\right)\right),\
\psi^{\varepsilon}\left[t\right]\coloneqq\psi\left(t,x^{\varepsilon}\left(t\right),\mathbf{E}x^{\varepsilon}\left(t\right),u^{\varepsilon}\left(t\right)\right);$
$\displaystyle\widehat{\psi}^{\varepsilon}\left[t\right]$
$\displaystyle\coloneqq\psi\left(t,\widehat{x}\left(t\right),\mathbf{E}\widehat{x}\left(t\right),u^{\varepsilon}\left(t\right)\right);$
$\displaystyle\widehat{\psi}_{u}\left[t;\varepsilon\right]$
$\displaystyle\coloneqq\int_{0}^{1}\psi_{u}\left(t,\widehat{x}\left(t\right),\mathbf{E}\widehat{x}\left(t\right),\widehat{u}\left(t\right)+\lambda\left(u^{\varepsilon}\left(t\right)-\widehat{u}\left(t\right)\right)\right)d\lambda;$
$\displaystyle\widehat{z}\left(t;\varepsilon\right)$
$\displaystyle:=\widehat{z}\left(t\right)+\lambda\left(z^{\varepsilon}\left(t\right)-\widehat{z}\left(t\right)\right),\quad\text{where}\quad
z=x,y,u;$ $\displaystyle\widehat{\psi}_{z}\left[t;\varepsilon\right]$
$\displaystyle:=\int_{0}^{1}\psi_{z}\left(t,\widehat{x}\left(t\right)+\lambda\left(x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)\right),\mathbf{E}\widehat{x}\left(t\right)+\lambda\left(\mathbf{E}x^{\varepsilon}\left(t\right)-\mathbf{E}\widehat{x}\left(t\right)\right),\widehat{u}\left(t\right)\right)d\lambda,$

where $z=x$ or $z=y:=\mathbf{E}x$ in the last definition.

######  Lemma 6.

Assume that assumption (A1) holds. Then we have for all $\varepsilon>0$ that

$\max_{0\leq k\leq
N}\mathbf{E}\left\|x^{\varepsilon}\left(t_{k}\right)-\widehat{x}\left(t_{k}\right)\right\|^{2}\leq
L\varepsilon^{2}\mathbf{E}\left\|\Delta v\right\|^{2}.$ (18)

###### Proof.

For $t=t_{0},...,\theta-h,$ it is clear that
$x^{\varepsilon}\left(t+h\right)=\widehat{x}\left(t+h\right),\
\mathbf{E}x^{\varepsilon}\left(t+h\right)=\mathbf{E}\widehat{x}\left(t+h\right)$.
By recursive iteration for $t=\theta,$ we have

$x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)=\widehat{f}^{\varepsilon}\left[\theta\right]-\widehat{f}\left[\theta\right]+\sum_{j=1}^{d}\left(\widehat{\sigma}^{\varepsilon
j}\left[\theta\right]-\widehat{\sigma}^{j}\left[\theta\right]\right)w_{h}^{j}\left(\theta\right).$

Then observe for all $\varepsilon>0$ that

$\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)\right\|^{2}\leq\mathbf{E}\left\|\widehat{f}^{\varepsilon}\left[\theta\right]-\widehat{f}\left[\theta\right]\right\|^{2}+\mathbf{E}\sum_{j=1}^{d}\left\|\widehat{\sigma}^{\varepsilon
j}\left[\theta\right]-\widehat{\sigma}^{j}\left[\theta\right]\right\|^{2}.$

By boundedness of $f_{u}$ and $\sigma_{u}^{j},$ and assumption (A1) we have
that

$\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)\right\|^{2}\leq
L\varepsilon^{2}\mathbf{E}\left\|\Delta v\right\|^{2}.$

For $t=\theta+h$ by boundedness of $f_{x}$ and $\sigma_{x}^{j}$ and assumption
(A1), we have that

$\displaystyle\mathbf{E}\left\|x^{\varepsilon}\left(\theta+2h\right)-\widehat{x}\left(\theta+2h\right)\right\|^{2}$
$\displaystyle\leq\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)\right\|^{2}$
$\displaystyle+\mathbf{E}\left\|f^{\varepsilon}\left(\theta+h,x^{\varepsilon}\left(\theta+h\right),\mathbf{E}x^{\varepsilon}\left(\theta+h\right),\widehat{u}\left(\theta+h\right)\right)-\widehat{f}\left[\theta+h\right]\right\|^{2}$
$\displaystyle+\mathbf{E}\sum_{j=1}^{d}\left\|\sigma^{\varepsilon
j}\left(\theta+h,x^{\varepsilon}\left(\theta+h\right),\mathbf{E}x^{\varepsilon}\left(\theta+h\right),\widehat{u}\left(\theta+h\right)\right)-\widehat{\sigma}^{j}\left[\theta+h\right]\right\|^{2}$
$\displaystyle\leq\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)\right\|^{2}+2L\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)\right\|^{2}$
$\displaystyle+2L\left\|\mathbf{E}x^{\varepsilon}\left(\theta+h\right)-\mathbf{E}\widehat{x}\left(\theta+h\right)\right\|^{2}+2L\varepsilon^{2}\mathbf{E}\left\|\Delta
v\right\|^{2}\leq L\varepsilon^{2}\mathbf{E}\left\|\Delta v\right\|^{2}.$

Similarly, by assumption (A1) we also have that

$\displaystyle\left\|\mathbf{E}x^{\varepsilon}\left(\theta+2h\right)-\mathbf{E}\widehat{x}\left(\theta+2h\right)\right\|^{2}$
$\displaystyle\leq\left\|\mathbf{E}x^{\varepsilon}\left(\theta+h\right)-\mathbf{E}\widehat{x}\left(\theta+h\right)\right\|^{2}$
$\displaystyle+\left\|\mathbf{E}f\left(\theta+h,x^{\varepsilon}\left(\theta+h\right),\mathbf{E}x^{\varepsilon}\left(\theta+h\right),\widehat{u}\left(\theta+h\right)\right)-\mathbf{E}\widehat{f}\left[\theta+h\right]\right\|^{2}$
$\displaystyle\leq\left\|\mathbf{E}x^{\varepsilon}\left(\theta+h\right)-\mathbf{E}\widehat{x}\left(\theta+h\right)\right\|^{2}+L\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)\right\|^{2}$
$\displaystyle+L\left\|\mathbf{E}x^{\varepsilon}\left(\theta+h\right)-\mathbf{E}\widehat{x}\left(\theta+h\right)\right\|^{2}\leq
L\varepsilon^{2}\mathbf{E}\left\|\Delta v\right\|^{2}.$

Therefore, by recursive iteration for $t=\theta+2h,\ldots,t_{N}$, we obtain
the desired result (18). The proof of Lemma 6 is thus complete. ∎

###  5.1 Duality analysis

Let $\xi=\left\\{\xi\left(t\right)\right\\}_{t=t_{0}}^{t_{N}}$ be the solution
to the following difference equation,

$\displaystyle\begin{cases}\xi\left(t+h\right)&=\xi\left(t\right)+\widehat{f}_{x}\left[t\right]\xi\left(t\right)+\widehat{f}_{y}\left[t\right]\textbf{E}\xi\left(t\right)+\delta_{t\theta}\widehat{f}_{u}\left[t\right]\varepsilon\Delta
v\\\
&+\sum_{j=1}^{d}\left(\widehat{\sigma}_{x}^{j}\left[t\right]\xi\left(t\right)+\widehat{\sigma}_{y}^{j}\left[t\right]\mathbf{E}\xi\left(t\right)+\delta_{t\theta}\widehat{\sigma}_{u}^{j}\left[t\right]\varepsilon\Delta
v\right)w_{h}^{j}\left(t\right),\\\ \xi\left(0\right)&=0.\end{cases}$ (19)

By Lemma 6, observe for all $\varepsilon>0$ that

$\max_{t\in\mathbb{T}}\mathbf{E}\left\|\xi\left(t\right)\right\|^{2}\leq
L\varepsilon^{2}\mathbf{E}\left\|\Delta v\right\|^{2}.$

######  Lemma 7.

Assume that assumption (A1) holds. Then we have for all $\varepsilon>0$ that

$\max_{t\in\mathbb{T}}\mathbf{E}\left\|x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)-\xi\left(t\right)\right\|^{2}=o\left(\varepsilon^{2}\right).$

###### Proof.

For $t=t_{0},...,\theta,$ observe that
$x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)-\xi\left(t\right)=0$
and
$\textbf{E}x^{\varepsilon}\left(t\right)-\mathbf{E}\widehat{x}\left(t\right)-\mathbf{E}\xi\left(t\right)=0$.
By recursive iteration, for $t=\theta$, we have

$\displaystyle
x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)-\xi\left(\theta+h\right)$
$\displaystyle=\left(\widehat{f}_{u}\left[\theta;\varepsilon\right]-\widehat{f}_{u}\left[\theta\right]\right)\varepsilon\Delta
v+\left(\widehat{f}_{v}\left[\theta;\varepsilon\right]-\widehat{f}_{v}\left[\theta\right]\right)\varepsilon\Delta\mathbf{E}v$
$\displaystyle+\sum_{j=1}^{d}\left(\widehat{\sigma}_{u}^{j}\left[\theta;\varepsilon\right]-\widehat{\sigma}_{u}^{j}\left[\theta\right]\right)\varepsilon\Delta
vw_{h}^{j}\left(\theta\right)+\sum_{j=1}^{d}\left(\widehat{\sigma}_{v}^{j}\left[\theta;\varepsilon\right]-\widehat{\sigma}_{v}^{j}\left[\theta\right]\right)\varepsilon\Delta\mathbf{E}vw_{h}^{j}\left(\theta\right).$

Then we have

$\displaystyle\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)-\xi\left(\theta+h\right)\right\|^{2}\leq
L\varepsilon^{2}\mathbf{E}\left\|\widehat{f}_{u}\left[\theta;\varepsilon\right]-\widehat{f}_{u}\left[\theta\right]\right\|^{2}\left\|\Delta
v\right\|^{2}$
$\displaystyle+L\varepsilon^{2}\mathbf{E}{\displaystyle\sum_{j=1}^{d}}\left\|\widehat{\sigma}_{u}^{j}\left[\theta;\varepsilon\right]-\widehat{\sigma}_{u}^{j}\left[\theta\right]\right\|^{2}\left\|\Delta
v\right\|^{2}.$

Similarly, we have

$\displaystyle\left\|\mathbf{E}x^{\varepsilon}\left(\theta+h\right)-\mathbf{E}\widehat{x}\left(\theta+h\right)-\mathbf{E}\xi\left(\theta+h\right)\right\|^{2}\leq
L\varepsilon^{2}\left\|\mathbf{E}\widehat{f}_{u}\left[\theta;\varepsilon\right]-\mathbf{E}\widehat{f}_{u}\left[\theta\right]\right\|^{2}\left\|\Delta
v\right\|^{2}.$

Since
$\left\|\widehat{f}_{u}\left[\theta;\varepsilon\right]-\widehat{f}_{u}\left[\theta\right]\right\|^{2},\left\|\widehat{\sigma}_{u}^{j}\left[\theta;\varepsilon\right]-\widehat{\sigma}_{u}^{j}\left[\theta\right]\right\|^{2}\rightarrow
0$ as $\varepsilon\rightarrow 0^{+},$ we get

$\displaystyle\lim_{\varepsilon\rightarrow
0}\frac{1}{\varepsilon^{2}}\mathbf{E}\left\|x^{\varepsilon}\left(\theta+h\right)-\widehat{x}\left(\theta+h\right)-\xi\left(\theta+h\right)\right\|^{2}$
$\displaystyle=0,$ $\displaystyle\lim_{\varepsilon\rightarrow
0}\frac{1}{\varepsilon^{2}}\left\|\mathbf{E}x^{\varepsilon}\left(\theta+h\right)-\mathbf{E}\widehat{x}\left(\theta+h\right)-\textbf{E}\xi\left(\theta+h\right)\right\|^{2}$
$\displaystyle=0.$

For $t=\theta+2h,...,t_{N},$

$\displaystyle
x^{\varepsilon}\left(t+h\right)-\widehat{x}\left(t+h\right)-\xi\left(t+h\right)$
$\displaystyle=\widehat{f}_{x}\left[t;\varepsilon\right]\left(x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)-\xi\left(t\right)\right)+\widehat{f}_{y}\left[\theta;\varepsilon\right]\left(\mathbf{E}x^{\varepsilon}\left(t\right)-\mathbf{E}\widehat{x}\left(t\right)-\mathbf{E}\xi\left(t\right)\right)$
$\displaystyle+\left(\widehat{f}_{x}\left[t;\varepsilon\right]-\widehat{f}_{x}\left[t\right]\right)\xi\left(t\right)+\left(\widehat{f}_{y}\left[t;\varepsilon\right]-\widehat{f}_{y}\left[t\right]\right)\mathbf{E}\xi\left(t\right)$
$\displaystyle+\sum_{j=1}^{d}\widehat{\sigma}_{x}^{j}\left[\theta;\varepsilon\right]\left(x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)-\xi\left(t\right)\right)w_{h}^{j}\left(t\right)+\sum_{j=1}^{d}\widehat{\sigma}_{y}^{j}\left[\theta;\varepsilon\right]\left(\mathbf{E}x^{\varepsilon}\left(t\right)-\mathbf{E}\widehat{x}\left(t\right)-\mathbf{E}\xi\left(t-h\right)\right)w_{h}^{j}\left(t\right)$
$\displaystyle+\sum_{j=1}^{d}\left(\widehat{\sigma}_{x}^{j}\left[t;\varepsilon\right]-\widehat{\sigma}_{x}^{j}\left[t\right]\xi\left(t\right)w_{h}^{j}\left(t\right)\right)+\sum_{j=1}^{d}\left(\widehat{\sigma}_{y}^{j}\left[t;\varepsilon\right]-\widehat{\sigma}_{y}^{j}\left[t\right]\mathbf{E}\xi\left(t\right)w_{h}^{j}\left(t\right)\right).$

Then we have

$\displaystyle\mathbf{E}\left\|x^{\varepsilon}\left(t+h\right)-\widehat{x}\left(t+h\right)-\xi\left(t+h\right)\right\|^{2}$
$\displaystyle\leq
L\mathbf{E}\left\|x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)-\xi\left(t\right)\right\|^{2}+L\left\|\mathbf{E}x^{\varepsilon}\left(t\right)-\mathbf{E}\widehat{x}\left(t\right)-\mathbf{E}\xi\left(t\right)\right\|^{2}$
$\displaystyle+\mathbf{E}\left(\left\|\widehat{f}_{x}\left[t;\varepsilon\right]-\widehat{f}_{x}\left[t\right]\right\|^{2}+\sum_{j=1}^{d}\left\|\widehat{\sigma}_{x}^{j}\left[t;\varepsilon\right]-\widehat{\sigma}_{x}^{j}\left[t\right]\right\|^{2}\right)\left\|\xi\left(t\right)\right\|^{2}$
$\displaystyle+\mathbf{E}\left(\left\|\widehat{f}_{y}\left[t;\varepsilon\right]-\widehat{f}_{y}\left[t\right]\right\|^{2}+\sum_{j=1}^{d}\left\|\widehat{\sigma}_{y}^{j}\left[t;\varepsilon\right]-\widehat{\sigma}_{y}^{j}\left[t\right]\right\|\right)\left\|\mathbf{E}\xi\left(t\right)\right\|^{2}.$

By assumption (A1), we obtain as $\varepsilon\rightarrow 0^{+}$ that

$\displaystyle\left\|\widehat{f}_{x}\left[t;\varepsilon\right]-\widehat{f}_{x}\left[t\right]\right\|^{2}+\sum_{j=1}^{d}\left\|\widehat{\sigma}_{x}^{j}\left[t;\varepsilon\right]-\widehat{\sigma}_{x}^{j}\left[t\right]\right\|^{2}$
$\displaystyle\rightarrow 0,$
$\displaystyle\left\|\widehat{f}_{y}\left[t;\varepsilon\right]-\widehat{f}_{y}\left[t\right]\right\|^{2}+\displaystyle\sum_{j=1}^{d}\left\|\widehat{\sigma}_{y}^{j}\left[t;\varepsilon\right]-\widehat{\sigma}_{y}^{j}\left[t\right]\right\|^{2}$
$\displaystyle\rightarrow 0,$

Therefore, we get that

$\displaystyle 0\leq\lim_{\varepsilon\rightarrow
0}\frac{1}{\varepsilon^{2}}\mathbf{E}\left\|x^{\varepsilon}\left(t+h\right)-\widehat{x}\left(t+h\right)-\xi\left(t+h\right)\right\|^{2}$
$\displaystyle\leq\lim_{\varepsilon\rightarrow
0}\frac{1}{\varepsilon^{2}}\mathbf{E}\left\|x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)-\xi\left(t\right)\right\|^{2}$
$\displaystyle+\lim_{\varepsilon\rightarrow
0}\frac{1}{\varepsilon^{2}}\mathbf{E}\left\|x^{\varepsilon}\left(t-h\right)-\widehat{x}\left(t-h\right)-\xi\left(t-h\right)\right\|^{2}=0.$

Finally, the conclusion is obtained by induction. ∎

######  Lemma 8.

We have the first-order increment of cost functional (2) as follows:

$\displaystyle J\left(u^{\varepsilon}\right)-J\left(\widehat{u}\right)$
$\displaystyle=\mathbf{E}\left\langle\varphi_{x}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),\xi\left(t_{N+1}\right)\right\rangle+\mathbf{E}\left\langle\varphi_{y}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),\mathbf{E}\xi\left(t_{N+1}\right)\right\rangle$
$\displaystyle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\Big{(}\left\langle\widehat{l}_{x}\left[t\right],\xi\left(t\right)\right\rangle+\textbf{E}\left\langle\widehat{l}_{y}\left[t\right],\mathbf{E}\xi\left(t\right)\right\rangle+\left\langle\delta_{t\theta}\widehat{l}_{u}\left[t\right],\varepsilon\Delta
v\right\rangle\Big{)}+o\left(\varepsilon\right),$ (20)

where $u^{\varepsilon}(\cdot)$ is the so-called a spike (or needle) variation
of $\widehat{u}(\cdot)$, defined as follows:

$\displaystyle
u^{\varepsilon}(t)\coloneqq\begin{cases*}\widehat{u}(t),\quad\text{for}\quad
t\neq\theta,\\\ v,\qquad\text{for}\quad t=\theta.\end{cases*}$

###### Proof.

The proof is trivial and based on Lemmas 6 and 7. By the first-order Taylor
expansion, we get

$\displaystyle J\left(u^{\varepsilon}\right)-J\left(\widehat{u}\right)$
$\displaystyle=\mathbf{E}\left(\varphi\left(x^{\varepsilon}\left(t_{N+1}\right),\mathbf{E}x^{\varepsilon}\left(t_{N+1}\right)\right)-\varphi\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right)\right)$
$\displaystyle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left(l^{\varepsilon}\left[t\right]-\widehat{l}\left[t\right]\right)$
$\displaystyle=\mathbf{E}\left\langle\varphi_{x}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),x^{\varepsilon}\left(t_{N+1}\right)-\widehat{x}\left(t_{N+1}\right)\right\rangle$
$\displaystyle+\left\langle\mathbf{E}\varphi_{y}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),\textbf{E}x^{\varepsilon}\left(t_{N+1}\right)-\textbf{E}\widehat{x}\left(t_{N+1}\right)\right\rangle$
$\displaystyle+\mathbf{E}\int_{0}^{1}\left\langle\varphi_{x}\left(\widehat{x}(t_{N+1})+\lambda(x^{\varepsilon}(t_{N+1})-\widehat{x}(t_{N+1})),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),x^{\varepsilon}\left(t_{N+1}\right)-\widehat{x}\left(t_{N+1}\right)\right\rangle
d\lambda$
$\displaystyle+\mathbf{E}\int_{0}^{1}\left\langle\textbf{E}\varphi_{y}(\widehat{x}(t_{N+1}),\textbf{E}\widehat{x}(t_{N+1})+\lambda(\textbf{E}x^{\varepsilon}(t_{N+1})-\textbf{E}\widehat{x}(t_{N+1}))),\textbf{E}x^{\varepsilon}\left(t_{N+1}\right)-\textbf{E}\widehat{x}\left(t_{N+1}\right)\right\rangle
d\lambda$
$\displaystyle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\Big{(}\left\langle\widehat{l}_{x}\left[t\right],x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)\right\rangle+\mathbf{E}\left\langle\widehat{l}_{y}\left[t\right],\textbf{E}x^{\varepsilon}\left(t\right)-\textbf{E}\widehat{x}\left(t\right)\right\rangle\Big{)}+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left\langle\widehat{l}_{u}\left[t\right],\delta_{t\theta}\varepsilon\Delta
v\right\rangle$
$\displaystyle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left\langle\widehat{l}_{x}\left[t;\varepsilon\right]+\mathbf{E}\widehat{l}_{y}\left[t;\varepsilon\right],x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)\right\rangle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left\langle\widehat{l}_{u}\left[t;\varepsilon\right],\delta_{t\theta}\varepsilon\Delta
v\right\rangle.$

Therefore, we get the desired result (8) with $o(\varepsilon)$ which is
defined as follows:

$\displaystyle o(\varepsilon)$
$\displaystyle\coloneqq\mathbf{E}\int_{0}^{1}\left\langle\varphi_{x}\left(\widehat{x}(t_{N+1})+\lambda(x^{\varepsilon}(t_{N+1})-\widehat{x}(t_{N+1})),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),x^{\varepsilon}\left(t_{N+1}\right)-\widehat{x}\left(t_{N+1}\right)\right\rangle
d\lambda$
$\displaystyle+\mathbf{E}\int_{0}^{1}\left\langle\textbf{E}\varphi_{y}(\widehat{x}(t_{N+1}),\textbf{E}\widehat{x}(t_{N+1})+\lambda(\textbf{E}x^{\varepsilon}(t_{N+1})-\textbf{E}\widehat{x}(t_{N+1}))),\textbf{E}x^{\varepsilon}\left(t_{N+1}\right)-\textbf{E}\widehat{x}\left(t_{N+1}\right)\right\rangle
d\lambda$
$\displaystyle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left\langle\widehat{l}_{x}\left[t;\varepsilon\right]+\mathbf{E}\widehat{l}_{y}\left[t;\varepsilon\right],x^{\varepsilon}\left(t\right)-\widehat{x}\left(t\right)\right\rangle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left\langle\widehat{l}_{u}\left[t;\varepsilon\right],\delta_{t\theta}\varepsilon\Delta
v\right\rangle.$

Hence, our conclusion follows. ∎

Now we are in a position to complete the proof of first-order necessary
conditions as stated in Theorem 2. As in the deterministic case the proof is
based on the following identity:

$\mathbf{E}\left\langle
p\left(t_{N+1}\right),\xi\left(t_{N+1}\right)\right\rangle=\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left\langle\Delta
p\left(t\right),\xi\left(t\right)\right\rangle+\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left\langle
p\left(t+h\right),\Delta\xi\left(t\right)\right\rangle.$ (21)

Observe from (3) and (19), respectively that

$\displaystyle\Delta p(t)$
$\displaystyle=-h\Big{(}\mathbf{E}\left\\{\widehat{f}_{x}^{\intercal}\left[t\right]p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\}+\mathbf{E}\left\\{\widehat{f}_{y}^{\intercal}\left[t\right]p\left(t+h\right)\right\\}\Big{)}$
$\displaystyle+\textbf{E}\widehat{l}_{x}[t]+\textbf{E}\left\\{\widehat{l}_{y}[t]\right\\}-\mathbf{E}\sum_{j=1}^{d}\Big{(}\left(\widehat{\sigma}_{x}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)+\mathbf{E}\left\\{\left(\widehat{\sigma}_{y}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)\right\\}\Big{)}$

and

$\displaystyle\Delta\xi(t)$
$\displaystyle=h\Big{(}\widehat{f}_{x}\left[t\right]\xi\left(t\right)+\widehat{f}_{y}\left[t\right]\mathbf{E}\xi\left(t\right)+\delta_{t\theta}\widehat{f}_{u}\left[t\right]\varepsilon\Delta
v\Big{)}$
$\displaystyle+\sum_{j=1}^{d}\mathbf{E}\Big{(}\widehat{\sigma}_{x}^{j}\left[t\right]\xi\left(t\right)+\widehat{\sigma}_{y}^{j}\left[t\right]\mathbf{E}\xi\left(t\right)+\delta_{t\theta}\widehat{\sigma}_{u}^{j}\left[t\right]\varepsilon
v\Big{)}(w_{h}^{j}(t))^{\intercal}$

Then we have that

$\displaystyle\mathbf{E}\left\langle\Delta
p\left(t\right),\xi\left(t\right)\right\rangle$
$\displaystyle=-h\textbf{E}\left\langle\mathbf{E}\left\\{\widehat{f}_{x}^{\intercal}\left[t\right]p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\}+\mathbf{E}\left\\{\widehat{f}_{y}^{\intercal}\left[t\right]p\left(t+h\right)\right\\},\xi\left(t\right)\right\rangle$
$\displaystyle+\textbf{E}\langle\widehat{l}_{x}[t]+\textbf{E}\left\\{\widehat{l}_{y}[t]\right\\},\xi(t)\rangle$
(22)
$\displaystyle-\mathbf{E}\sum_{j=1}^{d}\left\langle\left(\widehat{\sigma}_{x}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)+\mathbf{E}\left\\{\left(\widehat{\sigma}_{y}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)\right\\},\xi\left(t\right)\right\rangle.$

Also, we obtain that

$\displaystyle\mathbf{E}\left\langle
p\left(t+h\right),\Delta\xi\left(t\right)\right\rangle$
$\displaystyle=h\mathbf{E}\left\langle
p\left(t+h\right),\widehat{f}_{x}\left[t\right]\xi\left(t\right)+\widehat{f}_{y}\left[t\right]\mathbf{E}\xi\left(t\right)+\delta_{t\theta}\widehat{f}_{u}\left[t\right]\varepsilon\Delta
v\right\rangle$ (23) $\displaystyle+\sum_{j=1}^{d}\mathbf{E}\left\langle
p\left(t+h\right)w_{h}^{j}\left(t\right),\left(\widehat{\sigma}_{x}^{j}\left[t\right]\xi\left(t\right)+\widehat{\sigma}_{y}^{j}\left[t\right]\mathbf{E}\xi\left(t\right)+\delta_{t\theta}\widehat{\sigma}_{u}^{j}\left[t\right]\varepsilon
v\right)\right\rangle.$

Substituting (5.1) and (23) into (21) yields that

$\displaystyle\mathbf{E}\left\langle
p\left(t_{N+1}\right),\xi\left(t_{N+1}\right)\right\rangle$
$\displaystyle=\sum_{t=t_{0}}^{t_{N}}\textbf{E}\langle\widehat{l}_{x}[t]+\textbf{E}\left\\{\widehat{l}_{y}[t]\right\\},\xi(t)\rangle$
$\displaystyle+\sum_{t=t_{0}}^{t_{N}}\Big{(}h\mathbf{E}\left\langle\delta_{t\theta}p\left(t+h\right),\widehat{f}_{u}\left[t\right]\varepsilon\Delta
v\right\rangle+\sum_{j=1}^{d}\mathbf{E}\left\langle\delta_{t\theta}p\left(t+h\right)w_{h}^{j}\left(t\right),\widehat{\sigma}_{u}^{j}\left[t\right]\varepsilon\Delta
v\right\rangle\Big{)}.$

So, we obtain for all positive integer $N$ that that

$\displaystyle\mathbf{E}\left\langle
p\left(t_{N+1}\right),\xi\left(t_{N+1}\right)\right\rangle$
$\displaystyle=\sum_{t=t_{0}}^{t_{N}}\textbf{E}\langle\widehat{l}_{x}[t]+\textbf{E}\left\\{\widehat{l}_{y}[t]\right\\},\xi(t)\rangle$
$\displaystyle+\sum_{t=t_{0}}^{t_{N}}\varepsilon\mathbf{E}\left\langle
h\widehat{f}_{u}^{\intercal}\left[\theta\right]\mathbf{E}\left\\{p\left(\theta+h\right)\mid\mathfrak{F}_{\theta}\right\\}+\sum_{j=1}^{d}\left(\widehat{\sigma}_{u}^{j}\left[\theta\right]\right)^{\intercal}q^{j}\left(\theta\right),\Delta
v\right\rangle.$

Therefore, by taking into account of the expression of $p(t_{N+1})$ in (3), we
have

$\displaystyle\mathbf{E}\left\langle\varphi_{x}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),\xi\left(t_{N+1}\right)\right\rangle+\mathbf{E}\left\langle\varphi_{y}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),\mathbf{E}\xi\left(t_{N+1}\right)\right\rangle$
$\displaystyle=-\sum_{t=t_{0}}^{t_{N}}\Big{(}\textbf{E}\langle\widehat{l}_{x}[t]+\textbf{E}\left\\{\widehat{l}_{y}[t]\right\\},\xi(t)\rangle-\varepsilon\mathbf{E}\left\langle
h\widehat{f}_{u}^{\intercal}\left[\theta\right]p\left(\theta+h\right)+\left(\widehat{\sigma}_{u}^{j}\left[\theta\right]\right)^{\intercal}q^{j}\left(\theta\right),\Delta
v\right\rangle\Big{)}.$

Thus, according to Lemma 8, it follows from the fact that

$0\leq
J\left(u^{\varepsilon}\right)-J\left(\widehat{u}\right)=-\varepsilon\mathbf{E}\left\langle
H_{u}\left(\theta,\widehat{u}\left(\theta\right)\right),\Delta
v\right\rangle+o\left(\varepsilon\right),\quad\theta\in\mathbb{T},$ (24)

where

$\displaystyle H_{u}(\theta,\widehat{u}(\theta))=\varepsilon\left\langle
h\widehat{f}_{u}^{\intercal}\left[\theta\right]p\left(\theta+h\right)+\left(\widehat{\sigma}_{u}^{j}\left[\theta\right]\right)^{\intercal}q^{j}\left(\theta\right),\Delta
v\right\rangle-\varepsilon\left\langle\widehat{l}_{u}[\theta],\Delta
v\right\rangle.$

Finally, we deduce the variational inequality (4) from (24). The proof is
complete.

##  6 Sufficient conditions for optimality

In this section, we will show that the general controlled discrete-time
stochastic systems formulated earlier, the maximum condition in terms of $H$
function plus some convexity conditions constitute sufficient conditions for
optimality.

We first introduce additional assumptions:

(A5.1)

The function $\varphi$ is convex in $(x,y)$ and functions $f,\sigma,l$ are
convex;

(A5.2)

The Hamiltonian $H$ is convex in $(x,y,v)$;

(A5.3)

The functions $f_{y},\sigma_{y},\varphi_{y},l_{y}$ are non-negative.

######  Theorem 9.

Assume the conditions (A5.1)–(A5.3) are satisfied and let
$\widehat{u}\in\mathfrak{U}_{ad}$ with state trajectory $\widehat{x}$ be given
and such that there exist solutions $(p(\cdot),q^{j}(\cdot)),j=1,\ldots,d$ to
the adjoint equation. Then, if

$H(t,\widehat{x},\widehat{u},p,q)=\inf_{v\in\mathcal{U}(t)}H(t,\widehat{x},v,p,q),$
(25)

then for all $t\in\mathbb{T}$, $\mathbb{P}$-a.s., $\widehat{u}$ is an optimal
control.

######  Remark 10.

By assumption $(A5.2)$, the conditions (4) and (25) are equivalent.

###### Proof.

We use the same short-hand notations which are defined in Section 4. Moreover,
we denote

$\displaystyle H(t)$ $\displaystyle=H(t,x(t),u(t),p(t),q(t)),$
$\displaystyle\widehat{H}(t)$
$\displaystyle=H(t,\widehat{x}(t),\widehat{u}(t),p(t),q(t)).$

Since $\varphi$ is convex, the following inequality holds:

$\displaystyle J\left(\widehat{u}\right)-J\left(u\right)$ (26)
$\displaystyle\leq-\mathbf{E}\left\langle\varphi_{x}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right)-\mathbf{E}\varphi_{y}\left(\widehat{x}\left(t_{N+1}\right),\mathbf{E}\widehat{x}\left(t_{N+1}\right)\right),\widehat{x}\left(t_{N+1}\right)-x\left(t_{N+1}\right)\right\rangle$
$\displaystyle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left(\widehat{l}\left[t\right]-l\left[t\right]\right).$

Therefore, the first term of (26) will be as follows:

$\displaystyle\mathbf{E}\left\langle
p\left(t_{N+1}\right),\widehat{x}\left(t_{N+1}\right)-x\left(t_{N+1}\right)\right\rangle=\sum_{t=t_{0}}^{t_{N}}\Delta\mathbf{E}\left\langle
p\left(t\right),\widehat{x}\left(t\right)-x\left(t\right)\right\rangle$
$\displaystyle=\sum_{t=t_{0}}^{t_{N}}\left(\mathbf{E}\left\langle\Delta
p\left(t\right),\widehat{x}\left(t\right)-x\left(t\right)\right\rangle+\mathbf{E}\left\langle
p\left(t+h\right),\Delta\widehat{x}\left(t\right)-\Delta
x\left(t\right)\right\rangle\right)$
$\displaystyle=-\sum_{t=t_{0}}^{t_{N}}\Big{\langle}\left(h\widehat{f}_{x}^{\intercal}\left[t\right]\right)\mathbf{E}\left\\{p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\}+\mathbf{E}\left\\{\widehat{f}_{y}^{\intercal}\left[t\right]p\left(t+h\right)\right\\}$
$\displaystyle+\sum_{j=1}^{d}\left(\widehat{\sigma}_{x}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)+\mathbf{E}\left\\{\left(\widehat{\sigma}_{y}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)\right\\},\widehat{x}\left(t\right)-x\left(t\right)\Big{\rangle}$
$\displaystyle+\sum_{t=t_{0}}^{t_{N}}\left\langle\left(\widehat{l}_{x}[t]+\textbf{E}\widehat{l}_{y}[t]\right),\widehat{x}\left(t\right)-x\left(t\right)\right\rangle$
$\displaystyle+\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left\langle
p\left(t+h\right),h\left(\widehat{f}\left[t\right]-f\left[t\right]\right)+\sum_{j=1}^{d}\left(\widehat{\sigma}^{j}\left[t\right]-\sigma^{j}\left[t\right]\right)w_{h}^{j}\left(t\right)\right\rangle$
$\displaystyle=-\sum_{t=t_{0}}^{t_{N}}\Big{\langle}\left(h\widehat{f}_{x}^{\intercal}\left[t\right]\right)\mathbf{E}\left\\{p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\}+\mathbf{E}\left\\{\widehat{f}_{y}^{\intercal}\left[t\right]p\left(t+h\right)\right\\}$
$\displaystyle+\sum_{j=1}^{d}\left(\widehat{\sigma}_{x}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)+\mathbf{E}\left\\{\left(\widehat{\sigma}_{y}^{j}\left[t\right]\right)^{\intercal}q^{j}\left(t\right)\right\\},\widehat{x}\left(t\right)-x\left(t\right)\Big{\rangle}$
$\displaystyle+\sum_{t=t_{0}}^{t_{N}}\left\langle\left(\widehat{l}_{x}[t]+\textbf{E}\widehat{l}_{y}[t]\right),\widehat{x}\left(t\right)-x\left(t\right)\right\rangle+\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left(\widehat{H}\left[t\right]-H\left[t\right]\right)+\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left(\widehat{l}\left[t\right]-l\left[t\right]\right),$

where in the last step we have used the definition of Hamiltonian function.
Finally, we differentiate the Hamiltonian and use the convexity of functions
to get for all $t\in\mathbb{T}$, $\mathbb{P}-a.s.$,

$\displaystyle\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left(\widehat{H}\left[t\right]-H\left[t\right]\right)$
$\displaystyle\leq\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left\langle\widehat{H}_{x}\left[t\right]+\widehat{H}_{y}\left[t\right],\widehat{x}\left(t\right)-x\left(t\right)\right\rangle$
$\displaystyle+\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left\langle\widehat{H}_{u}\left[t\right],\left(\widehat{u}\left(t\right)-u\left(t\right)\right)\right\rangle$
$\displaystyle\leq\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left\langle\widehat{H}_{x}\left[t\right]+\widehat{H}_{y}\left[t\right],\widehat{x}\left(t\right)-x\left(t\right)\right\rangle,$

where we have applied
$\left\langle\widehat{H}_{u}\left[t\right],\left(\widehat{u}\left(t\right)-u\left(t\right)\right)\right\rangle\leq
0$ in the last step due to the minimum condition (25). Combining the above
inequalities gives us

$\displaystyle J\left(\widehat{u}\right)-J\left(u\right)$
$\displaystyle\leq\mathbf{E}\left\langle
p\left(t_{N+1}\right),\widehat{x}\left(t_{N+1}\right)-x\left(t_{N+1}\right)\right\rangle+\mathbf{E}\sum_{t=t_{0}}^{t_{N}}\left(\widehat{l}\left[t\right]-l\left[t\right]\right)$
$\displaystyle\leq\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left(\widehat{H}\left[t\right]-H\left[t\right]\right)-\sum_{t=t_{0}}^{t_{N}}\mathbf{E}\left\langle\widehat{H}_{x}\left[t\right]+\widehat{H}_{y}\left[t\right],\widehat{x}\left(t\right)-x\left(t\right)\right\rangle\leq
0$

and thus, $\widehat{u}$ is an optimal. ∎

##  7 Application

Portfolio selection is to seek a best allocation of wealth among a basket of
securities. This model is the foundation of modern finance theory and inspired
literally hundreds of extensions and applications.

Our application is concerned with a discrete-time portfolio selection model
that is formulated as production and consumption choice optimization problems.
The objective is to maximize the expected terminal return and minimize the
variance of the terminal wealth. By putting weights on the two criteria one
obtains a single objective stochastic control problem which is however not in
the standard form due to the variance term involved.

We suppose that an investor is able to invest his wealth to produce some
production, and he can get profit from the production. Denote by $x(t)$ the
capital of this investor at time $t$ and by $v(t)$ the rate of consumption.

Now, we consider some risk in the investment process

$\displaystyle\begin{cases}x\left(t+h\right)=x\left(t\right)+h\left(f\left(x\left(t\right)\right)-\delta
x\left(t\right)\right)-v\left(t\right)+\sigma\left(x\left(t\right)\right)w_{h}\left(t\right),\\\
x\left(0\right)=x_{0}\in\mathbb{R}^{+},\qquad t=t_{0},...,t_{N},\end{cases}$
(27)

where $f\left(x\right)$ is the income production, $\delta$ is the depreciation
rate of the capital, $\sigma\left(x\right)$ denotes the effect of influenced
by the exogenous environment and $w_{h}\left(t\right)$ are the $1$-dimentional
white noises. Our objective is to choose the optimal consumption rate
$v\left(t\right)\geq 0$ to maximize the following functional
$J\left(v\right):$

$J\left(v\right)=\textbf{E}x\left(t_{N+1}\right)+\textbf{E}\sum\limits_{k=0}^{N}l\left(v\left(t_{k}\right)\right),$
(28)

where $x\left(t_{N+1}\right)$ is the capital left over after consumption in
the last period $t_{N+1}$. $l$ is the utility function given by

$l\left(v\right)=\frac{\delta}{\delta-1}v^{1-\frac{1}{\delta}},\ \ \
0<\delta<1.$

Then Hamiltonian function is

$H\left(t,v\right):=l\left(v\right)+hp\left(t+h\right)\left(f\left(x\left(t\right)\right)-\delta
x\left(t\right)-v\right)+q\left(t\right)\sigma\left(x\left(t\right)\right)$

and the corresponding backward stochastic difference equation is

$\left\\{\begin{array}[c]{l}p\left(t\right)=\left(1+hf_{x}\left(x\left(t\right)\right)\right)\mathbf{E}\left\\{p\left(t+h\right)\mid\mathfrak{F}_{t}\right\\}+\sigma_{x}\left(x\left(t\right)\right)q\left(t\right)\\\
q\left(t\right)=\mathbf{E}\left\\{p\left(t+h\right)w_{h}\left(t\right)\mid\mathfrak{F}_{t}\right\\}\\\
p\left(t_{N+1}\right)=1,\ q\left(t_{N}\right)=0.\end{array}\right.$

Assume that $f\left(x\right)=x$ and $\sigma\left(x\right)=\frac{1}{2}x.$ Then
the Hamiltonian function is

$\displaystyle H\left(t,v\right)$
$\displaystyle=\frac{\delta}{\delta-1}v^{1-\frac{1}{\delta}}+hp\left(t+h\right)\left(x\left(t\right)-\delta
x\left(t\right)-v\right)+\frac{1}{2}q\left(t\right)x\left(t\right).$

Solving the equation with respect to $v$:

$\displaystyle
H_{v}\left(t,v\right)=v^{-\frac{1}{\delta}}-hp\left(t+h\right)=0,$

we obtain

$v\left(t\right)=h^{-\delta}p^{-\delta}\left(t+h\right),\quad t\in\mathbb{T}.$

By Theorem 2, we know that $\left\\{v(t)\right\\}$ is the optimal consumption
rate for the optimization problem (27)-(28).

If $N=5$ then we have

$\displaystyle p\left(6h\right)$ $\displaystyle=1,\ q\left(5h\right)=0,$
$\displaystyle p\left(5h\right)$
$\displaystyle=h\left(2-\delta\right)\mathbf{E}\left\\{p\left(6h\right)\mid\mathfrak{F}_{t}\right\\}+\frac{1}{2}0=h\left(2-\delta\right),$
$\displaystyle q\left(4h\right)$ $\displaystyle=0,$ $\displaystyle
p\left(4h\right)$
$\displaystyle=h\left(2-\delta\right)\mathbf{E}\left\\{h\left(2-\delta\right)\mid\mathfrak{F}_{t}\right\\}=h^{2}\left(2-\delta\right)^{2}.$

The trajectory of $v(t)$ with $\delta=0.5$ and $h=0.5$ will be illustrated by
the following figure:

Figure 1: Trajectories of optimal consumption rate $v(t)$

##  8 Conclusions

The core of this paper is to study a new version of maximum principle of
discrete-time mean-field stochastic systems and to provide a proof of main
result stated in Theorem 2. As a consequence of the main theorem constructing
stochastic maximum principle to the discrete-time stochastic systems, by
comparing these results with some existing results in the literature proved
from different point of view.

The main contributions of our work are described in detail as follows:

  * •

We established a new version of the maximum principle for discrete-time mean-
field stochastic optimal control problems and the first-order necessary and
sufficient optimality conditions for discrete-time stochastic optimal control
problems;

  * •

We introduced discrete-time backward (matrix) stochastic equation. Based on
the discrete-time backward stochastic equations, we have obtained necessary
first-order and second-order optimality conditions for the stochastic discrete
optimal control problem (1)-(2);

  * •

Finally, we considered an application as a kind of optimization problem for
production and consumption.

Other related research directions in maximum principle may include the various
relevant topics that may be useful in the optimal control theory of stochastic
differential equations with mean-field type, e.g., one can consider a maximum
principle of a fractional analogue of discrete and continuous time stochastic
differential equations. Our method can also be applied to more complicated
discrete-time stochastic optimal control problems, for example, problems with
delays, terminal constraint problems, and problems with neutral term.

Although there are many articles on the maximum principle of stochastic and
deterministic systems, there still remain many other interesting open problems
concerning their fractional analogues, which can be extended by methods
analogous to those used for fractional derivations of Caputo and Riemann-
Liouville type. To this end, one can consider the method given in [30] to
study the optimal control problem in which a dynamical system is controlled by
a nonlinear Caputo fractional state equation.

## References

  * [1] Dubovitskii, A.Y., Milyutin, A.A., Extremum problems with constraints, Dokl, Akad, Nauk, USSR, English transl. in Sovet Math. Dokl., 4(1963) 2, 452-455.
  * [2] Halkin, H., On the Necessary Conditions for the Optimal Control of Nonlinear Systems, J. Math. Anal. 12, 1–82 (1964)
  * [3] Holtzman, J.M., Convexity and the Maximum Principle for Discrete Systems, IEEE Trans. Automat. Control 11, 30–35 (1966)
  * [4] Holtzman, JM, Halkin, H, Directional convexity and the maximum principle for discrete systems. SIAM J. Control 4(2), 213-275 (1966).
  * [5] Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V., Mishchenko, E.F.: Matematicheskaya teoriya optimal’nykh protsessov (Mathematical Theory of Optimal Processes), Moscow: Nauka, (1969).
  * [6] Rozonoer, LI., Pontryagin’s maximum principle in theory of optimal systems. Avtomatika i Telemekhanika, 20(11), (1959), 1288–1302.
  * [7] Jordan, BW., Polak, E., Theory of a class of discrete optimal control systems, University of California, Berkeley, Electronics Research Lab., Internal Tech. Memo 1,147, 1961.
  * [8] Butkovskii, AG., On necessary and sufficient optimality conditions for impulse control systems. Avtom. Telemeh. 24, 1056-1064 (1963) (in Russian)
  * [9] Gabasov, RF., Kirillova, FM.,The qualitative theory of optimal processes (in Russian). Nauka, Moscow (1971)
  * [10] Song, T., Liu, B., A maximum principle for fully coupled controlled forward-backward stochastic difference systems of mean field-type, Adv. Diff. Equ. 2020, 188 (2020). https://doi.org/10.1186/s13662-020-02640-x.
  * [11] Vinter, RB., Optimality and Sensitivity of Discrete Time Processes, Control Cybernet. 17, 191–211 (1988)
  * [12] Gamkrelidze, RV, Time-optimal processes with restricted phase coordinates, Doklady Akad Nauk, USSR, 3 (1959) 475-478.
  * [13] Mahmudov, NI., Necessary first-order and second-order optimality conditions in discrete-time stochastic systems. J. Optim. Theory Appl. 182, 1001–1018 (2019)
  * [14] Mahmudov, NI., Bashirov, AE., First order and second order necessary conditions of optimality for stochastic systems. Statistics and control of stochastic processes (Moscow, 1995/1996), 283–295, World Sci. Publ., River Edge, NJ, (1997)
  * [15] Yong J., Zhou, X., Stochastic controls: Hamiltonian systems and HJB equations. Springer Verlag, (1999)
  * [16] Rami, M., Chen, X., Zhou, X., Discrete-time indefinite LQ control with state and control dependent noises, Journal of Global Optimization, 23, 245–265 (2002)
  * [17] Zhang, H., Duan, G., Xie, L., Linear quadratic regulation for linear time-varying systems with multiple input delays, Automatica, vol. 42, 1465–1476 (2006)
  * [18] Blot, J., An infinite-horizon stochastic discrete-time pontryagin principle, Nonlinear Analysis: Theory, Methods & Applications, 71, 999–1004 (2009)
  * [19] Wang, H., Zhang, H., Wang, X., Optimal control for stochastic discrete-time systems with multiple input-delays, Proceedings of the 10th World Congress on Intelligent Control and Automation, Beijing, China, July 6-8, 1529–1534 (2012)
  * [20] Mingaleeva, ZT., Shvartsman, IA., Second-order necessary optimality conditions for a discrete optimal control problem. Translation of Differ. Uravn. 50, 1640–1646 (2014)
  * [21] Lin, X., Zhang, W., A maximum principle for optimal control of discrete-time stochastic systems with multiplicative noise. IEEE Trans. Automat. Control 60, 1121–1126 (2015).
  * [22] Mahmudov NI., Maximum Principle for Stochastic Discrete-Time Itô Equations , Brownian Motion: Elements, Dynamics and Applications, Chapter 6, Nova (2015)
  * [23] Buckdahn, R., Djehiche, B., Li, J., A general stochastic maximum principle for SDEs of mean-field type, Applied Mathematics and Optimization, vol. 64, no. 2, pp. 197-216, October 2011.
  * [24] Buckdahn, R., Djehiche, B., Li, J., Peng, S., Mean-field backward stochastic differential equations: a limit approach, Annals of Probability, vol. 37, no. 4, pp. 1524-1565, 2009.
  * [25] Peng, S., A general stochastic maximum principle for optimal control problems, SIAM Journal on Control and Optimization 28 (1990) 966–979.
  * [26] Wu, Z., Zhang, F., Maximum principle for discrete-time stochastic optimal control problem and stochastic game. Mathematical Control & Related fields, 12(2) 2022 475-493.
  * [27] Hu, M., Ji, S., Xue, X., A global stochastic maximum principle for fully coupled forward-backward stochastic systems, SIAM J. Control Optim., 56 (2018), 4309-4335.
  * [28] Yong, J., Optimality variational principle for controlled forward-backward stochastic differential equations with mixed initial-terminal conditions, SIAM J. Control Optim., 48(2010), 4119-4156.
  * [29] Wu, Z., A general maximum principle for optimal control of forward–backward stochastic systems, Automatica, 49(2013), 1473-1480.
  * [30] Yusubov, Sh., Mahmudov, E., Optimality conditions of singular controls for systems with Caputo fractional derivatives, J. Ind. Manag. Optim., (2021) doi: 10.3934/jimo.2021182.