# On Novel Fixed-Point-Type Iterations with Structure-Preserving Doubling
Algorithms for Stochastic Continuous-time Algebraic Riccati equations

Tsung-Ming Huang  Department of Mathematics, National Taiwan Normal
University, Taipei 116, Taiwan (min@ntnu.edu.tw). Yueh-Cheng Kuo Department of
Applied Mathematics, National University of Kaohsiung, Kaohsiung 811, Taiwan
(yckuo@nuk.edu.tw). Ren-Cang Li  Department of Mathematics, University of
Texas at Arlington, Arlington, USA (rcli@uta.edu). Wen-Wei Lin  Nanjing Center
for Applied Mathematics, Nanjing, China; Department of Applied Mathematics,
National Yang Ming Chiao Tung University, Hsinchu 300, Taiwan.
(wwlin@math.nctu.edu.tw).

(today)

###### Abstract

In this paper we mainly propose efficient and reliable numerical algorithms
for solving stochastic continuous-time algebraic Riccati equations (SCARE)
typically arising from the differential state-dependent Riccati equation
technique from the 3D missile/target engagement, the F16 aircraft flight
control and the quadrotor optimal control etc. To this end, we develop a fixed
point (FP)-type iteration with solving a CARE by the structure-preserving
doubling algorithm (SDA) at each iterative step, called FP-CARE_SDA. We prove
that either the FP-CARE_SDA is monotonically nondecreasing or nonincreasing,
and is R-linearly convergent, with the zero initial matrix or a special
initial matrix satisfying some assumptions. The FP-CARE_SDA (FPC) algorithm
can be regarded as a robust initial step to produce a good initial matrix, and
then the modified Newton (mNT) method can be used by solving the corresponding
Lyapunov equation with SDA (FPC-mNT-Lyap_SDA). Numerical experiments show that
the FPC-mNT-Lyap_SDA algorithm outperforms the other existing algorithms.

##  1 Introduction

The nonlinear dynamics of the stochastic state-dependent control system in
continuous-time subject to multiplicative white noises can be described as

$\displaystyle
dx(t)=A(x)x+B(x)u+\sum_{i=1}^{r}(A_{0}^{i}(x)x+B_{0}^{i}(x)u)dw_{i}(t),$ (1.1)

wher $A(x)$, $A_{0}^{i}(x)\in\mathbb{R}^{n\times n}$, $B(x)$,
$B_{0}^{i}(x)\in\mathbb{R}^{n\times m}$ for $i=1,\ldots,r$, $x(t)$ and $u(t)$
are the state and the control input, $w(t)=[w_{1}(t),\cdots,w_{r}(t)]^{\top}$
is a standard Wiener process satisfying that each $w_{i}(t)$ is a standard
Brownian motion. Under the nonlinear dynamical system (1.1), we consider the
cost functional with respect to the control $u(t)$ with a given initial
$x_{0}$

$\displaystyle
J(t_{0},x_{0};u)=E\left\\{\int_{t_{0}}^{\infty}\begin{bmatrix}x\\\
u\end{bmatrix}^{\top}\begin{bmatrix}Q(x)&L(x)\\\
L^{\top}(x)&R(x)\end{bmatrix}\begin{bmatrix}x\\\ u\end{bmatrix}dt\right\\},$
(1.2)

in which $Q(x)\in\mathbb{R}^{n\times n}$, $L(x)\in\mathbb{R}^{m\times n}$, and
$R(x)>0\in\mathbb{R}^{m\times m}$ with $Q(x)-L(x)R(x)^{-1}L(x)^{\top}\geq 0$.

Recently, the state-dependent Riccati equation (SDRE) [3] generalizes the
well-known linear quadratic regulator (LQR) [23] and attacks broad attention
in nonlinear optimal controls [26, 16, 25]. The SDRE scheme manifests state
and control weighting functions to ameliorate the overall performance [9], as
well as, capabilities and potentials of other performance merits as global
asymptotic stability [21, 19, 4]. In practical applications such as the
differential SDRE with impact angle guidance strategies [20, 24] models a 3D
pursuer/target trajecting tracking or interception engagement, finite-time
SDRE for F16 aircraft flight controls [6], SDRE optimal control design for
quadrotors for enhancing the robustness against unmodeled disturbances [5],
and position/velocity controls for a high-speed vehicle [1]. However, these
application problems are real-world manipulation, the stochastic SDRE (SSDRE)
should indispensably be considered. That is, the goal in SSDRE control is to
minimize the cost function (1.2) and compute the optimal control $u$ at each
fixed/frozen state $x$.

Assume that, with $A_{0}\equiv A\equiv A(x)$, $\\{A_{0}^{i}\equiv
A_{0}^{i}(x),B_{0}^{i}\equiv B_{0}^{i}(x)\\}_{i=1}^{r}$, $Q\equiv Q(x)$,
$R\equiv R(x)$, $L\equiv L(x)$,

(c1) The pair $(\\{A_{0}^{i}\\}_{i=0}^{r},\\{B_{0}^{i}\\}_{i=0}^{r})$ is
stabilizable, i.e., there exists $F\in\mathbb{R}^{m\times n}$ such that the
linear differential equation

$\displaystyle\frac{d}{dt}Z(t)=\mathcal{L}_{F}Z(t)\equiv(A+BF)^{\top}Z+Z(A+BF)+\sum_{i=1}^{r}(A_{0}^{i}+B_{0}^{i}F)^{\top}Z(A_{0}^{i}+B_{0}^{i}F)$

is exponentially stable, i.e., $e^{\mathcal{L}_{F}Z(t-t_{0})}$ is
exponentially stable;

(c2) The pair $(\\{A_{0}^{i}\\}_{i=0}^{r},C)$ is detectable with
$C^{\top}C=Q-LR^{-1}L^{\top}$, that is,
$(\\{A_{0}^{i\top}\\}_{i=0}^{r},\\{C_{i}\\}_{i=0}^{r})$ is stabilizable with
$C_{0}=C$ and $\\{C_{i}\equiv 0\\}_{i=1}^{r}$.

For a fixed state $x$, if both (c1) and (c2) hold, then the stochastic
continuous-time algebraic Riccati equation (SCARE) corresponding to (1.1) and
(1.2)

$\displaystyle A^{\top}X+XA+\sum_{i=1}^{r}A_{0}^{i\top}XA_{0}^{i}+Q$
$\displaystyle-\left(XB+\sum_{i=1}^{r}A_{0}^{i\top}XB_{0}^{i}+L\right)\left(\sum_{i=1}^{r}B_{0}^{i\top}XB_{0}^{i}+R\right)^{-1}\left(B^{\top}X+\sum_{i=1}^{r}B_{0}^{i\top}XA_{0}^{i}+L^{\top}\right)=0,$
(1.3)

has a unique positive semi-definite (PSD) stabilizing solution $X_{*}\equiv
X_{*}(x)$ [11] such that the system
$(A(x)+B(x)F_{X_{*}},A_{0}^{1}(x)+B_{0}^{1}(x)F_{X_{*}},\cdots,A_{0}^{r}(x)+B_{0}^{r}(x)F_{X_{*}})$
is stable with

$\displaystyle
F_{X_{*}}=-\left(\sum_{i=1}^{r}B_{0}^{i}(x)^{\top}X_{*}B_{0}^{i}(x)+R(x)\right)^{-1}\left(B(x)^{\top}X_{*}+\sum_{i=1}^{r}B_{0}^{i}(x)^{\top}X_{*}A_{0}^{i}(x)+L(x)^{\top}\right).$
(1.4)

In fact, $X_{*}$ is a stabilizing solution if and only if the closed-loop
system

$\displaystyle
dx(t)=(A(x)+B(x)F_{X_{*}})x+\sum_{i=1}^{r}(A_{0}^{i}(x)x+B_{0}^{i}(x)F_{X_{*}})dw_{i}(t)$

is exponentially stable.

For convenience for discussion of convergence, the SCARE (1.3) can be
rewritten as

$A^{\top}X+XA+Q+\Pi_{11}(X)\\\
-(XB+L+\Pi_{12}(X))(R+\Pi_{22}(X))^{-1}(XB+L+\Pi_{12}(X))^{\top}=0,$ (1.5)

where

$\displaystyle\Pi(X)$
$\displaystyle\equiv\begin{bmatrix}\Pi_{11}(X)&\Pi_{12}(X)\\\
\Pi_{12}(X)^{\top}&\Pi_{22}(X)\end{bmatrix}\equiv\begin{bmatrix}\sum_{i=1}^{r}A_{0}^{i\top}XA_{0}^{i}&\sum_{i=1}^{r}A_{0}^{i\top}XB_{0}^{i}\\\
\sum_{i=1}^{r}B_{0}^{i\top}XA_{0}^{i}&\sum_{i=1}^{r}B_{0}^{i\top}XB_{0}^{i}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}A_{0}^{1\top}&\cdots&A_{0}^{r\top}\\\
B_{0}^{1\top}&\cdots&B_{0}^{r\top}\\\ \end{bmatrix}(I_{r}\otimes
X)\begin{bmatrix}A_{0}^{1}&B_{0}^{1}\\\ \vdots&\vdots\\\
A_{0}^{r}&B_{0}^{r}\\\ \end{bmatrix}.$ (1.6)

Note that if $X\geq 0$, then $\Pi(X)\geq 0$ and it also holds that

$\displaystyle\Pi(X)\geq\Pi(Y)\geq 0\ \mbox{ for }\ X\geq Y\geq 0.$ (1.7)

The main task of this paper is to develop an efficient, reliable, and robust
algorithm for solving SCARE in (1.3) or (1.5). Several numerical methods, such
as the fixed-point (FP) method [13], Newton (NT) method [10], modified Newton
(mNT) method [8, 12, 15], structure-preserving double algorithms (SDAs) [13],
have been proposed for solving the positive semi-definite (PSD) solutions for
SCARE in (1.3). In fact, based on the FP, NT, mNT and SDA methods, we can
combine these methods to propose various algorithms to efficiently solve SCARE
(1.3) or (1.5). Nevertheless, factors such as the conditions for convergence,
monotonic convergence, convergence speed, accuracy of residuals and
computational cost of iterative algorithms can ultimately determine the
effectiveness, reliability and robustness of a novel algorithm.

We consider the following four algorithms for solving SCARE.

(i) FP-CARE_SDA: We rewrite (1.5) as

$\displaystyle\mathcal{R}(X)\equiv
A_{c}(X)^{\top}X+XA_{c}(X)-XG_{c}(X)X+H_{c}(X)=0$ (1.8)

where $A_{c}(X),G_{c}(X)$ and $H_{c}(X)$ are defined in (2.2b), (2.2c) and
(2.2d), respectively. We frozen $X$ in $A_{c}(X),G_{c}(X)$ and $H_{c}(X)$ as a
fixed-point iteration and solve the CARE (1.8) by SDA [14, 22], called the FP-
CARE_SDA.

Let $\widehat{X}\geq 0$ be a PSD solution of SCARE (1.5). The sequence
$\\{X_{k}\\}_{k=0}^{\infty}$ generated by FP-CARE_SDA with $X_{0}=0$ is
monotonically nondecreasing and R-linearly convergent to a PSD
$0\leq\widehat{X}_{-}\leq\widehat{X}$ (see Theorem 3.3). On the other hand, if
$X_{0}\geq\widehat{X}$ such that $A_{c}(X_{0})-G_{c}(X_{0})X_{0}$ is stable,
and $\mathcal{R}(X_{0})\leq 0$, FP-CARE_SDA generates a monotonically
nonincreasing sequence which converges R-linearly to a PSD
$\widehat{X}_{+}\geq\widehat{X}\geq 0$ (see Theorem 3.4).

(ii) NT-FP-Lyap_SDA: The SCARE (1.5) can be regarded as a nonlinear equation
and solved by Newton’s method in (1.9) which has been derived by [10]

$\displaystyle A_{X_{k}}^{\top}X+XA_{X_{k}}+\Pi_{X_{k}}(X)=-M_{X_{k}},$ (1.9)

where $A_{X_{k}}$, $\Pi_{X_{k}}(X)$ and $M_{X_{k}}$ are given in (4.2b),
(4.2c) and (4.2d), respectively. To solve the Newton’s step (1.9), we frozen
$X$ in $\Pi_{X_{k}}(X)$ as a fixed-point iteration and solve the associated
Lyapunov equation by the special L-SDA in Algorithm 2.

Let $\widehat{X}\geq 0$ be a PSD solution of SCARE (1.5). If there is a
$X_{0}\geq\widehat{X}$ such that $\mathcal{R}_{X_{0}}^{\prime}$ is stable,
then $\\{X_{k}\\}_{k=0}^{\infty}$ generated by NT-FP-Lyap_SDA is monotonically
nonincreasing and quadratically convergent to a PSD
$\widehat{X}_{+}\geq\widehat{X}$ (see [10]). Here
$\mathcal{R}_{X_{0}}^{\prime}$ denotes the Frèchet partial derivative with
respect to $X_{0}$.

(iii) mNT-FP-Lyap_SDA: The method NT-FP-Lyap_SDA proposed in (ii) can be
accelerated by the modified NT (mNT) step [12]. That is, the Lyapunov equation
(1.9) with the frozen term $\Pi_{X_{k}}(X)$ only solved once for approximating
the modified Newton step:

$\displaystyle A_{X_{k}}^{\top}X+XA_{X_{k}}=-\Pi_{X_{k}}(X_{k})-M_{X_{k}}.$

If $X_{0}\geq\widehat{X}$ with $\mathcal{R}(X_{0})\leq 0$ such that
$A_{0}-G_{0}X_{0}$ is stable, and $X_{k}\geq 0$ for all $k\geq 0$ generated by
mNT-FP-Lyap_SDA algorithm, $\\{X_{k}\\}_{k=1}^{\infty}$ is monotonically
nonincreasing and converges to a PSD $\widehat{X}_{+}\geq\widehat{X}\geq 0$
(see Theorem 5.4).

(iv) FP: In [13], the authors used the Möbius transformation to transform the
SCARE (1.5) into a stochastic discrete-time algebraic Riccati equation (SDARE)
and proposed following fixed-point (FP) iteration to solve SDARE:

$\displaystyle X_{1}$ $\displaystyle=H_{\gamma},$ $\displaystyle X_{k+1}$
$\displaystyle=E_{\gamma}^{\top}(X_{k}\otimes
I_{r+1})\left(I_{n(r+1)}+G_{\gamma}(X_{k}\otimes
I_{r+1})\right)^{-1}E_{\gamma}+H_{\gamma},$

where $E_{\gamma}$, $G_{\gamma}$, and $H_{\gamma}$ are defined in (8.1). The
convergence of the fixed-point iteration has been proved in [13].

This paper is organized as follows. We propose the FP-CARE_SDA (FPC) algorithm
for solving SCARE (1.3) in Section 2 and prove its monotonically
nondecreasing/nonincreasing convergence with suitable initial matrices in
Section 3. We propose the mNT-FP-Lyap_SDA algorithm and prove its monotonic
convergences in Sections 4 and 5, respectively. Moreover, for the practical
applications, we proposed FPC-NT-FP-Lyap_SDA and FPC-mNT-FP-Lyap_SDA
algorithms, which are combined FP-CARE_SDA with NT-FP-Lyap_SDA and mNT-FP-
Lyap_SDA algorithms. In Section 6, numerical examples from benchmarks show the
convergence behavior of proposed algorithms. The real-world practical
applications demonstrate the robustness of FP-CARE_SDA and the efficiency of
the FPC-mNT-FP-Lyap_SDA.

##  2 Fixed-Point CARE_SDA method

In this section, we propose a FP-CARE_SDA algorithm to solve (1.5). Let

$\displaystyle L_{c}(X)=L+\Pi_{12}(X),\quad R_{c}(X)=R+\Pi_{22}(X),\quad
Q_{c}(X)=Q+\Pi_{11}(X).$ (2.1)

Then equation (1.5) is equivalent to

$\displaystyle
A^{\top}X+XA-(XB+L_{c}(X))R_{c}(X)^{-1}(XB+L_{c}(X))^{\top}+Q_{c}(X)=0$

or

$\displaystyle A_{c}(X)^{\top}X+XA_{c}(X)-XG_{c}(X)X+H_{c}(X)=0$ (2.2a) with
$\displaystyle A_{c}(X)$ $\displaystyle=A-BR_{c}(X)^{-1}L_{c}(X)^{\top},$
(2.2b) $\displaystyle G_{c}(X)$ $\displaystyle=BR_{c}(X)^{-1}B^{\top},$ (2.2c)
$\displaystyle H_{c}(X)$
$\displaystyle=Q_{c}(X)-L_{c}(X)R_{c}(X)^{-1}L_{c}(X)^{\top}.$ (2.2d)

For a given $X_{k}\geq 0$, we denote

$\displaystyle R_{k}=R_{c}(X_{k}),\quad G_{k}=G_{c}(X_{k}),\quad
A_{k}=A_{c}(X_{k}),\quad H_{k}=H_{c}(X_{k})$ (2.3)

and consider the following CARE

$\displaystyle A_{k}^{\top}X+XA_{k}-XG_{k}X+H_{k}=0,$ (2.4)

where $G_{k}\geq 0$ from (2.3).

######  Proposition 2.1.

If $(A,B)$ is stabilizable, $(A_{k},G_{k})$ is stabilizable.

###### Proof.

Let $y^{\top}(A-BR_{k}^{-1}L_{c}(X_{k}))=\lambda y^{\top}$ with ${\rm
Re}(\lambda)\geq 0$. Suppose $y^{\top}G_{k}=y^{\top}BR_{k}^{-1}B^{\top}=0$.
This implies that $y^{\top}BR_{k}^{-1/2}=0$, i.e., $y^{\top}B=0$. Then we have
$y^{\top}A=\lambda y^{\top}$, ${\rm Re}(\lambda)\geq 0$, which implies that
$y=0$ by stabilzability of $(A,B)$. ∎

######  Proposition 2.2.

Let $C^{\top}C=Q-LR^{-1}L^{\top}$ (a full rank decomposition). Suppose that
${\rm Ker}C\subseteq\left(\bigcap_{i=1}^{r}{\rm
Ker}A_{0}^{i}\right)\bigcap{\rm Ker}L$. If $(C,A)$ is detectable,
$(H_{k},A_{k})$ is detectable and $H_{k}\geq 0$ for $X_{k}\geq 0$.

###### Proof.

Let $A_{k}x=Ax-
BR_{k}^{-1}(L+\sum_{i=1}^{r}(B_{0}^{i})^{\top}X_{k}A_{0}^{i})x=\lambda x$ with
${\rm Re}(\lambda)\geq 0$. If $H_{k}x=0$, then from (2.3) and (2.2d) follows
that

$\displaystyle x^{\top}C^{\top}Cx+x^{\top}\Lambda(X_{k})x=0,$ (2.5)

where
$\Lambda(X)\equiv(Q-C^{\top}C+\Pi_{11}(X))-L_{c}(X)R_{c}(X)^{-1}L_{c}(X)^{\top}$.
From (1), it holds that

$\displaystyle 0\leq\Pi(X_{k})+\begin{bmatrix}Q-C^{\top}C&L\\\
L^{\top}&R\end{bmatrix}=\begin{bmatrix}I&L_{c}(X_{k})R_{c}(X_{k})^{-1}\\\
0&I\end{bmatrix}\begin{bmatrix}\Lambda(X_{k})&0\\\
0&R_{c}(X_{k})\end{bmatrix}\begin{bmatrix}I&0\\\
R_{c}(X_{k})^{-1}L_{c}(X_{k})^{\top}&I\end{bmatrix},$

which implies that $\Lambda(X_{k})\geq 0$ and
$H_{k}=C^{\top}C+\Lambda(X_{k})\geq 0$. Hence, we have $Cx=0$ and
$x^{\top}\Lambda(X_{k})x=0$. Since ${\rm
Ker}C\subseteq\left(\bigcap_{i=1}^{r}{\rm Ker}A_{0}^{i}\right)\bigcap{\rm
Ker}L$, we obtain that $A_{0}^{i}x=0$ for $i=1,\ldots,r$. This means that
$Ax=\lambda x$ with ${\rm Re}(\lambda)\geq 0$. Therefore, $x=0$ by the
detectability of $(C,A)$ and then $(H_{k},A_{k})$ is detectable. ∎

######  Theorem 2.1.

For a given $X_{k}\geq 0$, the CARE of (2.4) has a unique PSD stabilizing
solution $X_{k+1}\geq 0$.

###### Proof.

From Propositions 2.1 and 2.2, as well as the uniqueness and existence theorem
[17] follows that there is a unique stabilizing solution $X_{k+1}\geq 0$ for
(2.4). ∎

It has been shown in [14, 22] that the structure-preserving doubling algorithm
(SDA) can be applied to solve (2.4) with quadratic convergence. In this paper,
we consider using the SDA as an implicit fixed-point iteration, called FP-
CARE_SDA as stated in Algorithm 1, for solving the SCARE (1.5).

Algorithm 1 FP-CARE_SDA method for solving SCARE (1.5)

0: $A,Q\in\mathbb{R}^{n\times n}$, $B,L\in\mathbb{R}^{n\times m}$,
$R=R^{\top}\in\mathbb{R}^{m\times m}$, $A_{0}^{i}\in\mathbb{R}^{n\times n}$,
$B_{0}^{i}\in\mathbb{R}^{n\times m}$ for $i=1,\ldots,r$, and a tolerance
$\varepsilon$.

0: Solution $X$.

1: Set $k=0$ and $\delta=\infty$.

2: Choose an initial $X_{0}^{\top}=X_{0}$.

3: while $\delta>\varepsilon$ do

4: Set $L_{c}=L+\Pi_{12}(X_{k})$, $R_{c}=\Pi_{22}(X_{k})+R$ and
$Q_{c}=\Pi_{11}(X_{k})+Q$.

5: Set $A_{c}=A-BR_{c}^{-1}L_{c}^{\top}$, $G_{c}=BR_{c}^{-1}B^{\top}$ and
$H_{c}=Q_{c}-L_{c}R_{c}^{-1}L_{c}^{\top}$.

6: Use SDA to solve the stabilizing solution $X_{k+1}$ of the CARE
$\displaystyle A_{c}^{\top}X+XA_{c}-XG_{c}X+H_{c}=0.$

7: Compute the normalized residual $\delta$ of (1.5) with $X=X_{k+1}$.

8: Set $k=k+1$.

9: end while

##  3 Convergence Analysis for FP-CARE_SDA Method

In this section, we will study the convergence analysis of the FP-CARE_SDA
algorithm. Let $\widehat{X}$ be a PSD solution of SCARE (1.3). With different
initial conditions of $X_{0}$, the sequence $\\{X_{k}\\}_{k=0}^{\infty}$
generated by FP-CARE_SDA is monotonically nondecreasing/nonincreasing to
$\widehat{X}$ and R-linearly convergent to some PSD solution
$\widehat{X}_{-}/\widehat{X}_{+}$, respectively, with
$\widehat{X}_{+}\geq\widehat{X}\geq\widehat{X}_{-}\geq 0$.

Let

$\displaystyle\Gamma(X)$
$\displaystyle\equiv\left[\begin{array}[]{cc|c}0&-A&-B\\\
-A^{\top}&Q_{c}(X)&L_{c}(X)\\\
\hline\cr-B^{\top}&L_{c}(X)^{\top}&R_{c}(X)\end{array}\right]\equiv\left[\begin{array}[]{c|c}\Gamma_{11}(X)&\Gamma_{12}(X)\\\
\hline\cr\Gamma_{21}(X)^{\top}&R_{c}(X)\end{array}\right]\in\mathbb{R}^{(2n+m)\times(2n+m)},$
(3.6)

where $Q_{c}(X)$, $L_{c}(X)$, and $R_{c}(X)$ are given in (2.1). Then the
Schur complement of the block matrix $\Gamma(X)$ has the form

$\displaystyle\Omega(X)=$ $\displaystyle\begin{bmatrix}0&-A\\\
-A^{\top}&Q_{c}(X)\end{bmatrix}-\begin{bmatrix}-B\\\
L_{c}(X)\end{bmatrix}R_{c}(X)^{-1}\begin{bmatrix}-B^{\top}&L_{c}(X)^{\top}\end{bmatrix}$
$\displaystyle=$
$\displaystyle\begin{bmatrix}-BR_{c}(X)^{-1}B^{\top}&-(A-BR_{c}(X)^{-1}L_{c}(X)^{\top})\\\
-(A^{\top}-L_{c}(X)R_{c}(X)^{-1}B^{\top})&Q_{c}(X)-L_{c}(X)R_{c}(X)^{-1}L_{c}(X)^{\top}\end{bmatrix}$
$\displaystyle=$ $\displaystyle\begin{bmatrix}-G_{c}(X)&-A_{c}(X)\\\
-A_{c}(X)^{\top}&H_{c}(X)\end{bmatrix},$ (3.7)

where $A_{c}(X)$, $G_{c}(X)$ and $H_{c}(X)$ in (2.2) are the coefficient
matrices of CARE (2.2a).

######  Lemma 3.1.

Let $\Omega(X)$ be defined in (3). Then

$\displaystyle\Omega(X)\geq\Omega(Y)\text{ for }X\geq Y\geq 0.$ (3.8)

###### Proof.

Suppose that $X\geq Y\geq 0$. From (3.6), (2.1) and (1), we have

$\displaystyle\Gamma(X)$ $\displaystyle=\left[\begin{array}[]{c|cc}0&-A&-B\\\
\hline\cr-A^{\top}&Q&L\\\
-B^{\top}&L^{\top}&R\end{array}\right]+\left[\begin{array}[]{ccc}0&\cdots&0\\\
\hline\cr A_{0}^{1\top}&\cdots&A_{0}^{r\top}\\\
B_{0}^{1\top}&\cdots&B_{0}^{r\top}\end{array}\right](I_{r}\otimes
X)\left[\begin{array}[]{c|cc}0&A_{0}^{1}&B_{0}^{1}\\\ \vdots&\vdots&\vdots\\\
0&A_{0}^{r}&B_{0}^{r}\end{array}\right]$
$\displaystyle\geq\left[\begin{array}[]{c|cc}0&-A&-B\\\
\hline\cr-A^{\top}&Q&L\\\
-B^{\top}&L^{\top}&R\end{array}\right]+\left[\begin{array}[]{ccc}0&\cdots&0\\\
\hline\cr A_{0}^{1\top}&\cdots&A_{0}^{r\top}\\\
B_{0}^{1\top}&\cdots&B_{0}^{r\top}\end{array}\right](I_{r}\otimes
Y)\left[\begin{array}[]{c|cc}0&A_{0}^{1}&B_{0}^{1}\\\ \vdots&\vdots&\vdots\\\
0&A_{0}^{r}&B_{0}^{r}\end{array}\right]=\Gamma(Y).$

Now, we show that (3.8). For any $\textbf{w}\in\mathbb{R}^{2n}$, since
$\Gamma(X)\geq\Gamma(Y)$, we have

$\displaystyle\Gamma_{\textbf{w}}(X)\equiv\left[\begin{array}[]{cc}\textbf{w}^{\top}\Gamma_{11}(X)\textbf{w}&\textbf{w}^{\top}\Gamma_{12}(X)\\\
\Gamma_{12}(X)^{\top}\textbf{w}&R_{c}(X)\end{array}\right]\geq\left[\begin{array}[]{cc}\textbf{w}^{\top}\Gamma_{11}(Y)\textbf{w}&\textbf{w}^{\top}\Gamma_{12}(Y)\\\
\Gamma_{12}^{\top}(Y)\textbf{w}&R_{c}(Y)\end{array}\right]\equiv\Gamma_{\textbf{w}}(Y).$
(3.13)

It is easily seen that $\textbf{w}^{\top}\Omega(X)\textbf{w}$ and
$\textbf{w}^{\top}\Omega(Y)\textbf{w}$ are Schur complement of
$\Gamma_{\textbf{w}}(X)$ and $\Gamma_{\textbf{w}}(Y)$, respectively. Suppose
that $\alpha\in\mathbb{R}_{+}$ such that

$\displaystyle\textbf{w}^{\top}\Omega(Y)\textbf{w}+\alpha>0.$

Let

$\displaystyle\Gamma_{\textbf{w}}^{\alpha}(X)=\Gamma_{\textbf{w}}(X)+\left[\begin{array}[]{cc}\alpha&0\\\
0&0\end{array}\right],\ \
\Gamma_{\textbf{w}}^{\alpha}(Y)=\Gamma_{\textbf{w}}(Y)+\left[\begin{array}[]{cc}\alpha&0\\\
0&0\end{array}\right].$ (3.18)

Note that $\textbf{w}^{\top}\Omega(Y)\textbf{w}+\alpha$ is the Schur
complement of $\Gamma_{\textbf{w}}^{\alpha}(Y)$. Since $R_{c}(Y)>0$ and
$\textbf{w}^{\top}\Omega(Y)\textbf{w}+\alpha>0$, it follows from (3.13) that
$\Gamma_{\textbf{w}}^{\alpha}(X)\geq\Gamma_{\textbf{w}}^{\alpha}(Y)>0$. Hence,

$\displaystyle
0<\Gamma_{\textbf{w}}^{\alpha}(X)^{-1}\leq\Gamma_{\textbf{w}}^{\alpha}(Y)^{-1}.$
(3.19)

From (3.18), we obtain that the $(1,1)$ entries of
$\Gamma_{\textbf{w}}^{\alpha}(X)^{-1}$ and
$\Gamma_{\textbf{w}}^{\alpha}(Y)^{-1}$ are
$(\textbf{w}^{\top}\Omega(X)\textbf{w}+\alpha)^{-1}$ and
$(\textbf{w}^{\top}\Omega(Y)\textbf{w}+\alpha)^{-1}$, respectively. It follows
from (3.19) that
$0<(\textbf{w}^{\top}\Omega(X)\textbf{w}+\alpha)^{-1}\leq(\textbf{w}^{\top}\Omega(Y)\textbf{w}+\alpha)^{-1}$
for any $\textbf{w}\in\mathbb{R}^{2n}$. This implies that (3.8) holds. ∎

######  Remark 1.

Let $\Omega(X)$ be defined in (3). The SCARE in (2.2a) can be represented as

$\displaystyle\begin{bmatrix}X&-I\end{bmatrix}\Omega(X)\begin{bmatrix}X\\\
-I\end{bmatrix}=0.$ (3.20)

Let $\\{X_{k}\\}_{k=0}^{\infty}$ be a sequence generated by Algorithm 1.
Suppose that $X_{k-1}\geq 0$. From Theorem 2.1, it holds that the CARE

$\displaystyle A_{k-1}^{\top}X_{k}+X_{k}A_{k-1}-X_{k}G_{k-1}X_{k}+H_{k-1}=0,$
(3.21)

has a PSD solution $X_{k}\geq 0$. That is,

$\displaystyle\begin{bmatrix}X_{k}&-I\end{bmatrix}\Omega(X_{k-1})\begin{bmatrix}X_{k}\\\
-I\end{bmatrix}=0.$ (3.22)

We now state a well-known result [18] to show the monotonicity of
$\\{X_{k}\\}_{k=0}^{\infty}$.

######  Lemma 3.2.

[18] Let $A,Q\in\mathbb{R}^{n\times n}$ with $A$ stable and $Q\leq 0$. Then
the solution of the Lyapunov equation $A^{\top}X+XA=Q$ is PSD.

######  Theorem 3.3.

Let $\widehat{X}\geq 0$ be a solution of SCARE (1.5). Suppose that
$\\{X_{k}\\}_{k=0}^{\infty}$ is a sequence generated by Algorithm 1 with
$X_{0}=0$. Then

(i)

$X_{k}\leq X_{k+1}$, $X_{k}\leq\widehat{X}$, $\mathcal{R}(X_{k})\geq 0$, and
$A_{k}-G_{k}\widehat{X}$ is stable for each $k\geq 0$.

(ii)

$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{-}$, where
$\widehat{X}_{-}\leq\widehat{X}$ is a solution of SCARE (1.5).

(iii)

$\sigma\left(A_{c}(\widehat{X}_{-})-G_{c}(\widehat{X}_{-})\widehat{X}_{-}\right)\subset\mathbb{C}_{-}\bigcup
i\mathbb{R}$.

###### Proof.

(i). We prove the assertion by induction for each $k\geq 0$. For $k=0$, we
have $0=X_{0}\leq\widehat{X}$ and $\mathcal{R}(X_{0})=H_{0}\geq 0$. Since
$X_{1}$ is the stabilizing PSD solution of

$\displaystyle A_{0}^{\top}X+XA_{0}-XG_{0}X+H_{0}=0,$

where $(A_{0},G_{0})$ is stabilizable and $(H_{0},A_{0})$ is detectable, we
have $X_{1}\geq X_{0}=0$. Now, we claim that $A_{0}-G_{0}\widehat{X}$ is
stable. From Lemma 3.1 and using the fact that $X_{0}\leq\widehat{X}$, we have

$\displaystyle(A_{0}^{\top}-\widehat{X}G_{0})\widehat{X}+\widehat{X}(A_{0}-G_{0}\widehat{X})$
$\displaystyle=\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{0})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}-H_{0}-\widehat{X}^{\top}G_{0}\widehat{X}$ $\displaystyle\leq-
H_{0}-\widehat{X}^{\top}G_{0}\widehat{X}.$ (3.23)

Suppose that $A_{0}-G_{0}\widehat{X}$ is not stable. Then there exist $y\neq
0$ and Re$(\lambda)\geq 0$ such that $(A_{0}-G_{0}\widehat{X})y=\lambda y$.
From (3), we have

$\displaystyle 2{\rm
Re}(\lambda)y^{*}\widehat{X}y\leq-y^{*}H_{0}y-y^{*}\widehat{X}^{\top}G_{0}\widehat{X}y.$

Since $\widehat{X},H_{0},G_{0}\geq 0$ and Re$(\lambda)\geq 0$, we obtain that
$H_{0}y=0$ and $G_{0}\widehat{X}y=0$. This implies that $A_{0}y=\lambda y$
which is a contradiction because of the detectabity of $(H_{0},A_{0})$. Hence,
$A_{0}-G_{0}\widehat{X}$ is stable. We now assume that the assertion (i) is
true for $k\geq 0$ and $X_{k+2}$ is the stabilizing solution of the CARE

$\displaystyle\begin{bmatrix}X&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X\\\
-I\end{bmatrix}=0.$

Then we have $A_{k+1}-G_{k+1}X_{k+2}$ is stable and

$\displaystyle(A_{k+1}^{\top}-X_{k+2}G_{k+1})(X_{k+2}-X_{k+1})+(X_{k+2}-X_{k+1})(A_{k+1}-G_{k+1}X_{k+2})$
$\displaystyle=$ $\displaystyle
A_{k+1}^{\top}X_{k+2}+X_{k+2}A_{k+1}-X_{k+2}G_{k+1}X_{k+2}-(A_{k+1}^{\top}X_{k+1}+X_{k+1}A_{k+1}-X_{k+1}G_{k+1}X_{k+1})$
$\displaystyle-(X_{k+2}-X_{k+1})G_{k+1}(X_{k+2}-X_{k+1})$ $\displaystyle=$
$\displaystyle-\mathcal{R}(X_{k+1})-(X_{k+2}-X_{k+1})G_{k+1}(X_{k+2}-X_{k+1}).$

From Lemma 3.1 and using the fact that $X_{k}\leq X_{k+1}$, we have

$\displaystyle
0=\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}\leq\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}=\mathcal{R}(X_{k+1}).$

Hence, $X_{k+1}\leq X_{k+2}$ by Lemma 3.2. Using the fact that
$X_{k}\leq\widehat{X}$ and Lemma 3.1, we have

$\displaystyle(A_{k}^{\top}-\widehat{X}G_{k})(\widehat{X}-X_{k+1})+(\widehat{X}-X_{k+1})(A_{k}-G_{k}\widehat{X})$
$\displaystyle=\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}-(\widehat{X}-X_{k+1})G_{k}(\widehat{X}-X_{k+1})\leq-(\widehat{X}-X_{k+1})G_{k}(\widehat{X}-X_{k+1}).$

This leads to $X_{k+1}\leq\widehat{X}$ because $A_{k}-G_{k}\widehat{X}$ is
stable. Now, we claim that $A_{k+1}-G_{k+1}\widehat{X}$ is stable. After
calculation, it holds that

$\displaystyle(A_{k+1}^{\top}-\widehat{X}G_{k+1})(\widehat{X}-X_{k+1})+(\widehat{X}-X_{k+1})(A_{k+1}-G_{k+1}\widehat{X})$
$\displaystyle=\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}-\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}-(\widehat{X}-X_{k+1})G_{k+1}(\widehat{X}-X_{k+1}).$ (3.24a)
Similarly, using the fact that $X_{k}\leq X_{k+1}\leq\widehat{X}$ and Lemma
3.1, we have
$\displaystyle\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}\leq 0\ \text{ and }\
-\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}\leq 0.$ (3.24b)

Assume that $A_{k+1}-G_{k+1}\widehat{X}$ is not stable. Then there exist
$y\neq 0$ and $\mbox{Re}(\lambda)\geq 0$ such that
$(A_{k+1}-G_{k+1}\widehat{X})y=\lambda y$. From (3.24) and
$X_{k+1}\leq\widehat{X}$, we obtain that $0\leq
2\mbox{Re}(\lambda)y^{*}(\widehat{X}-X_{k+1})y$ and hence,

$\displaystyle 0$
$\displaystyle=y^{*}\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}y$
$\displaystyle=y^{*}\left[(A_{k+1}^{\top}-\widehat{X}G_{k+1})\widehat{X}+\widehat{X}(A_{k+1}-G_{k+1}\widehat{X})+\widehat{X}G_{k+1}\widehat{X}+H_{k+1}\right]y$
$\displaystyle=2\mbox{Re}(\lambda)y^{*}\widehat{X}y+y^{*}\widehat{X}G_{k+1}\widehat{X}y+y^{*}H_{k+1}y.$

Since $\widehat{X},H_{k+1},G_{k+1}\geq 0$, we obtain that $H_{k+1}y=0$ and
$G_{k+1}\widehat{X}y=0$. This means that

$\displaystyle\lambda y=(A_{k+1}-G_{k+1}\widehat{X})y=A_{k+1}y,\quad y\neq
0,\quad\mbox{Re}(\lambda)\geq 0,\quad\mbox{ and }\quad H_{k+1}y=0,$

which contradicts the detectability of $(H_{k+1},A_{k+1})$. Therefore,
$A_{k+1}-G_{k+1}\widehat{X}$ is stable. The induction process is complete.

(ii). Since the sequence $\\{X_{k}\\}$ is monotonically nondecreasing and
bounded above by $\widehat{X}$, we have
$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{-}$, where
$\widehat{X}_{-}\leq\widehat{X}$ is a solution of SCARE (1.5).

(iii). By (i), all $A_{k}-G_{k}\widehat{X}_{-}$ are stable, by continuity
$\sigma\left(A_{c}(\widehat{X}_{-})-G_{c}(\widehat{X}_{-})\widehat{X}_{-}\right)\subset\mathbb{C}_{-}\bigcup
i\mathbb{R}$. ∎

From another perspective, we can prove that $\\{X_{k}\\}_{k=0}^{\infty}$ is
monotonically nonincreasing if the initial $X_{0}$ satisfies some strict
conditions.

######  Theorem 3.4.

Let $\widehat{X}\geq 0$ be a solution of SCARE (1.5). Suppose that
$\\{X_{k}\\}_{k=0}^{\infty}$ is a sequence generated by Algorithm 1. Let
$X_{0}\geq\widehat{X}$ such that $A_{0}-G_{0}X_{0}$ is stable and
$\mathcal{R}(X_{0})\leq 0$. Then

(i)

$X_{k}\geq X_{k+1}$, $X_{k}\geq\widehat{X}$, $\mathcal{R}(X_{k})\leq 0$ for
each $k\geq 0$.

(ii)

$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{+}$, where
$\widehat{X}_{+}\geq\widehat{X}$ is a solution of SCARE (1.5).

(iii)

$\sigma\left(A_{c}(\widehat{X}_{+})-G_{c}(\widehat{X}_{+})\widehat{X}_{+}\right)\subset\mathbb{C}_{-}\bigcup
i\mathbb{R}$.

###### Proof.

(i). We prove the assertion by mathematical induction on $k\geq 0$. First, we
show that if $X_{k-1}\geq\widehat{X}$, $X_{k}\geq\widehat{X}$. From (3.22), we
have

$\displaystyle(A_{k-1}^{\top}-X_{k}G_{k-1})(\widehat{X}-X_{k})+(\widehat{X}-X_{k})(A_{k-1}-G_{k-1}X_{k})$
$\displaystyle=$ $\displaystyle\
(\widehat{X}-X_{k})G_{k-1}(\widehat{X}-X_{k})+\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{k-1})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}.$ (3.25)

Since $X_{k-1}\geq\widehat{X}$ and $\widehat{X}$ is a solution of SCARE (1.5),
it follows from Lemma 3.1 that

$\displaystyle\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{k-1})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}\geq\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(\widehat{X})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}=0.$

Since $G_{k-1}\geq 0$ and $A_{k-1}-G_{k-1}X_{k}$ is stable, it follows from
(3) and Lemma 3.2 that $X_{k}\geq\widehat{X}$. Since $X_{0}\geq\widehat{X}\geq
0$, we obtain that $X_{k}\geq\widehat{X}$ for each $k$.

Next, we show that the sequence $\\{X_{k}\\}$ is monotonically nonincreasing.
First, we show that $X_{0}\geq X_{1}$. Since $X_{1}$ is the stabilizing
solution of

$\displaystyle A_{0}^{\top}X+XA_{0}-XG_{0}X+H_{0}=0,$

we have

$\displaystyle(A_{0}-G_{0}X_{0})^{\top}(X_{1}-X_{0})+(X_{1}-X_{0})(A_{0}-G_{0}X_{0})=(X_{0}-X_{1})G_{0}(X_{0}-X_{1})-\mathcal{R}(X_{0})\geq
0.$

Since $(A_{0}-G_{0}X_{0})$ is stable, from Lemma 3.2, $X_{0}\geq X_{1}$.
Suppose that $X_{k-1}\geq X_{k}\geq 0$ and $X_{k+1}$ stabilizing solution of

$\displaystyle\begin{bmatrix}X&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}X\\\
-I\end{bmatrix}=0.$

Since $X_{k}$ is stabilizing solution satisfying (3.22), we have
$A_{k-1}-G_{k-1}X_{k}$ is stable. On the other hand,

$\displaystyle(A_{k-1}^{\top}-X_{k}G_{k-1})(X_{k+1}-X_{k})+(X_{k+1}-X_{k})(A_{k-1}-G_{k-1}X_{k})$
$\displaystyle=$
$\displaystyle(X_{k+1}-X_{k})G_{k-1}(X_{k+1}-X_{k})+\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k-1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}.$

From Lemma 3.1 and using the fact that $X_{k-1}\geq X_{k}$, we have

$\displaystyle\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k-1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}\geq\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}=0.$

Hence, $X_{k}\geq X_{k+1}$ by Lemma 3.2. Since $X_{0}\geq X_{1}$, we find that
$\\{X_{k}\\}_{k=0}^{\infty}$ is a monotonically nonincreasing sequence by
mathematical induction. Finally, we show that $\mathcal{R}(X_{k})\leq 0$ for
each $k$. Since $\\{X_{k}\\}_{k=0}^{\infty}$ is monotonically nonincreasing,
from Lemma 3.1 we have

$\displaystyle\mathcal{R}(X_{k})=\begin{bmatrix}X_{k}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}X_{k}\\\
-I\end{bmatrix}\leq\begin{bmatrix}X_{k}&-I\end{bmatrix}\Omega(X_{k-1})\begin{bmatrix}X_{k}\\\
-I\end{bmatrix}=0.$

(ii). Since the sequence $\\{X_{k}\\}_{k=0}^{\infty}$ is monotonically
nonincreasing and bounded below by $\widehat{X}$, we have
$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{+}$, where $\widehat{X}_{+}$ is a
solution of SCARE (1.5). (iii). Using the fact that $A_{k}-G_{k}X_{k+1}$ is
stable for each $k$, by continuity
$\sigma\left(A_{c}(\widehat{X}_{+})-G_{c}(\widehat{X}_{+})\widehat{X}_{+}\right)\subset\mathbb{C}_{-}\bigcup
i\mathbb{R}$. ∎

######  Remark 2.

From Theorems 3.3 and 3.4, we obtain that $\widehat{X}_{-}$ and
$\widehat{X}_{+}$ are minimal and maximal PSD solutions of SCARE (1.5),
respectively. In addition, we obtain that all eiganvalues of
$A_{c}(\widehat{X}_{-})-G_{c}(\widehat{X}_{-})\widehat{X}_{-}$ and
$A_{c}(\widehat{X}_{+})-G_{c}(\widehat{X}_{+})\widehat{X}_{+}$ are in
$\mathbb{C}_{-}\bigcup i\mathbb{R}$. If the PSD solution of SCARE (1.5) is
unique, we have $\widehat{X}_{-}=\widehat{X}_{+}$.

From Theorems 3.3 and 3.4, we proved that the sequence
$\\{X_{k}\\}_{k=1}^{\infty}$ generated by the FP-CARE_SDA converges
monotonically. Suppose that

$\displaystyle\lim_{k\to\infty}X_{k}=\widehat{X}\geq 0.$

From (2.4), we have the implicit relation between $X_{k}$ and $X_{k+1}$ as

$\displaystyle\mathcal{R}(X_{k},X_{k+1})$
$\displaystyle=A_{k}^{\top}X_{k+1}+X_{k+1}A_{k}-X_{k+1}G_{k}X_{k+1}+H_{k}=0,$

where $R_{k}=R+\Pi_{22}(X_{k})$, $G_{k}=BR_{k}^{-1}B^{\top}$,
$A_{k}=A-BR_{k}^{-1}(L+\Pi_{12}(X_{k}))^{\top}$ and
$H_{k}=Q+\Pi_{11}(X_{k})-(L+\Pi_{12}(X_{k}))R_{k}^{-1}(L+\Pi_{12}(X_{k}))^{\top}$.
Let $\mathcal{H}^{n}$ be the set of all Hermitian matrices and

$\displaystyle\mathcal{R}_{X_{k}}^{\prime}:\mathcal{H}^{n}\to\mathcal{H}^{n},\quad\mathcal{R}_{X_{k+1}}^{\prime}:\mathcal{H}^{n}\to\mathcal{H}^{n}$

be the Frèchet partial derivative with respect to $X_{k}$ and $X_{k+1}$,
respectively, i.e.,

$\displaystyle\mathcal{R}_{X_{k}}^{\prime}(Z)$ $\displaystyle=\lim_{t\to
0}\frac{1}{t}(\mathcal{R}(X_{k}+tZ)-\mathcal{R}(X_{k})),$ (3.26a)
$\displaystyle\mathcal{R}_{X_{k+1}}^{\prime}(Z)$ $\displaystyle=\lim_{t\to
0}\frac{1}{t}(\mathcal{R}(X_{k+1}+tZ)-\mathcal{R}(X_{k+1})).$ (3.26b)

From the implicit function theorem in Banach space, it holds that

$\displaystyle\mathcal{R}_{X_{k}}^{\prime}(Z)+\mathcal{R}_{X_{k+1}}^{\prime}(Z)\frac{dX_{k+1}}{dX_{k}}=0.$
(3.27)

From (3.27), we can derive the derivation of the explicit expression of
$X_{k+1}:=F(X_{k})$.

From (3.26b) follows that

$\displaystyle\mathcal{R}_{X_{k+1}}^{\prime}(Z)=(A_{k}^{\top}-X_{k+1}G_{k})Z+Z(A_{k}-G_{k}X_{k+1}).$
(3.28)

Furthermore, the Frèchet derivatives of $H_{c}(X)$, $G_{c}(X)$ and $A_{c}(X)$
at $X=X_{k}$ are

$\displaystyle H_{c}^{{}^{\prime}}(X_{k})(Z)=$
$\displaystyle\Pi_{11}(Z)-\Pi_{12}(Z)R_{c}(X_{k})^{-1}L_{c}(X_{k})^{\top}-L_{c}(X_{k})R_{c}(X_{k})^{-1}\Pi_{12}(Z)^{\top}$
$\displaystyle+L_{c}(X_{k})R_{c}(X_{k})^{-1}\Pi_{22}(Z)R_{c}(X_{k})^{-1}L_{c}(X_{k})^{\top}$
$\displaystyle=$
$\displaystyle\left[I,-L_{c}(X_{k})R_{c}(X_{k})^{-1}\right]\Pi(Z)\left[\begin{array}[]{c}I\\\
-R_{c}(X_{k})^{-1}L_{c}(X_{k})^{\top}\end{array}\right],$ (3.29c)
$\displaystyle G_{c}^{{}^{\prime}}(X_{k})(Z)=$ $\displaystyle-
BR_{c}(X_{k})^{-1}\Pi_{22}(Z)R_{c}(X_{k})^{-1}B^{\top},$ (3.29d)
$\displaystyle A_{c}^{{}^{\prime}}(X_{k})(Z)=$ $\displaystyle
BR_{c}(X_{k})^{-1}\Pi_{22}(Z)R_{c}(X_{k})^{-1}L_{c}(X_{k})^{\top}-BR_{c}(X_{k})^{-1}\Pi_{12}(Z)^{\top}.$
(3.29e)

From (3.26a) and (3.29), we have

$\displaystyle\mathcal{R}_{X_{k}}^{\prime}(Z)=\left(A_{c}^{{}^{\prime}}(X_{k})(Z)\right)^{\top}X_{k+1}+X_{k+1}A_{c}^{{}^{\prime}}(X_{k})(Z)-X_{k+1}G_{c}^{{}^{\prime}}(X_{k})(Z)X_{k+1}+H_{c}^{{}^{\prime}}(X_{k})(Z)$
(3.30)

Denote

$\displaystyle\widehat{A}=A_{c}(\widehat{X}),\ \
\widehat{G}=G_{c}(\widehat{X}),\ \ \widehat{L}=L_{c}(\widehat{X}),\ \
\widehat{R}=R_{c}(\widehat{X}).$

From (3.28), (3.29) and (3.30), we set

$\displaystyle\mathcal{L}_{\widehat{X}}(Z)\equiv$
$\displaystyle(\widehat{A}^{\top}-\widehat{X}\widehat{G})Z+Z(\widehat{A}-\widehat{G}\widehat{X}),$
$\displaystyle\Psi_{\widehat{X}}(Z)\equiv$
$\displaystyle\left(\widehat{L}\widehat{R}^{-1}\Pi_{22}(Z)\widehat{R}^{-1}B^{\top}-\Pi_{12}(Z)\widehat{R}^{-1}B^{\top}\right)\widehat{X}+\widehat{X}\left(B\widehat{R}^{-1}\Pi_{22}(Z)\widehat{R}^{-1}\widehat{L}^{\top}-B\widehat{R}^{-1}\Pi_{12}(Z)^{\top}\right)$
$\displaystyle+\widehat{X}B\widehat{R}^{-1}\Pi_{22}(Z)\widehat{R}^{-1}B^{\top}\widehat{X}+\left[I,-\widehat{L}\widehat{R}^{-1}\right]\Pi(Z)\left[\begin{array}[]{c}I\\\
-\widehat{R}^{-1}\widehat{L}^{\top}\end{array}\right].$

######  Theorem 3.5.

Suppose $\lim_{k\to\infty}X_{k}=\widehat{X}\geq 0$ monotonically by SDA-CARE.
If $\rho(\mathcal{L}_{\widehat{X}}^{-1}\Psi_{\widehat{X}})<1$, then
$\\{X_{k}\\}_{k=1}^{\infty}$ converges to $\widehat{X}$ R-linearly with

$\displaystyle\limsup_{k\to\infty}\sqrt[k]{\|X_{k}-\widehat{X}\|}\leq\widehat{\gamma}<1$

for some $0<\widehat{\gamma}<1$ and $\|\cdot\|$ being any matrix norm.

###### Proof.

From implicit function theorem, there is a function $F$ such that

$\displaystyle X_{k+1}=F(X_{k}),\quad\widehat{X}=F(\widehat{X}).$ (3.31)

Then

$\displaystyle
X_{k+1}-\widehat{X}=F(X_{k})-F(\widehat{X})=F^{\prime}(\widehat{X})(X_{k}-\widehat{X})+o(\|X_{k}-\widehat{X}\|).$

By implicit function Theorem, we have

$\displaystyle-F^{\prime}(\widehat{X})(Z)=\mathcal{L}_{\widehat{X}}^{-1}\Psi_{\widehat{X}}(Z).$
(3.32)

Since $\rho(\mathcal{L}_{\widehat{X}}^{-1}\Psi_{\widehat{X}})<1$, there exists
a matrix norm $\|\cdot\hat{\|}$ and a constant $M$ such that $\|\cdot\|\leq
M\|\cdot\hat{\|}$ and

$\displaystyle\|\mathcal{L}_{\widehat{X}}^{-1}\Psi_{\widehat{X}}\hat{\|}<1.$
(3.33)

Therefore, it implies

$\displaystyle\|X_{k+1}-\widehat{X}\|\leq\|F^{\prime}(\widehat{X})^{k}\|\|X_{0}-\widehat{X}\|\leq
M\|F^{\prime}(\widehat{X})^{k}\hat{\|}\|X_{0}-\widehat{X}\|.$

That is

$\displaystyle\sqrt[k]{\|X_{k+1}-\widehat{X}\|}\leq
M^{1/k}\|F^{\prime}(\widehat{X})\hat{\|}\|X_{0}-\widehat{X}\|^{1/k}.$

Hence,
$\limsup_{k\to\infty}\sqrt[k]{\|X_{k}-\widehat{X}\|}\leq\widehat{\gamma}<1$
for some $\widehat{\gamma}<1$. ∎

##  4 Newton method and Modified Newton method

Newton’s method for solving the SCARE (1.5) has been studied in [10]. The
associated Newton’s iteration is

$\displaystyle
X_{k+1}=X_{k}-(\mathcal{R}_{X_{k}}^{\prime})^{-1}\mathcal{R}(X_{k}),\ \mbox{
for }k=0,1,2,\ldots,$ (4.1)

provided that the Frèchet derivative $\mathcal{R}_{X_{k}}^{\prime}$ are all
invertible. The iteration in (4.1) is equivalent to solve $X_{k+1}$ of the
nonlinear matrix equation

$\displaystyle A_{X_{k}}^{\top}X+XA_{X_{k}}+\Pi_{X_{k}}(X)=-M_{X_{k}},$ (4.2a)
where $\displaystyle A_{X_{k}}$
$\displaystyle=A-BR_{k}^{-1}S_{X_{k}}^{\top}=A_{k}-G_{k}X_{k},$ (4.2b)
$\displaystyle\Pi_{X_{k}}(X)$
$\displaystyle=\begin{bmatrix}I&-S_{X_{k}}R_{k}^{-1}\end{bmatrix}\Pi(X)\begin{bmatrix}I\\\
-R_{k}^{-1}S_{X_{k}}^{\top}\end{bmatrix},$ (4.2c) $\displaystyle M_{X_{k}}$
$\displaystyle=\begin{bmatrix}I&-S_{X_{k}}R_{k}^{-1}\end{bmatrix}\begin{bmatrix}Q&L\\\
L^{\top}&R\end{bmatrix}\begin{bmatrix}I\\\
-R_{k}^{-1}S_{X_{k}}^{\top}\end{bmatrix}$ (4.2d) with $\displaystyle
S_{X_{k}}=X_{k}B+L_{c}(X_{k})\equiv X_{k}B+L_{k}.$ (4.2e)

The sequence $\\{X_{k}\\}_{k=0}^{\infty}$ computed by Newton’s method
converges to the maximal solution $\widehat{X}_{+}$ of (1.5) under mild
conditions.

######  Theorem 4.1 ([10]).

Assume that there exists a solution $\widehat{X}$ to $\mathcal{R}(X)\geq 0$
and a stabilizing matrix $X_{0}$ (i.e., $A_{X_{0}}$ is stable and
$\rho((A_{X_{0}}^{\top}X_{0}+X_{0}A_{X_{0}})^{-1}\Pi_{X_{0}})<1$). Then,

(i)

$X_{k}\geq X_{k+1}\geq\widehat{X}$, $\mathcal{R}(X_{k})\leq 0$ for all $k$.

(ii)

$\lim_{k\to\infty}X_{k}=\widehat{X}_{+}$ is the maximal solution of (1.5).

The authors in [10] provided the convergence analysis in Theorem 4.1. However,
it lacks the numerical methods for solving (4.2a). A method of solving (4.2a)
is the following fixed-point iteration

$\displaystyle
A_{X_{k}}^{\top}Y_{j+1}+Y_{j+1}A_{X_{k}}=-(\Pi_{X_{k}}(Y_{j})+M_{X_{k}})$
(4.3)

for $j=0,1,2,\ldots$ with $Y_{0}=X_{k}$.

The Lyapunov equation in (4.3) is a special case of CARE. As shown in [14], we
can find the Möbius transformation $g(z)=\frac{z+\alpha}{z-\alpha}$ with
$\alpha>0$ to transform (4.3) into a DARE

$\displaystyle Y=A_{d}^{\top}YA_{d}+H_{d}$ (4.4a) with $\displaystyle A_{d}$
$\displaystyle=I+2\alpha(A_{X_{k}}-\alpha I)^{-1},$ (4.4b) $\displaystyle
H_{d}$ $\displaystyle=2\alpha(A_{X_{k}}^{\top}-\alpha
I)^{-1}(\Pi_{X_{k}}(Y_{j})+M_{X_{k}})(A_{X_{k}}-\alpha I)^{-1}.$ (4.4c)

Therefore, the Lyapunov equation (4.3) can be solved efficiently by L-SDA in
Algorithm 2.

0: Coefficient matrices $C\geq 0$ and $E$ (stable), parameter $\alpha>0$, the
stopping tolerance $\varepsilon$.

0: The solution $Y$.

1: Compute $E_{0}=I+2\alpha(E-\alpha I)^{-1}$ and
$Y_{0}=2\alpha(E^{\top}-\alpha I)^{-1}C(E-\alpha I)^{-1}$.

2: Set $k=0$.

3: repeat

4: Compute $E_{k+1}=E_{k}^{2}$ and $Y_{k+1}=Y_{k}+E_{k}^{\top}Y_{k}E_{k}$.

5: Compute the normalized residual $\delta$ of Lyapunov equation.

6: Set $k=k+1$.

7: until $\delta<\varepsilon$.

8: Set $Y=Y_{k}$.

Algorithm 2 L-SDA for solving Lyapunov equation $E^{\top}Y+YE+C=0$ Algorithm 3
Newton fixed-point Lyapunov (NT-FP-Lyap_SDA) method for solving SCARE (1.5)

0: $A,Q\in\mathbb{R}^{n\times n}$, $B,L\in\mathbb{R}^{n\times m}$,
$R=R^{\top}\in\mathbb{R}^{m\times m}$, $A_{0}^{i}\in\mathbb{R}^{n\times n}$,
$B_{0}^{i}\in\mathbb{R}^{n\times m}$ for $i=1,\ldots,r$, initial
$X_{0}^{\top}=X_{0}$, and a tolerance $\varepsilon$.

0: Solution $X$.

1: Set $k=0$ and $\delta=\infty$.

2: while $\delta>\varepsilon$ do

3: Compute $S_{k}=X_{k}B+L+\Pi_{12}(X_{k})$, $R_{k}=\Pi_{22}(X_{k})+R$.

4: Compute $A_{X_{k}}=A-BR_{k}^{-1}S_{k}^{\top}$,
$P_{k}=[I,-S_{k}R_{k}^{-\top}]^{\top}$ and
$M_{k}=P_{k}^{\top}\begin{bmatrix}Q&L\\\ L^{\top}&R\end{bmatrix}P_{k}$.

5: Set $Y_{0}=X_{k}$, $j=0$ and $\delta_{f}=\infty$.

6: while $\delta_{f}>\varepsilon$ do

7: Compute $\Pi_{Y_{j}}=P_{k}^{\top}\Pi(Y_{j})P_{k}$ and
$C_{j}=\Pi_{Y_{j}}+M_{k}$.

8: Use L-SDA to solve the stabilizing solution $Y_{{j+1}}$ of the Lyapunov
equation  $\displaystyle A_{X_{k}}^{\top}Y+YA_{X_{k}}+C_{j}=0.$

9: Compute the normalized residual $\delta_{f}$ of (4.2a) with $X=Y_{{j+1}}$
and set $j=j+1$.

10: end while

11: Set $X_{k+1}=Y_{j}$ and compute the normalized residual $\delta$ of (1.5)
with $X=X_{k+1}$.

12: Set $k=k+1$.

13: end while

14: Set $X=X_{k}$.

Algorithm 3, called the NT-FP-Lyap_SDA algorithm, states the processes of
solving SCARE in (1.5) by combining Newton’s iteration in (4.2a) with fixed-
point iteration in (4.3). Steps 2-13 are the outer iterations of Newton’s
method. The inner iterations in Steps 6-10 are the fixed-point iteration to
solve the nonlinear matrix equation in (4.2a).

It is well-known that the fixed-point iteration in (4.3) has only a linear
convergence if it converges. Even if the convergence of Newton’s method is
quadratic, the overall iterations for solving (4.2a) by (4.3) can be large. To
reduce the iterations, as the proposed method in [12], we modify the Newton’s
iteration in (4.2a) as

$\displaystyle A_{X_{k}}^{\top}X+XA_{X_{k}}=-\Pi_{X_{k}}(X_{k})-M_{X_{k}}$
(4.5)

and solve it by L-SDA in Algorithm 2. Solving (4.5) is equival to solve (4.2a)
by using one iteration of fixed-point iteration in (4.3) with $Y_{0}=X_{k}$.
We state the modified Newton’s iteration in Algorithm 4 and call it as mNT-FP-
Lyap_SDA algorithm.

Algorithm 4 mNT-FP-Lyap_SDA method for solving SCARE (1.5)

0: $A,Q\in\mathbb{R}^{n\times n}$, $B,L\in\mathbb{R}^{n\times m}$,
$R=R^{\top}\in\mathbb{R}^{m\times m}$, $A_{0}^{i}\in\mathbb{R}^{n\times n}$,
$B_{0}^{i}\in\mathbb{R}^{n\times m}$ for $i=1,\ldots,r$, and a tolerance
$\varepsilon$.

0: Solution $X$.

1: Set $k=0$ and $\delta=\infty$.

2: Choose an initial $X_{0}^{\top}=X_{0}$.

3: while $\delta>\varepsilon$ do

4: Compute $S_{k}=X_{k}B+L+\Pi_{12}(X_{k})$, $R_{k}=\Pi_{22}(X_{k})+R$.

5: Compute $A_{X_{k}}=A-BR_{k}^{-1}S_{k}^{\top}$,
$P_{k}=[I,-S_{k}R_{k}^{-\top}]^{\top}$ and
$M_{k}=P_{k}^{\top}\begin{bmatrix}Q&L\\\ L^{\top}&R\end{bmatrix}P_{k}$.

6: Compute $C_{k}=P_{k}^{\top}\Pi(X_{k})P_{k}+M_{k}$.

7: Use L-SDA to solve the stabilizing solution $X_{{k+1}}$ of the Lyapunov
equation  $\displaystyle A_{X_{k}}^{\top}X+XA_{X_{k}}+C_{k}=0.$

8: Compute the normalized residual $\delta$ of (1.5) with $X=X_{k+1}$.

9: Set $k=k+1$.

10: end while

11: Set $X=X_{k}$.

In Theorem 4.1, a stabilizing matrix $X_{0}$ is needed to guarantee the
convergence of Newton’s method. However, the initial $X_{0}$ is generally not
easy to find. To tackle this drawback, we integrate FP-CARE_SDA (FPC) method
in Algorithm 1 into Algorithm 3. We take a few iterations of Algorithm 1, e.g.
$k=5$, to obtain $X_{k}$ and then use $X_{k}$ as an initial $X_{0}$ in
Algorithm 3. We summarize it in Algorithm 5, called the FPC-NT-FP-Lyap_SDA
method.

Algorithm 5 FPC-NT-FP-Lyap_SDA method for solving SCARE (1.5)

0: $A,Q\in\mathbb{R}^{n\times n}$, $B,L\in\mathbb{R}^{n\times m}$,
$R=R^{\top}\in\mathbb{R}^{m\times m}$, $A_{0}^{i}\in\mathbb{R}^{n\times n}$,
$B_{0}^{i}\in\mathbb{R}^{n\times m}$ for $i=1,\ldots,r$, and a tolerance
$\varepsilon$.

0: Solution $X$.

1: Set $k=0$ and $\delta=\infty$.

2: Choose an initial $X_{0}^{\top}=X_{0}$.

3: Use FP-CARE_SDA algorithm (Algorithm 1) with initial $X_{0}$ to compute an
approximated $X_{k}$ satisfied $\|X_{k}-X_{k-1}\|<0.01$.

4: Use NT-FP-Lyap_SDA algorithm (Algorithm 3) with initial matrix $X_{k}$ to
solve SCARE (1.5).

##  5 Convergence Analysis for Modified Newton Method

In this section, we study the convergence of the modified Newton’s method. For
the special case that $L=0$ and $B_{0}^{i}=0$, $i=1,\ldots,r$, Guo has shown
the following result in [12].

######  Theorem 5.1.

[12] Let $\widehat{X}\geq 0$ be the solution of SCARE (1.5) with $L=0$ and
$B_{0}^{i}=0$, $i=1,\ldots,r$. Assume $X_{0}\geq\widehat{X}$ such that
$A_{0}-G_{0}X_{0}$ is stable and $\mathcal{R}(X_{0})\leq 0$. Then
$\\{X_{k}\\}_{k=0}^{\infty}$ generated by mNT-FP-Lyap_SDA algorithm is
satisfied $\mathcal{R}(X_{k})\leq 0$, $X_{0}\geq X_{1}\geq\cdots\geq
X_{k}\geq\cdots\geq\widehat{X}$, and $\lim_{k\to\infty}X_{k}=\widehat{X}_{+}$,
where $\widehat{X}_{+}$ is the maximal solution of SCARE (1.5).

Now, we consider the general case, i.e., $B_{0}^{i}\neq 0$. Here, we assume
that the sequence $\\{X_{k}\\}_{k=0}^{\infty}$ exists and $X_{k}\geq 0$ for
all $k$. Under this posterior assumption, we have the monotonic convergence of
$\\{X_{k}\\}_{k=0}^{\infty}$.

######  Lemma 5.2.

Let $C^{\top}C=Q-LR^{-1}L^{\top}$. Assume that ${\rm
Ker}C\subseteq\left(\bigcap_{i=1}^{r}{\rm Ker}A_{0}^{i}\right)\bigcap{\rm
Ker}L$ and $(C,A)$ is detectable. Suppose that $\\{X_{k}\\}_{k=0}^{\infty}$ is
a sequence generated by mNT-FP-Lyap_SDA algorithm (Algorithm 4). Then
$A_{k}-G_{k}X_{k}$ is stable if $X_{k+1}\geq 0$.

###### Proof.

From Proposition 2.2, $(H_{k},A_{k})$ is detectable. Let
$H_{k}=C_{k}^{\top}C_{k}$. From (4.2b) and (5.1), we have

$\displaystyle(A_{k}-G_{k}X_{k})^{\top}X_{k+1}+X_{k+1}(A_{k}-G_{k}X_{k})=-(H_{k}+X_{k}G_{k}X_{k}).$

Assume that $A_{k}-G_{k}X_{k}$ is not stable. Then there exist $y\neq 0$ and
$\mbox{Re}(\lambda)\geq 0$ such that $(A_{k}-G_{k}X_{k})y=\lambda y$ which
implies that

$\displaystyle 0\leq
2\mbox{Re}(\lambda)y^{*}X_{k+1}y=-(y^{*}H_{k}y+y^{*}X_{k}G_{k}X_{k}y).$

Since $H_{k},G_{k}\geq 0$, we obtain the following $H_{k}y=0$ and
$G_{k}X_{k}y=0$. This means that

$\displaystyle\lambda y=(A_{k}-G_{k}X_{k})y=A_{k}y,\quad y\neq
0,\quad\mbox{Re}(\lambda)\geq 0,\quad\mbox{ and }\quad H_{k}y=0,$

which contradicts the detectability of $(H_{k},A_{k})$. ∎

From (2.1), (4.2c) and (4.2d), we have

$\displaystyle\Pi_{X_{k}}(X_{k})+M_{X_{k}}$
$\displaystyle=\begin{bmatrix}I&-S_{X_{k}}R_{c}(X_{k})^{-1}\end{bmatrix}\begin{bmatrix}Q_{c}(X_{k})&L_{c}(X_{k})\\\
L_{c}(X_{k})^{\top}&R_{c}(X_{k})\end{bmatrix}\begin{bmatrix}I\\\
-R_{c}(X_{k})^{-1}S_{X_{k}}^{\top}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}I&-L_{c}(X_{k})R_{c}(X_{k})^{-1}\end{bmatrix}\begin{bmatrix}Q_{c}(X_{k})&L_{c}(X_{k})\\\
L_{c}(X_{k})^{\top}&R_{c}(X_{k})\end{bmatrix}\begin{bmatrix}I\\\
-R_{c}(X_{k})^{-1}L_{c}(X_{k})^{\top}\end{bmatrix}$
$\displaystyle+\begin{bmatrix}I&-L_{c}(X_{k})R_{c}(X_{k})^{-1}\end{bmatrix}\begin{bmatrix}Q_{c}(X_{k})&L_{c}(X_{k})\\\
L_{c}(X_{k})^{\top}&R_{c}(X_{k})\end{bmatrix}\begin{bmatrix}0\\\
-R_{c}(X_{k})^{-1}B^{\top}X_{k}\end{bmatrix}$
$\displaystyle+\begin{bmatrix}0&-X_{k}BR_{c}(X_{k})^{-1}\end{bmatrix}\begin{bmatrix}Q_{c}(X_{k})&L_{c}(X_{k})\\\
L_{c}(X_{k})^{\top}&R_{c}(X_{k})\end{bmatrix}\begin{bmatrix}I\\\
-R_{c}(X_{k})^{-1}S_{X_{k}}^{\top}\end{bmatrix}$
$\displaystyle=H_{c}(X_{k})+X_{k}BR_{c}(X_{k})^{-1}B^{\top}X_{k}.$ (5.1)

Substituting (4.2b) and (5.1) into (4.5), the modified Newton’s iteration in
(4.5) is equivalent to solve the following equation

$\displaystyle A_{k}^{\top}X+XA_{k}-XG_{k}X+H_{k}+(X-X_{k})G_{k}(X-X_{k})=0.$
(5.2)

The matrix $H_{c}(X)$ in (2.2d) is PSD and satisfies $H_{c}(X)\geq H_{c}(Y)$
for $X\geq Y\geq 0$. Suppose that $\widehat{X}_{-}\geq 0$ in Theorem 3.3 is
the minimal PSD solution of SCARE (1.5). From (3.29c), the Frèchet derivative
$H_{c}^{{}^{\prime}}(\widehat{X}_{-}):\mathcal{H}^{n}\rightarrow\mathcal{H}^{n}$
of $H_{c}(X)$ at $\widehat{X}_{-}\geq 0$ is given by

$\displaystyle H_{c}^{{}^{\prime}}(\widehat{X}_{-})(Z)=$
$\displaystyle[I,-L_{c}R_{c}^{-1}]\Pi(Z)\left[\begin{array}[]{c}I\\\
-R_{c}^{-1}L_{c}^{\top}\end{array}\right],$ (5.5)

where $Z\in\mathcal{H}^{n}$ is Hermitian, $R_{c}=R_{c}(\widehat{X}_{-})$, and
$L_{c}=L_{c}(\widehat{X}_{-})$. Note that if $Z\geq 0$, then
$H_{c}^{{}^{\prime}}(\widehat{X}_{-})(Z)\geq 0$. In the following theorem, we
show the monotonically nondecreasing property if the initial $X_{0}$ is close
enough to the solution $\widehat{X}_{-}$.

######  Theorem 5.3.

Let $\widehat{X}_{-}\geq 0$ in Theorem 3.3 be the minimal PSD solution of
SCARE (1.5). Assume that there exists $\alpha>0$ such that the Frèchet
derivative of $H_{c}(X)$ at $X=\widehat{X}_{-}$ satisfies
$H_{c}^{{}^{\prime}}(\widehat{X}_{-})(Z)\geq\alpha Z$ for all $Z\geq 0$. Let
$X_{0}\leq\widehat{X}_{-}$ be close enough to the solution $\widehat{X}_{-}$
and $\mathcal{R}(X_{0})\geq 0$. Suppose that $\\{X_{k}\\}_{k=0}^{\infty}$ with
$X_{k}\geq 0$ for all $k$ is a sequence generated by mNT-FP-Lyap_SDA algorithm
(Algorithm 4). Then

(i)

$X_{k}\leq X_{k+1}$, $X_{k}\leq\widehat{X}_{-}$, $\mathcal{R}(X_{k})\geq 0$
for each $k\geq 0$.

(ii)

$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{-}$.

###### Proof.

Since $X_{k+1}\geq 0$, from Lemma 5.2, $A_{k}-G_{k}X_{k}$ is stable for
$k=0,1,\ldots$. We first show that

$\displaystyle\begin{bmatrix}Y&I\end{bmatrix}(\Omega(Y)-\Omega(X_{k}))\left[\begin{array}[]{c}Y\\\
I\end{array}\right]\geq(Y-X_{k})G_{k}(Y-X_{k}),$ (5.8)

for $X_{0}\leq X_{k}\leq Y\leq\widehat{X}_{-}$. From the first-order
expansion, we have

$\displaystyle
H_{c}(Y)-H_{c}(X_{k})=H_{c}^{{}^{\prime}}(X_{k})(Y-X_{k})+o(Y-X_{k}).$

Using the fact that $H_{c}^{{}^{\prime}}(\widehat{X}_{-})(Z)\geq\alpha Z$ for
$Z\geq 0$, $G_{c}(\widehat{X}_{-})\geq 0$ is bounded and $X_{0}$ is
sufficiently close to $\widehat{X}_{-}$, by continuity we have

$\displaystyle H_{c}(Y)-H_{c}(X_{k})\geq\frac{\alpha}{2}(Y-X_{k})\text{ and
}G_{k}=G_{c}(X_{k})\text{ is bounded.}$

Hence,

$\displaystyle H_{c}(Y)-H_{c}(X_{k})\geq(Y-X_{k})G_{k}(Y-X_{k}).$

From Lemma 3.1, we have $\Omega(Y)-\Omega(X_{k})\geq 0$, and then

$\displaystyle\begin{bmatrix}Y&I\end{bmatrix}(\Omega(Y)-\Omega(X_{k}))\left[\begin{array}[]{c}Y\\\
I\end{array}\right]\geq\begin{bmatrix}Y&I\end{bmatrix}\left[\begin{array}[]{cc}0&0\\\
0&H_{c}(Y)-H_{c}(X_{k})\end{array}\right]\left[\begin{array}[]{c}Y\\\
I\end{array}\right]\geq(Y-X_{k})G_{k}(Y-X_{k}).$

This shows that (5.8).

Now, we prove assertion (i) by induction. For $k=0$, we already have
$X_{0}\leq\widehat{X}_{-}$ and $\mathcal{R}(X_{0})\geq 0$. From (5.1), we have

$\displaystyle(A_{0}^{\top}-X_{0}G_{0})(X_{1}-X_{0})+(X_{1}-X_{0})(A_{0}-G_{0}X_{0})=-\mathcal{R}(X_{0})\leq
0.$

Since $A_{0}-G_{0}X_{0}$ is stable, by Lemma 3.2, we obtain that $X_{0}\leq
X_{1}$. We now assume that assertion (i) is true for $k\geq 0$. From (5.2), we
obtain that

$\displaystyle(A_{k}^{\top}-X_{k}G_{k})(\widehat{X}_{-}-X_{k+1})+(\widehat{X}_{-}-X_{k+1})(A_{k}-G_{k}X_{k})$
$\displaystyle=\begin{bmatrix}\widehat{X}_{-}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}\widehat{X}_{-}\\\
-I\end{bmatrix}+(\widehat{X}_{-}-X_{k})G_{k}(\widehat{X}_{-}-X_{k}).$

Using the fact that $\widehat{X}_{-}$ is a solution of SCARE, it follows from
(5.8) that

$\displaystyle\begin{bmatrix}\widehat{X}_{-}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}\widehat{X}_{-}\\\
-I\end{bmatrix}+(\widehat{X}_{-}-X_{k})G_{k}(\widehat{X}_{-}-X_{k})\leq 0.$

Therefore, $X_{k+1}\leq\widehat{X}_{-}$ by Lemma 3.2. Now, we show that
$X_{k+1}\leq X_{k+2}$ and $\mathcal{R}(X_{k+1})\geq 0$. Since

$\displaystyle[X_{k+1},-I]\Omega(X_{k})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}+(X_{k+1}-X_{k})G_{k}(X_{k+1}-X_{k})=0$

and $X_{k}\leq X_{k+1}\leq\widehat{X}_{-}$, we have

$\displaystyle\mathcal{R}(X_{k+1})$
$\displaystyle\equiv\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}$
$\displaystyle=[X_{k+1},-I](\Omega(X_{k+1})-\Omega(X_{k}))\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}-(X_{k+1}-X_{k})G_{k}(X_{k+1}-X_{k}).$

From (5.8), we have $\mathcal{R}(X_{k+1})\geq 0$. Then we obtain that

$\displaystyle(A_{k+1}^{\top}-X_{k+1}G_{k+1})(X_{k+2}-X_{k+1})+(X_{k+2}-X_{k+1})(A_{k+1}-G_{k+1}X_{k+1})$
$\displaystyle=-\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}=-\mathcal{R}(X_{k+1})\leq 0.$

Hence, $X_{k+1}\leq X_{k+2}.$ The induction process is complete.

(ii). From the assertion (i), we have the sequence
$\\{X_{k}\\}_{k=0}^{\infty}$ is monotonically nondecreasing and bounded above
by $\widehat{X}_{-}$. This implies that
$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{-}$. ∎

Next, we show the monotonically nonincreasing property and the convergence of
$\\{X_{k}\\}_{k=0}^{\infty}$ as follows.

######  Theorem 5.4.

Let $\widehat{X}\geq 0$ be the solution of SCARE (1.5). Let
$X_{0}\geq\widehat{X}$ such that $A_{0}-G_{0}X_{0}$ is stable and
$\mathcal{R}(X_{0})\leq 0$. Suppose that $\\{X_{k}\\}_{k=0}^{\infty}$ with
$X_{k}\geq 0$ for all $k$ is a sequence generated by mNT-FP-Lyap_SDA algorithm
(Algorithm 4). Then

(i)

$X_{k}\geq X_{k+1}$, $X_{k}\geq\widehat{X}$, $\mathcal{R}(X_{k})\leq 0$ for
each $k\geq 0$.

(ii)

$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{+}$, where
$\widehat{X}_{+}\geq\widehat{X}$ is a solution of SCARE (1.5).

(iii)

$\sigma\left(A_{c}(\widehat{X}_{+})-G_{c}(\widehat{X}_{+})\widehat{X}_{+}\right)\subset\mathbb{C}_{-}\bigcup
i\mathbb{R}$.

###### Proof.

(i). Since $X_{k+1}\geq 0$, from Lemma 5.2, $A_{k}-G_{k}X_{k}$ is stable for
$k=0,1,\ldots$. Now, we prove by induction that for each $k\geq 0$,

$\displaystyle X_{k}\geq X_{k+1},\ \ X_{k}\geq\widehat{X}\text{ and
}\mathcal{R}(X_{k})\leq 0.$ (5.9)

For $k=0$, we already have $X_{0}\geq\widehat{X}$ and $\mathcal{R}(X_{0})\leq
0$. From (5.1), we have

$\displaystyle(A_{0}^{\top}-X_{0}G_{0})(X_{1}-X_{0})+(X_{1}-X_{0})(A_{0}-G_{0}X_{0})=-\mathcal{R}(X_{0})\geq
0.$

Since $A_{0}-G_{0}X_{0}$ is stable, by Lemma 3.2, we obtain that $X_{0}\geq
X_{1}$. We now assume that (5.9) is true for $k\geq 0$. From (5.2) and Lemma
3.1, we obtain that

$\displaystyle(A_{k}^{\top}-X_{k}G_{k})(\widehat{X}-X_{k+1})+(\widehat{X}-X_{k+1})(A_{k}-G_{k}X_{k})=\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}+(\widehat{X}-X_{k})G_{k}(\widehat{X}-X_{k})$
$\displaystyle\geq\begin{bmatrix}\widehat{X}&-I\end{bmatrix}\Omega(\widehat{X})\begin{bmatrix}\widehat{X}\\\
-I\end{bmatrix}+(\widehat{X}-X_{k})G_{k}(\widehat{X}-X_{k})=(\widehat{X}-X_{k})G_{k}(\widehat{X}-X_{k})\geq
0.$

Therefore, $X_{k+1}\geq\widehat{X}$ by Lemma 3.2. Now, we show that
$X_{k+1}\geq X_{k+2}$ and $\mathcal{R}(X_{k+1})\leq 0$. From Lemma 3.1,
$X_{k}\geq X_{k+1}$ and $X_{k+1}$ satisfies (5.2), we have

$\displaystyle\mathcal{R}(X_{k+1})$ $\displaystyle\equiv
A_{k+1}^{\top}X_{k+1}+X_{k+1}A_{k+1}-X_{k+1}G_{k+1}X_{k+1}+H_{k+1}=\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}$
$\displaystyle\leq\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}=-(X_{k+1}-X_{k})G_{k}(X_{k+1}-X_{k})\leq 0.$

Then we obtain that

$\displaystyle(A_{k+1}^{\top}-X_{k+1}G_{k+1})(X_{k+2}-X_{k+1})+(X_{k+2}-X_{k+1})(A_{k+1}-G_{k+1}X_{k+1})$
$\displaystyle=\begin{bmatrix}X_{k+2}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+2}\\\
-I\end{bmatrix}-\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}+(X_{k+2}-X_{k+1})G_{k+1}(X_{k+2}-X_{k+1})$
$\displaystyle=-\begin{bmatrix}X_{k+1}&-I\end{bmatrix}\Omega(X_{k+1})\begin{bmatrix}X_{k+1}\\\
-I\end{bmatrix}=-\mathcal{R}(X_{k+1})\geq 0.$

Hence, $X_{k+1}\geq X_{k+2}$. The induction process is complete.

(ii). From the assertion (i), we have the sequence
$\\{X_{k}\\}_{k=0}^{\infty}$ is monotonically nonincreasing and bounded below
by $\widehat{X}$. This implies that
$\lim_{k\rightarrow\infty}X_{k}=\widehat{X}_{+}$, where $\widehat{X}_{+}$ is a
solution of SCARE (1.5). (iii). Since $A_{k}-G_{k}X_{k}$ is stable for each
$k$, by continuity
$\sigma\left(A_{c}(\widehat{X}_{+})-G_{c}(\widehat{X}_{+})\widehat{X}_{+}\right)\subset\mathbb{C}_{-}\bigcup
i\mathbb{R}$. ∎

In practice, solving Lyapunov equation by Algorithm 2 outperforms solving CARE
by SDA algorithm. However, it is difficult to find an initial matrix $X_{0}$
satisfying the conditions in Theorem 5.3 or Theorem 5.4. Moreover, comparing
the Lyapunov equation in (4.5) with the CARE in (2.4), the CARE is more
approximated the SCARE (1.5) than the Lyapunov equation. The convergence of
FP-CARE_SDA algorithm will be better than that of mNT-FP-Lyap_SDA algorithm as
shown in Figure 1.

To preserve the advantage and tackle the drawbacks of mNT-FP-Lyap_SDA
algorithm, we proposed a practical algorithm that integrates FP-CARE_SDA and
mNT-FP-Lyap_SDA algorithms. If we ignore the quadratic term
$(X-X_{k})G_{k}(X-X_{k})$ in (5.2), then (5.2) is equal the CARE in (2.4).
Assume that $X_{k}\geq 0$ in (5.2) is closed to the nonnegative solution of
the SCARE (1.5). The norm of this quadratic term will be small, which means
that FP-CARE_SDA and mNT-FP-Lyap_SDA algorithms will have the similar
convergence. Furthermore, from Theorem 3.3, $\\{X_{k}\\}_{k=0}^{\infty}$
produced by FP-CARE_SDA with $X_{0}=0$ is monotonically nondecreasing and
$\mathcal{R}(X_{k})\geq 0$ for all $k\geq 0$. Therefore, we can use FP-
CARE_SDA method in Algorithm 1 with $X_{0}=0$ to obtain $X_{k}$, which is
close enough to $\widehat{X}_{-}$ and $\mathcal{R}(X_{k})\geq 0$. Based on
Theorem 5.3, such $X_{k}$ can be used as an initial matrix of the mNT-FP-
Lyap_SDA algorithm to improve the convergence of the mNT-FP-Lyap_SDA
algorithm. We call this integrating method as the FPC-mNT-FP-Lyap_SDA method
and state in Algorithm 6.

Algorithm 6 FPC-mNT-FP-Lyap_SDA method for solving SCARE (1.5)

0: $A,Q\in\mathbb{R}^{n\times n}$, $B,L\in\mathbb{R}^{n\times m}$,
$R=R^{\top}\in\mathbb{R}^{m\times m}$, $A_{0}^{i}\in\mathbb{R}^{n\times n}$,
$B_{0}^{i}\in\mathbb{R}^{n\times m}$ for $i=1,\ldots,r$, and a tolerance
$\varepsilon$.

0: Solution $X$.

1: Set $k=0$ and $\delta=\infty$.

2: Choose an initial $X_{0}^{\top}=X_{0}$.

3: Use FP-CARE_SDA algorithm (Algorithm 1) with initial $X_{0}$ to compute an
approximated $X_{k}$ satisfied $\|X_{k}-X_{k-1}\|<0.01$.

4: Use mNT-FP-Lyap_SDA algorithm (Algorithm 4) with initial matrix $X_{k}$ to
solve SCARE (1.5).

##  6 Numerical experiments

In this section, we construct eight examples of SCAREs, in which $A$, $B$,
$Q$, $R$ come from the references and $A_{0}^{1},\ldots,A_{0}^{r}$,
$B_{0}^{1},\ldots,B_{0}^{r}$ are constructed in this paper. In the first four
examples, we find an initial $X_{0}$ satisfying the conditions in Theorems 3.4
and 5.4 to numerically validate $\\{X_{k}\\}_{k=0}^{\infty}$ produced by FP-
CARE_SDA and mNT-FP-Lyap_SDA algorithms to be monotonically nonincreasing and
compare the convergence of the FP-CARE_SDA, NT-FP-Lyap_SDA and mNT-FP-Lyap_SDA
algorithms. Such initial $X_{0}$ is difficult to find in Examples 5 \- 8, then
we use $X_{0}=0$ as the initial matrix. Theorem 3.3 shows
$\\{X_{k}\\}_{k=0}^{\infty}$ to be monotonically nondecreasing for the FP-
CARE_SDA algorithm. Numerical results show the failure of the convergence for
the NT-FP-Lyap_SDA and mNT-FP-Lyap_SDA algorithms. However, the FPC-NT-FP-
Lyap_SDA and FPC-mNT-FP-Lyap_SDA algorithms converge well to demonstrate the
robustness of FP-CARE_SDA algorithm.

We use the normalized residual

$\displaystyle\mbox{NRes}_{k}\equiv\frac{\|\mathcal{R}(X_{k})\|_{F}}{2\|AX_{k}\|_{F}+\|Q\|_{F}+\|\Pi_{11}(X_{k})\|_{F}+\|\mathcal{B}(X_{k})\|_{F}}$

to measure the precision of the calculated solutions $X_{k}$, where
$\mathcal{B}(X)=(XB+L+\Pi_{12}(X))(R+\Pi_{22}(X))^{-1}(XB+L+\Pi_{12}(X))^{\top}$.
In Subsections 6.1 and 6.2, we demonstrate the convergence behaviors of
$\mbox{NRes}_{k}$ for FP-CARE_SDA (Algorithm 1), NT-FP-Lyap_SDA (Algorithm 3),
mNT-FP-Lyap_SDA (Algorithm 4), FPC-NT-FP-Lyap_SDA (Algorithm 5), FPC-mNT-FP-
Lyap_SDA (Algorithm 6), and fixed-point iteration in (8.2) proposed by Guo and
Liang [13]. The efficiency of these algorithms is demonstrated in Subsection
6.3.

All computations in this section are performed in MATLAB 2023a with an Apple
M2 Pro CPU.

###  6.1 Numerical Validation of the convergence

######  Example 1.

Let

$\displaystyle A$ $\displaystyle=\begin{bmatrix}0.9512&0\\\
0&0.9048\end{bmatrix},\quad B=\begin{bmatrix}4.8770&4.8770\\\
-1.1895&3.5690\end{bmatrix},\quad Q=\begin{bmatrix}0.005&0\\\
0&0.020\end{bmatrix},\quad R=\begin{bmatrix}\frac{1}{3}&0\\\
0&3\end{bmatrix},\quad L=0,$ $\displaystyle A_{0}^{1}$
$\displaystyle=\begin{bmatrix}-0.1&0.1\\\ -0.2&0.2\end{bmatrix},\quad
A_{0}^{2}=\begin{bmatrix}1&-0.1\\\ 0.5&0\end{bmatrix},\quad
A_{0}^{3}=\begin{bmatrix}0&-0.2\\\ 0.2&0.5\end{bmatrix},$ $\displaystyle
B_{0}^{1}$ $\displaystyle=\begin{bmatrix}0&-0.1\\\ 0.1&0\end{bmatrix},\quad
B_{0}^{2}=\begin{bmatrix}0.5&1\\\ -0.1&0.2\end{bmatrix},\quad
B_{0}^{3}=\begin{bmatrix}1&-1\\\ -0.2&1\end{bmatrix},$

where $A$, $B$, $Q$ and $R$ are given in Example 2.2 of [2].

######  Example 2.

Let

$\displaystyle A$
$\displaystyle=\varepsilon\begin{bmatrix}\frac{7}{3}&\frac{2}{3}&0\\\
\frac{2}{3}&2&-\frac{2}{3}\\\ 0&-\frac{2}{3}&\frac{5}{3}\end{bmatrix},\quad
B=\frac{1}{\sqrt{\varepsilon}}I_{3},\quad
A_{0}^{1}=0.1\begin{bmatrix}0.1&-0.1&0.01\\\ -0.2&0.1&-0.1\\\
0.05&-0.01&0.3\end{bmatrix},\quad R=I_{3},\quad L=0,$ $\displaystyle Q$
$\displaystyle=\begin{bmatrix}(4\varepsilon+4+\varepsilon^{-1})/9&2(2\varepsilon-1-\varepsilon^{-1})/9&2(2-\varepsilon-\varepsilon^{-1})/9\\\
2(2\varepsilon-1-\varepsilon^{-1})/9&(1+4\varepsilon+4/\varepsilon)/9&2(-1-\varepsilon+2/\varepsilon)/9\\\
2(2-\varepsilon-\varepsilon^{-1})/9&2(-1-\varepsilon+2/\varepsilon)/9&(4+\varepsilon+4/\varepsilon)/9\end{bmatrix},\quad
B_{0}^{1}=0.1\begin{bmatrix}0&0&0.2\\\ 0.36&-0.6&0\\\
0&-0.95&-0.032\end{bmatrix}$

with $\varepsilon=0.01$, where $A$, $B$, $Q$ and $R$ are given in Example 12
of [7].

######  Example 3.

Let

$\displaystyle A$ $\displaystyle=\begin{bmatrix}0.9512&0\\\
0&0.9048\end{bmatrix},\quad B=\begin{bmatrix}4.8770&4.8770\\\
-1.1895&3.5690\end{bmatrix},\quad Q=\begin{bmatrix}0.0028&-0.0013\\\
-0.0013&0.0190\end{bmatrix},$ $\displaystyle R$
$\displaystyle=\begin{bmatrix}1/3&0\\\ 0&3\end{bmatrix},\quad
A_{0}^{1}=6.5\begin{bmatrix}0.1&0.2\\\ 0.2&0.1\end{bmatrix},\quad
B_{0}^{1}=6.5I_{2},\quad L=0,$

modified in [27].

######  Example 4.

Let

$\displaystyle A$ $\displaystyle=\begin{bmatrix}3-\varepsilon&1\\\
4&2-\varepsilon\end{bmatrix},\quad B=\begin{bmatrix}1\\\ 1\end{bmatrix},\quad
Q=\begin{bmatrix}4\varepsilon-11&2\varepsilon-5\\\
2\varepsilon-5&2\varepsilon-2\end{bmatrix},\quad
A_{0}^{1}=\begin{bmatrix}0.1&-0.1\\\ -0.2&0.1\end{bmatrix},$ $\displaystyle
B_{0}^{1}$ $\displaystyle=\begin{bmatrix}0.1\\\ 0\end{bmatrix},\quad
L=\begin{bmatrix}0\\\ 0\end{bmatrix},\quad R=1,\ \mbox{ with }\
\varepsilon=5,$

where $A$, $B$, $Q$ and $R$ are given Example 11 of [7].

(a) Example 1

(b) Example 2 with $\varepsilon=0.01$

(c) Example 3

(d) Example 4 with $\varepsilon=5$

Figure 1:  $\mbox{NRes}_{k}$ for FP-CARE_SDA, FPC-NT-FP-Lyap_SDA, FPC-mNT-FP-
Lyap_SDA and fixed-point method.

In Examples 1 \- 4, we firstly take $X_{0}=0$ as the initial matrix for FP-
CARE_SDA, NT-FP-Lyap_SDA and mNT-FP-Lyap_SDA algorithms. FP-CARE_SDA has a
monotonically nondecreasing convergence as shown in Theorem 3.3. The numerical
result shows that (i) the produced sequence converges to the solution
$\widehat{X}_{-}\geq 0$ of the SCARE; (ii) NT-FP-Lyap_SDA and mNT-FP-Lyap_SDA
algorithms fail to converge to $\widehat{X}_{-}$.

Next, we choose an initial $X_{0}$ with $\mathcal{R}(X_{0})\leq 0$,
$X_{0}\geq\widehat{X}$, $A_{0}-G_{0}X_{0}$ being stable and
$\rho(\mathcal{L}_{A_{X_{0}}}^{-1}\Pi_{X_{0}})<1$ so that
$\\{X_{k}\\}_{k=0}^{\infty}$ produced by the FP-CARE_SDA, NT-FP-Lyap SDA and
mNT-FP-Lyap SDA algorithms are monotonically nonincreasing and converge to the
solution $\widehat{X}_{+}$ of the SCARE. Using the same initial $X_{0}$ for
FP-CARE_SDA, NT-FP-Lyap_SDA and mNT-FP-Lyap_SDA algorithms, the normalized
residual $\mbox{NRes}_{k}$ of the algorithms are shown in Figure 1. Numerical
results show that $\widehat{X}_{+}=\widehat{X}{-}$ for these four examples.

The numerical results in Figure 1 show that the convergence of the FP-CARE_SDA
algorithm is obviously better than that of the mNT-FP-Lyap_SDA algorithm.
Furthermore, the normalized residuals at the first few iterations of the FP-
CARE_SDA algorithm are less than those of the NT-FP-Lyap_SDA and mNT-FP-
Lyap_SDA algorithms. These indicate that we can use a few iterations of FP-
CARE_SDA algorithm to produce a good initial matrix $X_{0}$ for the NT-FP-
Lyap_SDA and mNT-FP-Lyap_SDA algorithms, which are, called FPC-NT-FP-Lyap_SDA
and FPC-mNT-FP-Lyap_SDA, respectively, used to improve the convergence and
efficiency of the algorithms.

###  6.2 Robustness of FP-CARE_SDA algorithm

Due to the difficulty in selecting the initial matrix for Theorems 3.4 and
5.4, based on Theorem 3.3, we take $X_{0}=0$ as the initial matrix so that the
FP-CARE_SDA algorithm has monotonically nondecreasing convergence. Numerical
results show that the NT-FP-Lyap_SDA and mNT-FP-Lyap_SDA algorithms with
$X_{0}=0$ fail to converge. In the following, we introduce four typical
practical examples and discuss the convergence of the FPC-NT-FP-Lyap_SDA and
FPC-mNT-FP-Lyap_SDA algorithms for these examples to demonstrate the
robustness of FP-CARE_SDA algorithm.

######  Example 5.

The coefficient matrices, $A$, $B$, $Q$, $R$, and $L$ of a mathematical model
of position and velocity controls for a string of high-speed vehicles [1] can
be described as,

$\displaystyle A$ $\displaystyle=\begin{bmatrix}C&D\\\ &\ddots&\ddots\\\
&&C&D\\\ &&&C&-1\\\
&&&&-1\end{bmatrix}\in\mathbb{R}^{(2m-1)\times(2m-1)},\quad
C=\begin{bmatrix}-1&0\\\ 1&0\end{bmatrix},\quad D=\begin{bmatrix}0&0\\\
-1&0\end{bmatrix},$ $\displaystyle B$
$\displaystyle=[b_{ij}]\in\mathbb{R}^{(2m-1)\times m},\quad
b_{ij}=\begin{cases}1,&\mbox{ for }i=1,3,5,\ldots,2m-1,j=(i+1)/2,\\\ 0,&\mbox{
otherwise},\end{cases}$ $\displaystyle Q$
$\displaystyle=[q_{ij}]\in\mathbb{R}^{(2m-1)\times(2m-1)},\quad
q_{ij}=\begin{cases}10,&\mbox{ for }i=2,4,6,\ldots,2m-1,\ j=i,\\\ 0,&\mbox{
otherwise},\end{cases}$ $\displaystyle R$ $\displaystyle=I_{m},\quad
L=0_{(2m-1)\times m},$

where $m$ is the number of vehicles. Furthermore, the multiplicative white
noises are

$\displaystyle A_{0}^{i}$ $\displaystyle=0.1\times
i\times\frac{\|A\|_{\infty}}{\|\widehat{A}_{0}^{i}\|_{\infty}}\widehat{A}_{0}^{i}\text{
with }\widehat{A}_{0}^{i}=\mbox{{wgn}}(2m-1,2m-1,8\times i),\mbox{ for
}i=1,\ldots,5,$ $\displaystyle B_{0}^{i}$ $\displaystyle=0.15\times
i\times\frac{\|B\|_{\infty}}{\|\widehat{B}_{0}^{i}\|_{\infty}}\widehat{B}_{0}^{i}\text{
with }\widehat{B}_{0}^{i}=\mbox{{wgn}}(2m-1,m,3\times i),\mbox{ for
}i=1,\ldots,5,$

where wgn is a MATLAB built-in function used to generate white Gaussian noise.

######  Example 6.

The differential SDRE with impact angle guidance strategies models the 3D
missile/target interception engagement [20, 24]. The nonlinear dynamics can be
expressed as in (1.1), where the system, the control and the white noise
matrix are given by

$\displaystyle A(x)$ $\displaystyle=\begin{bmatrix}0&1&0&0&0\\\
0&\frac{-2\dot{r}}{r}&0&\frac{-1}{2}\dot{\psi}\sin(2x_{1}+2\theta_{f})&\frac{g_{1}}{z_{a}}\\\
0&0&0&1&0\\\
0&2\dot{\psi}\tan(x_{1}+\theta_{f})&0&\frac{-2\dot{r}}{r}&\frac{g_{2}}{z_{a}}\\\
0&0&0&0&-\eta\end{bmatrix}\ \ B(x)=\begin{bmatrix}0&0\\\
\frac{-\cos\theta_{M}}{r}&0\\\ 0&0\\\
\frac{\sin\theta_{M}\sin\psi_{M}}{r\cos\theta}&\frac{-\cos\psi_{M}}{r\cos\theta}\\\
0&0\end{bmatrix},$ $\displaystyle A_{0}^{i}$ $\displaystyle=0.2\times
i\times\frac{\|A\|_{\infty}}{\|\widehat{A}_{0}^{i}\|_{\infty}}\widehat{A}_{0}^{i}\
\text{ with }\ \widehat{A}_{0}^{i}=\mbox{{wgn}}(5,5,8\times i),\mbox{ for
}i=1,\ldots,4,$ (6.1a) $\displaystyle B_{0}^{i}$ $\displaystyle=0.1\times
i\times\frac{\|B\|_{\infty}}{\|\widehat{B}_{0}^{i}\|_{\infty}}\widehat{B}_{0}^{i}\
\text{ with }\ \widehat{B}_{0}^{i}=\mbox{{wgn}}(5,2,3\times i),\mbox{ for
}i=1,\ldots,4,$ (6.1b)

in which $r$ measures the distance between missile and target,
$\\{\psi,\theta\\}$ (respective to $\\{\psi_{M},\theta_{M}\\}$) are the
azimuth and elevation angles (respective to missile) corresponding to the
initial frame (respective the line-of-sight (LOS) frame). Furthermore, the
state vector $x=[x_{1},x_{2},x_{3},x_{4},x_{5}]^{\top}$ is defined by
$x_{1}=\theta-\theta_{f}$, $x_{2}=\dot{x}_{1}$, $x_{3}=\psi-\psi_{f}$,
$x_{4}=\dot{x}_{3}$ with $\theta_{f}$ and $\psi_{f}$ being prescribed final
angles, $x_{5}$ is a slow varying stable auxiliary variable governed by
$\dot{x}_{5}=-\eta x_{5}$, $\eta>0$, $g_{1}$ and $g_{2}$ are highly nonlinear
functions on $r$, $\theta$, $\psi$, $\theta_{f}$, $\psi_{f}$, $\theta_{T}$,
$\psi_{T}$, $a_{T}^{z}$ and $a_{T}^{y}$ (see [24] for details), where
$\\{\psi_{T},\theta_{T}\\}$ are the azimuth and elevation angles of the target
to the LOS frame, $a_{T}^{z}$ and $a_{T}^{y}$ are the lateral accelerations
for the target, and $(u_{1},u_{2})=(a_{M}^{z},a_{M}^{y})$ is the control
vector for the maneuverability of missile.

For a fixed state $x$ at the state-dependent technique of SDRE [20], the
coefficient matrices of the SCARE (1.3) associated with the optimization
problem (1.2) under the nonlinear dynamics for the missile/target engagement
have the forms.

$\displaystyle A$ $\displaystyle=\begin{bmatrix}0&1&0&0&0\\\
0&0.0696&0&-0.0307&-1.91\times 10^{-4}\\\ 0&0&0&1&0\\\
0&0.123&0&0.0696&6.13\times 10^{-4}\\\ 0&0&0&0&-0.1\end{bmatrix},\quad
B=\begin{bmatrix}0&0\\\ -9.13\times 10^{-5}&0\\\ 0&0\\\ 2.42\times
10^{-5}&-1.30\times 10^{-4}\\\ 0&0\end{bmatrix},$ $\displaystyle Q$
$\displaystyle={\rm diag}(1000,1000,1000,1000,0),\quad R=I_{2},\quad L=0,\quad
A_{0}^{i}\text{ and }B_{0}^{i}\text{ are as in \eqref{eq6.1}.}$

(a) Example 5 with $m=100$

(b) Example 6

(c) Example 7

(d) Example 8

Figure 2:  Normalized residual $\mbox{NRes}_{k}$ for FP-CARE_SDA, FPC-NT-FP-
Lyap_SDA, FPC-mNT-FP-Lyap_SDA and fixed-point method.

######  Example 7.

[6] Finite-time SDRE of F16 aircraft flight control system [6] can be
described as

$\displaystyle\dot{\mathbf{x}}\equiv\begin{bmatrix}\dot{u}\\\ \dot{v}\\\
\dot{w}\\\ \dot{p}\\\ \dot{q}\\\ \dot{r}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}(g\sin\theta)/u&0&0&0&-w&v\\\
(-g\sin\phi\cos\theta)/u&0&0&w&0&-u\\\ (-g\cos\phi\cos\theta)/u&0&0&-v&u&0\\\
0&0&0&c_{1}q/2&(c_{1}p+c_{2}r)/2&c_{2}q/2\\\
0&0&0&c_{3}p+c_{4}r/2&0&c_{4}p/2-c_{3}r\\\
0&0&0&c_{5}q/2&(c_{5}p+c_{6}r)/2&c_{6}q/2\end{bmatrix}\begin{bmatrix}u\\\ v\\\
w\\\ p\\\ q\\\ r\end{bmatrix}$ $\displaystyle+\begin{bmatrix}1/m&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&c_{lp}&0&c_{np}\\\ c_{mq}Z_{TP}&0&c_{mq}&0\\\
0&c_{lr}&0&c_{nr}\end{bmatrix}\begin{bmatrix}F_{T}\\\ L\\\ M\\\
N\end{bmatrix}\equiv A(\mathbf{x})\mathbf{x}+B(\mathbf{x})\mathbf{u},$

where $g$ is the gravity force, $m$ is the aircraft mass, $(u,v,w)$, $(p,q,r)$
are the aircraft velocity and angular velocity vectors, respectively, for roll
$\phi$, pitch $\theta$ and yaw $\psi$ angles, $c$’s parameters are the
suitable combinations of coefficients of the aircraft inertial matrix, the
control vector $\mathbf{u}=[F_{T},L,M,N]^{\top}$ consists of the thrust, the
resulted rolling, pitch and yawing moments from the layout of ailerons,
elevators and rudders, $Z_{TP}$ is the position of the thrust point.

For a fixed state $\mathbf{x}$, the coefficient matrices of the SCARE (1.3)
with the optimization problem (1.4) has the forms

$\displaystyle A$ $\displaystyle=\begin{bmatrix}3.958\times
10^{-5}&0&0&0&-5.866&-6.985\\\ 2.116\times 10^{-4}&0&0&5.866&0&-84.66\\\
-0.1158&0&0&6.985&84.66&0\\\ 0&0&0&1.791\times 10^{-4}&4.303\times
10^{-3}&-5.006\times 10^{-3}\\\ 0&0&0&-5.329\times 10^{-3}&0&-4.259\times
10^{-2}\\\ 0&0&0&-4.769\times 10^{-3}&3.253\times 10^{-2}&-1.791\times
10^{-4}\end{bmatrix},$ $\displaystyle B$
$\displaystyle=\begin{bmatrix}1.076\times 10^{-4}&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&7.780\times 10^{-5}&0&7.780\times 10^{-5}\\\ 3.964\times
10^{-6}&0&1.321\times 10^{-5}&0\\\ 0&1.211\times 10^{-6}&0&1.171\times
10^{-5}\end{bmatrix},$ $\displaystyle Q$ $\displaystyle=5000I_{6},\quad
R=2\times 10^{-4}I_{4},\quad L=0,$ $\displaystyle A_{0}^{i}$
$\displaystyle=0.012\times
i\times\frac{\|A\|_{\infty}}{\|\widehat{A}_{0}^{i}\|_{\infty}}\widehat{A}_{0}^{i}\
\text{ with }\ \widehat{A}_{0}^{i}=\mbox{{wgn}}(6,6,100\times i),\mbox{ for
}i=1,2,3,$ $\displaystyle B_{0}^{i}$ $\displaystyle=0.012\times
i\times\frac{\|B\|_{\infty}}{\|\widehat{B}_{0}^{i}\|_{\infty}}\widehat{B}_{0}^{i}\
\text{ with }\ \widehat{B}_{0}^{i}=\mbox{{wgn}}(6,4,40\times i),\mbox{ for
}i=1,2,3.$

######  Example 8.

[5] The SDRE optimal control design for quadrotors for enhancing robustness
against unmodeled disturbances [5] can be described as in (1.1) and (1.2) with
$\mathbf{x}=[u,v,w,p,q,r,\phi,\theta,z_{a}]^{\top}$,
$\mathbf{e}=\mathbf{x}-\mathbf{x}_{f}$,

$\displaystyle A(\mathbf{x})$
$\displaystyle=\begin{bmatrix}0&\frac{r}{2}&-\frac{q}{2}&0&-\frac{w}{2}&\frac{v}{2}&0&-\frac{g(\sin\theta)}{\theta}&0\\\
-\frac{r}{2}&0&\frac{p}{2}&\frac{w}{2}&0&-\frac{u}{2}&\frac{g(1+\cos\theta)\sin\phi}{2\phi}&-\frac{g(1-\cos\theta)\sin\phi}{2\theta}&0\\\
\frac{q}{2}&-\frac{p}{2}&0&-\frac{v}{2}&\frac{u}{2}&0&-\frac{2g}{\phi}\sin^{2}(\frac{\phi}{2})\cos^{2}(\frac{\theta}{2})&-\frac{2g}{\theta}\sin^{2}(\frac{\theta}{2})\cos^{2}(\frac{\phi}{2})&\frac{g}{z_{a}}\\\
0&0&0&0&\frac{c_{1}r}{2}&\frac{c_{1}q}{2}&0&0&0\\\
0&0&0&\frac{c_{2}r}{2}&0&\frac{c_{2}p}{2}&0&0&0\\\
0&0&0&\frac{c_{3}q}{2}&\frac{c_{3}p}{2}&0&0&0&0\\\
0&0&0&1&\frac{\sin\phi\tan\theta}{3}&\frac{\tan\theta(1+\cos\phi)}{3}&\frac{\alpha_{1}}{3\phi}&\frac{\alpha_{2}}{3\theta}&0\\\
0&0&0&0&\frac{1+\cos\phi}{2}&-\frac{\sin\phi}{2}&-\frac{\alpha_{3}}{2\phi}&0&0\\\
0&0&0&0&0&0&0&0&-\eta\end{bmatrix},$ $\displaystyle B$
$\displaystyle=\begin{bmatrix}0&0&-\frac{1}{m}&0&0&0&0&0&0\\\
0&0&0&\frac{1}{I_{x}}&0&0&0&0&0\\\ 0&0&0&0&\frac{1}{I_{y}}&0&0&0&0\\\
0&0&0&0&0&\frac{1}{I_{z}}&0&0&0\end{bmatrix}^{\top},$ $\displaystyle Q$
$\displaystyle=\mbox{\rm diag}(2000,2000,3000,10,10,100,0,0,0),\ R=I_{4},\
L=0,$ $\displaystyle A_{0}^{i}$ $\displaystyle=0.025\times
i\times\frac{\|A\|_{\infty}}{\|\widehat{A}_{0}^{i}\|_{\infty}}\widehat{A}_{0}^{i}\
\text{ with }\ \widehat{A}_{0}^{i}=\mbox{{wgn}}(9,9,10\times i),\mbox{ for
}i=1,2,3,$ $\displaystyle B_{0}^{i}$ $\displaystyle=0.01\times
i\times\frac{\|B\|_{\infty}}{\|\widehat{B}_{0}^{i}\|_{\infty}}\widehat{B}_{0}^{i}\
\text{ with }\ \widehat{B}_{0}^{i}=\mbox{{wgn}}(9,4,4\times i),\mbox{ for
}i=1,2,3,$

where $\alpha_{1}=q\tan\theta\sin\phi-r\tan\theta+r\tan\theta\cos\phi$,
$\alpha_{2}=q\tan\theta\sin\phi+r\tan\theta\cos\phi$, and
$\alpha_{3}=q(1-\cos\phi)+r\sin\phi$, $g$ is the gravity force, $m$ is the
quadrotor mass, $(u,v,w)$ and $(p,q,r)$ are the velocity and the angular
velocity on the body-fixed frame, respectively, for roll $\phi$, pitch
$\theta$ and yaw $\psi$ angles, $z_{a}$ is a slow varying stable auxiliary
variable governed by $\dot{z}_{a}=-\eta z_{a}$, $\eta>0$,
$c_{1}=\frac{I_{y}-I_{z}}{I_{x}}$, $c_{2}=\frac{I_{z}-I_{x}}{I_{y}}$,
$c_{3}=\frac{I_{x}-I_{y}}{I_{z}}$, and $I_{x}$, $I_{y}$, $I_{z}$, are initial
parameters.

For a fixed state $\mathbf{x}$ with $m=1$, $I_{x}=I_{y}=0.01466$, and
$I_{z}=0.02848$, the coefficient matrix $A\equiv A(\mathbf{x})$ has the form

$\displaystyle A$
$\displaystyle=\begin{bmatrix}0&-8.208\mbox{e-}4&-1.047\mbox{e-}2&0&-1.234\mbox{e-}4&1.178&0&-9.8000&0\\\
8.208\mbox{e-}4&0&-1.603\mbox{e-}3&1.234\mbox{e-}4&0&2.203\mbox{e-}2&9.800&-5.436\mbox{e-}4&0\\\
1.047\mbox{e-}2&1.603\mbox{e-}3&0&-1.178&-2.203\mbox{e-}2&0&0&0&9.820\mbox{e-}1\\\
0&0&0&0&7.738\mbox{e-}4&-9.871\mbox{e-}3&0&0&0\\\
0&0&0&-7.738\mbox{e-}4&0&-1.511\mbox{e-}3&0&0&0\\\ 0&0&0&0&0&0&0&0&0\\\
0&0&0&1.000&1.386\mbox{e-}8&2.499\mbox{e-}4&2.617\mbox{e-}6&-5.464\mbox{e-}4&0\\\
0&0&0&0&1.000&0&-9.650\mbox{e-}3&0&0\\\ 0&0&0&0&0&0&0&0&-0.100\end{bmatrix}.$

In Examples 5 \- 8, we use $X_{0}=0$ as the initial matrix so that
$\\{X_{k}\\}_{k=0}^{\infty}$ is monotonically non-decreasing for the FP-
CARE_SDA algorithm. The associated convergence of $\mbox{NRes}_{k}$ for these
examples is presented in Figure 2. Numerical results show that

  * •

It is well known that Newton’s method (NT-FP-Lyap_SDA algorithm) is highly
dependent on the initial choice $X_{0}$. When $X_{0}=0$ or $X_{0}\geq 0$ is
randomly constructed, NT-FP-Lyap_SDA does not converge to $\widehat{X}\geq 0$
for these four examples. However, in Figure 2, after a few $m$ iterations of
FP-CARE_SDA for obtaining $X_{m}$ with $\|X_{m}-X_{m-1}\|_{2}\leq 0.01$,
$X_{m}$ is closed to $\widehat{X}$ and Newton’s method with initial $X_{m}$
(i.e. FPC-NT-FP-Lyap_SDA algorithm) has quadratic convergence.

  * •

Numerical results in Figure 2 show that the fixed-point method, FP-CARE_SDA
and FPC-mNT-FP-Lyap_SDA algorithms have linear convergence. But, the iteration
number of the fixed-point method is obviously larger than that of the other
methods. Comparing the results in Figure 2 with those in Figure 1, the
iteration number of mNT-FP-Lyap_SDA algorithm in Figure 1 is obviously larger
than that of FP-CARE_SDA algorithm. As shown in (5.2), the quadratic term
$(X-X_{k})G_{k}(X-X_{k})$ is small when $X_{m}$ is used as an initial matrix
of mNT-FP-Lyap_SDA algorithm. We can expect that the mNT-FP-Lyap_SDA algorithm
with initial $X_{m}$ (i.e. FPC-mNT-FP-Lyap_SDA algorithm) has the same
convergent behavior of the FP-CARE_SDA algorithm as presented in Figure 2.
Without such initial $X_{m}$ (i.e., with $X_{0}=0$ or a randomly initial
$X_{0}$), mNT-FP-Lyap_SDA algorithm will not converge for Examples 5 \- 8.
This demonstrates the importance and robustness of the FP-CARE_SDA algorithm.

  * •

For Examples 6 \- 8, FPC-mNT-FP-Lyap_SDA algorithm produces monotonically non-
decreasing sequences and has a linear convergence. The sequence produced by
FPC-NT-FP-Lyap_SDA algorithm does not contain any monotonic property.

From the results in Figures 1 and 2, the FP-CARE_SDA algorithm has a linear
convergence. Moreover, the FP-CARE_SDA algorithm also provides a good initial
matrix for the NT-FP-Lyap_SDA algorithm and the modified Newton’s iteration in
(4.5) so that the FPC-NT-FP-Lyap_SDA and FPC-mNT-FP-Lyap_SDA algorithms have
quadratic and linear convergence, respectively. These two algorithms can not
compute the solution $\widehat{X}$ without such an initial matrix. Therefore,
our proposed FP-CARE_SDA algorithm is a reliable and robust algorithm.

###  6.3 Efficiency of algorithms

In each iteration of the FP-CARE_SDA (Algorithm 1), FPC-NT-FP-Lyap_SDA
(Algorithm 5), and FPC-mNT-FP-Lyap_SDA (Algorithm 6), solving the CARE and
Lyapunov equation are the most computationally intensive steps. To compare the
efficiency of these methods, we first show the total numbers of solving CARE
(2.4), Lyapunov equation (4.3)/(4.5), and fixed-point iteration (8.2) for each
method in Table 1. The third and fifth columns show that we only need a few
iterations of the FP-CARE_SDA algorithm to get a good initial matrix in FPC-
NT-FP-Lyap_SDA and FPC-mNT-FP-Lyap_SDA algorithms.

| FP-CARE_SDA | FPC-NT-FP-Lyap_SDA | FPC-mNT-FP-Lyap_SDA | FP
---|---|---|---|---
| solving (2.4) | solving (2.4) | solving (4.3) | solving (2.4) | solving (4.5) | solving (8.2)
Exa. 5 | 14 | 3 | 16 | 3 | 11 | 86
Exa. 6 | 83 | 10 | 111 | 10 | 74 | 185
Exa. 7 | 74 | 8 | 114 | 8 | 71 | 410
Exa. 8 | 85 | 10 | 107 | 10 | 74 | 300
Table 1: Number of solving CARE (2.4), Lyapunov equation (4.3)/(4.5), and fixed-point (FP) iteration (8.2). | CPU times (seconds)
---|---
| FP-CARE_SDA | FPC-NT-FP-Lyap_SDA | FPC-mNT-FP-Lyap_SDA | FP
Exa. 5 | 7.947e-1 | 8.112e-1 | 5.729e-1 | 5.790
Exa. 6 | 6.151e-3 | 4.632e-3 | 3.928e-3 | 9.634e-3
Exa. 7 | 5.527e-3 | 4.897e-3 | 3.816e-3 | 2.041e-2
Exa. 8 | 1.217e-2 | 7.759e-3 | 6.627e-3 | 3.804e-2
Table 2: CPU times for solving SCARE (1.5).

  * •

The fourth and sixth columns of Table 1 show that the FPC-NT-FP-Lyap_SDA
algorithm still has a quadratic convergence, the convergence of the fixed-
point iteration (4.3) for the FPC-NT-FP-Lyap_SDA algorithm is linear so that
the total number of solving Lyapunov equations is larger than that for the
FPC-mNT-FP-Lyap_SDA algorithm. As shown in Table 2, the performance of the
FPC-mNT-FP-Lyap_SDA algorithm is better than that of the FPC-NT-FP-Lyap_SDA
algorithm.

  * •

Sum of the numbers in the fifth and sixth columns of Table 1 is almost equal
to the number in the second column. However, the computational cost for
solving CARE (2.4) by SDA in [14] is larger than for solving Lyapunov equation
(4.5) by L-SDA in Algorithm 2. This indicates that the FPC-mNT-FP-Lyap_SDA
algorithm outperforms the FP-CARE_SDA algorithm as shown in the comparison of
the CPU time in Table 2.

  * •

The dimension of the coefficient matrix for the linear system of the fixed-
point iteration in (8.2) is $n(r+1)$. Compared to the other three algorithms,
in which matrix dimensions are $n$, the fixed-point iteration needs more
computational cost to solve the linear systems. On the other hand, as shown in
Figures 1-2 and the last column in Table 1, the fixed-point iteration needs a
large iteration number to get the solution. These lead to the performance of
the fixed-point iteration being worse than that of the other three algorithms
as shown in Table 2.

Summarized above results, we can see that the FPC-mNT-FP-Lyap_SDA algorithm
outperforms the FP-CARE_SDA, FPC-NT-FP-Lyap_SDA algorithms, and fixed-point
iteration.

##  7 Conclusions

We propose a reliable FP-CARE_SDA for solving the SCARE (1.3) and prove its
monotonically non-decreasing convergence with $X_{0}=0$, its monotonically
non-increasing convergence with $\mathcal{R}(X_{0})\leq 0$ and
$A_{0}-G_{0}X_{0}$ being stable. Furthermore, we propose the mNT-FP-Lyap_SDA
algorithm to be used to accelerate the convergence with the FP-CARE_SDA as a
robust initial step and prove its monotonically non-decreasing convergence
under the assumption that the resulting $X_{k}\geq 0$. Numerical experiments
of real-world practical applications from the 3D missile/target engagement,
the F16 aircraft flight control, and the quadrotor optimal control show the
robustness of our proposed FPC-mNT-FP-Lyap_SDA algorithm.

##  8 Appendix

In [13], the authors proposed a fixed-point iteration to solve stochastic
discrete-time algebraic Riccati equations (SDAREs) and proved the convergence
of the fixed-point iteration. Moreover, the Möbius transformation is applied
to transform the SCARE (1.5) into a SDARE so that the proposed fixed-point
iteration can be applied to solve the associated solution $X$. Now, we state
the fixed-point iteration for solving SCARE. Let

$\displaystyle\widehat{A}=A-BR^{-1}L^{\top},\quad\widehat{B}=BR^{-1/2},\quad\widehat{C}^{\top}\widehat{C}=Q-LR^{-1}L^{\top},$

$\Pi$ and $\widehat{\Pi}\in\mathbb{R}^{n(r+1)\times n(r+1)}$ be the
permutations satisfying

$\displaystyle\Pi^{\top}(X\otimes I_{r})\Pi$ $\displaystyle=I_{r}\otimes X,$
$\displaystyle\mbox{diag}\left(X,X\otimes I_{r}\right)$
$\displaystyle=\widehat{\Pi}^{\top}\left(X\otimes
I_{r+1}\right)\widehat{\Pi},$

and

$\displaystyle\mathcal{A}=\Pi\left(\begin{bmatrix}A_{0}^{1}\\\ \vdots\\\
A_{0}^{r}\end{bmatrix}-\begin{bmatrix}B_{0}^{1}\\\ \vdots\\\
B_{0}^{r}\end{bmatrix}R^{-1}L^{\top}\right),\quad\mathcal{B}=\Pi\begin{bmatrix}B_{0}^{1}\\\
\vdots\\\ B_{0}^{r}\end{bmatrix}R^{-/12}.$

Give $\gamma>0$ so that $\widehat{A}_{\gamma}\equiv\widehat{A}-\gamma I_{n}$
is nonsingular. Define
$Z_{\gamma}=\widehat{C}\widehat{A}_{\gamma}^{-1}\widehat{B}$. By Möbius
transformation, we have

$\displaystyle E_{\gamma}$
$\displaystyle=\widehat{\Pi}\begin{bmatrix}\widehat{A}_{\gamma}+2\gamma
I_{n}+\widehat{B}Z_{\gamma}^{\top}\widehat{C}\\\
\sqrt{2\gamma}(\mathcal{A}+\mathcal{B}Z_{\gamma}^{\top}\widehat{C})\end{bmatrix}(I_{n}+\widehat{A}_{\gamma}^{-1}\widehat{B}Z_{\gamma}^{\top}\widehat{C})^{-1}\widehat{A}_{\gamma}^{-1}\in\mathbb{R}^{n(r+1)\times
n},$ (8.1a) $\displaystyle H_{\gamma}$
$\displaystyle=2\gamma\widehat{A}_{\gamma}^{-\top}\widehat{C}^{\top}(I_{\ell}+Z_{\gamma}Z_{\gamma}^{\top})^{-1}\widehat{C}\widehat{A}_{\gamma}^{-1}\in\mathbb{R}^{n\times
n},$ (8.1b) $\displaystyle G_{\gamma}$
$\displaystyle=\widehat{\Pi}\begin{bmatrix}\sqrt{2\gamma}\widehat{A}_{\gamma}^{-1}\widehat{B}\\\
\mathcal{A}\widehat{A}_{\gamma}^{-1}\widehat{B}-\mathcal{B}\end{bmatrix}(I_{m}+Z_{\gamma}^{\top}Z_{\gamma})^{-1}\begin{bmatrix}\sqrt{2\gamma}\widehat{A}_{\gamma}^{-1}\widehat{B}\\\
\mathcal{A}\widehat{A}_{\gamma}^{-1}\widehat{B}-\mathcal{B}\end{bmatrix}^{\top}\widehat{\Pi}^{\top}\in\mathbb{R}^{n(r+1)\times
n(r+1)}.$ (8.1c)

Then, the fixed-point iteration in [13] is given as

$\displaystyle X_{1}$ $\displaystyle=H_{\gamma},$ $\displaystyle X_{k+1}$
$\displaystyle=E_{\gamma}^{\top}(X_{k}\otimes
I_{r+1})\left(I_{n(r+1)}+G_{\gamma}(X_{k}\otimes
I_{r+1})\right)^{-1}E_{\gamma}+H_{\gamma}$ (8.2)

for $k=1,2,\ldots$.

Each iteration in (8.2) only needs to compute $(X_{k}\otimes
I_{r+1})\left(I_{n(r+1)}+G_{\gamma}(X_{k}\otimes
I_{r+1})\right)^{-1}E_{\gamma}$ by the direct method. However, the dimension
of the coefficient matrix is enlarged to $n(r+1)\times n(r+1)$. It needs to
have more and more computational cost as $r$ is large. On the other hand, the
iteration can not preserve the symmetric property of the solution.

## Acknowledgments

This work was partially supported by the National Science and Technology
Council (NSTC), the National Center for Theoretical Sciences, and the Nanjing
Center for Applied Mathematics. T.-M. Huang, Y.-C. Kuo, and W.-W. Lin were
partially supported by NSTC 110-2115-M-003-012-MY3, 110-2115-M-390-002-MY3 and
112-2115-M-A49-010-, respectively.

## References

  * [1] J. Abels and P. Benner.  CAREX - a collection of benchmark examples for continuous-time algebraic Riccati equations.  SLICOT Working Note 1999-16 (Version 2.0), December 1999.  Available online at slicot.org/working-notes/.
  * [2] J. Abels and P. Benner.  DAREX – a collection of benchmark examples for discrete-time algebraic Riccati equations.  SLICOT Working Note 1999-14 (Version 2.0), December 1999.  Available online at slicot.org/working-notes/.
  * [3] R. Babazadeh and R. Selmic.  Distance-based multi-agent formation control with energy constraints using SDRE.  IEEE Trans. Aerosp. Electron. Syst., 56(1):41–56, 2020.
  * [4] P. Benner and J. Heiland.  Exponential stability and stabilization of extended linearizations via continuous updates of Riccati-based feedback.  Int. J. Robust Nonlinear Control, 28(4):1218–1232, 2018.
  * [5] S.-W. Cheng and H.-A. Hung.  Robust state-feedback ${H}_{\infty}$ control of quadrotor.  International Automatic Control Conference, 2022.
  * [6] M. Chodnicki, P. Pietruszewski, M. Weso ̵lowski, and S. Stȩpień.  Finite-time SDRE control of F16 aircraft dynamics.  Arch. Control Sci., 32:557–576, 2022.
  * [7] E. K.-W. Chu, H.-Y. Fan, and W.-W. Lin.  A structure-preserving doubling algorithm for continuous-time algebraic Riccati equations.  Linear Algebra Appl., 396:55–80, 2005.
  * [8] E. K.-W. Chu, T. Li, W.-W. Lin, and C.-Y. Weng.  A modified Newton’s method for rational Riccati equations arising in stochastic control.  In 2011 International Conference on Communications, Computing and Control Applications (CCCA), pages 1–6, 2011.
  * [9] T. Çimen.  Survey of state-dependent Riccati equation in nonlinear optimal feedback control synthesis.  J. Guid. Control Dyn., 35(4):1025–1047, 2012.
  * [10] T. Damm and D. Hinrichsen.  Newton’s method for a rational matrix equation occurring in stochastic control.  Linear Algebra Appl., 332-334:81–109, 2001.
  * [11] V. Dragan, T. Morozan, and A.-M. Stoica.  Mathematical Methods in Robust Control of Linear Stochastic Systems.  Springer-Verlag, 2013.
  * [12] C.-H. Guo.  Iterative solution of a matrix Riccati equation arising in stochastic control.  In Linear Operators and Matrices, pages 209–221. Birkhäuser Basel, 2002.
  * [13] Z.-C. Guo and X. Liang.  Stochastic algebraic Riccati equations are almost as easy as deterministic ones theoretically.  arXiv, 2207.11220v3, 2023.
  * [14] T.-M. Huang, R.-C. Li, and W.-W. Lin.  Structure-Preserving Doubling Algorithms for Nonlinear Matrix Equations.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2018\.
  * [15] I. G. Ivanov.  Iterations for solving a rational Riccati equation arising in stochastic control.  Comput. Math. Appl., 53:977–988, 2007.
  * [16] S. Kim and S. J. Kwon.  Nonlinear optimal control design for underactuated two-wheeled inverted pendulum mobile platform.  22(6):2803–2808, 2017.
  * [17] P. Lancaster and L. Rodman.  Existence and uniqueness theorems for the algebraic Riccati equation.  Int. J. Control, 32:285–309, 1980.
  * [18] P. Lancaster and L. Rodman.  Algebraic Riccati Equations.  Oxford Science Publications. Oxford University Press, 1995.
  * [19] L.-G. Lin, Y.-W. Liang, and L.-J. Cheng.  Control for a class of second-order systems via a state-dependent Riccati equation approach.  SIAM J. Control Optim., 56(1):1–18, 2018.
  * [20] L.-G. Lin, R.-S. Wu, P.-K. Huang, M. Xin, C.-T. Wu, and W.-W. Lin.  Fast SDDRE-based maneuvering-target interception at prespecified orientation.  pages 1–8, 2023.
  * [21] L.-G. Lin and M. Xin.  Alternative SDRE scheme for planar systems.  IEEE Trans. Circuits Syst. II-Express Briefs, 66(6):998–1002, 2019\.
  * [22] W.-W. Lin and S.-F. Xu.  Convergence analysis of structure-preserving doubling algorithms for Riccati-type matrix equations.  SIAM J. Matrix Anal. Appl., 28(1):26–39, 2006.
  * [23] Y. Meng, Q. Chen, X. Chu, and A. Rahmani.  Maneuver guidance and formation maintenance for control of leaderless space-robot teams.  IEEE Trans. Aerosp. Electron. Syst., 55(1):289–302, 2019.
  * [24] R. V. Nanavati, S. R. Kumar, and A. Maity.  Spatial nonlinear guidance strategies for target interception at pre-specified orientation.  Aerosp. Sci. Technol., 114:106735, 2021.
  * [25] S. R. Nekoo.  Tutorial and review on the state-dependent Riccati equation.  8(2):109–166, 2019.
  * [26] V. Polito.  Nonlinear business cycle and optimal policy: A VSTAR perspective.  CESifo Working Paper Series, (8060):1–54, 2020.
  * [27] P. C.-Y. Weng and F. K.-H. Phoa.  Perturbation analysis and condition numbers of rational Riccati equations.  Ann. Math. Sci. Appl., 6:25–49, 2021.