About this Journal Submit a Manuscript Table of Contents
Journal of Combustion
Volume 2012 (2012), Article ID 927345, 11 pages
Research Article

Robust Output Model Predictive Control of an Unstable Rijke Tube

Institute of Automatic Control, RWTH Aachen University, 52074 Aachen, Germany

Received 6 January 2012; Accepted 12 March 2012

Academic Editor: Xue-Song Bai

Copyright © 2012 Fabian Jarmolowitz et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


This work investigates the active control of an unstable Rijke tube using robust output model predictive control (RMPC). As internal model a polytopic linear system with constraints is assumed to account for uncertainties. For guaranteed stability, a linear state feedback controller is designed using linear matrix inequalities and used within a feedback formulation of the model predictive controller. For state estimation a robust gain-scheduled observer is developed. It is shown that the proposed RMPC ensures robust stability under constraints over the considered operating range.

1. Introduction

Modern gas turbines have to comply with increasingly stringent emission requirements for 𝑁𝑂𝑥. One of the most effective ways reducing these emissions is the development of lean premixed (LP) combustor systems. For gas turbines operated with natural gas thermal 𝑁𝑂𝑥 is the most relevant source, and its formation is highly temperature dependent. Both, a lean mixture and a premixing of fuel and oxidizer, which increases the mixture homogeneity and therefore avoids temperature peaks in the flame, reduce the combustion temperature and consequently the thermal 𝑁𝑂𝑥 emissions. Alongside the aforementioned advantages LP systems are more susceptible to combustion oscillations than conventional burner, because the heat release of lean-premixed flames is very sensitive to flow disturbances. Fluctuating heat release leads to fluctuating gas expansion. The gas expansion in the combustion zone acts as acoustic source and the emitted acoustic waves in turn influence the flame if reflected. Due to this coupling between acoustics of the combustion chamber and the heat release of the flame a feedback path is established which can give rise to thermoacoustic instabilities (see Figure 1). A well-known criterion for thermoacoustic instability is the Rayleigh-criterion. Originally, Rayleigh defined the criterion as 𝑇0𝑝(𝑡)̇𝑞(𝑡)𝑑𝑡>0 [1], but for example in [2] it is shown, that the acoustic losses 𝐿 have to be considered for checking instability. This leads to the original Rayleigh-criterion determines only if acoustic energy is fed into the system by the combustion process as energy source: 𝑇0𝑝(𝑡)̇𝑞(𝑡)𝑑𝑡>𝐿,(1) where 𝑝(𝑡) is the acoustic pressure, ̇𝑞 the fluctuating rate of heat release, and 𝐿 the acoustic losses. If the integral over one period of oscillation 𝑇 is positive, the heat release and acoustic pressure interfere constructively meaning that the heat release is in phase (up to 90) with the pressure waves. Thus, the flame feeds energy into the system. If this energy feed is higher than the acoustic losses 𝐿, the system is unstable and moves after an exponential growth into a limit cycle. The limit cycle is clearly a phenomenon of a nonlinear system and possibly a result of saturation effects in the heat release. It represents a stable trajectory, where losses and energy feed are in balance.

Figure 1: Active control of thermoacoustic instabilities (in parts taken from [5]).

As a consequence high-pressure oscillations take place in the combustor which are detrimental for performance, emissions and durability of the combustor components. To avoid these drawbacks, two main directions are possible, namely, passive and active control, to stabilize the thermoacoustic system [3, 4]. Passive control techniques include geometry modifications and integration of additional acoustic dampers like Helmholtz resonators. However, the operational range at which stabilization is achieved is limited. To overcome this limitation, active control can be applied to enlarge the region of stable operations. In addition, existing burners can be equipped with active control systems as a retrofit. This can be necessary not only for older systems but also for systems, where computational fluid dynamics (CFD) design tools failed to predict instabilities during development.

In this paper active control is investigated, that is, the closed-loop control of the thermoacoustic system. Closed-loop control uses a feedback path consisting of sensors, a control law and actuators to control a system, cf. Figure 1. The main goal for active control is to stabilize the system robustly possibly under constraints over a wide operating range. Robust stabilization means stabilization in the presence of uncertain system parameters. The increased demand for variable operations of gas turbines especially in power plants stems from the increased use of renewable energy sources. As a consequence, gas-fired power plants are operated more flexible to compensate the fluctuating input of renewable energy sources. The connected problem is two-fold. On the one hand, it may be the case that the optimal operating point in terms of efficiency, and emissions might be unstable itself and thus cannot be operated at. On the other hand, a change of the operating point can be delayed because unstable conditions have to be avoided during transition. Another demand for an active control system is to handle multiple inputs and multiple outputs (MIMO), since gas turbines are equipped with different types of actuators and sensors at different locations [3, 4], cf. Figure 1. Different sensors are presented in the literature like pressure transducers and chemiluminescence imaging techniques. For actuation different types of fuel modulation and acoustic forcing are applied. A MIMO controller has the advantage of using the combined information of all available sensors and the actuator operation of possibly different actuators is adapted to the actual machine state. In case multiple actuators are used, it is very likely that one actuator on its own has not the control authority to stabilize the system. Consequently, constraints on the control input have to be considered. In addition, the MIMO setup has another important benefit. The controller can be designed to be fault-tolerant, that is, if there exists some redundancy in the sensor and/or actuator setup the controller can maintain operation if some sensor and/or actuator fails. This fault-tolerance can be improved further by adapting the input constraints of the failed actuator.

In summary, the following demands have to be fulfilled by an active control system: (i)robust stabilization under constraints; (ii)handling MIMO; (iii) fault-tolerance.

To cope with these demands, robust output model predictive control (RMPC) is investigated in this work. A common physical modelling approach is applied to the Rijke tube, and from this a simplistic linear polytopic model is derived to model the system dynamics with (assumed) parameter uncertainties. The model is used to determine a robust linear state feedback controller by the use of linear matrix inequalities (LMI). This controller is incorporated in an approach presented in [6] to formulate a robust MPC and combined with a robust gain-scheduled observer. The proposed RMPC can have all the proposed characteristics needed for active control systems applied in combustion systems. In this work, however, the solution for a single input and single output (SISO) setup is shown, thus focusing on robust stabilization under constraints.

The paper is organised as follows. Section 2 presents the Rijke tube setup, the corresponding physical modelling, and its validation. Section 3 introduces the proposed robust output MPC, and Section 4 gives some control results for, varied tube length to demonstrate the robustness.

2. Modelling of the Rijke Tube

2.1. Setup

A schematic drawing of the Rijke tube setup is shown in Figure 2 with corresponding dimensions in Table 1. The test rig consists of a vertical glass tube that contains a smaller tube which can be drawn out vertically in order to vary the total length from the lower to upper end in a range of 1.09m up to 1.365m. This change of parameter is used to validate the robustness of the proposed active controller. The flame holder in the lower part of the tube is connected to a fuel feed line and holds a diffusion flame. A microphone is positioned a small distance beneath the flame holder. The test rig features both types of adjustment control, the variation of the fuel mass flow rate via valves at the feed line and antinoise. For the last mentioned purpose, a loudspeaker is positioned underneath the tube, which is used as actuator in this study.

Table 1: Dimensions of the Rijke tube setup.
Figure 2: Schematic drawing of the Rijke tube setup.
2.2. Acoustics

The acoustic model is a one-dimensional acoustic network consisting of simple geometric components, which are analytically tractable. The underlying assumption, which is common in thermoacoustic community [7, 8] is that the mean flow is homogeneous, the length to diameter ratio of the considered elements is sufficiently large such that only axial waves are relevant, the acoustic disturbances are linear, and that the flow is isentropic. Then, it is sufficient to consider only one-dimensional linear acoustics and the corresponding acoustic state vector (four-pole) is of order 𝑛=2. These assumptions are fulfilled for a Rijke tube [9]. The state can be defined, for example, by acoustic velocity 𝑢 and acoustic pressure 𝑝. The total pressure and velocity are defined as the sum of mean and acoustic value as 𝑝(𝑥,𝑡)=𝑝+𝑝(𝑥,𝑡), and 𝑢(𝑥,𝑡)=𝑢+𝑢(𝑥,𝑡). The solution of the governing conservation equations of mass, momentum, and energy [8] leads to the well-known wave equation whose solution under the aforementioned assumptions is 𝑝(𝑥,𝑡)=𝑥𝑝+𝑓𝑡𝑐+𝑢𝑥+𝑔𝑡+𝑐𝑢,𝑢(𝑥,𝑡)=1𝑢+𝜌𝑐𝑓𝑥𝑡𝑐+𝑢𝑥𝑔𝑡+𝑐𝑢,(2) with 𝑐 the speed of sound, 𝑇 the temperature and 𝜌 the density. The functions 𝑔 and 𝑓 are the Riemann invariants representing the up and down travelling waves. Instead of pressure and velocity, we use the Riemann invariants as acoustic states. The advantage is that the resulting transfer functions are always causal if incoming waves are taken as input and outgoing waves as output. As a consequence, every element in the one-dimensional network has to connect the Riemann invariants on both sides of the element. Straight duct elements can be represented by the time delays 𝜏𝑢 and 𝜏𝑑 which correspond to the travel times of the acoustic waves. As gas air under atmospheric conditions is assumed for the upstream part, for the downstream part with a mean temperature of 347K due to the heat release of the flame. The mean flow speed 𝑢 is neglected because of its negligible influence on the acoustic travel times compared to the speed of sound. Therefore, 𝜏𝑢 and 𝜏𝑑 can be calculated with 𝜏=𝑥𝑖/𝑐 with 𝑥𝑖 the considered distance, and 𝑐=𝛾𝑅𝑇 with 𝑅 the specific gas constant and 𝛾 the adiabatic exponent. At both ends of the tube, the model is closed by acoustic reflection coefficients 𝑅𝑢 and 𝑅𝑑 relating incoming with outgoing wave, cf. Figure 3. Acoustic reflection coefficients less than one represent the only acoustic losses (see (1)) in the model and can be interpreted as acoustic flux across the boundaries of the system.

Figure 3: Signal-oriented model of the Rijke tube in MATLAB/SIMULINK.
2.3. Flame

Since only low frequencies are of interest the flame zone is short compared to the acoustic wavelength. Thus, the flame zone represents a discontinuity for the acoustic waves. An approach derived in [7, 9] exploiting this simplification is used to model the flame zone. The flame zone can then be modelled with the conservation equations for mass, momentum, and energy: 𝐴1𝜌1𝑢1=𝐴2𝜌2𝑢2,𝐴1𝑝1+𝐴1𝜌1𝑢21=𝐴2𝑝2+𝐴2𝜌2𝑢22,𝑐𝑝𝑇1+12𝑢21𝜌1𝑢1𝐴1=𝑐𝑝𝑇2+12𝑢22𝜌2𝑢2𝐴2+𝑄,(3) with the subscripts 1 and 2 for up- and down-stream of the flame zone (cf. Figure 2), 𝑄 the heat release, 𝑐𝑝 the specific heat capacity, and 𝐴 the area. Since 𝐴1=𝐴2, the area is omitted in the following. Linearisation leads to the following system of equations:𝑋𝑓2𝑔2𝑓=𝑌1𝑔1+0𝑞,(4) with the matrices 𝑋 and 𝑌. To close the feedback between acoustics and combustion, a dynamic relation between acoustic disturbances and fluctuating heat release is needed. Because a large number of physical pathways exist [10], one often chooses to subsume these effects in one flame transfer function (the term “transfer function” is used for the Laplace transform of an ordinary differential equation. Note that the term “nonlinear transfer function” is not compatible with this definition), relating one acoustic variable to the integral heat release. In the flame model, we relate the acoustic velocity 𝑢1(𝑡) upstream of the flame to the rate of heat release 𝑞=𝑄(𝑡)𝑄. In [7], it is shown that𝑝𝑝𝛾𝑂𝑀𝑢𝑢,(5) so that at a low Mach number 𝑀 pressure fluctuations remain small even, when 𝑢/𝑢𝑂(1). The same conclusions apply for density and temperature fluctuations. Thus, it is reasonable to regard acoustic velocity fluctuations as the main excitation source for heat release fluctuations at low Mach numbers [11] as is the case for the present setup. As flame transfer function the 𝑛-𝜏-model of Crocco and Cheng [12] is often used, where 𝑛 is the so-called interaction index or combustion efficiency and 𝜏𝑓 a time delay. The time delay can be explained by convective transportation and mixing time. This model can be combined with a low-pass filter. Several researchers reported a low-pass behaviour of the flame transfer function at least for premixed flames [11, 13]. These considerations lead to the following transfer function: 𝑞𝑘(𝑠)=𝑓1+𝑇𝑓𝑠𝑒𝑠𝜏𝑓𝑢1(𝑠),(6) with 𝑘𝑓 the static gain and 𝑇𝑓 the time constant of the low-pass element, 𝜏𝑓 the travel respectively dead time and 𝑢1 the acoustic velocity upstream of the flame.

Due to the resulting interaction between acoustics and combustion, the thermoacoustic instability reaches a limit cycle after a period of exponential growth which is a clear indicator for nonlinear effects in the system. Since only linear acoustics are expected in the operational range of a combustor, the nonlinear effects are supposed to be in the combustion itself. The main non-linear effect of the combustion is possibly a saturation in the heat release. Besides sophisticated models relating the flame surface to the heat release via the G-equation [14, 15], there are empirical ways of predicting a limit cycle by using a parameter-dependent 𝑛-𝜏-model [11]. This non-linear behaviour, however, is not pursued further in this paper because the goal for active control is to keep the system near the steady state.

2.4. Model Validation

For simulation purpose, the presented model is implemented in signal oriented form in Matlab/SIMULINK as shown in Figure 3. In order to validate the model and the chosen parameters, the frequency response of the analytic model and from measurements are compared. The used modelling parameters are shown in Table 2. For the derivation of the complete set of parameters with a discussion of the diffusion flame in the setup, see [16].

Table 2: Model parameters for the physical model.

Figure 4 shows a comparison of the frequency response of the model and measurements from the test rig taken from closed-loop identification. The four illustrated positions represent measurements respectively simulation results at increasing length of the tube. Since a positive 180 phase shift indicates instability, it can be seen that at positions 2–4 the system is unstable, whereas at position 1 it is stable. Since all parameters are kept constant over the presented operating range, the change of dynamical behaviour depends only on the value of the time delays in the upper part of the tube that are determined by its length. In the presented diagrams, the length was varied from 1.09m to 1.365m resulting in a maximal difference of 1.5ms for the time delays. This explains the high sensitivity of the dynamical behaviour of the model concerning small changes in the parameters of the flame. The diagrams show that a parameter configuration could be found, that predicts instability or stability over the chosen operating range correct. The simulated resonant frequencies are in good agreement with the measurement results. There are bigger differences, however, in the amplitudes in these points especially at the unstable positions 2–4. Obviously, the amplification of acoustic fluctuations at the flame is too high at these frequencies, however there is only little room for improvement via parameter variation in the given setup of the model, as a reduction of the influence of the flame stabilises the system, assuming constant energy losses.

Figure 4: Comparison of frequency response of measurements (green) and analytical model (blue).

3. Robust Output Model Predictive Control

Model predictive control (MPC) is a control technique that utilizes modern optimization algorithms for control by using a predictive optimization formulation and algorithms for constrained optimization online. In a receding horizon policy the cost function:𝐽(𝑥(0),𝐮)=𝑁1𝑘=0𝑙(𝑥(𝑘),𝑢(𝑘))+𝑉𝑓(𝑥(𝑁))(7) is minimized over the prediction horizon 𝑁 iteratively in every time step with updated state, 𝑥(0) under the input, state and terminal constraints 𝐮𝕌, 𝑥(𝑘)𝕏 for 𝑘={0,,𝑁1} and 𝑥(𝑁)𝕏𝑓. In addition, the dynamic model of the system:𝑥(𝑘+1)=𝑓(𝑥(𝑘),𝑢(𝑘))(8) with the state 𝑥 and the input 𝑢 acts as a state constraint for 𝑘={0,,𝑁1}. The stage cost 𝑙(𝑥(𝑘),𝑢(𝑘)) is used to define the control objective, for example, by penalizing the deviation of the models future trajectory of some reference trajectory. In the following𝑙(𝑥(𝑘),𝑢(𝑘))=𝑥𝑇(𝑘)𝑄𝑥(𝑘)+𝑢𝑇(𝑘)𝑅𝑢(𝑘)(9) is assumed with the positive definite matrices 𝑄and𝑅. The terminal penalty 𝑉𝑓(𝑥(𝑁)) can be used to guarantee stability (see Section 3.3). The first entry of the minimizer 𝐮 of 𝐽() is used as the actual control input to the system. The advantage of the MPC approach over classical control approaches like PIDs is the capability of explicitly accounting for constraints in the actuating and/or state variables.

3.1. Polytopic System of the Rijke Tube

In order to explicitly take into account model uncertainties within the MPC framework, we consider a linear polytopic model of the Rijke tube in the following form:𝑥(𝑘+1)=𝐴𝑘𝑥(𝑘)+𝐵𝑘𝑢(𝑘)(10) with (𝐴𝑘,𝐵𝑘)𝐶𝑜{(𝐴1,𝐵1),,(𝐴𝑚,𝐵𝑚)}. In order to derive such a model, we follow a pragmatic approach for reduced modelling. We model the unstable first (fundamental) mode as an oscillatory second-order element and neglect higher-order dynamics:̇𝑥=01𝜔202𝐷𝜔00𝑥+𝑘𝜔20𝑢,(11) with the characteristic angular frequency 𝜔0 and the damping coefficient 𝐷. Of course this is only an approximation of the real system, but the unstable mode is by far the most dominant one in the system. Nevertheless, higher-order dynamics could destabilize the system if the controller excites them. Therefore, it is mandatory that the higher orders are filtered out. This can be accomplished using the reduced model within the state observer (see Section 3.5) with suitable tuning. Using the analytic model we found a parameter range for the damping of 𝐷={0.002,,0.012} and for the characteristic angular frequency of 𝜔0={815/s,,930/s} due to a change in the length of the Rijke tube between 𝑙={1.2m,,1.365m}. This range is chosen because it represents the unstable regime regarding the tube length in the setup. The static gain is constant with 𝑘=0.0106. As linear parameters 𝑝1=𝜔20 and 𝑝2=2𝐷𝜔0 are chosen. The resulting discretized linear parameter varying system (LPV) is 𝑥𝐴(𝑘+1)=0+𝑝1𝐴1+𝑝2𝐴2𝑥+𝐵(𝑘)0+𝑝1𝐵1+𝑝2𝐵2𝑦𝐶𝑢(𝑘),(𝑘)=0+𝑝1𝐶1+𝑝2𝐶2𝑥(𝑘).(12)

This LPV system can be transformed into a polytopic system (10) with 𝑚=4 vertices. As sampling time 𝑇𝑠=0.0006s is chosen in order to be able to resolve the frequency range of the unstable modes. This is also the time slot in which each optimization problem for the MPC has to be solved. For the proposed robust observer in Section 3.5 the actual model and therefore the parameters 𝑝1 and 𝑝2, respectively 𝐷 and 𝜔0, have to be known. Since we assume that the length of the tube can be measured online, the relation between the total tube length 𝑙 and these parameters is sought. For 𝜔0 it is a reasonable assumption to consider only the acoustics for determining the resonant frequency since the flame has negligible influence on it in the current setup. Thus, the following relation for the fundamental mode can be used:𝜔0=𝜋𝑙𝑢/𝑐𝑢+𝑙𝑑/𝑐𝑑(13) with the index 𝑢 for the lower part of the tube below the flame holder and 𝑑 the downstream part of the tube. To relate the damping 𝐷 to the tube length 𝑙, a data fit from the simulation is used. The LPV model (12) and its equivalent polytopic form (10) are used in the following for the robust state observer and for the determination of the constraints within the optimization problem.

3.2. Feedback MPC

Because an open-loop prediction of the systems state trajectory, as it is standard in most MPC formulations, cannot account for the reduced sensitivity to disturbances or modelling errors due to a feedback controller, we use the approach of incorporating an internal state feedback controller 𝑢=𝐾𝑓𝑥 in the prediction model. This technique is called closed-loop paradigm (CLP, [17]). In addition, the CLP has the advantage of a better numerical conditioning of the prediction matrices especially for an unstable system because the system with controller is stable. As a result the influence of uncertainties, in the present case represented by the polytopic system, is reduced over the prediction horizon due to the presence of the internal controller. Otherwise, the optimization problem would be far TOO conservative, because no controller action would be assumed over the prediction horizon to react on the uncertainties. The re-parametrisation starts with𝐽(𝑥(0),𝐮)𝑐𝑙𝑝=𝑘=0𝑥𝑇(𝑘)𝑄𝑥(𝑘)+𝑢𝑇(=𝑘)𝑅𝑢(𝑘)𝑁1𝑘=0𝑥𝑇(𝑘)𝑄𝑥(𝑘)+𝑢𝑇(𝑘)𝑅𝑢(𝑘)+𝑉𝑓(𝑥(𝑁))(14) with the inputs𝑢(𝑘)=𝐾𝑓𝑥(𝑘)+𝑐(𝑘),𝑘={0,,𝑁1},(15)𝑢(𝑘)=𝐾𝑓𝑥(𝑘),𝑘𝑁(16) and the terminal penalty𝑉𝑓(𝑥)=𝑘=𝑁𝑥𝑇(𝑘)𝑄𝑥(𝑘)+𝑢𝑇(𝑘)𝑅𝑢(𝑘).(17) Formulating the cost function 𝐽() in 𝐜 results in (see [17])𝐽(𝑥(0),𝐜)CLP=𝑁1𝑘=0𝑐𝑇(𝑘)𝑊𝑐(𝑘)(18) under the (virtual) input, state and terminal constraints. Since (18) is just a re-parametrisation of the original optimization problem, the optimum is identical. The optimization variable 𝐜 is the deviation of the nominal controller in (15) and is used for constrained handling only as is evident from (18). The minimizer 𝐜 is identical zero in case of no active constraints. A pragmatic view of designing a MPC in CLP formulation is to design a controller matrix 𝐾𝑓 and use 𝑊 in (18) for the weighting of 𝑐(𝑘) over the prediction horizon for constraint handling. This approach is used in the following.

3.3. Robust Stability of MPC

Formulations for guaranteed stability within the MPC framework can be considered as a mature topic nowadays. A standard approach is to utilize the cost function 𝐽(𝑥,𝐮) as a Lyapunov function by requiring that a minimizer 𝑢 formin𝑢𝕌𝑉𝑓(𝑓(𝑥,𝑢))+𝑙(𝑥,𝑢)𝑓(𝑥,𝑢)𝕏𝑓𝑉𝑓(𝑥).(19) exists for all 𝑥𝕏𝑓 (for details see [18, 19]). Note that (19) implicitly renders 𝕏𝑓 positive invariant Positive invariant means 𝑓(𝑥,𝑢)𝕏𝑓, for all 𝑥𝕏𝑓. In order to extend this requirement to the robust case, one has to consider the worst-case for guaranteed stability. This can be done by use of a min-max optimization: min𝑢𝕌maxΔ𝑁1𝑘=0𝑙(𝑥(𝑘),𝑢(𝑘),𝛿(𝑘))+𝑉𝑓(𝑥(𝑁))(20) with 𝛿(𝑘)Δ the model uncertainty. For the terminal penalty this results in the following requirement for robust stability [19]:min𝑢𝕌maxΔ𝑉𝑓(𝑓(𝑥,𝑢,𝛿))+𝑙(𝑥,𝑢)𝑓(𝑥,𝑢,𝛿)𝕏𝑓𝑉𝑓(𝑥),(21) for all 𝑥𝕏𝑓. Since Δ0 as 𝑥0 asymptotic stability of the origin can be established for the LPV model, and a solution to (21) exists in case of adequate constraints. Thus, the problem is to find a terminal penalty 𝑉𝑓() which fulfills this condition in the closed set 𝕏𝑓. In the unconstrained case and for the system (10) this problem can be cast as the following linear matrix inequality (LMI). Having (21) in mind one searches for a linear state feedback 𝑢=𝐾𝑓𝑥 and 𝑉𝑓(𝑥)=𝑥𝑇𝑃𝑓𝑥 with 𝑃𝑓>0, an upper bound of the worst-case costs, that fulfills:𝑥(𝑘)𝑇𝑄𝑥(𝑘)+𝑢(𝑘)𝑇𝑅𝑢(𝑘)+𝑥(𝑘+1)𝑇𝑃𝑓𝑥(𝑘+1)𝑥(𝑘)𝑇𝑃𝑓𝑥(𝑘)(22)for all (𝐴𝑘,𝐵𝑘)𝐶𝑜{(𝐴1,𝐵1),,(𝐴𝑚,𝐵𝑚)} and with 𝑄,𝑅>0 which makes 𝑉𝑓() a control-Lyapunov function for the complete model set. Inserting the linear state feedback 𝑢=𝐾𝑓𝑥 and the model leads to the following matrix inequality:𝐴𝑘𝐵𝑘𝐾𝑓𝑇𝑃𝑓𝐴𝑘𝐵𝑘𝐾𝑓𝑃𝑓𝑄𝐾𝑓𝑅𝐾𝑓.(23)A linear state feedback fulfilling this condition is robustly stabilizing the unconstrained LPV system. According to [20], this matrix inequality can be transformed into an LMI by substituting 𝑊=𝑃𝑓1 and 𝐿=𝐾𝑓𝑃𝑓1 and applying a Schur complement:𝑊𝐴𝑘𝑊+𝐵𝑘𝐿𝑇𝑊𝐿𝑇𝐴𝑘𝑊+𝐵𝑘𝐿𝑊00𝑊0𝑄10𝐿00𝑅10.(24)Since the model set (𝐴𝑘,𝐵𝑘) and 𝑉𝑓() are convex, it is sufficient to check the vertices of the set (𝐴𝑘,𝐵𝑘). If one searches in addition for a controller that minimizes the worst-case cost:𝑃𝑓maxΔ𝑘=0𝑥(𝑘)𝑇𝑄𝑥(𝑘)+𝑢(𝑘)𝑇𝑅𝑢(𝑘),(25)one has to solve the semi definite program:max𝑊,𝐿tr(𝑊),(26)subject to the LMI (24). When comparing (21) and (22) it can easily be seen that by choosing the stage cost in (7) to be equal to 𝑙(𝑥,𝑢)=𝑥(𝑘)𝑇𝑄𝑥(𝑘)+𝑢(𝑘)𝑇𝑅𝑢(𝑘), condition (21) is fulfilled, since𝑥(𝑘)𝑇𝑄𝑥(𝑘)+𝑢(𝑘)𝑇𝑅𝑢(𝑘)+𝑥(𝑘+1)𝑇𝑃𝑓𝑥(𝑘+1)𝑥(𝑘)𝑇𝑄𝑥(𝑘)+𝐾𝑓𝑥(𝑘)𝑇𝑅𝐾𝑓𝑥(𝑘)+𝑥(𝑘+1)𝑇𝑃𝑓𝑥(𝑘+1)𝑥(𝑘)𝑇𝑃𝑓𝑥(𝑘)(27) holds for the complete model set in the unconstrained case (if the global minimum is guaranteed to be found by the optimizer, then this condition holds. In this work a quadratic program is solved, so this condition holds). Therefore, the terminal state 𝑥(𝑁) has to reach the aforementioned terminal region 𝕏𝑓, here a region where 𝑢=𝐾𝑓𝑥𝕌 and 𝑥(𝑘+1)=(𝐴𝑘𝐵𝑘𝐾𝑓)𝑥(𝑘)𝕏 for all 𝑥𝕏𝑓 holds, that is, no active constraints. Furthermore, 𝕏𝑓 has to be a robust positive invariant set.

In the literature, this technique of guaranteed stability is sometimes called dual-mode control. The optimizer has to steer the system into the terminal region 𝕏𝑓 taking into account the constraints. Within the terminal region a virtual second controller 𝑢=𝐾𝑓𝑥 becomes active with guaranteed stability. Because the terminal region 𝕏𝑓 is designed in such a way that this controller never violates the constraints in this set, linear theory can be used to design 𝑢=𝐾𝑓𝑥, in this work by the use of LMIs. As already mentioned, (21) guarantees robust stability only in case a min-max optimization is used in standard MPC formulation. A min-max optimization is time consuming and can usually not be solved in real time, that is, within the sampling interval. In case of a re-parametrisation as feedback formulation using the derived 𝐾𝑓 though robust stability can be guaranteed for a minimization of (18) as quadratic program (QP) [6]. It can be shown that 𝐽() is Lyapunov, that is, strictly monotonically decreasing over time.

3.4. Robust Constraint Satisfaction

As terminal region 𝕏𝑓 the maximal robust positive invariant set (MAS) for polytopic systems can be constructed using the algorithm presented in [21] if it holds that𝐴𝑘𝐵𝑘𝐾𝑓𝑇𝑃𝑓𝐴𝑘𝐵𝑘𝐾𝑓<𝑃𝑓,(28) known as quadratic stability, for all 𝑘={1,,𝑚}. This is fulfilled because of (23). In order to reach the terminal set 𝕏𝑓 robustly, we follow the lines of [6] and use all permutations of model set vertices for a 𝑛-step ahead prediction:𝑥(𝑘+𝑛)=𝑛1𝑘=0𝐴𝑘𝐵𝑘𝐾𝑓+𝑥(𝑘)𝑛1𝑙=0𝑛𝑙2𝑚=0𝐴𝑛𝑙𝐵𝑛𝑙𝐾𝑓𝐵𝑙𝑐(𝑙),(29) to predict the future states. Applying the constraints to all possible future states leads to the constraint set:𝑁𝑥(𝑘)+𝑀𝐜𝑑,(30) which can be checked for redundancy using the (reformulated) constraint set𝐜(𝑁𝑀)𝑥(𝑘)𝑑.(31) Restricting the future trajectory to (30) guarantees robust constraint satisfaction for the complete model set and arbitrary fast model changes within the prediction horizon.

3.5. Robust Observer

Since the actual state 𝑥(0) is not measurable in the Rijke tube setup it has to be estimated by a state observer. In this work, the observer structure of a Luenberger observer is used ̂𝑥(𝑘+1)=𝐴𝑘̂𝑥(𝑘)+𝐿𝑜(𝑦(𝑘)̂𝑦(𝑘))+𝐵𝑘𝑢(𝑘),̂𝑦(𝑘)=𝐶𝑘̂𝑥(𝑘)+𝐷𝑘𝑢(𝑘),(32) where ̂𝑥(𝑘) and ̂𝑦(𝑘) are the estimated state and output, and 𝐿𝑜 is the observer gain. Defining the observer error as 𝑒(𝑘)=̂𝑥(𝑘)𝑥(𝑘) leads to𝑒𝐴(𝑘+1)=𝑘𝐿𝑜𝐶𝑘𝑒(𝑘).(33) Thus, the observer error converges to zero if (33) is asymptotically stable for all (𝐴𝑘,𝐵𝑘). To do so, one can make use of the duality between the control and observer problem. Since the eigenvalues, and therefore the stability properties of 𝐴𝐿𝑜𝐶 are identical with its transposed (𝐴𝑘𝐿𝑜𝐶𝑘)𝑇=𝐴𝑇𝑘𝐶𝑇𝑘𝐿𝑇𝑜 the observer system (32) can be cast into a control problem [22]: 𝑥𝑇(𝑘+1)=𝐴𝑇𝑘𝑥𝑇(𝑘)+𝐶𝑇𝑘𝑢𝑇𝑢(𝑘),𝑇(𝑘)=𝐿𝑇𝑜𝑥𝑇(𝑘)(34) and one can use the semi-definite program shown in Section 3.3 to find a positive definite 𝑃𝑜 and the corresponding observer gain 𝐿𝑜 fulfilling the (quadratic stability) condition 𝐴𝑇𝑘𝐶𝑇𝑘𝐿𝑇𝑜𝑇𝑃𝑜𝐴𝑇𝑘𝐶𝑇𝑘𝐿𝑇𝑜𝑃𝑜𝑄𝐿𝑜𝑅𝐿𝑇𝑜,(35)for all (𝐴𝑇𝑘,𝐶𝑇𝑘)𝐶𝑜{(𝐴𝑇1,𝐶𝑇1),,(𝐴𝑇𝑚,𝐶𝑇𝑚)} with 𝑄,𝑅>0. As a consequence, (33) is robustly asymptotically stable for the model set. Thus, the observer error converges to zero for the complete model set, if the actual model is known to the observer. Here, the model is known due to the online measurement of the tube length.

4. Results

The proposed robust MPC with its internal polytopic model and the robust state observer is able to stabilize the Rijke tube robustly by solving a quadratic program (QP) online. Critical for the application are the number of constraints in the QP as a result of the feedback MPC formulation. For a prediction horizon of 𝑁=4 and 𝑚=4 vertices in the polytopic model under the input constraints 5𝑢5 the problem formulation leads to a QP with 180 constraints. This QP can be solved online within the sampling time of 𝑇𝑠=0.0006s with modern hardware. In Figure 5 the control of the physical model is shown for four different lengths of the Rijke tube in simulation. The controller is activated after 2.1s. This delay can be interpreted as a temporary loss of actuation or a disturbance deflecting the system from steady state. The controller has the task to drive the system back to the origin. It can be seen that the RMPC is capable of robustly stabilizing the different configurations while respecting the input constraints. Especially for the length 𝑙=1.285m, the RMPC uses the full input range. The reason for this extensive use of the actuator can be seen in the fact that the RMPC has to guarantee stability for any arbitrarily fast model change (change in tube length) within the prediction horizon. This conservatism could be reduced by using only the actual model for constraint satisfaction while sacrificing some robustness concerning model changes. Then, however, the guaranteed stability would be lost during a parameter variation, that is, during transition from one tube length to another. Another source of conservatism is the formulation as feedback MPC. The robust optimization problem, that is, finding an input trajectory which guarantees constraint satisfaction in the presence of (known) uncertainties, can be solved exactly for example by dynamic programming [19]. Unfortunately, dynamic programming is not applicable for most systems due to its complexity. The feedback MPC formulation is a good compromise between conservatism and optimization complexity. Figure 6 presents the control results for the RMPC of the Rijke tube for 𝑙=1.2m and 𝑙=1.365m in the experiment. The successful control proofs the real-time capability of the control. As input and output constraints 7𝑢7 and 1𝑦1 are used leading to 176 constraints for the QP. The maximum calculation time on the dSPACE rapid control prototyping hardware DS1006 equipped with a 2.8GHz AMD Opteron CPU is below the sampling time of 𝑇𝑠=0.0006s. As QP solver qpOASES was used [23]. The robust state observer uses the presented LPV model which was derived from the analytical model. Using the estimated state, the RMPC can stabilize the Rijke tube while respecting the given constraints. For comparison the output estimated by the robust state observer is plotted. It can be seen that the frequency is in excellent agreement with the measurement from the microphone, and the estimated amplitude is a little bit to conservative, when the controller becomes active.

Figure 5: RMPC of the unstable Rijke tube with different length in simulation. The RMPC has a prediction horizon of 𝑁=4 and uses a robust state observer. The RMPC is activated at 𝑡=2.1s.
Figure 6: RMPC of the unstable Rijke tube with different length in the experiment. The RMPC has a prediction horizon of 𝑁=4 and uses a robust state observer. The maximum calculation time for 𝑙=1.2m is 3.35×105s and for 𝑙=1.365m1.9×104s.

5. Conclusions

An analytical model for a Rijke tube has been applied which is able to reproduce the stability map of the thermoacoustic setup as well as the dynamic behaviour over different lengths of the tube. Therefore, it is an ideal test bench for the robust control of unstable thermoacoustic systems, especially the real-time capability of MPC algorithms. It is used to derive a simplistic linear parameter varying system which represents the unstable modes. Using this system, it is shown that a robust output MPC can be designed that is capable of steering the system robustly to the origin under constraints. The proposed RMPC is a good compromise between conservatism and computational load when considering the very fast system dynamics of a thermoacoustic system and as a consequence the short calculation times needed for the application online. Thus, RMPC is a promising approach in order to fulfill the demands for an active control system in modern gas turbines. The robust stabilization and constraint handling are shown in this paper. In addition, a MIMO setup and fault-tolerant control can be incorporated quiet naturally into the MPC framework. Finally, the importance of physical understanding and modelling for the estimation of unavoidable uncertainties is demonstrated for the robust control of thermoacoustic instabilities.


The authors gratefully acknowledge the contribution of the Deutsche Forschungsgemeinschaft through the Collaborative Research Center 686 “Model-Based Control of Homogenized Low-Temperature Combustion”.


  1. J. S. W. Rayleigh, The Theory of Sound, vol. 2, Courier Dover Publications, 1945.
  2. F. Nicoud and T. Poinsot, “Thermoacoustic instabilities: should the Rayleigh criterion be extended to include entropy changes?” Combustion and Flame, vol. 142, no. 1-2, pp. 153–159, 2005. View at Publisher · View at Google Scholar · View at Scopus
  3. S. Candel, “Combustion dynamics and control: progress and challenges,” in Proceedings of the 29th International Symposium on Combustion Hokkaido University Sapporo Japan, vol. 29, pp. 1–28, July 2002. View at Scopus
  4. A. P. Dowling and A. S. Morgans, “Feedback control of combustion oscillations,” Annual Review of Fluid Mechanics, vol. 37, pp. 151–182, 2005. View at Publisher · View at Google Scholar · View at Scopus
  5. B. T. Zinn and T. C. Lieuwen, “Combustion instabilities: basic concepts,” in Combustion Instabilities in Gas Turbine Engines: Operational Experience, Fundamental Mechanisms, and Modeling, T. C. Lieuwen and V. Yang, Eds., vol. 210, pp. 3–26, American Institute of Aeronautics and Astronautics, 2005.
  6. B. Pluymers, J. A. Rossiter, J. Suykens, and B. De Moor, “A simple algorithm for robust MPC,” in Proceedings of the IFACWorld Congress, Prague, Czech Republic, 2005.
  7. A. P. Dowling, “Nonlinear self-excited oscillations of a ducted flame,” Journal of Fluid Mechanics, vol. 346, pp. 271–290, 1997. View at Scopus
  8. T. Poinsot and D. Veynante, Theoretical and Numerical Combustion, R. T. Edwards, 2nd edition, 2005.
  9. A. S. Morgans and A. P. Dowling, “Model-based control of combustion instabilities,” Journal of Sound and Vibration, vol. 299, no. 1-2, pp. 261–282, 2007. View at Publisher · View at Google Scholar · View at Scopus
  10. T. Lieuwen, “Modeling premixed combustion-acoustic wave interactions: a review,” Journal of Propulsion and Power, vol. 19, no. 5, pp. 765–781, 2003. View at Scopus
  11. P. A. Hield, M. J. Brear, and S. H. Jin, “Thermoacoustic limit cycles in a premixed laboratory combustor with open and choked exits,” Combustion and Flame, vol. 156, no. 9, pp. 1683–1697, 2009. View at Publisher · View at Google Scholar · View at Scopus
  12. L. Crocco and S. L. Cheng, “Theory of combustion instability in liquid propellant rocket motors,” AGARDOGRAPH 8, Butterworths Science Publication, 1956.
  13. T. Schuller, D. Durox, and S. Candel, “A unified model for the prediction of laminar flame transfer functions: comparisons between conical and V-flame dynamics,” Combustion and Flame, vol. 134, no. 1-2, pp. 21–34, 2003. View at Publisher · View at Google Scholar · View at Scopus
  14. A. P. Dowling, “A kinematic model of a ducted flame,” Journal of Fluid Mechanics, vol. 394, pp. 51–72, 1999. View at Scopus
  15. X. Yuan, K. Glover, and A. P. Dowling, “Modeling investigation for thermoacoustic oscillation control,” in Proceedings of the American Control Conference, pp. 3323–3328, Baltimore, Md, USA, 2010.
  16. F. Jarmolowitz, Ch. Gro-Weege, Th. Lammersen, and D. Abel, “Modelling and robust Model Predictive Control of an unstable thermoacoustic system with input constraints,” in Proceedings of the American Control Conference, Montreal, Canada, 2012.
  17. J. A. Rossiter, Model-Based Predictive Control: A Practical Approach, CRC Press, 2003.
  18. D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert, “Constrained model predictive control: stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000. View at Publisher · View at Google Scholar · View at Scopus
  19. J. B. Rawlings and D. Q. Mayne, Model Predictive Control: Theory and Design, Nob Hill, 2009.
  20. S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear Matrix Inequalities in System and Control Theory, SIAM Studies in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 1994.
  21. B. Pluymers, J. A. Rossiter, J. Suykens, and B. De Moor, “The effcient computation of polyhedral invariant sets for linear systems with polytopic uncertainty,” in Proceedings of the American Control Conference, pp. 804–809, Portland, Ore, USA, 2005.
  22. J. Lunze, Regelungstechnik 2: Mehrgroensysteme, Digitale Regelung, Springer, 2010.
  23. H. J. Ferreau, H. G. Bock, and M. Diehl, “An online active set strategy to overcome the limitations of explicit MPC,” International Journal of Robust and Nonlinear Control, vol. 18, no. 8, pp. 816–830, 2008. View at Publisher · View at Google Scholar · View at Scopus