Feedback control of combustion instabilities from within limit cycle oscillations using H∞ loop-shaping and the ν-gap metric

Combustion instabilities arise owing to a two-way coupling between acoustic waves and unsteady heat release. Oscillation amplitudes successively grow, until nonlinear effects cause saturation into limit cycle oscillations. Feedback control, in which an actuator modifies some combustor input in response to a sensor measurement, can suppress combustion instabilities. Linear feedback controllers are typically designed, using linear combustor models. However, when activated from within limit cycle, the linear model is invalid, and such controllers are not guaranteed to stabilize. This work develops a feedback control strategy guaranteed to stabilize from within limit cycle oscillations. A low-order model of a simple combustor, exhibiting the essential features of more complex systems, is presented. Linear plane acoustic wave modelling is combined with a weakly nonlinear describing function for the flame. The latter is determined numerically using a level set approach. Its implication is that the open-loop transfer function (OLTF) needed for controller design varies with oscillation level. The difference between the mean and the rest of the OLTFs is characterized using the ν-gap metric, providing the minimum required ‘robustness margin’ for an H∞ loop-shaping controller. Such controllers are designed and achieve stability both for linear fluctuations and from within limit cycle oscillations.


Introduction
It is desirable to operate modern industrial gas turbines and aeroengines under lean premixed combustion conditions in order to reduce NO x emissions. However, lean premixed systems are highly susceptible to combustion instabilities arising from the coupling between the heat release rate perturbations and the acoustic disturbances within the combustion chamber [1]. These instabilities are nearly always undesirable as they lead to large oscillation amplitudes and damaging structural vibrations of the combustion chamber [2].
The stability of a combustion chamber is determined by the balance between the energy gained from the flame/acoustic interactions and various dissipation processes [1,3,4]. Active feedback control can be used to interrupt the coupling between the acoustic waves and unsteady heat release and prevent or suppress instability. The design of most types of controller requires prior knowledge of how the combustor responds to actuation-the 'open-loop transfer function' (OLTF) between monitored thermodynamic properties within the combustion chamber (typically one or more pressure measurements) and the actuator signal(s) [4].
Although linear models may be accurate at low perturbation levels, the dynamics of real unstable combustors become dominated by nonlinear mechanisms once amplitudes grow sufficiently. In many combustors, including gas turbine and aeroengine combustors, it is notable that even during large amplitude oscillations the behaviour of the acoustic waves remains linear, with nonlinearity arising almost entirely from the way in which the flame's unsteady heat release responds to oncoming flow disturbances [5][6][7]. (Acoustic waves in rocket engines are, however, typically nonlinear [8].) This flame nonlinearity is 'weak' such that a 'describing function' provides a good model of the nonlinear flame response. This assumes that the flow forcing amplitude level is a quasi-steady parameter, and that the dominant frequency of the unsteady heat release rate response matches the incoming flow forcing frequency, but with a gain and phase shift that depend on the forcing level as well as the frequency. This gives rise to a 'family' of frequency responses that depend on forcing level. It has been shown that this weak nonlinearity can cause saturation into limit cycle either via saturation of the heat release rate amplitude or via a change in phase between the heat release rate and acoustic pressure as the modulation level increases [6,9]. This may cause the dominant unstable frequency to change, which, in turn, alters the OLTF and complicates the controller design. Although some recent work has identified nonlinear states other than limit cycle oscillations as resulting from combustion instability [10][11][12], the focus of this work will be limit cycle oscillations, which are by far the most common reported nonlinear state.
The design of linear feedback controllers is typically based on the OLTF corresponding to small, linear disturbances. Such feedback controllers are then activated from within nonlinear limit cycle oscillations, during which the linear OLTF will not be valid. It is typically observed that such controllers successfully stabilize oscillations despite this paradox [13][14][15]. However, their effectiveness is not guaranteed, and if they fail then there has, until now, been no framework for systematically improving their design.
In this work, the concept of a flame describing function (FDF) is combined with the ν-gap metric [16] in order to design controllers which are guaranteed to achieve stabilization across all oscillation levels, from small linear fluctuations to large limit cycle oscillations. In modern robust control systems theory, the concept of a ν-gap between open-loop plants represents the most useful measure of distance between systems when one is concerned specifically with applying feedback to those systems. The ν-gap metric fits very naturally into H ∞ loop-shaping controller design, providing a bound on the minimum required 'robustness stability margin' for an H ∞ loop-shaping controller [16,17]. Modern linear robust control, of which H ∞ loop-shaping is one methodology, has been successfully used in diverse fields, such as control of combustion instabilities [18][19][20], control of developing flow [21], voltage converter design [22] and vehicle oscillation damper optimization [23], to achieve control even in nonlinear systems with limit cycles. H ∞ loop-shaping combines classical loop-shaping, which can obtain tradeoffs between good performance and robust stability, with modern H ∞ optimization in order to guarantee closed-loop stability and a level of robust stability at all frequencies [16,17]. It can be applied to plants with high dimensions and is computationally inexpensive [19], at least for plants of low and moderate order. The resulting controllers generally have orders comparable to those of the plants.
It is typical to approximate the OLTF as a high-order rational transfer function (RTF) for the purposes of controller design, meaning that high-order controllers will be generated. These are less reliable to implement than low-order controllers. As a final step, this work proposes a low-order controller design method, which combines H ∞ loop-shaping, the ν-gap metric and a low-order fitting procedure.
The remainder of the paper is organized as follows. Section 2 presents the ducted laminar flame combustor model used as a test case in the paper, including the coupling of linear modelling for the acoustic waves with the nonlinear flame response. Section 3 demonstrates the determination of the nonlinear FDF from simulations using a level set approach (LSA), ensuring that our flame model is grounded in flow physics. Section 4 analyses the resulting nonlinear thermoacoustic instabilities for different flame positions. Section 5 presents the design of linear controllers which are guaranteed to be stabilizing from within limit cycle oscillations, using H ∞ loop-shaping and the ν-gap metric. This includes the design of standard highorder controllers and a method for obtaining low-order controllers. The viscothermal damping of acoustic waves in ducts is presented in appendix, and conclusions are drawn in the final section.

Combustor model
A geometrically simple model combustor is considered for this study. The geometrical simplicity permits us to embed the advanced modelling features that capture the important phenomena at play in real combustors. This includes linear acoustic wave behaviour with realistic boundary and viscothermal losses, and complex nonlinear flame kinematics for sufficiently high amplitude perturbations.
The model considered is a ducted laminar flame combustor, often known as a Rijke tube. It is shown schematically in figure 1, comprising a cylindrical Rijke tube with both ends open to the atmosphere. Denoting the distance along the duct axis by the vector x, the tube inlet and outlet are at x = x 0 = 0 and x = x 2 = l, respectively. The inner diameter of the tube is D, and a laminar premixed methane-air conical flame is located at x = x 1 = x f . The Bunsen-type burner has an axisymmetric geometry with outlet diameter d. A loudspeaker is located at the bottom of the tube as an actuator for control, with a pressure sensor located a distance x = x r downstream of the combustion zone. The modelling analysis makes the following assumptions: (i) the combustion zone is considered 'compact' compared with the acoustic wavelength, and we limit attention to low frequencies where we need to account only for longitudinal waves [24]. (ii) The fluids before and after the flame are assumed to be perfect gases [24]. (iii) There is no noise produced by the entropy waves formed during the unstable combustion process-these waves are assumed to leave the tube without interaction with the flow at the end of tube [25]. (iv) The flame remains stabilized at the burner outlet. Flame intrinsic instabilities [26] and irregular response to strong disturbances [27] are not accounted for.
The flow can be taken to be composed of a steady uniform mean flow (denoted¯) and small perturbations (denoted ). Harmonic time variations are considered for which all fluctuating variables have the form a =â e iωt , where ω is complex angular frequency. The regions upstream and downstream of the flame are indicated by the subscripts 1 and 2, respectively. Then, the pressure, p, and longitudinal velocity, u, in the two regions can be expressed in terms of the amplitudes A + n and A − n of the upward and downward propagating acoustic waves, respectively, their wavenumbers, k ± =k ± + k ± , wherek ± = ω/(c ±ū) and k ± are wavenumber corrections for the viscothermal damping effects (see appendix), and the speed of sound, c, and density, ρ [28] The acoustic wave strengths, A + n and A − n , are related by pressure reflection coefficients 1 at the tube boundaries-the inlet and outlet pressure reflection coefficients are characterized by R 1 and R 2 , respectively.
For an open duct, radiation of sound from the duct ends can be accounted for using frequency-dependent reflection coefficients that act as low-pass filters [29]. The tube also appears acoustically longer than its physical length, and this end correction can be accounted for. The reflection coefficient for an unflanged pipe end can be approximately expressed in the Laplace domain as [29,30] where τ d = D/c. It was confirmed in [31] that the gain of this reflection coefficient decreases with both duct diameter and frequency. The reflection coefficients at the inlet and outlet of the Rijke tube are denoted by R 1 and R 2 , respectively. In order to relate the acoustic wave strengths across the flame, the linearized mass, momentum and energy conservation equations across the flame zone (x = x f ) are combined with the perfect gas equation. Denoting the jump across the flame as [·] 2 1 , the amplitude of the flame's heat release rate perturbation asQ and the cross-sectional area of the duct as S, this yields two equations which are independent of the downstream density, ρ 2 , which contains an entropy component [32]  linking the acoustic wave strengths before and after the flame are obtained in terms of s, the Laplace transform variable [33] (where s n = s + Λ n √ s, n = 1, 2, the second term accounting for viscothermal damping with the coefficient Λ n defined in equation (A 2)), time delays τ ± n = (x n − x n−1 )/(c n ±ū n ) and Ξ = (S 1 /S 2 )(ρ 2c2 )/(ρ 1c1 ). The governing equation is expressed as ⎡ The matrix coefficients, M n , are given in appendix. In order to close equation (2.5) that contains five unknowns (A ± n andQ), a flame model relatingQ to the acoustic fluctuations is needed. In this work, we use an FDF, whose form is derived numerically in §3.
Although the combustor configuration contains only two acoustic modules (the two regions upstream and downstream of the flame), the system is complex owing to the frequencydependent pressure reflection coefficients, viscothermal damping terms and the presence of flame nonlinearity. A numerical solver is used to solve the nonlinear eigenproblem [34]. The solver is part of OSCILOS, an open source combustion instability low-order network simulation tool. It represents complicated combustor geometries as a network of simple connected modules, modelling the acoustic wave behaviour analytically using linear wave-based methods, and able to incorporate a wide range of linear and nonlinear flame models. It can predict thermoacoustic modal frequencies, growth rates, mode shapes and the time evolution of disturbances, and has been validated against experiments [35].
In this work, a relatively low cost numerical method, compared with large eddy simulations [35,36], is used to determine the FDF of the laminar conical premixed flame. This numerical method is described in §3.

Determining the flame describing function
In the model combustor, a laminar conical premixed flame within the duct generates a mean and fluctuating heat release rate. The heat release rate is taken to be directly proportional to the flame's surface area, which means that an accurate flame model for insertion into equation (2.5) can be found by capturing the kinematic response of the flame to oncoming acoustic disturbances.
The instantaneous location of a moving flame front subjected to oncoming flow perturbations can be captured using the LSA proposed by Markstein [37]. This kinematic model is also known as the G-equation model, and treats the flame as an extremely thin interface (corresponding to the iscontour G = 0), separating fresh reacting flow (denoted by G < 0) from the burned gases (G > 0). The flame front is assumed to propagate normal to itself with the velocity U · n − s L , where U is the fresh reacting flow velocity vector and in this work is two-dimensional, U = (u x , u r ), owing to the axisymmetric configuration. s L is the flame burning speed and n = VG/|VG| is the normal vector. The transport equation for G is [38] ∂G ∂t The mutual interactions of the flame front kinematics and the hydrodynamic flow field are not accounted for [39], meaning the fresh mixture flow field above the burner outlet and inside the flame can be treated using a simple model which is independent of the flame evolution. It has been shown that when the G-equation model is combined with a simple incompressible velocity perturbation model in the fresh mixture, complex flame front evolution phenomena, including the 'pinch-off' that occurs for strong flow perturbations, can be captured [40,41]. This model was recently validated using direct numerical simulations [42], and is expressed as The composition of the fresh reacting flow is considered constant, and the dependence of the flame burning speed on the flame front curvature is taken to be s L = s L0 (1 − Lκ), where s L0 is the burning speed of a laminar planar flame, L is the Markstein length, which depends on mixture equivalence ratio φ and thermal conditions [43] (and hence is constant in this work), whereas κ = V · n is the local flame curvature [37,43].
The G-equation (equation (3.1)) combined with the incompressible flow velocity model (equation (3.2)) and the flame curvature effect is solved numerically using a fifth-order weighted essentially non-oscillatory (WENO) scheme [44] for spatial discretization and a thirdorder total variation diminishing (TVD) Runge-Kutta scheme for time integration [45]. The spatial and temporal resolutions are fixed in all calculations: x = r = 5 × 10 −3 d and t = 5 × 10 −4 d/ū 1 , respectively, where d is the diameter of the burner outlet. To significantly reduce the computational cost, a local LSA is employed, only accounting for the grids around the flame front [46]. The boundary condition for the centreline of the flame is symmetry-flame axisymmetry is exploited to simulate only half of the flame. The flame base is considered always attached to the burner lip, with G at the flame base reinitialized to 0 at every simulation time step. No boundary condition is needed for other boundaries, because the local LSA is implemented.
For premixed flames, the heat release rate can be calculated usinġ whereρ 1 denotes the mean density of the fresh reacting flow,ω T indicates the heat release rate per volume from combustion and A indicates the whole space of the LSA simulation. δ(G) is the delta function, representing integration over the flame front, which is performed using a high-order scheme [47].
The weakly nonlinear flame model to be extracted from these simulations is in the form of an 'FDF', expressed as whereQ 1 (û 1 /ū 1 , s) is the Laplace transform of the fundamental term of heat release rate perturbations. Note that this can be thought of as an 'input amplitude dependent' transfer function or frequency response. To validate the LSA, a flame for which experimental data are available is firstly simulated. A laminar conical methane-air premixed flame with burner outlet diameter d = 20 mm, equivalence ratio φ = 1.0, mean velocityū 1 = 1.5 m s −1 , laminar flame speed s L0 = 0.368 m s −1 and Markstein length L = 1 mm [43] is considered. Figure 2 compares the simulated flame front location from the LSA with a four-colour Schlieren image from experiments, when the perturbation frequency and level are f p = 100 Hz andû 1 /ū 1 = 0.1, respectively. The flame front is seen to respond strongly even for a relatively weak perturbation level. The prediction matches the experimental result with the flame cusps accurately captured. The evolution of the flame front throughout an entire forcing cycle for modulation at f p = 50 Hz andû 1 /ū 1 = 0.15 is shown in figure 3a. At t/T = 1 2 , where T = 1/f p , pinch-off occurs and a flame pocket is produced. The frequency spectrum of the heat release rate perturbations is shown for 50 Hz forcing in figure 3b. Although the nonlinear flame response causes harmonics at high frequencies, the fact that the dominant flame response frequency exactly matches the forcing frequency validates the 'weakly' nonlinear assumption of the FDF.
The FDF, F , calculated using the LSA is compared with experimental results for a similar laminar conical premixed flame [50], with almost the same burner outlet diameter. The simulation and experimentally measured FDFs are shown in figure 4. The FDF has the shape of lowpass filter at constantû 1 /ū 1 , with nonlinearity evident at higher frequencies-the gain of FDF   decreases with increasingû 1 /ū 1 , but with the response linear at low frequencies even for a largê u 1 /ū 1 = 0.214. The phase dependence on forcing amplitude begins at a higher frequency than the gain dependence. The gain and phase lag of the FDF both generally decrease with modulation level, which is consistent with the experimental results. It can clearly be concluded that the LSA captures both the flame front motion and the main features of the FDF of laminar conical premixed flames.
Simulation is then performed for the combustor test-case flame. For a laminar conical premixed flame, the cut-off frequency of the FDF is inversely proportional to the diameter of the burner outlet [6,51,52]. To guarantee that a sufficiently large heat release feeds the unstable thermoacoustic mode, the diameter of the burner outlet is reduced to d = 6 mm (meaning that the cut-off frequency is approximately three times that of the FDF shown in figure 4), with the other flow parameters unchanged. Simulations are performed for modulation frequencies ranging from 10 to 500 Hz in intervals of 10 Hz, and with normalized velocity perturbation levels ranging from 0.025 to 0.4 in intervals of 0.025. The FDF is shown in figure 5. The shape remains similar, with the response generally falling off with both frequency and modulation level. Although |F| does not strictly decrease with increasingû 1 /ū 1 above 300 Hz, the accuracy of the flow velocity perturbation model (equation (3.2)) breaks down at high frequencies. With the dominant unstable frequency in the Rijke tube being 210 Hz, the FDF results can be considered accurate below 300 Hz. The FDF provided by the LSA simulations (denoted F LSA ) is discrete in both frequency and forcing amplitude and is limited to the frequency range [10, 500] Hz. The FDF can be considered         where Q =T 2 /T 1 − 1 andT denotes the mean flow temperature. The thermoacoustic modes are given by the values of s = iω = σ + iω i at constantû 1 /ū 1 for which det(M 0 + M 1 ) = 0, where det denotes the matrix determinant, ω i = 2π f is real angular frequency and σ is the growth rate. The stability of the mode is defined by the sign of the growth rate, with a positive value corresponding to an unstable mode. This calculation is performed within OSCILOS. The combustor parameters considered for the analysis are listed in table 1.

Simulation and analysis of the unstable combustor
It is found that only the mode with the eigenfrequency f ≈ 210 Hz has the potential for instability (this corresponds to the cold flow fundamental frequency of 173 Hz for perfect acoustic reflections from the boundaries). This is consistent with both the FDF gain fall-off and the acoustic losses increasing with frequency, having a stabilizing effect on higher-order modes. A growth rate map for this unstable mode as a function of both flame position within the duct and flame velocity perturbation level is shown in figure 6, with only positive (unstable) growth rates shown. The most dangerous unstable position is a quarter of the way along the duct, which is generally the  situation in a Rijke tube [53]. For each flame location, limit cycle oscillations will be established at the velocity fluctuation level corresponding to the uppermost contour, as this is where the growth rate has dropped to zero. The flame location, x f = 0.2 m, is chosen for feedback controller design, corresponding to the limit cycle being established at a relatively largeû 1 /ū 1 . Figure 7 shows the trajectories of the growth rate and eigenfrequency with increasingû 1 /ū 1 for this flame location. For weak perturbations, the growth rate is strongly positive, implying a rapidly growing exponential envelope, e.g. forû 1 /ū 1 = 0.025, the growth rate is 46 s −1 , implying exponentially growth with an envelope of exp (46t). This growth rate decreases withû 1 /ū 1 , becoming zero whenû 1 /ū 1 = 0.18, hence this is the velocity fluctuation amplitude at which we expect limit cycle oscillations to occur.

Feedback control from within limit cycle oscillations (a) Open-loop transfer function
Actuators for feedback control of combustion instabilities generally fall within two categories. The first seeks to modify the fuel-air composition, normally via a fuel valve, whereas the second uses acoustic actuators to modify the pressure waves inside the combustor directly [1,18]. This work assumes the latter-these are generally better suited to light duty combustors. A loudspeaker at the bottom end of the duct is assumed, injecting a pressure signal p L . A microphone mounted x r = 0.5 m downstream of the flame measures the pressure at that location. To simplify the analysis, we do not account for the instrument transfer functions of the loudspeaker and microphone [13,54]: we consider that the electrical driving signal to the loudspeaker I in equals p L , and the microphone electrical output signal I out equals the measured pressure signal p r . Furthermore, we assume that the fuel-air ratio is unaffected by the acoustic waves, maintaining the applicability of the FDF calculated in §3. For model-based feedback controller design (figure 8), it is necessary to have a model for the OLTF from p L to p r :  Each modal peak is associated with a phase change of π , implying a second-order transfer function contribution. The direction of the phase change provides stability information-a phase increase in π corresponds to an unstable conjugate pair of poles, whereas a phase decrease in π indicates a stable conjugate pair. From figure 9, normalized velocity perturbation amplitudes, u 1 /ū 1 , of 0-0.15 correspond to unstable OLTFs, whereas amplitudes of 0.2 and 0.25 correspond to stable OLTFs. This is consistent with the previous finding that limit cycle oscillations occur at u 1 /ū 1 = 0.18.
It is thus clear that a controller design based on just one of these linear OLTFs cannot guarantee the entire set of plants. In order to guarantee stability from within limit cycle oscillations, a robust controller must -reduce the size of perturbations when the system is in limit cycle, with the effect of reducingû 1 /ū 1  It involves designing a pre-or post-compensator to 'shape' the singular values (i.e. the gain for single-input-singleoutput (SISO) systems) of the open-loop system, so that they have a high gain where disturbance rejection is important, and a low gain where robustness to multiplicative plant uncertainty is required (see references [17,55] for an explanation of different types of uncertainty from a control systems perspective). A stabilizing controller is then synthesized using Matlab (e.g. using the ncfsyn command), which minimizes the H ∞ norm of the closed-loop transfer function matrix (see equation (5.5))-this norm also provides a measure of the system robustness. If it is suitably small, then a robust stabilizing controller has been found. The gain of the shaped transfer function will then not be significantly affected by the controller, but the phase will be altered so as to achieve stability [17]. If the norm is not suitably small then another iteration on the choice of the pre-or post-compensator must be performed. The underpinning mathematics can be found in [17]. The shape of the open-loop plant P is first modified using pre-and post-compensators W 1 and W 2 ( figure 8). For an SISO system, modifying either one of these has the same effect. The H ∞ controller K ∞ is then synthesized in Matlab, with the final feedback controller K, being K = The ν-gap metric provides a useful measure of the distance in the state space between two open-loop plants with respect to how they behave when connected in a unity feedback control loop [16]. For SISO systems, it can be computed from the frequency response: the ν-gap, δ ν , between plants P 1 and P 2 is given by δ ν (P 1 , P 2 ) = sup ω ψ(P 1 (iω), P 2 (iω)), provided the winding number condition is satisfied 4 The ν-gap metric fits very naturally into H ∞ loop-shaping control design, which returns both a feedback controller design, and the corresponding 'stability margin', b ∞ , for that controller/plant pairing. The stability margin b ∞ for a shaped plant P W = W 1 P 1 W 2 and the resulted controller K ∞ is represented as [17] where · ∞ is the infinity norm. If the controller is applied to a different plant (e.g. P 2 ), separated from the initial by ν-gap δ ν , as long as b ∞ > δ ν , then the H ∞ loop-shaping controller is also guaranteed to stabilize this second system.

(c) H ∞ loop-shaping controller design
H ∞ loop-shaping controller design will need to be designed based on a nominal linear plant. The differences between this nominal plant and all of the possible OLTFs will then be characterized using the ν-gap metric. The controller design can be schematically explained by the flowchart shown in figure 10. For combustion instabilities, nonlinearities are introduced by the FDF, which can be treated as a discrete set of FTFs for different velocity perturbation levelsû 1 /ū 1 . From figure 6, limit cycle oscillations are generally established forû 1 /ū 1 ≤ 0.2. We thus calculate the mean FTF, F 0 , by averaging over the FTFs for velocity perturbation levelsû 1 /ū 1 ≤ 0.25. The mean FTF is fitted to a state space form in order to construct the mean OLTF (indicated by P 0 ) for feedback controller design. The presence of time delays in the wave-based models (e.g. see equations (4.1), and (5.2)) means the system is infinite dimensional. In order to implement H ∞ synthesis, the mean OLTF needs to be approximated by a RTF which will be finite dimensional. The exponential terms could be replaced by Padé approximants [56], but the presence of complicated damping mechanisms means that this is not straightforward. In this work, a frequency-domain fitting procedure using the Matlab command fitfrd is applied to obtain a RTF approximation to the mean OLTF. As the FDF is only valid over the [10,500] Hz frequency range, we fit the OLTF over a restricted frequency range [0, 1000] Hz, and also add a weighted low-pass filter to the RTF to reduce its gain at high frequencies. The resulting RTF, represented by P 0F , typically has an order of 16. A comparison of the original infinite-dimensional and RTF-fitted plants is shown in figure 11a.
As P 0F is SISO, we let the post-compensator W 2 = I and design the pre-compensator W 1 . A first-order low-pass filter with a cut-off frequency 300 Hz is chosen as W 1 , to drop the gain at high frequencies to provide robustness to multiplicative plant uncertainty. The ν-gaps between the plant P 0F and the plants P for differentû 1 /ū 1 can be identified from the frequency variations of ψ (see equation (5.4)), shown in figure 12a. Large bumps occur at each mode, with the deviation reaching a maximum at the unstable frequency-210 Hz. The maximum value of the ν-gap across all differentû 1 /ū 1 levels is max(δ ν ) = 0.18.    A H ∞ loop-shaping synthesis is then performed, using the Matlab command ncfsyn in the Robust control toolbox. The designed controller has an order of 16, and its stability margin is b ∞ = 0.41, which is significantly larger than max(δ ν ) above. This suggests that the controller should guarantee closed loop stability of all possible plants.
The performance of the negative feedback controller can be assessed by calculating the poles of the closed-loop transfer function P(s)/(1 + K(s)P(s)), where P(s) is the original OLTF (equation (5.1)). Figure 13a shows the poles locations for differentû 1 /ū 1 in the presence of the feedback controller. The growth rates of all modes are now stable, confirming closed loop stability for allû 1 /ū 1 .
Of course, for sufficiently large values ofû 1 /ū 1 , the combustor exhibits limit cycle oscillations and the linear OLTF appears stable even without feedback control. For these conditions, the sensitivity transfer function provides an essential measure of the controller performance. reduced. Thus, assuming that the limit cycle frequency matches or is close to the frequency of the unstable mode, f , the requirement now becomes |S(i2π f )| < 1 5 . Figure 14a shows the plots of |S| versus frequency for differentû 1 /ū 1 . The unstable frequency is from 205 to 215 Hz in the amplitude range 0 ≤û 1 /ū 1 ≤ 0.2. The fact that |S(i2π f )| < 1 at these frequencies means that the current controller will reduce the limit cycle oscillation amplitude. For sufficient attenuation, the plant will appear unstable (see the yellow and orange colour lines in figure 9) and the feedback controller will then stabilize the system.  It has been shown that H ∞ loop-shaping combined with the ν-gap metric can be used to design a robust controller that guarantees stabilization from within the limit cycle. The designed controller typically has an order comparable to that of the shaped plant [55]. High-order controllers are difficult to implement and exhibit lower reliability than low-order controllers. Although controller order reduction methods, such as balanced truncation methods [57], can be used on the plant or the controller, this work proposes an alternative approach to directly reduce the order of plant, provided that the maximum ν-gap (max(δ ν )) is smaller than the stability margin b ∞ .
The low-order method swaps the order of the shaping and fitting procedures. The plant P 0 is firstly shaped with W 1 using an 'aggressive' low-pass filter, which does not significantly alter the gain at low frequencies near the unstable mode, but reduces the gain at frequencies above the unstable mode, to provide robustness to fitting errors at these frequencies. A low-order fitting procedure is then implemented on the shaped plant P 0W , with the fit over low frequencies near the unstable mode prioritized. A fourth-order low-pass filter is chosen as the pre-compensator W 1 to shape the original plant, given by where ω L = 2π × 200 s −1 and ξ L = √ 2/2. W 1 has a gain of 0.5 at 200 Hz, and so does not significantly alter the gain at the unstable frequency. The gain rapidly falls off to 0.025 at 500 Hz.
The controller designed by H ∞ loop-shaping synthesis generally has a dimension comparable to that of the shaped plant, and so a fitted plant with low dimension is desirable. The shaped system contains only one unstable mode and exhibits fast gain fall-off at higher frequencies.
The fitting procedure is performed on the shaped plant P 0W = P 0 W 1 -this allows the dimension of the fitted RTF to be significantly reduced and yields a RTF with an order of 4. Figure 11b compares the original, P 0W , and fitted, P 0WF , plants. The fitted RTF captures the shape of original plant at low frequencies and has the shape of low-pass filter at high frequencies. Figure 12b shows the evolution of ψ(P 0WF (iω), P W (iω)) with frequency. Compared with the previous approach (figure 12a), there are larger deviations at high frequencies ( f ≥ 250 Hz) where the low-order fitting is less accurate. However, owing to the shape of W 1 , relatively small values of ψ(P 0WF (iω), P W (iω)) are guaranteed for all frequencies and the maximum of ν-gap for allû 1 /ū 1 is max(δ ν ) = 0.34.
A controller with an order of 3 is then obtained from the H ∞ loop-shaping synthesis. The final controller K = W 1 K ∞ has an order of 7 in the denominator and 3 in the numerator of its RTF, which is straightforward to implement. The stability margin is b ∞ = 0.52, which is larger than max (δ ν ).
The locations of the closed-loop poles for differentû 1 /ū 1 are shown in figure 13b. The unstable mode is stabilized for allû 1 /ū 1 , and the modes at high frequencies are not disturbed. However, the growth rate reduction for the unstable mode is less than for the high-order controller (figure 13a). This is due to the larger ν-gap caused by the low-order fitting. The gain of the sensitivity function S(s) with frequency is shown in figure 14b and is seen to be smaller than unity near the unstable mode frequency. This, combined with the fact that b ∞ > max (δ ν ), means that the low-order controller is also guaranteed to be stabilizing from within limit cycle oscillations.
The robustness of the low-order controller design to changes in operating condition is investigated by moving the flame position along the combustor. Figure 15 shows that the growth rate of the dominant mode becomes negative for all unstable configurations. Thus, the controller designed stabilizes all possible flame locations within this combustor.
The controller implementation in time domain simulations for the design condition of x f = 0.2 m is shown in figure 16 (more information on the time domain simulation method is found in [31,34,58]). The normalized velocity perturbations before the flame, u 1 /ū 1 , and the actuation signal from the loudspeaker p L are shown both before feedback control is activated and after it has been switched on. Prior to control, small disturbances are seen to grow rapidly, with a limit cycle established approximately beyond t = 0.4 s. When the controller is switched on at t = 0.5 s,  the velocity perturbations decrease rapidly to the weak background white noise. The loudspeaker signal, p L , is initially large in order to attain control and decreases to a much smaller value in order to maintain control.

Conclusion
This article has presented a method of designing linear feedback controllers that are guaranteed to stabilize combustion instability, whenever they are activated. The guarantees are subject to some very reasonable assumptions on the nature of the flame nonlinearity that would be satisfied very widely in practice. A single feedback controller design then can guarantee stability whether control is activated within the regime of small linear disturbances and exponential growth or within the regime of limit cycle oscillations in which the amplitude has saturated. The method combines the FDF approach for weakly nonlinear flame responses with H ∞ loopshaping for controller design and the ν-gap metric for characterizing nonlinearity as deviations from the 'average' linear model. The FDF gives rise to a 'family' of linear OLTFs corresponding to different flame perturbation levels. The maximum ν-gap between these and the average or nominal plant transfer function sets the minimum robustness stability margin required of the controller. An H ∞ loop-shaping controller is then designed which has a stability margin exceeding this, and which also attenuates fluctuations within limit cycle.
The control strategy has been demonstrated on a representative model of a nonlinear combustion system comprising an unsteady laminar flame inside a Rijke tube. Linear plane waves within the duct are assumed to suffer damping owing to end radiation and viscothermal wall damping mechanisms. The nonlinear kinematic response of the flame to flow disturbances is modelled using an LSA, which combines the G-equation kinematic flame model with a twodimensional flow velocity perturbation model and a flame front curvature-dependent burning