Automatic tuning of model predictive controllers based on multiobjective optimization

M. Francisco and P.Vega

Dpto. de Informática y Automática, ETSII de Béjar. Universidad de Salamanca (Spain)
mfs@usal.es; pvega@usal.es

Abstract - In this work a general procedure for tuning multivariable model predictive controllers (MPC) with constraints is presented. Control system parameters are obtained by solving a multiobjective optimization problem. The set of objectives includes controllability aspects, in terms of the H_∞ norms of some closed loop transfer functions of the system, and others related to the range of manipulated and controlled variables, expressed using the l₁ norm. Moreover, the use of multiple linearized models for tuning, allows for the specification of robust performance criteria through a set of constraints. The mathematical optimization for tuning all controller parameters is tackled in two iterative steps. First, integer parameters are obtained using a specific random search, and secondly a sequential programming based method is used to tune the real parameters. As a validation example, the tuning of the control system for the activated sludge process of a wastewater treatment plant has been selected.

Keywords - Model Predictive Control. Activated Sludge Process. Mixed Sensitivity Problem. Robust Control Theory. l₁ Norm.

I. INTRODUCTION

Model Based Predictive Control (MPC or MBPC) has become the leading form of advanced multivariable control in the process industries. The popularity of MPC is due to the successful results, the natural way of incorporating constraints, and its simplicity for operators. Excellent reviews of MPC and comparisons of commercial MPC are given in Maciejowski (2002) and Qin and Badgwell (2003).

MPC controllers have been tuned traditionally by a trial and error procedure, determining prediction and control horizons, input and output weights in the objective function, and constraints. The tuning task can be particularly difficult if the system is multivariable, since the whole set represents a formidable array of possible tuning combinations and also because many of these parameters have overlapping effects on the closed-loop performance and robustness. In these cases the advantages of using automatic MPC tuning methods is clear.

In the literature there are many works dealing with automatic tuning of MPC, but due to their complexity and the difficulty to perform analytical studies, the development of a general method is still a challenge.  Some researchers have proposed tuning methods for specific types of MPC, or leaving apart some particular aspects of the problem, for instance not considering the horizons (Li and Du, 2002). In other works, tuning methods have been developed considering simplifications in some way. For instance, Sridhar and Cooper (1997) focused on a single input, single output case with first order plus dead time models and Al-Ghazzawi et al. (2001) developed a tuning method based on a linear approximation between the closed loop predicted output and the tuning controller parameters. Some other papers that also deal with the tuning problem of MPC are Abu-Ayyad and Dubai (2006), Liu and Wang (2000) and Van der Lee et al. (2008).

Many researches pay also attention to the tuning of MPC by solving different optimization problems. Ali and Zafiriou (1993) proposed an off-line procedure for the non linear MPC tuning specifying time performance criteria, and incorporating a random grid search for obtaining the horizons. Li and Du (2002) developed also a simulation-optimization method based on fuzzy decision criteria for tuning only the weighting factors of the control variables. In Francisco et al. (2005) and Francisco and Vega (2006) the authors describe a new methodology for the on-line automatic tuning of the whole set of parameters of linear MPC, carried out by minimizing dynamical indexes as performance measures. An important drawback of these works is the required dynamical simulations within the optimization algorithm, making the procedure extremely slow.

Frequency domain methods for tuning linear optimal controllers have been studied since the beginning of 1980's (Doyle et al., 1992), and they are a good alternative to speed up MPC automatic tuning procedures avoiding dynamical simulations.  Based on that, the authors propose in Vega et al. (2007) a tuning procedure by solving a mixed sensitivity problem with constraints using a nominal model of the process. However, due to the use of a single linearized model, some problems of stability and robustness were detected in the presence of nonlinearities and load disturbances acting on the plant. For these reasons, the objective of this work is to propose a general methodology in which predictive controllers are tuned systematically using multiple models and robust theory fundamentals to guarantee robust performance (Morari and Zafiriou, 1989; Sideris and Rotstein, 1993; Tadeo et al.,1998).

Another reason to use frequency domain techniques for tuning is that they provide the way to analyse some properties of the feedback system that time domain techniques do not. For example, disturbance rejection within a range of frequencies, not only considering a specific set of disturbances. On the other hand, the MPC time domain formulation is intuitive, flexible and with easy incorporation of constraints, so the mixture proposed seems to be promising.

The developed tuning procedure allows for the evaluation of the control system parameters by solving a constrained multiobjective MINLP/DAE (Mixed Integer Non linear Programming with Differential and Algebraic Equations) problem. The selected objectives include controllability indexes based on the H_∞ norms of some weighted closed loop transfer functions. Other indexes related to the range of manipulated and controlled variables expressed using the l₁ norm have been considered to keep the process away from physical and operational bounds. Nevertheless, the tuning method is suitable not only for the unconstrained MPC, but also when there is a fixed set of active constraints. Moreover, multiple models have been considered to ensure robust performance of the closed loop process, in the face of high non linear dynamics and load disturbances. 

The approach has been validated on a simulated nonlinear model representing a real wastewater treatment plant (Manresa, Spain). The interest in selecting this type of plant lies on the complexity and variability of the biological processes involved, making this case study very suitable to test the tuning method.

The paper is organized as follows. First, the MPC formulation and basic relationships are presented. Then, the method for MPC automatic tuning is posed and explained in detail. The robust MPC tuning, particularly with the multiple models approach, and some considerations about constraints are presented, to follow with the optimization procedure. Finally, the activated sludge process model selected for validation is described, to end up with the tuning method results, conclusions and future extensions.

II. MPC FORMULATION

A. Basic MPC formulation

The basic MPC formulation consists of the on-line calculation of the future control moves by solving the following constrained optimization problem subject to constraints on inputs, predicted outputs and changes in manipulated variables.

(1)

where k denotes the current sampling point,  is the predicted output vector at time k+i, depending of measurements up to time k, r(k+i) is the reference trajectory, Δu are the changes in the manipulated variables, H_p is the upper prediction horizon, H_w is the lower prediction horizon, H_c is the control horizon, R and Q are positive definite matrices representing the weights of the change of control variables and the weights of the set-point tracking errors respectively. In this work the matrices R and Q are diagonal but not time dependent, so the error vector  is penalized at every point in the prediction horizon and the changes in the control signal Δu(k+i) are penalized at every point in the control horizon.

The problem (1) is a Quadratic Programming (QP) problem that gives a sequence of control moves Δu(k+i). The first component of this sequence is applied to the system in time k+1, and the optimization problem (1) is repeated at the next sampling time (receding horizon strategy).

The MPC prediction model used in this paper is a linear discrete state space model of the plant obtained by linearizing the first-principles nonlinear model of the process (Maciejowski, 2002). We assume full state measurement without loss of generality, and then the equations are:

(2)

where x(k) is the state vector, u(k) is the input vector and d(k) the disturbance vector. It is assumed that d(t) ∈ L₂ [0,∞) where L₂ [0,∞) is the set of continuous signals on [0,∞) that have finite 2-norm. Matrices A, B and B_d are of adequate dimensions.

Solving the optimization problem (1) without constraints, the following MPC control law can be obtained. The first element of vector Δu(k) will be applied to the process:

(3)

where

	(4)
	(5)
	(6)
	(7)
	(8)
	(9)
	(10)
	(11)

As it is a common practice, we assume that future measured disturbances remain constant at the last measured value:

(12)

See Maciejowski (2002) for a more detailed description of the control algorithm.

B. MPC basic relationships

The MPC control law (3) can be also expressed in the Laplace domain as (see block diagram of Fig. 1):

(13)

where K_i  are the transfer functions between the control signal and the different inputs (r(s), y(s), d(s)) which depend on the control system tuning parameters (Q, R, H_p , H_w and H_c). Note that the proposed MPC structure is a combined feedforward-feedback system.

Fig. 1: Nominal closed loop system

Consequently, taking into account control law (13) and the transfer function of the open loop system, the closed loop response can be obtained from

(14)

where  are the filtered disturbances , and the nominal transfer functions are denoted by G₀ and Gd₀

Equation (14) can be expressed substituting the sensitivity function and complementary sensitivity functions:

(15)

where

	(16)
	(17)

In order to state the automatic tuning problem for the nominal case, the sensitivity function S₀(s)R_d0(s) between the load disturbances and the outputs will be considered, as well as the control sensitivity transfer function M₀(s) between the load disturbances and the control signals when the reference is set to zero. Its calculation is straightforward applying block algebra to diagram of Fig. 1:

(18)

For simplicity, only the SISO case has been considered, but when several inputs and outputs are present, the closed loop system matrices can be also calculated in the same way, and the subsequent work could be generalized.

III. AUTOMATIC TUNING OF MPC

This section discusses the development of an automatic tuning procedure for the previously described MPC control system, which combines feedforward and feedback actions. The control structure is tuned to achieve good disturbance rejection and the proposed method provides the optimum MPC parameters R and the horizons H_p and H_c by solving a constrained multiobjective optimization problem. An extension of the method is straightforward to tune H_w and Q if desired.

The criteria for the optimization are some disturbance rejection measures based on the H_∞ and l₁ norms of some closed loop transfer functions of the system. The calculation of these transfer functions is possible not only in the unconstrained case, but also when the set of active constraints is fixed, as it will be shown in point V. The combined feedforward-feedback control system is tuned simultaneously due to the particular characteristics of the MPC.

A. Mixed sensitivity optimization problem

The problem of finding an optimal MPC is stated as a mixed sensitivity optimization problem that takes into account both disturbance rejection and control effort objectives in the same tuning function. The problem definition is:

(19)

where

(20)

c = (H_p, H_c, R) are the tuning parameters, and it is subject to the set of constraints explained below.

The dependence on s of the transfer functions N₀, S₀ and M₀ has been omitted for brevity.

Wp(s) and Wesf(s) are suitable weights to achieve closed loop performance specifications and to reduce the control efforts respectively. Note that control efforts rather than magnitudes of control are included in the objective function by considering the derivative of the transfer function M₀.

B. Performance constraints

In order to ensure disturbance rejection we need (considering normalized disturbances) the following equation to be satisfied in the disturbances frequency range:

(21)

where S₀(jω) is the frequency response of the sensitivity function, d(ω) is the disturbance spectra, and (ω₁, ω₂) is the disturbances frequency range. By choosing a weight Wp(s) satisfying

(22)

disturbance rejection can be assured imposing the following constraint in the optimization tuning procedure:

(23)

A typical choice for the weight Wp(s) is a rational function with one zero and one pole.

C. Limits on control and output variables

Although the unconstrained MPC model is considered, the maximum value of the manipulated and controlled variables (for the worst case of disturbances) can be constrained to be less than u_max and y_max respectively, by means of the l₁ norm and the following conditions:

(24)

D. Multiobjective optimization approach

The optimization problem for automatic tuning can be defined as a multiobjective one, considering constraints (24) as objectives to avoid possible feasibility problems. Mathematically it is stated as:

and constraints (23)

where the objectives are:

(25)

with the respective goals  f₁^*, f₂^*, f₃^*. In order to keep satisfying constraints (24) when the solutions do not reach the objectives exactly, goals are chosen in the following way:

(26)

IV. MPC ROBUST TUNING

In order to give robustness to the tuned MPC, and to extend the procedure to more realistic cases, let us assume that there is available a nominal model of the process and the corresponding uncertainty limits to represent the real plant.

A.  Robust Performance constraints

The first step is the selection of the region around the nominal point in which performance specifications are going to be satisfied. Assuming the classic multiplicative uncertainty, two families of models П_u and П_d  are defined respectively by means of:

	(27)
	(28)

where lm_u(s) and lm_d(s) are some bounded uncertainty functions. The sensitivity and complementary sensitivity functions for П_u and П_d  families  are given by:

(29)

From (17), a new transfer function can be defined, representing the open loop residual disturbance remaining after the feedforward compensation and due to the uncertainty.

(30)

Recall from (23) that the nominal performance condition for the closed loop system is:

(31)

In the presence of uncertainty, the robust performance is attained if for every d(t) L₂ [0,∞), the desired performance condition is satisfied for every perturbed model G_p ∈ Π_u and Gd ∈ Π_d. This sufficient condition for robust disturbance rejection can be stated analogously to the nominal case:

(32)

This can be also formalized as follows (Adam and Marchetti, 2004)

(33)

The condition (33) is more conservative than the equivalent one derived for just feedback control loops. The point is that in robust control theory it is frequently accepted that nominal performance plus robust stability implies robust performance, but according with equation (33) this is not the case for the MPC combined feedforward-feedback control system, because R_d in (33) is different from the nominal R_d0 in equation (31). The equation (33) reduces to the traditional feedback condition in the absence of feedforward action.

B.  Robustness with multiple models

In this work, multiple models belonging to the families П_u and П_d have been considered for the MPC robust tuning, and the MPC obtained is optimal in the region that they define. Open loop transfer functions G(s) and Gd(s) are easily calculated for every model by linearization (variable s is omitted for brevity and N is the number of models considered)

(34)

Within the optimization process,  the following sensitivity functions S_i must be calculated with the current controller values (c) and fixed transfer functions (34), including the nominal model.

(35)

Then, for every model, the following constraint must be imposed to guarantee robust performance:

(36)

Note that multiple models could be also considered in objectives (25) or any other aspect of the tuning methodology.

V. CONSIDERATIONS FOR CONSTRAINED MPC TUNING

One of the major points of MPC controllers is their ability to deal with constraints when solving the optimization problem (1). In this point we show that the MPC tuning method proposed could be also used when the set of active constraints is fixed. In this case, the constrained optimization problem (1) can be posed as:

(37)

where

(38)

subject to constraints on inputs, outputs and control moves in the form:

(39)

where E, F and G are matrices of suitable dimensions,

	(40)
	(41)

Problem (37) has the form of a typical QP problem, since all constraints (39) can be expressed linearly in terms of ΔU(k)(Maciejowski, 2002). The equivalent QP problem is:

(42)

where E = [Ω β] if only constraints on control moves are considered to simplify the expression.

When the set of active constraints is fixed, we could pose the optimization problem (42) subject to:

(43)

which could be solved by the theory of Lagrange multipliers as the following unconstrained optimization problem (where λ is the Lagrange multiplier):

(44)

In this case, the constrained predictive control law is also linear and some analysis and design can be performed.

Note that in practice, the assumption of a fixed active set of constraints is a bit strong, because at the transient regime, different sets of control variables can be temporarily activated and, consequently, the set of active constraints changes. Therefore, the MPC control law cannot be reduced to a single linear controller and the proposed methodology is only fully applicable for the unconstrained case. For that reason, the constraints (24) are used to keep the variables within the feasibility region. However, sometimes it is necessary to activate some constraints in order to use the whole system capability. In those cases, the consideration of (24) may deteriorate the controller performance.

V. ALGORITHM DESCRIPTION AND IMPLEMENTATION

The main difficulty when solving the generated optimization problem is the existence of real and integer variables (control and prediction horizons). In this work a two iterative steps algorithm has been proposed, combining a random search based on Solis and Wets (1981) for tuning horizons, and the goal attainment multiobjective method for tuning weight R, implemented in MATLAB.  Similar two steps approaches are presented in Leyffer (2001).

The multiobjective algorithm for the real part is stated as a sequential quadratic programming (SQP) problem that minimizes parameter γ, which represents the deviation of objectives (f_j) from goals (f_j^*)

(45)

where w_j are the weights for every objective and c = (H_p, H_c, R) are the tuning parameters. In this work the values of these weights are such that the importance of all objectives is the same.

The random search basically generates new horizons by adding and subtracting random integers to the current point, selecting the candidate solution with the lowest cost (Francisco and Vega, 2006; Vega et al. 2007).

For solving the MPC optimization problem, MPC Toolbox of MATLAB has been used, with some specific modifications (Maciejowski, 2002), implementing an extended state space representation.

VI. ACTIVATED SLUDGE PROCESS DESCRIPTION WITH MPC CONTROL

A. Plant description

The purpose of wastewater treatment plants is to process sewage and return clean water to the river. Activated sludge process is a very important part of the cleaning procedure in those plants, and its control is the objective of this work. A detailed description of the whole plant considered in this work is given in Moreno et al. (1992). The water treatment comprises the following basic steps:

a) The primary treatment that is divided into three parts. The first is a screening process to remove the gross solids from the water. The second part takes place in a vessel where the sand settles and wastewater surface is skimmed to remove oil and grease. A primary sedimentation is the last step of this stage. This process removes up to 50 % of the total polluting sewage load.

b) The secondary treatment is the activated sludge process. The mixed outlet stream from the primary sedimentation tanks is now passed to a central feeder channel, where recycled activated sludge is introduced and the flow is divided between six aeration tanks that form the reactor. In these tanks, the aerobic action of a mixture of microorganisms is used to reduce the substrate concentration in the water. A bacterial culture degrades the organic substrate converting it into inorganic products, more biomass and water. The dissolved oxygen required is provided by a set of aeration turbines.

c) Clarification. The effluent is feed into clarification tanks, where the activated sludge and clean water are separated. After this, the water contains approximately 10 % of the waste material. The water is discharged to the river. Between 25 % and 100 % of the settled activated sludge is recycled to re-inoculate the reactors during the residence time of the water in the aeration tanks.

B. Mathematical model of the activated sludge process

This work focuses solely on the activated sludge and the clarification processes (Fig. 2). These two parts of the plant have been modeled based on mass balances.

Fig. 2: Aeration tank and clarifier layout with an MPC controller

· Aeration tanks

The rate of change of the biomass, organic substrate and dissolved oxygen concentrations are given by:

	(46)
	(47)
	(48)

where x, s and c are the biomass, the substrate (the Chemical Oxygen Demand (COD)) and the dissolved oxygen (DO) concentration at the output of the aeration tanks (mg/l); xir and sir are respectively the inlet biomass and substrate (mg/l). μ_max is the maximum growth rate of the microorganisms, q is the inlet flow (m³/h), K_d is the kinetic coefficient of biomass decay by endogenous metabolism (1/h), K_c is the kinetic coefficient of biomass decay by biological waste, V₁ is the total useful volume for the six aeration tanks (m³), y is the yield coefficient between cellular growth and substrate elimination, f_kd is the yield coefficient between biomass endogenous and substrate contribution to the medium, c_s is the DO concentration at saturation, K_la is the mass transfer coefficient, fk₁ is the aeration factor which depends on the number and speed of working turbines, OUR is the oxygen uptake rate and K_s is the saturation constant.

For the rate of change of the biomass (x), the first term describes the biomass growth following the Monod model, the second describes cell death (as in the Volterra-Leslie modified model), the third describes the biological waste, and the final term quantifies the dilution effects (Moreno et al., 1992). For the rate of consumption of organic substrate (s), the first term expresses the decrease of the substrate through the activity of the biomass (Monod model), the second and third ones describe the transformation part of the dead biomass and biological waste into organic substrate, and the last term is the difference between the input and output substrate mass flows. Finally, for the dissolved oxygen concentration (c), the equation follows the classic literature: the first term is the rate of oxygen transferred to the water, the second describes the rate of oxygen used by the microorganisms (uptake rate), and the final term quantifies the dilution effects.

Algebraic equations for xir and sir are expressed as mass balances:

(49)

where x_i , s_i are the biomass and substrate at the influent,  q_i is the input flow to the process, and s_r , x_r , q_r  are the recycled concentrations and flow.

Equation for oxygen uptake rate is:

(50)

where K₀₁ is the yield coefficient between the cellular growth and the oxygen consumption rate.

Note that the assumption that the reactions take place in one perfectly mixed tank is considered and the effluent from the primary treatment is symmetrically distributed into the six reactors.

· Secondary clarifiers (settlers)

The operation of these elements is described by mass balance equations and an expression for the settling of activated sludge. The model takes into account the difference in settling rates between layers of increasing biomass concentration. This model has been developed by Moreno et al. (1992), in an attempt to capture the dynamic behaviour of the clarifiers but in a simpler way than the partial derivative equations model that describes the evolution of biomass with respect to the depth in the clarifier in the classical way:

	(51)
	(52)
	(53)

where x_d is the biomass concentration at the surface of the clarifiers, leaving the plant, q_sal is the flow of clean water at the output of the clarifiers, x_b is the biomass concentration in the second layer, q₂ is the activated sludge flow at the output, x_r is the biomass concentration at the bottom of the clarifier, vs  is the settling rate of the activated sludge, A is the area of the clarifiers, and l_d, l_b, l_r are the height of the first, second and third layer, respectively (Fig. 2). Note that the clarifier input flow q enters that unit at the second layer level.

The settling rate is calculated experimentally, the parameters are evaluated to fit a curve defined by experimental points:

(54)

The relations between the different flows are:

(55)

where q_p is the purge flow.

Some of the parameters of this model are known, such as the volume of the aeration tanks of the area of the settlers, but the rest of them are calculated to minimize a function expressed as the difference between the model output and real data form the plant when the same inputs are applied to both systems (Moreno et al., 1992). The parameters calculated are described in Table 1.

Table 1: Parameters of the wastewater plant model

C. Control problem

The control of this process aims to keep the substrate at the output (s₁) below a legal value despite the large variations of the flow rate and the substrate concentration in the incoming water (q_i and s_i), which are the input disturbances and one of the main problems when trying to control the plant properly. Another control objective is to keep dissolved oxygen concentration (c₁) around 2 mg/l, concentration that is necessary for the proper working of activated sludge process. The set of disturbances used in dynamic simulations (Fig. 3) has been determined by COST 624 program and its benchmark (Copp, 2002).

Fig. 3:  Substrate disturbances at the influent (s_i, q_i)

As for the manipulated variables, the recycling flow qr₁ has been chosen (Moreno et al., 1992; Vega and Gutiérrez, 1999), with the possibility of considering also the purge flow (q_p). The biomass concentration is only a constrained variable for a good performance of the process and it is not controlled. In this work substrate control has been considered, although the methodology proposed is general and it could be also extended to oxygen control, by using the aeration factor (fk₁) as manipulated variable.

VII. TUNING RESULTS

The controller considered is a constrained linear MPC with sample period of T=0.5 hours, suitable for representing the process dynamics.  MPC constraints include the following limits on inputs and outputs:

(56)

The controlled variable is the substrate at the output (s₁), and the biomass concentration x₁ is required to be in a certain range for proper operation of the activated sludge process. The selected plant is fixed with dimensions V_1n=7668 m³ (reactor volume) and A_n=2970.88 m² (settler area) and a steady state point defined by s_1n=58.445 mg/l and qr_1n=220 m³/h (nominal values) Disturbances s_i and q_i are assumed to be measured and scaled to make methodology improvement clearer. Their nominal values are: s_in=340 mg/l and q_in=1150 m³/h.

The steps followed to obtain and compare all results in this paper are explained here:

• First, the MPC has been tuned as presented in this paper using some linearized models of the plant (state space models).
• Then the MPC has been tested with the linearized model of the plant, i.e. the prediction model and plant model are the same. In all figures and tables the corresponding variables are labelled with their name plus "linear". (e.g. s₁ linear, qr₁ linear)
• Finally, the MPC has been tested with the full nonlinear plant, using the differential equations of the activated sludge process. In all figures and tables the corresponding variables are labelled only with their names (e.g. s₁, qr₁)

In order to determine the multiple models for tuning the MPC, several criteria have been considered. The first one consists of changing the recycling flow, which in turn changes the steady state output substrate (s₁). The MPC obtained in this way is robust to this type of variations, satisfying the performance criteria inside the region.

The second criterion to obtain the multiple models consists of modifying the influent characteristics (input flow q_i and substrate concentration at the input s_i). This is very interesting because the plant influent has always a large variability and in this way the MPC tuned is robust to those variations.

Weights Wp(s) and Wesf(s) for the mixed sensitivity indexes are kept constant in order to evaluate the influence of multiple models. The weight Wp(s) has been chosen considering condition (22) with benchmark and real disturbance spectra, while Wesf(s) has been chosen to impose a certain penalty on control moves. Their values for disturbance s_i are the following:

	(57)
	(58)

The performance indexes used in this work to tune the controllers are f₁ = ||N₀||_∞; f₂ = ||M₀||₁, and the objective vector is:

(59)

The optimization procedure is iterative in two steps as explained in previous paragraphs. For the random search part a multiobjective function must be defined weighting adequately all objectives and considering constraints as penalty terms. The extension to include objective f₃ is straightforward.

A. Results changing output substrate (s₁)

The first tuning case of study includes performance constraints (36) for two new models calculated changing the operation point 20 mg/l around the nominal value s_1n=58.445 mg/l. The optimal MPC obtained produces better disturbance rejection than the controller obtained considering only one model for those plants working even on the edge of the region (s_1n+20 mg/l, s_1n-20 mg/l).

In Table 2 some results are shown for the nominal working point. In the column headed "multiple models", the MPC parameters and some performance indexes are presented, while column headed "single model" presents MPC tuning results with only one model, for comparison. Control efforts are larger with the MPC tuned with multiple models because the controller must satisfy performance conditions not only for nominal model, but also for the near ones.  Nevertheless, for this working point, the constraint over ||W_p · S₀ · R_d0||_∞ is still satisfied in both cases.

Table 2: Results for automatic tuning with single and multiple models changing s₁, in the nominal point

Then MPC performance in the edge of the region has been studied to evaluate controller robustness. Numerical results are shown in Table 3, for the case of the working point s_1n-20 mg/l. This case is the most difficult to control because we demand a substrate level much lower than the nominal one. In this case, when the MPC has been tuned with only one model, the constraint over ||W_p · S · R_d||_∞ is no longer satisfied, unlike when the MPC is tuned with multiple models.  For nonlinear plant simulations an inner point in the region has been considered (s_1n-10 mg/l) and the disturbance rejection is also better for closed loop system with the MPC tuned with multiple models.

Table 3: Results for automatic tuning with single and multiple models changing s₁

In Fig. 4 a comparison of the sensitivity functions for the nominal working point is presented. In this case both controllers satisfy disturbance rejection because magnitude of S(s)·R_d(s) is smaller than magnitude of Wp_si^-1(s). But in Fig. 5 a comparison of the sensitivity functions for both closed loop systems working in the point s_1n-20 mg/l is shown, and it is clear that only the MPC tuned with multiple models satisfies the constraint imposed by Wp_si^-1(s). In Fig. 6, results considering the linearized plant are presented, and in Fig. 7 results with the nonlinear plant around the point (s_1n -10 mg/l) have been performed. Finally, magnitude of Wesf_si^-1(s) and s^.M are presented in Fig. 8 for both controllers.

Fig. 4: Magnitude of Wp_si^-1 (s olid line) and sensitivity functions for the nominal model with MPC tuned using multiple models (dotted line) and single model (dash-dotted line)

Fig. 5: Magnitude of Wp_si^-1 (solid line) and sensitivity functions for the model located in s_1n-20, with MPC tuned using multiple models (dotted line) and single model (dash-dotted line)

Fig. 6: Substrate concentration (s1 linear) for linearized plant and MPC tuned with multiple models (solid line) or single model (dash-dotted line), in the point s_1n-20 mg/l.

Fig. 7: Substrate concentration (s₁) and recycling flow (qr₁) for nonlinear plant and MPC tuned with multiple models (solid line) or single model (dash-dotted line)

Fig. 8: Magnitude of Wesf_si^-1  (solid line) and sensitivity functions s^.M for the model located in s_1n-20, with MPC tuned with multiple models (dotted line) and with single model (dash-dotted line)

In Fig. 9 some results with arbitrary but reasonably tuning parameters are presented (R=0.012;H_p=7;H_c=2) showing that the optimum parameters produce better disturbance rejection.

Fig. 9: Substrate concentration (s₁) for nonlinear plant and MPC tuned with arbitrary parameters (dash-dotted line) or optimal parameters (solid line)

B. Results changing influent characteristics (s_i, q_i)

The second case considers two multiple models obtained changing the nominal plant influent (s_in+120 mg/l, q_in+230 m³/h) and (s_in-120 mg/l, q_in-230 m³/h).

The optimal MPC tuned considering multiple models produces also better disturbance rejection even for worst case working point in the region (s_in+120 mg/l, q_in+230m³/h). For this case numerical results are shown in Table 4. Substrate concentrations and comparison of sensitivity functions are similar to the previous case.

Table 4:  Results for automatic tuning with single and multiple models changing s_i, q_i

The results presented in the tables have been validated solving the optimization problem with different starting points. Some previous analysis of the optimization problem was performed, showing that for MPC tuning applied to the activated sludge process, local minima hardly appear. Nevertheless, there is no guarantee that the optimal solution will be found in a general case.

VIII. FUTURE EXTENSIONS

The tuning method proposed in this work can be extended to other MPC structure formulation with stability guarantee.

For instance, following the Kwon and Pearson (1979) results, a straightforward application of the methodology proposed in this paper could be carried out considering a MPC formulation with a terminal constraint over a coincidence horizon that may be included as an additional tuning parameter. The implementation can be performed by using Lagrange multipliers as described in section V.

Furthermore, for the constrained MPC, we are considering now the formulation proposed by Rawlings and Muske (1993). In this paper they formulate an infinite horizon problem but keeping a finite number of decision variables in the optimization problem so it can be solved on line as a quadratic program. The advantage of this formulation is that nominal stability is ensured for any selection of the tuning parameters, even in the presence of constraints, if they are feasible. As the authors indicate, the design for other performance goals that can be posed as optimization over the tuning parameters since it is essentially an unconstrained optimization problem. Consequently, the proposed tuning procedure can be applied with minor modifications.

IX. CONCLUSIONS

In this work a method for tuning robust model predictive controllers has been developed, based on some frequency domain performance norm based indexes. The use of multiple linear models also allows for the specification of many robust performance criteria conditions. This method has been tested successfully in MPC applied to the activated sludge process, and the closed loop responses for substrate concentration in the reactor show that obtained controllers are properly tuned, taking into account the large magnitude of influent disturbances.

The methodology proposed here is based on a state space representation, but it is a general one, and any other linear MPC representation and performance criteria can be considered. Particularly, the extension to other structure formulations with stability guarantee is straightforward.

Finally it is important to show that the developed method is particularly suitable for its inclusion in the resolution of the Integrated Design optimization problem, which determines the optimum controller and the optimum plant at the same time.

ACKNOWLEDGMENTS
The authors gratefully acknowledge the support of the Spanish Government through the MEC project DPI2006-15716-C02-01 and Samuel Solórzano Foundation Project of the University of Salamanca (Spain).

REFERENCES
1. Abu-Ayyad, M. and R. Dubai, "MIMO extended predictive control-implementation and robust stability analysis," ISA Transactions, 45, 545-561 (2006).         [ Links ]
2. Adam, E.J. and J.L. Marchetti, "Designing and tuning robust feedforward controllers," Computers and Chemical Engineering, 28, 1899-1911 (2004).         [ Links ]
3. Al-Ghazzawi, A., E. Ali, A. Nouh and E. Zafiriou, "On-line tuning strategy for model predictive controllers," Journal of Process Control, 11, 265-284 (2001).         [ Links ]
4. Ali, E. and E. Zafiriou, "On the tuning of Nonlinear Model Predictive Control Algorithms," American Control Conference, 786-790 (1993).         [ Links ]
5. Copp, J.B. The COST Simulation Benchmark: Description and Simulator Manual, Office for Official Publications of the European Community (2002).         [ Links ]
6. Doyle, J.C., B.A. Francis and A.R. Tannenbaum, Feedback Control Theory, MacMillan (1992).         [ Links ]
7. Francisco, M., P. Vega and O. Pérez, "Process Integrated Design within a Model Predictive Control framework," IFAC'05 World Congress (2005).         [ Links ]
8. Francisco, M. and P. Vega, "Diseño Integrado de procesos de depuración de aguas utilizando Control Predictivo Basado en Modelos," Rev. Iberoamericana de Automática e Informática Industrial, 3, 88-98 (2006).         [ Links ]
9. Kwon, W.H. and A.E. Pearson, "On feedback stabilization of time-varying discrete linear systems," IEEE Trans. On Automatic Control, 23, 479-481 (1979).         [ Links ]
10. Leyffer, S., "Integrating SQP and Branch and Bound for Mixed Integer Nonlinear Programming," Computational Optimization and Applications, 18, 295-309 (2001).         [ Links ]
11. Li, S. and G. Du, "On-line tuning scheme for generalized predictive control via simulation-optimization," IEEE International Conference on Fuzzy Systems, 1381-1386 (2002).         [ Links ]
12. Liu, W. and G. Wang, "Auto-tuning procedure for model-based predictive controller," Proceedings IEEE SMC, 5, 3421-3426 (2000).         [ Links ]
13. Maciejowski, J.M., Predictive Control with Constraints, Prentice Hall (2002).         [ Links ]
14. Morari, M. and E. Zafiriou, Robust process control, Prentice Hall (1989).         [ Links ]
15. Moreno, R., C. De Prada, J. Lafuente, M. Poch and G. Montague, "Non-linear predictive control of dissolved oxygen in the activated sludge process," IFAC BIO 2, 289-298 (1992).         [ Links ]
16. Qin, S.J. and T.A. Badgwell, "A survey of industrial model predictive control technology," Control Engineering Practice, 11, 733-764 (2003).         [ Links ]
17. Rawlings, J.B. and K.R. Muske, "The stability of constrained receding horizon control," IEEE transactions on automatic control, 38, 1512-1516 (1993).         [ Links ]
18. Sideris, A. and H. Rotstein, "Single-input-Single-output H_∞ control with Time Domain Constraints," Automatica, 29, 969-983 (1993).
19. Solis, F.J. and R.J-B. Wets, "Minimization by random search techniques," Mathematics of Operations Research, 6, 19-30. (1981)         [ Links ]
20. Sridhar, R. and D. Cooper, "A tuning strategy for unconstrained SISO model predictive control," Ind. Eng Chem Res, 36, 729-746 (1997).         [ Links ]
21. Tadeo, F., T. Hologan and P. Vega, "l₁-optimal regulation of a pH control plant," Computers and Chemical Engineering, 22, S459-S466 (1998).         [ Links ]
22. Van der Lee, J.H., W.Y. Svrcek and B.R. Young, "A tuning algorithm for model predictive controllers based on genetic algorithms and fuzzy decision making," ISA Transactions, 47, 53-59 (2008).         [ Links ]
23. Vega, P., M. Francisco and E. Sanz, "Norm based approach for automatic tuning of Model Predicitive Controllers," Proceedings of ECCE-6, Copenhague (2007).         [ Links ]
24. Vega, P. and G. Gutiérrez, "Optimal Design Control and Operation of wastewater treatment plants," European Control Conference, Germany (1999).         [ Links ]

Received: September 22, 2008
Accepted: July 1, 2009
Recommended by Subject Editor: Jorge Solsona