MetaHeuristic Algorithms for a MultiProduct Dynamic LotSizing Problem with a Freight Container Cost
 Author: Kim Byung Soo, Lee WoonSeek
 Organization: Kim Byung Soo; Lee WoonSeek
 Publish: Industrial Engineering and Management Systems Volume 11, Issue3, p288~298, 30 Sep 2012

ABSTRACT
Lot sizing and shipment scheduling are two interrelated decisions made by a manufacturing plant and a thirdparty logistics distribution center. This paper analyzes a dynamic inbound ordering problem and shipment problem with a freight container cost, in which the order size of multiple products and single container type are simultaneously considered. In the problem, each ordered product placed in a period is immediately shipped by some freight containers in the period, and the total freight cost is proportional to the number of containers employed. It is assumed that the load size of each product is equal and backlogging is not allowed. The objective of this study is to simultaneously determine the lotsizes and the shipment schedule that minimize the total costs, which consist of production cost, inventory holding cost, and freight cost. Because the problem is NPhard, we propose three metaheuristic algorithms: a simulated annealing algorithm, a genetic algorithm, and a new populationbased evolutionary metaheuristic called selfevolution algorithm. The performance of the metaheuristic algorithms is compared with a local search heuristic proposed by the previous paper in terms of the average deviation from the optimal solution in small size problems and the average deviation from the best one among the replications of the metaheuristic algorithms in large size problems.

KEYWORD
Dynamic LotSizing , Shipment Scheduling , MultiProducts , MetaHeuristics

1. INTRODUCTION
Most manufacturing companies do not have a specialty of logistics functions. In order to obtain the benefit of shipping and warehousing costs, they strategically ally with a specialized thirdparty logistics (TPL) company for offering a wide range of services including: order fulfillment, inbound freight, warehousing, freight consolidation, outbound dispatching. They outsource the logistics operations to TPL company and concentrate on their core manufacturing operation.
Among the services provided by a TPL company, this paper focuses on the most important process is replenishing the multiple products by a number of freight containers at right time from manufacturers to the TPL warehouse. This process leads to managerial decision problems in determining replenishment lotsizes and shipment schedule for each demand. For these problems the freight cost as well as the ordering and holding cost are considered. Since the replenishment orders are assumed to be shipped in containers to the warehouse, the freight cost is assumed to be proportional to the number of containers used. To solve this problem, this paper proposes an integrated optimization model of inbound ordering and outbound dispatching for single product in a TPL company.
The demand in this paper is assumed to be known, but dynamic over a discrete and finite time horizon. The productioninventory problem with the demand is called the dynamic lotsizing model (DLSM), which has stemmed from the work of Wagner and Whitin (1958). Though there are a number of research results dealing with the DLSM, the majority of DLSMs have not considered any productioninventory problem incorporating shipment problem as well. These days, the issue of shipment scheduling for ordered products (or delivering orders) by a proper shipping freight container mode at right time becomes significantly important in production and distribution management because of the increasing fuel price.
Several studies have focused on various general costs and the capacitated resources as the extended works of the classical DLSM. Lippman (1969) studied two deterministic multiperiod production planning models; monotone cost model and concave model. Florian
et al . (1980) and Bitran and Yanasse (1982) proved the general capacitated single product lot sizing problem is NPhard. Hwang and Sohn (1985) dealt with a DLSM in which the transportation mode and the order size for a deteriorating product are simultaneously considered. However, they considered no capacity restriction on the transportation mode. Lee (1989) considered a DLSM allowing multiple setup costs including a fixed charge cost and a freight cost, where a fixed single container type with limited carrying capacity is considered and the freight cost is proportional to the number of containers used. Fumero and Vercellis (1999) proposed an integrated optimization model for production and distribution planning considering such operational decisions as capacity management, inventory allocation, and vehicle routing. The solution of the integrated optimization model was obtained using the Lagrangian relaxation technique. Leeet al . (2003) extended the works of Lee (1989) by considering multiple heterogeneous vehicle types to immediately transport the finished product in the same period it is produced. It is also assumed that each vehicle has a typedependent carrying capacity and the unit freight cost for each vehicle type is dependent on the carrying capacity. Leeet al . (2003) considered a dynamic model for inventory lotsizing and outbound shipment scheduling in the thirdparty warehousing domain. They presented a polynomial time algorithm for computing the optimal solution. Jaruphongsaet al . (2005) analyzed a dynamic lotsizing model in which replenishment orders may be delivered by multiple shipment modes with different lead times and cost functions.Anily and Tzur (2005) considered a dynamic model of shipping multiple items by capacitated vehicles. They presented an algorithm based on a dynamic programming approach. Van Norden and van de Velde (2005) dealt with a multiple product problem of determining transportation lotsizes in which the transportation cost function has piecewise linear as to a transportation capacity reservation contract. They proposed a Lagrangian relaxation algorithm to compute lower and upper bounds. Lee
et al . (2005) proposed a heuristic algorithm for a dynamic lotsizing and shipping problem, in which the order size of multiple products and a single container type are simultaneously considered. Kim and Lee (2012) are proposed a tight lowerbound using the shortest reformulation model compared to the solution of the upperbound of the heuristic algorithm proposed by Leeet al . (2005). The limitation of their studies is that the performance is not successive in large size problems because the proposed solution approach is a local search heuristic.In this paper, we propose three metaheuristic algorithms to define the multiproduct lotsize problem and shipment scheduling that minimize the total costs, which consist of ordering cost, inventory holding cost, and freight cost by generalizing the basic model of Kim and Lee (2012). The proposed heuristic algorithms are the first study in simultaneously determining multiproduct lotsizing problem and the shipment scheduling.
This paper is organized as follows. In the next section, the mathematical model of the problem is described. In Section 3, the properties of the optimal solution are characterized. Three metaheuristic algorithms based on the optimal solution properties are described in Section 4. The computational results from a set of simulation experiments are presented in Section 5.
2. PROBLEM DESCRIPTION
The following notations are defined to formulate the problem:
Parameters TH = the set of discrete time horizon, (t = 1, 2, …, T)
PD = the set of products, (i = 1, 2, …, M)
dti = amount demanded for product i in period t,
wi = volume of product i,
W = carrying capacity of a container,
St = ordering cost of product i in period t,
hti = unit inventory holding cost of product i from period t to period t+1,
F = unit freight cost of container,
Variables xti = amount of product i ordered and shipped by container in period t,
zti = 1, if an order is incurred for product i in period t, and 0, otherwise,
yt = number of containers used in period t (nonnegative integer),
lti = inventory level of product i in period t
The objective of the problem is to determine (
x_{ti} ,y_{t} ) for ∀i ∈TH andt ∈TH so that all demands over the given horizon are satisfied at the minimum total cost. Therefore, theT period problem (P ) can be formulated in a mathematical programming as follows:Subject to If we assume that all products have the same volume, the model (
P ) can be specialized to (P1 ) proposed by Kim and Lee (2012).Subject to and (1), (3), (4), (5), and (6)
Constraint (2) ensures that the inventory level of each product in the current period balances the inventory level in the previous period, the ordering amount and the demand in the current period. Constraints (3) and (7) imply that the total ordering amount is restricted by the total carrying capacity associated with the number of containers used in the period. Constraint (4) ensures the relation between
x_{ti} andz_{ti} . The constraints in the above model (P ) and (P1 ) define a closed bounded convex set and the objective function is concave, so that the problem attains its minimum at an extreme point of the convex set. The extreme points will be further characterized in association with the optimal solution in the next section.3. OPTIMAL SOLUTION PROPERTIES
The mathematical model (
P ) can be represented by a network model as Figure 1. In the network, two flow types are defined as follows. 1) The aggregate flow is defined as the flow between node 0 and nodes (t = 1, 2, …,T ). 2) The individual flow is defined as the flow between nodes (t = 1, 2, …,T ) and nodes ((1, 1), (1, 2), …, (T ,M )).Here, the arcs in the aggregate flow are restricted by the capacities associated with the number of containers used whereas the arcs in the individual flow are not.
The optimal solution of the model (
P ) occurs at extreme points. In a network theory, such an extreme point can be interpreted as an extreme flow (Florinet al ., 1980; Zangwill, 1969). In a network without arc capacities, a feasible flow is an extreme flow if it does not have a positive loop. Also, in a network with arc capacities, a feasible flow is an extreme flow if and only if each loop has at least one saturated arc.In Figure 1, loops can be formed by two cases as follows. 1) Between the aggregate flow and the individual flow, for example, the loop can be formed by the sequences of nodes (0), (1), (1, 1), (2, 1), (2) and (0); 2) On the individual flows, for example, the loop can be formed by the sequences of nodes (1, 1), (1), (1, M), (2, M), (2), (2, 1), and (1, 1).
Ordering and shipment schedule satisfying the optimal solution property by Wagner and Whitin (1958),
x_{ti} ？I_{ti} = 0, is no longer the extreme flow. Such policies may form positive loops between the aggregate flows. The properties of Theorems 1 and 2 must be satisfied so as to have an extreme flow. To examine the properties, the production point, the partial aggregate point, and the inventory point are defined, respectively as follows. 1) Periodt is a ordering point for producti ifx_{ti} ≥ 0 ; 2) Container typej is a partial aggregate point in periodt if(
n is a nonnegative integer); 3) Periodt is an inventory point for product i ifI_{ti} = 0.Theorem 1: In the model P, the optimal solution has at most one partial aggregate point for product i between two consecutive inventory points for product i .Proof: Suppose that there exists the optimal solution which has two partial aggregate point between two consecutive inventory points for producti . In the network of Figure 1, this case can have the loop formed by the sequences of nodes (0), (1), (1, 1), (2, 1), (2), and (0). The condition that this feasible flow will be the extreme flow is that at least one of the arcs (0, 1) and (0, 2) must be saturated. This feasible flow is not an extreme flow. Therefore, the proof is completed.Theorem 2: In the model P, the optimal flow must do not form the positive loop in the individual flow .Proof: Suppose that there exists a feasible flow that satisfies the properties of Theorems 1 and 2, which has a loop formed by the sequences of nodes (1, 1), (1), (1, 2), (2, 2), (2), (2, 3), (3, 3), (3), (3, 1), (2, 1) and (1, 1) in Figure 1. Because arcs in the individual flow are not capacitated, there exists a positive loop formed by unsaturated arcs. This feasible flow is not an extreme flow. Therefore, the proof is completed.Van Norden and van de Velde (2005) proved that the single level noncapacitated multiproduct multiperiod lotsizing with transportation cost model is NPhard in the strong sense. Hence, the problem (
P ) is NPhard. So, it is not easy to make optimization viable for large problems. In the next section, three metaheuristic algorithms are presented based on the properties of Theorems 1 and 2.4. METAHEURISTIC ALGORITHMS
In this section, we propose three metaheuristics: a simulated annealing (SA) algorithm, a genetic algorithm (GA), and a selfevolution algorithm (SEA) based on the properties of Theorems 1 and 2.
4.1 Solution Representation
For the solution representation of three metaheuristic algorithms, a solution is defined as two dimensional array of (
T ,M ) of 01 genes as follows:In this representation,
e_{ti} can take a value 0, 1, 2 which is the key to find the value ofx_{ti} . In this situation, the value ofe_{ti} directly decides the value ofx_{ti} . Then depending variable setsz_{ti} ,l_{ti} , andy_{t} can be calculated, once the value ofx_{ti} is determined by the following procedure:Step 1: Seti = 0.Step 2: Letp be the earliest period of producti aftert that has nonzeroe_{ti} . If all the period aftert have zero values fore_{pi} , letp beT +1Step 3: SetStep 4: Ife_{ti} = 1,x_{ti} =X andz_{ti} = 1. On the other hands, else ife_{ti} = 2, /x_{ti} =W ？X /W ？ andz_{ti} = 1, otherwisex_{ti} =z_{ti} = 0Step 5: Ifi =M , determinethen STOP. Otherwise,
i =i +1Property 1: For the modelP1 , the variable setx_{ti} , are determined, the other decision variable setsz_{ti} ,l_{ti} , andy_{t} can be determined.In this representation, one can easily find that the first nonzero value of
e_{ti} ,i ∈PD must be no later than the first positive value ofd_{ti} to meet the constraint (4), i.e., the constraint of nobacklogging. Thus, a repairing process is required after genetic operations of GA and SEA and perturbations of SA.4.2 Genetic Algorithm
The GA, which has been widely used in solving various optimization problems for the last three decades, is a stochastic search algorithm based on the mechanism of natural selection and natural genetics. Being different from conventional search techniques, it starts with an initial set of (random) solutions called a population. Each individual in the population is called a chromosome, representing a solution to the problem at hand. The chromosomes are evaluated, using some measures of fitness. Generally speaking, the GA is applied to spaces that are too large to be exhaustively searched (Goldberg, 1989). In this paper, the chromosomes that have higher fitness value (lower objective function value) than the average fitness of current population make a potential parent pool and new chromosomes in the next generation are randomly selected from the parent pool or selected from the roulette wheel method.
Two kinds of crossover operators are randomly used: partiallymatched crossover (PMX) and uniform crossover (UX). Eight kinds of mutation operators are
randomly used: the eight kinds of mutation operators (pull operator, swap operator, insert operator, inner random operator, outer random operator, uniform random operator, single change operator, and rest change operator) to make a new chromosome. For the operators, one, two or multiple points in the selected original chromosome are randomly selected and mutated. The pull operator is illustrated in Figure 2a. The genes on right side of point 2 (including point 2) are pulled to the position of point 1, and the genes between point 1 and 2 (including point 1) are placed after. The two genes at the points are interchanged for swap operator as shown in Figure 2b. Insert operator simply insert the gene at point 2 into the position of point 1 as shown in Figure 2c. Inner random operator and outer random operator are illustrated in Figure 2d and e. The inner or outer genes of point 1 and 2 are randomly replaced for the operators. Uniform random operator randomly select multiple points and change as shown in Figure 2f. Figure 2g and h describe single change and other change in simple way. Using the crossover operators, the mutation operators, and the selection operators, the selected parents reproduce new chromosomes (i.e., children) to generate a population for the next generation. The various parameters of the GA heuristic are summarized in Table 1. These parameters were selected based on extensive preliminary experimentations to the best combination with highest performance. GA evaluates a total of 1,000×(
T +M +N ) of fitness function.The following procedure outlines the pseudo code of the GA:
Procedure GA for multiproduct lotsizing problem with a freight container
Begin
Set crossover rate and mutation rate as rc = 0.8 and rm = 0.2
Set a generation size as g = 0
Generate an initial population Pc
Evaluate f ( ？ ) all chromosomes in Pc
Find the best chromosome scbest in Pc
Let Pbest = scbest and fbest = f(scbest)
Repeat (the generation loop)
Repeat (the population loop)
Randomly select two parents sc1 and sc2 among the population Pc using roulette wheel selection
Generate two new children sn1 and sn2 by following genetic operations:
If (random [0, 1] < rc )
Then perform randomly PMX crossover or UX crossover
If (random [0, 1] < rm )
Then perform randomly a mutation among 8 mutations
Until (number of populations reaches a maximum of T+M+N)
Evaluate f ( ？ ) in Pn
Update Pc by selecting an elitist and the rest of chromosomes in the population using randomly roulette wheel selection or random selection within a good solution pool from Pn
If fbest > f (scbest)
Then
Pbest > scbest and fbest > f (scbest) in the population Pc
Until (number of generation reaches a maximum of 1,000 times)
End
4.3 Simulated Annealing
SA was first proposed by Kirkpatrick
et al . (1983) to solve combinatorial optimization problems. SA explores the solution space by successive moves from one potential solution to another that is a variation of its predecessor. At the start of the process, an initial solution is generated and evaluated based on the total cost (fitness) in Section 2. SA differs from local search algorithms in that the procedure uses random selection and will sometimes accept nonimproving moves (interchanges) hoping to expand the search space and ultimately reach a better overall solution. The nonimproving moves are probabilistically performed using the Boltzmann probability mass function as follows (Wolsey, 1998):where
t is the current temperature, Δf =f (P_{n} ) ？f (P_{c} ) andp_{n} andp_{c} are the candidate fitness and the current fitness values, respectively.In the implementation of this paper, the SA heuristic starts with a randomly generated chromosome and is controlled using two loops: an outer loop and an inner loop. The outer loop controls the temperature reduction according to the cooling schedule. In the inner loop, the temperature is held constant and a predetermined number of explorations are made. In the outer loop, explorations are made at each temperature. For fair comparison with other heuristics, the number of inner loops is 50 and the number of outer loops is 20, and the number of perturbations consists of (
T +M +N ) using random selection of the eight mutations in Figure 2 to increase the chances of obtaining a better solution. This leads a total of 1,000×(T +M +N ) evaluations. This high number of perturbations was motivated by Barbarosoglu and Ozdamar (2000) who indicated that the increase in the number of search moves is carried out by SA significantly improves its performance. The various parameters of the SA heuristic are summarized in Table 2. These parameters were selected based on extensive preliminary experimentations to determine the best combination that leads to the highest frequency of hitting the optimal solution.The following outlines the outlines pseudo code of the SA algorithm:
Procedure SA for multiproduct lotsizing problem with a freight container
Begin
t = 0
θ0 = 1,000
Randomly select a search point Pc
Evaluate f(Pc)
Let Pbest = Pc and fbest = f(Pc)
Repeat (the outer cooling loop)
Repeat (the inner cooling loop)
Repeat
Apply a specific perturbation on Pc to produce Pn
Evaluate f(Pn) with a local search heuristic
If f(Pn) < f(Pc)
Then
Pc = Pn
If fbest > f(Pn)
Then
Pbest > Pn and fbest > f(Pc)
Else if random [0, 1] < e(f (Pc )？ f (Pn )) / θt
Then
Pc = Pn
Until ((T+M+N) random trials are exhausted among 8 perturbations)
Until (number of iterations reaches a maximum of 50 times)
t = t+1
Until (t reaches a maximum of 20 times)
End
4.4 SelfEvolution Algorithm
SEA is a metaheuristic algorithm which has a population (a set of solutions) based mechanism that uses the evolution of a solution by itself (selfevolution). Similar to GA, the set of chromosomes forms a population. Initial population is generated randomly, and the chromosomes in the population are evaluated by the measure of the fitness introduced in Section 2. A chromosome from the population is randomly selected and executes a selfreproduction using a random selection of one of mutation operators described in Section 4.2. Then the new chromosome is evaluated and it replaces the original chromosome, if the fitness value of the new chromosome is better than that of the original chromosome. The algorithm continues until the number of selfreproductions becomes a predetermined stopping value.
We propose the operators used in the eight mutations in GA. SEA is running without providing any parameters for the algorithm, because all the selection processes in SEA, such as selection of chromosome from the population for selfevolution, selection of evolution operator, and selection of points for the operator are randomly executed. SEA showed good performance for a machine scheduling problem compared with any other metaheuristics (Joo and Kim, 2012).
For fair comparison with other heuristics, the number of generations is fixed as 1,000 for terminating, and the number of evolution operators consists of (
T +M +N ) using random selection of the eight mutations in Figure 2 to increase the chances of obtaining a better solution. This leads a total of 1,000×(T +M +N ) evaluations with the same number of evaluations of GA and SA. The various parameters of the SEA heuristic are summarized in Table 3. These parameters were selected based on extensive preliminary experimentations to determine the best combination that leads to the highest frequency of hitting the optimal solution.The following procedure outlines the pseudo code of the SEA:
Procedure SEA for multiproduct lotsizing problem with a freight container
Begin
Set a generation size as g = 0
Generate an initial population Pc
Evaluate f(Pc)
Find the best solution scbest in Pc
Let sbest = scbest and fbest = f(scbest)
Repeat (the generation loop)
Randomly select a solution sc among the population Pc using roulette wheel selection
Apply a specific perturbation sc in the incumbent population Pc to produce sn for the population Pn of the next generation using randomly selected the perturbation operations a mutation among 8 mutations
Evaluate f(sn)
If f(sn) < f(sc)
Then
sc = sn
If fbest = f(sn)
Then
Pbest = sn and fbest = f(sn)
Until (number of generation reaches a maximum of 1,000×(T+M+N) times)
End
5. COMPUTATIONAL RESULTS
To analyze the performance of the proposed heuristic algorithm, the following experimental conditions were designed.
Set M = 3, 6, 10 and T = 4 and 8 in small sized problems and T = 15, 18, and 24 in large sized problems
Demands were generated from a normal distribution,
Mean μi was generated from an uniform distribution, U(25, 100).
Standard deviation σi was equally likely selected from μi or μi/5.
Setup cost was selected by
and TSi = 1, 3, 6, where TSi denotes EOQ time supply.
hti = 1 was assumed without loss of generality.
Set wi were generated from an uniform distribution, U(0, 5)
W = 100, 200, 300 and respective unit freight cost was proportionally selected by W as follows:
F =i ？ W, i =1, 3, 6.
To evaluate the performance of the heuristic, the C# computer code for the proposed heuristic was run on an Intel Core™2 Duo CPU with 2.00 GHz RAM. Also, CPLEX package for finding the optimal solution was run on the same computer and run so as to find the best solution within 700,000 node limits.
For the performance of the heuristics, the relative percentage deviation (RPD) calculated with the expression (10) and the computing time by GA, SA, and SEA of 10 replications for each test problem.
Where
z_{B} = objective value of the best solution by CPLEX in small sized test or the objective value of the best solution among the solutions of replications of GA, SA, and SEA andz_{H} = the average objective value of the heuristic.In small sized problems, four instances were made for each combination of input parameters and 10 replicates for each instance are performed in the metaheuristics (GA, SA, and SEA). In Table 4, we compare SEA shown in the best performance among SA and GA with the local search heuristic (Heu) proposed by Kim and Lee (2012). In the second column in the Table 4, SEA shown in the best performance offers on the no more than 5% from the optimal solution with small computa
tional time for small sized test problems. The difference of RPDs between SEA and Heu increases as
T (the size of planning horizon) andM (the number of items) increases, because Heu is easy to converge local optima as the problem size increases. The third and fourth columns in Table 4 show the average computing times taken by CPLEX and SEA that are measured in seconds. From the table, it is shown that SEA offers a good solution recorded in 0.34 seconds in average sense. Meanwhile, the computing time of CPLEX dramatically increases asT (the size of planning horizon) andM (the number of items).In large sized problems, the absolute performance between CPLEX and the best solution of the heuristics (ARPD) is presented with the expression (11) and the relative performance of the heuristics (SA, GA, and SEA) is compared by the best value of the heuristics are shown in Table 5.
Where
z_{HB} = the best heuristic solution among SA, GA, and SEA andz_{Opt} = the optimal solution using CPLEX without node limits.Table 5 shows that the optimal solution using CPLEX was not obtained within 1 hour for many largesized test problems having more than 15 periods as shown in ‘NA’ because of the limitations of a computer performance. This result indicates that it is difficult to find the optimal solution as the problem size increases (
T andM ). For the relative comparison, SEA offers the best performance by showing 1.40% of RPD value in average sense compared to SA and GA. In order to validate the results, it is important to check if the observed differences in the RPD values of each metaheuristic algorithm are statistically significant. Figure 3 shows the mean plots and Tukey honest significant difference (HSD) intervals at the 95% confidence level of all problems in Table 5. The superiority of SEA can be explained by the graphs in the graph in Figure 3a. The RPD value of SEA is consistent with small variance asT (the size of planning horizon) andM (the number of items) increase, but it increases asW (container capacity) increases. The computing time of SEA presents no more than 5 seconds.6. CONCLUSION
In this paper, a dynamic lot sizing problem and shipment scheduling are simultaneously analyzed, where the order size of multiple products and a single container type are simultaneously considered. To address the problem, two different solution approaches are proposed. The first approach is based on a mixed integer programming model. Because the problem is NPhard, three metaheuristic algorithms are proposed to increase solution efficiency based on the optimal solution properties. SEA is a new metaheuristic algorithm which has a population (a set of solutions) based selfevolution mechanism. Two problem groups are tested to verify the performance of proposed metaheuristic algorithms. SEA offers the best performance by showing 1.40% of RPD value in average sense compared to SA and GA. Also, the RPD value of SEA is consistent with small variance as
T (the size of planning horizon) andM (the number of items) increase, but it increases asW (container capacity) increases. Overall, the test results indicate that SEA is very effective and efficient algorithm with low variation for the dynamic inbound ordering problem in calculating the transportation cost.Further study is required to assess the performance of SEA with other metaheuristics (tabusearch and antcolony optimization, etc.) in solving the inbound ordering problem.

[Figure 1.] Network representation of the model P1.

[Figure 2.] Mutation operators of genetic algorithm: (a) pull, (b) swap, (c) insert, (d) inner random, (e) outer random, (f) uniform random, (g) single change, and (h) rest change operator.

[Table 1.] Parameters of the genetic algorithm heuristics

[Table 2.] Parameters of the simulated annealing heuristics

[Table 3.] Parameters of the selfevolution algorithm heuristics

[Table 4.] RPD and CPU times in small sized problems

[Table 5.] Absolute RPDs between CPLEX and SEA and relative RPDs among SA, GA, and SEA in large sized problems

[Figure 3.] Mean plots and Tukey honest significant difference intervals at the 95% confidence level: (a) algorithm, (b) period (T), (c) products (M), and (d) container capacity (W). GA: genetic algorithm, SA: simulated annealing, SEA: selfevolution algorithm, RPD: relative percentage deviation.