- Email: [email protected]
Identification and adaptation of linear look-up table parameters using an efficient recursive least-squares technique
ISA Transactions 48 (2009) 476–483
Contents lists available at ScienceDirect
ISA Transactions journal homepage: www.elsevier.com/locate/isatrans
Identification and adaptation of linear look-up table parameters using an efficient recursive least-squares technique James C. Peyton Jones ∗ , Kenneth R. Muske 1 Center for Nonlinear Dynamics and Control, College of Engineering, Villanova University, 800 Lancaster Ave., Villanova, PA 19085, United States
article
info
Article history: Received 29 January 2009 Received in revised form 9 April 2009 Accepted 13 April 2009 Available online 21 May 2009 Keywords: Lookup table Piece-wise linear models Parameter estimation Least-squares estimation
abstract Linear look-up tables are widely used to approximate and characterize complex nonlinear functional relationships between system input and output. However, both the initial calibration and subsequent realtime adaptation of these tables can be time consuming and prone to error as a result of the large number of table parameters typically required to map the system and the uncertainties and noise in the experimental data on which the calibration is based. In this paper, a new method is presented for identifying or adapting the look-up table parameters using a recursive least-squares based approach. The new method differs from standard recursive least squares algorithms because it exploits the structure of the look-up table equations in order to perform the identification process in a way that is highly computationally and memory efficient. The technique can therefore be implemented within the constraints of typical embedded applications. In the present study, the technique is applied to the identification of the volumetric efficiency look-up table commonly used in gasoline engine fueling strategies. The technique is demonstrated on a Ford 2.0L I4 Duratec engine using time-delayed feedback from a sensor in the exhaust manifold in order to adapt the table parameters online. © 2009 ISA. Published by Elsevier Ltd. All rights reserved.
1. Introduction Model-based representations are concise and able to characterize a wide range of operation, including dynamic behavior, with a small set of parameters. Look-up tables or maps, on the other hand, are useful for characterizing systems where the functional relationship is not known or simply too complex to represent analytically. Linear look-up tables are widely used in function approximation applications ranging from machine learning, floating point numerical computation, nonlinear network analysis, and nonlinear adaptive control systems. Details concerning the theoretical development, construction, and implementation of these tables can be found in the text by Leenaerts and van Bokhoven [1]. In this paper, the primary interest is the approximation of nonlinear system response, as in the early work of Bellman and Roth [2], with particular emphasis on real-time control and monitoring applications. In these applications, the use of look-up tables enables the nonlinear response of the system, within some operating region, to be characterized almost arbitrarily by adjusting the relevant table parameter values. Generally only static relationships are represented in this way, although it is possible to represent system dynamics by
∗
Corresponding author. Tel.: +1 610 519 4216. E-mail addresses: [email protected] (J.C. Peyton Jones), [email protected] (K.R. Muske). 1 Tel.: +1 610 519 6195.
mapping the rate of change of a given signal of interest onto one of the axes of the table. Linear look-up tables are therefore highly flexible, but this flexibility is achieved at the expense of a large number of parameters required to define the grid and the commensurately high effort required to identify their values. For the characterization of multivariable functions, each additional independent variable adds a new axis or dimension to the table and the number of table parameters increases geometrically. For example, a 2-dimensional table with five points on each axis contains only 25 elements but the number of elements increases to 3125 elements for a 5dimensional table. Because of this geometric increase, linear lookup tables and maps are generally confined to 1 or 2-dimensional problems. Even in these cases, however, the issue of table parameter value identification remains. Not only do appropriate values have to be determined in order to populate the table initially, but in many cases there is also a need to adapt or calibrate the table parameters on-line in order to account for slowly time-varying behaviors of the system. In this paper, the look-up table input-output relationship is first cast in the form of an over-determined set of linear equations. A new recursive least squares estimation algorithm is then presented in order to identify the table parameters while minimizing the computational and data storage requirements relative to standard recursive least squares methods. These requirements can be particularly critical in real-time embedded applications where the computational resources are limited. The new method takes
0019-0578/$ – see front matter © 2009 ISA. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.isatra.2009.04.007
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
advantage of the particular sparse structure of the recursive least squares matrix equations for 1-dimensional and 2-dimensional look-up tables to obtain a very efficient numeric algorithm for computing the recursive least squares estimate. In order to apply the techniques presented in this work, it is assumed that the lookup table structure (dimension and grid point placement) has been previously established. The determination of the interpolation grid point locations can be carried out using techniques discussed in [1,3] and is not considered within the scope of this work. This paper is structured as follows. Section 2 presents the mathematical relationships for the linear look-up table parameter identification problem in one and two dimensions respectively. Least-squares identification is discussed in Section 3 along with an analysis of the computational and storage requirements necessary for these algorithms. An alternate recursive least squares approach that takes advantage of the sparse structure of the matrix equations representing the linear look-up table parameter identification problem is presented in Section 4 where a simple 1-dimensional look-up table example is used to illustrate the application of the proposed technique. On-line calibration of the 2-dimensional volumetric efficiency look-up table for an experimental Ford 2.0L I4 Duratec engine is then presented in Section 5 where timedelayed feedback from a wide-ranging exhaust gas oxygen sensor in the exhaust manifold is used to update the table parameters.
477
Fig. 1. 1-dimensional linear interpolation.
Eq. (5) is in a form whose solution is readily obtained, but it is unlikely in practice that the data will be constrained to a single table cell. A general analysis to account for data values falling within several cells is presented in the following section. 2.1. General structure of 1-dimensional tables To generalize the analysis as the input varies over multiple cells, the cell index i is redefined to be function of the time step k i(k) = max(j) | (xj ≤ x(k)) ⇒ xi(k) ≤ x(k) < xi(k)+1 .
(7)
This re-definition propagates to the offset equation 2. Linear look-up table mathematical description Consider the 1-dimensional look-up table illustrated in Fig. 1. The input space of the table is defined by the set of ordered (but not necessarily equi-spaced) points xi , i ∈ [1 . . . nx ]. The output values at each of these points are defined by the corresponding table parameters, θi . Intermediate output values are determined using linear interpolation. The index i of the desired table cell is identified such that the input x falls within the ith table cell i = max(j) | (xj ≤ x) ⇒ xi ≤ x < xi+1
(1)
and the normalized offset of x within the cell is then determined according to, r =
x − xi x i +1 − x i
.
(2)
The interpolated output is given by the relation, y = (1 − r )θi + r θi+1
(3)
which can also be written as the vector–vector dot product shown in Eq. (4). y = (1 − r )
θi r . θi+1
(4)
y = 8θ
(5)
where y is the vector of output values, 8 is a matrix of inputdependent interpolation factors, and θ is a 2-element vector of the table parameters associated with the cell of interest as shown in Eq. (6). y(1) y(2)
y= ..
. y(n) θi θ= . θi+1
(1 − r (1)) (1 − r (2)) 8= .. .
(1 − r (n))
x(k) − xi(k)
r (1) r (2)
.. . r (n)
(8)
xi(k)+1 − xi(k)
and to the interpolation equation. y(k) = (1 − r (k))
r (k)
θi(k) . θi(k)+1
(9)
Although Eq. (9) is similar in form to Eq. (4), it is no longer possible to re-write it according to the form presented in Eqs. (5) and (6) because the two parameters θi(k) , θi(k)+1 no longer refer exclusively to a single cell and may have different identities at different time steps. The vector θ is therefore re-defined to contain all nx table parameter values even though only two of these will be ‘active’ at any given time step.
θ = θ1
θ2
···
θnx
T
.
(10)
To select the appropriate pair of parameters, the matrix 8 (which now has nx columns) is defined with only two non-zero entries in any row ϕ T (k). These entries define the required interpolation factors at the given time step k as follows.
( (1 − r (k)) ϕ (k)j = r (k) T
0
If only a single observation of the output is available, then it is not possible to identify the two unknowns (θi and θi+1 ) using Eq. (4). If, however, there are many observations y(k) of the response to excitation within a single cell i, then Eq. (4) can be written in matrix-vector form
r (k) =
j = i(k) j = i(k) + 1 otherwise.
(11)
For example, if the input x(k) is in cells i(k) = {2, 1, 3} for the first three time-steps and the table itself contained nx = 5 unknown parameters, then the matrix 8 becomes 0 (1 − r (2)) 8= 0
.. .
(1 − r (1)) r (2) 0
.. .
r (1) 0 (1 − r (3))
.. .
0 0 r (3)
.. .
0 0 0 .
(12)
.. .
With these new definitions, Eq. (5) can be used to describe the lookup table response to an input whose values x(k) can vary across the entire input space of the table. 2.2. General structure of 2-dimensional tables
(6) The procedure used to compute the output of a 1-dimensional lookup table in the previous section can be extended to compute the output of a 2-dimensional table. In this case, the input space
478
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
(1 − r2 (1))(1 − r1 (1)) (1 − r2 (2))(1 − r1 (2)) 8= .. .
θ = 2[i1 ,i2 ]
2[i1 +1,i2 ]
(1 − r2 (1))r1 (1) (1 − r2 (2))r1 (2) .. .
2[i1 ,i2 +1]
2[i1 +1,i2 +1]
T
r2 (1)(1 − r1 (1)) r2 (2)(1 − r1 (2))
.. .
r2 (1)r1 (1) r2 (2)r1 (2)
.. .
.
(18) Box I.
stacking each column of 2 into a nx1 · nx2 × 1 dimensional vector as follows.
θi1 +(nx1 −1)i2 = 2i1 ,i2 .
(19)
The matrix of regressors 8 also must be expanded into a sparse matrix with nx1 · nx2 columns but only four non-zero column entries in each row defined by the analog of Eq. (11) for two dimensions.
(1 − r2 ) · (1 − r1 ) (1 − r2 ) · r1 ϕ T (k)j = r2 · (1 − r1 ) r2 · r1 0
Fig. 2. 2-dimensional linear interpolation.
of the table is defined by a 2-dimensional grid of points x1[i1 ] , i1 ∈ [1 . . . nx1 ] and x2[i2 ] , i2 ∈ [1 . . . nx2 ]. The matrix of table parameters 2i1 ,i2 contains the output values at each point on the grid and defines the mesh surface of the table. Intermediate values on this surface are obtained using 2-dimensional linear interpolation. The location of the desired table cell is first identified such that the inputs x1 (k), x2 (k) at time-step k fall within the cell whose ‘lower left’ corner has index i1 , i2 . i1 (k) = max(j) | (x1 [j] ≤ x1 (k)) i2 (k) = max(j) | (x2 [j] ≤ x2 (k)).
r2 (k) =
x1 (k) − x1[i1 (k)] x1[i1 (k)+1] − x1[i1 (k)] x2 (k) − x2[i2 (k)] x2[i2 (k)+1] − x2[i2 (k)]
(15)
ˆ t) = V (θ,
(16)
t 1X
2 k =1
(17)
with 8 and θ defined as shown in Eq. (18) (Box I). When the input crosses cell boundaries, however, the parameter vector θ must be expanded to incorporate every value in the table. Since these values are conceptually arranged as the 2-dimensional matrix 2, the parameter vector is obtained by
2 β t −k y(k) − ϕ T (k)θˆ
(21)
where the ‘forgetting factor’ β ≤ 1 diminishes the weight of older data relative to data that was obtained more recently. Typically β is set to a value between 0.95 and 1 depending on the desired rate of adaptation. In off-line applications, β is generally set equal to 1 and the closed form solution is given by
ˆ t) = θ(
t X k=1
that relates the output to the four values at each corner of the cell. Note that the time-step dependency of r1 , r2 , i1 and i2 is not shown explicitly in Eq. (16) for reasons of notational clarity. As long as the cell indices i1 (k) and i2 (k) do not change value between time-steps, then the vector of table outputs y can be obtained according to y = 8θ
3. Recursive least-squares parameter identification
(14)
(13)
and the table output is obtained by interpolating first along the two parallel cell edges in the x1 direction and then between these two points in the x2 direction as shown in Fig. 2. The structure of this interpolation can be seen in the matrix equation
(1 − r2 ) · (1 − r1 ) T 2[i1 ,i2 ] (1 − r2 ) · r1 2[i1 +1,i2 ] y(k) = r2 · (1 − r1 ) 2[i1 ,i2 +1] r2 · r1 2[i1 +1,i2 +1]
Note that the condition on j is obtained by mapping the 2 indices in Eq. (19) according to Eq. (20). With these definitions, Eq. (17) can be used to describe the look-up table response to an input x(k) that varies across the entire 2-dimensional input space of the table.
Eqs. (5) and (17) show how output values are obtained from the table parameters θ and inputs x. In this application, however, the table values θi are to be identified given observed values for the output y and input x. The set of over-determined linear equations represented by Eq. (5) for the 1-dimensional table or Eq. (17) for the 2-dimensional table, with the corresponding definitions for the parameter vector θ and the matrix of regressors 8, can be solved to obtain the unknown table parameters θi using least squares [4,5]. This method corresponds to finding the parameters that minimize the weighted mean squared error between observations and model predictions according to the cost function
The normalized offsets of x1 (k), x2 (k) in each cell dimension are given by the pair of equations r1 (k) =
j = i1 (k) + (nx1 − 1)i2 (k) j = i1 (k) + (nx1 − 1)i2 (k) + 1 j = i1 (k) + nx1 i2 (k) (20) j = i1 (k) + nx1 i2 (k) + 1 otherwise.
! −1 ϕ(k)ϕ T (k)
t X
! ϕ(k)y(k)
(22)
k =1
where the solution to Eq. (22) should yield estimated values for the table parameters of interest as long as the input is sufficiently exciting. The estimates will be unbiased providing the regressors are uncorrelated with the noise process and providing the lookup table model structure admits an exact description of the true system. It is worth noting that in most applications the piecewise linear look-up table model is only an approximation of the (generally smooth) nonlinear system that it is used to represent. The table parameters are therefore likely to be ‘biased’ relative to the exact values of the system at the table grid points, but the resultant ‘bias’ will tend to minimize the overall mean square prediction error between the true process output and that predicted from the linear look-up table approximation—at least
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
ˆ 0) = θinit , θ(
479
P (0) = δ I
(28)
where δ 1 is used indicate uncertainty in the initial parameter estimates if reasonable starting values are not known a priori. 4. An efficient recursive least squares algorithm
Fig. 3. Operating point bias in RLS estimation.
for the data set used during identification. Such behavior is often acceptable or even desirable but it is important to ensure that the data set used during identification is well distributed across and within the cells since the model fit will favor any dominant regions of operation. Fig. 3, for example, illustrates a 1-dimensional slice of a function with high curvature that is being approximated by a linear interpolation table. If the operating condition remains in a narrow range about point (a), then the estimated parameters will be as shown in the left hand plot of the figure. If instead the operating point were constrained to a narrow region around point (b), then the fit shown in the right hand plot of Fig. 3 would result. In order to avoid these issues, the cell grid should be designed to limit the maximum error between the ‘true’ function and its linear approximation and the RLS algorithm should also be fed with data distributed across the cell. If necessary, this distribution can be achieved by windowing the data so that the algorithm only executes when the operating condition has moved significantly from its previous position. In on-line applications, it is usually not possible or desirable to store large data records. It is much more efficient merely to update the table parameter estimates as the next sample of data becomes available. It is also useful to retain the forgetting factor so that the parameters are weighted towards the most recent data. The recursive form, which relates the estimates at time t to those at time t − 1, is obtained by extracting the uppermost elements of the summations in Eq. (22) and re-arranging to give
h
i ˆ t ) = θ( ˆ t − 1) + R (t )ϕ(t ) y(t ) − ϕ (t )θ( ˆ t − 1) θ( −1
T
(23)
where R(t ) =
t X
β t −k ϕ(k)ϕ T (k) = β R(t − 1) + ϕ(t )ϕ T (t ).
(24)
k=1
In order to avoid having to compute the inverse of matrix R(t ) in Eq. (23), it is convenient to define the covariance matrix P (t ) = R−1 (t ) and update this inverse directly. Applying the matrix inversion lemma to Eq. (24) enables P (t ) to be expressed as P (t ) =
1
β
I − L(t )ϕ T (t ) P (t − 1)
(25)
R(t )1θˆ (t ) = f (t )
(29)
where
1θˆ (t ) = θˆ (t ) − θˆ (t − 1) f (t ) = ϕ(t ) y(t ) − ϕ T (t )θˆ (t − 1) .
(30)
As noted above, the general disadvantage of this approach is that the linear system in Eq. (29) involving the R(t ) matrix must be solved at each time step. In this particular application, however, the special structure of the regression vectors (Eq. (11) for 1-D or Eq. (20) for 2-D tables) results in a sparse banded R(t ) matrix that has minimal memory requirements and for which a solution to Eq. (29) can be efficiently obtained as discussed in the sequel. These benefits do not hold for the standard RLS algorithm because the P (t ) matrix does not retain a sparse banded structure upon inversion of the R(t ) matrix in general. 4.1. Recursive least squares approach for 1-dimensional tables
where L(t ) =
One of the drawbacks of the recursive least-squares formulation is the requirement to store and update a P (t ) matrix of dimension (nx1 )2 for a 1-dimensional table, or of dimension (nx1 · nx2 )2 in the 2-dimensional case. This storage requirement therefore increases rapidly with table dimension and size. A 1-dimensional table with 8 data points requires only 64 storage locations for the P (t ) matrix, but a 2-dimensional table with 8 data points on each axis requires 3996. Memory storage requirements can therefore become a constraining issue in real-time embedded applications where computational resources are limited. One way to address this issue is to approximate the covariance matrix by its main and first upper and lower diagonals on the basis that only neighboring cells are strongly correlated to each other [6]. This approximation significantly reduces the data storage requirements to 2n elements as compared to the n2 elements required for the full matrix, where n represents the total number of parameters in the table, but the validity of the approximation has not been proven theoretically. In this paper a different approach is adopted which does not rely upon approximations, but instead exploits the sparse structure of the R(t ) matrix in order to achieve a memory-efficient and computationally-efficient implementation. The advantage of this approach is that the recursive least squares estimate is itself unchanged so all the convergence and stability properties carry over from the standard algorithm. It is only the numeric algorithm with which the estimate is computed that has been redesigned. The new method is based on computing the estimate by rearranging Eq. (23) and using the matrix R(t ), defined in Eq. (24), directly rather than by means of the matrix inversion lemma in Eqs. (25) and (26). Re-arranging Eq. (23) results in the following linear system of equations.
P (t − 1)ϕ(t ) . β + ϕ T (t )P (t − 1)ϕ(t )
(26)
Furthermore, using Eqs. (25) and (26), it can be shown that P (t )ϕ(t ) = L(t ). Eq. (23) may therefore be re-written as
h i ˆ t ) = θ( ˆ t − 1) + L(t ) y(t ) − ϕ T (t )θ( ˆ t − 1) . θ(
(27)
Eqs. (25)–(27) define the well known Recursive Least-Squares (RLS) algorithm. A more detailed derivation, together with a discussion of numerical and convergence issues can be found in [4,5]. The algorithm is initialized by setting
Consider the evaluation of Eq. (24) at time-step t, where ϕ (t) is defined by Eq. (11), in the case of a 1-dimensional look-up table. Inspection of Eq. (11) reveals that the update term for R(t ) in Eq. (24) is a sparse matrix whose only non-zero terms are the 2 × 2 symmetric square matrix shown in Eq. (31).
(1 − r (k))2 ϕ(k)j ϕ (k)j = r (k)(1 − r (k)) T
r (k)(1 − r (k)) . r 2 (k)
(31)
These non-zero elements are positioned starting at location j, j within the R(t ) matrix as shown in Fig. 4. The 4-element submatrix defined in Eq. (31) is positioned on the main diagonal of R(t )
480
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
(1 − r2 )2 (1 − r1 )2 (1 − r2 )2 (1 − r1 )r . 0..
ϕ(k)j ϕ (k)j = r2 (1 − r2 )(1 − r1 )2 r2 r1 (1 − r2 )(1 − r1 ) T
(1 − r2 )2 (1 − r1 )r1 (1 − r2 )2 r12 . 0..
r2 r1 (1 − r2 )(1 − r1 ) r2 r12 (1 − r2 )
...0... ...0... ...0... ...0...
r2 (1 − r2 )(1 − r1 )2 r2 r1 (1 − r2 )(1 − r1 )
.
0..
( − r1 ) ( − r1 )
r22 1 r22 r1 1
2
r2 r1 (1 − r2 )(1 − r1 ) r2 r12 (1 − r2 )
.
0.. 2 r2 r1 (1 − r1 ) r22 r12
.
(34)
Box II.
Fig. 4. Structure of the R(t ) matrix for 1-dimensional tables.
at the location specified by the index j = i(k) that represents the currently active cell. The matrix R(t ) is therefore symmetric and tri-diagonal as also illustrated in Fig. 4. The advantage of this matrix structure is two-fold. First, the R(t ) matrix requires only storage for the (2 ∗ n − 1) non-zero elements, as compared to n2 non-zero elements for the corresponding inverse matrix P (t ), and second, the solution to the linear system in Eq. (29) is efficiently obtained using the Thomas algorithm [7]. The Thomas algorithm is briefly summarized as follows. Let b(i), c (i) represent the elements on main and first upper and lower diagonals of the matrix R(t ), respectively, and let f (i) represent the elements of the vector f (t ). The forward elimination phase of the algorithm is defined as follows: b1(1) = 1; for k = 2 to n, b1(k) = b(k) − c (k) ∗ c (k − 1)/b1(k − 1); f 1(k) = f (k) − c (k) ∗ f 1(k − 1)/b(k − 1); end;
(32)
The backward substitution phase then yields the solution according to:
1θ(n) = f 1(n)/b1(n), for k = n − 1 downto 1, 1θ(k) = f 1(k) − c (k) ∗ f 1(k + 1)/b1(k); end;
Fig. 5. Structure of the R(t ) matrix for 2-dimensional tables.
necessary to accommodate differences in the form of the regression matrices. In the 2-dimensional case, inspection of Eq. (20) shows that the update component of the R(t ) matrix in Eq. (24) has the form as in Box II: In this case, the non-zero elements are grouped as four 2 × 2 symmetric sub-matrices separated from each other by an index offset of nx1 in either dimension. This group of four sub-matrices is positioned starting at element j, j of the R(t ) matrix where j = i1 (k) + (nx1 − 1)i2 (k) corresponds to the ‘lower left’ corner of the active cell as shown in Fig. 5. As the active cell changes, the position within the R(t ) matrix shifts up or down the leading diagonal resulting in the symmetric, banded block tri-diagonal matrix that is illustrated in Fig. 5. The banded block structure of the R(t ) matrix for 2-dimensional tables prevents the direct application of the Thomas algorithm as presented in Section 4.1. To address this issue, the R(t ) matrix is partitioned into a tri-diagonal portion RD (t ), similar to that shown in Fig. 4, and a remainder (R(t ) − RD (t )) that contains the upper and lower banded part. Eq. (29) can then be rewritten as RD (t )1θˆ (t ) = f (t ) − R(t ) − RD (t ) 1θˆ (t ).
(33)
The process for efficient estimation of 1-dimensional table parameters can therefore be summarized as follows:
ˆ 0) = θinit , R(0) = δ I where (i) Initialize algorithm by setting θ( δ < 1 is used indicate uncertainty in the initial parameter estimates if reasonable starting values are not known a priori. (ii) Update the sparse, tri-diagonal matrix R(t ) according to Eqs. (24) and (31). (iii) Compute 1θˆ (t ) and hence update the parameter estimate using Eqs. (29)–(33). (iv) Loop back to step (ii).
(35)
The solution to this equation can now be determined by applying the Thomas algorithm iteratively where 1θˆ (t ) in the right-hand side of the equation is approximated by the estimate obtained at the previous iteration, RD (t )1θˆm (t ) = f (t ) − R(t ) − RD (t ) 1θˆm−1 (t )
(36)
and where m represents the iteration number. The iteration is initiated by setting 1θˆ0 (t ) = 0. Practical application of the method indicates that the solution converges to within 2% of the ‘true’ values obtained by direct solution with the full R(t ) matrix within a few iterations. Apart from this additional iteration step, the method is essentially identical to that presented for the 1dimensional table case.
4.2. Recursive least squares approach for 2-dimensional tables
4.3. 1-dimensional look-up table illustrative example
A similar procedure can be used to estimate the parameters of a 2-dimensional lookup table, although some modifications are
To assess the performance of the various recursive leastsquares algorithms and approximations when applied to table
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
k
25
Known Table
y(k)
20 Estimated parameters, θ
x(k)
481
RLS Estimator k Fig. 6. Simulation block diagram.
15
10
5 25 Output, y(k)
0 2.5
20
3
3.5 4 4.5 Table x-axis values [-]
5
5.5
Fig. 8. Snapshots of estimated table parameters (circles); final estimate (diamond); ‘true’ values (solid).
15
5
0
0
0.2
0.4
0.6
0.8
1 1.2 Time [s]
1.4
1.6
1.8
2
Fig. 7. Measured table input and output data.
identification, simulation data generated from a known table with zero mean unity variance additive noise is used for parameter estimation as shown in Fig. 6. For the sake of illustration, only a small table with seven entries is used where the ‘true’ table values θ0 at nodes x are defined by x = 2.5
3.0
3.5
θ0 = 10
15
19
4.0 21
4.5 22
5.0 21.5
5.5
T
T
.
20
Estimated parameters, θ
10
25 20 15 10 5
3
0.1 3.5
0.05
4
4.5 Table x-axis values [-]
5
5.5
Time [s]
0
Fig. 9. Evolution of estimated parameters with time.
25
(37)
The input is a uniformly distributed random signal designed to span the full input range of the table and the time step is 0.01 s. The resultant input/output time histories, shown in Fig. 7, were then presented to the RLS identification algorithm. In the simulation study described here, θinit is set equal to zero and the R(0) matrix is started with δ = 0.1 on initialization. The value of β is set initially to 0.98. Various snapshots of the estimated table parameters are shown (as circles joined by dotted lines) in Fig. 8. The last of these, corresponding to the parameter values at t = 1 s, is indicated with a diamond symbol. Also shown in Fig. 8 (as a solid line) is the ‘true’ table given by Eq. (37). As shown in the figure, the two match quite closely indicating successful convergence of the algorithm. The time evolution of the estimated parameters is illustrated in Fig. 9. The form of the table can be seen emerging with time from the initial zero parameter values. In the early stages of the table identification, some points adapt more rapidly than others, presumably because the input was exciting these cells more than the others. A more quantitative assessment of the convergence process can be made by slicing Fig. 9 along lines parallel with the time axis that result in the traces shown in Fig. 10. Each trace represents the evolution of an individual parameter in the table, although only four such traces are shown for reasons of clarity. After an initial sharp transient the parameters gradually converge towards their ‘true’ values. The effects of noise are seen in occasional fluctuations away from this ideal.
0.2 0.15
0 2.5
20 Estimated parameters, θi
Measured table input and output [-]
Input, x(k)
15
10 θ7 ...=20 θ5 ...=22
5
θ2 ...=15 θ1 ...=10
0
0
0.2
0.4
0.6
0.8
1 1.2 Time [s]
1.4
1.6
1.8
2
Fig. 10. Evolution of estimated parameters (β = 0.98).
The tradeoff between noise immunity and convergence rate is controlled by the value of the forgetting factor β . The results shown in Figs. 8–10 were obtained with β equal to 0.98. Fig. 11 presents the identification results when the simulation is repeated with β reduced to 0.94. Although the same general trends are observed, the estimates are now more sensitive to noise but should adapt more quickly to genuine changes in the ‘true’ parameter values. The value of β = 0.98 seems more appropriate in this case, though
482
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
25
Estimated parameters, θi
φ d(k) 20
φ m(k)
Engine
rpm(k) load(k)
15 ηvol Table / Map
10 θ7 ...=20 θ5 ...=22
5
z-d
z-d
k
θ2 ...=15 θ1 ...=10
0
0
0.2
0.4
0.6
0.8 1 1.2 Time [s]
1.4
1.6
1.8
z-(d+1)
RLS Estimator
ηvol
2 Fig. 12. 2-D identification block diagram.
Fig. 11. Evolution of estimated parameters (β = 0.94).
care should be taken to tune the forgetting factor in different applications. 5. Identification of a 2-dimensional volumetric efficiency map Because of their computational efficiency and the highly nonlinear nature of automotive powertrain systems, look-up tables are used extensively in automotive engine control and monitoring applications. In order to identify the parameters, it is often necessary to run the engine in steady-state conditions at each point in the table grid, a process which is expensive and time-consuming. The parameter identification itself is also often heuristic, relying on operator judgment for recording the appropriate value from a noisy signal. Even when the recorded value is accurate, errors can also arise if the engine is not operated at exactly the point defined by the table grid. Automated table identification is therefore desirable, both off-line using engine testbed data, and online in order to adapt the table in-vehicle as the engine ages. In either case, both the table input and desired output need to be available during the identification process. In the offline case, such measurements are often available from the engine test-bed instrumentation. In the on-line case, however, the required output is generally only known ‘after the event’. If it were not, a table would hardly be required. Knock spark angle correction, applied by table lookup, is one such example. When knock is detected in a previous firing, the tabulated value can be adjusted to avoid knock in subsequent firings [8–10]. In many cases, the table updates are performed heuristically sometimes simply overwriting previous table values and thereby discarding the information contained in earlier measurements [9,10]. The calibration of engine fueling is another example where a table is used to predict the required fueling for the current cycle as part of the air/fuel ratio feedforward control strategy [11]. After the firing has taken place, and the gas has reached the exhaust gas oxygen sensor, the measured air/fuel ratio can be used to update the table for subsequent cycles. The identification procedure presented in this work is applied to this example. The algorithm is implemented on a rapid prototyping system attached to a Ford 2.0L I4 engine and then used to estimate the parameters of the volumetric efficiency map. A block diagram of the configuration is shown in Fig. 12. In this case the ‘true’ map is represented by the behavior of the engine itself. The inputs are speed and load and the output is volumetric efficiency. The output is not measured directly, but is computed assuming the injectors are well calibrated. The measured equivalence ratio φm , is related to desired equivalence ratio φd according to
φm (k) =
mbf (k − d) ma (k − d)
AFRs φd (k − d)
(38)
where mbf is the base fuel mass, ma the air mass flow, and AFRs is the stoichiometric air/fuel ratio by mass. The lag d represents the transport delay of the engine in engine revolutions. We assume φd is constant, and that the other variables change only slowly. The dynamics of the response are therefore ignored in Eq. (38). The base fuel mass is generally computed from the ratio of the estimated air mass to the stoichiometric air/fuel ratio. Air mass flows are also assumed proportional to volumetric efficiency ηv ol such that the following relationships hold.
φm (k) =
ˆ a (k − d + 1) m ma (k − d)
φd (k) =
ηˆ vol (k − d + 1) φd (k) ηvol (k − d)
φd (k) ηˆ vol (k − d + 1). (39) φm (k) Note that the delay on φd can be ignored since its value is held constant at φd = 1 in this example. ⇒ ηvol (k − d) =
The value of the measured volumetric efficiency is computed ‘after the event’ according to Eq. (39) and the appropriately delayed speed and load inputs are then fed to the RLS algorithm implemented on the rapid prototyping system as shown in Fig. 12. The algorithm is executed at each engine firing where a forgetting factor of β = 0.995 was used. For practical reasons, only four table parameters were monitored at any instant, these being the parameters associated with the current cell whose values were being adapted. Fig. 13 shows the results that were obtained as engine speed (Fig. 14) varied slowly across two cells. Initially the parameter values were constant at their default calibration value. At t = 2 s, however, the RLS algorithm was enabled. The process of algorithm initialization results in an undesirable transient, but the parameters then begin to adapt toward their new values. At t = 7 s, engine speed crosses the 1600 rpm level which moves the operating condition into the next cell of the table. The two lower speed parameters 25,4 , 25,5 represented by the dotted traces in Fig. 13 for 0 < t < 7 are therefore dropped from the monitored list, and replaced by the two higher speed parameters 27,4 , 27,5 associated with the new cell (dotted traces t > 7 s). The parameters, 26,4 , 26,5 , are common to both cells, so there is no change in identity for these traces. As the operating condition enters the new cell, the two new parameters 27,4 , 27,5 adapt rapidly from their initial values at t = 7 s to levels around 0.9, 0.95 respectively. However the degree of uncertainty or ‘noise’ in the estimated values, seems higher for these parameters than for the parameters 26,4 , 26,5 which now define the low speed end of the cell. This is believed to be because
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
1.3 Θ5,4 , Θ7,4
1.2 Estimated parameters, θi
6. Conclusions
Θ6,4
1.25
Θ6,5 Θ5,5 , Θ7,5
1.15 1.1 1.05 1 0.95 0.9 0.85 0.8
0
5
10
15
20
483
25
An efficient recursive least-squares identification algorithm in terms of both computation and storage is presented for the initial determination and adaptation of the parameters in linear look-up tables. This algorithm has been applied successfully to 1- and 2dimensional problems in simulation and in real time. The real-time application concerned the estimation of volumetric efficiency for engine fueling, but the method can be applied in any case where input and (generally delayed) output signals are available. For the method to be robust, care must be taken to ensure that the data is well distributed over the range of operation defined by the table and that the table grid is well designed relative to the curvature of the function being approximated. Further work into developing windowing and weighting strategies to address these issues is ongoing.
Time [s]
Acknowledgments
Fig. 13. Evolution of estimated ηv ol parameters.
The authors gratefully acknowledge financial and technical support from Ford Motor Company, Johnson-Matthey, ExxonMobil, and Mr. E. Barry, ChE ‘54.
1700
Engine Speed [rpm]
1650
References
1600
1550
1500
1450
0
5
10
15
20
25
Time [s] Fig. 14. Engine speed time history.
engine speed fell mainly in the lower end of the cell, resulting in better estimates of the nearby parameters (26,4 , 26,5 ), compared to the others which are further away.
[1] Leenaerts D, van Bokhoven W. Piecewise linear modeling and analysis. Berlin: Springer; 1998. [2] Bellman R, Roth R. Curve fitting by segmented straight lines. Amer Statist Assoc J 1969;64:1079–84. [3] Nelles O. Nonlinear system identification. New York: Springer-Verlag; 2000. [4] Ljung I. System identification: Theory for the user. 2nd ed. Upper Saddle River (NJ): PTR Prentice Hall; 1999. [5] Soderstrom T, Stoica P. System identification. UK: Prentice Hall International; 1989. [6] Peyton Jones J, Muske K. Automatic calibration of 1 and 2-D lookup tables using recursive least squares identification techniques. In: 2007 SAE world congress. Paper 2007-1333. 2007. [7] Press W, Teukolsky S, Vetterling W, Flannery B. Numerical recipes. 2nd ed. Cambridge (UK): Cambridge University Press; 1992. [8] Wu G. A table update method for adaptive knock control. In: 2006 SAE world congress. Paper 2006-0607. 2006. [9] Entenmann R, Unland S, Torno O, Haeming W, Rothhaar U, Surjadi I. et al. Method for the adaptive knock control of an internal combustion engine. US Patent No. 5645033. 1997. [10] Entenmann R, Unland S, Haeming W. Process for controlling knocking in internal combustion engine. US Patent No. 5243842. 1993. [11] Knorad JA, Freeland M, Czuhai BS. Method and system for controlling flexible fuel vehicle fueling. US Patent No. 5467755. 1995.
Contents lists available at ScienceDirect
ISA Transactions journal homepage: www.elsevier.com/locate/isatrans
Identification and adaptation of linear look-up table parameters using an efficient recursive least-squares technique James C. Peyton Jones ∗ , Kenneth R. Muske 1 Center for Nonlinear Dynamics and Control, College of Engineering, Villanova University, 800 Lancaster Ave., Villanova, PA 19085, United States
article
info
Article history: Received 29 January 2009 Received in revised form 9 April 2009 Accepted 13 April 2009 Available online 21 May 2009 Keywords: Lookup table Piece-wise linear models Parameter estimation Least-squares estimation
abstract Linear look-up tables are widely used to approximate and characterize complex nonlinear functional relationships between system input and output. However, both the initial calibration and subsequent realtime adaptation of these tables can be time consuming and prone to error as a result of the large number of table parameters typically required to map the system and the uncertainties and noise in the experimental data on which the calibration is based. In this paper, a new method is presented for identifying or adapting the look-up table parameters using a recursive least-squares based approach. The new method differs from standard recursive least squares algorithms because it exploits the structure of the look-up table equations in order to perform the identification process in a way that is highly computationally and memory efficient. The technique can therefore be implemented within the constraints of typical embedded applications. In the present study, the technique is applied to the identification of the volumetric efficiency look-up table commonly used in gasoline engine fueling strategies. The technique is demonstrated on a Ford 2.0L I4 Duratec engine using time-delayed feedback from a sensor in the exhaust manifold in order to adapt the table parameters online. © 2009 ISA. Published by Elsevier Ltd. All rights reserved.
1. Introduction Model-based representations are concise and able to characterize a wide range of operation, including dynamic behavior, with a small set of parameters. Look-up tables or maps, on the other hand, are useful for characterizing systems where the functional relationship is not known or simply too complex to represent analytically. Linear look-up tables are widely used in function approximation applications ranging from machine learning, floating point numerical computation, nonlinear network analysis, and nonlinear adaptive control systems. Details concerning the theoretical development, construction, and implementation of these tables can be found in the text by Leenaerts and van Bokhoven [1]. In this paper, the primary interest is the approximation of nonlinear system response, as in the early work of Bellman and Roth [2], with particular emphasis on real-time control and monitoring applications. In these applications, the use of look-up tables enables the nonlinear response of the system, within some operating region, to be characterized almost arbitrarily by adjusting the relevant table parameter values. Generally only static relationships are represented in this way, although it is possible to represent system dynamics by
∗
Corresponding author. Tel.: +1 610 519 4216. E-mail addresses: [email protected] (J.C. Peyton Jones), [email protected] (K.R. Muske). 1 Tel.: +1 610 519 6195.
mapping the rate of change of a given signal of interest onto one of the axes of the table. Linear look-up tables are therefore highly flexible, but this flexibility is achieved at the expense of a large number of parameters required to define the grid and the commensurately high effort required to identify their values. For the characterization of multivariable functions, each additional independent variable adds a new axis or dimension to the table and the number of table parameters increases geometrically. For example, a 2-dimensional table with five points on each axis contains only 25 elements but the number of elements increases to 3125 elements for a 5dimensional table. Because of this geometric increase, linear lookup tables and maps are generally confined to 1 or 2-dimensional problems. Even in these cases, however, the issue of table parameter value identification remains. Not only do appropriate values have to be determined in order to populate the table initially, but in many cases there is also a need to adapt or calibrate the table parameters on-line in order to account for slowly time-varying behaviors of the system. In this paper, the look-up table input-output relationship is first cast in the form of an over-determined set of linear equations. A new recursive least squares estimation algorithm is then presented in order to identify the table parameters while minimizing the computational and data storage requirements relative to standard recursive least squares methods. These requirements can be particularly critical in real-time embedded applications where the computational resources are limited. The new method takes
0019-0578/$ – see front matter © 2009 ISA. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.isatra.2009.04.007
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
advantage of the particular sparse structure of the recursive least squares matrix equations for 1-dimensional and 2-dimensional look-up tables to obtain a very efficient numeric algorithm for computing the recursive least squares estimate. In order to apply the techniques presented in this work, it is assumed that the lookup table structure (dimension and grid point placement) has been previously established. The determination of the interpolation grid point locations can be carried out using techniques discussed in [1,3] and is not considered within the scope of this work. This paper is structured as follows. Section 2 presents the mathematical relationships for the linear look-up table parameter identification problem in one and two dimensions respectively. Least-squares identification is discussed in Section 3 along with an analysis of the computational and storage requirements necessary for these algorithms. An alternate recursive least squares approach that takes advantage of the sparse structure of the matrix equations representing the linear look-up table parameter identification problem is presented in Section 4 where a simple 1-dimensional look-up table example is used to illustrate the application of the proposed technique. On-line calibration of the 2-dimensional volumetric efficiency look-up table for an experimental Ford 2.0L I4 Duratec engine is then presented in Section 5 where timedelayed feedback from a wide-ranging exhaust gas oxygen sensor in the exhaust manifold is used to update the table parameters.
477
Fig. 1. 1-dimensional linear interpolation.
Eq. (5) is in a form whose solution is readily obtained, but it is unlikely in practice that the data will be constrained to a single table cell. A general analysis to account for data values falling within several cells is presented in the following section. 2.1. General structure of 1-dimensional tables To generalize the analysis as the input varies over multiple cells, the cell index i is redefined to be function of the time step k i(k) = max(j) | (xj ≤ x(k)) ⇒ xi(k) ≤ x(k) < xi(k)+1 .
(7)
This re-definition propagates to the offset equation 2. Linear look-up table mathematical description Consider the 1-dimensional look-up table illustrated in Fig. 1. The input space of the table is defined by the set of ordered (but not necessarily equi-spaced) points xi , i ∈ [1 . . . nx ]. The output values at each of these points are defined by the corresponding table parameters, θi . Intermediate output values are determined using linear interpolation. The index i of the desired table cell is identified such that the input x falls within the ith table cell i = max(j) | (xj ≤ x) ⇒ xi ≤ x < xi+1
(1)
and the normalized offset of x within the cell is then determined according to, r =
x − xi x i +1 − x i
.
(2)
The interpolated output is given by the relation, y = (1 − r )θi + r θi+1
(3)
which can also be written as the vector–vector dot product shown in Eq. (4). y = (1 − r )
θi r . θi+1
(4)
y = 8θ
(5)
where y is the vector of output values, 8 is a matrix of inputdependent interpolation factors, and θ is a 2-element vector of the table parameters associated with the cell of interest as shown in Eq. (6). y(1) y(2)
y= ..
. y(n) θi θ= . θi+1
(1 − r (1)) (1 − r (2)) 8= .. .
(1 − r (n))
x(k) − xi(k)
r (1) r (2)
.. . r (n)
(8)
xi(k)+1 − xi(k)
and to the interpolation equation. y(k) = (1 − r (k))
r (k)
θi(k) . θi(k)+1
(9)
Although Eq. (9) is similar in form to Eq. (4), it is no longer possible to re-write it according to the form presented in Eqs. (5) and (6) because the two parameters θi(k) , θi(k)+1 no longer refer exclusively to a single cell and may have different identities at different time steps. The vector θ is therefore re-defined to contain all nx table parameter values even though only two of these will be ‘active’ at any given time step.
θ = θ1
θ2
···
θnx
T
.
(10)
To select the appropriate pair of parameters, the matrix 8 (which now has nx columns) is defined with only two non-zero entries in any row ϕ T (k). These entries define the required interpolation factors at the given time step k as follows.
( (1 − r (k)) ϕ (k)j = r (k) T
0
If only a single observation of the output is available, then it is not possible to identify the two unknowns (θi and θi+1 ) using Eq. (4). If, however, there are many observations y(k) of the response to excitation within a single cell i, then Eq. (4) can be written in matrix-vector form
r (k) =
j = i(k) j = i(k) + 1 otherwise.
(11)
For example, if the input x(k) is in cells i(k) = {2, 1, 3} for the first three time-steps and the table itself contained nx = 5 unknown parameters, then the matrix 8 becomes 0 (1 − r (2)) 8= 0
.. .
(1 − r (1)) r (2) 0
.. .
r (1) 0 (1 − r (3))
.. .
0 0 r (3)
.. .
0 0 0 .
(12)
.. .
With these new definitions, Eq. (5) can be used to describe the lookup table response to an input whose values x(k) can vary across the entire input space of the table. 2.2. General structure of 2-dimensional tables
(6) The procedure used to compute the output of a 1-dimensional lookup table in the previous section can be extended to compute the output of a 2-dimensional table. In this case, the input space
478
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
(1 − r2 (1))(1 − r1 (1)) (1 − r2 (2))(1 − r1 (2)) 8= .. .
θ = 2[i1 ,i2 ]
2[i1 +1,i2 ]
(1 − r2 (1))r1 (1) (1 − r2 (2))r1 (2) .. .
2[i1 ,i2 +1]
2[i1 +1,i2 +1]
T
r2 (1)(1 − r1 (1)) r2 (2)(1 − r1 (2))
.. .
r2 (1)r1 (1) r2 (2)r1 (2)
.. .
.
(18) Box I.
stacking each column of 2 into a nx1 · nx2 × 1 dimensional vector as follows.
θi1 +(nx1 −1)i2 = 2i1 ,i2 .
(19)
The matrix of regressors 8 also must be expanded into a sparse matrix with nx1 · nx2 columns but only four non-zero column entries in each row defined by the analog of Eq. (11) for two dimensions.
(1 − r2 ) · (1 − r1 ) (1 − r2 ) · r1 ϕ T (k)j = r2 · (1 − r1 ) r2 · r1 0
Fig. 2. 2-dimensional linear interpolation.
of the table is defined by a 2-dimensional grid of points x1[i1 ] , i1 ∈ [1 . . . nx1 ] and x2[i2 ] , i2 ∈ [1 . . . nx2 ]. The matrix of table parameters 2i1 ,i2 contains the output values at each point on the grid and defines the mesh surface of the table. Intermediate values on this surface are obtained using 2-dimensional linear interpolation. The location of the desired table cell is first identified such that the inputs x1 (k), x2 (k) at time-step k fall within the cell whose ‘lower left’ corner has index i1 , i2 . i1 (k) = max(j) | (x1 [j] ≤ x1 (k)) i2 (k) = max(j) | (x2 [j] ≤ x2 (k)).
r2 (k) =
x1 (k) − x1[i1 (k)] x1[i1 (k)+1] − x1[i1 (k)] x2 (k) − x2[i2 (k)] x2[i2 (k)+1] − x2[i2 (k)]
(15)
ˆ t) = V (θ,
(16)
t 1X
2 k =1
(17)
with 8 and θ defined as shown in Eq. (18) (Box I). When the input crosses cell boundaries, however, the parameter vector θ must be expanded to incorporate every value in the table. Since these values are conceptually arranged as the 2-dimensional matrix 2, the parameter vector is obtained by
2 β t −k y(k) − ϕ T (k)θˆ
(21)
where the ‘forgetting factor’ β ≤ 1 diminishes the weight of older data relative to data that was obtained more recently. Typically β is set to a value between 0.95 and 1 depending on the desired rate of adaptation. In off-line applications, β is generally set equal to 1 and the closed form solution is given by
ˆ t) = θ(
t X k=1
that relates the output to the four values at each corner of the cell. Note that the time-step dependency of r1 , r2 , i1 and i2 is not shown explicitly in Eq. (16) for reasons of notational clarity. As long as the cell indices i1 (k) and i2 (k) do not change value between time-steps, then the vector of table outputs y can be obtained according to y = 8θ
3. Recursive least-squares parameter identification
(14)
(13)
and the table output is obtained by interpolating first along the two parallel cell edges in the x1 direction and then between these two points in the x2 direction as shown in Fig. 2. The structure of this interpolation can be seen in the matrix equation
(1 − r2 ) · (1 − r1 ) T 2[i1 ,i2 ] (1 − r2 ) · r1 2[i1 +1,i2 ] y(k) = r2 · (1 − r1 ) 2[i1 ,i2 +1] r2 · r1 2[i1 +1,i2 +1]
Note that the condition on j is obtained by mapping the 2 indices in Eq. (19) according to Eq. (20). With these definitions, Eq. (17) can be used to describe the look-up table response to an input x(k) that varies across the entire 2-dimensional input space of the table.
Eqs. (5) and (17) show how output values are obtained from the table parameters θ and inputs x. In this application, however, the table values θi are to be identified given observed values for the output y and input x. The set of over-determined linear equations represented by Eq. (5) for the 1-dimensional table or Eq. (17) for the 2-dimensional table, with the corresponding definitions for the parameter vector θ and the matrix of regressors 8, can be solved to obtain the unknown table parameters θi using least squares [4,5]. This method corresponds to finding the parameters that minimize the weighted mean squared error between observations and model predictions according to the cost function
The normalized offsets of x1 (k), x2 (k) in each cell dimension are given by the pair of equations r1 (k) =
j = i1 (k) + (nx1 − 1)i2 (k) j = i1 (k) + (nx1 − 1)i2 (k) + 1 j = i1 (k) + nx1 i2 (k) (20) j = i1 (k) + nx1 i2 (k) + 1 otherwise.
! −1 ϕ(k)ϕ T (k)
t X
! ϕ(k)y(k)
(22)
k =1
where the solution to Eq. (22) should yield estimated values for the table parameters of interest as long as the input is sufficiently exciting. The estimates will be unbiased providing the regressors are uncorrelated with the noise process and providing the lookup table model structure admits an exact description of the true system. It is worth noting that in most applications the piecewise linear look-up table model is only an approximation of the (generally smooth) nonlinear system that it is used to represent. The table parameters are therefore likely to be ‘biased’ relative to the exact values of the system at the table grid points, but the resultant ‘bias’ will tend to minimize the overall mean square prediction error between the true process output and that predicted from the linear look-up table approximation—at least
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
ˆ 0) = θinit , θ(
479
P (0) = δ I
(28)
where δ 1 is used indicate uncertainty in the initial parameter estimates if reasonable starting values are not known a priori. 4. An efficient recursive least squares algorithm
Fig. 3. Operating point bias in RLS estimation.
for the data set used during identification. Such behavior is often acceptable or even desirable but it is important to ensure that the data set used during identification is well distributed across and within the cells since the model fit will favor any dominant regions of operation. Fig. 3, for example, illustrates a 1-dimensional slice of a function with high curvature that is being approximated by a linear interpolation table. If the operating condition remains in a narrow range about point (a), then the estimated parameters will be as shown in the left hand plot of the figure. If instead the operating point were constrained to a narrow region around point (b), then the fit shown in the right hand plot of Fig. 3 would result. In order to avoid these issues, the cell grid should be designed to limit the maximum error between the ‘true’ function and its linear approximation and the RLS algorithm should also be fed with data distributed across the cell. If necessary, this distribution can be achieved by windowing the data so that the algorithm only executes when the operating condition has moved significantly from its previous position. In on-line applications, it is usually not possible or desirable to store large data records. It is much more efficient merely to update the table parameter estimates as the next sample of data becomes available. It is also useful to retain the forgetting factor so that the parameters are weighted towards the most recent data. The recursive form, which relates the estimates at time t to those at time t − 1, is obtained by extracting the uppermost elements of the summations in Eq. (22) and re-arranging to give
h
i ˆ t ) = θ( ˆ t − 1) + R (t )ϕ(t ) y(t ) − ϕ (t )θ( ˆ t − 1) θ( −1
T
(23)
where R(t ) =
t X
β t −k ϕ(k)ϕ T (k) = β R(t − 1) + ϕ(t )ϕ T (t ).
(24)
k=1
In order to avoid having to compute the inverse of matrix R(t ) in Eq. (23), it is convenient to define the covariance matrix P (t ) = R−1 (t ) and update this inverse directly. Applying the matrix inversion lemma to Eq. (24) enables P (t ) to be expressed as P (t ) =
1
β
I − L(t )ϕ T (t ) P (t − 1)
(25)
R(t )1θˆ (t ) = f (t )
(29)
where
1θˆ (t ) = θˆ (t ) − θˆ (t − 1) f (t ) = ϕ(t ) y(t ) − ϕ T (t )θˆ (t − 1) .
(30)
As noted above, the general disadvantage of this approach is that the linear system in Eq. (29) involving the R(t ) matrix must be solved at each time step. In this particular application, however, the special structure of the regression vectors (Eq. (11) for 1-D or Eq. (20) for 2-D tables) results in a sparse banded R(t ) matrix that has minimal memory requirements and for which a solution to Eq. (29) can be efficiently obtained as discussed in the sequel. These benefits do not hold for the standard RLS algorithm because the P (t ) matrix does not retain a sparse banded structure upon inversion of the R(t ) matrix in general. 4.1. Recursive least squares approach for 1-dimensional tables
where L(t ) =
One of the drawbacks of the recursive least-squares formulation is the requirement to store and update a P (t ) matrix of dimension (nx1 )2 for a 1-dimensional table, or of dimension (nx1 · nx2 )2 in the 2-dimensional case. This storage requirement therefore increases rapidly with table dimension and size. A 1-dimensional table with 8 data points requires only 64 storage locations for the P (t ) matrix, but a 2-dimensional table with 8 data points on each axis requires 3996. Memory storage requirements can therefore become a constraining issue in real-time embedded applications where computational resources are limited. One way to address this issue is to approximate the covariance matrix by its main and first upper and lower diagonals on the basis that only neighboring cells are strongly correlated to each other [6]. This approximation significantly reduces the data storage requirements to 2n elements as compared to the n2 elements required for the full matrix, where n represents the total number of parameters in the table, but the validity of the approximation has not been proven theoretically. In this paper a different approach is adopted which does not rely upon approximations, but instead exploits the sparse structure of the R(t ) matrix in order to achieve a memory-efficient and computationally-efficient implementation. The advantage of this approach is that the recursive least squares estimate is itself unchanged so all the convergence and stability properties carry over from the standard algorithm. It is only the numeric algorithm with which the estimate is computed that has been redesigned. The new method is based on computing the estimate by rearranging Eq. (23) and using the matrix R(t ), defined in Eq. (24), directly rather than by means of the matrix inversion lemma in Eqs. (25) and (26). Re-arranging Eq. (23) results in the following linear system of equations.
P (t − 1)ϕ(t ) . β + ϕ T (t )P (t − 1)ϕ(t )
(26)
Furthermore, using Eqs. (25) and (26), it can be shown that P (t )ϕ(t ) = L(t ). Eq. (23) may therefore be re-written as
h i ˆ t ) = θ( ˆ t − 1) + L(t ) y(t ) − ϕ T (t )θ( ˆ t − 1) . θ(
(27)
Eqs. (25)–(27) define the well known Recursive Least-Squares (RLS) algorithm. A more detailed derivation, together with a discussion of numerical and convergence issues can be found in [4,5]. The algorithm is initialized by setting
Consider the evaluation of Eq. (24) at time-step t, where ϕ (t) is defined by Eq. (11), in the case of a 1-dimensional look-up table. Inspection of Eq. (11) reveals that the update term for R(t ) in Eq. (24) is a sparse matrix whose only non-zero terms are the 2 × 2 symmetric square matrix shown in Eq. (31).
(1 − r (k))2 ϕ(k)j ϕ (k)j = r (k)(1 − r (k)) T
r (k)(1 − r (k)) . r 2 (k)
(31)
These non-zero elements are positioned starting at location j, j within the R(t ) matrix as shown in Fig. 4. The 4-element submatrix defined in Eq. (31) is positioned on the main diagonal of R(t )
480
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
(1 − r2 )2 (1 − r1 )2 (1 − r2 )2 (1 − r1 )r . 0..
ϕ(k)j ϕ (k)j = r2 (1 − r2 )(1 − r1 )2 r2 r1 (1 − r2 )(1 − r1 ) T
(1 − r2 )2 (1 − r1 )r1 (1 − r2 )2 r12 . 0..
r2 r1 (1 − r2 )(1 − r1 ) r2 r12 (1 − r2 )
...0... ...0... ...0... ...0...
r2 (1 − r2 )(1 − r1 )2 r2 r1 (1 − r2 )(1 − r1 )
.
0..
( − r1 ) ( − r1 )
r22 1 r22 r1 1
2
r2 r1 (1 − r2 )(1 − r1 ) r2 r12 (1 − r2 )
.
0.. 2 r2 r1 (1 − r1 ) r22 r12
.
(34)
Box II.
Fig. 4. Structure of the R(t ) matrix for 1-dimensional tables.
at the location specified by the index j = i(k) that represents the currently active cell. The matrix R(t ) is therefore symmetric and tri-diagonal as also illustrated in Fig. 4. The advantage of this matrix structure is two-fold. First, the R(t ) matrix requires only storage for the (2 ∗ n − 1) non-zero elements, as compared to n2 non-zero elements for the corresponding inverse matrix P (t ), and second, the solution to the linear system in Eq. (29) is efficiently obtained using the Thomas algorithm [7]. The Thomas algorithm is briefly summarized as follows. Let b(i), c (i) represent the elements on main and first upper and lower diagonals of the matrix R(t ), respectively, and let f (i) represent the elements of the vector f (t ). The forward elimination phase of the algorithm is defined as follows: b1(1) = 1; for k = 2 to n, b1(k) = b(k) − c (k) ∗ c (k − 1)/b1(k − 1); f 1(k) = f (k) − c (k) ∗ f 1(k − 1)/b(k − 1); end;
(32)
The backward substitution phase then yields the solution according to:
1θ(n) = f 1(n)/b1(n), for k = n − 1 downto 1, 1θ(k) = f 1(k) − c (k) ∗ f 1(k + 1)/b1(k); end;
Fig. 5. Structure of the R(t ) matrix for 2-dimensional tables.
necessary to accommodate differences in the form of the regression matrices. In the 2-dimensional case, inspection of Eq. (20) shows that the update component of the R(t ) matrix in Eq. (24) has the form as in Box II: In this case, the non-zero elements are grouped as four 2 × 2 symmetric sub-matrices separated from each other by an index offset of nx1 in either dimension. This group of four sub-matrices is positioned starting at element j, j of the R(t ) matrix where j = i1 (k) + (nx1 − 1)i2 (k) corresponds to the ‘lower left’ corner of the active cell as shown in Fig. 5. As the active cell changes, the position within the R(t ) matrix shifts up or down the leading diagonal resulting in the symmetric, banded block tri-diagonal matrix that is illustrated in Fig. 5. The banded block structure of the R(t ) matrix for 2-dimensional tables prevents the direct application of the Thomas algorithm as presented in Section 4.1. To address this issue, the R(t ) matrix is partitioned into a tri-diagonal portion RD (t ), similar to that shown in Fig. 4, and a remainder (R(t ) − RD (t )) that contains the upper and lower banded part. Eq. (29) can then be rewritten as RD (t )1θˆ (t ) = f (t ) − R(t ) − RD (t ) 1θˆ (t ).
(33)
The process for efficient estimation of 1-dimensional table parameters can therefore be summarized as follows:
ˆ 0) = θinit , R(0) = δ I where (i) Initialize algorithm by setting θ( δ < 1 is used indicate uncertainty in the initial parameter estimates if reasonable starting values are not known a priori. (ii) Update the sparse, tri-diagonal matrix R(t ) according to Eqs. (24) and (31). (iii) Compute 1θˆ (t ) and hence update the parameter estimate using Eqs. (29)–(33). (iv) Loop back to step (ii).
(35)
The solution to this equation can now be determined by applying the Thomas algorithm iteratively where 1θˆ (t ) in the right-hand side of the equation is approximated by the estimate obtained at the previous iteration, RD (t )1θˆm (t ) = f (t ) − R(t ) − RD (t ) 1θˆm−1 (t )
(36)
and where m represents the iteration number. The iteration is initiated by setting 1θˆ0 (t ) = 0. Practical application of the method indicates that the solution converges to within 2% of the ‘true’ values obtained by direct solution with the full R(t ) matrix within a few iterations. Apart from this additional iteration step, the method is essentially identical to that presented for the 1dimensional table case.
4.2. Recursive least squares approach for 2-dimensional tables
4.3. 1-dimensional look-up table illustrative example
A similar procedure can be used to estimate the parameters of a 2-dimensional lookup table, although some modifications are
To assess the performance of the various recursive leastsquares algorithms and approximations when applied to table
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
k
25
Known Table
y(k)
20 Estimated parameters, θ
x(k)
481
RLS Estimator k Fig. 6. Simulation block diagram.
15
10
5 25 Output, y(k)
0 2.5
20
3
3.5 4 4.5 Table x-axis values [-]
5
5.5
Fig. 8. Snapshots of estimated table parameters (circles); final estimate (diamond); ‘true’ values (solid).
15
5
0
0
0.2
0.4
0.6
0.8
1 1.2 Time [s]
1.4
1.6
1.8
2
Fig. 7. Measured table input and output data.
identification, simulation data generated from a known table with zero mean unity variance additive noise is used for parameter estimation as shown in Fig. 6. For the sake of illustration, only a small table with seven entries is used where the ‘true’ table values θ0 at nodes x are defined by x = 2.5
3.0
3.5
θ0 = 10
15
19
4.0 21
4.5 22
5.0 21.5
5.5
T
T
.
20
Estimated parameters, θ
10
25 20 15 10 5
3
0.1 3.5
0.05
4
4.5 Table x-axis values [-]
5
5.5
Time [s]
0
Fig. 9. Evolution of estimated parameters with time.
25
(37)
The input is a uniformly distributed random signal designed to span the full input range of the table and the time step is 0.01 s. The resultant input/output time histories, shown in Fig. 7, were then presented to the RLS identification algorithm. In the simulation study described here, θinit is set equal to zero and the R(0) matrix is started with δ = 0.1 on initialization. The value of β is set initially to 0.98. Various snapshots of the estimated table parameters are shown (as circles joined by dotted lines) in Fig. 8. The last of these, corresponding to the parameter values at t = 1 s, is indicated with a diamond symbol. Also shown in Fig. 8 (as a solid line) is the ‘true’ table given by Eq. (37). As shown in the figure, the two match quite closely indicating successful convergence of the algorithm. The time evolution of the estimated parameters is illustrated in Fig. 9. The form of the table can be seen emerging with time from the initial zero parameter values. In the early stages of the table identification, some points adapt more rapidly than others, presumably because the input was exciting these cells more than the others. A more quantitative assessment of the convergence process can be made by slicing Fig. 9 along lines parallel with the time axis that result in the traces shown in Fig. 10. Each trace represents the evolution of an individual parameter in the table, although only four such traces are shown for reasons of clarity. After an initial sharp transient the parameters gradually converge towards their ‘true’ values. The effects of noise are seen in occasional fluctuations away from this ideal.
0.2 0.15
0 2.5
20 Estimated parameters, θi
Measured table input and output [-]
Input, x(k)
15
10 θ7 ...=20 θ5 ...=22
5
θ2 ...=15 θ1 ...=10
0
0
0.2
0.4
0.6
0.8
1 1.2 Time [s]
1.4
1.6
1.8
2
Fig. 10. Evolution of estimated parameters (β = 0.98).
The tradeoff between noise immunity and convergence rate is controlled by the value of the forgetting factor β . The results shown in Figs. 8–10 were obtained with β equal to 0.98. Fig. 11 presents the identification results when the simulation is repeated with β reduced to 0.94. Although the same general trends are observed, the estimates are now more sensitive to noise but should adapt more quickly to genuine changes in the ‘true’ parameter values. The value of β = 0.98 seems more appropriate in this case, though
482
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
25
Estimated parameters, θi
φ d(k) 20
φ m(k)
Engine
rpm(k) load(k)
15 ηvol Table / Map
10 θ7 ...=20 θ5 ...=22
5
z-d
z-d
k
θ2 ...=15 θ1 ...=10
0
0
0.2
0.4
0.6
0.8 1 1.2 Time [s]
1.4
1.6
1.8
z-(d+1)
RLS Estimator
ηvol
2 Fig. 12. 2-D identification block diagram.
Fig. 11. Evolution of estimated parameters (β = 0.94).
care should be taken to tune the forgetting factor in different applications. 5. Identification of a 2-dimensional volumetric efficiency map Because of their computational efficiency and the highly nonlinear nature of automotive powertrain systems, look-up tables are used extensively in automotive engine control and monitoring applications. In order to identify the parameters, it is often necessary to run the engine in steady-state conditions at each point in the table grid, a process which is expensive and time-consuming. The parameter identification itself is also often heuristic, relying on operator judgment for recording the appropriate value from a noisy signal. Even when the recorded value is accurate, errors can also arise if the engine is not operated at exactly the point defined by the table grid. Automated table identification is therefore desirable, both off-line using engine testbed data, and online in order to adapt the table in-vehicle as the engine ages. In either case, both the table input and desired output need to be available during the identification process. In the offline case, such measurements are often available from the engine test-bed instrumentation. In the on-line case, however, the required output is generally only known ‘after the event’. If it were not, a table would hardly be required. Knock spark angle correction, applied by table lookup, is one such example. When knock is detected in a previous firing, the tabulated value can be adjusted to avoid knock in subsequent firings [8–10]. In many cases, the table updates are performed heuristically sometimes simply overwriting previous table values and thereby discarding the information contained in earlier measurements [9,10]. The calibration of engine fueling is another example where a table is used to predict the required fueling for the current cycle as part of the air/fuel ratio feedforward control strategy [11]. After the firing has taken place, and the gas has reached the exhaust gas oxygen sensor, the measured air/fuel ratio can be used to update the table for subsequent cycles. The identification procedure presented in this work is applied to this example. The algorithm is implemented on a rapid prototyping system attached to a Ford 2.0L I4 engine and then used to estimate the parameters of the volumetric efficiency map. A block diagram of the configuration is shown in Fig. 12. In this case the ‘true’ map is represented by the behavior of the engine itself. The inputs are speed and load and the output is volumetric efficiency. The output is not measured directly, but is computed assuming the injectors are well calibrated. The measured equivalence ratio φm , is related to desired equivalence ratio φd according to
φm (k) =
mbf (k − d) ma (k − d)
AFRs φd (k − d)
(38)
where mbf is the base fuel mass, ma the air mass flow, and AFRs is the stoichiometric air/fuel ratio by mass. The lag d represents the transport delay of the engine in engine revolutions. We assume φd is constant, and that the other variables change only slowly. The dynamics of the response are therefore ignored in Eq. (38). The base fuel mass is generally computed from the ratio of the estimated air mass to the stoichiometric air/fuel ratio. Air mass flows are also assumed proportional to volumetric efficiency ηv ol such that the following relationships hold.
φm (k) =
ˆ a (k − d + 1) m ma (k − d)
φd (k) =
ηˆ vol (k − d + 1) φd (k) ηvol (k − d)
φd (k) ηˆ vol (k − d + 1). (39) φm (k) Note that the delay on φd can be ignored since its value is held constant at φd = 1 in this example. ⇒ ηvol (k − d) =
The value of the measured volumetric efficiency is computed ‘after the event’ according to Eq. (39) and the appropriately delayed speed and load inputs are then fed to the RLS algorithm implemented on the rapid prototyping system as shown in Fig. 12. The algorithm is executed at each engine firing where a forgetting factor of β = 0.995 was used. For practical reasons, only four table parameters were monitored at any instant, these being the parameters associated with the current cell whose values were being adapted. Fig. 13 shows the results that were obtained as engine speed (Fig. 14) varied slowly across two cells. Initially the parameter values were constant at their default calibration value. At t = 2 s, however, the RLS algorithm was enabled. The process of algorithm initialization results in an undesirable transient, but the parameters then begin to adapt toward their new values. At t = 7 s, engine speed crosses the 1600 rpm level which moves the operating condition into the next cell of the table. The two lower speed parameters 25,4 , 25,5 represented by the dotted traces in Fig. 13 for 0 < t < 7 are therefore dropped from the monitored list, and replaced by the two higher speed parameters 27,4 , 27,5 associated with the new cell (dotted traces t > 7 s). The parameters, 26,4 , 26,5 , are common to both cells, so there is no change in identity for these traces. As the operating condition enters the new cell, the two new parameters 27,4 , 27,5 adapt rapidly from their initial values at t = 7 s to levels around 0.9, 0.95 respectively. However the degree of uncertainty or ‘noise’ in the estimated values, seems higher for these parameters than for the parameters 26,4 , 26,5 which now define the low speed end of the cell. This is believed to be because
J.C. Peyton Jones, K.R. Muske / ISA Transactions 48 (2009) 476–483
1.3 Θ5,4 , Θ7,4
1.2 Estimated parameters, θi
6. Conclusions
Θ6,4
1.25
Θ6,5 Θ5,5 , Θ7,5
1.15 1.1 1.05 1 0.95 0.9 0.85 0.8
0
5
10
15
20
483
25
An efficient recursive least-squares identification algorithm in terms of both computation and storage is presented for the initial determination and adaptation of the parameters in linear look-up tables. This algorithm has been applied successfully to 1- and 2dimensional problems in simulation and in real time. The real-time application concerned the estimation of volumetric efficiency for engine fueling, but the method can be applied in any case where input and (generally delayed) output signals are available. For the method to be robust, care must be taken to ensure that the data is well distributed over the range of operation defined by the table and that the table grid is well designed relative to the curvature of the function being approximated. Further work into developing windowing and weighting strategies to address these issues is ongoing.
Time [s]
Acknowledgments
Fig. 13. Evolution of estimated ηv ol parameters.
The authors gratefully acknowledge financial and technical support from Ford Motor Company, Johnson-Matthey, ExxonMobil, and Mr. E. Barry, ChE ‘54.
1700
Engine Speed [rpm]
1650
References
1600
1550
1500
1450
0
5
10
15
20
25
Time [s] Fig. 14. Engine speed time history.
engine speed fell mainly in the lower end of the cell, resulting in better estimates of the nearby parameters (26,4 , 26,5 ), compared to the others which are further away.
[1] Leenaerts D, van Bokhoven W. Piecewise linear modeling and analysis. Berlin: Springer; 1998. [2] Bellman R, Roth R. Curve fitting by segmented straight lines. Amer Statist Assoc J 1969;64:1079–84. [3] Nelles O. Nonlinear system identification. New York: Springer-Verlag; 2000. [4] Ljung I. System identification: Theory for the user. 2nd ed. Upper Saddle River (NJ): PTR Prentice Hall; 1999. [5] Soderstrom T, Stoica P. System identification. UK: Prentice Hall International; 1989. [6] Peyton Jones J, Muske K. Automatic calibration of 1 and 2-D lookup tables using recursive least squares identification techniques. In: 2007 SAE world congress. Paper 2007-1333. 2007. [7] Press W, Teukolsky S, Vetterling W, Flannery B. Numerical recipes. 2nd ed. Cambridge (UK): Cambridge University Press; 1992. [8] Wu G. A table update method for adaptive knock control. In: 2006 SAE world congress. Paper 2006-0607. 2006. [9] Entenmann R, Unland S, Torno O, Haeming W, Rothhaar U, Surjadi I. et al. Method for the adaptive knock control of an internal combustion engine. US Patent No. 5645033. 1997. [10] Entenmann R, Unland S, Haeming W. Process for controlling knocking in internal combustion engine. US Patent No. 5243842. 1993. [11] Knorad JA, Freeland M, Czuhai BS. Method and system for controlling flexible fuel vehicle fueling. US Patent No. 5467755. 1995.