- Email: [email protected]
Optimization and Control of Combustion Process Products
Copyright @ IFAC Power Plants and Power Systems Control, Brussels, Belgium, 2000
OPTIMIZATION AND CONTROL OF COMBUSTION PROCESS PRODUCTS Daniel Pachner .,1 Vladimir Havlena •• ,1
• Trnka Laboratory of Automatic control, Faculty of Electrical Engineering, Czech Technical University e-mail: [email protected] .. Honeywell Technology Center Prague, Pod vodarenskou vezl 4, 182 08 Prague 8, Czech Republic e-mail: [email protected]
Abstract: A practical adaptive control application is reported. In coal-burned thermal plant, the air supply rate is controlled using a simple non-linear model of carbon monoxide generation. Parameters of the model are on-line estimated using the standard least-squares estimation technique. This model is used for the the optimal air/fuel ratio re-calculation, which should ensure that the steady state CO concentration would exceed a given limit with only a small probability (cautious optimization). Both the limit and probability threshold are the two principal input parameters. To avoid covariance wind-up, the regularized forgetting algorithm is used. To model the fact the parameter variations are caused mainly by the fuel quality fluctuations, the forgetting is directional. The output data are provided by CO concentration sensors, which may get saturated or clogged. Mathematical solutions to both problems is presented. Copyright @2000 IFAC Keywords: power plant; thermal plant; combustion control; cautious optimization techniques; carbon monoxide emission control; failure detection; least-squares estimation; regularized forgetting; sensor saturation; covariance wind-up
1. INTRODUCTION
2. MODELS OF COMBUSTION PRODUCTS GENERATION
When coal is burned, it is desirable to adjust the air supply rate (excess air) in an adaptive manner, so that the changes offuel quality would be compensated for. Adaptive air supply rate allows to stabilize the carbon monoxide concentration bellow a prescribed limit, which also automatically ensures good efficiency of the combustion process (Hanby 1994).
A model relating the air and fuel supply rates to the controlled concentrations is a necessary prerequisite for the control. The structure of the model cannot be however complicated, because its numerical parameters are thought to be estimated by a reliable on-line method. The models developed on the basis of deeper knowledge of the combustion processes and its physics always contradict this requirement. The solution presented is therefore based on a statistical
1 This research was partially supported by the Grant Agency of the Czech republic, grant 102/97/0466, and by the Ministry of Education of the Czech Republic under Project VS97/034
127
~ logO
of the model (1) would fluctuate near the steady-state characteristics expO
(2)
.!.=.!!.
(The symbols M and V denote the mean value and variance respectively.) Provided the function rj> can be approximated by its tangent, the steady state variance can be roughly approximated by
logO %-Q
FIR!
FIR2
V(y)
i
i
wet)
e(t)
u(t)
= V(e)
drj>(x) dx
2 1
I/l-1M(II)
(1 1-
0'~2.
(3)
0'
Thus, the steady state characteristics is a stationary random process and its statistical properties are functions of inputs uo, Wo. Apart from e(t), the other source of uncertainty are the model parameters, which must be also considered to be random variables.
~ yet)
Fig. 1. A model of the CO generation process: its concentration y(t) is a function of filtered the air w(t) to the filtered fuel u(t) supply rates.
The concentration of carbon monoxide exhibits a nonlinearity which can be hardly approximated by its tangent, but should be better approximated by exponential function rj>(x) = exp(cI x + C2)' The constants Cl, C2 must be calculated numerically. To predict mean value of the steady CO concentration, it is better to approximate the output nonlinearity by an exponential and use the well-known formula
model, which captures only statistical characteristics of the dependency instead of modeling it in detail. The remnant difference between the model and the true process is captured by a fictitious random signal. This is the standard way how to tackle overly complicated processes. As it is known that the steady state concentration of any of the combustion products (CO, O2 , NO:t) is a function, say rj>, of air to fuel ratio, the model should comply with this item of information. The following model was suggested:
This way of approximation decreases the bias of CO prediction and should be used rather than (2). (1)
3. PARAMETER ESTIMATION
The symbols have the following meaning:
The model (1) contains a number of unknown parameters: ak,bj,c,O',rj>. It is considered safer not to estimate 0' to avoid unstable (0' > 1) models. The value of 0' is therefore set equal to an off-line estimate. The function rj> can be set equal to a steady state characteristic, which can be obtained experimentally. With the parameters (J = (ak,bj,c) estimated, the model captures the following class of steady state characteristics
y(t) - Product concentration (e. g. CO) u(t) - Fuel supply rate w(t) - Air supply rate
e(t) - Fictitious noise rj> - Steady state characteristics 0' -
Autoregressive coefficient
ak, bj - Moving average coefficients
y(ak,bj,c)=rj> (
c - Shift of origin
C+WI:jb j ) I:kak u
.
(5)
It is believed that this class is sufficiently broad to adapt the model to slight variations in fuel quality. Note that w (air) has unknown origin, but the origin of u (fuel) can be fixed. This is due to the fact it is difficult to measure the total amount of air entering the combustion, but it can be done so with fuel.
The error signal e(t) is supposed to be a sequence of independent normally distributed random variables with zero mean and common variance V(e). If 0' E [0,1) and if constant input values uo, Wo would be applied sufficiently long period, the output
128
with covariance matrix or its inverse expressed in terms of LDL' product, where L is lower triangular with units on its diagonal and D is a diagonal matrix. Several variants of the algorithm are known and described in literature (Peterka 1986). For some reasons which will be made apparent in section (10) we used the inverse covariance matrix.
4. LINEARIZED MODEL In order to make the use of a conventional identification method (recursive least squares) possible, the model (1) should be linearized. Because both a and 4> are not estimated, the model can be simplified introducing the filtered output Q(t)
Q(t) = 4>-1 (y(t» - a4>-I(y(t-l». I-a
The advantages gained by expressing the inverse covariance matrix as the product are two: (1) the matrix is ensured to be positive semidefinite, (2) the algorithm is insensitive to its singularity.
(6)
This filtered output satisfies
Q(t)
b·w(t -)') = c +" LJj 1 + e(t) = Lk aku(t-k) = Qd(B) + e(t).
6. HEURISTIC WEIGHT-SHAPING (7)
Data obtained by means of a number of experiments showed that the variance V(e) varies with the mean value. This dependence suggests to use the weighted least squares algorithm with the weight function W (y) as follows: V(e) = a~W(y). Here the function WO is considered to be known and a e is an unknown constant, which is estimated.
It is possible to use the derivatives 8Qd/8B to approximate the non-linear function Qd(B). Nevertheless, it has been suggested by intuition and verified by experiments that it is advantageous to use the logarithm q(t) = log Q(t) instead of Q(t), because the logarithm of the CO concentration is a linear function of functions of both fuel and air supply rates
Because the purpose of the model is to predict the air supply rate which would produce an acceptable CO concentration, say YL, on the output for a given fuel supply rate or power demand, it is clear that the model is used mainly locally in the vicinity of YL. This argument suggests to use the following form of W (y)
q(t) = log(c + I)jw(t-j»j
-log(Laku(t-k»+e(t).
(8)
k
This equation may be approximated by
Wl
W(y) q(t)
+
~
~~ I
qd(M(B»+
={ Wl
(B - M(B» + e(t).
Y < Yl Yl ~ Y ~ Y2· >1 Y > Y2
>1
1
(12)
Thus, the model is rather a local approximation of the true process. The constants Wk,Yk are tuning knobs of the algorithm and must be adjusted empirically.
(9)
M(9)
The equation error is still denoted as e(t), though its p.d.f. has changed. .
7. TREATMENT OF SENSOR SATURATION 5. NUMERICALLY RELIABLE ESTIMATION
Bayesian parameter estimation is nothing but a chained conditioning of a prior distribution of probability
By means of linearization, the non-linear model is put to the simple form
y/(t)
= 1/JB + e(t).
Pt+! (B)
(10)
~~ I
M(B) - qd(M(B»
=
M(9)
=
~~ IM(9) B+ e(t).
Po(B) given
(13)
D(t) are the data acquired at the time instant t. With measurements of concentrations, one often encounters a problem due to sensor saturation. Instead of the item of information y(t) = v - output equals a value - one gets known just y(t) > Vsat - output is above a value. Here v is the upper bound of the sensor. Several solutions to this problem are possible. Among them: (a) To process the parameter update as if for y = v. (b) To wait until the saturation vanishes and ignore the data. (c) To perform the duly conditioning Pt(Bly(t + 1) > Vsat)·
That is
q(t) +
= Pt(BID(t + 1»,
(11)
Here 1/J is a vector related to data, and y/ is filtered, mapped and shifted output. The numerically reliable algorithm performing data-conditioning works either
129
let)
The first solution suffers the disadvantage the parameter estimate is biased if the sensor saturation occurs. The second solution can get the estimate trapped after the sensor would have been saturated - and it could be very dangerous. There are also problems with the third solution: the major one is caused by the fact pt(Oly(t + 1) > v) is not Gaussian even though Pt(O) would be so. It is therefore necessary to perform just an approximation of this solution.
uet) wet)
A way to perform the conditioning can be based on Monte-Carlo simulations. It is necessary to generate a sample On fromp(O). Then each sample is not rejected with probability proportional to p(y > vIOn ). Those On which passed through can be used, for example, to calculate the empirical mean and variance to set up Pt+1 (0)
Fig. 2. A model of the fuel quality variations: Ideal fuel is multiplied by a random walk process. Le. there is a large prediction error, the measurement y > Vsat is processed almost as if y = Vsat. If it is not surprising and prediction M(y) also (in the mean) satisfies M(y) > Vsat, the measurement is almost ignored. This observation could be perhaps used as a simpler solution to the problem.
The approximation can be also based on the following obvious idea: A Gaussian p.dJ. in the parameter space is determined by N = n + (n + n 2 )/2 numbers, where n is its dimension. One calculates N functionals (e.g. moments) of pt(Oly(t+ 1) > v) and replaces this p.d.f. with p(m, P), this time a Gaussian one, letting its functional values be the same. We used the following functionals determining the vector of mean values m and the covariance matrix P :
BPI BB
-0
8. FORGETTING To allow for tracking of slowly varying parameter values, it is necessary to perform a transformation of the probability distribution representing the knowledge about their true values. This transformation is called forgetting. There are two features which must be possessed by any forgetting algorithm used to solve the problem discussed in this document. These are treated in the next two sections.
(14)
(J=m -
2
B B02p I
= const. P -1 .
yet)
Boiler
(15)
(J=m
The Bayes formula
9. LINEAR DIRECTIONAL FORGETTING
p(Bly > v) = const.p(B)p(y > vIB).
(16)
The combustion process has two inputs: the fuel and air supply rates (u, w.) It is necessary to forget the parameters in a way reflecting mainly fuel quality variations. An extra parameter 9 can be introduced for the purpose. Then it is supposed that the quality of a fictitious fuel u is always the same, but the gain g(t) varies
is to be used. The term p(y > vlO) is related to the cumulative probability distribution function 4> :
p(y
=1 -
> vlB)
=
u(t) = g(t)u(t), (17)
a e 4>((4>B - v)/a).
g(t+l)
If the following equation is solved numerically by
Vsat IB)
=0
= g(t) + f(t).
(20)
Here f(t) is supposed to be a zero-meaned random variable. The model can be re-written (writing only two of a and b for simplicity):
means of linear approximations (Newton method applied to likelihood equation)
B log p( B)p(y > BB
(19)
(18)
y(t+ 1)
'
= (a1
a2
b1
~ )X
u(t - du)g(t - du) ) u(t-du-1)g(t-du -1) x w(t-dw) ( w(t-dw-1)
then matrix of second order derivatives must be evaluated at each iteration. This can be then directly used as an approximate inverse covariance matrix. The behaviour of the algorithm can be described as follows: If the fact the sensor is saturated is surprising,
+e
(t) .
(21)
The following derivation may seem somewhat confusing, which is due to the fact that the parameter
130
variation model must, rather than describe a real random mechanism, be simple and preserve analytic form of the probability distribution functions. Firstly it is reasonable to suppose the gain g(t) varies slowly compared to u(t). It is therefore possible to write g(t-du) ~ g(t-du -1) and
To allow a slow forgetting also in directions other than fuel quality, the following linear forgetting step can be used
P(t)
= P(t-1) + M(B')M(B)V(f) + -yI.
(32)
Here -yI is a matrix proportional to the unit matrix.
10. REGULARIZED FORGETTING The linear forgetting described in the preceding section performs satisfactorily only if the plant is sufficiently excited by its inputs. In practice this assumption is far from being true and simple linear forgetting must not be used. Equations from the preceding section should be understood only as an auxiliary material explaining an idea of parameter variation.
(22)
It is convenient to introduce time-varying parameters with the gain 9 already absorbed in
B(t) = Bg(t), B(t+ 1) = Bg(t+ 1) = Bg(t) B(t+1)
(23)
+ Bf(t),
= B(t) + Bf(t).
An alternative which proved itself to be better was regularized exponential forgetting (Kulhavy and Kraus 1996). The difference between linear and regularized exponential forgetting will be now explained. The forgetting step is an operator F acting on p( B). With linear forgetting one uses the explicit model of parameter variations to define F. Regularized exponential forgetting does not define F with the help of such a model, but directly as
(24) (25)
To make the calculation simple Bf(t) will be supposed to be a normal random variable. The true parameter value B will be approximated with B(t). Thus the linear forgetting step could be expressed as follows
P(t)
= P(t-1) + B'BV(f).
(26)
The true parameter value must be approximated with M(B) :
P(t)
= P(t-1) + M(B')M(B)V(f).
F(P(B)) =
mjn{)..D(f(B),Po(B))
(27)
= diag(ai,O,O,O),
(33)
The symbols have the following meaning:
It is also necessary to express this matrix M(B')M(B)V(f) as an LDL' product, where L is lower triangular with units on its diagonal and D is a diagonal matrix. This product has the form as follows: D
+ )"D(f(B),p(B))}.
(28)
D(p, q) = Kullback-Leibler information p(B) D(p, q) = p(B) log q(B) dB / Po(B) = worst-case information. The general solution has the form: (34)
(29)
The Kullback-Leibler (Cover and Thomas 1991) divergence will not be discussed here. The optimization problem (33) can be easily solved in the case p, Po are Gaussian probability distribution functions p( ml , PI) and p(mo, Po) respectively. The solution is also Gaussian p.d.f. p(m*, P*) as follows
The above definition of D and L gives sense only for al other than zero. Suppose that al is zero and a2 is not. The matrices will change to
D
= diag(O, a~, 0),
(30)
+ (1 - )..)PI- I = )..PO-Imo + (1 - )..)PI-Iml'
r - I = )..PO-
L= Similarly if al .. .ak
(1°° °° 0) 1
a3/a2
r-Im*
.
(31)
I
(35) (36)
The three parameters ).., mo (vector of mean values) and Po (covariance matrix) are design parameters. It is plain to see that the desired directional properties
1 all equal zero.
131
of the estimation algorithm can be achieved by the following values
Po = V(f)M(8')M(8) + "(1.
CO
(37)
To avoid a possibility of divergence of M(8) due to feedback via Po, a normalized version was suggested V(f)
Po
= 1 + "(2 M (8')M(8) MW)M(8) + "(1.
(38)
CO limit
The LDL' form given by (9) of this covariance matrix remains valid. The mean value 111.0 can be set to an off-line estimate, which is considered to be reliable. Another possibility sets 111.0 = M(8). The former possibility is considered to be safer. The major advantage of the regularized exponential forgetting is that covariance wind-up can occur. Optimal 11. Am/FUEL RATIO SCHEDULING Fig. 3. The optimal air/fuel ratio scheduling visualized.
The model whereof parameters are estimated is used to predict the air supply rate Wo for which the steadystate CO concentration would satisfy Prob(y(t)
> YL)
= E.
The equation (41) is solved via the 'regula falsi' method, which always converges.
(39)
I. e. the limit YL will be exceeded with a small probability E for any fixed t. This probability is independent of t, because yet) is supposed to be a stationary random process. This limit is mapped into qL, which satisfies (provided 4> is increasing)
qL
12. SENSOR FAULT DETECTION The sensors measuring the concentrations of combustion products are sensitive devices which can be subjects to failures. The failures cannot be detected easily, because the functionality of sensor degrades slowly. It is desirable the identification algorithm could detect sensor failures automatically.
= log 4>-1 (Yd,
P(y(t) > YL) = P(q(t) > qL)'
(40)
It is quite important to notice that a model, which is used for the purpose of fault detection may differ from a model for the purpose of prediction. For example, one has to model CO concentration as a system with two outputs in case two sensors are available, though the use of single output model (Y1 + Y2) /2 may be tolerable for the purpose of prediction. Multi sensor model improves sensor fault detection ability automatically. The same applies for other auxiliary items of information, which should be incorporated to the model as outputs. As an example, it is possible to add O2 concentration as an extra output to the model of CO generation. Such model then comprises the knowledge CO concentration is low if O2 concentration is high and vice versa. If the real data contradict this knowledge, it is possible to detect a failure even if the concentration is still within reasonable limits and would be quite possible for some coal types. Thus, the first remark concerning sensor fault detection is that a MIMO model (in case two sensors are available) parameterized with a
The sought air supply rate Wo can be calculated as solution to the following equation:
Prob(q(uo, w~) > qL) Here
E
= E.
(41)
is a value set by the personnel. The equation
(41) may be re-written
M(q(uo,w~)) + nu(Eh!V«uo, wo)) = qL. (42) The two statistics equal approximately
M(q(uo,wo)) d - Oqd - 8
~ qd(UO,W~),
(43)
I8=M(8)'
(44)
V«uo,w~)) ~ (1 +dJPd)V(e).
(45)
Function nu(E) is related to Student p.dJ. It proved to be satisfactory to let the personnel control nu directly, rather than E.
132
13. CONCLUSION
vector of parameters 0 should be used and p.d.f. p(O) calculated.
The ideas presented in this article were implemented to perform the optimal air/fuel ratio scheduling for the Advanced Combustion Control (by Honeywell). The ratio calculated by this module is used as a setpoint for the main controller which manipulates the two air and fuel supply rates. The controller is able to maintain their ratio constant even during a dynamic regime. The ratio must be constant at the combustion chamber. The experience gained during the working tests proved that this mode of control decreases the average air excess and the efficiency increases. The mean CO concentration could be increased, but its variance decreased - it is steadier.
For the purpose of fault detection it is necessary to calculate integral of the following type: 1=
Js
(46)
p(O)d8.
The correct behaviour of the fault detection algorithm is defined by means of a partition of the set of all parameter values and the true (unknown) value O· :
=> normal process, f/. S => faulty process.
O. { E S
(47)
The too simple model (1) loses its capability to predict the CO concentration during a highly dynamic regime. As a result, the estimated error variance increases and so does the excess air (as seen on the figure 3). Thus, during a highly dynamic regime the optimization is ineffective and the excess air has usual values. In contrast, during a steadier regime the excess air can be significantly decreased. It would be desirable to develop a more realistic model of the CO concentration, which would enable further decrease of the air excess.
For the process discussed in this document, elliptical sets S were used. The integral (46) is not easily evaluated even though p(O) is Gaussian. Nevertheless, its value can be approximated evaluating another integral
la. =
J
p(O)exp((O-v)'W(O-v))d8.
(48)
This integral approximates I if W and v are chosen as follows:
S = {x,exp((O - v)'W(O - v))
> 1/2}.
Because of safety considerations, the air supply rate has not been controlled by the CO emissions only in the final product. The optimizer contains also a dynamic model of O2 concentration similar to (1). Then the optimal air supply rate is chosen according to (39), provided the predicted O 2 concentration is not too small. If it is the case, a higher air supply rate is used.
(49)
The value of la. can be calculated analytically. It is convenient to introduce a few auxiliary symbols:
+ W, (P-IM(O) + Wv).
H = p-l h = H-
1
(50) (51)
Then the integral la. equals to 1
exp
(-~((h -
14. REFERENCES
Basseville, M. and I. Nikiforov (1993). Detection of abrupt changes - theory and applications. Prentice Hall. Cover, Thomas M. and Joy A. Thomas (1991). Elements of Information Theory. John Wiley. Hanby, V. I. (1994). Combustion and Pollution Control in Heating Systems. Springer-Verlag. Isermann, Rolf (1984). Process fault detection based on modelling and estimation methods - a survey. Automatica 20,387-404. Kulhavy, R. and F. J. Kraus (1996). On duality of regularized exponential and linear forgetting. Automatica 32, 1403-1415. Peterka, Vaclav (1986). Control of uncertain processes, applied theory and algorithms. The Journal of the Czechoslovak Association for Cybernetics. Supplement to the Journal.
v)'W(h - v)+
VdetPdetH 2 + (h - M(O))'P-l(h - M(O))).
(52)
It is possible to calculate both upper and lower limits as functions of this integral (53) If the exact value of I would be known, it would
be possible to generate alarm signals a(I < T) for a threshold T. Because only a lower limit can be calculated, the alarms must be defined as a(L(Ia.) < T). The way the upper and lower limits can be found will be published elsewhere. How to define the alarm signal thresholds is studied by the theory of statistical decisions, refer to Basseville and Nikiforov (1993).
133
["Q'ip'i'·"i'tN"'Q·:M'Mfflicti.iiiQ'im Ai!
35
-J"'""!"'" I Nt_
n
Fud.. der
lotode(de9<. 1
w_ le",.", Pocf-ol
CP;: ~ ~
-r·~.·.
_.~
{P~.¥ jOO1lD1li: lifihld 11l1200lXXl fuel~
_ien l1OlOOlooo .:;<:ii ]02rolOO
Hfih .....
All """'
10 OOlXXl': low ~
AlS'
101lXWJ
III V
1100
0,4
j40OOlXXl .
tDw_ 112~
.
)., 00lXXl
,. ll\ii....g/l~
t«U.C,
l'1gt"""i2..J~
.J'~~=~_~ ne:':?"M ..[~]_ . . .~.J.. ........•.:..,... :..,.,.
Fig. 4. The optimizer engineer's interfacer: The parameters related to model orders, forgetting and weighting can be changed from here. The ak and bk parameters are visuzlized. The steady-state characteristics is given in the tabular form.
,WtN
TOTGUS08 (R110.61
31 Hc f
99 15 :28 :20
8
~ OPT" ~oneywell 1 .:,. ·.1
{(:.
';.:H
~~~L,
IL:·1GM~
P:dIO
LJ~
;..
:
i:~
, l;
~.
-
L8g~
Ijp.
...' '!
·.'e
(·00:; [·OOG
~ .. 6~
';0(1( •
.3 S7.37
GOOD
~ l>;~;
',. y~
i&Ct~~.
s. g,'(
;.• ~~ ?,Q<;:
GOn(·
: Total AIR Correction ~:.(,\\
:·.tJ~
:'l: 4.0t
\l.~·.. ~e
'\.~;·e'
l.net
~.'.::;~~
O.' .. FEF R~rIC_U'~F
I~:..
~llf.:t
~;.,.,n
l~r..
s::~:.
:)·:~r.
lE:. -; H :
'.\1~
~.e::c
.:. . ~H.~
lH1L;~L:
::~~. 0~
i
l.. . . .
..........1 ..............., ••, , , . . .
---
.......,.
,
,
Lr..~.~.c
J
L..
# ••• _ , - ,
.
Fig. 5. The optimizer operator's interface: The parameters CO-MAX and N-SIGMA are changed from here. The actual values of emissions and input supply rates can be seen.
134
OPTIMIZATION AND CONTROL OF COMBUSTION PROCESS PRODUCTS Daniel Pachner .,1 Vladimir Havlena •• ,1
• Trnka Laboratory of Automatic control, Faculty of Electrical Engineering, Czech Technical University e-mail: [email protected] .. Honeywell Technology Center Prague, Pod vodarenskou vezl 4, 182 08 Prague 8, Czech Republic e-mail: [email protected]
Abstract: A practical adaptive control application is reported. In coal-burned thermal plant, the air supply rate is controlled using a simple non-linear model of carbon monoxide generation. Parameters of the model are on-line estimated using the standard least-squares estimation technique. This model is used for the the optimal air/fuel ratio re-calculation, which should ensure that the steady state CO concentration would exceed a given limit with only a small probability (cautious optimization). Both the limit and probability threshold are the two principal input parameters. To avoid covariance wind-up, the regularized forgetting algorithm is used. To model the fact the parameter variations are caused mainly by the fuel quality fluctuations, the forgetting is directional. The output data are provided by CO concentration sensors, which may get saturated or clogged. Mathematical solutions to both problems is presented. Copyright @2000 IFAC Keywords: power plant; thermal plant; combustion control; cautious optimization techniques; carbon monoxide emission control; failure detection; least-squares estimation; regularized forgetting; sensor saturation; covariance wind-up
1. INTRODUCTION
2. MODELS OF COMBUSTION PRODUCTS GENERATION
When coal is burned, it is desirable to adjust the air supply rate (excess air) in an adaptive manner, so that the changes offuel quality would be compensated for. Adaptive air supply rate allows to stabilize the carbon monoxide concentration bellow a prescribed limit, which also automatically ensures good efficiency of the combustion process (Hanby 1994).
A model relating the air and fuel supply rates to the controlled concentrations is a necessary prerequisite for the control. The structure of the model cannot be however complicated, because its numerical parameters are thought to be estimated by a reliable on-line method. The models developed on the basis of deeper knowledge of the combustion processes and its physics always contradict this requirement. The solution presented is therefore based on a statistical
1 This research was partially supported by the Grant Agency of the Czech republic, grant 102/97/0466, and by the Ministry of Education of the Czech Republic under Project VS97/034
127
~ logO
of the model (1) would fluctuate near the steady-state characteristics expO
(2)
.!.=.!!.
(The symbols M and V denote the mean value and variance respectively.) Provided the function rj> can be approximated by its tangent, the steady state variance can be roughly approximated by
logO %-Q
FIR!
FIR2
V(y)
i
i
wet)
e(t)
u(t)
= V(e)
drj>(x) dx
2 1
I/l-1M(II)
(1 1-
0'~2.
(3)
0'
Thus, the steady state characteristics is a stationary random process and its statistical properties are functions of inputs uo, Wo. Apart from e(t), the other source of uncertainty are the model parameters, which must be also considered to be random variables.
~ yet)
Fig. 1. A model of the CO generation process: its concentration y(t) is a function of filtered the air w(t) to the filtered fuel u(t) supply rates.
The concentration of carbon monoxide exhibits a nonlinearity which can be hardly approximated by its tangent, but should be better approximated by exponential function rj>(x) = exp(cI x + C2)' The constants Cl, C2 must be calculated numerically. To predict mean value of the steady CO concentration, it is better to approximate the output nonlinearity by an exponential and use the well-known formula
model, which captures only statistical characteristics of the dependency instead of modeling it in detail. The remnant difference between the model and the true process is captured by a fictitious random signal. This is the standard way how to tackle overly complicated processes. As it is known that the steady state concentration of any of the combustion products (CO, O2 , NO:t) is a function, say rj>, of air to fuel ratio, the model should comply with this item of information. The following model was suggested:
This way of approximation decreases the bias of CO prediction and should be used rather than (2). (1)
3. PARAMETER ESTIMATION
The symbols have the following meaning:
The model (1) contains a number of unknown parameters: ak,bj,c,O',rj>. It is considered safer not to estimate 0' to avoid unstable (0' > 1) models. The value of 0' is therefore set equal to an off-line estimate. The function rj> can be set equal to a steady state characteristic, which can be obtained experimentally. With the parameters (J = (ak,bj,c) estimated, the model captures the following class of steady state characteristics
y(t) - Product concentration (e. g. CO) u(t) - Fuel supply rate w(t) - Air supply rate
e(t) - Fictitious noise rj> - Steady state characteristics 0' -
Autoregressive coefficient
ak, bj - Moving average coefficients
y(ak,bj,c)=rj> (
c - Shift of origin
C+WI:jb j ) I:kak u
.
(5)
It is believed that this class is sufficiently broad to adapt the model to slight variations in fuel quality. Note that w (air) has unknown origin, but the origin of u (fuel) can be fixed. This is due to the fact it is difficult to measure the total amount of air entering the combustion, but it can be done so with fuel.
The error signal e(t) is supposed to be a sequence of independent normally distributed random variables with zero mean and common variance V(e). If 0' E [0,1) and if constant input values uo, Wo would be applied sufficiently long period, the output
128
with covariance matrix or its inverse expressed in terms of LDL' product, where L is lower triangular with units on its diagonal and D is a diagonal matrix. Several variants of the algorithm are known and described in literature (Peterka 1986). For some reasons which will be made apparent in section (10) we used the inverse covariance matrix.
4. LINEARIZED MODEL In order to make the use of a conventional identification method (recursive least squares) possible, the model (1) should be linearized. Because both a and 4> are not estimated, the model can be simplified introducing the filtered output Q(t)
Q(t) = 4>-1 (y(t» - a4>-I(y(t-l». I-a
The advantages gained by expressing the inverse covariance matrix as the product are two: (1) the matrix is ensured to be positive semidefinite, (2) the algorithm is insensitive to its singularity.
(6)
This filtered output satisfies
Q(t)
b·w(t -)') = c +" LJj 1 + e(t) = Lk aku(t-k) = Qd(B) + e(t).
6. HEURISTIC WEIGHT-SHAPING (7)
Data obtained by means of a number of experiments showed that the variance V(e) varies with the mean value. This dependence suggests to use the weighted least squares algorithm with the weight function W (y) as follows: V(e) = a~W(y). Here the function WO is considered to be known and a e is an unknown constant, which is estimated.
It is possible to use the derivatives 8Qd/8B to approximate the non-linear function Qd(B). Nevertheless, it has been suggested by intuition and verified by experiments that it is advantageous to use the logarithm q(t) = log Q(t) instead of Q(t), because the logarithm of the CO concentration is a linear function of functions of both fuel and air supply rates
Because the purpose of the model is to predict the air supply rate which would produce an acceptable CO concentration, say YL, on the output for a given fuel supply rate or power demand, it is clear that the model is used mainly locally in the vicinity of YL. This argument suggests to use the following form of W (y)
q(t) = log(c + I)jw(t-j»j
-log(Laku(t-k»+e(t).
(8)
k
This equation may be approximated by
Wl
W(y) q(t)
+
~
~~ I
qd(M(B»+
={ Wl
(B - M(B» + e(t).
Y < Yl Yl ~ Y ~ Y2· >1 Y > Y2
>1
1
(12)
Thus, the model is rather a local approximation of the true process. The constants Wk,Yk are tuning knobs of the algorithm and must be adjusted empirically.
(9)
M(9)
The equation error is still denoted as e(t), though its p.d.f. has changed. .
7. TREATMENT OF SENSOR SATURATION 5. NUMERICALLY RELIABLE ESTIMATION
Bayesian parameter estimation is nothing but a chained conditioning of a prior distribution of probability
By means of linearization, the non-linear model is put to the simple form
y/(t)
= 1/JB + e(t).
Pt+! (B)
(10)
~~ I
M(B) - qd(M(B»
=
M(9)
=
~~ IM(9) B+ e(t).
Po(B) given
(13)
D(t) are the data acquired at the time instant t. With measurements of concentrations, one often encounters a problem due to sensor saturation. Instead of the item of information y(t) = v - output equals a value - one gets known just y(t) > Vsat - output is above a value. Here v is the upper bound of the sensor. Several solutions to this problem are possible. Among them: (a) To process the parameter update as if for y = v. (b) To wait until the saturation vanishes and ignore the data. (c) To perform the duly conditioning Pt(Bly(t + 1) > Vsat)·
That is
q(t) +
= Pt(BID(t + 1»,
(11)
Here 1/J is a vector related to data, and y/ is filtered, mapped and shifted output. The numerically reliable algorithm performing data-conditioning works either
129
let)
The first solution suffers the disadvantage the parameter estimate is biased if the sensor saturation occurs. The second solution can get the estimate trapped after the sensor would have been saturated - and it could be very dangerous. There are also problems with the third solution: the major one is caused by the fact pt(Oly(t + 1) > v) is not Gaussian even though Pt(O) would be so. It is therefore necessary to perform just an approximation of this solution.
uet) wet)
A way to perform the conditioning can be based on Monte-Carlo simulations. It is necessary to generate a sample On fromp(O). Then each sample is not rejected with probability proportional to p(y > vIOn ). Those On which passed through can be used, for example, to calculate the empirical mean and variance to set up Pt+1 (0)
Fig. 2. A model of the fuel quality variations: Ideal fuel is multiplied by a random walk process. Le. there is a large prediction error, the measurement y > Vsat is processed almost as if y = Vsat. If it is not surprising and prediction M(y) also (in the mean) satisfies M(y) > Vsat, the measurement is almost ignored. This observation could be perhaps used as a simpler solution to the problem.
The approximation can be also based on the following obvious idea: A Gaussian p.dJ. in the parameter space is determined by N = n + (n + n 2 )/2 numbers, where n is its dimension. One calculates N functionals (e.g. moments) of pt(Oly(t+ 1) > v) and replaces this p.d.f. with p(m, P), this time a Gaussian one, letting its functional values be the same. We used the following functionals determining the vector of mean values m and the covariance matrix P :
BPI BB
-0
8. FORGETTING To allow for tracking of slowly varying parameter values, it is necessary to perform a transformation of the probability distribution representing the knowledge about their true values. This transformation is called forgetting. There are two features which must be possessed by any forgetting algorithm used to solve the problem discussed in this document. These are treated in the next two sections.
(14)
(J=m -
2
B B02p I
= const. P -1 .
yet)
Boiler
(15)
(J=m
The Bayes formula
9. LINEAR DIRECTIONAL FORGETTING
p(Bly > v) = const.p(B)p(y > vIB).
(16)
The combustion process has two inputs: the fuel and air supply rates (u, w.) It is necessary to forget the parameters in a way reflecting mainly fuel quality variations. An extra parameter 9 can be introduced for the purpose. Then it is supposed that the quality of a fictitious fuel u is always the same, but the gain g(t) varies
is to be used. The term p(y > vlO) is related to the cumulative probability distribution function 4> :
p(y
=1 -
> vlB)
=
u(t) = g(t)u(t), (17)
a e 4>((4>B - v)/a).
g(t+l)
If the following equation is solved numerically by
Vsat IB)
=0
= g(t) + f(t).
(20)
Here f(t) is supposed to be a zero-meaned random variable. The model can be re-written (writing only two of a and b for simplicity):
means of linear approximations (Newton method applied to likelihood equation)
B log p( B)p(y > BB
(19)
(18)
y(t+ 1)
'
= (a1
a2
b1
~ )X
u(t - du)g(t - du) ) u(t-du-1)g(t-du -1) x w(t-dw) ( w(t-dw-1)
then matrix of second order derivatives must be evaluated at each iteration. This can be then directly used as an approximate inverse covariance matrix. The behaviour of the algorithm can be described as follows: If the fact the sensor is saturated is surprising,
+e
(t) .
(21)
The following derivation may seem somewhat confusing, which is due to the fact that the parameter
130
variation model must, rather than describe a real random mechanism, be simple and preserve analytic form of the probability distribution functions. Firstly it is reasonable to suppose the gain g(t) varies slowly compared to u(t). It is therefore possible to write g(t-du) ~ g(t-du -1) and
To allow a slow forgetting also in directions other than fuel quality, the following linear forgetting step can be used
P(t)
= P(t-1) + M(B')M(B)V(f) + -yI.
(32)
Here -yI is a matrix proportional to the unit matrix.
10. REGULARIZED FORGETTING The linear forgetting described in the preceding section performs satisfactorily only if the plant is sufficiently excited by its inputs. In practice this assumption is far from being true and simple linear forgetting must not be used. Equations from the preceding section should be understood only as an auxiliary material explaining an idea of parameter variation.
(22)
It is convenient to introduce time-varying parameters with the gain 9 already absorbed in
B(t) = Bg(t), B(t+ 1) = Bg(t+ 1) = Bg(t) B(t+1)
(23)
+ Bf(t),
= B(t) + Bf(t).
An alternative which proved itself to be better was regularized exponential forgetting (Kulhavy and Kraus 1996). The difference between linear and regularized exponential forgetting will be now explained. The forgetting step is an operator F acting on p( B). With linear forgetting one uses the explicit model of parameter variations to define F. Regularized exponential forgetting does not define F with the help of such a model, but directly as
(24) (25)
To make the calculation simple Bf(t) will be supposed to be a normal random variable. The true parameter value B will be approximated with B(t). Thus the linear forgetting step could be expressed as follows
P(t)
= P(t-1) + B'BV(f).
(26)
The true parameter value must be approximated with M(B) :
P(t)
= P(t-1) + M(B')M(B)V(f).
F(P(B)) =
mjn{)..D(f(B),Po(B))
(27)
= diag(ai,O,O,O),
(33)
The symbols have the following meaning:
It is also necessary to express this matrix M(B')M(B)V(f) as an LDL' product, where L is lower triangular with units on its diagonal and D is a diagonal matrix. This product has the form as follows: D
+ )"D(f(B),p(B))}.
(28)
D(p, q) = Kullback-Leibler information p(B) D(p, q) = p(B) log q(B) dB / Po(B) = worst-case information. The general solution has the form: (34)
(29)
The Kullback-Leibler (Cover and Thomas 1991) divergence will not be discussed here. The optimization problem (33) can be easily solved in the case p, Po are Gaussian probability distribution functions p( ml , PI) and p(mo, Po) respectively. The solution is also Gaussian p.d.f. p(m*, P*) as follows
The above definition of D and L gives sense only for al other than zero. Suppose that al is zero and a2 is not. The matrices will change to
D
= diag(O, a~, 0),
(30)
+ (1 - )..)PI- I = )..PO-Imo + (1 - )..)PI-Iml'
r - I = )..PO-
L= Similarly if al .. .ak
(1°° °° 0) 1
a3/a2
r-Im*
.
(31)
I
(35) (36)
The three parameters ).., mo (vector of mean values) and Po (covariance matrix) are design parameters. It is plain to see that the desired directional properties
1 all equal zero.
131
of the estimation algorithm can be achieved by the following values
Po = V(f)M(8')M(8) + "(1.
CO
(37)
To avoid a possibility of divergence of M(8) due to feedback via Po, a normalized version was suggested V(f)
Po
= 1 + "(2 M (8')M(8) MW)M(8) + "(1.
(38)
CO limit
The LDL' form given by (9) of this covariance matrix remains valid. The mean value 111.0 can be set to an off-line estimate, which is considered to be reliable. Another possibility sets 111.0 = M(8). The former possibility is considered to be safer. The major advantage of the regularized exponential forgetting is that covariance wind-up can occur. Optimal 11. Am/FUEL RATIO SCHEDULING Fig. 3. The optimal air/fuel ratio scheduling visualized.
The model whereof parameters are estimated is used to predict the air supply rate Wo for which the steadystate CO concentration would satisfy Prob(y(t)
> YL)
= E.
The equation (41) is solved via the 'regula falsi' method, which always converges.
(39)
I. e. the limit YL will be exceeded with a small probability E for any fixed t. This probability is independent of t, because yet) is supposed to be a stationary random process. This limit is mapped into qL, which satisfies (provided 4> is increasing)
qL
12. SENSOR FAULT DETECTION The sensors measuring the concentrations of combustion products are sensitive devices which can be subjects to failures. The failures cannot be detected easily, because the functionality of sensor degrades slowly. It is desirable the identification algorithm could detect sensor failures automatically.
= log 4>-1 (Yd,
P(y(t) > YL) = P(q(t) > qL)'
(40)
It is quite important to notice that a model, which is used for the purpose of fault detection may differ from a model for the purpose of prediction. For example, one has to model CO concentration as a system with two outputs in case two sensors are available, though the use of single output model (Y1 + Y2) /2 may be tolerable for the purpose of prediction. Multi sensor model improves sensor fault detection ability automatically. The same applies for other auxiliary items of information, which should be incorporated to the model as outputs. As an example, it is possible to add O2 concentration as an extra output to the model of CO generation. Such model then comprises the knowledge CO concentration is low if O2 concentration is high and vice versa. If the real data contradict this knowledge, it is possible to detect a failure even if the concentration is still within reasonable limits and would be quite possible for some coal types. Thus, the first remark concerning sensor fault detection is that a MIMO model (in case two sensors are available) parameterized with a
The sought air supply rate Wo can be calculated as solution to the following equation:
Prob(q(uo, w~) > qL) Here
E
= E.
(41)
is a value set by the personnel. The equation
(41) may be re-written
M(q(uo,w~)) + nu(Eh!V«uo, wo)) = qL. (42) The two statistics equal approximately
M(q(uo,wo)) d - Oqd - 8
~ qd(UO,W~),
(43)
I8=M(8)'
(44)
V«uo,w~)) ~ (1 +dJPd)V(e).
(45)
Function nu(E) is related to Student p.dJ. It proved to be satisfactory to let the personnel control nu directly, rather than E.
132
13. CONCLUSION
vector of parameters 0 should be used and p.d.f. p(O) calculated.
The ideas presented in this article were implemented to perform the optimal air/fuel ratio scheduling for the Advanced Combustion Control (by Honeywell). The ratio calculated by this module is used as a setpoint for the main controller which manipulates the two air and fuel supply rates. The controller is able to maintain their ratio constant even during a dynamic regime. The ratio must be constant at the combustion chamber. The experience gained during the working tests proved that this mode of control decreases the average air excess and the efficiency increases. The mean CO concentration could be increased, but its variance decreased - it is steadier.
For the purpose of fault detection it is necessary to calculate integral of the following type: 1=
Js
(46)
p(O)d8.
The correct behaviour of the fault detection algorithm is defined by means of a partition of the set of all parameter values and the true (unknown) value O· :
=> normal process, f/. S => faulty process.
O. { E S
(47)
The too simple model (1) loses its capability to predict the CO concentration during a highly dynamic regime. As a result, the estimated error variance increases and so does the excess air (as seen on the figure 3). Thus, during a highly dynamic regime the optimization is ineffective and the excess air has usual values. In contrast, during a steadier regime the excess air can be significantly decreased. It would be desirable to develop a more realistic model of the CO concentration, which would enable further decrease of the air excess.
For the process discussed in this document, elliptical sets S were used. The integral (46) is not easily evaluated even though p(O) is Gaussian. Nevertheless, its value can be approximated evaluating another integral
la. =
J
p(O)exp((O-v)'W(O-v))d8.
(48)
This integral approximates I if W and v are chosen as follows:
S = {x,exp((O - v)'W(O - v))
> 1/2}.
Because of safety considerations, the air supply rate has not been controlled by the CO emissions only in the final product. The optimizer contains also a dynamic model of O2 concentration similar to (1). Then the optimal air supply rate is chosen according to (39), provided the predicted O 2 concentration is not too small. If it is the case, a higher air supply rate is used.
(49)
The value of la. can be calculated analytically. It is convenient to introduce a few auxiliary symbols:
+ W, (P-IM(O) + Wv).
H = p-l h = H-
1
(50) (51)
Then the integral la. equals to 1
exp
(-~((h -
14. REFERENCES
Basseville, M. and I. Nikiforov (1993). Detection of abrupt changes - theory and applications. Prentice Hall. Cover, Thomas M. and Joy A. Thomas (1991). Elements of Information Theory. John Wiley. Hanby, V. I. (1994). Combustion and Pollution Control in Heating Systems. Springer-Verlag. Isermann, Rolf (1984). Process fault detection based on modelling and estimation methods - a survey. Automatica 20,387-404. Kulhavy, R. and F. J. Kraus (1996). On duality of regularized exponential and linear forgetting. Automatica 32, 1403-1415. Peterka, Vaclav (1986). Control of uncertain processes, applied theory and algorithms. The Journal of the Czechoslovak Association for Cybernetics. Supplement to the Journal.
v)'W(h - v)+
VdetPdetH 2 + (h - M(O))'P-l(h - M(O))).
(52)
It is possible to calculate both upper and lower limits as functions of this integral (53) If the exact value of I would be known, it would
be possible to generate alarm signals a(I < T) for a threshold T. Because only a lower limit can be calculated, the alarms must be defined as a(L(Ia.) < T). The way the upper and lower limits can be found will be published elsewhere. How to define the alarm signal thresholds is studied by the theory of statistical decisions, refer to Basseville and Nikiforov (1993).
133
["Q'ip'i'·"i'tN"'Q·:M'Mfflicti.iiiQ'im Ai!
35
-J"'""!"'" I Nt_
n
Fud.. der
lotode(de9<. 1
w_ le",.", Pocf-ol
CP;: ~ ~
-r·~.·.
_.~
{P~.¥ jOO1lD1li: lifihld 11l1200lXXl fuel~
_ien l1OlOOlooo .:;<:ii ]02rolOO
Hfih .....
All """'
10 OOlXXl': low ~
AlS'
101lXWJ
III V
1100
0,4
j40OOlXXl .
tDw_ 112~
.
)., 00lXXl
,. ll\ii....g/l~
t«U.C,
l'1gt"""i2..J~
.J'~~=~_~ ne:':?"M ..[~]_ . . .~.J.. ........•.:..,... :..,.,.
Fig. 4. The optimizer engineer's interfacer: The parameters related to model orders, forgetting and weighting can be changed from here. The ak and bk parameters are visuzlized. The steady-state characteristics is given in the tabular form.
,WtN
TOTGUS08 (R110.61
31 Hc f
99 15 :28 :20
8
~ OPT" ~oneywell 1 .:,. ·.1
{(:.
';.:H
~~~L,
IL:·1GM~
P:dIO
LJ~
;..
:
i:~
, l;
~.
-
L8g~
Ijp.
...' '!
·.'e
(·00:; [·OOG
~ .. 6~
';0(1( •
.3 S7.37
GOOD
~ l>;~;
',. y~
i&Ct~~.
s. g,'(
;.• ~~ ?,Q<;:
GOn(·
: Total AIR Correction ~:.(,\\
:·.tJ~
:'l: 4.0t
\l.~·.. ~e
'\.~;·e'
l.net
~.'.::;~~
O.' .. FEF R~rIC_U'~F
I~:..
~llf.:t
~;.,.,n
l~r..
s::~:.
:)·:~r.
lE:. -; H :
'.\1~
~.e::c
.:. . ~H.~
lH1L;~L:
::~~. 0~
i
l.. . . .
..........1 ..............., ••, , , . . .
---
.......,.
,
,
Lr..~.~.c
J
L..
# ••• _ , - ,
.
Fig. 5. The optimizer operator's interface: The parameters CO-MAX and N-SIGMA are changed from here. The actual values of emissions and input supply rates can be seen.
134