A Bilevel Optimization Approach to Coordination of Distributed Model

Subscriber access provided by UCL Library Services

Article

A Bi-level Optimization Approach to Coordination of Distributed Model Predictive Control Systems Bardia Hassanzadeh, Jinfeng Liu, and J. Fraser Forbes Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b02414 • Publication Date (Web): 24 Dec 2017 Downloaded from http://pubs.acs.org on December 30, 2017

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 44 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

A Bi-level Optimization Approach to Coordination of Distributed Model Predictive Control Systems Bardia Hassanzadeh,∗ Jinfeng Liu, and J. Fraser Forbes Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada, T6G 1H9 E-mail: [email protected]

Phone: +1(780)964-5141. Fax: +1(780)492-2881

Abstract In this paper, a novel algorithm is proposed to coordinate distributed model predictive control (DMPC) systems. The main goal of coordinated distributed model predictive (CDMPC) is to improve the performance of an existing decentralized MPC system, which neglects the interactions among subsystems. In the first scenario, an analytic approach derives the global optimal solution to the plant-wide problem in the absence of inequality constraints in local controllers. In the second scenario, an iterative approach determines a local optimal solution to the general CDMPC problem based on the method of feasible directions. The proposed CDMPC algorithm is capable of stabilizing open-loop unstable dynamics.

1

Introduction

Traditionally, decentralized MPC and centralized MPC are two candidates for designing the control structure in networked process units. Throughout this paper, decentralized MPC is ∗

To whom correspondence should be addressed

1

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

referred to a specific type of network when the interactions among subsystems are ignored; and, centralized MPC is referred to a type of network in which a single controller considers all the interactions. While the decentralized MPC is straightforward to implement, it may lead to degradation of plant-wide performance or even loss of closed-loop stability. This is due to the fact that the interactions between sub-systems are neglected in the decentralized control fashion. On the other hand, centralized MPC is expected to give the highest performance; however, it may become too complicated to implement as the size of the control problem grows. These considerations have motivated significant efforts in the control of networked systems using distributed MPC (DMPC) strategies. One important class of DMPC systems is called cooperative DMPC; some of the important recent work in this category includes the metho cooperative DMPC approaches proposed by Rawlings and Stewart 1 and Stewart et al. 2 , as well as the Lyapunov-based sequential DMPC approaches proposed by Liu et al. 3 ,4 . Another important type of DMPC, pertaining to the present work, is coordinated DMPC (CDMPC), which strives to reach the optimal plant-wide performance by adding a coordination level to a currently installed decentralized network. In a CDMPC scheme, distributed MPCs communicate with a coordinator to achieve improved performance. Different algorithms have been developed for CDMPC including price-driven approach, 5–8 and prediction-driven approach. 9 In CDMPC schemes, the local MPC controllers communicate with a coordinator to achieve the centralized MPC performance. In this scheme, an existing network of decentralized MPC controllers are modified to exchange information with an upper-level coordinator. On the other hand, in cooperative DMPC schemes, the distributed controllers communicate with each other to improve the overall performance towards the centralized MPC optimal trajectory. In general, cooperative MPC schemes are not applied to an existing decentralized MPC network, because a centralized model is required inside each local controller. The scope of this paper is limited to CDMPC design based on price-driven coordination. This method was first introduced by Cheng et al. 5 , and then extended by Marcos et al. 7 .

2


Page 2 of 44

Page 3 of 44 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


The price-driven CDMPC is essentially a primal-dual problem which can be regarded as a bi-level programming (BLP) problem. The two levels of this BLP problem consist of local MPC controllers as the lower-level problem, and a coordinator as the upper-level problem. In practice, local optimization algorithms are used to solve the CDMPC problem, due to real-time implications imposed by global optimization algorithms. Conventionally, the pricedriven CDMPC scheme is solved using a nested approach. 7 In this approach, MPC controllers solve their local problems for a fixed value of price; then, the local controllers communicate with a coordinator to update the price vector via Newton’s method. This procedure continues until an optimal price is found. Then, each local controller calculates the receding horizon action, using the optimal price vector received from the coordinator. Nested approaches may exhibit convergence problems, especially when the open-loop dynamic system is unstable. 10 In this paper, a novel price-driven CDMPC algorithm is proposed to extend the idea of Cheng et al. 5 ,6 , Marcos et al. 7 , and Hassanzadeh et al. 11 . Correspondingly, two scenarios are considered based on the presence of inequality constraints in the local MPC problems. The main contributions of this paper are as follows: • In the first scenario, all the constraints are in the form of equalities. Accordingly, a closed-form solution is derived for the price vector and the local variables, which is the global optimal solution to the plant-wide problem. The analytical closed-form solution eliminates the need for iterative procedures 5,7 in every coordination cycle. Furthermore, using this algorithm reduces the coordination cycles to one cycle; and, unlike the method proposed by Marcos et al. 7 , the obtained solution is globally optimal to the corresponding CDMPC problem. • In the second scenario, a convex set of inequality constraints is allowed in each local MPC formulation. In this regard, the lower-level problem remains convex and several strategies 12–17 can be selected to solve the resulting BLP problem. Here, an iterative approach is proposed based on finding a feasible descent direction. 12 The methodology is an alternative to the nested primal-dual programming CDMPC approach used 3



Page 4 of 44

by Cheng et al. 5 , Marcos et al. 7 , and Mohseni 10 , which guarantees convergence to the optimal plant-wide solution without relying on any Newton-based 5–7 coordination systems. • The CDMPC methodologies provided in this paper intend to resolve the convergence issues with the nested CDMPC approaches and open-loop unstable systems reported by Mohseni 10 . In other words, the previously proposed CDMC approaches by Cheng et al. 5 , Marcos et al. 7 , and Mohseni 10 are unable to handle such unstable dynamical systems. Here, an on-line adaptive horizon selection scheme is embedded to the proposed algorithm to ensure closed-loop stability of the entire network. This feature does not change the classic formulation of the price-driven CDMPC, and does not impose any requirements for the sub-systems to be open-loop stable. In the next section, an analytic algorithm is proposed to derive a closed-form solution to the CDMPC problem when only equality constraints are present in the local MPC problems. Next, the generic CDMPC problem is considered when inequality constraints are also present and a feasible direction method is proposed to solve the resulting MPC problem. The control performance of the proposed CDMPC algorithms is illustrated via the applications to a twoCSTR process, as well as an alkylation benchmark in the simulation section of this paper.

2

Preliminaries

In this work, the entire plant is considered to be composed of m interconnected sub-systems with the following state-space representation:

x(k + 1) = Ψx(k) + Γu(k)

4


(1)

Page 5 of 44 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


where: x = [xT1 · · · xTi · · · xTm ]T ∈ Rnx denotes the vector of state variables, and u = [uT1 · · · uTi · · · uTm ]T ∈ Rnu denotes the vector of manipulated input variables. Here, Ψ and Γ are state and input coefficient matrices, respectively; these matrices can be defined by aggregating sub-systems i = 1, · · · , m as: 

 Ψ11 · · · Ψ1j · · · Ψ1m    ..  . . . .  .  . .      Ψ =  Ψi1 · · · Ψii · · · Ψim  ,   ..   .. .. . .   .   Ψm1 · · · Ψmi · · · Ψmm



 Γ11 · · · Γ1j · · · Γ1m    ..  . . . .  .  . .      Γ =  Γi1 · · · Γii · · · Γim     . .  . .. ..  .  .    Γm1 · · · Γmi · · · Γmm

(2)

In this model, the sub-systems are interconnected through states and manipulated input variables, (i.e., the manipulated variable of one sub-system might affect other sub-systems as they explicitly appear within the models). In addition, it is assumed that xi and ui belong to a local set of convex constraints Ci , defined as follows: Ci , {(xi , ui ) ∈ Rnxi × Rnui |gi (xi , ui ) ≤ 0, i = 1, . . . , m}

(3)

where: gi is a set of linear inequalities over xi and ui for sub-system ‘i’. Then, the finite-time centralized MPC formulation, at time instant tk , can be defined as:

min JC = X,U

1 (X(k) − Xset )T Q(X(k) − Xset ) + U (k)T RU (k) 2

s.t. xˆ(k + l + 1|k) = Ψˆ x(k + l|k) + Γˆ u(k + l|k), l = 0, · · · , N − 1 (ˆ xi , uî ) ∈ Ci

(4a) (4b) (4c)

where: N is the prediction horizon, also xˆ and uˆ are states and manipulated input variables inside the controller, respectively. In this paper, for the sake of simplicity, it is assumed

5



Page 6 of 44

the control horizon and the prediction horizon of the system are the same. Although it might lead to a higher computation load, longer control horizons can improve the performance and compensate for package dropouts in networked control systems. 18 Additionally: X(k) = [X1 (k)T , · · · , Xm (k)T ]T is the vector of the predicted state trajectory; Xi (k) = [ˆ xi (k +1|k)T , · · · , xî (k +N |k)T ]T is the vector of the predicted state trajectory for sub-system ‘i’; U (k) = [U1 (k)T , · · · , Um (k)T ]T is the vector of the calculated manipulated variable moves; Ui (k) = [ˆ ui (k|k)T , · · · , uî (k + N − 1|k)T ]T is the vector of the calculated manipulated variable moves for sub-system ‘i’; Q is a positive-definite block-diagonal weighting matrix for the states (i.e., Q = diag{Qi }); and R is a positive-definite block-diagonal weighting matrix for the manipulated variables of the overall system (i.e., R = diag{Ri }); for i = 1, . . . , m. The states of the m sub-systems are assumed to be sampled synchronously at time instants tk = k∆t (with k = 0, 1, . . .). In the remainder of this paper, k is used to denote tk in the discrete-time model. For sub-system ‘i’, i = 1, . . . , m, the prediction model used in the formulation of local MPC has the following form:

xî (k + l + 1|k) = Ψii xî (k + l|k) + Γii uî (k + l|k) +

m X

(1 − β)Ψij xˆj (k|k)

j6=i

+

m X

βΨij xˆj (k + l|k) + Γij uˆj (k + l|k) ,

j6=i

l = 0, · · · , N − 1

(5a)

xî (k|k) = xi (k)

(5b)

with:

β=

   0 l = 0

(5c)

  1 l = 1, · · · , N − 1 where: Ψii and Γii are matrices corresponding to sub-system ‘i’. In (5a), the term ‘Ψii xî (k + P l|k) + Γii uî (k + l|k)’ is the local information, and the term ’ m ˆj (k|k)’ is the j6=i (1 − β)Ψij x 6


Page 7 of 44 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


measured value of the state variables of other sub-systems ‘j 6= i’. This measurement is transmitted to an upper-level coordinator from the corresponding local MPCs and then shared with the local MPC ‘i’ at the beginning of the prediction horizon ‘l = 0’. Accordingly, the P term ’ m ˆj (k + l|k) + Γij uˆj (k + l|k) ’ characterizes the unknown interaction inforj6=i βΨij x mation between the other sub-systems and sub-system ‘i’. Accordingly, equation (5a) can be modified to include an additional variable vî which represents the interaction information from other subsystems as:

xî (k + l + 1|k) = Ψii xî (k + l|k) + Γii uî (k + l|k) +

m X

(1 − β)Ψij xˆj (k|k) + vî (k + l|k) (6)

j6=i

In the proposed control scheme, vî is a decision variable in each local CDMPC controller. Sections 3 and 4 show that the goal of CDMPC is to ensure that vî predicted by subsystem ‘i’ is equal to the actual effect of the interactions whose contributions are given by xˆj and uˆj . Correspondingly, the local interaction error ‘ei ’ is defined as follows:

ei (k + l|k) , vî (k + l|k) −

m X

βΨij xˆj (k + l|k) + Γij uˆj (k + l|k)

(7)

j6=i

where: ei (k + l|k) denotes the difference between the calculated value for vî and the actual other sub-systems with sub-system ‘i’ as captured by the plant-wide model. A specific objective of a coordinated DMPC system is to ensure that the optimal value for interaction term vî renders ei (k + l|k) = 0. In other words, as ei (k + l|k) → 0, the CDMPC approaches the performance of the corresponding centralized MPC. Thus, the overall interaction error over the prediction horizon can be described as follows:

E(k|k) =

m X

Ei (k|k)

i=1

7


(8a)


Page 8 of 44

where, based on (7), E(k|k) can be rewritten in the following form:   Xi (k) m X    E(k|k) = Θi (k)   Ui (k)    i=1 Vi (k)

(8b)

In this formulation, Vi (k) , [ˆ vi (k|k)T , · · · , vî (k+N −1|k)T ]T is the vector of predicted linking variables for sub-system ‘i’, and Θi (k) is the coefficient matrix for the linking constraints, defined as: T

T T T Θi (k) = θ1,i , · · · , θm,i , · · · , θi,i

where: θj,i =

    0N nx

, 0N nxi ×N nui , I i ×N nxi

   θΨj,i , θΓj,i , 0N nx

j ×N nxi

(9)

for j = i (10) for j 6= i

with ‘0’ denoting zero matrices of appropriate dimensions; I being an identity matrix of size N nxi × N nxi ; as well as, θΨj,i and θΓj,i being N nxj × N nxi and N nxj × N nui matrices, respectively: 

θΨj,i

0 ···  0  −Ψji 0 ···  = .. ..  0 . .   0 0 −Ψji

 0  0  , ..   .  0



θΓj,i



0 ··· 0  −Γji  ..  .  0 −Γji . . .    = .  .. ..  .. . . 0      0 ··· 0 −Γji

(11)

Matrix Θi can be written in terms of an augmented matrix of coefficients for Xi , Ui , and Vi as follows:

8


Page 9 of 44 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




θΨ1,i   ..  .   Θi =  0N nxi ×N nxi  ..  .   θΨm,i

θΓ1,i .. .



0 .. .

     IN nxi ×N nxi   = [ΘXi , ΘUi , ΘVi ]  ..  .   0

0N nxi ×N nui .. . GΓm,i

(12)

Therefore, the centralized MPC problem (4) can be reformulated as the following overall plant-wide MPC formulation:

min JP =

X,U,V

m X

JPi

(13a)

i=1

s.t. xî (k + l + 1|k) = Ψii xî (k + l|k) + Γii uî (k + l|k) +

X

(1 − β)Ψij xˆj (k|k) + vî (k + l|k),

j6=i

(ˆ xi , uî ) ∈ Ci   Xi (k) m X   =0 Θi  U (k) i     i=1 Vi (k)

(13b)

(13c)

where: JPi =

1 (Xi (k) − Xi,set )T Qi (Xi (k) − Xi,set ) + Ui (k)T Ri Ui (k) 2

(13d)

Optimization problem (13) provides the basis for the formulation of the sub-system MPCs and the coordinator. In price-driven coordinated distributed MPC, the idea is to penalize constraint (13c), which represents the interactions in the plant-wide problem (13), via a price vector into the objective function JPi . Specifically, the generic price-driven CDMPC problem is a primal-dual problem which includes two levels: the lower-level problem which consists of the distributed MPCs, and the upper-level problem which consists of a coordinator. The

9



Page 10 of 44

bi-level CDMPC problem can be written as the following bi-level optimization problem:

max JD (p, X ∗ (k), U ∗ (k), V ∗ (k)) p   ∗ X (k) m X    U ∗ (k)  = arg min{JD = J¯Di }     i=1 V ∗ (k)   Xi (k)    = bi s.t. Ai  U (k) i     Vi (k) (ˆ xi , uî ) ∈ Ci

(14a)

(14b)

(14c)

(14d)

where:

JDi

  Xi (k)   1  (Xi (k) − Xi,set )T Qi (Xi (k) − Xi,set ) + Ui (k)T Ri Ui (k) + pT Θi  = U (k)  i  2   Vi (k)

Ai = [AXi , AUi , AVi ],  0  Inxi  −Ψii In xi  A Xi =  .  0 ..   0 0

··· ··· .. .



0   0   , AUi ..  .   

−Ψii Inxi

(14e)

(14f)   0 ··· 0  −Γii  ..  ...  0 −Γ .  ii   =  .  , AVi = − IN nxi ×N nxi . .  .. .. .. 0      0 ··· 0 −Γii (14g)

  P  j6=i (1 − β)Ψij xˆj (k)     0   bi =   .   . .     0

(14h)

10


Page 11 of 44 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


for i = 1, · · · , m. In the optimization problem (14), ‘p’ is a price vector calculated by the coordinator.

3

Coordination Algorithm Scenario I: Equality Constraints

In this section, an analytic CDMPC algorithm is presented for linear systems, when no inequality constraint is present in the overall formulation (14). Specifically, the CDMPC problem is formulated as follows:

max JD (p, X ∗ (k), U ∗ (k), V ∗ (k)) p   ∗ X (k) m X    U ∗ (k)  = arg min{JD = J¯Di }     i=1 V ∗ (k)   Xi (k)    s.t. Ai  U (k)  i  = bi   Vi (k)

(15a)

(15b)

(15c)

Remark 3.1. The closed-form optimal solution of the lower-level problem (14) for the current value of the price vector ‘p’ can be derived as:

       T −1 QXset  Xset  −Q ΘX + AX  X (k; p)  − [ΘX + AX , ΘU + AU ]T p =  −   = Υ−1  T  p (16a) 0 0 −R−1 ΘU + AU U ∗ (k; p) 

∗

  ∗ X (k; p) V ∗ (k; p) = [AX , AU ]  −b U ∗ (k; p)

(16b)

Proof. See Appendix S.1 Remark 3.2. For the bi-level optimization problem (15), the closed-form solution for the 11



Page 12 of 44

optimal price vector can be derived as: −1 p∗ = [θX , θU ]Υ−1 [θX , θU ]T θX Xset − b

(17a)

θX = ΘX + AX

(17b)

θU = ΘU + AU

(17c)

where:

Proof. See Appendix S.2. Theorem 3.1. The closed-form solution of the optimal price vector stated in (17a) p∗ is not singular. Proof. See Appendix S.3. Thus, based on Remarks 3.1 and 3.2, the closed-form solution to the CDMPC problem (15) can be derived as:     T −1 −1 Xset  −Q ΘX + AX  −1 T Z ∗ (k; p) =  [θ , θ ]Υ [θ , θ ] θ X − b (18a) +  X U X U X set T 0 −R−1 ΘU + AU and:     ! T −1 −1 Xset  −Q ΘX + AX  ∗ −1 T V (k; p) = [AX , AU ]  θX Xset − b −b + T  × [θX , θU ]Υ [θX , θU ] 0 −R−1 ΘU + AU

(18b)   X  where: Z =  . U Theorem 3.2. The closed-form solution set (18) for ‘Z’ and ‘V ’ is globally optimal to the CDMPC problem (15). 12


Page 13 of 44 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


Proof. See Appendix S.4. The next remark summarizes the necessity of developing this algorithm w.r.t the classical price-driven CDMPC approaches. 5–7 The analytic CDMPC approach is summarized in Algorithm 1. Remark 3.3. In classical primal-dual CDMPC approaches 5–7 Newton-like methods are used. Given sensitivity information of the local controllers, i.e.

dZ dpT

, the coordinator uses an it-

erative gradient-based procedure such as Newton’s method to update the price vector. Thus, the coordinator and local controllers communicate until the norm of the overall interaction constraint (13c) is less than a predefined upper limit and the optimal plant-wide solution is achieved. This always results in a locally optimal solution to the CDMPC problem (15) and requires at least two coordination cycles and two iterations per cycle to converge. 7 In the proposed approach, the KKT condition is solved analytically to obtain the optimal price vector, and the closed-form solutions for Z ∗ , V ∗ , and p∗ reduces the coordination cycles to one and one iteration per cycle to obtain the globally optimal solution. Algorithm 1: The analytic CDMPC algorithm Initialization: Number of sub-systems m, initial state x(k), positive definite matrix Q, positive definite matrix R, initial prediction horizon N ; Local MPCs: Given x(k) and N , formulate local KKT conditions from (15); Local MPCs: Calculate θXi and θUi ; Local MPCs: Send θXi , θUi , Υi , Xspi , bi to Coordinator; Coordinator: Calculate the optimal price vector from (17a); Coordinator: Send the optimal price vector and N to the local controllers; Local MPCs: Calculate Zi∗ (k) and Vi∗ (k) using (18); Local MPCs: Apply the RHC action to the plant; Output: Z ∗ , V ∗ , and p∗ ;

13



4

Page 14 of 44

Coordination Algorithm Scenario II: Inequality Constraints

In this section, a feasible descent direction approach is presented to coordinate distributed MPCs subject to convex sets of equality and inequality constraints. In the proposed formulation, sub-system MPCs and the coordinator communicate and exchange information iteratively every sampling time to find the optimal price and state trajectories. Feasible direction methods often require first-order derivate information of the optimization problem, which defines the main content of communication between the coordinator and local controllers. In the proposed hybrid solution strategy, the algorithm is initiated with a relaxed version of KKT conditions of the lower-level problem and subsequently solves a single level optimization. Once the initial active set of constraints is identified, a feasible direction for the price vector is iteratively calculated through coordination cycles between local MPCs and the coordinator. Then, the optimal price is sent to local MPCs, and each local controller calculates the receding horizon action based on the fixed optimal value of price.

4Appendix S.1

Algorithm Description

In the rest of this paper, for the sake of simplicity, the CDMPC problem (14) is rewritten as the following BLP:

min F (p, Z¯ ∗ )

(19a)

p

¯ Z¯ ∗ = arg min f (p, Z)

(19b)

s.t. Hi,j (Z¯i ) = 0,

j ∈ JHi

(19c)

Gi,j (Z¯i ) ≤ 0,

j ∈ JGi

(19d)

14


Page 15 of 44 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


where:   Xi (k)    Z¯i =   Ui (k)  ,   Vi (k)

(19e)

F (p, Z¯ ∗ ) = − JD (p, Z¯ ∗ ), ¯ = f (p, Z¯i ) = JD (p, Z)

m X

(19f) JDi (p, Z¯i ),

(19g)

i=1

Hi,j (Z¯i ) = Ai Z¯i − bi ,

(19h)

Gi,j (Z¯i ) = Aineq Z¯i − bineq , i i

(19i)

ineq = [Aineq Aineq i Xi , AUi , 0]

(19j)

and JHi , JGi are finite sets of indexes corresponding to equality and inequality constraints in sub-system ‘i ’, respectively. In the BLP problem (19), the lower-level problem represents the distributed network of MPC controllers and the upper-level problem corresponds to the coordinator. In addition, for a fixed value of the price vector, Si (p) = {(Z¯i ) ∈ R2×nxi +nui |Hi,j (Z¯i ) = 0, Gi,j (Z¯i ) ≤ 0, j ∈ JHi ∪ JGi } is defined as the set of feasible solutions and JGi,act (p) = {j ∈ JGi |Gi,j (Z¯i ) = 0} is defined as the set of indexes of active constraints of sub-system ‘i ’. Assumption 4.1. To guarantee that there exists at least one set of solution to the BLP problem (19), it is assumed that the feasible solution set of the lower-level problem for a fixed value of the price vector (Si (p)) is unique and uniformly compact. In other words, there exists an open set C such that for any value of ‘p’, Si (p) ∈ C is satisfied. (See Appendix S.5 for more details.) Assumption 4.2. The vectors ∇Z¯i Hi,j (Z¯i∗ ) and ∇Z¯i Gi,j (Z¯i∗ ), for j ∈ JHi ∪ JGi,act , are linearly independent. (See Appendix S.6 for more details.) For any fixed value of ‘p’, the Lagrangian of the lower-level problem in (19) can be stated 15



Page 16 of 44

as:

¯ ν, λ) = L(p, Z,

m X

Li (p, Z¯i , νi , λi )

(20a)

i=1

where:

Li (p, Z¯i , νi , λi ) = fi (p, Z¯i ) +

X

νi,j Hi,j (Z¯i ) +

j∈JHi

X

λi,j Gi,j (Z¯i )

(20b)

j∈JGi,act

and (νi , λi ) are the Lagrange multiplier associated with equality and inequality constraints, respectively. Additionally, the following set of assumptions are required to ensure that the lower-level problem in (19) satisfies regularity of constraints at Z¯i∗ ∈ (Si (p)): Remark 4.1. Second order sufficiency condition (SOSC) holds at Z¯i∗ ∈ (Si (p)) for the lowerlevel problem in (19). SOSC holds if for some feasible point Z¯i∗ ∈ (Si (p)), there exist Lagrange multipliers (νi∗ , λ∗i ) such that the KKT conditions for the lower-level problem in (19) are satisfied and for all non-zero feasible directions d that wiT ∇2[Xi ,Ui ,Vi ] Li (p, Xi∗ , Ui∗ , Vi∗ , νi∗ , λ∗i )wi > 0. Here, the lower-level problem in (19) is a convex quadratic programming problem and the weighting matrices Qi and Ri are positive definite. Therefore, ∇2[Xi ,Ui ,Vi ] Li (p, Xi∗ , Ui∗ , Vi∗ , νi∗ , λ∗i ) is positive definite, which verifies this optimality condition. Remark 4.2. Assumptions 4.1 and 4.2 together with Remark 4.1 ensure that Z¯i is a Lipschitz function of the price vector ‘p’. 13

4Appendix S.1.1

The Initial Active-set

In order to determine the initial set of active constraints JGi,act (p) for the BLP problem (19), consider the following optimization:

¯ min F (p, Z)

(21a)

¯i ,νi ,λi p,Z

T T s.t. ∇Z¯i fi (p, Z¯i ) + ∇Z¯i Hi (Z¯i ) νi + ∇Z¯i Gi,act (Z¯i ) λi,act = 0, 16


(21b)

Page 17 of 44 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


Hi,j∈JHi (Z¯i ) = 0,

(21c)

Gi,j∈JGi (Z¯i ) ≤ 0,

(21d)

λi,j∈JGi Gi,j∈JGi (Z¯i ) = 0,

(21e)

λi,j∈JGi ≥ 0

(21f)

where:

Gi,act (Z¯i ) = Gi,j (Z¯i )

,

λi,act = λi,j

j∈JGi,act

(21g) j∈JGi,act

in which the lower-level problem in (19) is replaced by the corresponding KKT conditions to form the single-level optimization problem (21). Note that, based on Assumptions 4.1 and 4.2, the corresponding Lagrange multipliers νi and λi are bounded and can be uniquely determined. The resulting problem (21) can also be regarded as a mathematical program with equilibrium constraints (MPEC); however, this problem is neither differentiable nor regular 19 due to the complementarity constraint present in (21). Then, instead of (21e), it is replaced by the following relaxed equality constraint:

λi,j∈JGi Gi,j∈JGi (Z¯i ) = ηi2

(22)

where: ηi ≥ 0 is defined as a perturbation parameter for sub-system ‘i’. This parameter is decreased at each iteration until it approaches zero. However, the optimization problem remains non-differentiable, and it is required to replace (22) with a smoothed constraint. In this work, a CHKS smoothing 20 function ‘Π’ is used, which is defined as:

Π(a, b, η) = a + b −

p (a − b)2 + 4η 2

(23)

In other words, for any given value of ηi ≥ 0, Π(−Gi,j∈JGi (Z¯i ), λi,j∈JGi , ηi ) = 0 is equivalent 17



Page 18 of 44

to (22). This leads to the following smoothed nonlinear programming problem:

¯ min F (p, Z)

(24a)

¯i ,νi ,λi p,Z

T T s.t. ∇Z¯i fi (p, Z¯i ) + ∇Z¯i Hi (Z¯i ) νi + ∇Z¯i Gi,act (Z¯i ) λi,act = 0,

(24b)

Hi,j∈JHi (Z¯i ) = 0,

(24c)

Gi,j∈JGi (Z¯i ) ≤ 0,

(24d)

Π(−Gi,j∈JGi (Z¯i ), λi,j∈JGi , η) = 0,

(24e)

λi,j∈JGi ≥ 0

(24f)

for a decreasing sequence of η. Equivalently, problem (24) can also be written as:

¯ min − JD (p, Z)  Qi  s.t.  Ri  

(25a)

¯i ,νi ,λi p,Z

  T  Xi,set Qi      Z¯i −  0  + ΘTi p + ATi νi + AineqT λi,act = 0, i,act       0 0 

(25b)

Ai Z¯i − bi = 0,

(25c)

Aineq Z¯i − bineq ≤ 0, i i

(25d)

¯ Π(−Aineq i,j∈JG (Zi ), λi,j∈JGi , η) = 0,

(25e)

λi,j∈JGi ≥ 0

(25f)

i

Therefore, an initial guess for JGi,act of the BLP (19) can be identified by solving (25). The procedure is listed in Algorithm 2.

18


Page 19 of 44 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


Algorithm 2: Initial guess for JGi,act of the BLP problem (19) Input: 0 < ηmin < η 0 ≤ 1, l = 0, β ∈ (0, 1), Zi0 (k), p0 = 0; while (η l ≥ ηmin ) do Coordinator - Local MPCs: Solve the NLP problem (24); η l+1 = βη l ; l = l + 1; Coordinator - Local MPCs: update Z¯il (k) and pl ; end for j ∈ JGi do if Gi,j (pl , Z¯i ) = 0 then Coordinator: JG0i,act = {j ∪ JG0 i,act }; end end Output: JG0i,act , Z¯il (k), and pl ;

4Appendix S.1.2

The Feasible Directions Method

Once the initial stage is finished, an approximate local solution to the original BLP (19) is determined. In the next stage, this solution is improved using a parametric optimization ¯ stays in a certain tolerance. In approach, until the change in the objective function ‘F (p, Z)’ this approach, the lower-level problem in (19) can be thought of as a parametric optimization problem, where the main optimization variable ‘Z¯i ’ is parameterized by the price vector ‘p’. Accordingly, for a perturbation in ‘p’ along direction ‘d’, the impact on ‘Z¯i ’ in the lower-level problem can be measured along direction ‘w’. Following the work of Savard and Gauvin 12 , the steepest descent direction along ‘d’ can be derived by solving the following BLP:

¯ ¯ min ∇p F (p, Z(k))d + ∇Z¯ F (p, Z(k))w

(26a)

s.t. ||d||∞ ≤ 1

(26b)

d

19



min w

dT , w T

  d ¯ ∇2p,Z(k) L(p, Z(k), νLP , λLP )   ¯ w

¯ s.t. ∇p Hi (Z(k))d + ∇Z¯i Hi (Z¯i (k))wi = 0

m X

(26c)

(26d)

¯ ∇p Gi,act (Z(k))d + ∇Z¯i Gi,act (Z¯i (k))wi ≤ 0 ¯ ¯ ∇p f (p, Z(k))d + ∇Z¯ f (p, Z(k))w =

Page 20 of 44

∗ , λ∗LP,i )d ∇p Li (p, Z¯i (k; p), νLP,i

(26e) (26f)

i=1

where: Li (p, Z¯i (k; p), νLP,i , λLP,i ) is the Lagrangian of the lower-level problem in (19). In the optimization problem (26), the coordinator finds the optimal direction of the price vector, then the local controllers optimize over their corresponding feasible direction to minimize the ¯ Equivalently, the main idea behind the BLP problem plant-wide objective function F (p, Z). (26) is to find a feasible set of directions ‘d’ and ‘wi ’ that results in decreasing the upper-level ¯ ¯ objective function of the original BLP problem (19), i.e. ∇p F (p, Z(k))d + ∇Z¯ F (p, Z(k))w < 0. Similar to Section 4Appendix S.1.1, the lower-level problem in (26) can be replaced by its KKT conditions and then smoothed using a CHKS smoothing function to form the following smooth single-level optimization:

min

d,w,νQP ,λQP

¯ ¯ ∇p F (p, Z(k))d + ∇Z¯ F (p, Z(k))w

s.t. ||d||∞ ≤ 1

(27a) (27b)

¯ νLP , λLP )w + ΘT d)+ 2(∇Z¯ L(p, Z,  T ¯  ∇Z¯ H(Z(k))  T ¯ λQP = 0   νQP + Gact (Z(k)) ¯ ∇Z¯ f (p, Z(k)) ¯ ∇p Hi (Z(k))d + ∇Z¯i Hi (Z¯i (k))wi = 0 ¯ ¯ ∇p f (p, Z(k))d + ∇Z¯ f (p, Z(k))w −

m X

(27c)

(27d) ∗ ∇p Li (p, Z¯i (k), νLP,i , λ∗LP,i )d = 0

(27e)

i=1

¯ ∇p Gi,act (Z(k))d + ∇Z¯i Gi,act (Z¯i (k))wi ≤ 0 20


(27f)

Page 21 of 44 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


¯ Π − ∇p Gi,act (Z(k))d − ∇Z¯i Gi,act (Z¯i (k))wi , λQP,i , ηQP,i = 0

(27g)

λQP,i ≥ 0

(27h)

T T where: νQP = [νQP,1 , . . . , νQP,m ]T and λQP = [λTQP,1 , . . . , λTQP,m ]T are Lagrange multipliers

associated with equality and inequality constraints of the lower-level problem in (26), respectively. Once the optimal directions ‘ds ’ and ‘wis ’ are found at iteration ‘s’, the following feasibility problem is performed:

min α

(28a)

s.t. α > 0

(28b)

α

¯ ps + αds )) < F (ps , Z(k; ¯ ps )) F (ps + αds , Z(k;

(28c)

Z¯i (k; ps + αds ) ∈ Si (ps + αd∗ )

(28d)

to find the minimum positive step-size along direction ‘d’ which ensures feasibility of ‘Z¯is ’. Then, the price vector ‘p’ and local optimization variables ‘Z¯i ’ are updated along the optimal directions ‘ds ’ and ‘wis ’, respectively. In other words, ‘α∗s ’ is found such that: ¯ ps + α∗s ds )) < F (ps , Z(k; ¯ ps )) while the solution ‘Z¯i (k; ps + α∗s ds )’ remains F (ps + α∗s ds , Z(k; feasible. 12 Therefore, the price vector and the local optimization variables are updated, based on the optimal value of the step-size ‘αs ’:

ps+1 = ps + α∗s ds

(29)

¯ ps+1 ) Z¯i = Z(k;

(30)

Accordingly, the procedure to find a feasible descent direction for the BLP problem (19) is listed in Algorithm 3

21



Algorithm 3: Feasible descent direction algorithm to solve the BLP (19) Input: 0 < ηmin < η 0 ≤ 1, l = 0, s = 0, β ∈ (0, 1), Z¯il (k), pl , JG0i,act ; ¯ < 0) do ¯ while (∇p F (p, Z(p))d + ∇Z¯ F (p, Z)w while (η l ≥ ηmin ) do Coordinator - Local MPCs: Solve the NLP problem (27) ; η l+1 = βη l ; l = l + 1; Coordinator - Local MPCs: update dl and wl ; end s = s + 1; Set ds = dl and ws = wl ; Coordinator - Local MPCs: Solve problem (28); Coordinator: Update ps+1 = ps + αs ds ; ¯ ps + α s ds ) ; Local MPCs: Update Z¯is+1 = Z(k; for j ∈ JGi do if Gi,j (ps , Z¯is ) = 0 then Coordinator: JGs i,act = {j ∪ JGs i,act }; end end Set l = 0; end Output: p∗ , Z¯i∗ (k);

4Appendix S.2

Convergence Analysis

In this section, convergence property of the proposed CDMPC algorithm is studied. First, the initialization stage is considered, in which the initial set of active constraints JGi,act (p) are determined using problem (24), for a decreasing sequence of ηi . Theorem 4.1. For every solution derived from (24), there exists ηmax > 0 and a compact ¯ max ) such that: S(p, ¯ η) ⊆ C(η ¯ max ) holds for all η ∈ (0, ηmax ]. set C(η Proof. See Appendix S.7.

22


Page 22 of 44

Page 23 of 44 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


¯ η(k)) is non-empty and compact Remark 4.3. Based on Theorem 4.1, the solution set S(p, for all values of 0 < η ≤ ηmax = 1. 20 Therefore, for a decreasing sequence of η, a continuous path towards the solution of (21) can be generated. 20,21 Once the initial set of active constraints is identified, the BLP problem (26) along with step-size calculation (28) are solved iteratively to find a local solution to the original BLP problem (19). In the following theorem, it is proved that the CDMPC algorithm, i.e. Algorithm 3, is globally convergent to a locally optimal solution of the plant-wide problem (19). Theorem 4.2. If assumptions 4.1 and 4.2 hold, then: (i) For the decreasing sequence η → 0, the initialization stage (Algorithm 2) finds an initial set of active constraints for the BLP problem (19). (ii) If (p∗ , Z¯ ∗ ) denotes an optimal solution to the BLP problem (19), then:

¯ p∗ ))d + ∇Z¯ F (p∗ , Z(k, ¯ p∗ ))w(p∗ , d) ≥ 0 ∇p F (p∗ , Z(k,

(31)

(iii) The proposed CDMPC approach (Algorithm 3) is globally convergent to an optimal solution of (19). Proof. See Appendix S.8.

5

Closed-loop Stability

In this work, the idea proposed by Grüne and Pannek 22,23 is used to implement a stability criterion. The main idea is to ensure closed-loop stability without imposing any terminal constraints or cost functions into the local MPC problems. This criterion is applicable to MPC schemes subject to admissible sets of states and input constraints. 18

23



5Appendix S.1

Page 24 of 44

Stability of the CDMPC scheme

If X ∗ (k), U ∗ (k), and V ∗ (k) denote the optimal open-loop solution of the proposed CDMPC algorithms; and, u∗i,RHC (N, xˆ∗i (n)) denotes the receding horizon control (RHC) action calculated by solving either (14) or (19) with a prediction horizon length of N and an initial value of xî (n), a value function VND can be defined as:

VND

−1 X N D ∗ ∗ xˆ(n) = l xˆ (n), uRHC (N − n, xˆ(n))

(32a)

n=0

where: lD is the stage cost of the overall system, defined as:

l

D

xˆ

∗

(n), u∗RHC (N

m T X xˆ∗i (n) − xî,set (k) Qi (n) xˆ∗i (n) − xî,set (k) − n, xˆ(n))) = i=1

T + u∗i,RHC (N − n, xî (n)) Ri (n) u∗i,RHC (N − n, xî (n))

(32b)

and, the matrices Qi and Ri are defined over the prediction horizon as: 



  Qi (1)  Ri (1)      . .  , Ri =   .. .. Qi =          Qi (N ) Ri (N )

(32c)

The difference between the open-loop control and receding horizon control (RHC) trajectories for a dynamic programming problem at sampling time k and the initial prediction horizon N = 5 is depicted in Figure 1. In this example, the system performs five open-loop optimization problems. Once the open-loop trajectory for N = 5 is found, the next openloop optimization is performed with N = 4 from the initial move previously computed. This process continues until five RHC points are obtained and the value function VND is calculated.

In the following theorem, it is shown that under certain conditions VND can be a Lyapunov function for the closed-loop system (1) under the control of the proposed CDMPC schemes, 24


Page 25 of 44 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


N 5

x i 1  u i , RHC

N

, x i 0 

x i 0 

u i , RHC  N  1, x i 1 

x i 2  u i , RHC  N  2, x i 2 

x i 3  u i , RHC  N  3, x i 3 

x i 4  u i , RHC  N  4, x i 4 

Receding Horizon Control

x i 5 

Open-loop Control

Figure 1: Receding horizon control (RHC) vs open-loop control trajectories for sub-system i with a typical prediction horizon length of N = 5.

i.e. Algorithm 1 or Algorithm 3. Theorem 5.1. If the value function VND is calculated from the result of the proposed CDMPC Algorithm 1 or Algorithm 3, using a posteriori approach, 18 and the following inequality holds: VND (ˆ x(k)) ≥ VND (ˆ x(k + 1)) + lD xˆ(k), u∗RHC (N, xˆ(k))

(33)

then the closed-loop system (1) is stable around the origin. Proof. See Appendix S.9. Based on Theorem 5.1, an adaptive horizon selection procedure is implemented to ensure the closed-loop system (1) is stable. Since no terminal costs or constraints are added to the local MPC controllers, the formulation of the local MPCs remains unchanged. Furthermore, it is assumed that all sub-system controllers have the same length of prediction horizon ‘N ’. This procedure is explained in detail in Algorithm 4.

25



Algorithm 4: A posteriori adaptive horizon algorithm for the CDMPC schemes Initialization: n = 0, prediction horizon N , number of sub-systems m, initial state xˆ(n) = x(k), Qi , Ri ; Coordinator - Local MPCs: Calculate VND (ˆ x(n)) with N, m, xˆ∗ (n), Qi , Ri ; Local MPCs: Apply u∗RHC (N, xˆ(n)) to the internal model; Local MPCs: Calculate xˆ∗ (n + 1); Coordinator - Local MPCs: Calculate VND (ˆ x(n + 1)) with N, m, xˆ∗ (n + 1), Qi , Ri ; if VND (ˆ x(n)) ≤ VND (ˆ x(n + 1)) then Coordinator - Local MPCs: Perform Horizon Prolongation; else if

x(n))−VND (ˆ x(n+1)) VND (ˆ x(n)) lD x ˆ∗ (n),u∗RHC (N,ˆ

≥ 1 then

Coordinator - Local MPCs: Perform Horizon Shortening; else Coordinator - Local MPCs: Perform Horizon Prolongation; end end Output: X ∗ based on the accepted value of N ;

To ensure stability of the closed-loop system (1), at each sampling instant, ‘N ’ is iteratively increased until (33) is satisfied. This value for the horizon is used for the next sampling period as an initial guess. Thus, to ensure that iteratively increasing ‘N ’ will eventually lead to a horizon for which (33) holds, the following assumption is made. ¯’ Assumption 5.1. It is assumed that for all xˆ(k), there exists a finite horizon length ‘N ¯. such that (33) holds for N ≥ N Assumption 5.1 can be seen as a performance assumption which requires the existence of ¯ . If no such horizon exists, no prolongation strategy based on Theorem a horizon length N 5.1 can be designed which guarantees the closed-loop stability of system (1). The following remark summarizes the adaptive strategy to obtain a proper horizon length ‘N ’ at each sampling time ‘k’. Remark 5.1. The adaptive horizon scheme, described in Algorithm 4, is initiated with N = 2 at k = 0. Then, it finds the lowest possible value for ‘N ’ that satisfies (33) through a 26


Page 26 of 44

Page 27 of 44 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


horizon prolongation procedure. Next, at each sampling time k > 0, ‘N ’ is initiated from the previously calculated value at sampling time k − 1. Throughout this procedure, it could be required to shorten the length of the horizon as long as the inequality (33) is satisfied for a smaller value of ‘N ≥ 2’.

5Appendix S.2

The Overall CDMPC Algorithm

In this section, a summary of the CDMPC algorithms applied to the plant-wide problem is provided. The hierarchical dual decomposition of the hypothetical centralized controller defined in (13) is related to the modification applied to the corresponding existing network of local decentralized MPC controllers. The modification can be interpreted as penalizing local violations of the overall interaction equality constraints via the price vector into the objective function of local MPCs. The local MPC controllers formulate their own optimization problems and send the required local information to the coordination level. The coordinator decides whether to apply the analytic CDMPC or feasible descent direction CDMPC based on the existence of inequality constraints in sub-problem formulations. Then, based on the chosen CDMPC method, finds the minimum required prediction horizon to ensure the stability of the closed-loop system. Once the plant-wide optimum solution, i.e. Z ∗ (k), is found, the receding horizon control action is applied to the plant. The overall process is presented in Algorithm 5.

27



Algorithm 5: The overall CDMPC algorithm Initialization: Number of sub-systems m, initial state x(k), positive definite matrix Q, positive definite matrix R, initial prediction horizon N0 , 0 < ηmin < η 0 ≤ 1, and β ∈ (0, 1); if inequality constraints present in local MPC problems then Local MPCs: Given x(k) and N0 , formulate local KKT conditions of problem (19); else Local MPCs: Given x(k) and N0 , formulate local KKT conditions of problem (14); end Coordinator - Local MPCs: Perform Algorithm 4 with x(k), Q, R, N0 , ηmin , η 0 , β ; Coordinator: Send the optimal price vector and N to the local controllers; Local MPCs: Calculate the optimum local vector of variables Zi∗ (tk ) using p∗ (tk ) and N ; Local MPCs: Apply the RHC action to the plant; Local MPCs: Update N0 = N ;

6

Simulation Case Studies

In this section, two benchmarks are studied to illustrate the CDMPC algorithms discussed in this study. The first case study is a process composed of two interconnected CSTRs, 24 which represents an open-loop unstable system subject to bounds on the control actions. The third case study is an alkylation process 25 which is an open-loop stable system without any bounds defined over the decision variables. The optimization problems are formulated inside MATLAB using the YALMIP 26 interface. Simulations were performed on an Intel Core-i7 processor with 8 GB of memory. The performance of the proposed algorithms, in this paper, were compared to the corresponding centralized MPC and decentralized MPC systems. As mentioned above, the decentralized MPC scheme in this work ignores all the interactions among sub-systems; and,

28


Page 28 of 44

Page 29 of 44 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


the corresponding optimization problem is defined as:

min JDC =

X,U,V

m X

JDCi

(34a)

i=1

s.t. xî (k + l + 1|k) = Ψii xî (k + l|k) + Γii uî (k + l|k), (ˆ xi , uî ) ∈Ci

(34b)

where: JDCi =

1 (Xi (k) − Xi,set )T Qi (Xi (k) − Xi,set ) + Ui (k)T Ri Ui (k) 2

(34c)

for l = 0, . . . , N − 1, with JDC being the overall objective function of the decentralized MPC network. Remark 6.1. In the following case studies, the adaptive horizon algorithms for the decentralized MPC (34) and the centralized (plant-wide) MPC (13) formulations are implemented based on Algorithm 6. Algorithm 6: A posteriori adaptive horizon algorithm for the decentralized MPC Initialization: n = 0, N , m, xˆ(n) = x(k), Qi , Ri ; Calculate VND (ˆ x(n)) with N, m, xˆ∗ (n), Qi , Ri ; Apply u∗RHC (N, xˆ(n)) to the internal model; Calculate xˆ∗ (n + 1); x(n + 1)) with N, m, xˆ∗ (n + 1), Qi , Ri ; Calculate VND (ˆ x(n + 1)) then x(n)) ≤ VND (ˆ if VND (ˆ Perform Horizon Prolongation; else if

VND (ˆ x(n))−VND (ˆ x(n+1)) lD

x(n)) x ˆ∗ (n),u∗RHC (N,ˆ

≥ 1 then

Perform Horizon Shortening; else Perform Horizon Prolongation; end end Output: X ∗ based on the accepted value of N ;

29



Page 30 of 44

Remark 6.2. In this work, the decentralized MPC (34) and the centralized (plant-wide) MPC (13) performance benchmarks were solved using the Active-set algorithm for convex quadratic programming problems, available in the MATLAB’s Optimization toolbox. Remark 6.3. To quantify the difference between the performance improvement of the proposed CDMPC schemes vs the corresponding decentralized MPC, the mean absolute percentage error (MAPE) is used in this paper. The MAPE between the performance curve of the CDMPC and the decentralized MPC is calculated, over the simulation life-time, via the following formula :

M AP E = 100 × mean abs

6Appendix S.1

DC − JD JD

J

! (35)

Two-CSTR Process

In this section, the feasible direction scheme is implemented on a two reactive continuous stirred tank reactor (CSTR) benchmark, taken from. 24 The system is comprised of two non-isothermal CSTRs with interconnections and a recycle stream as depicted in Figure 2; and, the following set of exothermic reactions take place in these reactors with substances k

k

k

1 2 3 A and B: (i) A − → B, (ii) A − → U P , and (iii) A − → DP . Here, UP and DP stand for the

undesired product and the desired product, respectively. As in Figure 2, CSTR I has two feed streams (one has fresh stream of substance ‘A’ with molar concentration CA0 , flow-rate F0 and temperature T0 , and the other stream is from the output of CSRT II containing the recycle stream of unreacted substance A at flow-rate Fr , molar concentration CA2 and temperature T2 ); additionally, CSTR II is fed from the output of CSTR I as well as another fresh stream of substance A at flow-rate F3 , molar concentration CA03 , and temperature T03 . The two CSTRs are equipped with jackets to remove/provide heat, due to non-isothermal nature of reactions. The main focus of this simulation is on stabilizing the system around the open-loop unstable steady-state operating point, to avoid high temperatures, while simultaneously 30


Page 31 of 44 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


achieving reasonable conversion. The manipulated variables of the system are the heat input rates Q1 and Q2 , and the inlet concentrations CA0 and CA03 . The continuous-time plant model of this process can be derived based on linearization around the unstable steady-state operating condition, i.e.

s s (T1s , CA1 , T2s , CA2 ) = (457.9[K], 1.77[kmol/m3 ], 415.5[K], 1.75[kmol/m3 ])

(36a)

s s ) = (0[kJ/hr], 4[kmol/m3 ], 0[kJ/hr], 2[kmol/m3 ]) , Qs2 , CA03 (Qs1 , CA0

(36b)

as the following:   0 0 0 0  0.0043  25.2 1284.3 35      0  −0.3 −45.9  4.998 0 0 0 35     x + x˙ =   u    13.3  0 0.0014 0  0 −2.8 336.2   0       0 0 0 10 0 13.3 −0.1 −24.9 



(37)

The ODE system (37) can be decomposed into two subsystems: CSTR I and CSTR II, which is depicted in Figure 2. Accordingly, the two sub-systems are defined as as:

xT = ([x1,1 , x1,2 ]|[x2,1 , x2,2 ])T = ([T1 , CA1 ]|[T2 , CA2 ])T

(38)

uT = ([u1,1 , u1,2 ]|[u2,1 , u2,2 ])T = ([Q1 , CA0 ]|[Q2 , CA03 ])T

(39)

in addition, the corresponding discretized-time model can be derived as:  1.1357   0.0014  xˆ(k + 1) =   0.0704  0.0001

6.1467 0.7954 0.2556 0.0557

0.1852 0.6696





0.0000     0.0000 0.0002 0.1466    xˆ(k) +    0.9917 1.5870 0.0000   0.0005 0.8875 0.0000

0.0778 0.0223 0.0022 0.0007

0.0000 0.0113



  0.0000 0.0039   uˆ(k)  0.0000 0.0404  0.0000 0.0471 (40)

31



F0, T0, CA0

Page 32 of 44

Fr, T2, CA2 F1, T1, CA1

Coolant in

Coolant out CSTR I

Composition analyzer

F3, T03, CA03

F2, T2, CA2

Coolant in

Coolant out

Temperature sensor

CSTR II MPC I Composition analyzer

Temperature sensor

MPC II

Coordinator

Figure 2: Schematic of the Two CSTR Case Study

for a sampling period of Ts = 0.005 [s]. The CDMPC trajectory is calculated according to the approach discussed in Section 3, based on the adaptive horizon scheme explained in Section 5Appendix S.1. The summary of controller parameters and the variable bounds of these two sub-systems are listed in Table 1. Table 1: The Two-CSTR variable bounds and controller parameters Initial Conditions: Weighting Matrices: Upper bounds: Lower Bounds: Sampling Period

MPC I x1 (0) = 462.9[K], 0.27[kmol/m3 ] Q1 = 5I, R1 = 15I lbu1 = 5[kJ/hr], 8[kmol/m3 ] ubu1 = − 5[kJ/hr], 0[kmol/m3 ] 0.005[s]

MPC II x2 (0) = 410.5[K], 3.45[kmol/m3 ] Q2 = 5I, R2 = 15I lbu2 = 5[kJ/hr], 4[kmol/m3 ] ubu2 = − 5[kJ/hr], 0[kmol/m3 ] 0.005[s]

The comparison between the required length of predictions horizons (N ) between the decentralized MPC, the centralized MPC and the CDMPC trajectories is depicted in Figure 3. As expected from the CDMPC algorithm, this networked control system is acting very similar to an equivalent monolithic centralized MPC in terms of their stability criteria. The state trajectories and manipulated input variables of the Two CSTR system are depicted in Figures 4 and 5, respectively. As shown, the CDMPC trajectory follows the centralized 32


Page 33 of 44

20 Decentralized MPC Centralized MPC CDMPC

Prediction Horizon (N)

15

10

5

0

0

0.05

0.1

0.15

0.2

0.25

time[s]

Figure 3: Comparison of number predictions horizons needed (N) between the Decentralized MPC, the Centralized MPC and the CDMPC trajectories vs simulation time (Two-CSTR Process)

MPC trend perfectly while the decentralized MPC shows a lower performance. Therefore, issues due to the convergence of iterative approach, reported by Mohseni 10 , for the open-loop unstable case study is resolved using Algorithm 5. 2.5

455 Decentralized MPC Centralized MPC CDMPC

450 445 0

0.05

0.1

0.15

0.2

A1

T1[K]

460

C [kmol/m3]

465

2 1.5

0.5 0

0.25

Decentralized MPC Centralized MPC CDMPC

1

0.05

0.1

time[s]

0.2

0.25

A2

416

C [kmol/m3]

3.5 Decentralized MPC Centralized MPC CDMPC

418

414 412 0

0.15

time[s]

420

T2[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


0.05

0.1

0.15

0.2

2.5 2 1.5 1 0

0.25


3

0.05

time[s]

0.1

0.15

0.2

0.25

time[s]

Figure 4: Comparison of state variables trajectories between the Decentralized MPC, the Centralized MPC and the CDMPC trajectories vs simulation time (Two-CSTR Process)

33



In addition, the overall objective functions of the three distributed control schemes JDC , JC , and JD are compared in Figure 6. The MAPE between the objective function of the CDMPC and the objective function of the decentralized MPC is approximately %200, which shows the significance of benefits of applying the CDMPC approach over the length of simulation. These results show that the CDMPC scheme achieves the centralized MPC performance with the lowest possible overall objective function. In other words, it improves the decentralized MPC overall performance by applying modifications to the currently installed decentralized MPC network. 0.1

8

0 0

0.05

0.1

0.15

0.2


3

CA0[kmol/m ]

Q1[kJ/hr]


0 0

0.25

0.05

0.1

time[s]

− 0.005 0

0.15

0.2

0.25

time[s]

0

4

0.05

0.1

0.15

0.2


3

CA03[kmol/m ]


Q2[kJ/hr]

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 34 of 44

0 0

0.25

0.05

time[s]

0.1

0.15

0.2

0.25

time[s]

Figure 5: Comparison of manipulated variables trajectories between the Decentralized MPC, the Centralized MPC and the CDMPC trajectories vs simulation time (Two-CSTR Process)

6Appendix S.2

Alkylation Process

In this section, the price-driven CDMPC is applied to a benzene alkylation process. The process consists of four continuous stirred tank reactors (CSTRs) and a flash tank separator, as shown in Figure 7. Pure benzene is fed via stream F1 and pure ethylene is fed via streams F2 , F4 and F6 . Two catalytic reactions take place in CSTR-1, CSTR-2 and CSTR3. In reaction 1, benzene (A) reacts with ethylene (B) and produces the desired product 34


Page 35 of 44

1200 Decentralized MPC (J

)

DC

Overall Objective Function

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


Centralized MPC (JC) CDMPC (JD) 800

400

0 0

0.05

0.1

0.15

0.2

0.25

time[s]

Figure 6: Comparison of overall objective functions between the Decentralized MPC, the Centralized MPC and the CDMPC trajectories vs simulation time (Two-CSTR Process)

ethylbenzene (C); in reaction 2, ethylbenzene reacts with ethylene to form a byproduct, 1, 3-diethylbenzene (D). In the flash tank separator, most of benzene is separated, flows overhead and recycled back to the plant. A portion of the recycle stream Fr1 is fed to CSTR4 and the other portion Fr2 is fed back to CSTR-1 together with an additional feed stream F10 , which contains 1, 3-diethylbenzene. Moreover, reaction 2 and reaction 3, in which 1, 3diethylbenzene reacts with benzene to produce ethylbenzene, take place in CSTR-4. Finally, it is assumed that throughout the process, all the materials are in liquid phase and their molar volumes are constant. The definition of the process variables is given in Tables 2. A more detailed version of the process model is presented by Christofides et al.. 25 In the design of the proposed CDMPC, the process is divided into two sub-systems. The first sub-system includes CSTR-1 and CSTR-2. The second sub-system includes CSTR-3, the separator, and CSTR-4. Consequently, the following sub-system decomposition of state and input variables is defined as:    X1 = [CA1 , CB1 , CC1 , CD1 , T1 , CA2 , CB2 , CC2 , CD2 , T2 ]   U1 = [Q1 , Q2 , F4 ] 35


(41a)


F2,B

F4,B

Page 36 of 44

F6,B

F1,A

F3

Q1

CSTR-1

F4

Q2

CSTR-2

Q3

CSTR-3 F7

Fr2 Fr

FP

Fr1

F10,D

F9

F8 Q5

CSTR-4

Q4

Separator

Figure 7: Process flow diagram of alkylation of benzene Table 2: Alkylation process parameters Process Parameter Description CA1 ,CB1 ,CC1 ,CD1 Concentration of species in CSTR-1 CA2 ,CB2 ,CC2 ,CD2 Concentration of species in CSTR-2 CA3 ,CB3 ,CC3 ,CD3 Concentration of species in CSTR-3 CA4 ,CB4 ,CC4 ,CD4 Concentration of species in Separator CA5 ,CB5 ,CC5 ,CD5 Concentration of species in CSTR-4 CAr ,CBr ,CCr ,CDr Concentration of species in Fr , Fr1 , Fr2 T1 ,T2 ,T3 ,T4 ,T5 Temperature in each vessel F3 ,F5 ,F7 ,F8 ,F9 Effluent flow rates from each vessel F1 ,F2 ,F4 ,F6 ,F10 Feed flow rates to each vessel Fr , Fr1 , Fr2 Recycle flow rates Q1 ,Q2 ,Q3 ,Q4 ,Q5 External heat/coolant input to each vessel CpA ,CpA ,CpA ,CpA Heat capacity at liquid phase

and:    X2 = [CA3 , CB3 , CC3 , CD3 , T3 , CA4 , CB4 , CC4 , CD4 , T4 , CA5 , CB5 , CC5 , CD5 , T5 ]   U2 =

(41b)

[Q3 , Q4 , Q5 , F6 ]

In this paper, the nonlinear system is linearized around the equilibrium point reported by Christofides et al. 25 . Furthermore, the inequality constraints are removed from the original case study which makes it suitable for both algorithms proposed in this paper. The sampling

36


Page 37 of 44

period is selected to be T = 30[s]; also, the weighting matrices Q and R are chosen as I25 and 102 × I7 , respectively;. Here, ‘I’ represents the identity matrix of appropriate dimension.

480

T2 [K]

T1 [K]

480

460

450

Centralized MPC / CDMPC Decentralized MPC 440 0

500

1000

1500 time[s]

2000

2500


3000

500

450

1500 time[s]

2000

2500

3000

460


1000

480

T4 [K]

480

T3 [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


500

1000

1500 time[s]

2000

2500


3000

500

1000

1500 time[s]

2000

2500

3000

Figure 8: Comparison of the temperature trajectories for CSTR-1, CSTR-2, CSTR-3, and CSTR-4 between the Decentralized MPC, the Centralized MPC and the CDMPC trajectories vs simulation time (Alkylation Process)

The temperature profiles of the CSTR units 1 to 4 are depicted in Figure 8, which show that the CDMPC approach perfectly converged to the corresponding centralized MPC optimal trajectory. Furthermore, the overall objective functions of the three distributed control schemes JDC and JC /JD are compared in Figure 9. Moreover, the evolution of prediction horizon required for the decentralized MPC and the centralized MPC/CDMPC is shown in Figure 10. The MAPE between the objective function of the CDMPC and the objective function of the decentralized MPC is approximately %32, which shows the benefits of applying the CDMPC approach over the length of simulation in terms of the overall objective function of the plant-wide problem. These results show that the CDMPC scheme improves the plant-wide performance of the previously installed decentralized MPC 37



7

14

x 10

Centralized MPC / CDMPC Decentralized MPC

Overall Objective Function

12

10

8

6

4

2

0 0

500

1000

1500

2000

2500

3000

time [s]

Figure 9: Comparison of overall objective functions between the Decentralized MPC, the Centralized MPC and the CDMPC trajectories vs simulation time (Alkylation Process)

by applying modifications to the local MPC formulation and adding a coordination layer to the network. 14 CDMPC/Centralized MPC Decentralized MPC

12

Prediction Horizon (N)

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 38 of 44

8

4

0

0

500

1000

1500

2000

2500

3000

time[s]

Figure 10: Comparison of number predictions horizons needed (N) between the Decentralized MPC, the Centralized MPC and the CDMPC trajectories vs simulation time (Alkylation Process)

38


Page 39 of 44 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


7

Conclusions

Two novel solution methodologies to solve the price-driven CDMPC problem were proposed. The analytic CDMPC scheme presented is capable of coordinating MPC networks subject to equality constraints; and, it is able to reduce the coordination cycles of conventional CDMPC methods to one. The feasible descent direction approach to solve a CDMPC problem proposed is globally convergent to an optimal plant-wide solution for the general case when inequality constraints are also present in the local MPC formulations. Furthermore, the adaptive horizon scheme enables the network to stabilize open-loop unstable systems, and resolves the convergence problems reported by Marcos et al. 7 , in conventional Newton-based nested BLP approaches. Results are promising and CDMPC can perfectly track the centralized optimum trajectory. Therefore, the performance of the decentralized MPC system is improved with modifications applied to the objective function of local MPC controllers and adding a coordination level to the network, in the price-driven CDMPC architecture. Note that, the proposed feasible decent direction approach cannot be used directly to nonlinear systems, since it requires the lower-level problem to be convex. However, it can be extended to a successively linearized version of a nonlinear CDMPC problems when there are no inequality constraints on the state trajectories, and the sub-systems are open-loop stable. In other words, the proposed schemes are expandable to open-loop stable nonlinear systems when there are no restrictions on the state trajectories of the system.

39



8

Supporting Information

The reader is encourage to find the above-mentioned proofs and verifications in the supporting information document attached to this paper. The following list summarizes the sections included in the supporting information document:

Appendix S.1

: Proof of Remark 3.1

Appendix S.2

: Proof of Remark 3.2

Appendix S.3

: Proof of Theorem 3.1

Appendix S.4


Appendix S.5

: Verification of Assumption 4.1

Appendix S.6

: Verification of Assumption 4.2

Appendix S.7


Appendix S.8


Appendix S.9


9

Acknowledgments

The authors would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support of this work.

References (1) Rawlings, J. B.; Stewart, B. T. Coordinating multiple optimization-based controllers: new opportunities and challenges. J. Process Control 2008, 18, 839–845. 40


Page 40 of 44

Page 41 of 44 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


(2) Stewart, B. T.; Venkat, A. N.; Rawlings, J. B.; Wright, S. J.; Pannocchia, G. Cooperative Distributed Model Predictive Control. Syst. Control Lett. 2010, 59, 460–469. (3) Liu, J.; Muñoz de la Peña, D.; Christofides, P. D. Distributed model predictive control of nonlinear process systems. AIChE J. 2009, 55, 1171–1184. (4) Liu, J.; Chen, X.; Muñoz de la Peña, D.; Christofides, P. D. Sequential and Iterative Architectures for Distributed Model Predictive Control of Nonlinear Process Systems. Part I: Theory. AIChE J. Baltimore, Maryland, 2010; pp 2137–2149. (5) Cheng, R.; F. Forbes, J.; Yip, W. Price-driven Coordination Method for Solving Plantwide MPC Problems. J. Process Control 2007, 17, 429–438. (6) Cheng, R.; Forbes, J. F.; Yip, W. DantzigâĂŞWolfe Decomposition and Plant-wide MPC Coordination. Comput. Chem. Eng. 2008, 32, 1507–1522. (7) Marcos, N. I.; Forbes, J. F.; Guay, M. Price-Driven Coordination of Distributed MPC Controllers for Constrained Dynamic Systems. I EC Res. 2013, 52, 17451–17464. (8) Hassanzadeh, B.; Pakravesh, H.; Liu, J.; Forbes, J. F. Distributed Model Predictive Control of Nonlinear Systems Based on Price-Driven Coordination. Industrial & Engineering Chemistry Research 2016, 55, 9711–9724. (9) Marcos, N. I.; Forbes, J. F.; Guay, M. Prediction-driven coordination of distributed MPC controllers for linear unconstrained dynamic systems. International Journal of Control 2014, 87, 1496–1512. (10) Mohseni, P. G. Coordination techniques for distributed model predictive control. PhD Dissertation, University of Alberta, 2013. (11) Hassanzadeh, B.; Liu, J.; Forbes, J. F. An analytic price-driven coordination scheme for distributed model predictive control systems. 53rd IEEE Conference on Decision and Control. 2014; pp 2461–2466. 41



(12) Savard, G.; Gauvin, J. The steepest descent direction for the nonlinear bilevel programming problem. Oper. Res. Lett. 1994, 15, 265–272. (13) Bard, J. F. Practical Bilevel Optimization : Algorithms and Applications; Kluwer Academic Publishers, 1998; pp 332–339. (14) Dempe, S. Annotated Bibliography on Bilevel Programming and Mathematical Programs with Equilibrium Constraints. Optimization 2003, 52, 333–359. (15) Colson, B.; Marcotte, P.; Savard, G. Bilevel Programming: A Survey. 4OR A Q. J. Oper. Res. 2005, 3, 87–107. (16) Thoai, N. V.; Yamamoto, Y.; Yoshise, A. Global Optimization Method for Solving Mathematical Programs with Linear Complementarity Constraints. J. Optim. Theory Appl. 2005, 124, 467–490. (17) Wan, Z.; Wang, G.; Sun, B. A hybrid intelligent algorithm by combining particle swarm optimization with chaos searching technique for solving nonlinear bilevel programming problems. Swarm Evol. Comput. 2013, 8, 26–32. (18) Pannek, J. Receding Horizon Control: A Suboptimality-based Approach. PhD Dissertation, University of Bayreuth, 2009. (19) Chen, Y.; Florian, M. The Nonlinear Bilevel Programming Problem: Formulations, Regularity and Optimality Conditions. Optimization 1995, 32, 193–209. (20) Facchinei, F.; Jiang, H.; Qi, L. A Smoothing Method for Mathematical Programs with Equilibrium Constraints. Math. Program. 1999, 85, 107–134. (21) Jiang, Y.; Li, X.; Huang, C.; Wu, X. An augmented Lagrangian multiplier method based on a CHKS smoothing function for solving nonlinear bilevel programming problems. Knowledge-Based Syst. 2014, 55, 9–14.

42


Page 42 of 44

Page 43 of 44 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


(22) Grüne, L.; Pannek, J. Practical NMPC suboptimality estimates along trajectories. Syst. Control Lett. 2009, 58, 161–168. (23) Grüne, L.; Pannek, J. Nonlinear Model Predictive Control: Theory and Algorithm; Springer-Verlag London Limited, 2011; pp 113–209. (24) Sun, Y.; El-Farra, N. H. Quasi-decentralized model-based networked control of process systems. Comput. Chem. Eng. 2008, 32, 2016–2029. (25) Christofides, P. D.; Liu, J.; Muñoz De La Peña, D. Networked and distributed predictive control: methods and nonlinear process network applications; Advances in Industrial Control Series; Springer-Verlag, London, England, 2011. (26) Löfberg, J. YALMIP : A Toolbox for Modeling and Optimization in MATLAB. IEEE Int. Symp. Comput. Aided Control Syst. Des. Taipei, Taiwan, 2004; pp 284–289.

43



10

Table of Contents Graphic

Decentralized MPC: MPC 1

MPC i

sub-system 1

MPC m

sub-system i

sub-system m Interactions

Interactions

Process

Coordinated Distributed MPC: CDMPC Coordinator

CDMPC Controller 1

CDMPC Controller i

sub-system 1

CDMPC Controller m

sub-system i


Interactions

Process

Centralized MPC: Centralized MPC

sub-system 1

sub-system i


Interactions

Process

TOC: The abstract graphical representation of the decentralized MPC, the proposed coordinated distributed MPC, and the centralized MPC architectures

44


Page 44 of 44

A Bilevel Optimization Approach to Coordination of Distributed Model

Recommend Documents