Method and associated device for sensing the air/fuel ratio of an internal combustion engine

FIELD OF THE INVENTION

0001This invention relates to control systems for the operating parameters of internal combustion engines, and, more particularly, to a method and associated device for sensing the air/fuel ratio (briefly AFR) of an internal combustion engine, and an associated control system that uses this sensing device.

BACKGROUND OF THE INVENTION

0002In the last twenty years, fundamental goals of engine manufacturers are to achieve significant reductions of the amounts of pollutants emitted at the engine exhaust, and lower fuel consumption without compromising speed and torque performances. For these reasons, an efficient engine control based on a comprehensive monitoring of the many engine working parameters is desired.

0003To maintain a strict control of the engine working parameters, Engine Management Systems (EMS) or Engine Control Units (ECU) are used. The EMS implements control strategies which achieve the optimum trade-off between several contradictory objectives: high output power when required by the driver, low emission levels and low fuel consumption. At the same time, in a spark-ignition engine, the EMS brings and maintains the engine in a specified operating range such that the three-way catalytic converter can further reduce the undesired content of the exhaust gases. The EMS controls the amount of fuel injected in the engine combustion chamber (fuel pulse width), the point in the engine cycle at which the mixture air fuel is ignited (ignition timing) and in advanced engine designs, other parameters, such as the valve timing. The EMS determines values for these parameters from measured quantities such as speed, torque, air mass flow rate, inlet-manifold pressure, temperatures at several critical points and throttle angle.

0004<figref idref="DRAWINGS">FIG. 1</figref> illustrates the EMS function. The EMS determines values for Controlled Variables from knowledge of the Measured Variables to achieve the System Aims. EMS essentially includes three components: engine maps (look-up tables stored in a ROM), a controller and sensors, as schematically depicted in <figref idref="DRAWINGS">FIG. 2</figref>.

0005In addition to sensors for measuring quantities of interest, such as speed, manifold pressure, air mass flow rate, temperature (that is, the Measured Variables appearing in both <figref idref="DRAWINGS">FIGS. 1 and 2</figref>), in <figref idref="DRAWINGS">FIG. 2</figref> appear other sensors too. These additional devices monitor whether the engine is working according to the System Aims, or not. Therefore, they have an active part in the real time updating process of controlled variables and, eventually of the engine maps. For example, in a spark-ignition engine a sensor of this type is the so-called lambda sensor. The lambda sensor, mounted in the exhaust stream as schematically shown in the block diagram of <figref idref="DRAWINGS">FIG. 3</figref>, determines whether the lambda ratio (that is AFR/AFR_{stoichiometric}) is above or below unity from the amount of oxygen detected in the exhaust gas mixture. The EMS uses this information to adjust the fuel pulse width and/or the ignition timing to keep the lambda ratio as close as possible to unity.

0006To keep the air/fuel ratio (AFR) within such a narrow range, a lambda sensor is inserted in the outlet of exhaust gases for monitoring the amount of oxygen in the exhaust gases. The lambda sensor provides a signal representative of the value of the ratio

0007<maths id="MATH-US-00001" num="00001"><math overflow="scroll"><mrow><mi>λ</mi><mo>=</mo><mfrac><mrow><mi>Air</mi><mo>⁢</mo><mstyle><mtext>/</mtext></mstyle><mo>⁢</mo><mi>Fuel</mi></mrow><mrow><mi>Air</mi><mo>⁢</mo><mstyle><mtext>/</mtext></mstyle><mo>⁢</mo><msub><mi>Fuel</mi><mi>stoichiometric</mi></msub></mrow></mfrac></mrow></math></maths><br /> If λ<1 the mixture is rich of fuel, while if λ>1 the mixture is lean of fuel, as schematically shown in <figref idref="DRAWINGS">FIG. 4</figref>.

0008The signal generated by the lambda sensor is input to the controller of the engine that adjusts the injection times and thus the fuel injected during each cycle for reaching the condition λ=1.

0009Many lambda sensors actually available, the so-called on/off lambda sensors, do not evaluate the ratio of the mixture and thus the exact value of λ, but signal whether the mixture is reach or lean. Considering that the injection time should ideally be proportional to the air/fuel ratio, these on/off lambda sensors do not allow a precise regulation.

0010There are lambda sensors that generate a signal representative of the effective value of the air/fuel ratio, but these lambda sensors (the so-called “wide-band lambda sensors”) are either very expensive or not very accurate. The following table compares costs and accuracies of commercially available “wide-band lambda sensors”:

0011<tables id="TABLE-US-00001" num="00001"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="5"><colspec colname="offset" colwidth="56pt" align="left" /><colspec colname="1" colwidth="35pt" align="center" /><colspec colname="2" colwidth="49pt" align="center" /><colspec colname="3" colwidth="35pt" align="center" /><colspec colname="4" colwidth="42pt" align="center" /><thead><row><entry /><entry namest="offset" nameend="4" align="center" rowsep="1" /></row><row><entry /><entry>accuracy</entry><entry>accuracy for</entry><entry>accuracy</entry><entry /></row><row><entry /><entry>for lean</entry><entry>stoichiometric</entry><entry>for rich</entry><entry>cost</entry></row><row><entry /><entry>mixtures</entry><entry>mixtures</entry><entry>mixtures</entry><entry>(USD)</entry></row><row><entry /><entry namest="offset" nameend="4" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry /></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="5"><colspec colname="1" colwidth="56pt" align="left" /><colspec colname="2" colwidth="35pt" align="char" char="." /><colspec colname="3" colwidth="49pt" align="char" char="." /><colspec colname="4" colwidth="35pt" align="char" char="." /><colspec colname="5" colwidth="42pt" align="center" /><tbody valign="top"><row><entry>McLaren</entry><entry>1.7%</entry><entry>0.1%</entry><entry>1.7%</entry><entry>1500-1800</entry></row><row><entry>electronic systems</entry></row><row><entry>MoTeC</entry><entry>2.5%</entry><entry>1.75% </entry><entry>1.75% </entry><entry>800-900</entry></row><row><entry>Bosch LSM 11</entry><entry>1.5%</entry><entry>unknown</entry><entry>unknown</entry><entry>300-400</entry></row><row><entry>Horiba LD-700</entry><entry>8.0%</entry><entry>4.0%</entry><entry>8.0%</entry><entry>60-80</entry></row><row><entry namest="1" nameend="5" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

0012Engines manufacturers are generally reluctant to a proliferation of sensors unless they produce valuable improvements that could not otherwise be attained. Virtual-sensors techniques are generally welcome because of their comparably lower cost, reliability and sturdiness. Virtual-sensors allow estimates of quantities of interest without the necessity for sensors dedicated to the measurements. In this field, intelligent systems models, such as neural networks, are attractive because of their capabilities in pattern recognition and signal analysis problems [1].

0013An approach to realize a virtual lambda sensor uses neural networks to correlate certain features of spark plug voltage waveforms with specific values of air fuel ratio [2], [3]. The spark plug is in direct contact with the combustion processes which are occurring in the engine cylinder, hence analysis of the spark plug voltage waveforms seems to be potentially a suitable method of monitoring combustion in spark ignition engines.

0014There are essentially two methods of using a spark plug as a combustion sensor, namely: the Ionic-Current and Spark Voltage Characterization (SVC) methods. In the ionic-current system, the spark plug is used as a sensor during the “non-firing” phase of the pressure cycle, which is the part of the pressure cycle after the spark advance, that is, after the spark ignition. This is done by applying a small voltage of about 100 Volts to the spark plug and measuring the current. The current is supported by reactive ions in the flame that carry on ionic current across the spark plug gap. The type and the number of ions, formed during and after the combustion, depends on the combustion conditions. The Ionic-Current depends also on other parameters such as temperature, pressure and other. Recently, much work has been done on the use of Ionic-Current for monitoring combustion [4], [5], [6] [7].

0015The SVC method rests on the analysis of the time-varying voltage detected across the gap of the spark plug. Since the SVC method involves the analysis of the ignition voltage waveform itself, it does not require additional biasing means and associated high voltage switching circuitry.

0016<figref idref="DRAWINGS">FIG. 5</figref> illustrates a typical spark voltage waveform. The shape of spark voltage waveform has several predictable phases. When the EHT (Extra High Tension) pulse is generated, the potential difference across the gap rises up to 12 kV and breakdown occurs. Breakdown is a fall in voltage that produces a characteristic voltage spike of about 10 μs in duration. Thereafter, a glow-discharge tail region of the waveform of a few milliseconds duration appears. Tests have demonstrated that changes of engine working parameters lead to changes of the shape of certain features of the waveform. However, it is far from being easy to predict these variations as the engine parameters are varied. In fact, random variations occur between successive sparks even when engine working parameters are kept constant.

0017Interactions of parameters, such as combustion temperatures, compression, composition of the air-fuel gas mixture, affect the shape of the breakdown voltage spike in the spark voltage waveform. Changes of the lambda ratio lead to breakdown voltage changes and to subtle changes in the overall shape of the ignition spark waveform. Lambda ratio changes appear to affect both the shapes of the breakdown voltage spike and of the flow-discharge tail portion of the waveform. An analytic relationship between lambda values and instantaneous voltage values of the spark voltage waveforms has not been found yet. However, several articles ([8] and [9]) sustain a correlation between the vector formed through a periodic sampling of the spark plug voltage (spark-voltage vector) and lambda values.

0018The Spark Voltage Characterization (SVC) technique is based on setting up an effective neural network for associating the spark-voltage vector and lambda ratio.

0000AFR Estimation Using Spark Voltage Characterization by Neural Network

0019According to R. J. Howlett et al. in [8], [9], and [10] it is possible to design a Virtual Lambda Sensor, that is a device for sensing the air/fuel without analyzing the exhaust gases of the engine.

0020Such a virtual sensor is based on a neural network trained to find the best correlation between characteristic aspects of the spark voltage waveform and lambda values. The trained neural network determines, for a current vector of characteristic values of the spark voltage, whether the air/fuel ratio (lambda value) is in the stoichiometric mixture range or in lean or rich mixture ranges.

0021<figref idref="DRAWINGS">FIG. 6</figref> shows a typical experimental arrangement to acquire data for training of virtual lambda sensor models. The dynamometer, by which an engine “dummy” load may be varied as desired, is used to measure load-torque and to calculate the output power. Setting of throttle position and fuel pulse width allows changing the air-fuel ratio. In this way, a data set related to the whole range of lambda values may be established.

0022The blocks EMU, A-D converter and DSP are an Engine Management Unit, Analog-to-Digital converter and Digital Signal Processor, respectively.

0023Air-fuel ratio values are measured by an exhaust gas analyzer. To measure spark plug voltage the ignition system is modified by the addition of a high-voltage test-probe at the spark plug.

0024In these approaches, a MLP (Multiple Layer Perceptron) neural network, with a single hidden layer unit and sigmoidal activation unit, is used as a spark-voltage vector classifier.

0025In a supervised training paradigm, a back-propagation learning algorithm sets the MLP training. The training file contains N_tpairs input-output; model input is an instantaneous spark-voltage vector of the form V_i=(ν₁, ν₂, . . . , ν_m), with i=1, . . . , N_tand m equal to the length of the spark-voltage vector; model output is a desired output vector of the form D_r=(0,0,1), D_stoi=(0,1,0) and D_l,=(1,0,0), depending on whether the lambda value, associated to the current spark-voltage vector, is rich (<1), stoichiometric (≈1) or lean (>1).

0026Three sets of spark-voltage vectors and their associated desired output vectors build the training file. Similar files, built by data not to be used for training, are created for validation and test purposes. In this case, during the testing phase, to estimate the model forecast capability it is sufficient to count the number of times in which model output doesn't match the desired output value. The ratio between this number and the total number of estimates represents the model classification error. An alternative quantity for describing the model forecast capability can be simply obtained as difference between 1 and the classification error. This alternative quantity is usually called correct classification rate.

0027R. J. Howlett et al. [8], [9] carried out a multi-speed test with the same 92 cc single-cylinder four-stroke engine. In this case, they used a more closely-spaced range of lambda values, i.e. 0.9, 1.0 and 1.1. <figref idref="DRAWINGS">FIG. 7</figref> shows the trend of the correct classification rate of the virtual lambda sensor model versus engine speed for various model training file sizes. The normalized size of the training file σ, used in this test, satisfies the following relationship N_t=N_wσ, where N_tis the size of training file and N_wis the number of weights of the MLP neural network modeling virtual lambda sensor.

0028These approaches have important drawbacks. The above virtual lambda sensors are unable to indicate the actual AFR but only if the AFR is in one or the other range. In other words, they cannot confirm lambda values approximately equal to 0.95 or 1.05 as illustrated by the rectangles in <figref idref="DRAWINGS">FIG. 8</figref>.

0029The number of cycles of integration, according to the approach aimed at reducing the effect of random variations observed in successive spark waveforms, is not specified. However, this would be an important parameter when realizing a fast gasoline engine injection control system.

0030The forecast capability of the system of R. J. Howlett et al. [8-9] has a strong dependence on engine speed.

0031It has been shown [20] that at an MBT condition (Maximum spark advance, evaluated in respect to the TDC for the Best Torque) the pressure peak in a cylinder during combustion is correlated with the air/fuel ratio, while the location of the pressure peak at a fixed air/fuel ratio value is correlated with the spark advance. Therefore, it is possible to regulate the air/fuel ratio at stoichiometric conditions simply by correcting the fuel injection to keep constant the position of the crank at which the pressure peak is attained, and keeping the pressure peak at a certain value.

0032The so-called MBT condition is the operating condition of the engine when the spark advance takes on the maximum value before bringing the engine toward the knocking phenomena. Normally, this condition is not often verified during the functioning of the engine.

0033In [20], a neural network for sensing the position of the crank when the pressure peak occurs (that is the Location of the Pressure Peak, or briefly the LPP parameter) and the pressure peak value (briefly, the PP parameter) is also disclosed. This neural network is embodied in an air/fuel ratio feedback regulator, and provides to a control system of the engine, signals representing the LPP and the PP parameters. This control system drives the motor in order to keep constant the LPP parameter and to keep constant the air/fuel ratio by regulating the pressure peak in the cylinders.

0034Unfortunately, this document though establishing that there is only a relationship between the air/fuel ratio and the pressure peak if the LPP parameter of the motor is constant (in particular, if the LPP parameter corresponds to the value for MBT condition), is silent about any possibility of assessing the actual air/fuel ratio as a function of the pressure peak without employing a classic lambda sensor under any condition of operation of the engine.

0035As a matter of fact, the correlation between the pressure peak and the air/fuel ratio has been demonstrated only in steady-states at certain operating conditions, that is, at MBT conditions, at 2000 rpm and MAP of 0.5 and 0.8 bar.

0036The system disclosed in that document does not lend itself for sensing the air/fuel ratio, that is for generating a signal that represents at each instant the current value of the air/fuel ratio of the engine.

0037Therefore, the need remains for a low cost manner of sensing the air/fuel ratio with a sufficient accuracy under any condition of operation of the engine.

SUMMARY OF THE INVENTION

0038It has been found a method of sensing the air/fuel ratio in a combustion chamber of an internal combustion engine that may be easily implemented by a respective low-cost device.

0039The device of the invention has a pressure sensor and a learning machine that generates a sensing signal representing the air/fuel ratio by processing the waveform of the pressure in at least one cylinder of the engine. In practice, the learning machine extracts characteristic parameters of the waveform of the pressure and as a function of a certain number of them generates the sensing signal.

0040Surprisingly, the device of this invention is even more accurate than the classic lambda sensors in any operating condition of the engine.

0041The characteristic parameters to be used for sensing the air/fuel ratio are preferably averaged on a certain number of pressure cycles, for reducing noise and improving the accuracy of the sensing. This certain number of pressure cycles is established by using a clustering algorithm on a data set comprising various moving averages of these parameters carried out on different number of pressure cycles.

0042According to an innovative aspect, the learning machine of the device for sensing the air/fuel ratio is based on a kind of neural network, herein referred as MultiSpread-PNN. An appropriate method of training such a neural network is also disclosed.

0043The device of the invention is conveniently inserted in a feedforward-and-feedback control system of an engine for regulating its air/fuel ratio at the stoichiometric value. All the methods of this invention may be implemented by a software computer program.

BRIEF DESCRIPTION OF THE DRAWINGS

0044The different aspects and advantages of the invention will become even clearer through a detailed description of practical embodiments referring to the attached drawings, wherein:

0045<figref idref="DRAWINGS">FIG. 1</figref> is a block diagram that illustrates schematically a control of an internal combustion engine by means of an Engine Management System as in the prior art;

0046<figref idref="DRAWINGS">FIG. 2</figref> is a block diagram of an Engine Management System as in the prior art;

0047<figref idref="DRAWINGS">FIG. 3</figref> illustrates an Engine Management System for a spark ignition engine as in the prior art;

0048<figref idref="DRAWINGS">FIG. 4</figref> is a graph that illustrates the conversion percentages of exhaust gases of the engine as a function of the air/fuel ratio as in the prior art;

0049<figref idref="DRAWINGS">FIG. 5</figref> depicts a simple spark voltage waveform of an internal combustion engine as in the prior art;

0050<figref idref="DRAWINGS">FIG. 6</figref> depicts an arrangement for training a lambda sensor based on a neural network as in the prior art;

0051<figref idref="DRAWINGS">FIG. 7</figref> illustrates how the correct classification rate of the lambda sensor of <figref idref="DRAWINGS">FIG. 6</figref> varies as a function of the engine speed as in the prior art;

0052<figref idref="DRAWINGS">FIG. 8</figref> shows schematically how a Multi-Layer Perceptron neural network of <figref idref="DRAWINGS">FIG. 6</figref> determines the air/fuel ratio of an internal combustion engine as in the prior art;

0053<figref idref="DRAWINGS">FIG. 9</figref> depicts an arrangement of this invention for training a lambda sensor based on a learning machine;

0054<figref idref="DRAWINGS">FIGS. 10</figref><i>a </i>to <b>10</b><i>d </i>are graphs that show results of a clustering algorithm on the values of the parameter Pratio<b>40</b> sensed at each pressure cycle, averaged on five, ten and fifteen pressure cycles, respectively, as in the invention;

0055<figref idref="DRAWINGS">FIG. 11</figref> is a sample graph of the performance index of a chosen clustering algorithm as a function of the number of clusters in which data are grouped as in the invention;

0056<figref idref="DRAWINGS">FIG. 12</figref> is a sample graph of the clustering factor of a chosen clustering algorithm as a function of the number of samples used for calculating the moving averages on which the clustering algorithm is applied, for random initializations of the centers of the clusters, as in the invention;

0057<figref idref="DRAWINGS">FIG. 13</figref> is a sample graph of the clustering factor of a chosen clustering algorithm as a function of the number of samples used for calculating the moving averages on which the clustering algorithm is applied, for deterministic initializations of the centers of the clusters, as in the invention;

0058<figref idref="DRAWINGS">FIG. 14</figref> is a graph showing a sample result of a chosen clustering algorithm on moving averages of the values of the parameters Pratio<b>40</b>, Pratio<b>50</b> and on the value lambda of the air/fuel ratio, as in the invention;

0059<figref idref="DRAWINGS">FIG. 15</figref> is a graph normalized in the range [0; 1] showing a sample result of a chosen clustering algorithm on moving averages of the values of the parameters Pratio<b>40</b>, Pratio<b>50</b> and on the value lambda of the air/fuel ratio, as in the invention;

0060<figref idref="DRAWINGS">FIG. 16</figref> shows schematically how a device of this invention based on a neural network of <figref idref="DRAWINGS">FIG. 9</figref> for sensing the air/fuel ratio of an internal combustion engine works, as in the invention;

0061<figref idref="DRAWINGS">FIG. 17</figref> highlights the incertitude range of the device of the invention;

0062<figref idref="DRAWINGS">FIG. 18</figref> depicts a typical architecture of a classic RBF neural network, as in the prior art;

0063<figref idref="DRAWINGS">FIG. 19</figref> depicts a typical architecture of a classic RBF-PNN neural network as in the prior art;

0064<figref idref="DRAWINGS">FIG. 20</figref> compares a classic RBF-PNN neural network with a MultiSpread-PNN neural network of this invention;

0065<figref idref="DRAWINGS">FIGS. 21 and 22</figref> show test results obtained with a known RBF-PNN neural network and a MultiSpread-PNN neural network of this invention for sensing the air/fuel ratio of an engine;

0066<figref idref="DRAWINGS">FIG. 23</figref> compares the mean test errors of the neural networks of <figref idref="DRAWINGS">FIGS. 21 and 22</figref> as a function of the number of times in which the neural networks RBF-PNN and MultiSpread-PNN have been trained, validated and tested on different training, validation and testing data sets, of this invention;

0067<figref idref="DRAWINGS">FIG. 24</figref> shows a feedforward-and-feedback control system of this invention that uses a device of this invention for sensing the air/fuel ratio of the engine;

0068<figref idref="DRAWINGS">FIG. 25</figref> shows a sample embodiment of the error evaluation subsystem B<b>1</b> of <figref idref="DRAWINGS">FIG. 24</figref>;

0069<figref idref="DRAWINGS">FIG. 26</figref> shows a sample embodiment of the correction subsystem B<b>2</b> of <figref idref="DRAWINGS">FIG. 24</figref>;

0070<figref idref="DRAWINGS">FIGS. 27 to 29</figref> show sample membership functions for the inputs and output of the correction unit CONTROLLER of <figref idref="DRAWINGS">FIG. 26</figref>;

0071<figref idref="DRAWINGS">FIG. 30</figref> shows the nine fuzzy rules of the correction unit CONTROLLER of <figref idref="DRAWINGS">FIG. 26</figref>;

0072<figref idref="DRAWINGS">FIGS. 31 and 32</figref> are graphs that describe the fuzzy rules of the correction unit CONTROLLER of <figref idref="DRAWINGS">FIG. 26</figref>;

0073<figref idref="DRAWINGS">FIG. 33</figref> illustrates schematically how the parameters for training the learning machine of the device of this invention are chosen;

0074<figref idref="DRAWINGS">FIG. 34</figref> illustrates schematically embodiments of the device of this invention for sensing the air/fuel ratio;

0075<figref idref="DRAWINGS">FIG. 35</figref> shows an embodiment of the device of this invention based on a neural network;

0076<figref idref="DRAWINGS">FIGS. 36</figref><i>a </i>and <b>36</b><i>b </i>compare sample graphs of the output generated by the device of <figref idref="DRAWINGS">FIG. 35</figref> without and with the pre-processor and the post-processor, respectively, of this invention;

0077<figref idref="DRAWINGS">FIG. 37</figref> shows an embodiment of the device of this invention based on three fuzzy subsystems;

0078<figref idref="DRAWINGS">FIG. 38</figref> is a graph of sensed values of the air/fuel ratio during a test in which an engine in steady state at 4600 rpm and WOT condition was controlled by the feedforward-and-feedback control system of this invention;

0079<figref idref="DRAWINGS">FIG. 39</figref> is a graph of sensed values of the air/fuel ratio during a test in which an engine in steady state at less than 4600 rpm and throttle not in WOT condition was controlled by the feedforward-and-feedback control system of this invention;

0080<figref idref="DRAWINGS">FIG. 40</figref> is a graph of sensed values of the air/fuel ratio during a test in which an engine in steady state at 4600 rpm and torque of 1.5 Nm was controlled by the feedforward-and-feedback control system of this invention;

0081<figref idref="DRAWINGS">FIG. 41</figref> is a graph of sensed values of the air/fuel ratio during tests in which an engine in transient conditions after variations of the throttle position from 38% to 100% was controlled by the feedforward-and-feedback control system of this invention;

0082<figref idref="DRAWINGS">FIG. 42</figref> is a graph of sensed values of the air/fuel ratio during a test in which an engine in transient conditions at 4600 rpm and WOT condition was controlled by the feedforward-and-feedback control system of this invention;

0083<figref idref="DRAWINGS">FIG. 43</figref> is a graph of sensed values of the air/fuel ratio during a test in which an engine initially functioning with lean mixtures was controlled by the feedforward-and-feedback control system of this invention;

0084<figref idref="DRAWINGS">FIG. 44</figref> displays the software console of the used test system during a test in which an engine in a steady state at 4600 rpm and WOT condition was controlled by the feedforward-and-feedback control system of this invention; and

0085<figref idref="DRAWINGS">FIGS. 45 and 46</figref> display the software console of the used test system during a test in which an engine in transient conditions at 5596 rpm and 4605 rpm, respectively, was controlled by the feedforward-and-feedback control system of this invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

0086According to this invention, the inputs for modeling a virtual lambda sensor are obtained from an engine cylinder pressure signal generated by a pressure sensor, such as for instance the integrated pressure sensor disclosed in [21]. According to the invention, a virtual device capable of sensing the air/fuel ratio is based on a learning machine trained according to the scheme of <figref idref="DRAWINGS">FIG. 9</figref>.

0087Of course, it is possible to train the learning machine using also characteristics of other signals (the speed of the engine, for instance) in addition to the characteristics of the waveform of the pressure in a cylinder, but surprisingly it has been found that the pressure waveform features alone permit to achieve an outstandingly accurate assessment of the air/fuel ratio.

0088Indeed, a wealth of operating parameters of the engine could be extracted from the waveform of pressure in a cylinder. Of course, it is not convenient to consider all of them because the computational load of the learning machine would become excessive.

0089sample set of characteristic parameters that are correlated with the air/fuel ratio are resumed in the following table

0090<tables id="TABLE-US-00002" num="00002"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="3"><colspec colname="offset" colwidth="56pt" align="left" /><colspec colname="1" colwidth="42pt" align="left" /><colspec colname="2" colwidth="119pt" align="left" /><thead><row><entry /><entry namest="offset" nameend="2" align="center" rowsep="1" /></row><row><entry /><entry>unit of</entry><entry /></row><row><entry /><entry>measure</entry><entry>element description</entry></row><row><entry /><entry namest="offset" nameend="2" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry /></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="3"><colspec colname="1" colwidth="56pt" align="left" /><colspec colname="2" colwidth="42pt" align="left" /><colspec colname="3" colwidth="119pt" align="left" /><tbody valign="top"><row><entry>Speed</entry><entry>rpm</entry><entry>engine speed</entry></row><row><entry>lambda</entry><entry>[ ]</entry><entry>lambda values</entry></row><row><entry>Aircycle</entry><entry>mg/cycle</entry><entry>air massive flow</entry></row><row><entry>BstMap</entry><entry>bar</entry><entry>intake manifold pressure</entry></row><row><entry>BurDur</entry><entry>deg</entry><entry>combustion duration</entry></row><row><entry>pEVC</entry><entry>bar</entry><entry>pressure at exhaust valve closure</entry></row><row><entry>pEVO</entry><entry>bar</entry><entry>pressure at exhaust valve opening</entry></row><row><entry>pIVC</entry><entry>bar</entry><entry>pressure at intake valve closure</entry></row><row><entry>pIVO</entry><entry>bar</entry><entry>pressure at intake valve opening</entry></row><row><entry>Pratio40</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 40 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pratio50</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 50 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pratio60</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 60 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pratio70</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 70 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pratio80</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 80 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pratio90</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 90 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pratio100</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 100 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pratio110</entry><entry>[ ]</entry><entry>pressure ratio between pressures</entry></row><row><entry /><entry /><entry>at 110 crank angles before and</entry></row><row><entry /><entry /><entry>after TDC</entry></row><row><entry>Pmax</entry><entry>bar</entry><entry>maximum pressure</entry></row><row><entry>PcompMax</entry><entry>[ ]</entry><entry>ratio between maximum of pressure</entry></row><row><entry /><entry /><entry>cycle and maximum of pressure</entry></row><row><entry /><entry /><entry>cycle without combustion</entry></row><row><entry namest="1" nameend="3" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

0091These parameters have been identified as relevant for estimating the air/fuel ratio during extensive tests carried out on the commercial scooter engine: Yamaha YP125 (four stroke spark ignition engine with a displacement of 125 cc). The tests have been performed at several engine speeds, throttle positions, and spark advances, for considering all possible functioning conditions of the engine.

0092Given that a learning machine processing detected relevant parameters for estimating the air/flow ratio would be relatively slow, a small number of parameters to be used has been chosen.

0093A data pre-processing campaign was carried out to identify the parameters more correlated to the air/flow ratio (lambda value). During each pressure cycle all the above parameters were measured and the corresponding lambda value was sensed by a lambda sensor. The correlation of each parameter with the sensed air/fuel ratio was calculated, and only the three parameters that resulted most correlated with the air/fuel ratio as directly measured by the sensor were chosen as inputs of the learning machine.

0094Of course, it is possible to choose more than three parameters or even two or only one parameter as inputs of the learning machine, but the choice of three parameters appeared a good compromise. While choosing to use a larger number of parameters will increase the computational load, a too small number of parameters may impair the accuracy under varying functioning conditions of the engine.

0095The three parameters most correlated with the air/fuel ratio resulted to be Pratio<b>40</b>, Pratio<b>50</b> and Pmax and, according to a preferred embodiment of this invention, these three parameters were used as the inputs of the learning machine.

0096In view of the fact that the values of these parameters as detectable may be corrupted by noise, it is advisable to use a moving average of the parameter calculated over a certain number of pressure cycles for estimating the air/fuel ratio.

0097A problem faced by the Applicants consisted in determining the number of detections values to be taken for calculating the moving average value to be used for the air/fuel ratio calculation. The larger is the number of successive samples, the more filtered from noise are the inputs of the learning machine, but the less prompt is the tracking of a time-changing air/fuel ratio.

0098In order to find the most effective approach, numerous moving averages with different numbers of samples of the three parameters have been calculated, and the moving average that resulted most correlated with the air/fuel ratio was chosen through a clustering analysis, that will be described in detail below.

0099<figref idref="DRAWINGS">FIGS. 10</figref><i>a </i>to <b>10</b><i>d </i>show the results of a clustering process carried out among moving averages of the parameter Pratio<b>40</b> calculated over different numbers of samples, and the corresponding λ value. It was experimentally determined that the larger the number of samples over which the moving average is carried out, the larger is the correlation between Pratio<b>40</b> and the corresponding λ value.

0100In the graphs, each correlation is represented with a circle of pre-established radius and the circles of each cluster have the same color. From a mathematical point of view, this is equivalent to assume that small input variations generate small output variations, that is equivalent to assume that the air/fuel ratio depends on this parameter through a well-posed mathematical problem ([11], [12) and [13]).

0101A more detailed description of the Yamaha engine data set clustering analysis is presented in detail below. A novel factor, called “clustering factor”, to compare how different data sets fit an user requested clusters number has been used and it has been found that the moving averages should be carried out on 16 successive samples for obtaining the best trade-off between noise filtering and tracking speed.

0000Clustering Factor of M-Dimensional Data

0102Clustering is an important processing tool used in different scientific fields, e.g. the compression of audio and video signals, pattern recognition, computer vision, recognition of medical images, etc. Clustering may be used for data pre-processing too. For example, it can be exploited for optimizing the learning of neural network models or/and to verify data pre-classifications.

0103There are two kinds of data distribution which can be clustered: data sequences and data ensembles. Data sequence means data come from a temporal sampling, e.g. a signal temporal sampling, or also the temporal sequence of engine pressure cycles. On the other hand, data ensemble means data (in an M-dimensional space) that are not temporally linked.

0104There are several clustering algorithms [22], for each of them there is a different approach to find the “optimal data clustering”. But what does the word “clustering” means?

0105Let a data set of N points in the space R^m, X=x₁, . . . , x_N, clustering of X involves a “natural” partition of the ensemble in 1&lt;c&lt;N sub-structures. There are three different ways to obtain these c sub-structures in a given data set. Each way defines suitably the membership matrix (c×N matrix) of X elements versus c clusters. The membership matrix is built by the c×N matrix U_ikbeing k=1, . . . , N and i=1, . . . , c.

0106U_ikelements are suitably linked to the distances of X points from c temporary cluster centers.

0107<maths id="MATH-US-00002" num="00002"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mi>M</mi><mi>P</mi></msub><mo>=</mo><mrow><mo>{</mo><mrow><mrow><mn>0</mn><mo>≤</mo><msub><mi>U</mi><mi>ik</mi></msub><mo>≤</mo><mrow><mn>1</mn><mo>⁢</mo><mrow><mo>∀</mo><mi>i</mi></mrow></mrow></mrow><mo>,</mo><mrow><mi>k</mi><mo>⁢</mo><mrow><mo>∀</mo><mrow><mi>k</mi><mo>⁢</mo><mrow><mo>∃</mo><mrow><mi>i</mi><mo>⁢</mo><mrow><mo></mo><mrow><msub><mi>U</mi><mi>ik</mi></msub><mo>></mo><mn>0</mn></mrow><mo>}</mo></mrow></mrow></mrow></mrow></mrow></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>M</mi><mi>F</mi></msub><mo>=</mo><mrow><mo>{</mo><mrow><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>c</mi></munderover><mo>⁢</mo><msub><mi>U</mi><mi>ik</mi></msub></mrow><mo>=</mo><mrow><mn>1</mn><mo>⁢</mo><mrow><mo>∀</mo><mrow><mrow><mi>k</mi><mo>⁢</mo><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><msub><mi>U</mi><mi>ik</mi></msub></mrow></mrow><mo>></mo><mrow><mn>0</mn><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>⁢</mo><mrow><mo>∀</mo><mi>i</mi></mrow></mrow></mrow></mrow></mrow></mrow><mo>}</mo></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>M</mi><mi>C</mi></msub><mo>=</mo><mrow><mo>{</mo><mrow><mrow><msub><mi>U</mi><mi>ik</mi></msub><mo>=</mo><mrow><mn>0</mn><mo>⁢</mo><mstyle><mspace width="0.8em" height="0.8ex" /></mstyle><mo>⁢</mo><mi>or</mi><mo>⁢</mo><mstyle><mspace width="0.8em" height="0.8ex" /></mstyle><mo>⁢</mo><mn>1</mn><mo>⁢</mo><mrow><mo>∀</mo><mi>i</mi></mrow></mrow></mrow><mo>,</mo><mi>k</mi></mrow><mo>}</mo></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0108The equations (1), (2) e (3) describe the three feasible membership matrices. In M_P=M_{possibilistic}, U_ikelement is the possibility (typicality) that x_kpoint belongs to the i-th sub-structure. In M_F=M_fuzzy, U_ikelement corresponds to the membership probability of x_kpoint to the i-th sub-structure; these probabilities satisfy for each k a normalization condition. Finally, M_C=M_crispis a Boolean matrix in which each element U_ik=1 if and only if x_kbelongs to the current sub-structure. The three matrices are related by the following relation: <br /><i>M</i>_crisp<i>⊂M</i>_fuzzy<i>⊂M</i>_{possibilistic}  (4)

0109Finding the “optimal partition” of a data set X means to find the matrix U_ikwhich better represents the unknown sub-structures of X in comparison to the clustering model that the algorithm induces.

0110Given a data set clustering, we need to introduce some criterions to appraise it.

0111Moreover, other criterions are needed in order to find what strategies to be followed for improving it. Concerning the first of these requests, there are not objective criterions to appraise a data set clustering; in fact, the current criterions depend on the application.

0112To improve a data set clustering, the more used approach is based on an iterative solution searching. An exhaustive solution searching, in the space of the possible solutions, could be too much onerous from a computational viewpoint. Indeed, the total number of partitions in c classes of a data set with N elements is c^N/c!. A sub-optimal approach allows for each iteration to improve the solution going to optimize the selected criterion.

0113Even if the approach does not guarantee the attainment of the absolute optimal solution, it often is used for its low computational complexity. However, a relevant problem of similar approaches is the sensibility to the initial choosing of the clusters.

0114Two clustering algorithm (FCM and FPCM) implementations will be described. The performances of the two algorithms have been compared on a data set obtained from the pressure cycles of the Yamaha YP125 gasoline engine. After having identified the algorithm features, the algorithms are implemented on a space of M-dimensional data. Besides, the influence of the cluster center vector initialization has been analyzed keeping in mind the theoretical results in [22]. At last, a clustering degree measure of a data set X, called “clustering factor” is proposed. This factor seems able to compare the induced clusterings on several data sets X_(i)for a fixed choosing of the number of sub-structures to be found.

0000Fuzzy C-Means Algorithm (FCM)

0115The FCM algorithm is based on Fuzzy System Theory which is used as a precious mathematical tool in several application fields. A fuzzy set is an element set with a “blurred” membership concept. FCM is an iterative procedure based on the idea that clusters can be handled as fuzzy sets. Each point x_k(with k=1, . . . , N) may belong at the same time to different clusters with membership degrees U_kj(with j=1, . . . , c) which change during the procedure. There is only a constraint: for each element x_kand for each algorithm step, the sum of membership degrees must be equal to 1.

0116From a mathematical perspective, the FCM clustering model can be described as an optimization problem with constraints. The objective function to be minimized is:

0117<maths id="MATH-US-00003" num="00003"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><msub><mi>J</mi><mi>m</mi></msub><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>U</mi><mo>,</mo><mi>V</mi></mrow><mo>)</mo></mrow></mrow><mo>=</mo><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>c</mi></munderover><mo>⁢</mo><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><mrow><msubsup><mi>U</mi><mi>ik</mi><mi>m</mi></msubsup><mo>⁢</mo><msubsup><mi>D</mi><mi>ik</mi><mn>2</mn></msubsup></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0118The following constraints are associated to the function (5):

0119<maths id="MATH-US-00004" num="00004"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>c</mi></munderover><mo>⁢</mo><msub><mi>U</mi><mi>ik</mi></msub></mrow><mo>=</mo><mrow><mrow><mn>1</mn><mo>⁢</mo><mstyle><mspace width="1.7em" height="1.7ex" /></mstyle><mo>⁢</mo><mrow><mo>∀</mo><mi>k</mi></mrow></mrow><mo>=</mo><mn>1</mn></mrow></mrow><mo>,</mo><mi>…</mi><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>,</mo><mi>N</mi></mrow></mtd><mtd><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0120In eq. (5) m is the system fuzzyness degree while D_ikmatrix represents the distances between distribution points (x_k) and cluster centers (v_i). For m=0 the fuzzy clusters become classical clusters, that is, each sample belongs only to a cluster. For m>>0 the system fuzzyness level grows. If m→∞ we can observe that the membership degrees of data set points approach to l/c and cluster centers approach to the distribution center. The FCM algorithm optimizes a criterion which is the “fuzzy” version of the “trace criterion” [23]. Since the algorithm depends by the initial cluster centers, we implemented it considering this. In fact, user can choice the initialization way (stochastic and deterministic) of the cluster centers.

0121The FCM procedure can be executed for several clustering strategies of data set X, that is, for different number of clusters (from c_minto c_max). The final result will be a sequence of J_min(c) (with c=c_min, . . . , c_max), each of theme is the function (5) minimum.

0122There is a performance index P(c), given by eq. (7), by which it is possible to find the “optimal” number of clusters.

0123<maths id="MATH-US-00005" num="00005"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>P</mi><mo>⁡</mo><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mrow><msub><mover><mi>J</mi><mo>~</mo></mover><mi>min</mi></msub><mo>⁡</mo><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow></mrow><mo>-</mo><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>c</mi></munderover><mo>⁢</mo><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><mrow><msubsup><mi>U</mi><mi>ik</mi><mi>m</mi></msubsup><mo>·</mo><mrow><mo></mo><mrow><mover><mi>x</mi><mi>_</mi></mover><mo>-</mo><msub><mi>v</mi><mi>i</mi></msub></mrow><mo></mo></mrow></mrow></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0124The “optimal” number of clusters c_opt: is one minimizes the performance index P(c).

0125P(c) has a minimum when data set clustering has a minimum intra-cluster variance (i.e. small values of D_ikin e {tilde over (J)}_min(c)) and a maximum inter-cluster variance (i.e. maximum cluster center distances v_ifrom data set centers <o ostyle="single">x</o>). Hence, a graph of the performance index versus the number of clusters may be plot. <figref idref="DRAWINGS">FIG. 11</figref> shows a given data set grouped with different numbers of clusters (c=2, . . . , 8).

0126Even looking at <figref idref="DRAWINGS">FIG. 11</figref> it is possible to see that the best data set clustering is obtained with four clusters.

0000Fuzzy Possibilistic C-Means Algorithm (FPCM)

0127The FCM is of course a useful data pre-processing tool, but it is burdened by the drawback of being sensible to noisy data and to the outlier problem (see [22] and [24]). There is a lot of clustering algorithms which try to solve these problems, as in [25] and [26].

0128In this FPCM implementation, the Bezdek approach disclosed in [22] has been used. The main feature of this approach consists in using another grouping strategy called typicality. Basically, while the membership degree is a “local” grouping feature, that is, x_khas a probability of belonging to c clusters normalized to 1, the typicality is a grouping feature involved by the same clusters. In other words, it is supposed that the clustering ways of a data set X are established by an observer.

0129For the FCM, the observer is, for every time, in x_kpoint. In this case, the observer thinks to set his membership to c temporary sub-structures with a probability inversely proportional to his distances from cluster centers. There is only a constraint: the membership degrees are normalized to 1.

0130For the FPCM, the observer is not, for every time, only in x_kpoint but also in the i-th cluster center; in this last case, he thinks to set the membership of all X points to the current cluster with a probability inversely proportional to the distances of X points from the observer. There is only a constraint: the typicality degrees are normalized to 1 according to 8.

0131<maths id="MATH-US-00006" num="00006"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><msub><mi>T</mi><mi>ik</mi></msub></mrow><mo>=</mo><mrow><mrow><mn>1</mn><mo>⁢</mo><mstyle><mspace width="1.7em" height="1.7ex" /></mstyle><mo>∀</mo></mrow><mo>=</mo><mn>1</mn></mrow></mrow><mo>,</mo><mi>…</mi><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>,</mo><mi>c</mi></mrow></mtd><mtd><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0132From a mathematical viewpoint, the FPCM clustering model can be described as an optimization problem with constraints. The objective function to be minimized is:

0133<maths id="MATH-US-00007" num="00007"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><msub><mi>J</mi><mrow><mi>m</mi><mo>,</mo><mi>η</mi></mrow></msub><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>U</mi><mo>,</mo><mi>T</mi><mo>,</mo><mi>V</mi></mrow><mo>)</mo></mrow></mrow><mo>=</mo><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>c</mi></munderover><mo>⁢</mo><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><mrow><mrow><mo>(</mo><mrow><msubsup><mi>u</mi><mi>ik</mi><mi>m</mi></msubsup><mo>+</mo><msubsup><mi>t</mi><mi>ik</mi><mi>η</mi></msubsup></mrow><mo>)</mo></mrow><mo>⁢</mo><msubsup><mi>D</mi><mi>ik</mi><mn>2</mn></msubsup></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0134The following constraints are associated to the function (9):

0135<maths id="MATH-US-00008" num="00008"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>c</mi></munderover><mo>⁢</mo><msub><mi>U</mi><mi>ik</mi></msub></mrow><mo>=</mo><mrow><mn>1</mn><mo>⁢</mo><mrow><mo>∀</mo><mi>k</mi></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><msub><mi>T</mi><mi>ik</mi></msub></mrow><mo>=</mo><mrow><mn>1</mn><mo>⁢</mo><mrow><mo>∀</mo><mi>i</mi></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0136In eq. (9) m and η are the system fuzzyness degree. Since the algorithm depends on the initial centers of clusters, we implemented it considering this. In fact, an user can choose the initialization way (stochastic or deterministic) of the cluster centers.

0000Data Set Clustering Factor

0137The innovation consists in the introduction of a measure of the clustering degree of a given data set X, called “clustering factor”. This factor is useful to compare the same clustering on several data sets. The idea is to divide the performance index P(c) given by eq. (7) by its asymptotic behavior P_asym.(c). P_asym.(c) is P(c) estimated when data set clustering is “ideal”. “Ideal” clustering means that the grouping of data set points is in sub-structures with the minimum intra-cluster variances (i.e. data set points falling on the cluster centers) and the maximum inter-cluster variances (i.e. maximum difference of the features, that is, maximum distances of the cluster centers from data set centers).

0138For an “ideal” clustering, supposing that n_ipoints, amongst data set elements, fall on i-th cluster (with i=1, . . . , c), the membership and typicality matrices (U_ikand T_ik) take the following form:

0139<maths id="MATH-US-00009" num="00009"><math overflow="scroll"><mtable><mtr><mtd><mrow><mo>(</mo><mtable><mtr><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mi>⋯</mi></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mi>⋯</mi></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mi>⋯</mi></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋯</mi></mtd><mtd><mi>⋮</mi></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mi>⋯</mi></mtd><mtd><mn>0</mn></mtd></mtr></mtable><mo>)</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>(</mo><mtable><mtr><mtd><mn>0</mn></mtd><mtd><mfrac><mn>1</mn><msub><mi>n</mi><mn>1</mn></msub></mfrac></mtd><mtd><mn>0</mn></mtd><mtd><mfrac><mn>1</mn><msub><mi>n</mi><mn>1</mn></msub></mfrac></mtd><mtd><mn>0</mn></mtd><mtd><mi>⋯</mi></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mfrac><mn>1</mn><msub><mi>n</mi><mn>2</mn></msub></mfrac></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mi>⋯</mi></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mfrac><mn>1</mn><msub><mi>n</mi><mn>3</mn></msub></mfrac></mtd><mtd><mi>⋯</mi></mtd><mtd><mfrac><mn>1</mn><msub><mi>n</mi><mn>3</mn></msub></mfrac></mtd></mtr><mtr><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mfrac><mn>1</mn><msub><mi>n</mi><mi>c</mi></msub></mfrac></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mi>⋯</mi></mtd><mtd><mn>0</mn></mtd></mtr></mtable><mo>)</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0140Data set points n_ifalling on c clusters must satisfy the constraint:

0141<maths id="MATH-US-00010" num="00010"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>c</mi></munderover><mo>⁢</mo><msub><mi>n</mi><mi>i</mi></msub></mrow><mo>=</mo><mi>N</mi></mrow></mtd><mtd><mrow><mo>(</mo><mn>14</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0142The N elements of a data set can fall in c clusters in c^N/c! different ways. Each falling way must satisfy the constraint in eq. (14). This means that the number of the “ideal” partitions of data set in c sub-structures is equal to c^N/c!. To build P_asint.(c) amongst the c^N/c! “ideal” partitions, we choose that whereby

0143<maths id="MATH-US-00011" num="00011"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mi>n</mi><mn>1</mn></msub><mo>=</mo><mrow><msub><mi>n</mi><mn>2</mn></msub><mo>=</mo><mrow><mi>…</mi><mo>=</mo><mrow><msub><mi>n</mi><mi>c</mi></msub><mo>=</mo><mfrac><mi>N</mi><mi>c</mi></mfrac></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>15</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0144Considering eqs. (12), (13), (14) and (15) it is simple to obtain the asymptotic performance index of the FCM and FPCM algorithms (P_{a sin t.}^FCM(c) and P_{a sin t.}^FPCM(c)):

0145<maths id="MATH-US-00012" num="00012"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><msubsup><mi>P</mi><mrow><mi>asint</mi><mo>.</mo></mrow><mi>FCM</mi></msubsup><mo>⁡</mo><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mo>-</mo><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><msup><mrow><mo></mo><mrow><mover><mi>x</mi><mi>_</mi></mover><mo>-</mo><msub><mi>v</mi><mi>i</mi></msub></mrow><mo></mo></mrow><mn>2</mn></msup></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mrow><msubsup><mi>P</mi><mi>asint</mi><mi>FPCM</mi></msubsup><mo>⁡</mo><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mrow><mo>-</mo><mrow><mo>(</mo><mrow><mfrac><mi>N</mi><mi>c</mi></mfrac><mo>+</mo><msup><mrow><mo>(</mo><mfrac><mi>c</mi><mi>N</mi></mfrac><mo>)</mo></mrow><mrow><mi>η</mi><mo>-</mo><mn>1</mn></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo>·</mo><mrow><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><msup><mrow><mo></mo><mrow><mover><mi>x</mi><mi>_</mi></mover><mo>-</mo><msub><mi>v</mi><mi>i</mi></msub></mrow><mo></mo></mrow><mn>2</mn></msup></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0146To obtain the “clustering factor” of a data set it is necessary to divide the performance index by its asymptotic form. “Clustering factor” is always in [0, 1]. It is able to recognize amongst several data sets X_i, which have been clustered by the same clustering algorithm and according to the same user requested number of clusters, the one which better fits the “ideal” clustering in c sub-structures.

0000Data Set Pre-Processing

0147The considered data set was a data ensemble extracted from the pressure cycles of the test engine. Pre-processing of the data set, by FCM and FPCM clustering algorithms, found the inputs most correlated to the model output (lambda values) but also found the number of the pressure cycle instantaneous values that were to be averaged to obtain the best correlation between VLS inputs and output.

0148The clustering process between Pratio<b>40</b> (a possible input of the VLS model) and λ versus the number of the pressure cycle instantaneous values which we averaged on was analyzed by increasing the number of cycles, the correlation between Pratio<b>40</b> and λ increases strongly.

0149The correlation between Pratio<b>40</b> and λ does not increase indefinitely but has a maximum. The maximum was found when the number of successive pressure cycles (averaged samples) was 16. This was established with the “clustering factor”.

0150Having fixed the model input (as Pratio<b>40</b>) and the cluster number to be induced in the data set (in this sample case, a partition of data set in 3 clusters is to be induced), the different data sets have been labeled with the respective number of pressure cycles on which each input parameter has been averaged.

0151Several cluster center initializations have been used. <figref idref="DRAWINGS">FIGS. 12 and 13</figref> are graphic representations obtained by simulations of values of clustering factor as a function of the number of samples on which parameters are averaged, for random and deterministic initialization of cluster centers, respectively. The maximum value of the “clustering factor” is attained for a number of the pressure cycles, on which each input parameter is to be averaged, equal to 16. Therefore, the characteristic parameters are best averaged on 16 pressure cycles, because this number better fits the clustering requested by the user.

0152Preferably, evolutionary algorithms are used to search the optimal design of a neural network model, which emulates an on-off lambda sensor.

0153In general, the clustering process of data sets depends on the scaling of data values.

0154Therefore, the algorithm input data should be previously normalized in [0, 1]. <figref idref="DRAWINGS">FIGS. 14 and 15</figref> show the difference of the results that were obtained by the same algorithm executed on two different data sets. In the first data set, data (λ, Pratio<b>40</b> and Pratio<b>50</b> values) have values within different ranges, in fact Pratio<b>40</b> and Pratio<b>50</b> values are in [2, 4.5] while λ values are in [0.8, 1.2]. In the second data set, the same data are normalized in [0, 1]. The six clusters, selected by the clustering algorithm for the first data set, are arranged along Pratio<b>40</b> axis. This indicates a data set grouping polarization induced by data component which has higher values. Instead, the four clusters selected by the algorithm in the second data set are arranged along the λ axis and they correspond to the three lambda values ranges useful to the setting-up of a VLS model: rich mixture range (orange and green cluster), stoichiometric range (brown cluster) and lean mixture range (blue cluster). Indeed, a model, which emulates effectively a real lambda sensor, must distinguish the three regions above mentioned in the input space. In this way, the model is almost able to forecast if the engine is above or below the stoichiometric region (i.e. emulation of an on-off lambda sensor).

0155There is a definitive difference between the computational complexity of the FCM and FPCM algorithms. Accordingly, the FCM could be chosen for decreasing the computing time. Moreover, the cluster centers initialization could slow down the algorithm convergence velocity.

0156In practice, there are several critical factors that limit the use of the clustering algorithms for data pre-processing in real time systems. On the other side, clustering algorithms remain an important tool for pre-processing data for the off-line learning of models, e.g. models having applications in the automotive field.

0000A Learning Machine

0157According to an embodiment of this invention, the learning machine is based on a new kind of working logic, substantially different from that of models proposed by R. J. Howlett et al.: such a learning machine is herein referred as Multi-Spread Probabilistic Neural Network.

0158As depicted in <figref idref="DRAWINGS">FIGS. 16 and 17</figref>, the gist of the model of virtual lambda sensor of this invention is to group lambda values in two classes: the class λ=1 in which lambda values are above the stoichiometric range and the class λ=0 in which lambda values are below the stoichiometric, according to inputs values selected by an initial data pre-processing phase.

0159Differently from the known models, incertitude region of the model of this invention is an extremely narrow range around the stoichiometric lambda (λ=1.0). According to the embodiment analyzed, this range corresponds to lambda values within 0.98 and 1.02.

0160In the incertitude region of the model of virtual lambda sensor of this invention, the model forecast capability is significantly lower than the forecast capability that the prior art models have in their working regions. By contrast, in the working regions (red rectangles in <figref idref="DRAWINGS">FIG. 17</figref>) the model of this invention has an outstanding forecast capability that is not tied bounded to a particular speed of engine nor to any particular position of the throttle.

0161From a mathematical viewpoint, a neural network with a scalar output can be described as an hyper-surface Γ in R^m+1space, that is a map s: R^m→R¹. In this formalism, the index m represents design space dimension. Neural network design can be described as a multivariable interpolation in high dimensional space. Given a set of N different points {x⁽ⁱ⁾εR^m|i=1,2, . . . ,N} and a corresponding set of N real numbers {d_iεR¹|i=1,2, . . . ,N}, it is necessary to find a function F: R^m→R¹that satisfies the interpolation conditions: <br /><i>F</i>(<i>x</i>⁽ⁱ⁾)=<i>d</i>_i<i>, i=</i>1, 2<i>, . . . , N</i>  (17)

0162For RBF (Radial Basis Function) neural networks, the map F(x) has the following form, [14] and [15]:

0163<maths id="MATH-US-00013" num="00013"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>F</mi><mo>⁡</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><mrow><msub><mi>w</mi><mi>i</mi></msub><mo>⁢</mo><mrow><mi>ϕ</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mo></mo><mrow><mi>x</mi><mo>-</mo><msup><mi>x</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></msup></mrow><mo></mo></mrow><mo>)</mo></mrow></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>18</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths><br /> where {φ(∥x−x⁽ⁱ⁾∥)|i=1,2, . . . ,N} is a set of N arbitrary nonlinear functions, known as radial basis function, the symbol ∥.∥ denotes a norm, that is usually Euclidean, and the points x⁽ⁱ⁾are the centers of these functions. <figref idref="DRAWINGS">FIG. 18</figref> depicts the architecture of the neural network RBF described by M. J. D. Powell equation (18).

0164Optimal values of the coefficients w_iare determined by interpolation conditions, that is:

0165<maths id="MATH-US-00014" num="00014"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mrow><mo>[</mo><mtable><mtr><mtd><msub><mi>ϕ</mi><mn>11</mn></msub></mtd><mtd><msub><mi>ϕ</mi><mn>12</mn></msub></mtd><mtd><mi>…</mi></mtd><mtd><msub><mi>ϕ</mi><mrow><mn>1</mn><mo>⁢</mo><mi>N</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mi>ϕ</mi><mn>21</mn></msub></mtd><mtd><msub><mi>ϕ</mi><mn>22</mn></msub></mtd><mtd><mi>…</mi></mtd><mtd><msub><mi>ϕ</mi><mrow><mn>2</mn><mo>⁢</mo><mi>N</mi></mrow></msub></mtd></mtr><mtr><mtd><mi>⋮</mi></mtd><mtd><mi>⋮</mi></mtd><mtd><mi>⋰</mi></mtd><mtd><mi>⋮</mi></mtd></mtr><mtr><mtd><msub><mi>ϕ</mi><mrow><mi>N</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mn>1</mn></mrow></msub></mtd><mtd><msub><mi>ϕ</mi><mrow><mi>N</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mn>2</mn></mrow></msub></mtd><mtd><mi>…</mi></mtd><mtd><msub><mi>ϕ</mi><mi>NN</mi></msub></mtd></mtr></mtable><mo>]</mo></mrow><mo>·</mo><mrow><mo>[</mo><mtable><mtr><mtd><msub><mi>w</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><msub><mi>w</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mi>⋮</mi></mtd></mtr><mtr><mtd><msub><mi>w</mi><mi>N</mi></msub></mtd></mtr></mtable><mo>]</mo></mrow></mrow><mo>=</mo><mrow><mo>[</mo><mtable><mtr><mtd><msub><mi>d</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><msub><mi>d</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mi>⋮</mi></mtd></mtr><mtr><mtd><msub><mi>d</mi><mi>N</mi></msub></mtd></mtr></mtable><mo>]</mo></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths><br /> where φ_ij=φ(∥x^(j)−x⁽ⁱ⁾∥). We may rewrite in matrix form the previous equation: <br />Φ·<i>w=d</i>  (20)

0166Assuming that interpolation matrix Φ is nonsingular, we have that: <br /><i>w=d·Φ</i>⁻¹  (21)

0167In literature, there are several classes of functions for which interpolation matrix Φ is always invertible: <ul id="ul0001" list-style="none"><li id="ul0001-0001" num="0000"><ul id="ul0002" list-style="none"><li id="ul0002-0001" num="0168">Multiquadrics: φ(r)=√{square root over (r²+c²)} for some c>0 and rεR</li><li id="ul0002-0002" num="0169">Inverse multiquadrics:</li></ul></li></ul>

0170<maths id="MATH-US-00015" num="00015"><math overflow="scroll"><mrow><mrow><mi>ϕ</mi><mo>⁡</mo><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow><mo>=</mo><mfrac><mn>1</mn><msqrt><mrow><msup><mi>r</mi><mn>2</mn></msup><mo>+</mo><msup><mi>c</mi><mn>2</mn></msup></mrow></msqrt></mfrac></mrow></math></maths><ul id="ul0003" list-style="none"><li id="ul0003-0001" num="0000"><ul id="ul0004" list-style="none"><li id="ul0004-0001" num="0171"> for some c>0 and rεR</li><li id="ul0004-0002" num="0172">Gaussian functions:</li></ul></li></ul>

0173<maths id="MATH-US-00016" num="00016"><math overflow="scroll"><mrow><mrow><mi>ϕ</mi><mo>⁡</mo><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mi>exp</mi><mo>(</mo><mrow><mo>-</mo><mfrac><msup><mi>r</mi><mn>2</mn></msup><mrow><mn>2</mn><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><msup><mi>σ</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></math></maths><ul id="ul0005" list-style="none"><li id="ul0005-0001" num="0000"><ul id="ul0006" list-style="none"><li id="ul0006-0001" num="0174"> for some σ>0 and rεR</li></ul></li></ul>

0175For a more detailed description of the functions classes, as previously mentioned, the interested reader is addressed to the work [16].

0176There is a special class of RBF neural networks, known as RBF-PNN, where the acronym PNN means Probabilistic Neural Network. These networks are used to solve classification problems, that is they work as classifiers.

0177From mathematics viewpoint, a classification problem can be formalized in the following way. Given a set of points {X⊂R^m|x⁽ⁱ⁾εR^m∀i=1,2, . . . ,N}, a clustering process induces a partition of x in I&lt;c&lt;N sub-structures. Membership of x points to c sub-structures, determined by clustering process, is fixed by membership matrix U_ik, where k=1, . . . , c and i=1, . . . , N.

0178Matrix element U_ikrepresents the probability that the i-th point belongs to c-th cluster. Usually, matrix elements U_iksatisfy some normalization conditions.

0179These conditions and the different ways by which a clustering process can be performed distinguish several clustering algorithms known in literature.

0180<figref idref="DRAWINGS">FIG. 19</figref> shows the typical architecture of a RBF neural network. In this scheme, after the hidden layer there is a block which computes membership matrix elements corresponding to input vector. There are several paradigms used to build U_ikmatrix. For RBF-PNN networks, the most common is “the winner takes all neurons”. In other words, we simply evaluate hidden layer neurons outputs, hence we assign to the input vector the membership class of hidden layer neurons which has the largest output. Hidden layer neurons outputs depend on the distance of the input vector from the centers of radial basis functions and on spreads of the same functions. A radial basis function is a Gaussian function having the following form:

0181<maths id="MATH-US-00017" num="00017"><math overflow="scroll"><mtable><mtr><mtd><mtable><mtr><mtd><mrow><msub><mi>ϕ</mi><mi>k</mi></msub><mo>=</mo><mrow><mi>exp</mi><mo>(</mo><mrow><mo>-</mo><mfrac><msup><mrow><mo></mo><mrow><msup><mi>x</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></msup><mo>-</mo><msup><mi>x</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></msup></mrow><mo></mo></mrow><mn>2</mn></msup><msub><mi>S</mi><mi>k</mi></msub></mfrac></mrow><mo>)</mo></mrow></mrow></mtd><mtd><mrow><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mo>,</mo><mi>…</mi><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>,</mo><mrow><mrow><mi>N</mi><mo>⁢</mo><mstyle><mspace width="0.8em" height="0.8ex" /></mstyle><mo>⁢</mo><mi>and</mi><mo>⁢</mo><mstyle><mspace width="0.8em" height="0.8ex" /></mstyle><mo>⁢</mo><mi>k</mi></mrow><mo>=</mo><mn>1</mn></mrow><mo>,</mo><mi>…</mi><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>,</mo><msub><mi>N</mi><mn>1</mn></msub></mrow></mtd></mtr></mtable></mtd><mtd><mrow><mo>(</mo><mn>22</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0182In eq. (22), N is the number of vectors used for testing of neural network model while N₁is the number of neurons of the hidden layer; usually, the last matches the number of samples used for neural network training. In eq. (2) S_k, which is related to Gaussian function variance, represents the so-called “spread” factor.

0183Its value is in [0, 1] range and it modulates the neuronal activation function sensitive. The smaller the parameter S_k, the more sensitive the neuron. For a better comprehension of the concept of “neurons sensitiveness”, it must kept in mind that: <br />φ_k≧0.5 <i>∀x</i>^(*)<i>εR</i>^m<i>|∥x</i>^(k)<i>−x</i>^(*)∥≦0.8326·√{square root over (<i>S</i>_k)}  (23)

0184The points x^(*)satisfying eq. (23) describe a hypersphere, having the center in x^(k), whose radius increases as S_kvalues decrease. In brief, little values of S_kinduce, for a fixed threshold of the membership probability (in the example of eq. (7) this threshold is equal to 0.5) of testing vectors to the class of the current hidden neuron (in this case it is the k-th), a larger hypersphere.

0185Known RBF-PNN neural networks have two limitations: they use the same spread factor for each neuron of hidden layer and they have not an explicit and definite procedure to determine the optimal value of S_kaccording to the current classification problem.

0186The neural network model developed by the applicants, called MultiSpread-PNN, overcomes the above noted two limitations of known models.

0187First, the hidden layer coupling each neurons is built with different spread factors. Second, an explicit and definite procedure to determine the optimal string of N₁spread factors is established. In this last phase, EA (Evolutionary Algorithms) are used.

0188<figref idref="DRAWINGS">FIG. 20</figref> compares graphically a RBF-PNN and a MultiSpread-PNN neural network of this invention. The spheres represent a three-dimensional section of hyperspheres, having unitary ray, by which the R^msubspace can be described from which the model inputs are extracted. The N₁spread factors of MultiSpread-PNN are endogenous parameters of the model. Research space of these parameters is R_[0,1]^N^{<sub2>1</sub2>}, where the symbol R_[0,1] represents real numbers in the compact range [0, 1]. An exhaustive research in a similar space is not possible, thus for this research EA, such as ES−(λ+μ) (Evolution Strategy), PSOA (Particle Swarm Optimization Algorithm) and different variants of DE (Differential Evolution algorithm) were used.

0189A trivial choice of the fitness function could be the classification error on testing data set of the MultiSpread-PNN model. In so doing, a “generalized” estimate of the endogenous parameters of model cannot be obtained. “Generalized” estimate means that the choice of the endogenous parameters of a model is made to increase model generalization capability, that is model “generalized forecast capability” [17].

0190The shape of the fitness function that was used is the following one:

0191<maths id="MATH-US-00018" num="00018"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><msub><mi>V</mi><mn>0</mn></msub><mo>⁡</mo><mrow><mo>(</mo><mrow><msub><mi>S</mi><mn>1</mn></msub><mo>,</mo><msub><mi>S</mi><mn>2</mn></msub><mo>,</mo><mi>…</mi><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>,</mo><msub><mi>S</mi><msub><mi>N</mi><mn>1</mn></msub></msub></mrow><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mfrac><mn>1</mn><msup><mi>N</mi><mo>*</mo></msup></mfrac><mo>⁢</mo><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><msup><mi>N</mi><mo>*</mo></msup></munderover><mo>⁢</mo><mrow><mfrac><mn>1</mn><mi>N</mi></mfrac><mo>⁢</mo><mrow><munderover><mo>∑</mo><mrow><msub><mi>k</mi><mi>i</mi></msub><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>⁢</mo><mrow><mo>[</mo><mrow><msub><mi>λ</mi><msub><mi>k</mi><mi>i</mi></msub></msub><mo>-</mo><mrow><msubsup><mi>F</mi><mrow><mo>(</mo><mrow><msub><mi>S</mi><mn>1</mn></msub><mo>,</mo><msub><mi>S</mi><mn>2</mn></msub><mo>,</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>…</mi><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>,</mo><msub><mi>S</mi><msub><mi>N</mi><mn>1</mn></msub></msub></mrow><mo>)</mo></mrow><mi>i</mi></msubsup><mo>⁡</mo><mrow><mo>(</mo><msup><mi>x</mi><mrow><mo>(</mo><msub><mi>k</mi><mi>i</mi></msub><mo>)</mo></mrow></msup><mo>)</mo></mrow></mrow></mrow><mo>]</mo></mrow></mrow></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>24</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mi>where</mi></mtd><mtd><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle></mtd></mtr><mtr><mtd><mrow><msup><mi>N</mi><mo>*</mo></msup><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mrow><mi>N</mi><mo>+</mo><msub><mi>N</mi><mn>1</mn></msub></mrow></mtd></mtr><mtr><mtd><mi>N</mi></mtd></mtr></mtable><mo>)</mo></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>25</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

0192The formula (24) derives from generalization of the “ordinary cross-validation estimate” of endogenous parameters of a neural network (chapter 5 of [18], and [19]). The parameter N* is the number of possible choices of a testing set with N samples in a data set composed by N+N₁input output couples. In eq. (24), k_ilabels the N elements of testing data set selected with the i-th choice. MultiSpread-PNN output is described by the symbol Fⁱ(_{S1,S2, . . . ,S}_N_{<sub2>1</sub2>}).

0193The optimal string of spread factors is the one which minimizes the functional V₀(S₁,S₂, . . . , S_N_{<sub2>1</sub2>}). To search this minimum, the above mentioned evolutionary algorithms were used.

0194It should be remarked that notwithstanding the differences between the MultiSpread-PNN and the RBF-PNN, the time spent for setting them up (that is for determining the values of their parameters) is substantially identical.

0195<figref idref="DRAWINGS">FIGS. 21 and 22</figref> show the results obtained by the best RBF-PNN and MultiSpread-PNN, respectively. The models were trained and tested on the Yamaha engine data set. These preliminary outcomes suggest that in the sensible lambda values ranges, where a satisfactory real lambda sensor is practically infallible, virtual emulator of an on-off lambda sensor according to the model of this invention fails about just one time over 100 estimates. In other words, this means that in the sensible lambda values ranges, the neural network models that have been found have a classification rate of about 99%.

0196<figref idref="DRAWINGS">FIG. 23</figref> shows a performance comparison between RBF-PNN and MultiSpread-PNN on 25 random permutations of the Yamaha data set. For each trial, model classification error was compared for, on the whole range of lambda values, related only to the testing data set. For the novel MultiSpread-PNN, the classification error averages to 4.22% while for the classical RBF-PNN, the classification error averages to 5.13%.

0197By comparing the virtual lambda sensor model of this invention with the models described in literature the following remarks can be made. The novel model of the applicants has a different working logic, it shows only one an incertitude region of relatively small width (blue rectangle in <figref idref="DRAWINGS">FIG. 17</figref>).

0198By having as inputs engine speed and inlet manifold pressure, the novel model has a forecast capability that is not limited to a single engine speed and/or a single throttle position.

0199The novel model is defined through a data pre-processing that establishes the optimal number of instantaneous cycles to be averaged for maximizing the correlation between MultiSpread-PNN model inputs and outputs.

0200Unlike neural network models known in literature to solve classification problems, the novel MultiSpread-PNN model has on average a larger forecast capability for the same set-up computational complexity.

0201A neural network model as the novel MultiSpread-PNN can be simply implemented with a low cost micro-controller. The model can be downloaded on the micro-controller memory as a sequence of matrices. Computational cost for a real-time application of the MultiSpread-PNN would be equal to the time that micro-controller spends to perform simple matrices products.

0202The only limitation of the MultiSpread-PNN model for real-time applications is related to the number of successive pressure cycles (16) over which input parameters must be averaged for maximizing the correlation between inputs and outputs. For example, in order to set-up a fuel injection control system, it must be accounted the fact that it will take 16 cycles to obtain the model inputs and the successive strategy for updating the fuel injection law as determined by the controller. In other words, the control system has a delay equivalent to 16 pressure cycles.

0203Of course, this limitation can be overcome by storing in a queue the parameters values during 16 consecutive pressure cycles. The MultiSpread-PNN model inputs would be obtained by averaging the queued values; the last value is updated by a FIFO (first in first out) strategy. By this expedient an injection law control system with a delay equal to 1 pressure cycle can be realized.

0204<figref idref="DRAWINGS">FIG. 24</figref> is a block diagram scheme of a feedforward-and-feedback control system employing a virtual lambda sensor of this invention.

0205Basically, the system is composed of a feedforward controller A, a pressure sensor, a device of this invention C for sensing the air/fuel ratio and a feedback controller B.

0206The feedforward controller A is input with signals representative of the speed and of the load of the engine, and outputs a signal DI_FF, that represents a duration of the fuel injection of the engine, and a spark generation control signal SA, that determines the spark-advance of the engine. The levels of these signals DI_FFand SA are calculated by the feedforward controller A as a function of the current speed and load of the engine using a pre-determined model of the engine.

0207The feedback controller B generates a feedback signal DI_FBfor correcting the calculated fuel injection duration, represented by the signal DI_FF, as a function of the difference between the signal generated by the virtual lambda sensor C of this invention and a reference value REF.

0000Feedforward Controller

0208This block generates the signals SA and DI_FFas a function of the speed and load of the engine by using control maps of the engine. In practice, according to a common control technique, a mathematical model of the functioning of the engine is determined during a test phase in order to determine, for any pair of values of speed and load of the engine, the best values of the duration of the fuel injection and the spark-advance of the engine.

0209The optimal duration of the fuel injection is that corresponding to the condition λ=1. In practice, the feedforward controller A compares the input values of speed and load with those stored in a look-up table generated during a test phase of the engine, and outputs the signals DI_FFand SA of corresponding values. When the input values of the speed and load do not correspond to any pair of the look-up table, the feedfoward controller A calculates the levels of the signals DI_FFand SA by linear interpolation.

0210Notably: <ul id="ul0007" list-style="none"><li id="ul0007-0001" num="0000"><ul id="ul0008" list-style="none"><li id="ul0008-0001" num="0211">the best fuel injection duration, for a pair of input values of the speed and load not contemplated in this look-up table, could differ from the value obtained by linear interpolation because an engine is a nonlinear system;</li><li id="ul0008-0002" num="0212">the test phase is carried out under pressure and temperature conditions, may differ sensibly from real working conditions of the engine.</li></ul></li></ul>

0213Therefore, the calculated duration of fuel injection represented by the signal DI_FFis corrected by a feedback signal DI_FBgenerated as a function of the sensed air/fuel ratio of the engine.

0000Feedback Controller

0214The feedback controller B generates a feedback signal for correcting the duration of fuel injection calculated by the feedforward controller A, as a function of the difference between the signal output by the virtual lambda sensor of this invention λ and a reference value REF.

0215The feedback controller B comprises an error evaluation subsystem B<b>1</b> and a correction subsystem B<b>2</b>. The error evaluation subsystem B<b>1</b> generates an error signal E<smallcaps>RROR </smallcaps>and an error variation signal C_<smallcaps>ERROR</smallcaps>, the values of which at a generic instant T are determined according to the following equations:

0216<maths id="MATH-US-00019" num="00019"><math overflow="scroll"><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mrow><mi>Error</mi><mo>⁡</mo><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><msub><mi>N</mi><mn>1</mn></msub><mo>·</mo><mrow><mo>(</mo><mrow><mi>REF</mi><mo>-</mo><mi>λ</mi></mrow><mo>)</mo></mrow></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mrow><mi>C_error</mi><mo>⁢</mo><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><msub><mi>N</mi><mn>2</mn></msub><mo>·</mo><mrow><mo>(</mo><mrow><mrow><mi>Error</mi><mo>⁡</mo><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow></mrow><mo>-</mo><mrow><mi>Error</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>T</mi><mo>-</mo><mrow><mi>Δ</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><msub><mi>T</mi><mn>1</mn></msub></mrow></mrow><mo>)</mo></mrow></mrow></mrow><mo>)</mo></mrow></mrow></mrow></mtd></mtr></mtable><mo> </mo></mrow></mrow></math></maths><br /> wherein N₁and N₂are normalization constants and ΔT₁is a time delay.

0217A sample embodiment of the error evaluation subsystem B<b>1</b> is shown in <figref idref="DRAWINGS">FIG. 25</figref>.

0218The correction subsystem B<b>2</b> is preferably composed of a correction unit CONTROLLER and an output stage, as shown in <figref idref="DRAWINGS">FIG. 26</figref>.

0219The correction unit CONTROLLER is input with the signals E<smallcaps>RROR </smallcaps>and C_<smallcaps>ERROR </smallcaps>and generates, as a function thereof, a correction signal Δ_DI of the feedback signal DI_FBin order to nullify its input signals E<smallcaps>RROR </smallcaps>and C_<smallcaps>ERROR. </smallcaps>

0220Preferably, the correction unit CONTROLLER is a fuzzy logic unit with two antecedents, that are the normalized values <smallcaps>EN </smallcaps>and <smallcaps>CEN </smallcaps>of the signals E<smallcaps>RROR </smallcaps>and C_<smallcaps>ERROR</smallcaps>, respectively, and one consequent, that is the normalized value <smallcaps>OUTPUT</smallcaps><b>1</b> of the correction signal Δ_DI, and preferably it is defined by three membership functions for each antecedent and consequent.

0221<figref idref="DRAWINGS">FIG. 27</figref> shows sample membership functions E_N, E_Z and E_P when the antecedent <smallcaps>EN </smallcaps>is negative, null (in a fuzzy sense) or is positive, respectively. Similarly, <figref idref="DRAWINGS">FIGS. 28 and 29</figref> show sample membership functions for the antecedent <smallcaps>CEN </smallcaps>and the consequent <smallcaps>OUTPUT</smallcaps><b>1</b>, respectively.

0222The fuzzy correction unit CONTROLLER generates the correction signal Δ_DI according to the nine fuzzy rules shown in <figref idref="DRAWINGS">FIG. 30</figref>. These fuzzy rules are graphically illustrated in <figref idref="DRAWINGS">FIG. 31</figref> and in the tree-dimensional graph of <figref idref="DRAWINGS">FIG. 32</figref>.

0223The output stage of the correction subsystem B<b>2</b> is preferably composed of an amplifier N₃of the correction signal Δ_DI, and of a positive feedback loop that generates the feedback signal DI_FBby adding to the amplified correction signal Δ_DI a delayed replica thereof Δ_DI(T−ΔT₂) of a certain time delay ΔT₂: <br /><i>ID</i>_FB<i>=N</i>₃·Δ_—<i>DI</i>(<i>T</i>)+Δ_—<i>DI</i>(<i>T−ΔT</i>)<br /> Learning Machine

0224The learning machine C includes an identification subsystem C<b>2</b> that chooses the smallest number of characteristic parameters of the detected pressure signal sufficient for estimating the value of the air/fuel ratio. The subsystem C<b>2</b> is input with the pressure signal generated by the pressure sensor in contact with at least a cylinder of the engine and implements a clustering algorithm for choosing moving averages of these characteristic parameters.

0225It is important that the number of characteristic parameters to be considered be as small as possible for reducing memory requirements and the number of calculations to be performed for estimating the air/fuel ratio. By contrast, a too small number of parameters may degrade accuracy.

0226In practice, the identification subsystem C<b>2</b> generates data sets composed of values of moving averages of characteristic parameters of a certain initial set of parameters, that are potentially useful for evaluating the air/fuel ratio, each for a respective number of pressure cycles. As a function of the desired number of clusters, it groups in clusters the moving averages of each data set with a respective execution of a clustering algorithm. Then, the number of pressure cycles on which these moving averages are calculated is chosen as the number corresponding to the execution of the clustering algorithm for which the ratio between the clustering performance index and the ideal clustering performance index is maximum.

0227For these operations any clustering method available in literature may be used.

0228It may be remarked that the clustering method thus allows to choose the best number of pressure cycles for averaging these characteristic parameters in order to evaluate the lambda factor, as schematically illustrated in <figref idref="DRAWINGS">FIG. 33</figref>.

0229The core C<b>1</b> of the lambda sensor of this invention is input with the parameters chosen by the identification subsystem C<b>2</b>. The selected characteristic parameters of the pressure signal are the pressure Pratio<b>40</b> in a cylinder for a 40° rotation of the crank, the pressure Pratio<b>50</b> in a cylinder for a 50° rotation of the crank, and the pressure peak Pmax independently from the position of the crank at which it is attained.

0230Indeed, the core of the lambda sensor may either be a neural network, a stochastic machine, a support vector machine, a committee machine or a hybrid learning machine.

0231However, as schematically illustrated in <figref idref="DRAWINGS">FIG. 34</figref>, the most performing cores are composed of a neural network or a hybrid learning machine.

0000Virtual Lambda Sensor Core Based on a Neural Network

0232<figref idref="DRAWINGS">FIG. 35</figref> shows a block scheme of a lambda sensor core based on a neural network.

0233As already remarked, the preferred inputs of the core are the parameters Pratio<b>40</b>, Pratio<b>50</b> and Pmax, that are the characteristic parameters of the pressure signal that more often have been chosen by the identification subsystem C<b>2</b> during the extensive tests that were carried out.

0234Signals representative of these parameters are processed by an input stage E<smallcaps>DGE </smallcaps>D<smallcaps>ETECTOR</smallcaps>, for sampling the input signals and generating a synchronization pulse for the subsequent blocks in cascade.

0235According to a preferred embodiment of this invention, the instant at which the synchronization pulse is generated is determined as a function of the instant at which a pressure peak is detected. Indeed, it has been experimentally verified that this instant is sufficiently stable and relatively free of spurious variations.

0236The block P<smallcaps>RE</smallcaps>-P<smallcaps>ROC </smallcaps>is a pre-processor that generates a moving average of the data output by the input stage E<smallcaps>DGE </smallcaps>D<smallcaps>ETECTOR</smallcaps>, for filtering these data from noise.

0237The number of samples to be considered for calculating the moving average, which is a number of pressure cycles of the engine, must be established. This number should not be too large, because this would make the virtual lambda sensor less prompt in reacting to functioning variations of the engine, nor too small, otherwise the level of noise corrupting the moving average would be too high.

0238The block P<smallcaps>OST</smallcaps>-P<smallcaps>ROC </smallcaps>is a post-processor analogous to the pre-processor P<smallcaps>RE</smallcaps>-P<smallcaps>ROC. </smallcaps>

0239The pre-processor and the post-processor effectively reduce the noise corrupting the output signal of the virtual lambda sensor core, as may be easily inferred by comparing the graphs in <figref idref="DRAWINGS">FIGS. 36</figref><i>a </i>and <b>36</b><i>b </i>of the output signal with and without the blocks P<smallcaps>RE</smallcaps>-P<smallcaps>ROC </smallcaps>and P<smallcaps>OST</smallcaps>-P<smallcaps>ROC</smallcaps>, respectively.

0240The learning part N<smallcaps>EURAL </smallcaps>N<smallcaps>ETWORK </smallcaps>of the virtual lambda sensor core is a neural network, preferably a MLP (Multi Layer Perceptron) with three inputs, thirty neurons in the hidden layer and a single output. The number of neurons has been chosen by a so-called “ordinary cross-validation” procedure with single-objective and multi-objectives optimization.

0241The neural network may be trained with classic learning algorithms, such as the “resilient back propagation” algorithm and/or the Levenberg-Marquardt algorithm, and stochastic search algorithm, such as the Particle Swarm Optimization Algorithm (PSOA).

0000Virtual Lambda Sensor Core Based on a Hybrid Learning Machine

0242An embodiment of the core of the virtual lambda sensor that includes a so-called “committee machine” is shown in <figref idref="DRAWINGS">FIG. 37</figref>. It includes an input stage E<smallcaps>DGE </smallcaps>D<smallcaps>ETECTOR</smallcaps>, a pre-processor P<smallcaps>RE</smallcaps>-P<smallcaps>ROC </smallcaps>and three post-processors P<smallcaps>OST</smallcaps>-P<smallcaps>ROC </smallcaps>that are identical to the homologous blocks of <figref idref="DRAWINGS">FIG. 35</figref>.

0243Differently from the embodiment of <figref idref="DRAWINGS">FIG. 35</figref>, the embodiment of <figref idref="DRAWINGS">FIG. 37</figref> comprises three fuzzy subsystems L<smallcaps>EAN </smallcaps>FIS, <smallcaps>ON</smallcaps>-<smallcaps>OFF </smallcaps>FIS and R<smallcaps>ICH </smallcaps>FIS, that are typically “TSK Singleton” fuzzy inference subsystems, that are defined by different parameters and generate respective signals related to the air/fuel ratio of the engine for λ<1, λ>1 and λ≈1, respectively.

0244In practice, the fuzzy subsystem <smallcaps>ON</smallcaps>-<smallcaps>OFF </smallcaps>FIS works like an on/off lambda sensor. An output multiplexer S<smallcaps>WITCH </smallcaps>generate an output signal <smallcaps>OUTPUT</smallcaps><b>2</b> that is the signal generated by the subsystem L<smallcaps>EAN </smallcaps>FIS or R<smallcaps>ICH </smallcaps>FIS depending on the value of the output of the <smallcaps>ON</smallcaps>-<smallcaps>OFF </smallcaps>FIS subsystem. Even this last signal is made available (<smallcaps>OUTPUT</smallcaps><b>1</b>), thus the virtual lambda sensor core of <figref idref="DRAWINGS">FIG. 37</figref> may function either as a linear lambda sensor (<smallcaps>OUTPUT</smallcaps><b>2</b>) or as an on/off lambda sensor (<smallcaps>OUTPUT</smallcaps><b>1</b>)

0245Each fuzzy subsystem has three antecedents and one consequent and preferably is defined by three membership functions for each antecedent or consequent.

0246The fuzzy subsystems are trained using experimental data with a “supervised training” procedure implementing a stochastic search algorithm (such as a PSOA) for calculating optimal values to be assigned to the parameters of the membership functions, the crisp values and the like.

0247The fuzzy subsystem R<smallcaps>ICH </smallcaps>FIS is trained such that the output signal preferably ranges between 0.85 and 0.98, while the fuzzy subsystem L<smallcaps>EAN </smallcaps>FIS outputs a signal that ranges between 1.02 and 1.15. The fuzzy subsystem <smallcaps>ON</smallcaps>-<smallcaps>OFF </smallcaps>FIS is trained for estimating very accurately the air/fuel ratio for 0.98&lt;λ&lt;1.02. An eventual error in this range would make the output multiplexer S<smallcaps>WITCH </smallcaps>select the output of the subsystem L<smallcaps>EAN </smallcaps>FIS when the air/fuel mixture is reach or vice versa.

0248Anyway, such an error does not worsen significantly the accuracy of the device of this invention because the subsystems L<smallcaps>EAN </smallcaps>FIS and R<smallcaps>ICH </smallcaps>FIS have almost the same accuracy in the range 0.98&lt;λ&lt;1.02. Obviously, it is important that the subsystem <smallcaps>ON</smallcaps>-<smallcaps>OFF </smallcaps>FIS do not make an error larger than 0.2 in determining λ, otherwise the accuracy could be negatively affected.

0249The pre-processor and the post-processors filter the noise that may corrupt signals input to and output from the fuzzy subsystems. They are useful because spurious spikes corrupting the signal output by the subsystems could degrade relevantly the accuracy of the sensing of the air/fuel ratio. In particular, spurious spikes output by the <smallcaps>ON</smallcaps>-<smallcaps>OFF </smallcaps>FIS subsystem must be filtered because they could induce spurious switching of the multiplexer.

0250Compared with classic lambda sensors the virtual device of this invention for sensing the air/fuel ratio has numerous advantages. The device of this invention need not warming up for functioning correctly, has a relatively low cost and tracks the air/fuel ratio very quickly.

0251By contrast, certain lambda sensors (such as the lambda sensor HEGO) must attain a temperature of about 300° C. before starting to function accurately, they are relatively expensive and subjected to wear. Moreover their accuracy is limited by the fact that they are installed in the exhaust gas pipe of the engine, which means that they cannot generate a signal representing the air/fuel ratio before the exhaust gases reach the sensor. Therefore, classic lambda sensors are generally sluggish in responding to rapid changes of the functioning conditions of the engine.

0252Formula 1 cars often damage their exhaust gases pipes where the lambda sensors are installed. In these situations, a classic lambda sensor may not sense correctly the air/fuel ratio and the engine risks to be miscontrolled.

0253The feedforward-and-feedback control system of this invention has been real-time tested on the Yamaha engine YP125.

0254In order to determine the relationship between air/fuel ratio (lambda value) and parameters from the cylinder pressure cycle, STMicroelectronics and Yamaha agreed on the following conditions for the experimental tests:

02551. 4600 rpm, torque=1.5 Nm

02562. 5600 rpm, torque=4.4 Nm

02573. 4600 rpm, WOT(Wide Open Throttle) condition

0258Yamaha constraints on these engine conditions were respectively to control the maintain engine close the stoichiometric condition with a maximum 1% error and to have a response time of the control system equal or of less than 100 milliseconds from the moment the engine reaches the desired steady state.

0259The tests have been conducted maintaining, for each condition, fixed spark advance, throttle position, injection timing and modifying only the duration of the fuel injection.

0260The goal of realizing an efficient injection control system for the Yamaha YP125 mono-cylinder gasoline engine meant to maintain the YP125 engine close to stoichiometric combustion in all the above mentioned three conditions of operation.

0261A closed lop injection control system based on soft computing models was realized. The loop included a Virtual Lambda Sensor and control system of this invention.

0262<figref idref="DRAWINGS">FIG. 38</figref> describes the time-changing trend of the lambda values as sensed by a real lambda sensor when the control system of this invention was activated, in order to test its performances. The engine is effectively maintained close to the stoichiometric condition with an error smaller than 1%.

0263During the tests, the engine working conditions were changed with satisfactory results even at lower engine speeds than 4600 rpm (down to 3600 rpm) and even for different throttle positions (down to 33% opening). These results are reported in <figref idref="DRAWINGS">FIGS. 39</figref> and <b>40</b>.

0264Finally, the control system of this invention was tested under different transient conditions. After few seconds, the control system brought back the engine to stoichiometric combustion conditions within 1% error.

0265<figref idref="DRAWINGS">FIGS. 41</figref>, <b>42</b> and <b>43</b> illustrate the results of these tests under transient conditions.

0266<figref idref="DRAWINGS">FIG. 44</figref> depicts the test console for simulating the functioning of the engine controlled by the system of this invention at 4600 rpm and WOT condition.

0267<figref idref="DRAWINGS">FIGS. 45 and 46</figref> depict the test console when the feedforward-and-feedback control system of the engine had restored stoichiometric conditions (λ≈1) starting from a condition in which the engine in working with a rich mixture (λ<1) and with a lean mixture (λ>1), respectively.

0268The meaning of the labels in <figref idref="DRAWINGS">FIGS. 44 to 46</figref> is resumed in the following table:

0269<tables id="TABLE-US-00003" num="00003"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="3"><colspec colname="offset" colwidth="21pt" align="left" /><colspec colname="1" colwidth="91pt" align="left" /><colspec colname="2" colwidth="105pt" align="left" /><thead><row><entry /><entry namest="offset" nameend="2" align="center" rowsep="1" /></row></thead><tbody valign="top"><row><entry /><entry>Coppia (Nm)</entry><entry>Torque</entry></row><row><entry /><entry>P (kW)</entry><entry>Power</entry></row><row><entry /><entry>Bst_Map (bar)</entry><entry>Intake manifold pressure</entry></row><row><entry /><entry>ALPHA (%)</entry><entry>Throttle position</entry></row><row><entry /><entry>DEG_DGMT (kg/h)</entry><entry>Intake manifold air flow</entry></row><row><entry /><entry>T_AIR (° C.)</entry><entry>Air temperature</entry></row><row><entry /><entry>FB_TEMP (° C.)</entry><entry>Fuel balance temperature</entry></row><row><entry /><entry>T_ACQ_US (° C.)</entry><entry>Cooling system water</entry></row><row><entry /><entry /><entry>temperature</entry></row><row><entry /><entry>T_ASP1 (° C.)</entry><entry>Intake manifold</entry></row><row><entry /><entry /><entry>temperature</entry></row><row><entry /><entry>T_ASP2 (° C.)</entry><entry>Intake manifold</entry></row><row><entry /><entry /><entry>temperature</entry></row><row><entry /><entry>T_SCARI1 (° C.)</entry><entry>Exhaust gases temperature</entry></row><row><entry /><entry>T_SCARI2 (° C.)</entry><entry>Exhaust gases temperature</entry></row><row><entry /><entry>LAMBDA 1</entry><entry>Lambda value</entry></row><row><entry /><entry>FB_VAL (kg/h)</entry><entry>Fuel balance value</entry></row><row><entry /><entry namest="offset" nameend="2" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

REFERENCES

0000<ul id="ul0009" list-style="none"><li id="ul0009-0001" num="0270">[1] T. W. Long, C. M. Atkinson and E. L. Hanzevack. Virtual sensing: a neural-network-based intelligent performance and emissions prediction system for on-board diagnostics and engine control. In Proceedings of the 1998 SAE International Congress &amp; Exposition, vol. 1357, 208, pages 39 {51, 1998.</li><li id="ul0009-0002" num="0271">[2] P. A. Howson, R. J. Howlett, S. D. Walters and I. Park. Air-fuel ratio measurement in an internal combustion engine using a neural network. In International Conference on Advances in Vehicle Control and Safety, AVCS'98, Amiens, France, July 1998.</li><li id="ul0009-0003" num="0272">[3] S. D. Walters, R. J. Howlett, M. M. de Zoysa and P. A. Howson. Neural network techniques for monitoring and control of internal combustion engines. In International ICSC Symposium on Intelligent Industrial Automation, Genoa, Italy, June 1999.</li><li id="ul0009-0004" num="0273">[4] A. M. Wahl, K. Ratton, E. N. Balles, E. A. VanDyne and M. C. Lai. In cylinder air/fuel ratio approximation using spark gap ionization sensing. In Proceedings of the 1998 SAE International Congress &amp; Exposition, vol. 1356, pages 39-44, 1998.</li><li id="ul0009-0005" num="0274">[5] F. Tanabe, Y. Ohashi, W. Fukui and A. Ueda. The application of ionic current detection system for the combustion limit control. In Proceedings of the 1998 SAE International Congress &amp; Exposition, vol. 1356, pages 79-85, 1998.</li><li id="ul0009-0006" num="0275">[6] T. Rognvaldsson, N. Wickstrom, C. Carlsson, M. Larsson, M. Hellring, T. Munther and J. Nytomt. Spark advance control using the ion current and neural soft sensor. SAE Paper, 99P-78.</li><li id="ul0009-0007" num="0276">[7] T. Rognvaldsson, N. Wickstrom, C. Carlsson, M. Larsson, M. Hellring, T. Munther and J. Nytomt. Robust afr estimation using the ion current and the neural networks. SAE Paper, 99P-76.</li><li id="ul0009-0008" num="0277">[8] R. J. Howlett and S. D. Walters. A multi-computer neural network architecture. IEEE Electronics Letters, 35.</li><li id="ul0009-0009" num="0278">[9] M. M. de Zoysa, R. J. Howlett and S. D. Walters. Monitoring internal combustion engines by neural network based on virtual sensing. In Recent Advances in Artificial Neural Networks Design and Applications, pages 291-318. Jain L. and Fanelli A. M., CRC Press, 2000.</li><li id="ul0009-0010" num="0279">[10] R. J. Howlett, M. M. de Zoysa, S. D. Walters and P. A. Howson, Neural Network Techniques for Monitoring and Control of Internal Combustion Engines, International Symposium on Intelligent Industrial Automation 1999.</li><li id="ul0009-0011" num="0280">[11] A. N. Tikhonov. On solving incorrectly posed problems and method of regularization. <i>Doklady Akademii Nauk USSR, </i>151.</li><li id="ul0009-0012" num="0281">[12] A. N. Tikhonov. On regularization of ill-posed problems. <i>Doklady Akademii Nauk USSR, </i>153.</li><li id="ul0009-0013" num="0282">[13] A. N. Tikhonov and V.Y. Arsenin. Solution of Ill-posed Problems. W. H. Winston.</li><li id="ul0009-0014" num="0283">[14] M. J. D. Powell. Radial basis function for multivariate interpolation: A review. IMA <i>Conference on Algorithms for the Approximations of Functions and Data. </i></li><li id="ul0009-0015" num="0284">[15] M. J. D. Powell. The theory of radial basis function approximation in 1990. <i>Advances in Numerical Analysis </i>Vol. II: Wavelets, Subdivision Algorithms, and Radial Basis Functions.</li><li id="ul0009-0016" num="0285">[16] M. J. D. Powell. Radial basis function approximations to polynomials. In <i>Numerical Analysis </i>1987 <i>Proceedings</i>, Dundee, UK, pages 223-241, 1988.</li><li id="ul0009-0017" num="0286">[17] V. N. Vapnik. An overview of statistical learning theory. <i>IEEE Trans. On Neural Networks, </i>9, 1999.</li><li id="ul0009-0018" num="0287">[18] S. Haykin. <i>Neural Networks: A comprehensive Foundation</i>. Prentice Hall International Editions, Piscataway, 1999.</li><li id="ul0009-0019" num="0288">[19] P. Craven and G. Wahba. Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation. <i>Numerische Mathematik, </i>31.</li><li id="ul0009-0020" num="0289">[20] P. Yoon, S. Park, M. Sunwoo, I. Ohm and K. J. Yoon, <i>Closed Loop Control of Spark Advance and Air</i>-<i>Fuel Ratio in SI Engines Using Cylinder Pressure</i>, SAE 2000-01-0933</li><li id="ul0009-0021" num="0290">[21] S. Coffa and L. Occhipinti, Monolithically integrated pressure sensor and related process for manufacturing, U.S. Pat. No. 6,743,654 and EP1215476,</li><li id="ul0009-0022" num="0291">[22] N. R. Pal, K. Pal and J. C. Bezdek, <i>A Mixed c</i>-<i>Means Clustering Model</i>, Proceedings of IEEE International Conference on Fuzzy Systems, pp. 11-21, 1997.</li><li id="ul0009-0023" num="0292">[23] R. Maitra, <i>Clustering Massive Datasets with Applications in Software Metrics and Tomography</i>, J. Amer. Stat. Assoc., 1998.</li><li id="ul0009-0024" num="0293">[24] H. Frigui and R. Krishnapuram, <i>A Possibilistic approach to Clustering</i>, IEEE Transactions on Fuzzy Systems, vol. 1, No. 2, pp. 98-110, 1993.</li><li id="ul0009-0025" num="0294">[25] J. M. Leski, <i>An</i>-<i>Insensitive Fuzzy C</i>-<i>Means Clustering</i>, Int. J. Appl. Math. Comp. Sci., vol. 11, No. 4, pp. 993-1007, 2001.</li><li id="ul0009-0026" num="0295">[26] H. Frigui and R. Krishnapuram, <i>A Robust Competitive Clustering Algorithm With Applications in Computer Vision</i>, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 21, No. 5, pp. 450-465, 1999.</li></ul>