3944
Ind. Eng. Chem. Res. 2009, 48, 3944–3954
MPC Constraint AnalysissBayesian Approach via a Continuous-Valued Profit Function Seyi Akande, Biao Huang,* and Kwan Ho Lee Department of Chemical and Materials Engineering, UniVersity of Alberta, Edmonton, AB, Canada, T6G 2G6
Model predictive control (MPC) is one of the most studied modern control technologies. Among the various subjects investigated, controller performance assessment of MPC has received considerable attention in recent time. Various approaches and algorithms have been proposed for the assessment of MPCs. In this work, we propose a novel approach to MPC constraint analysis by considering the economic objective function as a continuous-valued function within a Bayesian probabilistic framework. The analysis involves inference of the effect of a decision to adjust the limits of the constrained variables with regards to the achievable profits (decision evaluation) as well as inference of constraint limits that should be adjusted so as to achieve a specified profit value (decision making). The benefits of this approach include a more generalized definition of quality variables, the development of a more rigorous formulation of the problem to address linear and quadratic objective functions and thereby obtaining closed form solutions, and maximum-likelihood location determination of the quality variables in the decision making process. The approach is illustrated with the use of simulations and a pilot-scale experiment. 1. Introduction The current academic and industrial interest in performance monitoring and assessment of industrial process controllers is motivated by the use of control systems to achieve goals related to quality, safety, and asset utilization.1-4 As a technology for asset-management, performance assessment helps industrial automation systems in maintaining high efficiency in production operations. Certain unexpected events can lead to disruptions in output and/or quality of processes and eventually, the process suffers in terms of loss of quality control of products, reduced machine efficiency, and ultimately, increased costs. A common example is the malfunctioning of sensors and/or actuators in an industrial process. It becomes necessary to quickly detect and correct such malfunctions and reduce associated variability. Consequently, it is beneficial to provide an automated procedure that would enable plant operators to effectively monitor and evaluate industrial processes and thereby increase profit margins. This is the main objective of controller performance monitoring/ assessment.3 A 0.1% improvement in an oilsands extraction/ upgrading process, for example, could result in savings of millions of dollars per year (given a yearly output of about 108 barrels) in a typical oil industry. The use of advanced process control (APC) strategies in industries has become crucial to achieving economic performance objectives. Among them, model predictive control (MPC) stands out as a key strategy for dealing with multivariable systems. “MPC refers to a class of algorithms that compute a sequence of manipulated variable adjustments in order to optimize the future behaviour of a plant.”5 MPC is a multivariable model-based controller and therefore it inherently considers the interaction between the process variables. One major advantage of MPC is its ability to handle input and output constraints. Tuning of an MPC is highly dependent on process knowledge and the particular control philosophy chosen for the process.5 Most approaches to tuning MPCs involve one or more of the following: (1) the penalty matrix, (2) the prediction horizon, and (3) the control horizon. All these parameters are important to the performance of the MPC, but the * To whom correspondence should be addressed. E-mail: biao.huang@ ualberta.ca. Tel.: +1-780-492-9016. Fax: +1-780-492-2881.
constraint limits and the variability of the process variables also play critical roles in controller performance. In general, however, there is no specific procedure or framework for choosing or tuning MPC constraints. Most approaches could be said to be largely subjective.6,7 In an MPC controller, the controlled variables (CVs) are often said to indicate the desired product qualities which are to be considered for optimization and the manipulated variables (MVs) are generally used to control and optimize the process within a constraint set. CV and MV constraints are often set conservatively in practice;8 therefore, the probability that adjusting constraints could provide more degrees of freedom for the controller and improving profit margins by doing so, typically exists.6,7 A method to estimate the effect of the constraints on the performance of an MPC will provide more insight and thus reduce conservativeness when setting up the constraints during commissioning or otherwise. In trying to achieve optimal MPC economic performance, the following questions arise: which constraints should be changed, and if constraint change is indeed possible, what profits can be obtained by doing so? These questions are answered by studying the relationship between variability and constraints within a Bayesian probabilistic framework. The details are described as follows:6,7 (A) Decision evaluation: To infer the expected return (or profit) that could be obtained, if certain decisions regarding constraint limits change are made. (B) Decision making: To obtain the combination of constraint changes of relevant variables, based on the maximum-aposteriori explanation, that will most likely achieve a userspecified target profit. This work considers the case where the expected return is obtained from a continuous-valued function. To make Bayesian inference and for the decision making/evaluation, the plant routine operating data, the plant steady state gains, which CVs and MVs are possibly allowed to change their limits, and the preference to make such changes, are required. Xu et al. (2007)9 formulated a linear matrix inequality based performance assessment (LMIPA) algorithm and that approach
10.1021/ie8011566 CCC: $40.75 2009 American Chemical Society Published on Web 03/16/2009
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3945
considers only the mean operating points of the quality variables in its analysis. It is more realistic, however, to consider the distribution (the mean and standard deviation or variance for a Gaussian distribution) of the quality variables, rather than just the mean operating points. This gives a more practical understanding of the process. Agarwal et al. (2007)6,7 developed an algorithm by considering the distribution of the quality variables and their effect on the potential profit to be obtained from the process. This previous work considered the controlled variables in the economic objective function as the quality variables, and made an assumption that the continuous-valued function for obtaining the expected returns could be discretized into zones of profit related to the process being above, within, or below the desired product specifications. While this was a commendable first step, the resolution from this kind of discretization is rather low, the lost information due to the discretization could have been valuable for the analysis, and the dimension of the problem can grow quickly, making its application in large-scale process impossible. In this work, we show that the continuous-valued function can be used without discretization and thus we avoid loss of information due to poor resolution, and the results are compact. The main contributions of this paper can be summarized as follows: (1) The use of the original continuous-valued profit function for the MPC constraint analysis. The derivation of the results is more mathematically elegant and, as a result, closed forms of solutions are obtained. (2) Maximum-likelihood location determination of the quality variable(s) in the decision making process. (3) Using the continuous-valued function allows for a natural extension to the consideration of statistical dependence between quality variables. (4) Detailed illustrations based on simulations and pilot-scale experiments are presented. The paper is organized as follows: Section 2 provides a brief background for the analysis by revisiting the Bayesian-based LMIPA problem formulation6,7,9 and introducing the continuousvalued profit function. This function is used to develop the basis for the decision evaluation and decision making formulations in section 3, in which cases of both dependence and independence between quality variables are considered. In section 4, pilot-scale experiments are used to evaluate the proposed algorithm, and finally, conclusions are made in section 5. 2. Discussion Bayesian-Based LMIPA Revisited. Linear matrix inequality based performance assessment (LMIPA) is a method for MPC performance assessment developed by Xu et al. (2007).9 For a particular process, the assessment of yield is calculated for various cases as follows:9,10 (1) Assessment of ideal yield. (2) Assessment of optimal yield without tuning the controller. (3) Assessment of improved yield by variability reduction. (4) Assessment of improved yield by relaxing constraints. (5) Constraint tuning for desired yield. These aspects have been extensively described and discussed in Xu et al. (2007)9 and Agarwal et al. (2007).6,7 In this work, we will mainly consider the effect of constraint adjustment (or tuning) of the MPC on the potential yield of the plant and make inference on constraint analysis within a Bayesian statistics framework. Consider an MPC application with n inputs and m outputs, with K as the steady state gain matrix and (yji0;ujj0) as the current mean (or base case) operating points for the jth input and ith
output. This state of the process based on the current steady state operation, prior to the profit analysis, is called the base case operation. Also let the number of controlled variables (CVs) and the number of manipulated variables (MVs) for which constraint limit changes are allowed (by relaxation) be a and b, respectively. This makes a total of N ) a + b variables for which we can possibly make limit changes. With yes and no as the options for applying the limit change to these N variables, there will be 2N combinations for applying the constraint changes. Each combination of limits change will have a specific optimal return and the optimum values of the mean operating point can be obtained through optimization, which will affect the MPC performance. Also, let (Lyi;Hyi) be the low and the high limits for yi, and (Luj;Huj) be the low and the high limits for uj, respectively. The optimization is carried out for each of the possible combinations, with the routine operating CV/MV data collected. The objective function is the economic objective function of the MPC controller and the constraints are the CV and MV constraint limits, while also taking into account the variability and the steady state gain relations.6 The quality variables are defined as all the variables (CVs and MVs) in the MPC economic profit function which have nonzero linear and/or quadratic coefficients. The optimization problem can be defined as a linear or a quadratic function of the following form: m
J)
∑
n
(εi yi + θ2i (yi - ηi)2) +
i)1
∑ (δ u
j j
+ γ2j (uj - νj)2)
j)1
(1) where, (yi,uj) and (ηi,νj) are the mean operating points and target values for the ith CV (yi) and jth MV (uj), respectively. εi, θi are the linear and quadratic coefficients for yi (outputs) and δj, γj are the linear and quadratic coefficients for uj (inputs), respectively. We also assume that within the constraints of the MVs and CVs, the derivative of the objective function is nonzero, that is, the optimum does not occur inside the constraint limits.6,7 With (yji0;ujj0) as the base case mean operating point, (yi;uj) as the optimum operating point, when the base case operating points are moved by (∆yi;∆uj), the equality and inequality constraints that need to be satisfied for economic optimization are shown in eqs 2-8:9 n
∆yi )
∑ (K
ij
× ∆uj)
(2)
j)1
yi ) jyi0 + ∆yi
(3)
uj ) ujj0 + ∆uj
(4)
Lyi e yi e Hyi
(5)
Luj e uj e Huj
(6)
Lyi - yholi × ryi e yi e Hyi + yholi × ryi
(7)
Luj - uholj × ruj e uj e Huj + uholj × ruj
(8)
where i ) 1, 2,..., m, and j ) 1, 2,..., n. We allow up to 5% constraint violation11,12 for the constrained variables in the desire to minimize the cost or maximize the profit potential. Kij represents the elements of the gain matrix of the process. Lyi and Hyi are the low and the high limits for yi, respectively. Lyj and Huj are the low and the high limits for uj, respectively,
3946 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009
ryi and ruj are the allowable change (percentage of the range) in the constraint limits for the process variables, in this paper 5% for process variables with changeable constraints and 0% for others; yholi, uholj are half of the constraint limits for yi and uj, respectively. Thus, the economic objective function for each constraint tuning can be specified as9 min J
The following results can be shown by using the expectation operation: E(F) )
[∑ ] q
∫
p(z1, ..., zq) z ,...,zq 1
Fi(zi) dz1, ..., dzq )
i)1
Fi(zi) (13)
i)1
(9)
jyi,ujj
[∑ ] q
E For the linear case,
subject to eqs 2-8. The optimum number of operating points obtained for all combinations of possible constraint changes is 2N. Superimposing the 2N optimum operating points with the base case variability and assuming the data to be Gaussian distributed, the probability distribution for the data under optimal operations can be obtained.6,7 The Bayesian network is subsequently created with N parent nodes, q child nodes (where q refers to the number of quality variables, that is, all variables in the economic objective function with nonzero linear and/or quadratic coefficients). For each of the cases of limit change considerations, the parent nodes have two states (yes, no) where yes implies change the limits and no implies do not change the limits. By changing the constraint limits, the quality variables are optimized to operate as close as possible to their maximum-return or minimum cost values. At the same time, other nonquality variables are also moved because of the interaction between the process variables. Each variable can take any value within the operating range and is continuous. The expected cost from the process can be calculated using eq 10. E(F) )
∫
z1,...,zq
p(z1, ..., zq) F(z1, ..., zq) dz1, ..., dzq
(10)
Fi(Zi) ) RiZi therefore,
[∑ ] q
E(F) ) E
Fi(zi) )
i)1
q
∑
q
RiE(zi) )
i)1
∑ R jz
i i
(14)
i)1
For the quadratic case, Fi(zi) ) Ri zi + βi2(zi - µi)2 therefore, q
E(F) )
∑ (R jz
i i
i)1
+ βi2(σzi2 + jzi2 - 2zji µi + µi2))
(15)
where zi represents the q quality variables, jzi the mean value of zi, Ri, and βi the linear and quadratic coefficients, σzi2 the variances of the quality variables and µi the target values. The above results (for the linear and quadratic cases) hold whether the quality variables are statistically correlated or uncorrelated. The above formulation can be applied to eq 1 to give the general form of the expected cost, as follows: m
where we introduce a variable z that represents all the quality variables (both CVs and MVs). E(F) is the expected cost, q represents the number of quality variables in the economic objective function and F(z1,..., zq) is an economic cost function which is continuous-valued in nature. The profit is estimated as the difference between the value of the base case cost function and the value of the cost function obtained after constraint limits have been adjusted. For this reason, the continuous-valued economic objective function is also alternatively referred to as the continuous-valued profit function in this paper, namely the two terms “cost” or “profit” are exchangeable, and we will not distinguish them unless where confusion arises. Continuous-Valued Cost Function. From the preliminary information on LMIPA previously provided in section 2 and eq 10, and with F(z1,..., zq) having the following additive form:6 q
F(z1, ..., zq) )
∑ F (z ) i
(11)
i
i)1
we can rewrite eq 10 as E(F) )
∫
z1,...,zq
q
p(z1, ..., zq)[
∑ F (z )] dz , ..., dz i
i
1
q
(12)
i)1
Previously this problem was solved using a discretization of the above function into finite zones of costs based on the probabilities of the quality variables (specifically the products) to be above, within, or below desired specifications.6 However, we propose a solution to this problem by considering the continuous-valued function directly and thereby derive a closed form of solutions.
E(F) )
∑ (ε jy
i i
i)1
+ θi2(σy2i + jyi2 - 2yjiηi + ηi2)) + n
∑ (δ uj
j j
j)1
+ γj2(σuj2 + ujj2 - 2ujjνj + νj2)) (16)
Bayesian Approach for MPC Constraint Analysis. For the cases of constraint change only, the covariance remains the same as that of the base case operation since for this case we assume the dynamic controller is not being tuned and the profit comes from simple adjustment of steady-state operating point. Optimizing the objective function subject to the new constraints provides the optimum operating points that are used as the new mean operating points for the benefit analysis. The parent nodes are defined as the variables that have the possibility of having their limits changed. The prior probability or the a priori for the parent nodes can be user defined or obtained from historical data. It indicates the preference to change or not to change the limits. For example, if a parent node has a priori of 0.7 for making a change to the limits, this means that the constraint limits for this variable has 70% tendency to change and 30% not to change. In this work, we use equal probability to change or not change the constraints (i.e., 50% to change or not to change the limits). In Agarwal et al. (2007),6,7 the selection of equal prior probabilities has been shown to give the maximumlikelihood solution. Figure 1 shows the DAGs (directed acyclic graphs) of two Bayesian networks. Figure 1a shows the structure for two parent nodes and two child nodes, while Figure 1b shows the structure for a network with multiple parent and child nodes and how they affect the obtainable profit. We refer to the work by Agarwal et al. (2007)6 for the procedure for building Bayesian networks for this type of analysis.
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3947
Figure 1. Illustrative directed acyclic graphs (DAG) of Bayesian networks.
The child nodes are defined as the quality variables, that is, the variables which affect the cost function. These are the variables with nonzero linear and/or quadratic coefficients in the MPC economic objective function. The conditional probability distributions of the child nodes are specified on the basis of the optimum operating points obtained from the optimization step and the base-case standard deviations from the data. The Bayesian analysis is based on the probability distribution of the data. For this work, the distribution is assumed to be Gaussian. The Bayesian analysis for decision evaluation is carried out using the Bayesian Network toolbox (BNT) developed by Kevin Murphy13-15 and the decision making is directly solved by optimization. Decision Evaluation. This refers to inferring the achievable profit, if certain decisions regarding limits change are made. For decision evaluation, the decision whether to change the limits is provided. This is equivalent to the evidence for the Bayesian network conditional on which, the probabilities of the locations of the quality variables are then estimated. Thus the cost or profit can be evaluated based on the relation specified in eq 15. Decision Making. This refers to obtaining the maximum-aposteriori explanation that will help to achieve a target value of the profit or cost. For decision making purposes the target profit or cost is provided and the corresponding optimum values for the quality variables affecting the cost function are estimated. On the basis of the optimum values for the quality variables obtained, the algorithm suggests a combination of variables for which changing their constraint limits will most likely achieve the specified profit. It is not necessary however that, given a profit or cost value, there is a unique set of quality CVs/MVs location. Multiple sets of locations may exist to yield the same desired profit or cost. This is shown in Figure 2 and will be elaborated shortly. Based on the probability distribution assumption for the base case variability, we use an optimization approach to find the locations of the quality variables (i.e., their values) that will have the maximum probability of occurring while also yielding the desired profit. As a result, the maximum-a-posteriori estimate of the states of the parent nodes (i.e., change the limits or not) can then be obtained. Based on the formulas for the continuous-valued cost function, we proceed to show how the decision making can be
Figure 2. Schematic illustrating choice of combination of quality variables to achieve a target return.
carried out. First consider the case when the quality variables are independent: q
p(z1, ..., zq) )
∏ p(z )
(17)
i
i)1
Thus: q
q
∏ p(z ) ) ∏ i
i)1
i)1
[
1
√2πσi2
(
exp -
(zi - jzi)2 2σi2
)]
(18)
where Gaussian probability distributions of the variables have been assumed for each of the q quality variables. jzi is a set of optimal operating points for a specific constraint tuning, obtained from the decision evaluation step. To obtain the most likely locations of the quality variables that will give the desired return, we maximize the probability
3948 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009
density function: q
∏
p ) max
zi,...,zq
i)1
using log-likelihood: q
max
zi,...,zq
∏ p(z ) ) max i
zi,...,zq
i)1
[
1
√2πσi
2
{ ( [√ q
log
∏ i)1
thereby obtaining q
max
zi,...,zq
(
exp -
1
zi,...,zq
2σi
(
2πσi2
)]
(19)
(zi - jzi)
2
2σi2
)])}
(20)
( (
∏ p(z ) ) max ∑ i
2
exp -
q
i)1
(zi - jzi)2
Thus:
-
(zi - jzi)2 2σi2
i)1
))
p(z1) p(z2) )
max
zi,...,zq
( ( q
∏ p(z ) ) min ∑ i
i)1
zi,...,zq
i)1
(zi - jzi)2 2
2σi
i)1
+ βi2(σzi2 + zi2 - 2ziµi + µi2)) ) RT
))
(21)
max {p(z1, ..., zq)1, ..., p(z1, ..., zq)2N}
1
×
2πσ22
(
(z2 - jz2)2 2σ22
)
(25)
(26)
max p(z1) p(z2) ) max{log(p(z1) p(z2))}
(27)
This gives
[(
z1,z2
(22)
(23)
We will illustrate the above formulations using an example where q ) 2, and z1 and z2 are independent, that is, p(z1, z2) ) p(z1) p(z2)
2σ1
)√
p ) max p(z1, z2)
max -
and also subject to the hard constraints on the quality variables. As illustrated in Figure 2, all points falling on the circle in the z1 and z2 plane satisfy the specified cost (alternatively return) constraints. Thus the solution to this problem is nonunique and the above optimization is justified. In Figure 2 the locations A, B, C, and D are examples of these points. These locations could all satisfy the specified target cost (or return) as shown in eq 22. The probability associated with all the possible locations is determined by the joint distribution of z1 and z2 and the optimization is carried out to find the location that has the highest probability of occurring. Equation 21 is an optimization problem with a quadratic objective function, a quadratic equality constraint (eq 22), and inequality constraints (eqs 5 and 6). This optimization problem can be solved using YALMIP,16 an open-source optimization toolbox for MATLAB.17 YALMIP is an interface for SeDuMi, SDPT3, and many other optimization algorithms. It is used for rapid prototyping of optimization problems and although it initially focused on semidefinite optimization problems, its scope has been significantly extended to include various other optimization problems. Other optimization toolboxes used are the MATLAB optimization toolbox,18 SeDuMi toolbox,19 and SDPT3,20 which are also useful for semidefinite programming and optimization over symmetric cones. The optimization is performed 2N times, once for each of the combinations of limits change, where N is the number of variables chosen for limits change, and a set of optimum values are obtained in each case. This results in 2N sets of optimum values. The decision making process is then completed by using eq 23 to determine which set of optimum values out of the 2N cases has the highest probability of occurring. (z1, ..., zq)*
2
Using log-likelihood:
q
i i
(z1 - jz1)2
As stated above, to obtain the most likely locations of the quality variables that will give the desired cost, we maximize the probability density function:
subject to the specified cost constraint:
∑ (R z
√2πσ12
(
exp -
exp -
or equivalently q
1
(24)
or equivalently max z1,z2
[(
(z1 - jz1)2 2σ12
(z1 - jz1)2 2σ21
) ( -
) ( +
(z2 - jz2)2 2σ22
(z2 - jz2)2 2σ22
)]
)]
(28)
subject to 2
∑ (R z
i i
i)1
+ βi2(σzi2 + zi2 - 2ziµi + µi2)) ) RT
(29)
and subject to eq 5 and 6. The values of jz1 and jz2 in eq 28 are obtained from the optimization of operating points previously described. In this illustration, since we have 22 ) 4 combinations of limits change, we will have four pairs of mean values (zj1 and jz2) for each case of constraints change. For each pair of means, the optimization described in eq 28 and its associated constraints will give the corresponding optimum values of z1 and z2. Overall, for all the cases considered, four pairs of optimum values for z1 and z2 are obtained (i.e., (z1,z2)1, (z1,z2)2, (z1,z2)3 and (z1,z2)4). Finally, using eq 30, the pair of optimum values (i.e., (z1,z2)1, (z1,z2)2, (z1,z2)3 or (z1,z2)4) that have the maximum probability of occurring is chosen, that is, max {p((z1, z2)1), p((z1, z2)2), p((z1, z2)3), p((z1, z2)4)} (30)
(z1,z2)*
where (z1,z2)* is the pair of optimum values of all four cases, with the maximum probability of occurrence. This solution is also the maximum a posteriori with the specified cost value given as the evidence. When the quality variables are dependent, we consider the joint probability distribution of the variables and the following formulation applies: p(z1, ..., zq) )
1 1 j )TΣ-1(Z - Z j ) (31) exp - (Z - Z 1/2 2 2π|Σ|
(
)
j is where Z is a vector of quality variables for optimization, Z the vector of corresponding mean values, and Σ is the corresponding covariance matrix.
( )
z1 - jz1 j) ) l (Z - Z zq - jzq
and Σ )
[
σ12 · · · σ1q σ22 · · · σ2n
σ21 σ21
· ·· l l σq1 σq2 · · ·
l σ2q
]
(32)
where Σ can be estimated from the base case operation data. The objective function for this case can also be formulated as shown previously using log-likelihood: 1 j )TΣ-1(Z - Z j) max - (Z - Z z1, · · · ,zq 2
(
)
or equivalently min
z1, · · · ,zq
( 21 (Z - Zj) Σ
T -1
)
j) (Z - Z
(33)
subject to eq 5 and 6 and subject to the specified cost: q
∑ (R z
i i
i)1
+ βi2(σzi2 + jzi2 - 2zji µi + µi2)) ) RT
(34)
Subsequently, eq 33 is used 2N times for the possible combinations of limits change; eq 23 is finally used to complete the decision making process. Choice of Variables for Limits Change. The fundamental assumption of the LMIPA algorithm is that the optimum values of the quality variables lie close to their constraint limits. Relaxing the constraints of other process variables that influence the quality variables will result in moving the quality variables closer to their limits and thereby obtain more profit. These process variables are believed to be under hard constraints that require some relaxation. It may be determined by the engineers which CV/MV may allow its constraint limits to be changed. This choice is either by the understanding of the process or by an intention to test how sensitive the economic performance to this CV/MV. The latter does not really change the constraint limits but tests sensitivity through the algorithm. The choice of variables for limits change could also be made based on existence of relations between quality variables and other variables according to the steady state gain matrix. They could also be chosen by identifying constrained variables via visual inspection and making such decisions from experience or using heuristics. Illustrative Example of a Bayesian Network for MPC Constraint Analysis. The following example illustrates the procedure for building a Bayesian network to carry out the procedures for MPC constraint analysis as described previously. We consider a 2 input-2 output multivariable system with a steady state gain matrix as follows: K)
[
k11 k12 k21 k22
]
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3949 Table 1. Linear (L) and Quadratic (Q) Coefficients of Variables and Limits Change Specifications variable
L coeff
Q coeff
limits change
CV1 CV2 MV1 MV2
0 RCV2 RMV1 0
0 0 0 0
yes no no yes
Table 2. Specifications of Parent and Child Nodes parent nodes
CV1*
MV2*
states child nodes states
(yes, no) CV2 (µ,σ2)
(yes, no) MV1 (µ,σ2)
Since we have specified two variables for limits change, the Bayesian network will have two parent nodes (CV1* and MV2*). Also, since we have two variables with nonzero linear and quadratic coefficients, we have two child nodes for the network (CV2 and MV1). Table 2 shows the specifications for the parent and child nodes. We assume that the prior probability of each of the parent nodes is 0.5, implying that there is an equal probability that we will make a change in the limits or leave the limits unchanged. For this system, we will have 22 ) 4 combinations of limits change for which we can carry out optimizations to obtain the optimum values of the quality variables based on our choice of limits change. The conditional probability distribution (CPD) of the child nodes, being continuous, is defined by mean values and variances for each node. The mean values of these child nodes are the optimum values obtained from the optimization of the operating points and the variances or covariances are the base case variances or covariances of the quality variables, obtained from data. Since CV2 and MV1 are identified as quality variables and their CPD is specified as discussed above, with the conditional probability table (CPT) for CV1* and MV2* as parent nodes, each being specified probability of 0.5 each, we have a Bayesian network as shown in Figure 3, which can subsequently be used for decision making and decision evaluation. Readers are referred to Agarwal et al.6,7 for illustrative examples. A Case Study of a Binary Distillation Column. The binary distillation column is a system with 10 CVs (outputs, y) and 4 MVs (inputs, u) (see Table 3). A detailed description of this system is available in the literature.10,21 A study was carried out by Agarwal et al. (2007)6,7 for this system based on the control objective described in Volk et al. (2005).21 The linear quadratic objective function coefficients required to carry out the analysis are provided in Table 4. This case study is included here to compare the proposed approach with the previous work. The steady state gain matrix10 for the process is shown in eq 36.
(35)
The input variable is also referred to as the manipulated variable (MV) and the output variable is also known as the controlled variable (CV). In the following discussion we will differentiate between the process variable (CV or MV) and the binary variable (the decision whether to change limits) by labeling them as CV or MV (process variables) and CV* or MV* (binary variables), respectively. Let CV2 and MV1 be quality variables (i.e., variables in the MPC economic objective functions with nonzero linear and/or quadratic coefficients) while CV1 and MV2 are constraint variables, see Table 1.
Figure 3. Bayesian network for illustrative example.
[ ]
3950 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009
0.9500 -0.0469 -0.0600 0 -0.0170 K) 0 0 0 0 0
0 0 0 -4.9852 0.1341 0.3324 -6.3750 0.1715 0.4250 -25.000 0 0 -8.8125 0.1930 0.5100 0.3750 0 0 0 0.2680 0 0 -0.0700 0 0 0 0.1700 0 0 0.6650
Table 3. Description of Controlled and Manipulated Variables Controlled Variables
(36)
As an illustration, MV1 and MV3 are chosen for limits change since they have some influence on one or both quality variables, CV2 and CV8. MV1 has some influence on CV2 but no influence on CV8, while MV3 has influence on both CV2 and CV8. CV5 was chosen because it is one of the CVs that have interactions with all the MVs, and CV9 is chosen without any particular basis. Since CV9 does not have any interaction with MV1 and MV3, we will subsequently show that such an apparently irrelevant choice does not result in a problem, although by itself, it does not contribute to the profits achieved. The implication of these choices will be discussed later. Thus the variables that are chosen for limits changes are MV1, MV3, CV5, and CV9. The objective function can be defined using the information in Table 4 and substituting them into eq 1. εi and θi are the first 10 values of the linear and quadratic rows in Table 4, respectively. Likewise, δj and γj are the last four values of the linear and quadratic rows in Table 4, respectively. 3. Results Linear Objective Function. Decision Evaluation. The potential profit that can be obtained from the change of limits is calculated as follows: profit ) base case cost - cost from limits change (37) where the base case cost is the cost from operating at current mean values. In this case study, base case cost ) -69.6293. The cost from limits change is the cost from operating at optimum values due to limits change. Equation 37 is used to obtain Table 5, which shows the profits calculated from the proposed Bayesian-based approach (continuous), the profits calculated from the discrete Bayesian-based approach (discrete), and the profits based on LMIPA (RLMIPA). The notation in Table 5 indicates the combination of limits change for each case. For example, “YNYN” indicates “yes, no, yes, no” for each of CV5, CV9, MV1, and MV3, respectively. From Equations 14 and 37, the profits are calculated, and the results shown in the first three columns of Table 5 indicate that the continuous Bayesian approach produce the same results as those of the deterministic LMIPA for the linear case. For the quadratic case, we would expect a difference since the formulation in eq 15 includes the variance (or covariance) term(s) as well as the mean values. From Table 5, comparing case 1 and case 16 (and the cases in between) of the “continuous” results column, we observe that there is an overall trend between case 1 (where only one variable has its constraint limits changed) and case 16 (where all four process variables have their constraint limits changed). This trend indicates that the more the number of process variables that are chosen for constraint change, the more the profits that can be obtained from
reflux flow light petrol FBP top PCT pressure valve position column bottom PCT column pressure feed temperature reboiler furnace duty duty bypass valve position Manipulated Variables reflux flow controller SP column pressure controller SP feed temp control valve position duty valve position Table 4. Linear and Quadratic Optimization CoefficientssLinear Case CV number 1
2
MV number
3 4 5 6 7
8
9 10
linear 0 0.2364 0 0 0 0 0 0.1714 0 quadratic 0 0 0 0 0 0 0 0 0
0 0
1
2
3
4
0 0
0 0
0 0
0 0
Table 5. Comparison of Decision Evaluation Results for Discrete and Continuous Algorithmsa case
R-LMIPA
continuous
discrete
1 (NNNN) 2 (YNNN) 3 (NYNN) 4 (YYNN) 5 (NNYN) 6 (YNYN) 7 (NYNN) 8 (NYYN) 9 (YYYN) 10 (NNNY) 11 (NYNY) 12 (YYNY) 13 (NNYY) 14 (YNYY) 15 (NYYY) 16 (YYYY)
1.7513 2.0257 1.7513 2.0303 1.9212 2.2025 1.9212 2.2071 1.7613 2.0309 1.7613 2.0364 1.9312 2.2086 1.9312 2.2132
1.7513 2.0257 1.7513 2.0303 1.9212 2.2025 1.9212 2.2071 1.7613 2.0309 1.7613 2.0364 1.9312 2.2086 1.9312 2.2132
2.2809 2.2809 2.2809 2.2809 2.2809 2.2809 2.2809 2.2809 2.5511 2.5117 2.5511 2.5511 2.5511 2.5511 2.5511 2.5511
a R-LMIPA: Profits based on LMIPA. Continuous: Profits calculated from the proposed Bayesian-based approach. Discrete: Profits calculated from the discrete Bayesian-based approach. The notation indicates the combination of limits change for each case. For example, “YNYN” indicates ”yes, no, yes, no” for each of CV5, CV9, MV1, and MV3, respectively.
Table 6. Decision Making Results for Linear Case: Locations of Quality Variables quality variable
base case mean
operating point (decision making)
CV2 CV8
67.9986 -500.025
60.5912 -500.0181
Table 7. Quality Variables and Specified Values
CV2 CV5 CV6 CV8
description
mean
min
max
light petrol FBP column bottom PCT column pressure reboiler furnace duty
67.9986 85.81 0.16 -500.03
40.35 81.15 0.21 -993.73
69.75 118.85 0.29 -6.27
the process. This observation would intuitively be expected because we get more degrees of freedom as more variables have their limits relaxed. The overall trend in Table 5 for the “continuous” results column shows that the values of the profits for the even-numbered cases are generally larger than the values of the profits for the odd-
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3951 Table 8. Linear and Quadratic Optimization Coefficients and Quadratic Targets CV number 1
2
3 4
5
6 7
MV number 8
9 10 1
linear 0 0.2364 0 0 1 1.5 0 0.1714 0 0 quadratic 0 0.5 0 0 2 0 0 0.1 0 0 target - 61.23 - - 80.65 - - -499.58 - -
0 0 -
2
3
4
0 0 -
0 0 -
0 0 -
Table 9. Decision Evaluation: Profits for Quadratic Case case
R-LMIPA
continuous
1 2 3 4
4014.04 4983.54 4014.72 5075.68
4013.32 4982.83 4014.01 5074.97
Table 10. Decision Making Results for Quadratic Case: Locations of Quality Variables quality variable
base case mean
operating point (decision making)
CV2 CV5 CV6 CV8
67.97 86.32 0.157 -500.11
61.11 79.93 0.209 -497.49
numbered cases. While this may not always be so, the explanation for this observation is due to the fact that some variables contribute more to the total profit when their constraints are adjusted, than others, and whenever these variables are included in the combination of variables having limits change as yes, they affect the value of the profit more. For example, by comparing the results from cases 1, 2, and 3 with respect to the combination of limits change, we note that changing the limits for CV5 produces additional profit and so whenever CV5 is included in the combination as a yes, additional profit would be expected. Whereas, in case 3, changing CV9 alone does not appear to produce any extra profit (same profit as in case 1). It is important to note, however, that the choice of CV9 as a candidate for limits change does not result in a loss; the profit just remains the same as it would have been if CV9 was not changed. In case 5, changing the limits of only MV1 produces some profit, but it is not up to the profit produced when the limits of both CV5 and CV9 are changed. Hence, although changing only CV9 does not produce any profit, we find that changing the limits for both CV5 and CV9 produces more profit than from only CV5 due to the interaction of CV5 and CV9. Therefore, although we find an overall trend indicating an increase in profits as more variables have constraint limits changed, the profit obtainable from each case depends on the combination of variables chosen and the individual contribution of each variable to the total profit. Comparing this proposed approach using the continuous-valued function with the approach proposed by Agarwal et al. (2007)6,7 for the discrete calculation and LMIPA by Xu et al. (2007)9 for the deterministic calculation, Table 5 shows that the LMIPA and the proposed continuous version of the Bayesian approach agree perfectly, whereas the results obtained from using the previous discrete approach6 does not clearly differentiate between the cases of limits change due to less resolution of the profit zones. Decision Making. Here, we seek to determine the process variables which should have their constraints changed on the basis of a user specified profit. Equation 21 together with associated constraints is used for the optimization to obtain the locations of the quality variables that will give the desired profit and then the algorithm determines the most probable combinations of variables for which changing their limits will give the desired results. It is important to note that the decision making algorithm involves the determination of a specific location of the quality
Figure 4. Schematic of multitank system configuration.
variable as evidence and not just the mean values as in the case of decision evaluation. Thus, the profit used in decision making is a point-wise profit while the profit used in decision evaluation is the averaged one. We illustrate a comparison between the decision making results for the discrete and continuous algorithms, with the following example. To achieve a target profit of 1.75 units, the discrete algorithm suggests a change in limits of both CV5 and CV9, while continuous algorithm suggests a change in limits for CV9 only. In fact it appears that the similar profit could be expected even without limits change according to Table 5. It is likely that this profit is achieved by simply moving the operating point closer to the profit limits without relaxing any constraint. For the given target profit value of 1.75 units, the calculated locations for the quality variables are shown in Table 6. The mean values, as obtained from the base case data, are compared with the suggested operating points obtained from the decision making analysis for the continuous algorithm. We observe that the operating point of CV2 is lowered, compared to its mean value, thereby creating room for potential profit. CV2 is the light petrol FBP, and it has its upper and lower limits as 69.75 and 40.35 °C, respectively, as shown in Table 7. Since the new suggested operating point of CV2 does not violate its constraints, but is still closer to the upper limits than the lower limit, this change is acceptable. Quadratic Objective Function. The linear and quadratic coefficients used for this case study are shown in Table 8. From eq 15, the expected cost can be calculated on the basis of the mean operating points and the variances of each quality variable, and the expected profit can be drawn. CV5 and MV3 are considered for limits change in this illustration. Decision Evaluation. There is a difference between the results obtained by the deterministic LMIPA and the Bayesian-based methods. This difference can be accounted for by the product of the quadratic coefficients of the objective function and the variance of the corresponding quality variable, as shown in eq 38: q
∑ (β
2
i
σi2)
i)1
The decision evaluation results are shown in Table 9.
(38)
()(
3952 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009
CV1
CV2
CV3
-0.04967 0 0 z - 0.8904 -0.04411z + 0.04201 0.03994z - 0.04013 0 ) 2 z - 1.764z + 0.7789 z2 - 1.841z + 0.847 0.03782z - 0.03758 -0.02095z + 0.0211 0.0139z - 0.01395 3 2 2 z - 1.897z + 0.9023z z - 1.778z + 0.79 z2 - 1.891z + 0.8908
Decision Making. The results using the decision making indicate that changing the limits of CV5 will achieve a desired profit of 4500 units. The suggested locations of the quality variables from decision making are shown in Table 10. These values are most likely achieved when the limits of CV5 have been changed by 5%. Comparing to Table 7, the new suggested operating points of the quality variables are well within the specified limits but are closer to the quadratic profit targets, thereby increasing the profits. 4. Experimental Evaluation Multitank Process: Process Description and Controller Design. The multitank system22 consists of three tanks (upper, middle, and lower tanks) arranged above each other, each equipped with drain valves. A pump fills the upper tank from a storage tank located at the bottom of the third tank as shown in Figure 4. Because of the vertical arrangement of the tanks, the liquid outflows all tanks due to gravity. The control objective is to maintain desired liquid levels in each of the tanks. Level sensors are mounted in the tanks to measure these levels. Each tank has a manual and an automatic valve. These valves can be used to stabilize the desired levels. The horizontal cross section of the first tank is constant while those of the second and third tanks are variable (prismatic and spherical, respectively) as shown in Figure 4. The system can be used as a single input single output (SISO), single input multiple output (SIMO), multiple input single output (MISO), or multiple input multiple output (MIMO) system, depending on the choice and/or combination of input and output variables chosen. For our purpose in this work, we consider the system as a 3 × 3 MIMO system with all three tank levels as the outputs (or controlled variables) and the corresponding automated valves are operated as the inputs (or manipulated variables). The pump flowrate is used as an unmeasured disturbance. A model was identified using the Matlab System Identification toolbox.23 The input signal was chosen to be a random binary sequence (RBS). The experimental conditions during the system identification are similar to those of the normal process operation such that a realistic model which is representative of the process can be obtained. The upper and lower limits of the input signals based on these considerations are then used in the design of the RBS signal. The following operating conditions were chosen: pump inflow, 34% (7.0353 × 10-5 m3/s); valve 1, 75% open; valve 2, 75% open; valve 3, 75% open. The valves only operate within the range of 50% to 100%. The estimated model is shown in eq 39. This model is subsequently used in the design of the MPC. The sampling time is 8 s. Linear Objective Function. The control objective in this case is to maintain desired liquid levels in the tanks, and we further assume our economic objective function to be the maximization of the liquid level in the third tank (i.e., max (100CV3) or min (-100CV3), which implies that CV3 has a linear coefficient value of -100. All other conditions are left unchanged unless further specified.
)( ) MV1
MV2
(39)
MV3
MPC Controller Parameters. The constraints are set as 0.5 e MVi e 1.0, and 0.05 e CVi e 0.25, where i ) 1,2,3. The input rate weights are given as [0.5, 0.5, 0.5], and the output weights are [1, 0.8, 0.9]. The control interval is 8 s, and the prediction and control horizons are specified as 15 and 2, respectively. By applying the designed MPC controller, and ensuring that all of the controlled variables and manipulated variables are running at steady state, the experimental data is collected and regarded as the base case operation data. The variables chosen for limits change are CV1 and MV3. This implies that we will have four (22 ) 4) possible combinations of limits change (see Table 11). The limits are allowed to be changed by 10%. Equation 40 shows the steady state gain matrix for this process. Since CV3 is the quality variable, it is obvious that from eq 40 that MV3 alone can be used to influence its operating point. Hence MV3 is chosen for limits change. Also, to further evaluate the choice of CVs being considered for limits change, CV1 is also chosen.
[
-0.45319 0 0 -0.01275 -0.35 0 K) -0.00943 0.02 -0.75
]
(40)
Decision Evaluation. Using the Bayesian method developed in this work, we obtain the decision evaluation shown in Table 12. The results are shown to be the same as the LMIPA results as previously discussed for this linear objective function case (R-case). The decision evaluation results obtained and the associated optimum operating values for each of the four cases were then used to conduct four validation experiments to ascertain the practicability of the suggested constraint tunings and determine if the potential profits can be actualized in reality. As shown in Figure 5, there is a good agreement between the calculated and experimentally achieved profits. Decision Making. When target profits of 1.7 and 2.5 are specified respectively, the algorithm determines the most probable combination of variables for limits change that will satisfy the target profits. The results indicate that changing the limits of CV1 will most likely achieve the target profit of 1.7 units, while changing the limits of MV3 will likely achieve a target profit of 2.5 units. The suggested operating points and constraint changes were then implemented in the experiment. Figure 6 shows that the decision making results have been experimentally evaluated. Note that there are some disagreements since the experimental profits are calculated as an averaged value while the decision making gives a point-wise profit. In addition, the maximum a posterior decision making only gives the most likely constraint tuning among a given set of tunings to achieve the desired profit target. Quadratic Objective Function. Following the above specifications for the linear case, we can also illustrate the case of a quadratic economic objective function by specifying that CV3 has a linear coefficient value of -100 and quadratic coefficient value of 50. We also choose MV1 and MV3 for limits change,
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3953 Table 11. Four Combinations of Limits Change parent nodes
case 1
case 2
case 3
case 4
CV1* MV3*
no no
yes no
no yes
yes yes
Table 12. Decision Evaluation: Profits for Pilot-Scale ExperimentsLinear Objective Function case
R-LMIPA
R-case
1 2 3 4
1.6837 1.7046 2.8795 2.8795
1.6837 1.7046 2.8795 2.8795 Figure 7. Experimental validation of decision evaluation resultssquadratic case.
Table 13. Decision Evaluation: Profits for Pilot-Scale ExperimentsQuadratic Objective Function case
R-LMIPA
R-case
1 2 3 4
23.6001 23.7438 34.4680 34.5975
23.8619 24.0056 34.7297 34.8592
indicating four possible cases of potential profits. The choice of MV1 and MV3 follows from the previous discussion based on the steady state gain matrix. Decision Evaluation. The results for decision evaluation are shown in Table 13 to compare the LMIPA and Bayesian-based approaches for the quadratic case. As shown in Figure 7, there is indeed profit to be obtained from constraint relaxation, even though the experimentally achieved profit is slightly less than what was calculated. Decision Making. When target profits of 25 and 31 are specified, respectively, the algorithm determines that changing the limits of CV1 and MV3 individually will most likely satisfy the respective target profits. As shown in Figure 8, this result is reasonable given that the experimental profits are calculated as an averaged value while the decision making gives a pointwise profit, as previously noted.
Figure 8. Experimental validation of decision making resultssquadratic case. Table 14. Comparing Features of Discrete and Continuous Algorithms discrete
continuous
compared with LMIPA quality variables
approximation mathematically accurate only CVs considered for all variables in economic computational simplicity objective function decision making approach discrete solution from direct solution through Bayesian software optimization
5. Conclusion
Figure 5. Experimental validation of decision evaluation resultsslinear objective function.
Figure 6. Experimental validation of decision making resultsslinear case.
A method based on the continuous-valued linear and quadratic economic objective functions has been developed. It provides the constraint analysis for the process variables of the MPC controller based on a Bayesian framework and the profit that can be achieved based on limits change (decision evaluation) is calculated. For decision making, the profits are specified and corresponding decisions concerning variables for limits change are inferred from the optimization. We have compared the results obtained using the approach proposed by Agarwal et al. (2007)6,7 and the results obtained using the proposed continuous-valued approach. The comparisons show that with this method, we can achieve calculations without loss of information due to discretization of profits, as in the case of the discrete algorithm. Table 14 summarizes the improvements of the continuous algorithm over the discrete version. The benefit of using the Bayesian approach, as shown in this paper, includes the ability to handle uncertainties by dealing with the probability distributions associated with the quality variables. The illustrations provided indicate the usefulness of this method for constraint analysis and tuning during process operation and maintenance of MPC controllers. It is also applicable for accessing controller performance and evaluating the opportunity for improvement and increasing profit margins. Finally, the utility of this method has been evaluated using pilotscale experiments.
3954 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009
Nomenclature σi ) standard deviation for yi σi0 ) base case standard deviation for yi E(F) ) expected cost J ) objective function K ) steady state process gain matrix Kij ) steady state gain for yi and uj uj ) jth input variable (MV) for MPC yi ) ith output variable (CV) for MPC Luj ) low limit for uj Lyi ) low limit for yi m ) number of output variables n ) number of input variables N ) number of constraint variables chosen for limits change P(¬A) ) probability for A not to be true P(A|B) ) probability for A to be true given B P(A) ) probability for A to be true P(B|A) ) likelihood of B to be true given A q ) number of quality variables ruj ) change in limits (as percentage of original range) for uj ryi ) change in limits (as percentage of original range) for yi ujj ) optimized mean operating point for uj jyi ) optimized mean operating point for yi ujj0 ) base case mean operating point for uj jyi0 ) base case mean operating point for yi uholj ) half of limits for uj yholi ) half of limits for yi
Acknowledgment This work is supported in part by Natural Sciences and Engineering Research Council of Canada. Literature Cited (1) Huang, B.; Shah, S. Performance Assessment of Control Loops: Theory and Applications; Springer-Verlag: London, 1999. (2) Bauer, M.; Craig, I. K. Economic assessment of advanced process controlsA survey and framework. J. Process Control 2007, 18, 2–8. (3) Jelali, M. An overview of control performance assessment technology and industrial applications. Control Eng. Pract. 2006, 14, 441–466. (4) Harris, T.; Yu, W. Analysis of multivariable controllers using degree of freedom data. Int. J. Adapt. Control Signal Process. 2003, 17, 569–588.
(5) Qin, S.; Badgwell, T. An overview of industrial model predictive control technology. Int. Conf. Chem. Process Control 5th 1997, 93, 232– 256. (6) Agarwal, N.; Huang, B.; Tamayo, E. Assessing model prediction control (MPC) performance. 2. Bayesian approach for constraint tuning. Ind. Eng. Chem. Res. 2007, 46, 8112–8119. (7) Agarwal, N.; Huang, B.; Tamayo, E. Assessing model prediction control (MPC) performance. 1. Probabilistic approach for constraint analysis. Ind. Eng. Chem. Res. 2007, 46, 8101–8111. (8) Singh, P.; Seto, K. Analyzing APC performance. Chem. Eng. Progress 2002, 98, 60–66. (9) 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. (10) Agarwal, N. M. Sc. Thesis. University of Alberta, Edmonton, Canada, 2007. (11) Latour, P.; Sharpe, J.; Delaney, M. Estimating benefits from advanced control. ISA Trans. 1986, 25, 13–21. (12) Martin, G.; Turpin, L.; Cline, R. Estimating control function benefits. Hydrocarbon Process. 1991, 68–73. (13) Murphy, K. P. How to use Bayes Net Toolbox. How to use Bayes Net Toolbox, 2004; http://www.cs.ubc.ca/∼murphyk/Software/BNT/bnt. html. (14) Murphy, K. P. The bayes net toolbox for MATLAB; Technical Report, 2001; www.cs.ubc.ca/∼murphyk/Papers/bnt.ps.pz. (15) Murphy, K. The Bayes Net Toolbox for Matlab. Comput. Sci. Stat. 2001, 33, 331–351. (16) Lofberg, J. YALMIP: A Toolbox for Modeling and Optimization in MATLAB; Proceedings of the CACSD Conference, Taipei, Taiwan, 2004. (17) MathWorks, Matlab User’s Guide; The MathWorks Inc.: Natick, MA, 2005. (18) MathWorks, Optimization Toolbox User’s Guide; The MathWorks Inc.: Natick, MA, 2007. (19) Sturm, J. Using SeDuMi 1.02, A Matlab toolbox for optimization over symmetric cones. Optimiz. Methods Software 1999, 11, 625–653. (20) Toh, K.; Todd, M.; Tutuncu, R. SDPT3° a Matlab software package for semidefinite programming. Optimiz. Methods Software 1999, 11, 545– 581. (21) 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. (22) Inteco, Multitank System User’s Manual; Inteco Ltd: Poland, 2006; www.inteco.com.pl. (23) Ljung, L. MATLAB System Identification Toolbox: User’s Guide; The Mathworks, Inc.: Nantucket, MA, 2007.
ReceiVed for reView July 27, 2008 ReVised manuscript receiVed February 5, 2009 Accepted February 6, 2009 IE8011566
