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

[0001] This 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.

[0002] In 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.

[0003] To 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] FIG. 1 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 FIG. 2.

[0005] In 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 FIGS. 1 and 2), in FIG. 2 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 FIG. 3, 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.

[0006] To 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] $\lambda =\frac{\mathrm{Air}\text{/}\mathrm{Fuel}}{\mathrm{Air}\text{/}{\mathrm{Fuel}}_{\mathrm{stoichiometric}}}$
If λ<1 the mixture is rich of fuel, while if λ>1 the mixture is lean of fuel, as schematically shown in FIG. 4.

[0008] The 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.

[0009] Many 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.

[0010] There 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]

	accuracy	accuracy for	accuracy
	for lean	stoichiometric	for rich	cost
	mixtures	mixtures	mixtures	(USD)
McLaren	1.7%	0.1%	1.7%	1500-1800
electronic systems
MoTeC	2.5%	1.75%	1.75%	800-900
Bosch LSM 11	1.5%	unknown	unknown	300-400
Horiba LD-700	8.0%	4.0%	8.0%	60-80

[0012] Engines 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].

[0013] An 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.

[0014] There 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].

[0015] The 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] FIG. 5 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.

[0017] Interactions 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.

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

AFR Estimation Using Spark Voltage Characterization by Neural Network

[0019] According 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.

[0020] Such 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] FIG. 6 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.

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

[0023] Air-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.

[0024] In 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.

[0025] In a supervised training paradigm, a back-propagation learning algorithm sets the MLP training. The training file contains N_t pairs input-output; model input is an instantaneous spark-voltage vector of the form V_i=(ν₁, ν₂, . . . , ν_m), with i=1, . . . , N_t and 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).

[0026] Three 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.

[0027] R. 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. FIG. 7 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_t is the size of training file and N_w is the number of weights of the MLP neural network modeling virtual lambda sensor.

[0028] These 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 FIG. 8.

[0029] The 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.

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

[0031] It 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.

[0032] The 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.

[0033] In [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.

[0034] Unfortunately, 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.

[0035] As 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.

[0036] The 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.

[0037] Therefore, 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.

[0038] It 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.

[0039] The 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.

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

[0041] The 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.

[0042] According 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.

[0043] The 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.

[0086] According 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 FIG. 9.

[0087] Of 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.

[0088] Indeed, 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.

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

[0090]

	unit of
	measure	element description
Speed	rpm	engine speed
lambda	[ ]	lambda values
Aircycle	mg/cycle	air massive flow
BstMap	bar	intake manifold pressure
BurDur	deg	combustion duration
pEVC	bar	pressure at exhaust valve closure
pEVO	bar	pressure at exhaust valve opening
pIVC	bar	pressure at intake valve closure
pIVO	bar	pressure at intake valve opening
Pratio40	[ ]	pressure ratio between pressures
		at 40 crank angles before and
		after TDC
Pratio50	[ ]	pressure ratio between pressures
		at 50 crank angles before and
		after TDC
Pratio60	[ ]	pressure ratio between pressures
		at 60 crank angles before and
		after TDC
Pratio70	[ ]	pressure ratio between pressures
		at 70 crank angles before and
		after TDC
Pratio80	[ ]	pressure ratio between pressures
		at 80 crank angles before and
		after TDC
Pratio90	[ ]	pressure ratio between pressures
		at 90 crank angles before and
		after TDC
Pratio100	[ ]	pressure ratio between pressures
		at 100 crank angles before and
		after TDC
Pratio110	[ ]	pressure ratio between pressures
		at 110 crank angles before and
		after TDC
Pmax	bar	maximum pressure
PcompMax	[ ]	ratio between maximum of pressure
		cycle and maximum of pressure
		cycle without combustion

[0091] These 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.

[0092] Given 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.

[0093] A 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.

[0094] Of 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.

[0095] The three parameters most correlated with the air/fuel ratio resulted to be Pratio40, Pratio50 and Pmax and, according to a preferred embodiment of this invention, these three parameters were used as the inputs of the learning machine.

[0096] In 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.

[0097] A 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.

[0098] In 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] FIGS. 10a to 10d show the results of a clustering process carried out among moving averages of the parameter Pratio40 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 Pratio40 and the corresponding λ value.

[0100] In 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]).

[0101] A 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.

Clustering Factor of M-Dimensional Data

[0102] Clustering 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.

[0103] There 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.

[0104] There 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?

[0105] Let 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<c<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_ik being k=1, . . . , N and i=1, . . . , c.

[0106] U_ik elements are suitably linked to the distances of X points from c temporary cluster centers.

[0107] $\begin{array}{cc}{M}_P=\{0\le U_{\mathrm{ik}}\le 1\forall i,k\forall k\exists iU_{\mathrm{ik}}>0\}& \left(1\right)\\ M_F=\left\{\sum _{i=1}^c{U}_{\mathrm{ik}}=1\forall k\sum _{k=1}^N{U}_{\mathrm{ik}}>0\forall i\right\}& \left(2\right)\\ M_C=\left\{U_{\mathrm{ik}}=0\mathrm{or}1\forall i,k\right\}& \left(3\right)\end{array}$

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

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

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

[0111] Moreover, 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.

[0112] To 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.

[0113] Even 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.

[0114] Two 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.

Fuzzy C-Means Algorithm (FCM)

[0115] The 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_k and for each algorithm step, the sum of membership degrees must be equal to 1.

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

[0117] $\begin{array}{cc}{J}_m\left(U,V\right)=\sum _{i=1}^c\sum _{k=1}^N{U}_{\mathrm{ik}}^m{D}_{\mathrm{ik}}^2& \left(5\right)\end{array}$

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

[0119] $\begin{array}{cc}\sum _{i=1}^c{U}_{\mathrm{ik}}=1\forall k=1,\dots ,N& \left(6\right)\end{array}$

[0120] In eq. (5) m is the system fuzzyness degree while D_ik matrix 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.

[0121] The FCM procedure can be executed for several clustering strategies of data set X, that is, for different number of clusters (from c_min to 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.

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

[0123] $\begin{array}{cc}P\left(c\right)={\tilde{J}}_{\min }\left(c\right)-\sum _{i=1}^c\sum _{k=1}^N{U}_{\mathrm{ik}}^m·\stackrel{\_}{x}-v_i& \left(7\right)\end{array}$

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

[0125] P(c) has a minimum when data set clustering has a minimum intra-cluster variance (i.e. small values of D_ik in e {tilde over (J)}_min(c)) and a maximum inter-cluster variance (i.e. maximum cluster center distances v_i from data set centers x). Hence, a graph of the performance index versus the number of clusters may be plot. FIG. 11 shows a given data set grouped with different numbers of clusters (c=2, . . . , 8).

[0126] Even looking at FIG. 11 it is possible to see that the best data set clustering is obtained with four clusters.

Fuzzy Possibilistic C-Means Algorithm (FPCM)

[0127] The 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].

[0128] In 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_k has 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.

[0129] For the FCM, the observer is, for every time, in x_k point. 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.

[0130] For the FPCM, the observer is not, for every time, only in x_k point 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] $\begin{array}{cc}\sum _{k=1}^N{T}_{\mathrm{ik}}=1\forall =1,\dots ,c& \left(8\right)\end{array}$

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

[0133] $\begin{array}{cc}{J}_{m,\eta }\left(U,T,V\right)=\sum _{i=1}^c\sum _{k=1}^N\left(u_{\mathrm{ik}}^m+t_{\mathrm{ik}}^{\eta }\right)D_{\mathrm{ik}}^2& \left(9\right)\end{array}$

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

[0135] $\begin{array}{cc}\sum _{i=1}^c{U}_{\mathrm{ik}}=1\forall k& \left(10\right)\\ \sum _{k=1}^N{T}_{\mathrm{ik}}=1\forall i& \left(11\right)\end{array}$

[0136] In 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.

Data Set Clustering Factor

[0137] The 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).

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

[0139] $\begin{array}{cc}\left(\begin{array}{ccccccc}0& 1& 0& 1& 0& \cdots & 0\\ 1& 0& 0& 0& 0& \cdots & 0\\ 0& 0& 0& 0& 1& \cdots & 1\\ ⋮& ⋮& ⋮& ⋮& ⋮& \cdots & ⋮\\ 0& 0& 1& 0& 0& \cdots & 0\end{array}\right)& \left(12\right)\\ \left(\begin{array}{ccccccc}0& \frac{1}{n_1}& 0& \frac{1}{n_1}& 0& \cdots & 0\\ \frac{1}{n_2}& 0& 0& 0& 0& \cdots & 0\\ 0& 0& 0& 0& \frac{1}{n_3}& \cdots & \frac{1}{n_3}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \frac{1}{n_c}& 0& 0& \cdots & 0\end{array}\right)& \left(13\right)\end{array}$

[0140] Data set points n_i falling on c clusters must satisfy the constraint:

[0141] $\begin{array}{cc}\sum _{i=1}^c{n}_i=N& \left(14\right)\end{array}$

[0142] The 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] $\begin{array}{cc}{n}_1=n_2=\dots =n_c=\frac{N}{c}& \left(15\right)\end{array}$

[0144] Considering 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] $\begin{array}{cc}{P}_{\mathrm{asint}.}^{\mathrm{FCM}}\left(c\right)=-\sum _{k=1}^N{\stackrel{\_}{x}-v_i}^2& \left(16\right)\\ P_{\mathrm{asint}}^{\mathrm{FPCM}}\left(c\right)=-\left(\frac{N}{c}+{\left(\frac{c}{N}\right)}^{\eta -1}\right)·\sum _{k=1}^N{\stackrel{\_}{x}-v_i}^2& \left(16\right)\end{array}$

[0146] To 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.

Data Set Pre-Processing

[0147] The 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.

[0148] The clustering process between Pratio40 (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 Pratio40 and λ increases strongly.

[0149] The correlation between Pratio40 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”.

[0150] Having fixed the model input (as Pratio40) 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.

[0151] Several cluster center initializations have been used. FIGS. 12 and 13 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.

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

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

[0154] Therefore, the algorithm input data should be previously normalized in [0, 1]. FIGS. 14 and 15 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 (λ, Pratio40 and Pratio50 values) have values within different ranges, in fact Pratio40 and Pratio50 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 Pratio40 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).

[0155] There 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.

[0156] In 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.

A Learning Machine

[0157] According 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.

[0158] As depicted in FIGS. 16 and 17, 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.

[0159] Differently 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.

[0160] In 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 FIG. 17) 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.

[0161] From a mathematical viewpoint, a neural network with a scalar output can be described as an hyper-surface Γ in R^m+1 space, 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:
F(x⁽ⁱ⁾)=d_i, i=1, 2, . . . , N  (17)

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

[0163] $\begin{array}{cc}F\left(x\right)=\sum _{i=1}^N{w}_i\varphi \left(x-x^{\left(i\right)}\right)& \left(18\right)\end{array}$
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. FIG. 18 depicts the architecture of the neural network RBF described by M. J. D. Powell equation (18).

[0164] Optimal values of the coefficients w_i are determined by interpolation conditions, that is:

[0165] $\begin{array}{cc}\left[\begin{array}{cccc}{\varphi }_{11}& {\varphi }_{12}& \dots & {\varphi }_{1N}\\ {\varphi }_{21}& {\varphi }_{22}& \dots & {\varphi }_{2N}\\ ⋮& ⋮& ⋰& ⋮\\ {\varphi }_{N1}& {\varphi }_{N2}& \dots & {\varphi }_{\mathrm{NN}}\end{array}\right]·\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_N\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\\ ⋮\\ d_N\end{array}\right]& \left(16\right)\end{array}$
where φ_ij=φ(∥x^(j)−x⁽ⁱ⁾∥). We may rewrite in matrix form the previous equation:
Φ·w=d  (20)

[0166] Assuming that interpolation matrix Φ is nonsingular, we have that:
w=d·Φ⁻¹  (21)

[0167] In literature, there are several classes of functions for which interpolation matrix Φ is always invertible:

• • [0168] Multiquadrics: φ(r)=√{square root over (r²+c²)} for some c>0 and rεR
  • [0169] Inverse multiquadrics:

[0170] $\varphi \left(r\right)=\frac{1}{\sqrt{r^2+c^2}}$

• • [0171]  for some c>0 and rεR
  • [0172] Gaussian functions:

[0173] $\varphi \left(r\right)=\exp \left(-\frac{r^2}{2{\sigma }^2}\right)$

• • [0174]  for some σ>0 and rεR

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

[0176] There 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.

[0177] From 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<c<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.

[0178] Matrix element U_ik represents the probability that the i-th point belongs to c-th cluster. Usually, matrix elements U_ik satisfy some normalization conditions.

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

[0180] FIG. 19 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_ik matrix. 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] $\begin{array}{cc}\begin{array}{cc}{\varphi }_k=\exp \left(-\frac{{x^{\left(k\right)}-x^{\left(i\right)}}^2}{S_k}\right)& i=1,\dots ,N\mathrm{and}k=1,\dots ,N_1\end{array}& \left(22\right)\end{array}$

[0182] In 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.

[0183] Its 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:
φ_k≧0.5 ∀x^(*)εR^m|∥x^(k)−x^(*)∥≦0.8326·√{square root over (S_k)}  (23)

[0184] The points x^(*) satisfying eq. (23) describe a hypersphere, having the center in x^(k), whose radius increases as S_k values decrease. In brief, little values of S_k induce, 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.

[0185] Known 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_k according to the current classification problem.

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

[0187] First, 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] FIG. 20 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^m subspace 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]^N1, 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.

[0189] A 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].

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

[0191] $\begin{array}{cc}{V}_0\left(S_1,S_2,\dots ,S_{N_1}\right)=\frac{1}{N^{*}}\sum _{i=1}^{N^{*}}\frac{1}{N}\sum _{k_i=1}^N\left[{\lambda }_{k_i}-F_{\left(S_1,S_2,\dots ,S_{N_1}\right)}^i\left(x^{\left(k_i\right)}\right)\right]& \left(24\right)\\ \mathrm{where}& \\ N^{*}=\left(\begin{array}{c}N+N_1\\ N\end{array}\right)& \left(25\right)\end{array}$

[0192] The 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_i labels the N elements of testing data set selected with the i-th choice. MultiSpread-PNN output is described by the symbol Fⁱ(_{S1,S2, . . . ,SN1}).

[0193] The optimal string of spread factors is the one which minimizes the functional V₀(S₁,S₂, . . . , S_N1). To search this minimum, the above mentioned evolutionary algorithms were used.

[0194] It 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] FIGS. 21 and 22 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] FIG. 23 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%.

[0197] By 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 FIG. 17).

[0198] By 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.

[0199] The 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.

[0200] Unlike 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.

[0201] A 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.

[0202] The 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.

[0203] Of 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] FIG. 24 is a block diagram scheme of a feedforward-and-feedback control system employing a virtual lambda sensor of this invention.

[0205] Basically, 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.

[0206] The 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_FF and 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.

[0207] The feedback controller B generates a feedback signal DI_FB for 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.

Feedforward Controller

[0208] This block generates the signals SA and DI_FF as 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.

[0209] The 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_FF and 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_FF and SA by linear interpolation.

[0210] Notably:

• • [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;
  • [0212] the test phase is carried out under pressure and temperature conditions, may differ sensibly from real working conditions of the engine.

[0213] Therefore, the calculated duration of fuel injection represented by the signal DI_FF is corrected by a feedback signal DI_FB generated as a function of the sensed air/fuel ratio of the engine.

Feedback Controller

[0214] The 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.

[0215] The feedback controller B comprises an error evaluation subsystem B1 and a correction subsystem B2. The error evaluation subsystem B1 generates an error signal ERROR and an error variation signal C_ERROR, the values of which at a generic instant T are determined according to the following equations:

[0216] $\{\begin{array}{c}\mathrm{Error}\left(T\right)=N_1·\left(\mathrm{REF}-\lambda \right)\\ \mathrm{C\_error}\left(T\right)=N_2·\left(\mathrm{Error}\left(T\right)-\mathrm{Error}\left(T-\Delta T_1\right)\right)\end{array}\hspace{1em}$
wherein N₁ and N₂ are normalization constants and ΔT₁ is a time delay.

[0217] A sample embodiment of the error evaluation subsystem B1 is shown in FIG. 25.

[0218] The correction subsystem B2 is preferably composed of a correction unit CONTROLLER and an output stage, as shown in FIG. 26.

[0219] The correction unit CONTROLLER is input with the signals ERROR and C_ERROR and generates, as a function thereof, a correction signal Δ_DI of the feedback signal DI_FB in order to nullify its input signals ERROR and C_ERROR.

[0220] Preferably, the correction unit CONTROLLER is a fuzzy logic unit with two antecedents, that are the normalized values EN and CEN of the signals ERROR and C_ERROR, respectively, and one consequent, that is the normalized value OUTPUT1 of the correction signal Δ_DI, and preferably it is defined by three membership functions for each antecedent and consequent.

[0221] FIG. 27 shows sample membership functions E_N, E_Z and E_P when the antecedent EN is negative, null (in a fuzzy sense) or is positive, respectively. Similarly, FIGS. 28 and 29 show sample membership functions for the antecedent CEN and the consequent OUTPUT1, respectively.

[0222] The fuzzy correction unit CONTROLLER generates the correction signal Δ_DI according to the nine fuzzy rules shown in FIG. 30. These fuzzy rules are graphically illustrated in FIG. 31 and in the tree-dimensional graph of FIG. 32.

[0223] The output stage of the correction subsystem B2 is preferably composed of an amplifier N₃ of the correction signal Δ_DI, and of a positive feedback loop that generates the feedback signal DI_FB by adding to the amplified correction signal Δ_DI a delayed replica thereof Δ_DI(T−ΔT₂) of a certain time delay ΔT₂:
ID_FB=N₃·Δ_—DI(T)+Δ_—DI(T−ΔT)
Learning Machine

[0224] The learning machine C includes an identification subsystem C2 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 C2 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.

[0225] It 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.

[0226] In practice, the identification subsystem C2 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.

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

[0228] It 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 FIG. 33.

[0229] The core C1 of the lambda sensor of this invention is input with the parameters chosen by the identification subsystem C2. The selected characteristic parameters of the pressure signal are the pressure Pratio40 in a cylinder for a 40° rotation of the crank, the pressure Pratio50 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.

[0230] Indeed, 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.

[0231] However, as schematically illustrated in FIG. 34, the most performing cores are composed of a neural network or a hybrid learning machine.

Virtual Lambda Sensor Core Based on a Neural Network

[0232] FIG. 35 shows a block scheme of a lambda sensor core based on a neural network.

[0233] As already remarked, the preferred inputs of the core are the parameters Pratio40, Pratio50 and Pmax, that are the characteristic parameters of the pressure signal that more often have been chosen by the identification subsystem C2 during the extensive tests that were carried out.

[0234] Signals representative of these parameters are processed by an input stage EDGE DETECTOR, for sampling the input signals and generating a synchronization pulse for the subsequent blocks in cascade.

[0235] According 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.

[0236] The block PRE-PROC is a pre-processor that generates a moving average of the data output by the input stage EDGE DETECTOR, for filtering these data from noise.

[0237] The 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.

[0238] The block POST-PROC is a post-processor analogous to the pre-processor PRE-PROC.

[0239] The 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 FIGS. 36a and 36b of the output signal with and without the blocks PRE-PROC and POST-PROC, respectively.

[0240] The learning part NEURAL NETWORK 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.

[0241] The 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).

Virtual Lambda Sensor Core Based on a Hybrid Learning Machine

[0242] An embodiment of the core of the virtual lambda sensor that includes a so-called “committee machine” is shown in FIG. 37. It includes an input stage EDGE DETECTOR, a pre-processor PRE-PROC and three post-processors POST-PROC that are identical to the homologous blocks of FIG. 35.

[0243] Differently from the embodiment of FIG. 35, the embodiment of FIG. 37 comprises three fuzzy subsystems LEAN FIS, ON-OFF FIS and RICH 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.

[0244] In practice, the fuzzy subsystem ON-OFF FIS works like an on/off lambda sensor. An output multiplexer SWITCH generate an output signal OUTPUT2 that is the signal generated by the subsystem LEAN FIS or RICH FIS depending on the value of the output of the ON-OFF FIS subsystem. Even this last signal is made available (OUTPUT1), thus the virtual lambda sensor core of FIG. 37 may function either as a linear lambda sensor (OUTPUT2) or as an on/off lambda sensor (OUTPUT1)

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

[0246] The 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.

[0247] The fuzzy subsystem RICH FIS is trained such that the output signal preferably ranges between 0.85 and 0.98, while the fuzzy subsystem LEAN FIS outputs a signal that ranges between 1.02 and 1.15. The fuzzy subsystem ON-OFF FIS is trained for estimating very accurately the air/fuel ratio for 0.98<λ<1.02. An eventual error in this range would make the output multiplexer SWITCH select the output of the subsystem LEAN FIS when the air/fuel mixture is reach or vice versa.

[0248] Anyway, such an error does not worsen significantly the accuracy of the device of this invention because the subsystems LEAN FIS and RICH FIS have almost the same accuracy in the range 0.98<λ<1.02. Obviously, it is important that the subsystem ON-OFF FIS do not make an error larger than 0.2 in determining λ, otherwise the accuracy could be negatively affected.

[0249] The 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 ON-OFF FIS subsystem must be filtered because they could induce spurious switching of the multiplexer.

[0250] Compared 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.

[0251] By 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.

[0252] Formula 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.

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

[0254] In 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:

[0255] 1. 4600 rpm, torque=1.5 Nm

[0256] 2. 5600 rpm, torque=4.4 Nm

[0257] 3. 4600 rpm, WOT(Wide Open Throttle) condition

[0258] Yamaha 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.

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

[0260] The 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.

[0261] A 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] FIG. 38 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%.

[0263] During 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 FIGS. 39 and 40.

[0264] Finally, 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] FIGS. 41, 42 and 43 illustrate the results of these tests under transient conditions.

[0266] FIG. 44 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] FIGS. 45 and 46 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.

[0268] The meaning of the labels in FIGS. 44 to 46 is resumed in the following table:

[0269]

	Coppia (Nm)	Torque
	P (kW)	Power
	Bst_Map (bar)	Intake manifold pressure
	ALPHA (%)	Throttle position
	DEG_DGMT (kg/h)	Intake manifold air flow
	T_AIR (° C.)	Air temperature
	FB_TEMP (° C.)	Fuel balance temperature
	T_ACQ_US (° C.)	Cooling system water
		temperature
	T_ASP1 (° C.)	Intake manifold
		temperature
	T_ASP2 (° C.)	Intake manifold
		temperature
	T_SCARI1 (° C.)	Exhaust gases temperature
	T_SCARI2 (° C.)	Exhaust gases temperature
	LAMBDA 1	Lambda value
	FB_VAL (kg/h)	Fuel balance value

• [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 & Exposition, vol. 1357, 208, pages 39 {51, 1998.
• [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.
• [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.
• [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 & Exposition, vol. 1356, pages 39-44, 1998.
• [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 & Exposition, vol. 1356, pages 79-85, 1998.
• [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.
• [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.
• [0277] [8] R. J. Howlett and S. D. Walters. A multi-computer neural network architecture. IEEE Electronics Letters, 35.
• [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.
• [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.
• [0280] [11] A. N. Tikhonov. On solving incorrectly posed problems and method of regularization. Doklady Akademii Nauk USSR, 151.
• [0281] [12] A. N. Tikhonov. On regularization of ill-posed problems. Doklady Akademii Nauk USSR, 153.
• [0282] [13] A. N. Tikhonov and V.Y. Arsenin. Solution of Ill-posed Problems. W. H. Winston.
• [0283] [14] M. J. D. Powell. Radial basis function for multivariate interpolation: A review. IMA Conference on Algorithms for the Approximations of Functions and Data.
• [0284] [15] M. J. D. Powell. The theory of radial basis function approximation in 1990. Advances in Numerical Analysis Vol. II: Wavelets, Subdivision Algorithms, and Radial Basis Functions.
• [0285] [16] M. J. D. Powell. Radial basis function approximations to polynomials. In Numerical Analysis 1987 Proceedings, Dundee, UK, pages 223-241, 1988.
• [0286] [17] V. N. Vapnik. An overview of statistical learning theory. IEEE Trans. On Neural Networks, 9, 1999.
• [0287] [18] S. Haykin. Neural Networks: A comprehensive Foundation. Prentice Hall International Editions, Piscataway, 1999.
• [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. Numerische Mathematik, 31.
• [0289] [20] P. Yoon, S. Park, M. Sunwoo, I. Ohm and K. J. Yoon, Closed Loop Control of Spark Advance and Air-Fuel Ratio in SI Engines Using Cylinder Pressure, SAE 2000-01-0933
• [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,
• [0291] [22] N. R. Pal, K. Pal and J. C. Bezdek, A Mixed c-Means Clustering Model, Proceedings of IEEE International Conference on Fuzzy Systems, pp. 11-21, 1997.
• [0292] [23] R. Maitra, Clustering Massive Datasets with Applications in Software Metrics and Tomography, J. Amer. Stat. Assoc., 1998.
• [0293] [24] H. Frigui and R. Krishnapuram, A Possibilistic approach to Clustering, IEEE Transactions on Fuzzy Systems, vol. 1, No. 2, pp. 98-110, 1993.
• [0294] [25] J. M. Leski, An-Insensitive Fuzzy C-Means Clustering, Int. J. Appl. Math. Comp. Sci., vol. 11, No. 4, pp. 993-1007, 2001.
• [0295] [26] H. Frigui and R. Krishnapuram, A Robust Competitive Clustering Algorithm With Applications in Computer Vision, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 21, No. 5, pp. 450-465, 1999.