A low-complexity explicit MPC controller for AFR control

Control Engineering Practice 42 (2015) 118–127

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

A low-complexity explicit MPC controller for AFR control Marek Honek a,n, Michal Kvasnica b, Alexander Szűcs b, Peter Šimončič a, Miroslav Fikar b, Boris Rohal’-Ilkiv a a Institute of Automation, Measurement and Applied Informatics, Slovak University of Technology in Bratislava, Faculty of Mechanical Engineering, Námestie slobody 17, 812 31 Bratislava, Slovakia b Institute of Automation, Information Engineering and Mathematics, Slovak University of Technology in Bratislava, Faculty of Chemical and Food Technology, Radlinského 9, 812 31 Bratislava, Slovakia

art ic l e i nf o

a b s t r a c t

Article history: Received 28 November 2014 Accepted 16 May 2015

Control of the air–fuel ratio in combustion engines is of imminent importance when aiming at reduction of the fuel consumption and mitigation of emissions. In this paper this is achieved by employing Model Predictive Control (MPC), which keeps the air-fuel ratio close to the stoichiometric value. The real-time implementation of MPC is accelerated by pre-computing the optimal solution and using the explicit MPC approach. The memory footprint of such a controller is subsequently reduced by employing a twolayered compression scheme. Quality of the regulation is demonstrated using real engine data and is compared to traditional control approach based on PID. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Gasoline engine Air/Fuel ratio control Identification Explicit Model Predictive Control

1. Introduction The problem of Air/Fuel Ratio (AFR) control plays an important role in the complex problem of emission reduction for combustion engines. The mixture quality is essential for the efficiency of a three-way catalytic converter and therefore suitable control techniques are needed to fulfill emission legislations. A comprehensive description of problems occurring in the field of AFR control, such as requirements on AFR performance as well as modeling of fuel and air path, is provided in Guzzella and Onder (2010). Besides basic issues and experimental results, the reference also describes various AFR control strategies, including feedforward controllers, conventional PI controllers with narrow-band λ sensor, H 1 feedback controllers, or an advanced feedback internal-model control approach utilizing narrow-band and wide-band λ sensors. The authors in Hendricks (1991) and Hendricks, Jensen, Kaidantzis, Rasmussen, and Vesterholm (1993) address the improvement of transient performance of conventional AFR control strategies. These strategies are based on combination of statically measured engine maps that ensure relatively fast AFR transient responses and PI controller used for proper operation in AFR steady-states. However, tuning of such controllers is cumbersome and boils down to a trial and error approach, which is time consuming.

n

Corresponding author. Tel.: þ 421 911 391 197; fax: þ421 5249 5315. E-mail address: [email protected] (M. Honek).

http://dx.doi.org/10.1016/j.conengprac.2015.05.009 0967-0661/& 2015 Elsevier Ltd. All rights reserved.

More sophisticated approaches are required to eliminate drawbacks stemming from manual tuning and to achieve better performance. Improvements can be found e.g. in Onder and Geering (1993), where the authors proposed to use an LQR controller. In Choi and Hendricks (1998) the authors introduced an observer-based fuel injection control algorithm to improve AFR operation, utilizing sliding mode control. Other approaches to AFR control include, but are not limited to, use of gain-scheduling (Alfieri, Amstutz, & Guzzella, 2009), fuzzy (Khiar et al., 2007), and multi-model (Zhang, Shen, & Marino, 2010) techniques. Even when such approaches represent a certain enhancement compared to the traditional PID-based strategies, they are in general not able to guarantee the best achievable performance while simultaneously guaranteeing satisfaction of all physical constraints. Recently, Model Predictive Control (MPC) has become a very popular optimization-based control strategy that is able to guarantee optimal performance while taking constraints into account. Several applications of MPC to control AFR are reported in the literature. For instance, Zhai, Yu, Tafreshi, and Al-Hamidi (2011, 2007) propose nonlinear MPC schemes for control of Air/Fuel Ratio for spark ignition combustion engines. A nonlinear MPC scheme for control of air management in the case of diesel engines is discussed in Garcia-Nieto, Martinez, Blasco, and Sanchis (2008). Although performing well on simulations, such a nonlinear approach is cumbersome from practical implementation point of view. In particular, complexity of the resulting nonlinear optimization prohibits its implementation on a typical electronic control unit (ECU). This concern is addressed, to some extent, in Wong, Wong, and Vong (2012), where MPC is implemented using online

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

sequential relevance vector machines. The algorithm was subsequently implemented as embedded Matlab code and executed via the LabVIEW interface and shown to outperform traditional control approaches. The main reason why MPC is not yet wide-spread in massproduced ECUs is mainly due to high complexity of closed-loop implementation of MPC. In particular, introducing feedback requires repetitively solving a given optimization problem at each sampling instant, inducing the necessity to use powerful computing units and complex, iterative-based algorithms. This wellknown drawback of MPC can be abolished by employing the concept of explicit MPC. Since pioneered by Bemporad, Morari, Dua, and Pistikopoulos (2002), the concept has garnered much popularity among theoreticians and practitioners because it allows MPC to be implemented effectively even on cheap hardware with low computational power. This is achieved by pre-computing, offline, the explicit solution to a given optimization problem. For a wide range of MPC setups, such an explicit solution takes a form of a Piecewise Affine (PWA) function which maps state measurements onto the optimal control inputs. The closed-loop implementation of explicit MPC hence boils down to a mere function evaluation. Not only is the function evaluation faster than solving an optimization problem, but it can also be performed in a division-free manner, i.e., only using additions and multiplications. Hence the implementation code is very simple and easy to certify even for mission-critical applications. Successful practical applications of explicit MPC range from power electronics, through mechanical and automotive systems, up to control of autonomous vehicles, see Kvasnica (2009) and references therein. The main drawback of explicit MPC, however, is that simplicity of implementation is traded for storage space. In particular, an explicit PWA feedback function often consumes hundreds of kilobytes to several megabytes of memory. Such an amount can easily exceed capabilities of a given hardware implementation platform, especially when the hardware has to accommodate multiple feedback loops. In order to implement an MPC algorithm for AFR control on a typical ECU, it is therefore of imminent importance to reduce the memory footprint to an acceptable level. The reader is referred to Kvasnica and Fikar (2012) and references therein for an overview of available methods which accomplish reduction of explicit MPC solutions. This work proposes an explicit MPC approach to AFR control in a spark ignition (SI) combustion engine that can be implemented even in ECUs with limited amount of memory storage. A low memory footprint of the controller is achieved by employing a data compression technique, preliminary reported in Szũcs, Kvasnica, and Fikar (2011). The compression is performed in two steps. First, redundant floating-point data is removed from the definition of the explicit feedback and replaced by integer pointers to the stack of non-redundant entries. Subsequently, these integer pointers are further compressed using Huffman encoding which replaces them by short bit-wise codes. As will be illustrated in Section 4.4, the proposed scheme reduces the memory consumption by a factor of 25 and hence enables the implementation of explicit MPC even in memory-constrained control platforms. A further advantage is that the compression is lossless and thus preserves control performance. The paper is organized as follows. First, Section 2 describes the physical setup and presents results of an experimental identification procedure which yields a mathematical model of the AFR loop suitable for control design. The MPC-based control synthesis is then presented in Section 3. Then, Section 4 describes the compression algorithm which reduces the memory footprint of the explicit MPC controller to a desirable level. Finally, Section 5 presents results of an experimental case study which compares the performance of the proposed MPC-based strategy against a

conventional setup feedforward part.

119

consisting

of

a

PID

controller

with

2. Physical setup and AFR modeling Combustion engines belong to a class of discrete event based systems. From digital AFR control point of view it would be logical to sample the signals in crankshaft domain instead of time domain and in this sense use the event based control strategy. Due to mixing dynamics in the exhaust manifold, mixture formation can be transformed from the discrete event process to continuous changes of AFR information. Therefore, classical controllers in time domain can be used. AFR in a SI engines is defined as a mass ratio of the air and fuel. Normalized AFR is defined as λ¼

ma 1 ; mf Lth

ð1Þ

where ma and mf represent, respectively, the mass of air and fuel. Moreover, Lth  14:64 is the theoretical mass of air necessary for the ideal combustion of a unit mass of fuel. In other words the value of Lth represents stoichiometric AFR. Particular values of ma and mf can be indirectly measured with a delay at the confluence point, shown in Fig. 1. To model the dynamical behavior of λ, the authors in Polóni, Johansen, and Rohal’-Ilkiv (2008) suggest to employ two subsystems: one for the air-path branch with the input variable tr representing the throttle position, and one for the fuel-path with the fuel pulse width fp as the manipulated variable, as proposed in Polóni et al. (2008). The engine speed ne together with the throttle position tr define the operating point ϕ. From the AFR control point of view, ne and tr represent measured disturbances, while fp is the control input which can be manipulated by the control algorithm. Depending on the value of λ, three possible scenarios can occur. In the first case with λ ¼ 1, a stoichiometric mixture is achieved. If the value of λ exceeds 1, than the mixture is lean, and on the contrary, in cases with λ o1, the mixture is rich. The control objective is to keep values of λ as close as possible to the stoichiometric level in transient and steady-state operating conditions as well. Transient operating conditions can be interpreted as changes in the engine speed ne or throttle position tr. Naturally, all these variables are kept at constant values in steady-state operating conditions. During transient, however, the mean value of inlet air flow into the cylinder is changing, therefore the controller is required to change the mass of injected fuel accordingly to keep λ  1. The dynamics of λ is influenced by the dynamics of air mass flow, fuel evaporation, and by the gas mixing effect. All these tr air filter injectors

air flow meter

f p ms dynamic brake connection ne min− 1 catalyst

λ− (confluence point)

Fig. 1. Engine setup with input/output relations; dashed arrows – inputs, solid arrows – outputs.

120

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

sample time and band values, were chosen to respect the engine specifications as elaborated in more details in Polóni et al. (2008). Standard identification algorithms (Söderström & Stoica, 1989) were then applied to estimate the parameters of a linear ARX model for the λ dynamics. In the chosen operating point with N observations of input-output data, the prediction error for the multi-input single-output autoregressive exogenous (MISO ARX) model is defined as

effects are nonlinear and depend on the current operating regime. To measure λ and to estimate its dynamical behavior in the proposed experimental workbench, the LSU 4.9 Bosch lambda probe was used. This probe is capable of measuring λ in the interval from 0.65 to 1. Furthermore, a multi-point injection system with two PS-DINJ 2/1 modules was used to drive four low impedance electromagnetic valves. Because two different injectors were driven by a single PS-DINJ 2/1 module, it was necessary to introduce an upper bound on the injection time. Existence of this upper bound stems from the fact that two injectors cannot be driven to open state at the same time, because both power outputs are controlled by one load current controller. For example, for the engine speed 2500 RPM, when two injectors driven by one module are ordered for cylinders for which the offset of strokes equals to 3601, the crankshaft rotates over this angle once every 24 ms. This value hence constitutes the upper bound on the fuel injection duration fp. Fuel is injected to the intake manifold, where constant pressure difference between the fuel pressure in rail and air pressure in intake manifold is maintained by a pneumatic fuel pressure controller on the level of 3 bars (1 bar ¼ 105 kg m-1 s-2). In this setup, the mass of injected fuel is directly proportional to injection time fp through the constant flow capacity determined by the constant pressure difference.

εðk; θÞ ¼ ΔλðkÞ  ΓT ðkÞθ;

ð2Þ

where


ΓðkÞ ¼  Δλðk  1Þ; …;  Δλðk  nyÞ; Δt r ðk 1  da Þ; …; Δt r ðk  nua da Þ; iT Δf p ðk 1  df Þ; …; Δf p ðk  nuf  df Þ ;

ð3aÞ

is the vector of past measurements, θ ¼ ½a1 ; …; any ; bt;1 þ da ; …; bt;nua þ da ; bf ;1 þ df ; …; bf ;nuf þ df T

is the vector of current model ARX parameters to be identified, and

2.1. Identification procedure

ΔλðkÞ ¼ λðkÞ  1;

ð5aÞ

Δf p ðkÞ ¼ f p ðkÞ  f p;e ;

ð5bÞ

Δt r ðkÞ ¼ t r ðkÞ  t r;e ;

ð5cÞ

are the deviation variables. Estimation of ARX parameters θ was achieved by minimizing the prediction error in a least squares sense, i.e.,

To derive a mathematical model of AFR, including the air and fuel paths in an engine, an ARX-based parameter estimation scheme was employed. Parameters of the model were estimated using real engine data as measured by the exhaust gas oxygen (EGO) sensor. An open loop identification experiment was carried out, exciting of both the air and the fuel subsystems simultaneously. The experiment always starts from the stoichiometric steady-state value in the given operating point ϕ, where t r;e and f p;e represent the throttle position and the fuel pulse width in this point. The speed of the engine was kept constant during the experiment. To excite the air and the fuel path dynamics, pseudorandom binary excitation signals (PRBS) were applied to the throttle and fuel pulse width signals. The time profile of experimentally measured values of λ is shown in Fig. 2, along with profiles of the throttle position tr and the fuel pulse width fp. The PRBS excitation signals were generated by the idinput function of the System Identification Toolbox for Matlab. Properties of the PRBS signals, such as the amplitude,

θ^ ¼ argmin θ

N X

ε2 ðk; θÞ:

ð6Þ

k¼1

Using the arx function of the System Identification Toolbox for Matlab, the following ARX AFR model was obtained: aðqÞΔλðkÞ ¼ bf ðqÞΔf p ðkÞ þ bt ðqÞΔt r ðkÞ þ eλ ðkÞ;

aðqÞ ¼ 1 0:423355q  1 þ 0:066642q  2  0:019692q  3 ; bf ðqÞ ¼  0:010509q  1 0:131001q  2 þ 0:005034q  3 ; bt ðqÞ ¼ 0:002580q  1 þ 0:029514q  2 þ 0:006737q  3 :

t r[ ◦ ] f p [s]

20

30

40

50

60

70

80

20

30

40

50

60

70

80

20

30

40

50

60

70

80

−3

x 10

4.38 4.13 10

λ [–]

1.15 1 0.85 10

ð8Þ

The corresponding model was obtained for the operating point ϕ, represented by the throttle position t r;e ¼ 241, fuel pulse width f p;e ¼ 4:38  10  3 s, and engine speed ne ¼ 2500 RPM, respectively.

24

4.63

ð7Þ

with

25

23 10

ð4Þ

time[s] Fig. 2. Time profile of experimentally measured values of λ.

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

Since the MPC approach proposed in Section 3 is based on a state-space model, the ARX representation in (7) was converted into xðk þ 1Þ ¼ AxðkÞ þ Bf Δf p ðkÞ þ Bt Δt r ðkÞ;

ð9aÞ

yðkÞ ¼ CxðkÞ;

ð9bÞ

with 2

0

6 A ¼ 4 0:250000 0 2

0:107789

0:157535

0

3

0

0:133285 7 5;

0:500000

0:423355

3

6 7 Bf ¼ 4 0:118056 5; 0:005159

2

 0:080540

ð10aÞ

3

6 7 Bt ¼ 4  0:524004 5;  0:021018

C ¼ ½0 0 0:5;

ð10bÞ

ð10cÞ

where Δf p ðkÞ and Δt r ðkÞ represent deviation variables defined per (5b) and (5c), respectively. The system's output y(k) represents ΔλðkÞ, which in turn denotes the deviation of λ from the desired stoichiometric value per (5a). With the aim to achieve closed loop behavior with zero steadystate error, state-space representation for MPC controller design was augmented with integral state variable q(k) to reject nonmodeled effects. Then the model takes the form: "

# #   "     xðk þ 1Þ xðkÞ A 0 Bf Bt ¼ þ Δf p ðkÞ þ Δt r ðkÞ; qðk þ 1Þ qðkÞ C 1 0 0 |fflffl{zfflffl}B~ t |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}zðk þ 1Þ |fflfflfflfflffl{zfflfflfflfflffl}A~ |fflfflfflffl{zfflfflfflffl}zðkÞ |fflffl{zfflffl}B~ f

ð11aÞ
 yðkÞ ¼ C 0 zðkÞ: |fflfflfflffl{zfflfflfflffl}C~

ð11bÞ

In what follows Δf p ðkÞ is considered as the manipulated variable, and Δt r ðkÞ is treated as a measured disturbance. 3. Explicit model Predictive Control of AFR The control objective is to manipulate the fuel pulse width fp such that λ  1 is maintained despite varying position of the ⋆ throttle tr. The optimal selection of the control input f p ðkÞ at the k-th discrete time step can be achieved by solving the following optimal control problem: min

N 1 X

(12e) is always feasible. The manipulated variables u(k) are selected such that the actual fuel pulse width, given by f p ðkÞ ¼ uðkÞ þ f p;e , respects physical constraints per (12e). The cost function (12a) employs three tuning parameters: Q z , Q u , and Q ϵ , which respectively penalize the augmented states, increments of the control inputs, and the slack variables ϵ. The trade-off between achieving λ  1 and obtaining a smooth profile of the manipulated variable can be thus adjusted by changing the Q ϵ =Q u ratio with higher values leading to a better tracking and lower values corresponding to a more smooth control profile. Problem (12a) is initialized by measurements (or estimates) of xðkÞ, q(k), and Δt r ðkÞ, which represents the measured disturbance. It is assumed that its value is constant throughout the prediction. For given measurements of the extended states zðkÞ and the ⋆ disturbance Δt r ðkÞ, the optimal fuel pulse width f p ðkÞ can be obtained by solving (12a) as a convex quadratic program (QP). Its solution yields the open-loop optimal sequence of control actions fu⋆ ðkÞ; …; u⋆ ðk þN  1Þg, from which only the first component, i.e., ⋆ f p ðkÞ ¼ u⋆ ðkÞ þ f p;e is employed for feedback control according to the receding horizon principle. However, to minimize the increments of control inputs Δu in (12a), it is also required to provide the value of uðk  1Þ, i.e., the control action from the previous time instant. Therefore the QP (12a) is parameterized in zðkÞ (which consists of xðkÞ and q(k) per (11a)), Δt r ðkÞ, and uðk  1Þ. In what follows these parameters are aggregated into a single vector of parameters ξ A R6 : 2 3 xðkÞ 6 7 6 qðkÞ 7 7 ð13Þ ξðkÞ ¼ 6 6 Δt r ðkÞ 7: 4 5 uðk  1Þ To avoid having to re-solve the QP (12a) every time the measurements aggregated in ξðkÞ change, it is proposed to precompute, off-line, the explicit solution to (12a) for all feasible values of the vector of parameters ξ using parametric programming (Bemporad et al., 2002). By doing so an explicit relation between ⋆ the parameters ξ and the optimal control inputs f p ¼ u⋆ þ f p;e in the form of a piecewise affine (PWA) function is obtained. Theorem 3.1 (Bemporad et al., 2002). The receding horizon feedn back f p ¼ κðξÞ for the problem in (12a) is given by 8 > < F1 ξ þ g 1 if ξ A R1 n ⋮ ð14Þ f p ¼ κðξÞ≔ > : F ξþg if ξ A R ; R

zðk þt þ 1ÞT Q z zðk þ t þ 1Þ þ ⋯ 2

⋯ þQ u ðΔuðk þ tÞÞ þ Q ϵ ϵðk þ tÞ

R

R

where

t¼0 2

121

ð12aÞ

s.t. ~ zðk þ t þ 1Þ ¼ Azðk þtÞ þ B~ f uðk þ tÞ þ B~ t Δt r ðkÞ;

ð12bÞ

~ δ  εðk þtÞ r Czðk þtÞ r δ þ εðk þ tÞ;

ð12cÞ

εðk þ tÞ Z0;

ð12dÞ

0 r uðk þ tÞ þf p;e r 24 ms;

ð12eÞ

Δuðk þ tÞ ¼ uðk þ tÞ  uðk þt  1Þ;

ð12f Þ

where constraints in (12b)–(12f) are imposed for all t ¼ 0; …; N  1 with N denoting the prediction horizon. Predictions of the extended state vector z are updated according to the linear dynamics (11a) where uðkÞ  Δf p ðkÞ. Therefore, f p ðkÞ ¼ uðkÞ þ f p;e by (5b). Moreover, (12c) forces λ to stay in the 7 δ neighborhood of the desired value of λ ¼ 1 via (11b) and (5a). This neighborhood is softened using non-negative slack variables ϵ to guarantee that

Ri ¼ fξj Hi ξ r Ki g;

i ¼ 1; …; R;

ð15Þ

are polytopic critical regions and R denotes the total number of critical regions. The advantage of explicit MPC feedbacks in (14) is that they provide a simple and fast implementation of MPC. In particular, for n given measurements of ξ, computing the optimal control action f p reduces to a mere evaluation of the PWA function κðÞ. Even for a n large number of critical regions, computing f p by evaluating κðξÞ via (14) can be performed much faster compared to solving (12a) at each sampling instant. The main drawback, however, is that the amount of memory consumed by the explicit feedback law κðÞ in (14) due to storing the critical regions Ri and the corresponding matrices Fi , gi is often large and can reach several hundreds of kilobytes. Thus it can exceed capabilities of the implementation hardware, especially when the control platform has to accommodate multiple control loops and/or when not the whole storage space is reserved just for a single controller. Therefore the next section shows how to

122

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

alleviate this limitation by reducing the memory footprint of the explicit MPC feedback law (14) by using data compression techniques.

4. Data compression technique for explicit MPC Memory footprint of the explicit MPC feedback law κðÞ is determined by two factors: by the number of polytopes R, and by the number of constraints defining individual polytopes Ri , cf. (15). The proposed data compression scheme applies two layers of compression to decrease the number of constraints of each polytope, hence reducing the total memory footprint of κðÞ. 4.1. Data de-duplication layer This layer reduces memory footprint of κðÞ in (14) by eliminating redundant entries from the definition of regions Ri in (15). First, introduce two sets of stacks 2 3 2 3 H1 K1 6 ⋮ 7 6 ⋮ 7 H¼4 ð16Þ 5 5; K ¼ 4 HR KR ; and let I dH denote the indices of rows of H which are duplicate under multiplication by 1 or by  1, i.e., I dH ¼ fij (j a i s:t: ½Hi ¼ ½Hj or ½Hi ¼  ½Hj g;

ð17Þ I uH

denote the where ½i represents the i-th row of the matrix. Let index set of unique rows, i.e., I uH ¼ f1; …; Mg⧹I dH , where M is the total number of rows of H in (16). Index set I uK , is constructed in a similar fashion. Note that, in practice, it is recommended to use a small numerical tolerance to detect unique rows. In this manuscript, all entries in (16) are first rounded to 5 decimal places before identifying unique rows. If the cardinality of I uH is smaller than the number of rows of H (which always appears to be true in practical cases), then the amount of memory required to represent the explicit feedback κðÞ in (14) can be decreased as follows. Let Hu ¼ f½Hi j i A I uH g be the unique rows of H under multiplication by 71, and denote by Ku the unique rows of K. Assume these sets of unique data are stored in the memory. Then, instead of storing matrices Hi and Ki in full, storing the i-th polytope Ri from (15) only requires recording the indices to the stacks Hu and Ku . The advantage being that such indices are signed integers, and therefore consume less memory compared to storing full floating point numbers. Typically, integers require 2 bytes of memory, while floating point numbers need to be represented either by 4 or 8 bytes, depending on the chosen precision. As an example, consider three regions given in their respective half-space representation by 2 3 2 3 2 3 0 1 1 0:5 0 1 6 1 6 1 6 0  0:5 7 1 7 0:5 7 6 7 6 7 6 7 H1 ¼ 6 7; H 2 ¼ 6 7; H 3 ¼ 6 7; 4 0 4 1 4 0 1 5 1 5 0:5 5 1 2

0:5 3

2:4 6 3:1 7 6 7 K1 ¼ 6 7; 4  1:5 5 0

0 3

2

2:4 6 5:0 7 6 7 K2 ¼ 6 7; 4  3:1 5 0

1 2

1

 0:5

3

0 6 0 7 6 7 K3 ¼ 6 7: 4 1:5 5 3:1

The stacks of rows that are unique under the 7 1 scaling are then given by Hu ¼ f½0 1; ½1  0:5g; Ku ¼ f  3:1;  1:5; 0; 2:4; 5:0g: Then the (signed) index set representation of these three regions

becomes I H1 ¼ f1; 2;  1;  2g; I K 1 ¼ f4;  1; 2; 3g; I H2 ¼ f1; 2;  2;  1g; I K 2 ¼ f4; 5; 1; 3g; I H3 ¼ f 2;  1; 1; 2g; I K 3 ¼ f3; 3;  2;  1g;

ð18Þ

where each element of I Hi and I K i points to the corresponding entry in Hu and Ku and hence allows to re-construct the i-th polytope on-the-fly. In this simple example, memory footprint of regions Ri was reduced from 36 floating point numbers representing matrices Hi , Ki to 9 floats for Hu , Ku , and 24 integer pointers. Assuming that one float is represented by 8 bytes and one integer by 2 bytes, the proposed de-duplication layer reduces required memory from 288 bytes to 120 bytes. In practical cases, though, the achievable reduction is even higher, as will be documented in Section 5.2. Needless to say, the same de-duplication procedure can also be applied to reduce memory footprint of the local feedback laws Fi ξ þ g i in (14). 4.2. Huffman encoding layer In this layer, which is applied on top of the de-duplication procedure of the previous section, the signed-integer representation is further simplified by devising shorter bit-wise codes to represent most abundant integers. Given are index set representations I H ¼ [ i I Hi and I K ¼ [ i I K i , whose entries point to corresponding rows in the set of unique elements Hu and Ku . In traditional implementation, each element of I H and I K would need to be represented as a (signed) integer, i. e. by 16 bits (provided that cardinality of H and K does not exceed 216). A more efficient representation can be obtained by using a prefix-free variable-length encoding (Knuth, 1985) where bit-wise codewords are assigned to each element of the index sets. Length of a codeword is inverse-proportional to its abundance, allowing to compress I K and I H as much as possible. As an illustration, consider the index sets in (18) and let I K ¼ I K 1 [ I K2 [ I K 3 . Then the integers to encode appear with frequencies reported in Table 1. The corresponding Huffman tree is shown in Fig. 3. Here, the optimal codewords are CðI K Þ ¼ f½  2 : 111; ½  1 : 001; ½1 : 110; ½2 : 101; ½3 : 01; ½4 : 000; ½5 : 100g. It is easy to verify that such an encoding is prefix-free, i.e., that no codeword is a prefix of another codeword. Moreover, the most abundant integer 3 is encoded using fewest number of bits as to minimize the total length of binary representation of I K . Hence, instead of storing I K as an array of 12 integers (i.e., 24 bytes), it suffices to store the tree (7 integers or 14 bytes) and 12 codewords of a total size 32 bits, or 4 bytes. Constructing the Huffman tree from a given set of integer indices can be performed using priority queue, as shown in Dasgupta, Papadimitriou, and Vazirani (2006). The tree can be represented as a linked list of nodes that represent the unique integer pointers. Decompression of codewords is performed by traversing the tree, starting from the root node. Depending on the value of the first bit of a particular codeword, the tree is either traversed to the left (if bit equals to 0), or to the right. The processed bit is then discarded and the procedure repeated with the next bit until a leaf node is reached. Value of this node provides the decoded integer value of the binary codeword, which is then used to extract the corresponding row from the stacks of unique data Hu and Ku . Note that both of these stacks have their separate decoding tree associated to them. Table 1 Frequencies of integers to encode. Integer Frequency

2 1

1 2

1 1

2 1

3 4

4 2

5 1

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

0

1

3 4

-1

5

2

1

-2

Fig. 3. Huffman tree for I K ¼ f4;  1; 2; 3; 4; 5; 1; 3; 3; 3;  2;  1g.

4.3. Decompression To reiterate, the de-duplication layer of Section 4.1 replaces individual expressions of matrices Hi and Ki in (15) by integer pointers I Hi and I K i to the stacks Hu and Ku of unique rows. Subsequently, the Huffman encoding layer of Section 4.2 replaces the integer pointers by binary codewords CðI Hi Þ and CðI K i Þ of minimal length such that the prefix-free property holds. Finally, only the stacks Hu , Ku are stored, together with the codewords CðI Hi Þ and CðI K i Þ for each critical region, along with the associated Huffman trees. ⋆ To obtain the value of f p for a particular value of ξ from (14) at each sampling instant, it is necessary to identify which critical region contains the vector ξ. This can be performed by searching through all critical regions sequentially as follows: 1. For the i-th critical region, convert the codewords CðI Hi Þ and CðI K i Þ into the corresponding integer pointers I Hi and I K i by traversing the corresponding Huffman trees top to bottom. 2. Recover Hi and Ki from the rows of Hu and Ku indexed by I Hi and I K i . 3. Check whether Hi ξ r Ki . If all inequalities are satisfied, set i⋆ ¼ i and abort. 4. Otherwise set i ¼ i þ1 and continue from step 1. Upon termination, the algorithm yields the index i⋆ of the critical region that is “active” for the particular value of the vector of parameters. Then the optimal fuel injection pulse is given by ⋆ f p ¼ F i⋆ ξ þ g i⋆ . Notice that such an implementation algorithm is very simple, as it only performs multiplications, additions, and comparisons. Moreover, it can be straightforwardly implemented in the C-language using basic data structures such as linked lists. Finally, it is important to note that the control action obtained via decompression is identical to that of the uncompressed controller. This follows from the fact that de-duplication and Huffman encoding layers merely encode the feedback law in (14) in a more memory-efficient way. 4.4. Compression of the AFR MPC controller The explicit representation of the optimal feedback law in (14) was first computed by formulating (12a) using YALMIP (Löfberg, 2004) and solving it according to Theorem 3.1 by the MultiParametric toolbox (Herceg, Kvasnica, Jones, & Morari, 2013). With N ¼25, δ ¼ 0:01, Q z ¼ diagð0; 0; 1; 1Þ, Qu ¼0.05, and Q ϵ ¼ 1  104 , the PWA function in (14) consists of 735 regions in the 6dimensional parametric space, cf. (13). The total amount of space consumed by the controller's critical regions Ri was 63,196 floating point numbers (which corresponds to a total of 9028 rows in the stacks H and K in (16)). When each floating point number is represented by a 8 bytes according to the doubleprecision standard, the controller's regions would consume 493 kB ⋆ in the ECU. The evaluation of f p for a specific value of ξ per (14) per the procedure of Section 4.3 requires, at most, 117,364 floating point operations (FLOPs) using the sequential search approach, i.e., when Steps 1 and 2 of the implementation algorithm are omitted.

123

Since only a portion of the ECU storage capacity is reserved for the air/fuel controller, the controller's memory footprint was subsequently reduced by using the compression techniques described above. First, the stacks H and K were formed per (16). As reported above, initially, all controller regions are defined using 9028 inequalities. However, 8641 rows in H and 8565 entries in K are redundant under (17). Therefore the unique stacks Hu and Ku are represented by 387 and 463 rows, respectively. The duplicate entries required to recover each critical region Ri in (15) are then recorded using integer pointers into Hu and Ku . In total, storing the unique stacks requires 11 kB of memory, and the integer pointers contribute with additional 35 kB. Therefore applying the de-duplication compression layer allows to decrease the memory footprint of the critical regions in (14) from 493 kB to just 46 kB, more than a tenfold reduction. However, when fp is evaluated from (14) compressed using de-duplication, the critical regions need to be de-compressed, one by one, on-the-fly to check which critical regions contain ξ. The additional decompression effort, which corresponds to Step 2 of the algorithm in Section 4.3, is just 850 FLOPs, an increase by mere 0.7% over the runtime figure reported above. The integer pointers obtained during de-duplication can be further compressed using the Huffman encoding layer. Specifically, designing prefix-free bit-wise codes for the most abundant integer pointers allowed to further decrease the memory required to store such pointers from 35 kB to just 9 kB. Therefore the total footprint of the controller is just 20 kB, a reduction by a factor of 2.3 versus de-duplication, or by a factor of 25 compared to using the “raw”, uncompressed solution. The price one has to pay is that the decompression of Huffman codes induces additional computational effort. Specifically, it requires additional 26,291 FLOPs to be performed at Step 1 of the algorithm of Section 4.3 at each sampling ⋆ instant when f p is evaluated from (14) for a particular value of ξ. Aggregated results for individual layers of the proposed compression technique are reported in Table 2. It is important to note that both compression layers preserve optimality of the compressed controller, i.e., no amount of control performance is lost due to compression/decompression. From a computational point of view, the construction of the explicit feedback law in (14) took 2 min on a 2.7 GHz machine using MPT3. The subsequent compression took under one second. Runtimes of these magnitudes allow for re-tuning iterations to be performed in a reasonable time.

5. Case study. This section compares the performance of the proposed explicit MPC controller, compressed using the algorithms of Section 4 against a conventional control strategy and discusses implementation details. 5.1. Hardware setup All experiments presented in this paper, including the identification procedure of Section 2.1, were performed on a spark ignition engine that is mounted in VW Polo and Škoda Fabia cars, Table 2 Efficiency of the compression. Compression level

Memory footprint (kB) Implementation effort (kFLOPs)

Original solution (14) 493 De-duplication 46 Huffman 20

117 118 144

124

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

marked with the product code AUA. The picture of the engine is shown in Fig. 4 and its parameters are in Table 3, reported in the Appendix. All measurements and control algorithms used dSpace (dSPACE, 2009) Rapid Control Prototyping (RCP) platform, which was connected to the engine. The Master/Slave (dSpace/RapidPro) configuration depicted in Fig. 5 was employed.

Fig. 4. The spark ignition engine on the experimental workbench.

Table 3 Engine specifications. Car Code Cubature Number of cylinders Number of inlet valves Borehole  throw Compression ratio Power Torque Jetting Ignition Fuel

Škoda Fabia or VW Polo AUA 1390 cm3 4 16 76:5  75:6 mm 10.5 55 kW (75 hp) at 5000 RPM 126 Nm at 3800 RPM Multi point, 1 injection valve/cylinder Electronic Natural 95

The modular master system is based on the microprocessor board DS1005, which serves as the main processing and host interface unit. The DS1005-09 revision of the board is equipped with the IBM PowerPC 750GX processor which runs at 1 GHz and provides 128 MB of memory storage. Host service specific code as a part of overall real-time application provides interface for user handling, data acquisition, and monitoring of real time application with the use of a host computer. Physical layer of communication between the host computer and the master system is made by the connection between host interface module DS814 on the side of master system and the ExpressCard-based module DS821 on the side of the host computer. The slave modular system is based on the DS1602 microcontroller module, equipped with the Motorola MPC565 microcontroller running at 56 MHz. The slave system consists of two units: the control unit and the power unit. The DS1602 microcontroller module provides the control functionality along with signal conditioning and communication modules and power modules for driving the actuators, as summarized in Table 4. The real-time control algorithm is executed on the master system in the Master/Slave arrangement. Control commands generated by the algorithm are then communicated to the slave system. The slave system performs signal conditioning and provides to the master system measurements of the crankshaft angle, engine speed, λ, temperature of the lambda probe's heating element, and information about the knocking sensor amplitude at user-specified frequencies. The bi-directional information exchange between the master and the slave system is achieved by high-speed LVDS (Low Voltage Differential Signaling) connection with the rate of 250 Mbit/s, with the use of communication modules DS1606 and DS4121, see Tables 4 and 5. To measure other signals and to drive the power electronics for electronic throttle positioning and load torque control, the DS2202 I/O interface board was used as a module in master system. The technical specifications are provided in Table 5. The dSpace RCP system is fully programable from Simulink. This allows the MPC-based control algorithm to be programmed in the C language and be provided to dSpace in form of a C mex Sfunction. However, besides the MPC algorithm, the control software performs other tasks as well, such as monitoring of process variables and communication. Therefore just a portion of the overall dSpace capacity is available to accommodate the MPCbased control algorithm.

RapidPro Power Unit

Master

Actuators RapidPro SC Unit

Unit Connection Bus (UCB)

RapidPro Control Unit

Unit Connection Bus (UCB)

RCP System

Sensors HostPC with Configuration Desk

Slave

LVDSL ink USB for configuration

dSPACE Link Board (DS81x) or slot CPU Fig. 5. The dSpace rapid control prototyping setup.

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

125

Table 4 List of RapidPro modules (RapidPro as a Slave system). Type

Description

MC-MPC565 1/1 module (DS1602) SC-CCDI 6/1 module (DS1637) SC-EGOS 2/1 module (DS1634)

Microcontroller module, based on the MPC565 microcontroller Conditioning of crankshaft and camshaft signals Conditioning of exhaust gas oxygen sensor signal (lambda probe LSU4.9 Bosch) Conditioning of knock sensor signal Configuration from Host PC Communication with Master dSpace system Low side driver, for driving of lambda probe heater element and ignition modules Driving of electromagnetic injection valves

SC-KNOCK 4/1 module (DS1635) COM-USB-CI 1/1 module (DS1609) COM-LVDS 1/1 module (DS1606) PS-LSD 6/1 (DS1662) 2  PS-DINJ 2/1 module (DS1664)

Table 5 List of dSpace modules. Type

Description

DS1005 module

Processor board, based on real-time processor PowerPC 750GX, it is a main processing unit and host interface I/O board, analog inputs and outputs for: control of electronic throttle position, MAF sensor, MAP sensor, oil and cooling water temperatures, braking torque and battery voltage Interface board, makes communication link between Master and Slave systems Host computer communication interface

DS2202 module

DS4121 module DS814 module

To enable assessment of performance of different control algorithms, the stock SI engine described previously was modified as follows. First, the Bosch HFM5 air mass flow sensor was installed in order to directly measure the air mass flow upstream of the throttle. Next, the 036133035A pneumatic fuel pressure regulator was added to maintain constant difference between the intake manifold and fuel pressures. This ensures that the fuel flow rate is constant. Moreover, to enable a continuous measurement of λ, the wideband Bosch LSU 4.9 λ probe was installed to replace the narrow-band exhaust gas oxygen sensor which only provides binary information. To be able to inject higher mass of fuel, the existing fuel injectors were replaced by more capable ones. Finally, since the SC-CCDI 6/1 module for crankshaft signal conditioning requires a zero crossing signal, a device based on an analog comparator was inserted between the hall sensor and the SC-CCDI 6/1 module. All these modifications allow to control the engine in the whole range of operating regimes with a high level of accuracy. 5.2. Experimental results To assess the quality of the proposed MPC-based control of the fuel injection times, two different experiments were performed. First, to obtain a relevant basis for comparison, AFR was controlled by a conventional PID controller with a feed-forward part. Obtained experimental profiles are shown in Fig. 6. The engine speed was kept constant at ne ¼ 2500 RPM during the whole experiment. The feed-forward part of this controller computes the injection time fp proportionally to the air mass into the cylinder per one suction stroke, estimated by a numerical integration of the air mass flow signal. The discrete-time PID controller then serves to correct the actions of the feed-forward controller to reject offsets due to uncertainties and external unmeasured disturbances. It is worth emphasizing that the controller operates in the crankshaft angle domain. The air mass flow signal is therefore sampled each 3:61 of the crankshaft angle and the λ signal is sampled each 1801. Therefore the control commands are computed and applied to the fuel injector valves each 1801 of the crankshaft angle. In the second experiment the fuel injection times fp were controlled by the explicit MPC controller (14) whose memory

footprint was compressed using the approaches described in Section 4. For particular measurements of the vector of parameters ⋆ ξ (cf. (13)), the controller's code computes f p by evaluating the PWA function (14) per the algorithm reported in Section 4.3. The algorithm was implemented in the C language and exported to dSpace via the Real-Time Workshop. The performance of the MPC-based control policy is shown in Fig. 7. Again, the engine speed was kept constant at ne ¼ 2500 RPM and the same profile of tr as in Fig. 6 was considered. As can be seen, the MPC controller manages to keep λ within the 7 1% zone around the desired stoichiometric value. However, unlike the PID controller in Fig. 6, the MPC-based policy provides a significantly smoother profile of the injection times. This is a consequence of minimizing increments of the control action in (12a). Note that the constraints on the λ zone are often active in the middle plot in Fig. 7, cf. (12c). The main difference, however, is in the way the two control strategies achieve the required tracking of λ  1. In the PID case, good performance was achieved by a proper tuning, which is a trial-and-error and time-consuming procedure. The MPC strategy, on the other hand, directly enforces λ A ½1  δ; 1 þ δ via (12c). Therefore the MPC policy represents a much more systematic approach to achieving desired performance and safety characteristics. Moreover, the MPC-based controller does not require information about the air mass flow, hence it does not rely on the intake manifold pressure or the air mass flow sensors. However, they typically still need to be present for engine diagnostics. Finally, even though the MPC controller runs at sampling time of 0.1 s, which is about 8 times slower than the PID controller, it provides at least comparable, if not better, performance. Less frequent execution of the algorithm allows the ECU to perform other tasks which otherwise would be suppressed by a more frequent execution of the conventional control algorithm.

6. Conclusion This paper has presented the design, implementation, and experimental validation of a memory-efficient explicit Model

126

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

tr[◦ ]

25 24 23

0

10

20

30

40

λ [–]

1.05

λ

1.01

λ

0.99

1 0.95

4.9091

f p [s]

50

0

10

20

10

20

30

40

50

30

40

50

−3

x 10

4.5091 4.1091

0

time [s] Fig. 6. Air/Fuel Ratio controlled by the feed-forward þ PID controllers.

tr [ ◦ ]

25 24 23

0

10

20

30

40

λ [–]

1.05

50

λ

1.01

λ

0.99

1 0.95

0

10

20

10

20

30

40

50

30

40

50

−3

f p [s]

5.0049

x 10

4.6049 4.2049

0

time[s] Fig. 7. Air/Fuel Ratio controlled by the explicit MPC controller.

Predictive Control strategy for controlling the Air/Fuel Ratio in spark ignition combustion engines. First, the dynamical model of the Air/Fuel Ratio was obtained by system identification techniques from real engine data. Then an MPC optimization problem was set up to force λ  1 by using soft constraints. The MPC problem was solved using parametric programming, which gave rise to an explicit relation between measured states and the optimal control inputs. However, although such an explicit solution features a fast and simple implementation, it often requires too much storage space, which might not be available in standard ECUs (especially when they have to accommodate multiple feedback loops). Therefore the paper has proposed a compression strategy to reduce the required memory footprint. The performance of the suggested MPC strategy was compared against a conventionally used control approach in an experimental study. The results show that the MPC controller achieves desired performance while providing a much more smooth profile of control actions. The main advantage of the presented methodology is that it allows to include performance specifications and constraints into control design in a systematic fashion. Specifically, the satisfaction of the λ zone is achieved via constraints, which means it is independent of tuning parameters. This allows for more rapid re-tuning iterations since change of the tuning parameters does not compromise constraint satisfaction.

Acknowledgments The work has been supported by the Slovak Research and Development Agency under grant APVV 0551-11, APVV 0090-10 and APVV 0015-12. This support is very gratefully acknowledged.

References Alfieri, E., Amstutz, A., & Guzzella, L. (2009). Gain-scheduled model-based feedback control of the air/fuel ratio in diesel engines. Control Engineering Practice, 17, 1417–1425. Bemporad, A., Morari, M., Dua, V., & Pistikopoulos, E. N. (2002a). The explicit linear quadratic regulator for constrained systems. Automatica, 38, 3–20. Choi, S. B., & Hendricks, J. K. (1998). An observer-based controller design method for improving air–fuel characteristics of spark ignition engines. IEEE Transactions on Control Systems Technology, 6, 325–334. Dasgupta, S., Papadimitriou, C., & Vazirani, U. (2006). Algorithms (1st edition). New York, NY, Columbus: McGraw-Hill Science/Engineering/Math. dSPACE GmbH. HelpDesk Application, 2009. Garcia-Nieto, S., Martinez, M., Blasco, X., & Sanchis, J. (2008). Nonlinear predictive control based on local model networks for air management in diesel engines. Control Engineering Practice, 16, 1399–1413. Guzzella, L., & Onder, C. H. (2010). Introduction to modeling and control of internal combustion engine systems. Berlin, Heidelberg: Springer-Verlag. Hendricks, E. (1991). Transient errors in classical SI engine controllers. In Proceedings of the 1991 American control conference (ACC), paper no. TP14-15:00.

M. Honek et al. / Control Engineering Practice 42 (2015) 118–127

Hendricks, E., Jensen, M., Kaidantzis, P., Rasmussen, P., & Vesterholm, T. (1993). Transient A/F ratio errors in conventional SI engine controllers. SAE technical paper no. 930856. Herceg, M., Kvasnica, M., Jones, C.,& Morari, M. (2013). Multi-parametric toolbox 3.0. In 2013 European Control Conference (pp. 502–510). Khiar, D., Lauber, J., Floquet, T., Colin, G., Guerra, T., & Chamaillard, Y. (2007). Robust Takagi–Sugeno fuzzy control of a spark ignition engine. Control Engineering Practice, 15, 1446–1456. Knuth, D. E. (1985). Dynamic Huffman coding. Journal of Algorithms, 6, 163–180. Kvasnica, M. (2009). Real-time model predictive control via multi-parametric programming: Theory and tools. Saarbruecken: VDM Verlag. Kvasnica, M., & Fikar, M. (2012). Clipping-based complexity reduction in explicit MPC. IEEE Transactions on Automatic Control, 57, 1878–1883. Löfberg, J. (2004). YALMIP. Available from: 〈http://users.isy.liu.se/johanl/yalmip/〉. Onder, C. H., & Geering, H. P. (1993). Model-based multivariable speed and air-to-fuel ratio control of an SI engine. SAE technical paper no. 930859. Polóni, T., Johansen, T. A., & Rohal’-Ilkiv, B. (2008). Modeling of air–fuel ratio dynamics of gasoline combustion engine with arx network. Journal of Dynamic

127

Systems, Measurement and Control-Transactions of the ASME, 130, 061009/ 1–061009/10. Söderström, T., & Stoica, P. (1989). System identification. New York: Prentice Hall. Szũcs, A., Kvasnica, M., & Fikar, M. (2011). A memory-efficient representation of explicit MPC solutions. In Proceedings of the 50th CDC and ECC, Orlando, Florida (pp. 1916–1921). Wong, H. C., Wong, P. K., & Vong, C. M. (2012). Model predictive engine air-ratio control using online sequential relevance vector machine. Journal of Control Science and Engineering, 15. Zhai, Y. J., Yu, D. L., & Wang, K. L. (2007). Comparison of single-dimensional and multi-dimensional optimization approaches in adaptive model predictive control for air–fuel ratio of SI engines. International Journal of Information and Systems Sciences, 3, 129–149. Zhai, Y. J., Yu, D. L., Tafreshi, R., & Al-Hamidi, Y. (2011). Fast predictive control for air–fuel ratio of SI engines using a nonlinear internal model. International Journal of Engineering, Science and Technology, 3, 1–17. Zhang, J., Shen, T., & Marino, R. (2010). Model-based cold-start speed control scheme for spark ignition engines. Control Engineering Practice, 18, 1285–1294.