Multilinear Model Decomposition and Predictive Control of MIMO Two

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX-XXX

pubs.acs.org/IECR

Multilinear Model Decomposition and Predictive Control of MIMO Two-Block Cascade Systems Jian Zhang,*,† Kwai Sang Chin,† and Maciej Ławryńczuk‡ †

Department of System Engineering and Engineering Management, City University of Hong Kong, 83 Tatchee Avenue, Kowloon, Hong Kong ‡ Institute of Control and Computation Engineering, Faculty of Electronics and Information Technology, Warsaw University of Technology, ul. Nowowiejska 15/19, 00-665 Warsaw, Poland ABSTRACT: This paper presents a new multilinear model decomposition method for multiple-input multiple-output (MIMO) two-block cascade systems and a model predictive control (MPC) algorithm for the resulting representation. First, a normal vector included angle division method is developed to decompose the operating space and determine the minimum linear model bank through evaluating the nonlinearity of the steady-state I/O surfaces. For a prescribed angle threshold, the minimum linear model bank can be constructed to approximate the original two-block cascade system sufficiently closely. Next, a multilinear MPC algorithm is designed with the proposed trajectory scheduling technique, which can reduce output oscillations caused by hard switching and avoid the difficulty of calculating/tuning complex weighting functions/parameters used in soft switching. A benchmark chemical reactor process is studied to illustrate the effectiveness and advantages of the proposed decomposition method and the predictive control algorithm.

1. INTRODUCTION The multilinear model control approach is an effective control technique for nonlinear processes with wide operating ranges and large set-point changes. Its key concept is to first approximate a nonlinear process as a combination of multiple linear models and then design the system controller based on classical control techniques. Various multilinear model control algorithms have been proposed and applied in practical systems.1,2 The most crucial problem of multilinear model control is how to decompose a nonlinear system, namely, how many and which linear models are required to span the expected operating space of a nonlinear system.2−4 This problem has not been completely solved and is still under study. Determination of the number and the positions of linear submodels is usually performed experimentally, but the trial and error approach is really timeconsuming and costly. In ref 5, Galán et al. first proposed using the gap metric concept to compare the candidate models and reduce the number of linear submodels in a given model bank. Tan et al.6 extended the method in Galán et al.,5 where operating point selection is integrated with local controller design via the loop-shaping H∞ approach. However, the operating points were selected from an existing point set, and the selection requires advance knowledge of the achievable closed-loop performance. Based on gap metric, Du et al.4 first designed the operating range division method for a nonlinear system with one scheduling variable, which does not need the assumption of existing model banks or operating point sets. Then, Du et al.7 generalized the method in ref 4 to accommodate nonlinear systems with multiple © XXXX American Chemical Society

scheduling variables, together with a gridding algorithm for reducing gridding points and decreasing computational burden. To accurately measure the distance between candidate linear models in all ranges of frequencies, Shaghaghi et al.8 studied the H-gap metric based division method. All the aforementioned methods divide the nonlinear system by comparing local linear models. Such methods depend on the system dynamic information around operating points and are just effective in the case where a set of local linear models is available before dividing, such as the cases in refs 5 and 6. Once local linear models are not available, system nonlinear models or a large number of local linear models would have to be identified from measured data, such as the cases in refs 4, 7, and 8. Because identifying an accurate nonlinear model or a large number of local linear models is generally time-consuming and costly, the gap metric based division methods are strictly restricted by the requirement of system dynamic information. Nevertheless, for a general nonlinear system, there seems to be no better alternative than the gap metric based methods. In contrast, for the Hammerstein system (a nonlinear steadystate block followed by a linear dynamic one) and the Wiener system (a linear dynamic block followed by a nonlinear steadystate one), their static I/O mappings (curves or surfaces) can sufficiently represent the nonlinearity and can be used for Received: Revised: Accepted: Published: A

August 7, 2017 November 7, 2017 November 8, 2017 November 8, 2017 DOI: 10.1021/acs.iecr.7b03265 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research multilinear model decomposition. As the two systems both consist of a nonlinear static block and a linear dynamic subsystem, they are called a ”two-block cascade system” in this article. In ref 9, Du et al. propose an included angle division method for single-input-single-output (SISO) two-block cascade systems. The difference between two slope angles at two operating points in the static I/O curve is employed to represent the change of system linearity. Thus, this method just requires the static I/O curve instead of entire system dynamic information. This method is simple and intuitive, but it is not applicable for MIMO systems because, at each operating point of a MIMO system, there are an infinite number of slopes in different directions on its I/O surface. On the other hand, MIMO two-block cascade systems with more than one scheduling variable exist widely in industry, such as the MIMO continuous stirred tank reactor,10,11 the polymerization reactor,10 the binary distillation column,12 the laser-aided powder deposition process,13 and so on. Therefore, a low-cost but effective multilinear model division method is urgently needed for MIMO two-block cascade systems. Once the multilinear model is obtained, the local linear controllers can be initially designed by using many types of control methods (e.g., MPC, PID, and LQR). MPC is recommended here because of its natural ability to deal with MIMO systems and constraints.14 The main problem in controller design is how to integrate these local MPC controllers as a global one. The simplest method is the hard switching method,2 in which only one local controller is selected to calculate the control law and local controllers are switched between different operating subregions. Hard switching of controllers usually leads to unacceptable output oscillations, e.g., the unexpected impulsive phenomena. To reduce output oscillations, many soft switching methods have been proposed based on different weighting functions, such as the Gaussian function based method,15 the Bayesian weighting function based method,16 the gap metric based method,17,18 and so on. However, designing a proper weighting function is generally not easy, because usually the weighting parameters can only be determined by engineering experience and calculation of gap metric requires complete system dynamic information. For these reasons, a simple but smooth switching method is expected for local MPC controller combination. This study focuses on the multilinear model decomposition and control of MIMO two-block cascade systems. The main contributions of this paper are summarized as follows: 1. The proposed normal vector included angle division method. This method extends the idea in ref 9 to MIMO systems by using the normal vector included angle concept. It initially partitions the system I/O surface into a certain number of polytopes and then employs the included angles between the normal vectors to these polytopes to evaluate system nonlinearity and decompose the system. It is usually more practical than the gap metric based methods,4−8 because only static I/O mappings instead of system dynamic information are required for the system decomposition. 2. The designed multilinear MPC algorithm with trajectory scheduling (MMPC-TS). At each sampling instant, by using a scheduling variable trajectory, the multilinear predictive model is transformed to be a certain local linear model within the prediction horizon. The proposed trajectory scheduling technique enables MMPC-TS to

reduce output oscillation in comparison with the control algorithms with hard switching, and meanwhile it avoids the complex weighting functions used in soft switching. A benchmark MIMO continuous stirred tank reactor system is studied to illustrate the effectiveness of the proposed division method and the developed control algorithm. The rest of this paper is organized as follows. Section 2 introduces the two-block cascade system and its multilinear model. Section 3 discusses the normal vector included angle division method. Section 4 details development of the nonlinear MPC based on the multilinear model. Section 5 gives the simulation results for the MIMO CSTR benchmark system. Section 6 concludes the paper.

2. MIMO TWO-BLOCK CASCADE SYSTEM AND ITS MULTILINEAR MODEL In this article, the two-block cascade system is a general term for the Hammerstein system ⎧ vk = f (uk) ⎪ n ⎨ yk = −∑ Fy ⎪ i̅ k−i + ⎪ ⎩ i=1

n

∑ H̅ ivk− i (1)

i=1

and the Wiener system n ⎧ ⎪ i̅ k−i + ⎪ wk = −∑ Fw ⎨ i=1 ⎪ ⎪ y = g(wk) ⎩ k

n

∑ H̅ iuk− i i=1

(2)

nu

where uk ∈  is the system input vector (the vector of manipulated variables), vk ∈ nv and wk ∈ nw are intermediate variable vectors, yk ∈ ny is the system output vector (the vector of controlled variables), k is the sample time, F̅i and H̅ i are the l i n e a r c o e ffi c i e n t m a t r i c e s , f ( ·): nu → nv a n d g( ·): nw → ny are static nonlinear functions. Note the static nonlinear functions are not necessarily invertible for this study. The only assumption on the nonlinear parts is the functions f(·) and g(·) are Lipschitz-continuous.19 This assumption is not restrictive because the static description of numerous technological processes10−13 satisfies the Lipschitzcontinuous condition. The multilinear model of a two-block system is a collection of multiple local linear submodels, with designated operating subregions and operating points scheduled by the scheduling variables. Typically, the scheduling variables can be system inputs, outputs, states, and disturbances. Let s denote the vector of scheduling variables. The variation range of s is the operating space of the system, denoted as Φ. Suppose Φ is divided into Nm subregions Φi (i = 1, ..., Nm). Φi should satisfy Φi ⊆ Φ, Φ1 ∪ ... ∪ ΦNm = Φ and Φi ∩ Φj = Ø for i ≠ j.7 For each subregion Φi (i = 1, ..., Nm), an operating point is assigned, denoted as OPi (uio,yio). To facilitate the following operating space division, the system inputs are selected as the scheduling variables, i.e., s = u, and thereby the operating space Φ becomes the input space. In the following, dividing the operating space is equivalent to dividing the input space. In the ith subregion Φi, a local linear submodel can be acquired around OPi (uio,yio) to approximate the studied two-block cascade system, which is expressed by the following difference equation: B

DOI: 10.1021/acs.iecr.7b03265 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research δyk + F1iδyk − 1 + ... + Fniδyk − n= H1iδuk − 1 + H2i δuk − 2 + ... + Hniδuk − n + δεk , s ∈ Φi

Thus, we can assess the nonlinearity of the I/O surface by indirectly analyzing the nonlinearity between these polytopes, and we finally divide the input space (i.e., the operating space). 2. The normal vector included angle of two hyperplanes can effectively assess the nonlinearity between the two hyperplanes. Suppose the system I/O surface has been well approximated by some polytopes, each of which is contained within a hyperplane of dimensionality nu and corresponds to a unique input region (a simplex). If the normal vector included angle of two hyperplanes is sufficiently small (e.g., θ12 in Figure 2), we can believe

(3)

where δuk = uk − uio, δyk = yk − yio, Fi1, ..., Fin, and Hi1, ..., Hin are model coefficients. The error term δεk = εk − εio, where δεk denotes the general modeling error including model mismatch, external disturbances, etc., and εio is the steady-state modeling error at the operating point OPi. Equation 3 is the multilinear model of a two-block cascade system, which is also called the linear model bank because it contains Nm linear submodels. The two-block cascade system has the important property that its nonlinearities are “primarily” static whereas the dynamic characteristics are “basically” linear. That means only the static system gain changes with operating points whereas the system dynamics almost remain the same in the full operating space. Thus, for a MIMO two-block cascade system, its static I/O surfaces can provide enough information for operating space division and multilinear model construction.

3. NORMAL VECTOR INCLUDED ANGLE DIVISION METHOD In this section, a division method based on system I/O surfaces is developed for MIMO tow-block cascade systems. The steadystate I/O mappings are always assumed to be already available. In practice, they can be obtained by system identification or by first principle modeling. For most systems, such as chemical processes10−13 and large-scale wind tunnels,20 the static I/O mappings are generally more readily available and much easier to identify than detailed dynamic models, especially for the systems with available controllers (even if the controllers may not have perfect performance, e.g., simple PID controllers). Therefore, this method is more practical and lower-cost than dynamic information based division methods.4−8 3.1. Normal vector included angle concept. To divide the operating space, some metric must be found to evaluate the nonlinearity of the system I/O surfaces. Based on the following two properties, the normal vector included angle is proposed to act as the metric: 1. Any Lipschitz-continuous functions def ined on a compact domain can be approximated arbitrarily closely by a PieceWise Linear (PWL) function based on the so-called simplicial partition.21,22 This means a I/O surface characterized by a Lipschitz-continuous nonlinear function can be sufficiently approximated by a finite number of polytopes, as shown in Figure 1. Each polytope corresponds to a simplex region in the input space.

Figure 2. Normal vector included angles.

there is not system nonlinearity and the studied system has a same static gain in the two input regions. Otherwise, if the included angle is very large (e.g., θ13 in Figure 2), we can believe the studied system has significant nonlinearity in the two input regions. Calculation of the normal vector included angle consists of two steps. First, calculate the normal vector to a hyperplane. This can be achieved by using the hyperplane equation or by nu linear independent points on this hyperplane. Because the nu vertexes of a polytope, donated as pt1, pt2, ..., ptnu+1, are naturally linear independent, they are used to compute the corresponding normal vector n = [n1 n2... nnu+1]T, namely ⎧ ⎪ n · (pti − ptj ) = 0, 1 ≤ i < j ≤ nu + 1 ⎨ ⎪ ⎩ nnu + 1 = 1

(4)

In eq 4, pti − ptj (1 ≤ i < j ≤ nu + 1) are actually nu noncollinear vectors on the hyperplane. The last equation “nnu+1 = 1” indicates all the calculated normal vector is assigned to be the positive direction of the yiy (1 ≤ iy ≤ ny) axis. The constant 1 can be replaced by other nonzero constants, which would not influence the final included angle. Once two normal vectors (e.g., ni and nj) are obtained, the included angle between them can be easily computed as ⎛ n ·n ⎞ i j ⎟ θ(ni , nj) = arccos⎜⎜ ⎟ ⎝ |ni nj| ⎠

Figure 1. Piecewise linear function approximation. C

(5) DOI: 10.1021/acs.iecr.7b03265 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

cally combine those simplices to form the final operating subregions. To simplify the following subregion combination and facilitate the description of the final operating subregion, this study regards a hypercube subregion as the basic division unit and finally the entire input space is divided into several hypercuboid regions. The issue of combining two hypercube subregions is studied in subsection 3.3.1, based on which subsection 3.3.2 discusses the method to combine hypercuboid subregions. 3.3.1. Hypercube subregions combination. According to the simplicial partition in section 3.2, in the input space each hypercube H h contains n! simplices and each simplex corresponds to a polytope approximating the I/O surface. Thus, for a given I/O surface, we can calculate n! normal vectors for a hypercube input subregion Hh. Consider two hypercubes Hhi and Hhj in the input space (where hi and hj are two index vectors). For the iyth (iy = 1, ..., ny) I/O surface, two sets of normal vectors can be calculated for Hhi and Hhj:

where | · | represents the magnitude of the vector. The range of included angle is [0, π]. If the included angle is very small (e.g., 0), the system has similar static gains and behaves similarly at the two regions; otherwise, the static gains and the system behaviors change significantly from one region to the other. In this study, the normal vector included angle is used to evaluate the difference of system behavior at the two operating regions and further to develop the division method. 3.2. Subprocedure1: simplicial partition of the input space. Figure 3 shows the division scheme. Prior to introducing

i

i ,1

i ,2

i , nu !

Nhyi = {n hyi , n hyi ..., n hyi

i

i ,1

i ,2

i ,n !

}, Nhyj = {n hyj , n hyi ..., n hyj u } (7)

These two sets of normal vectors contain the system static information, which is used to evaluate system nonlinearity between the two hypercube regions. Through pairwise comparison on the two sets of normal vectors, an included angle vector with nu! × nu! elements is formed as the follows:

Figure 3. Normal vector included angle based division of a operating space in 2 .

i

θhyi , hj =

the main procedure of the division method, two subprocedures are initially proposed according to the theoretical basis in section 3.1, based on which main procedure is developed in section 3.4. The first subprocedure of the division method is to partition the input space. Consider the following input space

i ,1

i ,1

i ,1

i ,2

i , nu !

[θ(n hyi , n hyj ), θ(n hyi , n hyj ), ..., θ(n hyi

i ,n !

, n hyj u )](Tnu !× nu !) × 1 (8)

The largest included angle in θihyi,hj, defined as ΔΘihyi,hj = max{θihyi,hj}, indicates the maximal difference of system behavior in terms of i the iyth system output between Hhi and Hhj. That means, if ΔΘhyi,hj is sufficiently small, the system can be regarded as a linear one in the two subregions and thus the two hypercube regions can be combined together as a hypercuboid region in terms of the iyth system output, and vice versa. Considering the system has ny outputs, the total ny largest included angles can be computed for the ny I/O surfaces, which are written in the following vector form:

D = {u ∈ nu | u̲ iu ≤ uiu ≤ u ̅ iu , iu = 1, 2, ..., nu} (6) iu iu where u and u are the lower bound and the upper bound of uiu, respectively. The simplicial partition of the input space consists of two steps.21,22 First, by specifying the number of divisions (miu) associated with the uiu axis (i.e., partitioning the interval [u iu, uiu] into miu subintervals), the domain D is evenly partitioned into ∏in=u 1miu hypercubes. Each hypercube, denoted as Hh, is indexed by a nu-dimension vector h = [h1 h2, ..., hnu]T, where hiu (1 ≤ iu ≤ nu) indicates the hypercube locates at the iuth subinterval on the uiu axis. Accordingly, the grid step in the axis ui is δiu = (uiu − u iu)/ miu. Second, each hypercube is subdivided into nu! simplices. By using the two-step simplicial partition, the input space D is finally partitioned into ∏in=u 1miu × nu! simplices. More details about the simplicial partition are found in refs 21 and 22. In this subprocedure, the number of divisions are designer-assigned parameters, which are written in the vector m = [m1, m2, ..., mnu]T for notational simplicity. Based on the above simplicial partition, a unique continuous high dimensional surface composed of polytopes can be found to approximate the original I/O surface, as shown in Figure 1. Namely, the system I/O surface is approximated as a set of certain polytopes. 3.3. Subprocedure2: combining subregions. In the subprocedure1, the input space is initially partitioned into hypercube subregions and then is partitioned into different simplex subregions. The subprocedure2 will discuss how to systemati-

i

n

ΔΘhi , hj = [ΔΘ1hi , hj , ΔΘ2hi , hj , ..., ΔΘhyi , hj , ..., ΔΘhyi , hj]nTy × 1 (9)

Given a designer-assigned angle threshold γ = [γ , γ , ···, γ ] where γiy (1 ≤ iy ≤ ny) is the threshold for the iyth steady-state I/O surface, if 1

i

ΔΘhyi , hj < γ iy

for all iy , 1 ≤ iy ≤ ny

2

ny T

(10)

then the two hypercube regions, Hhi and Hhj, can be merged into a large subregion in which the system can be regarded to be linear under the criteria γ. The inequality 10 is the basic condition of combining two hypercube regions, which is also the basis of the following hypercuboid subregion combination. 3.3.2. Hypercuboid subregion combination. As shown in Figure 3, suppose the current divided region, which is also called the “original region” for the current division, starts at hs = [h1s , h2s , n n i i ..., hs u]T and ends at he = [h1e , h2e , ..., he u]T, where 1 ≤ hsu, heu ≤ mnu (1 c ≤ iu ≤ nu). The current “considered region” is in the u direction D


Article

Industrial & Engineering Chemistry Research (c can be 1,2, ..., nu), starts at hc = [h1s , ..., hce + 1, ..., hns u]T, and ends at hc = [h1e , ..., hce + 1, ..., hne u]T. The original region and the considered region are both hypercuboid regions and contain No i i n i i = ∏nuiu=1(heu − hsu + 1) and Nc = ∏iuu=1,iu≠c(heu − hsu + 1) hypercubes, respectively. The following two conditions determine whether the considered region can be combined with the original region: 1. Any two different hypercubes in the considered region satisfy the combination condition (10), namely i

ΔΘhyic , hjc < γ iy


in which the subscripts ic and jc mean

i ΔΘhyic,hjc

3.

(11)

are

computed by using the hypercubes in the considered region. 2. Any hypercube in the original region and any one in the considered region satisfy the combination condition in eq 10, namely i

ΔΘhyio , hjc < γ iy


4.

(12)

in which the subscript io means one hypercube comes from the original region. Writing Θihyic,hjc and Θihyio,hjc in vector form of eq 9, we can finally construct the following two included angle matrices Λcc = [ΔΘhic , hjc ]ny × (Nc(Nc − 1)/2) ,

Λoc = [ΔΘhio , hjc ]ny × (NoNc) (13)

where ΔΘhio,hjc and ΔΘhic,hjc are formed according to the method to compute ΔΘhi,hj in eq 9. Define Λ = [Λoo , Λoc]ny × (Nc(Nc − 1)/2 + NoNc)

5.

(14)

Calculating the maximal values of Λ in a row, we can get the included angle vector 1

2

n

ΔΘ̅ = [ΔΘ̅ , ΔΘ̅ , ..., ΔΘ̅ y]nTy × 1

(15)

Thus, similar to the combination condition in eq 10, the criterion to combine two hypercuboid regions is i

ΔΘ̅ y < γ iy

for all iy ,

1 ≤ i y ≤ ny

(16)

This criterion (eq 16) provides the hypercuboid region combination condition for the main procedure of the division method. 3.4. Main procedure: division of operating space. Without loss of generality, we introduce the division method for a two-input two-output (TITO) two-block cascade system. If there exist more system inputs and outputs, the method can be directly extended to accommodate such situations. The main procedure of the normal vector included angle division method consists of the following steps. 1. Parameter setting. Choose the division number vertor m = [m1, m2]T and the threshold value vector γ = [γ1, γ2]T. Partition the input space into m1 × m2 × 2 simplices by subprocedure1. 2. Initialization. (a) Set hs = [1, 1]T, he = hs. (b) Define the block direction flag vector flagbd = [f lag1bd, f lag2bd]T and set f lagbd = [0, 0]T. The vector

6.

E

f lagbd is used to block the division direction. flagibdu = 1 means the uiu direction is blocked, and vice versa. (c) Define the updating direction flag vector f lagud = [f lag1ud, f lag2ud]T and set f lagud = [0, 0]T. The vector f lagud indicates the division direction in the previous steps. f lagbd and f lagud are used together to determine the following division direction. Determine the division direction. (a) If all directions are blocked (f lagbd = [1, 1]T), go to Step 5. (b) Otherwise, check each direction from u1 to unu in sequence. If f lagibdu = 0 and f lagiudu = 0, then go to Step 4 to conduct division in the direction uiu. Division in the direction uiu. (a) If hieu = miu or some hypercube in the considered region has been classified, set f lagibdu = 1 and f lagiudu = 1. If the vector f lagud = [1, 1]T, set flagud = [0, 0]T. Go to Step 3. (b) Otherwise, calculate the included angle vector ΔΘ̅ for the original region and the considered region. i. If ΔΘ̅ satisfies the hypercuboid combination condition (16), set hieu = hieu + 1. If hieu = miu, set f lagibdu = 1. ii. If ΔΘ̅ does not satisfy the hypercuboid i combination condition (eq 16), set f lagbdu = 1. iii. To update the division direction, set f lagiudu = 1. If the vector flagud = [1, 1]T, set flagud = [0, 0]T. Go to Step 3. Subregion determination and reinitialization. (a) At this step, total (h1e − h1s + 1) × (h2e − h2s + 1) hypercubes have been classified to one group and the corresponding hypercuboid region {u | u 1 + (h1s − 1)δ1 ≤ u1 ≤ u 1 + h1e δ1, u 2 + (h2s − 1)δ2 ≤ u2 ≤ u 2 + h2e δ2} forms an independent subregion. The center of the subregion is chosen as the operating point, around which a linear submodel can be identified. (b) If there are remainder hypercubes to be classified, find the index of unclassified hypercube h = [h1,h2]T with the minimum coordinate h1 in the u1 direction (if h1 = m1, consider the minimum coordinate h2 in the u2 direction, and so forth). Update hs = h = [h1,h2]T, he = hs, and f lagbd = [0, 0]T. Go to Step 3. (c) If all hypercubes have been classified, the subregion division of the two-block cascade system is completed. Go to Step 6. Remove redundant submodels. Suppose Nm subregions are classified by using the above steps. Possibly, not all corresponding submodels are necessary. This means, although the Nm subregions are classified into different parts, the studied system are likely to have similar static gain in some of these subregions. Namely, some submodels can be exempt. Suppose the operating points of two subregions belong to the two hypercubes Hhi and Hhj in the input space, respectively. If the corresponding included angle vector ΔΘhi,hj satisfies the hypercube combination criterion in eq 10, the system static gains in the two subregions are close to one another and only one of the two corresponding submodels is to be contained in the model bank for control. By pairwise comparing the previously obtained subregions, the minimal model bank is finally constructed. DOI: 10.1021/acs.iecr.7b03265 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research ⎧ Δx = A i Δx + B i Δu + B Δε ⎪ k+1 N k N k N ,ε k+1 ⎨ , i ⎪ Δy = CN Δxk ⎩ k

Remark 3.1: This division method is naturally applicable to SISO systems. In such a system, the normal vector included angle has the same effect with the slope included angle. This means that the slope included angle division method proposed by Du et al.9 is a special case of the discussed method. Remark 3.2: The above steps start the division procedure from the hypercube region with the minimum values of system input. The users can change the starting region by initializing hs as other values in step 2(a). Different starting regions may lead to different dividing results. These results are all “optimal” in the sense of the prescribed angle threshold γ. Remark 3.3: The computational burden of this division method is significant. It strongly depends on the number of system inputs nu, because the number of simplices (∏i n=u 1miu × nu!) is highly related to nu. The method discussed is quite computationally efficient when there are 2 or 3 system inputs. This is not a very limiting constraint as numerous technological processes, e.g. chemical reactors10,11,13 and distillation columns,12 satisfy that condition. On the other hand, for systems with a large number of inputs, the division parameters miu (iu = 1,..., nu) can be employed to balance the computational burden and the finer division. Decreasing miu will reduce the computational burden at the cost of rougher divising, and vice versa. Remark 3.4: The above method is efficient for many common industrial processes, but it may become very time-consuming when the I/O surfaces in some outputs are extremely complicated. In such a case, one can first carry out the above procedure for each I/O surface, respectively. In particular, for the I/O surface in which the output changes periodically with the system inputs, one just needs to carry out the division procedure for some one-period operating subregion and then map the dividing results to the remaining operating space accordingly. Next, the resulting subregion boundaries from all the I/O surfaces can be together employed to finally divide the operating space.

(19)

where ANi = ⎡− F i ⎢ 1 ⎢ In × n ⎢ y y ⎢⋮ ⎢ ⎢0 ⎢ ⎢0 ⎢0 ⎢ ⎢⋮ ⎢ ⎣0

BNi =

+

+ ... +

HniΔuk − n

+ Δεk ,

0

0

⋱ ⋮

⋮

⋮

⋯ Iny × ny

0

0

⋯ 0 ⋯ 0

0 0

0 Inu × nu

⋱ ⋮

⋮

⋮

⋯ 0

0

⋯

x

[(H1i)T

x

0 ⋯ 0 Inu × nu 0 ⋯ 0 ]nTx × nu

0 and I are a zero matrix and a unit matrix with appropriate dimensions, respectively. Define the output tracking error as ek = yk − rk (20) where rk is the set-point at the instant k. By combining eqs 19 and 20, we derive ek+1 as follows: ek + 1 = ek + CNi ANi Δxk + CNi BNi Δuk + Δεk + 1 − Δrk + 1 , s ∈ Φi

(21)

By combining the state vector Δxk and the output tracking error ek, the extended nonminimum state vector is defined as24,25

⎡ Δxk ⎤ zk = ⎢ ⎥ ⎣ ek ⎦

(22)

where zk ∈  and nz = ny × (n + 1) + nu × (n − 1). Based on eqs 19, 21, and 22, we get the MENMSS model: nz

zk + 1 = Ai zk + Bi Δuk + C iΔrk + 1 + BεΔεk + 1 ,

s ∈ Φi (23)

s ∈ Φi

where ⎡ Ai ⎡0 ⎤ ⎡ i ⎤ 0nx × ny ⎤ ⎢ N ⎥ i ⎢ BN ⎥ i ⎢ nx × ny ⎥ ;B A =⎢ ⎥ ; B = ⎢ i i ⎥ ; C = ⎢− I i i ⎥ ε ⎣ BN AN ⎦ ⎣ ny × ny ⎦ ⎢⎣CN AN Iny × ny ⎥⎦

(17)

i

Based on the works in Wang23 and Zhang,24 the nonminimum state vector Δxk is chosen as

⎡ BN , ε ⎤ ⎥ =⎢ ⎢⎣ Iny × ny ⎥⎦

Δxk = [ΔykT , ΔykT− 1 , ..., ΔykT− n + 1 , ΔukT− 1 , ΔukT− 2 , ..., ΔukT− n + 1]T

⋯ 0

⋯ Hni − 1 Hni ⎤ ⎥ 0 ⎥ ⋯ 0 ⎥ ⋱ ⋮ ⋮ ⎥ ⎥ 0 ⎥ ⋯ 0 ⎥ 0 ⎥ ⋯ 0 0 ⎥ ⋯ 0 ⎥ ⋱ ⋮ ⋮ ⎥ ⎥ 0 Inu × nu 0 ⎦ n ×n

CNi = [ Iny × ny 0 ⋯ 0 0 0 ⋯ 0 ]ny × nx

Δyk + F1iΔyk − 1 + ... + FniΔyk − n= H2i Δuk − 2

⋯ −Fni − 1 −Fni H2i

BN , ε = [ Iny × ny 0 ⋯ 0 0 0 ⋯ 0 ]nTx × ny

4. MULTILINEAR MODEL PREDICTIVE CONTROL In this part, the multilinear extended nonminimal state space (MENMSS) model is initially formed to reject process disturbance and deal with model mismatch. Then, an improved multilinear model predictive control algorithm is designed with the proposed trajectory scheduling technique. 4.1. Multilinear extended nonminimal state space model. Consider the multilinear model in the form of eq 3. By adding the back shift operator Δ to eq 3, eq 3 is expressed as

H1iΔuk − 1

s ∈ Φi

(18)

(24)

Note, the frequently used ”DMC type” disturbance model is adopted in the following control algorithm. Over the prediction horizon in the following control algorithm, the error term εk and its increment Δεk are assumed to be a constant and zero, 11,16

where Δxk ∈ nx and nx = ny × n + nu × (n − 1). Then, eq 17 can be transformed into the following state space model: F


Article


Suppose a set of local linear controllers are designed based on the minimum model bank. To simplify design of the following global controller, we choose the same predictive horizon P and control horizon M for these local controllers, whereas the weights i i are separately tuned, donated as wz,k + j|k and wΔu,k + j|k (i = 1, ..., Nm), respectively. Then, the proposed MMPC-TS algorithm is developed as follows. Define the predicted scheduling variable trajectory at sampling instant k as

respectively. Therefore, the effect of the error term, the term BεΔεk+1 in eq 23, is exempt from the following algorithm. 4.2. Predictive control based on the MENMSS model. 4.2.1. MPC optimization problem. At each sampling instant, MPC determines the future sequence of input increments ΔUk = [ΔukT| k ΔukT+ 1 | k ... ΔukT+ Hc − 1 | k]T

by solving the following problem with constraints P

∑

min

Δuk | k , ···,

j=1

Δuk + M − 1 | k

s. t.

Sk = [skT+ 1 | kskT+ 2 | k ... skT+ P | k]T

M−1

zk + j | k

2 wz , k + j | k

+

∑

Δuk + j | k

j=0

2 wΔu , k + j | k

and define ⎡ zk + 1 | k ⎤ ⎢z ⎥ k+2|k ⎥ Zk = ⎢ , ⎢ ⋮ ⎥ ⎢ ⎥ ⎣ zk + P | k ⎦

zk + j | k = Ai zk + j − 1 | k + Bi Δuk + j − 1 | k + C iΔrk + j | k , if sk + j | k ∈ Φi(i = 1, ..., Nm) Δuk + j | k = 0,

(26)

if j ≥ M

umin ≤ uk + j | k ≤ umax

⎡ Δrk + 1 | k ⎤ ⎥ ⎢ ⎢ Δrk + 2 | k ⎥ ΔR k = ⎢ ⎥ ⎢ ⋮ ⎥ ⎥ ⎢ Δr ⎣ k+P|k ⎦

(27)

Taking into account the MENMSS model (eq 23), the predicted state trajectory can be derived as

Δumin ≤ Δuk + j | k ≤ Δumax

Zk = Fkzk + Gk ΔUk + ΨkΔR k

(25)

(28)

where Fk, Gk, and Ψk are parameter matrices

where the norm is defined as ∥x∥2W ≜ xTWx; P is the prediction horizon; M is the control horizon; wz,k+j|k ≥ 0nz×nz, wΔu,k+j|k ≥ 0nu×nu are the weights of the minimized cost-function; umin ∈ nu , umax ∈ nu , Δumin ∈ nu , and Δumax ∈ nu are the constraints imposed on the magnitude and on the increments of system inputs, respectively. Here, we ignore output constraints for simplicity. Interested readers can refer to Wang,23 which shows output constraints can be easily integrated into the problem (eq 25). At each sampling instant, M future control increments are calculated from the MPC optimization problem (eq 25); then only the increments for the instant k are actually applied to the process, namely uk = uk−1 + Δuk. At the next sampling instant, the whole optimization procedure is repeated. Because different local linear models are used to characterize the process in different subregions Φi, the predicted states zk+j|k are essentially nonlinear functions of the optimized online control increments ΔUk. Directly solving the complex nonlinear optimization problem (eq 25) for ΔUk is computationally difficult and will lead to a heavy computational burden. The more practical approach is first to design local MPC controllers from the local linear models and then to combine these local controllers into a global one.2,17 Based on this idea, a computational efficient multilinear MPC algorithm with trajectory scheduling (MMPC-TS) is developed in the following subsection. 4.2.2. Multilinear MPC algorithm with trajectory scheduling. In the proposed MMPC-TS algorithm, the multilinear model is scheduled as different linear models over the predictive horizon by using a predicted scheduling variable trajectory. Then, the predicted states zk+j|k become the linear functions of the future control increments ΔUk, and the nonlinear programming problem (eq 25) is simplified into a quadratic programming (QP) problem. So, ΔUk can be easily computed (the computational burden of quadratic optimization is much lower than that of general nonlinear optimization and the global optimal solution is always found). For accurate scheduling, the aforementioned trajectory scheduling and the resulting QP optimization may be repeated several times in internal iterations at each sampling instant.

⎡ Ak + 1 | k ⎤ ⎢ ⎥ ⎢ Ak + 2 | k Ak + 1 | k ⎥ ⎢ ⎥ ⋮ Fk = ⎢ ⎥ ⎢ 1 ⎥ ⎢ ∏ Ak + j | k ⎥ ⎢⎣ j = P ⎥⎦ n P×n z

⎡ Bk + 1 | k ⎢ ⎢ Ak + 2 | k Bk + 1 | k ⎢ ⋮ ⎢ ⎢ 2 ⎢ ∏ A ⎢ j=M+1 k+j|k Gk = ⎢ ⎢ Bk + 1 | k ⎢ ⋮ ⎢ ⎢ 2 ⎢ ∏ Ak + j | k ⎢ j=P ⎢ Bk + 1 | k ⎣

(29)

z

0

⋯

Bk + 2 | k

⋯

⋮

⋱

3

∏

Ak + j | k ⋯

j=M+1

Bk + 2 | k ⋮

⋱

∏ Ak + j | k

⋯

3 j=P

Bk + 2 | k

⎤ ⎥ 0 ⎥ ⎥ ⋮ ⎥ ⎥ Ak + M + 1 | k Bk + M | k ⎥ ⎥ ⎥ ⎥ ⎥ ⋮ ⎥ M+1 ⎥ ∏ Ak + j | k ⎥ ⎥ j=P ⎥ Bk + M | k ⎦n P × n M 0

z

u

(30)

Ψk = ⎡Ck + 1 | k 0 ⎢ Ck + 2 | k ⎢ Ak + 2 | k Ck + 1 | k ⎢ ⋮ ⎢⋮ ⎢ 2 3 ⎢∏ A C A C ∏ ⎢⎣ j = P k + j | k k + 1 | k j = P k + j | k k + 2 | k

⎤ ⎥ ⋯ 0 ⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ Ck + P | k ⎥ ⎥⎦ n P×n P ⋯ 0

z

y

(31)

Based on the scheduling variable sk+j|k(j = 1, ..., P), the model matrices Ak+j|k, Bk+j|k, and Ck+j|k are chosen from the obtained MENMSS model. For example, if sk+j|k ∈ Φi (i = 1, ..., Nm), then Ak+j|k = Ai, Bk+j|k = Bi, and Ck+j|k = Ci. G


Article


3. Calculate the predicted state trajectory Ztk, i.e., the matrices Ftk, Gtk and Ψtk are found, as well as the weighting matrices Wtz,k and WtΔu,k. 4. The MPC QP problem (eq 32) is solved to find ΔUtk and the vector Stk is updated by using ΔUtk. 5. Decide whether internal iterations should be terminated. If the trajectory ΔUtk is close to the trajectory at the previous internal iteration, i.e.

Taking into account the derived state trajectory (eq 28) and considering that Uk = LΔUk + Uk−1 where ⎡ In × n 0 ⎢ u u ⎢ In × n In × n L=⎢ u u u u ⋮ ⎢ ⋮ ⎢I ⎣ nu × nu Inu × nu ⎡ uk − 1⎤ ⎢u ⎥ k − 1⎥ Uk − 1 = ⎢ ⎢ ⋮ ⎥ ⎢ ⎥ ⎣ uk − 1⎦

0 ⎤ ⎥ ⋯ 0 ⎥ , ⎥ ⋱ ⋮ ⎥ ⋯ Inu × nu ⎥⎦ n M×n M ⋯

u

u

ΔUkt − ΔUkt− 1

the nonlinear optimization problem (eq 25) is simplified as the following QP problem ΔUk

2 Wz , k

+ ΔUk

2 WΔu , k

s . t . Ul ≤ Uk − 1 + LΔUk ≤ Uu ΔUl ≤ ΔUk ≤ ΔUu

(32)

where the vectors of length nuM are ⎡ umax ⎤ ⎢u ⎥ max ⎥ , Uu = ⎢ ⎢ ⋮ ⎥ ⎢ ⎥ ⎣ umax ⎦ ⎡ Δumax ⎤ ⎢ ⎥ ⎢ Δumax ⎥ ΔUu = ⎢ ⎥ ⎢ ⋮ ⎥ ⎢⎣ Δu ⎥⎦ max

⎡ umin ⎤ ⎢u ⎥ min ⎥ Ul = ⎢ , ⎢ ⋮ ⎥ ⎢u ⎥ ⎣ min ⎦

< σu

(33)

or t > tmax, then go to Step 6 (σu is a threshold value to be tuned). Otherwise, update t: = t + 1 and go to Step 3. 6. Apply the first nu elements of the vector ΔUtk (i.e., the vector Δutk|k) to the process. At the next sampling instant, the algorithm starts from Step 1 and the above six-step procedure is repeated. Remark 4.1: If tmax is chosen as 1 and sk1 + j|k (1 ≤ j ≤ Nm) are set to the same value, the proposed MMPC-TS reduces to the conventional MPC algorithm with hard switching. This means MPC with hard switching is a special case of this developed MMPC-TS algorithm. Remark 4.2: Compared with the MPC algorithms with soft switching, the proposed MMPC-TS algorithm does not require any tuning parameters for switching, or the complete process dynamic information which is difficult to be obtained in nonlinear systems. Due to the trajectory scheduling technique, the MMPC-TS algorithm still can alleviate output oscillations in comparison with hard switching.

nuM × 1

min Fkzk + Gk ΔUk + ΨkΔR k

2

⎡ Δumin ⎤ ⎢ ⎥ ⎢ Δumin ⎥ ΔUl = ⎢ ⎥, ⎢ ⋮ ⎥ ⎢⎣ Δu ⎥⎦ min

5. CASE STUDY: A MIMO CONTINUOUS STIRRED TANK REACTOR PROCESS 5.1. Process description. Consider a MIMO continuous stirred tank reactor (CSTR) process, which has been modeled by a Hammerstein structure or a Wiener structure in previous studies.10,11 It consists of an irreversible, exothermic reaction, A → B, and takes place in a constant volume reactor cooled by a single coolant stream. It can be characterized by the following continuous-time nonlinear equations10,11 q ⎧ ̇ −E / RT ⎪CA = V (CA 0 − CA) − k 0CAe ⎪ ⎪ q ΔHk 0 −E / RT CAe ⎨Ṫ = (T0 − T ) − ρCp V ⎪ ρc Cpc ⎪ + qc(1 − e−hA / ρc Cpcqc)(Tc 0 − T ) ⎪ ρ C V p ⎩ (34)

the matrix Wz,k = blockdiag {wz,k+1|k, wz,k+2|k, ..., wz,k+P|k} and the matrix WΔu,k = blockdiag{wΔu,k|k,wΔu,k+1|k, ..., wΔu,k+M−1|k}. Similar with the model matrices (Ak+j|k, Bk+j|k, Ck+j|k), wz,k+j|k (1 ≤ j ≤ P) and wΔu,k+j|k (0 ≤ j ≤ M − 1) are selected from the designed local controllers based on the scheduling variable sk+j|k (if sk+j|k ∈ Φi, i i wz,k+j|k = wz,k + j|k, and wΔu,k+j|k = wΔu,k + j|k). If the scheduling variable trajectory Sk is known at sampling instant k, the nonlinear MPC optimization problem (eq 25) can be accurately transformed into the QP problem (eq 32). Unfortunately, the trajectories Sk are generally unknown because they are usually functions of optimized variables. Inaccurate trajectory Sk will deteriorate control performance, which has a similar effect with model mismatch. To predict Sk as accurately as possible, the trajectory scheduling (eq 28) and the resulting QP optimization are repeated several times in internal iterations. For convenience, index the related vectors and matrices at the tth iteration with superscript t, such as the vector Stk, the matrix Ftk, and so on. The proposed MMPC-TS algorithm is summarized as the following steps: 1. Measure the system output yk. Form the ENMSS state space vector zk according to eqs 18, 20, and 22. 2. Set the index t to 1. Initialize the scheduling variable trajectory S1k = [(s1k+1|k)T(s1k+2|k)T ... (sbk1 + P|k)T]T (sk1 + j|k, j = 1, ..., P, are uniformly chosen as uk−1 because system inputs are used as the scheduling variables in this study).

The nominal model parameters are listed in Table 1. The sampling period is Ts = 0.1 min. The process output variables are the concentration of A (CA) and the temperature (T). The process input variables are the process flow rate (q) and the coolant flow rate (qc). The control objective is to regulate CA and T for set-point tracking by manipulating q and qc. The variation ranges of the inputs are q ∈ [95, 150] and qc ∈ [60, 110]. The variation ranges of the outputs are CA ∈ [0.02, 0.15] and T ∈ [430, 490]. Figure 4 and Figure 5 show the steadystate surfaces of the concentration CA and the temperature T, respectively. This CSTR shows apparent nonlinearity, and thus a single linear controller is not sufficient for the wide operating range.10,11 5.2. Division of the operating space. The proposed normal vector included angle method is employed to decompose H


Article


Case 1. m = [10, 10]T, γ = [0.003, π/10]T. By using the first five steps of the proposed division method, the operating space is initially decomposed into three subregions, donated as Φ1, Φ2, and Φ3. The corresponding three operating points are OP1 (122.5, 77.5, 0.0303, 468.79), OP2 (106, 102.5, 0.0855, 443.04) and OP3 (133.5, 102.5, 0.0562, 456.83). Next, according to Step 6 of the division method, the operating points OP1 and OP3 are close in the normal vector included angle sense, where the CSTR has similar static gains. Thus, subregions Φ1 and Φ3 can share a common submodel (which is selected as the submodel at OP1 in this study). The final model bank should only include the two linear submodels at OP1 and OP2. Case 2. m = [30, 30]T, γ = [0.003, π/10]T. In this case, the operating space is also divided into three subregions and the final model bank should also include two submodels at the two operating points OP1(122.5, 76.67, 0.0295, 469.39) and OP2 (104.17, 101.67, 0.0860, 442.55). Case 3. m = [60, 60]T, γ = [0.003, π/10]T. As shown in Figure 6, this case has the same dividing results with Case 2. In the above three cases, the designer-assigned parameters h = [m1, m2]T are selected to be different values, whereas the threshold value vector γ is chosen as the same [0.003, π/10]T. The three cases have almost the same dividing results: the same number of subregions, nearly the same subregion ranges, and nearly the same operating points. In particular, Case 2 has exactly the same dividing results as Case 3. The larger the numbers of divisions m1 and m2 are, the smaller simplices the input region is partitioned into and the more system characteristics the dividing procedure considers. Compared with Case 1, Case 2 and Case 3 obtain the narrower subregion Φ1 in the qc direction and the narrower subregion Φ2 in the q direction. This is because larger m1 and m2 lead to finer and more accurate dividing in Case 2 and Case 3, while smaller m1 and m2 make the division in Case 1 rougher. Therefore, sufficienly large numbers of divisions miu (iu = 1, 2, ..., nu) are generally suggested for accurate dividing (e.g., m = [30, 30]T for the CSTR) if the off-line computation time is allowable. Case 4. m = [30, 30]T, γ = [0.0028, π/30]T. In this case, the operating space is divided into five subregions as depicted in Figure 6. The subregions Φ2 and Φ3 share the same submodel with the subregion Φ1. The minimal model bank resulted from the proposed division method includes the three submodels at the three operating points: OP1 (113.33, 68.33, 0.0247, 471.69), OP3 (104.17, 101.67, 0.0860, 442.55), and OP5 (140.83, 85, 0.0318, 470.77), respectively. Case 4 and Case 2 choose the same numbers of divisions (m = [30, 30]T), whereas Case 4 uses a smaller threshold values γiy (iy = 1, 2, ..., ny) than Case 2. The smaller the threshold value γiy is, the smaller nonlinearity the dividing procedure allows in each subregion for the iyth I/O surface, and vice versa. For this reason, the operating space is decomposed into more subregions in Case 4 than in Case 2. In practice, proper threshold values should be chosen according to the prior knowledge of the system and the control performance. In this study, the dividing result of Case 2 is employed to construct the linear model bank. Two second-order TITO linear submodels in the form of eq 3 are respectively identified at OP1 and OP2, whose parameters are given as follows:

Table 1. Nominal Model Parameters of CSTR Measured product concentration Reactor temperature Coolant flow rate Process flow rate Feed concentration Feed temperature Inlet coolant temperature CSTR Volume Heat transfer term Reaction rate constant Activation energy term Heat of reaction Liquid densities Specific Heats

CA T qc q CA0 T0 Tc0 V hA k0 E/R ΔH ρ, ρc Cp,Cpc

0.1 mol L−1 438.5 K 103.41 L min−1 100 L min−1 1 mol L−1 350 K 350 K 100 L 7 × 105 cal min−1 K−1 7.2 × 1010 min−1 1 × 104 K −2 × 105 cal min−1 1 × 103 g L−1 1 cal g−1 K−1

Figure 4. Steady-state surface of the concentration CA.

Figure 5. Steady-state surface of the concentration T.

the CSTR system. The operating points (OP) for the CSTR are expressed in the form of (q, qc, CA, T). Four cases of division are studied to show the effect of the designer-assigned parameters. The dividing results are depicted in Figure 6. I


Article


Figure 6. Dividing results of the four cases for the CSTR.

⎤ ⎡−0.84371 0 F11 = ⎢ ⎥, ⎣0 −0.84371⎦ ⎤ ⎡ 3.8045 × 10−2 0 ⎥, F21 = ⎢ ⎢⎣ 0 3.8045 × 10−2 ⎥⎦ ⎡ 2.148 × 10−4 1.3291 × 10−4 ⎤ ⎥, H11 = ⎢ ⎢⎣ 2.7358 × 10−2 − 1.3392 × 10−1⎥⎦ ⎡−2.7355 × 10−4 4.7247 × 10−5 ⎤ ⎥, H21 = ⎢ ⎣⎢ 5.277 × 10−2 − 4.6077 × 10−3 ⎦⎥ ⎤ ⎡−1.5235 0 F12 = ⎢ ⎥, ⎣0 −1.5235⎦ ⎡ 6.4138 × 10−1 0 ⎤ ⎥, F22 = ⎢ ⎣⎢ 0 6.4138 × 10−1⎦⎥ ⎡ 6.1264 × 10−4 1.9117 × 10−4 ⎤ ⎥, H12 = ⎢ ⎢⎣−3.4271 × 10−2 − 1.1993 × 10−1⎥⎦ ⎡−8.1253 × 10−4 1.6477 × 10−4 ⎤ ⎥ H22 = ⎢ ⎢⎣1.0622 × 10−1 3.1286 × 10−2 ⎥⎦

algorithm is applied to the CSTR for set-point tracking. For comparison, two simple Linear MPC (LMPC) controllers (LMPC1 and LMPC2) are designed based on the linear submodel 1 and submodel 2, respectively. They are employed to demonstrate the necessity of utilizing the multilinear model technique and the effectiveness of the proposed division method. Moreover, a multilinear model predictive controller with hard switching (MMPC-HS) is also designed based on the same multilinear model. It is also compared with MMPC-TS to illustrate the advantages of the trajectory scheduling technique introduced in section 4.2.2. Initially, the two LMPC controllers are tuned with the 6 6 6 1 parameters P = 20, M = 5, wz,k + j|k = diag{10 , 1, 10 , 1, 0, 0, 10 , 1 2 5 5 1}, wΔu, k + j|k = diag{1, 1}, wz,k + j|k = diag{10 , 1, 10 , 1, 0, 0, 105, 2 1}, wΔu, k + j|k = diag{1, 1}. Constraints on the manipulated variables are 95 ≤ q ≤ 150, 60 ≤ qc ≤ 110, − 10 ≤ Δq ≤ 10, and −10 ≤ Δqc ≤ 10. Next, these parameters are directly used as the control parameters of the MMPC-HS controller and the proposed MMPC-TS controller. Additional parameters of the MPC-TS controller are σu = 1, tmax = 4. Figure 8 depicts the closed-loop performance of the MMPCTS controller. The process outputs CA and T follow the setpoints quickly and precisely, and the process inputs vary exactly in their predetermined ranges, q ∈ [95, 150] and qc ∈ [60, 110]. More importantly, the closed-loop tracking is very smooth, during which neither output chattering nor input chattering happen. The MMPC-TS controller based on the two linear submodels performs quite satisfactorily in the wide operating range. This illustrates two linear submodels and the corresponding divided subregions provide sufficient information for this setpoint tracking control, and the combination of the parameters m = [30, 30]T and γ = [0.003, π/10]T is a good choice for the division of the MIMO CSTR. The comparisons between the MMPC-TS controller and the other two linear MPC controllers (LMPC1, LMPC2) are

Figure 7 compares the modeling performance of the two linear submodels and the multilinear model. The linear submodel 1 cannot track the sampled data with high concentration (CA) and low temperature (T), e.g., the trajectories at around the 1000th sampling point. The submodel 2 even gives negative concentration, which is physically not possible. By comparison, the global multilinear model is very precise, whose trajectories are very close to the sampled data. To sum up, the multilinear model can closely approximate the nonlinear MIMO CSTR process in the entire operating space, with much higher accuracy than a single linear submodel. 5.3. Predictive control of the CSTR. In this part, based on the obtained linear model bank, the proposed MMPC-TS J


Article


Figure 7. Open-loop model validation of the MIMO CSTR: the linear submodels vs the multilinear model.

Figure 8. Simulation results of the MIMO CSTR using the MMPC-TS controller.

the subregions Φ1 and Φ3. This illustrates a single submodel is sufficient to characterize the system dynamic in the two

displayed in Figures 9−12. In Figure 9, the LMPC1 controller performs as well as the MMPC-TS controller at the set-points in K


Article


Figure 9. Closed-loop responses: the MMPC-TS controller vs the LMPC1 controller.

Figure 11. Closed-loop responses: the MMPC-TS controller vs the LMPC2 controller.

Figure 10. Control inputs: the MMPC-TS controller vs the LMPC1 controller.

Figure 12. Control inputs: the MMPC-TS controller vs LMPC2 controller.

subregions. However, the LMPC1 controller leads to serious output oscillations at the set-points in the subregion Φ2. Similarly, in Figure 11, the LMPC2 controller has the same satisfactory performance as the MMPC-TS controller at the setpoints in the subregion Φ2, whereas large overshoots and slow tracking responses happen in other subregions (Φ1 and Φ3). As a whole, the linear MPC controllers based on the two linear submodels just perform well in their own subregions. Therefore, two linear submodels are necessary for this nonlinear MIMO CSTR process. This also confirms the effectiveness of the proposed normal vector included angle division method. Figures 13−14 compare the MMPC-TS controller with the MMPC-HS controller. The MMPC-HS controller can regulate the process outputs to follow their set-points, and it performs as well as the MMPC-TS controller when no submodel switching is activated, such as the responses from k = 0 to k = 50 and from k = 201 to k = 300. Nevertheless, during switching submodels, such as the responses from k = 51 to k = 200 and from k = 301 to k = 350, the MMPC-HS controller causes unexpected large output oscillations. This is due to the rough and inaccurate switching

caused by the hard switching technique in the MMPC-HS controller. By comparison, the proposed trajectory scheduling technique can achieve smooth switching, and its internal iteration mechanism makes the predicted trajectory and the resulting switching in the MMPC-TS controller more accurate than those in the MMPC-HS controller. The total computational time of the MMPC-TS controller is 1.28× that of the MMPC-HS controller. Therefore, the trajectory scheduling technique enables the MMPC-TS controller to give better performance at a slightly higher computational time cost than the MMPC-HS controller. Moreover, it does not require extra complicated weighting parameters/functions used in soft switching. Thanks to the proposed trajectory scheduling technique, the resulting MMPC-TS controller achieves satisfactory balance between smooth control and difficulty in computing control laws and designing controllers. To compare these controllers more systematically, two kinds of criteria are calculated for the above simulations: the settling time Ts26 and the overshoot Mp.27 The settling time Ts refers to L


Article


times), the average settling time (Ts), the average overshoot (Mp), the maximal settling time (Tmax s ), and the maximal overshoot (Mmax p ) of these stages are calculated. Table 2 gives the values of Ts and Tmax (in terms of the three s most common prescribed percentages27) and the values of Mp and Mmax resulted from these MPC controllers. The proposed p MMPC-TS has the minimal values in terms of all these criteria among these controllers. It contributes to faster closed-loop responses and higher control precision than other controllers. Therefore, the effectiveness of the proposed division method and the designed MMPC-TS control algorithm is confirmed.

6. CONCLUSION This study provides a complete multilinear control scheme for MIMO two-block cascade systems. First, a normal vector included angle division method is proposed to decompose the system operating space and determine the minimum linear model bank. It works for the systems with multiple inputs and is an extension of the method proposed by De et al.9 This method just requires steady-state I/O mappings and thus is more convenient and lower-cost than dynamic information based division methods. Moreover, a novel multilinear MPC algorithm (MMPC-TS) is developed with the trajectory scheduling technique for local controller combination. The trajectory scheduling technique enables MMPC-TS to reduce output oscillations at a slightly higher computational cost compared with the control algorithms with hard switching, and meanwhile it avoids the complex weighting parameters/functions used in soft switching. The presented division method has a great potential when applied to the classical technological processes with 2 or 3 inputs. For systems with a large number of inputs, users can use the parameters miu to balance the computational load and the division precision. Further research should focus on reducing computational load without losing precision.

Figure 13. Closed-loop responses: the MMPC-TS controller vs the MMPC-HS controller.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jian Zhang: 0000-0002-1530-6323

Figure 14. Control inputs: the MMPC-TS controller vs the MMPC-HS controller.

Notes

The authors declare no competing financial interest.

■

the time from changing set-points to the instant all the controlled variables first reach and thereafter remain within a prescribed percentage ± ϵ of their corresponding set-points, which reflects the response speed of the closed-loop system. The overshoot Mp is the percentage between the maximum amount the system overshoots its final value and the final value of the system output, which is an important stability indicator for control. As the simulation consists of seven stages (the set-point changes seven

REFERENCES

(1) Porfírio, C.; Neto, E. A.; Odloak, D. Multi-model predictive control of an industrial C3/C4 splitter. Control Eng. Pract. 2003, 11, 765−779. (2) Murray-Smith, R.; Johansen, T. Multiple model approaches to nonlinear modelling and control; Taylor & Francis: London, England, 1997. (3) Galan, O.; Romagnoli, J. A.; Palazoglu, A. Real-time implementation of multi-linear model-based control strategies−an application to a

Table 2. Control Performance of the MPC Controllersa Tmax s /min

Ts/min

a

Mmax p

Mp

Algorithm

ϵ = 5%

ϵ = 2%

ϵ = 1%

ϵ = 5%

ϵ = 2%

ϵ = 1%

CA

T

CA

T

MMPC-TS LMPC1 LMPC2 MMPC-HS

1.214

1.614

1.929

1.6

2.2

2.7

2.59%

1.07%

6.39%

5.36%

1.357 1.386

1.914 1.800

2.357 2.086

2.0 2.1

2.3 2.8

2.8 2.9

4.82% 4.41%

3.08% 2.63%

12.99% 10.00%

6.08% 9.22%

Because LMPC1 leads to serious output oscillations in some set-points, these criteria are not applicable to LMPC1. M


Article

Industrial & Engineering Chemistry Research bench-scale pH neutralization reactor. J. Process Control 2004, 14, 571− 579. (4) Du, J.; Song, C.; Li, P. Application of gap metric to model bank determination in multilinear model approach. J. Process Control 2009, 19, 231−240. (5) Galán, O.; Romagnoli, J. A.; PalazoÇ ğlu, A.; Arkun, Y. Gap metric concept and implications for multilinear model-based controller design. Ind. Eng. Chem. Res. 2003, 42, 2189−2197. (6) Tan, W.; Marquez, H. J.; Chen, T.; Liu, J. Multimodel analysis and controller design for nonlinear processes. Comput. Chem. Eng. 2004, 28, 2667−2675. (7) Du, J.; Song, C.; Yao, Y.; Li, P. Multilinear model decomposition of MIMO nonlinear systems and its implication for multilinear modelbased control. J. Process Control 2013, 23, 271−281. (8) Shaghaghi, D.; Fatehi, A.; Khaki-Sedigh, A. Multi-linear model set design based on the nonlinearity measure and H-gap metric. ISA Trans. 2017, 68, 1−13. (9) Du, J.; Song, C.; Li, P. Multilinear model control of Hammersteinlike systems based on an included angle dividing method and the MLDMPC strategy. Ind. Eng. Chem. Res. 2009, 48, 3934−3943. (10) Cervantes, A. L.; Agamennoni, O. E.; Figueroa, J. L. A nonlinear model predictive control system based on Wiener piecewise linear models. J. Process Control 2003, 13, 655−666. (11) Ławryńczuk, M. Suboptimal nonlinear predictive control based on multivariable neural Hammerstein models. Appl. Intell. 2010, 32, 173−192. (12) Fruzzetti, K.; Palazoğlu, A.; McDonald, K. Nolinear model predictive control using Hammerstein models. J. Process Control 1997, 7, 31−41. (13) Cao, X.; Ayalew, B. Control-oriented mimo modeling of laser-aided powder deposition processes. American Control Conference. 2015; pp 3637−3642. (14) Darby, M. L.; Nikolaou, M. MPC: Current practice and challenges. Control Eng. Pract. 2012, 20, 328−342. (15) Foss, B. A.; Johansen, T. A.; Sørensen, A. V. Nonlinear predictive control using local models - applied to a batch fermentation process. Control Eng. Pract. 1995, 3, 389−396. (16) Aufderheide, B.; Bequette, B. W. Extension of dynamic matrix control to multiple models. Comput. Chem. Eng. 2003, 27, 1079−1096. (17) Du, J.; Johansen, T. A. A gap metric based weighting method for multimodel predictive control of MIMO nonlinear systems. J. Process Control 2014, 24, 1346−1357. (18) Du, J.; Johansen, T. A. Integrated multimodel control of nonlinear systems based on gap metric and stability margin. Ind. Eng. Chem. Res. 2014, 53, 10206−10215. (19) Polak, E. Optimization: algorithms and consistent approximations; Springer-Verlag: New York, 2012; Vol. 124. (20) Wang, X.; Yuan, P.; Mao, Z. The modified feature subsets ensemble applied for the mach number prediction in wind tunnel. IEEE Trans. Aerosp. Electron. Syst. 2016, 52, 863−874. (21) Julian, P.; Desages, A.; Agamennoni, O. High-level canonical piecewise linear representation using a simplicial partition. IEEE Trans. Circuits Syst. I, Fundam. Theory 1999, 46, 463−480. (22) Julian, P. A high level canonical piecewise linear representation: theory and applications. Ph.D. thesis, Universidad Nacional del Sur, Bahia Blanca, Argentina, 1999. (23) Wang, L.; Young, P. C. An improved structure for model predictive control using non-minimal state space realisation. J. Process Control 2006, 16, 355−371. (24) Zhang, R.; Wu, S.; Gao, F. State Space Model Predictive Control for Advanced Process Operation: A Review of Recent Development, New Results, and Insight. Ind. Eng. Chem. Res. 2017, 56, 5360−5394. (25) Zhang, R.; Xue, A.; Wang, S.; Ren, Z. An improved model predictive control approach based on extended non-minimal state space formulation. J. Process Control 2011, 21, 1183−1192. (26) Yepes, A. G.; Vidal, A.; Malvar, J.; López, O.; Doval-Gandoy, J. Tuning method aimed at optimized settling time and overshoot for synchronous proportional-integral current control in electric machines. IEEE Trans. Power Electron. 2014, 29, 3041−3054.

(27) Aström, K. J.; Murray, R. M. Feedback systems: an introduction for scientists and engineers; Princeton University Press: 2010.

N


Multilinear Model Decomposition and Predictive Control of MIMO Two

Recommend Documents