Adaptive Nonlinear Regulation Control of Thermoacoustic Oscillations in Rijke-Type Systems

Abstract

Adaptive nonlinear control of self-excited oscillations in Rijke-type thermoacoustic systems is considered. To demonstrate the methodology, a well-accepted thermoacoustic dynamic model is introduced, which includes arrays of sensors and monopole-like actuators. To facilitate the derivation of the adaptive control law, the dynamic model is recast as a set of nonlinear ordinary differential equations, which are amenable to control design. The control-oriented nonlinear model includes unknown, unmeasurable, nonvanishing disturbances in addition to parametric uncertainty in both the thermoacoustic dynamic model and the actuator dynamic model. To compensate for the unmodeled disturbances in the dynamic model, a robust nonlinear feedback term is included in the control law. One of the primary challenges in the control design is the presence of input-multiplicative parametric uncertainty in the dynamic model for the control actuator. This challenge is mitigated through innovative algebraic manipulation in the regulation error system derivation along with a Lyapunov-based adaptive control law. To address practical implementation considerations, where sensor measurements of the complete state are not available for feedback, a detailed analysis is provided to demonstrate that system observability can be ensured through judicious placement of pressure (and/or velocity) sensors. Based on this observability condition, a sliding-mode observer design is presented, which is shown to estimate the unmeasurable states using only the available sensor measurements. A detailed Lyapunov-based stability analysis is provided to prove that the proposed closed-loop active thermoacoustic control system achieves asymptotic (zero steady-state error) regulation of multiple thermoacoustic modes in the presence of the aforementioned model uncertainty. Numerical Monte Carlo-type simulation results are also provided, which demonstrate the performance of the proposed closed-loop control system under various sets of operating conditions.

1. Introduction

Rijke-type instability is a widely investigated example of a thermoacoustic phenomenon, which describes the generation of potentially unstable pressure oscillations that results from the dynamic coupling between unsteady heat transfer and acoustics [1, 2]. The resulting oscillations in Rijke-type systems can degrade performance and even cause structural damage in combustion systems. Based on this fact, thermoacoustic instability is a primary challenge that must be addressed in the design and manufacture of land-based gas turbines and aircraft engines [3–11]. Other applications for which thermoacoustic oscillations are a concern include boilers, furnaces, ramjet engines, and rocket motors. The myriad practical engineering applications impacted by Rijke-type instability necessitate the design of reliable control systems to regulate the potentially catastrophic effects of thermoacoustic oscillations.

Control design methods for thermoacoustic oscillation suppression systems can be separated into two main categories: passive control and active control approaches. Passive control methods [12–18] can employ acoustic dampers, such as Helmholtz resonators [13] or acoustic liners [12], or they can be achieved by physically redesigning the system by changing the location of the heat source, for example. Passive approaches offer the virtues of simplicity and inexpensive maintenance; however, the performance of passive control methods can only be ensured over a relatively narrow range of operating conditions [4]. To expand the usable range of operating conditions, active control methods offer the capability to automatically adjust the level of control actuation in response to sensor stimuli.

Active control methods are usually implemented in closed-loop configurations, where sensor measurements are utilized in a feedback loop to automatically drive the input signal to the actuators. Figure 1 provides an example of functional schematic of a closed-loop thermoacoustic oscillation suppression system. The two primary strategies for achieving closed-loop active control of thermoacoustic oscillations include (1) using a monopole-like acoustic source such as a loudspeaker to control the acoustic field [19] or (2) using a secondary fuel injector to control the unsteady heat release rate [20, 21]. Several active control approaches to suppress thermoacoustic oscillations have been presented in recent research literature.

Figure 1.

A functional implementation diagram of a thermoacoustic oscillation control system, including a microphone for sensing and a loudspeaker for actuation.

Standard linear control systems for thermoacoustic oscillation suppression are based on stabilizing the closed-loop system through causing the dominant eigenmodes to exponentially decay. However, for realistic thermoacoustic systems where the eigenmodes are nonorthogonal, controlling only the dominant eigenmode can result in the excitation of other modes as a result of the coupling between the acoustic modes. To address this challenge, a transient growth controller is presented in [22–24], which achieves strict dissipativity. Experimental or numerical empirical methods for thermoacoustic oscillation control have been widely considered, but more systematic approaches such as robust and adaptive control have gained popularity in more recent research. Active linear control methods have been widely investigated for applications considering simplified thermoacoustic dynamic models [1, 25]. However, by leveraging the tools of nonlinear control, effective suppression control of thermoacoustic oscillations can be achieved over a wide range of operating conditions and dynamic model uncertainty.

Physically speaking, thermoacoustic oscillation suppression can be achieved by disrupting the inherent dynamic coupling between the unsteady heat release and the acoustic waves. By designing active control systems to alter the interaction between the acoustic waves and the unsteady heat release, the amplitude of the thermoacoustic oscillations can be forced to decrease, instead of increase. Additional challenges in designing control systems for thermoacoustic oscillations can be incurred as a result of parametric model uncertainty and unmodeled operating conditions. The recent result in [26] presents a nonlinear active control method, which is proven to asymptotically regulate thermoacoustic oscillations in Rijke-type systems that do not include parametric model uncertainty and unmodeled nonlinearities. The design of active closed-loop control systems for thermoacoustic oscillation suppression that achieve reliable performance over a wide range of operating conditions and model uncertainty remains very much an open problem.

In this chapter, an observer-based nonlinear active closed-loop control method is presented, which achieves asymptotic suppression of self-excited thermoacoustic oscillations in a Rijke-type system, where the system dynamic model includes unmodeled nonlinearities and parametric uncertainty in the system dynamics and actuator dynamics. To achieve the result, a well-accepted thermoacoustic model is utilized, which employs arrays of sensors and monopole-like actuators. To facilitate the control design, the original dynamic equations are recast in a control-amenable form, which explicitly includes the effects of unmodeled, nonvanishing external disturbances and linear time delay. A sliding-mode observer-based nonlinear control law is then derived to regulate oscillations in the thermoacoustic system. A primary challenge in the control design is the presence of input-multiplicative parametric uncertainty in the control-oriented model. This challenge is handled through innovative algebraic manipulation in the regulation error system derivation along with a Lyapunov-based adaptive law. A rigorous Lyapunov-based stability analysis is used to prove that the closed-loop system achieves asymptotic regulation of a thermoacoustic system consisting of multiple modes. Numerical Monte Carlo-type simulation results are also provided, which demonstrate the performance of the proposed closed-loop active thermoacoustic oscillation suppression system.

2. Thermoacoustic system model

The thermoacoustic system model that will be utilized in this chapter consists of a horizontal Rijke tube with multiple actuators. The model is identical to that studied in our previous work in [22–24, 27].

Consider the system shown in Figure 2 , where the actuators are modeled as multiple monopole-like moving pistons. It will be assumed that K ≥ 1 actuators are available for control purposes. To facilitate the following observer and control design and analysis, a block diagram is also provided in Figure 3 .

Figure 2.

A control-oriented schematic of a combustion system with actuators modeled as monopole-like moving pistons.

Figure 3.

A block diagram illustrating the main components of the proposed robust and adaptive thermoacoustic oscillation control system.

To facilitate the subsequent model development, nondimensional system variables are defined as

$u=\frac{\tilde{u}}{u_0},p=\frac{\tilde{p}}{\gamma M_0{p}_0},{{\displaystyle \dot{Q}}}_s=\frac{{{\displaystyle \dot{{\displaystyle \tilde{Q}}}}}_s}{\gamma p_0{u}_0},$ E1

$x=\frac{\tilde{x}}{L_0},t=\frac{\tilde{t}{c}_0}{L_0},\frac{\delta \left(x-x_f\right)}{L_0}=\tilde{\delta }\left(\tilde{x}-{\tilde{x}}_f\right),$ E2

where the above tilde notation denotes the dimensional quantities and the subscript 0 denotes the mean values. In Eqs. (1) and (2), x ∈  $\mathbb{R}$ denotes the location along the duct (the actuators are located at x_ak  ∈  $\mathbb{R}$ , for k = 1 ,  …  , K, and the heat source is located at x_f ), t ∈  $\mathbb{R}$ _≥0 denotes nondimensional time, p(x, t) ∈  $\mathbb{R}$ is the acoustic pressure, u(x, t) ∈  $\mathbb{R}$ denotes the velocity, M ∈  $\mathbb{R}$ is the Mach number, c ∈  $\mathbb{R}$ is the speed of sound, L ₀ ∈  $\mathbb{R}$ is the length of the duct, and γ ∈  $\mathbb{R}$ is the ratio of specific heats.

By using the nondimensionalized variables defined in Eqs. (1) and (2), the thermoacoustic system with K actuators can be expressed as

$\frac{\partial u}{\partial t}+\frac{\partial p}{\partial x}=0,$ E3

$\frac{\partial p}{\partial t}+\mathit{\zeta p}+\frac{\partial u}{\partial x}=\left(\gamma -1\right){\dot{Q}}_s\delta \left(x-x_f\right)+\gamma \sum _{k=1}^K‍{\alpha }_{\mathit{ak}}{v}_{\mathit{ak}}\delta \left(x-x_{\mathit{ak}}\right),$ E4

In the expressions in Eqs. (3) and (4), ζ ∈  $\mathbb{R}$ denotes a damping coefficient; thus, the term ζp physically expresses losses resulting from the effects of friction and thermo-viscous damping, and α_ak  ∈  $\mathbb{R}$ represents a dimensionless area ratio that can be explicitly defined as α_ak  = S_ak /S for k = 1 ,  …  , K. The nondimensional heat release rate ${\dot{Q}}_s\in \mathbb{R}$ is explicitly defined as [19].

${\dot{Q}}_s=\mathcal{K}\left[\sqrt{\mid \frac{1}{3}+u_f\left(t-\tau \right)\mid }-\sqrt{\frac{1}{3}}\right],$ E5

where

$\mathcal{K}=\frac{2L_w\left(T_w-{\overline{T}}_0\right)}{\sqrt{3u_0}\mathit{S\gamma }{p}_0}\sqrt{\mathit{\pi \lambda }{c}_v{\rho }_0\frac{d_w}{2}}.$ E6

In Eq. (6), d_w , L_w , and T_w  ∈  $\mathbb{R}$ denote the diameter, length, and temperature of the heated wire, respectively; ρ ∈  $\mathbb{R}$ denotes air density; T ∈  $\mathbb{R}$ is temperature; λ is thermal conductivity; c_v  ∈  $\mathbb{R}$ denotes the specific heat capacity at constant volume; and τ ∈  $\mathbb{R}$ represents the time delay between the initial velocity field perturbation (i.e., the actuation) and the resulting effect on the heat release. Readers are referred to [22] for details on the numerical values of the physical parameters used for the thermoacoustic model being considered in this chapter. The gas is assumed to be inviscid, perfect, and nonconductive.

The acoustic pressure p and velocity u inside the duct can be expressed as a superposition of the duct natural modes as

$p\left(x,t\right)=-\sum _{j=1}^N‍\frac{sin\left(\mathit{j\pi x}\right)}{\mathit{j\pi }}{\eta }_j\left(t\right),$ E7

$u\left(x,t\right)=\sum _{j=1}^N‍cos\left(\mathit{j\pi x}\right){\eta }_j\left(t\right),$ E8

where N ∈ ℕ denotes the number of modes considered in the numerical discretization.

The actuation signal v_ak  ∈  $\mathbb{R}$ of the k^th monopole-like source (e.g., a loudspeaker) [28] can now be expressed as

$v_{\mathit{ak}}={\mathcal{R}}_{k}u\left(x_{\mathit{ak}}\right)+{\mathcal{S}}_{k}p\left(x_{\mathit{ak}}\right),$ E9

where $\mathcal{R}$ _k and ${\mathcal{S}}_k\in \mathbb{R}$ are dimensionless control parameters of the actuators. After using the expressions in Eqs. (4)–(8), the discretized governing equations are obtained as

$\frac{{\ddot{\eta }}_j}{\mathit{j\pi }}+\mathit{j\pi }{\eta }_j+{\zeta }_j\frac{{\dot{\eta }}_j}{\mathit{j\pi }}=-2\left(\gamma -1\right){\dot{Q}}_s\left(x_f,t-\tau \right)sin\left(\mathit{j\pi }{x}_f\right)-2\gamma \sum _{k=1}^K‍{\alpha }_{\mathit{ak}}{v}_a\left(x_{\mathit{ak}},t\right)sin\left(\mathit{j\pi }{x}_{\mathit{ak}}\right).$ E10

In Eq. (10), ζ represents the overall damping in the system [29, 30]. For the model under consideration here, the terms p and $\frac{\text{∂}u}{\text{∂}x}$ are taken to be zero at the ends of the duct. Moreover, it is assumed that no acoustic energy is dissipated in the thermal and viscous boundary layers at the duct walls [29–31].

3. Control-oriented model derivation

To facilitate the presentation of the main ideas, we consider a thermoacoustic system with two modes (i.e., N = 2). However, the theoretical development presented here can be directly extended to address N modes.

In the following discussion, the vector of modes (i.e., the state vector) will be annotated as η(t) ≜ [η ₁(t), η ₂(t)] ^T , and define

$\mathbf{\Psi }\left(x\right)={\left[cos\left(\mathit{\pi x}\right),cos\left(2\mathit{\pi x}\right)\right]}^T,$ E11

$\mathbf{\Phi }\left(x\right)={\left[sin\left(\mathit{\pi x}\right),sin\left(2\mathit{\pi x}\right)\right]}^T.$ E12

Assuming that ∣u_f  (t − τ)∣ < 1/3, the heat release rate can be approximated as

${\dot{Q}}_s\left(x_f,t-\tau \right)\approx \frac{\sqrt{3}\mathcal{K}}{2}{\mathbf{\Psi }}_f^T\eta \left(t-\tau \right),$ E13

where Ψ _f  ≜ Ψ(x_f ).

By following a derivation procedure similar to that presented in [24], the dynamics of the duct natural modes can be expressed as

$\mathbf{M}\ddot{\eta }+\left(\mathbf{D}-\left(\gamma -1\right)\mathcal{K}\sqrt{3}\tau {\mathbf{\Phi }}_f{\mathbf{\Psi }}_f^T\right)\dot{\eta }+\left({\mathbf{M}}^{-1}+\left(\gamma -1\right)\mathcal{K}\sqrt{3}{\mathbf{\Phi }}_f{\mathbf{\Psi }}_f^T\right)\eta -\mathbf{h}\left(\eta ,\dot{\eta }\right)-{\tau }_d=\mathcal{B}\mathbf{v},$ E14

where

$\mathbf{M}=diag\left\{\frac{1}{\pi },\frac{1}{2\pi }\right\},\mathbf{D}=diag\left\{\frac{{\zeta }_1}{\pi },\frac{{\zeta }_2}{2\pi }\right\}.$ E15

In Eq. (14), η(t) = [η ₁(t), η ₂(t)] ^T  ∈  $\mathbb{R}$ ² is a vector containing the natural modes, $\mathbf{h}\left(\eta ,\dot{\eta }\right)\in {\mathbb{R}}^2$ is an unknown nonlinear function, and τ_d (t) ∈  $\mathbb{R}$ ² is a general unknown bounded disturbance. To facilitate the control development in the following analysis, the dynamic equation in Eq. (14) is rewritten in the control-oriented form:

$\mathbf{M}\ddot{\eta }+{\mathbf{C}}_{\mathbf{d}}\dot{\eta }+\mathbf{K}\eta -\mathbf{h}\left(\eta ,\dot{\eta }\right)-{\tau }_d=\mathcal{B}\mathbf{v},$ E16

where the uncertain constant terms C _d and K ∈  $\mathbb{R}$ ^2 × 2 are defined as

${\mathbf{C}}_{\mathbf{d}}=\mathbf{D}-\left(\gamma -1\right)\mathcal{K}\sqrt{3}\tau {\mathbf{\Phi }}_f{\mathbf{\Psi }}_f^T,$ E17

$\mathbf{K}={\mathbf{M}}^{-1}+\left(\gamma -1\right)\mathcal{K}\sqrt{3}{\mathbf{\Phi }}_f{\mathbf{\Psi }}_f^T.$ E18

Also, in Eq. (16), the uncertain constant control input gain matrix $\mathcal{B}\in R^{2\times 2}$ is defined via the relationship:

$\mathcal{B}\mathbf{v}=-2\gamma \left[\sum _{k=1}^K‍{\alpha }_{\mathit{ak}}{\mathcal{R}}_k\mathbf{\Phi }\left(x_{\mathit{ak}}\right){\mathbf{\Psi }}^T\left(x_{\mathit{ak}}\right)\right]\eta +2\gamma \left[\sum _{k=1}^K‍{\alpha }_{\mathit{ak}}{\mathcal{S}}_k\mathbf{\Phi }\left(x_{\mathit{ak}}\right){\mathbf{\Phi }}^T\left(x_{\mathit{ak}}\right)\mathbf{M}\right]\dot{\eta }$ E19

where v(t) ∈  $\mathbb{R}$ ² is a subsequently defined auxiliary control signal.

Remark 1 Note that Eq. (19) highlights one of the primary challenges in the control design presented in this chapter. Specifically, the input-multiplicative parametric uncertainty in $\mathcal{B}$ presents a nontrivial control design challenge, which will be mitigated in the proposed control method through the use of a novel Lyapunov-based adaptive law.

Assumption 1 The unknown nonlinear disturbance τ_d (t) satisfies

$‖{\tau }_d\left(t\right)‖\le \zeta ,\forall t\ge 0,$ E20

where ζ ∈  $\mathbb{R}$ denotes a positive bounding parameter.

4. Control development

In this section, a rigorous regulation error system development will be utilized to develop a nonlinear control system, which will be proven to effectively compensate for the inherent parametric uncertainty in the dynamic model of the thermoacoustic system in addition to the uncertain actuator model. Moreover, the proposed controller compensates for unmodeled, norm-bounded disturbances present in the dynamic model (e.g., the disturbances could represent unmodeled nonlinearities resulting from time delays due to the finite heat release rate).

4.1. Open-loop error system

The robust and adaptive nonlinear control design presented here is motivated by the desire to eliminate the transient growth of acoustical energy in a thermoacoustic dynamic system. To present the control design methodology, we consider a simplified N = 2 mode scenario, which will be shown to regulate the modes η ₁(t) and η ₂(t) to zero in the sense that

$‖\eta \left(t\right)‖,‖\dot{\eta }\left(t\right)‖\to 0.$ E21

To mathematically describe the regulation control objective, an auxiliary regulation error signal r(t) ∈  $\mathbb{R}$ ² is defined as

$\mathbf{r}=\dot{\eta }+\mathit{\alpha \eta },$ E22

where α ∈  $\mathbb{R}$ denotes a positive, constant control gain. After taking the time derivative of Eq. (22), multiplying the result by M, and using Eq. (16), the regulation error dynamics can be expressed as

$\mathbf{M}\dot{\mathbf{r}}=-{\mathbf{D}}_{\mathbf{d}}\dot{\eta }-\mathbf{K}\eta +\mathbf{h}\left(\eta ,\dot{\eta }\right)+\mathcal{B}\mathbf{v}+{\tau }_d.$ E23

To address the case where the constant matrices D _d  = C _d  + M α, K, and $\mathcal{B}$ are uncertain, the dynamics can be linearly parameterized as

${\mathbf{Y}}_1{\theta }_1=-{\mathbf{D}}_{\mathbf{d}}\dot{\eta }-\mathbf{K}\eta ,$ E24

${\mathbf{Y}}_2{\theta }_2=\mathcal{B}\mathbf{v}.$ E25

In Eqs. (24) and (25), ${\mathbf{Y}}_1\left(\eta ,\dot{\eta }\right)\in {\mathbb{R}}^{2\times p_1}$ and Y ₂(v) ∈  $\mathbb{R}$ ^{2 × p ₂} are measurable regression matrices, and θ ₁ ∈  $\mathbb{R}$ ^{p ₁} and θ ₂ ∈  $\mathbb{R}$ ^{p ₂} are vectors containing the uncertain constant parameters in D _d , K, and $\mathcal{B}$ . The constants p ₁ and p ₂ ∈ ℕ denote the number of uncertain parameters in the vectors θ ₁ and θ ₂, respectively.

To facilitate the subsequent Lyapunov-based adaptive control law development to compensate for the input-multiplicative uncertain matrix $\mathcal{B}$ , an estimate ${\widehat{\theta }}_2\left(t\right)\in {\mathbb{R}}^{p_2}$ of the uncertain vector θ ₂ is defined via the linear parameterization:

${\mathbf{Y}}_2{\widehat{\theta }}_2=\widehat{\mathcal{B}}\mathbf{v}.$ E26

In Eq. (26), $\widehat{\mathcal{B}}\left(t\right)\in {\mathrm{ℝ}}^{2\times 2}$ denotes a time-varying estimate of the uncertain constant matrix $\mathcal{B}$ . By adding and subtracting the term $\widehat{\mathcal{B}}\left(t\right)\mathbf{v}\left(t\right)$ in Eq. (23) and using Eqs. (24) and (26), the open-loop error dynamics can be expressed as

$\mathbf{M}\dot{\mathbf{r}}={\mathbf{Y}}_1{\theta }_1+\mathbf{h}\left(\eta ,\dot{\eta }\right)+{\mathbf{Y}}_2{\tilde{\theta }}_2+\widehat{\mathcal{B}}\mathbf{v}+{\tau }_d,$ E27

where ${\tilde{\theta }}_2\left(t\right)\in {\mathbb{R}}^{p_2}$ denotes the parameter estimate mismatch, which is defined as

${\tilde{\theta }}_2\triangleq {\theta }_2-{\widehat{\theta }}_2.$ E28

The error dynamics in Eq. (27) are now in a form amenable for the design of a robust and adaptive control law, which compensates for the parametric uncertainty and unmodeled nonlinearities present in the system dynamics.

Assumption 2 The unknown nonlinear term $\mathbf{h}\left(\eta ,\dot{\eta }\right)$ can be upper bounded as

$‖\mathbf{h}\left(\eta ,\dot{\eta }\right)‖\le \rho \left(‖\mathbf{z}‖\right)‖\mathbf{z}‖,$ E29

where ρ(·) ∈  $\mathbb{R}$ is a positive, globally invertible nondecreasing function and z(t) ∈  $\mathbb{R}$ ⁴ is defined as

$\mathbf{z}\left(t\right)={\left[{\eta }^T\left(t\right){\mathbf{r}}^T\left(t\right)\right]}^T.$ E30

In Eq. (29), ‖·‖ denotes the standard Euclidean norm of the vector argument.

Assumption 2 is mild in the sense that inequality (29) is satisfied for a wide range of nonlinear function $\mathbf{h}\left(\eta ,\dot{\eta }\right)$ .

4.2. Closed-loop error system

Based on the open-loop error system in Eq. (27), the control input v(t) is designed as

$\mathbf{v}={\widehat{\mathcal{B}}}^{-1}\left(-{\mathbf{Y}}_1{\widehat{\theta }}_1-\left(k_s+1\right)\mathbf{r}-\beta sgn\left(\mathbf{r}\right)-\eta \right),$ E31

where k_s  ∈  $\mathbb{R}$ denotes a positive, constant control gain and β ∈  $\mathbb{R}$ ^2 × 2 is a positive-definite, diagonal control gain matrix. In Eq. (31), sgn(·) denotes a vector form of the standard signum function.

After substituting the control input expression in Eq. (34) into the open-loop dynamics in Eq. (27), the closed-loop error system is obtained as

$\mathbf{M}\dot{\mathbf{r}}={\mathbf{Y}}_1{\tilde{\theta }}_1+{\mathbf{Y}}_2{\tilde{\theta }}_2+\mathbf{h}\left(\eta ,\dot{\eta }\right)-\eta -\left(k_s+1\right)\mathbf{r}-\beta sgn\left(\mathbf{r}\right)+{\tau }_d,$ E32

where ${\tilde{\theta }}_1\left(t\right)\in {\mathbb{R}}^{p_1}$ is the parameter estimate mismatch defined as

${\tilde{\theta }}_1\triangleq {\theta }_1-{\widehat{\theta }}_1.$ E33

Based on Eq. (32) and the subsequent stability analysis, the parameter estimates ${\widehat{\theta }}_1\left(t\right)$ and ${\widehat{\theta }}_2\left(t\right)$ are generated online according to the adaptive laws:

${{\displaystyle \dot{{\displaystyle \widehat{\theta }}}}}_1=\text{proj}\left({\mathrm{\Gamma }}_1{\mathbf{Y}}^T\mathbf{r}\right),{{\displaystyle \dot{{\displaystyle \widehat{\theta }}}}}_2=\text{proj}\left({\mathrm{\Gamma }}_2{\mathbf{Y}}_2^T\mathbf{r}\right),$ E34

where Γ₁ ∈  $\mathbb{R}$ ^{p ₁ ×  p ₁} and Γ₂ ∈  $\mathbb{R}$ ^{p ₂ ×  p ₂} are positive-definite adaptation gains.

Remark 2 The function proj(·) in Eq. (34) denotes a normal projection algorithm, which ensures that the following inequalities are satisfied:

${\underset{¯}{\theta }}_1\le {\widehat{\theta }}_1\le {\overline{\theta }}_1,{\underset{¯}{\theta }}_2\le {\widehat{\theta }}_2\le {\overline{\theta }}_2,$ E35

where ${\underset{¯}{\theta }}_1,{\overline{\theta }}_1,{\underset{¯}{\theta }}_2\mathrm{and}{\overline{\theta }}_2\in \mathbb{R}$ represent known, constant lower and upper bounds of the elements of ${\widehat{\theta }}_1\left(t\right)$ and ${\widehat{\theta }}_2\left(t\right)$ , respectively. In the current result, the use of the proj(·) function is primarily motivated by the desire to avoid singularities in the matrix estimate and facilitate the matrix inverse calculation in Eq. (31).

To facilitate the following stability analysis, the control gain matrix β will be selected to satisfy the sufficient condition:

${\lambda }_{\min }\left\{\beta \right\}>\zeta ,$ E36

where ζ is introduced in Eq. (20) and λ _min{·} denotes the minimum eigenvalue of the argument.

5. Stability analysis

Theorem 1 The control law in Eq. (31) with adaptive laws defined as in Eq. (34) ensures asymptotic regulation of the thermoacoustic modes η ₁(t) and η ₂(t) in the sense that

$‖\eta \left(t\right)‖\to 0\mathrm{as}t\to \infty$ E37

provided that k_s is selected as sufficiently large (see the subsequent stability proof) and β is selected to satisfy inequality (36).

Proof. Let $V\left(\eta ,\mathbf{r},{\widehat{\theta }}_1,{\widehat{\theta }}_2,t\right)\in \mathbb{R}$ be defined as the nonnegative function:

$V\left(t\right)\triangleq \frac{1}{2}{\eta }^T\eta +\frac{1}{2}{\mathbf{r}}^T\mathbf{Mr}+\frac{1}{2}{\tilde{\theta }}_1^T{\mathrm{\Gamma }}_1^{-1}{\tilde{\theta }}_1+\frac{1}{2}{\tilde{\theta }}_2^T{\mathrm{\Gamma }}_2^{-1}{\tilde{\theta }}_2.$ E38

After taking the time derivative of Eq. (38) and using Eq. (32), $\dot{V}\left(t\right)$ can be expressed as

$\dot{V}\left(t\right)={\eta }^T\left(\mathbf{r}-\mathit{\alpha \eta }\right)+{\mathbf{r}}^T\left(\mathbf{h}\left(\eta ,\dot{\eta }\right)-\eta -\left(k_s+1\right)\mathbf{r}-\beta sgn\left(\mathbf{r}\right)+{\tau }_d\right)+{\mathbf{r}}^T\left({\mathbf{Y}}_1{\tilde{\theta }}_1+{\mathbf{Y}}_2{\tilde{\theta }}_2\right)-{\tilde{\theta }}_1^T{\mathrm{\Gamma }}_1^{-1}{{\displaystyle \dot{{\displaystyle \widehat{\theta }}}}}_1-{\tilde{\theta }}_2^T{\mathrm{\Gamma }}_2^{-1}{{\displaystyle \dot{{\displaystyle \widehat{\theta }}}}}_2,$ E39

where Eq. (22) was utilized. After substituting the adaptive laws in Eq. (34) and canceling common terms, $\dot{V}\left(t\right)$ can be expressed as

$\dot{V}\left(t\right)=-\alpha {\eta }^T\eta +{\mathbf{r}}^T\left(\mathbf{h}\left(\eta ,\dot{\eta }\right)-\left(k_s+1\right)\mathbf{r}-\beta sgn\left(\mathbf{r}\right)+{\tau }_d\right).$ E40

By using inequalities of Eqs. (20) and (29), the expression in Eq. (40) can be upper bounded as

$\dot{V}\left(t\right)\le -\alpha {‖\eta ‖}^2-\left(k_s{‖\mathbf{r}‖}^2-\rho \left(‖\mathbf{z}‖\right)‖\mathbf{z}‖‖\mathbf{r}‖\right)-{‖\mathbf{r}‖}^2-\beta {\mathbf{r}}^{T}sgn\left(\mathbf{r}\right)+\zeta ‖\mathbf{r}‖.$ E41

After completing the squares for the parenthetic terms in Eq. (41), the upper bound on $\dot{V}\left(t\right)$ can be expressed as

$\dot{V}\left(t\right)\le -\alpha {‖\eta ‖}^2-{‖\mathbf{r}‖}^2-\beta \left|\mathbf{r}\right|+\zeta ‖\mathbf{r}‖-k_s{\left(‖\mathbf{r}‖-\frac{\rho \left(‖\mathbf{z}‖\right)}{2k_s}‖\mathbf{z}‖\right)}^2+\frac{{\rho }^2\left(‖\mathbf{z}‖\right)}{4k_s}{‖\mathbf{z}‖}^2,$ E42

where the fact that r ^T  sgn(r) = |r| was utilized. After using inequality (36), the upper bound in Eq. (42) can be expressed as

$\dot{V}\left(t\right)\le -\left({\lambda }_0-\frac{{\rho }^2\left(‖\mathbf{z}‖\right)}{4k_s}\right){‖\mathbf{z}‖}^2,$ E43

where λ ₀ ≜ min {α, 1} and the triangle inequality (i.e., |r| ≥ ‖r‖ ∀r ∈  $\mathbb{R}$ ⁿ ) was utilized. Based on Eq. (43), $\dot{V}\le -c{‖\mathbf{z}‖}^2$ , for some positive constant c, inside the set $\mathcal{R}$ , where $\mathcal{R}$ is defined as

$\mathcal{R}\triangleq \left\{\mathbf{z}|\mathbf{z}<{\rho }^{-1}\left(2\sqrt{{\lambda }_0{k}_s}\right)\right\}.$ E44

The expressions in Eqs. (38) and (43) can be used to prove that η(t), r(t), ${\tilde{\theta }}_1\left(t\right)$ , and ${\tilde{\theta }}_2\left(t\right)\in {\mathcal{L}}_{\infty }$ in $\mathcal{R}$ and Eq. (22) can then be used to prove that $\dot{\eta }\left(t\right)\in {\mathcal{L}}_{\infty }$ in $\mathcal{R}$ . Given that η(t) and r(t) ∈  $\mathcal{L}$ _∞ in $\mathcal{R}$ , Eq. (31) can be used along with Eq. (35) to prove that v(t) ∈  $\mathcal{L}$ _∞ in $\mathcal{R}$ . Since η(t), $\dot{\eta }\left(t\right)$ , and v(t) ∈  $\mathcal{L}$ _∞ in $\mathcal{R}$ , ${\mathbf{Y}}_1\left(\eta ,\dot{\eta }\right)$ and Y ₂(v) ∈  $\mathcal{L}$ _∞ in $\mathcal{R}$ . Given that η(t), r(t), $\dot{\eta }\left(t\right)$ , ${\mathbf{Y}}_1\left(\eta ,\dot{\eta }\right)$ , and Y ₂(v) ∈  $\mathcal{L}$ _∞ in $\mathcal{R}$ , Eq. (32) can be used along with Eq. (35) to show that $\dot{\mathbf{r}}\left(t\right)\in {\mathcal{L}}_{\infty }$ in $\mathcal{R}$ . Since $\dot{\eta }\left(t\right)\mathrm{and}$ $\dot{\mathbf{r}}\left(t\right)\in {\mathcal{L}}_{\infty }$ in $\mathcal{R}$ , η(t) and r(t) are uniformly continuous in $\mathcal{R}$ . It then follows from Eq. (30) that z(t) is uniformly continuous in $\mathcal{R}$ . Given that η(t), r(t), ${\tilde{\theta }}_1\left(t\right)$ , and ${\tilde{\theta }}_2\left(t\right)\in {\mathcal{L}}_{\infty }$ in $\mathcal{R}$ , V(t) ∈  $\mathcal{L}$ _∞ in $\mathcal{R}$ , and Eq. (43) can be integrated to prove that ${\int }_0^{\infty }‍{‖\mathbf{z}\left(t\right)‖}^2\mathit{dt}\in {\mathcal{L}}_{\infty }$ in $\mathcal{R}$ . Thus, z(t) ∈  $\mathcal{L}$ _∞ ∩  $\mathcal{L}$ ₂ in $\mathcal{R}$ . Barbalat’s lemma can now be invoked to prove that ‖z(t)‖→0 as t→∞. Hence, ‖η(t)‖→0 as t→∞ in $\mathcal{R}$ , where the set $\mathcal{R}$ can be made arbitrarily large by increasing the control gain k_s —a semiglobal result.

6. Sliding-mode observer design

In practical thermoacoustic systems, the full state of the dynamic system is not directly measurable, and so it must be estimated through direct sensor measurements of velocity and pressure. This section presents an observer design, which is utilized to estimate the complete state of the system. The necessary observability condition can easily be satisfied through judicious sensor placement.

Let $\mathbf{x}={\left[{\eta }^T,{\dot{\eta }}^T\mathbf{M}\right]}^T$ denote the state. For simplicity in the subsequent development, only one sensor is assumed to be available at location x_s . Then, the system can be rewritten as

$\dot{\mathbf{x}}\left(t\right)=\mathbf{Ax}\left(t\right)+\mathbf{Bx}\left(t-\tau \right)+\mathbf{Gu}\left(t\right),$ E45

$y=\mathbf{Cx},$ E46

where

$\mathbf{A}=\left[\begin{array}{ll}\mathbf{0}& {\mathbf{M}}^{-1}\\ -{\mathbf{M}}^{-1}& -\mathbf{D}\end{array}\right],\mathbf{B}=\left[\begin{array}{ll}\mathbf{0}& \mathbf{0}\\ -\mathbf{W}& \mathbf{0}\end{array}\right],\mathbf{G}=\left[\begin{array}{l}\mathbf{0}\\ \mathbf{I}\end{array}\right],$ E47

where

$\mathbf{W}=\sqrt{3}\left(\gamma -1\right)\mathcal{K}{\mathbf{\Phi }}_f{\mathbf{\Psi }}_f^T,$ E48

and the output matrix C is determined by the sensor choice and its location. The output equation for the velocity sensor case is given by

$y=u\left(x_s,t\right)={\mathbf{\Psi }}_s^T\eta ,$ E49

and thus

$\mathbf{C}=\left[\begin{array}{ll}{\mathbf{\Psi }}_s^T& {\mathbf{0}}_{1\times 2}\end{array}\right].$ E50

For the pressure sensor case, we have

$y=p\left(x_s,t\right)=-{\mathbf{\Phi }}_s^T\dot{\eta },$ E51

and hence

$\mathbf{C}=\left[\begin{array}{ll}{\mathbf{0}}_{1\times 2}& -{\mathbf{\Phi }}_s^T\mathbf{M}.\end{array}\right].$ E52

It is assumed that the sensor location x_s is chosen such that the systems (45) and (46) subject to τ = 0 are observable, i.e.,

$rank\left[{\mathbf{C}}^T{\mathbf{E}}^T{\mathbf{C}}^T{\left({\mathbf{E}}^T\right)}^2{\mathbf{C}}^T{\left({\mathbf{E}}^T\right)}^3{\mathbf{C}}^T\right]=4,$ E53

where E = A + B.

It can be shown that there exists a coordinate transformation of the forms z = [z ₁, z ₂] ^T  = Tx, z ₁ ∈  $\mathbb{R}$ , and z ₂ ∈  $\mathbb{R}$ ³, such that the systems (45) and (46) in the new variables take the following form:

$\dot{\mathbf{z}}\left(t\right)=\tilde{\mathbf{A}}\mathbf{z}\left(t\right)+\tilde{\mathbf{B}}\mathbf{z}\left(t-\tau \right)+\tilde{\mathbf{G}}\mathbf{u}\left(t\right),$ E54

$y\left(t\right)=z_1\left(t\right),$ E55

where $\tilde{\mathbf{A}}=\mathbf{TA}{\mathbf{T}}^{-1}$ , $\tilde{\mathbf{B}}=\mathbf{TB}{\mathbf{T}}^{-1}$ , and $\tilde{\mathbf{G}}=\mathbf{TG}$ .

The estimate $\widehat{\mathbf{z}}$ of the the state z will now be generated via the observer equation

$\dot{\widehat{\mathbf{z}}}\left(t\right)=\tilde{\mathbf{A}}\widehat{\mathbf{z}}\left(t\right)+\tilde{\mathbf{B}}\widehat{\mathbf{z}}\left(t-\tau \right)+\tilde{\mathbf{G}}\mathbf{u}\left(t\right)+\mathbf{L}sgn\left(y\left(t\right)-{\widehat{z}}_1\left(t\right)\right).$ E56

Then, the error dynamics $\overline{\mathbf{z}}=\mathbf{z}-\widehat{\mathbf{z}}$ can be written as

$\dot{\overline{\mathbf{z}}}\left(t\right)=\tilde{\mathbf{A}}\overline{\mathbf{z}}\left(t\right)+\tilde{\mathbf{B}}\overline{\mathbf{z}}\left(t-\tau \right)-\mathbf{L}sgn\left({\overline{z}}_1\left(t\right)\right).$ E57

Partition the system above as

${\dot{\overline{z}}}_1\left(t\right)={\tilde{A}}_{11}{\overline{z}}_1\left(t\right)+{\tilde{\mathbf{A}}}_{12}{\overline{\mathbf{z}}}_2\left(t\right)+{\tilde{B}}_{11}{\overline{z}}_1\left(t-\tau \right)+{\tilde{\mathbf{B}}}_{12}{\overline{\mathbf{z}}}_2\left(t-\tau \right)-L_{1}sgn\left({\overline{z}}_1\left(t\right)\right),$ E58

${\dot{\overline{\mathbf{z}}}}_2\left(t\right)={\tilde{\mathbf{A}}}_{21}{\overline{z}}_1\left(t\right)+{\tilde{\mathbf{A}}}_{22}{\overline{\mathbf{z}}}_2\left(t\right)+{\tilde{\mathbf{B}}}_{21}{\overline{z}}_1\left(t-\tau \right)+{\tilde{\mathbf{B}}}_{22}{\overline{\mathbf{z}}}_2\left(t-\tau \right)-{\mathbf{L}}_{2}sgn\left({\overline{z}}_1\left(t\right)\right).$ E59

The following result can now be stated:

Theorem 2 Let ${\tilde{\mathbf{E}}}_{12}={\tilde{\mathbf{A}}}_{12}+{\tilde{\mathbf{B}}}_{12}$ and ${\tilde{\mathbf{E}}}_{22}={\tilde{\mathbf{A}}}_{22}+{\tilde{\mathbf{B}}}_{22}$ , and assume that the pair $\left[{\tilde{\mathbf{E}}}_{22},{\tilde{\mathbf{E}}}_{12}\right]$ is observable. Then, the observer gain L ₁ can be chosen such that the system (57) is asymptotically stable for all τ ∈ [0, τ_max ] with τ_max defined by:

${\tau }_{\max }=\frac{1}{2\sqrt{{\lambda }_{\max }\left({\mathbf{Q}}_1^{-T}{\mathbf{F}}^T\mathbf{PH}{\mathbf{Q}}^{-1}{\mathbf{H}}^T\mathbf{PF}{\mathbf{Q}}_1^{-1}\right)}},$ E60

where $\mathbf{F}={\tilde{\mathbf{E}}}_{22}-{\mathbf{L}}_2{\tilde{\mathbf{E}}}_{12}/L_1$ , $\mathbf{H}={\tilde{\mathbf{B}}}_{22}-{\mathbf{L}}_2{\tilde{\mathbf{B}}}_{12}/L_1$ , and Q is any symmetric, positive-definite matrix such that P is a symmetric, positive-definite matrix P solution of the Lyapunov matrix equation:

${\mathbf{F}}^T\mathbf{P}+\mathbf{PF}=-\mathbf{Q},$ E61

and Q ₁ is the square root of the matrix Q, i.e.,

${\mathbf{Q}}_1^T{\mathbf{Q}}_1=\mathbf{Q}.$ E62

Proof. From Eq. (58), sliding mode exists in an area:

$L_1>\mid {\tilde{A}}_{11}{\overline{z}}_1\left(t\right)+{\tilde{\mathbf{A}}}_{12}{\overline{\mathbf{z}}}_2\left(t\right)+{\tilde{B}}_{11}{\overline{z}}_1\left(t-\tau \right)+{\tilde{\mathbf{B}}}_{12}{\overline{\mathbf{z}}}_2\left(t-\tau \right)\mid .$ E63

Condition (63) guarantees sliding in Eq. (58) along the manifold ${\overline{z}}_1=0$ ; thus $\widehat{\eta }\to \eta$ . According to the equivalent control method [32], the system in sliding mode behaves as if $L_{1}sgn\left({\overline{z}}_1\right)$ is replaced by its equivalent value ${\left(L_{1}sgn\left({\overline{z}}_1\right)\right)}_{\mathrm{eq}}$ which can be calculated from subsystem (58) assuming ${\overline{z}}_1=0$ and ${\dot{\overline{z}}}_1=0$ . Hence,

${\left(L_{1}sgn\left({\overline{z}}_1\right)\right)}_{\mathrm{eq}}={\tilde{\mathbf{A}}}_{12}{\overline{\mathbf{z}}}_2\left(t\right)+{\tilde{\mathbf{B}}}_{12}{\overline{\mathbf{z}}}_2\left(t-\tau \right).$ E64

Substitution of Eq. (64) into Eq. (59) yields

${\dot{\overline{\mathbf{z}}}}_2=\left({\tilde{\mathbf{A}}}_{22}-{\mathbf{L}}_2{\tilde{\mathbf{A}}}_{12}/L_1\right){\overline{\mathbf{z}}}_2\left(t\right)+\left({\tilde{\mathbf{B}}}_{22}-{\mathbf{L}}_2{\tilde{\mathbf{B}}}_{12}/L_1\right){\overline{\mathbf{z}}}_2\left(t-\tau \right).$ E65

Using the fact that the pair $\left[{\tilde{\mathbf{E}}}_{22},{\tilde{\mathbf{E}}}_{12}/L_1\right]$ is observable, the observer gain L ₂ can be chosen such that the eigenvalues of the matrices ${\tilde{\mathbf{E}}}_{22}-{\mathbf{L}}_2{\tilde{\mathbf{E}}}_{12}/L_1$ have negative real parts. Thus, the subsystem (65) is asymptotically stable for τ ≤ τ _max. This implies that $\overline{\mathbf{z}}\to$ 0 and $\widehat{\mathbf{z}}\left(t\right)=\mathbf{z}\left(t\right)$ , and hence, $\widehat{\mathbf{x}}=\mathbf{x}$ .

7. Simulation results

A numerical simulation was created for two modes and two actuators (i.e., N = K = 2) to demonstrate the performance of the control law described by Eqs. (26), (31), and (34). The simulation utilizes the dynamics described in Eq. (16) for a system with two modes. The physical parameters used in the simulation are given in Table 1 .

ρ	1.025 kg/m^3	λ	0.0328 W/m K
cv	719 J/kg K	γ	1.4
L 0	1 m	Lw	2.5 m
c	344 m/s	u 0	0.3 m/s
T 0	295 K	Tw	1680 K
dw	0.5 × 10^−3 m	S	1.56 × 10^−3 m
P 0	8.69 × 10^4 Pa	α a 1  = α a 2	0.01
ζ 1	0.0440	ζ 2	0.1657

Table 1.

Physical parameters.

The initial conditions for the modes were selected as η ₁(0)  =  0.07 and η ₂(0) = 0. The initial conditions of the velocity and pressure parameters at the location of the heat source are u(x_f , 0) = − 0.1646 and p(x_f , 0) = 0 for x_f = 0.7. Figure 4 shows the open-loop velocity and pressure perturbations at the heat source in the absence of control actuation.

Figure 4.

Time response of the velocity u(t) and pressure p(t) at the heat source location during open-loop (uncontrolled) operation.

In closed-loop operation, the adaptation gain matrices used in the simulation were selected as Γ ₁ = 0.1I ₈ and Γ ₂ = 0.001I ₄. The value for the gain α was selected as 1, and the control gain matrices β and k _s were selected as β = diag {2.43, 2.1} ,  k _s  = diag {4.5, 2.4}. The results of 20 Monte Carlo-type simulations for the closed-loop operation are shown in Figures 5 – 10 . Figure 5 shows the time evolution of the velocity and pressure values at the heat source location during closed-loop controller operation. Figure 6 shows the time history of the modes η ₁(t) and η ₂(t) during closed-loop operation. The results clearly show the capability of the proposed robust and adaptive control law to drive the states to zero. The commanded control signals are shown in Figure 7 . The control actuation remains within the reasonable limits throughout closed-loop operation. Figures 8 – 10 show the time responses of the elements of the parameter estimate vectors ${\widehat{\theta }}_1\left(t\right)$ and ${\widehat{\theta }}_2\left(t\right)$ during closed-loop controller operation.

Figure 5.

Time response of the velocity u(t) and pressure p(t) during closed-loop controller operation at the heat source location.

Figure 6.

Time response of the oscillation modes η ₁(t) and η ₂(t) during closed-loop controller operation.

Figure 7.

Commanded control inputs ν ₁(t) and ν ₂(t) during closed-loop controller operation.

Figure 8.

Time response of the adaptive parameter estimates ${\widehat{\theta }}_{11}\left(t\right)$ , ${\widehat{\theta }}_{12}\left(t\right)$ , ${\widehat{\theta }}_{13}\left(t\right)$ , and ${\widehat{\theta }}_{14}\left(t\right)$ during closed-loop controller operation.

Figure 9.

Time response of the adaptive parameter estimates ${\widehat{\theta }}_{15}\left(t\right)$ , ${\widehat{\theta }}_{16}\left(t\right)$ , ${\widehat{\theta }}_{17}\left(t\right)$ , and ${\widehat{\theta }}_{18}\left(t\right)$ during closed-loop controller operation.

Figure 10.

Time response of the adaptive parameter estimates ${\widehat{\theta }}_{21}\left(t\right)$ , ${\widehat{\theta }}_{22}\left(t\right)$ , ${\widehat{\theta }}_{23}\left(t\right)$ , and ${\widehat{\theta }}_{24}\left(t\right)$ during closed-loop controller operation.

8. Conclusion

A robust and adaptive nonlinear control method is presented, which asymptotically regulates thermoacoustic oscillations in a Rijke-type system in the presence of dynamic model uncertainty and unknown disturbances. To demonstrate the methodology, a well-accepted thermoacoustic dynamic model is introduced, which includes arrays of sensors and monopole-like actuators. To facilitate the derivation of the adaptive control law, the dynamic model is recast as a set of nonlinear ordinary differential equations, which are amenable to control design. To compensate for the unmodeled disturbances in the dynamic model, a robust nonlinear feedback term is included in the control law. One of the primary challenges in the control design is the presence of input-multiplicative parametric uncertainty in the dynamic model for the control actuator. This challenge is mitigated through innovative algebraic manipulation in the regulation error system derivation along with a Lyapunov-based adaptive control law. To address practical implementation considerations, where sensor measurements of the complete state are not available for feedback, a detailed analysis is provided to demonstrate that system observability can be ensured through judicious placement of pressure (and/or velocity) sensors. A sliding-mode observer design is developed, which is shown to estimate the unmeasurable states using only the available sensor measurements. A detailed Lyapunov-based stability analysis is provided to prove that the proposed closed-loop active thermoacoustic control system achieves asymptotic (zero steady-state error) regulation of multiple thermoacoustic modes in the presence of the aforementioned model uncertainty. Numerical Monte Carlo-type simulation results are also provided, which demonstrate the performance of the proposed closed-loop control system under 20 different sets of operating conditions.