Robust Output Model Predictive Control of an Unstable Rijke Tube

Abstract

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: 𝐴₁𝐴₁𝜌₁𝑢₁=𝐴₂𝜌₂𝑢₂,𝑐_𝑝𝑝₁+𝐴₁𝜌₁𝑢²¹=𝐴₂𝑝₂+𝐴₂𝜌₂𝑢²²,𝑇₁+12𝑢²¹𝜌¹𝑢₁𝐴₁=𝑐_𝑝𝑇₂+12𝑢²²𝜌²𝑢₂𝐴₂+𝑄,(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 𝐴₁=𝐴₂, the area is omitted in the following. Linearisation leads to the following system of equations:𝑋⎛⎜⎜⎝𝑓²𝑔₂⎞⎟⎟⎠=𝑌⎛⎜⎜⎝𝑓₁𝑔₁⎞⎟⎟⎠+⎛⎜⎜⎝𝑞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 𝑢₁ 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 (𝐴_𝑘,𝐵_𝑘)∈𝐶𝑜{(𝐴₁,𝐵₁),…,(𝐴_𝑚,𝐵_𝑚)}. 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:𝑥̇=⎛⎜⎜⎝−0𝜔₂₀1−2𝐷𝜔₀⎞⎟⎟⎠𝑥+⎛⎜⎜⎝𝑘𝜔0₂₀⎞⎟⎟⎠𝑢,(11) with the characteristic angular frequency 𝜔₀ 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 𝜔₀={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 𝑝₁=−𝜔²⁰ and 𝑝₂=−2𝐷𝜔₀ are chosen. The resulting discretized linear parameter varying system (LPV) is 𝑥(𝑘+1)=+𝐴𝐵₀+𝑝₁𝐴₁+𝑝₂𝐴₂𝑥(𝑘)₀𝑦(𝑘)=𝐶₀+𝑝₁𝐵₁+𝑝₂𝐵₂𝑢(𝑘),+𝑝₁𝐶₁+𝑝₂𝐶₂𝑥(𝑘).(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 𝑝₁ and 𝑝₂, respectively 𝐷 and 𝜔₀, 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 𝜔₀ 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:𝜔₀=𝜋𝑙_𝑢/𝑐_𝑢+𝑙_𝑑/𝑐_𝑑(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_𝑢∈𝕌≤ma𝑉x_Δ𝑓𝑉_𝑓(𝑓(𝑥,𝑢,𝛿))+𝑙(𝑥,𝑢)∣𝑓(𝑥,𝑢,𝛿)∈𝕏_𝑓(𝑥),(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 (𝐴_𝑘,𝐵_𝑘)∈𝐶𝑜{(𝐴₁,𝐵₁),…,(𝐴_𝑚,𝐵_𝑚)} 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:⎛⎜⎜⎜⎜⎜⎜⎝𝑊𝐴_𝑘𝑊+𝐵_𝑘𝐿^𝑇𝑊𝐿^𝑇𝐴_𝑘𝑊𝑊+𝐵_𝑘𝐿𝑊0𝑄0₋₁0𝐿00𝑅0₋₁⎞⎟⎟⎟⎟⎟⎟⎠≥0.(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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

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.

(a)

(b)

(c)

(d)

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×10−5s and for 𝑙=1.365m1.9×10−4s.

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.

Acknowledgments

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