An Improved Shuffled Frog-Leaping Algorithm for Solving the Dynamic and Continuous Berth Allocation Problem (DCBAP)

This research deals with the dynamic and continuous berth allocation problem (DCBAP) in which both arrived and incoming ships are considered and a quay is used as a continuous line to accommodate as many ships as possible at one time. The DCBAP is solved by a two-stage procedure. In the first stage a heuristic/metaheuristic is used to generate alternative ship placement sequences while in the second stage a specific heuristic is employed to place ships and resolve overlaps of ships for the development of a feasible solution. Different methods, including FCFS (First Come First Served), SFLA (Shuffled Frog-Leaping Algorithm), and ISFLA (Improved Shuffled Frog-Leaping Algorithm), were employed in the first stage for comparison. The experimental results show that the ISFLA outperforms the others in terms of solution quality, implying that the ISFLA has the potential to deal with the DCBAP in a container terminal.


Introduction
International trade is the exchange of goods or services across nations. There are different kinds of transportations for exchanging goods, including air, land, and maritime. Among these transportations, maritime transport is essential due to its having a relatively lower transport cost as a result of mass transport. Among various kinds of maritime transport, container transport is especially important, as the number of global container shipments is increasing continuously. To serve calling ships, many of the busy global ports, such as Hamburg, Rotterdam, and Antwerp in Europe, and Busan, Shanghai, and Hong Kong, employ multi-user terminals [1]. To improve productivity at this kind of terminal, a better solution for the berth allocation problem (BAP) is essential.
There are many approaches available for improving the productivity of a container terminal. Basically, these approaches can be classified into the two categories: hardware and software approaches. For example, increasing the number of quays or quay cranes (QCs) is an example of a hardware approach, whereas improving the utilization of quays or QCs is an example of a software approach. In this research, the software approach is used, as it is less expensive and time-consuming.
In a container terminal, the operations can mainly be classified into three areas: seaside, yard, and landside. In the seaside area, there are three famous seaside operational problems: the berth allocation in comparison to FCFS and SFLA, which were employed in the first stage. The results showed that the ISFLA outperformed the other two in terms of solution quality.
The rest of this research is organized as follows: Section 2 presents a literature review on the CBAP and basic SFLA; Section 3 presents a formulation of the Mixed Integer Linear Programming (MILP) model for the DCBAP; Section 4 develops the ISFLA for solving the DCBAP; Section 5 includes numerical examples; Finally, Section 6 presents the conclusion and future research directions. Figure 1 illustrates two ships berthing alongside the quay of a container terminal. The container terminal is separated into the following three areas: seaside, yard and landside [3].

The Operations in a Seaport Container Terminal
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 28 search the solution space for a possible solution. We conducted experiments investigating the effectiveness of the ISFLA in comparison to FCFS and SFLA, which were employed in the first stage. The results showed that the ISFLA outperformed the other two in terms of solution quality. The rest of this research is organized as follows: Section 2 presents a literature review on the CBAP and basic SFLA; Section 3 presents a formulation of the Mixed Integer Linear Programming (MILP) model for the DCBAP; Section 4 develops the ISFLA for solving the DCBAP; Section 5 includes numerical examples; Finally, Section 6 presents the conclusion and future research directions. Figure 1 illustrates two ships berthing alongside the quay of a container terminal. The container terminal is separated into the following three areas: seaside, yard and landside [3]. The quay in the container terminal can be configured as either a discrete or continuous quay. With a discrete configuration, the quay is separated into several fixed sections, and each section is allowed to accommodate one ship at a time. With a continuous configuration, the entire quay is used as a continuous line that can accommodate as many ships as possible at one time. Along the quayside, rail-mounted quay cranes are used to load/unload containers to/from trucks. The trucks will transport containers between the berthed ships and the storage areas in the yard area.

The Operations in a Seaport Container Terminal
The yard area serves as a buffer between the sea and land sides, which is an area of temporary storage for containers. The yard area is divided into many blocks, and each block is composed of several rows for stacking containers. For each block, there is a yard crane used to load and unload containers on or off the trucks. The containers in the yard area will be transshipped to other ships or transported to customers via road or rail. Figure 2a shows the flows of intra-terminal transshipment between ship 1 and ship 2. The quay in the container terminal can be configured as either a discrete or continuous quay. With a discrete configuration, the quay is separated into several fixed sections, and each section is allowed to accommodate one ship at a time. With a continuous configuration, the entire quay is used as a continuous line that can accommodate as many ships as possible at one time. Along the quayside, rail-mounted quay cranes are used to load/unload containers to/from trucks. The trucks will transport containers between the berthed ships and the storage areas in the yard area.
The yard area serves as a buffer between the sea and land sides, which is an area of temporary storage for containers. The yard area is divided into many blocks, and each block is composed of several rows for stacking containers. For each block, there is a yard crane used to load and unload containers on or off the trucks. The containers in the yard area will be transshipped to other ships or transported to customers via road or rail. Figure 2a shows the flows of intra-terminal transshipment between ship 1 and ship 2.
The land side contains operations consisting of transporting containers to and from the container terminal, by trucks, road or railway. Figure 2b shows two types of intermodal transportation: ship-to-road and ship-to-rail. In this yard area, both intra-terminal transshipment and intermodal transshipments operate simultaneously.
The yard area serves as a buffer between the sea and land sides, which is an area of temporary storage for containers. The yard area is divided into many blocks, and each block is composed of several rows for stacking containers. For each block, there is a yard crane used to load and unload containers on or off the trucks. The containers in the yard area will be transshipped to other ships or transported to customers via road or rail. Figure 2a shows the flows of intra-terminal transshipment between ship 1 and ship 2. In a container terminal, there are different operational planning problems. These can be listed as follows: • Allocation of berths for calling ships (BAP).

•
Assigning of QCs for berthed ships for loading and unloading containers (QCAP).

•
Assigning of yard trucks.

•
Scheduling of yard trucks.
The BAP, QCAP and QCSP are widely recognized seaside operational problems. In this research, we will focus on the BAP of the continuous quay configuration. Relevant CBAP studies are reviewed in the next section.
[1], Wang and Lim [32], Lee and Chen [33], Zhen et al. [34] carried out studies focused on the CBAP. Due to the focus on quays with continuous configuration, we detail these CBAP studies below.
The study by Lim [4] is an early study focusing on the CBAP. In that study, the CBAP was considered as a two-dimensional packing problem, but transformed into a graph theoretical representation once the authors had captured the CBAP characteristics. Then, the author proposed a heuristic for solving this problem with the objective of minimizing the quay length required to accommodate all incoming ships under the assumption that these ships would not further change their berthing positions. On the other hand, Li et al. [28] solved the CBAP both with and without this restriction in ships' movements. Their objective was to minimize the makespan of the schedule. Guan et al. [29] developed a heuristic for solving the CBAP. Their objective was to find the solution with the minimum total weighted completion time of ships. Park and Kim [30] proposed a subgradient optimization method for dealing with the CBAP with the objective of minimizing costs arising from the delayed departures of ships due to undesirable service order and the additional complexity of handling containers when ships are served at non-optimal mooring locations in the port. This study appears to be more practical than the above-mentioned studies, in that it took the actual berthing positions of ships into consideration. Kim and Moon [5] addressed the same CBAP discussed in Park and Kim [30]. The authors proposed a simulated annealing method to deal with this problem, finally concluding that this approach had the capacity to find an optimal solution. However, the CBAP studies mentioned above neglected the impact of berthing position on the handling time of a ship. In a later study, Park and Kim [31] studied a CBAP with a similar objective function to those used in Park and Kim [30] and Kim and Moon [5]. A major difference of the objective function used in [31] was that it accounted for additional costs resulting from the early or late start of ship handling with respect to the ETA of a ship. Their approach generated the optimal starting times of ship services and berthing positions, as well as QC assignments to ships. However, in that study, the handling times of ships were still independent from the ships' mooring positions. Having completed the DBAP [24][25][26], Imai et al.
[1] moved on to deal with the CBAP, which was treated as a two-dimensional cutting stock problem. The CBAP was solved by using a two-stage heuristic. The first stage generates an initial solution to the DBAP, and then this solution is further transformed into a feasible one by repositioning overlapped ships as well as sparsely located ships at the second stage. In that study, the impacts of the mooring positioning of a ship on its handling time were taken into consideration. Wang and Lim [32] solved the CBAP by proposing a stochastic beam search. Their experiments showed that the proposed approach was able to find a near-optimal solution, with the solution being better than those found by state-of-art metaheuristics and traditional deterministic beam search in terms of accuracy and efficiency. Lee and Chen [33] proposed a neighborhood-search optimization heuristic to solve the CBAP. In that study, factors including the First Come First Served (FCFS) rule, clearance distance between ships, and possible ship shifting were taken into consideration. The experimental results showed that a better decision could be achieved when taking these factors into consideration. In Zhen et al. [34], the authors proposed a two-stage decision model for dealing with the CBAP. The arrival times and handling times of ships in the model are treated as uncertainties. In the first stage, this model first initiates a baseline schedule that is further adjusted by using a reactive recovery strategy in the second stage. The objective of that approach was to minimize penalty costs incurred by deviating from the initial schedule. For the CBAP in a large-scale environment, the authors proposed a metaheuristic approach based on the simulated annealing approach. Our literature review shows there is still a lack of implementation of SFLA-based approaches for dealing with the BAP, especially the DCBAP.
The CBAP is not completely equivalent to the cutting-stock problem (CSP), because not all rectangles (placing items)-which correspond to calling ships in the CBAP-are already available in the CSP, as a result of the incoming ships in the CBAP. In addition, the constraints in the CBAP can be graphically portrayed. For example, some rectangles in a berth plan cannot be placed before corresponding predetermined vertical lines that represent the estimated arrival times (ETAs) of calling ships. This restriction is simply based on the fact that ships cannot be berthed before they have arrived at the port. Thus, the static version of the CBAP is identical to the fixed-orientation CSP but the dynamic BAP is not [1].
In many CBAP studies, such as those by Lim [4], Li et al. [28], and Guan et al. [29], the authors only considered arrived ships and assumed that a ship's handling time was independent of its berthing location. Although incoming ships were taken into consideration in some studies, such as those by Park and Kim [30,31] and Kim and Moon [5], and fixed handling times had been assumed for those ships, Ref. [31] can only be regarded as a static version of CBAP, as ships were able to be served before their ETAs, while penalty costs were then imposed in the objective function for the use of such services. This relaxation downgrades the dynamic version of BAP to a static one [5]. Consequently, its constraint is also equivalent to that of the CSP with fixed orientation. In this present research, the CBAP to be solved is quite different from the conventional CSP, as the handling time of a ship is dependent on the ship's berthing location, and incoming ships are to be taken into consideration.
In addition, heuristics have been widely used to deal with the CBAP. However, high-level heuristics, i.e., metaheuristics, have seldom been used to deal with the CBAP. In particular, there is still a lack of using SFLA to solve the BAP. Such applications are still at the emerging stage.

The Basic SFLA
As a kind of population and evolutionary-based metaheuristic, the basic SFLA (B-SFLA) employs a swarm of frogs to find the position (solution) with the greatest amount of food. The goodness of a position can be indicated by an objective function value. In a D-dimensional space, the current position of a frog f is denoted as X f (t) = (X f 1 , . . . , X f D ). During memetic evolution, frogs are separated into m groups, with the best frog being assigned to the memeplex 1; the 2nd best to the memeplex 2; the mth best to the memeplex m; the m + 1th best then back to the memeplex 1 again, and so on. As a result, each group has n = F/m frogs, and these frogs are ordered decreasingly according to their objective function values [35]. Following this, triangular distribution sampling, as in Equation (1), is used to form sub-memeplexes.
With Equation (1), the best frog in each memeplex has the highest probability 2/(n + 1), while the worst frog has the lowest probability 2/n(n + 1), of being selected into a corresponding sub-memeplex. Then, Roulette Wheel Selection (RWS) can be used to select q frogs from the original memeplex into the sub-memeplex for evolution and local searching. At the current time t, the leaping distance of the frog f is denoted as D f (t) and determined by Equation (2).
The R is a random number in [0, 1]; X b (t) is the position of the best frog in the sub-memeplex; S m is the maximum step allowed. The next position of the frog f is determined by Equation (3).
If feasible and better, this next position will be accepted; otherwise, the X b (t) in Equation (2) will be replaced with X g (t) and the next position will be found again, where X g (t) is the position of the global best frog. If this next position is once again not better, a random jump will be made for the frog. Once all frogs in a sub-memeplex have completed several local searches, they then move on to process the next sub-memeplex until all sub-memeplexes have been completed. After completing one iterative run, the frogs are reshuffled and the next iterative run is started. This is carried out until the termination criteria are reached. For more details about the B-SFLA, one can refer to Eusuff et al. [35].
Our literature review shows that SFLA and its variants have been used in various areas, including water distribution network design [19], project management [36], knapsack problem [37], component pick-and-place problem [38], routing problem [39], gateway loading balance [40], local dimming optimization [41], underwater sonar image detention [42], job scheduling [43], big data [44], feature selection [45], etc. However, to our best knowledge, SFLA has never been used to deal with the BAP. the length of a quay; H: the planning horizon; V: a finite set of calling vessels; V = {1, 2, . . . , N}, where N is the total number of calling ships; Ç: a finite set of berthing constraints; P: a set of plans with assignments of berthing positions and start berthing times for ships; ∀ p i ∈ P, ∃ p i = ∪ j=N j=1 B j , S j 0 ≤ B j ≤ L − l j , 0 ≤ S j , where B j , S j is a pair of assignments, where B j is the berthing position, S j is the start berthing time and l j is length of the calling ship j. Z: an objective function mapping p i to a time/cot value Z. The objective of the DCBAP is to find a p i or p * (p i , p * ∈ P) that minimizes the objective function value Z, in which p i is a feasible solution, while p * is the optimal solution subjecting to Ç. For each ship j a berthing position B j , S j on P is assigned. Finding the B j , S j for each ship j (j = 1, . . . , n, where n is total number of calling ships) is an NP-hard problem.

Berthing Plan Representation for the DCBAP
A berth plan can represent a solution to the CBAP. Figure 3 illustrates a berth plan in which the X-axis represents the time dimension (planning horizon), while the Y-axis represents the space dimension (quay length). This plan includes one overlap of two ships j and k that are respectively located at positions at a (x j 0 , y j 0 ) and a x k 0 , y k 0 , where x j 0 and x k 0 are berthing times and y j 0 and y k 0 are berthing positions. However, this overlap makes this plan unfeasible. To be feasible, it needs to resolve this overlap. A ship being selected and moved to resolve an overlap is termed "target ship". It is noted that the arrival time (ETA j ) and the desired berthing position (d j ) are the best coordinates in a berthing plan for a ship j, as there are no increased waiting or handling times. Any deviation from these coordinates will increase the cost for the ship.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 28 are berthing positions. However, this overlap makes this plan unfeasible. To be feasible, it needs to resolve this overlap. A ship being selected and moved to resolve an overlap is termed "target ship". It is noted that the arrival time (ETAj) and the desired berthing position (dj) are the best coordinates in a berthing plan for a ship j, as there are no increased waiting or handling times. Any deviation from these coordinates will increase the cost for the ship.

Estimation of Increased Handling Times for a Ship
In this research, it is assumed that the berthing position of a ship will affect the ship's handling time. In Figure 3, we give each ship j or k an initial handling time based on their respective desired berthing positions. Thus, if any of the ships deviates from its desired berthing position, then the ship's handling time will increase due to the greater distance that a truck will need to move a container.
Equation (4) is the formula for calculating the distance (∆ ) for a ship j deviating from its desired berthing position.
where is the actual berthing position of ship j alongside the quay; is the desired berthing position of ship j alongside the quay.
Then, the increased handling times (∆ ) for the ship j is determined by Equation (5) where is the berth deviation factor; a rate estimating the increase of handling time due to the deviation from the desired berthing position ( ).
Equation (6) is used to estimate , the actual handling times of a ship j.
where Figure 3. The coordinates of two corners of the ships j and k.

Estimation of Increased Handling Times for a Ship
In this research, it is assumed that the berthing position of a ship will affect the ship's handling time. In Figure 3, we give each ship j or k an initial handling time based on their respective desired berthing positions. Thus, if any of the ships deviates from its desired berthing position, then the ship's handling time will increase due to the greater distance that a truck will need to move a container.
Equation (4) is the formula for calculating the distance (∆b j ) for a ship j deviating from its desired berthing position.
where B j is the actual berthing position of ship j alongside the quay; d j is the desired berthing position of ship j alongside the quay.
Then, the increased handling times (∆H j ) for the ship j is determined by Equation (5) where β is the berth deviation factor; a rate estimating the increase of handling time due to the deviation from the desired berthing position (d j ).
Appl. Sci. 2019, 9, 4682 8 of 28 Equation (6) is used to estimate T j , the actual handling times of a ship j.
where h j is the original handling time for ship j (this is not affected by berthing position). Finally, given the actual start berthing time of ship j (S j ), the completion time of a ship j can be estimated by Equation (7).

Estimation of Increased Waiting for a Ship
In Figure 1, one way of resolving the overlap is to move ship j in the +X direction. This will not increase the handling times, but rather the waiting times of the ship. The increased waiting times (∆W j ) for a ship j can be estimated by Equation (8), where S j is the start berthing time of the ship j.

Mathematical Model
This section formulates a mathematical model for the DCBAP. Firstly, the assumptions, parameters, indices, and decision variables for this model are introduced.

Assumptions
(1) Each ship is handled continuously until it is completed.
Each ship has an estimated processing time.
Each ship has a desired berthing position. (4) Inter-ship clearance distance is included in ship length. X jxy = 1, if the ship j is assigned with the berthing postion (x, y) 0, otherwise B j the berthing position of ship j ( j ∈ J) S j the actual start berthing time of ship j ( j ∈ J) C j the completion time of ship j (j ∈ J) ∆W j the increased waiting time of ship j ( j ∈ J) ∆H j the increased handling time of ship j ( j ∈ J).
The Mixed Integer Programing (MIP) model of the DCBAP can be formulated as follows.
Equation (9) is the objective function for minimizing the total cost (Z), including the increased waiting cost, as well as the increased handling cost. Equation (10) stipulates that one berthing position can be assigned to only one ship at a time. Equation (11) is a requirement for ∆W j stipulating that S j should not start before ETA j . Equation (12) defines a requirement for ∆H j . Equations (13) and (14) are conditions for the assigned berthing position for a ship. Equations (15) and (16) are physical conditions for a ship. Equations (17) comprise a binary constraint for decision variables X jxy .

The ISFLA for Solving the Simultaneous DCBAP
Due to the NP-hard problem, the MILP formulated in Section 3.3 will become computationally intractable when used to deal with problems of practical size. This section details the development of an ISFLA for dealing with the DCBAP.

Position Representation
Equation (9) is the objective function for minimizing the total cost (Z), including the increased waiting cost, as well as the increased handling cost. Equation (10) stipulates that one berthing position can be assigned to only one ship at a time. Equation (11) is a requirement for ∆ stipulating that should not start before . Equation (12) defines a requirement for ∆ . Equations (13) and (14) are conditions for the assigned berthing position for a ship. Equations (15) and (16) are physical conditions for a ship. Equations (17) comprise a binary constraint for decision variables .

The ISFLA for Solving the Simultaneous DCBAP
Due to the NP-hard problem, the MILP formulated in Section 3.3 will become computationally intractable when used to deal with problems of practical size. This section details the development of an ISFLA for dealing with the DCBAP.  The ROV technique used in Hsu (2016) [7] is used to transform real values into ranking numbers that represent a placement sequence of ships. For instance, the sequence [0.5,0.6,0.7,0.4] will be transformed into the ranking set [2,3,4,1], corresponding to the placement sequence for ships 1 to 4.

The ISFLA
The ISFLA possesses the following features: (1) multiple groups of frogs, (2) shuffled mechanism, (3) discrete operators with self-adaptive jump, and (4) self-adaptive mutation The ROV technique used in Hsu (2016) [7] is used to transform real values into ranking numbers that represent a placement sequence of ships. For instance, the sequence [0.5,0.6,0.7,0.4] will be transformed into the ranking set [2,3,4,1], corresponding to the placement sequence for ships 1 to 4.

Multiple Groups of Frogs
For a given solution space, the ISFLA uses multiple groups of frogs and forms multiple search areas (focuses) that can be concurrently searched by the swarm of frogs. Frogs in the same group will search among them for the best frog in that group. Figure 5 shows the evolution of three groups of frogs from the iteration i to the iteration i + 1. The ISFLA does not use sub-memeplex or triangular distribution sampling. This allows all frogs in a memeplex to attend to evolution directly, which results in population advantage.

Shuffled Mechanism
The shuffled mechanism is used to regroup frogs belonging to the same swarm, which can initiate other opportunities for these frogs to be affected by other elites (i.e., the best frog in each group). This helps diversify frogs for the exploration of the solution space. In a D-dimensional solution space, the position of a frog i is denoted as = ( ,… , ), and the frogs of this swarm are separated into m groups according to the following rules:  the best is assigned to group 1;  the 2nd best is assigned to group 2;  the mth best is assigned to group m;  the m + 1th best is assigned to group 1 again, and so on.
As a result, each group then contains F/m frogs (assuming that there are a total number of F frogs), and the frogs in each group are ordered in decreasing order based on their objective function values.

Discrete Operators with Self-Adaptive Jump
In this research, the ISFLA uses discrete operators to determine the next position for a frog based on the position of the frog relative to that of the target frog. An example is given in order to specify the self-adaptive jump of frogs. Figure 6 illustrates three frogs, , i, and j, with positions of ( ) = [4,3,2,1,5,6], ( ) = [1,2,3,4,5,6], and ( ) = [4,1,3,2,5,6], respectively. Among the three frogs, the frog o, termed the "target frog", is in the best position, attracting frogs i and j to search around it. The frogs and are considered to be closer, due to sharing four of the same position elements (while frogs and only shar two of the same elements). A self-adaptive mechanism is required to cause the non-target frogs to adaptively approach the target frog, i.e., to enable a bigger jump for more distant frogs to be able to approach the target frog quickly, while enabling smaller jumps for closer frogs to prevent them from jumping over the optima. However, non-target frogs should not jump directly to the target frog, as this will waste a local search.

Shuffled Mechanism
The shuffled mechanism is used to regroup frogs belonging to the same swarm, which can initiate other opportunities for these frogs to be affected by other elites (i.e., the best frog in each group). This helps diversify frogs for the exploration of the solution space. In a D-dimensional solution space, the position of a frog i is denoted as X i = (X i 1, . . . ,X iD ), and the frogs of this swarm are separated into m groups according to the following rules: • the best is assigned to group 1; • the 2nd best is assigned to group 2; • the mth best is assigned to group m; • the m + 1th best is assigned to group 1 again, and so on.
As a result, each group then contains F/m frogs (assuming that there are a total number of F frogs), and the frogs in each group are ordered in decreasing order based on their objective function values.

Discrete Operators with Self-Adaptive Jump
In this research, the ISFLA uses discrete operators to determine the next position for a frog based on the position of the frog relative to that of the target frog. An example is given in order to specify the self-adaptive jump of frogs. Figure 6 illustrates three frogs, o, i, and j, with positions of X o (t) = [4,3,2,1,5,6], X i (t) = [1,2,3,4,5,6], and X j (t) = [4,1,3,2,5,6], respectively. Among the three frogs, the frog o, termed the "target frog", is in the best position, attracting frogs i and j to search around it. The frogs i and o are considered to be closer, due to sharing four of the same position elements (while frogs j and o only shar two of the same elements). A self-adaptive mechanism is required to cause the non-target frogs to adaptively approach the target frog, i.e., to enable a bigger jump for more distant frogs to be able to approach the target frog quickly, while enabling smaller jumps for closer frogs to prevent them from jumping over the optima. However, non-target frogs should not jump directly to the target frog, as this will waste a local search. Figure 6 illustrates three frogs, , i, and j, with positions of ( ) = [4,3,2,1,5,6], ( ) = [1,2,3,4,5,6], and ( ) = [4,1,3,2,5,6], respectively. Among the three frogs, the frog o, termed the "target frog", is in the best position, attracting frogs i and j to search around it. The frogs and are considered to be closer, due to sharing four of the same position elements (while frogs and only shar two of the same elements). A self-adaptive mechanism is required to cause the non-target frogs to adaptively approach the target frog, i.e., to enable a bigger jump for more distant frogs to be able to approach the target frog quickly, while enabling smaller jumps for closer frogs to prevent them from jumping over the optima. However, non-target frogs should not jump directly to the target frog, as this will waste a local search.  To enable self-adaptive jumps, we change Equations (2)- (18).
RV i (t) refers to the adaptive binary vector, and is used to adjust the moving distance D i (t) for a frog i at the time t.
The operator "~" refers to the total distance measuring operator, which measures the distance between frogs o and i in cases where frog i jumps towards the target frog o. The operator works as follows.
In Equation (19), k indicates the kth elements in the position vectors of frog o and i. Given In Equation (18), the operator " " refers to the binary step multiply operator. It determines the result of the multiplication of X o (t)~X i (t) with RV i (t) = [RV i,1 , . . . , RV i,D ], in which each RV i,k (t) (k ∈ [1,D]) is a binary value. This operator works as shown in Equation (20).
For example, given RV i (t) = [0,1,0,0,0,0], we can derive D i (t) as follows. The above example shows that more binary values of 1 in RV i (t) result in a larger leaping distance. Thus, we let each RV i,k (t) be determined by Equation (21).
R1 is a random number within the range [0,1], while BR1 i (t) is a threshold controlling the generation of the binary values 0 or 1 in RV i (t) for the frog i. To enable a bigger jump for a more distant frog, Equations (22)- (24) are designed to initiate a bigger BR1 i (t). First, Equation (22) counts the total number of shared elements between X i (t) and X o (t), denoted as NSE i,o (t). Then, Equation (23) calculates ω i (t), the maximum number of binary values of 1 currently allowed for the RV i (t) of frog i. Finally, Equation (24) is used to determine the value of BR1 i (t).
In Equation (22), the symbol "≈" is the comparing operator, and compares two elements with the same position in X o (t) and X i (t), and this operator works as follows.
In Equation (26), the symbol " " is as the leaping operator, which is used to determine the next position of a frog. It works on the basis of the following 4 steps: (1) the operator takes the first non-zero value out of the D i (t) and replaces the value at the same position in X i (t), (2) the replaced value takes the position of the non-zero value in X i (t), (3) repeat the first and second steps until there are no non-zero values in D i (t), (4) the operator copies the current X i (t) as X i (t + 1).
Take Equation (27) as an example: the operator first takes the non-zero value "3" from D i (t) and replaces the "2" in X i (t). Then, the replaced value "2" takes the place of the "3" in X i (t). Because there are no non-zero elements in D i (t), the next position of the frog i will be X i (t + 1) = [1, 3,2,4,5,6].
This example shows that a non-zero element in D i (t) can exchange two elements in X i (t). In terms of distance, the two frogs i and o are now considered to be closer due to their having four shared position elements, and the frog i possesses the following data: In addition to the self-adaptive leap mechanism, the ISFLA also makes it possible to include the following novel features: • Direct-jump prevention: The direct-jump prevention will stop the generation of binary values of 1 for the RV i (t) of frog i in order to prevent it jumping directly to the target frog, as this would waste one local search. This will happen if the frog i is satisfied with the condition NSE i,o (t) ≥ D − 2.
In Equation (22), the term −2 refers to the direct-jump prevention, which prevents frog i from jumping directly to the target frog o.

•
Neighborhood jump: However, if BR1 i (t) 0, there is always a chance for a frog to jump to the target frog directly. Thus, we let the ISFLA count the total number of binary values of 1 generated so far for the RV i (t) of frog i; if the threshold ω i is reached, instead of allowing a direct jump, the swap(p1, p2) operator, which switches two randomly selected elements in a frog's position vector, will be used to perform a neighborhood search. The neighborhood search avoids wasting local searches for frogs due to the imposition of direct-jump prevention.

•
Push jump: For a frog i freed from direct-jump prevention but idling at its current position, which happens when D k=1 RV i,k (t) = 0, the ISFLA will introduce one random binary value 1 into RV i (t) in order to push jump this frog.

The Self-Adaptive Mutation Mechanism
The ISFLA uses a self-adaptive mutation mechanism to vary the frogs in order to prevent them from being trapped in local optima. Two mutation operators used in the ISFLA are detailed below. TM has a higher degree of mutation than SM due to the greater number of mutated genes.
The trigger for the mutation of a frog depends on a random number R2 and a pre-defined mutation rate θ; if R2 < θ, then the frog is mutated. Furthermore, the use of SM or TM depends on Equation (28), which measures the degree of similarity between X i (t) and X g (t).
The higher the SD i,g (t), the higher the similarity. We assume that frogs with less similarity require greater variations so as to be able to escape the unfavorable position. Thus, if SD i,g (t) < δ, then TM is used in order to produce a bigger variation; otherwise, SM is used in order to generate a smaller variation. The δ, with its value within the range [0, 1], is a parameter controlling the use of TM or SM.

The Two Stage Procedure
In addition to the ROV technique, a two-stage approach was used to deal with the DCBAP.

The First Stage
This first stage focuses on generating a ship placement sequence.

The Second Stage
The second stage focuses on developing a feasible DCBAP solution after resolving overlaps of ships. This process contains three steps: place ships one by one into a berth plan, identify each overlap of ships, and resolve the overlap. These three tasks are repeated until all ships have been placed into the berth plan without any overlap.
Step 1: Place ships one by one into the berth plan From the objective function, we know that if every ship is assigned its best position (ETA j , d j ), this will lead to the optimal berth plan (i.e., Z = 0) due to there being no increased waiting or handling costs for the ships. However, if there are too many ships calling a container terminal at the same time, then this ideal arrangement becomes impossible due to the limited available resources in terms of both space and time. Repositioning overlapped ships is necessary, and repositioning around the best position will lead to the generation of optimal/near-optimal solutions.
Step 2: Identify overlaps Figure 7 shows a berth plan with two ships that are overlapping in the dimensions of "space" and "time". Given that a = (x j 0 , y j 0 ) and c = (x j 1 , y j 1 ) are the coordinates of the lower-lef and upper-right corners of the ship j, while a = (x k 0 , y k 0 ) and c = (x k 1 , y k 1 ) are the coordinates of the lower-left and upper-right corners of the ship k, respectively, then Equation (29) is used to check whether the two ships are overlapping. If the following coordinate relationships hold simultaneously, then the two ships are concluded to be overlapping.
From the objective function, we know that if every ship is assigned its best position (ETAj, dj), this will lead to the optimal berth plan (i.e., Z = 0) due to there being no increased waiting or handling costs for the ships. However, if there are too many ships calling a container terminal at the same time, then this ideal arrangement becomes impossible due to the limited available resources in terms of both space and time. Repositioning overlapped ships is necessary, and repositioning around the best position will lead to the generation of optimal/near-optimal solutions.
Step 2: Identify overlaps Figure 7 shows a berth plan with two ships that are overlapping in the dimensions of "space" and "time". Given that a = ( , ) and c = ( , ) are the coordinates of the lower-lef and upper-right corners of the ship j, while a' = ( , ) and c' = ( , ) are the coordinates of the lower-left and upper-right corners of the ship k, respectively, then Equation (29) is used to check whether the two ships are overlapping. If the following coordinate relationships hold simultaneously, then the two ships are concluded to be overlapping. Step 3: Resolve the overlap Resolving all overlaps of ships is necessary in order to make a feasible plan. In this research, a target ship is permitted to move in one of the following three directions: +Y, −Y, and +X, in order to resolve an overlap. However, while moving towards the +Y/−Y direction, handling times will be increased for the target ship due to the longer moving distances that trucks will have to move containers ( Figure 8). When moving towards the +X direction, waiting time (cost) will be increased for the target ship ( Figure 9). The lower the moving cost, the higher the priority.
For a target ship j moving toward the +Y/−Y direction, Equation (30) is used for estimating the increased distance ∆Y j .
Equation (31) is the cost of increased handling times, where c 2 is the cost rate of handling time, β = 2/(100 * 60), and ∆b j = ∆Y j .
After repositioning, the coordinates of the lower-left and upper-right corners of the target ship j become a = (x j 0 , y j 0 + ∆Y j ) and c = (x j 1 , y j For a target ship j moving towards the +X direction, Equation (32) is available for estimating the increased waiting time.
Equation (33) is the increased waiting cost, where c 1 is the cost rate of waiting time.
After repositioning, the coordinates of the lower-left and upper-right corners of the target ship j become a = (x  Step 3: Resolve the overlap Resolving all overlaps of ships is necessary in order to make a feasible plan. In this research, a target ship is permitted to move in one of the following three directions: +Y, −Y, and +X, in order to resolve an overlap. However, while moving towards the +Y/−Y direction, handling times will be increased for the target ship due to the longer moving distances that trucks will have to move containers ( Figure 8). When moving towards the +X direction, waiting time (cost) will be increased for the target ship ( Figure 9). The lower the moving cost, the higher the priority.
For a target ship j moving toward the +Y/−Y direction, Equation (30) is used for estimating the increased distance ∆ .  Step 3: Resolve the overlap Resolving all overlaps of ships is necessary in order to make a feasible plan. In this research, a target ship is permitted to move in one of the following three directions: +Y, −Y, and +X, in order to resolve an overlap. However, while moving towards the +Y/−Y direction, handling times will be increased for the target ship due to the longer moving distances that trucks will have to move containers ( Figure 8). When moving towards the +X direction, waiting time (cost) will be increased for the target ship ( Figure 9). The lower the moving cost, the higher the priority.
For a target ship j moving toward the +Y/−Y direction, Equation (30) is used for estimating the increased distance ∆ .  Equation (31) is the cost of increased handling times, where is the cost rate of handling time, = 2/(100 * 60), and ∆ = |∆ |. ∆ Figure 9. The movement of a target ship toward the +X direction. Figure 10 shows the main flow of the two-stage approach for dealing with the DCBAP. The first stage (steps 1 and 2) focuses on generating a ship placement sequence by using different approaches. The second stage (steps 3 to 10) is a heuristic used to place ships into a berth plan and resolve overlaps of ships. Each step is detailed as follows.

The Main Flow Chart of the Two-Stage Approach
Step 1: Generate data (including a j , d j , l j and h j ) for all calling ships (j = 1, . . . , n; where n is the total number of calling ships).
Step 2: Develop a placement sequence using the ISFLA (or FCFS, GA that are to be compared); set s = 1, where s indicates the placement order of a ship.
Step 3: Place the ship j (at the placement order s) into the berth plan, with the lower-left and upper-right corners being located at the coordinates (x j 0 , y j 0 ) = (a j , c j ) and (x j 1 , y j 1 ) = (a j + h j , c j + l j ), respectively; Set k = 1.
Step 4: Check whether the ship j with the placement order s has overlapped with the ship with the placement order s-k by using Equation (29). If "No", then go to Step 5; otherwise, go to Step 8.
Step 5: Check whether s = k. If "Yes", go to Step 6; otherwise k = k + 1 and go to Step 4 to check the next ship.
Step 6: Check whether the ship j remains with overlap(s) after repositioning. If "Yes", go to Step 10; otherwise, go to Step 7.
Step 7: Check whether s = N. If "Yes", go to Step 11; otherwise s = s + 1 and go to Step 3.
Step 8: Estimate the moving cost of each moving direction for repositioning the target ship with the placement order s using Equations (31) and (33). In addition, let i = 1 index the least-cost movement direction as the first priority; i = 2 index the second least-cost movement direction as the second priority; i = 3 index the most-cost movement direction as the 3rd priority.
Step 9: Repositioning the target ship towards the direction with the priority index i. Update the coordinates of the lower-left and upper-right corners of the target ship using Equation (30) or (32). Go to Step 5.
Step 10: Check whether the ship j has been ever overlapped with this ship or had an overlap resolved.
If "Yes", set i = i + 1 (choose the next repositioning direction) and then go to Step 9; otherwise, go to Step 8.
Step 11: Check whether there is another iteration to perform. If "Yes", go to Step 2; otherwise go to Step 12.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 16 of 28 Figure 10. The main flow of the two-stage procedure. Figure 10. The main flow of the two-stage procedure.

Numerical Examples
Java is used as the programming language to implement the different approaches for the purpose of comparison. In the two-stage procedure, the FCFS, SFLA, and ISFLA were each employed in the first stage, whereas a specific heuristic was employed in the second stage. For each comparison, a set of ship data was randomly generated by a computer, and experiments were conducted on a computer equipped with an Intel PENTIUM CPU (64 bit and 1.8 GHz) and 4G DRAM. Section 5.1 details the parameter values set for the different approaches. Section 5.2 shows the results obtained from a small-sized experiment with N = 10. Section 5.3 shows the results obtained from a large-sized experiment with N = 50. Section 5.4 presents an analysis and discussion on the obtained experimental results.

Parameter Settings for Different Approaches
The following parameter settings were used for the experiments: the planning horizon (H) was set to one week, so that the ETAs of ships would be within the interval [0, 168] (hours); the quay length (L) was set to 1000 m; the length of a ship j (l j ) was a random value within the range [50, 200] (meters); the desired berthing position of a ship j (d j ) was a random value within the interval [0, L-l j ]; the cost rates c 1 and c 2 were set to USD 1000/h; and the β = 2/(100 * 60) (per hour), which means the deviation of 100 m from the desired berthing position of a ship would increase handling time by 2 min.
The parameter values for the three different methods used in the first stage of the two-stage procedure are listed in Table 1. The parameter F indicates the total number of frogs in the swarm; the parameter m indicates the total number of memeplexes; the parameter l_iter specifies the number of local search for each frog; the parameter iterations indicates the total number of iterations; the parameter q indicates the total number of frogs in each memeplex; the parameter δ indicates the rate controlling the use of TM; the parameter θ indicates the mutation rate used; the symbol "-:" indicates not used.
F: total number of frogs; m: total number of memeplexes; l_iter: number of local search; iterations: total number of iterations; q: number of frogs in each memeplex; δ: the rate of using TM; θ the mutation rate; -: not used; N: total number of ships. Table 2 shows the original and repositioned ship data of a small-sized experiment (N = 10). The field SID indicates the ID of a ship; the l j shows the length of ship j; the ETA j is Estimated Time of Arrival of ship j; d j is the desired berthing position of ship j; ∆W j shows the increased waiting time of ship j; ∆H j shows increased handling time of ship j; (X o , Y o ) and (X 1 , Y 1 ) are coordinates of left-lower corner and right-upper corners of a ship j, respectively. Two ships (7 and 2) were repositioned to resolve overlaps. Ship 7 is moving towards the +Y direction to avoid overlapping with ship 9, while ship 2 is repositioning towards the −Y direction to avoid overlapping with ship 1. Figure 11 shows the berth plan with overlapping ships, while Figure 12 shows the feasible berth plan after resolving all overlaps of ships. Table 2. Original and repositioning ship data after resolving overlaps (N = 10). while ship 2 is repositioning towards the −Y direction to avoid overlapping with ship 1. Figure 11 shows the berth plan with overlapping ships, while Figure 12 shows the feasible berth plan after resolving all overlaps of ships.     Table 3 shows the original and repositioned data of a large-sized experiment with N = 50 ships. The 1st column shows the ship ID (SID); the 2nd column shows the ship length; the 3rd column shows the Estimated Time of Arrival of ship j ( ); the 4th column shows the desired berthing position of  Table 3 shows the original and repositioned data of a large-sized experiment with N = 50 ships. The 1st column shows the ship ID (SID); the 2nd column shows the ship length; the 3rd column shows the Estimated Time of Arrival of ship j (ETA j ); the 4th column shows the desired berthing position of ship j (d j ); the 5th column shows the original coordinates of ship j (X 0 , Y 0 , X 1 , Y 1 ); the 6th column shows the repositioned coordinates of ship j (X 0 , Y 0 , X 1 , Y 1 ); the 7th column shows the increased waiting time of ship j (c 1 ∆W j ); the 8th column shows the increased handling time of ship j (c 2 ∆H j ). Figure 13 shows a berth plan that is unfeasible due to the overlapping of the ships, based on the original coordinates of ships. To derive a feasible solution, it is necessary to resolve all overlaps of ships. Figure 14 shows the feasible berth plan obtained from the ISFLA after resolving all overlaps among the ships. In Table 3, the bolded figures indicate updated coordinates for those target ships that were repositioned to avoid overlapping.   Table 4 shows the results obtained from different approaches with different problem sizes. The Z indicates the objective function value; the T indicates the computational times; the GAP indicates the gap in percentage compared to the ISFLA. The experimental results show that, on average, the ISFLA exhibits better performance in terms of solution quality. Figure 15 gives a graph with the average Zs obtained from different approaches under different problem sizes (N = 10, 20, 30, 40, 50 and 60).   Table 4 shows the results obtained from different approaches with different problem sizes. The Z indicates the objective function value; the T indicates the computational times; the GAP indicates the gap in percentage compared to the ISFLA. The experimental results show that, on average, the ISFLA exhibits better performance in terms of solution quality. Figure 15 gives a graph with the average Zs obtained from different approaches under different problem sizes (N = 10, 20, 30, 40, 50 and 60).    Table 4 shows the results obtained from different approaches with different problem sizes. The Z indicates the objective function value; the T indicates the computational times; the GAP indicates the gap in percentage compared to the ISFLA. The experimental results show that, on average, the ISFLA exhibits better performance in terms of solution quality. Figure 15 gives a graph with the average Zs obtained from different approaches under different problem sizes (N = 10, 20, 30, 40, 50 and 60).

Analysis of Results
(1) Different methods used in the first stage of the two-stage procedure can lead to different planning results. (2) Among the FCFS, SFLA, and ISFLA, the experimental results showed the ISFLA outperformed the two others, but at the cost of longer computational times. (3) The FCFS rule is simple and fast in finding a solution. However, with an expected lower waiting time (cost), the FCFS rule lacks the capacity to improve solution quality continuously. (4) The increase in iterative runs enables ISFLA to explore a wider solution space. This helps improve solution quality, but at the cost of longer computational times. (5) When the ISFLA repositions a target ship towards the +Y/−Y direction, the handling time for the target ship will increase; when repositioning a target ships towards the +X direction, the waiting time for the target ship will increase. Since the movement towards the +Y/−Y direction results in a smaller increase in handling cost for a target ship, these two directions will have a higher priority, and this can result in better utilization of quay space. However, due to the limited quay space, movement towards the +X direction is sometimes necessary for a target ship. Such repositioning of ships leads to the finding of a near or the optimal solution. (6) As the average computational times required for the ISFLA to solve a problem size of 60 calling ships is about 1.4 h, the ISFLA is thus concluded to be applicable in practice. (7) For the ISFLA, increasing the total number of iterative runs usually leads to the finding of a higher-quality solution, due to the wider exploration of the solution space. (8) Both cost factors c1 and c2 are found able to affect the selection of repositioning direction for a target ship. They can affect the choices of moving direction toward +Y, −Y, or +X for a target ship. In this research, the two cost rates are set to USD 1000/h.

Discussion
(1) To our best knowledge, this is the first research employing SFLA for dealing with the BAP in a seaport container terminal. The SFLA also shows potential in dealing with problems in the yard and landside areas, in addition to the seaside area. (2) Some early CBAP studies, such as those by Lim [4], Li et al. [28], and Guan et al. [29], assumed static berth allocation and independent handling times of ships in their berthing positions. In such studies, the problems to be solved are equivalent to the CSP with fixed orientation of placed items. In contrast, this present research considers the dynamic nature of ship arrivals and assumes that the handling times of ships are dependent on their berthing positions. (3) This present research also differs from the studies of Park and Kim [30,31] and Kim and Moon [5], as those studies assumed fixed handling times for ships. In addition, in a strict sense, the BAP solved in the study of Park and Kim [31] was downgraded to a static version of the problem due to ships being able to be served before their ETAs, albeit with some penalty cost imposed in the objective function [1]. (4) The CBAP solved in this study is quite different from the conventional CSP due to the following three differences: first, in the CBAP a variable handling time is considered; second, the CBAP considers the dynamic nature of ship arrivals; third, the CBAP considers the orientation of the ship placed in a berth plan. (5) Ref.
[1] also proposed a two-stage approach for dealing with the CBAP. The first stage generates an initial solution for the DBAP, and then this solution is further transformed into a feasible one by repositioning overlapped and sparsely located ships at the second stage. However, this approach appears to be rigorous, because to find the final solution to the CBAP, it is first necessary to find the solution to the DBAP. In addition, in the first stage, it is necessary to set a minimum and maximum berth length for the DBAP, and to prepare several intermediate berth lengths for this approach. Such calculations for the setting of parameters take a long time. In contrast, our approach appears to be simpler, as it finds the solution to the CBAP directly. (6) In [1], the proposed heuristic arranges ships in ascending order of their handling time, which is similar to the FCFS rule, before resolving overlaps of ships. This differs from our approach, which creates alternative placement sequences of ships by using the ISFLA so as to explore more alternative solutions. (7) Based on a simulated annealing approach, the metaheuristic proposed in [34] is capable of exploring alternative solutions by means of an R parameter that denotes the number of sequences of the neighborhood search. However, the sequence approach proposed in [34] is different from the ISFLA proposed in this research. In addition, [34] did not describe the resolution of overlapping ships. (8) Compared with [1,34], our approach appears to be simpler and more practically applicable.
However, it is too early to say that our approach is better than the approaches proposed in [1,34], as more experiments are required for comparison. Nevertheless, according to the objective defined and used in this research, we know that the optimal berthing position and time for a ship in a berth plan are the d j and ETA j of the ship. As the ISFLA will assign a ship to or around the optimal berthing position (d j ) and start working time (ETA j ) for each ship, we can conclude that the resulting berth plan will be the optimal solution (when Z ≥ 0), or a near-optimal solution (when Z > 0). As the average computational times for our approach when solving a problem with 60 ships was about 1.4 h; this shows the applicability of the proposed method in practice. (9) This ISFLA is able to improve the simplicity of simple heuristics while avoiding the computational intractability of exact approaches by tuning the total number of iterative runs.

Conclusions and Future Research Direction
The dynamic and continuous berth allocation problem is a seaside operational problem routinely faced by container terminal planners. As this problem can considerably affect the productivity of a container terminal, a better solution for this problem is needed. In this research, we have proposed a two-stage approach for dealing with this problem. In the first stage, the First Come First Served rule (FCFS), Shuffled Frog-Leaping Algorithm (SFLA), and Improved Shuffled Frog-Leaping Algorithm (ISFLA) were each employed in order to generate alternative ship placement sequences that could serve as inputs for the second stage. In the second stage, a specific heuristic was used to assign berthing positions and resolve overlaps of ships in order to develop a feasible solution. The experimental results showed that the Improved Shuffled Frog-Leaping Algorithm outperformed the other two methods in terms of solution quality. The contributions of this research are highlighted below: (1) A novel Shuffled Frog-Leaping Algorithm-based approach (i.e., the ISFLA) was developed to deal with the dynamic and continuous berth allocation problem. To our best knowledge, this kind of approach has never been used for this purpose in a container terminal. (2) In addition to the Improved Shuffled Frog-Leaping Algorithm, a heuristic was developed, in the second stage of a two-stage procedure, for placing ships and resolving overlaps of ships for the development of a feasible solution. (3) The Java programing language was used to implement the two-stage procedure, as it facilitates the generation of BAP solutions for practical use. Our small-sized experiments showed that the proposed approach was capable of finding optimal/near-optimal solutions in terms of the objective function defined in this research. In addition, our experiments demonstrated the feasibility of the ISFLA with respect to computational time for solving a large-sized problem of up to 60 ships.
In this research, we achieved the first step in dealing with the dynamic and continuous berth allocation problem (DCBAP) by using an Improved Shuffled Frog-Leaping Algorithm (ISFLA). The next step could consist of applying the Improved Shuffled Frog-Leaping Algorithm to deal with the dynamic and continuous berth allocation problem and quay crane assignment problem (QCAP) simultaneously. Using ISFLA to deal with the quay crane scheduling problem (QCSP) is another research direction. In addition, a comparison between the Improved Shuffled Frog-Leaping Algorithm and other approaches proposed in past research (such as the [1,34]) could further be conducted. In addition, improving the heuristic used in the second stage of the two-stage procedure is another direction of research.

Conflicts of Interest:
The authors declare no conflict of interest.