Abstract:
The present invention relates to a method of estimating the fuel/air ratio in each cylinder of an injection internal-combustion engine comprising an exhaust circuit on which a detector measures the fuel/air ratio of the exhaust gas. An estimator based on a Kalman filter is coupled with a physical model representing the expulsion of the gases from the cylinders and their travel in the exhaust circuit to the detector. The method has application to engine controls.

Description:
CROSS REFERENCE TO RELATED APPLICATION 
     Reference is made to commonly assigned related application Serial No. U.S. Ser. No. 11/437,702, filed in the U.S. Patent and Trademark Office on May 22, 2006, entitled “Method of Estimating the Fuel/Air Ratio in a Cylinder of an Internal Combustion Engine by Means of an Adaptive NonLinear Filter 
     BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The present invention relates to a method of estimating the fuel/air ratio of each cylinder of a fuel injected internal-combustion engine from a fuel/air ratio measurement downstream from the manifold and an extended Kalman filter. 
     2. Description of the Prior Art 
     Knowledge of the fuel/air ratio, characterized by the ratio of the mass of fuel to the mass of air, is important for all vehicles, whether equipped with gasoline engines, since it conditions good combustion of the mixture when it is close to 1, or with diesel engines, in which case the interest of knowing the fuel/air ratio is different insofar as diesel engines run under lean mixture conditions (ratio below 1). In particular, catalysts using a NOx trap lose efficiency in the course of time. In order to recover optimum efficiency, the fuel/air ratio has to be maintained close to 1 for some seconds, prior to returning to normal running conditions with a lean mixture. Depollution by DeNox catalysis therefore requires precise control of the fuel/air ratio cylinder by cylinder. 
     A probe arranged at the turbine outlet (supercharged engine) and upstream from the NOx trap therefore gives a measurement of the mean fuel/air ratio as a result of the exhaust process. This measurement being highly filtered and noise-affected is used for control of the masses injected into the cylinders during fuel/air ratio phases of 1, each cylinder receiving then the same mass of fuel. 
     In order to control more precisely, and in particular individually, injection of the fuel masses into the cylinders, reconstruction of the fuel/air ratio in each cylinder is necessary. Since installing fuel/air ratio probes at the outlet of each cylinder of a vehicle cannot be done considering their cost price, setting an estimator working from the measurements provided by a single probe advantageously allows to separately know the fuel/air ratios in each cylinder. 
     An engine control can thus, from the reconstructed fuel/air ratios, adjust the fuel masses injected into each cylinder so that the fuel/air ratios are balanced in all the cylinders. 
     French Patent 2,834,314 describes the definition of a model, based upon observation and filtering by means of a Kalman filter. This model contains no physical description of the mixture in the manifold and does not take into account highly pulsating flow rate phenomena. 
     Estimation of the fuel/air ratio in the cylinders is only conditioned by the coefficients of a matrix, coefficients that can be identified off-line by means of an optimization algorithm. Furthermore, a different adjustment of the matrix, therefore an identification of its parameters, corresponds to each working point (engine speed/load). This estimator thus requires heavy acquisition test means (with 5 fuel/air ratio probes) and has no robustness in case of engine change. 
     SUMMARY OF THE INVENTION 
     The present invention allows finer modelling of the exhaust process so as to, on the one hand, do without the identification stage and, on the other hand, provide the fuel/air ratio estimation model with more robustness, for all the engine working points. The invention furthermore allows performing of measurements every 6° of crankshaft rotation, and therefore to obtain high-frequency fuel/air ratio measurement information without being disturbed by the measurement noise. 
     The present invention thus relates to a method of estimating the fuel/air ratio in each cylinder of an internal-combustion engine comprising a gas exhaust circuit including at least cylinders connected to a manifold and a fuel/air ratio detector (λ) downstream from the manifold. The method is characterized in that it comprises the following steps: 
     defining an estimation of the fuel/air ratio (λ) measured by the detector from at least one variable of the model; 
     performing a modelling of the transfer function of said detector wherein the measured fuel/air ratio estimation is taken into account; 
     establishing a physical model representing in real time the expulsion of the gases from each one of the cylinders and their travel in the exhaust circuit up to the detector, wherein the modelling of the transfer function is taken into account; 
     coupling said model with an extended Kalman type nonlinear estimator; and 
     performing a real-time estimation of the fuel/air ratio value in each cylinder from the extended Kalman type nonlinear estimator. 
     According to the method, the transfer function can be modelled from a first order filter. 
     A lag time due to the gas transit time and to the detector response time can also be evaluated by carrying out a test disturbance in a determined cylinder and by measuring its effect on the detector. 
     According to an embodiment, the physical model can comprise at least the following four variables: the total mass of gas in the exhaust manifold (M T ), the mass of fresh air in the exhaust manifold (M air ), the fuel/air ratio measured by said detector (λ) and the fuel/air ratios in each cylinder (λ i ). This embodiment can also include at least the following two output data: the total mass of gas in the exhaust manifold (M T ) and the mass flow rates leaving the cylinders (d i ). 
     The measured fuel/air ratio (λ) can be estimated as a function of the total mass of gas in the exhaust manifold (M T ) and the mass of fresh air in the exhaust manifold (M air ). Estimation of the fuel/air ratio value in each cylinder can then involve real-time correction of the estimation of the total mass of gas in the exhaust manifold (M T ), of the estimation of the mass of fresh air in the exhaust manifold (M air ) and of the estimation of the fuel/air ratio value in each cylinder (λ i ). 
     Finally, the method can be applied to an engine control for adapting the fuel masses injected into each cylinder in order to adjust the fuel/air ratio in all the cylinders. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Other features and advantages of the present invention will be clear from reading the description hereafter of a non-limitative embodiment example, illustrated by the accompanying figures, wherein: 
         FIG. 1  diagrammatically shows the descriptive elements of the exhaust process; 
         FIG. 2  illustrates the reference fuel/air ratios (λ i   ref ) as a function of time (T) and the results of the estimator according to the invention ({circumflex over (λ)} i ) as a function of time (T), for each one of the four cylinders; 
         FIG. 3  shows the structure of the estimator; 
         FIGS. 4A to 4C  illustrate the reference fuel/air ratios (λ i   ref ) as a function of time (T) and the results of the estimator taking account of the lag according to the invention ({circumflex over (λ)} i ) on a first working point at 2000 rpm at high load (9 bar) as a function of time (T), for each one of the four cylinders; and 
         FIGS. 4D to 4F  illustrate the reference fuel/air ratios (λ i   ref ) as a function of time (T) and the results of the estimator taking account of the lag according to the invention ({circumflex over (λ)} i ) on a second working point at 2500 rpm at low load (3 bar) as a function of time (T), for each one of the four cylinders. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Description of the Exhaust Process 
     The exhaust process comprises the path traveled by the gases from the exhaust valve to the open air, at the exhaust muffler outlet. The engine in the present embodiment example is a 2200-cm 3  4-cylinder engine. It is equipped with a variable-geometry turbo/supercharger. The diagram of  FIG. 1  shows the descriptive elements of the exhaust process, wherein: 
     λ 1  to λ 4  are the fuel/air ratios in each one of the four cylinders, 
     SR is the fuel/air ratio probe, 
     CE corresponds to the exhaust manifold, 
     T corresponds to the turbine of the turbo/supercharger, 
     DS 1  to DS 4  represent the flow rates at the cylinder outlets. 
     Fuel/air ratio probe (SR) is arranged just after turbine (T). The gases, after combustion in the cylinder, undergo the following actions: 
     passage through the exhaust valve. The latter is controlled by a camshaft, with the lift law being bell-shaped. The flow rates will go from a high value, when the valve opens, to a lower value when the cylinder and manifold pressures become equal, and they will eventually increase again when the piston starts to slide up again to expel the exhaust gases; 
     passage through a short pipe connecting the manifold to the cylinder head outlet; 
     a mixing phase in exhaust manifold (CE) where the flows (DS 1  to DS 4 ) from the four cylinders meet. It is here that the individual cylinder outputs of exhaust mix, depending on the manifold type (symmetrical or asymmetrical), on the EEO (Early Exhaust Opening) and on the LEC (Late Exhaust Closing) which will determine the flow overlap proportion; 
     passage through the turbine which supplies the compressor arranged upstream from the intake with the required torque. Although its action on the flow rates is not well known, one may consider that it is going to mix even more of the individual cylinder outputs coming from the various cylinders; and measurement by the UEGO type probe. 
     The composition of the exhaust gases depends on the amounts of fuel and of air fed into the combustion chamber, on the composition of the fuel and on the development of the combustion. 
     In practice, the fuel/air ratio probe measures the O 2  concentration inside a diffusion chamber connected to the exhaust pipe by a diffusion barrier made of a porous material. This configuration can induce differences depending on the location selected for the probe, notably because of the temperature and/or pressure variations near the fuel/air ratio probe. 
     This fuel/air ratio variation phenomenon depending on the pressure or on the temperature has however been disregarded since what is sought is detection of fuel/air ratio disparities between the cylinders, the mean value being normally kept by the estimator. 
     In the real-time physical model used by the estimator according to the invention, the measured fuel/air ratio (A) is related to the mass of air (or flow of air) around the probe and to the total mass (or total flow rate). The model is based on a three-gas approach: air, fuel and burnt gases. It is thus considered that, with a lean mixture, all of the gas remaining after combustion is a mixture of air and of burnt gases. For a rich mixture, the fuel being in excess, unburnt fuel and burnt gases are present after combustion, whereas all of the air has disappeared. In reality, the combustion is never 100% complete, but for the estimator it is considered to be complete. 
     A formulation relating the fuel/air ratio to the masses of the three species mentioned is defined. In the case of a lean mixture: the air is in excess, and no fuel is left after combustion. Before combustion, the following masses are assumed to be present in the cylinder: M air  the mass of air, M carb  the mass of fuel and M gazB  the mass of burnt gases:
         M air =; M carb =Y; M gazB =0       

     Knowing that 14.7 times as much air as fuel is required to reach stoichiometric conditions, the table hereunder giving the masses of each species before and after combustion can be drawn up: 
     
       
         
               
               
               
               
             
               
               
               
               
             
           
               
                   
                   
               
               
                   
                   
                   
                 Mass of burnt 
               
               
                   
                 Mass of air 
                 Mass of fuel 
                 gas 
               
               
                   
                   
               
             
             
               
                   
               
             
          
           
               
                 Before combustion 
                 x 
                 y 
                 0 
               
               
                 After combustion 
                 x − 14.7 × y 
                 0 
                 y + 14.7 × y 
               
               
                   
               
             
          
         
       
     
     The fuel/air ratio λ representing ratio 
                 M   carb     /     M   air       ,         
after calculation the following formulation is obtained, valid only if the mixture is a lean mixture:
 
             λ   =         M   gazB     ·   PCO           M   air     ·     (     1   +   PCO     )       +       M   gazB     ·   PCO               
where PCO corresponds to ratio
 
               M   air     /     M     carb   ⁢                     
when the mixture is stoichiometric. PCO is the calorific value of the fuel.
 
     For a rich mixture, the formula is as follows 
     
       
         
           
             λ 
             ⁢ 
             
                 
             
             = 
             
               
                 
                   
                     M 
                     carb 
                   
                   · 
                   
                     ( 
                     
                       1 
                       + 
                       PCO 
                     
                     ) 
                   
                 
                 
                   M 
                   gazB 
                 
               
               + 
               1 
             
           
         
       
     
     However, these formulas are valid in the case where the mixture contains no EGR since the presence of burnt gases at the intake changes the concentrations of the three gases at the exhaust. 
     In the present embodiment, only the fuel/air ratio formula for lean mixtures is used in the estimator, for integration of the fuel/air ratio in equation (7), a very small part of the air (&lt;3%) being disregarded. However, the invention is not limited to this embodiment; in fact, the formula is continuous in the vicinity of a fuel/air ratio of 1, and its inversion poses no problems for rich mixtures. 
     In order to better apprehend the way the gases mix in the exhaust pipes, a diesel engine model was used with the AMESim software of the IMAGINE Company (France). This model, which cannot be inverted, will be used as a reference to validate the model according to the invention. 
     AMESim is a 0D modelling software, particularly well-suited for thermal and hydraulic phenomena. It notably allows to model volumes, pipes or restrictions. P The exhaust model comprises: 
     the exhaust pipes represented by a volume and a tube, 
     the exhaust manifold with thermal exchanges, 
     the turbine and the bypass valve, 
     a volume at the confluence of the turbine and valve flow rates, 
     a tube between the turbine and the measuring probe, 
     a volume and a tube for the exhaust line. 
     The elementary blocks for modelling the pipes, restrictions and volumes are described in the AMESim instruction manual “Thermal Pneumatic Library”. The standard equations are used to calculate a flow rate through a restriction and the mass and energy conservation. Furthermore, the model takes account of the inertias of the gases, which is important to study the gas composition dynamics. 
     Since it is a 0D model, the time dimension is not taken into account, and it is not possible to model a lag time with a physical approach. If an input variable is changed, the output is immediately changed. The transport time is thus disregarded. This limitation is important when trying to work on real-time acquisitions. 
     According to the invention, a single real-time physical model is defined to model the global system, that is the entire path traveled by the exhaust gases, from the cylinders through the manifold up to the exhaust downstream from the turbine. 
     I—Definition of a Real-Time Physical Model 
     In the present embodiment, the temperature variation is considered low over an engine cycle, and that its action is limited on the flow rate variations. The pressure variations are in fact essential in the process since they are directly related to the flow rates. A fixed temperature is thus set for each element: cylinders, manifold and turbine. The heat exchanges are therefore not modelled either. This simplification hypothesis does not have much impact. 
     In a first approach, two gases are considered: fresh air and burnt gases. The conventional equations describe the evolution of the total mass of the gases in the volumes, and of the mass of fresh air. The burnt gases can then be deduced therefrom. This procedure is valid in the case of lean mixture conditions, but similar equations can be written for the fuel and the burnt gases, in the case of a rich mixture. 
     A) Physical Model of the Exhaust Manifold 
     The exhaust manifold is modelled according to a volume in which there is mass conservation. The temperature is assumed to be substantially constant and determined from a chart as a function of the engine speed and load. 
     According to the invention, it has been chosen to relate the measured fuel/air ratio to the mass of air around the probe and to the total mass. Thus, conservation of the total mass in the manifold expresses the fact that the exhaust gas mass in the manifold is equal to the exhaust gas mass entering the manifold (cylinder outlet flow rate) decreased by the mass leaving the manifold. The composition of the flow in the turbine is assumed to be the same as at the manifold outlet. Thus, the mass leaving the manifold is equal to the flow passing through the turbine. The following formula is for the total mass: 
                       N   e     ⁢       ⅆ     M   T         ⅆ   α         =         ∑     i   =   1       n   cyl       ⁢       d   i     ⁡     (   α   )         -       d   T     ⁡     (     M   T     )                 (   1   )               
with:
 
     N e : engine speed 
     α: crankshaft angle 
     M T : total mass in the exhaust manifold 
     d i : mass flow rate leaving cylinder i 
     d T : total flow rate passing through the turbine. 
     Similarly, for the air mass conservation, the relationship is: 
                       N   e     ⁢       ⅆ     M   air         ⅆ   α         =         ∑     i   =   1       n   cyl       ⁢       (     1   -     λ   i       )     ·       d   i     ⁡     (   α   )           -       d   air     ⁡     (     M   air     )                 (   2   )               
with:
 
     N e : engine speed 
     α: crankshaft angle 
     M air : fresh air mass in the exhaust manifold 
     λ i : fuel/air ratio in each cylinder 
     d i : mass flow rate leaving cylinder i 
     d air : air flow rate passing through the turbine. 
     The physical models are described allowing determination of the flow rate at the cylinder outlet and the flow rates passing through the turbine. 
     Model Allowing Determination of the Flow Rate at the Cylinder Outlet: Gas Expulsion 
     The gas flow rate at the cylinder outlet can be modelled by means of a physical model describing the flow rate at the outlet of the exhaust valves. Three variables are used for this expulsion model of the gases through the valves: 
     crankshaft angle (α), 
     the flow inducted through the cylinder d asp  (variable estimated by the upstream engine control), 
     the mean fuel/air ratio value measured by the probe  λ  over a cycle. 
     The mean outgoing flow is known from the intake flow and from the injected gasoline flow rate. The instantaneous value of the outgoing flow is based on a template depending on the sucked flow. This template is a physical model (curve) based on an empirical law allowing estimation of a mean flow rate for a cylinder as a function of the crankshaft angle from the engine speed, the crankshaft angle, the intake flow by the cylinder and the mean fuel/air ratio value being measured by the probe over a cycle. The only constraint of this physical law is to respect the mean outgoing flow (curve area) and to provide a curve accounting for the two phenomena as follows: 
     the cylinder/exhaust pressure balance expressed by a flow rate peak as a function of the crankshaft angle; 
     a flow rate that depends on the cross section of flow of the exhaust valve, expressed by a second flow rate peak of lower amplitude. 
     This template (  d ) provides at the output a course (curve) of the mass flow rate at the exhaust valves outlet d i , that is a common estimation (  d ) of the flow rate for all the cylinders. It is obtained in correlation with the engine test bench measurements. Depending on the intake flow and on the mean fuel/air ratio measured by the probe  λ , then a homothetic transformation of the template and a phase shift for each cylinder as a function of the crankshaft angle is carried out, so as to deduce the flow rate of the gases at the outlet of each cylinder:
 
 d   i (α)=   d   (α 1 +α)·α 0  
 
with:
 
     d i (α) the gas outlet flow rate at the outlet of cylinder i, 
       d (α) the template, that is an estimation of the flow rate at the cylinder outlet, 
     
       
         
           
             
               α 
               0 
             
             = 
             
               
                 d 
                 asp 
               
               · 
               
                 ( 
                 
                   1 
                   + 
                   
                     
                       λ 
                       _ 
                     
                     PCO 
                   
                 
                 ) 
               
             
           
         
       
     
     α i  the phase shift angle for cylinder i. 
     The phase shift of the template curve can be diagrammatically seen in  FIG. 1  (DS 1  to DS 4 ). 
     The physical models allowing determination of the flow rate passing through the turbine are described hereafter. 
     Model Allowing Determination of the Flow Rate Passing Through the Turbine: Turbine Model 
     The turbine is modelled according to a flow passing through a flow rate restriction. The flow rate in the turbine is generally given by mapping (chart) as a function of the turbine speed and of the pressure ratio upstream/downstream from the turbine. 
     The flow rate passing through the turbine d T  is a function of the total mass (M T ) in the exhaust manifold, of the temperature in the exhaust manifold, of the turbo/supercharger speed and of the turbo/supercharger geometry. The input data of this model thus are: 
     the total mass of exhaust gas (M T ), 
     the mass of air (M air ), 
     the engine speed (N e ), and 
     the (turbo/supercharger) turbine speed. 
     This flow rate can be estimated from a concave function of the total mass M T . This function is denoted by p. The flow rate in the turbine is then written as follows: d T (M T )=M T ·p(M T ). 
     Function p is a root type function that is expressed as a function of the turbine speed on the one hand and of the ratio of the total mass in the exhaust manifold (M T ) to the mass in the manifold under atmospheric conditions (M 0 ) on the other hand. Thus, the mapping gives p(M T ) as a function of ratio 
               M   T       M   0           
and of the (turbo/supercharger) turbine speed. The formula used by this mapping is:
 
               p   ⁡     (     M   T     )       =       f   ⁡     (     turbine   ⁢             ⁢             ⁢   speed     )       ·         2   ·   g       g   -   1         ·     (         (       M   T       M   0       )         -   2     g       -       (       M   T       M   0       )         -     (     g   +   1     )       g         )             
where:
 
     ƒ is a polynomial function, 
     g is a constant. 
     The parameters of function ƒ are optimized by correlation with the turbine mapping. 
     Furthermore, the air composition is assumed to be the same as in the exhaust manifold. The flow of air passing through the turbine thus is: 
     
       
         
           
             
               
                 d 
                 air 
               
               ⁡ 
               
                 ( 
                 
                   M 
                   air 
                 
                 ) 
               
             
             = 
             
               
                 
                   
                     M 
                     air 
                   
                   
                     M 
                     T 
                   
                 
                 ⁢ 
                 
                   
                     d 
                     T 
                   
                   ⁡ 
                   
                     ( 
                     
                       M 
                       T 
                     
                     ) 
                   
                 
               
               = 
               
                 
                   M 
                   air 
                 
                 · 
                 
                   p 
                   ⁡ 
                   
                     ( 
                     
                       M 
                       T 
                     
                     ) 
                   
                 
               
             
           
         
       
     
     Thus, by means of the physical gas expulsion and turbine models, equations (1) and (2) are written as follows: 
     
       
         
           
             
               
                 
                   { 
                   
                     
                       
                         
                           
                             
                               N 
                               e 
                             
                             ⁢ 
                             
                               
                                 ⅆ 
                                 
                                   M 
                                   T 
                                 
                               
                               
                                 ⅆ 
                                 α 
                               
                             
                           
                           = 
                           
                             
                               
                                 ∑ 
                                 
                                   i 
                                   = 
                                   1 
                                 
                                 
                                   n 
                                   cyl 
                                 
                               
                               ⁢ 
                               
                                 
                                   d 
                                   i 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   α 
                                   ) 
                                 
                               
                             
                             - 
                             
                               
                                 M 
                                 T 
                               
                               · 
                               
                                 p 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     M 
                                     T 
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                     
                     
                       
                         
                           
                             
                               N 
                               e 
                             
                             ⁢ 
                             
                               
                                 ⅆ 
                                 
                                   M 
                                   air 
                                 
                               
                               
                                 ⅆ 
                                 α 
                               
                             
                           
                           = 
                           
                             
                               
                                 ∑ 
                                 
                                   i 
                                   = 
                                   1 
                                 
                                 
                                   n 
                                   cyl 
                                 
                               
                               ⁢ 
                               
                                 
                                   ( 
                                   
                                     1 
                                     - 
                                     
                                       λ 
                                       i 
                                     
                                   
                                   ) 
                                 
                                 · 
                                 
                                   
                                     d 
                                     i 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     α 
                                     ) 
                                   
                                 
                               
                             
                             - 
                             
                               
                                 M 
                                 air 
                               
                               · 
                               
                                 p 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     M 
                                     T 
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     This system of equations (3) is the physical model of the exhaust manifold. 
     The input data of this model are: 
     N e : engine speed 
     α: crankshaft angle 
     d i : mass flow rate leaving cylinder i 
     d T : total flow rate passing through the turbine 
     d air : air flow rate passing through the turbine 
     λ i : fuel/air ratio in each cylinder and the unknowns of the system are: 
     M T : total mass in the exhaust manifold 
     M air : mass of fresh air in the exhaust manifold. 
     The first equation contains one unknown: M T . The second one contains two unknowns: M air  and λ i . This leads to the additional hypotheses described hereafter. 
     B) Hypothesis on the Cylinder Outlet Fuel/Air Ratios 
     In order to complete the real-time physical model (RTM) of the exhaust manifold, it is assumed that the fuel/air ratios at the cylinder outlet are constant over a working point, therefore: 
     
       
         
           
             
               
                 
                   
                     
                       N 
                       e 
                     
                     · 
                     
                       
                         ⅆ 
                         
                           λ 
                           i 
                         
                       
                       
                         ⅆ 
                         α 
                       
                     
                   
                   = 
                   0 
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
     In fact, since calculation is carried out in real time, constants λ i  are estimated. 
     C) Hypothesis on the Detector Dynamics 
     In the previous equations, the fuel/air ratio measured at the detector is calculated from the fuel/air ratio in the cylinders, the flow of air at the cylinder outlet and the total flow of gas. This structure is difficult to use in a Kalman filter because the inputs of the model have to be estimated. The model is therefore completed by addition of the inputs M T  and λ (Mohinder S. Grewal: “Kalman Filtering Theory and Practice”, Prentice Hall, 1993). 
     Therefore the response dynamics of the fuel/air ratio detector is taken into account and the transfer function of the measuring probe (of UEGO type) is modelled according to a first order filter. The fuel/air ratio (λ) given by the model downstream from the turbine is equal to the fuel/air ratio in the manifold and, by using the fuel/air ratio equation described above, it is possible to estimate the measured fuel/air ratio from M T  and M air : 
     
       
         
           
             
               
                 
                   λ 
                   = 
                   
                     
                       
                         
                           M 
                           gazB 
                         
                         · 
                         PCO 
                       
                       
                         
                           
                             M 
                             air 
                           
                           · 
                           
                             ( 
                             
                               1 
                               + 
                               PCO 
                             
                             ) 
                           
                         
                         + 
                         
                           
                             M 
                             gazB 
                           
                           · 
                           PCO 
                         
                       
                     
                     = 
                     
                       1 
                       - 
                       
                         
                           M 
                           air 
                         
                         
                           M 
                           T 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     The probe response can thus be modelled by the following relation: 
                       N   e     ⁢       ⅆ   λ       ⅆ   α         =       1   τ     ⁢     (           M   T     -     M   air         M   T       -   λ     )               (   5   )               
with: τtime constant of the fuel/air ratio detector considered to be of first order.
 
     D) Expression of the Real-Time Physical Model 
     Finally, the real-time physical model RTM can be expressed in the matrical form as follows from equations (3), (5) and (6): 
     
       
         
           
             
               
                 
                   X 
                   = 
                   
                     
                       
                         [ 
                         
                           
                             
                               
                                 M 
                                 T 
                               
                             
                           
                           
                             
                               
                                 M 
                                 air 
                               
                             
                           
                           
                             
                               
                                 λ 
                                 1 
                               
                             
                           
                           
                             
                               
                                 λ 
                                 2 
                               
                             
                           
                           
                             
                               
                                 λ 
                                 3 
                               
                             
                           
                           
                             
                               
                                 λ 
                                 4 
                               
                             
                           
                         
                         ] 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           N 
                           e 
                         
                         · 
                         
                           
                             ⅆ 
                             X 
                           
                           
                             ⅆ 
                             α 
                           
                         
                       
                     
                     = 
                     
                       [ 
                       
                         
                           
                             
                               
                                 
                                   ∑ 
                                   
                                     i 
                                     = 
                                     1 
                                   
                                   
                                     n 
                                     cyl 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     d 
                                     i 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     α 
                                     ) 
                                   
                                 
                               
                               - 
                               
                                 
                                   M 
                                   T 
                                 
                                 · 
                                 
                                   p 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       M 
                                       T 
                                     
                                     ) 
                                   
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 
                                   ∑ 
                                   
                                     i 
                                     = 
                                     1 
                                   
                                   
                                     n 
                                     cyl 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         λ 
                                         i 
                                       
                                     
                                     ) 
                                   
                                   · 
                                   
                                     
                                       d 
                                       i 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       α 
                                       ) 
                                     
                                   
                                 
                               
                               - 
                               
                                 
                                   M 
                                   air 
                                 
                                 · 
                                 
                                   p 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       M 
                                       T 
                                     
                                     ) 
                                   
                                 
                               
                             
                           
                         
                         
                           
                             0 
                           
                         
                         
                           
                             0 
                           
                         
                         
                           
                             0 
                           
                         
                         
                           
                             0 
                           
                         
                       
                       ] 
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     The unknowns of the physical model are eventually M T , M air , λ and the □ i . 
     The output data of the physical model are M T  and d i . 
     E) Exhaust Lag Time 
     The lag times due to the transportation of gas in the pipes and the various volumes, and the “idle time” of the measuring probe, are not taken into account in the physical model described above (system of equations 7). However, the model is constructed linearly in relation to these lag times because transportation in the pipes is disregarded. They can therefore be compiled into a single lag time for all of the exhaust process, and the physical model can be inverted as it is, since the influence of the lag time can be considered later, as explained hereafter. 
     II—Fuel/Air Ratio Estimator 
     The above physical model (7) describes that the fuel/air ratio downstream from the turbine (considered to be identical to the fuel/air ratio in the manifold) is expressed as a function of the composition of the gas flow at the exhaust manifold inlet. 
     The measured data are: 
     Fuel/air ratio measured by the probe: λ 
     The other known data of the system are: 
     Engine speed: N e    
     Crankshaft angle: α 
     (Supercharger) turbine speed 
     Intake flow inducted by the cylinder. 
     The modelled data of the system are: 
     Mass flow rate coming from cylinder i: d i    
     Total flow rate passing through the turbine: d T    
     Flow of air passing through the turbine: d air    
     Total mass in the exhaust manifold: M T . 
     The unknowns thus are: 
     Fuel/air ratio in each one of the four cylinders: λ i    
     Mass of fresh air in the exhaust manifold: M air . 
     Physical model (7) is nonlinear, and such a system cannot be solved in real time. It is therefore necessary to use an estimator rather than try to directly calculate the unknowns of the system. Selection of the estimator according to the invention is based on the fact that this model has a structure that can be used in an extended Kalman filter. Thus, in order to estimate the unknowns from the physical model RTM, the method according to the invention uses an estimator based on an extended Kalman filter. Such a filter is well known in the art and it is described in the following document: 
     Greg Welch and Gary Bishop: “An Introduction to the Kalman Filter”, University of North Carolina—Chapel Hill TR95-041. May 23, 2003. 
     The structure of an extended Kalman filter is described hereafter. The extended Kalman filter (EKF) allows estimation of the state vector of a process in cases where the latter, or the measuring process, is nonlinear. 
     It is assumed that a process (x) is governed by a nonlinear stochastic equation (f). At the time increment k, it can be written: 
     
       
         
           
             
               x 
               k 
             
             = 
             
               
                 
                   f 
                   ⁡ 
                   
                     ( 
                     
                       
                         x 
                         
                           k 
                           - 
                           1 
                         
                       
                       , 
                       
                         w 
                         
                           k 
                           - 
                           1 
                         
                       
                     
                     ) 
                   
                 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 with 
                 ⁢ 
                 
                   : 
                 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   x 
                   k 
                 
               
               = 
               
                 [ 
                 
                   
                     
                       
                         M 
                         
                           T 
                           ⁡ 
                           
                             ( 
                             k 
                             ) 
                           
                         
                       
                     
                   
                   
                     
                       
                         M 
                         
                           air 
                           ⁡ 
                           
                             ( 
                             k 
                             ) 
                           
                         
                       
                     
                   
                   
                     
                       
                         λ 
                         
                           ( 
                           k 
                           ) 
                         
                       
                     
                   
                   
                     
                       
                         λ 
                         
                           1 
                           ⁢ 
                           
                             ( 
                             k 
                             ) 
                           
                         
                       
                     
                   
                   
                     
                       
                         λ 
                         
                           2 
                           ⁢ 
                           
                             ( 
                             k 
                             ) 
                           
                         
                       
                     
                   
                   
                     
                       
                         λ 
                         
                           3 
                           ⁢ 
                           
                             ( 
                             k 
                             ) 
                           
                         
                       
                     
                   
                   
                     
                       
                         λ 
                         
                           4 
                           ⁢ 
                           
                             ( 
                             k 
                             ) 
                           
                         
                       
                     
                   
                 
                 ] 
               
             
           
         
       
     
     A measurement (y) is given by a nonlinear observation equation h. At the time increment k, it can be written:
 
 y   k   =h ( x   k   ,v   k )
 
where the random variables w k  and v k  respectively represent the model noises and the measurement noises.
 
     The estimation method comprises an estimator based on a prediction/correction technique, that is a prediction of the variable is performed, then a correction is applied thereto at each time increment. These prediction/correction stages are described below within the general context of an extended Kalman filter: 
     It is denoted: 
     by {circumflex over (x)}   k    the prediction of x at time increment k, 
     by {circumflex over (x)} k  the estimation of x at time increment k. 
     Prediction stage: 
     
       
         
           
               
             
               
                 
                   
                     
                       
                         x 
                         ^ 
                       
                       k 
                       - 
                     
                     = 
                     
                       f 
                       ⁡ 
                       
                         ( 
                         
                           
                             
                               x 
                               ^ 
                             
                             
                               k 
                               - 
                               1 
                             
                           
                           , 
                           0 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   
                     
                       P 
                       k 
                       - 
                     
                     = 
                     
                       
                         
                           A 
                           k 
                         
                         ⁢ 
                         
                           P 
                           
                             k 
                             - 
                             1 
                           
                         
                         ⁢ 
                         
                           A 
                           k 
                           T 
                         
                       
                       + 
                       
                         
                           W 
                           k 
                         
                         ⁢ 
                         
                           Q 
                           
                             k 
                             - 
                             1 
                           
                         
                         ⁢ 
                         
                           W 
                           k 
                           T 
                         
                       
                     
                   
                 
               
             
           
         
       
     
     In the above equations, the following variables are used:
         P k  is the error dispersion positive definite symmetric matrix,   Q k  is the Gaussian white noise w k  dispersion positive definite symmetric matrix,   R k  is Gaussian white noise v k  dispersion positive definite symmetric matrix,   P   k    is the error propagation positive definite symmetric matrix.       

     The various elements are initialized by means of the values obtained during simulation with AmeSim. 
     Correction Stage: 
     Estimation of x at time increment k, {circumflex over (x)} k , is obtained by correction of the prediction {circumflex over (x)}   k    as follows: 
                           K   k     =       P   k   -     ⁢         H   k   T     ⁡     (         H   k     ⁢     P   k   -     ⁢     H   k   T       +       V   k     ⁢     R   k     ⁢     V   k   T         )         -   1                         x   ^     k     =         x   ^     k   -     +       K   k     ⁡     (       y   k     -     h   ⁡     (         x   ^     k   -     ,   0     )         )                       P   k     =       (     I   -       K   k     ⁢     H   k         )     ⁢     P   k   -                     
where:
 
     A is the Jacobian matrix of the partial derivatives of f with respect to x: 
     
       
         
           
             
               A 
               
                 [ 
                 
                   i 
                   , 
                   j 
                 
                 ] 
               
             
             = 
             
               
                 
                   ∂ 
                   
                     f 
                     
                       [ 
                       i 
                       ] 
                     
                   
                 
                 
                   ∂ 
                   
                     x 
                     
                       [ 
                       j 
                       ] 
                     
                   
                 
               
               ⁢ 
               
                 ( 
                 
                   
                     
                       x 
                       ^ 
                     
                     k 
                   
                   , 
                   0 
                 
                 ) 
               
             
           
         
       
     
     W is the Jacobian matrix of the partial derivatives of f with respect to w: 
     
       
         
           
             
               W 
               
                 [ 
                 
                   i 
                   , 
                   j 
                 
                 ] 
               
             
             = 
             
               
                 
                   ∂ 
                   
                     f 
                     
                       [ 
                       i 
                       ] 
                     
                   
                 
                 
                   ∂ 
                   
                     w 
                     
                       [ 
                       j 
                       ] 
                     
                   
                 
               
               ⁢ 
               
                 ( 
                 
                   
                     
                       x 
                       ^ 
                     
                     k 
                   
                   , 
                   0 
                 
                 ) 
               
             
           
         
       
     
     H is the Jacobian matrix of the partial derivatives of h with respect to x: 
     
       
         
           
             
               H 
               
                 [ 
                 
                   i 
                   , 
                   j 
                 
                 ] 
               
             
             = 
             
               
                 
                   ∂ 
                   
                     h 
                     
                       [ 
                       i 
                       ] 
                     
                   
                 
                 
                   ∂ 
                   
                     x 
                     
                       [ 
                       j 
                       ] 
                     
                   
                 
               
               ⁢ 
               
                 ( 
                 
                   
                     
                       x 
                       ^ 
                     
                     k 
                     - 
                   
                   , 
                   0 
                 
                 ) 
               
             
           
         
       
     
     V is the Jacobian matrix of the partial derivatives of h with respect to v: 
               V     [     i   ,   j     ]       =         ∂     h     [   i   ]           ∂     x     [   j   ]           ⁢     (         x   ^     k   -     ,   0     )             
It can be noted that, in order to lighten the notations, the index of time increment k is not given, even though these matrices are in fact different at each increment.
 
     At the input of the Kalman filter, the fuel/air ratio λ downstream from the turbine and the total mass of gas in the manifold (M T ) are necessary. The measured or modelled input parameters of the estimator thus are: 
     
       
         
           
             
               
                 
                   Y 
                   = 
                   
                     [ 
                     
                       
                         
                           
                             M 
                             T 
                           
                         
                       
                       
                         
                           λ 
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
           
         
       
     
     The fuel/air ratio is measured and the total mass of gas is the result of the calculation of model (7) in parallel with the Kalman filter. 
     This estimator based on a Kalman filter finally allows estimation of the fuel/air ratio cylinder by cylinder from the fuel/air ratio measurement provided by the detector arranged behind the turbine. 
     The estimator thus constructed allows real-time correction of M T , M air , λ i  and λ from a first value of M T  provided by the real-time model (7) and from the fuel/air ratio measured by the probe. 
     The Kalman filter is numerically solved in real time, the calculator using an explicit Euler discretization known in the art. 
     Simulation Results: Estimator Test (9) 
     From the known individual fuel/air ratios, it is estimated by means of the reference modelling AMESim a fuel/air ratio at the probe (λ). This fuel/air ratio value (λ) is used at the estimator input. The dynamics of the probe has not been taken into account. An injection unbalance is applied and the cylinder by cylinder estimation of fuel/air ratio (λ i ) from the fuel/air ratio measured behind the turbine (λ) is observed. 
     For simulation, the 4 cylinders by introducing 80 μs more injection on the cylinder, then cylinder  1  and  4  are unbalanced similarly.  FIGS. 2A and 2B  show, at the bottom, the fuel/air ratios (λ i   ref ) given by AmeSim as a function of time (T) and, at the top, the results of the estimator ({circumflex over (λ)} i ) as a function of time (T). The four curves correspond to each one of the four cylinders. The performance of the estimator based on the Kalman filter is very good. A slight phase difference, due to the inertia of the gas that is not taken into account in the present model, is however noted. It is therefore decided to complete the model and the estimator by an exhaust lag time estimator. 
     Exhaust Lag Time Estimator 
     The estimator implemented as described above does not allow the estimation method to account for the lag time between the cylinder exhaust and the signal acquired by the probe. In reality, the lag time is due to several sources: transport time in the pipes and through the volumes, idle time of the measuring probe. 
     By applying a lag time D at the estimator input to the variables from the model, the estimator can be synchronized with the fuel/air ratio measurements. The structure of the estimator with a lag time is illustrated in  FIG. 3 , wherein: 
     N e  and α are the input data of the real-time model RTM described by equations (7); 
     MMBO is the Open Loop Mass Model (model RTM); 
     D is the lag time applied to the output variables of model RTM (MMBO); this lag time is obtained from equation (10); 
     SR is the probe measuring the fuel/air ratio downstream from the turbine used; 
     ERFK is the Fuel/Air Ratio Estimator based on a Kalman Filter and described by equation (9); and 
     λ i  is the fuel/air ratio in cylinder i estimated by estimator ERFK. 
     The lag time depends on the operating conditions: engine speed, load, exhaust manifold pressure, etc. Since the lag time is difficult to model, an identification method was developed to calculate in real time the lag time between the estimator and the measurements without using an additional instrument. The principle is to apply a small increment in the vicinity of the injection point of cylinder  1 , and in calculating the estimated fuel/air ratio variations for each cylinder. Then, an identification criterion J k  is constructed so as to penalize the variations of cylinders  2 ,  3  and  4 . 
                   {           β   =     [     0   ,   1   ,     -   1     ,   2     ]                   J   k     =     β   *     (       Λ   k     -     Λ   0       )                       (   10   )               
with:
 
     
       
         
           
             
               Λ 
               k 
             
             = 
             
               
                 [ 
                 
                   
                     
                       
                         
                           λ 
                           ^ 
                         
                         1 
                       
                     
                   
                   
                     
                       
                         
                           λ 
                           ^ 
                         
                         2 
                       
                     
                   
                   
                     
                       
                         
                           λ 
                           ^ 
                         
                         3 
                       
                     
                   
                   
                     
                       
                         
                           λ 
                           ^ 
                         
                         4 
                       
                     
                   
                 
                 ] 
               
               ⁢ 
               
                 : 
               
               ⁢ 
               
                   
               
               ⁢ 
               composition 
               ⁢ 
               
                   
               
               ⁢ 
               at 
               ⁢ 
               
                   
               
               ⁢ 
               increment 
               ⁢ 
               
                 
                     
                 
                 ⁢ 
                 
                     
                 
               
               ⁢ 
               k 
             
           
         
       
     
     
       
         
           
             
               Λ 
               0 
             
             = 
             
               
                 [ 
                 
                   
                     
                       
                         λ 
                         1 
                         ref 
                       
                     
                   
                   
                     
                       
                         λ 
                         2 
                         ref 
                       
                     
                   
                   
                     
                       
                         λ 
                         3 
                         ref 
                       
                     
                   
                   
                     
                       
                         λ 
                         4 
                         ref 
                       
                     
                   
                 
                 ] 
               
               ⁢ 
               
                 : 
               
               ⁢ 
               
                   
               
               ⁢ 
               reference 
               ⁢ 
               
                   
               
               ⁢ 
               
                 composition 
                 . 
               
             
           
         
       
     
     The penalization is given by β. If there is a positive variation of the fuel/air ratio value estimated for cylinder  2 , the lag time between the estimator and the measurements is positive. If there is a variation on cylinder  3 , the lag time is negative and the penalization is negative. A variation of cylinder  4  can be considered to be a consequence of a positive or negative lag time. Lag time D applied to the output variables of model RTM is an additive delay, it is calculated by least squares by minimizing J k . 
     Criterion J k  is controlled at zero by a controller PI (Integral Proportional) on the estimator lag time. When the controller is stabilized, the estimated fuel/air ratio variation is maximum on cylinder  1  and minimum on cylinder  4 . The estimator is then in phase with the measurements. 
     Results 
       FIGS. 4A to 4F  illustrate the estimation of the fuel/air ratio cylinder by cylinder by means of the estimator described above for two working points, 2000 rpm at high load, 9 bar ( FIGS. 4A to 4C ) and 2500 rpm at low load, 3 bar ( FIGS. 4D to 4F ). These figures show, at the top, the reference fuel/air ratios (λ i   ref ) as a function of time (T) and, at the bottom, the results of the estimator ({circumflex over (λ)} i ) as a function of time (T). The four curves correspond to each one of the four cylinders. 
     The present invention relates to an estimation method comprising construction of an estimator allowing, from the fuel/air ratio measured by the probe (λ) and the information on the total mass of gas inside the manifold (M T ), to estimate the fuel/air ratios at the outlet of the four cylinders (λ i ). The estimator thus achieved is efficient and, above all, it requires no additional adjustment in case of a working point change. No identification stage is necessary, a single measurement noise and model adjustment only has to be performed. 
     In order to make the estimation according to the invention more robust, whatever the operating conditions, a lag time controller is used in parallel with the estimator, allowing to re-adjust the lag time after an injection time increment on a cylinder. This allows optimum calibration of the estimator, for example before a fuel/air ratio 1 phase. 
     The invention also allows performing a measurement every 6° of crankshaft rotation and thus to have high-frequency information of the fuel/air ratio measurement without however being affected by the measurement noise. Furthermore, the high-frequency representation allows to accounting for the pulsating effect of the system. The modelled system is periodic and allows obtaining an estimator with better dynamics: the exhaust pulsation is anticipated. 
     Besides, the invention allows the calculating time to be reduced by approximately a factor of 80 in relation to prior methods.