Inverse NARMA a robust control method applied to SI engine

*Corresponding author.Tel.:+44-151-794-4828;fax:+44-151-794-4892.E-mail address:shenton@(A.T.Shenton).0967-0661/02/$-see front matter r2002Elsevier Science Ltd.All rights reserved. PII:S0967-0661(02)00111-9potentially more cost effective than feedforward due to the high cost of load sensor hardware.In any case,feedforward control is usually implemented with an element of feedback to deal with the uncertainty present in the disturbance loads and the varying or unknown system dynamics.Due to the increasingly stringent legislation on emission levels and fuel economy together with the very competitive performance and cost requirements of the automotive industry,the idle-speed control problem has received a great deal of attention in the literature (see survey paper by Hrovat and Sun (1997)and the references therein).In this paper a NARMA filter is proposed for the direct inverse filter.Previous engine modelling based on NARMA/NARMAX/NARX models includes the very first NARMAX application carried out by Billings,Chen,and Backhouse (1989),in which a turbocharged automotive diesel engine was modelled.Another inter-esting approach to using NARMAX modelling in engine control was presented by Glass and Franchek,1999.A NARX model was also applied to the idle-speed control problem by De Nicolao,Rossi,Scattolini,and Suffritti (1999),using both spark and air channels for feedback control.After an overview of the inverse control scheme in Section 2,the structure of a suitable NARMA filter is considered in Section 3and then the selection of a suitable identification method for the filter is discussed in Section 4.The experimental set-up and the physical processes involved in the idle dynamics are discussed in Section 5.To aid with the determination of the inverse-NARMA (INARMA)filter,a simple but practical heuristic structure detection method is presented in Section 6to enable fast elimination of inconsequential terms.In Section 7the effectiveness of the linearisation using the resulting inverse is determined by linear identification of the system,pre-compensated with the identified INARMA filter.Linear controllers for the linearised system are accordingly presented in Section 8.As a final validation of the proposed technique,in Section 9the engine system,pre-compen-sated by its inverse,is compared experimentally onan electric dynamometer to the system without inverse compensation.2.Direct inverse controlDirect inverse control (DIC)was originated by Horowitz (1981)in the context of non-linear QFT design.The idea of DIC is to pre-compensate or post-compensate a non-linear plant G directly by its exact or approximate inverse dynamics L :This is done with a view to obtaining the advantage of linearising any non-linear dynamics and also to providing a unity path.A feedback controller K ;round the inverse compensated plant is then implemented using linear design methods.Fig.1shows the pre-compensated DIC set-up for the constant engine-speed regulation disturbance-rejection problem considered in this paper,in which torque disturbances add to the system speed output via the load dynamics G d :In the case of the pre-compensated system considered here,the unity path between pre-filter input and plant output is used for implementing output regulation and tracking systems.The pre-filter is then required to constitute a right-inverse (Kotta,1995)for the plant.In the alternative post-compensated scheme the unity path between the plant input and a post-filter output can be used for control effort limiting (saturation avoidance)systems.The post-filter must then constitute a left-inverse (Kotta,1995)for the plant.DIC is contrasted to inversion by feedback linearisa-tion (Marino &Tomei,1995;Hunt &Meyer,1997),in which state or output feedback is used to linearise the system,and in which there then exists an inner linearising feedback loop.By contrast in the DIC approach,the single feedback loop design allows direct consideration of the complete system robust-perfor-mance and the robustness-performance trade-off (Hor-owitz,1982).In DIC design,this design trade-off may be directly addressed either by a linear QFT approach (Horowitz,1982;Borghesani,Chait,&Yaniv,1994)or by any alternative linear robust performance techniquesFig.1.The idle-speed Direct Inverse Control scheme.A.P.Petridis,A.T.Shenton /Control Engineering Practice 11(2003)279–290280(Chiang&Safonov,1992;Balas,Doyle,Glover, Packard,&Smith,1991;Besson&Shenton,2000).In its simplest conception DIC may be used to linearise a(i.e.uncertain)certain plant and thereby eliminate the amount of uncertainty necessary to account for the non-linearity when using a linear control method.In a more realistic setting DIC is used to reduce the amount of uncertainty necessary not only to account for the non-linearity but also for any plant or disturbance uncertainty.A DIC method is consequently required to determine not only the approximate inverse but also the corresponding reduced uncertainty.In Horowitz’s non-linear QFT technique(Horowitz,1981), a non-linear plant model in input–output format is inverted and the QFT plant templates of the inverse compensated system are then found based on parametric uncertainty about a‘known’nominal non-linear model. Because the nominal model is somehow‘centred’1in the uncertainty,the resulting template sizes are reduced and the robustness-performance trade-off is then less tight. The resulting non-causal inverse is combined with a linear controller with an excess of poles over zeros to provide an overall causal control scheme.This approach requires a minimum-phase(MP)plant and a detailed parametric model,which accurately represents the system uncertainty.Unfortunately,both these assump-tions are often not valid in practical situations such as in engine control.Previous work by the authors(Petridis&Shenton, 1999,2000,2002;Petridis,2000;Shenton&Petridis, 2002)proposed related direct-inverse techniques where-by identified continuous time non-linear input–output or state-space plant models were mathematically in-verted resulting in a non-causal inverse.Such inversions can result in highly complex inverses for simple forward systems(Petridis,2000).As mentioned,implementing such continuous methods eventually requires discretisa-tion.In any case,the physical processes in the engine although broadly understood are actually complex and difficult to model without resorting to experimentation on the candidate engine.One is consequently led to consider the use of NARMA models.The discrete time technique proposed in this paper has the advantage that the system’s inverse is modelled directly hence bypassing the need to develop a forward model,be it phenomen-ologically based or otherwise,and bypassing subsequent mathematical inversion.A further motivation for this approach has been the fact that the analytic inverses of NARMA forward models have been found to be unstable.The proposed inverse technique identifies a non-linear inverse plant model directly from reversed input–output (i.e.output–input)data.Approximate LTIE plant templates(uncertainty regions)which neglect the time-varying nature of the linear spot point models are then found by means of linear system identification.This results in collections of linear spot point models of the open loop direct-inverse compensated plant,for a range of different operating conditions.These templates are conveniently considered as circles in the Nyquist chart (unstructured uncertainty),representing the uncertainty about the nominal linear model.In the subsequent linear controller design approach presented here,the time varying nature of the switching between the linear spot point models is(conservatively)taken into account by considering the circle theorem of Sandberg(1964).3.Filter structure and relative degreeTo implement the inverse control strategy,we require to identify the inverse system as a pre-compensating controlfilter.Since we propose the use of parameter estimation for this identification we must establish candidatefilter structures including the required relative degree.If the relative degree of a plant is unknown,then that of an INARMAfilter may be determined on a bestfit basis by comparing thefit of trialfilters with different relative degrees.If the relative degree of a certain plant is d P and the number of sample-times in the pure time-delay is d D;then the relative degree of the INARMA filter,d I must be d I¼Àd Pþd D to exactly invert the plant.Thus if this difference is known,for the most accurate results the trialfits should be with INARMA filters of this relative degree.Since real plants require d PÀd D>0this means such an INARMAfilter is non-causal.Nevertheless,in simulation(during identification of the inverse),this may be readily handled by advancing the plant output data d P sample times to make it causal.Then once this causal inversefilter with shifted data is identified,thefilter may itself be advanced d PÀd D sample times to account for this and recover the true non-causalfilter equations.To use such a non-causal inversefilter in thefinal control implementation, the linear controller is then designed with sufficient excess poles over zeros to produce an overall causal combined controller-filter.Several robust linear control design techniques which directly address the possibility of designing for such acausal processes are currently available,including QFT loop-shaping(Horowitz,1982) and robustfixed-orderfixed-structure(includingfiltered PID)methods(Besson&Shenton,1997,1999,2000; Shenton&Besson,2000).Although increasing(making less negative)the relative degree of thefilter beyond any true relative degree is liable to result in loss of accuracy and overfit(Ljung,1987)if thefilter relative degree is increased to or beyond zero,thefilter is then causal and this,although not necessary,may sometimes have obvious convenience.1Currently a trial-and-error process.A.P.Petridis,A.T.Shenton/Control Engineering Practice11(2003)279–290281In order to demonstrate the capability of the method, and to contrast to the approaches in Petridis(2000)and Petridis and Shenton(2002),for this paper a(less exact) compensator of zero relative degree was chosen.In the case of the direct-inversefilter for the idle-speed control scheme,the input is the desired engine speed,here denoted u:The inversefilter output is the duty cycle to the ABV,denoted y:A non-linear NARMA model structure consisting of35parameters was initially selected to represent the inverse system without time delay:y t¼p1þp2u tþp3u tÀ1þp4u tÀ2þp5u2t þp6u2tÀ1þp7u2tÀ2þp8y tÀ1þp9y tÀ2þp10y2tÀ1þp11y2tÀ2þp12u t y tÀ1þp13u t y tÀ2þp14u t y2tÀ1þp15u t y2tÀ2þp16u tÀ1y tÀ1þp17u tÀ1y tÀ2þp18u tÀ1y2tÀ1þp19u tÀ1y2tÀ2þp20u tÀ2y tÀ1þp21u tÀ2y tÀ2þp22u tÀ2y2tÀ1þp23u tÀ2y2tÀ2þp24u2ty tÀ1þp25u2t y tÀ2þp26u2ty2tÀ1þp27u2ty2tÀ2þp28u2tÀ1y tÀ1þp29u2tÀ1y tÀ2þp30u2tÀ1y2tÀ1þp31u2tÀ1y2tÀ2þp32u2tÀ2y tÀ1þp33u2tÀ2y tÀ2þp34u2tÀ2y2tÀ1þp35u2tÀ2y2tÀ2:ð1ÞThe35terms were chosen to include all combinations up to second order bilinear expressions.This was considered reasonable since in the forward case, previous successfully identified non-linear state-space and NARMA engine models(Billings,Chen,&Back-house,1989;Hariri,Shenton,&Dorey,1998;Petridis& Shenton,2000)have consisted of up to second order bilinear terms.4.INARMA identification methodThough least-squares(LS)identification methods (Billings,Korenberg,&Chen,1988)are highly efficient in producing a model when the parameters appear linearly as in(1),it is possible that the inverse-systems which are identified by LS are unstable.In fact a primary motivation for the current work was the discovery that forward NARMA models developed for the engine had locally unstable inverses.Of course,an exactly identified inversefilter for a non-minimum phase (NMP)process should be unstable,but this will be expected to produce uncontrollable unstable modes with consequent unbounded control inputs when used as a pre-compensator.Accordingly,it is required in such cases to approximate the system inverse by a stable system.It is thereby required to have an identification method which can be guaranteed to produce a stable model of the inverse system when LS does not deliver one. Conveniently however,in view of the fact that the experimental test signals and measured values are always bounded,the input–output data sets for the reverse identification(output–input)always represent bounded-input–bounded-output(BIBO)data sets and so the process can be readily identified by parameter estimation methods.A parameter estimation method capable of imposing stability constraints is thus required to identify stable approximate systems.With such a method it follows that in the case of the non-linear plant displaying unstable zero dynamics,it is thereby possible to ensure the stability of the INARMA filter at the expense of loss of accuracy in the inversion and consequent loss of advantageous effect from the pre-compensator.Increased difficulties in control are of course always to be expected in NMP systems,and some loss in benefit from the inverse control will consequently be unavoidable.Whether there is increased degradation over and above that inherent in the NMP characteristics and how to quantify a priori any such increased degradation,if it exists,is an open question.Never-theless,an assessment of the net benefit of the INARMAfilter as measured by the reduction in linear uncertainty to account for the system non-linearity is always an outcome of the proposed scheme.In the current paper the stable identification method is implemented by prediction error minimisation using a genetic algorithm(GA)based search technique in which model simulation is used to evaluate an RMS error performance function.The result of an additional local stability test(here a simple step response)is combined in the cost function to penalise unstable simulations,so as to eliminate them at each generation.The search method was implemented using the Matlabgenetic algorithm optimisation toolbox(GAOT)(Houck,Joines,&Kay, 1995).An additional feature of the GA technique used is an heuristic approach to structural identification.This was to enable afilter structure of reasonable quality but without excessive complexity to be found.Such lower orderfilters can represent significantly less software and hardware implementation overhead in practical engine management systems.5.Experimental set-up and engine processesThe experimental work described here was carried out on the low inertia electric engine dynamometer at the University of Liverpool(Dorey,Maclay,Shenton,& Shafiei,1995).The subject engine was a4-cylinder Ford Zetec1:6L16valve engine with sequential port fuel injection.The Matlabreal time workshop(RTW)was used to interface the dynamometer using dSpace-Auto-box to enable a Simulink representation of the inverse system to be applied in real time as a pre-filter to the engine.A.P.Petridis,A.T.Shenton/Control Engineering Practice11(2003)279–290 282The engine process considered is a single-input–single-output(SISO)system.The component subsystems are comprised of the idle-speed control valve,the manifold plenum,the engine pumping dynamics,combustion and rotational dynamics.The control input is the(non-dimensional)duty cycle to the air-bleed valve(ABV), and the system output is the engine speed(rpm).There is known to be a pure time delay T d(s)between the ABV input and the engine speed output.This delay represents the transport-delay between the mass-airflow into the inlet manifold and the indicated torque during the intake-to-combustion stroke.In a simple state-space representation of the process,the mass-airflow rate into the intake manifold can be modelled as a function of ABV opening.The rate of change of intake manifold pressure at constant intake manifold temperature can be modelled as the difference between the mass-airflow rate into the manifold and the mass-airflow rate out of the manifold.The mass-airflow rate out of the intake manifold can be modelled as a function of intake manifold pressure and engine speed.The rate of change of engine speed can be modelled in terms of the mass-air flow rate out of the manifold and engine speed.Refer to (Petridis&Shenton,2000)for a mathematical descrip-tion of this process,and an identified phenomenologi-cally based forward model.6.INARMA model identificationTo ease the identification process the intake-to-combustion delay T d was identified in the forward sense using step response tests.The value of T d was thus found to be5T s s;where T s represents the nominal sample time of33ms corresponding to1801crank rotation at a nominal engine idle speed of880rpm.In thefirst stage of the identification the effect of T d was then removed by shifting the input–output data byfive samples.Initially,the GAOT-based search method was used to identify the INARMA model parameters p1–p35:The excitation signal applied to the ABV was designed to excite the expected envelope of engine idle speed operating conditions ranging from between approxi-mately500and1200rpm.The shape of the excitation signal can be seen in Fig.3in the lower subplot as the recorded ABV duty signal.Once identified,a simplified structure was then detected using a technique to select the most significant terms.The survey paper by Haber and Unbehauen (1990)summarises a number of structure identification methods based on step and impulse tests,frequency response measurements,correlation analysis and time response data.The method employed in this paper is based on simulated time response data,in which all possible regressions are atfirst considered.Each parameter was excited in turn by a percentage of the identified parameter coefficient and the effect on the system’s output rms error was monitored.A sinusoidal signal was used to excite each parameter,the magnitude of which was30%of each parameter value.Fig.2a shows each response superimposed,where the effect of each parameter on the inverse output y can be seen. Fig.2b shows a histogram for each parameter p1–p35 showing the relative effect each has on the rms error between the measured data and the inversemodel Fig.2.(a)Effect of parameter excitation on system output.(b)Parameter significance.A.P.Petridis,A.T.Shenton/Control Engineering Practice11(2003)279–290283prediction.The original rms error for the 35-parameter model was 0.01479over 68s :The 10terms (above the line y ¼0:014795in Fig.2b )which had the most influence on the system output due to parameter excitation were selected to form the new model structure.After re-identification of the new structure using GAOT,the best rms error achieved was 0.0153over 68s :The loss in model quality from the 35-parameter structure was considered acceptable when balanced with the reduction in model complexity.The simplified 10parameter model is y t ¼p 1þp 2u t þp 3u t À1þp 5u 2tþp 6u 2t À1þp 7u 2t À2þp 16u t À1y t À1þp 17u t À1y t À2þp 18u t À1y 2t À1þp 26u 2t y 2t À1;ð2Þfor which the associated identified parameter coefficients are given in Table 1.The model was validated using previously unseen experimental output–input data by inputting the mea-sured speed data to the inverse system and comparing the resultant ABV output to the recorded ABV from the data set as shown in Fig.3.7.Linear model identification and compensator assessmentLinear identification was carried out,firstly for the uncompensated engine between the ABV duty cycle as input and the engine speed as output,and secondly for the inverse compensated system for which the input signal is the engine speed demand and the output is the actual engine speed.Pseudo-random binary sequence (PRBS)signals were used to excite a wide range of frequencies.The magnitude of the PRBS ABV duty cycle signals was 0.015,which was added to and subtracted from the nominal ABV duty,which at the desired idle speed was 0.34.The period of the signal was designed based on the step response behaviour of the system yielding a 6-shift register PRBS expanded 5times,based on the sampling time T s ¼33ms :The MatlabSystem Identification toolb ox (Ljung,1991)was used to implement the ARMAX identification algo-rithm.By this means,collections of linear models were identified at different operating conditions characterised by coefficients of rational transfer functions with time delay.The normalised frequency response of these models,consisting of eight such identified models in each case,Table 1Identified parameter coefficients Parameter ValueUnits p 10.22521p 20.2633Â10À31=rpm p 3À38.303Â10À61=rpm p 5 5.364Â10À91=rpm 2p 695.131Â10À91=rpm 2p 718.196Â10À91=rpm 2p 16À0.3078Â10À31=rpm p 17À0.2473Â10À31=rpm p 1821.138Â10À61=rpm p 26À40.408Â10À91=rpm 21020304050607040060080010001200I n p u t u (R P M )Time (s)102030405060700.250.30.350.40.45O u t p u t y (A B V D u t y )Time (s)Fig.3.Inverse model validation (a)measured speed (b)recorded vs.predicted ABV duty.A.P.Petridis,A.T.Shenton /Control Engineering Practice 11(2003)279–290284are shown in Fig.4.It can be seen that a reduction in the spread of the inverse compensated models is apparent about the desired crossover frequency range of about 1–4rad =s :However,visually comparing the uncertainty in a Bode plot is not a reliable or easily quantifiable method.In fact,any direct comparison is not a straightforward task and requires some engineering judgement.For a fair assessment,the system’s gain and phase uncertainty must be compared at each frequency but particularly around crossover frequencies.In the engine problem the system uncertainty is unstructured around these frequencies.Accordingly,the uncertainty is represented by fitting minimal disks to the transfer function loci at discrete frequencies in the complex plane (Nyquist diagram).Then to enable a direct comparison of the two systems at similar frequencies,we reason it is acceptable to apply any certain MP linear compensation to the open loop (uncertain)system since this can be exactly removed by its inverse as part of the controller.Accordingly,such compensation was applied to the uncertain Nyquist loci to cause the vector margin (VM)of each system to meet the uncertainty disks at the same approximate crossover frequency.The VM is then an indication of the relative amount of robust stability of the open loop models which is of course related to the achievable robust performance of the closed loop system.Fig.5shows disks spanning the frequency range 1.3–3:9rad =s :The VM of both systems intersects the disks at 3:027rad =s :It can be clearly seen that at these frequencies,the spread of the uncertainty of the inverse compensated system has been significantly reduced when compared to the system without inverse compen-sation.The VM of the inverse compensated system is 0.6133compared to the uncompensated system VM of 0.412.8.Controller designA robust idle-speed controller is required to achieve a number of performance targets.These are specified in terms of minimal speed deviation,fast settling time,small overshoot and zero steady-state error in the face of disturbance torque loading,sensor noise,environmen-tally produced parameter variation and engine-to-engine variation.In order to account for the uncertainties in engine characteristics,the closed loop system is required to have significant gain,phase and vector margins.To enable a fair comparison between the systems both with and without the inverse compensation,the con-troller design in this paper is optimised for each system using the same technique.A weighted sensitivity parameter space method (Besson &Shenton,1997,1999,2000;Shenton &Besson,2000)was used to carry out controller design.This method allows explicit design for robust performance in a non-conservative manner for the unstructured non-parametric uncertainties and pure time delays which are present.To quantify the plant uncertainty in the controller design process,disks are fitted to the Nyquist plant loci as described in the previous section.One hundred and fourty-one such disks were fitted at discrete frequencies ranging from 0.1to 100rad =s :Following Jayasuriya and Franchek (1994),the desired time domain response was used to derive the frequency weighting functions.10-1100101Frequency (rad/s)1011010Frequency (rad/s)10-1100101Frequency (rad/s)10-1100101M a g n i t u d e10-1100101M a g n i t u d e−200−150 −100 −50P h a s e (d e g )−200 −150 −100 −500P h a s e (d e g )Fig.4.Normalised ABV to speed frequency response for the uncompensated engine (lhs),and the inverse compensated engine (rhs).A.P.Petridis,A.T.Shenton /Control Engineering Practice 11(2003)279–290285The designed PID controller gains are as follows for the system without and with non-linear compensation,respectively:K 1¼2:068Â10À4s 2þ4:417Â10À4s þ4:781Â10À4s 2þ10s;ð3ÞK 2¼5:587s 2þ31:520s þ40:770s þ10s:ð4ÞThe controllers were applied to the SISO control loop for each case as shown in Fig.1.For regulation about a fixed set-point without non-linear compensation,the plant G is represented by the collection of linear transfer function models,controlled by controller K with L ¼I :When inverse pre-compensation is included in the system,the plant becomes G L where L represents the inverse pre-filter.The systems are subject to disturbance loads as shown,where G d represents the input torque dynamics.Attenuation of sensor noise before the combustion cycle frequency is also required.The disturbance load dynamics G d were identified using the dynamometer to excite the engine speed with a PRBS torque input signal.The torque loading is represented by a step input through the identified (Petridis,2000)second order transfer function:G d ðs Þ¼À0:0046s þ0:281s þ9:2671s þ18:9697:ð5ÞThe magnitude of the input step disturbance to be applied to the above transfer function is 24,889,which is equivalent to 18:4Nm of applied torque.This was chosen to represent the following (actual)loads appliedto the engine during controller testing:2:6Nm cooling fan,1:3Nm headlights,2:6Nm heater/air conditioning fan and 11:7Nm alternator loading via a resistor load bank.The critical disk of Sandberg’s circle theorem (Sand-berg,1964)is introduced at this stage.The critical disk is an extra (conservative)robust stability constraint which takes into account the rate and degree of change from one linear-time-invariant (LTI)model to another within the uncertainty disks.The critical disk is located in the complex plane and must be avoided by the uncertain loop locus.The size of the disk is related to the change in the gain of the nominal transfer function model G nom ;i.e.,a ¼j G j minj G j nom;b ¼j G j maxj G j nom;ð6Þwhere the minimum and maximum gains are j G j min and j G j max ;respectively,and j G j nom is the nominal (mean)gain.a and b translate into the real axis points À1=a and À1=b which determine the diameter of the critical disk,which is centred on the real axis.Constructing the critical disk based solely on the change of the gain of the LTI models is justified from an engineering point of view because excepting for very low frequency and DC gain,the dynamics of the spot point models have similar shaped loci (see Fig.4).Fig.6shows the unperturbed loop function L ¼GK for the system without inverse compensation.The critical disk is then determined by a ¼0:605and b ¼1:546as shown.The loop function for the inverse compensated system is shown in Fig.7.The critical disk−2−1.5−1−0.500.5−1−0.8−0.6−0.4−0.2Real axisI m a g i n a r y a x i s−2−1.5−1−0.500.5−1−0.8−0.6−0.4−0.2Real axisI m a g i n a r y a x i sFig.5.VM of uncertainty disks (a)without and (b)with non-linear inverse compensation.A.P.Petridis,A.T.Shenton /Control Engineering Practice 11(2003)279–290286。