Gap-Metric-Based Multiple-Model Predictive Control with a Polyhedral

Aug 26, 2015 - Department of Automation and Key Laboratory of System Control and Information Processing, Ministry of Education of China, Shanghai Jiao...
0 downloads 6 Views 984KB Size
Subscriber access provided by UNIV OF ARIZONA

Article

Gap Metric Based Multiple Model Predictive Control with Polyhedral Stability Region xiangyuan tao, Dewen Li, Yishi Wang, Ning LI, and Shaoyuan Li Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/ie5042758 • Publication Date (Web): 26 Aug 2015 Downloaded from http://pubs.acs.org on September 18, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Gap Metric Based Multiple Model Predictive Control with Polyhedral Stability Region Xiangyuan Tao, Dewen Li, Yishi Wang, Ning Li*, Shaoyuan Li Department of Automation, Shanghai Jiao Tong University, and Key Laboratory of System Control and Information Processing, Ministry of Education of China, Shanghai 200240 Corresponding author*: [email protected]

Abstract— This paper proposes a gap metric based multiple model predictive control (MMPC) method for nonlinear system with wide operating range. The gap metric theory, integrated into a neighborhood estimation algorithm, is utilized to partition the whole operating range into sub-regions corresponding to operating points. A local controller is then designed in each sub-region, which is composed of a constant feedback control and a receding infinite horizon control with estimated polyhedral stable region. To guarantee the global stability of the whole system, we design a novel switching rule activating between local controllers. The simulation study on continuous stirred tank reactor (CSTR) is presented to validate the proposed methods. Index Terms—Multiple model predictive control; Gap metric; Switching stability; Polyhedral stable regions

I. INTRODUCTION

M

Predictive Control (MPC) has been studied for decades, while it remains open problems in both theory and application[1, 2], such as when facing nonlinear systems [3]. Nonlinear Model Predictive Control (NMPC) is feasible to solve the nonlinear system control problem [1, 3]. However, the high on-line computational burden hinders its application. Distributed MPC was proposed to reduce the computational scale and complexity of on-line nonlinear programming in NMPC, using the idea of decomposition and coordination in large scale system theory [4, 5]. Bemporad [6] presented the Explicit MPC (EMPC) algorithm that moves the on-line computational burden to the off-line. Linearization technique is also taken to design predictive controller for nonlinear systems [7], where the linear MPC theory can be directly adopted. Since the linear model obtained by linearization around a certain operating point cannot reflect the dynamic/static property of nonlinear system in the whole operating range, the controller cannot reach satisfied performance; even ensure the stability of closed-loop system. Another feasible way, to control nonlinear system with wide operating range, is the multi-model method, which utilizes a set of linear models to approximate the nonlinear system [8]. Among the multi-model based control methods, Multiple Model Predictive Control (MMPC) is widely used to control the strong nonlinear processes that are constrained, multivariable and uncertain. In [9], a multi-model predictive control framework was proposed for distributed parameter systems, where it selects and updates the model online. Using ODEL

the techniques of linear matrix inequalities (LMIs), a MMPC algorithm with a quasi-infinite horizon objective function was designed for a solution copolymerization reactor, which is modeled as a piecewise linear system [8]. Using the parallel distribution compensation method, MMPC for T-S fuzzy systems was presented in [10], where the global controller output is the fuzzy weighted integration of the local ones. Other applications of MMPC in hypersonic vehicle, nuclear stream generator, and wind turbines can be found in [11-13]. MMPC methods can greatly decrease the high on-line computational burden of NMPC and improve the approximation of whole nonlinear system using a bank of linear models instead of only one linearized model. However, the obstacles for MMPC design are obvious. 1. The selection of operating point and linear model set. The first step to design multi-model controller is to obtain the linear model set to approximate the original nonlinear system. As a model based control method, MMPC employs predictive model to predict the future system variation. The more accurate the linear model set, the more useful the predictive information the controller will have in the receding optimization phase, which brings better control performance. How to select the proper linear model set is always an open problem for MMPC. 2. System stability. Closed-loop stability is of great importance in any control strategy. However, to design a MMPC with stability guarantee is difficult. On the one hand, the receding horizon optimization of MPC leads the difficulties in stability analysis. On the other hand, MMPC system usually has more than one local controller, and the control signal is obtained by a certain kind of switching rule of local controllers. Thus the closed-loop stability not only relates to the stability of local controllers but also associates with the switching rule. These issues have been discussed in literatures. Tan et al. [14] presented a method to select operating points via gap metric in case the prior knowledge of the achievable closed-loop performance is obtained. In [15], a linear bank determination algorithm based on gap metric was proposed to construct the linear model set without a prior linear model bank. Galan et al. [16] studied a bench-scale pH neutralization system of the gap metric based control strategy where local MPC is used as one of the alternative controller designs. The multi-linear modeling philosophy is complemented with an adaptation mechanism to explain the

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

nonlinear plant behavior in the whole operating range. The multi-linear controllers exhibit good closed-loop behavior. In [17], a gap metric based closed-loop performance is defined and used to assign weights to the local controllers of a global controller. Kothare [18] analyzed the stability of a MMPC algorithm by employing multiple Lyapunov functions with strong assumptions. In [20], the authors proposed an efficient scheduled MPC for constrained nonlinear systems. It locally represents the nonlinear system as a linear time varying (LTV) system. In each sub-region, the robust MPC (RMPC) is developed. The stability of the switching control is ensured by utilizing the property of asymptotically stable invariant ellipsoid of RMPC. In this paper, we propose a gap metric based neighborhood estimation algorithm for building the linear model set. Local controller is then developed for each linear model. It is composed of a constant feedback control and a receding infinite horizon control with estimated polyhedral stable region. Then we design a MMPC algorithm for nonlinear system. We extend the controller switching rule presented in [20], which makes the proposed MMPC algorithm guarantee both control performance and stability. The paper is organized as follows. In Section II, we focus on the problem formulation. Section III presents a gap metric based neighborhood estimation algorithm. Section IV presents the steps for designing the local controller, which is composed of a constant feedback control and a modified predictive control with additional unsymmetrical state constraint obtained by gap metric theory. The complete design steps of MMPC are proposed in section V and its stability analysis is also given. This strategy is validated through CSTR model. Section VI concludes the paper.

problem. Suppose we decompose X into p sub-regions i , i  1,2,..., p satisfing i  X,1 2  ...  p  X . In each

sub-region  i , a linearization operating point OPi  ( xi , ui ) can be assigned, which satisfies equation (1). Linearize the nonlinear system (1) around operating point OPi . The following linear model  i can be obtained (2) xi (k  1)  Ai xi (k )  Bi ui (k ), x  i where xi  x  xi , ui  u  ui . Then a linear model predictive controller can be designed for each linear model  i . For MMPC strategy, the control input of whole system is acquired by either the weighting sum of each sub-controller [21] or a switching to some sub-controller [8, 18]. This paper studies the switching structure of MMPC and focuses on the following questions: 1. How to decompose the whole operating range X using the property of gap metric theory? 2. How to design local predictive controllers with infinite horizon cost function in the existence of additional unsymmetrical state constraints? 3. How to select operating points and design the controller switching rule to guarantee the whole system stable? The proposed MMPC strategy is shown in Figure 1. There are mainly 3 parts corresponding to the above 3 questions: linear model set, local MPC, and switching rule. The following 3 sections will concern them respectively. Local MPC

r

Reference Trajectory

i-th sub-operating region in X

OPi

i-th system operating point in  i

i

i-th linear model at OPi

Xi

the covered region of  i

Xi

symmetrical region in X i

i

elliptical stable region in X i

i

polyhedral invariant set in X i

Linear Model 1

MIHMPC 2

MIHMPC m

X i

y

Linear Model Set

Nomenclature Physical meaning

Controlled Plant

MIHMPC 1

.. .

Symbol

Page 2 of 12

Switching Function

...

_

e11

+

Linear Model  2

_

.. .

Linear Model  m

+

_

e2

+

em

full operating range Figure 1 The structure of proposed MMPC system

III. GAP METRIC BASED NEIGHBORHOOD ESTIMATION ALGORITHM

II. PROBLEM FORMULATION Consider a nonlinear system described as following: x(k  1)  f ( x(k ), u(k ))

(1)

where x  is the state vector, u U   is the control input vector. Let X be the full operating range of the system (1). Suppose X is large, thus a single linear model, obtained by linearization techniques around an equilibrium point, cannot fully describe the system in the whole operating range. Then we can take the multi-model strategy into account for this nx

nu

Gap metric was firstly introduced to control theory as a tool to quantify the tolerable uncertainties that preserve closed loop stability in robust control [22, 23]. Subsequently, it is applied to select operating points in multi-model methods [14, 15, 24, and 25]. In this section, we propose a neighborhood estimation algorithm based gap metric to help build the linear model set, which is the basis of designing local predictive controller and MMPC. Let P be a finite dimensional linear system. Its transfer function P( s) takes on a normalized right coprime factorization P(s)  N (s) D1 (s), with NN  DD  I where N (s)  N (s)T , D( s)  D(s)T . The gap between two systems P1 and P2 is defined as

ACS Paragon Plus Environment

Page 3 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

 ( P1 , P2 )  max{ ( P1 , P2 ),  (P2 , P1 )}

(3)

where

 D   Dj 

 ( Pi , Pj )  inf  i     Q , i, j  {1, 2}, i  j QH  Ni   N j   

The properties of gap metric are listed as follows [15]: (1) The gap metric can be regarded as a measure of distance between two linear systems. (2) The gap metric lies between 0 and 1. (3) If gap metric between two systems is small, there exists at least one feedback controller that can stabilize both systems. Considering these properties, we propose an algorithm to estimate a neighborhood region around the operating point OPi  ( xi , ui ) in which the linearization model can approximate the nonlinear system. Algorithm 1. (Neighborhood estimation) Step 1. Linearize the nonlinear system (1) around the chosen operating point OPi , and get linear model  i . Step 2. Choose a sufficient large state space X s ,i  X around the operating point OPi and distribute N ( N is a large number) reference points in X s ,i . Compute the Jaccobi matrices  Ai , j Bi , j  , j  1, 2,..., N of system (1) in these reference points and get N linear model Si , j , j  1,2,..., N satisfying Si , j : xi (k  1)  Ai , j xi (k )  Bi , j ui (k ) . Step 3. Compute the gap metric between Si , j and  i

 j  gap(Si , j , i ), j  1,2,..., N Step 4. Prescribe a threshold level  . Cluster the reference points that satisfy  j   , then the sub-region X i can be constructed. Therefore, we can decompose the whole operating range X into sub-regions X i using Algorithm 1. Remark 1. The obtained sub-region X i is an irregular shape that can be approximated by a polyhedron as X i ={x M i ( x  xi )  di } . Remark 2. In most cases, the sub-region X i constructed by Algorithm 1 is a connected domain that includes the operating point OPi . If more than one connected domain is obtained, only the region that includes the operating point is needed. Remark 3. The selected state space X s ,i is large enough to cover the neighborhood of the operating point OPi . It is suggested to be selected using the prior knowledge of the system, which is a big range to include the operating point OPi . After X s ,i is given, N could be selected since it is only related to the discretized sampling interval. The threshold level  should be set appropriately. A high threshold will lead a large

sub-region where a single linear model may not be sufficient to span. On the other hand, conservative approximated sub-region will be obtained if the threshold is set too low. Therefore, the whole system will be divided to more sub-regions, and will be described by more submodels, which leads to more frequent switching. Therefore the selection of the threshold is definitely a trade-off. We suggest determining its value case by case. Referring to [15, 24], the empirical value or the initial value could be selected as 0.4    0.6 . These issues will be explained quantitatively combined with local controller design in section IV through a simulation example. IV.GAP METRIC BASED LOCAL PREDICTIVE CONTROLLER DESIGN

A. MPC with Infinite Horizon cost function (IHMPC) Consider the tracking problem that drives any state in region X i obtained by Algorithm 1 to operating point OPi . The aim of the predictive control algorithm is to find a state feedback control law ui (k  j k )  Ki xi (k  j k ) that can converge xi to origin. The MPC is designed as the following, with infinite horizon cost. min J i (k ) ui ( k  p | k ), p  0,1,..., 

 x (k  j k )   0   xi (k  j k )  J i (k )    i     i  0 ui ( k  j k )   0 R   ui (k  j k )  where   0, R  0 and

(4)

ui ,r (k  i k )  ui ,r ,max , i  0, r  1,2,..., nu

(5)

xi ,r (k  i k )  xi ,r ,max , i  0, r  1,2,..., nx

(6)

T



Kothare [26] solved the above optimization problem using LMIs to minimize the upper bound  of the objective function, and proved the closed-loop stability with an estimated elliptical stable region. Lemma 1. [26] (MPC with Infinite Horizon cost function, IHMPC) For the linear system (2), the state feedback control law that minimizes the upper bound  of the objective function J i (k ) and robustly stabilize the closed system with an estimated elliptical stable region T 1 i  {x ( x  xi ) Qi ( x  xi )  1} given by

ui (k  i k )  Fi (k ) xi (k  i k ), Fi (k )  Yi (k )Qi1 (k ) , where Yi and Qi are obtained by solving the following problem:

min

 i ,Qi ,Yi , X i , Zi

i

(7)

s.t.  1 xi (k |k )T     0, Qi  0 Qi   xi (k |k)

ACS Paragon Plus Environment

(8)

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Qi  * *   *

Qi AiT  Y T BiT Qi * *

Qi 1/ 2 0 I *

Y T R1/ 2   0   0, 0    I 

Zi  Qi  0, Zi ,rr  xi2,r ,max , r  1,..., nx

Algorithm 2. (Modified-IHMPC, MIHPC) (9)

 X i Yi  2 (11)   Q   0, X i , rr  ui , r ,max , r  1,..., nu i  The matrix inequalities (10) and (11) help handle the symmetrical state constraint (5) and the input constraint (6). Here we need to design the controller within the neighborhood region X i obtained by Algorithm 1, yet the state constraint is most likely asymmetrical. Thus, to apply Lemma 1, the state constraint needs to be preprocessed. We take the following method to deal with the asymmetrical constraints. 1. Find a symmetrical region X i ={x M i ( x  xi )  di } in X i . 2. Define virtual output yi ,r  M i ,r xi , r  1, 2..., ny , where M i , r is the r th row of M i , and n y is the row of matrix M i . 3. Re-express the symmetrical region X i as symmetrical output constraints (12)

th

where di , r is the r row of d i . 4. Based on the derivation in [27], (12) can be satisfied by adding the following LMI into optimization problem (7). Wi M i ( Ai Qi  BiYi )  2    0,Wi , rr  yi , r ,max , r  1,..., ny (13)  Q i   where Wi is the optimization variable. The above steps turn the asymmetrical state constraints into the symmetrical output constraints that can be handled by LMI. However, this inevitably leads to conservation since the symmetrical region X i sometimes is much smaller than X i .

where the linear model (2) can approximate nonlinear system (1) using Algorithm 1. Step 2. Estimate the maximum asymptotically stable ellipsoidal region i  {x | ( x  xi )T Gi 1 ( x  xi )  1} of IHMPC and obtain the corresponding state feedback gain Fi  Yt ,i Gi1 , where Gi , Yt ,i are the optimal solution of Qi and Yi respectively by solving the following maximization problem [18], max log det(Qi )  i , Qi ,Yi , X i , Z i ,Wi (14) s.t. (9)  (11), (13) Step 3. For state feedback gain Fi obtained from Step 2, the polyhedral invariant set i  {x Si ( x  xi )  hi } can be constructed by the following steps [29]: (3.1) Set m  1 and initialize Si  [M i , FiT ,  FiT ]T , T T hi  [di , umax  uiT , (umin  uiT )] .

(3.2) Select row m of (Si , hi ) labeled as ( Si , m , hi , m ) and solve the following problem to check if the constraint is satisfied at next sampling time max Wm x

(15)

s.t

 Wm  Si , m ( Ai  Bi Fi )( x  xi )  hi ,m  Si ( x  xi )  di   If Wm  0 , add the constraint to (Si , hi ) by assigning T

Si  [SiT ,(Si , m ( Ai  Bi Fi ))T ]T and hi  [hT , hiT, m ]T . (3.3) Let m  m  1 and go to step (3.2) until m is larger than the number of rows in (Si , hi ) . On-line:

B. Modified MPC with Infinite Horizon cost function (MIHMPC) To compensate for the performance degradation caused by handling asymmetrical constraint, this section proposes a modified local controller design strategy using polyhedral invariant sets. Definition 3.1. [28] For autonomous system x(k  1)  f ( x(k )), x  X , if set    x Sx  h  X

Off-line: Step 1. Estimate the neighborhood region X i ={x M i ( x  xi )  di } around the operating point OPi ,

(10)

yi ,r (k  i k )  di ,r , i  0, r  1,2,..., ny

Page 4 of 12

satisfies:

whenever

x( k )  

,

x(k  1)   , i  1, 2,...,  , then  is a polyhedral invariant set.

At every sampling time, the control law is  F u (k ), x(k )  i   i ui (k )   i i  IHMPCi , x(k )   i

(16)

where ( i  i ) represents the area that is within region  i but outside region  i , IHMPCi denotes the controller designed by Lemma 1 for operating point OPi . The control law (16) can be proved stabilizing the closed-loop system by the following theorem. Theorem 1 For linear system (2), the MIHMPC, expressed as Algorithm 2, can asymptotically stabilize the closed-loop system with an estimated polyhedral region of stability  i .

ACS Paragon Plus Environment

Page 5 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Proof. Case 1: If x(k ) i  i , the constant feedback control law ui (k )  Fi xi (k ) will work. Consider Lyapunov function as following (17) Vi (k )  xi (k )T Px Pi   i Qi 1 i i (k ), where  i and Qi are obtained by solving (7) for linear system

(2). The satisfaction of (9) for constant feedback matrix Fi ensures [ Ai  Bi Fi ]T Pi [ Ai  Bi Fi ]  Pi  Gi  FiT RFi  0 (18)



q E (C Af  C A )  k0 exp( )C A V RT (21)  q H E UA T  (T f  T )  k0 exp( )C A  (Tc  T ) V C p RT V C p The input constraints is (22) umin  Tc  umax . where umin  250, umax  500 . The physical meanings of these parameters are listed in TABLE I. TABLE I CA 

Multiply xi (k )T and xi (k ) at the left and right side of (18) respectively and substitute xi (k +1) for [ Ai  Bi Fi ]xi (k ) . We have:

xi (k  1) Pi xi (k  1)  xi (k ) Px i i (k )   xi (k ) [Gi  Fi RFi ]xi (k ) (19) The inequality (19) reveals (20) Vi (k  1)  Vi (k ) So Vi (k ) is a strictly decreasing Lyapunov function. Besides, T

T

T

T

the set of initial state in polyhedral region  i constructed by Algorithm 1 guarantees all future states to stay within  i . Thus for initial state in ( i  i ), the constant feedback control law can asymptotically stabilize the closed-loop system. Thus after finite sampling time, the system state will definitely enter  i to approach xi . Case 2: If x(k )  i , the predictive controller IHMPCi will be activated. Based on Lemma 1, the closed-loop system is asymptotically stable. Remark 4. Theorem 1 can be proved for linear system (2). According to Algorithm 1, the linear model is a good representative of the nonlinear system (1) in the polyhedral region i  X i , so

MODEL PARAMETERS AND CONSTRAINTS Symbol

Physical meaning

CA

Concentration of A

0  CA  1 mol / l

T

250K  T  500K

C Af

Reactor temperature Temperature of coolant stream Flow Volume of reactor Feed concentration

Tc q V

Example I. CSTR The following simulation on a benchmark continuous stirred tank reactor (CSTR) shows the detailed steps to apply the Algorithm 1&2. The dynamic model of the CSTR system, shown as Figure 2, can be described by the following differential equations (21) based on material and energy balance [18].

CAf , T f , q

250K  TC  500K 100 l / min 100 l 1 mol / l

k0

Time constant

4.71108 min 1

H

Reaction heat

-2 105 J / mol

UA

Heat exchange coefficient Density of liquid

105 J / min/ K

C

Thermal capacity

1J / g / K

E/R

Reaction activated energy

8000 K



103 g / l

Let state vector x  CA T  , u  Tc . Discretize (21) with T

sampling time 18s . The control aim is to drive the current state to the operating point OP1  [ x1  [0.522,398.52]T , u1  302] whose linearization model 1 is

 i can also be regarded as the stable region for the

closed-loop system composed of the nonlinear system (1) and the control law (16). C. Simulation Examples

Value/constraint

x(k  1)  A1 x(k )  B1u(k )

(23)

 1.095e  5 0.942 0.7337e  3 where A1    , B1   0.03129  . 5.615 1.086    

Next we will illustrate the implementation of the proposed Algorithm 1&2 step by step. Step 1: Estimate the neighborhood region X 1 (1) Choose a large region X s ,1 around the OP1 within the range of state constraints, shown in Table I. For this example, We select X s,1  x 0.12  CA  0.92,348.8  T  448.8 . Within X s ,1 , discretize the state CA and T using sampling interval 0.1

T

Tc

CA , q Tc

and 5 respectively. Then we have N=9*21=189 reference points, which are the 189 intersections of the grid shown in Figure 4. Computation of the Jaccobi matrix around these reference points results in 189 linear models S1, j , j  1,2,...,189 . (2) Compute the gap metric  j  gap(S j , 1 ), j  1,2,.., N shown

Figure 2. CSTR system

in Figure 3.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

In other words, the sub-region X1 can be described by one linearization model 1 (23). Step 2: Estimate the maximum asymptotically stable ellipsoidal region 1 of IHMPC.

1 0.8

(1) Find a symmetrical polyhedron X i1 ={x M1 ( x  x1 )  d1} to

0.6 Gap

approximate neighborhood region X 1 , where

0.4

 1.32 0.02   4.22 0.135  , M1    1.73 0.035   18.16 0.15 

0.2 0 450 1

400

0.8 0.6

350 T

0.4 300

0.2 0

CA

Figure 3. The gap between the 189 reference models and 1

(3) Choose threshold level. The smaller the threshold  , the smaller the neighborhood of OP is. Therefore, the whole system will be described by more submodels, which leads to more frequent switching. But if the threshold  is too big, the approximation of submodel will be worse, which leads to worse control performance. In this example, we tried several thresholod to compare the neighborhood. When   0.3 , the neighborhood of OP is composed of 10 reference points; when   0.4 , the neighborhood includes 17 reference points; when   0.5 , the number becomes 27, and when   0.6 , there are 40 reference points in the neighborhood. In pursuit of balance between approximation and switching effect, we choose threshold level   0.4 in this example. X1 ={x M1 ( x  x1 )  d1} is then constructed as the green region around the red point in Figure 4.

(2)

Define

virtual

0.2256 1.5228   d1    1.528     1.857  y1   y11 y12

output

y13 y14 

T

where y1  M1 x1 , then y1,max  d1 . (3) Estimate the maximum asymptotically stable ellipsoidal region 1  {x | ( x  x1 )T G11 ( x  x1 )  1} of IHMPC by solving optimization problem (7) and get the corresponding F1   406.93 6.58 feedback gain , where 0.015 0.9  G1   .  0.9 116.74 

Step 3. Estimate polyhedral 1  {x S1 ( x  x1 )  h1} of F1 .

invariant

set

T T (1) Initialize matrices S1 and h1 with S1  [ M1, F1 ,  F1 ],

T T h1  [d1 , umax  u1T , (umin  u1T )] .

(2) The polyhedral invariant 1 can be obtained after 43 iterations based on Algorithm 2 which is showed as the region included by blue lines in Figure 5. 440 430 X: 0.375 Y: 411.1

420

X

410 400



T

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 12

390

1

1



1

380 370 360 350 0.1

0.2

0.3

0.4

0.5 CA

0.6

0.7

0.8

0.9

Figure 5. The regions of X1 , 1 , 1 for OP1

Figure 4. The approximate region of linear model 1

where

 11 0.05   2.3068  19 0.13  2.2038     M 1   13 0.42  , d1  5.0733     0   10  3.770   33  5.5791 0.5 

Hereto, the design of MIHMPC for operating point OP1 is finished. To show the effectiveness of the designed controller, we make the controller drive the initial state x0  0.375 411.1 to x1 . The state transition trajectory is presented in Figure 5. Figure 6 gives the changes of state variables and control input. The simulation results verify the good control performance in estimated polyhedral stable region 1 . Since the initial state x0 is not in the elliptical stable region 1 of IHMPC, thus IHMPC cannot drive it

ACS Paragon Plus Environment

Page 7 of 12



to x1 . This manifests the advantage of MIHMPC because of an enlarging feasible region.

.

x1  



vt vt x1  u Lt0 lt0

(24)



vt x2  x1 Lt0 .

412 410



.

x3 

408 406



vt vt sin[ x2  x1 ] t0 2L

T(K)

where

404

x1 : angle difference between truck and trailer

402

x2 : horizontal position of rear end of trailer x3 : vertical position of rear end of trailer u : angle of driver The parameters of this model are given as v  1 ,

400 398 0

50

100

150 Sampling time

200

250

300

(a) State variable -- T

_

l  2.8 , L  5.5 , t0  0.5 , t  2 ,   u   . Assume state

0.55

vector x  [ x1 , x2 , x3 ]T . Discretize (24) with sampling time 0.03s. Let the initial state be x0  [3,0.5,3]T .The aim is to drive the initial state to the equilibrium state xeq  [0, 0, 0]T , ueq  0 .The linearization model  is

C (mol/L)

0.5

A

0.45

x(k  1)  Ax(k )  Bu(k )

0.4

(25)

Where 0.35

0

50

100

150 Sampling time

200

250

0 1.022  A   -0.008823 1  0.02259 -0.12 

300

(b) State variable-- C A 305

0  0 , 1 

 -0.04333    B   0.0001884   -0.0004785   

As mentioned above, we will illustrate the implementation of Algorithm 1&2 step by step. Step 1: Estimate the neighborhood region X (1)Choose a large region

300

295 C

T (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

290

X s  {3  x1  3,  3  x2  3,  4  x3  4} around xeq . Within

285

X s , discretize the state x1 , x2 and x3 using sampling interval

280

0

50

100

150 Sampling time

200

250

300

1,0.5 and 1 respectively. Thus we have N=7*13*9=819 reference points. Computation of the Jaccobi matrix around

(c) Control variable-- Tc Figure 6. Simulation results of MIHMPC on CSTR

Example II. Truck-trailer In order to illustrate the proposed methods, we apply the above technique to another example-truck-trailer with more states [19]. The model can be described as following differential equations.

these reference points results in 819 linear models Si , i=1,2,...,819. (2)Compute the gap metric  j  gap(S j , ), j  1,2,.., N (3)Choose

threshold

level   0.4 .

The

sub-region

X ={x M ( x  xeq )  d} is showed as the blue region in Figure 7.

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

feedback

F  [1.3061, 3.0789,0.4646]

gain

,

where

invariant

set

 9.0730 1.7039 -0.2039    G   1.7039 0.9045 0.6429  .  -0.2039 0.6429 7.2723   

5

x3

Step

polyhedral   {x S ( x  xeq )  h} of F .

0

3. Estimate

(1) Initialize matrices S and h with S  [M , F T , F T ], -5 2 1

T T h  [d , umax  uT , (umin  uT )] .

4 2

0

(2) The polyhedral invariant  can be obtained based on

0

-1 x2

-2 -2

-4

x1

Algorithm 2 which is showed as the yellow region in Figure

Figure 7.The approximate region of linear model



8.

After uniformization of d , we have

 0   0  0.2000   0  -0.2000 M   0  0.1429   0.3333   -0.1429  -0.3333 

1   1 1   1 1 d    1 1   1 1   1  

0 0.2500   0.6667 0  -0.8000 0   -0.6667 0  0.8000 0   0 -0.2500  -0.8571 0   0 0   0.8571 0  0 0 

5

x3

0

-5 2 1

4 2

0

x2

x1

Figure 8.The regions of



The initial state x0  [3,0.5,3] is then driven to xeq .The

(1) Find a symmetrical polyhedron X1 ={x M ( x  xeq )  d } to

1 0 -0.2500     -1.0000 0  1 _ 1 , 0 0  d    1.0000 0  1 1  0 0    0 0.2500  1 _

_

_

_

state trajectory is presented in Figure 9. The change of state variables and control input are presented in Figure 10.

approximate neighborhood region X 1 , where

_

-4

T

ellipsoidal region  of IHMPC.

_

-2 -2

Step 2: Estimate the maximum asymptotically stable

0   0.1667  _  0.3333 M   -0.1667  -0.3333  0 

0

-1

4 3 2 1 x3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 12

0 -1

_

(2) Define virtual output y  [ y1 , y2 , y3 , y4 , y5 , y 6 ]T ,

-2 -3

_

where y  M xeq , then ymax  d .

-4

2 1

0

0 -1

(3) Estimate the maximum asymptotically stable ellipsoidal x2

region

1

  {x | ( x  xeq ) G ( x  xeq )  1} of T

IHMPC

by

solving optimization problem (7) and get the corresponding

ACS Paragon Plus Environment

-2 x1

(a) Front view

Page 9 of 12

controller. Here we extend this idea and propose the gap metric based MMPC with stability guaranteed.

-3

-2

Algorithm 3. (Gap Metric based MMPC, GM-MMPC) Step 1. Choose the desired state xt as the first operating

x1

-1

point OP1 and set i  1 . Step 2. Design local controller MIHMPCi, with the polyhedral invariant set  i of feedback gain Fi and

0

1

asymptotical stable region  i of IHMPCi .

2

Step 3. If the initial state satisfies x0 i , choose the next 3

-1

0 x2

operating point OPi 1  ( xi 1 , ui 1 ) satisfying xi 1 i - i , which has the smallest Euclidean distance to the desired state among all reference points in (i - i ) . Set i  i  1 and go to step 2. If

1

(b)Top view Figure 9.The state transition trajectory

the initial state satisfies x0 i , set p  i and terminate. Through Algorithm 3, several operating points are selected. The polyhedral stable regions for the corresponding local controllers overlap with each other. Employing these adjacent polyhedral stable regions, the controller switching rule can be designed as follows, which guarantees the controller switching stability:

4 x1 x2 x3

3 2

state

1 0

Controller switching rule. For x(k ) i 1  i , i  1,2..., p 1 ,

-1 -2

-4

100

200

300

400 500 Sampling time

600

700

800

3

2

1

0

-1

-2

0

i +1  i

region  i , the local controller MIHMPCi+1 will drive the state 0

(a)State variables

-3

where

represents the area that is within region i +1 but outside

-3

u

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

200

400

600

800 1000 1200 Sampling time

1400

1600

1800

2000

(b)Control variable u Figure 10.Simulation results of MIHMPC on truck-trailer

V. GAP METRIC BASED MMPC FOR NONLINEAR SYSTEM A. Gap metric based MMPC Using the above MIHMPC for sub-region, this section presents MMPC strategy for a nonlinear system with wide operating range. [20] proposed a multi-model control method, where the stable region of sub-controllers nested with each other. The switching stability can be guaranteed by selecting the operating point in the stable region of adjacent local

to region  i . Once the state reaches the boundary of  i , the controller immediately switches to MIHMPCi. The flow chart of the Algorithm 3 can be explained in Figure 11. As can be seen, the control system is a two-layer switching structure. The switching in outer layer determines which local controller should be activated based on current state. In the inner layer, switching control law (16) is applied to the nonlinear system. For nonlinear system (1), it can be proved x(k ) ip0 i , the proposed MMPC can drive the state to desired one xt in finite time. Theorem 2. For the nonlinear system (1), the GM-MMPC in Algorithm 3 with two-layer switching has an estimated stable region ip1 i around the desired state xt . Proof The linear model set constructed in Algorithm 3 can be combined into a piecewise linear system: : {1 , 2 ,...,  p } (23) where i (i  1,.., p) is given as (2).

For system (23),

if x(k ) i 1  i, i {1,..., p 1} , system state will gradually close to xi 1 using local controller MIHMPCi+1according to Theorem 1. Since xi 1 i (See Step 3 in Algorithm 3), the state transition trajectory will enter  i in finite time. The process will continue along with the switching of controller until the state is driven into the stable region 1 where local

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

controller MIHMPC1 is responsible to take state to the desired state xt . According to Remark 4, the piecewise linear system (23) well describes the nonlinear system (1) in the region ip1 i , so ip1 i is also the stable region for the closed-loop system composed of the nonlinear system (1) and proposed MMPC. x(0)  x0 , k  0

No

x(k )  1?

Yes

Outer layer

x(k )   2?

i=1

... ...

No

x(k )   p?

No

Yes

i=p

x(k )  i ?

No Inner layer

corresponding elliptical regions 1 ~ 4 constructed by Algorithm 3. Figure 12 (b) provides the information of local controller switching of GM-MMPC. The initial state x0  [0.189, 432.083]T was moved by IHMPC4 and entered the stable sub-region  3 after 7 samples. Similarly, local Controller MIHMPC3 and MIHMPC2 drove the state to stable sub-region  2 and 1 , respectively. Local controller MIHMPC1 ensures the convergence of the state to the desired one xt  [0.52,398.77]T . The selected operating points of GM-MMPC are as follows: [0.375, 411.07]T  xt  [0.522,398.77]T . The operating points,

except the initial and the desired states, are selected within (i - i ) . To avoid confusion, we did not plot the 6 elliptical regions of EMMPC in Figure 12 (a). The selected operating points of EMMPC are as follows: x0  [0.189, 432.083]T  [0.208, 429.29]T  [0.253, 423.48]T  [0.32, 416.26]T  [0.357, 412.7]T  [0.412, 407.86]T 

Yes u (k )  IHMPCi u (k )  Fi x(k )

while EMMPC and SMPC have 6 and 10 ones respectively. Figure 12 (a) shows the nested stable regions 1 ~4 and

x0  [0.189, 432.083]T  [0.215, 428.3]T  [0.266, 422]T 

Yes

i=2

Controller # I

(MIHMPC)

x(k  1)  f ( x(k ), u (k )) k  k 1

xt  [0.522,398.77]T . Here the operating points, except the

initial and the desired states, are selected within  i . For SMPC, the operating points changes as: x0  [0.189, 432.083]T  [0.2, 430.43]T  [0.22, 427.66]T 

[0.243, 424.72]T  [0.272, 421.32]T  [0.305, 417.78]T  [0.34, 414.31]T

No

Page 10 of 12

x(k )  xt ?



[0.38, 410.6]T



[0.425, 406.71]T



[0.478, 402.3]  xt  [0.522,398.77] . As shown in Table II T

T

Yes Stop Figure 11. Online implementation of GM-MMPC algorithm

B. Simulation Example In this part, we will implement Algorithm 3 into the nonlinear state transition process of CSTR, as shown in Figure 2. The initial state is x0  [0.189, 432.083]T and the control objective is to regulate the initial state to the desired state xt  [0.522,398.77]T by manipulating Tc . Four local MIHMPCs are designed since the termination condition of Algorithm 3 satisfies x0 4 . Here the proposed algorithm GM-MMPC is compared with the scheduled model predictive control (SMPC) proposed in [20] and the MMPC using IHMPC with the elliptical invariant set only (EMMPC). For SMPC, with the polyhedral invariant set the LTV system was utilized to represent the nonlinear system in the neighborhood of an operating point. For EMMPC, the algorithm is almost the same as GM-MMPC except it takes IHMPC as the local controller and the next operating point is not selected within (i - i ) , but within the elliptical region  i . In this simulation, GM-MMPC uses 4 local controllers,

and Figure 12 (c)(d)(e), the proposed GM-MMPC strategy has the least local controllers, but has better control performance in terms of the rapidity and smoothness of state transition among all three methods. In Figure 12 (e), we illustrate the control inputs by the stair figure with a partial enlarged sub-figure, which helps show the input change clearly. The comparison of the control input shown also verifies the superiority of the proposed MMPC with less fluctuation in manipulating variable. TABLE II COMPARISION OF SIMULATION RESULTS

Algorithm Algorithm 3 EMMPC SMPC

ACS Paragon Plus Environment

Number of local controllers 4 6 10

Settling time (sampling) 231 400 448

Page 11 of 12

Tc(K)

350

305

340

300

330

295

320

290

310

285 20

GM-MPC EMMPC SMPC

25

30

35

40

300 290 280 270 260

(a). State transfer trajectory

250 4

0

50

100

Sub-region number 3.5

Region/Controller number

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Sub-controller number

150

200 250 300 Sampling time

350

400

450

500

(e) Change of Control Variable -- TC

3

Figure 12. Simulation results of GM-MMPC on CSTR

2.5

VI.CONCLUSIONS

2

MIHMPC2

IHMPC4

1.5

F3

1 0.5

IHMPC3

F2

0

20

F1

IHMPC2

0 40

60

80

time(min)

(b) Switching sequence

(c).Response of state variable-- C A

IHMPC1 100

120

This work presents a stabilizing multi-model predictive control algorithm for nonlinear systems with wide operating region. Using the property of gap metric, the whole operating range is decomposed into multiple sub-regions, where linear models can then be built. The local controller MIHMPC for linear model is designed as a composition of a constant feedback control and a MPC with infinite cost function, which reduces the conservation derived from turning the asymmetrical state constraints to symmetrical ones. The GM-MMPC with switching between overlapping polyhedral invariant sets guarantees the closed-loop stability. As for the online computational burden, the proposed MMPC is also efficient since only a linear optimization problem with a few constraints needs solve online. As pointed in [18],a common limitation of multiple model approaches is the stability results are only valid provided that the nonlinear system can be described by linear piecewise systems. However, the basis of our work is gap metric theory that can help select linear models to well represent the nonlinear system, so the stability results of the proposed algorithm are also reliable to nonlinear systems. Our work did not consider any uncertainties. However, the infinite horizon MPC can be extended to handle uncertainties modeled as a polytopic uncertainty description [26], so it is possible to extend our work to uncertain systems which is also the focus of future study. ACKNOWLEDGMENT

(d).Response of state variable-- T

This work is supported by the National Nature Science Foundation of China (61374109, 61304078, 61473184), the National Basic Research Program of China (973 Program-2013CB035500), and National High Technology R&D Program of China (863 Program - 2015AA043102).

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

REFERENCES [1]

[2]

[3]

[4]

[5]

[6]

[7] [8]

[9]

[10]

[11] [12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22] [23]

[24]

Mayne, D.Q. Nonlinear model predictive control: Challenges and opportunities, Progress in systems and control theory. Springer Basel AG, vol. 26, pp. 23-44, 2006 Camacho, E.F.; Bordons, C. Nonlinear model predictive control: An Introductory Review, Assessment and future directions of nonlinear model predictive control. Lecture Notes in Control and Inform, Springer, Berlin, vol. 358, pp. 1-16, 2007 Qin, S.J.; Badgwell, T.A. A survey of industrial model predictive control technology. Control engineering practice, vol. 34, pp. 469-475, 2003 Du, X.N. ; Xi, Y.G. ; Li, S.Y. Distributed model predictive control for large-scale systems. Proceedings of the 2001 American Control Conference. Arlington, VA: IEEE, pp. 3142-3143, 2001 Camponogara, E. D.; Jia, B.H.; Krogh, S. Talukdar. Distributed model predictive control. IEEE Control Systems Magazine, vol. 22(1), pp. 44-52, 2002, Bemporad, A.; Morari, M. ; Dua, V. ; Pistikopoulos, E.N. The explicit linear quadratic regulator for constrained systems. Automatica, vol. 38(1), pp. 3-20, 2002 Murray-Smith, R.; Johansen (Eds), T.A. Multiple model approaches to modeling and control. London: Taylor and Francis. 1997 Özkan, L.; Kothare, M.V.; Georgakis, C. Control of a solution copolymerization reactor using multi-model predictive control. Chemical Engineering Science, vol. 58, pp. 1207-1221, 2003 Garcia, M.R.; Vilas, C.; Santos, L.O.; Alonso, A.A. A robust multi-model predictive controller for distributed parameter systems. Journal of Process Control, vol. 22, pp. 60-71, 2012 Li, N.; Li, S.Y.; Xi, Y.G. Multi-model predictive control based on the Takagi-Sugeno fuzzy models: a case study. Information Science, vol. 175, pp.247-263, 2004 Hua, C.; Li, N.; Li, S.Y. Switching multi-model predictive control for hypersonic vehicle, 8th Asian Control Conference, pp. 677-681,2011 Hu, K.; Yuan, J.Q. Multi-model predictive control method for nuclear steam generator water level. Energy conversion and management, vol. 49, pp. 1177-1184, 2008 Mostafa, S.; Malik, O.P.; Westwick, D.T. Multiple model predictive control for wind turbines with doubly fed induction generators. IEEE Transactions on Sustainable Energy, vol. 2, pp. 215-225, 2011 Tan, W.; Marquez, H.J.; Chen, T.W.; Liu, J.Z. Multi-model analysis and controller design for nonlinear processes. Computers and Chemical Engineering, vol. 28, pp. 2667-2675, 2004 Du, J.; Song, C.; Li, P. Application of gap metric to model bank determination in multi-linear model approach. Journal of Process Control, vol. 19, pp. 231–240, 2009 Galan, O.; Romagnoli, J.A.; Palazoglu, A. Real-Time Implementation of Multi-Linear Model-Based Control Strategies – An Application to a Bench-Scale pH Neutralization Reactor. Journal of Process Control, 2004, 14(5): 571-579 Arslan, E.; Çamurdan, M.C. ; Palazoglu, A.; et al. Multi-model Scheduling Control of Nonlinear Systems Using Gap Metric. Industrial & Engineering Chemistry Research, 2004, 43 (26): 8275–8283 Özkan, L.; Kothare, M.V. Stability analysis of a multi-model predictive control algorithm with application to control of chemical reactor. Journal of Process Control, vol. 17, pp. 81-90, 2006 Cao Y Y Frank PM.Stability analysis and synthesis of nonlinear time—delay systems via linear Takagi-Sugeno fuzzy models[J],Fuzzy Sets and Systems,200 1,1 24:2 1 3—229. Wan, Z.; Kothare, M.V. Efficient scheduled stabilizing model predictive control for constrained nonlinear systems. International Journal of robust and nonlinear control, vol. 13, pp. 331-346, 2003 He, M.; Cai, W.J.; Li, S.Y. Multiple fuzzy model-based temperature predictive control for HVAC systems. Information Sciences, vol. 169, pp. 155-174, 2005 Zames, G.; El-Sakkary, A.K. Unstable systems and feedback: The gap metric. Proceeding of Allerton Conference, pp. 380-385, 1980 El-Sakkary, A.K. The gap metric: robustness of stabilization of feedback systems. IEEE Transactions on Automatic Control, vol. 30, pp. 240–247, 1985 Hariprasad, K.; Bhartiya, S.; Gudi, R.D. A gap metric based multiple model approach for nonlinear switched systems. Journal of Process Control, vol. 22. pp. 1743-1754, 2012

Page 12 of 12

[25] Galán, O.; Romagnoli, J.A. ; Palazoglu, A. ; Arkun, Y. Gap metric concept and implications for multi-linear model-based controller design. Industrial & Engineering Chemistry Research, vol. 42, pp. 2189-2197, 2003 [26] Kothare, M.V.; Balakrishnan, V.; Morari, M. Robust constrained model predictive control using linear matrix inequalities. Automatica, vol. 32, pp. 1361–1379, 1996 [27] Boyd, S.; El Ghaoui, L.; Feron, E.; Balakrishnan, V. Linear matrix inequalities in system and control theory. SIAM, Philadelphia,1994 [28] Pluymers, B.; Rossiter, J.A.; Suykens, J.A.K; De Moor, B. The efficient computation of polyhedral invariant sets for linear systems with polytopic uncertainty. Proceedings of the American Control Conference, Portland, USA, pp.804-809, 2005 [29] Bumroongsri, B.; Kheawhom, S. An off-line robust MPC algorithm for uncertain polytopic discrete-time systems using polyhedral invariant sets. Journal of Process Control, vol. 22, pp. 975-983, 2012

ACS Paragon Plus Environment