A Methodology To Determine Controllability Indices in the Frequency

Jun 18, 2008 - can be transferred into Matlab that is used to calculate the controllability indices ... variables but they can be extended to the freq...
2 downloads 2 Views 4MB Size
Ind. Eng. Chem. Res. 2008, 47, 4807–4816

4807

A Methodology To Determine Controllability Indices in the Frequency Domain Mate Gabor* and Peter Mizsey Department of Chemical and EnVironmental Process Engineering, Budapest UniVersity of Technology and Economics, Budafoki str. 8, H-1521 Budapest, Hungary

The control structure design is an important step of process synthesis. This activity can be supported with controllability indices of the process to be designed. Existing methods to calculate these indices are very complicated and time consuming. In this article, a methodology is presented for the simple determination of controllability indices in the frequency domain. With this methodology, the process of making open loop simulations to approximate transfer function matrices can be avoided. The results obtained can be used to evaluate different control structures. The methodology uses the Control Design Interface module of Aspen Dynamics to describe the system with the linear state space representation. The state space representation can be transferred into Matlab that is used to calculate the controllability indices in the function of frequency from the matrices transferred. The methodology is introduced on two thermally coupled distillation systems (side-stream stripper, SS, and side-stream rectifier, SR). The results show that the methodology is well applicable. The SS has better controllability properties than the SR if the BC separation is easier than AB. In the other cases the SR is better. 1. Introduction Control structure synthesis has become the integral part of process synthesis tasks. This activity should be completed in the early stage of the process design so that a complex evaluation of the design alternatives can be completed. The integrated treatment of equipment and control system design is not new; Ziegler and Nichols1 have already called attention to this problem. They have also introduced the concept of controllability, that is, the ability of a process to achieve and maintain the desired equilibrium value. Since then, several strategies have been developed for control system synthesis, e.g., Luyben, Douglas, Skogestad, and Postlethwaite.2–4 Besides such general strategies for the control system synthesis, there are several studies for complex, integrated investigation of different chemical processes. Mizsey et al.5 have studied energy integrated distillation alternatives; Engelien and Skogestad6 have investigated heat-integrated systems with prefractionator. The different MIMO (multiple input multiple output) control system synthesis strategies are in agreement that the proper control system design can be supported with different analysis tools such as controllability indices. These tools help to get information in the early stage of process design about the controllability behavior of the process, e.g., stability, resiliency, robustness features, and interaction among control loops. Since the process synthesis task is currently supported with professional flowsheeting software tools, it is an obvious intention that such software tools should offer the possibility to investigate the controllability features of the process design alternatives during the synthesis activity. Determining the process model is always a critical part of such an activity. Weitz and Lewin7 presented a procedure for assessing the controllability and resiliency of a process flowsheet using PRO/II for the steady-state calculations. For the dynamic investigation, however, they could not determine the exact dynamic model of the system and they used approximated time constants and delays based on either mixed tank or plug flow assumptions as appropriate. This simplification might be misleading since the * To whom correspondence should be addressed. Email: mate.gabor@ gmail.com.

Figure 1. Connection procedure between Aspen Dynamics and Matlab.

processes are usually more complex than a hydrodynamic approximation, e.g., composition changes. In contrast to Weitz and Lewin,7 in our methodology, we use a more accurate method that involves the linearized dynamic process model determined by professional flowsheet simulator. The Aspen Dynamics flowsheet simulator offers an option, the so-called Control Design Interface (CDI) to get the dynamic model of the process. The CDI determines the differential equation system in the form of state space representation around the operating point of the system to be investigated. Our methodology is based on this option. Moreover, the results obtained with CDI can be directly processed and evaluated with the help of the Matlab software package.

10.1021/ie070952e CCC: $40.75  2008 American Chemical Society Published on Web 06/18/2008

4808 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008 Table 1. Example Script in Aspen Dynamics To Calculate Matrices of the State Space Representation Set Doc ) Simulation set CDI ) Doc.CDI CDI.Reset CDI.AddInputVariable “BLOCKS(““1””).RefluxRatio” CDI.AddInputVariable “BLOCKS(““1””).QRebR” CDI.AddInputVariable “STREAMS(““D2””).FmR” CDI.AddOutputVariable “STREAMS(““D1””).In_F.z(““N-BUT-01””)” CDI.AddOutputVariable “STREAMS(““B1””).In_F.z(““N-PEN-01””)” CDI.AddOutputVariable “STREAMS(““D2””).In_F.z(““2-MET-01””)” CDI.MatricesRequired “A”,“B”,“C”,“D” CDI.Calculate

2. Description of the New Methodology If MIMO systems are to be controlled, it can be completed traditionally in two ways. One possibility is to use controllers having nondiagonal transfer functions. This involves using pre- and postcompensators that counteract the interactions in the plant. As one might expect, this method is very sensitive to modeling errors. In this article, the other method is used and explained in more detail. The so-called decentralized control structure uses diagonal transfer function matrices for the controller. This means that we try to treat an n×n MIMO system as n separate SISO systems. The most important question in this case is how to take the interactions among the control loops into account and then how to pair the controlled and manipulated variables to have good control performance and robustness. To choose the best control structure, the controllability indices can be applied to give information about the behavior of the system. The singular value decomposition provides a useful way of quantifying multivariable directionality. The smallest singular value and the ratio of the largest and the smallest one play an important role when examining a specific control structure.3 This will be explained later in more detail. The relative gain array (RGA)8 gives us information about the interactions among the control loops. These indices were first defined as stationary variables but they can be extended to the frequency domain, for example, as shown by McAvoy et al.9 Calculating frequency-dependent controllability indices is a question of great importance when examining control properties because of their ability to show potential control problems. In this way we can get a deeper insight into the process dynamics. An often-used method can be done in two steps (Hernandez and Jimenez,10 Segovia-Herna´ndez and Herna´ndez11). First, the process identification is carried out by doing open loop simulations and approximating the

Figure 3. Sidestream rectifier (SR). Table 2. Mixtures Investigated

M1 M2 M3

A

B

C

RAC

RBC

RAB

pentane izopentane butane

hexane pentane izopentane

heptane hexane pentane

7.38 3.62 2.95

2.67 2.78 1.3

2.76 1.3 2.26

Table 3. Estimated Total Annual Costs (TAC) of the Design Alternatives TAC (105$/year) product purity (mol%)

SS

SR

4.97 4.55 4.17

4.7 4.16 3.67

12 8.33 7

15.1 10.4 7.06

19.7 11.5 8.83

11.9 9.51 6.71

M1 99 95 90 M2 99 95 90 M3 99 95 90

Table 4. Groups of the Manipulated Variables Examineda SR xA-xB-xC

SS xA-xB-xC

D1-D2-Q1 R1-D2-Q1 D1-R2-Q1 D1-L2-Q1

D1-Q2-Q1 R1-Q2-Q1 L1-Q2-Q1

a The first variable in each group controls the composition of product A, the second controls product B, and the third controls product C. D ) distillate flow rate, L ) reflux flow rate, R ) reflux ratio, Q ) reboiler duty; 1 ) main column, 2 ) side column; xA ) composition of product A, xB ) composition of product B, xC ) composition of product C.

Figure 2. Sidestream stripper (SS).

elements of the transfer function based on the response curves. Second, the indices can be calculated with the use of the transfer function matrix. This method can be time consuming for larger systems and, furthermore, the result of the process identification is usually a first- or second-order unit with dead time element. The methodology introduced here is based on using the state space model of the process that can be calculated with the Control Design Interface module of Aspen Dynamics. With

Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008 4809

Figure 4. Controllability indices for mixture M1 for different groups of manipulated variables (the best group is underlined).

the use of the state space model of any system, more accurate dynamic behavior can be determined than with the use of approximated transfer functions obtained from step responses. The controllability indices can be calculated in the function of frequency in Matlab using the state space model obtained in Aspen Dynamics. The connection procedure between Aspen and Matlab is illustrated in Figure 1. 2.1. SVD and the Controllability Indices Used. Every matrix G can be transformed into its singular value decomposition (SVD) in the following way3 G ) U Σ VH (1) where G is an l×m matrix, Σ is an l×m matrix with k ) min(l,m) nonnegative singular values, σi in its main diagonal;

U is an l×l matrix of output singular vectors, ui; V is an m×m matrix of input singular vectors Vi; and H denotes the complex conjugate transpose of the corresponding matrix. It is important to note that the column vectors of V and similarly the column vectors of U are orthonormal (orthogonal and of unit length). This means that the vectors of V form a basis in the input space and likewise the vectors of U form a basis in the output space. The singular values of the open loop frequency function matrix of a process at a given frequency are the gains of the process at this frequency in the directions of the corresponding input singular vectors. (As the input singular vectors form a basis in the input space the gain can be calculated in every direction.) These gains play an important role when doing

4810 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008

Figure 5. Controllability indices for mixture M2 for different groups of manipulated variables (the best group is underlined).

controllability analysis of a process and for a complex analysis they have to be evaluated at different frequencies. Morari Resiliency Index (MRI). MRI is the smallest singular value of the process open loop frequency function matrix.3 The larger its value the more controllable the process is. If it is zero it means that there is an input direction where the gain is zero and the matrix is not invertible.

Condition Number (CN). CN is the ratio of the largest and smallest singular values of the process open loop frequency function matrix. If it is large then the matrix has strong directionality which means that the gains vary strongly depending on the input directions. Such a matrix is said to be ill-conditioned. A large CN means that the system is sensitive to input and model uncertainty and therefore the process is less controllable.

Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008 4811

Relative Gain Array Number (RGAno). The RGA of a nonsingular square matrix G is a square matrix defined as follows: RGA(G) ) G X (G-1)T (2) where X denotes the element by element multiplication and not the traditional matrix multiplication. T denotes the transpose of the corresponding matrix. The RGA gives information about the interactions in the process and helps us to choose good pairings of controlled and manipulated variables. To measure the interactions with one number the RGAno can be introduced as follows:

RGAno ) |RGA - I|sum

(3)

where | RGA - I |sum is the norm that gives the sum of the absolute values of the matrix and I is the unit matrix. Of course to use the RGAno, first we have to rearrange the inputs and outputs to make the plant diagonally dominant. Pairings with low RGAno are preferred due to weak interactions (RGA is close to unit matrix). 2.2. Determination of the State Space Representation of the Process with Aspen. The linearized state space representation of a process consists of the following system of equations:

Figure 6. Controllability indices for mixture M3 for different groups of manipulated variables (the best group is underlined).

4812 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008

˙ x ) Ax + Bu

(4)

y ) Cx + Du

(5)

where u is the vector of input variables, y is the vector of output variables, x is the vector of state variables, and A, B, C, D are real matrices. Determination of the state space representation is made with the help of the Control Design Interface module of Aspen Dynamics. To do this, a short script is needed to be written in Aspen in which we give the input variables and output variables, and call the appropriate functions to calculate the matrices A, B, C, D of the state space model. In the case of three input and three output variables, a simple script can be used (see Table 1.). It is important that the specification of the input variables always has to be fixed and the output variables have to be free.

In order to use the script, the simulation has to be started and then stopped when steady state is reached. (The system will be linearized around this point. If a system is considered with a set point change in a given range, then this investigation should be repeated at several possible set points. From this investigation, we can conclude how linear the system is and how much the results can be generalized.) When the simulation is stopped, the script can be run. There are five output files of the script. One contains basic information about the results: name, number, and steady state value of the input, output and state variables as well as the time of linearization and the number of nonzero elements of the A, B, C, D matrices. The four other files contain the A, B, C, D matrices in a sparse matrix representation. 2.3. Determination of the Controllability Indices with Matlab. The files containing the A, B, C, D matrices can be loaded and used by Matlab. Then the matrices have to be

Figure 7. Controllability indices for the best manipulated variables for SR and SS (contains one curve from each graph in Figures 3-5, the curve that corresponds to the best group of manipulated variables)

Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008 4813

Figure 8. Closed loop response curves for mixture M1 for the closed loop best manipulated variables (SS: D1-Q2-Q1; SR: R1-D2-Q1).

transformed to full matrices in order to be able to make calculations with them. After that the input and output variables have to be scaled because the MRI and CN depend strongly on the units of measurement. This can be done for example as proposed by Skogestad and Postlethwaite4 by dividing each input and output variable by its desired value (setpoint). This transformation can be described mathematically with diagonal scaling matrices as follows: State space representation after scaling is x˙ ) Ax + BNuus

(6)

ys ) Ny-1Cx + Ny-1DNuus

(7)

where us is the vector of scaled input variables, ys is the vector of scaled output variables, Nu is a diagonal scaling matrix for the input variables containing the steady state values of the input variables in its diagonal, and Ny is a diagonal scaling matrix for the output variables containing the desired values of the output variables in its diagonal.

In Matlab we can give the new B,C,D matrices with the old ones and create the state space model using the ss (state space) function of Matlab: B)B*Nu; C)inv(Ny)*C; D)inv(Ny)*D*Nu; sys)ss(A,B,C,D); The frequency function matrix and the RGAno now can be computed at discrete frequencies for example as shown by Skogestad and Postlethwaite:4 omega)logspace(-5,2,100); for i)1:length(omega) Gf)freqresp(sys,omega(i)); RGAw(:,:,i))Gf.*pinv(Gf).′; RGAno(i))sum(sum(abs(RGAw(:,:,i)-eye(3)))); end Finally the singular values at discrete frequencies can be calculated: H)sigma(sys,omega);

4814 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008

Figure 9. Closed loop response curves for mixture M2 for the closed loop best manipulated variables (SS: D1-Q2-Q1; SR: R1-D2-Q1).

m)max(H); MRI)min(H); for i)1:100 CN(i))m(i)/MRI(i); end 3. Case Studies Our methodology is tested on thermally coupled distillation systems. Such energy-integrated systems have received considerable attention because they often have more favorable energy consumption than the conventional sequences (Annakou and Mizsey,12 Mizsey et al.,5 Hernandez et al.,13,14 and Skogestad15). In this work, two partially thermally coupled distillation systems are investigated, the sidestream stripper (SS, see Figure 2) and the sidestream rectifier (SR, see Figure 3) for the separation of three ternary (ABC) mixtures. SR can be considered as the alternative of the traditional direct sequence and

SS as the alternative of the traditional indirect sequence. With the help of our methodology, the controllability indices are determined to choose the best manipulated variables to control the three product compositions. 3.1. Process Design and Economic Investigation. Since the aim of any process design is the complex investigation of the design alternatives, the optimal design parameters (number of trays, feed tray location, flow rates of the interconnecting streams) are determined for the distillation systems. First, short cut column models are used to give an estimation for the design parameters. Second, these parameters are verified with rigorous simulations and their optimal values are determined where the objective function is the total annual cost (TAC). TAC is calculated for the possible cases using an Excel subroutine called by Aspen Dynamics. For the cost estimation the philosophy of Douglas is followed.2 Three mixtures are examined (Table 2). In terms of the ease of separation, the first mixture is symmetrical; in the second one the BC separation is easier while in the third case the AB

Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008 4815

Figure 10. Closed loop response curves for mixture M3 for the closed loop best manipulated variables (SS: D1-Q2-Q1; SR: D1-D2-Q1).

separation is easier. The feed contains always 33 kmol/h A, 34 kmol/h B, and 33 kmol/h C. The steady-state results show that the SR is economically favorable if the AB separation is easier than BC that is typical for the direct sequence. The SS is better than SR if the BC separation is easier that is typical for the indirect sequence (see Table 3.). This is in accordance with the heuristic rule of leaving the hardest separation at the end. Also the SR is cheaper than the SS if the separation is symmetrical. This is in agreement with the separation heuristics that in such cases the direct sequence is favored to the indirect one. 3.2. Controllability Investigation. Controllability investigation is an integral part of the process design and analysis. The controlled variables are the three product compositions. To accomplish this goal, seven groups of manipulated variables are selected for examination (Table 4). The first manipulated variable of each group controls the composition of product A, the second controls the composition of product B, and the third controls product C. The controllability indices for each mixture for the different manipulated variables can be seen in Figures 4–6. Figure 7 shows the indices for the open loop best

manipulated variables (that are underlined on Figures 4–6) for every mixture. Figures 4–6 show the MRI, CN, and RGAno as the function of frequency. These open loop controllability indices give indication for the proper control structure design. The higher the MRI, the better the manipulated variable set is. The CN shows how the different control loops are conditioned. This value should be close to 1. The RGAno shows the interaction among the control loops and in the optimal case it should be close to zero. These indices should be evaluated in the whole frequency range in order to help the designer to determine the best closed loop control structure. For example, in Figure 4 for the SS, the best control loop pairing is the D1-Q2-Q1 because it has the highest MRI value, its CN is the closest to 1, and its RGAno is also among the best ones. As a consequence, this structure is selected for closed loop investigation. The other figures are evaluated in the same manner. The closed loop simulations are made for three types of disturbances: step change in the feed flow rate (from 100 kmol/h to 101 kmol/h), feed composition (from 0.3300/0.3400/0.3300 to 0.3333/0.3434/0.3233), and set point (from 0.95 to 0.96). PI

4816 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008

controllers are tuned with the Ziegler Nichols method using the automatic tuning tool of Aspen Dynamics in decentralized control structure to maintain the three product compositions at the set point. The quality of the control for different cases is compared with the use of integral absolute error (IAE). Closed loop simulation results are on Figures 8–10. It can be seen in Figures 4–6 that the best group of manipulated variables for the SS is the D1-Q2-Q1 group for every mixture. The best group of manipulated variables for the SR is the D1-L2-Q1 for mixture 1 but the D1-D2-Q1 for mixture 2 and mixture 3. In the case of mixture 1 for the SR the R1-D2-Q1 group has the best RGAno but in the same time it has got the worst MRI and CN values. In every case there is a product stream among the manipulated variables in the best group. In all cases the indices are almost constant until 0.1 rad/s and after they begin to change. It can be seen in Figure 7 that the SR has better controllability properties for the first and third mixtures while the SS is better in the case of the second mixture. This is also in accordance with the heuristic rule of leaving the hardest separation at the end (although this heuristic rule does not take controllability properties into account). In most cases the group of manipulated variables that is chosen using the controllability indices in the frequency domain (Figures 4–6) shows the best controllability properties (IAE values) also in the time domain during closed loop simulations (Figures 8–10). Closed loop results (Figures 8 and 10) show that SR has better controllability properties for mixtures 1 and 3 while SS is better for mixture 2. This is in accordance with the results obtained with our methodology (Figure 7). 4. Conclusions A methodology has been introduced to calculate controllability indices for an Aspen Dynamics simulation in the frequency domain without doing open loop simulations to approximate the transfer function matrix of the process. The methodology works in the ASPEN flowsheeting environment using the statespace representation of the system to be investigated. The MATLAB software supports the evaluation of the controllability indices. The methodology is successfully applied for two partially thermally coupled distillation systems, the sidestream stripper and the sidestream rectifier. The main advantage of this method is that it is fast and simple to apply for different processes. It can be concluded from the results that the controllability indices used are able to show the control properties of the systems for the different manipulated variable sets in most cases. Therefore, the number of the time-consuming closed loop simulations can be reduced and they have to be completed only for the best control structures. Nomenclature SS ) sidestream stripper SR ) sidestream rectifier

MIMO ) multiple input multiple output SISO ) single input single output SVD ) singular value decomposition MRI ) Morari resiliency index CN ) condition number RGA ) relative gain array RGAno ) relative gain array number A, B, C, D ) matrices of the state space representation G ) transfer function matrix of the process IAE ) integral absolute error D1 ) distillate flow rate of the first column D2 ) distillate flow rate of the second column Q1 ) reboiler duty of the first column Q2 ) reboiler duty of the second column R1 ) reflux ratio of the first column

Literature Cited (1) Ziegler, J. G.; Nichols, N. B. Process lags in automatic control circuits. Trans. ASME 1943, 65, 433–444. (2) Douglas, J. Conceptual design of chemical processes; McGraw Hill: New York, 1988. (3) Luyben, L. W. Process modeling, simulation, and control for chemical engineers; McGraw-Hill: Singapore, 1990. (4) Skogestad, S.; Postlethwaite, I. MultiVariable feedback control; Wiley: Chippenham, Wiltshire, UK, 2005. (5) Mizsey, P.; Hau, N. T.; Benko, N.; Kalmar, I.; Fonyo, Z. Process control for energy integrated distillation schemes. Comput. Chem. Eng. 1998, 22, 427–434. (6) Engelien, H. K.; Skogestad, S. Selecting appropriate control variables for a heat-integrated distillation system with prefractionator. Comput. Chem. Eng. 2004, 28, 683–691. (7) Weitz, O.; Lewin, D. R. Dynamic controllability and resiliency diagnosis using steady state process flowsheet data. Comput. Chem. Eng. 1996, 20, 325–335. (8) Bristol, E. On a new mesure of interaction for multivariable process control. IEEE Trans. Autom. Control 1966, AC-11, 133. (9) Mc Avoy, T.; Yaman, A.; Rong, C.; Robinson, D.; Schnelle, P. D. A new approach to defining a dynamic relative gain. Control Eng. Pract. 2003, 11, 907–914. (10) Herna´ndez, S.; Jimenez, A. Controllability Analysis of Thermally Coupled Distillation Systems. Ind. Eng. Chem. Res. 1999, 38, 3957–3963. (11) Herna´ndez, J. G.; Herna´ndez, S. Dynamic Behaviour of Thermally Coupled Distillation Configurations for the Separation of Multicomponent Mixtures. Chem. Biochem. Eng. Q. 2006, 20 (2), 125–133. (12) Annakou, O.; Mizsey, P. Rigorous Comparative Study of EnergyIntegrated Distillation Schemes. Ind. Eng. Chem. Res. 1996, 35, 1877– 1885. (13) Herna´ndez, J. G. S.; Herna´ndez, S.; Rico-Ramı´rez, V.; Jime´nez, A. A comparison of the feedback control behaviour between thermally coupled and conventional distillation schemes. Comput. Chem. Eng. 2004, 28, 811–819. (14) Herna´ndez, J. G. S.; Herna´ndez, S.; Jime´nez, A. Analysis of dynamic properties of alternative sequences to the Petlyuk column. Comput. Chem. Eng. 2005, 29, 1389–1399. (15) Engelien, H. K.; Skogestad, S. Multi-effect distillation applied to an industrial case study. Chem. Eng. Process. 2005, 44, 819–826.

ReceiVed for reView July 12, 2007 ReVised manuscript receiVed April 25, 2008 Accepted May 2, 2008 IE070952E