Full-order optimal compensators for flow control: the multi-input case

TL;DR

This paper develops an iterative algorithm for designing full-order optimal flow controllers with multiple actuators and sensors, avoiding complex Riccati equations, and demonstrates its effectiveness and scalability in large dynamical systems.

Contribution

It extends the ADA algorithm to multi-input, coupled, and robust H-infinity controllers, enabling direct design without model reduction for large flow systems.

Findings

Algorithm outperforms existing methods in convergence speed.

Scalable with the number of actuators in large systems.

Effective in complex flow control scenarios like Kuramoto-Sivashinsky.

Abstract

Flow control has been the subject of numerous experimental and theoretical works. In this numerical study, we analyse full-order, optimal controllers for large dynamical systems in presence of multiple actuators and sensors. We start from the original technique proposed by Bewley et al. (2016), the adjoint of the direct-adjoint (ADA) algorithm. The algorithm is iterative and allows bypassing the solution of the algebraic Riccati equation associated with the optimal control problems, typically unfeasible for large systems. We extend ADA into a more generalized framework that includes the design of multi-input, coupled controllers and robust controllers based on the H_{\infty} framework. The full-order controllers do not require any preliminary step of model reduction or low-order approximation: this feature allows to pre-assess the optimal performances of an actuated flow without relying on any estimation process or further hypothesis. We show that the algorithm outperforms analogous technique, in terms of convergence performances considering two numerical cases: a distributed system and the linearized Kuramoto-Sivashinsky equation, mimicking a full three-dimensional control setup. For the ADA algorithm we find excellent scalability with the number of inputs (actuators) in terms of convergence to the solution, making the method a viable way for full-order controller design in complex settings.

Tables

Table 1. Table 1: The table summarizes the main features of the algorithms compared in Sec.  4.3 : the column IC indicates the initial conditions used; the third column Grad. indicates the dimensions of the gradients used during the process of update; in the column Loops , the number of iteration-loops required for computing the control gains of dimension n × m 𝑛 𝑚 n\times m is indicated. In the last column, the characterizing features of the methods are highlighted.

Algorithm	IC	Grad.	Loops	Notes
DA	Random at each loop	$m\times n_t$	$n$	$𝑲$ is solution of a linear system.
ADA	$-{𝐑}_{i,i}^{-1}{𝑩}_i$	$m\times n_t$	$m$	Each row of $𝑲$ is solution of one iteration-loop.
MR	Random at each iter.	$m\times n$	$1$	$𝑲$ is solution of one iteration-loop, $\forall m$.

Equations

$$\mathcal{J}=\frac{1}{2}{\int }_0^T\left({\text{z}}^H\text{z}+{\text{u}}^H\mathbf{R}\text{u}\right)dt.$$

$$\dfrac{d\boldsymbol{q}}{dt}&=\boldsymbol{A}\boldsymbol{q}+\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\boldsymbol{\textsf{u}},\qquad with\qquad\boldsymbol{q}(0)=\boldsymbol{q}_{0},\\ \boldsymbol{\textsf{z}}=\boldsymbol{C}_{\boldsymbol{\textsf{z}}}\boldsymbol{q},$$

$$\mathcal{J}=\frac{1}{2}{\int }_0^T\left({\mathbf{q}}^H\mathbf{Q}\mathbf{q}+{\text{u}}^H\mathbf{R}\text{u}\right)dt+\frac{1}{2}\mathbf{q}(T{)}^H{\mathbf{Q}}_T\mathbf{q}(T)$$

$$\tilde{\mathcal{J}}$$

$$\dfrac{d\boldsymbol{q}}{dt}&=\boldsymbol{A}\boldsymbol{q}+\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\boldsymbol{\textsf{u}},\quad\qquad with\quad\boldsymbol{q}(0)=\boldsymbol{q}_{0},\\ \dfrac{d\boldsymbol{p}}{dt}=-\boldsymbol{A}^{H}\boldsymbol{p}-\boldsymbol{Q}\boldsymbol{q},\quad\,\,with\quad\boldsymbol{p}(T)=0,\\ \dfrac{\partial\mathcal{\tilde{J}}}{\partial\boldsymbol{\textsf{u}}}=\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}\boldsymbol{p}+\boldsymbol{R}\boldsymbol{\textsf{u}}.$$

$$\text{u}(t{)}^{i+1}=\text{u}(t{)}^i-{\beta }^i{\left(\frac{\partial \tilde{\mathcal{J}}}{\partial \text{u}}\right)}^i,$$

$$\dot{\left(\begin{array}[]{c}\boldsymbol{q}\\ \boldsymbol{p}\end{array}\right)}=\left[\begin{array}[]{cc}\boldsymbol{A}&-\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}\\ -\boldsymbol{Q}&-\boldsymbol{A}^{H}\end{array}\right]\left(\begin{array}[]{c}\boldsymbol{q}\\ \boldsymbol{p}\end{array}\right).$$

$${\mathbf{A}}^H\mathbf{X}+\mathbf{X}\mathbf{A}-\mathbf{X}{\mathbf{B}}_{\text{u}}{\mathbf{R}}^{-1}{\mathbf{B}}_{\text{u}}^H\mathbf{X}+\mathbf{Q}=0,$$

$$\mathbf{K}=-{\mathbf{R}}^{-1}{\mathbf{B}}_{\text{u}}^H\mathbf{X}.$$

$$\left[{\text{u}}_0^1{\text{u}}_0^2\dots {\text{u}}_0^n\right]={\mathbf{K}}_{(1\times n)},{\left[{\mathbf{q}}_0^1{\mathbf{q}}_0^2\dots {\mathbf{q}}_0^n\right]}_{(n\times n)},$$

$$\boldsymbol{J}=\left[\begin{array}[]{rcccr}\boldsymbol{0}&&&\boldsymbol{I}\\ -\boldsymbol{I}&&&\boldsymbol{0}\end{array}\right],$$

$${\mathbf{T}}^H\mathbf{J}+\mathbf{J}\mathbf{T}=0.$$

$$\boldsymbol{T}=\left[\begin{array}[]{cc}\boldsymbol{A}&-\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}\\ -\boldsymbol{Q}&-\boldsymbol{A}^{H}\end{array}\right].$$

$$\Omega (t)=\tilde{\mathbf{x}}(t{)}^H\mathbf{J}\mathbf{x}(t)=\tilde{\mathbf{p}}(t{)}^H\mathbf{q}(t)-\tilde{\mathbf{q}}(t{)}^H\mathbf{p}(t).$$

$$\tilde{\boldsymbol{p}}_{0}^{H}\boldsymbol{q}_{0}&=\tilde{\boldsymbol{q}}_{0}^{H}\boldsymbol{p}_{0},\\ \boldsymbol{\textsf{u}}=\boldsymbol{K}\boldsymbol{q}_{0}=(-\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H})\boldsymbol{p}_{0}.$$

$$\frac{d\mathbf{q}}{dt}$$

$$\mathcal{J}$$

$$\dot{\hat{\boldsymbol{q}}}&=\boldsymbol{A}\hat{\boldsymbol{q}}+\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\boldsymbol{\textsf{u}}-\boldsymbol{L}(\boldsymbol{\textsf{y}}-\hat{\boldsymbol{\textsf{y}}}),\\ \hat{\boldsymbol{\textsf{y}}}=\boldsymbol{C}_{\boldsymbol{\textsf{y}}}{\hat{\boldsymbol{q}}}.$$

$$&\boldsymbol{A}\boldsymbol{Y}+\boldsymbol{Y}\boldsymbol{A}^{H}-\boldsymbol{Y}\boldsymbol{C}_{\boldsymbol{\textsf{y}}}^{H}\boldsymbol{G}^{-1}\boldsymbol{C}_{\boldsymbol{\textsf{y}}}\boldsymbol{Y}+\boldsymbol{B}_{\mathsf{\textsf{d}}}\boldsymbol{W}\boldsymbol{B}_{\mathsf{\textsf{d}}}^{H}=0,\\ \boldsymbol{L}=-\boldsymbol{Y}\boldsymbol{C}_{\boldsymbol{\textsf{y}}}^{H}\boldsymbol{G}^{-1}.$$

$$\mathcal{N}\left(\mathbf{p}\left(\tilde{\text{y}}\right),\tilde{\text{y}}\right)=\frac{1}{2}{\int }_0^T\left({\mathbf{p}}^H\mathbf{W}\mathbf{p}+{\tilde{\text{y}}}^H\mathbf{G}\tilde{\text{y}}\right)dt$$

$$\mathcal{J}=\frac{1}{2}{\int }_0^T\left({\mathbf{q}}^H\mathbf{Q}\mathbf{q}+{\text{u}}^H\mathbf{R}\text{u}-{\gamma }^{-2}{\text{d}}^H\mathbf{W}\text{d}\right)dt,$$

$$\frac{d\mathbf{q}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$\text{u}=-{\mathbf{R}}^{-1}{\mathbf{B}}_{\text{u}}^H\mathbf{p};$$

$$\text{d}={\gamma }^2{\mathbf{W}}^{-1}{\mathbf{B}}_{\text{d}}^H\mathbf{p}.$$

$$\begin{array}[]{c}\boldsymbol{T}\end{array}=\left[\begin{array}[]{ccc}\boldsymbol{A}&&\gamma^{2}\boldsymbol{B}_{\mathsf{\textsf{d}}}\boldsymbol{W}^{-1}\boldsymbol{B}^{H}_{\mathsf{\textsf{d}}}-\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}\\ -\boldsymbol{Q}&&-\boldsymbol{A}^{H}\end{array}\right].$$

$$&\boldsymbol{A}^{H}\boldsymbol{X}+\boldsymbol{X}\boldsymbol{A}-\boldsymbol{X}\left(\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}-\gamma^{2}\boldsymbol{B}_{\mathsf{\textsf{d}}}\boldsymbol{W}^{-1}\boldsymbol{B}_{\mathsf{\textsf{d}}}^{H}\right)\boldsymbol{X}+\boldsymbol{Q}=0,\\ \boldsymbol{K}=-\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}\boldsymbol{X},\\ \boldsymbol{Y}=\gamma^{2}\boldsymbol{W}^{-1}\boldsymbol{B}_{\mathsf{\textsf{d}}}^{H}\boldsymbol{X}.$$

$$&\boldsymbol{q}_{\boldsymbol{K},0}=-\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H},\\ \boldsymbol{q}_{\boldsymbol{Y},0}=\gamma^{2}\boldsymbol{W}^{-1}\boldsymbol{B}_{\mathsf{\textsf{d}}}^{H}.$$

$${\mathbf{K}}^{i+1}={\mathbf{K}}^i-{\beta }^i{\left({\nabla }_{\mathbf{K}}\tilde{\mathcal{J}}\right)}^i,$$

$$\beta =\frac{{\left(d{\mathcal{J}}^i\right)}^{T}d{\mathcal{J}}^i}{{\left(d{\mathcal{J}}^i\right)}^{T}d{\mathcal{J}}^{i+1}}.$$

$$\begin{array}[]{c}\mathbf{A}\end{array}=\left[\begin{array}[]{rrr rrr rrr}1&-4&0&&&&&&\\ &-1&-4&1&&&&&\\ &&1&-4&-2&&&&\\ &&&\ddots&\ddots&\ddots&&&\\ &&&&1&-4&-2&&\\ &&&&&8&-8&-1&\\ &&&&&&0&-2&1\end{array}\right].$$

Full text

We extend ADA into a more generalized framework that includes the design of multi-input, coupled controllers and robust controllers based on the $\mathcal{H}_{\infty}$ framework. The full-order controllers do not require any preliminary step of model reduction or low-order approximation: this feature allows to pre-assess the optimal performances of an actuated flow without relying on any estimation process or further hypothesis.

We show that the algorithm outperforms analogous technique, in terms of convergence performances considering two numerical cases: a distributed system and the linearized Kuramoto-Sivashinsky equation, mimicking a full three-dimensional control setup. For the ADA algorithm we find excellent scalability with the number of inputs (actuators) in terms of convergence to the solution, making the method a viable way for full-order controller design in complex settings.

††journal: Theoretical and Computational Fluid Dynamics

1 Introduction

Flow control based on linear strategies has been applied to a large variety of flows in the last decades, see 1990:aber:tema , arfm2007-kim-bewley , fabbiane2014adaptive , sipp2016linear . The interest in such applications has been invigorated by the possible, numerous outcomes, ranging from drag reduction to acoustic emission mitigation. From a physical point of view, one of the main assumptions is that the modification of the coherent structures can be achieved by properly acting on the flow, both in the presence of small amplitude perturbations or large coherent structures at higher Reynolds number. If a linear approximation is valid, rigorous methods for linear control can be applied.

In this paper we focus on active control. By definition active control is characterized by energy input in the system to be controlled, by means of actuators. The system to be controlled will often be referred to as plant. The presence of actuators allows more flexibility in the control design: sensors can be introduced such that the actuation is optimized with respect to an objective function.

An example of a plant to be controlled using active control strategies is given in Fig 1, where the sketch of a boundary layer flow developing on a flat plate is shown, including $m$ actuators $\boldsymbol{B}_{\boldsymbol{\textsf{u}}}$ and $p$ sensors $\boldsymbol{C}_{z}$. We seek for a control law $\boldsymbol{\textsf{u}}(t)$, feeding the $m$ actuators. A possible way, it is to define a control kernel $\boldsymbol{K}\in\mathbb{R}^{m\times n}$ such that the control signal $\boldsymbol{\textsf{u}}(t)$ – defined in a time interval $t=[0,\,T]$ – is proportional to the state vector $\boldsymbol{q}\in\mathbb{R}^{n}$, i.e. $\boldsymbol{\textsf{u}}(t)=\boldsymbol{K}\boldsymbol{q}$. The controller is designed for fulfilling a target; in linear quadratic regulators (LQR) the following objective function to be minimized can be defined as follow

[TABLE]

The first term includes the signals recorded at the location of the sensors $\boldsymbol{C}_{\boldsymbol{\textsf{z}}}$, as $\boldsymbol{\textsf{z}}=\boldsymbol{C}_{\boldsymbol{\textsf{z}}}\boldsymbol{q}$. The matrix $\boldsymbol{R}$ contains the control penalties for tuning the control effort (see for instance semeraro2011feedback , semeraro2013transition and reference therein).

A difficulty that arises in fluid mechanics is related to the dimensions of the dynamical systems: although it is relatively common to deal with numerical simulations with a number of degrees of freedom $n>10^{6}$, control design tools become infeasible for much smaller dimensions. A classical approach to circumvent this limitation consists in replacing the full-order system with reduced-order models capturing the essential dynamics of the system. This methodology is sometimes called reduce-then-design.

1.1 Full-order control design

Despite the necessity of reduced-order modelling for real application, the design of full-order controllers can be of great interest. Indeed, the reduced-order modelling is often performed in open-loop, thus discarding states of the system that might be reachable in the presence of a controller; the un-modelled states can lead to inconsistencies in the control design. Moreover, when designing a controller based on reduced-order models, an estimator is often necessary: this leads to a quite larger parameters space to be explored, that includes numerous choices such as the location and the spatial distribution of the sensors/actuators pairs. In this sense, it is not always straightforward to assess whether a system is controllable and – when possible – if the achieved performances are the best possible due to these numerous choices that must be done during the control design. For this reason, we believe that it is of interest to pre-assess the performance of a controller by minimizing the number of these choices using optimal full-order controllers that do not require any preliminary step of model reduction or estimation processes.

Within the optimal control framework, the control gain $\boldsymbol{K}$ can be obtained as solution of the algebraic Riccati equation (ARE), lewis ; however, already for systems of dimensions $n\approx 10^{3}-10^{4}$ the direct method is computationally intractable 2008:benn:li:penz . Neglecting the low-order design, alternative methods for the full-order design were proposed in the last decades, based on the solution of surrogate systems of equations, iterative procedures or algorithms exploiting the inherent sparsity of the considered systems. A common trait of these techniques is their feasibility when the number of actuators $m\ll n$. A classic technique is the Chandrasekhar method: the full-order ARE is replaced by a set of partial differential equations. Alternatively, the projection on low-order Krylov spaces by means of Arnoldi-type processes (see 2004:benn and references therein) and Newton methods have been proposed (1991:bank:ito ). An efficient solution method for this set of equations was proposed in akhtar2015using , where long-time integrators are used in combination with reduced-order models based on proper orthogonal decomposition (POD) for the control of the Karman street developing behind a cylinder. Full-order kernels were computed for the control of the plane channel by hogberg:bewley:henning:03:a , based on the works by bamieh2002 ; in this method, optimal control kernels with a compact support in physical space are designed in the limit of parallel flow, by solving an optimal control problem for each wavenumber pair in Fourier space, independently. A Fourier anti-transform allows reconstructing the kernel in physical domain. The procedure was applied to weakly spatially developing flows by chevalier:hoepffner:akervik:henning:07 and monokrousos2008dns , although the main drawback is the distribution of sensors/actuators pair that cannot be localized in this approach. Methods based on the pole placement are discussed, as the minimum-energy control (MCE), already introduced in 2016:bewley:luchini:pralits and bewley2007minimal ; in this limit, the control kernel is defined starting from the unstable adjoint model of the system. Examples are given by carini2015feedback for the control of vortex shedding behind a cylinder, and carini2016global for the control of the instabilities rising in the wake of a thick plate at higher Reynolds number. In the latter contribution, the linearization is performed around different mean flows, computed as RANS and U-RANS solutions.

1.2 Iterative methods and present investigation

Finally, iterative methods can be considered. The standard approach consists of an iteration based on the direct equation describing the system dynamics and its adjoint (see luchini2014adjoint and citations therein). The aim of this technique is the identification of the optimal control law $\boldsymbol{\textsf{u}}(t)$, and it can be generalized to nonlinear settings as well as receding horizon control problem (bewley2001dns , cherubini2013nonlinear , kim2014adjoint ).

However, the control signal in the standard direct-adjoint iteration (hereafter indicated with DA) depends on the initial condition from which the trajectory emanates. This limitation is tackled in Bewley et al. 2016:bewley:luchini:pralits , where the adjoint of the direct-adjoint (ADA) algorithm is introduced, a method capable of identifying the kernel gain $\boldsymbol{K}$. Examples of this application can be found for the optimal control of the von Kármán street developing past a cylinder in pralits2010riccati or the control of Tollmien-Schlichting waves developing in a two-dimensional boundary layer flow in semeraro2013riccati .

From the theoretical point of view, the ADA algorithm replaces the original direct-adjoint optimization with the corresponding sensitivity analysis by considering the adjoint of the entire problem: this procedure changes an optimization problem of size $n$, the number of states, into a problem of size $m$, the number of inputs. The algorithm is detailed in Sec. 2, where the original formulation is revisited. In this article we focus on two aspects of the algorithm. First, we thoroughly reconsider the design of a multi-input controller based on ADA; in Sec. 3, we extend the algorithm to the multi-input case, by considering the centralized control design glad2000control . This allows also the application of the same framework to the robust optimal control, $\mathcal{H}_{\infty}$ zhou:doyle:glover:02 . Moreover, we assess the performance of the algorithm from the convergence point of view; we compare the multi-input ADA algorithm is compared to an analogous iterative methodology introduced in maartensson2012gradient , based on stochastic gradient (Sec. 4.3). The comparison makes use of a simple distributed system. In Sec. 5, we apply the multi-input ADA to a two-dimensional version of the Kuramoto-Sivashinsky (KS) equation. The work finalizes with conclusions in Sec. 6.

2 Linear Quadratic Regulator (LQR): iterative solution and Riccati equation

In this section, we concisely state the optimal control problem and derive the time-continuous algebraic Riccati equation (CARE). For a deeper discussion, we refer to the specialized literature lewis . We introduce the input-output equations

[TABLE]

where $\boldsymbol{A}\in\mathbb{Rb}^{n\times n}$ is the system matrix of dimensions $n$, the degrees of freedom. In this work, we consider time-continuous, spatially-discretized systems. The spatial distribution of the actuators is described by the matrix $\boldsymbol{B}_{\boldsymbol{\textsf{u}}}\in\mathbb{R}^{n\times m}$. The variable z is scalar and represent a time signal. If the system arises from the discretization of the Navier-Stokes equations linearized around a baseflow or a mean flow, the vector $\boldsymbol{q}(t)\in\mathbb{R}^{n}$ represents the state of the fluid system. We want to identify a control signal $\boldsymbol{\textsf{u}}(t)$ such that

[TABLE]

is minimized. The matrices $\boldsymbol{Q}\geq 0\in\mathbb{R}^{n\times n}$, $\boldsymbol{R}>0\in\mathbb{R}^{m\times m}$ and $\boldsymbol{Q}_{T}\geq 0\in\mathbb{R}^{n\times n}$ contain weights in the entries. The matrix $\boldsymbol{Q}$ can be low-rank, for instance by defining it based on the sensors $\boldsymbol{C}_{\boldsymbol{\textsf{z}}}$ (see Fig. 1 and Eq. 1). In the following, the matrices $\boldsymbol{Q}$ and $\boldsymbol{R}$ are assumed diagonal, and the final condition $\boldsymbol{Q}_{T}$ is set null. These assumptions do not lead to any loss of generality.

The solution of the control problem in Eq. 2–3 can be obtained by defining an augmented cost function $\mathcal{\tilde{J}}$

[TABLE]

Applying integration by parts, the following system of equations is cast

[TABLE]

The equation for the adjoint state $\boldsymbol{p}(t)\in\mathbb{R}^{n}$ is obtained by zeroing the gradient ${\partial\mathcal{\tilde{J}}}/{\partial\boldsymbol{q}}$ and is integrated backward in the interval $t\in[T,0]$. The matrix $\boldsymbol{A}^{H}$ denotes the adjoint operator, satisfying the inner-product $\langle\boldsymbol{A}\boldsymbol{q},\boldsymbol{p}\rangle=\langle\boldsymbol{q},\boldsymbol{A}^{H}\boldsymbol{p}\rangle$. The unknown of the system is the control signal $\boldsymbol{\textsf{u}}(t)$ in the time interval $t\in[0,T]$. The solution can be approximated by iteration: at each step of the direct-adjoint iteration the control signal is updated as

[TABLE]

using 5$(c)$. The resulting control signal is optimal for a given initial condition $\boldsymbol{q}_{0}$. A gradient descent algorithm can be used for defining the step $\beta$ 2007:numericalrecipes . We will refer to this technique as direct-adjoint (DA) iteration in what follows.

A direct solution of the control problem is obtained by solving the associated Riccati equation. The system in Eq. 5 can be written in matrix form as

[TABLE]

Assuming the relation $\boldsymbol{p}=\boldsymbol{X}\boldsymbol{q}$, the following Riccati equation is obtained

[TABLE]

for linear time-invariant system, in the steady case limit $\dot{\boldsymbol{X}}=0$. The matrix $\boldsymbol{X}\in\mathbb{R}^{n\times n}$ is the solution of the Riccati equation and it is positive-definite and symmetric. The control signal $\boldsymbol{\textsf{u}}(t)$ is proportional to the state $\boldsymbol{q}$ as $\boldsymbol{\textsf{u}}(t)=\boldsymbol{K}\boldsymbol{q}$; the constant control kernel $\boldsymbol{K}\in\mathbb{R}^{m\times n}$ is

[TABLE]

For large-system $n>10^{3}$ the direct solution of the Riccati equation is not feasible, due to computational costs of order $O(n^{3})$ – regardless of the structure of the system matrix $\boldsymbol{A}$ – and storage requirements which are at least of order $O(n^{2})$, see 2004:benn . A viable alternative is represented by iterative methods for the computation of the control kernel $\boldsymbol{K}$.

2.1 Adjoint of the Direct-Adjoint (ADA) algorithm for the solution of the LQR problem

The DA iteration does not allow the direct computation of the optimal control kernel $\boldsymbol{K}\in\mathbb{R}^{n}$; indeed, the unknowns of the problem for a given initial condition $\boldsymbol{q}_{0}^{i}$ is the control signal $\boldsymbol{\textsf{u}}(t)^{i}$ defined in $t\in[0,T]$, with $T$ the final time of optimization. However, as observed by 2016:bewley:luchini:pralits , the following linear system can be formed

[TABLE]

where $n$ different solutions $\boldsymbol{\textsf{u}}(t)$ of the optimal control problem emanating from $n$ linearly independent initial conditions $\boldsymbol{q}_{0}^{i}$ are used. Thus, $n$ iterations need to be solved, each with a different initial condition. In particular, the known vector on the left-hand side can be formed by taking the corresponding control signal at $t=0$. The columns of the matrix on the right-hand side are represented by the initial conditions of each of the $n$ iterative loops. The last step consists of the solution of a linear system of dimensions $n\times n$; so, for large system this iterative scheme is again unfeasible.

However, as elucidated in pralits2010riccati and 2016:bewley:luchini:pralits , one may drastically reduce the computational costs of the problem by analysing the sensitivity with respect to the initial condition using the adjoint of the DA system. The original demonstration of the algorithm makes use of integration by-parts and it is reported in the appendix A for sake of completeness. Here, we propose an alternative version based on the properties of the Hamiltonian systems. Introducing the symplectic matrix

[TABLE]

the following property is fulfilled

[TABLE]

The matrix $\boldsymbol{T}\in\mathbb{R}^{2n\times 2n}$ for the control problem in Eq 2 is defined as

[TABLE]

Let us introduce the state {\boldsymbol{x}}=\left(\begin{array}[]{c}{\boldsymbol{q}}\\ {\boldsymbol{p}}\end{array}\right) as solution of the system $\dot{{\boldsymbol{x}}}=\boldsymbol{T}{\boldsymbol{x}}$, with initial condition {\boldsymbol{x}_{0}}=\left(\begin{array}[]{c}{\boldsymbol{q}}_{0}\\ {\boldsymbol{p}}_{T}\end{array}\right) and its adjoint state \tilde{\boldsymbol{x}}=\left(\begin{array}[]{c}\tilde{\boldsymbol{q}}\\ \tilde{\boldsymbol{p}}\end{array}\right), solution of the system $\dot{\tilde{\boldsymbol{x}}}=\boldsymbol{T}^{H}\tilde{\boldsymbol{x}}$, with initial condition {\tilde{\boldsymbol{x}}_{0}}=\left(\begin{array}[]{c}\tilde{\boldsymbol{q}}_{0}\\ \tilde{\boldsymbol{p}}_{T}\end{array}\right). Note that the state ${\tilde{\boldsymbol{x}}}$ is the adjoint of the direct state ${{\boldsymbol{x}}}$ with respect to the symplectic product $\Omega(t)$, defined as

[TABLE]

In Hamiltonian systems, this product is constant $\forall t$; including the boundary conditions, we obtain $\Omega(t)=0$. The relation 11 can be compared with the optimality condition at $t=0$

[TABLE]

We can notice that, if we introduce as initial condition of the dual system $\tilde{\boldsymbol{q}}_{0}$ one row of $-\boldsymbol{R}^{-1}\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}$, the resulting adjoint solution $\tilde{\boldsymbol{p}}_{0}^{H}$ corresponds to one row of $\boldsymbol{K}$. Thus, the solution of the Riccati problem is obtained without solving the algebraic equation Eq. 7. In principle, the exact solution is obtained for $T\rightarrow\infty$. In practice, a sufficiently long time-window for the optimization guarantees convergence towards the optimal solution. Strictly speaking, the resulting problem is not anymore an optimal control problem, despite it makes use of the same equations.

3 The multi-input version of ADA

In this section we discuss how to extend the ADA algorithm to multivariable systems, i.e. the systems characterized by the presence of multiple inputs and outputs. Within the context of flow control, this is typically the case of full, three-dimensional setups extended in the spanwise direction. An example is sketched in Fig. 1. When only one actuator/sensor pair is considered, we usually refer to a Single-Input-Single-Output (SISO) system; vice versa, the opposite case is when multiple inputs/outputs are introduced, i.e. the Multi-Input-Multi-Output (MIMO) case. In the forthcoming, we consider the MIMO case. The main issues for the MIMO design are introduced in Sec. 3.1. The extension of the ADA algorithm for the solution of the optimal control is discussed in section Sec. 3.2; the dual iteration for the estimation is briefly introduced in Sec. 3.3. Finally, the robust, optimal control design is revisited in Sec. 3.4, where the multivariable extension of ADA is adapted for the $\mathcal{H}_{\infty}$ design (skogestad05 ).

3.1 Multivariable systems

From the design point-of-view, a multivariable system is potentially characterized by cross coupling between inputs and outputs. This cross-coupling is the root of difficulties in multivariable control glad2000control : indeed, a change in one input can affect multiple outputs.

A first, simplistic approach for the MIMO design consists of designing a number of SISO closed-loop that equals the number of sensor/actuator pairs; considering the example in Fig. 1, $m=p$ closed loop can be designed for each of the actuators, based on one sensor located at the same upstream location $z_{0}$. This approach is called decentralized approach. In this case the number of actuators must equal the number of sensors, unless neglecting some of the elements; more importantly, the cross-couplings are disregarded. The cross-couplings may of course affect both stability and performance of the closed loop system. In general, the stability of the closed-loop is not guaranteed. However, if the decentralized controller is stable in each SISO loop and weak couplings characterize the original plant, then the closed loop is also stable. In terms of performance, the resulting controller will be sub-optimal.

The opposite approach consists of accounting all the possible couplings. This methodology is called centralized control design. This choice guarantees optimal performances and stability of the closed loop, although it is less easy to implement in practical situations.

In the following, we discuss only these two limits. In general, linking the sensor signals with the controller that have the strongest interactions allow for hybrid solutions where a number of sensors are wired with a number of actuators; this is the so-called pairing problem.

3.2 Centralized controllers using ADA

Here, we analyse the centralized and decentralized design with respect to the optimal problem. The two cases can be summarized as follows

Decentralized control: $m$ control gains of dimension $n$ are designed independently from each other. Thus, $m$ Riccati equations are solved. 2. 2.

Centralized control: one control gain of dimensions $m\times n$ is designed, by solving one Riccati equation.

The application of the ADA algorithm for the decentralized case is straightforward and it is the one de facto discussed by 2016:bewley:luchini:pralits . Each of the $m$ SISO Riccati equation is replaced by the corresponding ADA iteration, with one input and one output. As already mentioned, this case is only suboptimal and does not guarantee the stability of the closed loop.

The ADA algorithm does not allow identifying a low-rank approximation for the matrix $\mathbf{X}$; thus, the design of the centralized controller needs to be performed row-by-row, using a number of iteration-loops that equals the number of rows. The design of a control kernel of dimensions $\boldsymbol{K}\in\mathbb{R}^{m\times n}$ requires $m$ iteration loops. The cross-coupling are accounted for by choosing the actuators and the sensors to be coupled in the input system and optimization cost function, respectively,

[TABLE]

where we indicate with $\boldsymbol{B}_{\boldsymbol{\textsf{u}},m}$ the matrix including all the $m$ actuators and $\boldsymbol{C}_{\boldsymbol{\textsf{z}},p}$ all the sensors of the original plant. The initial condition of each of the iterations is chosen taking $\tilde{\boldsymbol{q}_{0}}_{i}^{H}=-\boldsymbol{R}_{i,i}^{-1}\boldsymbol{B}_{i}^{H}$, such that each of the $i$-th rows of the matrix $\boldsymbol{K}$ is obtained as adjoint solution of the iteration process.

This design strategy guarantees the coupling among the multiple inputs and corresponds to the solution of a centralized problem. In principle, following the pairing problem, one can select the pairs of sensors/actuators that share the strongest interactions a-priori. In the limit where only one actuator is considered in each of the $m$ loops, we recover the decentralized case where each of the actuator is designed independently from each other.

3.3 Estimation problem

The technique can be extended to any problem based on the solution of algebraic Riccati equations. A first example is provided by the estimation problem as shown in semeraro2013riccati for the single-input-single-output setting. An estimator is defined as the following dynamical system

[TABLE]

The estimator allows to reconstruct the original state $\boldsymbol{q}$, starting from local measurements y; the state $\hat{\boldsymbol{q}}$ is defined such that the error $\boldsymbol{\textsf{y}}-\hat{\boldsymbol{\textsf{y}}}$ is minimized. The estimator is driven by the error term via the matrix $\boldsymbol{L}$, referred to as estimation gain. The estimation gain $\boldsymbol{L}$ is the unknown of the associated estimation problem, obtained as solution of the following Riccati equation

[TABLE]

The $\mathbf{G}$ contains in each of the diagonal entries the estimation penalty, similarly to the control case. From the mathematical point of view, a fictitious adjoint problem is cast such that the cost function

[TABLE]

is mimized. In this dual system, the output $\boldsymbol{C}_{\boldsymbol{\textsf{y}}}^{H}$ is an input of the system, while the adjoint of the inputs $\boldsymbol{B}^{H}_{\mathsf{\textsf{d}}}$ and $\boldsymbol{B}_{\boldsymbol{\textsf{u}}}^{H}$ play the role of the outputs bagheri2009input . The analogy is completed by observing that the “feedback law” is now represented by $\tilde{\boldsymbol{\textsf{y}}}(t)=\boldsymbol{L}^{H}\boldsymbol{p}(t)$. Thus, the first term of the cost function expresses the energy of the dual state $\boldsymbol{p}$, while the second term minimizes the energy input of the feedback law.

A justification of this deterministic approach is given by arfm2007-kim-bewley . The physical interpretation of the method is given within the stochastic framework amr2009-bagheri-et-al . In this case, the unknown of the Riccati equation (16) is the expected energy of the estimation error $\boldsymbol{e}(t)$, while $\mathbf{W}$ is the covariance of the forcing.

The optimization machinery employed for the control problem can be adopted in an analogous manner for computing the full-dimensional estimation gain $\boldsymbol{L}$. The full formulation is given in semeraro2013riccati , where the algorithm is referred to as Adjoint of the Adjoint-Direct (AAD). Also the AAD algorithm can be generalized to the centralized/decentralized case. Thus, in presence of $p$ estimation sensors $\boldsymbol{C}_{\boldsymbol{\textsf{y}}}$, $p$ estimators can be designed based only on one sensor, independently from each other, or a centralized gain can be designed by running an iteration for each of the column of $\boldsymbol{L}$, based on all the sensors $\boldsymbol{C}_{\boldsymbol{\textsf{y}}}$.

3.4 Riccati solution for $\mathcal{H}_{\infty}$ problem

The optimal control framework can be extended to the solution of robust optimization problems, introducing the worst disturbance scenario. The cost function associated with the problem is the following

[TABLE]

where the extra-term to maximize is the unknown, worst disturbance d. The optimization process consists of a simultaneous optimization problem: the signal u is computed such that the a worst-case disturbance is minimized (min-max optimization). The parameter $\gamma$ is chosen by the user. In this sense, the resulting control is sub-optimal. For values of $\gamma\rightarrow\infty$, we approach the limit of the full-order LQR problem already considered. The direct equation of the problem reads

[TABLE]

By introducing the corresponding augmented Lagrangian, it is possible to define the adjoint equation

[TABLE]

Two optimal conditions are defined, by zeroing the respective gradients. The first condition corresponds to the optimal control signal u,

[TABLE]

the second condition is the worst disturbance d

[TABLE]

An ARE can be obtained also for this robust, optimal control problem by imposing the equivalence $\boldsymbol{p}=\boldsymbol{X}\boldsymbol{q}$, into the Hamiltonian system

[TABLE]

The solution is obtained by solving the following problem, based on a ARE equation

[TABLE]

The control gain $\boldsymbol{K}$ and the worst disturbance $\boldsymbol{Y}$ are computed based on the solution $\boldsymbol{X}$ of the ARE equation.

The ADA algorithm can be applied for the full-order approximation of this problem. Considering the single-input-single-output setting, two iterations are necessary for computing, as adjoint solution at final time $T$, $\boldsymbol{K}$ and $\boldsymbol{Y}$ introducing as initial conditions

[TABLE]

Indeed, the robust problem can be seen as a generalized multi-input Riccati equation to be solved, where all the inputs are coupled together. From the algorithm point of view, a relevant difference between the standard $\mathcal{H}_{2}$ problem and the $\mathcal{H}_{\infty}$ is related to the optimization process: due to the simultaneous minimization/maximization underlying the $\mathcal{H}_{\infty}$, we do not seek for a minimum but for a saddle point of the objective function (see for instance gumussoy2009multiobjective ). Finally, note that the both the disturbances $\boldsymbol{B}_{\mathsf{\textsf{d}}}$ and the $\boldsymbol{B}_{\boldsymbol{\textsf{u}}}$ can be multiple; also in this case, a centralized approach needs to be used.

4 Comparisons among algorithms

This section provides numerical results for assessing the performance of the ADA algorithm in presence of multiple inputs and centralized control. We consider a toy-problem and compare ADA with a recent algorithm proposed by maartensson2012scalable . The test-bed is mainly meant at verifying the speed of the convergence for the different algorithms.

4.1 Full-order controllers using stochastic gradients

The algorithm introduced by maartensson2012scalable , and here indicated as MR, is used for comparison. The algorithm computes the control gain by updating at each step the solution as

[TABLE]

where $\nabla_{\boldsymbol{K}}\mathcal{\tilde{J}}$ is the gradient for the update at each step of the iteration. The basic formulation is proposed in appendix B. In principle, the algorithm suffers of the same limitations of the DA iteration: the identified solution depends by the initial condition. The problem is circumvented by choosing a new initial condition randomly at each iteration. The methods can be interpreted as a stochastic gradient descent method; at each step of the iteration, an initial state is chosen for which the gradient is computed. Using different initial states at every iteration, the final feedback matrix will be not depend from any of the initial conditions.

4.2 Implementation

The implemented routines are in prototypical form. The main idea is to analyze the scaling of the computational costs as a function of the number of actuators $m$ and compare the resulting trends. Following this rational, a simple steepest descent algorithm, where the step along the gradient direction is computed as

[TABLE]

The variable $d\mathcal{J}^{i+1}$ indicates the gradient at the current step of the iteration. In principle, acceleration techniques might be applied (for instance, implicit restarting or Nesterov’s accelerated methods) for further improving the performance of the algorithms. Two stop-criteria are chosen: the difference of the control cost between two iterations and the norm of the gradient. The iterations stop when one of the two criteria is fulfilled below a chosen tolerance.

In table 1, we highlight the main features of the algorithms. For all the cases the centralized/decentralized synthesis needs to be imposed when choosing the actuators of each single-input sub-system. In this section we only consider the centralized version of the algorithm that is the most expensive from the computational point of view. Note that for the MR algorithm the basic form is analysed, so without considering a pre-determined structure for the kernel; due to the random selection of the initial condition, the assessment of the performance is done by considering $10$ simulations for each of the analysed cases.

4.3 Numerical example: a distributed system

A toy-problem mimicking a distributed system of dimension $n=50$ is considered. The system matrix $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ is in tri-diagonal form

[TABLE]

The matrix $\boldsymbol{B}\in\mathbb{R}^{n\times m}$ contains non-null, unitary entries only along the main diagonal. The chosen system resembles the one used in maartensson2012scalable . The number of actuators is indicated by $m$; 10 cases are analysed, where the number of actuators is chosen as $m=[5,10,15,\dots,50]$. For the last case, $n=m$. The system is marched in time until the final time $T=20$, with $n_{t}=501$. The value of the final time was deemed sufficient for properly approximating the solution.

The solution of the optimal problem is computed for three different penalty kernels, namely $l=[25,50,100]$ and benchmarked against the Riccati solution, obtained in Matlab using care.m.

In Fig. 2, the number of iterations necessary for convergence is shown as a function of the number of actuators $m$. The chosen tolerance, given as the difference between the cost at $i$ and $i-1$, is $\varepsilon=1.0\times 10^{-6}$. All the algorithms scale with the number of actuators $m$. In particular, it is possible to observe a remarkable regularity for the ADA algorithm, that scales linearly with $m$; different control penalties require different number of iteration. In particular, decreasing $l$ requires a greater number of iterations. The behaviour was already observed in ADA, see semeraro2013riccati . This same behaviour is observed in the DA algorithm; in fact, only at $l=100$ is observed a constant number of iterations with $m$, while a larger number of iterations is required with the increasing number of actuators for the solution of the problem at lower $l$.

The DA algorithm requires on average more iterations when compared to the other iterative methods. This is not surprising, as $n$ loops are required. However, the main bottleneck of DA is the final solution of a linear system, together with the storage of a large matrix containing along the columns the $n$ initial conditions; in that sense, the use of the DA algorithm is unfeasible for larger system and is considered here only for the sake of completeness of our discussion.

The MR algorithm has a less clear behaviour, due to the choice of the initial conditions; an increase of the number of iterations is observed with $m$ also for this case. The choice of the penalty does not affect the number of iterations required. This is clearer when observing Fig. 3, where the iteration number as a function of $m$ is shown for two different penalty gains.

In conclusion, we can summarize the performance of three algorithms as follows

• •

For all the cases, an increase of the total number of iterations with $m$ is observed. This is expected in ADA, where $m$ iteration-loop are required, while it is less obvious for the MR case, where only one iteration-loop is required regardless of $m$.

• •

ADA scales linearly with $m$. This is not the case for the MR algorithm, due to the stochastic gradient application.

• •

For the MR algorithm the averaged number of iterations is not affected by the control effort applied (i.e., the choice of the control penalty $l$).

5 Numerical example: the two-dimensional Kuramoto-Sivashinsky

In this section, a modified version of the Kuramoto-Sivashinsky equation is used for testing the multi-input case, mimicking the three-dimensional setup sketched in Fig. 1; this proof-test provides the design feasibility in larger computational domains.

The flow field is approximated in a two-dimensional domain $x$-$z$ plane, and governed by the equation

[TABLE]

With respect to the original equation, the dynamics is linearized around the convective velocity, $V=0.4$. The non-dimensional terms $R$ and $P$ are defined as

[TABLE]

with $\eta$ representing an energy production term and $\mu$ the dissipation; $L$ is the reference length of the system. The parameters can be chosen ad-hoc such that the dispersion relation of the system mimics the evolution of travelling packets of Tollmien-Schlichting (TS) waves. Introducing the maximum streamwise and spanwise wavenumbers as $\alpha_{max},\beta_{max}\in\mathbb{R}$ respectively, and the complex temporal frequency $\omega$, the non-dimensional numbers can be expressed as

[TABLE]

with $S$ providing the modulation along the spanwise direction of the wave. The numbers are chosen such that the dispersion relation closely resembles the behaviour in the horizontal plane of a TS wave evolving on a flat-plate at $Re_{\delta^{*}}=1000$ (see semeraro2011feedback , semeraro2013transition ). More details are reported by johan2015controlling .

The computational box extends along the streamwise direction in the interval $x\in[0,500]$ and $z\in[-90,90]$ along the spanwise direction; a grid with $N_{x}=256$ and $N_{z}=96$ is chosen. The spatial discretization is performed by means of a pseudo-spectral method, that includes a fringe region extending between $x=400$ and $x=500$; periodic boundary conditions are imposed along the spanwise direction. The time marching is performed using a three-steps Runge-Kutta scheme; the basic implementation can be found in the repository https://github.com/nfabbiane/ks2D.

5.1 A multi-input full-order controller

The resulting dynamics is depicted in Fig. 4$a$-$c$-$e$, where the impulse response of the system is shown at three different instants. The disturbance is introduced at $\mathbf{x}_{0}=\left(x_{0},z_{0}\right)=\left(2.5,0\right)$, and its spatial distribution is modelled as

[TABLE]

with $\sigma=4$. The dynamics of the wavepacket mimics the evolution of a TS wave, growing as it propagates downstream along the streamwise direction while extending along the spanwise direction with a backward bending. The corresponding energy growth associated with the evolving wavepacket is shown in Fig. 5$b$ (black, solid line); the root-mean-square energy is obtained from the statistics associated with a stochastically driven simulation, averaged over a span of $T=10000$ time units.

A LQR controller is designed using the ADA algorithm. The controller is centralized: all the gains are designed such that the set of sensors and actuators is coupled. The setup closely resembles the one analysed by semeraro2011feedback , semeraro2013transition and reproduced schematically in Fig. 1: a row composed by $n=9$ localized actuators ($\boldsymbol{B}_{\boldsymbol{\textsf{u}}}$) is placed along the spanwise direction at $x=200$, equi-spaced of $\Delta z=20$. The same setup along the spanwise direction is chosen for the sensors $\boldsymbol{C}_{z}$, placed at $x=300$. All the chosen elements are modelled as Gaussian distributions, Eq. 30. Nine iterations are set, one for each control gain. An example is given in Fig. 5$a$, where the control gain placed at $\mathbf{x}=\left(200,0\right)$, with $l=100$; the controller shows the typical signature of the adjoint solution, as it is bent backward with respect to the propagation of the direct solution. The performances of two controllers are shown in Fig. 5$b$, where the clean configuration is compared with two cases: $l=100$ (red solid line), and $l=500$ (blue solid line). The results are in qualitative agreement with semeraro2011feedback . More interestingly, it shows that the technique enables the computation of multi-input, full order controllers without using a preliminary model reduction step.

6 Conclusions

We extended the framework of the Adjoint of the adjoint algorithm to multivariate large systems, by highlighting how decentralized, centralized and coupled controllers can be computed for full-order, optimal control. These observations apply straightforwardly to the dual problem, the estimation problem, and to robust controllers belonging to the $\mathcal{H}_{\infty}$ framework. The main advantage of the ADA algorithm (and its extensions) is the independence of the final solution from the initial conditions. From the physical point of view, this implies that it is not required any knowledge of the disturbances active in the flow. This feature makes the algorithm particularly appealing for the pre-assessment of the optimal performance of a controller, based on linear (or linearized) plant. Moreover, one can assess the controllability of the dynamical system by such techniques and estimate the efficiency of the resulting controller.

Two test cases are considered. A toy-problem mimicking a distributed system is analysed for assessing the convergence performance of ADA against an analogous algorithm by Mårtensson et al. maartensson2012scalable , indicated as MR. The results show a remarkable scalability of ADA, with respect the number of actuators; with respect to the MR algorithm, the number of iterations required is on average three times smaller and - more importantly - it is not affected by the stochastic approach required by the MR algorithm for identifying the global solution. The second example is based on a modified version of the Kuramoto-Sivashinsky. This example reproduces a control setup based on multiple, localized sensors and actuators placed along the cross-flow direction. The feasibility of the approach is demonstrated suggesting that the technique might applied also to full three-dimensional cases in more complex fluid mechanics settings.

Future work will be devoted to the application of these methods in combination with sparsity promoting algorithms (see 2013:lin:fard:jova ). In our application, the coupling among the actuators is pre-determined and not optimized. However, it is possible to introduce the sensor/actuator pairing within the optimization process (pairing problem). Moreover, sensor and actuator placement for large-scale systems might be analysed starting from this application (chen2011h , N. Fabbiane, private communication, 2016).

Acknowledgments Anders Rantzer and his collaborators (Lund University) are acknoweldged for interesting discussions and for bringing up their work on optimal control synthesis in large systems. Nicolò Fabbiane (DAFE, Onera) provided the basic implementation of the two-dimensional Kuramoto-Sivashinsky equation. The first author thanks the support by the Agence Nationale de la Recherche (ANR) under the CoolJazz project, grant number ANR-$12$-BS$09$-$0024$.

Appendix A Derivation of ADA algorithm using integration by parts

In this section we briefly summarize the original derivation of the ADA algorithm, as first proposed by 2016:bewley:luchini:pralits . We start from the direct-adjoint system

[TABLE]

where the optimality condition has already been included. By introducing two variables, namely the state vectors $\boldsymbol{a}(t)$ and $\boldsymbol{b}(t)$, we can recombine the system as

[TABLE]

Integration by parts allows rewriting the relation as

[TABLE]

Gathering the resulting state equations for $\boldsymbol{a}(t)$ and $\boldsymbol{b}(t)$, and changing the variables $\boldsymbol{a}\rightarrow\tilde{\boldsymbol{p}}$, $\boldsymbol{b}\rightarrow-\tilde{\boldsymbol{q}}$, allow us to write the following system

[TABLE]

Equations 32 correspond to the ones of the original optimal control problem 31. The last relation 32$(c)$ is exactly the condition recovered by symplectic product in Eq. (11); thus, the observations done in Eqs. (11–12) are valid also for this case.

Appendix B Full-order controllers using stochastic gradients

In this appendix, the algorithm discussed by Mårtensson and Rantzer in maartensson2012scalable is derived following the Lagrangian approach. For sake of conciseness, we will indicate the algorithm with the acronym MR. The state equation is re-written as

[TABLE]

The problem is rewritten now assuming the control kernel $\boldsymbol{K}$ as unknown of the problem, replacing $\boldsymbol{\textsf{u}}(t)$. The cost function is now

[TABLE]

and the resulting augmented Lagrangian is

[TABLE]

The resulting system is obtained by considering the gradients with respect to $\boldsymbol{p}$, $\boldsymbol{q}$ and $\boldsymbol{K}$

[TABLE]

Zeroing the gradients, the system can be arranged in matrix form as

[TABLE]

The solution is updated by using $\ref{eq:sysK}c$ as

[TABLE]

The iteration is usually initialized with a guess for $\boldsymbol{K}$, unless the system under consideration is asymptotically stable; in the latter case, the controller is designed for modifying the transient dynamics and the starting guess is simply a null vector. In this formulation, the problem depends on the initial conditions. The solution is not guaranteed to be the one corresponding to the global minimum for a given initial condition. This drawback can be circumvented by using ad-hoc strategies as the stochastic gradient approach described in maartensson2012scalable . Moreover, with respect to the standard algorithm, the MR algorithm does not require the solution of Lyapunov equations, replaced with an approximation based on the adjoint and direct solutions.

B.1 Relation with the ADA algorithm: generalizing the sensitivity approach

The symplectic product adopted in ADA can be introduced also for the system in Eq. B. It is possible to show that the relation

[TABLE]

is fulfilled also in this case. Considering again the gradient,

[TABLE]

it is null $\forall\boldsymbol{q}$ if

[TABLE]

By rescaling the relation via $\boldsymbol{R}^{-1}$, we recover the equivalence between the product and the optimality condition at $t=0$

[TABLE]

Thus, also for the basic equation of the MR algorithm, introducing as initial conditions of the system $\boldsymbol{\tilde{q}}=-\boldsymbol{R}_{i,i}^{-1}\boldsymbol{B}_{i}$, the iteration enables to compute the $i$-th kernel $\boldsymbol{K}$ as solution of the adjoint equation at the final time $T$. All the observations done for the multi-input, coupled version of the ADA algorithm are valid also for this algorithm. Interestingly, this algorithm produces the control kernels in two different ways: i) as adjoint solution; ii) as result of the gradient-based iteration

[TABLE]

The sensitivity with respect to the initial conditions applied to the optimal control problem leads to the same results with different choice of gradients. The gradient with respect to $\boldsymbol{K}$ can be rewritten, using the chain-rule, as

[TABLE]

The underlying system of equations is equivalent in the two cases.