Structured real radius of controllability for higher order LTI systems

TL;DR

This paper introduces a method to compute the radius of controllability for higher order LTI systems by relating it to the low rank approximation of a Toeplitz structured matrix, enhancing robustness analysis.

Contribution

It establishes a novel link between controllability radius and structured low rank approximation, providing an effective computational approach for higher order systems.

Findings

Method effectively computes controllability radius.

Approach outperforms benchmark numerical methods.

Numerical examples validate the approach.

Abstract

In this paper, we consider the problem of computing the nearest uncontrollable (C-uncontrollable) system to a given higher order system. The distance to the nearest uncontrollable system, also termed as the radius of controllability, is a good measure of gauging the numerical robustness of the given system with respect to controllability. Here, we invoke the equivalence of C-controllability of a higher order system with full rank property of a certain Toeplitz structured matrix. This enables us to pose the problem of computing the radius of controllability as equivalent to the problem of computing the nearest structured low rank approximation of this Toeplitz structured matrix. Through several numerical examples and comparison with the benchmark numerical problem, we illustrate that our approach works well.

Tables

$\delta$	1	0.6	0.4	0.2	0.1	0.01	0
$S{R}_c$	0.3193	0.3820	0.4132	0.2000	0.1	0.01	0

Table 2. TABLE I :

$(C_1,C_2,L,R)$	$\min \{C_1,C_2,L\}$	$S{R}_c$
(1, 1, 1, 1)	1	0.9997
(2, 1.5, 3, 1)	1.5	1.4998
(2, 3.5, 1.2, 4)	1.2	1.2000
(0.0001, 0.1, 10, 3)	0.0001	0.0001
(8, 0.01, 0.1, 4)	0.01	0.0100

Table 3. TABLE II :

$(n,m)$	Avg. Time (Seconds)	Avg. No. of iterations
(5, 1)	0.0166	8.3
(7, 2)	0.0957	18.5
(10, 3)	0.6242	23.9
(15, 4)	2.4168	22.6
(20, 5)	11.1736	20.5

Table 4. TABLE III :

$\mu$	$S{R}_c$
	Using SLRA	According to [9]
0.05	0.0587	0.0587
0.1	0.1031	0.1031
0.15	0.1470	0.1470
0.2	0.1901	0.1901
0.5	0.4227	0.4227
1	0.6813	0.6811
10	0.9959	0.9959
100	1.0000	1.0000
1000	1.0000	1.0000

Equations

$$P_{\mathtt{d}}\frac{d^{\mathtt{d}}x(t)}{d{t}^{\mathtt{d}}}+P_{\mathtt{d}-1}\frac{d^{\mathtt{d}-1}x(t)}{d{t}^{\mathtt{d}-1}}+\cdots +P_1\frac{dx(t)}{dt}+P_{0}x(t)=bu(t)$$

$$E\stackrel{.}{z}(t)=Az(t)+Bu(t);E,A\in {\mathbb{R}}^{n\times n},B\in {\mathbb{R}}^{n\times m}$$

$$R_c=\underset{{\Delta }_{P_i}\in {\mathbb{R}}^{N\times N},{\Delta }_b\in {\mathbb{R}}^{N\times M}}{\min }\{\Vert \left[{\Delta }_{P_{\mathtt{d}}}\cdots {\Delta }_{P_0}{\Delta }_b\right]{\Vert }_F|$$

$$(P_{\mathtt{d}}+{\Delta }_{P_{\mathtt{d}}},\dots P_0+{\Delta }_{P_0},b+{\Delta }_b)\text{ is uncontrollable}\}$$

$$R_c=\underset{{\Delta }_E,{\Delta }_A\in {\mathbb{R}}^{n\times n},{\Delta }_B\in {\mathbb{R}}^{n\times m}}{\min }\{\Vert \left[{\Delta }_E{\Delta }_A{\Delta }_B\right]{\Vert }_F|$$

$$(E+{\Delta }_E,A+{\Delta }_A,B+{\Delta }_B)\text{ is uncontrollable}\}.$$

$$E\stackrel{.}{z}(t)=Az(t)+Bu(t),$$

$$\mathcal{R}(s)=U(s)[s{E}_r-A_r]V(s)$$

$$\mathcal{Q}(s)\equiv\left[\begin{array}[]{cccc|cccc}I_{N}&0&\dots&0&0&0&\dots&0\\ 0&I_{N}&\dots&0&0&0&\dots&0\\ &&\ddots&&&&&\\ 0&0&\dots&P(s)&0&0&\dots&b\\ \hline\cr 0&0&\dots&0&I_{M}&0&\dots&0\\ 0&0&\dots&0&0&I_{M}&\dots&0\\ &&\ddots&&&&\ddots&\\ 0&0&\dots&0&0&0&\dots&I_{M}\end{array}\right]$$

$$\left[\begin{array}[]{c|ccccc}[sE-A~{}B]&0&0&0&\dots&0\\ \hline\cr 0&-I_{M}&sI_{M}&0&\dots&0\\ 0&0&-I_{M}&sI_{M}&\dots&0\\ &&&&\ddots&\\ 0&0&0&\dots&-I_{M}&sI_{M}\end{array}\right]$$

$$[P(s)b]\text{ has full rank }\iff$$

$$\iff$$

$$r_1=$$

$$\text{rank}[s(E+{\Delta }_E)-(A+{\Delta }_A)B+{\Delta }_B]<n$$

$$\text{ for some }s\in \mathbb{C}\}$$

$$r_2=\underset{{\Delta }_E,{\Delta }_B}{\min }\{\Vert [{\Delta }_E{\Delta }_B]{\Vert }_F,$$

$$\text{ such that, }\text{rank}[(E+{\Delta }_E)(B+{\Delta }_B)]<n\}$$

$$\mathcal{C}(E,A,B)$$

$$R_c={\Vert \left[\Delta E\Delta A\Delta B\right]\Vert }_F.$$

$$\mathcal{C}(E+\Delta E,A+\Delta A,B+\Delta B)$$

$$=\left[\begin{array}{cccccccc}-\hat{A}& \hat{B}& & & & & & \\ \hat{E}& 0& -\hat{A}& \hat{B}& & & & \\ & & \hat{E}& 0& & & & \\ & & & & \ddots & & & \\ & & & & & -\hat{A}& \hat{B}& \\ & & & & & \hat{E}& 0& \hat{B}\end{array}\right]$$

$${\Vert \left[\Delta E\Delta A\Delta B\right]\Vert }_F=\underset{{\Delta }_E,{\Delta }_A,{\Delta }_B}{\min }\{{\Vert \left[{\Delta }_E{\Delta }_A{\Delta }_B\right]\Vert }_F|$$

$$\mathcal{C}(E+{\Delta }_E,A+{\Delta }_A,B+{\Delta }_B)\text{ is not full row rank}\}$$

$$S{R}_c=\underset{{\Delta }_E\in {\mathcal{S}}_E,{\Delta }_A\in {\mathcal{S}}_A,{\Delta }_B\in {\mathcal{S}}_B}{\min }\{\Vert \left[{\Delta }_E{\Delta }_A{\Delta }_B\right]{\Vert }_F|$$

$$(E+{\Delta }_E,A+{\Delta }_A,B+{\Delta }_B)\text{ is uncontrollable}\}.$$

$$S{R}_c=\underset{{\Delta }_E\in {\mathcal{S}}_E,{\Delta }_A\in {\mathcal{S}}_A,{\Delta }_B\in {\mathcal{S}}_B}{\min }\{{\Vert \left[{\Delta }_E{\Delta }_A{\Delta }_B\right]\Vert }_F|$$

$$\mathcal{C}(E+{\Delta }_E,A+{\Delta }_A,B+{\Delta }_B)\text{ is not full row rank}\}.$$

$$\underset{X\in \mathcal{X}\cap \Omega }{\min }\Vert X-X_{\ast }{\Vert }_F,$$

$$\underset{\alpha ,z}{\min }{\Vert \left[\begin{array}{c}\omega r(\alpha ,z)\\ D\alpha \end{array}\right]\Vert }_2$$

$$r(\alpha +\Delta$$

$$\approx (y+f_1)-(Y+E_1)z-(Y+E_1)\Delta z-\Delta E_{1}z$$

$$\approx r(\alpha ,z)+(f_1-\Delta E_{1}z)-(Y+E_1)\Delta z$$

$$\underset{\Delta \alpha ,\Delta z}{\min }{\Vert \left[\begin{array}{cc}\omega (S-P)& Y+E_1\\ I& 0\end{array}\right]\left[\begin{array}{c}\Delta \alpha \\ \Delta z\end{array}\right]+\left[\begin{array}{c}-\omega r\\ \alpha \end{array}\right]\Vert }_2$$

$$E\stackrel{.}{x}(t)=Ax(t)+Bu(t)$$

$$\left[\begin{array}{cccc}{C}_1& 0& 0& 0\\ 0& C_2& 0& 0\\ 0& 0& -L& 0\\ 0& 0& 0& 0\end{array}\right]\stackrel{.}{x}=\left[\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ -1& 1& 0& 0\\ 1& 0& 0& R\end{array}\right]x+\left[\begin{array}{c}0\\ 0\\ 0\\ -1\end{array}\right]u_e$$

$$M\stackrel{..}{x}(t)+K(\mu )x(t)=f(t)$$

Full text

I Introduction

Many real life phenomena can be modeled as linear time-invariant (LTI) higher order systems as defined in the following equation:

[TABLE]

Here $P_{0},P_{1},\dots,P_{{\tt d}-1},P_{\tt d}\in\mathbb{R}^{N\times N},~{}b\in\mathbb{R}^{N\times M}$. It should be noted that the case when ${\tt d}=1$, the system in (1) becomes an LTI descriptor system. In addition if we consider the leading coefficient matrix $P_{\tt d}$ to be non-singular, the system in (1) can easily be transformed into the usual state space system. However, in many real life applications, the matrix $P_{\tt d}$ turns out to be a singular matrix for ${\tt d}\geq 1$. (see for instance [1, 2]).

Controllability of the system in (1) plays a key role in the analysis of such systems. However, the notion of controllability gives only a yes/no type of answer. It may so happen that small perturbations to system parameters may render the system uncontrollable. In such a case, it is not sufficient to know only whether the given system is controllable, but it is advisable to know how far this system is from being uncontrollable. In other words, one may be interested in knowing the nearest uncontrollable system to the given system. The distance to the nearest uncontrollable system, also called the radius of controllability, is a robust measure of controllability. Several studies in the past have been conducted in the area of computation of the radius of controllability for state space systems [3, 4, 5, 6] as well as for descriptor systems and higher order systems [7, 8, 9].

In [10], the author first introduced the notion of real radius of controllability for state space systems of the form $\overset{.}{x}(t)=Ax(t)+Bu(t)$ and many researchers have followed this work as evident in the literature (see for instance [11, 12, 3, 4, 5, 6] and references therein). In [12], the radius of controllability for state space system was computed by minimizing the minimum singular value function of a matrix pencil, which is not numerically efficient. A similar concept was implemented in [4, 3] to compute the radius of controllability. In [6], the problem of computing the radius of controllability for state space system was shown to be equivalent to a problem of computing the nearest Structured Low Rank Approximation (SLRA) of a Toeplitz structured matrix.

In [9], a higher order linear time invariant system was considered. The radius of controllability was computed using generalized singular value decomposition of a matrix pair, though finding generalized singular value decomposition may not be numerically efficient in general. Basically, a generalized real perturbation value function is minimized in the complex plane, which is an one dimensional optimization problem. The structured radius of controllability of higher order system was also discussed in [7].

In this paper we will illustrate an approach to compute the radius of controllability which involves computation of the nearest SLRA of a structured matrix. Due to some physical constraints, it may happen that some of the entries of $P_{i}(i=0,1,\dots,{\tt d})$ and $b$ are not allowed to be perturbed while computing the radius of controllability. In such cases, the proposed approach is shown to be easily applied to compute the nearest uncontrollable system. The proposed technique is computationally efficient as illustrated subsequently and avoids the computation of generalized singular value decomposition as done in [9].

Criterion for Controllability of Higher order Systems

Here, we state various types of controllability concepts introduced in the literature for LTI higher order, descriptor systems. For a comprehensive study of the descriptor systems the reader is referred to [13, 1, 14]. Controllability for linear descriptor systems is analyzed by some graph-theoretic approach in [15]. The controllability of the system in (1) is discussed in the following result from [7].

Theorem I.1

The Higher order system in (1) is $C^{\tt d}$-controllable if and only if

rank$([P(s)~{}b])=N,$* for all finite $s\in\mathbb{C}$, and* 2. 2.

rank$([P_{\tt d}~{}b])=N$**

where, $P(s):=P_{\tt d}s^{\tt d}+\dots+P_{1}s+P_{0}$.

Descriptor systems

Descriptor system is a particular case of higher order system with ${\tt d}=1$. Consider the following descriptor system

[TABLE]

Controllability of the descriptor system in (2) plays an important role in our paper. Now we state the theorem regarding the criterion of controllability of descriptor systems, which directly follows from Theorem I.1.

Theorem I.2

The descriptor system in (2) is C-controllable if and only if

rank$[sE-A~{}~{}B]=n$, for all finite $s\in\mathbb{C}$, and 2. 2.

rank$[E~{}~{}B]=n$.

II Radius of controllability

As discussed in the previous section, we now formally define the real radius of controllability, denoted as $R_{c}$, for the system in (1). In order to simplify the notation we write $(P_{\tt d},\dots,P_{1},P_{0},b)$ is either C-controllable/uncontrollable to mean the system in (1) is C-controllable/uncontrollable respectively.

[TABLE]

Before we proceed further, there are several remarks in order.

Remark II.1

The perturbations $\Delta P_{i}\in\mathbb{R}^{N\times N},~{}\Delta b\in\mathbb{R}^{N\times M}$ are allowed to be only real and hence the name real radius of controllability. The case where complex perturbations are allowed is not considered in this paper, though the theory proposed in this paper can be extended very easily to accommodate for the complex case. **

Remark II.2

The norm of the perturbation matrices considered here is the Frobenius norm as opposed to the spectral norm as in previous works. The advantage of Frobenius norm is that one can quantify the change in individual entries of $(P_{\tt d},\dots,P_{1},P_{0},b)$ which has great physical advantage. **

Remark II.3

From the definition of real radius of controllability, it is clear that there is no restriction on the structure of the perturbations. In many cases, such structure restriction on the perturbation matrices can arise naturally due to structure constraints of the systems. In such a case, when the perturbations $\Delta P_{i}~{}(i=0,1,\dots,{\tt d}),~{}\Delta b$ are structured perturbations, the radius of controllability, for obvious reasons, is called as structured real radius of controllability. We note here, that the theory that we develop to compute radius of controllability can be used to obtain the structured radius of controllability as well by slightly modifying certain equations. This will be illustrated by several numerical examples. A particular case of interest that perturbation to the leading coefficient matrix $P_{\tt d}$ is not permitted is explored in the numerical examples presented later. **

In this paper, we transform the higher order system into a descriptor system, then we compute the radius of controllability for the descriptor system only. The radius of controllability for descriptor system in (2) can be defined in similar way.

[TABLE]

The remainder of the paper is organized as follows: in the next section, we show the equivalence of computing the radius of controllability with the computation of the nearest SLRA of a certain Toeplitz matrix. In Section IV, we discuss an efficient algorithm (Structured Total Least Norm) to solve the SLRA problem. We illustrate the advantages of our approach by demonstrating several numerical case studies in Section V before we conclude in Section VI.

III Mathematical Formulation

In this section, we reformulate the problem of computing the radius of controllability of a linear time invariant higher order system as the problem of computing the nearest SLRA of a certain Toeplitz matrix. We first transform the higher order LTI system (1) into its canonical form. The canonical form of equation (1) can be written as a descriptor system

[TABLE]

where $E=\begin{bmatrix}P_{\tt d}&0&\dots&0\\ 0&I_{N}&\dots&0\\ &&\ddots&\\ 0&0&\dots&I_{N}\end{bmatrix}\in\mathbb{R}^{N{\tt d}\times N{\tt d}}$,

$A=\begin{bmatrix}-P_{{\tt d}-1}&-P_{{\tt d}-2}&\dots&-P_{1}&-P_{0}\\ I_{N}&0&\dots&0&0\\ 0&I_{N}&\dots&0&0\\ &&\ddots&&\\ 0&0&\dots&I_{N}&0\end{bmatrix}\in\mathbb{R}^{N{\tt d}\times N{\tt d}}$, $B=\begin{bmatrix}b^{T}&0^{T}&\cdots&0^{T}\end{bmatrix}\in\mathbb{R}^{N{\tt d}\times M}$

Therefore System (5) is same as System (2) with $n=N{\tt d}$ and $m=M$. Now we will state and prove a theorem which connects the controllability of the higher order system in (1) and its canonical form in (5). Before proving the theorem, we state the following necessary lemma.

Lemma III.1

[16*,  6.3-20]**

Consider the $\tt r^{\rm{th}}$ order polynomial matrix $R(s)=R_{\tt r}s^{\tt r}+\dots+R_{1}s+R_{0}$, where $R_{i}\in\mathbb{R}^{r_{1}\times r_{2}}$ for $i=0,1,\dots,\tt r$. Then*

[TABLE]

*where $\mathcal{R}(s)=\text{block diag}\{\underbrace{I_{r_{2}},\dots,I_{r_{2}},R(s)}_{\tt r\text{-terms}}\},$

$E_{r}=\begin{bmatrix}R_{\tt r}&0&\dots&0\\ 0&I_{r_{2}}&\dots&0\\ &&\ddots&\\ 0&0&\dots&I_{r_{2}}\end{bmatrix}$,

$A_{r}=\begin{bmatrix}-R_{\tt r-1}&-R_{\tt r-2}&\dots&-R_{1}&-R_{0}\\ I_{r_{2}}&0&\dots&0&0\\ 0&I_{r_{2}}&\dots&0&0\\ &&\ddots&&\\ 0&0&\dots&I_{r_{2}}&0\end{bmatrix}$ and $U(s),V(s)$ are unimodular111A matrix $U(s)$ is called unimodular if Det$[U(s)]$ is a non-zero constant. matrices.*

Theorem III.1

The higher order system in (1) is $C^{\tt d}$-controllable if and only if its canonical form in (5) is C-controllable.

Proof: It is observed that if we can prove the equivalence of the conditions of Theorem I.1 and Theorem I.2, then the proof is done. Now look at the structures of the matrices $E,A$ and $B$ shown in (5). From the structure of $[E~{}B]$, it is observed that $[E~{}B]$ has full rank if and only if $[P_{\tt d}~{}b]$ has full rank. Therefore condition (2) of Theorem I.1 and Theorem I.2 are equivalent. Now it is sufficient to prove only that rank$([P(s)~{}b])=N$ if and only if rank$([sE-A~{}B])=n$.

Let $Q(s)=[P(s)~{}b]\in\mathbb{R}^{N\times(N+M)}$. It can be written that $Q(s)=Q_{\tt d}s^{\tt d}+\dots+Q_{1}s+Q_{0}$, where $Q_{0}=[P_{0}~{}b],~{}Q_{i}=[P_{i}~{}0]~{}\text{for}~{}i=1,2,\dots,{\tt d}.$ Now from Lemma III.1, we get there exist unimodular matrices $U(s)$ and $V(s)$ such that $\mathcal{Q}(s)=U(s)[sE_{1}-A_{1}]V(s)$ where $\mathcal{Q}(s)=\text{block diag}\{\underbrace{I_{(N+M)},\dots,I_{(N+M)},Q(s)}_{{\tt d}\text{-terms}}\},$

$E_{1}=\textrm{diag}\begin{bmatrix}Q_{\tt d}&&I_{(N+M)}&\cdots&I_{(N+M)}\end{bmatrix}$, and

$A_{1}=\begin{bmatrix}-Q_{{\tt d}-1}&-Q_{{\tt d}-2}&\dots&-Q_{1}&-Q_{0}\\ I_{(N+M)}&0&\dots&0&0\\ 0&I_{(N+M)}&\dots&0&0\\ &&\ddots&&\\ 0&0&\dots&I_{(N+M)}&0\end{bmatrix}$

Therefore $\mathcal{Q}(s)$ has full rank if and only if $[sE_{1}-A_{1}]$ has full rank for all finite $s\in\mathbb{C}$. Now substitute $Q(s)=[P(s)~{}b]$ in $\mathcal{Q}(s)$, then after some elementary row/column operations we get

[TABLE]

Now from this structure it is clear that $\mathcal{Q}(s)$ has full rank if and only if $[P(s)~{}b]$ has full rank for all finite $s\in\mathbb{C}$. Again substitute $Q_{0}=[P_{0}~{}b],~{}Q_{i}=[P_{i}~{}0]~{}\text{for}~{}i=1,2,\dots,{\tt d}$ in $sE_{1}-A_{1}$. Similarly after some elementary row/column operations we get $sE_{1}-A_{1}\equiv$

[TABLE]

From the structure of the matrix, it is also clear that $sE_{1}-A_{1}$ has full rank if and only if $[sE-A~{}B]$ has full rank for all finite $s\in\mathbb{C}.$ Therefore combining the results we get for all finite $s\in\mathbb{C}$,

[TABLE]

Hence the result is proved. $\square$

Remark III.1

From Theorem III.1, it directly follows that nearest uncontrollable system to the higher order system (1) can be obtained from the nearest uncontrollable system to the descriptor system (5). Suppose $(\tilde{E},\tilde{A},\tilde{B})$ is the nearest uncontrollable system to System (5) satisfying the additional constraint that $(\tilde{E},\tilde{A},\tilde{B})$ and $(E,A,B)$ have the same structures. Therefore from (5) it can be observed that the perturbations are allowed only in the blocks $(P_{\tt d},\dots,P_{1},P_{0},b)$ in the perturbed $(E,A,B)$. Suppose the perturbed matrices $(\tilde{E},\tilde{A},\tilde{B})$ follows the perturbed matrices $(\tilde{P}_{\tt d},\dots,\tilde{P}_{1},\tilde{P}_{0},\tilde{b})$. Then we can conclude that the corresponding perturbed higher order system $(\tilde{P}_{\tt d},\dots,\tilde{P}_{1},\tilde{P}_{0},\tilde{b})$ is the nearest uncontrollable system to (1). Therefore, in order to compute the real radius of controllability of higher order system, it is sufficient to compute the structured real radius of controllability (the radius of controllability obtained by respecting the structure of perturbation in the parameter matrices) of the equivalent descriptor system. **

In light of condition (2) of Theorem I.2, the problem of computing the real radius of controllability of System (2) can be split into two following subproblems:

Problem III.1

[TABLE]

Problem III.2

[TABLE]

Therefore it can be observed that radius of controllability of the system (2) can be found by choosing $R_{c}=\min\{r_{1},r_{2}\}$ with the understanding that the system will be uncontrollable if any of the matrices among $[sE-A~{}~{}B]$ and $[E~{}~{}B]$ loses its rank. In [9], these two problems were solved individually to find $r_{1}$ and $r_{2}$. We now present a result which states that controllability of System (2) can be shown to be equivalent to the full rank property of a certain Toeplitz structured matrix.

Theorem III.2

[13, Theorem 2-2.1]* Consider the descriptor system given in (2). Then the following statements are equivalent.*

*The system in (2) is C-controllable. * 2. 2.

$\mathcal{C}(E,A,B)\in\mathbb{R}^{n^{2}\times n(n+m-1)}$* is full row-rank matrix, where,*

[TABLE]

The equivalence of C-controllability in Statement 2 of Theorem III.2 is the key to formulate the computation of the real radius of controllability as a structured low rank approximation problem. From Theorem III.2, we immediately conclude the following theorem which is the main result of the paper.

Theorem III.3

Let an LTI descriptor system be given as in (2). Then computing the radius of controllability, denoted by $R_{c}$, of the system in (2), is equivalent to computing nearest SLRA $~{}\mathcal{C}(\hat{E},\hat{A},\hat{B})$ to $\mathcal{C}(E,A,B)$, introduced in Statement 2 of Theorem III.2. In fact,

[TABLE]

where, $\hat{E}=E+\Delta E,~{}\hat{A}=A+\Delta A,~{}\hat{B}=B+\Delta B$.

Proof: Suppose that the given system $(E,A,B)$ is C-controllable, then the statement 2 states that the Toeplitz structured matrix $\mathcal{C}(E,A,B)\in\mathbb{R}^{n^{2}\times n(n+m-1)}$ is a full row rank matrix. In order to compute the real radius of controllability, we seek the perturbations $\Delta E,\Delta A,\Delta B$ of appropriate dimensions such that the perturbed system $(E+\Delta E,A+\Delta A,B+\Delta B)$ is uncontrollable. Then Statement 2 in Theorem III.2 would imply that the Toeplitz structured matrix associated with the perturbed system, denoted as $\mathcal{C}(E+\Delta E,A+\Delta A,B+\Delta B)$ as in (8), is not a full row rank matrix.

[TABLE]

where $\hat{A}=A+\Delta A$, $\hat{E}=E+\Delta E$ and $\hat{B}=B+\Delta B$. Clearly, $\mathcal{C}(E+\Delta E,A+\Delta A,B+\Delta B)$ is a structured perturbation of $\mathcal{C}(E,A,B)$.

Conversely, given a system $(E,A,B)$, construct the Toeplitz structured matrix as in (6). If we are able to compute some perturbations $(\Delta E,\Delta A,\Delta B)$ of appropriate sizes such that the perturbed Toeplitz matrix in (8) is not full row rank, then the corresponding system $(E+\Delta E,A+\Delta A,B+\Delta B)$ is clearly uncontrollable. Let these perturbations, in particular, have the following property.

[TABLE]

In such a case, the matrix $\mathcal{C}(E+\Delta E,A+\Delta A,B+\Delta B)$ is the nearest structured low rank approximation to the matrix $\mathcal{C}(E,A,B)$. Clearly, from (9) both these matrices have the same structure. Further, not only the matrix $\mathcal{C}(E+\Delta E,A+\Delta A,B+\Delta B)$ is of lower rank than $\mathcal{C}(E,A,B)$, but is also the nearest with the low rank property. This justifies why this perturbed matrix is called as the nearest Structured Low Rank Approximation (SLRA) of the given matrix. Thus, from (4) and (9), it is clear that the problem of computing the real radius of controllability of the system in (2) is equivalent to computation of the nearest SLRA $~{}\mathcal{C}(\hat{E},\hat{A},\hat{B})$ of $\mathcal{C}(E,A,B)$. Also clearly $R_{c}=\left\|\left[\Delta E~{}~{}\Delta A~{}~{}\Delta B\right]\right\|_{F},$ where, $\hat{E}=E+\Delta E,~{}\hat{A}=A+\Delta A,~{}\hat{B}=B+\Delta B$. $\square$

Remark III.2

Note that we are using Statement 2 in Theorem III.2 to compute the radius of controllability. This enables us to compute the radius of controllability directly as opposed to the other approaches in the literature where $r_{1}$ and $r_{2}$ are computed as in Problem III.1 and Problem III.2 and the minimum of these is declared to be $R_{c}$. **

Remark III.3

From the definition and Theorem III.3, we can understand the radius of controllability in another way. Theorem III.3, in fact, says that there a ball around the system $(E,A,B)$ with radius equal to $R_{c}$ such that every system inside this ball is C-controllable. The larger radius of this ball ensures higher robustness with respect to controllability. On the other hand a small value of $R_{c}$ indicates that an uncontrollable system is nearby indicating the system may be susceptible to small perturbations. **

In the literature, the problem of computing the nearest SLRA to a given matrix is studied extensively. We will explain the problem of SLRA formally and a procedure to solve this problem in the next section.

As pointed out in Remark II.3 and Remark III.1, it may so happen that certain physical constraints may impose structural constraints on the perturbations while computing the nearest uncontrollable system to (2). We now define Structured Real Radius of Controllability for descriptor systems under the assumption that these structure constraints on the perturbations are linear or affine. Let $\mathcal{S}_{E}\subseteq\mathbb{R}^{n\times n}$ be a linear subspace (or affine subset) which characterizes the perturbations allowed in the system matrix $E$. Similarly, let $\mathcal{S}_{A}\subseteq\mathbb{R}^{n\times n}$ and $\mathcal{S}_{B}\subseteq\mathbb{R}^{n\times m}$ be the linear subspaces (or affine subset) which characterize the perturbations to be given to the system matrices $A$ and $B$ respectively. Then the Structured Real Radius of Controllability is defined as follows.

[TABLE]

It should be clear that by performing the same calculations as in the case of real radius of controllability, we can show that the problem of computing the structured real radius of controllability can be posed as an SLRA problem. For the sake of completeness we state the following theorem.

Theorem III.4

Let an LTI descriptor system be given as in (2). Let $SR_{c}$ denote the structured real radius of controllability of the system in (2). Then $SR_{c}$ can be computed by computing the nearest SLRA $\mathcal{C}(E+\Delta E,A+\Delta A,B+\Delta B)$ to $\mathcal{C}(E,A,B)$ by solving the following optimization problem.

[TABLE]

Proof: In the Theorem III.3, the nearest SLRA was computed to $\mathcal{C}(E,A,B)$ for any $E,A\in\mathbb{R}^{n\times n},~{}B\in\mathbb{R}^{n\times m}$, but here the same thing has to be done for the structured matrices $E\in\mathcal{S}_{E},~{}A\in\mathcal{S}_{A},~{}B\in\mathcal{S}_{B}$. Therefore, the proof is similar as the proof of Theorem III.3, we omit this. $\square$

We now give the algorithm to compute the radius of controllability for a given higher order system.

Algorithm** III.1**

Structured real radius of controllability

Input:* Higher order system $(P_{\tt d},\dots,P_{1},P_{0},b)$

Output: The structured real radius of controllability, $SR_{c}$

Step 1: Construct structured matrices $(E,A,B)$ given in Eq.(5)

Step 2: Construct the Toeplitz structured matrix $\mathcal{C}(E,A,B)$ from the triple $(E,A,B)$ as in (6).

Step 3: Compute nearest SLRA $\mathcal{C}(E+\Delta E,A+\Delta A,B+\Delta B)$ to $\mathcal{C}(E,A,B)$ as in (9) (See Algorithm 4.1).

Step 3: Compute $SR_{c}=\left\|\left[\Delta E~{}~{}\Delta A~{}~{}\Delta B\right]\right\|_{F}$.*

IV Structured Low Rank Approximation (SLRA)

In this section, we state the problem of computing the nearest SLRA to a given matrix and explain the Structured Total Least Norm (STLN) algorithm to compute the nearest SLRA to the given matrix. Let $X_{*}\in\mathbb{R}^{M\times N}$ be a given matrix with a particular affine structure. We want to compute the nearest low rank approximation $X$ which inherits the same structure as $X_{*}$. Therefore, the SLRA problem is an optimization problem where the distance between $X_{*}$ and $X$ has to be minimized. Mathematically, the SLRA problem can be formulated as

[TABLE]

where,

$\begin{aligned} &\Omega\in\mathbb{R}^{M\times N}\text{ is a space (or set) of particular linear (or affine)}\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{ structured matrices},\\ &\mathcal{X}:=\left\{X\in\mathbb{R}^{M\times N}~{}~{}|~{}~{}\text{rank}(X)<\text{rank}(X_{*})\right\}.\end{aligned}$ In literature, many techniques are reported to solve SLRA problem (see for instance [17, 18, 19] and the references therein). In this paper, we will use Structured Total Least Norm (STLN) algorithm to compute the nearest SLRA [17].

IV-A Structured total least norm (STLN)

Suppose the given matrix $X_{*}$ is partitioned like $X_{*}=[Y~{}~{}y]$, where $y$ is the last column of $X_{*}$. The target is to solve the system $Yz\approx y$ approximately. In other words, we aim to find the minimum perturbations in $Y$ and $y$ in such a way that $(Y+E_{1})z=(y+f_{1}),$ where $[E_{1}~{}~{}f_{1}]$ belongs to the linear structure of the space $\Omega$. Let $\{\Phi_{1},\Phi_{2},\ldots\Phi_{\ell}\}$ be an ordered basis of the linear space of $\Omega$ and $\alpha\in\mathbb{R}^{\ell}$ be the co-ordinate vector of $X_{*}$ with respect to this basis. Then the STLN problem can be stated as:

[TABLE]

where $r(\alpha,z)=y-(Y+E_{1})z$ is the residual of the system, $\omega$ is any large weight and $D$ is a suitably chosen weight matrix. In our case, we assume $D=I$. Let a small change in $\alpha$ be denoted as $\Delta\alpha$ and a small change in $z$ be $\Delta z$, and a small perturbation in the error matrix $E_{1}$ be $\Delta E_{1}$. Then the perturbed residual becomes

[TABLE]

Compute a matrix $P$ such that $f_{1}=P\Delta\alpha$, and also compute a structured matrix $S$ dependent on $z$ which satisfies $S\Delta\alpha=\Delta E_{1}z$. Formation of the matrix $S$ is explained in [20]. Then $r(\alpha+\Delta\alpha,z+\Delta z)$ becomes $r(\alpha,z)-(S-P)\Delta\alpha-(A_{1}+E_{1})\Delta z$. In this notation, the STLN problem is reformulated as

[TABLE]

Therefore the Algorithm for STLN can be written as:

Algorithm** IV.2**

STLN

Input: Matrices $Y,~{}y,~{}\text{error tolerance }\epsilon>0$

Output: Error matrix $[E_{1}~{}~{}f_{1}]$, with least norm.

Step 1: Set $\omega$ to be a large number, $E_{1}=0,f_{1}=0$ and find $z$ from $\min_{z}\|y-Yz\|$, and construct $S$ from $z$.

Step 2: Set $r=y-Yz$

Step 3: Repeat

           a. Solve the minimization problem in (13).

           b. Set $z:=z+\Delta z,\alpha:=\alpha+\Delta\alpha$

           c. Construct $[E_{1}~{}~{}f_{1}]$ from $\alpha$ and $S$ from $z$

           d. $r=(y+f_{1})-(Y+E_{1})z$

until ($\|\alpha\|,\|z\|<\epsilon$)*

Remark IV.1

In every iteration of Algorithm IV.2, the optimization problem in (13) is solved. Note that, this optimization problem is a least square problem. There are many efficient methods like $QR$ decomposition to solve the least square problem efficiently. Further, it is clear from (13) that matrices involved have a particular structure of zeros which can be exploited to solve this least square problem much more efficiently. **

Remark IV.2

It should be noticed that the optimization problem in SLRA formulation has many local solutions (as explained in [21]). STLN algorithm does not guarantee the global solution to this optimization problem. Also, STLN may be very sensitive with respect to the partition of the matrix we are considering at the beginning of the formulation. In the formulation, we have partitioned the matrix by its last column as $X_{*}=[Y~{}~{}y]$, where $y$ is the last column. However, instead of chosing the last column, any other column of the matrix $X_{*}$ can be chosen at the time of partition. Though theoretically chosing any partition of the matrix $X_{*}$ is equivalent and should lead to the same result, in practice it is observed that this affects the numerical results. Therefore care must be exercised while choosing the partition. Another important factor in STLN is the choice of the constant $\omega$. This choice may also affect the rate of convergence as well the local minimum to which STLN converges. **

Remark IV.3

As stated earlier, we now explain a way in which we can compute the structured radius of controllability. Observe that the structure constraint imposed on the perturbation matrices can be characterized by a linear subspace. In the STLN formulation, we then chose a basis for this linear structure space accordingly. Once this is done, then the rest of the calculation remains same. This shows that, in this formulation, any linear structure restriction on perturbations can be very easily incorporated. This is a great advantage of the proposed formulation. **

V Numerical Examples

In this section, we study some numerical examples to compute the structured radius of controllability of higher order, descriptor systems.

Example V.1

This is an example taken from the paper [9]. The descriptor system is given by

[TABLE]

*where, the system matrices are given by

$E=\begin{bmatrix}1.8&0&0\\ 0&0.34&0\\ 0&0&0\\ \end{bmatrix},A=\begin{bmatrix}2&-0.91&-0.088\\ 0.19&0.25&0.51\\ 0.64&0.31&-0.59\end{bmatrix},\\ B=\begin{bmatrix}-0.63\\ 0.53\\ -0.58\end{bmatrix}.$ The system (14) is controllable, our interest is to compute the nearest uncontrollable system. In this case it is assumed that no perturbation in the matrix $E$ is allowed i.e., $\Delta E=0$. According to the theory, STLN technique (seting, $\omega=10^{13}$) is used to compute the SLRA of the matrix $\mathcal{C}(E,A,B)\in\mathbb{R}^{9\times 9}$, constructed from the given $(E,A,B)$. In this case, we get the structured radius of controllability $SR_{c}=0.3436$, (we have taken 2-norm of the perturbation matrix in order to facilitate the comparison with the results in [9]) and the corresponding nearest uncontrollable system is $\hat{E}\overset{.}{x}(t)=\hat{A}x(t)+\hat{B}u(t)$, where

$\hat{A}=A+\Delta A=\begin{bmatrix}2.0000&-0.9100&-0.3333\\ 0.1900&0.2500&0.4362\\ 0.6400&0.3100&-0.4887\end{bmatrix}\\ \hat{B}=B+\Delta B=\begin{bmatrix}-0.4471\\ 0.5850\\ -0.6555\end{bmatrix}$. It can be observed that the result exactly matches with the result shown in [9].*

Example V.2

*A parametric family of descriptor systems is taken as $E\overset{.}{x}=Ax+Bu$, where the system matrices are given as

$E=\begin{bmatrix}0&2.1&0\\ 1&0&0\\ 0&0&0\end{bmatrix}A=\begin{bmatrix}1&3&0\\ 2&1&1\\ 3&1&5\end{bmatrix}B=\begin{bmatrix}1\\ 0\\ \delta\end{bmatrix},\delta\in\mathbb{R}$

We have computed the radius of controllability for different values of $\delta$ by our approach, considering no perturbation in $E$. According to the structure of the matrix $[E~{}~{}B]=\begin{bmatrix}0&2.1&0&1\\ 1&0&0&0\\ 0&0&0&\delta\end{bmatrix}$, it is clear that rank$[E~{}~{}B]$ is determined by the parameter $\delta$. When the value of $\delta$ is very close to zero, the matrix $[E~{}~{}B]$ is going to loose its rank, equivalently the system is going to become uncontrollable. Therefore according to our approach, using the matrices $E,A,B$ we constructed the matrix $\mathcal{C}(E,A,B)\in\mathbb{R}^{9\times 9}$, and computed the SLRA using STLN technique. The final results are shown in the following table.*

Example V.3

The following example, based on an electric circuit system, is adopted from [13]. Consider the circuit given in the Figure 1.

Choose the state variable as $x=[u_{C_{1}}~{}~{}u_{C_{2}}~{}~{}I_{2}~{}~{}I_{1}]^{T}$, where $u_{C_{1}},u_{C_{2}},I_{2},I_{1}$ are the voltages of capacitors: $C_{1},~{}C_{2}$ and the currents flowing over them. According to Kirchoff’s laws, the following linear descriptor system can be obtained:

[TABLE]

For this system, we have, rank$[sE-A~{}~{}B]=\text{rank}\begin{bmatrix}sC_{1}&0&0&-1&0\\ 0&sC_{2}&-1&0&0\\ 1&-1&-sL&0&0\\ -1&0&0&-R&-1\end{bmatrix}=4$, $\forall s\in\mathbb{C}$. Also rank$[E~{}~{}B]$=rank$\begin{bmatrix}C_{1}&0&0&0&0\\ 0&C_{2}&0&0&0\\ 0&0&-L&0&0\\ 0&0&0&0&-1\end{bmatrix}=4$ for any non-zero values of $C_{1},C_{2}$ and $L$. Therefore the system (15) is controllable for any values of $R$ and for any non-zero values of $C_{1},C_{2}$ and $L$. It is clear from the structure of the matrices that $[sE-A~{}~{}B]$ can not lose rank for any values of $C_{1},C_{2},L,R$, however the only way for the matrix $[E~{}~{}B]$ to lose its rank is by setting any of the entries $C_{1},C_{2},L$ to zero. Therefore, the nearest uncontrollable system can be obtained by setting the smallest entry of $C_{1},C_{2},L$ to zero. The radius of controllability is then $\min\{C_{1},C_{2},L\}$.

We can verify the result by computing the structured radius of controllability $SR_{c}$ in the approach we proposed here. Using the matrices $E,A,B$ we obtain $\mathcal{C}(E,A,B)\in\mathbb{R}^{16\times 16}$, whose nearest SLRA needs to be computed. However it should be remembered that the perturbations are permitted in $C_{1},C_{2},L$ and $R$ only, the other elements should remain unchanged. Therefore, this is the case of finding structured radius of controllability, where not all the entries system matrices are permitted to be perturbed. For different values of $C_{1},C_{2},L$ and $R$, our results are given in Table I, from which the accuracy and efficiency of our technique is proved.

Example V.4

In this example, we have considered ten randomly generated linear descriptor systems $(E,A,B)$ of different dimensions. In each of these cases, we have computed the structred real radius of controllability where we do not allow perturbations in $E$. The purpose of this example is to show the computational efforts required to compute the radius of controllability in terms of time and number of iterations required in STLN algorithm. We note the average number of iterations in STLN (taking $\omega=10^{8}$, $\epsilon=10^{-3}$) and average time required to compute the nearest SLRA. These results are summarized in Table II.

It should also be noted that in the case of multi-input descriptor systems, the matrix $C(E,A,B)$ is a wide matrix. Therefore, we use transpose of the $C(E,A,B)$ and implement STLN as discussed.

Example V.5

(A quadratic brake model) In [9, 22], the vibrations of a drum brake system are studied. The model consists of a rotating disc with mass $m$, radius $r$, and two contact points subject to friction and is described by the quadratic equation

[TABLE]

where the mass matrix $M$ is given by

[TABLE]

and the stiffness matrix $K(\mu)$ is

[TABLE]

$\mu$* is the friction coefficient, and $k$ and $\gamma$ are other parameters specific to the model. As in [9, 22], suppose the force on the brake system has just the vertical component determined by the input*

[TABLE]

System (16) is a second order system, which can be transformed into its first order canonical form

[TABLE]

where $B=\begin{bmatrix}0\\ 1\end{bmatrix}$. For the parameters $m=5,~{}k=1,~{}\gamma=\frac{\pi}{100}$, we computed structured radius of controllability for different values of $\mu$. It is assumed that the mass matrix $M$ is not supposed to get perturbed, only perturbations are allowed in the matrices $K$ and $B$. The results are compared with the corresponding results from [9] (see Table III). It is noticed that the results are exactly same for most of the cases, which proves the efficiency of our algorithm.

VI Concluding Remarks

In this paper, we considered the problem of computing the structured radius of controllability for higher order systems. It was shown that this problem is equivalent to the problem of computing the nearest SLRA of a certain Toeplitz structured matrix obtained from the system matrices. It is also shown that the structure constraint on the perturbations can be inherently modeled in the SLRA formulation making the computation of structured radius of controllability on par with the radius of controllability. Several numerical case studies were presented to highlight advantages of the proposed algorithm.