# A Hybrid Algorithm of Ant Colony and Benders Decomposition for Large-Scale Mixed Integer Linear Programming

Authors: Tingting Shan; Zhaoxuan Qiu
Journal: Computational Intelligence and Neuroscience
 (2022)

---

## Abstract

The mixed integer linear programming (MILP) has been widely applied in many fields such as supply chain management and robot control, while how to develop a more efficient algorithm to solve large-scale MILP is still in discussion. This study addresses a hybrid algorithm of the ant colony and Benders decomposition to improve the efficiency. We firstly introduce the design of our algorithm, in which the Benders algorithm decomposes the MILP into a master problem and a slack problem, the ant colony algorithm generates initial solutions for the master problem, and heuristic rules obtain feasible solutions for the slack problem. Then, the computational experiments are carried out to verify efficiency, with a benchmark test and some medium-large scale examples. Compared with other algorithms like CPLEX, GUROBI, and traditional ACA, our algorithm shows a better performance with a 0.3%–4.0% optimality gap, as well as a significant decrease of 54.3% and 33.6% on average in the CPU time and iterations, respectively. Our contribution is to provide a low-workload, time-saving, and high-accuracy hybrid algorithm to solve MILP problems with a large amount of variables, which can be widely used in more commercial solvers and promote the utilization of the artificial intelligence.


---

## Body

## 1. Introduction

Mixed integer linear programming (MILP) is a highly complex combinatorial optimization problem that contains both variables and continuous variables [1], as well as an NP-complete problem [2]. Due to the multiple types of its decision variables, it has broad applications in solving managerial planning problems, such as supply chain management [3], robot control problems [4], and traffic flow problems [5, 6]. A general MILP model is shown as follows:(1)PMinz=cx+fy,Ax≥b,Bx+Gy≥d,x∈Z+n1,y∈R+n2.In Model 1,z denotes the value of the objective function of the MILP, which often stands for some minimized planning goal like the total cost or the distance of path; x and y are the decision vectors subject to linear inequality constraints and the requirements in a mixed integer linear programming model, and all the variables are positive; n is the total amount of the decision variables, including n1 integer variables and n2 continuous variables, n=n1+n2; m is the amount of the constraints, including m1 constraints with only integer variables and m2 constraints with both integer and continuous variables, m=m1+m2; c and f is the vector of the continuous type of coefficients of the integer decision vector x and the continuous decision vector y, respectively, and the number of the elements in each vector equals to the amount of each kind of variables, c∈Rn1 and f∈Rn2; the matrices A, B, and G contain the coefficients in the constraints, which is related to all the variables in pure integer constraints, the integer variables and the continuous ones in mixed constraints, respectively. A∈Rm1×n1, B∈Rm2×n1,G∈Rm2×n2; b and d are the vectors including the constraint values, b∈Rm1 and d∈Rm2. If n1=n (n2=0), the problem is a pure integer linear program; if xi∈0,1, ∀I, the problem is a binary (0-1) integer program; if n2 ≥ 1, the problem is a MILP problem.Though MILP has a great significance in several managerial and computational contexts, there is a problem: the size of the MILP expands exponentially with the dimension of the decision variables; that is, the more the decision variables in the model, the more difficult it for the algorithm to run and work in CPU efficiently. When solving large-scale MILP problems, traditional algorithms will be hard to match the requirement of running in commercial solvers, which can be an obstacle of spreading the MILP model to the daily practices. Therefore, how to develop a more efficient algorithm to solve large-scale MILP problems is still in discussion in the academic community.To deal with this difficulty, previous research focused on three algorithms, i.e., decomposition algorithm, evolutionary algorithm, and hybrid algorithm. The decomposition algorithm is to decompose a MILP into an original problem and a slack problem, such as Benders decomposition [7–10], column generation [11–15], and Lagrangian decomposition [16–18]. The evolutionary algorithm solves a MILP in accordance with evolutionary rules [19–25]. However, both decomposition and evolutionary algorithm have their unique shortcomings in the calculation. To make use of the advantages of more than two algorithms and avoid all the drawbacks, the hybrid algorithm is proposed to solve MILP problems [26, 27].Our contributions include: (1) combining the both sides of advantages, we come up with a new hybrid algorithm based on the Benders decomposition with the ant colony algorithm and heuristic algorithm, which gives a new perspective to solving MILP in a special context of large-scale variables in an efficient way, trying to fill a research gap in developing a more efficient hybrid algorithm for large-scale MILP; (2) as our hybrid algorithm has excellent performances compared to other traditional algorithms, it can provide a low-workload, time-saving, and high-accuracy way to solve large-scale MILP problems in several managerial and computational contexts such as the logistics network management. As the requirement for the function of CPU is reduced, the cost paid in high-quality computing equipment will decrease significantly and thus the MILP can be widely used in more commercial solvers and significantly promote the utilization of the artificial intelligence in more and more parts of the life and production process.Our paper is organized as follows: in Section2, we review the literature relevant to our study. In Section 3, we introduce the solution procedure of the proposed hybrid algorithm. In Section 4, we carry out numerical experiments and report numerical results. In Section 5, we summarize the conclusions.


## 2. Literature Review

This study focuses on proposing a solution algorithm for MILP. The related research problems include the Benders decomposition algorithm, the hybrid algorithm of Benders decomposition, and evolutionary algorithms.The Benders decomposition algorithm refers to decomposing MILP into an original problem and a slack problem to improve the computational speed. Existing research has analysed the performance of Benders decomposition to solve MILP through solution time, iteration number, and degree of acceleration. Poojari and Beasley compare Benders decomposition with a standalone solver (CPLEX) and benders decomposition using a genetic algorithm [9]. Fischetti et al. prove that Benders decomposition allows for a significant boost in the performance of a mixed integer programming solver [28]. Víctor et al. use Benders decomposition to solve the coverage problem of network-oriented design [29]. François et al. use benders decomposition to solve the two-dimensional packing problem [30]. For the robot cell scheduling problem, Komari Alaei et al. prove that the Benders decomposition algorithm can improve the efficiency of scheduling [31]. Kergosiena et al. propose a Benders decomposition-based heuristic that makes it possible to find feasible solutions and lower bounds in a production and outbound distribution scheduling problem with strict delivery constraints [32]. However, the general Benders decomposition algorithm has a drawback as well. It takes a long time to calculate the optimal solution even in an infinite loop. Based on this situation, Boland et al. present a heuristic approach to two-stage mixed integer linear stochastic programming models with continuous second-stage variables [33]. In addition, to improve computation efficiency, the authors investigate the use of proximity search as a tactical tool to drive Benders decomposition. The Benders decomposition algorithm can effectively solve medium-large scale problems. Lin et al. propose a heuristic implementation of the Benders decomposition method that routes additional single and multiple flows without resolving the routing problem [34]. Previous studies have shown that the proposed hybrid algorithm is superior in convergence rate, CPU time, and iteration number.The hybrid algorithm incorporates evolutionary algorithm and heuristic rules into the benders decomposition framework. For the former, Awad et al. propose a mixed Benders decomposition approach with the ant colony optimization technique to solve a transfer line balancing problem [35]. This hybrid algorithm can get an optimal solution in a short time for the transfer line balancing medium-large scale problems. Ma and Zhang use multiobjective evolutionary algorithms to reveal the cost-efficiency balance of human brain networks [36]. Su et al. propose a multiobjective evolutionary algorithm for the detection of large-scale complex network communities [37]. For the letter, Osman and Baki propose a heuristic algorithm based on column generation to solve the vehicle routing problem with time windows [38]. Comparing the general benders decomposition with benders decomposition using a genetic algorithm, the hybrid algorithm is better than the general benders decomposition in terms of computational efficiency [9]. Liu and Dessouky propose a decomposition-based hybrid heuristic algorithm to solve the joint passenger and freight train scheduling problem [39].


## 3. Hybrid Algorithm of Ant Colony and Benders Decomposition

### 3.1. Basic Benders Decomposition

Benders decomposition algorithm provides a basic framework to solve MILP through decomposing the original complex problem into two problems, i.e., an original problem and a slack problem. The original problem and the slack problem could be solved in turn through the transmission of information [40]. Model 1 contains two decision problems, defined by integer decision variables x and continuous variables y, respectively. In this paper, the optimization on x is defined as the original problem or “master problem” (denoted as PM), and y is defined as a slack problem (denoted as PS). The slack problem can be written as follows:(2)PSMinz=cx^+fy,Gy≥d−Bx^,y∈R+n2.In Model 2,z is the objective function value of the PS; x^ denotes a feasible solution of x. When solving the slack problem PS, x^ can be treated as a constant vector temporarily. Therefore, the dual problem of model 2, denoted as PD, can be obtained as follows:(3)PDMaxz1=π¯Td−Bx^,π¯TG≤f,π≥0.In Model 3,z1 denotes the objective function value of the PD, which will equal to z only when both of them reach the optimal value; π denotes the dual variables of the PS. That means we can transfer PS into a standard form of a general linear programming problem like PD and then the optimal y will be easily solved.There are three possible solutions with respect toPS: (1) infeasible, then exit; (2) unbounded, in which case choose any unbounded extreme ray (denoted as π¯T) and add a feasibility cut π¯Td−Bx≤0 into the original problem. Since we need the minimum of z in the PS, the feasible cut can be seen as the upper limit of the objective function value of the PS; (3) bounded, in which case take an optimal solution (denoted as π˜T) and add an optimality cut π˜Td−Bx≤θ into the original problem. Due to the same reason in the unbounded solution, the optimality cut will be treated as the lower limit of the objective function value of the PS.Therefore, the original problem of Model 4, denoted asPM, can be written as follows:(4)PMMinw=cx+θ,Ax≥b,0≥π˜OTd−Bx,∀f=1,…,F,θ≥π˜OTd−Bx,∀o=1,…,O,x∈Z+n1,θ∈R.In thePM, w denotes the final objective function value of MILP. θ is a variable of the master problem, θ=fy and the y is the optimal solution vector for all the continuous variables we gain by solving the PS. F and O denote the collection of feasibility cut and optimality cut, respectively. The procedures of Benders decomposition are shown as follows (Algorithm 1):Algorithm 1: Procedures of benders decomposition.
Initialization: k: = 0, f = 0, o = 0; LB: = −∞, UB: = ∞; set an initial feasible solution x^, and algorithm stop condition ε = 0.while UB − LB > ε doSolve the PD.if Infeasible thenExit.end ifif Unbounded thenObtain an extreme rayπ∗;letf = f + 1;letπ¯=π∗;Add a feasibility cut 0≥π¯fTd−Bx to PM.end ifif Bounded thenObtain an extreme rayπ∗;leto = o + 1;letπ˜=π∗;Add an optimality cutθ≥π˜fTd−Bx to PM;Obtain a upper bound,UB=minUB,c⋅x^+π∗d−Bx^.end ifSolve PM.Obtain a solutionx∗,θ∗;Update the lower bound,LB=maxLB,cx∗+θ∗;k: = k + 1;x^: = x∗;end whileReturn the optimum result.Algorithm 2: Procedures of hybrid algorithm of ant colony algorithm and benders decomposition.
Initialization: k: = 0, f = 0, o = 0; LB: = −∞, UB: = ∞; set an initial feasible solution x^, and algorithm stop condition ε = 0.Generate the initial solution by ant colony algorithm:Letρ = 0.7, index = 1, Q = 104, m = 20, NCmax = 30;Setτmn: = ones (UBn − LBn + 1, n1 + n2)·(UBn − LBn + 1);fori: = 1 ⟶ mdoRandomly obtainedm solutions.end forFind local minimum solutionfmin;Updateτmn and pmn.while index < NCmax − 1 dofori: = 1 ⟶ m − 1 doRandomly obtainedm − 1 new solutions.end forAdd local minimum solution to the new solutions;Updatefmin, τmn and pmn;index = index + 1.end whileObtain the initial solutionx¯.Obtain upper bound by Benders decomposition:while UB − LB > εdoSolve the PD.if Infeasible thenExit.end ifif Unbounded thenObtain an extreme rayπ∗;letf = f + 1;letπ¯=π∗;Add a feasibility cut 0≥π¯fTd−Bx to PFM.end ifif Bounded thenObtain an extreme rayπ∗;leto = o + 1;letπ˜=π∗;Add an optimality cutθ≥π˜fTd−Bx to PFM;Obtain a upper bound, UB: = Min{UB,c⋅x^+π∗d−Bx^}.end ifObtain feasible solution and lower bound by feasible adjustment rule:Let maxiter: = 25,x0 = x¯, α0 = 1, μ = 0.9;Set UBFM: = ∞, and set the loop stop conditionλ = 10 − 4;Solve PFM and obtain the solution (x∗, θ∗);Letx∗ = round(x∗), α1 = µα0;Update UBFM: = Min{UBFM, ||x∗ − x0||2};fori: = 1 ⟶ maxiter − 1 doif ||xi + 1 − xi||2 > λthenLetxi = x∗;Solve PFM;Update the solution (x∗, θ∗);if ||xi + 1 − xi||2 < UBFM thenLet UBFM: = ||xi + 1 − xi||2;x∗ = round (xi + 1);end ifend ifUpdateαi + 1: = µαi.end forUpdate the lower bound, LB: = Max{LB,cx∗+θ∗};k: = k + 1;x^: = x∗;end whileReturn the optimum result


### 3.2. Hybrid Algorithm of Ant Colony Algorithm and Benders Decomposition

Medium-large scale MILP indicates it has a tremendous number of decision variables, which makes MILP more complicated and difficult to solve. Benders decomposition gives a framework to reduce the problem size, but it cannot guarantee the solution speed and convergence rate. To improve the solution speed of medium-large scale MILP, we propose a hybrid algorithm of the ant colony algorithm and Benders decomposition. It utilizes the ant colony algorithm to generate the initial solution for the master problem and uses a feasibility heuristic to get feasible solutions for the subproblem in each iteration. The hybrid algorithm flow chart is shown in Figure1.Figure 1
The flow chart of hybrid algorithm.The ant colony algorithm is a kind of heuristic algorithm imitating ants’ behavior. The idea is that a group of cooperating ants can find the shortest path between the nest and the food source. Ants achieve cooperation by leaving a certain amount of pheromone on the road. The pheromone left by an ant can be detected by other ants, and the more pheromones left on a road, the higher the probability of other ants will follow this path, and the path of the track will be more strengthened. The ant colony algorithm had been proposed for the MILP problems [41–43].For the MILP in Model 1, there aren1+n2, decision variables. We regard each solution of each decision variable as a node. The ants select a node which is a current value for each variable. Therefore, after n1+n2 times, it can obtain a solution for MILP. The solution procedures of the hybrid algorithm of the ant colony algorithm and Benders decomposition are explained as follows:Step 1:Initialization. Let m denote the number of the ant colonies. Therefore, it can generate m solutions according to Model 1. In each iteration, the newly updated pheromone concentration can be calculated as follows:(5)τmnnew=ρτmnold+Qfmin,whereτmnnew denotes the pheromone concentration of mth ant in nth variable in the current iteration,τmnold denotes the pheromone concentration of mth ant in nth variable in the previous iteration, ρ represents the disappearance rate of pheromone, Q is a constant which equals to 104 in this paper, fmin denotes the optimal value in the current iteration.Step 2:Fitness function. The fitness function is defined by the object function and the number of solution violating the constraints according to Model 1. It is defined as follows:(6)Ff=λf+1−λκ,whereFf denotes the fitness function, λ represents the proportion of each solution object value, f denotes the object value of each solution, κ denotes the number of solution violating the constraints.Step 3:Selection. Choose an optimal solution for Model 1 in the current iteration and update this optimal solution for the next iteration. For the ant choosing which solution node for each decision variable, the selection probability of mth ant in nth variable pmn is given by the following equation:(7)pmn=τmn∑n∈Nτmn∑,∀m=1,2,…,m,0,otherwise.Step 4:Feasibility adjustment. The generated initial solutions using the ant colony algorithm may not satisfy all the constraints in Model 4. There are two cases: the first one is feasible solutions if initial solutions satisfy all the constraints; the second one is infeasible solutions if initial solutions do not satisfy all the constraints. To guarantee the feasibility of the master problem, we propose a feasibility adjustment rule to solve Model 4. The solution idea is illustrated as follows.Letξ denotes the corresponding polyhedron of its LP relaxation for Model 4. The feasibility adjustment rule that obtains the feasible solution of the master problem can be constructed in the following:(8)PFMMinαIxI+1−xI2+1−αIcxI+1+θ,xI,θ∈ξ,whereαI is geometrically decreased with a fixed factor μ∈0,1, i.e., αI+1 = μαI and α0 = [0, 1], xI denotes the feasible integral solution in iteration I, and I=0 denotes the initial solution from ant colony algorithm output, and I>0 denotes the feasible solution from the feasibility heuristic in iteration I–1. We set the original α0 to 1. In order to generate a feasible solution, the feasibility heuristic method is also used by other scholars [1, 44–47].In general, the procedure of the hybrid algorithm of the ant colony algorithm and Benders decomposition is described as follows (Algorithm 2).


## 3.1. Basic Benders Decomposition

Benders decomposition algorithm provides a basic framework to solve MILP through decomposing the original complex problem into two problems, i.e., an original problem and a slack problem. The original problem and the slack problem could be solved in turn through the transmission of information [40]. Model 1 contains two decision problems, defined by integer decision variables x and continuous variables y, respectively. In this paper, the optimization on x is defined as the original problem or “master problem” (denoted as PM), and y is defined as a slack problem (denoted as PS). The slack problem can be written as follows:(2)PSMinz=cx^+fy,Gy≥d−Bx^,y∈R+n2.In Model 2,z is the objective function value of the PS; x^ denotes a feasible solution of x. When solving the slack problem PS, x^ can be treated as a constant vector temporarily. Therefore, the dual problem of model 2, denoted as PD, can be obtained as follows:(3)PDMaxz1=π¯Td−Bx^,π¯TG≤f,π≥0.In Model 3,z1 denotes the objective function value of the PD, which will equal to z only when both of them reach the optimal value; π denotes the dual variables of the PS. That means we can transfer PS into a standard form of a general linear programming problem like PD and then the optimal y will be easily solved.There are three possible solutions with respect toPS: (1) infeasible, then exit; (2) unbounded, in which case choose any unbounded extreme ray (denoted as π¯T) and add a feasibility cut π¯Td−Bx≤0 into the original problem. Since we need the minimum of z in the PS, the feasible cut can be seen as the upper limit of the objective function value of the PS; (3) bounded, in which case take an optimal solution (denoted as π˜T) and add an optimality cut π˜Td−Bx≤θ into the original problem. Due to the same reason in the unbounded solution, the optimality cut will be treated as the lower limit of the objective function value of the PS.Therefore, the original problem of Model 4, denoted asPM, can be written as follows:(4)PMMinw=cx+θ,Ax≥b,0≥π˜OTd−Bx,∀f=1,…,F,θ≥π˜OTd−Bx,∀o=1,…,O,x∈Z+n1,θ∈R.In thePM, w denotes the final objective function value of MILP. θ is a variable of the master problem, θ=fy and the y is the optimal solution vector for all the continuous variables we gain by solving the PS. F and O denote the collection of feasibility cut and optimality cut, respectively. The procedures of Benders decomposition are shown as follows (Algorithm 1):Algorithm 1: Procedures of benders decomposition.
Initialization: k: = 0, f = 0, o = 0; LB: = −∞, UB: = ∞; set an initial feasible solution x^, and algorithm stop condition ε = 0.while UB − LB > ε doSolve the PD.if Infeasible thenExit.end ifif Unbounded thenObtain an extreme rayπ∗;letf = f + 1;letπ¯=π∗;Add a feasibility cut 0≥π¯fTd−Bx to PM.end ifif Bounded thenObtain an extreme rayπ∗;leto = o + 1;letπ˜=π∗;Add an optimality cutθ≥π˜fTd−Bx to PM;Obtain a upper bound,UB=minUB,c⋅x^+π∗d−Bx^.end ifSolve PM.Obtain a solutionx∗,θ∗;Update the lower bound,LB=maxLB,cx∗+θ∗;k: = k + 1;x^: = x∗;end whileReturn the optimum result.Algorithm 2: Procedures of hybrid algorithm of ant colony algorithm and benders decomposition.
Initialization: k: = 0, f = 0, o = 0; LB: = −∞, UB: = ∞; set an initial feasible solution x^, and algorithm stop condition ε = 0.Generate the initial solution by ant colony algorithm:Letρ = 0.7, index = 1, Q = 104, m = 20, NCmax = 30;Setτmn: = ones (UBn − LBn + 1, n1 + n2)·(UBn − LBn + 1);fori: = 1 ⟶ mdoRandomly obtainedm solutions.end forFind local minimum solutionfmin;Updateτmn and pmn.while index < NCmax − 1 dofori: = 1 ⟶ m − 1 doRandomly obtainedm − 1 new solutions.end forAdd local minimum solution to the new solutions;Updatefmin, τmn and pmn;index = index + 1.end whileObtain the initial solutionx¯.Obtain upper bound by Benders decomposition:while UB − LB > εdoSolve the PD.if Infeasible thenExit.end ifif Unbounded thenObtain an extreme rayπ∗;letf = f + 1;letπ¯=π∗;Add a feasibility cut 0≥π¯fTd−Bx to PFM.end ifif Bounded thenObtain an extreme rayπ∗;leto = o + 1;letπ˜=π∗;Add an optimality cutθ≥π˜fTd−Bx to PFM;Obtain a upper bound, UB: = Min{UB,c⋅x^+π∗d−Bx^}.end ifObtain feasible solution and lower bound by feasible adjustment rule:Let maxiter: = 25,x0 = x¯, α0 = 1, μ = 0.9;Set UBFM: = ∞, and set the loop stop conditionλ = 10 − 4;Solve PFM and obtain the solution (x∗, θ∗);Letx∗ = round(x∗), α1 = µα0;Update UBFM: = Min{UBFM, ||x∗ − x0||2};fori: = 1 ⟶ maxiter − 1 doif ||xi + 1 − xi||2 > λthenLetxi = x∗;Solve PFM;Update the solution (x∗, θ∗);if ||xi + 1 − xi||2 < UBFM thenLet UBFM: = ||xi + 1 − xi||2;x∗ = round (xi + 1);end ifend ifUpdateαi + 1: = µαi.end forUpdate the lower bound, LB: = Max{LB,cx∗+θ∗};k: = k + 1;x^: = x∗;end whileReturn the optimum result


## 3.2. Hybrid Algorithm of Ant Colony Algorithm and Benders Decomposition

Medium-large scale MILP indicates it has a tremendous number of decision variables, which makes MILP more complicated and difficult to solve. Benders decomposition gives a framework to reduce the problem size, but it cannot guarantee the solution speed and convergence rate. To improve the solution speed of medium-large scale MILP, we propose a hybrid algorithm of the ant colony algorithm and Benders decomposition. It utilizes the ant colony algorithm to generate the initial solution for the master problem and uses a feasibility heuristic to get feasible solutions for the subproblem in each iteration. The hybrid algorithm flow chart is shown in Figure1.Figure 1
The flow chart of hybrid algorithm.The ant colony algorithm is a kind of heuristic algorithm imitating ants’ behavior. The idea is that a group of cooperating ants can find the shortest path between the nest and the food source. Ants achieve cooperation by leaving a certain amount of pheromone on the road. The pheromone left by an ant can be detected by other ants, and the more pheromones left on a road, the higher the probability of other ants will follow this path, and the path of the track will be more strengthened. The ant colony algorithm had been proposed for the MILP problems [41–43].For the MILP in Model 1, there aren1+n2, decision variables. We regard each solution of each decision variable as a node. The ants select a node which is a current value for each variable. Therefore, after n1+n2 times, it can obtain a solution for MILP. The solution procedures of the hybrid algorithm of the ant colony algorithm and Benders decomposition are explained as follows:Step 1:Initialization. Let m denote the number of the ant colonies. Therefore, it can generate m solutions according to Model 1. In each iteration, the newly updated pheromone concentration can be calculated as follows:(5)τmnnew=ρτmnold+Qfmin,whereτmnnew denotes the pheromone concentration of mth ant in nth variable in the current iteration,τmnold denotes the pheromone concentration of mth ant in nth variable in the previous iteration, ρ represents the disappearance rate of pheromone, Q is a constant which equals to 104 in this paper, fmin denotes the optimal value in the current iteration.Step 2:Fitness function. The fitness function is defined by the object function and the number of solution violating the constraints according to Model 1. It is defined as follows:(6)Ff=λf+1−λκ,whereFf denotes the fitness function, λ represents the proportion of each solution object value, f denotes the object value of each solution, κ denotes the number of solution violating the constraints.Step 3:Selection. Choose an optimal solution for Model 1 in the current iteration and update this optimal solution for the next iteration. For the ant choosing which solution node for each decision variable, the selection probability of mth ant in nth variable pmn is given by the following equation:(7)pmn=τmn∑n∈Nτmn∑,∀m=1,2,…,m,0,otherwise.Step 4:Feasibility adjustment. The generated initial solutions using the ant colony algorithm may not satisfy all the constraints in Model 4. There are two cases: the first one is feasible solutions if initial solutions satisfy all the constraints; the second one is infeasible solutions if initial solutions do not satisfy all the constraints. To guarantee the feasibility of the master problem, we propose a feasibility adjustment rule to solve Model 4. The solution idea is illustrated as follows.Letξ denotes the corresponding polyhedron of its LP relaxation for Model 4. The feasibility adjustment rule that obtains the feasible solution of the master problem can be constructed in the following:(8)PFMMinαIxI+1−xI2+1−αIcxI+1+θ,xI,θ∈ξ,whereαI is geometrically decreased with a fixed factor μ∈0,1, i.e., αI+1 = μαI and α0 = [0, 1], xI denotes the feasible integral solution in iteration I, and I=0 denotes the initial solution from ant colony algorithm output, and I>0 denotes the feasible solution from the feasibility heuristic in iteration I–1. We set the original α0 to 1. In order to generate a feasible solution, the feasibility heuristic method is also used by other scholars [1, 44–47].In general, the procedure of the hybrid algorithm of the ant colony algorithm and Benders decomposition is described as follows (Algorithm 2).


## 4. Computational Experiments

In this section, we present computational experiments to compare the efficiency and effectiveness of the proposed hybrid algorithm with the Benders decomposition algorithm and CPLEX solver. We adopt a mathematical model for the vehicle assignment and distribution problem in short supply to describe the general mixed integer linear programming. We use this model to test our proposed algorithm and report the numerical results from one small scale and six medium-large scale experiments, respectively. All the tests in this section were tested on a Lenovo Y400 with Intel Core i5-3230 M CPU, 2.60 GHz frequency, and 4 GB memory.


### 4.1. Vehicle Assignment and Distribution Problem

In our model formulation, we useS to denote the collection of all supply points, D denotes the collection of all demand points, and N denotes the collection of nodes. If node i, j∈S, node i and j represent the distribution center, otherwise, node i and j represent the demand point. Let cij represent the time cost from note i to note j. Let k represents the vehicle number, respectively. Let K denote the collection of vehicles. We use qk to represent the maximum loading capacity of the vehicle k. Let Qj denote the material demand at demand point j. There are two types of decision variables in the whole rescue vehicle dispatch period. Let xijk represent the flow originating at node i and arriving at node j. xijk is a binary decision variable. Let yjk represent the resource allocated to node j by vehicle k. The assumption that the total amount of rescue materials is less than the total amount of rescue materials at the demand point are used in this study and each vehicle cannot travel between supply points.The mathematical model for the vehicle assignment and distribution problem in short supply can be rewritten by the following equation:(9)Min∑i∈N∑j∈N∑k∈Kcij⋅xijk+∑j∈NQj−∑k∈Kyjk,s.t.∑i∈S∑j∈Dxijk=1,∀k∈K,∑i∈Dxijk−∑i∈Dxjik=0,∀j∈D,∀k∈K,∑j∈Nxijk≤1,∀i∈D,∀k∈K,∑k∈Kyjk≤Qj,∀j∈D,∑j∈Dyjk≤qk,∀k∈K,yjk≤Qj⋅∑i∈Nxijk,∀j∈D,∀k∈K,xijk∈0,1,yjk∈Z+.where ∑i∈N∑j∈N∑k∈Kcij⋅xijk denotes the total time cost of vehicles and ∑j∈NQj−∑k∈Kyjk denotes the total unmet material at demand points. In this optimization model, constraints 10∼12 are the flow conservation equations. Constraint 13 states that the quantity of materials released is smaller or equal to the quantity demand at each demand point. Constraint 14 ensures that the quantity of materials released is smaller and equal to its maximum load by each vehicle. Constraint 15 states the relationship between two variables. Finally, Constraint 16 ensures that xijk is a binary variable and yjk is the positive integer variable.


### 4.2. Small Scale Experiment: A Benchmark Test

All parameter values are shown in this section. In the example, we assume that the number of supply pointS is equal to 2, the number of demand point D is equal to 5, the number of vehicle K is equal to 3, and the material demand at demand points 1, 2, 3, 4, and 5 are 400, 500, 200, 600, and 300 separately, the maximum loading capacity of the vehicle 1, 2, and 3 are 800, 400, and 300, respectively, the time cost cij is shown in Table 1. Due to the assumption that each vehicle cannot travel between supply points, we cite a very large number M in Table 1.Table 1
Temperature and wildlife count in the three areas covered by the study.
S = 1S = 2D = 1D = 2D = 3D = 4D = 5S = 1MM12159624S = 2MM11812134D = 11211M186169D = 215818M152011D = 3912615M1714D = 4613162017M11D = 52449111411MThe results using CPLEX (version 12.6), GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm are shown in Table2 and Figure 2, respectively. According to the results shown in Table 2, the CPU times of the small scale experiment using CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm are 1.58 seconds, 1.21 seconds, 32.3 seconds, 0.39 seconds, 0.40 seconds, and 0.35 seconds, respectively. For one thing, our hybrid algorithm can save 99.8% of the CPU time in maximum compared with general Benders decomposition and 10.3% of the CPU time in minimum compared with ADMM, with an average decrease of 54.3%; for another, this algorithm can save 50% of the number of iterations in maximum compared with CPLEX and 14.3% of the number of iterations in minimum compared with ADMM, with an average decrease of 33.6%. Besides, we compare the upper bound, lower bound, CPU time, iteration number, and gap of six algorithms. The trend of the upper and lower bound for CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm is shown in Figure 2. Obviously, the upper and lower bounds of six algorithms are equal in the 12th, 11th, ninth, eighth, seventh, and sixth iteration, respectively. Therefore, this result proves that the efficiency of the proposed hybrid algorithm is better than the other five algorithms. The accuracy of CPLEX, GUROBI, and general Benders decomposition algorithm are better than the proposed hybrid algorithm in this small scale experiment. The accuracy of the proposed hybrid algorithm is better than ADMM and ACA in this small scale experiment.Table 2
Numerical results of benchmark example.
ApproachUpper boundLower boundCPU (s)IterationGap (%)CPLEX5445441.58120GUROBI5445441.21110Benders54454432.390ADMM5905310.3988.2ACA6025270.40710.6Hybrid algorithm5865370.3564Figure 2
The trend of the upper and lower bound.


### 4.3. Medium-Large Scale Experiments: Further Comparison

In order to better compare the advantages and disadvantages of six algorithms in this paper, we expand the number of supply point, the number of vehicle, and the number of demand point based on a benchmark example. Meanwhile, we expand other parameters in the model. There are 9 cases of medium-large scale experiments in the paper. We set the small scale experiment as Case 1. The numerical results of medium-large scale numerical examples that are solved by CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm are recorded in Tables3 and 4. The total physical nodes of the vehicle assignment and distribution network in the nine medium-large scale cases are 20, 25, 30, 35, 35, 40, 45, 50, and 60, respectively. Therefore, the total variable size of the vehicle assignment and distribution network in the nine medium-large scale cases are 1175, 2075, 3225, 4600, 6450, 9200, 12500, 16300, and 30825, respectively.Table 3
Numerical results of medium-large scale numerical examples.
ParametersVariable sizeCPLEXGUROBIBenders decompositionSDKxyUBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)253147155445441.5805445441.21054454432.305105112550325732573.200325732572.8103257325710205155200075431343134.600431343133.98043134313134.0052053125100677767775.660677767774.23067776777156.00102054500100705170516.050705170515.58070517051196.00520106250200788678866.840788678865.92078867886275.001020109000200800880087.320800880086.59080088008365.0010251012250250811781173510811781172350811781171221.001030101600030010307103072524010307103072012014270100013600.03815301530375450————————3440798533600.0145Table 4
Numerical results of medium-large scale numerical examples.
ParametersVariable sizeADMMACAHybrid algorithmSDKxyUBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)253147155905310.398.26025270.4010.65865370.354.05105112550330232141.801.4331432011.801.7329832311.700.75155200075439542872.501.9442042392.602.4432743062.400.352053125100689066514.202.5699966083.503.0681466993.601.2102054500100723469013.906.2798166813.908.3709269903.800.5520106250200801676294.203.3821075234.204.6790178424.100.51020109000200829978935.204.6846577685.605.7810479624.901.2102510122502508207801416.701.68287799517.302.18186808715.500.31030101600030011416996021.308.411409989320.3010.6104231016719.701.415301530375450145901302746.706.3154091258049.3010.1144071366443.202.6We can obtain upper bound, lower bound, CPU time, and gap by CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm. We report CPU time and iteration number amongst six algorithms for each case in Figures3 and 4, respectively. In order to better reflect the relationship between CPU time data, we do the followings: (1) divide the CPU time greater than 300 by 10 in Figure 3; (2) we set the CPU time for Case 10 by CPLEX and GUROBI equal to 5000 in Figure 3; (3) due to the CPU time of Benders decomposition is too long, we set a maximum time point to 3600 seconds. In order to better reflect the relationship between iteration number data, we do the following: set the iteration number for Case 10 by CPLEX and GUROBI equal to 150 in Figure 4. From the results of nine medium-large scale experiments, the calculation accuracy and calculation time are the same as that of small scale experiment. That is because Benders decomposition solves the master problem by using a function named CPLEXMILP in each iteration. The master problem has more and more cuts. Therefore, the running time of general Benders decomposition is longer than other algorithms for each example, which is not an appropriate technique if we solve the MILP problem with a commercial solver.Figure 3
CPU time amongst three algorithms for each case.Figure 4
Iteration number amongst three algorithms for each case.Summarily, compared to the Benders decomposition, ADMM, and ACA, the hybrid algorithm has better computational efficiency. The hybrid algorithm has better computational accuracy than ADMM and ACA in this paper. Obviously, the hybrid algorithm is more suitable for the MILP problem than CPLEX, GUROBI, Benders decomposition, ADMM, and ACA.


## 4.1. Vehicle Assignment and Distribution Problem

In our model formulation, we useS to denote the collection of all supply points, D denotes the collection of all demand points, and N denotes the collection of nodes. If node i, j∈S, node i and j represent the distribution center, otherwise, node i and j represent the demand point. Let cij represent the time cost from note i to note j. Let k represents the vehicle number, respectively. Let K denote the collection of vehicles. We use qk to represent the maximum loading capacity of the vehicle k. Let Qj denote the material demand at demand point j. There are two types of decision variables in the whole rescue vehicle dispatch period. Let xijk represent the flow originating at node i and arriving at node j. xijk is a binary decision variable. Let yjk represent the resource allocated to node j by vehicle k. The assumption that the total amount of rescue materials is less than the total amount of rescue materials at the demand point are used in this study and each vehicle cannot travel between supply points.The mathematical model for the vehicle assignment and distribution problem in short supply can be rewritten by the following equation:(9)Min∑i∈N∑j∈N∑k∈Kcij⋅xijk+∑j∈NQj−∑k∈Kyjk,s.t.∑i∈S∑j∈Dxijk=1,∀k∈K,∑i∈Dxijk−∑i∈Dxjik=0,∀j∈D,∀k∈K,∑j∈Nxijk≤1,∀i∈D,∀k∈K,∑k∈Kyjk≤Qj,∀j∈D,∑j∈Dyjk≤qk,∀k∈K,yjk≤Qj⋅∑i∈Nxijk,∀j∈D,∀k∈K,xijk∈0,1,yjk∈Z+.where ∑i∈N∑j∈N∑k∈Kcij⋅xijk denotes the total time cost of vehicles and ∑j∈NQj−∑k∈Kyjk denotes the total unmet material at demand points. In this optimization model, constraints 10∼12 are the flow conservation equations. Constraint 13 states that the quantity of materials released is smaller or equal to the quantity demand at each demand point. Constraint 14 ensures that the quantity of materials released is smaller and equal to its maximum load by each vehicle. Constraint 15 states the relationship between two variables. Finally, Constraint 16 ensures that xijk is a binary variable and yjk is the positive integer variable.


## 4.2. Small Scale Experiment: A Benchmark Test

All parameter values are shown in this section. In the example, we assume that the number of supply pointS is equal to 2, the number of demand point D is equal to 5, the number of vehicle K is equal to 3, and the material demand at demand points 1, 2, 3, 4, and 5 are 400, 500, 200, 600, and 300 separately, the maximum loading capacity of the vehicle 1, 2, and 3 are 800, 400, and 300, respectively, the time cost cij is shown in Table 1. Due to the assumption that each vehicle cannot travel between supply points, we cite a very large number M in Table 1.Table 1
Temperature and wildlife count in the three areas covered by the study.
S = 1S = 2D = 1D = 2D = 3D = 4D = 5S = 1MM12159624S = 2MM11812134D = 11211M186169D = 215818M152011D = 3912615M1714D = 4613162017M11D = 52449111411MThe results using CPLEX (version 12.6), GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm are shown in Table2 and Figure 2, respectively. According to the results shown in Table 2, the CPU times of the small scale experiment using CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm are 1.58 seconds, 1.21 seconds, 32.3 seconds, 0.39 seconds, 0.40 seconds, and 0.35 seconds, respectively. For one thing, our hybrid algorithm can save 99.8% of the CPU time in maximum compared with general Benders decomposition and 10.3% of the CPU time in minimum compared with ADMM, with an average decrease of 54.3%; for another, this algorithm can save 50% of the number of iterations in maximum compared with CPLEX and 14.3% of the number of iterations in minimum compared with ADMM, with an average decrease of 33.6%. Besides, we compare the upper bound, lower bound, CPU time, iteration number, and gap of six algorithms. The trend of the upper and lower bound for CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm is shown in Figure 2. Obviously, the upper and lower bounds of six algorithms are equal in the 12th, 11th, ninth, eighth, seventh, and sixth iteration, respectively. Therefore, this result proves that the efficiency of the proposed hybrid algorithm is better than the other five algorithms. The accuracy of CPLEX, GUROBI, and general Benders decomposition algorithm are better than the proposed hybrid algorithm in this small scale experiment. The accuracy of the proposed hybrid algorithm is better than ADMM and ACA in this small scale experiment.Table 2
Numerical results of benchmark example.
ApproachUpper boundLower boundCPU (s)IterationGap (%)CPLEX5445441.58120GUROBI5445441.21110Benders54454432.390ADMM5905310.3988.2ACA6025270.40710.6Hybrid algorithm5865370.3564Figure 2
The trend of the upper and lower bound.


## 4.3. Medium-Large Scale Experiments: Further Comparison

In order to better compare the advantages and disadvantages of six algorithms in this paper, we expand the number of supply point, the number of vehicle, and the number of demand point based on a benchmark example. Meanwhile, we expand other parameters in the model. There are 9 cases of medium-large scale experiments in the paper. We set the small scale experiment as Case 1. The numerical results of medium-large scale numerical examples that are solved by CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm are recorded in Tables3 and 4. The total physical nodes of the vehicle assignment and distribution network in the nine medium-large scale cases are 20, 25, 30, 35, 35, 40, 45, 50, and 60, respectively. Therefore, the total variable size of the vehicle assignment and distribution network in the nine medium-large scale cases are 1175, 2075, 3225, 4600, 6450, 9200, 12500, 16300, and 30825, respectively.Table 3
Numerical results of medium-large scale numerical examples.
ParametersVariable sizeCPLEXGUROBIBenders decompositionSDKxyUBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)253147155445441.5805445441.21054454432.305105112550325732573.200325732572.8103257325710205155200075431343134.600431343133.98043134313134.0052053125100677767775.660677767774.23067776777156.00102054500100705170516.050705170515.58070517051196.00520106250200788678866.840788678865.92078867886275.001020109000200800880087.320800880086.59080088008365.0010251012250250811781173510811781172350811781171221.001030101600030010307103072524010307103072012014270100013600.03815301530375450————————3440798533600.0145Table 4
Numerical results of medium-large scale numerical examples.
ParametersVariable sizeADMMACAHybrid algorithmSDKxyUBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)UBLBCPU (s)Gap (%)253147155905310.398.26025270.4010.65865370.354.05105112550330232141.801.4331432011.801.7329832311.700.75155200075439542872.501.9442042392.602.4432743062.400.352053125100689066514.202.5699966083.503.0681466993.601.2102054500100723469013.906.2798166813.908.3709269903.800.5520106250200801676294.203.3821075234.204.6790178424.100.51020109000200829978935.204.6846577685.605.7810479624.901.2102510122502508207801416.701.68287799517.302.18186808715.500.31030101600030011416996021.308.411409989320.3010.6104231016719.701.415301530375450145901302746.706.3154091258049.3010.1144071366443.202.6We can obtain upper bound, lower bound, CPU time, and gap by CPLEX, GUROBI, Benders decomposition, ADMM, ACA, and hybrid algorithm. We report CPU time and iteration number amongst six algorithms for each case in Figures3 and 4, respectively. In order to better reflect the relationship between CPU time data, we do the followings: (1) divide the CPU time greater than 300 by 10 in Figure 3; (2) we set the CPU time for Case 10 by CPLEX and GUROBI equal to 5000 in Figure 3; (3) due to the CPU time of Benders decomposition is too long, we set a maximum time point to 3600 seconds. In order to better reflect the relationship between iteration number data, we do the following: set the iteration number for Case 10 by CPLEX and GUROBI equal to 150 in Figure 4. From the results of nine medium-large scale experiments, the calculation accuracy and calculation time are the same as that of small scale experiment. That is because Benders decomposition solves the master problem by using a function named CPLEXMILP in each iteration. The master problem has more and more cuts. Therefore, the running time of general Benders decomposition is longer than other algorithms for each example, which is not an appropriate technique if we solve the MILP problem with a commercial solver.Figure 3
CPU time amongst three algorithms for each case.Figure 4
Iteration number amongst three algorithms for each case.Summarily, compared to the Benders decomposition, ADMM, and ACA, the hybrid algorithm has better computational efficiency. The hybrid algorithm has better computational accuracy than ADMM and ACA in this paper. Obviously, the hybrid algorithm is more suitable for the MILP problem than CPLEX, GUROBI, Benders decomposition, ADMM, and ACA.


## 5. Conclusions

This study focuses on proposing a new hybrid algorithm that contains ant colony and feasible heuristic decomposition with consideration of solving the large-scale MILP problems with higher computing efficiency. Firstly, we have used the ant colony algorithm and heuristic rules to calculate the initial and feasible solution for the master and slack problem, respectively, gaining the initial elements that are more approaching to the optimal solution of the MILP to speed up the Benders computing process. Subsequently, using vehicle assignment and distribution problem as a background, our hybrid algorithm has been applied in a series of computational experiments to verify its efficiency. The numerical results in the benchmark test have demonstrated that our hybrid algorithm significantly saves about 54.3% and 33.6% on average in the CPU time and iterations, respectively. Compared with other algorithms, our hybrid algorithm shows the superiority of computational performance with only 43.20s of CPU time, a 2.6% gap even in the largest 60-node example, which means our algorithm is always the fastest one and required the least iterations with a high robustness.Our research also has some limitations, which could pave the fruitful path for the future works. Firstly, the major limitation is that we only focus on extending the application of the Benders decomposition with ant colony in solving large-scale MILP problems in the reality due to the advantages indicated by previous studies. Besides, inefficiency may occur in coincidence when dealing with some data with very special structure due to the limitation of the random search technique applied in many heuristic algorithm. For future studies, more traditional algorithms such as CPLEX or GUROBI, etc., can be taken into consideration to improve the efficiency of other large-scale programming problems. Scholars can explore more new hybrid algorithms with other excellent metaheuristic algorithms, especially to avoid the coincident inefficiency when processing some special data, or develop new ways to deal with large-scale problems in the context of the Big Data epoch in the future.



---

*Source: 1634995-2022-08-31.xml*