Controllability and Resiliency Analysis for Homogeneous Azeotropic Distillation Columns

Controllability and Resiliency Analysis for Homogeneous Azeotropic Distillation Columns


1MB Sizes 0 Downloads 156 Views

Copyright @ IFAC Advanced Control of Chemical Processes, Pisa, Italy, 2000


Boris M. Solovyev and Daniel R. Lewin

Wo/fson Department of Chemical Engineering Technion, Israel Institute of Technology, Haifa 32000, Israel

Abstract: Flowsheet controllability and resiliency analysis is a desirable capability in the early stages of process design. When implemented, it is commonly carried out using low-order approximations for the process dynamics, generated on the basis of preliminary flowsheet information. However, processing units that conditionally feature undesirable dynamics such as inverse response and overshoots require more complex models to account for these phenomena. Such models can be generated using a new method that relies on appropriate partitioning of the process into a number of component parts, from which an approximation of the linear dynamics is derived. In this paper, an azeotropic distillation column is analyzed using the new approach, leading to results that concur with and extend those obtained using rigorous dynamic simulation. Copyright © 2000 IFAC Keywords: Controllability and resiliency analysis; process dynamics; inverse response; homogeneous azeotropic distillation.


methods; (c) carrying out simultaneous design and control synthesis. Clearly, the unconstrained sequential approach does not guarantee a design that will meet its specifications when challenged by external disturbances, since the design may be devoid of the necessary degrees of freedom. The commonly accepted practice to subsequently test controllability by rigorous dynamic simulation should be seen only as a necessary verification step, and is too costly in terms of required engineering effort to implement as a screening method. In contrast, the full integration of design and control by simultaneous synthesis of a process and its control system is an attractive and truly optimal alternative, but one which cannot deal with the complexities of a large-scale process. Apart for anything else, many design decisions do not lend themselves to full automation, and will invariably involve significant intervention of the designer. Consequently, the inclusion of C&R screening

Traditionally, the chemical industry has been driven largely by the economy of scale, with increased profitability attained through increased production. Following the energy crisis of the early 70's, and with the resulting sharp increase in the cost of energy, profitability became associated with decreased production costs. Process synthesis since then has involved a balance between the search for designs that are attractive economically and those that can be operated safely and to meet specifications. This balance inevitably involves the interaction of design and control. As pointed out by Lewin (1999), approaches to integrated design adopted in practice are either: (a) sequential design in which the plant is synthesized first and controllability issues are then addressed; (b) sequential design constrained using short-cut controllability and resiliency (C&R) screening


The desired degree of separation is achieved by employing several possible operating conditions, depicted schematically in Fig. 2. The operation curve defines the limits of feasible operation as a function of the entrainer flowrate, E, and the column boil up, V. Andersen et al. (1991) defme three distinct operating regions, and point out that the operating line identifies possible input multiplicities when running in Regions Il and Ill, since for the same value of E, the desirable separation can be achieved by two distinct values of V.

methods as an integral part of any modem design approach seems to be the most practical alternative (Seider et al., 1999). The impact of the process structure on the plant's dynamic properties has been studied since the 1960's (e.g., GiIIiland et aI., 1964; Marlin, 1981). In cases where the units can be reliably described by perfect mixing-tank assumptions, the over-all process dynamics can be obtained directly from information available from steady-state material and energy balances. The dominant time constants are generated on the basis of designed flow rates and the results of preliminary equipment sizing (Weitz and Lewin, 1996; Lewin et al., 1996).

Table 1. Azeotropic Column Data Value 40 5 25 1 kmoUmin 0.9:0.1 10-7 : 10-4

Item Number of ideal trays Entrainer feed tray Azeotrope feed tray Azeotrope feed rate Azeotrope - Acetone:Heptane Entrainer - Acetone:Heptane Acetone top specification AcetonelHeptane bottom spec.

Since a more complex dynamic response can be generated by combining those of simple lags, a processing unit whose response trajectory is potential non-monotonic can be mode led by appropriate partitioning. In this approach, each partitioned segment of the unit is characterized by a single dominant driving force. The time constant associated with this driving force is obtained using information from the flowsheet material and energy balance. In this paper, this partitioning approach to dynamic mode ling is demonstrated according to its impact on the C&R analysis of a homogeneous azeotropic distillation column, which potentially features significant inverse response and overshoots, depending on the selected operating point.

0.998 0.005

Domain of improved separation Region I



A typical process design for the separation of a binary azeotrope of acetone and heptane using toluene as the entrainer consists of a sequence of two distillation columns as shown in Fig. 1. The entrainer recovery column presents no special difficulties, since it separates an almost ideal binary mixture, and so the analysis is focused on the azeotropic column (Andersen et al., 1991). Molar basis data for the column is summarized in Table 1.

Fig. 2. Operation curve as derived by Andersen et al. (1991) Inverse response occurs when a change in a process input affects an output in such a way that the initial trajectory moves in the opposite direction to its fmal, steady-state response. Mathematically, this occurs when the transfer function between the input and the output contains a right half plane (RHP) zero. It is most significant when the characteristic time constant of the zero is of comparable magnitude to that of the dominant time constant of the process. Andersen et al. (1991) indicate that significant inverse response with respect to E is to be expected in Region Il, whereas in Region Ill, this occurs with respect to both E and V. In Fig. 2, expected inverse response phenomena are indicated by dotted arrows aligned in the direction of the respective input variable. Thus, Region I is expected to be devoid of inverse response, and in practice, the column is usually operated in this region, since one would intuitively except improved controllability there. In contrast, the analysis of Jacobsen et al. (1991) suggests that


Fig. 1. Flowsheet for the separation of a binary azeotrope.


overloading the column with entrainer makes it no easier to control, which provides an incentive to run the column in the more profitable Region H.

3.3 Overall process dynamics.

The standard form is used: .rts) = ~(s) 1!(s) + Us) g(s)



where .rts) is a vector of process outputs, 1!(s) is a vector of process inputs, g(s) is a vector of process disturbances, and E..{s) and Us) are the process transfer matrices that relate output variables and manipulated variables and disturbances, respectively. It is noted that the linearizing output transformation suggested by Andersen et al. (1991) is used:

3.1 Mass and energy balances. A steady-state model for the azeotropic column was developed using the HYSYS commercial simulator from AEA Technology Software Engineering. The Wilson activity model was used for thermodynamic property predictions using data from Horsley (1952) for verification. The operation curve obtained using HYSYS is shown in Fig. 3, which is close to that provided by Andersen et al. (1991). In the following, approximate dynamic models are derived for three distinct operating points, as defined by the molar flow rates of E and V (in kmol/min), as indicated in Fig. 3: Region I (E = 1.05, V = 2.85), Region H (E = 0.85, V = 2.87), and Region III (E = 0.85 V= 3.33).


The manipulated variable and disturbance vectors are: 1!(s) = [D, L, L R, B, V, VR, E]T and g(s) = [E, F, ZF]T. It is noted that E appears alternatively in both 1! and g to study the impact of its allocation on the C&R of the column. The over-all dynamic models can be obtained by algebraic manipulation of the Trays and Utilities partition matrices following Lewin et al. (1996). The process outputs are elements of the vector LV, which is expressed in terms of VD :

3.2 Partitioning the process.

According to the proposed approach, the dominant driving forces need to be isolated in at least two partitions, to provide a potential mechanism for inverse response. For the azeotropic column, as illustrated in Fig. 4, the two partitions selected were the Trays and the Utilities sections. The former consists of the column tray internals, while the latter accounts for the reboiler and condenser. This partitioning is motivated by the observation that rapid transients can be expected from the tray hydrodynamics, while slower responses are expected from the interaction between the trays and the utility units. In Fig. 4, the vectors, TB, LV and VD consist of the following elements:


~ ~

I.s3.4 ...... ;::.~

TB = [Yt. Y2, Y3, Xt. X 2, x 3f , where Y, and X, are the component molar flowrates at the top and bottom, respectively. Note that the component ordering selected for i (acetone, heptane and toluene) depends on the operating region. LV = [log(1 - YD), L, Yt. Y2,V, Xt. X2, log(xR)]T, where L and V are the molar flowrates of reflux and boilup respectively, Yi and Xi are component molar ratios at top and bottom, respectively and YD and XR are molar ratios of acetone at top and acetonelheptane at bottoms, respectively. VD = [D, E, F, B, ZF, L R, VR]T, are the molar flow rates of top product, entrainer, azeotrope feed and bottom product respectively, the azeotrope feed concentration, and reflux and boilup ratios (LR =LlD and VR =VIB), respectively.

egion of improved separation

El 3.2


Region I









E, [-.-]


Fig. 3. Operation Curve obtained by steady-state simulation using HYSYS.

Fig. 4. Block diagram for partitioned process model.







only in Region Ill, with the response being minimum phase in Region 11 (and also in Region I, not shown in the figure). The model suggests weak inverse response with respect to E in regions II and Ill. These results are in good agreement with those of Andersen et al. (1991). It is noted that the overall dynamics are dominated by a characteristic time constant in the range 500-700 min. The large material recycle in the column is responsible for this phenomenon, which has been widely reported (e.g., Gilliland et al., 1964; Denn and Lavie, 1982; Kapoor et a/., 1986).

~ -10

"2 .~

ii -15





.S! -25

-30 0.97



I 1.01 F [kmole • min J



Fig. 5. Effect of F on 10g(xR) 3.4 Estimation of static gains.

These can be obtained from local linearization of the y - u curves, obtained by running parametric sensitivity analyses using the flowsheet simulator. For extremely small perturbations (±0.05% of full-scale input values), the y - u curves are almost linear. However, for larger perturbations (see Jacobsen et a/., 1991), they may be highly nonlinear, which will be reflected in significant uncertainty associated with the static gains. Fig.5 shows the y - u plot for 10g(xR)/F, showing the linear region. Jacobsen et a/. (1991) indicate that the column is resilient up to ± 15% disturbances in F, but Fig. 5 shows that this may be problematic for positive disturbances.




~ -1




Full model Reduced model.

-2 -;0-'






w .(radlmin)

(a) Region 11

3.5 Estimation of dynamic components.

Perfect mixing tank models are assumed for the Trays and Utilities partition matrix elements, with components of the and /i;j blocks in Eq. (3) being transfer functions of the form:




y(s) = _e_u{s) 13' + 1


~ 0

a: -2 In Eq. (4), the static gains, k, are computed using the procedure described in the previous section. The hydrodynamic time constants are computed as the ratio of the liquid holdup, M, to the liquid supply rate, 1. In contrast, the tray time delays, (Shinskey, 1984) and tray composition characteristic time (Moczek et a/., 1965) can be obtained using: 'tc,a = dM/dL.


Full model Reduced model




w .(radlmin)

(b) Region III Fig. 6. Bode plots for 10g(S)/V



3.6 Existence of Inverse Response.

Since all of the blocks and lk are proper, asymptotically stable, transfer function matrices, open-loop stability is determined by the roots of the characteristic equation:

To test the reliability of the partitioned model, Bode plots are generated to describe the effect of V and E on 10g(S), where S = y D /{1- y D)x R (degree of separation). As illustrated in Fig. 6, the model predicts severe inverse response with respect to V



The Nyquist plot in Fig. 7, generated for Region 11, accounts for the known uncertainty in the static gain coefficients and a range of feasible characteristic time constants. The absence of encirclements of the origin indicates open-loop stability. In the same way, the other two operating regions are shown to be open loop stable.

very sensitive to the uncertainty associated with the process nonlinearities. The example in Fig. 8 shows the effect of uncertainty in the static gain log(xR)/F on the static DC for the L/rB configuration, and indicates that steady-state disturbance rejection is not possible for larger disturbances in F. This range of uncertainty is well within the range associated with perturbations considered by Jacobsen et al. (1991). 1.6

.:: 0.8

r ~


0 .6

.... ~ .... 0.4


1. 2





O~~~==============~==~ -0.25 0 0.25 0.5 0.75 1.25 Re(ll-ull ·cllD

0.8 0.6

-~3000~---~~~OO~--2000~~--1~500~--~I~OOO~-- 5~OO~--0~

Fig. 7. Nyquist plot for Region 11

. 1 Gain between log(xo ) and F. (kmolc / min

4.2 Disturbance Cost.

Fig. 8. Effect of uncertainty in 10g(XR)/F on static DC for LR-B configuration

From Eq. (1) the required control effort to attain perfect disturbance rejection is: ~(s) = - f.-l£afl(s)


4.3 Dynamic Resiliency.


To differentiate between the feasible configurations, Eq. (7) is computed as a function of frequency, as shown in Fig. 9. The analysis indicates that the largest bandwidth for disturbance rejection is possible when configuration L/rB is selected for both regions. Although disturbance rejection in Region I appears to be superior to that in Region 11 (for the L/rB configuration), both provide acceptable bandwidth. As a fmal check, the frequency-dependent Relative Gain Array (RGA, Bristol, 1966) is plotted in Fig. 10. This indicates a low degree of loop interaction suggesting that decentralized diagonal control can be implemented to achieve good tracking performance. Similar performance is predicted for Region I with the same control configuration.

The disturbance cost, DC (Lewin, 1996), is computed by scaling the norm of the required control effort by the maximum possible perturbation in u, lul max and for the worst possible disturbance, Idlmax:

Following scaling, values of the steady-state DC in excess of unity are an indication of unfeasible designs or control configurations. Similarly, violations of DC(ooc) ~ 1 are indicative of bandwidth limitations to perfect disturbance rejection. Decentralized control of the column is desirable due to its robustness and simplicity (Morari and Zafiriou, 1989). Since there are two output variables, and given the number of possible manipulated variables, the total number of candidate 2 x2 control configurations is 32. Direct screening of these candidates is carried by computing the static DC, that is, by computing the worst-case value of Eq. (7) for both manipulated variables for each of the 32 candidates, in the steady state. The conclusions of the static DC screening are that E should be controlled, and the feasible candidate configurations are reduced to only L-B, L-V, L-VI?, L/rB, L/rV, L/rVR for Region I and only L/rB for Region 11. The above analysis is carried out using static gains computed with ±O.05% perturbations in inputs. The static DC is

s. CONCLUSIONS The proposed mode ling approach has been demonstrated in the analysis ofC&R for a relatively complex process. The results obtained for the homogeneous azeotropic column analyzed are comparable to those obtained by rigorous dynamic simulation. The analysis has indicated that uncertainties associated with the linear gains may lead to potential disturbance rejection difficulties, and has highlighted the need to handle this issue. It would appear that operation in the most economic operating region is possible without


Seider W. D., Seader, 1. D. and Lewin. D. R. (1999). Process Design Principles: Synthesis, Analysis and Evaluation, John Wiley and Sons Shinskey, F.G. (1984). "Distillation Control For Productivity and Energy Conservation," McGraw- Hill, 56, 82-89 Weitz, O. and Lewin, D. R. (1996). "Dynamic Controllability and Resiliency Diagnosis using Steady State Process Flowsheet Data," Comp. Chem. Eng., 20(4), 325-335.

having to sacrifice dynamic performance. However, the control of the entrainer flow rate is crucial for successful regulation of the column, due the minimum constraint on E as well as the poor disturbance rejection attainable when E is left uncontrolled.

ACKNOWLEDGEMENT This research was supported by a grant from the Fund for the Promotion of Research at the Technion.

1 . 2r-~---~---~-'-~-~~

-L-B __ L- E ._._ L-V



_ ~-B __ ~-E I---_~


Andersen, H.W., Laroche, L. and Morari M. (1991). "Dynamics of Homogeneous Azeotropic Distillation Columns", 1. & E. C. Res. , 30,1846-1855 Bristol, E. H. (1966). "On a New Measure oflnteractions for Multivariable Process Control," IEEE Trans. Automatic Control, AC -11,133-134 Denn, M. M. and Lavie, R. (1982). "Dynamics of Plants with Recycle", Chem. Eng. Jnl, 24, 55-59 Gilliland, E. R., Gould, L. A. and Boyle, T. 1. (1964). "Dynamic Effect of Material Recycle," JACC, Stanford, 140. Horsley, L. H. (1952). Azeotropic Data, ACS. Jacobsen, E.W., Laroche, L., Morari M., Skogestad S., Andersen, H.W. (1991). "Robust Control of Homogeneous Azeotropic Distillation Columns," AIChEJ., 37, 1810-1823 Kapoor, N., McAvoy, T.J., Marlin T.E. (1986). "Effect of Recycle Structure on Distillation Tower Time Constants," AIChE J., 32, 411-418 Lewin, D. R. (1996). "A Simple Tool for Disturbance Resiliency Diagnosis and Feedforward Control Design," Comp. Chem. Eng. , 20, 13-25 Lewin, D. R., J.-P. Gong, and R. Gani (1996). "Optimal Design and Controllability Assessment for a Plant-wide Benchmark," Proc. of the 13th IFAC World Congress, M, 79-84 Lewin, D. R. (1999). "The Interaction of Design and Control," Proc. of the 7th IEEE Mediterranean Con! on Control and Automation (MED '99), 1384-1391. Marlin, T. E. (1981). "The effect on Control Systems of Viewing the Plant as a System of Integrated Units," Proc. Chem. Process Control JI, Engineering Foundation, Sea Island, GA Moczek, J. S., Otto, R.E., Williams, T. 1. (1965). "Approximation Model for the Dynamic Response of Large Distillation Columns," Chem. Eng. Prog. Symp. Ser., 61,135 Morari, M. and Zafrriou, E. (1989): Robust Process Control, Prentice-Hall, Englewood Cliffs N1.


...... ........ L- V


~- VR





10- 2



(I), (min- ' ]

(a) Region I 1 . 2r-~---~---~..---ro-,

1 .... - - .. - - . - - -... .- ...._ ._;,;~



0.8 ~--------------'

go.6 0.4




10" w, [min- ']

(b) Region 11 Fig. 9. Frequency-dependent Disturbance Cost for Regions I and II











0.6 0.4 02 10"


Fig. 10. Frequency-dependent RGA for LJ