- Email: [email protected]
Optimal setpoint control of complex heat exchanger systems
83
Applications
Optimal Setpoint Control of Complex Heat Exchanger Systems N. Toomarian, Z.J. Palmor and J. Dayan * Faculty of Mechanical Engineering, Teehnion, Israel Institute of Technology, Haifa 32000, Israel This paper presents an approach for calculating the best way of distributing the streams following through a certain class of complex heat exchanger systems in order to achieve maximum heat recovery within the system. A computer code has been developed by which the described method is demonstrated, off-line, for two real cases. This program can be readily integrated into an over-all, on-line computer control system for any complex process consisting of an exchanger system of this class. Using an accurate and detailed heat exchanger model, the exit temperatures of each exchanger are calculated by a simple mathematical procedure based on Gilmour's design method. This procedure has been included in a general model for the complete scheme of the system. The scheme is made up of a series of heat exchanger groups with parallel paths in each group. The optimal distribution of the streams within a group is found by the direct search method of Hooke and Jeeves, modified to include constraints; while the overall optimization of the system is achieved via dynamic programming.
1. Introduction Heat exchanger systems are quite common engineering systems in various process industries. An example of such a system is the crude-oil preheat-
Dr. Zalman J. Palmor has been Senior Lecturer at the Faculty of Mechanical Engineering, Technion, Israel, since 1979. He received the B.Sc. (1966) and M.Sc. (1972) in Mechanical Engineering from the Technion and the Ph.D. (1976) in Chemical Engineering from the City University of New York. In 1977 the City University of New York presented him with the Stanley Katz honor award in Chemical Engineering. His professional career includes project supervision and coordination as a Navy (Israel Defence Force) officer and the development of advanced digital process control algorithms and tuning procedures as Head of the R & D group for advanced control strategies at Taylor Instrument Co. and as Senior Consultant to that company. Dr. Palmor is a member of the Israel Association for Automatic Control (member of IFAC). His current research interests are in the practical instabilities in control systems, the application of control theory to process control and self tuning controllers. He is the author or co-author of numerous reports and published technical papers in these fields.
Keywords: Heat Exchangers, Heat Recovery, Optimization Techniques, Computers, Controls
N. Toomarian is a Ph.D. student in the Nuclear Engineering department of the Technion, Israel. He received the B.Sc. (1978) from the Mechanical Engineering faculty of the Arya-Meher University of Technology in Iran, and the M.Sc. (1982) from the Mechanical Engineering faculty of the Technion, Israel. He has authored three journal papers and 15 reports. His research interests include process control, the thermal-hydraulics of nuclear reactors, and the sensitivity of mathematical models. * [Correspondence should be addressed to Prof. Dayan.] North-Holland Computers in Industry 6 (1985) 83-98
Dr. Joshua Dayan is Associate Professor of Mechanical Engineering at the Technion, Israel Institute of Technology, Haifa, Israel. He received the B.Sc. (1961) and M.Sc. (1964) in Chemical Engineering from the Technion and the Ph.D. (1967) from IIT - Illinois Institute of Technology, Chicago, lllinois. Prior to joining the faculty of Mechanical Engineering at the Technion he worked for the American Oil Co., Whiting, Indiana, and the Institute of Automation, Tel-Aviv, Israel, implementing process computer control systems for various refinery and other industrial processes. Along with his academic career Professor Dayan is an active consultant to many industrial firms in both Israel and the United States in the area of process dynamics and control. Professor Dayan is a member of AIChE and Israeli professional societies. He is the author of 30 published technical papers and numerous reports. His research interests and activities includes process dynamics, computer control systems, mathematical modelling and simulation and the development of various energy related processes.
0166-3615/85/$3.30 © 1985 Elsevier Science Publishers B.V. (North-Holland)
84
Applications
ing system in a crude distillation column, Fig. 1. The system is designed to utilize the sensible heat of the various distilled products leaving the column to heat up the feed before it enters the heating furnace and the distillation column itself. This particular system is used in the present paper to demonstrate the proposed method for calculating the optimal distribution of flow streams among the individual exchangers. Operating the system of heat exchangers under strictly controlled conditions is crucial in saving a substantial amount of energy. Much of the work dealing with the optimization of heat exchanger systems can be divided into two groups: 1. Optimal design of a serial heat exchanger system for a given schematic configuration [1,6,29]. In this work the main goal is to arrive at an optimal area for each individual heat exchanger in the system such that the capital and maintenance cost is minimized. The heat exchanger model used in these works is the simplest one possible, i.e., Q = UA(AT)~nFT. 2. Detailed design of the optimal configuration of the heat exchanger system. Some of the workers in this category [13,16,23] have used the heat exchanger system merely as an example to demonstrate the applicability of their methods for system analysis. Others have dealt specifically with heat exchanger systems [7,8,10,11,14,18,21,22,25]. Minimizing the system costs is still the main goal, and the model used is quite trivial. Unlike the above mentioned two classes of problems, the present work does not deal with the design problem but with the operational problem, i.e. the optimal distribution of the flows among the various heat excangers where the latter are arranged in a fixed and known scheme and their geometry is given. The optimal control is designed to redistribute the flows among the heat exchangers such that the overall heat recovery within the whole system is maximized throughout the operation time. Little attention has been given, so far, to this problem, although in several work [5,12,24] its importance in crude units performance and in other refinery operations has been pointed out. For example, in 1976, Ryskamp et al. [24] claimed that raising the crude preheat temperature by 1 °F for a 100,000 B / D flow rate was worth $10,000 per year in fuel savings (based on 3 5 5 / b b l crude). Dibiano [5] also checked the economic feasibility
(om/mter~ m ],ulu.vlr~
of installing a control system for optimizing exchanger operation and has shown the possibility for much higher savings. The simplicity of the exchanger model used in optimizing the design of the exchangers in the system or in optimizing the configuration of the system enables one to employ methods such as dynamic programming or Pontryagin's maximum principle in a rather straight forward manner. However, to find the optimal distribution of flow in a control application, where the calculations are based on relatively small temperature variations, the needed exchanger model must be more precise than those models used for the design. This latter model should be able to account for the temperature dependency of the heat transfer coefficients, for the thermodynamic properties of the fluids and for the geometrical dimensions of each heat exchanger. A suitable exchanger model as employed in the present work is based on Gilmour's design method as modified by Palmor et al. [20], and further improved by Toomarian [27] (see Appendix A). Along with the use of more complex models for the exchangers the adoption of the required optimization techniques also becomes more involved. Dibiano [5] has pointed out that the heat exchanger trains for the preheat of the crude feed present a complex nonlinear mathematical relationship. Supplying maximum heat to the crude oil while operating within all of the system constraints represents a dynamic control problem which cannot be attempted to any degree of complexity with conventional control. On the other hand, Ryskamp et al. [24] claimed that the optimal distribution always dictates equal temperatures of the parallel flows at the exit of the system. This latter claim essentially contradicts Dibiano's claim, since it is quite simple to drive the flow toward its equalization of the output temperatures. In the present work it is shown that the Ryskamp et al. claim is true only for complete symmetry of the system (geometry, fouling factors, etc.). In all other cases, an accurate search for the optimum is rather complex and Dibiano seems to be correct.
2. Problem Formulation and Method of Solution
In general, a heat exchanger system is composed of several groups of heat exchangers, where
each group includes several parallel lines (trains) of serially connected exchangers. Fig. 1 represents such a system containing K groups. The purpose of building a model for a group of exchangers is to enable the calculation of the overall heat, (Q~)k, exchanged within group k for a given distribution of the streams among the group's lines as well as the calculation of the exit temperature of the cold stream, Tk_l. These heat quantities and temperatures can be computed if the exact group scheme is known and if all the geometrical parameters of the exchangers as well as the operational data (rates of all flows, their inlet temperatures and their physical properties) are given. The overall heat exchanged in group k is the sum of all the heat quantities exchanged through each individual exchanger in the group. (J)4 (1),.~
(Oc,)k = ~_~ E j=
85
N. Toomarian et al. / Optimal Setpoint Control
Computers in Industry
(1)
(O),,j,k,
i=1
where (J)k
= number of parallel lines in group k; k = l , 2. . . . . K (I)j,k = number of exchangers connected in series in line j of group k; j = 1, 2 . . . . ( J ) , Qi,j.k = heat exchanged in heat exchanger i of line j and group k; i = 1, 2 . . . . (I)j,k The calculation of the temperature of the combined cold stream, Tk_ 1, leaving group k is based upon the assumption that the total cold stream
enthalpy, h k_ 1, is the sum of the enthalpies of the individual cold streams leaving each line j in group k.
hk_l = C~[Tk_I]" Tk_I (J)~
=± 2 ((Wc)j.,'Cc[(Too),.j.k] Wcj= 1 •
(2)
For the simple case, where each hot stream passes through a single heat exchanger in any particular group, its inlet temperature is known. Also, the cold inlet flow temperature for an exchanger is known as =
(Tco)i-l.j.,
(3)
except for the first exchanger in each line where (T~i),,j.k = T k.
(4)
For this simple case each of the exchangers is calculated by the Palmor [20] and Toomarian [27] method. The overall group calculations (Fig. 2) are conducted, systematically, starting with the first exchanger of the first line to the next on the same line through all the exchangers in this line. The calculation for the next parallel line is conducted in the same manner and so on to the end. After completing these individual calculations, the overall rate of heat exchanged can be obtained from Eq. 1, and the exit temperature of the cold stream, T k_ 1, is obtainable from Eq. 2. For the more complex case where the hot stream flows through more than one exchanger in the
Hot Streams-Distillates
and pumparounds
[ I
~-----7
F ~-
TI( Cold Crude Feed
i
r
I
I
I
-
Group K
~J
I
L~ Group k
Fig. 1. A Typical Schematic of a Preheating System in a Crude Distillation Process.
Furnace Group I
Distillation Column
86
('omputer~ m
Applications
group (see Case 2) the general scheme of Fig. 2 has to be modified accordingly and tailored to the specific arrangement. Such systems, however, are not very common. Once the procedure for computing the exchanged heat has been established, the optimization method, including the search procedure used to calculate the optimal flow distribution, has to be selected. It is clear that the optimal distribution of flow results in the maximum heat exchanged within the system; it also results in the hottest cold stream leaving the system. Thus, the "return function" which has to be maximized is: max
K
(J)k (/)/.~
E E
E
/,=1 j = l
i=1
max
:
K
w,.k E
(5)
k=l
Routine Entry Obtaining
data/
or e x c h a n g e r s / and streams /
©
¢5 Computi ng
(Q'Tco'Tho)i ,j ,k from Palmor-Toomarian model
1 (Thi)i
Set
,j,k=(Tho)i-l,j,k
> ,),o Yes
lndmtr~
where the flows, W~,~, in all trains and groups are the problem variables. The solution must account for several constraints imposed on the various streams. First, the sum of all parallel flows in each group must be equal to the flow, W,,,,I: ~,/),~
E
=
(6)
i-!
Second, the maximum flow through any exchanger is limited by the maximum allowable pressure drop characteristics for that particular exchanger (for both tube side and shell side).
__<(ae, / ? =< (as>,,>),.,.,)
(7,8)
where the subscripts s and t stand for "shell" and " t u b e " respectively and the additional subscript p stands for "permissible". This is a multistage problem with at least F,['= l I( J)k - l ] variables. The multistage characteristics of the problem together with its large number of variables suggests the use of Discrete Dynamic Programming. For the present problem an exchanger group has been chosen as a stage and the flow rates (W~) in each group as the decision variables. Since the relationships among the rates of heat exchanged in each group, the exit temperature of each group and the input variables to the group, are quite complex functions which cannot be analytically represented or solved without iteration, it is necessary to incorporate a search method for the optimal values along the various computational stages of the dynamic programming. The selected method is a modified version of the direct search method of Hooke and Jeeves [9]. Its details, and the way it is integrated into the dynamic programming along with the problem's constraints, are given in Appendix B.
Yes
3. Case Studies
Calculating (QG)N , TN_1
(Return
)
Fig. 2. Principles of the Computations for a Group of Heat Exchangers.
The optimization method has been applied to the heat exchanger systems of two different refinery crude units. The method has been adapted to each case and the program code provided for calculating the optimal set points for the flows through each line.
N. Toomarian et al. / Optimal Setpoint Control
Computers in Industry
products and intermediate cooling pump arounds are all fixed and known. Table 2 summarizes the various stream data. In all exchangers of this example it has been assumed that the heat exchange coefficient through the fouling layer, HF, is 300 btu/hr-ft2-°F and the permitted pressure drop in all exchangers, both in the tube side and the shell side, has been taken as (AP)s p = ( A P ) t p = 15 psi. The computations of the heat rates exchanged in each of the system's groups have been carried out according to Fig. 2. The program has also been used to study the effect of the initial size of the search step (in M l b / h r ) on the final values of the various streams as well as the effect of the initial starting point (assumed initial values of the parallel flow rates)
3.1. Case 1
-
Consider the crude preheat system described by Fig. 3. The system is comprised of fourteen heat exchangers arranged in three groups whose geometrical data are given in Table 1. The purpose of the optimal control program is to find the best way of distributing the main cold stream (crude feed to an atmospheric distillation tower) among the parallel lines in each group so that the total heat recovery from the system is maximized. The temperature of the crude leaving the system is the highest possible and therefore the extra additional heat needed to bring the crude to the desired feed temperature in the furnace is minimized. For this particular case, the flow rates and temperatures of the hot streams - the distillation column side
Table 1 Geometrical Data for the Heat Exchangers of the System in Case 1 (Fig. 3).
H.E. Numberi, j , k
1.1.1 2.1.1
1.2.1
1.3.1
2.3.1 3.3.1
1.4.1
1.1.2 2.1.2
1.2.2
1.3.2
1.1.3
1.2.3 2.2.3
Inside shell diameter d i (in)
33.0
33.0
39.0
39.0
21.25
33.0
25.0
33.0
34.0
34.0
Length of shell L S (ft.)
15
15
15
15
15
15
15
15
15
15
Tube inside diameter d i (in)
0.56
0.56
0.76
0.76
0.76
0.56
0.56
0.56
0.56
0.56
Tube outside diameter d0(in )
0.75
0.75
1.0
1.0
1.0
0.75
0.75
0.75
0.75
0.75
Num b er of shells NS
2
1
1
2
1
2
1
1
1
2
Num b er of tube passes per shell Nr,r
2
4
8
8
2
4
2
4
4
4
711
178
75
75
80
356
262
178
172
214
Baffle spacing PB (in)
8
8
11
11
11
8
6
8
15
15
Baffle Cut B o
0.09
0.09
0.25
0.25
0.40
0.09
0.02
0.09
0.45
0.25
Place of Cold stream flow 0 --- tubes 1 --- shell
0
0
1
1
1
0
0
0
0
0
Type of flow 0 = parallel 1 --- cross
1
1
1
1
1
1
1
1
1
1
Arrangement of tubes 0 --- square 1 -= triangle
0
0
0
0
1
0
1
0
0
1
Types of tubes 0 --- straight 1 --- u tubes
0
0
0
0
0
0
0
0
0
0
N u m b e r of tubes per pass N
87
88
Applications
('omputerLv ttl
Group 3
Group 2
]tldlL~li~l
Group 1
Tp
L . _
_
__1
Fig. 3. A Schematic Description of the Heat Exchanger System of Case 1.
on the final optimal distribution. The initial starting points for the three groups are given in Fig. 4. Table 3 summarizes the results of this study. These effects can be best recognized by observing the value of the optimal return function of the system, f~[Mbtu/hr] = Q~' + Q~ + Q~. It can be seen that the optimal distribution of the main stream among the parallel paths is independent of the initial point of the direct search, and the size of the initial step. Thus, it supports the assumption that there is a single maximum only. The differences between the sets of optimal distributions calculated from the different initial starting points or by using a different size of initial steps are, at the most, of the order of the final search step (usually less than 10M l b / h r ) . These small differences also support Latour's claim [12] that the obtained maximum of heat duty is quite smooth and flat.
3.2. Case 2
Another crude unit preheating system consists of one group of twelve heat-exchangers (whose geometrical data are given in Table 4) arranged according to the scheme of Fig. 5. It is different from the first case in two major points: 1. The system is a complex single group arrangement of exchangers compared to the simple three group arrangement of the first example. For the optimization of a single group there is no need for the dynamic programming procedure; however, the direct search procedure is still needed. 2. In this system the search for an optimal distribution of the flows applies not only to the main cold stream (the crude feed) but also to the hot HAGO (Heavy Atmospheric Gas Oil) secondary stream.
Table 2 Physical Data for the Streams of Case 1. H.E. Number
HotStreams 1.1.1 2.1.1
Inlet temperature to H.E. T(°F) 600 Flow Rate W (Mlb/hr) Specific Gravity s (60°/60°F) Characterization Factor K Viscosities at 100°F z(cp) Viscosities at 210°F z(cp)
100
0.8611
10.6
ColdStream
1.2.1
1.3.1
2.3.1 3.3.1
1.4.1
1.1.2 2.1.2
1.2.2
1.3.2
1.1.3
1.2.3 2.2.3
1.1.3 1.2.3
600
490
490
660
435
410
400
450
400
90
40
150
150
80
160
72
43
53
132
600
0.8611
10.6
3.72
3.72
1.39
1.39
0.9590
11.8
118 11.8
0.9590
11.8
118 11.8
0.8853
11.8
33 4.4
0.8274
10.8
0.8274
10.8
0.8617
10.6
0.7902
10.8
0.7178
11.2
0.8686
11.6
1.44
1.44
3.72
0.953
0.44
4.7
0.67
0.67
1.39
0.51
0.361
1.63
Optimal return from system f3 MM B t u / h r
108.25 108.3 108.4 108.3
108.2 108.2 108.4 108.2
108.25 108.3 108.4 108.3
108.2 108.2 108.4 108.2
108.25 108.4 108.4 108.4
Starting step size A M lb/hr
150 100 80 50
150 100 80 50
150 100 80 50
150 100 80 50
150 100 80 50
Starting point no.
1
2
3
4
5
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
Optimal heat exchanged (group 3) O~ MM B t u / h r
216 212.5 210 212.5
216 210 210 210
216 212.5 210 212.5
208 217.5 210 217.5
216 212.5 210 212.5
W31
384 387.5 390 387.5
384 390 390 390
384 387.5 390 387.5
392 382.5 390 382.5
384 387.5 390 387.5
W32
Optimal flowrates vector (group 3) W~' M lb/hr
192.8 192.7 192.7 192.7
192.8 192.7 192.7 192.7
192.8 192.7 192.7 192.7
192.7 192.8 192.7 192.8
192.8 192.7 192.7 192.7
Optimal outlet temp. from group 3 7"l* F
28.5 28.5 28.5 28.5
28.5 28.5 28.5 28.5
28.5 28.5 28.5 28.5
28.5 28.45 28.5 28.45
28.5 28.5 28.5 28.5
Optimal heat exchanged (group 2) Q~ MMBtu/hr
Table 3 Optimal Policy for the System of Case 1 (Fig. 3) for Crude Inlet Temperature of 90°F.
W22
187.5 198.5 190 187.5 187.5 187.5 190 187.5 187.5 187.5 190 187.5 188 185 190 185 188 185 190 185
W21
319 325 330 325 321 317.5 330 317.5 319 325 330 325 319 335 330 335 319 325 330 325
93 90 80 90
93 80 80 80
93.5 87.5 80 87.5
91.5 95 80 95
93.5 87.5 80 87.5
W23
Optimal flowrates vector (group 2) W~' MM B t u / h r
284.4 284.3 284.3 284.3
284.4 284.3 284.3 284.3
284.4 284.3 284.3 284.3
284.3 284.2 284.3 284.2
284.4 284.3 284.3 284.3
T~' °F
Optimal outlet temp. from group 2
50.55 50.7 50.7 50.7
50.5 50.5 50.7 50.5
50.55 50.6 50.7 50.6
50.5 50.55 50.5 50.55
50.55 50.5 50.7 50.6
O~' MM B t u / h r
Optimal heat exchanged in group 1
66 62.5 70 62.5 56 62.5 70 62.5 66 62.5 70 62.5 57 60 70 60 66 62.5 70 60
150 162.5 150 162.5 161 155 150 155 150 162.5 150 162.5 160 160 150 160 150 162.5 150 162.5
M lb/hr
w?
300 287.5 290 287.5
301 297.5 290 297.5
300 287.5 290 287.5
300 300 290 300
300 287.5 290 287.5
Optimal flowrates vector (group 1)
84 87.5 90 90
82 82.5 90 82.5
84 87.5 90 87.5
83 82.5 90 82.5
84 87.5 90 87.5
430.7 431.0 431.0 431.0
430.6 430.6 431.0 430.6
430.7 430.8 431.0 430.8
430.5 430.6 431.0 430.6
430.7 430.8 431.0 430.8
oF
To*
Maximum outlet temp. from systems
e~
90
Applications
('omputer~ m lnduMrv
_~_~ .~, ~
(
4
3
*3,, -
2
Group 3
3~
~ IllO
o
16o
300
Group
~
600-
2
Fig. 4. Starting Points for the Calculations of Case 1.
In addition to the above differences, an improvement in the model calculations has been employed in Case 2. The heat exchange coefficients of the exchangers have not been taken as
constants, but were actually calculated from the measured flow rates, and the input and the output temperatures. For this example, the input temperatures of the
H .A.G.O.
fuel oil
J=l
L
F crude
~/_
j=z
+ Fig. 5. A Schematic of the Crude Preheating System of Case 2.
N. Toomarian et a L / Optimal Setpoint Control
Computers in Industry
91
Table 4 Geometrical Data for the Heat Exchangers of Case 2 (Fig. 5)
Heat Exchanger Number, E(i, j)
E(1, j)
E(2, j)
E(3, j)
E(4, j)
Inside shell diameter D i (inch) Length of shell L s (ft) Tube inside diameter, d i (in) Tube outside diameter, d o (in) Number of shells Ns Number of passes for shell Np Number of tubes per pass N Baffle spacing p a(in) Baffle cut Bo Place of cold stream flow 0 -= tubes; 1 --- shell Type of flow 0 --- parallel; 1 -= cross Arrangement of tubes 0 ~- square; 1 ~ triangle Type of tubes 0 -= straight tubes; 1 --- u tubes
33 20 0.56 0.75 1 2 361 7.3 0.2 0 1 0 0
29 20 0.56 0.75 1 2 345 13.2 0.25 0 1 0 0
33 20 0.81 1.00 1 8 54 12.7 0.25 1 1 1 0
33 20 0.81 1.00 1 8 54 12.7 0.25 1 1 1 0
hot streams, H A G O and FO (Fuel Oil), to exchangers E(2, j ) and E(4, j ) are known for all flow rates. However, the input temperatures to exchangers E(1, j ) and E(3, j ) are not known. To calculate these unknown temperatures, as well as the cold stream temperature and the heat exchanged in the system, a special routine [27] whose flowsheet is given in Fig. 6, has been developed. The routine requires a "first guess" for the hot stream inlet temperature to E(1, j). The prediction of this first guess has been based on an observed, almost linear (for small variations of flow rates), dependency of the hot stream temperature entering E(1, j ) on flow rate:
~
Retrivingdata / for exchangersl andstreams/
J! I-
Calculate: E(I,j) andE(2,j) i
CalculateThi(3,j) usingDNh (Eq. 11) Cal culate E(3,j)
+
+ b[ ( Wh)2,j -- ( I~Vh)2.j] .
(9)
.
Thi(1,j):Tha(3,J)
The parameters a and b were found experimentally to be constants whose values for this particular example are: a = 2[°F/M
lb/hr];
I
Ym [CalculateE(4,J~ N~ Thi(3,j)=Tho(4,j) CalculateE(3,j)
b = ~[°F/M lb/hr].
a and b seem to be effected by changes in temperature levels or physical properties of the streams. Thus, these have to be checked and recalculated whenever the above mentioned factors are changed. The temperature of the input hot stream to E(3, j ) (whose exit temperature is the input temperature to E(1, j ) ) can be estimated rapidly by a nondimensional number DN. This number is based on Deliquit's principle [3]: For any heat exchanger with constant flow, this number is a nondimensional factor which stays constant when a change
Calculate Qj =i~iQi
F-N°F
Yes
d
o
QG=~!IQj Fig. 6. Principlesof the Computation for the Complex Group of Case 2.
92
Applications
('omputcr.sm
lndusl O'
Table 5 Flow Data for the Heat Exchanger System of Case 2. Run No.
Hot/ Cold Stream
Specific Gravity gr/cc
Viscosity c.p.
Fuel oil
HAGO
0.9996
0.8880 0.8765
Crude & HAGO Flow Rates (gpm)
Fuel oil
HAG()
1
2
3
4328
8.3 2.13
349 205
740 225
597 154
761
676
759
1820
14.0 4.6
1
Cold Hot
2
Cold Hot
0.9906
0.8908 0.8648
Cold Hot
0.8917 0.8673
1706
13.8 4.64
766
0.9905
Cold Hot
0.8948 0.8667
612
13.8 4.64
540
0.9852
Cold Hot
0.8945 0.8704
663
13.8 4.64
471
0.9875
Cold Hot
0.8935 0.8715
1807
13.8 4.64
514
0.9892
3
4
5
6
Inlet Temp.
in input temperature for one of the streams occurs. For a change in the temperature of the cold stream, this factor is given by:
673
465
464
499
Fuel oil
HAGO
659
309 560
665
292 543
663
295 542
646
289 550
652
290 554
660
288 549
755
468 532 600
600
(10)
the values obtained without the optimal control. Given the flow rates for a particular refinery and the cost of fuel that would have been required to supply this extra AT, the annual saving can be estimated.
and for a change in the hot stream temperature, it is given by:
4. Concluding Remarks and Recommendations for Control Applications
DN~.
T ~ - Too _ Tc'i - Tc'o
T ~ - Th~
DN h -
E!i - Th~
T h ° - Thi-- T h ° - Thi T~i - Thi Tci- T~i
(11)
where (') denotes the temperature after the change. The initial values of the dimensionless factors DN~ and D N h are calculated from the measured temperatures when the computer control program is called in for the first time. Then, the values of the D N ' s are updated within the control program (each time the heat exchanger model is called) and a new set of flow rates is used. The results of the optimization procedure for Case 2 are presented in the following tables: Table 5 is a summary of the data of the various streams as measured at different time intervals and fed to the optimization program. Table 6 summarizes the calculated results of the optimal values for the stream distributions. It can be seen that the search for the optimal distribution of the flows leads to an average increase (improvement) of A ( TG)co = 3.45°F for the cold stream exit temperature over
The problem of optimal operation of heat exchanger systems has been dealt with, employing dynamic programming, a direct search method, and a proper model for the heat exchangers behavior. With the aid of two case studies it has been shown that a proper on-line computer control program can obtain the necessary measured data from the system (at predetermined time intervals) and calculate an optimal distribution of the flows for this particular set of input (measured) data. Of course, an on-line control scheme requires control valves, or set-point stations for local control, manipulated by the computer, in order to transfer the new settings of the flows to the system. However, the program can also be used off-line as an "advisory tool" which then requires the manual intervention of the operator to transfer the results to the system. An economic feasibility study weighing the investment in the control system versus the expected return has to be conducted for each particular case. Some analysis of the expected
Computers in Industry
N. Toomarian et al. / Optimal Setpoint Control
Outlet Temp. of the Flows °F 1
2
3
Final Temp. °F
Heat Duty MM b t u / h r
~o
1
2
3
4
1
2
3
4
l
2
3
4
~o~co
319 399
391 455
418 500
435 592
319 349
362 440
388 434
432 522
220 407
363 451
399 476
449 546
438.69
78.40
334 363
371 481
415 472
456 579
324 339
373 450
423 430
464 553
336 358
369 476
415 459
467 558
462.28
136.01
332 358
369 482
410 466
449 575
327 341
376 455
425 434
465 555
338 361
372 479
415 463
465 559
459.46
131.23
316 328
360 450
393 422
436 552
323 328
375 426
422 423
470 530
340 358
380 456
433 453
489 546
463.95
93.62
321 336
371 458
406 433
454 534
324 336
376 429
423 425
471 535
332 350
368 453
416 440
470 537
465.22
93.57
316 331
367 471
401 434
447 536
332 336
377 444
420 429
468 537
327 346
362 471
407 437
460 536
488.39
99.44
changes in the system (with respect to time) that affect the frequency of the calculations must be made, as the behavior of the changes affects the type of control system needed and its cost. A typical computer run for Case Studies 1 and 2 on an IBM 370/168 computer took about 30 seconds. This affects both the expense of the run and the frequency at which it could be run (since a different computer will perform differently). For the actual systems used in Case Studies 1 and 2, the frequency of changes on both the temperature and the rate of crude entering the systems was quite low. Thus, running the optimization program only once per shift (every 8 hours) was sufficient. It has also been estimated that a supervisory control system using local control with computer manipulated set point stations (a local control computer
performing other control tasks will serve for the necessary data input and output and as a link between the system and the main computer on which the optimization program will run) can pay itself back (through savings) in about one year! It has been pointed out that the program developed for Cases 1 and 2 is capable of finding the optimal distribution regardless of the chosen initial starting points or the size of the initial steps in the direct search portion of the optimization routine. It is clear, though, that a large initial step size results in a large number of calculations while a too small size causes the program to converge very slowly. Thus, in either case the computer time becomes expensive. Palmor has shown [20] that a "good" estimate for the initial step size is about 10% of the total flow. As far as the initial starting
Table 6 Optimal Flow Rates for the System of Case 2 as Calculated by the Computer Program. Run
Trains Flow Rates (gpm)
No.
Crude
1 2 3 4 5 6
HAGO
Fuel Oil
Outlet Temp. °F
Heat Duty MM B t u / h r
1
2
3
1
2
3
1
2
3
(TG)~o
A(TG)co
QG
AQG
383 740 717 421 457 500
622 704 673 521 527 555
681 752 804 531 483 558
204 303 310 177 206 236
256 441 431 190 216 331
151 258 323 171 110 179
71 240 217 134 128 118
197 202 197 144 152 155
280 262 264 211 211 200
441.27 465.42 461.78 469.34 469.03 461.94
2.58 3.14 2.32 5.39 3.71 3.55
80.09 138.81 133.29 96.90 95.84 101.76
1.69 2.80 2.06 3.28 2.28 2.32
where: A(TG)co = (TG)co - (TG)co and AQG = QH - Q~
93
Apphcations
94
p o i n t is concerned, T o o m a r i a n showed [27] that starting with the present actual m e a s u r e d flow rate values is the best e d u c a t e d guess one can make. Nomenclature A a
b C DNE
--
f
-
F~
-
HF
-
h 1 i J
-
j
K k p
-
Q
-
q
-
T
--
U
-
W
-
heat exchanger area [ft 2] constant in Eq. 9 [ ° F / M l b / h r ] constant in Eq. 9 [ ° F / M l b / h r ] specific heat Deliquits ratio of t e m p e r a t u r e differences (Eqs. 10,11) heat exchanger in Case 2 return function correction factor for (AT)~n heat transfer coefficient in the fouling layer [btu/hr-ft2-°F] e n t h a l p y rate [ b t u / l b - h r ] n u m b e r of exchangers in a train exchanger n u m b e r [ 1 . . . I ) n u m b e r of parallel lines line n u m b e r ( 1 . . . J ) n u m b e r of groups in the system group n u m b e r ( 1 . . . K ) pressure heat load or exchanged heat rate [ b t u / h r ] a function representing the d e p e n d e n c y of heat rate on flow rate a n d inlet t e m p e r a ture t e m p e r a t u r e [°F] overall heat transfer coefficient [ b t u / h r ft2-°F] flow rate [ l b / h r ]
Subscripts c cal. ci co G h hi ho i j k
-
o
-
p s
-
t
-
cold calculated cold in cold out ofagroup hot hot in hot out of exchanger i of line j of group k out permissible at shell side at tube side
('omputer~ m Imlustr)
Superscripts *
o p t i m a l value measured value (n) the n'th discrete value ( )' - measured value of t e m p e r a t u r e after a control change.
Appendix Model
A. A S h e l l - a n d - T u b e
Heat
Exchanger
The heat exchanger m o d e l e m p l o y e d in the p r e s e n t work is based on G i l m o u r ' s design method. This m e t h o d c o m b i n e d the classical, empirical, equations for heat transfer coefficients, the heat b a l a n c e a n d the geometrical d i m e n s i o n of the exc h a n g e r into a simple overall equation (see T a b l e I in [19] a n d T a b l e II in [15]). T h e r e are four d i m e n s i o n l e s s factors in G i l m o u r ' s equation which account for the tempera t u r e d r o p across each of the resistance elements in the system. A s s u m i n g that the true t e m p e r a t u r e difference in the exchanger, A Tt, is given by
A T t = ( A T ) , , . FT,
[A-l]
where ( A T ) , . = (Thi - T~°) - ( T h ° - T~) = L M T D . In Thi - To° Z h o - - Tci
The true t e m p e r a t u r e difference can be considered to be the sum of the t e m p e r a t u r e d r o p s on the four resistance elements in the system
aTt = A T s + a T s + A T w + A T F.
[A-2]
T h e four dimensionless factors can then be defined as
T u b e factor:
F T = A T 1 / A Tt ]
Shell factor:
FS = A T s / A ~
W a l l factor: F o u l i n g factor:
~ FW = ATw/AT ti F F = A Tv/za T,
[A-3]
)
These factors depend, as m e n t i o n e d above, on the physical p r o p e r t i e s of the fluid, on the geometrical d i m e n s i o n s of the exchanger, a n d on its performance. Such a p r e s e n t a t i o n enables a clear view of the cause a n d the effect of each factor in the heat exchanger. F o r a p r o p e r design the G i l m o u r m e t h o d requires that the sum of p r o d u c t s (SOP) - the sum
Computers in Industry
N. Toomarian et al. / Optimal Setpoint Control
of dimensionless factors - should be equal to one SOP = F T + FS + F W + F F = 1.
[A-4]
The design parameters are changed until such equality is obtained. In the same way, SOP = 1 has to be satisfied when the exit temperatures are the free parameters. Obtaining SOP = 1 means that the exit temperatures are found. For this purpose the dependency of the SOP on any of the exit temperatures has been studied (for example, the cold stream temperature Too). Large numbers of different types of industrial shell and tubes heat exchangers were checked with hundreds of possible flow rate samples. It was found, in all these examples, that as long as 0.5 < SOP < 1.5 the plot of In(SOP) versus Too lay, very closely, on a straight line [19]. In order to calculate the exact heat transfer coefficient in the fouling layer, for a given working heat exchanger (SOP = 1), one can substitute the measured value of the flow rate along with the inlet and outlet temperatures of both shell and tube sides into Equation (A-4). This c o m p u t a t i o n will be done only once at the beginning of each run whenever it is assumed that this value of heat transfer coefficient remains constant for small changes in the inlet flow rates or in inlet temperatures during the optimization calculations.
A p p e n d i x B. T h e M o d i f i e d H o o k e and J e e v e s Direct Search M e t h o d
As mentioned, it is impossible to present qk and fk in an explicit analytical form and thus, it would be impossible to use classical analytical methods to calculate the m a x i m u m return functions max{ qk(T~ "), Wk) + fk_,[tk(T(k ",, Wk)] }. Instead the direct search method of H o o k e & Jeeves, which does not require an analytical representation of the function and its derivatives, is used. The m e t h o d is based solely on the comparison of the function's values. It is a compromise between methods that are easy to apply but slow to converge and those which are fast to converge but quite difficult to handle. (Sometimes, particularly when there are discontinuities in the return function, it is even advantageous not to use derivatives).
95
Certain improvements and modifications have been added to the original Hooke and Jeeves method. Bell and Pike have introduced [2] the size of the search step as a function of the variable value, A((h,). They have also been able to reduce the n u m b e r of calculations by suggesting the continuation of the local search around each point in the same direction in which improvement was last detected. First, a large step size is employed and t h e n , when no further improvement of the maximized function is detected, one step size is decreased and the search is continued (still in the same direction). De Vogelaere [4] and Tolmin and Smith [26] have suggested other modifications, that apply mostly to some programming problems. These latter are (1), neutralizing the effect of computational errors that can cause, in some cases, a slowing down of the calculations, and (2) limiting the n u m b e r of the c o m p u t e d points to a predetermined maximum. The flowchart of the improved and modified algorithm as merged into the discrete dynamic p r o g r a m m i n g is given in Fig. B-1. The "local search", which is part of the direct search method, Table B-1 Notation for the Flowsheets in Figs. B-1 and B-2. The variables ~k, 0 and q, are points in a J dimensional space; the rest of the variables are unidimensional 0 the previous base point ~p the current base point q, the base point resulting from the current move S~k the functional value at the base point (return function f at ~/,) T~p the output variable value at 4, (exit temperature rk) Sq, the functional value for this move (return function f at (h) ' Tq, the output variable value at ,~ (exti temperature TK) ss the functional value before this move (the largest value so far of the return function attained by the set of exploratory moves) TS the best value so far of the output variable A current step size as a function of the variable value (`p -=-w) 6 "minimum" step size as a function of the variable 4, (,p =- w) p reduction factor for step size, p < 1, q~j one of the coordinate values for q~, j = 1, 2..... J. J number of coordinates for the points (number of parallel paths in the group)
MAXEVA Allowable maximum number of calculations for the return function.
Applications
96
Computct'~ ~Iz /nduYtr~
, ' j , ! j ,, ,MA×EV~ Set: EVA = 1 Calculate: TI,=T k I Obtain:
I
[ S~= Xk(Tk(n),~) .L
©
Local Search
Set:
;j=:j+,j EVA=EVA+I
Calculate:
T,~=T k 1
S~=Xk(TK(n))
4
-
the Neighborhood0f / No
i
ye s
y~e~
~
No
No
Yes
Yes Set:
=~',j = - / j ; 'ij
EVA=EVA+I
= jj+2,,'j
Calculate:
T;
= Tk_ I
S,=Xk(Tk(n),i)
yes
fk(T~])=S,
nl ' ='° I Calculate: T~,=Tk_1 SJ =Xk(T~n)- ,~) Set: SS=S; , TS=T; EVA = EVA+I LOCAL SEARCH Search f o r Maximum
Return Function at the Neiahberhood o f
Fig. B-1. Flowsheet of the Modified Hooke and Jeeves Search Method as Used in the Discrete Dynamic Programming Procedure for the Heat Exchanger Problem.
is given in Fig. B-2. Table B-1 (based on H o o k e and Jeeves' Table III) summarizes the variables used in the program. An example of the progress of the search method is illustrated by Fig. B-3 which represents the way an optimum has been reached for group 2 in Case 1. As with other search methods it is impossible to ensure that the computed maximum is indeed the absolute maximum and not a local one. In most cases, starting the search at remotely located initial points and converging to the same "maximum" is considered as an accepted proof that the maximum is the absolute maximum. Converging to such a maximum from two different initial starting
Fig. B-2. Flowsheet of the "Local Search" Subroutine in the Direct Search Method of Hooke and Jeeves.
points (point numbers 2 and 5) is demonstrated in Fig. B-3. Table 3 shows the optima obtained for all five initial starting points for this group as well as the other optima and initial starting points used in Case 1. The direct search method has been originally developed and used to find the extremum of unconstrained functions. However, the heat exchanger systems contain two kinds of constraints, the constant sum of the flows in parallel paths and the limited allowable pressure drop in each heat exchanger. These constraints have to be integrated into the search procedure. In order to force the procedure to take into account the limit of total flow, a "zero gain" has been assigned to any point representing total flow greater than Wtota~. Thus, the search would "see" such a point as a non improving direction and move away from it. Similarly, the limitation on pressure drops have been integrated into the search procedure. If during the
Computers in Industry
N. Toomarian et al. / Optimal Setpoint Control
97
60C
Stqrtina ooint no. 2 -Direction of computations progres •
Points computed with steps of 50 Optimum with steps of
5
; 40C
50
•
Points computed with steps of 2,5
V
Optimum with steps of
25
Starting point no. 5 - - - D i r e c t i o n of computations progress \
\ \
\
I I
o
Points computed with steps of 50
~)
Optimum with steps of
V
Points computed with steps of 25
V
Optimum computed with steps of 25
50
300
200
IOC
0
I00
200
300
400
500
600
Fig. B-3. "Direct Search" for/'2 with Initial Step Size of 50 for Group 2. Example 1. calculation of a certain heat exchanger train the m a x i m u m pressure d r o p in a n y of the exchangers has been exceeded (too high a flow rate in either the shell or the tubes) the c o n t r i b u t i o n of that p a r t i c u l a r e x c h a n g e r to the overall return function is forced to be zero, thus forcing the p r o g r a m to r e d u c e the flow in o r d e r to gain s o m e t h i n g out of this exchanger. A l t h o u g h not a i m e d for that, such a p r o c e d u r e will also indicate faulty design a r r a n g e m e n t s where in a certain train of exchangers several units are m a d e to transfer high flow rates while others d o not. In other words, it is possible that the o p t i m a l d i s t r i b u t i o n dictates a high flow rate through a certain p a t h in which a small heat exchanger, which c a n n o t take this high flow, is located. T h e p r o g r a m w o u l d then suggest that the " p r o b l e m a t i c " heat exchanger should be bypassed, a n d the o p t i m a l gain is c a l c u l a t e d without it.
References [1] M.A. Ait-Ali and D.J. Wilde, 1980, "Optimal Area Allocation in Multistage Heat-Exchanger System", Transaction of the ASME, Journal of Heat Transfer, Vol. 102, p. 199. [2] M. Bell and M.C. Pike, 1966, "Remark on algorithm 178 (E4) direct search". Comm. ACM, Vol. 9, No. 9, p. 685. [3] C. Deliquit, April 1979, "Simple Method Predicts HeatExchanger Outlet Temperature", Hydr. Proc., p. 213. [4] R. De Vogelaere, 1968, "Remark on algorithm 178 (E4) direct search". Comm. ACM. Vol. 11, No. 7, p. 498. [5] R. Dibiano, Feb. 1981, "Importance versatile control strategy on crude unit", Chem. Eng. Pro., p. 86. [6] L.T. Fan, C.L. Hwang and C.S. Wang, 1965, "Optimization of multistage heat-exchanger system by the discrete maximum principle". Chem. Eng. Prog. Sym. Ser., Vol. 61, No. 59, p. 243. [7] J.R. Flower and B. Linnhoff, 1980, "A thermodynamic combinatorial approach to the design of optimum heat exchanger networks". AIChE Journal, Vol. 26, No. 1, p. 1. [8] C.S. Hwa, 1965. "Mathematical formulation and optimi-
98
[9]
[10]
[11]
[12] [13]
[14]
[15]
[16]
[17]
[18]
Applications
zation of heat exchanger networks using separable programming", AIChE-I. Chem. E. Meeting London, 1965, Sym. Ser., No. 4-Application of math. models in chem. eng. r&d, pp. 101-106. R. Hooke and T.A. Jeeves, 1961, "Direct Search Solution of Numerical and Statistical Problems". J. ACM. Vol. 8, No. 2, p. 212. M.G. Kesler and R.O. Parker, 1969, "Optimal networks of exchanger". Chem. Eng. Prog. Sym. Ser., Vol. 65, No. 92, p. 111. S. Kobayashi, T. Umeda and A. lchikawa, 1971, "Synthesis of optimal heat exchanger systems - an approach by the optimal assignment problem in linear programming". Chem. Eng. Sci., Vol. 26, No. 9, p. 1367. P.R. Latour, June 1979, "On-line computer optimization 1; what it is and where to do it". Hydr. Proc., p. 73. K.F. Lee, D.F. Rudd and A.H. Masso, 1970, "Branch and Bound Synthesis of Integrated Process Designs", Ind. Eng. Chem. Fund., Vol. 9, No. 1, p. 48. B. Linnhoff and J.R. Flower, 1978, "Synthesis of heat exchanger networks", AIChE Journal, Vol. 24, No. 4, p. 633. R.C. Lord, P.E. Minton and R.P. Slusser, Jan. 26, 1970, "Design of heat exchangers", Chem. Eng., Vol. 77, No. 2, p. 96. A.H. Masso and D.F. Rudd, 1969, "The synthesis of system design: Heuristic structuring", AIChE Journal, Vol. 15, No. 1, p. 10. N. Nishida, S. Kobayashi and A. lchikawa, 1971, "Optimal synthesis of heat exchanger system. Necessary conditions for minimum heat transfer area and their application to system synthesis". Chem. Eng. Sci., Vol. 26, p. 1841. N. Nishida, Y.A. Liu and L. Lapidus, 1977, "A Simple and practical approach to the optimal synthesis of heat exchanger networks". AIChE Journal. Vol. 23, No. 1, p. 77.
Computers m lndu.strv
[19] Z. Palmor, J. Dayan and M. Avriel, 1973, ""Evaluating thc performance of shell- and -tube heat exchangers". Israel Journal of Technology. Vol. 11, 4, p. 273. [20] Z. Palmor, 1972, "Determination of optimal set-points for the control of flow distributed through heat-exchangers system". M.Sc. Thesis, Department of Mechanical Engineering, Technion, Haifa~ Israel. [21] T.K. Pho and L. Lapidus, 1973, "'Topics in computer aided design", AIChE Journal, Vol. 19, No. 6, p. 1182. [22] J.W. Ponton and R.A. Donaldson, 1980, "A fast method for the synthesis of optimal heat exchanger networks". AIChE Journal, Vol. 26, No. 1, p. 1. [23] D.F. Rudd and C.C. Watson, 1968, "'Strategy of process engineering", 1st. ed., John Wiley and Sons, p. 212. [24] C.J. Ryskamp, H.L. Wade and R. Britton, May 1976, "Improve crude unit operation". Hydro. Proc., p. 81. [25] J.V. Shah and A.W. Westerberg, 1975, "'Evolutionary Synthesis of Heat Exchanger Networks". Paper No. 60C. AIChE Annual Meeting, Los Angeles, California. [26] F.K. Tolmin and L.B. Smith, 1969, "Direct search. Remark on algorithm 178(E4)." Comm. ACM, Vol. 12, No. 11, p. 637. {27] N. Toomarian, 1981, "'Design of digital computer control system for a crude distillation tower". M.Sc. Thesis, Department of Mechanical Engineering, Technion, Haifa, Israel. [28] N. Toomarian, Z.J. Palmor and J. Dayan, 1982, "Optimal control of parallel streams following through complex heat-exchanger systems". TME-413, Department of Mechanical Engineering, Technion, Haifa, Israel. [29] G.T. Westbrook, 1961, "Use this method to size each stage for best operation". Hydro. Proc. Pet. Ref., Vol. 40, No. 9. p. 201.
Applications
Optimal Setpoint Control of Complex Heat Exchanger Systems N. Toomarian, Z.J. Palmor and J. Dayan * Faculty of Mechanical Engineering, Teehnion, Israel Institute of Technology, Haifa 32000, Israel This paper presents an approach for calculating the best way of distributing the streams following through a certain class of complex heat exchanger systems in order to achieve maximum heat recovery within the system. A computer code has been developed by which the described method is demonstrated, off-line, for two real cases. This program can be readily integrated into an over-all, on-line computer control system for any complex process consisting of an exchanger system of this class. Using an accurate and detailed heat exchanger model, the exit temperatures of each exchanger are calculated by a simple mathematical procedure based on Gilmour's design method. This procedure has been included in a general model for the complete scheme of the system. The scheme is made up of a series of heat exchanger groups with parallel paths in each group. The optimal distribution of the streams within a group is found by the direct search method of Hooke and Jeeves, modified to include constraints; while the overall optimization of the system is achieved via dynamic programming.
1. Introduction Heat exchanger systems are quite common engineering systems in various process industries. An example of such a system is the crude-oil preheat-
Dr. Zalman J. Palmor has been Senior Lecturer at the Faculty of Mechanical Engineering, Technion, Israel, since 1979. He received the B.Sc. (1966) and M.Sc. (1972) in Mechanical Engineering from the Technion and the Ph.D. (1976) in Chemical Engineering from the City University of New York. In 1977 the City University of New York presented him with the Stanley Katz honor award in Chemical Engineering. His professional career includes project supervision and coordination as a Navy (Israel Defence Force) officer and the development of advanced digital process control algorithms and tuning procedures as Head of the R & D group for advanced control strategies at Taylor Instrument Co. and as Senior Consultant to that company. Dr. Palmor is a member of the Israel Association for Automatic Control (member of IFAC). His current research interests are in the practical instabilities in control systems, the application of control theory to process control and self tuning controllers. He is the author or co-author of numerous reports and published technical papers in these fields.
Keywords: Heat Exchangers, Heat Recovery, Optimization Techniques, Computers, Controls
N. Toomarian is a Ph.D. student in the Nuclear Engineering department of the Technion, Israel. He received the B.Sc. (1978) from the Mechanical Engineering faculty of the Arya-Meher University of Technology in Iran, and the M.Sc. (1982) from the Mechanical Engineering faculty of the Technion, Israel. He has authored three journal papers and 15 reports. His research interests include process control, the thermal-hydraulics of nuclear reactors, and the sensitivity of mathematical models. * [Correspondence should be addressed to Prof. Dayan.] North-Holland Computers in Industry 6 (1985) 83-98
Dr. Joshua Dayan is Associate Professor of Mechanical Engineering at the Technion, Israel Institute of Technology, Haifa, Israel. He received the B.Sc. (1961) and M.Sc. (1964) in Chemical Engineering from the Technion and the Ph.D. (1967) from IIT - Illinois Institute of Technology, Chicago, lllinois. Prior to joining the faculty of Mechanical Engineering at the Technion he worked for the American Oil Co., Whiting, Indiana, and the Institute of Automation, Tel-Aviv, Israel, implementing process computer control systems for various refinery and other industrial processes. Along with his academic career Professor Dayan is an active consultant to many industrial firms in both Israel and the United States in the area of process dynamics and control. Professor Dayan is a member of AIChE and Israeli professional societies. He is the author of 30 published technical papers and numerous reports. His research interests and activities includes process dynamics, computer control systems, mathematical modelling and simulation and the development of various energy related processes.
0166-3615/85/$3.30 © 1985 Elsevier Science Publishers B.V. (North-Holland)
84
Applications
ing system in a crude distillation column, Fig. 1. The system is designed to utilize the sensible heat of the various distilled products leaving the column to heat up the feed before it enters the heating furnace and the distillation column itself. This particular system is used in the present paper to demonstrate the proposed method for calculating the optimal distribution of flow streams among the individual exchangers. Operating the system of heat exchangers under strictly controlled conditions is crucial in saving a substantial amount of energy. Much of the work dealing with the optimization of heat exchanger systems can be divided into two groups: 1. Optimal design of a serial heat exchanger system for a given schematic configuration [1,6,29]. In this work the main goal is to arrive at an optimal area for each individual heat exchanger in the system such that the capital and maintenance cost is minimized. The heat exchanger model used in these works is the simplest one possible, i.e., Q = UA(AT)~nFT. 2. Detailed design of the optimal configuration of the heat exchanger system. Some of the workers in this category [13,16,23] have used the heat exchanger system merely as an example to demonstrate the applicability of their methods for system analysis. Others have dealt specifically with heat exchanger systems [7,8,10,11,14,18,21,22,25]. Minimizing the system costs is still the main goal, and the model used is quite trivial. Unlike the above mentioned two classes of problems, the present work does not deal with the design problem but with the operational problem, i.e. the optimal distribution of the flows among the various heat excangers where the latter are arranged in a fixed and known scheme and their geometry is given. The optimal control is designed to redistribute the flows among the heat exchangers such that the overall heat recovery within the whole system is maximized throughout the operation time. Little attention has been given, so far, to this problem, although in several work [5,12,24] its importance in crude units performance and in other refinery operations has been pointed out. For example, in 1976, Ryskamp et al. [24] claimed that raising the crude preheat temperature by 1 °F for a 100,000 B / D flow rate was worth $10,000 per year in fuel savings (based on 3 5 5 / b b l crude). Dibiano [5] also checked the economic feasibility
(om/mter~ m ],ulu.vlr~
of installing a control system for optimizing exchanger operation and has shown the possibility for much higher savings. The simplicity of the exchanger model used in optimizing the design of the exchangers in the system or in optimizing the configuration of the system enables one to employ methods such as dynamic programming or Pontryagin's maximum principle in a rather straight forward manner. However, to find the optimal distribution of flow in a control application, where the calculations are based on relatively small temperature variations, the needed exchanger model must be more precise than those models used for the design. This latter model should be able to account for the temperature dependency of the heat transfer coefficients, for the thermodynamic properties of the fluids and for the geometrical dimensions of each heat exchanger. A suitable exchanger model as employed in the present work is based on Gilmour's design method as modified by Palmor et al. [20], and further improved by Toomarian [27] (see Appendix A). Along with the use of more complex models for the exchangers the adoption of the required optimization techniques also becomes more involved. Dibiano [5] has pointed out that the heat exchanger trains for the preheat of the crude feed present a complex nonlinear mathematical relationship. Supplying maximum heat to the crude oil while operating within all of the system constraints represents a dynamic control problem which cannot be attempted to any degree of complexity with conventional control. On the other hand, Ryskamp et al. [24] claimed that the optimal distribution always dictates equal temperatures of the parallel flows at the exit of the system. This latter claim essentially contradicts Dibiano's claim, since it is quite simple to drive the flow toward its equalization of the output temperatures. In the present work it is shown that the Ryskamp et al. claim is true only for complete symmetry of the system (geometry, fouling factors, etc.). In all other cases, an accurate search for the optimum is rather complex and Dibiano seems to be correct.
2. Problem Formulation and Method of Solution
In general, a heat exchanger system is composed of several groups of heat exchangers, where
each group includes several parallel lines (trains) of serially connected exchangers. Fig. 1 represents such a system containing K groups. The purpose of building a model for a group of exchangers is to enable the calculation of the overall heat, (Q~)k, exchanged within group k for a given distribution of the streams among the group's lines as well as the calculation of the exit temperature of the cold stream, Tk_l. These heat quantities and temperatures can be computed if the exact group scheme is known and if all the geometrical parameters of the exchangers as well as the operational data (rates of all flows, their inlet temperatures and their physical properties) are given. The overall heat exchanged in group k is the sum of all the heat quantities exchanged through each individual exchanger in the group. (J)4 (1),.~
(Oc,)k = ~_~ E j=
85
N. Toomarian et al. / Optimal Setpoint Control
Computers in Industry
(1)
(O),,j,k,
i=1
where (J)k
= number of parallel lines in group k; k = l , 2. . . . . K (I)j,k = number of exchangers connected in series in line j of group k; j = 1, 2 . . . . ( J ) , Qi,j.k = heat exchanged in heat exchanger i of line j and group k; i = 1, 2 . . . . (I)j,k The calculation of the temperature of the combined cold stream, Tk_ 1, leaving group k is based upon the assumption that the total cold stream
enthalpy, h k_ 1, is the sum of the enthalpies of the individual cold streams leaving each line j in group k.
hk_l = C~[Tk_I]" Tk_I (J)~
=± 2 ((Wc)j.,'Cc[(Too),.j.k] Wcj= 1 •
(2)
For the simple case, where each hot stream passes through a single heat exchanger in any particular group, its inlet temperature is known. Also, the cold inlet flow temperature for an exchanger is known as =
(Tco)i-l.j.,
(3)
except for the first exchanger in each line where (T~i),,j.k = T k.
(4)
For this simple case each of the exchangers is calculated by the Palmor [20] and Toomarian [27] method. The overall group calculations (Fig. 2) are conducted, systematically, starting with the first exchanger of the first line to the next on the same line through all the exchangers in this line. The calculation for the next parallel line is conducted in the same manner and so on to the end. After completing these individual calculations, the overall rate of heat exchanged can be obtained from Eq. 1, and the exit temperature of the cold stream, T k_ 1, is obtainable from Eq. 2. For the more complex case where the hot stream flows through more than one exchanger in the
Hot Streams-Distillates
and pumparounds
[ I
~-----7
F ~-
TI( Cold Crude Feed
i
r
I
I
I
-
Group K
~J
I
L~ Group k
Fig. 1. A Typical Schematic of a Preheating System in a Crude Distillation Process.
Furnace Group I
Distillation Column
86
('omputer~ m
Applications
group (see Case 2) the general scheme of Fig. 2 has to be modified accordingly and tailored to the specific arrangement. Such systems, however, are not very common. Once the procedure for computing the exchanged heat has been established, the optimization method, including the search procedure used to calculate the optimal flow distribution, has to be selected. It is clear that the optimal distribution of flow results in the maximum heat exchanged within the system; it also results in the hottest cold stream leaving the system. Thus, the "return function" which has to be maximized is: max
K
(J)k (/)/.~
E E
E
/,=1 j = l
i=1
max
:
K
w,.k E
(5)
k=l
Routine Entry Obtaining
data/
or e x c h a n g e r s / and streams /
©
¢5 Computi ng
(Q'Tco'Tho)i ,j ,k from Palmor-Toomarian model
1 (Thi)i
Set
,j,k=(Tho)i-l,j,k
> ,),o Yes
lndmtr~
where the flows, W~,~, in all trains and groups are the problem variables. The solution must account for several constraints imposed on the various streams. First, the sum of all parallel flows in each group must be equal to the flow, W,,,,I: ~,/),~
E
=
(6)
i-!
Second, the maximum flow through any exchanger is limited by the maximum allowable pressure drop characteristics for that particular exchanger (for both tube side and shell side).
__<(ae, / ? =< (as>,,>),.,.,)
(7,8)
where the subscripts s and t stand for "shell" and " t u b e " respectively and the additional subscript p stands for "permissible". This is a multistage problem with at least F,['= l I( J)k - l ] variables. The multistage characteristics of the problem together with its large number of variables suggests the use of Discrete Dynamic Programming. For the present problem an exchanger group has been chosen as a stage and the flow rates (W~) in each group as the decision variables. Since the relationships among the rates of heat exchanged in each group, the exit temperature of each group and the input variables to the group, are quite complex functions which cannot be analytically represented or solved without iteration, it is necessary to incorporate a search method for the optimal values along the various computational stages of the dynamic programming. The selected method is a modified version of the direct search method of Hooke and Jeeves [9]. Its details, and the way it is integrated into the dynamic programming along with the problem's constraints, are given in Appendix B.
Yes
3. Case Studies
Calculating (QG)N , TN_1
(Return
)
Fig. 2. Principles of the Computations for a Group of Heat Exchangers.
The optimization method has been applied to the heat exchanger systems of two different refinery crude units. The method has been adapted to each case and the program code provided for calculating the optimal set points for the flows through each line.
N. Toomarian et al. / Optimal Setpoint Control
Computers in Industry
products and intermediate cooling pump arounds are all fixed and known. Table 2 summarizes the various stream data. In all exchangers of this example it has been assumed that the heat exchange coefficient through the fouling layer, HF, is 300 btu/hr-ft2-°F and the permitted pressure drop in all exchangers, both in the tube side and the shell side, has been taken as (AP)s p = ( A P ) t p = 15 psi. The computations of the heat rates exchanged in each of the system's groups have been carried out according to Fig. 2. The program has also been used to study the effect of the initial size of the search step (in M l b / h r ) on the final values of the various streams as well as the effect of the initial starting point (assumed initial values of the parallel flow rates)
3.1. Case 1
-
Consider the crude preheat system described by Fig. 3. The system is comprised of fourteen heat exchangers arranged in three groups whose geometrical data are given in Table 1. The purpose of the optimal control program is to find the best way of distributing the main cold stream (crude feed to an atmospheric distillation tower) among the parallel lines in each group so that the total heat recovery from the system is maximized. The temperature of the crude leaving the system is the highest possible and therefore the extra additional heat needed to bring the crude to the desired feed temperature in the furnace is minimized. For this particular case, the flow rates and temperatures of the hot streams - the distillation column side
Table 1 Geometrical Data for the Heat Exchangers of the System in Case 1 (Fig. 3).
H.E. Numberi, j , k
1.1.1 2.1.1
1.2.1
1.3.1
2.3.1 3.3.1
1.4.1
1.1.2 2.1.2
1.2.2
1.3.2
1.1.3
1.2.3 2.2.3
Inside shell diameter d i (in)
33.0
33.0
39.0
39.0
21.25
33.0
25.0
33.0
34.0
34.0
Length of shell L S (ft.)
15
15
15
15
15
15
15
15
15
15
Tube inside diameter d i (in)
0.56
0.56
0.76
0.76
0.76
0.56
0.56
0.56
0.56
0.56
Tube outside diameter d0(in )
0.75
0.75
1.0
1.0
1.0
0.75
0.75
0.75
0.75
0.75
Num b er of shells NS
2
1
1
2
1
2
1
1
1
2
Num b er of tube passes per shell Nr,r
2
4
8
8
2
4
2
4
4
4
711
178
75
75
80
356
262
178
172
214
Baffle spacing PB (in)
8
8
11
11
11
8
6
8
15
15
Baffle Cut B o
0.09
0.09
0.25
0.25
0.40
0.09
0.02
0.09
0.45
0.25
Place of Cold stream flow 0 --- tubes 1 --- shell
0
0
1
1
1
0
0
0
0
0
Type of flow 0 = parallel 1 --- cross
1
1
1
1
1
1
1
1
1
1
Arrangement of tubes 0 --- square 1 -= triangle
0
0
0
0
1
0
1
0
0
1
Types of tubes 0 --- straight 1 --- u tubes
0
0
0
0
0
0
0
0
0
0
N u m b e r of tubes per pass N
87
88
Applications
('omputerLv ttl
Group 3
Group 2
]tldlL~li~l
Group 1
Tp
L . _
_
__1
Fig. 3. A Schematic Description of the Heat Exchanger System of Case 1.
on the final optimal distribution. The initial starting points for the three groups are given in Fig. 4. Table 3 summarizes the results of this study. These effects can be best recognized by observing the value of the optimal return function of the system, f~[Mbtu/hr] = Q~' + Q~ + Q~. It can be seen that the optimal distribution of the main stream among the parallel paths is independent of the initial point of the direct search, and the size of the initial step. Thus, it supports the assumption that there is a single maximum only. The differences between the sets of optimal distributions calculated from the different initial starting points or by using a different size of initial steps are, at the most, of the order of the final search step (usually less than 10M l b / h r ) . These small differences also support Latour's claim [12] that the obtained maximum of heat duty is quite smooth and flat.
3.2. Case 2
Another crude unit preheating system consists of one group of twelve heat-exchangers (whose geometrical data are given in Table 4) arranged according to the scheme of Fig. 5. It is different from the first case in two major points: 1. The system is a complex single group arrangement of exchangers compared to the simple three group arrangement of the first example. For the optimization of a single group there is no need for the dynamic programming procedure; however, the direct search procedure is still needed. 2. In this system the search for an optimal distribution of the flows applies not only to the main cold stream (the crude feed) but also to the hot HAGO (Heavy Atmospheric Gas Oil) secondary stream.
Table 2 Physical Data for the Streams of Case 1. H.E. Number
HotStreams 1.1.1 2.1.1
Inlet temperature to H.E. T(°F) 600 Flow Rate W (Mlb/hr) Specific Gravity s (60°/60°F) Characterization Factor K Viscosities at 100°F z(cp) Viscosities at 210°F z(cp)
100
0.8611
10.6
ColdStream
1.2.1
1.3.1
2.3.1 3.3.1
1.4.1
1.1.2 2.1.2
1.2.2
1.3.2
1.1.3
1.2.3 2.2.3
1.1.3 1.2.3
600
490
490
660
435
410
400
450
400
90
40
150
150
80
160
72
43
53
132
600
0.8611
10.6
3.72
3.72
1.39
1.39
0.9590
11.8
118 11.8
0.9590
11.8
118 11.8
0.8853
11.8
33 4.4
0.8274
10.8
0.8274
10.8
0.8617
10.6
0.7902
10.8
0.7178
11.2
0.8686
11.6
1.44
1.44
3.72
0.953
0.44
4.7
0.67
0.67
1.39
0.51
0.361
1.63
Optimal return from system f3 MM B t u / h r
108.25 108.3 108.4 108.3
108.2 108.2 108.4 108.2
108.25 108.3 108.4 108.3
108.2 108.2 108.4 108.2
108.25 108.4 108.4 108.4
Starting step size A M lb/hr
150 100 80 50
150 100 80 50
150 100 80 50
150 100 80 50
150 100 80 50
Starting point no.
1
2
3
4
5
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
29.2 29.2 29.2 29.2
Optimal heat exchanged (group 3) O~ MM B t u / h r
216 212.5 210 212.5
216 210 210 210
216 212.5 210 212.5
208 217.5 210 217.5
216 212.5 210 212.5
W31
384 387.5 390 387.5
384 390 390 390
384 387.5 390 387.5
392 382.5 390 382.5
384 387.5 390 387.5
W32
Optimal flowrates vector (group 3) W~' M lb/hr
192.8 192.7 192.7 192.7
192.8 192.7 192.7 192.7
192.8 192.7 192.7 192.7
192.7 192.8 192.7 192.8
192.8 192.7 192.7 192.7
Optimal outlet temp. from group 3 7"l* F
28.5 28.5 28.5 28.5
28.5 28.5 28.5 28.5
28.5 28.5 28.5 28.5
28.5 28.45 28.5 28.45
28.5 28.5 28.5 28.5
Optimal heat exchanged (group 2) Q~ MMBtu/hr
Table 3 Optimal Policy for the System of Case 1 (Fig. 3) for Crude Inlet Temperature of 90°F.
W22
187.5 198.5 190 187.5 187.5 187.5 190 187.5 187.5 187.5 190 187.5 188 185 190 185 188 185 190 185
W21
319 325 330 325 321 317.5 330 317.5 319 325 330 325 319 335 330 335 319 325 330 325
93 90 80 90
93 80 80 80
93.5 87.5 80 87.5
91.5 95 80 95
93.5 87.5 80 87.5
W23
Optimal flowrates vector (group 2) W~' MM B t u / h r
284.4 284.3 284.3 284.3
284.4 284.3 284.3 284.3
284.4 284.3 284.3 284.3
284.3 284.2 284.3 284.2
284.4 284.3 284.3 284.3
T~' °F
Optimal outlet temp. from group 2
50.55 50.7 50.7 50.7
50.5 50.5 50.7 50.5
50.55 50.6 50.7 50.6
50.5 50.55 50.5 50.55
50.55 50.5 50.7 50.6
O~' MM B t u / h r
Optimal heat exchanged in group 1
66 62.5 70 62.5 56 62.5 70 62.5 66 62.5 70 62.5 57 60 70 60 66 62.5 70 60
150 162.5 150 162.5 161 155 150 155 150 162.5 150 162.5 160 160 150 160 150 162.5 150 162.5
M lb/hr
w?
300 287.5 290 287.5
301 297.5 290 297.5
300 287.5 290 287.5
300 300 290 300
300 287.5 290 287.5
Optimal flowrates vector (group 1)
84 87.5 90 90
82 82.5 90 82.5
84 87.5 90 87.5
83 82.5 90 82.5
84 87.5 90 87.5
430.7 431.0 431.0 431.0
430.6 430.6 431.0 430.6
430.7 430.8 431.0 430.8
430.5 430.6 431.0 430.6
430.7 430.8 431.0 430.8
oF
To*
Maximum outlet temp. from systems
e~
90
Applications
('omputer~ m lnduMrv
_~_~ .~, ~
(
4
3
*3,, -
2
Group 3
3~
~ IllO
o
16o
300
Group
~
600-
2
Fig. 4. Starting Points for the Calculations of Case 1.
In addition to the above differences, an improvement in the model calculations has been employed in Case 2. The heat exchange coefficients of the exchangers have not been taken as
constants, but were actually calculated from the measured flow rates, and the input and the output temperatures. For this example, the input temperatures of the
H .A.G.O.
fuel oil
J=l
L
F crude
~/_
j=z
+ Fig. 5. A Schematic of the Crude Preheating System of Case 2.
N. Toomarian et a L / Optimal Setpoint Control
Computers in Industry
91
Table 4 Geometrical Data for the Heat Exchangers of Case 2 (Fig. 5)
Heat Exchanger Number, E(i, j)
E(1, j)
E(2, j)
E(3, j)
E(4, j)
Inside shell diameter D i (inch) Length of shell L s (ft) Tube inside diameter, d i (in) Tube outside diameter, d o (in) Number of shells Ns Number of passes for shell Np Number of tubes per pass N Baffle spacing p a(in) Baffle cut Bo Place of cold stream flow 0 -= tubes; 1 --- shell Type of flow 0 --- parallel; 1 -= cross Arrangement of tubes 0 ~- square; 1 ~ triangle Type of tubes 0 -= straight tubes; 1 --- u tubes
33 20 0.56 0.75 1 2 361 7.3 0.2 0 1 0 0
29 20 0.56 0.75 1 2 345 13.2 0.25 0 1 0 0
33 20 0.81 1.00 1 8 54 12.7 0.25 1 1 1 0
33 20 0.81 1.00 1 8 54 12.7 0.25 1 1 1 0
hot streams, H A G O and FO (Fuel Oil), to exchangers E(2, j ) and E(4, j ) are known for all flow rates. However, the input temperatures to exchangers E(1, j ) and E(3, j ) are not known. To calculate these unknown temperatures, as well as the cold stream temperature and the heat exchanged in the system, a special routine [27] whose flowsheet is given in Fig. 6, has been developed. The routine requires a "first guess" for the hot stream inlet temperature to E(1, j). The prediction of this first guess has been based on an observed, almost linear (for small variations of flow rates), dependency of the hot stream temperature entering E(1, j ) on flow rate:
~
Retrivingdata / for exchangersl andstreams/
J! I-
Calculate: E(I,j) andE(2,j) i
CalculateThi(3,j) usingDNh (Eq. 11) Cal culate E(3,j)
+
+ b[ ( Wh)2,j -- ( I~Vh)2.j] .
(9)
.
Thi(1,j):Tha(3,J)
The parameters a and b were found experimentally to be constants whose values for this particular example are: a = 2[°F/M
lb/hr];
I
Ym [CalculateE(4,J~ N~ Thi(3,j)=Tho(4,j) CalculateE(3,j)
b = ~[°F/M lb/hr].
a and b seem to be effected by changes in temperature levels or physical properties of the streams. Thus, these have to be checked and recalculated whenever the above mentioned factors are changed. The temperature of the input hot stream to E(3, j ) (whose exit temperature is the input temperature to E(1, j ) ) can be estimated rapidly by a nondimensional number DN. This number is based on Deliquit's principle [3]: For any heat exchanger with constant flow, this number is a nondimensional factor which stays constant when a change
Calculate Qj =i~iQi
F-N°F
Yes
d
o
QG=~!IQj Fig. 6. Principlesof the Computation for the Complex Group of Case 2.
92
Applications
('omputcr.sm
lndusl O'
Table 5 Flow Data for the Heat Exchanger System of Case 2. Run No.
Hot/ Cold Stream
Specific Gravity gr/cc
Viscosity c.p.
Fuel oil
HAGO
0.9996
0.8880 0.8765
Crude & HAGO Flow Rates (gpm)
Fuel oil
HAG()
1
2
3
4328
8.3 2.13
349 205
740 225
597 154
761
676
759
1820
14.0 4.6
1
Cold Hot
2
Cold Hot
0.9906
0.8908 0.8648
Cold Hot
0.8917 0.8673
1706
13.8 4.64
766
0.9905
Cold Hot
0.8948 0.8667
612
13.8 4.64
540
0.9852
Cold Hot
0.8945 0.8704
663
13.8 4.64
471
0.9875
Cold Hot
0.8935 0.8715
1807
13.8 4.64
514
0.9892
3
4
5
6
Inlet Temp.
in input temperature for one of the streams occurs. For a change in the temperature of the cold stream, this factor is given by:
673
465
464
499
Fuel oil
HAGO
659
309 560
665
292 543
663
295 542
646
289 550
652
290 554
660
288 549
755
468 532 600
600
(10)
the values obtained without the optimal control. Given the flow rates for a particular refinery and the cost of fuel that would have been required to supply this extra AT, the annual saving can be estimated.
and for a change in the hot stream temperature, it is given by:
4. Concluding Remarks and Recommendations for Control Applications
DN~.
T ~ - Too _ Tc'i - Tc'o
T ~ - Th~
DN h -
E!i - Th~
T h ° - Thi-- T h ° - Thi T~i - Thi Tci- T~i
(11)
where (') denotes the temperature after the change. The initial values of the dimensionless factors DN~ and D N h are calculated from the measured temperatures when the computer control program is called in for the first time. Then, the values of the D N ' s are updated within the control program (each time the heat exchanger model is called) and a new set of flow rates is used. The results of the optimization procedure for Case 2 are presented in the following tables: Table 5 is a summary of the data of the various streams as measured at different time intervals and fed to the optimization program. Table 6 summarizes the calculated results of the optimal values for the stream distributions. It can be seen that the search for the optimal distribution of the flows leads to an average increase (improvement) of A ( TG)co = 3.45°F for the cold stream exit temperature over
The problem of optimal operation of heat exchanger systems has been dealt with, employing dynamic programming, a direct search method, and a proper model for the heat exchangers behavior. With the aid of two case studies it has been shown that a proper on-line computer control program can obtain the necessary measured data from the system (at predetermined time intervals) and calculate an optimal distribution of the flows for this particular set of input (measured) data. Of course, an on-line control scheme requires control valves, or set-point stations for local control, manipulated by the computer, in order to transfer the new settings of the flows to the system. However, the program can also be used off-line as an "advisory tool" which then requires the manual intervention of the operator to transfer the results to the system. An economic feasibility study weighing the investment in the control system versus the expected return has to be conducted for each particular case. Some analysis of the expected
Computers in Industry
N. Toomarian et al. / Optimal Setpoint Control
Outlet Temp. of the Flows °F 1
2
3
Final Temp. °F
Heat Duty MM b t u / h r
~o
1
2
3
4
1
2
3
4
l
2
3
4
~o~co
319 399
391 455
418 500
435 592
319 349
362 440
388 434
432 522
220 407
363 451
399 476
449 546
438.69
78.40
334 363
371 481
415 472
456 579
324 339
373 450
423 430
464 553
336 358
369 476
415 459
467 558
462.28
136.01
332 358
369 482
410 466
449 575
327 341
376 455
425 434
465 555
338 361
372 479
415 463
465 559
459.46
131.23
316 328
360 450
393 422
436 552
323 328
375 426
422 423
470 530
340 358
380 456
433 453
489 546
463.95
93.62
321 336
371 458
406 433
454 534
324 336
376 429
423 425
471 535
332 350
368 453
416 440
470 537
465.22
93.57
316 331
367 471
401 434
447 536
332 336
377 444
420 429
468 537
327 346
362 471
407 437
460 536
488.39
99.44
changes in the system (with respect to time) that affect the frequency of the calculations must be made, as the behavior of the changes affects the type of control system needed and its cost. A typical computer run for Case Studies 1 and 2 on an IBM 370/168 computer took about 30 seconds. This affects both the expense of the run and the frequency at which it could be run (since a different computer will perform differently). For the actual systems used in Case Studies 1 and 2, the frequency of changes on both the temperature and the rate of crude entering the systems was quite low. Thus, running the optimization program only once per shift (every 8 hours) was sufficient. It has also been estimated that a supervisory control system using local control with computer manipulated set point stations (a local control computer
performing other control tasks will serve for the necessary data input and output and as a link between the system and the main computer on which the optimization program will run) can pay itself back (through savings) in about one year! It has been pointed out that the program developed for Cases 1 and 2 is capable of finding the optimal distribution regardless of the chosen initial starting points or the size of the initial steps in the direct search portion of the optimization routine. It is clear, though, that a large initial step size results in a large number of calculations while a too small size causes the program to converge very slowly. Thus, in either case the computer time becomes expensive. Palmor has shown [20] that a "good" estimate for the initial step size is about 10% of the total flow. As far as the initial starting
Table 6 Optimal Flow Rates for the System of Case 2 as Calculated by the Computer Program. Run
Trains Flow Rates (gpm)
No.
Crude
1 2 3 4 5 6
HAGO
Fuel Oil
Outlet Temp. °F
Heat Duty MM B t u / h r
1
2
3
1
2
3
1
2
3
(TG)~o
A(TG)co
QG
AQG
383 740 717 421 457 500
622 704 673 521 527 555
681 752 804 531 483 558
204 303 310 177 206 236
256 441 431 190 216 331
151 258 323 171 110 179
71 240 217 134 128 118
197 202 197 144 152 155
280 262 264 211 211 200
441.27 465.42 461.78 469.34 469.03 461.94
2.58 3.14 2.32 5.39 3.71 3.55
80.09 138.81 133.29 96.90 95.84 101.76
1.69 2.80 2.06 3.28 2.28 2.32
where: A(TG)co = (TG)co - (TG)co and AQG = QH - Q~
93
Apphcations
94
p o i n t is concerned, T o o m a r i a n showed [27] that starting with the present actual m e a s u r e d flow rate values is the best e d u c a t e d guess one can make. Nomenclature A a
b C DNE
--
f
-
F~
-
HF
-
h 1 i J
-
j
K k p
-
Q
-
q
-
T
--
U
-
W
-
heat exchanger area [ft 2] constant in Eq. 9 [ ° F / M l b / h r ] constant in Eq. 9 [ ° F / M l b / h r ] specific heat Deliquits ratio of t e m p e r a t u r e differences (Eqs. 10,11) heat exchanger in Case 2 return function correction factor for (AT)~n heat transfer coefficient in the fouling layer [btu/hr-ft2-°F] e n t h a l p y rate [ b t u / l b - h r ] n u m b e r of exchangers in a train exchanger n u m b e r [ 1 . . . I ) n u m b e r of parallel lines line n u m b e r ( 1 . . . J ) n u m b e r of groups in the system group n u m b e r ( 1 . . . K ) pressure heat load or exchanged heat rate [ b t u / h r ] a function representing the d e p e n d e n c y of heat rate on flow rate a n d inlet t e m p e r a ture t e m p e r a t u r e [°F] overall heat transfer coefficient [ b t u / h r ft2-°F] flow rate [ l b / h r ]
Subscripts c cal. ci co G h hi ho i j k
-
o
-
p s
-
t
-
cold calculated cold in cold out ofagroup hot hot in hot out of exchanger i of line j of group k out permissible at shell side at tube side
('omputer~ m Imlustr)
Superscripts *
o p t i m a l value measured value (n) the n'th discrete value ( )' - measured value of t e m p e r a t u r e after a control change.
Appendix Model
A. A S h e l l - a n d - T u b e
Heat
Exchanger
The heat exchanger m o d e l e m p l o y e d in the p r e s e n t work is based on G i l m o u r ' s design method. This m e t h o d c o m b i n e d the classical, empirical, equations for heat transfer coefficients, the heat b a l a n c e a n d the geometrical d i m e n s i o n of the exc h a n g e r into a simple overall equation (see T a b l e I in [19] a n d T a b l e II in [15]). T h e r e are four d i m e n s i o n l e s s factors in G i l m o u r ' s equation which account for the tempera t u r e d r o p across each of the resistance elements in the system. A s s u m i n g that the true t e m p e r a t u r e difference in the exchanger, A Tt, is given by
A T t = ( A T ) , , . FT,
[A-l]
where ( A T ) , . = (Thi - T~°) - ( T h ° - T~) = L M T D . In Thi - To° Z h o - - Tci
The true t e m p e r a t u r e difference can be considered to be the sum of the t e m p e r a t u r e d r o p s on the four resistance elements in the system
aTt = A T s + a T s + A T w + A T F.
[A-2]
T h e four dimensionless factors can then be defined as
T u b e factor:
F T = A T 1 / A Tt ]
Shell factor:
FS = A T s / A ~
W a l l factor: F o u l i n g factor:
~ FW = ATw/AT ti F F = A Tv/za T,
[A-3]
)
These factors depend, as m e n t i o n e d above, on the physical p r o p e r t i e s of the fluid, on the geometrical d i m e n s i o n s of the exchanger, a n d on its performance. Such a p r e s e n t a t i o n enables a clear view of the cause a n d the effect of each factor in the heat exchanger. F o r a p r o p e r design the G i l m o u r m e t h o d requires that the sum of p r o d u c t s (SOP) - the sum
Computers in Industry
N. Toomarian et al. / Optimal Setpoint Control
of dimensionless factors - should be equal to one SOP = F T + FS + F W + F F = 1.
[A-4]
The design parameters are changed until such equality is obtained. In the same way, SOP = 1 has to be satisfied when the exit temperatures are the free parameters. Obtaining SOP = 1 means that the exit temperatures are found. For this purpose the dependency of the SOP on any of the exit temperatures has been studied (for example, the cold stream temperature Too). Large numbers of different types of industrial shell and tubes heat exchangers were checked with hundreds of possible flow rate samples. It was found, in all these examples, that as long as 0.5 < SOP < 1.5 the plot of In(SOP) versus Too lay, very closely, on a straight line [19]. In order to calculate the exact heat transfer coefficient in the fouling layer, for a given working heat exchanger (SOP = 1), one can substitute the measured value of the flow rate along with the inlet and outlet temperatures of both shell and tube sides into Equation (A-4). This c o m p u t a t i o n will be done only once at the beginning of each run whenever it is assumed that this value of heat transfer coefficient remains constant for small changes in the inlet flow rates or in inlet temperatures during the optimization calculations.
A p p e n d i x B. T h e M o d i f i e d H o o k e and J e e v e s Direct Search M e t h o d
As mentioned, it is impossible to present qk and fk in an explicit analytical form and thus, it would be impossible to use classical analytical methods to calculate the m a x i m u m return functions max{ qk(T~ "), Wk) + fk_,[tk(T(k ",, Wk)] }. Instead the direct search method of H o o k e & Jeeves, which does not require an analytical representation of the function and its derivatives, is used. The m e t h o d is based solely on the comparison of the function's values. It is a compromise between methods that are easy to apply but slow to converge and those which are fast to converge but quite difficult to handle. (Sometimes, particularly when there are discontinuities in the return function, it is even advantageous not to use derivatives).
95
Certain improvements and modifications have been added to the original Hooke and Jeeves method. Bell and Pike have introduced [2] the size of the search step as a function of the variable value, A((h,). They have also been able to reduce the n u m b e r of calculations by suggesting the continuation of the local search around each point in the same direction in which improvement was last detected. First, a large step size is employed and t h e n , when no further improvement of the maximized function is detected, one step size is decreased and the search is continued (still in the same direction). De Vogelaere [4] and Tolmin and Smith [26] have suggested other modifications, that apply mostly to some programming problems. These latter are (1), neutralizing the effect of computational errors that can cause, in some cases, a slowing down of the calculations, and (2) limiting the n u m b e r of the c o m p u t e d points to a predetermined maximum. The flowchart of the improved and modified algorithm as merged into the discrete dynamic p r o g r a m m i n g is given in Fig. B-1. The "local search", which is part of the direct search method, Table B-1 Notation for the Flowsheets in Figs. B-1 and B-2. The variables ~k, 0 and q, are points in a J dimensional space; the rest of the variables are unidimensional 0 the previous base point ~p the current base point q, the base point resulting from the current move S~k the functional value at the base point (return function f at ~/,) T~p the output variable value at 4, (exit temperature rk) Sq, the functional value for this move (return function f at (h) ' Tq, the output variable value at ,~ (exti temperature TK) ss the functional value before this move (the largest value so far of the return function attained by the set of exploratory moves) TS the best value so far of the output variable A current step size as a function of the variable value (`p -=-w) 6 "minimum" step size as a function of the variable 4, (,p =- w) p reduction factor for step size, p < 1, q~j one of the coordinate values for q~, j = 1, 2..... J. J number of coordinates for the points (number of parallel paths in the group)
MAXEVA Allowable maximum number of calculations for the return function.
Applications
96
Computct'~ ~Iz /nduYtr~
, ' j , ! j ,, ,MA×EV~ Set: EVA = 1 Calculate: TI,=T k I Obtain:
I
[ S~= Xk(Tk(n),~) .L
©
Local Search
Set:
;j=:j+,j EVA=EVA+I
Calculate:
T,~=T k 1
S~=Xk(TK(n))
4
-
the Neighborhood0f / No
i
ye s
y~e~
~
No
No
Yes
Yes Set:
=~',j = - / j ; 'ij
EVA=EVA+I
= jj+2,,'j
Calculate:
T;
= Tk_ I
S,=Xk(Tk(n),i)
yes
fk(T~])=S,
nl ' ='° I Calculate: T~,=Tk_1 SJ =Xk(T~n)- ,~) Set: SS=S; , TS=T; EVA = EVA+I LOCAL SEARCH Search f o r Maximum
Return Function at the Neiahberhood o f
Fig. B-1. Flowsheet of the Modified Hooke and Jeeves Search Method as Used in the Discrete Dynamic Programming Procedure for the Heat Exchanger Problem.
is given in Fig. B-2. Table B-1 (based on H o o k e and Jeeves' Table III) summarizes the variables used in the program. An example of the progress of the search method is illustrated by Fig. B-3 which represents the way an optimum has been reached for group 2 in Case 1. As with other search methods it is impossible to ensure that the computed maximum is indeed the absolute maximum and not a local one. In most cases, starting the search at remotely located initial points and converging to the same "maximum" is considered as an accepted proof that the maximum is the absolute maximum. Converging to such a maximum from two different initial starting
Fig. B-2. Flowsheet of the "Local Search" Subroutine in the Direct Search Method of Hooke and Jeeves.
points (point numbers 2 and 5) is demonstrated in Fig. B-3. Table 3 shows the optima obtained for all five initial starting points for this group as well as the other optima and initial starting points used in Case 1. The direct search method has been originally developed and used to find the extremum of unconstrained functions. However, the heat exchanger systems contain two kinds of constraints, the constant sum of the flows in parallel paths and the limited allowable pressure drop in each heat exchanger. These constraints have to be integrated into the search procedure. In order to force the procedure to take into account the limit of total flow, a "zero gain" has been assigned to any point representing total flow greater than Wtota~. Thus, the search would "see" such a point as a non improving direction and move away from it. Similarly, the limitation on pressure drops have been integrated into the search procedure. If during the
Computers in Industry
N. Toomarian et al. / Optimal Setpoint Control
97
60C
Stqrtina ooint no. 2 -Direction of computations progres •
Points computed with steps of 50 Optimum with steps of
5
; 40C
50
•
Points computed with steps of 2,5
V
Optimum with steps of
25
Starting point no. 5 - - - D i r e c t i o n of computations progress \
\ \
\
I I
o
Points computed with steps of 50
~)
Optimum with steps of
V
Points computed with steps of 25
V
Optimum computed with steps of 25
50
300
200
IOC
0
I00
200
300
400
500
600
Fig. B-3. "Direct Search" for/'2 with Initial Step Size of 50 for Group 2. Example 1. calculation of a certain heat exchanger train the m a x i m u m pressure d r o p in a n y of the exchangers has been exceeded (too high a flow rate in either the shell or the tubes) the c o n t r i b u t i o n of that p a r t i c u l a r e x c h a n g e r to the overall return function is forced to be zero, thus forcing the p r o g r a m to r e d u c e the flow in o r d e r to gain s o m e t h i n g out of this exchanger. A l t h o u g h not a i m e d for that, such a p r o c e d u r e will also indicate faulty design a r r a n g e m e n t s where in a certain train of exchangers several units are m a d e to transfer high flow rates while others d o not. In other words, it is possible that the o p t i m a l d i s t r i b u t i o n dictates a high flow rate through a certain p a t h in which a small heat exchanger, which c a n n o t take this high flow, is located. T h e p r o g r a m w o u l d then suggest that the " p r o b l e m a t i c " heat exchanger should be bypassed, a n d the o p t i m a l gain is c a l c u l a t e d without it.
References [1] M.A. Ait-Ali and D.J. Wilde, 1980, "Optimal Area Allocation in Multistage Heat-Exchanger System", Transaction of the ASME, Journal of Heat Transfer, Vol. 102, p. 199. [2] M. Bell and M.C. Pike, 1966, "Remark on algorithm 178 (E4) direct search". Comm. ACM, Vol. 9, No. 9, p. 685. [3] C. Deliquit, April 1979, "Simple Method Predicts HeatExchanger Outlet Temperature", Hydr. Proc., p. 213. [4] R. De Vogelaere, 1968, "Remark on algorithm 178 (E4) direct search". Comm. ACM. Vol. 11, No. 7, p. 498. [5] R. Dibiano, Feb. 1981, "Importance versatile control strategy on crude unit", Chem. Eng. Pro., p. 86. [6] L.T. Fan, C.L. Hwang and C.S. Wang, 1965, "Optimization of multistage heat-exchanger system by the discrete maximum principle". Chem. Eng. Prog. Sym. Ser., Vol. 61, No. 59, p. 243. [7] J.R. Flower and B. Linnhoff, 1980, "A thermodynamic combinatorial approach to the design of optimum heat exchanger networks". AIChE Journal, Vol. 26, No. 1, p. 1. [8] C.S. Hwa, 1965. "Mathematical formulation and optimi-
98
[9]
[10]
[11]
[12] [13]
[14]
[15]
[16]
[17]
[18]
Applications
zation of heat exchanger networks using separable programming", AIChE-I. Chem. E. Meeting London, 1965, Sym. Ser., No. 4-Application of math. models in chem. eng. r&d, pp. 101-106. R. Hooke and T.A. Jeeves, 1961, "Direct Search Solution of Numerical and Statistical Problems". J. ACM. Vol. 8, No. 2, p. 212. M.G. Kesler and R.O. Parker, 1969, "Optimal networks of exchanger". Chem. Eng. Prog. Sym. Ser., Vol. 65, No. 92, p. 111. S. Kobayashi, T. Umeda and A. lchikawa, 1971, "Synthesis of optimal heat exchanger systems - an approach by the optimal assignment problem in linear programming". Chem. Eng. Sci., Vol. 26, No. 9, p. 1367. P.R. Latour, June 1979, "On-line computer optimization 1; what it is and where to do it". Hydr. Proc., p. 73. K.F. Lee, D.F. Rudd and A.H. Masso, 1970, "Branch and Bound Synthesis of Integrated Process Designs", Ind. Eng. Chem. Fund., Vol. 9, No. 1, p. 48. B. Linnhoff and J.R. Flower, 1978, "Synthesis of heat exchanger networks", AIChE Journal, Vol. 24, No. 4, p. 633. R.C. Lord, P.E. Minton and R.P. Slusser, Jan. 26, 1970, "Design of heat exchangers", Chem. Eng., Vol. 77, No. 2, p. 96. A.H. Masso and D.F. Rudd, 1969, "The synthesis of system design: Heuristic structuring", AIChE Journal, Vol. 15, No. 1, p. 10. N. Nishida, S. Kobayashi and A. lchikawa, 1971, "Optimal synthesis of heat exchanger system. Necessary conditions for minimum heat transfer area and their application to system synthesis". Chem. Eng. Sci., Vol. 26, p. 1841. N. Nishida, Y.A. Liu and L. Lapidus, 1977, "A Simple and practical approach to the optimal synthesis of heat exchanger networks". AIChE Journal. Vol. 23, No. 1, p. 77.
Computers m lndu.strv
[19] Z. Palmor, J. Dayan and M. Avriel, 1973, ""Evaluating thc performance of shell- and -tube heat exchangers". Israel Journal of Technology. Vol. 11, 4, p. 273. [20] Z. Palmor, 1972, "Determination of optimal set-points for the control of flow distributed through heat-exchangers system". M.Sc. Thesis, Department of Mechanical Engineering, Technion, Haifa~ Israel. [21] T.K. Pho and L. Lapidus, 1973, "'Topics in computer aided design", AIChE Journal, Vol. 19, No. 6, p. 1182. [22] J.W. Ponton and R.A. Donaldson, 1980, "A fast method for the synthesis of optimal heat exchanger networks". AIChE Journal, Vol. 26, No. 1, p. 1. [23] D.F. Rudd and C.C. Watson, 1968, "'Strategy of process engineering", 1st. ed., John Wiley and Sons, p. 212. [24] C.J. Ryskamp, H.L. Wade and R. Britton, May 1976, "Improve crude unit operation". Hydro. Proc., p. 81. [25] J.V. Shah and A.W. Westerberg, 1975, "'Evolutionary Synthesis of Heat Exchanger Networks". Paper No. 60C. AIChE Annual Meeting, Los Angeles, California. [26] F.K. Tolmin and L.B. Smith, 1969, "Direct search. Remark on algorithm 178(E4)." Comm. ACM, Vol. 12, No. 11, p. 637. {27] N. Toomarian, 1981, "'Design of digital computer control system for a crude distillation tower". M.Sc. Thesis, Department of Mechanical Engineering, Technion, Haifa, Israel. [28] N. Toomarian, Z.J. Palmor and J. Dayan, 1982, "Optimal control of parallel streams following through complex heat-exchanger systems". TME-413, Department of Mechanical Engineering, Technion, Haifa, Israel. [29] G.T. Westbrook, 1961, "Use this method to size each stage for best operation". Hydro. Proc. Pet. Ref., Vol. 40, No. 9. p. 201.