Network Design for Controllability Metrics

TL;DR

This paper develops methods for tuning edge weights in fixed-topology networks to satisfy controllability metrics, using convex relaxations and optimization, with applications demonstrated on power systems.

Contribution

It introduces a convex relaxation approach for controllability-based edge weight tuning and proposes a sparsity-promoting cost function for network design.

Findings

Convex relaxations effectively solve controllability feasibility problems.

Sparsity-promoting costs reduce the number of modified edges.

Numerical simulations validate the proposed methods on power system models.

Abstract

In this paper, we consider the problem of tuning the edge weights of a networked system described by linear time-invariant dynamics. We assume that the topology of the underlying network is fixed and that the set of feasible edge weights is a given polytope. In this setting, we first consider a feasibility problem consisting of tuning the edge weights such that certain controllability properties are satisfied. The particular controllability properties under consideration are (i) a lower bound on the smallest eigenvalue of the controllability Gramian, which is related to the worst-case energy needed to control the system, and (ii) an upper bound on the trace of the Gramian inverse, which is related to the average control energy. In both cases, the edge-tuning problem can be stated as a feasibility problem involving bilinear matrix equalities, which we approach using a sequence of convex relaxations. Furthermore, we also address a design problem consisting of finding edge weights able to satisfy the aforementioned controllability constraints while seeking to minimize a cost function of the edge weights, which we assume to be convex. In particular, we consider a sparsity-promoting cost function aiming to penalize the number of edges whose weights are modified. Finally, we verify our results with numerical simulations over many random network realizations as well as with an IEEE 14-bus power system topology.

Figures

Equations

$$x(k+1)=A(\mathcal{G})x(k)+Bu(k),$$

$$J(T,x_T):=x_T^{{}^{⊺}}(W_{r,T}{)}^{-1}{x}_T,$$

$$A(\mathcal{G})W_r^{\infty }A(\mathcal{G}{)}^{{}^{⊺}}-W_r^{\infty }+B{B}^{{}^{⊺}}=0.$$

$$W_r^{\infty }-\tilde{\lambda }{I}_n⪰0.$$

$$n\tilde{\tau }-\mathrm{tr}\{(W_r^{\infty }{)}^{-1}\}\ge 0,$$

$$W_r^{\infty }\in {\mathcal{W}}_{\theta },$$

$$x(k+1)=[A(\mathcal{G})+\Delta (\mathcal{G})]x(t)+Bu(t).$$

$$\mathrm{find}$$

$$\mathrm{subjectto}$$

$$\Delta \in \mathcal{D},$$

$$(A+\Delta )$$

$$|{\lambda }_i(A+\Delta )|<1,i=1,\dots ,n,$$

$$\underset{\begin{array}{c}\Delta \in {\mathbb{R}}^{n\times n}\\ W\in {\mathbb{S}}_{++}^n\end{array}}{\mathrm{minimize}}$$

$$\mathrm{subjectto}$$

$$M(W,H)N(\Delta )=Q,$$

$$M(W,H):=\left[\begin{array}{cc}{H}^{{}^{⊺}}& -W\\ -W& H\end{array}\right],N(\Delta ):=\left[\begin{array}{c}(A+\Delta {)}^{{}^{⊺}}\\ I_n\end{array}\right],Q:=\left[\begin{array}{c}-B{B}^{{}^{⊺}}\\ 0\end{array}\right].$$

$$(A+\Delta )H-W$$

$$H-W(A+\Delta {)}^{{}^{⊺}}$$

$$\mathcal{Z}(W,H,\Delta )$$

$$=\left[\begin{array}{ccc}{I}_n& 0& (A+\Delta {)}^{{}^{⊺}}\\ 0& I_n& I_n\\ H^{{}^{⊺}}& -W& -B{B}^{{}^{⊺}}\\ -W& H& 0\end{array}\right].$$

$$\mathrm{rank}[Z]=\mathrm{rank}[Z_{11}]+\mathrm{rank}[Z/Z_{11}].$$

$${\eta }_r(X):=\sum _{i=r+1}^{\min \{m,n\}}{\sigma }_i(X),\text{vspace}-3mm$$

$${\eta }_r(X)=\Vert X{\Vert }_{\ast }-{\Vert X\Vert }_{\lceil r\rceil },$$

$${\eta }_r(X)=\Vert X{\Vert }_{\ast }-\underset{\begin{array}{c}L{L}^{{}^{⊺}}=I_r\\ R{R}^{{}^{⊺}}=I_r\end{array}}{\sup }\mathrm{tr}\{LX{R}^{{}^{⊺}}\},$$

$$\underset{W,H,\Delta }{\mathrm{minimize}}$$

$$\mathrm{subjectto}$$

$$\underset{W,H,\Delta }{\mathrm{minimize}}$$

$$\mathrm{subjectto}$$

$$\underset{W,H,\Delta }{\mathrm{minimize}}$$

$$\mathrm{subjectto}$$

$$=\underset{W,H,\Delta }{\mathrm{minimize}}$$

$$-\underset{\begin{array}{c}L{L}^{{}^{⊺}}=I_{2n}\\ R{R}^{{}^{⊺}}=I_{2n}\end{array}}{\sup }\mathrm{tr}\{L\mathcal{Z}(W,H,\Delta )R^{{}^{⊺}}\}$$

$$\mathrm{subjectto}$$

$$\underset{W,H,\Delta }{\mathrm{minimize}}$$

$$\mathrm{subjectto}$$

$$A=$$

$$\displaystyle\left[\setcounter{MaxMatrixCols}{14}\begin{smallmatrix}\,\,\cdot\,\,&0.06&\,\,\cdot\,\,&\,\,\cdot\,\,&0.22&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ 0.06&\,\,\cdot\,\,&0.20&0.18&0.17&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ \,\,\cdot\,\,&0.20&\,\,\cdot\,\,&0.17&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ \,\,\cdot\,\,&0.18&0.17&\,\,\cdot\,\,&0.04&\,\,\cdot\,\,&0.21&\,\,\cdot\,\,&0.56&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ 0.22&0.17&\,\,\cdot\,\,&0.04&\,\,\cdot\,\,&0.25&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.25&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.20&0.26&0.13&\,\,\cdot\,\,\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.21&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.18&0.11&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.18&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.56&\,\,\cdot\,\,&\,\,\cdot\,\,&0.11&\,\,\cdot\,\,&\,\,\cdot\,\,&0.08&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.27\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.08&\,\,\cdot\,\,&0.19&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.20&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.19&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.26&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.20&\,\,\cdot\,\,\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.13&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.20&\,\,\cdot\,\,&0.35\\ \,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.27&\,\,\cdot\,\,&\,\,\cdot\,\,&\,\,\cdot\,\,&0.35&\,\,\cdot\,\,\\ \end{smallmatrix}\right]$$

Full text

I Introduction

Many technological, biological, chemical, and social systems can be modeled as large ensembles of dynamical units connected via an intricate pattern of interactions [1]. From an engineering perspective, we are interested in efficiently steering the dynamics of these complex systems via external actuation. In this direction, control theory provides us with the notion of controllability to decide whether a given system can be steered towards an arbitrary state [2]. Furthermore, the so-called controllability Gramian of a system, which implicitly depends on the system’s dynamics and the configuration of its actuators, can be used to quantify the energy required to steer the system, assuming the system is controllable [2]. Leveraging these notions, several papers have recently focused on the problem of optimally allocating actuators throughout the network under several performance metrics [3, 4, 5, 6, 7, 8, 9, 10, 11, 12].

In some scenarios, instead of designing the location of external actuators, one may consider the alternative problem of modifying the network’s dynamics given a fixed configuration of actuators. For example, in power systems, one can tune the electrical parameters of the transmission lines using, for example, flexible AC transmission system (FACTS) devices [13, 14]. Also, in multi-agents networks, the interactions between agents can usually be modified to achieve a particular objective [15]. For instance, in leader-follower multi-agent networks, one may consider the scenario where both the communication topology and the location of the external actuators are fixed. Then, one can seek a set of edge weights (e.g., the agents’ update rules) such that the average and/or worst-case energy required to drive the state of the network satisfies certain bounds. In this regard, the present work first considers the feasibility problem of finding the edge weights of a linear networked system such that certain bounds on controllability metrics are satisfied. Secondly, we address the design problem of finding edge weights able to satisfy the aforementioned bounds while seeking to minimize a cost function of the edge weights, which we assume to be convex. In particular, we consider a $1$-norm sparsity-promoting cost function aiming to penalize the number of edges whose weights are modified in the resulting design.

I-1 Related Work

In recent years, the problem of designing systems to satisfy certain controllability metrics has mostly focused on finding optimal actuator configurations, i.e., the location of those nodes to be externally actuated by control inputs [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. In addition, a considerable amount of research has been dedicated to understanding how the network topology impacts control performance [16, 17, 18, 19, 20, 21, 22, 23, 7, 24, 25]. In particular, [25] establishes necessary and sufficient graph-theoretical conditions for a discrete-time networked system to exhibit a diagonal controllability Gramian. In [26], the authors characterize the minimum input energy required to transfer a discrete-time dynamical system with bilinear dynamics from the origin to a desired state. The work in [27] proposes the notion of observability radius, which measures how much the parameters of a dynamical system can be perturbed before the system becomes unobservable. In a similar direction, the work in [28] investigates the effect of adding network edges to improve spectral performance metrics for the case of consensus dynamics over networks.

More generally, the works in [29, 30, 31, 32, 33] investigate design problems that seek to optimize network dynamical properties such as the dominant eigenvalue of the system matrix, with applications to virus spread and wireless control networks appearing in [34, 35, 36, 37].

The present paper extends previous work by the authors in [38] through several contributions. Specifically, in this paper we: (i) address the discrete-time case, in which the discrete Lyapunov equation introduces higher-degree products in its decision variables and requires new transformation steps for its treatment; (ii) provide an analysis of the conditions under which stability of the designed system is assured; (iii) consider cost functions over edge weights, which can be used to promote solutions with higher sparsity in edge modifications; (iv) propose a convex relaxation approach, which enables a more detailed analysis of convergence; (v) consider average controllability as an additional controllability metric; and (vi) present comprehensive computational experiments to illustrate the above aspects.

I-2 Structure and contributions of the paper

The rest of the paper is organized as follows. In Section II, we formalize both the network feasibility and the network design problems, in which we are tasked with tuning the weights of the edges in a given network in order to satisfy certain controllability metrics. Specifically, we consider two metrics: (i) the worst-case control energy, which is related to the smallest eigenvalue of the Gramian, and (ii) the average energy required to drive the system, which is related to the trace of the Gramian inverse. In Section III, we provide a detailed description of the strategy followed to solve both problems. In particular, we cast both the feasibility and the design problems into nonlinear optimization programs with quadratic bilinear terms, which are, in general, computationally hard to solve. We approach these optimization problems by lifting the space of variables and adding a rank constraint on a matrix whose entries depend affinely on the decision variables. We then propose a sequence of convex problems to relax this rank constraint using a truncated nuclear norm. In Section IV, we illustrate the validity of our results via computational experiments on random graphs, as well as a $1$-norm sparsity-promoting design problem considering the IEEE 14-bus system. We conclude and enumerate some possibilities for future work in Section V.

Notation

We denote by $[X]_{i,j}$ the entry at the $i$-th row and $j$-th column of the matrix $X\in\mathbb{R}^{m\times n}$. The transpose of X is written as $X^{{}^{\intercal}}$. The $n\times n$ identity matrix is denoted by $I_{n}$. The operator $\operatorname{diag}(a_{1},\ldots,a_{n})$ returns a diagonal matrix having $a_{1},\ldots,a_{n}$ as entries in its diagonal. The inner product between two matrices $X,Y\in\mathbb{R}^{m\times n}$ is given by $\left<X,Y\right>=\operatorname{tr}\{X^{{}^{\intercal}}Y\}$,

where $\operatorname{tr}\{X^{{}^{\intercal}}Y\}=\sum_{i=1}^{n}[X^{{}^{\intercal}}Y]_{i,i}$ denotes the trace operator.

The $1$-norm of a matrix $X\in\mathbb{R}^{m\times n}$ is defined as the $\ell_{1}$-norm of its vectorization, i.e., $\|X\|_{1}=\|\!\operatorname{vec}(X)\|_{1}$.

Likewise, the [math]-norm of a matrix is defined as the $\ell_{0}$-quasi-norm of its vectorization, i.e., the number of nonzero entries. The infinity norm of $X$ is defined as $\|X\|_{\infty}=\max_{i,j}[X]_{i,j}$.

The nuclear norm of $X$ is defined in terms of its singular values $\sigma_{i}(X)$, $i=1,\ldots,\min\{m,n\}$, as $\|X\|_{\ast}=\sum_{i=1}^{\min\{m,n\}}\sigma_{i}(X)$.

The operator norm of $X$ is denoted by $\|X\|$ and computed as $\|X\|=\sigma_{1}(X)$, the largest singular value of $X$.

We denote by $\mathbb{S}^{n}$ the set of symmetric matrices of dimension $n$. Likewise, $\mathbb{S}_{+}^{n}$ (resp., $\mathbb{S}_{++}^{n}$) is the set of symmetric positive semidefinite (resp., definite) matrices. Correspondingly, the semidefinite partial ordering is denoted $X\succeq Y$ (resp., $X\succ Y$) when $X-Y\succeq 0$ (resp., $X-Y\succ 0$).

A set $\mathcal{S}\subset\mathbb{R}^{m}$ is a spectrahedron [39, Def. 2.6] if it can be represented in the form $\mathcal{S}=\{(x_{1},\ldots,x_{m})\in\mathbb{R}^{m}:Q_{0}+\sum_{i=1}^{m}Q_{i}x_{i}\succeq 0\}$, for $Q_{0},\ldots,Q_{m}\in\mathbb{S}^{n}$.

A proper algebraic variety $\mathcal{V}\subset\mathbb{R}^{n}$ is the set of common zeros of a finite number of nonzero polynomials in $n$ variables.

II Problem Formulation

Consider a networked system following a discrete-time linear time-invariant dynamics, described by

[TABLE]

where $x(k)\in\mathbb{R}^{n}$ denotes the vector of states and $u(k)\in\mathbb{R}^{m}$ is the vector of inputs at instant $k$. The sparsity pattern of the state matrix $A(\mathcal{G})\in\mathbb{R}^{n\times n}$ is constrained by a directed interdependency graph $\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$ defined by a set of nodes $\mathcal{V}=\{1,\ldots,n\}$ and a set of edges $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$, such that $[A(\mathcal{G})]_{i,j}\in\mathbb{R}$ if the edge $(j,i)\in\mathcal{E}$, and $[A(\mathcal{G})]_{i,j}=0$ if $(j,i)\notin\mathcal{E}$. Also, the input matrix $B\in\mathbb{R}^{n\times m}$ is such that $[B]_{i,l}\neq 0$ if the external input signal $[u(k)]_{l}$ directly influences $[x(k+1)]_{i}$, and $[B]_{i,l}=0$ otherwise.

Next, consider the problem of driving the state of the network from a given initial state $x_{0}\equiv x(0)$ to a desired target state $x_{T}\equiv x(T)$ within a time horizon $T>0$, by designing a sequence of inputs $u(k)$ for $k\in\{0,1,\ldots,T-1\}$. If any $x_{T}\in\mathbb{R}^{n}$ is attainable from $x_{0}=0_{n}$ within a time horizon  $T$, then the system (1) is said to be reachable, which we refer to $(A(\mathcal{G}),B)$ being reachable. Furthermore, it is known that the minimum input control energy to steer the system to a desired final state $x_{T}$ from $x_{0}=0_{n}$ is given by [2]

[TABLE]

where ${W_{r,T}}$ is called the finite-horizon reachability Gramian, defined as ${W_{r,T}}\coloneqq\sum_{k=0}^{T-1}A(\mathcal{G})^{k}BB^{{}^{\intercal}}(A(\mathcal{G})^{{}^{\intercal}})^{k}.$ The infinite-horizon reachability Gramian is then obtained as the limit $W_{r}^{\infty}\coloneqq\lim_{T\rightarrow\infty}{W_{r,T}}$. This Gramian is positive definite, and can be computed as the (unique) solution to the discrete-time Lyapunov equation

[TABLE]

when the system is reachable and $A(\mathcal{G})$ is stable [2].

II-A Reachability Metrics

We focus on two metrics related to the reachability Gramian to quantify the minimum input energy to drive the system [40, 16, 26].

Worst-case minimum input energy

Because $W_{r}^{\infty}$ is (symmetric) positive definite when the system is reachable, its eigenvalues $\lambda_{1}\leq\ldots\leq\lambda_{n}$ are positive real numbers, with corresponding eigenvectors $v_{i}$ for $i=1,\ldots,n$. It turns out that the final state $x_{T}$ satisfying $\|x_{T}\|_{2}=1$ requiring the largest minimum input energy to be reached from $x_{0}=0_{n}$ is given by the (normalized) eigenvector $v_{1}$. The energy required to drive the state from the origin towards $v_{1}$ within an infinite horizon is equal to $\lambda_{1}^{-1}$, which we call the worst-case minimum input energy. Therefore, if we require the worst-case minimum input energy to be less than or equal to a desired value $\tilde{\lambda}^{-1}{\,>0}$, then the reachability Gramian must satisfy the following semidefinite constraint:

[TABLE]

Average minimum input energy

The expected energy required to steer the system from the origin towards a random final state uniformly distributed over the unit sphere is equal to $\frac{1}{n}\operatorname{tr}\{(W_{r}^{\infty})^{-1}\}$ [40], which we call the average minimum input energy. In a manner similar to the worst-case minimum input energy metric, we can constrain the average minimum input energy to be upper-bounded by a target value $\tilde{\tau}{\;<\infty}$ via the condition

[TABLE]

which is also representable by a semidefinite constraint over $W_{r}^{\infty}$ (see Lemma A.3 in the Appendix).

In what follows, we will refer to the aforementioned reachability constraints on $W_{r}^{\infty}$ by the set membership condition

[TABLE]

where $\mathcal{W}_{\theta}$ is a convex set (more precisely, a spectrahedron) defined by constraints (4) and/or (5), and indexed by the parameters in $\theta=(\tilde{\lambda},\tilde{\tau})$.

II-B Network Design for Reachability

As previously mentioned, we consider the problem of tuning the edge weights of a given network in order to satisfy certain minimum control energy requirements (either in worst-case or in average). In particular, we assume that we are able to add a matrix $\Delta(\mathcal{G})\in\mathbb{R}^{n\times n}$ to the state matrix $A(\mathcal{G})$, such that $\Delta(\mathcal{G})$ presents the same sparsity pattern as the interdependency graph, i.e., $[\Delta(\mathcal{G})]_{i,j}=0$ for $(j,i)\notin\mathcal{E}$. After this addition, the dynamics of the network becomes

[TABLE]

Furthermore, we may require that $\Delta(\mathcal{G})$ be contained in a given polytope $\mathcal{D}$ encoding acceptable limits for its entries. For example, we can impose upper and lower bounds of the form $[\Delta(\mathcal{G})]_{i,j}\in[\iota_{i,j},\upsilon_{i,j}]$ for $(j,i)\in\mathcal{E}$ in the design problem. Subsequently, we consider the model described by (7) and address the following two problems111For compactness of notation, we will denote $A(\mathcal{G})$, $W_{r}^{\infty}$, and $\Delta(\mathcal{G})$ simply by $A$, $W$, and $\Delta$, respectively, in the rest of the paper..

II-B1 Feasible Design for Reachability Metrics

We seek an addition $\Delta\in\mathcal{D}$ such that the resulting reachability Gramian $W\in\mathbb{S}_{++}^{n}$ satisfies $W\in\mathcal{W}_{\theta}$. This can be posed as the following feasibility problem:

$\operatorname{\mathcal{P}_{1}}$Feasible Design for Reachability Metrics:

Given the interdependency graph $\mathcal{G}$, with $(A,B)$ reachable, we would like to

[TABLE]

where constraint (10) arises from the discrete-time Lyapunov equation associated with (7), and constraint (11) enforces the stability of the designed system.

Remark 1*:*

Partial design, allowing only a subset of the edge weights to be modified, can be performed by imposing additional constraints $[\Delta]_{i,j}=0$ for the edges $(j,i)$ that cannot be affected by the design procedure.

As we will show in the next section, this feasibility problem can be addressed using a sequence of convex relaxations. This problem also lays the foundation to our second problem, described next.

II-B2 Design for Reachability with Structural Penalties

In this formulation, we introduce an optimization objective that penalizes entries of $\Delta$ with large magnitudes, while meeting the reachability requirements on $W$ and structural constraints on $\Delta$. In particular, aiming at penalizing the number of edges modified, we consider the $1$-norm penalty over the entries of $\Delta$ as our cost function.

The $1$-norm behaves as a convex envelope to the [math]-norm (i.e., the number of non-zero entries in the matrix), and has found wide use in the signal processing and optimization literature [41, 42, 43]. In control systems problems, it has been successfully applied to promote sparsity in control architectures, for instance, in [44, 45].

$\operatorname{\mathcal{P}_{2}}$Design for Reachability with Structural Penalties:

Given an interdependency graph $\mathcal{G}$ and a reachable system $(A,B)$, find a structural addition $\Delta$ seeking to

[TABLE]

As will be described in Section III-D, this problem can be addressed by a sequence of convex relaxations involving an additive penalty term over the 1-norm of $\Delta$, whose limiting value is obtained by a procedure called regularization path [46].

Remark 2*:*

More generally, in $\operatorname{\mathcal{P}_{2}}$, we could consider a cost function having individual weights over the entries of $\Delta$. For simplicity, in this paper we consider all entries to have unit weight.

III Design for Reachability Algorithm

In this section, we propose a computational procedure to address $\operatorname{\mathcal{P}_{1}}$ and $\operatorname{\mathcal{P}_{2}}$. We begin by providing preliminary analyses of the Lyapunov equation (10) and of the stability constraint (11). We show that the Lyapunov equation constraint can be transformed into a rank constraint, and that its solution will imply the stability of $A+\Delta$ almost surely. Then, we solve $\operatorname{\mathcal{P}_{1}}$ by handling the rank constraint through a sequence of convex problems with guaranteed convergence. Subsequently, we address $\operatorname{\mathcal{P}_{2}}$ by computing a regularization path over a weight parameter that controls the sparsity of the generated solutions.

III-A Stability from a positive solution to the Lyapunov Equation

In this section, we show that constraint (11) is satisfied almost surely by all $\Delta\in\mathcal{D}$ that satisfy the Lyapunov constraint in (10). Following methodologies similar to [47, 48, 49, 50], we formalize this result in the next theorem.

Theorem 1* (Stability of the designed system):*

For a solution $(W,\Delta)$ to (10) with $W\succ 0$, if the original system $(A,B)$ is reachable, then the system $A+\Delta$ will be stable for any $\Delta\in\mathcal{D}\setminus\mathcal{V}$, where $\mathcal{V}$ is a set with Lebesgue measure zero.

Proof.

Applying Lemma A.1 from the Appendix for the matrix $A+\Delta$, we have that a solution $W$ to (10) exists and is unique for all $\Delta\in\mathcal{D}\setminus\mathcal{V}_{0}$, where $\mathcal{V}_{0}$ is a proper algebraic variety with Lebesgue measure zero. Further, since the pair $(A,B)$ is reachable and $\Delta$ is restricted to the structure of $A$ by $\mathcal{D}$, from [48, Proposition 2], the pair $(A+\Delta,B)$ is also reachable for $\Delta\in\mathcal{D}\setminus\mathcal{V}_{1}$, where $\mathcal{V}_{1}$ is a proper algebraic variety with Lebesgue measure zero. Therefore, since a finite union of proper algebraic varieties is a proper algebraic variety, we have that the system $A+\Delta$ will be reachable and will have a unique solution $W\succ 0$ to (10) for any $\Delta\in\mathcal{D}\setminus\mathcal{V}$, where $\mathcal{V}\coloneqq\mathcal{V}_{0}\cup\mathcal{V}_{1}$ is a proper algebraic variety with zero Lebesgue measure. Thus, applying Lemma A.2, we have that $A+\Delta$ will be stable for all $\Delta\in\mathcal{D}\setminus\mathcal{V}$. ∎

Therefore, seeking a tractable computational strategy for $\operatorname{\mathcal{P}_{1}}$, we consider constraint (11) to be implicitly satisfied by all points satisfying (8) and (10) which do not lie in $\mathcal{V}$. Consequently, if the solution to $\operatorname{\mathcal{P}_{1}}$, as determined by specific constraint sets $\mathcal{W}_{\theta}$ and $\mathcal{D}$, is such that $\Delta\in\mathcal{V}$, then, we declare $\operatorname{\mathcal{P}_{1}}$ to be infeasible for the parameters defining those sets. The same considerations apply to $\operatorname{\mathcal{P}_{2}}$.

III-B Discrete-time Lyapunov Equation as a Rank Condition

Notice that, for both problems $\operatorname{\mathcal{P}_{1}}$ and $\operatorname{\mathcal{P}_{2}}$, the discrete-time Lyapunov constraint (10) induces double and triple products between the decision matrices $\Delta$ and $W$. To address this issue, we first show that (10) can be alternatively satisfied by the solution of a lifted bilinear matrix equation (BME). Then, we approximate the solution of the resulting BME-constrained problem using a sequence of convex problems. We begin by lifting the constraint in (10) into a BME using the following lemma.

Lemma 1*:*

The discrete-time Lyapunov equation (10) is satisfied by $W$ and $\Delta$ when the following BME is satisfied by the variables $W\in\mathbb{S}_{++}^{n}$, $H\in\mathbb{R}^{n\times n}$, and $\Delta\in\mathbb{R}^{n\times n}$:

[TABLE]

where

[TABLE]

Proof.

The equation in (12) is equivalent to the following system of matrix equations:

[TABLE]

From (13b), we have that $H=W(A+\Delta)^{{}^{\intercal}}$. Substituting this $H$ in (13a), we obtain the Lyapunov equation in (10), as desired. ∎

We now rewrite the BME in (12) as an equivalent rank constraint over a matrix with a specific block structure, as stated in the next theorem.

Theorem 2* (Rank condition for Lyapunov equation):*

Let $\mathcal{Z}(W,H,\Delta)\in\mathbb{R}^{4n\times 3n}$ be the structured matrix defined as

[TABLE]

If $\operatorname{rank}[\mathcal{Z}(W^{\star},H^{\star},\Delta^{\star})]=2n$, then $W^{\star}$ and $\Delta^{\star}$ satisfy the discrete-time Lyapunov equation in (10).

Proof.

Consider the Schur complement of $Z_{11}$ in $Z\equiv\mathcal{Z}(W^{\star},H^{\star},\Delta^{\star})$, given by $Z/Z_{11}=Z_{22}-Z_{21}Z_{11}^{-1}Z_{12}$. From (14), we have that $Z/Z_{11}=Q-M^{\star}N^{\star}$, where $M^{\star}\equiv M(W^{\star},H^{\star})$ and $N^{\star}\equiv N(\Delta^{\star})$. According to Guttman’s rank additivity formula [51], the following holds:

[TABLE]

Since $\operatorname{rank}(Z_{11})=2n$, we have that $\operatorname{rank}(Z)=2n$ if and only if $\operatorname{rank}[Z/Z_{11}]=0=\operatorname{rank}[Q-M^{\star}N^{\star}]$, or equivalently, $Q=M^{\star}N^{\star}$. Thus, by Lemma 1, it follows that $W^{\star}$ and $\Delta^{\star}$ satisfy the discrete-time Lyapunov equation in (10). ∎

Equipped with the above result, we can replace the constraint in (10) by the rank constraint $\operatorname{rank}[\mathcal{Z}(W,H,\Delta)]=2n$ in both problems $\operatorname{\mathcal{P}_{1}}$ and $\operatorname{\mathcal{P}_{2}}$. Importantly, notice that the blocks of $\mathcal{Z}(W,H,\Delta)$ depend affinely on the problem decision matrices $W$ and $\Delta$. Next, we show that this reformulation can be approached using a sequence of convex programs.

III-C Design for Reachability via Sequential Optimization

As introduced in Theorem 2, a solution $(W^{\star},\Delta^{\star})$ to (7) will be obtained when the rank of $\mathcal{Z}(W^{\star},H^{\star},\Delta^{\star})$ equals $2n$. To achieve this condition, one would in principle seek to minimize the rank of $\mathcal{Z}(W,H,\Delta)$, which is a non-convex and discontinuous function. Alternatively, problems having the rank as an objective function have been approached by considering the nuclear norm (i.e., the sum of a matrix’s singular values) as a relaxation[42]. Further, from Theorem 2, we have a-priori information on the specific optimal value (equal to $2n$) for the rank of $Z$. In this case, alternative functions related to the nuclear norm have been shown to produce better approximations to the rank function [52]. In particular, the truncated nuclear norm function, defined next, uses the rank as an index restricting the number of (ordered) singular values considered in its computation.

Definition 1* (Truncated nuclear norm function):*

The truncated nuclear norm function (TNN) of a matrix $X\in\mathbb{R}^{m\times n}$ with respect to an integer parameter $r$ satisfying $r<\min\{m,n\}$ is defined as

[TABLE]

where $\sigma_{i}$ takes values over the set of singular values of $X$ sorted in descending order.

Using this definition, we can re-state the conditions in Theorem 2 in terms of the TNN, as described below.

Corollary 1* (TNN sufficient condition for Lyapunov equation):*

If the tuple $(W^{\star}\in\mathbb{S}_{++},H^{\star}\in\mathbb{R}^{n\times n},\Delta^{\star}\in\mathbb{R}^{n\times n})$ satisfies $\eta_{2n}(\mathcal{Z}(W^{\star},H^{\star},\Delta^{\star}))=0$, then $(W^{\star},\Delta^{\star})$ satisfies the discrete-time Lyapunov equation (10).

Proof.

The value $\eta_{2n}(\mathcal{Z}(W^{\star},H^{\star},\Delta^{\star}))=0$ implies $\sigma_{i}=0$ for $i=2n+1,\ldots,{3n}$. This, in turn, implies that $\operatorname{rank}[\mathcal{Z}(W^{\star},H^{\star},\Delta^{\star}]=2n$ in (14), and subsequently (10) is satisfied by invoking Theorem 2. ∎

The next lemma establishes a useful fact associated with Definition 1.

Lemma 2* (TNN via Von Neumann’s inequality [52]):*

Let $\left\|X\right\|_{\left\lceil r\right\rceil}\coloneqq\sum_{i=1}^{r}\sigma_{i}(X)$ denote the Ky Fan norm of a matrix $X\in\mathbb{R}^{m\times n}$ with respect to an integer $r\leq\min\{m,n\}$. Then, the TNN can be written as

[TABLE]

which is a difference-of-convex function of $X$. Moreover, the TNN is equivalently given by

[TABLE]

for $L\in\mathbb{R}^{r\times m}$ and $R\in\mathbb{R}^{r\times n}$.

Proof.

We have $\|X\|_{\ast}-\left\|X\right\|_{\left\lceil r\right\rceil}=\sum_{i=1}^{\min\{m,n\}}\sigma_{i}(X)-\sum_{i=1}^{r}\sigma_{i}(X)=\sum_{i=r+1}^{\min\{m,n\}}\sigma_{i}(X)=\eta_{r}(X)$. This form is clearly a difference of convex functions, since it is a difference between the nuclear and Ky Fan norms of $X$. Equation (16) is proved by observing the equivalence of $\left\|X\right\|_{\left\lceil r\right\rceil}$ with ${\sup_{LL^{{}^{\intercal}}=I_{r},RR^{{}^{\intercal}}=I_{r}}}\operatorname{tr}\{LXR^{{}^{\intercal}}\}$, as established by Lemma A.4 in the Appendix. The supremum term is defined over a family of affine functions parameterized by the matrices $L$ and $R$; hence, it is convex. ∎

Using Corollary 1, we can reformulate $\operatorname{\mathcal{P}_{1}}$ by seeking to minimize $\eta_{2n}(\mathcal{Z}(W,H,\Delta))$ subject to the reachability requirements in (8) and structural constraints in (9). Using Lemma 2, a solution to $\operatorname{\mathcal{P}_{1}}$ can be found by solving the following problem.

$\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$Difference-of-norms problem:

[TABLE]

As established in Theorem 1, a solution to $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ will fulfill the stability constraint in (11) almost surely. Further, despite its non-convexity, $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ has a known global optimal value when $\operatorname{\mathcal{P}_{1}}$ is feasible. From Corollary 1, this optimal value is equal to $\eta_{2n}(\mathcal{Z}(W,H,\Delta))=0$.

Next, taking inspiration from related problems in the literature [52], we employ a specific strategy consisting of solving a sequence of convex problems. More specifically, a convex relaxation of $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ is obtained by replacing the supremum over parameters $L$ and $R$ in (16) by fixed values $\check{L}$ and $\check{R}$, respectively, as formalized next.

$\operatorname{\mathcal{P}_{1-\mathrm{SUB}}}$Convex sub-problem for $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$:

For fixed $\check{L}\in\mathbb{R}^{2n\times 4n}$ and $\check{R}\in\mathbb{R}^{2n\times 3n}$, we define the convex problem $\mathcal{C}(\check{L},\check{R};\theta)$ as

[TABLE]

Subsequently, using Von Neumann’s trace inequality in Lemma A.4, a sequence of convex problems can be defined by iteratively solving $\operatorname{\mathcal{P}_{1-\mathrm{SUB}}}$ according to the following rule: At each iteration $k$, the parameters $L^{(k)}$ and $R^{(k)}$ are fixed, and convex sub-problem $\mathcal{C}(L^{(k)},R^{(k)};\theta)$ is solved. Then, the left- and right-singular vectors of the current solution $\mathcal{Z}^{(k)}(W,H,\Delta)=\operatorname{argmin}_{W,H,\Delta}\mathcal{C}(L^{(k)},R^{(k)};\theta)$ are used, respectively, to update parameters $L^{(k+1)}$ and $R^{(k+1)}$ for the next iteration. Such procedure, summarized in Algorithm 1, generates a monotonically convergent sequence of objective function values, as shown in the next theorem.

Theorem 3* (Convergence of Algorithm 1):*

Let $\alpha_{k}\coloneqq\eta_{2n}(\mathcal{Z}(W^{(k)},H^{(k)},\Delta^{(k)}))$. Then, the sequence $\{\alpha_{k}\}$ generated by $(W^{(k)},H^{(k)},\Delta^{(k)})=\operatorname{argmin}\,\mathcal{C}(L^{(k)},R^{(k)};\theta)$, according to Algorithm 1, is monotonically non-increasing.

Proof.

We assume that the sets $\mathcal{D}$ and $\mathcal{W}_{\theta}$ are non-empty, i.e., there exists at least one feasible solution $(W^{(0)},H^{(0)},\Delta^{(0)})$ to the relaxed problem $\mathcal{C}(L^{(0)},R^{(0)};\theta)$. For example, for the worst-case minimum energy design, a feasible solution can be constructed by letting any $\Delta^{(0)}\in\mathcal{D}$, $W^{(0)}=\tilde{\lambda}I_{n}$, and $H^{(0)}=W^{(0)}(A+\Delta)^{{}^{\intercal}}$. Because Step A (in Algorithm 1) does not affect feasibility of the initial feasible solution $(W^{(0)},H^{(0)},\Delta^{(0)})$, this solution will remain feasible for Step B, which will also retain feasibility, by construction. Therefore, a solution $(W^{(k)},H^{(k)},\Delta^{(k)})$ will remain feasible at any iteration $k$. Let $\phi(Z,L,R)\coloneqq\|\mathcal{Z}(W,H,\Delta)\|_{\ast}-\operatorname{tr}\{L\,\mathcal{Z}(W,H,\Delta)\,R^{{}^{\intercal}}\}$ be the value of the objective function of $\mathcal{C}(L,R;\theta)$ evaluated at $Z$, for $Z\equiv\mathcal{Z}(W,H,\Delta)$. We now analyze the behavior of the objective function at any iteration $k$. Denote by $p_{A}^{(k)}\coloneqq\phi(Z^{(k)},L^{(k)},R^{(k)})$ the objective function value returned after execution of Step A in Algorithm 1. Likewise, denote by $p_{B}^{(k)}\coloneqq\phi(Z^{(k+1)},L^{(k)},R^{(k)})$ the objective function value returned after execution of Step B. Because Step B involves the solution of a (feasible) convex optimization problem, we have $p_{B}^{(k)}\leq p_{A}^{(k)}$. Further, by invoking Lemma 2, we have that $p_{A}^{(k+1)}\leq p_{B}^{(k)}$. Therefore, we have $p_{A}^{(k+1)}\leq p_{A}^{(k)}$ for any $k$, and $\alpha_{k}=p_{A}^{(k)}$. Thus, for any $\epsilon_{\eta}>0$, there exists an iteration number $k$ such that $|\alpha_{k+1}-\alpha_{k}|\leq\epsilon_{\eta}$, and the sequence $\{\alpha_{k}\}$ is monotonically non-increasing. ∎

III-D Design for Reachability with Structural Penalties

We now build on the results obtained for the feasibility problem $\operatorname{\mathcal{P}_{1}}$ to address the more challenging problem $\operatorname{\mathcal{P}_{2}}$, which seeks to penalize large magnitudes in the entries of $\Delta$. First, we observe that using the definition of the truncated nuclear norm introduced in the previous section, $\operatorname{\mathcal{P}_{2}}$ can be approximated by solving the following problem for increasing values of the positive weight $\gamma$.

$\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$Penalized difference-of-norms problem:

For $\gamma$ a positive scalar, a relaxation of $\operatorname{\mathcal{P}_{2}}$ can be written as

[TABLE]

[TABLE]

where we have removed the explicit stability constraint (11) based on the results presented in Theorem 1. Besides using a relaxation strategy similar to the one previously used for $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ (i.e., replacing the supremum operator with fixed values for $L$ and $R$), we associate with $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ the following convex sub-problem.

$\operatorname{\mathcal{P}_{2-\mathrm{SUB}}}$Convex sub-problem for $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$:

For $\gamma>0$ with fixed $\check{L}\in\mathbb{R}^{2n\times m}$ and $\check{R}\in\mathbb{R}^{2n\times n}$, we define the convex sub-problem $\mathcal{C}_{\gamma}(\check{L},\check{R};\theta)$ as

[TABLE]

Note that $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ presents two competing objectives with relative importance balanced by the weight $\gamma$. On one hand, we have the truncated nuclear norm term, associated with the residual of the Lyapunov equation (10). On the other hand, we have the 1-norm penalty aiming to promote sparsity on the design variable $\Delta$. As a result, a sequential optimization strategy similar to the one applied for $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ can introduce an unwanted side-effect: depending on the magnitude of $\gamma$, convergence in terms of the truncated nuclear norm is not guaranteed. More specifically, while the overall cost of $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ can be still assured to be monotonically non-increasing (using similar arguments from Theorem 3), higher values of $\gamma$ might promote iterations where a decrease in the overall objective function (including the penalty term $\gamma\|\Delta\|_{1}$) will be obtained at the expense of an increase in the term associated with the truncated nuclear norm $\|\mathcal{Z}(W,H,\Delta)\|_{\ast}-\operatorname{tr}\{\check{L}\,\mathcal{Z}(W,H,\Delta)\,\check{R}^{{}^{\intercal}}\}$.

To control this effect, we propose an iterative procedure that seeks an approximation for the largest value of $\gamma$ for which $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ can be solved. The proposed procedure begins by solving $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$$(\gamma)$ with $\gamma=0$. In this configuration, $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$$(\gamma)$ is equivalent to the unpenalized problem $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$. Therefore, Algorithm 1 can be applied to achieve convergence as established in Theorem 3. Then, we attempt to solve $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$$(\gamma)$ for increasing values of $\gamma$, using the solution of the current problem as an initialization for the next problem, until a stopping criterion is met. This type of strategy is commonly referred to as regularization path, and has been applied to control problems, for instance, in [53, 46].

Formally, we consider a sequence $\{\gamma_{t}\}_{t=1}^{N}$ of increasing positive weights, and begin by applying Algorithm 1 to solve \operatorname{\mathcal{P}_{2-\mathrm{DN}}}$$(\gamma_{0}) with a preliminary weight $\gamma_{0}=0$. If Algorithm 1 fails to produce a feasible solution at convergence, we declare $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ infeasible. Otherwise, if it produces a solution $\mathcal{Z}(\bar{W},\bar{H},\bar{\Delta})$ with $\eta_{2n}(\mathcal{Z}(\bar{W},\bar{H},\bar{\Delta}))<\epsilon_{\eta}$, we make $Z^{(0)}\equiv\mathcal{Z}(\bar{W},\bar{H},\bar{\Delta})$ and use $L^{(0)}=[u_{1}^{(0)},\ldots,u_{2n}^{(0)}]^{{}^{\intercal}}$ and $R^{(0)}=[v_{1}^{(0)},\ldots,v_{2n}^{(0)}]^{{}^{\intercal}}$ from $\operatorname{svd}\{Z^{(0)}\}$ as initial parameters for \operatorname{\mathcal{P}_{2-\mathrm{DN}}}$$(\gamma_{1}). Then, for each $\gamma_{t}$, we seek to solve \operatorname{\mathcal{P}_{2-\mathrm{DN}}}$$(\gamma_{t}) by a sequence of convex subproblems $\{\mathcal{C}_{\gamma_{t}}(L^{(k)},R^{(k)};\theta)\}_{k}$ and evaluate the stopping condition in terms of the inner-loop solution $Z^{(k)}\equiv\mathcal{Z}(W^{(k)},H^{(k)},\Delta^{(k)})$ to each $\mathcal{C}_{\gamma_{t}}(L^{(k)},R^{(k)};\theta)$, as follows. If $\eta_{2n}(Z^{(k)})<\epsilon_{\eta}$, we consider the algorithm to have converged for the current weight $\gamma_{t}$, and move on to the next weight in the sequence. Otherwise, we choose to stop the sequence if $\eta_{2n}(Z^{(k)})\geq\eta_{2n}(Z^{(k-1)})$ holds for $K>1$ successive iterations of $\mathcal{C}_{\gamma_{t}}(L^{(k)},R^{(k)};\theta)$, where $K$ is a parameter of choice. For this purpose, we define the function $\mathrm{stop}_{K}(Z^{({\min\{0,k-K+1\}})},\ldots,Z^{(k)})$, which returns true if $\eta_{2n}(Z^{(k)})\geq\eta_{2n}(Z^{(k-1)})$ for $k-K+2,\ldots,k$ when $k\geq K$, and false otherwise. The proposed procedure is summarized in Algorithm 2.

IV Computational Experiments

To illustrate the effectiveness of our proposed approaches, in this section we perform several computational experiments considering both worst-case and average reachability designs. In the first set of experiments, we analyze random networks generated by the directed Erdős-Rényi model. The main goal is to verify the convergence of our algorithm for different random system realizations and different reachability objectives. As we will illustrate, our algorithm typically reaches solutions characterized by a very low value (i.e., below a pre-specified tolerance) of the truncated nuclear norm after a relatively small number of iterations.

In the second set of experiments, we examine a networked system with the topology of the IEEE 14-bus system[54].

We take inspiration from [6], which considers the problem of improving transient stability properties of power grids to damp frequency oscillations and prevent rotor angle instability. In this setting, the physical design variables are associated with the placement of high voltage direct current (HVDC) links, which are modeled as ideal AC current sources on the terminal buses [55]. Further, in their problem formulation, the nonlinear swing equations of system are linearized, and the HVDC placements are evaluated using controllability Gramian metrics. Our presentation consists of a simplification of the aforementioned experiment, with the goal of illustrating the effects of sparsity obtained by applying the procedure for design with structural penalties described in Section III-D. Further, as described in our problem statement, we restrict our edge design variables to follow the existing network topology. The code and data generated for both sets of experiments are available in [56].

IV-A Erdős-Rényi

We generate $L=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.p.l}{thevalue}$ random realizations of directed Erdős-Rényi (ER) systems, with state dimension $n=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.p.n}{thevalue}$ and input dimension $m=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.p.m}{thevalue}$. Each system $l=1,\ldots,L$ is defined by a pair $(A^{(l)},B^{(l)})$ that is generated as follows: The sparsity pattern encoded by the set $\{(i,j):i,j=1,\ldots,n;(i,j)\in\mathcal{G}\}$, is obtained by following the ER process until the resulting density of nonzero entries, i.e., $\|A^{(l)}\|_{0}/n^{2}$, reaches a value of $0.5$. The weights of the edges in the network are sampled from a standard uniform distribution, i.e., $[A^{(l)}]_{i,j}\sim\text{uniform}(0,1)$, for all $(i,j)\in\mathcal{G}$, with self-loops being allowed. To assure stability, the entries of each matrix $A^{(l)}$ were simultaneously scaled such that the absolute value of the largest eigenvalue of the matrix was less than one. The entries of the input matrices $B^{(l)}=[b_{1}^{(l)}|\ldots|b_{m}^{(l)}]$ were selected to have each column $b_{j}$ ($j=1\ldots,m$) defined as a canonical indicator vector $e_{\pi_{j}(n)}$, where $\pi_{j}(n)$ denotes the index of the entry equal to $1$ and is obtained as a random permutation of the $1,\ldots,n$ possible indices. Each pair was tested for reachability by assuring that $\operatorname{rank}[\,\mathcal{C}(A^{(l)},B^{(l)})\,]=n$, where $\mathcal{C}(A,B)=[\,B\,|\,AB\,|\cdots|\,A^{n-1}B\,]$.

We consider two types of design problems: (i) design for worst-case reachabililty, associated with the minimum eigenvalue $\lambda_{1}(W)$, and (ii) design for average reachability, associated with $\tau=\frac{1}{n}\operatorname{tr}\{W^{-1}\}$. For each objective, we explore two cases: one with a low target improvement value, and one with a high target improvement value. For the case of design for worst-case reachabililty, we define the ratio of improvement ${}_{\lambda}=\tilde{\lambda}_{1}/\lambda_{1}$ and fix target values $\tilde{\ratio}_{\lambda}^{\text{\,low}}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.lam_{m}in_{g}ains(1)}{thevalue}$ and $\tilde{\ratio}_{\lambda}^{\text{\,high}}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.lam_{m}in_{g}ains(2)}{thevalue}$. For the case of design for average reachability, we define the ratio of improvement ${}_{\tau}=\tilde{\tau}/\tau$ and fix target values $\tilde{\ratio}_{\tau}^{\text{\,low}}=\frac{1}{\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.tr_{i}nv_{g}ains_{r}ec(1)}{thevalue}}$ and $\tilde{\ratio}_{\tau}^{\text{\,high}}=\frac{1}{\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.tr_{i}nv_{g}ains_{r}ec(2)}{thevalue}}$. The maximum and minimum allowed perturbation magnitudes $[\Delta]_{i,j}$ were set to $\upsilon_{i,j}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.amax}{thevalue}$ and $\iota_{i,j}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.amin}{thevalue}$, respectively, for all $i$ and $j$. We then observe the evolution of the truncated nuclear norm $\eta_{2n}(Z^{(k)})$ as a function of the iteration $k$ for each system realization, until a stopping criterion is met. In particular, this criterion was set to $\epsilon_{\eta}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.m.tol_{t}nn}{thevalue}$, i.e., the algorithm stops when it reaches an iteration $k^{\star}$ for which $\eta_{2n}(Z^{(k^{\star})})\leq\epsilon_{\eta}$. The results from the execution of the algorithm are presented in Figure 1. It can be seen that $\eta_{2n}(Z^{(k)})$ reached the threshold $\epsilon_{\eta}$ for all cases considered, indicating that the desired reachability improvement, as captured by the constraint $W\in\mathcal{W}_{\theta}$, was feasible in relation to the structural constraints imposed by $\Delta\in\mathcal{D}$. Further, the median iteration value $k^{\star}$ for which such threshold was achieved is below 100 for the four scenarios considered. Finally, it can observed that the iteration for which the desired improvement in reachability is achieved typically coincides with the iteration at which the truncated nuclear norm reaches the lowest point.

IV-B IEEE Electric Power Network

We generate a network following the topology of the IEEE 14-bus system [54], with state dimension $n=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.p.n}{thevalue}$ and input dimension $m=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.p.m}{thevalue}$. The maximum and minimum allowable perturbation magnitudes $[\Delta]_{i,j}$ are set to $\upsilon_{i,j}=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.s.amax}{thevalue}$ and $\iota_{i,j}=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.s.amin}{thevalue}$, respectively, for all $i$ and $j$. As a simplification of the experiments presented in [6], the initial weights of the network were symmetrically associated with the resistance values of the transmission lines, with particular numerical values set to those available in [57]. The resulting matrix $A$ has sparsity pattern and weights as displayed next, with values rounded for compactness.

[TABLE]

In the above matrix, the symbol ‘$\cdot$’ denotes an absence of interconnection, corresponding to an entry with numerical value [math]. In particular, the network represented by $A$ has a total of $\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.card_{A}0}{thevalue}$ edges out of $\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.num_{A}0}{thevalue}$ possible, resulting in a density of $0.204$ nonzero entries.

In a similar fashion to the previous experiment, we consider two types of design: (i) design for worst-case reachabililty, associated with the minimum eigenvalue $\lambda_{1}(W)$, and (ii) design for average reachability, associated with $\tau=\frac{1}{n}\operatorname{tr}\{W^{-1}\}$. For each objective, we explore two cases: one with a low target improvement value, and one with a high target improvement value. For case of design for worst-case reachabililty, we define the ratio of improvement ${}_{\lambda}=\tilde{\lambda}_{1}/\lambda_{1}$ and set target values $\tilde{\ratio}_{\lambda}^{\text{\,low}}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.lam_{m}in_{g}ains(1)}{thevalue}$ and $\tilde{\ratio}_{\lambda}^{\text{\,high}}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.lam_{m}in_{g}ains(2)}{thevalue}$. For the case of design for average reachability, we define the ratio of improvement ${}_{\tau}=\tilde{\tau}/\tau$ and set target values $\tilde{\ratio}_{\tau}^{\text{\,low}}=\frac{1}{\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.tr_{i}nv_{g}ains_{r}ec(1)}{thevalue}}$ and $\tilde{\ratio}_{\tau}^{\text{\,high}}=\frac{1}{\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.tr_{i}nv_{g}ains_{r}ec(2)}{thevalue}}$.

To evaluate the effect of the sparsity inducing penalty, we define the cardinality index $\alpha(\Delta)\coloneqq\|\Delta\|_{0}/\|A\|_{0}$, which aims at computing the density of nonzero entries of $\Delta$ in terms of the available system entries, as induced by the sparsity pattern of the original system matrix $A$. We solve $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ using Algorithm 2 for $\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.s.num_{g}am}{thevalue}$ different values of the penalization parameter $\gamma$, whose logarithm values are set to be uniformly spaced in the pre-specified interval $\log_{10}\gamma\in[\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.s.min_{g}am}{thevalue},\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.s.max_{g}am}{thevalue}]$.

In practice, this range just needs to be chosen wide enough such that its lower limit allows $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ to be solved within the prescribed tolerance, and, conversely, its upper limit causes $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ not to be solved (i.e, the $\mathrm{stop}_{K}$ function returns true at some iteration $k^{\star}$).

In particular, Algorithm 2 is set to stop at iteration $k^{\star}$ if $\eta_{2n}(Z^{(k)})\geq\eta_{2n}(Z^{(k-1)})$ holds for $K=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.m.max_{n}o_{d}ecrease}{thevalue}$ successive iterations preceding $k^{\star}$. The results from the execution of the algorithm are presented in Figure 2. We notice the decrease of the penalty term $\|\Delta\|_{1}$ associated with a decrease in the cardinality index $\alpha(\Delta)$, for all the four cases studied.

The total number of iterations (i.e., convex subproblems solved) for the worst-case controllability metric was of $47$ and $61$ for the low and high improvement ratios, respectively. Likewise, the total number of iterations for the average controllability metric was of $49$ and $60$, respectively, for the low and high improvement ratios.

Further, for concreteness, we display the specific values of $\Delta$ for the initial and final values of the penalization weight $\gamma$, considering the scenario where we seek the design for average reachability with a high target value of improvement $\tilde{\ratio}_{\tau}^{\text{\,high}}=\DTLifdbexists{res_{e}rdos}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{e}rdos}{res_{e}rdos.txt}}\DTLfetch{res_{e}rdos}{thekey}{par.s.lam_{m}in_{g}ains(2)}{thevalue}$ (c.f. panel (h) in Figure 2). The entries of the perturbation matrix obtained for the initial value of the penalization parameter $\gamma_{\text{\,first}}=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.gammas.first}{thevalue}$ were

[TABLE]

Here, the symbol ‘$\ast$’ means that the specific entry had a value approximately zero (i.e., within a tolerance $\epsilon_{\text{s}}=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.m.tol_{s}parsity}{thevalue}$), even though the original network topology and sparsity constraints allowed a non-zero intervention value. More specifically, $\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.card_{f}irst}{thevalue}$ out of $\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.card_{A}0}{thevalue}$ non-zero possible entries were used. The algorithm was executed for increasing values of $\gamma$ until the stopping criterion was met, in particular, occurring for $\gamma_{\text{\,last}}=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.gammas.stopped}{thevalue}$. The penalized values obtained in this case were given by

[TABLE]

Here, the symbol ‘$\circledast$’ indicates that the corresponding entry resulted in an approximately zero value (i.e., within a tolerance $\epsilon_{\text{s}}=\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{par.m.tol_{s}parsity}{thevalue}$) for this value of $\gamma_{\text{\,last}}$, whereas the same entry took a nonzero value when the penalization weight $\gamma_{\text{\,first}}$ was considered. In particular, while $\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.card_{f}irst}{thevalue}$ nonzero entries were used for $\gamma_{\text{\,first}}$, this number was reduced to $\DTLifdbexists{res_{i}eeebus}{}{\DTLloaddb[noheader,keys={thekey,thevalue}]{res_{i}eeebus}{res_{i}eeebus.txt}}\DTLfetch{res_{i}eeebus}{thekey}{hst.card_{l}ast}{thevalue}$ for $\gamma_{\text{\,last}}$, as a result of effect of the structural penalty.

V Conclusion

In this paper, we have formulated and solved two problems involving the tuning of edge weights in a given discrete-time networked dynamical system such that certain reachability requirements, defined in terms of the reachability Gramian, are satisfied. In our first problem, we aimed at finding a feasible tuning of the edge weights. A direct formulation of this problems results in highly nonlinear optimization program. In order to overcome this challenge, we proposed a chain of transformations allowing us to reformulate this problem as an optimization program involving a rank constraint over a structured matrix presenting an affine dependence on the decision variables. We then relax this rank constraint using a truncated nuclear norm and proposed a sequence of convex programs to solve this relaxation. Furthermore, we have also considered a second problem in which we aimed at finding edge-weights in order to satisfy certain reachability requirements while tuning a small number of edges. Our computational approach to solve these problems has been illustrated with several numerical experiments. As future work, we plan to examine a more comprehensive class of systems, including bilinear and stochastic systems, through their corresponding reachability Gramians. Another interesting avenue of investigation would be to provide insights on the graph-theoretic characteristics of optimal designs produced for different network topologies.

Appendix A Additional Lemmas

Lemma A.1*:*

(Uniqueness for the Lyapunov equation) A solution $W\in\mathbb{S}^{n}$ to

[TABLE]

exists and is unique for any matrices $A\equiv A(\mathcal{G})\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times m}$ except for a proper algebraic variety $\mathcal{V}_{0}\subset\mathbb{R}^{|\mathcal{E}|}$, where $|\mathcal{E}|$ is the number of free entries in $A$.

Proof.

Existence and uniqueness of a solution $W\in\mathbb{S}^{n}$ to (17) can determined by examining the result of applying the vectorization operator on both sides to get

[TABLE]

where the symbol $\otimes$ denotes the Kronecker product, and the function $\mathrm{vec}(\cdot)$ is the vectorization operator. Equation (18) will have a unique solution whenever the coefficient matrix $(A\otimes A-I_{n^{2}})$ is nonsingular. Following [50], we let $a_{\mathcal{E}}\coloneqq\left([A]_{i,j}:(j,i)\in\mathcal{E}\right)$ represent an ordered set containing the entries of $A$ in lexicographic order. Next, we define a correspondence between $a_{\mathcal{E}}$ and a vector $z\in\mathbb{R}^{d},\,d=|\mathcal{E}|$, and notice that $\varphi(z)\coloneqq\det(A\otimes A-I_{n^{2}})$ is a polynomial function of the components of $z$. Then, we observe that the set $\mathcal{V}_{0}\coloneqq\{z\in\mathbb{R}^{d}:\varphi(z)=0\}$ defines a proper algebraic variety of $\mathbb{R}^{d}$ [58] where the matrix $(A\otimes A-I_{n^{2}})$ is singular. Therefore, for any matrix $A$ having entries from the correspondence between $a_{\mathcal{E}}$ and $z$ such that $z\in\mathbb{R}^{d}\setminus\mathcal{V}_{0}$, the matrix $(A\otimes A-I_{n^{2}})$ will be nonsingular, and (17) will have a unique solution $\mathrm{vec}(W)=-(A\otimes A-I_{n^{2}})^{-1}\cdot\mathrm{vec}(BB^{{{}^{\intercal}}})$. ∎

Lemma A.2*:*

(Stability from the Lyapunov equation)

Consider the discrete-time Lyapunov equation (17) with a unique solution $W$. If $W\succ 0$ and the pair $(A,B)$ is reachable, then the matrix $A$ is Schur stable.

Proof.

The proof is a trivial extension to discrete-time systems of the proof to Theorem 12.5 in [2, p.103]. To begin, we pick a left eigenvector $v$ of $A$ such that $A^{{{}^{\intercal}}}v=\lambda v$. Then, we compare the quadratic forms for $v$ at both sides of (17):

[TABLE]

where $v^{\ast}$ denotes the conjugate-transpose of $v$. Because we assumed that $W\succ 0$, it is the case that $v^{\ast}Wv>0$. Then, since $(A,B)$ is reachable by assumption, from the Popov-Belevitch-Hautus (PBH) test for controllability [2, c.f. Theorem 12.3, p.101], there is no eigenvector $v$ of $A^{{{}^{\intercal}}}$ such that $B^{{{}^{\intercal}}}v=0$. Therefore, we have that $\|B^{{{}^{\intercal}}}v\|^{2}>0$, which implies $|\lambda|<1$ in (19). Hence, the matrix $A$ is Schur stable. ∎

Lemma A.3* (Trace-inverse as semidefinite constraint):*

The condition $n\tau-\operatorname{tr}\{\left[W\right]^{-1}\}\geq 0$ for $W\in\mathbb{S}_{++}^{n}$ can be formulated as a semidefinite constraint requiring the existence of a variable $P\in\mathbb{R}^{n\times n}$ such that

[TABLE]

Proof.

Note that $P-W^{-1}\succeq 0\Rightarrow\operatorname{tr}\{P\}-\operatorname{tr}\{W^{-1}\}\geq 0$. Then, applying the Schur complement on $P-W^{-1}\succeq 0$ yields the relationship in terms of the inverse of $W$. ∎

Lemma A.4* (Von Neumann’s Trace Inequality):*

For any $X\in\mathbb{R}^{m\times n}$ and pair $(L,R)\in\{L\in\mathbb{R}^{r\times m},R\in\mathbb{R}^{r\times n}:LL^{{}^{\intercal}}=I_{r},RR^{{}^{\intercal}}=I_{r}\}$, where $1\leq r\leq\min\{m,n\}$, we have

[TABLE]

Further, consider the singular value decomposition $X=U\Sigma V^{{}^{\intercal}}$, where $U=[u_{1},\ldots,u_{m}]$ and $V=[v_{1},\ldots,v_{n}]$. Then, (20) holds with equality if $L=[u_{1},\ldots,u_{r}]^{{}^{\intercal}}$ and $R=[v_{1},\ldots,v_{r}]^{{}^{\intercal}}$.

Proof.

See Theorem 3.1 [52] and Theorem 7.4.1.1 [59, p. 458]. ∎

Appendix B Additional Conditions for Optimality of $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ and $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$

The results established in Theorem 3 guarantee that Algorithm 1 will converge to a limit value in terms of $\eta_{2n}(\mathcal{Z}(W^{(k)},H^{(k)},\Delta^{(k)}))$. However, because of the non-convexity of $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$, such a limit value does not need to correspond to its optimal value $\eta_{2n}(\mathcal{Z}(W^{(k)},H^{(k)},\Delta^{(k)}))=0$, attainable when $\operatorname{\mathcal{P}_{1}}$ is feasible. This fact motivates us to seek additional conditions for optimality of $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ by examining limit points associated with the limit values attained by Algorithm 1 in terms of their Karush-Kuhn-Tucker (KKT) conditions. Next, with this intent, we introduce a standardized version for $\operatorname{\mathcal{P}_{1-\mathrm{SUB}}}$.

$\operatorname{\mathcal{P}_{1-\mathrm{STD}}}$Standard form for $\operatorname{\mathcal{P}_{1-\mathrm{SUB}}}$:

This form consists of expressing $\operatorname{\mathcal{P}_{1-\mathrm{SUB}}}$ in terms of a single unstructured matrix variable $X\in\mathbb{R}^{4n\times 3n}$, along with affine and semidefinite constraints. Specifically, we introduce the equality constraint $X=\mathcal{Z}(W,H,\Delta)$, along with the reachability constraint $W\in\mathcal{W}_{\theta}$ and the structural constraint $\Delta\in\mathcal{D}$. Then, we jointly encode these three constraints by an equality constraint $\mathcal{A}(X)=a_{0}$ and a semidefinite constraint $\mathcal{B}(X)\succeq B_{0}$. Here, $\mathcal{A}:\mathbb{R}^{4n\times 3n}\rightarrow\mathbb{R}^{d_{\mathcal{A}}}$ and $\mathcal{B}:\mathbb{R}^{4n\times 3n}\rightarrow\mathbb{S}^{d_{\mathcal{B}}}$ are linear operators222

The operator $\mathcal{A}(X):\mathbb{R}^{4n\times 3n}\rightarrow\mathbb{R}^{d_{\mathcal{A}}}$ can be concretely expressed as $\mathcal{A}(X)=M\operatorname{vec}(X)$ for some matrix $M\in\mathbb{R}^{d_{\mathcal{A}}\times 4n\cdot 3n}$. The operator $\mathcal{B}(X)$ can be expressed as $\mathcal{B}(X)=\sum_{i=1}^{m}\sum_{j=1}^{n}Q_{i,j}[X]_{i,j}$ for symmetric matrices $\{Q_{i,j}\in\mathbb{S}^{d_{\mathcal{B}}}\}_{i,j=1}^{m,n}$.

with ${d_{\mathcal{A}}}$ and ${d_{\mathcal{B}}}$ depending on specific forms of $\mathcal{W}_{\theta}$ and $\mathcal{D}$. Further, we denote the term $\operatorname{tr}\{LXR^{{}^{\intercal}}\}$ by its inner product representation $\left\langle C,X\right\rangle$, where $C\coloneqq L^{{}^{\intercal}}R$. Therefore, the standard form representation of $\operatorname{\mathcal{P}_{1-\mathrm{SUB}}}$ is described as

[TABLE]

Lemma 3* (Optimality conditions for $\operatorname{\mathcal{P}_{1-\mathrm{STD}}}$):*

Consider a point $X^{\star}$ with rank $q$ and singular value decomposition $U\Sigma V^{{{}^{\intercal}}}=\operatorname{svd}\{X^{\star}\}$, where $\Sigma=\operatorname{diag}(\sigma_{1},\ldots,\sigma_{q},0,\ldots,0)$, $U\in\mathbb{R}^{4n\times 3n}$ with $U=\left[U_{q}|U_{y}\right]$, $U_{q}=\left[u_{1}|\ldots|u_{q}\right]$ and $U_{y}=\left[u_{q+1}|\ldots|u_{3n}\right]$, and $V\in\mathbb{R}^{3n\times 3n}$ with $V=\left[V_{q}|V_{y}\right]$, $V_{q}=\left[v_{1}|\ldots|v_{q}\right]$ and $V_{y}=\left[v_{q+1}|\ldots|v_{3n}\right]$. Also, consider the following set, associated with the subdifferential of the nuclear norm of $X$ at $X^{\star}$:

[TABLE]

Further, let $\mu\in\mathbb{R}^{d_{\mathcal{A}}}$ and $\Gamma\in\mathbb{S}^{d_{\mathcal{B}}}$, be Lagrange multipliers for the constraints associated with operators $\mathcal{A}$ and $\mathcal{B}$, respectively, and define the mapping $G:\mathbb{R}^{d_{\mathcal{A}}}\times\mathbb{S}^{d_{\mathcal{B}}}\rightarrow\mathbb{R}^{4n\times 3n}$ as

[TABLE]

Here, $\mathcal{A}^{\ast}$ and $\mathcal{B}^{\ast}$ denote the adjoint333

An adjoint operator $\mathcal{A}^{\ast}(X)$ with respect to an operator $\mathcal{A}(X)$ and inner product $\left\langle\cdot,\cdot\right\rangle$ is such that $\left\langle\mathcal{A}(X),a_{0}\right\rangle=\left\langle X,\mathcal{A}^{\ast}(a_{0})\right\rangle$.

of their respective operators. Then, for a primal-dual feasible point $X^{\star},(\mu^{\star},\Gamma^{\star})$ to be optimal for $\operatorname{\mathcal{P}_{1-\mathrm{STD}}}$ it needs to satisfy the complementary slackness, and, additionally, the Lagrangian stationary condition

[TABLE]

for some $Y\in\mathcal{Y}|_{X^{\star}}$.

Proof.

Applying the KKT conditions to the convex problem $\operatorname{\mathcal{P}_{1-\mathrm{STD}}}$, we have that Lagrangian stationarity requires

[TABLE]

where $\partial\|X^{\star}\|_{\ast}$ denotes the subdifferential of the nuclear norm at $X^{\star}$. Using the conjugacy property of linear operators and evaluating the gradient of the above equation implies that

[TABLE]

Then, using Lemma B.2 (in the Appendix) for the subdifferential of the nuclear norm, condition (24) is obtained. ∎

Next, we use the conditions specified above to analyze stationary points of $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$, further characterizing such solutions in terms of their optimality.

Theorem 4* *(Optimality conditions

for $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$):

Consider a primal feasible stationary limit point $\bar{X}$ for $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$ along with its singular value decomposition and subdifferential $\mathcal{Y}|_{\bar{X}}$ (as defined in Lemma 3). Further, consider its dual-feasible point $(\bar{\mu},\bar{\Gamma})$, for which the corresponding complementary slackness conditions hold. Then, if the Lagrange multipliers $(\bar{\mu},\bar{\Gamma})$ are such that

[TABLE]

we have that $\bar{X},(\bar{\mu},\bar{\Gamma})$ attains $\eta_{2n}(\bar{X})=0$.

Proof.

The proof is by contradiction. We assume that $q>2n$, implying $\eta_{2n}(\bar{X})>0$. A limit point $\bar{X}$ will satisfy $X^{(k+1)}=X^{(k)}=\bar{X}$ for $k\rightarrow\infty$. Considering the updates performed by Algorithm 1, we have $L^{(k+1)}=L^{(k)}=\bar{L}$ and $R^{(k+1)}=R^{(k)}=\bar{R}$ such that $\bar{L}=[\bar{u}_{1},\ldots,\bar{u}_{2n}]^{{}^{\intercal}}$, and $\bar{R}=[\bar{v}_{1},\ldots,\bar{v}_{2n}]^{{}^{\intercal}}$. From Lemma 3, the Lagrange stationarity condition for $\bar{X}$ requires

[TABLE]

We now split the term $\bar{U}\bar{V}^{{}^{\intercal}}$ as the sum $\bar{U}\bar{V}^{{}^{\intercal}}=\sum_{i=1}^{2n}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}+\sum_{i=2n+1}^{q}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}$ and compare it with the product $\bar{L}^{{}^{\intercal}}\bar{R}=\sum_{i=1}^{2n}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}$, as implied by the stationarity condition. This allows us to rewrite (26) as $G(\bar{\mu},\bar{\Gamma})=Y+\sum_{i=2n+1}^{q}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}},$ where the common terms between $\bar{U}\bar{V}^{{}^{\intercal}}$ and $\bar{L}^{{}^{\intercal}}\bar{R}$ have been canceled. Since, from Lemma 3, the optimality conditions require that $Y\in\mathcal{Y}_{\bar{X}}$, we note that any right- or left-singular vectors of $Y$ must be orthogonal to the right- and left-singular vectors appearing in $\sum_{i=2n+1}^{q}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}$. Subsequently, as $G(\bar{\mu},\bar{\Gamma})$ lies in $\mathcal{Y}_{\bar{X}}$ by (25), we must have $Y=G(\bar{\mu},\bar{\Gamma})$. This implies $\sum_{i=2n+1}^{q}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}=0$, which is impossible for $q>2n$. Therefore, it is the case that $q=2n$, since the minimum rank of $\bar{X}$ is $2n$, by Theorem 2. This fact implies that $\sigma_{i}=0$ for $i=2n+1,\ldots,3n$, and, consequently, $\eta_{2n}(\bar{X})=0$. ∎

Remark 3*:*

Conditions under which sequences of points generated by updates as performed by Algorithm 1 will produce limit points that converge to stationary points can be found in Section 3 of [60].

We can interpret Theorem 4 to get some intuition on the conditions for optimality of $\operatorname{\mathcal{P}_{1-\mathrm{STD}}}$. First, we recall that primal feasibility requires constraints (21) and (22) to be satisfied at $\bar{X}$. Next, we notice from (23) that, for optimality, the subdifferential set $\mathcal{Y}|_{\bar{X}}$ is required to be orthogonal to the row and column spaces of primal-feasible point $\bar{X}$. Thus, Theorem 4 shows us how the set $\mathcal{Y}|_{\bar{X}}$ restricts the space for dual feasibility of $\operatorname{\mathcal{P}_{1-\mathrm{STD}}}$. Specifically, it requires the existence of Lagrange multipliers $\bar{\mu}$ and $\bar{\Gamma}$ such that $G(\bar{\mu},\bar{\Gamma})$ becomes contained in the low-dimensional space defined by $\mathcal{Y}|_{\bar{X}}$. Simply stated, the higher the dimension of the space required for primal feasibility – as induced by the structural and reachability constraints – the more restricted becomes the space $\mathcal{Y}_{\bar{X}}$ available for dual feasibility (and hence, for optimality). A geometrical representation of the stationary conditions for optimality obtained in Theorem 4 is presented in Figure 3.

Next, in line with the result presented in Theorem 4 for $\operatorname{\mathcal{P}_{1-\mathrm{DN}}}$, we derive a Lagrangian stationarity condition for a limit point $\bar{X}$ now associated with a particular value of $\gamma$, to satisfy $\eta_{2n}(\bar{X})\!=\!0$ for $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$.

Theorem 5* *(Optimality conditions

for $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$):

Consider a primal feasible stationary limit point $\bar{X}=\mathcal{Z}(\bar{W},\bar{H},\bar{\Delta})$ for some $\gamma$, with $\operatorname{rank}[\bar{X}]=q$ and singular value decomposition as described in Theorem 4. Also, consider its dual-feasible point $(\bar{\mu},\bar{\Gamma})$, for which the corresponding complementary slackness conditions hold. Further, define the matrix

[TABLE]

If the Lagrange multipliers $(\bar{\mu},\bar{\Gamma})$ associated with the primal-dual feasible limit point $\bar{X},(\bar{\mu},\bar{\Gamma})$ are such that

[TABLE]

where $Y\in\mathcal{Y}|_{\bar{X}}$ and $F$ is in the set $\mathcal{F}|_{\bar{X}}\!\coloneqq\{F\in\mathbb{R}^{4n\times 3n}\!:[F]_{i,j}\!=\!0\text{ if }[\bar{X}]_{i,j}\!\neq 0,\|F\|_{\infty}\leq\gamma\},$ then, we have that $\bar{X}$ attains $\eta_{2n}(\bar{X})\!=\!0$ for $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$. Here, $\operatorname{sgn}(\Delta)$ applies the signum function over each entry of $\Delta$, evaluating to $[\operatorname{sgn}(\Delta)]_{i,j}=1$ if $[\Delta]_{i,j}>0$, $[\operatorname{sgn}(\Delta)]_{i,j}=-1$ if $[\Delta]_{i,j}<0$, and $[\operatorname{sgn}(\Delta)]_{i,j}=0$, otherwise.

Proof.

We define the linear operator $\mathcal{L}_{{}^{{}_{\Delta}}}:\mathbb{R}^{4n\times 3n}\rightarrow\mathbb{R}^{4n\times 3n}$, which extracts $\Delta$ from the upper-right block of $\bar{X}$, such that $\mathcal{L}_{{}^{{}_{\Delta}}}(X;A)=\Delta^{\!\llcorner}$. In particular, we have $\|\mathcal{L}_{{}^{{}_{\Delta}}}(X;A)\|_{1}=\|\Delta\|_{1}$. This operator allows $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$ to be expressed in the standard form, with objective function written in terms of the single variable $X\in\mathbb{R}^{4n\times 3n}$ and convex constraints imposed by the linear operators $\mathcal{A}(X)=a_{0}$ and $\mathcal{B}(0)\succeq B_{0}$, to which we associate the function $G(\mu,\Gamma)\coloneqq\mathcal{A}^{\ast}(\mu)+\mathcal{B}^{\ast}(\Gamma)$. By Lemmas B.1 and B.2 in this section of the Appendix, the Lagrangian stationarity condition for the standard form associated with $\operatorname{\mathcal{P}_{2-\mathrm{DN}}}$, at the primal-dual feasible point $\bar{X},(\bar{\mu},\bar{\Gamma})$, requires

[TABLE]

where $Y\in\mathcal{Y}|_{\bar{X}}$ and $F\in\mathcal{F}|_{\bar{X}}$. By splitting $\bar{U}\bar{V}^{{}^{\intercal}}$ as the sum $\sum_{i=1}^{r}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}+\sum_{i=r+1}^{q}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}$, we have

[TABLE]

where the common terms between $\bar{U}\bar{V}^{{}^{\intercal}}$ and $\bar{L}^{{}^{\intercal}}\bar{R}$ have been canceled. Therefore, $\eta_{2n}(\bar{X})=0$ will be attained if $\sum_{i=r+1}^{q}\bar{u}_{i}\bar{v}_{i}^{{}^{\intercal}}=0$. By letting $Y$ and $F$ such that

[TABLE]

the desired condition is achieved. ∎

B-A Additional Lemmas

Lemma B.1* (Subdifferential of matrix $1$-norm):*

Let $X\in\mathbb{R}^{m\times n}$, and denote by $\operatorname{sgn}(X)$ the matrix containing the result of the signum function applied at each entry of $X$. Then, the subdifferential of $\|X\|_{1}$ is given by $\partial\|X\|_{1}=\{\operatorname{sgn}(X)+F:F\in\mathbb{R}^{m\times n},[F]_{i,j}=0\text{ if }[X]_{i,j}\neq 0,\|F\|_{\infty}\leq 1\},$ where $i=1,\ldots,m$, $j=1,\ldots,n$, and $\|F\|_{\infty}=\max_{i,j}[F]_{i,j}$.

Proof.

See [61, p.244]. ∎

Lemma B.2* (Subdifferential of matrix nuclear norm):*

Let $X\in\mathbb{R}^{m\times n}$ with rank $q$, and singular value decomposition

$X=U_{q}S_{q}V_{q}^{{}^{\intercal}}$, with $S_{q}=\operatorname{diag}(\sigma_{1},\ldots,\sigma_{q})$, $U_{q}\in\mathbb{R}^{m\times q}$, and $V_{q}\in\mathbb{R}^{n\times q}$.

Then, the subdifferential of $\|X\|_{\ast}$ is given by $\partial\|X\|_{\ast}=\{U_{q}V_{q}^{{}^{\intercal}}+Y:Y\in\mathbb{R}^{m\times n},YV_{q}=0,U_{q}^{{}^{\intercal}}Y=0,\|Y\|\leq 1\}$, where $\|Y\|$ denotes the operator norm of $Y$.

Proof.

See [62] and [42, p.481]. ∎