No Subject

[Mike Klopfer]

I'm interested in coming up with a fuel control algorithm. Perhaps eventually
even a semiautomated means of arriving at a good algorithm given an arbitrary
SI engine. Unfortunately I don't have access to a fuel injected engine. So I
was wondering if there is enough interest by people who do have the capability
of performaing measurements on an engine to perhaps collaborate on a control
algorithm. In an attempt to figure out what this might involve I've been
looking at some of the SAE "Electronics Engine Controls" books. I'm trying to
see if I can use the information from the
paper SAE 940373 "Adaptive Air-Fuel Ratio Control of a Spark-Ignition
Engine" in the book "Electronic Engine Controls 1994" to derive a control
algorithm. This was for a single cyliner engine so I'm not sure how much
difficulty would be involved in extending it to a multicylinder engine.
Perhaps someone could suggest a better set of equations. Anyway the these
researchers acheived <.5% rms deviation of the air-fuel ration for
step changes of throttle angle of 10 degree. This engine model tells nothing
about the power of efficiency of the engine. Perhaps that would be of more
interest but perhaps it is also more difficult to arrive at an accurate
model. Below is a description of the
equations and some random ideas on how to begin to use them to come up with
a control algorithm.

The fuel dynamics equations are:

fp(k+1)= (1-f_b) fp(k) + (1-f_a) on(k)
phie(k+1)= (f_b fp(k) + f_a on(k))/ (injSt(k)- t0)       /* equation A */
phid(k+1)= phie(k)
phim(k+1)= g2 phie(k) + g1 phid(k) + g- phim(k)

on= inj - t0
g0= exp(-T/te)
g1= exp(-m T/te)-exp(-T/te)
g2= 1-exp(-m T/te)
m= 2-td/T

fp = fuel puddle mass equivalent time
phie= unit cycle delay of equivalence ratio of cylinder air-fuel ratio
phid= unit cycle delay of phie
phim= UEGO sensor measurement of equivalence ratio
inj= injector pulse width
injSt= steady state pulse width yielding air-fuel equivalence ratio of 1
t0= injector no flow time
f_a= fraction of injected fuel entering cylinder
f_b= fraction of puddle mass entering cylinder
T= engine cycle period
td= transport delay
te= UEGO sensor time constant

The system parameters are f_a, f_b, t0, te, td, and injSt. These discrete
time equations appear to have been derived from a continuous time equations 
much like the following from John S Gwynne <jsg at coulomb.eng.ohio-state.edu>:


>                                      ----> C2 --->|
>                                      |            |
>desired F/A --> + --(error)--> C1 --->|            + ---> C4--->exhaust F/A
>                                      |            |         |
>                ^ (-)                 ----> C3 --->|         |
>                |                                            |
>                |                                            |
>                 ----- [sensor] <-------- C5 <---------------
> 
>where C1 = Control law
>      C2 = Inlet manifold (fast fuel) = 1/(t_1 s + 1)
>      C3 = Inlet manifold (slow fuel) = 1/(t_2 s + 1)
>      C4 = time delay (engine rotation...) = e^(-s T)
>      C5 = Sensor lag (exhaust gases mixing in the manifold)
>         = 1/(t s + 1)
>      sensor = sensor characteristic; Vout as a function of F/A
 
The authors seem
to mainly consider the case where injSt, f_a and f_b are functions of throttle 
angle and RPM. It seems possible that air temperature, coolant temperature and
air pressure could affect these as well.

The following control law is derived:

on(k)= (injSt(k)-t0) phiDes - 
       (f_b(k)/f_a(k)) {efp(k)-[(1-f_a(k))/f_b(k)](injSt(k)-t0) phiDes}

     = [(injSt(k)-t0) phiDes/f_a(k)]-[(f_b(k)efp(k))/f_a(k)]

Where efp(k) is the estimated fuel puddle mass and phiDes(k+1) is the desired
value of phie . The authors estimated efp(k)
using an extended Kalman filter. I didn't follow there description of how
they derived the above equation (something aoubt a Ricatti equation) but since
it just seem to be equation A above solved for on(k) I didn't pursue it. The 
following are the equations I have found
for the Kalman filter which seems to apply to the discrete nonstationary
equations given above. However I'm not sure if this is the "extended Kalman 
filter". For a system  of the following form where y(k) is the measured 
variable:

x(k+1)= A(k)x(k)+ B(k)u(k) + v(k)
y(k)= C(k)x(k)+w(k)

where v(k) and w(k) are zero mean random variables with variance
E(v(k)v(k))= Q(k) and E(w(k)w(k))= R(k)

we can estimate x(k+1|k) (i.e. x at time k+1 given measurements up to and
includeing time k) using the following equations:

x(k+1|k)= A(x)x(k|k) + B(k)u(k)
x(k|k)= x(k|k-1)+K(k)[y(k)-C(k)x(k|k-1)]
P(k+1|k)= A(k)P(k|k)A'(k) + Q(k)
K(k)= P(k|k-1)C'(k) inverse[C(k)P(k|k-1)C'(k)+R(k)]
P(k|k)= P(k|k-1)-K(k)C(k)P(k|k-1)= [I-K(k)C(k)]P(k|k-1)

Letting x'= [fp phie phid phim] and y= [phim] the following matrices are 
obtained from the above equations:

          
          1-f_b(k)             0     0     0
          f_b(k)/[injSt(k)-t0] 0     0     0
A(k)=     0                    1     0     0
          0                    g2(k) g1(k) g0(k)

          1-f_a(k)
B(k)=     f_a(k)/[injSt(k)-t0]
          0
          0

C(k)= [0 0 0 1]


The discrete equation for the UEGO sensor output was obtained from the 
following continuous time equation:

 te phim'(t) + phim(t) = phie(t-td)

The continuous time solution to this equation is the following:

phim(t)= exp(-(t-t_0)/te) phim(t_0) + 
         integral{d,t_0,t,{exp(-(t-d)/te) phie(d-td)/te}]

where integral[u,v,x,y] is the integral with respect to variable u of the
expression y from v to x.

Evaluating this as discrete times t_0,t_1, ... ,t_k,t_k+1, ... the following
equations can be derived:

phim(t_k)= exp(-t_k/te) { exp(t_0/te) phim(t_0) +
                          integral[d,t_0,t_k,{exp(d/te) phie(d-td)/te}]}

phim(t_k+1)= exp(-t_k+1/te) { exp(t_0/te) phim(t_0) +
                              integral[d,t_0,t_k+1,{exp(d/te) phie(d-td)/te}]}

           = exp(-(t_k+1 - t_k)/te) phim(t_k) +
             exp(-t_k+1 /te) integral[d,t_k,t_k+1,{exp(d/te) phie(d-td)/te}]

Assuming that phie(t)= phie(t_k) for t_j =< t < t_j+1 for all j (i.e phie is
constant between sampling times, and assuming that td < t_k+1 - t_k the 
following can be derived:

exp(-t_k+1 /te) integral[d,t_k,t_k+1,{exp(d/te) phie(d-td)/te}] =
   phie(t_k-1)[exp(-(t_k+1 - t_k - td)/te) - exp(-(t_k+1 - t_k)/te)] +
   phie(t_k)[1- exp(-(t_k+1 - t_k - td)/te)]


These equations for phim (letting T= t_k+1- t_k) are the same as those from the 
paper except that in this case m= 1- td/T instead of m= 2- td/T. I'm not sure
what is the source of this discrepancy but since they don't show the derivation
of their equations its hard to figure out. Unless someone could point out an
error in my derivation.

The transport delay (td) was estimated by controlling the throttle plate
angle with square wave pattern. The UEGO sensor output was then used to
estimate td using some nonlinear equation fitting method.

Perhaps a similar approach could be applied for a multicylinder engine. For
example one cylinder could be given a rich air fuel mixture relative to the 
other cylinders. If anyone is interested in making the measurements of
EGO, throttle angle, RPM, injection pulse width and frequency and any other
parameters of a running engine available I would like to try to figure out how 
to extract the model parameters from such a data log.

One part of the transient behavior that isn't modeled in the above equations
is the manifold dynamics. It would be interesting to see whether and under
what circumstances this portion of the engine behavior significantly impacts
the air-fuel ratio. The SAE papers 930856 and 910258 suggest that the manifold 
dynamics does significantly effect the air-fuel ratio during transients. The
latter is a nice introduction to what is apparently a fairly simple and 
accurate mathematical model of spark ignition engines called the Mean Value
Engine Model (MVEM). I didn't see anything
in these models that would help in choosing the ignition timing though.