SOLVING THE NONLINEAR TRANSPORTATION PROBLEM BY GLOBAL OPTIMIZATION

. The aim of this paper is to present the suitability of three different global optimization methods for specifically the exact optimum solution of the nonlinear transportation problem (NTP). The evaluated global optimization methods include the branch and reduce method, the branch and cut method and the combination of global and local search strategies. The considered global optimization methods were applied to solve NTPs with reference to literature. NTPs were formulated as nonlinear programming (NLP) optimization problems. The obtained optimal results were compared with those got from literature. A comparative evaluation of global optimization methods is presented at the end of the paper to show their suitability for solving NTPs.


Introduction
The transportation problem (TP) is a well-known network optimization problem that was originally introduced by Hitchcock (1941). The objective of the TP is to determine the minimum cost distribution plan for the shipment of a single commodity from a number of sources to a number of destinations subject to the capacities of each source and each destination, Çakmak and Ersöz (2007). When transportation cost on a given route is nonlinearly dependent on the number of the units transported, the transportation problem becomes the nonlinear optimization problem. Determining the optimal solution of the nonlinear transportation problem (NTP) has been the subject of intensive research on logistics management. Various different approximate heuristic and exact mathematical programming methods have been proposed to solve the NTP.
As regards approximate heuristic optimization methods, genetic algorithms (GA) by Holland (1975), tabu search (TS) by Glover (1977), simulated annealing (SA) by Kirkpatrick et al. (1983), particle swarm optimization (PSO) by Kennedy and Eberhart (1995) and ant colony optimization (ACO) by Dorigo et al. (1996) have been proposed to solve the NTP in most of the cases. In this way, Michalewicz et al. (1991), Michalewicz (1994), Tsujimura et al. (2002) and Jo et al. (2007) have presented GA methods for the optimum solution of the NTP. Ilich and Simonovic (2001) have developed an evolution program (EP) for the NTP. TS optimization techniques for the NTP have been suggested by Cao and Uebe (1995) and Sun (1998). Yan and Luo (1999) have introduced the SA approach for the NTP. Li and Wang (2007) have devised the PSO technique to solve the NTP. Recently, the ACO method, the hybrid GA-ACO method and the hybrid TA-SA method for the NTP have been worked out in the papers by Altiparmak and Karaoglan (2006, 2007. Considering the exact mathematical programming methods, the optimum solution of the NTP has been frequently handled by different linear programming (LP) methods, see e.g. Cao (1992), Dangalchev (1996), Bell et al. (1999), Kuno and Utsunomiya (2000), Dangalchev (2000) and Nagai and Kuno (2005). However, research effort has been also devoted to nonlinear programming (NLP) techniques for the optimum solution of the NTP. For instance, Michalewicz et al. (1991) have applied the reduced gradient (RG) method to obtain the optimum solution of the NTP. Youssef and Mahmoud (1996) have described an iterative tangent line approximation procedure for solving the NTP under concave cost function. Tuy et al. (1996) and Puri and Puri (2006) have put forward a strongly polynomial algorithm for a concave NTP.
The aim of this paper is to present the suitability of three different global optimization methods for specifically the exact optimum solution of the NTP. The evaluated global optimization methods include the branch and reduce method, the branch and cut method and the combination of global and local search strategies. The considered optimization methods were applied to solve NTPs with reference to literature. NTPs were formulated as NLP optimization problems. The obtained optimal results were compared and a comparative evaluation of global optimization methods is presented at the end of the paper to show their suitability for solving NTPs.

Formulation of the NLP Optimization Problem
The general NLP optimization problem may be formulated in the following form: subject to: where: x is a vector of continuous variables defined within compact set X. Functions f(x), h(x) and g(x) are (non)linear functions involved in objective function z, equality and inequality constraints respectively. All functions f(x), h(x) and g(x) must be continuous and differentiable.
In the context of the NTP, continuous variables define the amounts of shipments (i.e. transporting flows) from sources to destinations. The objective function determines the total transportation cost. Equality and inequality constraints as well as the bounds of continuous variables represent a rigorous system of supply, demand and transporting flow constraints.

Formulation of the NLP Optimization Model
Taking into account the general formulation of the NLP optimization problem, the formulation of the latter model for the NTP is more specific, particularly in terms of variables and constraints. The objective of the NTP is to minimize the total nonlinear transportation cost while meeting supply, demand and transporting flow constraints. The formulation of the optimization model for the NTP is given in the following way: , 0, i j x ≥ for i = 1, 2,..., m and j = 1, 2,..., n, where: objective variable CT represents the total transportation cost for the shipment of a single commodity from sources m to destinations n, f i,j (x i,j ), denotes the cost function of transporting flow x i,j from source i to destination j, s i and d j are the capacities of each source i and each destination j respectively.
The total cost objective function defined by Eq. (1) represents the sum of the individual cost contributions of transporting flows. Eq. (2) and Eq. (3) determine supply and demand constraints respectively. While the set of constraints presented in Eq.
(2) ensures that the sum of shipments from a source equals its supply, the set of constraints given by Eq. (3) requires that the sum of shipments to a destination must satisfy its demand. Finally, the required non-negativity condition for transportation flows is handled by the set of constraints shown in Eq. (4). Only when all cost functions f i,j (x i,j ) are linear, the optimization model is linear. Otherwise, the above-defined optimization model determines the NTP.
The presented formulation of the optimization model determines the balanced transportation problem. The structure of the optimization model implies that total supply

Optimization Methods
The defined NLP optimization problem may be solved applying a suitable optimization method. The NLP class of optimization problems can, in principle, be solved using several classical local search algorithms and their extensions such as the reduced gradient method (RG) by Wolfe (1963), the generalized reduced gradient method (GRG) by Abadie and Carpentier (1969), augmented Lagrangian (AL) by Powell (1969) and Hestenes (1969), sequential quadratic programming (SQP) by Powell (1978) and the interior point method (IP) by Karmarkar (1984). In this paper, the suitability of the exact global optimization techniques for solving NTPs was investigated. The considered optimization techniques are briefly presented in the following sections.

Branch and Reduce Method
The first exact global optimization method applied is the branch and reduce method (BR) by Ryoo and Sahinidis (1996) implemented in the computational system BARON (Branch and Reduce Optimization Navigator) by Sahinidis and Tawarmalani (2008). The BR global optimization approach integrates the conventional branch and bound method (BB) with a wide variety of range reduction tests applied to every sub-problem of the search tree in pre-and post-processing steps to contract search space and reduce the relaxation gap. Many of reduction tests are based on duality and applied when relaxation is convex and solved by an algorithm that provides the dual, in addition to the primal, solution of the relaxed problem.
An important element of the computational system is the implementation of heuristic techniques for the approximate solution of the problem that yield improved bounds for the variables (feasibility-based tightening). It also includes a number of compound branching schemes that accelerate the convergence of standard branching strategies. Currently, the BR algorithm can handle nonlinear functions that involve exp(x), ln(x), x α for real α, β x for real β, x y and |x|. On the other hand, there is no support for other functions, including trigonometric ones such as sin(x), cos(x), arctan(x) etc. lems include 7×7 and 10×10 node problems with six different nonlinear continuous cost functions. The capacities of sources s i , the capacities of destinations d j and the matrix of the cost parameters c i,j of the 7×7 problem are given in Table 1.
It should be noted that 7×7 cost matrix represents a symmetrical matrix with zero cost parameters on the diagonal and six cost parameters with a large value of 1000 in a relative comparison to the rest of cost parameters. The input data of the 10×10 node problem is given in Table 2.
The cost functions applied in the tests were defined as proposed by Michalewicz et al. (1991) and are listed in Table 3.

Branch and Cut Method
The second exact global optimization technique applied is the branch and cut method (BC) implemented in solver LINDOGlobal by LINDO Systems, Inc. (2008).
The BC method is used to break an NLP optimization model down to a list of sub-problems each of which is analyzed and either: -is shown not to have a feasible or optimal solution; -an optimal solution to the sub-problem is found; -the sub-problem is further split into two or more sub-problems, which are then put on the list. The optimization algorithm automatically linearizes a number of nonlinear relationships through the addition of constraints and variables, so the transformed linearized model is mathematically equivalent to the original nonlinear model.
The optimization procedure has a multi-start feature that restarts the standard (non-global) nonlinear algorithm from a number of the intelligently generated points which allows the solver to find a number of locally optimal points and report the best one found.

Combination of Global and Local Search Strategies
The last global optimization approach employed is based on the seamless combination of global and local search strategies proposed by Pintér (2007) and implemented in optimization software LGO (Lipschitz Global Optimizer) by Pintér (2008). The optimizer integrates the following global scope algorithms: -branch and bound (adaptive partition and sampling) based global search (BB); -adaptive global random search (GARS); -adaptive multi-start global random search (MS) and the following local search (LS) strategies: -bound-constrained local search based on the use of an exact penalty function (EPF); -constrained local search based on a generalized reduced gradient approach (GRG). The global search algorithms are used to generate a close approximation of the global solution points (initial solution for subsequent local search). For this purpose, any of global algorithms BB, GARS or MS can be selected within a given solver run. Upon the generation of the global solution points, optimization switches over to LS using both the EPF and the GRG.
While BB + LS and GARS + LS search strategies launch a single local search from the best point found in the global search phase, MS + LS search strategy applies several local searches. Although MS + LS search strategy requires the most computational effort (due to its multiple local searches), it usually finds the best solutions. In this way, MS + LS search strategy was used to solve NTPs.

Input data
The considered set of test problems was originally presented by Michalewicz et al. (1991). The set of test prob-   Both values of parameters PA and PB were set to 1000. For function A, S was set to 2 while for functions B, E, and F, the value of 5 was used as proposed by Michalewicz et al. (1991).
The total transportation cost objective function for both above-defined 7×7 and 10×10 node problems was formulated as follows: where: f(x i,j ) is the cost function. While the shape of cost function f(x i,j ) is the same on all arcs, variation in arcs is obtained using cost parameters c i,j . Ilich and Simonovic (2001) have also applied the above-presented set of test problems to investigate the efficiency of a strongly feasible evolution program (SFEP) for solving NTPs. The gained optimum solutions presented in their research work were used in this research as the initial points for optimization.

Optimization
The formulation of the NLP optimization model for the NTP was applied to solve test problems. A high-level language GAMS (General Algebraic Modelling System) by Brooke et al. (1988) was used for modelling and data inputs/outputs. Test problems were solved using personal computer Intel Core2 Duo T8100, 2.10 GHz, 4 GB RAM DDR2 and 250 GB hard disc.
It should be noted that the n×n node optimization problem includes n 2 variables (e.g. transporting flows x i,j ) + an objective variable (e.g. total transportation cost CT), 2 n constraints (e.g. supply and demand constraints) and a nonlinear objective function, e.g. see Eq. (5). In this way, the resulting formulation of the NLP optimization model for the 7×7 node problem included 50 variables, 14 constraints and the objective function while the optimization model for the 10×10 node problem included 101 variables, 20 constraints and the objective function.

Computational Results for Function A
Function A represents a nonlinear continuous and differentiable arc-tangent approximation of a five-step piece-wise linear function, see Fig. 1.
Since BARON cannot handle trigonometric function arctan(x) within the optimization model, the optimum solutions of NTPs with function A were found only by LINDOGlobal and LGO. Both solvers found identical optimum solutions for the 7×7 node problem with function A, see Table 4. The gained minimum objective function value was 4.264. Table 5 and Table 6 represent optimum solutions for the 10 × 10 node problem with function A obtained using LINDOGlobal and LGO respectively. The algorithms calculated different optimum solutions with an identical minimum objective function value of 174.067.

Computational Results for Function B
Function B represents a nonlinear arc-tangent approximation of a piece-wise linear function with three gradients, see Fig. 2. Similarly as before, the optimum solutions of NTPs with function B were gained only by LINDOGlobal and LGO while BARON could not handle arc-tangent functions within the optimization model. Both LINDOGlobal and LGO found identical solutions for the 7×7 node problem and for the 10×10 node problem respectively. The minimum objective function value of 183.586 was obtained for the 7×7 node problem while the objective function value of the optimal solution for the 10×10 node problem was 146.987, see Table 7 and Table 8.

Computational Results for Function C
Function C is a regular quadratic function, see Fig. 3. The NTP with a quadratic cost function represents the convex optimization problem, for which classical gradient based algorithms usually easily find the global optimum.
As expected, BARON, LINDOGlobal and LGO solvers calculated the identical optimum solutions of the convex 7×7 and 10×10 node problems with function C respectively. Table 9 demonstrates the optimum solution of the 7×7 node problem with the objective function value of 2535.293 while Table 10 represents the optimum solution of the 10×10 node problem with the objective function value of 4401.649.

Computational Results for Function D
Function D is a square root function, see Fig. 4. This concave function gives large contribution to the overall objective function value even for small values of decision variables.
All three global optimization methods calculated the identical optimum solutions also in cases of 7×7 and 10×10 NTPs with function D respectively. For the 7×7 test problem, the applied global search algorithms found the optimum solution with the objective function value of 480.164, see Table 11.
The obtained minimum cost objective function value for the 10×10 node optimization problem was 388.910, see Table 12.

Computational Results for Function E
Function E infrequently appears within the cost functions of practical NTPs, see Fig. 5. However, this nonconvex function was proposed by Michalewicz et al. (1991) and included into NTPs as a severe test for optimization algorithms. Also, in cases of NTPs with function E, BARON, LINDOGlobal and LGO found the identical optimum solutions for 7×7 and 10×10 test problems respectively. Table 13 demonstrates the minimum cost solution of 204.842 for the 7×7 test problem.
For the 10×10 node problem, global optimization methods calculated the optimum solution with the cost objective function value of 71.657, see Table 14

Computational Results for Function F
Function F represents a non-convex function with multiple valleys and peaks, see Fig. 6. This function with multiple sub-optimal points usually causes difficulties in classical gradient based search techniques to find the global optimum solution.
Optimum solutions to NTPs with function F were found only using LINDOGlobal and LGO. BARON could not handle sinus functions within the optimization models of NTPs.
For the 7×7 test optimization problem, LINDOGlobal calculated a solution with the total transportation cost of 70.248 while LGO found the best solution with the minimum cost objective function value of 54.562, see Table 15 and Table 16. It should be noted that even after several hours of spent overall solver time, LINDOGlobal could not improve the minimum cost solution presented in Table 15.
On the other hand, both LINDOGlobal and LGO found identical optimum solutions for the 10 × 10 node optimization problem with function F, see Table 17. The calculated minimum value of the cost objective function for the optimum solution of the 10 × 10 node problem was 153.493.   Table 18 and Table 19 demonstrate a comparison between the obtained optimum results from this research and the previously reported optimum results obtained by the following optimization algorithms: -MINOS (RG) by Michalewicz et al. (1991); -GENETIC-2 (GA) by Michalewicz et al. (1991); -GENOCOP (GA) by Michalewicz (1994); -SFEP (EP) by Ilich and Simonovic (2001). Note, however, that the obtained optimum results for NTPs with functions A, B and F in Table 18 and Table  19 are presented only for global solvers LINDOGlobal and LGO since BARON could not handle trigonometric functions within the optimization model. Based on the above presented comparisons of the obtained results, the global optimization methods, except a few cases discussed bellow, found smaller or (almost) the same objective function values for the optimum solutions of NTPs in comparison with the optimization methods presented by Michalewicz et al. (1991), Michalewicz (1994 and Ilich and Simonovic (2001).

Comparison of Results
For instance, Table 18 shows that LINDOGlobal and the LGO calculated a higher value of the objective function for the 7×7 test problem with function A than GENETIC-2 and SFEP. The main reason for such results lies in the fact that the presented exact global optimization techniques can handle only the continuous and differentiable approximation of original piece-wise linear function A while both GENETIC-2 and SFEP can handle the original one. Contrary to the original piece-wise linear function A, arc-tangent approximation has no zero function values (e.g. see Fig. 1 for 0 ≤ x ≤ 2). Consequently, GENETIC-2 and SFEP calculated the zero minimum value of the objective function while LINDOGlobal and LGO calculated a very small non-zero minimum value of the objective function for the 7×7 test problem with function A. Note that the optimum solution of the 7×7 test problem calculated by LINDOGlobal and LGO indicates zero value for the objective function with original piece-wise linear function A, see cost parameters c i,j in Table 1, function A for 0 ≤ x ≤ 2 in Fig. 1 and the obtained optimum solution in Table 4.  Alongside the above mentioned differences in function formulations, conceptual differences between the considered optimization techniques also indicated some influence on the objective function values of the optimum solutions. Namely, GENETIC-2, GENOCOP and SFEP represent heuristic optimization methods while MINOS, BARON, LINDOGlobal and LGO represent the exact optimization methods. While the heuristic methods calculate the approximate optimum solutions, the exact optimization methods calculate the exact ones. The exact optimization methods are expected to find a global optimum of the convex optimization problems (e.g. test problems with quadratic function C) and the heuristic methods are usually gauged by measuring their closeness to global optimum found by the exact method.
In this way, the approximate optimum solution calculated by the heuristic method can sometimes represent an infeasible solution for the exact optimization technique. This, in turn, may result in some differences in the objective function value between both the approximate and exact solutions. For example, the approximate optimum solution obtained by SFEP indicates a 0.04% smaller objective function value for the convex 7×7 test problem with quadratic function C than the exact global optimum solution, see Table 18. Similarly, SFEP also found an approximate optimum solution for the 10 × 10 test problem with function A with a 0.62% smaller objective function value than LINDOGlobal and LGO, see Table 19.
On the other hand, LINDOGlobal and LGO calculated optimum solutions with a significantly smaller objective function value for NTPs with function A in comparison with GENOCOP (-82.36%) and MINOS (-38.05%). Further, the results for testing NTPs with functions B and F obtained by LINDOGlobal and LGO were found to be better in comparison with the best previously reported results, i.e. the optimum results obtained by SFEP.
For the 7×7 test problem with function B, the objective function value of the optimum solution obtained by LINDOGlobal and LGO was 9.89% smaller than the one obtained by SFEP. The minimum objective function value for the 10 × 10 test problem with function B found by LINDOGlobal and LGO was 8.01% smaller than that gained by SFEP.
Considering the test on NTPs with function F, the optimum solution for the 7×7 node problem obtained by LGO shows a 22.33% smaller objective function value than the solution found by LINDOGlobal. The minimum objective function value for the 7×7 test problem obtained by LGO was 30.77% smaller in comparison with the one found by SFEP. For the same test problem, LINDOGlobal calculated the optimum solution with a 10.86% smaller objective function value than that calculated by SFEP. As regards the 10×10 test problem with function F, both LINDOGlobal and LGO found the identical optimum solution with a 11.41% smaller objective function value than the solution calculated by SFEP. The optimum solutions for NTPs with functions C, D and E calculated by the considered global optimization techniques were found to be (nearly) the same as the best previously reported results.

Conclusions
All three global optimization methods find acceptable the exact solutions for the considered test on NTPs by reasonably low consumption of the total solver time. The gained solutions to the test on NTPs calculated applying global optimization techniques were found to be better or at least almost the same as the best previously reported results.
The BR algorithm inside computational system BARON appears to be fast and more robust of global optimization techniques. The convergence of the optimum solutions for the test on NTPs with functions C, D and E was achieved by the BR algorithm in less than a second. A lack of support for trigonometric functions within the optimization model reduces the applicability of the BR algorithm for solving NTPs.
The BC method implemented in solver LINDOGlobal, is an applicable tool for solving NTPs that can include a wide variety of nonlinear functions. However, non-convex functions with multiple valleys and peaks (i.e. such as function F) defined inside the optimization model for the NTP may cause difficulties with the BC method to find the global optimum solution to reasonable total solver time.
A combination of the global and local search strategies of LGO software (i.e. MS + LS search strategy) solves both convex and non-convex NTPs fast and efficiently. Most of the considered test problems were solved employing LGO within a few seconds. Only 10×10 test on the NTP with function A (i.e. nonlinear continuous and differentiable arc-tangent approximation of a five-step piece-wise linear function) was found to be the one that required the most consumption of the total solver time for which both LINDOGlobal LGO solvers required approximately one minute.