# A Directed Continuous Time Random Walk Model with Jump Length Depending on
Waiting Time

Long Shi1,2, Zuguo Yu1,3, Zhi Mao1, Aiguo Xiao1
1Hunan Key Laboratory for Computation and Simulation in Science and
Engineering and Key
Laboratory of Intelligent Computing and Information Processing of Ministry of
Education,
Xiangtan University, Xiangtan, Hunan 411105, China.
2Institute of Mathematics and Physics, Central South University of Forest and
Technology, Changsha,
Hunan 410004, China.
3School of Mathematical Sciences, Queensland University of Technology, GPO Box
2434,
Brisbane, Q4001, Australia. Corresponding author, e-mail<EMAIL_ADDRESS>or<EMAIL_ADDRESS>
###### Abstract

In continuum one-dimensional space, a coupled directed continuous time random
walk model is proposed, where the random walker jumps toward one direction and
the waiting time between jumps affects the subsequent jump. In the proposed
model, the Laplace-Laplace transform of the probability density function
$P(x,t)$ of finding the walker at position $x$ at time $t$ is completely
determined by the Laplace transform of the probability density function
$\varphi(t)$ of the waiting time. In terms of the probability density function
of the waiting time in the Laplace domain, the limit distribution of the
random process and the corresponding evolving equations are derived.

## 1\. Introduction

The continuous time random walk (CTRW) theory, which was introduced by
Montroll and Weiss [1] to study random walks on a lattice, has been applied
successfully in many fields (see, e.g., the reviews [2-4] and references
therein).

In continuum one-dimensional space, a CTRW process is generated by a sequence
of independent identically distributed (IID) positive waiting times
$T_{1},T_{2},T_{3},\cdots$, and a sequence of IID random jump lengths
$X_{1},X_{2},X_{3},\cdots$. Each waiting time has the same probability density
function (PDF) $\varphi(t),t\geq 0$, and each jump length has the same PDF
$\lambda(x)$ (usually chosen to be symmetric $\lambda(x)=\lambda(-x)$).
Setting $t_{0}=0,t_{n}=T_{1}+T_{2}+\cdots+T_{n}$ for $n\in N$ and
$x_{0}=0,x_{n}=X_{1}+X_{2}+\cdots+X_{n},x(t)=x_{n}$ for $t_{n}\leq t<t_{n+1}$,
we get a microscopic description of the diffusion process [5]. If
$\\{X_{n}\\}$ and $\\{T_{n}\\}$ are independent, the CTRW is called decoupled.
Otherwise it is called coupled CTRW [6]. The decoupled CTRW, which is
completely determined by mutually independent random jump length and random
waiting time, has been widely studied in recent years [3-20].

In some applications it becomes important to consider coupled CTRW [7-8]. The
coupled CTRW can be described by the joint PDF $\phi(x,t)$ of jump length and
waiting time. Because $\phi(x,t)dxdt$ is the probability of a jump to be in
the interval $(x,x+dx)$ in the time interval $(t,t+dt)$, the waiting time PDF
$\varphi(t)=\int_{-\infty}^{+\infty}\phi(x,t)dx$ and the jump length PDF
$\lambda(x)=\int_{0}^{+\infty}\phi(x,t)dt$ can be deduced. Some kinds of
couplings and correlations were proposed in [21-25], where the symmetric jump
length PDF is chosen. For the coupled CTRW, there exist two coupled forms:
$\phi(x,t)=\lambda(x)\varphi(t|x)$ and $\phi(x,t)=\varphi(t)\lambda(x|t)$. The
first coupled form has been studied sufficiently in many literatures [8,
21-23]. The famous model is Lévy walk. Recently, we considered the second
coupled form, discussed the asymptotic behaviors of the coupled jump
probability density function in the Fourier-Laplace domain, and derived the
corresponding fractional diffusion equations from the given asymptotic
behaviors [25].

In this work, we introduce a directed CTRW model with jump length depending on
waiting time (i.e. $\phi(x,t)=\varphi(t)\lambda(x|t),x>0,t>0$). In our model,
the Laplace-Laplace transform [26] of $P(x,t)$ of finding the walker at
position $x$ at time $t$ is completely determined by the Laplace transform of
$\varphi(t)$. Generally, CTRW processes can be categorised by the mean waiting
time $T=\int_{0}^{+\infty}t\varphi(t)dt$ being finite or infinite. Here we
find that the long-time limit distributions of the PDF $P(x,t)$ are a Dirac
delta function for finite $T$ and a beta-like density for infinite $T$, the
corresponding evolving equations are a standard advection equation for finite
$T$ and a pseudo-differential equation with fractional power of coupled space
and time derivative for infinite $T$.

This paper is organized as follows. In section 2, we introduce the basic
concepts of the coupled CTRW. In section 3, a coupled directed CTRW model is
introduced. In section 4, the limit distributions and the corresponding
evolving equations of the coupled directed CTRW model are derived. The
conclusions are given in section 5.

## 2\. The coupled continuous time random walk

Now we recall briefly the general theory of CTRW [3]. Let $\eta(x,t)$ is the
PDF of just having arrived at position $x$ at time $t$. It can be expressed by
$\eta(x^{\prime},t^{\prime})$ (the PDF of just having arrived at position
$x^{\prime}$ at time $t^{\prime}<t$) as:

$\eta(x,t)=\int_{-\infty}^{+\infty}dx^{\prime}\int_{0}^{+\infty}dt^{\prime}\eta(x^{\prime},t^{\prime})\phi(x-x^{\prime},t-t^{\prime})+\delta(x)\delta(t).$
(1)

Then, the PDF $P(x,t)$ with the initial condition $P(x,0)=\delta(x)$ can be
described by the following integral equation [3]

$P(x,t)=\int_{0}^{t}\eta(x,t^{\prime})\omega(t-t^{\prime})dt^{\prime},$ (2)

where $\omega(t)=1-\int_{0}^{t}\varphi(\tau)d\tau$ is the probability of not
having made a jump until time $t$.

Let $\widehat{f}(k)$ and $\widetilde{g}(s)$ be the transforms of Fourier and
Laplace of sufficiently well-behaved (generalized) functions $f(x)$ and $g(t)$
respectively, defined by

$\widehat{f}(k)={\cal
F}\\{f(x);k\\}=\int_{-\infty}^{+\infty}f(x)e^{ikx}dx,\hskip 14.22636ptk\in R,$
(3)

$\widetilde{g}(s)={\cal L}\\{g(t);s\\}=\int_{0}^{+\infty}g(t)e^{-st}dt,\hskip
14.22636pts>s_{0}.$ (4)

After using the Fourier-Laplace transforms and the convolution theorems for
integral equation (2), one can obtain the following famous algebraic relation
[3]

$\widehat{\widetilde{P}}(k,s)=\frac{1-\widetilde{\varphi}(s)}{s}\cdot\frac{1}{1-\widehat{\widetilde{\phi}}(k,s)}.$
(5)

## 3\. A coupled directed CTRW model

In Ref. [23], the author considered a CTRW model with waiting time depending
on the preceding jump length, where the author supposed that the PDF of the
waiting time is a function of a preceding jump length. In that model, the
author introduced a natural ”physiological” analogy: after making a jump one
needs time to rest and recover. The longer the jump distance is, the recovery
and the waiting time needed are longer. This is an interesting hypothetical
physiological example. Motivated by this, we consider a directed CTRW model
with jump length depending on the waiting time and give an analogue
physiological explanation.

A directed CTRW model with jump length depending on the waiting time can be
generated by a sequence of IID positive waiting times
$T_{1},T_{2},T_{3},\cdots$, and a sequence of jumps
$X_{1},X_{2},X_{3},\cdots$. each waiting time has the same PDF
$\varphi(t),t\geq 0$. Every time jump has the same direction and each jump
length has the same conditional PDF $\lambda(x|t),x\geq 0$, which is the PDF
of the random walker making a jump of length $x$ following a waiting time $t$.

A natural assumption is that the jump length is proportional to the waiting
time. So we can take the simplest jump length PDF as
$\lambda(x|t)=\delta(x-vt),v>0$. Without the loss of generality, we take $v=1$
in the following discussion. Setting $t_{0}=0,t_{n}=T_{1}+T_{2}+\cdots+T_{n}$
for $n\in N$ and $x_{0}=0,x_{n}=X_{1}+X_{2}+\cdots+X_{n},x(t)=x_{n}$ for
$t_{n}\leq t<t_{n+1}$, we get a directed CTRW process, where the joint PDF
$\phi(x,t)$ can be expressed by $\phi(x,t)=\varphi(t)\delta(x-t)$. A
physiological explanation can be made as follows: the walker has a random time
for a rest to supplement energy, then makes a jump. The longer the rest time
is, the jump length can be longer.

Since the variable $x$ takes positive values in proposed directed CTRW model,
it is convenient to replace the Fourier transform for variable $x$ in the
formula (5) by the Laplace transform (i.e. $\widetilde{f}(k)={\cal
L}\\{f(x);k\\}=\int_{0}^{+\infty}f(x)e^{-kx}dx$) to obtain the following
Laplace-Laplace relation [26]:

$\widetilde{\widetilde{P}}(k,s)=\frac{1-\widetilde{\varphi}(s)}{s}\cdot\frac{1}{1-\widetilde{\widetilde{\phi}}(k,s)}.$
(6)

Since

$\begin{array}[]{lll}\widetilde{\widetilde{\phi}}(k,s)&=&\int_{0}^{+\infty}dt\int_{0}^{+\infty}\phi(x,t)e^{-kx-
st}dx\\\ \\\
&=&\int_{0}^{+\infty}dt\int_{0}^{+\infty}\varphi(t)\delta(x-t)e^{-kx-st}dx\\\
\\\ &=&\int_{0}^{+\infty}\varphi(t)e^{-(s+k)t}dt\\\ \\\
&=&\widetilde{\varphi}(s+k),\end{array}$ (7)

Eq.(6) is recast into

$\widetilde{\widetilde{P}}(k,s)=\frac{1-\widetilde{\varphi}(s)}{s}\cdot\frac{1}{1-\widetilde{\varphi}(s+k)}.$
(8)

The $n$th ($n=1,2$) moment of $P(x,t)$ is given by

$\begin{array}[]{lll}\langle
x^{n}\rangle(t)&=&\int_{0}^{+\infty}x^{n}(t)P(x,t)dx\\\ \\\
&=&(-1)^{n}\frac{\partial^{n}}{\partial k^{n}}\widetilde{P}(k,t)\mid_{k=0}\\\
\\\ &=&{\cal
L}^{-1}\\{\frac{1-\widetilde{\varphi}(s)}{s}\cdot(-1)^{n}\frac{\partial^{n}}{\partial
k^{n}}\frac{1}{1-\widetilde{\varphi}(s+k)}\mid_{k=0}\\}.\end{array}$ (9)

In the following section, we will study the possible behaviors of $P(x,t)$ and
its $n$th ($n=1,2$) moment.

## 4\. The limit distributions of the coupled directed CTRW model

From Eq.(8), we can see that the Laplace-Laplace transform of PDF $P(x,t)$ is
completely determined by the Laplace transform of the waiting time PDF
$\varphi(t)$. Usually, the random waiting time is characterized by its mean
value $T$. It may be finite or infinite.

For finite mean waiting time $T$, the Laplace transform of $\varphi(t)$ is of
the form

$\widetilde{\varphi}(s)=1-sT+o(s),\hskip 14.22636pts\rightarrow 0.$ (10)

Substituting Eq.(10) into Eq.(8), in the limit $(k,s)\rightarrow(0,0)$, we get
the asymptotic relation

$\widetilde{\widetilde{P}}(k,s)\sim\frac{1-(1-sT)}{s}\cdot\frac{1}{1-(1-(s+k)T)}=\frac{1}{s+k}.$
(11)

After taking the inverse Laplace transforms for Eq.(11) about $k$ and $s$, we
have

$P(x,t)=\delta(x-t).$ (12)

For long times

$\langle x\rangle(t)=t,$ (13) $\langle x^{2}\rangle(t)=t^{2}.$ (14)

From Eq.(11), we get

$s\widetilde{\widetilde{P}}(k,s)-1+k\widetilde{\widetilde{P}}(k,s)=0.$ (15)

Using ${\cal L}\\{\frac{\partial P(x,t)}{\partial
t};s\\}=s\widetilde{P}(x,s)-P(x,0)$, ${\cal L}\\{\frac{\partial
P(x,t)}{\partial x};k\\}=k\widetilde{P}(k,t)-P(0,t)$, initial condition
$P(x,0)=\delta(x)$ and natural boundary conditions, we obtain the partial
differential equation

$\frac{\partial P(x,t)}{\partial t}+\frac{\partial P(x,t)}{\partial x}=0,$
(16)

which is the standard advection equation.

In many applications, one needs to consider a long waiting time (i.e. $T$ is
infinite), it is natural to generalize Eq.(10) to the following form:

$\widetilde{\varphi}(s)=1-s^{\beta}+o(s^{\beta}),\hskip 14.22636pts\rightarrow
0,0<\beta\leq 1.$ (17)

Inserting Eq.(17) into Eq.(8), in the limit $(k,s)\rightarrow(0,0)$, we get
the asymptotic relation

$\widetilde{\widetilde{P}}(k,s)\sim\frac{1-(1-s^{\beta})}{s}\cdot\frac{1}{1-(1-(s+k)^{\beta})}=\frac{s^{\beta-1}}{(s+k)^{\beta}}.$
(18)

After taking the Laplace inverse transform for Eq.(18) about $s$, one has

$\begin{array}[]{lll}\widetilde{P}(k,t)&=&\frac{t^{-\beta}}{\Gamma(1-\beta)}\ast[e^{-kt}\frac{t^{\beta-1}}{\Gamma(\beta)}]\\\
\\\
&=&\int_{0}^{t}e^{-k\tau}\frac{\tau^{\beta-1}(t-\tau)^{-\beta}}{\Gamma(\beta)\Gamma(1-\beta)}d\tau,\end{array}$
(19)

where we use the formulas ${\cal
L}\\{t^{\beta-1};s\\}=\frac{\Gamma(\beta)}{s^{\beta}}$ for $\beta>0$, ${\cal
L}\\{e^{-at}g(t);s\\}=\widetilde{g}(s+a)$ and ${\cal L}\\{(f\ast
g)(t);s\\}=\widetilde{f}(s)\widetilde{g}(s)$.

According to the formula (9) and Eq.(19), for long times, one gets

$\langle x\rangle(t)=\beta t,$ (20) $\langle
x^{2}\rangle(t)=\frac{\beta(\beta+1)}{2}t^{2}.$ (21)

Then taking the Laplace inverse transform for Eq.(19) about $k$, the following
form is obtained

$\begin{array}[]{lll}P(x,t)&=&\int_{0}^{t}\delta(x-\tau)\frac{\tau^{\beta-1}(t-\tau)^{-\beta}}{\Gamma(\beta)\Gamma(1-\beta)}d\tau\\\
\\\
&=&\frac{x^{\beta-1}(t-x)^{-\beta}}{\Gamma(\beta)\Gamma(1-\beta)},\end{array}$
(22)

which is the density of a random variable $tB$, where $B$ has a Beta
distribution with parameters $\beta$ and $1-\beta$.

From Eq.(19), we can also obtain

$(s+k)^{\beta}\widetilde{\widetilde{P}}(k,s)=s^{\beta-1},$ (23)

which leads to the pseudo-differential equation [27-28]

$(\frac{\partial}{\partial t}+\frac{\partial}{\partial
x})^{\beta}P(x,t)=\delta(x)\frac{t^{-\beta}}{\Gamma(1-\beta)}$ (24)

with a coupled space-time fractional derivative operator on the left-hand
side.

Eq.(24) is useful to model flow in porous media and other physical systems
characterized by a link between the waiting time and the jump length.

## 5\. Conclusions

In this work, we introduce a directed CTRW model with jump lengths depending
on waiting times. By the Laplace-Laplace transform technique, we find that the
PDF $P(x,t)$ is determined only by the waiting times PDF $\varphi(t)$. For
finite and infinite mean waiting time, we deduce the limit distributions of
$P(x,t)$ from the asymptotic behaviors of $\varphi(t)$ in the Laplace domain
respectively. The corresponding evolving equations are also derived. For
finite mean waiting time, the limit behavior of the PDF $P(x,t)$ is governed
by a standard advection equation. For infinite mean waiting time, the limit
behavior of the PDF $P(x,t)$ is governed by a pseudo-differential equation
with coupled space-time fractional derivative. We also calculate the first
order moment $\langle x\rangle(t)$ and the second order moment $\langle
x^{2}\rangle(t)$ of $P(x,t)$. An interesting phenomenon is obtained: there
exist the relations $\langle x\rangle(t)\sim t$, $\langle x^{2}\rangle(t)\sim
t^{2}$, whether the mean waiting time is finite or not.

## Acknowledgements

This project was supported by the Natural Science Foundation of China (Grant
Nos. 11371016 and 11271311), the Chinese Program for Changjiang Scholars and
Innovative Research Team in University (PCSIRT) (Grant No. IRT1179), the
Research Foundation of Education Commission of Hunan Province of China (grant
no. 11A122), the Lotus Scholars Program of Hunan province of China.

## References

  * [1] E.W. Montroll, G.H. Weiss, Random walks on lattices II, J. Math. Phys. 6 (1965), 167-181.
  * [2] J.P. Bouchaud, A. Georges, Anomalous diffusion in disordered media: statistical mechanisms models and physical applications. Phys. Rep. 195 (1990), 127-293.
  * [3] R. Metzler, J. Klafter, The random walk s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 339 (2000), 1-77.
  * [4] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A: Math. Gen. 37 (2004), R161-R208.
  * [5] R. Gorenflo, A. Vivoli, F. Mainardi, Discrete and continuous random walk models for space-time fractional diffusion. Nonlinear Dynamics 38 (2004), 101-116.
  * [6] M.M. Meerschaert, E. Nane, Y. Xiao, Correlated continuous time random walks. Stat. Prob. Lett. 79 (2009), 1194-1202.
  * [7] M.F. Shlesinger, J. Klafter, Y.M. Wong, Random walks with infinite spatial and temporal moments. J. Stat. Phys. 27 (1982), 499-512.
  * [8] J. Klafter, A. Blumen, M.F. Shlesinger, Stochastic pathway to anomalous diffusion. Phys. Rev. A 35 (1987), 3081-3085.
  * [9] G.H. Weiss, Aspects and Applications of the Random Walk. North Holland, Amsterdam, 1994.
  * [10] H.E. Roman, P.A. Alemany, Continuous-time random walks and the fractional diffusion equation. J. Phys. A: Math. Gen. 27 (1994), 3407-3410.
  * [11] R. Hilfer, L. Anton, Fractional master equations and fractal time random walks. Phys. Rev. E 51 (1995), R848-R851.
  * [12] E. Scalas, R. Gorenflo, F. Mainardi, Fractional calculus and continuous-time finance. Physica A 284 (2000), 376-384.
  * [13] E. Barkai, CTRW pathways to the fractional diffusion equation. Chem. Phys. 284 (2002), 13-27.
  * [14] R. Hilfer, On fractional diffusion and continuous time random walks. Physica A 329 (2003), 35-40.
  * [15] E. Scalas, R. Gorenflo, F. Mainardi, Uncoupled continuous-time random walks: Solution and limiting behavior of the master equation. Phys. Rev. E 69 (2004), 011107.
  * [16] E. Scalas, The application of continuous-time random walks in finance and economics. Physica A 362 (2006), 225-239.
  * [17] R. Gorenflo, F. Mainardi, A. Vivoli, Continuous-time random walk and parametric subordination in fractional diffusion. Chaos, Solitons and Fractals 34 (2007), 87-103.
  * [18] A.V. Chechkin, M. Hofmann I.M. Sokolov, Continuous-time random walk with correlated waiting time. Phys. Rev. E 80 (2009), 031112.
  * [19] V. Tejedor, R. Metzler, Anomalous diffusion in correlated continuous time random walks. J. Phys. A: Math. Theor. 43 (2010), 082002.
  * [20] K.S. Fa, Uncoupled continuous-time random walk: finite jump length probability density function. J. Phys. A: Math. Theor. 45 (2012), 195002.
  * [21] A. Blumen, G. Zumofen, J. Klafter, Transport aspects in anomalous diffusion: Levy walks. Phys. Rev. A 40 (1989), 3964-3973.
  * [22] G. Zumofen, J. Klafter, Scale-invariant motion in intermittent chaotic systems. Phys. Rev. E 47 (1993), 851-863.
  * [23] V.Y. Zaburdaev, Random walk model with waiting times depending on the preceding jump length. J. Stat. Phys. 123 (2006), 871-881.
  * [24] J. Liu, J.D. Bao, Continuous time random walk with jump length correlated with waiting time. Physica A 392 (2013), 612-617.
  * [25] L. Shi, Z. Yu, Z. Mao, A. Xiao, H. Huang, Space-time fractional diffusion equations and asymptotic behaviors of a coupled continuous time random walk model. Physica A 392 (2013), 5801-5807.
  * [26] R. Gorenflo, F. Mainardi, Laplace-Laplace analysis of the fractional Poisson process, in Analytical Methods of Analysis and Differential Equations (Editor: S.Rogosin),Minsk, 2012, 43-58.
  * [27] P. Becher-Kern, M.M. Meerschaert, H.P. Scheffler, Limit theorems for coupled continuous time random walks. The Annals of Probability 32 (2004), 730-756.
  * [28] A. Jurlewicz, P. Kern, M.M. Meerschaert, H.P. Scheffler, Fractional governing equations for coupled random walks. Comput. Math. Appl. 64 (2012), 3021-3036.


# Solving the Model Unavailable MARE using Q-Learning Algorithm

Fei Yan<EMAIL_ADDRESS>Jie Gao<EMAIL_ADDRESS>Tao Feng
<EMAIL_ADDRESS>Jianxing Liu<EMAIL_ADDRESS>School of Information
Science and Technology, Southwest Jiaotong University, Chengdu, Sichuan, PR
China  National Engineering Laboratory of Integrated Transportation Big Data
Application Technology, Chengdu, China  Department of Control Science and
Engineering, Harbin Institute of Technology, Harbin 150001, PR China

###### Abstract

In this paper, the discrete-time modified algebraic Riccati equation (MARE) is
solved when the system model is completely unavailable. To achieve this,
firstly a brand new iterative method based on the standard discrete-time
algebraic Riccati equation (DARE) and its input weighting matrix is proposed
to solve the MARE. For the single-input case, the iteration can be initialized
by an arbitrary positive input weighting if and only if the MARE has a
stabilizing solution; nevertheless a pre-given input weighting matrix of a
sufficiently large magnitude is used to perform the iteration for the multi-
input case when the characteristic parameter belongs to a specified subset.
Benefit from the developed specific iteration structure, the Q-learning (QL)
algorithm can be employed to subtly solve the MARE where only the system
input/output data is used thus the system model is not required. Finally, a
numerical simulation example is given to verify the effectiveness of the
theoretical results and the algorithm.

###### keywords:

Modified algebraic Riccati equation (MARE); model unavailable; Q-learning (QL)
algorithm.

, , ,

##  1 Introduction

The Riccati equation constitutes an essential component within the framework
of Linear Quadratic Regulator (LQR) and has been extensively investigated [1,
2, 3]. Nevertheless, the practical considerations of imprecise detection,
packet losses and stochastic disturbances cannot be disregarded, which will
introduce significant uncertainties into estimation of target state. Kalman
filter and optimal estimation are two typical techniques to address this
problem, which are closely linked to a modified algebraic Riccati equation
(MARE) that was initially derived decades ago [4]. Since then, the MARE has
gained compelling attention owing to its broad applications in Kalman
filtering [5], intermittent observations [6], network synchronization [7] and
optimal estimator [8], etc. Additionally, the MARE has been further broadened
to systems that are impacted by external disturbances constrained by input
saturation [9] and discrete-time mean-field systems characterized by input
delays[10]. It has been established that the stability characteristics of such
systems are intimately linked to the stabilizing solution of MARE.

Considerable amount of researches have concentrated on the solution of MARE.
In the context where delectability is not assumed, the uniqueness of the
almost stabilizing solution has been discussed in [11]. In an effort to
further tackle this problem, the framework of cone-invariant operators has
been employed, yielding an explicit, necessary, and sufficient condition for
the existence of a mean-square stabilizing solution [12]. Subsequent
contributions have focused on the explicit characterization of the solution.
An analytic solution for the homogeneous MARE for single-input systems [13]
and a closed-form solution in terms of closed-loop poles locations [14] have
been deduced, respectively. However, inherent requirement of solving a set of
Linear Matrix Inequalities (LMIs), which requires accurate system matrices,
fails to adapt these methods to model unavailable scenarios. In scenarios
where the system model is unavailable, QL algorithm [15] has been employed to
resolve such problems. Subsequently, it has been implemented in multi-agent
discrete-time graphical games [16], optimal tracking control [17] and optimal
output regulation [18]. In [19], the QL algorithm is developed to solve a
single DARE, thereby facilitating the computation of the optimal feedback
gain. Then the conditions for the consensus of multi-agent system is
determined and the consensus is achieved when system model is completely
unavailable. However, to best of our knowledge, the development of a model
unavailable algorithm to solve the MARE when system model is completely
unavailable remains an unsolved problem.

In this paper, we aim to develop a new iteration method to solve the MARE, and
then use it to propose a model unavailable algorithm based on the QL
algorithm. In [20], the authors rewrite the MARE as a standard DARE with
specific constraint and thus facilitating an iterative method. However, the
initial conditions to the algorithm still hinge on the LMI problem and the
iteration can only be started by trial and error, which is not suitable to the
model unavailable scenario. In current study, we will address this problem by
proposing a new iterative algorithm to solve the MARE. Specially, in the
single-input case, we show that iteration can be started by an arbitrary
positive input weighting; in the multi-input case, we reduce the selection of
the input weighting matrix within the entire matrix space to a single
sufficient large scalar. Based on this, the QL algorithm can be employed to
solve the MARE when system matrices are completely unavailable.

The structure of this paper is as follows. The problem is established in
section 2. New iterative method to solve MARE is presented in section 3. In
section 4, a model unavailable method based on QL algorithm is developed.
Finally, simulation is given in section 5 and the conclusion is drawn in
section 6.

Notations: ${\lambda_{i}^{u}{\rm{(}}A{\rm{)}}}$ denotes the unstable
eigenvalues of matrix $A$. ${{\lambda_{\max}}{\rm{(}}A{\rm{)}}}$ denotes the
maximum eigenvalue of matrix $A$. ${\mathbb{R}^{n\times m}}$ denotes matrix
space. $I_{n}$ denotes $n$-dimensional identity matrix. Matrix
$A\succ(\succeq)B$ means matrix $A-B$ is positive (semi-) definite.
$A\prec(\preceq)B$ means matrix $A-B$ is negative (semi-) definite.

##  2 Problem formulation

Consider the following discrete-time MARE

$\displaystyle\ X={A^{T}}XA-\gamma{A^{T}}XB{({B^{T}}XB+R)^{-1}}{B^{T}}XA+Q$
(1)

where $A$ and $B$ are system matrices, the characteristic parameter
$\gamma\in(0,1)$, the input weighting matrix $R$ is positive definite and the
state weighting matrix $Q$ is semi-positive definite. Throughout this paper,
it is assumed that $(A,B)$ is stabilizable and $(A,\sqrt{Q})$ is detectable.
Then, it is well known that the MARE (1) has a stabilizing positive definite
solution when the characteristic parameter $\gamma$ is larger than a critical
value $\gamma_{c}$ [21] which is bounded by

$\displaystyle
1-\frac{1}{{{{\max}_{i}}{\rm{|}}\lambda_{i}^{u}{\rm{(}}A{\rm{)}}{{\rm{|}}^{2}}}}\leq{\gamma_{c}}\leq
1-\frac{1}{{\prod_{i}{\rm{|}}\lambda_{i}^{u}{\rm{(}}A{\rm{)}}{{\rm{|}}^{2}}}}.$

In literature [20], the MARE can be simply solved by the LMI technique and
iterative method only when the system matrices are precisely obtained.
However, this paper aims to solve the MARE when the system matrices $A$ and
$B$ are completely unavailable. Therefore, a new iterative method should be
proposed which is especially suitable for developing a model unavailable
algorithm only using the input/output date of the dynamic system.

##  3 Solving MARE by a new iterative method

In this section, we will propose a new iterative method to solve the MARE. For
the single-input case, the iterative method can be initialized by an arbitrary
positive input weighting $R=r>0$ if and only if the characteristic parameter
$\gamma\in(1-({{\prod_{i}{\rm{|}}\lambda_{i}^{u}{\rm{(}}A{\rm{)}}{{\rm{|}}^{2}}}})^{-1},1)$,
i.e., the MARE has a stabilizing solution. For the multi-input case, we
restrict to the case of $\gamma\in(\overline{\gamma},1)$ with
$\overline{\gamma}$ being an estimate of ${\gamma}_{c}$, which facilitates us
to propose the model free algorithm in the next section. It should be pointed
out that the estimate
$\overline{\gamma}=1-({{\prod_{i}{\rm{|}}\lambda_{i}^{u}{\rm{(}}A{\rm{)}}{{\rm{|}}^{2}}}})^{-1}$
for the single-input case.

###  3.1 Single-Input Case

For the single-input case, the control input matrix $B=b\in\mathbb{R}^{n}$,
and the input weighting matrix $R=r$ reduced to a positive scalar. The
following lemma is required for the theoretical development.

######  Lemma 1.

Consider the following DARE

$\displaystyle X={A^{T}}XA-{A^{T}}Xb{({b^{T}}Xb+\sigma)^{-1}}{b^{T}}XA+Q,$ (2)

where $(A,B)$ is stabilizable and $(A,\sqrt{Q})$ is detectable. Defining that
$\beta={b^{T}}Xb/\sigma$, then the scalar $\beta$ converges to its minimum as
$\sigma\to+\infty$ and the infimum is given by

$\displaystyle{\beta_{c}}={\prod_{i}{\rm{|}}\lambda_{i}^{u}{\rm{(}}A{\rm{)}}{{\rm{|}}^{2}}}-1.$

Proof. First, we show that $X/\sigma$ is non-increasing as $\sigma\to\infty$.
Consider the following quadratic performance index

$\displaystyle{J_{\sigma}}(x(0))=x_{0}^{T}{X_{\sigma}}{x_{0}}=\sum\limits_{k=0}^{\infty}{[{x^{T}}(k)Qx(k)+{u^{T}}(k){\sigma}u(k)]}.$

Suppose $0<{\sigma_{1}}<{\sigma_{2}}$, then we will obtain

$\displaystyle J_{{\sigma_{1}}}^{*}({\sigma_{2}}x(0))$ $\displaystyle=$
$\displaystyle x_{0}^{T}({\sigma_{2}}{X_{{\sigma_{1}}}}){x_{0}}$
$\displaystyle=$
$\displaystyle\sum\limits_{k=0}^{\infty}{[{\sigma_{2}}{x^{T}}(k)Qx(k)+{\sigma_{2}}{u^{T}}(k){\sigma_{1}}u(k)]}$
$\displaystyle>$
$\displaystyle\sum\limits_{k=0}^{\infty}{[{\sigma_{1}}{x^{T}}(k)Qx(k)+{\sigma_{1}}{u^{T}}(k){\sigma_{2}}u(k)]}$
$\displaystyle\geq$ $\displaystyle
x_{0}^{T}({\sigma_{1}}{X_{{\sigma_{2}}}}){x_{0}}=J_{{\sigma_{2}}}^{*}({\sigma_{1}}x(0)),$

where ${X_{{\sigma_{i}}}}(i=1,2)$ is the solution of the following DARE

$\displaystyle{X_{{\sigma_{i}}}}={A^{T}}{X_{{\sigma_{i}}}}A-{A^{T}}{X_{{\sigma_{i}}}}b{({b^{T}}{X_{{\sigma_{i}}}}b+{\sigma_{i}})^{-1}}{b^{T}}{X_{{\sigma_{i}}}}A+Q.$

Since this is true for any ${x_{0}}\in{\mathbb{R}^{n}}$, the inequality
${\sigma_{2}}{X_{{\sigma_{1}}}}>{\sigma_{1}}{X_{{\sigma_{2}}}}$ holds, that
is, $({X_{{\sigma_{1}}}}/{\sigma_{1}})>({X_{{\sigma_{2}}}}/{\sigma_{2}})$.
This indicates that ${b^{T}}{X_{\sigma}}b/\sigma$ decreases with respect to
$\sigma$. Due to the fact that ${b^{T}}{X_{\sigma}}b/\sigma>0$, the lower
bound ${\beta_{c}}$ exists as $\sigma\to\infty$.

Then, consider the following scalar

$\displaystyle\delta=\sqrt{\frac{\sigma}{{\sigma+{b^{T}}Xb}}}=\sqrt{\frac{1}{{1+{b^{T}}\frac{X}{\sigma}b}}},$

now we can see that the value $\delta$ reaches its maximum as
${b^{T}}Xb/\sigma$ reaches its minimum, i.e., $\sigma\to\infty$. In the view
of Lemma 5 in [19], the scalar $\delta$ has a supremum

$\displaystyle{\delta_{c}}=\mathop{\lim}\limits_{\sigma\to\infty}\sqrt{\frac{\sigma}{{\sigma+{b^{T}}Xb}}}=\frac{1}{{\prod_{i}{\rm{|}}\lambda_{i}^{u}{\rm{(}}A{\rm{)|}}}}.$
(3)

Solving ${\beta_{c}}$ from equation (3), we will have

$\displaystyle{\beta_{c}}=\mathop{\lim}\limits_{\sigma\to\infty}\frac{{{b^{T}}Xb}}{\sigma}={\prod_{i}{\rm{|}}\lambda_{i}^{u}{\rm{(}}A{\rm{)}}{{\rm{|}}^{2}}}-1.$
(4)

######  Theorem 1.

For the single-input MARE

$\displaystyle X={A^{T}}XA-\gamma{A^{T}}Xb{({b^{T}}Xb+r)^{-1}}{b^{T}}XA+Q,$
(5)

there must exist a positive value ${\omega_{\gamma}}$ and a positive definite
matrix ${X_{\gamma}}$ such that

$\displaystyle{X_{\gamma}}={A^{T}}{X_{\gamma}}A-{A^{T}}{X_{\gamma}}b{({b^{T}}{X_{\gamma}}b+{\omega_{\gamma}})^{-1}}{b^{T}}{X_{\gamma}}A+Q,$
(6a)
$\displaystyle{\omega_{\gamma}}=\frac{1}{\gamma}r+\frac{{1-\gamma}}{\gamma}{b^{T}}{X_{\gamma}}b.$
(6b)

Proof. For any given positive initial value ${\omega_{0}}$, we substitute
$\omega_{0}$ into the following standard algebraic Riccati equation

$\displaystyle X={A^{T}}XA-{A^{T}}Xb{({b^{T}}Xb+{\omega})^{-1}}{b^{T}}XA+Q,$
(7)

which yields a stabilizing solution $X_{0}$. If the pair
$({\omega_{0}},{X_{0}})$ satisfies equation (6a) and (6b) already, then take
${\omega_{\gamma}}={\omega_{0}}$ and ${X_{\gamma}}={X_{0}}$, the issue is
addressed. Otherwise, we consider the following two situations, respectively.

1) The pair $({\omega_{0}},{X_{0}})$ satisfies equation (6a) such that

$\displaystyle{\omega_{0}}>\frac{1}{\gamma}r+\frac{{1-\gamma}}{\gamma}{b^{T}}{X_{0}}b.$
(8)

Then, we construct $\omega_{1}$ which satisfies

$\displaystyle{\omega_{1}}=\frac{1}{\gamma}r+\frac{{1-\gamma}}{\gamma}{b^{T}}{X_{0}}b.$

Obviously, we will have ${\omega_{0}}\geq{\omega_{1}}>0$. Similarly,
substituting $\omega_{1}$ into DARE (7) yields the second stabilizing solution
$X_{1}$. It follows [22] that ${X_{0}}\geq{X_{1}}>0$ when $\omega_{0}$ is
larger than $\omega_{1}$. Then,
$\omega_{2}=r/\gamma+[(1-\gamma){b^{T}}{X_{0}}b]/\gamma$ can also be
constructed and we will have two non-increasing sequences
${\omega_{0}}\geq{\omega_{1}}\geq{\omega_{2}}>0$ and
${X_{0}}\succeq{X_{1}}\succeq{X_{2}}\succ 0$ for the same reason.
Subsequently, two non-increasing sequences
$\\{{\omega_{0}},{\omega_{1}},{\omega_{2}},\ldots\\}$ and
$\\{{X_{0}},{X_{1}},{X_{2}},\ldots\\}$ are easily obtained, where
${\omega_{i}}>0$ and ${X_{i}}=X_{i}^{T}\succ 0$. As a result, these two
sequences converge to ${\omega_{\gamma}}$ and ${X_{\gamma}}$ that satisfies
(6).

2) The pair $({\omega_{0}},{X_{0}})$ satisfies equation (6a) such that

$\displaystyle{\omega_{0}}<\frac{1}{\gamma}r+\frac{{1-\gamma}}{\gamma}{b^{T}}{X_{0}}b.$
(9)

Following the same development way in situation 1), we obtain two non-
decreasing sequences $\\{{\omega_{0}},{\omega_{1}},{\omega_{2}},\ldots\\}$,
$\\{{X_{0}},{X_{1}},{X_{2}},\ldots\\}$, where ${\omega_{i}}>0$ and
${X_{i}}=X_{i}^{T}\succ 0$. Suppose sequence $\omega_{i}$ is bounded, these
two sequences will definitely converge to ${\omega_{\gamma}}$ and
${X_{\gamma}}$. In the converse, i.e., the sequence $\omega_{i}$ is unbounded,
we will have

$\displaystyle\mathop{\lim}\limits_{\omega_{i}\to+\infty}\frac{1}{\gamma}\frac{r}{\omega_{i}}+\frac{{1-\gamma}}{\gamma}\frac{{{b^{T}}Xb}}{\omega_{i}}$
$\displaystyle=$
$\displaystyle\mathop{\lim}\limits_{\omega_{i}\to+\infty}\frac{{1-\gamma}}{\gamma}\frac{{{\gamma_{c}}}}{{1-{\gamma_{c}}}}$
$\displaystyle\leq$
$\displaystyle\frac{{1-{\gamma_{c}}}}{{{\gamma_{c}}}}\frac{{{\gamma_{c}}}}{{1-{\gamma_{c}}}}$
$\displaystyle=$ $\displaystyle 1.$

due to the fact that $\gamma\in({\gamma_{c}},1]$. This indicates that with the
proposed iterative method, there exists $\omega_{i}$ large enough such (8)
holds, i.e., the stabilizing solution ${X_{\gamma}}$ to MARE (5) is certainly
to be obtained.

######  Remark 1.

In [20], the authors have proposed a similar iterative method to solve the
MARE. However, the algorithm must started by an inequality coupling with a
DARE thus the iteration can only be start by trial and error. On the contrary,
Theorem 1 shows that the initial conditions for the input weighting $\omega$
can be an arbitrary given positive scalar to start the iteration. Therefore,
the proposed iterative method is more ease of use. It needs to be emphasized
that the crucial significance of such an improvement of the initialization of
the iteration will establish the solid foundation for proposing the model free
algorithm in the following section.

###  3.2 Multi-Input Case

For the multi-input case, the input matrix $B\in\mathbb{R}^{n\times m}$ is
assumed to be full column rank. The results in Lemma 1 are first extended to
multi-input case and we will give an estimate of the critical value
$\gamma_{c}$.

For simplicity, consider the following multi-input DARE with the input
weighting matrix $\Sigma=\sigma I_{m}$

$\displaystyle X={A^{T}}XA-{A^{T}}XB{({B^{T}}XB+\sigma
I_{m})^{-1}}{B^{T}}XA+Q,$ (10)

then the scalar ${\beta}={\lambda_{\max}}({B^{T}}XB/{\sigma})$ which is
obviously non-increasing with respect to $\sigma$ in the view of Lemma 1. Thus
the infimum of $\beta$ exists which is denoted by
${{\beta}_{c}}=\mathop{\lim}\limits_{\sigma\to\infty}{\lambda_{\max}}({{B^{T}}XB}/{\sigma}).$
Defining that

$\displaystyle\bar{\gamma}=\frac{{{\beta}_{c}}}{1+{\beta}_{c}},$

then we have the following conclusion.

######  Theorem 2.

For the DARE

$\displaystyle X={A^{T}}XA-{A^{T}}XB{({B^{T}}XB+\omega
I_{m})^{-1}}{B^{T}}XA+Q$ (11)

and any given $\gamma\in(\bar{\gamma},1)$, there exists a finite sufficient
large value $\omega_{t}$ such that when $\omega>\omega_{t}$, it holds that

$\displaystyle\omega
I_{m}\succ\frac{1}{\gamma}R+\frac{{1-\gamma}}{\gamma}{B^{T}}XB.$ (12)

Furthermore, the MARE (1) has a stabilizing positive definite solution.

Proof. First we proof that a finite sufficient value $\omega_{t}$ exists for
any given $\gamma\in(\bar{\gamma},1)$ such that equation (12) holds when
$\omega>\omega_{t}$. Then an iterative method similar to single-input case can
be developed when $\omega>\omega_{t}$.

We know that $\beta={\lambda_{\max}}({B^{T}}XB/{\omega})$ is monotonically
decreasing with respect to $\omega$ and the positive definite matrices
$\Omega=\omega I_{m}$ and $X$ satisfy (11). As $\omega\to\infty$,
$\tilde{\gamma}$ converges to its infimum which is
$\bar{\gamma}=\beta_{c}/(1+\beta_{c})$. Then, for any given
$\gamma\in(\bar{\gamma},1)$, there exists a finite sufficient large value
$\omega_{t}$ such that $\gamma_{t}\in(\bar{\gamma},\gamma)$. When taking
$\omega>\omega_{t}$, we will obtain $1>[(1-\gamma)/\gamma]({B^{T}}XB/\omega)$.
Note that $\omega_{t}$ is sufficient large and $\omega>\omega_{t}$, which
means $R/\omega\to 0$. Then we will have

$\displaystyle
1>\frac{1}{\gamma}{\lambda_{\max}}[\frac{R}{\omega}+\frac{(1-\gamma)}{\gamma}\frac{{{B^{T}}XB}}{\omega}],$

which ensures the inequality (12). Then the iterative method is developed as
follows.

Construct matrix $\Omega_{1}$ as in 1 and we have
$\Omega\succeq{\Omega_{1}}\succ 0$. Substituting $\Omega_{1}$ into DARE (11),
we will have the stabilizing solution $X_{1}$. It follows [22] that
$X\succeq{X_{1}}\succ 0$. Then $\Omega_{2}$, $\Omega_{3}$,… can be constructed
in the same way and as a result, two non-increasing sequences
$\\{{\Omega_{0}},{\Omega_{1}},{\Omega_{2}},\ldots\\}$ and
$\\{{X_{0}},{X_{1}},{X_{2}},\ldots\\}$, where ${\Omega_{i}}\succ 0$ and
${X_{i}}=X_{i}^{T}\succ 0$, will be obtained. Obviously, these two sequences
will converge to ${\Omega_{\gamma}}$ and $X_{\gamma}$ such that

$\displaystyle{X_{\gamma}}$
$\displaystyle={A^{T}}{X_{\gamma}}A-{A^{T}}{X_{\gamma}}B{({B^{T}}{X_{\gamma}}B+{\Omega_{\gamma}})^{-1}}{B^{T}}{X_{\gamma}}A+Q,$
$\displaystyle{\Omega_{\gamma}}$
$\displaystyle=\frac{1}{\gamma}R+\frac{{1-\gamma}}{\gamma}{B^{T}}{X_{\gamma}}B.$

Then, the positive definite solution $X_{\gamma}$ is obtained.

Theorem 2 is analogous to Theorem 1 in the previous subsection but with some
constraints specifically for multi-input case. In a special case where matrix
$\Omega$ be with a specialized structure given by $\Omega={\omega}I_{m}$, the
MARE (1) admits an positive definite solution if and only if $\gamma$ is
larger than the estimate value $\bar{\gamma}$ which can be determined via
Q-Learning algorithm when input weighting matrix $R$ is characterized by a
sufficient large magnitude. Given that $\omega_{t}$ is a finite value for any
given $\gamma$ within the interval $(\bar{\gamma},1]$, one can infer that a
scalar $\omega$ of sufficient large magnitude must exists, which permits the
start of iteration. From this perspective, the QL algorithm can be adapted to
solve MARE when the system model is completely unavailable.

##  4 Solving MARE without system model

In this section, the QL algorithm will be employed to derive the solution to
MARE (5) when system matrices $A$ and $b$ are unavailable. In Theorem 1, it is
established that a single standard DARE needs to be considered in each
iteration. This suggests that addressing MARE equates to solve a series of
DAREs. Then, a model free algorithm can be developed as follows.

Define the $\mathbb{Q}$-function as

$\displaystyle\mathbb{Q}(x(k),u(k))$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}x(k)\\\
u(k)\end{array}\right]^{T}\mathbb{H}\left[\begin{array}[]{c}x(k)\\\
u(k)\end{array}\right]$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}x(k)\\\
u(k)\end{array}\right]^{T}\left[\begin{array}[]{cc}\mathbb{H}_{xx}&{\mathbb{H}_{ux}^{T}}\\\
\mathbb{H}_{ux}&\mathbb{H}_{uu}\end{array}\right]\left[\begin{array}[]{c}x(k)\\\
u(k)\end{array}\right],$

where the kernel matrix $\mathbb{H}$ is partitioned into

$\displaystyle{\mathbb{H}_{xx}}=Q+{A^{T}}XA,{\mathbb{H}_{uu}}=r+{b^{T}}Xb,{\mathbb{H}_{ux}}={b^{T}}XA.$

Then, we know that the optimal control is given by

$\displaystyle u(k)=-{\mathbb{H}_{uu}^{-1}\mathbb{H}_{uu}x(k)}.$

And $\gamma_{c}$ can be computed by

$\displaystyle{\gamma_{c}}=1-{\delta_{c}}^{2}=\mathop{\lim}\limits_{r\to\infty}{\sqrt{r{({\mathbb{H}_{uu}})}^{-1}}}.$
(16)

Therefore, the critical value $\gamma_{c}$ can be obtained by means of
$\delta_{c}$. Besides, the solution $X_{\gamma}$ and positive value
$\omega_{\gamma}$ in equation (6) can be expressed by

$\displaystyle{X_{\gamma}}={\mathbb{H}_{xx}}-\mathbb{H}_{ux}^{T}{{({\mathbb{H}_{uu}})}^{-1}}{\mathbb{H}_{ux}},$
(17a)
$\displaystyle{\Omega_{\gamma}}=\frac{{1-{\gamma}}}{\gamma}{\mathbb{H}_{uu}}+r.$
(17b)

In view of this, QL algorithm emerges as a way to address a series of DAREs
when system model is not accessible. The variables $\omega$ and $X$ can be
computed iteratively until conditions in equations (17) are satisfied, i.e.,
the solution has been acquired.

Algorithm 1

Input: A control policy $u(k)=-Kx(k)+n(k)$, where $K$ is optimal feedback gain
stabilizes $(A,b)$, $n(k)$ is an exploring noise; a predefined small threshold
$\varepsilon>0$; an any given initial value $\omega_{0}>0$.

Step 1: Let $Z(k)={[{x^{T}}(k){u^{T}}(k)]^{T}}$ and solve $\mathbb{H}^{t}$
from

$\displaystyle Z(k)$ $\displaystyle=$
$\displaystyle[{x^{T}}(k)Qx(k)]+r{[{u^{t}}(k)]^{T}}{u^{t}}(k)$
$\displaystyle+$ $\displaystyle{Z^{T}}(k+1){{}^{t}}Z(k+1).$

Step 2: Obtain the iterative optimal control ${u^{t+1}}(k)$ as

$\displaystyle{u^{t+1}}(k)=-{(\mathbb{H}_{uu}^{t})^{-1}}\mathbb{H}_{ux}^{t}x(k).$

Step 3: Let $t\leftarrow t+1$, repeat Step 1 and Step 2 until

$\displaystyle||{\mathbb{H}_{t+1}}-{\mathbb{H}_{t}}||<\varepsilon.$

Step 4: Compute $\gamma_{c}$ by

$\displaystyle{\gamma_{c}}=1-{\delta_{c}}^{2}=1-r{{(\mathbb{H}_{uu})}^{-1}}.$

Step 5: If modified parameter $\gamma>\gamma_{c}$, take
$\omega_{n}=\omega_{0}$ as the initialized input weighting matrix. Otherwise,
the algorithm ends.

Step 6: Repeat Step 1 to Step 3 and compute $X_{n}$ by

$\displaystyle
X_{n}={\mathbb{H}_{xx}}-\mathbb{H}_{ux}^{T}{{({\mathbb{H}_{uu}})}^{-1}}{\mathbb{H}_{ux}}.$

Step 7: Obtain the iterative input weighting matrix $\omega_{n+1}$ as

$\displaystyle\omega_{n+1}=\frac{{1-\gamma}}{\gamma}{{{{(\mathbb{H}_{uu}}-\omega_{n}+r)}}+r}.$

Step 8: Let $n\leftarrow n+1$, repeat Step 6 to Step 7 until

$\displaystyle||\omega_{n}-\frac{{1-\gamma}}{\gamma}{\mathbb{H}_{uu}}+r||<\varepsilon.$

Step 9: Obtain the solution $x_{\gamma}=x_{n}$.

End.

######  Remark 2.

As $t\to+\infty$, it is established in [15] that the QL algorithm converges.
And the iterative kernel matrix $\mathbb{H}$ will converge to the true value.
Assuming the input weighting matrix $r$ be sufficiently large at first, the
critical value $\gamma_{c}$ in single-input case can be approximated with an
arbitrary degree of precision. Therefore, the existence of stabilizing
solution for MARE (5) can be verified when system matrices $A$ and $b$ are
entirely unknown.

######  Remark 3.

It is evident that an error is introduced at each iteration step. When
implementing QL algorithm, the error is accumulated iteratively, potentially
leading to the algorithm’s failure to converge particularly in high-
dimensional problems. In light of Theorem 1, it is observed that if the pairs
$({\omega_{l}},{X_{l}}{\rm{)}}$ and $({\omega_{r}},{X_{r}}{\rm{)}}$ satisfy
inequality (9) and (8) respectively, the two sequences
${X_{l}}\prec{\rm{}}{X_{\gamma}}\prec{X_{r}}$,
${\omega_{l}}\leq{\rm{}}{{\rm{\omega}}_{\gamma}}\leq{\omega_{r}}$ holds.
Hence, it suffices to consider the appropriately selected parameter $\omega$.
Only a single DARE (6a) needs to be solved until the condition (6b) is met.
The dichotomy method based on QL algorithm can also be developed which ensures
a higher accuracy.

##  5 Simulation

The system matrices are given by

$\displaystyle A=\left[\begin{array}[]{cc}1&0.5\\\ 0&1\\\ \end{array}\right],\
\ b=\left[\begin{array}[]{c}0\\\ 1\\\ \end{array}\right].$

Obviously, matrix $A$ is neutrally unstable. Set $Q={I_{2}}$, $r=10$,
$\gamma=0.8$. Take ${\omega_{0}}=10$ as the initial value, then we compute
${\omega_{\gamma}}$ and ${X_{\gamma}}$ by the iterative method in Theorem 1
and after 7 iterations, we obtain $X_{\gamma}$ and $\omega_{\gamma}$ as
follows

$\displaystyle{X_{\gamma}}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}5.5114&5.2161\\\ 5.2161&11.7659\\\
\end{array}\right],$ $\displaystyle{\Omega_{\gamma}}$ $\displaystyle=$
$\displaystyle 15.4415$

within the error of ${10^{-7}}$.

When the system matrices are unavailable, we run Algorithm 1 and obtain the
solution to MARE (5). Fig.1 shows the convergence process of
${\gamma_{c}}^{t}$, where ${\gamma_{c}}^{120}=0.1653$. As we can see in Fig.2,
the kernel matrix $\mathbb{H}$ converges to the true value with the error less
than $10^{-10}$ at each iteration. For each iteration, we compute ${X_{n}}$
and the norm of ${\omega_{n}}-[(1-\gamma)/\gamma]{\mathbb{H}_{uu}}-r$ which is
presented in Fig.3. As shown in Fig.3, the error is less than $10^{-12}$ after
60 iterations. Let
$\Delta(X)={A^{T}}XA-\gamma{A^{T}}Xb{({b^{T}}Xb+r)^{-1}}{b^{T}}XA+Q$, the norm
of $X_{n}-\Delta_{n}(X_{n})$ is presented for the $n$th iteration in Fig.4
which can be used to determine the solution. When the difference is $0$,
$X_{n}=\Delta_{n}(X_{n})$ and MARE (5) holds, which means the solution is
obtained by $X_{\gamma}=X_{n}$. After 60 iterations, the solution
${X_{\gamma}}$ and the positive value ${\omega_{\gamma}}$ are obtained as

$\displaystyle{X_{\gamma}}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}5.5114&5.2161\\\ 5.2161&11.7659\\\
\end{array}\right],$ $\displaystyle{\omega_{\gamma}}$ $\displaystyle=$
$\displaystyle 15.4415,$

with the error $\varepsilon={10^{-12}}$ achieved, which means the MARE is
solved by Algorithm 1.

Figure 1: Convergence process of $\gamma_{c}$ Figure 2: Error between matrix
$\mathbb{H}^{t}$ and the true value of $\mathbb{H}$ Figure 3: The norm of
${\omega_{n}}-\frac{{1-\gamma}}{\gamma}{\mathbb{H}_{uu}}-r$ Figure 4: The norm
of $X_{n}-\Delta_{n}(X_{n})$

##  6 Conclusion

In current study, the modified algebraic Ricatti equation (MARE) is solved
when system model is absolutely unavailable. A novel iterative approach is
developed and explicated for both single-input and multiple-input case.
Furthermore, an estimate of the critical value $\gamma_{c}$ is deduced by our
proposed method, which is attainable when system model is unavailable via QL
algorithm. For single-input case, the positive definite solution of MARE can
be obtained from arbitrarily specified input weighting value. Moreover, a
dichotomy method is introduced to guarantee the accuracy of the solution. For
multi-input case, the results are presented for a particular instance,
$\Omega=\omega I_{m}$. Upon the pre-given positive definite matrix $\Omega$
with sufficiently large parameter $\omega$, the problem is tackled when the
characteristic parameter $\gamma$ is larger than the pre-determined estimate
bound $\bar{\gamma}$. Based on this, a model-free algorithm is developed and
by which the MARE is solved.

## References

  * [1] Huaguang Zhang, Tao Feng, Hongjing Liang, and Yanhong Luo.  Lqr-based optimal distributed cooperative design for linear discrete-time multiagent systems.  IEEE Transactions on Neural Networks and Learning Systems, 28(3):599–611, 2017.
  * [2] Junfeng Zhang, Fei Yan, Tao Feng, Tao Deng, and Yue Zhao.  Fastest containment control of discrete-time multi-agent systems using static linear feedback protocol.  Information Sciences, 614:362–373, 2022.
  * [3] Aiguo Wu, Huijie Sun, and Ying Zhang.  Two iterative algorithms for stochastic algebraic riccati matrix equations.  Applied Mathematics and Computation, 339:410–421, 2018.
  * [4] Thomas E. Fortmann, Yaakov Bar-Shalom, Molly Scheffe, and Saul Gelfand.  Detection thresholds for tracking in clutter—a connection between estimation and signal processing.  IEEE Transactions on Automatic Control, 30(3):221–229, 1985.
  * [5] Jianming Zhou, Guoxiang Gu, and Xiang Chen.  Distributed kalman filtering over wireless sensor networks in the presence of data packet drops.  IEEE Transactions on Automatic Control, 64(4):1603–1610, 2019.
  * [6] Bruno Sinopoli, Luca Schenato, Massimo Franceschetti, Kameshwar Poolla, Michael I. Jordan, and Shankar S. Sastry.  Kalman filtering with intermittent observations.  IEEE Transactions on Automatic Control, 49(9):1453–1464, 2004\.
  * [7] Michael Z. Q. Chen, Liangyin Zhang, Housheng Su, and Guanrong Chen.  Stabilizing solution and parameter dependence of modified algebraic riccati equation with application to discrete-time network synchronization.  IEEE Transactions on Automatic Control, 61(1):228–233, 2016.
  * [8] Hong Lin, Hongye Su, Peng Shi, Zhan Shu, Renquan Lu, and Zheng-Guang Wu.  Optimal estimation and control for lossy network: Stability, convergence, and performance.  IEEE Transactions on Automatic Control, 62(9):4564–4579, 2017.
  * [9] Housheng Su, Yanyan Ye, Yuan Qiu, Yang Cao, and Michael Z. Q. Chen.  Semi-global output consensus for discrete-time switching networked systems subject to input saturation and external disturbances.  IEEE Transactions on Cybernetics, 49(11):3934–3945, 2019.
  * [10] Qingyuan Qi, Lihua Xie, and Huanshui Zhang.  Linear–quadratic optimal control for discrete-time mean-field systems with input delay.  IEEE Transactions on Automatic Control, 67(8):3806–3821, 2022.
  * [11] Liangyin Zhang, Michael Z. Q. Chen, and Chanying Li.  On the uniqueness of the almost stabilizing solution of the modified algebraic riccati equation.  In 2016 35th Chinese Control Conference (CCC), pages 7321–7326, 2016.
  * [12] Jianying Zheng and Li Qiu.  Existence of a mean-square stabilizing solution to a modified algebraic riccati equation.  Siam Journal on Control and Optimization, 56(1):367–387, 2018.
  * [13] Liangyin Zhang, Michael Z. Q. Chen, and Chanying Li.  The analytic solutions of the homogeneous modified algebraic riccati equation.  IEEE Transactions on Automatic Control, 62(9):4834–4839, 2017.
  * [14] Alejandro J. Rojas.  Modified algebraic riccati equation closed-form stabilizing solution.  IEEE Access, 9:140667–140675, 2021.
  * [15] Asma Al-Tamimi, Frank L. Lewis, and Murad Abu-Khalaf.  Model-free q-learning designs for linear discrete-time zero-sum games with application to h-infinity control.  Automatica, 43(3):473–481, 2007.
  * [16] Mohammed I. Abouheaf, Frank L. Lewis, Kyriakos G. Vamvoudakis, Sofie Haesaert, and Robert Babuska.  Multi-agent discrete-time graphical games and reinforcement learning solutions.  Automatica, 50(12):3038–3053, 2014.
  * [17] Weijie Sun, Guanyue Zhao, and Yunjian Peng.  Adaptive optimal output feedback tracking control for unknown discrete-time linear systems using a combined reinforcement q-learning and internal model method.  IET Control Theory and Applications., 13(18):3075–3086, 2019.
  * [18] Jialu Fan, Qian Wu, Yi Jiang, Tianyou Chai, and Frank L. Lewis.  Model-free optimal output regulation for linear discrete-time lossy networked control systems.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 50(11):4033–4042, 2020.
  * [19] Tao Feng, Jilie Zhang, Yin Tong, and Huaguang Zhang.  Q-learning algorithm in solving consensusability problem of discrete-time multi-agent systems.  Automatica, 128:109576, 2021.
  * [20] Francisco J. Vargas and Rodrigo A. González.  On the existence of a stabilizing solution of modified algebraic riccati equations in terms of standard algebraic riccati equations and linear matrix inequalities.  IEEE Control Systems Letters, 4(1):91–96, 2020.
  * [21] Luca Schenato, Bruno Sinopoli, Massimo Franceschetti, Kameshwar. Poolla, and S. Shankar Sastry.  Foundations of control and estimation over lossy networks.  Proceedings of the IEEE, 95(1):163–187, 2007.
  * [22] P. Lancaster and L. Rodman.  Algebraic Riccati Equations.  Clarendon Press, Oxford, U.K., 1995.