Dev Chem. Eng. Mineral Process, lO(l/2), p p 19-31,2002 Time Optimal Control of a Binary Distillation Column Rein Luus* Department of Chemical Engineering, University of Toronto, Toronto, Ontario M5S 3E5, Canada lime optimal control problems involving high dimensional systems are generally very dificult to solve. Instead of using suboptimal control based on reduced system, iterative dynamic programming is used directly to establish the time optimal control policy of a binary distillation column described by 11 ordinary differential equations. The resulting optimal control policy that is bang-bang in nature was obtained very rapidly. & minimum time of 2236 seconds obtained here is substantially better than 4120 seconds obtained by a suboptimal control policy. Simulations show that the procedure may be used for on-line optimal control of the distillation column, since the optimal control policy can be established very fast on a personal computer. Introduction Time optimal control problems, where the desired state is to be reached in minimum time, are important in the operation of engineering systems For example, in the startup of a chemical reactor we may wish to reach the desired state of operation as fast as possible, or in the case of some upset in the operation we may wish to recover in minimum time In the operation of a distillation column, some upset in the system may occur, and we may want to recover the desired state of operation as fast as possible. Therefore time optimal control is important in industry and has been studied by numerous researchers. For linear systems one approach that has been used successfully is to convert the continuous system to a system of difference equations by integrating the system over a sampling time, and then to find the minimum number of sampling periods necessary to reach the desired state. Such an approach with linear programming was used for linear systems by Lapidus and Luus [l], Bashein [2], and by Rosen et a1 [3]. For nonlinear systems a similar approach with sequential quadratic programming (SQP) was used by Rosen and Luus [3, 41. Other approaches include the use of dynamic programming and penalty functions [6] and the use of a suboptimal approach [7] For time optimal control of high-dimensional systems the use of model reduction was used by Wong and Luus [8, 91. More recently, Bojkov and Luus [lo, 111 used iterative dynamic programming to solve time optimal control problems. By using a quadratic penalty function with shifting terms, Luus [ 12, 131 showed that the time optimal control problem can be solved successfully for several nonlinear systems that had been difficult to solve by other methods ~~ *Authorfor correspondence (email luus@ecf utoronto cn) 19 The goal here is to test the viability of iterative dynamic programming for the solution of the time optimal control problem involved in transferring the operation of a binary distillation column to its desired state of operation. The distillation column chosen is the one modelled by Davison [14] and used for the study of interaction of control systems by Davison [ 151, and by Davison and Man [161. This model was used for model reduction by Bosley and Lees [ 171, and for time suboptimal control based on reduced models by Wong and Luus [8, 9). The model consists of 8 plates, plus a condenser and a reboiler The control is achieved by changing the heat input to the reboiler and the heat withdrawn from the condenser. There are therefore two control variables and 10 concentration variables, plus the pressure in the column, giving a total of 11 state variables. The detailed development of the equations for the distillation column with pressure variation is given in the appendix of the paper [ 141. Time Optimal Control Problem The model of the binary distillation column with pressure variation, as modelled by Davison [14], is given by the differential equations: 3 = -1.4 x 10-22~+ 4.3 x 10-3z2 dt h 2 dt = 9.5 x 1 0 - 3 ~ 1- 1.38 x 10-222 3.0 x 1 0 - ~ 2~ 3.0~x + 4.6 X 1 0 - 3 ~ 3+ 3.0 x ~ O - ~ U ~ dx3 = 9.5 X 10-3x2 - 1.41 x 10-223 + 6.3 x 1 0 - 3 ~ 4- dt dt &6 dt 5.0 x 1 0 - ~ 2~ 5.0~x + 5.0 x 10-5u2 8.0 x 1 0 - ~ 2-~5.0 ~ x 10-6u1 + 5.0 x ~ o - ~ u ~ - 3.12 x 10-2z5 + 1.5 x 10-2z48.0 x 1 0 - 4 ~ 1 1- 5.0 x ~ O - ~+U5.0~ x ~ O - ~ U = 9.5 x 10-3z4 - 3.52 x 1 0 - 2 ~ 6+ 2.2 x 10-2~3 8.0 x 10-4~11 - 5.0 x 10-%~+ 5.0 x 10-5u2 = 2.02 X 10-225 dx7 - = 2.02 x 10-2x6- 4.22 x 10-2x7+ 2.8 x 1 0 - 2 ~ 3dt 8.0 x 1 0 - ~ 2-~5.0 ~ x 10-%~+ 5.0 x 20 ~ Time Optimal Control of a Binary Distillation Column dZ8 - = 2.02 x 10-2x7 - 4.82 x 1 0 - 2 ~ s 3.7 x 10-2x3 - - = 2.02 x 10-2z8 dt dX9 dt + 6.0 x 10-~211- 2.0 x 10% ~ + 5.0 x 1 0 % ~ - 5.72 x 1 0 - 2 ~ 9+ 4.2 x 3.0 x 1 0 - 4 ~ 1 1 4.0 x + ~ O - ~ 4.0 U ~x 10-2z3 10%~ - - - 2.02 x ~ o - -~4.83 &lo z ~x ~ o - ~ 2.0 x x~ ~ dt kll - dt (9) + 2.0 x 10- %~ (10) + + = 2.55 x 1 0 - 2 ~ 1 2.55 x 10-2~10- 1.85 x 10-2~11 4.6 x 10-4u1 + 4.6 x ~ o - ~ u ~ . (1 1) By using vector-matrix notation, these 11 equations can be written in a more compact form as: dx - = AX + Bu. (12) dt The initial state is: ~ ( 0 ) = [-0.0836 -0.8901 -0.2730 -0.6102 - 0.6135 - 0.3105 -0.9104 0.1298 -0.7639 - 1.0 - 0.2946IT. (13) - The state variables zi, i = 1, ..,10 represent normalized concentration deviations from the desired states of the distillate, the 8 plates, and the bottoms, respectively The state variable 211 is the normalized deviation of the pressure in the column from the pressure at the desired state of operation. Normalization is carried out so that the desired state is the origin The control variables 261 and zb2 are the normalized change of heat input to the reboiler and the normalized heat output from the condenser, respectively. These control variables are bounded by: -1 5 uj 5 1, j = 1,2 (14) and are normalized such that they are zero at the desired state of operation. The time optimal control problem is to find the control policy u(t) in the time interval 0 5 t < t f that transfers the system from the given initial state to the origin in minimum time. The performance index to be minimized is therefore chosen as: I =tf. (15) Numerically, it is required that in a minimum value of the final time t f ,we satisfy the constraint: [Zi(tf)l5 €i, i = 1,2,* . *, 11, (16) 21 where Q is some tolerance, such as the measurement error. For time suboptimal studies for each state variable. of this system, Wong and Luus [8] chose Q to be Direct Approach to Time Optimal Control To solve the optimal control problem, according to the suggestion of Luus [12, 131, we choose the augmented performance index to be minimized as: n J = tf + B C(Zi(tf) - siy, i=1 and update the shifting terms si after every pass of the IDP optimization procedure according to: sf" = si9 - zi(tf)i (18) where the superscript 9 is used to indicate the pass number, and z i ( t f )is the value of the state variable at the end of the pass The product 28si provides sensitivity of the change of the performance index with respect to the equality constraint. Unless sensitivity information is available, the initial choice for si is taken to be zero. The effectiveness of this type of penalty function to handle an equality constraint in the optimal control of a fed-batch reactor was shown by Luus and Hennessy [18]. To measure the closeness of the final states to the origin, we use the norm of the sum of absolute deviations, namely: i= 1 In order to solve the time optimal control problem, we transform the continuous control policy into a piecewise constant control problem by dividing the time interval [0,t f ]into P stages, each of variable length, such that: v(k) = t& - t k - 1 , k = 1,2,...,P. (20) Computationally. we impose the condition: v(k) > 0, k = 1,2, * . * , P , (21) so that we do not deal with negative stage lengths. Let us now introduce a normalized time variable T,by defining d~ in the time interval t k - 1 t < t k through the relationship: < dt = v(k)Pd.r (22) so that: t&- t k - 1 = V ( k ) P ( T k - 71f-1). Therefore: 22 (23) Time Optimal Control of a Binary Distillation Column and in the normalized time, the stages are of equal length and the final time is T = 1,so that iterative dynamic programming (IDP) can be used easily. The only change is that the differential equation now becomes: dx = v ( k ) P ( A x+ Bu) dr in the kth time interval. Algorithm for IDP We consider the stage length at each stage as an additional control variable and augment the control vector to a (3 x 1)vector, but for clarity of presentation we keep it separate We use KDP in a multi-pass fashion, where after a pass the region is restored to a fraction 1 of its value at the beginning of the previous pass. The algorithm can be presented in eight steps as follows 1. Choose the number of time stages P, the number of grid points N, the number of allowable values for control (including stage lengths) R at each grid point, the region contraction factor 7,the region restoration factor 77, initial values for the control and the stage lengths, the initial region sizes, the number of iterations to be used in every pass, and the number of passes. 2. By choosing N values for control and the stage length, evenly distributed inside the specified regions, integrate Eq. (25) from 7 = 0 to T = 1to generate N trajectories. The N values for x at the beginning of each time stage constitute the N grid points at each stage 3 Starting at stage P,corresponding to the normalized time 7 = (P - l)/P,for each grid point generate R sets of values for control and stage length: + D$(P) (26) v ( P ) = v'j(P) + w w j ( P ) (27) u(P) = u'j(P) where D is a 2 x 2 diagonal matrix with different random numbers between -1 and 1 along the diagonal and w is another random number between -1 and 1; u'j(P) is the best value for control, and v*j(P) is the best value for the stage length, both obtained for that particular grid point in the previous iteration; w j is the stage length at iteration j . Integrate Eq (25) from (P - 1 ) / P to T = 1 once with each of the R allowable values for control and stage length to yield x ( t f ) so that the performance index can be evaluated Compare the R values of the augmented performance index and choose the control and stage length which give the minimum value. The corresponding control and stage length are stored for use in step 4. 4. Step back to stage P - 1, corresponding to time r = (P - 2 ) / P . For each grid point generate R allowable sets of control and stage lengths. Integrate Eq. (25) from (P- 2 ) / Pto (P- 1 ) / Ponce with each of the R sets. To continue integration, choose the control and the stage length from step 3 that corresponds to the grid point that is closest to the x at r = (P l ) / P .Now compare the R values of the augmented performance index and store the control policy and the stage length that yield the minimum value. - 23 R Luus 5. Step back to stage P - 2, and continue the procedure in the previous step. Continue until stage 1 corresponding to T = 0 with the given initial state as the grid point is reached. Make the comparison of the R values of the performance index to give the best control and the stage length for this stage We now have the best control policy and the stage length for each stage in the sense of minimizing the augmented performance index from the allowable choices. 6. In preparation for the next iteration, reduce the size of the allowable regions: .i+'(k)=+(k), Wj+'(k) k=1,2,-**,P = p J j ( k ) , k = 1,2,**.,P (28) (29) where 7 is the region reduction factor and j is the iteration index. Use the best control policy and the stage lengths from step 5 as the midpoint for the next iteration. 7. Increment the iteration index j by 1 and go to step 2 to generate another set of grid points. Continue for the specified number of iterations. 8 Increment the pass number index by 1, set the iteration index j to 1 and go to step 2. Continue for the specified number of passes, and examine the results. Computational Results All computations were done in double precision. For integration of the differential equations we used the subroutine DVERK of Hull et a1 1191 with a local error tolerance of to give reliable results. The listing of this subroutine is given in [13]. The problem was solved in FORTRAN on a PentiumIIY600 personal computer. To solve this time optimal control problem by IDP, we chose the region contraction factor used after every iteration 7 = 0.70, region restoration factor 7) = 0.98, a single grid point (N = l), R = 10 randomly chosen points, P = 10 stages and allowed 200 passes, each consisting of 20 iterations. The initial stage length for each stage was taken to be 250 s, and an initial control policy of u1 = 0, u2 = 0 was taken for each stage. The initial region sizes for control and stage length were taken to be 1.0. The first run was performed with 0 = 2 x 10'. After 200 passes requiring 348 s of computation time on a PentiumIIU600we obtained a final time of t f = 2305 with the sum of absolute deviations S = 0.00396. All the state variables were within of the origin. In the resulting control policy, the first 9 stages were at u1 = -1, u2 = 1, for a time of 2061.4 s. The last stage had u1 = 0.17391,~2= -0.20644 and a length of 243 5 s For the next run, the penalty function factor was increased to B = 2 x lo6. At the end of 200 passes, S was decreased to 0.00036, but the final time was increased to t f = 2410.8s. The first 8 stages had the control policy ul = -1,212 = 1. The gth stage of length 241.02 had u1 = 0.33925,u2 = 1, and the last stage of length 246.91 had the policy u1 = 0.02375,~2= -0.04342. The third run with 0 = 2 x lo4 did not yield a final state that was within the required tolerance. The convergence profiles for these three values of the penalty function factor are shown in Figure 1. It is seen that B = 2 x lo5 is a good choice for the penalty function factor. 24 Time Optimal Control of a Binary Distillation Coliunn 1E4 i , 50 PASS 100 150 200 NUMBER Figure 1. Convergence prof2es for different values of the penalqfunctionfactor 8. To refine the results, we now took as the starting condition the first stage of length 2000 s with '111 = -1,912 = 1, and the remaining 9 stages of length 33.3 s. The initial control policy of u = 0 was chosen for these stages, and we used the penalty function factor 8 = 2 x lo5. The run of 200 passes, each consisting of 20 iterations, took a computation time of 135 s, yielding t i = 2236.4 and S = 0.001124 The final state was x ( t / ) = [0.00030 0.00006 0.00008 -0.00009 -0.00021 0.00017 0.00005 -0.00011 0.00003 0.00002 O.OOOOO]T. The resulting optimal control policy is given in Table 1. It can be seen that the time optimal control policy is really a four-stagebang-bang policy, where the first stage is of length 2016.98, the second stage is of length 174.20, the third stage is of length 12 41 and the fourth stage is of length 32.84 The state trajectories corresponding to this control policy are given in Figures 2 4 With the exception of zll that exhibits oscillatory behavior, the states approach the origin in a very smooth manner. We have brought the system closer to the origin than had been required. Therefore upon examination of the output, it was observed that already at the end of the 3gfhpass, all the state variables are within loe3 of the origin. At the end of pass 39 requiring 22.7 s of computation time, the final time t i = 2167.4, S = 0.004485, and the final state is x ( t / ) = [0.00016 -0.00031 -0.00070 -0.00094 -0.00043 0.00048 0.00072 0.00040 0.00020 0.00005 O.OOO1l]T. Here the control policy, as shown in Table 2, is not bang-bang. It is interesting to note that the t / obtained here is almost one-half of the t f = 4120 25 Table 1. lime optimal control policy for the binary distillation column giving t f = 2236.4 s with S = 0 001124. Stage number u1 1 -1.ooooo 2 3 4 5 6 7 8 9 10 -1 .ooooo - 1.00000 1.00000 1.ooooo 1.ooooo 1.00000 1.ooooo 1.ooooo -1 .ooooo Stage length 2000.00 1.00000 10.43 6.55 1.00000 1.00000 34.85 35.59 1.ooooo 38.18 1.00000 37.31 1.00000 28.27 1.00000 12.41 -1.ooooo -1.ooooo 32.84 112 1.ooooo Table 2. Control policyfor the binary distillation column giving tf = 2167.4 s with S = 0.004485. Stage number 1 2 3 4 5 6 7 8 9 10 211 -1.00000 -1.00000 0.98650 1.ooooo 1.00000 1.00000 1.ooooo 1.ooooo -1.00000 -1.ooooo Stage length 2000 00 1.ooooo 23.66 1.00000 28.04 1.ooooo 24 40 1.00000 22.34 1.00000 24.56 1.ooooo 2 1.98 1.ooooo 24.06 -0.24856 27.72 -1.ooooo 31.19 u2 1.00000 obtained for this system by Wong and Luus [9] by the use of time suboptimal control obtained from a reduced model. On-line Time Optimal Control For on-line time optimal control we are faced with the problem of how to control the system while computations are being carried out. Since a good approximation to the optimal control policy can be obtained by using a small number of passes and a small number of random points, we can specify control action for the first time stage and then use IDP to find the best control for the other stages. The first stage is taken as short as possible, but it must still be long enough to enable time optimal calculations to be carried out. For this system we chose a total of 11 stages. We let the first stage be 60 s long and kept u = 0 throughout the run. Thus, no optimization was carried out on this stage. For the other 10 stages we chose each of them initially of length 240 s and the initial control u = 0. Since we wanted the calculations to be executed in less than one minute, we 26 Time Optimal Control of a Binary Distillation Column Figure 2. State trajectoriesfor X I ,22,x3. and x4. -1 2 0 500 lcm 1500 2ooo n m 1 2500 Figure 3. State trajectoriesfor 25,x6, xl,and xg. 27 R Luus Figure 4. State trajectories for xg, 210, and 211. chose R = 5, and allowed only 40 passes with the region restoration factor 7 = 0.95. The other parameters in IDP were the same as before. The computation time for this run was 52 s. After 40 passes, we obtained the control policy in Table 3 with t f = 2273.4, and final state x(t,) = [-0.00009 -0.00030 -0.00071 -0.00086 -0.00025 0.00015 0.00047 0.00051 0.00033 0.00016 0.00008IT. Therefore, we should switch to u1 = -1,112 = 1 after 60 s. We could continue with the control policy in Table 3 and bring the system to the vicinity of the origin (such that all variables are within in magnitude) in a total time of 2273.4 s. However, we can improve the result by carrying out another optimization run where we have for the initial control policy u1 = -1,212 = 1 and a new initial state obtained after the system has been running for 60 s This new initial state is x(0) = [-0.08192 -0.25390 -0.57503 -0.86554 -0.72484 -0.95424 -0.84941 -0.58713 -0.30062 -0.12744 -0.90243jT. We therefore carried out the run with 11 stages, keeping the first stage control policy at 211 = -1,212 = 1, and choosing its initial length to be 1741.66. The remaining 10 stages were chosen of initial length 45.8 and initial control policy u = 0. For this run we chose R = 10, 7) = 0.98, and allowed 200 passes After 200 passes, requiring a computation time of 180 s, we obtained the control policy in Table 4, giving tf = 2182.5 and S = 0.001091. The final state is x ( t f ) = (0.00030 0.00008 0.00010 -0.00008 -0.00024 0.00013 0.00008 -0.00006 0.00000 0.00001 O.OOOO1]T, so all the states are within 0.0003 of the origin. The time required to take the system from the original initial state to the origin is therefore 60 + 2182.5 = 2242.5 s. It is interesting to note that this value of final time is only 28 Time Optimal Control of a Binary Distillation Column 6 s more than the final time obtained by the optimal control policy in Table 1 . This simulation run supports the idea of using optimal control for batch reactors as outlined by Luus and Okongwu [20]. Table 3. Control policy for the binary distillation column where the control policy u = 0 is imposed f o r thejrst stage. Stage number 1 2 3 4 5 6 7 8 9 10 11 u1 0.00000 -1.00000 -1.ooooo -1.ooooo -1 .ooooo -1.ooooo -1.ooooo -1 .ooooo -1.ooooo -0.84869 0.04203 u2 0.00000 1.00000 1.ooooo 1.00000 1.00000 1.ooooo 1.00000 1.00000 1.ooooo 1.ooooo -0.03460 Stage length 60 00 211 02 210 42 216 88 216 20 217.61 222 59 224 01 222 93 232.72 238.02 Table 4. On-line rime optimal control policy for the binary distillation column a f e r the control policy for thejrst stage is obtained in 60 s Stage number 1 2 3 4 5 6 7 8 9 10 11 211 212 -1.ooooo -1.ooooo 1.ooooo 1.ooooo -1.00000 1.00000 1.00000 1.ooooo 1.00000 1.00000 1.00000 1.ooooo -1.00000 -1.00000 -1.00000 -1.ooooo 1 .ooooo 1.ooooo 1.00000 -0.60848 1.00000 -1.00000 Stage length 1742.00 34.34 37.20 40.97 33.99 31.51 36 16 52.64 59 57 68.91 45.21 Conclusions The time optimal control of a binary distillation column was easily established by iterative dynamic programming Since the computation time is small, such optimization can be done on-line and the method should find industrial application. The minimum time to reach the desired state is substantially less than has been reported in the literature. The state trajectories show that the origin is approached very smoothly in spite of the bang-bang control, so there is no concem about the stability of the system once the desired state has been reached. 29 R. Luus Acknowledgement Financial support from the Natural Sciences and Engineering Research Council of Canada is gratefully appreciated. Nomenclature A B D I J L n N P 9 r R si S t t/ % U V 20 xi X State coefficient matrix (11 x 11) Control coefficient matrix (11 x 2) Diagonal matrix of random numbers between -1 and 1 Performance index Augmented performance index Length of a time stage Number of state variables Number of grid points Number of time stages Pass number Region vector over which allowable values for control are chosen Number of randomly chosen values for control Shifting term corresponding to constraint i Sum of absolute values of deviations Time Final time of operation j t h element of control vector u Control vector (2 x 1) Variable stage length Region over which allowable stage lengths are chosen ith state variable State vector (n x 1) Subscripts f Final time i Index j Index Index pertaining to time stage k Superscripts Best value obtained from previous iteration j Iteration step 4 Passnumber Transpose Greek letters Region contraction factor by which the region is reduced after every iteration 7 e Tolerance 71 Region restoration factor 0 Penalty function factor 7 Normalized time variable Random number between -1 and 1 w 30 Time Optimal Control of a Binary Distillation Colwnn References 1 Lapidus, L and LUUS.R 1967 Optimal Control of Engineering Processes Blaisdell, Waltham, M A , 184-189 2 Bashein, G A 1971 A simplex algorithm for on-line computation of time-optimal controls IEEE Trans Autom Control, AC-16,479482 3 Rosen, 0 ,Imanudin, and Luus, R 1987 Sensitivity of the final state specification for time-optimal control problems Int J Control, 45, 1371-1381 4 Rosen, 0 and Luus, R 1989 Sensitivity of optimal control to final state specification by a combined continuation and nonlinear programming approach Chem Eng Sci ,44,2527-2534 5 Rosen, 0 and Luus, R 1991 Evaluation of gradients for piecewise constant optimal control Comput Chem Eng , 15,273-281 6 Edgar, T F and Lapidus, L 1972 The computation of optimal singular control bang-bang control 11: Nonlinear systems AIChE J , 18,780-785 7 Nishida, N ,Liu, Y A ,Lapidus, L ,and Hiratsuka, S 1976 An effective computational algorithm for suboptimal singular ancVor bang-bang control AIChE J ,22,505-523 8 Wong, K T and Luus, R 1983 Time suboptimal feedback control of high order linear systems Int J Control, 37,89-109 9 Wong. K T and LUIS R 1983 Time suboptimal feedback control design through coordinate transformation Int J Control, 37, 723-739 10 Bojkov, B and Luus,R 1994 Time-optimal control by iterative dynamic programming Ind Eng Chem Res ,33, 1486-1492 11 Bojkov, B and Luus,R 1995 Time optimal control of high dimensional systems by iterative dynamic programming Can J Chem Eng ,73,380-390 12 Luus, R 1998 Direct approach to time optimal control by iterative dynamic programming In: Proc fASTED Inrl Conf Intelligent Systems and Control, June 1-4, 1998, Halifax. Nova Scotia, Canada, pp 121-125 13 Luus,R 2000 Iterative Dynamic Programming, Chapman & HalYCRC, London, UK, Chapter 12 14 Davison, E J 1967 Control of a distillation column with pressure variation Trans lnstn Chem Engrs ,45,229-250 15 Davison, E J 1970 The interaction of control systems in a binary distillation column Automatica, 6,447461 16 Davison, E J and Man, F T 1970 Interaction index for multivariable control systems Proc IEE, 117,45942 M J and Lees, F P 1973 The reduction of the order of state variable models using the method of moments Chem Eng Sci .28,2071-2077 17 Bosky, 18 Luus, R and Hennessy, D 1999 Optimization of fed-batch reactors by the Luus-Jaakola optimization procedure Ind Eng Chem Res ,38,1948-1955 19 Hull, T E ,Enright, W D ,and Jackson, K R 1976 User Guide to DVERK - a Subroutine for Solving Nonstiff ODE’S, Report 100. Dept of Comp Sci , Univ of Toronto, Toronto, Canada 20 Luus,R and Okongwu, 0 N 1999 Towards practical optimal conml of batch reactors Chem Eng J .75, 1-9 Received: 4 December 2000; Accepted afier revision: 20 April 2001. 31

1/--страниц