- Email: [email protected]
Black-Box Modeling of an Industrial Combustion Engine Plant
Copyright © IFAC System Identification, Kitakyushu , Fukuoka, Japan, 1997
BLACK-BOX MODELING OF AN INDUSTRIAL COMBUSTION ENGINE PLANT
Reinhard Korb, Peter Skorjanz, H. Peter Jorgl
Institute for Machine and Process Automation Vienna University afTechnology, Vienna, Austria e-mail.·[email protected]
Abstract: The linear black-box modeling of a combustion engine is described . As the plant is unstable and the variance of the disturbing noise in the system output is high, standard LS-identification yields insufficient results. Various other well established identification methods such as extended least squares, least squares with data filtering, least squares with order reduction in the frequency domain, weighting sequence estimation, instrumental variable method and correlation analysis with least squares estimation are discussed for a two input single output system. The applicability of these methods for identification of combustion engines is discussed and model validation plots are presented. Keywords: Linear estimation, Internal combustion engines, Identification, MISO
Furthermore it is possible to burn gases, which cannot be utilized by other processes.
I. INTRODUCTION In black-box modeling, a dynamic description of the input-output mapping of a system, using no inner structure information, but measured in- and output data is obtained. This is very useful for plants, for which physical modeling is too time consuming and not effective because of unknown structure and parameters. Plenty of literature is available and surveyed in Isermann (1992a and I 992b). We examined data sets measured at a combustion engine plant. This plant has two characteristics which make identification difficult: Firstly the plant is unstable and secondly the measured signals are corrupted by strong. unknown noise. Only few suggestions for this case are available in literature.
' t
_ift._
I " ,' - -: synchronous i , . --" generator i . :
gas
,
~
'"
:
't . ,__ .. ___._.__ . _ ... __ ..1 engine speed sensor
cooling water IllIXer , A
air
2. DESCRIPTION OF THE PLANT
Fig. I: Flow diagram of combustion engine plant. The numbers indicate measurement locations
The combustion engines investigated in this work, are used for the cogeneration of heat and electrical power. Therefore. these engines are distinguished by a high utilization of primary energy of almost 90 O/C .
The stand alone operation mode has been investigated , which is the more challenging one,
191
The difference
because in this case the frequency must be kept constant by a speed controller. Figure I shows a simple flow diagram of the plant.
(4 )
£(k)=y(k)-,}'(k)
is called the output error. Unfortunately this error is not linear in the unknown parameters. Furthermore in the case of unstable plants this error definition is not useful , as y(k) and y(k) diverge.
As the speed control loop is to be modeled, the control signal is the throttle position, the controlled signal is the engine speed. The electrical power is the disturbance variable . The gas mixer is driven by a separate controller, which guarantees low NO, emissions in the exhaust gas. Figure 2 shows the inand outputs of the two input, single output system .
Instead of eqn . (I) equ. (5) can be used : A(q -' ) y(k)= B(q -' )q -d" u(k)+
(5)
+ B, (q -' )q -d, z(k )+b+ A(q -' )n( k). ath throttle
engine generator
Notice, that A(q· I), B(q·l) and B:(q·l) are different from those in equation (I), as A(q· l) in eqn. (5) is the common denominator of A(q·l) and A:(q· l) in eqn . (I). The model output eqn . (3) is substituted by
engme speed n rnOI
Fig. 2: Black box model of the plant
A(q-I ).v(k)= =B(q - I )q -J" u(k)+B: (q-1 )q - d: z(k)+b
Everyone is familiar with the stable behavior of combustion engines in cars. They are stable, because the load increases strongly as the engine speed increases. In the case of the investigated combustion engines the situation is completely different: The load is independent of the engine speed. Therefore the plant is unstable and measurements are not possible in open loop.
(6)
and for a given set of estimates A, B, B: the so called equation error is defined by e(k)
=A(q -' ) y(k)- B(q -' )q -J" u(k)-
(7)
_B, (q -')q -d, z(k)+b.
3.2 Standard least squares (LS) identification
3. IDENTIFICA TlON OF LINEAR MODELS The 3.1 Basic definitions
LS-Algorithm
determines
the
polynomials
A(q· I) , B(q· l) and Biq· l) by minimizing
(J(A'B'B»=(I~(e(k» ' ) =>A,8.B,
It is assumed , that the inputs are related to the output of the system by B(q-I) _ B _(q - I) _ y(k)=---q d" u(k)+' q d: z(k)+b+n(k) . A(q-I) A: (q - I)
Using matrix abbreviations
calculations
and
the
(8) following
(I)
A more detailed representation is
(2)
where y(k) is the output, u(k) the control variable, z(k) the di sturbance signal and n(k) some additional noise. q.1 is the backward shift operator. d" respectively d: are time delays. 11 and n: are the system orders. The mean value of u(k) , y(k) and z(k) need not be zero as the term b is introduced . The identification algorithms computes estimates h, hi ,,,,, hn , al'" '' all' h:I , . .. , h:II _ , a:I....' aZII , for the parameters. The model output is given by
y=[y(n+d) . . . y(N)]T
(10)
e=[e(n+d) .. . e(N)]T
(I I)
x,(k)=[-y(k-I) ... -y(k-n)]
(12)
x"Ck)=[u(k-dll-l) . . . uCk-d,,-Il)]
(13)
x:(k)=[z(k-d: -I) . .. z(k-d: -n)]
(14)
x(k)=[1 x,(k) x,,(k) x:(k)]
( 15)
x(n+d)
X= X(Il+;d+l) [
1
(16)
x(N) (17)
the well known LS estimate equation is obtained : 8=(X T X)-IX T y .
192
(18)
This equation is solvable, if the system is persistently excited. Furthermore the parameter estimates are unbiased (£(9)=9), if I . the order of the controller is higher or equal than the order of the model or the reference signal yJk) is persistently exciting 2. the noise signal can be generated by IIA-filtering of white noise. This means that white noise can be generated by FIR-filtering of n(k) with the filter polynomial A and filter order n . 3. the disturbing variable z(t) changes its value only at sampling instants. In this case
(21 ) is estimated by the equations obtained from the real part and the imaginary part of eqn. (20) CRe(jw)+ +( CRe (jw )COS( W)+ Clm (j w )sin( W) )a l + ... + +( CRe (jw)cos(n w)+C lm (jw)sin(n w) )an =
(22)
=~ cos( w)+b2 cos(2w)+ .. .+bn cos(n w) Clm(jw)+ + (C lm (jw )cos( W)-C Re (jw )sin( w) )a l + .. .+ +( Clm (jw )cos(n W)-C Re (jW )sin(n W) )an =
(23)
-bl sine w)-b2 sin(2w)-.. .-bn sin(n w)
holds .
From the real and imaginary part of equations (22) and (23) for m (m>n) frequencies w, a matrix equation for the unknown parameters of the low order approximation is derived . It has to be mentioned, that the upper bound for the chosen w has to be designed carefully.
ad I. During the experiments the reference signal was not kept constant to fulfill this condition. ad 2. This is not fulfilled in our case. The bias of the parameters is proportional to the variance of the noise. Due to bad signal to noise (SIR) ratio, the simulation results are not satisfying. ad 3. This requirement is not fulfilled. If the sampling time is small enough, the difference between the left and the right side of eqn . (19) is also small.
3.5 Standard least squares (LS) idelltification with data filtering An considerable part of the noise is caused by the torsional oscillation of the crank shaft. Therefore the in- and output data sets were filtered forwards and backwards by Butterworth filter of lOth order. The cut-off frequency was chosen to be 1/20 of the sampling frequency - 1.67 Hz. The parameters of the filter must be chosen very carefully.
The simulation results using the standard LS-method are not satisfying.
3.3 Standard least squares (LS) identification with high model orders Using very high model orders (e.g. 40) yields good results. Obviously the plant order is not so high, especially not in the discrete time domain. The transfer functions between u(k) respectively z(k) and y(k) are overdetermined by models of such high order. But the higher order also leads to a longer FIRFilter between n(k) and the white noise v(k).
3.6 Weighting sequence estimation
Consider the convolution sum:
y(k)= Lg(l) u(k-I-dl/)+ /=0
(24)
~
3.4 Standard least squares (LS) identification with high model orders and order reduction in the frequency domain
+ Lg:(I) z(k-l-dJ+b+n(k) /=0
This equation can be utilized for identification, if there is a bound p, such that
In order to reduce the high order a model reduction in the frequency domain has been used. The frequency response function of the discrete time model is given by C( '
)W
) C (.
=
) 'C (. ) B(Z-') ",. )W +) I", )W = A(Z -') z=e )WT
Ig(l)I<£ V/>p
system
(25)
for all arbitrary £ holds. The advantage of identifying weighting sequences is that no noise filter specifications are demanded.
(20)
O:5:w:5:rr
Equation (25) does not hold for unstable systems. Kouvaritakis and Rossiter (1992) suggested to replace the convolution sum of the causal weighting sequence by the convolution sum of the so called bicausal weighting sequence, which describes a
C(jw) can be determined by the high order model. A new low order model
193
As E(e(k» is zero, the following cross correlation functions are also zero :
unstable system by a non causal weighting sequence
g h(l), I E )-00,00[ .
(0)
y(k)= 'Ig ,,(l)u(k-l-d.)+ (26)
/:--
(11 ) f:-oo
Using the abbreviations
Kouvaritakis and Rossiter (1992) state that even for unstable systems a bound p can be found that eqn. (25) holds. In this application, no such a bound p could be found, probably because of the noise. Therefore no further investigations were made into this field .
w(k +d -d ll + I) w(k)= : [ w(k -d + N +tmin - t max ) ] ll z (k+d-dll+l) z(k)= : ] [ z(k -d + N +tmin - t max ) ll
T
(2)
T
(1)
3.7 Extellded least squares (ELS) idelltification (4)
The extended least squares (ELS) - method gives unbiased parameter estimates for noise filters of the type D( lI(k)
N= N +tmin -t max + I equation (0) and (11) can be written notation for all t from t/llin to t"IUX
- I
q}v(k) A(q )
(27)
05) In
matrix
where v(k) is white noise . In the recursive version of the ELS algorithm it can be seen that the parameters, especially those of D(q· l) do not converge within the length (N: 1000-4000 points ) of the available data sets. Therefore we decided, that ELS is not the right tool to solve our problem .
3.8 Instrumental Variable (IV) method This method is suggested in literature (e.g. Isermann (1992a» for the applications, which do not fulfill special noise filter conditions. However, in this application the iterative algorithm did not converge because of the unstable system identification in closed loop and correlation of the auxilary model and the equation error.
if the assumptions 1/11I1X ~ -d and 1",in ~ 0 hold . E( . ) is the symbol for the expected value operator. The suitable choice of 1/11in and 1/11I1X is described at the end of this chapter. Otherwise, the dimensions of the matrices must be corrected because of the finite length of the data vectors. Using the equation (38)
y=[y(1 max +d) .. . y(N+1min)]T
(40)
with
3.9 Minimalisation of the cross correlation function between the input and the equation error (eORLS) Another approach to parameter estimation is suggested by Isermann (I 992b ): The model parameters are chosen in the best way, if and only if the equation error does not include components, which can be traced back to the system inputs . In the case of close loop identification the cross covariance function of the external exciting signal w(k) (e.g. ylk» and e(k) resp. of z(k) and e(k) must be zero for all t :
and the definitions of x and e in the equations (15) and (17), equations (6) and (37) yield:
(28)
Rj t)=O 'lit
e=y-xe
(29)
194
It is necessary to use more rows than parameters in eqn . (50) in order to decrease the influence of defective estimates of the first step. The bounds "tmin and "t max are chosen such that
From the definition of the correlation function, we get for equation (41) =[<1>"'\ <1> ...\. .
('t m:in +d,,) : :. <1>...\ ('t min:+ d " -n) .
.
.
1
;/"t)""O \i('t"t max)
Notice, that this methods yields unbiased parameter estimates for all stationary noise signals. No special noise filter is required . Using an additional closed loop input, this method is also very well suited for identification of unstable systems.
(43)
<1>".,. ('t max +d,,) . .. <1>...,. ('t max +d" -n)
...// ('t min ) ...w//('tmin-n)] :
. ..
(51 )
i,jE{U,w,y,z}
(44)
:
[ <1>".// ( 't max) .. . \\"// ( 't max -n)
4. SIMULATION RESULTS As the plant is unstable, validation of the model is possible in closed loop only. The measured signals .vJ.k) (desired engine speed) and PetCk) (electrical power output) are used in simulation as inputs and the simulated signals atl,(k) (throttle position) and nff/m(k) (engine speed) are compared to the measured signals CJ.'h(k) and nnlll,(k). The first 30 sec of the data sets are used for training the model , the rest is used for validation.
(46)
Using similar abhreviations we get for equation (42)
Pel
(47)
Yr+
Equs. (46) and (47) yield a set of 2('tmllx-'tm;n+l) linear equations for the unknown parameter vector 8. Additionally the stationary hias-equation should be fulfilled :
In order to compare the applicability of the different algorithms for the combustion engine plant, simulation studies were carried out. Four typical simulation results are shown in figure 4 - 7. For each identification the same data set is used . The system order is assumed to be 4. It is demanded, that the model does not reproduce the torsional oscillation of the crank shaft.
or with lll xm being a (n X m)-matrix of ones (49)
Now equation (49), (46) and (47) are written in matrix notation: E(y)
I
$" . \ (, min )
E(w)
:
$..., (, ma, )
E(w)
$" . ('m ;n)
E( z )
:
$" (, ma, )
"nmot
Fig 3: Simulation in the closed loop
ECv)=b-a l E(y)- .. .-aIlE( y )+
E (y )=[I\-E(y) · ll xlliE( y) ll xn iE( y )llxlIle
model of the " controller CJ.th
model of engine and generator
-E(y) · l lxn E(u) · l lxn E( z) · l lxn
1520
-G>wy
- $:,
w"
$,,,
'"i:
$w:
'E 1500
--
e
~ o
1480
0-
$::
Cl)
~ 1460 'Oi;
E(z)
(50)
51440 '
The following application of these equations is suggested in order to avoid large matrices: I . Estimates of the expected values and the correlation functions are computed. 2. The matrix eqn. (50) is build. 3. The parameter vector 8 is computed using the Least Squares-solution of eqn . (50).
... . ... · desired - - - - - - measured ---simulated
1420 L~~~~_~=======:.J 34 36 38 40 42 44 46 48 50 time [s1 Fig. 4 : Measured and simulated signals of the engine speed. The 4th order model was generated using the standard LS-method .
195
In figure 4 it can be seen that the simple least squares model reproduces the plant's behavior roughly, but at the peaks the model's output is poor. Figure 5 shows, that filtering of the data vectors improves the simulation result considerably. This means, that the part of noise, which cannot be represented by the I lA-filter, is located mainly in the high frequency domain. Be aware of the fact, that the filter order, the cut off frequency and the system order have to be chosen carefully. Therefore some trials are necessary for getting good results with a new type of plant.
1520
c
] 1500
13 1480
...Co.
~ 1460
t:
'on
~ 1440 1420 34
desired - - - - - - mt:asured - - - simulated 36
38
40
42
44
46
48
50
time [5] Fig. 7: Measured and Simulated signals of the engine speed. The 4lh order model was generated using the minimalization of the cross correlation function between the input and the equation error
1520 t:
] 1500
Finally. figure 7 shows that the cross correlation minimalization method yields simulation results, which are comparable to those shown in figure 5 and 6. The great advantage of this method is that the only parameters to be chosen (tmin and t max) are insensitive to changes. Their choice over a wide range influences the simulation result negligibly .
13
...Co. 14RO
~ 1460
t: '8)
i 3 1440 1420 L~~~~_~=======.J 34 36 40 42 44 46 48 50
time [s] Fig 5: Measured and simulated signals of the engine speed. The 4lh order model was generated using the standard LS-method with filtering the data vectors.
5. CONCLUSION
In this paper different linear identification methods are compared to the standard least squares method for their applicability in modeling an industrial combustion engine plant. As the plant is unstable, and there are many unmodeled dynamics and noise in the data sets, numerous methods suggested in literature fail. However, three algorithms are presented, which yield a drastic improvement compared to the result obtained by the standard least squares method .
1520
c
.- 1500 ~ :; 1480
......
0-
~
1460
t:
'on ~
. ... . desired - - - - - - measured - - - simulated
1440
142~4
36
38
40
42
44
46
48
ACKNOWLEDGMENTS 50
Many thanks to lenbacher Energiesysteme AG , who made the data sets available to us. The project was supported by the Austrian Forschungs-ForderungsFonds der gewerblichen Wirtschaft (FFF).
time [s1 Fig. 6: Measured and simulated signals of the engine speed. The 4lh order model was generated by frequency domain order reduction of a 40 lh order LS-model
REFERENCES Figure 6 visualizes, that order reduction in the frequency domain also yields good approximation. Therefore, high order standard least squares models are able to reproduce the system behavior even if the noise filter specification is violated . Furthermore, it is shown that the order can be reduced in the frequency domain without deteriorating the simulation result. As in the case of data filtering, careful tuning of some parameters is necessary .
Isermann, R. (1992a). Identifikation dynamischer Systeme . Vol. I, 2 nd print, Springer, Berlin . Isermann, R. (1992b). Identifikation dynamischer S.vsteme . Vol. 2, 2 nd print, Springer, Berlin. Kay S.M. (1988) . Modern spectral estimation : theory and application. Prentice-Hall, London . Kouvaritakis, B.; Rossiter, I.A . (1992) . Use of bicausal weighting sequences in least squares identification of open loop unstable dynamic systems. lEE Proceedings 139, No. 3, pp328-336
196
BLACK-BOX MODELING OF AN INDUSTRIAL COMBUSTION ENGINE PLANT
Reinhard Korb, Peter Skorjanz, H. Peter Jorgl
Institute for Machine and Process Automation Vienna University afTechnology, Vienna, Austria e-mail.·[email protected]
Abstract: The linear black-box modeling of a combustion engine is described . As the plant is unstable and the variance of the disturbing noise in the system output is high, standard LS-identification yields insufficient results. Various other well established identification methods such as extended least squares, least squares with data filtering, least squares with order reduction in the frequency domain, weighting sequence estimation, instrumental variable method and correlation analysis with least squares estimation are discussed for a two input single output system. The applicability of these methods for identification of combustion engines is discussed and model validation plots are presented. Keywords: Linear estimation, Internal combustion engines, Identification, MISO
Furthermore it is possible to burn gases, which cannot be utilized by other processes.
I. INTRODUCTION In black-box modeling, a dynamic description of the input-output mapping of a system, using no inner structure information, but measured in- and output data is obtained. This is very useful for plants, for which physical modeling is too time consuming and not effective because of unknown structure and parameters. Plenty of literature is available and surveyed in Isermann (1992a and I 992b). We examined data sets measured at a combustion engine plant. This plant has two characteristics which make identification difficult: Firstly the plant is unstable and secondly the measured signals are corrupted by strong. unknown noise. Only few suggestions for this case are available in literature.
' t
_ift._
I " ,' - -: synchronous i , . --" generator i . :
gas
,
~
'"
:
't . ,__ .. ___._.__ . _ ... __ ..1 engine speed sensor
cooling water IllIXer , A
air
2. DESCRIPTION OF THE PLANT
Fig. I: Flow diagram of combustion engine plant. The numbers indicate measurement locations
The combustion engines investigated in this work, are used for the cogeneration of heat and electrical power. Therefore. these engines are distinguished by a high utilization of primary energy of almost 90 O/C .
The stand alone operation mode has been investigated , which is the more challenging one,
191
The difference
because in this case the frequency must be kept constant by a speed controller. Figure I shows a simple flow diagram of the plant.
(4 )
£(k)=y(k)-,}'(k)
is called the output error. Unfortunately this error is not linear in the unknown parameters. Furthermore in the case of unstable plants this error definition is not useful , as y(k) and y(k) diverge.
As the speed control loop is to be modeled, the control signal is the throttle position, the controlled signal is the engine speed. The electrical power is the disturbance variable . The gas mixer is driven by a separate controller, which guarantees low NO, emissions in the exhaust gas. Figure 2 shows the inand outputs of the two input, single output system .
Instead of eqn . (I) equ. (5) can be used : A(q -' ) y(k)= B(q -' )q -d" u(k)+
(5)
+ B, (q -' )q -d, z(k )+b+ A(q -' )n( k). ath throttle
engine generator
Notice, that A(q· I), B(q·l) and B:(q·l) are different from those in equation (I), as A(q· l) in eqn. (5) is the common denominator of A(q·l) and A:(q· l) in eqn . (I). The model output eqn . (3) is substituted by
engme speed n rnOI
Fig. 2: Black box model of the plant
A(q-I ).v(k)= =B(q - I )q -J" u(k)+B: (q-1 )q - d: z(k)+b
Everyone is familiar with the stable behavior of combustion engines in cars. They are stable, because the load increases strongly as the engine speed increases. In the case of the investigated combustion engines the situation is completely different: The load is independent of the engine speed. Therefore the plant is unstable and measurements are not possible in open loop.
(6)
and for a given set of estimates A, B, B: the so called equation error is defined by e(k)
=A(q -' ) y(k)- B(q -' )q -J" u(k)-
(7)
_B, (q -')q -d, z(k)+b.
3.2 Standard least squares (LS) identification
3. IDENTIFICA TlON OF LINEAR MODELS The 3.1 Basic definitions
LS-Algorithm
determines
the
polynomials
A(q· I) , B(q· l) and Biq· l) by minimizing
(J(A'B'B»=(I~(e(k» ' ) =>A,8.B,
It is assumed , that the inputs are related to the output of the system by B(q-I) _ B _(q - I) _ y(k)=---q d" u(k)+' q d: z(k)+b+n(k) . A(q-I) A: (q - I)
Using matrix abbreviations
calculations
and
the
(8) following
(I)
A more detailed representation is
(2)
where y(k) is the output, u(k) the control variable, z(k) the di sturbance signal and n(k) some additional noise. q.1 is the backward shift operator. d" respectively d: are time delays. 11 and n: are the system orders. The mean value of u(k) , y(k) and z(k) need not be zero as the term b is introduced . The identification algorithms computes estimates h, hi ,,,,, hn , al'" '' all' h:I , . .. , h:II _ , a:I....' aZII , for the parameters. The model output is given by
y=[y(n+d) . . . y(N)]T
(10)
e=[e(n+d) .. . e(N)]T
(I I)
x,(k)=[-y(k-I) ... -y(k-n)]
(12)
x"Ck)=[u(k-dll-l) . . . uCk-d,,-Il)]
(13)
x:(k)=[z(k-d: -I) . .. z(k-d: -n)]
(14)
x(k)=[1 x,(k) x,,(k) x:(k)]
( 15)
x(n+d)
X= X(Il+;d+l) [
1
(16)
x(N) (17)
the well known LS estimate equation is obtained : 8=(X T X)-IX T y .
192
(18)
This equation is solvable, if the system is persistently excited. Furthermore the parameter estimates are unbiased (£(9)=9), if I . the order of the controller is higher or equal than the order of the model or the reference signal yJk) is persistently exciting 2. the noise signal can be generated by IIA-filtering of white noise. This means that white noise can be generated by FIR-filtering of n(k) with the filter polynomial A and filter order n . 3. the disturbing variable z(t) changes its value only at sampling instants. In this case
(21 ) is estimated by the equations obtained from the real part and the imaginary part of eqn. (20) CRe(jw)+ +( CRe (jw )COS( W)+ Clm (j w )sin( W) )a l + ... + +( CRe (jw)cos(n w)+C lm (jw)sin(n w) )an =
(22)
=~ cos( w)+b2 cos(2w)+ .. .+bn cos(n w) Clm(jw)+ + (C lm (jw )cos( W)-C Re (jw )sin( w) )a l + .. .+ +( Clm (jw )cos(n W)-C Re (jW )sin(n W) )an =
(23)
-bl sine w)-b2 sin(2w)-.. .-bn sin(n w)
holds .
From the real and imaginary part of equations (22) and (23) for m (m>n) frequencies w, a matrix equation for the unknown parameters of the low order approximation is derived . It has to be mentioned, that the upper bound for the chosen w has to be designed carefully.
ad I. During the experiments the reference signal was not kept constant to fulfill this condition. ad 2. This is not fulfilled in our case. The bias of the parameters is proportional to the variance of the noise. Due to bad signal to noise (SIR) ratio, the simulation results are not satisfying. ad 3. This requirement is not fulfilled. If the sampling time is small enough, the difference between the left and the right side of eqn . (19) is also small.
3.5 Standard least squares (LS) idelltification with data filtering An considerable part of the noise is caused by the torsional oscillation of the crank shaft. Therefore the in- and output data sets were filtered forwards and backwards by Butterworth filter of lOth order. The cut-off frequency was chosen to be 1/20 of the sampling frequency - 1.67 Hz. The parameters of the filter must be chosen very carefully.
The simulation results using the standard LS-method are not satisfying.
3.3 Standard least squares (LS) identification with high model orders Using very high model orders (e.g. 40) yields good results. Obviously the plant order is not so high, especially not in the discrete time domain. The transfer functions between u(k) respectively z(k) and y(k) are overdetermined by models of such high order. But the higher order also leads to a longer FIRFilter between n(k) and the white noise v(k).
3.6 Weighting sequence estimation
Consider the convolution sum:
y(k)= Lg(l) u(k-I-dl/)+ /=0
(24)
~
3.4 Standard least squares (LS) identification with high model orders and order reduction in the frequency domain
+ Lg:(I) z(k-l-dJ+b+n(k) /=0
This equation can be utilized for identification, if there is a bound p, such that
In order to reduce the high order a model reduction in the frequency domain has been used. The frequency response function of the discrete time model is given by C( '
)W
) C (.
=
) 'C (. ) B(Z-') ",. )W +) I", )W = A(Z -') z=e )WT
Ig(l)I<£ V/>p
system
(25)
for all arbitrary £ holds. The advantage of identifying weighting sequences is that no noise filter specifications are demanded.
(20)
O:5:w:5:rr
Equation (25) does not hold for unstable systems. Kouvaritakis and Rossiter (1992) suggested to replace the convolution sum of the causal weighting sequence by the convolution sum of the so called bicausal weighting sequence, which describes a
C(jw) can be determined by the high order model. A new low order model
193
As E(e(k» is zero, the following cross correlation functions are also zero :
unstable system by a non causal weighting sequence
g h(l), I E )-00,00[ .
(0)
y(k)= 'Ig ,,(l)u(k-l-d.)+ (26)
/:--
(11 ) f:-oo
Using the abbreviations
Kouvaritakis and Rossiter (1992) state that even for unstable systems a bound p can be found that eqn. (25) holds. In this application, no such a bound p could be found, probably because of the noise. Therefore no further investigations were made into this field .
w(k +d -d ll + I) w(k)= : [ w(k -d + N +tmin - t max ) ] ll z (k+d-dll+l) z(k)= : ] [ z(k -d + N +tmin - t max ) ll
T
(2)
T
(1)
3.7 Extellded least squares (ELS) idelltification (4)
The extended least squares (ELS) - method gives unbiased parameter estimates for noise filters of the type D( lI(k)
N= N +tmin -t max + I equation (0) and (11) can be written notation for all t from t/llin to t"IUX
- I
q}v(k) A(q )
(27)
05) In
matrix
where v(k) is white noise . In the recursive version of the ELS algorithm it can be seen that the parameters, especially those of D(q· l) do not converge within the length (N: 1000-4000 points ) of the available data sets. Therefore we decided, that ELS is not the right tool to solve our problem .
3.8 Instrumental Variable (IV) method This method is suggested in literature (e.g. Isermann (1992a» for the applications, which do not fulfill special noise filter conditions. However, in this application the iterative algorithm did not converge because of the unstable system identification in closed loop and correlation of the auxilary model and the equation error.
if the assumptions 1/11I1X ~ -d and 1",in ~ 0 hold . E( . ) is the symbol for the expected value operator. The suitable choice of 1/11in and 1/11I1X is described at the end of this chapter. Otherwise, the dimensions of the matrices must be corrected because of the finite length of the data vectors. Using the equation (38)
y=[y(1 max +d) .. . y(N+1min)]T
(40)
with
3.9 Minimalisation of the cross correlation function between the input and the equation error (eORLS) Another approach to parameter estimation is suggested by Isermann (I 992b ): The model parameters are chosen in the best way, if and only if the equation error does not include components, which can be traced back to the system inputs . In the case of close loop identification the cross covariance function of the external exciting signal w(k) (e.g. ylk» and e(k) resp. of z(k) and e(k) must be zero for all t :
and the definitions of x and e in the equations (15) and (17), equations (6) and (37) yield:
(28)
Rj t)=O 'lit
e=y-xe
(29)
194
It is necessary to use more rows than parameters in eqn . (50) in order to decrease the influence of defective estimates of the first step. The bounds "tmin and "t max are chosen such that
From the definition of the correlation function, we get for equation (41) =[<1>"'\ <1> ...\. .
('t m:in +d,,) : :. <1>...\ ('t min:+ d " -n) .
.
.
1
;/"t)""O \i('t
Notice, that this methods yields unbiased parameter estimates for all stationary noise signals. No special noise filter is required . Using an additional closed loop input, this method is also very well suited for identification of unstable systems.
(43)
<1>".,. ('t max +d,,) . .. <1>...,. ('t max +d" -n)
...// ('t min ) ...
. ..
(51 )
i,jE{U,w,y,z}
(44)
:
[ <1>".// ( 't max) .. . \\"// ( 't max -n)
4. SIMULATION RESULTS As the plant is unstable, validation of the model is possible in closed loop only. The measured signals .vJ.k) (desired engine speed) and PetCk) (electrical power output) are used in simulation as inputs and the simulated signals atl,(k) (throttle position) and nff/m(k) (engine speed) are compared to the measured signals CJ.'h(k) and nnlll,(k). The first 30 sec of the data sets are used for training the model , the rest is used for validation.
(46)
Using similar abhreviations we get for equation (42)
Pel
(47)
Yr+
Equs. (46) and (47) yield a set of 2('tmllx-'tm;n+l) linear equations for the unknown parameter vector 8. Additionally the stationary hias-equation should be fulfilled :
In order to compare the applicability of the different algorithms for the combustion engine plant, simulation studies were carried out. Four typical simulation results are shown in figure 4 - 7. For each identification the same data set is used . The system order is assumed to be 4. It is demanded, that the model does not reproduce the torsional oscillation of the crank shaft.
or with lll xm being a (n X m)-matrix of ones (49)
Now equation (49), (46) and (47) are written in matrix notation: E(y)
I
$" . \ (, min )
E(w)
:
$..., (, ma, )
E(w)
$" . ('m ;n)
E( z )
:
$" (, ma, )
"nmot
Fig 3: Simulation in the closed loop
ECv)=b-a l E(y)- .. .-aIlE( y )+
E (y )=[I\-E(y) · ll xlliE( y) ll xn iE( y )llxlIle
model of the " controller CJ.th
model of engine and generator
-E(y) · l lxn E(u) · l lxn E( z) · l lxn
1520
-G>wy
- $:,
$,,,
'"i:
$w:
'E 1500
--
e
~ o
1480
0-
$::
Cl)
~ 1460 'Oi;
E(z)
(50)
51440 '
The following application of these equations is suggested in order to avoid large matrices: I . Estimates of the expected values and the correlation functions are computed. 2. The matrix eqn. (50) is build. 3. The parameter vector 8 is computed using the Least Squares-solution of eqn . (50).
... . ... · desired - - - - - - measured ---simulated
1420 L~~~~_~=======:.J 34 36 38 40 42 44 46 48 50 time [s1 Fig. 4 : Measured and simulated signals of the engine speed. The 4th order model was generated using the standard LS-method .
195
In figure 4 it can be seen that the simple least squares model reproduces the plant's behavior roughly, but at the peaks the model's output is poor. Figure 5 shows, that filtering of the data vectors improves the simulation result considerably. This means, that the part of noise, which cannot be represented by the I lA-filter, is located mainly in the high frequency domain. Be aware of the fact, that the filter order, the cut off frequency and the system order have to be chosen carefully. Therefore some trials are necessary for getting good results with a new type of plant.
1520
c
] 1500
13 1480
...Co.
~ 1460
t:
'on
~ 1440 1420 34
desired - - - - - - mt:asured - - - simulated 36
38
40
42
44
46
48
50
time [5] Fig. 7: Measured and Simulated signals of the engine speed. The 4lh order model was generated using the minimalization of the cross correlation function between the input and the equation error
1520 t:
] 1500
Finally. figure 7 shows that the cross correlation minimalization method yields simulation results, which are comparable to those shown in figure 5 and 6. The great advantage of this method is that the only parameters to be chosen (tmin and t max) are insensitive to changes. Their choice over a wide range influences the simulation result negligibly .
13
...Co. 14RO
~ 1460
t: '8)
i 3 1440 1420 L~~~~_~=======.J 34 36 40 42 44 46 48 50
time [s] Fig 5: Measured and simulated signals of the engine speed. The 4lh order model was generated using the standard LS-method with filtering the data vectors.
5. CONCLUSION
In this paper different linear identification methods are compared to the standard least squares method for their applicability in modeling an industrial combustion engine plant. As the plant is unstable, and there are many unmodeled dynamics and noise in the data sets, numerous methods suggested in literature fail. However, three algorithms are presented, which yield a drastic improvement compared to the result obtained by the standard least squares method .
1520
c
.- 1500 ~ :; 1480
......
0-
~
1460
t:
'on ~
. ... . desired - - - - - - measured - - - simulated
1440
142~4
36
38
40
42
44
46
48
ACKNOWLEDGMENTS 50
Many thanks to lenbacher Energiesysteme AG , who made the data sets available to us. The project was supported by the Austrian Forschungs-ForderungsFonds der gewerblichen Wirtschaft (FFF).
time [s1 Fig. 6: Measured and simulated signals of the engine speed. The 4lh order model was generated by frequency domain order reduction of a 40 lh order LS-model
REFERENCES Figure 6 visualizes, that order reduction in the frequency domain also yields good approximation. Therefore, high order standard least squares models are able to reproduce the system behavior even if the noise filter specification is violated . Furthermore, it is shown that the order can be reduced in the frequency domain without deteriorating the simulation result. As in the case of data filtering, careful tuning of some parameters is necessary .
Isermann, R. (1992a). Identifikation dynamischer Systeme . Vol. I, 2 nd print, Springer, Berlin . Isermann, R. (1992b). Identifikation dynamischer S.vsteme . Vol. 2, 2 nd print, Springer, Berlin. Kay S.M. (1988) . Modern spectral estimation : theory and application. Prentice-Hall, London . Kouvaritakis, B.; Rossiter, I.A . (1992) . Use of bicausal weighting sequences in least squares identification of open loop unstable dynamic systems. lEE Proceedings 139, No. 3, pp328-336
196