Model Predictive Control of A Gas Turbine

Model Predictive Control of a Gas Turbine Installation
Bas Vroemen WFW report 97.002
B.G. VROEMEN Faculty of Mechanical Engineering Eindhoven University of Technology January 1997
I would like to thank all people who helped m e fulfill this project, but Harm van Essen in particular, whose pleasant support I couldnt have managed without
Furthermore, I wish to express m y gratitude towards m y parents, who have made this all possible for m e
To the memory of Tante Liesbeth
Model Predictive Control of a Gas Turbine Installation
B. G. Vroemen
Faculty of Mechanical Engineering Eindhoven University of Technology January 15, 1997
Final project Supervisor: Harm van Essen Professor: J. J. Kok
Nomenclat ure
Summary
1 Introduction
11
1.1 Gas turbines
....................................
11
12 13 13
14
1.2 Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Goals of this project
................................
1.4 Projectsetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Organization of report . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Theory: MPC
2.1
15
Unconstrained Model Predictive Control

2.1.1
.....................
15
16
The optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2
Constrained Model Predictive Control
......................
17
17
2.3 Tuning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 2.3.2 2.3.3 The sample time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The prediction horizon . . . . . . . . . . . . . . . . . . . . . . . . . . . The control horizon
18 18
20
............................
2
2.3.4 2.3.5 2.3.6 2.3.7
Input- and output weightings . . . . . . . . . . . . . . . . . . . . . . . The filter gain
20
21
...............................
A constraint window . . . . . . . . . . . . . . . . . . . . . . . . . . . . Control blocks
23 24 25 25 27 27
............................... .............................
...........................
2.4 Alternative MPC methods 2.4.1 2.4.2
Repeated linearization Nonlinear MPC
............................... .................
2.5 MPC implementations: Matlab and PRIMACS
3 The gas turbine installation
30
3.1 Construction and operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The compressor
30 31 35 36 37 39 39 40 42
...................................
3.3 The combustion chamber 3.4 The turbine
.............................
.....................................
3.5 The overall setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Operation and control 3.6.1 3.6.2 3.6.3
...............................
Inputs and outputs Constraints
.............................
.................................
Control objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 Modeling the gas turbine
44
4.1 The compressor model
...............................
44
4.2 Extensions t o the basic compressor model
....................
46
48 54
4.3 Modeling the complete gas turbine . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
M P C applied on the gas turbine: simulations
57
57 57 59 59 60 61 62 62 64 66 67 67 68 69 69 70 70 71 71 72
5.1 The MPC implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 5.1.2 5.1.3 5.1.4 5.1.5 5.2 Scaling
...................................
Discretization Reconstruction Constraints
................................
...............................
................................. .............................
Parameter settings
Simulation results 5.2.1 5.2.2 5.2.3 5.2.4 5.2.5 5.2.6 5.2.7 5.2.8 5.2.9 Simulation 1 Simulation 2
.................................
................................
................................
Simulation 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 6
................................
Simulation 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.10 Simulation 10
................................
5.2.11 Simulation 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.12 Unfortunately . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Experiments with PRIMACS
75
7 Conclusions & Recommendat ions
78
7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Recommendations
78 79
81
.................................
Bibliography
A MPC
A .1 Unconstrained Model Predictive Control . . . . . . . . . . . . . . . . . . . . . A.2 Constrained Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . A.3 The filter gain
83
83
85
88
88
.................................... .............................
A.4 Matlab versus PRIMACS
B Gas turbine model specifics

B.l The Garrett configuration B.2 The BBC configuration
90
.............................
90 93
..............................
C MPC implementation specifics

D MPC parameter settings
100
103
Only those symbols are included that are frequently used. We divided them into two groups being MPC (Model Predictive Control) related symbols on the one hand, and gas turbine related symbols on the other. If any symbol appears in both categories, it is assumed not to be too easily confused.
MPC
IM IOM IP M QP
k
U
Internal Model Internal Optimization Model Internal Prediction Model Quadratic Program sample counter inputs (or: manipulated variables) states state estimates outputs output estimate at k 1, based on measurement information up to k output measurements outputs before the addition of measurement noise output references unmeasured disturbances measured disturbances measurement noise offset between process and model manipulated variable moves prediction horizon control horizon number of samples after which the system settles MPC controller gain first-order filter gain Kalman filter gain projected error vector (z - y for zero future inputs) T output- and input weights step response matrices accounting for effects of U and d, respectively
X
X
Y Y@
Ym
+W )
Y
r d
20
z 4
AU
P m n
KMPC K p
JcKalman
E P
ry, r u
SU, Sd
corrected vector of future outputs, defined by (A.8) state-space realization discrete matrix connecting z ( k ) to z ( k 1) discrete matrix connecting u ( k ) t o z ( k i),not t o be confused with input weightings r" discrete matrix connecting d ( k ) to z ( k -t 1) discrete matrix connecting w ( k ) to z ( k 1) discrete matrix connecting u ( k ) to y(k) discrete matrix connecting d ( k ) t o y(k) discrete matrix connecting w ( k ) to y(k) output prediction vector, truncated at k p output reference vector, truncated at k p matrix connecting Y with Y manipulated variable moves vector, truncated at k p low and high constraint on u, similar definitions for y unitary matrix and matrix with unitary matrices on the lower triangular
+ +
+
+ +
Gas turbine
P
m
T V A L N
rl
pressure [Pa] mass-flow [kg/s] temperature [KI (plenum) volume [m3] (duct) area [m'] (duct) length [m] rotational speed [rev/s] efficiency [-I
subscripts
ss comp bl t hr tb in out rPm
T fuel bar max Pl cc amplitude offset
steady state compressor blow-off throttle turbine at inlet at outlet revolutions per second at high temperature fuel property pressure in [bar] maximum value plenum combustion chamber amplitude of signal offset of signal
S
SB ST
sv
DEVSurge
c,
s
F
C p , Cv
valve position E [O, 1 [-I 1 blow-off valve (position) throttle valve (position) fuel valve (position) deviation from the surge line [kg/s] steady state compressor pressure rise [Pa] pressure drop over valve [Pa] specific heat capacities at constant pressure and temperature, respectively [J/kg K]
CP/CU
Y
R
[-I
gas constant [N m/kg KI P(n) (normalized) air density [kg/m3] a speed of sound [m/s] U impeller tip speed [m/s] WH Helmholtz resonance frequency [Hz] A/AComp : duct area relative t o compressor duct area [-I 4 L/LcOmp: length relative t o compressor duct length duct G n p ,A P pressure drops [Pa] 7time lag [SI I inertia [kg m2] B Greitzer B parameter [-I e n in-going power [W] WP requested power [W] (overall) friction factor [-I Scale scaling factor, defined in Appendix B Acomp, Beonp, Ceonp coefficients of compressor characteristic scaling factors for state i , output i and input i ei, cy,%,e u , i scaled x dimensionless x
[-I
2
x
Summary
In this project, we investigated the use of Model Predictive Control (MPC) on a gas turbine laboratory installation. This is part of the research into modeling and control of compressor and turbine systems, done by the sections Energy Technology and Systems & Control from Mechanical Engineering (TUE), in co-operation with TNO-TPD (TNOs Technisch Physische Dienst) at Delft. In a gas turbine, a compressor pressurizes gas, which after being heated in a combustion chamber, enters a turbine where pressure ratio is converted t o rotational speed of the shaft that turbine and compressor share. This way, the turbine drives the compressor by providing shaft power, and the compressor in its turn drives the turbine by delivering gas under pressure. The addition of heat in the combustion chamber enables the gas turbine to supply power to some external load. However, the system under consideration, which is a scale model of a gas turbine, serves as a test-stand for research purposes-not as a power supply. The gas turbine is limited in its operation by several undesirable effects, imposed by the (combination of) components it comprises. The two most important limitations are swge, and the maximum turbine inlet temperature. Surge denotes the flow-instability associated with a sudden drop in compressor delivery pressure and with violent dynamic pulsation which is transmitted throughout the whole machine. It is t o be prevented at all times, since it can be most destructive. Next, the turbine inlet temperature is not allowed to exceed a certain maximum, beyond which the stress exerted on the turbine rotor-blades becomes undesirably high. These limitations are generally referred to as constraints. Other than these restrictions, actuators impose constraints as well. For instance, valves cannot be opened beyond fully opened, or-just as obviously-close beyond fully closed. Furthermore, they can not move at arbitrarily high speed. This is why move constraints will have t o be specified as well.
It is especially these constraints, which have resulted in choosing MPC for a control strategy. Unlike most other control algorithms, MPC can explicitly incorporate constraints in its computations. Shortly put, MPC computes its controller moves by minimizing the difference between desired and predicted output response, while satisfying all constraints. Outputs are predicted by using a linear model of the process.
In this project we first obtained a model of the gas turbine installation. This model has t o be
9
reasonably accurate on the one hand, since it largely determines the performance of MPC. On the other hand, computational demands require it t o be relatively simple. Not without reason has MPC so far been applied on comparatively slow (e.g. chemical) systems, mostly. Controller execution time will therefore be of crucial importance. Next, we simulated this model using Matlabs MPC Toolbox, but especially TNOs PRIMACS. Eventualiy, we sought t o apply MPC in a real-time implementation, which unfortunately wasnt really accomplished, since time did not permit to implement MPC in a closed-loop situation. However, we did get a chance t o use MPC open-loop, which means that we applied inputs, obtained from a certain simulation, to the laboratory installation and compared the process response to the model response. The results of this seem t o be rather promising, as will be seen in Chapter 6, although the model we used might still be improved. Furthermore, it is difficult t o say what will happen if MPC is eventually used in closed-loop, since especially measurement noise is something which must not be underestimated. Still, based on the results thus far, we believe that computational demands will not render real-time implementation impossible.
Chapter 1
Introduction
In co-operation with TNO-TPD (TNOs Technisch Physische Dienst) at Delft, the sections Energy Technology and Systems & Control from Mechanical Engineering (TUE) are researching the modeling and control of compressor and turbine systems. Compressors are often to be found in process-industrial applications. and essentially provide a certain mass-flow under some specified pressure. The compressors operation is unfortunately limited by some undesirable effects. First of all, surge is the flow-instability associated with a sudden drop in delivery pressure and with violent dynamic pulsation which is transmitted throughout the whole machine. It is t o be prevented at all times, since it can be most destructive. Secondly, choke is the phenomenon that limits the mass-flow delivered by the compressor for some constant rotational speed. In this project, we consider the compressor as a component of the gas turbine.
1.1
Gas turbines
In a gas turbine, a turbine expands pressurized gas, basically converting pressure ratio t o rotational speed. The shaft power so obtained then drives a compressor, which in its turn provides the pressure ratio needed for the turbine. If the installation were without losses of any kind, it would be able to drive itself, but of course no more than that. In practice, the losses that are always present have t o be overcome which is why energy has t o be added one way or another. In our case, this is accomplished by the incorporation of a combustion chamber, which heats the pressurized gas coming out of the compressor before it enters the turbine (see Figure 1.1). If more energy is added than necessary t o overcome losses, the gas turbine can drive some external load, for instance, a generator. In the system under consideration, we do not make use of this possibility. Instead, the installation we investigate, which is a scale model of a gas turbine, primarily serves as a test-stand for research into dynamic behaviour of compressors
11
12
C H A P T E R i. INTRODUCTION
Figure 1.1: The gas turbine in its simplest form and turbines and their interaction-not as a power plant. This laboratory installation uses air for working fluid and natural gas for fuel. Next to surge and choke, gas turbine operation is limited by the turbine inlet temperature not being allowed t o exceed a certain maximum, beyond which the stress exerted on the turbine rotor- blades becomes undesirably high. Apart from these output constraints, limits on actuator operation form input constraints. For instance, valves cannot open beyond fully opened or, just as obviously, close beyond fully closed. Since we wish to control gas turbine operation, for instance, operate at constant rotational speed for varying loads without exceeding any of these constraints, the choice of a control strategy will very much be determined by the strategy's capability in handling constraints.
1.2
Model Predictive Control
It is these constraints which have resulted in choosing Model Fredictive Control (MPC) or a control strategy. MPC is a control methodology, which explicitly includes constraints in its optimization. The vast majority of widely used control algorithms does not exhibit this important feature, and is therefore hardly equipped to handle constraints1. MPC uses a model of the process which is to be controlled, to predict system response over a future (prediction) horizon. Through minimizing the difference between desired and predicted output response, the inputs (or manipulated variables) can be determined. The use of a prediction horizon also has the benefit that future setpoint (desired output) changes or known disturbances can be
'Most control algorithms can only indirectly influence-not guarantee-variables to meet their constraints, by using ad hoc methods, mostly. To some extend (problem must be feasible N solvable), MPC can in fact guarantee variables not to violate their limits.
1.3. GOALS OF THIS PROJECT
13
anticipated for. Apart from this MPC can control multivariable systems, which may exhibit nonlinear behaviour, interaction as well as dead-times. All in all, MPC seems t o offer benefits especially useful to the control of the gas turbine, which indeed exhibits highly interactive and nonlinear behaviour. So far, though, MPC has been applied mostly on comparatively slow systems, such as chemical processes-and not without reason. Mechanical systems, such as the gas turbine, are characterized by much smaller timeconstafits, causing computing-time related problems that must not be underestimated. MPC is a control algorithm, which requires relatively large execution times, especially if many constraints are present. Furthermore, it remains t o be seen if linear MPC can cope with the highly nonlinear behaviour we know is present.
1.3
Goals of this project
This constitutes the main goal of this project. We attempt t o assess the feasibility of MPC on the gas turbine installation. We do expect problems concerning execution times and nonlinearity and therefore, the problem statement can be posed more specifically as: Find out what can be achieved with MPC, using a model as simple as needed to allow for real-time implementation, and, to what extend must we compromise in order to do so? Having done this, do the advantages we claimed MPC t o possess over other control strategies (for this specific case, that is), still hold? The final question will be if the resulting closed-loop behaviour is acceptable, and furthermore, if it can compete with other, more common, control strategies3. Most probably, MPC will have to give in somewhat on performance, in exchange of the desired constraint handling. Secondary t o this main problem, we will concentrate on testing (and possibly improving) TNOs MPC-implementation PRIMACS, which still is subject of constant development.
1.4
Project setup
In order t o arrive at these goals, we first developed a simulation model as a compromise between simplicity (minimizing computing-time) and accuracy (maximizing usefulness of prediction). Next, control-goals and -strategy were specified after which the MPC-setup could be developed and implemented in standard packages being Matlabs MPC-Toolbox (by Morari and Ricker) and TNOs PRIMACS. The following step was to tune MPC- and filter-parameters. Before the final experiments could take place, further development of the simulation-model was carried out t o approximate reality as good as practically possible (use only measured outputs for feedback, introduce measurement noise, deliberately create differences between simulation-model and prediction-model to test for robustness, etc.). Finally, open-loop ex2This is the system which combines process and controller. 3The comparison with other control strategies actually is not part of the present assignment, but definitely of the total project.
14
C H A P T E R 1. INTRODUCTION
periments took place using PRIMACS and LabVIEW, where inputs were computed from simulations off-line. Unfortunately, time did not permit for the real-time closed-loop experiments t o be carried out. This, and the comparison with other control algorithms will have to be carried out by a successor. Next, some related previous research using PRIMACS should be mentioned. First of all, we mention the application of MPC under PRIMACS on the FCC (Fluidized Catalytic Cracker) process by Peeters [Peeters 951. This was the first project to take place within the co-operation of TNO-TPD and TUE. Secondly, an application of MPC using PRIMACS on a Thermal Hydraulic Process was examined by Frits van der Meulen [Van der Meulen 951, whose report we will refer t o more than once. More closely related to this project is the research into an MPC-application on a compressor-station of two radial compressors in parallel by Inge Satter [Satter 961. Such a configuration serves to increase capacity, whereas a serial construction would allow for larger pressure ratios. Combinations of these configurations can be imagined as well (see, e.g. [Van Essen 95bl). An additional reason for employing multiple compressors may be t o improve reliability. Further attention to the parallel configuration is paid by Gert Leenheers [Leenheers 971 in an ongoing project which should be finished somewhere in March 1997. The latter two projects, and the present project as well, all arose from the PhD-project by Harm van Essen, who is researching the modeling and control of compressor related systems, primarily making use of model based control strategies.
1.5
Organization of report
The report is organized as follows: In Chapter 2 the theory behind MPC will (to some extend) be explained, along with a (brief) description of two practical implementations being Matlabs MPC-Toolbox and TNOs PRIMACS. Chapter 3 describes the gas turbine installation, whereas in Chapter 4 a lumped-parameter model is developed. Then, in Chapter 5, simulations in PRIMACS are presented followed by a discussion of the results thus far. Results of real-time experiments on the laboratory installation are the subject of Chapter 6 and, finally, Chapter 7 closes with overall conclusions and some recommendations for future work in this area. Finally, we should mention the fact that in the gas turbine installation, we used two different compressor/turbine combinations, i.e. Garrett and BBC, respectively. The reason for this is that, although originally the installation was designed for the BBC configuration, at some time it was used as a test-stand for a turbo charger which can be found in trucks. During this project, the installation was reconfigured t o its original shape. Although some attention will be paid t o the different characteristics, this will not be vital t o our problem statement.
Chapter 2
Theory: MPC
2.1
Unconstrained Model Predictive Control
We will start with describing MPC in its simplest form which is linear unconstrained MPC. Only in Section 2.4.2 we will slightly go into nonlinear MPC algorithms, since our interest lies in linear MPC, being the least time-consuming method of the two, and execution time will be of crucial importance. Model Predictive Control (MPC) utilizes a model of the process which is to be controlled to predict the outputs up to a certain time-instant, based on the inputs t o the system and the most recent feedback values. It uses this model prediction t o compute the moves in the manipulated variables that will make the system perform in some desirable manner over a future time horizon. The model used by MPC for predicting output response and computing manipulated variable moves will be referred t o as Internal Model (IM). For example, consider Figure 2.1. At the present time k , the response of the output y(k) t o changes in the manipulated variables a ( k ) is predicted over the prediction horizon p using the IM. The manipulated variables are allowed t o vary over the control horizon 'M. and can be computed such that future deviations of the controlled variables from a desired target r ( k ) (setpoint) are minimized. Of the computed optimal control moves, only the values for the first sample are actually implemented. This way the most recent feedback values can be used t o calculate a new sequence of control moves of which again only the values corresponding to the first sample are used. This mechanism is also known as a moving (or receding) horizon. Normally, MPC uses the IM both for predicting outputs and computing inputs. However, in some cases, these models may differ. For instance, linear MPC always uses a linear(ized) model for computation of moves, but the model used for prediction is sometimes seen t o be nonlinear. We will pay more attention t o this particular approach in Section 2.4.2. In such a
15
16
CHAPTER 2. THEORY: MPC
Figure 2.1: The basic concept of MPC case, we will refer to the model used for prediction as Internal Prediction Model (IPM) and the model used for optimization as Internal Optimization Model (IOM).
2.1.1
The optimization
The wish to minimize future deviations of the controlled variables from their setpoints must now be expressed in mathematical terms so that a control law can be obtained in algorithmic form. A commonly used criterion in MPC is the quadratic objective, which, when combined with the wish t o prevent the inputs from becoming inadmissibly large, can be stated as
where weights are included to express the relative importance of outputs following their reference trajectory on the one hand and trading this off with reducing the action of manipulated variables on the other. In this notation : and I represent the output- and input-weightings, ' I ' ? respectively. Furthermore, y(k -k 1 I K ) denotes the estimate of y(7 I) obtained at k , taking into account al measurement information up to 7. Au are the manipulated variable moves. l
The solution to this optimization is obtained in Appendix A and will only be stated here as the resulting control law
Au(-) = K M P C E P ( k
+ 117)
where Ep(7 117) is the measurement corrected vector of future output deviations from the reference trajectory (or projected error vector), assuming all future moves are zero. Simply Put, EP = reference T - output y (for zero future moves).
2.2. CONSTRAINED MODEL PREDICTIVE CONTROL
17
K M p c can be computed off-line when EP(- i l k ) is the only time varying element, which see normally is the case. Only when weights (I?, r.) or model (SU, Appendix A) need t o be updated, K M p C must be recomputed. The latter is the case, for instance, when repeated linearization is used, which will be discussed in Section 2.4.1.
2.2
Constraints are always present in any real life process control situation. To include these in the optimization problem of MPC, the control formulation is expanded, to arrive at an algorithm which finds the best solution satisfying all constraints. There are three types of process constraints: Manipulated variable constraints: hard limits on the inputs u ( k ) t o account for, e.g. valve saturation constraints. Manipulated variable rate constraints: hard limits on the size of manipulated variable moves A u ( k ) to directly influence the rate of change of the manipulated variables. Output variable constraints: hard or soft' constraints on system outputs that may not exceed some minimum or maximum value; these outputs can, but needn't be, controlled variables with setpoint t r aject ories .
Incorporating these constraints eventually leads t o a Quadratic Program which must be solved every controller executing time. More details can be found in Appendix A and [Morari et al. 911.
2.3
Tuning parameters
In this section the effects of tuning parameter settings are discussed. These tuning parameters are:
1. the sample time

2. the prediction horizon p
3. the control horizon m
4. the input weightings I?"

'The exact meaning of hard and soft constraints will be defined in Section 2.5. For now, we will say that soft constraints may be (temporarily) released in case of an infeasible optimization problem, whereas hard constraints may never be released.
18

5. the output weightings Y
6. the filter gain K 7. a constraint window
8. control blocks
2.3.1
The sample time
Probably the first tuning parameter that must be chosen is the sample time. The sample time is defined as the time interval between two successive moments at which MPC computes manipulated variable moves. Since the same sample intervals are used for discretization and computation, the sample time can not be chosen too large, since doing so would result in a highly inaccurate (if not unstable) model response. On the other hand, choosing the sample time too small greatly increases execution time, while the discretization may not get considerably more accurate. Therefore, we have t o choose a compromise between accuracy and execution time.
If fast system dynamics force the discretization sample time to be rather small, real-time implementation may become an impossibility. When this is the case, one might decide t o only consider slower system dynamics, if this is indeed acceptable. As will be seen in Chapter 4, we will in fact choose t o ignore faster system dynamics, being mass-flow response, by using static instead of differential equations to compute mass-flows. Apart from the sampling time, execution time is also determined by some of the tuning parameters that will follow.
2.3.2
The prediction horizon
The prediction horizon is defined as the number of samples into the future for which output variables are predicted. When choosing this parameter, several issues must be taken into consideration. Roughly taken, the following rules of thumb exist:
o
the prediction horizon must exceed the inverse response period, when this type of behavior is indeed present [Satter 961; the prediction horizon must exceed any dead-time periods present [Satter 961; the prediction horizon must be at least as large as the largest time-constant corresponding to quantities t o be controlled; increasing the prediction horizon results in
-
o o
less aggressive control actions, slower system response, and a more robust closed-loop behaviour,
2.3. TUNING PARAMETERS and vice versa;

o
19
the prediction horizon would ideally be chosen as large as possible, thus being more likely t o yield nominal stability of the closed-loop system; the prediction horizon should on the other hand be chosen as small as necessary for real-time implementation.
Some of these points need further explanation. The first of these is rather obvious: if p is smaller than the inverse response period, this inverse response is all MPC will ever see in its prediction. So the moves MPC computes will always be the inverse of what eventually would make the system respond in the desired way. Next, choosing p smaller than the dead-time period will result in MPC having absolutely no idea of what t o do. It will simply see no response whatsoever (in the extreme case of all outputs exhibiting dead-times) and might therefore conclude not to do anything, since the deviation of outputs from their references can not be reduced within one prediction horizon2. Similarly, it can be imagined that p must exceed the largest time-constant that is associated with outputs we wish t o control. Choosing p much smaller than this will cause MPC t o overlook long-term effects of its moves. Since choosing p smaller results in MPC to become rather short-sighted, it will behave more aggressive; it just doesnt see the long-term effects of its actions. This will then cause the system t o respond faster, although the effects of model mismatch and disturbances will become ever more evident-in other words the closed-loop behaviour is less robust. For open-loop stable processes, nominal stability of the closed-loop system only depends on KM,, (for the unconstrained case) which in turn is affected by p , m, and rp. No y precise conditions on p , m, and ;l exist which guarantee closed-loop stability. Increasing y p relative t o m tends t o stabilize the system, apart from making the controller less aggressive. For p = co closed-loop nominal stability is guaranteed for any finite m and time-invariant input- and output weights [Van der Meulen 951. Since the control scheme will have t o be implemented real-time eventually, computing-time becomes an important issue. Therefore, it might be necessary t o take p equal to the largest time-constant present, if in fact (and this is usually the case) it is the time-constant which imposes the minimum value on p . When acceptable, it might even be necessary t o only consider the system behavior related to smaller time-constants, if the main interest lies here. Alternatively, if slower system dynamics are considered more important, one can decide t o take a larger sampling period. Of course, this has its limits as well, since the discretized linear model of the process (which is used as Internal Model) might become highly inaccurate or even unstable for too large a discretization sampling period, as was discussed in the previous Subsection.
2So any moves
> O only
increase the optimization criterion.
20
C H A P T E R 2. THEORY: MPC
2.3.3
The control horizon
The control horizon is the number of samples into the future for which optimal control moves are computed. As for the prediction horizon, this tuning parameter can be chosen using some basic rules:
the control horizon must be at least one sampling period and can never exceed the prediction horizon, or 1 I m 5 p ; : in practice, the control horizon is chosen to be 1 / 6 t o 1/3 of the prediction horizon [Van der Meulen 951; increasing the control horizon results in
-
more aggressive control actions, faster system response, and less robust closed-loop behaviour,
- a
and vice versa;

o
the control horizon should preferably be chosen as small as possible in order to reduce computing time.
The first of these points is rather trivial: of course, choosing m = O cancels out the controller and, just as obvious, taking m > p is useless, since MPC cant see beyond the prediction horizon. Varying the size of m has effects inverse to those of varying p . If we choose m quite small, say m = 2, MPC will become rather tame since it cannot afford t o mess up with only 2 input moments per prediction horizon available. However, if m is much larger, the first input moves can be rather big, since the inputs that follow can take care of slowing the system down again. Again, one should weigh these recommendations for each specific case. If, for instance, faster system response is desirable, enough computation-time must be available or otherwise the control horizon can no longer be enlarged. It is, in fact, especially the control horizon that influences execution time because with it the number of degrees of freedom in the Quadratic Program (which has to be solved each controller execution time) increases proportionally.
2.3.4 Input- and output weightings

We treat both weightings in one subsection, since it is only the ratio between the two that actually matters. Once output weightings are determined such that relative importance of the several outputs is reflected, input weightings can be chosen freely in order to effect whatever closed-loop system behaviour is desired. Decreasing input weightings basically has the same
2.3. TUNING PARAMETERS
21
effects as increasing the control horizon, with the notable difference that input weightings do not affect execution time. : is most commonly used as tuning parameter. Furthermore, ' I when ; increases, the influence of m on closed-loop system behaviour decreases. ' I
2.3.5
The filter gain
A filter is used t o correct model information for process information. Model and process differ because of model mismatch on the one hand and system- and measurement noise on the other. System noise denotes disturbances (noise) that directly influence the process, i.e. actuator signals and states. Measurement noise is the noise that is present on the output measurements.
To correct for these undesirable effects, we use a filter. Depending on the fact if the filter is model-based or not, we can divide filters into two groups, starting with the non-model-based type1. Non-model-based filtering
This type of filter only uses model outputs and process output measurements-not the model itself. An example is the first-order filter K F , used as the general filter in standard linear MPC as described in [Morari and Ricker 931, and stated here in its simplest form
where y,(k) is the measurement at time k and y ( k l k - 1) denotes the estimate of y(L) obtained at k - 1, taking into account all measurement information up t o k - 1. For more information of the filter as used in [Morari and Ricker 931, the reader is referred t o Appendix A.3. Choosing the filter gain relatively large would in theory best compensate for the differences between model and measurements. In fact, if we choose the filter constants equal t o 1, we have a zero-order filter, where the corrected model output exactly equals the process output measurement. This also means that any noise present in the process measurements directly enters the model, which is why filter settings should not be too tight. Furthermore, if the r e d system only j n s t exceeded some constraint, filter related actions can, so t o speak, 'pull' the Internal Model over the constraint, beyond which the QP problem might become infeasible, whereas the uncorrected model most probably wouldn't have experienced any severe problems. For example, if in case of the gas turbine, the turbine inlet temperature would exceed its maximum when already near surge, filter corrections can cause the model temperature to exceed its constraint as well, after which the QP problem may become infeasible if the controller has no options (inputs) left t o immediately return t o a feasible situation. Situations like this have been seen to occur more than once. If filter gains, on the other hand, become too small, the control algorithm barely makes use of measurements and open-loop control is almost what we get.
22 2. Model-based filtering
This type of filtering uses process measurements, model outputs and the model itself to determine corrective actions. An example of such a filter is the Kalman filter. This filter can be applied only if the IM is linear and secondly, if noise can be assumed to be zero mean white noise '. Let's consider a (discrete) state-space representation of the process (in its simplest form)
where z ( k ) denotes the model states and u ( k ) and y(k) the inputs and outputs. Assuming that system noise s ( k ) and measurement noise z ( k ) are present4, and furthermore that D = O, which is rather common, the observed system can be represented by
Reconstructing the states can then be done using the optimal Kalman filter
with IC(k1k - 1) the estimate of .(IC), t o k - 1.
obtained at k
1, and based on measurements up
The optimal Kalman gain for the stationary filter follows from the discrete algebraic Riccati equation
Q = V,
KKalman
+ Q [Q-' + C'K-'C]-' QCT [V,+ CQCT]-'
a>*
where V, and V, are intensity-matrices for s and z respectively. A stabilizing, unambiguous solution to this problem exists if and only if Q = QT 2 O. Still, if offsets between process and model exist, the assumption that the system noise s has zero mean is no longer valid. In that case, the standard Kalman filter will not suffice, which may lead to offsets in the controlled outputs if the filter is not modified in some way. Instead, it is recommended to extend the standard Kalman algorithm with a disturbance-model, in order t o get an estimate of the offsets. Such an extension can be an output-disturbance model, an input-disturbance model or, even, a combination of both (apart from variations like the 'velocity-form', see [Van der Meulen 951).
~~
3This is white noise with a mean value of zero. White noise is a mathematical fiction, which denotes a signal that is totally unpredictable, i.e. one could never predict a future value based on information up to the present time sample. Also see [Kok 901. 41n [Kok 901 system noise is denoted by w ( k ) , whereas u ( k ) is the measurement noise. Since we will use these symbols in other places, we have chosen t o use the definitions stated in the above.
2.3. TUNING PARAMETERS
23
The first case of these denotes the situation where the states are augmented with stepwise disturbances q ( z in [Van der Meulen 95]), which are then simply added to the output y. This way we get
[Van der Meulen 951 takes a closer look at this. The input-disturbance model assumes that q ( k ) has a direct effect on z ( k ) instead of y(k). In state-space notation this yields
x(kk + 1)) ] *( 1
[A o
G ] I
[ ;[:i]+[
+ 1).
::]u@)+[
where G is introduced to connect q ( k ) to z ( k
Most of what was treated here was taken from [Van der Meulen 95, Muske and Rawlings 931.
2.3.6
A constraint window
The constraint window denotes a time interval within the prediction horizon, for which output constraints are included in the optimization. Outside of this window, output constraints are simply ignored. Figure 2.2 depicts this situation. In PRIMACS, the constraint window is assumed t o start at a pre-specified sample E [ l , p ] and stop at the end of the prediction horizon5. If the constraint window starts at sample 1, it equals the prediction horizon. This will be referred t o as not using a constraint window, although this may sound somewhat paradoxical, since such a constraint window actually is the largest one that can be used. In Matlab, one can simply specify constraints t o be valid for each sample independently. For instance, if p = 15, one might decide to specify constraints for samples 3 , 4, 5, 11 and 15, or any other combination. Not all combinations seem useful though, which is why PRIMACS assumes the constraint window t o be continuous. Still, it might be useful t o leave outputs unconstrained for some samples in between of starting and ending point of the constraint window in order t o reduce computation effort, which highly depends on the total number of inequality constraints. Of course, this can only be done if the system can be assumed not t o oscillate violently within this unconstrained time interval.
51n the following we will use PRIMACSs definition of a constraint window.
24
CHAPTER 2 THEORY: MPC .
p e d i c(e2 out p u t
I
k
k+ P
Figure 2.2: The constraint window. Here a situation is depicted, where the predicted outputs do in fact exceed their limits outside of the constraint window. Still, the prediction at sample k , which can be assumed to be reasonably accurate, does not exceed the constraint Especially for processes exhibiting dead time behaviour, the use of constraint windows is very much recommended, if not obligatory. The starting point of the constraint window should correspond to a number of intervals larger than the dead time. But also in other cases a constraint window can be very useful, as will be seen in Chapter 5. Without a constraint window, the QP problem becomes infeasible if during operation outputs temporarily exceed their constraints. Depending on the implementation in question, constraints are either released or, even worse, the controller stops altogether because of infeasibility. If, however, a constraint window is used, MPC is allowed to ignore constraints during the first few samples of the prediction horizon. Although this can lead t o some overshoot, it is preferable t o the two situations mentioned earlier, being constraint-release and controller-stop. Too big a constraint window might still leave the modified QP problem infeasible with the same consequences as before. Decreasing the size of the constraint window can then solve the problem, although eventually this might lead to a situation where constraints are never (guaranteed to be) satisfied, which is no better than releasing the constraints altogether. Still, when the situation occurs that constraints are about to be released, but wouldnt be with a slightly smaller constraint window, it might be recommendable to decrease the constraint window. Releasing constraints is in fact almost the same as using a constraint window of zero size. Frits van der Meulen recommends that this procediire is automated, so that the maximum is done to keep a constraint from being released. More on this subject will be seen in Chapter 5 and also can be found in [Van der Meulen 951.
2.3.7
Control blocks
Finally, control blocks can be used, meaning that the manipulated variables can no longer vary from one sample to another, but only from one control block to the other. For this purpose, the control horizon is divided into blocks of possibly unequal size, the size varying from 1
2.4. ALTERNATIVE MPC METHODS
25
t o m samples6. Thus, during one block, the manipulated variable moves Au are taken t o be zero during the optimization, thereby greatly reducing execution times. Actually intended to prevent the controller from ringing (or oscillating)) control blocking is often more effective in doing so than the alternative of choosing the prediction horizon significantly larger than the control horizon. This of course at the expense of a more sluggish system response. Control blocking is certainly not the same as choosing a larger sampling period for MPC, since doing so would also lead t o a far less accurate (if not unstable7) output prediction. Another reason for incorporating control blocks might be t o improve anticipation abilities) without increasing aggressiveness of the controller. For example, if the same number of blocks is used but with blocks twice as large, the controller can (directly) influence the process twice as far into the future without becoming more aggressive and without increasing controller execution time. Summarizing) the size of control blocks should be chosen as a compromise between sluggishness and aggressiveness, in conjunction with possible anticipation demands.
2.4
Alternative MPC methods
In this section we summarize some approaches other than standard linear MPC, which was discussed so far. We will conclude with a brief discussion of nonlinear methods, but first describe approaches which lie somewhere between linear and nonlinear MPC, where nonlinear MPC denotes the form that uses nonlinear models both for prediction and optimization.
2.4.1
Repeated linearization
We will start with the intermediate form that still uses a linear model both for prediction and optimization, but this time obtained at more than one sample, which normally is the case (linearizing only once, that is). This means that a nonlinear model is linearized once every few samples, or ultimately) every sample-but still not within the prediction horizon. Repeated linearization can be very useful, especially if the actual process exhibits strong nonlinearities-which indeed the gas turbine does-since process and linear model will deviate more and more as we move away from the original linearization point. Therefore) it seems worthwhile t o repeat linearization every time the difference between process (which is thus far still represented by a nonlinear model) and linear model exceeds some prespecified bound. This already describes the first and most clever) one of the four linearization criteria that
~
Not using control blocking in this sense is equivalent to using m blocks one sample wide. When speaking of control blocks, at least one block should be two samples in size. 7See the remark on sampling periods in Subsection 2.3.2. 81f we linearize only once, we will do this around the starting conditions. Of course, one could just as l) well choose to do this around the desired ending position (if such a position can be distinguished at al, or somewhere between beginning and end, even.
26 will follow.
CHAPTER 2 THEORY: MPC .
But before that, it must be mentioned that repeated linearization also takes its time t o be executed. Linearizing the model means that K M p C has t o be recomputed, since the entire system has changed. Without repeated linearization, K M p c had to be computed only once (unless weightings were changed) which could take place off-line, so laborious computations were not really a problem. When repeated linearization is used, it becomes eminently clear that recomputation of K i ~ P C forms the major part of computation efforts. However, we wonldnt use it if we didnt expect it to offer great improvements in closed-loop performance. Furthermore, repeated linearization requires a reasonably accurate knowledge of all states in order t o be effective. This means that states that arent available as measurements will have t o be reconstructed one way or another. Usually this is done using a Kalman filter (see Section 2.3.5), but different methods can be imagined as well. In Subsection 5.1.3, we wio pay attention t o this problem once more. Four possible moments to linearize are:
1. Linearize whenever the deviation between process and linear model reaches some prespecified bound S. We measure the deviation by taking the maximum absolute difference between model output (or, actually, state) and measurement. Of course, when not all states are measured, some will have to be reconstructed from those states which are measured.
2. Linearize every sample: this is the most time-consuming option of all. On the other hand, this is also the simplest option t o implement [Haarsma 951.
3. Linearize every 2 samples: the second simplest option, where execution time can be reduced (compared to option 2, that is) by choosing 2 > 1.
4. Linearize whenever new setpoints are specified [Haarsma 951. This option may result in only one linearization if setpoints are given at the beginning of the movement only, for instance, in a point-to-point trajectory.
To perform the linearization, information on partial derivatives of the differential equations will have t o be available to MPC in one form or another. Partial derivatives can either be obtained numerically or analytically. Although the first option has the small benefit that no software besides PRIMACS will have to be used, analytically obtained symbolic expressions can simply be evaluated at each new operating point. Such expressions can be obtained either by hand, or, for more complex cases, with symbolic computation packages such as Maple [Maple 941.
If repeated linearization is used within the prediction horizon, one can still distinguish between using this newly obtained model both as IPM and IOM, or just one of the two. It is not known whether such approaches exist, or are even possible at all, since optimization with a changing (different for each sample) model can be expected to pose considerable problems.
2.5. MPC IMPLEMENTATIONS: M A T L A B A N D PRIMACS
27
2.4.2
Nonlinear MPC
In this subsection we mention methods which use nonlinear models, whether this occurs during prediction, optimization or both. The first t o be mentioned is the approach that uses a nonlinear IFM combined with a linear IOM. The underlying premise is that the superposition theorem-which is valid for linear systems-still holds to a reasonable extend. Thus reasoning, nonlinearly predicted effects of initial conditions can be added t o the effects of manipulated variable moves obtained from the optimization, much like is done in Appendix A.l, albeit in this appendix both effects are computed with the same linear model, of course. For instance, [Satter 961 uses this intermediate approach and does so with what appear t o be rather promising results. It must however be investigated if this method can be guaranteed not to cause any unforeseeable problems. It would be better still if the method somehow could be guaranteed t o improve controller performance at all times, since on the other hand it will increase controller execution time at all times. Combinations of nonlinear prediction and repeated linearization-whether or not within the prediction horizon-can be imagined as well; all of this t o stay away from nonlinear optimization, which requires computationally very demanding techniques. Some of these techniques use iterative optimization until all constraints are met, while others substitute constraints into the nonlinear differential equations. One can imagine that both approaches require considerable computational efforts, which can more or less be said for all nonlinear methods. Apart from that, it remains to be seen if convergence is guaranteed. In the area of nonlinear MPC, no clean solution has been found yet, which is why it remains the subject of continuous development. For more detailed information on nonlinear MPC, the reader is referred to, e.g. [Bequette 91, Van Essen 95~1.
2.5
MPC implementations: Matlab and PRIMACS
In this section we will take a closer look at two MPC implementations, being Matlabs MPCTools [Morari and Ricker 931 and TNOs PRIMACS. That is, no exhaustive description of these implementations is pursued, but typical features and striking differences will be discussed. First of all, it should be mentioned that PRIMACS was developed t o improve real-time implementation on practical applications, over other implementations such as Matlabs MPC Toolbox. PRIMACS can receive measurement data from, and send input signals to any process, whether this is a linear model, nonlinear model or real-life process. This is right away one of the most significant differences with Matlab: Matlab uses linear models both for the Internal Model and for the Simulation Model. These two linear models can, by the way, still be different t o simulate plant/controller model mismatch. However, if nonlinear simulation models
28
C H A P T E R 2. THEORY: MPC
are t o be used in Matlab, no straightforward methods are available to do so9. Let us consider the following discrete-time linear time-invariant representation-which ally is the one used by Matlab: actu-
In this representation d and w denote measured- and unmeasured disturbances respectively, and t represents measurement noise. The symbol i simply denotes the outputs without measurements. Matlab requires that neither the manipulated variables u nor the measured disturbances d have an immediate effect on the outputs, in other words, the model is restricted t o be strictly causal. In terms of the above representation, this means that both D, and D must consist d of nothing but zeros. In PRIMACS, on the other hand, an Internal Model is utilized where at time k the inputs obtained at k - 1 are used, thus creating a strictly causal model. Appendix A.4 takes a closer look at this. Next, PRIMACS offers the possibility to specify output constraints as soft or hard. This is done using priorities which can vary from 1t o 100. Priority 100 indicates the output constraint is hard, meaning that under no circumstance is the hard constrained output allowed to exceed its limits. All other priorities indicate the relative importance'" of soft constrained outputs not t o exceed their bounds. The constraint handling strategy in PRIMACS is now such that whenever outputs exceed their limits, soft constraints are released in order of increasing priority. When, after this, the QP problem remains infeasible, the whole control problem is infeasible, since any hard constraints that might be present, are never released. This directly constitutes the problem that may arise when too many (or worse, all) output constraints are specified as hard, often rendering the QP problem infeasible leading t o a complete controller stop. In fact, Matlab doesn't use any type other than hard constraints, with all consequences of which. With respect t o releasing constraints, it should be mentioned that once certain constraints are released, the algorithm doesn't 'care' how big the violation is. The QP problem could not be solved with this particular constraint, so the controller acts as if the constraint is nonexistent until the QP problem with constraint is feasible again. Therefore, it could be useful to penalize the violation in some suitably weighted criterion, for instance as an extra quadratic term in the optimization. This way, while the constraint is being released, the controller will be forced, for example, t o under-perform somewhat in the setpoint tracking in order to reduce constraint violations. The question remains whether this can be realized without the algorithm becoming overly complex and, consequently, computationally demanding.
'This would involve programming new algorithms, which in fact is known to have been done. "The actual values of the priorities are irrelevant, as long as they reflect the order in which constraints can be released. So, priorities serve as 'tags' rather than weights.
2.5. M P C IMPLEMENTATIONS: M A T L A B A N D PRIMACS
29
If repeated linearization is t o be used, as will be done in Chapter 5, it must be possible t o update the Internal Model every number of samples. This can be done in PRIMACS, but not (to our knowledge) in Matlab. More specifics about the repeated linearization will be given in Chapter 5. Finally, a minute feature of PRIMACS is the possibility t o disable anticipation. The advantage of MPC t o anticipate for future setpoint changes as well as known future disturbances, can now be disabled t o imitate the effect that nothing is known about future changes. Actually, this possibility was added, more than anything else, to visualize the value of anticipation in comparison with no anticipation (which is characteristic for most other control strategies).
Chapter 3
The gas turbine installation

In this chapter a general description of the gas turbine installation will be given. Construction and operation of the complete gas turbine and the several elements it comprises will be discussed. Modeling of the gas turbine will be treated in Chapter 4.
3.1
Construction and operation
A gas turbine produces mechanical (shaft) power by expansion of compressed gas-usually

air-through a turbine. To produce an expansion in the turbine, a pressure ratio must be provided which explains the need for a compressor. To develop useful power, energy has to be added to raise the temperature of the gas prior to expansion. This is accomplished in a combustion chamber which is positioned between compressor and turbine. Basically, this describes the gas turbine in its simplest form. Figure 3.1 schematically depicts this setup. Furthermore, Figure 3.2 shows the p - V and T - s diagrams of the (idealized) gas turbine
Fuel
Combustion chamber
Compressor
Turbine
Figure 3.1: The gas turbine in its simplest form cycle.
30
3.2. THE COMPRESSOR
31
I
P
I
T
V-
s-
Figure 3.2: p
-V
and T - s diagrams of the gas turbine cycle
In practice, part of the energy added is used to overcome losses which occur both in compressor and turbine. These days, gas turbine engines are getting more and more competitive thanks to increased component efficiencies, to name one factor. The second factor determining the performance of gas turbines is the turbine working temperature. The higher this temperature is, the better the turbine efficiency. Naturally, this temperature can not rise unlimitedly, since at some point material-dependent physical bounds are exceeded. Development of ever better turbine materials have increased critical stress levels, thus allowing higher gas temperatures. In the following, we will describe the gas turbine component-wise, starting with the compressor.
3.2
The compressor
Compressors can be divided into centrifugal compressors on the one hand and axial compressors on the other, both having their specific advantages'. Their concepts of operation, though, are basically the same. Therefore, we will confine ourselves t o describing the centrifugal compressor only, partly because the laboratory installation uses this type of compressor. The centrifugal compressor essentially consists of a stationary casing containing a rotating impeller which accelerates the air, and several fixed diverging passages in which the air is decelerated, causing the static pressure to rise. This accounts for part of the pressure rise
'For instance, axial compressors can deliver both higher pressure ratios and larger flow rates for a given frontal area, which can be of great advantage, especially to jet engines. Apart from that, they operate at potentially higher efficiencies than their counterparts. On the other hand, centrifugal compressors are characterized by shorter lengths than equivalent axial compressors, better resistance to foreign object damage, less susceptibility to loss of performance by build-up of deposits on the blade surfaces and the ability to operate over a wide range of mass flows at a particular rotational speed.
32
C H A P T E R 3. THE G A S TURBINE INSTALLATION
only; the remainder of the pressure rise is produced in the impeller itself, by centrifugal effects. Figure 3.3 is a diagrammatic sketch of the centrifugal compressor.
lesvina mt0r
Resultant V?ioatY
Outward velocity
Tangential
Volute
-Outlet
Figure 3.3: The basic construction of a centrifugal compressor. For explanation of the terminology, the reader is referred t o [Cohen et al. 871
More of interest to us are the inherent restrictions which compressors display. To explain the compressors operation, we consider Figure 3.4. The characteristic depicted here, can be described by the standard cubic characteristic representation
where +c denotes the compressor pressure rise, and 4 the (mass)flow through the compressor. W and H are parameters determining the shape of the characteristic, as is depicted. Departing from point A, if we move to the left along the Characteristic, it is seen that the compressor can deliver a higher pressure ratio only at the expense of mass-flow. When moving past the top of the characteristic, we meet one of the compressors restrictions. Being the most severe restriction of all, surge denotes the phenomenon associated with a sudden drop in delivery pressure, and with violent dynamic pulsation which is transmitted throughout the whole machine.
+co,
Surge is likely2 to happen when the mass-flow through the compressor decreases enough to enter the part of the compressor characteristic where slopes are positive. A positive slope
2The transition from stable to unstable operation may not happen immediately after the surge point (extremum of characteristic) is trespassed, because the pressure downstream may drop at a higher rate than the delivery pressure, to name one factor. So, the compressor might still operate at a mass-flow @ in Figure 3.4, without entering surge.
3.2. THE COMPRESSOR
33
?r
%
-W
"i t -
I
3
-@-
Figure 3.4: The compressor characteristic in standard cubic representation means that a decrease in mass-flow will be accompanied by a fall of delivery pressure. If the pressure downstream does not fall fast enough, the air tends to reverse its direction and flow back in the direction of the resulting pressure gradient. When this happens, the pressure ratio over the compressor drops dramatically, only to rise again when the pressure downstream of the compressor has fallen also. Unless the compressor is somehow forced to operate in the stable zone again (i.e. where slopes are negative), the surge instability is bound to persist, repeating the cycle of events at high frequency. Eventually, this can get highly destructive, which is why it is to be prevented at all times. The most effective way to get back into the stable zone, is to immediately open a recycle- or blow-off valve (see Figures 3.5 and 3.6). Without a turbine present, this simply causes the mass-flow through the compressor to
Compressor
:r
Blow-off
Figure 3.5: Compressor and plenum instantly rise, which could also be accomplished by opening the throttle-valve. With a turbine present, as is the case in the gas turbine installation, opening the blow-off valve results in a
34
C H A P T E R 3. T H E G A S TURBINE INSTALLATION
Figure 6: Preventing surge by opening the blow-off valve. Here mcomp- ti,, denotes the dimensionless mass-flow through the compressor and Delta- P- tilde the dimensionless pressure rise delivered by it
rapid pressure drop, again causing the compressor t o enter safety zone (this will be seen in Subsection 3.6.1). Another important cause of instability and poor performance is rotating staZZ. Although it can exist in the nominal operating range, it is known t o potentially contribute t o surge. When there is any non-uniformity in the flow or geometry of the channels between vanes or blades, breakdown in the flow in one channel, say B in Figure 3.7, causes the air to be deflected in such a way that C receives fluid at a reduced angle of incidence and channel A at an increased incidence. Channel A then stalls, resulting in a reduction of incidence t o channel B enabling the flow in that channel to recover. Thus the stall passes from channel to channel: at the impeller eye it would rotate in a direction opposite t o the direction of rotation of the impeller. Although ways exist to model rotating stall (see e.g. [Greitzer 761 and references therein), this will not be pursued in the present report. The final limitation to the operating range is known as choke. When moving to the right in the compressor characteristic, mass-flow rises simultaneously with pressure drop. In short, the compressor moves out of the working area it was designed for, thus causing losses which grow considerably with increasing velocities. At some point the position is reached where no further increase in mass-flow can be obtained, and choking is said to have occurred. Other factors restricting the compressor in its operation are imposed by other components in the installation and their limitations, the next of which to be treated is the combustion chamber.
3.3. THE COMBUSTION CHAMBER
35
Figure 3.7: The concept of rotating stall But before we do so, a little more can be said about the compressor characteristics we will actually use. Throughout this project, we first used a Garrett and then a BBC compressor. The Garrett compressor characteristic is essentially similar to that of the BBC, and is depicted in Figure D.l, which can (rather inconveniently) be found in the back of this report. As can be seen, several lines are drawn which correspond to different rotational speeds. The extrema of the various characteristics are connected by the surge line. Also depicted are lines of constant compressor efficiency (which in fact we will not use in our model-we will assume the efficiency to be constant and equal to 70 % , in case of the Garrett).
3.3
The combustion chamber
In the laboratory installation, a choice has been made for a tubular combustion chamber, with natural gas for fuel. For the gas to be injected into the combustion chamber, it must be under sufficient pressure. For this purpose, an mxiliary (constant speed, electrical) compressor is installed, with (intercooled) by-pass control and overflow valve. The tubular combustion chamber, which is fairly standard for conventional gas turbine configurations, basically consists of two concentrically mounted cylinders. These are the liner and the casing, as depicted in Figure 3.8. The combustion process takes place in, consecutively, the primary zone, the secondary or intermediate zone and the tertiary or dilution zone. In the primary zone air and fuel are mixed, by swirling mostly. This is achieved by flat vanes in between concentrical rings (also known as the air swirler) through which the air flow passes. The secondary zone mainly serves as a place for the flow to reside and be cooled to an intermediate temperature by means of the secondary liner holes, thus allowing the combustion to proceed t o completion3. Finally, in the tertiary zone the remaining air entering through the tertiary liner holes, mixes with the combustion gases and further reduces gas temperature, which is necessary in order to reach a turbine inlet temperature that doesnt exceed the maximally allowed. This basically
3At extreme temperatures (ca 2200 [KI), dissociation of the main combustion products prevents the combustion from reaching completion. However, at somewhat reduced temperatures, say 1800 [KI, combustion can pretty much proceed as normal.
36
C H A P T E R 3. THE G A S TURBINE INSTALLATION
Figure 3.8: The tubular combustion chamber describes the operation of the tubular combustion chamber. For detailed design specifications, the reader is referred t o [Van Essen 95al. It is important t o realize that fuel flow cannot increase arbitrarily quickly, since lack of air will cause excessive fuel to leave the combustion chamber unburnt . Alternatively, opening the throttle-valve all too quickly will result in very lean burn and possibly even flame-out. Furthermore, as was mentioned earlier, the maximum fuel/air ratio that may be used is governed by the working temperature of the highly stressed turbine blades. Moreover, opening the fuel valve rather quickly results in high temperature peaks, which may exceed the maximum turbine inlet temperature. This is due to the inertia which characterizes the gas turbine installation. Opening the fuel valve will result in almost immediate (turbine inlet) temperature rise, followed by an acceleration of the turbine axle which in its turn drives the compressor. This causes the mass-flow to increase, finally resulting in a temperature drop that renders the turbine inlet temperature to have increased only slightly compared t o its starting point. This effect will be seen in Figure 4.7 (where actually, the fuel valve is closed, with inverse effects).
3.4
The turbine
Although much can be said about turbines, we will restrict ourselves t o a few remarks only, since basically turbines are the inverse of compressors: expanders. While compressors need shaft power t o deliver a pressure ratio, turbines produce shaft power through expanding a pressurized mass-flow. Like compressors, turbines either are of the axial or radial type. The vast majority of gas turbines employs the axial flow turbine, usually having the highest efficiency of the two. However, when mounted back-to-back with a centrifugal compressor, the radial turbine offers the benefit of a very short and rigid rotor, which can be seen, for instance, in the Garrett configuration. Contrary to the compressor characteristic, for the turbine no dependency on rotational speed is present, see Figure 3.9. This is due to the fact that in a turbine the flow is effectively forced
3.5. T H E O V E R A L L SETUP
37
Figure 3.9: The experimentally obtained Garrett turbine characteristic t o converge, whereas in a compressor it has to diverge. Figure 3.10 shows this mechanism, where it can be seen that in the compressor the flow tends t o break-away while in the turbine flow separation will never occur. Instead, the flow at the smallest area will always be near to, or in a choking situation. The result of this is that, contrary to the compressor situation, losses in the turbine are very much independent of rotational speed. In principle, the outgoing flow of the turbine can be recycled) in which case it must be precooled before re-entering the compressor. In the system under consideration) though, the combustion products are exhausted into the atmosphere, and heat is not re-used4. Furthermore, we do not use this installation for producing external shaft power.
35 .
The overall setup
*In a cogeneration plant, the gas turbine drives a generator and the exhaust gases are used as a source of low-grade heat. Heat at a relatively low temperature is required for heating buildings and operating airconditioning systems, for instance. Other conceivable applications are paper drying and chemical processes.
38
C H A P T E R 3. THE GAS TURBINE INSTALLATION
Figure 3.10: The conceptual difference between compressor (a) and turbine (b)
Figure 3.11: The laboratory gas turbine installation diagrammatically depicted Figure 3.11 schematically depicts the complete gas turbine as employed in the laboratory installation. In this schematic view the three elements mentioned earlier can be recognized, viz. compressor, combustion chamber and turbine. Furthermore a buifer-tank (or plenum) is positioned in the section between compressor and combustion chamber. Its main function is t o decouple the flows out of the compressor and into the turbine. Other components worth mentioning are the blow-off valve, the return valve, the compressed air valve, the throttle valve and the fuel valve. The blow-off valve mainly serves as an antisurge valve, as was explained earlier. Apart from this, the blow-off valve can also be used to reach certain points in the compressor characteristic (setpoints) even more quickly than already possible with the other control valves. More about control goals will be seen in Section 3.6. The return-valve (between blow-off and buffer tank) and the compressed air valve both are incorporated for start-up purposes. A start-up facility must be present to provide sufficient rotational speed of the turbine axle for the compressor to compress enough air to drive the turbine. This can be accomplished either by connecting an external (electrical) motor to the turbine axle-which is elegant but very expensive-or by injecting (heated) compressed air
3.6. OPERATION A N D CONTROL
39
from an external compressor. This last option has been selected for various practical reasons (see [Van Essen 95al). During start-up the return-valve prevents instigation of the compressor. The function of the throttle valve (between buffer tank and combustion chamber) is t o control air flow through the gas turbine. Apart from delivering some desired mass-flow, it can also serve as a means t o quickly decrease turbine inlet temperature, by opening the throttle valve. Finally, the fuel valve controls the fuel flow into the conbustion chamber and thereby the power supplied t o the gas turbine.
3.6
Operation and control
We conclude this chapter with a discussion of control related issues. We will not yet go into model specific details which will be done in the following chapter.
3.6.1
Inputs and outputs
We will start with specifying in- and outputs for the gas turbine installation. For inputs we have the blow-off valve, the throttle valve and the fuel valve. More inputs than this seem unnecessary for our control purposes, since we don't seek t o control start-up and shut-down of the gas turbine. These procedures should still be carried out by hand. Any inputs less than this is hardly preferable, which we will explain by describing the valves' functions. An important valve for mass-flow control mostly, is the throttle valve ( S T ) ,which can be seen from Figure 3.12, where we used our (Garrett) model t o simulate the response t o closing the throttle valve from 30 % to 20 %. Apart from decreasing the mass-flow through the the compressor meomp, pressure pcomp drops somewhat as well. The same thing is seen in Most interesting is the rotational speed N and mass-flow through the throttle vaive mthr. turbine inlet temperature Ttbin which rises due t o smaller mass-flows. , The fuel valve mainly serves to control compressor pressure ratio, although not without influencing mass-flow through the gas turbine. This can be seen in Figure 3.13, where we opened the fuel valve (SV) from 23.7 % to 30.0 %. The effect described in Section 3.3 can already be seen here: when opening SV, the turbine inlet temperature quickly rises from 760 [KI t o 830 [KI, after which it drops again t o 790 [KI, which is only 30 [KI higher than the starting point. Still, trajectories exist for which the turbine inlet temperature at beginning and end do not differ at all. This is because of the fact that, when we vary SV only, turbine inlet temperature reaches a minimum value somewhere in the center of the compressor characteristic. Surge avoidance is the primary function of the blow-off valve. Apart from this, the blow-off valve can also be used t o reach certain setpoints faster, but only if these setpoints lie in a direction which can be reached by opening the valve, since this valve should preferably be entirely closed in a stationary situation. The response t o opening the blow-off valve is depicted
40
C H A P T E R 3. THE GAS TURBINE INSTALLATION

1
lo5. PcomP 9
~ 0
.1
3 mcomp r 6
1;.
io4
0.34
1.85
0.32
1.% 150
5OTtbiJO0
0.3 O
5a mthr
loa
150
6.6
50
100
150
-3 ;d
0
Q
"
'
'
'
'
E2 O o Q 1
O
.. ..
.: .
...
.:. . . . . .i.
. ..
. . . . .:. . . . .: .
...
'..
.. .
0.2 0.4 0.6 compressor flow [kg/s]
Figure 3.12: The response to closing the throttle valve ST from 30 % t o 20 % ( S B = 0.00; SV = 0.237) in Figure 3.14. It is clear that the blow-off valve is highly effective in reducing the pressure delivered by the compressor, since we opened S B by 5 % only. Other than that, turbine inlet temperature increases considerably, which is why using this valve must be done with caution. All in all, we conclude that these three inputs are the minimum number we should have if we want t o reach more or less every point within the compressor characteristic. Still, some points just cannot be reached, simply because certain combinations of (low) mass-flow and (high) rotational speed cause temperatures t o reach excessively high values. The outputs we (may) wish t o control are formed by the three variables which determine the position in the compressor characteristic. These are the pressure ratio over the compressor, the mass-flow through the compressor and the rotational speed. Of course, we will never seek to control all these three variables at the same time, because this would result in a n inherently over-determined problem: two out of these three variables fully determine the position in the compressor charact eristic. Apart from this, we may wish to define some variables as outputs in order t o specify constraints in MPC. These will be seen in the next subsection.
3.6.2
Constraints
We already mentioned the most important reason for choosing MPC as a control strategy, which was the possibility of explicitly incorporating constraints on in- and outputs in the optimization algorithm.
41
~ 2
lo5 I PcomP
mcomp
7~~~
1.8 O
5OTtbin00
150
0.35
100 50 mthr
150
50
100
150
0
Q
;a3
E2 o
750A
50
100
150
0.3 O
50
100
150
1 O 0.2 0.4 0.6 compressor flow [kg/s]
Figure 3.13: The response t o opening the fuel valve SV from 23.7 % t o 30 % ( S B = 0.00; ST = 0.30) First of all, valves should be constrained not to operate beyond fully opened and fully closed (saturations). This can be easily realized in MPC, whereas other control strategies do not offer a solution as straightforward as this. Just as importantly, maximum moving rates should be specified, since real-life actuators cannot simply move with arbitrarily high speed. In MPC this is realized by specifying the maximum move size per sample period. For instance, the throttle valve ST requires 120 [SI to move from completely opened t o completely closed. When using a sample time of 1 [SI, this means that the move per sample may not exceed &. Furthermore, some outputs may not exceed certain (physical) limits. In our case, the turbine inlet temperature is not allowed t o (continuously) exceed the bound of 950 [KI. Beyond this temperature, the stress exerted on the turbine blades becomes unacceptably large. Actually, this temperature is measured at some distance before the turbine, which is why in our model (where we did take the turbine inlet temperature t o be measured at the actual turbine inlet) this temperature is always somewhat lower than that actually measured. Therefore, we must lower the constraint proportionally. The actual bounds differ for the two configurations used (Garrett and BBC) and will be stated in Appendix B. The second output constraint stems from the need to avoid surge. For this purpose we will define a dummy output DEVSurge which measures the distance from the surge-line (in terms of mass-flow). Keeping this distance greater than zero then is the same as staying out of the zone left of the surge-line, which is exactly what we want. In principle, such a dummy variable could also be used t o detect choke. In practice, choke is not really a problem, not being an instability like surge. Still, we may want t o confine the area of operation by defining something of a choke-line, which the model is not allowed t o pass, since the real process will not either. Similarly, minimum and maximum rotational speeds may be specified, although especially the maximum rotational speed will probably never be reached before the turbine
42
~ i 2
;3 :q
C H A P T E R 3. T H E GAS TURBINE INSTALLATION
mcomp
0.32 i00 50Ttbin
1.8
1.6
1.4
150
750 800
~~~~
0.3 O
50 m b l O O
150
~
I
50
i00
150
50
100
150
O O
50
100
150
O 0.2 0.4 0.6 compressor flow [kg/s]
Figure 3.14: The response to opening the blow-off valve S B from O % t o 5 % ( S T = 0.30; SV = 0.237) inlet temperature passes its maximum.
3.6.3
Control objectives
In a practical gas turbine application, objectives one can think of are operating at constant rotational speed (for generator purposes), or, if the gas turbine is t o be used as a power supply for some external load (for instance in a jet engine) it might be important t o deliver constant power. Such demands are generally referred t o as peak-shaving, which basically means keeping either one of these quantities constant under various circumstances and when exposed t o all kinds of disturbances. Furthermore, operating at an optimal point can be of significant importance. This means operating at maximum efficiency with minimal by-pass (blow-off). Such an objective can be very hard t o implement, since optimal efficiency cannot easily be incorporated in the optimization. One could think of including compressor or turbine efficiency (or even better: the product of these) as an output with a constant setpoint at maximum efficiency. Likewise, one could define the blow-off valve position as an output with a setpoint at fully closed position in order t o minimize by-pass. The danger in such strategies is that if weightings are too stiff, some setpoints can never be reached. However, in this project, we will not primarily seek t o meet these real-life control objectives. Instead we will mostly pay attention to what we will call setpoint and it trajectory control. This means that we will specify objectives that are interesting from the control-oriented point of view, mostly.
43
For instance, we could try t o reach different operating points while keeping constant one of the following quantities:
-
mass-flow through the compressor pressure delivered by the compressor rotational speed
Furthermore, we might specify trajectories in which anticipation abilities are shown, or alternatively, interaction or nonlinearities are investigated.
In Chapter 5 we will mostly pay attention to trajectories which in way or another challenge the controller. That is, we will try to reach setpoints where one or more constraints are bound t o be exceeded. If the controller can keep these constraints from being exceeded, probably through ignoring the desired setpoint, it is said to perform acceptably.
Chapter 4
Modeling the gas turbine

Most of the system equations that are presented here, will be stated without explanation. For information on how t o obtain these equations, the reader is referred to [Van Essen 961.
4.1
The compressor model
First a model of the compressor and plenum' will be obtained. This model is similar to the one proposed by Greitzer in [Greitzer 761, which is a lumped parameter model where the compressor itself is modeled as an actuator disc. Figure 4.1 depicts this setup where the following definitions were used mass-flows through compressor, throttle and blow-off, respectively [kg/s] steady-state pressure rise according t o compressor characteristic [Pa] pressure drop at throttle valve and blow-off valve, respectively [Pa] volume [m3] of- and temperature [KI and pressure [Pa] at plenum temperature [KI and pressure [Pa] at inlet rotational speed (kept constant) [rev/s] pressure at outlet [Pa] (4.1) compressor efficiency [-I duct areas at Compressor, throttle and blow-off [m2] duct lengths at compressor, throttle and blow-off [m] cp/cV [-I, the gas constant [N m/kg KI and the density of air (at atmospheric conditions) [kg/m3]
'In general, a plenum is a chamber in which air (gas) pressure is greater than that in the outside atmosphere.
It can be a (fictitious) representation of 'volumes' from al over the system (i.e. a lumped p a r a m e t e r ) , or it can l be a real physical chamber such as a buffer tank. In our models, it will usually be a combination of both.
44
4.1. T E E COMPRESSOR MODEL
45
Figure 4.1: The basic compressor model This system can be described by the following equations
Pin
= Pout
where
In many cases it can be helpful t o nondimensionalize these equations'. We nondimensionalize 1 by multiplying: mass-flows X"A pressure differences x spu2 time where
WH
XWH
=a
d.PfEx,is the Helmholtz resonance frequency and U the tipspeed.

~
This yields the following equations
'Mostly for ease of comparison; different compressors with varying dimensions can be characterized in terms of a few dimensionless parameters.
46
CHAPTER 4. MODELING THE GAS TURBINE
The dimensionless parameters B and influence the dynamic behaviour of the compressor system. The first of these is by far the most important parameter in determining the system dynamics. Apparently, the bigger B , the more likely the compressor system is t o enter surge once the surge line is passed. This corresponds to running the compressor at relatively high speed or with a large volume behind it or some combination of the two. The second parameter3 only has minor effect on compressor operation and therefore usually is kept constant in parameter analysis types of experiments. qthr and qbl effectively serve as measuring-rods for the valve positions. Figure 4.2 shows the response of this simple dimensionless model to changes in q t h r , which is basically the same as the response to changes in qbl, apart from the fact that the mass-flow out of the compressor is divided differently over blow-off and throttle duct. Furthermore,
mcomp-tilde
[-I
Figure 4.2: The response of the dimensionless model to an increase in Figure 4.3 depicts surge cycles for B = 2.
qth,.
4.2
Extensions to the basic compressor model
Before moving on to modeling the complete gas turbine installation, some extensions to the thus far obtained most simple form of compressor model should be mentioned. When rotating stall needs to be described, a quasi-steady approximation of the compressor response will not be adequate, since there is a definite time lag between the onset of the instability and the establishment of the fully developed rotating stall pattern. In this case, a simple,
31n literature, is usually denoted by just G, and equal to
~~~~~~t~~
since no blow-off is present there,
see e.g. [Greitzer 761.
4.2. EXTENSIONS T O T H E BASIC COMPRESSOR MODEL
47
i.....
.........................................................
.............................................................. ......................
0 3 ...................... .
0.2 ..........................
01
.I
....................................................................
O
:.... ..
...........\i
..........
i
-0.5
rncornp-tilde
[-I
05 .
Figure 4.3: Surge simulated with the dimensionless model first order transient response is adopted to simulate this time lag in compressor response. In dimensionless form this can be written as
Note that, for a given compressor (R, N and Leomp constant), 7 is proportional to are
$.
Furthermore, the model can be augmented with differential equations for rotational speed N and plenum temperature Tpl.Especially the first extension seems very useful since our aim is t o model a gas turbine, which is hardly characterized by a constant rotational speed. Only the use of a constant speed (electrical) motor would result in such behaviour. Including rotational speed yields the following (dimension-full) equations dN IN= Pi, - Wp dt with
where I is inertia, Pi, the in-going power and W p the power requested by the compressor. The incorporation of dynamic temperatures improves realism, but can be argued, since this type of low-frequency dynamics could also be chosen to be ignored, just as slow varying ambient conditions will not be accounted for. In dimension-full form the corresponding differential equations take the following shape
48
C H A P T E R 4. MODELING THE GAS TURBINE
where Ts, = Ti,

css
(2)
1 3.577comp
)
with Pss the steady-state compressor (absolute) pressure (=
+Pin).
Finally, one can choose to replace the differential equations for the mass-flows by static relations) simply because the dynamics involved here are generally so fast that mass flows change instantly with other system variables. By doing so, computing times for the model can be somewhat reduced. On the other hand, defining mass-flows as outputs instead of states can cause major practical problems when trying to implement the model in the Matlab MPC Toolbox, i.e. MPC-Tools does not accept any models other than strictly causal, so outputs cannot be explicit functions of inputs (as was discussed in Section 2.5 and Appendix A.4).
4.3
Modeling the complete gas turbine
To obtain a model that combines compressor) plenum, combustion chamber and turbine) first some things need to be said about the latter two components.
The combustion chamber can be modeled in various ways. We wiU present three methods, based on the following scheme:
1. N o plenum for the combustion chamber: in this case the temperature at the combustion chamber Te,will have to be determined statically.
2. A plenum to represent the combustion chamber:
(a) Teecan still be determined statically;
(b) Teecan also be determined dynamically. Other variations can be imagined, depending on how much detail is needed for in the description. For instance, temperatures can be obtained in a physical way, but just as well through a 'black-box' function of the fiJel-valve position SV.
ad i: The simplest way is to represent the combustion chamber as a temperature rise only, taking place somewhere between buffer tank and turbine, in which case the turbine inlet pressure (ptbin) is related to the pressure in the buffer tank (or plenum: p P l )through a static equation. This would take the following form:
TCc -p l = T
CPT
Cpcomp
Power
mthrcpT
4.3. MODELING T H E COMPLETE GAS TURBINE

= P p l - AP
49
&bin
where C p T denotes the specific heat corresponding t o the (heated) air/uel mixture (T indicates temperature) and Ap the pressure drop caused by pipe lines and throttle-valve, mostly. T t b i n and Tpl are temperatures at turbine inlet and buffer tank. Power serves as the input used t o obtain some desired Tec. so desired, Power can be a function of If the fuel-valve position SV. The temperature rise can be obtained in a more detailed way, by actually determining the fuel flow into the combustion chamber from the fuel valve position and the pressure ratio over this valve. Next , the temperature is computed using specific heat capacities and combustion energy. This method-which assumes pressure and temperature of the natural gas t o be fixed-yields the following equations
Tcc
mfuel3.85 l o 7 f mfuelCpfuel313
CpT(ljZfue1
mthrCpcompTpl
%thr)
Pccbar
Ptbin
105
(this is p t b i n in [bar])
where mfuel the fuel flow, e p f u e the specific heat for fuel (usually natural gas), and is l ICvfuer (momentary4) fuel valve capacity [m3/hr]. The factor 3.85. lo7 is the heat of the combustion for natural gas.
o
ad 2a: If it is desirable t o decouple p t b i n from pp1, an extra plenum at the combustion chamber is required5. With this extra plenum representing the combustion chamber, the resulting equations (for the combustion chamber) will look something like this:
is where Vee the (combustion chamber) plenum volume, m t h r the mass-flow through the throttle valve and l j Z t b the mass-flow through the turbine.
e
ad 2b: In our final model, the combustion chamber is represented, not by its temperature, but by a Power term added t o the differential equations of &bin and T t b i n , which corresponds t o the power actually being supplied at the combustion chamber plenum, instead of somewhere between buffer tank and combustion chamber (in the form of a
* K W f u e i determined from some valve characteristic which describes the relationship between hvfuei is and the fuel valve position, SV. This characteristic can be proportional or, just as common, exponential (Kvfuel Kvfueimaz = . where SV = 1.0 corresponds to a fully opened valve). More about valve characteristics can be found in Appendix B. 5This also causes these variables to become states instead of outputs, which can be especially useful if too many variables (outputs) are algebraically interconnected, which causes (symbolic) linearization to get overly complex, and, eventually, the nonlinear model will have to be linearized in order to produce a model which (linear) MPC can use. We will pay more attention to this in Section 4.4.
50

temperature rise). This power can either be determined in a physical way, or alternatively, simply as a fitted (black box) function of the fuel valve position SV. The latter option will be used in our model.
To model the turbine, we use its characteristic, which was depicted in Figure 3.9. As was
mentioned before, we can either determine the mass-flow dynamically or statically, i.e.
or
mtb
= fturbinetptbini Ptbout, T t b i n )
where Psstbis the steady state turbine inlet pressure needed for a certain turbine mass-ow, and fzurbine the inverse of Psstb6. We will now state our model of the complete gas turbine installation, where we have chosen to use static equations for all mass-flows. Furthermore, yet another plenum-this time between compressor and buffer tank-is added, t o reflect the fact that in reality the blow-off is positioned near the compressor and not at the buffer tank. Apart from this, it enables an overall pressure drop t o be modeled apart from the throttle valve pressure drop. This setup is depicted in Figure 4.4.
x
Ptbin Ttbin mthr
4
qt b
Power( SV)
N ,I
Figure 4.4: The final gas turbine model setup
For the actual fits, the reader is referred t o Appendix B.
4.3. MODELING THE COMPLETE GAS TURBINE
51
DEVSurge =
mcomp
N) 2AComp (N)Scale
Bcomp (
Scale is defined in Appendix B, see the explanation below, also.
S B ST
SVJ
52
The various fits (f...) in this model can be found in Appendix B, both for the Garrett- and the BBC-configuration. The static equation for m,, is simply obtained from the pressure drop over a duct with area A, inlet pressure and -temperature pl and T I ,outlet pressure p,, flow speed u and friction factor [, which is given by
where determines the size of the overall pressure drop, which-in and located between plenum Vcomp the buffer tank (Vpt).
our model-is
(fictitiously)
Furthermore, the additional variable DEVSurge is introduced t o measure the distance from the surge-line, i.e. the difference between the actual mass-flow through the compressor and the mass-flow for which the surge-line is crossed (at the same rotational speed). Since the compressor characteristic is fitted with second order polynomials (pressure ratio = Acomp(N)(~compScale)2 Bcomp ) ~ c o m p S c a l eC e o m p ( N ) , see Appendix B) this simply (N boils down t o
DEVSurge =
-
mcomp -
camp,
d pressure ratio-O
dmcomp
mcomp
-Beomp ( N ) 2AcO,,(N)Scale
where Scale mainly scales the mass-flow obtained from the characteristic t o the right dimensions (see Appendix B). DEVSurge can (and will) also be used t o (try and) prevent choke. Although in real-life choke doesnt constitute much of a problem, it is important with regard t o modeling the gas turbine. The fit is simply not equipped t o describe the actual process (far) beyond this line, where the mass-flow should-but in the fit doesnt really-automatically reach the point which it cannot pass (see Figure 4.5, the most of which will be explained in the following). The choke line is sometimes? considered t o be similar t o the surge line, only t o be shifted t o the right (towards higher mass-flows). In reality, this is certainly not exactly true, but it seems useful as a first approximation. However, it should be noted that in our case, even the just mentioned approximation is violated in some way. If we define the choke-line t o lie at, say, DEVSurge=0.2, this is not the same as shifting the surge-line over a distance of 0.2, because the choke-point is determined from the surge-point ( f 0 . 2 ) at the same rotational speed and, more importantly, for higher N , the curves in the characteristic fall off faster for increasing mass-flows. The effect of this is best shown in Figure 4.5, where our choke-line is plotted together with the earlier suggested approximation and the real choke-line. As can be seen, the choke-he is really riot very wel! described by either of the suggested approhirr,atior,s. However, this will not be of crucial importance.
7See e.g. [Satter 961.
4.3. MODELING THE COMPLETE GAS TURBINE
53
t
O
0.1
02 .
03 .
0.4 05 . compressor flow [ks/s]
0.6
07 .
08 .
Figure 4.5: Choke-lines in the (Garrett) compressor characteristic; -: - . - : suggested
real, --: ours9,
Also stated were 2,y and g-the states, outputs and inputs of the model, respectively. Note that in our final mo(ie1 that was presented here, some of the outputs, for instance m,,, are explicit functions of the inputs, in this case S B . We mentioned earlier that this is not allowed in Matlabs MPC-Tools. The model presented here is the model used in PRIMACS, whereas for Matlab simulation models, dynamic mass-flows were used, but since these models can be obtained rather easily from what was discussed so far and since the PRIMACS results are much more valuable to us, the Matlab models are omitted here. In Figures 4.6 and 4.7 we compared our (Garrett) model to measurements of the real installation. Figure 4.6 shows stationary points in the compressor characteristic for different valve positions. The left line corresponds to ST = 0.3, while the right line depicts points for which ST = 1.0. As can be seen the model quite well describes stationary behavior, both for several fuel valve positions and for different throttle valve positions. In Figure 4.7, where the model responses are pointed to by arrows, we compared modeland actual response to a ramp-wise change in the fuel valve position from 35.3 % t o 8.4 %. Actually, in the simulation, SV was closed until 11.3 % only, since the fit of the compressor characteristic is no longer valid for stationary points below the point we stopped at. This will cause the model to reach a point somewhat higher up in the characteristic, which is indeed the case, judging from Figure 4.7. Most model responses quite well simulate the actual response, although the turbine inlet temperature drops somewhat further than it should have. Still, it is hard t o tell how big the error really is, since turbine inlet temperature was already said to reach a minimum somewhere in the center of the compressor characteristic, t o rise again when moving away in either direction along the line of constant ST and S B (only SV is varied). This means that indeed this temperature would be expected t o be somewhat smaller
54
C H A P T E R 4 . MODELING T H E G A S TURBINE
0.1
02 .
03 .
0.4 0.5 compressor flow [kgis]
06 .
0.7
08 .
Figure 4.6: Stationary points in the (Garrett) compressor characteristic; experimental (*) and simulated (o) for SV = 11.3 % than for SV = 8.4 %. Other than this, rotational speed seems t o be a little too slow in its response. This could not be counteracted by decreasing inertia I, since the time-constant displayed in the transient response of the (model) rotational speed is fully determined by the slower system dynamics, being the plenum response.
4.4
Linearization
In order t o obtain a linear model which (linear) MPC can use for its computations, we must linearize the nonlinear differential equations. With the states i, outputs y and the inputs the - as defined in the previous section, we can represent the (nonlinear) system-equations by u
where
dx 2= x.
Linearizing this system would result in the linear system

-
a = Ai!$)
+ &(t)
- = C i ( t ) Du(t) y
where A ( i , j ) or instance contains the partial derivative af (see, e.g. [Kok 50, Section 3.41).
2
%EJ
Substituting the point (states and inputs) around which is linearized (this is usually the initial condition) yields the (quantified) matrices A , B , C and D.
4.4. LINEARIZATION
55
"PI-!
Figure 4.7: Comparison of model and reality: transient response to change in S V from 35.3 % to 8.4 % /11.3 %, where (a): S V (lower 2 lines); (b): N (lower 2 lines); (c): h t h r ; (d): Tcomp; (e): Ttbin (upper 2 lines and 'arrow-line'); (f): peomp(upper 2 lines)
If the outputs contain mutually coupled (algebraic) expressions, linearization can get rather complex. For instance, if we have two outputs
yi = g,(-,U,y,)
-2
= g2(-,u,y,)
it might be impossible t o explicitly express both the outputs in terms of g and U, only. In that case it may be possible to obtain the partial derivatives by solving a set of algebraic equations. If, however, the number of mutually coupled outputs gets considerably larger, the linearization gets rather complex'. To avoid this, certain outputs can be defined as states either by using some physically justifiable differential equation or by introducing a first-order time lag approximation. For this very reason it was convenient t o decouple ptbin and pp by introducing a third plenum, which was described in the previous section.
As was seen in Subsection 2.4.1, the linearization can be obtained analytically or numerically. We chose the first option and used Maple to obtain the symbolic expressions for the partial derivatives.
This linear model next was evaluated at the initial conditions, the result of which is shown in Figure 4.8, where we compared linear and nonlinear model (Garrett) for a step in S V from 23.7 % to 28.0 %. For this (not too big a) step, the deviations between the two models seem rather small, maximally about 15 %. Of course, the linear model will be accurate only within the proximity of the original linearization point. Repeated linearization (evaluation of the linear matrices at multiple points) is expected to improve accuracy.
More about the implementation of the model in MPC will be seen in the following chapter.
'For q algebraically coupled outputs, this would involve the symbolic inversion of a matrix with dimensions
4 x Q.
56
C H A P T E R 4. MODELING T H E G A S TURBINE
4rnc
14 . 0 i
rnthr
rntb 4
0.38 0.35 0.36 0.34
0.38 0.36
50
PP1
100
150
0.3
* :::f-q F l
O
50ptbi,,100
150
0.34 O
50
100
150
1200
O
400
390
50 TDllOO
150
I
i .a O
a20i
50Ttbi,100
150
1
1150
50
100
150
380
n -7 if
Figure 4.8: Linear and nonlinear model compared; -: nonlinear, --: linear
Chapter 5
MPC applied on the gas turbine:

simulations
In this chapter we will state results of MPC applied on the gas turbine, obtained with PRIMACS. Although we did run numerous simulations in Matlab, we will pay no attention t o them, since these simulations-not incorporating a nonlinear model-seem t o be of little value compared t o those run in PRIMACS. We will discuss these results and draw some conclusions. Furthermore, we will design an experiment which is t o be carried out on the real laboratory installation. Out of safety considerations, this will be a reasonably mild trajectory. But before we do so, we must say some things about the implementation. Specifics can be found in Appendix C.
5.1
The MPC implementation
Most of what will be treated in this section, stems from handling numerical problems. After that, we will state and motivate the constraints that were implemented, and our final MPC parameter settings.
5.1.1
Scaling
The first of these issues is the need for scaling. Basically, this is a rather trivial problem, but nevertheless it is essential t o solving the control problem.
The results of this will be discussed in the next chapter.
57
58
C H A P T E R 5. MPC APPLIED ON THE GAS TURBINE: SIMULATIONS
Systems that are characterized by widely varying time-constants, yield stiff differential equations that cause optimizations (but even simulations) to become ilZ-conditioned2. Proper scaling can3 alleviate these problems.
The concept of scaling is quite simple. For example, lets consider mcomp its scaled version and - mcomp where cl = 0.5 ( w order of magnitude of hcomp). mcomp Likewise, we
CCLE
define scded versions of
states xi j % = outputs yi = inputs u =+- = i
2 5 2
With this we get
% = - c U ( i , j )= $A(i,j) 2 CC .
yielding
Similarly,
define scaled versions of B , C en D :

-
B(i,j) =
Eventually, E, y and g can be obtained from xi = ei . %, etc. This basically describes the scaling procedure. In PRIMACS the same procedure is used, but next t o that PRIMACS sometimes (internally) uses scaled signals as well. The definition of these scaled signals is slightly different from the earlier described scaling method. By subtracting an estimated mean of the signal, we get a signal with (approximately) zero mean. Next, this signal is scaled like before. This yields the following expression for e.g. x i
where
xi,offset and xi,amplito&
denote offset and amplitude of xi.
2Definitions of these terms can be found in various works on mathematics, such as [Borowski and Borwein 891. 3For very stiff systems, scaling may not be sufficient t o prevent numerical problems.
5.1. THE MPC ILMPLEMENTATION
59
5.1.2
Discretization
Another important issue is the discretization of the matrices A , B , C and D.We obtained the discrete matrices by using
cd
Dd
CeA(At-8)
[iAt-e +
eATdr]B
c d
where At is the sampling interval and the 'data procession' time (see [Kok 90, Section 6.11). In our case 6 = At, since computed moves (based on data from sample k ) are implemented at the following sample ( k i), so Cd and D take the following form d
c
D
D d
Ad and
B d
are often seen t o be referred t o as Q> and I', respectively.
We implemented this by including the Matlab function c2d (converted t o C++) in the PRIMACS source-code. This function computes the exponential functions by using either Taylor or Pad approximations, preceded by scaling of the matrices A and B and followed by repeated squaring t o undo the scaling (see [Matlab 921). The corresponding PRIMACS source-code can be found in Appendix C.
5.1.3
Reconstruction
If repeated linearization is to be used, all model states will have t o be known reasonably accurate (see Subsection 2.4.1). This means that states that aren't available as measurements will have t o be reconstructed one way or another. Section 2.3.5 proposed the use of a Kalman filter for this purpose, possibly even combined with a disturbance-model. However, in the end we did not apply this filter, because time did not permit t o get the implementation t o work properly. Instead, we 'reconstructed' unmeasured states in a somewhat 'ad hoc' way. With the following measurements available
Peomp Teomp
Ptbin
Ttbin
mthr
we can compute p p l and T,l, if we assume that Tplevolves exactly like Teomp (thus differing only in their initial conditions4), by using the pressure drop over the throttle valve (for subcritical flow, see Appendix B)
- 4 S 0 we assume that T,l = Tcomp.
60
CHAPTER 5. MPC APPLIED ON THE GAS TURBINE: SIMULATIONS
*
and
Ppl
= ?:2pt:r
(STKuT,,,)
1
mthr
7.0
4 -
Tpl from Tcomp using equation (5.1), but modified such that we get by
Tl = TcompTpl,amplitude p
Tpl,ojjset.
The exact corresponding source-code can be found in Appendix C.
5.1.4
Constraints
In the following simulations we specified the following constraints: DEVSurge low constraint DEVsurge high constraint Ttbin constraint high
: 0.01 [kg/s] : 0.2 [kg/s] : 950, 900 and 750 [KI consecutively
(5-2)
The three values for the maximum temperature constraint correspond t o different models or different safety margins. At first, when we were still using the Garrett model, we took this constraint to be 950 [KI. After a while, we felt that-out of safety precautions (since the constraint was slightly exceeded)-the constraint would best be lowered t o 900 [KI. When finally, we switched to the BBC model, we changed the constraint t o 750 [KI, because this model continuously resulted in a far lower Ttbin. Next, we used the following input constraints and input move constraints:
S B low constraint high constraint move constraint ST low constraint high constraint move constraint SV low constraint high constraint move constraint
1.0 [-I 1/72 [l/s] 0.01 [-I 1.0 [-I 1/120 [1/s] : 0.01 [-] : 1.0 : 1/48 [1/s]
: : : : :
: 0.0
[-]
(5.3)
[-I
Most of these constraints are rather obvious. We specified low constraints on ST and SV slightly larger than zero, since it seems rather useless (and maybe dangerous) to completely close these valves. The move constraints were taken from [Van Essen 95a]. For instance, S B needs 72 [s] to move from fully closed to fully opened. In the laboratory installation, these valves always move at maximum rate, but can of course be made t o move during part of a sampling interval only. Other than this, we did not specify any constraints, mainly t o keep computational efforts (which grow considerably with increasing numbers of constraints) within practical limits.
5.i. THE M P C IMPLEMENTATION

5.1.5
61
Parameter settings
In the simulations, we used the following MPC parameter settings (unless stated otherwise): sample-time (= MPC execution-time) discretization time-interval integration sample-time prediction horizon p control horizon m constraint window starting-sample filter constants
=
= = = = =
1.0 [SI 1.0 [SI 0.0001 [s] 15 [samples] 4 [samples] 5 [samples]
(5.4)
Teornp, P t b i n ,
= 0.1 for
iV;-
peornp,
Ttbin,
eomp7 m t h r and DEVSurge; all other filter constants are zero.
Consecutively, we will motivate all of these choices. In Subsection 2.3.1, we mentioned the fact that the (MPC) sampling period can not be chosen arbitrarily large, since it is also used as discretization sample time. This is why we chose the first two sampling intervals t o be equal to 1 [SI. From Matlab simulations, it was concluded that such a sample time still allowed for reasonably accurate modeling, while on the other hand it caused no considerable computing-time related problems5. In PRIMACS, this sampling period proved t o be rather large, although with the use of somewhat faster computing hardware, this problem is expected not to be insuperable. Other than that, the integration routine used in PRIMACS t o simulate the nonlinear model, might quite easily be improved. At the moment, this still is an Euler forward integration scheme. Of course, the execution time also largely depends on the tuning parameters that are t o follow. The integration sample time again was chosen as large as practically possible. Choosing a period larger than this, resulted in highly inaccurate (or even unstable) model response. The prediction horizon p was determined as a compromise between minimizing computational efforts and predicting slow system dynamics properly. Although the gas turbine displays maximum time-constants of about 30 [SI, we decided on choosing p = 15, in order to reduce execution time. Most of the other rules of thumb that were presented in Subsection 2.3.2 do not apply, since the gas turbine does not exhibit any dead times, and furthermore inverse responses do not really6 occur in the quantities we wish to control. The turbine inlet temperature does respond inversely t o changes in S V , but we only seek t o keep this output beneath its maximum. Next, we chose the control horizon m = 4, in accordance with the recommendation t o choose
5Actually, in Matlab, we did not check if MPCs execution time never exceeded the limit of 1 [s] per sample. We merely found that for a simulation run of 150 [SI, MPC needed about 20 [SI (depending on the number of constraints that were active). This might s t 2 mean that, say, for the first sample, MPC needed (say) 3 [s], which of course is unacceptable for real-time implementation. The mass-flow through the compressor is sometimes seen to display an inverse response typically lasting for about 1 [SI so this is of RO importance to us.
62
C H A P T E R 5. MPC APPLIED ON T H E GAS TURBINE: SIMULATIONS
it somewhere between 1/6 and 1 / 3 of the prediction horizon, which in our case corresponds t o 2 t o 5 samples. This control horizon was found t o result in acceptable controller behaviour. We decided t o use a constraint window starting from the 5th sample, because in several simulations we ran with the Garrett model, not using a constraint window resulted in early release of constraints, after which (due to the effects of, for instance, turbine inlet temperatures of almost 2000 [KI), the controller could not manage t o get back to the desired situation. In case of the BBC mode!, we did aot yet experience these problems, which is why simdation 9 (see Section 5 2 , where we didnt use a constraint window, does not show any of these .) problems. Finally, we defined filter constants # O for the outputs that are in fact measured. Furthermore, mcOmp DEVSurge are filtered as well, where we assumed that these outputs can and be reconstructed from peomp and N via the compressor characteristic and the definition of DEVSurge, respectively. In order not t o take filter settings too tight on the one hand, but still make use of process measurements on the other, we chose the filter constants equal t o 0.1. These parameter settings were used in all simulations, unless explicitly stated otherwise. Parameter variations will be stated and motivated in the next section. All other parameter settings can be found in Appendix D which should accompany the simulations treated in the next section.
5.2
Simulation results
This section contains results obtained with the Garrett model (the first two simulations), followed by those obtained with the BBC model (all other simulations). In all figures, we used the following convention:
5.2.1
Simulation 1
We start with an example of setpoint control, depicted in Figure 5.1 and referred t o as Simulation 1 (Garrett). Here we specified setpoints for mcOmp from 0.35232 t o 0.2 [kg/s] simultaneously with a compressor pressure (pcomp) drop from 187899 to 150000 [Pa].It is seen that the turbine inlet temperature reaches its constraint, but doesnt exceed it (not even the linear model does), which is exactly what we expected from MPC. Furthermore, if one compares the valve positions t o the temperature response, it becomes clear that at 10 t o 15 seconds before the temperature constraint is about t o be exceeded, MPC decides t o open the fuel valve less quickly in combination with closing the throttle valve
5.2. SIMULATION RESULTS
63
rncomp
100 200
1O0
200
1O0
200
Ttbin
DEVsurge
O a 2
pcornp vs. mcornp
100
200
1O 0
200
0.5
S5
ST
SV
0.2
0.2
0.15
0.15
100
200
1O0
200
o. 1
1O0
200
Figure 5.1: Simulation 1
slower, both actions of which can be easily seen (see e.g. Figures 3.12 and 3.13) to slow down temperature rise. Next to these two valves, in the beginning the blow-off valve is opened at maximum rate for about 5 seconds and then closed at a slower rate for another 15 seconds. This is especially useful to decrease pcomp quickly (see Figure 3.14), after which the blow-off must be closed again to prevent this pressure from becoming too low again. It can also be noted that DEVSurge is especially sensitive to differences between linear and nonlinear model. If DEVSurge is close to a constraint this may cause the controller to become overcautious (which isnt such a bad thing), since the linear model has the tendency of (always) responding more violently than the nonlinear model. In this simulation, we linearized on each sample, which is basically why this rather difficult setpoint change (since a constraint is reached) is performed reasonably well. Without repeated linearization, the linear model might have become too inaccurate to keep this control problem from getting infeasible.
64
5.2.2
C H A P T E R 5. MPC APPLIED ON T H E GAS TURBINE: SIMULATIONS
Simulation 2
This is definitely true in Simulation 2, where we tried to reach a point where the turbine inlet temperature would exceed its maximum7. This was done by specifying a setpoint change in ?-hcomp from 0.35232 t o 0.2 [kg/s], while maintaining a constant pcomp.Figure 5.2 shows that the gas turbine displays highly interactive behaviour, because when we try to decrease ?-hcomp, pcomp drops as well, whereas it was supposed to be kept constant.
mcomp
N 1200 7
1000.
i, .
100
200
-0
~ : ~ ~
DEVsurge
100
200
800 I o 100 200
!,--.-------
500'
100
200
100 200
ST
SV
100
200
100 200
- 0
100 200
Figure 5.2: Simulation 2 Figure 5.2 shows that in order to keep the temperature from violating its constraint, MPC decides to (more or less) ignore the setpoint on pcomps. At first, MPC still tries to get this pressure back t o its initial value (by opening S V ) , but pretty soon comes t o the conclusion that Tibin exceed its maximum, if no counteraction will
~~
71n this simulation (which actually preceeded Simulation i) we still used a high constraint of 950 [KI for the tempeiaiwe, wheeas in the sbsecperit Gaiett simdations this constrint was lowered to 900 [KI out of safety reasons. 8 0 f course, this 'choice' heavily depends on the weights used for hcomp pcomp. With different weightings, and MPC might just as well decide to let go of the other setpoint, or-for that matter-meet both setpoints halfway.
65
is undertaken. Immediately, ST is being opened again in conjunction with closing S V , which results in the linear temperature not exceeding the maximum. Unfortunately, the nonlinear temperature quite severely does, but only for a short period of time. Apart from all this, (linear) DEVSurge reaches its constraint as well, but is kept from violating it by the very actions that cause T b nto stay beneath its maximum. It is also seen that in ti this simulation, where we did not (yet) filter DEVSurge, linear and nonlinear version display a huge diEerence. This causes the c o n t d e r to Secorne wedy caiitim, were it m t f the m fact that control actions were undertaken anyway, since Ttbin nearing its constraint. was When looking at what actually happens, it is not surprising that MPC more or less releases the setpoint on peomp,since the number of setpoints (2) added to the number of active constraints (ior 2) exceeds the number of effectiveginputs (2). With effectively 1setpoint and 1 constraint active, the problem becomes feasible again. This especially hard control problem could not be 'solved' if we linearized only once, which is shown in Figure 5.3. Soon the model becomes too inaccurate to keep the problem from
mcomp
1.6
120
0 . . . . . .. . . . . . . . . . . . . . o,2 . 3
0.1
40 Ttbin
20
40
~izl
1000 .............
,
DEVsurge
- 20
'~'..."'.
40
4.1'
20
40
ST
0 . 4 7 1
0.41
sv
I
O ' O
20
40
Figure 5.3: Same simulation as Simulation 2, only now we linearized only once becoming infeasible. Still, we must mention the fact that in the repeated linearization we used in Simulation 2, the model was (incorrectly) not updated for new inputs-only for the outputs. The more surprising it is, that this faulty repeated linearization still yielded better results than the simii!ation where we Enearized only once.
' S B just cannot contribute in increasing pcomp.
66
C H A P T E R 5. MPC APPLIED ON THE GAS TURBINE: SIMULATIONS
5.2.3
Simulation 3
In this simulation (BBC) we tried to do something similar t o what was done in Simulation 1, but unfortunately stumbled upon some inconvenient effects. Here we wanted t o decrease pcomp from 164171 t o 135000 [Pa] together with rheomp decreasing from 0.4919 t o 0.33 [kgls].
The results (Figure 5.4) show that for octpxt weightings 5 fm pcompand 1 for meompi the desired pressure is reached, but not the desired mass-flow. Although this seemed unexplainable at first, eventually it was found to be caused by the rather inconvenient character of the throttle valve.
mcomp
,_
Ttbin
- - .- . . . . - - -
0.4
0.35O O
m
~ ~
.. 1O0 _ _ - a -4. .
~ ~
600
0-05 oF i i
O
1O0
1O0
SB
ST
SV
O O
1O0
0.95
1O0
1O0
This valve has such a large capacity, that when it is operated somewhere near fully opened, it is almost without effect. The pressure drop over this valve is just too small to be of any contribution to the process. Only when the valve is closed beyond 50 % the effect becomes noticeable. Since we chose an initial condition where ST = 1.0 (fully opened), this behaviour is something we will have to live with, unless of course we choose to take a different initial condition. This is something we recommend for future simulations, but didn't do ourselves, again out of lack of time. The effects on the optimization in MPC can be explained as follows. Consider the optimization criterion, which was represented in Chapter 2
67
When minimizing this criterion, MPC will find that-for the input weighting we used for S T , which is 2-(large) changes in ST (Au) add more t o the optimization criterion, than they contribute t o reducing this cost by reduction of (weighted) output deviations from their setpoints (y - T ) . This is why MPC will not generate any changes in ST other than very small.
pcow
1 . 6 j
1.4
o;;.%]
0.35 ........... ...........
............... ...............
1000
2000
1000
2000
1000
2000
SB
ST
O. 06
r0.25t-l l
0.3
0.2
sv
i000
2000
1000
2000
1000
2000
Figure 5.5: Simulation 4 And indeed, Simulation 4 (Figure 5.5) shows that if we wait long enough, the throttle valve very slowly reaches the zone where it becomes effective again, and it is seen that in the end both pressure ratio and mass-flow reach their desired setpoint-it just takes about 10 times as long as it normally would have! The easiest solution t o this problem is to lower the weight on ST considerably, such that the throttle valve can move at maximum rate in order to reach effectiveness as soon as possible. The only danger in this is the fact that when setpoints are (almost) reached, the controller can become rather nervous (as will be seen in Simulation 6), since changes in ST are barely penalized.
5.2.5
Simulation 5
This simulation shows (Figure 5.6) that if we lower the weight on S T , the desired setpoint from the previous two simulations is reached within an acceptable time interval. We eventually chose such a weight on S T , that the throttle valve is not too restricted in its moves on the one hand, and will not act too nervous on the other (see Simulation 6). This resulted in a weight of 0.3, whereas other weights were unchanged.
does reach its setpoint, Apart from the fact that, thanks to a far more active S T , mcomp nothing else has changed. The turbine inlet temperature still doesnt exceed its maximum.
This simulation was actually designed for use in the next chapter, where it will serve as the example which the laboratory installation should follow, being actuated with the inputs we
68
CHAPTER 5. MPC APPLIED ON THE GAS TURBINE: SIMULATIONS
mcomp
Ttbin
100
200
100
200
100
200
SB
ST
SV
0.35I
0 OO .
100 0
200 5
~
0
100
200
1O0
200
computed here. Next, we will show the disadvantage of choosing an input weight as small as the one we chose for ST.
5.2.6
Simulation 6
In this simulation, we tried to decrease m, ,p from 0.4919 to 0.33 [kg/s], while maintaining a constant pcomp.As in Simulation 2 this describes a point for which Tbin by far exceeds its maximum. Therefore, it is not surprising that neither of the setpoints are met (see Figure 5.Q since the (linear) Ttain pretty soon reaches its constraint, once the throttle valve gets effective. At about 80 [SI, the linear Ttbin almost exceeds its maximum, after which it suddenly drops again, corresponding to the point where S B stops closing, since by then it is completely closed. Next, the throttle valve starts to oscillate with a rather small amplitude, which still results in quite severe oscillations in the linear turbine inlet temperature. This is caused by the small weight on ST. The weights on S B and SV in this simulation were chosen to be 10 and 100 respectively, in order t o prevent these valves (especially S V ) from oscillating as well, since this was seen (not shown here) to happen in earlier simulations. The problem depicted here, forced us to (from here on) choose setpoints which can be reached by varying the fuel valve, mostly. Therefore, in the following 4 simulations, where we wish to depict the effects of the various tuning parameters, we leave off from a simulation where the throttle valve is hardly used.
69
SB
ST
o::
sv
0.3 200
O
100
1O0
200

5.2.7
Simulation 7
In this simulation, we specified new setpoints for mcomp pcomp.The mass-flow had to and increase from 0.4919 to 0.64 [kg/s], the pressure from 164171 to 200000 [Pa]. Figure 5.8 shows
that the controller manages t o reach the desired setpoints quickly, although Ttbin momentarily (but rather severely) exceeds its constraint. This is again caused by the use of a constraint window.
As was expected, the controller mainly uses SV to reach the setpoints, but next t o that-at about 10 [SI-it opens the blow-off somewhat, in order t o stop the increase of pcomp,mainly. As was seen in Figure 3.14, the blow-off is highly effective in doing so.
Next, we vary the various tuning parameters. We compared the results from that to Simulation 7, but unfortunately didn't use different line types for the two responses. However, the accompanying descriptions should be able to make clear which is which.
5.2.8
Simulation 8
Simulation 8 (Figure 5.9) shows the effect of choosing m = 1 instead of m = 4. As was t o be expected, the controller is far less aggressive in its actions. All three valves are less active, and pcomp. Furthermore, Ttbin still exceeds its leading t o a slower response in both mcomp constraint, but not by as much as in the previous simulation.
70
C H A P T E R 5. MPC APPLIED ON T H E G A S TURBINE: SIMULATIONS
i I .6 .8
tj
O
SB
50
PCO~P
rncornp
50
Ttbin
O
ST
50
SV
50
50
50
5.2.9
Simulation 9
Next, we depict the effect of not using a constraint window (Figure 5.10). In this simulation, at sample 6, Ttbinwas about to exceed its maximum, after which constraints were being released from sample 7 t o 12. In this interval, Ttbin is seen to immediately rise, until (luckily) the setpoints are reached, which causes the fuel valve to close again. If it wasnt for this fortunate situation, Ttai, would most probably have exceeded its maximum for a much longer l period of time. Because of a l this, the three valves display something of a delay compared to Simulation 7, which causes setpoints to be reached slightly slower. As we mentioned in Subsection 5.1.5, while running the Garrett model, several times we experienced big problems to occur when a constraint window wasnt used. Introducing a constraint window was then often t o be found to solve the problem. When constraints are released, it can simply not be guaranteed that the controller will manage to meet the constraints again. It should be investigated, exactly how using a constraint window changes this situation.
5.2.10
Simulation 10
Instead of a control horizon of 4 samples, we now use 2 control blocks, each measuring 2 samples (see 5.11). Although the differences are rather small, it can be seen that indeed when we use these control blocks, all controller actions are less severe, resulting in a somewhat slower response in Ij2comp, especially.
71
rncornp
Ttbin
50
o
ST
50
50
SB
SV
04 R .5
04 .
09 .5
03 .5 50 0
50
5.2.11
Simulation 11
Finally, we wish to show the effect of choosing a prediction horizon of 30 samples instead of 15. When we did this starting from Simulation 7, no differences were to be seen. Therefore, we chose to use the situation from Simulation 5, and depicted the response of Ttbinin Figure 5.12.
Differences are still rather small, but it can be seen that when we take p = 30, the controller manages to keep both the linear and the nonlinear Ttbin exceeding the maximum, whereas from for p = 15 the Linear Ttbin in fact slightly exceed its constraint. did When p = 30 the controller can in fact see the point where Ttbin reach its constraint will right from the very start, since this happens at 19 [SI.
5.2.12
Unfortunately
...
When trying to show that repeated linearization possessed great benefits over Linearizing only once, something quite unfortunate was stumbled upon. This is depicted in Figure 5.13, where we ran the same simulation as in Simulation 5, apart from linearizing only once. We expected Ttai, to severely exceed its maximum because Simulation 5 was already close t o doing this. Quite to the contrary, the controller performed better than in Simulation 5, keeping Ttbina safe distance from its maximum, while reaching the setpoints even faster. at
72
C H A P T E R 5. MPC APPLIED ON T H E G A S TURBINE: SIMULATIONS
mcomp
Ttbln
50
Of= l
0.5 O
O
ST
50
50
Si3
sv
SO
50
Figure 5.10: Simulation 9 This is something we are still trying to explain. Although we might be dealing with just an unfortunate example of a case in which repeated linearization doesn't work too well, it suggests that something is wrong in the repeated linearization procedure, especially since linear and nonlinear model display smuller differences here than in Simulation 5. Lack of time kept us from finding out what exactly went wrong here, but we definitely recommend to look into this until the problem is solved. Still, in the Garrett model we experienced repeated linearization t o do nothing but good, which is why we feel that it should be able t o improve performance at (practically) all times.
5.3
Conclusions
Looking back at the simulations we presented in the foregoing, we can conclude that
o
Under normal circumstances, MPC is well capable of reaching setpoints, without exceeding constraints, which is probably the most important feature of MPC. Although MPC uses a linear model, it can still handle highly nonlinear processes, such as the gas turbine, especially if repeated linearization is used. This must of course partly be attributed t o the filter actions. However, repeated linearization was shown not always to improve performance. Since we do not (yet) really know what caused this, it doesn't seem appropriate t o draw too specific a conclusion out of this.
5.3. CONCL USIONS
73
mcomp
SB
0.005
0.5
o O
50 6
O *.O1 ~ O
50
SV
ST
50
50

o
If setpoints are specified that simply cannot be reached, without violating one or more
constraints, MPC partly ignores the setpoints but keeps the constrained outputs from (permanently) violating their limits.
Inputs that are hardly effective cause MPC to either become very slow in reaching the setpoints, if the corresponding input weight is comparatively large - become rather nervous when setpoints are almost reached, if the corresponding input weight is relatively small.
-
Such situations should therefore be avoided.

o
The rules of thumb for choosing tuning parameters provide a useful basis t o start from. However, the tuning process still involves trial and error procedures. Apart from that, we found the basic effects of the various tuning parameters to apply to the gas turbine as well (although not always as evident). A constraint window can be very useful, t o keep constraints from being permanently released. Our experience is that, with a constraint window, constraints are violated sometimes as well, but (almost) never for long periods of time. Finally, we mention that while running these simulations, MPC neeaed about 5 Es] to execute its computations (on a 486DX66 machine), if repeated linearization was used. Without repeated linearization, it could do with about 1.3 [SI, which is almost within the limit required for real-time implementation. However, execution time does not seem t o be much of a problem, since faster computers are quite commonly available.
74
C H A P T E R 5. M P C APPLIED ON T H E G A S TURBINE: SIMULATIONS
UV"
10
20
30
40
50
60
70
80
90
100
Figure 5.12: Simulation 11; - .-: p = 15, nonlinear;

--: p = 30, linear
: p = 15, linear; -: p = 30, nonlinear;
mcomp
Tbn ti
1O0
200
1O0
200
100
200
SB
ST
SV
1O0
200
100
200
0.26 0.24 0.22
1O0
200
Figure 5.13: The simulation that turned everything 'upside down'
Chapter 6
Experiments with PRIMACS

In this chapter we state the results of applying the inputs, obtained in Simulation 5, on the laboratory gas turbine installation. Out of safety reasons we chose a reasonably mild trajectory for the controller to perform. To actually apply off-line computed inputs t o the installation, it was necessary t o design a program which, via the data acquisition package LabVIEW, could drive the valves. These valves normally operate at maximum moving rate only, which is why moves have t o be translated to operating-time intervals. Through the manipulated variable move constraints, it is guaranteed that the valves will not have to move for any time longer than one sampling interval (which is approximately1 1 second). Still, when moves are smaller than some minimal value (determined by the resolution of the driving unit), they will have t o be accumulated until their total at least equals the minimum move size. All of this (and more) was implemented by Twan Althuizen, whose report should appear somewhere in February 1997 [Althuizen 971. Looking at the actual implementation of the inputs we obtained in Simulation 5, it can be seen that it is reasonably accurate (Figure 6.1). The resulting response is depicted in Figures 6.2, 6.3 and 6.4. In these plots, we filtered out the spikes, by simply not connecting these points. First of all, it can be seen that, a fair amount of noise is present. The measurement from which mass-flow somewhere near the throttle valve should be obtained (this is a pressure drop measurement), in fact still seems t o be too unreliable t o use. Instead, we reconstructed the mass-flow through the compressor by using the compressor characteristic. The result of this is depicted in (the left figure of) Figure 6.2. Furthermore, we mention that the turbine inlet temperature is measured by two thermo-couples. The smallest thermo-couple exhibits more noise, but responds much faster than the thicker one (Figure 6.4). Also depicted is the temperature response our model predicted, compensated with the 200 [KI we subtracted from
In the actual application of the inputs, the sample interval might have been somewhat longer than 1 second, due to (small) practical problems, which we will not describe.
75
76
C H A P T E R 6. EXPERIMENTS W I T H PRIMACS
Figure 6.1: Computed (- .-) and actually implemented (-)
input signals
2Ri.
Figure 6.2: meomp (left) and N (right); - .-: simulated; -: measured (.
in case of heomp)
the constraint (see Subsection 3.6.2 and Appendix B.2). The model reasonably well simulates the temperature response, apart from a much larger peak at about 20 seconds. Of course, it would have been worse if it were the other way around, since in our case, the controller just becomes somewhat too cautious. Finally, it is seen that, although the system responds almost the way we wanted, peompand meomp not drop far enough, which seems to indicate that in do reality, the throttle valve becomes effective even later than the one we used in our model.
77
Figure 6.3: peomp (left) and ptbin (right); -. -: simulated; - . - : measured
Ttbin 3mm en Ttbin, TEMPTURBIN
8ao -
..
720 700
v,
50
1O0 150
Figure 6.4: Turbine inlet temperature: - .-: simulated; -: measured, where the upper signal is generated by the smallest thermo-couple
Chapter ?
Conclusions & Recommendations

7.1
Conclusions
We have tried t o assess the feasibility of MPC as a control strategy for the gas turbine installation.
As we expected, controller execution time is an important issue, although not insuperable. With the computing hardware we used, it was not yet possible t o keep the execution time per sample beneath one second, if repeated linearization was used. In that case, the controller needed about 5 [SI to execute its computations. Without repeated linearization, however, this goal was almost achieved, with execution times of approximately 1.3 [SI. Still, when faster computers are available, this issue doesnt seem to be much of a problem anymore.
Next, we expected nonlinearities t o cause problems, since MPC uses linearized models for its computations. Even when we didnt linearize repeatedly, MPC performed rather well under most circumstances. On one occasion though, the model which was linearized only once performed even better than the model obtined by repeated linearization. We could not (yet) figure out why this happened, since we definitely expected the inverse t o have happened. However, we did experience the repeatedly linearized model t o improve performance in the simulations we did with the Garrett model. This causes us t o believe still should be worth its while. Still, this is something which needs t o be further investigated. We derived a lumped parameter model of the gas turbine which reasonably well described the real installation. Still, this model (for the BBC configuration) should be improved t o better describe the response t o variations in the throttle valve. This same valve caused some problems, due t o the fact that it is hardly effective when operated somewhere near completely opened. The model which described the Garrett configuration on the other hand described the corresponding process surprisingly well. The reason for the difference in model match between 78
7.2. RECOMMENDATIONS
79
the two still isnt clear. Looking back at what was caused by all this, we feel that the reconfiguration from Garrett t o BBC was (unwillingly) perhaps a little more significant than it was initially expected t o be. All in all, the models we obtained seem t o be a fine compromise between accuracy and simplicity, which was needed for our purposes. The specific advantages MPC claims t o possess over other control strategies can be said t o hold up t o some extend. In most cases, MPC manages t o reach certain setpoints while keeping coristrairits from Seir,g violated. Howeve, sometimes when constraints are aibmt t o be exceeded the problem may become infeasible, if hard constraints are present. These constraints may never be released, whereas soft constraints can in fact be (temporarily) released, t o prevent the optimization problem from becoming infeasible. It can still not be guaranteed that this release of constraints will in fact be temporary. The use of a constraint window may cause constraints t o be exceeded as well, but our experience is that these violations more often than not last for only a short period of time. One of our recommendations stems from this problem. In Chapter 6 we applied an off-line computed sequence of inputs t o the laboratory installation, and compared this response to the corresponding simulation. The results were about as good as we expected them to be, since the model did not describe throttle valve response very well. From what was seen in this one experiment, we feel that the implementation of MPC on the real installation seems rather promising, also because input constraints and input move constraints were correctly incorporated. Secondary t o this, we implemented the model in PRIMACS, which in itself was rather laborious. A new discretization routine had to be programmed in order t o yield a stable model. Apart from that, we had t o reconstruct the states that werent available as measurements. Understanding the source-code used in PRIMACS seems rather important, if the various operations are t o be carried out properly.
9.2
Recommendations
Finally, we wish t o recommend a few things for future work in this area. Most of these recommendations leave off from the conclusions we stated in the above. Starting with the first (and probably most trivial), we recommend the use of computing hardware which is about five times as fast. Since we used a 486DX66, this does not seem t o be much of a problem. Secondly, the model we used for the BBC configuration should be improved (if possible), t o better describe throttle valve response. This seems to be easier said than done, since all sorts of things have already been tried t o achieve this. Our third recommendation actually was done by Frits van der Meulen in his report [Van der Meulen 951 and regards the use of an automated process which decreases the
80
C H A P T E R 7. CONCLUSIONS & RECOMMENDATIONS
size of a constraint window, if constraints are in fact about t o be released. However, if our fourth recommendation is complied with, the third one may no longer be necessary at all. This recommendation, which is in fact already taken into consideration (having been stated earlier in [Satter 96]), regards the penalization of constraint violations. If violations are somehow incorporated into the quadratic optimization criterion (2.1),the controller will be forced t o reduce them, by under-performing in trying t o reach the setpoints. The fifth recommendation we wish t o make, suggests the use of more efficient integration routines in the PRIMACS source-code. At the time of writing, an Euler forward integration scheme is made use of. We finish with mentioning some of the work that should preferably be carried out by a successor. Apart from the recommendations stated here, a successor should look into the implementation of a(n extended) Kalman filter, once more. If noise is t o be coped with, a real-time implementation can probably hardly do without something like it. First of all though, attention should be paid to the problem we encountered in Chapter 5, which indicated that the procedure of repeated linearization as it is currently implemented, is certainly not without faults. Next, it would be very interesting t o compare the performance of MPC-in simulations as well as in experiments-with other control strategies. Meanwhile, it will be useful t o stay in touch with other work done in this area, especially if this also involves the use of PRIMACS. The soonest publication in this context will most probably be the work done by Gert Leenheers [Leenheers 971.
Bibliography
[Althuizen 971 T. A. M. Althuizen, Klepaansturing van een StIG-installatie (?), Technical Report, Eindhoven University of Technology, 1997 ( t o appear) [Bequette 911 B. W. Bequette, Nonlinear control of chemical processes: A review, Industrial Engineering Chemical Research, vol. 30, no. 7, pp. 1391-1413,1991 [Borowski and Borwein 891 E. J. Borowski and J. M. Borwein, Dictionary of mathematics, Collins (Reference), 1989 [Cohen et al. 871 H. Cohen, G. F. C. Rogers and H. I. H. Saravanamuttoo, Gas turbine theory, Longman, Scientific k Technical, 1987 3 [Econosto 911 Algemene Econosto catalogus, Cappelle aan de IJssel, 1991 [Van Essen 95a] H. A. van Essen, Design of a laboratory gas turbine installation, Technical Report WOC-WET 95.012, Eindhoven University of Technology, March 1995 [Van Essen 95b] H. A. van Essen, Centrifugal compressor; principles of operation and control, Technical Report, August 1995 [Van Essen 95c] H. A. van Essen, Model Predictive Control (MPC)-Report ciples and algorithms, Technical Report, November 1995 on basic prin-
[Van Essen 961 H. A. van Essen, Compressor modeling and control, Technical Report, 1996 [Greitzer 761 E. M. Greitzer, Surge and rotating stall in axial flow compressors-Part I: Theoretical compression system model, Journal of Engineering for Power, vol. 98, pp. 190198, 1976 [Haarsma 951 H .Haarsma, Robust model predictive oxygen control: The development of a robust model predictive dissolved oxygen controller for an activated sludge process, Final Project, Vakgroep Meet-, Regel- en Systeemtechniek, Landbouw Universiteit Wageningen, 1995 [Kok 901 J. J. Kok, Werktuigkundige regeltechniek 11, collegedictaat werktuigbouwkunde, Eindhoven University of Technology, 1990 [Lahoye 961 M. Lahoye, Experimentele verificatie van de karakteristieken van een radiale uitlaatgasturbine, Technical Report WOC-WET 96.028, Eindhoven University of Technology, 1996
81
82
BIBLIOGRAPHY
[Van der Meulen 951 F. van der Meulen, Toepassing van lineaire model predictive control op een thermisch hydraulisch proefproces, Technical Report NR-1928 (1995-11-6), Eindhoven University of Technology, 1995 [Muske and Rawlings 931 E(. R. Muske and J. B. Rawlings, Model predictive control with linear models, AIChE Journal, vol. 39, no. 2, pp. 262-287, 1993 [Leenheers 971 G. J. C. C. M. Leenheers, Model predictive control of a compressor station (?), Technical Report, Eindhoven University of Technology, 1997 ( t o appear) [Maple 941 Maple V Release 3 for DOS and Windows-Getting Software 3, Springer- Verlag, 1994 [Matlab 921 Matlab Reference Guide, The Math Works, inc., 1992 [Morari and Ricker 931 M. Morari and N.L. Ricker, Model predictive control toolbox (MPCTools)-Matlab functions for the analysis and design of model predictive control syst e m ~ , 1993 ~ [Morari et al. 911 M. Morari, C. E. Garcia, J. H. Lee and D. M. Prett, Model predictive control, Preniice-Ball, 1991 (still to appear) [Peeters 951 J. J. F. A. Peeters, The behavior of a MPC controller and the the implementation in PRIMACS, Technical Report, Eindhoven University of Technology, 1995 [Satter 961 I. H. G. Satter, Model Predictive Control (MPC) in de PRIMACS-omgeving, Technical Report 72-01-28-728-020, Landbouw Universiteit Wageningen, 1996 Started, Waterloo Maple
Appendix A - -
MPC
A. 1 Unconstrained Model Predictive Control
The objective was stated as
The control problem can now be expressed as the following optimization problem
AU(k)
min { 11 rY[y(k 1Ik) - ~
+ k11 112 + II ruAu(k)It'} 1
(A4
such that
Y ( k t 1Ik) = h Y ( k I k ) + SdAd(L) S u A U ( k ) t
where the following definitions are used
I" = diag{I'y,. . . ,r;} '
ry= diag{I'l,. . . ,I?;}

Y(k
+ lp) =
R ( k f 1) =
83
84
APPENDIX A. MPC
AU(k
M=
0 9
1 6
.
I
.
............ 0 C .........
I
o ...... o
. .
a
. .
... O
I I/
'I
In equation (A.2) the first two' terms are completely defined by past control actions and present measurements and the last term describes the effect of future manipulated variable moves. This equation originates from the superposition theorem stating that response = effect of initial condition
+ effect of manipulated variable move
So the first term in equation (A.2) corresponds to effects which can be attributed t o initial conditions without moving the input at/after k - 1. Since it is assumed that the system will settle after some time, say n time intervals, M has n columns, whereas the prediction only looks p time intervals into the future corresponding t o the p rows in M . Actually, Y ( k l k ) is the corrected vector of future outputs, where KF is the filter gain and y) , ( k the measurement at time 5 Y ( k l k ) = Y ( k l k - 1) K F ( y m ( k ) - y(klk - 1) )
Furthermore, A d ( k ) is the measured disturbance vector. S d and SU are step-response coefficient matrices, completely describing how a particular input (or disturbance) affects the system output when the system is at rest. Note that the first column of SU shifts downwards in the subsequent columns t o account for the fact that the system is assumed t o be strictly causal, i.e. the output at a certain time instant only depends on the input up t o that time instant, exciuding the time instant itseif.
'Ad()) is assumed not to change in the future, for in general its exact future behaviour wili not be known a priori.
A.2. CONSTRAINED MODEL PREDICTIVE CONTROL
85
To solve the optimization problem (A.l)-(A.2), the following set of linear matrix equations is formed
(A.3)
or
with
E p ( k l [ k ) = R(k t 1) - MY(k1k) - S d A d ( k )
The solution that solves the optimization problem equals the least squares solution to matrix equation (A.4) and is given by
au(k)= (suTryTrysu r u T r u ) - i s u T r y T r y q kt i l k ) +
and the resulting control law becomes
au(k) =
[I o
... o
(SuTryTrys"+r21Sru)-1SuTrySryEp(k+ilk)
=: I c M P C q , ( k + i l q
where only the first move (at time k ) is computed. Normally spoken, K M p c can be computed off-line, since the only time-varying element is the projected error vector Ep(k I l k ) . Only when y , I" and/or SU need to be updated, K M p C must be recomputed. SU would have to ' be updated, for instance, if repeated linearization was used, as will be done in Chapter 5.
A.2
The three types of process constraints will have to be formulated as linear inequalities in order to solve the constrained control problem as a Quadratic Program (QP) which can generally be stated as min x T H x - gTx
X
such that where
Cz 2 c
H is denoted as the Hessian matrix,

g is the gradient vector,
C is the inequality constraint equation matrix, and
86
c is the inequality constraint equation vector.
APPENDIX A. MPC
So, a manipulated variable constraint at a future interval of time k

Ulow(k
+ Z of the form
+ Z) 5
AU(k
j=O
+j ) t
U(k -
1) 5
uhigh(k
+I )
I = o, 1,.. . ,m - 1
can-when
augmented with all other projections-be
expressed in matrix inequality form as
u@ - 1) - % h i g h ( k )
u(k - 1) - u h i g h ( k m - i) urow(k) u(k. - 1) UlOW(k
tm - 1) - u ( k - 1)
where
I, =
I o ... O 7 I i ... O . . . . ... . . . . I I

* * '
I can be
Rather similarly, manipulated variable rate constraints of the form /Au/ Au,,, 5 expressed in matrix inequality form as
-Au,,,(k where Au,,,
+ m - 1)
> O is the bound on the magnitude of the moves.
Finally, when using equation (A.2), output variable constraints can be formulated as
where
Ylow(k
+ 1) =
A.2. CONSTRAINED MODEL PREDICTIVE CONTROL With this, the Quadratic Program becomes
AU(k)
87
min { 11 rY[S"nu(k) such that
+ Ilk)] 112 + II r"[au(ql112}

C " A U ( k ) 2 C(k
+ lik) +
where (A.5), (A.6) and (A.7) are combined into one matrix inequality. C" and C ( k i l k ) collect all left hand sides and right hand sides of the inequalities, respectively, yielding
and
C(k
+ llk) =
utow(k m - 1)- u(k - 1) -A&iXX(k)
Converting this to the standard QP formulation, we get min -AU(k)T'FI"AU(k) 1

AU(k)
2
C"AU(k) 2 C ( k
T p " r y p
such that with Hessian
+ lik)
3-1" = p
and gradient vector
+ yuTy
q k + ilk) = s u T r y T r w P ( k ilk). +
This Quadratic Program will be solved at each controller execution time. In principle, the Hessian 'FI" can be computed off-line as was the case for the computation of K M P Cwith , the same exception when weights or model need t o be updated. In all other cases, the only time-varying elements are Ep(k i l k ) and C ( k i l k ) . Then a parametric QP algorithm which employs the pre-inverted Hessian in its computations is preferable in order t o reduce on-line computation effort.
88
APPENDIX A. MPC
A.3
The filter gain

K F
In [Morari and Ricker 931, the first-order filter gain vector of future outputs Y ( k I k ) as
is used t o determine the corrected

- y(klk -
Y ( k I k ) = Y ( k l k - 1)
+ IC&&@)
1))
(A4
where Y(K1K - 1) denotes the estimate of Y ( k )ebtained at k - 1, and bases o2 rileasfirenlent information up t o k - 1. Here Y ( k )is the vector of future outputs when at time k - 1 we stop varying all inputs AV (consisting of manipulated variables disturbances) t o a process
Y ( k ) :=
=
[ ~ ( k ) ' ,y(k
[ ~ ( k ) ' y(- ,) . . . , y(k n -, , 'I ' y(k ) 2 for Av(k i ) = O ; i 2 O
+,')I
... , y ( k + n -, ' ) L
+n
-
i)']'
y(k
+ n ''2 ])
since al outputs after k l settles).
+ n - 2 are assumed to be constant (after n time intervals the system
A.4
Matlab versus PRIMACS
Matlab uses the discrete-time (linear time invariant) state-space representation:
.(k)
= az(k - i)
=
Y(k)
m+m C.(k) + D,u(k) + Ddd(k) + D,w(k) + .(k)
+rUu(k
i)
+ rdd(k - i) + r,w(k
i)
In this representation d and w denote measured- and unmeasured disturbances respectively, and z represents measurement noise. The outputs without measurement noise are denoted by
YAlthough D, and D directly connect u ( k ) and d ( k ) (respectively) with y(k), these matrices d may not contain any elements other than zero. So, in fact, this is the same as saying that
Y) = C.@) @
+ D,w(k) + .(k)
-
with the exception that for some reason D, and Dd have to exist and be specified in the right dimensions. The result is that y(k) only depends on inputs up to, and including, sample ( k from w(k) and z ( k ) , of course-yielding a strictly causal model. 1)-apart
In PRIMACS, the same result is obtained, but in a different way. This time, the D-matrix does not have t o consist of zeros, but the discrete-time representation is defined differently instead:
x(q =
y(k)
=
az(s - 1) + r u ( k - 1)
C.(k)
+ Du(k - 1)
A.4. MATLAB VERSUS PRIMACS
89
In this representation, where we left out2 d, w and z , its clear that y(k) does not depend on inputs from sample k but only on inputs older than that. Still, one can specify a D-matrix containing elements other than zero, which can be very useful.
In [Van der Meulen 95, Appendix A], more on this subject can be found. In that appendix, Van der Meulen also pays attention to implementation related differences in disturbancemodeling, as well as in state-estimation. Apart from this, Matlab and PRIMACS also slightly differ in their definitions of weighting factors. We will not $0 into any of this any Farther.
2This is done in analogy with the description in [Satter 961, from which this was taken. This does not mean that d , w and z arent present in the PRIMACS model.
Appendix B
Gas turbine model specifics

In this appendix all fits that were used for both the Garrett and the BBC configuration are listed, starting with the Garrett.
B.l
The Garrett configuration
Starting with the compressor characteristic, we used the following second order fit for the compressor pressure ratio as a function of mass-flow mcGmp rotational speed N ) : (and
where Scale = "Pt," 43::778 0.4536 scales the mass-flow t o the appropriate dimensions' and 6o Acomp, BcGmp Ceomp third order fits on the rotational speed: and are
AcGmp= ((-2.788188. 10-'sNrpm 3.64473. 10-14)N,pm
+ 1.879789. 10-8)NTpm
9.490669 = 2.654884 * 10-6)NrpmBcGmp ((6.263607. 10-'6N,pm - 7.269489 * l0-")Nrpm 0.0229777141 Ccomp = (( -2.310107. 10-14NTpm 3.663545 * 10-9)N,pm - 1.723704 10-4)NTp,
3.7372517755 with Nrpmin [rev/min]. N as used in our model is defined in [rev/s], which is why we have bothered to use the subscript here-so, Nrpm= 60.N. The resulting compressor characteristic is depicted in Figure B.1. More specifics about this fit can be found in [Lahoye 961. The real compressor characteristic (Figure B.7) displays lines of constant compressor efficiency. We did not pursue t o fit this efficiency, but assumed it to be constant and equal t o 0.70.
*Scale corrects for deviations from the standard conditions as well as converts from [kg/s] to pbs/min]. The latter (non-metric) dimension is used in the Garrett compressor characteristic.
90
B.I. T H E G A R R E T T CONFIGURATION
91
10
20
30
40
60
70
niascallw [l&min]
Figure B.l: The (Garrett) compressor characteristic fit (from [Lahoye 96]), unfortunately accompanied by Dutch captions If mcomp to be obtained from the static equation is
m c o m p - fcompressor(N,Peomp,Pin,
Tin)
fcompressor
will have the following form (where N is omitted for clarity)
fcompressor
..
2AcOmp Scale
. .
which simply is one of the two2 solutions to quadratic equation (BA). Next, to determine mass-flow through a valve, the following two-empirically obtainedrelationships for gaseous fluids are used [Econosto 911. They consecutively describe sub- and supercritical flow, the last case of which describes the situation when fluid velocity has reached its maximum value (at the smallest area): subcritical flow:
E 2 0.53
supercritical flow:
m = 0.0711i-~pl&
E L: 0.53
JB:,,,
B + ~ ~-4 A~~~~ ( cComp P i n 1 ~ ~ - com mg 2The other solution, fcompresaor = corresponds t o mass-flows in the 2AcompScale area left of the surge-line. Although we mean to stay out of this area, this second solution is used to describe system behaviour left of the surge-line to some extend.
92
APPENDIX B. GAS TURBINE MODEL SPECIFICS
where Tiand pi are temperature and pressure before the valve and p2 the pressure after, A p the pressure difference ( p i - p 2 ) , ICv is the valve capacity [m3/hr] and pn the normalized density (at 1.013 [bar] and 273 [KI). Pressures should be entered in [bar] and temperatures in [KI to get k in [kg/s]. The valve capacity ICv is the momentary value which means that it varies depending on the valve position S. This dependency can, for instance, be proportional
with S E [O, i]and ICvma, the maximum valve capacity. Also rather common is an exponential valve charact eristic ICv = ICvm,, but other functions can be imagined as well. Figure B.2 depicts the two characteristics suggested here.
Figure B.2: Proportional and exponential compressor characteristic (base = 20)
In our model we only used proportional valve characteristics, both for the blow-off and the throttle valve, where for the Garrett we have
Kvblmaz 160 [m3/hr] =
and
ICvthrmaz 400 [m3/hr]. =
For the turbine, we have the characteristic depicted in Figure 3.9. We used the following fits for mass-flow through and efficiency of the turbine
ktb =
((((((-8.345\ \ \ \
litbout
+ 110.542) 1
- 604.202) -+ 1748.819)
Ptbout Ptbout
Ptbin-
1 Ptbout
2840.761) - 2475.591) -- 870.560) -/88.30.00756 tbin Ptbin ptbin 101321 Ttbin Ptbout Ptbout
B.2. THE BBC CONFIGURATION

and qtb == ((o.0414- Ptbin &bout
- 0.3260)
93
Ptbin + 0.7500)
Ptbout
*
Ptbout
f 0.1658
Next, we fitted Power t o the fuel valve position S V . We did this by determining the power needed t o approximate a measured stationary point (for each point separately) as good as practically possible, fdowed by 5tting these v2lues GE the corresponding (measured) fuel valve positions S V . Of course, this fit largely depends on the choices made for the overall friction factor (=1.0), the valve capacities K v b l m a s (=160) and K v t h r m a z (=400), and the efficiencies vcomp (=0.70) and v t b (fit). Anyway, after some iterations, we arrived at the following fit Power = e13,421718-0.756721 where SV E [O, i]and Power has the dimension of [W]. As could be seen in Figures 4.6 and 4.7, this fit works rather well. Based on simulations, we decided t o take a maximum turbine inlet temperature (for the Garrett) of 900 [KI (see Subsection 3.6.2 for motivation). Fitting the transient response was done by varying Vcomp, Veeand I. We ended up with Vpl, the following choices
Vcomp 0.5 [m3] =
V,,
vee
= 5.5 [m3] = 2.0 [m3]

=
(J3.3)
0.001 [kgm]
All other parameter choices will be stated together with those used for the BBC configuration in the table at the end of this appendix.
B.2
The BBC configuration
Here, we state the fits used for the BBC model, where we will stay with the order used in the previous section.
So, we start with the compressor fit, which is similar t o the previous one, apart from the fact that coefficients are different:
pressure ratio = Peomp
-= Pin
Aeomp(N)(meompSCale)
+ Beomp(N)&mpScale + Ccomp(N)
(3.4)
94 speed
Acomp = 2.77309428403675 - lo1 - 1.15180363512399 10-2Nrpm 1.57304616982506- 10-6N,pm - 1.058110041784 10-loN,pm

4.36553159400385 10-25Nfpm
+
-t
3.72374699463208 ~dotlO-'~N,P, - 6.51912048495668 10-20N,!pm
Beomp = -1.45197329U98387 i 2 4.04918186433896 - iO-21?rp, 4.52872638270511 - 10-6N:pm 2.60954693821889. 10-loN,pm 8.11642530137263 - 10-15N,pm 1.28243058048326* 10-19N~'m7.88466994246646 - 10-25Nfpm C e o m p = 7.25717433841597 - 10' - 1.86190895383627 10-2Nrp, 1.95883720952405* 10-6N,?pm- 1.06175287530094* 10-loN~p, 3.10922060166141 10-15N:pm - 4.61132963429711- 10-20N,?pm 2.63995234246694 10-25N,6,m
+ +
+ +
Again Nrpm denotes the rotational speed in [rev/min]. This fits of these coefficients are depicted in Figure B.3. The resulting characteristic is compared to the experimentally determined compressor characteristic3 in Figure B.4. The compressor efficiency is again assumed t o be constant: qcomp 0.77. =
(:a2), Bcomp (:al) and Ccomp (:aO) Figure B.3: The fits of coefficients Aeomp
If meomp to be determined statically, equation (B.2) should be used. is

3The real BBC characteristic is SO depicted in Figure B.5 at the end of this appendix.
B.2. T E E BBC CONFIGURATION
95
0.1
02 .
03 .
0.4
0.5
0.6
07 .
0.8
09 .
Figure B.4: The (BBC) compressor characteristic fit; experimental measurements are denoted by + Next , we still used proportional valve characteristics both for blow-off and throttle valve, 3 where l<&[mac = 160 [m /hr] and IiiJthrmax = 550 [m3/hr]. The turbine characteristic is depicted in Figure B.6 at the end of this appendix. We used the following fits for mass-flow through and efficiency of the turbine
mtb
= 10
((((-0.006181-
ptbout
+ 0.088298)
Ptbout ptbin
0.476851) Ptbout @!E
+ 1.167342) Ptbout &+
0.082815)21.1-and
?)tb =
Ptbin Ptbout Ptbout
Ptbout
lo5
TinRT
((-0.00048-
Ptbin Ptbout
- 0.003704) @!E - 0.009924)

Ptbout
ptbout
Ptbin
+ 0.788052
The maximum turbine inlet temperature we decided on was 700 [KI, based on the considerations mentioned in Subsection 3.6.2. Finally, we fitted stationary points and transient responses as best as possible in the same manner as before. This yielded the following fit
where SV E [O, 1 and Power has the dimension of [W]. 1
96
Furthermore, we have
Vcomp = 0.3 [m3] V,, = 2.5 [m3] v,, = 0.5 [m3] 1 = 0.01 [kgm2]
Other parameters used in this (BBC) model can be found in the following table, together with the parameters from the Garrett model. E
1.0
Vcomp
v p ~ vcc
I
0.001 0.01
YT
1.4 1.375
RT
287 291.5
Kvblrnaz
Kvthrmaz
Tcornp
Cpcomp
Garrett
BBC
0.5 0.3
5.5 2.5
2.0 0.5
160 160
400 550
41.7
0.70 0.77
1000 1010
1040 1060
We conclude this appendix with remarking that the Garrett model approximates reality far better than the BBC model. For the BBC model we did not really succeed in fitting the stationary points where the throttle valve was partly closed, although we most certainly tried. We investigated the separate influences of Kuthrmaz,, qcOmp,b and Power(SV), but E qt found that neither combination yielded the desired results. It still isnt clear why this was so hard compared to the Garrett situation. It is therefore recommended that future extensions to this project look into that once more.
B.2. T H E BBC CONFIGURATION
97
Figure B.5: The experimentally obtained BBC compressor characteristic
98
1.35
i3 0
:.25
1.21
L.IS
:.Lo
ai:
I.
:.ia
Figure B.6: The experimentally obtained BBC turbine characteristic
B.2. THE BBC CONFIGURATION
99
Figure B.7: The experimentally obtained Garrett compressor characteristic
Appendix C
MPC implementation specifics

In this appendix we state PRIMACS source-code associated with implementation specific problems, treated in Chapter 5. We start by listing the source-code corresponding t o the discretization routine.
// --- Discretisation with Taylor (plus scaling and unscaling by // --- squaring)
// --- Discretisation of matrix A; first a matrix // --- A = [ L BLI*t; zeros(NU,NX+NU)I is constructed B[A
NT = NX + NU; AB.init(PDIM(NT,NT>>; for (i=l;i<=NX;i++) for (j=;j<=NX; j++> AB(i,j) = msample*AL(i,j); for (i=l;i<=NX;i++) for (j=NX+I ;j<=NT; j++) AB(i,j) = msample*BL(i,j-NX);
// //
--- Scale AB by 2 --- 1s < 1 / 2 . .
so
that its infinity-norm (largest row-sum)
infnorm=O; sum = o; for (i=l;i<=NT;i++)
c
sumold = infnorm; for (j=,sum=O;j<=IT;j++) sum += fabs(AB(i,j)); infnorm = --max ( sumold, sum) ;
>
100
101
mantissa = frexp(infnorm, &exponent); s = --max( O ,exponent) ; AB = AB*(l/pow(2,s));
//
--- Taylor series or exp(A)
E.init(PDIM(NT, NT) ) ; F = eye(NT);

k = I; RES = (E + F) - E;
onenorm = O; sum = o; for (-j=i; j<=NT;j++)
c
sumold = onenorm; for (i=l,sum=O;i<=NT;i++) sum += fabs(RES(i,j)); onenorm = --max(sumold,sum);
1
while (onenorm > O)
c
E = E + F ; F = (AB*F); kreal = k ; F = (i/kreal)*F; k = k + l ; RES = (E + F) - E; onenorm=O;
sum = o;
for (j=l;j<=NT; j++)
c
sumold = onenorm; for (i=i,sum=O;i<=NT;i++) sum += fabs(RES(i, j)); onenorm = --max(sumold,sum);
3 3
// --- Undo scaling by repeated squaring
for (k=i;kc=s;k++) E = E*E; S = E;
// --- Get AL and BL from S; CL and DL are unchanged
102
f o r (i=i;i<=NX;i++) f o r (j=i; j<=NX;j++> AL(i,j) = S(i,j); f o r (i=i;i<=NX;i++) for (j=NX+i; j<=NT;j++) BL(i,(j-WX)) = S(i,j);
APPENDIX C. MPC IMPLEMENTATION SPECIFICS
DL &= DL;
In this code we used the following definitions:
AL, BL, CL, DL NX, NU, NT msample AB

S
matrices A , B , C and D obtained by linearizing number of states, inputs and their total, respectively discretization sample interval the NT x NT matrix which combines AL and BL (augmented with zeros) the Ciiscrete version of AB denote PRIMACS functions (such
Other undefined names-if not self-explanatory-usually as PDIM), or C++ functions (e.g. fabs).
Next, we state the procedure which reconstructs the states, needed for linearization.
where yp (=[licomp Teompppl Tpl l ) t b i n Ttbin N I ) denotes the scaled output(measurement)s and x-old the states that will be used to compute the linearization. Furthermore, Y . scale and Y .unscale are (PRIMACS) functions that scale and unscale the outputs, respectively, by using equation (5.1). In this notation, Y .unscale(yp(L) ,4), instance, means for the unscaling of yp(2) with amplitude and offset of yp(4).
Appendix D
MPC parameter settings

In this appendix, we represent a table with the parameters as they were used in Chapter 5. The simulation numbers we used here, can also be found in that chapter. Unless stated otherwise in the colilma Remarks, we linearized on each sample. Other standard settings were stated in Chapter 5.
103

Model Predictive Control of A Gas Turbine

Uploaded by

Copyright:

Available Formats

Model Predictive Control of A Gas Turbine

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Model Predictive Control of A Gas Turbine

Uploaded by

Copyright:

Available Formats

Model Predictive Control of a Gas Turbine Installation

Bas Vroemen WFW report 97.002

To the memory of Tante Liesbeth

Model Predictive Control of a Gas Turbine Installation

Final project Supervisor: Harm van Essen Professor: J. J. Kok

1.1 Gas turbines

1.2 Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3 Goals of this project

1.5 Organization of report . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Unconstrained Model Predictive Control

Constrained Model Predictive Control

2.3.4 2.3.5 2.3.6 2.3.7

Input- and output weightings . . . . . . . . . . . . . . . . . . . . . . . The filter gain

A constraint window . . . . . . . . . . . . . . . . . . . . . . . . . . . . Control blocks

2.4 Alternative MPC methods 2.4.1 2.4.2

Repeated linearization Nonlinear MPC

2.5 MPC implementations: Matlab and PRIMACS

3 The gas turbine installation

3.1 Construction and operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The compressor

3.3 The combustion chamber 3.4 The turbine

Inputs and outputs Constraints

4 Modeling the gas turbine

4.1 The compressor model

4.2 Extensions t o the basic compressor model

4.3 Modeling the complete gas turbine . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

M P C applied on the gas turbine: simulations

Discretization Reconstruction Constraints

Simulation 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 6

Simulation 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2.11 Simulation 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.12 Unfortunately . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Experiments with PRIMACS

7 Conclusions & Recommendat ions

7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Recommendations

A.4 Matlab versus PRIMACS

B Gas turbine model specifics

C MPC implementation specifics

T fuel bar max Pl cc amplitude offset

Model Predictive Control

1.3. GOALS OF THIS PROJECT

Goals of this project

Unconstrained Model Predictive Control

CHAPTER 2. THEORY: MPC

2.2. CONSTRAINED MODEL PREDICTIVE CONTROL

Constrained Model Predictive Control

1. the sample time

3. the control horizon m

4. the input weightings I?"

CHAPTER 2. THEORY: MPC

The sample time

The prediction horizon

2.3. TUNING PARAMETERS and vice versa;

increase the optimization criterion.

The control horizon

and vice versa;

2.3.4 Input- and output weightings

2.3. TUNING PARAMETERS

The filter gain

CHAPTER 2. THEORY: MPC

with IC(k1k - 1) the estimate of .(IC), t o k - 1.