Robust MPC for nonlinear multivariable systems

7 3URFHHGLQJVRIWKHWK0HGLWHUUDQHDQ&RQIHUHQFHRQ &RQWURO $XWRPDWLRQ-XO\$WKHQV*UHHFH Robust MPC for nonlinear multivariable systems Bouzouita Badreddine∗ , Bouani Faouzi§ , Ksouri Mekki∗ ∗ École Nationale d’Ingénieurs de Tunis, Tunisia. National des Sciences Appliquées et de Technologie, Tunis, Tunisia. E-mail : Badreddine.Bouzouita@enit.rnu.tn, Faouzi.Bouani@insat.rnu.tn, Mekki.ksouri@insat.rnu.tn § Institut Abstract— In this work, a robust predictive controller of uncertain nonlinear multivariable systems is developed. The control design is based on Multi-Input MultiOutput (MIMO) Nonlinear Auto Regressive Moving Average (NARMA) model. To cope with uncertain dynamic behavior of the system, the structured uncertainty is adopted. In fact, the main limitation of the robust predictive controllers is the computational burden leading to a lack of on line implementation. In this work, an efficient method is proposed. This method is based on transformation variables which reduce the initial non-convex problem to a convex programming. The efficiency of the proposed method is tested and compared with LMI and genetic algorithms optimizers on benchmark functions. The robustness of the proposed control law is experimented on three tanks system. I. P ROBLEM FORMULATION Linear MPC approach resorts to linear model to predict the future behavior of the system to be controlled. Although this controller have found successful applications [1], its success is restricted to given operating point. In fact, linear models are not able to describe the global behavior of the system over the whole operating range. This motivates the employ of nonlinear system description that leads to nonlinear model predictive control (NMPC) [2]. In this paper, the nonlinear MIMO NARMA model is adopted. The proposed representation has known considerable interest in applications due the fact that it deals with nonlinearity on input and output signals [3]–[5]. In addition to nonlinear behavior, many industrial process models are characterized by the presence of timevarying uncertainties such as unknown process parameters and external disturbances which, if not accounted for in the controller design, may cause performance deterioration and even closedloop instability. and inputs, respectively; ny and nu are their associated maximum lags; fj (.) are unknown nonlinear functions. Leontaritis and Billings have been demonstrated that polynomial representation of NARMA model work well in practical application [7]. Thus, fj (.) is expressed as a polynomial of degree L (where L is the degree of the nonlinearity): fj (x1 , x2 , ..., xn ) = N
θji i q1  q2  xp1 p1 =1 xp2 · · · p2 =1 qr  xpr pr =1 (2) with q1 , q2 , ..., qr ≥ 0 and q1 + q2 + ... + qr ≤ L Then, from equations (1) and (2), the outputs of the system can be rearranged in a compact form: Y (k) = θφ(k) (3) where θ is a matrix of scalar parameters and φ(k) represents the data vector which includes the past values of inputs and outputs:   θ10 · · · θ1N   .. (4) θ= , . θn0 ··· θnN φ(k) =[1, u1 (k − 1), ..., up (k − nu), y1 (k − 1), ..., yn (k − ny), ..., r uri i (k − p) × · · · × uj j (k − q), ..., ylrl (k − h) × · · · × urt t (k − s), ..., rm (k − h) × · · · × yprp (k − s)]T ym (5) Note that the polynomial model is nonlinear in the output and input variables but linear in the parameters. A. MIMO NARMA representation Therefore, the set of coefficients can be estimated by a Least-Squares (LS) algorithm. The presence of a numeric model is a necessary The robustification of the predictive controller consists condition for the development of the predictive control. to take into account, in an explicit manner, the uncertainSince, it permits to predict the future behavior of the ties at the time of calculation of the control law. Most process. To cope with nonlinear multi-variable systems predictive controllers consider additive uncertainties [8]– with n outputs and p inputs, the MIMO-NARMA model [11]. However, this type of uncertainties is limited to is adopted in this work. Indeed, this last provides a unified measure errors on output signals and it is not suitable to representation for a wide class of non-linear systems [6], describe the uncertain behavior of the physical system. In [7]. The outputs of the system are given by: the present paper, structured uncertainty is adopted. This yj (k) = fj (Y (k−1), ..., Y (k−ny), U (k−1), ..., U (k−nu)) last is adequate to model the uncertain dynamic behavior of the system [12], [13]. (1) where j = 1, ..., n; Y (k) = [y1 (k), ..., yn (k)] ∈ Rn To deal with uncertain behavior of physical system, and U (k) = [u1 (k), ..., up (k)] ∈ Rp are system outputs structured uncertainty is considered. Using a polytope 7 3URFHHGLQJVRIWKHWK0HGLWHUUDQHDQ&RQIHUHQFHRQ &RQWURO $XWRPDWLRQ-XO\$WKHQV*UHHFH description of uncertainties, the future outputs of the process are given by: r
αi θi )φ(k + j) (6) Y (k + j) = ( i=1 r αi = 1. where αi ≥ 0, i=1 B. Control law To indicate how well the process follows the desired trajectory, predictive control employs a cost function. This last depends on the future error between output signals and setpoints, and future increment controls to be optimized. In this work, the cost function to be optimized is as follows: Ny T [Y (k + i) − W (k + i)] [Y (k + i) − W (k + i)] J= i=1 N u−1 T [∆U (k + j)] R [∆U (k + j)] + j=0 (7) y1 (k + i) · · · yn (k + i) where Y (k + i) = w1 (k + i) · · · wn (k + i) , W (k + i) = and ∆U (k+j) = ∆u1 (k + j) · · · ∆up (k + j) are respectively the vectors of output predictions, setpoints and future increment of inputs; N y is the output prediction horizon, N u is the control horizon and R is the control weighting diagonal matrix. Robust predictive control is based on worst case strategy. The control represents the best solution for the worst case defined by the set of uncertain models [10], [14]. Hence, the input control is obtained by the resolution of the following min-max optimization problem: min − max J(∆U, θ) ∆U ∈Ω αj=1,...,r (8) where Ω is the set of constraints on the input/output signals. Robust predictive control algorithm suffers from a great computational burden leading to a lack of on line implementation [8], [15]. In fact, many authors have proposed the linear matrix inequality (LMI) optimization technique to solve the min-max problem [16]–[19]. This method requires a substantial computing time which increases exponentially with control horizon and the number of variables to be optimized [19]. To overcome this problem, Generalized Geometric Programming (GGP) can be used [20]. The adopted optimization method is addressed to problems of minimizing or maximizing a multivariate polynomial under polynomial constraints. Indeed, this kind of problem is encountered in a wide variety of applications in production planning, engineering design, risk management, etc [21], [22]. This paper is organized as follows: Section II presents the GGP optimization algorithm. The effectiveness of this algorithm opposite LMI and genetic algorithms is also exhibited in this section. In section III, the implementation and robustness of the proposed control law are demonstrated through simulations using three tanks system example. Finally in section IV, the conclusions are presented. II. G LOBAL SOLUTION In the present work, the Generalized Geometric Programming optimization method is adopted. The proposed method is addressed to solve non-convex problems of which the objective function and constraints are polynomials. The mathematical formulation of GGP is defined as follows [23]: T0 min cj zj j=1 subjet to            Tk hkq zkq ≤ lk , k = 1, . . . , K q=1 (9) α α α zj = x1 j1 x2 j2 . . . xnjn β β β zkq = x1 kq1 x2 kq2 . . . xnkqn xi ≤ xi ≤ xi where cj , hkq , lk , αpi and βkqj ∈ ℜ. zj and zkq are called signomial term. Usually the domain is xi ∈ ℜ+ . This is, however, no essential restriction since simple translations of the variables can often be used to fulfill the requirement for variables originally taking negative values. A. Convexification strategies Several methods have been proposed for solving this kind of problem. These methods are based on variable transformations and some other techniques [23]–[27]. Lemma 1: For a twice-differentiable function f (X) =  n i c i=1 xα i , X = (x1 , ..., xn ), xi ≥ 0, c, xi , αi ∈ ℜ, ∀i, let H(X) be the Hessian matrix of f (X). The determinant of H(X) can be expressed as [27]: n  
n det H(X) = (−c)n αi xinαi −2 1 − αi i=1 i=1 (10) Remark 1: if c ≥ 0, xi ≥ 0, and αi ≤ 0 (for all i), then det H(X) ≥ 0 [27]. Remark 2: if c ≤ 0, xi ≥ 0, αi ≥ 0 (for all i) and n (1 − i=1 αi ) ≥ 0, then det H(X) ≥ 0 [27]. Using remarks 1 and 2, we give the following propositions: Proposition 1: A twice-differentiable function f (X) =  n i c i=1 xα i is convex for c ≥ 0, xi ≥ 0, and αi ≤ 0 (for i = 1, . . . , n) [27]. Proposition 2: A twice-differentiable function f (X) =  n i is convex for c ≤ 0, xi ≥ 0, αi ≥ 0 (for c i=1 xα i n i = 1, . . . , n) and (1 − i=1 αi ) ≥ 0 [27]. The convexification strategy consists to transform each non-convex monomial term of the problem (9) to convex one. For instance, considering the following function: f (x1 , ..., xn ) = cxr11 xr22 ...xrnn , (11) then the convexification of this function depends of the sign of c and the values of ri : • if c > 0, and ri ≥ 0 then new variables Xi are introduced according to xi = exp(Xi ). Thus, the function given by (11) can be rewritten as: cxr11 xr22 ...xrnn = c exp(r1 X1 + r2 X2 + ... + rn Xn ), (12) 7 3URFHHGLQJVRIWKHWK0HGLWHUUDQHDQ&RQIHUHQFHRQ &RQWURO $XWRPDWLRQ-XO\$WKHQV*UHHFH • which is the exponentiel of affine function. Therefore, it is convex. n if c < 0, and ”ri < 0 or i=1 ri ≥ 1” then let 1 n r ≤ R. So, the function xi = XiR where i=1 i (11) is given by: cxr11 xr22 ...xrnn r1 R r1 R rn R = cX1 X2 ...Xn (13) the equation (8) can also be expressed by the equations (14) and (15) as follows: min − max J(∆U, θ) = min J ∗ (∆U ) ∆U αj=1,...,r ∆U (14) with J ∗ (∆U ) = max J(∆U, θ) αj=1...r which is convex by proposition 2. = min − J(∆U, θ) αj=1...r B. Evaluation of GGP In order to evaluate the effectiveness of generalized geometric programming, the performance of this last is compared with Linear Matrix Inequality optimization algorithm and Genetic Algorithm (GA) through solving a set of benchmark problems listed in the appendix. To avoid misinterpretation of the optimization results related to the choice of any particular initial points, each of the algorithms was run 100 times from random initial points. The following criteria summarize the results from 100 times minimization per function: • Errors: it is the sum of errors between the reached solution and the global minima given in the appendix. • The CPU time: is the total time (in second) put for 100 times minimization per function. In fact, GA optimization algorithm requires some parametrization. For these benchmark functions, it is configured as follows: • The real codification, arithmetic crossover and arithmetic mutation are used. • Number of individuals in initial population is equal to 50. • The algorithm is stopped when the maximal number of generation is reached, which equal to 100. Table I presents the computational results obtained by GGP presented in this work, ’GloptiPoly’ which based on LMI optimization algorithm [28], [29] and genetic algorithm. From the errors obtained by GGP and LMI, we can conclude that both algorithms converge to global minimum whatever the starting points. Whereas, it is not the case for GA. Indeed, the error in the case of ’Colville’ and ’Rosenbrok’ functions is great. This due to the fact that GA is based on stochastic rules and decisions. The CPU times illustrated in the table I show that the proposed GGP is more faster than LMI and GA. In fact, the computational burden of LMI increases exponentially with number of variables. However, the CPU time of GA depends on the number of generation and the number of individuals. Therefore, the proposed optimization method (GGP) is an alternative for control fast systems and solve nonconvex optimization problem. C. Implementation of control law The min-max optimization problem presented by equation (8) is bilevel [30]. It gives the solution of the best design in terms of future increments of control ∆U for the worst case defined by the uncertain model. Therefore, (15) Equation (15) maximizes the objective function with respect to the uncertain parameter θ, and after, minimizes it with respect to ∆U (equation (14)). From relation (2), we can prove that the output prediction yi (k + j) is non-convex and under polynomial shape. Therefore, the criterion J given by (7) is non-convex and it can be transformed to a convex function by using on each signomial term the correspond convexification rule as given in section II-A. After the transformation of the criterion, we can use a standard optimization technique to solve it. Consequently, the computation time is reduced and the global solution is reached. III. S IMULATION STUDY A. System description In this section the performance of the developed controller is tested on interconnected tank system depicted in figure (1). The process is composed of three cylindric tanks numbered from 1 to 3 which are connected through valves µ13 and µ23 . The valves µ10 , µ20 and µ30 are the emptying valves to the main tank. Tanks 1 and 2 of section equal to 0.049 m2 are fed into water respectively by respectively pump 1 and pump 2. The section of tank 3 is 0.0638 m2 . The level for each tank depends on the sum of the water flowing into and flowing off the tank that can be adjusted by the flow rate of the pump 1 and pump 2. Then, the system can be conveniently represented by: dh1 = q1 − q10 − q13 dt dh2 = q2 − q20 − q23 S2 dt dh3 = q13 + q23 − q30 S3 dt S1 (16) (17) (18) where hi is the tank level, q1 and q2 are the input flows, Sj is the section of tank j and qij represents the water flow rate from tank i to j ( i, j = 1, 2, 3), which, according to Torricellis rule, is given by:  (19) qij = sij µij sign(hi − hj ) 2g |hi − hj | with sij = 6.36 10−5 m2 is the section of valve, g is the gravity coefficient and µij ∈ {0, 1} (where µij = 0 denotes that the valve is close and µij = 1 indicates that the valve is open). Notice that qi0 (i = 1, 2, 3) represents the outflow rate with:  (20) qi0 = si0 µi0 2ghi 7 3URFHHGLQJVRIWKHWK0HGLWHUUDQHDQ&RQIHUHQFHRQ &RQWURO $XWRPDWLRQ-XO\$WKHQV*UHHFH TABLE I C OMPARISON RESULTS : GGP, LMI AND GA Functions Price Colville Booth Schwefel Rosenbrok Nb of variables 2 4 2 3 7 GGP 6.1506e − 06 1.6743e − 05 3.3001e − 09 1.8771e − 05 5.3232e − 04 Errors LMI 9.0222e − 11 2.0127e − 08 −2.8248e − 10 4.6966e − 10 8.2584e − 07 where si0 = 6.36 10−5 m2 is the section of corresponding valve and µi0 ∈ {0, 1} (0 for close and 1 for open). We aim to control the water levels of tanks 1 and 2 by adjusting the flow rate of pump 1 and 2. B. Modeling and identification The process, although non-differentiable, may be regarded as a hybrid system. Indeed, it has many possible state locations (h1 > h3 or h1 < h3 or h2 > h3 or h2 < h3 ). Furthermore, the dynamic of the system depends on the state of the valves (µij ). In this simulation, we assume that the state of valves µi0 (i = 1, 2, 3) can be modified at any time but the other valves are always open. A general inspection reveals that a linear second order system is a good representation for small variations of the inputs and of the outputs. This means that the global nonlinear model, after linearization, should provide a second order discrete time system. Hence, a MIMONARMA model with ny = 2, nu = 2 and L = 2 is identified for seven combinations of valves µi0 (i = 1, 2, 3). The parameters of these models are obtained off line using LS identification algorithm with sample time equals to 5s. Consequently, the uncertain MIMONARMA representing the global dynamic behavior of the system is constructed from these models as equation (6). GA 0.2880 233.1446 0.0112 1.7657 6.5038e + 03 GGP 8.291 8.322 5.347 6.840 22.122 CPU time (s) LMI 174.310 80.937 34.690 90.390 265.492 GA 32.1970 33.2680 32.3060 32.8870 34.5800 estimated model, we can conclude that the adopted representation allows to model this multivariable nonlinear system with a small modeling error. C. Results The simulation experiments have been performed in order to emphasize the robustness property of the proposed control scheme opposite the uncertain dynamic behavior of the process. Hence, during the simulations the valves µi0 have modified as mentioned in table (II). The references are also generated, as shown in figures (3) and (4), of a manner to test the robustness of the proposed controller opposite the state locations (h1 > h2 , h1 = h2 and h1 < h2 ) and opposite the tracking problem. TABLE II S TATES OF VALVES µi0 (i = 1, 2, 3) Valves µ10 µ20 µ30 0-150 1 1 1 151-300 0 0 1 Sample time 301-500 501-700 0 1 1 0 0 1 701-800 1 1 0 To attain the control gaols, the controller is designed as follows: • The sample time is 5s. • The cost function to be optimized is defined by the quadratic norm criterion given by the relation (7) with: Ny = 4, Nu = 1 and R = diag(2.2970 106 ) • Fig. 1. Three tanks system. Figure 2 depicts the evolution of water levels h1 and h2 obtained using random inputs in the case µi0 = 1. Comparing the water levels of the true system and the Since N u = 1, the constraints on the input signals are: 0 ≤ qj (k) ≤ 5 10−4 m3 s−1 , −5 3 −1 −4 10 m s ≤ ∆qj (k) ≤ 4 10−5 m3 s−1 ; j = 1, 2; Figure 3 plots the evolution of water levels (h1 and h2 ) and the flow rates (q1 and q2 ). As shown in this simulation, the robust nonlinear multivariable predictive control exhibit good performances. Indeed, although the change of the states of valves which accompany at iterations 151, 301, 501 and 701, the output signals arrive to reach the desired setpoints. In the second simulation, we aim to test the robustness of the developed controller opposite the tracking capability. From figure 4, we note that the water levels of tank 1 7 3URFHHGLQJVRIWKHWK0HGLWHUUDQHDQ&RQIHUHQFHRQ &RQWURO $XWRPDWLRQ-XO\$WKHQV*UHHFH −4 Input signals x 10 6 q1 q2 4 2 0 0 500 1000 1500 2000 2500 3000 3500 k 4000 True system Identified model 1 h1 0.8 0.6 0.4 0.2 0 0 500 1000 1500 2000 2500 3000 3500 k 4000 True system Identified model 1 h2 0.8 0.6 0.4 0.2 0 0 Fig. 2. 500 1000 1500 2000 2500 3000 3500 k 4000 Fig. 4. Closed loop responses of the system in the case of tracking setpoints. Identification of water levels h1 and h2 for µi0 = 1. A PPENDIX L IST OF BENCHMARK FUNCTIONS [31] Price (2 variables): P rice(x) = (2x31 x2 − x32 )2 + (6x1 − x22 + x2 )2 s.t. −10 ≤ xi ≤ 10 3 global minimums: x∗ = (0, 0); x∗ = x∗ = (1.464352, −2.506012); P rice(x∗ ) = 0 Fig. 3. Closed loop responses of the system. and tank 2 flow the desired reference despite the change of the dynamic behavior of the process. In these simulations, the average time required to compute the control inputs q1 and q2 for each sample time is 0.11s. Therefore, we can conclude that the proposed optimization method, generalized geometric programming, presents an alternative to control fast systems. IV. CONCLUSIONS This paper has proposed the robust nonlinear multivariable predictive control based on multivariable NARMA model with structured uncertainties. This kind of uncertainty, we allow to deal with uncertain dynamic behavior of the system. The control law is formulated as nonconvex min-max problem. An efficient optimization technique is presented to overcome this problem. The efficiency of the proposed algorithm is tested on benchmark functions and compared with LMI and genetic algorithms optimizers. The obtained results show the superiority of the proposed algorithm. The robustness of the proposed controller has been tested on three tanks system. The obtained simulation results have showed that the developed controller can deal with the uncertain physical behavior of the process. (2, 4); Colville4 (4 variables): Colville(x) = 100(x2 − x22 )2 + (1 − x1 )2 + 90(x4 − x23 )2 + (1 − x3 )2 + 10.1((x2 − 1)2 + (x4 − 1)2 ) + 19.8(x2 − 1)(x4 − 1); s.t. −10 ≤ xi ≤ 10 1 global minimum: (x)∗ = (1, 1, 1, 1); Colville(x∗ ) = 0 Booth (2 variables): Booth(x) = (x1 + 2x2 − 7)2 + (2x1 + x2 − 5)2 ; s.t. −10 ≤ xi ≤ 10 1 global minimum: (x1 , x2 )∗ = (1, 3); Booth((x1 , x2 )∗ ) = 0 Schwefel 3.2 (3 variables): Schwef el(x) = (x1 − x22 )2 + (1 − x2 )2 + (x1 − x32 )2 + (1 − x3 )2 ; s.t. −10 ≤ xi ≤ 10 1 global minimum: (x1 , x2 , x3 )∗ = (1, 1, 1); Scwef el((x1 , x2 , x3 )∗ ) = 0 Extended Rosenbrok (7 variables): 7 Rosenbrok(x) = 100(xi − xi−1 )2 + (1 − xi−1 )2 ; i=2 s.t. −10 ≤ xi ≤ 10 1 global minimum: (x1 , ..., x7 )∗ Rosenbrok((x1 , ..., x7 )∗ ) = 0 = (1, ..., 1); R EFERENCES [1] S. Qin and T. Badgwell, “An overview of industrial model predictive control technology,” Fifth International Conference on Chemical Process Control, vol. 93, pp. 232–256, 1997. 3URFHHGLQJVRIWKHWK0HGLWHUUDQHDQ&RQIHUHQFHRQ 7 &RQWURO $XWRPDWLRQ-XO\$WKHQV*UHHFH [2] S. L. Oliveira and M. Morari, “Contractive model predictive control for constrained nonlinear systems,” IEEE Transactions on Automatic Control, vol. 45, no. 6, pp. 1053–1071, June 2000. [3] C. Loh and J. Duh, “Analysis of nonlinear system using narma models,” JSCE - Structural Eng./Earthquake Eng, vol. 13, pp. 11– 21, 1996. [4] N. Chiras, C. Evans, and D. Rees, “Nonlinear gas turbine modeling using narmax structures,” IEEE Trans. on Instrumentation and Measurement, vol. 50, no. 4, pp. 893–898, Aug 2001. [5] A. Leva and L. Piroddi, “Narx-based technique for the modeling of magnetorheological damping devices,” Smart Materials and Structures, vol. 11, no. 2, pp. 79–88, 2002. [6] I. Leontaritis and S. Billings, “Input-output parametric models for non-linear systems. part 2: Stochastic non-linear systems,” International Journal of Control, vol. 41, no. 2, pp. 329–344, 1985. [7] I. J. Leontaris and S. A. Billings, “Input-output parametric models for non-linear systems. part 1: Deterministic non-linear systems,” International Journal of Control, vol. 41, no. 2, pp. 303–328, 1985. [8] P. O. M. Scokaert and D. Q. Mayne, “Min-max feedback model predictive control for constrained linear systems,” IEEE Trans. Automat. Contr., vol. 43, pp. 1136–1142, 1998. [9] P. Grieder, P. A. Parrilo, and M. Morari, “Robust receding horizon control - analysis and synthesis,” Proceedings of the 42nd IEEE Conference on Decision and Control, Maui, Hawaii USA, vol. 42, pp. 941–946, December 2003. [10] D. Ramirez, M. R. Arahal, and E. F. Camacho, “Min-max predicitve control of a heat exchanger using a neural network solver,” IEEE transactions on control systems technology, vol. 12, no. 05, pp. 776–786, 2004. [11] D. M. de la Pena, T. Alamo, D. Ramı́rez, and E. Camacho, “Minmax model predictive control as a quadratic program,” Proceedings of 16th IFAC World Congress, Praga, Repblica Checa, 2005. [12] S. Ploix, “Diagnostic des systèmes incertains : l´approche bornante.” PhD thesis, Institut National Polytechnique de Lorraine, 23 décembre 1998. [13] O. Adrot, “Diagnostic à base de modèles incertains utilisant l´analyse par intervalles : l´approche bornante,” PhD thesis, Institut National Polytechnique de Lorraine, 4 décembre 2000. [14] D. Ramirez, T. Alamo, and E. Camacho, “Efficient implementation of constrained min-max model predictive control with bounded uncertainties,” Proceedings of the 41th IEEE Conference on Decision and Control, Las Vegas, Nevada USA, 2002. [15] J. H. Lee and Z. Yu, “Worst-case formulations of model predictive control for systems with bounded parameters,” Automatica, vol. 33, no. 5, pp. 763–781, 1997. [16] M. V. Kothare, V. Balakishnan, and M. Morari, “Robut constrained model predictive control using linear matrix inequalities,” Automatica, vol. 10, no. 32, pp. 1361–1379, Octobre 1996. [17] F. A. Cuzzola, J. C. Geromel, and M. Morari, “An improved approach for constrained robust model predictive control,” Automatica, vol. 38, pp. 1183–1189, 2002. [18] Z. Wan and M. Kothare, “An efficient off-line formulation of robust model predictive control using linear matrix inequalities,” Automatica, vol. 39, pp. 837–846, 2003. [19] A. Casavola, D. Famularo, and G. Franzé, “Robust constrained predictive control of uncertain norm-bounded linear systems,” Automatica, vol. 40, pp. 1865–1876, 2004. [20] B. Bouzouita, F. Bouani, and M. Ksouri, “Efficient implementation of mpc with parametric uncertainties,” accepted in European Control Conference, Kos, Greece, 2-5 July 2007. [21] K. Das, T. Roy, and M. Maiti, “Multi-item inventory model with under imprecise objective and restrictions: a geometric programming approach,” Production Planning and Control, vol. 11, no. 8, pp. 781–788, 2000. [22] C. F. C.D. Maranas, “Global optimization in generalized geometric programming,” Computers and Chemical Engineering, vol. 21, no. 4, p. 351369, 1997. [23] J.-F. Tsai and M.-H. Lin, “An optimization approach for solving signomial discrete programming problems with free variables,” Computers and Chemical Engineering, vol. 30, pp. 1256–1263, 2006. [24] K.-M. Björk, P. O. Lindberg, and T. Westerlund, “Some convexifications in global optimization of problems containing signomial terms,” Computers and Chemical Engineering, vol. 27, pp. 669– 679, 2003. [25] P. Shen and K. Zhang, “Global optimization of signomial geometric programming using linear relaxation,” Applied Mathematics and Computation, vol. 150, pp. 99–114, 2004. [26] Y. Wang and Z. Liang, “A deterministic global optimization algorithm for generalized geometric programming,” Applied Mathematics and Computation, vol. 168, pp. 722–737, 2005. [27] J.-F. Tsai, M.-H. Lin, and Y.-C. Hu, “On generalized geometric programming problems with non-positive variables,” European Journal of Operational Research, vol. 178, pp. 10–19, 2007. [28] D. Henrion and J.-B. Lasserre, GloptiPoly: Global Optimization over Polynomials with Matlab and SeDuMi, 11 March 2005, version 2.3.0, www.laas.fr@∼henrion, www.laas.fr@∼lasserre. [29] J. B. Lasserre, “Global optimization with polynomials and the problem of moments,” Society for Industrial and Applied Mathematics, vol. 11, no. 3, pp. 796–817, 2001. [30] R. F. Alfredo, L. Cardozo, and M. Ijar, “On the structural optimization for space applications using the min-max technique,” 55th International Astronautical Congress, Vancouver, Canada, 2004. [31] R. J. V. Iwaarden, “An improve unconstrained global optimization algorithm,” Doctor of Philosophy, University of Colorado at Denver, 1996.