Cruise control using model predictive control with constraints

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/compag

Cruise control using model predictive control with constraints T. Coen ∗ , J. Anthonis, J. De Baerdemaeker BIOSYST – MeBioS, Faculty of Bioscience Engineering, K.U. Leuven, Kasteelpark Arenberg 30, B-3001 Leuven, Belgium

a r t i c l e

i n f o

a b s t r a c t

Article history:

Model Predictive Control (MPC) originated in process industry, but in recent years it has been

Received 1 October 2007

used in many applications beyond this sector. MPC implies solving a quadratic optimisation

Received in revised form

problem with constraints online. For a practical prototype implementation, it is often neces-

17 February 2008

sary to implement the optimisation method yourself. First, the application is presented; then

Accepted 5 March 2008

an overview of the most important optimisation methodologies is given. A modified Active Set Method (ASM) is derived. The advantage of this method is its straightforward implemen-

Keywords:

tation and good timing behaviour. This method is compared to a commercially implemented

Cruise control

Interior Point Method (IPM) on a numerical example and on a real-life implementation.

Optimisation problems

© 2008 Elsevier B.V. All rights reserved.

Model-based Predictive Control Quadratic programming

1.

Introduction

The objective is to design a cruise control system for a combine harvester. The propulsion system of modern combine harvester consists of a diesel engine (with variable engine speed), which drives a hydrostatic pump (with variable pump setting). The machine speed thus depends on the rotation speed of the diesel engine and on the pump setting. The relation between both input signals and the actual speed is nonlinear, which can be translated into a time-varying linear model. The input signals have to be restricted to guarantee operator comfort and safe machine operation. The time-varying nature of the model, as well as the need for constraints, call for an online calculated control law. Model Predictive Control (MPC) is very well suited for this task. MPC was launched in the process industry (Wang et al., 1997), but in recent years MPC is being used in all kinds of sectors (Blasco et al., 2007; Bellemans et al., 2006; Pinon et al., 2005). MPC implies that an optimisation problem is solved online. In most cases the optimisation problem consists of a

∗

Corresponding author. Tel.: +32 16 32 16 13; fax: +32 16 32 85 90. E-mail address: [email protected] (T. Coen). 0168-1699/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.compag.2008.03.003

quadratic objective function and linear equality and inequality constraints. Typical of the earlier applications in process industry, is the slow update rate of the controller. An update rate of once a minute or even less is no exception. Applications in other sectors often are slow processes or simulations, which do not call for a fast calculation either. The cruise control system of a harvester however calls for a sampling rate of 20 Hz. Of course, a lot of methodologies exist to solve these optimisation problems very fast. However, when building a prototype, or a commercial implementation with a low volume, not only computational efficiency (number of operations necessary to solve the problem) but also the time required to implement the method is important. After all, often these methods are not yet available on the implementation platform of choice. For instance, LabView (National Instruments, Austin, TX, USA), a leading software platform for real-time prototype implementation, only offers an Active Set Method (ASM) since mid 2006 (an Interior Point Method (IPM) is already long before available).

228

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

Nomenclature S PS RPM ˇ(PS) xk yk PSk RPMk S uk

Q S R X Y U

travel speed (km h−1 ) pump setting (%) engine speed of the diesel engine (rpm) local amplification factor between PS and S state vector of the system model at time step (k) output vector of the system at time step k (travel speed, km h−1 ) change of pump setting at time step k (%) change of engine speed at time step k (rpm) Change of travel speed (km h−1 ) input vector, containing the change of pump setting and the change of engine speed at time step (k) weighting matrix of the state error weighting matrix of the output error weighting matrix of the input feasible region for xk feasible region for yk feasible region for uk

Computational efficiency determines which hardware is needed for the application. The more computing power is needed, the higher the cost. On the other hand, for low volume applications, the engineering cost to implement the method is also significant. For a lot of applications, the fastest, most complicated method thus is not the most economical solution. In the first part of this paper the cruise control application is presented. In the next section some background to quadratic programming is presented. Subsequently a variation on the general Active Set Method is derived, and contrasted with the standard implementation. This method is benchmarked with an Interior Point Method on a problem from literature. Then the derived method is extended to solve linear programs (LPs) as well. Finally the results of the method on the implementation of a cruise control system on a combine harvester are shown.

2. Application: cruise control for a combine harvester 2.1.

System description

The propulsion system of a combine harvester consists of a diesel engine, which drives a variable hydrostatic pump (Fig. 1). The hydrostatic flow thus depends on the rotation speed of the diesel engine and on the pump setting. The hydrostatic flow as well as the engine speed can be set independently by the operator. The engine speed can vary between a minimum and a maximum level. For the combine harvester used for the field tests in this paper (a New Holland CR) this is respectively 1300 and 2100 rpm. Of course the exact values depend on the combine type. Note however, that different combinations of engine speed and pump setting result in the same machine speed. During road transport, the operator can

Fig. 1 – A scheme of the propulsion system of a combine harvester.

freely vary the engine speed and pump setting. To reduce the noise caused by high engine speed and the exhaust emissions of the machine, the engine speed should be at all time as low as possible. This is necessary to comply with the latest legislation (EC, 2004, 2005, 2006). This cruise control system is implemented using Model Predictive Control. MPC controllers are designed based on a dynamical model of the system that has to be controlled (i.e. the plant), and apply mathematical optimisation techniques in order to obtain the optimal inputs to be applied to the plant (Maciejowski, 2002; Camacho and Bordons, 2004). In MPC an objective function is optimised over all possible input sequences, subject to equality and inequality constraints. The objective is expressed as a function of the states, outputs and inputs of the system. The states and outputs are predicted over a given prediction horizon, the inputs can be varied over the control horizon. In each time step only the first element of the calculated input sequence is applied to the system. The main equality constraint is the dynamic model. Inequality constraints are actuator limits (input constraints), state constraints or output constraints. The inputs of the model are engine speed and pump setting, the output is the machine speed. The model of the system is nonlinear, and is derived by Coen et al. (2008). The local static model of the machine speed S in function of the engine speed and pump setting is: S = ˇ(PS) × (RPM × PS)

(1)

where RPM is the rotation speed of the diesel engine (rounds per minute, rpm), PS is the pump setting (%, between −100 and +100) and ˇ is the local amplification which depends on the pump setting.  indicates the change of the engine speed–pump setting product. In other words, the change of machine speed equals the local amplification times the change of the engine speed–pump setting product. The dynamic model of the system is a Hammerstein model: the static nonlinearity of Eq. (1) is followed by second order dynamics (Coen et al., 2008). This model structure can easily be explained: the hydrostatic flow depends nonlinearly on the pump setting and the engine speed, whilst the machine dynamics can be situated between the hydrostatic flow and the machine speed. The nonlinearity between the actual hydrostatic flow and the pump setting can be dealt with by online estimating the amplification between both signals. This is illustrated further on in Eq. (3). In order to linearise the system, the input is transformed as follows: the input of the model is taken to be PSk and

229

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

RPMk , which is respectively the change of the pump setting and the change of the engine speed, at time instant k. The system is discretised at a 20 Hz sampling rate. The sampling rate is selected such that all the dynamics of the system can be captured. Since pump setting and engine speed both need to be constrained to a physical range, an integrator state for pump setting and engine speed needs to be added. That way the controller can restrict the control signals (otherwise the actual pump setting and engine speed are not known to the controller, since the input signals are the change of those variables). The linearised system thus becomes:



xk+1 = AL xk + BL uk

(2)

yk = CL xk + DL uk

where uk = [PSk ; RPMk ] and with the matrices AL , BL , CL and DL :

⎡

1

0

0

0

0

⎤

⎢0 1 0 0 0⎥ ⎢ ⎥ ⎢ ⎥ AL = ⎢ 0 0 1.9105 −0.9141 1 ⎥ ⎢ ⎥ ⎣0 0 1 0 0⎦ ⎡

0

0

0

0

1

⎤

0

ˇRPMk−1

CL = 0



0

min

(3)

7.7E − 5

DL = 0.0009ˇRPMk−1

xk ,yk ,uk

+

T

(yk − ry,k ) S(yk − ry,k ) +

k=1

M−1 

N 

T

(xk − rx,k ) Q(xk − rx,k )

k=1

uT k Ruk

(4)

k=0

ˇPSk−1

0.0035

N 

1

⎢ ⎥ 0 1 ⎢ ⎥ ⎢ ⎥ BL = ⎢ ˇRPMk−1 ˇPSk−1 ⎥ ⎢ ⎥ ⎣ ⎦ 0 0

is in this case defined as the change of engine speed and pump setting. This cost is chosen because it is related to the comfort of the operator. One of the goals of this cruise control system is to reduce the average engine speed, which would result in a lower noise emission. In order to minimise the engine speed an extra penalty term in the engine speed state is added to the objective function, next to the error term and the input cost term. The objective now mathematically states to follow the speed set point as accurately as possible, with a minimal engine speed, whilst taking into account the input cost (for operator comfort). The equality constraints are formed by the dynamic model of the system of Eq. (2). The inequality constraints are input constraints (thus on change of pump setting and engine speed) and state constraints (on the pump setting and engine speed state). The input constraints are imposed to guarantee operator comfort, the state constraints confine the pump setting and engine speed to their physical region. The setpoint error, the engine speed minimisation and the input cost are all written as quadratic functions. The optimisation problem at every time step thus becomes:

0.0009

0.0009ˇPSk−1





subject to:



xk+1 = AL xk + BL uk y = CL x + DL u

,

k = 0, . . . , N − 1

with ˇ the local amplification, and RPMk−1 and PSk−1 the previous value of respectively the engine speed and the pump setting. Note that the amplification between each of the inputs and the machine speed (situated in matrix BL ) depends on the machine state (RPMk−1 and PSk−1 ).

k k ⎧ k x ∈ X, k = 1, . . .,N ⎪ ⎨ k

2.2.

where ry,k is the desired machine speed, S ∈ R1×1 is the weight on the setpoint error, rx,k contains the desired state values, Q ∈ R5×5 is the weight on the state errors and R ∈ R2×2 is the input cost. The matrices Q, R and S have the following structure:

Objective

Any controller makes a trade-off between control effort and control performance. Control performance is how well the controller can satisfy the control specifications, that is how fast it can correct disturbances and how long it takes to reach the set point. Control effort on the other hand defines how much input or input variation is required to obtain the current performance. Input effort is of course especially important when the input is costly, such as a heating input, hence it is also called input cost. An important aspect of a cruise control system is that the operator is on board. In other words, the operator feels all the speed variations applied by the controller. Thus the ideal controller could be described as one that maximizes the control performance with a minimal input cost. Thus, the objective of the controller is to minimise the setpoint error (between actual machine speed and desired machine speed) with minimal input variation. The input cost

y ∈ Y,

(5)

k = 0, . . . , N − 1

k ⎪ ⎩ u ∈ U, k = 0, . . . , M − 1 k

⎡0 ⎢0 ⎢ ⎣0

Q=⎢0 0

⎤

0

0

0

0

wRPM

0

0

0⎥

0

0

0

0

0

0

⎥ ⎦ 0

0

0

0

0

0⎥,

 R=

wPS

0

0

wRPM

 ,

S = [werror ]

(6)

where wRPM is the weight of the deviation of the engine speed from the minimum, wPS is the cost of the change of the pump setting, wRPM is the cost of the change of the engine speed and werror is the cost of the difference between the desired and the actual machine speed. X, Y and U represent the feasible region for respectively the vectors xk , yk and uk . These regions

230

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

are defined by the following inequality constraints:

⎡ ⎤ ⎡ ⎤⎫  1 0 0 0 0 PSmin ⎪ ⎪  ⎢ RPM ⎥⎪ ⎬ ⎢ 0 1 0 0 0 ⎥ min ⎢ ⎥ ⎢ ⎥  X = x ∈ R5×1 ⎢ ⎥x ≥ ⎢ ⎥ ⎪ ⎣ −1 0 0 0 0 ⎦ ⎣ −PSmax ⎦ ⎪ ⎪ ⎪ ⎪ ⎪  ⎭ ⎩  0 −1 0 0 0 −RPMmax       1 speedmin Y = y ∈ R1×1  (7) y≥ −speedmax −1 ⎧ ⎡ ⎤ ⎡ ⎤⎫  1 0 PSmin ⎪ ⎪ ⎪ ⎪  ⎪ ⎪ ⎢ RPM ⎥⎬ ⎨ ⎢ 0 1 ⎥ min ⎢ ⎥ ⎢ ⎥  U = u ∈ R2×1 ⎢ ⎥u ≥ ⎢ ⎥ ⎪ ⎣ −1 0 ⎦ ⎣ −PSmax ⎦ ⎪ ⎪ ⎪ ⎪ ⎪  ⎭ ⎩  0 −1 −RPMmax ⎧ ⎪ ⎪ ⎪ ⎨

with PSmin the minimal pump setting (expressed in %), PSmax the maximal pump setting, RPMmin the minimal engine speed (expressed in rpm), RPMmax the maximal engine speed, speedmin the minimal machine speed, speedmax the maximal machine speed, PSmin the minimal change of pump setting (expressed in %), PSmax the maximal change of pump setting, RPMmin the minimal change of engine speed (expressed in rpm) and RPMmax the maximal change of engine speed. This description assumes that the feasible region is constant over the horizon. Of course this formulation can easily be extended to time-varying constraints.

3.

Theory of quadratic programming

f (x) =

subject to

x aT i aT x i

x

1 T T 2 x Gx + x d

= bi ,

i∈E

≥ bi ,

i∈I

(8)

where G is a symmetric n × n matrix. d, x and {ai }, i ∈ E ∪ I are vectors with n elements. The optimisation problem described by Eqs. (4) and (5) can be brought into this form by stacking all the optimisation variables xk , yk and uk (for all k) in one vector x. The optimisation problem of Eq. (8) is equivalent with: min maxL(x, ␭) = f (x) − x

␭



i (aT i x − bi )

(9)

i ∈ E∪I

where L is called the Langrangian. In the following sections some methods for solving equality-constrained and inequality-constrained problems are treated. More details can be found in Nocedal and Wright (1999) and Boyd and Vandenberghe (2004).

3.1.

min 21 xT Gx + xT d x

Equality-constrained quadratic programs

Denote the number of constraints m, the number of optimisation variables n, assume that the constraints are independent, and that m ≤ n. After all, if the number of constraints is larger than the number of optimisation variables, there is no exact

(10)

subject to Ax = b, where A is the Jacobian of constraints defined by A = [ai ]T i ∈ E.

(11)

Note that since the constraints are independent, A is of full row rank. Based on the Langrangian as formulated in Eq. (9), some optimality conditions can be derived. The Karush–Kuhn–Tucker (KKT) (Nocedal and Wright, 1999) conditions state that for a given solution x* , a vector ␭* must exist such that the following system of equations is satisfied:



G

−AT

A

0



x∗



∗

 =

−d b

 .

(12)

The system (12) can be rewritten in a form that is more suited for computation by expressing x* as x* = x + p, where x is some initial estimate of the solution and p is the desired step. By introducing this notation and rearranging the equations, we obtain:



An optimisation problem with a quadratic objective function and linear constraints is called a quadratic program. The general quadratic program (QP) can be stated as: min

solution possible. Moreover, the constraints are assumed to be independent. Based on (8), the quadratic program can be written as:

G

AT

A

0



−p ∗



  =

g c

(13)

,

where c = Ax − b,

g = d + Gx,

p = x∗ − x.

(14)

The matrix in (13) is called the KKT matrix. Note that this matrix is always indefinite due to its special structure. This severely complicates solving (13). Several methods to obtain p and ␭* exist, such as the range-space and the null-space method (Nocedal and Wright, 1999). These methods are all based on a decomposition of the KKT matrix. The structure of the matrix can be exploited to reduce the computational complexity. This requires however implementing the correct numerical solvers for the decomposition. These numerical aspects render implementing the decomposition a tedious task.

3.2.

Inequality-constrained quadratic programs

The methods presented in this section can easily be extended to the general case of a quadratic program with equalities and inequalities. First an Interior Point Method is discussed, then a short introduction to Active Set Methods is given.

3.2.1.

Interior Point Method

An IPM is basically an iterative method to solve the nonlinear KKT conditions. Based on the Langrangian formulation of Eq.

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

231

(9), the KKT conditions are: ∇x L(x∗ , ␭∗ ) = 0, x∗ = bi , aT i

∀i ∈ E,

aT x∗ ≥ i ∗ i ≥ 0, x∗ ∗i aT i

∀i ∈ I,

bi ,

(15)

∀i ∈ I, = 0,

∀i ∈ E ∪ I

x − bi , i ∈ I, the KKT conBy introducing a slack variable yi = aT i ditions can be rewritten as: Gx − A∗T  + d = 0, A∗ x − y − b = 0, yi i = 0,

(16)

i = 1, 2, . . . , m,

Fig. 2 – The flow chart of an Active Set Method.

(y, ) ≥ 0. with: A∗ = [ai ]T i ∈ I.

(17)

Given a current iterate (x,y,␭) that satisfies (y,␭) ≥ 0, a duality measure ␮ can be defined as: yT  1 yi i = m m m

=

(18)

i=1

The algorithm now calculates steps (x, y, ␭) to satisfy the first two equations of (16), and iteratively lower the duality measure  of Eq. (18). By eliminating y, this comes down to solving the following system:



G

−AT

A

−1 Y



x 



 =





−rd

−rb + −y + −1 e



(19)

where Y = diag(y1 , y2 , . . . , ym ), T

e = (1, 1, . . . , 1) ,

 = diag(1 , 2 , . . . , m ),

r d = Gx − AT  + d,

rb = Ax − y − b

(20)

 is the duality measure,  ∈ (0,1) is a tuning parameter which determines the convergence speed of the algorithm. The step size ˛ is defined such that the last equation of (16) remains satisfied. The next intermediate solution is then defined as: (x+ , y+ , + ) = (x, y, ) + ˛(x, y, ),

(21)

where ˛ is the previously defined step size and (x, y, ␭) is the calculated optimal step from the starting value (x,y,␭). More details about IPMs can be found in Nocedal and Wright (1999) and in Boyd and Vandenberghe (2004). Note that choosing an initial value is much easier for the IPM than for the ASM. For an IPM the feasible region is only limited by (y,␭) ≥ 0.

3.3.

Active Set Method

The scheme of an ASM is presented in Fig. 2. Only a brief outline of the method is given here. More details can be found in

the next section. First of all a feasible starting point x0 is calculated. A feasible point is a point that satisfies all the inequality constraints (and equality constraints if present). In every iteration k an intermediate solution xk is available. In the first iteration this is the starting point x0 . The active set is defined as the set of inequality constraints that would be violated if the method proceeds from the current solution xk to the unconstrained solution x* . In order to avoid this, the equalityconstrained problem (with the active set equations as equality constraints) is solved instead of the unconstrained problem. The next solution xk+1 is defined as xk+1 = xk + pk . If the step pk is different from zero, a step size ˛k is calculated. In the case that ˛k is smaller than 1, the blocking constraint is added to the active set. The blocking constraint is defined as the constraint that would be violated (first) if the step pk is taken. If the step pk is zero, the gap between the unconstrained solution x* and the current solution xk cannot be further decreased without violating one or more of the active set constraints. The Lagrange multipliers are then used to check whether all constraints are necessary. If the active set contains a constraint that would not be violated by a step towards the unconstrained optimum x* , this constraint is dropped from the active set. If all constraints are necessary, the current solution is the optimal solution.

3.4.

Final remarks

In literature, it is shown that IPMs scale better with increasing system sizes than ASMs. The number of iterations of ASMs increases with the system dimension, whereas the number of iterations of IPMs is known to increase with the root of the system dimension (Boyd and Vandenberghe, 2004). For small to mid-size systems however, ASMs generally perform faster. Both methods come down to solving an equalityconstrained problem, which requires factorisations (which are very hard to implement). But, since the matrix structure (location of zeros) remains constant for IPMs, standard sparse factorisation techniques can be used. For ASMs, the structure of the matrix depends on the current active set. This renders an efficient implementation very difficult. In the next section an ASM algorithm that avoids solving an equality-constrained problem will be presented.

232

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

nificantly. The dimension of the optimisation vector v only depends on the control horizon (often much smaller than the prediction horizon) and the input dimension. We assume that a feasible starting point v* (that satisfies all the constraints) is available. The problem of finding such a starting point is treated in the next section. In a standard ASM the equality constraints are not eliminated. However, in the case of MPC, the elimination is very straightforward, and the dimension reduction (of the optimisation vector) speeds up the online optimisation. Moreover, the elimination of these equality constraints can be performed offline.

4.3.

Fig. 3 – The flow chart of the Active Set Method presented in Section 4.

4.

Implementation of an Active Set Method

4.1.

Motivation

When solving QPs online, speed is essential. For small to mid-size problems ASMs are clearly faster than IPMs (from experiments, see further on). ASMs also have a nice geometrical interpretation, which makes them easier to understand. The major disadvantage of implementing an ASM is the complicated sparse decomposition, which is very hard to implement efficiently. In this paper an alternative implementation of ASM is presented, that avoids this problem. This method is a combination of the previously described ASM and the gradient-projection method (Nocedal and Wright, 1999). The following sections describe the algorithm, starting from the general optimisation problem defined in (8). An overview of the algorithm is given in Fig. 3.

4.2.

Eliminate equality constraints

First of all the equality constraints from (8) are removed. In the case of MPC, the equality constraints are the dynamic model, and give the relation between state and output on one hand, and initial state (current state of the system) and input on the other hand. Using this model, the objective can be written as a function of the initial state and the input. We then have an optimisation problem with only inequality constraints. The problem can be written as: min x

1 T T 2 v Hv + v r

subject to ␪T i v ≥ bi ,

i∈I

(22)

Note that the elimination of the constraints comes down to a transformation of the optimisation variable from x to v. This changes the objective function (G and d become respectively H and r) as well as the inequality constraints (vector ai becomes ␪i ). Eliminating the equality constraints reduces the problem dimension (and thus the computation time) sig-

Calculate the optimal step

The solution of (22) in iteration k + 1 can be written as vk+1 = vk + qk where v0 = v* and vk is the solution in the previous iteration. The solution of the unconstrained problem is vunconstr = vk + pk where pk is the optimal step for the k+1 unconstrained problem, and is calculated by solving (22) without constraints. This is completely analogous to a standard ASM. The algorithm starts with an empty active set. The next section explains how this active set is updated. Each member of the active set is a hyperplane, to which the solution vk+1 has to belong. Such a hyperplane is defined by a vector ␪i , perpendicular to all points in the hyperplane. In other words, if the step pk has a component along one of these vectors ␪i , this means the solution will move away from the hyperplane. The active set thus forms a basis for the infeasible region. The active set can be written as a matrix A† , where each row contains an orthonormal basis vector. To avoid moving into the infeasible region, all components of pk along any of the vectors contained in the rows of A† need to be removed. Assuming that the basis is orthonormal, the optimal step qk can be defined as: qk = pk − AT A† pk ,

(23)

where each row of A† contains a vector of the orthonomal basis defined by the active set. The step size ˛k is determined just as for the standard ASM:

 ˛k ∈ [0, 1],

˛k = min

1,

min

i∈W / k ,␪T qk <0 i

4.4.

bi − ␪T i vk ␪T i qk

 (24)

Transform the active set

The update of the active set is analogous to the standard ASM. A constraint is added to the active set, if the step size ˛k is smaller than 1. If the optimal step qk is zero, the Lagrangian multipliers are used to check whether all constraints are necessary. If not, the superfluous constraint is removed. The difference with the standard method is that the active set is transformed into an orthonormal basis. The matrix A (see Eq. (11)) is thus transformed into matrix A† using Gram-Schmidt orthonormalisation (Friedberg et al., 1989). Depending on the problem size, it is recommended to perform

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

4.6.

233

Stopping criterion

Different stopping criteria can be used. It should be taking into account that in a real-life application terminating the algorithm may also be necessary to avoid numerical problems. The criteria described below are implemented in the abovepresented method. A standard ASM will have very similar stopping criteria. Theoretical (if the impact of numerical errors is not considered) stopping criteria are:

Fig. 4 – Performance of the implemented ASM if the objective function is not rotational symmetric. A is the solution found by the algorithm. B is the actual solution. The horizontal line is a constraint.

• The number of constraints equals the degrees of freedom. In this case the solution is exactly defined, and no more iterations are necessary. • The maximal step (step size 1) is taken. The optimal solution is reached, and the optimiser can be terminated.

Following criteria are necessary because of time or numerical constraints: this orthonormalisation recursively (Huang et al., 2001) or to recalculate the basis at each update.

4.5.

Rotational symmetry

So far it was assumed implicitly that the objective function is rotational symmetric relative to its minimum. In fact the solver minimises the distance between the unconstrained solution and the intermediate solution. However, this is only correct if the objective is rotational symmetric. Fig. 4 illustrates what happens if the objective function is not rotational symmetric. The algorithm ends up in point A, which would be the solution if the objective had been rotational symmetric. Point B is the actual solution. This can be solved by performing another transformation of variables. Assume that the objective function is xT Hx. It can be transformed into a rotational symmetric problem by √ defining a new variable w = Hx such that xT Hx = wT Iw with I the identity matrix. The transformation performed on the √ −1 objective function and on the constraints then is x = ( H) w. Note that this transformation is only possible if H is well conditioned. Due to the origin of H (the objective function of the MPC after eliminating the equality constraints), this condition is satisfied in most situations. Moreover, bad condition of H should always be avoided with any method, since it will lead to very slow convergence. This can be illustrated by an alternative implementation of the algorithm. If instead of the optimal step, the direction of steepest descent is followed, an objective function that is not rotational symmetric is no issue. This approach implies some other problems (the step size is no longer contained in the interval [0,1] for instance), and the convergence speed depends on the condition of the problem. Experiments have shown that this gradient approach has more problems with a badly conditioned objective function then the earlier presented method with the transformation in order to obtain a rotational symmetric system. Since a standard ASM solves the equality-constrained problem exactly, the above-described transformation is not necessary for the standard method.

• If the intermediate solution no longer satisfies the constraints, the optimiser should be stopped. Rounding errors may cause the system to leave the feasible region. Typically the first constraint error is only of the size of the machine precision. If the optimiser is allowed to continue, the error will grow exponentially due to the fact that an invariant condition of the ASM (the method always needs to stay in the feasible region), is violated. • If the step in an iteration is very small, and no blocking constraint is present (which caused the algorithm to take a small step), the optimiser should be stopped as well. • The optimiser may also be stopped if time has run out. A possible method to implement this is to compare the available time left to the time needed for the last iteration. • Depending on the application and the circumstances, some of these criteria may be left out, and others may be added.

4.7.

Final remarks

Some special numerical precautions need to be taken to implement this algorithm (or any algorithm for that matter) in a real-life application. Some of the precautions necessary for this method when implemented in LabView are presented below. The calculated step size and the calculated step is truncated to a lower precision than the machine precision. This keeps small numerical errors from spreading throughout the calculation. The feasible region is made slightly smaller for the calculation of the step size ˛k by increasing the right hand side of the equation. This reduces the probability of exiting the feasible region due to numerical errors.

4.8.

Computation speed

In this section the presented ASM method is compared to an IPM method as presented in Section 3. Both methods are implemented in LabView. The comparison is made on a simple

234

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

Fig. 5 – The optimisation problem as described in Nocedal and Wright (1999). The ‘x’ marks the unconstrained optimum, the diamond marks the constrained optimum, the circles are the different initialisation values, the thick lines are the constraints, and the dashed lines are the contour lines of the objective function.

problem (Nocedal and Wright, 1999) (Fig. 5): 2

minq(x) = (x1 − 1) + (x2 − 2.5) x

subject to

tems is very difficult, since the allocation of the tasks to the different cores of the processor depends on the type of task. A possible speed criterion is the number of multiplications per second the processor can execute within a LabView for-loop (since the algorithm is also implemented in LabView). Both systems can perform about 7 × 107 multiplications per second. One optimisation is too short to get a correct measurement of the execution time by the system clock (especially in Windows). Therefore each optimisation is performed 1000 times, after which the average is calculated. The results are shown in Table 1. For one initialisation value the calculation time for a Windows operating system is also given. The ASM outperforms the IPM by at least a factor of 2 on this simple problem. In the application section (Section 6) it will be shown that for the real-life problem the difference between both methods is even larger. The ASM is hardly influenced by the choice of the initialisation value, unlike the IPM. LabView RT is clearly faster than Windows. Note however that IPM is almost a factor 10 times slower than the ASM on Windows. It is not clear what causes this difference with LabView RT.

5.

Initialisation of an Active Set Method

5.1.

Motivation

2

x1 − 2x2 + 2 ≥ 0, −x1 − 2x2 + 6 ≥ 0, x1 + 2x2 + 2 ≥ 0, x1 ≥ 0, x2 ≥ 0,

(25)

The calculation time of both methods depends on the initialisation value. Five different start positions are defined on Fig. 5. For a real-life implementation, as in Section 6, a realtime operating system (LabView RT) is used on a PXI system (National Instruments industrial PC). This system is better suited for control because of its predictable timing behaviour. The timing behaviour of a Windows system depends on the active processes and the user interaction (Craessaerts et al., 2005). The computation speed of both methods is compared on LabView RT as well as on Windows. The real-time operating system runs on a 2.0 GHz Pentium M 760 processor with 512 MB RAM. The Windows operating system runs on an Intel Core Duo processor (T2500) at 2 GHz with 1 GB RAM. Comparing the computation speed of both sys-

An ASM requires an initial point in the feasible region. For some problems this initial point may be available from physical insights, for other problems it needs to be calculated. This section shows how this initial value can be calculated, and how the above-presented ASM can be adapted to solve the resulting linear program. Since the initialisation has to be done in every time step, computation speed is equally important for the LP as for the QP. Experience has shown that for many applications it may be necessary to implement the LP-solver yourself to reach the necessary computation speed.

5.2.

Slack variables

Assume that the equality constraints have been eliminated, such that there are only inequality constraints left. The procedure presented below can be extended to the more general case of a problem with both equality and inequality constraints. To find a feasible initial point, slack variables are introduced on the constraints. The following linear program

Table 1 – The time required to solve the optimisation problem of Eq. (21) for IPM and ASM, for different initialisation values and different operating systems (LabView RT and Windows) Initialisation value

(0,0) (1,0) (2,0) (1,0.5) (3,1)

LabView RT

Windows

IPM

ASM

IPM

ASM

6.278E−04 5.432E−04 6.434E−04 3.887E−04 3.935E−04

2.068E−04 2.077E−04 2.073E−04 2.072E−04 2.090E−04

2.312E−03

3.440E−04

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

235

is then defined: mineT z subject to

(x,z) aT x + i zi i

≥ bi ,

i∈I

(26)

z≥0

zi is the slack variable of the ith inequality. The LP needs to be initialised from a feasible point (x,z) The slack vector z can however be calculated from x, thus initialisation is no issue here. The LP in (26) minimises the H1 norm of the slack vector. It is also possible to minimise the H∞ norm of the slack vector (by using one value z for all constraints) or to minimise the H2 norm (by solving a QP). Since a LP can be solved faster than a QP, the H1 norm and the H∞ norm are used the most.

5.3.

Linear programming

LPs are generally solved with the simplex method. This method is nothing more than a simplified version of an ASM. Thus with some small modifications, the above-presented ASM can be used for solving LPs as well. There are two main differences between the method for QPs and the method for LPs: • The optimal step of the unconstrained LP would be infinitely large. Thus instead of this optimal step, the steepest descent direction (the gradient) of the LP is used. • Closely related to the previous remark, the optimal step size is no longer defined for a LP (for a QP the optimal step size was 1). Thus the step size is allowed to be infinitely large. The stopping criterion then has to be modified: step size is no longer relevant for the termination of the algorithm. Concluding, by using a different calculation for the optimal step, and by removing the upper limit on the step size, the ASM can be used for solving LPs.

6.

Fig. 6 – (a) The machine speed response with the MPC controller. The full line is the response, the dashed line is the reference speed. (b) The inputs applied by the MPC controller. The full line is the engine speed (left axis), the dashed line is the pump setting (right axis). During the first part of the acceleration both inputs are increased, in the second part the pump setting is further increased to compensate for the decreasing engine speed.

Results

The MPC controller as described in Section 2 is implemented using the ASM of Section 4 and the initialisation method of Section 5. The objective of this system is to regulate the machine speed towards a desired value, whilst both input signals, namely the engine speed of the diesel engine and the pump setting, remain in the feasible region. To reduce the exhaust and noise emission, the engine speed should also be kept as low as possible. The optimisation problem is described in Section 2. In order to test the above-proposed optimisation methodology, the approach was implemented on the application example of Section 2. The cruise control system was evaluated for several steps in reference speed. A step from 0 to 6 km h−1 is shown in Fig. 6. This figure shows the machine speed response and the applied inputs. Note that the controller varies these inputs just like a human would do, given the same controller objective. In the first phase both the pump setting and the engine speed are increased to guarantee maximal acceleration. Once the machine speed approaches the set point, the

engine speed returns to its minimal value, whilst the pump setting further increases to compensate for the decrease in engine speed. In this application the operator comfort is guaranteed by limiting the change of rotation speed of the diesel engine and the pump setting separately. For other applications, this may not suffice. A possible additional constraint is to limit the sum of the change of rotation speed of the engine and the pump setting. It is clear that the maximal machine speed can only be reached if the engine speed is allowed to be larger than the minimum in steady state; adding some special precautions can easily accommodate for this situation. More details about the system identification and the controller itself can be found in Coen et al. (2008). However, the objective of this paper is to illustrate the advantages of the ASM implementation. The optimisation problem for the cruise control system is a system with five states and two inputs. The prediction horizon is taken to be

236

c o m p u t e r s a n d e l e c t r o n i c s i n a g r i c u l t u r e 6 3 ( 2 0 0 8 ) 227–236

ten samples. The control horizon is set to five samples. Since there are four state constraints per time sample (upper and lower limit on pump setting and engine speed) and four input constraints per time sample (upper and lower limit on the change of pump setting and engine speed), the optimisation problem has 40 constraints (in this case the state constraints remain constant once the input is zero). Both methods (ASM and IPM) are compared offline on an optimisation problem that occurred during real-time operation of the cruise control system. The ASM (from Section 4) solves the problem in 1.3 × 10−3 s whilst the IPM needs 1 × 10−1 s (both with LabView RT). Note that the controller in practice only has 5 × 10−2 s available to come up with the solution (20 Hz sampling rate), which means the IPM is useless for this application. The increasing gap between the ASM and the IPM can partly be attributed to the fact that the IPM uses a standard solver instead of a sparse solver (which becomes more important for larger problems). Note that LabView RT is a real-time operating system, which is ideally suited for rapid prototyping. On one hand relatively high-level programming is still possible, which enables a swift implementation of the above-presented method. On the other hand the code is compiled and run on a real-time operating system, which can guarantee that the code is executed on time. Contrary to a Windows system, there is no influence of external devices or uncontrolled processes (moving the mouse actually influences the time required to solve the optimisation problem). The ASM is also well suited for implementation on other lower level controllers. The main advantages of the method are: • Easy and transparent implementation: You do not have to be an optimisation expert to understand the algorithm because of its graphical interpretation. • Intermediate results: Contrary to an IPM (except for some of the most recent academic versions), an ASM can be configured to return an intermediate solution if time is up. • Fast execution: This paper shows that the method is very fast, and even faster than the IPM present in LabView.

7.

Conclusion

To implement a Model Predictive Controller (MPC) a quadratic programming (QP) problem needs to be solved in every time step. This paper is concerned with an easy and fast implementation of a QP solver. Control engineers often have to build a prototype application using programming environments, which often do not contain the suited optimisation algorithms. First the application, cruise control on a combine harvester, is presented. Then a short introduction to quadratic programming is given. The computation speed of these methods however depends on the implementation of good sparse solvers. To make the implementation more accessible (and quicker), an alternative Active Set Method is proposed, which does not require advanced sparse solvers. This method is compared to a standard Interior Point Method on a simple problem from literature. The ASM is at least 2 times faster than the IPM.

The ASM is then extended to solve linear programs, which are needed to determine a correct initialisation value of the ASM. Finally a real-life implementation on a combine harvester is presented. A cruise control system is implemented using the proposed ASM. The computation time of the ASM is almost 100 times smaller than that of the standard IPM. In conclusion, this algorithm enables a swift implementation, and is faster than other methods with a similar implementation complexity. Above all it is transparent, even to people from other fields.

Acknowledgements The authors gratefully acknowledge the cooperation of National Instruments, and especially Jim Nagle (National Instruments R&D, Austin, TX, USA) and Wouter Van Hoof (National Instruments support, Zaventem, Belgium). Tom Coen is supported by a scholarship of I.W.T.-Flanders (Institute for the Promotion of Innovation through Science and Technology). Jan Anthonis is funded as a Postdoctoral Fellow of the Research Foundation – Flanders (FWO). Josse De Baerdemaeker is a full professor at K.U. Leuven.

references

Bellemans, T., De Schutter, B., De Moor, B., 2006. Model predictive control for ramp metering of motorway traffic: a case study. Control Engineering Practice 14 (5), 441–466. Boyd, S., Vandenberghe, L., 2004. Convex Optimization. Cambridge University Press, Cambridge. Blasco, X., Martinez, M., Herrero, J.M., Ramos, C., Sanchis, J., 2007. Model-based predictive control of greenhouse climate for reducing energy and water consumption. Computers and Electronics in Agriculture 55 (1), 49–70. Camacho, Bordons, 2004. Model Predictive Control. Springer. Coen, T., Saeys, W., Missotten, B., De Baerdemaeker, J., 2008. Cruise control on a combine harvester. Biosystems Engineering 99 (1), 47–55. Craessaerts, G., Maertens, K., De Baerdemaeker, J., 2005. A Windows-based design environment for combine automation via CANbus. Computers and Electronics in Agriculture 49 (2), 233–245. EC (European Commission), 2004. Emissions for non-road mobile machinery, Directive 2004/26/EC. EC (European Commission), 2005. Noise emissions for outdoor equipment, Directive 2005/88/EC. EC (European Commission), 2006. Machinery, Directive 2006/42/EC. Friedberg, S., Insel, A., Spence, L., 1989. Linear Algebra. Prentice Hall, London. Huang, X., Caron, M., Hindson, D., 2001. Recursive Gram-Schmidt orthonormalization procedure and its application to communications. In: Proceedings of the Third IEEE Signal Processing Workshop on Signal Processing Advances in Wireless Communications, Taoyuan, Taiwan, pp. 340–343. Maciejowski, J.M., 2002. Predictive Control with Constraints. Prentice Hall, London. Nocedal, J., Wright, S.J., 1999. Numerical Optimization. Springer, New York. Pinon, S., Camacho, E.F., Kuchen, B., Pena, M., 2005. Constrained predictive control of a greenhouse. Computers and Electronics in Agriculture 49 (3), 317–329. Wang, Q., Chalaye, G., Thomas, G., Gilles, G., 1997. Predictive control of a glass process. Control Engineering Practice 5 (2), 167–173.