Sparsity-Promoting Optimal Control of Cyber-Physical Systems over Shared Communication Networks

TL;DR

This paper develops sparse control strategies for cyber-physical systems sharing communication networks, addressing delays and fairness, and introduces algorithms that optimize performance while managing bandwidth and sparsity trade-offs.

Contribution

It introduces novel sparse H2 control algorithms that incorporate network delays and fairness in shared communication environments, with proven convergence and practical validation.

Findings

Algorithms effectively balance sparsity and performance.

Delay can be a function of sparsity, affecting control trade-offs.

Fair bandwidth sharing improves multi-user control performance.

Abstract

Recent years have seen several new directions in the design of sparse control of cyber-physical systems (CPSs) driven by the objective of reducing communication cost. One common assumption made in these designs is that the communication happens over a dedicated network. For many practical applications, however, communication must occur over shared networks, leading to two critical design challenges, namely - time-delays in the feedback and fair sharing of bandwidth among users. In this paper, we present a set of sparse H2 control designs under these two design constraints. An important aspect of our design is that the delay itself can be a function of sparsity, which leads to an interesting pattern of trade-offs in the H2 performance. We present three distinct algorithms. The first algorithm preconditions the assignable bandwidth to the network and produces an initial guess for a stabilizing controller. This is followed by our second algorithm, which sparsifies this controller while simultaneously adapting the feedback delay and optimizing the H2 performance using alternating directions method of multipliers (ADMM). The third algorithm extends this approach to a multiple user scenario where optimal number of communication links, whose total sum is fixed, is distributed fairly among users by minimizing the variance of their H2 performances. The problem is cast as a difference-of-convex (DC) program with mixed-integer linear program (MILP) constraints. We provide theorems to prove convergence of these algorithms, followed by validation through numerical simulations.

Tables

Table 1. Table 1: Optimization variables for MILP reformulation of 𝓞 4 subscript 𝓞 4 \boldsymbol{\mathcal{O}}_{4}

${\delta }_{i,j}\in \{0,1\}$,	${\delta }_{i,j}=1\iff s_i\ge s_{i,j+1}$
$j\in {\mathbb{N}}_{p_i-2}$
$z_{i,1}\in ℝ$	$=\{\begin{array}{cc}{a}_{i,2}{s}_i+b_{i,2},\text{if}{\delta }_{i,1}=1& \\ a_{i,1}{s}_i+b_{i,1},\text{else}& \end{array}$
$z_{i,j}\in ℝ$,	$=\{\begin{array}{cc}\left(a_{i,j+1}-a_{i,j}\right)s_i+\left(b_{i,j+1}-b_{i,j}\right),\text{if}{\delta }_{i,j}=1& \\ 0,\text{otherwise}.& \end{array}$
$j\in {\mathbb{N}}_{p_i-2}\backslash 1$
${\tilde{z}}_i\in ℝ$	$={\sum }_{j=1}^{pi-2}{z}_{i,j}$
$d_i^{max}\in {\mathbb{Z}}_{+}$	$=s_{i,p_i}-s_{i,1}$
$d_i^{min}\in {\mathbb{Z}}_{+}$	$=s_{i,2}-s_{i,p_i}$
$t_{i,j}\in ℝ$, $j=1,2$	$\underset{s_i\in {ℛ}_i}{\text{max}}{a}_{i,j}{s}_i-b_{i,j}$
$t_{i,j}\in ℝ$, $j=3,\dots ,p_i-1$	$\underset{s_i\in {ℛ}_i}{\text{max}}\left(a_{i,j}-a_{i,j-1}\right)s_i-\left(b_{i,j}-b_{i,j-1}\right)$

Table 2. Table 2: Reformulating the constraints of 𝓞 4 subscript 𝓞 4 \boldsymbol{\mathcal{O}}_{4}

$j=1,\dots ,p_i-1$	$d_i^{max}{\delta }_{i,j}-s_i\le -s_{i,j+1}+d_i^{max}$
	$\left(d_i^{min}-ϵ\right){\delta }_{i,j}+s_i\le s_{i,j+1}-ϵ$
$j=2,\dots ,p_i-1$	${\delta }_{i,j}-{\delta }_{i,l}\le 0$,
$j=3,\dots ,p_i$	$\mp t_{i,j}{\delta }_{i,j-1}\pm z_{i,j-1}\le 0,$
	$\mp t_{i,j}{\delta }_{i,j-1}\pm z_{i,j-1}\mp \left(a_{i,j}-a_{i,j-1}\right)s_i\mp \left(b_{i,j}\mp b_{i,j-1}\right)\le -t_{i,j}$
$j=3,\dots ,p_i$	$\pm \left(t_{i,2}-t_{i,1}\right){\delta }_{i,1}\mp z_{i,1}\pm a_{i,2}{s}_i\le \pm \left(t_{i,2}-t_{i,1}\right)\mp b_{i,2}$
$j=3,\dots ,p_i$	$\pm \left(t_{i,2}-t_{i,1}\right){\delta }_{i,1}\mp z_{i,1}\pm a_{i,1}{s}_i\le \mp b_{i,1}$
$j=3,\dots ,p_i$	$\pm {\tilde{z}}_i\mp {\sum }_{j=1}^{{\overline{p}}_i}{z}_{i,j}\le 0$

Equations

$$\dot{x}(t)=Ax(t)+Bu(t)+B_{w}w(t),$$

$$\dot{x}(t)=Ax(t)-BKx(t-\tau )+B_{w}w(t),$$

$$z(t)=C_{1}x(t)+D_{2}u(t),C_1=\left[\begin{array}{c}{Q}^{\frac{1}{2}}\\ 0\end{array}\right],D_2=\left[\begin{array}{c}0\\ R^{\frac{1}{2}}\end{array}\right],$$

$$\tau=\tau_{t}+\tau_{p}=\underbrace{\kappa\left(\frac{\texttt{Card}(K)}{c}\right)+\tau_{p}}_{\mathcal{Z}\big{(}\texttt{Card}(K),c,\tau_{p}\big{)}},$$

$$\displaystyle\boldsymbol{\mathcal{O}}_{1}:\$$

$$\text{subject to}K\in {\mathcal{V}}_{\tau },$$

$$\displaystyle\hskip 51.21504pt\ \tau=\mathcal{Z}\big{(}\texttt{Card}(K),c,\tau_{p}\big{)},$$

$$(K^{\prime},\tau^{\prime}):K^{\prime}\in\mathcal{V}_{\tau^{\prime}},\ \tau^{\prime}=\mathcal{Z}\big{(}\texttt{Card}(K^{\prime}),c,\tau_{p}\big{)},$$

$${\Phi }_{\text{L}}=\left(\begin{array}{cccc}{\Delta }_0& {\Delta }_1& -D_1^a& -D_1^b\\ \star & S_1& -D_2^a& -D_2^b\\ \star & \star & {\Delta }_2& 0\\ \star & \star & \star & {\Delta }_3\end{array}\right)\succ 0,$$

$${\Psi }_{\text{L}}=\left(\begin{array}{ccc}P& Q_0& Q_1\\ \star & R_{00}+\frac{1}{\tau }{S}_0& R_{01}\\ \star & \star & R_{11}+\frac{1}{\tau }{S}_1\end{array}\right)\succ 0,$$

$${\Delta }_0=-PA-A^{T}P-Q_0-Q_0^T-S_0,$$

$${\Delta }_1=Q_1-PBK,{\Delta }_3=3(S_0-S_1),$$

$${\Delta }_2=\tau (R_{00}-R_{11})+(S_0-S_1),$$

$$D_1^a=\frac{\tau }{2}{A}^T(Q_0+Q_1)+\frac{\tau }{2}(R_{00}+R_{01})-(Q_0-Q_1),$$

$$D_2^a=\frac{\tau }{2}{K}^T{B}^T(Q_0+Q_1)-\frac{\tau }{2}(R_{01}^T+R_{11}),$$

$$D_1^b=-\frac{\tau }{2}{A}^T(Q_0-Q_1)-\frac{\tau }{2}(R_{00}-R_{01}),$$

$$D_2^b=-\frac{\tau }{2}{K}^T{B}^T(Q_0-Q_1)+\frac{\tau }{2}(R_{01}^T-R_{11}).$$

$${\Phi }_0⪰$$

$$\alpha \ge$$

$$\gamma \ge$$

$$\eta \ge$$

$${\Psi }_0⪰$$

$$\beta ⪰$$

$$\displaystyle\boldsymbol{\mathcal{O}}_{2}:\$$

$$+\Vert \tau -{\tau }_k{\Vert }^2$$

$$\text{subject to}\Vert \tau -{\tau }_k\Vert \le \Delta ,$$

$$(K,\tau )\in {\mathcal{X}}_{(K_k,{\tau }_k)},$$

$${\left[\begin{array}{ccc}{\theta }_1=-\tau ,& \cdots ,& {\theta }_{\mathcal{N}}=0\end{array}\right]}^T,{\mathcal{M}}_i=\cos \left(\pi -\frac{i\pi }{\mathcal{N}-1}\right),$$

$${\theta }_{i+1}=\frac{\tau }{2}({\mathcal{M}}_i-1),i\in \{0,1,\dots ,\mathcal{N}-1\}.$$

$$\zeta ={\left[\underset{tail}{\underbrace{{\varphi }^T(-\tau ),{\varphi }^T({\theta }_2),\cdots ,{\varphi }^T({\theta }_{N-1})}},\underset{head}{\underbrace{x^T(t)}}\right]}^T.$$

$$\mathcal{A}\zeta =\mathcal{A}\left[\begin{array}{c}\varphi ({\theta }_1=-\tau )\\ ⋮\\ x(t)\end{array}\right]=\left[\begin{array}{c}{\partial }_{{\theta }_1}\varphi ({\theta }_1=-\tau )\\ ⋮\\ Ax(t)+K\varphi ({\theta }_1=-\tau )\end{array}\right].$$

$$M_i={\left[I_n{e}_1,\cdots ,I_n{e}_{j=i},\cdots ,I_n{e}_{\mathcal{N}}\right]}^T,$$

$$\dot{\zeta }(t)=\mathcal{A}\zeta (t)+\mathcal{B}w(t),$$

$$\mathcal{A}=\tilde{A}(\tau )-M_{\mathcal{N}}BK{M}_1^T,$$

$${\tilde{A}}_{ij}=\{\begin{array}{l}{q}_j({\theta }_i)I_n,j\in {\mathbb{N}}_{\mathcal{N}},i\in {\mathbb{N}}_{\mathcal{N}-1}\\ A,i=j=\mathcal{N}\\ 0,\text{else},\end{array}$$

$$q_i({\theta }_i)=\frac{1}{{\theta }_{ji}}\prod _{m=1,m\ne i,j}^{\mathcal{N}}\frac{{\theta }_{im}}{{\theta }_{jm}},q_j({\theta }_i)=\sum _{m=1,m\ne j}^{\mathcal{N}}\frac{1}{{\theta }_{jm}},$$

$$J_{\tau }(K,\tau )=\text{Trace}(\mathcal{CL}{\mathcal{C}}^T)=\text{Trace}({\mathcal{B}}^T\mathcal{P}\mathcal{B})$$

$$\mathcal{A}(K,\tau )\mathcal{L}+\mathcal{L}{\mathcal{A}}^T(K,\tau )=-\mathcal{B}{\mathcal{B}}^T,$$

$${\mathcal{A}}^T(K,\tau )\mathcal{P}+\mathcal{P}\mathcal{A}(K,\tau )=-{\mathcal{C}}^T\mathcal{C},$$

$${\mathcal{C}}^T\mathcal{C}=\tilde{Q}+M_1{K}^{T}RK{M}_1^T,$$

$$\displaystyle\boldsymbol{\mathcal{O}}_{3}:\$$

$$\text{s.t.}K-G=0,$$

$$\tau =\mathcal{Z}\left(\text{Card}(K),c_f,{\tau }_p\right),$$

Full text

1 Introduction

Recent achievements in engineering, computer science and communication technology have led to the development of large-scale cyber-physical systems (CPSs) that comprise of multiple dynamical systems interacting over robust communication networks. In the past two decades, significant amount of research has been done on the control-theoretic aspects of CPS, for example, see Antsaklis et al., (2013), Hespanha et al., (2007) and Sinopoli et al., (2003). Conventional CPS control designs usually rest on two idealistic assumptions, namely, that there is no limit on communication cost, and that the communication network is dedicated for control. In reality, however, network operators prefer sparse communication topologies for reducing the communication cost, and multiple users share the same communication network leading to feedback delay. Both factors give rise to important concerns on closed-loop stability and the fairness of bandwidth allocation among users.

Motivated by this problem, in this paper, we propose a network control design that lies at the intersection of three design challenges, namely - delays, sparsity of communication network, and optimal $\mathcal{H}_{2}$ performance for a LTI system. Results addressing pairs of these three are known in the literature. For example, at the junction of sparsity and optimal control, sparse linear quadratic regulators (LQR) and $\mathcal{H}_{2}$ controllers have been proposed by Lin et al., (2013) using alternating directions method of multipliers (ADMM), by Wytock and Kolter, (2013) using proximal Newton method, and by Arastoo et al., (2014) using rank-constrained convex optimization, among others. The approach has been extended to sparse output feedback in Lin and Bopardikar, (2017) using a variant of Kalman filtering.

A vast literature exists on the stability and control of LTI systems with time-delays, starting from the seminal works of Kolmanovskii and Richard, (1999), Gu et al., (2003), Krstic, (2009), and others. These papers use various analytical tools such as Riccati matrix equations, Lyapunov-Krasovskii functionals, and small-gain theorem to derive sufficient conditions for stability of time-delayed systems. Optimal control of such systems has been studied in Kushner and Barnea, (1970), Ross, (1971), Esfahani et al., (1998), and more recently in Vanbiervliet et al., (2011) and Hale et al., (2015) using spectral discretization of delays.

However, when all the three design challenges are combined together, the resulting problem becomes far more difficult to solve due to its non-convex nature as well as due to the complicated interplay between the design factors, constrained by nonlinear matrix inequalities (NMIs) arising from delay stability requirements. The problem becomes more complicated when the delay itself is a function of sparsity as eliminating communication links can free up bandwidth and reduce the per-link delay. More importantly, the interpretation of these challenges when multiple users share the same communication network is also an open question as in this case the network operator must divide the bandwidth in a fair way so that each user has a guaranteed relative $\mathcal{H}_{2}$ performance while maintaining closed-loop stability. Classical results on bandwidth allocation such as by Kelly et al., (1998) only provide protocols defined for static utility functions. Our problem, in contrast, requires bandwidth sharing for dynamic and sparse $\mathcal{H}_{2}$ control.

Summarizing the above points, the main contributions of this paper are as follows:

We present two algorithms by which one can design sparse optimal controllers for LTI systems in the presence of feedback delays. Starting from an infinite bandwidth assumption, Algorithm 1 designs a practical finite value of the network bandwidth in Section 3 by co-designing the delays and the corresponding stabilizing controller. We derive a set of convex relaxations of NMI based sufficient conditions on delay stability to perform this co-design. Thereafter, we present Algorithm 2 in Section 4, which uses ADMM to sparsify this controller while minimizing the $\mathcal{H}_{2}$ norm.

• -

We next extend the sparse controller of Section 4 to the case when multiple sets of such controllers need to be designed for $N$ decoupled systems belonging to $N$ different users, sharing the same communication network, while ensuring that each user gets a fair share of communication links. The problem is formulated in terms of minimizing the variance of the $\mathcal{H}_{2}$ performance of each individual user. Algorithm 3 is presented in Section 6 to solve this as a difference-of-convex program with MILP constraints.

• -

We validate our results using numerical simulations in Sections 5 and 7. The former illustrates the inter-dependencies between delay, sparsity and optimal control. The latter shows the efficiency of our bandwidth sharing algorithm in reaching its global minimum.

Some preliminary results on this topic have been reported in our recent conference paper Negi and Chakrabortty, (2018). The results in this paper, however, are significantly extensive in comparison. For example, Negi and Chakrabortty, (2018) only minimized a lower bound for the $\mathcal{H}_{2}$ norm of a time-delayed system using convex relaxations of bilinear matrix inequalities. In contrast, in this paper, we take a completely different approach by minimizing the exact $\mathcal{H}_{2}$ norm of the delayed system using spectral discretization. The preconditioning algorithm (Algorithm 1) in Section 3 and the bandwidth sharing algorithm (Algorithm 3) in Section 6 are also added as new contributions.

Notations: $\mathbb{R}$, $\mathbb{R}_{+}$ and $\mathbb{R}_{++}$ represent the set of real, non-negative real and positive real numbers respectively. $\mathbb{Z}_{+}$ represents the set of non-negative integers. The set $\mathbb{N}_{n}$ represents the set of natural numbers till $n$. The set $B\backslash A$ denotes the set of elements in $B$ but not in $A$. The Euclidean norm and Frobenius norm are represented by $\|\cdot\|$ and $\|\cdot\|_{\text{F}}$, respectively. Unless otherwise stated, $I_{n}$ denotes an identity matrix of order $n$. The symbol $\star$ in the lower triangular part of a matrix represents a symmetric block.

2 Problem Formulation

2.1 CPS Modeling

We consider an LTI system represented by the following state-space equation:

[TABLE]

where $x\in\mathbb{R}^{n}$, $u\in\mathbb{R}^{m}$ and $w\in\mathbb{R}^{p}$ are the state, control input and exogenous input, respectively. $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$ and $B_{w}\in\mathbb{R}^{n\times p}$ are the corresponding matrices. Assuming $(A,B_{w})$ to be stabilizable and $(A,B)$ to be controllable, our objective is to design a static state feedback controller $u(t)=-Kx(t)$, with $K\in\mathbb{R}^{m\times n}$, that minimizes the $\mathcal{H}_{2}$ norm of (1). This controller will be implemented in a distributed way using a communication network as in Fig. 1, which shows the CPS representation of the closed-loop system.

The communication of state $x_{j}$ to compute the input $u_{i}$ will suffer a delay $\tau_{ij}$ (in seconds), $i\in\mathbb{N}_{m}$, $j\in\mathbb{N}_{n}$. This delay consists of two parts, namely, $\tau_{ij}=\tau_{p_{ij}}+\tau_{t_{ij}}$, where $\tau_{p_{ij}}$ is the propagation delay and $\tau_{t_{ij}}$ is the transmission delay. $\tau_{p_{ij}}$ is defined as the ratio of link length to the speed of light, which is assumed to be a constant denoted as $\tau_{p}$ for every pair $i,j$. We assume that the communication link connecting any $j$-th sensor to any $i$-th actuator is allotted equal bandwidth, which implies that $\tau_{t_{ij}}$ is also constant for every $i,j$ (Bertsekas et al.,, 1992). We, therefore, denote $\tau_{ij}$ simply as a constant $\tau$. In reality, even if $\tau_{ij}$ deviates from $\tau$ due to traffic fluctuations and network uncertainties, our design will still hold as long as this deviation lies within the stability radius of the plant (See Proposition 1.14 and Theorem 1.16 of Michiels and Niculescu, (2014)). The $i$-th control input is written as $u_{i}(t)=-\sum_{j=1}^{n}K_{ij}x_{j}(t-\tau)$, and accordingly the closed-loop system is written as:

[TABLE]

where $Q\succeq 0$, $R\succ 0$. $K$ follows from the minimization of the $\mathcal{H}_{2}$ norm of the transfer function of (2) from input $w(t)$ to output $z(t)$. In general, $K$ is a dense matrix, which incurs high communication and computation costs. To reduce these costs, we impose a sparsity constraint on $K$. As $K$ is sparsified, the bandwidth $c$ is equally redistributed among the remaining links. From Bertsekas et al., (1992, Ch. 3), $\tau$ can be written as

[TABLE]

where $\texttt{Card}(K)$ denotes the cardinality or the number of non-zero elements in $K$ and $\kappa>0$ is a proportionality constant.The sparsity level $s$ of $K\in\mathbb{R}^{m\times n}$ is defined as $s=mn-\texttt{Card}(K)$. Equation (3) implies that $\tau$ will change as the sparsity level of $K$ changes. This change is captured by the function $\mathcal{Z}(\cdot)$.

2.2 Problem Statement

Let $\mathcal{V}_{\tau}$ be the set of $K$ that stabilize (2) for a delay $\tau$. Our objective is to design a sparse controller $K$ for the closed-loop stability of (2) while minimizing its $\mathcal{H}_{2}$ norm $J_{\tau}(K)$. The controller design problem is stated as:

[TABLE]

where $\lambda>0$ is a regularization parameter to emphasize the sparsity of $K$. We assume that the set $\mathcal{V}_{\tau}$ is non-empty. The closed form expression of $J_{\tau}(K)$ will be derived shortly in Section 4.1, where we present our main algorithm for solving $\boldsymbol{\mathcal{O}}_{1}$. Before stating that, we first develop a pre-conditioning algorithm by which we can construct an appropriate initial condition for it.

3 Preconditioning Algorithm : Co-design of $\tau$ and $K$

Given $A$, $B$, $B_{w}$ in (2), our sparsity-promoting algorithm for solving $\boldsymbol{\mathcal{O}}_{1}$ in Section 4 (namely, Algorithm 2) will require a pair

[TABLE]

as an initial condition, i.e., $(K^{\prime},\tau^{\prime})$ must stabilize (2), where \tau^{\prime}=\mathcal{Z}\big{(}\texttt{Card}(K),c,\tau_{p}\big{)}; $c$ is a nominal allotted bandwidth and $\tau_{p}$ is a known constant. The choice of $c$ plays an important role in finding this pair. For example, let $K_{o}$ be any matrix such that $(A-BK_{o})$ is Hurwitz. Then from Gu et al., (2003, Chapter 2) it follows that there exists a $\tau^{m}_{o}\geq 0$ such that $x=0$ is an asymptotically stable equilibrium of (2a) with $K=K_{o}$, $w=0$ and $\tau\in[0,\tau^{m}_{o}]$; $\tau^{m}_{o}$ is referred to as the delay margin of $K_{o}$. Thus, if $c$ is such that $\tau^{*}_{o}:=\mathcal{Z}(\texttt{Card}(K_{o}),c,\tau_{p})$ is less than or equal to $\tau^{m}_{o}$, then one can simply use $(K^{\prime},\tau^{\prime})\equiv(K_{o},\tau^{*}_{o})$ as an initial condition for Algorithm 2. However, if $c$ is too small so that $\tau^{m}_{o}<\tau^{*}_{o}$, then a different initial guess needs to be found. In this section we develop an algorithm for finding $(K^{\prime},\tau^{\prime})$ for such conservative values of $c$. In case the pair $(K^{\prime},\tau^{\prime})$ does not exist for the given $c$, the steps of this algorithm will adapt $c$ to $c^{\prime}$ so that such a pair can be found. We refer to the interval $[0,\tau^{m}_{o}]$ as a stable delay interval of $(A,B,K_{o})$. Depending on these three matrices, a system may have multiple non-overlapping stable delay intervals. This fact will also be accounted for in the algorithm.

3.1 Delay Stability Conditions

As $(K^{\prime},\tau^{\prime})$ in (5) must stabilize (2), we first recall an LMI based sufficient condition that guarantees the stability of (2) for a given $K$ and $\tau$.

Theorem 1**.**

(Gu et al.,, 2003)* The system in (2a) is asymptotically stable for a given $K$, delay ${\tau}>0$ and $w=0$ if there exist $n\times n$ matrices $P=P^{T},\ Q_{0},\ Q_{1},\ S_{0}=S^{T}_{0},\ S_{1}=S^{T}_{1},\ R_{00}=R^{T}_{00},\ R_{01},R_{11}=R_{11}^{T}$ such that the following LMIs are feasible:*

[TABLE]

For our problem, however, the matrices $P$, $Q_{0}$, $Q_{1}$, $R_{00}$, $R_{01}$, $R_{11}$, $S_{0}$, $S_{1}$ are unknown in addition to the design variables $K$ and $\tau$. As a result, the affine conditions in (6) change to nonlinear matrix inequalities (NMIs). In the following subsections, we will replace the LMI matrices $\Phi_{\text{L}}$ and $\Psi_{\text{L}}$ in (6) by $\Phi_{\text{N}}$ and $\Psi_{\text{N}}$ to signify their conversion to NMIs.

3.2 Relaxing the NMIs : Design of $K^{\prime}$ and $\tau^{\prime}$

If the NMIs in (6) could be solved then one would obtain a $(K,\tau)$ with $K\in\mathcal{V}_{\tau}$, thereby at least partly satisfying (5)111That means that the delay-stability condition is satisfied, but $\tau=\mathcal{Z}(\texttt{Card}(K),c,\tau_{p})$ is not ensured from the solution $(K,\tau)$ of the NMIs (6a)-(6b).. However, solving these NMIs is not tractable. Therefore, we first relax these NMIs to a set of convex constraints. We carry out these relaxations in the form of an iterative algorithm with initial condition $(K_{o},\tau^{m}_{o})$ where $K_{o}$ stabilizes the delay-free system (i.e., $A-BK_{o}$ is Hurwitz) and $\tau^{m}_{o}$ is the delay margin of $(A,B,K_{o})$. $\tau^{m}_{o}$ can be found from either a frequency-sweeping test as shown in Gu et al., (2003, Theorem 2.6; Example 2.9), or by checking the stability of a higher-order Pade approximation of the delayed model. The following theorem relaxes the NMIs (6a)-(6b) into a set of convex constraints to obtain a pair $(K,\tau)$ with $K\in\mathcal{V}_{\tau}$ from $(K_{o},\tau_{o}^{m})$ such that $\tau>\tau_{o}^{m}$.

Theorem 2**.**

Let the pair $(K_{o},\tau^{m}_{o})$ satisfy the LMIs of Theorem 1 with the auxiliary variables $P_{o}$, $Q_{0o}$, $Q_{1o}$, $R_{00o}$, $R_{01o}$, $R_{11o}$, $S_{0o}$ and $S_{1o}$. Let us consider perturbations $(\Delta K,\Delta\tau)$ on $(K_{o},\tau^{m}_{o})$ and $\Delta P$, $\Delta Q_{0}$, $\Delta Q_{1}$, $\Delta R_{00}$, $\Delta R_{{01}}$, $\Delta R_{{11}}$, $\Delta S_{0}$ and $\Delta S_{1}$ on the auxiliary variables as the solution of the following constraints:

[TABLE]

where $\Delta\tau\in\mathbb{R}_{+}$, $\mathcal{B}_{i}(K_{o},{\tau}^{m}_{o}),$ $i\in\mathbb{N}_{4}$, and $\mathcal{C}_{1}(K_{o},{\tau}^{m}_{o})$ are given in Appendix (Section 9.1), $\alpha,\beta\in\mathbb{R}_{++}$ are optimization variables, and $\epsilon,\gamma,\eta\in\mathbb{R}_{++}$ are scalar design variables. Then, the pair $(K=K_{o}+\Delta K,\tau=\tau^{m}_{o}+\Delta\tau)$ satisfies Theorem 1 with the auxiliary variables $P=P_{o}+\Delta P$, $Q_{0}=Q_{0o}+\Delta Q_{0}$, $Q_{1}=Q_{1o}+\Delta Q_{1}$, $R_{00}=R_{00o}+\Delta R_{00}$, $R_{01}=R_{01o}+\Delta R_{01}$, $R_{11}=R_{11o}+\Delta R_{{11}}$, $S_{0}=S_{0o}+\Delta S_{0}$ and $S_{1}=S_{1o}+\Delta S_{1}$. We use the shorthand notation of $(K,\tau)\in\mathcal{X}_{(K_{o},\tau^{m}_{o})}$ to denote that $(K,\tau)$ is obtained from a given pair $(K_{o},\tau^{m}_{o})$ using (7).

Proof.

The proof is given in Appendix (Section 9.1). ∎

We next use Theorem 2 to develop an algorithm starting from $(K_{o},\tau^{m}_{o})$, where $\Delta K$ and $\Delta\tau$ are evaluated iteratively such that we eventually obtain a stabilizing pair $(K^{\prime},\tau^{\prime})$ that fully satisfies (5), i.e., $\tau^{\prime}=\mathcal{Z}(\texttt{Card}(K^{\prime}),c,\tau_{p})$. Equations (7a)-(7f) in Theorem 2 provide the set of constraints that are required to search for a new stabilizing pair $(K_{k+1}=K_{k}+\Delta K_{k},$ $\tau_{k+1}=\tau_{k}+\Delta\tau_{k})$ in the $(k+1)$-th iteration of this algorithm. We incorporate these constraints using the conjugate gradient Steighaug’s method (Steihaug,, 1983) as shown next.

Let us define $f(\tau)=(T-\tau)^{2}$, where $T$ is a constant. Starting from the known stabilizing pair $(K_{k},\tau_{k})$ for (2), $T$ can be chosen as $T\geq\mathcal{Z}(\texttt{Card}(K_{k}),c,\tau_{p})$ and the $(k+1)$-th iteration of the algorithm can be obtained as the solution of the following program:

[TABLE]

where $\Delta$ is the radius of the trust-region, and $\mathcal{X}_{(K_{k},{\tau}_{k})}$ is a subset of the stabilizing pairs from Theorem 2. Algorithm 1 describes the iterative solution for (8). Details on the practical choice of $\epsilon$, $\eta$ and $\gamma$ can be found in Negi and Chakrabortty, (2018). The following Theorem provides the conditions under which Algorithm 1 converges.

Theorem 3**.**

Using Algorithm 1, a sequence $\{\Delta{\tau}_{k}\}\subset\mathbb{R}_{+}$ is established in $k$ iterations with $\lim_{k\to\infty}\Delta{\tau}_{k}=0$ if $\eta$ and $\gamma$ are chosen such that the convex program in Theorem 2 is feasible for iterations $\{1,2,\ldots,k\}$.

Proof.

Constraint set in (8) satisfies linear constraint qualification. Additionally, $(K^{\prime}_{i+1},{\tau}^{\prime}_{i+1})\in\mathcal{X}_{(K^{\prime}_{i},{\tau}^{\prime}_{i})}$ and $\Delta\tau_{i}\in\mathbb{R}_{+}$ is ensured for every $i$-th iteration through (7) and the assumption on the values of $\eta$ and $\gamma$. Therefore, using Steihaug, (1983, Theorem 3.3), global convergence of (1) is guaranteed. ∎

If at any $k$-th iteration, we obtain a pair $(K_{k},\tau_{k})$ such that $K_{k}$ stabilizes (2) for $\tau^{*}_{k}=\mathcal{Z}(\texttt{Card}(K_{k}),c,$ $\tau_{p})$, then the test in (5) is passed. In that case, we stop Algorithm 1, set $(K^{\prime},\tau^{\prime})\equiv(K_{k},\tau^{*}_{k})$ and allocate a final bandwidth as $c_{f}=c$. However, depending on where we start from, Algorithm 1 may converge to a pair $\lim_{k\to\infty}(K_{k},\tau_{k})$ such that $K_{k}$ stabilizes (2) for a stable delay interval $[\tau_{k-},\tau_{k+}]$ ($\tau_{k}$ belongs to this interval) but does not stabilize it for $\tau^{*}_{k}=\mathcal{Z}(\texttt{Card}(K_{k}),c,\tau_{p})$. In that case, instead of re-running the algorithm with a different initial pair $(K_{o},\tau^{m}_{o})$, an easier alternative is to set $(K^{\prime},\tau^{\prime})\equiv(K_{k},\tau_{k+})$ if $\tau^{*}_{k}>\tau_{k+}$ or $(K^{\prime},\tau^{\prime})\equiv(K_{k},\tau_{k-})$ otherwise, and update $c$ to $c^{\prime}=\kappa\frac{\texttt{Card}(K^{\prime})}{\tau^{\prime}-\tau_{p}}$. The final bandwidth is then updated as $c_{f}=c^{\prime}$. We next proceed to Algorithm 2 with $(K^{\prime},\tau^{\prime})$ as the initial condition satisfying (5) with bandwidth $c_{f}$.

4 Main Algorithm

4.1 $\mathcal{H}_{2}$ norm of the Delayed Model

We use the method of spectral discretization for approximating the time-delayed system (2) to a finite dimensional LTI system, as derived in (Vanbiervliet et al.,, 2011). $\mathcal{H}_{2}$ norm of the original time-delayed system is then approximated as that of this LTI system. Considering any stabilizing tuple $(K,\tau)$, the interval $[-{\tau},0]$ is first divided into a grid of $\mathcal{N}$ Chebyshev extremal points

[TABLE]

We define function segments $\phi(\theta)=x(t+\theta),\ -{\tau}\leq\theta\leq 0$, which are used to represent both the past (tail) state $\phi(\theta_{i})$, $i\in\mathbb{N}_{\mathcal{N}-1}$ and the present (head) state $x(t)$ in the augmented state vector $\zeta\in\mathbb{R}^{\mathcal{N}n}$ as

[TABLE]

Let $\mathcal{W}:=\mathcal{L}_{2}([-\tau,0],\mathbb{R}^{n})\times\mathbb{R}^{n})$ be defined as the vector space containing $\phi(\theta)$, $-\tau\leq\theta\leq 0$. The discretization in (9) replaces the infinite dimensional space $\mathcal{W}$ by $\mathbb{W}$, which is the space of discrete functions over the $\mathcal{N}$ points in (9). A linear operator $\mathcal{A}$ is defined to act on $\zeta=(\phi,x)\in\mathbb{W}$ such that it mimics the behavior of the delayed system (2) in the following way:

[TABLE]

Let us define matrices $M_{i}$, $i\in\mathbb{N}_{\mathcal{N}}$ as follows:

[TABLE]

where $e_{j}=1$ if $j=i$ and [math] otherwise. The approximated LTI system can then be written as

[TABLE]

where $\mathcal{B}=M_{\mathcal{N}}B_{w}\in\mathbb{R}^{\mathcal{N}n\times p}$ is the exogenous input matrix for this augmented system. $\mathcal{A}\in\mathbb{R}^{\mathcal{N}n\times\mathcal{N}n}$ represents the closed-loop state matrix, which can be further decomposed as

[TABLE]

where $\tilde{A}(\tau)$ is defined using (9) as

[TABLE]

and $\theta_{pq}=\theta_{p}-\theta_{q}$. The $\mathcal{H}_{2}$ norm of the time-delayed model can finally be approximated by the $\mathcal{H}_{2}$ norm of (13) which is defined as follows:

[TABLE]

where $\tilde{Q}=\texttt{BlkDiag}\left(\boldsymbol{0},Q\right)$, and $\mathcal{P},\mathcal{L}\in\mathbb{R}^{\mathcal{N}n\times\mathcal{N}n}$ are positive-definite Lyapunov matrices. The $\mathcal{H}_{2}$ norm $J_{\tau}$ is finite and positive for a stable delayed system. We next use this expression as the objective function in the ADMM formulation for solving $\boldsymbol{\mathcal{O}}_{1}$.

4.2 Design Formulation using ADMM

As $\boldsymbol{\mathcal{O}}_{1}$ has two conflicting objectives, namely, sparsity promotion and minimization of the $\mathcal{H}_{2}$ norm, it fits the structure of ADMM-based optimization. Lin et al., (2013) formulated an ADMM based design for the delay-free LTI system. In our problem, we include the added objective of maintaining stability of (2a) while updating the delay as a function of sparsity. We relax $\boldsymbol{\mathcal{O}}_{1}$ in the ADMM form as follows:

[TABLE]

where $c_{f}$ and $\tau_{p}$ are known constants, and $\lambda$ is a known regularization parameter. Due to the non-convex nature of $\texttt{Card}(K)$, sparsity is promoted through its closest convex approximation, which is the weighted $l_{1}$ norm $g(G)$. The norm is computed iteratively as

[TABLE]

where $G_{i}$ is known from the $i$-th iteration and $0<\epsilon_{1}\ll 1$. $\mathbf{\mathcal{O}}_{3}$ is solved using two loops. The inner ADMM loop will update $K$ and $G$ assuming $\tau$ to be constant while the outer loop will update $\tau$ according to (17c), as shown shortly in Section 4.3. The augmented Lagrangian for the inner ADMM loop can be written as:

[TABLE]

where $\Lambda\in\mathbb{R}^{m\times n}$ is the dual variable and $\rho$ is a constant. $\lambda$ is increased in steps to allow gradual sparsity promotion. For the rest of this subsection, we continue with the formulation of the inner-loop. The outer-loop for $\tau$-update is presented in the next subsection.

At any $(k+1)$-th iteration of the inner loop, the algorithm alternates between the following three minimization problems:

$\bullet$ $K$-min : Minimizes the $\mathcal{H}_{2}$ norm of (2) as

[TABLE]

$\bullet$ $G$-min : Minimizes the weighted $l_{1}$ norm as:

[TABLE]

$\bullet$ $\Lambda$-Update : $\Lambda_{k+1}=\Lambda_{k}+(K_{k+1}-G_{k+1})$.

The $K$-min step can be solved using Anderson-Moore method (Rautert and Sachs,, 1997). This method requires the gradient of $J_{\tau}$, which is derived as follows.

Theorem 4**.**

The gradient of $J_{\tau}(K)$ is given as:

[TABLE]

where $\mathcal{L}(K)$ and $\mathcal{P}(K)$ are obtained from (16b) and (16c), respectively.

Proof.

Taking the partial derivative of (16c), we get

[TABLE]

Partial derivative of the $\mathcal{H}_{2}$ norm in (16a) is

[TABLE]

Multiplying (16c) by $\partial\mathcal{P}$ and taking its trace, we get

[TABLE]

Using (24) and (25), we can write

[TABLE]

We next multiply (23) by $\mathcal{L}$ and take its trace. Using (26), $\text{Tr}(dK^{T}\nabla J_{\tau}(K))$ is equivalent to

[TABLE]

The final expression for the gradient (22) follows from the above equation. ∎

We next use Theorem 4 to derive the solution of the $K$-min step, followed by the solution of the $G$-min step.

4.2.1 Solving the $K$-min Step

To solve the $K$-min step, we employ the Anderson-Moore method using the gradient of the objective function in (20). Denoting $V_{k}=G_{k}-\frac{1}{\rho}\Lambda_{k}$ and holding $\tau$ constant for this step, we use Theorem 4 to derive the necessary optimality condition for (20) as:

[TABLE]

where $K^{+}$ is the descent direction, and $\mathcal{L}(K_{k})$ and $\mathcal{P}(K_{k})$ are the solutions of (16b) and (16c), respectively. The Armijo-Goldstein rule for determining the step-size guarantees $\lim_{k\to\infty}\|\nabla J_{\tau}(K_{k})\|=0$ (Rautert and Sachs,, 1997, Theorem 4.2).

4.2.2 Solving the $G$-min Step

The solution of the $G$-min step is well-known in the literature, and is given by the soft-thresholding operator (Boyd et al.,, 2011). Let $K_{pq}$ be the $(p,q)$-th element of the solution of the $K$-min step. Denoting $U_{pq}=\Lambda_{pq}+\rho K_{pq}$, the solution of the $G$-min step is given as

[TABLE]

A high value of $\lambda$ favors the promotion of zero entries in $G$. The ADMM loop terminates in $K\in\mathcal{I}_{s}$ upon reaching the convergence criteria given in Boyd et al., (2011, Section 3.3), where $\mathcal{I}_{s}$ represents a set of matrices with a specific sparsity structure. A polishing step will be introduced in Section 4.4 for further improving $K$.

4.3 Update of $\tau$

The delay $\tau$ is assumed to be constant in the inner loop since it is a function of $\texttt{Card}(K)$ (and not $K$). Moreover, Theorem 4 can guarantee a solution of $K$-min only if $\tau$ remains constant. For these reasons, $\tau$ is updated with respect to $\texttt{Card}(K)$ (i.e., according to (17c)) at the beginning of every re-weighting step (18), which we refer to as the outer loop. Thus, at every outer loop iteration, sparsity level of $K$ decides the value of $\tau$ and, in turn $\tau$ affects the future value of sparse $K$ through the three steps of ADMM. This is another difference between our proposed algorithm and that in Lin et al., (2013).

Before proceeding with the outer loop design, we first recall a fundamental property from the stability theory of delayed systems following Gu et al., (2003, Chapter 2). The characteristic polynomial for (2a) with $w=0$ is

[TABLE]

If $\tau^{m}$ is the delay margin of $(A,B,K)$ in (2a), then the frequency $\omega$ at which $p(j\omega;e^{-j\tau^{m}\omega})=0$ becomes the first crossing of the characteristic roots of (29) from the stable left hand complex plane to the unstable right half complex plane, also known as a zero crossing. For any given system $(A,B,K)$, multiple but finite zero crossings exist, which indicate that (2a) is also stable for some delay intervals other than $[0,\tau^{m}]$. We assume that $\mathfrak{n}$ such crossings exist, and denote the crossings as $\omega_{k}$, $k\in\mathbb{N}_{\mathfrak{n}}$. Also, we define a dimensionless quantity $\hat{\omega}_{k}\in[0,2\pi]$ such that $p(j\omega;e^{-j\hat{\omega}_{k}})=0$ $\forall\ \omega>0$. Defining $\nu_{k}=\nicefrac{{\hat{\omega}_{k}}}{{\omega_{k}}}$, we assume without loss of generality that

[TABLE]

where $\nu_{k}$, $k\in\mathbb{N}_{\mathfrak{n}}$. The system in (2a) is stable for any $\tau\in(\nu_{k},\nu_{k+1})$ if it is known to be stable for some $\tau^{*}$ in that interval. As mentioned in Section 3, such continuous intervals are called the stable delay intervals of $(A,B,K)$. We next derive the update law for $\tau$ using the concept of stable delay intervals.

Let the $i$-th re-weighting step start with the stabilizing pair $(K_{i},\tau_{i})$ and end with the stabilizing pair $(K_{i+1},\tau_{i})$, where $\texttt{Card}(K_{i+1})\leq\texttt{Card}(K_{i})$. In the next iteration, $\tau_{i}$ should be updated to $\mathcal{Z}(\texttt{Card}(K_{i+1}),c_{f},\tau_{p})$, which is denoted by $\tau^{*}_{i+1}$. The following scenarios can arise depending on whether the pair $(K_{i+1},\tau^{*}_{i+1})$ is stabilizing or not:

4-1

If $(K_{i+1},\tau^{*}_{i+1})$ stabilizes (2a), we set $\tau_{i+1}=\tau^{*}_{i+1}$ and start the ADMM loop. We illustrate this case in 4-1 panel of Fig. 2.

4-2

If $(K_{i+1},\tau^{*}_{i+1})$ does not stabilize (2a), then $\tau_{i+1}$ cannot be updated to $\tau^{*}_{i+1}$ since it does not belong to any stable delay interval of $(A,B,K_{i+1})$. To ensure that the bandwidth does not increase, we look for the closest stable delay intervals towards the right of $\tau^{*}_{i+1}$ as shown in panels 4-2a and 4-2b of Fig. 2. Let us denote the stable interval in which $\tau_{i}$ belongs to as $[\nu_{k},\nu_{k+1}]$, where $k\in\mathbb{N}_{\mathfrak{n}}$. Two situations can arise in this case. The first situation (panel 4-2a) is when the closest stable delay interval on the right of $\tau^{*}_{i+1}$ is $[\nu_{k},\nu_{k+1}]$. We denote the start of this interval by $\nu^{*}$, and update $\tau_{i+1}=\nu^{*}$. In the second situation (panel 4-2b), the closest stable delay interval to the right of $\tau^{*}_{i+1}$ is $[\nu_{r},\nu_{r+1}]$, where $r<k$. In this case, assuming that one can find this stable interval, we denote $\nu^{*}=\nu_{r}$ and set $\tau_{i+1}=\nu^{*}$. In reality, however, for this case it may be difficult to find the exact location of the nearest interval $[\nu_{r},\nu_{r+1}]$ numerically. Since the location of $[\nu_{k},\nu_{k+1}]$ is already known from the $i$-th iteration, $\tau_{i+1}$ can be updated to $\nu^{*}=\nu_{k}$ for this case as well.

Let $(K_{f},\tau_{f})$ denote the output of Algorithm 2. If 4-1 is encountered then $\tau_{f}=\mathcal{Z}(\texttt{Card}(K_{f}),c_{f},\tau_{p})$. If 4-2 is encountered, we may obtain $(K_{f},\tau_{f})$ such that $\tau_{f}>\mathcal{Z}(\texttt{Card}(K_{f}),c_{f},\tau_{p})$, which can be accommodated by decreasing the allotted bandwidth $c_{f}$. While this saves us bandwidth, it may cost us $\mathcal{H}_{2}$ performance as we obtain a controller for a higher feedback delay than necessary. Also, in this case it is easy to verify that the difference between $\tau_{f}$ and $\tau_{f}^{*}$ (where, $\tau_{f}^{*}=\mathcal{Z}(\texttt{Card}(K_{f}),c_{f},\tau_{p})$) is always bounded as

[TABLE]

where $\epsilon=\max_{i}(\tau^{*}_{i+1}-\tau_{i+1})$, and $r_{max}$ is the total number of re-weighting steps.

4.4 Polishing Step

The final step of the inner ADMM loop is to refine its output $K\in\mathcal{I}_{s}$ to $K^{\prime}\in\mathcal{I}_{s}$ such that the latter produces the least $J_{\tau}$ attainable for the sparsity set $\mathcal{I}_{s}$. The gradient of $J_{\tau}(K)$ for fixed $\mathcal{I}_{s}$ is given by

[TABLE]

The ADMM loop of the $(i+1)$-th re-weighting step produces a stabilizing tuple $(K_{i+1}\in\mathcal{I}_{s},\tau_{i})$, as shown in line-17 of Algorithm 2. Newton’s method can be used to construct the polishing step with a quadratic convergence rate. Let the output of the polishing step be

[TABLE]

where $\tilde{K}\in\mathcal{I}_{s}$ is the Newton’s descent direction and $\mathbf{s}$ is the appropriately chosen step-size. Using Rautert and Sachs, (1997), (34) is written in the form of the following necessary optimality condition:

[TABLE]

where $\mathcal{L}(K_{i+1})$ and $\mathcal{P}(K_{i+1})$ are computed from (16b) and (16c), respectively. The matrices $Z_{1}=\partial\mathcal{L}$ and $Z_{2}=\partial\mathcal{P}$ are computed by differentiating (16b) and (16c) to obtain the following equations:

[TABLE]

where $\displaystyle G_{1}=M_{\mathcal{N}}B\tilde{K}M^{T}_{1}\mathcal{L}$, $G_{2}=\mathcal{P}M_{\mathcal{N}}B\tilde{K}M^{T}_{1}-M_{1}K$ $R\tilde{K}M_{1}^{T}$ and $\mathcal{A}$ follows from (14) with $K_{i+1}$ as the feedback gain. We first compute the gradient $\nabla J_{\tau}(K_{i+1})|_{\mathcal{I}_{s}}$ and then solve (35), (36) and (37) for the unknown variables $Z_{1}$, $Z_{2}$ and $\tilde{K}$. Beginning from $K_{i+1}$, a sequence of decreasing $\{J_{\tau}(K^{\prime}_{i+1})\}$ is generated by using Armijo-Goldstein rule to choose step-size $\mathbf{s}$ in (33). Conjugate gradient method can be used to compute the Newton direction similar to Fardad et al., (2009, Section IV-B).

5 Examples

We illustrate the effectiveness of algorithms 1 and 2 using numerical examples. We refer to the two algorithms jointly as a delay-aware algorithm, as compared to the delay-agnostic algorithm of Lin et al., (2013) where no feedback delay was considered. We present three case studies, each for the following two conditions:

S1 $\tau$ is constant.

S2 $\tau$ changes with $\texttt{Card}(K)$ following (3).

5.1 Case-1 : Element-wise Sparsity

We consider a $50$-th order LTI system with a randomly generated state matrix $A$. The numerical value of $A$ can be found in the data set (Git Hub) in Negi, (2019). We assume $B=B_{w}=I_{50}$, $\tau_{p}=9.83$ ms, $c=956$ and $\kappa=0.01$. Algorithm 1 starts with a full matrix $K_{o}$ and corresponding $\tau^{m}=19.7$ ms, and outputs a full $K^{\prime}$ with the corresponding $\tau^{\prime}=36$ ms such that $c_{f}=c$. Algorithm 2 then sparsifies $K^{\prime}$ while maintaining optimality. Fig. 4(a) and 4(a) compare the dynamic performance of the closed-loop states between the delay-aware and delay-agnostic algorithms for conditions S1 (i.e., when $\tau$ is kept constant at $36$ ms), and S2 (i.e, when $\tau$ decreases from $36$ ms for $s=0$ to $12.42$ ms for $s=2252$222$s$ is the sparsity level of $K$ defined as $s=mn-\texttt{Card}(K)$.), respectively. Fig. 4(b) and 4(b) show the comparison of $J_{\tau}$ versus $\texttt{Card}(K)$ between the delay-agnostic and delay-aware algorithms for S1 and S2, respectively. For S1, we can see in Fig. 4(b) that while $K$ obtained by both algorithms are of comparable sparsity levels, the one resulting from the delay-agnostic algorithm easily destabilizes the closed-loop system. In contrast, the delay-aware algorithm consistently maintains closed-loop stability and $\mathcal{H}_{2}$ optimality of the delayed system. Similar conclusions can be drawn for S2 as seen in Fig. 4(b). Since $J_{\tau}$ is not defined for unstable systems, a nominal value of $J_{\tau}=-1000$ is allotted for representation.

The upper panel of Fig. 4(b) for S1 shows an initial decrease in $J_{\tau}$ till it reaches a minimum, followed by a gradual increase as sparsity further increases. A similar trend is observed for S2 in Fig. 4(b) but here the minimum value of $J_{\tau}$ is reached at a much higher sparsity level compared to S1. This is because S1 forces $\tau$ to remain at the initial high value of $36$ ms, while in S2 $\tau$ decreases with increasing sparsity, which aids in achieving a higher sparsity level for a lower $J_{\tau}$. The initial detrimental effect of increasing sparsity on $J_{\tau}$ is weaker compared to the aiding effect of the decreasing delay. After the minimum $J_{\tau}$ is reached, the effect of increasing sparsity takes over and $J_{\tau}$ starts increasing monotonically. Algorithm 2 successfully captures this minimum.

5.2 Case-2 : Block Sparsity

Next, we consider a $10$-th order LTI system. The numerical values of the randomly generated state matrix $A$ and the control input matrix $B=\texttt{BlkDiag}\left(B_{1},B_{2},B_{3}\right)$, with $B_{1}\in\mathbb{R}^{2\times 2}$, $B_{2}\in\mathbb{R}^{3\times 3}$ and $B_{3}\in\mathbb{R}^{5\times 5}$ are provided in Negi, (2019). We assume $B_{w}=I_{10}$. This case presents the extension of our design to a block sparse $K$. Again, as shown in Fig. 6(b) and Fig. 6(b) respectively, our algorithm guarantees closed-loop stability while the delay-agnostic algorithm fails to do so for both conditions S1 and S2. Also, delay-aware algorithm for S2 achieves higher sparsity for a lower $J_{\tau}$ compared to S1 due to the flexibility in changing $\tau$.

5.3 Case-3 : $K$ is not stabilizing for the delay-free system

One benefit of Algorithm 1 is that we can obtain a stabilizing pair $(K^{\prime},\tau^{\prime})$ where $\tau^{\prime}=\mathcal{Z}(\texttt{Card}(K^{\prime}),c,$ $\tau_{p})$ lies beyond the delay margin interval $[0,\tau^{\prime m}]$ of $(A,B,K^{\prime})$, as explained in Section 4.3. This means that even for a relatively low value of the bandwidth $c$ and a high value of $\tau_{p}$, Algorithm 1 can successfully find a stable pair $(K^{\prime},\tau^{\prime})$, thereby fulfilling the condition in (5), and preventing the case where one may need to revise $c$ to a higher bandwidth $c^{\prime}$. The trade-off is that such a starting pair for Algorithm 2 may result in the 4-2 scenario of Section 4.3, where we lose on $\mathcal{H}_{2}$ optimality. This case describes such a scenario. We consider a $10$-th order LTI system with randomly generated state matrix $A$ (Negi,, 2019) and $B=B_{w}=I_{10}$. For $c=10.5$, $\tau_{p}=28.34$ ms, and $\kappa=0.01$ (compared to Case-1, $c$ is lower and $\tau_{p}$ is higher). Algorithm 1 starts with a full matrix $K_{o}$ and the corresponding $\tau^{m}_{o}=56.5$ ms. It converges to the desired $\tau^{\prime}=123.5$ ms corresponding to a full matrix $K^{\prime}$. However, the pair $(K^{\prime},\tau^{\prime})$ is such that $K^{\prime}$ is not stabilizing for $\tau=0$. Therefore, a comparison with the delay-agnostic algorithm is not possible. The sparsity behavior is seen to be similar to Case-1, as shown in Fig. 8(b) and 8(b) for S1 and S2, respectively. For S1, $\tau$ is kept constant at $123.5$ ms, and the minimum $J_{\tau}=375.5$ is achieved at $s=12$. Thereafter, $J_{\tau}$ increases steeply with sparsity. In S2, $\tau$ decreases from $123.5$ ms to $97.87$ ms as we reach $s=58$. We see a steady decline in $J_{\tau}$ such that minimum $J_{\tau}=273.4$ is reached at $s=49$, thereby confirming that for S2, a lower $\mathcal{H}_{2}$ norm is obtained at a much higher sparsity level compared to S1.

6 Cooperative Coupled Sparsification

Algorithms 1 and 2 have provided the designs for sparse $K$ when there is a single user consuming the entire available bandwidth. We next extend this design to a multi-user system, which consists of multiple dynamically decoupled LTI systems that share the same set of communication links for their individual closed-loop control. Thus, the closed-loop performance of the users are indirectly dependent on each other even though they are physically decoupled. This section presents a characterization of the equivalent $\mathcal{H}_{2}$ performance between the users followed by an algorithm that enables sparse controller designs by fairly allocating the number of available communication links.

6.1 Problem Formulation

Let there be ${N}$ users sharing $\ell$ communication links. Let the $i$-th user’s closed-loop state-space model be represented as

[TABLE]

where $\tau_{i}$ is the feedback delay associated with the $i$-th system, $x_{i}\in\mathbb{R}^{n_{i}}$ are the states, $u_{i}\in\mathbb{R}^{m_{i}}$ are the control inputs and $K_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ are the feedback gains. Ideally, one would need $\sum_{i=1}^{N}m_{i}n_{i}$ number of links to accommodate an optimal state-feedback controller for every user. If $\ell$ is less than this number then one must promote at least $\mathfrak{s}$ number of zero entries in $K_{1},\ldots,K_{N}$, where $\mathfrak{s}=\sum_{i=1}^{N}m_{i}n_{i}-\ell$. Let $s_{i}\in\mathbb{Z}_{+}$ be the sparsity level of the $i$-th user, defined as the number of zero entries promoted in $K_{i}$. This means that $\sum_{i=1}^{N}s_{i}=\mathfrak{s}$. $K_{i}$ can be designed for different sparsity levels $s_{i}$ using Algorithm 1 and 2. Since multiple structures of $K_{i}$ are possible for each such $s_{i}$, we choose that $K_{i}$ which results in the lowest $J_{\tau_{i}}$. This generates a $J_{\tau_{i}}$ vs $s_{i}$ characteristic curve for every user $i$. An example of this characteristic curve is seen in Fig. 4(b). To normalize these curves over $N$ users, let $J^{o}_{i}=\text{min}\ J_{\tau_{i}}$ be called as the nominal performance of the $i$-th user. We define the performance cost ratio $r_{i}$ of the $i$-th user as

[TABLE]

The ratio curve $r_{i}(s_{i})\geq 1$ is the normalized characteristic of the $i$-th user; see Fig. 9(a). Let $s=[s_{1},\ldots,s_{N}]^{T}\in\mathbb{R}^{N}$ be the vector of sparsity levels that will be allotted to the users, with the corresponding performance ratios being $r(s)=[r_{1}(s_{1}),\ldots,r_{N}(s_{N})]^{T}\in\mathbb{R}^{N}$. We denote $H(r_{1}(s_{1}),\ldots,$ $r_{N}(s_{N}))$ as the variance of the $\mathcal{H}_{2}$ performance ratio between the users defined as

[TABLE]

The variance $H\geq 0$ measures the fairness of link allotment. We state our problem as:

• •

Find $s=[s_{1},\ldots,s_{N}]$ to minimize the variance $H(r(s))$ of $\mathcal{H}_{2}$ performance ratio between the users such that $\mathfrak{s}=\sum_{i=1}^{N}s_{i}$, where $\mathfrak{s}$ is a known number.

To present the above objective as an optimization problem, we define the following sets:

[TABLE]

$\mathcal{S}_{i}$ is the set of possible sparsity levels of the $i$-th user. $\mathcal{S}$ is the feasible set of our problem. The optimization problem can be stated as:

[TABLE]

where $\sigma$ is a positive scalar parameter, $\mathfrak{s}\in\mathcal{S}_{M}$ and $h(r(s))=\sum_{i=1}^{N}r_{i}(s_{i})$. The term $h(\cdot)$ is added to ensure that the sum of the ratios $r_{i}$ do not become too large. For instance, between $r=[1.33,1.33]^{T}$ and $r=[1.32,1.31]^{T}$ for the same $\mathfrak{s}$, the latter would be preferable even though the former has a lower $H$. We choose $0<\sigma\ll 1$ to ensure that $\boldsymbol{\mathcal{O}}_{4}$ places a higher relative importance on minimizing the variance.

The optimization problem in (45) is defined over $\mathcal{S}\subset\mathbb{Z}^{N}$, and is therefore, combinatorial. Moreover, our objective in (45a) is defined as a functional on the desired domain $\mathcal{S}$, in contrast to a standard optimization problem. We first propose a framework for relaxing the domain of $\boldsymbol{\mathcal{O}}_{4}$ from $\mathbb{Z}^{N}$ to $\mathbb{R}^{N}$, followed by introducing a convex-concave programming (CCP) based algorithm that solves the relaxed program.

6.2 Relaxing the Optimization Constraints

Since $\mathcal{S}_{i}$ is a discrete set $\forall$ $i\in\mathbb{N}_{N}$, $r_{i}(s_{i})$ can be represented as a piece-wise affine function with discontinuities at every point in $\mathcal{S}_{i}$. Without loss of generality, we can assume that the discrete set $\mathcal{S}_{i}$ for the $i$-th user has $p_{i}$ number of elements ordered as

[TABLE]

We define the continuous intervals $\mathcal{R}_{i,j}=[s_{i,j},s_{i,j+1}]$, $j\in\mathbb{N}_{p_{i}-1}$ such that $\mathcal{R}_{i}:=\bigcup^{p_{i}-1}_{j=1}\mathcal{R}_{i,j}$ for the $i$-th user.

Example (a) - If user $1$ has $p_{1}=4$ sparsity levels $\mathcal{S}_{1}=\{0,3,4,5\}$, then $\mathcal{R}_{1,1}=[0,3]$, $\mathcal{R}_{1,2}=[3,4]$, $\mathcal{R}_{1,2}=[4,5]$, and $\mathcal{R}_{1}=[0,5]$. If user $2$ has $p_{2}=4$ with $s_{2}=\{0,2,4,5\}$ and $\mathfrak{s}=8$, then the feasible set of $\boldsymbol{\mathcal{O}}_{4}$ is $\{(4,4),(3,5)\}$.

We can represent $r_{i}(s_{i})$ for the $i$-th user as follows:

[TABLE]

Since the piecewise-affine characteristics of $r_{i}(s_{i})$ is defined over finite real set $\mathcal{R}_{i}$ $\forall$ $i$, $\boldsymbol{\mathcal{O}}_{4}$ will search over closed and bounded sets. We next reformulate the integer program of $\boldsymbol{\mathcal{O}}_{4}$ as an MILP with the non-convex objective of (45a). Table 1 defines the optimization variables for an $i$-th user with the sparsity level $s_{i}\in\mathcal{R}_{i}$ and the corresponding $r_{i}(s_{i})$ piece-wise defined over $p_{i}-1$ intervals. The $p_{i}-2$ binary variables $\delta_{i,j}\in\{0,1\}$ determine the interval $\mathcal{R}_{i,j}$ for $s_{i}$. $\tilde{z}_{i}\in\mathbb{R}$ represents the corresponding piece of $r_{i}(s_{i})$. Using Table 1, we define vector $v_{i}$ for every $i\in\mathbb{N}_{N}$, which represents the $i$-th user’s decision variable in the reformulated MILP:

[TABLE]

Table 2 translates the definitions of Table 1 into mixed integer linear inequalities, thereby, relaxing the combinatorial constraints of $\boldsymbol{\mathcal{O}}_{4}$. The inequalities in Table 2 can be written in the form of compact LMIs $G_{i}v_{i}\leq g_{i}$, $i\in\mathbb{N}_{N}$. We reformulate $\boldsymbol{\mathcal{O}}_{4}$ as

[TABLE]

Let $\mathbb{S}_{\mathfrak{s}}$ define the solution set of $\boldsymbol{\mathcal{O}}_{5}$ for $\mathfrak{s}\in\mathcal{S}_{M}$ as:

[TABLE]

Since $\mathbb{S}_{\mathfrak{s}}$ is the Cartesian product of closed and bounded sets $\mathcal{R}_{i}$, $i\in\mathbb{N}_{N}$, it is also closed, bounded and non-empty. The objective function $F(r_{1}(s_{1}),\ldots,r_{N}(s_{N}))$ is continuous over the sets $\mathcal{R}_{i}$, which are, in turn, continuous over $s_{i}$ from (47). Therefore, $F$ is continuous over $s_{i}$, $i\in\mathbb{N}_{N}$ and a global minima $s^{*}=(s^{*}_{1},\ldots,s^{*}_{N})$ of $F$ with $\sum_{i=1}^{N}s^{*}_{i}=\mathfrak{s}$ exists in $\mathbb{S}_{\mathfrak{s}}$ from Weierstrass’s theorem.

While $\boldsymbol{\mathcal{O}}_{4}$ is defined over $\mathcal{S}\subset\mathbb{Z}^{N}_{+}$, $\boldsymbol{\mathcal{O}}_{5}$ is defined over $\mathbb{S}\subset\mathbb{R}^{N}_{+}$, where

[TABLE]

This implies that $\mathcal{S}\subset\mathbb{S}$. We next modify the $r_{i}(s_{i})$ characteristics of all users to ensure that the local minima of $\boldsymbol{\mathcal{O}}_{5}$ always lie in $\mathcal{S}$ despite being defined over the much larger set $\mathbb{S}$.

Example (a) [Contd.] The local minima for $\boldsymbol{\mathcal{O}}_{5}$ should belong to $\{(4,4),(3,5)\}$. However, since $\mathcal{R}_{1}=\mathcal{R}_{2}=[0,5]$ are continuous real sets on which the piece-wise affine characteristics of (47) are defined, the local minimum can very well be at some point $(4.5,3.5)$, which is unacceptable.

6.3 Modifying $r_{i}(s_{i})$ Characteristics to Accommodate Integer Sparsity Levels

To ensure that any local minima of $\boldsymbol{\mathcal{O}}_{5}$ belongs to $\mathcal{S}$ instead of $\mathbb{S}\backslash\mathcal{S}$, the following theorem modifies the piece-wise affine $r_{i}(s_{i})$ at $s_{i}\notin\mathcal{S}_{i}$ $\forall$ $i\in\mathbb{N}_{N}$.

Theorem 5**.**

Let $r_{i}(s_{i})$ be defined on intervals $\mathcal{R}_{i,j}=[s_{i,j},s_{i,j+1}]$, $j\in\mathbb{N}_{p_{i}-1}$, as given in (46) and (47). With $\epsilon>0$, we partition each interval $\mathcal{R}_{i,j}$ for the $i$-th user as

[TABLE]

such that $\mathcal{R}_{i,j}=\bigcup_{q=1}^{3}R^{(q)}_{i,j}$. Let $D\in\mathbb{R}$ be a scalar parameter. We modify $r_{i}(s_{i})$ to $\tilde{r}_{i}(s_{i})$ as

[TABLE]

which is represented by the piecewise affine function

[TABLE]

Denoting $r^{(m)}=\text{max}\{r^{(m)}_{i}\}$, $r^{(m)}_{i}=\text{max}_{\mathcal{S}_{i}}\{r_{i}(s_{i})\}$, if $D\geq\frac{2N^{3}(r^{(m))^{2}}}{\epsilon}$ and $0<\epsilon<\text{min}(\sigma,\nicefrac{{1}}{{N}})$, then all the local minima for $\boldsymbol{\mathcal{O}}_{5}$ are contained inside the set $\mathcal{S}$.

Proof.

The proof is stated in Appendix (Section 9.2). ∎

The piecewise affine function $\tilde{r}_{i}(s_{i})$ retain $r_{i}$ values for all $s_{i}\in\mathcal{S}_{i}$ following (54). The corresponding modified intervals for the $i$-th user $\mathcal{R}^{(k)}_{i,j}$, $k\in\mathbb{N}_{3}$, $j\in\mathbb{N}_{p_{i}-1}$ are used to determine the optimization variables in Table 1, followed by the linear inequalities in Table 2 in the form of the modified LMIs $\tilde{G}_{i}\tilde{v}_{i}\leq\tilde{g}_{i}$.

We next use convex-concave procedure (CCP) to relax the non-convex objective function in (49a) (Lipp and Boyd,, 2016). This is done by considering any point $s_{(0)}=[s_{1_{0}},\ldots,s_{N_{0}}]^{T}\in\mathcal{S}_{(\mathfrak{s})}$, where $\mathcal{S}_{(\mathfrak{s})}\subset\mathcal{S}$ is defined as

[TABLE]

and bounding the concave part of the DC objective $-g(s)$ by its convex relaxation around this point as

[TABLE]

In that case, the relaxed objective function for $\boldsymbol{\mathcal{O}}_{5}$ can be written as

[TABLE]

Algorithm 3 shows the steps for solving $\boldsymbol{\mathcal{O}}_{5}$. The following proposition provides the conditions under which Algorithm 3 converges.

Proposition 1**.**

For a given $\mathfrak{s}\in\mathcal{S}_{M}$, Algorithm 3 iteratively solves $\boldsymbol{\mathcal{O}}_{5}$ such that

(a) the difference between $F(s_{(j)})$ and its convex upper bound $\hat{F}(s_{(j)};s_{(j-1)})$ with respect to the known $s_{(j-1)}$ is given as:

[TABLE]

where $r_{i}(s_{i_{j}})$ is the performance ratio of the $i$-th user for the $j$-th iteration.

(b) Beginning from an initial guess $s_{(0)}\in\mathcal{S}_{(\mathfrak{s})}$, Algorithm 3 will converge to a local minima $s_{(k)}\in\mathcal{S}_{(\mathfrak{s})}$ such that

[TABLE]

Proof.

The proof is given in Appendix (Section 9.3). ∎

7 Examples

To test the efficiency of Algorithm 3, simulations are carried out for $N=2$ and $N=3$ users. For each case, we assign a ranking system to each $s\in\mathcal{S}$, where $\mathcal{S}$ is the solution set of our original problem $\boldsymbol{\mathcal{O}}_{4}$. The reason for this is the reformulation of $\boldsymbol{\mathcal{O}}_{4}$ to $\boldsymbol{\mathcal{O}}_{5}$ due to which all the discrete points in $\mathcal{S}$ now become the local minima for Algorithm 3. Ranking of these local minima is needed to evaluate the performance of the algorithm. Separate ranking is created for each individual $\mathfrak{s}\in\mathcal{S}_{M}$. Increasing ranks are allotted according to the increasing $F(r(s))$.

For $N=2$ case, the state matrices for both users are randomly chosen as $A_{i}\in\mathbb{R}^{10\times 10}$ with $B_{i}=B_{w_{i}}=I_{10}$, $i=1,2$. For $N=3$ case, the three systems have $n_{1}=5$, $n_{2}=4$ and $n_{3}=5$ number of states. The state matrices for these systems are randomly generated $A_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$, with $B_{i}=B_{w_{i}}=I_{n_{i}}$. Fig. 11 and 11 show the simulation results for $N=2$ and $N=3$, respectively, where the x-axis represents the given $\mathfrak{s}$, and, the y-axis shows the absolute difference between the values of $F(r(s))$ for rank-1 and rank-2 solutions, represented by $\Delta F=|F_{\text{rank 1}}-F_{\text{rank 2}}|$. The color of the scatter plot represents the rank of the solution obtained from Algorithm 3 with ‘high’ referring to rank-1 (high rank) and ‘low’ referring to the worst-ranked solution.

For $N=2$ case, while user-1 achieves $p_{1}=41$ sparsity levels $0<2<\cdots<76<77$, user-2 achieves $p_{2}=36$ sparsity levels ranging from $0<4<\cdots<67<69$. Fig. 11 shows that Algorithm 3 achieves high ranking solutions for most of the $\mathfrak{s}$ values. However, low ranking solutions can be obtained when $\Delta F=|F_{rank1}-F_{rank2}|$ is too small due to the steps of the branch and bound method that is used to solve the MILPs. Similar trend of low ranking solutions are observed for $N=3$ as well, as shown in Fig. 11.

A three-dimensional plot of the solution set for $\mathfrak{s}=89$ in the case of $N=2$ users is shown in Fig. 12, where the blue dots denote the solution set of $\boldsymbol{\mathcal{O}}_{4}$, while the connecting blue curve creates the solution set for $\boldsymbol{\mathcal{O}}_{5}$. The global minima is represented by the red dot. The $r_{i}(s_{i})$ graphs for each user are shown in Fig. 9. In summary, this example shows the case where a total of $111$ links are divided, beginning from an initial allocation of $80$ and $31$ links to user-1 and 2 respectively, which corresponds to a high discrepancy between their performance ratios $r_{1}=1.00$ and $r_{2}=1.45$. Starting from this initial guess, Algorithm 3 reaches the global minima in the second iteration, which decreases the number of links allotted to user-1 from $80$ to $53$, and, increases the number of links of user-2 from $31$ to $58$, resulting in $r_{1}\approx r_{2}=1.05$, which ensures the lowest possible disparity in their performances.

8 Conclusion

We develop algorithms to design sparse $\mathcal{H}_{2}$ optimal controllers for an LTI system in the presence of feedback delay. Our algorithms bring out the interplay between delays, sparsity and $\mathcal{H}_{2}$ objective. Both sparsity and $\mathcal{H}_{2}$ norm are facilitated when we include the flexibility in controlling the delay as compared to when the delay is kept constant. While traditional sparsity promoting algorithms destabilize the delayed system in most cases, ours preserves stability all throughout. In the second half of the paper, we extend these algorithms to a multiple user scenario, and assign a fair number of communication links to each user with the aim of reducing the variance of their $\mathcal{H}_{2}$ performances. Numerical simulations are presented to showcase the efficiency of all algorithms. Our future research will include extension of these strategies to a CPS with multiple delays, and analyzing their sensitivity to variations in network traffic using ideas from sparse stochastic optimal control.

9 Appendix

9.1 Proof of Theorem 2

Proof.

Substituting $K=K_{o}+\Delta K$, $\tau=\tau_{o}+\Delta\tau$, $P=P_{o}+\Delta P$, $\cdots$, $S_{1}=S_{1o}+\Delta S_{1}$ in (6a) and assuming $\Delta\mathcal{F}:=\texttt{blkdiag}(B\Delta K,B\Delta K,$ $B\Delta K,B\Delta K)$, the NMI $\Phi_{\text{N}}\succ 0$ can be written as

[TABLE]

where $\Phi_{0}$ is given as

[TABLE]

and $\Phi_{i}$, $i=1,2,3,4$ are given as

[TABLE]

with $\mathcal{B}_{i}$, $i=1,2,3,4$ defined as follows:

[TABLE]

The matrices $\Phi_{1}$, $\Phi_{3}$ and $\Phi_{4}$ are symmetric Hamiltonian (Negi and Chakrabortty,, 2018, Theorem 4.3.1), while $\Phi_{2}$ is symmetric, but not Hamiltonian. Therefore, the following relations can be derived using (Negi and Chakrabortty,, 2018, Lemma 1, Theorem 2):

[TABLE]

If the NMI $\Phi_{\text{N}}=\sum_{i=0}^{4}\Phi_{i}\succ 0$, then $\sum_{i=1}^{4}\lambda_{\text{min}}\left(\Phi_{i}\right)\succ 0$ as well (Horn and Johnson,, 2012, Theorem 1.2). Therefore, (65a)-(65d) can be used to derive the relaxed convex constraint for ensuring $\Phi_{\text{N}}\succ 0$. However, the bounds in (65a)-(65d) are still bilinear. Hence, two additional convex constraints are introduced next as:

[TABLE]

where $\eta>0$ and $\gamma>0$ are scalar design variables. Substituting (66) in (65a)-(65d), one can obtain (7a) to (7d) as the final relaxed convex constraint set for $\Phi_{\text{N}}\succ 0$.

Next, $\Psi_{N}\succ 0$ is relaxed. Assuming $0<\gamma\ll 1$, which implies $\Delta{\tau}\ll 1$, we can write:

[TABLE]

Using (67), $\Psi_{\text{N}}\succ 0$ can be rewritten as:

[TABLE]

Using (66), (68e)-(68i) and the same procedure as in (65) and (66), $\Psi_{\text{N}}\succ 0$ can be equivalently written as (7e)-(7f) which imply $\Phi_{\text{N}}\succeq\epsilon I\implies\Phi_{\text{N}}\succ 0$ and $\Psi_{\text{N}}\succeq\epsilon I\implies\Psi_{\text{N}}\succ 0$. ∎

9.2 Proof of Theorem 5

Proof.

Let the solution of $\boldsymbol{\mathcal{O}}_{5}$ be denoted as $s=(s_{1},\ldots,s_{N})\in\mathbb{S}_{\mathfrak{s}}$ for a given $\mathfrak{s}\in\mathcal{S}_{M}$, i.e., $\sum_{i=1}^{N}s_{i}=\mathfrak{s}$. Denoting $r(s)=[r_{1}(s_{1}),\ldots,r_{N}(s_{N})]^{T}$, $F(r(s))=H(r(s))+\sigma h(r(s))$ is denoted as $F(s)=H(s)+\sigma h(s)$. Since $\boldsymbol{\mathcal{O}}_{5}$ is solved using $\tilde{r}(s)$, the objective function $F(\tilde{r}(s))$ is denoted as $\tilde{F}(s)$.

We define the set $\mathcal{\tilde{S}}:=\mathcal{\tilde{R}}_{1}\times\cdots\times\mathcal{\tilde{R}}_{N}$, where $\mathcal{\tilde{R}}_{i}=\bigcup_{j=1}^{p_{i}-1}\bigcup_{q=1,3}R^{(q)}_{i,j}$ for the $i$-th user. We further define its subset $\mathcal{\tilde{S}}_{(\mathfrak{s})}$ $:=\{s\in\mathcal{\tilde{S}}:$ $\sum_{i=1}^{N}s_{i}=\mathfrak{s}\}$. Note that the desired solution set $\mathcal{S}_{(\mathfrak{s})}$ defined in (43) is such that $\mathcal{S}_{(\mathfrak{s})}\subset\mathcal{\tilde{S}}_{(\mathfrak{s})}\subset\mathbb{S}_{\mathfrak{s}}$. Since $\boldsymbol{\mathcal{O}}_{5}$ is evaluated over $\mathbb{S}_{\mathfrak{s}}$, we first prove that with the given bounds on $D$ and $\epsilon$, the local minima of $\tilde{F}(s)$ cannot belong to $\mathbb{S}_{\mathfrak{s}}/\mathcal{\tilde{S}}_{(\mathfrak{s})}$. We further prove that the local minima only belong to the subset $\mathcal{S}_{(\mathfrak{s})}$ of $\mathcal{\tilde{S}}_{(\mathfrak{s})}$. For any $s\in\mathbb{S}_{\mathfrak{s}}$, $\tilde{F}(s)$ is defined as:

[TABLE]

where $s_{i}=s^{*}_{i}\pm\epsilon_{i}$ is an $\epsilon_{i}>0$ perturbation on the nearest sparsity level $s^{*}_{i}\in\mathcal{S}_{i}$ such that $\sum_{i=1}^{N}s^{*}_{i}\pm\epsilon_{i}=\mathfrak{s}$. The modified curve function of (54) can be written as:

[TABLE]

where $\bar{\epsilon}_{i}=({\epsilon_{i}}/{\epsilon})\in[0,1]$ $\forall\ i$. If $s\in\mathbb{S}_{\mathfrak{s}}/\mathcal{\tilde{S}}_{(\mathfrak{s})}$, then for $1\leq k\leq N$ users, $\bar{\epsilon}_{i}=1$; for the remaining $N-K$ users $\bar{\epsilon}_{i}\in[0,\epsilon]$. The objective function value for $s\in\mathbb{S}_{\mathfrak{s}}/\mathcal{\tilde{S}}_{(\mathfrak{s})}$ is evaluated as (dropping arguments from $\tilde{r}_{i}$ and $r_{i}$):

[TABLE]

Using (70), the given bounds on $D$ and $\epsilon$, we obtain:

[TABLE]

which shows that local minima cannot belong to $\mathbb{S}_{\mathfrak{s}}/\mathcal{\tilde{S}}_{(\mathfrak{s})}$. We next evaluate the objective over the set ${\tilde{S}}_{(\mathfrak{s})}$. For any $s\in\mathcal{\tilde{S}}_{(\mathfrak{s})}$, $s_{i}$ is an $\epsilon_{i}\in[0,\epsilon]$ perturbation on the nearest sparsity level $s^{*}_{i}\in\mathcal{S}_{(\mathfrak{s})}$, i.e., $s_{i}=s^{*}_{i}\pm\epsilon_{i}$ $\forall$ $i$, where $s$ is such that $\sum_{i=1}^{N}s_{i}=\mathfrak{s}$. It can be shown easily that $\|s-s^{*}\|<\sqrt{N}\epsilon$ $\forall$ $s\in\mathcal{\tilde{S}}_{(\mathfrak{s})}$. The modified ratio curves $\tilde{r}_{i}(s_{i})$ follow (70) with $\bar{\epsilon}_{i}=\nicefrac{{\epsilon_{i}}}{{\epsilon}}\in[0,1]$ $\forall$ $i$. To show that the local minima will only exist in the subset $\mathcal{S}_{(\mathfrak{s})}$ of $\mathcal{\tilde{S}}_{(\mathfrak{s})}$, we next prove that $\tilde{F}(s)\geq\tilde{F}(s^{*})$ for $s\in\mathcal{\tilde{S}}_{(\mathfrak{s})}$.

Since $\epsilon_{i}$ follow $0<\epsilon_{i}<\frac{1}{N}$ $\forall$ $i$, the perturbations are such that $\sum_{i=1}^{N}\pm\epsilon_{i}=0$ to ensure $\sum_{i=1}^{N}s=\mathfrak{s}$. Using (71), $\tilde{F}(s)$ over the set $\mathcal{\tilde{S}}_{(\mathfrak{s})}$ is evaluated as $\tilde{F}(s)|_{s\in\mathcal{\tilde{S}}_{(\mathfrak{s})}}$ $=\tilde{F}(s^{*})+$ $\Delta\tilde{F}(s)|_{s\in\mathcal{\tilde{S}}_{(\mathfrak{s})}}$, where $\Delta\tilde{F}(s)|_{s\in\mathcal{\tilde{S}}_{(\mathfrak{s})}}$ follows the condition:

[TABLE]

where $\Delta H(s)>0$. Moreover, it can be shown easily that as $\epsilon_{i}\to 0$ $\forall\ i$, $\Delta\tilde{F}(s)|_{s\in\mathcal{\tilde{S}}_{(\mathfrak{s})}}\to 0$. Using (72) and (73), there exists a positive $\varepsilon>0$ such that if $\|s-s^{*}\|<\varepsilon$, $\tilde{F}(s)\geq\tilde{F}(s^{*})$ for all $s\in\mathbb{S}_{\mathfrak{s}}$ and $s^{*}\in\mathcal{S}_{(\mathfrak{s})}$. This proves that all the local minima belong to $\mathcal{S}_{(\mathfrak{s})}$. ∎

9.3 Proof of Proposition 1

Proof.

Let $r_{(j)}=[r_{1_{j}},\ldots,r_{N_{j}}]^{T}$ represent the vector of performance ratios obtained in the $j$-th iteration. Let $F_{j}=f_{j}-g_{j}+\sigma h_{j}$ represent the objective value for $s_{(j)}$ and $\hat{F}_{j,j-1}=f_{j}-g_{j-1}-\nabla g^{T}_{j-1}(r_{(j)}-r_{(j-1)})+\sigma h_{j}$ represent the relaxed objective value for $s_{(j)}$ with respect to $s_{(j-1)}$. It follows that $\hat{F}_{j,j-1}-F_{j}$ is given as:

[TABLE]

We next derive the convergence proof. Let $y_{j-1}=f_{j-1}-g_{j-1}+\sigma h_{j-1}=F_{j-1}$ for the $(j-1)$-th iteration. Then, $y_{j-1}=F_{j-1}=\hat{F}_{j-1,j-1}\geq\hat{F}_{j,j-1}$ as CCP minimizes $\hat{F}(j,j-1)$ at every iteration $j$. Therefore, one can write $y_{j-1}\geq f_{j}-\hat{g}_{j,j-1}\geq y_{j}$. Thus, a non-increasing sequence $\{y_{j}\}$ is obtained, which converges to a local minima of $\boldsymbol{\mathcal{O}}_{5}$. Moreover, since $\hat{F}_{j,j-1}-F_{j}=\delta_{j,j-1}$, where $\delta_{j,j-1}\geq 0$ $\forall$ $j$, beginning from the [math]-th iteration, we have $\hat{F}_{1,0}=F_{1}+\delta_{1,0},\cdots,\hat{F}_{k,k-1}=F_{k}+\delta_{k,k-1}$. Also, since $F_{j}\geq\hat{F}_{j,j-1}$ for any $j$-th iteration, it follows that $F_{0}\geq F_{1}+\delta_{1,0},\cdots,F_{k-1}\geq F_{k}+\delta_{k,k-1}$. Since $F_{j}\geq 0$ $\forall$ $j$, we directly obtain (60). ∎