Robustness of a Class of Bias Update Controllers for Blending

Aug 15, 2002 - The aim of this paper is to study the robustness of a class of control strategies with bias update for blending operation. The class of...
0 downloads 0 Views 129KB Size
Download PDF
4786

Ind. Eng. Chem. Res. 2002, 41, 4786-4793

Robustness of a Class of Bias Update Controllers for Blending Systems Jose Alvarez-Ramirez,*,‡ America Morales,† and Rodolfo Suarez† Programa de Matematicas Aplicadas y Computacion, Instituto Mexicano del Petroleo, and Division de Ciencias Basicas e Ingenieria, Universidad Autonoma Metropolitana-Iztapalapa, Apartado Postal 55-534, Mexico D.F., 09340 Mexico

The aim of this paper is to study the robustness of a class of control strategies with bias update for blending operation. The class of blender control systems consists of a linear programming endowed with a bias updating scheme to reduce (nonlinear) model/plant mismatch within a real-time operation procedure. It is proven that these control strategies are robust against fluctuations in the stream properties and blending nonlinearities as long as deviations from the ideal mixing rules are sufficiently small. Discussions are concluded with a demonstration of the results on simulations of a simple gasoline blending control problem. 1. Introduction Blending operations are widely used in industrial plants including oil refineries, chemical plants, and cement plants. The operation mixes two or more streams with different properties to a given specification (e.g., temperatures, quality, etc.). For quality assurance, feedback control is often employed.1,2 Typically, a commercial blender control system consist of a linear programming based controller where cost function is optimized under certain flow rate availability and desired blended properties constraints.3-5 However, it can be difficult to achieve the desired blend characteristics if the compositions of the inlet streams vary in an unknown manner6 and if the blending rules contain uncertain/unknown nonlinearities.5 In fact, it has been shown7 that uncertainties can lead to a substantial loss in blender profitability given even small variation in the feedstock qualities. This is because uncertainties yield an inappropriate process model that cause the optimizing controller to drive the manipulated variables far from the process optimum. To address the blending operation control under model uncertainties, real-time optimization (RTO) procedures have been proposed.1,2,7 A real-time optimization controller (RTOC) consists of a set of integrated subsystems for (i) validating process measurements, (ii) updating the process model using these measurements, (iii) performing the optimization using the updated model, and (iv) validating the optimizer commands.7 Success of the optimizing controller in converging to the plant optimum depends on the ability of the model to accurately represent key aspects of the process behavior. Hence, the updating stage of the RTOC is of prime importance in practical situations. Bias updating (BU) schemes have been proposed to reduce the effects of model/plant mismatch in the performance of blending controllers. The idea behind the BU schemes is to estimate model uncertainties via a measurement-driven balancing of the constraints. More specifically, the plant * Corresponding author. E-mail: [email protected]. Fax: +52-5-8044900. † Instituto Mexicano del Petroleo. ‡ Universidad Autonoma Metropolitana-Iztapalapa.

outputs are compared with those predicted by the process model, the difference between the actual and the predicted plant outputs is used to update a bias or constant term appearing in the constraints of the model used by the RTO.7 Forbes and Marlin7 provided a criterion for determining sufficient conditions for the model-based optimization embedded within the RTOC, using BU scheme, is capable of finding the plant optimum despite model/plant mismatch. Recently, Singh et al.8 proposed an improved formulation for the gasoline blend optimization problem that incorporates a stochastic model of disturbances into the RTOC problem. By use of a nonlinear programming approach, improvements over linear programming approaches were reported. Such results motivates to explore more elaborated RTOC strategies to overcome the drawbacks in traditional (i.e., linear programming-based) blending strategies. Despite of the importance of blending processes in the process industry, to the best of our knowledge, a rigorous convergence and stability analysis of the iteration procedure induced by BU schemes has not been addressed. Since BU schemes are widely used by commercial blending systems,1,2 it would be desirable to dispose of a stability analysis to evaluate their viability and performance limits in the face of nonlinear mixing effects, uncertainties, etc. The aim of this paper is study the stability of a class of BU methods for RTO of blending processes. More specifically, we study the stability of linear programming plus BU scheme procedures to solve nonlinear blending problems. In this way, we constraint ourselves to the case where blending rules are given as an ideal (i.e., linear) mixing rule plus a nonlinear correction term due to high order interaction among chemical species. Regarding existing results in the literature7 for this class of RTOCs, our contribution can be summarized as follows: • Using feedback control and contracting map ideas, a sufficient stability condition in terms of the “size” of uncertainties (including nonlinearities or deviations from ideal mixing rules) is provided. From a practical viewpoint, such condition implies that convergence of the iteration procedure based on linear programming plus BU schemes will converge to the optimum of the

10.1021/ie0109455 CCC: $22.00 © 2002 American Chemical Society Published on Web 08/15/2002

Ind. Eng. Chem. Res., Vol. 41, No. 19, 2002 4787

actual nonlinear program as long as deviations from ideal mixing rule are sufficiently small. • It is shown that the BU scheme can be interpreted as a linear integral feedback controller acting on the “optimum regulation error”. This opens some possible extensions of BU schemes to the nonlinear programming case. The paper is organized as follows. Section 2 describes the class of blending control problem that will be considered in this paper. Section 3 describes RTOC and the corresponding BU scheme. Section 4 presents the stability of the BU schemes in the face of blending nonlinearities and inlet stream disturbances. Section 5 illustrates the results by means of a numerical simulation example. Section 6 closes the paper with some conclusions. Acronyms: BU, bias update; RTO, real-time optimization; RTOC, real-time optimizer controller. Notation: xT denotes the transpose of x; In denotes the n-dimensional identity matrix; the symbol 1n denotes the vector

second virial coefficient model retains the main nonlinear interactions of homogeneous mixtures. Summarizing, a mixture rule can be described either (a) as a pure nonlinear function pb ) fMR(u, p), or (b) as an ideal mixing rule corrected by a nonlinear function pb(u) ) n pb,id(u) + pb,nl(u), where pb,id ) ∑i)1 uipi is the ideal mixing rule and pb,nl is the (nonlinear) excess property or deviation from ideality. The ith stream is defined as Si ) {qi,min, qi,max, Pi, ci} where qi,min and qi,max are respectively the minimum and maximum flow rate availabilities, Pi ) (pi,1, ..., pi,m)T ∈ Rn is the vector of quality properties, and ci is the cost of the stream. Let Sb ) {qb,min, qb,max, Pb,min, Pb,max, cb} be a desired mixture. The blending problem consists of finding a suitable fraction composition u, such that n ∑i)1 uiSi f Sb at minimum cost cb. Let q ) (q1, ..., qn)T ∈ n n R be the vector of flow rates, and let qT ) ∑i)1 qi be the total flow rate. Without loss of generality, let us consider the case qb,min ) qb,max ) qT, where qT is the total required flow rate of the blended mixture. Then, ui ) qi/qT or u ) q/qT. Essentially, the blending control problem is a cost optimization problem on the manipulated variables or control input vector u ) (u1, ..., un)T. By assuming negligible storage dynamical effects, this program can be posed as follows:

the symbol 0n denotes the zero vector

mincTu u

(1)

subjected to the following constraints. • Minimum and maximum availability: 2. Control Problem Statement Consider a homogeneous mixture of n pure compon nents, with a fraction proportion u ) (u1, ..., un), ∑i)1 ui ) 1. The fraction ui’s can be given in any basis, e.g., volumetric, molar, etc. Let p ) (p1, ..., pn) be a given (quality) property vector. If the mixture was ideal, the property of the (blended) mixture pb is simply the average of the properties of the pure components at the same temperature and pressure, each weighted according to a given fraction base (e.g., molar, volumetric). If this were true in general for all the extensive properties, n uipi. In general, this equation one could write pb ) ∑i)1 is not correct since the mixing rule is given as a nonlinear functionality pb ) fMR(u, p). An example is given by the Reid vapor pressure model adopted in n gasoline blending process:9 fMR(u, p) ) [∑i)1 uipi1.25]0.8, where ui represents volumetric fraction and pi is the component Reid vapor pressure. In other cases, an alternative is to add a correction term ∆pb to the ideal n uipi + ∆pb, where the term ∆pb mixing rule: pb ) ∑i)1 is called the excess property or deviation from ideality. For ideal solutions, ∆pb ) 0. For real solutions, ∆pb is a nonlinear function of the composition u. Since an adequate equation of state for liquids is rarely known, modeling of ∆pb relies on empirical equations. Following ideas from the virial equation for gases, the following second virial coefficient equation has been proposed:10 n-1 n ∆pb ) ∑j)1 ∏i)j+1 uiujπij, where πij is the interaction coefficient for components i and j. The parameter πij depends mainly on temperature and pressure. The above virial-type equation can be generalized to account for third virial coefficients; however, few data are available for the third virial coefficients of pure materials, and still less is known about third cross coefficients. From a practical viewpoint, it is expected that the

umin e u e umax

(2)

u T 1n ) 1

(3)

• Total demand:

• Quality requirements:

Pb,min e Pb e Pb,max

(4)

where umin ) qmin/qT g 0 and umax ) qmax/qT > umin. Equations 1-4 become a nonlinear programming problem in the n-dimensional decision vector u with 2n + m + 2 constraints. In the sequel, it will be tacitly assumed that the nonlinear programming problem for (1)-(4) is feasible. Although the general control problem should consider nonlinear mixing rules of the form Pb ) FMR(u, p), where FMR(u, p) ) (fMR(u, p)1, ..., fMR(u, p)m), in this paper we restrict ourselves to the simpler case where the mixing rule can be written as Pb ) Pb,id(u) + Pb,nl(u). That is, we shall consider the case where the mixing rule is given as an ideal mixing rule Pb,id(u) corrected by a nonlinear term Pb,nl(u), which can be seen also as a deviation from ideality. This is not restrictive since even a pure nonlinear mixing rule fMR(u, p) can be expressed as fMR(u, p) ) pb,id(u) + pb,nl(u). This can be justified via the Taylor theorem11 when the mixing rule fMR(u, p) is an analytical function of its arguments. In such case, the nonlinear component pb,nl(u) is the residual part after the Taylor expansion. For instance, if fMR(u, p) ) n [∑i)1 uipi1.25]0.8 as in the Reid vapor pressure, then the deviation from ideality is the nonlinear function pb,nl(u) ) fMR(u, p) - pb,id(u).

4788

Ind. Eng. Chem. Res., Vol. 41, No. 19, 2002

In this way, the quality constraint (eq 4) can be rewritten as

Pb,min e Πu + N(u) e Pb,max

(5)

where Π ) (P1, ..., Pn) ∈ Rm×n is the matrix of stream quality properties and N(u) is a nonlinear function describing nonlinear mixing effects or deviations from ideality. Since most commercial optimizing controllers for blending systems are based on linear programming with a BU scheme,1-4 from a practical viewpoint, it would be desirable to dispose of an analysis framework to evaluate their robustness capabilities against nonlinearities and unmeasured stream qualities. 3. Optimizing Controllers with Bias Update Scheme In this section, we provide a brief description of RTOC with BU scheme7 from a modeling error compensation framework. Many commercial blender automation systems are based on a RTO layer that is usually based on linear programming with bias updating12 (see for instance the Internet page www.foxboro.com/nmr/prdctblndng.html). Such systems work by (i) taking a measurement of blended qualities, (ii) calculating biases (or modeling errors) as the difference between measured blended qualities and those predicted by the linear model, (iii) solving a linear programming problem using the updated biases to compute new blender feedstock recipes uk, and (iv) implementing the blend recipe. Basically, steps i-iv define a feedback control loop. On the basis of this information, the following assumptions will be considered. A.1. The nonlinearity N(u) is unknown. This means that the nonlinear function N(u) is not available for the solution of the optimization problem in (1)-(4). A.2. The blended qualities Pb(t) are available from online measurements, with a sampling period Ts > 0. A.3. Nominal values Π h of the stream qualities Π(t) are known. In many commercial blender automation systems, stream qualities Π(t) are measured, which introduces a feedforward path into the RTOC. In our analysis, we shall consider the worst case situation where Π(t) is not completely measured. This case arises when the measurement equipment is so expensive that measurements are made only at the blender output. Consider the blending process as a continuous-time dynamical system, with t denoting the time variable. Define the modeling error function B(q(t)) as

B(q) ) [Π - Π h ]u + N(u)

(6)

Then, the programming problem in (1)-(5) can be written as a RTO as follows:

min

reconstructed from measurements of the blend qualities Pb(t) (assumption A.2) and the availability of the stream molar proportions u(t). We know that Pb(t) ) Πu(t) + N(u(t)) for all t g 0. In this way, one can write (see eq h u(t) + B(u(t)) Hence, the modeling error B(u) 6) Pb(t) ) Π is trivially computed as

B(u(t)) ) Pb(t) - Π h u(t)

(8)

for all t g 0. This means that the modeling error signal B(u(t)) is uniformly observable for all t g 0. Physically, eq 8 means that the modeling error signal B(u(t)) closes the material balance associated with the mixing rule (2). In principle, the modeling error estimator (8) can be used in the RTOC (7) to compensate for uncertainties in the mixing rules. However, this procedure leads to an improper optimizing blending controller in the sense that u(t) is required to compute the modeling error signal B(u(t)) in eq 8, which is subsequently used in the RTOC (7) to compute u(t). That is, u(t) is required to compute itself. From the control theory viewpoint, this problem can be overcame with the use of a suitable estimation procedure. To this end, advantage is taken from the fact that the blending qualities Pb are measured with sampling period Ts > 0 (assumption A.2). The basic idea is to use delayed estimation to obtain a proper RTO. Let tk ) kTs, uk ) u(tk), Bk ) B(u(tk)), and so on. The corresponding discrete-time RTOC with modeling error estimation can be described in the following manner:

min cTuk, for k g 1 uk

subjected to umin e uk e umax uTk 1n ) 1 Pb, min e Π h uk + Bk e Pb,max

(9)

where Bk is computed by

h uk-1, with B0 Bk ) Pb,k-1 - Π

(10)

Posed as above, the discrete-time RTOC becomes a linear programming problem in the n-dimensional decision vector uk, which can be easily solved with traditional simplex methods. Besides, the modeling error estimator (10) is the so-called BU scheme,7 which becomes the adaptive part of the RTOC. Every value uk is pointwise optimum in the sense that it is the best value of the control vector u on the base of the available measurements Pb,k-1 and estimated information Bk. On the other hand, the aim of the modeling error estimator or BU scheme in (10) is to compensate the effects of the nonlinearity N(u).

u(t)

4. Stability Analysis

subjected to umin e u(t) e umax uT1n ) 1 Pb,min e Π h u(t) + B(u(t)) e Pb,max

(7)

The modeling error function B(u) represents the uncertain component of the nonlinear program (eq 7). Let us show that the modeling error trajectory B(u(t)) can be

The procedure using eqs 9 and 10 generates a sequence of vector values {uk ) 1, uk ) 2, ....}. Therefore, the sequence of linear programming problems (9) and (10) defines a dynamical system. From the practical viewpoint, a natural question arises; namely, under what conditions the sequence of linear programming problems (9) and (10) will converge to the optimum uop of the nonlinear programming problem (1)-(5)? In other words, how will the presence of the nonlinearity N(u)

Ind. Eng. Chem. Res., Vol. 41, No. 19, 2002 4789

affect the convergence of sequence {uk)1, uk)2, ...} to the optimum uop? In the following, we will provide a sufficient condition for convergence. The stability analysis of the general RTOC with BU scheme, even in the simple case considered in this paper where Pb ) Πu + N(u), is a very challenging theoretical problem because we add to the complexity of controlling an nonlinear uncertain system the fact that the number of active constraints can change at every iteration step. Hence, we are confronted with a bona fide hybrid system.13 Despite growing interest in hybrid systems the existing results are applicable to very simple and very specific (typically linear) dynamics.13 In this way, to handle more complex dynamics, like the application at hand, it is reasonable to try to develop instead a problem-specific analysis. For convenience in the stability analysis, let us rewrite the nonlinear programming problem (1)-(5) in a general form:

min cTu(t) u(t)

subjected to D1Tu(t) + b1 ) 0

(11)

D2u(t) + b2 + N(u(t)) e 0p where D1 ∈ Rn, D2 ∈ RR2p×n, b1 ∈ R, and b2 ∈ R2p, p ) 2(n + m), and N(u) is a nonlinear function corresponding to mixing nonlinearities. Let us use first-order sufficient conditions to convert the nonlinear programming problem (11) into an algebraic system.14 For each time t > 0, first-order conditions for problem (11) state that, if uop(t) is a regular point for all t g 0, uop(t) is a extremum relative to the problem (11) if there exists λ ∈ R and µ ∈ Rp, such that

µ g 0p µT[D2uop(t) + b2 + N(uop(t))] ) 0p

(12)

cT + λD1T + µT[D2 + ∇N(uop(t))] ) 0n The above system becomes a (n + p)-dimensional algebraic system with p constraints (µ g 0p) and (n + p) unknowns. A values µi > 0 defines an active inequalities. Then, for all active constraints, system (12) defines a set of r ) n + 1 algebraic equations with r unknowns, denoted by x ∈ Rr (vector x contains the control vector u and the Lagrange multiplier λ). According to the structure of system (12), and without losing generality, such set of equations can be written as

Ax(t) + M(x(t)) ) ysp

(13)

where y ∈ Rr is the measured output, which can correspond to some qualities Pb and some entries of the proportion vector u ∈ Rn, and ysp is a desired output value (set point) corresponding to active constraint values. Besides, A ∈ Rr×r is a given matrix depending on the stream qualities matrix Π, and M(x) ∈ Rn is a nonlinear function of the manipulated variables u ∈ Rn depending only on the blending nonlinearity N(u). For instance, if n ) 3 and m ) 1, and the constraint u1 e umax and Πu + N(u) g Pb,min are active, the correspond-

ing algebraic nonlinear system has

[

1 0 0 1 1 A) 1 p1,1 p2,1 p3,1

]

and M(u) ) (0, 0, N(u))T. Notice that the second row in A is related to the constraint u1 + u2 + u3 ) 1. It can be shown that det(A) ) p3,1 - p2,1, so that det(A) * 0 as long as the second and the third streams have not the same quality value. In this case, the measured output is y ) (u1, 1, pb,1)T and the corresponding set point value is ysp ) (umax, 1, pb,1,min)T. Notice that some entries of the vector M(u) can be zero, as those associated with availability and total demand constraints. In this way, for the sake of stability analysis only (and not to find the optimum), solving the nonlinear programming problem (11) is equivalent to solving the algebraic system (13) (see Luenberger14). To simplify the analysis, consider the following assumptions: A.4. xop is the unique solution of the algebraic system (13) for all t g 0; i.e., there is only one value xop such that Axop + M(xop) ) ysp. This means that uop is the unique extremum of the program (1)-(5). If uop is not the unique solution, the convergence results described below will only hold only locally (i.e., in a neighborhood of xop). A.5. If A h is an estimate of A, both A and A h are invertible matrices. This means that regularity of the value uop, relative to the problem (11), is assumed. In this way, according to assumptions A.1-A.5, the solution of the nonlinear programming problem (11) is equivalent to solve a set of algebraic equations (13) where the nonlinearity M(x) is unknown. Besides, the “output” y(t) ) Ax(t) + M(x(t)) is available for measurements. From a control theory framework, this problem can be seen as that of designing a robust (against nonlinear uncertainties) “control strategy” u(t) such that y(t) f ysp asymptotically. Let us describe the real-time control strategy with BU scheme corresponding to the system (13). As before, let B(x) ) [A - A h ]x + M(x) be the modeling error associated with the algebraic system (13), where A h is an invertible estimate of the matrix A. Then, one can write

A h x(t) + B(x(t)) ) ysp Analogously, the corresponding iteration is given by

xk ) A h -1(ysp - Bk)

(14)

h xk-1 Bk ) yk-1 - A

(15)

where

By virtue of assumption A.4, the aim of the iteration given by eqs 14 and 15 is to compute the control input trajectory xk such that yk f ysp asymptotically. Equations 14 and 15 can be combined to give

h -1(ysp - yk-1 + A h xk-1) ) xk-1 + A h -1(ysp - yk-1) xk ) A (16) The measurement yk introduces the correction A h -1(ysp - yk-1), by virtue of which xk is updated. It is then clear that, in the event the iteration (16) be stable,

4790

Ind. Eng. Chem. Res., Vol. 41, No. 19, 2002

yk f ysp asymptotically. By virtue of assumption A.4, this implies that xk f xop (hence, uk f uop) as desired. One has that yk ) Axk + M(xk), for all k g 1. Then, from (16) one can get

(17)

(c) It has been assumed that the stream qualities matrix Π is not measured. In practice, the matrix Π can be measured to update the corresponding estimate A h. In this way, eq 18 can be generalized to consider fluctuations in the matrix A and periodic quality measurements to give

Let ek ) xk - xop. Since ysp ) A h xop + B(xop) ) Axop + M(xop), eq 18 can be used to obtain

ek ) (In - A h k-1-1Ak-1)ek-1 A h k-1[M(ek-1 + xop) - M(xop)]

xk ) xk-1 + A h -1(ysp - Axk-1 - M(xk-1))

h -1A)ek-1 - A h -1[M(ek-1 + uop) - M(xop)] ek ) (In - A (18) By virtue of assumption A.4, eq 18 has a unique equilibrium point eeq ) 0r at least in a neighborhood of the origin. Notice that, if M(x) ) 0 and A h ) A, then ek ) 0r. That is, if there is no uncertainty, the BU scheme converges in one step. Hence, in the presence of “small” uncertainty, the BU scheme will converge asymptotically to the optimum value xop. This is stated in the following proposition, which is straightforward consequence of the Contraction Mapping Theorem.11 Proposition 1. Consider the BU scheme described before. Then, uk f uop if -1

h (A + JM(x))|| < 1 ||In - A for all x ∈ Rr, where JM(x) is the Jacobian matrix of the nonlinearity M(x). In regard to this result, the following comments are in order. (a) The main conclusion from proposition 1 is that, compared to the ideal (linear) mixing rule, nonlinearities should be sufficiently small in order to guarantee convergence of the linear programming plus BU scheme. Conversely, if the nonlinear mixing effects are large, instabilities are expected to be found. This fact limits seriously the applicability of linear programming plus BU schemes to solve blending control problems with significant nonlinear effects. Hence, alternative RTOC schemes should be explored. In this way, recent results by Singh et al.8 and Zhang et al.15 provide a suitable alternative. Specifically, Singh et al.8 have proposed a RTOC for gasoline blending that incorporates a stochastic model of disturbances within a nonlinear programming problem. Important improvements over linear programming approaches have been reported. Subsequently, Zhang et al.15 have studied some implementation and convergence characteristics of real-time optimization under parametric uncertainty. This optimization approach should be explored in commercial blending systems to overcome the limitation of linear programming approaches. (b) In some practical situations, the nonlinearity N(u) has a low contribution to the blended quality Pb. If In A h -1A is small enough, it is expected that JM(x) acts as a perturbation to the matrix A, so that ||In - A h -1(A + JM(x))|| < 1. An example of this case is the mixing rule n fMR(u, p) ) [∑i)1 uipiβ]1/β for small values of |β - 1|. In such case, the nonlinear component or departure from n n ideality is given by pb,nl(u) ) [∑i)1 uipiβ]1/β - ∑i)1 u i pi . Notice that pb,nl(u) f 0 as β f 1, so that, continuity arguments imply that |pb,nl(u)| is small for small values of |β - 1|. This would support the successful application of BU schemes to the RTO of blending controllers with marginal nonlinear mixing effects.3

h k is the where Ak is the value of A at time t ) kTs and A corresponding available estimate. Notice that in commercial blending systems, measurements of component quality variations are used to introduce a feedforward control scheme. In this way, the impact of component quality variations on the product quality is drastically reduced via a precompensation (feedforward) path. (d) Interestingly, the BU scheme can be seen as an integral action acting on the “regulation error” ysp - y. In fact, from eq 16 it is easy to get

h -1 xk ) x0 + A

k-1

∑ (ysp - yj)

j)1

where A h -1 plays the role of an integral control gain. In the continuous case, such integral action looks like x(t) h -1∫t0(ysp - y(t′)) dt′. In this way, the aim of the ) x0 + A BU scheme is to add integral control action into the RTOC to compensate for model/plat mistmaches. Given existing stability results of nonlinear control systems (see for instance Desoer and Lin16), it seems that BU scheme can be also used to enhance the stability and convergence properties of nonlinear programming approaches for real-time optimization. Although implementation results have been reported recently (see for instance Zhang et al.15), a rigorous analysis should be made to evaluate stability limitations. 5. Examples In this section, two examples are presented to illustrate the stability features of linear programming plus BU schemes. In the first example, a mixing model of the form pb(u) ) pb,id(u) + pb,nl(u) is taken. In the second example, a pure nonlinear mixing rule fMR(u, p) is considered, so that the “correction” part to an ideal mixing model is given by pb,nl(u) ) fMR(u, p) - pb,id(u). Although these examples are related to realistic physical situations and are used only to illustration purposes of the stability/instabilities features, it should be stressed that they do not describe a precise application situation. Example 1. Consider the blending of five feedstocks to obtain a gasoline with minimum (research) octane quality.4,9 Although it is a simple blending system, it retains the main characteristics of blending systems; namely, more than two streams with nonlinear mixing rules.5,8 The feedstocks are, for instance, (1) reformate, (2) LSR naphtha, (3) n-butane, (4) FCC gasoline, and (5) alkylate (Table 1). The parameters of the feedstock are the following: c ) (34.00, 26.00, 10.30, 31.30, 37.00)T in $/bbl, umin ) 05, and umax ) (1, 0.8, 0.38, 0.55, 0.6)T. The stream qualities have (10% fluctuation around the nominal value Π h ) (90.1, 75.7, 93.8, 82.9, 95.0). The nonlinear part of the mixing rule is modeled as a second n-1 n ∏i)j+1 uiujπij virial coefficient nonlinearity ∆pb ) ∑j)1 with π12 ) -6.0, π13 ) -8.25, π14 ) 8.25, π15 ) -6.0, π23 ) -7.8, π24 ) 9.0, π25 ) 11.0, π34 ) 11.0, π35 ) -

Ind. Eng. Chem. Res., Vol. 41, No. 19, 2002 4791 Table 1. Components and Properties component cost ($/bbl) umax octane RVP (psia)

(1) reformate 34.0 1.0 90.1 3.8

(2) LSR naphtha 26.0 0.8 75.7 12

(3) n-butane 10.3 0.38 93.8 138.0

(4), FCC gasoline 31.3 0.55 82.9 5.3

(5) alkylate 37.0 0.6 95.0 6.6

Figure 1. Performance of the RTOC for the sampling period Ts ) 1 min.

Figure 2. Performance of the RTOC for the following set of disturbances: a step set point change from Pb,max ) 90.0 to Pb,max ) 91.0 at t ) 25 min and a random (10% step change in the stream qualities matrix Π at t ) 50 min.

9.0, and π45 ) 11.5. These values are arbitrary, although they are of the order of the interaction coefficients reported by Morris et al.9 In this case, the contribution of the excess property is about (15% to the total blended octane. Figure 1 presents the performance of the RTOC for Π ) (92.1, 73.7, 91.8, 87.9, 92.0), Pb,min ) 85.0, Pb,max ) 90.0, and sampling period Ts ) 1 min. This figure clearly shows that the RTOC adjusted the manipulated variables to quickly and successfully satisfy product quality constraint. As it has been assumed, the number of active constraints is equal to the dimension of the decision vector u. In fact, it is noted that the five active constraints are u1 ) u1,min ) 0, u3 ) u3,min ) 0, u5 ) 5 u5,max ) 0.6, ∑j)1 uj ) 1, and Πu + N(u) ) Pb,max. This set of active constraints can be easily reduced to a twodimensional nonlinear algebraic system in the variables u2 and u4. Let us show that the RTOC is able to damp the effects of external disturbances. Figure 2 presents the performance of the RTOC for the following set of disturbances: a step set point change from Pb,max ) 90.0 to Pb,max ) 91.0 at t ) 25 min and a random (10% step change in the stream qualities matrix Π at t ) 50 min. It is observed that the RTOC is able to reject external disturbances while maintaining asymptotically the desired product qualities. This figure also illustrate the complexity of the blending system. The disturbance can change the structure of the active constraints. Every 5 uj ) 1, and Πu + N(u) ) time, the constraints ∑j)1

Pb,max are active. After t ) 50 min, the constraints u1 ) u1,min ) 0, u3 ) u3,min ) 0, and u5 ) u5,max ) 0.6 are active. After the disturbance at t ) 50 min, the constraint u4 ) u4,min ) 0 becomes active and the constraint u5 ) u5,max ) 0.6 becomes inactive. This shows that the optimization of blending system is a dynamical system with complex switching of active/ inactive constraints. The effects of the sampling period Ts are illustrated as follows. Figure 3 presents the evolution of the blending variables for two different values of Ts. By reducing the sampling period, the RTOC drives faster the product qualities Pb to the prescribed value. This means that an RTOC using BU scheme can be only successful if the sampling period is sufficiently small as compared to the characteristic dynamics of the disturbances in the stream qualities matrix Π. Finally, let us show that the RTOC based on a BU scheme can undergo instabilities. Figure 4 presents the RTOC performance when π45 ) 11.5 is changed to π45 ) -11.5. In this case, the BU scheme is unable to achieve the optimum value uop. Instabilities of the RTOC equipped with a BU scheme can be either to (i) the fail in the fulfillment of the condition of proposition 1, or (ii) to a more structural problem where the program has not a feasible solution. In case i, BU schemes cannot lead to the actual optimal solution, so that another model update schemes should be explored.

4792

Ind. Eng. Chem. Res., Vol. 41, No. 19, 2002

Figure 3. Effect of the sampling period Ts on the performance of the RTOC.

Figure 5. Dynamics of the linear programming plus BU scheme under a pure nonlinear mixing rule.

Figure 4. Instability in the RTOC due to a sudden change in the interaction parameter π45.

Example 2. Under the same cost and feedstock conditions as in Example 1, consider the case of the n nonlinear mixing rule fMR(u, p) ) [∑i)1 uipi1.25]0.8. In n this case, only the ideal mixing rule pb,id(u) ) ∑i)1 u i pi is available for RTOC, so that the nonlinear uncertainty n (or deviation from ideality) is pb, nl(u) ) [∑i)1 uipi1.25]0.8 n - ∑i)1uipi. Property vector p ) (3.8, 12.0, 138.0, 5.3, 6.6)T, which corresponds to Reid vapor pressure (psia).5,8

The following property thresholds were chosen: Pb,max ) 15.0 and Pb,min ) 12.5 for t < 10 min, and Pb,max ) 15.0 and Pb,min ) 10.0 for t > 10 min. Figure 5 presents the performance of the linear programming plus BU scheme. In the step change from Pb,min ) 12.5 to Pb,min ) 10.0 at t ) 10 min, there are switches in the set of active constraints. For instance, for t < 10, the quality constraint Pb,min ) 12.5 is active, while for t > 10 min, such constraint is inactive. This example shows that a linear programming plus BU scheme is able to handle pure nonlinear mixing rules as long as departures from n the ideal mixing rule pid(u) ) ∑i)1 uipi are sufficiently small. In this case, departures from pid(u) are not longer than 7%. This shows that the most important feature to ensure stability of linear programming plus BU schemes is the size of deviations from the ideality pid(u) and not the particular form of the mixing rule. Finally, the switching structure of active/inactive shows the hybrid nature of RTOC when BU schemes are used. This is an interesting feature from both theoretical and application viewpoints, and should be explored further. 6. Conclusions The robust stability of linear programming plus BU schemes for a class of blending control systems was addressed in this paper. The main result establishes that BU schemes converge to the optimum value as long as the nonlinear part (or deviation from ideality) of mixing rules and feedstock changes have a sufficiently

Ind. Eng. Chem. Res., Vol. 41, No. 19, 2002 4793

small contribution to the blended qualities. This stability property restricts seriously the applicability of linear programming-based blending systems in case where, e.g., feedstock variations are significant. The results in this paper should be seen as a first step toward a complete characterization of stability and convergence properties of RTOC for blending systems. More sophisticated approaches8,15 deserve a detailed stability analysis to evaluate their capabilities and limitations oriented to industrial applications. Studies on these issues will be addressed in a future work. Literature Cited (1) Serpemen, Y.; Wenzel, F. W.; Huubel, A. Blending technology key to making new gasolines. Oil Gas J. 1991, 89, 62. (2) White, J.; Hall, F. Gasoline blending optimization cuts use of expensive components. Oil Gas J. 1991, 89, 62. (3) Rigby, A.; Ladon, L. S.; Warren, A. D. The evolution of Texaco’s blending systems: from OMEGA to StarBlend. Interfaces 1995, 25, 64. (4) Diaz, A.; Barsamian, J. A. Meet changing fuel requirements with online blend optimization. Hydrocarbon Process. 1996, Feb, 71. (5) Chang, D.-M.; Yu, C.-C. Chien, I.-L. Coordinated control of blending systems. IEEE Trans. Control Syst. Technol. 1998, 6, 495. (6) Murakami, K.; Seborg, D. E. Constrained parameter estimation with application to blending operations. J. Process Control 2000, 10, 195. (7) Forbes, J. F.; Marlin, T. E. Model accuracy for economic optimizing controllers: the bias update case. Ind. Eng. Chem. Res. 1994, 33, 1919.

(8) Singh, A.; Forbes, J. F.; Vermeer, P.; Woo, S. S. Model-based real-time optimization of automotive gasoline blending operations. J. Process Control 2000, 10, 143. (9) Gary, J. H.; Handwerk, G. E. Petroleum Refining Technology and Economics, 3rd ed.; Marcer Dekker: New York, 1994. (10) Morris, W. E. Optimum blending gives best pool octane. Oil Gas J. 1986, 85, 63. (11) Kolmogorov, A. N.; Fomin, S. V. Elements of the Theory of Functional Analysis-I: Metric and Normed Spaces; Graylock Press: Rochester, NY, 1957. (12) Margoulas, K.; Marinos-Kouris, D.; Lygeros, A. Instructions are given for building a gasoline blending, LP. Oil Gas J. 1988, July 4, 32. (13) van der Schaft, A. J.; Schumacher, J. M. An Introduction to Hybrid Dynamical Systems; Springer Lectures Notes in Control and Information Sciences 251; Springer-Verlag: London, U.K., 2000; p 174. (14) Luenberger, D. G. Linear and Nonlinear Programming, Addison-Wesley: Reading, MA, 1984. (15) Zhang, Y.; Monder, D.; J. F. Forbes, Real-time optimization under parametric uncertainty: a probability constrained approach. J. Process Control 2002, 12, 373. (16) Desoer, Ch. A.; Lin, Ch.-A. Tracking and disturbance rejection of MIMO nonlinear systems with PI controller. IEEE Trans. Automat. Contr. 1985, 30, 861.

Received for review November 26, 2001 Revised manuscript received April 26, 2002 Accepted July 1, 2002 IE0109455