# Robust Output Model Predictive Control of an Unstable Rijke Tube

Authors: Fabian Jarmolowitz; Christopher Groß-Weege; Thomas Lammersen; Dirk Abel
Journal: Journal of Combustion
 (2012)

---

## 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.


---

## Body

## 1. Introduction

Modern gas turbines have to comply with increasingly stringent emission requirements forNOx. 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 NOx 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  NOx 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 ∫0Tp′(t)·q̇′(t)dt>0 [1], but for example in [2] it is shown, that the acoustic losses L 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: (1)∫0Tp′(t)⋅q̇′(t)dt>L,
where p′(t) is the acoustic pressure, q̇′ the fluctuating rate of heat release, and L the acoustic losses. If the integral over one period of oscillation T 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 L, 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. Figure1. 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. Section2 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 Figure2 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.
SizeValuelu0.275mlm0.195mls0.1md0.044mlPos. 1:1.090mPos. 2:1.180mPos. 3:1.285mPos. 4:1.365mFigure 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 n=2. These assumptions are fulfilled for a Rijke tube [9]. The state can be defined, for example, by acoustic velocity u′ and acoustic pressure p′. The total pressure and velocity are defined as the sum of mean and acoustic value as  p(x,t)=p¯+p′(x,t), and  u(x,t)=u¯+u′(x,t). 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 (2)p(x,t)=p¯+f(t-xc¯+u¯)+g(t+xc¯-u¯),u(x,t)=u¯+1ρ¯c¯[f(t-xc¯+u¯)-g(t+xc¯-u¯)],
with c the speed of sound,  T the temperature and ρ the density. The functions g and f 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 τu and τd 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 u¯ is neglected because of its negligible influence on the acoustic travel times compared to the speed of sound. Therefore, τu and τd can be calculated with τ=xi/c with xi  the considered distance, and c=γRT with R the specific gas constant and γ the adiabatic exponent. At both ends of the tube, the model is closed by acoustic reflection coefficients Ru and Rd 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:
(3)A1ρ1u1=A2ρ2u2,A1p1+A1ρ1u12=A2p2+A2ρ2u22,(cpT1+12u12)ρ1u1A1=(cpT2+12u22)ρ2u2A2+Q,
with the subscripts 1 and 2 for up- and down-stream of the flame zone (cf. Figure 2), Q the heat release, cp the specific heat capacity, and A the area. Since A1=A2, the area is omitted in the following. Linearisation leads to the following system of equations:(4)X(f2g2)=Y(f1g1)+(0q′),
with the matrices X and Y. 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 u1′(t) upstream of the flame to the rate of heat release q′=Q(t)-Q¯. In [7], it is shown that(5)p′p¯∝O(γM¯)u′u¯,
so that at a low Mach number M¯ pressure fluctuations remain small even, when u′/u¯∝O(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 n-τ-model of Crocco and Cheng [12] is often used, where n is the so-called interaction index or combustion efficiency and τf 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:
(6)q′(s)=kf1+Tfse-sτfu1′(s),
with kf the static gain and Tf the time constant of the low-pass element, τf the travel respectively dead time and u1 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 theG-equation [14, 15], there are empirical ways of predicting a limit cycle by using a parameter-dependent n-τ-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 Figure3. 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.
ParameterValuec1343m/sc2374m/sRu−0.89Rd−0.89T1=Tu293KT2=Td347KTf0.00391sτf0.0018skf918Figure4 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).


## 2.1. Setup

A schematic drawing of the Rijke tube setup is shown in Figure2 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.
SizeValuelu0.275mlm0.195mls0.1md0.044mlPos. 1:1.090mPos. 2:1.180mPos. 3:1.285mPos. 4:1.365mFigure 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 n=2. These assumptions are fulfilled for a Rijke tube [9]. The state can be defined, for example, by acoustic velocity u′ and acoustic pressure p′. The total pressure and velocity are defined as the sum of mean and acoustic value as  p(x,t)=p¯+p′(x,t), and  u(x,t)=u¯+u′(x,t). 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 (2)p(x,t)=p¯+f(t-xc¯+u¯)+g(t+xc¯-u¯),u(x,t)=u¯+1ρ¯c¯[f(t-xc¯+u¯)-g(t+xc¯-u¯)],
with c the speed of sound,  T the temperature and ρ the density. The functions g and f 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 τu and τd 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 u¯ is neglected because of its negligible influence on the acoustic travel times compared to the speed of sound. Therefore, τu and τd can be calculated with τ=xi/c with xi  the considered distance, and c=γRT with R the specific gas constant and γ the adiabatic exponent. At both ends of the tube, the model is closed by acoustic reflection coefficients Ru and Rd 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:
(3)A1ρ1u1=A2ρ2u2,A1p1+A1ρ1u12=A2p2+A2ρ2u22,(cpT1+12u12)ρ1u1A1=(cpT2+12u22)ρ2u2A2+Q,
with the subscripts 1 and 2 for up- and down-stream of the flame zone (cf. Figure 2), Q the heat release, cp the specific heat capacity, and A the area. Since A1=A2, the area is omitted in the following. Linearisation leads to the following system of equations:(4)X(f2g2)=Y(f1g1)+(0q′),
with the matrices X and Y. 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 u1′(t) upstream of the flame to the rate of heat release q′=Q(t)-Q¯. In [7], it is shown that(5)p′p¯∝O(γM¯)u′u¯,
so that at a low Mach number M¯ pressure fluctuations remain small even, when u′/u¯∝O(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 n-τ-model of Crocco and Cheng [12] is often used, where n is the so-called interaction index or combustion efficiency and τf 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:
(6)q′(s)=kf1+Tfse-sτfu1′(s),
with kf the static gain and Tf the time constant of the low-pass element, τf the travel respectively dead time and u1 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 theG-equation [14, 15], there are empirical ways of predicting a limit cycle by using a parameter-dependent n-τ-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 Figure3. 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.
ParameterValuec1343m/sc2374m/sRu−0.89Rd−0.89T1=Tu293KT2=Td347KTf0.00391sτf0.0018skf918Figure4 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:(7)J(x(0),u):=∑k=0N-1l(x(k),u(k))+Vf(x(N))
is minimized over the prediction horizon N iteratively in every time step with updated state, x(0) under the input, state and terminal constraints u∈𝕌, x(k)∈𝕏 for k={0,…,N-1} and x(N)∈𝕏f. In addition, the dynamic model of the system:(8)x(k+1)=f(x(k),u(k))
with the state x and the input u acts as a state constraint for k={0,…,N-1}. The stage cost l(x(k),u(k)) 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)l(x(k),u(k)):=xT(k)Qx(k)+uT(k)Ru(k)
is assumed with the positive definite matrices QandR. The terminal penalty Vf(x(N)) can be used to guarantee stability (see Section 3.3). The first entry of the minimizer u* of  J(·) 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:(10)x(k+1)=Akx(k)+Bku(k)
with (Ak,Bk)∈Co{(A1,B1),…,(Am,Bm)}. 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:(11)ẋ=(01-ω02-2Dω0)x+(0kω02)u,
with the characteristic angular frequency ω0 and the damping coefficient  D. 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  D={-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  l={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 k=0.0106. As linear parameters p1=-ω02  and  p2=-2Dω0 are chosen. The resulting discretized linear parameter varying system (LPV) is (12)x(k+1)=(A0+p1A1+p2A2)x(k)+(B0+p1B1+p2B2)u(k),y(k)=(C0+p1C1+p2C2)x(k).This LPV system can be transformed into a polytopic system (10) with m=4 vertices. As sampling time Ts=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  p1 and p2, respectively D 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 l 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:(13)ω0=πlu/cu+ld/cd
with the index u for the lower part of the tube below the flame holder and d the downstream part of the tube. To relate the damping D to the tube length l, 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 controlleru=-Kfx 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(14)J(x(0),u)clp:=∑k=0∞xT(k)Qx(k)+uT(k)Ru(k)=∑k=0N-1xT(k)Qx(k)+uT(k)Ru(k)+Vf(x(N))
with the inputs(15)u(k)=-Kfx(k)+c(k),k={0,…,N-1},(16)u(k)=-Kfx(k),k≥N
and the terminal penalty(17)Vf(x)=∑k=N∞xT(k)Qx(k)+uT(k)Ru(k).
Formulating the cost function J(·) in c results in (see [17])(18)J(x(0),c)CLP:=∑k=0N-1cT(k)Wc(k)
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 c is the deviation of the nominal controller in (15) and is used for constrained handling only as is evident from (18). The minimizer c* 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 Kf and use W in (18) for the weighting of c(k) 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 functionJ(x,u) as a Lyapunov function by requiring that a minimizer u* for(19)min⁡u∈U{Vf(f(x,u))+l(x,u)∣f(x,u)∈Xf}≤Vf(x).
exists for all x∈𝕏f (for details see [18, 19]). Note that (19) implicitly renders 𝕏f positive invariant Positive invariant means f(x,u)∈𝕏f, for all x∈𝕏f. 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:
(20)min⁡u∈Umax⁡Δ∑k=0N-1l(x(k),u(k),δ(k))+Vf(x(N))
with δ(k)∈Δ the model uncertainty. For the terminal penalty this results in the following requirement for robust stability [19]:(21)min⁡u∈Umax⁡Δ{Vf(f(x,u,δ))+l(x,u)∣f(x,u,δ)∈Xf}≤Vf(x),
for all  x∈𝕏f. Since  Δ↦0  as  x↦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  Vf(·) which fulfills this condition in the closed set  𝕏f. 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 u=-Kfx  and  Vf(x)=xTPfx with Pf>0, an upper bound of the worst-case costs, that fulfills:(22)x(k)TQx(k)+u(k)TRu(k)+x(k+1)TPfx(k+1)≤x(k)TPfx(k)for all (Ak,Bk)∈Co{(A1,B1),…,(Am,Bm)}  and with Q,R>0  which makes Vf(·) a control-Lyapunov function for the complete model set. Inserting the linear state feedback u=-Kfx and the model leads to the following matrix inequality:(23)(Ak-BkKf)TPf(Ak-BkKf)-Pf≤-Q-Kf′RKf.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  W=Pf-1  and  L=KfPf-1 and applying a Schur complement:(24)(W(AkW+BkL)TWLTAkW+BkLW00W0Q-10L00R-1)≥0.Since the model set (Ak,Bk) and Vf(·) are convex, it is sufficient to check the vertices of the set (Ak,Bk). If one searches in addition for a controller that minimizes the worst-case cost:(25)Pf≥max⁡Δ∑k=0∞x(k)TQx(k)+u(k)TRu(k),one has to solve the semi definite program:(26)max⁡W,Ltr⁡(W),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 l(x,u)=x(k)TQx(k)+u(k)TRu(k), condition (21) is fulfilled, since(27)x(k)TQx(k)+u*(k)TRu*(k)+x(k+1)TPfx(k+1)≤x(k)TQx(k)+(-Kfx(k))TR(-Kfx(k))+x(k+1)TPfx(k+1)≤x(k)TPfx(k)
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 x(N) has to reach the aforementioned terminal region 𝕏f, here a region where u=-Kfx∈𝕌 and x(k+1)=(Ak-BkKf)x(k)∈𝕏 for all x∈𝕏f holds, that is, no active constraints. Furthermore, 𝕏f 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𝕏f taking into account the constraints. Within the terminal region a virtual second controller u=-Kfx becomes active with guaranteed stability. Because the terminal region 𝕏f is designed in such a way that this controller never violates the constraints in this set, linear theory can be used to design u=-Kfx, 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 Kf though robust stability can be guaranteed for a minimization of (18) as quadratic program (QP) [6]. It can be shown that J(·) is Lyapunov, that is,  strictly monotonically decreasing over time.


### 3.4. Robust Constraint Satisfaction

As terminal region𝕏f the maximal robust positive invariant set (MAS) for polytopic systems can be constructed using the algorithm presented in [21] if it holds that(28)(Ak-BkKf)TPf(Ak-BkKf)<Pf,
known as quadratic stability, for all  k={1,…,m}. This is fulfilled because of (23). In order to reach the terminal set 𝕏f robustly, we follow the lines of [6] and use all permutations of model set vertices for a n-step ahead prediction:(29)x(k+n)=∏k=0n-1(Ak-BkKf)x(k)+∑l=0n-1{∏m=0n-l-2(An-l-Bn-lKf)}Blc(l),
to predict the future states. Applying the constraints to all possible future states leads to the constraint set:(30)Nx(k)+Mc≤d,
which can be checked for redundancy using the (reformulated) constraint set(31)(NM)(x(k)c)≤d.
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 statex(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 (32)x̂(k+1)=Akx̂(k)+Lo(y(k)-ŷ(k))+Bku(k),ŷ(k)=Ckx̂(k)+Dku(k),
where x̂(k) and ŷ(k) are the estimated state and output, and Lo is the observer gain. Defining the observer error as e(k)=x̂(k)-x(k) leads to(33)e(k+1)=(Ak-LoCk)e(k).
Thus, the observer error converges to zero if (33) is asymptotically stable for all (Ak,Bk). 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 A-LoC are identical with its transposed (Ak-LoCk)T=AkT-CkTLoT the observer system (32) can be cast into a control problem [22]: (34)xT(k+1)=AkTxT(k)+CkTuT(k),uT(k)=-LoTxT(k)
and one can use the semi-definite program shown in Section 3.3  to find a positive definite Po and the corresponding observer gain Lo fulfilling the (quadratic stability) condition
(35)(AkT-CkTLoT)TPo(AkT-CkTLoT)-Po≤-Q-LoRLoT,for all (AkT,CkT)∈Co{(A1T,C1T),…,(AmT,CmT)} with Q,R>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.


## 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:(10)x(k+1)=Akx(k)+Bku(k)
with (Ak,Bk)∈Co{(A1,B1),…,(Am,Bm)}. 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:(11)ẋ=(01-ω02-2Dω0)x+(0kω02)u,
with the characteristic angular frequency ω0 and the damping coefficient  D. 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  D={-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  l={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 k=0.0106. As linear parameters p1=-ω02  and  p2=-2Dω0 are chosen. The resulting discretized linear parameter varying system (LPV) is (12)x(k+1)=(A0+p1A1+p2A2)x(k)+(B0+p1B1+p2B2)u(k),y(k)=(C0+p1C1+p2C2)x(k).This LPV system can be transformed into a polytopic system (10) with m=4 vertices. As sampling time Ts=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  p1 and p2, respectively D 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 l 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:(13)ω0=πlu/cu+ld/cd
with the index u for the lower part of the tube below the flame holder and d the downstream part of the tube. To relate the damping D to the tube length l, 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 controlleru=-Kfx 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(14)J(x(0),u)clp:=∑k=0∞xT(k)Qx(k)+uT(k)Ru(k)=∑k=0N-1xT(k)Qx(k)+uT(k)Ru(k)+Vf(x(N))
with the inputs(15)u(k)=-Kfx(k)+c(k),k={0,…,N-1},(16)u(k)=-Kfx(k),k≥N
and the terminal penalty(17)Vf(x)=∑k=N∞xT(k)Qx(k)+uT(k)Ru(k).
Formulating the cost function J(·) in c results in (see [17])(18)J(x(0),c)CLP:=∑k=0N-1cT(k)Wc(k)
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 c is the deviation of the nominal controller in (15) and is used for constrained handling only as is evident from (18). The minimizer c* 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 Kf and use W in (18) for the weighting of c(k) 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 functionJ(x,u) as a Lyapunov function by requiring that a minimizer u* for(19)min⁡u∈U{Vf(f(x,u))+l(x,u)∣f(x,u)∈Xf}≤Vf(x).
exists for all x∈𝕏f (for details see [18, 19]). Note that (19) implicitly renders 𝕏f positive invariant Positive invariant means f(x,u)∈𝕏f, for all x∈𝕏f. 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:
(20)min⁡u∈Umax⁡Δ∑k=0N-1l(x(k),u(k),δ(k))+Vf(x(N))
with δ(k)∈Δ the model uncertainty. For the terminal penalty this results in the following requirement for robust stability [19]:(21)min⁡u∈Umax⁡Δ{Vf(f(x,u,δ))+l(x,u)∣f(x,u,δ)∈Xf}≤Vf(x),
for all  x∈𝕏f. Since  Δ↦0  as  x↦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  Vf(·) which fulfills this condition in the closed set  𝕏f. 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 u=-Kfx  and  Vf(x)=xTPfx with Pf>0, an upper bound of the worst-case costs, that fulfills:(22)x(k)TQx(k)+u(k)TRu(k)+x(k+1)TPfx(k+1)≤x(k)TPfx(k)for all (Ak,Bk)∈Co{(A1,B1),…,(Am,Bm)}  and with Q,R>0  which makes Vf(·) a control-Lyapunov function for the complete model set. Inserting the linear state feedback u=-Kfx and the model leads to the following matrix inequality:(23)(Ak-BkKf)TPf(Ak-BkKf)-Pf≤-Q-Kf′RKf.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  W=Pf-1  and  L=KfPf-1 and applying a Schur complement:(24)(W(AkW+BkL)TWLTAkW+BkLW00W0Q-10L00R-1)≥0.Since the model set (Ak,Bk) and Vf(·) are convex, it is sufficient to check the vertices of the set (Ak,Bk). If one searches in addition for a controller that minimizes the worst-case cost:(25)Pf≥max⁡Δ∑k=0∞x(k)TQx(k)+u(k)TRu(k),one has to solve the semi definite program:(26)max⁡W,Ltr⁡(W),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 l(x,u)=x(k)TQx(k)+u(k)TRu(k), condition (21) is fulfilled, since(27)x(k)TQx(k)+u*(k)TRu*(k)+x(k+1)TPfx(k+1)≤x(k)TQx(k)+(-Kfx(k))TR(-Kfx(k))+x(k+1)TPfx(k+1)≤x(k)TPfx(k)
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 x(N) has to reach the aforementioned terminal region 𝕏f, here a region where u=-Kfx∈𝕌 and x(k+1)=(Ak-BkKf)x(k)∈𝕏 for all x∈𝕏f holds, that is, no active constraints. Furthermore, 𝕏f 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𝕏f taking into account the constraints. Within the terminal region a virtual second controller u=-Kfx becomes active with guaranteed stability. Because the terminal region 𝕏f is designed in such a way that this controller never violates the constraints in this set, linear theory can be used to design u=-Kfx, 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 Kf though robust stability can be guaranteed for a minimization of (18) as quadratic program (QP) [6]. It can be shown that J(·) is Lyapunov, that is,  strictly monotonically decreasing over time.


## 3.4. Robust Constraint Satisfaction

As terminal region𝕏f the maximal robust positive invariant set (MAS) for polytopic systems can be constructed using the algorithm presented in [21] if it holds that(28)(Ak-BkKf)TPf(Ak-BkKf)<Pf,
known as quadratic stability, for all  k={1,…,m}. This is fulfilled because of (23). In order to reach the terminal set 𝕏f robustly, we follow the lines of [6] and use all permutations of model set vertices for a n-step ahead prediction:(29)x(k+n)=∏k=0n-1(Ak-BkKf)x(k)+∑l=0n-1{∏m=0n-l-2(An-l-Bn-lKf)}Blc(l),
to predict the future states. Applying the constraints to all possible future states leads to the constraint set:(30)Nx(k)+Mc≤d,
which can be checked for redundancy using the (reformulated) constraint set(31)(NM)(x(k)c)≤d.
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 statex(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 (32)x̂(k+1)=Akx̂(k)+Lo(y(k)-ŷ(k))+Bku(k),ŷ(k)=Ckx̂(k)+Dku(k),
where x̂(k) and ŷ(k) are the estimated state and output, and Lo is the observer gain. Defining the observer error as e(k)=x̂(k)-x(k) leads to(33)e(k+1)=(Ak-LoCk)e(k).
Thus, the observer error converges to zero if (33) is asymptotically stable for all (Ak,Bk). 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 A-LoC are identical with its transposed (Ak-LoCk)T=AkT-CkTLoT the observer system (32) can be cast into a control problem [22]: (34)xT(k+1)=AkTxT(k)+CkTuT(k),uT(k)=-LoTxT(k)
and one can use the semi-definite program shown in Section 3.3  to find a positive definite Po and the corresponding observer gain Lo fulfilling the (quadratic stability) condition
(35)(AkT-CkTLoT)TPo(AkT-CkTLoT)-Po≤-Q-LoRLoT,for all (AkT,CkT)∈Co{(A1T,C1T),…,(AmT,CmT)} with Q,R>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 ofN=4 and m=4 vertices in the polytopic model under the input constraints  -5≤u≤5  the problem formulation leads to a QP with 180 constraints. This QP can be solved online within the sampling time of Ts=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 l=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 l=1.2m and l=1.365m in the experiment. The successful control proofs the real-time capability of the control. As input and output constraints -7≤u≤7  and  -1≤y≤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 Ts=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.RMPC of the unstable Rijke tube with different length in simulation. The RMPC has a prediction horizon ofN=4 and uses a robust state observer. The RMPC is activated at t=2.1s.
(a)(b)(c)(d)(e)(f)(g)(h)RMPC of the unstable Rijke tube with different length in the experiment. The RMPC has a prediction horizon ofN=4 and uses a robust state observer. The maximum calculation time for l=1.2m is 3.35×10-5s and for l=1.365m1.9×10-4s.
(a)(b)(c)(d)


## 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.



---

*Source: 927345-2012-05-27.xml*