1. INTRODUCTION
The world fleet expanded to more than nine billion dollars during 2012, reaching more than 1.5 billion deadweight tons (DWT) in 2012, an increase of over 37 percent in just four years. Container transport is regarded as a crucial part of the world's truly global supply chain. Numerous factors should be taken into account for ship route design, one of the most important of which is port service accessibility. In addition, timeline design plays a role in air pollution due to its impact on ship fuel consumption (Lei, 2022). Containers are moved by container shipping companies over scheduled routes. Several products such as manufactured goods, food, clothing and apparel are shipped in containers. Generally, shipping services operate within a fixed schedule and has an appropriate port arrangement, in a way that arrival and departure time at ports are specified and planned. Port schedules are posted on the APL website, and customers can set the delivery time of their shipments based on the dates announced on the website and determine the time required for reaching the port (Salido et al., 2011;Alharbi et al., 2014;Chang et al., 2010;Qi and Song, 2012). Therefore, container shipping is significant part of the world’s global supply chain.
The network of container shipping lines includes many shipping lines that must be assessed for the route of each ship. Scheduling shipping lines is a planned decision at technical level made every three to six months. Port service accessibility is the first factor considered for ship route design, without which the ship route schedule design is not possible. The rotating time decreases when the ship’s speed is higher in container shipping lines, which leads to a lower number of ships required for transport. The present study aims to develop a model for scheduling shipping lines in each network of container shipping lines in order to minimize the overall ship costs, including the costs of charging ships and the cost of maintaining the inventory of empty containers in the port by considering the time window of the port and the amount of fuel. The remainder of the article is structured, as follows: Section 2 presents the research background, and Section 3 addresses the proposed mathematical model. Section 4 explains the solution method, and Section 5 focuses on the generation of initial solution. Section 6 presents computational results and section seven concludes and makes suggestions for future studies.
2. RESEARCH BACKGROUND
There are limited studies on the design of ship route schedule design. The first group of studies is related to designing a schedule at the technical planning level. In a research, Wang and Meng (2012) focused on a tacticallevel liner ship route schedule design in a container shipping company with a high number of ports and different routes and a fixed sailing speed on each voyage leg (Wang et al., 2014). In another research, X and Song presented an optimal vessel schedule to minimize the total expected fuel consumption. In addition, the rotation time was randomly considered, meaning that ships arrived at a port of call no later than the announced time (Trappey et al., 2011). Wang and Meng (2012) considered a schedule for the ships existing on the route when uncertainty in port operations and scheduling were considered while assuming that the ships could make up for the delay by sailing across the ocean (Yan et al., 2009). In their previous research, Wang and Meng (2012) improved uncertainty at sea and ports by incorporating speed optimization. These scholars adopted a dynamic planning approach to a practical tactical liner ship route schedule design problem with a port time window. They assumed that each port on a ship route is visited only once in a roundtrip journey. In reality, however, there are many routes where a port can be visited twice (Yin et al., 2011). Wang et al. (2015) developed their previous research with the assumption that each port on the route of the ship can be visited twice. In addition, they focused on the route of one ship instead of a container network with several shipping routes (Sun et al., 2012).
Yan et al. (2009) developed a container routing model aiming at maximizing operating profit in order to set a schedule at operational level. In another research, Brouer et al. (2013) presented the Vessel Schedule Recovery Problem (VSRP) to evaluate a given disruption scenario and to select a recovery action balancing the tradeoff between increased fuel consumption and the impact on cargo in the remaining network and the customer service level. Other studies, such as those developed by Chang et al. (2010, 2011), He et al. (2010), Du et al. (2011), have focused on port operations. Sun et al. (2012) created a general simulation platform, which aims to provide a flexible modeling system. Yan et al. (2009) introduced a distributed agent system for dynamic port planning and scheduling and examined their hypothesis with a case study. Salido et al. (2011) developed a planning technique for solving the container stacking problem and a set of optimized allocation algorithms for solving the berth allocation problem independently. Moreover, He et al. (2010) postulated a strategy for resolve the issue of sharing internal porters across multiple container terminals and selected a nearoptimal solution using an optimization method.
Alharbi et al. (2014) examined a practical liner shipping schedule design problem with port time windows for container supply chain networks. They proposed a mixedinteger nonlinear nonconvex model that incorporates the availability of ports to minimize the sum of ship cost and fuel cost while considering the port time window. In the end, the model was reformulated as an integer linear optimization model and was solved applying an iterative optimization approach.
3. PROPOSED MATHEMATICAL MODEL
In order to better design the mathematical model, the symbols, variables and parameters are defined first, followed by presenting the objective function and related constraints. Moreover, the necessary explanations of the details of the mathematical model are provided below.
3.1 Symbols and Variables
The symbols and variables used in the mathematical model are presented below:
R shows the series of shipping routes and P indicates the set of ports, in a way that the route of a ship r ∈ R is shown, as follows:
P_{ri} ∈ P shows the port corresponding to the ith port of call, and N_{r} is the number of ports existing on the route of the ship. The set of ports of call on the r route are defined, as follows:

t_{rij} : the time (day) required for a ship to move from the ith port to the jth port on the r route.

W: total days of the week, W = {0,1,2,3,4,5,6}, where zero is indicative of Saturday, one shows Sunday, and…

${t}_{ri}^{arr}$ : time of arrival at the ith port of call on the r route, ${t}_{ri}^{arr}\in W$.

${t}_{r,{N}_{r}+1}^{arr}$ : the time (day) required by a ship to return to the first port of call on the r route.

m_{r} : the number of ships existing on the r route.

${x}_{ri}^{h}$ : the number of empty containers reserved at the ith port on the r route.

${z}_{ri}^{w}:\{\begin{array}{c}1\\ 0\end{array}$ : a binary variable; 1, if the ship arrives at the ith port of call on the r route in the w day of the week. Otherwise, 0.

${z}_{ri}^{bw}:\{\begin{array}{c}1\\ 0\end{array}$: a binary variable; 1, if the ship arrives at the bth harbor of the ith port of call on the r route in the w day of the week. Otherwise, 0.
In addition, fuel consumption (tons/miles) is considered as a function shown below:
3.2 Parameters
Several parameters are used in the proposed mathematical model, as shown below:

α : Ship fuel prices (unit/ton)

${\delta}_{b}^{w}$ : 1, if the b ∈ B_{p} port is free in the w day; otherwise, 0.

β and a_{ri} : coefficients in the fuel consumption function (a_{ri} > 0)

B_{p} : the set of docks of the P port

${C}_{r}^{s}$ : the costs of charging ships existing on the r route

${C}_{r}^{k}$ : the cost of maintaining the inventory of empty containers at the ith port in the r route

L_{ri} : the distance between the ith port and the jth port in the r route

${t}_{ri}^{port}$ : the time passed by a ship at the ith route in the r route

${t}_{ri}^{min}$ : the shortest time expected for shipping to the ith port in the r route

${m}_{r}^{max}$ : the maximum number of ships existing in the r route

${x}_{ri}^{hmax}$ : the maximum number of empty containers reserved at the ith port in the r route

${v}_{r}^{max}$ : the maximum ship speed in the r route

${z}^{+}$ : the set of positive integers
According to the mentioned symbols, parameters and variables, the mathematical model proposed encompasses a fiveexpression objective function and 14 sets of constraints, as described below:
Accordingly, the objective is to minimize the total costs of recharging ships, costs of inventory of empty containers in the port, and costs of fuel.
Constraints (5) guarantee solution symmetry. Constraints (6) indicates that ship cannot travel at speeds in excess of the speed limit. Constraints (7) defines the relationship between a time of a round trip. Constraints (8) shows the time of return to the first port of call on the return path. Constraints (9) demonstrates that the number of ships does not exceed the maximum number of ships and is also a positive integer. The time of arrival at each port of call is a nonnegative integer shown in Constraints (11). The second set of constraints show the time of week for arrival at each port of call in the network:
The day of arrival at the port of call on the route is shown by constraints (12), whereas constraints (13) demonstrate that each ship must arrive at a port of call in each week. Constraints (14) define the binary variable, and the k_{ri} auxiliary variable is defined as a nonnegative integer. The third set of constraints points out the accessibility of ports:
According to Constraints (16), an available port cannot serve more than one ship in a day. Moreover, Constraints (17) and (18) indicate that each ship can use each port of call within a week.
4. THE PROPOSED GA SOLUTION METHOD
Since the problem in the present study is an NPhard problem, it is impossible to use precise methods for solving the problem in reasonable time. Therefore, heuristic or metaheuristic methods should be used to solve the problem. In this research, the genetics algorithm (GA) is applied to solve the problem. In addition, GAMS and GA are used to exhibit the solvability of the model at different dimensions. The GA implementation stages and the solution process are explained below.
Algorithm 1. Pseudocode of the genetic algorithm

1: t →0;

2: Forming an initial population [P(t)];

3: [P(t)]; population member assessment and sorting

4: to be performed when the termination conditions are not realized

5: P'(t); formation of new solutions

6: P'(t); evaluation of new solutions

7: $\left[P\text{'}(t){{\displaystyle \cup}}^{\text{}}\hspace{0.17em}Q\right]$ applying genetics operators → P(t + 1); the next generation population

8: T ← t + 1

9: end (until)
5. INITIAL SOLUTION GENERATION
An initial population of chromosomes must be generated after determining the coding system and specifying the method of turning each solution into chromosome. In most cases, the initial population is created randomly. However, heuristic methods are sometimes used to increase algorithm quality and speed to generate the initial population. The size of the initial population often depends on the other coded string. For instance, the initial population selected must be definitely larger when there are 32bit chromosomes in the problem, compared to 16 bit chromosomes. The cutoff point is often between 80 to 95%, whereas the probability of mutation is considered between 0.5 and 1 percent and the population size regarded in the range of 2030. Based on the fit function, the selected chromosomes are allocated a real amount, which shows their value, and the GA stages continue. If there is a very low number of chromosomes, the GA will be less able to perform less combination operations, and only a small part of the search space will be discovered. On the other hand, the process of the GA algorithm will be slow if there is a very high number of chromosomes. According to the results, the use of the population will not be extremely effective as a result of some constraints that mostly depend on coding and the problem itself.
5.1 Coding
This is probably the most difficult problem solution stage of the algorithm method. Binary coding is one of the simplest coding techniques and the best conversion for genetic operators. In binary conversion, the members of the population become strings of zeros and ones. Figure 1 shows the chromosome coding structure.
5.2 Population
In the GA, the concept of population is similar to what actually exists in normal life. As the first stage of the GA, it is required to generate a set of doable solutions as the initial population. The members of this set are often selected randomly. However, some conditions are used in the optimization algorithms so that there is not an excessively scattered population. In addition, the number of population members depend on the type of problem. In fact, the number of members is a parameter that can be changed by improving the accuracy of solutions and search convergence speed. While an eightmember population might be completely suitable for some problems, a population with more than 100 members might not be sufficient for other problems. Therefore, it is best to specify the population within the range of 10160 members.
5.3 Population Size
Equation (19) is used to calculate the population size:
For instance, if each chromosome has a length of 25, then we will have: ${N}_{pop}=\frac{1}{65}\times {2}^{(0/21\times 25)}$, or if the length of each chromosome is equal to 25, the population size will be equal to 270.
5.4 Roulette Wheel Selection
The roulette wheel is one of the most suitable random selection methods, in which the probability of selection corresponding to each chromosome is estimated based on its fitness. In this regard, if fk shows the fit value of the kth chromosome, the probability of survival corresponding to that chromosome will be as shown in Equation (20):
Afterwards, the chromosomes are sorted based on P_{k}, and q_{k}, which is the cumulative values of P_{k}, are obtained using Equation (21):
In fact, the roulette wheel creates a random number between zero and one for selection of each chromosome, and the chromosome corresponding to the mentioned number is selected whatever the range of the number is.
5.5 Mutation
In nature, some factors, such as ultraviolet radiation, cause unpredictable changes in chromosomes. Since GAs follow the law of evolution, these algorithms use the lowprobability mutation operator. The mutation causes a search in the intact spaces of the problem, and the most important responsibility of mutation is to avoid convergence to local optimization. Figure 2 shows mutation and its function.
In general, mutation prevents the GA from locating itself in local extremes. Each member is a value determined by the user that depends on the mutation possibility (P_{m}).
5.6 GA Parameter Tuning
The efficiency of metaheuristic algorithms is directly related to tuning its parameters so that the correct choice of parameter values increases the efficiency of the algorithm. The improper tuning of parameters could lead to improper solutions for the problem under study. While various statistical methods have been proposed for designing test, increasing the number of factors studied, complex and timeconsuming calculations are performed in using a comprehensive approach such as complete factor experiments. Taguchi introduced a set of fractional factor experiments, which specifically reduce the number of tests required for assessment while maintaining the required information (1). In this project, control factors of Taguchi method include parameters of GA, population size and number of iterations, probability of crossover performance and probability of mutation performance. Three levels are used for analyzing the optimal values of parameters in Taguchi method, as shown in Table 1:
For the factors considered from the standard tables of orthogonal arrays, the most suitable array is L9, which is presented in Table 2. The integer values in each column shows its level in each test. The goal of the method is to find optimal levels of significant and controllable factors and minimize the effect of confounding factors. The qualitative features of the values measured from the tests are turned into signal to noise (S/N). This rate shows the level of deviations in the solution variable. In this study, the lower the value of objective function, the better. Therefore, the (S/N) ratio has the feature of the larger the better. In the Taguchi method, this index is defined as Equation (22):
Five sample problems are considered for the implementation of the tests. Each of the nine different tests designed in the orthogonal array are implemented five times for each problem. Figure 3 shows the value of (S/N) for different GA levels. As observed, decrease of algorithm deviations occurs when the parameters of the problem are tuned for the algorithm, as follows: firstlevel crossover rate, thirdlevel mutation rate, combination of the number of initial population and thirdlevel generation, minimum value of secondlevel external weight and combination of number of particles and thirdlevel generation.
6. COMPUTATIONAL RESULTS
After tuning the optimal levels of GA parameters, the GA results can be compared with exact solutions of GAMS 24.1.2 in terms of the value of objective function and solution time using coding in MATLAB. Equation 23 is applied to measure the relative gap between GAMS and MATLAB solutions:
Table 3 shows the input values of the most important parameters of the problem model:
According to the presented model and its solution methods, some problems are generated as example and the problem is solved for each example. The results obtained from comparing the solution of 24 example problems with the results obtained from GA are presented in Table 4. A comparison of GA results with the optimal solution of GAMS 24.1.2 shows that the software is able to achieve an optimal solution in solving largescale problems based on GA. On the other hand, the solution time of problems in two sections of GA and GAMS 24.1.2 indicates less solution time and is associated with a very slight increase in large scales.
As observed in Figure 4, GAMS is unable to solve all problems. meanwhile, the GA could solve all the problems in reasonable time. There is a significant difference between GAMS and GA in terms of solution time with an increase in the problem dimension, which shows that the inability of GAMS to operate when there is an increase in the scale of the problem.
According to Figure 5, the error of the metaheuristic method is very small among the problems that have provided the ability to compare the GA with GAMS software, in a way that the maximum error rate is less than 1%. This confirms the high efficiency of the GA regarding the optimization of the mathematical model proposed in the current research.
7. CONCLUSION
This study modeled a scheduling problem for a container supply chain in order to minimize costs of recharging ships and costs of inventory of empty containers at the port, which is actually a planned decision technical level. The results of the present study could be used for realworld problems with a slight change. The results showed the effectiveness of displacement inside ports and fuel price on the overall costs, the optimal number of ships used, and the optimal scheduling table. The model presented in the current study was solved by GA for GAMS 24.1.2 at large scales. According to Table 4, the error percentage of the GA objective function was less than two percent in all solved problems, compared to the software, which demonstrated the efficiency of GA. It is not possible to solve the model at a reasonable time with available processing facilities for scales larger than 28 ports of call, which justified the use of the metaheuristic approach. According to Table 4, the time required for problem solving with GAMS 24.1.2 was extremely higher, compared to the GA. Therefore, with regard to the solutions obtained and very short time required for problem solving, the GA was reported to have high efficiency in this regard. It is recommended that scheduling problems be assessed and compared to container routing and the type of cargo transported between ports.