Algorithms for LQR via Static Output Feedback for Discrete-Time LTI Systems

Abstract

Randomized and deterministic algorithms for the problem of LQR optimal control via static-output-feedback (SOF) for discrete-time systems are suggested in this chapter. The randomized algorithm is based on a recently introduced randomized optimization method named the Ray-Shooting Method that efficiently solves the global minimization problem of continuous functions over compact non-convex unconnected regions. The randomized algorithm presented here has a proof of convergence in probability to the global optimum. The suggested deterministic algorithm is based on the gradient method and thus can be proved to converge to local optimum only. A comparison between the algorithms is provided as well as the performance of the hybrid algorithm.

1. Introduction

The application of static-output-feedbacks (SOFs) for linear-quadratic regulators (LQR) is very attractive, since they are cheap and reliable and their implementation is simple and direct, because their components has direct physical interpretation in terms of sensors amplification rates and actuator activation power. Moreover, the long-term memory of dynamic feedbacks is useless for systems subject to random disturbances, to fast dynamic loadings or to random bursts and impulses, and the application of state feedbacks is not always possible due to unavailability of full-state measurements (see, e.g., [1]). Also, the use of SOF avoids the need to reconstruct the state by Kalman filter or by any other state reconstructor.

On the other hand, in practical applications, the entries of the needed SOFs are bounded, and since the problem of SOFs with interval constrained entries is NP-hard (see [2, 3]), one cannot expect the existence of a deterministic efficient (i.e., polynomial-time) algorithm to solve the problem. Randomized algorithms are thus natural solutions to the problem. The probabilistic and randomized methods for the constrained SOF problem and robust stabilization via SOFs (among other hard problems) are discussed in [4, 5, 6, 7]. For a survey of the SOF problem see [8], and for a recent survey of the robust SOF problem see [9].

The Ray-Shooting Method was recently introduced in [10], where it was used to derive the Ray-Shooting (RS) randomized algorithm for the minimal-gain SOF problem, with regional pole assignment, where the region can be non-convex and unconnected. The Ray-Shooting Method was successfully applied recently also to the following hard complexity control problems for continuous-time systems:

• Structured and structured-sparse SOFs (see [11])

• LQR via SOF for continuous-time LTI systems (see [12])

• LQR optimal, $H_{\infty }$ -optimal, and $H_2$ -optimal proportional-integral-differential (PID) controllers (see [13])

• Robust control via SOF (see [14])

The contribution of the research presented in the current chapter is as follows:

1. The randomized algorithm presented here (which we call the RS algorithm) is based on the Ray-Shooting Method (see [10]), which opposed to smooth optimization methods, and has the potential of finding a global optimum of continuous functions over compact non-convex and unconnected regions.

2. The RS algorithm has a proof of convergence in probability and explicit complexity.

3. Experience with the algorithm shows good quality of controllers (in terms of reduction of the LQR functional value with relatively small controller norms), high percent of success, and good run-time, for real-life systems. Thus, the suggested practical algorithm efficiently solves the problem of LQR via SOF for discrete-time systems.

4. The RS algorithm does not need to solve any Riccati or quadratic matrix equations (QMES) and thus can be applied to large systems.

5. The RS algorithm is one of the few, dealing with the problem of LQR via SOF for discrete-time systems.

6. A deterministic algorithm for the problem that generalizes the algorithm of Moerder and Calise [15], for discrete-time systems, is given (we call it the MC algorithm). The MC algorithm has a proof of convergence to a local optimum only, and it needs other algorithms for computing initial stabilizing SOF.

7. A comparison between the RS and the MC algorithms, as well as the performance of the hybrid algorithm, for real-life systems, is provided.

The reminder of the chapter is organized as follows:

In Section 2 we formulate the problem and give some useful lemmas (without a proof). In Section 3, we introduce the randomized algorithm for the problem of LQR via SOF for discrete-time LTI systems. Section 4 is devoted to the deterministic algorithm for the problem. In Section 5, we give the results of the algorithms for some real-life systems. Finally, in Section 6 we conclude with some remarks.

2. Preliminaries

Let a discrete-time system be given by

$$\left\{\begin{array}{l}{x}_{k+1}={\mathit{Ax}}_k+{\mathit{Bu}}_k,k=0,1,\dots \\ y_k={\mathit{Cx}}_k\end{array}\right.$$ E1

where $A\in {\mathbb{R}}^{p\times p},B\in {\mathbb{R}}^{p\times q}$ , $C\in {\mathbb{R}}^{r\times p}$ , and $x_0\in {\mathbb{R}}^p$ . Let the LQR cost functional be defined by

$$J\left(x_0,u\right)≔\sum _{k=0}^{\infty }\left(x_k^T{\mathit{Qx}}_k+u_k^T{\mathit{Ru}}_k\right),$$ E2

where $Q>0$ and $R>0$ . Let $u_k=-{\mathit{Ky}}_k$ be the SOF, and let $A_{c\mathrm{\ell }}\left(K\right)≔A-\mathit{BKC}$ denote the closed-loop matrix. Let $\mathbb{D}$ denote the open unit disk, let $0<\alpha <1$ , and let ${\mathbb{D}}_{\alpha }$ denote the set of all $z\in \mathbb{D}$ with $\left|z\right|<1-\alpha$ (where $\left|z\right|$ is the absolute value of $z$ ). For a square matrix $Z$ , we denote by $\sigma \left(Z\right)$ its spectrum. For any rectangular matrix $Z$ , we denote by $Z^{+}$ its Moore-Penrose pseudo-inverse. By ${‖Z‖}_F=\mathit{\text{trace}}{\left(Z^{T}Z\right)}^{\frac{1}{2}}$ we denote the Frobenius norm of $Z$ , and by $‖Z‖={\left(max\left(\sigma \left(Z^{T}Z\right)\right)\right)}^{\frac{1}{2}}$ we denote the spectral norm of $Z$ . By $L_Z$ and $R_Z$ , we denote the (left and right) orthogonal projections $I-Z^{+}Z$ and $I-{\mathit{ZZ}}^{+}$ on the spaces $\mathit{Ker}\left(Z\right)$ and $\mathit{Ker}\left(Z^{+}\right)$ , respectively. For a topological space $\mathcal{X}$ and a subset $\mathcal{U}\subset \mathcal{X}$ , we denote by $\mathit{int}\left(\mathcal{U}\right)$ the interior of $\mathcal{U}$ , i.e., the largest open set included in $\mathcal{U}$ . By $\overline{\mathcal{U}}$ we denote the closure of $\mathcal{U}$ , i.e., the smallest closed set containing $\mathcal{U}$ , and by $\partial \mathcal{U}=\overline{\mathcal{U}}-\mathit{int}\left(\mathcal{U}\right)$ we denote the boundary of $\mathcal{U}$ . Let ${\mathcal{S}}^{q\times r}$ denote the set of all matrices $K\in {\mathbb{R}}^{q\times r}$ such that $\sigma \left(A_{c\mathrm{\ell }}\right)\subset \mathbb{D}$ (i.e., stable in the discrete-time sense), and let ${\mathcal{S}}^{q\times r}$ denote the set of all matrices $K\in {\mathbb{R}}^{q\times r}$ such that $\sigma \left(A_{c\mathrm{\ell }}\right)\subset {\mathbb{D}}_{\alpha }$ . If the last is nonempty, we say that $A_{c\mathrm{\ell }}$ is $\alpha$ -stable and we call $\alpha$ the degree of stability. Let $K\in {\mathcal{S}}_{\alpha }^{q\times r}$ be given. Substitution of the SOF $u_k=-{\mathit{Ky}}_k=-{\mathit{KCx}}_k$ into (2) yields:

$$J\left(x_0,K\right)=\sum _{k=0}^{\infty }{x}_k^T\left(Q+C^T{K}^T\mathit{RKC}\right)x_k.$$ E3

Since $Q+C^T{K}^T\mathit{RKC}>0$ and $A_{c\mathrm{\ell }}\left(K\right)$ is stable, it follows that the Stein equation

$$P-A_{c\mathrm{\ell }}{\left(K\right)}^T{\mathit{PA}}_{c\mathrm{\ell }}\left(K\right)=Q+C^T{K}^T\mathit{RKC}$$ E4

has a unique solution $P>0$ , given by

$$P\left(K\right)=\mathit{mat}\left({\left(I_p\otimes I_p-A_{c\mathrm{\ell }}{\left(K\right)}^T\otimes A_{c\mathrm{\ell }}{\left(K\right)}^T\right)}^{-1}\cdot \mathit{vec}\left(Q+C^T{K}^T\mathit{RKC}\right)\right).$$ E5

Substitution of (4) into (3) and using $x_k=A_{c\mathrm{\ell }}{\left(K\right)}^k{x}_0$ with the fact that $A_{c\mathrm{\ell }}\left(K\right)$ is stable leads to

$$\begin{array}{ll}J\left(x_0,K\right)& =\sum _{k=0}^{\infty }{x}_k^T\left(P-A_{c\mathrm{\ell }}{\left(K\right)}^T{\mathit{PA}}_{c\mathrm{\ell }}\left(K\right)\right)x_k\\ & =\sum _{k=0}^{\infty }{x}_0^T{A}_{c\mathrm{\ell }}{\left(K\right)}^{\mathit{Tk}}\left(P-A_{c\mathrm{\ell }}{\left(K\right)}^T{\mathit{PA}}_{c\mathrm{\ell }}\left(K\right)\right)A_{c\mathrm{\ell }}{\left(K\right)}^k{x}_0\\ & =x_0^{T}P\left(K\right)x_0={‖P{\left(K\right)}^{\frac{1}{2}}{x}_0‖}^2.\end{array}$$

Thus, when $x_0$ is known, we search for $K\in {\mathcal{S}}_{\alpha }^{q\times r}$ that minimizes the functional

$$J\left(x_0,K\right)=x_0^{T}P\left(K\right)x_0.$$ E6

Let

$${\sigma }_{\mathit{max}}\left(K\right)≔max\left(\sigma \left(P\left(K\right)\right)\right).$$ E7

Now, since $\frac{J\left(x_0,K\right)}{{‖x_0‖}^2}\le {\sigma }_{\mathit{max}}\left(K\right)$ for any $x_0\ne 0$ , with equality in the worst case, therefore

$$\underset{x_0\ne 0}{sup}\left(\frac{J\left(x_0,K\right)}{{‖x_0‖}^2}\right)=\underset{x_0\ne 0}{sup}\left(\frac{{‖P{\left(K\right)}^{\frac{1}{2}}{x}_0‖}^2}{{‖x_0‖}^2}\right)={‖P{\left(K\right)}^{\frac{1}{2}}‖}^2=‖P\left(K\right)‖={\sigma }_{\mathit{max}}\left(K\right).$$

Thus, when $x_0$ is unknown, we search for $K\in {\mathcal{S}}_{\alpha }^{q\times r}$ , such that ${\sigma }_{\mathit{max}}\left(K\right)=‖P\left(K\right)‖$ is minimal. Note that if $\lambda$ is an eigenvalue of $A_{c\mathrm{\ell }}\left(K\right)$ and $v$ is a corresponding eigenvector, then (4) yields $1-{\left|\lambda \right|}^2=\frac{v^{\ast }\left(Q+C^T{K}^T\mathit{RKC}\right)v}{v^{\ast }P\left(K\right)v}\ge \frac{v^{\ast }\mathit{Qv}}{v^{\ast }P\left(K\right)v}>0$ . Therefore, ${\left|\lambda \right|}^2\le 1-\frac{v^{\ast }\mathit{Qv}}{v^{\ast }P\left(K\right)v}<1$ , and thus, minimizing ${\sigma }_{\mathit{max}}\left(K\right)$ results in eigenvalues that are getting closer to the boundary of $\mathbb{D}$ . Since $\alpha$ , the degree of stability, is important to get satisfactory decay rate of the state to $0$ , and for disturbance rejection, we allow the user of the algorithms to determine $\alpha$ . Note that too high value of $\alpha$ might result in nonexistence of any SOF for the system or in complicating the search for a starting SOF. Higher values of $\alpha$ result in higher values of the optimal value of the LQR functional, i.e., higher energy consumption for decaying the disturbance $x_0$ to $0$ .

The functionals $J\left(x_0,K\right)$ and ${\sigma }_{\mathit{max}}\left(K\right)$ are generally not convex since their domain of definition ${\mathcal{S}}_{\alpha }^{q\times r}$ (and therefore ${\mathcal{S}}_{\alpha }^{q\times r}$ ) is generally non-convex. Necessary conditions for optimality for continuous-time systems were given as three QMEs in [15, 16, 17, 18]. Necessary and sufficient conditions for optimality for continuous-time systems, based on linear matrix inequalities (LMI), were given in [19, 20, 21]. However, algorithms based on these formulations are generally not guaranteed to converge, seemingly because of the non-convexity of the coupled matrix equations or inequalities, and when they converge, it is to a local optimum only.

In the sequel, we will use the following lemmas, given here without proofs.

Lemma 2.1. We have:

1. The equation $\mathit{AX}=B$ has solutions if and only if ${\mathit{AA}}^{+}B=B$ or equivalently, if and only if $R_{A}B=0$ . In this case, the set of all solutions is given by

$$X=A^{+}B+L_{A}Z,$$

where $Z$ is arbitrary.

1. The equation $\mathit{XA}=B$ has solutions if and only if ${\mathit{BA}}^{+}A=B$ or equivalently, if and only if ${\mathit{BL}}_A=0$ . In this case, the set of all solutions is given by

$$X={\mathit{BA}}^{+}+{\mathit{ZR}}_A,$$

where $Z$ is arbitrary.

Lemma 2.2. We have:

1. Let $A,B,X$ be matrices with sizes $p\times q,r\times p,q\times r$ , respectively. Then

$$\frac{\partial }{\partial X}\mathit{\text{trace}}\left(\mathit{AXB}\right)=A^T{B}^T.$$

1. Let $A,B,C,X$ be matrices with sizes $p\times q,r\times r,q\times p,r\times q$ , respectively. Then

$$\frac{\partial }{\partial X}\mathit{\text{trace}}\left({\mathit{AX}}^T\mathit{BXC}\right)=B^T{\mathit{XA}}^T{C}^T+\mathit{\text{BXCA}}.$$

3. The randomized Ray-Shooting Method-based algorithm

The Ray-Shooting Method works as follows, for general function minimization: let $f\left(x\right)\ge 0$ be a continuous function defined over some compact set $X\subset {\mathbb{R}}^n$ . Let $\mathrm{ϵ}>0$ be given and assume that we want to compute $x_{\ast }\in X$ such that $y_{\ast }≔f\left(x_{\ast }\right)={min}_{x\in X}f\left(x\right)$ up to $\mathrm{ϵ}$ , i.e., to find $\left(x,y\right)$ in the set $\mathcal{S}\left(\mathrm{ϵ}\right)=\left\{\left(x,y\right)\left|x\in X,f\left(x_{\ast }\right)\le y=f\left(x\right)\le f\left(x_{\ast }\right)+\mathrm{ϵ}\right.\right\}$ . Let $x_0\in \mathcal{X}$ be given, let $y_0≔f\left(x_0\right)$ and let ${\mathcal{S}}^{\left(0\right)}=\left\{\left(x,y\right)\left|x\in X,f\left(x\right)\le y\le y_0\right.\right\}$ denote the search space, which is a subset of the epigraph of $f$ . Let ${\mathcal{D}}^{\left(0\right)}=\left\{\left(x,y\right)\left|x\in \mathcal{X},0\le y\le y_0\right.\right\}$ denote the cylinder enclosed between $\mathcal{X}$ and the level $y_0$ . Let ${\mathcal{L}}^{\left(0\right)}=\left\{\left(x,y\right)\left|x\in \mathcal{X},y=0\text{or}x\in \partial \mathcal{X},0\le y\le y_0\right.\right\}$ . Let $z_0≔\left(x_0,y_0\right)$ and note that $z_0\in {\mathcal{S}}^{\left(0\right)}$ . Then, we choose $w_0$ in ${\mathcal{L}}^{\left(0\right)}$ randomly, according to some distribution, and we define the ray as $z\left(t\right)≔\left(1-t\right)z_0+{\mathit{tw}}_0,0\le t\le 1$ . We scan the ray and choose the largest $0\le t_0\le 1$ such that $\left(1-t_0\right)z_0+t_0{w}_0\in {\mathcal{S}}^{\left(0\right)}$ (actually, we scan the ray from $t=1$ in equal-spaced points and take the first $t$ for which this happens). We define $z_1≔\left(1-t_0\right)z_0+t_0{w}_0$ and update sets ${\mathcal{S}}^{\left(0\right)}$ , ${\mathcal{D}}^{\left(0\right)}$ , and ${\mathcal{L}}^{\left(0\right)}$ by replacing $y_0$ with $y_1$ , where $\left(x_1,y_1\right)=z_1$ . Let ${\mathcal{S}}^{\left(1\right)}$ , ${\mathcal{D}}^{\left(1\right)}$ , and ${\mathcal{L}}^{\left(1\right)}$ denote the updated sets. We continue the process similarly from $z_1\in {\mathcal{S}}^{\left(1\right)}$ , and we define a sequence $z_n\in {\mathcal{S}}^{\left(n\right)},n=0,1,\dots$ . Note that $\mathcal{S}\left(\mathrm{ϵ}\right)\subset {\mathcal{S}}^{\left(n+1\right)}\subset {\mathcal{S}}^{\left(n\right)}$ for any $n=0,1,\dots$ , unless we have $z_n\in \mathcal{S}\left(\mathrm{ϵ}\right)$ for some $n$ (in which the process is ceased). One can show that the sequence ${\left\{z_n\right\}}_{n=0}^{\infty }$ converges (in probability) to a point in $\mathcal{S}\left(\epsilon \right)$ . Note that shooting rays from the points of local minimum have positive probability to hit $\mathcal{S}\left(\epsilon \right)$ (under the following mild assumption), because any global minimum is visible from any local minimum. Moreover, for a given level of certainty, we hit $\mathcal{S}\left(\epsilon \right)$ in a finite number of iterations (see Remark 3.1 below). Practically, we may stop the algorithm if no improvement is detected within a window of $20\%$ of the allowed number of iterations. The function need not be smooth or even continuous. It only needs to be well defined and measurable over the compact domain $\mathcal{X}$ , and $\mathcal{S}\left(\epsilon \right)$ should have non-negligible measure (i.e., should have some positive volume). Obviously, global minimum points belong to the boundary of the search space ${\mathcal{S}}^{\left(0\right)}$ , and actually such points are where the distance between the compact sets $\mathcal{X}\times \left\{0\right\}$ and ${\mathcal{S}}^{\left(0\right)}$ in ${\mathbb{R}}^{n+1}$ is accepted. This is essential for the efficiency of the Ray-Shooting Method, although we raised the search space dimension from $n$ to $n+1$ .

In order to apply the Ray-Shooting Method for the LQR via SOF problem, we need the following definitions: assume that $K^{\left(0\right)}\in \mathit{int}\left({\mathcal{S}}_{\alpha }\right)$ was found by the RS algorithm (see [10]) or by any other method (see [22, 23, 24]). Let $h>0$ and let ${\mathcal{U}}^{\left(0\right)}$ be a unit vector (actually a matrix, but we consider here the space of matrices as a normed vector space) with respect to the Frobenius norm, i.e., ${‖{\mathcal{U}}^{\left(0\right)}‖}_F=1$ . Let $L^{\left(0\right)}=K^{\left(0\right)}+h\cdot {\mathcal{U}}^{\left(0\right)}$ and let $\mathcal{L}$ be the hyperplane defined by $L^{\left(0\right)}+V$ , where ${⟨V,{\mathcal{U}}^{\left(0\right)}⟩}_F=0$ . Here $\mathcal{L}$ is the tangent space at $L^{\left(0\right)}$ to the closed ball $\overline{\mathbb{B}\left(K^{\left(0\right)},h\right)}$ centered at $K^{\left(0\right)}$ with radius $h$ , with respect to the Frobenius norm on ${\mathbb{R}}^{q\times r}$ . Let $r_{\infty }>0$ and let ${\mathcal{R}}_{\infty }$ denote the set of all $F\in \mathcal{L}$ , such that ${‖F-L^{\left(0\right)}‖}_F\le r_{\infty }$ . Let ${\mathcal{R}}_{\infty }\left(\mathrm{ϵ}\right)=R_{\infty }+\overline{\mathbb{B}\left(0,\mathrm{ϵ}\right)}$ , where $\overline{\mathbb{B}\left(0,\mathrm{ϵ}\right)}$ denotes the closed ball centered at $0$ with radius $\mathrm{ϵ}$ ( $0<\mathrm{ϵ}\le \frac{1}{2}$ ). Let ${\mathcal{D}}^{\left(0\right)}=\mathit{\text{chull}}\left(K^{\left(0\right)},{\mathcal{R}}_{\infty }\left(\mathrm{ϵ}\right)\right)$ denote the convex hull of the vertex $K^{\left(0\right)}$ with the basis ${\mathcal{R}}_{\infty }\left(\mathrm{ϵ}\right)$ . Let ${\mathcal{S}}_{\alpha }^{\left(0\right)}={\mathcal{S}}_{\alpha }\cap {\mathcal{D}}^{\left(0\right)}$ and note that ${\mathcal{S}}_{\alpha }^{\left(0\right)}$ is compact (but generally not convex). We wish to minimize the continuous function ${\sigma }_{\mathit{max}}\left(K\right)$ (or the continuous function $J\left(x_0,K\right)$ , when $x_0$ is known) over the compact set ${\mathcal{S}}_{\alpha }\cap \overline{\mathbb{B}\left(K^{\left(0\right)},h\right)}$ . Let $K_{\ast }$ denote a point in ${\mathcal{S}}_{\alpha }\cap \overline{\mathbb{B}\left(K^{\left(0\right)},h\right)}$ where the minimum of ${\sigma }_{\mathit{max}}\left(K\right)$ is accepted. Obviously, $K_{\ast }\in {\mathcal{D}}^{\left(0\right)}$ , for some direction $U^{\left(0\right)}$ from $K^{\left(0\right)}$ .

The Ray-Shooting Algorithm 1 for the LQR via SOF problem, works as follows: we start with a point $K^{\left(0\right)}\in \mathit{int}\left({\mathcal{S}}_{\alpha }\right)$ , found by the RS algorithm (see [10]). Assuming that $K_{\ast }\in {\mathcal{D}}^{\left(0\right)}$ , the inner loop ( $j=1,\dots ,n$ ) uses the Ray-Shooting Method in order to find an approximation of the global minimum of the function ${\sigma }_{\mathit{max}}\left(K\right)$ over ${\mathcal{S}}_{\alpha }^{\left(0\right)}$ —the portion of ${\mathcal{S}}_{\alpha }$ bounded in the cone $D^{\left(0\right)}$ . The proof of convergence in probability of the inner loop and its complexity (under the above-mentioned assumption) can be found in [10] (see also [11]). In the inner loop, we choose a search direction by choosing a point $F$ in ${\mathcal{R}}_{\infty }\left(\mathrm{ϵ}\right)$ —the base of the cone ${\mathcal{D}}^{\left(0\right)}$ . Next, in the most inner loop ( $k=0,\dots ,s$ ), we scan the ray $K\left(t\right)≔\left(1-t\right)K^{\left(0\right)}+\mathit{tF}$ and record the best controller on it. Repeating this sufficiently many times, we reach $K_{\ast }$ (or an $\mathrm{ϵ}$ neighborhood of it) with high probability, under the assumption that $K_{\ast }\in {\mathcal{D}}^{\left(0\right)}$ (see Remark 3.1).

The reasoning of the Ray-Shooting Method is that sampling the whole search space will lead to the probabilistic method that is doomed to the “curse of dimensionality,” which the method tries to avoid. This is achieved by slicing the search space into covering cones ( $m$ is the number of cones allowed), because any point in the cone is visible from its vertex. At each cone we shoot rays ( $n$ is the number of rays per cone) from its node toward its basis, where each ray is sampled from its head toward its tail, while updating the best point found so far. Note that the global minimum of ${\sigma }_{\max }\left(K\right)$ over any compact subset of ${\mathcal{S}}_{\alpha }$ is achieved on the boundary of the related portion of the epigraph of ${\sigma }_{\max }\left(K\right)$ . Therefore, we can break the most inner loop; in the first moment, we find an improvement in ${\sigma }_{\max }\left(K\right)$ . This bypasses the need to sample the whole search space (although we raise by $1$ the search space dimension) and explains the efficiency of the Ray-Shooting Method in finding global optimum. Another advantage of the Ray-Shooting Method which is specific to the problem of LQR via SOF is that the search is concentrated to the parameter space (the $\mathit{qr}$ -dimension space where the $K$ rests) and not to the certificate space (the $p^2$ -dimension space where the Lyapunov matrices $P$ rests). Thus, the method avoids the need to solve any Riccati, LMI, and BMI equations, which might make crucial difference for large-scale systems (i.e., where $p^2>>\mathit{qr}$ ).

Algorithm 1. The Ray-Shooting Algorithm for LQR via SOF for discrete-time systems.

Require: An algorithm for deciding $\alpha$ -stability, an algorithm for computing ${\sigma }_{\max }\left(K\right)$ and algorithms for general linear algebra operations.

Input: $0<\mathrm{ϵ}\le \frac{1}{2},0<\alpha <1,h>0,r_{\infty }>0$ , integers: $m,n,s>0$ , controllable pairs $\left(A,B\right)$ and $\left(A^T,C^T\right)$ , matrices $Q>0,R>0$ and $K^{\left(0\right)}\in \mathit{int}\left({\mathcal{S}}_{\alpha }\right)$ .

Output: $K\in {\mathcal{S}}_{\alpha }$ close as possible to $K_{\ast }$ .

1. compute $P\left(K^{\left(0\right)}\right)$ as in (5)

2. $P^{\left(\text{best}\right)}\leftarrow P\left(K^{\left(0\right)}\right)$

3. ${\sigma }_{\mathit{max}}^{\left(\text{best}\right)}\leftarrow max\left(\sigma \left(P^{\left(\text{best}\right)}\right)\right)$

4. $\upsilon \leftarrow 1$

5. for $i=1$ to $m$ do

6. choose ${\mathcal{U}}^{\left(0\right)}$ such that ${‖{\mathcal{U}}^{\left(0\right)}‖}_F=1$ , uniformly at random

7. $L^{\left(0\right)}\leftarrow K^{\left(0\right)}+h\cdot {\mathcal{U}}^{\left(0\right)}$

8. for $j=1$ to $n$ do

9. choose $F\in {\mathcal{R}}_{\infty }\left(\mathrm{ϵ}\right)$ , uniformly at random

10. for $k=0$ downto $s$ do

11. $t\leftarrow \frac{k}{s}$

12. $K\left(t\right)\leftarrow \left(1-t\right)K^{\left(0\right)}+\mathit{tF}$

13. if $K\left(t\right)\in {\mathcal{S}}_{\alpha }$ then

14. compute $P\left(K\left(t\right)\right)$ as in (5)

15. ${\sigma }_{\mathit{max}}\left(K\left(t\right)\right)\leftarrow max\left(\sigma \left(P\left(K\left(t\right)\right)\right)\right)$

16. if $\left({\sigma }_{\mathit{max}}\left(K\left(t\right)\right)<{\sigma }_{\mathit{max}}^{\left(\text{best}\right)}\right)$ then

17. $K^{\left(\text{best}\right)}\leftarrow K\left(t\right)$

18. $P^{\left(\text{best}\right)}\leftarrow P\left(K\left(t\right)\right)$

19. ${\sigma }_{\mathit{max}}^{\left(\text{best}\right)}\leftarrow {\sigma }_{\mathit{max}}\left(K\left(t\right)\right)$

20. end if

21. end if

22. end for

23. end for

24. if $i>\upsilon \cdot ⌈e\sqrt{2\pi \mathit{qr}}⌉$ then

25. $K^{\left(0\right)}\leftarrow K^{\left(\text{best}\right)}$

26. $\upsilon \leftarrow \upsilon +1$

27. end for

28. end for

29. return $K^{\left(\text{best}\right)},P^{\left(\text{best}\right)},{\sigma }_{\mathit{max}}^{\left(\text{best}\right)}$

Remark 3.1. In [12] it is shown that by taking $m=⌈e\cdot \sqrt{2\mathit{\pi qr}}⌉$ iterations in the outer loop, we have $K_{\ast }\in D^{\left(0\right)}$ , for some direction ${\mathcal{U}}^{\left(0\right)}$ , almost surely. Let ${\mathcal{S}}_{\alpha }^{\left(0\right)}\left(\mathrm{ϵ}\right)$ denote the set $\left\{K\in {\mathcal{S}}_{\alpha }^{\left(0\right)}\left|{\sigma }_{\max }\left(K\right)\le {\sigma }_{\max }\left(K_{\ast }\right)+\mathrm{ϵ}\right.\right\}$ . Then, the total number of arithmetic operations of the RS algorithm that guarantees a probability of at least $1-\beta$ to hit ${\mathcal{S}}_{\alpha }^{\left(0\right)}\left(\mathrm{ϵ}\right)$ is given by $O\left(\frac{\left|ln\left(\beta \right)\right|h}{\mathrm{ϵ}}{\left(\frac{r_{\infty }}{r_{\mathrm{ϵ}}}\right)}^{q_0{r}_0}\left(max{\left(q,r\right)}^3+p^6\right)\right)$ , for systems with $q\le q_0,r\le r_0$ for fixed $q_0,r_0$ , where $r_{\mathrm{ϵ}}$ is the radius of the basis of a cone with height $\mathrm{ϵ}$ that has the same volume as of ${\mathcal{S}}_{\alpha }^{\left(0\right)}\left(\mathrm{ϵ}\right)$ ; see [10, 11, 12]. This is a polynomial-time algorithm by restricting the input and by regarding $\left(\frac{r_{\infty }}{r_{\mathrm{ϵ}}}\right)$ as the size of the problem.

4. The deterministic algorithm

The deterministic algorithm we introduce here as Algorithm 2 (which we call the MC algorithm) generalizes the algorithm of Daniel D. Moerder and Anthony A. Calise (see [15]) to the case of discrete-time systems. To the best of our knowledge, this is the best algorithm for LQR via SOF published so far, in terms of rate of convergence (to local minimum).

Here, we wish to minimize the LQR functional

$$J\left(x_0,P\right)=x_0^T{\mathit{Px}}_0,$$ E8

under the constraints

$$Y\left(K,P\right)≔Q+C^T{K}^T\mathit{RKC}-P+A_{c\mathrm{\ell }}{\left(K\right)}^T{\mathit{PA}}_{c\mathrm{\ell }}\left(K\right)=0,P>0.$$ E9

Since $Y^T=Y$ , there exist orthogonal matrix $\mathcal{U}$ such that $\widehat{Y}={\mathcal{U}}^{T}Y\mathcal{U}$ is diagonal. Now, minimizing (8) under the constraints (9) is equivalent to minimizing

$$\mathcal{L}\left(K,P,S\right)=\mathit{\text{trace}}\left(x_0^T{\mathit{Px}}_0\right)+\sum _{i=1}^p{\widehat{S}}_{i,i}{\widehat{Y}}_{i,i},$$

under the constraint $P>0$ , where ${\widehat{S}}_{i,i}$ are the Lagrange multipliers. We have

$$\begin{array}{cc}\mathcal{L}\left(K,P,S\right)& =\mathit{\text{trace}}\left(x_0^T{\mathit{Px}}_0\right)+\sum _{i=1}^p{\widehat{S}}_{i,i}{\widehat{Y}}_{i,i}\\ & =\mathit{\text{trace}}\left(x_0^T{\mathit{Px}}_0\right)+\mathit{\text{trace}}\left(\widehat{S}\widehat{Y}\right)\\ & =\mathit{\text{trace}}\left(x_0^T{\mathit{Px}}_0\right)+\mathit{\text{trace}}\left(\widehat{S}{\mathcal{U}}^{T}Y\mathcal{U}\right)\\ & =\mathit{\text{trace}}\left(x_0^T{\mathit{Px}}_0\right)+\mathit{\text{trace}}\left(\mathcal{U}\widehat{S}{\mathcal{U}}^{T}Y\right)\\ & =\mathit{\text{trace}}\left(x_0^T{\mathit{Px}}_0\right)+\mathit{\text{trace}}\left(\mathit{SY}\right)\end{array}$$

where $S=\mathcal{U}\widehat{S}{\mathcal{U}}^T$ . Note that $S^T=S$ . Let the Lagrangian be defined by

$$\mathcal{L}\left(K,P,S\right)=\mathit{\text{trace}}\left(x_0^T{\mathit{Px}}_0\right)+\mathit{\text{trace}}\left(\mathit{SY}\left(K,P\right)\right),$$ E10

for any $K$ any $P>0$ and any $S$ such that $S^T=S$ . The necessary conditions for optimality are $\frac{\partial \mathcal{L}}{\partial K}=0$ , $\frac{\partial \mathcal{L}}{\partial P}=0$ , and $\frac{\partial \mathcal{L}}{\partial S}=Y^T=Y=0$ .

Now, using Lemma 2.2, we have

$$\begin{array}{ll}& \frac{\partial \mathcal{L}}{\partial P}=0\\ & \iff x_0{x}_0^T-S^T+A_{c\mathrm{\ell }}{S}^T{A}_{c\mathrm{\ell }}^T=0\\ & \iff x_0{x}_0^T-S+A_{c\mathrm{\ell }}{\mathit{SA}}_{c\mathrm{\ell }}^T=0\\ & \iff S-A_{c\mathrm{\ell }}{\mathit{SA}}_{c\mathrm{\ell }}^T=x_0{x}_0^T\\ & \iff \left(I_p\otimes I_p-A_{c\mathrm{\ell }}\otimes A_{c\mathrm{\ell }}\right)\mathit{vec}\left(S\right)=\mathit{vec}\left(x_0{x}_0^T\right)\\ & \iff S=\mathit{mat}\left({\left(I_p\otimes I_p-A_{c\mathrm{\ell }}\otimes A_{c\mathrm{\ell }}\right)}^{-1}\mathit{vec}\left(x_0{x}_0^T\right)\right),\end{array}$$

where the last passage is affordable because $\sigma \left(A_{c\mathrm{\ell }}\right)\subset \mathbb{D}$ . Note that the last with the stability of $A_{c\mathrm{\ell }}$ implies that $S\ge 0$ .

We also have

$$\begin{array}{ll}\frac{\partial \mathcal{L}}{\partial K}& =\frac{\partial }{\partial K}\mathit{\text{trace}}\left(\mathit{SY}\right)\\ & =\frac{\partial }{\partial K}\mathit{\text{trace}}\left(S\left(Q+C^T{K}^T\mathit{RKC}-P+\left(A^T-C^T{K}^T{B}^T\right)P\left(A-\mathit{BKC}\right)\right)\right)\\ & =\frac{\partial }{\partial K}\mathit{\text{trace}}\left({\mathit{SC}}^T{K}^T\mathit{RKC}-{\mathit{SA}}^T\mathit{\text{PBKC}}-{\mathit{SC}}^T{K}^T{B}^T\mathit{PA}+{\mathit{SC}}^T{K}^T{B}^T\mathit{\text{PBKC}}\right)\\ & =\frac{\partial }{\partial K}\mathit{\text{trace}}\left({\mathit{SC}}^T{K}^T\mathit{RKC}-{\mathit{SA}}^T\mathit{\text{PBKC}}-A^T{P}^T{\mathit{\text{BKCS}}}^T+{\mathit{SC}}^T{K}^T{B}^T\mathit{\text{PBKC}}\right)\\ & =R^T{\mathit{KCS}}^T{C}^T+{\mathit{\text{RKCSC}}}^T-B^T{P}^T{\mathit{AS}}^T{C}^T-B^T{\mathit{\text{PASC}}}^T\\ & +B^T{P}^T{\mathit{\text{BKCS}}}^T{C}^T+B^T{\mathit{\text{PBKCSC}}}^T\\ & =2{\mathit{\text{RKCSC}}}^T-2B^T{\mathit{\text{PASC}}}^T+2B^T{\mathit{\text{PBKCSC}}}^T.\end{array}$$

Therefore,

$$\begin{array}{ll}& \frac{\partial \mathcal{L}}{\partial K}=0\\ & \iff {\mathit{\text{RKCSC}}}^T-B^T{\mathit{\text{PASC}}}^T+B^T{\mathit{\text{PBKCSC}}}^T=0\\ & \iff \left(R+B^T\mathit{PB}\right){\mathit{\text{KCSC}}}^T=B^T{\mathit{\text{PASC}}}^T\\ & \iff {\mathit{\text{KCSC}}}^T={\left(R+B^T\mathit{PB}\right)}^{-1}{B}^T{\mathit{\text{PASC}}}^T.\end{array}$$

Thus, if ${\mathit{CSC}}^T$ is invertible, then

$$K={\left(R+B^T\mathit{PB}\right)}^{-1}{B}^T{\mathit{\text{PASC}}}^T{\left({\mathit{CSC}}^T\right)}^{-1}.$$ E11

Otherwise, if

$${\left(R+B^T\mathit{PB}\right)}^{-1}{B}^T{\mathit{\text{PASC}}}^T\cdot L_{{\mathit{CSC}}^T}=0,$$

which is equivalent to

$$B^T{\mathit{\text{PASC}}}^T\cdot L_{{\mathit{CSC}}^T}=0,$$ E12

then

$$K={\left(R+B^T\mathit{PB}\right)}^{-1}{B}^T{\mathit{\text{PASC}}}^T{\left({\mathit{CSC}}^T\right)}^{+}+Z\cdot R_{{\mathit{CSC}}^T},$$ E13

where $Z$ is arbitrary $q\times r$ matrix (and we may take $Z=0$ , unless some other constraints on $K$ are needed). Note that if condition (12) does not happen, then $\frac{\partial \mathcal{L}}{\partial K}\ne 0$ . We conclude with the following theorem:

Theorem 4.1. Assume that $\mathcal{L}\left(K,P,S\right)$ given by (10) is minimized locally at some point $K_{\ast }$ , $P_{\ast }>0$ , and $S_{\ast }$ such that $S_{\ast }^T=S_{\ast }$ . Then

$$\left\{\begin{array}{l}{K}_{\ast }={\left(R+B^T{P}_{\ast }B\right)}^{-1}{B}^T{P}_{\ast }{\mathit{AS}}_{\ast }{C}^T{\left({\mathit{CS}}_{\ast }{C}^T\right)}^{+}+Z_{\ast }\cdot R_{{\mathit{CSC}}^T},\mathit{\text{for\hspace{0.17em}some}}q\times r\mathit{\text{matrix}}{Z}_{\ast }\\ P_{\ast }=\mathit{mat}\left({\left(I_p\otimes I_p-A_{c\mathrm{\ell }}{\left(K_{\ast }\right)}^T\otimes A_{c\mathrm{\ell }}{\left(K_{\ast }\right)}^T\right)}^{-1}\cdot \mathit{vec}\left(Q+C^T{K}_{\ast }^T{\mathit{RK}}_{\ast }C\right)\right)\\ S_{\ast }=\mathit{mat}\left({\left(I_p\otimes I_p-A_{c\mathrm{\ell }}\left(K_{\ast }\right)\otimes A_{c\mathrm{\ell }}\left(K_{\ast }\right)\right)}^{-1}\mathit{vec}\left(x_0{x}_0^T\right)\right),\end{array}\right.$$ E14

where $A_{c\mathrm{\ell }}\left(K_{\ast }\right)=A-{\mathit{BK}}_{\ast }C$ .

Proof:

Since $\mathcal{L}\left(K_{\ast },P_{\ast },S_{\ast }\right)$ is minimal in some neighborhood of $\left(K_{\ast },P_{\ast },S_{\ast }\right)$ , it follows that $\frac{\partial \mathcal{L}}{\partial K}\left(K_{\ast },P_{\ast },S_{\ast }\right)=0$ , $\frac{\partial \mathcal{L}}{\partial P}\left(K_{\ast },P_{\ast },S_{\ast }\right)=0$ , and $\frac{\partial \mathcal{L}}{\partial S}\left(K_{\ast },P_{\ast },S_{\ast }\right)=Y^T\left(K_{\ast },P_{\ast }\right)=Y\left(K_{\ast },P_{\ast }\right)=0$ .

The condition $\frac{\partial \mathcal{L}}{\partial S}\left(K_{\ast },P_{\ast },S_{\ast }\right)=Y^T\left(K_{\ast },P_{\ast }\right)=Y\left(K_{\ast },P_{\ast }\right)=0$ is just

$$P_{\ast }-A_{c\mathrm{\ell }}{\left(K_{\ast }\right)}^T{P}_{\ast }{A}_{c\mathrm{\ell }}\left(K_{\ast }\right)=Q+C^T{K}_{\ast }^T{\mathit{RK}}_{\ast }C$$

which with $P_{\ast }>0$ and $Q>0,R>0$ implies that $A_{c\mathrm{\ell }}\left(K_{\ast }\right)=A-{\mathit{BK}}_{\ast }C$ is stable. Now, since $\sigma \left(A_{c\mathrm{\ell }}\left(K_{\ast }\right)\right)\subset \mathbb{D}$ , it follows that $I_p\otimes I_p-A_{c\mathrm{\ell }}{\left(K_{\ast }\right)}^T\otimes A_{c\mathrm{\ell }}{\left(K_{\ast }\right)}^T$ is invertible, and therefore,

$$P_{\ast }=\mathit{mat}\left({\left(I_p\otimes I_p-A_{c\mathrm{\ell }}{\left(K_{\ast }\right)}^T\otimes A_{c\mathrm{\ell }}{\left(K_{\ast }\right)}^T\right)}^{-1}\cdot \mathit{vec}\left(Q+C^T{K}_{\ast }^T{\mathit{RK}}_{\ast }C\right)\right).$$

Since $I_p\otimes I_p-A_{c\mathrm{\ell }}\left(K_{\ast }\right)\otimes A_{c\mathrm{\ell }}\left(K_{\ast }\right)$ is invertible, $\frac{\partial \mathcal{L}}{\partial P}\left(K_{\ast },P_{\ast },S_{\ast }\right)=0$ implies that

$$S_{\ast }=\mathit{mat}\left({\left(I_p\otimes I_p-A_{c\mathrm{\ell }}\left(K_{\ast }\right)\otimes A_{c\mathrm{\ell }}\left(K_{\ast }\right)\right)}^{-1}\mathit{vec}\left(x_0{x}_0^T\right)\right).$$

Finally, $\frac{\partial \mathcal{L}}{\partial P}\left(K_{\ast },P_{\ast },S_{\ast }\right)=0$ implies that $K_{\ast }{\mathit{CS}}_{\ast }{C}^T={\left(R+B^T{P}_{\ast }B\right)}^{-1}{B}^T{P}_{\ast }{\mathit{AS}}_{\ast }{C}^T$ , which in view of Lemma 2.1 implies ${\left(R+B^T{P}_{\ast }B\right)}^{-1}{B}^T{P}_{\ast }{\mathit{AS}}_{\ast }{C}^T\cdot L_{{\mathit{CS}}_{\ast }{C}^T}=0$ and

$$K_{\ast }={\left(R+B^T{P}_{\ast }B\right)}^{-1}{B}^T{P}_{\ast }{\mathit{AS}}_{\ast }{C}^T{\left({\mathit{CS}}_{\ast }{C}^T\right)}^{+}+Z_{\ast }\cdot R_{{\mathit{CS}}_{\ast }{C}^T},$$

where $Z_{\ast }$ is some $q\times r$ matrix. $■$

Note that the equations are coupled tightly, in the sense that $P_{\ast }$ and $S_{\ast }$ need $K_{\ast }$ , while $K_{\ast }$ needs $P_{\ast }$ and $S_{\ast }$ . Note also the cubic dependencies (that can be made quadratic by introducing new variables). These make the related QMEs non-convex and, therefore, hard to compute.

Remark 4.1. When $x_0$ is unknown, it is customary to assume that $x_0$ is uniformly distributed on the unit sphere, which implies that $E\left[x_0{x}_0^T\right]=I_p$ , where $E\left[•\right]$ is the expectation operator. Thus, changing the problem to that of minimizing $E\left[J\left(x_0,P\right)\right]$ amounts to replacing $S_{\ast }$ with

$$E\left[S_{\ast }\right]=\mathit{mat}\left({\left(I_p\otimes I_p-A_{c\mathrm{\ell }}\left(K_{\ast }\right)\otimes A_{c\mathrm{\ell }}\left(K_{\ast }\right)\right)}^{-1}\mathit{vec}\left(I_p\right)\right)>0.$$

Therefore, there is change in Algorithm 2.

Remark 4.2. The convergence of Algorithm 2 to local minimum can be proved similarly to the proof appearing in [15], under the assumptions that ${\mathcal{S}}_{\alpha }^{q\times r}$ is nonempty and that $\left(Q^{\frac{1}{2}},A\right)$ is detectable (here, we do not need this condition because of the assumption that $Q>0$ ). The convergence can actually be proved for the more general problem that adds ${‖K‖}_F^2=\mathit{\text{trace}}\left(K^{T}K\right)$ to the LQR functional, thus minimizing also the Frobenius norm of $K$ . In this context, note that by adding ${‖K‖}^2=max\left(\sigma \left(K^{T}K\right)\right)$ to the LQR functional will lose the argument, and there will be a need to more general proof, because in the proof appearing in [15], the demand is for ${\mathbf{C}}^1$ smooth function of $K$ , while ${‖K‖}^2={\sigma }_{\max }\left(K^{T}K\right)$ is continuous but not Lipschitz continuous. The RS algorithm can use any continuous function of $K$ and can deal also with sparse SOFs for LQR and with regional pole-placement SOFs for LQR.

Example 4.1. In the following simple example, we illustrate the notions appearing in the definition of the RS algorithm, and we demonstrate the operation of the RS algorithm. Consider the unstable system

$$\left\{\begin{array}{l}{x}_{k+1}=\left[\begin{array}{cc}2& 1\\ 0& -\frac{1}{2}\end{array}\right]x_k+\left[\begin{array}{c}1\\ 1\end{array}\right]u_k,k=0,1,\dots \\ y_k=x_k,\end{array}\right.$$

where we look for SOF $K$ stabilizing the system while reducing the LQR functional (2) with $Q=I,R=1$ . Let $K=\left[k_1{k}_2\right]$ then

$$A_{c\mathrm{\ell }}\left(K\right)=A-\mathit{BK}=\left[\begin{array}{cc}2-k_1& 1-k_2\\ -k_1& -\frac{1}{2}-k_2\end{array}\right],$$

with characteristic polynomial $z^2+z\left(k_1+k_2-\frac{3}{2}\right)+\frac{3}{2}{k}_1-2k_2-1$ . Applying the Y. Bistritz stability criterion (see [25]), we have

$$\upsilon =\mathit{Var}\left\{\frac{5}{2}{k}_1-k_2-\frac{3}{2},-\frac{3}{2}{k}_1+2k_2+2,\frac{1}{2}{k}_1-3k_2+\frac{3}{2}\right\},$$

where $\upsilon$ is the number of sign variations in the set. According to the Bistritz criterion, the system is stable if and only if $\upsilon =0$ . We conclude that $\mathcal{S}$ is the set of all $K$ such that $\frac{5}{2}{k}_1-k_2-\frac{3}{2}>0,-\frac{3}{2}{k}_1+2k_2+2>0,\frac{1}{2}{k}_1-3k_2+\frac{3}{2}>0$ or $\frac{5}{2}{k}_1-k_2-\frac{3}{2}<0,-\frac{3}{2}{k}_1+2k_2+2<0,\frac{1}{2}{k}_1-3k_2+\frac{3}{2}<0$ , where the last branch is empty (which could have make the set non-convex). The set $\mathcal{S}$ appears in Figure 1 as the blue region, where the golden star is the analytic global optimal solution $K_{\ast }=\left[\begin{array}{ll}1.09473459& 0.36138828\end{array}\right]$ (computed by the related discrete algebraic Riccati equation).

Figure 1.

The stability region $S$ .

Algorithm 2 The MC algorithm for LQR via SOF for discrete-time systems.

Require: An algorithm for deciding $\alpha$ -stability, an algorithm for computing ${\sigma }_{\max }\left(K\right)$ and algorithms for general linear algebra operations.

Input: $0<\mathrm{ϵ}\le \frac{1}{2},0<\alpha <1$ , integers: $m,s>0$ , controllable pairs $\left(A,B\right)$ and $\left(A^T,C^T\right)$ , matrices $Q>0,R>0$ and $K^{\left(0\right)}\in \mathit{int}\left({\mathcal{S}}_{\alpha }\right)$ .

Output: $K\in {\mathcal{S}}_{\alpha }$ that locally minimizes ${\sigma }_{\mathit{max}}\left(K\right)$ .

1. $j\leftarrow 0;A_0\leftarrow A-{\mathit{BK}}_{0}C$

2. $P_0\leftarrow \mathit{mat}\left(\begin{array}{c}{\left(I_p\otimes I_p-A_0^T\otimes A_0^T\right)}^{-1}\cdot \mathit{vec}\left(Q+C^T{K}_0^T{\mathit{RK}}_{0}C\right)\end{array}\right)$

3. $S_0\leftarrow \mathit{mat}\left({\left(I_p\otimes I_p-A_0\otimes A_0\right)}^{-1}\cdot \mathit{vec}\left(I_p\right)\right);{\sigma }_{\mathit{max}}\left(K_0\right)\leftarrow max\left(\sigma \left(P_0\right)\right)$

4. $\mathrm{\Delta }{K}_0\leftarrow {\left(R+B^T{P}_{0}B\right)}^{-1}{B}^T{P}_0{\mathit{AS}}_0{C}^T{\left({\mathit{CS}}_{0}C\right)}^{+}-K_0$

5. $\mathit{\text{flag}}\leftarrow 0$

6. for $k=0$ to $s$ do

7. $t\leftarrow \frac{k}{s}$

8. $K\left(t\right)\leftarrow \left(1-t\right)K_0+t\mathrm{\Delta }{K}_0$

9. if $K\left(t\right)\in {\mathcal{S}}_{\alpha }$ then

10. $A\left(t\right)\leftarrow A-\mathit{BK}\left(t\right)C$

11. $P\left(t\right)\leftarrow \mathit{mat}\left(\begin{array}{c}{\left(I_p\otimes I_p-A{\left(t\right)}^T\otimes A{\left(t\right)}^T\right)}^{-1}\cdot \mathit{vec}\left(Q+C^{T}K{\left(t\right)}^T\mathit{RK}\left(t\right)C\right)\end{array}\right)$

12. $S\left(t\right)\leftarrow \mathit{mat}\left({\left(I_p\otimes I_p-A\left(t\right)\otimes A\left(t\right)\right)}^{-1}\cdot \mathit{vec}\left(I_p\right)\right);{\sigma }_{\mathit{max}}\left(K\left(t\right)\right)\leftarrow max\left(\sigma \left(P\left(t\right)\right)\right)$

13. if ${\sigma }_{\mathit{max}}\left(K\left(t\right)\right)<{\sigma }_{\mathit{max}}\left(K_0\right)$ then

14. $K_1\leftarrow K\left(t\right);A_1\leftarrow A-{\mathit{BK}}_{1}C;P_1\leftarrow P\left(t\right);S_1\leftarrow S\left(t\right);{\sigma }_{\mathit{max}}\left(K_1\right)\leftarrow {\sigma }_{\mathit{max}}\left(K\left(t\right)\right)$

15. $\mathit{\text{flag}}\leftarrow 1$

16. end if

17. end if

18. end for

19. if $\mathit{\text{flag}}==1$ then

20. while $\left(\left|{\sigma }_{\mathit{max}}\left(K_{j+1}\right)-{\sigma }_{\mathit{max}}\left(K_j\right)\right|\ge \mathrm{ϵ}\right)$ and $\left(j<m\right)$ do

21. $\mathrm{\Delta }{K}_j\leftarrow {\left(R+B^T{P}_{j}B\right)}^{-1}{B}^T{P}_j{\mathit{AS}}_j{C}^T{\left({\mathit{CS}}_j{C}^T\right)}^{+}-K_j$

22. for $k=0$ to $s$ do

23. $t\leftarrow \frac{k}{s}$

24. $K\left(t\right)\leftarrow \left(1-t\right)K_j+t\mathrm{\Delta }{K}_j$

25. if $K\left(t\right)\in {\mathcal{S}}_{\alpha }$ then

21. $A\left(t\right)\leftarrow A-\mathit{BK}\left(t\right)C$

22. $P\left(t\right)\leftarrow \mathit{mat}\left(\begin{array}{c}{\left(I_p\otimes I_p-A{\left(t\right)}^T\otimes A{\left(t\right)}^T\right)}^{-1}\cdot \mathit{vec}\left(Q+C^{T}K{\left(t\right)}^T\mathit{RK}\left(t\right)C\right)\end{array}\right)$

23. $S\left(t\right)\leftarrow \mathit{mat}\left({\left(I_p\otimes I_p-A\left(t\right)\otimes A\left(t\right)\right)}^{-1}\cdot \mathit{vec}\left(I_p\right)\right);{\sigma }_{\mathit{max}}\left(K\left(t\right)\right)\leftarrow max\left(\sigma \left(P\left(t\right)\right)\right)$

29. if ${\sigma }_{\mathit{max}}\left(K\left(t\right)\right)<{\sigma }_{\mathit{max}}\left(K_j\right)$ then

30. $K_{j+1}\leftarrow K\left(t\right);A_{j+1}\leftarrow A-{\mathit{BK}}_{j+1}C;P_{j+1}\leftarrow P\left(t\right);S_{j+1}\leftarrow S\left(t\right);{\sigma }_{\mathit{max}}\left(K_{j+1}\right)\leftarrow {\sigma }_{\mathit{max}}\left(K\left(t\right)\right)$

31. end if

32. end if

33. end for

34. $j\leftarrow j+1$

35. end while

36. end if

37. return $K^{\left(\text{best}\right)}\leftarrow K_j,P^{\left(\text{best}\right)}\leftarrow P_j,{\sigma }_{\mathit{max}}^{\left(\text{best}\right)}\leftarrow {\sigma }_{\mathit{max}}\left(K_j\right)$

In Figure 1, we can see how the RS algorithm works: we fix $\alpha ={10}^{-3},\mathrm{ϵ}={10}^{-16},r_{\infty }=2,h=2$ , and we set $m=1,n=5,s=10$ for a single iteration, where the single cone is sampled along $5$ rays and each ray is sampled $10$ times. The sampled points are circled, where red circles indicate infeasible or non-improving points and the black circles indicate improving points. The green star point is the initial point $K^{\left(0\right)}$ found by the Ray-Shooting algorithm for minimal-norm SOF. The bold black circle is the boundary of the closed circle $\overline{\mathbb{B}\left(K^{\left(0\right)},h\right)}$ . We choose ${\mathcal{U}}^{\left(0\right)}$ randomly to define the search direction, and we set $L^{\left(0\right)}=K^{\left(0\right)}+h\cdot {\mathcal{U}}^{\left(0\right)}$ to be the point where the direction meets the boundary of the circle. $L$ is the tangent line at $L^{\left(0\right)}$ to the circle, and ${\mathcal{R}}_{\infty }\left(\mathrm{ϵ}\right)$ is the $2r_{\infty }$ width segment on the line, inflated by $\mathrm{ϵ}$ . The search cone ${\mathcal{D}}^{\left(0\right)}=\mathit{\text{chull}}\left(K^{\left(0\right)},R_{\infty }\left(\mathrm{ϵ}\right)\right)$ is the related black triangle. Here ${\mathcal{S}}_{\alpha }^{\left(0\right)}={\mathcal{S}}_{\alpha }\cap {\mathcal{D}}^{\left(0\right)}$ is just the portion of the blue region inside the triangle, and we can see that the assumption that $K_{\ast }\in {\mathcal{D}}^{\left(0\right)}$ is in force. For the current problem $⌈e\sqrt{2\pi \mathit{qr}}⌉=10$ , and therefore, by making $10$ iterations, $K_{\ast }$ will be inside some triangle almost surely.

The algorithm chooses $F$ in the basis of the triangle and defines $K\left(t\right)$ to be the ray from $K_0$ to $F$ . The ray is sampled at $10$ equally spaced points, and the best feasible point is recorded.

In Figure 2, we can see that $5$ iterations suffice to include $K_{\ast }$ in some triangle and to find improving points very close to $K_{\ast }$ . In Figure 3, we can see that when we allow $20$ iterations, after $10$ iterations, the center $K^{\left(0\right)}$ is switched to the best point found so far (see lines $24-26$ in Algorithm 1). This is done in order to raise the probability to hit $K_{\ast }$ or its $\mathrm{ϵ}$ -neighborhood, and as we can see, the final best point (green star) is very close to $K_{\ast }$ (Figure 4).

Figure 2.

Single iteration of the RS algorithm.

Figure 3.

Five iterations of the RS algorithm.

Figure 4.

Twenty iterations of the RS algorithm.

The results of the algorithm for $1,5$ and $20$ iterations are the following. Note that ${\sigma }_{\max }\left(K_{\ast }\right)=5.9551$ , and note the “huge variations” the function ${\sigma }_{\mathit{max}}\left(K\right)$ has.

For $m=1$ we had

$$\begin{array}{ll}& K^{\left(0\right)}=\left[\begin{array}{rr}0.58739333& -0.15823016\end{array}\right],{\sigma }_{\max }^{\left(0\right)}=25.7307,\\ & \mathbf{RS}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}1.17786349& 0.35034398\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=6.1391,\\ & \mathbf{MC}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}0.58739333& -0.15823016\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=25.7307,\\ & \mathbf{RS}\mathbf{+}\mathbf{MC}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}1.05244278& 0.31681948\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=6.0001.\\ & \end{array}$$

Note that in this case, the $\mathit{MC}$ algorithm makes no improvement, while the RS and RS + MC have very close values to the global optimal value, with slightly better value for the RS + MC, over the RS algorithm.

For $m=5$ we had

$$\begin{array}{ll}& K^{\left(0\right)}=\left[\begin{array}{rr}0.60478870& -0.06023828\end{array}\right],{\sigma }_{\max }^{\left(0\right)}=36.4583,\\ & \mathbf{RS}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}1.04166520& 0.40826562\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=6.1655,\\ & \mathbf{MC}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}0.60478870& -0.06023828\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=36.4583,\\ & \mathbf{RS}\mathbf{+}\mathbf{MC}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}1.04166520& 0.40826562\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=6.16557843.\end{array}$$

For $m=20$ we had

$$\begin{array}{ll}& K^{\left(0\right)}=\left[\begin{array}{rr}0.51029365& -0.22521376\end{array}\right],{\sigma }_{\max }^{\left(0\right)}=3198.8196,\\ & \mathbf{RS}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}1.11453066& 0.33955607\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=5.9893,\\ & \mathbf{MC}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}0.51029365& -0.22521376\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=3198.8196,\\ & \mathbf{RS}+\mathbf{MC}:K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{rr}1.11453066& 0.33955607\end{array}\right],{\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=5.9893.\end{array}$$

In Figure 5, the initial condition response of the open-loop system is given. One can see the unstable mode related to the unstable eigenvalue $2$ . In Figures 6–8, the initial condition responses of the closed-loop systems with the SOFs for $m=20$ , with $x_0={\left[\begin{array}{rr}3& 1\end{array}\right]}^T$ and sampling time $T_s=0.01$ , are given. One can see that the responses of the closed-loop systems with the SOFs computed by RS and RS + MC are very close to the global optimal response, while the response of the closed-loop system with the SOF computed by the MC algorithm (actually with the initial SOF), although stable, is unacceptable.

Figure 5.

The initial condition response of the open-loop system.

Figure 6.

The initial condition response of the closed-loop system with the SOF computed by the RS algorithm (blue) compared with the global optimal response (red).

Figure 7.

The initial condition response of the closed-loop system with the SOF computed by the MC algorithm (blue) compared with the global optimal response (red).

Figure 8.

The initial condition response of the closed-loop system with the SOF computed by the RS + MC algorithm (blue) compared with the global optimal response (red).

5. Experiments

In the following experiments, we applied Algorithms 1 and 2, for systems taken from the libraries [26, 27, 28]. The systems given in these libraries are real-life continuous-time systems. In order to get related discrete-time systems, we sampled the systems using the Tustin method with sampling rate $T_s=0.01\left[sec\right]$ . We took only the systems for which the RS algorithm succeeded in finding SOFs for the continuous-time systems (see [10], Table 8, p. 231). In order to initialize the MC Algorithm, we also used the RS algorithm to find a starting $\alpha$ -stabilizing SOF. In all the experiments, we used $m=2⌈e\sqrt{2\pi \mathit{qr}}⌉,n=100,s=100$ ; $h=100,r_{\infty }=100,\mathrm{ϵ}={10}^{-16}$ , for the RS algorithm; and $m=200⌈e\sqrt{2\pi \mathit{qr}}⌉,s=100$ , for the MC Algorithm, in order to get the same number of total iterations and the same number $s=100$ of iterations for the local search. We took $Q=I_p,R=I_q$ in all the cases.

The stability margin column of Table 1 relates to $0<\alpha <1$ for which the absolute value of any eigenvalue of the closed loop is less than $1-\alpha$ . The values of $\alpha$ in Table 1 relates to the largest $0<\alpha <1$ for which the RS algorithm succeeded in finding a starting SOF $K^{\left(0\right)}$ . As we saw above, it is worth searching for a starting point $K^{\left(0\right)}$ that maximizes $0<\alpha <1$ . This can be achieved efficiently by running a binary search on the $0<\alpha <1$ and using the RS algorithm as an oracle. Note that the RS CPU time appearing in the fourth column of Table 1 relates to running the RS algorithm for known optimal value of $0<\alpha <1$ . The RS algorithm is sufficiently fast also for this purpose, but other algorithms such as the HIFOO (see [24]) and HINFSTRUCT (see [29]) can be applied in order to get a starting SOF. The advantage of the use of the RS is of finding starting SOF with relatively small norm.

Let ${\sigma }_{\mathit{max}}\left(F\right)$ denote the functional (7) for the system $\left(A,B,I_p\right)$ , where $A-\mathit{BF}$ is stable, i.e., $F\in {\mathcal{S}}^{q\times p}$ . Let $P\left(F\right)$ denote the Lyapunov matrix (5) for the system $\left(A,B,I_p\right)$ with $F$ as above. Let ${\sigma }_{\mathit{max}}\left(K\right)$ denote the functional (7) for the system $\left(A,B,C\right)$ with $K\in {\mathcal{S}}^{q\times r}$ and related Lyapunov matrix $P=P\left(K\right)$ given by (5). Now, if $A-\mathit{BKC}$ is stable for some $K$ , then $A-\mathit{BF}$ is stable for $F=\mathit{KC}$ (but there might exist $F$ such that $A-\mathit{BF}$ is stable, which cannot be defined as $\mathit{KC}$ for some $q\times r$ matrix $K$ ). Therefore,

$${\sigma }_{\mathit{max}}\left(F_{\ast }\right)=\underset{F\in {\mathcal{S}}^{q\times p}}{min}{\sigma }_{\mathit{max}}\left(F\right)\le \underset{K\in {\mathcal{S}}_{\alpha }^{q\times r}\cap \overline{\mathbb{B}\left(K^{\left(0\right)},h\right)}}{min}{\sigma }_{\mathit{max}}\left(K\right)={\sigma }_{\mathit{max}}\left(K_{\ast }\right),$$ E15

where $F_{\ast }$ is an optimal LQR state-feedback (SF) controller for the system $\left(A,B,I_p\right)$ . We conclude that ${\sigma }_{\mathit{max}}\left(F_{\ast }\right)\le {\sigma }_{\mathit{max}}\left(K_{\ast }\right)\le {\sigma }_{\mathit{max}}\left(K^{\left(\text{best}\right)}\right)$ . Thus, ${\sigma }_{\mathit{max}}\left(F_{\ast }\right)$ is a lower bound for ${\sigma }_{\mathit{max}}\left(K^{\left(\mathit{\text{best}}\right)}\right)$ and can serve as a good estimator for it, in order to quantify the convergence of the algorithm to the global minimum (as is evidently seen from Table 1 in many cases) and in order to stop the algorithm earlier, since ${\sigma }_{\mathit{max}}\left(F_{\ast }\right)$ can be calculated in advance. The lower bound appears in the last column of Table 1.

For all the systems, we had $\left(A,B\right),\left(A^T,C^T\right)$ controllable, except for ROC1 and ROC2. All the systems are unstable, except for AC6, AC15, and NN16 which are stable, but not $\alpha$ -stable, for $\alpha$ given in the stability margin column.

The experiments were executed on:

Computer: LAPTOP-GULIHG OV, ASUSeK COMPUTER, INC.

TUF GAMING FX504GM-FX80GM.

Processor: Intel(R) Core(TM) i7-8750H CPU@2.20GHz.

Platform: MATLAB, Version R2018b Win 64.

5.1 Conclusions from the experiments

Regarding the experimental results in Table 2 and the comparison between the RS algorithm and the MC algorithm, we conclude:

1. The RS algorithm performs in magnitude better than the MC algorithm for the systems: AC1, AC11, HE1, HE4, ROC1, ROC4, TF1, and NN5.

2. The MC algorithm performs in magnitude better than the RS algorithm for the systems AC5 and DIS5.

3. The MC algorithm performs slightly better than the RS algorithm for systems HE3 and NN16.

System	$\begin{array}{l}\text{Size}\\ \left(p,q,r\right)\end{array}$	$\begin{array}{l}\text{Stab}.\\ \mathrm{Mgn}.\end{array}$	$\begin{array}{l}\mathrm{RS}\\ \mathrm{CPU}\\ \text{time}\\ \left[sec\right]\end{array}$	$\begin{array}{l}{\sigma }_{\max }^{\left(0\right)}\\ \text{for}\\ \left(A,B,C\right)\end{array}$	$\begin{array}{l}{\sigma }_{\max }\left(F_{\ast }\right)\\ \text{for}\\ \left(A,B\right)\end{array}$
AC1	$\left(\mathrm{5,3,3}\right)$	$0.01$	$2.6226$	$1.0701\cdot {10}^4$	$1.3073\cdot {10}^3$
AC5	$\left(\mathrm{4,2,2}\right)$	$0.001$	$1.5468$	$1.5888\cdot {10}^9$	$8.4264\cdot {10}^7$
AC6	$\left(\mathrm{7,2,4}\right)$	$0.001$	$0.7094$	$3.1767\cdot {10}^3$	$5.9783\cdot {10}^2$
AC11	$\left(\mathrm{5,2,4}\right)$	$0.01$	$1.0575$	$1.2968\cdot {10}^4$	$5.8777\cdot {10}^2$
HE1	$\left(\mathrm{4,2,1}\right)$	$0.001$	$0.0872$	$1.5040\cdot {10}^3$	$3.0013\cdot {10}^2$
HE3	$\left(\mathrm{8,4,6}\right)$	$0.001$	$2.6845$	$5.4064\cdot {10}^6$	$6.1185\cdot {10}^4$
HE4	$\left(\mathrm{8,4,6}\right)$	$0.001$	$2.5633$	$4.1660\cdot {10}^6$	$2.2992\cdot {10}^4$
ROC1	$\left(\mathrm{9,2,2}\right)$	${10}^{-5}$	$0.5279$	$1.5906\cdot {10}^7$	$1.1207\cdot {10}^5$
ROC4	$\left(\mathrm{9,2,2}\right)$	${10}^{-5}$	$0.4677$	$1.2273\cdot {10}^6$	$8.5460\cdot {10}^4$
DIS4	$\left(\mathrm{8,4,6}\right)$	$0.01$	$2.5074$	$4.5133\cdot {10}^3$	$1.7556\cdot {10}^2$
DIS5	$\left(\mathrm{4,2,2}\right)$	$0.001$	$1.2187$	$2.8686\cdot {10}^8$	$9.0756\cdot {10}^6$
TF1	$\left(\mathrm{7,2,4}\right)$	${10}^{-4}$	$0.8011$	$7.9884\cdot {10}^5$	$5.8134\cdot {10}^3$
NN5	$\left(\mathrm{7,1,2}\right)$	${10}^{-4}$	$0.4138$	$5.4066\cdot {10}^6$	$2.8789\cdot {10}^5$
NN13	$\left(\mathrm{6,2,2}\right)$	$0.01$	$0.4876$	$7.8402\cdot {10}^2$	$63.5366$
NN16	$\left(\mathrm{8,4,4}\right)$	${10}^{-4}$	$3.5530$	$1.9688\cdot {10}^3$	$2.3327\cdot {10}^2$
NN17	$\left(\mathrm{3,2,1}\right)$	$0.001$	$0.0925$	$3.2733\cdot {10}^4$	$3.1358\cdot {10}^2$

Table 1.

General information of the systems and initial values.

System	$\begin{array}{l}{\sigma }_{\max }^{\left(\text{best}\right)}\\ \text{for}\\ \left(A,B,C\right)\\ \mathrm{RS}\\ \text{Algo}.\end{array}$	$\begin{array}{l}{\sigma }_{\max }^{\left(\text{best}\right)}\\ \text{for}\\ \left(A,B,C\right)\\ \mathrm{MC}\\ \text{Algo}.\end{array}$	$\begin{array}{l}{\sigma }_{\max }^{\left(\text{best}\right)}\\ \text{for}\\ \left(A,B,C\right)\\ \mathrm{RS}+\mathrm{MC}\\ \text{Algo}.\end{array}$	$\begin{array}{l}\mathrm{RS}\\ \text{Algo}.\\ \mathrm{CPU}\\ \text{time}\\ \left[sec\right]\end{array}$	$\begin{array}{l}\mathrm{MC}\\ \text{Algo}.\\ \mathrm{CPU}\\ \text{time}\\ \left[sec\right]\end{array}$	$\begin{array}{l}\mathrm{RS}+\mathrm{MC}\\ \text{Algo}.\\ \mathrm{CPU}\\ \text{time}\\ \left[sec\right]\end{array}$
AC1	$1.9207\cdot {10}^3$	$1.0701\cdot {10}^4$	$1.9207\cdot {10}^3$	$2.9843$	$0.0468$	$3.0312$
AC5	$1.5888\cdot {10}^9$	$2.5905\cdot {10}^8$	$2.5905\cdot {10}^8$	$2.0156$	$0.4062$	$2.2500$
AC6	$6.1449\cdot {10}^2$	$6.5913\cdot {10}^2$	$6.1389\cdot {10}^2$	$5.1250$	$0.2500$	$5.1718$
AC11	$2.4234\cdot {10}^3$	$1.2968\cdot {10}^4$	$2.4234\cdot {10}^3$	$2.9531$	$0.0468$	$3.0000$
HE1	$9.1253\cdot {10}^2$	$1.0968\cdot {10}^3$	$9.1253\cdot {10}^2$	$2.2812$	$0.0468$	$2.3437$
HE3	$8.6808\cdot {10}^4$	$7.1816\cdot {10}^4$	$8.1737\cdot {10}^4$	$13.4843$	$0.2343$	$13.7656$
HE4	$5.2247\cdot {10}^4$	$1.1817\cdot {10}^6$	$3.1783\cdot {10}^4$	$10.9687$	$0.0468$	$11.2656$
ROC1	$6.6239\cdot {10}^5$	$1.5906\cdot {10}^7$	$6.6239\cdot {10}^5$	$7.6250$	$0.1250$	$7.7500$
ROC4	$5.9923\cdot {10}^5$	$1.2273\cdot {10}^6$	$5.9923\cdot {10}^5$	$5.3750$	$0.0625$	$5.4218$
DIS4	$1.7590\cdot {10}^2$	$2.0376\cdot {10}^2$	$1.7590\cdot {10}^2$	$7.2187$	$0.1250$	$7.2656$
DIS5	$2.8686\cdot {10}^8$	$3.2079\cdot {10}^7$	$3.2079\cdot {10}^7$	$1.8593$	$0.1562$	$1.9687$
TF1	$1.9289\cdot {10}^4$	$1.1230\cdot {10}^5$	$1.9270\cdot {10}^4$	$5.0000$	$0.0625$	$5.0468$
NN5	$9.6780\cdot {10}^5$	$1.3372\cdot {10}^6$	$9.6780\cdot {10}^5$	$1.2968$	$0.1718$	$1.3593$
NN13	$2.0521\cdot {10}^2$	$2.6553\cdot {10}^2$	$1.7953\cdot {10}^2$	$2.3281$	$0.2656$	$2.4843$
NN16	$6.4416\cdot {10}^2$	$6.0032\cdot {10}^2$	$6.0030\cdot {10}^2$	$6.4062$	$0.5468$	$6.7343$
NN17	$3.6805\cdot {10}^3$	$4.9674\cdot {10}^3$	$3.6787\cdot {10}^3$	$0.9375$	$0.2812$	$1.1718$

Table 2.

Experimental results.

Regarding the experimental results in Table 2 and the performance of the RS + MC algorithm, we conclude:

1. The RS + MC algorithm performs better than each algorithm separately, for systems AC6, HE4, TF1, NN13, NN16, and NN17.

2. The RS + MC algorithm performs better than the RS algorithm for systems AC5, HE3, and DIS5.

3. The RS + MC algorithm performs exactly as the RS algorithm for systems AC1, AC11, DIS4, HE1, ROC1, ROC4, and NN5. This observation assesses the claim for convergence of the RS algorithm to global optimum.

4. The RS + MC algorithm performs exactly as the MC algorithm for systems AC5 and DIS5.

Regarding improvements over the starting point, we had:

1. The RS algorithm failed in finding any improvement over ${\sigma }_{\max }\left(K^{\left(0\right)}\right)$ for systems AC5 and DIS5.

2. The MC algorithm failed in finding any improvement over ${\sigma }_{\max }\left(K^{\left(0\right)}\right)$ for systems AC1, AC11, ROC1, and ROC4. This observation assesses the heuristic that it is better to start with a SOF that brings the poles of the closed loop as close as possible to the boundary of the disk ${\mathbb{D}}_{\alpha }$ .

3. The RS + MC algorithm improved ${\sigma }_{\max }\left(K^{\left(0\right)}\right)$ in any case.

Regarding the assessment of convergence to a global minimum, we had the following results:

1. The RS algorithm and the RS + MC algorithm had very close values of ${\sigma }_{\max }\left(K^{\left(\mathit{\text{best}}\right)}\right)$ (or exactly the same value) which is very close to the lower bound, for systems AC1, AC6, HE1, HE3, HE4, ROC1, DIS4, NN5, and NN16.

2. The MC algorithm achieved a very close value of ${\sigma }_{\max }\left(K^{\left(\mathit{\text{best}}\right)}\right)$ to the lower bound, for the systems AC6, HE3, DIS5, NN5, and NN16.

As was expected, the MC algorithm seems to perform better locally, while the RS algorithm seems to perform better globally. Thus, practically, the best approach is to apply the RS algorithm in order to find a close neighborhood of a global minimum and then to apply the MC algorithm on the result, for the local optimization, as is evidently seen from the performance of the RS + MC algorithm.

5.2 Some specific results

Example 5.1. For the HE4 system with

$$\begin{array}{l}\begin{array}{ll}& A=\left[\begin{array}{cc}{A}_1& A_2\end{array}\right],\text{where}\\ & A_1=\left[\begin{array}{rrrr}0.99999985& 0.00000014& 0.00001538& 0.00988556\\ -0.00000082& 0.99999927& 0.00944714& -0.00015179\\ -0.00016276& -0.00014358& 0.89058607& -0.02380208\\ -0.00002661& 0.00002887& 0.00410983& 0.98016538\\ -0.00002930& -0.00000576& -0.01923117& -0.00428369\\ -0.32100350& -0.00003189& -0.00471073& 0.02122199\\ 0.00098904& 0.32051910& -0.02076100& -0.00474084\\ -0.01910442& 0.01711196& 0.00004146& -0.00066123\end{array}\right],\\ & A_2=\left[\begin{array}{rrrr}0.00053188& 0.00000088& 0.00000089& -0.00000005\\ 0.00059005& 0.00000508& -0.00000451& 0.00000034\\ -0.00060285& 0.00100715& -0.00089937& 0.00006730\\ -0.00000054& 0.00016706& 0.00018078& -0.00001158\\ 0.99268256& 0.00018120& -0.00003706& 0.00002045\\ -0.00008474& 0.99978709& -0.00020742& 0.00015757\\ 0.00841925& 0.00020153& 0.99963047& 0.00000293\\ -0.00000005& 0.00013964& -0.00000914& 0.99709909\end{array}\right],\\ & \end{array}\\ \begin{array}{ll}& B=\left[\begin{array}{rrrr}-0.00000097& 0.00002352& -0.00000082& -0.00000055\\ 0.00000679& 0.00000357& -0.00013147& -0.00000146\\ 0.00117865& 0.00072316& -0.02601681& -0.00016953\\ -0.00035692& 0.00470509& 0.00008442& -0.00000015\\ 0.00302236& 0.00013040& -0.00468257& -0.00205817\\ 0.00286635& -0.00539604& -0.00009665& 0.00000026\\ -0.00018970& 0.00014386& -0.00517989& 0.00234114\\ -0.04813606& -0.00000574& -0.00000060& -0.00000001\end{array}\right],\\ & \end{array}\end{array}$$

$$\begin{array}{ll}& C=\left[\begin{array}{cc}{C}_1& C_2\end{array}\right],\text{where}\\ & C_1=\left[\begin{array}{rrrr}-0.00000185& 0.00001067& -0.00071398& 0.00083459\\ 0.99999992& 0.00000007& 0.00000769& 0.00494278\\ -0.00000041& 0.99999963& 0.00472357& -0.00007589\\ -0.00001393& -0.00000365& -0.00972548& -0.05509146\\ -0.00008138& -0.00007179& 0.94529303& -0.01190104\\ -0.00001330& 0.00001443& 0.00205491& 0.99008269\end{array}\right],\\ & C_2=\left[\begin{array}{rrrr}0.00022183& 0.05942943& 0.05327854& -0.99534942\\ 0.00026594& 0.00000044& 0.00000044& -0.00000002\\ 0.00029502& 0.00000254& -0.00000225& 0.00000017\\ 0.99634129& 0.00008613& -0.00002336& 0.00001053\\ -0.00030142& 0.00050357& -0.00044968& 0.00003365\\ -0.00000027& 0.00008353& 0.00009039& -0.00000579\end{array}\right],\end{array}$$

with

$$\sigma \left(A\right)=\left\{\begin{array}{l}0.638773652186517,0.847768449750652\\ 1.002501752901569,0.960047795833900\\ 0.990353602254223\end{array}\right\},$$

we had the following results:

by the RS algorithm for minimal-gain SOF (see [10])

$$\left\{\begin{array}{ll}& K^{\left(0\right)}=\left[\begin{array}{cc}{K}_1^{\left(0\right)}& K_2^{\left(0\right)}\end{array}\right],\text{where}\\ & K_1^{\left(0\right)}=\left[\begin{array}{lll}-0.05281866& 0.30558099& -0.04123125\\ -0.74370605& -0.07272045& 0.16180699\\ -0.34989799& -0.70937255& -0.03438071\\ 0.36682921& -0.55329174& -0.42930790\end{array}\right],\\ & K_2^{\left(0\right)}=\left[\begin{array}{lll}-0.07894070& -0.83106320& 0.63513665\\ -0.57049060& 0.02944824& 0.03985277\\ 0.01822942& -0.40097959& -0.32739026\\ -0.32035712& 0.23550532& 0.55239497\end{array}\right],\\ & \sigma \left(A-{\mathit{BK}}^{\left(0\right)}C\right)=\left\{\begin{array}{l}0.88529956,0.98303140,\\ 0.99890558\pm 0.00709027i,\\ 0.99895917\pm 0.00548801i,\\ 0.99648905,0.99228590\end{array}\right\},\\ & {\sigma }_{\max }^{\left(0\right)}=4.1660\cdot {10}^6,‖K^{\left(0\right)}‖=1.1052,\end{array}\right.$$

by the RS Algorithm1

$$\left\{\begin{array}{ll}& K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{cc}{K}_1^{\left(\mathit{\text{best}}\right)}& K_2^{\left(\mathit{\text{best}}\right)}\end{array}\right],\text{where}\\ & K_1^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{lll}2.60115550& 0.71670943& 1.38518242\\ -1.21472623& 7.41955425& -2.28748737\\ 0.08032866& -1.01678491& -3.50944968\\ 0.88639191& -0.92895747& 0.77107501\end{array}\right],\\ & K_2^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{lll}-0.43606038& -0.52797919& -1.91992615\\ -0.37181168& -0.49434738& -1.97765275\\ -0.43167652& 2.40298411& -0.14274690\\ -1.34628983& -2.73288946& -0.19682963\end{array}\right],\\ & \sigma \left(A-{\mathit{BK}}^{\left(\mathit{\text{best}}\right)}C\right)=\left\{\begin{array}{l}0.86726666,0.98939915,\\ 0.99600356\pm 0.01886288i,\\ 0.99670041\pm 0.00199237i,\\ 0.97791267\pm 0.01292326i,\\ \end{array}\right\},\\ & {\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=5.2247\cdot {10}^4,‖K^{\left(\mathit{\text{best}}\right)}‖=8.1964,\end{array}\right.$$

by the MC Algorithm 2

$$\left\{\begin{array}{ll}& K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{cc}{K}_1^{\left(\mathit{\text{best}}\right)}& K_2^{\left(\mathit{\text{best}}\right)}\end{array}\right],\text{where}\\ & K_1^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{lll}0.16219293& 0.22685502& -1.41842788\\ 0.10558055& 0.04299252& -1.12924145\\ -0.01638169& -0.08184317& -0.14133958\\ -0.03859594& 0.07947463& 0.12534983\end{array}\right],\\ & K_2^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{lll}-0.23992446& -0.17536269& -0.44875450\\ -0.14481130& -0.18150136& -0.27764220\\ -0.04123859& 0.04014410& 0.05631322\\ -0.06158264& 0.05048672& 0.16208075\end{array}\right],\\ & \sigma \left(A-{\mathit{BK}}^{\left(\mathit{\text{best}}\right)}C\right)=\left\{\begin{array}{l}0.53091472,0.95364142,\\ 0.98741833,0.99897830,\\ 0.99717581\pm 0.00564301i,\\ 0.95033459\pm 0.08956374i,\\ \end{array}\right\},\\ & {\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=1.1817\cdot {10}^6,‖K^{\left(\mathit{\text{best}}\right)}‖=195.3621,\end{array}\right.$$

and by RS + MC Algorithms 1 and 2:

$$\left\{\begin{array}{ll}& K^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{cc}{K}_1^{\left(\mathit{\text{best}}\right)}& K_2^{\left(\mathit{\text{best}}\right)}\end{array}\right],\text{where}\\ & K_1^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{lll}0.16101573& 0.52768736& -1.99112247\\ -1.8453102& 11.28218548& -3.45382160\\ 0.45368113& -0.97316671& -6.96467764\\ 0.57106266& -0.53606197& -0.45209883\end{array}\right],\\ & K_2^{\left(\mathit{\text{best}}\right)}=\left[\begin{array}{lll}-0.67588625& -0.18705486& -0.28941025\\ -0.69074912& -0.26655181& 0.20004806\\ -0.83280356& 3.20877857& -0.31860927\\ -3.09450298& -0.73194640& -0.07799223\end{array}\right],\\ & \sigma \left(A-{\mathit{BK}}^{\left(\mathit{\text{best}}\right)}C\right)=\left\{\begin{array}{l}0.98195937\pm 0.03551664i,\\ 0.99264852\pm 0.01953917i,\\ 0.99354901\pm 0.00357529i,\\ 0.99828371,0.98490431\end{array}\right\}\\ & {\sigma }_{\max }^{\left(\mathit{\text{best}}\right)}=3.1783\cdot {10}^4,‖K^{\left(\mathit{\text{best}}\right)}‖=12.1029.\end{array}\right.$$

6. Concluding remarks

The Ray-Shooting Method is a powerful tool, since it practically solves the problem of LQR via SOF, for real-life discrete-time LTI systems. The proposed hybrid algorithm RS + MC has good performance in terms of run-time, in terms of the quality of controllers (by reducing the starting point LQR functional value and by reducing the controller norm) and in terms of the success rate in finding a starting point feasible with respect to the needed $\alpha$ -stability. The RS + MC algorithm has a proof of convergence in probability to a global minimum (as is evidently seen from the experiments). This enlarges the practicality and scope of the Ray-Shooting Method in solving hard complexity control problems, and we expect to receive more results in this direction.