# Data-driven computation of adjoint sensitivities without adjoint solvers: An
application to thermoacoustics

Defne E. Ozan  Imperial College London, Department of Aeronautics, Exhibition
Road, London SW7 2BX, UK  Luca Magri  Imperial College London, Department of
Aeronautics, Exhibition Road, London SW7 2BX, UK  The Alan Turing Institute,
London NW1 2DB, UK  Politecnico di Torino, DIMEAS, Corso Duca degli Abruzzi,
24 10129 Torino, Italy

###### Abstract

Adjoint methods have been the pillar of gradient-based optimization for
decades. They enable the accurate computation of a gradient (sensitivity) of a
quantity of interest with respect to all system’s parameters in one
calculation. When the gradient is embedded in an optimization routine, the
quantity of interest can be optimized for the system to have the desired
behaviour. Adjoint methods, however, require the system’s Jacobian, whose
computation can be cumbersome, and is problem dependent. In this paper, we
propose a computational strategy to infer the adjoint sensitivities from data
(observables), which bypasses the need of the Jacobian of the physical system.
The key component of this strategy is an echo state network, which learns the
dynamics of nonlinear regimes with varying parameters, and it evolves
dynamically via a hidden state. Although the framework does not make
assumptions on the dynamical system, we focus on thermoacoustics, which are
governed by nonlinear and time-delayed systems. First, we show that a
parameter-aware Echo State Network (ESN) infers the parameterized dynamics.
Second, we derive the adjoint of the ESN to compute the sensitivity of time-
averaged cost functionals. Third, we propose the Thermoacoustic Echo State
Network (T-ESN), which hard constrains the physical knowledge in the network
architecture. Fourth, we apply the framework to a variety of nonlinear
thermoacoustic regimes of a prototypical system. We show that the T-ESN
accurately infers the correct adjoint sensitivities of the time-average
acoustic energy with respect to the flame parameters. The results are robust
to noisy data, from periodic, through quasiperiodic, to chaotic regimes. The
inferred adjoint sensitivities are employed to suppress an instability via
steepest descent. We show that a single network predicts the nonlinear
bifurcations on unseen scenarios, which allows it to converge to the minimum
of the acoustic energy. This work opens new possibilities for gradient-based
data-driven design optimization.

##  I Introduction

Engineers and scientists strive to design flow and acoustic systems that
behave in an optimal way [e.g., 1, 2, 3, 4, 5]. Finding the global optimum is
challenging because the temporal dynamics can be nonlinear, time-delayed, and
chaotic, which makes the objective function non-convex [e.g, 6, 7]. A
practical approach is to solve the optimization problem by gradient-descent,
in which the objective function is optimized along the gradient direction
until a local optimum is found. The gradient (also known as sensitivity)
offers the information about how the quantity of interest changes with the
system’s parameters. Central to gradient-based optimization methods are
adjoint methods. Adjoint methods enable the accurate computation of the
gradient with a computational cost that is not affected by the number of
design parameters. For example, if the system has $N$ parameters, adjoint
methods compute the gradient of a objective function with only one
calculation. This is is contrast with tangent linear or finite difference
methods, which would require $N$ computations [e.g., 8, 9, 10, 2, 3, 5].

The optimization of a quantity of interest of a nonlinear system with adjoint
methods requires, in general, six ingredients: (i) The system’s nonlinear
governing equations are integrated forward in time to obtain the baseline
(direct) solution; (ii) the nonlinear equations are linearized to obtain the
Jacobian of the system; (iii) the adjoint equations are derived from the
Jacobian, (iv) the adjoint equations are initialized at the end of the direct
solution, and are integrated backwards in time (because in nonlinear systems
the Jacobian depends on the direct solution, strategies such as checkpointing
become often necessary to decrease the memory requirements [11]); (v) at the
end of the adjoint integration, the gradient of the objective function to all
the design parameters of interest can be evaluated by an inner product; and
(iv) the gradient is embedded in a gradient-based optimization scheme to
iteratively compute the optimum of the objective function. This process is
also known as the direct-adjoint looping. On the one hand, in finite
difference or tangent linear approaches, the computational cost increases
linearly with the number of parameters because the sensitivity calculation
needs to be performed for each parameter. On the other hand, the adjoint
method computes the gradient of the objective with respect to all parameters
with a single operation, which makes it computationally cheaper when the
number of parameters is larger than the number of objective functionals [e.g.,
12, 13, 14, 15, 8]. For this reason, adjoint methods have found extensive
application in design optimization, and both linear and nonlinear sensitivity
analysis of fluid dynamics applications.

Although adjoint methods are excellent tools, they are problem dependent,
i.e., they depend on the system’s equations. First, the adjoint sensitivity is
only as accurate as the model assumptions. Low-order models are
computationally cheap to evaluate, however the uncertainties in the model can
lead to errors in the sensitivity computations, especially in nonlinear
systems that are sensitive to perturbations. Second, the system must be
linearized, i.e., by computing the Jacobian. This means the computation of the
adjoint solver needs to be performed for different fluid problems, numerical
schemes, and systems. Automatic differentiation [16, 17] and surrogate models
[18] help the automation of the process, but they still require a code that
may be problem dependent.

To bypass the need of system’s equations, a data-driven model can be employed.
The key idea of this paper is to learn an approximation of the dynamical
systems from observables, and then compute the adjoint of it. Recurrent neural
networks are nonlinear dynamical systems, which infer the governing dynamics
from observables (data). We employ a type of reservoir computer, the Echo
State Network (ESN) [19], which is a universal approximator of the dynamics
[20]. The ESN benefits from straightforward training because it does not
require backpropagation. Its Jacobian can also be computed in a
straightforward way due to the shallow architecture [21]. ESNs have been shown
to successfully (i) make time-accurate short-term forecasts [22, 23, 24], (ii)
reconstruct unmeasured variables [25, 26], (iii) learn ergodic properties,
such as Lyapunov exponents, and long-term statistics of chaotic flows [27, 7,
21], and (iv) infer model errors in thermoacoustics [28]. In recent advances,
low-order representations of turbulence were learned from high-dimensional
data, and the dynamics of this low-dimensional system was propagated with the
ESN [29]. In this paper, we take an additional step. Because in order to
perform design optimization of system’s parameters, we require a parameter-
aware extension of the Echo State Network [30], which learns the parameterized
dynamics and can operate across parameter regimes.

We apply the data-driven sensitivity framework to nonlinear thermoacoustic
oscillations. These are phenomena that can be detrimental in the manufacturing
of gas turbines [31], and rocket engines [32]. In thermoacoustic systems, if
the heat release and acoustic pressure are sufficiently in phase, self-
sustained oscillations can occur [33, 34]. Therefore, their prediction,
control, and suppression are important both during design and during
operation. This is especially challenging because thermoacoustics exhibit
nonlinear behaviour and bifurcations to periodic, quasiperiodic, and chaotic
regimes [35], which motivates the development of methods for nonlinear
analysis. Preliminary design of thermoacoustic systems is performed by linear
stability [31, 36]. Linear stability analysis determines the growth rates of
infinitesimal perturbations to the mean flow, i.e., the eigenvalues of the
linearized system of equations around the mean flow [37]. To name a few
examples, adjoint-based linear stability analysis enabled calculation of
sensitivity to parameters and passive control mechanisms [38], optimal
placement and tuning of acoustic dampers [39] and stabilization of
thermoacoustic oscillations by geometry changes [40]. An in-depth review of
development of adjoint methods for thermoacoustics can be found in [5, 37].
Beyond linear stability analysis, adjoint methods were used for weakly
nonlinear analysis to predict subcritical Hopf bifurcations [41], and Floquet
analysis to calculate the stability and the period of limit cycle oscillations
[5]. The stability of chaotic thermoacoustics was studied by Huhn and Magri
[6] using covariant Lyapunov analysis, who analytically showed that this is
equivalent to eigenvalue and Floquet analysis for fixed point and limit cycle
solutions, respectively. As a result, Huhn and Magri [6] achieved gradient-
based minimization of acoustic energy in chaotic thermoacoustics using
shadowing methods. The application of this method was challenged by the
discontinuities observed in the acoustic energy when bifurcations occur,
computational cost of shadowing methods, and non-hyperbolicity of
thermoacoustics for some regimes. To circumvent these issues, a gradient-free
optimization method based on Bayesian optimization was proposed [7]. Among
applications of data-driven methods in thermoacoustics, the ESN was previously
used to predict the nonlinear dynamics [7] and calculate the Lyapunov
exponents [42] of a chaotic regime, and infer the model error of a low-order
model [28].

The overarching objective of this paper is to address the challenge of
developing code-specific adjoint solvers, which requires the system’s
equations and the Jacobian to be known. We propose data-driven nonlinear
direct-adjoint looping to compute parametric sensitivities, and apply this to
a prototypical thermoacoustic system. The specific objectives are to build a
robust parameter-aware ESN architecture for this purpose, derive its adjoint,
and demonstrate the data-driven inference of parameterized dynamics and
associated adjoint sensitivities. We refer to this architecture as the
thermoacoustic ESN and show that (i) it learns the parameterized dynamics of
nonlinear, time-delayed thermoacoustics, (ii) predicts bifurcations, (iii)
accurately obtains the adjoint sensitivities, (iv) the acoustic energy is
minimized in a gradient-based optimization framework, and (v) demanding
scenarios such as training with noisy or chaotic data can be handled.

This paper is organized as follows. In Sec. II, we give a background on the
mathematical framework for sensitivity analysis and the parameter-aware ESN.
In Sec. III, we derive the adjoint of the ESN. In Sec. IV, we present the
nonlinear and time-delayed thermoacoustic model, namely the Rijke tube. In
Sec. V, we introduce the time-delayed ESN designed for thermoacoustics and
derive its adjoint. In Sec. VI, we show the performance of the ESN in terms of
prediction and accuracy of the learned adjoint sensitivities, and employ it
for the gradient-based design optimization of thermoacoustics in Sec. VII. We
address training with noisy data in Sec. VIII, and learning chaotic regimes in
Sec. IX. Section X concludes the paper.

##  II Background

###  II.1 Sensitivity Analysis

We consider a dynamical system in the form of

$\dot{\mathbf{x}}-\mathbf{F}(\mathbf{x},\mathbf{p})=0,$ (1)

where $\mathbf{x}\in\mathbb{R}^{n_{x}}$ is the state vector,
$\dot{\mathbf{x}}$ is its time-derivative, $\mathbf{p}\in\mathbb{R}^{n_{p}}$
is the parameters’ vector, and
$\mathbf{F}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{p}}\rightarrow\mathbb{R}^{n_{x}}$
is the operator that describes the system’s evolution. We assume that
$\mathbf{F}$ is smooth so that the Jacobian of the system, i.e.,
$\partial\mathbf{F}/\partial\mathbf{x}$, is defined everywhere, which is a
requirement for gradient-based sensitivity analysis. We wish to optimize an
objective functional, $\mathcal{J}$, which is a time-average of
$\tilde{\mathcal{J}}$ [5]

$\mathcal{J}=\frac{1}{T}\int_{0}^{T}\tilde{\mathcal{J}}(\mathbf{x}(t),\mathbf{p})dt,$
(2)

subjected to the system dynamics (1). A straightforward approach to obtain
sensitivity information is to apply a small perturbation, $h$, to the
parameters of interest, $p_{i}\in\mathbf{p}$ and numerically compute the
sensitivity using finite differences

$\frac{d\mathcal{J}}{d\mathbf{p}}\bigg{\rvert}_{\mathbf{p}=\bar{{\mathbf{p}}}}=\frac{d\mathcal{J}}{dp_{i}}\bigg{\rvert}_{p_{i}=\bar{{p}}_{i}}\approx\frac{\mathcal{J}\big{\rvert}_{p_{i}=\bar{p}_{i}+h}-\mathcal{J}\big{\rvert}_{p_{i}=\bar{p}_{i}}}{h}+O(h),\quad
i=1,2,\dots,n_{p},$ (3)

which needs to be repeated for each $p_{i}\in\mathbf{p}$ to determine
$d\mathcal{J}/d\mathbf{p}$. This means that the computational cost of the
sensitivity increases linearly with the number of parameters. A more efficient
approach is to use the adjoint equation. Analytically, we can express the
gradient $d\mathcal{J}/d\mathbf{p}$ as

$\displaystyle\frac{d\mathcal{J}}{d\mathbf{p}}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\left(\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{p}}+\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{x}}\mathbf{Q}\right)dt,$
(4a) $\displaystyle\dot{\mathbf{Q}}$
$\displaystyle=\frac{\partial\mathbf{F}}{\partial\mathbf{p}}+\frac{\partial\mathbf{F}}{\partial\mathbf{x}}\mathbf{Q},$
(4b) $\displaystyle\mathbf{Q}(0)$ $\displaystyle=\mathbf{0},$ (4c)

where we define $\mathbf{Q}\equiv
d\mathbf{x}/d\mathbf{p},\;\mathbf{Q}\in\mathbb{R}^{n_{x}\times n_{p}}$. We
refer to solving the system of equations (4) as the tangent linear problem. In
this case, the sensitivity can be computed by formulating the tangent linear
system around a base trajectory and solving for $\mathbf{Q}$ in time. Whilst
this approach eliminates the finite difference numerical error, it still
suffers from a computational cost that is linear with the number of parameters
as the dimension of $\mathbf{Q}$ grows with the dimension of $\mathbf{p}$.
Adjoint methods circumvent this issue. By introducing the Lagrangian of the
objective (2) and the adjoint variables $\mathbf{q}^{+}\in\mathbb{R}^{n_{x}}$
as the Lagrange multipliers of the constraint (1), the gradient
$d\mathcal{J}/d\mathbf{p}$ can be rewritten as

$\displaystyle\frac{d\mathcal{J}}{d\mathbf{p}}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\left(\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{p}}+\frac{\partial\mathbf{F}}{\partial\mathbf{p}}^{\top}\mathbf{q}^{+}\right)dt,$
(5a) $\displaystyle\dot{\mathbf{q}}^{+}$
$\displaystyle=-\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{x}}^{\top}-\frac{\partial\mathbf{F}}{\partial\mathbf{x}}^{\top}\mathbf{q},$
(5b) $\displaystyle\mathbf{q}^{+}(T)$ $\displaystyle=0,$ (5c)

where the adjoint variables $\mathbf{q}^{+}$ need to solved by integrating the
adjoint system backwards in time. The details of the derivation of the adjoint
problem from a Lagrangian perspective for a continuous-time system can be
found in [5]. In this paper, we derive the tangent linear and adjoint
equations for the proposed data-driven algorithm in Appendix A.

###  II.2 Parameter-aware Echo State Networks

Given the evolution equations of a dynamical system (1), we can compute its
adjoint (5). Developing the adjoint solver however requires the system
Jacobian, $\partial\mathbf{F}/\partial\mathbf{x}$. This Jacobian is code
specific. The key idea of this paper is to eliminate code-specific adjoint
solvers by having a machine that learns the dynamics from data, and is itself
a dynamical system, the Jacobian of which is readily available independent of
the system in question. Recurrent neural networks (RNNs) are well-suited for
this task, as they are neural network architectures that can process
sequential data and make forecasts of timeseries from dynamical systems.
Although the dynamical systems we consider (1) are continuous in time, RNNs
accept data as a sequence in discrete time. RNNs possess memory through a
hidden state and thus they can learn temporal relationships. The hidden state
enables the RNNs to be framed as discrete dynamical systems [21]. In this
work, we utilize the Echo State Network (ESN) [43], which is a type of RNN
that is straightforward to train because it does not require backpropagation
during training. The ESN is a universal approximator [20] and it follows this
dynamic equation

$\mathbf{r}(i+1)=(1-\alpha)\mathbf{r}(i)+\alpha\tanh(\mathbf{W}_{in}\mathbf{y}_{in}(i)+\mathbf{W}\mathbf{r}(i)),$
(6a) where $\mathbf{y}_{in}\in\mathbb{R}^{n_{y}}$ is the input vector,
$\mathbf{r}(i)\in\mathbb{R}^{n_{r}}$ is the reservoir state,
$\mathbf{W}_{in}\in\mathbb{R}^{n_{r}\times n_{y}}$ is the input matrix, and
$\mathbf{W}\in\mathbb{R}^{n_{r}\times n_{r}}$ is the state matrix. Equation
(6a) maps the reservoir state at time $t_{i}$ to the reservoir state at time
$t_{i+1}=t_{i}+\Delta t$. The prediction of the output is a linear combination
of the reservoir state
$\mathbf{\hat{y}}(i+1)=\mathbf{W}_{out}\mathbf{r}(i+1),$ (6b)

where $\mathbf{y}\in\mathbb{R}^{n_{y}}$ is the target output vector, the
dynamics of which we aim to model, $\mathbf{\hat{y}}\in\mathbb{R}^{n_{y}}$ is
the predicted output vector, and $\mathbf{W}_{out}\in\mathbb{R}^{n_{y}\times
n_{r}}$ is the output matrix that maps the reservoir state to the output
vector. The matrices $\mathbf{W}_{in}$ and $\mathbf{W}$ are sparse, randomly
generated, and not trained, whilst $\mathbf{W}_{out}$ is trained via ridge
regression [19]. The input matrix is scaled as
$\mathbf{W}_{in}=\sigma_{in}\tilde{\mathbf{W}}_{in}$, where
$\tilde{\mathbf{W}}_{in}$ is constructed such that each row contains only one
nonzero element drawn from a uniform distribution $\sim\mathcal{U}(-1,1)$,
i.e., the inputs are activated individually. The state matrix $\mathbf{W}$
describes the reservoir connections, and is scaled as
$\mathbf{W}=\rho\tilde{\mathbf{W}}$, where $\tilde{\mathbf{W}}$ is an Erdős-
Renyi matrix with a given average number of nonzero elements per row
(connectivity) drawn from a uniform distribution $\sim\mathcal{U}(-1,1)$. The
scalar $\rho$ is the spectral radius of $\mathbf{W}$. In this architecture,
the reservoir state variables “echo” the input, and so the reservoir holds
memory. During training $\mathbf{y}_{in}$ is fed into the network in open-loop
and the evolution of the reservoir state is saved. The reservoir state is
initialized to a zero vector. An initial time window is used for the washout
stage, which eliminates the effects of the initial condition and is discarded
in training. The training is performed by solving for $\mathbf{W}_{out}$ [19]

$\mathbf{W}_{out}^{*}=\underset{\mathbf{W}_{out}}{\mathrm{arg\,min}}\;\frac{1}{n_{y}}\sum_{j=1}^{n_{y}}\left(\sum_{i=1}^{N_{t}}(y_{j}(i)-\hat{y}_{j}(i))^{2}+\lambda||\mathbf{w}_{out,j}||^{2}\right),$
(7)

where $N_{t}$ is the number of training steps, $\mathbf{w}_{out,j}$ is the
$j^{th}$ row of $\mathbf{W}_{out}$, and $\lambda$ is the Tikhonov coefficient.
This is a ridge regression problem, which minimizes both the error between
network predictions for the next step and the norm of the output weights for
regularization. The solution is [19]

$\mathbf{W}_{out}^{*}=\mathbf{Y}\mathbf{R}^{\top}(\mathbf{R}\mathbf{R}^{\top}+\lambda\mathbf{I})^{-1},$
(8)

where
$\mathbf{R}=[\mathbf{r}(1)\;\mathbf{r}(2)\;\dots\;\mathbf{r}(N_{t})],\;\mathbf{R}\in\mathbb{R}^{n_{r}\times
N_{t}}$,
$\mathbf{Y}=[\mathbf{y}(1)\;\mathbf{y}(2)\;\dots\;\mathbf{y}(N_{t})],\;\mathbf{Y}\in\mathbb{R}^{n_{y}\times
N_{t}}$, and $\mathbf{I}$ is the identity matrix. Tuning an ESN involves
searching for optimal hyperparameters that consist of input matrix scaling
$\sigma_{in}$, the spectral radius of the state matrix $\rho$ (typically $\leq
1$ to satisfy the echo state property), the leak rate $\alpha$, and the
Tikhonov regularizer $\lambda$. These hyperparameters are determined by
validation. A detailed explanation of training and robust validation of Echo
State Networks can be found in [24]. After training, the ESN can run
autonomously in closed-loop by setting
$\mathbf{y}_{in}(i)=\mathbf{\hat{y}}(i)$.

In the formulation (6a), the ESN is unaware of the system parameters, i.e., it
can only operate in one regime. For our purposes, we require an ESN to learn
the parameter-dependent dynamics of a system such that we can capture the
variations in different regimes. A parameter-aware extension of the ESN has
been employed to predict when the system bifurcates to fixed point solutions,
a phenomenon known as amplitude death [30], and bifurcations in multi-stable
regimes including chaos [44]. It can be formulated as [30]

$\mathbf{r}(i+1)=(1-\alpha)\mathbf{r}(i)+\alpha\tanh(\mathbf{W}_{in}[\mathbf{y}_{in}(i);\mathrm{diag}(\mathbf{\bm{\sigma}}_{p})(\mathbf{p}-\mathbf{k}_{p})]+\mathbf{W}\mathbf{r}(i)),$
(9)

where $(;)$ denotes vertical concatenation, $\mathbf{p}\in\mathbb{R}^{n_{p}}$
denotes the system’s parameters and are fed to the network through the input
channel, which we assume to be constant in time for each system. The input
matrix
$\mathbf{W}_{in}=[\mathbf{W}_{in}^{y}\;\mathbf{W}_{in}^{p}],\;\mathbf{W}_{in}^{y}\in\mathbb{R}^{n_{r}\times
n_{y}},\;\mathbf{W}_{in}^{p}\in\mathbb{R}^{n_{r}\times n_{p}}$ is constructed
such that the parameters are fully connected with all the reservoir nodes,
i.e., the weights associated with the parameters $\mathbf{W}_{in}^{p}$ are
nonzero for all rows. The parameter-aware ESN adds two additional
hyperparameters per physical parameter $p_{i}$. A scalar $k_{p_{i}}$ shifts
the parameter $p_{i}$ as $p_{i}-k_{p_{i}}$, and a scalar $\sigma_{p_{i}}$
scales it as $\sigma_{p_{i}}(p_{i}-k_{p_{i}})$. Hence, we apply the input
scaling $\sigma_{in}$ on $\mathbf{W}_{in}^{y}$ only. The training of the
parameter-aware ESN is performed as the standard ESN (8). The dataset contains
data from different parameter regimes, which means that we need to save the
evolution of the reservoir state for each training regime, and concatenate as
$\mathbf{R}=[\mathbf{R}^{1}\;\mathbf{R}^{2}\;\dots,\mathbf{R}^{N_{tr}}]$
(likewise, the outputs are concatenated as
$\mathbf{Y}=[\mathbf{Y}^{1}\;\mathbf{Y}^{2}\;\dots\mathbf{Y}^{N_{tr}}]$, where
$N_{tr}$ is the number of training regimes.)

##  III Adjoint of the Echo State Network

In this section, we derive the adjoint of the parameter-aware ESN. When the
ESN runs in closed-loop, it acts as an autonomous dynamical system. This is
the configuration we use for the computation of the adjoint. We consider the
time-averaged objective functional (2), which we wish to optimize. In discrete
time, we reformulate it as a sum over a discrete number of time steps, denoted
by $N$

$\mathcal{J}=\frac{1}{N}\sum_{i=1}^{N}\tilde{\mathcal{J}}(\mathbf{r}(i)).$
(10)

Without loss of generality, we assume that the objective functional does not
explicitly depend on the parameters, i.e.,
$\partial\tilde{\mathcal{J}}/\partial\mathbf{p}=0$. For the parameter-aware
ESN, the tangent linear problem is

$\displaystyle\frac{d\mathcal{J}}{d\mathbf{p}}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\frac{d\tilde{\mathcal{J}}(\mathbf{r}(i))}{d\mathbf{r}(i)}\mathbf{Q}(i),$
(11a) $\displaystyle\mathbf{Q}(i+1)$
$\displaystyle=\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{p}}+\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{r}(i)}\mathbf{Q}(i),$
(11b) $\displaystyle\mathbf{Q}(0)$ $\displaystyle=\mathbf{0},$ (11c)

where $\mathbf{Q}(i)\equiv
d\mathbf{r}(i)/d\mathbf{p},\;\mathbf{Q}(i)\in\mathbb{R}^{n_{r}\times n_{p}}$.
Physically, Eq. (11c) means that the initial condition is independent of the
parameter. Because the dimension of the tangent linear system,
$\mathbf{Q}(i)$, grows with the number of parameters, adjoint methods enable a
computationally cheap and accurate way to obtain sensitivities [5]. Using the
Lagrangian formulation of the objective and constraints given by the evolution
equations (Appendix A), the adjoint problem is formulated as

$\displaystyle\frac{d\mathcal{J}}{d\mathbf{p}}$
$\displaystyle=\sum_{i=1}^{N}\mathbf{q}^{+\top}(i)\frac{\partial\mathbf{r}(i)}{\partial\mathbf{p}},$
(12a) $\displaystyle\mathbf{q}^{+}(i)$
$\displaystyle=\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{r}(i))}{\partial\mathbf{r}(i)}^{\top}+\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{r}(i)}^{\top}\mathbf{q}^{+}(i+1),$
(12b) $\displaystyle\mathbf{q}^{+}(N)$
$\displaystyle=\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{r}(N))}{\partial\mathbf{r}(N)}^{\top}.$
(12c)

Equation (12) is in a form similar to the continuous time equations (5),
however, we have replaced the system’s Jacobian with the ESN Jacobian, thus,
the adjoint variables are the adjoint of the ESN’s reservoir state.
Algorithmically, we first obtain the direct solution, $\mathbf{r}(i)$, as a
baseline trajectory by running the ESN in closed-loop for a time window. Then,
the adjoint equations (12b) are solved backwards in time starting from the
terminal condition $\mathbf{q}^{+}(N)$ (12c). This procedure requires the
computation of the Jacobian, $\partial\mathbf{r}(i+1)/\mathbf{r}(i)$, and the
partial derivative with respect to the parameters
$\partial\mathbf{r}(i)/\partial\mathbf{p}$, associated with the ESN. These
gradients are provided in Appendix B.

The objective functional we are interested in depends on the output vector,
e.g.,

$\mathcal{J}=\frac{1}{N}\sum_{i=1}^{N}||\hat{\mathbf{y}}(i)||^{2}_{2}=\frac{1}{N}\sum_{i=1}^{N}||\mathbf{W}_{out}\mathbf{r}(i)||^{2}_{2}.$
(13)

Thus, the gradient of the objective with respect to the reservoir state is

$\frac{\partial\tilde{\mathcal{J}}(\mathbf{r}(i))}{\partial\mathbf{r}(i)}=2\mathbf{W}_{out}^{\top}\mathbf{W}_{out}\mathbf{r}(i).$
(14)

Minimization of equation (13) defines the optimization problem where the ESN
models the physical dynamical system.

##  IV Thermoacoustics

We consider a prototypical thermoacoustic system known as the Rijke tube,
which consists of a duct with a heat source located inside. The Rijke tube is
a setup used to propose fundamental methods of thermoacoustic instabilities
because it captures the qualitative nonlinear dynamics and bifurcations
observed in real-life applications [e.g., 45, 35, 6]. In this model, we make
the following assumptions; (i) the acoustics are one-dimensional, (ii) the
mean-flow has a low Mach number and the acoustic variables are modelled as
perturbations on top of a mean flow with uniform density across the tube,
(iii) the heat release is compact, and (iv) the boundary conditions are ideal,
i.e., fully reflective. The governing partial differential equations (PDEs)
are derived from conservation of momentum, energy, and mass [e.g., 5]

$\displaystyle\frac{\partial u}{\partial t}+\frac{\partial p}{\partial x}=0,$
(15a) $\displaystyle\frac{\partial p}{\partial t}+\frac{\partial u}{\partial
x}+\zeta p-\dot{q}\delta(x-x_{f})=0,$ (15b)

where $u$ and $p$ are the acoustic velocity and pressure fluctuations,
respectively, $\zeta$ is the frequency dependent damping, which accounts for
the acoustic dissipation, $\dot{q}$ is the nonlinear heat release model, and
$\delta(x-x_{f})$ is a Dirac delta at the heat source location $x_{f}$. The
solution to the governing PDEs is expressed by a Galerkin decomposition on
$n_{g}$ acoustic eigenfunctions as [46]

$u(x,t)=\sum_{j=1}^{n_{g}}\eta_{j}(t)\cos(j\pi x),\quad
p(x,t)=-\sum_{j=1}^{n_{g}}\mu_{j}(t)\sin(j\pi x),$ (16)

where $j\pi$ are the non-dimensional angular frequencies of the purely
acoustic modes [47]. We denote the vectors of velocity and pressure Galerkin
amplitudes by
$\mathbf{\bm{\eta}}=[\eta_{1};\;\eta_{2};\;\dots;\;\eta_{n_{g}}]$ and
$\mathbf{\bm{\mu}}=[\mu_{1};\;\mu_{2};\;\dots;\;\mu_{n_{g}}]$, respectively.
By substituting pressure and velocity variables in the PDEs (15) with their
Galerkin decompositions in (16) and projecting the dynamics onto the Galerkin
modes, the dynamics of the Galerkin variables $\eta_{j}$ and $\mu_{j}$ are
described by a $2n_{g}$ dimensional system of ordinary differential equations
(ODEs)

$\displaystyle\dot{\eta}_{j}-\mu_{j}j\pi=0,$ (17a)
$\displaystyle\dot{\mu}_{j}+\eta_{j}j\pi+\zeta_{j}\mu_{j}+2\dot{q}\sin(j\pi
x_{f})=0,$ (17b)

where the modal damping is modelled by $\zeta_{j}=c_{1}j^{2}+c_{2}j^{1/2}$,
and $\dot{q}$ is the heat release rate described by a nonlinear time-delay
model [48]

$\dot{q}=\beta\left(\sqrt{|1+u(x_{f},t-\tau)|}-1\right),$ (18)

where $\beta$ is the heat release strength, and $\tau$ is the time delay
between the velocity at the heat source and heat release. The system exhibits
different nonlinear solutions, e.g., bifurcation from fixed point to limit
cycle, when the parameters $x_{f}$, $\beta$ and $\tau$ are varied. In the
remainder of this paper, we fix $x_{f}=0.2$, which causes thermoacoustic
instabilities to occur, and analyse the dynamics as functions $\beta$ and
$\tau$. The derivative of the nonlinear heat release model (18) with respect
to $u(x_{f},t-\tau)$ is not defined for $u(x_{f},t-\tau)=-1$, which means that
the Jacobian and consequently the sensitivity of the system can not be
calculated around this value. Therefore, the heat-release model (18) is
modified to have smooth dynamics by approximating a small-neighbourhood of
$u(x_{f},t-\tau)=-1$ with a fourth-degree polynomial [6]. Further, the time-
delayed problem is reformulated as an advection equation via the dummy
variable $v$ to make the problem Markovian [6]

$\frac{\partial v}{\partial t}+\frac{1}{\tau}\frac{\partial v}{\partial
X}=0,\quad 0\leq X\leq 1,\quad v(X=0,t)=u(x_{f},t).$ (19)

The PDE (19) is discretized using a Chebyshev spectral method with $n_{c}$
points such that the dimension of the system is $2n_{g}+n_{c}$ and thus, can
be solved with a standard time-marching scheme. For the thermoacoustic system,
the objective functional $\mathcal{J}$ is the time-averaged acoustic energy,
which relates to the average energy of the undesired oscillations, given by
[e.g. 6]

$E_{ac}(t)=\int_{0}^{1}\frac{1}{2}(u^{2}(x,t)+p^{2}(x,t))dx=\frac{1}{4}\sum_{j=1}^{n_{g}}(\eta_{j}^{2}(t)+\mu_{j}^{2}(t)).$
(20)

##  V Thermoacoustic Echo State Network

Figure 1: Schematic of the proposed thermoacoustic Echo State Network (T-ESN)
in closed-loop configuration. T-ESN hard constrains the physical knowledge
about the nonlinear and time-delayed thermoacoustics.

We propose a parameter-aware ESN with hard-constrained physics for
thermoacoustics [49]. In the context of machine learning, hard constrainining
means embedding prior knowledge into the architecture of the neural network,
as opposed to soft-constraining, which refers to imposing prior knowledge as a
penalty during training [23]. This approach makes neural networks more robust
and generalizable on unseen data [49]. First of all, our aim is to develop a
parameter-aware ESN that learns the dynamics of the Galerkin amplitudes (16)
from the thermoacoustic system introduced in Sec. IV. Initially, we consider
two cases without any hard constraining: (i) the Galerkin amplitudes as the
input and output vectors
$\mathbf{y}_{in}(i)=\hat{\mathbf{y}}(i)=[\hat{\mathbf{\bm{\eta}}}(i);\;\hat{\mathbf{\bm{\mu}}}(i)]$,
and (ii) the Galerkin amplitudes and the advection variables as the input and
output vectors
$\mathbf{y}_{in}(i)=\hat{\mathbf{y}}(i)=[\hat{\mathbf{\bm{\eta}}}(i);\;\hat{\mathbf{\bm{\mu}}}(i);\;\hat{\mathbf{v}}(i)]$.
In these cases, we feed both $\beta$ and $\tau$ as parameters through the
input channel, $\mathbf{p}=[\beta;\;\tau]$. However, neither of these options
is inherently time delayed. Therefore, we design an ESN for time-delayed
systems and refer to it as thermoacoustic ESN (T-ESN). During closed-loop
prediction, we compute the delayed velocity at the flame location from past
predictions of the velocity Galerkin amplitudes. We then feed the delayed
velocity into the network as an input for the prediction of the next time
step. One step of T-ESN is

$\displaystyle u_{f}(i-N_{\tau})$ $\displaystyle=\sum_{j=1}^{n_{g}}\cos(j\pi
x_{f})\hat{\eta}_{j}(i-N_{\tau}),$ (21a) $\displaystyle\mathbf{r}(i+1)$
$\displaystyle=(1-\alpha)\mathbf{r}(i)+\alpha\tanh(\mathbf{W}_{in}[\mathbf{y}_{in}(i);u_{f}(i-N_{\tau});\sigma_{\beta}(\beta-
k_{\beta})]+\mathbf{W}\mathbf{r}(i)),$ (21b)
$\displaystyle\hat{\mathbf{y}}(i+1)$
$\displaystyle=[\hat{\mathbf{\bm{\eta}}}(i+1);\;\hat{\mathbf{\bm{\mu}}}(i+1)]=\mathbf{W}_{out}\mathbf{r}(i+1),$
(21c)

where $N_{\tau}=\mathrm{int}(\tau/\Delta t)$, which implies that the ESN
handles delays that are integer multiples of the time step, as it is a
discrete-time map. (To extend this, one can use interpolation to approximate
the velocity at the desired time delay, however this is beyond the scope of
this paper.) A schematic of the T-ESN in closed-loop configuration is shown in
Fig. 1.

The tangent linear problem for computing the gradient with respect to the
parameters becomes

$\displaystyle\frac{d\mathcal{J}}{d\mathbf{p}}$
$\displaystyle=\sum_{i=1}^{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{r}(i))}{\partial\mathbf{r}(i)}\mathbf{Q}(i),$
(22a) $\displaystyle\mathbf{Q}(i+1)$
$\displaystyle=\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{p}}+\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{r}(i)}\mathbf{Q}(i)+\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{r}(i-N_{\tau})}\mathbf{Q}(i-N_{\tau}),$
(22b) $\displaystyle\mathbf{Q}(\leq 0)$ $\displaystyle=\mathbf{0}.$ (22c)

The adjoint of the time-delayed system is time-advanced

$\displaystyle\frac{d\mathcal{J}}{d\mathbf{p}}$
$\displaystyle=\sum_{i=1}^{N}\mathbf{q}^{+\top}(i)\frac{\partial\mathbf{r}(i)}{\partial\mathbf{p}},$
(23a) $\displaystyle\mathbf{q}^{+}(i)$
$\displaystyle=\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{r}(i))}{\partial\mathbf{r}(i)}^{\top}+\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{r}(i)}^{\top}\mathbf{q}^{+}(i+1)+\frac{\partial\mathbf{r}(i+1+N_{\tau})}{\partial\mathbf{r}(i)}^{\top}\mathbf{q}^{+}(i+1+N_{\tau}),$
(23b) $\displaystyle\mathbf{q}^{+}(N)$
$\displaystyle=\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{r}(N))}{\partial\mathbf{r}(N)}^{\top},$
(23c) $\displaystyle\mathbf{q}^{+}(\geq N)$ $\displaystyle=\mathbf{0}.$ (23d)

The derivation of the time-delayed adjoint is provided in Appendix A. The
time-delayed Jacobian is provided in Appendix B.

To obtain the sensitivity to time delay, we need

$\frac{\partial\mathbf{r}(i+1)}{\partial\tau}=\frac{\partial\mathbf{r}(i+1)}{\partial
u_{f}(i-\tau)}\frac{\partial u_{f}(i-\tau)}{\partial\tau}$ (24)

In the delayed formulation of the T-ESN, the delayed velocity $u_{f}(t-\tau)$
(21a) does not continuously depend on $\tau$. We approximate the gradient
$\partial u_{f}(t-\tau)/\partial\tau$ with a central finite difference

$\frac{\partial
u_{f}(i-\tau)}{\partial\tau}\approx\frac{u_{f}(i-(N_{\tau}+1))-u_{f}(i-(N_{\tau}-1))}{2\Delta
t}.$ (25)

Motivated by thermoacoustic physics, we design the input weights matrix based
on the following observations; (i) the acoustics are linear and the modes are
sparsely connected, (ii) the heat release is nonlinear, (iii) the heat release
strength $\beta$ is only coupled with $u_{f}(t-\tau)$

$\mathbf{W}_{in}=\begin{bmatrix}\mathbf{w}_{1,1}&0&0&0&0&0\\\
0&\mathbf{w}_{2,2}&0&0&0&0\\\ 0&0&\ddots&0&0&0\\\
0&0&0&\mathbf{w}_{2n_{g},2n_{g}}&0&0\\\
0&0&0&0&\mathbf{w}_{f}&\mathbf{w}_{\beta}\end{bmatrix},$ (26)

where
$\mathbf{w}_{j,j}\in\mathbb{R}^{\mathrm{int}(n_{r}/(2n_{g}+1))},\;j=1,\dots,2n_{g},\mathbf{w}_{f}\in\mathbb{R}^{\mathrm{int}(n_{r}/(2n_{g}+1))},\mathbf{w}_{\beta}\in\mathbb{R}^{\mathrm{int}(n_{r}/(2n_{g}+1))}$
The $\tanh$ activation acts approximately linearly near zero, so the terms
that linearly contribute to the dynamics have a smaller scaling, whereas the
nonlinear terms have a larger scaling. Therefore, we scale the input weights
that multiply the linear acoustics with one scaling i.e.,
$\mathbf{w}_{j,j}=\sigma_{in}\tilde{\mathbf{w}}_{j,j}$, and the input weights
that multiply the delayed velocity that contributes to the nonlinear heat
release term with another scaling, i.e.,
$\mathbf{w}_{f}=\sigma_{f}\tilde{\mathbf{w}}_{f}$, because we expect that
$\sigma_{f}>\sigma_{in}$ ($\tilde{\mathbf{w}}_{j,j}$, $\tilde{\mathbf{w}}_{f}$
and $\mathbf{w}_{\beta}$ are drawn from the uniform distribution
$\sim\mathcal{U}(-1,1)$ as before). While doing so we do not impose the
linearity, but rather allow the machine to discover this via the optimization
of the scaling terms. Separately, the parameter $\beta$ is shifted and scaled
by $\sigma_{\beta}(\beta-k_{\beta})$. The scalings are determined by the
hyperparameter search. In summary, we hard constrain the physical knowledge
about the thermoacoustics in the network architecture in two ways: (i) by
making the network explicitly time-delayed, and (ii) by designing the input
weights matrix.

##  VI Learning thermoacoustics and parametric sensitivities

In this section, we learn the dynamics of the Rijke tube parameterized by the
heat release strength, $\beta$, and the time delay, $\tau$. First, we focus on
limit cycle regimes and find that the thermoacoustic ESN (T-ESN) outperforms
other parameter-aware ESN architectures in terms of generalizability over
unseen regimes. Then, using the trained T-ESN, we infer the adjoint
sensitivities of the acoustic energy to these parameters. We compare the time
series predictions and inferred adjoint sensitivities to the ground truth. The
ground truth is obtained by integrating the original system described by the
ODEs (17), and its adjoint equations (5), which requires the computation of
the Jacobian of (17).

###  VI.1 Training and validation

For training and validation, we create a dataset from time-delays,
$\tau=\\{0.1,0.15,0.2,0.25,0.3\\}$ and heat release strengths,
$\beta=\\{1.0,2.0,3.0,4.0,5.0\\}$, giving a total of 25 thermoacoustic
regimes. These parameters generate limit cycle solutions except for
$(\beta=5.0,\tau=0.15)$, which produces a quasiperiodic solution. We fix the
flame location $x_{f}=0.2$, modal damping coefficients $c_{1}=0.1$ and
$c_{2}=0.06$, the number of Galerkin modes $n_{g}=4$, and the number of
Chebyshev discretization $n_{c}=10$ [6]. With four Galerkin modes, we observe
the rich nonlinear behaviour associated with thermoacoustics. We generate the
data by time-marching the ODEs (17) using the scipy.odeint tool with a time
step of 0.001, which we upsample to 0.01 for the ESN. We remove a transient of
200 time units, select the next 4 time units as washout and 8 time units as
training data. (We found that increasing the training length further does not
improve the performance whereas it increases the computational effort.) We
randomly split this dataset into 20 regimes for training and 5 regimes for
validation of the hyperparameters. During validation, we evaluate the closed-
loop performance of 5 random realisations of the ESN, i.e., realisations of
the randomly generated input and state matrices, over 5 trajectories, or
folds, of length 4 time units (approximately 2 periods) starting from random
time steps. As error metric, we compute relative $\ell_{2}$-error between the
prediction, $\hat{\mathbf{y}}$ and the ground truth $\mathbf{y}$

$\epsilon=\frac{\sqrt{\sum_{i=1}^{N_{v}}\sum_{j=1}^{n_{y}}(y_{j}(i)-\hat{y}_{j}(i))^{2}}}{\sqrt{\sum_{i=1}^{N_{v}}\sum_{j=1}^{n_{y}}y_{j}^{2}(i)}},$
(27)

where $N_{v}$ is the number of validation steps, which we average over the
number of folds, number of validation regimes, and number of realisations. We
employ a Bayesian optimization scheme using the skopt.gp_minimize tool to
determine the set of hyperparameters that have the minimal validation error
[24].

We consider three networks: (i) the standard parameter-aware ESN (9) with the
Galerkin amplitudes as the input and output vectors
$\mathbf{y}_{in}(i)=\hat{\mathbf{y}}(i)=[\hat{\mathbf{\bm{\eta}}}(i);\;\hat{\mathbf{\bm{\mu}}}(i)]$,
(ii) we augment the input and output vectors with the advection variable
$\mathbf{y}_{in}(i)=\hat{\mathbf{y}}(i)=[\hat{\mathbf{\bm{\eta}}}(i);\;\hat{\mathbf{\bm{\mu}}}(i);\;\hat{\mathbf{v}}(i)]$,
and (iii) thermoacoustic ESN (T-ESN) (21). The standard and the advection-
augmented parameter-aware ESNs have both $\beta$ and $\tau$ as explicit
parameters in the input channel, i.e., in Eq. (9) $\mathbf{p}=[\beta;\;\tau]$.
The T-ESN has $\beta$ as a parameter in the input channel and $\tau$ is
modelled as a time-delay in the delayed velocity term in the input channel
$u_{f}(t-\tau)$ (21a). For comparison, we fix the reservoir size to 1200 and
connectivity of the reservoir state variables to 20. The optimal
hyperparameters are provided in Table 1. The first row of Table 1 shows the
ranges from which the hyperparameters are sampled during the validation;
$\rho$, $\sigma_{in}$, $\sigma_{f}$, $\sigma_{\beta}$, $\sigma_{\tau}$,
$\alpha$ and $\lambda$ are sampled uniformly on a $\log_{10}$ scale, while
$k_{\beta}$ and $k_{\tau}$ are sampled uniformly.

Table 1: Optimal hyperparameters of (i) standard parameter-aware ESN (9), (ii) advection-augmented parameter-aware ESN (9), (iii) thermoacoustic ESN (21) trained on $\beta=\\{1.0,2.0,3.0,4.0,5.0\\}$ and $\tau=\\{0.1,0.15,0.2,0.25,0.3\\}$ regimes. Refer to the Eqs. (9), and (21) for the hyperparameters. | $\rho$ | $\sigma_{in}$ | $\sigma_{f}$ | $\sigma_{\beta}$ | $k_{\beta}$ | $\sigma_{\tau}$ | $k_{\tau}$ | $\alpha$ | $\lambda$
---|---|---|---|---|---|---|---|---|---
| 0.01 – 1.0 | 0.01 – 2.0 | 0.01 – 2.0 | 0.01 – 2.0 | -10.0 – 10.0 | 0.01 – 2.0 | -1.0 – 1.0 | 0.01 – 1.0 |  $10^{-6}$ – $10^{-1}$
ESN (i) | 0.0363 | 0.3476 | - | 0.0152 | -5.3759 | 0.7741 | -0.2063 | 0.0597 | $10^{-5}$
ESN (ii) | 0.0367 | 0.3293 | - | 0.0162 | -5.4293 | 0.7319 | 0.72734 | 0.0627 | $10^{-5}$
ESN (iii) | 0.0125 | 0.0790 | 1.0501 | 0.1092 | -10.0 | - | - | 0.0363 | $10^{-6}$

###  VI.2 Prediction performance

Figure 2: Short-term prediction performance of (a) standard parameter-aware
ESN (9), (b) advection-augmented parameter-aware ESN (9), (c) thermoacoustic
ESN (21). Relative $\ell_{2}$ error between prediction and ground truth of
$[\mathbf{\bm{\eta}};\;\mathbf{\bm{\mu}}]$ over 20 time units on a $\log_{10}$
scale.

The computed adjoint sensitivities are only as accurate as the model we have.
Therefore, the first task of the ESN is to accurately emulate the dynamics of
regimes both seen and unseen during training. This includes two subtasks: (i)
time-accurate prediction from a given initial condition, and (ii) convergence
to the correct attractor. The latter is necessary for long-term performance
because we are interested in the sensitivity of long-time averaged quantities.
Figure 2 shows the short-term prediction performance of the standard,
advection-augmented and thermoacoustic ESNs over a grid of $\beta$ and $\tau$
values. For a given regime, after a washout stage, the ESN is run in closed-
loop for prediction. This means that the network is forecasting the dynamics.
The relative $\ell_{2}$ error (27) between prediction and ground truth of
$[\mathbf{\bm{\eta}};\;\mathbf{\bm{\mu}}]$ over 20 time units is averaged over
5 folds starting from random time steps, and 5 realisations of the ESN.
Standard parameter-aware ESN is not capable of sustaining the limit cycle
solutions for approximately 10 periods even for the training regimes, whereas
the augmentation of the advection variables $\mathbf{v}$ in the input and
output significantly improves the performance. Notably, the T-ESN, generalizes
to unseen regimes. Owing to the nonlinearity in the system, the prediction
performance decreases as we move away from the training regimes.

In Fig. 3, we show the performance of the T-ESN on (a) a limit cycle regime
seen during training ($\beta=2.0,\;\tau=0.25$), and (b) an unseen limit cycle
regime ($\beta=4.5,\;\tau=0.12$). Based on Fig. 2, we choose regime (b) as an
example of an unseen regime with higher relative $\ell_{2}$ error compared to
the other limit cycle regimes. In these figures, we assess the performance by
(from left to right) short-term time-series prediction for 20 time units after
washout, the probability density function (PDF), i.e., statistics, and
amplitude spectrum of the first pressure mode $\mu_{1}$, and the acoustic
energy $E_{ac}$ (20). For the PDF and amplitude spectrum, we feed the initial
condition of zeros except $\eta_{1}(t_{0})=1.0$ in the washout stage, let the
ESN evolve in closed-loop, discard a transient of 200 time units and use the
next 1000 time units. This procedure and results verify that the ESN has
learned the correct limit cycle attractor because it converges to this set of
states from an initial condition that is not in the training set.

Figure 3: Thermoacoustic Echo State Network learns limit cycle regimes both
seen and unseen during training. Short-term prediction, probability
distribution function (PDF), and amplitude spectrum of the first pressure mode
and acoustic energy of (a-f) a limit cycle regime seen during training
($\beta=2.0,\;\tau=0.25$) and (g-l) an unseen limit cycle regime
($\beta=4.5,\;\tau=0.12$). In the PDFs, the mean of an ensemble of 5
realisations of the ESN is shown, one standard deviation from the mean is
plotted in lighter colour.

###  VI.3 Computation of adjoint sensitivities

Figure 4: Adjoint sensitivities of the true system and adjoint sensitivities
inferred by the thermoacoustic Echo State Network over an integration time of
100 time units. Top row, sensitivity of the time-averaged acoustic energy to
heat release strength $\beta$, $dJ/d\beta$ at fixed (a) $\tau=0.07$, (b)
$\tau=0.12$, (c) $\tau=0.22$, (d) $\tau=0.32$. Bottom row, sensitivity of the
time-averaged acoustic energy to time-delay $\tau$, $dJ/d\tau$ at fixed (e)
$\beta=1.25$, (f), $\beta=2.5$, (g)$\beta=3.75$, (h) $\beta=4.5$. The mean of
the sensitivity of an ensemble of 5 realisations of the ESN is shown, one
standard deviation from the mean is plotted in lighter colour.

We compare the adjoint sensitivities computed by the T-ESN following the
procedure described in Sec. III with the adjoint sensitivities obtained from
the original system described by Eq. (17). We refer to the latter as the
ground truth, which is computed by deriving the Jacobian of the ODEs (17) and
integrating the adjoint equations (5). The sensitivities are shown for varying
$\beta$ and $\tau$ values in Fig. 4. For ease of results, we show the
sensitivity to $\beta$ when $\beta$ is varied and $\tau$ is fixed, and vice
versa. In the washout stage, we feed the initial condition of zeros except
$\eta_{1}(t_{0})=1.5$, assuming that no data is available from the original
system for the washout. This value of $\eta_{1}(t_{0})$ makes sure that the
trajectories converge to the stable limit cycle and not the stable fixed point
in the multistable regimes. After the washout and the transient, we solve the
adjoint equations around a trajectory of 100 time units, for which we find
that the adjoint sensitivity converges to a steady-state value. The learned
adjoint sensitivities match closely with the ground truth, while the results
are also robust to different realisations of the network. Even though the
training data was generated from discrete values of $(\beta,\tau)$ on a grid
(Fig. 2), we are able to obtain sensitivity information from a range of
$(\beta,\tau)$ values that are unseen during the training using a single
instance of the T-ESN. The smoothness of the sensitivity curve also indicates
that the dynamics of the T-ESN varies smoothly with the parameters. The
accuracy decreases as we move away from the training regimes, e.g., beyond
$\beta=5.0$ in the$dJ/d\beta$ plot at $\tau=0.32$, but the ESN captures the
trend of the sensitivities.

Fig. 4 also shows that the T-ESN infers the existence of bifurcations to fixed
point solutions. The fixed point of the Rijke tube is the zero vector, i.e.,
the mean flow, which means that the sensitivity is zero for the fixed point
regimes, e.g., at fixed $\tau=0.22$ (Fig. 4(c)) for $\beta<0.45$. Albeit the
parameter value at which the bifurcation occurs might have a small error
compared to the true system. For example, at fixed $\tau=0.22$ the Hopf
bifurcation occurs at $\beta=0.45$, whereas 4 out of 5 realisations of the
T-ESN predict it as $\beta\approx 0.35$, and the remaining realisation
predicts it as $\beta\approx 0.30$.

##  VII Gradient-based optimization

Figure 5: Gradient-based optimization of the time-averaged acoustic energy
starting from different initial $(\beta,\tau)$ using adjoint sensitivities
computed by a single instance of the thermoacoustic Echo State Network. In
(a), example optimization paths starting from different initial $(\beta,\tau)$
are plotted on top of the landscape of the objective functional $\mathcal{J}$.
The values $\mathcal{J}$ takes during the iterations of optimization are
plotted in (b).

Figure 6: Gradient-based optimization of time-averaged acoustic energy marked
as Path 3 in Fig. 5. Shown (a) iteration 0; $(\beta=4.0,\tau=0.25)$, (b)
iteration 2; $(\beta=3.57,\tau=0.17)$, (c) iteration 4;
$(\beta=2.85,\tau=0.1)$, and (d) iteration 8; $(\beta=1.64,\tau=0.07)$. From
left to right; the first two pressure modes in the phase space, short-term
timeseries prediction of the acoustic energy, and long-term amplitude spectrum
of the acoustic energy. The optimization follows a path along which the
amplitude of the oscillations and the acoustic energy decreases, while the ESN
accurately predicts the changing dynamics.

We apply the inferred adjoint sensitivities to a gradient-based optimization
framework. One realisation of the T-ESN suffices for the computation of the
gradients and optimization for all parameters. In this work, the goal of the
optimization is to determine the set of parameters $(\beta,\tau)$ that
minimizes the time-averaged acoustic energy (20). We restrict the domain to
$\beta=[0.5,5.5],\;\tau=[0.05,0.35]$, which prevents the parameters to take
negative (non-physical) values, and the system to bifurcate to quasiperiodic
or chaotic regimes. We wish to avoid this during gradient-based optimization
because first, the acoustic energy tends to display discontinuous jumps at the
parameter values where the system bifurcates [6], which can lead to undefined
gradients between different nonlinear regimes. Second, while the adjoint
sensitivity can be computed for limit-cycle and quasiperiodic regimes for
arbitrary integration times, for chaotic regimes the integration time is
limited before the adjoint system becomes unstable, and the sensitivities
diverge [50, 6]. Another consideration is the respective orders of magnitudes
of the parameters and the sensitivities to them. From Fig. 4, $\beta$ and
$\tau$ as well as $dJ/d\beta$ and $dJ/d\beta$ have different orders of
magnitudes, i.e., $\beta\sim O(1),\;\tau\sim O(10^{-1})$. The difference in
the orders of magnitudes requires the parameters and their respective
gradients to be normalized, such that
$\breve{\beta}=\beta/\bar{\beta},\;\breve{\tau}=\tau/\bar{\tau},\;dJ/d\breve{\beta}=\bar{\beta}dJ/d\beta,\;dJ/d\breve{\tau}=\bar{\tau}dJ/d\tau$,
where $\bar{(\cdot)}$ denotes the mean of the parameter in the domain
considered [6]. The optimization is performed with the steepest-descent method
using the normalized gradients with a step size of 0.2. (Since $\tau$ can only
take discrete values, we round it according to the time-step size of the ESN
before the next step of the optimization.) The optimization stops when (i) the
acoustic energy is near zero, $J_{i}<10^{-4}$, or (ii) the relative change in
acoustic energy is below a threshold, $(J_{i}-J_{i-1})/(J_{i-1})<10^{-4}$
where subscript $i$ denotes the $i^{th}$ iteration step. Iterations of
optimization starting from different initial $(\beta,\tau)$ are shown in Fig.
5. The optimal $(\beta,\tau)$ values are found to be in the fixed point
regimes, i.e., when the instability is suppressed. For the Rijke tube, the
fixed point is the zero vector (mean flow), and as such the acoustic energy,
i.e., the acoustic energy of the fluctuations from the mean flow, is equal to
zero at these values. This result shows the significance of the ability of the
ESN to infer fixed points. The process of optimization is illustrated for Path
3 in Fig 6. The optimization follows a path where the amplitude of the
oscillations and so the acoustic energy decreases, for which the ESN
accurately predicts the changing dynamics. Similar behaviour is observed for
other paths. An exceptional case is Path 4, which passes through a saddle in
the landscape of the objective function near $(\beta=4.5,\tau=0.15)$. The
change in the sign of the gradient, as shown in Fig. 4 in the sensitivity plot
for $\beta=4.5$ around $\tau=0.13$, causes the “zigzagging” path in this
region.

##  VIII Robutsness to noise

If the data is acquired from real-life experiments, the measurements may be
corrupted by noise. This is the case when the system equations are unknown and
we wish to learn sensitivities from experimental data. To test the performance
of the T-ESN under such a scenario, we consider zero-mean additive Gaussian
noise with a standard deviation that is equal to 5% of the standard deviation
of the true data (this corresponds to 26 dB signal-to-noise ratio). We repeat
the training on the noisy data using the hyperparameters chosen before (Tab.
1) and three Tikhonov coefficients, $\lambda=10^{-6}$ as before,
$\lambda=10^{-5}$ and $\lambda=10^{-3}$ to observe the effect of increasing
regularization. We show the performance of closed-loop prediction on a
parameter regime $(\beta=2.0,\tau=0.25)$, which is in the training set, in
Fig. 7. The T-ESN successfully learns the limit cycle attractor from noisy
data and does not become unstable when the training data is noisy. We compute
the adjoint sensitivities for the same test cases in Sec.VI.3, these are shown
in Fig. 8. The sensitivity computations are affected by the noisy training
especially in small $\tau$ regimes ($\tau\leq 0.1$) that are unseen during
training. These errors are alleviated by increasing the regularization.
However, using a very large regularization might decrease the performance, for
example in higher $\tau$ regimes ($\tau\geq 0.25$) as shown in Fig. 7.

Figure 7: Thermoacoustic Echo State Network learns the true limit cycle
attractors of the original system from noisy data. Short-term closed-loop
prediction on the training data from $(\beta=2.0,\tau=0.25)$ regime is shown
for the first pressure mode as timeseries (a) and for the first and second
pressure modes in the phase space (b). The amplitude spectrum (c) is computed
for the same time window as the training data (8 time units) for consistency
in frequency resolution and shown in logarithmic scale. Short-term prediction,
probability distribution function (PDF), and amplitude spectrum of the
acoustic energy are shown as in Fig. 3. Figure 8: Adjoint sensitivities of the
true system and adjoint sensitivities inferred by the thermoacoustic Echo
State Network when the training is done on noisy data. Top row, sensitivity of
the time-averaged acoustic energy to heat release strength $\beta$,
$dJ/d\beta$ at fixed (a) $\tau=0.07$, (b) $\tau=0.12$, (c) $\tau=0.22$, (d)
$\tau=0.32$. Bottom row, sensitivity of the time-averaged acoustic energy to
time-delay $\tau$, $dJ/d\tau$ at fixed (e) $\beta=1.25$, (f), $\beta=2.5$,
(g)$\beta=3.75$, (h) $\beta=4.5$. The mean of the sensitivity of an ensemble
of 5 realisations of the ESN is shown, one standard deviation from the mean is
plotted in lighter colour. Increasing Tikhonov coefficient regularizes the
training.

##  IX Bifurcations to nonlinear regimes

Table 2: Optimal hyperparameters of the thermoacoustic ESN (21) trained on data from $\beta=\\{1.0,2.0,3.0,4.0,5.0\\}$ and $\beta=\\{6.0,6.5,7.0,7.5,8.0\\}$ regimes ($\tau=\\{0.1,0.15,0.2,0.25,0.3\\}$). Refer to (21) for the hyperparameters. | $\rho$ | $\sigma_{in}$ | $\sigma_{f}$ | $\sigma_{\beta}$ | $k_{\beta}$ | $\alpha$ | $\lambda$
---|---|---|---|---|---|---|---
| 0.01 – 1.0 | 0.01 – 2.0 | 0.01 – 2.0 | 0.01 – 2.0 | -10.0 – 10.0 | 0.01 – 1.0 |  $10^{-6}$ – $10^{-1}$
ESN $\beta_{train}=[1,5]$ | 0.0125 | 0.0790 | 1.0501 | 0.1092 | -10.0 | 0.0363 | $10^{-6}$
ESN $\beta_{train}=[6,8]$ | 0.1602 | 0.0122 | 1.296 | 0.3800 | -4.008 | 0.6369 | $10^{-6}$

Figure 9: Thermoacoustic Echo State Network captures bifurcations to period-
doubling and quasi-periodic regimes. Prediction, statistics, and amplitude
spectrum of (a-f) an unseen limit cycle regime where period-doubling occurs
($\beta=7.5,\;\tau=0.3$), (g-l) an unseen quasi-periodic regime
($\beta=6.1,\;\tau=0.2$)

Figure 10: Thermoacoustic Echo State Network captures bifurcations to chaotic
regimes. Prediction, statistics, and amplitude spectrum of (a-f) an unseen
chaotic regime ($\beta=7.6,\;\tau=0.22$) with a Lyapunov time of 8.5, and
(g-l) an unseen chaotic regime ($\beta=8.7,\;\tau=0.23$) with a Lyapunov time
of 3.9. Figure 11: Bifurcation diagrams obtained from the true system (top
row), thermoacoustic ESN trained on $\beta=[1,5]$ (middle row), and trained on
$\beta=[6,8]$ (bottom row). The peaks of the first pressure mode $\mu_{1}$ and
the acoustic energy are plotted. Left half (a-f) shows varying $\beta$ for
fixed $\tau=0.22$, and right half (g-l) shows varying $\tau$ for fixed beta =
7.6.

In Sec. VI, we discuss that the ESN predicts bifurcations from limit cycle
oscillations to fixed points when the heat release strength $\beta$ is
decreased. On the other hand, if $\beta$ is increased, the thermoacoustic
system displays bifurcations from limit cycle oscillations to chaos via
period-doubling and quasiperiodic routes [6]. In this section, we investigate
whether the ESN can also predict these bifurcations and attractors. We train
and validate a separate T-ESN on data generated by larger values
$\beta=\\{6.0,6.5,7.0,7.5,8.0\\}$ (for brevity, we refer to this dataset as
$\beta_{train}=[6.0-8.0]$). The reservoir size and connectivity are 1200 and
20, respectively as before. The optimal hyperparameters found for this case
are provided in Table 2 along with the hyperparameters for the previous case
when the training was done with data from $\beta=\\{1.0,2.0,3.0,4.0,5.0\\}$
(for brevity, we refer to this dataset as $\beta_{train}=[1.0-5.0]$). The
faster changing dynamics, e.g., chaotic behaviour, in the larger $\beta$
regimes might be responsible for the increase in the spectral radius and leak
rate.

In Fig. 9(a-f) we show a limit cycle regime where period-doubling occurs
($\beta=7.5,\;\tau=0.3$), and in Fig. 9(g-l) a quasi-periodic regime
($\beta=6.1,\;\tau=0.2$). Chaotic regimes are characterized by the time at
which two nearby trajectories diverge exponentially. This characteristic
timescale is known as the Lyapunov time (LT) and is an indicator of the
predictability of a system. We examine two chaotic thermoacoustic regimes with
different Lyapunov times; in Fig. 10(a-f) a regime ($\beta=7.6,\;\tau=0.22$)
with a Lyapunov time of 8.5 time units and in Fig. 10(g-l) a regime
($\beta=8.7,\;\tau=0.23$) with a Lyapunov time of 3.9 time units. The ability
of the T-ESNs trained on different $\beta$ regimes to predict bifurcations is
further shown in Fig. 11 in the bifurcation diagrams for fixed $\tau=0.22$ and
$\beta=7.6$ cases. The ESN trained on $\beta_{train}=[1.0-5.0]$ can
qualitatively predict the bifurcations and capture the statistics and dominant
frequencies correctly, especially for the non-chaotic regimes. This ability of
the parameter-aware ESN to predict bifurcations to quasi-periodic or chaotic
regimes at unseen parameter values is in itself not an obvious result as the
network only had access to limit cycle regimes at lower $\beta$ values for
training and validation. For example, when we increase $\beta$ and a
bifurcation occurs, the acoustic energy increases with a discontinuous jump,
which indicates that the amplitude of the oscillations abruptly change.
Therefore, the prediction of these regimes becomes an extrapolation problem
because the range of data we are interested in is different from the range of
data the network was trained on. In this case, incorporating physical
knowledge via the T-ESN architecture enables extrapolation performance.
Training the T-ESN on $\beta_{train}=[6.0-8.0]$, which encompasses the regimes
that we are testing, increases accuracy of both short- and long-term
predictions. The bifurcation diagrams confirm that in the higher $\beta$
regimes, the T-ESN trained on $\beta_{train}=[1.0-5.0]$ can qualitatively
capture the bifurcations to quasiperiodic and chaotic regimes, however is not
as quantitatively accurate as the T-ESN trained on $\beta_{train}=[6.0-8.0]$.
(The reverse is true for the lower $\beta$ regimes.) When $\beta$ is further
increased beyond $\beta=10.0$, the quantitatively accurate performance of the
T-ESN trained on $\beta_{train}=[6.0-8.0]$ also starts diminishing.

In chaotic systems, due to the positive Lyapunov exponent, the adjoint system
is unstable and the adjoint sensitivities diverge [50, 6]. This means that we
can only integrate adjoint equations for short time windows. Using the T-ESN
trained on $\beta_{train}=[6.0-8.0]$, we compute the adjoint sensitivities of
short-time trajectories for the ($\beta=7.6,\;\tau=0.22$) regime. Figures
12(a,b) show the evolution of the adjoint sensitivities with integration time
for one trajectory for an ensemble of 5 T-ESNs. For short integration times,
the T-ESN predicts the chaotic sensitivities close to the original system, but
diverges as the integration time increases due to the chaotic nature of the
system. One approach to estimate the sensitivity of chaotic systems is by
taking the mean of an ensemble trajectories [50]. In Figs. 12(c,d) and (e,f),
we show the probability distributions and the mean of the adjoint
sensitivities for 2000 trajectories for integration times of 0.5 LT and 1 LT,
respectively. In order to obtain the PDFs, we let the T-ESN run autonomously
without any washout data from the original system. For an integration time of
0.5 LT, we can recover the PDF as well as the mean, while for 1.0 LT, low
statistical moments are matched.

Figure 12: Computation of chaotic adjoint sensitivities using thermoacoustic
ESN. The adjoint sensitivity for a single trajectory at
$(\beta=7.6,\;\tau=0.22)$ is shown on the left. The probability distribution
of the sensitivity is computed from 2000 short trajectories with integration
time of 0.5 LT (c,d) and 1.0 LT (e,f).

##  X Conclusion

Adjoint methods are established tools for gradient-based optimization.
Although adjoint methods offer an accurate and efficient gradient computation,
they are code specific. This means that an adjoint solver needs to be
derived/implemented for new problems of interest. In this work, we propose a
strategy to infer adjoint sensitivities from data (observables). We showcase
the method on a prototypical nonlinear time-delayed thermoacoustic system.
This system is described by a nonlinear wave equation. The parameters of
interest are the heat release strength and the time-delay at the heat source
location. By varying these parameters, the system displays fixed-point, limit-
cycle, quasiperiodic, and chaotic solutions. The cost functional is the time-
average acoustic energy, which we wish to minimize to suppress an instability.

First, we propose a parameter-aware Echo State Network (ESN). This network
learns the parameterized dynamics, i.e., it infers how the dynamics change
with varying parameters, thereby making predictions for parameters that are
not used during training. Second, we derive the adjoint parameter-aware Echo
State Network, which is used to compute the adjoint sensitivities of the
physical system. Third, we design a parameter-aware ESN motivated by the
physics of thermoacoustics (thermoacoustic ESN, or T-ESN in short) by (i)
making the ESN explicitly time-delayed, and (ii) designing the input matrix in
consideration of acoustic physics. We extend the derivation of the adjoint to
the time-delayed case. We find that the T-ESN accurately predicts the dynamics
of different nonlinear regimes, and successfully infers the adjoint
sensitivities. These results are validated against the ground truth generated
from integration of the direct and adjoint equations of the original system.
We optimize the acoustic energy with steepest gradient-descent, in which the
gradients are computed by only one T-ESN. Once the T-ESN has been trained, the
algorithm runs autonomously, i.e., we do not run any additional simulations of
the original model. As the optimization proceeds with the iterations, the
acoustic energy and the amplitude of the oscillations decrease, finally
reaching a fixed-point solution with zero acoustic energy. Fourth, we
investigate the robustness to noisy data and chaotic regimes. The T-ESN can
qualitatively predict bifurcations to nonlinear regimes beyond the training
range of parameter values. Quantitative accuracy increases when the range of
training parameters covers the regime in which we are interested. This work
opens up possibilities for data-driven optimization without developing code-
specific adjoint solvers.

###### Acknowledgements.

This research has received financial support from the ERC Starting Grant No.
PhyCo 949388 and UKRI AI for Net Zero grant EP/Y005619/1. L.M. is also
grateful for the support from the grant EU-PNRR YoungResearcher TWIN ERC-
PI_0000005.

##  Appendix A Adjoints of discrete-time dynamical systems with and without
time delay

###  A.1 Adjoint of a discrete-time system

The following derivation of the adjoint sensitivity applies to any discrete
map of choice that can be expressed as

$\mathbf{x}(i)=\mathbf{f}(\mathbf{x}(i-1),\mathbf{p}),$ (28)

where $\mathbf{x}(i)\in\mathbb{R}^{n_{x}}$ is the state vector at time step
$i$. Equation (28) is rewritten as a constraint

$\mathbf{F}(\mathbf{x}(i),\mathbf{x}(i-1),\mathbf{p})\coloneqq\mathbf{x}(i)-\mathbf{f}(\mathbf{x}(i-1),\mathbf{p})=0.$
(29)

We consider a time-averaged quantity as the objective functional

$\mathcal{J}(\mathbf{x},\mathbf{p})\coloneqq\frac{1}{N}\sum_{i=1}^{N}\tilde{\mathcal{J}}(\mathbf{x}(i),\mathbf{p}).$
(30)

Following [5], the problem of minimizing the objective functional (30) under
the constraints (29) can be formulated as a Lagrangian optimization problem

$\mathcal{L}\coloneqq\mathcal{J}(\mathbf{x},\mathbf{p})-\langle\mathbf{q}^{+},\mathbf{F}\rangle,$
(31)

where $\mathbf{q}^{+}$ are the Lagrange multipliers. The inner product
$\langle\mathbf{q}^{+},\mathbf{F}\rangle$ is given by

$\langle\mathbf{q}^{+},\mathbf{F}\rangle\coloneqq\sum_{i=1}^{N}\mathbf{q}^{+\top}(i)\mathbf{F}(i),$
(32)

where we denote
$\mathbf{F}(i)=\mathbf{F}(\mathbf{x}(i),\mathbf{x}(i-1),\mathbf{p})$ for
brevity. Because $d\mathcal{J}/d\mathbf{p}=d\mathcal{L}/d\mathbf{p}$, we seek
the gradient $d\mathcal{L}/d\mathbf{p}$

$\frac{d\mathcal{L}}{d\mathbf{p}}=\frac{1}{N}\sum_{i=1}^{N}\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{p}}+\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{x}(i)}\frac{d\mathbf{x}(i)}{d\mathbf{p}}-\sum_{i=1}^{N}\mathbf{q}^{+\top}(i)\left(\frac{\partial\mathbf{F}(i)}{\partial\mathbf{p}}+\frac{\partial\mathbf{F}(i)}{\partial\mathbf{x}(i)}\frac{d\mathbf{x}(i)}{d\mathbf{p}}+\frac{\partial\mathbf{F}(i)}{\partial\mathbf{x}(i-1)}\frac{d\mathbf{x}(i-1)}{d\mathbf{p}}\right).$
(33)

Without loss of generality, we assume that the objective function does not
explicitly depend on $\mathbf{p}$ and therefore
$\partial\tilde{\mathcal{J}}/\partial\mathbf{p}$ is zero. Writing out the
terms of the summation, we have

$\begin{split}\frac{d\mathcal{L}}{d\mathbf{p}}=\sum_{i=1}^{N}&-\mathbf{q}^{+\top}(i)\frac{\partial\mathbf{F}(i)}{\partial\mathbf{p}}\\\
&+\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(1))}{\partial\mathbf{x}(1)}\frac{d\mathbf{x}(1)}{d\mathbf{p}}-\mathbf{q}^{+\top}(1)\frac{\partial\mathbf{F}(1)}{\partial\mathbf{x}(1)}\frac{d\mathbf{x}(1)}{d\mathbf{p}}-\mathbf{q}^{+\top}(1)\frac{\partial\mathbf{F}(1)}{\partial\mathbf{x}(0)}\frac{d\mathbf{x}(0)}{d\mathbf{p}}\\\
&+\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(2))}{\partial\mathbf{x}(2)}\frac{d\mathbf{x}(2)}{d\mathbf{p}}-\mathbf{q}^{+\top}(2)\frac{\partial\mathbf{F}(2)}{\partial\mathbf{x}(2)}\frac{d\mathbf{x}(2)}{d\mathbf{p}}-\mathbf{q}^{+\top}(2)\frac{\partial\mathbf{F}(2)}{\partial\mathbf{x}(1)}\frac{d\mathbf{x}(1)}{d\mathbf{p}}\\\
&\vdots\\\
&+\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(N))}{\partial\mathbf{x}(N)}\frac{d\mathbf{x}(N)}{d\mathbf{p}}-\mathbf{q}^{+\top}(N)\frac{\partial\mathbf{F}(N)}{\partial\mathbf{x}(N)}\frac{d\mathbf{x}(N)}{d\mathbf{p}}-\mathbf{q}^{+\top}(N)\frac{\partial\mathbf{F}(N)}{\partial\mathbf{x}(N-1)}\frac{d\mathbf{x}(N-1)}{d\mathbf{p}}.\end{split}$
(34)

Because of Eq. (28), the terms $\partial\mathbf{F}(i)/\partial\mathbf{x}(i)$
are equal to the identity matrix, and $d\mathbf{x}(0)/d\mathbf{p}$ is equal to
zero. We choose the Lagrange multipliers to eliminate the terms
$d\mathbf{x}(i)/d\mathbf{p}$, the size of which grow with the dimension of
$\mathbf{p}$. For this purpose, we gather the terms multiplying
$d\mathbf{x}(i)/d\mathbf{p}$, thereby obtaining the evolution equations (12)
for the Lagrange multipliers, i.e., the adjoint variables,

$\begin{split}&\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(1))}{\partial\mathbf{x}(1)}-\mathbf{q}^{+\top}(1)-\mathbf{q}^{+\top}(2)\frac{\partial\mathbf{F}(2)}{\partial\mathbf{x}(1)}=0\\\
&\vdots\\\
&\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(N))}{\partial\mathbf{x}(N)}-\mathbf{q}^{+\top}(N)=0,\end{split}$
(35)

where from Eq. (28),
$\partial\mathbf{F}(i)/\partial\mathbf{x}(i-1)=-\partial\mathbf{f}(i)/\partial\mathbf{x}(i-1)$.

###  A.2 Adjoint of a time-delayed discrete-time system

The discrete map for the time-delayed can be expressed as

$\mathbf{F}(\mathbf{x}(i),\mathbf{x}(i-1),\mathbf{x}(i-1-N_{\tau}),\mathbf{p})=0,$
(36)

where $N_{\tau}$ is the time-delay in discrete time steps. The Lagrangian
$\mathcal{L}$ is formulated similarly to Eq. (31). The gradient
$d\mathcal{L}/d\mathbf{p}$ is given by

$\frac{d\mathcal{L}}{d\mathbf{p}}=\frac{1}{N}\sum_{i=1}^{N}\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{p}}+\frac{\partial\tilde{\mathcal{J}}}{\partial\mathbf{x}(i)}\frac{d\mathbf{x}(i)}{d\mathbf{p}}-\sum_{i=1}^{N}\mathbf{q}^{+\top}(i)\left(\frac{\partial\mathbf{F}(i)}{\partial\mathbf{p}}+\frac{\partial\mathbf{F}(i)}{\partial\mathbf{x}(i)}\frac{d\mathbf{x}(i)}{d\mathbf{p}}+\frac{\partial\mathbf{F}(i)}{\partial\mathbf{x}(i-1)}\frac{d\mathbf{x}(i-1)}{d\mathbf{p}}+\frac{\partial\mathbf{F}(i)}{\partial\mathbf{x}(i-1-N_{\tau})}\frac{d\mathbf{x}(i-1-N_{\tau})}{d\mathbf{p}}\right).$
(37)

Writing out the terms of the summation, we have

$\begin{split}\frac{d\mathcal{L}}{dp}=\sum_{i=1}^{N}&-\mathbf{q}^{+\top}(i)\frac{\partial\mathbf{F}(i)}{\partial\mathbf{p}}\\\
&+\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(1))}{\partial\mathbf{x}(1)}\frac{d\mathbf{x}(1)}{d\mathbf{p}}-\mathbf{q}^{+\top}(1)\frac{\partial\mathbf{F}(1)}{\partial\mathbf{x}(1)}\frac{d\mathbf{x}(1)}{d\mathbf{p}}-\mathbf{q}^{+\top}(1)\frac{\partial\mathbf{F}(1)}{\partial\mathbf{x}(0)}\frac{d\mathbf{x}(0)}{d\mathbf{p}}-\mathbf{q}^{+\top}(1)\frac{\partial\mathbf{F}(1)}{\partial\mathbf{x}(-N_{\tau})}\frac{d\mathbf{x}(-N_{\tau})}{d\mathbf{p}}\\\
&+\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(2))}{\partial\mathbf{x}(2)}\frac{d\mathbf{x}(2)}{d\mathbf{p}}-\mathbf{q}^{+\top}(2)\frac{\partial\mathbf{F}(2)}{\partial\mathbf{x}(2)}\frac{d\mathbf{x}(2)}{d\mathbf{p}}-\mathbf{q}^{+\top}(2)\frac{\partial\mathbf{F}(2)}{\partial\mathbf{x}(1)}\frac{d\mathbf{x}(1)}{d\mathbf{p}}-\mathbf{q}^{+\top}(2)\frac{\partial\mathbf{F}(2)}{\partial\mathbf{x}(1-N_{\tau})}\frac{d\mathbf{x}(1-N_{\tau})}{d\mathbf{p}}\\\
&\vdots\\\
&+\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(2+N_{\tau}))}{\partial\mathbf{x}(2+N_{\tau})}\frac{d\mathbf{x}(2+N_{\tau})}{d\mathbf{p}}-\mathbf{q}^{+\top}(2+N_{\tau})\frac{\partial\mathbf{F}(2+N_{\tau})}{\partial\mathbf{x}(2+N_{\tau})}\frac{d\mathbf{x}(2+N_{\tau})}{d\mathbf{p}}-\mathbf{q}^{+\top}(2+N_{\tau})\frac{\partial\mathbf{F}(2+N_{\tau})}{\partial\mathbf{x}(1+N_{\tau})}\frac{d\mathbf{x}(1+N_{\tau})}{d\mathbf{p}}\\\
&-\mathbf{q}^{+\top}(2+N_{\tau})\frac{\partial\mathbf{F}(2+N_{\tau})}{\partial\mathbf{x}(1)}\frac{d\mathbf{x}(1)}{d\mathbf{p}}\\\
&\vdots\\\
&+\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(N))}{\partial\mathbf{x}(N)}\frac{d\mathbf{x}(N)}{d\mathbf{p}}-\mathbf{q}^{+\top}(N)\frac{\partial\mathbf{F}(N)}{\partial\mathbf{x}(N)}\frac{d\mathbf{x}(N)}{d\mathbf{p}}-\mathbf{q}^{+\top}(N)\frac{\partial\mathbf{F}(N)}{\partial\mathbf{x}(N-1)}\frac{d\mathbf{x}(N-1)}{d\mathbf{p}}\\\
&-\mathbf{q}^{+\top}(N)\frac{\partial\mathbf{F}(N)}{\partial\mathbf{x}(N-1-N_{\tau})}\frac{d\mathbf{x}(N-1-N_{\tau})}{d\mathbf{p}}.\end{split}$
(38)

Because of Eq. (28), the terms $\partial\mathbf{F}(i)/\partial\mathbf{x}(i)$
are equal to the identity matrix, and due to the initial conditions
$d\mathbf{x}(i\leq 0)/d\mathbf{p}$ are equal to zero. Eliminating
$d\mathbf{x}/d\mathbf{p}$ terms, the evolution of the adjoint variables is
given by

$\begin{split}&\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(1))}{\partial\mathbf{x}(1)}-\mathbf{q}^{+\top}(1)-\mathbf{q}^{+\top}(2)\frac{\partial\mathbf{F}(2)}{\partial\mathbf{x}(1)}-\mathbf{q}^{+\top}(2+N_{\tau})\frac{\partial\mathbf{F}(2+N_{\tau})}{\partial\mathbf{x}(1)}=0,\\\
&\vdots\\\
&\frac{1}{N}\frac{\partial\tilde{\mathcal{J}}(\mathbf{x}(N))}{\partial\mathbf{x}(N)}-\mathbf{q}^{+\top}(N)=0.\end{split}$
(39)

We set $\mathbf{q}^{+}(i>N)=0$ because there are no terms with
$d\mathbf{x}(i>N)/d\mathbf{p}$.

##  Appendix B Gradients of the Echo State Network

The gradients associated with the ESN are provided here. We first define,

$\tilde{\mathbf{r}}(i)=\frac{\mathbf{r}(i+1)-(1-\alpha)\mathbf{r}(i)}{\alpha},$
(40)

which is used to compute the derivative of the $\tanh$ term in the ESN
equations, e.g., (6a). The Jacobian and the partial derivative with respect to
the system’s parameters are

$\displaystyle\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{r}(i)}$
$\displaystyle=(1-\alpha)\mathbf{I}_{n_{r}\times
n_{r}}+\alpha\mathrm{diag}(1-\tilde{\mathbf{r}}^{2}(i))(\mathbf{W}_{in}^{y}\mathbf{W}_{out}+\mathbf{W}),$
(41) $\displaystyle\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{p}}$
$\displaystyle=\alpha\mathrm{diag}(1-\tilde{\mathbf{r}}^{2}(i))\mathbf{W}_{in}^{p}\mathrm{diag}(\mathbf{\bm{\sigma}}_{p}),$
(42)

where $\mathbf{I}$ denotes the identity matrix, and $\mathrm{diag}(\cdot)$
denotes a diagonal matrix that has $(\cdot)$ as its diagonal. For the
thermoacoustic ESN, we have the following gradients

$\displaystyle\frac{\partial\mathbf{r}(i+1)}{\partial u_{f}(i-N_{\tau})}$
$\displaystyle=\alpha\mathrm{diag}(1-\tilde{\mathbf{r}}^{2}(i))\mathbf{W}_{in}^{u_{f}}$
(43)
$\displaystyle\frac{\partial\mathbf{r}(i+1)}{\partial\mathbf{r}(i-N_{\tau})}$
$\displaystyle=\frac{\partial\mathbf{r}(i+1)}{\partial
u_{f}(i-N_{\tau})}\Psi\mathbf{W}_{out},$ (44)

where $\Psi=[\cos(\pi x_{f});\;\cos(2\pi x_{f});\;\dots;\;\cos(n_{g}\pi
x_{f});\;\mathbf{0}_{n_{g}}]^{\top}$ contains the Galerkin modes associated
with velocity and masks the Galerkin modes associated with pressure.

## References

  * Peter and Dwight [2010] J. E. Peter and R. P. Dwight, Numerical sensitivity analysis for aerodynamic optimization: A survey of approaches, Computers & Fluids 39, 373 (2010).
  * Sipp  _et al._ [2010] D. Sipp, O. Marquet, P. Meliga, and A. Barbagallo, Dynamics and Control of Global Instabilities in Open-Flows: A Linearized Approach, Applied Mechanics Reviews 63, 030801 (2010).
  * Luchini and Bottaro [2014] P. Luchini and A. Bottaro, Adjoint Equations in Stability Analysis, Annual Review of Fluid Mechanics 46, 493 (2014).
  * Camarri [2015] S. Camarri, Flow control design inspired by linear stability analysis, Acta Mechanica 226, 979 (2015).
  * Magri [2019] L. Magri, Adjoint Methods as Design Tools in Thermoacoustics, Applied Mechanics Reviews 71, 020801 (2019).
  * Huhn and Magri [2020] F. Huhn and L. Magri, Stability, sensitivity and optimisation of chaotic acoustic oscillations, Journal of Fluid Mechanics 882, A24 (2020).
  * Huhn and Magri [2022] F. Huhn and L. Magri, Gradient-free optimization of chaotic acoustics with reservoir computing, Physical Review Fluids 7, 014402 (2022).
  * Giles and Pierce [2000] M. B. Giles and N. A. Pierce, An Introduction to the Adjoint Approach to Design, Flow, Turbulence and Combustion 65, 393 (2000).
  * Gunzburger [2002] M. D. Gunzburger,  _Perspectives in Flow Control and Optimization_ (Society for Industrial and Applied Mathematics, 2002).
  * Giannetti and Luchini [2007] F. Giannetti and P. Luchini, Structural sensitivity of the first instability of the cylinder wake, Journal of Fluid Mechanics 581, 167 (2007).
  * Griewank and Walther [2000] A. Griewank and A. Walther, Algorithm 799: Revolve: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation, ACM Transactions on Mathematical Software 26, 19 (2000).
  * Pironneau [1974] O. Pironneau, On optimum design in fluid mechanics, Journal of Fluid Mechanics 64, 97 (1974).
  * Jameson [1988] A. Jameson, Aerodynamic design via control theory, Journal of Scientific Computing 3, 233 (1988).
  * Jameson  _et al._ [1998] A. Jameson, L. Martinelli, and N. Pierce, Optimum Aerodynamic Design Using the Navier-Stokes Equations, Theoretical and Computational Fluid Dynamics 10, 213 (1998).
  * Jameson [1999] A. Jameson, Re-Engineering the Design Process Through Computation, Journal of Aircraft 36, 36 (1999).
  * Mohammadi and Mohammadi [1997] B. Mohammadi and B. Mohammadi, Optimal shape design, reverse mode of automatic differentiation and turbulence, in  _35th Aerospace Sciences Meeting and Exhibit_ (American Institute of Aeronautics and Astronautics, Reno,NV,U.S.A., 1997).
  * Biava  _et al._ [2016] M. Biava, M. Woodgate, and G. N. Barakos, Fully Implicit Discrete-Adjoint Methods for Rotorcraft Applications, AIAA Journal 54, 735 (2016).
  * Xiao  _et al._ [2021a] C. Xiao, Y. Deng, and G. Wang, Deep-Learning-Based Adjoint State Method: Methodology and Preliminary Application to Inverse Modeling, Water Resources Research 57, e2020WR027400 (2021a).
  * Lukoševičius [2012] M. Lukoševičius, A Practical Guide to Applying Echo State Networks, in  _Neural Networks: Tricks of the Trade_, Vol. 7700, edited by G. Montavon, G. B. Orr, and K.-R. Müller (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012) pp. 659–686.
  * Grigoryeva and Ortega [2018] L. Grigoryeva and J.-P. Ortega, Echo state networks are universal, Neural Networks 108, 495 (2018).
  * Margazoglou and Magri [2023a] G. Margazoglou and L. Magri, Stability analysis of chaotic systems from data, Nonlinear Dynamics 111, 8799 (2023a).
  * Pathak  _et al._ [2018] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Model-Free Prediction of Large Spatiotemporally Chaotic Systems from Data: A Reservoir Computing Approach, Physical Review Letters 120, 024102 (2018).
  * Doan  _et al._ [2020a] N. Doan, W. Polifke, and L. Magri, Physics-informed echo state networks, Journal of Computational Science 47, 101237 (2020a).
  * Racca and Magri [2021] A. Racca and L. Magri, Robust Optimization and Validation of Echo State Networks for learning chaotic dynamics, Neural Networks 142, 252 (2021).
  * Lu  _et al._ [2017] Z. Lu, J. Pathak, B. Hunt, M. Girvan, R. Brockett, and E. Ott, Reservoir observers: Model-free inference of unmeasured variables in chaotic systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 041102 (2017).
  * Doan  _et al._ [2020b] N. A. K. Doan, W. Polifke, and L. Magri, Learning Hidden States in a Chaotic System: A Physics-Informed Echo State Network Approach, in  _Computational Science – ICCS 2020_, Vol. 12142, edited by V. V. Krzhizhanovskaya, G. Závodszky, M. H. Lees, J. J. Dongarra, P. M. A. Sloot, S. Brissos, and J. Teixeira (Springer International Publishing, Cham, 2020) pp. 117–123.
  * Pathak  _et al._ [2017] J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data, Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 121102 (2017).
  * Nóvoa  _et al._ [2024] A. Nóvoa, A. Racca, and L. Magri, Inferring unknown unknowns: Regularized bias-aware ensemble Kalman filter, Computer Methods in Applied Mechanics and Engineering 418, 116502 (2024).
  * Racca  _et al._ [2023] A. Racca, N. A. K. Doan, and L. Magri, Predicting turbulent dynamics with the convolutional autoencoder echo state network, Journal of Fluid Mechanics 975, A2 (2023).
  * Xiao  _et al._ [2021b] R. Xiao, L.-W. Kong, Z.-K. Sun, and Y.-C. Lai, Predicting amplitude death with machine learning, Physical Review E 104, 014205 (2021b).
  * Lieuwen and Yang [2006] T. C. Lieuwen and V. Yang,  _Combustion Instabilities In Gas Turbine Engines: Operational Experience, Fundamental Mechanisms, and Modeling_ (American Institute of Aeronautics and Astronautics, Reston ,VA, 2006).
  * Culick and Kuentzmann [2006] F. Culick and P. Kuentzmann, Unsteady Motions in Combustion Chambers for Propulsion Systems, NATO RTO-AG-AVT-039, AGARDograph (2006).
  * Rayleigh [1878] Rayleigh, THE EXPLANATION OF CERTAIN ACOUSTICAL PHENOMENA 1, Nature 18, 319 (1878).
  * Magri  _et al._ [2020] L. Magri, M. P. Juniper, and J. P. Moeck, Sensitivity of the Rayleigh criterion in thermoacoustics, Journal of Fluid Mechanics 882, R1 (2020).
  * Kabiraj  _et al._ [2012] L. Kabiraj, R. I. Sujith, and P. Wahi, Bifurcations of Self-Excited Ducted Laminar Premixed Flames, Journal of Engineering for Gas Turbines and Power 134, 031502 (2012).
  * Juniper and Sujith [2018] M. P. Juniper and R. Sujith, Sensitivity and Nonlinearity of Thermoacoustic Oscillations, Annual Review of Fluid Mechanics 50, 661 (2018).
  * Magri  _et al._ [2023] L. Magri, P. J. Schmid, and J. P. Moeck, Linear Flow Analysis Inspired by Mathematical Methods from Quantum Mechanics, Annual Review of Fluid Mechanics 55, 541 (2023).
  * Magri and Juniper [2013] L. Magri and M. P. Juniper, Sensitivity analysis of a time-delayed thermo-acoustic system via an adjoint-based approach, Journal of Fluid Mechanics 719, 183 (2013).
  * Mensah and Moeck [2017] G. A. Mensah and J. P. Moeck, Acoustic Damper Placement and Tuning for Annular Combustors: An Adjoint-Based Optimization Study, Journal of Engineering for Gas Turbines and Power 139, 061501 (2017).
  * Aguilar and Juniper [2018] J. G. Aguilar and M. P. Juniper, Adjoint Methods for Elimination of Thermoacoustic Oscillations in a Model Annular Combustor via Small Geometry Modifications, in  _Volume 4A: Combustion, Fuels, and Emissions_ (American Society of Mechanical Engineers, Oslo, Norway, 2018) p. V04AT04A054.
  * Orchini  _et al._ [2016] A. Orchini, G. Rigas, and M. P. Juniper, Weakly nonlinear analysis of thermoacoustic bifurcations in the Rijke tube, Journal of Fluid Mechanics 805, 523 (2016).
  * Margazoglou and Magri [2023b] G. Margazoglou and L. Magri, Data-Driven Stability Analysis of a Chaotic Time-Delayed System, in  _Computational Science – ICCS 2023_, Vol. 10476, edited by J. Mikyška, C. De Mulatier, M. Paszynski, V. V. Krzhizhanovskaya, J. J. Dongarra, and P. M. Sloot (Springer Nature Switzerland, Cham, 2023) pp. 406–413.
  * Jaeger and Haas [2004] H. Jaeger and H. Haas, Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication, Science 304, 78 (2004).
  * Roy  _et al._ [2022] M. Roy, S. Mandal, C. Hens, A. Prasad, N. V. Kuznetsov, and M. Dev Shrimali, Model-free prediction of multistability using echo state network, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 101104 (2022).
  * Juniper [2011] M. P. Juniper, Triggering in the horizontal Rijke tube: Non-normality, transient growth and bypass transition, Journal of Fluid Mechanics 667, 272 (2011).
  * Zinn and Lores [1971] B. T. Zinn and M. E. Lores, Application of the Galerkin Method in the Solution of Non-linear Axial Combustion Instability Problems in Liquid Rockets, Combustion Science and Technology 4, 269 (1971).
  * Magri and Juniper [2014] L. Magri and M. P. Juniper, Global modes, receptivity, and sensitivity analysis of diffusion flames coupled with duct acoustics, Journal of Fluid Mechanics 752, 237 (2014).
  * Heckl [1990] M. A. Heckl, Nonlinear Acoustic Effects in the Rijke Tube, Acustica 72, 63 (1990).
  * Ozan and Magri [2023] D. E. Ozan and L. Magri, Hard-constrained neural networks for modeling nonlinear acoustics, Physical Review Fluids 8, 103201 (2023).
  * Lea  _et al._ [2000] D. J. Lea, M. R. Allen, and T. W. N. Haine, Sensitivity analysis of the climate of a chaotic system, Tellus A 52, 523 (2000).