A Self-Organizing Extreme-Point Tabu-Search Algorithm for Fixed Charge Network Problems with Extensions
Richard S. Barr, Fred Glover, Toby Huskinson, Gary Kochenberger
11 A Self-Organizing Extreme-Point Tabu-Search Algorithm for Fixed Charge Network Problems with Extensions
Richard S. Barr , Fred Glover , Toby Huskinson , Gary Kochenberger Abstract
We propose a new self-organizing algorithm for fixed-charge network flow problems based on ghost image (GI) processes as proposed in Glover (1994) and adapted to fixed-charge transportation problems in Glover, Amini and Kochenberger (2005). Our self-organizing GI algorithm iteratively modifies an idealized representation of the problem embodied in a parametric ghost image, enabling all steps to be performed with a primal network flow algorithm operating on the parametric GI. Computational tests are carried out on an extensive set of benchmark problems which includes the previous largest set in the literature, comparing our algorithm to the best methods previously proposed for fixed-charge transportation problems, though our algorithm is not specialized to this class. We also provide comparisons for additional more general fixed-charge network flow problems against Cplex 12.8 to demonstrate that the new self-organizing GI algorithm is effective on large problem instances, finding solutions with statistically equivalent objective values at least 700 times faster. The attractive outcomes produced by the current GI/TS implementation provide a significant advance in our ability to solve fixed-cost network problems efficiently and invites its use for larger instances from a variety of application domains. Introduction: Problem Definition and Background
We define the network fixed charge problem as
NetFC : Minimize x o [FC] = cx + F(x) Subject to: Ax = b U ≥ x ≥
0 where x is the vector given by x = (x j : j ∈ N = {1, …, n}) and the matrix A is a node-arc incidence matrix, so that the equation Ax = b constitutes a classical network representation of the flow equations defining a pure network problem and the variables x j correspond to arcs of the network. The fixed charge function is given by F(x) = ∑ (F j y j : j ∈ N) where each fixed charge EMIS Department, Lyle School of Engineering, Southern Methodist University, Dallas, TX 75275, [email protected] ECEE, College of Engineering and Applied Science, University of Colorado – Boulder, Boulder, CO, 80309 USA [email protected] Computer Science Department, Lyle School of Engineering, Southern Methodist University, Dallas, TX 75275, [email protected] Business School, University of Colorado at Denver, Denver, CO 80217 USA, [email protected] coefficient F j is nonnegative and the y j variables take on binary values that satisfy y j = 1 if x j > 0 and y j = 0 otherwise. F(x) may be equivalently written as
F(x) = ∑ (F j y j : j ∈ N(FC)), where N(FC) = {j ∈ N: F j > 0} and we call N(FC) the set of (effective) fixed charge coefficients. Applications of the problem NetFC arise in many areas, including facility location, network design, logistics and supply chain, and specific problems, such as lot-sizing, course scheduling, and others. Location problems include the uncapacitated and capacitated facility or plant location problems as described in Fernández and Landete (2015) and Eiselt, Marianov, and Bhadury, (2015). In-depth coverage in Daskin (2013) provides in-depth coverage of the area and extensive list of application papers are available in Eiselt, Marianov, and Bhadury (2015). Network design applications arise in telecommunications (Forsgren and Prytz (2006) and Pioro and Medhi (2004)), including related location problems (Fortz (2015), regional wastewater system design (Jarvis,. Rardin, Unger, Moore and Schimpeler (1978)), and electrical smartgrid data network design, including equipment placement, described in Barr, Jones, and Klinkert (2018). NetFC problems also have useful applications in supply chain optimization (Alizadeh (2009)), logistics (Alumur, Kara, and Melo (2015)), vanpool assignment (Kaan and Olinick (2013)), and distribution networks (Mateus and Patrocinio (2006)). In addition, they emerge in multi-level lot-sizing within an MRP (Steinberg and H. Albert Napier (1980)) and scheduling training courses (Glover, Klingman, and Phillips (1992)). See other applications enumerated in Nicholson and Zhang, (2016) In the following, we assume the reader has a basic acquaintance with formulations and solution algorithms for pure networks and is familiar at a rudimentary level with primal simplex algorithms for pure networks. (For references containing useful background information, see for example Ahuja, Magnanti, and Orlin (1993); Bazaraa, Jarvis, and Sherali (2010); and Murty (1992).) The remainder of this paper is organized as follows. Section 2 introduces our Self-Organizing GI Approach FixNetGI for the network fixed-charge problem and gives a pseudocode for its main algorithm, followed by an explanation of the procedure. Section 3 provides the pseudocodes for the routines invoked by the main algorithm, together with a description of their functions. The design for testing our algorithm and the computational results, together with a comparison involving outcomes obtained by applying the Cplex MIP code [15], are presented in Section 4. As shown, the outcomes demonstrate significant advantages for our algorithm both in solution time and solution quality in solving large and challenging
NetFC problems. Finally, Section 5 concludes the paper, highlighting the implications of the computational results and identify directions for future research. The Self-Organizing GI Approach
The general form of the self-organizing ghost image (GI) approach derives from a collection of problem-solving principles detailed in Glover (1994). Our focus on applying the GI framework to fixed-charge network problems builds on the work of Glover, Amini and Kochenberger (2005) that studies an earlier version of the approach applied to the special case of fixed-charge transportation networks.
Within the pure network setting of
NetFC , our method exploits the problem structure by introducing a non-negative penalty vector p = ( p j : j ∈ N ) and an associated penalized cost vector given by c(p) = ( c j + p j : j ∈ N ). The penalty vector p is determined by a self-organizing parameterization to give the following parametric network linear programming relaxation of the fixed-charge problem LP( p ) Minimize x o (p) = c(p)x Subject to: Ax = b U ≥ x ≥
0 The parameterization defining p j occurs by setting p j = F j /v j where v j denotes a quantity that is systematically updated throughout the algorithm. Hence p j allocates the fraction j of the fixed cost F j to the total cost of x j . We apply the convention that a denominator v j close to 0 (smaller than a chosen ε value) translates into setting p j = BigM provided F j > 0, where BigM is a large positive number, and similarly a denominator v j that exceeds BigM translates into setting p j = 0. However, we will in several instances identify the p j values directly without bothering to refer to v j . (Note, if F j = 0 then automatically p j = 0 regardless of the value of v j , since F j = 0 expresses the fact that x j is not a fixed-charge variable. We also interpret the value of x j to be 0 if this value is less than ε .) In the case p = 0 (where all p j = 0), we have c(p) = c, and obtain the simple network linear programming relaxation LP:
Minimize x o = cx Subject to: Ax = b U ≥ x ≥
0 The method sketched in Glover (1994) begins by solving LP, and then solves a succession of problems LP( p ) produced by progressively modifying and updating p j in alternation with applying an improvement method for enhancing the solution to LP( p ), utilizing adaptive memory strategies from tabu search. In our adaptation of the self-organizing GI method to the present context, for simplicity we use the convention of identifying the value of the (nonlinear) fixed-charge objective function x o [FC] for a given trial solution vector x (e.g., x = x' , x" , etc.) as x o (hence, x o ' = x o [FC: x' ], x o " = x o [FC: x" ], etc.) It is important to keep in mind that in such cases x o includes reference to the fixed charge component of the objective function, with the sole exception of explicitly referring to the problem LP. The values U o and U jo defined below are used as proxy bounds for x j that will be introduced to replace the original bound U j in certain calculations of the algorithm. Apart from trial solution vectors, we maintain a locally optimal solution vector x* and an overall (“global”) best solution vector x G , i.e., x oG (= x o [FC: x G ]) is the minimum of the x o * (= x o [FC: x* ]) values. We first give a pseudo-code for the main routines of our Ghost Image Tabu Search (GI/TS) method embodied in our FixNetGI code and then describe the rationale that explains the key steps.
GI/TS Algorithm
The algorithm requires setting the following user input parameters:
Search limits: MaxIter: maximum inside loop iterations per invocation 2.
MaxPass: number of diversification invocations required to terminate algorithm 3.
MaxInsideImprove: number of consecutive non-improving inside loop iterations that will trigger an exit from the inside loop 4.
BadLuck: number of consecutive x*-improvement failures that will trigger a diversification 5.
OutOfLuck: number of consecutive non-improving outside loop iterations that will trigger an exit from the outside loop.
Updating (cid:1874) : Alpha(i), i = 1 to 3: weighting factors, summing to 1, for updating v j values. Weights: Alpha(1) for current x j *, Alpha(2) for current v j value, and Alpha(3) for the historical mean j plus U jo as adjusted by Beta 7. Beta: weight for historical average associated with Alpha(3) and v j update MaxSol: when updating v j , the maximum number of previous x j * values used to calculate mean j for the Alpha(3) term Tabu control: TabuTenure: pivots required before a leaving arc can reenter the network tree (LP basis)
Duplicate solutions
LimMatch: limits the number of times a solution duplication occurs before triggering diversification 11. sLim: number of solutions saved for duplicate-solution checking 12.
ZeroRefresh: number of diversifications performed that will trigger refreshing the duplicate-check solution list with all counts equal 0
The GI/TS algorithm as defined here is supported by several subsidiary procedures to update v , control the descent and tabu phases, perform moves/pivots, check for duplicate solutions, and diversity the search. These components are defined and discussed separately. GI/TS Main Routine Step 0: solve LP and create initial v, p, and locally best solution x* A. Initialize parameters: i. JIter= 0, v iter = MaxIter/4, Pass = 0, LastInsideImprove = 0, Zero(s) = (0, … 0) for s = 1 to sLim (i.e., Zero(s,j) = 0 for j = 1 to n), nMatch = 0, Recover = 0, DoTabu = True, NumSol = 0, GbestIter = 0, NoLuck = 0, BigM = large positive number, AscentTenure = DescentTenure = TabuTenure, ii. Set x oG = BigM B. Solve LP, save the solution as the first locally best solution x* and identify the fixed charge objective function value x o * = x o [FC: x* ]. C. Save the scalar U o as the largest flow value x j , j ∈ N(FC) in the solution to LP. D. Save individual values U jo ( ≤ U j ) = x j as the max flow (so far) for each arc j ∈ N(FC). E. Set v j = U j so that initially p j = F j /U j , mean j = U j for all j ∈ N(FC). F. Step 1: create and solve LP(p) to get first test solution x ′ i. Solve LP( p ) by re-optimization to get x ′ and identify the fixed charge objective function value x o ' = x o [FC: x' ]. NumSol = 1 ii. Update U jo = max(U jo , x j ′ ), for each j ∈ N(FC). iii. If x o ′ < x o * then x o * = x o ' , x* = x' , set Descent = True and perform V_UPDATE iv. Create the n -vector ZeroØ, where ZeroØ(j) = 1 if x j ′ = 0 and F j > 0 (j ∈ N(FC)), else ZeroØ(j) = 0. v. Set First = 1 and Zero(1) = SumZeroØ = ZeroØ vi.
OutsideOK = True While (OutsideOK): (Execute the
Outside loop ) A. Step 2:
Improve the current solution x ′ , move to local optimum x ″ , and then to TS improvement i. Phase I:
Refine x ′ by LP Restriction: a. Set p : p j = BigM if ZeroØ(j) = 1, else p j = 0 b. Solve LP(p) by re-optimization to get x ″ (and x o ″ ) c. If x o ″ < x o * then set Descent = True (Recording of x* = x" will be handled later) d. If JIter < v iter , update U jo = max(U jo , x j ″ ), for each j ∈ N(FC) ii.
Phase II: a. Initialize parameters: a. Set InsideIter = BestIter = TSImprove = DescentImprove = LastInsideImprove = 0, Descent = True, Improve = False, TabuTenure = DescentTenure b. Set Tabu(j) = 0 for each j ∈ N c. Set Aspire = Min(x o ″ , x o *). InsideOK = True b. Repeat while (InsideIter < MaxIter and InsideOK) (Execute the
Inside Loop ) a. ++InsideIter, j*=k*=0 b. For every NB arc j ∈ N: Compute x oj , the change in the objective function x o " (= x o [ FC:x" ]) if x j is pivoted into the basis (and one or more variables x k are driven to their lower or upper bounds to become candidates to leave the basis). Restrict consideration to j ∈ N satisfying Tabu(j) < InsideIter or satisfying the aspiration criterion x oj < Aspire – x o ″ c. Save the best arc j* = arg min( x oj : for j subject to the restriction in b.), and identify a leaving arc k*. (k* = j* if there is a “bound flip” where x j* leaves the basis at its opposite bound) d. Perform DESCEND to carry out the pivot and associated update for the choice of j* and k*. e. If (InsideIter – LastInsideImprove > MaxInsideImprove) then InsideOK = False (
Exit the Inside Loop ) EndWhile (for the
Inside Loop ) c. If ++JIter > MaxIter, OutsideOK = False (Exit the Outside Loop) d. If (Improve) NoLuck = 0 Else a. ++NoLuck b. If NoLuck = OutOfLuck, OutsideOK = False, BREAK (Exit the Outside Loop, to conclude at Step 3)
ElseIf NoLuck = BadLuck, then a. v j = Max(U o – v j , 1) for each j ∈ N(FC) (mini-diversification) b. if x o * < x oG then update x oG = x o * and x G = x* c. x o * = BigM (to assure LP(p) starts over to make a new local optimum x* ) B. Create and solve LP(p) to get new test solution x ′ and Check for Duplications i. Set p j = F j / v j for each j ∈ N(FC) ii. Solve LP(p) by post-optimization to get x ′ and x o ′ iii. Update U o = max(U o , x j ' ) for each j ∈ N(FC) iv. If x o ′ < x o * then a. Update x o * = x o ′ and x* = x ′ b. Perform V_UPDATE v. Create the n-vector ZeroØ, where ZeroØ(j) = 1 if F j > 0 and x j ′ = 0, else ZeroØ(j) = 0 vi. Perform DUPCHECK ( which may include DIVERSIFY ) Endwhile ( Outside Loop ) Conclusion, after exit Outside Loop
A. If x o * < x oG then x oG = x o * and x G = x* and set BestPass = Pass B. STOP
Discussion of the GI/TS Main Routine
In the initialization step, Step 0, the original linear programming relaxation LP is solved, and its solution is saved as the first locally optimal solution x* . Also, to initiate alternative formulas for updating the parameter vector v , the constant U o is initialized to be the largest x j value obtained in solving LP. In addition, the solution value for each variable x j is recorded in U jo . In Step 1, the problem LP( p ) is solved for the first time by re-optimizing the solution obtained in Step 0 for the modified objective function of LP( p ), to obtain an LP optimum solution, x'. The fixed-charge objective function value x o ' ( = x o [FC: x' ]) for x' is calculated and x' replaces the locally best solution x* if x o ' < x o * (= x o [FC: x* ]). We continue to update the values U jo designated to maintain the maximum value attained by x j for the first v iter iterations. To investigate the potential for further improvement to the current solution, x' , in Step 2-Phase I the objective function coefficients of the variables with nonzero and zero values are set to their variable costs c j and c j + Big M , respectively (as a result of setting p j = 0 and p j = BigM in these two cases), resulting in the specified form of LP( p ), which is then solved by post-optimization, yielding x" . The main purpose of setting the cost of variables with the zero values in the trial solution to Big M is to maintain their values at zero during the current post-optimization process, and these variables alternatively could simply be handled by temporarily setting their upper bounds to 0 during this step. Remaining variables that were positive in the solution to the previous LP(p) problem receive their original costs c j so that the solution will be evaluated relative to the original variable costs. Following the calculation of the fixed charge objective function value for the resulting solution x" , the current locally best solution x* is replaced by x" if this new solution turns out to be better. Also, in Phase I the value U j , identifying the maximum value for each x j throughout the first v iter iterations, is updated. Next, the Inside Loop is initiated within Phase II that executes a tentative pivot exploration process, where each nonbasic variable x j , j ∈ N, is considered as a potential entering variable, and the candidates for the leaving variable, x k , are identified, to determine the change x oj in the fixed charge objective function that would result if x j were selected to enter the basis tree. The process is guided by a simple tabu search approach, where attention is restricted to j ∈ N satisfying Tabu(j) < InsideIter or satisfying the aspiration criterion x oj < Aspire – x o ″ , conditions that are irrelevant initially but that become relevant based on updates in the DESCEND routine. The pseudocodes for the DESCEND procedure and other procedures invoked by the main algorithm appear below, followed by a description of the functions of these procedures. At the completion of the tentative pivot explorations within the main algorithm, the variable x j* that yields the greatest reduction in the fixed-charge objective function, is selected for pivoting to bring it into the basis. To further improve the current solution, the process returns to the tentative pivot exploration phase, using the current basis representation. The Inside Loop ends once the current iteration, InsideIter, exceeds the maximum allowed number of iterations, MaxInsideImprove, beyond the last improvement of the locally best solution x* . At the conclusion of the Inside Loop the Outside Loop continues by setting the counter NoLuck to 0 if the Inside Loop had succeeded in improving the locally best solution x* . Otherwise NoLuck is incremented and if NoLuck = OutOfLuck the Outside Loop terminates to record the final global best solution x G at Step 3. Barring this, if NoLuck = BadLuck, a “mini-diversification” step is initiated. Phase II proceeds to generate the current p vector based on the vector v , and then solves LP(p) by post-optimization to obtain x' . If the fixed charge objective function value x o ' = x o [FC: x' ] improves on x o * then x* is updated and the V_UPDATE routine is executed. Finally, the DUPCHECK routine is executed, which may involve executing the DIVERSIFY procedure, to lay the foundation for the next iteration of the Outside Loop. Supporting Procedures
We first give the pseudocode for the supporting procedures used within the main routine, in the order in which they first appear in the main routine and in other supporting procedures, and then explain their functions.
Procedure DESCEND If Descent = True then A. If x oj * < 0 then ( the Decent Phase continues to improve ) i. Perform PIVOTJSTAR ( to pivot in j* and remove k*, … ) ii. Update x o ″ and set Aspire = Min(x o *, x o ″ ). iii. ++DescentImprove B. Else i. Descent = False ( happens the first time that leave Descent Phase ) ii. TabuTenure = AscentTenure iii.
If x o ″ < x o * then a. Improve = True b. BestIter = LastInsideImprove = InsideIter – 1 and BestIterG = JIter c. Update x o * = x o ″ and x* = x ″ d. Perform V_UPDATE iv.
If DoTabu = False, then InsideOK = False, RETURN (EXIT
Inside loop) v. Perform PIVOTJSTAR Else (
Descent = False and we are now in the TS phase ) A. Perform PIVOTJSTAR B. If x oj * < 0 then i. TabuTenure = DescentTenure ii.
If x o ″ < x o * then a. Improve = True b. BestIter = LastInsideImprove = InsideIter, and BestIterG = JIter c. Update x o * = x o ″ and x* = x ″ d. ++TSImprove and ++AllTSImprove e. Aspire = x o * f. Perform V_UPDATE C. Else TabuTenure=AscentTenure Update: Tabu(k*) = InsideIter + TabuTenure Return END DESCEND
Procedure V_UPDATE`
1. ++NumSol 2.
Y = Min(NumSol,MaxSol), X = 1/Y 3.
For each arc j ∈ N(FC): a.
Mean j = (X) x j * + (1-X)Mean j b. UMean = Beta·Mean j + (1-Beta)U o c. v j = Alpha(1)· x j * + Alpha(2)· v j + Alpha(3)·UMean 4. If x o * < x oG then x oG = x o * and x G = x* , GbestIter = JIter END V_UPDATE Procedure PIVOTJSTAR Pivot in j* and remove k* from the tree (or perform a bound flip) yielding a new x ″ and updating x o ″ . As x" is created, set U jo = max(U jo , x j " ) along the basis equivalent path. Return END PIVOTJSTAR 0
Procedure DUPCHECK Set s = First and Match = False (
Match will change to True if some Zero(s) = ZeroØ ) For CheckDup = 1 to sLim and Match = False (
DupCheck loop ) A. If Zero(CheckDup) = ZeroØ then Match = True, Exit this loop B. Else If ++s > sLim then s = 1 If Match = True then ( new change ) A. If ++nMatch > LimMatch then i. sMax = Max(sMax, CheckDup) ( Records how far we had to go to find a match. ) ii. Execute DIVERSIFY iii. nMatch = 0 Else A. If nMatch > 0 then i. ++Recover ii.
MaxRecover = Max(Recover,MaxRecover) iii. nMatch = 0 B. SumZeroØ = SumZeroØ + ZeroØ C. If First > 1 then Last = First – 1, else Last = sLim. D. Replace Zero(Last) by setting Zero(Last) = ZeroØ and set First = Last Return END DUPCHECK PROCEDURE DIVERSIFY If x o * < x oG then x G = x * and x oG = x o* and set BestPass = Pass If Pass = MaxPass then STOP, else ++Pass Let Max = Max(SumZeroØ(j), over j ∈ N(FC)) For all j ∈ N(FC) A. Let f j = SumZeroØ(j)/Max B. If SumZeroØ(j) > Max/2, v j = ⌈ f j ·U j ⌉ C. Else v j = Max( ⌈ f j ·U jo ⌉ , 1) D. p j = F j / v j Create and solve LP(p) to get new “first” test solution x’ A. Solve LP(p) by post-optimization to get x ′ and x o ′ B. Begin x* again from scratch to set x* = x ′ and x o * = x o ' C. Update U jo = max(U jo , x j ') for each j ∈ N(FC) D. Create the n-vector ZeroØ, where ZeroØ(j) = 1 if F j > 0 and x j ' = 0, else ZeroØ(j) = 0 1 E. Perform V_UPDATE Set First = 1 and Zero(1) = ZeroØ Set Zero(s) = (0, … 0) for s = 2 to sLim If Pass is a multiple of ZeroRefresh, then also re-initialize SumZeroØ = (0, … 0), but otherwise let SumZeroØ continue to accumulate END DIVERSIFY
Discussion of the Supporting Procedures
The DESCEND routine is the first supporting procedure invoked by the main routine, to implement the choice of x j* as the incoming pivot variable and the associated x k* as the leaving variable. If the algorithm is in a descent phase (Descent = True and TabuTenure = DescentTenure), and if the value x oj* continues the descent (x oj* < 0), then the routine simply performs the PIVOTJSTAR procedure which pivots in x j* and removes x k* from the basis tree, to produce the updated solution x" and its fixed charge objective x o ", and updates U jo for variables along the basis exchange path. Once the descent ends, Descent is set to False, TabuTenure is set to AscentTenure, and a check is performed to see if the solution x" (before updating by the basis exchange of x j* and x k* ) improves on x* ( x o " < x o * ). In this case, x* is updated as customary and the routine performs V_UPDATE, which updates the v j values as a foundation for subsequently determining the p j values that define the problem LP( p ). PIVOTJSTAR is likewise performed now that the descent ends. When the DESCEND routine is invoked and Descent = False, the PIVOTJSTAR routine is immediately performed and if x o " < x o * , then x* is updated as before. (The value x oj* can be improving after the initial descent has concluded. Instead of bouncing in and out of successive descent and ascent phases, once the initial descent has concluded, all subsequent steps are treated as an “ascent tabu phase.” However, TabuTenure is set to DescentTenure whenever an improving step occurs, and to AssetTenure otherwise.) Finally, Tabu(k*) = InsideIter + TabuTenure for the variable x k* that leaves the basis tree and becomes non-basic. Having discussed V_CHECK and PIVOTJSTAR in the explanation of DESCEND, it remains to discuss the supporting procedure DUPCHECK and the DIVERSIFY procedure that is invoked within it. The DUPCHECK routine is designed to check whether there are any duplications among the most recent ZeroØ vectors stored in Zero(s) for s = 1 to sLim. Since each ZeroØ vector identifies the variables x j that equal 0 in a given solution (by setting ZeroØ(j) = 1), and setting these variables to 0 automatically determines the network solution that sets remaining variables to 1, a duplication in these vectors implies that the associated fixed charge solutions are duplicated. DUPCHECK carries out a check for duplications (matches) by recording Zero(s) as a wraparound list, where the most recent ZeroØ vector is stored in Zero(First) and Zero(Last) is the ZeroØ vector recorded sLim iterations ago. The Zero(s) array starts from s= First until 2 reaching s = sLim, and then continues at s = 1 until reaching s = First – 1. Then the new (now most recent) ZeroØ vector is recorded by writing over the oldest one in the location s = First – 1 and then First is updated by setting First = First – 1. (Special case: If First = 1 then the location First – 1 is sLim.) This device avoids having to write the vectors into a temporary array and then write them back into Zero(s) to allow Zero(s) to always go from s = 1 to sLim. If the number of matches nMatch is found to exceed the limit LimMatch, the DIVERSIFY routine is executed which updates x G if the current x* improves upon it and if the DIVERSIFY routine has been invoked MaxPass times the algorithm stops. Otherwise the diversification proceeds by generating new f j values based on the formula f j = SumZeroØ(j)/Max, where SumZeroØ(j) counts the number of times x j = 0 in a solution that produced a ZeroØ vector in the DUPCHECK routine, and Max is the maximum of these SumZeroØ(j) values. The new v j values are then determined by setting v j = ⌈ f j ·U j ⌉ if SumZeroØ(j) > Max/2 and otherwise setting v j = Max( ⌈ f j ·U jo ⌉ , 1). From this, the p j values are determined by the usual formula p j = F j / v j as a basis for creating the problem LP(p) which is then solved by post-optimization to obtain a solution x'. The locally optimal solution x* starts again “from scratch” by setting x* = x', and the bounds U jo are updated in the customary way, along with establishing the ZeroØ vector as in the first step of the main algorithm. Finally, the V_UPDATE routine is executed, and the arrays associated with ZeroØ are likewise re-initialized, to conclude the DIVERSIFY procedure. In the event than Match is not True in the DUPCHECK procedure (and hence nMatch is not checked for exceeding LimMatch, and DIVERSIFY is not executed), then the DUPCHECK procedure updates values for tracking the algorithm’s performance, assures that nMatch = 0, and updates the Zero(s) array in accordance with the explanation above. In conjunction with the main routine, these supporting procedures complete the GI/TS algorithm.
4. GI/TS Computational Testing
An implementation of the above GI/TS algorithm, our code FixNetGI, was built using the alternating–path primal network simplex methods and data structures described in [1, 2, 3]. This solver is implemented in Fortran, compiled with gfortran –O3 , and tested under the Centos 6.10 version of the Linux operating system at Southern Methodist University. The test hardware is a Dell R720 with a Dual Six Core Intel Xeon @ 3.5GHz with all runs executed in single-thread mode. To assess the performance of FixNetGI, computational comparisons in terms of solution quality and speed are made with the IBM commercial optimization software Cplex 12.8, running with default parameters except for specifying single-threaded execution mode and a time limit per problem. Since Cplex is a general-purpose optimizer for linear and mixed-integer problems, the special-purpose heuristic approach of FixNetGI gives it major advantages. This comparison, however, is valuable because: no comparable solver for NetFC is available, Cplex is widely used 3 and respected by practitioners and researchers, and the comparison will indicate the heuristic’s efficiency and solution quality for use on real-world industry problems of this type. To test the effectiveness of the new solution approach, two problem test sets are used for benchmarking. The first is a collection of known problems from the literature and the second is a new suite of larger problems generated to explore the effects of problem characteristics on performance. Since there are over a dozen tuning parameters for the heuristic, we performed preliminary testing to identify a single set of parameters to use for all computational results reported herein. Randomly selected values from assigned ranges were run on the test sets, giving quite varied results, but providing guidance as to what value ranges seemed appropriate. The following parameter settings are employed for all runs reported: MaxIter = 50, MaxPass = 10, MaxInsideImprove = 40, BadLuck = 5, OutOfLuck = 20, Alpha(1) = 0.3, Alpha(2) = 0.45, Alpha(3) = 0.25, Beta = 0.4, MaxSolLimit = 1000, TabuTenure = 10, LimMatch = 10, sLim = 10, and ZeroRefresh = 30.
Test Set 1: Description
This first set of studied problems is drawn from the comprehensive FCTP testbed of Sun, et al [4] with a variety of problem dimensions and characteristics. The problems were originally created with a version of the well-known NETGEN random problem generator [6, 7], modified to include fixed costs on arcs. These Test Set 1 problems have seven problem dimensions, eight fixed-cost ranges (or types, labeled A-H), and 17 randomly generated instances of each combination. See Table 1 for definitions of these characteristics.
Table 1. Test Set 1 problem characteristics: (a) dimensions, (b) fixed cost range [4] Each test problem is a totally dense capacitated fixed-charge transportation problem with randomly distributed supplies and demands per Table 1(a) and with each arc randomly assigned a discrete variable cost between 3 and 8 plus a fixed cost in the associated range from Table 1(b). A subset of the 896 original testbed problems were selected for computational experiments with the GI2 code, following the choices of [5]. For the six smallest problem sizes, two instances of type A were used for this experimentation. For the largest and most difficult 50x100 size, all 15 instances of each fixed-charge type (A-H) were included, for a total of 132 problems. Hence the focus is on mixed-integer programs with 50,000 binary variables. Test Set 1: Computational Results and Analysis
Table 2 describes the solution results for the 12 smaller problems tested. Shown are the dimensions of the transportation problem, the problem identifier, the best solution value found and CPU solution time for Cplex 12.8 (run with a 7200-second time limit) and the FixNetGI code, the ratio of the two solvers’ solution values (Z-ratio = FixNetGI’s x oG /CPLEX (cid:1878) ∗ ) and the Cplex time as a multiple of the FixNetGI solution time (Time-X). Table 2 Test Set 1 s olution results for small problems, type A
With these smaller problems, the heuristic’s x oG solution values are within 0.1% of the Cplex optimal, on average, and were identified an average of three orders of magnitude faster. One third of FixNetGI’s solutions were optimal and its solution times averaged half a second. The bulk of the testing was focused on the more-difficult totally dense fixed-charge transportation problems with 50 source and 100 sink nodes, 50,000 arcs, supply of 50,000, and all fixed charge ranges as described in Table 1(b). Table 3 summarizes the results from solving 15 problem instances from each of the eight fixed-charge ranges (A-H). Detailed computational results from these 120 problems are found in Tables 4-11. The results on the larger problems underscore the effectiveness of the GI/TS algorithm. In every case, Cplex did not run to completion and exited at the 7,200-second time limit, while FixNetGI used an average of 1.11 seconds of CPU time. Although FixNetGI’s solution values averaged 9% higher, these were identified 6,000 times faster. To evaluate these solvers’ ability to handle even more challenging problems, such as is found in industrial applications, a new problem set was created. The problems are not only larger, but the suite is structured to facilitate statistical analysis of problem characteristics. Table 3. Test Set 1: summary of difficult, large 50x100 problems, averages of 15 problems per fixed-charge type
Table 4. Test Set 1: solution results for larger, difficult problems, type A fixed costs in range [50, 200] Table 5. Test Set 1: s olution results for large r, difficult problems, type B fixed costs in range [100,
Table 6. Test Set 1: solution results for larger, difficult problems, type C fixed costs in range [200, 800] Table 7. Test Set 1: solution results for larger, difficult problems, type D fixed costs in range [400, 1600]
Table 8. Test Set 1: solution results for larger, difficult problems, type E fixed costs in range [800, 3200] Table 9. Test Set 1: solution results for larger, difficult problems, type F fixed costs in range [1600, 6400]
Table 10. Test Set 1: solution results for larger, difficult problems, type G fixed costs in range [3200, 12800] Table 11. Test Set 1: solution results for larger, difficult problems, type H fixed costs in range [6400, 25600]
Test Set 2: Overview and Experimental design
To explore still larger problems and the possible effects of problem structure on solution time and quality, an experimental design using randomly generated test problems was established. For this, the NETGEN problem generator [9], modified to include fixed charges, created a new structured suite of transportation and transshipment problems with up to 33 times as many nodes, 100,000 binary variables, and a variety of problem characteristics. Test Set 2 consists of 96 problems, each generated with a different seed value, and with problem characteristics varied to enable a full-factorial experimental design. All combinations of five factors are used: number of problem nodes (500, 1000, 3000, and 5000), percentage of source and sink nodes (30% / 70%, transportation, and 20% / 20%, transshipment), number of arcs (10,000, 50,000, and 100,000), total supply (100,000 and 500,000), and fixed-cost range (20-200 and 1600-6400). All arcs have a fixed cost, a variable cost between 3-8, and an arc capacity between 200 and 1500 units. Transshipment sources and sinks are not used. Tables 12 and 13 display Test Set 2’s problem characteristics and solution results from the FixNetGI code and Cplex 12.8, run with a one-hour time limit and a single CPU thread. Problem characteristics shown are problem identifier and the number of nodes, sources and sinks, arcs, total supply, and fixed-cost range. Solution results are: the best solution value found (Best Z) for each application, the ratio of FixNetGI Z to Cplex Z (Z-ratio), the solution time using FixNet, and the Cplex time (3600 seconds in all instances) as a multiple of the FixNet solution time (Cplex Time-X). 1
Table 12 Test Set 2, 500- and 1000-node problem characteristics and solution results for FixNetGI and Cplex 12.8 Table 13 Test Set 2, 3000- and 5000-node problem characteristics and solution results for FixNetGI and Cplex 12.8 Table 14 Problem group and overall average Z-ratio, FixNetGI time, Cplex Time Multiple
Summary performance statistics by problem size and structure are given in Table 14. In terms of solution quality between the two solvers, The FixNetGI solution values average 1.2% larger than Cplex’s, but for 13 of the 96 problems FixNetGI solutions are superior (Z-ratio less than 1), including some larger instances where Cplex’s Best Z is 30 times larger. Based on average Z-ratio, the heuristic’s solution quality tends to be superior for transportation problems when compared to transshipment problems with the same number of nodes. In terms of solution speed, Cplex runs to the one-hour time limit in all cases. FixNetGI averages 10.1 seconds per problem, or 700 times faster than Cplex’s 3600-seconds, as shown in the Cplex Time-X column of Table 14. These multiples are better for the smaller problems, but all multiples would be much larger if Cplex had been allowed to run to optimality.
Test Set 2: Computational Results and Statistical Analysis
The structure of the test set enables rigorous statistical analysis of the relative performance of Cplex and FixNetGI solvers in terms of solution values and solution time, and the effect of the five factors described above. SAS 9.2’s analysis of variance procedure (ANOVA) and comparisons of means using Tukey’s Significant Difference (TSD) test are employed to determine whether the average results differed by solution method and whether factors affected the average results. The TSD procedure compares and ranks solver performance under the effect of different single-factor levels and treatment combinations. Specifically, we test hypotheses that the mean solution times and solution values are the same for both solvers and under different factor levels. Based on the problem solution times and values in Tables 12 and 13, ANOVA shows a statistically significant difference in mean solution times between the Cplex and FixNetGI codes. Hence, as expected, the mean solution speeds of the two solvers are statistically different, with FixNetGI being the faster. Statistical differences in time are also found between the four levels of problem node count, the two fixed-charge ranges, transportation and transshipment network 4 structures, the three levels of number of problem arcs, and two levels of total supply and demand. Hence, all hypotheses of equivalent means are rejected when runtime is the performance metric. However, when comparing solvers based on problem solution values (Z), the TSD test finds no statistically significant difference between the solvers. Therefore, while the mean Z-ratio for FixNetGI is slightly higher than Cplex’s, ANOVA shows that the mean solution values are not statistically different and the hypothesis of equality of mean solution values is not rejected. The two fixed-charge ranges do produce statistically different average solution values, as expected, but transportation and transshipment problems do not demonstrate statistically different values, nor do the numbers of problem arcs. Problems with 5000 nodes had mean solution values that are statistically different from those with 500 and 1000 nodes, but not those with 3000 nodes. This combination of hypothesis outcomes validates the effectiveness and speed of the GI/TS algorithm as implemented in FixNetGI for these larger and more challenging problem types. With solution times three orders of magnitude faster than Cplex while producing comparable objective function values, this approach advances the state of the art for fixed-charge network problems and renders solvable large practical instances from industrial settings.
5. Conclusions
Statistical testing reveals that the FixNetGI code is not only dramatically faster than Cplex in identifying its best solutions, but its mean solution quality is statistically equivalent to that of Cplex. This implementation of the GI/TS algorithm makes it appropriate for applications requiring high-quality results quickly, as in time-critical logistics, military response, airline re-scheduling, telecommunications and content-delivery network reconfiguration for demand fluctuations, and other near-real-time decision-making situations. There are a variety of opportunities to improve the GI/TS algorithm in the future. The tabu search procedure currently employed in the method is exceedingly simple, and a more advanced version may well enhance overall performance. Another conspicuous opportunity for future improvement will be to determine better parameters settings (for example, based on problem size and network class). A related possibility for investigation is to shortcut the Inside Loop operation and solve LP(p) more often, with the option of updating the solution each time by solving the restricted LP problem. Within the DUPCHECK procedure, the trade-offs between the sLim and the LimMatch values likewise invite examination, as do the values of the “alpha parameters” in V_UPDATE. The attractive outcomes produced by the current version of GI/TS embodied in FixNetGI provides a significant advance in our ability to solve fixed cost network problems efficiently and motivates a study devoted to the solution of practical problems in multiple areas. 5
References [1] Barr, R., F. Glover, and D. Klingman. “The Generalized Alternating Path Algorithm for Transportation Problems,”
European Journal of Operational Research
2, 137-144, 1978. [2] Barr, R., F. Glover, and D. Klingman. “Enhancements to Spanning Tree Labelling Procedures for Network Optimization,”
INFOR
Extremal Methods and Systems Analysis , A. V. Fiacco and K. O. Kortanek, eds., Springer-Verlag, New York, 250-274. [4] Sun, M., J.E. Aronson, P.G. McKeown, and D. Drinka, 1998. “A Tabu Search Heuristic Procedure for the Fixed Charge Transportation Problem.”
European Journal of Operational Research
Journal of Heuristics , 11, 307-336. [6] Klingman, D., Napier, A., Stutz, J., 1974. “NETGEN—A Program For Generating Large Scale (Un)Capacitated Assignment, Transportation, and Minimum Cost Flow Network Problems.”
Management Science
20 (5), 814–822. [7] Barr. R.S., Glover. F.. Klingman, D.. 1981. A New Optimization Method for Large-Scale Fixed Charge Transportation Problems,”
Operations Research
29 (3) 448--463. [8] Ahuja, R. K., T. L. Magnanti, and J. B. Orlin, 1993.
Network Flows: Theory, Algorithms, and Applications . Englewood Cliffs, N.J: Prentice Hall. [9] Bazaraa, M. S., John J. Jarvis, and Hanif D. Sherali, 2010.
Linear Programming and Network Flows, th ed. Hoboken, N.J: John Wiley & Sons. [10] Murty, K. G., 1992. Network Programming . Englewood Cliffs, N.J: Prentice Hall. [11] Glover, F., D Klingman, and N. V. Phillips, 1992.
Network Models in Optimization and Their Applications in Practice . New York: Wiley. [12] Glover, F. and M. Laguna, 1997.
Tabu Search . Kluwer Academic Publishers. Hingham, MA. [13] Balinski, M.L. (1961). “Fixed Cost Transportation Problem.”
Naval Research Logistics Quarterly
8, 41–54. 6 [14] Walker,W.E., 1976. “A Heuristic Adjacent Extreme Point Algorithm for the Fixed Charge Problem,”
Management Science
Network and Discrete Location : Models, Algorithms, and Applications , 2 nd ed. New York: Wiley. [17] Pioro, M., and D. Medhi, 2004. Routing, Flow, and Capacity Design in Communication and Computer Networks . San Francisco: Elsevier. [18] Mateus, G.R., and Z.K.G. Patrocinio, Jr., 2006. “Optimization Issues in Distribution Network Design,” in Resende, M.G.C. and F. M. Pardalos, eds.
Handbook of Optimization in Telecommunications , New York: Springer. [19] Hewitt, M., G. Nemhauser, and M. Savelsbergh, 2010. “Combining Exact and Heuristic Approaches for the Capacitated Fixed-Charge Network Flow Problem.”
INFORMS Journal on Computing
Computers & Industrial Engineering
99, 260–268. [22] Barr, R.S., R. Jones, and A. Klinkert, 2018. “An Efficient Optimization Approach to Designing Large-Scale Hierarchical Smart-Grid Data Networks,” technical report. [23] Alizadeh, M., 2009. “Facility Location in Supply Chain,” in R. Z., Farahani and M. Hekmatfar, eds.,
Facility Location: Concepts, Models, Algorithms and Case Studies , Berlin: Springer-Verlag. [24] Laport, G., S. Nickel, and F. Saldanha da Gama, eds, 2015.
Location Science . New York: Springer. [25] Fernández, E., and M. Landete, 2015. “Fixed-Charge Facility Location Problems,” in [24], 47—77. [26] Alumur, S. A., B. Y. Kara, and M. T.a Melo, 2015. “Location and Logistics.” In [24]. [27] Fortz, B., 2015. “Location Problems in Telecommunications,” in [24]. [28] Earle Steinberg and H. Albert Napier, 1980. “Optimal Multi-Level Lot Sizing for Requirements Planning Systems.”
Management Science
Operations Research [30] Eiselt, H.A., V. Marianov, and J. Bhadury, 2015. “Location Analysis in Practice,” in Eiselt, H.A. and V. Marianov, eds.,
Applications of Location Analysis . New York: Springer. [31] Kaan, L., and E. V. Olinick, 2013. “The Vanpool Assignment Problem: Optimization Models and Solution Algorithms,”
Computers & Industrial Engineering
66, 24–40. [32] Forsgren, A. and M. Prytz, 2006. “Telecommunications Network Design,” in Resende, M.G.C. and F. M. Pardalos, eds.