Distributed Optimal Control of Graph Symmetric Systems via Graph Filters

As noted in §4 of [anderson_system_2019], we can compactly write the system dynamics (2) in the frequency domain as

where $\mathbf{x}=\sum_{t=0}^{\infty}z^{-t}x(t)$ is the signal $x(t)$ in the $z$-domain, and idem for $\mathbf{u}$ and $\mathbf{w}$. For a (dynamic) linear state-feedback controller $\mathbf{u}=\mathbf{K}\mathbf{x}$, it follows immediately that

where $\mathbf{\Phi}_{x}(z)\in\mathbb{C}^{N\times N}$ and $\mathbf{\Phi}_{u}(z)\in\mathbb{C}^{N\times N}$ are system responses that map the disturbance $\mathbf{w}$ to state $\mathbf{x}$ and control input $\mathbf{u}$, respectively. The following SLS theorem states that all achievable responses lie in an affine subspace of strictly proper stable rational transfer functions $\frac{1}{z}\mathcal{R}\mathcal{H}_{\infty}$.

For disturbance $w(t)\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,I)$, one can recast the optimization problem (4) in terms of system responses $\mathbf{\Phi}_{x}$ and $\mathbf{\Phi}_{u}$ as

With a slight abuse of notation, we define $J(\mathbf{\Phi}_{x},\mathbf{\Phi}_{u})$ to be the LQR cost achieved by $\mathbf{\Phi}_{x}$ and $\mathbf{\Phi}_{u}$ in the objective of (11) and $J(\mathbf{K})$ as the cost (3) achieved by applying controller $\mathbf{K}$.

We now proceed to show that under Assumption 1, the optimal system response for a graph symmetric system that solves the LQR problem (11) can be written as a graph filter.

The main implication of Theorem 2 is that the optimal linear state-feedback controller for graph symmetric systems under Assumption 1 is a graph filter and can thus be implemented via distributed message passing. We note, however, that the above result implies that the resulting optimal system response $\mathbf{\Phi}^{\star}:=(\mathbf{\Phi}^{\star}_{x},\mathbf{\Phi}^{\star}_{u})$ could be dense, as $S^{N-1}$ is dense if $S$ defines a connected graph. This can be undesirable in practice, as it requires $N-1$ communication exchanges with one-hop neighbors, potentially causing significant delays if the size $N$ of the graph is large. In the next section, we leverage a robust variant of the SLS parameterization given in Theorem 1 to restrict the optimal system responses to only the first $F\ll N-1$ filter taps while guaranteeing stability and near optimal performance.

We end this section by noting that in the graph signal processing literature \citepSandryhaila2013-DSPG, a controller of the form of (12) is known as a linear, shift-invariant (LSI) graph filter, and is analogous to an LTI filter. Note that $S\mathbf{\Phi}^{\star}\mathbf{x}=\mathbf{\Phi}^{\star}S\mathbf{x}$, hence the name. In particular, Equation (12) is a spatially finite impulse response (FIR) graph filter \citepSegarra2017-GraphFilterDesign that is completely characterized by a finite set of $N$ filter taps that can be conveniently described by a collection of transfer functions $\boldsymbol{\phi}^{\star}=[\boldsymbol{\phi}_{0}^{\star},\ldots,\boldsymbol{\phi}_{N-1}^{\star}]^{\mathsf{T}}\in\mathbb{R}^{N}$. We emphasize that the transfer functions themselves are not restricted to be temporally FIR. Spatially FIR graph filters are also known as convolutional graph filters \citepRuiz2021-GNN due to their sum-and-shift nature, understanding that the effect of the operation $S\mathbf{x}$ is to shift the signal around the graph (thus, oftentimes, the GMD $S$ is referred to as the graph shift operator). Furthermore, spatially FIR graph filters satisfy the convolution theorem that indicates that a convolution in the vertex domain can be computed by means of an elementwise multiplication in the spectrum domain \citepSandryhaila2013-DSPG. Finally, in the context of finite graphs, it is observed that the space of FIR graph filters of the form (12), characterized by $N$ filter taps, is equivalent to the space of spatially infinite impulse response (IIR) graph filters as well as autoregressive, moving average (ARMA) graph filters \citepIsufi2017-ARMA.

We propose an approach inspired by the graph signal processing literature, wherein we exploit the graph filter structure of the optimal system responses \citepSegarra2017-GraphFilterDesign. More specifically, we project the optimal system responses $\mathbf{\Phi}^{\star}$ onto graph filters of order $F$ in the $\mathcal{H}_{2}$ norm by solving the following optimization problem

Here $\boldsymbol{\phi}:=[\boldsymbol{\phi}_{0}(z),\ldots,\boldsymbol{\phi}_{F-1}(z)]^{\top}\in\mathbb{C}^{F}$ collects the $F$ transfer function filter taps defining $\mathbf{\Phi}:=\sum_{k=0}^{F-1}\boldsymbol{\phi}_{k}(z)S^{k}$. If we further restrict each transfer function filter tap $\phi_{k}(z)$ to be FIR of order $n$, i.e., if we write $\boldsymbol{\phi}_{k}(z)=\sum_{i=1}^{n}z^{-i}\phi_{k}[i]$, then this reduces to solving the following unconstrained quadratic program

Prop. 3 determines in closed-form how to compute the filter of order $F$ that best approximates the optimal linear distributed controller in the $\mathcal{H}_{2}$ norm. It also shows that each filter tap transfer function $\boldsymbol{\phi}_{k}(z)$ is obtained as the optimal tap $\boldsymbol{\phi}_{k}^{\star}(z)$ with an added corrective term that accounts for the $N-F$ taps that could not be included.

This approach is easy to implement computationally, as it only requires solving a least squares problem to minimize the projection cost. However, the resulting controller cannot be guaranteed to be stabilizing. As we show later via numerical simulation, approximations with a small number of filter taps $F$ are often unstable. This motivates an approach that takes into account the stability of the resulting controller.

Robust SLS \citepmatni2017scalable, anderson_system_2019 offers a systematic way to reason about approximate system responses, i.e., system responses that do not exactly satisfy the achievability constraint (10). In particular, as shown in the following result, robust SLS allows for an explicit characterization of the effects of using approximate system responses for controller design.

We leverage this result to provide an upper-bound on the amount of truncation that can be applied to $\mathbf{\Phi}^{\star}$ without destabilizing the system.

Corollary 5 offers a way to synthesize stable truncated system responses. Specifically, to synthesize a system response that uses only $F<N-1$ filter taps while guaranteeing both stability and performance, we propose the following optimization problem:

By Corollary 5, the solution to (18) defines a controller that is stabilizing. We further show in the next result that it enjoys guaranteed suboptimality bounds relative to the optimal controller defined by $\mathbf{\Phi}^{\star}$. We first introduce the following notation: for a system response of the form $\mathbf{\Phi}=\sum_{k=0}^{N-1}\boldsymbol{\phi}_{k}(z)S^{k}$, we define the $F$-truncation $P_{F}(\mathbf{\Phi})$ and the $F$-tail $P_{F\perp}(\mathbf{\Phi})$ as

This optimization problem is jointly quasi-convex and can be solved efficiently using bisection. Further, feasibility provides a stability certificate in the form of $\left\lVert\mathbf{\Delta}\right\rVert_{\mathcal{H}_{\infty}}<1$.

Lastly, we can combine the robustness constraints used in robust SLS with the signal-processing-based projection method. Specifically, we solve the following optimization problem

We note that solving this problem does not give an upper bound on the cost of the resulting controller, but the robustness constraint ensures that the resulting controller is stabilizing.

We end this section by detailing practical implementation details for optimization problems (18) and (20). First, we note that computationally, one cannot directly optimize for the IIR system responses as is written in (18), (20). In practice, we use an FIR approximation of the strictly proper transfer functions $\boldsymbol{\phi}_{x}(z)$ and $\boldsymbol{\phi}_{u}(z)$, i.e., $\boldsymbol{\phi}(z)$ is parameterized as

for some given FIR order $n$. As shown in [anderson_system_2019], the suboptimality incurred by such an FIR approximation decays exponentially in the horizon $n$.

The $\mathcal{H}_{\infty}$-norm constraints on the (also FIR) transfer matrix $\mathbf{\Delta}$ can then be enforced via semidefinite programming (see Theorem 5.8 in [dumitrescu2007positive]), potentially introducing a nontrivial computational burden. However, we note that one can replace the $\mathcal{H}_{\infty}$-norm constraint in optimization problems (18) and (20) with any induced norm constraint. A particularly appealing option is the $\ell_{1}\to\ell_{1}$ induced norm (which defines the $\mathcal{L}_{1}$-norm of the transpose system), as this norm decomposes columnwise. As shown in [anderson_system_2019, wang_separable_2018], the resulting robustness constraints are linear and embarrassingly parallelizable. We defer this extension to future work.

In the following experiments, we consider the distributed linear quadratic regulator (LQR) problem (4) over $N=10$ scalar subsystems. We generate the GMD $S$ and dynamic matrices $A$ and $B$ using a process similar to that in [Gama2022-DistributedLQR]. To generate a problem instance, we start by creating the communication network $\mathcal{G}$ by randomly sampling $N$ numbers $\{u_{i}\}_{i=1}^{N}\sim{}U[0,1]$, and creating a bi-directional link between $v_{i}$ and each of its $3$ nearest points as defined by the topology on the interval $[0,1]$ under the metric $d(v_{i},v_{j})=|u_{i}-u_{j}|$. We then take $S$ to be a symmetric matrix that shares the sparsity pattern of the Laplacian of $\mathcal{G}$, with its entry values sampled independently from $\mathcal{N}(0,1)$. The GMD $S$ is then normalized to have a spectral radius of $1$. We generate the dynamics matrices $A$ and $B$ to share the same eigenvectors as $S$, and sample their eigenvalues i.i.d. from the standard normal distribution – hence both $A$ and $B$ are symmetric matrices. We take the cost matrices $Q=R=I_{N}$. For both of the following experiments, we randomly generate $50$ problem instances using this process. We end by noting that $\mathcal{G}$ generated this way have, on average, a diameter of $5.92$ hops.

For the implementation of the optimization problems, we approximate the transfer functions $\boldsymbol{\phi}(z)$ with an FIR horizon of $n=10$. Further, for the robust SLS problem (18), instead of using bisection to determine the best value of $\gamma$, we fix $\gamma=0.98$, as empirically the value of $\gamma$ does not significantly affect the cost achieved by the controllers.

In this experiment, we demonstrate the importance of the robust SLS-based stability constraints in synthesizing stable distributed controllers. We vary the number of allowed filter taps and apply both the naive projection (13) and robust SLS (20) methods to $50$ randomly generated problem instances. For the naive projection method, we report the percentage of resulting controllers that are stable. For robust SLS, we report the percentage of optimization problems (18) being feasible, as feasible solutions optimization problem (18) are guaranteed to be stabilizing. The results are shown in Figure 1.

First, we observe that as expected, a higher number of filter taps result in a higher probability of synthesizing stable responses for both methods. However, the naive projection method has nonzero probability of resulting in an unstable controllers even with a large number of filter taps i.e., even when the projection error between the stable optimal system responses and the projected responses is small. On the other hand, the percentage of stable solutions resulting from the robust SLS problem increases monotonically with the number of filter taps, which is expected for a principled way of synthesizing stable controllers. Robust SLS also generally achieves a higher percentage of certifiably stable responses than that of naive projection, except for in the extremely sparse case of $F=3$. This suggests that the robust constraint might be too restrictive for synthesizing extremely sparse responses. Combined with the low computation cost of naive projection, this suggests a potential benefit of applying both methods in the sparse regime.

In this experiment, we compare the performance of robust SLS and robust projection for different filter tap numbers $F$ on the same $50$ randomly generated problem instances. We show the median (solid lines), $25$-th and $75$-th percentile (shaded regions) of the costs achieved by both methods in Figure 2.

First, we note that the median costs decreases monotonically for both methods as the number of hops increase. This shows that that the optimization problems can leverage the increase in expressivity of the graph filters to achieve better performance, which matches our intuition. Second, we note that the robust SLS-based method achieves a lower cost than robust Projection over for all numbers of filter taps considered. We also note that to the left of $4$ hops, the upper boundary of the shaded region, which represents the $75$-th percentile of the cost, is infinite, indicating that at least $25\%$ of the robust synthesis problems are infeasible. This again suggests a need to develop more flexible methods in the sparse regime.

We show that for any optimal system response $\mathbf{\Phi}^{\star}_{x}$, $\mathbf{\Phi}^{\star}_{u}$ of the optimization problem (11) that is not diagonalizable with $V$, i.e.,

are not diagonal, we can construct a simultaneously diagonalizable system response $\mathbf{\Phi}_{x}^{\prime}$, $\mathbf{\Phi}_{u}^{\prime}$ that is equally optimal. In particular, we construct such a system response as follows:

We first show that $\mathbf{\Phi}_{x}^{\prime}$, $\mathbf{\Phi}_{u}^{\prime}$ are feasible solutions of (11). From the achievability condition on $\mathbf{\Phi}^{\star}_{x},\mathbf{\Phi}^{\star}_{u}$, we have that

Using the simultaneous diagonalizability of $A$ and $B$, we have

Since the matrices $(zI-\Lambda_{A})$, $\Lambda_{B}$ and $I$ are diagonal, we have

where $\mathsf{diag}[\mathbf{M}]$ projects a matrix $\mathbf{M}$ onto its diagonal elements. Therefore, $\mathbf{\Phi}_{x}^{\prime}$, $\mathbf{\Phi}_{u}^{\prime}$ is also feasible.

Now, we show that $\mathbf{\Phi}_{x}^{\prime}$, $\mathbf{\Phi}_{u}^{\prime}$ gives a cost at least as good as that of $\mathbf{\Phi}^{\star}_{x},\mathbf{\Phi}^{\star}_{u}$. By the simultaneous diagonalizability of the matrices $Q$ and $R$, and the fact that the $\mathcal{H}_{2}$-norm is invariant under unitary transformations, we have that

Denoting the $i$-th eigenvalue of $\mathbf{Q}$ and $\mathbf{R}$ with $\lambda_{Q,i},\lambda_{R,i}$, respectively, we have the inequality

which follows from the definition of $\mathbf{\Lambda}_{x}^{\prime}$ and $\mathbf{\Lambda}_{u}^{\prime}$ in equation (22). Reversing the steps in (26), we see that $\mathbf{\Phi}_{x}^{\prime},\mathbf{\Phi}_{u}^{\prime}$ achieves a cost at least as good as that of $\mathbf{\Phi}^{\star}_{x},\mathbf{\Phi}^{\star}_{u}$. We can thus conclude that there always exists an optimal simultaneously diagonalizable system response to the LQR problem (11). The controller $\mathbf{K}^{\prime}=(\mathbf{\Phi}_{u}^{\prime})(\mathbf{\Phi}_{x}^{\prime})^{-1}$ is thus optimal and simultaneously diagonalizable by $V$.

Here we show that optimal control problems over graph symmetric systems and controllers satisfy quadratic invariance \citeprotkowitz2005characterization. Before proceeding, we remark that the analysis of LQR optimal control problem over GSSs does not require quadratic invariance. In particular, in Theorem 2 we analyze the unconstrained optimal control problem and show that the resulting unconstrained optimal controller satisfies a corresponding notion of graph symmetry. However, in the interest of completeness, we show that if such a constraint were imposed on the controller during synthesis, the resulting problem satisfies quadratic invariance.

To that end, the corresponding constrained controller synthesis problem can be stated as

where $\mathcal{S}:=\{S\in\mathcal{RH}_{\infty}\ |\ S\ \text{is diagonalizable by }V\}.$ Denoting the plant input-output transfer function as

we have the following proposition.