Optimal Sensor and Actuator Selection using Balanced Model Reduction

TL;DR

This paper introduces an efficient method for optimal sensor and actuator placement in high-dimensional systems using balanced model reduction and greedy QR pivoting, significantly reducing computational complexity.

Contribution

It develops a scalable greedy algorithm for sensor and actuator selection that approximates optimal placements with linear runtime scaling.

Findings

Algorithm efficiently approximates optimal placements.

Runtime scales linearly with system dimension.

Demonstrated on Ginzburg-Landau system, matching costly gradient descent results.

Abstract

Optimal sensor and actuator selection is a central challenge in high-dimensional estimation and control. Nearly all subsequent control decisions are affected by these sensor/actuator locations, and optimal placement amounts to an intractable brute-force search among the combinatorial possibilities. In this work, we exploit balanced model reduction and greedy optimization to efficiently determine sensor and actuator selections that optimize observability and controllability. In particular, we determine locations that optimize scalar measures of observability and controllability via greedy matrix QR pivoting on the dominant modes of the direct and adjoint balancing transformations. Pivoting runtime scales linearly with the state dimension, making this method tractable for high-dimensional systems. The results are demonstrated on the linearized Ginzburg-Landau system, for which our algorithm approximates known optimal placements computed using costly gradient descent methods.

Equations

$$\dot{\mathbf{x}}=\mathbf{Ax}+\mathbf{Bu}$$

$$\mathbf{y}=\mathbf{Cx},$$

$${\mathbb{S}}_C$$

$${\mathbb{S}}_B$$

$$\mathbf{y}=\hat{\mathbf{C}}\mathbf{x}=[x_{{\gamma }_1}{x}_{{\gamma }_2}\dots x_{{\gamma }_r}{]}^T.$$

$${\mathbf{W}}_o={\int }_0^{\infty }{e}^{{\mathbf{A}}^{\ast }t}{\mathbf{C}}^{\ast }\mathbf{C}{e}^{\mathbf{A}t}dt,{\mathbf{W}}_c={\int }_0^{\infty }{e}^{\mathbf{A}t}\mathbf{B}{\mathbf{B}}^{\ast }{e}^{{\mathbf{A}}^{\ast }t}dt,$$

$${\mathcal{E}}_c=\{{\mathbf{W}}_c^{1/2}\mathbf{x}\mid \Vert \mathbf{x}{\Vert }_2\le 1\},$$

$${\mathcal{E}}_o=\{{\mathbf{W}}_o^{1/2}\mathbf{x}\mid \Vert \mathbf{x}{\Vert }_2\le 1\}.$$

$$\Vert G{\Vert }_2^2=\frac{1}{4{\pi }^2}{\int }_0^{\infty }\mathrm{tr}(G(j\omega {)}^{\ast }G(j\omega ))d\omega .$$

$$\Vert G{\Vert }_2^2$$

$$={\int }_0^{\infty }\mathrm{tr}({\mathbf{B}}^{\ast }{e}^{{\mathbf{A}}^{\ast }t}{\mathbf{C}}^{\ast }\mathbf{C}{e}^{\mathbf{A}t}\mathbf{B})dt=\mathrm{tr}({\mathbf{B}}^{\ast }{\mathbf{W}}_o\mathbf{B})$$

$$\log |\mathbf{C}{\mathbf{W}}_c{\mathbf{C}}^{\ast }|,\log |{\mathbf{B}}^{\ast }{\mathbf{W}}_o\mathbf{B}|,$$

$${\mathbf{B}}_{\star }=\underset{\mathbf{B}}{\mathrm{argmax}}\log |\mathbf{C}{\mathbf{W}}_c{\mathbf{C}}^{\ast }|.$$

$${{\mathbb{S}}_C}_{\star }$$

$${{\mathbb{S}}_B}_{\star }$$

$$\dot{\mathbf{a}}={\mathbf{\Phi }}^{\ast }\mathbf{A\Psi a}+{\mathbf{\Phi }}^{\ast }\mathbf{Bu}$$

$$\mathbf{y}=\mathbf{C\Psi }\mathbf{a}.$$

$${\tilde{\mathbf{W}}}_c={\mathbf{\Phi }}^{\ast }{\mathbf{W}}_c\mathbf{\Phi },{\tilde{\mathbf{W}}}_o={\mathbf{\Psi }}^{\ast }{\mathbf{W}}_o\mathbf{\Psi }.$$

$${\tilde{\mathbf{W}}}_c={\tilde{\mathbf{W}}}_o=\mathbf{\Sigma },$$

$${\tilde{\mathbf{W}}}_c{\tilde{\mathbf{W}}}_o={\mathbf{\Phi }}^{\ast }{\mathbf{W}}_c{\mathbf{W}}_o\mathbf{\Psi }={\mathbf{\Sigma }}^2⟹{\mathbf{W}}_c{\mathbf{W}}_o\mathbf{\Psi }=\mathbf{\Psi }{\mathbf{\Sigma }}^2.$$

$$\Vert G-G_r{\Vert }_{\infty }\le 2\sum _{i=r+1}^n{\sigma }_i.$$

$${\mathbf{C}}_{\star }$$

$$=\underset{{\mathbb{S}}_C}{\mathrm{argmax}}|{\mathbb{S}}_C\mathbf{C}{\mathbf{\Psi }}_r{|}^2\cdot |{\mathbf{\Sigma }}_r|$$

$$=\underset{{\mathbb{S}}_C}{\mathrm{argmax}}|{\mathbb{S}}_C\mathbf{C}{\mathbf{\Psi }}_r|.$$

$${\mathbf{B}}_{\star }=\underset{{\mathbb{S}}_B}{\mathrm{argmax}}|{\mathbf{\Phi }}^{\ast }\mathbf{B}{\mathbb{S}}_B|.$$

$$\mathbf{VP}=\mathbf{QR}.$$

$$\left[\begin{array}{cc}{\hat{\mathbf{V}}}_P\mid & \ast \end{array}\right]=[\mathbf{Q}]\left[\begin{array}{cc}\mathbf{T}\mid & \ast \end{array}\right].$$

$$|R_{ii}{|}^2\ge \sum _{j=i}^k|R_{jk}{|}^2;1\le i\le k\le p.$$

$$|{\hat{\mathbf{V}}}_P|=|\mathbf{Q}||\mathbf{T}|=\prod _{i=1}^r|T_{ii}|,$$

$${\mathbb{S}}_C\triangleq ({\mathbf{P}}_{.,j}{)}^T,\text{mbox}wherej:1\to r.$$

$$({\mathbf{\Phi }}_r^{\ast }\mathbf{B})\tilde{\mathbf{P}}=\tilde{\mathbf{Q}}\tilde{\mathbf{R}}.$$

$$\hat{\mathbf{C}}={\mathbb{S}}_C\mathbf{C},\hat{\mathbf{B}}=\mathbf{B}{\mathbb{S}}_B$$

$$\mathbf{y}=\hat{\mathbf{C}}\mathbf{x}\approx \hat{\mathbf{C}}{\mathbf{\Psi }}_r{\mathbf{a}}_r,$$

$$\hat{\mathbf{x}}={\mathbf{\Psi }}_r(\hat{\mathbf{C}}{\mathbf{\Psi }}_r{)}^{-1}\mathbf{y}={\mathbf{\Psi }}_r(\hat{\mathbf{C}}{\mathbf{\Psi }}_r{)}^{-1}\hat{\mathbf{C}}\mathbf{x}.$$

Full text

I Introduction

Optimizing the selection of sensors and actuators is one of the foremost challenges in feedback control [1]. For high-dimensional systems it is impractical to monitor or actuate every state, hence a few sensors and actuators must be carefully positioned for effective estimation and control. Determining optimal selections with respect to a desired objective is an NP-hard selection problem, and in general can only be solved by enumerating all possible configurations. This combinatorial growth in complexity is intractable; therefore, the placement of sensors and actuators are typically chosen according to heuristics and intuition. In this paper, we propose a greedy algorithm for sensor and actuator selection based on jointly maximizing observability and controllability in linear time-invariant systems. Our approach (see Fig. 1) exploits low-rank transformations that balance the observability and controllability gramians to bypass the combinatorial search, enabling favorable scaling for high-dimensional systems.

To understand the challenges of sensor and actuator placement for estimation and control, we will first consider optimal sensor placement, which has mostly been used to reconstruct static signals. The primary challenge of sensor selection is that given $n$ possible locations and a budget of $r$ sensors, there are combinatorially many, $n\choose r$, configurations to evaluate in a brute-force search. Fortunately, there are heuristics that employ greedy selection of sensors based on maximizing mutual information [2] and information theoretic criteria [3]. Another popular approach relaxes sensor selection to a weighted convex combination of possible sensors [4, 5, 6], typically solved using semidefinite programming. Both heuristic approaches optimize submodular objective functions [7], which bound the distance between heuristic and optimal placement. Some objectives, such as those based on the quality of a Kalman filter, are not submodular [8]. Alternatively, sparsity-promoting optimization can be used to determine sensors and actuators [9, 10, 11], although non-differentiability of sparsity promoting terms motivates other optimization techniques [12].

Even such heuristics cannot accommodate the large dimension of many physical models, such as in fluid dynamics. Fortunately, high-dimensional systems often evolve according to relatively few intrinsic degrees of freedom. Thus, it is possible to leverage dimensionality reduction to strategically select sensors. One approach to place point sensors  [13] computes the empirical interpolation points via EIM [14] corresponding to the proper orthogonal decomposition (POD) [15] of data, to determine important locations in state space.

For systems with actuation, it is necessary to simultaneously consider the placement of sensors and actuators, since the most observable and most controllable subspaces are often different. Sensors and actuators for optimal feedback control are generally placed along the most observable and controllable directions, respectively [16, 17, 18, 19, 7], using objective functions based on the associated observability or controllability gramians. Standard metrics for evaluating a certain sensor/actuator configuration include the $H_{2}$ norm [20, 16], a measure of the average impulse response, and the $H_{\infty}$ norm to measure the worst case performance. A chief drawback is the need to recompute the controller with each new configuration of sensors and actuators given by either the gradient minimization computation or brute-force searches. Moreover, these methods do not exploit the state-of-the-art in model reduction to optimize sensor and actuator placement.

Contribution. This work develops a scalable sensor and actuator selection algorithm based on balanced truncation [21], in which modes are hierarchically ordered by their observability and controllability. We use empirical interpolation of the low-rank balanced representation to find maximally observable and controllable states. The resulting locations correspond to near-optimal point sensor and actuator configurations. The quality of our optimized configurations are evaluated using the $H_{2}$ norm of the resulting system, which is an average measure of its output energy. The closed loop $H_{2}$ norm is more relevant than open loop metrics for control performance, given a specific $H_{2}$ cost function. Our approach, when used to optimize the open loop $H_{2}$ norm, is agnostic to the specific choice of controller weight matrices, and instead maximizes the input–output energy of the reduced order model. We also show that it is possible to apply our framework to closed loop systems, demonstrating near optimal sensor and actuator selection in comparison with more expensive iterative closed loop $H_{2}$ optimization. The runtime scales linearly with the number of state variables, after a one-time offline computation of the balancing transformation, which is less expensive than iterative alternatives. The resulting sensor and actuator configurations reproduce known optimal locations at a fraction of the cost associated with competing gradient descent methods.

II Problem Setup

Consider the following linear time-invariant system with a given state-space realization

[TABLE]

with large state dimension, i.e., $n\gg 1$. It is assumed that the system is stable, and $\mathbf{B}$ and $\mathbf{C}$ are linear actuation and measurement operators that make the system observable and controllable. Our objective is to choose a minimal subset of these sensors and actuators to obtain a system that is most jointly controllable and observable. For illustration we begin with $\mathbf{B}=\mathbf{C}=\mathbb{I}$ which correspond to pointwise sensing and actuation, but in general the subset selection can be adapted for arbitrary $\mathbf{B}$ and $\mathbf{C}$. This subset selection corresponds to multiplying inputs and outputs by the selection matrices

[TABLE]

Here $\mathbf{e}_{j}$ are the canonical basis vectors for $\mathbb{R}^{n}$ with a unit entry at the selected index $j$ and zeros elsewhere, where $\boldsymbol{\gamma}=\{\gamma_{1},\dots,\gamma_{r}\}\subset\{1,\dots,p\}$ denotes the index set of sensor locations with $\mbox{card}(\boldsymbol{\gamma})=r$. Similarly, actuator selection indices are given by $\boldsymbol{\beta}=\{\beta_{1},\dots,\beta_{r}\}$. The new measurement and actuation operators are $\hat{\mathbf{C}}=\mathbb{S}_{C}\mathbf{C}$ and $\hat{\mathbf{B}}=\mathbf{B}\mathbb{S}_{B}$ respectively. In the special case $\mathbf{B}=\mathbf{C}=\mathbb{I}$, the new operators $\hat{\mathbf{C}}=\mathbb{S}_{C}\mathbb{I}$ and $\mathbf{B}_{\star}=\mathbb{I}\mathbb{S}_{B}$ select subsets of state inputs and outputs, and the output would consist of $r$ components of $\mathbf{x}$

[TABLE]

Problem statement: *What are the best $r$-subsets of a given set of $p$ sensors and $q$ actuators, where $r\ll n$? *

To answer this question, we first quantify the degree of observability and controllability for a given set of sensors and actuators, i.e. for a given choice of $\mathbf{C}$ and $\mathbf{B}$. Optimizing over these directly involves a combinatorial search, and thus a heuristic approach is necessary for high-dimensional systems.

II-A Observability and controllability

The degrees of observability and controllability for the state-space system (1) are quantified by the observability gramian $\mathbf{W}_{o}$ and controllability gramian $\mathbf{W}_{c}$

[TABLE]

which may be visualized as controllable and observable ellipsoids (Fig. 2). These depend on the actuation and measurement operators, which consist of all states reachable from a bounded initial state

[TABLE]

and all states that may be observed

[TABLE]

Because the gramians depend on $\mathbf{B}$ and $\mathbf{C}$, they are often used to evaluate the observability/controllability of a given sensor and actuator placement. One important evaluation metric is the $H_{2}$ norm of a system. It measures the average output gain over all frequencies of the input, or the output energy. For the state-space system (1) with transfer function $G(s)=\mathbf{C}(s\mathbb{I}-\mathbf{A})^{-1}\mathbf{B}$, it is given by

[TABLE]

By the Plancherel theorem, it is also defined in the time domain by the impulse response ${y_{ij}(t)=\mathbf{C}_{i}e^{\mathbf{A}t}\mathbf{B}_{j}}$ - the output in component $i$ given an impulse in input $j$,

[TABLE]

which explicitly relate each gramian to both $\mathbf{B}$ and $\mathbf{C}$. A related alternative to the average output energy metric is given by the volumetric measure, the log determinant, denoted

[TABLE]

which are the logarithms of the geometric mean of the axes of the ellipsoid skewed by $\mathbf{B}$ or $\mathbf{C}$, by comparison the trace is the arithmetic mean. This metric is introduced by Summers et al [7] to place actuators using a greedy optimization scheme for the submodular objective function

[TABLE]

For $H_{2}$ optimal control it is desirable to minimize the average gain from stochastic disturbance $\hat{w}$ to control output $\hat{z}(s)=\hat{G}(s)\hat{w}(s)$, namely, minimizing $\|\hat{G}\|_{2}$. Several strategies seek to build the controller and choose actuators simultaneously, using expensive gradient optimization schemes. The drawback of such closed loop metrics is having to recompute the gramians - an $O(n^{3})$ operation - for every iteration that selects the next best actuator. This cubic scaling may be intractable for high-dimensional systems with large $n$.

There are cases where optimizing sensors and actuators using the closed loop $H_{2}$ norm is more relevant for control [20, 16]. By contrast, our approach reverses the strategy by instead starting from a maximally actuated and sensed optimal controller, then seeks a subset of these sensors/actuators to preserve (maximize) the geometric control measure, namely

[TABLE]

Now, the gramians no longer depend on the optimization variable and need only be computed once, and both objectives are still fundamentally linked to the $H_{2}$ norm of the system. Critically, we will extract the dominant controllable and observable subspaces from a balanced coordinate transformation of the gramians.

III Balanced Model Reduction

Many systems of interest are exceedingly high dimensional, making them difficult to characterize and limiting controller robustness due to significant computational time-delays. However, even if the ambient dimension is large, there may still be a few dominant coherent structures that characterize the system. Thus, significant effort has gone into obtaining efficient reduced-order models that capture the most relevant mechanisms for use in real-time feedback control [1].

The goal of balanced model reduction is to find a transformation $\mathbf{T}$ from state-space (leaving inputs and outputs unchanged), \begin{bmatrix}\begin{array}[]{c|c}\mathbf{A}&\mathbf{B}\\ \hline\cr\mathbf{C}&0\end{array}\end{bmatrix} to \begin{bmatrix}\begin{array}[]{c|c}\mathbf{T}\mathbf{A}\mathbf{T}^{-1}&\mathbf{T}\mathbf{B}\\ \hline\cr\mathbf{C}\mathbf{T}^{-1}&0\end{array}\end{bmatrix}, such that the transformed coordinates $\mathbf{a}=\mathbf{T}^{-1}\mathbf{x}$ are hierarchically ordered by their joint observability and controllability. This permits an $r$-dimensional representation made possible by truncating the $n-r$ least observable and controllable states.

The seminal work of Moore in 1981 [21] showed it is possible to compute this coordinate system $\mathbf{\Psi}$ where the controllability and observability gramians are equal and diagonal, denoted by the balanced model

[TABLE]

Here $\mathbf{T}^{-1}\triangleq\mathbf{\Psi}$ are direct modes and $\mathbf{T}\triangleq\mathbf{\Phi}^{*}$, the adjoint modes.

The balanced state $\mathbf{a}$ is then truncated, keeping only the first $r\ll n$ most jointly controllable and observable states in $\mathbf{a}_{r}$, so that $\mathbf{x}\approx\mathbf{\Psi}_{r}\mathbf{a}_{r}$. This results in the balanced truncation model [21] G_{r}=\begin{bmatrix}\begin{array}[]{c|c}\mathbf{\Phi}_{r}^{*}\mathbf{A}\mathbf{\Psi}_{r}&\mathbf{\Phi}_{r}^{*}\mathbf{B}\\ \hline\cr\mathbf{C}\mathbf{\Psi}_{r}&0\end{array}\end{bmatrix}. Since gramians depend on the particular choice of coordinate system, they will transform under a change of coordinates. The controllability and observability gramians for the balanced truncated system are

[TABLE]

The coordinate transformation $\mathbf{\Psi}$ that makes the controllability and observability gramians equal and diagonal,

[TABLE]

is given by the matrix of eigenvectors of the product of the gramians $\mathbf{W}_{c}\mathbf{W}_{o}$ in the original coordinates:

[TABLE]

The resulting balanced system is quantifiably close to the original system in the $H_{\infty}$ norm in terms of the Hankel singular values or diagonal entries of $\mathbf{\Sigma}$

[TABLE]

In practice, computing the gramians $\mathbf{W}_{c}$ and $\mathbf{W}_{o}$ and the eigendecomposition of the product $\mathbf{W}_{c}\mathbf{W}_{o}$ in (15) may be prohibitively expensive for high-dimensional systems. Instead, the balancing transformation may be approximated with data from impulse responses of the direct and adjoint systems, utilizing the singular value decomposition for efficient extraction of the relevant subspaces. The method of empirical gramians is quite efficient and is widely used [21, 22, 23, 24]. Moore’s approach computes the entire $n\times n$ balancing transformation, which is not suitable for exceedingly high-dimensional systems. In 2002, Willcox and Peraire [23] generalized the method to high-dimensional systems, introducing a variant based on the rank-$r$ decompositions of $\mathbf{W}_{c}$ and $\mathbf{W}_{o}$ obtained from snapshots of direct and adjoint simulations. It is then possible to compute the eigendecomposition of $\mathbf{W}_{c}\mathbf{W}_{o}$ using efficient eigenvalue solvers. This approach requires as many adjoint impulse-response simulations as the number of output equations, which may be prohibitively large for full-state measurements. In 2005, Rowley [24] addressed this issue by introducing output projection, which limits the number of adjoint simulations to the number of relevant POD modes in the data. It is particularly advantageous to use these data-driven methods or low-rank alternating direction methods [25] to approximate the gramians when there are fewer than full measurements and actuation of the state.

IV Sensor & actuator optimization via QR pivoting

We now describe an efficient matrix pivoting algorithm to optimize the log determinant over the choices of sensors and actuators. The representation of the gramians in balanced truncation coordinates plays a crucial role.

IV-A Matrix volume objective

Recall the goal of optimizing a set of $r$ sensors and actuators out of a fixed set $p$ and $q$ possible choices. The budget $r$ determines the balancing rank truncation, which necessarily must be less than both $p$ and $q$. Our sensor-actuator selection can be regarding as interpolating this rank-$r$ representation, that is, choosing locations or interpolation points that are heavily weighted in the dominant $r$ balanced modes.

Summers et al [7] show that it suffices to only consider controllable or observable subspaces for selecting sensors and actuators using the log determinant objective. Thus, we can substitute rank-$r$ balanced approximation of the gramians, $\hat{\mathbf{W}}_{c}$ and $\hat{\mathbf{W}}_{o}$, into the log determinant objective

[TABLE]

This result follows from the monotonicity of logarithms and the product property of determinants, then omitting the term that is independent of the sensors, $\det\mathbf{\Sigma}_{r}$. Likewise, in the actuator case, the objective $\operatorname*{arg\rm{}max}_{\mathbb{S}_{B}}\log|\hat{\mathbf{B}}^{T}\hat{\mathbf{W}}_{o}\hat{\mathbf{B}}|$ simplifies

[TABLE]

Consider for now the case of sensor placement. The absolute determinant is a measure of matrix volume, and $\mathbb{S}_{C}$ is a row selection matrix. The transformed objectives may be viewed as a submatrix volume maximization problem, which involves choosing the optimal $r$-row selection of $\mathbf{C}\mathbf{\Psi}_{r}$ with the largest possible determinant. Finding this optimum is an NP-hard, intractable combinatorial search over all possible $r$-row submatrices of $\mathbf{C}\mathbf{\Psi}_{r}$. However, it can be optimized greedily and efficiently via one-time matrix QR factorization requiring $\mathcal{O}(pr^{2})\mbox{ and }O(qr^{2})$ operations, as described next.

IV-B * QR pivoting algorithm*

The QR factorization with column pivoting is a greedy submatrix volume optimization scheme that we will use to construct $\mathbf{C}$ and $\mathbf{B}$, given $\mathbf{\Psi}_{r}$ and $\mathbf{\Phi}_{r}$. The pivoted QR factors any input matrix $\mathbf{V}\in\mathbb{R}^{r\times p}$ into a unitary matrix $\mathbf{Q}$, and upper-triangular matrix $\mathbf{R}$, and column permutation matrix $\mathbf{P}$ so that the permuted matrix $\mathbf{VP}$ is better conditioned than $\mathbf{V}$

[TABLE]

However, we seek a well-conditioned row permutation of $\mathbf{C}\mathbf{\Psi}_{r}$. Consider the input $\mathbf{V}=(\mathbf{C}\mathbf{\Psi}_{r})^{*}$ to the QR factorization, and the leading $r\times r$ square submatrices of the permuted input on both sides of (24),$\mathbf{\hat{V}}_{P}$ and $\mathbf{T}$

[TABLE]

Each iteration of pivoting works by applying orthogonal projections to successive columns of $\mathbf{V}$ to introduce subdiagonal zeros in $\mathbf{R}$. For our purposes, $\mathbf{P}$ plays the crucial role: at each step $\mathbf{P}$ stores the column “pivot” index of the column selected at each iteration to guarantee the following diagonally dominant structure in $\mathbf{R}$

[TABLE]

Observe that the quantity of interest, the determinant of the row-selected submatrix $\mathbf{\hat{V}}_{P}$ corresponding to the subset selection of measurements, now satisfies

[TABLE]

since $\mathbf{Q}$ is unitary and $\mathbf{T}$ is upper-triangular.

Because the determinant is the product of these diagonal entries, it can be seen that diagonal dominance guaranteed by the pivoting implicitly optimizes the desired submatrix determinant. Thus $\mathbb{S}_{C}$ is constructed from the first $r$ columns of $\mathbf{P}$ transposed

[TABLE]

Actuator selection proceeds similarly to construct a submatrix of $r$ columns of $\mathbf{B}^{*}\mathbf{\Phi}_{r}$ with maximal determinant, using one additional QR factorization

[TABLE]

The solution $\mathbb{S}_{B}$ is precisely the leading $r$ columns of $\tilde{\mathbf{P}}$, $\mathbb{S}_{B}\triangleq\tilde{\mathbf{P}}_{.,j}$, and we denote by

[TABLE]

the new measurement and actuation operators obtained in this manner.

The QR pivoting routine is a standard tool in scientific computing for matrix decomposition and linear least-squares problems. We use a block accelerated implementation of classical Businger-Golub pivoting [26] in MATLAB. Recently QR pivoting was used for interpolating nonlinear terms in EIMs [14], which would otherwise require the evaluation of high-dimensional inner products. In this setting, the interpolation point selection operator is analogous to our selection operator $\mathbb{S}_{C}$ used with pointwise measurements ($\mathbf{C}=\mathbb{I}$). The algorithm can be analyzed in terms of the error between the full state and the interpolant approximation at QR pivot interpolation points. The interpolation points can now be written

[TABLE]

where $\mathbf{\Psi}_{r}$ are the POD modes of the reduced model, and $\mathbf{a}_{r}$ are the modal coefficients. Recovering the state using the interpolant in the POD basis is accomplished with standard least-squares approximation

[TABLE]

This can be expressed as a projection $\mathbb{P}_{C}\triangleq\mathbf{\Psi}_{r}(\hat{\mathbf{C}}\mathbf{\Psi}_{r})^{-1}\hat{\mathbf{C}}$ of the true state $\mathbf{x}$ into the observable subspace.

As we shall see, the upper bound on the approximation error

[TABLE]

is given by $\|(\hat{\mathbf{C}}\mathbf{\Psi}_{r})^{-1}\|_{2}=1/|T_{rr}|$. The connection between the latter and maximizing the submatrix determinant can be made explicit in terms of the Hankel singular values of $G$.

V Analysis

The best approximation to the state in the span of the direct modes is given by $\mathbf{x}_{\star}\triangleq\mathbf{\Psi}_{r}\mathbf{\Phi}_{r}^{*}\mathbf{x}$ in the ideal measurement scenario $\mathbf{y}=\mathbf{x}$, i.e. $\hat{\mathbf{C}}=\mathbb{I}$. Here the approximation is bounded by the well-known balanced truncation error

[TABLE]

where $\sigma_{k}$ are the Hankel singular values, the diagonal entries of the balanced gramian $\boldsymbol{\Sigma}$ (14).

The analysis of empirical QR interpolation in the balanced modes begins with an established result for measurements selected using QR, which states that $\|(\hat{\mathbf{C}}\mathbf{\Psi}_{r})^{-1}\|_{2}$ at most grows as $\sqrt{p}\mathcal{O}(2^{r})$.

Lemma 1** (Drmac & Gugercin [14])**

The spectral norm of $(\mathbb{S}\mathbf{U})^{-1}$ where $\mathbb{S}$ is computed from the QR factorization (23) of the full-rank matrix $\mathbf{U}\in\mathbb{R}^{p\times r}$ is bounded above

[TABLE]

We generalize this result to the setting of arbitrary linear measurements and actuation, by analyzing the residual between the state and its interpolation in balanced coordinates. Note that the residual between the state and its projection into balanced modes $\mathbf{v}=\mathbf{x}-\mathbf{x}_{\star}$ satisfies

[TABLE]

The interpolation error from QR pivot selection satisfies

[TABLE]

Substituting (29),(30) above yields the following result.

Theorem 2

The approximation error from interpolating QR-selected observations (23) in balanced truncated modes is controlled by the discarded Hankel singular values and the norms of the given measurements and direct modes

[TABLE]

The term $\|\mathbf{C}\|_{2}=\|\hat{\mathbf{C}}\|_{2}$ results from information loss when $\mathbf{C}\neq\mathbb{I}$. An analogous result is obtained for actuator selection by considering the dual problem of estimating the adjoint state from actuation matrix $\hat{\mathbf{B}}$ - which is now the measurement operator of the adjoint system. The resulting projection operator, $\mathbb{P}_{B}\triangleq\mathbf{\Phi}_{r}(\hat{\mathbf{B}}^{*}\mathbf{\Phi}_{r})^{-1}\hat{\mathbf{B}}^{*}$, now projects on the span of the adjoint modes $\mathbf{\Phi}_{r}$. Making appropriate substitutions of $\mathbb{P}_{B}$ in the above results yields the following.

Corollary 1

The approximation error from interpolating QR-selected observations (25) of the adjoint state in balanced truncated modes is controlled by the discarded Hankel singular values and the norms of the given actuators and adjoint modes

[TABLE]

We now relate the approximation error bounds using QR pivot sensors and actuators to the log determinant objectives.

Theorem 3

*Given direct modes $\mathbf{\Psi}_{r}$, QR pivot sensors $\hat{\mathbf{C}}$ guarantee the following lower bound for the log determinant *

[TABLE]

Proof:

Noting the relationship between the singular values of a matrix and its QR factorization, we can express $|\hat{\mathbf{C}}\mathbf{\Psi}_{r}|$ in terms of the diagonal entries of its $\mathbf{R}$ factor

[TABLE]

due to nondecreasing $\sigma_{i}(\hat{\mathbf{C}}\mathbf{\Psi}_{r})$ for increasing $i$. By squaring the inequality and multiplying by $|\mathbf{\Sigma}_{r}|$ we obtain

[TABLE]

where taking logarithms yields

[TABLE]

Because $\|(\hat{\mathbf{C}}\mathbf{\Psi}_{r})^{-1}\|_{2}=1/|T_{rr}|$, the upper bound (30) in Lemma 2 is the inverse lower bound for $|T_{rr}|$, which can now be substituted above to obtain the final result. ∎

An analogous lower bound can be obtained for the objective using QR pivot actuators by appropriately substituting $\hat{\mathbf{B}},\tilde{\mathbf{R}}$ and adjoint modes $\mathbf{\Phi}_{r}$ in the above proof.

Corollary 2

*Given adjoint modes $\mathbf{\Phi}_{r}$, $\hat{\mathbf{B}}$ satisfies the following lower bound for the log determinant *

[TABLE]

VI Results

We evaluate the selection algorithm in two settings. The first compares QR pivot selections with all possible sensor subset selections in a random state-space model of tractable size. Next we consider an application to closed-loop flow control using LQG control to stabilize unstable Ginzburg-Landau dynamics. The LQG controller with full actuation and sensing is also tractable, and we approximate the $H_{2}$ optimal placements computed using gradient descent [16] with our QR scheme.

VI-A Discrete random state space

Our first example investigates sensor and actuator selection for random state-space systems with randomized $\mathbf{A},\mathbf{B},\mathbf{C}$. First, we compare the results of QR sensor placement against a brute-force search across all possible sensor selections for a system with $n=25$ states and $r=7$ randomized measurements. The log determinant objective (11) is evaluated for all possible choices of 7 sensors, since the system is small enough to explicitly compute the full gramian for all ${n\choose r}=480,700$ choices of $\hat{\mathbf{C}}$. These results are binned in Fig. 3, and compared with the value resulting from our method (red line). The input to the QR scheme, the balancing modes, are computed only once from the full system. The sensors resulting from our method are observed to be near optimal for the log determinant, exceeding $99.99\%$ of all others, and also good substitutes for $H_{2}$ optimal sensors. On average, our method surpasses $99.8\%$ of possible outcomes with a standard deviation of $0.85\%$, over a randomly generated ensemble of 500 model realizations. Therefore, QR sensors are closer to optimal than the analysis suggests.

We now investigate performance on a larger random state-space model with $n=100$ states, and likewise initialize the model with randomized actuation and sensing such that $p=q=100$. Figure 4 shows the log determinant objective that is being optimized for various sensor and actuator configurations. The log determinant of the gramian volume is plotted for the truncated model with QR-optimized sensor and actuator configurations (red circles) and with random configurations (blue violin plots). The truncation level $r$ for the balanced truncation is chosen to match the sensor and actuator budget on the $x$-axis. The QR-optimized configurations dramatically outperform random configurations. As more modes are retained, the chosen sensors and actuators better characterize the input–output dynamics, and their performance gap over random placement increases over all random ensembles, giving empirical validation of our approach.

Because the system is randomly generated and the dynamics do not evolve according to broad, non-localized features in state-space, many sensors and actuators are required to characterize the system. In particular, this is reflected in the slow decay of Hankel singular values. By contrast, the next example is generated by a physical fluid flow model, and has coherent structure that allow for a more physical interpretation of sensor and actuator placements with enhanced sparsity.

VI-B Linearized Ginzburg-Landau with stochastic disturbances

We consider the closed-loop linearized Ginzburg-Landau model evolving velocity perturbations in a flow, given a controller with full actuation and sensing, which is often not feasible in practice. The equations modeling the plant dynamics are unstable because the system matrix has eigenvalues in the right half plane. The dynamical system matrix $\mathbf{A}$ is formed from Hermite pseudospectral discretization of the linearized Ginzburg-Landau operator

[TABLE]

The spatial grid $\boldsymbol{\xi}\in\mathbb{R}^{n}$ is discretized at the $n=100$ roots of Hermite polynomials, and $\nu,\beta,\mbox{ and }\mu(\xi)$ are advection, diffusion and wave amplification parameters. Each $i$th sensor $\xi_{s}$ (row of $\mathbf{C}_{2}$) and actuator at $\xi_{a}$ (column of $\mathbf{B}_{2}$) are weighted by Gaussian kernels and the trapezoidal integration weights $\mathbf{M}$

[TABLE]

The linear quadratic Gaussian (LQG) controller stabilizes the dynamics by minimizing the $H_{2}$ optimal cost function $J(\mathbf{x},\mathbf{u})=\mathbf{x}^{T}\hat{\mathbf{Q}}\mathbf{x}+\mathbf{u}^{T}\hat{\mathbf{R}}\mathbf{u}$, where $\hat{\mathbf{Q}}$ and $\hat{\mathbf{R}}$ are user-specified weight matrices. The output $\mathbf{u}$ of the LQG controller, given by

[TABLE]

stabilizes the dynamics given white noise stochastic disturbance $\mathbf{d}$ and noise $\mathbf{n}$ at all sensors and actuators, with covariances $\mathbf{V}=4\cdot 10^{-8}\mathbb{I}$ and $\mathbf{W}=\mathbb{I}$. Since every state is observed and actuated, each spatial gridpoint corresponds to one $\xi_{a}$ and $\xi_{s}$. The idea is to preserve as much of this “ideal” controller as possible using a subset of the original sensors and actuators

[TABLE]

which is similar to our original problem formulation. Hence we can perform balanced model reduction on either impulse responses or the controller directly, and then QR pivoting to optimize placements. In this formulation, $\mathbf{u}$ is the output and $\mathbf{y}$ is the input, which encapsulates the notion of maximizing the gain from feedback to $\mathbf{u}$ to stabilize the dynamics, which translates to the observabillity of (36). Thus we compute gramians and adjoint, direct modes of the LQG matrices $\mathbf{A}\triangleq\mathbf{A}-\mathbf{B}_{2}\mathbf{F}-\mathbf{LC}_{2},\mathbf{B}\triangleq\mathbf{L},\mathbf{C}\triangleq-\mathbf{F}$.

We compare our approach to established gradient descent techniques for computing the $H_{2}$ optimal controller and sensor-actuator placements simultaneously. The particular algorithm for comparison is the optimal placement for this model determined using the gradient descent scheme of Chen and Rowley [16].

Their $H_{2}$ norm optimization scheme permits placement of sensors and actuators at locations that may not be grid points. The major drawback is that each Newton iteration requires solving $2r$ $n\times n$ Lyapunov equations until convergence, although recent work simplifies this to 2 equations per iteration [27]. Furthermore, the procedure requires an ensemble of random initial conditions to avoid converging to a local minimum. In [16], the optimal placement is computed using conjugate gradient optimization for the same spatial discretization $n=100$, which becomes computationally expensive as the grid resolution increases. In this case, gradient descent is more costly than balancing the fully actuated and observed system, which comes at a one-time cost of solving 2 Lyapunov equations for the gramians, and 2 Riccati equations for the LQG gain matrices ($O(n^{3})$ each). Therefore, our algorithm is sensible when the grid discretization is sufficiently fine. Furthermore, our solution is a good starting point for the convergence of the gradient descent scheme, thus eliminating the need for optimization over a large ensemble of randomized starting points. QR pivoting runtime scales as $\mathcal{O}(nr^{2})$ and the deviation of the resulting placement from the $H_{2}$ optimum (fig. 5) decreases with increasing $r$.

Figure 5 plots sensor and actuator configurations from the QR algorithm and $H_{2}$ gradient optimization, which are compared with the $H_{2}$ optimal placements in [16]. The resulting placements for the cases $r=1$ to $r=5$ sensors and actuators are plotted vertically, and the horizontal axis is the spatial domain $\xi\in[-12,12]$ with a shaded wave amplification region in which fluid perturbations are amplified. For each value of $r$, we apply QR pivoting to the rank $r$ truncated balanced modes. QR pivoting collocates sensors and actuators, indicating that $\mathbf{A}$ is approximately symmetric and hence the direct and adjoint modes (pictured in Fig. 1) are identical up to a scaling factor. In practice, sensors are often slightly downstream to account for time delays, so we enforce via the pivoting procedure that sensors are not placed at previously chosen actuators. The $H_{2}$ norms of the resulting placement on the $y$-axis indicate that the QR selections closely approximate the optimal placements. The $H_{2}$ optimal placement [16, Fig. 4] of five sensors and actuators, with $H_{2}$ norm 27.4, agrees exactly with the $H_{2}$ optimum and is closely approximated by the QR pivoted placement (27.8).

Figure 6 compares controller performance between QR pivoting and the $H_{2}$ optimum via the LQG gain of a given signal from each sensor to each actuator. The LQG gains are identical to those produced by the $H_{2}$ optimal method of Chen and Rowley [16, Fig. 5]. The diffusive nature of the dynamics favors nearly collocating the sensors and actuators, since the high-frequency oscillations mostly propagate to the nearest actuator. This confirms that our framework is useful for optimizing sensors and actuators. Balanced truncation applied to the closed loop system is critical to achieving this, since the open loop dynamics are unstable and it is shown in [16] that the dominant eigenmodes of the dynamics lead to vastly suboptimal placements.

VII Discussion and Outlook

In this work we develop scalable sensor and actuator selection whose runtime scales linearly with the number of state variables, after a one-time offline computation of the balanced modes. Our approach relies on balanced model reduction [21, 23, 24], which hierarchically orders modes by their observability and controllability. We extend EIMs to interpolate the low-rank balancing modes of the system and determine maximally observable and controllable locations (sensor & actuators) in state space. The performance of this algorithm is demonstrated on random state-space systems, and optimal $H_{2}$ control of the linearized Ginzburg-Landau model. Our optimized placements vastly exceed the performance of random placements, and closely approximate $H_{2}$ optimal placements computed by costly gradient minimization schemes, but achieved at a fraction of the runtime.

Sensors and actuators are critical for feedback control of large high-dimensional complex systems. This work advocates sensor and actuator selection using QR pivots of the direct and adjoint modes of a system’s balancing transformation. The resulting placement is empirically shown to preserve the dynamics of the full system. The method has deep connections to system observability, controllability, modal sampling methods and classical experimental design criteria. Furthermore, QR pivoting is more computationally efficient than leading greedy and convex optimization methods, and thus critically enlarges the search space of possible selections. This is particularly valuable in spatiotemporal models where high-resolution grids generate a large number of states, and balanced modes and QR method exploit the spatial structures.

This work opens a variety of future directions in pivoting sensor and actuator optimization. Rapid advances in data collection yield extremely large search spaces, for which the computation of balanced modes and QR pivoting may be accelerated using randomized linear algebra. Our method relies on a known model of the dynamics, but it would also be interesting to generalize the method to data-driven system identification models. In addition, point sensors and actuators are simplifications of constrained or nonlinear sensing and actuation that may occur in practice. Nonlinear sensing constraints remain an open challenge.