This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

1

Robust Model Predictive Control of Nonlinear Systems With Unmodeled Dynamics and Bounded Uncertainties Based on Neural Networks Zheng Yan, Student Member, IEEE, and Jun Wang, Fellow, IEEE

Abstract— This paper presents a neural network approach to robust model predictive control (MPC) for constrained discretetime nonlinear systems with unmodeled dynamics affected by bounded uncertainties. The exact nonlinear model of underlying process is not precisely known, but a partially known nominal model is available. This partially known nonlinear model is first decomposed to an affine term plus an unknown highorder term via Jacobian linearization. The linearization residue combined with unmodeled dynamics is then modeled using an extreme learning machine via supervised learning. The minimax methodology is exploited to deal with bounded uncertainties. The minimax optimization problem is reformulated as a convex minimization problem and is iteratively solved by a two-layer recurrent neural network. The proposed neurodynamic approach to nonlinear MPC improves the computational efficiency and sheds a light for real-time implementability of MPC technology. Simulation results are provided to substantiate the effectiveness and characteristics of the proposed approach. Index Terms— Extreme learning machine (ELM), real-time optimization, recurrent neural networks (RNNs), robust model predictive control (MPC), unmodeled dynamics.

I. I NTRODUCTION

M

ODEL predictive control (MPC) is a powerful modelbased control technique, which explicitly optimizes the overall performance of a system to be controlled. In MPC, the control laws are obtained by solving a finite-horizon constrained optimization problem during each sampling interval, using the current state as an initial state. The optimization yields an optimal control sequence to the plant. There are several attractive features of MPC; e.g., it handles multivariable control problems naturally, it considers input and output constraints, and it adapts structural changes [1]. Manuscript received June 29, 2012; revised February 22, 2013, April 13, 2013, and June 25, 2013; accepted July 26, 2013. This work was supported in part by the Research Grants Council of the Hong Kong Special Administrative Region under Grant CUHK417209E, Grant CUHK416811E, and Grant CUHK416812E, and in part by the National Natural Science Foundation of China under Grant 61273307. Z. Yan is with the Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, Hong Kong (e-mail: [email protected]). J. Wang is with the Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, Hong Kong, and also with the School of Control Science and Engineering, Dalian University of Technology, Dalian 116024, China (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TNNLS.2013.2275948

A key issue for MPC implementation lies in the effectiveness and efficiency of real-time optimization. The computational efficiency determines the reliability of any MPC approach. In the literature, several approaches are available with the aim of reducing the computational burden of nonlinear MPC. The first category is MPC with online linearization. As a result of linearization, the MPC optimization task becomes a quadratic programming problem, which can be solved very efficiently [2]–[5]. The second category is MPC, in which the neural network approximates the time-varying feedback law. The control signal is calculated explicitly without any online optimization [6]. The third category is MPC, in which the neural network directly calculates the control signal without any optimization [7], [8]. The fourth category is MPC with an explicit piecewise linear state feedback approximator, which can be found offline using a multiparametric nonlinear programming. The controller is realized online by binary tree search [9]. One common limitation for most of the existing methods is that they may not be competent for problems with high dimensionality and stringent computation time requirement. In the past two decades, recurrent neural networks (RNNs) emerged as promising computational tools for solving various optimization problems. The essence of neurodynamic optimization lies in its inherent nature of parallel and distributed information processing and the availability of hardware implementation. Based on the duality and projection methods, many neural network models have been proposed for convex optimization (e.g., [10]–[16]) and pseudoconvex optimization (e.g., [17], [18]). These neural networks have shown good performance in terms of global convergence, low computational complexity, and good robustness [19]. The recent advance of RNNs offers a new optimization paradigm for real-time applications, such as MPC. MPC based on the neural networks was probably the first time proposed in [20] where the use of quadratic programming to implement nonlinear MPC with RNNs was described. Many studies on incorporating the neural networks with MPC have been carried out in the last two decades. Generally speaking, the applications of neural networks for MPC fall into three categories: 1) using neural networks for system modeling [21]–[26]; 2) using neural networks for solving online optimization problems [27]–[31]; and 3) using neural networks to approximate offline the MPC law [8]. In these works, the use of neural networks in MPC design

2162-237X © 2013 IEEE

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 2

greatly improved the computational efficiency and control performance. As MPC relies on the explicit use of system models, the representation of the underlying process is of vital importance. A nonlinear model generally results in a nonconvex optimization problem, which is computationally intractable. Jacobian linearization enables the MPC to be synthesized by solving a sequence of convex optimization problems. Because the Jacobian linearization approximates a nonlinear system locally, it would result in an unsatisfactory performance, particularly in the presence of model uncertainties. Another important, but a challenging issue in MPC is achieving robustness against uncertainties. Two types of uncertainties often exist, one is parametric uncertainties of the system structure and the other one is persistent disturbances. Some formulations have been proposed to address the robustness issues. The first method to achieve robustness is to minimize a nominal performance index with tightened state and terminal constraints [32], [33]. The second method is to solve a minimax optimization problem. Specifically, in an open-loop formulation, a sequence of control actions is obtained by minimizing the performance index under the worst case effect of uncertainties [34]–[37]. In a closed-loop formulation, a vector of feedback control policies is computed [38], [39]. The third method is to calculate the disturbance invariant sets such that the evolution of closed-loop trajectory lie in a invariant tube [40]–[42]. In this paper, a minimax MPC approach for nonlinear systems with unmodeled dynamics and bounded uncertainties is proposed. Considering Taylor expansion at an operating point, a nonlinear model can be decomposed to an affine model plus an unknown high-order term. The unknown higher-order term, together with the unmodeled dynamics are modeled by a feedforward neural network called extreme learning machine (ELM) [43] via supervised learning. The minimax approach is adopted to deal with bounded uncertainties. The minimax optimization problem is reformulated as a convex nonlinear minimization problem. A two-layer RNN is then applied for solving the nonlinear optimization problem in real time. The contribution of the proposed MPC approach is twofold. First, it can significantly reduce computational burden. Using successive linearization, a local linear model is obtained at each sampling instant. The nonlinear part and unmodeled dynamics are approximated using an ELM. Therefore, the entailing optimization problem can be formulated as a convex optimization. The issue of heavy computational burden results from minimax formulation is addressed. Second, the proposed robust MPC scheme is real-time implementable. Despite of the progress made in the research of nonlinear MPC, there are few approaches that are suitable for a real-time implementation in practice [3]. Because of the inherent parallelism of RNNs, the proposed scheme would not result in lower computational efficiency even if the optimization problem is of large size. In contrast, most of the existing MPC algorithms are computationally demanding when dealing with large-scale multivariable industrial processes. The rest of this paper is organized as follows. In Section II, the minimax MPC is formulated. In Section III, neural

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

networks are applied for modeling the unknown parameters and solving the optimization problems. In Section IV, simulation results are provided. Finally, Section V concludes this paper. II. P ROBLEM F ORMULATION A. Minimax Problem Formulation Consider a discrete-time nonlinear time-invariant model in the form of x(k + 1) = f (x(k), u(k)) + w(x(k), u(k)) + v(k) y(k) = C x(k)

(1)

m

is the state vector, u ∈ is the input where x ∈ vector, y ∈  p is the output vector, and w ∈ n is a vector of unmodeled dynamics, which can represent unidentified structures and unknown parameters of the plant, which may be state and input dependent. For the sake of brevity, we will not point out the functional dependence of w(x, u) except where strictly needed. v ∈ n is a vector of additive bounded uncertainties, f (·) is a nonlinear function, and C ∈  p×n is the output matrix. The following assumptions will be used throughout this paper. Assumption 1: All state variables are measurable or can be accurately estimated. Assumption 2: f (·) is Lipschitz continuous and differentiable for all (x, u), with f (0, 0) = 0. Assumption 3: w(x, u) is continuous at the origin with w(0, 0) = 0. Assumption 4: The additive uncertainty v is limited in a compact set, which contains the origin with known bound, i.e., v ∈ V, v2 ≤ ζ . In addition, its mean is assumed to be zero. Let L f be the Lipschitz constant of f . The difference between the state evolution without an additive uncertainty v and real evolution at a future time instant k + j is bounded j by ζ(L f − 1)/(L f − 1) [33]. To apply MPC on (1), a constrained finite-horizon optimal control problem is solved over a prediction horizon N during each sampling period, using the current state as the initial state. When bounded uncertainties are explicitly considered, the optimal control actions can be obtained by solving the following minimax optimization problem repeatedly over k: n

min max

u(k) v(k)

=

N  j =1

+

J (u(k), v(k), x(k), y(k)) y(k + j |k) − r (k + j )2Q j

N u −1 j =0

u(k + j |k)2R j + F(x(k + N|k))

subject to u min ≤ u(k + j |k) ≤ u max , j = 0, 1, . . . , Nu − 1 v min ≤ v(k + j |k) ≤ v max , j = 1, 2, . . . , N u min ≤ u(k + j |k) ≤ u max , j = 0, 1, . . . , Nu − 1 x min ≤ x(k + j |k) ≤ x max , j = 1, 2, . . . , N ymin ≤ y(k + j |k) ≤ ymax , j = 1, 2, . . . , N

(2)

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. YAN AND WANG: ROBUST MPC OF NONLINEAR SYSTEMS

3

where u(k) is the input increment vector defined as u(k) = u(k) − u(k − 1), r (k + j ) is a known output reference vector, y(k + j |k) is the predicted output vector, N and Nu are, respectively, prediction horizon (1 ≤ N) and control horizon (0 < Nu ≤ N), || · || denotes Euclidean norm, Q and R and weight matrices with compatible dimension (tuning methods are discussed in [44]), F is a terminal cost for stability, x(k + N|k) is the predicted terminal state within the prediction horizon, u min ≤ u max , u min ≤ u max , x min ≤ x max , and ymin ≤ ymax are lower and upper bounds of system constraints, and v min ≤ v max are the lower and upper bounds of the additive uncertainty. Because of the nonlinear relation between the predicted output vector y(k + j |k) and control increment vector, the optimization problem (2) generally becomes nonconvex, which is computationally intractable. It is desirable to reformulate (2) as a convex optimization. Via Jacobian linearization, the nonlinear system (1) can be decomposed, around an operating point [x o (k), u o (k)], to an affine system plus unknown terms x(k + 1) = f (x o (k), u o (k)) ∂ f (x o (k), u o (k)) (x(k) − x o (k)) + ∂x ∂ f (x o (k), u o (k)) (u(k) − u o (k)) + ∂u + ε f (k) + w(k) + v(k)

(3)

where ε f (k) is the high-order term of Taylor series of f . The operating point [x o (k), u o (k)] is chosen as the state and input vectors at the current time k, so the operating point at each time instant is known beforehand. Denote A(k) = ∂ f (x o (k), u o (k))/∂ x, B(k) = ∂ f (x o (k), u o (k))/∂u, δx(k) = x(k) − x o (k), δu(k) = u(k) − u o (k), ε(k) = ε f (k) + w(k) + f (x o (k), u o (k)) − x o (k), the system (1) can be reformulated as δx(k + 1) = A(k)δx(k) + B(k)δu(k) + ε(k) + v(k) δy(k) = Cδx(k). (4) Denote δr (k + j ) = r (k + j ) − C x o (k), then the MPC of the transformed model (4) can be formulated as the following optimization problem: min max

u(k) v(k)

=

N  j =1

+

J (u(k), v(k), δx(k), δy(k)) δy(k + j |k) − δr (k + j )2Q j

N u −1 j =0

u(k + j |k)2R j + F(δx(k + N|k))

where δu min = u min − u o (k), δu max = u max − u o (k) δx min = x min − x o (k), δx max = x max − x o (k) δymin = ymin − C x o (k), δymax = ymax − C x o (k). It is assumed that the problem (5) is feasible at the initial time k = 0. B. Input-to-State Practical Stability In recent years, some results on the stability of minimax MPC were presented [45]– [47]. It was shown in general that only input-to-state practical stability (ISpS) can be ensured for minimax MPC of nonlinear systems with persistent disturbances [45]. Compared with input-to-state stability, ISpS is a weaker property, which does not guarantee asymptotic stability when disturbances vanish. The ISpS definitions, assumptions, and conditions for MPC of system (1) are stated as follows. Definition 1: A continuous function αˆ : + → + is a K-function if α(0) ˆ = 0, α(s) ˆ > 0 for all s > 0, and α(s) ˆ is strictly increasing. A continuous function βˆ : + × Z + → ˆ 1 , s2 ) is a K-function in s1 for all + is a KL-function if β(s ˆ s2 ≥ 0, β(s1 , s2 ) is strictly decreasing in s2 for all s1 > 0, and ˆ 1 , s2 ) → 0 as s2 → 0. β(s Definition 2: A set that contains the origin in its interior is a robust positively invariant set for system (1) if for all x(k) ∈ and all v(k) ∈ V it holds that x(k + 1) ∈ . Definition 3: A system is ISpS if for all x 0 ∈ , and all v ∈V x(k) ≤ ρ(x 0 , k) + (v(k − 1)) + c0 ∀k ≥ 1 where is a robust positively invariant set, ρ is a KLfunction, is a K-function, and c0 ≥ 0. Consider the nonlinear system (1) with constraints x ∈ X , u ∈ U, and v ∈ V, let L(·, ·) be a stage cost, F(·) be a terminal cost, and be a terminal set such that the following assumption [46] holds: Assumption 5: There exist positive parameters a L , a F , b F (b F ≥ a L ), sˆ ; nonnegative parameters e1 , e2 ; a function h : n → m with h(0) = 0; and a K function ϑ such that 1) ⊆ X × U and 0 ∈ 2) is a robust positively invariant set for (1) in closedloop with u = h(x) 3) L(x, u) ≤ a L xsˆ , ∀ x ∈ X and u ∈ U 4) a F xsˆ ≤ F(x) ≤ b F xsˆ + e1 , ∀x ∈ 5) F( f (x, h(x))+w(x, h(x))+v)− F(x) ≤ −L(x, h(x))+ ϑ(v), ∀x ∈ and v ∈ V.

u min ≤ u(k + j |k) ≤ u max , j = 0, 1, . . . , Nu − 1

ISpS condition [47]: The nonlinear system (1) in closedloop with the minimax MPC control obtained via (2) is ISpS in X f (N) if there exists a number bv ≥ b F such that

v min ≤ v(k + j |k) ≤ v max , j = 1, 2, . . . , N δu min ≤ δu(k + j |k)δ ≤ u max , j = 0, 1, . . . , Nu − 1

F(x) + L(x, h(x)) ≤ bv xsˆ ∀x ∈ X f (N) \

subject to

δx min ≤ δx(k + j |k)δ ≤ x max , j = 1, 2, . . . , N δymin ≤ δy(k + j |k)δ ≤ ymax , j = 1, 2, . . . , N

(5)

(6)

where X f (N) is the set of all states that can be robustly controlled into in N steps.

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 4

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

S˜ = A N ∈ n×n , V˜ = (A N−1 + · · · + I )B ∈ n×m M˜ = [(A N−1 + · · · + I )B · · · (A N−Nu + · · · + I )B]

In this paper, the terminal cost and the terminal set in the optimization problem (5) are defined as follows: F(δx(k + N|k)) = ||δx(k + N|k)||2P  

= δx(k + N|k) ∈ n | ||δx(k + N|k)||2P ≤ α

∈ n×Nu m H˜ = [ A N−1 A N−2 · · · I ] ∈ n×Nn .

(7) (8)

where P is a symmetric positive semidefinite matrix and α is a nonnegative number. The ISpS condition (6) can be fulfilled by selecting suitable P and α and a detailed discussion on the design strategy can be found in [45].

Therefore, the original optimization problem (5) becomes ˜ ||C(Sδx(k) + V δu(k − 1) + Mu(k) ¯

min max

u(k) ¯ v(k) ¯

2 ¯ + H ε¯ (k) + H v(k)) ¯ − δr¯ (k)||2Q +u(k) R ˜ ˜ u(k) + || Sδx(k) + V˜ δu(k −1) + M ¯ 2 + H˜ ε¯ (k) + H˜ v(k)|| ¯ P,

C. Minimax Problem Reformulation s.t.

Denote the following vectors:

u¯ min ≤ u(k) ¯ ≤ u¯ max

δ y¯ (k) = [δy(k + 1|k) . . . δy(k + N|k)]T ∈  N p

v¯min ≤ v(k) ¯ ≤ v¯max

δr¯ (k) = [δr (k + 1) . . . δr (k + N)]T ∈  N p δ u(k) ¯ = [δu(k|k) . . . δu(k + Nu − 1|k)]T ∈  Nu m

δ u¯ min ≤ δ u(k ¯ − 1) + I˜u(k) ¯ ≤ δ u¯ max δ x¯min ≤ δ x(k) ¯ ≤ δ x¯max

δ x(k) ¯ = [δx(k + 1|k) . . . δx(k + N|k)]T ∈  Nn u(k) ¯ = [u(k|k) . . . u(k + Nu − 1|k)]T ∈  Nu m v(k) ¯ = [v(k|k) . . . v(k + N − 1|k)]T ∈  Nn ε¯ (k) = [ε(k|k) . . . ε(k + N − 1|k)]T ∈  Nn .

δ y¯min ≤ δ y¯ (k) ≤ δ y¯max ||δx(k + N)||2P ≤ α (9)

At each time instant k, the linearization is performed once and the same linearized model is then used for predicting the future system behaviors within the entire prediction horizon. Therefore, A(k) and B(k) are written as A and B for brevity. Using (4) as the prediction model. The predicted output δ y¯ (k) is then expressed as follows: ˜ δ y¯ (k) = C(Sδx(k) + V δu(k − 1) + Mu(k) ¯ +H ε¯ (k) + H v(k)) ¯

(10)

where C˜ = block-diag[C, . . . , C]

where Q = block-diag[Q 1, Q 2 , . . . , Q N p ], R = block-diag [R1 , R2 , · · · , R Nu m ] ⎡ ⎤ I 0 ... 0 ⎢ I I . . . 0⎥ ⎢ ⎥ (13) I˜ = ⎢ . . . . ⎥ ∈  Nu m × Nu m . ⎣ .. .. . . .. ⎦ I I ... I Write u(k) ¯ and v(k) ¯ as u¯ and v¯ for brevity, the problem (12) at time k can be rewritten as the following minimax problem:  T   T  u¯ u¯ c u¯ min max W1 + 1 c2 v¯ v¯ v¯ u¯ v¯ s.t.

⎡

⎤ ⎤ ⎡ A B ⎢ A2 ⎥ ⎥ ⎢ (A + I )B ⎢ ⎥ ⎥ ⎢ S=⎢ . ⎥∈  Nn×n , V=⎢ ⎥∈  Nn×m . . . ⎣ . ⎦ ⎦ ⎣ . (A N−1 + · · · + I )B AN ⎤ ⎡ B ... 0 ⎥ ⎢ (A + I )B ... 0 ⎥ ⎢ M =⎢ ⎥ .. .. .. ⎦ ⎣ . . . (A N−1

⎡ ⎢ ⎢ H =⎢ ⎣

I A .. .

+ · · · + I )B · · · 0 I .. .

... ... .. .

(A N−Nu

⎤ 0 0⎥ ⎥ Nn×Nn . .. ⎥ ∈  .⎦

u¯ min ≤ u¯ ≤ u¯ max v¯min ≤ v¯ ≤ v¯max E 1 u¯ ≤ b1  T   T  u¯ u¯ γ u¯ 1 + 1 ≤ α˜ γ2 v¯ v¯ v¯

(14)

where W1

+ · · · + I )B)

∈  Nn×Nu m

A N−1 A N−2 . . . I

(12)

 (C˜ M)T Q C˜ M + M˜ T P M˜ + R (C˜ M)T Q C˜ H + M˜ TP H˜ = (C˜ H )T Q C˜ M + H˜ T P M˜ (C˜ H )T Q C˜ H + H˜ TP H˜ ∈ (Nu m+Nn)×(Nu m+Nn)

˜ c1 = 2(C˜ M)TQ(C(Sδx(k) + V δu(k − 1) + H ε¯ (k)) T ¯ ˜ ˜ + V˜ δu(k − 1) + H˜ ε¯ (k)) −δ (r )(k)) + 2 M P( Sδx(k) ∈  Nu m ˜ c2 = 2(C˜ H ) Q(C(Sδx(k) + V δu(k − 1) + H ε¯ (k)) T ¯ )(k)) + 2 H˜ P( Sδx(k) ˜ −δ (r + V˜ δu(k − 1) + H˜ ε¯ (k)) T

The predicted terminal state within the prediction horizon N can be expressed as ˜ ˜ u(k) δx(k + N) = Sδx(k) + V˜ δu(k − 1) + M ¯ ˜ ˜ + H (k)¯ε(k) + H v(k) ¯ (11)

E 1 = − I˜ I˜ − MM − C˜ M C˜ M

where

b1

T

∈  Nn ∈ 2(Nu m + Nn + N p)×Nu m

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. YAN AND WANG: ROBUST MPC OF NONLINEAR SYSTEMS

5

⎡

⎤ −δ u¯ min + δ u(k ¯ − 1) ⎢ ⎥ ¯ − 1) δ u¯ max − δ u(k ⎢ ⎥ ⎢−δ x¯ min + Sδx(k) + V δu(k − 1) + H ε¯ (k) + H v¯ min ⎥ ⎢ ⎥ = ⎢ δ x¯ ⎥ ⎢ max − Sδx(k) − V δu(k − 1) − H ε¯ (k) − H v¯ max ⎥ ˜ ⎣−δ y¯min + C(Sδx(k)+V δu(k −1)+ H ε¯ (k)+ H v¯ min) ⎦ ˜ δ y¯max − C(Sδx(k)−V δu(k − 1) − H ε¯ (k) − H v¯ max )

where η = [0, 0, . . . , 0, 1]T ∈  Nu m + Nn + 1  z s= ∈  Nu m + Nn + 1 θ  W1 0 W = ∈ (Nu m + Nn + 1)×(Nu m + Nn + 1) 0 0  c˜ c= ∈  Nu m + Nn + 1 −1

∈ 2Nu m+2Nn+2N p  M˜ T P M˜ M˜ T P H˜ 1 = ˜ T ˜ ˜ T ˜ ∈ (Nu m+Nn)×(Nu m+Nn) H PM H PH

E = [E 2 0] ∈ (4Nu m + 4Nn + 2N p)×(Nu m + Nn + 1)  1 0 = ∈ (Nu m + Nn + 1)×(Nu m + Nn + 1) 0 0  γ˜ γ = ∈  Nu m + Nn + 1 . (17) 0

˜ γ1 = 2 M˜ T P( Sδx(k) + V˜ δu(k − 1) + H˜ ε¯ (k)) ∈  Nu m T ˜ γ2 = 2 H˜ P( Sδx(k) + V˜ δu(k − 1) + H˜ ε¯ (k)) ∈  Nn ˜ α˜ = α − || Sδx(k) + V˜ δu(k − 1) + H˜ ε¯ (k)||2P .

Let θ = max z T W1 z + c˜ T z v¯

 T T  where z = u¯ T , v¯ T and c˜ = c1T , c2T the optimization problem (14) is then equivalent to min

III. N EURAL N ETWORK A PPROACH A. Extreme Learning Machine

θ

s.t.  T    T  z c˜ W1 0 z z + ≤ 0 θ 0 0 θ −1 θ   z  E2 0 ≤ b2 θ  T    T  z 1 0 z z γ˜ ≤ α˜ + 0 0 θ θ θ 0

The solution to the constrained minimization problem (16) gives the optimal control input increment vector u(k) ¯ that minimizes the worst case performance index at time instant k over the prediction horizon.

(15)

The ELM ( [43], [48]–[50]), is a single hidden-layer feedforward neural network. The connection weights from inputs to hidden neurons are randomly generated, whereas the connection weights from hidden neurons to outputs are computed analytically. The hidden layer of ELM does not need tuning. For N arbitrary distinct samples (xi , di ), where xi =  T xi1 , xi2 , . . . , xip ∈ p and di = [di1 , di2 , . . . , dim ]T ∈ m , ELM with L hidden nodes and activation function G(x) are mathematically modeled as L 

where ⎡

⎤

E1 0 ⎢ I 0 ⎥ ⎢ ⎥ (4Nu m + 4Nn + 2N p) × (Nu m + Nn) ⎢ E2 = ⎢ 0 I ⎥ ⎥∈ ⎣−I 0 ⎦ 0 −I ⎡ ⎤ b1 ⎢ u¯ max ⎥ ⎢ ⎥ 4Nu m + 4Nn + 2N p ⎥ b=⎢ ⎢ v¯max ⎥ ∈  ⎣−u¯ min ⎦ −v¯min  γ1 γ˜ = ∈  Nu m+Nn . γ2

i=1

L 

βi G(aiT x j + bi ) = o j , j = 1, . . . , N

i=1

T



(18)

where ai = ai1 , ai2 , . . . , aip is the weight vector connecting the i th hidden node and the input nodes, βi = [βi1 , βi2 , . . . , βim ]T is the weight vector connecting the i th hidden node and the output nodes, and bi is the threshold of the i th hidden node. aiT x j is the inner product of ai and x j . G(x) can be any infinitely differential activation function.  T   ∈ N ×m , β = β1T , . . . , β LT ∈ Let d = d1T , . . . , d NT  L×m and ⎤ ⎡ G(a1 · x1 + b1 ) · · · G(a L · x1 + b L ) ⎥ ⎢ .. .. N ×L H=⎣ ⎦∈ . ... .

Furthermore, the minimization problem (15) can be rewritten as the following linear minimization problem subjected to quadratic and linear inequality constraints: min η T s s.t.

G(a1 · xN + b1 ) · · · G(a L · xN + b L )

where d is a target vector, β is a output weight matrix, and H is called hidden layer output matrix. The learning algorithm of ELM is to find a least-squares solution β˜ of the linear system Hβ = d    ˜ (19) Hβ  = min Hβ − d . β

s T W s + cT s ≤ 0 Es ≤ b s T s + γ T s ≤ α˜

βi G i (x j ) =

(16)

The minimized norm least-squares solution of the above linear system is β˜ = H+ d (20)

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 6

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

where H+ is the Moore–Penrose generalized inverse of matrix H. Thus, ELM learning can be summarized as follows: Given a training set ℵ = {(xi , di ) | xi ∈ p , di ∈ m , i = 1, . . . , N }, hidden node activation function G and hidden node number L ELM Learn Algorithm: 1) Randomly generate hidden node weights ai and thresholds bi , i = 1, . . . , L. 2) Calculate the hidden layer output matrix H. 3) Calculate the output weight vector β: β = H+ d. B. Two-Layer RNN In [51], a RNN is presented for nonlinear convex optimization, which is capable of handling nonlinear inequality constraints. Consider the nonlinear convex program in the following form: minimize f (s) subject to c(s) ≤ 0

(21)

where s is a vector of decision  variables, f isT a real-valued vector function, and c(s) = c1 (s), . . . , cq (s) is a vectorvalued continuous function. The functions f , c1 , . . . , cq are assumed to be convex and twice differentiable. The Karush– Kuhn–Tucker conditions for (21) can be written as λ≥0 c(s) ≤ 0 ∇ f (s) + ∇c(s)λ ≥ 0 λT c(s) = 0

(22)

where λ is a dual vector for the corresponding  T Lagrangian  function, c(s)  = c1 (s), · · · , cq (s) , ∇c(s) = ∇c1 (s), · · · , ∇cq (s) is the gradient of c(s), and ∇ f (s) is the gradient of f (s). With the well-known saddle point theorem [52] and the projection theorem [53], the two-layer RNN model for solving (21) is presented in [51] with its dynamic equations as follows:     d s ∇ f (s) + ∇c(s)λ (23)  =− λ − (λ + c(s))+ dt λ T  where  is a positive constant, λ+ = (λ1 )+ , . . . , (λq )+ and [λi ]+ = max {0, λi }. According to the theoretical analysis in [51], the two-layer RNN (23) is Lyapunov stable and globally convergent to the optimal solutions to the convex program (21). C. Overall MPC Approach Rewrite (16) in a general form

As W1 and 1 in (15) are the symmetric and positive definite, W and  in (24) are the symmetric and positive semidefinite as well. Thus, c1 , c2 , and c3 are all the convex and twice differentiable. Problem (24) is a convex programming and is suitable for applying the two-layer RNN (23). However, it is worth noting that the term ε¯ (k) is still unknown so far. Therefore, the parameters c, b, γ , and α˜ in (24) are unknown. The optimization problem (24) cannot be solved unless ε¯ (k) can be estimated. ε¯ (k) represents the unmodeled dynamics and the linearization residue, and it does not depend only on the current state and input, but also the future evolution of the system and future inputs. Explicit computation of this term is demanding or even intractable. It is worth noting that due to the terminal constraint (8), the terminal state δx(k + N) must be bounded in . Thus, ε¯ (k) is bounded during the control process. Because of the universal approximation capability, ELM can be used in a black-box identification of nonlinear functions from input and output data [54]. Suppose that the state x is measurable. From a pair of input-state sample data randomly generated within the operational domain, according to (1), the state response xˆ can be measured. The errors resulting from linearization and unmodeled dynamics can be computed as ε(x, u, x o , u o ) = xˆ −

∂ f (x o ) ∂ f (u o ) (x − x o ) − (u − u o ). ∂x ∂u

An ELM can be applied for mapping the relation between the two data sets by treating [x u x o u o ]T as the input vector and ε as the target vector. A well-trained ELM from sample data sets is capable of estimating the numerical value of ε(k + j |k) j = 1, . . . , N − 1 at each time instant k. The minimax MPC for partially known nonlinear systems with bounded additive uncertainties based on the neural networks is summarized as follows. 1) Let k = 1. Set control time terminal T , prediction horizon N, control horizon Nu , sampling period τ , and weight matrices Q and R. Calculate the terminal cost matrix P and the terminal region . 2) Model ε via supervised learning using an ELM. 3) Estimate ε¯ (k) via the trained ELM in Step 2. Calculate the optimization problem coefficients, such as W , c, E, b, , and γ . 4) Solve the convex optimization problem (24) using the two-layer neural network (23) to obtain the optimal control increment vector u(k). ¯ 5) Calculate the optimal control vector u(k) ¯ and implement u(k|k). 6) If k < T , set k = k + 1, go to Step 3; otherwise end.

min f (s) = η T s IV. S IMULATION R ESULTS

s.t. c1 (s) = s W s + c s ≤ 0 T

T

c2 (s) = Es − b ≤ 0 c3 (s) = s T s + γ T s − α˜ ≤ 0.

(24)

In this section, simulation results are provided to demonstrate the performance of the proposed minimax MPC approach.

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. YAN AND WANG: ROBUST MPC OF NONLINEAR SYSTEMS

7

TABLE I ELM L EARNING R ESULTS IN E XAMPLE 1

TABLE II MLP-BP L EARNING R ESULTS IN E XAMPLE 1

Fig. 1.

L 1 -norm of training errors in Example 1.

Fig. 2.

L 1 -norm of testing errors in Example 1.

Fig. 3.

Controlled output in Example 1.

TABLE III LS-SVM L EARNING R ESULTS IN E XAMPLE 1

Example 1: Consider a nonlinear plant model presented in [55]: x˙1 = x 2 + w1 + v 1   x˙ 2 = x 12 + 0.15u 3 + 0.1 1 + x 22 u + w2 + v 2 y = x1.

(25)

Assume w = [w1 , w2 = is the unmodeled dynamics of the plant, v = [v 1 , v 2 ]T is assumed to be −0.1 ≤ v ≤ 0.1. The tracking objective is to force the output y to follow a reference signal r = sin t + cos(0.5t). The initial condition is x(0) = [0.5, 0]T . Set N = 5, Nu = 2, Q = 50I , R = 0.1I , and τ = 0.1s. The input and output constraints are −4 ≤ u ≤ 3 and −2 ≤ y ≤ 2. The strategy presented in [46] is used to compute P and α. The calculation results are P = diag(0.067, 0.148) and α = 3.4. To apply the proposed MPC approach, we first model ε using an ELM via supervised learning. The learning results based on 2000 training data and 2000 testing data are shown in Table I. To compare the performance, we also applied classical multilayer perceptron backpropagation (MLP-BP ) [56] and least-squares supporting vector machines (LS-SVM) [57] to model ε whose results are shown in Tables II and III. As ELM randomly generates hidden weights and MLP randomly generates initial weights, we run 10 times for each fixed number of hidden nodes and take the mean square errors. All experiments are performed using MATLAB toolboxes on the same computer. The L 1 -norm of training errors and testing errors corresponding to the highlighted cases are shown in Figs. 1 and 2. From the learning results, we observe that ELM ]T

[0, sin(0.1u)]T

can learn thousands of times faster and provide comparable accuracy. It is reasonable to choose ELM as the modeling tool for MPC design. The reformulated optimization problem (24) is then repeatedly solved using the RNN (23). The controlled output is shown in Fig. 3. The control input and input increments are shown in Fig. 4. To illustrate the performance, we also applied

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 8

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

TABLE V MLP-BP L EARNING R ESULTS IN E XAMPLE 2

TABLE VI LS-SVM L EARNING R ESULTS IN E XAMPLE 2

Fig. 4.

Control input in Example 1. TABLE IV ELM L EARNING R ESULTS IN E XAMPLE 2

the proposed minimax MPC on the exact model of the plant, the methods presented in [58] and [59]. The results show that the MPC approach proposed herein results in a desirable tracking performance in the presence of unmodeled dynamics and bounded uncertainties. In contrast, as the method presented in [58] does not exploit minimax strategy, it cannot constantly resist uncertainties. The method presented in [59] does not consider the higher-order term associated with linearization, inaccurate prediction model results in an undesirable tracking performance. Example 2: Consider a liquid-level system model in [60]

Fig. 5.

L 1 -norm of training errors in Example 2.

Fig. 6.

L 1 -norm of testing errors in Example 2.

y(k) = 0.9722y(k − 1) + 0.3578u(k − 1) −0.1259u(k − 2) − 0.3103y(k − 1)u(k − 1) −0.04228y 2(k − 2) + 0.1663y(k − 2)u(k − 2) −0.03259y 2(k − 1)y(k − 2) −0.3153y 2(k − 1)u(k − 2) +0.3804y(k − 1)y(k − 2)u(k − 2) +w + v(k − 2) + v(k − 1) + v(k).

(26)

Assume that w = 0.1087y(k − 2)u(k − 1)u(k − 2) is the unknown model dynamics, the bounds of the uncertainty v are ±0.1. The reference set-points are assumed to be ranging within ±0.5. Set N = 5, Nu = 2, Q = I , and R = 2I . The terminal cost ant terminal region are designed to be P = 0.8I and α = 1.2. The initial condition of [y(k − 1), y(k − 2), u(k − 2)] is assumed to be [0.1125, 0.1677, 0.6133]T . The modeling results of ε based on 2000 data using ELM is

shown in Table IV. MLP and LS-SVM are also applied for comparison and the results are shown in Tables V and VI. The training and testing errors are shown in Figs. 5 and 6. The minimax strategy is then applied to control the plant after modeling ε. The reformulated optimization problem (24) is repeatedly solved by the RNN (23) at each

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. YAN AND WANG: ROBUST MPC OF NONLINEAR SYSTEMS

9

TABLE VII ELM L EARNING R ESULTS IN E XAMPLE 3

TABLE VIII MLP-BP L EARNING R ESULTS IN E XAMPLE 3

Fig. 7.

Controlled output in Example 2.

TABLE IX LS-SVM L EARNING R ESULTS IN E XAMPLE 3

The nonlinear function α is Fig. 8.

Control input and input increment in Example 2.

sampling instant. The simulation results are shown in Figs. 7 and 8. The results show that the approach herein results in asymptotic tracking for the controlled nonlinear system in the presence of unmodeled dynamics and bounded uncertainties. Example 3: Consider the cement milling circuit model in [61] ˆ d)) + w1 + v 1 0.3 y˙ f = −y f + (1 − α(l, v s , d)ϕ(l, ˙l = −ϕ(l, ˆ d) + u f + yr + w2 + v 2 y˙r = −yr + α(l, v s , d)ϕ(l, ˆ d) + w3 + v 3 (27) where y f is the product flow rate (tons/h), l is the load in the mill (tons), yr is the tailing flow rate (tons/h), ϕ(l, d) is the output flow rate of the mill (tons/h), u f is the feed flow rate (tons/h), v s is the classifier speed (rpm), d = 1 is the hardness of the material inside the mill, and v 1 , v 2 , v 3 are the bounded uncertainties. The estimated nonlinear function ϕˆ is   ϕ(l, ˆ d) = max 0; −0.1dl 2 + 16l whereas the actual ϕ is

  ϕ(l, d) = max 0; −0.1116dl 2 + 16.5l .

α(l, v s , d) =

ϕˆ 0.8 v s4 . 3.56 × 1010 + ϕˆ 0.8 v s4

T w = [w  theTfunctional  errors.TDefine  1 , w2 , wT3 ] represent x = y f , l, yr , u = u f , v s , y = y f , yr , and v = [v 1 , v 2 , v 3 ]T . v is assumed to be −2 ≤ v ≤ 2. The reference set point at t < 3h is r = [120, 450]T , and at t ≥ 3h is r = [110, 425]T . The initial state is x(0) = [115, 60, 455]T , and the initial input is u(0) = [130, 170]T . Set N = 5, Nu = 3, Q = I , R = 0.1I and sampling period be 2 min. The constraints are: 80 ≤ u f ≤ 150 and 165 ≤ v s ≤ 180. The region are ⎤ ⎡ terminal cost and terminal 0.33 −0.09 −0.03 chosen to be P = ⎣ − 0.09 1.31 0.17 ⎦ and α = 120. − 0.03 0.17 0.81 The learning results of ε based on 2000 data using ELM is shown in Table VII. MLP and LS-SVM are also applied for comparison purpose and the results are shown in Tables VIII and IX. Training and testing errors are shown in Figs. 9 and 10. After obtaining the model of ε, the reformulated optimization problem (24) is solved by the RNN (23) at each sampling time k. The simulation results are shown in Figs. 11 and 12. The controlled outputs achieved robust tracking in the presence of bounded uncertainties.

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 10

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

Fig. 9.

L 1 -norm of training errors in Example 3. Fig. 12.

Control inputs in Example 3.

kr 2 sin x 3 + w2 + v 2 4 j1 x + w3 + v 3  4 sat(u 2 ) mˆ 2 gr kr 2 kr sin x 3 + − (l − b) + α2 j2 4 j2 2 j2 j2 kr 2 + sin x 1 + w4 + v 4 4 j2 x1 x3 (28) +

x˙3 = x˙4 =

y1 = y2 =

Fig. 10.

L 1 -norm of testing errors in Example 3.

where x 1 and x 3 are the vertical angular displacements of the pendulums, x 2 and x 4 are the angular velocities, u 1 and u 2 are the input torques generated by the servomotor, mˆ 1 = 1.8 and mˆ 2 = 2 are the estimated pendulum masses, j1 = 0.5 and j2 = 0.625 are the moments of inertia, k = 100 is the spring constant, r = 0.5 is the pendulum height, l = 0.5 is the natural length of the spring, b = 0.4 is the distance between pendulum hinges, α1 = 25 and α2 = 25 are the control input gains, and g = 9.81 is the gravitational acceleration. All the parameters are in SI unit. The actual masses are assumed to be m 1 = 2 and m 2 = 2.5. w = [w1 , w2 , w3 , w4 ]T represents the parameter variations. The function sat (·) represents the nonlinearity of actuators, which is implemented by tanh(·). The bounds of uncertainties are −0.02 ≤ v ≤ 0.02. The control objective is to track the reference r1 (t) = 0.55 cos (6.28t) r2 (t) = 0.35 cos (9.42t).

Fig. 11.

Controlled outputs in Example 3.

Example 4: Consider the double-parallel invert pendulum system model in [62] x˙1 = x 2 + w1 + v 1   sat(u 1 ) mˆ 1 gr kr 2 kr sin x 1 + x˙2 = − (l − b) + α1 j1 4 j1 2 j1 j1

(29)

The initial conditions are x(0) = [0.1, 0, −0.1, 0]T and u(0) = [0, 0]T . Let N = 5, Nu = 2, Q = 20I , R = 0.1I , and sampling periods be 0.05s. P and α are designed to be 0 for this example. We show that without enforcing terminal constraint, the controlled system can be stable by selecting a suitable Q, R, and N. The constraint is −π/2 ≤ y ≤ π/2. The modeling results of ε based on 2000 sample data using ELM is shown in Table X. The learning results of MLP and LS-SVM are shown in Tables XI and XII. Training and testing errors are shown in Figs. 13 and 14.

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. YAN AND WANG: ROBUST MPC OF NONLINEAR SYSTEMS

11

TABLE X E LM L EARNING R ESULTS IN E XAMPLE 4

TABLE XI MLP-BP L EARNING R ESULTS IN E XAMPLE 4

Fig. 14.

L 1 -norm of testing errors in Example 4.

Fig. 15.

Controlled outputs in Example 4.

Fig. 16.

Control inputs in Example 4.

TABLE XII LS-SVM L EARNING R ESULTS IN E XAMPLE 4

Fig. 13.

L 1 -norm of training errors in Example 4.

After modeling ε, the reformulated optimization problem (24) is repeatedly solved. The simulation results are shown in Figs. 15 and 16. The method herein results in robust tracking in the presence of bounded uncertainties. As shown by the four simulation examples, we observe that ELM generally trains the network thousands of times faster than MLP-BP and LS-SVM with the same number of hidden neurons/features. ELM needs more hidden neurons than MLP-BP to deliver comparable accurate outputs. The number of support vectors of LS-SVM is approximately

equal to half of the size of training data. We also find that MLP-BP is more prone to suffer over fitting problem as we increase the number of hidden neurons. But for ELM, it is capable of resisting over fitting for many cases [54]. Even if the number of hidden neurons is larger than the size of the training data, the network may still perform quite well. In addition, unlike MLP-BP and LS-SVM, the hidden layer of ELM need not be tuned. In other words, ELM is a

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 12

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS

suitable tool for modeling ε in the proposed nonlinear MPC approach. V. C ONCLUSION This paper presents a minimax MPC approach to nonlinear systems with unmodeled dynamics and bounded uncertainties based on an ELM and a two-layer RNN. The ELM is used to model the combined term of residue errors from Jacobian linearization and unmodeled dynamics. Well-trained ELM after validation is used to estimate the unknown parameters in the reformulated convex nonlinear optimization problem. A two-layer RNN is used for realizing real-time optimization. Simulation results demonstrated the superior performance of the proposed approach. R EFERENCES [1] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000. [2] J. Mu, D. Rees, and G. P. Liu, “Advanced controller design for aircraft gas turbine engines,” Control Eng. Pract., vol. 13, no. 8, pp. 1001–1015, 2005. [3] P. Tatjewski, Advanced Control of Industrial Processes, Structures and Algorithms.. New York, NY, USA: Springer-Verlag, 2003. [4] G. Colin, Y. Chamaillard, G. Bloch, and G. Corde, “Neural control of fast nonlinear systems-application to a turbocharged SI engine with VCT,” IEEE Trans. Neural Netw., vol. 18, no. 4, pp. 1101–1114, Jul. 2007. [5] M. Lawrynczuk and P. Tatjewski, “Nonlinear predictive control based on neural multi-models,” Int. J. Appl. Math. Comput. Sci., vol. 20, no. 1, pp. 7–21, 2010. [6] M. Lawrynczuk, “Explicit neural network-based nonlinear predictive control with low computational complexity,” Rough Sets Current Trends Comput., LNCS, vol. 6086, pp. 649–658, Jan. 2010. [7] L. Cavagnari, L. Magni, and R. Scattolini, “Neural network implementation of nonlinear receding horizon control,” Neural Comput. Appl., vol. 8, no. 1, pp. 86–92, 1999. [8] B. M. Akesson and H. T. Toivonen, “A neural network predictive controller,” J. Process Control, vol. 16, no. 9, pp. 937–946, 2006. [9] T. A. Johansen, “Approximate explicit receding horizon control of constrained nonlinear systems,” Automatica, vol. 40, no. 2, pp. 293–300, 2004. [10] Y. Xia and J. Wang, “A general projection neural network for solvingoptimization and related problems,” IEEE Trans. Neural Netw., vol. 15, no. 2, pp. 318–328, Mar. 2004. [11] Y. Xia, G. Feng, and J. Wang, “A primal-dual neural network for on-line resolving constrained kinematic redundancy in robot motion control,” IEEE Trans. Syst., Man, Cybern. B, Cybern., vol. 35, no. 1, pp. 54–64, Feb. 2005. [12] S. Liu and J. Wang, “A simplified dual neural network for quadratic programming with its KWTA application,” IEEE Trans. Neural Netw., vol. 17, no. 6, pp. 1500–1510, Nov. 2006. [13] X. Hu and J. Wang, “A recurrent neural network for solving a class of general variational inequalities,” IEEE Trans. Syst. Man. Cybern. B, Cybern., vol. 37, no. 3, pp. 528–539, Jun. 2007. [14] X. Hu and J. Wang, “An improved dual neural network for solving a class of quadratic programming problems and its k-winnerstake-all application,” IEEE Trans. Neural Netw., vol. 19, no. 12, pp. 2022–2031, Dec. 2008. [15] Y. Xia, G. Feng, and J. Wang, “A novel neural network for solving nonlinear optimization problems with inequality constraints,” IEEE Trans. Neural Netw., vol. 19, no. 8, pp. 1340–1353, Aug. 2008. [16] Q. Liu and J. Wang, “A one-layer recurrent neural network with a discontinuous hard-limiting activation function for quadratic programming,” IEEE Trans. Neural Netw., vol. 19, no. 4, pp. 558–570, Apr. 2008. [17] Z. Guo, Q. Liu, and J. Wang, “A one-layer recurrent neural network for pseudoconvex optimization subject to linear equality constraints,” IEEE Trans. Neural Netw., vol. 22, no. 12, pp. 1892–1900, Dec. 2011.

[18] Q. Liu and J. Wang, “A one-layer projection neural network for nonsmooth optimization subject to linear equalities and bound constraints,” IEEE Trans. Neural Netw. Learn. Syst., vol. 24, no. 5, pp. 812–824, May 2013. [19] Y. Shen and J. Wang, “Robustness analysis of global exponential stability of recurrent neural networks in the presence of time delays and random disturbances,” IEEE Trans. Neural Netw. Learn. Syst., vol. 23, no. 1, pp. 87–96, Jan. 2012. [20] P. J. Werbos, T. McAvoy, and T. Su, “Neural networks, system identification, and control in the chemical process industries,” in Handbook Intelligent Control: Neural, Fuzzy, and Adaptive Approaches, D. White and D. Sofge, Eds. New York, NY, USA: Van Nostrand Reinhold, 1992. [21] A. Draeger, S. Engell, and H. Ranke, “Model predictive control using neural networks,” IEEE Control Syst. Mag., vol. 15, no. 5, pp. 61–66, Oct. 1995. [22] G. Liu and V. Kadirkamanathan, “Predictive control for non-linear systems using neural networks,” Int. J. Control, vol. 71, no. 6, pp. 1119–1132, 1998. [23] S. Piche, B. Rodsari, D. Johnson, and M. Gerules, “Nonlinear model predictive control using neural networks,” IEEE Control Syst. Mag., vol. 20, no. 3, pp. 53–62, Feb. 2002. [24] J. Huang and F. Lewis, “Neural network predictive control for nonlinear dynamic systems with time-delay,” IEEE Trans. Neural Netw., vol. 14, no. 2, pp. 377–389, Mar. 2003. [25] Y. Song, Z. Chen, and Z. Yuan, “New chaotic PSO-base neural network predictive control for nonlinear process,” IEEE Trans. Neural Netw., vol. 18, no. 2, pp. 595–601, Mar. 2007. [26] M. Han, J. Fan, and J. Wan, “A dynamic feedforward neural network based on gaussian particle swarm optimization and its application for predictive control,” IEEE Trans. Neural Netw., vol. 22, no. 9, pp. 1457–1468, Sep. 2011. [27] L. Wang and F. Wan, “Structured neural networks for constrained model predictive control,” Automatica, vol. 37, no. 8, pp. 1235–1243, 2001. [28] L. Cheng, Z. Hou, and M. Tan, “Constrained multi-variable generalized predictive control using a dual neural network,” Neural Comput. Appl., vol. 16, no. 6, pp. 505–512, 2007. [29] Y. Pan and J. Wang, “Model predictive control of unknown nonlinear dynamical systems based on recurrent neural networks,” IEEE Trans. Ind. Electron., vol. 59, no. 8, pp. 3089–3101, Aug. 2012. [30] Z. Yan and J. Wang, “Model predictive control of nonlinear systems with unmodeled dynamics based on feedforward and recurrent neural networks,” IEEE Trans. Ind. Inf., vol. 8, no. 4, pp. 746–756, Nov. 2012. [31] Z. Yan and J. Wang, “Model predictive control for tracking of underactuated vessels based on recurrent neural networks,” IEEE J. Ocean. Eng., vol. 37, no. 4, pp. 717–726, Oct. 2012. [32] L. Chisci, J. A. Rossiter, and G. Zappa, “Systems with persistent disturbances: Predictive control with restricted constraints,” Automatica, vol. 37, no. 7, pp. 1019–1028, 2001. [33] D. L. Marruedo, T. Alamo, and E. F. Camacho, “Input-to-state stable MPC for constrained discrete-time nonlinear systems with bounded additive uncertainties,” in Proc. 41st IEEE Conf. Decision Control, vol. 4. Dec. 2002, pp. 4619–4624. [34] J. Lee and Z. Yu, “Worst-case formulation of model predictive control for systems with bounded parameters,” Automatica, vol. 33, no. 5, pp. 763–781, 1997. [35] V. Kothare, V. Balakrishnan, and M. Morari, “Robust constrained model predictive control using linear matrix inequalities,” Automatica, vol. 32, no. 10, pp. 1361–1379, 1996. [36] F. A. Cuzzola, J. C. Geromel, and M. Morari, “An improved approach for constrained robust model predictive control,” Automatica, vol. 38, no. 7, pp. 1183–1189, 2002. [37] A. Bemporad, F. Borrelli, and M. Morari, “Min-max control of constrained uncertain discrete-time linear systems,” IEEE Trans. Autom. Control, vol. 48, no. 9, pp. 1600–1606, Sep. 2003. [38] P. O. M. Scokaert and D. Q. Mayne, “Min-max feedback model predictive control for constrained linear systems,” IEEE Trans. Autom. Control, vol. 43, no. 8, pp. 1136–1142, Aug. 1998. [39] E. C. Kerrigan and J. M. Maciejowski, “Feedback min-max model predictive control using a single linear program: Robust stability and the explicit solution,” Int. J. Robust Nonlinear, vol. 14, no. 4, pp. 395–413, 2004. [40] W. Langson, I. Chryssochoos, S. V. Rakovic, and D. Q. Mayne, “Robust model predictive control using tubes,” Automatica, vol. 40, no. 1, pp. 125–133, 2004. [41] D. Q. Mayne, M. M. Seron, and S. V. Rakovic, “Robust model predictive control of constrained linear systems with bounded disturbances,” Automatica, vol. 41, no. 2, pp. 219–224, 2005.

This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. YAN AND WANG: ROBUST MPC OF NONLINEAR SYSTEMS

[42] G. Pin, D. M. Raimondo, L. Magni, and T. Parisini, “Robust model predictive control of nonlinear systems with bounded and statedependent uncertainties,” IEEE Trans. Autom. Control, vol. 54, no. 7, pp. 1681–1687, Jul. 2009. [43] G. B. Huang, Q. Zhu, and C. Siew, “Extreme learning machine: Theory and application,” Neurocomputing, vol. 70, nos. 1–3, pp. 489–501, 2005. [44] J. Garriga and M. Soroush, “Model predictive control tuning methods: A review,” Ind. Eng. Chem. Res., vol. 49, no. 8, pp. 3505–3515, 2010. [45] D. Limon, T. Alano, F. Salas, and E. F. Camacho, “Input to state stability of min-max MPC controllers for nonlinear systems with bounded uncertainties,” Automatica, vol. 42, no. 4, pp. 797–803, 2006. [46] M. Lazar, D. Pena, W. Heemels, and T. Alamo, “On input-to-state stability of min-max nonlinear model predictive control,” Syst. Control Lett., vol. 57, no. 1, pp. 39–48, 2008. [47] D. M. Raimondo, D. Limon, M. Lazar, L. Magni, and E. F. Camacho, “Min-max model predictive control of nonlinear systems: A unifying overview on stability,” Eur. J. Control, vol. 15, no. 1, pp. 5–21, 2009. [48] N. Liang, G. B. Huang, P. Saratchandran, and N. Sundrarajan, “A fast and accurate online sequential learning algorithm for feedforward networks,” IEEE Trans. Neural Netw., vol. 17, no. 6, pp. 1411–1423, Nov. 2006. [49] R. Zhang, Y. Lan, G. B. Huang, and Z. Xu, “Universal approximation of extreme learning machine with adaptive growth of hidden nodes,” IEEE Trans. Neural Netw. Learn. Syst., vol. 23, no. 2, pp. 365–371, Feb. 2012. [50] Y. Yang, Y. Wang, and X. Yuan, “Bidirectional extreme learning machine for regression problem and its learning effectiveness,” IEEE Trans. Neural Netw. Learn. Syst., vol. 23, no. 9, pp. 1498–1505, Sep. 2012. [51] Y. Xia and J. Wang, “A recurrent neural network for nonlinear convex optimization subject to nonlinear inequality constraints,” IEEE Trans. Circuit Syst. I, Regular Papers, vol. 51, no. 7, pp. 1385–1394, Jul. 2004. [52] M. Bazaraa, H. Sherali, and C. Shetty, Nonlinear Programming: Theory and Algorithms. New York, NY, USA: Wiley, 1993. [53] D. Kinderlehrer and G. Stampacchia, An Introduction to Variational Inequalities and Their Applications. San Diego, CA, USA: Academic, 1980. [54] G. B. Huang, H. Zhou, X. Ding, and R. Zhang, “Extreme learning machine for regression and multiclass classification,” IEEE Trans. Syst. Man. Cybern. B, Cybern., vol. 42, no. 2, pp. 513–529, Apr. 2012. [55] S. S. Ge, C. C. Hang, and T. Zhang, “Adaptive neural network control of nonlinear systems by state and output feedback,” IEEE Trans. Syst. Man. Cybern. B Cybern., vol. 29, no. 6, pp. 818–828, Dec. 1999. [56] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural Netw., vol. 2, no. 5, pp. 359–366, 1989. [57] J. Suykens, T. Gestel, J. Brabanter, B. DeMoor, and J. Vandewalle, Least Squares Support Vector Machines. Singapore: World Scientific, 2002. [58] Y. Pan and J. Wang, “A neurodynamic optimization approach to nonlinear model predictive control,” in Proc. IEEE Int. Conf. Syst., Man, Cybern., Oct. 2010, pp. 1597–1602. [59] Y. Pan and J. Wang, “Model predictive control for nonlinear affine systems based on the simplified dual neural network,” in Proc. IEEE Int. Symp. Intell. Control, Jul. 2009, pp. 683–688. [60] M. S. Ahmed, “Neural-net-based direct adaptive control for a class of nonlinear plants,” IEEE Trans. Autom. Control, vol. 45, no. 1, pp. 119–124, Jan. 2000. [61] L. Magni, G. Nicolao, and R. Scattolini, “Output feedback and tracking of nonlinear systems with model predictive control,” Automatica, vol. 37, no. 10, pp. 1601–1607, 2001. [62] Y. Tang, M. Tomizuka, G. Guerrero, and G. Montemayor, “Decentralized robust control of mechanical systems,” IEEE Trans. Autom. Control, vol. 45, no. 4, pp. 771–776, Apr. 2000.

13

Zheng Yan (S’11) received the B.Eng. degree in automation and computer-aided engineering from The Chinese University of Hong Kong, Shatin, Hong Kong, in 2010, where he is currently pursuing the Ph.D. degree with the Department of Mechanical and Automation Engineering. His current research interests include computational intelligence and model predictive control.

Jun Wang (S’89–M’90–SM’93–F’07) received the B.S. degree in electrical engineering and the M.S. degree in systems engineering from the Dalian University of Technology, Dalian, China, in 1982 and 1985, respectively, and the Ph.D. degree in systems engineering from Case Western Reserve University, Cleveland, OH, USA, in 1991. He is a Professor with the Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, Shatin, Hong Kong. He held various academic positions with the Dalian University of Technology, Case Western Reserve University, and University of North Dakota, Grand Forks, ND, USA. He held various short-term or parttime visiting positions with the U.S. Air Force Armstrong Laboratory in 1995; RIKEN Brain Science Institute, Wako, Japan, in 2001, Universite Catholique de Louvain, Louvain, Belgium, in 2001, Chinese Academy of Sciences, Beijing, China, in 2002, Huazhong University of Science and Technology, Wuhan, China, from 2006 to 2007, Shanghai Jiao Tong University, Shanghai, China, as a Cheung Kong Chair Professor from 2008 to 2011, and Dalian University of Technology as a National Thousand-Talent Chair Professor since 2011. His current research interests include neural networks and their applications. Dr. Wang has served as an Associate Editor of the IEEE T RANSACTIONS ON C YBERNETICS and its predecessor since 2003 and an Editorial Board Member of Neural Networks since 2012. He served as an Associate Editor of the IEEE T RANSACTIONS ON N EURAL N ETWORKS from 1999 to 2009 and the IEEE T RANSACTIONS ON S YSTEMS , M AN , AND C YBERNETICS PART C from 2002 to 2005, and an Editorial Advisory Board Member of the International Journal of Neural Systems, from 2006 to 2012. He was a Guest Editor of special issues of the European Journal of Operational Research in 1996, the International Journal of Neural Systems in 2007, and Neurocomputing in 2008. He served as the President of the Asia Pacific Neural Network Assembly in 2006, the General Chair of 13th International Conference on Neural Information Processing in 2006, and the IEEE World Congress on Computational Intelligence in 2008. He has served on many committees, such as the IEEE Fellow Committee. He was an IEEE Computational Intelligence Society Distinguished Lecturer from 2010 to 2012, the recipient of the Research Excellence Award from The Chinese University of Hong Kong from 2008 to 2009, two Natural Science Awards (first class) respectively from Shanghai Municipal Government in 2009 and the Ministry of Education of China in 2011, the Outstanding Achievement Award from Asia Pacific Neural Network Assembly, the IEEE T RANSACTIONS ON N EURAL N ETWORKS Outstanding Paper Award (with Qingshan Liu) from the IEEE Computational Intelligence Society in 2011.