Optimal expansion of a hydrodesulfurization process

Compnfm and Chemic~ Enghwhg, Vol. 4. pp. 101-112 @ Pergam Press Ltd., 1980. Printed in GreatBritain

OPTIMAL EXPANSION OF A HYDRODESULFURIZATION PROCESS D. M. HIMMELBLAU* Department of Chemical Engineering, The University of Texas, Austin, TX 78712,U.S.A.

and T. c.

BICKELt

Sandia Laboratories, Albuquerque, New Mexico,U.S.A. (Received 20 July 1979)

Abstract-The problem solved is: Given deterministic sales forecasts and production requirements over a future time horizon, and a finite set of possible plant expansions to an existing process, when and how should the chemical plant be expanded in order to satisfy the production requirements and maximize the present worth of net profit? A branch and bound procedure combined with nonlinear programming was developed, and applied to a hydrodesulfurization plant to illustrate the technique, its advantages, and its limitations. Scope-Expansion of a chemical plant is one aspect of the more general problem of determining the optimal configuration (flowsheet) for an entire chemical plant. Expansion problems differ from the more general problem in that a portion (perhaps most) of the configuration of equipment is already in existence and thus is fixed with respect to arrangement. A major limitation of previous work on the chemical plant expansion problem is the assumption of independence between possible expansion projects. This assumption leads to the conclusion that the addition of an expansion to the chemical process can increase production, but the assumption holds only if all possible expansions are complete within themselves and can operate independently of all other equipment in the plant. For individual pieces of equipment in a process, this assumption does not generally hold true. As a result, a new procedure has been developed to determine the optimal expansion schedule of an existing chemical plant for a prespecified set of n possible expansions. The objective was to maximize the present worth of net profit, subject to the constraint that production must satisfy a prespecified demand which is increasing over the time horizon. Four significant questions answered were: (1) What kind of equipment should be installed in the.plant? (2) When should the expansion equipment be installed? (3) What would be the optimal operating conditions of the chemical process during each interval over the time horizon?; and (4) When should existing equipment that is obsolete be phased out of production? The expansion time horizon was discretized into p expansion intervals. For each pseudo-steady state interval, linear and nonlinear equality and inequality constraints were developed to model the example HDS process. To avoid a complete enumeration of all of the (p t 1)” possible discrete solutions to the problem, the resulting nonlinear integer programming problem was solved using a branch and bound algorithm that included the generalized reduced gradient constrained nonlinear optimization algorithm, with numerical evaluation of derivatives. The branch and bound algorithm successively partitioned the discrete solution space into smaller sets. An upper bound was calculated for the capacity expansidn objective function for each set. A model of a fuel oil hydrodesulfurization (HDS) process was developed to test and evaluate the chemical process plant expansion algorithm. Although the modelling is somewhat simplified, it served to illustrate the major features that have to be accomodated in any plant model, and had the advantage that the process data were available in the open literature. For the three possible expansions and four expansion intervals (over a 10yr time period) considered, the problem comprised 57 variables, 16 linear equality constraints, 21 nonlinear equality constraints, 6 linear inequality constraints, and 114 bounds on the variables. Conclusions and SigrWanee-Because of the high cost of capital investment, coupled with increased feedstock costs and operation expenses, capacity expansion decisions are subject to considerable risk. Any quantitative ways of reducing the risk such as described here are useful tools for an engineer. tFormerIy, Assistant Instructor, Departmentof Chemical Engineering, The University of Texas at Austin, Austin, TX 78712,U.S.A. *Correspondence should be addressedto this author. 101

102

D. M. HIMMELBLAUand T. C.

BICKEL

Figure 2 shows the example HDS plant and Fig. 3 shows the branch and bound enumeration tree generated during the solution of the example expansion problem. Only 47 problems (vertices) of the total possible 125, in the complete enumeration tree had to be solved, i.e. an upper bound to the objective function was calculated with the nonlinear optimization code. The search terminated at 8 vertices because the relaxed problem at the vertex was found to be infeasible and at 18 vertices because of the upper bounding criterion developed in this work. Of the three possible expansions, only two were utilized. While the algorithm required considerable computation time, it was concluded that it was an effective tool when planning for the future. As the size and complexity of the chemical plant increases, the calculation time of the nonlinear optimization problem will increase at a greater than linear rate, and therefore the total calculation time of the capacity expansion problem will increase correspondingly. This increased cost of calculation time, and the stochastic nature of future demand, must be weighed against the potential costs of poor decisions in capacity expansion implementation.

1. INTRODUCTION

Decision making for chemical plant expansion requires that an engineer maintain an extensive and complex set of skills. The high cost of capital investment, coupled with increased feedstock costs and operating expenses, make capacity expansion decisions subject to considerable risk. Chemical plant expansion problems at a single location form a subset of the more general problem in synthesis, that of determining the optimal configuration for an entire chemical plant. Expansion problems differ from this more general problem in that a portion (perhaps most) of the configuration of equipment is already in existence, and, thus, is essentially fixed with respect to arrangement. Only a fraction of the plant is modified or augmented. Given a deterministic sales forecast or production requirements for the future over a fixed time horizon of, say 5-10 yr, typical questions to be answered are: (1) What kind of equipment should be installed? (2) When should the expansion equipment be installed? (3) What are the optimal operating conditions for the new and old units; and (4) When should units that are obsolete be phased out of production? Such decision problems are among the most difficult to treat because of their combinatorial properties. Plant expansions represent discrete entities, and as a consequence, the capacity expansion problem can be viewed as a discrete decision problem in which a decision is required for each possible expansion project during each capacity expansion time interval. A decision is made in each future time interval as to whether or not capacity expansion project t should be built during the time interval /c, and the decision itself is a yes-no (O-l) variable. It is well known that the total number of solutions for consideration in the capacity expansion problem will be (p + l)“, where p is the number of capacity expansion intervals, and n is the number of possible capacity expansion projects. We describe a method that employs a branch and bound technique to determine the optimal solution of the overall discrete capacity expansion problem with nonlinear programming used to solve each subproblem of the overall problem. The method is illustrated by application to a hydrodesulfurization (HDS) plant. The important difference between the chemical plant capacity expansion problem and the typical discrete problem is that in the former each subproblem in the decision tree is a nonlinear optimization problem, whereas in the latter, a linear subproblem is normally encountered. In developing the algorithm, the objectives were (1) to develop a

quantitative approach that could augment “rules of thumb” in decision making, (2) to make the algorithm general so as to accommodate most steady state plant models, (3) to make it easy to implement and modify, and (4) to keep the computational cost as low as possible. Figure 1 illustrates the main components of the overall algorithm. We first summarize the HDS process model, and then the features of the branch and bound algorithm used to solve the overall capacity expansion problem. Problems in meshing the chemical plant model with the optimization code are described together with some typical computational results. 2. PRRVIOUSWORK

Previous investigators have treated the expansion problem on a grander scale than is the subject here, namely the synthesis of the configuration of an entire new plant, but, as a consequence, they had to assume each expansion was independent of other expansions and of existing equipment as well. For the smaller scale expansions considered here involving the addition of such equipment as reactors, pumps, and distillation columns that affect the whole process when introduced into the flowsheet, the assumption of independence often is not justified. No one has as yet reported how to treat quantitatively this type of chemical process plant expansion problem. Integer linear programming is not applicable because of the extreme nonlinearity of the typical chemical process plant model. Dynamic programming has been applied to chemical plant expansions in the past, but becomes ineffective when information (process) recycle streams are involved. Dynamic programming is also confined to quite small problems due to storage limitations. 3. PROCESSDESCRlPTlON To illustrate the application of the plant expansion algorithm, a model of the RCD Isomax fuel oil desulfurization process (HDS .process) was developed. Expansion of the HDS process over a ten year time horizon was considered. The production schedule of the HDS process was determined from data of the U.S. Bureau of Mines [ 11on the proportion of desulfurized oil products. 3.1 RCD isomax hydrodesulfwization (HDS) process Figure 2 is a flow diagram of the HDS process together with three possible expansions. In the existing process, the feed stream (stream 1) of fuel oil, containing sulfur compounds in the boiling range of trimethylbenzothiophene and dibenzothiophene, is characterized using

103

Optimal expansion of a hydrodesulfurization process

I I

.__-___J

r-l NONLINEAR

0

PROGRAMMING ALGORITHM

r_____-_--_____-____

___-______ ___-----f

_______

I

CHEMICAL

FLOWSHEET

PROCESS

-

SIMULATION

ASSIGNMENT

OBJECTIVE

PROCESS OBJECTIVE

-

PROCESS

FUNCTION

_

CONSTRAINT JACOBIAN

CONSTRAINTS

GRADIENT

FUNCTION

a

VARIABLE

I

THERMODYNAMIC 5 PROPERTIES

I ‘$+lVERBITY

OF

IEXAS

OPTIMIZING

lo I

ELANT

EXPANSION

FGORITHY

Fig. 1. Information flowsheet of capacity expansion algorithm.

a 2

4

3

2

3

4

I

Fig. 2. Hydrodesulfurization process (dashed lines enclose all possible expansions).

D. hf.

104

HIMMELBLAUand

boiling weight fractions. After the feed stream is mixed in equipment No. 3 with make up hydrogen (stream 4) and recycle hydrogen (stream ll), the temperature of the mixture is increased in the reactor preheater (No. 4). The removal of sulfur compounds from the oil occurs in the vapor-liquid catalytic trickle bed reactor (HDS reactor). The product stream of the reactor (stream 7) contains fuel oil with a reduced sulfur content in the liquid phase, and hydrogen and hydrogen sulfide in the vapor phase. Subsequently, the hydrogen and hydrogen sulfide are removed from the fuel oil in a vapor-liquid separator (No. 6). The hydrogen sulfide is removed from the hydrogen recycle stream (stream 10) in a monoethanolamine scrubbing unit (No. 7). The recycle hydrogen is compressed (No. 8) and mixed with the feed stream fuel oil. The thermodynamic and physical property data for the fuel oil, hydrogen, and hydrogen sulfide are listed in Bickel[t]. The desulfurized fuel oil (stream 8) is split in a distillation column into light end products (boiling point 450°F and below) which are sold as product and the bottoms product (boiling point 550°F and above) which is used internally elsewhere in the plant. Figure 2 also shows the possible compatible expan-

T. C.

BICKEL

sions to the HDS process; the expansion projects are enclosed by dashed lines. Expansions 1 and 2 represent two additional trickle bed reactors with support equipment, and expansion 3 is an additional distillation column. Bickel lists all of the sizes of the major pieces of expansion equipment as well as the equipment costs. Table 1 lists all of the process operating variables for one time interval as if all of the expansions existed in the flowsheet. As equipment was removed from ‘the flowsheet, certain variables (and associated constraints) would be correspondingly excluded. For example, if expansion No. 2 did not exist in the time interval being solved, variables 13-24 would be excluded from the nonlinear optimization problem. Removing equipment from the flowsheet posed no hazard to the internal consistency of the problem except in the case in which downstream equipment from a stream splitter was removed from the flowsheet and only one exit stream existed for the stream splitter. The resolution of this difficulty is discussed below subsequent to the development of the process stream splitter constraints. Certain special characteristics of the process model and some of the problems encountered with the model are described next.

Table I. HDS process operating variables Number I 2 3 4 5 6 7 8 9 IO II I2 I3 I4 I5 I6 I7 I8 I9 20 21 22 23 24 25 26 21 28

Variable T5

T6 T7 TS T,, P6 p7 Y1.4 Y1.6 Y1.7 Y2.7 F2 TX5 TM T17 TM T21 PI6 PI7

Y 1.14 Y 1.16 Y 1.17 Y 2.17 62

T25 T26 T27 T28 T3,

:‘o 31 32 33 34 35 36 37 38 39 40

&6 p27

Y 1.24 Y 1.26 Y 1.27 y2.27 F22 T32 62 y3.34

X 4.34

Description Mixer No. 3 exit temperature Reactor No. 5 entrance temperature Reactor No. 5 exit temperature Vapor-liquid separator No. 6 exit temperature Recycle compressor No. 8 exit temperature Reactor No. 5 entrance pressure Reactor No. 5 exit pressure Make up hydrogen to reactor No. 5 Hydrogen feed to reactor No. 5 Exit hydrogen from reactor No. 5 Exit hydrogen sulfide from reactor No. 5 Oil tlowrate to reactor No. 5 Mixer No. IOexit temperature Reactor No. I2 entrance temperature Reactor No. I2 exit temperature Vapor-liquid separator No. I3 exit temperature Recycle compressor No. I5 exit temperature Reactor No. I2 entrance pressure Reactor No. I2 exit pressure Make up hydrogen to reactor No. I2 Hydrogen feed to reactor No. I2 Exit hydrogen from reactor No. I2 Exit hydrogen sulfide from reactor No. I2 Oil flowrate to reactor No. I2 Mixer No. I7 exit temperature Reactor No. I9 entrance temperature Reactor No. I9 exit temperature Vapor-liquid separator No. 20 exit temperature Recycle compressor No. 22 exit temperature Reactor No. I9 entrance pressure Reactor No. I9 exit pressure Make up hydrogen from reactor No. I9 Hydrogen feed to reactor No. 19 Exit hydrogenfrom reactor No. I9 Exit hydrogen sulfide from reactor No. 19, Oil Rowrate to reactor No. I9 Exit temperature of splitter No. 23 Oil flowrate to distillation column No. 25 35O”Fboiling weight fraction flowrate in distillation stream No. 34 450°F boiling weight fraction flowrate in distillate stream No. 34

105

Optimal expansion of a hydrodesulfurization process Table 1. (Confd.) Number 41

X5.34

42

X6.34

43

x3.35

44

X4.35

45

X 5.35

46

x6.35

47 48 49

Description

Variable

21 F36 X3.38

50

X5.38

51

X5.38

52

X6.38

53

x3.39

54

X4.39

55

X5.39

56

X6.39

57

&

550°F boiling weight fraction flowrate in distillate stream No. 34 650°F boiling weight fraction flowrate in distillate stream No. 34 350°F boiling weight fraction flowrate in bottoms stream No. 35 450°F boiling weight fraction flowrate in bottoms stream No. 35 650°F boiling weight fraction flowrate in bottoms stream No. 35 650°F boiling weight fraction flowrate in bottoms stream No. 35 Reflux rate of distillation column No. 25 Oil flowrate of distillation column No. 27 350°F boiling weight fraction tlowrate in distillate stream No. 39 450°F boiling weight fraction Aowrate in distillate stream No. 38 550°F boiling weight fraction flowrate in distillate stream No. 38 650°F boiling weight fraction tlowrate in distillate stream No. 38 350°F boiling weight fraction flowrate in bottoms stream No. 39 450°F boiling weight fraction flowrate in bottoms stream No. 39 55O”Fboiling weight fraction flowrate in bottoms stream No. 39 650°F boiling weight fraction flowrate in bottoms stream No. 39 Reflux rate of distillation column No. 27

3.1.1 Process stream splitters. The mole fractions of the components of the fuel oil are constant throughout the HDS process until the fuel oil is split in the distillation columns. In an effort to keep the dimensionality of the problem as small as possible, the total molar flowrate of the fuel oil was used as a process stream variable, replacing the individual molar flowrates of the fuel oil components. The molar flowrates of the fuel oil components were computed as the total molar flowrate of the fuel oil times the constant mole fraction of the respective component. The thermodynamic properties of the process stream were always calculated from the individual fuel oil components throughout the HDS process. Since the mole fractions of the fuel oil in streams Nos. 2, 12, and 22 are constant parameters, the

average mass =

(

(mass in) - (mass out) = o average mass

(2)

where:

(mass in t mass out)/2 if (mass in t mass out)/2 2 1.0 1.0 if (mass in t mass out)/2 < 1.0.

only equation (model constraint) required for stream splitter No. 1 is an overall mass balance: F,-Fz-F,z-Fzz=O

For example, the value of a mass balance constraint might be of the order of 103,while the value of an energy balance constraint might be of the order of lob. To rectify this difIIculty we scaled the value of each type of constraint to the order of unity, and then used a permultiplier in front of the constraint to provide the required number of digits of accuracy. For mass balance constraints, the scaling technique used was to divide the mass balance by a pseudo average of the mass entering and exiting the equipment

(1)

where F, is the molar flowrate of fuel oil to the HDS process. An energy balance constraint was not required for stream splitter No. 1 since it operates isothermally at 77°F. Numerical difficulties were encountered during the nonlinear optimization phase of the algorithm because of the different magnitudes of the values of the constraints.

For an equality constraint to be satisfied, its absolute value had to be less than 5 x 10e4. To achieve three decimal places of accuracy when the mass balance equality constraints were satisfied, a premultiplier of 0.5 was used in front of all mass balance constraints. Thus, Eq. 1 became O.S(F,- Fz- F,z- F22) maximum 1 o (FI + F2 + Kz + F22 . , 2

(3) =

‘-

Stream splitter No. 23 had an analogous mass balance.

106

D. M. HIMMELBLAU and T. C. BICKEL

In addition, an energy balance was required for stream splitter No. 23 since the HDS reactor exit temperatures were not equal. The energy balance (another model constraint) for stream splitter No. 23 was F&s+ F,sH,8+&Hz- FdL-F~~H,~ =0

(4)

where H8, H18, Hz*, H3*, and H3a were the molar enthalpies of streams Nos. 8, 18, 28, 32, and 36, respectively. The values of the energy balance constraints were scaled similarly to the mass balance constraints with average energy defined analogously to average mass. For Eq. (4)

Two decimal places of accuracy were desired for the values of the energy balance constraints (since the accuracy of the process stream enthalpies was approximately two decimal places), and, as a result, a premultiplier of 0.05 was used in front of all energy balance constraints. During the testing of the HDS process model, when expansion project 3 (distillation column No. 27) was excluded from the flowsheet, the resulting nonlinear optimization problem became degenerate. The degeneracy occurred because stream splitter No. 23 had only one exit stream (process stream No. 32) and, as a result, the mass balances for splitters Nos. 1 and 23 were identical (since Fl and FS2 had to be equal from an overall mass balance on the HDS process). The degeneracy was overcome by excluding variable F3* and the mass balance for splitter No. 23 from the problem. Subsequently, prior to optimization, all stream splitters were tested to assure that they had more than one exit stream. If a stream splitter had only one exit stream, the stream splitter mass balance and mole fraction constraints were excluded from the problem, and through user supplied input data, the variables and constants associated with the exit stream were replaced (reducing the dimensionality of the problem) with variables and constants associated with other process streams. If a stream splitter had only one entrance stream and had an energy balance constraint, the exit stream temperatures were set to the entrance stream temperature and the energy balance was excluded. 3.1.2 Process stream mixers. A hydrogen mass balance and an energy balance were used to model stream mixers Nos. 3, 10, and 17. The mixer exit stream flowrate of the fuel oil was replaced by the mixer entrance stream flowrate of the fuel oil to avoid introducing a fuel oil mass balance. The hydrogen mass balance for stream mixer No. 3 was

WY,,+ maxhum

( 1

Y,.,, - Y,,d Yl,rl+ Y1.5) ., 2 o

(Y1.4+

=O*

>

and the balances for mixers Nos. 10 and 17 were similar. 3.1.3 Reactor preheaters. No mass balances were used to model the reactor preheaters Nos. 4, 11 and 18. The exit stream molar flowrates of all components (hydrogen and fuel oil) were replaced by the entrance stream molar flowrates, eliminating the need for mass balance constraints. The energy balance constraints were solved explicitly for the heat absorbed by the preheaters, Q, and thus did not have to be included as constraints in the model. 3.1.4 HDS reactors. The catalyst beds of each of the HDS reactors Nos. 5, 12, and 19 were made up of l/8 in diameter cobalt molydate pellets with a bed bulk density

of 441b/ft3. The HDS kinetics to calculate the H2S formed for this type of catalyst were developed by Frye 8z Mosby [3] for an isothermal reactor using a LangmuirHinselwood concept of a surface controlling reaction. If the temperature of the reactor was assumed to be the average of the entrance and exit stream temperatures and the pressure of the reactor was assumed to be the average of the entrance and exit stream pressures, the results of an example adiabatic flash calculation showed that 99% of the hydrogen was in the vapor phase, 94% of the hydrogen sulfide was in the vapor phase and 66% of the fuel oil was in the liquid phase. If the vapor phase components of the fuel oil had been made variables in the problem, it would have increased the dimensionality of the problem by 30 variables. Because the size of the nonlinear problem was already approaching the capacity of the CDC-6600 used for the computations, it was assumed that all the fuel oil remained in the liquid phase. Scaled component mass balances on the hydrogen and hydrogen sulfide, a scaled energy balance, and a pressure drop constraint were used to model each HDS reactor. The equation for the reactor kinetics could not be used in determining the consumption of hydrogen in the reactor because some cracking of the fuel oil occurs in the reactor, a factor which was not taken into account in developing the reaction kinetics. The consumption of hydrogen in the HDS reactors was estimated from the data of Schuman 8t Shaht[41 to be 135 SCFHJbbl fuel oil. Heat generation by the HDS reaction was also that estimated from the data of Shuman and Shalit. The pressure drop across the HDS reactor packed bed was calculated by combining the Ergun equation[5], the correlation developed by Sato et a1.[6] for two phase flow through a packed bed, and the viscous friction factor formula of Tallmadge [7]:

(6) where AP = total pressure drop across packed bed; APv, APL = pressure drop across packed bed if vapor and

and the other two balances were similar. The energy balance for process mixer No. 3 was

O.O~(Y,,~H~+FSH~+ Y,.,, HI,-FsH5 maximum 10 tY,.4H4+F,%+ Y,,IIHII+ FsH,) =O . , 2

(7)

107

Optimal expansion of a hydrodesulfurization process

liquid flowed separately; A = (APJAPv)“L. The scaled pressure drop for reactor No. 5, for example, was thus (9)

A premultiplier of 0.01 was used in front of the constraint because the accuracy of the pressure drop correlation was approximately two decimal places. Similar equations were used for reactors 12 and 19. 3.1.5 Vapor liquid separators. Only energy balances were needed for the vapor liquid separators Nos. 6, 13, and 20 as the gas was hydrogen and hydrogen sulfide while the liquid was fuel oil. To mix process streams Nos. 8, 18, and 28 the pressures of the exit streams of the vapor liquid separators were set to 300 psia. The reduction in pressure in the separators caused a reduction of the temperature of the exit streams of the separator, and, as a result, a scaled energy constraint was required for each separator similar to Eq. 7. 3.1.6 H2S scrubbers. The amine scrubbers were in the flowsheet to remove the hydrogen sulfide from the recycle hydrogen (since H,S inhibits the HDS reaction). It was assumed that all the H2S was removed in the scrubber, and, as a result, no mass balance or energy balance constraints were required. Also, no operating costs were assumed for the scrubbers since the H2S removed in the units was small, and the rich liquor could be used elsewhere in the chemical plant without treatment. A fixed capital cost was incurred for each scrubbing unit itself. 3.1.7 Hydrogen recycle compressors. No mass balances were used to model compressors Nos. 8, 15, and 22, but an energy balance was required because the pressure of the exit stream was increased over the entrance stream pressure. To keep the dimensionality of the nonlinear optimization problem as small as possible, the compression was assumed to be isentropic. The scaled balance for compressor 8, for example, was (1o)

0.05 [T,, - (Z’,ot 459.69)(P,,/P,o)k-“k -459.691= o

the saturated liquid bottoms product stream temperature was determined to be 641°F. The heat duty of the overhead condenser was determined from the difference of the enthalpies of the saturated distillate vapor at its dewpoint temperature and the saturated liquid distillate at 536°F to be 17,27OBtu/lb mole while the heat duty of the reboiler was determined from the difference of the enthalpies of the saturated bottoms vapor at its dewpoint and the saturated liquid bottoms at 641°F to be 19,825 Btu/lb mole. Because the conditions of the exit streams of the distillation column were specified and the values of the condenser and reboilers duties determined, no energy balance was required for the distillation columns. The mass balance distillation columns No. 25 was 0.5(X.w- x.35 - X.36) maximum 10 Wi,34+X.35+X,36) . , 2 (

=O

i=

3"

* * '6

(11)

I

and a similar equation was used for column No. 27. The reflux rate constraint for distillation column No. 25 was: 001

I,,_(~~D)m-G

=o

(12)

G-l

.

where (L/O), is the minimum reflux rate and G = ((L/D) - (WD),)/((UD) - 1). A similar equation was used for column No. 27. Equation (12) did not have to be scaled, but a premultiplier of 0.01 was used to obtain two decimal places of accuracy in the value of the constraints. All the equations of the HDS process have now been described; the next section deals with the inequality constraints. 3.2 HDS process inequality constraints Three types of inequality constraints were used in the modeling of the HDS process. The first type guaranteed removal of sulfur to at least 0.8 Wt/% of fuel oil such as inequality (13) for reactor No. 5:

’ (13) A premultiplier of 0.05 was used in order to obtain two decimal places of accuracy in the value of the constraints. 3.1.8 Distihztion columns. In the distillation columns, the fuel oil components were used as the stream variables (rather than the total molar flowrate of the oil) so that component mass balances for the fuel oil were required. The feed streams to the distillation columns comprised a saturated liquid at 4Opsia, and the conditions of the distillation feed streams were set in order to fix the feed tray location. From a bubble point calculation, the temperature of the distillation column feed streams was determined to be 562°F. The pressure of the fuel oil was reduced from 300 psia to 40 psia in a pressure reduction valve (equipment Nos. 24 and 26). Streams Nos. 32 and 36 were either heated or cooled to 562PF after the pressure reduction. The overhead condenser on the distillation columns was a total condenser in which the overhead distillate vapor was condensed to saturated liquid. From a bubble point temperature calculation, the overhead product stream temperature was determined to be 536°F. From a bubble point calculation,

The second type of inequality constraint guaranteed the purity of the distillate stream products such as inequality (14) for distillation column No. 25: x3.34 o.95

-

( x3.34+x4.34+

+

x4.34 X5.34-tX6.34

(14)

I So.

The third inequality constraint guaranteed the distillate product molar flowrate of both distillation columns to be at least one half of the feed stream molar flowrate:

-

(~3,38+~4.38+~5.38+~6.38)


(15)

4. OBJECTIVE FUNCTION

Among the many alternatives that might be chosen for the objective function, a simple criterion was selected, namely the maximization of the discounted profits over

108

D. M. HIMMELBLAU and T. C. BICKEL

the entire time horizon. The profit of the chemical plant can be expressed as: Objective function = S - CF - Co.

each expansion interval and each piece of equipment or stream. Specific cost data for the HDS plant are listed in Bickel.

(16)

specifically:

5. THE BRANCH AND BOUND PROCEDURE time horizon over which the expansions are to be considered is discretized into p expansion intervals. /3, is the capacity expansion interval in which project t is first utilized. The construction of expansion t is assumed to be started and completed prior to interval p,, but fixed and operating charges do not commence until pt. Let (I~ be the interval after which project t can be utilized, and let yt be the interval such that project t must be utilized (built) during or before interval ‘y,. The variable PI, its lower bound a,, and its upper bound yt, are all integers: The

where S, is the sales revenue in time interval k, k = 1, . . . , p, s(k,r) = l/(1 t r)’ is the present worth factor, k is the expansion interval index, and r is the modified discount rate, i.e. an interval rate, not necessarily an annual rate. Sales revenues are computed from the product of the production rate and the sales prices of the respective products. Prices are assumed to be increased periodically. It is also assumed that demand is growing in each period, and that production is scheduled so as to meet the specified demand. Expansion is carried out at appropriate times so that demand is met in every time period. The fixed capital cost term, C,, in the objective function is the present value of total fixed charges of all the expansion projects. Before /11,the (integer) interval at which the project t can first be utilized, no costs are incurred by project f because it does not yet exist. At the start of expansion interval &, the capital cost of project t is represented as &C,. The discounted capital cost attributable to project t over the time horizon considered is:

where C, = the capital cost of project t; 4t = capital recovery factor, the fraction of C, charged per interval; s&)=Oif k<& and 1 if krp,. Thus, the total fixed charges of all the expansion I’ projects will be a function of /3.

The operating costs of the chemical plant, Co, can be expressed as the sum of the discounted costs of the raw materials and the plant costs for utilities, etc. for the initially existing plant (the first term on the right hand side of Bq. 19) plus the operating costs for the expansion projects implemented (the second term on the right hand side of Eq. 19). Co(B,T,PX,Y,Z) = = $, &(T,P,X,Y,Z)s(k,r) + 2 5 D:(T,P,X,Y,Z)lr(k,r) ,=I k-l

(19)

where I&’ = the operating cost of expansion project t during interval k, if it were to exist in time interval k (I& is zero otherwise); Bk = the operating cost of equipment existing at time t = 0, during interval k. Equation (16), after introduction of Eqs. (17)-(19), is the objective function for the capacity expansion problem. Note that the objective function involves the vector of discrete variables fl and the continuous variables T,P,X,Y,Z, all of which have upper and lower bounds for

01a,
t=l,.*.,

n

(20)

where n is the number of possible capacity expansion projects. The upper bound of p t 1 rather than p permits an expansion project to be excluded from being utilized, i.e. not built. The strategy of branch and bound optimization is simply the coordinated search over all possible combinations of the discrete variables carried out in such a manner as to eliminate certain solutions to the discrete problem. Thus, an optimal solution to a problem can be obtained with less that a complete enumeration of all the possible solutions. A number of branch and bound algorithms have been proposed to solve a wide variety of combinatorial problems (Golamb t Baumert[8], Little et al. [9]). Hillier & Lieberman[lO] provide a brief, simple summary of the technique. Each individual problem, however, requires specialized features. In essence, one repeatedly partitions the discrete problem solution space into smaller and smaller subspaces. A lower bound in the case of minimization, or an upper bound in the case of maximization, is computed for the value of the objective function for each subspace. A problem can be represented as a decision tree which illustrates the partitioning of the solution space. The tree consists of vertices, denoted by circles, and edges which connect the vertices. Each vertex represents a set of constraints associated with the objective function of the capacity expansion problem. Each edge represents one additional integer constraint added to the constraint set of the preceeding vertex. Thus, the edges of the enumeration tree represent the addition of increasingly restrictive construction date bounds on the expansion projects. Initially (vertex 1), the bounds on the construction dates are: a, =o

t=1,...,n

yt=ptl

t=1 ,...,

II.

Either two edges, called a partition, or no edges will issue from a vertex. When a vertex originates two edges, the left branch below the vertex increases the lower bound on the construction date of an expansion project, while the right branch below the vertex decreases the upper bound on the construction date of an expansion project. In the branch and bound strategy, one systematically applies more restrictive upper and lower bounds to the expansion project construction dates, /3, ... untu:

Optimal expansion of a hydrodesulfurization process

109

(2) The upper bound to the capacity expansion objective function at the vertex is less than the best feasible solution found so far in the search. at which stage the branch is completed. (3) The relaxed problem at the node is found to be A partition cannot occur from a vertex where (Y(t 1 = ‘yt,r = 1;. . . , R because doing so would violate the con- infeasible, i.e. the condition that “production must meet demand” cannot be satisfied. dition that a < y. The object of the branch and bound When one of the three above criteria occur, the search procedure is not to solve a problem at each vertex of the entire tree, but to exclude sections of the tree prior to the backs up (“backtracks”) to continue from the first precondition that (Yt 1 = y. This objective is achieved by decessor vertex from which only one branch has been the calculation of an upper bound to the capacity expansearched. When all vertices in the tree have been fathomed using the above criteria, then the best feasible sion objective function at the vertices in the tree. solution found during the search is the optimal solution The conclusion that a partition should not be introduced at a particular vertex, t)i, (called fathoming a to the discrete problem. The computer code to execute the entire optimization vertex) is based upon the fact that all problems below t+ in the decision tree are restrictions of the problem at Vi. is documented in Bickel[2]. Therefore, an upper bound to the value of the capacity expansion objective function at ai is alsoan upper bound 6.RESULTSOF THE CAPACITYEXPANSIONOF THE HYDRO. DESIJLFURIZATlONPROCESS to the value of the objective function of all problem The HDS process model developed above was solutions found below ui. If the upper bound found at vi is less than the best value of the objective function of the analyzed over a 10yr time horizon from the viewpoint of capacity expansion problem found at another vertex in possible expansions. With the time horizon discretized the tree where a t 1 = y (i.e. at a feasible solution to the into four equal expansion intervals, and with three possible expansions of the HDS process, there existed 125 capacity expansion problem), then vertex tri is fathomed, since no solutions exist to the capacity expansion prob- discrete solutions to the capacity expansion problem. lem below ui that are better than the best feasible solu- With all the expansions included for each expansion interval, the nonlinear optimization problem was comtion found so far in the tree. prised of 57 variables, 16 linear equality constraints, 21 The upper bound to the capacity expansion objective nonlinear equality constraints, 6 linear inequality confunction is obtained at a vertex by solving a nonlinear programming problem which is a relaxation of the straints and 114 bounds on the variables. After some capacity expansion problem described in Sets. 3 and 4. experimentation to select a robust nonlinear programThe relaxation employed is in the terms Co and CF. In Co ming code, it was decided to employ the Generalized we let Reduced Gradient (GRG) code of Abadie & Guigou [ 11,121 modified so as to employ numerical 0:’ = 0 if k 5 0, derivatives. (22) Figure 3 shows the branch and bound enumeration tree Dz[’= Dkl if k > a,. generated during the solution of the capacity expansion problem. Only 47 problems (vertices) in the enumeration What Eq. (22) means is that during the interval (Y,to /3,, tree were solved, i.e. a relaxation of the discrete capacity expansion project t is allowed to operate even though it expansion problem was solved 47 times using the nondoes not yet exist. In the term C,, the fixed capital linear programming algorithm. Eight vertices were charges for the expansion equipment in the objective fathomed prior to starting the solution of the relaxation function are modified by substituting S*(&) for S(/%), of the discrete capacity expansion problem because the where: constraints were not satisfied. These vertices are shown in Fig. 3 as unnumbered vertices. Eighteen vertices were fathomed using the upper bound fathoming criterion. The starting point at vertex 1 of the enumeration tree was feasible with respect to the HDS process modeling Equation (23) means that no fixed capital costs are constraints. A starting point was determined for each charged for expansion project t until after the upper expansion interval by: bound on the construction date on project t. The proof (1) Setting the molar flowrate of fuel oil to each HDS that this relaxation does indeed yield an upper bound to reactor at a value such that the mass balance constraint the objective function has been demonstrated by of stream splitter 1 was satisfied. Bickel[2]. (2) Setting the inlet temperature and pressure of each. Based upon the value of the objective function at a HDS reactor and support equipment (the molar flowrate vertex, the decision is made either to fathom the vertex, into each reactor had been previously set in step l), and or to generate two new edges. The generation of two new solving the constraints associated with each reactor and edges requires the selection of the expansion project for support equipment simultaneously. which to modify the construction date bounds, and the (3) Sending 50% of the fuel oil to each distillation expansion interval on which to partition. The selection column from stream splitter No. 23; the exit temperature criteria for these choices are too extensive to describe of the splitter was calculated interactively by assuming an here but are explained in Bickel. exit temperature, and evaluating Eq. (5). This procedure The search for an improved value of the objective was repeated until Eq. (5) was satisfied. function of the capacity expansion problem continues (4) Satisfying the component mole balances on each down a branch of the tree until a vertex is fathomed, i.e. distillation column until Eq. (4) and its equivalent for no edges are generated from the vertex. The criteria for column No. 27 were satisfied, and then using Underfathoming are summarized below: wood’s technique to determine the reflux ratio of each (1) a, + 1 = .y,, t = 1,. . .) n distillation column. a,+l=y,

t=1,...,n

110

D. M. HIMMELBLAU and T. C. o=Io,o,o

BICKEL

1

V=lS,S,Sl 1

A

v3= 1

2 v,=r

a2=4

d -

20 v,=i

a3=4 V3’4

21 v2=2

a2=2 22

24 v2=3

a2=3 f\

23 -

-

-

23 v2=2

a2=2 26 V3’2

a3=2 21

32 v3=3

a3= 3 8

a,=3

-

30

v,=s

-

29 -

Fig. 3. HDS processexpansiontree (the underscoreon the circle designatesa fathomedvertex).

As the capacity expansion problem became more restricted (proceeding down the tree), the starting point for the non-linear optimization problem at a vertex became infeasible because of the removal of expansion equipment from the HDS process flowsheet. In an effort to maintain as close to a feasible starting point as possible to save computer time, an option in the code allowed the addition of the component molar flowrates from one process stream to another. For example, if expansion No. 3 did not exist during a particular expansion interval, the starting point for the nonlinear optimization algorithm for the expansion interval would be modified in that the molar flowrate of process stream No. 37 would

be added to process stream No. 32 and all subsequent streams. At the conclusion of the calculations for the example problem, the solution was found to be /I = (4,1,5), i.e. expansion No. 1 is built in expansion interval four, while expansion No. 2 is built in interval one. Expansion No. 3 is never built. The value of the capacity expansion objective function (excluding charges for equipment in place at the initial time) at the optimal solution was $98.94~ 106, a value not much higher than that of most of the feasible fathomed vertices. While the relative difference between the values of the objective function is small, the absolute difference

111

Optimalexpansion of a hydrodesulfmizationprocess

exceeds one million dollars. Of course, many nodes were fathomed using the upper bound criterion prior to reaching a unique solution to the capacity expansion problem, and the value of the capacity expansion objective function at such nodes could be considerably less than $97.91x lo”. The small relative differences in the values of the objective function are representative of the low growth in utilities costs used in the simulation. The computation time for the HDS process expansion on a WC-6600 computer amounted to 9.7 hr. Figure 4 presents the percent of computation time spent in the major phases of the capacity expansion calculations. Clearly, the largest consumption of time was spent in the calculation of the thermodynamic properties. While this does not necessarily indicate that the CADC CONCEPT THEKMOPAK computer routines employed were inefficient, it does indicate that any small improvement in the calculation time of the thermodynamic properties would decrease the computation time of the overall algorithm. Table 2 lists the major elements of computer time for the various individual calculations performed during the solution of the capacity expansion problem. Clearly, the evaluations of process stream enthalpies and phase equilibria K-values were the most time consuming together with the inversion of the Jacobian matrix of the constraints in the nonlinear programming algorithm. The large time required for the calculation of process stream mole fractions is not surprising considering the number of times the mole fractions are needed, such as in the calculation of the enthalpies and K-values. While the overall capacity expansion algorithm is time consuming, the results obtained by the algorithm need to be weighed against the results achieved by the“rules of thumb” calculations currently in use today or by case studies. 7. DISCUSSIONAND CONCLUSIONS

The solution procedure developed in this work can

accommodate any number of discretized time intervals and any number of possible equipment expansions keeping in mind that maximum number of problems that have to be solved is given by the expression (p t 1)“. Some fraction of the maximum number will actually be solved. Any number of process operating variables can be accommodated, limited only by the ability of the nonlinear programming algorithm to solve the nonlinear optimization problem. The process streams in the example were

1

I 4

Fig. 4. Histogramof computationtime for the capacity expansion algorithm.1= Branch and bound algorithm,input, output, storage; 2 = Chemical process simulation; 3 = Thermodynamic calculations;4 = Nonlinearprogramming alborithm;5 = CDC6600 system routines. limited to liquid and vapor phases because of the inability of the CADC CONCEPT THEKMOPAK

package[l3] to calculate phase equilibria values for streams containing solids. The process operating variables are limited to continuous variables. The algorithm can treat processes with multiple recycle and bypass streams. Most equipment indigenous to chemical processes can be incorporated into the flowsheet, but the kinetics of all reactions occurring in the process must be supplied by the user. Sizing the expansions is made possible by allowing multiple expansions of different sizes comprised of individual pieces of equipment with the overall result being the optimal size of a specific expansion. Obsolescence of equipment can be accommodated by examining the values of the continuous variables at an optimal solution. If they were zero for all intervals subsequent to interval k, then the equipment no longer interacted with the rest of the process, and could be removed. (This condition did not occur for the HDS process.) In contrast to decisions based on rules of thumb, the results of the algorithm are quantitative. In the example of the expansion of the HDS plant, the higher capital cost expansion project No. 2 (HDS reactor No. 19) was implemented before the lower capital cost expansion project No. 1 (HDS reactor No. 12) because of the difference in operating costs (economies of scale). If the lower cost expansion project No. 2 (such as at vertex

Table 2. Majoruses of computertime of the capacity expansionalgorithm

Computeroperation Evaluationof phase equilibriaK-values Evaluationof process streammolarenthalpies Determinationof process streammole fractions Inversionof constraintJacobianmatrix Evaluationof chemicalprocess constraints Evaluationof constraintJacobianmatrix Evaluationof capacity expansionobjective function Updatingprocess streampropertiesduringevaluation of numericalderivatives Evaluationof reducedgradientand reentering constraintfeasible region Evaluationof gradientof capacity expansion objective function CACE

Vol. 4. No. 2-C

5

Percent computer time

Histogram section Fig. 4

23.4 17.3 14.2 11.5 13.4 1.1 1.0 1.8

; 2 4 2 2 2 2

1.8

2

0.1

2

D. M. HIMMELBLMI and T. C. BICKEL

112

No. 23 of the enumerative tree shown in Fig. 3), a reduction in profit of nearly one-half million dollars would have occurred. Naturally, the proposed algorithm has some limitations. First the cost of computation time to carry out the analysis (approximately $2200 for the expansion of the HDS plant using a CDC-6600 computer at a cost of $230/hr of CPU time) may seem high, but most likely the computation cost will be a very minor element in the overall cost of the design, particularly the cost relative to the cost of collecting and validating the data for the calculation of the thermodynamic properties that may run into tens of thousands of dollars. The cost of computer time is even more insignificant in comparison with the savings that can be achieved via the optimization procedure compared with a case study approach to the same problem. If one assumes that the sizes of the nonlinear optimization problems solved during the branch and bound enumeration is represented by a binomial distribution, an approximation to the computation time, T, in hours for the capacity expansion algorithm will be:

7 = 0.2 z0 &

eNi’s

where Ni equals the number of elements in the Jacobian matrix of the constraints when n expansions exist in the process diagram. The results for the HDS plant indicate that the capacity expansion objective function was “flat”, i.e. relatively minor changes in the value of the objective function occur between single solutions to the discrete problem. If this situation were true in other potential problems, a near optimal solution could be achieved to the discrete problem much quicker by modifying the upper bound fathoming criterion, in that a vertex would be fathomed if its upper bound were within, say, one half of 1% of the best feasible solution found so far to the discrete problem. A second limitation exists because of the stochastic nature of some of the data needed for the optimization. An engineer would have to reevaluate future expansion policy as the data change. Third, to run the capacity expansion algorithm, an elaborate model of the chemical plant must be developed. On the other hand, once the model has been developed (as long as similar equipment is added to the plant) modification of the model would be trivial. In any case, the model would probably be necessary in the design and operation of the plant. Finally, the expansions to a chemical plant must be prespecified, and special care must be used in the choice of possible expansions. Too few expansions may yield a suboptimal result, while too many expansions greatly increase the computation time for the algorithm. NOMENCLATURE

CF present value of total fixed charges of all expansion projects

Gl present value of operating costs of existing plant plus expansion implements capital cost of project t operating cost of expansion project t during interval k if it was implemented 0:’ relaxed operating cost given by Eq. (22) Ek operating cost of equipment existing at I = 0 during interval k

2

total molar flowrate of fuel oil (subscript denotes stream) factor in Gilliland correlation used Eq. (12) molar enthalpy of stream designated by subscript heat capacity ratio (C,/C,) = 1.288 reflux ratio in distillation columns minimum reflux rate number of possible capacity expansion projects number of elements in the Jacobian matrix of the constraints number of capacity expansion projects pressure of stream designated by subscript total pressure drop across a packed bed pressure drop across a packed bed with liquid flowing pressure drop across a packed bed with vapor flowing heat absorbed modified discount rate (interval rate) present value of sales revenue for all products sales revenue in time interval k index for project number temperature of stream designated by subscript vertex i in decision tree molar flowrate for liquid component (subscripts designate component and stream) molar flowrate for gaseous component (subscripts designate component and stream) equipment variable, in general (design variable specific to one piece of equipment) reflux rate of distillation column No. 25 reflux rate of distillation column No. 27 Greek symbols

interval after which project t can be utilized integer valued interval at which project t can first be utilized interval such that project t must be utilized before

interval y, _ integer factor for equipment utilization modified factor for equipment utilization given by Eq. (23) present worth factor computation time in hours capital recovery factor, i.e. fraction of C, charged per interval

I. U.S. Departmentof the Interior, Bureau of Mines Mineral Industry Surveys, Crude Petroleum, Petroleum Products. and Natural-Gas Liquids, Pi&burg, Pennsylvania (1977).

2. T. C. Bike]. The o&ma1 cauacitv exoansion of a chemical plant via n&linear ;nteger p~ogra~mmi~g. Ph.D. Dissertation, The University of Texas at Austin, Austin, Texas (1978). 3. C. G. Frye and J. F. Mosby, Chem. Engng Prog., 63, 66 (1%7). 4. S. C. Shuman and H. Shalit, Catal. Reu. 4, 245 (1970). 5. S. Ergun, Chem. Engng Prog., 48, 89 (1952). 6. Y. Sate, T. Hirose, F. Takahashi, and M. Toda, J. Chem. Engng Japan 6, 147 (1973). 7. J. A. Tallmadge, A.LCh.E.J. 16, 1092(1970). 8. R. E. Golamb and L. D. Baument, J. Assoc. Comput. Much. 12,516 (I%$

9. J. D. C. Little, K. G. Murty, D. W. Sweeney, and C. Karel, Op. Res. 11,972 (1%3). IO. F. S. Hillier and G. J. Lieberman, Introduction to Operations Research, Holden-Day, San Francisco, California, (1967). pp. 565-570. 11. J. Abadie and J. Guigou, Gradient Reduit Generalise, Note E.D.F. HIO69/2 du Avril, Electricitt de France, Paris (1%9). 12. I. Guigou, Presentation et Utilisation du Code GREG, Note HI 582/2,25 Mai, Electric& de France, Paris (1971). 13. Computer Aided Design Centre, CONCEPT Mark III User’s Manual, Cambridge, England (1973).