- Email: [email protected]
An approach to the optimal design of networks
Chemical Engineering Science, 1972, Vol. 21, pp. 113 1-I 141,
Pergamon
Press.
Printed
in Great
Britain
An approach to the optimal design of networks B. A. MURTAGH Computer Sciences Department, Fluor Corporation, Los Angeles, U.S.A. (First received 22 November 1970; accepted 9 September 197 1) Abstract-The design problem considered in this paper is the optimal distribution of pressure drop over a network in which flow-rates are specified, but pipe diameters may vary. By using elements of the duality theory of nonlinear programming, the problem is cast in a form which is computationally attractive. The approach presented is of general applicability to networks in which the equations of flow are nonlinear. 1. INTRODUCTION
has been much work recently on the application of mathematical programming to problems which arise in determining the optimal distribution of flows in networks[l, 21. The algorithms obtain optimal flow patterns by minimizing some cost function, while satisfying the laws of conservation around the junctions of the network. The problem treated in this paper is not the usual type of network flow problem in that the flows throughout the network are considered to be fixed. The maximum total drop in potential energy in each path from source to sink is also specified. The problem is to adjust the resistance in each arc of a path to achieve a minimum total cost of the network, considered to be a function of resistance. The treatment which is given is directed to a specific problem which frequently arises in chemical engineering design: one of compressible fluid flow through a network of pipes. However, there are many analogous problems which arise in other contexts - some generalizations will be discussed in Section 7. The optimization problem which arises enjoys a particularly simple structure, which is exploited by formulating a second (“dual”) problem, the solution of which also solves the original (“primal”) problem. The concept of duality in nonlinear programming has often been passed off as a mathematical nicety of little practical value. (For example, Wolfe[3] remarked in a recent review that, “I have yet THERE
to see a case in which the dual problem is truly easier to solve than the primal”.) However, in this particular instance, the dual problem is not only theoretically pleasing, but also much more easily solved, as computer results show. 2. A CHEMICAL RELIEF
ENGINEERING PROBLEM: HEADER DESIGN
Consider a network of pipes as shown in Fig. 1. These pipes lead from relief valves on various refinery vessels to a discharge zone, or flare. A pipe section is defined as the length of pipe between two junctions or between a valve or the flare and a junction. A path is defined as the sequence of pipe sections between any one valve and the flare. The configuration of the pipe network is assumed to be already determined by the layout of the refinery vessels. Thus the lengths of all the pipe sections are specified and the only design variables remaining are their diameters. In an emergency, the vessels must discharge their contents through the relief valve at maximum (sonic) velocity for an initial transitory period. The problem is to select the pipe diameters in each path from valve to flare so that the total pressure drop over the path remains within a pre-specified limit. Assuming that the flare pressure is specified (usually slightly above atmospheric pressure), the limit is imposed so that the back-pressure in the pipe at the valve is below the critical value, thus ensuring sonic velocity in the throat of the valve.
1131
B. A. MURTAGH A
t
wheref=
Flare
friction factor = 0.0475 (Rc+-O.‘~~
+---I I5
13
=
o.a475
(This gives factor.)
14
(3)
(!E3)-0’186
a conservative
value
of friction
and
Valve
I
2
3
4
Fig. 1. Schematic representation
5
6
7
6
of a piping network.
There are obviously many feasible solutions to this problem. Note also that some pipe sections are common to many paths, so the problem cannot be solved for each path individually. In formulating this as an optimization problem, we impose all the pressure drop limitations as a set of constraints, and choose as the ‘best’ feasible solution that set of pipe diameters which minimizes the total weight of pipe. The weight per unit length of pipe is assumed to be a linear function of diameter; (this is a valid assumption when the diameter is much larger than the thickness): Wi= (Y+ pdi lb/ft.
(1)
where di is the diameter in inches and the subscript i refers to the number of the pipe section. The coefficients a and fl are obtained by linear regression. Following API-RP520 design standards [4], the appropriate equation is for isothermal flow of a compressible fluid. This equation takes the form [5] : Pz_P2=64.f.L.W2.R.T. A B r2.gc.M.D5
x log10EpA
(2)
PA = Ps = D = L= W= R = T= M = gc = p =
inlet pressure (abs) (lbF/ft2) outlet pressure (abs) (lb,/ft2) dia. (ft) length of pipe (ft) weight rate of flow (lblsec) gas constant = 1546 (ft-1bJlb mole “R) absolute temperature (“R) molecular weight (lb/lb mole) 32-17 (lb ft/lbF sec2) viscosity (centipoise).
For any section of pipe, the first three terms above vary, while those below are constant. The properties at the valve are specified, and the properties downstream are calculated by simple mixing rules. (There are many assumptions worthy of contention used here, especially that of isothermal flow. However, we will not take issue on these since they are either set forth by API-RP520, the usual standard in contracts, or used in practice with safety as the justification.) The pressure drop over each section of pipe is given by an equation of the form of (2). Representing all constant terms by one constant, Ki, for each pipe section, i, and leaving the presence of the subscript as understood, Eq. (2) with (3) becomes: P,z_P:_I=gg where C, represents the term in brackets in Eq. (2), the kinetic energy correction factor, and Pi and Pi-1 are the inlet and outlet pressures respectively. Note that Ci is a function of Pi and Pi_l; we will, however, ignore this fact
1132
An approach to the optimal design of networks
for the present and treat it as a constant. (This point will be discussed later.) Consider r pipes in series numbered from 1 to r consecutively; the individual pressure drops are calculated by recursion:
(5) P,“-P,2_,=f-g$.
r
Hence, summing:
boundaries on the feasible region in which the point x may lie. A constraint j is deemed to be active if the current point satisfies (8) exactly, so that the equality sign holds. By the introduction of some convenient notation, it is possible to express the problem described in the previous section in the form of (7) and (8). Let n be the number of pipe sections and m be the number of paths. Assign a unique number between 1 and n to each pipe section, and a unique number between 1 and m to each valve (see Fig. 1). It is evident from Fig. 1 that this last set of numbers may be used to refer to both the valve itself and the path from the valve to the flare since there is a one-to-one correspondence. Let the subscript i represent the number of the pipe section, and the subscript j represent the number of the path. Let Sj be the set of pipe section numbers in path j. Then, for example, in Fig. 1:
where P,. and P, are the inlet and exit pressures respectively. The problem is to adjust each dia., diy to minimize the sum of LiWi over all pipe sections, i, while obeying pressure drop limitations of the form given by Eq. (6).
S, = (1,9,13,15) S, = (2,9,13,15)
SS = (8, 12, 14, 15). 3. OPTIMIZATION:
PRIMAL
FORMULATION
The objective function may be expressed simply as the sum of the weights of each pipe section: From Eq. (1):
Linearly constrained nonlinear optimization problems may be expressed in the form: minimizef(x)
= _&x1---x,)
subject to 5 a x. d bi i=l 15t
j=
(7)
l-_-m
.f(dl ---dn) = (8)
where n is the number of variables, and m is the number of constraints. Ax) represents the objective function, which usually expresses some form of cost associated with the values of the independent variables, xt, i = 1 ---n. ati is the coefficient of variable i in constraint j, and b, is the corresponding constant term. Geometrically, the equations
represent
$, au.Xi-bj=O
j=
n-dimensional
hyperplanes
l---m
(8a)
which form
5 &(a
i=l
+ pdi).
(9)
The constraints may also be written down easily with the aid of the above notation. There will be m pressure drop equations of the form of Eq. (6), one for each path j, j = 1 ---m. They may be expressed as inequality constraints by using the maximum allowable inlet pressure in place of the actual inlet pressure: Let bJ = Pj” spec. - PO2
j=
l---m
(10)
where Pj spec. is the specified maximum allowable back pressure at valve j, and P,, is the flare pressure.
1133
B. A. MURTAGH
Let ati =
I
&C,, 0
i E si ,
i= l---n,
l---m
i$Ssj.
(11)
In order to make the constraints introduce the transformation: xi=
j=
l/di4’R14 i= l-_-n.
(12)
minimize f(x) = E L*((Yf /3x,-0.20”)
(13)
i=l
$l a+~~ 5 bj
j=
4. REVIEW
linear, we
We then have the problem:
subject to
the dimensionality from n to m.
1 ---m
(14)
which is the desired form. Gradient projection methods of nonlinear programming [6] are particularly well suited to problems posed in this form. There are also auxiliary constraints in the form of lower bounds on the pipe diameters: di 2 dimin. > 0 + xi 5 ximal. < ~4 i= 1 ---II
of the optimization
OF DUALITY
problem
THEORY
In this section we shall briefly describe those essential elements of the duality theory of nonlinear programming which are required. Duality theory has been put forward in various guises, but the most easily visualized concept is that of a saddle point of an associated Lagrangian function. The prima1 problem involves minimizing in one space of variables, while the dual problem involves maximizing in another space of variables. Under certain regularity assumptions, the two solutions are equal and correspond to a saddle point (not necessarily a stationary point) of the Lagrangian in the combined space of both variables. Consider the problem expressed by (7) and (8). Using Lagrange multipliers, hj, to adjoin the constraint equations to the objective function we obtain a Lagrangian function, which takes the form: m
~(x,h)=f(x)-
(15)
where dimi”. is the minimum section i, and
diameter
of pipe
xi may.= (l/dimin.)4+14 i= l-_-n.
(l2a)
These constraints are formed by the requirement that the diameter of all pipe sections on any path j must be greater than or equal to that of the relief valve j. Hence, each dimin’,i = 1 ---n, is found by considering all paths j on which section i lies, and taking the largest of all the corresponding valve diameters. The presence of rounding errors, and of course computer storage limitations, limit the size of linearly constrained nonlinear programming problems which can be solved by gradientprojection methods to approximately 100 variables, in this case 100 pipe sections. Since the number of pipe sections, IZ, is always very much greater than the number of paths, m, it is attractive to consider the possibility of reducing
2 Ai
(16)
i=l
Kuhn and Tucker[7] have presented a set of necessary and sufficient conditions for the existence of a solution to (7) and (8). In terms of the Lagrangian, (16), they may be stated as fol1ows: If X solves (7) and (8), then, (assuming f(X) is differentiable), there exists ax such that: 1. s(X,h)=Oi.e.F=
z
2m ahii 1 i=
l---n. (17)
This means that the gradient of the objective function must be orthogonal to the surface formed by the boundaries of the active linear so no further progress can be constraints, made. 2. g
1134
(Z,i)
2 Oi.e.iatixr-bJ
s 0 j=
l---m.
i=l
(18)
An approach to the optimal design of networks
This merely states that the point X must be feasible, i.e. satisfy the constraints given by (8). 3. XS Oi.e.X, 5 Oj= l---m. (19) This, in conjunction with (17), means that the negative gradient (steepest descent) direction must be pointing orthogonally outwards from the feasible region. (If it were pointing inwards with respect to any given constraint, further progress could obviously be made by moving off that constraint boundary.)
min.
max.
x E En
A 5 0.
b(x, A) (22)
Consider the function defined as follows: Y(A) = min
4(x, A)
x E En, A fixed.
(23)
We shall call this the dual function. Provided that this function exists for all A % 0 and that 4(x, A) attains a finite minimum, the solution of the saddle-point problem is obtained by solving the dual problem:
4. XT !$ (jE,h) = 0 max. y(A) i.e. either& = 0 or 2 a. R - bi = 0 <=I *( j=
l-_-m.
AS0
(20)
This means that if X is not on the boundary of any constraint, j; then its corresponding Lagrange multiplier, X, , is zero. These four conditions are necessary, they are also sufficient if the optimization problem is convex.* The inequalities (18) and (19) reflect the fact that the constraints (8) are inequality constraints; for equality constraints (18) would be an equality and the corresponding Lagrange multipliers would be unrestricted in sign. By considering the left hand relations in (17)(20), it is evident that, under convexity assumptions, the Kuhn-Tucker conditions are also necessary and sufficient for the existence of a saddle-point of 4 (x, A) at (X, A) , i.e.
($(%A) 5 (%,A) 5 4(X,X)
for all ASO,XEE~.
or, alternatively,
(21)
Thus the solution of (7) and (8) may be obtained by solving the saddle-point problem: *A convex programming problem is defined as a problem where the constraints define a feasible region which is convex, i.e. the line segment between every two points in the feasible region is also in the feasible region, and the objective function being minimized (maximized) is convex (concave).
min. - y(A) A 5 0.
(24)
It is, in fact, possible to generalize the definition of the dual function to the form: y(A) = min $(x, A) x E R, A fixed.
(23a)
where R is some subset of En. This amounts to a partitioning of the original set of constraints into two sets: one set is adjoined to the objective function with Lagrange multipliers to form the Lagrangian function, and the other set remains to provide constraints on the associated minimization problem implied in Eq. (23). The dual problem, (24) looks deceptively simple; it is an m-dimensional problem and the only constraints are simple upper bounds. However, note from (23) that each evaluation of the dual function involves a complete minimization. Unless the structure of the problem makes this minimization relatively simple, there is little advantage to be gained by considering the dual problem in preference to the primal problem. There are other difficulties in being assured that the dual problem is a convex programming problem, even though the primal problem may be convex[8]. A sufficient condition which ensures that the dual problem is convex is that the primal constraints be convex, and the primal objective function be strictly convex.
1135
B. A. MURTAGH
It will be evident that there are some analogies between duality theory and Pontryagin’s Maximum Principle and its extensions, on which there has been much written in recent Chemical Engineering literature [9- 1 I]. In particular, the Maximum Principle is related to the right-hand inequality in (21). (In this case it is a minimum principle since we are concerned with minimizing, not maximizing, the Lagrangian function.) For an excellent discussion of these analogies, the reader is referred to Mangasarian’s recent book [ 121. Duality theory also forms the basis of various decomposition techniques in linear and nonlinear programming; an application to recycle processes has been proposed by Brosilow and Lasdon[ 131. Much of this related work is concerned with equality constraints, so the problem becomes one of determining a stationary saddle point of the associated Lagrangian function. For nonlinear equality constraints, it is not a trivial matter to determine that a stationary point satisfies Eq. (2 1) and, hence, solves the primal problem. 5. DUAL
FORMULATION
Consider the structure of the network optimization problem formulated in Section 3. The total cost of the system is the sum of the costs of each pipe section, so that regardless of how this cost is expressed, the objective function is always separable into a sum of functions of the individual variables xi, i = 1 ---n. i.e.
f(x) = 2 j&q). i=l
(25)
Fig. 2. Form of objective function for an individual pipe section.
, on each variable would alone prevent this aiymptotic value being approached. In the region 0 < xi s ximas., i= 1 ---n, the second derivative of each h(xi) is greater than zero, which is a sufficient condition for J(xi) to be strictly con-
X,max.
vex, and hencef(x) zigIh(Xi) to be strictly convex. The linear constraints represented by Eq. (14) do not form a closed feasible region; however, the coefficients ati and bj are all positive and so provide an upper limit on each Xi, which corresponds to the direction of decreasing function value. Considering only the linear constraints given by (14), the Lagrangian function corresponding to (13) and ( 14) takes the form:
Note, also, that by expressing the cost as the weight of each pipe section, the form of each individual function is the same for each variable: fi(xJ = Li(~+~xi-o.20’7)
i= 1 ---n.
(26)
The graph of this function is shown in Fig. 2. Observe that in the region of interest, Xi > 0, fi(xJ is monotonically decreasing and bounded below by an asymptotic value of Lia. Regardless of the other linear constraints, the upper bound, 1136
=
gl [4(Di +flXi-"'2077)
An approach to the optimal design of networks
minimize - y (A)
This, in turn, is strictly convex in x in the range 0 < x S xmax.for all A I 0. Define the dual functions as follows: y(A) = minimum 4(x, A) over x E R = (xl0 < xi 5 ximax., i= l---n} A fixed.
subjectto&
(28)
The solution of the minimization problem implied in Eq. (28) is obtained analytically. For a stationary point, we have: v
=
0 at xi
=
t
A&L&
(-5
0.2077)
s 0 j=
(34) I---m.
(35)
The dimensionality of the problem to be solved has been reduced from n to m, and the task of calculating the n optimal pipe diameters, dii = 1 ---n, has been reduced to a trivial sub-problem. A flow diagram of the procedure for solving the dual problem is shown in Fig. 3. The steps which are taken are as follows~ (1) Assume
an initial set of values of Ajj =
1- - - ~ltall less than zero.
-l/1.2077 Enter
5=1
(29) Set
= xi*(A), say.
(30)
initial
“OIWS
for A
8
This must be a minimum if xc(h) s Ximax., because of the strict convexity of the Lagrangian. Note that this is finite unless hj(S 0) is zero for all paths j such that i E Sj. The solution point of the subproblem in (28) is:
Calculate x*(x, from
Eqn.(30) t
Calculate min x”(X),
B(h)= 1
I
%(A) =min.Of
($(A)pXi”““.),i=
l---n.
(31)
Calculate
r(X)
from Eqn.i32)
t
This exists, and is unique and finite, for all hj S 0. The dual function for our particular problem is thus:
Use gradient pmjectbn method to ge value
‘Fnw L -_-I
Test KuhnTucker condns for optimality
Y(A) = 46(A), A) = gI &(a + &%(A)-0’2077) - fI Aj ( :I atixi(A) - S) (32) with i(A) given by (3 1). The gradient of the dual function, which exists and is continuous for all AI S 0,j = 1 ---m, is:
ayw_ -_aAj
2 au&(A)-bj)
j=
l---m.
Calculate dioms. from Eqn (I 2). Recolculote C, i= Ibn frc.n Eqnwmd aquiv.lgth.CweCtsj I
(33)
We are therefore able to use a gradient projection method of minimization to solve the linearly constrained dual problem: 1137
Fig. 3. Flowchart
of calculation
procedure.
B. A. MURTAGH
(2) Calculate the optimal values of xii= 1 ---n, for the given A from Eqs. (30) and (3 1). (3) Calculate the value of the dual objective function from Eq. (32), and its gradient from Eq. (33). (Note that we wish to maximize the dual, so the objective function in a minimization program is the negative of the dual.) (4) Test for optimality. If the assumed value of A is optimal, stop. Otherwise, use the gradientprojection minimization program to assign a new value of A 5 0. Return to step 2. Note that the strict convexity in x of the Lagrangian function 4(x, A) in the range 0 < xi 5 xtmaK.for all values of A 5 0 ensures the uniqueness of x(A). Mere convexity of the Lagrangian may, in general, give rise to a non-unique point x(A) which solves the dual sub-problem, and it is evident that the complication of non-uniqueness is transmitted to the evaluation of the gradient of the dual function. For computational purposes, we may avoid difficulty with the strict inequality in the definition of the region R by placing a lower bound E on Xi, i = l---n, where E is a small positive constant chosen to ensure that no difficulty arises through rounding errors. In the region R, sufficient conditions for a convex dual problem are satisfied so that any optimum found in step 4 of the above procedure will be global over R. The constraints which were used to define R still enable i(A) to be written down term by term. It is possible to incorporate any form of lower and upper bound in each variable without losing this ability, but this ability is immediately lost if a constraint used to define R contains terms in more than one variable. Such constraints must be adjoined to the Lagrangian, each requiring an additional Lagrange multiplier, and it is evident that any advantage in dimensionality of the dual problem is lost if there are a large number of these. Note that the upper bound, ximax., keeps each ii (A) finite and thus the corresponding diameter greater than zero. Note, also from Fig. 1, that for each path j, there is at least one pipe section i (the section adjacent to the valve) which is on that path
only, i.e. there is at least one i for which
j=l
has only the one non-zero ati term. If the bounds x.~~~., i= l---n were absent, then at the optimum all hj, j = 1 ---m, would be non-zero, which means that all the constraints of the primal problem given by Eq. (14) would be active. The physical interpretation of this is that the pipe diameters would be squeezed down to utilize the maximum pressure drop available in each in ridiculously small path, often resulting diameters in the sections adjacent to the valve. We have thus far considered the kinetic energy correction factor, Cii = 1 ---n, to be constant. While this is true for incompressible fluid flow, it will be seen from Eq. (2) that for compressible fluid flow each C.depends in a nonlinear manner on both the diameter and terminal pressures of each section. It is therefore necessary in this particular application to construct an outer iterative cycle (see Fig. 3) in which the values of Ci are recalculated using the optimal diameters and the resulting pressure drops. As the reviewer has pointed out, it is not possible to prove that such a cycle will converge to the correct solution, or that it will converge at all. This difficulty, which is present in both the primal and dual formulations, is not just a technical matter of being unable to provide a convergence proof; the nonlinearity of the correction factor Ci makes even observed convergence of this outer iterative cycle subject to doubt as to the validity of the converged solution, and justification for accepting it must be supplied from other considerations. Along this line, it is possible to obtain an accurate measure of the actual discrepancy in the constraints at each iteration; given the diameters of each section and starting with the known flare pressure, P,, it is possible to show that the sequential calculation of pressure by the iteration P.r+l 1 =f(P,‘) for each section, where
1138
(36)
An approach to the optimal design of networks
X [
Pi’ 1 + Kidi”“14 log10 PI_, -
(
l/2 )I)
(37)
and where Ki = 4’61/4fiLidi”‘*14 will converge for all values of Pi-1 and di in the range of interest. Thus it is possible to calculate the actual back pressure at each valve, say Pj, and determine the discrepancy, pj-Pj spec., in the total pressure drop for those constraints j which are active. It is also possible to obtain an estimate of the penalty incurred in the objective function for each of these. Let & = I’j’- PO’, and Ab = &jbj. A well-known economic interpretation of Lagrange multipliers (see for example [2]), which applies equally to linear and nonlinear inequality constraints, is that hj is a measure of the extent to which the optimum value of the objective function, say f, could be further decreased by increasing the bound bj on the constraint, i.e. hj = (a! (E)/abj). We may thus get an approximate estimate of the effect of each discrepancy on the objective function by the term XjAbj. (This is only an approximate estimate, since we are assuming that x approximates the solution to the problem with the bound 6, on each constraint j, and that Abj is sufficiently small for a first order estimate to be adequate.) Since in practice the correction caused by Ci is itself small[5], it may be satisfactory merely to know the actual discrepancy in pressure drop. 6. COMPUTATIONAL
EXPERIENCE
Although emphasis has been placed on the approach rather than the application, this paper would not be complete without some description of computational experience. The optimization method used to solve the problems is a recently developed ‘variable-metric’ projection method, so called because it accumulates and utilizes an approximation to the local inverse matrix of second partial-derivatives of the objective function. This method has been used very success-
fully in unconstrained minimization ([ 14, 151); the basis of the approach in the presence of linear constraints is the use of projection operations which confine the search to the subspace formed by the active constraints while still retaining complete curvature information. A full description and theoretical justification of the method may be found in Refs. [6] and [ 161. Two practical points should be mentioned at this stage. The first is that, because of the various classifications of pipes which are made in practice, the number of pipe sections n is much greater in comparison to the number of paths m than that implied by the definition in Section 2. Also, it is possible to include an ‘equivalent length’ correction for elbows and bends by including this in the term Ci; the same remarks made in the previous section concerning convergence of the iterative adjustment of Ci apply equally here. The primal problem formulated in Section 3 is quite tractable; the main motivation for continuing on to formulate the dual problem was the fact that relief header networks of well over a hundred pipe sections are becoming more common as refineries increase in size and complexity. In order to make a comparative evaluation of the primal and dual formulations a number of medium sized problems taken from previous contracts were solved using both. The dual problem was consistently the most easily solved, both in terms of number of function calculations required and total execution time. For example, in one problem consisting of 3 1 pipe sections and 9 valves, the primal formulation required 35 function calculations and took 8.15 set computing time on a CDC 6400, whereas the dual formulation required 10 function calculations and took 1.31 sec. The maximum discrepancy in pressure drop resulting from the assumed and calculated values of Ci (including an equivalent length term) was 1 per cent after two iterations for this problem. The solution obtained is subjected to a postoptimal treatment, which selects the appropriate commercial pipe diameters by a partially enumerative search and also checks for critical
1139
B. A. MURTAGH
velocity in each pipe section. The program is now in full use as a design tool and has already achieved considerable savings in capital cost on a number of contracts. The largest system which has been solved consisted of a total of 225 pipe sections. 7. DISCUSSION AND SOME GENERALIZATIONS
The two reasons why the dual formulation turns out to be so successful computationally are the fact that the constrained sub-problem is solved analytically, and that the number of nontrivial constraints in the primal problem is much less than the number of variables. Thus we not only reduce the number of variables by formulating the dual problem, but also are able to avoid the -usual handicap of having a “minimization within a minimization.” The approach presented in this paper is not restricted to fluid-flow systems, but can be applied to any system whose configuration can be represented by a network. By observing the analogies in other contexts of the three basic variables: pressure drop or change in potential, flowrate, and diameter or resistance to flow, similar optimization problems may be formulated in such diverse areas as transportation, scheduleconomics, and of course engineering im systems. The approach can be applied in a straightforward manner to equations of the form: (Potential Drop)‘] = (Flow)j2 X Resistance)j3. Any departure from this form by the use of correction factors, as in the example presented here of compressible fluid flow, would require further treatment and physical reassurance that the contribution was small. Note also that the equations which describe the system being optimized are assumed to be nonlinear; in the case of linear systems the use of linear programming is already well-established. One may often be, as in this example, dependent on the nonlinearity of the flow equations with respect to the independent variables in order to obtain a strictly convex primal objective function. Although the network considered here consisted of multiple sources but only one sink,
there is no more difficulty in considering multiple sources and sinks as long as each path between source and sink is not ambiguously defined and the flowrate in each section is either specified or calculable. The number of sections will always be greater than the number of paths, regardless of the configuration of the network, so it would seem that the dual formulation would always be an attractive alternative in such applications. Acknowledgements-The author wishes to thank Fluor Corporation for permission to publish this paper. Although mentioned only briefly in the text, a considerable effort was expended by Manuel de la Puente and Chan-shan Fan in developing a program for selecting the appropriate commercial pipe sizes and checking the resultant pressure drops rigorously. The author also wishes to thank the reviewer for helpful comments and constructive criticism.
NOTATION
a,
coefficients of constraint equations (see Eq. 8) bj constant term of constraint equations (see Eq. 8) kinetic energy correction factor (see ci Eqs. 2 and 4) di diameter of pipe section i, in. dimin. minimum diameter of pipe section i, in. En euclidean n-space i subscript to denote pipe section j subscript to denote path from source to sink (and, in the problem presented here, the corresponding valve) & constant term in pressure-drop equation for pipe section i (see Eq. 2) ~~ length of pipe section i, ft p,, flare pressure, psia exit pressure of section i, psia pi P5spec. specified maximum pressure at valve j, psia R region in En (see Eq. 28) sj the set of pipe section numbers in path j weight per unit length of pipe section i, wi lb/ft independent variable in primal minimXi ization problem X;(A) unconstrained optimal value of x,, given X (see Eqs. 29 and 30)
1140
An approach to the optimal design of networks X,max.
I
maximum
value
of
xi :ximaX.=
~/(~imin.)4+314
ii(X) f
optimal value of Eq. 3 1) optimal value of x
Xi
for a given A (see
Greek symbols
(Y,p
coefficients
(see Eq. 1) dual function (see Eq. 23) 4 Lagrangian function (see Eq. 16) Xj Langrange multiplier associated with constraint j of primal minimization problem. Independent variable in dual problem. X optimal value of A y
relating weight of pipe to dia. REFERENCES
VI FULKERSON
PI HU T. C.,
18. and Network Flows. Addison-Wesley,
D. R., .I. Sec. Ind.
Integer Programming
Appl. Math. 19619
Mass. 1969. P., Notes to accompany lectures, NATO Summer School on Integer and Nonlinear Programming, Bandol, France 1969. API-RP520, Recommended Practice for the Design and Installation of Pressure-Relieving Systems in Refineries: Part I-Design. p. 22, Third Edn., American Petroleum Institute 1967. PERRY J. H., ChemicalEngineer’s Handbook, 4th Edn., pp. 5-24. McGraw-Hill, New York 1963. MURTAGH B. A. and SARGENT R. W. H., Optimization, pp. 215-246, (Ed. R. FLETCHER). Academic Press, London 1969. KUHN H. W. and TUCKER A. W., Proc. Second Berkeley Symposium on Mathematical Statistics and Probability, pp. 481-492, (Ed. J. NEYMAN), University of California Press, Berkeley 1951. FALK J. E.,/. Math. Anal. Appt. 1967 19 141. HORN F. and JACKSON R., Ind. Engng Chem. Fund/s. 1965 4 110; 1965 4 239; 196.5 4 387. JACKSON R., Chem. Engng Sci. 1964 19 19; 1964 19 253; 1965 20 405. KATZ S.,A.I. Ch.E.-Inst. Chem. Engrs. Symp. Ser. 1965 4 69. MANGASARIAN 0. L. Nonlinear Programming. McGraw-Hill, New York 1969. BROSILOW C. and LASDON L.,A.Z. Ch.E.-Inst. Chem. Engrs. Symp. Ser. 1965 4 75. FLETCHER R. and POWELL M. J. D., C0mpt.J. 1963 6 163. MURTAGH B. A. and SARGENT R. W. H., C0mpt.J. 1970 13 185. SARGENT R. W. H. and MURTAGH B. A., paper presented at the 7th International Symposium on Mathematical Programming, The Hague, Sept. 1970.
131 WOLFE [41 [51 [61 [71 [81 [91 [lOI u 11 t121 [I31 [I41 t151 [I61
R&urn&-
Le probleme de conception trait& dans cette etude se rapporte 21la distribution optimale de la chute de pression dans un reseau oh le debit est d&ermine mais oti les diametres des tubes peuvent varier. Si l’on utilise des elements de la theorie de dualite de programmation non-lineaire, le probleme semble pouvoir etre resolu a I’aide dun calculateur. Les auteurs presentent une approche qui a une applicabilite generale aux reseaux dans lesquels les equations de I’ecoulement sont non-lineaires. Zusammenfassung-
Das in diesem Artikel behandelte Konstruktionsproblem betrifft die optimale Verteilung des Druckabfalls iiber ein Netzwerk, in welchem die Stromungsgeschwindigkeiten festgelegt sind, w&end die Rohrdurchmesser variieren kiinnen. Unter Verwendung der Dualitatstheorie nicht-linearer Programmierung kann das Problem in eine computermassig zugangliche Form gebracht werden. Die dargelegte Methode ist von allgeneiner Anwendbarkeit auf Netzwerke, in welchen die Stromungsgleichungen nicht-linearer Art sind.
1141
Pergamon
Press.
Printed
in Great
Britain
An approach to the optimal design of networks B. A. MURTAGH Computer Sciences Department, Fluor Corporation, Los Angeles, U.S.A. (First received 22 November 1970; accepted 9 September 197 1) Abstract-The design problem considered in this paper is the optimal distribution of pressure drop over a network in which flow-rates are specified, but pipe diameters may vary. By using elements of the duality theory of nonlinear programming, the problem is cast in a form which is computationally attractive. The approach presented is of general applicability to networks in which the equations of flow are nonlinear. 1. INTRODUCTION
has been much work recently on the application of mathematical programming to problems which arise in determining the optimal distribution of flows in networks[l, 21. The algorithms obtain optimal flow patterns by minimizing some cost function, while satisfying the laws of conservation around the junctions of the network. The problem treated in this paper is not the usual type of network flow problem in that the flows throughout the network are considered to be fixed. The maximum total drop in potential energy in each path from source to sink is also specified. The problem is to adjust the resistance in each arc of a path to achieve a minimum total cost of the network, considered to be a function of resistance. The treatment which is given is directed to a specific problem which frequently arises in chemical engineering design: one of compressible fluid flow through a network of pipes. However, there are many analogous problems which arise in other contexts - some generalizations will be discussed in Section 7. The optimization problem which arises enjoys a particularly simple structure, which is exploited by formulating a second (“dual”) problem, the solution of which also solves the original (“primal”) problem. The concept of duality in nonlinear programming has often been passed off as a mathematical nicety of little practical value. (For example, Wolfe[3] remarked in a recent review that, “I have yet THERE
to see a case in which the dual problem is truly easier to solve than the primal”.) However, in this particular instance, the dual problem is not only theoretically pleasing, but also much more easily solved, as computer results show. 2. A CHEMICAL RELIEF
ENGINEERING PROBLEM: HEADER DESIGN
Consider a network of pipes as shown in Fig. 1. These pipes lead from relief valves on various refinery vessels to a discharge zone, or flare. A pipe section is defined as the length of pipe between two junctions or between a valve or the flare and a junction. A path is defined as the sequence of pipe sections between any one valve and the flare. The configuration of the pipe network is assumed to be already determined by the layout of the refinery vessels. Thus the lengths of all the pipe sections are specified and the only design variables remaining are their diameters. In an emergency, the vessels must discharge their contents through the relief valve at maximum (sonic) velocity for an initial transitory period. The problem is to select the pipe diameters in each path from valve to flare so that the total pressure drop over the path remains within a pre-specified limit. Assuming that the flare pressure is specified (usually slightly above atmospheric pressure), the limit is imposed so that the back-pressure in the pipe at the valve is below the critical value, thus ensuring sonic velocity in the throat of the valve.
1131
B. A. MURTAGH A
t
wheref=
Flare
friction factor = 0.0475 (Rc+-O.‘~~
+---I I5
13
=
o.a475
(This gives factor.)
14
(3)
(!E3)-0’186
a conservative
value
of friction
and
Valve
I
2
3
4
Fig. 1. Schematic representation
5
6
7
6
of a piping network.
There are obviously many feasible solutions to this problem. Note also that some pipe sections are common to many paths, so the problem cannot be solved for each path individually. In formulating this as an optimization problem, we impose all the pressure drop limitations as a set of constraints, and choose as the ‘best’ feasible solution that set of pipe diameters which minimizes the total weight of pipe. The weight per unit length of pipe is assumed to be a linear function of diameter; (this is a valid assumption when the diameter is much larger than the thickness): Wi= (Y+ pdi lb/ft.
(1)
where di is the diameter in inches and the subscript i refers to the number of the pipe section. The coefficients a and fl are obtained by linear regression. Following API-RP520 design standards [4], the appropriate equation is for isothermal flow of a compressible fluid. This equation takes the form [5] : Pz_P2=64.f.L.W2.R.T. A B r2.gc.M.D5
x log10EpA
(2)
PA = Ps = D = L= W= R = T= M = gc = p =
inlet pressure (abs) (lbF/ft2) outlet pressure (abs) (lb,/ft2) dia. (ft) length of pipe (ft) weight rate of flow (lblsec) gas constant = 1546 (ft-1bJlb mole “R) absolute temperature (“R) molecular weight (lb/lb mole) 32-17 (lb ft/lbF sec2) viscosity (centipoise).
For any section of pipe, the first three terms above vary, while those below are constant. The properties at the valve are specified, and the properties downstream are calculated by simple mixing rules. (There are many assumptions worthy of contention used here, especially that of isothermal flow. However, we will not take issue on these since they are either set forth by API-RP520, the usual standard in contracts, or used in practice with safety as the justification.) The pressure drop over each section of pipe is given by an equation of the form of (2). Representing all constant terms by one constant, Ki, for each pipe section, i, and leaving the presence of the subscript as understood, Eq. (2) with (3) becomes: P,z_P:_I=gg where C, represents the term in brackets in Eq. (2), the kinetic energy correction factor, and Pi and Pi-1 are the inlet and outlet pressures respectively. Note that Ci is a function of Pi and Pi_l; we will, however, ignore this fact
1132
An approach to the optimal design of networks
for the present and treat it as a constant. (This point will be discussed later.) Consider r pipes in series numbered from 1 to r consecutively; the individual pressure drops are calculated by recursion:
(5) P,“-P,2_,=f-g$.
r
Hence, summing:
boundaries on the feasible region in which the point x may lie. A constraint j is deemed to be active if the current point satisfies (8) exactly, so that the equality sign holds. By the introduction of some convenient notation, it is possible to express the problem described in the previous section in the form of (7) and (8). Let n be the number of pipe sections and m be the number of paths. Assign a unique number between 1 and n to each pipe section, and a unique number between 1 and m to each valve (see Fig. 1). It is evident from Fig. 1 that this last set of numbers may be used to refer to both the valve itself and the path from the valve to the flare since there is a one-to-one correspondence. Let the subscript i represent the number of the pipe section, and the subscript j represent the number of the path. Let Sj be the set of pipe section numbers in path j. Then, for example, in Fig. 1:
where P,. and P, are the inlet and exit pressures respectively. The problem is to adjust each dia., diy to minimize the sum of LiWi over all pipe sections, i, while obeying pressure drop limitations of the form given by Eq. (6).
S, = (1,9,13,15) S, = (2,9,13,15)
SS = (8, 12, 14, 15). 3. OPTIMIZATION:
PRIMAL
FORMULATION
The objective function may be expressed simply as the sum of the weights of each pipe section: From Eq. (1):
Linearly constrained nonlinear optimization problems may be expressed in the form: minimizef(x)
= _&x1---x,)
subject to 5 a x. d bi i=l 15t
j=
(7)
l-_-m
.f(dl ---dn) = (8)
where n is the number of variables, and m is the number of constraints. Ax) represents the objective function, which usually expresses some form of cost associated with the values of the independent variables, xt, i = 1 ---n. ati is the coefficient of variable i in constraint j, and b, is the corresponding constant term. Geometrically, the equations
represent
$, au.Xi-bj=O
j=
n-dimensional
hyperplanes
l---m
(8a)
which form
5 &(a
i=l
+ pdi).
(9)
The constraints may also be written down easily with the aid of the above notation. There will be m pressure drop equations of the form of Eq. (6), one for each path j, j = 1 ---m. They may be expressed as inequality constraints by using the maximum allowable inlet pressure in place of the actual inlet pressure: Let bJ = Pj” spec. - PO2
j=
l---m
(10)
where Pj spec. is the specified maximum allowable back pressure at valve j, and P,, is the flare pressure.
1133
B. A. MURTAGH
Let ati =
I
&C,, 0
i E si ,
i= l---n,
l---m
i$Ssj.
(11)
In order to make the constraints introduce the transformation: xi=
j=
l/di4’R14 i= l-_-n.
(12)
minimize f(x) = E L*((Yf /3x,-0.20”)
(13)
i=l
$l a+~~ 5 bj
j=
4. REVIEW
linear, we
We then have the problem:
subject to
the dimensionality from n to m.
1 ---m
(14)
which is the desired form. Gradient projection methods of nonlinear programming [6] are particularly well suited to problems posed in this form. There are also auxiliary constraints in the form of lower bounds on the pipe diameters: di 2 dimin. > 0 + xi 5 ximal. < ~4 i= 1 ---II
of the optimization
OF DUALITY
problem
THEORY
In this section we shall briefly describe those essential elements of the duality theory of nonlinear programming which are required. Duality theory has been put forward in various guises, but the most easily visualized concept is that of a saddle point of an associated Lagrangian function. The prima1 problem involves minimizing in one space of variables, while the dual problem involves maximizing in another space of variables. Under certain regularity assumptions, the two solutions are equal and correspond to a saddle point (not necessarily a stationary point) of the Lagrangian in the combined space of both variables. Consider the problem expressed by (7) and (8). Using Lagrange multipliers, hj, to adjoin the constraint equations to the objective function we obtain a Lagrangian function, which takes the form: m
~(x,h)=f(x)-
(15)
where dimi”. is the minimum section i, and
diameter
of pipe
xi may.= (l/dimin.)4+14 i= l-_-n.
(l2a)
These constraints are formed by the requirement that the diameter of all pipe sections on any path j must be greater than or equal to that of the relief valve j. Hence, each dimin’,i = 1 ---n, is found by considering all paths j on which section i lies, and taking the largest of all the corresponding valve diameters. The presence of rounding errors, and of course computer storage limitations, limit the size of linearly constrained nonlinear programming problems which can be solved by gradientprojection methods to approximately 100 variables, in this case 100 pipe sections. Since the number of pipe sections, IZ, is always very much greater than the number of paths, m, it is attractive to consider the possibility of reducing
2 Ai
(16)
i=l
Kuhn and Tucker[7] have presented a set of necessary and sufficient conditions for the existence of a solution to (7) and (8). In terms of the Lagrangian, (16), they may be stated as fol1ows: If X solves (7) and (8), then, (assuming f(X) is differentiable), there exists ax such that: 1. s(X,h)=Oi.e.F=
z
2m ahii 1 i=
l---n. (17)
This means that the gradient of the objective function must be orthogonal to the surface formed by the boundaries of the active linear so no further progress can be constraints, made. 2. g
1134
(Z,i)
2 Oi.e.iatixr-bJ
s 0 j=
l---m.
i=l
(18)
An approach to the optimal design of networks
This merely states that the point X must be feasible, i.e. satisfy the constraints given by (8). 3. XS Oi.e.X, 5 Oj= l---m. (19) This, in conjunction with (17), means that the negative gradient (steepest descent) direction must be pointing orthogonally outwards from the feasible region. (If it were pointing inwards with respect to any given constraint, further progress could obviously be made by moving off that constraint boundary.)
min.
max.
x E En
A 5 0.
b(x, A) (22)
Consider the function defined as follows: Y(A) = min
4(x, A)
x E En, A fixed.
(23)
We shall call this the dual function. Provided that this function exists for all A % 0 and that 4(x, A) attains a finite minimum, the solution of the saddle-point problem is obtained by solving the dual problem:
4. XT !$ (jE,h) = 0 max. y(A) i.e. either& = 0 or 2 a. R - bi = 0 <=I *( j=
l-_-m.
AS0
(20)
This means that if X is not on the boundary of any constraint, j; then its corresponding Lagrange multiplier, X, , is zero. These four conditions are necessary, they are also sufficient if the optimization problem is convex.* The inequalities (18) and (19) reflect the fact that the constraints (8) are inequality constraints; for equality constraints (18) would be an equality and the corresponding Lagrange multipliers would be unrestricted in sign. By considering the left hand relations in (17)(20), it is evident that, under convexity assumptions, the Kuhn-Tucker conditions are also necessary and sufficient for the existence of a saddle-point of 4 (x, A) at (X, A) , i.e.
($(%A) 5 (%,A) 5 4(X,X)
for all ASO,XEE~.
or, alternatively,
(21)
Thus the solution of (7) and (8) may be obtained by solving the saddle-point problem: *A convex programming problem is defined as a problem where the constraints define a feasible region which is convex, i.e. the line segment between every two points in the feasible region is also in the feasible region, and the objective function being minimized (maximized) is convex (concave).
min. - y(A) A 5 0.
(24)
It is, in fact, possible to generalize the definition of the dual function to the form: y(A) = min $(x, A) x E R, A fixed.
(23a)
where R is some subset of En. This amounts to a partitioning of the original set of constraints into two sets: one set is adjoined to the objective function with Lagrange multipliers to form the Lagrangian function, and the other set remains to provide constraints on the associated minimization problem implied in Eq. (23). The dual problem, (24) looks deceptively simple; it is an m-dimensional problem and the only constraints are simple upper bounds. However, note from (23) that each evaluation of the dual function involves a complete minimization. Unless the structure of the problem makes this minimization relatively simple, there is little advantage to be gained by considering the dual problem in preference to the primal problem. There are other difficulties in being assured that the dual problem is a convex programming problem, even though the primal problem may be convex[8]. A sufficient condition which ensures that the dual problem is convex is that the primal constraints be convex, and the primal objective function be strictly convex.
1135
B. A. MURTAGH
It will be evident that there are some analogies between duality theory and Pontryagin’s Maximum Principle and its extensions, on which there has been much written in recent Chemical Engineering literature [9- 1 I]. In particular, the Maximum Principle is related to the right-hand inequality in (21). (In this case it is a minimum principle since we are concerned with minimizing, not maximizing, the Lagrangian function.) For an excellent discussion of these analogies, the reader is referred to Mangasarian’s recent book [ 121. Duality theory also forms the basis of various decomposition techniques in linear and nonlinear programming; an application to recycle processes has been proposed by Brosilow and Lasdon[ 131. Much of this related work is concerned with equality constraints, so the problem becomes one of determining a stationary saddle point of the associated Lagrangian function. For nonlinear equality constraints, it is not a trivial matter to determine that a stationary point satisfies Eq. (2 1) and, hence, solves the primal problem. 5. DUAL
FORMULATION
Consider the structure of the network optimization problem formulated in Section 3. The total cost of the system is the sum of the costs of each pipe section, so that regardless of how this cost is expressed, the objective function is always separable into a sum of functions of the individual variables xi, i = 1 ---n. i.e.
f(x) = 2 j&q). i=l
(25)
Fig. 2. Form of objective function for an individual pipe section.
, on each variable would alone prevent this aiymptotic value being approached. In the region 0 < xi s ximas., i= 1 ---n, the second derivative of each h(xi) is greater than zero, which is a sufficient condition for J(xi) to be strictly con-
X,max.
vex, and hencef(x) zigIh(Xi) to be strictly convex. The linear constraints represented by Eq. (14) do not form a closed feasible region; however, the coefficients ati and bj are all positive and so provide an upper limit on each Xi, which corresponds to the direction of decreasing function value. Considering only the linear constraints given by (14), the Lagrangian function corresponding to (13) and ( 14) takes the form:
Note, also, that by expressing the cost as the weight of each pipe section, the form of each individual function is the same for each variable: fi(xJ = Li(~+~xi-o.20’7)
i= 1 ---n.
(26)
The graph of this function is shown in Fig. 2. Observe that in the region of interest, Xi > 0, fi(xJ is monotonically decreasing and bounded below by an asymptotic value of Lia. Regardless of the other linear constraints, the upper bound, 1136
=
gl [4(Di +flXi-"'2077)
An approach to the optimal design of networks
minimize - y (A)
This, in turn, is strictly convex in x in the range 0 < x S xmax.for all A I 0. Define the dual functions as follows: y(A) = minimum 4(x, A) over x E R = (xl0 < xi 5 ximax., i= l---n} A fixed.
subjectto&
(28)
The solution of the minimization problem implied in Eq. (28) is obtained analytically. For a stationary point, we have: v
=
0 at xi
=
t
A&L&
(-5
0.2077)
s 0 j=
(34) I---m.
(35)
The dimensionality of the problem to be solved has been reduced from n to m, and the task of calculating the n optimal pipe diameters, dii = 1 ---n, has been reduced to a trivial sub-problem. A flow diagram of the procedure for solving the dual problem is shown in Fig. 3. The steps which are taken are as follows~ (1) Assume
an initial set of values of Ajj =
1- - - ~ltall less than zero.
-l/1.2077 Enter
5=1
(29) Set
= xi*(A), say.
(30)
initial
“OIWS
for A
8
This must be a minimum if xc(h) s Ximax., because of the strict convexity of the Lagrangian. Note that this is finite unless hj(S 0) is zero for all paths j such that i E Sj. The solution point of the subproblem in (28) is:
Calculate x*(x, from
Eqn.(30) t
Calculate min x”(X),
B(h)= 1
I
%(A) =min.Of
($(A)pXi”““.),i=
l---n.
(31)
Calculate
r(X)
from Eqn.i32)
t
This exists, and is unique and finite, for all hj S 0. The dual function for our particular problem is thus:
Use gradient pmjectbn method to ge value
‘Fnw L -_-I
Test KuhnTucker condns for optimality
Y(A) = 46(A), A) = gI &(a + &%(A)-0’2077) - fI Aj ( :I atixi(A) - S) (32) with i(A) given by (3 1). The gradient of the dual function, which exists and is continuous for all AI S 0,j = 1 ---m, is:
ayw_ -_aAj
2 au&(A)-bj)
j=
l---m.
Calculate dioms. from Eqn (I 2). Recolculote C, i= Ibn frc.n Eqnwmd aquiv.lgth.CweCtsj I
(33)
We are therefore able to use a gradient projection method of minimization to solve the linearly constrained dual problem: 1137
Fig. 3. Flowchart
of calculation
procedure.
B. A. MURTAGH
(2) Calculate the optimal values of xii= 1 ---n, for the given A from Eqs. (30) and (3 1). (3) Calculate the value of the dual objective function from Eq. (32), and its gradient from Eq. (33). (Note that we wish to maximize the dual, so the objective function in a minimization program is the negative of the dual.) (4) Test for optimality. If the assumed value of A is optimal, stop. Otherwise, use the gradientprojection minimization program to assign a new value of A 5 0. Return to step 2. Note that the strict convexity in x of the Lagrangian function 4(x, A) in the range 0 < xi 5 xtmaK.for all values of A 5 0 ensures the uniqueness of x(A). Mere convexity of the Lagrangian may, in general, give rise to a non-unique point x(A) which solves the dual sub-problem, and it is evident that the complication of non-uniqueness is transmitted to the evaluation of the gradient of the dual function. For computational purposes, we may avoid difficulty with the strict inequality in the definition of the region R by placing a lower bound E on Xi, i = l---n, where E is a small positive constant chosen to ensure that no difficulty arises through rounding errors. In the region R, sufficient conditions for a convex dual problem are satisfied so that any optimum found in step 4 of the above procedure will be global over R. The constraints which were used to define R still enable i(A) to be written down term by term. It is possible to incorporate any form of lower and upper bound in each variable without losing this ability, but this ability is immediately lost if a constraint used to define R contains terms in more than one variable. Such constraints must be adjoined to the Lagrangian, each requiring an additional Lagrange multiplier, and it is evident that any advantage in dimensionality of the dual problem is lost if there are a large number of these. Note that the upper bound, ximax., keeps each ii (A) finite and thus the corresponding diameter greater than zero. Note, also from Fig. 1, that for each path j, there is at least one pipe section i (the section adjacent to the valve) which is on that path
only, i.e. there is at least one i for which
j=l
has only the one non-zero ati term. If the bounds x.~~~., i= l---n were absent, then at the optimum all hj, j = 1 ---m, would be non-zero, which means that all the constraints of the primal problem given by Eq. (14) would be active. The physical interpretation of this is that the pipe diameters would be squeezed down to utilize the maximum pressure drop available in each in ridiculously small path, often resulting diameters in the sections adjacent to the valve. We have thus far considered the kinetic energy correction factor, Cii = 1 ---n, to be constant. While this is true for incompressible fluid flow, it will be seen from Eq. (2) that for compressible fluid flow each C.depends in a nonlinear manner on both the diameter and terminal pressures of each section. It is therefore necessary in this particular application to construct an outer iterative cycle (see Fig. 3) in which the values of Ci are recalculated using the optimal diameters and the resulting pressure drops. As the reviewer has pointed out, it is not possible to prove that such a cycle will converge to the correct solution, or that it will converge at all. This difficulty, which is present in both the primal and dual formulations, is not just a technical matter of being unable to provide a convergence proof; the nonlinearity of the correction factor Ci makes even observed convergence of this outer iterative cycle subject to doubt as to the validity of the converged solution, and justification for accepting it must be supplied from other considerations. Along this line, it is possible to obtain an accurate measure of the actual discrepancy in the constraints at each iteration; given the diameters of each section and starting with the known flare pressure, P,, it is possible to show that the sequential calculation of pressure by the iteration P.r+l 1 =f(P,‘) for each section, where
1138
(36)
An approach to the optimal design of networks
X [
Pi’ 1 + Kidi”“14 log10 PI_, -
(
l/2 )I)
(37)
and where Ki = 4’61/4fiLidi”‘*14 will converge for all values of Pi-1 and di in the range of interest. Thus it is possible to calculate the actual back pressure at each valve, say Pj, and determine the discrepancy, pj-Pj spec., in the total pressure drop for those constraints j which are active. It is also possible to obtain an estimate of the penalty incurred in the objective function for each of these. Let & = I’j’- PO’, and Ab = &jbj. A well-known economic interpretation of Lagrange multipliers (see for example [2]), which applies equally to linear and nonlinear inequality constraints, is that hj is a measure of the extent to which the optimum value of the objective function, say f, could be further decreased by increasing the bound bj on the constraint, i.e. hj = (a! (E)/abj). We may thus get an approximate estimate of the effect of each discrepancy on the objective function by the term XjAbj. (This is only an approximate estimate, since we are assuming that x approximates the solution to the problem with the bound 6, on each constraint j, and that Abj is sufficiently small for a first order estimate to be adequate.) Since in practice the correction caused by Ci is itself small[5], it may be satisfactory merely to know the actual discrepancy in pressure drop. 6. COMPUTATIONAL
EXPERIENCE
Although emphasis has been placed on the approach rather than the application, this paper would not be complete without some description of computational experience. The optimization method used to solve the problems is a recently developed ‘variable-metric’ projection method, so called because it accumulates and utilizes an approximation to the local inverse matrix of second partial-derivatives of the objective function. This method has been used very success-
fully in unconstrained minimization ([ 14, 151); the basis of the approach in the presence of linear constraints is the use of projection operations which confine the search to the subspace formed by the active constraints while still retaining complete curvature information. A full description and theoretical justification of the method may be found in Refs. [6] and [ 161. Two practical points should be mentioned at this stage. The first is that, because of the various classifications of pipes which are made in practice, the number of pipe sections n is much greater in comparison to the number of paths m than that implied by the definition in Section 2. Also, it is possible to include an ‘equivalent length’ correction for elbows and bends by including this in the term Ci; the same remarks made in the previous section concerning convergence of the iterative adjustment of Ci apply equally here. The primal problem formulated in Section 3 is quite tractable; the main motivation for continuing on to formulate the dual problem was the fact that relief header networks of well over a hundred pipe sections are becoming more common as refineries increase in size and complexity. In order to make a comparative evaluation of the primal and dual formulations a number of medium sized problems taken from previous contracts were solved using both. The dual problem was consistently the most easily solved, both in terms of number of function calculations required and total execution time. For example, in one problem consisting of 3 1 pipe sections and 9 valves, the primal formulation required 35 function calculations and took 8.15 set computing time on a CDC 6400, whereas the dual formulation required 10 function calculations and took 1.31 sec. The maximum discrepancy in pressure drop resulting from the assumed and calculated values of Ci (including an equivalent length term) was 1 per cent after two iterations for this problem. The solution obtained is subjected to a postoptimal treatment, which selects the appropriate commercial pipe diameters by a partially enumerative search and also checks for critical
1139
B. A. MURTAGH
velocity in each pipe section. The program is now in full use as a design tool and has already achieved considerable savings in capital cost on a number of contracts. The largest system which has been solved consisted of a total of 225 pipe sections. 7. DISCUSSION AND SOME GENERALIZATIONS
The two reasons why the dual formulation turns out to be so successful computationally are the fact that the constrained sub-problem is solved analytically, and that the number of nontrivial constraints in the primal problem is much less than the number of variables. Thus we not only reduce the number of variables by formulating the dual problem, but also are able to avoid the -usual handicap of having a “minimization within a minimization.” The approach presented in this paper is not restricted to fluid-flow systems, but can be applied to any system whose configuration can be represented by a network. By observing the analogies in other contexts of the three basic variables: pressure drop or change in potential, flowrate, and diameter or resistance to flow, similar optimization problems may be formulated in such diverse areas as transportation, scheduleconomics, and of course engineering im systems. The approach can be applied in a straightforward manner to equations of the form: (Potential Drop)‘] = (Flow)j2 X Resistance)j3. Any departure from this form by the use of correction factors, as in the example presented here of compressible fluid flow, would require further treatment and physical reassurance that the contribution was small. Note also that the equations which describe the system being optimized are assumed to be nonlinear; in the case of linear systems the use of linear programming is already well-established. One may often be, as in this example, dependent on the nonlinearity of the flow equations with respect to the independent variables in order to obtain a strictly convex primal objective function. Although the network considered here consisted of multiple sources but only one sink,
there is no more difficulty in considering multiple sources and sinks as long as each path between source and sink is not ambiguously defined and the flowrate in each section is either specified or calculable. The number of sections will always be greater than the number of paths, regardless of the configuration of the network, so it would seem that the dual formulation would always be an attractive alternative in such applications. Acknowledgements-The author wishes to thank Fluor Corporation for permission to publish this paper. Although mentioned only briefly in the text, a considerable effort was expended by Manuel de la Puente and Chan-shan Fan in developing a program for selecting the appropriate commercial pipe sizes and checking the resultant pressure drops rigorously. The author also wishes to thank the reviewer for helpful comments and constructive criticism.
NOTATION
a,
coefficients of constraint equations (see Eq. 8) bj constant term of constraint equations (see Eq. 8) kinetic energy correction factor (see ci Eqs. 2 and 4) di diameter of pipe section i, in. dimin. minimum diameter of pipe section i, in. En euclidean n-space i subscript to denote pipe section j subscript to denote path from source to sink (and, in the problem presented here, the corresponding valve) & constant term in pressure-drop equation for pipe section i (see Eq. 2) ~~ length of pipe section i, ft p,, flare pressure, psia exit pressure of section i, psia pi P5spec. specified maximum pressure at valve j, psia R region in En (see Eq. 28) sj the set of pipe section numbers in path j weight per unit length of pipe section i, wi lb/ft independent variable in primal minimXi ization problem X;(A) unconstrained optimal value of x,, given X (see Eqs. 29 and 30)
1140
An approach to the optimal design of networks X,max.
I
maximum
value
of
xi :ximaX.=
~/(~imin.)4+314
ii(X) f
optimal value of Eq. 3 1) optimal value of x
Xi
for a given A (see
Greek symbols
(Y,p
coefficients
(see Eq. 1) dual function (see Eq. 23) 4 Lagrangian function (see Eq. 16) Xj Langrange multiplier associated with constraint j of primal minimization problem. Independent variable in dual problem. X optimal value of A y
relating weight of pipe to dia. REFERENCES
VI FULKERSON
PI HU T. C.,
18. and Network Flows. Addison-Wesley,
D. R., .I. Sec. Ind.
Integer Programming
Appl. Math. 19619
Mass. 1969. P., Notes to accompany lectures, NATO Summer School on Integer and Nonlinear Programming, Bandol, France 1969. API-RP520, Recommended Practice for the Design and Installation of Pressure-Relieving Systems in Refineries: Part I-Design. p. 22, Third Edn., American Petroleum Institute 1967. PERRY J. H., ChemicalEngineer’s Handbook, 4th Edn., pp. 5-24. McGraw-Hill, New York 1963. MURTAGH B. A. and SARGENT R. W. H., Optimization, pp. 215-246, (Ed. R. FLETCHER). Academic Press, London 1969. KUHN H. W. and TUCKER A. W., Proc. Second Berkeley Symposium on Mathematical Statistics and Probability, pp. 481-492, (Ed. J. NEYMAN), University of California Press, Berkeley 1951. FALK J. E.,/. Math. Anal. Appt. 1967 19 141. HORN F. and JACKSON R., Ind. Engng Chem. Fund/s. 1965 4 110; 1965 4 239; 196.5 4 387. JACKSON R., Chem. Engng Sci. 1964 19 19; 1964 19 253; 1965 20 405. KATZ S.,A.I. Ch.E.-Inst. Chem. Engrs. Symp. Ser. 1965 4 69. MANGASARIAN 0. L. Nonlinear Programming. McGraw-Hill, New York 1969. BROSILOW C. and LASDON L.,A.Z. Ch.E.-Inst. Chem. Engrs. Symp. Ser. 1965 4 75. FLETCHER R. and POWELL M. J. D., C0mpt.J. 1963 6 163. MURTAGH B. A. and SARGENT R. W. H., C0mpt.J. 1970 13 185. SARGENT R. W. H. and MURTAGH B. A., paper presented at the 7th International Symposium on Mathematical Programming, The Hague, Sept. 1970.
131 WOLFE [41 [51 [61 [71 [81 [91 [lOI u 11 t121 [I31 [I41 t151 [I61
R&urn&-
Le probleme de conception trait& dans cette etude se rapporte 21la distribution optimale de la chute de pression dans un reseau oh le debit est d&ermine mais oti les diametres des tubes peuvent varier. Si l’on utilise des elements de la theorie de dualite de programmation non-lineaire, le probleme semble pouvoir etre resolu a I’aide dun calculateur. Les auteurs presentent une approche qui a une applicabilite generale aux reseaux dans lesquels les equations de I’ecoulement sont non-lineaires. Zusammenfassung-
Das in diesem Artikel behandelte Konstruktionsproblem betrifft die optimale Verteilung des Druckabfalls iiber ein Netzwerk, in welchem die Stromungsgeschwindigkeiten festgelegt sind, w&end die Rohrdurchmesser variieren kiinnen. Unter Verwendung der Dualitatstheorie nicht-linearer Programmierung kann das Problem in eine computermassig zugangliche Form gebracht werden. Die dargelegte Methode ist von allgeneiner Anwendbarkeit auf Netzwerke, in welchen die Stromungsgleichungen nicht-linearer Art sind.
1141