- Email: [email protected]
Nonlinear programming for multiperiod capacity planning in a manufacturing system
EUROPEAN JOURNAL OF OPERATIONAL RESEARCH
ELSEVIER
European Journal of Operational Research 96 (1996) 167-179
Theory and Methodology
Nonlinear programming for multiperiod capacity planning in a manufacturing system Kurt M. Bretthauer *, Murray J. Ctt6 Department of Business Analysis and Research, Texas A & M University, College Station, TX 77 843, USA
Received July 1995
Abstract In this paper, we present nonlinear programming methods for capacity planning in a manufacturing system that consists of a set of machines or work stations producing multiple products. We model the facility as an open network of queues where capacity at each work station in the system may be changed in each of a finite number of time periods. To determine the timing and size of capacity changes, we present two nonlinear programming models and methods for solving the resulting problems. One model involves minimizing total capacity costs such that plant congestion is controlled via upper limits on work-in-process. The other model involves minimizing a weighted sum of product lead times subject to budget constraints on capacity costs. We present solution methods for continuous and discrete capacity options and convex and nonconvex (e.g., economies of scale) capacity cost functions. We use branch and bound and outer approximation techniques to determine globally optimal solutions to the nonconvex problems. Computational testing of the algorithms is reported. Keywords: Capacity planning; Nonlinear programming; Branch and bound; Nonconvex optimization
I. Introduction In this paper, we present nonlinear programming methods for capacity planning in a manufacturing system. Consider a facility consisting of a set of machines or work stations producing multiple products where the sequence of stations visited depends on the processing required by a particular product type. We are interested in determining the timing and size of capacity changes at each work station over a finite number of time periods. The capacity planning decision must take into account the product demand forecasts, existing capacity levels, the costs of changing capacity, the impact of capacity changes on the * Corresponding author.
performance of the manufacturing system, and the firm's manufacturing strategy. During periods of increasing demand it may become necessary to add capacity to some or all of the work stations to satisfy the demand and maintain stable operating conditions in the plant. If demand is expected to decrease for a long period of time, then it may be cost effective to reduce capacity. We model the manufacturing system as an open network of queues and present an optimization framework for capacity planning over a multiperiod planning horizon. The decision variables are the service rates (capacity) at each work station in each time period. Capacity can be controlled at the work stations via the number of machines, modernizing or updating equipment, additional maintenance, the
0377-2217/96/$15.00 @) 1996 Elsevier Science B.V. All fights reserved P// S0377-2217(96)00061-6
168
K.M. Bretthauer, M.I. C~t~ / European Journal of Operational Research 96 (1996) 167-179
number of workers, the number of shifts, the use of overtime, providing additional worker training, etc. We present two nonlinear programming models for capacity planning. Solution methods are developed for continuous and discrete capacity options and for convex and nonconvex (e.g., economies of scale) capacity cost functions. By incorporating stochastic queueing behavior at the work stations within the optimization framework, the models not only identify capacity levels such that forecasted demand in each period will be satisfied, but they also control performance measures such as product lead times, work-in-process, and capacity costs. Specifically, one model involves the minimization of capacity expansion costs minus capacity reduction savings (or plus capacity reduction costs, if appropriate) subject to upper limits on work-in-process. The other model involves the minimization of a weighted sum of product lead times subject to budget constraints on capacity costs. Also note that the capacity planning models presented here can be applied to other queueing network applications, such as computer/communication systems [14] and health care systems [15]. Bitran and Tirupati [7,8] consider a similar problem involving a one time, or single period, capacity change decision in a manufacturing queueing network. They address only the size of capacity changes at the work stations; they do not consider the timing aspect of the decision. Bitran and Tirupati present a model that minimizes capacity costs subject to a single work-in-process constraint. In [7], they treat capacity as a continuous variable and present a heuristic solution technique for expanding capacity. They briefly discuss how their method can be modified to also handle capacity reductions. In [8], they present a heuristic for the capacity expansion problem when capacity is available in discrete options. See Bitran and Sarkar [5,6], Boxma, Rinnooy Kan and Van Vliet [9], Bretthauer [10,11], Buss, Lawrence and Kropp [13], and Van Vliet and Rinnooy Kan [24] for related single period models. Karmarkar, Kekre and Kekre [19] and Solberg [23] also model the manufacturing system as a queueing network and discuss capacity issues, but do not present an optimization model for making capacity change decisions. Although it is relatively straightforward to extend
the single period models of Bitran and Tirupati [7,8] to the multiperiod case, solving the resulting nonlinear optimization problems can be quite difficult. As will be discussed in later sections of the paper, these nonlinear optimization models often will be nonconvex. Therefore, standard nonlinear programming methods may converge to only a locally optimal solution, not necessarily a globally optimal solution. Solving the capacity planning models may be further complicated by situations where capacity is available only in discrete options. In this paper, we present methods for determining globally optimal solutions to these types of nonconvex multiperiod problems, including the case of discrete capacity options. Bitran and Tirupati [7,8] present heuristics for the single period models. See [1,2,3,12,17,18,21] for a detailed discussion of global optimization algorithms and the related literature. Much of the capacity planning literature considers the long term problem, such as when and where to build additional facilities [20]. Here, we focus on a medium term planning horizon where capacity may be changed at the individual work stations in an existing facility. Other medium term capacity planning research considers methods for estimating how much capacity will be necessary at each work station to satisfy a given master production schedule [4]. However, these approaches do not present a decision model for adjusting capacity levels if a mismatch between capacity required and capacity available is identified. Also, these approaches do not take into account the relationship between capacity levels and system performance measures such as work-in-process and product lead times. Our models and solution methods provide recommendations for the size and timing of capacity changes that account for capacity costs and the impact of capacity levels on the performance of the system. The rest of the paper is organized as follows. In Section 2, we present two versions of a capacity planning model that involves minimizing capacity costs subject to work-in-process constraints. In Section 3, we present capacity planning models that involve minimizing product lead times subject to budget constraints on capacity costs. In Section 4, we discuss solution techniques for the continuous capacity, convex cost version of the models. Section 5 provides a branch and bound/outer approximation
K.M. Bretthauer, M.I. C~t~ / European Journal o f Operational Research 96 (1996) 167-179
algorithm for solving the discrete capacity, nonconvex cost version of the problem. In Sections 6 and 7, we provide computational results for continuous capacity, convex cost problems and discrete capacity, nonconvex cost problems. Section 8 contains some concluding remarks.
2. The minimum cost capacity planning model We consider two versions of a minimum cost capacity planning model that minimizes capacity costs subject to work-in-process constraints. Which version is appropriate depends on the manufacturing system under consideration, as discussed in the following subsections.
2.1. Simultaneous capacity increases and decreases permitted Problem (C1) permits the addition of new capacity and the elimination of old capacity in the same time period at a given work station. This option may be attractive if existing capacity at a station can be replaced with new cheaper or more efficient capacity (possibly a new technology). (C1) P
Min
P
~',fit(xit) - ~ iffil t=l
s.t.
i=1
(1)
Y'~gi,(Yit) t=l
~vitLit(lxit)<-..WlPt,
t~T,
(2)
i=1
I'Zit = [£i,t- 1 ~- Xit -- Yit,
i ~ N,
t~
T,
(3) (4)
O<~xlt<~cit,
i ~ N, t ~ T, i~N, t~T,
0 <~Yit <~ dit,
i ~ N,
(6)
air <~ ~it <~ bit,
t ~ T.
(5)
Subscripts i t N T
= = = =
Work station subscript, i = 1. . . . . n. Time period subscript, t = 1. . . . . p. { 1, 2 . . . . . n} = Index set for work stations. {1, 2 . . . . . p} = Index set for time periods.
Decision variables l,Zit Capacity (service rate) at work station i in ~--
period t.
169
xit = Increase in capacity at work station i available at the beginning of period t. Yit = Decrease in capacity at work station i at the beginning of period t.
Functions fit(x~t) = Cost of increasing capacity at work station i in period t by xit. ggt(Yit) = Savings from decreasing capacity at work station i in period t by Yir Lit(I,l~it) = M e a n number of jobs at work station i in period t as a function of /zit (discussed further below).
Parameters v,.t = Average dollar work-in-process value per job at work station i in period t. WlP t = An upper limit on total work-in-process allowed in the facility in period t. /xi0 = Initial capacity at work station i. Ais = Arrival rate (demand) at work station i in period t (assume A,.t > 0 for all (i, t)). ai,, bit = Lower and upper bounds on capacity at work station i in period t ( a , > A,). cit = Upper bound on capacity increase at work station i in period t. dit = Upper bound on capacity decrease at work station i in period t. ca,. t = Squared coefficient of variation (variance of a random variable divided by the square of its mean) of the interarrival times at work station i in period t. cs, = Squared coefficient of variation of service times at work station i in period t. Note that for all models presented in this paper, the capacity reduction savings function git(Yi,) can easily be modified to handle settings where capacity reductions result in costs instead of savings. Problem (C1) requires the ability to calculate system performance measures in a network of queues (in particular, the number of jobs at each station). To make the model as flexible as possible, we consider general interarrival time distributions and general service time distributions at each single server station. We use the approximate decomposition approach to analyze this class of queueing networks
K.M. Bretthauer, M.J. C~td/ European Journal of OperationalResearch 96 (1996) 167-179
170
(see [7,8] and the references therein). The parameters c a . and cs/t are assumed to be constant with respect to capacity. See Bitran and Tirupati [7,8] for justification of this assumption and a discussion of the system of equations to be solved to determine the parameters A., ca~t, and cs/t for all (i, t). We also use a constant arrival rate Ait within a particular period t. This assumption is appropriate when changes in demand occur at discrete points in time, such as when a contract is received from a new customer. The queueing results can also be expected to provide good approximations when demand changes continuously at a slow rate (see Hall [16, pp.178-180] for further discussion of this issue). Constraint set (2) imposes an upper limit on total work-in-process in the facility in each time period. Therefore, we must be able to calculate the average number of jobs Li,(jzit) at each work station as a function of the station capacity. Bitran and Tirupati [7] use the following approximation to calculate Lit(/~,) and show that it is a convex function of/z~:
t.,,( it.) = X,t/u. +
( c a i t + csit ) A~t
21.~it( I~it- Air )
h( X., #i,. ca.. cs.),
where
h('Ait, ~£it' ca/t,
CSit)
[ - 2(1
=, exp[
-
cait)(
Problem (C2): Same as problem (C1) except the following constraints are added:
xitYit=O,
iEN,
(7)
tET.
Constraint set (7) eliminates the possibility of increasing and decreasing capacity at a given work station in the same time period. However, adding these constraints can make problem (C2) much more difficult to solve than (C1), as is discussed later in the paper.
3. The minimum lead time capacity planning model The next capacity planning model that we consider involves minimizing a weighted sum of product lead times rather than minimizing capacity costs. To control capacity costs, we include constraints that limit the amount of money that can be spent on capacity. Minimizing lead times is attractive when management wants to focus on increasing manufacturing flexibility and improving customer service.
I~it- A/t )
3a~(-~a. + cs/~ 1
a given work station are not permitted. For example, a company may not wish to, or be able to, replace experienced workers with new inexperienced workers, even though this strategy may reduce salary expenses. Problem (C2) is as follows:
if c a , < 1,
3.1. Simultaneous capacity increases and decreases permitted
if c a i t > 1,
To handle multiple servers at a station, the expression for Lit given in Bitran and Tirupati [8, pp.134135] can be used in place of the single server equation given above. In this case,the number of servers would be the decision variable instead of the service rate, Constraint (3) is a set of 'capacity flow balance' equations. Constraint sets (4), (5), and (6) set lower and upper bounds on the capacity variables.
The minimum lead time model where increases and decreases in capacity are allowed at a given station in the same time period is formulated as follows. (L1)
(8)
Min s.t.
~L,(x.)i=1
2.2. Simultaneous capacity increases and decreases not p e r m i t t e d
Next we consider environments where increases and decreases in capacity in the same time period at
g.(yi,)~B,,
~it = l~ia - I + x . - )'it,
air <~ I,~it <~bit,
ter,
i=l
i ~ N,
i ~ N, t ~ T,
O <~ xit <~ cit,
i ~ N,
t ~ T,
O <~ Yit <~ dit,
i ~ N,
t E T.
t E T,
(9)
K.M. Bretthauer, M.J. C~M/ European Journal of Operational Research 96 (1996) 167-179
The objective function W(/~) represents a weighted sum of product lead times as a function of the vector/x of work station capacities. Lead times are calculated from Little's equation using the expression for Lit(i~it) previously described. A particular product's total lead time through the system is calculated by summing the time spent at each work station on its route. As in Karmarkar, Kekre and Kekre [19], we weight each product's lead time by its demand when calculating W(/.~). The lead times could also be weighted by cost, to allow for low demand but expensive to hold products. Constraint set (9) sets a capacity budget of B t dollars in period t, t = l . . . . . p. The remaining constraints are the same as in problem (C 1).
3.2. Simultaneous capacity increases and decreases not permitted Next we modify the minimum lead time model (L1) to handle environments where capacity cannot be increased and decreased at a given station in the same time period. Problem (L2) is as follows: Problem (L2): Same as problem (L1) except the following constraints are added:
)fit Yit =
0,
i ~ N,
t ~ T.
4. Continuous capacity, convex cost solution methods Here we address settings where capacity can be treated as a continuous variable and make the following assumption concerning the capacity cost and savings functions. Assumption 1. The cost function fit(x.) is convex with respect to xit for all (i, t) and the savings function g.(Y~t) is concave with respect to Yit for all (i, t). This assumption is appropriate for f , t ( x , ) when the cheapest capacity expansion alternatives are selected first and each subsequent addition may become more expensive on a per unit basis. The as-
171
sumption is appropriate for git(Yit) when the capacity reductions yielding the greatest savings are selected first and each subsequent reduction may yield less savings on a per unit basis. If Assumption 1 is not satisfied (e.g., a concave capacity cost function to handle economies of scale), then the solution obtained may be only locally optimal instead of globally optimal. This issue is discussed later in the paper. When Assumption 1 is satisfied and capacity is a continuous variable, models (CI) and (L1) are convex optimization problems and, therefore, can be solved with a commercially available nonlinear programming package.
4.1. Solving a special case of model (C2) Consider model (C2) for the case where capacity is a continuous variable and Assumption 1 is satisfied. Because of the constraints xitYit = 0 for all (i, t), problem (C2) is a nonconvex optimization problem. Therefore, standard convex programming techniques may converge only to a locally optimal solution instead of a globally optimal solution. Fortunately, we are able to show that with some very reasonable additional assumptions on the cost functions, the constraints x . y . = O for all (i. t) will automatically be satisfied and, therefore, can be eliminated from the formulation. Consider the following assumption: Assumption 2. f/t(0)= 0; git(O) = 0; f/,(x/t) is an increasing convex function of xit; git (Yit) is a nondecreasing concave function of yi,; and f / t ( f l ) > g . ( f l ) for all /3 > 0. In addition to f.(xit) being convex and g . ( y . ) being concave, we are also assuming that there are no costs or savings if capacity is not changed, that additional capacity always costs something, that additional capacity reductions always yield a savings or salvage value greater than or equal to zero, and that the cost of obtaining a given amount of new capacity is more than the savings or salvage value from reducing capacity by the same amount at a given station in a given time period.
K.M. Bretthauer, M.I. COtd/ European Journal of Operational Research 96 (1996) 167-179
172
R e m a r k 1. Assumption 2 implies
Then,
fit(xit) - f i t ( x i t - 3/) > g i t ( Y ) - git(0)
z*--~.=[f(x*)-g(y*)]-
for all xi, >1 y > O.
=
t Xit "i,t~R
R e m a r k 2. Assumption 2 implies
for all Yit ~ "Y > O.
Proof. Let Q = {(i, t) [ x:t Yi; = 0}, let
--
fit Xit--Yit i,
R
( E g i t ( Y i : ) -- E git( 0 ) } i,tER i,t~R
L t ( ' Y ) --fit(O) > git( Yit) -- git( Y i t - "Y)
Proposition 1. If Assumption 2 is satisfied, then the optimal solution Ix:t, x:t, Yi7 for all (i, t) to the convex problem (C 1) is a globally optimal solution to the nonconvex problem (C2).
[f( 2) -- g( ~)]
idES
idES
idES idES Because we have assumed R U S ¢ 0, the above expression for z * - ~ along with Remarks 1 and 2 imply z * - z > 0. The inequality z * - ~ > 0 contradicts the fact that xi,, Yit, Ixit for all (i, t) is the optimal solution to problem (C1). Therefore, we must have R U S = I ~ . []
and let
Proposition 1 implies that we can solve the convex problem (C1) instead of the much more difficult nonconvex problem (C2) when Assumption 2 is satisfied.
S = {(i, t) I x:t Yi7 ~ 0 and xit < Yit }.
4.2. Solving a special case of model (L2)
R = {( i, t) I get Yi7 ~ 0 and x:t >1YeT },
We must show that x:t Yi~ = 0 for all (i, t) to verify that Ixit, xit, Yit for all (i, t) is the globally optimal solution to problem (C2), i.e., that R U S = 0. The proof will be by contradiction. Assume, to the contrary, that R t,J S 4: 0. Consider the solution Txit = Ix:t for all (i, t), 2it = x:t and S:it = Yi~ for all (i, t) ~ Q, 2it = Xit -- Yi~ and Pit ~" 0 for all (i, t) ~ R, and 2i, = 0 and Yi, = Yi~ - x:t for all (i, t) E S. Note that the solution ~i,, Ycit, Pit is feasible to problems (C1) and (C2). Let f ( 2 ) = Y'-fit(xit) i,t
and g ( Y ) = E g i , ( ~Yi,). i,t Let f ( x* ) and g( y * ) be defined similarly. Also, let =f(2)
-- g ( y )
Here we consider model (L2) for the case where capacity is continuous and Assumption 2 is satisfied. Let IX¼, x:t, and Yi~ for all (i, t) now denote the optimal solution to the convex problem (L1) and let Q = {(i, t) I x:t Yi; = 0}, g = {(i, t) I x:t Yi7 ~ 0 and xit >~Yi, },
and S = {(i, t) lx:t Yi7 ~ 0 and xi, < Yi~ }. Similar to Proposition 1, consider the following solution: T-~it~-~lbi; for all (i, t), 2it = Xit and ~:it = Yi, for all (i, t ) ~ Q, 2it = x:t -Yi~ and Pit = 0 for all ( i , t ) E R , and 2 , = 0 and P i t = Y i ~ - x : t for all
(i, t ) E S . Proposition 2. If Assumption 2 is satisfied, then Ixit, xit, and Pit for all (i, t) is a globally optimal solution to the nonconvex problem (L2).
and
Proof. If R U S = 0, then Pit ~-~£it' 2it ~-Xit' and Pit = Yi~ for all (i, t) is clearly a globally optimal
Z* = f ( x * ) - - g ( y * ) .
solution to problem (L2). If R tJ S 4= ¢ and Assump-
K.M. Bretthauer, M.L C~t~/ European Journal of Operational Research 96 (1996) 167-179
tion 2 is satisfied, then, as in the proof of Proposition 1, we can show that [f(x')-g(y*)]-
[ f ( 3:) - g( ~)] > 0 .
This inequality implies that the solution izit, xit, and ~i, for all (i, t) is feasible to problem (L2). We also know W ( ~ ) = W(/z* ) because ~i, =/zit for all (i, t). Therefore, izit, xit, and 3it for all (i, t) is a globally optimal solution to problem (L2). [] Proposition 2 shows that if Assumption 2 is true, then a globally optimal solution to the nonconvex problem (L2) is easily constructed from the solution to the convex problem (L1). 4.3. A branch and bound algorithm f o r models (C2) and (L2)
Next we discuss a branch and bound method for finding a globally optimal solution to the nonconvex problems (C2) and (L2) for the case where capacity is continuous, Assumption 1 is satisfied, and Assumption 2 is not satisfied. The algorithm relaxes the constraints xitYit = 0 for all (i, t) and solves a convex nonlinear subproblem of the form of problem (C1) (or (L1)) at each node of the search tree. Branching is performed on a selected subproblem when its solution does not satisfy constraint set (7), i.e., when xitYit 4:0 for some (i, t). Consider such a subproblem where xitYit ~ 0. A new left branch subproblem is constructed by adding the upper bound xit ~< 0 to the selected subproblem. A new right branch subproblem is constructed by adding Yit <~0 to the selected subproblem. These two new subproblems are solved and added to the candidate list of unfathomed subproblems. Otherwise, standard branch and bound logic and fathoming rules are used. The algorithm terminates after a finite number of iterations with the globally optimal solution to problem (C2) or (L2).
5. Discrete capacity, nonconvex cost solution method In this section, we present a method for solving models (C1) and (C2) for the case of a nonconvex capacity cost (or savings) objective function and
173
discrete capacity options. Examples of discrete capacity variables include the number of machines, shifts, or workers. If costs exhibit economies of scale, then a nonconvex cost function is appropriate. The solution method combines branch and bound with outer approximation techniques. In Bretthauer [11], we presented a version of this algorithm for solving single period problems with a concave objective function, a set of convex nonlinear constraints, and discrete variables (see [12] for the special case where all constraints are linear). Here we extend the approach in [11] to handle a more general nonconvex objective function as well as the nonconvex constraints xitYit = 0 for all (i, t). The algorithm presented here also takes advantage of the fact that some of the constraints in (C 1) and (C2) are linear. In addition, we report computational results with the algorithm, which was not done in [11]. See Benson [2] for further background on branch and bound/outer approximation methods for concave cost, continuous variable problems. First we will present the algorithm for solving problem (C1), which does not include the constraints xit Yit = 0 for all (i, t). Then we will show how to modify the algorithm to handle these nonconvex constraints for problem (C2). To simplify the notation and presentation, we describe the algorithm for a general optimization problem with the structure of model (C 1). Consider the following problem:
(P) n
Min
E f/(Xi) i~l
s.t.
x~GNPNH,
x discrete. In problem (P), x = ( x I . . . . . Xn) is an n-vector of decision variables, f,.(Xi) , for i = 1. . . . . n, is a realvalued cost function (potentially nonconvex), G is a convex set defined by m convex nonlinear constraints gj(x)~< 0, for j = 1. . . . . m, P is a polyhedral set defined by a set of linear constraints, and H-
{ x ~ Rn I / i < x , < u i, i = 1 . . . . . n}
K.M. Bretthauer, M.I. COt~/ European Journal of Operational Research 96 (1996) 167-179
174
is a set of finite and integer valued lower and upper bounds on x. Also, let f(x)
= ~ f~(xi)
mate f ( x ) , leading to tighter lower bounds. The outer approximation cuts help force the subproblem solutions to satisfy the relaxed nonlinear constraint set G.
i=1
and, for notational convenience, let g ( x ) = m a x { g j ( x ) l j = 1 . . . . . m}. Note that g ( x ) is a convex function. Then G can be written as
G= { x ~ R"lg(x) <~0}. We assume integer decision variables. The following discussion remains valid for other discrete sets of variable values. To write the general problem (P) in terms of model (C1), the vector x would include lZit, x u, and Yit for all (i, t), G would include the nonlinear work-in-process constraints (2), P would include the linear capacity flow balance constraints (3), and H would include the bounds on the variables (4), (5), and (6). As in standard branch and bound methods, the algorithm generates a search tree that requires the solution of one subproblem at each of its nodes. To allow the calculation of lower bounds on the optimal objective value of the original problem (P), the subproblems are constructed using three types of 'relaxations': 1) the discrete variable requirements are relaxed; 2) a convex underestimating function to the nonconvex objective is used; and 3) the convex nonlinear constraints in the set G are replaced by a series of linear outer approximation cutting planes. These three relaxations result in convex programming subproblems with all linear constraints, making their solution much simpler than if the nonconvex objective was maintained. We are particularly interested in cost functions that exhibit economies of scale, i.e., concave cost functions. In this case, the underestimating function is linear, resulting in linear programming subproblems. In every iteration of the algorithm, either a subrectangle of H is divided into two smaller subrectangles (rectangular partitioning, or branching), yielding two new subproblems, or an outer approximation cutting plane is added, yielding one new subproblem. Rectangular partitioning helps force the variables to take on discrete values and also causes the convex underestimating function to more accurately approxi-
5.1. Rectangular partitioning Let H={X~nlii
<~Xi <~Ui, i =
1..... n}
denote any rectangular subset of the initial rectangle H. Also, let 2 denote a point in H. First consider the case where 2 has a fractional component 2q. We divide /~ into the two subrectangles,
-' = (x, filx,,
t2,,l)
and
,-,+ = (x +,+ x,, >- l+,,l), where [. ] and [. ] represent the integer round-down and round-up operations, respectively. Now consider the case where all elements of 2 are integer valued. We select acomponent 2q with lq < fiq. In this case we divide H into
H1= ( x E H I Xq<<.2q-1)
H2=(xEHIXq>/2q}. (Note that if 2 o = lq, then instead we divide H into
nl={xE-Hlxq<~2q}
and n 2 = { x ~ H I X q ~ 2 q
+ 1}.) If H is a single point, then no partitioning will be performed and the subproblem is fathomed.
5.2. Convex envelopes As previously mentioned, in the subproblems we replace the nonconvex objective f ( x ) with a convex underestimating function. In particular, we replace f ( x ) with its convex envelope over a given rectangular partition element. Definition 1 (Convex envelope [18]). The convex envelope of a function f ( x ) taken over a nonempty
K.M. Brenhauer, M.J. C~td / European Journal of Operational Research 96 (1996) 167-179
subset M of the domain of f ( x ) is a function c(x) such that: l) c(x) is a convex function defined over the convex hull of the set M. 2) c(x)<~f(x) for all x ~ m . 3) If h(x) is any convex function defined over the convex hull of M and if h(x) <~f(x) for all x E M, then h(x) ~< c(x) for all x in the convex hull of M. The convex envelope is particularly attractive as an underestimating function to a separable concave function because then it is easy to construct and linear. We will construct convex envelopes ~(x) of f ( x ) over H, which satisfies ~(x)<<,f(x) for all zEH.
5.3. The outer approximation cutting planes As the algorithm proceeds, a series of linear outer approximation cutting planes are added in an attempt to force the subproblem solutions to satisfy the relaxed nonlinear constraints. A new cut is added only when a subproblem is selected whose solution is integer and does not satisfy the nonlinear constraints in the set G. Let ~ denote a subproblem solution that is integer, but not in the set G. In this case, we construct a linear outer approximation cutting plane r(x) ~ 0 that 'cuts off' ~ from G. Consider the line segment connecting the point ~ with some point y in the interior of G. The first step in constructing the cut requires determining with a line search the unique point w on this line segment that intersects the surface of G. Let q denote a subgradient of g(x) at the point w. Then, the affine function r( x) = ( x - w, q) + g( w) yields the outer approximation cut r ( x ) ~< 0 that cuts off the point 2 from G, where ( a , b) denotes the inner product between the vectors a and b.
5.4. The continuous subproblems Given a rectangular subset H of the initial rectangle H, the convex envelope ~(x) of f ( x ) over H, and the set of s outer approximation cuts rw(x) ~< O, w = 1 , . . . , s, constructed up to the current iteration of the algorithm, the convex subproblems will be of the following form:
175
( P3 Min S.t.
rw( X ) <~0,
w=
1.....
s,
x~PnH. Note that when f ( x ) is concave these will be linear programming subproblems. Because of the relaxations used, the optimal objective value of subproblem ~ provides a lower bound on the minimum of f(_x) over the integer values of x in the set G f~ P N H.
5.5. Algorithm description Let LB k denote the lower bound on the optimal objective value of problem (P) at iteration k of the algorithm, let UB k denote the best upper bound found up to iteration k, and let x c denote the incumbent feasible solution to problem (P) providing upper bound LIB k. In iteration 0, we solve subproblem (SP °) which is of the form of (gP) where H = H and no outer approximation cuts have yet been added. The optimal objective value of subproblem (SP °) yields the initial lower bound LB ° on the optimal objective value of problem (P). If the solution to (SP °) is integer and satisfies the relaxed nonlinear constraints in the set G, then it is stored as the incumbent solution x c and evaluated in f ( x ) to yield upper bound UB ° on the optimal objective value of problem (P). Otherwise, upper bound UB ° is initialized to LIB° = + ~ and no x c is available. I f L B ° = UB °, then the algorithm terminates with the globally optimal solution x c and optimal objective value f(x¢). Otherwise, subproblem (SP °) is added to the candidate list of unfathomed subproblems and the algorithm proceeds to iteration 1. At the beginning of a general iteration k, the optimal objective value, solution, and rectangular partition element of every subproblem in the candidate list is available. Also, LB k- l, UB k- 1, x ¢, and the set of s outer approximation cuts generated so far are available. Any subproblem that is infeasible, has an optimal objective value greater than or equal to UB k- t, or has a corresponding rectangle H that is a single point cannot provide a lower cost solution than the
176
K.M. Bretthauer, M.I. COtE/ European Journal of Operational Research 96 (1996) 167-179
incumbent and, therefore, is eliminated from the candidate list (or, fathomed). If the candidate list is now empty, then x c is the globally optimal solution with optimal objective value f ( x ~). Otherwise, from the subproblems remaining in the candidate list after fathoming, one is selected for further consideration. Let this subproblem be denoted as (-S--P-)with solution 2 and rectangular partition element H. If 2 is fractional or if ~ ~ G and integer, then is divided into two subrectangles H l and H 2 as previously described. This yields the following two new subproblems.
(SP v)
v = 1,2
Min
c°(x)
s.t.
r , ( x) <<.O, x~POH
w = l . . . . . s,
updated to be the minimum of UB k- ~ and any upper bounds obtained in iteration k. The incumbent solution x c is updated to be the solution providing UB k. If LB k = UB k, then the algorithm terminates with the globally optimal solution x ¢. Otherwise, the algorithm proceeds to iteration k -I- 1.
Proposition 3. The branch and bound/outer approximation algorithm described above provides a globally optimal solution to problem (P) in a finite number of iterations. Proof. The algorithm belongs to the class of prototype branch and bound methods for discrete global optimization developed in [3]. Therefore, by Proposition 1 in [3], the result is true. []
°.
In subproblem (SP°), c°(x) is the convex envelope of f ( x ) over H o, v = 1, 2. These two subproblems are solved and added to the candidate list. The other possibility is that 2 ~ G and is integer. Then, instead of subdividing H we construct and add a new outer approximation cutting plane to subproblem ~ that cuts off the point 2 from G. The one new subproblem has the following form:
5.6. Extending the algorithm to solve model (C2) Extending the algorithm to solve model (C2) requires an additional branching mechanism to handle the nonconvex constraints x , yi, = 0, for all (i, t). This branching mechanism is as previously described in the branch and bound algorithm for the continuous capacity, convex cost version of model (C2).
(Sp cut) Min
~(x)
s.t.
rw(x)<~O,
6. Computational testing for continuous capacity, convex cost problems w = l . . . . . s,
x~PAH. In subproblem (SPOUt), s has been updated to s = s + 1 and rw(x) ~< 0, w = 1. . . . . s, includes the set of previously constructed cuts as well as the new cut that separates 2 from G. Subproblem (SP cut) is solved and added to the candidate list. After constructing and solving the new subproblem(s), lower bound LB k is set equal to the minimum of all the optimal objective values of the subproblems currently in the candidate list. All solutions obtained in iteration k that are integer and satisfy the relaxed nonlinear constraints are evaluated in f ( x ) to obtain upper bounds on the optimal objective value of problem (P). Upper bound U B k is
In this section, we report computational results for continuous capacity, convex cost problems solved with Smith and Lasdon's [22] LSGRG code on an IBM 3090-600E computer. As a source of test problems, we used the 13 station, 10 product semiconductor production facility described in [7,8] (see [7,8] for the product routing sequences, the service and interarrival time distribution data, the average dollar work-in-process values per job, and the capacity cost function data). We used a l0 period planning horizon and solved the capacity planning models for a variety of demand patterns, work-in-process levels, and capacity budgets. We solved the continuous capacity, convex cost version of model (C1) for the following two scenar-
K.M. Bretthauer, M.I. C~t~ / European Journal of Operational Research 96 (1996) 167-179
ios: 1) demand is constant from period to period and the firm wishes to gradually decrease work-in-process levels over time; and 2) demand is increasing over time and the finn wishes to keep work-in-process constant from period to period. Average solution time for the scenario 1) problems was 203.8 CPU seconds and average solution time for the scenario 2) problems was 251.7 CPU seconds. We also solved the continuous capacity, convex cost version of model (L1) for a series of different capacity budgets for the case where demand is steady from period to period. Solution time for these problems averaged 73.2 CPU seconds. These results indicate that when capacity can be treated as a continuous variable and costs are convex, commercially available nonlinear programming software can be used to solve models (C1) and (L1) for realistic sized problems.
7. Computational testing of the discrete capacity, nonconvex cost algorithm Next we report computational results with the discrete capacity, nonconvex cost algorithm for model (C1). We focus our computational testing on a six station system where management is interested only in adjusting capacity at two of the stations. To model a nonconvex cost function exhibiting economies of scale, a concave square root function was used. The algorithm was implemented in FORTRAN 77 on an IBM 3090-600E computer. We solved the problems to within 0.01% of optimality. That is, the algorithm terminates when UB k - LB k ~< 0.0001.
UB k
Table 1 Computational results Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU secs.)
1 2 3 4 5 6 7 8 9 10
6 6 6 6 6 6 6 6 6 6
36 15 85 35 32 42 3 10 32 3
4 3 3 5 2 4 1 2 2 1
0.120 0.058 0.254 0.135 0.098 0.134 0.025 0.043 O.O94 0.024
29.3
2.7
0.099
Average: Table 2 Computational results Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU secs.)
11 12 13 14 15 16 17 18 19 20
18 18 I8 18 18 18 18 18 18 18
62 306 10 300 1295 1791 101 1151 33 23
6 8 2 6 37 13 7 23 3 3
0.37 2.19 0.07 1.88 38.67 27.88 0.74 20.17 0.18 0.13
Average:
507.2
177
10.8
9.23
178
K.M. Bretthauer, M.I. C~t# / European Journal of Operational Research 96 (1996) 167-179
Table 1 provides the number o f variables, the number o f subproblems in the branch and bound tree, the number o f cuts added, and the solution time in CPU seconds for a set o f six variable problems. The six variable problems in Table 1 correspond to a single period capacity planning model. The results in Table 1 indicate that the algorithm is able to solve these single period, discrete capacity, nonconvex cost problems in very little computing time. The number o f subproblems in the branch and bound tree varied from 3 to 85, while the number of cuts required remained small for all problems, never exceeding five. To further test the algorithm, we also solved a set of larger problems. Tables 2 and 3 report results for 18 and 36 variable problems. The 18 and 36 variable
problems correspond to 3 and 6 periods in the planning horizon. The results in Tables 2 and 3 suggest that the algorithm can solve moderate sized discrete capacity, nonconvex cost problems in reasonable computing time. Note that one o f the 36 variable problems required 1084.3 CPU seconds to solve, much longer than the other 36 variable problems. This behavior is typical for integer programming because, in the worst case, branch and bound can require complete enumeration of all feasible discrete solutions. However, an advantage of the algorithm is that even if it is terminated before reaching the specified optimality tolerance, it provides a feasible solution within a known percentage of optimality. W e were also interested in determining how the
Table 3 Computational results Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU sees.)
21 22 23 24 25 26 27 28 29 30
36 36 36 36 36 36 36 36 36 36
7283 860 127 245 26369 421 227 87 134 192
11 12 13 7 23 15 7 7 8 8
202.3 12.1 1.6 2.7 1084.3 6.4 2.6 0.9 1.4 1.9
Average:
3 594.5
11.1
131.6
Table 4 Computational results: optimality tolerance = 0.01 Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU sees.)
21 22 23 24 25 26 27 28 29 30
36 36 36 36 36 36 36 36 36 36
432 15 15 15 52 9 14 12 9 8
4 5 5 5 4 3 4 4 3 2
5.77 0.18 0.18 0.18 0.54 0.11 0.16 0.15 0.11 0.09
3.9
0.75
Average:
58.1
K.M. Bretthauer, M.J. C~t~ / European Journal of Operational Research 96 (1996) 167-179
value of the optimality tolerance affects solution time. To address this issue, we solved the 36 variable problems from Table 3 to within 1% of optimality, instead of 0.01% of optimality as in Table 3. The results in Table 4 indicate, as expected, that increasing the value of the optimality tolerance can provide large reductions in solution time. Average solution time decreased from 131.6 CPU seconds for 0.01% of optimality to 0.75 CPU seconds for 1% of optimality.
[6]
[7]
[8]
[9]
8. Conclusions [10]
We have presented models and nonlinear programming methods for multiperiod capacity planning in manufacturing systems that can be modeled as an open network of queues. The optimization framework guides the timing and size of capacity changes and allows performance measures such as capacity costs, work-in-process, and product lead times to be controlled. Solution methods were presented for continuous and discrete capacity options and convex and nonconvex cost functions. Also, note that the capacity planning models presented here can be applied to other queueing network applications, such as computer and communication networks [14] and health care systems [15]. A useful extension to the work presented in this paper would be to modify the models and solution techniques to also treat demand (the arrival rate) as a decision variable that can be controlled by the Marketing function, as in [13].
[11]
[12]
[13]
[14]
[15]
[16] [17] [18]
References [1] Androulakis, I.P., Maranas, C.D., and Floudas, C.A., " a B B : A global optimization method for general constrained nonconvex problems", Journal of Global Optimization 7 (1995) 337-363. [2] Benson, H.P., "Separable concave minimization via partial outer approximation and branch and bound", Operations Research Letters 9 (1990) 389-394. [3] Benson, H.P., Erenguc, S.S., and Horst, R., " A note on adapting methods for continuous global optimization to the discrete case", Annals of Operations Research 25 (1990) 243-252. [4] Berry, W.L., Schmitt, T.G., and Volimann, T.E., "Capacity planning techniques for manufacturing control systems: Information requirements and operational features", Journal of Operations Management 3 (1982) 13-25. [5] Bitran, G.R., and Sarkar, D., "Throughput analysis in manu-
[19]
[20]
[21]
[22]
[23] [24]
179
facturing networks", European Journal of Operational Research 74 (1994) 448-465. Bitran, G.R., and Sarkar, D., "Targeting problems in manufacturing queueing networks - An iterative scheme and convergence", European Journal of Operational Research 76 (1994) 501-510. Bitran, G.R., and Tirupati, D., "Tradeoff curves, targeting and balancing in manufacturing queueing networks", Operations Research 37 (1989) 547-564. Bitran, G.R., and Tirupati, D., "Capacity planning in manufacturing networks with discrete options", Annals of Operations Research 17 (1989) 119-135. Boxma, O.J., Rinnooy Kan, A.H.G., and Van Vliet, M., "Machine allocation problems in manufacturing networks", European Journal of Operational Research 45 (1990) 47-54. Bretthauer, K.M., "Capacity planning in manufacturing and computer networks", European Journal of Operational Research 91 (1996) 386-394. Bretthauer, K.M., "Capacity planning in networks of queues with manufacturing applications", Mathematical and Computer Modelling 21 (1995) 35-46. Bretthauer, K.M., Cabot, A.V., and Venkataramanan, M.A., "An algorithm and new penalties for concave integer minimization over a polyhedron", Naval Research Logistics 41 (1994) 435-454. Buss, A.H., Lawrence, S.R., and Kropp, D.H., "Volume and capacity interaction in facility design", lIE Transactions 26 (1994) 36-49. Gerla, M., and Kleinrock, L., "On the topological design of distributed computer networks", IEEE Transactions on Communications 25 (1977) 48-60. Glenn, J.K., and Roberts, S.D., "The relationship between resource and utilization factors in an outpatient care system", AIIE Transactions 5 (1973) 24-32. Hall, R.W., Queueing Methods for Services and Manufacturing, Prentice-Hall, Englewood Cliffs, NJ, 1991. Horst, R., and Pardalos, P.M. (eds.), Handbook of Global Optimization, Kluwer Academic Publishers, Dordrecht, 1995. Horst, R., and Tuy, H., Global Optimization: Deterministic Approaches, Springer-Verlag, Berlin, 1990. Karmarkar, U.S., Kekre, S., and Kekre, S., "Capacity analysis of a manufacturing cell", Journal of Manufacturing Systems 6 (1987) 165-175. Luss, H., "Operations Research and capacity expansion problems: A survey", Operations Research 30 (1982) 907947. Pardalos. P.M., and Rosen, J.B., Constrained Global Optimization: Algorithms and Applications, Lecture Notes in Computer Science, Vol. 268, Springer-Verlag, Berlin, 1987. Smith, S., and Lasdon, L., "Solving large sparse nonlinear programs using GRG", ORSA Journal on Computing 4 (1992) 2-15. Solberg, J.J., "Capacity planning with a stochastic workflow model", AIIE Transactions 13 (1981) 116-122. Van Vliet, M., and Rinnooy Kan, A.H.G, "Machine allocation algorithms for job shop manufacturing", Journal of Intelligent Manufacturing 2 (1991) 83-94.
ELSEVIER
European Journal of Operational Research 96 (1996) 167-179
Theory and Methodology
Nonlinear programming for multiperiod capacity planning in a manufacturing system Kurt M. Bretthauer *, Murray J. Ctt6 Department of Business Analysis and Research, Texas A & M University, College Station, TX 77 843, USA
Received July 1995
Abstract In this paper, we present nonlinear programming methods for capacity planning in a manufacturing system that consists of a set of machines or work stations producing multiple products. We model the facility as an open network of queues where capacity at each work station in the system may be changed in each of a finite number of time periods. To determine the timing and size of capacity changes, we present two nonlinear programming models and methods for solving the resulting problems. One model involves minimizing total capacity costs such that plant congestion is controlled via upper limits on work-in-process. The other model involves minimizing a weighted sum of product lead times subject to budget constraints on capacity costs. We present solution methods for continuous and discrete capacity options and convex and nonconvex (e.g., economies of scale) capacity cost functions. We use branch and bound and outer approximation techniques to determine globally optimal solutions to the nonconvex problems. Computational testing of the algorithms is reported. Keywords: Capacity planning; Nonlinear programming; Branch and bound; Nonconvex optimization
I. Introduction In this paper, we present nonlinear programming methods for capacity planning in a manufacturing system. Consider a facility consisting of a set of machines or work stations producing multiple products where the sequence of stations visited depends on the processing required by a particular product type. We are interested in determining the timing and size of capacity changes at each work station over a finite number of time periods. The capacity planning decision must take into account the product demand forecasts, existing capacity levels, the costs of changing capacity, the impact of capacity changes on the * Corresponding author.
performance of the manufacturing system, and the firm's manufacturing strategy. During periods of increasing demand it may become necessary to add capacity to some or all of the work stations to satisfy the demand and maintain stable operating conditions in the plant. If demand is expected to decrease for a long period of time, then it may be cost effective to reduce capacity. We model the manufacturing system as an open network of queues and present an optimization framework for capacity planning over a multiperiod planning horizon. The decision variables are the service rates (capacity) at each work station in each time period. Capacity can be controlled at the work stations via the number of machines, modernizing or updating equipment, additional maintenance, the
0377-2217/96/$15.00 @) 1996 Elsevier Science B.V. All fights reserved P// S0377-2217(96)00061-6
168
K.M. Bretthauer, M.I. C~t~ / European Journal of Operational Research 96 (1996) 167-179
number of workers, the number of shifts, the use of overtime, providing additional worker training, etc. We present two nonlinear programming models for capacity planning. Solution methods are developed for continuous and discrete capacity options and for convex and nonconvex (e.g., economies of scale) capacity cost functions. By incorporating stochastic queueing behavior at the work stations within the optimization framework, the models not only identify capacity levels such that forecasted demand in each period will be satisfied, but they also control performance measures such as product lead times, work-in-process, and capacity costs. Specifically, one model involves the minimization of capacity expansion costs minus capacity reduction savings (or plus capacity reduction costs, if appropriate) subject to upper limits on work-in-process. The other model involves the minimization of a weighted sum of product lead times subject to budget constraints on capacity costs. Also note that the capacity planning models presented here can be applied to other queueing network applications, such as computer/communication systems [14] and health care systems [15]. Bitran and Tirupati [7,8] consider a similar problem involving a one time, or single period, capacity change decision in a manufacturing queueing network. They address only the size of capacity changes at the work stations; they do not consider the timing aspect of the decision. Bitran and Tirupati present a model that minimizes capacity costs subject to a single work-in-process constraint. In [7], they treat capacity as a continuous variable and present a heuristic solution technique for expanding capacity. They briefly discuss how their method can be modified to also handle capacity reductions. In [8], they present a heuristic for the capacity expansion problem when capacity is available in discrete options. See Bitran and Sarkar [5,6], Boxma, Rinnooy Kan and Van Vliet [9], Bretthauer [10,11], Buss, Lawrence and Kropp [13], and Van Vliet and Rinnooy Kan [24] for related single period models. Karmarkar, Kekre and Kekre [19] and Solberg [23] also model the manufacturing system as a queueing network and discuss capacity issues, but do not present an optimization model for making capacity change decisions. Although it is relatively straightforward to extend
the single period models of Bitran and Tirupati [7,8] to the multiperiod case, solving the resulting nonlinear optimization problems can be quite difficult. As will be discussed in later sections of the paper, these nonlinear optimization models often will be nonconvex. Therefore, standard nonlinear programming methods may converge to only a locally optimal solution, not necessarily a globally optimal solution. Solving the capacity planning models may be further complicated by situations where capacity is available only in discrete options. In this paper, we present methods for determining globally optimal solutions to these types of nonconvex multiperiod problems, including the case of discrete capacity options. Bitran and Tirupati [7,8] present heuristics for the single period models. See [1,2,3,12,17,18,21] for a detailed discussion of global optimization algorithms and the related literature. Much of the capacity planning literature considers the long term problem, such as when and where to build additional facilities [20]. Here, we focus on a medium term planning horizon where capacity may be changed at the individual work stations in an existing facility. Other medium term capacity planning research considers methods for estimating how much capacity will be necessary at each work station to satisfy a given master production schedule [4]. However, these approaches do not present a decision model for adjusting capacity levels if a mismatch between capacity required and capacity available is identified. Also, these approaches do not take into account the relationship between capacity levels and system performance measures such as work-in-process and product lead times. Our models and solution methods provide recommendations for the size and timing of capacity changes that account for capacity costs and the impact of capacity levels on the performance of the system. The rest of the paper is organized as follows. In Section 2, we present two versions of a capacity planning model that involves minimizing capacity costs subject to work-in-process constraints. In Section 3, we present capacity planning models that involve minimizing product lead times subject to budget constraints on capacity costs. In Section 4, we discuss solution techniques for the continuous capacity, convex cost version of the models. Section 5 provides a branch and bound/outer approximation
K.M. Bretthauer, M.I. C~t~ / European Journal o f Operational Research 96 (1996) 167-179
algorithm for solving the discrete capacity, nonconvex cost version of the problem. In Sections 6 and 7, we provide computational results for continuous capacity, convex cost problems and discrete capacity, nonconvex cost problems. Section 8 contains some concluding remarks.
2. The minimum cost capacity planning model We consider two versions of a minimum cost capacity planning model that minimizes capacity costs subject to work-in-process constraints. Which version is appropriate depends on the manufacturing system under consideration, as discussed in the following subsections.
2.1. Simultaneous capacity increases and decreases permitted Problem (C1) permits the addition of new capacity and the elimination of old capacity in the same time period at a given work station. This option may be attractive if existing capacity at a station can be replaced with new cheaper or more efficient capacity (possibly a new technology). (C1) P
Min
P
~',fit(xit) - ~ iffil t=l
s.t.
i=1
(1)
Y'~gi,(Yit) t=l
~vitLit(lxit)<-..WlPt,
t~T,
(2)
i=1
I'Zit = [£i,t- 1 ~- Xit -- Yit,
i ~ N,
t~
T,
(3) (4)
O<~xlt<~cit,
i ~ N, t ~ T, i~N, t~T,
0 <~Yit <~ dit,
i ~ N,
(6)
air <~ ~it <~ bit,
t ~ T.
(5)
Subscripts i t N T
= = = =
Work station subscript, i = 1. . . . . n. Time period subscript, t = 1. . . . . p. { 1, 2 . . . . . n} = Index set for work stations. {1, 2 . . . . . p} = Index set for time periods.
Decision variables l,Zit Capacity (service rate) at work station i in ~--
period t.
169
xit = Increase in capacity at work station i available at the beginning of period t. Yit = Decrease in capacity at work station i at the beginning of period t.
Functions fit(x~t) = Cost of increasing capacity at work station i in period t by xit. ggt(Yit) = Savings from decreasing capacity at work station i in period t by Yir Lit(I,l~it) = M e a n number of jobs at work station i in period t as a function of /zit (discussed further below).
Parameters v,.t = Average dollar work-in-process value per job at work station i in period t. WlP t = An upper limit on total work-in-process allowed in the facility in period t. /xi0 = Initial capacity at work station i. Ais = Arrival rate (demand) at work station i in period t (assume A,.t > 0 for all (i, t)). ai,, bit = Lower and upper bounds on capacity at work station i in period t ( a , > A,). cit = Upper bound on capacity increase at work station i in period t. dit = Upper bound on capacity decrease at work station i in period t. ca,. t = Squared coefficient of variation (variance of a random variable divided by the square of its mean) of the interarrival times at work station i in period t. cs, = Squared coefficient of variation of service times at work station i in period t. Note that for all models presented in this paper, the capacity reduction savings function git(Yi,) can easily be modified to handle settings where capacity reductions result in costs instead of savings. Problem (C1) requires the ability to calculate system performance measures in a network of queues (in particular, the number of jobs at each station). To make the model as flexible as possible, we consider general interarrival time distributions and general service time distributions at each single server station. We use the approximate decomposition approach to analyze this class of queueing networks
K.M. Bretthauer, M.J. C~td/ European Journal of OperationalResearch 96 (1996) 167-179
170
(see [7,8] and the references therein). The parameters c a . and cs/t are assumed to be constant with respect to capacity. See Bitran and Tirupati [7,8] for justification of this assumption and a discussion of the system of equations to be solved to determine the parameters A., ca~t, and cs/t for all (i, t). We also use a constant arrival rate Ait within a particular period t. This assumption is appropriate when changes in demand occur at discrete points in time, such as when a contract is received from a new customer. The queueing results can also be expected to provide good approximations when demand changes continuously at a slow rate (see Hall [16, pp.178-180] for further discussion of this issue). Constraint set (2) imposes an upper limit on total work-in-process in the facility in each time period. Therefore, we must be able to calculate the average number of jobs Li,(jzit) at each work station as a function of the station capacity. Bitran and Tirupati [7] use the following approximation to calculate Lit(/~,) and show that it is a convex function of/z~:
t.,,( it.) = X,t/u. +
( c a i t + csit ) A~t
21.~it( I~it- Air )
h( X., #i,. ca.. cs.),
where
h('Ait, ~£it' ca/t,
CSit)
[ - 2(1
=, exp[
-
cait)(
Problem (C2): Same as problem (C1) except the following constraints are added:
xitYit=O,
iEN,
(7)
tET.
Constraint set (7) eliminates the possibility of increasing and decreasing capacity at a given work station in the same time period. However, adding these constraints can make problem (C2) much more difficult to solve than (C1), as is discussed later in the paper.
3. The minimum lead time capacity planning model The next capacity planning model that we consider involves minimizing a weighted sum of product lead times rather than minimizing capacity costs. To control capacity costs, we include constraints that limit the amount of money that can be spent on capacity. Minimizing lead times is attractive when management wants to focus on increasing manufacturing flexibility and improving customer service.
I~it- A/t )
3a~(-~a. + cs/~ 1
a given work station are not permitted. For example, a company may not wish to, or be able to, replace experienced workers with new inexperienced workers, even though this strategy may reduce salary expenses. Problem (C2) is as follows:
if c a , < 1,
3.1. Simultaneous capacity increases and decreases permitted
if c a i t > 1,
To handle multiple servers at a station, the expression for Lit given in Bitran and Tirupati [8, pp.134135] can be used in place of the single server equation given above. In this case,the number of servers would be the decision variable instead of the service rate, Constraint (3) is a set of 'capacity flow balance' equations. Constraint sets (4), (5), and (6) set lower and upper bounds on the capacity variables.
The minimum lead time model where increases and decreases in capacity are allowed at a given station in the same time period is formulated as follows. (L1)
(8)
Min s.t.
~L,(x.)i=1
2.2. Simultaneous capacity increases and decreases not p e r m i t t e d
Next we consider environments where increases and decreases in capacity in the same time period at
g.(yi,)~B,,
~it = l~ia - I + x . - )'it,
air <~ I,~it <~bit,
ter,
i=l
i ~ N,
i ~ N, t ~ T,
O <~ xit <~ cit,
i ~ N,
t ~ T,
O <~ Yit <~ dit,
i ~ N,
t E T.
t E T,
(9)
K.M. Bretthauer, M.J. C~M/ European Journal of Operational Research 96 (1996) 167-179
The objective function W(/~) represents a weighted sum of product lead times as a function of the vector/x of work station capacities. Lead times are calculated from Little's equation using the expression for Lit(i~it) previously described. A particular product's total lead time through the system is calculated by summing the time spent at each work station on its route. As in Karmarkar, Kekre and Kekre [19], we weight each product's lead time by its demand when calculating W(/.~). The lead times could also be weighted by cost, to allow for low demand but expensive to hold products. Constraint set (9) sets a capacity budget of B t dollars in period t, t = l . . . . . p. The remaining constraints are the same as in problem (C 1).
3.2. Simultaneous capacity increases and decreases not permitted Next we modify the minimum lead time model (L1) to handle environments where capacity cannot be increased and decreased at a given station in the same time period. Problem (L2) is as follows: Problem (L2): Same as problem (L1) except the following constraints are added:
)fit Yit =
0,
i ~ N,
t ~ T.
4. Continuous capacity, convex cost solution methods Here we address settings where capacity can be treated as a continuous variable and make the following assumption concerning the capacity cost and savings functions. Assumption 1. The cost function fit(x.) is convex with respect to xit for all (i, t) and the savings function g.(Y~t) is concave with respect to Yit for all (i, t). This assumption is appropriate for f , t ( x , ) when the cheapest capacity expansion alternatives are selected first and each subsequent addition may become more expensive on a per unit basis. The as-
171
sumption is appropriate for git(Yit) when the capacity reductions yielding the greatest savings are selected first and each subsequent reduction may yield less savings on a per unit basis. If Assumption 1 is not satisfied (e.g., a concave capacity cost function to handle economies of scale), then the solution obtained may be only locally optimal instead of globally optimal. This issue is discussed later in the paper. When Assumption 1 is satisfied and capacity is a continuous variable, models (CI) and (L1) are convex optimization problems and, therefore, can be solved with a commercially available nonlinear programming package.
4.1. Solving a special case of model (C2) Consider model (C2) for the case where capacity is a continuous variable and Assumption 1 is satisfied. Because of the constraints xitYit = 0 for all (i, t), problem (C2) is a nonconvex optimization problem. Therefore, standard convex programming techniques may converge only to a locally optimal solution instead of a globally optimal solution. Fortunately, we are able to show that with some very reasonable additional assumptions on the cost functions, the constraints x . y . = O for all (i. t) will automatically be satisfied and, therefore, can be eliminated from the formulation. Consider the following assumption: Assumption 2. f/t(0)= 0; git(O) = 0; f/,(x/t) is an increasing convex function of xit; git (Yit) is a nondecreasing concave function of yi,; and f / t ( f l ) > g . ( f l ) for all /3 > 0. In addition to f.(xit) being convex and g . ( y . ) being concave, we are also assuming that there are no costs or savings if capacity is not changed, that additional capacity always costs something, that additional capacity reductions always yield a savings or salvage value greater than or equal to zero, and that the cost of obtaining a given amount of new capacity is more than the savings or salvage value from reducing capacity by the same amount at a given station in a given time period.
K.M. Bretthauer, M.I. COtd/ European Journal of Operational Research 96 (1996) 167-179
172
R e m a r k 1. Assumption 2 implies
Then,
fit(xit) - f i t ( x i t - 3/) > g i t ( Y ) - git(0)
z*--~.=[f(x*)-g(y*)]-
for all xi, >1 y > O.
=
t Xit "i,t~R
R e m a r k 2. Assumption 2 implies
for all Yit ~ "Y > O.
Proof. Let Q = {(i, t) [ x:t Yi; = 0}, let
--
fit Xit--Yit i,
R
( E g i t ( Y i : ) -- E git( 0 ) } i,tER i,t~R
L t ( ' Y ) --fit(O) > git( Yit) -- git( Y i t - "Y)
Proposition 1. If Assumption 2 is satisfied, then the optimal solution Ix:t, x:t, Yi7 for all (i, t) to the convex problem (C 1) is a globally optimal solution to the nonconvex problem (C2).
[f( 2) -- g( ~)]
idES
idES
idES idES Because we have assumed R U S ¢ 0, the above expression for z * - ~ along with Remarks 1 and 2 imply z * - z > 0. The inequality z * - ~ > 0 contradicts the fact that xi,, Yit, Ixit for all (i, t) is the optimal solution to problem (C1). Therefore, we must have R U S = I ~ . []
and let
Proposition 1 implies that we can solve the convex problem (C1) instead of the much more difficult nonconvex problem (C2) when Assumption 2 is satisfied.
S = {(i, t) I x:t Yi7 ~ 0 and xit < Yit }.
4.2. Solving a special case of model (L2)
R = {( i, t) I get Yi7 ~ 0 and x:t >1YeT },
We must show that x:t Yi~ = 0 for all (i, t) to verify that Ixit, xit, Yit for all (i, t) is the globally optimal solution to problem (C2), i.e., that R U S = 0. The proof will be by contradiction. Assume, to the contrary, that R t,J S 4: 0. Consider the solution Txit = Ix:t for all (i, t), 2it = x:t and S:it = Yi~ for all (i, t) ~ Q, 2it = Xit -- Yi~ and Pit ~" 0 for all (i, t) ~ R, and 2i, = 0 and Yi, = Yi~ - x:t for all (i, t) E S. Note that the solution ~i,, Ycit, Pit is feasible to problems (C1) and (C2). Let f ( 2 ) = Y'-fit(xit) i,t
and g ( Y ) = E g i , ( ~Yi,). i,t Let f ( x* ) and g( y * ) be defined similarly. Also, let =f(2)
-- g ( y )
Here we consider model (L2) for the case where capacity is continuous and Assumption 2 is satisfied. Let IX¼, x:t, and Yi~ for all (i, t) now denote the optimal solution to the convex problem (L1) and let Q = {(i, t) I x:t Yi; = 0}, g = {(i, t) I x:t Yi7 ~ 0 and xit >~Yi, },
and S = {(i, t) lx:t Yi7 ~ 0 and xi, < Yi~ }. Similar to Proposition 1, consider the following solution: T-~it~-~lbi; for all (i, t), 2it = Xit and ~:it = Yi, for all (i, t ) ~ Q, 2it = x:t -Yi~ and Pit = 0 for all ( i , t ) E R , and 2 , = 0 and P i t = Y i ~ - x : t for all
(i, t ) E S . Proposition 2. If Assumption 2 is satisfied, then Ixit, xit, and Pit for all (i, t) is a globally optimal solution to the nonconvex problem (L2).
and
Proof. If R U S = 0, then Pit ~-~£it' 2it ~-Xit' and Pit = Yi~ for all (i, t) is clearly a globally optimal
Z* = f ( x * ) - - g ( y * ) .
solution to problem (L2). If R tJ S 4= ¢ and Assump-
K.M. Bretthauer, M.L C~t~/ European Journal of Operational Research 96 (1996) 167-179
tion 2 is satisfied, then, as in the proof of Proposition 1, we can show that [f(x')-g(y*)]-
[ f ( 3:) - g( ~)] > 0 .
This inequality implies that the solution izit, xit, and ~i, for all (i, t) is feasible to problem (L2). We also know W ( ~ ) = W(/z* ) because ~i, =/zit for all (i, t). Therefore, izit, xit, and 3it for all (i, t) is a globally optimal solution to problem (L2). [] Proposition 2 shows that if Assumption 2 is true, then a globally optimal solution to the nonconvex problem (L2) is easily constructed from the solution to the convex problem (L1). 4.3. A branch and bound algorithm f o r models (C2) and (L2)
Next we discuss a branch and bound method for finding a globally optimal solution to the nonconvex problems (C2) and (L2) for the case where capacity is continuous, Assumption 1 is satisfied, and Assumption 2 is not satisfied. The algorithm relaxes the constraints xitYit = 0 for all (i, t) and solves a convex nonlinear subproblem of the form of problem (C1) (or (L1)) at each node of the search tree. Branching is performed on a selected subproblem when its solution does not satisfy constraint set (7), i.e., when xitYit 4:0 for some (i, t). Consider such a subproblem where xitYit ~ 0. A new left branch subproblem is constructed by adding the upper bound xit ~< 0 to the selected subproblem. A new right branch subproblem is constructed by adding Yit <~0 to the selected subproblem. These two new subproblems are solved and added to the candidate list of unfathomed subproblems. Otherwise, standard branch and bound logic and fathoming rules are used. The algorithm terminates after a finite number of iterations with the globally optimal solution to problem (C2) or (L2).
5. Discrete capacity, nonconvex cost solution method In this section, we present a method for solving models (C1) and (C2) for the case of a nonconvex capacity cost (or savings) objective function and
173
discrete capacity options. Examples of discrete capacity variables include the number of machines, shifts, or workers. If costs exhibit economies of scale, then a nonconvex cost function is appropriate. The solution method combines branch and bound with outer approximation techniques. In Bretthauer [11], we presented a version of this algorithm for solving single period problems with a concave objective function, a set of convex nonlinear constraints, and discrete variables (see [12] for the special case where all constraints are linear). Here we extend the approach in [11] to handle a more general nonconvex objective function as well as the nonconvex constraints xitYit = 0 for all (i, t). The algorithm presented here also takes advantage of the fact that some of the constraints in (C 1) and (C2) are linear. In addition, we report computational results with the algorithm, which was not done in [11]. See Benson [2] for further background on branch and bound/outer approximation methods for concave cost, continuous variable problems. First we will present the algorithm for solving problem (C1), which does not include the constraints xit Yit = 0 for all (i, t). Then we will show how to modify the algorithm to handle these nonconvex constraints for problem (C2). To simplify the notation and presentation, we describe the algorithm for a general optimization problem with the structure of model (C 1). Consider the following problem:
(P) n
Min
E f/(Xi) i~l
s.t.
x~GNPNH,
x discrete. In problem (P), x = ( x I . . . . . Xn) is an n-vector of decision variables, f,.(Xi) , for i = 1. . . . . n, is a realvalued cost function (potentially nonconvex), G is a convex set defined by m convex nonlinear constraints gj(x)~< 0, for j = 1. . . . . m, P is a polyhedral set defined by a set of linear constraints, and H-
{ x ~ Rn I / i < x , < u i, i = 1 . . . . . n}
K.M. Bretthauer, M.I. COt~/ European Journal of Operational Research 96 (1996) 167-179
174
is a set of finite and integer valued lower and upper bounds on x. Also, let f(x)
= ~ f~(xi)
mate f ( x ) , leading to tighter lower bounds. The outer approximation cuts help force the subproblem solutions to satisfy the relaxed nonlinear constraint set G.
i=1
and, for notational convenience, let g ( x ) = m a x { g j ( x ) l j = 1 . . . . . m}. Note that g ( x ) is a convex function. Then G can be written as
G= { x ~ R"lg(x) <~0}. We assume integer decision variables. The following discussion remains valid for other discrete sets of variable values. To write the general problem (P) in terms of model (C1), the vector x would include lZit, x u, and Yit for all (i, t), G would include the nonlinear work-in-process constraints (2), P would include the linear capacity flow balance constraints (3), and H would include the bounds on the variables (4), (5), and (6). As in standard branch and bound methods, the algorithm generates a search tree that requires the solution of one subproblem at each of its nodes. To allow the calculation of lower bounds on the optimal objective value of the original problem (P), the subproblems are constructed using three types of 'relaxations': 1) the discrete variable requirements are relaxed; 2) a convex underestimating function to the nonconvex objective is used; and 3) the convex nonlinear constraints in the set G are replaced by a series of linear outer approximation cutting planes. These three relaxations result in convex programming subproblems with all linear constraints, making their solution much simpler than if the nonconvex objective was maintained. We are particularly interested in cost functions that exhibit economies of scale, i.e., concave cost functions. In this case, the underestimating function is linear, resulting in linear programming subproblems. In every iteration of the algorithm, either a subrectangle of H is divided into two smaller subrectangles (rectangular partitioning, or branching), yielding two new subproblems, or an outer approximation cutting plane is added, yielding one new subproblem. Rectangular partitioning helps force the variables to take on discrete values and also causes the convex underestimating function to more accurately approxi-
5.1. Rectangular partitioning Let H={X~nlii
<~Xi <~Ui, i =
1..... n}
denote any rectangular subset of the initial rectangle H. Also, let 2 denote a point in H. First consider the case where 2 has a fractional component 2q. We divide /~ into the two subrectangles,
-' = (x, filx,,
t2,,l)
and
,-,+ = (x +,+ x,, >- l+,,l), where [. ] and [. ] represent the integer round-down and round-up operations, respectively. Now consider the case where all elements of 2 are integer valued. We select acomponent 2q with lq < fiq. In this case we divide H into
H1= ( x E H I Xq<<.2q-1)
H2=(xEHIXq>/2q}. (Note that if 2 o = lq, then instead we divide H into
nl={xE-Hlxq<~2q}
and n 2 = { x ~ H I X q ~ 2 q
+ 1}.) If H is a single point, then no partitioning will be performed and the subproblem is fathomed.
5.2. Convex envelopes As previously mentioned, in the subproblems we replace the nonconvex objective f ( x ) with a convex underestimating function. In particular, we replace f ( x ) with its convex envelope over a given rectangular partition element. Definition 1 (Convex envelope [18]). The convex envelope of a function f ( x ) taken over a nonempty
K.M. Brenhauer, M.J. C~td / European Journal of Operational Research 96 (1996) 167-179
subset M of the domain of f ( x ) is a function c(x) such that: l) c(x) is a convex function defined over the convex hull of the set M. 2) c(x)<~f(x) for all x ~ m . 3) If h(x) is any convex function defined over the convex hull of M and if h(x) <~f(x) for all x E M, then h(x) ~< c(x) for all x in the convex hull of M. The convex envelope is particularly attractive as an underestimating function to a separable concave function because then it is easy to construct and linear. We will construct convex envelopes ~(x) of f ( x ) over H, which satisfies ~(x)<<,f(x) for all zEH.
5.3. The outer approximation cutting planes As the algorithm proceeds, a series of linear outer approximation cutting planes are added in an attempt to force the subproblem solutions to satisfy the relaxed nonlinear constraints. A new cut is added only when a subproblem is selected whose solution is integer and does not satisfy the nonlinear constraints in the set G. Let ~ denote a subproblem solution that is integer, but not in the set G. In this case, we construct a linear outer approximation cutting plane r(x) ~ 0 that 'cuts off' ~ from G. Consider the line segment connecting the point ~ with some point y in the interior of G. The first step in constructing the cut requires determining with a line search the unique point w on this line segment that intersects the surface of G. Let q denote a subgradient of g(x) at the point w. Then, the affine function r( x) = ( x - w, q) + g( w) yields the outer approximation cut r ( x ) ~< 0 that cuts off the point 2 from G, where ( a , b) denotes the inner product between the vectors a and b.
5.4. The continuous subproblems Given a rectangular subset H of the initial rectangle H, the convex envelope ~(x) of f ( x ) over H, and the set of s outer approximation cuts rw(x) ~< O, w = 1 , . . . , s, constructed up to the current iteration of the algorithm, the convex subproblems will be of the following form:
175
( P3 Min S.t.
rw( X ) <~0,
w=
1.....
s,
x~PnH. Note that when f ( x ) is concave these will be linear programming subproblems. Because of the relaxations used, the optimal objective value of subproblem ~ provides a lower bound on the minimum of f(_x) over the integer values of x in the set G f~ P N H.
5.5. Algorithm description Let LB k denote the lower bound on the optimal objective value of problem (P) at iteration k of the algorithm, let UB k denote the best upper bound found up to iteration k, and let x c denote the incumbent feasible solution to problem (P) providing upper bound LIB k. In iteration 0, we solve subproblem (SP °) which is of the form of (gP) where H = H and no outer approximation cuts have yet been added. The optimal objective value of subproblem (SP °) yields the initial lower bound LB ° on the optimal objective value of problem (P). If the solution to (SP °) is integer and satisfies the relaxed nonlinear constraints in the set G, then it is stored as the incumbent solution x c and evaluated in f ( x ) to yield upper bound UB ° on the optimal objective value of problem (P). Otherwise, upper bound UB ° is initialized to LIB° = + ~ and no x c is available. I f L B ° = UB °, then the algorithm terminates with the globally optimal solution x c and optimal objective value f(x¢). Otherwise, subproblem (SP °) is added to the candidate list of unfathomed subproblems and the algorithm proceeds to iteration 1. At the beginning of a general iteration k, the optimal objective value, solution, and rectangular partition element of every subproblem in the candidate list is available. Also, LB k- l, UB k- 1, x ¢, and the set of s outer approximation cuts generated so far are available. Any subproblem that is infeasible, has an optimal objective value greater than or equal to UB k- t, or has a corresponding rectangle H that is a single point cannot provide a lower cost solution than the
176
K.M. Bretthauer, M.I. COtE/ European Journal of Operational Research 96 (1996) 167-179
incumbent and, therefore, is eliminated from the candidate list (or, fathomed). If the candidate list is now empty, then x c is the globally optimal solution with optimal objective value f ( x ~). Otherwise, from the subproblems remaining in the candidate list after fathoming, one is selected for further consideration. Let this subproblem be denoted as (-S--P-)with solution 2 and rectangular partition element H. If 2 is fractional or if ~ ~ G and integer, then is divided into two subrectangles H l and H 2 as previously described. This yields the following two new subproblems.
(SP v)
v = 1,2
Min
c°(x)
s.t.
r , ( x) <<.O, x~POH
w = l . . . . . s,
updated to be the minimum of UB k- ~ and any upper bounds obtained in iteration k. The incumbent solution x c is updated to be the solution providing UB k. If LB k = UB k, then the algorithm terminates with the globally optimal solution x ¢. Otherwise, the algorithm proceeds to iteration k -I- 1.
Proposition 3. The branch and bound/outer approximation algorithm described above provides a globally optimal solution to problem (P) in a finite number of iterations. Proof. The algorithm belongs to the class of prototype branch and bound methods for discrete global optimization developed in [3]. Therefore, by Proposition 1 in [3], the result is true. []
°.
In subproblem (SP°), c°(x) is the convex envelope of f ( x ) over H o, v = 1, 2. These two subproblems are solved and added to the candidate list. The other possibility is that 2 ~ G and is integer. Then, instead of subdividing H we construct and add a new outer approximation cutting plane to subproblem ~ that cuts off the point 2 from G. The one new subproblem has the following form:
5.6. Extending the algorithm to solve model (C2) Extending the algorithm to solve model (C2) requires an additional branching mechanism to handle the nonconvex constraints x , yi, = 0, for all (i, t). This branching mechanism is as previously described in the branch and bound algorithm for the continuous capacity, convex cost version of model (C2).
(Sp cut) Min
~(x)
s.t.
rw(x)<~O,
6. Computational testing for continuous capacity, convex cost problems w = l . . . . . s,
x~PAH. In subproblem (SPOUt), s has been updated to s = s + 1 and rw(x) ~< 0, w = 1. . . . . s, includes the set of previously constructed cuts as well as the new cut that separates 2 from G. Subproblem (SP cut) is solved and added to the candidate list. After constructing and solving the new subproblem(s), lower bound LB k is set equal to the minimum of all the optimal objective values of the subproblems currently in the candidate list. All solutions obtained in iteration k that are integer and satisfy the relaxed nonlinear constraints are evaluated in f ( x ) to obtain upper bounds on the optimal objective value of problem (P). Upper bound U B k is
In this section, we report computational results for continuous capacity, convex cost problems solved with Smith and Lasdon's [22] LSGRG code on an IBM 3090-600E computer. As a source of test problems, we used the 13 station, 10 product semiconductor production facility described in [7,8] (see [7,8] for the product routing sequences, the service and interarrival time distribution data, the average dollar work-in-process values per job, and the capacity cost function data). We used a l0 period planning horizon and solved the capacity planning models for a variety of demand patterns, work-in-process levels, and capacity budgets. We solved the continuous capacity, convex cost version of model (C1) for the following two scenar-
K.M. Bretthauer, M.I. C~t~ / European Journal of Operational Research 96 (1996) 167-179
ios: 1) demand is constant from period to period and the firm wishes to gradually decrease work-in-process levels over time; and 2) demand is increasing over time and the finn wishes to keep work-in-process constant from period to period. Average solution time for the scenario 1) problems was 203.8 CPU seconds and average solution time for the scenario 2) problems was 251.7 CPU seconds. We also solved the continuous capacity, convex cost version of model (L1) for a series of different capacity budgets for the case where demand is steady from period to period. Solution time for these problems averaged 73.2 CPU seconds. These results indicate that when capacity can be treated as a continuous variable and costs are convex, commercially available nonlinear programming software can be used to solve models (C1) and (L1) for realistic sized problems.
7. Computational testing of the discrete capacity, nonconvex cost algorithm Next we report computational results with the discrete capacity, nonconvex cost algorithm for model (C1). We focus our computational testing on a six station system where management is interested only in adjusting capacity at two of the stations. To model a nonconvex cost function exhibiting economies of scale, a concave square root function was used. The algorithm was implemented in FORTRAN 77 on an IBM 3090-600E computer. We solved the problems to within 0.01% of optimality. That is, the algorithm terminates when UB k - LB k ~< 0.0001.
UB k
Table 1 Computational results Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU secs.)
1 2 3 4 5 6 7 8 9 10
6 6 6 6 6 6 6 6 6 6
36 15 85 35 32 42 3 10 32 3
4 3 3 5 2 4 1 2 2 1
0.120 0.058 0.254 0.135 0.098 0.134 0.025 0.043 O.O94 0.024
29.3
2.7
0.099
Average: Table 2 Computational results Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU secs.)
11 12 13 14 15 16 17 18 19 20
18 18 I8 18 18 18 18 18 18 18
62 306 10 300 1295 1791 101 1151 33 23
6 8 2 6 37 13 7 23 3 3
0.37 2.19 0.07 1.88 38.67 27.88 0.74 20.17 0.18 0.13
Average:
507.2
177
10.8
9.23
178
K.M. Bretthauer, M.I. C~t# / European Journal of Operational Research 96 (1996) 167-179
Table 1 provides the number o f variables, the number o f subproblems in the branch and bound tree, the number o f cuts added, and the solution time in CPU seconds for a set o f six variable problems. The six variable problems in Table 1 correspond to a single period capacity planning model. The results in Table 1 indicate that the algorithm is able to solve these single period, discrete capacity, nonconvex cost problems in very little computing time. The number o f subproblems in the branch and bound tree varied from 3 to 85, while the number of cuts required remained small for all problems, never exceeding five. To further test the algorithm, we also solved a set of larger problems. Tables 2 and 3 report results for 18 and 36 variable problems. The 18 and 36 variable
problems correspond to 3 and 6 periods in the planning horizon. The results in Tables 2 and 3 suggest that the algorithm can solve moderate sized discrete capacity, nonconvex cost problems in reasonable computing time. Note that one o f the 36 variable problems required 1084.3 CPU seconds to solve, much longer than the other 36 variable problems. This behavior is typical for integer programming because, in the worst case, branch and bound can require complete enumeration of all feasible discrete solutions. However, an advantage of the algorithm is that even if it is terminated before reaching the specified optimality tolerance, it provides a feasible solution within a known percentage of optimality. W e were also interested in determining how the
Table 3 Computational results Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU sees.)
21 22 23 24 25 26 27 28 29 30
36 36 36 36 36 36 36 36 36 36
7283 860 127 245 26369 421 227 87 134 192
11 12 13 7 23 15 7 7 8 8
202.3 12.1 1.6 2.7 1084.3 6.4 2.6 0.9 1.4 1.9
Average:
3 594.5
11.1
131.6
Table 4 Computational results: optimality tolerance = 0.01 Problem number
Number of variables
Number of subprbs. in the tree
Number of cuts added
Solu. time (CPU sees.)
21 22 23 24 25 26 27 28 29 30
36 36 36 36 36 36 36 36 36 36
432 15 15 15 52 9 14 12 9 8
4 5 5 5 4 3 4 4 3 2
5.77 0.18 0.18 0.18 0.54 0.11 0.16 0.15 0.11 0.09
3.9
0.75
Average:
58.1
K.M. Bretthauer, M.J. C~t~ / European Journal of Operational Research 96 (1996) 167-179
value of the optimality tolerance affects solution time. To address this issue, we solved the 36 variable problems from Table 3 to within 1% of optimality, instead of 0.01% of optimality as in Table 3. The results in Table 4 indicate, as expected, that increasing the value of the optimality tolerance can provide large reductions in solution time. Average solution time decreased from 131.6 CPU seconds for 0.01% of optimality to 0.75 CPU seconds for 1% of optimality.
[6]
[7]
[8]
[9]
8. Conclusions [10]
We have presented models and nonlinear programming methods for multiperiod capacity planning in manufacturing systems that can be modeled as an open network of queues. The optimization framework guides the timing and size of capacity changes and allows performance measures such as capacity costs, work-in-process, and product lead times to be controlled. Solution methods were presented for continuous and discrete capacity options and convex and nonconvex cost functions. Also, note that the capacity planning models presented here can be applied to other queueing network applications, such as computer and communication networks [14] and health care systems [15]. A useful extension to the work presented in this paper would be to modify the models and solution techniques to also treat demand (the arrival rate) as a decision variable that can be controlled by the Marketing function, as in [13].
[11]
[12]
[13]
[14]
[15]
[16] [17] [18]
References [1] Androulakis, I.P., Maranas, C.D., and Floudas, C.A., " a B B : A global optimization method for general constrained nonconvex problems", Journal of Global Optimization 7 (1995) 337-363. [2] Benson, H.P., "Separable concave minimization via partial outer approximation and branch and bound", Operations Research Letters 9 (1990) 389-394. [3] Benson, H.P., Erenguc, S.S., and Horst, R., " A note on adapting methods for continuous global optimization to the discrete case", Annals of Operations Research 25 (1990) 243-252. [4] Berry, W.L., Schmitt, T.G., and Volimann, T.E., "Capacity planning techniques for manufacturing control systems: Information requirements and operational features", Journal of Operations Management 3 (1982) 13-25. [5] Bitran, G.R., and Sarkar, D., "Throughput analysis in manu-
[19]
[20]
[21]
[22]
[23] [24]
179
facturing networks", European Journal of Operational Research 74 (1994) 448-465. Bitran, G.R., and Sarkar, D., "Targeting problems in manufacturing queueing networks - An iterative scheme and convergence", European Journal of Operational Research 76 (1994) 501-510. Bitran, G.R., and Tirupati, D., "Tradeoff curves, targeting and balancing in manufacturing queueing networks", Operations Research 37 (1989) 547-564. Bitran, G.R., and Tirupati, D., "Capacity planning in manufacturing networks with discrete options", Annals of Operations Research 17 (1989) 119-135. Boxma, O.J., Rinnooy Kan, A.H.G., and Van Vliet, M., "Machine allocation problems in manufacturing networks", European Journal of Operational Research 45 (1990) 47-54. Bretthauer, K.M., "Capacity planning in manufacturing and computer networks", European Journal of Operational Research 91 (1996) 386-394. Bretthauer, K.M., "Capacity planning in networks of queues with manufacturing applications", Mathematical and Computer Modelling 21 (1995) 35-46. Bretthauer, K.M., Cabot, A.V., and Venkataramanan, M.A., "An algorithm and new penalties for concave integer minimization over a polyhedron", Naval Research Logistics 41 (1994) 435-454. Buss, A.H., Lawrence, S.R., and Kropp, D.H., "Volume and capacity interaction in facility design", lIE Transactions 26 (1994) 36-49. Gerla, M., and Kleinrock, L., "On the topological design of distributed computer networks", IEEE Transactions on Communications 25 (1977) 48-60. Glenn, J.K., and Roberts, S.D., "The relationship between resource and utilization factors in an outpatient care system", AIIE Transactions 5 (1973) 24-32. Hall, R.W., Queueing Methods for Services and Manufacturing, Prentice-Hall, Englewood Cliffs, NJ, 1991. Horst, R., and Pardalos, P.M. (eds.), Handbook of Global Optimization, Kluwer Academic Publishers, Dordrecht, 1995. Horst, R., and Tuy, H., Global Optimization: Deterministic Approaches, Springer-Verlag, Berlin, 1990. Karmarkar, U.S., Kekre, S., and Kekre, S., "Capacity analysis of a manufacturing cell", Journal of Manufacturing Systems 6 (1987) 165-175. Luss, H., "Operations Research and capacity expansion problems: A survey", Operations Research 30 (1982) 907947. Pardalos. P.M., and Rosen, J.B., Constrained Global Optimization: Algorithms and Applications, Lecture Notes in Computer Science, Vol. 268, Springer-Verlag, Berlin, 1987. Smith, S., and Lasdon, L., "Solving large sparse nonlinear programs using GRG", ORSA Journal on Computing 4 (1992) 2-15. Solberg, J.J., "Capacity planning with a stochastic workflow model", AIIE Transactions 13 (1981) 116-122. Van Vliet, M., and Rinnooy Kan, A.H.G, "Machine allocation algorithms for job shop manufacturing", Journal of Intelligent Manufacturing 2 (1991) 83-94.