# Model and Algorithm for Container Allocation Problem with Random Freight Demands in Synchromodal Transportation

Authors: Yan Xu; Chengxuan Cao; Bin Jia; Guangzhi Zang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper aims to investigate container allocation problem with random freight demands in synchromodal transportation network from container carriers’ perspective. Firstly, the problem is formulated as a stochastic integer programming model where the overall objective is to determine a container capacity allocation plan at operational level, so that the expected total transportation profit is maximized. Furthermore, by integrating simulated annealing with genetic algorithm, a problem-oriented hybrid algorithm with a novel gene encode method is designed to solve the optimization model. Some numerical experiments are carried out to demonstrate the effectiveness and efficiency of the proposed model and algorithm.


---

## Body

## 1. Introduction

As stated in [1], synchromodal freight transportation involves a structured, efficient, and synchronized combination of two or more transportation modes. Through synchromodal transportation, the carriers or customers select independently at any time the best transportation mode based on the operational circumstances and/or customer requirements. While multimodal freight transportation is defined as the transportation of goods by a sequence of at least two different modes of transportation, intermodal freight transportation is defined as a particular type of multimodal transportation where the load is transported from an origin to a destination in one and the same intermodal transportation unit. There obviously exist some differences in synchromodal freight transportation, multimodal freight transportation, and intermodal freight transportation as stated above. Besides, resource management, which is on how and when to optimally utilize the limited available resources, is a typical problem at operational level in freight transportation network. Some examples of these resources are containers, vehicles (e.g., planes and trucks), trailers, rail cars, locomotives, equipment, crew, power, and so forth. Once a resource is allocated to an activity, it is no longer available for certain duration.Containerization may improve safety, reduction of handling costs, standardization, and accessibility to multiple modes of transportation [2]. It plays a more and more important role in freight transportation network. While container central station is a large container handling station set up in accordance with the construction of modern logistics hub requirements, it is also a terminal station of loading and unloading containerized freight to maritime, hinterland, or air transportation. At present, many carriers have to determine the most profitable allocation of available capacity for utilizing limited container more efficiently. As pointed in [3], it is still a challenging task to optimize transportation profit by optimizing resource management. These reasons inspire us to research container capacity allocation problem. Since containerization transportation system is a considerable complex system, including container allocation, reposition, and retention, it is sometimes difficult to optimize carriers’ transportation profit by determining the scheme of container capacity allocation for container carries in strategic stage [4]. Therefore, we mainly focus on container capacity allocation problem for carriers at operational level in this paper, that is, how to allocate the limited container to maximize transportation profit in synchromodal freight transportation network.So far, many scholars have researched on capacity problem in transportation network. Bostel and Dejax [5] defined the capacity allocation problem as allocating capacity to match each demand. Demirag and Swann [6] considered capacity allocation problem for central headquarters by formulating it as a mixed integer program, and some simple heuristic algorithms are present to solve it. Bell et al. [7] focused on container assignment problem by taking empty containers flow into consideration and introduced a frequency based model for a maritime network with port rotations services. Arnone et al. [8] extended Bell et al.’s frequency based approach to the multimodal network, by adding services capacity constraints, loading/unloading, and duty check times while empty containers repositioning was not treated. Besides, Cao et al. [9] formulated a stochastic integer programming model for the capacity allocation problem taking into account supply matching with random demands in rail container transportation network, excluding the transport by ship or truck (significantly reducing the complexity of the problem). In addition, they also proposed a large-scale 0-1 integer programming model, subjected to limited container capacities, to maximize the total transportation profit generated by all freight bookings would be accepted in a multiperiod planning horizon [10]. Similarly, Jacko [11] considered the optimal dynamic allocation problem on limited knapsack for perishable or nonperishable items and formulated it in the framework of Markov decision processes.In this paper, we extend the container capacity allocation problem introduced by Cao et al. [9] to a synchromodal transportation network, including rail, river, and road transportation mode, while the empty containers repositioning problem is out of scope. Since the freight demands for container transportation are uncertain in synchromodal transportation network, this paper provides a stochastic integer programming model for the container capacity allocation problem with a single origin and multiple destinations. The aim of the model is to maximize the total transportation profit by allocating container capacity appropriately in planning horizon. Furthermore, a problem-oriented hybrid algorithm with a novel gene encoding method is proposed to address the large-scale stochastic integer programming formulation.The rest of the paper is organized as follows. The container capacity allocation problem with random freight demands in synchromodal transportation network is formulated as a stochastic integer programming model in Section2. According to the characteristics of the proposed model, a problem-oriented hybrid algorithm is designed to solve it in Section 3. And then, we perform some examples to indicate the efficiency of the proposed algorithm in Section 4. Section 5 concludes the paper with a summary of the capacity allocation problem with random demands.


## 2. Mathematical Formulation

This section aims to develop a mathematical model for container capacity allocation problem with random freight demands in synchromodal transportation network.As mentioned above, the container capacity allocation decision is not only a tactical level decision, but also an operational level optimization decision when more freight demands or other pieces of information are available. The overall objective is to determine the limited container capacity allocation, to maximize the potential profit for container carriers. In this paper, we assume that the container capacity allocation plan has to be made with limited knowledge of future scenarios of parameters. Therefore, to incorporate the uncertainty in parameters, these parameters can be realized as one ofS scenarios and the possible realizations probability of scenario s is denoted by Ps, as the scenario-based methods of dealing with random variables in [12].The mathematical model is subjected to several technical constraints. The obvious constraints are that volume of freight should not exceed the total available capacity and the crew, truck, locomotive, and ship at the container central station are assumed sufficient. Additionally, there is at most one container transportation line to each destination for reducing complexity of the loading work. We also assume that the container transport lines of a certain transportation mode, which have not enough demands, would be abandoned.In order to present the mathematical formulation of the container capacity allocation problem clearly, some notations are interpreted in Table1.Table 1
List of notations.
Symbol
Description
(1)Index sets
I
~
Set of freight types1,2,…,i,…,I.
J
~
Set of destinations1,2,…,j,…,J.
M
~
Set of transportation modes1,2,…,m,…,M.
S
~
Set of scenarios1,2,…,s,…,S.
K
~
m
Set of freight transported by modem.
(2)Random data
D
i
j
s
Demands for transportation of freight typei to destination j under scenario s.
α
m
i
j
s
Penalty per unit shortage in supplied capacity for freight typei transported by mode m to destination j under scenario s.
β
m
i
j
s
Penalty per unit overage in supplied capacity for freight typei transported by mode m to destination j under scenario s.
(3)Deterministic parameters
r
m
i
j
Profit for carrying a unit freight of typei to destination j by transportation mode m.
C
m
j
Requested bottom supply of capacity for transportation modem to destination j.
H
Amount of available container capacity at the originO.
(4)Decision variables
x
m
i
j
Container capacity for transportation modem of freight type i to destination j, which will be supplied for shipping cargoes in a planning horizon. In this paper, we assume that the unit of capacity is TEU (i.e., twenty-foot equivalent unit).
y
m
i
j
s
Amount of shortage capacity for freight typei transported by mode m to destination j under scenario s.
z
m
i
j
s
Amount of overage in capacity for freight typei transported by mode m to destination j under scenario s.
u
m
j
Binary variable: that is,umj=1, transportation mode m to destination j will be operated; umj=0, otherwise.As stated in Section1, a container capacity allocation problem is mainly for the carriers who are in the container central station. Containers can be shipped to the place where freight demands occurred, but the available container capacity to each demand is certain after transportation planning is determined. Therefore, there is only one origin in order to simplify. While I freight types will be received at the origin to be shipped to J destinations during a planning horizon and M transportation modes, implying the means of transportation, are available between origin and each destination, besides, S future scenarios indicating the uncertain demands are under consideration. Figure 1 shows the network structure of the container capacity allocation with I freight types, S scenarios, M transportation modes, one origin, and J destinations.Figure 1
Network structure of the container capacity allocation problem.Consequently, the container capacity allocation problem with random freight demands in synchromodal transportation network can be formulated as the following two-stage stochastic integer programming model (P1):(1)MaximizeZ=∑m∈M~∑i∈I~∑j∈J~rmij·xmij-∑s∈S~PsQsxmij(2)subject to∑m∈M~∑i∈I~∑j∈J~xmij≤H,(3)∑i∈I~xmij≥Cmjumj,∀m∈M~,∀j∈J~,(4)∑m∈M~umj≤1,∀j∈J~,(5)umj∈0,1,∀m∈M~,∀j∈J~,(6)xmij∈z+,∀m,i,j∈K~m,where, for all s,(7)Qsxmij=Minimize∑m∈M~∑i∈I~∑j∈J~αmijsymijs+rmijzmijs+βmijszmijs(8)subject toymijs-zmijs=Dijs·umj-xmij,∀i,j∈K~,∀m,i,j∈K~m,s∈S~,(9)ymijs,zmijs∈Z+,s∈S~,∀m,i,j∈K~m,where Z+ denotes the set of nonnegative integers. The objective function of first stage program (1) maximizes the total transportation profit by allocating capacity reasonably in operational planning horizon. The model is subjected to the following constraints. As the available capacity is limited, constraint (2) imposes that the aggregate supply of container capacity cannot exceed the available capacity at the origin O in the planning horizon. Constraint (3) assures the actual container capacity supplied by transportation mode m to destination j cannot be less than a certain bottom value, which may guarantee maximizing the utilization of limited resource. If less, cancel the transportation line of mode m to destination j; that is, umj=0. Constraint (4) assures that there is at most one transportation mode to destination j during the planning horizon to reduce the complexity of loading work. Constraint (5) indicates that the decision variable umj is a binary variable, meaning that transportation mode m to destination j will be operated or not. The objective function of second stage objective function (7) minimizes the overall penalties for overage and shortage in capacity under scenario s, excepted item rmijzmijs. Because zmijs indicates an overage in supplied capacity, it means that the container carriers lose the opportunities for making more transportation profit from the unutilized containers, that is, rmijzmijs. Thus, the total profit generated by carrying freight of type i to destination j by transportation mode m is rmijxmij-rmijzmijs under scenario s, whereas the profit deviation rmijzmijs is moved to (7). Constraint (8) implies the amount of overage or shortage in the supplied capacity under scenario s: for example, umj=1; if Dijs>xmij, it means the supplied capacity is in shortage, then ymijs=Dijs-xmij, and zmijs=0; if Dijs<xmij, it means the supplied capacity is in overage, Zmijs=xmij-Dijs, and ymijs=0; and if Dijs=xmij, it means the supplied capacity is just right, Zmijs=ymijs=0. While umj=0, constraint (8) is meaningless. As the definition of shortage in capacity of ymijs and overage in capacity of zmijs in constraint (9), ymijs=(Dijs·umj-xmij)+, zmijs=(xmij-Dijs·umj)+, where (a)+ denotes maxa,0 for the real number a.Since it is not easy to solve above two-stage stochastic integer program (P1), we reformulate the two-stage programming (P1) into one-stage stochastic integer programming problem (OSIPP) (P2) by integrating objective functions. The objective function of program (P2) is as the following:(10)MaximizeZ=∑m∈M~∑i∈I~∑j∈J~rmij·xmij-∑s∈S~PsQsxmij=∑m∈M~∑i∈I~∑j∈J~rmij·xmij-∑s∈S~Ps∑i∈I~∑j∈J~αmijsDijs·umj-xmij++rmij+βmijsxmij-Dijs·umj+=∑m∈M~∑i∈I~∑j∈J~rmij·xmij-∑s∈S~Ps∑m∈M~∑i∈I~∑j∈J~αmijsymijs+rmijzmijs+βmijszmijs.The constraints in program (P2) consist of constraint (2)–constraint (9), except the second stage objective function (7), for the objective functions of the two-stage objective functions (1) and (7) are integrated in the objective function (10).To the best of our knowledge, there is no effective method to deal with the large-scale stochastic integer programming model (P2) even when there is only one scenario in the problem [12]. Moreover, the large-scale stochastic integer programming model (P2) is NP-hard problem [13]. It is very unlikely to devise a polynomial-time algorithm for it. Before presenting an effective heuristic algorithm, we analyze some special properties of the model (P2) below.Proposition 1.
There is at least one feasible solution for the SIPP (P2).Proof.
Actually,umj=0, xmij=0, ymijs=zmijs=0, ∀(m,i,j)∈K~m, s∈S~ is a feasible solution for (P2), and Z=0.Proposition 2.
The large-scale SIPP (P2) has a finite optimal solution.Proof.
LetR-=max(m,i,j)∈K~mrmij; for any feasible solution of large-scale SIPP (P2), we have(11)Z=∑m∈M~∑i∈I~∑j∈J~rmij·xmij-∑s∈S~Ps∑i∈I~∑j∈J~αmijsymijs+rmijzmijs+βmijszmijs≤∑m∈M~∑i∈I~∑j∈J~rmij·xmij≤R-H.
So, the optimum of objective function of the large-scale SIPP (P2) has a finite upper boundR-H.


## 3. Simulated-Annealing-Based Genetic Algorithms (SAGA) for Capacity Allocation Problem

As stated above, for a NP-hard problem, it is widely accepted that there is no algorithm that can solve the problem in polynomial time. Generally, there are tens of thousands of integer decision variables in the real world capacity allocation problem with random demands for the container carrier and we cannot solve it by any known methods in reasonable time. So, we attempt to design a problem-oriented algorithm to solve the proposed mathematical formulation (P2).The genetic algorithm (GA) was first introduced by Holland in 1975 [14]. It is a kind of global search heuristic used to find exact or high quality solutions for various optimization problems [15–17]. However, simulated annealing (SA), firstly introduced by Kirkpatrick et al. in 1983 to deal with combination optimization problem [18], is effective in local optimum search [19]. At present, the focus of algorithm research has changed from being rather algorithm-oriented to being more problem-oriented [20–23]. Therefore, according to the characteristic of the proposed model, a hybrid algorithm, called simulated-annealing-based genetic algorithms (SAGA) in this paper, is designed for it by uniting the respective advantages of GA and SA. The structure of SAGA is as illustrated in Figure 2.Figure 2
The structure of SAGA.As illustrated in Figure2, the genetic operators are responsible for global research, mainly to produce traditional population A, while the simulated annealing operator is responsible for local research, dedicating to produce elite population B, where Pa is the splitting factor, which governs the size of populations A and B. In early evolution, populations mainly rely on genetic operators for global search. At this time, the size of population A would be larger with smaller Pa. In later period of evolution, populations primarily depend on the individual advantage to search solution near optimal solutions. In this case, the size of population B would be larger with larger Pa. In this paper, the regulation process is divided into two stages: that is, P=0.4Pa when g∈(0,Gmax/2]; P=Pa when g∈(Gmax/2,Gmax], where Pa is the maximum per centum of populations to dominant populations, Gmax is the maximum iterations, and g is the current iterations.


### 3.1. Coding and Decoding

Since most real life optimization problems are with multiple constraints, it is not easy to find a feasible solution. The methodologies discussed in previous literatures mostly relied on some extra penalty functions to revise the infeasible solutions. But, it is not easy to revise the infeasible solutions for most real world problem [13]. As the genes on chromosome are the deciding attributes for individual genetic trait, it is attractive to design a novel method of coding/encoding genes on chromosome without revision of infeasible solutions.Specifically, the genes on the chromosomes are encoded by uniformly distributed random number between 0 and 1 in this paper. Then, the individuals directly participate in crossover, mutation, and simulated annealing operator. Before selection, the genes on each chromosome are decoded into a feasible solution according to the constraints listed in (P2). The decoding process is illustrated in Figure3, where St denotes a certain chromosome t, Ru is the utilization ratio of the total available capacity H, and Xt is the phenotype of chromosome St, indicating a feasible solution. The pseudocode of decoding process is presented in Algorithm 1, where Stl (l∈L~) denotes the lth gene on chromosome St, L~ is a set of discrete locations 1,2,…,l,…,L and L=M×J+1+M×J×I, and M, I, and J represent the maximum number of transportation mode, freight type, and destination, respectively, as stated in Table 1. Cn (n∈N~) denotes the requested bottom supply of capacity in set C~, which is generated according to Cmj; that is, C11=C1, C12=C2,…,CMJ=CN, where N~ is a set 1,2,…,n,…,N, N=M×J. Xtl denotes decision variables umj or xmij according to the specific location on the string of integer; that is, if l<N, Xtl denotes decision variables umj; else if l>N+1, Xtl denotes decision variables xmij.Figure 3
The mapping of decoding.Algorithm 1.
Input parameters:H, St, N, C~, M, I, and J.
Step  1 (handling constraint (2) and (3)). Set total utilization of container capacity V0=ceil(H·RU), a1=St1+St2+⋯+StN. For l=1:N, set b1l=∑i=1i=lSti/a1, AVtl=ceil(b1l·V0); then xtl=AVtl-AVtl-1; if xtl≥Cl, set Xtl=1. Otherwise, Xtl=0, where xtl denotes the total allocating container capacity to a certain destination shipped by a certain transportation mode, corresponding to xtl=∑i∈I~xmij; Xtl corresponds to a certain decision variable umj.
Step  2 (handling constraint (4)). For j=1:J, set Uj=Stj+StJ+j+⋯+St(M-1)J+j, if Uj>1, meaning that there is more than one transportation mode to destination j, which violates constraint (4) in (P2). Thus, set Xtl0=0, l0∈j,J+j,…,(M-1)·J+j. Find a correct location l1=Location_max(Stj,StJ+j,…,St(M-1)·J+j), which denotes the location of the maximum value in set Stj,StJ+j,…,StM-1·J+j on chromosome St; set Xtl1=1.
Step  3 (allocating the container capacity)
Step  3.1. Set l=1, a2=St1·Xt1+St2·Xt2+⋯+StN·XtN. For n=1:N, b2n=∑i=1i=nSti·Xti/a2, AVtn=ceil(b2n·V0), the final container capacity allocated to lth location V0n=AVtn-AVtn-1.
Step  3.2. If V0l>0, go to Step  3.3; else, go to Step  3.4.
Step  3.3. Allocate capacity to each item. Set a3=St(N+1)+(l-1)I+1+St(N+1)+(l-1)I+2+⋯+St(N+1)+lI. For k=1:I, b3k=∑i=1i=kSt(N+1)+(l-1)I+k/a3, AVItk=ceil(b3k·V0l), then Xt(N+1)+(l-1)I+k=AVItk-AVItk-1, where Xt(N+1)+(l-1)I+k corresponds to a certain decision variable xmij.
Step  3.4. Set l=l+1. If l>N, stop; output the results. Otherwise, return to Step  3.2, where b1,2,3l is accumulated allocation proportion, ceil(x) rounds x to the nearest integer greater than or equal to x, AVtl, AVItk represents, respectively, the accumulated allocation capacity to a certain destination shipped by a certain transportation mode and to a certain item, AVt0=0, AVI0=0.According to Algorithm1, all chromosomes can be decoded into feasible solutions satisfying all of the constraints without revision or elimination of invalid solution. For instance, for xt(M-1)∗J+1<C(M-1)∗J+1, u(M-1)·J+1=0; contrarily, for xN≥CN, uN=1. Since the allocated container capacity to destination J by transportation mode M is not empty, that is, uN=1, it needs to further allocate the container capacity to each item; that is, xM1J=Xt(N+1)+(N-1)I+1=3800, xM2J=2700,…,xMIJ=3500, as shown in Figure 3. Additionally, the encoding method can increase the flexibility in crossover, mutation, and simulated annealing operators since there is no concern about whether the potential solutions are feasible or not.Furthermore, it seems quite natural to define its objective value as its fitness function value for maximum optimal problem; that is,f(S)=Z(x). Thus, in this paper, the relationships of search space, solution space, the objective function value, and the value of fitness function are illustrated in Figure 4.Figure 4
The relationship of search space, feasible solution space, the objective function value, and the fitness function value.


### 3.2. Operators

#### 3.2.1. Reproduction and Selection

Up to now, various reproduction methods have been proposed and considered, for example, ranking selection, elitist ranking selection, expected value selection, elitist expected value selection, roulette wheel selection, and elitist roulette wheel selection [24]. In this paper, we further develop the elitism preserve strategy, called developed elitism strategy (DES), and then adopt improved elitism roulette wheel selection to reproduce.In the existing literature, elitism preserve strategy is preserving the elite and directly putting it into the population of following generation. Obviously, it ignores that elite may generate a more excellent individual in crossover operator which also disappears easily in mutation operator. Therefore, the new elitism strategy is improved by preserving the newborn eliteEg_crossover which may be generated in crossover operator and then replacing the worst individual after mutation operator. The DES process is illustrated in Figure 5, where Pg′/Pg+1′, Pg_selection, Pg_crossover, and Pg_mutation present the current gth/(g+1)th populations, the current populations after selection, crossover, and mutation operator, respectively. Pg_mutation′ presents the current populations after replacing the worst individual by the newborn elite, Pg+1 presents the (g+1)th populations after genetic operators, and Eg/Eg+1 represents the elite in the gth/(g+1)th populations. From Figure 5, it can be easily seen that the elites both in initial populations and in crossover populations can not only be preserved but also have chance to produce a preferable individual in DES.Figure 5
The process of DES.


#### 3.2.2. Crossover

In traditional GA, the crossover position is mostly fixed. Employing the same crossover position may be not conducive to evolution, especially in later period of evolution that the crossover area should be larger gradually for individual similarity. For the advantage of code technique, the crossover position can easily be revised in this paper. Therefore, a nonuniform arithmetic crossover with adaptive crossover position is adopted. With two individualsS1g and S2g crossing, the new individuals are obtained as follows:(12)S1g+1=λS2g+1-λS1g,S2g+1=λS1g+1-λS2g.In (12), λ=exp(-λ0·Gmax/g), where λ0 is a constant, λ0=0.9 in this paper, Gmax⁡ is the maximum iterations, and g is the current iterations.


#### 3.2.3. Mutation

For increasing the population diversity, a nonuniform mutation rate is adopted. The individualSig mutate into a new individual according to the following equation:(13)Sig+1=0.95Sig+0.05θ,ifc1<0.5Pm,0.9Sig+0.1θ,ifc1<0.9Pm,θ,ifc1<Pm,where θ is 1×n dimensional matrix formed by uniformly distributed random number between 0 and 1, n is the number of genes on each chromosome, Pm is fixed mutation rate, and c1 is a uniformly distributed random number between 0 and 1. Equation (13) implies that the mutation rate is variable according to the uniformly distributed random number c1.


#### 3.2.4. Simulated Annealing Operator (SAO)

Simulated annealing is set as an operator in simulated-annealing-based genetic algorithms (SAGA). According to the characteristic of simulated annealing, the out layer loop temperature of simulated annealing operator in each generation is gradually decreased with evolution forward; that is,Tg←Tg-1·β, where β is the decay scale and is 0.9 in this paper. Then, set Tg as initial temperature T in simulated annealing operator of current generation g.However, the temperature is lowering byt←t·αIt in simulated annealing operator, where t is the current temperature and It is the current iteration. It is obvious that the cooling process is related not only to the decay scale α, α∈[0,1], but also to the current iteration It. Furthermore, a diverse method of perturbation is adopted for the stability is not the same at different temperatures. The detailed procedure of simulated annealing operator is described in Algorithm 2.Algorithm 2.
Input the out layer loop temperatureTg, Si, L, and Imin.
Step  1 (initialization). Set T=Tg, t=Tg, It=0, and S=Si. Compute fitness f(S) and PreBestS←S, PreBestZ←f(S), BestS←PreBestS, BestZ←PreBestZ, where PreBestS is the preoptimal solution, PreBestZ is the preoptimal the value of object function, BestS is the current optimal solution, and BestZ is the current value of object function.
Step  2. For l=1:L, where L is the Markov length.
Step  2.1 (movements). Perturb the placement through a defined movement and the specific perturbation model is as follows:(14)
S
′
=
0.7
S
+
0.3
θ
t
∈
T
2
,
T
,
0.85
S
+
0.15
θ
t
∈
T
10
,
T
2
,
0.95
S
+
0.05
θ
t
∈
0
,
T
10
,where θ is defined as in Section 3.2.3, S is the previous best solution, and S′ is a newborn solution after perturbing; t is the current temperature in the simulated annealing operator.
Step  2.2 (calculate score). Calculate the change in the score for perturbation. Calculate the increment of objectives by Δf=f(S′)-f(S), where f(S′) is the fitness of the newborn solution S′.
Step  2.3 (Metropolis Criterion). If Δf>0, accept the newborn solution S′, S←S′. Otherwise, accept the newborn S′ with a certain probability, S←PS′, where P=exp(Δf/(kt)), where k is the Boltzmann constant, k=0.85 in this paper.
Step  2.4 (estimate
S
). If f(S)>BestZ, set PreBestS←BestS, PreBestZ←BestZ, BestS←S, BestZ←f(S).
Step  3 (termination conditions). If (BestZ-PreBestZ)/PreBestZ<e (parameter tolerance e = 10−4) or It>Imin, where Imin is the least iterations, stop and output the result. Otherwise, go to Step  4.
Step  4 (lowering temperature). Set It←It+1, t←t·αIt; go to Step  2.


### 3.3. Termination Conditions

When applying the proposed hybrid algorithms to solve the capacity allocation problem (P2), we must obtain an approximate solution of desirable precision within a proper time. For this reason, two parametersGmax and Gmin are introduced, which denote the maximum and the minimum iterative generation, respectively. Then the following termination conditions are imposed:(i)
IfGt>Gmin and (fmax-fmean)/fmax<ε, stop.(ii)
IfGt>Gmax, stop. Otherwise, return to iterative process.ε is a tolerance parameter (ε=0.0005 in this paper), Gt is the current generation, and fmax, fmean are the maximum fitness and the mean fitness in the current populations, respectively.


## 3.1. Coding and Decoding

Since most real life optimization problems are with multiple constraints, it is not easy to find a feasible solution. The methodologies discussed in previous literatures mostly relied on some extra penalty functions to revise the infeasible solutions. But, it is not easy to revise the infeasible solutions for most real world problem [13]. As the genes on chromosome are the deciding attributes for individual genetic trait, it is attractive to design a novel method of coding/encoding genes on chromosome without revision of infeasible solutions.Specifically, the genes on the chromosomes are encoded by uniformly distributed random number between 0 and 1 in this paper. Then, the individuals directly participate in crossover, mutation, and simulated annealing operator. Before selection, the genes on each chromosome are decoded into a feasible solution according to the constraints listed in (P2). The decoding process is illustrated in Figure3, where St denotes a certain chromosome t, Ru is the utilization ratio of the total available capacity H, and Xt is the phenotype of chromosome St, indicating a feasible solution. The pseudocode of decoding process is presented in Algorithm 1, where Stl (l∈L~) denotes the lth gene on chromosome St, L~ is a set of discrete locations 1,2,…,l,…,L and L=M×J+1+M×J×I, and M, I, and J represent the maximum number of transportation mode, freight type, and destination, respectively, as stated in Table 1. Cn (n∈N~) denotes the requested bottom supply of capacity in set C~, which is generated according to Cmj; that is, C11=C1, C12=C2,…,CMJ=CN, where N~ is a set 1,2,…,n,…,N, N=M×J. Xtl denotes decision variables umj or xmij according to the specific location on the string of integer; that is, if l<N, Xtl denotes decision variables umj; else if l>N+1, Xtl denotes decision variables xmij.Figure 3
The mapping of decoding.Algorithm 1.
Input parameters:H, St, N, C~, M, I, and J.
Step  1 (handling constraint (2) and (3)). Set total utilization of container capacity V0=ceil(H·RU), a1=St1+St2+⋯+StN. For l=1:N, set b1l=∑i=1i=lSti/a1, AVtl=ceil(b1l·V0); then xtl=AVtl-AVtl-1; if xtl≥Cl, set Xtl=1. Otherwise, Xtl=0, where xtl denotes the total allocating container capacity to a certain destination shipped by a certain transportation mode, corresponding to xtl=∑i∈I~xmij; Xtl corresponds to a certain decision variable umj.
Step  2 (handling constraint (4)). For j=1:J, set Uj=Stj+StJ+j+⋯+St(M-1)J+j, if Uj>1, meaning that there is more than one transportation mode to destination j, which violates constraint (4) in (P2). Thus, set Xtl0=0, l0∈j,J+j,…,(M-1)·J+j. Find a correct location l1=Location_max(Stj,StJ+j,…,St(M-1)·J+j), which denotes the location of the maximum value in set Stj,StJ+j,…,StM-1·J+j on chromosome St; set Xtl1=1.
Step  3 (allocating the container capacity)
Step  3.1. Set l=1, a2=St1·Xt1+St2·Xt2+⋯+StN·XtN. For n=1:N, b2n=∑i=1i=nSti·Xti/a2, AVtn=ceil(b2n·V0), the final container capacity allocated to lth location V0n=AVtn-AVtn-1.
Step  3.2. If V0l>0, go to Step  3.3; else, go to Step  3.4.
Step  3.3. Allocate capacity to each item. Set a3=St(N+1)+(l-1)I+1+St(N+1)+(l-1)I+2+⋯+St(N+1)+lI. For k=1:I, b3k=∑i=1i=kSt(N+1)+(l-1)I+k/a3, AVItk=ceil(b3k·V0l), then Xt(N+1)+(l-1)I+k=AVItk-AVItk-1, where Xt(N+1)+(l-1)I+k corresponds to a certain decision variable xmij.
Step  3.4. Set l=l+1. If l>N, stop; output the results. Otherwise, return to Step  3.2, where b1,2,3l is accumulated allocation proportion, ceil(x) rounds x to the nearest integer greater than or equal to x, AVtl, AVItk represents, respectively, the accumulated allocation capacity to a certain destination shipped by a certain transportation mode and to a certain item, AVt0=0, AVI0=0.According to Algorithm1, all chromosomes can be decoded into feasible solutions satisfying all of the constraints without revision or elimination of invalid solution. For instance, for xt(M-1)∗J+1<C(M-1)∗J+1, u(M-1)·J+1=0; contrarily, for xN≥CN, uN=1. Since the allocated container capacity to destination J by transportation mode M is not empty, that is, uN=1, it needs to further allocate the container capacity to each item; that is, xM1J=Xt(N+1)+(N-1)I+1=3800, xM2J=2700,…,xMIJ=3500, as shown in Figure 3. Additionally, the encoding method can increase the flexibility in crossover, mutation, and simulated annealing operators since there is no concern about whether the potential solutions are feasible or not.Furthermore, it seems quite natural to define its objective value as its fitness function value for maximum optimal problem; that is,f(S)=Z(x). Thus, in this paper, the relationships of search space, solution space, the objective function value, and the value of fitness function are illustrated in Figure 4.Figure 4
The relationship of search space, feasible solution space, the objective function value, and the fitness function value.


## 3.2. Operators

### 3.2.1. Reproduction and Selection

Up to now, various reproduction methods have been proposed and considered, for example, ranking selection, elitist ranking selection, expected value selection, elitist expected value selection, roulette wheel selection, and elitist roulette wheel selection [24]. In this paper, we further develop the elitism preserve strategy, called developed elitism strategy (DES), and then adopt improved elitism roulette wheel selection to reproduce.In the existing literature, elitism preserve strategy is preserving the elite and directly putting it into the population of following generation. Obviously, it ignores that elite may generate a more excellent individual in crossover operator which also disappears easily in mutation operator. Therefore, the new elitism strategy is improved by preserving the newborn eliteEg_crossover which may be generated in crossover operator and then replacing the worst individual after mutation operator. The DES process is illustrated in Figure 5, where Pg′/Pg+1′, Pg_selection, Pg_crossover, and Pg_mutation present the current gth/(g+1)th populations, the current populations after selection, crossover, and mutation operator, respectively. Pg_mutation′ presents the current populations after replacing the worst individual by the newborn elite, Pg+1 presents the (g+1)th populations after genetic operators, and Eg/Eg+1 represents the elite in the gth/(g+1)th populations. From Figure 5, it can be easily seen that the elites both in initial populations and in crossover populations can not only be preserved but also have chance to produce a preferable individual in DES.Figure 5
The process of DES.


### 3.2.2. Crossover

In traditional GA, the crossover position is mostly fixed. Employing the same crossover position may be not conducive to evolution, especially in later period of evolution that the crossover area should be larger gradually for individual similarity. For the advantage of code technique, the crossover position can easily be revised in this paper. Therefore, a nonuniform arithmetic crossover with adaptive crossover position is adopted. With two individualsS1g and S2g crossing, the new individuals are obtained as follows:(12)S1g+1=λS2g+1-λS1g,S2g+1=λS1g+1-λS2g.In (12), λ=exp(-λ0·Gmax/g), where λ0 is a constant, λ0=0.9 in this paper, Gmax⁡ is the maximum iterations, and g is the current iterations.


### 3.2.3. Mutation

For increasing the population diversity, a nonuniform mutation rate is adopted. The individualSig mutate into a new individual according to the following equation:(13)Sig+1=0.95Sig+0.05θ,ifc1<0.5Pm,0.9Sig+0.1θ,ifc1<0.9Pm,θ,ifc1<Pm,where θ is 1×n dimensional matrix formed by uniformly distributed random number between 0 and 1, n is the number of genes on each chromosome, Pm is fixed mutation rate, and c1 is a uniformly distributed random number between 0 and 1. Equation (13) implies that the mutation rate is variable according to the uniformly distributed random number c1.


### 3.2.4. Simulated Annealing Operator (SAO)

Simulated annealing is set as an operator in simulated-annealing-based genetic algorithms (SAGA). According to the characteristic of simulated annealing, the out layer loop temperature of simulated annealing operator in each generation is gradually decreased with evolution forward; that is,Tg←Tg-1·β, where β is the decay scale and is 0.9 in this paper. Then, set Tg as initial temperature T in simulated annealing operator of current generation g.However, the temperature is lowering byt←t·αIt in simulated annealing operator, where t is the current temperature and It is the current iteration. It is obvious that the cooling process is related not only to the decay scale α, α∈[0,1], but also to the current iteration It. Furthermore, a diverse method of perturbation is adopted for the stability is not the same at different temperatures. The detailed procedure of simulated annealing operator is described in Algorithm 2.Algorithm 2.
Input the out layer loop temperatureTg, Si, L, and Imin.
Step  1 (initialization). Set T=Tg, t=Tg, It=0, and S=Si. Compute fitness f(S) and PreBestS←S, PreBestZ←f(S), BestS←PreBestS, BestZ←PreBestZ, where PreBestS is the preoptimal solution, PreBestZ is the preoptimal the value of object function, BestS is the current optimal solution, and BestZ is the current value of object function.
Step  2. For l=1:L, where L is the Markov length.
Step  2.1 (movements). Perturb the placement through a defined movement and the specific perturbation model is as follows:(14)
S
′
=
0.7
S
+
0.3
θ
t
∈
T
2
,
T
,
0.85
S
+
0.15
θ
t
∈
T
10
,
T
2
,
0.95
S
+
0.05
θ
t
∈
0
,
T
10
,where θ is defined as in Section 3.2.3, S is the previous best solution, and S′ is a newborn solution after perturbing; t is the current temperature in the simulated annealing operator.
Step  2.2 (calculate score). Calculate the change in the score for perturbation. Calculate the increment of objectives by Δf=f(S′)-f(S), where f(S′) is the fitness of the newborn solution S′.
Step  2.3 (Metropolis Criterion). If Δf>0, accept the newborn solution S′, S←S′. Otherwise, accept the newborn S′ with a certain probability, S←PS′, where P=exp(Δf/(kt)), where k is the Boltzmann constant, k=0.85 in this paper.
Step  2.4 (estimate
S
). If f(S)>BestZ, set PreBestS←BestS, PreBestZ←BestZ, BestS←S, BestZ←f(S).
Step  3 (termination conditions). If (BestZ-PreBestZ)/PreBestZ<e (parameter tolerance e = 10−4) or It>Imin, where Imin is the least iterations, stop and output the result. Otherwise, go to Step  4.
Step  4 (lowering temperature). Set It←It+1, t←t·αIt; go to Step  2.


## 3.2.1. Reproduction and Selection

Up to now, various reproduction methods have been proposed and considered, for example, ranking selection, elitist ranking selection, expected value selection, elitist expected value selection, roulette wheel selection, and elitist roulette wheel selection [24]. In this paper, we further develop the elitism preserve strategy, called developed elitism strategy (DES), and then adopt improved elitism roulette wheel selection to reproduce.In the existing literature, elitism preserve strategy is preserving the elite and directly putting it into the population of following generation. Obviously, it ignores that elite may generate a more excellent individual in crossover operator which also disappears easily in mutation operator. Therefore, the new elitism strategy is improved by preserving the newborn eliteEg_crossover which may be generated in crossover operator and then replacing the worst individual after mutation operator. The DES process is illustrated in Figure 5, where Pg′/Pg+1′, Pg_selection, Pg_crossover, and Pg_mutation present the current gth/(g+1)th populations, the current populations after selection, crossover, and mutation operator, respectively. Pg_mutation′ presents the current populations after replacing the worst individual by the newborn elite, Pg+1 presents the (g+1)th populations after genetic operators, and Eg/Eg+1 represents the elite in the gth/(g+1)th populations. From Figure 5, it can be easily seen that the elites both in initial populations and in crossover populations can not only be preserved but also have chance to produce a preferable individual in DES.Figure 5
The process of DES.


## 3.2.2. Crossover

In traditional GA, the crossover position is mostly fixed. Employing the same crossover position may be not conducive to evolution, especially in later period of evolution that the crossover area should be larger gradually for individual similarity. For the advantage of code technique, the crossover position can easily be revised in this paper. Therefore, a nonuniform arithmetic crossover with adaptive crossover position is adopted. With two individualsS1g and S2g crossing, the new individuals are obtained as follows:(12)S1g+1=λS2g+1-λS1g,S2g+1=λS1g+1-λS2g.In (12), λ=exp(-λ0·Gmax/g), where λ0 is a constant, λ0=0.9 in this paper, Gmax⁡ is the maximum iterations, and g is the current iterations.


## 3.2.3. Mutation

For increasing the population diversity, a nonuniform mutation rate is adopted. The individualSig mutate into a new individual according to the following equation:(13)Sig+1=0.95Sig+0.05θ,ifc1<0.5Pm,0.9Sig+0.1θ,ifc1<0.9Pm,θ,ifc1<Pm,where θ is 1×n dimensional matrix formed by uniformly distributed random number between 0 and 1, n is the number of genes on each chromosome, Pm is fixed mutation rate, and c1 is a uniformly distributed random number between 0 and 1. Equation (13) implies that the mutation rate is variable according to the uniformly distributed random number c1.


## 3.2.4. Simulated Annealing Operator (SAO)

Simulated annealing is set as an operator in simulated-annealing-based genetic algorithms (SAGA). According to the characteristic of simulated annealing, the out layer loop temperature of simulated annealing operator in each generation is gradually decreased with evolution forward; that is,Tg←Tg-1·β, where β is the decay scale and is 0.9 in this paper. Then, set Tg as initial temperature T in simulated annealing operator of current generation g.However, the temperature is lowering byt←t·αIt in simulated annealing operator, where t is the current temperature and It is the current iteration. It is obvious that the cooling process is related not only to the decay scale α, α∈[0,1], but also to the current iteration It. Furthermore, a diverse method of perturbation is adopted for the stability is not the same at different temperatures. The detailed procedure of simulated annealing operator is described in Algorithm 2.Algorithm 2.
Input the out layer loop temperatureTg, Si, L, and Imin.
Step  1 (initialization). Set T=Tg, t=Tg, It=0, and S=Si. Compute fitness f(S) and PreBestS←S, PreBestZ←f(S), BestS←PreBestS, BestZ←PreBestZ, where PreBestS is the preoptimal solution, PreBestZ is the preoptimal the value of object function, BestS is the current optimal solution, and BestZ is the current value of object function.
Step  2. For l=1:L, where L is the Markov length.
Step  2.1 (movements). Perturb the placement through a defined movement and the specific perturbation model is as follows:(14)
S
′
=
0.7
S
+
0.3
θ
t
∈
T
2
,
T
,
0.85
S
+
0.15
θ
t
∈
T
10
,
T
2
,
0.95
S
+
0.05
θ
t
∈
0
,
T
10
,where θ is defined as in Section 3.2.3, S is the previous best solution, and S′ is a newborn solution after perturbing; t is the current temperature in the simulated annealing operator.
Step  2.2 (calculate score). Calculate the change in the score for perturbation. Calculate the increment of objectives by Δf=f(S′)-f(S), where f(S′) is the fitness of the newborn solution S′.
Step  2.3 (Metropolis Criterion). If Δf>0, accept the newborn solution S′, S←S′. Otherwise, accept the newborn S′ with a certain probability, S←PS′, where P=exp(Δf/(kt)), where k is the Boltzmann constant, k=0.85 in this paper.
Step  2.4 (estimate
S
). If f(S)>BestZ, set PreBestS←BestS, PreBestZ←BestZ, BestS←S, BestZ←f(S).
Step  3 (termination conditions). If (BestZ-PreBestZ)/PreBestZ<e (parameter tolerance e = 10−4) or It>Imin, where Imin is the least iterations, stop and output the result. Otherwise, go to Step  4.
Step  4 (lowering temperature). Set It←It+1, t←t·αIt; go to Step  2.


## 3.3. Termination Conditions

When applying the proposed hybrid algorithms to solve the capacity allocation problem (P2), we must obtain an approximate solution of desirable precision within a proper time. For this reason, two parametersGmax and Gmin are introduced, which denote the maximum and the minimum iterative generation, respectively. Then the following termination conditions are imposed:(i)
IfGt>Gmin and (fmax-fmean)/fmax<ε, stop.(ii)
IfGt>Gmax, stop. Otherwise, return to iterative process.ε is a tolerance parameter (ε=0.0005 in this paper), Gt is the current generation, and fmax, fmean are the maximum fitness and the mean fitness in the current populations, respectively.


## 4. Numerical Experiments

Some numerical experiments are carried out to demonstrate the applications of the mathematical formulation and efficiency of the proposed algorithm. The algorithm has been coded in C++ and run under Microsoft Windows 7 Professional (Intel Core 2 Duo CPU 2.26 GHz and 2.0 GB of RAM). CPU times were obtained through the C++ function clock().Example 1.
As shown in Table2, we firstly suppose a small-scale instance of the optimization capacity allocation problem for container carrier in synchromodal transportation only with one origin: Shanghai and 3 destinations: Dalian, Zhengzhou, and Chongqing, which are container central stations. In the case, there are 3 transportation modes: River, Rail, and Road, 4 freight types: Electronic products, Mechanical equipment, Pharmaceuticals, and Textiles, and 2 scenarios, that is, M=3, I=4, J=3, and S=2 with 227 variables, 86 constraints. Since shortage in container capacity means that carriers lose the chance to make transportation profit, the penalty coefficient of shortage is much more than that of overage, almost equaling unit delivery profit. On the other side, the overage is relatively less as more penalty is included in rmijzmijs which indicates that the container carriers lose the opportunities for making more transportation profit from the unutilized containers zmijs. These penalty coefficients are consistent with [9]. The relevant experimental data are tabulated in Tables 3–5, including capacity and demands information, capacity and demands information, profit and penalties information.
Before verifying the effectiveness and efficiency of the proposed SAGA, the actual optimal solution of the above numerical example is obtained through the branch and bound method (B-and-B) in GLPK 3.2.2 (see the website of free software GLPK, 2002).
u
22
=
1
,
u
31
=
1, u33=1, x212=3800, x222=4500, x132=3200, x142=1900, x311=2400, x321=3200, x331=2500, x341=4400, x313=3800, x323=2700, x333=3400, and x343=3500. Other decision variables equal zero and the actual optimal objective value opt_obj=19428550.
Firstly, we test the sensitivity and feasibility of the proposed algorithm with respect to different parameters settings. Randomly generate 10 sets of parameters in SAGA, and for each set of experiments, the above small-scale instances were performed 20 times. The computational results are listed in Table6, where Pc, Pm, Psize, Gmin, and Gmax represent the population size, crossover probability, mutation probability, and maximum and minimum iterative generation in genetic operators, respectively; T0, α, L, and Imin represent the initial temperature, the factor of step length, Markov length, and the minimum iteration in simulated annealing operator, respectively. Besides, Best, Average, and Time represent the approximate optimal objectives, the average objectives, and the average processing time in 20 times, respectively. The relative gap between best approximate optimal objectives in each set of experiments and exact optimal objectives is denoted by δ and calculated by the following formulation:(15)δ=opt_obj-Bestopt_obj×100%.
It is easy to see from these results that the total average objectives and variance are 18901431 and 2.71%, respectively, which illustrates the robustness and steadiness of the proposed solution methodology. Additionally, it is interesting to note that the processing time is sensitive to the setting of parameters and varies in the numerical experiments for different parameters. For instance, the average processing time in seventh sets is just 39.1 s since relatively smaller parameters are adopted. In comparison, the tenth sets cost the longest processing time due to their much larger parameter settings. Actually, most of processing time in each experiment is wasted and useless. With the parameters in sets 7, the obtained best objective in this set of experiments is 19119450 and the relative gap is only 1.59%. Therefore, with consideration of processing time and objectives precision, the parameters in sets 7 are adopted in the rest of experiments; that is,Pc=0.1, Pm=0.8, Psize=50, Gmin=50, Gmax=70, Pa=0.12, T0=3000, α=0.9, L=120, and Imin=10.Table 2
Basic information.
The number of scenarios in a planning horizon: 2The container central station of origin: Shanghai
Freight type
Probability of each scenario
Destination
(1) Electronic products
(1) 1/2
(1) Dalian
(2) Mechanical equipment
(2) 1/2
(2) Zhengzhou
(3) Pharmaceuticals
(3) Chongqing
(4) TextilesTable 3
Capacity and demands information.
Freight type
Total available capacity (TEU) in the planning horizon: 40000
Dalian
Zhengzhou
Chongqing
Scenario 1
Scenario 2
Scenario 1
Scenario 2
Scenario 1
Scenario 2
1
2400
2200
3800
3500
3800
3200
2
3200
3100
4500
4100
2700
2300
3
2500
2300
3200
2900
3400
3100
4
4400
4200
1900
1500
3500
3000
Transportation mode
The requested bottom supply of capacity for each transportation mode to each destination
Road
3000
4000
3000
Rail
4000
5000
4000
River
3000
3000
4000Table 4
Profit (CNY¥) information.
Freight type
The unit delivery profit for each freight type shipped to each destination by each transportation mode
Dalian
Zhengzhou
Chongqing
Road
Rail
River
Average
Road
Rail
River
Average
Road
Rail
River
Average
1
500
490
580
523
440
510
460
470
460
490
480
477
2
550
540
600
563
450
520
440
470
490
470
490
483
3
530
490
570
530
470
530
430
477
480
500
510
497
4
440
480
560
493
450
470
450
457
500
450
460
470Table 5
Penalties (CNY¥) for per overage and shortage in supplied capacity under each scenario.
Transportation mode
Scenario
Dalian
Zhengzhou
Chongqing
1
2
3
4
1
2
3
4
1
2
3
4
Penalty per overage in supplied capacity for delivery to each destination by each mode under each scenario
Road
1
86
74
78
74
87
87
76
67
75
54
66
86
2
76
87
70
61
86
59
87
80
55
72
83
79
Rail
1
87
68
85
83
72
81
57
59
82
61
65
70
2
82
90
85
56
81
60
65
78
54
82
71
67
River
1
89
78
68
88
83
73
80
75
66
73
82
66
2
81
80
88
87
77
81
59
64
67
67
74
73
Penalty per shortage in supplied capacity for delivery to each destination by each mode under each scenario
Road
1
496
449
600
523
578
587
487
494
511
592
599
466
2
562
584
577
600
442
526
553
587
514
597
507
463
Rail
1
590
504
482
542
519
520
592
520
571
506
465
505
2
501
469
450
583
454
469
444
534
576
556
438
477
River
1
571
515
438
594
480
484
590
585
583
485
502
507
2
528
498
450
459
571
550
472
489
496
595
597
468Table 6
The computational results.
Sets
P
c
P
m
P
size
G
m
i
n
G
m
a
x
P
a
T
0
α
L
I
m
i
n
Time (×101)
Best (×101)
Average (×107)
δ
1
0.3
0.7
120
60
80
0.09
5000
0.7
90
12
7.64
1.9200
1.9100
1.18%
2
0.1
0.9
50
30
50
0.13
2000
0.9
100
10
2.36
1.9088
1.8837
1.75%
3
0.2
0.8
110
60
70
0.12
5000
0.9
110
12
11.07
1.9138
1.9031
1.50%
4
0.3
0.7
60
30
50
0.14
5000
0.9
80
13
2.63
1.9020
1.8655
2.10%
5
0.2
0.7
120
60
80
0.11
5000
0.8
110
10
8.96
1.9127
1.9082
1.55%
6
0.1
0.8
80
40
60
0.14
3000
0.9
100
14
6.58
1.9133
1.8772
1.52%
7
0.1
0.8
50
50
70
0.12
3000
0.9
120
10
3.91
1.9119
1.8991
1.59%
8
0.1
0.8
110
30
50
0.07
4000
0.8
120
12
3.59
1.9092
1.8718
1.73%
9
0.4
0.9
80
60
80
0.14
5000
0.7
80
15
8.31
1.9160
1.8923
1.38%
10
0.4
0.8
80
60
80
0.14
4000
0.8
110
15
11.31
1.9164
1.9010
1.36%
Average
1.9124
1.8901Example 2.
In this example, some experiments are implemented to demonstrate the efficiency of SAGA by comparing the results with these obtained by the traditional genetic algorithm (GA) and simple simulated annealing algorithm (SA). Besides, the parameter values used in SAGA were set as the same as mentioned in Example1. In consideration of fairness, Pc, Pm, T0, and α are set as the same in SAGA, while Psize, Gmax in GA and L, Imin in SA are adjusted to try to guarantee that the processing time of the three algorithms is more or less the same. In general, these parameters were adopted in GA: Pc=0.1, Pm=0.8, Psize=500, and G=500, while T0=3000, α=0.9, L=600, and Imin=300 in SA. Additionally, the code/decode techniques and other relative parameters are the same as SAGA.
For the small-problem in Example1, 20 times were performed through SAGA, GA, and SA, respectively. The results and processing time of SAGA, GA, and SA are shown in Figure 6. One of the evolutionary performances of these algorithms is illustrated in Figure 7.
It is easily seen from Figure6 that the results obtained through GA are relatively stable, while there are big fluctuations both in objectives and in processing time through SA. Moreover, in SA, the longer processing time does not imply more outstanding results. Obviously, SA has the initial solution-dependence disease during the evolutionary procession. It means that if the initial optimum has high quality, it can be quickly close to global optimal solutions. From Figure 7(a), we can find that GA can quickly find a relatively optimal solution in about 20 iterations. Unfortunately, it falls into the local optimal after 220 iterations. However, from Figure 7(b), it is easily seen that the simulated annealing procession can hardly produce a better solution after 45 accept points, which implies that SA is limited to local optimal research and incapable to break through local scope. It is worthwhile to note that it just needs 5 iterations to find a relatively optimal solution through SAGA, which is even better than that obtained through GA in its whole procession, as shown in Figure 7(c). Meanwhile, after 20 iterations, the objective is better than 1.9 × 107. By contrasting Figures 7(a)–7(c), we can draw a conclusion that SAGA has significant advantages, such as higher efficiency and shorter processing time since it takes full advantages of the remarkable characteristic of GA in global search and that of SA in local search. The comparative results are shown in Table 7.
For the small-scale instance, it can be seen from Table7 that the relative gap and the average processing time through SAGA are much smaller than through GA and SA. Although the minimum relative gap obtained through SA is 0.07%, which is less than 1.59% obtained through SAGA, the average processing time is much more than SAGA. Moreover, the worst objective is much smaller than SAGA. Referring to Figure 6, it is easily seen that SA is quite unsteady. However, the quantity of results obtained through GA is in between SA and SAGA, but it may not be unsatisfactory. Meanwhile, compared with B-and-B, the average processing time is smaller, and the relative gap is within the allowed relative gap range.Table 7
Experimental results for the small problem: 227 variables and 86 constraints (20 times).
Methods
Best (×107)
Average(×107)
Worst (×107)
Time (×101 s)
δ
GA
1.8831
1.8417
1.7672
9.07
3.08%
SA
1.9293
1.1727
0.5576
16.88
0.07%
SAGA
1.9119
1.8991
1.8842
3.91
1.59%
B-and-B
1.9429
13.58Figure 6
The results and processing time of SAGA, GA, and SA.
(a)
(b)Figure 7
Evolutionary performances of these algorithms, where (a) is the optimal solution of each generation through GA; (b) is the optimal solution of each accepted point through SA; (c) is the optimal solution of each generation through SAGA.
(a)
(b)
(c)Example 3.
In this example, we consider a larger-scale instance of the optimization capacity allocation problem.
In the medium-scale instance:M=3, I=8, J=5, and S=9, with 2338 variables and 1105 constraints. It is difficult to input real data manually to the model with tens of thousands of parameters and variables. Therefore, we randomly generated an instance from uniform distribution to test our algorithm for solution quality and computational time. Let [x] denote the largest integer that is less than or equal to x. In the section, Dijs, H, and Cmj are uniformly generated in the range [1800,4500], 1800×I×J,4500×I×J, [1800×I,4500×I], respectively. The parameters rmij, αmijs, and βmijs are randomly generated from uniform distribution in the range [450,650], [400,550], [60,90], respectively. And then, SAGA, GA, SA, and B-and-B for the instance are performed 20 times, respectively. The results are tabulated in Table 8, where “#” represents the exact optimal solution which cannot be obtained through the B-and-B method. Thus, there is no relative gap between the exact optimal solution and approximate optimal solutions obtained through SAGA, GA, and SA, indicated by “—” in the table.
In Table8, it proves once more that SA is very unsteady, regardless of the results or processing time. However, the results through GA are relatively stable but unsatisfactory. Contrarily, SAGA has obvious advantages in average objectives and processing time, showing promising computational results.
From what has been discussed in Examples1–3, we can conclude that the numerical solutions of the proposed SAGA are very close to the exact solution, and the relative gaps are all in the acceptable relative gap range. The main reason is that the proposed oriented-problem algorithm combined the GA with SA, so that it is good not only at local search but also in global search. The excellent individuals may have more opportunity to produce better individuals by simulated annealing operator. Meanwhile, it can overcome the unstable characteristic of SA by genetic operators. To sum up, the proposed algorithm is an effective algorithm for the model (P2) of the container capacity allocation problem with random demands for the container carrier in synchromodal transportation.Table 8
Experimental results for large-scale problem: 2338 variables and 1105 constraints (20 times).
Methods
Best(×107)
Average (×107)
Worst(×107)
Time(×101 s)
δ
GA
5.1493
5.0848
4.9471
13.42
—
SA
5.3107
3.7473
1.4473
20.93
—
SAGA
5.4416
5.4158
5.3831
7.58
—
B-and-B
#
#
# represents the exact optimal solution which cannot be obtained through the B-and-B method.


## 5. Conclusions

On the basis of analyzing the research status of capacity allocation in synchromodal transportation, this paper presents a stochastic integer programming model for the container capacity allocation problem with random freight demands in synchromodal transportation network from container carriers’ perspective. In this paper, the purpose of the proposed model is to maximize the expected total transportation profit over optimal decisions. Since there is no effective polynomial-time algorithm to solve the proposed formulation in the real world, a problem-oriented intelligence algorithm is designed for the proposed formulation while some numerical examples are carried out to demonstrate the effectiveness and efficiency of the hybrid algorithm.As the above analysis in Section4, four types of algorithms, including GA, SA, SAGA, and B-and-B, are tested on two instances, respectively. The proposed SAGA, by integrating simulated annealing method with genetic algorithm, shows promising computational results. The SAGA with a novel encoding method and a developed elitism strategy, comparing to simple GA and SA, has obvious advantage in processing time and results. To sum up, with the problem-oriented algorithm, the proposed programming formulation has important practical significance for the container carrier.For the container capacity allocation problem in this research without considering the repositioning of empty containers, one can consider the capacity allocation problem by taking the repositioning of empty containers into account to minimizing the operational cost in future research.



---

*Source: 986152-2015-11-03.xml*