State space model predictive control for advanced process operation

Article pubs.acs.org/IECR

State Space Model Predictive Control for Advanced Process Operation: A Review of Recent Development, New Results, and Insight Ridong Zhang,* Sheng Wu, and Furong Gao Department of Chemical and Biomolecular Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong ABSTRACT: Model predictive control (MPC) has acquired lots of developments and extensive applications in various industries during the past 40 years. For the early version of basic MPC strategies, there are challenges in both control performance and applications. On the basis of such backgrounds, state space models that contain more system information on the controlled system than conventional models and make the controller design straightforward under MIMO situations were introduced to the MPC framework, and many academic results have been put forward. One of the challenging missions for current MPC research is still focused on the optimal control of systems with improved (or robust) control performance. Therefore, it is of significance to investigate new developments of MPC. In order to improve the control performance of MPC, we introduced a new comprehensive state space model formulation strategy in which state variables and the output tracking errors are combined and can be regulated separately. On the basis of this model formulation, more degrees of freedom can be offered during the design of controllers for improved control performance. The corresponding MPC algorithms have been introduced to deal with singleinput single-output (SISO)/multiple-input multiple-output (MIMO) systems under process/model uncertainties to optimize PID controllers and to be implemented in typical industrial coke equipment in our recent work. These results demonstrate the effectiveness of the proposed models-based MPC algorithms. A brief review of the recent developments of state space MPC strategies, a comprehensive summary of our recent work, and further insights based on a comparison with the typical work by J. B. Rawlings in the literature (Pannocchia, G., Rawlings, J. B. Disturbance models for offset-free model-predictive control. AIChE J., 2003 49 (2), 426−437) are presented in this paper.

1. INTRODUCTION For industrial processes with various constraints, the applications of the classic control strategies, such as PID control,1−4 may not always achieve the ideal control performance. In order to handle such processes that exist widely in practice, MPC was put forward as a promising advanced control algorithm and had acquired extensive applications.5−7 As is known to all, the early basic MPC strategies, i.e., dynamic matrix control (DMC), model algorithm control (MAC), and generalized predictive control (GPC), encounter the challenges of both performance and applications under the increasing demands. How to enhance their control performance, extend their applications, reduce their computation time, etc. are all topics that need paid attention to. State space models that contain more system information promote the development of MPC schemes. On the basis of such models, the controller design for industrial processes is straightforward, and relevant system analysis is simplified. There are numerous fruits of state space MPC both on theory and applications during the past decades.8 Complex high-dimensional processes exist widely in various industries, and there are always dozens or even hundreds of variables and constraints. In order to solve the control issue of such processes online, the requirements on the computers in © 2017 American Chemical Society

terms of performance and computation speed are even stricter. Focusing on this problem, distributed MPC obtains a lot of attention. A large-scale online optimization problem can be decomposed into many small-scale distributed optimizations of intelligent agents by distributed MPC strategies, then the requirement on the computer performance will be reduced. At the same time, the performance reduction caused by information loss of the conventional decentralized MPC will be overcome in the distributed MPC scheme because there are information communications between these intelligent agents. Many researchers presented their significant results. On the basis of a cooperative game, a distributed MPC strategy was developed for the control of two constrained linear systems.9 In ref 10, the tutorial review of the recent developments of the distributed MPC scheme was presented. Liu et al.11 handled a class of nonlinear control problems by adopting the distributed MPC approach. In ref 12, an improved distributed MPC strategy was proposed for the discrete-time systems. Conte et al.13 presented a new distributed MPC method for the networks of linear systems. However, the condition of the Received: March 31, 2017 Accepted: April 2, 2017 Published: April 3, 2017 5360

DOI: 10.1021/acs.iecr.7b01319 Ind. Eng. Chem. Res. 2017, 56, 5360−5394

Article

Industrial & Engineering Chemistry Research

based robust economic MPC was presented for the nonlinear systems with unknown disturbances. Farima et al.34 reviewed the development of the stochastic robust MPC algorithms in the past 10 years. On the basis of the Gaussian processes, a stochastic MPC method was developed for the propagating system disturbances in ref 35. Although great developments are gained in theory for the robust MPC algorithm, its concrete implementations are limited by the amount of online calculation in practice. How to improve the ensemble control performance of MPC algorithms is still a theme needed to be explored. At the same time, the extension of the applications of MPC schemes is the other important research field. The studies about MPC on hybrid systems have received a lot of attention. Many researchers presented their important results on the adoption of the MPC approach to piecewise affine systems, switch systems, and so on. In ref 36, an economic MPC strategy was developed to reduce the energy cost for the hybrid heating systems. A novel MPC scheme using slope information was designed for the hybrid electric vehicle platooning to improve fuel economy by Yu et al.37 Hybrid MPC methods that switch between two predictor functions were described in ref 38. In ref 39, the MPC strategy was developed to optimize the fuel economy and reduce tailpipe emissions simultaneously for hybrid electric vehicles. The hybrid MPC was presented to obtain the optimal performance and balance the handling capacity and energy consumption for an automated container terminal in ref 40. There are usually several models in the hybrid systems, and the models or dynamics are switched during the running processes of the controlled system. It is inevitable that oscillations may appear and deteriorate the relevant system performance. Therefore, an improved MPC algorithm with stronger robustness is needed. Besides the conventional applications of MPC strategies in industrial processes that are known as the slow processes, applications of an MPC scheme in the embedded systems that are small and fast are also a hot topic. In ref 41, the quadratic programming problems in the embedded constrained control applications were solved by a dual fast gradient-projection method. The computational complexity of the embedded MPC using a dual gradient method for constrained systems was analyzed by Necoara et al.42 A hardware architecture for embedded real-time MPC was discussed in ref 43. Besselmann et al.44 applied the nonlinear MPC strategy to an embedded system. In ref 45, the real-time explicit MPC was implemented on a microcontroller unit. In order to realize the applications on embedded systems, the MPC algorithms need to be simplified to meet the computation power and memory capacity of these small processors. For these uncertainties in practice, better ensemble performance of the modified MPC scheme is also the point needed to be considered. In most cases, the profits are the principal objective. On the aspect of improving profits, the economic MPC shows its superiority. In ref 46, the economic MPC scheme was presented for a nonlinear chemical reactor with an arbitrary economic objective. Heidarinejad et al.47 presented MPC schemes for a broad class of nonlinear process systems with general economic considerations. In ref 48, a number of developments in receding horizon control for the nonlinear systems with economic objectives were reviewed. A terminal constraint economic MPC strategy was developed by Amrit et al.49 In ref 50, the economic MPC method was presented for the chemical process with input rate-of-change constraints. An

information communications between the intelligent agents affects the control performance of the distributed MPC strategies greatly. It is a fact that disturbances and uncertainties exist inevitably in practice, which implies that good information exchanges between these agents can hardly be guaranteed. Besides the idea of adopting the distributed control systems, how to reduce the computation time of MPC directly was also taken into account. Many researchers contribute to the development of fast computation of MPC. In ref 14, a fast nonlinear MPC algorithm using the particle swarm optimization algorithm was implemented on a field-programmable gate array. A series of approaches to improve the speed of MPC were discussed by Wang et al.15 In ref 16, a fast numerical algorithm was presented for nonlinear receding horizon control. Barros et al.17 presented a fast MPC method for multilevel converters. In ref 18, a fast, large-scale MPC scheme based on partial enumeration was addressed. Among these improved algorithms, the explicit MPC is a significant contribution to fast computation of MPC. By employing the multiparameter quadratic programming, the repetitive computation can be done offline, which reduces the amount of online calculation largely. On the basis of wavelet multiresolution analysis, an explicit MPC strategy with low complexity was presented by Summers et al.19 Through approximation of the nonlinear constrained finite-time optimal control, an explicit nonlinear MPC was developed by Ulbig et al.20 In ref 21, an explicit MPC scheme was proposed on the basis of alternative models. An online active set strategy was introduced to construct the fast MPC by Ferreau et al.22 In ref 23, the method to acquire the explicit form of min−max MPC using a constructive algorithm was given. However, the explicit MPC scheme is just suitable for the small-scale optimization problem. For complex large-scale control systems, the demands on high speed and real-time are not always fulfilled. As mentioned above, how to simplify the computation process of the MPC algorithm for complex systems is still an open issue. In the presence of these unknown disturbances and uncertainties in practice, how to maintain the acceptable control performance of MPC strategies is also a vital topic. In this branch, robust MPC was listed as the research object. Combining with the Lyapunov function, linear matrix inequality (LMI), state feedback control law, min−max problem, etc., various results of robust MPC algorithms are developed. In ref 24, an improved robust MPC scheme was addressed for the constrained linear systems. In order to cope with the inputsaturated polytopic linear parameter varying discrete-time systems, a new robust MPC approach was presented by Casavola et al.25 A robust hybrid MPC structure was proposed for nonlinear systems with input constraints and uncertain variables by Mhaskar et al.26 In ref 27, a novel robust MPC strategy was developed for input-saturated systems with polytopic model uncertainties. A robust MPC scheme was put forward for nonlinear systems subject to constraints, uncertainty, and actuator faults by Mhaskar.28 An overview of robust MPC was presented in ref 29. Raffo et al.30 proposed a nonlinear robust MPC algorithm to handle the path following problem of a quadrotor helicopter. At the same time, tubebased robust MPC schemes and stochastic robust MPC strategies also draw a lot of attention. In ref 31, a collective neurodynamic approach to robust MPC for nonlinear systems affected by bounded uncertainties was proposed. The tubedbased MPC scheme for a class of constrained linear switch systems was considered by Hariprasad et al.32 In ref 33, a tube5361


Article

Industrial & Engineering Chemistry Research economic MPC scheme without terminal constraints was discussed in ref 51. By considering the economic performance index, more benefits are obtained, but the complexity of the whole system is increased simultaneously. In order to make better use of the economic performance index, the improved MPC algorithm with lower complexity and stronger control performance is required. Besides the main topics mentioned above, there are also other research topics about the state space MPC strategies. How to develop an improved MPC algorithm with better ensemble control performance, acceptable amount of calculation, and lower complexity is still an open topic, which will promote the applications of MPC approaches to other industries simultaneously. Note that, in most cases, the process information is not fully used in the state variables, and the degree of freedom for controller design is limited. In the past few years, we developed approaches to improve the control performance of MPC using the extended state space models and the extended nonminimal state space (ENMSS) models in which the state variables and tracking error are combined and regulated separately.52−56 Compared with traditional state space MPC, the choice of the corresponding weighting matrices can be more flexible in the controller design based on such models. At the same time, the proposed strategy can also be used to optimize the PID controllers.57−61 We introduced the control law of PID control into the framework of such models based MPC and obtained PID controllers that inherit the simple structure of conventional PID controllers and the same control performance as the proposed MPC algorithms. In the following sections, a summary of our recent work on such models-based MPC and the corresponding industrial applications62,63 are presented. The corresponding results are compared with typical PID and MPC methods in the literature.64,65 Finally, further insights into the proposed MPC are verified by comparison with the typical work by Rawlings.66 In order to maintain consistency of the whole article, only single-input single-output (SISO) modelbased MPC strategies are introduced in detail, and note that extension to the multi-input multi-output (MIMO) cases is straightforward.

Δxm(k + 1) = A mΔxm(k) + BmΔu(k − d) Δym (k + 1) = CmΔxm(k + 1)

By taking ym(k) as the state variables, we can obtain a new state space model x(k + 1) = Ax(k) + BΔu(k − d) ym (k + 1) = Cx(k + 1)

(3)

where ⎡ Δx (k)⎤ ⎡ A ⎡ Bm ⎤ 0⎤ m ⎥; A = ⎢ m ⎥; B = ⎢ ⎥; x(k) = ⎢ ⎢⎣ y (k) ⎥⎦ ⎢⎣CmA m 1 ⎥⎦ ⎢⎣CmBm ⎥⎦ m C = [0 1]

Here, 0 in these system matrices are zero matrices with appropriate dimensions. The reference trajectory can be set as follows when the smoothing factor α is introduced s(k) = y(k) s(k + i) = α iy(k) + (1 − α i)ys

i = 1, 2, ···, P , ···, P + d (4)

where y(k) is the actual process output, ys is the set-point, and P is the prediction horizon. Here, choose M as the control horizon. The future predictions from sampling time instant k can be acquired based on eq 3 X = Fx(k + d) + ΦΔU Y = Tx(k + d) + GΔU

(5)

where ⎤ ⎡ x(k + d + 1) ⎤ ⎡ Δu(k) ⎥ ⎢ ⎥ ⎢ ⎢ x(k + d + 2) ⎥ ⎢ Δu(k + 1) ⎥ X=⎢ ⎥ ⎥ ; ΔU = ⎢ ⋮ ⋮ ⎥ ⎢ ⎥ ⎢ ⎢ x(k + d + P)⎥ ⎢ Δu(k + M − 1)⎥ ⎦ ⎣ ⎦ ⎣

2. STATE SPACE MODEL PREDICTIVE CONTROL In this section, the conventional state space MPC is first described. Second, we propose an improved version of MPC and PFC based on the extended state space model. Finally, the proposed control strategy is tested via the control of an inverseresponse process. 2.1. Conventional State Space Model Predictive Control. Here, the traditional idea of using state space MPC for set-point tracking will be briefly discussed. For simplicity, we consider the SISO controlled process with the following state space model

⎡ y (k + d + 1) ⎤ ⎡ s (k + 1 + d ) ⎤ ⎢ m ⎥ ⎢ ⎥ ⎢ y (k + d + 2) ⎥ ⎢ ⎥ s k d + + ( 2 ) ⎥; S = ⎢ Y=⎢ m ⎥ ⎢ ⎥ ⋮ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ s(k + P + d)⎥ ⎣ ⎦ y k d P + + ( ) ⎢⎣ m ⎥⎦ ⎡ B 0 0 ⎡A⎤ ⎢ AB B 0 ⎢ 2⎥ ⎢ ⎢A ⎥ F = ⎢ ⎥ ; Φ = ⎢ A2 B AB B ⎢ ⋮ ⎢ ⎥ ⋮ ⋮ ⋮ ⎢ ⎣⎢ AP ⎦⎥ ⎢ P−1 ⎣ A B AP − 2 B AP − 3B

xm(k + 1) = A mxm(k) + Bmu(k − d) ym (k + 1) = Cmxm(k + 1)

(2)

(1)

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ AP − M B ⎦ ⋯ ⋯

0 0

⎡ CA ⎤ ⎢ 2⎥ ⎢ CA ⎥ T=⎢ ; ⋮ ⎥ ⎢ ⎥ ⎢⎣CAP ⎥⎦

where xm(k), u(k), and ym(k) are the state, input, and output of the model at time instant k, respectively. Here, d is the dead time, and Am, Bm, Cm are the system matrices with appropriate dimensions. By adding the difference operator Δ to eq 1, we can obtain the following model 5362


Article

Industrial & Engineering Chemistry Research ⎡ CB 0 0 ⎢ CAB CB 0 ⎢ 2 G = ⎢ CA B CAB CB ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ P−1 P − 2 P ⎣CA B CA B CA − 3B

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ CAP − M B ⎦ ⋯ ⋯

⎡ Am ⎢ ⎢ 0 A1 = ⎢ 0 ⎢ ⎢⋮ ⎢⎣ 0

0 0

Δx1(k + 1) = A1Δx1(k) + B1Δu(k)

Yc = Y + Epre

Δym (k + 1) = C1Δx1(k + 1)

where

(12)

The reference trajectory can be set as follows:

⎡ y (k + d + 1) ⎤ ⎥ ⎢ mc ⎢ y (k + d + 2) ⎥ ⎥; Epre = [ epre(k) epre(k) ⋯ epre(k)]Τ Yc = ⎢ mc ⎥ ⎢ ⋮ ⎥ ⎢ ⎢⎣ ymc (k + d + P)⎥⎦

r(k) = y(k) r(k + i) = α iy(k) + (1 − α i)ys

i = 1, 2, ···, P

(13)

where y(k) is the actual process output, ys is the set-point, α is the smoothing factor, and P is the prediction horizon. Then, the output tracking error at time instant k can be described as

epre(k) = y(k) − ym (k)

e(k) = ym (k) − r(k)

The cost function is chosen as

(14)

The formulation of e(k + 1) can be estimated by combining eqs 12−14

(7)

where Q and L are the corresponding weight matrices with appropriate dimensions. By taking a derivative of J, the optimal control law can be obtained as follows: ΔU = −(GΤQG + L)−1GΤQ (Tx(k + d) + Epre − S)

⎡Cm ⎤Τ ⎢ ⎥ ⎢0⎥ ⎢⋮⎥ ⎢ ⎥ ⎢0⎥ ⎢⎣ 0 ⎥⎦

By adding the difference operator Δ to eq 11, the following holds

Consider the prediction error between actual process output and model output, the corrected model prediction output can be expressed as

J = (S − Yc)Τ Q (S − Yc) + ΔUΤLΔU

0 ⋯ 0 Bm ⎤ ⎡0⎤ ⎥ ⎢ ⎥ 0 ⋯ 0 0⎥ ⎢1 ⎥ ⎥ 1 0 ⋯ 0 ; B1 = ⎢ 0 ⎥ ; C1 = ⎥ ⎢ ⎥ ⋮ ⋱ ⋮ ⋮⎥ ⎢⋮⎥ ⎥ ⎣0⎦ 0 ⋯ 1 0⎦

e(k + 1) = e(k) + C1A1Δx1(k) + C1B1Δu(k) − Δr(k + 1) (15) T

Constructing a new state variable as z(k) = [Δx1(k) e(k)]T, we can obtain the extended state space (ESS) model z(k + 1) = Az(k) + BΔu(k) + C Δr(k + 1)

(8)

(16)

where

Let

⎡ A1 0 ⎤ ⎡ B1 ⎤ ⎡0⎤ ⎥; B = ⎢ ⎥; C = ⎢ ⎥ A=⎢ ⎣−1⎦ ⎣C1A1 1 ⎦ ⎣C1B1⎦

KF = (GΤQG + L)−1GΤQT KS = (GΤQG + L)−1GΤQ

(9)

Here, the 0 in eq 17 are zero matrices with appropriate dimensions. Choose M as the control horizon and define

Then, the control input at time instant k is Δu(k) = −kFx(k + d) + kS(S − Epre)

(10)

⎡ z(k + 1) ⎤ ⎡ ⎤ ⎡ Δr(k + 1) ⎤ Δu(k) ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ z(k + 2) ⎥ ⎢ Δu(k + 1) ⎥ ⎢ Δr(k + 2) ⎥ Z=⎢ ⎥ ; ΔU = ⎢ ⎥ ; ΔR = ⎢ ⎥ ⋮ ⋮ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ z(k + P)⎦ ⎣ Δu(k + M − 1)⎦ ⎣ Δr(k + P)⎦

where kF and kS are the first rows of KF and KS, respectively. Remark 1. In eq 8, Δu(k − d),···,Δu(k − 1) are the past inputs, and the computations of these items are meaningless. In order to simplify the computation process, we can also take Δu(k − d),···,Δu(k − 1) as the elements of the new state variables; then, eq 3 can be transformed into another state space model in which Δu(k) remains as the input, which is shown in the next section. 2.2. Extended State Space Model Predictive Control. 2.2.1. Extended State Space Model Based MPC [52]. Consider the same state space model in eq 1 and choose the new state variable as x1(k) = [xm(k), u(k − 1), ···, u(k − d)]T; then, eq 1 can be rewritten as

(18)

The future state variables form sampling time instant k based on eq 16 can be described as Z = Sz(k) + F ΔU + θ ΔR

(19)

where ⎡ B 0 0 ⎡A⎤ ⎢ B 0 ⎢ 2⎥ ⎢ AB ⎢A ⎥ 2 ⎢ S = ⎢ ⎥; F = A B AB B ⎢ ⋮ ⎢ ⎥ ⋮ ⋮ ⎢ ⋮ ⎢⎣ AP ⎥⎦ ⎢ P−1 P 2 P − ⎣ A B A B A − 3B

x1(k + 1) = A1x1(k) + B1u(k) ym (k + 1) = C1x1(k + 1)

(17)

(11)

where 5363

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ AP − M B ⎦ ⋯ ⋯

0 0


Article

Industrial & Engineering Chemistry Research ⎡ C 0 0 ⎢ AC C 0 ⎢ 2 θ=⎢ AC AC C ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ P−1 P − 2 P ⎣ A C A C A − 3C

0 0⎤ ⎥ 0 0⎥ 0 0⎥ ⎥ ⋱ ⋮⎥ ⎥ ⋯ C⎦

⎡ z(k + 1) ⎤ ⎡ Δr(k + 1) ⎤ ⎥ ⎢ ⎥ ⎢ ⎢ z(k + 2) ⎥ ⎢ Δr(k + 2) ⎥ Z=⎢ ⎥ ⎥ ; ΔR = ⎢ ⋮ ⋮ ⎥ ⎢ ⎥ ⎢ ⎢ z(k + P)⎥ ⎢ Δr(k + P)⎥ ⎦ ⎣ ⎦ ⎣

and combine eq 24 with eq 27; the future state variables of PFC from sampling time instant k can be described as

Considering the prediction error, the cost function for ESS model-based MPC is chosen as J(k) = (Z + Epre)Τ Q (Z + Epre) + ΔUΤLΔU

(28)

Z = Sz(k) + φη − Fu(k − 1) + θ ΔR

(20)

(29)

where

where

⎡ ⎤ BT0 ⎢ ⎥ ⎢ ⎥ (AB − B)T0 + BT1 ⎡A⎤ ⎡ B ⎤ ⎢ ⎥ ⎢ 2⎥ ⎥ ⎢ AB ⎥ ⎢(A2B − AB)T0 + (AB − B)T1 + BT2 ⎥ ⎢A ⎥ ⎢ S = ⎢ ⎥; F = ⎢ ;ϕ=⎢ ⎥ ⋮ ⎥ ⋮ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ P−1 ⎥ ⎢ P−1 ⎥ ⎢⎣ AP ⎥⎦ ⎣ A B⎦ ⎢ ⎥ k k−1 ( A B A B ) T BT − + ∑ 1 1 P k P − − − ⎢ ⎥ ⎣ k=1 ⎦

Epre = [E , E , ···, E]Τ ; E = [ 0 epre(k)]

Here, epre(k) = y(k) − ym(k) is the prediction error between the actual process output and the model output, and 0 in E is the zero vector with appropriate dimension. Note that Q and L are the weighting matrices for the state variables and control inputs, respectively, and it is further verified that the proposed model shown in eq 16 will facilitate the weighting on the state variables and the tracking error. Due to the above merits, more degrees of freedom are provided, and improved control performance may be anticipated. By combining eq 19 with eq 20 and taking a derivative of J(k), the optimal control law can be obtained as follows: ΔU = −(FΤQF + L)−1FΤQ (Sz(k) + θ ΔR + Epre)

⎡ C 0 0 ⎢ C 0 ⎢ AC 2 ⎢ θ= AC AC C ⎢ ⋮ ⋮ ⋮ ⎢ ⎢ P−1 ⎣ A C AP − 2 C AP − 3C

(21)

0 0⎤ ⎥ 0 0⎥ 0 0⎥ ⎥ ⋱ ⋮⎥ ⎥ ⋯ C⎦

Let Taking into account the prediction error, the cost function for the ESS model-based PFC is chosen as

KF = (FΤQF + L)−1FΤQS KS = (FΤQF + L)−1FΤQ

J(k) = (Z + Epre)Τ Q (Z + Epre)

(22)

Then, the control input at time instant k is Δu(k) = −kFz(k) − kS(θ ΔR + Epre)

where (23)

Epre = [E , E , ···, E]Τ ; E = [ 0 epre(k)]

where kF and kS are the first rows of KF and KS, respectively. 2.2.2. Extension to Predictive Functional Control (PFC).53 Adopting the same model formulation procedure shown in Section 2.2.1, we can obtain the following equation

Here, epre(k) = y(k) − ym(k) is the prediction error between the actual process output and the model output. Also, Q is the weighting matrix for the state variables, and 0 in E is the zero vector with appropriate dimension. By synthesizing eq 29 with eq 30 and taking a derivative of J(k), the optimal control law can be described as follows:

z(k + 1) = Az(k) + Bu(k) − Bu(k − 1) + C Δr(k + 1) (24)

It is known that PFC action is associated with the process character and base functions of the set-point as follows:

η = −(φ ΤQφ)−1φ ΤQ (Sz(k) + θ ΔR − Fu(k − 1) + Epre)

N

u(k + i) =

∑ μj f j (i) j=1

(31)

By further defining

(25)

where u(k + i) is the control input at time instant k + i, μj is the weight coefficient, f j(i) is the value of the base function at time instant k + i, and N is the degree of the base functions. Define

μ1 = −( 1, 0, ⋯ , 0)(φΤQφ)−1φΤQ (Sz(k) − Fu(k − 1) + θ ΔR + Epre) = − h1z(k) + hu1u(k − 1) − m1(θ ΔR + Epre) μ2 = −( 0, 1, ⋯ , 0)(φΤQφ)−1φΤQ (Sz(k) − Fu(k − 1) + θ ΔR + Epre)

Ti = [f1 (i), f2 (i), ⋯ , fN (i)], (i = 0, 1, ⋯ , P − 1) η = [μ1 , μ2 , ⋯ , μN ]Τ

= − h2z(k) + hu2u(k − 1) − m2(θ ΔR + Epre) ⋮

(26)

μN = −( 0, 0, ⋯ , 1)(φΤQφ)−1φΤQ (Sz(k) − Fu(k − 1) + θ ΔR + Epre)

then eq 25 can be rewritten as

u(k + i) = Tiη

(30)

= − hN z(k) + huN u(k − 1) − mN (θ ΔR + Epre)

(32)

(27)

the control input at time instant k can be described as

Denote 5364


Article


Figure 1. Closed-loop responses under case 1.

Figure 2. Closed-loop responses under case 2. N

u(k) =

at time instant k = 0, and the load disturbance with amplitude −0.2 is added to the process output at time instant k = 150. For the two MPC methods, the control parameters are the same with P = 21, M = 1, L = 0.001, and α = 0.65. For the proposed ESS-based MPC, the weighting matrix is Q = diag(Q1,Q2,···,Q21), and Qj = diag(0.02,0,0,0,1)(j = 1,2,···,21) is chosen. For the conventional state space model MPC, since its model is not expanded as an ESS model, it can only adopt Qj = 1 for the output error weighting. In practice, disturbance and uncertainty exist inevitably and cause model/plant mismatches. It is more meaningful for us to verify the ensemble performance of controllers under such situations. Here, we suppose the process and its model are mismatched, where the matrices A and B of the actual process are changed to the following

∑ μj f j (0) j=1

= −Hz(k) + Huu(k − 1) − M(θ ΔR + Epre)

(33)

where N

H=

∑ f j (0)hk , (k = 1, 2, ⋯, N ) j=1 N

Hu =

∑ f j (0)huk ,(k = 1, 2, ⋯, N ) j=1 N

M=

∑ f j (0)mk , (k = 1, 2, ⋯, N ) j=1

Case 1:

2.3. Simulation Examples. Consider an inverse-response process with unstable open loop dynamics and dead time

⎡ 0.95 ⎤ ⎡ 1.1053 0 ⎤ ;B=⎢ A=⎢ ⎥ ⎣−0.012 0.8745⎥⎦ ⎣ 0.0836 ⎦

x(k + 1) = Ax(k) + Bu(k − 2) y(k + 1) = Cx(k + 1)

(34)

(35)

Case 2:

where

⎡ 1.15 ⎤ ⎡ 1.1053 0 ⎤ ;B=⎢ A=⎢ ⎣ 0.0873⎥⎦ ⎣−0.011 0.8804 ⎥⎦

⎡1.1053 ⎡ 1 ⎤ 0 ⎤ A=⎢ ⎥ ; B = ⎢⎣ ⎥ ; C = [0 1] ⎣ −0.01 0.8186 ⎦ 0.0858 ⎦

(36)

In Figure 1, we can easily find that the two methods achieve the set-point tracking and disturbance rejection successfully. From an overall perspective, the control performance of ESS model-based MPC is better than that of conventional state space MPC. Overshoot and oscillations in ESS model-based MPC method are smaller than those of the traditional state space MPC method. In summary, the tracking error of the

Here, conventional state space MPC is introduced to be compared with ESS model-based MPC, and the simulation procedure is as follows. It is assumed that the process variables are at zeros when the control system is started, and this applied to the simulations in the other sections. A unit step change is added to the set-point 5365


Article

Industrial & Engineering Chemistry Research conventional state space MPC method is bigger than that of ESS model-based MPC. The situation in Figure 2 is the same as that in Figure 1, and the ensemble control performance using ESS model-based MPC is improved compared with conventional state space MPC.

x(k + 1) = Ax(k) + BΔu(k) y(k + 1) = Cx(k + 1)

where ⎡ Δx (k)⎤ ⎡ A ⎡ Bm ⎤ 0⎤ m ⎥; A = ⎢ m ⎥; B = ⎢ ⎥; C = [0 1] x(k) = ⎢ ⎢⎣CmA m 1 ⎥⎦ ⎢⎣CmBm ⎥⎦ ⎢⎣ y(k) ⎥⎦

3. INPUT−OUTPUT MODEL-BASED ENMSS MODEL PREDICTIVE CONTROL An input−output model is common in practice, and we need to transform it into a form of state space model, and then the state space MPC can be employed. In this section, the conventional nonminimal state space MPC (NMSSMPC) and extended nonminimal state space MPC (ENMSSMPC) based on input− output models are presented, and the control performance of these methods is compared through the case studies on a double integrating plant and a glasshouse microclimate process. At the same time, the relevant constrained MPC scheme is also developed. 3.1. Conventional (Nonminimal State Space) NMSS Model Predictive Control.64 Suppose the model of the controlled process can be described as follows:

and 0 is a zero matrix with appropriate dimension. The reference trajectory can be set as follows when the smoothing factor α is introduced s(k) = y(k) i

s(k + i) = α y(k) + (1 − α i)ys

(37)

⎡ s(k + 1) ⎤ ⎢ ⎥ ⎢ s(k + 2) ⎥ S=⎢ ⎥ ⎢ ⋮ ⎥ ⎢ ⎥ ⎣ s(k + P)⎦

Δy(k + 1) + H1Δy(k) + H2Δy(k − 1) + ··· + HmΔy(k − m + 1)

(38)

X = Fx(k) + ΦΔU

Τ

Δxm(k) = [Δy(k), Δy(k − 1), ···, Δy(k − m + 1),

Y = Tx(k) + GΔU

Δu(k − 1), Δu(k − 2), ···, Δu(k − n + 1)]

(44)

where (39)

Then, eq 38 can be rewritten as Δxm(k + 1) = A mΔxm(k) + BmΔu(k) (40)

where ⎡− H1 − H2 ⎢ 0 ⎢ 1 ⎢ 0 1 ⎢ ⋮ ⎢ ⋮ Am = ⎢ 0 0 ⎢ 0 ⎢ 0 ⎢ 0 0 ⎢ ⋮ ⎢ ⋮ ⎢⎣ 0 0

(43)

then the following holds

The NMSS vector Δxm(k) can be chosen as

Δy(k + 1) = CmΔxm(k + 1)

(42)

⎡ y(k + 1) ⎤ ⎤ ⎡ ⎡ x(k + 1) ⎤ Δu(k) ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ y(k + 2) ⎥ ⎢ Δu(k + 1) ⎥ ⎢ x(k + 2) ⎥ X=⎢ Δ = = ; U ; Y ⎥; ⎢ ⎥ ⎢ ⎥ ⋮ ⋮ ⋮ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ + Δ + − x ( k P ) u ( k M 1) + y ( k P ) ⎦ ⎣ ⎦ ⎣ ⎦ ⎣

where y(k) and u(k) are the output and input of the process at time instant k, respectively. By adding the difference operator Δ to eq 37, the following holds = L1Δu(k) + L 2Δu(k − 1) + ··· + LnΔu(k − n + 1)

i = 1, 2, ···, P

where ys is the set-point, and P is the prediction horizon. Here, denote M as the control horizon. The future state variables from sampling time instant k are based on eq 41 and define

y(k + 1) + H1y(k) + H2y(k − 1) + ··· + Hmy(k − m + 1) = L1u(k) + L 2u(k − 1) + ··· + Lnu(k − n + 1)

(41)

⋯ − Hm − 1 − Hm L 2 ⋯ Ln − 1 Ln ⎤ ⎥ ⋯ 0 0 0 ⋯ 0 0⎥ ⋯ 0 0 0 ⋯ 0 0⎥ ⎥ ⋯ ⋮ ⋮ ⋮ ⋯ ⋮ ⋮⎥ ⋯ 1 0 0 ⋯ 0 0⎥ ⎥ ⋯ 0 0 0 ⋯ 0 0⎥ ⋯ 0 0 1 ⋯ 0 0⎥ ⎥ ⋯ ⋮ ⋮ ⋮ ⋯ ⋮ ⋮⎥ ⋯ 0 0 0 ⋯ 1 0 ⎥⎦

⎡ B 0 0 ⎡A⎤ ⎢ B 0 ⎢ 2⎥ ⎢ AB ⎢A ⎥ 2 ⎢ F = ⎢ ⎥; Φ = A B AB B ⎢ ⋮ ⎢ ⎥ ⋮ ⋮ ⋮ ⎢ ⎢⎣ AP ⎥⎦ ⎢ P−1 ⎣ A B AP − 2 B AP − 3B

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ AP − M B ⎦

⎡ CB 0 0 ⎡ CA ⎤ ⎢ 0 CAB CB ⎢ 2⎥ ⎢ ⎢ CA ⎥ 2 ⎢ T=⎢ ; G = CA B CAB CB ⎢ ⋮ ⎥ ⎥ ⎢ ⋮ ⋮ ⋮ ⎢ ⎣⎢CAP ⎦⎥ ⎢ P−1 − P 2 P ⎣CA B CA B CA − 3B

⋯ ⋯

⋯ ⋯

0 0

⎤ ⎥ ⎥ 0 ⎥ ⋯ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ CAP − M B ⎦ 0 0

The cost function is chosen as J = (S − Y )Τ Q (S − Y ) + ΔUΤLΔU

(45)

where Q and L are the corresponding weighting matrices with appropriate dimensions. By taking a derivative of J, the optimal control law can be presented as follows:

Bm = [ L1 0 0 ⋯ 1 0 ⋯ 0 ]Τ

Cm = [1 0 0 ⋯ 0 0 0 0]

ΔU = −(GΤQG + L)−1GΤQ (Tx(k) − S)

Combine the NMSS vector with y(k), a new state space model can be presented as

(46)

Let 5366


Article

Industrial & Engineering Chemistry Research ⎡ Cm 0 0 ⎢ 0 Cm ⎢ A mCm ⎢ θ = ⎢ A m 2 Cm A mCm Cm ⎢ ⋮ ⋮ ⋮ ⎢ ⎢ A P − 1C A P − 2 C A P − 3C ⎣ m m m m m m

KF = (GΤQG + L)−1GΤQT KS = (GΤQG + L)−1GΤQ

(47)

Then, the control vector at time instant k is Δu(k) = −kFx(k) + ksS

(48)

where kF and ks are the first rows of KF and KS, respectively. 3.2. Extended NMSS (ENMSS) Model-Based Control.54 3.2.1. ENMSS Model-Based MPC. Here, we consider the same difference equation model in eq 37, and we can obtain the state space model in eq 40 first. The reference trajectory can be chosen as follows:

i = 1, 2, ···, P

(57)

ΔU = −(FΤQF + L)−1FΤQ (Sz(k) + θ ΔR ))

(58)

Let KF = (FΤQF + L)−1FΤQS KS = (FΤQF + L)−1FΤQθ

Combining eq 40 with eq 50, we can obtain the formulation of e(k + 1)

(59)

Then the control vector at time instant k is

e(k + 1) = e(k) + CAΔx(k) + CBΔu(k) − Δr(k + 1)

Δu(k) = −kFz(k) − ksΔR

(51)

(60)

where kF and ks are the first rows of KF and KS, respectively. 3.2.2. Extension to PFC.55 Here, we directly adopt eq 53 in Section 3.2.1 as the process model, and it can be further written as

The new state vector is chosen as

(52)

z(k + 1) = A mz(k) + Bmu(k) − Bmu(k − 1) + CmΔr(k + 1) (61)

Then the ENMSS model can be described as

It is known that PFC action is associated with the process character and the base functions of set-point as follows:

(53)

where

N

⎡ A 0⎤ ⎡B⎤ ⎡0⎤ ;B = ;C = Am = ⎢ ⎣CA 1 ⎥⎦ m ⎢⎣CB ⎥⎦ m ⎢⎣−1⎥⎦

u(k + i) =

Ti = [f1 (i), f2 (i), ⋯ , fN (i)], (i = 0, 1, ⋯ , P − 1) η = [μ1 , μ2 , ⋯ , μN ]Τ

u(k + i) = Tiη

The future state variables from sampling time instant k based on eq 53 can be described as


where

⋯ ⋱ ⋯

(64)

Denote

(56)

⋯

(63)

Then, eq 62 can be rewritten as

(55)

⋯

(62)

where u(k + i) is the control input at time instant k + i, μj is the weight coefficient, f j(i) is the value of the base function at time instant k + i, and N is the degree of the base functions. Define

⎡ z(k + 1) ⎤ ⎡ ⎤ ⎡ Δr(k + 1) ⎤ Δu(k) ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ z(k + 2) ⎥ ⎢ Δu(k + 1) ⎥ ⎢ Δr(k + 2) ⎥ Δ = Δ = Z=⎢ U R ; ; ⎥ ⎢ ⎥ ⎢ ⎥ ⋮ ⋮ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ z(k + P)⎦ ⎣ Δu(k + M − 1)⎦ ⎣ Δr(k + P)⎦

Z = Sz(k) + F ΔU + θ ΔR

∑ μj f j (i) j=1

(54)

and the 0 in eq 54 are a zero matrices with appropriate dimensions. Define

⎡ Bm 0 0 ⎡ Am ⎤ ⎢ Bm 0 ⎢ ⎥ ⎢ A mBm ⎢ A m2 ⎥ ⎢ 2 ⎥ ; F = ⎢ A m Bm S=⎢ A mBm Bm ⎢ ⋮ ⎥ ⎢ ⋮ ⋮ ⋮ ⎢ P⎥ ⎢ ⎣ Am ⎦ ⎢ A P − 1B A P − 2 B A P − 3B ⎣ m m m m m m

⋯

where Q and L are the weighting matrices for the state variables and control inputs, respectively. By combining eq 56 with eq 57 and taking a derivative of J(k), the optimal control law can be obtained as follows:

(50)

z(k + 1) = A mz(k) + BmΔu(k) + CmΔr(k + 1)

⋱

J(k) = ZΤQZ + ΔUΤLΔU

where ys is the set-point, α is the smoothing factor, and P is the prediction horizon. The tracking error can be further presented as

⎡ Δx(k)⎤ ⎥ z(k ) = ⎢ ⎢⎣ e(k) ⎥⎦

0

The cost function for ENMSS model-based MPC is chosen

(49)

e(k) = y(k) − r(k)

0

0⎤ ⎥ 0⎥ ⎥ 0⎥ ⎥ ⋮⎥ Cm ⎥⎦

as

r(k) = y(k) r(k + i) = α iy(k) + (1 − α i)ys

0

⎤ ⎥ 0 ⎥ ⎥ 0 ⎥ ⎥ ⋮ ⎥ P−M Am Bm ⎥⎦ 0

(65)

and combine eq 61 with eq 64; we can obtain the future state variables of PFC from sampling time instant k as Z = Sz(k) + φη − Fu(k − 1) + θ ΔR 5367

(66)


Article

Industrial & Engineering Chemistry Research N

where H=

⎡ Am ⎤ ⎡ Bm ⎤ ⎢ ⎥ ⎢ ⎥ 2 ⎢ Am ⎥ ⎢ A mBm ⎥ S=⎢ ⎥; F = ⎢ ⎥; ⋮ ⎢ ⋮ ⎥ ⎢ ⎥ ⎢ P⎥ ⎢ A P − 1B ⎥ ⎣ m m⎦ ⎣ Am ⎦

j=1 N

Hu =

∑ f j (0)huk ,(k = 1, 2, ⋯, N ) j=1 N

M=

0 0 0 ⋱ ⋯

∑ f j (0)mk , (k = 1, 2, ⋯, N ) j=1

⎡ ⎤ BmT0 ⎢ ⎥ ⎢ ⎥ (A mBm − Bm)T0 + BmT1 ⎢ ⎥ ⎢ (A m 2 Bm − A mBm)T0 + (A mBm − Bm)T1 ⎥ ⎢ ⎥ + BmT2 ϕ=⎢ ⎥ ⎢ ⎥ ⋮ ⎢ ⎥ ⎢ P−1 ⎥ k k−1 ⎢ ∑ (A m Bm − A m Bm)TP − 1 − k + BmTP − 1⎥ ⎢⎣ k = 1 ⎥⎦

⎡ Cm 0 0 ⎢ 0 Cm ⎢ A mCm ⎢ θ = ⎢ A m 2 Cm A mCm Cm ⎢ ⋮ ⋮ ⋮ ⎢ ⎢ A P − 1C A P − 2 C A P − 3C ⎣ m m m m m m

∑ f j (0)hk , (k = 1, 2, ⋯, N )

3.3. Constrained Model Predictive Control.56 Constraint dealing is also an important issue in MPC. In this section, the constraint dealing is discussed using ENMSS model-based MPC in Section 3.2.1 as an example. For simplicity, the control horizon and the prediction horizon are both chosen as P. For unconstrained MPC, the optimal control law can be presented as ΔU = −(FΤQF + L)−1FΤQ (Sz(k) + θ ΔR ))

(71)

The following constraints on the inputs and outputs are introduced

0⎤ ⎥ 0⎥ ⎥ 0⎥ ⎥ ⋮⎥ Cm ⎥⎦

ymin ≤ y(k) ≤ ymax , umin ≤ u(k) ≤ umax

(72)

where ymin, ymax, umin, and umax are the constraints. It is obvious that eq 71 will not always be optimal under constrained cases. Here, we introduce a new vector Z = [Z1T, Z2T,···,ZPT]T as follows:

The cost function for the ENMSS model-based PFC is chosen as

Z = (FT QF + L)ΔU + FT Q (Sz(k) + θ ΔR )

(73)

The new cost function is chosen as follows: J(k) = ZΤQZ

(67)

P

J[ΔU ] = min ∑ |Zi|

where Q is the weighting matrix for the state variables. By synthesizing eq 66 and eq 67, and taking a derivative of J(k), the optimal control law can be described as follows: η = −(φ ΤQφ)−1φ ΤQ (Sz(k) + θ ΔR − Fu(k − 1))

The goal of eq 74 is to drive the process output to the setpoint and ensure the steady process input simultaneously. Define

(68)

Zi = Zi + − Zi−

By further defining

= −h1z(k) + hu1u(k − 1) − m1ΔR Τ

(75)

where Zi+ > 0,Zi− > 0. Then, eq 74 can be rewritten as

μ1 = −(1, 0, ⋯ , 0)(φΤQφ)−1φΤQ (Sz(k) − Fu(k − 1) + θ ΔR )

P

J[ΔU ] = min ∑ |Zi + − Zi−|

−1 Τ

μ2 = −(0, 1, ⋯ , 0)(φ Qφ) φ Q (Sz(k) − Fu(k − 1) + θ ΔR )

(76)

i=1

= −h2z(k) + hu2u(k − 1) − m2ΔR

Finally, the constrained MPC can be converted to the following linear programming problem

⋮ μN = −(0, 0, ⋯ , 1)(φΤQφ)−1φΤQ (Sz(k) − Fu(k − 1) + θ ΔR )

P

= −hN z(k) + huN u(k − 1) − mN ΔR

min ∑ |Zi + − Zi−|

(69)

(77)

i=1

the control input at time instant k can be described as

subject to

N

u(k) =

(74)

i=1

ymin ≤ y(k + i) ≤ ymax (i = 1, 2, ⋯ , P)

∑ μj f j (0) = −Hz(k) + Huu(k − 1) − M ΔR

umin ≤ u(k + i − 1) ≤ umax (i = 1, 2, ⋯ , P)

j=1

(70)

(78)

Express all constraints in terms of the manipulated vector ΔU, the following holds

where 5368


Article



⎡ ⎤ u(k) ⎢ ⎥ ⎡I ⎤ ⎢ u(k + 1) ⎥ ⎢ I ⎥ ⎢ ⎥ = ⎢ ⎥u(k − 1) ⋮ ⎢ ⎥ ⎢I ⎥ ⎢ u(k + P − 1)⎥ ⎢⎣ I ⎥⎦ ⎣ ⎦ ⎡I ⎢ I +⎢ ⎢⋮ ⎢ ⎣I

0 ⋯ I ⋯ ⋮ ⋱ I

I

⎡ ⎤ Δu(k) 0 ⎤⎢ ⎥ ⎥ 0 ⎥⎢ Δu(k + 1) ⎥ ⎢ ⎥ 0 ⎥⎢ ⋮ ⎥ ⎥ I ⎦⎢⎣ Δu(k + P − 1)⎥⎦

3.4. Simulation Examples. 3.4.1. SISO Case Study. Here, the double integrating plant is introduced,64 and its discretetime transfer function is G (z ) =

(79)

C1umin − C2u(k − 1) ≤ C ΔU ≤ C1umax − C2u(k − 1) (80)

where C, C1, and C2 are the corresponding matrices. In terms of output constraints, the following equations can be obtained based on eq 40

(81)

where C̅ = [0 Iq]. Then, the following holds ⎡ y(k + 1) ⎤ ⎢ ⎥ ⎡ Δx (k)⎤ ⎢ y(k + 2) ⎥ m ̂ ΔU ⎥ + CF Y=⎢ ⎥ = Ĉ ⎢ ⋮ ⎣⎢ y(k) ⎥⎦ ⎢ ⎥ ⎢ ⎥ ⎣ y(k + P)⎦

(84)

In this section, conventional NMSS model-based MPC is adopted to be compared with ENMSS model-based MPC, and the simulation procedure is as follows. A unit step change is added to the set-point at time instant k = 0, and a step change of load disturbance with amplitude of −0.1 is added to the process at time instant k = 150. For the two MPC methods, the control parameters are the same with P = 20, M = 3, and α = 0.65. For the ENMSS model-based MPC, it has more degrees of freedom to improve the control performance by choosing its weighting matrices as Q = diag(Q1,Q2,···,Q20) and L = diag(L1, L2, L3), and here, Qj = diag(1,1,0,1) (j = 1,2,···,20) and Li = 1 (i = 1,2,3) are adopted. For the conventional NMSS model-based MPC, since its model is only formulated for the cost function to include the output tracking error and the control input, only Qj = 1 (j = 1,2,···,20) and Li = 1 (i = 1,2,3) can be adopted. Note that disturbance and uncertainty exist inevitably in practice and cause the model/plant mismatch, which shows that it is more meaningful to verify the ensemble performance of controllers under such situations. Supposing the process and its model are mismatched, the discrete-time transfer function of actual process is changed to

Combining eq 72 and eq 79, the constraints of manipulated variables ΔU can be presented as

⎡ Δx (k)⎤ m ⎥ y(k) = C̅ ⎢ ⎢⎣ y(k) ⎥⎦

0.5(z + 1) (z − 1)2

Case 3: G (z ) =

(82)

where Ĉ is a P-diagonal matrix with the C̅ matrix sitting on its diagonal. Combining eq 72 and eq 82, the constraints on the output vector can be presented as

0.41(z − 0.86) (z − 1.12)2

(85)

0.53(z − 0.82) (z − 1.19)2

(86)

Case 4: G (z ) =

⎡ Δx (k)⎤ ⎡ Δx (k)⎤ m ̂ ΔU ≤ C4y − Ĉ ⎢ m ⎥ ⎥ ≤ CF C3ymin − Ĉ ⎢ max ⎢⎣ y(k) ⎥⎦ ⎢⎣ y(k) ⎥⎦

In Figure 3, we can easily find that the two methods achieve the set-point tracking and disturbance rejection successfully. From an overall perspective, the control performance of ENMSS model-based MPC is better than that of conventional NMSS model-based MPC. The overshoot and oscillations in the NMSS-based model are bigger than those in ENMSS model-based MPC, and at the same time, the tracking error in

(83)

where C3 and C4 are the corresponding matrices. By synthesizing eq 73, eq 80, and eq 83, the constraints on Z can be calculated. Then, the linear programming problem in eq 77 can be solved with simplex method. 5369


Article




The set-points of two MPC methods are both 1, and output disturbances with amplitudes of −0.05 and −0.1 are added to the two outputs, respectively. The control parameters of the two methods are the same with P = 8, M = 1, and α = 0.6. Here, ENMSS-based MPC has more degrees of freedom to choose its weighting matrices Q = diag(Q 1 ,Q 2 ,···Q 8 ), and Q j = diag(2,2,1,1,0,0,0,0,0,0,1,1)(j = 1,2,···8), L = diag(0.1,0.1) is chosen. For conventional NMSS model-based MPC, since its model is not expanded, it can only choose Qj = diag(1,1)(j = 1,2,···8) and L = diag(0.1,0.1). Here, model/plant mismatch cases are taken into account. Suppose the discrete process model of the actual process is as follows:

the ENMSS model-based MPC method is smaller than that of NMSS model-based MPC. The situation results in Figure 4 are the same as those in Figure 3, and the ensemble control performance based on ENMSS model-based MPC is improved. 3.4.2. MIMO Case Study. Here, the following discrete model of a glasshouse microclimate process is introduced.64 ⎡ y (k)⎤ ⎢ 1 ⎥ ⎢ y (k)⎥ ⎣ 2 ⎦ ⎡ 0.1237z−1 + 0.04935z−2 ⎤ − 0.2929z−2 ⎥⎡ ⎢ −1 ⎤ 1 − 0.8669z 1 − 0.8669z−1 ⎥⎢ u1(k) ⎥ ⎢ =⎢ −2 −3 −1 −2 ⎥⎢ u (k)⎥ 0.2933z + 0.1496z ⎥⎣ 2 ⎦ ⎢ − 0.05833z − 0.2214z ⎥⎦ ⎢⎣ 1 − 0.9001z−1 1 − 0.897z−1

(87)

The simulation procedure is as follows. 5370


Article

Industrial & Engineering Chemistry Research Case 5:

Then the incremental reference in eq 55 is rewritten as ⎡ (1 − α)(y − y(k)) ⎤ s ⎢ ⎥ ⎢ α(1 − α)(y − y(k)) ⎥ s ⎥ ΔR = ⎢ ⎢ ⎥ ⋮ ⎢ ⎥ ⎢⎣ α P − 1(1 − α)(ys − y(k))⎥⎦

⎡ y (k)⎤ ⎢ 1 ⎥ ⎢ y (k)⎥ ⎣ 2 ⎦ ⎡ − 0.3175z−2 0.1049z−1 + 0.05614z−2 ⎤ ⎥⎡ ⎢ ⎤ −1 1 − 0.7912z 1 − 0.7417z−1 ⎥⎢ u1(k) ⎥ ⎢ =⎢ −2 −3 −1 −2 ⎥⎢ u (k)⎥ 0.3028z + 0.1601z ⎥⎣ 2 ⎦ ⎢ − 0.05671z − 0.2553z ⎥⎦ ⎢⎣ 1 − 0.8123z−1 1 − 0.7735z−1

Then, from eq 60, eq 52, eq 90, and eq 93 define

(88)

In Figure 5, we can easily find that the ensemble control performance of ENMSS model-based MPC is better than that of NMSS model-based MPC. The overshoot and oscillations in the NMSS model-based MPC method are bigger than those of the ENMSS model-based MPC method. The tracking error in the ENMSS model-based MPC method is smaller.

kR′ = kR1(1 − α) + kR 2α(1 − α) + ··· + kRPα P − 1(1 − α) (94)

The control law is derived as (1 − z −1) L(z)U (z) = −(1 − z −1) P(z)Y (z) + kR′ (Ys(z) − Y (z))

4. INTERPRETATIONS OF ENMSS-BASED MPC A brief interpretation of ENMSSMPC and the relationship between ENMSSMPC and NMSSMPC are described in this section. 4.1. Transfer Functional Interpretations.54 The ENMSS model-based MPC (ENMSSMPC) design can be cast in terms of transfer functions as NMSS model-based MPC (NMSSMPC) in ref 64. The following will give a brief interpretation of ENMSSMPC. The transfer function is B (z) assumed to be n , and the model order for the input and

(95)

where Ys(z) and Y(z) are the z-transforms of y(k) and ys, respectively. The transfer function from the set-point to output is derived as T ̅ (z ) =

kR′ Bn(z) −1

(1 − z )( L(z)Ad (z) + P(z)Bn(z)) + kR′ Bn(z) (96)

with lim T̅(z) = 1

Ad (z)

(97)

z→1

output transfer function is n. Then, the dimensionality of z(k) is 2n. Write the vectors ks and kR in eq 60 as

showing that the closed-loop system tracks the constant setpoint without steady error. The transfer function from the output disturbance to output is derived as

ks = [k1k 2⋯knkn + 1kn + 2⋯k 2n] kR = [kR1kR 2⋯kRP]

(93)

(89)

Define two polynomial functions as follows:

S ̅ (z ) =

P(z) = k1 + k 2z −1 + k 3z −2 + ··· + knz −(n − 1)

(1 − z −1) L(z)Ad (z) −1

(1 − z )( L(z)Ad (z) + P(z)Bn(z)) + kR′ Bn(z) (98)

L(z) = 1 + kn + 1z −1 + kn + 2z −2 + ··· + k 2n − 1z −(n − 1)

with (90)

lim S ̅ (z) = 0

(99)

z→1

The reference trajectory is set as

The transfer function from the input disturbance to output response is derived as

r(k) = y(k) r(k + i) = α iy(k) + (1 − α i)ys i = 1, 2, ···, P

Si̅ (z) =

(91)

(1 − z −1) L(z)Bn(z) (1 − z −1)( L(z)Ad (z) + P(z)Bn(z)) + kR′ Bn(z) (100)

where ys is the set-point vector, and α is the smoothing factor. From eq 91, note that

with lim Si̅ (z) = 0

r(k) = y(k)

(101)

z→1

Δr(k + 1) = r(k + 1) − r(k)

It is seen that the closed-loop control system can reject constant output disturbances and input disturbance totally without steady error. Equation 60 can be further formulated as

= (1 − α)(ys − y(k)) Δr(k + 2) = r(k + 2) − r(k + 1)

n

= α(1 − α)(ys − y(k)) ⋮

j=1

∑

kjΔu(k − j + n)

j=n+1 P

Δr(k + P) = r(k + P) − r(k + P − 1) = α P − 1(1 − α)(ys − y(k))

2n − 1

Δu(k) = −∑ kjΔy(k − j + 1)− − k 2ne(k) − (92)

∑ kRjΔr(k + j) j=1

5371

(102) DOI: 10.1021/acs.iecr.7b01319 Ind. Eng. Chem. Res. 2017, 56, 5360−5394

Article

Industrial & Engineering Chemistry Research The control law includes the integral action of output tracking error, state feedback, and P -step set-points feed forward to improve the system’s dynamic performance. 4.2. Relationship between ENMSSMPC and NMSSMPC.52 Overall, the NMSSMPC method can be regarded as the special case of the ENMSSMPC method. The ENMSSMPC method will reduce to the NMSSMPC method by choosing Qj = diag{0,0,···,0,0,0,···0,qje2,qje2,···,qjeq} 1 ≤ j ≤ P; i.e., only the output errors are considered in the closed-loop system. At the same time, we use some variables notations from ref 64 to make the proof clear. If we only consider the outputs errors in the controller design, the cost function can be rewritten as P

J=

Through eq 50 and eq 52, it is easily seen that the following holds ⎡ Δx(k) ⎤ ⎡ Δx(k)⎤ ̑ (k) = CS̑ ⎢ ⎥ ⎥ = CS̑ ⎢ CSz ⎢⎣ e(k) ⎥⎦ ⎣⎢ y(k) − r(k)⎥⎦ ⎡ Δx(k)⎤ ⎡ 0 ⎤ ⎥ − CS̑ ⎢ = CS̑ ⎢ ⎥ ⎢⎣ y(k) ⎥⎦ ⎣ r(k)⎦

Considering the definition of C̑ in eq 106 and S in eq 107, we obtain the detailed elements in C̑ S as follows:

⎡ *1 ⎢ ⎢ *2 CS̑ = ⎢ ⎢⋮ ⎢ ⎢⎣*P

M

∑ e Τ(k + j)Q̅ je(k + j) + ∑ Δu Τ(k + j)LjΔu(k + j) j=1

j=1

(103)

where P is the prediction horizon, M is the control horizon, and Q̅ j = diag{qje1,qje2,···,qjeq} 1 ≤ j ≤ P, Lj(j = 1,2,···,M) are the weighting matrices for outputs errors and the control input. On the basis of eqs 51 and 52, the following holds e(k) = [ 02

⎡ Δx(k)⎤ Iq ]⎢ ⎥ ⎢⎣ e(k) ⎥⎦

⎡−Iq 0 ⎢ ⎢−Iq −Iq C̑θ = ⎢ ⎢ ⋮ ⋮ ⎢ ⎣−Iq −Iq

Combining eq 55 with eq 104, the future error vector can be predicted as follows: ̑ = CSz ̑ (k) + CF ̑ ΔU + C̑θ ΔR E = CZ (106) where C̑ is the block diagonal matrix with the matrix C̃ sitting on its diagonal, and

⎡ C 0 0 ⎢ C 0 ⎢ AC 2 ⎢ θ= AC AC C ⎢ ⋮ ⋮ ⋮ ⎢ ⎢ P−1 ⎣ A C AP − 2 C AP − 3C

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ AP − M B ⎦

0 0⎤ ⎥ 0 0⎥ 0 0⎥ ⎥ ⋱ ⋮⎥ ⎥ ⋯ C⎦

0 0

Τ

J = E Q̅ E + ΔU LΔU

0 ⎤ ⎥ 0 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ −Iq ⎦ 0

(113)

which can further obtain ⎡ ⎤ −Δr(k + 1) ⎢ ⎥ ⎢−Δr(k + 1) − Δr(k + 2)⎥ ⎢ ⎥ Cθ̑ ΔR = ⎢ ⋮ ⎥ ⎢ ⎥ P ⎢ ⎥ −∑ Δr(k + i) ⎢⎣ ⎥⎦ i=1

(114)

Refering to eq 49, the following equations can be easily proved r(k + 1) = r(k) + Δr(k + 1)

(107)

r(k + 2) = r(k + 1) + Δr(k + 2) = r(k) + Δr(k + 1) + Δr(k + 2)

The cost function in eq 103 can be presented as Τ

(112)

Combining eq 106 with eq 107, it is easily obtained that

(105)

⋯ ⋯

(111)

⎡ r(k)⎤ ⎢ ⎥ ⎡ 0 ⎤ ⎢ r(k)⎥ CS̑ ⎢ ⎥=⎢ ⎥ ⎣ r(k)⎦ ⎢ ⋮ ⎥ ⎢ r(k)⎥ ⎣ ⎦

In order to simplify the derivation process, we denote C̃ = [02 Iq] and the future error vector as

⎡ B 0 0 ⎡A⎤ ⎢ B 0 ⎢ 2⎥ ⎢ AB ⎢A ⎥ 2 ⎢ S = ⎢ ⎥; F = A B AB B ⎢ ⋮ ⎢ ⎥ ⋮ ⋮ ⎢ ⋮ ⎢⎣ AP ⎥⎦ ⎢ P−1 P 2 P − ⎣ A B A B A − 3B

Iq ⎤ ⎥ Iq ⎥ ⎥ ⋮⎥ ⎥ Iq ⎥⎦

where ∗i, (i = 1,2,···,P) are irrelevant elements for the proof. Then the following can be obtained

(104)

E = [e(k + 1), e(k + 2), ···, e(k + P)]Τ

(110)

⋮ P

(108)

r(k + P) = r(k) +

where Q̅ is the block diagonal matrix with the matrix Q̅ j sitting on its diagonal, and L is the block diagonal matrix with the Lj sitting on its diagonal. By substituting eq 106 into eq 108, the optimal control law can be derived as

∑ Δr(k + i) i=1

(115)

On the basis of eqs 110, 112, and 114 with eq 115, the following equation can be derived ̑ (k) + C̑θ ΔR = CSx ̑ (k ) − R CSz

T ̑ + L)‐1FT C̑ TQ̅ (CSz ̑ (k) + C̑θ ΔR ) ΔU = −(FT C̑ Q̅ CF

(116)

where

(109) 5372


Article

Industrial & Engineering Chemistry Research ⎡ C 0 0 ⎢ AC C 0 ⎢ θ = ⎢ A2 C AC C ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ P−1 P − 2 P ⎣ A C A C A − 3C

⎡ r(k + 1) ⎤ ⎢ ⎥ ⎡ Δx (k)⎤ ⎢ r(k + 2) ⎥ m ⎥; R = ⎢ x(k) = ⎢ ⎥ ⎢⎣ y(k) ⎥⎦ ⋮ ⎢ ⎥ ⎢ r(k + P)⎥ ⎣ ⎦

are the variables notations defined in ref 64. Substituting eq 116 into eq 109, the optimal control trajectory under the situation that only the outputs errors are considered is described as follows: T ̑T

‐1 T ̑ T

̑ + L) F C Q̅ (CSx ̑ (k ) − R ) ΔU = −(F C Q̅ CF

The cost function for ENMSSMPC is chosen as J(k) = ZΤQZ + ΔUΤLΔU

(117)

P

∑ [y(k + j) − r(k + j)]Τ Q̅ j[y(k + j) − r(k + j)]

u(k) = u(k − 1) + K p(k)(e1(k) − e1(k − 1))

j=1 M

+

(121)

where Q and L are the weighting matrices for the state variables and control inputs, respectively. Here, the incremental PID controller is chosen to combine with the ENMSSMPC method, and its formulation can be presented as follows:

which is just the optimal control law presented in ref 64 by choosing the following cost function J=

0 0⎤ ⎥ 0 0⎥ 0 0⎥ ⎥ ⋱ ⋮⎥ ⎥ ⋯ C⎦

+ K i(k)e1(k) + Kd(k)(e1(k) − 2e1(k − 1) + e1(k − 2))

∑ Δu Τ(k + j)LjΔu(k + j) (118)

j=1

e1(k) = yr (k) − y(k)

This concludes the proof.

where Kp(k), Ki(k), and Kd(k) are the proportional coefficient, integral coefficient, and differential coefficient at time instant k, respectively, and e1(k) is the error between the set-point and the actual output value at time instant k. In order to simplify the computation process, eq 122 can be rewritten as

5. PID CONTROL OPTIMIZED BY ENMSS MODEL PREDICTIVE CONTROL Compared with the MPC strategy, there is still room for the improvement of control performance of the PID control. Meanwhile, PID controllers possess a simpler structure. In this section, the PID control and ENMSSMPC are combined, and a novel PID controller that inherits the simple structure of conventional PID controllers and the same control performance as the proposed MPC algorithms is developed. The effectiveness of the proposed PID controller is examined via case studies on the temperature control of an industrial surfactant reactor and the control of a glasshouse microclimate process. 5.1. SISO PID Control Optimized by ENMSS Model Predictive Control. 5.1.1. PID Control Optimized by ENMSSMPC.57 Consider the difference equation model of the controlled system in eq 37

u(k) = u(k − 1) + w(k)Τ E(k) w(k) = [w1(k), w2(k), w3(k)]Τ w1(k) = K p(k) + K i(k) + Kd(k), w2(k) = −K p(k) − 2Kd(k), w3(k) = Kd(k) E(k) = [e1(k), e1(k − 1), e1(k − 2)]Τ

= L1u(k) + L 2u(k − 1) + ··· + Lnu(k − n + 1) (119)

w(k) = (− ((FΤQF + L)E(k)Τ E(k))−1FΤQ (Sz(k) + θ ΔR ))E(k) (124)

where y(k) and u(k) are the output and input of the process at time instant k, respectively. The following derivation process is the same as in Section 3.2.1, and then, we can obtain the future state variables from sampling time instant k

Then, Kp(k), Ki(k), and Kd(k) can be further obtained as K p(k) = −w2(k) − 2Kd(k) K i(k) = w1(k) − KP(k) − Kd(k)

(120)

Kd(k) = w3(k)

where ⎡ B 0 0 ⎡A⎤ ⎢ B 0 ⎢ 2⎥ ⎢ AB ⎢A ⎥ 2 ⎢ S = ⎢ ⎥; F = A B AB B ⎢ ⋮ ⎢ ⎥ ⋮ ⋮ ⎢ ⋮ ⎢⎣ AP ⎥⎦ ⎢ P−1 P 2 P − ⎣ A B A B A − 3B

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ AP − M B ⎦ ⋯ ⋯

(123)

For simplicity, the control horizon here is taken as 1, and by taking a derivative of the cost function J(k), the optimal control law is presented as

y(k + 1) = H1y(k) + H2y(k − 1) + ··· + Hmy(k − m + 1)


(122)

(125)

From eq 124, it is easily seen that w(k) will be infinite if E(k)TE(k) is approximating zero, showing that Kp(k), Ki(k), and Kd(k) will be infinite and is unrealistic in practice. Therefore, it is necessary to set a small permissible error limitation δ in which the parameters of the PID controller are unchanged and kept as the same values as the previous sampling instant. The details are shown as follows:

0 0

5373


Article

Industrial & Engineering Chemistry Research ⎧ K p(k) = K p(k − 1) ⎪ ⎪ ⎨ K i(k) = K i(k − 1)··· ···|e1(k)| ≤ δ ⎪ ⎪ K (k) = K (k − 1) ⎩ d d

u(k) = u(k − 1) + K p(k)(e1(k) − e1(k − 1)) + K i(k)e1(k) + Kd(k)(e1(k) − 2e1(k − 1) + e1(k − 2))

⎧ K p(k) = − w2(k) − 2Kd(k) ⎪ ⎪ ⎨ K i(k) = w1(k) − KP(k) − Kd(k)··· ···|e1(k)| > δ ⎪ ⎪ K (k) = w (k) ⎩ d 3

e1(k) = yr (k) − y(k)

where Kp(k), Ki(k), and Kd(k) are the proportional coefficient, integral coefficient, and differential coefficient at time instant k, respectively, and e1(k) is the error between the set-point value and actual output value at time instant k. To simplify the following derivation process, eq 132 can be reconstructed as

(126)

Then the optimal control input at sampling instant k can be calculated as u(k) = u(k − 1) + K p(k)(e1(k) − e1(k − 1))

u(k) = u(k − 1) + w(k)Τ E(k)

+ K i(k)e1(k) + Kd(k)(e1(k) − 2e1(k − 1) + e1(k − 2))

w(k) = [w1(k), w2(k), w3(k)]Τ (127)

w1(k) = K p(k) + K i(k) + Kd(k), w2(k)

58

5.1.2. PID Control Optimized by ENMSSPFC. Consider the same difference equation model as Section 5.1.1 and the same derivation procedure as Section 3.2.2, we can obtain the following equation first.

= −K p(k) − 2Kd(k), w3(k) = Kd(k) E(k) = [e1(k), e1(k − 1), e1(k − 2)]Τ

z(k + 1) = A mz(k) + Bmu(k) − Bmu(k − 1) + CmΔr(k + 1)

w(k) = ( −(FΤQFE(k)Τ E(k))−1FΤQ (Sz(k) + θ ΔR ))E(k) (134)


then it is easy to obtain the control parameters Kp(k), Ki(k), and Kd(k) K p(k) = −w2(k) − 2Kd(k) K i(k) = w1(k) − KP(k) − Kd(k)

(129)

On the basis of eq 128, it is easy to acquire the future state variables of PFC from sampling time instant k Z = Sz(k) + Fu(k) − Fu(k − 1) + θ ΔR

Kd(k) = w3(k)

(135)

Here, the permissible error limitation δ is also necessary, just like the situation in Section 5.1.1. Finally, the optimal control input at sampling instant k can be obtained as

(130)

where ⎡ Am ⎤ ⎡ Bm ⎤ ⎢ ⎥ ⎥ ⎢ ⎢ A m2 ⎥ ⎢ A mBm ⎥ S=⎢ ⎥; F = ⎢ ⎥ ⋮ ⎢ ⋮ ⎥ ⎥ ⎢ ⎢ P⎥ ⎢ A P − 1B ⎥ ⎣ m m⎦ ⎣ Am ⎦

u(k) = u(k − 1) + K p(k)(e1(k) − e1(k − 1)) + K i(k)e1(k) + Kd(k)(e1(k) − 2e1(k − 1) + e1(k − 2))

0 0 0 ⋱ ⋯

(136)

5.2. MIMO PID Control Optimized by ENMSS Model Predictive Control.59,60 5.2.1. PID Control Optimized by ENMSSMPC. Consider the difference equation process model with p inputs and q outputs and referring to the derivation procedure in Section 3.2.1, we can obtain the following future state variables from sampling time instant k first

0⎤ ⎥ 0⎥ ⎥ 0⎥ ⎥ ⋮⎥ Cm ⎥⎦


(137)

where

We choose the cost function of the ENMSSPFC method as follows: J(k) = ZΤQZ

(133)

From eqs 130, 131, and 133, the following holds (128)

Define

⎡ Cm 0 0 ⎢ 0 Cm ⎢ A mCm ⎢ 2 θ = ⎢ A m Cm A mCm Cm ⎢ ⋮ ⋮ ⋮ ⎢ ⎢ A P − 1C A P − 2 C A P − 3C ⎣ m m m m m m

(132)

⎡ B 0 0 ⎡A⎤ ⎢ B 0 ⎢ 2⎥ ⎢ AB ⎢A ⎥ 2 ⎢ S = ⎢ ⎥; F = A B AB B ⎢ ⋮ ⎢ ⎥ ⋮ ⋮ ⎢ ⋮ ⎢⎣ AP ⎥⎦ ⎢ P−1 P 2 P − ⎣ A B A B A − 3B

(131)

where Q is the weighting matrix for the state variables. Similar to Section 5.1.1, the incremental PID controller is introduced 5374

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ AP − M B ⎦ ⋯ ⋯

0 0


Article

Industrial & Engineering Chemistry Research ⎡ C 0 0 ⎢ AC C 0 ⎢ 2 θ=⎢ AC AC C ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ P−1 P − 2 P ⎣ A C A C A − 3C

0 0⎤ ⎥ 0 0⎥ 0 0⎥ ⎥ ⋱ ⋮⎥ ⎥ ⋯ C⎦

wi1(k) = kpi(k) + kii(k) + kdi(k) wi2(k) = −kpi(k) − 2kdi(k) wi3(k) = kdi(k)

On the basis of eqs 137−140 and taking a derivative of J(k), the following equation holds

The corresponding cost function is chosen as J(k) = ZΤQZ + ΔUΤLΔU

w(k) = E(k)(− ((FΤQF + L)E(k)Τ E(k))−1FΤQ (Sz(k) + θ ΔR )) (141)

(138)

where Q and L are the weighting matrices for the state variables and control inputs, respectively. The PID controller introduced is the incremental PID controller, and its control law is

then the coefficients for the proposed PID controller are calculated as kpi(k) = −wi2(k) − 2kdi(k) kii(k) = wi1(k) − kpi(k) − kdi(k)

u(k) = u(k − 1) + kp(k)(e1(k) − e1(k − 1)) + ki(k)e1(k) + kd(k)(e1(k) − 2e1(k − 1) + e1(k − 2))

kdi(k) = wi3(k)

It is easy to find that the w(k) will be infinite when the actual process output is close enough to the reference trajectory, and then, E(k) is close to 0. Such a situation will cause the coefficients of PID controller to reach infinity and is unrealistic in practice. Here, we need to set an error permission limitation δ to prevent such situation, and the rule for the setting of δ is

e1(k) = c(k) − y(k) e1(k) = [ e11(k) e12(k) ⋯ e1q(k)]Τ (139a)

wherekp(k), ki(k), and kd(k) are the proportional coefficient matrix, integral coefficient matrix, and differential coefficient matrix at time instant k, respectively, and e1(k) is the error vector between the set-point vector and the actual output vector at time instant k. Here, we adopt the case that p is equal or greater than q as an example. We assume that the first q inputs are paired with the q outputs for control, and the other p−q inputs are zeros. Also kp(k), ki(k), and kd(k) can be adopted as follows: ⎡ kn1(k) 0 ⎢ ⎢ 0 kn2(k) ⎢ 0 0 ⎢ kn(k) = ⎢ ⋮ ⋮ ⎢ ⎢ ⋯ ⎢ 0 ⎢ ⋮ ⋮ ⎢ ⎣ 0 0

0 0 ⋱ 0 0 ⋯ ⋯

⎧ kpi(k) = kpi(k − 1) ⎪ ⎪ ⎨ kii(k) = kii(k − 1)··· ···|e1i(k)| ≤ δ ⎪ ⎪ k (k) = k (k − 1) ⎩ di di ⎧ kpi(k) = − wi2(k) − 2kdi(k) ⎪ ⎪ ⎨ kii(k) = wi1(k) − kpi(k) − kdi(k)··· ···|e1i(k)| > δ ⎪ ⎪ k (k) = w (k) ⎩ di i3

⋯ ⋯ 0⎤ ⎥ ⋮ ⋯ ⋮ 0⎥ ⎥ 0 ⋱ ⋮ ⋮⎥ ⎥ knq(k) 0 ⋮ ⋮ ⎥ , (n = p , i , d) ⎥ 0 0 ⋯ 0⎥ ⋱ ⋮ ⋱ 0⎥ ⎥ ⋯ ⋯ ⋯ 0⎦ (139b) ⋯

z(k + 1) = A mz(k) + Bmu(k) − Bmu(k − 1) + CmΔr(k + 1)

(144)

Define ⎡ z(k + 1) ⎤ ⎡ Δr(k + 1) ⎤ ⎥ ⎢ ⎥ ⎢ ⎢ z(k + 2) ⎥ ⎢ Δr(k + 2) ⎥ Z=⎢ ⎥ ⎥ ; ΔR = ⎢ ⋮ ⋮ ⎥ ⎢ ⎥ ⎢ ⎢ z(k + P)⎥ ⎢ Δr(k + P)⎥ ⎦ ⎣ ⎦ ⎣

(140)

where ⎡ E (k ) 0 0 ⎢ 1 ⎢ 0 E 2 (k ) 0 E (k ) = ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ 0 Eq(k) ⎣ 0

(143)

5.2.2. PID Control Optimized by ENMSSPFC. The idea in Section 5.2.1 can be extended to the PFC case, and we can obtain the following new state space model first

The control law of the corresponding PID controller in eq 139a can be simplified further as follows: u(k) = u(k − 1) + E(k)Τ w(k)

(142)

⋯ 0⎤ ⎥ ⋯ 0⎥ ⎥ ⋮ ⋮⎥ ⎥ ⋯ 0⎦ 3q × p

(145)

Then, the future state variables of PFC from sampling time instant k can be described as Z = Sz(k) + Fu(k) − Fu(k − 1) + θ ΔR

(146)

where ⎡ Am ⎤ ⎡ Bm ⎤ ⎢ ⎥ ⎥ ⎢ ⎢ A m2 ⎥ ⎢ A mBm ⎥ S=⎢ ⎥; F = ⎢ ⎥ ⋮ ⎢ ⋮ ⎥ ⎥ ⎢ ⎢ P⎥ ⎢ A P − 1B ⎥ ⎣ A m m⎦ ⎣ m ⎦

Ei(k) = [ e1i(k) e1i(k − 1) e1i(k − 2)]Τ w(k) = [ w1(k) w2(k) ⋯ wq(k)]Τ

wi(k) = [ wi1(k) wi2(k) wi3(k)] 5375


Article



⎡ Cm 0 0 ⎢ 0 Cm ⎢ A mCm ⎢ θ = ⎢ A m 2 Cm A mCm Cm ⎢ ⋮ ⋮ ⋮ ⎢ ⎢ A P − 1C A P − 2 C A P − 3C ⎣ m m m m m m

0 0 0 ⋱ ⋯

0⎤ ⎥ 0⎥ ⎥ 0⎥ ⎥ ⋮⎥ Cm ⎥⎦

where ⎡ E (k ) 0 0 ⎢ 1 ⎢ 0 E 2 (k ) 0 E (k ) = ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ 0 Eq(k) ⎣ 0

The cost function is selected as Τ

J(k) = Z QZ

Ei(k) = [ e1i(k) e1i(k − 1) e1i(k − 2)]Τ (147)

w(k) = [ w1(k) w2(k) ⋯ wq(k)]Τ

where Q is the weighting matrix for the state variables. The conventional PID controller is chosen as the incremental PID controller, and its control law is presented as follows:

wi(k) = [ wi1(k) wi2(k) wi3(k)] wi1(k) = kpi(k) + kii(k) + kdi(k)

u(k) = u(k − 1) + kp(k)(e1(k) − e1(k − 1))

wi2(k) = −kpi(k) − 2kdi(k)

+ ki(k)e1(k) + kd(k)(e1(k) − 2e1(k − 1) + e1(k − 2))

wi3(k) = kdi(k)

By combining eqs 146−149 and taking a derivative of J(k), the following equation holds

e1(k) = c(k) − y(k) e1(k) = [ e11(k) e12(k) ⋯ e1q(k)]Τ

w(k) = E(k)( −(FΤQFE(k)Τ E(k))−1FΤQ (Sz(k) + θ ΔR ))

(148)

(150)

where kp(k), ki(k), and kd(k) are the block diagonal matrices in which the proportional coefficients, integral coefficients, and differential coefficients are sitting on their diagonals, respectively. In order to simplify the following computation process, eq 148 can be rewritten as u(k) = u(k − 1) + E(k)Τ w(k)

⋯ 0⎤ ⎥ ⋯ 0⎥ ⎥ ⋮ ⋮⎥ ⎥ ⋯ 0⎦ 3q × p

Then the corresponding PID control parameters are calculated as kpi(k) = −wi2(k) − 2kdi(k) kii(k) = wi1(k) − kpi(k) − kdi(k) kdi(k) = wi3(k)

(149) 5376


Article



Here, we set an error permission limitation ε to avoid the same situations mentioned in Section 5.2.1, and the details are described as follows:

(j = 1,2,···,15). The control parameters of the PID controller are Kp = 9.93, Ti = 1324, and Td = 13.85. Considering that the disturbances and uncertainties exist in

⎧ kpi(k) = kpi(k − 1) ⎪ ⎪ ⎨ kii(k) = kii(k − 1)··· ···|e1i(k)| ≤ ε ⎪ ⎪ k (k) = k (k − 1) ⎩ di di ⎧ kpi(k) = − wi2(k) − 2kdi(k) ⎪ ⎪ ⎨ kii(k) = wi1(k) − kpi(k) − kdi(k)··· ···|e1i(k)| > ε ⎪ ⎪ k (k) = w (k) ⎩ di i3

practice, here we suppose two cases of the actual process Case 6: G (s ) =

(152)

2.8e−28s 1310s + 1

(154)

2.5e−27s 1153s + 1

(155)

Case 7:

5.3. Simulation Examples. 5.3.1. SISO Case Study. Here, the temperature model of an industrial surfactant reactor is introduced,57 and its transfer function is as follows: G (s ) =

3.0e−30s 1441s + 1

G (s ) =

In Figure 6, it is easily seen that the responses of the ENMSSMPC method are smoother, and the IMC-PID method

(153)

shows the poor recovery ability under disturbance. From an

In this section, the PID controller tuned by the inner model control (IMC) method65 is introduced to be compared with the PID controller optimized by the ENMSSMPC method and the simulation procedure is as follows. The set-point of the controlled system is changed from 0 to 1 at time instant k = 0. After the temperature control system has reached its steady state, output disturbance with amplitude of −0.05 is added to the output for both methods. Here, the permissible error limitation δ is set as 10−5. The control parameters of the ENMSSMPC method are P = 15, M = 1, α = 0, L = 0.01, Q = diag(Q1,Q2,···,Q15), and Qj = diag(0,0,0,0,0,1),

overall perspective, the ensemble performance of the ENMSSMPC method is better than that of the IMC-PID method, although these methods can adapt to mismatched situations. In Figure 7, the situations are the same as those in Figure 6. In a word, the ENMSSMPC method provides improved control performance. 5.3.2. MIMO Case Study. Here, the same discrete model of a glasshouse microclimate process in Section 3.4.2 is introduced. 5377


Article


Figure 8. Closed-loop responses under case 8. ⎡ y (k)⎤ ⎢ 1 ⎥= ⎢ y (k)⎥ ⎣ 2 ⎦

In Figure 8, it is obvious that the ensemble control performance of the PID controller optimized by the ENMSSMPC method is better than that of the PID controller optimized by the NMSSMPC method. The overshoot and oscillations in the NMSSMPC tuning method are bigger, and the tracking error of the ENMSSMPC-based method is smaller. In a word, the PID controller optimized by the ENMSSMPC method provides improved control performance, although both methods achieve the set-point tracking and disturbance rejection successfully.

⎡ − 0.2929z −2 0.1237z −1 + 0.04935z−2 ⎤ ⎥⎡ ⎢ −1 ⎤ 1 − 0.8669z 1 − 0.8669z −1 ⎥⎢ u1(k) ⎥ ⎢ ⎥ ⎢ −2 −3 −1 −2 ⎢ ⎥ ⎢ − 0.05833z − 0.2214z 0.2933z + 0.1496z ⎥⎣ u 2(k)⎦ −1 −1 ⎥ ⎢⎣ ⎦ 1 − 0.9001z 1 − 0.897z (156)

In this section, the PID controller optimized by the NMSSMPC method is introduced to be compared with the PID controller optimized by the ENMSSMPC method, and the simulation procedure is the same as in Section 3.4.2. The control parameters of both methods are the same with P = 8, M = 1, L = diag(0.1,0.1), α = 0.6, Q = diag(Q1,Q2,···,Q8), Qj = diag(2,2,1,1,0,0,0,0,0,0,1,1)(j = 1,2,···8) for ENMSSMPC, and Qj = diag(1,1)(j = 1,2,···,8) for NMSSMPC. The permissible error limitation δ is set as 10−6. Here, we also consider the case of model/plant mismatch, and the actual discrete process model is as follows:

6. GENETIC ALGORITHM (GA) OPTIMIZATION61 In the previous section, it is obvious that the weighting matrix of state variables Qj plays a vital role on the control performance. After we choose the cost function of the controller system, we can use GA to optimize the selection of Qj. Here, our previous work61 is briefly mentioned. For the process described by eq 1, specify Qj as Q j = diag {qjx1, qjx 2 , ··· , qjxn , qju1, qju2 , ··· , qjud , qje}1 ≤ j ≤ P (158)

Case 8:

In order to obtain the desired control performance, it is important to make a compromise between the weighting factors in Qj. It is noted that the weighting factor on the tracking error qje can be set as a fixed value; thus, the optimization on the weighting factors of the state variables is the point needed to be focused on. Here, qje is chosen as 1. Also, the weighting factors on the past control input qju1, qju2,···,qjud are not taken into account. Thus, the main task is to optimize the weighting factors on output changes qj x1, qj x2,···,qj xn.

⎡ y (k)⎤ ⎢ 1 ⎥= ⎢ y (k)⎥ ⎣ 2 ⎦ ⎡ 0.1035z −1 + 0.05521z−2 ⎤ − 0.3056z −2 ⎥⎡ ⎢ ⎤ 1 − 0.7899z −1 1 − 0.7502z −1 ⎥⎢ u1(k) ⎥ ⎢ ⎢ −2 −3 −1 −2 ⎥⎢ u (k)⎥ ⎢ − 0.05611z − 0.2536z 0.3019z + 0.1587z ⎥⎣ 2 ⎦ ⎥⎦ ⎢⎣ 1 − 0.8232z −1 1 − 0.7813z −1 (157) 5378


Article


Figure 9. Process flow of coking furnace (101/3).

reproduction, which depends on their fitness values. Individuals with bigger fitness values have a bigger possibility to survive. In order to obtain the probability of an individual being selected, the Roulette wheel method is adopted, and the probability P(cl) is

6.1. Encoding Method. In this section, each element can be obtained by the binary encoding method. A 10-bit binary code is introduced as example as follows: |1|0|0|1|1|1|0|1|1|1|

(159)

The elements qj x1, qj x2,···,qj xn in eq 158 can be acquired through the following equations qji = max(qji) × b/210 (i = 1, 2, ···, n)

P(cl) =

(160)

ο(t ) + tr(t )

(161)

where ο(t) is the overshoot and tr(t) is the rise time. To simplify the optimization of GA, the following equation of the fitness of GA is defined. Max

1/[c + ο(t ) + tr(t )]

(163)

where f(cl) is the fitness value of individual cl, and N is the number of the individuals. 6.3.2. Crossover and Mutation Operators. As for the current individual cu and cv, the crossover operation and mutation probability pc and pm are carried out, and then, the offspring chromosomes c′u, c′v will be produced. Through the above several steps, a proper choice of Qj can be obtained.

6.2. Objective of Qj Optimization. In order to obtain a small overshoot and short rise time simultaneously, the following objective can be introduced Min

f (c l ) N ∑l = 1 f (cl)

7. INDUSTRIAL APPLICATION OF ENMSSMPC The fuels and petrochemical needs are essential parts of life. Nowadays, the acquirement of these necessities is greatly dependent on coking. In the coking process, the outlet temperature of the coke furnace influences the coking rate

(162)

where c is a constant. 6.3. Operators. 6.3.1. Selection Operator. Among the current population, a set of individuals should be selected for 5379


Article


change. Also, T and τ are calculated through several methods such as “time for 63.2% approach to new steady state”, “maximum slope method”, and “two-point method”. In this section, the “two-point method” is selected for deriving T and τ. Unlike the method in ref 68 in which two times are selected in terms of response times, here the two time points are selected based on the process output values, which will result in the simple formulas for the parameter calculation. The calculation process is as follows. Define y(t) as the outlet temperature variable with a steady state value y(∞) and denote U0 as the step input. Then, the steady state process gain is K = y(∞)/U0. Here, y(t) can be calculated from the FOPDT model in eq 164 as

and the quality of the final products directly. During the operation process, process dynamics, uncertainties of the coking kinetics, and other situations may impact the control performance of the relevant controllers. The employment of the conventional PID control on such a process may be unsatisfactory. The proposed MPC has been applied to the temperature control and oxygen content control, and improved results are obtained.62,63 In this section, the industrial application of the ENMSSMPC strategy on the temperature control of the coke furnace is presented as an example. Through the input−output process data, the extended state space model is obtained, and then, the relevant control scheme is derived. The application results show the effectiveness of the proposed MPC approach. 7.1. Plant Description. The detailed framework of the coke furnace (F101/3) is described in Figure 9. The residual oil is divided into two flows (FRC8103/4) and sent to the convection room for preheating at first. Then, they mix together and flow into the fractionating tower (T102) to exchange heat with the gas oil in it. After the heat exchange process, the heavy part of them (called circulating oil) is pumped into two flows (FRC8107/8) and sent to the radiation room. In the radiation room, they will be controlled at 495 °C, and then, they mix together and go to the coke towers (T101/ 5,6) for coke removing. During the operation, one of the coke towers works until it is full, and then, it will be replaced by the other. This is known as coke towers switches. In the process, the radiation room outlet temperature variables (TRC8103/5) are the controlled variables, and the fuel volumes supplied to the furnace are the manipulated variables. Due to the uncertainties and disturbances, the relationship between them is difficult to be acquired. The temperature is influenced by the heat exchange in the fractionating tower (T102) because the circulating oil volume changes significantly. The gas oil volume is affected by coking of residues in the coke towers, and then, there is an impact on the temperature. Besides, change of the circulating oil volume, which is influenced by the coke towers switches and the random switch time of the coke towers, also affects the radiation room outlet temperature. For conventional methods, it is hard to obtain the desired control performance. Through the modules in the control station based on the CENTUM CS3000 control system, the original PID controllers are designed. It is a fact that the PID control can reject the abrupt disturbance rapidly; thus, the master−slave architecture will be adopted in the process. The set-point of the PID controller is provided by MPC. The goal of the master−slave architecture is to maintain the outlet temperature at a fixed value. 7.2. Modeling. Considering all the aforementioned uncertainties and disturbances, the model will be complicated, and the relevant controller design will be difficult. Hence, a simple model including the main process dynamics is necessary. Here, to be consistent with the experimental test, the first-order plus dead time (FOPDT) model is chosen Ke−τs G (s ) = Ts + 1

⎧0 t t1 > τ, we can obtain the following equations y(t1) = K − Ke−t1− τ / T t ≥ τ y(t 2) = K − Ke−t2 − τ / T t ≥ τ

(166)

Then, the following holds T=

t 2 − t1 ln[1 − y(t1)/K ] − ln[1 − y(t 2)/K ]

τ=

t 2 ln[1 − y(t1)/K ] − t1ln[1 − y(t 2)/K ] ln[1 − y(t1)/K ] − ln[1 − y(t 2)/K ]

(167)

Here, y(t1) = 0.39K and y(t2) = 0.63K are chosen; thus, T and τ are obtained as T = 2(t 2 − t1) τ = 2t1 − t 2

(168)

7.3. Controller Design. Under sampling time Ts, the discrete form of eq 164 is y(k + 1) + Fy(k) = Hu(k − L)

(169)

where y(k+1), y(k), and u(k−L) are the output variables and input variable at corresponding time instants, respectively. Also, k is the current time instant, and F, H, and L are model coefficients F = −exp( −Ts/T ) H = K (1 + F ) L = τ /Ts

(170)

By adding a difference operator Δ = 1 − z−1 to eq 169, the following formula is obtained Δy(k + 1) + F Δy(k) = H Δu(k − L)

(171)

Define a new state space vector as Δxm(k)T = [Δy(k), Δu(k − 1), Δu(k − 2), ··· , Δu(k − L)]

(164)

(172)

where K is the steady state process gain, T is the time constant of the process, and τ is time delay of the process. The identification of the three parameters has been mentioned in many text books.67,68 Here, K is obtained by the division of the steady value of process output with the input

Then, the following state space model is obtained Δxm(k + 1) = A mΔxm(k) + BmΔu(k) Δy(k + 1) = CmΔxm(k + 1) 5380


Article

Industrial & Engineering Chemistry Research ⎡ Δr(k + 1) ⎤ ⎥ ⎢ ⎢ Δr(k + 2) ⎥ ΔR = ⎢ ⎥ ⋮ ⎥ ⎢ ⎢ Δr(k + N )⎥ y ⎦ ⎣

where ⎡− F ⎢ ⎢ 0 Am = ⎢ 0 ⎢ ⎢ ⋮ ⎢⎣ 0

0 ⋯ 1 ⋮

⎡0⎤ H⎤ ⎥ ⎢ ⎥ 0⎥ ⎢1 ⎥ 0 ⎥ , Bm = ⎢ 0 ⎥ , Cm = [1 0 0 ⋯ 0] ⎥ ⎢ ⎥ ⋮⎥ ⎢⋮⎥ ⎢⎣ 0 ⎥⎦ ⎥ 1 0⎦ (174)

⋯ ⋯ 0 ⋱

0 ⋯ ⋯ ⋱

⋯ 0

On the basis of eq 178, the process output predictions can be described as

Denote the output error as

Z = Fz(k) + ΦΔU + SΔR

e(k) = y(k) − r(k)

(175)

⎡ B 0 0 ⎢ 0 B ⎢ AB Φ = ⎢ A2 B AB B ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ Ny − 1 Ny − 2 Ny − 3 ⎣A B A B A B ⎡ A⎤ ⎢ 2⎥ ⎢A ⎥ =⎢ ⎥ ⋮ ⎢ ⎥ ⎢⎣ ANy ⎥⎦

e(k + 1) = e(k) + CmA mΔxm(k) + CmBmΔu(k) − Δr(k + 1)

(176)

Define the new state space vector as ⎡ Δx (k)⎤ m ⎥ z(k ) = ⎢ ⎢⎣ e(k) ⎥⎦

(177)

Then, the ENMSS model is obtained z(k + 1) = Az(k) + BΔu(k) + C Δr(k + 1)

(179)

(185)

Q = block diag {Q 1 , Q 2 , ⋯ , Q N } y

L = block diag {L1 , L 2 , ⋯ , L Nu}

j=1

(186)

7.4. Experiment Results. Experiments are done using CENTUM CS3000 Distributed Control System (DCS) shown in Figure 10. Industrial field control stations (FCS), an engineering station (EWS), two human interface stations (HIS), a monitor computer (PC), a printer (PRT), and many I/O modules compose the DCS control system. Using the process data, the generalized process can be modeled by a SISO model, where the outlet temperature set-

s.t. j ≥ Nu

(180)

where Ny and Nu are the prediction horizon and the control horizon, respectively, and Ny ≥ Nu. Here, z(k + j) is the state prediction for time k + j made at time k, Lj ≥ 0 is the weighting factor on control action, and Qj is the symmetrical weighting matrix Q j = diag {qjy , qju1 , ···, qjuL , qje}, 1 ≤ j ≤ Ny

(184)

where

∑ z T(k + j)Q jz(k + j) + ∑ ΔuT(k + j − 1)LjΔu(k + j − 1)

Δu(k + j) = 0

0 0⎤ ⎥ 0 0⎥ 0 0⎥ ⎥ ⋱ ⋮⎥ ⎥ ⋯ C⎦

ΔU = −(ΦTQ Φ + L)−1ΦTQ (Fz(k) + SΔR )

Nu

j=1

0 0

Synthesizing eq 180 and eq 183, the control law can be obtained by the minimization of J

Here, 0 in eq 179 is a zero vector with appropriate dimensions. Remark 2. In traditional state space MPC, the future output prediction is obtained based on eq 173. The proposed MPC adopts eq 177 to acquire the future prediction and design the controller. The cost function is chosen as J=

⎤ ⎥ ⎥ 0 ⎥; F ⋯ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ ANy − Nu B ⎦ ⋯ ⋯

⎡ C 0 0 ⎢ 0 AC C ⎢ 2 ⎢ S= AC AC C ⎢ ⋮ ⋮ ⎢ ⋮ ⎢ Ny − 1 N 2 N − ⎣A C A y C A y − 3C

(178)

where ⎡ Am 0⎤ ⎡ Bm ⎤ ⎡0⎤ ⎥, B = ⎢ ⎥, C = ⎢ ⎥ A=⎢ ⎣−1⎦ ⎢⎣CmBm ⎥⎦ ⎢⎣Cm A m 1 ⎥⎦

(183)

where

Combining eq 173 and eq 175, we obtain

Ny

(182)

(181)

Remark 3. It is obvious that the output error and the state variables can both be regulated by tuning qje and qjy1, qj u1,···,qjuL, respectively. Define ⎡ z(k + 1) ⎤ ⎡ ⎤ Δu(k) ⎢ ⎥ ⎢ ⎥ ⎢ z(k + 2) ⎥ ⎢ Δu(k + 1) ⎥ Z=⎢ ⎥ , ΔU = ⎢ ⎥, ⋮ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ z(k + N )⎥ ⎢ Δu(k + N − 1)⎥ ⎣ ⎦ y ⎦ ⎣ u

Figure 10. Automation system of the outlet temperature system. 5381


Article

Industrial & Engineering Chemistry Research point of the PID is the input, and the measured outlet temperature is the output. By testing the step response of the generalized process around an operation point, the SISO model is obtained. At the same time, the PID parameters Kc = 70, TI = 120, and TD = 20 are taken into account. In Figure 11, we can see that there is a data database (π database) to store the process data in DCS. Here, the outlet

Table 1. Transfer Function Models group 1

group 2

A: G1(s) =

1.15e−120s 400s + 1

C: G3(s) =

0.96e−140s 500s + 1

B: G2(s) =

0.85e−150s 600s + 1

D: G4(s) =

1.1e−150s 300s + 1

Since the time delays of G2(s) and G4(s) are the largest, the model can be chosen from one of them to include the delays. Meanwhile, the time constants of G1(s), G3(s), and G4(s) are closer to one another; thus, G4(s) is chosen as the process model. Under sampling time Ts = 30 s, G4(s) is first transformed into discrete model as 0.1047 GTs(z) = 6 (188) z − 0.9048z 5 Figure 11. Data acquisition configuration for outlet temperature system.

which can then be transformed into the input−output formulation shown in eq 169. The proposed MPC has been applied to the outlet radiation temperature control of the coke furnace, where G4(s) and a sampling time of Ts = 30 s are used. Figure 13 shows the performance of set-point tracking from 493 to 495 °C. The lower graph shows the corresponding input

temperature Tset, measured outlet temperature Tm, and input fuel flow Fflow are collected. The response of the process is acquired by both increasing and decreasing the set-point of the PID. Here, the sampling time is 10 s. The filtered and modeling results with the setpoint change cycle (495 °C → 497 °C → 495 °C → 493 °C → 495 °C) and the real-time responses are shown in Figure 12. It

Figure 13. Performance of the closed-loop system for tracking.

is obvious that the dynamic response of the process is complex, and it is difficult to obtain an accurate model. Considering the noise effect, the process data were first filtered, and the low-pass filter is chosen as 0.15 Gf (z) = (187) z − 0.85

fuel flow. Also, the controller shows a fast response, which is approximately 25 min, whereas Figure 12 shows that the response of the original PID is more than 33 min. To test the load disturbance rejection, the residual oil flow into the coke furnace was changed from 46 t/h to about 44 t/h at t = 10 min. In Figure 14, it shows that the proposed MPC successfully maintained the required set-point. Two groups of regulatory control results are shown in Figures 15 and 16. The temperature is more stable compared with before. The corresponding statistical results concerning the average value, the maximum, minimum, standard deviation, and the extreme deviation of the outlet temperature are shown in Table 2, showing a more steady control by the proposed approach.

In Table 1, the transfer functions are listed. We can see that the tests result in different process models. Even for the two forward step-response tests “A” and “D” or the reverse stepresponse tests “B” and “C”, the dynamics differ a lot.

8. NEW INSIGHTS INTO THE PROPOSED METHOD In order to explore the essence of the proposed method further, the extended disturbance model-based MPC66 in which the

Figure 12. Real-time forward and reverse step-response tests.

5382


Article

Industrial & Engineering Chemistry Research Table 2. Statistical Results group

item

average

max

min

standard deviation

extreme deviation

1 1 2 2

PID proposed PID proposed

495.6 495.2 495.5 495.2

497.6 497.1 497.4 496.1

492.1 493.7 492.1 493.1

0.8474 0.4806 0.7661 0.4176

5.5 3.4 5.3 3.0

extended state space model-based PFC schemes. Here, we assume that the state is measurable and there is no noise. The relevant state space model is chosen as xm(k + 1) = A mxm(k) + Bmu(k − 2) ym (k + 1) = Cmxm(k + 1) Figure 14. Residual oil flow load change test.

(189)

where ⎡1.1053 ⎡ 1 ⎤ 0 ⎤ Am = ⎢ ⎥ ; Bm = ⎢⎣ ⎥ ; C = [0 1] ⎣ −0.01 0.8186 ⎦ 0.0858 ⎦ m

Select the new state variable as x(k) = [xm(k)T, u(k − 1), u(k − 2)]T, then the corresponding state space model is x(k + 1) = Ax(k) + Bu(k) ym (k + 1) = Cx(k + 1)

(190)

where ⎡ A m 0 Bm ⎤ ⎡0⎤ ⎢ ⎥ ⎢ ⎥ A = ⎢ 0 0 0 ⎥ ; B = ⎢ 1 ⎥ ; C = [Cm 0 0 ] ⎢⎣ 0 1 0 ⎥⎦ ⎣0⎦

On the basis of the aforementioned assumption, the disturbance integrated in the state variable is estimated as follows:66

Figure 15. Experimental results of PID and the proposed MPC: Group 1.

d(k) = d(k − 1) + y(k) − ym (k)

(191)

where d(k) is the integrating disturbance at time instant k, and y(k) and ym(k) are the actual process output and the model output at time instant k, respectively. On the basis of eqs 190 and 191, the disturbance model in which the output disturbance is considered can be expressed as66 ⎡ x(k + 1)⎤ ⎡ A 0 ⎤⎡ x(k)⎤ ⎡ B ⎤ ⎥ + ⎢ ⎥u(k) ⎢ ⎥=⎢ ⎥⎢ ⎢⎣ d(k + 1)⎥⎦ ⎣ 0 1 ⎦⎢⎣ d(k)⎥⎦ ⎣ 0 ⎦ ⎡ x(k + 1)⎤ ⎥ ym (k + 1) = [C 1]⎢ ⎢⎣ d(k + 1)⎥⎦

(192)

Refering to the procedure in Section 2.1, the predicted output can be derived based on the disturbance model in eq 192. The relevant cost function is selected as

Figure 16. Experimental results of PID and the proposed MPC: Group 2.

J(k) = (Y − S)Τ Q (Y − S)

disturbance is integrated in the state variable to acquire the offset-free control is introduced as the comparison. In the following, PFC is designed based on the proposed extended state space models and disturbance models to show the relationships between the two strategies. 8.1. Case Study on the State Space Model. In this section, the state space model is introduced to test the two

(193)

where Y is the predictive output vector, and S is the corresponding reference trajectory. By minimizing the objective function J(k), we can obtain the optimal control law u(k). Our proposed PFC strategy based on the proposed extended state space model is shown in Section 2.2.2. 5383


Article





For both approaches, the prediction horizon P is chosen as 20, and the smoothing factor of the reference trajectory is 0. Under the model/plant matched case, the weighting matrix Q for the proposed strategy is Q = diag(Q1,Q2,···,Qp), Qj = diag(0,0,0,0,1), (j = 1,2,···,P), i.e., only the tracking error is weighted. In the model/plant mismatched cases, Q is selected as Q = diag(Q1,Q2,···,Qp), Qj = diag(0,20,0,0,1), (j = 1,2,···,P) for the proposed scheme. Compared with the proposed approach, the design degrees of freedom of the disturbance model-based PFC strategy is fewer, and the relevant Q can only be adopted as Q = diag(Q1,Q2,···,Qp), Qj = 1,(j = 1,2,···,P) for all cases. In both cases, the set-point is 1, and a step disturbance with amplitude of −0.2 is added to the process output at time instant k = 150. Besides, a random white noise sequence with standard deviation 0.2 is added to the process output under these mismatched cases to test the control performance of both methods further. Figure 17 shows the responses for both methods under the model/plant match case. It is obvious that their responses are the same, which implies that the two approaches are the same when only the tracking error is weighted in the proposed

In order to verify the control performance of both methods, the model/plant match case (here referred to as Case 9) and three mismatched cases generated by the Monte Carlo simulation are introduced. The three model/plant mismatched cases are as follows: Case 10: ⎡ 1.1053 ⎡ 1.1 ⎤ 0 ⎤ Am = ⎢ ⎥ ; Bm = ⎢ ⎥ ; C = [0 1] ⎣−0.011 0.9173⎦ ⎣ 0.0736 ⎦ m

Case 11: ⎡ 1.1053 ⎡ 1.2 ⎤ 0 ⎤ Am = ⎢ ⎥ ; Bm = ⎢⎣ ⎥ ; C = [0 1] ⎣−0.009 0.8729 ⎦ 0.0911⎦ m

Case 12: ⎡ 0.9 ⎤ ⎡ 1.1053 0 ⎤ ; Bm = ⎢ Am = ⎢ ⎥ ; C = [0 1] ⎥ ⎣−0.012 0.8845⎦ ⎣ 0.0943⎦ m 5384


Article



Figure 21. Closed-loop responses under white noise for case 10.



can easily find that the oscillations of the proposed method are

strategy. Figures 18−20 show the responses of both schemes under model/plant mismatched cases. It can be easily seen that the response of the proposed method is smoother. The overshoot and oscillations of the disturbance model-based PFC strategy is bigger. From an overall perspective, the ensemble control performance of the proposed scheme is better. The responses under model/plant mismatched cases and white noise for the two approaches are shown in Figures 21−23, and the relevant standard deviation values are listed in Table 3. We

smaller, which proves the aforementioned viewpoints further. 8.2. Case Study on the Input−Output Model Based Process. In this section, we introduce an input−output model based process to evaluate the control performance of the two strategies further. The corresponding input−output model is chosen as follows: 5385


Article


test the control performance of the two PFC strategies. The three model/plant mismatched cases are as follows:

Table 3. Standard Deviation Values under Model/Plant Mismatched for Case of Section 8.1 items y u y u y u

case 1 case 2 case 3

G (z ) =

proposed

disturbance model

0.2175 0.1356 0.2154 0.1348 0.2163 0.1350

0.2245 0.2264 0.2214 0.2245 0.2232 0.2252

−0.002497 z 3 − 0.8752z 2

Case 14: G (z ) =

−0.001653z − 0.001878 z4 − 0.8642z 3

Case 15: G (z ) =

(194)

−0.002192z − 0.001999 z 3 − 0.8324z 2

Choose the state vector as xm(k)=[y(k) and u(k−1),u(k− 2)]T, then the relevant state space model can be described as Case 16:

xm(k + 1) = A mxm(k) + Bmu(k) y(k + 1) = Cmxm(k + 1)

G (z ) =

(195)

where

For both methods, the prediction horizon is selected as 20, and the smoothing factor of the reference trajectory is 0.65. Under the model/plant matched case, only the tracking error is weighted for the proposed scheme, and the relevant weighting matrix is Q = diag(Q1,Q2,···,QP), Qj = diag(0,0,0,1), (j = 1,2,···,P). For the model/plant mismatched cases, Q is chosen as Q = diag(Q1,Q2,···,QP), Qj = diag(500,0,0,1), (j = 1,2,···,P) for the proposed strategy. As to the disturbance model-based PFC approach, the weighting matrix Q is selected as Q = diag(Q1,Q2,···,QP), Qj = 1,(j = 1,2,···,P) for all cases. In all cases, the set-point is 1, and a step disturbance with amplitude of −0.1 is added to the output at time instant k = 100. Further, a random white noise sequence with standard deviation 0.2 is added to the process output to verify the relevant control performance of both strategies under the model/plant mismatched cases. Figure 24 shows the responses of both methods under model/plant match case. It can be easily found that their responses are the same. In other words, when only the tracking error is weighted in the proposed strategy, it is equivalent to the disturbance model-based PFC scheme. Figures 25−27 show the responses of two approaches under model/plant mismatched cases. Overall, the ensemble control performance of the proposed PFC method is better. The overshoot and oscillations of the disturbance model-based PFC approach is bigger. The responses of the proposed strategy are smoother. Figures 28−30 show the responses under mismatched cases and random white noise for both methods. The corresponding standard deviation values are listed in Table 4. It is obvious that

⎡ 0.8752 0 −0.002497 ⎤ ⎡0⎤ ⎡1 ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ; ; B = C = Am = ⎢ 0 1 0 0 ⎥ m ⎢ ⎥ m ⎢0⎥ ⎣0⎦ ⎣0⎦ ⎣ 0 ⎦ 1 0

Here, we assume that there are no noises, and the calculation formulas of the integrating disturbance are the same as eq 191. Then, the relevant disturbance model in which just the output disturbance is considered can be obtained. x(k + 1) = Ax(k) + Bu(k) y(k + 1) = Cx(k + 1)

(196)

where x(k) = [y(k), u(k − 1), u(k − 2), d(k)]Τ ⎡ 0.8752 ⎢ 0 A=⎢ ⎢ 0 ⎢⎣ 0

0 −0.002497 0 0 1 0 0 0

−0.0006451z − 0.000916 z4 − 0.8959z 3

⎡0⎤ ⎡1 ⎤ 1⎤ ⎥ ⎢ ⎥ ⎢ ⎥ 0 1 0⎥ ; B = ⎢ ⎥; C = ⎢ ⎥ ⎢0⎥ ⎢0⎥ 0⎥ ⎥ ⎢⎣ ⎥⎦ ⎢⎣ ⎥⎦ 0 0 1⎦

On the basis of eq 196, the relevant predicted output can be derived. Choosing the cost function as eq 193, the corresponding optimal control law u(k) can be acquired finally. Note that more degrees of freedom can be provided for the proposed controller design, and Section 3.2.2 shows the details. Similar to Section 8.1, the model/plant match case (referred to as Case 13) and three mismatched cases are introduced to

Figure 24. Closed-loop responses under case 13. 5386


Article






the oscillations of the proposed scheme are smaller, which proves the effectiveness of the proposed approach further. 8.3. Further Case Study on the State Space Model with the Role of Kalman Filter. In order to explore the role of Kalman filter and the relationships between the proposed method and that with Kalman filter, here the process and the Kalman filter-based MPC in ref 66 are introduced. The corresponding state space model is

x(k + 1) = Ax(k) + Bu(k) y(k + 1) = Cx(k + 1)

(197)

where 5387


Article




⎡1 0⎤ ⎢ ⎥ Cd = ⎢ 0 0 ⎥ ⎣0 1⎦

Table 4. Standard Deviation Values under Model/Plant Mismatched for Case of Section 8.2 items case 1 case 2 case 3

y u y u y u

proposed

disturbance model

0.4266 198.68 0.5629 237.57 0.3791 379.16

0.4861 329.60 0.8656 395.29 0.5633 1069.1

The gain of the Kalman filter is chosen as in ref 66 ⎡ 0 1.734 × 10−2 ⎢ 0.9934 ⎢0 ⎢ K g = ⎢ 0 −8.770 × 10−2 ⎢ 1 −1.734 × 10−2 ⎢ ⎢⎣ 0 8.770 × 10−2

⎡ 0.2511 −3.368 × 10−3 −7.056 × 10−4 ⎤ ⎢ ⎥ A = ⎢ 11.06 −2.545 0.3296 ⎥ ⎢⎣ 0 ⎥⎦ 0 1

(199)

Thus, the state and disturbance are estimated by the following formula

⎡−5.246 × 10−3 1.530 × 10−5 ⎤ ⎡1 0 0⎤ ⎢ ⎥ ⎢ ⎥ B=⎢ 1.297 0.1218 ⎥; C = ⎢ 0 1 0 ⎥ ⎢ ⎥ ⎣0 0 1⎦ ⎣ −6.592 × 10−2 ⎦ 0

⎡ x(̂ k|k)⎤ ⎡ x(̂ k|k − 1)⎤ ⎢ ⎥=⎢ ⎥ + K g (y(k) − Cx(̂ k|k − 1) ⎢⎣ d(̂ k|k)⎥⎦ ⎢⎣ d(̂ k|k − 1)⎥⎦ − Cdd(̂ k|k − 1))

According to ref 66, only the first and third outputs are

(200)

Finally, the state prediction is

required to be controlled at the set-point, then the relevant

⎡ x(̂ k + 1|k)⎤ ⎡ ⎤⎡ x(̂ k|k)⎤ ⎡ B ⎤ ⎢ ⎥ = ⎢ A 0 ⎥⎢ ⎥ + ⎢ ⎥u(k) ⎢⎣ d(̂ k + 1|k)⎥⎦ ⎣ 0 I ⎦⎢⎣ d(̂ k|k)⎥⎦ ⎣ 0 ⎦

disturbance model in which just the output disturbance is considered is obtained. ⎡ x(k + 1)⎤ ⎡ A 0 ⎤⎡ x(k)⎤ ⎡ B ⎤ ⎥ + ⎢ ⎥u(k) ⎢ ⎥=⎢ ⎥⎢ ⎢⎣ d(k + 1)⎥⎦ ⎣ 0 I ⎦⎢⎣ d(k)⎥⎦ ⎣ 0 ⎦ ⎡ x(k + 1)⎤ ⎥ ym (k + 1) = [C Cd ]⎢ ⎢⎣ d(k + 1)⎥⎦

0⎤ ⎥ 0⎥ ⎥ 0⎥ 0⎥ ⎥ 1 ⎥⎦

(201)

Select the same cost function as eq 193, then the optimal control law u(k) of the disturbance model-based PFC approach can be obtained. Refering to Section 2.2.2, the relevant optimal control law of the proposed strategy is acquired (the extended state vector is xe(k) = [x(k)T,e1(k),e2(k),e3(k)]T). In this section, the model/plant match case (referred to as case 17) and three mismatched cases (referred to as cases 18, 19, and 20) are generated to verify the control performance of both methods.

(198)

where 5388


Article




5389


Article




5390


Article



Case 18:

Case 20: ⎡ 0.3192 −2.949 × 10−3 −6.519 × 10−4 ⎤ ⎢ ⎥ A = ⎢ 12.79 −2.496 0.3455 ⎥ ⎢⎣ 0 ⎥⎦ 0 1

⎡ 0.2965 −3.079 × 10−3 −7.423 × 10−4 ⎤ ⎢ ⎥ A = ⎢ 13.34 −2.274 0.3081 ⎥ ⎢⎣ 0 ⎥⎦ 0 1

⎡−4.728 × 10−3 1.103 × 10−5 ⎤ ⎡1 0 0⎤ ⎢ ⎥ ⎢ ⎥ = B=⎢ ; C 1.006 0.1471 ⎥ ⎢0 1 0⎥ ⎢ ⎥ ⎣0 0 1⎦ ⎣ −5.879 × 10−2 ⎦ 0

⎡−6.314 × 10−3 1.189 × 10−5 ⎤ ⎡1 0 0⎤ ⎢ ⎥ ⎢ ⎥ B=⎢ 1.021 0.1304 ⎥; C = ⎢ 0 1 0 ⎥ ⎢ ⎥ ⎣0 0 1⎦ ⎣ −6.834 × 10−2 ⎦ 0

For both approaches, the prediction horizon is chosen as 10, and the smoothing factor of the reference trajectory is 0. As to the model/plant matched case, there is no need to use the extra freedom, so the weighting matrix Q = diag(0,0,0,1,0,1) is chosen for the proposed approach. For the three mismatched cases, the extra freedom can be employed to improve the control performance, and the relevant weighting matrix is Q = diag(50,0,0,1,0,1) in the proposed method. At the same time, Q = diag(1,0,1) is selected for the disturbance model-based PFC method. The set-point for the two outputs are 1, and a step disturbance with amplitude of −0.1 is added to the controlled outputs at time instant k = 80 in all cases. Besides, a random white noise sequence with standard deviation 0.01 is added to the controlled outputs to test the control performance of both strategies further. The responses of the two approaches under the model/plant match case are shown in Figure 31. It is obvious that these responses are nearly the same, and their control performances are both satisfactory. Figures 32−34 show the closed-loop responses of both methods for three model/plant mismatched cases. We can easily find that both approaches achieve the setpoint tracking and disturbance rejection successfully. In Figure

Case 19:

⎡ 0.2480 −3.812 × 10−3 −8.191 × 10−4 ⎤ ⎢ ⎥ A = ⎢ 9.21 −2.869 0.2920 ⎥ ⎢⎣ 0 ⎥⎦ 0 1 ⎡−6.314 × 10−3 1.727 × 10−5 ⎤ ⎡1 0 0⎤ ⎢ ⎥ ⎢ ⎥ = B=⎢ ; C 1.014 0.1163 ⎥ ⎢0 1 0⎥ ⎢ ⎥ ⎣0 0 1⎦ ⎣ −7.186 × 10−2 ⎦ 0

5391


Article




32, the overshoot of the first output of the disturbance modelbased PFC strategy is bigger than that of the proposed method. Note that the coefficients in the third row of matrix A are

special (i.e., 0, 0, and 1); thus, there are no mismatches in A for the third output. Thus, responses of the third output of both methods are almost the same. In Figures 33 and 34, the 5392


Article

Industrial & Engineering Chemistry Research Table 5. Statistical Results of y1 for Cases under White Noise item case 1 case 2 case 3

proposed disturbance model proposed disturbance model proposed disturbance model

standard deviation

max

min

extreme deviation

10−5 10−5 10−5 10−5 10−5 10−5

0.0312 0.0345 0.0321 0.0359 0.0318 0.0342

−0.0388 −0.0402 −0.0389 −0.0405 −0.0381 −0.0393

0.0700 0.0747 0.0710 0.0764 0.0699 0.0735

standard deviation

mean

max

min

extreme deviation

0.0103 0.0103 0.0103 0.0103 0.0102 0.0102

−6.1781 × 10−8 −6.2410 × 10−8 −3.6512 × 10−7 −3.6583 × 10−7 6.6635 × 10−7 6.6592 × 10−7

0.0297 0.0297 0.0297 0.0297 0.0299 0.0299

−0.0338 −0.0338 −0.0340 −0.0340 −0.0333 −0.0333

0.0635 0.0635 0.0637 0.0637 0.0632 0.0632

0.0113 0.0121 0.0114 0.0122 0.0111 0.0119

mean −2.1845 −3.4671 −2.2649 −3.5395 −2.0603 −3.1653

× × × × × ×

Table 6. Statistical Results of y2 for Cases under White Noise item case 1 case 2 case 3

proposed disturbance model proposed disturbance model proposed disturbance model

(8) Mayne, D. Q. Model predictive control: Recent developments and future promise. Automatica 2014, 50 (12), 2967−2986. (9) Maestre, J. M.; Munoz de la Pena, D.; Camacho, E. F. Distributed model predictive control based on a cooperative game. Optimal Control Appl. Methods 2011, 32 (2), 153−176. (10) Christofides, P. D.; Scattolini, R.; Munoz de la Pena, D.; Liu, J. F. Distributed model predictive control: a tutorial review and future research directions. Comput. Chem. Eng. 2013, 51, 21−41. (11) Liu, J. F.; Munoz de la Pena, D.; Christofides, P. D. Distributed model predictive control of nonlinear process systems. AIChE J. 2009, 55 (5), 1171−1184. (12) Farina, M.; Scattolini, R. Distributed predictive control: a noncooperative algorithm with neighbor-to-neighbor communication for linear systems. Automatica 2012, 48 (6), 1088−1096. (13) Conte, C.; Jones, C. N.; Morari, M.; Zeilinger, M. N. Distributed synthesis and stability of cooperative distributed model predictive control for linear systems. Automatica 2016, 69, 117−125. (14) Xu, F.; Chen, H.; Gong, X.; Mei, Q. Fast nonlinear model predictive control on FPGA using particle swarm optimization. IEEE Trans. Ind. Electronics 2016, 63 (1), 310−321. (15) Wang, Y.; Boyd, S. Fast model predictive control using online optimization. IEEE Trans. Control Syst. Technol. 2010, 18 (2), 267− 278. (16) Ohtsuka, T. A continuous/GMRES method for fast computation of nonlinear receding horizon control. Automatica 2004, 40 (4), 563−574. (17) Barros, J. D.; Silva, J. F. A.; Jesus, E. G. A. Fast-predictive optimal control of NPC multilevel converters. IEEE Trans. Ind. Electronics 2013, 60 (2), 619−627. (18) Pannocchia, G.; Rawlings, J. B.; Wright, S. J. Fast, large-scale model predictive control by partial enumeration. Automatica 2007, 43 (5), 852−860. (19) Summers, S.; Jones, C. N.; Lygeros, J.; Morari, M. A multiresolution approximation method for fast explicit model predictive control. IEEE Trans. Autom. Control 2011, 56 (11), 2530−2541. (20) Ulbig, A.; Olaru, S.; Dumur, D.; Boucher, P. Explicit nonlinear predictive control for a magnetic levitation system. Asian J. Control 2010, 12 (3), 434−442. (21) Mariethoz, S.; Morari, M. Explicit model-predictive control of a PWM inverter with an LCL filter. IEEE Trans. Ind. Electronics 2009, 56 (2), 389−399. (22) Ferreau, H. J.; Bock, H. G.; Diehl, M. An online active set strategy to overcome the limitations of explicit MPC. Int. J. Robust Nonlinear Control 2008, 18 (8), 816−830. (23) Munoz de la Pena, D.; Ramirez, D. R.; Camacho, E. F.; Alamo, T. Explicit solution of min-max MPC with additive uncertainties and quadratic criterion. Syst. Control Lett. 2006, 55 (4), 266−274.

situations are similar to that of Figure 32. The overshoot of the first output of the proposed approach is smaller, and the responses of the third output of both strategies are nearly the same. Figures 35−37 show the responses of both methods for all cases under white noises, and the relevant statistical results of the controlled outputs are listed in Tables 5 and 6. These statistical results verify the aforementioned viewpoints further.

9. CONCLUSION In this paper, a summary and new insight of our recent work about the ENMSS model predictive control is presented in detail. On the basis of the ENMSS model, there are more degrees of freedom for the controller design, and the corresponding controllers provide improved control performance. Issues of SISO MPC controller designs as an example, relationship between ENMSSMPC and traditional state space MPC, tuning of the PID controller control, parameter optimization, and industrial applications are discussed.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Ridong Zhang: 0000-0003-4184-3955 Furong Gao: 0000-0002-5900-1353 Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Skogestad, S. Simple analytic rules for model reduction and PID controller tuning. J. Process Control 2003, 13 (4), 291−309. (2) Ziegler, J.; Nichols, N. Optimum settings for automatic controllers. Trans. ASME 1942, 1942 (64), 759−768. (3) Cohen, G. H.; Coon, G. A. Theoretical investigation of retarded control. Trans. ASME 1953, 75, 827−834. (4) Tyreus, B. D.; Luyben, W. L. Tuning PI controllers for integrator/dead time processes. Ind. Eng. Chem. Res. 1992, 31 (11), 2625−2628. (5) Richalet, J. Industrial applications of model-based predictive control. Automatica 1993, 29 (5), 1251−1274. (6) Garcia, C. E.; Prett, D. M.; Morari, M. Model predictive control theory and practice - a survey. Automatica 1989, 25 (3), 335−348. (7) Qin, S. J.; Badgwell, T. A. A survey of industrial model predictive control technology. Control Eng. Practice 2003, 11 (7), 733−764. 5393


Article

Industrial & Engineering Chemistry Research (24) Badgwell, T. A. Robust model predictive control of stable linear systems. Int. J. Control 1997, 68 (4), 797−818. (25) Casavola, A.; Famularo, D.; Franze, G. A feedback min-max MPC algorithm for LPV systems subject to bounded rates of change of parameters. IEEE Trans. Autom. Control 2002, 47 (7), 1147−1153. (26) Mhaskar, P.; El-Farra, N. H.; Christofides, P. D. Robust hybrid predictive control of nonlinear systems. Automatica 2005, 41 (2), 209−217. (27) Casavola, A.; Giannelli, M.; Mosca, E. Min-max predictive control strategies for input-saturated polytopic uncertain systems. Automatica 2000, 36 (1), 125−133. (28) Mhaskar, P. Robust model predictive control design for faulttolerant control of process systems. Ind. Eng. Chem. Res. 2006, 45 (25), 8565−8574. (29) Bemporad, A.; Morari, M. Robust model predictive control: a survey. Robustness in Identification and Control 1999, 245, 207−226. (30) Raffo, G. V.; Ortega, M. G.; Rubio, F. R. An integral predictive/ nonlinear H(infinity) control structure for a quadrotor helicopter. Automatica 2010, 46 (1), 29−39. (31) Yan, Z.; Le, X. Y.; Wang, J. Tube-based robust model predictive control of nonlineat systems via collective neurodynamic optimization. IEEE Trans. Ind. Electronics 2016, 63 (7), 4377−4386. (32) Hariprasad, K.; Bhartiya, S. A computationally efficient robust tube based MPC for linear switched systems. Nonlinear Analysis-Hybrid Systems 2016, 19, 60−76. (33) Bayer, F. A.; Muller, M. A.; Allgower, F. Tube-based robust economic model predictive control. J. Process Control 2014, 24 (8), 1237−1246. (34) Farina, M.; Giulioni, L.; Scattolini, R. Stochastic linear model predictive control with chance constraints - A review. J. Process Control 2016, 44, 53−67. (35) Ocampo-Martinez, C.; Wang, Y.; Puig, V. Stochastic model predictive control based on Gaussian processes applied to drinking water networks. IET Control Theory Appl. 2016, 10 (8), 947−955. (36) Khanmirza, E.; Esmaeilzadeh, A.; Markazi, A. H. D. Predictive control of a building hybrid heating system for energy cost reduction. Appl. Soft Comput. 2016, 46, 407−423. (37) Yu, K. J.; Yang, H. Z.; Tan, X. G.; Kawabe, T.; Guo, Y. A.; Liang, Q.; Fu, Z. Y.; Zheng, Z. Model predictive control for hybrid electric vehicle platooning using slope information. IEEE Trans. Intelligent Transportation Syst. 2016, 17 (7), 1894−1909. (38) Zhang, K.; Sprinkle, J.; Sanfelice, R. G. Computationally aware switching criteria for hybrid model predictive control of cyber-physical systems. IEEE Trans. Automation Sci. Eng. 2016, 13 (2), 479−490. (39) Zhao, J. F.; Wang, J. M. Integrated model predictive control of hybrid electric vehicles coupled with aftertreatment systems. IEEE Trans. Veh. Technol. 2016, 65 (3), 1199−1211. (40) Xin, J. B.; Negenborn, R. R.; Lodewijks, G. Energy-efficient container handling using hybrid model predictive control. Int. J. Control 2015, 88 (11), 2327−2346. (41) Patrinos, P.; Bemporad, A. An accelerate dual gradientprojection algorithm for embedded linear model predictive control. IEEE Trans. Autom. Control 2014, 59 (1), 18−33. (42) Necoara, I. Computational complexity certification for dual gradient method: Application to embedded MPC. Syst. Control Lett. 2015, 81, 49−56. (43) Vouzis, P. D.; Bleris, L. G.; Arnold, M. G.; Kothare, M. V. A system-on-a-chip implementation for embedded real-time model predictive control. IEEE Trans. Control Syst. Technol. 2009, 17 (5), 1006−1017. (44) Besselmann, T. J.; Van de Moortel, S.; Almer, S.; Jorg, P.; Ferreau, H. J. Model predictive control in multi-megawatt range. IEEE Trans. Ind. Electronics 2016, 63 (7), 4641−4648. (45) Takacs, G.; Batista, G.; Gulan, M.; Rohal’-Ilkiv, B. Embedded explicit model predictive vibration control. Mechatronics 2016, 36, 54− 62. (46) Diehl, M.; Amrit, R.; Rawlings, J. B. A lyapunov function for economic optimizing model predictive control. IEEE Trans. Autom. Control 2011, 56 (3), 703−707.

(47) Heidarinejad, M.; Liu, J. F.; Christofides, P. D. Economic model predictive control of nonlinear process systems using lyapunov techniques. AIChE J. 2012, 58 (3), 855−870. (48) Angeli, D.; Amrit, R.; Rawlings, J. B. On average performance and stability of economic model predictive control. IEEE Trans. Autom. Control 2012, 57 (7), 1615−1626. (49) Amrit, R.; Rawlings, J. B.; Angeli, D. Economic optimization using model predictive control with a terminal cost. Ann. Rev. Control 2011, 35 (2), 178−186. (50) Durand, H.; Ellis, M.; Christofides, P. D. Economic model predictive control designs for input rate-of-change constraint handling and guaranteed economic performance. Comput. Chem. Eng. 2016, 92, 18−36. (51) Grune, L. Economic receding horizon control without terminal constraints. Automatica 2013, 49 (3), 725−734. (52) Zhang, R. D.; Lu, J. Y.; Qu, H. Y.; Gao, F. R. State space model predictive fault-tolerant control for batch processes with partial actuator failure. J. Process Control 2014, 24 (5), 613−620. (53) Zhang, R. D.; Lu, R. Q.; Xue, A. K.; Gao, F. Predictive functional control for linear systems under partial actuator faults and application on an injection molding batch process. Ind. Eng. Chem. Res. 2014, 53 (2), 723−731. (54) Zhang, R. D.; Xue, A. K.; Wang, S. Q.; Ren, Z. Y. An improved model predictive control approach based on extended non-minimal state space formulation. J. Process Control 2011, 21 (8), 1183−1192. (55) Zhang, R. D.; Gao, F. R. Multivariable decoupling predictive functional control with non-zero-pole cancellation and state weighting: Application on chamber pressure in a coke furnace. Chem. Eng. Sci. 2013, 94, 30−43. (56) Zhang, R. D.; Xue, A. K.; Wang, S. Q.; Zhang, J. M. An improved state space model structure and a corresponding predictive functional control design with improved control performance. Int. J. Control 2012, 85 (8), 1146−1161. (57) Zhang, R. D.; Wu, S.; Lu, R. Q.; Gao, F. R. Predictive control optimization based PID control for temperature in an industrial surfactant reactor. Chemom. Intell. Lab. Syst. 2014, 135 (15), 48−62. (58) Zhang, R. D.; Cao, Z. X.; Bo, C. M.; Li, P.; Gao, F. R. New PID controller design using extended non-minimal state space model based predictive functional control structure. Ind. Eng. Chem. Res. 2014, 53 (8), 3283−3292. (59) Wu, S. State space predictive functional control optimization based new PID design for multivariable processes. Chemom. Intell. Lab. Syst. 2015, 143 (15), 16−27. (60) Wu, S. Multivariable PID control using improved state space model predictive control optimization. Ind. Eng. Chem. Res. 2015, 54 (20), 5505−5513. (61) Zhang, R. D.; Zou, H. B.; Xue, A. K.; Gao, F. R. GA based predictive functional control for batch processes under actuator faults. Chemom. Intell. Lab. Syst. 2014, 137 (15), 67−73. (62) Zhang, R. D.; Xue, A. K.; Gao, F. R. Temperature control of industrial coke furnace using novel state space model predictive control. IEEE Tran. Ind. Informatics 2014, 10 (4), 2084−2092. (63) Zhang, R. D.; Xue, A. K.; Lu, R. Q.; Li, P.; Gao, F. R. Real-time implementation of improved state-space MPC for air supply in a coke furnace. IEEE Trans. Ind. Electronics 2014, 61 (7), 3532−3539. (64) Wang, L. P.; Young, P. C. An improved structure for model predictive control using non-minimal state space realisation. J. Process Control 2006, 16 (4), 355−371. (65) Rivera, D.; Morari, M.; Skogestad, S. Internal model control 4. PID controller design. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 252−265. (66) Pannocchia, G.; Rawlings, J. B. Disturbance models for offsetfree model-predictive control. AIChE J. 2003, 49 (2), 426−437. (67) Zhang, R. D.; Li, P.; Xue, A. K.; Jiang, A. P.; Wang, S. Q. A simplified linear iterative predictive functional control approach for chamber pressure of industrial coke furnace. J. Process Control 2010, 20 (4), 464−471. (68) Bequette, B. W. Process Control: Modeling, Design and Simulation; Prentice Hall: Upper Saddle River, NJ, 2003. 5394


State space model predictive control for advanced process operation

Recommend Documents