Uncertainty quantification of model predictive control for nonlinear systems with parametric uncertainty using hybrid pseudo-spectral method

Abstract In this paper, a hybrid pseudo-spectral (hPS) method is utilized to analyze the uncertainty of the model predictive control (MPC) for nonlinear systems with stochastic parameter uncertainty. The hPS method incorporates Generalized polynomial expansion (gPC) and pseudo-spectral (PS) optimal control in a hybrid format. To quantify the effect of uncertainty on the MPC states and control inputs, we use the gPC method, and for time discretization of the MPC problem, the pseudo-spectral method is employed. The proposed method will be compared with the well-known Monte Carlo (MC) simulation method in terms of computational speed. Two dynamical systems are tested. First, an industrial process, continuous stirred reactor tank, which is a highly nonlinear system, with Gaussian distributed parameter uncertainty and second, a typical second-order dynamical system with uniform and Weibull distributed parameter uncertainty. Simulation results demonstrate that the proposed method can approximate the two first moments of the states and control inputs under MPC much faster and computationally more efficient than MC simulations.


Introduction
MPC is one of the most practical methods of advanced process control in the industry. It is an optimal control strategy that is implemented in a receding horizon manner. By implementing the MPC strategy, Ali Namadchian ABOUT THE AUTHOR Ali Namadchian received the M.Sc. degree in control theory and control engineering from the Islamic Azad University of Mashhad, Mashhad, Iran, in 2013; he is currently pursuing the Ph.D. degree with the Department of Electrical and Control Engineering, Tafresh University, Tafresh, Iran. His current research interests include stochastic control, nonlinear control systems and optimal control.

PUBLIC INTEREST STATEMENT
Uncertainty quantification is the science of quantitative characterization of uncertainties in real-world applications. It tries to determine how likely certain outcomes are if some aspects of the system are not exactly known. Most practical systems are exposed to uncertainty and engineers need to know how the uncertainty will affect the system. In this paper, we proposed a hybrid pseudo-spectral algorithm to quantify the effect of uncertainty on nonlinear systems under a model predictive control strategy.
certain constraints can be imposed on outputs, inputs, and states of the system. Furthermore, it can be implemented on both linear and nonlinear dynamical systems. Many applications of linear MPC have been reported in the literature (Garcia, Prett, & Morari, 1989;Muske & Rawlings, 1993;Schwickart, Voos, Hadji-Minaglou, & Darouach, 2016;Widd, Liao, Gerdes, Tunestål, & Johansson, 2014). For some highly nonlinear systems, the approximated linear model is not sufficient. Therefore, the MPC should be designed directly based on the nonlinear model (Grüne & Pannek, 2011;Mayne, Rawlings, Rao, & Scokaert, 2000). Several theories have been proposed regarding the stability and feasibility of the MPC, both in linear and nonlinear context. In nonlinear MPC (NMPC), the system or constraints are nonlinear. The efforts to prove the stability of NMPC have led to the well-known four axioms of stability (Mayne et al., 2000). Such as nonlinearity, which is the inherent characteristic of real-world systems, parameter uncertainties occur in many practical systems. Uncertainties in the parameters of a system could either be bounded or have a stochastic nature.
When the uncertainty is bounded, the well-known robust MPC strategies are applied. Tube-based MPC is a dominant approach in this field (Mayne, Kerrigan, Van Wyk, & Falugi, 2011). For systems with stochastic parameter uncertainty, where the probability distribution of uncertainty is known, the stochastic MPC strategies are employed. In Visintini, Glover, Lygeros, and Maciejowski (2006), the Markov-chain Monte Carlo technique was proposed for solving constrained nonlinear stochastic optimization problems. In Bavdekar and Mesbah (2016), at each sampling time of the MPC, the gPC coefficients were calculated and then the MPC was performed on the gPC coefficients. The Fokker-Planck equation, which described the probability density of the states of the system, was utilized in Annunziato and Borzì (2013) to solve the MPC for the stochastic differential equation. In Annunziato and Borzì (2013), the MPC for the system of stochastic differential equations was transformed into a deterministic MPC with PDE constraint. The stochastic Lyapunov function was considered as an auxiliary constraint in the stochastic MPC framework in Mahmood and Mhaskar (2012) to guarantee the stochastic stability of the SDE under MPC strategy.
Control engineers should always have a general sight of the effects of uncertainty on the system under the MPC strategy. By analyzing the uncertainty effects on the system, one can choose the system hardware and software accordingly. For instance, the effect of the uncertainty on the control actuation can be evaluated and the size of the actuator can be determined. Uncertainty quantification is the science of quantitative characterization of uncertainties in real-world applications. It tries to determine how likely certain outcomes are if some aspects of the system are not exactly known. This paper is the continuation of our previous work (Namadchian & Ramezani, 2017) in which we used the pseudospectral method for time discretization of the MPC for a CSTR and showed the computational efficiency of the proposed method compared with RK4 time discretization method. In this paper, we have not only used the PS method for time discretization but also applied gPC, for uncertainty quantification and propose a hybrid pseudo-spectral method for uncertainty quantification of nonlinear MPC.
Pseudo-spectral methods referred to a class of spectral techniques which are solved by collocation points (Fornberg, 1998). In this paper, for uncertainty propagation in nonlinear MPC, the Generalized Polynomial Chaos (Xiu & Karniadakis, 2002) has been utilized. GPC is a very efficient spectral method for uncertainty evaluation. In the sample-based version of the gPC or collocation-based gPC, the uncertainty can be analyzed with much fewer samples than Monte Carlo (MC) simulations.
The proposed hybrid algorithm selects samples from the random space using collocation nodes, after that, inserts those sample values into the differential equations of the system, and then uses a pseudo-spectral based deterministic solver to generate solutions for each of the resulting deterministic problems.
The set of deterministic solutions is then used to construct a polynomial representation of the solution of the stochastic optimal control problem using the gPC method. Since both gPC and PS optimal control can be categorized as PS methods and are incorporated in a hybrid format, we call this approach, the hybrid pseudo-spectral or the hPS method.
PS optimal control is an efficient computational method for solving optimal control problems. The most advanced optimal control softwares such as GPOPS, PSOPT, and RPROPT use PS optimal control to discretize the nonlinear optimal control problem and transform it into a nonlinear programming problem (Becerra, 2010;Patterson & Rao, 2014. The hybrid method, involving a combination of the PS optimal control and gPC, has been appeared in few publications. Harmon (2017) utilized the PS optimal control with gPC to quantify the effect of uncertainty on the optimal states of the system. Matsuno, Tsuchiya, Wei, Hwang, and Matayoshi (2015) solved the aircraft conflict resolution by the hPS. The system was a 3D aircraft dynamic under wind uncertainty. The main goal of the optimal control problem was to determine the optimal trajectories in a way that the confliction of the aircraft is avoided. The conflict was modeled as the distance between two aircraft and as a constraint in the optimal control problem, the distance must not exceed a predefined threshold value. Motivated by his previous work, Matsuno developed his approach for aircraft conflict resolution for online applications . By Polynomial Chaos Kriging method, he constructed a surrogate model for the solution of stochastic optimal control method based on the hPS. Then, the surrogate model was used to circumvent the need for solving a stochastic optimal control problem on-line. Although the mentioned papers utilized both PS optimal control and gPC in a hybrid manner, they did not address the optimal control problem in a receding horizon approach, i.e. MPC framework. What distinguishes our work from these papers is considering the uncertainty quantification of optimal control in a receding horizon manner. The main contribution of our work is to quantify the effect of parameter uncertainty on the states and controls of the nonlinear system under MPC strategy with less computational effort than MC simulations.
The paper is structured as follows: Section 2 gives details on the problem formulation and notations. The PS fundamentals are described in Section 3. The polynomial chaos theory is introduced in Section 4. The hybrid PS-MPC uncertainty propagation algorithm is presented in Section 5. In Section 6, uncertainty evaluation is performed on two examples under MPC strategy and compared with MC simulations.

Problem formulation
Consider the following continuous-time system in the state space form: Since in MPC literature, the control is normally piecewise constant, without loss of generality, we can always express the system (1) with a difference equation by an appropriate sampling time Δt could be expressed as follows: where xðkÞ 2 R nx is the state of the system at time index k, uðkÞ 2 R nu is the manipulated input, θ 2 R nθ is the vector of unknown parameters. θ is stochastic uncertainty. θ i , the component of θ are assumed to be independent and identically distributed random variables with pdf f θ i such that θ i 2 K & L 2 ðΩ; F; PÞ, ðΩ; F; PÞ is the probability space, Ω is the sample space, F is the σ À a lg ebra of the events, and P is the probability measure of events. The first moment or expectation of a random variable θ is defined as EðθÞ ¼ ð Ω θðwÞdPðwÞ. L 2 ðΩ; F; PÞ is the Hilbert space of all random variables θ whose L 2 À norm is finite, i.e. jjθjj 2 ¼ E½jθj 2 1=2 <1. Also, the second moment of θ, called variance, is defined as VarðθÞ ¼ E½ðθ À E½θÞ 2 ¼ σ 2 θ , where σ θ is the standard deviation. We assume that f is continuously differentiable with respect to its arguments and the solution of (2) exists uniquely and almost surely.
The MPC formulation for Equation (2) is as follows: subject to : xðk þ 1jtÞ ¼ f ðxðkjtÞ; uðkjtÞ; θÞ (4) xðkjtÞ 2 X k ¼ 1; :::; N p (5) uðkjtÞ 2 U k ¼ 0; :::; N c (6) xðNjtÞ 2 X f (7) xð0jtÞ ¼ xðtÞ (8) where xðkjtÞ is the predicted value of x, k step ahead of t, based on the information available at time t. N p and N c are prediction and control horizons, respectively. X and U are state and control constraints, respectively, and X f is the terminal constraint. V N ðx; uÞ, the cost function, consists of the stage cost lðx; uÞ and the terminal cost V f ðxÞ.
It is assumed that at each step, the states of the system xðkÞ are measurable. Due to the randomness of θ, xðkÞ and uðkÞ are stochastic processes. Our goal is to find the mean and variance of the states and control inputs at each time index k to find how the randomness of θ will affect the distribution of control inputs and the states of the uncertain system under the MPC strategy.
The terminal constraint (7) imposes an implicit constraint on the control sequence of the form: where the control constraint set U N ðx; θÞ is the set of control sequence u :¼ fuð0Þ; uð1Þ; :::; uðN À 1Þg satisfying the state and control constraints: ¼ fðx; uÞjuðkÞ 2 U; xðkÞ 2 X; " k 2 0; 1; . . . ; N À 1 f g and xðNÞ 2 X f ðθÞ g Therefore, the optimal control problem P N ðx; θÞ is defined as: The following assumptions are considered for the MPC controller: Assumption 1.
(a) Continuity of the system: for each θ i 2 K & L 2 ðΩ; F; PÞ, the function f : and lð:; :Þ is strictly positive function.
(b) Properties of the constraint set: for each θ i 2 K & L 2 ðΩ; F; PÞ, the set X and X f ðθ i Þ are closed, X f ðθ i Þ 2 X and U are compact, each set contains the origin.
Theorem.1 (axioms of stability for MPC) (Mayne et al., 2000): If for any θ i 2 K & L 2 ðΩ; F; PÞ, the following assumptions are satisfied, then the states of the system f ðx; u; θ i Þ converge asymptotically to the origin under MPC formulation P N ðx; θ i Þ: A1: State constraint satisfied in X f ðθ i Þ: There exists a local controller κ f ðx; Proof: In the axioms of stability for the MPC, we consider X f ðθ i Þ, V f ðx; θ i Þ, κ f ðx; θ i Þ such that the axioms are satisfied for each θ i , so without loss of generality, to prove the theorem, we consider θ i as a constant and omit it from the equations.
The stability of the MPC will be proved based on using the value function as a decreasing Lyapunov like function. In the first part, we prove that initial feasibility implies recursive feasibility and in the second part, the convergence of the system to the origin will be guaranteed under MPC strategy.

Part2-Asymptotic stability
The most popular idea in the proof of the stability of the MPC is to show that optimal objective For the optimal control sequence u ¼ fu Ã ð0Þ; u Ã ð1Þ; :::; u Ã ðN À 1Þg calculated at xðtÞ, the optimal objective function is: Also, for sub-optimal control sequenceũ ¼ fu Ã ð1Þ; :::; u Ã ðN À 1Þ; κ f ðx Ã ðNÞÞg the suboptimal objective function is as follows: is a sub-optimal solution, we have: can be rewritten as follows: By adding and subtracting V f ðx Ã ðNÞÞ to the above expression, we obtain: It can be seen that the first two terms in the above equation is equal to V Ã N ðxÞ, thus we can rewrite it as: From A4 (MPC axioms of stability), for the last tree term in the above expression we have: Leaving us with the expression: Since it is assumed that lð0; 0Þ ¼ 0 and lð:Þ is a strictly positive function, we have: Indicating that the objective function V Ã N ðxðtÞ is strictly decreasing along closed-loop trajectories under MPC control, thus the trajectories converge to the origin. To design an MPC for the nonlinear system (2), we choose lðx; uÞ ¼ ð1=2Þðjxj 2 Q þ juj 2 R Þ, in which Q and R are positive definite. In order to select a proper V f ðx; θ i Þ and X f ðθ i Þ, the nonlinear system (1) is linearized at the origin to obtain a linear model terminal cost function is chosen to be: And X f ðθ i Þ is the sublevel set WðaÞ ¼ fxjV f ðx; θ i Þ ag for some suitably chosen constant a. Choosing such a controller, it can be proved that there exists an a ! 0 such that. V f ðf ðx; Kðθ i ÞxÞÞ þ lðx; Kðθ i ÞxÞ À V f ðx; θ i Þ 0; "x 2 X f :¼ WðaÞ In order to satisfy the state and control constraints, we choose a such that X f ðθ i Þ X and KX f ðθ i Þ U.
It can be demonstrated that by the designed controller, the closed-loop system under the MPC strategy is asymptotically stable (Rawlings & Mayne, 2012).
Remark 1: Most of the system variables are functions of θ. In fact, the whole MPC problem P N ðx; θÞ, the linearized system AðθÞ and BðθÞ, the local controller k f ðx; θÞ, the designed terminal constraint X f ðθÞ and also the states xðk; θÞ and uðk; θÞ are all functions of the random variable θ. The hard constraints X and U are fixed sets and are not the functions of θ.
Remark 2: This paper investigates that how the system variables must have changed under parameter uncertainty to have a stable system under the model predictive controller. Although we can quantify the effect of uncertainty on all of the parameters mentioned in Remark1, in this paper, we just focus on the uncertainty quantification of the states and controls input of the nonlinear system with parameter uncertainty under MPC.

Pseudo-spectral fundamentals
There are two approaches to solve a continuous-time optimal control problem, direct methods, and indirect methods (Von Stryk & Bulirsch, 1992). Indirect methods use the calculus of variation such as the Pontryagin minimum principle. Indirect methods are usually difficult to solve numerically. Direct methods are known as discretize-then-optimize methods. In direct methods, the state and control functions are approximated by polynomials.
When an optimal control problem solved by the direct method, having fewer discretization points is an advantage. A typical optimal control problem consists of three mathematical equations: (1) integration in the cost function, (2) differential equation governing the system state evolution, and (3) algebraic equations appear as the constraints. PS method is suitable for all three objects. There are a lot of theoretical proofs regarding the convergence, consistency, and optimality of the solution of PS optimal control (Gong, Ross, Kang, & Fahroo, 2008;Ross, 2005).
The PS approximation was initially developed for solving PDEs. In this method, continuous functions are approximated in a finite set of nodes called collocation nodes. These nodes are selected in a way that guaranty the high accuracy of the approximation. The main exclusivity of the PS method is the smart selection of collocation nodes. Consider the approximation of f ðtÞ on ½À1; 1 by some interpolants. The simplest choices of nodes are equal distance nodes, but if a special set of nodes such as roots of derivative of Legendre polynomial, known as Legendre-Gauss-Lobatto (LGL) nodes is employed, the higher accuracy with fewer nodes can be achieved. This inequality shows the impressive convergence speed of the PS method (Gong et al., 2007).
where I N f ðtÞ is the polynomial interpolation of f ðtÞ. C is a constant. N is the number of collocation points and m is the order of differentiability of f ðtÞ. It can be seen that if f ðtÞ 2 C 1 ðm ¼ 1Þ, the polynomial converges at a spectral rate.
In this paper, we chose the Gauss-PS method with Lagrange interpolant. Since the prediction horizon of optimal control problem in the MPC framework is finite, the LGL nodes are adapted for discretization.
Given t ¼ ½τ 0 ; τ 1 ; :::; τ N set of distinct points in ½À1; 1, a standard Lagrange interpolant for control function uðtÞ and xðtÞ can be written as: where u N ðtÞ and x N ðtÞ are a unique polynomial of degree N and L N i ðτÞ is the N-th order Lagrange polynomial which satisfies the Kronecker property L N i ðτ j Þ ¼ δ ij . This implies that uðτ j Þ ¼ u N ðτ j Þ and xðτ j Þ ¼ x N ðτ j Þ.
The first derivative of f ðxÞ can be written as where _ L N i ðtÞ is expressed as: From these two equations, the derivative of f ðxÞ at the LGL node τ k 2 ½À1; 1 is approximated by (14) can be represented in the matrix form: where D Nþ1 is a ðN þ 1Þ Â ðN þ 1Þ matrix defined by: It is obvious that the approximation of _ f ðtÞ only depends on f ðtÞ. For nodes that are not scaled, i.e. τ k 2 ½a; b, the derivative of f ðtÞ at the LGL nodes is as follows: In the PS method, the integration is approximated by the Gaussian quadrature rule. The rule is as follows: where in Gauss-PS method τ i ; i ¼ 1; :::; N are LGL nodes and w i ; i ¼ 1; :::; N are predetermined quadrature weight correspond to LGL nodes. For the arbitrary domain of integration over ½a; b, the Gauss quadrature rule can be rewritten as

Polynomial chaos uncertainty propagation
gPC is a method for approximation of random variables with finite second-order moments (Cameron & Martin, 1947;Wiener, 1938). The random variable ψ θ ð Þ 2 L 2 Ω; F; P ð Þ has the L 2 À convergent expansion (Xiu & Karniadakis, 2002): where a k are expansion coefficient, Φ a k θ ¼ Q n i¼1 Φ α i;k θ i are k À th order multivariate polynomials and Φ a i;k θ i ð Þ are univariate polynomials of degree a i;k . The optimal form of Φ a i;k θ i is based on the probability distribution of θ i and is chosen according to the Wiener-Askey scheme (Table 1). For instance, for Gaussian distributed θ i , the Hermit polynomials are the best choice. These polynomials satisfy the orthogonality property For practical reason, (22) is truncated, the L term of the series is considered, where L depends on the number of random variables. n and m are selected such So the random variable expansion can be reformulated as: The coefficients a i are calculated as follows: To calculate (27), the collocation method is adapted. A proper choice for collocation points is the roots of ðm þ 1Þ À th order polynomial that is chosen based on the Askey scheme. Employing orthogonality property, the first and second moments of the random variable ψ θ can be derived by gPC coefficients as follows:

Hybrid pseudo-spectral MPC uncertainty propagation
The states and inputs of the system with parameter uncertainty are stochastic processes. xðt; θÞ and uðt; θÞ are the functions of random variable as well as the function of time.
In Section 3, it has been shown that a function of time can be approximated by a polynomial. Just like PS optimal control method which approximates xðtÞ and uðtÞ, with polynomials, gPC approximates xðθÞ and uðθÞ with the same approach. At each sample time t 0 , xðt 0 ; θÞ ¼ xðθÞ and uðt 0 ; θÞ ¼ uðθÞ are random variables. The objective of gPC is to approximate these random variables at each sample time. The main difference between PS optimal control and gPC is that in the PS optimal control, the approximation is performed over the time domain, whereas in gPC the approximation is carried out over the random variable probability space. Therefore, the main goal is to approximate the distribution of xðt; θÞ and uðt; θÞ.
The procedure of the hybrid algorithm is as follows: (1) Initial parameters: • Pick fq i ; α i g N i¼1 , N collocation points and quadrature weights based on the distribution of each random variable and Wiener-Askey scheme. Create the tensor grid of this point where M is the number of random variables. Find the values of Φ a k p i ð Þ in (22).
• Choose fτ i ; w i g, the LGL collocation points and weights for PS optimal control.
• Choose overall simulation time T, MPC parameters such as prediction horizon, control horizon.
(2) For each p i , design X f ðp i Þ, V f ðx; p i Þ, κ f ðx; p i Þ according to Section 2.1 and solve the MPC problem by PS optimal control method: In this step, by utilizing PS optimal control, the problem 1 is transformed into the following nonlinear programming (NLP): w k lðx N ðτ k Þ; u N ðτ k ÞÞ s:t: In this paper, the sequential quadratic programming is employed to solve the above optimization problem. At this step, x N ðp i ; τ i Þ and u N ðp i ; τ i Þ are calculated during simulation time T at both p i and τ i collocation points.
(3) Calculate the gPC expansion coefficients For the values of x N ð:; τ i Þ and u N ð:; τ i Þ at each collocation points τ i , calculate the coefficients of the gPC using quadrature formula for Equation (27) (4) Derive the statistics of input and state of the system by Equations (28) and (29).

Numerical examples
In this section, the uncertainty propagation for two dynamical systems under MPC control is analyzed. The first one is an industrial process, continuous stirred reactor tank (CSTR). The CSTR has exposed to Gaussian distributed parameter uncertainty. The second one is a second-order dynamic system with uniform and Weibull distributed parameter uncertainty. All the simulations are performed on a Core i7 laptop with 6G Ram.
Example 1. CSTR as a highly nonlinear second-order dynamical system is a benchmark process for advanced control strategies. In a CSTR an irreversible first-order exothermic reaction of the form A ! B is taking place. The dynamics of a CSTR tank is as follows (Mahmood & Mhaskar, 2012): where C A is the concentration of species A, T and V R are the temperature and volume of the reactor, respectively. Q is the amount of heat input to the reactor. k 0 , E, and ΔH are the preexponential constant, the activation energy, and enthalpy of the reaction, respectively. c p and ρ denote heat capacity and density of the fluid in the reactor. The states of the system are x 1 ¼ C A and x 2 ¼ T. The manipulated variables are u 1 ¼ C A0 À C A0s and u 2 ¼ Q. The MPC aims to reach to the desired steady-state values of C As and T Rs . The uncertainty in the process is considered due to the uncertainty in the two physical parameters k 0 and ΔH. This uncertainty is inevitable because of the experimental measurement complexity of thermal and kinetic phenomena. The values of all the parameters are listed in Table 2.
As it can be seen from Table 2, the distributions of the uncertain parameters ðk 0 ; ΔHÞ are Gaussian, so for gPC expansion, the Hermite polynomial is selected based on the Wiener-Askey scheme. The nominal values of uncertain parameters are k 0 ¼ 72 Â 10 9 and ΔH ¼ À4:78 Â 10 4 . Our goal is to steer the state of the system to the equilibrium point x ss ¼ ðC Rs ; T Rs Þ. Without loss of generality, we shifted the equilibrium point to the origin.
The cost function is lðx; uÞ ¼ ð1=2Þðjxj 2 Q þ juj 2 R Þ with Q ¼ ½1 0; 0 1 and R ¼ ½0:1 0; 0 0:1. The constraints of the control are ju 1 j 1 kmol=m 3 and ju 2 j 92 kj=s. The state constraints are jC A j 1 kmol=m 3 and À 100 K < T A < 100 K. The initial condition is x ¼ ½À0:27; 5. We use Legendre polynomial for PS optimal control and Hermite polynomial for gPC uncertainty propagation. The prediction and control horizons are N p ¼ 20s and N c ¼ 0:25s, respectively. The system is simulated for 3 s. Also, 10 PS optimal control collocation points and 6 Â 6 ¼ 36 gPC collocation points are considered. As discussed in Section 2.1, (30) was linearized around the equilibrium point to find X f ðθ i Þ, V f ðx; θ i Þ and κ f ðx; θ i Þ. For the terminal constraint we choose X f ðθ i Þ ¼ fxjV f ðx; θ i Þ 10g and for the linearized feedback controller we choose κ f ðx; θ i Þ such that the closed-loop linearized system A K ðθ i Þ has two stable poles at −15 and −5 for each θ i . Figure 1 shows how the parameter uncertainty affects the terminal constraint X f . For different gPC collocation points, the terminal constraints are calculated and sketched in Figure 1. Figure 2 shows how the uncertainty in the parameters of the CSTR affects the optimal control inputs of the system. In other words, it demonstrates how the control input should change to cope with the uncertainty of the parameters in order to have an optimal controller under stable MPC. When the CSTR is exposed to uncertainty, the controls and states will be a stochastic process; thus in Figure 2, the control first and second moments are indicated. The first moment is the mean of the control sequence which is represented by red lines and the gray interval illustrates the threesigma area, i.e. μ AE 3σ which is a representation of the second moment. We choose to show the Table 2. Parameters of CSTR Based on the three-sigma rule for a Gaussian distributed random variable, the probability of values generated by Gaussian pdf lying in the three-sigma area is 0.9974. We cannot claim that for a nonlinear system, with Gaussian uncertainty parameter distribution, the distribution of the states and the controls of the system remains Gaussian for all time. Generally, with the three-sigma rule, we cannot identify the probability of system inputs but it can visualize the two first moments of the system inputs. Also, by the two first moments, we can apply Markov inequality to characterize the probability of the system states and controls in some bounded domain. Fortunately, in contrast to the three-sigma rule, the Markov inequality is not limited to just Gaussian distribution, it can be used for a wide group of the distributions. To highlight the numerical validation of the proposed method, it is compared with MC simulations. In order to be sure that MC simulations are convergent for our system, we already checked the convergence of MC simulations for both mean and variance of the states and the control of the system with parameter uncertainty under MPC. Figures 3 and 4 demonstrate the comparison between mean and the variance of the states of the system for both MC and gPC uncertainty propagation. Figures 3 and 4 show that the propose method results are consistent with MC simulations. Also, from Figures 3 and 4, it can be seen that the MPC objective, i.e. stabilization, is satisfied and the states converge to the origin.
The advantage of gPC method compared to MC simulations is the calculation of the mean and the variance of the system with much fewer samples and less computational time. Table 3 compares the computational time and complexity of the gPC and MC simulations for the CSTR tank. Figure 2. Effect of the parameter uncertainty on the control input for CSTR optimal control problem. Figure 3. Mean of the optimal states of the CSTR with MC and gPC uncertainty propagation. Table 3 shows that in order to find the effect of the parameter uncertainty for model predictive control of a CSTR system, 30,000 MC simulations must be performed to have a convergent first and second moment for the states and controls of the system. On the other hand, when gPC is used, only 36 samples are needed to determine the first and second moments of the states and the controls of the CSTR tank. The MPC computational time takes 35 h, whereas gPC method only took 11 min of computational time. Thus, it can be found that gPC is a much more efficient method for uncertainty propagation than the MC simulation method.

Example 2.
In Example 1, the validity of the hPS method has been confirmed for uncertainty propagation when the uncertain parameters both have Gaussian distributions. In this example, the ability of the hPS method will be demonstrated when the uncertain parameters of the system have uniform distributions. Consider the following second-order system with parameter uncertainties A and B: _ x 1 ðtÞ ¼ x 2 ðtÞ þ Bx 2 ðtÞ _ x 2 ðtÞ ¼ Àðx 1 ðtÞ þ Ax 1 ðtÞÞ þ uðtÞ þ ½1 À ðx 1 ðtÞ þ Ax 1 ðtÞÞ 2 ððx 2 ðtÞ þ Bx 2 ðtÞÞ The nonlinear dynamic is under MPC strategy with the cost function lðx; uÞ ¼ ð1=2Þðjxj 2 Q þ juj 2 R Þ with Q ¼ ½1000 0; 0 1000 and R ¼ ½0:0001 0; 0 0:0001, state constraint À 100 < x 1 ; x 2 < 100, control constraint À 100 < u < 100 and the initial condition x ¼ ½3; 3. The goal of the MPC is to stabilize the system states to the origin.
The parameter A has uniform distribution U½À0:01; 0:01 and parameter B has Weibull distribution WeðBÞ ¼ αβB βÀ1 e ÀαB β with α ¼ :1 and β ¼ 1. According to Wiener-Askey scheme for the uniform distribution, the Legendre polynomial should be selected, but no type of polynomial is determined for Weibull distribution. To use Legendre polynomial, fundamental transformation law of probabilities is utilized to map a random variable with uniform distribution to the B as follows: WeðBÞ ¼ ðÀ1=αÞ lnð1 À bÞ ½ 1 β ; b ¼ U½0; 1  Therefore, in the dynamics of the system in the place of B, the term ðÀ1=αÞ lnð1 À bÞ ½ 1 β should be inserted. We transfer the Weibull distribution to a uniform random distribution to be able to use a type of polynomial according to standard distribution in Wiener-Askey scheme (Table 1). For the terminal constraint we choose X f ðθ i Þ ¼ fxjV f ðx; θ i Þ 10g and for the linearized feedback controller, κ f ðx; θ i Þ is chosen such that the closed-loop linearized system A K ðθ i Þ has two stable poles at −2 and −1 for each θ i . For this example, Legendre polynomial for both PS optimal control and gPC uncertainty propagation is applied. The prediction horizon is N p ¼ :6 s, control horizon is N c ¼ 0:1 s, and the system is simulated for 11 s. Also, the number of PS collocation points in the time direction considered to be 4 and 20 Â 20 ¼ 400 gPC collocation points are considered. Figure  5 highlights how the parameter uncertainty affects the optimal control solution of the system. Figures 6 and 7 indicate the consistency of the gPC method with MC simulations in MPC uncertainty propagation. Also, the MPC objective, i.e. stabilization, is satisfied. It shows that gPC uncertainty propagation for Example 2 is much more efficient than MC simulations. Figures 6 and 7 show that gPC results are in consistent with MC simulation results. Figure 5. Effect of the parameter uncertainty on the control for optimal control problem of Example 2. Figure 6. Mean of the optimal states of Example 2 with MC and gPC uncertainty propagation.

Conclusion
In this paper, an efficient method is proposed to evaluate the effect of parameter uncertainty on the inputs and the states of a nonlinear dynamical system under the MPC strategy. The presented method implements a hybrid pseudo-spectral method to increase the speed of the process of uncertainty propagation. Based on the proposed method, the first and second moments of the control sequence of a nonlinear system with uncertain parameters have been found. Two simulation results validate the effectiveness of the proposed method. It is shown that the gPC results are consistent with MC simulations, but gPC method is much superior to MC simulations regarding the computational time. In this paper, the hybrid PS method evaluates the first and second moments of the distributions. For future works, the other moments of the distribution can be evaluated by gPC coefficients to have a better analysis of the non-Gaussian distributions. Figure 7. Variance of the optimal states of Example 2 with MC and gPC uncertainty propagation.