Assessing Model Prediction Control (MPC) - American Chemical Society

Nov 1, 2007 - Assessing Model Prediction Control (MPC) Performance. 1. Probabilistic Approach for Constraint Analysis. Nikhil Agarwal and Biao Huang*...
0 downloads 8 Views 319KB Size
Ind. Eng. Chem. Res. 2007, 46, 8101-8111

8101

Assessing Model Prediction Control (MPC) Performance. 1. Probabilistic Approach for Constraint Analysis Nikhil Agarwal and Biao Huang* Department of Chemical and Materials Engineering, UniVersity of Alberta, Edmonton, AB, Canada, T6G 2G6

Edgar C. Tamayo Syncrude Canada Ltd., Fort McMurray, Alberta, Canada, T9H 3L1

Advanced process control (APC)sin particular, model predictive control (MPC)shas emerged as the most effective control strategy in process industry, and numerous applications have been reported. Nevertheless, there are several factors that limit the achievable performance of MPC. One of the limiting factors considered in this paper is the presence of constraints. To exploit optimal control performance, continuous performance assessment, with respect to the constraints of MPC, is necessary. MPC performance assessment has received increasing interest, both in academia and in industry. This paper is concerned with a practical aspect of performance assessment of industrial MPC by investigating the relationship among process variability, constraints, and probabilistic economic performance of MPC. The proposed approach considers the uncertainties induced by process variability and evaluates the economic performance through probabilistic calculations. It also provides a guideline for the constraint tuning, to improve MPC performance. 1. Introduction Advanced process control (APC) applications such as model predictive control (MPC) emerged in mid-1970s. Since then, they have been regarded as the most popular industrial advanced multivariable control strategy. The main objective of the MPC controller is to calculate the control signals and to minimize the sum of the squares of the error between the reference signal and the output signal within a given future horizon.2 In addition to the dynamic control objective, an economic objective function is utilized to optimize steady-state economic performance, subject to process constraints. The economic optimization drives the controlled variables (CVs) and manipulated variables (MVs) to their steady-state optimum, thus generating economic benefits. The economic benefit analysis of implemented process control projects has attracted much interest.20 Despite the great success of MPC, there are certain limitations on the achievable performance. For example, the performance of practical MPC applications can be limited by model uncertainty,4,12 as well as possible conservative tuning such as conservativeness in setting up constraints.16 The latter is not uncommon and has been the motivation of this work. Many commercial controllers are now available that are able to improve MPC controller performance by resolving one or another of the practical issues.14 With most of the research being focused on the areas of developing MPC controllers, relatively less work has been reported in the field of evaluation and the practical tuning of MPC controllers. Conventional tuning for the MPC controller involves tuning the controller, with respect to the prediction and control horizon and the linear and quadratic objective functions.9,11 Guidelines are available in the literature, which provide qualitative relationships between these parameters; however, the quantitative relationships differ from one system to another. Tuning of the controllers through these parameters is generally done at the design stage and is not advised to be changed on a day-to-day basis. To change these parameters, it is essential to have thorough process knowledge, * To whom correspondence should be addressed. Tel.: +1-780-4929016. Fax: +1-780-492-2881. E-mail address: [email protected].

a deep understanding about the MPC control algorithm and the effect of these changes on the MPC controller performance. In daily operations, on the other hand, it is not uncommon to change some of the CV and MV constraint limits. This action is simple but affects the operating points of the CVs and the MVs. Sometimes it can improve the controller performance, and other times may have little effect. Thus, understanding the relationship between constraints, process variability, and economic objective functions is important to guide the constraint tuning. Change of constraints is a decision-making process in the economic optimization layer. In the words of Hummel et al.,6 “Plant optimization can be defined as the maximization of the profit function describing the economics of the plant. This function contains terms with product values, feedstock process, operational costs, etc”. Prett and Garcia13 states that the decision-making is a multilayered process, involving measurement, controls, optimization, and logistics. Measurement is the gathering of information of the process measurements. Control is the manipulation of the process degrees of freedom for satisfying the operation criteria. Optimization is defined as plant optimization, which is a technique of manipulating the process degrees of freedom to satisfy the plant economic objectives. Logistics is the highlevel scheduling to respond to external market changes for profit maximization. The optimization layer contributes the most to the profit improvement. This layer receives the inputs of economic targets from the logistics layer and then passes the optimum operating target to the control layer for realization. The previous work in the direction of performance assessment, with respective to constraints analysis, has focused primarily on evaluating the deterministic linear and quadratic objective functions,19 subject to process variability. An economic performance assessment of MPC based on a linear matrix inequality approach (LMIPA) has been developed.19 In this early work, the economic function was evaluated by considering the mean operating points of the process variables, subject to the

10.1021/ie070417e CCC: $37.00 © 2007 American Chemical Society Published on Web 11/01/2007

8102

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Figure 3. High probability of violating the lower limit.

Figure 1. True process data.

Figure 4. Lower probability of violating the lower limit.

Figure 2. Process data distribution.

constraints on the CVs/MVs, the process variance, and the steady-state process gains. However, because of the presence of the disturbances, there is a variability associated with the process operating points, and, consequently, process variables do not always operate at the mean operating points. It is necessary to consider the probabilistic objective function, which involves the uncertain operating points of the CVs and MVs, the profits associated with each operating points, and the probabilities for CVs and MVs to be in each of the operating points. The problem thus formulated will be nonlinear, and the approach taken in this paper is a step toward a more-realistic assessment of MPC control performance. The contributions of this paper can be summarized as follows: (1) A probabilistic approach for MPC performance assessment via constraint analysis is proposed. (2) Guidelines for optimal constraint tuning are derived. (3) Detailed simulations, as well as industrial case studies, are presented. The paper is organized as follows. Section 2 introduces the problem formulation and solutions. The formulation of the probabilistic objective function and the algorithms used for the probabilistic performance assessment are discussed in section 3. Section 4 illustrates an application for a simulated binary distillation column, and application to an industrial distillation column is discussed in section 5, followed by the conclusion in section 6. 2. Probabilistic Performance Analysis 2.1. Mean Operating Point. The previous MPC performance assessment, with respect to constraint analysis, calculates the economic optimum based on the mean operating point.19 Because the plant seldom operates precisely on the mean operating point, it is essential to consider the data distribution while evaluating the economic objective function. Thus, the probabilistic performance assessment method must be developed, which takes into consideration the probability distribution. To illustrate this, consider two outputs (y1 and y2) of a system. Let y1 be the quality variable and y2 be a constraint variable. The current operating data, which are defined as the base case operating data for y1 and y2, are shown in Figure 1. Figure 2 shows the probability distribution for the base case data when the Gaussian distribution is to be approximated.

As can be seen from Figures 1 and 2, although the mean operating points for y1 and y2 are within the constraint limits throughout the period for which the data are collected, the data trend shows that most of time they are not on the mean operating points, and sometime they are even outside the constraint limits. Because y1 is a quality CV, having its value outside the set limits is highly undesirable, because it can render the product unmarketable. Considering only the mean operating point for performance assessment, the performance of the controller for the aforementioned data set can be called satisfactory, but when also considering the data distribution, the aforementioned performance does not necessarily need to be satisfactory. There exists an optimum operating point at which the value of the expected economic objective function can be optimized, by taking into account the data distribution. This can be explained as follows. Assume that the target operating point for y1 is at its low limit. Now moving its mean operating point closer to the low limit will not necessarily increase the profit margins. This is true because, due to the presence of variance, moving the operating point closer to the low limit will also increase the probability of violating the low limits and, thus, a tradeoff must be considered while moving closer to the low limit (see Figures 3 and 4). 2.2. Effect of Changing Constraints. The limits for the CVs and the MVs and their variability determine the optimum operating point for the MPC controller. Relaxing the limits for one or more process variables provides the controller with an increased number of degrees of freedom and, thus, may help to improve the expected return, even if there is no reduction in the variability of the variables. This is indicated in Figure 5. Figure 5 shows that, although y1 is far from the optimum operating point, the high limit, it may not be further raised to an optimum value, because y2 is hitting the constraints at the high limit. Relaxing the constraint for y2 on the high limits may provide the controller with additional degrees of freedom and y1 may be raised to a point closer to the optimum operating point, leading to improvement of the controller performance. Therefore, the objective of this study is to investigate the potential of increasing the expected return by changing the operating points or changing constraints on selected CVs/MVs. The optimal operating conditions are identified by optimizing the probabilistic economic objective function, subject to the constraints and variability.8,10 Thus, this paper will answer the following questions: (1) Given the existing process variability and constraints, is the process operating at its economic optimum? If not, what is its potential to improve and which variables must be adjusted, relative to their operating points, to achieve the optimum?

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8103

Figure 6. Division of zones for a child node. Figure 5. Effect of change in limits.

(2) Given the existing process variability and possible adjustment of constraints for some variables, how much additional performance can be achieved? Which variables must be adjusted by their constraints, and by how much to achieve the performance improvements? All performance potentials are calculated by taking into account of the process variability through probability distribution functions. 2.3. Probabilistic Optimization. 2.3.1. Optimization Formulation. For an m × n system with n inputs and m outputs, let K be the steady-state gain matrix and (yji0,ujj0) be the current mean operating points for the ith CV and jth MV, i.e., the base case mean operating points. Let (yji,ujj) be the corresponding optimal operating points to be determined. For simplicity of the presentation, we assume that the profit is dependent only on the CVs (i.e., the CVs are quality variables). The results can be easily extended to include MVs as quality variables. Without any loss of generality, let the first q CVs be the quality variables. Given the mean operating points and the variance of the process data, the expected return, for a process with q quality variables, can be calculated as shown in eq 1:

E(R) )

∫y ,...,y 1

q

p(y1,...,yq)F(y1,...,yq) dy1 ... dyq

(1)

where p(y1,...,yq) is a probability density function, and F(y1,...,yq) is a profit function. Generally, quality CVs directly reflect the economic returns of the process and are very selective; thus, q is typically a small number (for example, purity of the product, which reflects economic return, whereas others are controlled and/or constraint variables). Economic profit functions are often expressed as a linear/quadratic function, such as q

profit )

(Riyi + βiyi2 + µi) ∑ i)1

(2)

where µi is the nominal return for the ith CV, and Ri and βi are the linear and quadratic coefficients for yi, respectively. In this paper, we assume that the optimum of the objective function occurs either at the lower or upper constraint limit; namely, within the constraint limits, the quadratic function is either monotonically increasing or decreasing. With a linear or

quadratic economic objective function, F(y1,...,yq) has the following additive form: q

F(i)(yi) ∑ i)1

F(y1,...,yq) )

(3)

where

F(i)(yi) ) Riyi + βiyi2 + µi Therefore,

E(R) )

∫y ,...,y 1

q

(

p(y1,...,yq)

q

)

F(i)(yi) ∑ i)1

dy1 ... dyq

(4)

The economic return has often been classified as in-spec, underspec, and over-spec;15 therefore, only three operating zones are needed to specify the profit. To increase the resolution, each CV can be discretized into a finite number of operating zones. For illustration purpose, six zones (Zone 1, Zone 2, Zone 3, Zone 4, Zone 5, Zone 6) are considered, which define the range in which the value of the output variables will lie, as illustrated in Figure 6. The state space of each CV (e.g., the ith CV) then can be written as

()( )

Ω1 Zone 1 Ω2 Zone 2 Ω3 Zone 3 i Ω ) Ω ) Zone 4 4 Ω5 Zone 5 Ω6 Zone 6

Namely, yi will take one of the six discrete value in Ωi. The number of zones can certainly be increased, depending on the resolution that is required. Zones 1 and 6 represent the region below and above the limits, respectively, and Zones 2-5 represent the four zones that are defined within the limits. Thus, if Lyi and Hyi are the low and the high limits for a particular CV (yi), then ∆i is the span of the range in which its value is to be maintained, i.e.,

∆i ) Hyi - Lyi The six zones for the states of the CV then can be defined as follows:

Zone 6 ) Hyi to ∞

(5)

8104

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Zone 1 ) -∞ to Lyi

(

Zone 2 ) Lyi to Lyi +

( ( (

)

∆i 4

) ( ) ( )

)

Zone 3 ) Lyi +

∆i ∆i to Lyi + 4 2

Zone 4 ) Lyi +

∆i 3∆i to Lyi + 2 4

Zone 5 ) Lyi +

3∆i to Hyi 4

) Figure 7. Probability for yi to be in the six zones.

By discretization, a simple discrete version of eq 1 can then be written as

E(R) )



P(y1,...,yq)F(y1,...,yq)

(6)

y1,...,yq

where P(y1,...,yq) is now a probability function. Furthermore, if the uncertainties associated with each of y1,...,yq are mutually independent, and the profit function is additive shown in eq 3, then eq 6 can be further simplified to q

E(R) )

6

∑ ∑ P(yi ∈ Ωk) × F(yi ∈ Ωk) i)1 k)1

(7)

(The independence of y1,...,yq, in terms of their uncertainties, is often observed, because the selected economic quality variables reflect distinct aspects of process operations.) As indicated by Figure 7, the probability for the quality variables to be in each of the six zones is estimated by assuming the data to be Gaussian-distributed with a mean yji and a covariance that is calculated from base operation data. Because we change only the constraints, the covariance remains the same as that of the base-case operation, because the dynamic control tuning is not being changed. Optimizing the economic objective function, subject to the original/new constraints (if constraints are to be adjusted), provides the optimum operating points that will provide the new mean operating points yji for the distribution function. The new mean operating points and base-case covariance, based on the multivariate Gaussian distribution, are used to obtain the probabilities for the CVs to be in each of the six zones defined for them. For the simpler case where the uncertainties associated with the quality variables are mutually independent or there is only one economic quality variable, eq 8 can be applied:

P(y1,...,yq) ) P(y1)...P(yq) where

P(yi ∈ Ωk) )

1 x σi 2Π 0

∫L

Hkyi kyi

(

exp -

(8)

(yi - yji)2 2

2σi0

)

The return for yi to be obtained in each of the six zones can be specified according to economic data. The following is an example of this calculation. The maximum return is obtained when the process operates at the maximum-return operating points (zones), and the maximum-return points typically lie in one of the constraint limits. Thus, the maximum return can be either at the high constraint limit or at the low constraint limit. The minimum return will lie in the other end of the constraint or in the zone of under-spec/over-spec. Given the maximum return and the minimum return, the profit for yi in the in-spec zones can be estimated through interpolation and be calculated as given below. Let Fk denote the return for the process variable yi to operate in the kth zone. If yi is to be maximized (i.e., the maximum return is at the high constraint limit, i.e., Fmax ) F5, then Fmin ) F1):

1 ∆F ) (F5 - F1) 4 F2 ) F1 + ∆F F3 ) F1 + 2∆F F4 ) F1 + 3∆F

(10)

This calculation is illustrated in the left panel of Figure 8. If yi is to be minimized (i.e., the maximum return is at the low constraint limit, i.e., Fmax ) F2, Fmin ) F6), then

1 ∆F ) (F2 - F6) 4

(11a)

F5 ) F6 + ∆F

(11b)

F4 ) F6 + 2 × ∆F

(11c)

F3 ) F6 + 3 × ∆F

(11d)

This calculation is illustrated in the middle panel of Figure 8. However, if yi is neither to be maximized or minimized, then

dyi

(9)

where yji is the mean operating point, σi0 is the standard deviation of yi in base-case operation, i ) 1, 2, ..., q, and k ) 1, 2, ..., 6. Lkyi and Hkyi are the low and the high limits for yi in the kth zone and can be identified as described in eq 5. In the rest of this paper, we focus on the case of either q ) 1 or the quality CVs being mutually independent. If the quality CVs are dependent, then eq 6 should be used for the profit calculation, which may increase the complexity, in the sense of numerical computations.

F 2 ) F3 ) F4 ) F5 ) F

(12)

This calculation is illustrated in the right panel of Figure 8. Note that, in some cases, F1 or F6 may simply be a loss and may not reflect the minimum return. Thus, F1 and F6 are shown as dotted lines in Figure 8, which indicates that these may have values less than, greater than, or equal to the returns in Zones 2 and 5, respectively, depending on the specific processes. 2.3.2. Optimization by Changing Operating Points. Not all processes operate at their optimum, even if there is potential to do so. Thus, it is possible to increase the profit by simply adjusting the operating points without doing any other tunings.

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8105

Figure 8. Return for yi to be in the six zones.

Figure 9. Binary distillation column.

There is a need for assessment of performance potential. The potential for expected return, E(Rp), by simply adjusting the operating points, without any other tuning such as constraint limits relaxation, is then estimated through the maximization of the expected economic objective function (eq 7). Thus, the objective function for the optimization can be specified as q

J)

6

∑ ∑ P(yi ∈ Ωk) × F(yi ∈ Ωk) i)1 k)1

(13)

With (yji0, uji0) as the base-case mean operating point, and (yji, uji) as the new mean operating point when the base-case operating points are adjusted by (∆yi, ∆uj), then the aforementioned maximization should be subject to the following equality constraints:

∑ j)1

Lyi + 2σi0 e yji e Hyi - 2σi0

(17)

Luj + 2Rj0 e ujj e Huj - 2Rj0

(18)

where i ) 1, 2, ..., m and j ) 1, 2, ..., n. σi0 is the standard deviation for the base operation data yi, Rj0 is the quarter of the existing operating range for the base operation data uj, and Luj and Huj are the low and the high constraint limits for uj, respectively. Lyi and Hyi are the low and the high constraint limits for yi, respectively. Thus, the performance assessment problem can be specified as

max J

n

∆yi )

limits for yi and uj can be identified. The following inequalities19 are required to be satisfied while optimizing the objective function defined in eq 13:

[Kij × ∆uj]

(14)

yji ) yji0 + ∆yi

(15)

ujj ) ujj0 + ∆uj

(16)

where, i ) 1, 2, ..., m and j ) 1, 2, ..., n, and K is the gain matrix of the process. Considering that up to 5% constraint violation is acceptable for the output variables,8,10 a set of inequalities that defines the

jy1,...,yjm uj1,...,ujn

(19)

s.t. eqs 14-18 The value of J thus estimated through the objective function defined in eq 13, subject to the constraints defined by eqs 1418, gives the potential to improve returns by simply adjusting the operating points without adjusting constraints or doing any other tuning. 2.3.3. Optimization by Tuning Constraints. In practice, because of some conservativeness or lack of information on

8106

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

the potential economic impacts of certain constraints when setting up their limits, the constraint limits for some CV and MV are often adjusted during process operations. Therefore, there is an incentive to determine the best adjustment of constraint limits, according to the base-case operation data. The results obtained can serve as a powerful tool set for process control engineers when they make adjustments of the constraints. On the other hand, even if there is no intention to change the constraints, the results can also provide an understanding of the sensitivity of the constraint limits of each variable in terms of economic return. These problems are considered next. Let B be the maximum allowable change in the constraint limits that can be made for a CV or MV, and let r (r e B) be the actual adjustment of the limits needed for that CV or MV to yield the maximum returns. Under this scenario, the expected return E(Rc) is given by the objective function described by eq 13, subject to the equalities specified in eqs 14-16; however, the inequalities given by eqs 17 and 18 change to eqs 20 and 21, respectively.

Lyi + 2σi0 - yholi × rLyi e yji e Hyi - 2σi0 + yholirHyi

(20)

Luj + 2Rj0 - uholj ×

(21)

rLuj

e ujj e Huj - 2Rj0 +

uholjruHj

The additional constraints are given by eqs 22-25.

0 e rLyi e BLyi

(22)

0 e ryHi e BHyi

(23)

0 e rLuj e BLuj

(24)

0 e ruHj e BuHj

(25)

where i ) 1, 2, ..., m and j ) 1, 2, ..., n and yholi, uholj are half of the range of the constraint limits for yi and uj; ByLi, ByHi , BuLj, and BuHj are the maximum allowable change (given as a percentage of the range specified) in the limits for the CVs and MVs respectively; ryLi, ryHi , ruLj, ruHj are the actual limit changes for the CVs and MVs, respectively. The subscripts “yi” and “uj” represent the ith output and jth input, respectively; the superscripts “L” and “H” represent the low and high limits, respectively. Thus, the performance assessment problem can be specified as

max J

jy1,...,yjm uj1,...,ujn ryL1,...,ryLm

(26)

Table 1. Assumed Control Objective and the Profits for the CVs Profit CV

description

objective

u/s

max

o/s

1 2 3 4 5 6 7 8 9 10

reflux flow light petrol FBP top percentage Pr Vlv OP bottoms percentage column pressure feed temperature reboiler duty duty bypass Vlv OP

constraint minimize constraint constraint constraint constraint constraint minimize constraint constraint

0 65 0 0 0 0 0 200 0 0

0 65 0 0 0 0 0 200 0 0

0 0 0 0 0 0 0 0 0 0

Table 2. List of MVs MV

description

1 2 3 4

reflux SP column pressure SP feed temperature Vlv OP duty Vlv OP

Table 3. Maximum Constraint Change Allowed Maximum Constraint Allowed in low limit

in high limit

0 10 0 0 10 10 10 0 0 0 10 10 10 0

0 10 0 0 0 0 0 0 0 0 10 10 50 0

CV1 CV2 CV3 CV4 CV5 CV6 CV7 CV8 CV9 CV10 MV1 MV2 MV3 MV4

For the purpose of explaining the calculations involved for the optimization, the terms involved in the objective function (eq 13) are elaborated below: The terms F(yi ∈ Ωk) are calculated using eqs 10, 11, or 12, and the probabilities are given by

P(yi ∈ Ω1) ) P(yi < Lyi) P(yi ∈ Ω6) ) P(yi g Hyi) P(yi ∈ Ωk) ) P(Lkyi e yi < Hkyi) where k ) 2, 3, 4, 5; Lyi and Hyi are the low and high constraint limits for yi; Lkyi and Hkyi are the low and high limits for the kth zone of yi. Because the process data have a Gaussian distribution, the aforementioned probabilities, in terms of the cumulative probability function, can be written as

ryH1,...,ryHm rLu1,...,rLun H rH u1,...,run

s.t. eqs 14-16 and 20-25 3. Algorithm The procedure discussed in the previous section provides the solutions to the problem using constraint limits relaxation for the CVs and MVs. The objective function that has been defined for this purpose is nonlinear, whereas the constraint equalities and the inequalities are linear.

P(yi < Lyi) )

[ ( )] [ ( )]

Lyi - yji 1 1 + erf 2 σi0x2

P(yi g Hyi) ) 1 -

Hyi - yji 1 1 + erf 2 σi0x2

P(Lkyi e yi < Hkyi) ) P(yi e Hkyi) - P(yi < Lkyi)

(27a)

(27b) (27c)

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8107

Figure 10. Comparison of the expected return. Table 4. Comparison of the Three Scenarios expected return POa

scenario base case base-case potential constraint adjustment a

145.36 190.64 192.38

Figure 11. Constraint tuning guidelines (expressed as a percentage).

Probabilistic optimization.

The inequalities specified by eqs 20-25 can be reformulated, in the Ax e b format, as

Table 5. Suggested Constraint Relaxtion Suggested Constraint Relaxation (%) in low limit

in high limit

0 5.00 0 0 3.20 10.00 9.99 0 0 0 10.00 10.00 0 0

0 10.00 0 0 0 0 0 0 0 0 10.00 10.00 50.00 0

CV1 CV2 CV3 CV4 CV5 CV6 CV7 CV8 CV9 CV10 MV1 MV2 MV3 MV4

Table 6. Control Objective and the Assumed Profits for the CVs Profit CV

description

objective

u/s

Max

o/s

1 2 3 4 5 6 7 8

SSP Flow naphtha 95%Pt SSP 5% Pt OVHD Boot Lvl OP column flood PA return temperature steam stream temperature column top temperature

constraint minimize constraint constraint constraint constraint constraint constraint

0 115.00 0 0 0 0 0 0

0 124.00 0 0 0 0 0 0

0 0 0 0 0 0 0 0

where erf(x) is defined as

2 erf(x) ) xΠ

∫0

x

2

exp(-t ) dt

n

Kij × ∆uj ∑ j)1 n

) yji0 +

Kij × (ujj - ujj ) ∑ j)1 0

n

) (yji0 -

n

Kij × ujj ) + ∑ Kij × ujj ∑ j)1 j)1

where i ) 1, 2, ..., m.

0

(30)

yholi × ryHi + yji e Hyi - 2σi0

(31)

-uholj × ruLj - ujj e -Luj - 2Rj0

(32)

uholj × ruHj + ujj e Huj - 2Rj0

(33)

-ryLi e 0

(34)

ryLi e ByLi

(35)

-ryHi e 0

(36)

ryHi e ByHi

(37)

-ruLj e 0

(38)

ruLj e BuLj

(39)

-ruHj e 0

(40)

ruHj e BHuj

(41)

(28)

Also, the equality constraints for the problem, as defined by eqs 14-16, can be simplified as

yji ) yji0 +

-yholi × ryLi - yji e -Lyi - 2σi0

(29)

where, in this algorithm, routine operating process data are used to calculate the base-case mean operating point (yji, ujj) and the base-case variance (σi02). The optimization function is now defined as given in eq 13, the equality constraint (eq 29), and the inequalities (eqs 3041) for the constraint limits relaxation case. The decision variables for the objective function are [ryLi, ryHi , ruLj, ruHj, jyi, ujj], where i ) 1, 2, ..., m and j ) 1, 2, ..., n. The values of these decision variables can now be calculated to optimize the objective function while satisfying the constraints. For the case that there is no constraint to be adjusted, the optimization algorithm can be formulated in a similar way. The objective function defined is nonlinear, and the constraints are linear. Sequential quadratic programming (SQP) is applied to solve the problem. In the SQP method, the nonlinear

8108

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Figure 12. Industrial distillation column.

optimization function is quadratically approximated, by its Taylor series, and then the approximated function is optimized.3,5,7,18 4. Simulation Study of a Binary Distillation Column The column (see Figure 9) is a binary distillation column that is used to separate light petrol and the heavy petrol from the feed obtained from an upstream desulfurization unit.17 The light petrol is comprised of components in the boiling range of 30-65 °C, and the heavy petrol has hydrocarbon components in the boiling range of 65-180 °C. The feed to the column is heated by steam and flashed into the column. The lighter fractions of the hot feed vaporize and are collected as the top product in the overhead vessel, and the heavier components are obtained at the bottoms of the column. The vapors from the top of the column are cooled and condensed by the air fin coolers. A portion of the condensed overhead vapors are sent back into the column as the reflux and the remainder is withdrawn, under a level control, as the top product (i.e., the light petrol). The reflux helps to maintain the column top temperature and also the vapor liquid traffic in the reflux zone of the column, which assists in better fractionation. The column operating pressure is maintained by a split range control, which adjusts the cooling of the overhead vapors and, if required, vents off the incondensable components from the column (as the off gas). The heavier components in the feed do not vaporize and are obtained from the column bottoms. A portion of the bottom stream is sent to the reboiler, where it is reboiled by a heating medium using a duty controller. The heating medium is heated in the reboiler furnace, whose duty is required to be maintained constant. The vapors from the reboiler are sent back into the column, and they assist in stripping off any lighter components that could not flash into the flash zone of the column and are carried to the bottom of the column. The balance of the heavier components is withdrawn from the column as the heavy petrol, under a bottom level controller. The controller designed for the column has 10 CVs and 4 MVs. The list of CVs, along with the assumed profit associated with them (identified as under-spec (u/s), in-spec (maximum), and over-spec (o/s)) and their control objective are listed in Table 1. The list of MVs is given in Table 2. As CV2 and CV8 govern the overall economics of the operations of the column, they are

chosen as the quality CVs here. The correlation coefficient between these two CVs is 0.0049, according to the base operation simulation data, and, thus, these CVs can be considered to be independent. Table 1 shows that CV2 and CV8 are the quality variables to be minimized, and the profits that are associated with these CVs, which are in the six zones, can be identified using eq 11. Note that fewer than 2.5% points of CV2 and CV8 would be operating under-spec under the optimal operating conditions owing to the constrain specification in eq 17 or eq 20. To optimize the process, with regard to the constraint limits change, subject to the maximum constraint change allowed for the process variables listed in Table 3, the optimization problem can be defined as described in eq 42. 6

max J )

jy1,...,yj10 uj1,...,uj4 ryL1,...,ryL10

∑ ∑ P(yi ∈ Ωk) × F(yi ∈ Ωk)

(42)

i)2,8 k)1

ryH1,...,ryH10 rLu1,...,rLu4 H rH u1,...,ru4

s.t. eqs 14-16 and eqs 20-25 The expected return estimated for the base-case operation, the potential by adjusting operating points, and the potential by constraints relaxation are calculated using the proposed probabilistic optimization (PO) method. The results obtained are listed in Table 4. The results for the three scenarios of the PO method are also shown in Figure 10. Figure 10 shows that the expected return of base-case operation is 145.36 units, whereas the potential return by adjusting the operating points is 190.64 units, and the potential return by adjusting constraints is 192.38 units. To achieve the expected returns of 192.38 units, the optimizer has suggested constraint relaxations for certain variables, as listed in Table 5 and shown in Figure 11. 5. Case Study of an Industrial Distillation Column The feed to the column of interest is dilute bitumen, which is to be fractionated to obtain the top product (naphtha), the middle distillate (called side striper product (SSP)), and the column bottoms (bitumen). The feed is heated by steam and

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8109 Table 7. List of MVs MV

description

1 2 3 4 5

column pressure SP reflux flow SP PA return temperature column stream steam side stream steam

Table 8. Maximum Constraint Change Allowed Maximum Constraint Change Allowed CV1 CV2 CV3 CV4 CV5 CV6 CV7 CV8 MV1 MV2 MV3 MV4 MV5

in low limit

in high limit

0 0 0 10 10 10 10 10 10 10 10 10 10

10 0 0 10 10 10 10 10 10 10 10 10 10

Figure 13. Comparison of the expected return.

Table 9. Comparison of the Three Scenarios scenario

expected return

base case base-case potential constraint adjustment

123.70 123.95 123.99

Table 10. Suggested Constraint Relaxation

Figure 14. Constraint tuning guidelines (expressed as a percentage).

Suggested Constraint Relaxation (%) CV1 CV2 CV3 CV4 CV5 CV6 CV7 CV8 MV1 MV2 MV3 MV4 MV5

in low limit

in high limit

0 0 0 0 4.08 1.82 3.33 6.99 5.00 3.72 6.32 5.00 5.00

9.99 0 0 9.99 4.77 6.47 3.60 6.01 4.99 4.99 5.00 3.87 6.36

routed into a primary flash drum. The vapors from the vessel are cooled and collected as one of the streams to the naphtha pool, and the liquid stream is further heated with steam, followed by heat recovery from the pump-around from the column. The heated feed is then flashed in the secondary pre-flash vessel for recovery of the vaporized naphtha, and the liquid stream from the drum is again heated by the hot SSP product, followed by further heating in the furnace. The hot stream from the furnace is then flashed in the flash zone of the distillation column from which the balance naphtha is recovered as overhead stream, SSP is drawn as the side draw, and the bitumen is obtained as the column bottoms. The stripping steam is also fed into the stripping zone of the column to recover any SSP that is carried away with the heavier components. The SSP from the draw tray is fed into the side stripper, where it is stripped off any lighter components by the side stripper stream. The pumparound stream of SSP is also drawn from the column to remove any excess amount of heat from the column. The typical flow diagram for the distillation plant previously described is shown in figure 12. The MPC controller for the column has five inputs (or MVs) and eight outputs (or CVs).

Table 6 lists the CVs for the controller, along with their control objectives and the profits associated with under-spec (u/s), maximum, and over-spec (o/s) conditions. Table 7 gives the list of MVs. Because CV2 reflects the product quality and is considered for economic assessment in industry, it is used as the quality variable. The process is subject to various measured and unmeasured disturbances (such as fuel gas pressure fluctuation, changes in feed composition, ambient temperature changes, fluctuation in steam header temperature and pressure) which have a tendency to move the quality variable away from the optimum operating point and sometimes also outside the constraint limits. The manipulated variables are column pressure, column reflux, pump-around flow, column steam, and side stripper steam. The maximum allowable constraint changes for the process variables are listed in Table 8, and the optimization problem is defined in eq 43. 6

max J )

jy1,...,yj8 uj1,...,uj5 ryL1,...,ryL8

∑ P(y2∈ Ωk) × F(y2∈ Ωk)

(43)

k)1

ryH1,...,ryH8 rLu1,...,rLu5 H rH u1,...,ru5

s.t. eqs 14-16 and eqs 20-25 The expected return estimated for the base-case operation, potential return by adjusting operating points, and potential return by constraints relaxation are calculated by the proposed PO method, and the results are listed in Table 9. Figure 13 compares the results for the three scenarios of the PO method.

8110

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Figure 13 shows that the expected return of base-case operation is 123.70 units, whereas the potential return by adjusting the operating points (E(Rp)) is 123.95 units, and the potential return by adjusting constraints (E(Rc)) is 123.99 units. To achieve a return of 123.99 units, the optimizer has suggested constraint relaxation for some variables. The constraint relaxations suggested by the algorithm are listed in Table 10 and are shown in Figure 14. Because the three values of the expected return presented for the analysis are very similar to each other, it can be concluded that, with the constraint limits and allowable changes provided, the controller has already been operating very close to the optimal operations. There is little benefit to do any adjustment. 6. Conclusion A performance assessment and constraint tuning method for the model prediction control (MPC) systems, which takes into account the distribution of the process data and the profits associated with the product quality, has been developed in this paper. The real processes always involve uncertainties; therefore, it is pertinent to develop a probabilistic approach for MPC performance assessment. The developed method calculates the expected return of base-case operation, and the potential improvement of the expected return through adjusting the operating points, as well as through adjusting constraints. It furthermore provides tuning guidelines and determines which variables must be adjusted and by how much. Two case studiessone for a simulated binary distillation column and the other for an industrial distillation columnshave been provided that illustrate the industrial utility of the proposed algorithm. The results from the studies show the feasibility of the proposed algorithm and also bring forth its utility for the process control engineers for the day-to-day maintenance of MPC controllers in a plant and for economic assessment of the performance of the controller. Notation BuHj ) maximum limits change (as percentage of original range) for uj high limit BuLj ) maximum limits change (as percentage of original range) for uj low limit ByHi ) maximum limits change (as percentage of original range) for yi high limit ByLi ) maximum limits change, as percentage of original range, for yi low limit E(Rc) ) expected returns potential by adjusting the constraints E(Rp) ) expected returns potential for the base case condition F(yi ∈ Ωk) ) Return associated with yi in Zone k Hkyi ) high limit of the kth zone for yi Huj ) high limit for uj Hyi ) high limit for yi J ) objective function K ) steady-state process gain matrix Lkyi ) low limit of the kth zone for yi Luj ) low limit for uj Lyi ) low limit for yi m ) number of output variables n ) number of input variables N ) total number of changeable CVs and MVs Pci(x) ) cumulative probability for yi e x P(yi ∈ Ωk) ) probability for yi to be in Zone k q ) number of quality variables ruHj ) optimum limits change (as percentage of original range) for uj high limit

ruLj ) optimum limits change (as percentage of original range) for uj low limit H ryi ) optimum limits change (as percentage of original range) for yi high limit ryLi ) optimum limits change (as percentage of original range) for yi low limit Rj ) quarter of range for uj Rj0 ) base-case quarter of range for uj ujj ) optimized mean operating point for uj yji ) optimized mean operating point for yi ujj0 ) base-case mean operating point for uj yji0 ) base-case mean operating point for yi uholj ) half of limits for uj yholi ) half of limits for yi uj ) jth input variable (MV) for MPC yi ) ith output variable (CV) for MPC Greek Symbols ∆i ) range of operation for yi ∆uj ) change in uj operating point ∆yi ) change in yi operating point Ω ) yi state matrix Acknowledgment This work is supported in part by Natural Sciences and Engineering Research Council of Canada and Syncrude Canada, Ltd. Literature Cited (1) Agarwal, N.; Huang, B.; Tamayo, E. C. Bayesian approach for constraint analysis of MPC and industrial application. Proceedings of DYCOPS, 2007. (2) Bars, R.; Haber, R. Predictive control applied for linear and nonlinear plants. Technical Report, Budapest University of Technology and Economics, Budapest, Hungary, 2006. (3) Berghen, F. V. CONDORsA constrained, non-linear, derivativefree parallel optimizer for continuous, high computing load, noisy objective functions, Ph.D. Thesis, Free University of Brussels, Belgium, 2004. (4) Clarke, D. W. Application of generalized predictive control to industrial processes. IEEE Control Syst. Mag. 1998, (April), 49-55. (5) Dixon, L. C. W. Conjugate gradient algorithm quadratic termination without linear searches. IMA J. Appl. Math. 1975, 15, 9-18. (6) Hummel, H. K.; de Wit, G. B. C.; Maarleveld, A. Optimization of EB plant by constraint control. Hydrocarbon Process. 1991, (March), 6771. (7) Jacobs, D. A. H. A generalization of conjugate-gradient method to solve complex systems. IMA J. Numer. Anal. 1986, 6, 447-452. (8) Latour, P. H.; Sharpe, J. H.; Delaney, M. C. Estimating benefits from advanced control. ISA Trans. 1986, 25 (4), 13-21. (9) Madar, J.; Abonyi, J.; Szeifert, F. Interactive evolutionary computation in process engineering. Comput. Chem. Eng. 2005, 29 (7), 1591-1597. (10) Martin, G. D.; Turpin, L. E.; Cline, R. P. Estimating control function benefits. Hydrocarbon Process. 1991, (June), 68-73. (11) Maurath, P. R.; Mellichamp, D. A.; Seborg, D. E. Predictive control design for SISO system. Ind. Eng. Chem. Res. 1988, 27, 956-963. (12) Morari, M.; Lee, J. H. Model predictive control: past, present and future. Comput. Chem. Eng. 1999, 23, 667-682. (13) Prett, D. M.; Garcı´a, C. E. Fundamental Process Control; Butterworths Series in Chemical Engineering; Butterworths: Boston, MA, 1988. (14) Qin, J.; Badgwell, T. A survey of industrial model predictive control technology. Control Eng. Pract. 2003, 11, 733-764. (15) Rahim, M. A.; Shaibu, A. B. Economic selection of optimal target values. Process Control Quality 2000, 11 (5), 369-381. (16) Singh, P.; Seto, K. Analyzing APC performance. Chem. Eng. Progress 2002, 98 (8), 60-66. (17) Volk, U.; Kniese, D. W.; Hahn, R.; Haber, R.; Schmitz, U. Optimized multivariable predictive control of an industrial distillation column considering hard and soft constraints. Control Eng. Pract. 2005, 13, 913-927.

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8111 (18) Wismer, D. A.; Chattergy, R. Introduction to Non-linear OptimizationsA Problem SolVing Approach; North-Holland Series in System Science and Engineering; North-Holland: New York, 1978. (19) Xu, F.; Huang, B.; Akande, S. Performance assessment of model predictive control for variability and constraint tuning. Ind. Eng. Chem. Res. 2007, 46, 1208-1219. (20) Zhou, Y.; Forbes, J. F. Determining controller benefits via probabilistic optimization. Int. J. AdaptiVe Control Signal Process. 2003, 17, 553-568.

(21) Agarwal, N.; Huang, B.; Tamayo, E. C. Assessing MPC Performance. 2. Bayesian Approach for Constraint Tuning. Ind. Eng. Chem. Res. 2007, 46, 8112-8119.

ReceiVed for reView March 21, 2007 ReVised manuscript receiVed July 25, 2007 Accepted August 29, 2007 IE070417E