Progress Report n°1
This is the first progress report for the ANR JCJC project CONMAN, covering the period from January 1st, 2023 to June 30th, 2023.
1 Identification
| Project acronym | CONMAN |
| Project title | Data-driven Control and Manifold Interpolation |
| Funding program | Programme Jeunes Chercheuses Jeunes Chercheurs |
| Agreement number | ANR-22-CE46-0008 |
| Principal Investigator | Jean-Christophe Loiseau |
| Starting date | 01/01/2023 |
| Duration | 48 months |
| Author of report | Jean-Christophe Loiseau |
| Date | June 2023 |
| Period covered | 01/01/2023 - 30/06/2023 |
2 Progress Report
Applications of closed-loop flow control have epic proportions: drag reduction, lift increase, mixing enhancement, or noise mitigation. It is however challenged by strong nonlinearities, partial state information, parametric dependencies, or time delays. By focusing on linear time-invariant (LTI) dynamical systems of the form
\[ \begin{aligned} \dfrac{d\mathbf{x}}{dt} & = \mathbf{Ax} + \mathbf{Bu} \\ \mathbf{y} & = \mathbf{Cx} + \mathbf{Du}, \end{aligned} \]
the present project aims at tackling some of these problems, most notably the parametric dependencies for reduced-order models. Likewise, optimal actuator and sensor placement (i.e. design of \(\mathbf{B}\) and \(\mathbf{C}\), respectively) still is an unsolved combinatorial optimization problem whose solution impacts all downstream tasks and applications.
To tackle these issues, the CONMAN project follows two parallel, albeit complementary routes :
- The development of a unified framework for data-driven linear modeling of parameterized systems based on the joint use of classical low-rank approximation techniques for dimensionality reduction purposes combined with high-order interpolation on matrix manifolds to capture the parametric dependencies.
- The formulation of a combinatorial optimization problem for actuators and sensors placement and its subsequent relaxation into a tractable form. Particular attention is being paid to highly efficient yet relevent (from an experimental point of view) algorithms with certains optimality guarantees.
In agreement with the recommendations put forth by Prof. Lorena Barba (George Washington University, USA) during her recent invited lecture Anti Patterns of Scientific Machine Learning to Fool the Masses at the Platform for Advanced Scientific Computing Conference (PASC23), extra care is taken in choosing the baseline models against which the rest of the techniques developped during the project will be compared in order to ensure fair and meaningful benchmarks. As a consequence, the first six months of this project have been dedicated to a thorough literature review as well as implementing some of the building blocks needed for the baseline reference models.
At the time of acceptance by ANR, plan was to hire a postdoctoral researcher in the first quarter of 2023. Unfortunately, no candidate was found before May 2023. The successful hiree, Simon Kern, has defended his Ph. D. work at KTH early June 2023 and will join DynFluid at the begining of the Fall semester. These human ressources issues forced us to postpone some tasks in Work Package #1 by a few months. This should not impact however the overall calendar of the project.
This is the first technical report (D0.1) of the CONMAN project. It is organized as follows: Section 2.1 provides a thorough presentation of the baseline techniques, namely Balanced Proper Orthogonal Decomposition for reduced-order modeling and the POD-QR method for near-optimal actuators/sensors selection. Section 2.2 describes the theoretical tools already available to the project members or under active development, along with their numerical implementations. Finally, Section 2.3 discusses some of the problems faced by members of the project as well as presenting some of the actions taken for the scientific communication and dissemination while Section 2.4 teases some of the scientific advances to be expected for the next technical report (D0.2, due Dec. 2023).
2.1 Baseline models
2.1.1 Reduced-order modeling
Direct design of controllers and estimators for high-dimensional dynamical systems resulting from the discretization of partial differential equations such as the Navier-Stokes equations is often faced with the curse of dimensionality. Model order reduction is thus a necessary step before control synthesis. Numerous techniques have been proposed in the literature. In this project, we restrict our attention to Balanced Truncation (see Section 2.1.1.1) and its computationally efficient (albeit equivalent) alternative Balanced Proper Orthogonal Decomposition (see Section 2.1.1.2).
2.1.1.1 Balanced Truncation
Given a multiple input multiple output (MIMO) linear time-invariant (LTI) dynamical system
\[ \begin{aligned} \dfrac{d\mathbf{x}}{dt} & = \mathbf{Ax} + \mathbf{Bu} \\ \mathbf{y} & = \mathbf{Cx} + \mathbf{Du} \end{aligned} \]
where \(\mathbf{x} \in \mathbb{R}^n\) is the state vector of the system, Moore (1981) proposed the Balanced Truncation (BT) procedure. Starting from the observability gramian \(\mathbf{W}_o\) and controllability gramian \(\mathbf{W}_c\) solution to the Lyapunov equations
\[ \mathbf{A}^* \mathbf{W}_o + \mathbf{W}_o \mathbf{A} = - \mathbf{C}^* \mathbf{C} \quad \text{and} \quad \mathbf{A} \mathbf{W}_c + \mathbf{W}_c \mathbf{A}^* = - \mathbf{BB}^*, \]
Balanced Truncation first applies a balancing transformation, i.e. finding an invertible linear transformation \(\hat{\mathbf{x}} = \mathbf{T} \mathbf{x}\) such that
\[ \mathbf{T} \mathbf{W}_c \mathbf{T}^{*} = \mathbf{T}^{-*} \mathbf{W}_o \mathbf{T}^{-1} = \boldsymbol{\Sigma} \]
where \(\boldsymbol{\Sigma} \in \mathbb{R}^{n \times n}\) is a diagonal matrix formed from the Hankel Singular Values of the system. The system in the original coordinates
\[ \begin{pmatrix} \begin{array}{c|c} \mathbf{A} & \mathbf{B} \\ \hline \mathbf{C} & \mathbf{0} \end{array} \end{pmatrix} \]
gets transformed to
\[ \begin{pmatrix} \begin{array}{c|c} \mathbf{TAT}^{-1} & \mathbf{TB} \\ \hline \mathbf{CT}^{-1} & \mathbf{0} \end{array} \end{pmatrix} \]
in the balancing coordinate system. Although in this new reference frame the observability and controllability gramians are now diagonal and equal, the system remains \(n\)-dimensional. Assuming the balancing transformation \(\mathbf{T}\) has been computed such that the Hankel singular values in \(\boldsymbol{\Sigma}\) are sorted in decreasing order, Moore (1981) proposed to simply truncate this transformation by considering only the most jointly controllable and observable directions. Denoting by \(\boldsymbol{\Psi}_r\) and \(\boldsymbol{\Phi}_r^*\) the \(r\) leading columns of \(\mathbf{T}^{-1}\) and rows of \(\mathbf{T}\), respectively, we arrive at the balanced truncation model
\[ \begin{aligned} \dfrac{d \hat{\mathbf{x}}}{dt} & = \boldsymbol{\Phi}_r^* \mathbf{A} \boldsymbol{\Psi}_r \hat{\mathbf{x}} + \boldsymbol{\Phi}_r^* \mathbf{Bu} \\ \mathbf{y} & \simeq \mathbf{C} \boldsymbol{\Psi}_r \hat{\mathbf{x}}. \end{aligned} \]
with the original state vector being recovered from \(\mathbf{x} \simeq \boldsymbol{\Psi}_r \hat{\mathbf{x}}\). Although balanced truncation is not optimal in any norm, the \(\mathcal{H}_{\infty}\) norm of the error resulting from balanced truncation is upper bounded by twice the sum of the neglected Hankel singular values. Upper bounds using different norms (e.g. \(\mathcal{H}_2\) norm) and their mathematical proofs can be found in the text book by Antoulas (2005, see chap. 7).
Despite its desirable mathematical properties, balanced truncation suffers from a crucial limitation: its computational cost. It requires having access to both the observability and controllability gramians of the system. These are \(n \times n\) symmetric positive definite solutions to Lyapunov equations, requiring \(\mathcal{O}(n^6)\) operations to compute with standard algorithms. Balanced truncation thus becomes rapidly computationally intractable as the dimension of the system increases. This observation motivated Rowley (2005) to develop the Balanced Proper Orthogonal Decomposition (BPOD), a reformulation of balanced truncation computable solely from simulation data.
2.1.1.2 Balanced Proper Orthogonal Decomposition
Balanced Proper Orthogonal Decomposition (BPOD) was proposed by Clarence Rowley from Princepton University (USA) in his seminal paper (Rowley 2005) as a computable alternative to balanced truncation for large-scale LTI systems. While the controllability and observability gramians are solutions to Lyapunov equations, they can also be defined as
\[ \mathbf{W}_c = \int_{0}^{\infty} \exp \left( \tau \mathbf{A} \right) \mathbf{BB}^T \exp \left( \tau \mathbf{A}^T \right) \ \mathrm{d} \tau \]
and
\[ \mathbf{W}_o = \int_{0}^{\infty} \exp \left( \tau \mathbf{A}^T \right) \mathbf{C}^T \mathbf{C} \exp \left( \tau \mathbf{A} \right) \ \mathrm{d} \tau, \]
respectively.
Empirical controllability gramian – Provided the matrix \(\mathbf{A}\) is Hurwitz, the controllability gramian \(\mathbf{W}_c\) can be approximated to an arbitrary tolerance using a simple quadrature scheme. For a system with \(p\) inputs, let the input-to-state matrix be
\[ \mathbf{B} = \begin{bmatrix} \mathbf{b}_1 & \mathbf{b}_2 & \cdots & \mathbf{b}_p \end{bmatrix} \in \mathbb{R}^{n \times p}. \]
The state responses to unit impulses are given by
\[ \begin{aligned} \mathbf{x}_1(t) & = \exp \left( t \mathbf{A} \right) \mathbf{b}_1 \quad \text{Response to impulsive input } u_1(t) = \delta(t) \\ \mathbf{x}_2(t) & = \exp \left( t \mathbf{A} \right) \mathbf{b}_2 \quad \text{Response to impulsive input } u_2(t) = \delta(t) \\ \vdots \\ \mathbf{x}_p(t) & = \exp \left( t \mathbf{A} \right) \mathbf{b}_p \quad \text{Response to impulsive input } u_p(t) = \delta(t). \end{aligned} \]
The controllability gramian \(\mathbf{W}_c\) can be expressed as
\[ \mathbf{W}_c = \int_{0}^{\infty} \left( \mathbf{x}_1(\tau) \mathbf{x}_1^T(\tau) + \mathbf{x}_2(\tau) \mathbf{x}_2^T(\tau) + \cdots + \mathbf{x}_p(\tau) \mathbf{x}_p^T(\tau) \right) \ \mathrm{d} \tau. \]
If the impulse responses are given at discrete times \(t_1, t_2, \cdots, t_m\), as is typical in numerical simulations, this integral can be approximated with a simple quadrature sum. Denoting the data matrix formed from the impulse responses as
\[ \mathbf{X} = \begin{bmatrix} \mathbf{x}_1(t_1) \sqrt{w_1} & \cdots & \mathbf{x}_1(t_m) \sqrt{w_m} & \cdots & \mathbf{x}_p(t_1)\sqrt{w_1} & \cdots & \mathbf{x}_p(t_m) \sqrt{w_m} \end{bmatrix}, \]
where \(w_i\) are quadrature coefficients, the empirical controllability gramian \(\hat{\mathbf{W}_c}\) is computed as
\[ \hat{\mathbf{W}_c} = \mathbf{XX}^T. \]
If \(\mathbf{A}\) is Hurwitz, this empirical approximation converges to the true controllability gramian as \(t_m\) tends to infinity. For dissipative systems, the data matrix \(\mathbf{X}\) often tends to be low-rank, a property leveraged later on in the BPOD algorithm.
Empirical observability gramian – For LTI systems, observability is the dual of controllability. Hence, similar techniques can be employed to approximate the observability gramian \(\mathbf{W}_o\). Let the measurement matrix \(\mathbf{C}\) be
\[ \mathbf{C} = \begin{bmatrix} \mathbf{c}_1^T \\ \mathbf{c}_2^T \\ \vdots \\ \mathbf{c}_q^T \end{bmatrix}. \]
The impulse responses of the adjoint system
\[ \dfrac{d\mathbf{z}}{dt} = \mathbf{A}^T \mathbf{z} + \mathbf{C}^T \mathbf{v} \]
are defined as
\[ \begin{aligned} \mathbf{z}_1(t) & = \exp \left( t \mathbf{A}^T \right) \mathbf{c}_1 \quad \text{Response to impulsive input } v_1(t) = \delta(t) \\ \mathbf{z}_2(t) & = \exp \left( t \mathbf{A}^T \right) \mathbf{c}_2 \quad \text{Response to impulsive input } v_2(t) = \delta(t) \\ \vdots \\ \mathbf{z}_q(t) & = \exp \left( t \mathbf{A}^T \right) \mathbf{c}_q \quad \text{Response to impulsive input } v_q(t) = \delta(t). \end{aligned} \]
The observability gramian \(\mathbf{W}_o\) can be expressed as
\[ \mathbf{W}_o = \int_{0}^{\infty} \left( \mathbf{z}_1(t) \mathbf{z}_1^T(t) + \mathbf{z}_2(t)\mathbf{z}_2^T(t) + \cdots + \mathbf{z}_q(t) \mathbf{z}_q^T(t) \right) \ \mathrm{d}t. \]
Denoting the data matrix formed from the adjoint impulse responses as
\[ \mathbf{Y} = \begin{bmatrix} \mathbf{z}_1(t_1) \sqrt{w_1} & \cdots & \mathbf{z}_1(t_m) \sqrt{w_m} & \cdots & \mathbf{z}_q(t_1)\sqrt{w_1} & \cdots & \mathbf{z}_q(t_m) \sqrt{w_m} \end{bmatrix}, \]
the empirical observability gramian \(\hat{\mathbf{W}_o}\) is computed as
\[ \hat{\mathbf{W}_o} = \mathbf{YY}^T. \]
Once again, if \(\mathbf{A}\) is Hurwitz, this empirical approximation converges to the true observability gramian as \(t_m\) tends to infinity. Likewise, for dissipative systems, the data matrix \(\mathbf{Y}\) often tends to be low-rank.
The BPOD procedure – Given the data matrices \(\mathbf{X}\) and \(\mathbf{Y}\), BPOD computes a convergent approximation of balanced truncation in two steps. First, the singular value decomposition of the inner-product matrix \(\mathbf{Y}^T \mathbf{X}\) is computed, leading to
\[ \mathbf{Y}^T \mathbf{X} = \begin{bmatrix} \mathbf{U}_1 & \mathbf{U}_2 \end{bmatrix} \begin{bmatrix} \boldsymbol{\Sigma}_1 & 0 \\ 0 & \boldsymbol{\Sigma}_2 \end{bmatrix} \begin{bmatrix} \mathbf{V}_1^T \\ \mathbf{V}_2^T \end{bmatrix}. \]
The balancing modes are then computed as
\[ \boldsymbol{\Psi} = \mathbf{XV}_1 \boldsymbol{\Sigma}_1^{-\frac12} \quad \text{and} \quad \boldsymbol{\Phi} = \boldsymbol{\Sigma}^{-\frac12} \mathbf{U}_1^T \mathbf{Y}^T. \]
A proposition proved by Rowley (2005) establishes that if \(r = n\) (i.e. the gramians are full rank), the matrix \(\boldsymbol{\Sigma}\) contains the Hankel singular values of the systems, while \(\boldsymbol{\Psi}\) determines the balancing transformation. Otherwise, \(\boldsymbol{\Phi}\) and \(\boldsymbol{\Psi}\) approximate the matrices obtained from balanced truncation.
Although BPOD formally converges to balanced truncation only in the limit \(t_m \to \infty\), a wealth of empirical evidence in the literature highlights that the error incurred by BPOD benefits from the same upper bounds as the error resulting from balanced truncation. The BPOD reduced order model
\[ \begin{aligned} \dfrac{d \mathbf{x}_r}{dt} & = \mathbf{A}_r \mathbf{x}_r + \mathbf{B}_r \mathbf{u} \\ \mathbf{y}_r & = \mathbf{C}_r \mathbf{x}_r + \mathbf{Du} \end{aligned} \]
is constructed in an analogous fashion as the one resulting from balanced truncation. One major limitation of BPOD is however the need for the adjoint system. Yet, the solid mathematical grounds of balanced truncation combined with the computational efficiency of BPOD makes this procedure an excellent baseline model against which to benchmark the performances of data-driven models developed subsequently in the CONMAN project.
nekStab library
Balanced Proper Orthogonal Decomposition is currently being implemented in the development branch of our open-source library nekStab. Being based on the high-performance CFD solver Nek5000, the BPOD implementation in nekStab will enable us to apply this procedure to design reduced-order models of a wide variety of flow configurations including fully three-dimensional flows, a feast never achieved (to the best of our knowledge) in the fluid dynamics literature.
2.1.2 Sensor and actuator placement
Although reduced-order modeling is a key enabler for analysis and feedback control, actuators and sensors placement (i.e. designing \(\mathbf{B}\) and \(\mathbf{C}\)) is an upstream task having important downstream consequences. Consider the situation where \(m\) possible sensors are available albeit we only have the budget for \(p\) of them. There are
\[ \binom{m}{p} = \dfrac{m!}{\left(m-p \right)! \times p!} \]
combinations to choose from. No matter the performance index considered, sensors selection thus fundamentally is a rapidly intractable combinatorial optimization problem.
Many techniques have been proposed to approximately solve this problem greedily as an alternative to the combinatorial approach. These include Gappy POD (Willcox 2006), the Discrete Empirical Interpolation Method (Drmac and Gugercin 2016) or through pivoted QR decomposition of the POD matrix (Manohar et al. 2018). All these techniques are however primarily targetting the state reconstruction problem. Given a matrix of possible measurements \(\mathbf{C}\) such that
\[ \mathbf{y} \simeq \mathbf{SCx}, \]
these techniques aim at optimizing the selection matrix \(\mathbf{S} \in \mathbb{R}^{p \times n}\) such that the state vector \(\mathbf{x}\) is best estimated from the sensor measurements \(\mathbf{y}\) in the least-squares sense. Although of practical interest for linear inverse problems, these techniques cannot be directly applied to select sensors or actuators for an LTI dynamical system. For systems with actuation, it is necessary to consider simultaneously the placement of sensors and actuators since the most observable and controllable subspaces might be quite different. Current state-of-the-art techniques include Chen and Rowley (2011), Burns and Herdman (2021), or Manohar, Kutz, and Brunton (2022) for instance. Although we only present the ideas put forth by Manohar, Kutz, and Brunton (2022) in this technical report for the sake of conciseness, all three of the aforementionned techniques will be considered as possible baseline reference in the rest of the CONMAN project.
Reminder about LTI systems – Given a linear time-invariant dynamical system
\[ \begin{aligned} \dfrac{d\mathbf{x}}{dt} & = \mathbf{Ax} + \mathbf{Bu} \\ \mathbf{y} & = \mathbf{Cx}. \end{aligned} \]
Its observability gramian \(\mathbf{W}_o\) and controllability gramian \(\mathbf{W}_c\) are defined as
\[ \mathbf{W}_o = \int_{0}^{\infty} \exp \left( t \mathbf{A}^T \right) \mathbf{C}^T \mathbf{C} \exp \left( t \mathbf{A} \right) \ \mathrm{d}t \]
and
\[ \mathbf{W}_c = \int_{0}^{\infty} \exp \left( t \mathbf{A} \right) \mathbf{BB}^T \exp \left( t \mathbf{A}^T \right) \ \mathrm{d}t. \]
Because they depend explicitly on \(\mathbf{B}\) and \(\mathbf{C}\), the gramians are often used to evaluate the observability/controllability of a given sensor and actuator placement. One standard performance index is the \(\mathcal{H}_2\) norm of the system, measuring the average output gain over all frequencies of the input. Given the transfer function \(\mathbf{G}(s) = \mathbf{C} \left( s \mathbf{I} - \mathbf{A} \right)^{-1} \mathbf{B}\), the \(\mathcal{H}_2\) norm of the system is defined as
\[ \| \mathbf{G} \|_{\mathcal{H}_2}^2 = \dfrac{1}{4\pi^2} \int_{0}^{\infty} \mathrm{trace} \left( \mathbf{G}^H(i\omega) \mathbf{G}(i\omega) \right) \ \mathrm{d}\omega. \]
By the Plancherel theorem, it can also be defined in the time-domain using the impulse response \(y_{ij}(t) = \mathbf{C}_i \mathrm{e}^{t\mathbf{A}} \mathbf{B}_j\) as
\[ \begin{aligned} \| \mathbf{G} \|_{\mathcal{H}_2}^2 & = \int_{0}^{\infty} \mathrm{trace} \left( \mathbf{C} \mathrm{e}^{t \mathbf{A}} \mathbf{BB}^T \mathrm{e}^{t \mathbf{A}^T} \mathbf{C}^T \right) \ \mathrm{d}t = \mathrm{trace} \left( \mathbf{C} \mathbf{W}_c \mathbf{C}^T \right) \\ & = \int_{0}^{\infty} \mathrm{trace} \left( \mathbf{B}^T \mathrm{e}^{t \mathbf{A}^T} \mathbf{C}^T \mathbf{C} \mathrm{e}^{t\mathbf{A}} \mathbf{B} \right) \ \mathrm{d}t = \mathrm{trace} \left( \mathbf{B}^T \mathbf{W}_o \mathbf{B} \right). \end{aligned} \]
The \(\mathcal{H}_2\) norm of a system thus relates explicitly the measurement matrix \(\mathbf{C}\) and input-to-state matrix \(\mathbf{B}\) to each of its gramians. An alternative to the average output energy metric are the volumetric measures given by
\[ \log\det \mathbf{C} \mathbf{W}_c \mathbf{C}^T \quad \text{and} \quad \log\det \mathbf{B}^T \mathbf{W}_o \mathbf{B}, \]
corresponding to the geoemtric mean of the axes of the ellipsoid skewed by \(\mathbf{B}\) or \(\mathbf{C}\). In comparison, the trace used to defined the \(\mathcal{H}_2\) norm corresponds to the arithmetic mean.
Optimal sensor and actuator selection using balanced model reduction – As baseline model, we consider the framework put forth by Manohar, Kutz, and Brunton (2022), combining balanced truncation and greedy maximization of the \(\log\det\) objective function corresponding to a D-optimal design strategy. Given the full-order LTI system
\[ \begin{aligned} \dfrac{d \mathbf{x}}{dt} & = \mathbf{Ax} + \mathbf{Bu} \\ \mathbf{y} & = \mathbf{Cx} \end{aligned} \]
with \(\mathbf{x} \in \mathbb{R}^n\) the (high-dimensional) state vector, they first compute the (truncated) balancing transformation \(\boldsymbol{\Psi}_r\) and \(\boldsymbol{\Phi}_r\) (either using the standard balanced truncation algorithm or BPOD). This step corresponds to a model order reduction assuming the system is fully actuated and sensed.
In a second step, Manohar, Kutz, and Brunton (2022) proposed a greedy approach such that the \(\mathcal{H}_2\) norm of the model gracefully degrades as we reduce the number of sensors and/or actuators. Denoting by \(\mathbb{S}_{C}\) the sensor selection operator, they define it as being the solution to the following optimization problem
\[ \max_{\mathbb{S}_{C}} \quad \log\det \mathbb{S}_{C} \mathbf{C} \mathbf{W}_c \mathbf{C}^* \mathbb{S}_{C}^* \]
Note that low-rank approximations of the gramians of the original system can be obtained as
\[ \mathbf{W}_c \simeq \boldsymbol{\Psi}_r \boldsymbol{\Sigma}_r \boldsymbol{\Psi}_r^* \quad \text{and} \quad \mathbf{W}_o \simeq \boldsymbol{\Phi}_r \boldsymbol{\Sigma}_r \boldsymbol{\Phi}_r^*. \]
Introducing this ansatz into the optimization problem leads to
\[ \max_{\mathbb{S}_C} \quad \vert \det \left( \mathbb{S}_C \mathbf{C} \boldsymbol{\Psi}_r \right) \vert. \]
A similar problem can be formulated for the actuator selection operator \(\mathbb{S}_B\). This problem remains a combinatorial problem. In order to approximately solve it, Manohar, Kutz, and Brunton (2022) rely on their prior works (Manohar et al. 2018) and suggest to greedily construct \(\mathbb{S}_C\) by applying the pivoted QR decomposition to the matrix \(\left( \mathbf{C} \boldsymbol{\Psi}_r \right) ^*\). Their approach has been shown to scale favourably as the dimension \(n\) of the full order system increases, while its performances are comparable to direct optimization of the \(\mathcal{H}_2\) norm of the system (Chen and Rowley 2011). Finally, while the authors considered mostly open-loop sensor and actuator placements, they have also illustrated how to extend their methodology to selection sensors and actuators for closed-loop systems.
To the best of our knowledge, the code base developped by Manohar, Kutz, and Brunton (2022) is not available. Note however that the algorithm is incredibly simple as it relies on the pivoted QR decomposition. This matrix factorization is readily available in most scientific programming languages, including MATLAB, SciPy, SciLab, or Octave for instance. In the rest of this project, this baseline model will be re-implemented using the Julia programming language.
2.2 Theoretical framework and numerical tools
2.2.1 The fundamental equation of component analysis
Many practical models in engineering sciences belong to the class of generalized linear models
\[ \mathbf{y} = \mathbf{Kx} + \boldsymbol{\varepsilon} \]
with \(\mathbf{y} \in \mathbb{R}^m\), \(\mathbf{x} \in \mathbb{R}^n\) and \(\mathbf{K} \in \mathbb{R}^{m \times n}\) a linear mapping from \(\mathbf{x}\) to \(\mathbf{y}\). In system identification, it includes for instance OKID where \(\mathbf{y}\) is the response of the unknown linear system, \(\mathbf{K}\) a Toeplitz matrix constructed from the input sequence, and \(\mathbf{x}\) is the vector of unknown Markov parameters of the system. Similarly, linear stochastic estimation or the linear deconvolution problem can also be cast as generalized linear models.
Of interest to us are situations where the operator \(\mathbf{K}\) is unknown. Given training pairs \(\left( \mathbf{x}_i, \mathbf{y}_i \right)\), an ordinary least-squares regression can be formulated to identify \(\mathbf{K}\). yet, for typical engineering problems, \(\mathbf{x}\) and/or \(\mathbf{y}\) are high-dimensional vectors. Hence, we are unlikely to have sufficient data to obtain a good statistical estimate of \(\mathbf{K}\). It can however be assumed to be a low-rank operator, an assumption often verified for high-dimensional dissipative systems. A good estimate can then be obtained by solving the following rank-constrained problem
\[ \begin{aligned} \min & \quad \| \mathbf{M}^\frac12 \left( \mathbf{Y} - \mathbf{PQ}^T \mathbf{X} \right) \mathbf{W}^\frac12 \|_F^2 \\ \mathrm{sub.~to} & \quad \mathrm{rank}~\mathbf{P} = r. \end{aligned} \]
Here, \(\mathbf{X}\) and \(\mathbf{Y}\) are data matrices, $ and \(\mathbf{W}\) are symmetric positive (semi-) definite matrices and \(r\) is the desired rank of the approximation. Although non-convex, this rank-constrained problem admits a closed-form solution (Héas and Herzet 2022; Loiseau 2020). In the statistics community, this problem is also known as the fundamental equation of component analysis (De la Torre 2012).
Many tools classically used in data-driven engineering fall in this class of problem. As an example, POD analysis can be recovered from this problem by setting \(\mathbf{X} = \mathbf{Y}\) and \(\mathbf{P} = \mathbf{Q}\) leading to the celebrated low-rank matrix approximation problem whose solution is given by the truncated singular value decomposition of the data matrix \(\mathbf{X}\). Likewise, Héas and Herzet (2022) and Loiseau (2020) have shown how Dynamic Mode Decomposition (Schmid 2022) can be cast in this form leading to a more robust and well-posed DMD problem.
Along with the aforementionned modal decompositions, Jean-Christophe Loiseau and Steven L. Brunton have been able over the course of the past year to cast numerous linear system identification techniques into this framework. It includes the EigenRealization Algorithm (Juang and Pappa 1985), N4SID (Van Overschee and Moor 1993), MOESP (Verhaegen and Dewilde 1992) and CVA (Larimore 1983). Reduced-rank regression thus provides a unifying framework for a large class of data-driven linear modeling techniques, simplyfing both the implementation of the algorithm as well as the analysis of their statistical properties. A review article by Loiseau & Brunton is currently in preparation.
nekStab library
Currently, efforts are being devoted to the implementation of the reduced-rank regression algorithm in the open-source library nekStab. This implementation will enable us to fit these models using high-dimensional fluid dynamics datasets by benefitting from the high-performance capabilities of the spectral element solver Nek5000.
DataDrivenLinearModels.jl
In parallel to the implementation in nekStab, specialized for incompressible fluid dynamics simulation data obtained from Nek5000, a stand-alone Julia package DataDrivenLinearModels.jl will be developed. This package will have a large usage base than thenekStab implementation as data will simply need to be stored in standard matrix formats.
2.2.2 Submodular set function maximization
Optimal sensor and actuator selection is an important unsolved problem in control theory. Nearly every downstream control decision is affected by the selected sensors and actuators. Yet, determining optimal locations is a brute-force search through the combinatorial possibilities. Given \(m\) available sensors but a budget for only \(p\) of them, there are
\[ \binom{m}{p} = \dfrac{m!}{\left( m - p \right)! \times p!} \]
possible combinations. No matter the performance index considered, this is a rapidly intractable combinatorial problem. Various heuristics have been proposed in the literature for LTI systems, including Chen and Rowley (2011), Summers, Cortesi, and Lygeros (2016), Zhang, Ayoub, and Sundaram (2017), Manohar et al. (2018), Burns and Herdman (2021), or Manohar, Kutz, and Brunton (2022) for instance. Most of these techniques are however targetting sensor and actuator placement for open-loop systems.
Although we are eventually interested in applications of sensor and actuator placement for closed-loop systems, let us illustrate the optimization methodology developped in this project on a simpler problem. For that purpose, consider the following linear inverse problem
\[ \mathbf{y} = \mathbf{Cx} + \mathbf{v}, \]
where \(\mathbf{x} \in \mathbb{R}^n\) is the state vector of the system, and \(\mathbf{C} \in \mathbb{R}^{m \times n}\) is the observation matrix with each of its rows describing a different sensor. The vector \(\mathbf{y} \in \mathbb{R}^m\) is the measurement vector whose i-th entry corresponds to a scalar observation of the i-th sensor, i.e. \(y_i = \mathbf{c}_i^T \mathbf{x} + v_i\). These measurements are being contaminated by the sensor noise \(\mathbf{v} \sim \mathcal{N}(\mathbf{0}, \mathbf{V})\). Our objective is to determine the best subset of sensors to activate as to minimize the estimation error, subject to a constraint on the number of activations. Let \(\mathbf{z} \in \left\{ 0, 1 \right\}^{m}\) where \(z_i\) indicates whether the i-th sensor has been selected (\(z_i = 1\)) or not (\(z_i = 0\)). The effective measurement equation can be expressed as
\[ \mathbf{y}_z = \mathbf{S}_z \mathbf{Cx} + \mathbf{S}_z \mathbf{v}, \]
where \(\mathbf{S}_z \in \left\{ 0, 1 \right\}^{\vert z \vert \times m}\) is the selection operator constructed from the binary vector \(\mathbf{z}\). Adopting a Bayesian point of view, the minimum mean-squared error estimate is given by
\[ \mathbf{x} = \bar{\mathbf{x}} + \left( \boldsymbol{\Sigma}^{-1} + \mathbf{C}^T \mathbf{S}_z^T \mathbf{V}_z^{-1} \mathbf{S}_z \mathbf{C} \right)^{-1} \mathbf{C}^T \mathbf{S}_z^T \left( \mathbf{y}_z - \mathbf{S}_z \mathbf{C} \bar{\mathbf{x}} \right) \]
with \(\boldsymbol{\Sigma}\) the state covariance matrix. The matrix
\[ \mathbf{F}_z = \boldsymbol{\Sigma}^{-1} + \mathbf{C}^T \mathbf{S}_z^T \mathbf{V}_z^{-1} \mathbf{S}_z \mathbf{C} \]
is the Fisher information matrix, and corresponds to the inverse of the estimator’s error covariance matrix. Our objective is thus to maximize the spectral content of this information matrix to ensure it as well conditionned as possible, thus reducing the estimation error.
Different objective functions can be used. These are classically chosen as \(f_1(\mathbf{z}) = \mathrm{trace} \left( \mathbf{F}_z \right)\), \(f_2(\mathbf{z}) = \log\det \left( \mathbf{F}_z \right)\), or \(f_3(\mathbf{z}) = \lambda_{\min}\left( \mathbf{F}_z \right)\), corresponding to different optimal design strategies. In either case, we aim to solve the following optimization problem
\[ \begin{aligned} \max_{\mathbf{z} \in \left\{ 0, 1 \right\}^{m}} & \quad f_i(\mathbf{z}) \\ \mathrm{subjecto~to} & \quad \vert \mathbf{z} \vert = p. \end{aligned} \]
Despite being a combinatorially intractable problem, one can leverage some properties of the objective functions to design greedy algorithms with approximation guarantees (Shamaiah, Banerjee, and Vikalo 2010). In particular, if the sensor noise covariance matrix is diagonal (i.e. uncorrelated noise), the function \(f_2(\mathbf{F}_z) = \log\det(\mathbf{F}_z)\) is a monotone submodular set function:
- monotone – Adding a new sensor can only improve the quality of our estimate.
- submodular – As the number of selected sensors increases, the benefit of adding yet another sensor decreases. It is the diminishing return property.
Maximizing a monotone submodular set function under a cardinality constraint can be done using extremely simple greedy algorithms (Minoux, 1978) with approximation guarantees. Current efforts by the Principal Investigator of the project are devoted to generalizing the results of Shamaiah, Banerjee, and Vikalo (2010) to arbitrary noise covariance matrix as well as implementing state-of-the-art submodular function maximization algorithms in the open-source library SubmodularOpt.jl.
SubmodularOpt.jl
SubmodularOpt.jl is a Julia package currently under active development by members of the CONMAN project. It provides efficient implementations of standard submodular function maximization algorithms, including the Lazy Greedy (Minoux, 1978) or the Lazier than Lazy Greedy (Mirzasoleiman et al., 2015) algorithms for maximizing a monotone submodular set function with cardinality constraint. Inspired from the concave-convex procedure for minimizing a difference of convex functions, SubmodularOpt.jl also provides an implementation of the SubSup procedure (Iyer & Bilmes, 2013) for the minimization/maximization of the difference of two submodular functions.
2.3 Problems encountered and actions taken for scientific communication and dissemination
2.3.1 Human ressources
Postdocotral researcher – Recruitment of the postdoctoral researcher was initially planned for the first quarter of 2023. Unfortunately, despite the numerous applications, no applicant had the sufficient background in numerical simulation, reduced-order modeling or linear dynamical system theory to ensure rapid progress of the project. Decision has been made to postpone this recruitment to Fall 2023. Doing so gave the opportunity to Simon Kern to join the project. Simon defended his Ph. D. work entitled Linear and non-linear dynamics of non-autonomous flows at the Royal Institute of Technology (KTH, Sweden) early June 2023 under the joint supervision of Ardeshir Hanifi and Dan Henningson. He will join us in September.
2.3.2 Scientific communication and dissemination
This section summarizes some of the actions taken for scientific communication and dissemination of the results related (work package WP5). So far, these actions are quite limited. This results from the delay in the recruitment of the postdoctoral researcher (initially expected to start during the first quarter of 2023 but effectively starting early september) and the fact the project itself only started a few months ago.
2.3.2.1 Journal Articles
Although the project officially started only early January 2023, three papers are already in preparation. A brief overview of each is given below.
Unifying framework for data-driven linear modeling with rank constraints
| Authors | J.-Ch. Loiseau (Arts et Métiers, France) & S. L. Brunton (UW, USA) |
| Tentative journal | SIAM Journal on Applied Dynamical Systems |
| Expected submission date | Fall 2023 |
As discussed in Section 2.2.1, many data-driven linear modeling frameworks for high-dimensional systems can be recast as the following reduced-rank regression problem
\[ \begin{aligned} \min & \quad \dfrac12 \| \mathbf{M}^\frac12 \left( \mathbf{Y} - \mathbf{PQ}^T \mathbf{X} \right) \mathbf{W}^{\frac12} \|_F^2 \\ \mathrm{subject~to} & \quad \mathrm{rank}(\mathbf{P}) = r. \end{aligned} \]
This includes Proper Orthogonal Decomposition (Berkooz, Holmes, and Lumley 1993) for which \(\mathbf{X} = \mathbf{Y}\) and \(\mathbf{M} = \mathbf{W} = \mathbf{I}\), Dynamic Mode Decomposition (Schmid 2022) where \(\mathbf{y}_t = \mathbf{x}_{t+1}\), or Linear Stochastic Estimation to name only the most popular in the fluid dynamics community. Likewise, numerous linear system identification techniques used in the control community can be derived from this rank-constrained optimization problem, e.g. EigenRealization Algorithm (Juang and Pappa 1985), N4SID (Van Overschee and Moor 1993) or Canonical Variate Analysis (Larimore 1983).
This article delves into this fundamental equation of component analysis, its derivation and the proof of its solution, and illustrates how classical techniques in mechanical engineering and control systems can be derived from it as special cases. In the process, we show how the performances and robustness of DMD-like techniques can be improved by correctly formulating these problems into this framework. Moreover, as a consequence of these techniques being special instances of a more general problem, we can also draw parallels and connections between them. In particular, we elucidate the connection between the modes resulting from Proper Orthogonal Decomposition and those from Dynamic Mode Decomposition, a long standing and unsatisfactrily answered question in the fluid dynamics literature.
SubmodularOpt.jl : a Julia package for greedy maximization of submodular set functions
| Author | J.-Ch. Loiseau (Arts et Métiers, France) |
| Tentative journal | Journal of Open Source Softwares |
| Expected submission date | Fall 2023 |
SubmodularOpt.jl is one of the main packages of this project. It is being openly developped using GitHub. This journal article, to be submitted to the Journal of Open Source Softwares, presents its core components as well as an example of sensor placement problem it can solve. Following the FAIR principles (Free, Accessible, Inclusive, Reproducible) for open and reproducible science, the corresponding release of SubmodularOpt.jl will be archived for long-term preservation on platforms such as Zenodo, granting it with a citable DOI.
2.3.2.2 Conferences and Seminars
Since CONMAN started only six months ago, members of the project did not have the occasion (yet) to present their work in conferences. Plans are being made however to send two contributions to the APD Division of Fluid Dynamics 76th Annual Meeting to be held in Washington DC, USA late November 2023. These contributions include:
- A unifying framework for data-driven linear modeling with rank constraints. J.-Ch. Loiseau and S. L. Brunton.
- D-optimal sensor placement with submodular optimization. J.-Ch. Loiseau.
Additionally, two seminars in French research laboratories are currently planned for Fall 2023:
- Tips and tricks for Koopman analysis. J.-Ch. Loiseau. ENSTA Paris, institut Polytechnique de Paris, France.
- Data-driven reduced-order modeling. J-Ch. Loiseau. Laboratoire de Physique, ENS Lyon, France.
As the project continues, great care will be taken to provide the postdoctoral research and Ph. D. candidate with opportunities to present their works at international conferences and seminars.
2.4 What’s next?
During the next period (July-December 2023), efforts will be focused on wrapping up the implementations of the baseline models in our different libraries. A particular attention will be devoted to validation & verification of the results for these baseline models as well as to their reproducibility. This includes:
- Validation and verification of the BPOD implementation in
nekStabagainst applications available in the literature. - Setting up the continuous integration feature for all the packages developed on GitHub.
- Designing informative unit tests for all the packages to minimize the chances of an unwanted behaviour slipping in.
- Setting up our reproducibility pipeline for results presented in journal articles using the Quarto publishing system and showyourwork.
From a purely scientific point of view, we expect the following:
- Putting the final piece to the unifying framework for data-driven linear modeling with rank constraint, including publication of the journal article.
- Formalising the proof of the submodularity (or submodular relaxation) of the sensor/actuator placement problem, both in the open-loop and closed-loop setups.
- Obtaining the first benchmark results for the submodular sensor placement vs. the approach proposed by Manohar, Kutz, and Brunton (2022) on the linearized Ginzburg-Landau equation.
It needs to be emphasized finally that a large portion of the time of J.-Ch. Loiseau will be dedicated to the education of the Ph. D. candidate arriving early Fall 2023. This includes familiarizing the Ph. D. candidate with both the theoretical and numerical tools on which this project relies as well as to the French Higher Education and Research ecosystem.