- Email: [email protected]
A simple method for robust control design, application on a non-linear and delayed system: engine torque control
ARTICLE IN PRESS
Control Engineering Practice 12 (2004) 417–429
A simple method for robust control design, application on a non-linear and delayed system: engine torque control Y. Chamaillard*, P. Higelin, A. Charlet L.M.E., University of Orleans, 8 rue l!eonard de Vinci, Orleans 45072, France Received 2 July 2002; accepted 1 May 2003
Abstract This paper presents a simple method for designing a robust controller which can be used on uncertain and non-linear systems. As an extension of the method, the case of delayed systems is examined. The successive steps of the method are: non-linearity analysis and description, system linearization and identification, operating model design by minimization of a weighted frequencies criteria, model uncertainty design, PI and LQ controller design with robust property verification. A Smith predictor is included in the global control scheme to take the delay of the system into account. The application is a design of a robust controller for engine torque control. r 2003 Elsevier Ltd. All rights reserved. Keywords: Robust control; Delay compensation; Uncertainty model; Engine torque control
1. Introduction Today, more and more stringent pollutant emission regulations as well as user demand for better fuel economy drives the automotive industry towards the development of internal combustion engines based on low pollutant emissions and high-efficiency concepts. These new generations of engines (i.e. GDI ) will have to switch between several operating modes (i.e. stratified, homogeneous charge) to achieve their best performance in terms of emissions and efficiency. To be able to switch operating mode (or change gear in a robotized gearbox) without compromising the driver’s comfort, the switching must be performed at constant torque. To achieve this goal, a torquebased engine control system must be developed. Current solutions are not optimal (Ericksson, 1999; Pestana, 1989) or do not use today’s best perform sensors (Kawamura, Shinshi, Sato, Takahashi, & Iriyama, 1988).
*Corresponding author. Tel.: +33-238-417-050; fax: +33-238417-383. E-mail address: [email protected] (Y. Chamaillard). 0967-0661/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0967-0661(03)00113-8
The main characteristics of an engine are that it is a highly non-linear system, a delayed system with variable time delay and a system with uncertainties; all of them are problematic for the control. The main objectives of a torque control are to guarantee robust stability in the whole operating range, to ensure a small transient time with low overshoot of the engine torque but with the command as smooth as possible (because of actuator considerations). Torque-based engine controls heavily rely on accurate and high-quality feedback data. Cylinder pressure is used to extract IMEP, which is very useful for indicated torque determination. But in this case cycle to cycle variability could be problematic for the stability of the closed loop scheme. Authors have demonstrated (Jaine, Chamaillard, Charlet, & Higelin, 2002) that a moving horizon filtering method is an appropriate method to reduce the IMEP fluctuation around its mean value. Low resolution in the crank angle degree (angular coder) can also be problematic for the stability of the closed loop scheme. In a first time, a high-precision sensor (like a perfect sensor) is used to validate the proposed approach. But in Labreuche et al. (2001) authors have already studied these problems and propose solutions. It has been demonstrated that a curve fitting (Bezier type) method is
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
418
Nomenclature GDI SA SI EGR IMEP BMEP TFMEP Wi We Wf Vd Te IMEP MAP ATDC BTDC ABDC BBDC A
gazoline direct injection; spark advance; spark ignition; exhaust gas recirculation; indicated mean effective pressure; break mean effective pressure; total friction mean effective pressure; indicated work; effective work; total friction work; volume displacement; effective torque; indicated mean effective pressure; mean admission pressure; after top dead center; before top dead center; after bottom dead center; before bottom dead center; area;
an appropriate method to reconstruct precisely the IMEP. Unfortunately, indicated torque is not the global engine output (in torque control objective). But effective torque on the crankshaft is, so the friction losses are needed. Authors have investigated in Labreuche et al. (2001) on two total friction mean effective pressure models. In this paper a method to synthesize a controller is proposed. This method will be as simple as possible to have a large domain of application but ensure that the controller is robust (Guzzella, 2000). This method makes it possible to take up the challenge of torque control. In a first step, an analysis of the non-linear behavior of the system is conducted, the non-linearity is described. In a second step, the system will be linearized and identified. The nominal model is designed by minimization of a weighted frequencies criterion in the third step. The fourth step models the design uncertainty. Finally, the nominal model is used to design a model-based controller (PID, LQ, etc.) but with a robust property verification.
Cd cp cv h m N P r Q t T U V W g
Subscripts i in; o out; u upstream; d downstream
dynamic cycle’’; as it can be seen on the cylinder pressure graph (Fig. 1). The indicated work per engine cycle represents the area enclosed by the pressure–volume cycle (Heywood, 1988): I Wi ¼ p dV : ð1Þ The effective work is determined from the brake effective torque. For a four stroke engine: We ¼ 4pTe :
ð2Þ
The indicated work produced by the gas volume enclosed inside the combustion chamber is not fully available on the crankshaft. Total friction losses (piston/liner, bearings, and accessories) are defined as Wf ¼ Wi We :
ð3Þ
To be able to compare engines with different displacements, IMEP, BMEP and TFMEP are defined as follows (Heywood, 1988): Wi IMEP ¼ ; Vd BMEP ¼
2. System description—engine model In a first step, before validation on a real engine, a very complete and representative Engine Model is used (see appendix for details). This model is very accurate, because of its high-frequency resolution ‘‘in the thermo-
discharge coefficient; specific heat at constant pressure; specific heat at constant volume; specific enthalpy; mass; engine speed; pressure; ideal gas constant; heat; time; temperature; internal energy; volume; work; specific heat ratio;
We ; Vd
TFMEP ¼ BMEP ¼
Wf ; Vd
Wi TFMEP: Vd
ð4Þ
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
way of a throttle actuator with its own dynamic and non-linearity and is also subject to the filling and emptying dynamics of the intake plenum. Fortunately for each operating point, defined by an engine speed and an admitted air mass, there is an optimal SA. A map of the optimal SA is defined and used in the control. In terms of control, the engine torque becomes a single input system: the throttle angle. Obviously, the output of the system is the effective torque but in a first step the IMEP could be considered as feedback data.
12
Pressure bar
10 8 6 4 2 0 0
90 180 270 360 450 540 630 720 Crank Angle Degree Fig. 1. Cylinder pressure.
With a well-calibrated pressure transducer, an accurate angular reference (optical encoder) and a highspeed data acquisition system, IMEP can be precisely computed from the pressure data. To compute BMEP from IMEP, a friction model must be used that gives TFMEP as a function of operating parameters (Labreuche et al., 2001). Two models are proposed. First of it is a black box model that estimates friction losses in a global way. The only variable parameter is the engine rotational speed (N): " # 1 N N 2 Wf ¼ Vbdc 1 þc aþb 105 : ð5Þ s 1000 1000 Only two parameters are used to characterize the engine: s ¼ Vbdc =Vtdc volumetric efficiency, Vbdc volume at bottom dead center. Three parameters (a; b; c) have to be calibrated on engine measurements. This model is the simplest available, only one input parameter is needed (N) but the accuracy of this model can be disputed. As the previous model, the second one is a black box model that estimates the effective mean torque but it is a first-order dynamic model:
3. Non-linearity description In the next step of the method, a linearization around the operating point will be performed, so the global nonlinear behavior of the system must be analyzed to define the maximum size step between each operating point. In this case only a non-linear gain description is necessary. So using a throttle angle input from 5 to 80 with a step of 5 (or 10 ), and by gain calculus (final value theorem) the non-linear description is obtained (see Fig. 2). In Fig. 2, it can be seen that the hardest part of the non-linearity is between 10 and 50 . So the operating points, for identification purpose, are chosen as follows: * *
5–60 with a step of 5 , 60–80 with a step of 10 .
For different engine speeds the shape of the nonlinearity is almost the same, so these operating points will be taken for all engine speeds. The local excitation has been chosen to be sufficient to obtain a significant variation of the output for all operating points.
dTe dMair þ Te ¼ C1 þ C2 dt dt þ C3 l þ C4 l2 þ C5 d þ C6 d2 2
þ C7 dN þ C8 N þ C9 N :
14
ð6Þ
The input measurement variables are: N engine rotational speed (rpm), d spark advance (SA, deg.), dMair =dt air mass flow (kg/h), l equivalence ratios. Nine parameters have to be identified on the engine. This model is quite simple and takes a representative set of engine variable parameters into account. This study has been performed for an SI engine, without EGR and at stoichoimetry. An SI engine (without EGR) as a torque producer is mainly a two coupling inputs system: the admitted air mass and the SA. The SA is directly commendable, quasi instantaneously, but the air mass is only commendable by the
12 IMEP in bar
te ðNÞ
419
3500 rpm 2500 rpm 1000 rpm
10 8 6 4 2 0
0
10
20
30 40 50 60 Throttle angle in °
70
80
Fig. 2. non-linear description of IMEP for throttle angle and speed engine variations.
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
420
4. Linearization and identification In general, controller design is easier to perform using linear models and normalized variables. In the case of non-normalized variables, numerical problems can be obtained and it is hard to compare controllers to each other. For these reasons, the model has to be normalized and linearized. For each defined operating point, a model has to be identified. A classical mean square method is used. But to ensure good results the throttle input signal must have a sufficiently large spectrum.
The new variables yðtÞ and uðtÞ are dimensionless and around 1 in magnitude. During linearization only small deviations from set points are analyzed and two new variables are introduced: uðtÞ ¼ 1 þ duðtÞ; yðtÞ ¼ 1 þ dyðtÞ:
ð8Þ
The system can also be identified as a linear transfer function between duðtÞ and dyðtÞ (see Fig. 4). 4.3. Identification The general transfer function is defined as
To investigate the transfer function between the IMEP and the throttle angle a PRBS is added. PRBS is a two-level preprogrammed sequence; it has a flat spectrum, which allows exciting the system at all frequencies in the desired range for identification purposes. Register length and update rate has to be optimized. PRBS is added to the base value (operating point definition). The sequence is repeated three times: first to bring the system to its operating point and the last two for identification purposes (see Fig. 3). 4.2. Normalization, linearization The following description of the non-linearity is used for the normalization step: for each throttle angle ðThrottle anglen Þ; a corresponding IMEP is defined ðIMEPn Þ: Normalized operation can be defined as follow: Throttle angleðtÞ ¼ Throttle anglen uðtÞ;
100
200
300 400 500 sample indice
600
700
0
100
200
300 400 500 sample indice
600
700
0.05
0
-0.05
12
20 0
2
4
6 8 time in s
10
12
15 10
d in sample period
throttle angle in°
0
14
40
1000 rpm 2500 rpm 3500 rpm
10 8 6 4 2
5 0
0.1 0.05 0 -0.05 -0.1
Fig. 4. Normalized variables.
60
0
ð9Þ
First time delay must be determined as multiple of sample time. In this case the problem is in the operating point dependency. Fig. 5 shows the variations of the time delay in function of the engine speed and in function of the engine load. It appears that the time delay is more sensitive to the engine speed rather than to the engine load.
ð7Þ
IMEPðtÞ ¼ IMEPn yðtÞ:
IMEP in bar
dyðz1 Þ b0 þ b1 z1 þ ? þ bnb znb ¼ zd : 1 duðz Þ 1 þ a1 z1 þ ? þ ana zna
y (t) - normalised output u (t) - normalised input
4.1. Pseudo-random binary sequence (PRBS)
0 0
2
4
6 time in s
8
10
Fig. 3. Pseudo-random binary sequence.
12
0
10
20
30 40 50 60 Throttle angle in°
Fig. 5. Time delay variations.
70
80
ARTICLE IN PRESS Table 1 Time delay choice 2500 rpm
3500 rpm
12
5
4
normalised output
d
1000 rpm
1 0.8 0.6 Cxx
normalised input
Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
0.4
421
0.1 0.05 0 -0.05 -0.1 0
100 200 300 400 500 600 sample indice
700 measured estimated
0.05 0 -0.05 0
0.2
100 200 300 400 500 600 sample indice
700
Fig. 7. Normalized identification results.
0
200
300 400 500 sample indice
600
700
Fig. 6. Error correlation.
So the time delay can be considered as constant for all throttle input angles but different for each engine speed. The following mean values (Table 1) have been chosen: Identification is computed by a classical recursive least-squares (RLS) method. According to the criterion method and time delay assumption, the least-squares transfer function structure has to be defined. For each engine speed and each operating point na will vary from 1 to 10, nb from 1 to 10 (with respect to na) and d from 1 to 15. The optimal structure obtained is for na ¼ 4 or 5 and nb ¼ 4 or 5. But after a precise analysis, it appears that a structure with na ¼ 1 and nb ¼ 1 is not significantly different (quadratic error estimation criterion consideration) but is really more simplistic. Identification results, for the operating point defined by a 2500 rpm engine speed and a 45 throttle angle, can be seen in Figs. 6–8. At the same engine speed (2500 rpm), but at two additional throttle angles (5 and 80 ) the results are in Figs. 9 and 10. As it can be seen the identification has been carried out successfully. The transfer function is of the following form: dyðz1 Þ b0 þ b1 z1 ¼ zd : 1 duðz Þ 1 þ a1 z1
60 40 20 0
0
200
400 600 800 sample indice
1000
1200
15 IMEP
100
10 measured estimated
5 0
0
200
400 600 800 sample indice
1000
1200
Fig. 8. Un-normalized comparison.
normalised input
0
normalised output
-0.4
Throttle angle
-0.2
0.1 0.05 0 -0.05 -0.1 0
100 200 300 400 500 600 700 sample indice
0.04 0.02 0 -0.02 -0.04
measured estimated
0
100
200 300 400 500 600 700 sample indice
Fig. 9. Normalized identification results (5 ).
ð10Þ
The pole and time constant variations are illustrated by Fig. 11. Pole variation clearly confirms the non-linear behavior of the system, and the necessity of take into account the time delay into control is illustrated by Fig. 5.
5. Nominal model The first assumption with the method is that a nonlinear and parameter varying system can be described by a set of linear time invariant systems. This assumption can be justified only if the parameters change slower
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
0
0.1 0.05 0 -0.05 -0.1
-0.1 -0.2
normalised output
0
100
200
300 400 500 sample indice
600
imaginary
normalised input
422
700 measured estimated
0.02
-0.3
-0.5
0.01 0
-0.6
-0.01
-0.7
-0.02
0
100
200
300 400 500 sample indice
600
Nominal Plant
-0.4
0
0.2
0.4
700
0.6
0.8 real
1
1.2
1.4
Fig. 12. Nominal diagram.
Fig. 10. Normalized identification results (80 ).
Pole value
0 -50
-150
constant time
3500 2500 1000
-100 0
10
20
30 40 50 60 Throttle angle in°
0.2
70
80
3500 2500 1000
0.15 0.1 0.05 0
0
10
20
30 40 50 60 Throttle angle in°
70
80
Fig. 11. Pole value and constant time variations as function of throttle angle and engine speed.
It must be noted that systems are considered as free time delay models. If necessary the time delay, for the nominal model, will be determined as a mean value of Table 1. The nominal model will be used for the control design and should be defined as a continuous transfer function. This transfer function can be found using an optimization method (from Matlab optimization toolbox) with an appropriate criterion to minimize. The criterion will be chosen to minimize the magnitude error and the phase error between the nominal plant and the nominal model at each frequency. First, a gain is defined to ‘‘normalize’’ magnitude and phase. Second, two weighted gains are defined, for two different frequency domains, to adapt more or less the magnitude or the phase error: crit ¼
than the relevant dynamic effects that are to be controlled. If this assumption is verified, the non-linear system can be considered as an uncertain system, and described as a family of linear systems centered on a nominal model. In a Nyquist diagram, families of linear systems and nominal plant can be shown (Fig. 12). The nominal plant is defined in the Nyquist diagram by finding a medium real part (Rn ) and a medium imaginary part (In ) for each frequency o: Rn and In are defined as the average of the minimum and the maximum values of all the real and imaginary parts taken into the families of linear systems for each frequency o: Rn ðoÞ ðmaxðrealðsystems; oÞ þ minðrealðsystems; oÞÞ ; ¼ 2 In ðoÞ ðmaxðimagðsystems; oÞ þ minðimagðsystems; oÞÞ : ð11Þ ¼ 2
L X
ðKm Kaðem ðkÞÞ2 þ Kph ðeph ðkÞÞ2 Þ=L
ð12Þ
k¼1
with em ¼ Magnominal eph ¼ Phanominal
# nominal model ; M ag # plant Phanominal model ; plant
$ ðk ¼ 1Þ ¼ 104 rd=s; $ ðk ¼ LÞ ¼ 104 rd=s and: Ka ¼
maxðMagnominal plant Þ ; maxðPhanominal plant Þ
Km ¼ 50’koL=3; Km ¼ 30’k > 2L=3; Kph ¼ 10’Phanominal
plant ðkÞo
20:
A hypothesis concerning the nominal model structure is necessary. First- or second-order nominal models have been used nominal model ¼
b0 þ b1 S 1 þ a1 S
ARTICLE IN PRESS magnitude in db
Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
423
1 0.5
magnitude error indb
0 10-4
nominal plant 1storder nominal model 2nd order nominal model
10-3
10-2
10-1
100 w in rd/s
101
102
103
104
100 w in rd/s
101
102
103
104
100 w in rd/s
101
102
103
104
100 w in rd/s
101
102
103
104
0.1 0 -0.1 10-4
1st order nominal model 2nd order nominal model
10-3
10-2
10-1
phase in°
0 -50
phase error in°
-100 10-4
nominal plant 1st order nominal model 2nd order nominal model
10-3
10-2
10-1
20 1st order nominal model 2nd order nominal model
10 0 -10 10-4
10-3
10-2
10-1
Fig. 13. First- and second-order nominal model results.
Table 2 Nominal models Nominal model
Roots
6:52 þ 0:037S 8:82 þ S 192:4 þ 11:89S þ 0:043S2 258:9 þ 48:61S þ S2
Zeros
8.82
172.9
6.08 42.52
17.28 255.64
or nominal model ¼
b0 þ b1 S þ b2 S 2 : 1 þ a1 S þ a2 S 2
ð13Þ
For a 2500 rpm engine speed, Fig. 13 shows the optimization results and Table 2 resumes the parameter definitions and root values obtained for the two different structures of the nominal model. Further developments will demonstrate that the simple first-order model is sufficient for the controller.
6. Model uncertainty A simple description of the model uncertainty is necessary to ease, after controller design, the analysis of
the closed loop system, especially the robustness property. In this work the simplest method is used. The simple idea is to parameterize the uncertainty by the set: 9 8 > = < PðsÞ ð1 þ W ðsÞDðsÞÞ > ; ð14Þ P ¼ jDð j $ Þjp1 > > ; : ArgðDð joÞÞAfp; pg where PðsÞ is the nominal model, centered in uncertain set, W ðsÞ is a weighting function (upper bound on the modulus of the uncertainty) and DðsÞ a variable. The set P is best visualized in a Nyquist diagram. For each frequency o the set Pð joÞ is a circle of radius jPð joÞW ð joÞj and centered at Pð joÞ: W ðsÞ has to be simple as possible but chosen to cover all the uncertain sets, while minimizing the overlap. Note that only the magnitude of W ðsÞ is important, and that in general, systems can be considered as low pass, W ðsÞ can be design as follows: *
*
First, consider W ðsÞ as a simple gain and determine the static gain of W ðsÞ to cover the uncertainty at static. Second, consider W ðsÞ as a simple first-order lowpass filter, where its constant time is determined for a specific frequency chosen in the middle of the frequency domain.
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
424 Table 3 Weighting function computation
Optimal SA
Desired
Gain
o-0 rd/s
o ¼ 4 rd/s
Min Max Nominal Delta (max) Kw t
0.1503 1.3463 0.7387 0.6086 0.8238
0.1498 1.1289 0.6724 0.5226 0.0879
0.8 0.6 0.4 0.2
SA
Controller Throttle Engine angle
Torque
Engine
Effective Torque
Fig. 15. Torque engine control scheme.
In this case the engine torque control is a single-input single-output control problem. The nominal model will be used to design a controller in a classical way for a linear system. The uncertainty model will be used for robustness verification. But two more problems have not been taken into account: the saturation of the command and the time delay in the system.
0 -0.2
7.1. PI(D) controller
-0.4 -0.6 -0.8 -0.2
0
0.2
0.4
0.6 0.8 real
1
1.2
1.4
Fig. 14. Uncertainty model.
For example, in the case of engine torque control, choices for W ðsÞ are defined by Table 3. With delta ¼ maxðnominal min; max nominalÞ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !!# u" u 1 Kw2 t : t¼ 1 þ o2 delta=nominal 2
Kp ¼ 2:66; Ti ¼ 0:14s;
So the weighting function is defined by W ðsÞ ¼
0:8238 : 1 þ 0:0879S
In the case of the PID controller the system is considered as a simple linear model defined by the nominal model (see 0). So the basic principle of Ziegler and Nichols is used to tune a PI controller. The PI control scheme (Fig. 16) includes normalization (Fig. 17), error construction, PI controller, Anti Reset Windup technique (because of the saturation of the command between 0 and p=2), and the ‘‘denormalization’’ step (Fig. 18). Note that the non-linearity description is important for this last step. The parameters obtained are the following:
ð15Þ
Fig. 14 shows the uncertainty, the nominal model, and the uncertainty model (conservative bound). As it can be seen, the uncertainty model completely covers all possible plants.
7. Torque engine control The control scheme adopted for engine torque control is illustrated by Fig. 15: It must be underlined that the SA is defined by the optimal advance map. In this case the efficiency (and consumption) is privileged. Pollutant problems are essentially solved by the air fuel ratio control in the case of a three-way catalytic exhaust. So the SA directly depends on the throttle angle and on the engine speed and that solution avoids to solve a multi-variable control problem.
Tarw ¼ 1=Ti: If control does not take into account of the time delay, it is sure that an instability problem occurs (see Fig. 20). The most classical way to compensate the problem caused by the delay of the system is used: a Smith predictor is included (Fig. 19) in the feedback. But it is necessary to have a model of the system. The nominal model is used again; this is another advantage of the proposed method. In a first approach, the delay is considered as a constant equal to five sample periods. In a second approach, it will be considered as variable with the engine speed (Fig. 5). In a last step, robustness must be verified. The idea is to use the bound W ðsÞ (Section 6) to test the uncertain system for closed loop stability. Assuming a nominal plant PðsÞ; a conservative bound on its uncertainty W ðsÞ and a controller CðsÞ; assuming the nominal closed loop is asymptotically stable; then the controller CðsÞ will produce an asymptotically stable closed loop system for all the plants in the set P if the following inequality is
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
Non linearity Descripation
Torque Nominal Normalization Measured
Error Construction
425
Throttle angle PID with ARW
Denormalization
Non Saturated Command
Saturated Saturation Command
Fig. 16. PID control scheme.
This figure shows explicitly the effectiveness of the smith predictor. Figs. 21 and 22 present results, with Smith predictor, for the two opposite operating points (10 and 80 ). It appears that the dynamic is different but stable and preferment in all the cases (10 , 35 , 80 ). In the case of a saturation of the command, responses can be compared with and without ARW technique. Fig. 22 clearly demonstrates the utility of the ARW feedback. These results must be confirmed based on the very accurate non-linear model. Fig. 23 presents closed loop response with the PI controller design using the nominal model, with or without Smith predictor. Using the non-linear model, the delay problem does not appear as critical as in the linear case because of identification procedure. But an instability tendency can also be seen if the Smith predictor is omitted.
/ + * -
Fig. 17. Normalization.
* + * +
Fig. 18. Denormalization.
verified for all frequencies o: j1 þ CðjoÞjXjCðjoÞPðjoÞW ðjoÞj:
ð16Þ
This method to verify robustness is very useful and very simple because it involves only one test. If the test of robustness fails, the problem occurs because of a family of two non-homogenous linear systems. In this case, a solution is to cut the family in two (or three) subfamilies and to define two (or three) suitable linear controllers. Of course robustness is verified in the case of the PI controller for engine torque control. 7.2. Engine torque control with PI controller results In all cases the PI controller is designed with the nominal model. Opposite linear models chosen in the family set of linear models have been used to verify the robustness of the closed loop system. The necessity of the ARW technique and of the Smith predictor is also demonstrated. Measurement noise is added. Three linear models are used, defined for the 10 , 35 and 80 throttle angle operating point. First results are for the 35 operating point which is the nearest of the nominal system. Fig. 20 presents results obtained with and without the Smith predictor:
7.3. Feedforward control with PI controller Another exploitable way is to use the nominal model in a feedforward control scheme (see Fig. 24). By inversion (if possible) of the nominal model, a direct command is defined and the PI controller is used only for relative correction. The results obtained in this case are not satisfactory because of the zero in the nominal transfer function. A low-pass filter must be associated in the direct way (the time constant is about 0.01 s). The results compared to feedback control data are shown in Fig. 25. It appears that the nominal model is too distant from the real system for some operating points. This is because only one set (and only one nominal model) has been considered for all the operating points. If the family set of linear models is cut in 2 or 3 subsets, perhaps the results could be improved.
7.4. LQ controller Considering that the nominal model is a simple firstorder model, a direct state space representation is
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
426
Denormalization Free delay Nominal model
+ +
Normalization
-
Denormalization
Nominal model with delay
+
Fig. 19. Smith predictor.
11
desired with out Smith predictor with Smith preditor
With Smith predictor Desired With out Smith predictor
10 9 pressure in bar
13 12.5 12 11.5 11 10.5 10 9.5 9 8.5 8
8 7 6 5 4 3
0
1
2 3 time in s
4
2
5
1
Fig. 20. 35 with and without Smith predictor.
2
3
4 time in s
5
6
7
Fig. 23. Non-linear with Smith predictor.
1.25 desired with Smith preditor
Nominal model Inverse
1.15
+ Saturation
1.1
+
1.05
PID with ARW
1 0.95 0.9
Smith predictor
0.85 0
2
4
6 time in s
8
10 Fig. 24. Feedforward control.
16
desired with out ARW with ARW
15 14
pressure in bar
Fig. 21. 10 with Smith predictor. feedback control Desired feedforward control
10 8 6 4 1
2
3
4
5
6
7
13 12 11 10
0
2
4 6 time in s
8
Fig. 22. 80 with and without ARW.
10
throttle angle
pressure in bar
pressure in bar
1.2
feedback control feedforward control
100 50 0
0
1
2
3
4
5
6
7
Fig. 25. Results with feedforward control.
8
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
427
Table 4 Comparative criterion results
defined as X’ ¼ AX þ BU ; Y ¼ CX :
2500 rpm
ð17Þ
To obtain an optimal command as LQ definition for minimizing the energy associated with the error signal and to the command, the criterion J is used: Z N J¼ ðxT Qx uT RuÞ ð18Þ
With Smith Without Smith Feedforward LQ
3500 rpm
Qe
Sd
Qe
Sd
0.29 0.3 0.41 0.38
7.78 7.80 9.73 8.26
0.2 0.25 0.27 0.25
8.47 9.14 9.93 8.57
0
with
"
x¼
# Torque R ; Torque
u ¼ ½Throttle angle : The state space representation is extended considering the integral of the error on the torque to compensate static error. The objective is also to find a feedback gain, which minimizes the J criterion. K is obtained by solving the Riccati equation: AT P þ PA þ Q PBR1 BT P ¼ 0 and K ¼ R1 BT P with
"
Q¼
0:01 0
ð19Þ
# 0 ; 1000
8. Conclusion
R ¼ ½0:5 : The results comparative to feedback control with PI are shown in Fig. 26: For a better comparison two comparative criteria have been computed:
pressure in bar
quadratic error ðQeÞ between desired value and output value; standard deviation ðSd Þ on the throttle angle command. PI & Smith control Desired LQ & Smith control
10 8 6 4 1
2
3
4
5
throttle angle
6
7
PI & Smith control LQ & Smith control
80 60 40 20 0
0
1
2
The first is a synthetic criterion to analyze the global performances of the regulation, where the second is a good criterion to analyze the solicitations of the actuator (results are in Table 4). The PI controller with ARW scheme associated to a Smith predictor appears as a good compromise for the engine torque control. The LQ controller, as it is known, ensures a good regulation too but with a smoother command. The complete effectiveness of the LQ controller can only be evaluated comparatively to the PI controller in a real case with noisy feedback data or in a solution taking into account other inputs like EGR.
3
4
5
6
Fig. 26. Results with LQ control.
7
8
In this paper a method to synthesize a robust controller is proposed. This method is as simple as possible but ensures an efficient control for non-linear systems and time delayed systems. By an inequality verification, the robustness of the closed loop uncertain system is tested. The problems of non-linearity have been solved by a normalization step, the time delay problem by a Smith predictor. An important result of the method is that a linear nominal model of the system is established with uncertainty description. This linear model is a multi-function model; it is very helpful for the controller definition, for delay compensation (Smith predictor) or for a feedforward control (if the model is invertible). First the global behavior of the system has to be described, second the system has to be linearized and identified. This linear nominal model is also obtained by minimizing a weighted frequencies criterion. The method has been illustrated on an explicit example: engine torque control. A very precise model is used to simulate the engine. Two controllers have been defined: a PI controller and an LQ controller. In the case of the PI controller two different schemes have been studied: only feedback control or feedforward and feedback control.
ARTICLE IN PRESS 428
Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
It must be underlined that all controllers have been determined with and for the linear nominal model. It is exactly the same controller that has been used with the non-linear system. The best results have been obtained for a PI controller associated to a Smith predictor. In the case of the feedforward control using the nominal inverse model, the non-linearity is too large (the nominal model is in some cases too distant from the real system). If the family set of linear models is cut as 2 or 3 subsets, perhaps the results could be improved. The LQ controller is not very explicit in this case (of course) because it will operate better in a multi-input configuration (Throttle, EGR, etc.) or in the case of signals with high-frequency noise. It can be noted that it will be the case on real engine. In all the cases, robust stability is guaranteed in the whole range of engine torque control. In closed loop small transient time and very low overshoot are achieved. In future work, the same method will be used on a real engine to prove its effectiveness. In this case multi-inputs must be considered like EGR or turbocharger. To analyze more finely the effectiveness of this method, it should be compared, on the same application, to a robust control design or to a non-linear control like sliding mode control.
Appendix A. The engine model This model represents a one cylinder, half-liter engine. All dimensions are those of a J4S Renault research engine (Table 5). Initial model validations have been performed by comparison with test bed results. The model which has been developed is a quasidimensional model (Heywood, 1988) where volumes (tank, intake, etc.) are connected together via flow singularities (valve, throttle, etc.). Due to the nature of the model, pipe acoustics as well as gas inertia have been
Table 5 Specifications of the engine Number of cylinders Cycle Number of intake valves Number of exhaust valves Bore Stroke Combustion chamber design Displacement Connecting rod length Compression ratio Piston Intake valve O/C Exhaust valve O/C
1 4 stroke 2 2 88.00 mm (3.45 in) 82.00 mm (3.23 in) Pentroof 498.7 cm3 136.9 mm (5.39 in) 9.5:1 Plate 37 BTDC/37 ABDC 37 BBDC/37 ATDC
neglected (Heywood, 1988). The model has been split into several components: valve lift (intake, exhaust, charging valve), pressure drop (intake, exhaust, charging valve), cylinder (variable volume vessel), air tank (constant volume vessel). A.1. Valve lift The valve lift determines the effective section of the valve flow restriction. It could represent a conventional camshaft (lift as a function of crankshaft angle) or a camless valve actuator (lift as a result of the valvetrain force balance). A.2. Pressure drop The pressure drop is calculated with Barre! de Saint Venant equation. Four cases must be distinguished. The flow can be in the forward or reverse direction (back flow) and it can be sub or supersonic. Eq (1) presents the critical ratio which determines the type of the flow. If the pressure drop is lower than the critical ratio, the flow is supersonic, otherwise the flow is subsonic.
2 gþ1
g=g1 :
ðA:1Þ
Eq. (19) represents the mass flow in the supersonic case. Eq. (A.1) represents the mass flow in the subsonic case. Cd is the discharge coefficient and Ar is the area of the singularity. This area is calculated from the valve lift. gþ1=2ðgþ1Þ dm Pu 2 ¼ Cd Ar pffiffiffiffiffiffiffiffi g0:5 ; dt gþ1 rTu ! dm Pu Pd 1=g ¼ Cd Ar pffiffiffiffiffiffiffi dt Pd rTu vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u g1=g ! u 2g Pd : t 1 g1 Pd
ðA:2Þ
ðA:3Þ
A.3. Cylinder The cylinder is the central part of the system. It exchanges work, heat and mass with outside. The first law of thermodynamics can by written as follows: dU ¼ dW þ dQ þ hi dmi þ ho dmo ;
ðA:4Þ
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429 1e5 0
Open
200
Close
296.5 Area
In take Valve P in P out T in T out Area
Intake Valve Area Open
90
Close
Area
Air Tank Valve Area Open
-45
Close
1e 5 29 6.5
Area
P in P out T int T out Area
Exhaust Valve Area
dm
P
dU V Cylinder
Volume
10 Advance 1.4
O
m
T
Cylinder Volume
A
dm in dm out dh in dh out Sonic ?
Exhaust Valve
Angla
Angle
dm in dm out dh in dh out Sonic ?
Air Tank Valve
Angle
-250
dm in dm out dh int dh out Sonic ?
P in P out T in T out Area
Angle
0
429
A Av m lamda O
dQ
Vibé 2 phases
Crankshaft Angle
N 1200
Engine Speed
Fig. 27. Engine model.
where U is the internal energy, W is the work generated by the moving walls, Q is the heat, m the mass and h the enthalpy. i indicates the intake flow, o indicates the outlet flow. For an ideal gas, the internal energy is dU ¼ mcv dT þ cv T dm
from the input mass and enthalpy fluxes. Conversely, the valve blocks compute enthalpy and mass fluxes from upstream and downstream pressure and temperature. The heat and mass balance are performed by the two sum blocks.
ðA:5Þ
with Eqs. (A.2) and (A.3) the following expression can be obtained: dW þ dQ þ dhi þ dho cv T dm dT ¼ : ðA:6Þ mcv The cylinder pressure is calculated with the ideal gas law (see Eq. (A.5)). The volume is determined with the engine geometry. The mass inside the cylinder is determined by integrating the mass flow through each valve and the temperature by integrating Eq. (A.4). mrT P¼ : ðA:7Þ V A.4. Global model The global model has been realized (Higelin & Charlet, 2001) with Matlab/Simulink. Fig. 27 shows a schematic view of the model. It can be seen that the volume blocks compute pressure, temperature and mass
References Eriksson, L. (1999). Spark advance modeling and control. Thesis of the . Linkoping University. Guzzella, L. (2000). Model-based and robust control systems development. Control System Workshop, Austin, Texas. Heywood, J. B. (1988). Internal combustion engine fundamentals. New York: McGraw-Hill. Higelin, P., & Charlet, A. (2001). Thermodynamic cycles for a new hybrid pneumatic—Combustion engine concept, Fifth international conference on internal combustion engines (ICE 2001), Naples. Jaine, T., Chamaillard, Y., Charlet, A., & Higelin, P. (2002). High frequency IMEP estimation, filtering for torque based SI engine control, SAE paper 2002-01-1276. Kawamura, Y., Shinshi, M., Sato, H., Takahashi, N., & Iriyama, M. (1988). MBT control through individual cylinder pressure detection, SAE paper 881779. Labreuche, G., Da Costa, A., Chamaillard, Y., Charlet, A., Higelin, P., & Perrier, C. (2001). Total friction effective pressure and torque estimator, MECA’O1, Italy. Pestana, G. W. (1989). Engine control methods using combustion pressure feedback, SAE paper 890758.
Control Engineering Practice 12 (2004) 417–429
A simple method for robust control design, application on a non-linear and delayed system: engine torque control Y. Chamaillard*, P. Higelin, A. Charlet L.M.E., University of Orleans, 8 rue l!eonard de Vinci, Orleans 45072, France Received 2 July 2002; accepted 1 May 2003
Abstract This paper presents a simple method for designing a robust controller which can be used on uncertain and non-linear systems. As an extension of the method, the case of delayed systems is examined. The successive steps of the method are: non-linearity analysis and description, system linearization and identification, operating model design by minimization of a weighted frequencies criteria, model uncertainty design, PI and LQ controller design with robust property verification. A Smith predictor is included in the global control scheme to take the delay of the system into account. The application is a design of a robust controller for engine torque control. r 2003 Elsevier Ltd. All rights reserved. Keywords: Robust control; Delay compensation; Uncertainty model; Engine torque control
1. Introduction Today, more and more stringent pollutant emission regulations as well as user demand for better fuel economy drives the automotive industry towards the development of internal combustion engines based on low pollutant emissions and high-efficiency concepts. These new generations of engines (i.e. GDI ) will have to switch between several operating modes (i.e. stratified, homogeneous charge) to achieve their best performance in terms of emissions and efficiency. To be able to switch operating mode (or change gear in a robotized gearbox) without compromising the driver’s comfort, the switching must be performed at constant torque. To achieve this goal, a torquebased engine control system must be developed. Current solutions are not optimal (Ericksson, 1999; Pestana, 1989) or do not use today’s best perform sensors (Kawamura, Shinshi, Sato, Takahashi, & Iriyama, 1988).
*Corresponding author. Tel.: +33-238-417-050; fax: +33-238417-383. E-mail address: [email protected] (Y. Chamaillard). 0967-0661/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0967-0661(03)00113-8
The main characteristics of an engine are that it is a highly non-linear system, a delayed system with variable time delay and a system with uncertainties; all of them are problematic for the control. The main objectives of a torque control are to guarantee robust stability in the whole operating range, to ensure a small transient time with low overshoot of the engine torque but with the command as smooth as possible (because of actuator considerations). Torque-based engine controls heavily rely on accurate and high-quality feedback data. Cylinder pressure is used to extract IMEP, which is very useful for indicated torque determination. But in this case cycle to cycle variability could be problematic for the stability of the closed loop scheme. Authors have demonstrated (Jaine, Chamaillard, Charlet, & Higelin, 2002) that a moving horizon filtering method is an appropriate method to reduce the IMEP fluctuation around its mean value. Low resolution in the crank angle degree (angular coder) can also be problematic for the stability of the closed loop scheme. In a first time, a high-precision sensor (like a perfect sensor) is used to validate the proposed approach. But in Labreuche et al. (2001) authors have already studied these problems and propose solutions. It has been demonstrated that a curve fitting (Bezier type) method is
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
418
Nomenclature GDI SA SI EGR IMEP BMEP TFMEP Wi We Wf Vd Te IMEP MAP ATDC BTDC ABDC BBDC A
gazoline direct injection; spark advance; spark ignition; exhaust gas recirculation; indicated mean effective pressure; break mean effective pressure; total friction mean effective pressure; indicated work; effective work; total friction work; volume displacement; effective torque; indicated mean effective pressure; mean admission pressure; after top dead center; before top dead center; after bottom dead center; before bottom dead center; area;
an appropriate method to reconstruct precisely the IMEP. Unfortunately, indicated torque is not the global engine output (in torque control objective). But effective torque on the crankshaft is, so the friction losses are needed. Authors have investigated in Labreuche et al. (2001) on two total friction mean effective pressure models. In this paper a method to synthesize a controller is proposed. This method will be as simple as possible to have a large domain of application but ensure that the controller is robust (Guzzella, 2000). This method makes it possible to take up the challenge of torque control. In a first step, an analysis of the non-linear behavior of the system is conducted, the non-linearity is described. In a second step, the system will be linearized and identified. The nominal model is designed by minimization of a weighted frequencies criterion in the third step. The fourth step models the design uncertainty. Finally, the nominal model is used to design a model-based controller (PID, LQ, etc.) but with a robust property verification.
Cd cp cv h m N P r Q t T U V W g
Subscripts i in; o out; u upstream; d downstream
dynamic cycle’’; as it can be seen on the cylinder pressure graph (Fig. 1). The indicated work per engine cycle represents the area enclosed by the pressure–volume cycle (Heywood, 1988): I Wi ¼ p dV : ð1Þ The effective work is determined from the brake effective torque. For a four stroke engine: We ¼ 4pTe :
ð2Þ
The indicated work produced by the gas volume enclosed inside the combustion chamber is not fully available on the crankshaft. Total friction losses (piston/liner, bearings, and accessories) are defined as Wf ¼ Wi We :
ð3Þ
To be able to compare engines with different displacements, IMEP, BMEP and TFMEP are defined as follows (Heywood, 1988): Wi IMEP ¼ ; Vd BMEP ¼
2. System description—engine model In a first step, before validation on a real engine, a very complete and representative Engine Model is used (see appendix for details). This model is very accurate, because of its high-frequency resolution ‘‘in the thermo-
discharge coefficient; specific heat at constant pressure; specific heat at constant volume; specific enthalpy; mass; engine speed; pressure; ideal gas constant; heat; time; temperature; internal energy; volume; work; specific heat ratio;
We ; Vd
TFMEP ¼ BMEP ¼
Wf ; Vd
Wi TFMEP: Vd
ð4Þ
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
way of a throttle actuator with its own dynamic and non-linearity and is also subject to the filling and emptying dynamics of the intake plenum. Fortunately for each operating point, defined by an engine speed and an admitted air mass, there is an optimal SA. A map of the optimal SA is defined and used in the control. In terms of control, the engine torque becomes a single input system: the throttle angle. Obviously, the output of the system is the effective torque but in a first step the IMEP could be considered as feedback data.
12
Pressure bar
10 8 6 4 2 0 0
90 180 270 360 450 540 630 720 Crank Angle Degree Fig. 1. Cylinder pressure.
With a well-calibrated pressure transducer, an accurate angular reference (optical encoder) and a highspeed data acquisition system, IMEP can be precisely computed from the pressure data. To compute BMEP from IMEP, a friction model must be used that gives TFMEP as a function of operating parameters (Labreuche et al., 2001). Two models are proposed. First of it is a black box model that estimates friction losses in a global way. The only variable parameter is the engine rotational speed (N): " # 1 N N 2 Wf ¼ Vbdc 1 þc aþb 105 : ð5Þ s 1000 1000 Only two parameters are used to characterize the engine: s ¼ Vbdc =Vtdc volumetric efficiency, Vbdc volume at bottom dead center. Three parameters (a; b; c) have to be calibrated on engine measurements. This model is the simplest available, only one input parameter is needed (N) but the accuracy of this model can be disputed. As the previous model, the second one is a black box model that estimates the effective mean torque but it is a first-order dynamic model:
3. Non-linearity description In the next step of the method, a linearization around the operating point will be performed, so the global nonlinear behavior of the system must be analyzed to define the maximum size step between each operating point. In this case only a non-linear gain description is necessary. So using a throttle angle input from 5 to 80 with a step of 5 (or 10 ), and by gain calculus (final value theorem) the non-linear description is obtained (see Fig. 2). In Fig. 2, it can be seen that the hardest part of the non-linearity is between 10 and 50 . So the operating points, for identification purpose, are chosen as follows: * *
5–60 with a step of 5 , 60–80 with a step of 10 .
For different engine speeds the shape of the nonlinearity is almost the same, so these operating points will be taken for all engine speeds. The local excitation has been chosen to be sufficient to obtain a significant variation of the output for all operating points.
dTe dMair þ Te ¼ C1 þ C2 dt dt þ C3 l þ C4 l2 þ C5 d þ C6 d2 2
þ C7 dN þ C8 N þ C9 N :
14
ð6Þ
The input measurement variables are: N engine rotational speed (rpm), d spark advance (SA, deg.), dMair =dt air mass flow (kg/h), l equivalence ratios. Nine parameters have to be identified on the engine. This model is quite simple and takes a representative set of engine variable parameters into account. This study has been performed for an SI engine, without EGR and at stoichoimetry. An SI engine (without EGR) as a torque producer is mainly a two coupling inputs system: the admitted air mass and the SA. The SA is directly commendable, quasi instantaneously, but the air mass is only commendable by the
12 IMEP in bar
te ðNÞ
419
3500 rpm 2500 rpm 1000 rpm
10 8 6 4 2 0
0
10
20
30 40 50 60 Throttle angle in °
70
80
Fig. 2. non-linear description of IMEP for throttle angle and speed engine variations.
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
420
4. Linearization and identification In general, controller design is easier to perform using linear models and normalized variables. In the case of non-normalized variables, numerical problems can be obtained and it is hard to compare controllers to each other. For these reasons, the model has to be normalized and linearized. For each defined operating point, a model has to be identified. A classical mean square method is used. But to ensure good results the throttle input signal must have a sufficiently large spectrum.
The new variables yðtÞ and uðtÞ are dimensionless and around 1 in magnitude. During linearization only small deviations from set points are analyzed and two new variables are introduced: uðtÞ ¼ 1 þ duðtÞ; yðtÞ ¼ 1 þ dyðtÞ:
ð8Þ
The system can also be identified as a linear transfer function between duðtÞ and dyðtÞ (see Fig. 4). 4.3. Identification The general transfer function is defined as
To investigate the transfer function between the IMEP and the throttle angle a PRBS is added. PRBS is a two-level preprogrammed sequence; it has a flat spectrum, which allows exciting the system at all frequencies in the desired range for identification purposes. Register length and update rate has to be optimized. PRBS is added to the base value (operating point definition). The sequence is repeated three times: first to bring the system to its operating point and the last two for identification purposes (see Fig. 3). 4.2. Normalization, linearization The following description of the non-linearity is used for the normalization step: for each throttle angle ðThrottle anglen Þ; a corresponding IMEP is defined ðIMEPn Þ: Normalized operation can be defined as follow: Throttle angleðtÞ ¼ Throttle anglen uðtÞ;
100
200
300 400 500 sample indice
600
700
0
100
200
300 400 500 sample indice
600
700
0.05
0
-0.05
12
20 0
2
4
6 8 time in s
10
12
15 10
d in sample period
throttle angle in°
0
14
40
1000 rpm 2500 rpm 3500 rpm
10 8 6 4 2
5 0
0.1 0.05 0 -0.05 -0.1
Fig. 4. Normalized variables.
60
0
ð9Þ
First time delay must be determined as multiple of sample time. In this case the problem is in the operating point dependency. Fig. 5 shows the variations of the time delay in function of the engine speed and in function of the engine load. It appears that the time delay is more sensitive to the engine speed rather than to the engine load.
ð7Þ
IMEPðtÞ ¼ IMEPn yðtÞ:
IMEP in bar
dyðz1 Þ b0 þ b1 z1 þ ? þ bnb znb ¼ zd : 1 duðz Þ 1 þ a1 z1 þ ? þ ana zna
y (t) - normalised output u (t) - normalised input
4.1. Pseudo-random binary sequence (PRBS)
0 0
2
4
6 time in s
8
10
Fig. 3. Pseudo-random binary sequence.
12
0
10
20
30 40 50 60 Throttle angle in°
Fig. 5. Time delay variations.
70
80
ARTICLE IN PRESS Table 1 Time delay choice 2500 rpm
3500 rpm
12
5
4
normalised output
d
1000 rpm
1 0.8 0.6 Cxx
normalised input
Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
0.4
421
0.1 0.05 0 -0.05 -0.1 0
100 200 300 400 500 600 sample indice
700 measured estimated
0.05 0 -0.05 0
0.2
100 200 300 400 500 600 sample indice
700
Fig. 7. Normalized identification results.
0
200
300 400 500 sample indice
600
700
Fig. 6. Error correlation.
So the time delay can be considered as constant for all throttle input angles but different for each engine speed. The following mean values (Table 1) have been chosen: Identification is computed by a classical recursive least-squares (RLS) method. According to the criterion method and time delay assumption, the least-squares transfer function structure has to be defined. For each engine speed and each operating point na will vary from 1 to 10, nb from 1 to 10 (with respect to na) and d from 1 to 15. The optimal structure obtained is for na ¼ 4 or 5 and nb ¼ 4 or 5. But after a precise analysis, it appears that a structure with na ¼ 1 and nb ¼ 1 is not significantly different (quadratic error estimation criterion consideration) but is really more simplistic. Identification results, for the operating point defined by a 2500 rpm engine speed and a 45 throttle angle, can be seen in Figs. 6–8. At the same engine speed (2500 rpm), but at two additional throttle angles (5 and 80 ) the results are in Figs. 9 and 10. As it can be seen the identification has been carried out successfully. The transfer function is of the following form: dyðz1 Þ b0 þ b1 z1 ¼ zd : 1 duðz Þ 1 þ a1 z1
60 40 20 0
0
200
400 600 800 sample indice
1000
1200
15 IMEP
100
10 measured estimated
5 0
0
200
400 600 800 sample indice
1000
1200
Fig. 8. Un-normalized comparison.
normalised input
0
normalised output
-0.4
Throttle angle
-0.2
0.1 0.05 0 -0.05 -0.1 0
100 200 300 400 500 600 700 sample indice
0.04 0.02 0 -0.02 -0.04
measured estimated
0
100
200 300 400 500 600 700 sample indice
Fig. 9. Normalized identification results (5 ).
ð10Þ
The pole and time constant variations are illustrated by Fig. 11. Pole variation clearly confirms the non-linear behavior of the system, and the necessity of take into account the time delay into control is illustrated by Fig. 5.
5. Nominal model The first assumption with the method is that a nonlinear and parameter varying system can be described by a set of linear time invariant systems. This assumption can be justified only if the parameters change slower
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
0
0.1 0.05 0 -0.05 -0.1
-0.1 -0.2
normalised output
0
100
200
300 400 500 sample indice
600
imaginary
normalised input
422
700 measured estimated
0.02
-0.3
-0.5
0.01 0
-0.6
-0.01
-0.7
-0.02
0
100
200
300 400 500 sample indice
600
Nominal Plant
-0.4
0
0.2
0.4
700
0.6
0.8 real
1
1.2
1.4
Fig. 12. Nominal diagram.
Fig. 10. Normalized identification results (80 ).
Pole value
0 -50
-150
constant time
3500 2500 1000
-100 0
10
20
30 40 50 60 Throttle angle in°
0.2
70
80
3500 2500 1000
0.15 0.1 0.05 0
0
10
20
30 40 50 60 Throttle angle in°
70
80
Fig. 11. Pole value and constant time variations as function of throttle angle and engine speed.
It must be noted that systems are considered as free time delay models. If necessary the time delay, for the nominal model, will be determined as a mean value of Table 1. The nominal model will be used for the control design and should be defined as a continuous transfer function. This transfer function can be found using an optimization method (from Matlab optimization toolbox) with an appropriate criterion to minimize. The criterion will be chosen to minimize the magnitude error and the phase error between the nominal plant and the nominal model at each frequency. First, a gain is defined to ‘‘normalize’’ magnitude and phase. Second, two weighted gains are defined, for two different frequency domains, to adapt more or less the magnitude or the phase error: crit ¼
than the relevant dynamic effects that are to be controlled. If this assumption is verified, the non-linear system can be considered as an uncertain system, and described as a family of linear systems centered on a nominal model. In a Nyquist diagram, families of linear systems and nominal plant can be shown (Fig. 12). The nominal plant is defined in the Nyquist diagram by finding a medium real part (Rn ) and a medium imaginary part (In ) for each frequency o: Rn and In are defined as the average of the minimum and the maximum values of all the real and imaginary parts taken into the families of linear systems for each frequency o: Rn ðoÞ ðmaxðrealðsystems; oÞ þ minðrealðsystems; oÞÞ ; ¼ 2 In ðoÞ ðmaxðimagðsystems; oÞ þ minðimagðsystems; oÞÞ : ð11Þ ¼ 2
L X
ðKm Kaðem ðkÞÞ2 þ Kph ðeph ðkÞÞ2 Þ=L
ð12Þ
k¼1
with em ¼ Magnominal eph ¼ Phanominal
# nominal model ; M ag # plant Phanominal model ; plant
$ ðk ¼ 1Þ ¼ 104 rd=s; $ ðk ¼ LÞ ¼ 104 rd=s and: Ka ¼
maxðMagnominal plant Þ ; maxðPhanominal plant Þ
Km ¼ 50’koL=3; Km ¼ 30’k > 2L=3; Kph ¼ 10’Phanominal
plant ðkÞo
20:
A hypothesis concerning the nominal model structure is necessary. First- or second-order nominal models have been used nominal model ¼
b0 þ b1 S 1 þ a1 S
ARTICLE IN PRESS magnitude in db
Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
423
1 0.5
magnitude error indb
0 10-4
nominal plant 1storder nominal model 2nd order nominal model
10-3
10-2
10-1
100 w in rd/s
101
102
103
104
100 w in rd/s
101
102
103
104
100 w in rd/s
101
102
103
104
100 w in rd/s
101
102
103
104
0.1 0 -0.1 10-4
1st order nominal model 2nd order nominal model
10-3
10-2
10-1
phase in°
0 -50
phase error in°
-100 10-4
nominal plant 1st order nominal model 2nd order nominal model
10-3
10-2
10-1
20 1st order nominal model 2nd order nominal model
10 0 -10 10-4
10-3
10-2
10-1
Fig. 13. First- and second-order nominal model results.
Table 2 Nominal models Nominal model
Roots
6:52 þ 0:037S 8:82 þ S 192:4 þ 11:89S þ 0:043S2 258:9 þ 48:61S þ S2
Zeros
8.82
172.9
6.08 42.52
17.28 255.64
or nominal model ¼
b0 þ b1 S þ b2 S 2 : 1 þ a1 S þ a2 S 2
ð13Þ
For a 2500 rpm engine speed, Fig. 13 shows the optimization results and Table 2 resumes the parameter definitions and root values obtained for the two different structures of the nominal model. Further developments will demonstrate that the simple first-order model is sufficient for the controller.
6. Model uncertainty A simple description of the model uncertainty is necessary to ease, after controller design, the analysis of
the closed loop system, especially the robustness property. In this work the simplest method is used. The simple idea is to parameterize the uncertainty by the set: 9 8 > = < PðsÞ ð1 þ W ðsÞDðsÞÞ > ; ð14Þ P ¼ jDð j $ Þjp1 > > ; : ArgðDð joÞÞAfp; pg where PðsÞ is the nominal model, centered in uncertain set, W ðsÞ is a weighting function (upper bound on the modulus of the uncertainty) and DðsÞ a variable. The set P is best visualized in a Nyquist diagram. For each frequency o the set Pð joÞ is a circle of radius jPð joÞW ð joÞj and centered at Pð joÞ: W ðsÞ has to be simple as possible but chosen to cover all the uncertain sets, while minimizing the overlap. Note that only the magnitude of W ðsÞ is important, and that in general, systems can be considered as low pass, W ðsÞ can be design as follows: *
*
First, consider W ðsÞ as a simple gain and determine the static gain of W ðsÞ to cover the uncertainty at static. Second, consider W ðsÞ as a simple first-order lowpass filter, where its constant time is determined for a specific frequency chosen in the middle of the frequency domain.
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
424 Table 3 Weighting function computation
Optimal SA
Desired
Gain
o-0 rd/s
o ¼ 4 rd/s
Min Max Nominal Delta (max) Kw t
0.1503 1.3463 0.7387 0.6086 0.8238
0.1498 1.1289 0.6724 0.5226 0.0879
0.8 0.6 0.4 0.2
SA
Controller Throttle Engine angle
Torque
Engine
Effective Torque
Fig. 15. Torque engine control scheme.
In this case the engine torque control is a single-input single-output control problem. The nominal model will be used to design a controller in a classical way for a linear system. The uncertainty model will be used for robustness verification. But two more problems have not been taken into account: the saturation of the command and the time delay in the system.
0 -0.2
7.1. PI(D) controller
-0.4 -0.6 -0.8 -0.2
0
0.2
0.4
0.6 0.8 real
1
1.2
1.4
Fig. 14. Uncertainty model.
For example, in the case of engine torque control, choices for W ðsÞ are defined by Table 3. With delta ¼ maxðnominal min; max nominalÞ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !!# u" u 1 Kw2 t : t¼ 1 þ o2 delta=nominal 2
Kp ¼ 2:66; Ti ¼ 0:14s;
So the weighting function is defined by W ðsÞ ¼
0:8238 : 1 þ 0:0879S
In the case of the PID controller the system is considered as a simple linear model defined by the nominal model (see 0). So the basic principle of Ziegler and Nichols is used to tune a PI controller. The PI control scheme (Fig. 16) includes normalization (Fig. 17), error construction, PI controller, Anti Reset Windup technique (because of the saturation of the command between 0 and p=2), and the ‘‘denormalization’’ step (Fig. 18). Note that the non-linearity description is important for this last step. The parameters obtained are the following:
ð15Þ
Fig. 14 shows the uncertainty, the nominal model, and the uncertainty model (conservative bound). As it can be seen, the uncertainty model completely covers all possible plants.
7. Torque engine control The control scheme adopted for engine torque control is illustrated by Fig. 15: It must be underlined that the SA is defined by the optimal advance map. In this case the efficiency (and consumption) is privileged. Pollutant problems are essentially solved by the air fuel ratio control in the case of a three-way catalytic exhaust. So the SA directly depends on the throttle angle and on the engine speed and that solution avoids to solve a multi-variable control problem.
Tarw ¼ 1=Ti: If control does not take into account of the time delay, it is sure that an instability problem occurs (see Fig. 20). The most classical way to compensate the problem caused by the delay of the system is used: a Smith predictor is included (Fig. 19) in the feedback. But it is necessary to have a model of the system. The nominal model is used again; this is another advantage of the proposed method. In a first approach, the delay is considered as a constant equal to five sample periods. In a second approach, it will be considered as variable with the engine speed (Fig. 5). In a last step, robustness must be verified. The idea is to use the bound W ðsÞ (Section 6) to test the uncertain system for closed loop stability. Assuming a nominal plant PðsÞ; a conservative bound on its uncertainty W ðsÞ and a controller CðsÞ; assuming the nominal closed loop is asymptotically stable; then the controller CðsÞ will produce an asymptotically stable closed loop system for all the plants in the set P if the following inequality is
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
Non linearity Descripation
Torque Nominal Normalization Measured
Error Construction
425
Throttle angle PID with ARW
Denormalization
Non Saturated Command
Saturated Saturation Command
Fig. 16. PID control scheme.
This figure shows explicitly the effectiveness of the smith predictor. Figs. 21 and 22 present results, with Smith predictor, for the two opposite operating points (10 and 80 ). It appears that the dynamic is different but stable and preferment in all the cases (10 , 35 , 80 ). In the case of a saturation of the command, responses can be compared with and without ARW technique. Fig. 22 clearly demonstrates the utility of the ARW feedback. These results must be confirmed based on the very accurate non-linear model. Fig. 23 presents closed loop response with the PI controller design using the nominal model, with or without Smith predictor. Using the non-linear model, the delay problem does not appear as critical as in the linear case because of identification procedure. But an instability tendency can also be seen if the Smith predictor is omitted.
/ + * -
Fig. 17. Normalization.
* + * +
Fig. 18. Denormalization.
verified for all frequencies o: j1 þ CðjoÞjXjCðjoÞPðjoÞW ðjoÞj:
ð16Þ
This method to verify robustness is very useful and very simple because it involves only one test. If the test of robustness fails, the problem occurs because of a family of two non-homogenous linear systems. In this case, a solution is to cut the family in two (or three) subfamilies and to define two (or three) suitable linear controllers. Of course robustness is verified in the case of the PI controller for engine torque control. 7.2. Engine torque control with PI controller results In all cases the PI controller is designed with the nominal model. Opposite linear models chosen in the family set of linear models have been used to verify the robustness of the closed loop system. The necessity of the ARW technique and of the Smith predictor is also demonstrated. Measurement noise is added. Three linear models are used, defined for the 10 , 35 and 80 throttle angle operating point. First results are for the 35 operating point which is the nearest of the nominal system. Fig. 20 presents results obtained with and without the Smith predictor:
7.3. Feedforward control with PI controller Another exploitable way is to use the nominal model in a feedforward control scheme (see Fig. 24). By inversion (if possible) of the nominal model, a direct command is defined and the PI controller is used only for relative correction. The results obtained in this case are not satisfactory because of the zero in the nominal transfer function. A low-pass filter must be associated in the direct way (the time constant is about 0.01 s). The results compared to feedback control data are shown in Fig. 25. It appears that the nominal model is too distant from the real system for some operating points. This is because only one set (and only one nominal model) has been considered for all the operating points. If the family set of linear models is cut in 2 or 3 subsets, perhaps the results could be improved.
7.4. LQ controller Considering that the nominal model is a simple firstorder model, a direct state space representation is
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
426
Denormalization Free delay Nominal model
+ +
Normalization
-
Denormalization
Nominal model with delay
+
Fig. 19. Smith predictor.
11
desired with out Smith predictor with Smith preditor
With Smith predictor Desired With out Smith predictor
10 9 pressure in bar
13 12.5 12 11.5 11 10.5 10 9.5 9 8.5 8
8 7 6 5 4 3
0
1
2 3 time in s
4
2
5
1
Fig. 20. 35 with and without Smith predictor.
2
3
4 time in s
5
6
7
Fig. 23. Non-linear with Smith predictor.
1.25 desired with Smith preditor
Nominal model Inverse
1.15
+ Saturation
1.1
+
1.05
PID with ARW
1 0.95 0.9
Smith predictor
0.85 0
2
4
6 time in s
8
10 Fig. 24. Feedforward control.
16
desired with out ARW with ARW
15 14
pressure in bar
Fig. 21. 10 with Smith predictor. feedback control Desired feedforward control
10 8 6 4 1
2
3
4
5
6
7
13 12 11 10
0
2
4 6 time in s
8
Fig. 22. 80 with and without ARW.
10
throttle angle
pressure in bar
pressure in bar
1.2
feedback control feedforward control
100 50 0
0
1
2
3
4
5
6
7
Fig. 25. Results with feedforward control.
8
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
427
Table 4 Comparative criterion results
defined as X’ ¼ AX þ BU ; Y ¼ CX :
2500 rpm
ð17Þ
To obtain an optimal command as LQ definition for minimizing the energy associated with the error signal and to the command, the criterion J is used: Z N J¼ ðxT Qx uT RuÞ ð18Þ
With Smith Without Smith Feedforward LQ
3500 rpm
Qe
Sd
Qe
Sd
0.29 0.3 0.41 0.38
7.78 7.80 9.73 8.26
0.2 0.25 0.27 0.25
8.47 9.14 9.93 8.57
0
with
"
x¼
# Torque R ; Torque
u ¼ ½Throttle angle : The state space representation is extended considering the integral of the error on the torque to compensate static error. The objective is also to find a feedback gain, which minimizes the J criterion. K is obtained by solving the Riccati equation: AT P þ PA þ Q PBR1 BT P ¼ 0 and K ¼ R1 BT P with
"
Q¼
0:01 0
ð19Þ
# 0 ; 1000
8. Conclusion
R ¼ ½0:5 : The results comparative to feedback control with PI are shown in Fig. 26: For a better comparison two comparative criteria have been computed:
pressure in bar
quadratic error ðQeÞ between desired value and output value; standard deviation ðSd Þ on the throttle angle command. PI & Smith control Desired LQ & Smith control
10 8 6 4 1
2
3
4
5
throttle angle
6
7
PI & Smith control LQ & Smith control
80 60 40 20 0
0
1
2
The first is a synthetic criterion to analyze the global performances of the regulation, where the second is a good criterion to analyze the solicitations of the actuator (results are in Table 4). The PI controller with ARW scheme associated to a Smith predictor appears as a good compromise for the engine torque control. The LQ controller, as it is known, ensures a good regulation too but with a smoother command. The complete effectiveness of the LQ controller can only be evaluated comparatively to the PI controller in a real case with noisy feedback data or in a solution taking into account other inputs like EGR.
3
4
5
6
Fig. 26. Results with LQ control.
7
8
In this paper a method to synthesize a robust controller is proposed. This method is as simple as possible but ensures an efficient control for non-linear systems and time delayed systems. By an inequality verification, the robustness of the closed loop uncertain system is tested. The problems of non-linearity have been solved by a normalization step, the time delay problem by a Smith predictor. An important result of the method is that a linear nominal model of the system is established with uncertainty description. This linear model is a multi-function model; it is very helpful for the controller definition, for delay compensation (Smith predictor) or for a feedforward control (if the model is invertible). First the global behavior of the system has to be described, second the system has to be linearized and identified. This linear nominal model is also obtained by minimizing a weighted frequencies criterion. The method has been illustrated on an explicit example: engine torque control. A very precise model is used to simulate the engine. Two controllers have been defined: a PI controller and an LQ controller. In the case of the PI controller two different schemes have been studied: only feedback control or feedforward and feedback control.
ARTICLE IN PRESS 428
Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429
It must be underlined that all controllers have been determined with and for the linear nominal model. It is exactly the same controller that has been used with the non-linear system. The best results have been obtained for a PI controller associated to a Smith predictor. In the case of the feedforward control using the nominal inverse model, the non-linearity is too large (the nominal model is in some cases too distant from the real system). If the family set of linear models is cut as 2 or 3 subsets, perhaps the results could be improved. The LQ controller is not very explicit in this case (of course) because it will operate better in a multi-input configuration (Throttle, EGR, etc.) or in the case of signals with high-frequency noise. It can be noted that it will be the case on real engine. In all the cases, robust stability is guaranteed in the whole range of engine torque control. In closed loop small transient time and very low overshoot are achieved. In future work, the same method will be used on a real engine to prove its effectiveness. In this case multi-inputs must be considered like EGR or turbocharger. To analyze more finely the effectiveness of this method, it should be compared, on the same application, to a robust control design or to a non-linear control like sliding mode control.
Appendix A. The engine model This model represents a one cylinder, half-liter engine. All dimensions are those of a J4S Renault research engine (Table 5). Initial model validations have been performed by comparison with test bed results. The model which has been developed is a quasidimensional model (Heywood, 1988) where volumes (tank, intake, etc.) are connected together via flow singularities (valve, throttle, etc.). Due to the nature of the model, pipe acoustics as well as gas inertia have been
Table 5 Specifications of the engine Number of cylinders Cycle Number of intake valves Number of exhaust valves Bore Stroke Combustion chamber design Displacement Connecting rod length Compression ratio Piston Intake valve O/C Exhaust valve O/C
1 4 stroke 2 2 88.00 mm (3.45 in) 82.00 mm (3.23 in) Pentroof 498.7 cm3 136.9 mm (5.39 in) 9.5:1 Plate 37 BTDC/37 ABDC 37 BBDC/37 ATDC
neglected (Heywood, 1988). The model has been split into several components: valve lift (intake, exhaust, charging valve), pressure drop (intake, exhaust, charging valve), cylinder (variable volume vessel), air tank (constant volume vessel). A.1. Valve lift The valve lift determines the effective section of the valve flow restriction. It could represent a conventional camshaft (lift as a function of crankshaft angle) or a camless valve actuator (lift as a result of the valvetrain force balance). A.2. Pressure drop The pressure drop is calculated with Barre! de Saint Venant equation. Four cases must be distinguished. The flow can be in the forward or reverse direction (back flow) and it can be sub or supersonic. Eq (1) presents the critical ratio which determines the type of the flow. If the pressure drop is lower than the critical ratio, the flow is supersonic, otherwise the flow is subsonic.
2 gþ1
g=g1 :
ðA:1Þ
Eq. (19) represents the mass flow in the supersonic case. Eq. (A.1) represents the mass flow in the subsonic case. Cd is the discharge coefficient and Ar is the area of the singularity. This area is calculated from the valve lift. gþ1=2ðgþ1Þ dm Pu 2 ¼ Cd Ar pffiffiffiffiffiffiffiffi g0:5 ; dt gþ1 rTu ! dm Pu Pd 1=g ¼ Cd Ar pffiffiffiffiffiffiffi dt Pd rTu vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u g1=g ! u 2g Pd : t 1 g1 Pd
ðA:2Þ
ðA:3Þ
A.3. Cylinder The cylinder is the central part of the system. It exchanges work, heat and mass with outside. The first law of thermodynamics can by written as follows: dU ¼ dW þ dQ þ hi dmi þ ho dmo ;
ðA:4Þ
ARTICLE IN PRESS Y. Chamaillard et al. / Control Engineering Practice 12 (2004) 417–429 1e5 0
Open
200
Close
296.5 Area
In take Valve P in P out T in T out Area
Intake Valve Area Open
90
Close
Area
Air Tank Valve Area Open
-45
Close
1e 5 29 6.5
Area
P in P out T int T out Area
Exhaust Valve Area
dm
P
dU V Cylinder
Volume
10 Advance 1.4
O
m
T
Cylinder Volume
A
dm in dm out dh in dh out Sonic ?
Exhaust Valve
Angla
Angle
dm in dm out dh in dh out Sonic ?
Air Tank Valve
Angle
-250
dm in dm out dh int dh out Sonic ?
P in P out T in T out Area
Angle
0
429
A Av m lamda O
dQ
Vibé 2 phases
Crankshaft Angle
N 1200
Engine Speed
Fig. 27. Engine model.
where U is the internal energy, W is the work generated by the moving walls, Q is the heat, m the mass and h the enthalpy. i indicates the intake flow, o indicates the outlet flow. For an ideal gas, the internal energy is dU ¼ mcv dT þ cv T dm
from the input mass and enthalpy fluxes. Conversely, the valve blocks compute enthalpy and mass fluxes from upstream and downstream pressure and temperature. The heat and mass balance are performed by the two sum blocks.
ðA:5Þ
with Eqs. (A.2) and (A.3) the following expression can be obtained: dW þ dQ þ dhi þ dho cv T dm dT ¼ : ðA:6Þ mcv The cylinder pressure is calculated with the ideal gas law (see Eq. (A.5)). The volume is determined with the engine geometry. The mass inside the cylinder is determined by integrating the mass flow through each valve and the temperature by integrating Eq. (A.4). mrT P¼ : ðA:7Þ V A.4. Global model The global model has been realized (Higelin & Charlet, 2001) with Matlab/Simulink. Fig. 27 shows a schematic view of the model. It can be seen that the volume blocks compute pressure, temperature and mass
References Eriksson, L. (1999). Spark advance modeling and control. Thesis of the . Linkoping University. Guzzella, L. (2000). Model-based and robust control systems development. Control System Workshop, Austin, Texas. Heywood, J. B. (1988). Internal combustion engine fundamentals. New York: McGraw-Hill. Higelin, P., & Charlet, A. (2001). Thermodynamic cycles for a new hybrid pneumatic—Combustion engine concept, Fifth international conference on internal combustion engines (ICE 2001), Naples. Jaine, T., Chamaillard, Y., Charlet, A., & Higelin, P. (2002). High frequency IMEP estimation, filtering for torque based SI engine control, SAE paper 2002-01-1276. Kawamura, Y., Shinshi, M., Sato, H., Takahashi, N., & Iriyama, M. (1988). MBT control through individual cylinder pressure detection, SAE paper 881779. Labreuche, G., Da Costa, A., Chamaillard, Y., Charlet, A., Higelin, P., & Perrier, C. (2001). Total friction effective pressure and torque estimator, MECA’O1, Italy. Pestana, G. W. (1989). Engine control methods using combustion pressure feedback, SAE paper 890758.