# A Hybrid Approach to the Optimization of Multiechelon Systems

Authors: Paweł Sitek; Jarosław Wikarek
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

In freight transportation there are two main distribution strategies: direct shipping and multiechelon distribution. In the direct shipping, vehicles, starting from a depot, bring their freight directly to the destination, while in the multiechelon systems, freight is delivered from the depot to the customers through an intermediate points. Multiechelon systems are particularly useful for logistic issues in a competitive environment. The paper presents a concept and application of a hybrid approach to modeling and optimization of the Multi-Echelon Capacitated Vehicle Routing Problem. Two ways of mathematical programming (MP) and constraint logic programming (CLP) are integrated in one environment. The strengths of MP and CLP in which constraints are treated in a different way and different methods are implemented and combined to use the strengths of both. The proposed approach is particularly important for the discrete decision models with an objective function and many discrete decision variables added up in multiple constraints. An implementation of hybrid approach in theECLiPSe system using Eplex library is presented. The Two-Echelon Capacitated Vehicle Routing Problem (2E-CVRP) and its variants are shown as an illustrative example of the hybrid approach. The presented hybrid approach will be compared with classical mathematical programming on the same benchmark data sets.


---

## Body

## 1. Introduction

In the modern freight transportation there are two main distribution strategies: direct shipping and multiechelon distribution. In the direct shipping, vehicles, starting from a depot, bring their freight directly to the destination, while in the multiechelon systems, freight is delivered from the depot to the customers through an intermediate point.The majority of multiechelon systems presented in the literature usually explicitly consider the routing problem at the last level of the transportation system, while a simplified routing problem is considered at higher levels [1].In recent years multiechelon systems have been introduced in different areas:(i)
logistics enterprises and express delivery service companies under competitions;(ii)
hypermarkets and supermarkets products distribution;(iii)
multimodal freight transportation;(iv)
supply chains;(v)
delivery in logistic competition;(vi)
E-commerce and home delivery services under competitions;(vii)
city and public logistics.The vast majority of models of optimization in freight transportation and logistics industry have been formulated as the mixed integer programming (MIP) or mixed integer linear programming (MILP) problems and solved using the operations research (OR) methods [2]. Their structures are similar and proceed from the principles and requirements of mathematical programming (MP) [2, 3].Unfortunately, high complexity of decision-making models and their integer nature contribute to the poor efficiency of OR methods. Therefore a new approach to solving these problems was proposed. As the best structure for the implementation of this approach, a declarative environment was chosen [4, 5].It seems that better results will be obtained by the use of the declarative constraint programming paradigms (CP/CLP) especially in modeling. The CP-based environments have advantage over traditional methods of mathematical modeling in that they work with a much broader variety of interrelated constraints and allow producing “natural” solutions for highly combinatorial problems.The main contribution of this paper is hybrid approach (mixed CP with MP paradigms) to modeling and optimization of the Multi-Echelon Capacitated Vehicle Routing Problems or the similar vehicle routing problems. In addition, some extensions and modifications to the standard Two-Echelon Capacitated Vehicle Routing Problems (2E-CVRP) are presented.The paper is organized as follows. In Section2 the literature related to Multi-Echelon Vehicle Routing Problems has been reviewed. Next section is about our motivation and contribution. In Section 4 the concept of hybrid approach to modeling and solving and the solution hybrid framework have been presented. Then, the general description of Multi-Echelon Vehicle Routing Problems and mathematical model of 2E-CVRP has been discussed. Finally test instances for 2E-CVRP with extension variants and some computational results were discussed in Section 6.


## 2. Literature Review

The Vehicle Routing Problem (VRP) is used to design an optimal route for a fleet of vehicles to serve a set of customers’ orders (known in advance), given a set of constraints. The VRP is used in supply chain management in the physical delivery of goods and services. The VRP is of the NP-hard type.Nowadays, the VRP literature offers a wealth of heuristic and metaheuristic approaches, which are surveyed in the papers of [6, 7] because exact VRP methods have a size limit of 50–100 orders depending on the VRP variant and the time-response requirements.There are several variants and classes of VRP like the capacitated VRP (CVRP), VRP with Time Windows (VRPTW), and Dynamic Vehicle Routing Problems (DVRP), sometimes referred to as Online Vehicle Routing Problems and so forth [6].Different distribution strategies are used in freight transportation. The most developed strategy is based on the direct shipping: freight starts from a depot and arrives directly to customers. In many applications and real situations, this strategy is not the best one and the usage of a multiechelon and particular two-echelon distribution system can optimize several features as the number of the vehicles, the transportation costs, loading factor, and timing.In the literature the multiechelon system and the two-echelon system in particular refer mainly to supply chain and inventory problems [1]. These problems do not use an explicit routing approach for the different levels, focusing more on the production and supply chain management issues. The first real application of a two-tier distribution network optimizing the global transportation costs is due to [8] and is related to the city logistics area. They developed a two-tier freight distribution system for congested urban areas, using small intermediate platforms, called satellites (intermediate points for the freight distribution). This system is developed for a specific situation and a generalization of such a system has not already been formulated. The complete mathematic model of The Two-Echelon Capacitated Vehicle Routing Problem (2E-CVRP) with the solution for sample test data in the classical approach has been proposed by [7], complemented with the method for boosting the computing efficiency (see Section 5).The increasing role of supply chains and their urban parts evokes a need to focus greater attention on this issue in modeling and efficient optimization methods, in particular.


## 3. Motivation and Contribution

Based on [2, 5–7, 9–11] and our previous work [3, 4, 12] we observed some advantages and disadvantages of both (CP/MP) paradigms.An integrated approach of constraint programming/constraint logic programming (CP/CLP) and mixed integer programming/mixed integer linear programming (MIP/MILP) can help to solve optimization problems that are intractable with either of the two methods alone [13–15]. Although Operations Research (OR) and Constraint Programming (CP) have different roots, the links between the two environments have grown stronger in recent years.Approaches known from the literature are based mostly on the division of the main problem into sub-problems and iteratively solving each of them in the proper CP/CLP or MP/MILP technique. This is usually a collection of many local optimization points of feasible solutions. Other approaches are based on a “blind” transformation for the CLP to the MILP model. In most cases, this results in an explosion of the number of constraints and variables, which has a negative impact on the effectiveness of optimization. In the proposed hybrid approach, a very important element is the transformation of the initial problem and its solution in the field of domains, which takes place in CP/CLP environment. Then the converted and “slimmed down” problem is solved in the MILP environment, thus creating a global approach to optimization [14, 16].Both MIP/MILP and finite domain CP/CLP involve variables and constraints. However, the types of the variables and constraints that are used, and the way the constraints are solved, are different in the two approaches [13, 15].MIP/MILP relies completely on linear equations and inequalities in integer variables; that is, there are only two types of constraints: linear arithmetic (linear equations or inequalities) and integer (stating that the variables have to take their values in the integer numbers). In finite domain CP/CLP, the constraint language is richer. In addition to linear equations and inequalities, there are various other constraints such disequalities, nonlinear and symbolic (alldifferent,disjunctive,cumulative, etc.) constraints. In both MIP/ MILP and CP/CLP, there is a group of constraints that can be solved with ease and a group of constraints that are difficult to solve. The easily solved constraints in MIP/MILP are linear equations and inequalities over rational numbers.Integer constraints are difficult to solve using mathematical programming methods and often the real problems of MIP/MILP make them NP-hard.In CP/CLP, domain constraints with integers and equations between two variables are easy to solve. The system of such constraints can be solved over integer variables in polynomial time. The inequalities between two variables, general linear constraints (more than two variables), and symbolic constraints are difficult to solve, which makes real problems in CP/CLP NP-hard. This type of constraints reduces the strength of constraint propagation. As a result, CP/CLP is incapable of finding even the first feasible solution.Both environments use various layers of the problem (methods, the structure of the problem, data) in different ways. The approach based on mathematical programming (MIP/MILP) focuses mainly on the methods of optimization and, to a lesser degree, on the structure of the problem. However, the data is completely outside the model. The same model without any changes can be solved for multiple instances of data. In the approach based on constraint programming (CP/CLP), due to its declarative nature, the methods are already built-in. The data and structure of the problem are used for its modelling in a significantly greater extent.To use so much different environments and a variety of functionalities such as modeling, optimization, and transformation, the declarative approach was adopted.The motivation and contribution behind this work were to create a hybrid method for constrained decision problems modelling and optimization instead of using mathematical programming or constraint programming separately.It follows from the above that what is difficult to solve in one environment can be easy to solve in the other.Moreover, such a hybrid approach allows the use of all layers of the problem to solve it. In our approach, to modelling and optimisation, we proposed the environment, where:(i)
knowledge related to the problem can be expressed as linear, logical, and symbolic constraints;(ii)
the optimization models solved using the proposed approach can be formulated as a pure model of MIP/MILP or of CP/CLP, or it can also be a hybrid model;(iii)
the problem is modelled in the constraint programming environment by CLP-based predicates, which is far more flexible than the mathematical programming environment/very important for decision-making problems under competitions;(iv)
transforming the decision model to explore its structure has been introduced by CLP-based predicates;(v)
constrained domains of decision variables, new constraints, and values for some variables are transferred from CP/CLP into MILP/MIP/IP by CLP-based predicates;(vi)
optimization is performed by MP-based environment.As a result, a more effective hybrid solution environment for a certain class of decision and optimization problems (2E-CVRP or similar) was obtained.


## 4. A Hybrid Solution Framework for Capacitated Vehicle Routing Problems (HSFCVRP)

Both environments have advantages and disadvantages. Environments based on the constraints such as CLPs are declarative and ensure a very simple modeling of decision problems, even those with poor structures if any. In the CLP a problem is described by a set of logical predicates. The constraints can be of different types (linear, nonlinear, logical, binary, etc.). The CLP does not require any search algorithms. This feature is characteristic of all declarative backgrounds, in which modeling of the problem is also a solution, just as it is in Prolog, SQL, and so on. The CLP seems perfect for modeling any decision problem.Numerous MP models of decision-making have been developed and tested, particularly in the area of decision optimization. Constantly improved methods and mathematical programming algorithms, such as the simplex algorithm, branch and bound, and branch-and-cost, have become classics now [2].The proposed method’s strength lies in high efficiency of optimization algorithms and a substantial number of tested models. Traditional methods when used alone to solve complex problems provide unsatisfactory results. This is related directly to different treatment of variables and constraints in those approaches (Section3).This schema of the hybrid solution framework for Capacitated Vehicle Routing Problems (HSFCVRP) and the concept of this framework with its predicates (P1–P7) are presented in Figure1. The names and descriptions of the CLP predicates and the implementation environment are shown in Table 1.Table 1
Description of CLP predicates.
Predicate
Description
P1CLP environment
The implementation of the model in CLP, the term representation of the problem in the form of predicates.
P2CLP environment
The transformation of the original problem aimed at extending the scope of constraint propagation. The transformation uses the structure of the problem. The most common effect is a change in the representation of the problem by reducing the number of decision variables and the introduction of additional constraints and variables, changing the nature of the variables, and so forth.
P3CLP environment
Constraint propagation for the model: constraint propagation is one of the basic methods of CLP. As a result, the variable domains are narrowed, and in some cases, the values of variables are set, or even the solution can be found.
P4CLP environment
Generation by the AG:(i) the model for mathematical programming: generation performed automatically using CLP predicate;(ii)  additional constraints on the basis of the results obtained by predicate P3; (iii) domains for different decision variables and other parameters based on the propagation of constraints. Transmission of this information in the form of fixed value of certain variables and/or additional constraints to the MP.
Merging files generated by predicate AG into one file. It is a model file format in MP format.
P5EPLEX environment
Finding the consistent area based on information from the CLP.
P6EPLEX environment
The solution of the model from the P4 by MP solver.
P7EPLEX environment
Solution transfer from EPLEX to CLP (predicateeplex_get(vars,Zm))Figure 1
The scheme of the hybrid solution framework for Capacitated Vehicle Routing Problems (HSFCVRP).From a variety of tools for the implementation of the CP/CLP,ECLiPSe software [4, 12] of constraint programming applications.  ECLiPSe contains several constraint solver libraries, a high-level modelling and control language, interfaces to third-party solvers, an integrated development environment, and interfaces for embedding into host environments. ECLiPSe was used to model the problem, transform it, and search for a domain solution by constraint propagation. This solution was then the basis for the final MP model, developed in the Eplex library [9] of the ECLiPSe environment. Since ECLiPSe version 5.7, standalone Eplex have become the standard. The previous lib(eplex), which loads Eplex with the range bounds keeper and the IC variant have now been phased out, so users of these old variants must now move to using standalone Eplex. The Eplex library allows MP/MIP/MILP problems to be modelled in ECLiPSe and solved (optimized) by an external MP solver.


## 5. Two-Echelon Capacitated Vehicle Routing Problem as an Illustrative Example

The Two-Echelon Capacitated Vehicle Routing Problem (2E-CVRP) is an extension of the classical Capacitated Vehicle Routing Problem (CVRP) where the delivery depot-customers pass through intermediate depots (called satellites). As in CVRP, the goal is to deliver goods to customers with known demands, minimizing the total delivery cost in the respect of vehicle capacity constraints. Multiechelon systems presented in the literature usually explicitly consider the routing problem at the last level of the transportation system, while a simplified routing problem is considered at higher levels [7, 8].In 2E-CVRP, the freight delivery from the depot to the customers is managed by shipping the freight through intermediate depots. Thus, the transportation network is decomposed into two levels (Figure2): the 1st level connecting the depot (d) to intermediate depots (s) and the 2nd one connecting the intermediate depots (s) to the customers (c). The objective is to minimize the total transportation cost of the vehicles involved in both levels. Constraints on the maximum capacity of the vehicles and the intermediate depots are considered, while the timing of the deliveries is ignored.Figure 2
Example of 2E-CVRP transportation network.From a practical point of view, a 2E-CVRP system operates as follows (Figure2):(i)
freight arrives at an external/first/base zone, the depot, where it is consolidated into the 1st-level vehicles, unless it is already carried into a fully loaded 1st-level vehicles;(ii)
each 1st-level vehicle travels to a subset of satellites that will be determined by the model and then it will return to the depot;(iii)
at a satellite, freight is transferred from 1st-level vehicles to 2nd-level vehicles.


### 5.1. Mathematical Model

The formal mathematical model (MILP) was taken from [7]. Table 2 shows the parameters and decision variables of 2E-CVRP. Figure 2 shows an example of the 2E-CVRP transportation network for this model:
(1)min⁡∑i,j∈V0∪Vsci,j·Xi,j+∑k∈Vs∑i,j∈Vs∪Vcci,j·Yk,i,j+∑k∈Vssk·Dsk(2)∑i∈VsX0,i≤M1(3)∑j∈Vs∪V0,j≠kXj,k=∑i∈Vs∪V0,i≠kXk,ifork∈Vs∪V0(4)∑k∈Vs∑j∈VcYk,k,j≤M2(5)∑i∈Vc,j∈VcYk,i,j=∑i∈Vc,j∈VcYk,j,ifork∈Vs(6)∑i∈V0∪Vs,i≠jQi,j1-∑i∈Vs,i≠jQj,i1=Dsjjisnotthedeport∑i∈Vc-diotherwiseforj∈Vs∪V0(7)Qi,j1≤k1·Xi,jfori,j∈Vs∪V0,i≠j(8)∑i∈Vs∪Vc,i≠jQk,i,j2-∑i∈Vc,i≠jQk,j,i2=Zk,jdjjisnotasatelite-Djotherwiseforj∈Vc∪VS,k∈Vs(9)Qk,i,j2≤k2·Yk,i,jfori,j∈Vs∪Vc,i≠j,k∈Vs(10)∑i∈VsQi,V01=0(11)∑j∈VcQk,j,k2=0fork∈Vs(12)Yk,i,j≤Zk,jfori∈Vs∪Vc,j∈Vc,k∈Vs(13)Yk,i,j≤Zk,jfori∈Vs,j∈Vc,k∈Vs(14)∑i∈Vs∪VcYk,i,j=Zk,jfork∈Vs,j∈Vc,i≠k(15)∑i∈VsYk,j,k=Zk,jfork∈Vs,j∈Vc,i≠k(16)∑i∈VsZi,j=1forj∈Vc(17)Yk,i,j≤∑l∈Vs∪V0Xk,lfork∈Vs,i,j∈Vc(18)Yk,i,j∈0,1,Zk,l∈0,1fork∈Vs,i,j∈Vs∪Vc,l∈Vc(19)Xk,j∈Ζ+fork,j∈Vs∪V0(20)Qi,j1≥0fori,j∈Vs∪V0;Qk,i,j2≥0fori,j∈Vs∪Vc,k∈Vs(21)∑i,j∈ScYk,i,j≤Sc-1forSc⊂Vc,2≤Sc≤Vc-2(22)Qk,i,j2≤k2-dj·Yk,i,jfori,j∈Vc,k∈Vs(23)Qk,i,j2-∑l∈VsQk,j,l2≤k2-dj·Yk,i,jfori,j∈Vc,k∈Vs(24)Dsk=∑l∈Vcdj·Zk,jfork∈Vs.
The objective function minimizes the sum of the routing and handling operations costs. Constraints (3) ensure, for k=V0, that each 1st-level route begins and ends at the depot, while when k is a satellite, impose the balance of vehicles entering and leaving that satellite. Constraints (5) force each 2nd-level route to begin and end to one satellite and the balance of vehicles entering and leaving each customer. The number of the routes in each level must not exceed the number of vehicles for that level, as imposed by constraints (2) and (4). The flows balance on each network node is equal to the demand of this node, except for the depot, where the exit flow is equal to the total demand of the customers and for the satellites at the 2nd-level, where the flow is equal to the demand (unknown) assigned to the satellites which provide constraints (6) and (8). Moreover, constraints (6) and (8) forbid the presence of subtours not containing the depot or a satellite, respectively. In fact, each node receives an amount of flow equal to its demand, preventing the presence of subtours. Consider, for example, that a subtour is present between the nodes i,j, and k at the 1st level. It is easy to check that, in such a case, any value does not exist for the variables Qi,j1, Qj,k1, and Qk,i1, satisfying the constraints (6) and (8). The capacity constraints are formulated in (7) and (9), for the 1st-level and the 2nd-level, respectively. Constraints (10) and (11) do not allow residual flows in the routes, making the returning flow of each route to the depot (1st-level) and to each satellite (2nd-level) equal to 0. Constraints (12) and (13) indicate that a customer j is served by a satellite k(Zk,j=1) only if it receives freight from that satellite (Yk,i,j=1). Constraint (16) assigns each customer to one and only one satellite, while constraints (14) and (15) indicate that there is only one 2nd-level route passing through each customer and connect the two levels. Constraint (17) allows the start of a 2nd-level route from a satellite k only if a 1st-level route has served it. Constraints from (17) to (20) result from the character of the MILP-formulated problem. Additional constraints were introduced by [7] to increase the solution search efficiency. They strengthen the continuous relaxation of the flow model. In particular, authors in [7] used two families of cuts, one applied to the assignment variables derived from the subtour elimination constraints (edge cuts) and the other based on the flows. The edge-cuts explicitly introduce the well-known subtours elimination constraints derived from the TSP (Traveling Sales Problem). They can be expressed as constraint (21). The inequalities explicitly forbid the presence in the solution of subtours not containing the depot, already forbidden by constraint (8). The number of potential valid inequalities are exponential, so that each customer reduces the flow of an amount equal to its demand di-constraints (22) and (23).Table 2
Summary indices, parameters, and decision variables.
Symbol
Description
Indices
n
s
Number of satellites
n
c
Number of customers
V
0
-
[
v
0
]
Deport
V
s
=
{
v
s
1
,
v
s
2
,
…
,
v
sn
}
Set of satellites
V
c
=
{
v
c
1
,
v
c
2
,
…
,
v
cn
}
Set of customers
Parameters
M
1
Number of the 1st-level vehicles
M
2
Number of the 2nd-level vehicles
K
1
Capacity of the vehicles for the 1st level
K
2
Capacity of the vehicles for the 2nd level
d
i
Demand required by customeri
c
i
,
j
Cost of the arc(i,j)
s
k
Cost of loading/unloading operations of a unit of freight in satellitek
Decision variables
X
i
,
j
An integer variable of the 1st-level routing is equal to the number of 1st-level vehicles using arc(i,j)
Y
k
,
i
,
j
A binary variable of the 2nd-level routing is equal to 1 if a 2nd-level vehicle makes a route start from satellitek and go from node i to node j and 0 otherwise
Q
i
,
j
1
The freight flow arc(i,j) for the first level
Q
k
,
i
,
j
2
The freight arc(i,j) where k represents the satellite where the freight is passing through.
Z
k
,
j
A binary variable that is equal to 1 if the freight to be delivered to customerj is consolidated in satellite k and 0 otherwise


### 5.2. Model Transformation

One of the most important features that characterize the hybrid approach is the ease of modeling and transformation of the problem. The transformation is usually used to reduce the size of the problem and increase the efficiency of the search for a solution. In this case the transformation is based on the transition from arc to the route notation. During the transformation in the CLP the TSP, traveling salesman problem, is repeatedly solved and only the best routes in terms of costs are generated. In the process of transformation, the capacity vehicles constraints and those resulting from the set of orders are taken into account at both first and second level. For 2E-CVRP variants, time and logic constraints are also included.The obtained optimization model after the transformation(T1)–(T9) has different decision variables (Table 3) and different constraints than those in the MILP (1)–(24). Some of the decision variables are redundant; other variables are subject to aggregation. This results in a very large reduction in their number. Decision variables before and after the transformation are shown in Table 4. The transformation also reduces or eliminates some of the constraints of the model (Table 4):
(T1)min⁡∑k=1TcYkT·Kck+∑lTsXlT·Ksl(T2)∑k=1TcYkT≤M2(T3)∑k=1TcYkT·Wk,j=1forj=1,…,nc(T4)∑k=1TcXk·Wk,i·Pck=∑lTsPsl·Ul,ifori=1,…,ns(T5)∑i=1ns∑k=1TcYkT·Wk,i·Pck=∑lTsPsl(T6)XlT·K1≥Pslforl=1,…,Ts(T7)∑l=1TsXlT≤M1(T8)YkT∈0,1fork=1,…,Tc(T9)XlT∈Cforl=1,…,Ts.Table 3
Summary indices, parameters, and decision variables for transformed model.
Symbol
Description
Indices
n
s
Number of satellites
n
c
Number of customers
T
s
Number of possible routes from depot to satellites (CLP-determined)
T
c
Number of possible routes from satellites to customers (CLP-determined)
i
Satellite index
l
Depot-satellite route index
j
Customer index
k
Satellite-customer route index
M
1
Number of the 1st-level vehicles
M
2
Number of the 2nd-level vehicles
Input parameters
s
s
Cost of loading/unloading operations of a unit of freight in satellites
D
j
Demand required by customerj
P
c
k
Total demand for routek (CLP-determined)
K
s
l
Routel cost (CLP-determined)
K
c
k
Routek cost (CLP-determined)
U
l
,
i
Ifi is located on route l  Ul,i=1, otherwise Ul,i=0
W
k
,
s
If satellite or receipientsis located on route k  Wk,s=1, otherwise Wk,s=0
K
1
Capacity of the vehicles for the 1st level
Decision variables
Y
l
If the tour takes place along the routel from the route set generated for level 1, then Yl=1, otherwise Yl=0
X
k
If the tour takes place along the routek from the route set generated for level 2, then Xk=1, otherwise Xk=0
Computed quantities
P
s
l
Total demand for routelTable 4
Decision variables and constraints beforei after transformation.
Before transformation
After transformation
Description
Decision variables
X
i
,
j
X
l
T
Transformation of decision variables level 1 from the arc model arc(i,j) to the route model (l).
Q
i
,
j
1
Y
k
,
i
,
j
Y
k
T
Transformation of decision variables level 2 from the arc model arc(i,j) to the route model (k).
Q
k
,
i
,
j
2
Z
k
,
j
Constraints
(1)
(
T
1
)
Objective function after transformation, different decision variables, the same in terms of the essence and functionality.
(2)
(
T
7
)
Number of 1-type resources (CLP-determined)
(3)
—
Supply balance equation for 1-level nodes is unnecessary after transformation. This is a result of the route model to which particular vehicles are allocated.
(4)
(
T
2
)
Number of 2-type resources (CLP-determined)
(5)
—
Vehicle balance equation for level 2 is unnecessary after transformation. This is a result of the route model to which particular vehicles are allocated.
(6)
(
T
4
)
Supply balance for satellites.
(7)
(
T
6
)
Number of tours for level 1 resulting from the capacity of vehicles.
(8)
—
Supply balance constraint for recipients is not required. In the route model, the supply volume is calculated for the route.
(9)
—
Supply volume constraint resulting from the vehicle capacity is unnecessary for level 2. The routes are generated only for the allowable capacities.
(10)
—
No return loads from satellite to depot (10). The routes are generated so as to automatically ensure this.
(11)
—
No return loads from the customer to satellite (11). The routes are generated so as to automatically ensure this.
(12), (13)
—
Nozk,j variable after transformation.
(14)–(16)
(
T
3
)
No overlapping deliveries to customers.
(17)
—
This is ensured by the route model.
(18)–(20)
(
T
8
), (T9)
Integer and binary
(21)–(23)
—
Additional constraints are not necessary in the model with routes.


## 5.1. Mathematical Model

The formal mathematical model (MILP) was taken from [7]. Table 2 shows the parameters and decision variables of 2E-CVRP. Figure 2 shows an example of the 2E-CVRP transportation network for this model:
(1)min⁡∑i,j∈V0∪Vsci,j·Xi,j+∑k∈Vs∑i,j∈Vs∪Vcci,j·Yk,i,j+∑k∈Vssk·Dsk(2)∑i∈VsX0,i≤M1(3)∑j∈Vs∪V0,j≠kXj,k=∑i∈Vs∪V0,i≠kXk,ifork∈Vs∪V0(4)∑k∈Vs∑j∈VcYk,k,j≤M2(5)∑i∈Vc,j∈VcYk,i,j=∑i∈Vc,j∈VcYk,j,ifork∈Vs(6)∑i∈V0∪Vs,i≠jQi,j1-∑i∈Vs,i≠jQj,i1=Dsjjisnotthedeport∑i∈Vc-diotherwiseforj∈Vs∪V0(7)Qi,j1≤k1·Xi,jfori,j∈Vs∪V0,i≠j(8)∑i∈Vs∪Vc,i≠jQk,i,j2-∑i∈Vc,i≠jQk,j,i2=Zk,jdjjisnotasatelite-Djotherwiseforj∈Vc∪VS,k∈Vs(9)Qk,i,j2≤k2·Yk,i,jfori,j∈Vs∪Vc,i≠j,k∈Vs(10)∑i∈VsQi,V01=0(11)∑j∈VcQk,j,k2=0fork∈Vs(12)Yk,i,j≤Zk,jfori∈Vs∪Vc,j∈Vc,k∈Vs(13)Yk,i,j≤Zk,jfori∈Vs,j∈Vc,k∈Vs(14)∑i∈Vs∪VcYk,i,j=Zk,jfork∈Vs,j∈Vc,i≠k(15)∑i∈VsYk,j,k=Zk,jfork∈Vs,j∈Vc,i≠k(16)∑i∈VsZi,j=1forj∈Vc(17)Yk,i,j≤∑l∈Vs∪V0Xk,lfork∈Vs,i,j∈Vc(18)Yk,i,j∈0,1,Zk,l∈0,1fork∈Vs,i,j∈Vs∪Vc,l∈Vc(19)Xk,j∈Ζ+fork,j∈Vs∪V0(20)Qi,j1≥0fori,j∈Vs∪V0;Qk,i,j2≥0fori,j∈Vs∪Vc,k∈Vs(21)∑i,j∈ScYk,i,j≤Sc-1forSc⊂Vc,2≤Sc≤Vc-2(22)Qk,i,j2≤k2-dj·Yk,i,jfori,j∈Vc,k∈Vs(23)Qk,i,j2-∑l∈VsQk,j,l2≤k2-dj·Yk,i,jfori,j∈Vc,k∈Vs(24)Dsk=∑l∈Vcdj·Zk,jfork∈Vs.
The objective function minimizes the sum of the routing and handling operations costs. Constraints (3) ensure, for k=V0, that each 1st-level route begins and ends at the depot, while when k is a satellite, impose the balance of vehicles entering and leaving that satellite. Constraints (5) force each 2nd-level route to begin and end to one satellite and the balance of vehicles entering and leaving each customer. The number of the routes in each level must not exceed the number of vehicles for that level, as imposed by constraints (2) and (4). The flows balance on each network node is equal to the demand of this node, except for the depot, where the exit flow is equal to the total demand of the customers and for the satellites at the 2nd-level, where the flow is equal to the demand (unknown) assigned to the satellites which provide constraints (6) and (8). Moreover, constraints (6) and (8) forbid the presence of subtours not containing the depot or a satellite, respectively. In fact, each node receives an amount of flow equal to its demand, preventing the presence of subtours. Consider, for example, that a subtour is present between the nodes i,j, and k at the 1st level. It is easy to check that, in such a case, any value does not exist for the variables Qi,j1, Qj,k1, and Qk,i1, satisfying the constraints (6) and (8). The capacity constraints are formulated in (7) and (9), for the 1st-level and the 2nd-level, respectively. Constraints (10) and (11) do not allow residual flows in the routes, making the returning flow of each route to the depot (1st-level) and to each satellite (2nd-level) equal to 0. Constraints (12) and (13) indicate that a customer j is served by a satellite k(Zk,j=1) only if it receives freight from that satellite (Yk,i,j=1). Constraint (16) assigns each customer to one and only one satellite, while constraints (14) and (15) indicate that there is only one 2nd-level route passing through each customer and connect the two levels. Constraint (17) allows the start of a 2nd-level route from a satellite k only if a 1st-level route has served it. Constraints from (17) to (20) result from the character of the MILP-formulated problem. Additional constraints were introduced by [7] to increase the solution search efficiency. They strengthen the continuous relaxation of the flow model. In particular, authors in [7] used two families of cuts, one applied to the assignment variables derived from the subtour elimination constraints (edge cuts) and the other based on the flows. The edge-cuts explicitly introduce the well-known subtours elimination constraints derived from the TSP (Traveling Sales Problem). They can be expressed as constraint (21). The inequalities explicitly forbid the presence in the solution of subtours not containing the depot, already forbidden by constraint (8). The number of potential valid inequalities are exponential, so that each customer reduces the flow of an amount equal to its demand di-constraints (22) and (23).Table 2
Summary indices, parameters, and decision variables.
Symbol
Description
Indices
n
s
Number of satellites
n
c
Number of customers
V
0
-
[
v
0
]
Deport
V
s
=
{
v
s
1
,
v
s
2
,
…
,
v
sn
}
Set of satellites
V
c
=
{
v
c
1
,
v
c
2
,
…
,
v
cn
}
Set of customers
Parameters
M
1
Number of the 1st-level vehicles
M
2
Number of the 2nd-level vehicles
K
1
Capacity of the vehicles for the 1st level
K
2
Capacity of the vehicles for the 2nd level
d
i
Demand required by customeri
c
i
,
j
Cost of the arc(i,j)
s
k
Cost of loading/unloading operations of a unit of freight in satellitek
Decision variables
X
i
,
j
An integer variable of the 1st-level routing is equal to the number of 1st-level vehicles using arc(i,j)
Y
k
,
i
,
j
A binary variable of the 2nd-level routing is equal to 1 if a 2nd-level vehicle makes a route start from satellitek and go from node i to node j and 0 otherwise
Q
i
,
j
1
The freight flow arc(i,j) for the first level
Q
k
,
i
,
j
2
The freight arc(i,j) where k represents the satellite where the freight is passing through.
Z
k
,
j
A binary variable that is equal to 1 if the freight to be delivered to customerj is consolidated in satellite k and 0 otherwise


## 5.2. Model Transformation

One of the most important features that characterize the hybrid approach is the ease of modeling and transformation of the problem. The transformation is usually used to reduce the size of the problem and increase the efficiency of the search for a solution. In this case the transformation is based on the transition from arc to the route notation. During the transformation in the CLP the TSP, traveling salesman problem, is repeatedly solved and only the best routes in terms of costs are generated. In the process of transformation, the capacity vehicles constraints and those resulting from the set of orders are taken into account at both first and second level. For 2E-CVRP variants, time and logic constraints are also included.The obtained optimization model after the transformation(T1)–(T9) has different decision variables (Table 3) and different constraints than those in the MILP (1)–(24). Some of the decision variables are redundant; other variables are subject to aggregation. This results in a very large reduction in their number. Decision variables before and after the transformation are shown in Table 4. The transformation also reduces or eliminates some of the constraints of the model (Table 4):
(T1)min⁡∑k=1TcYkT·Kck+∑lTsXlT·Ksl(T2)∑k=1TcYkT≤M2(T3)∑k=1TcYkT·Wk,j=1forj=1,…,nc(T4)∑k=1TcXk·Wk,i·Pck=∑lTsPsl·Ul,ifori=1,…,ns(T5)∑i=1ns∑k=1TcYkT·Wk,i·Pck=∑lTsPsl(T6)XlT·K1≥Pslforl=1,…,Ts(T7)∑l=1TsXlT≤M1(T8)YkT∈0,1fork=1,…,Tc(T9)XlT∈Cforl=1,…,Ts.Table 3
Summary indices, parameters, and decision variables for transformed model.
Symbol
Description
Indices
n
s
Number of satellites
n
c
Number of customers
T
s
Number of possible routes from depot to satellites (CLP-determined)
T
c
Number of possible routes from satellites to customers (CLP-determined)
i
Satellite index
l
Depot-satellite route index
j
Customer index
k
Satellite-customer route index
M
1
Number of the 1st-level vehicles
M
2
Number of the 2nd-level vehicles
Input parameters
s
s
Cost of loading/unloading operations of a unit of freight in satellites
D
j
Demand required by customerj
P
c
k
Total demand for routek (CLP-determined)
K
s
l
Routel cost (CLP-determined)
K
c
k
Routek cost (CLP-determined)
U
l
,
i
Ifi is located on route l  Ul,i=1, otherwise Ul,i=0
W
k
,
s
If satellite or receipientsis located on route k  Wk,s=1, otherwise Wk,s=0
K
1
Capacity of the vehicles for the 1st level
Decision variables
Y
l
If the tour takes place along the routel from the route set generated for level 1, then Yl=1, otherwise Yl=0
X
k
If the tour takes place along the routek from the route set generated for level 2, then Xk=1, otherwise Xk=0
Computed quantities
P
s
l
Total demand for routelTable 4
Decision variables and constraints beforei after transformation.
Before transformation
After transformation
Description
Decision variables
X
i
,
j
X
l
T
Transformation of decision variables level 1 from the arc model arc(i,j) to the route model (l).
Q
i
,
j
1
Y
k
,
i
,
j
Y
k
T
Transformation of decision variables level 2 from the arc model arc(i,j) to the route model (k).
Q
k
,
i
,
j
2
Z
k
,
j
Constraints
(1)
(
T
1
)
Objective function after transformation, different decision variables, the same in terms of the essence and functionality.
(2)
(
T
7
)
Number of 1-type resources (CLP-determined)
(3)
—
Supply balance equation for 1-level nodes is unnecessary after transformation. This is a result of the route model to which particular vehicles are allocated.
(4)
(
T
2
)
Number of 2-type resources (CLP-determined)
(5)
—
Vehicle balance equation for level 2 is unnecessary after transformation. This is a result of the route model to which particular vehicles are allocated.
(6)
(
T
4
)
Supply balance for satellites.
(7)
(
T
6
)
Number of tours for level 1 resulting from the capacity of vehicles.
(8)
—
Supply balance constraint for recipients is not required. In the route model, the supply volume is calculated for the route.
(9)
—
Supply volume constraint resulting from the vehicle capacity is unnecessary for level 2. The routes are generated only for the allowable capacities.
(10)
—
No return loads from satellite to depot (10). The routes are generated so as to automatically ensure this.
(11)
—
No return loads from the customer to satellite (11). The routes are generated so as to automatically ensure this.
(12), (13)
—
Nozk,j variable after transformation.
(14)–(16)
(
T
3
)
No overlapping deliveries to customers.
(17)
—
This is ensured by the route model.
(18)–(20)
(
T
8
), (T9)
Integer and binary
(21)–(23)
—
Additional constraints are not necessary in the model with routes.


## 6. Computational Tests: Two-Echelon Capacitated Vehicle Routing Problem

For the final validation of the proposed hybrid approach, the benchmark data for 2E-CVRP was selected. 2E-CVRP, a well described and widely discussed problem, corresponded to the issues to which our approach was applied.The instances for computational examples were built from the existing instances for CVRP [17] denoted as E-n13-k4. All the instance sets can be downloaded from the website [18]. The instance set was composed of 5 instances with 1 depot, 12 customers, and 2 satellites. The full instance consisted of 66 instances because the two satellites were placed over twelve customers in all 66 possible ways (number of combinations: 2 out of 12). All the instances had the same position for depot and customers, whose coordinates were the same as those of instance E-n13-k4. The instances differed in the choice of two customers who were also satellites (En13-k4-4, En13-k4-5, En13-k4-6, En13-k4-12, etc.).Numerical experiments were conducted for the same data in three runs. The first run was a classical implementation of models (1)–(20) and its solution in the MP-based environment (P1). The second run used the same environment for models (1)–(24) with additional edge-cuts (P2). In the next run the models (1)–(20) were transformed (T1)–(T9) and solved in the proposed hybrid solution framework (P3). The calculations were performed using a computer with the following specifications: Intel(R) Core(TM) 2, 2×2, 40 GHZ RAM 2 GB. The analysis of the results for the benchmark instances demonstrates that the hybrid approach may be a superior approach to the classical mathematical programming. For all the examples, the solutions were found 4–40 times faster than they are in the classical approach (Table 5). In many cases the calculations ended after 600 s as they failed to indicate that the solution was optimal. The number of constraints (C) and decision variables/integer (V/int V) was, for example, P1, P2, and P3, respectively, C = 1262, V = 744/368 for P1, C = 1982, V = 744/368 for P2, and C = 21, V = 1082/1079 for P3. Thus, the combinatorial spaces (V×C) for illustrative examples were(i)
P1 ≈ 1 000 000;(ii)
P2 ≈ 1 400 000;(iii)
P3 ≈ 22 700.Table 5
The results of numerical examples for 2E-CVRP.
E-n13-k4
HSFCVRP (P3)
MP + Edge-Cuts (P2)
MP (P1)
T
Fc
T
Fc
T
Fc
E-n13-k4-01
17,36
280
600*
280
600*
280
E-n13-k4-02
17,22
286
600*
286
600*
286
E-n13-k4-03
15,39
284
600*
284
600*
284
E-n13-k4-04
10,09
218
44
218
65
218
E-n13-k4-05
9,58
218
48
218
108
218
E-n13-k4-06
11,05
230
78
230
154
230
E-n13-k4-07
9,16
224
39
224
64
224
E-n13-k4-08
13,03
236
46
236
75
236
E-n13-k4-09
13,22
244
67
244
93
244
E-n13-k4-10
14,08
268
107
268
183
268
E-n13-k4-11
18,91
276
159
276
600*
276
E-n13-k4-12
20,38
290
600*
290
600*
290
E-n13-k4-13
15,14
288
600*
288
600*
288
E-n13-k4-14
9,53
228
29
228
67
228
E-n13-k4-15
9,38
228
42
228
86
228
E-n13-k4-16
11,48
238
61
238
90
238
E-n13-k4-17
10,38
234
40
234
64
234
E-n13-k4-18
10,28
246
52
246
79
246
E-n13-k4-19
11,30
254
78
254
126
254
E-n13-k4-20
12,14
276
76
276
487
276
E-n13-k4-21
15,11
286
600*
286
600*
286
E-n13-k4-22
9,97
312
600*
312
600*
312
E-n13-k4-23
15,36
242
51
242
50
242
E-n13-k4-24
14,39
242
54
242
92
242
E-n13-k4-25
10,38
252
67
252
121
252
E-n13-k4-26
12,19
248
36
248
67
248
E-n13-k4-27
12,02
260
51
260
69
260
E-n13-k4-28
24,09
268
53
268
65
268
E-n13-k4-29
17,11
290
83
290
94
290
E-n13-k4-30
15,00
300
104
300
136
290
E-n13-k4-31
16,27
246
61
246
84
246
E-n13-k4-32
10,28
246
100
246
600*
246
E-n13-k4-33
15,17
258
93
258
123
258
E-n13-k4-34
11,00
252
48
252
55
252
E-n13-k4-35
8,92
264
40
264
52
264
E-n13-k4-36
11,11
272
97
272
138
272
E-n13-k4-37
16,06
296
109
296
213
296
E-n13-k4-38
16,69
304
124
304
600*
304
E-n13-k4-39
12,58
248
58
248
65
248
E-n13-k4-40
11,50
254
27
254
38
254
E-n13-k4-41
16,19
256
58
256
79
256
E-n13-k4-42
14,20
262
58
262
74
262
E-n13-k4-43
14,34
262
62
262
64
262
E-n13-k4-44
15,28
262
40
262
41
262
E-n13-k4-45
15,14
262
32
262
55
262
E-n13-k4-46
11,42
280
135
280
600*
280
E-n13-k4-47
12,20
274
95
274
142
274
E-n13-k4-48
13,17
280
76
280
257
280
E-n13-k4-49
11,16
280
79
280
117
280
E-n13-k4-50
12,30
280
63
280
83
280
E-n13-k4-51
14,97
280
48
280
62
280
E-n13-k4-52
15,30
292
63
292
98
292
E-n13-k4-53
12,33
300
66
300
150
300
E-n13-k4-54
14,28
304
94
304
600*
304
E-n13-k4-55
14,19
310
216
310
600*
310
E-n13-k4-56
17,05
310
60
310
162
310
E-n13-k4-57
14,13
326
221
326
600*
326
E-n13-k4-58
9,17
326
78
326
600*
326
E-n13-k4-59
12,02
326
56
326
112
326
E-n13-k4-60
13,91
326
42
326
68
326
E-n13-k4-61
12,20
338
600*
338
600*
338
E-n13-k4-62
10,05
350
79
350
365
350
E-n13-k4-63
11,92
350
83
350
239
350
E-n13-k4-64
10,13
358
122
358
600*
358
E-n13-k4-65
12,94
358
219
358
600*
358
E-n13-k4-66
11,91
400
600*
400
600*
400
*Calculations stopped after 600 s, the feasible value of the objective function.
Fc: the optimal value of the objective function.The logical relationship between mutually exclusive variables was taken into account, which in real-world distribution systems means that the same vehicle cannot transport two types of selected goods or two points cannot be handled at the same time.Those constraints result from technological, marketing, sales safety or competitive reasons. Only declarative application environments based on constraint satisfaction problem (CSP) make it possible to implement of this type of constraint. Table6 presents the results of the numerical experiments conducted for 2E-CVRPs with logical constraints relating to the situation where two delivery points (customers) can be handled separately but not together in one route.Table 6
The results of numerical examples for 2E-CVRP with logical constraints.
E-n13-k4
Fc
T
C
V
exCustomer*
E-n13-k4-01
284
15,36
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-07
240
7,16
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-11
290
16,91
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-20
280
13,14
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-26
270
10,72
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-32
270
10,88
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-33
276
14,124
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-40
284
11,23
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-46
308
11,12
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
E-n13-k4-54
334
14,28
21
788
2,3; 2,4; 2,6; 2,7; 1,8; 1,9
*Pairs of customers that cannot be served on one route.The final stage of the research was to optimize Two-Echelon Capacitated VRP with Time Windows (2E-CVRP-TW). This problem is the extension of 2E-CVRP where time windows on the arrival or departure time at the satellites and/or at the customers are considered. The time windows can be hard or soft. This variant of the 2E-CVRP is extremely important in a competitive environment.In the first case the time windows cannot be violated, while in the second one if they are violated a penalty cost is paid. 2E-CVRP-TW has been implemented in a hybrid environment. This was followed by the optimization problem under the time constraints (time windows). There have been experiments with both windows hard and windows soft. The results are shown in Table7(a) and Table 7(b). The impact of these constraints on the value of the objective function can be clearly seen. For instance E-n13-k4-20 also shows graphically the optimal way of delivery (see Figure 3(a)) as well as the impact on the route of logical constraints (Figure 3(b)) and the time window (Figure 3(c)).Table 7
(a) The results of numerical examples for 2E-CVRP-TW (hard windows). (b) The results of numerical examples for 2E-CVRP-TW (soft windows, penalty = 30).
(a)
E-n13-k4
T
40
50
60
70
80
90
100
110
130
150
160
E-n13-k4-01
—
—
—
—
—
—
280
280
280
280
280
E-n13-k4-07
—
224
224
224
224
224
224
224
224
224
224
E-n13-k4-11
—
—
304
276
276
276
276
276
276
276
276
E-n13-k4-20
—
294
280
276
276
276
276
276
276
276
276
E-n13-k4-26
—
248
248
248
248
248
248
248
248
248
248
E-n13-k4-32
—
—
262
246
246
246
246
246
246
246
246
E-n13-k4-33
—
258
258
258
258
258
258
258
258
258
258
E-n13-k4-40
—
284
284
254
254
254
254
254
254
254
254
E-n13-k4-46
—
—
308
308
280
280
280
280
280
280
280
E-n13-k4-54
—
—
—
324
304
304
304
304
304
304
304
(b)
E-n13-k4
T
40
50
60
70
80
90
100
110
130
150
160
E-n13-k4-01
358
354
346
346
310
310
280
280
280
280
280
E-n13-k4-07
270
224
224
224
224
224
224
224
224
224
224
E-n13-k4-11
306
306
304
276
276
276
276
276
276
276
276
E-n13-k4-20
366
294
280
276
276
276
276
276
276
276
276
E-n13-k4-26
292
248
248
248
248
248
248
248
248
248
248
E-n13-k4-32
322
278
262
246
246
246
246
246
246
246
246
E-n13-k4-33
336
258
258
258
258
258
258
258
258
258
258
E-n13-k4-40
344
284
284
254
254
254
254
254
254
254
254
E-n13-k4-46
344
310
308
308
280
280
280
280
280
280
280
E-n13-k4-54
342
334
334
324
304
304
304
304
304
304
304(a) Example of 2E-CVRP transportation network for E-n13-k4-20 instance. (b) Example of 2E-CVRP transportation network for E-n13-k4-20 instance with logic constraints. (c) Example of 2E-CVRP-TW transportation network for E-n13-k4-20 instance.
(a)
(b)
(c)


## 7. Conclusion and Discussion on Possible Extension

The efficiency of the proposed approach is based on the reduction of the combinatorial problem and using the best properties of both environments. The hybrid approach (Table3) makes it possible to find solutions in the shorter time.In addition to solving larger problems faster, the proposed approach provides virtually unlimited modeling options with many types of constraints. Therefore, the proposed solution is recommended for decision-making problems under competitions and that has a structure similar to the presented models (Section5). This structure is characterized by the constraints and objective function in which the decision variables are added together.Further work will focus on running the optimization models with nonlinear and other logical constraints, multiobjective, uncertainty, and so on, in the hybrid optimization framework. The planned experiments will employ HSFCVRP for Two-Echelon Capacitated VRP with Satellites Synchronization, 2E-CVRP with Pickup and Deliveries, and other VRP issues in Supply Chain Sustainability [19] and other routing problems [20].



---

*Source: 925675-2015-01-26.xml*