A robust multiobjective model for the integrated berth and quay crane scheduling problem at seaside container terminals

The ever increasing demand for container transportation has led to the congestion of maritime container terminals in the world. In this work, the two interrelated problems of berth and quay crane scheduling are considered in an integrated multiobjective mathematical model. A special character of this model is that the arrival times of vessels and the failure (working) times of quay cranes are not deterministic and can vary based on some scenarios. Hence, a robust model is devised for the problem having three objectives of minimising the deviations from target berthing locations and times as well as departure delays of all vessels. This robust optimisation seeks to minimise the value of the objectives regarding all the scenarios. An exact solution approach based on the 𝜖-constraint method by the Gurobi software is applied. Moreover, regarding the complexity of the problem, two Simulated Annealing (SA) based metaheuristics, namely a Multi-Objective Simulated Annealing (MOSA) and a Pareto Simulated Annealing (PSA) approach are adapted with a novel solution encoding scheme. The three methods are compared based on some multiobjective metrics and a statistical test. The advantage of the integration of berth and quay crane scheduling is examined as well.

shipping as well as handling costs [47]. As the United Nations conferences in 2018 and 2017 [2] report, the global container port throughput in the world has raised from about 100 million 20-foot equivalent units (TEUs) in 1990 to over 750 million in 2019. Furthermore, the increase in the worldwide containerised trade volume was more than 40% in the decade 2007-2016 [11]. Regarding this considerable growth, container terminals have to deal with larger quantities of throughput and are becoming more congested. Hence, efficient planning and scheduling of operations at a seaside container terminal are getting more significance. According to [46] and [45], operations in a container terminal can be grouped based on the location into: seaside operations, yard operations, and land-side operations. A comprehensive review of Operations Research problems at container terminals is presented by [41] as well as [36].
At the seaside, firstly, arriving vessels have to be assigned to berths, then the Quay Cranes (QCs) are scheduled to serve them. Afterwards, in the yard, the vehicles and transportation equipments are planned to execute the intra-terminal movements of containers, e.g. from the quay to the storage area and vice versa. In the next step, an appropriate storage plan for containers in specific locations in the yard area is needed. Finally, at the land-side, transportation modes from outside such as trucks and trains are scheduled at gates of the terminal to the hinterland.
The most important operational problem at a seaside area is to schedule the existing berth length or slots to serve the arriving vessels. This is known as the berth scheduling problem (BSP) and if it only deals with the assignment of berths to vessels rather than giving berthing schedules, it is called the berth allocation problem (BAP) [11]. The assignment and scheduling of the available quay cranes, which work on the vessels, strongly influences the process time of vessels at the berth. This task is related to another problem known as the quay crane assignment problem (QCAP). Its more complete version is the quay crane scheduling problem (QCSP), which determines the start and finish time of each crane on each vessel.
An important goal is assigning the vessels to berths which are less distant from their desired berthing locations. The desired or target berthing location for a vessel depends mainly on the average distance from the storage sites of its inbound and outbound containers. Therefore, considerable costs and efforts which are related to the container movements can be saved. Likewise, there are target berthing times for vessels which are based on their arrival time at the terminal. Hence, berthing of each vessel with the shortest possible temporal deviation from its target berthing time is another significant objective. If the cumulative deviation from the target berthing times is reduced, the large costs of delay are decreased. A major goal is to reduce the process time of vessels at the berth. This strongly depends on the assignment and scheduling of the available quay cranes, which work on the vessels. By this achievement, the vessels can be served and leave the terminal sooner, thus more available time can be provided for both the terminal and the vessels. These three objectives are considered together in our berth and quay crane scheduling problem.
The significance of the three considered objectives, which may conflict with each other in practice, for a real case becomes evident as we pay attention to the fact that a container terminal must provide the best possible service to a large number of vessels and handle a huge volume of containers within a planning horizon. If it only focuses on the first objective or minimisation of the total deviations from the desired berthing locations, then the vessels must wait for their target locations and it results in more deviation from their target berthing times and larger delays. This decays the second and third objective. Therefore, a rational balance between them is vital. It is an important reason that our work presents a multiobjective optimisation framework of these goals to provide the decision maker with a bunch of high quality non-dominated solutions to choose from.
The BSP and the QCSP are very interrelated; therefore, they can be integrated as one problem called the berth and quay crane scheduling problem (BQCSP). The advantage of this integration is verified in our work by means of some numerical comparisons. A real aspect in such problems at container terminals is the variability of inputs, which has not been taken into account enough in the previous works. Therefore, this work focuses on a version of the BQCSP which includes a number of scenarios for the arrival times of vessels and the availability periods of some quay cranes. The first is because of the fact that some incidents can change the time that vessels can arrive at the terminal and the second is related to the breakdown possibility of the cranes, which is a common issue in the real world.
Our approach to deal with these uncertainties is to develop a robust comprehensive mathematical model with the components of both problems. This novel model includes three objectives of minimising the weighted sum of distances from the desired (target) berthing locations, the weighted sum of temporal deviations from the target berthing times and the weighted sum of the departure (completion time) delays of all vessels. A number of scenarios are considered for two problem inputs which can easily vary in the reality, namely arrival times of vessels and availability periods of quay cranes. All of these scenarios are respected in the constraints of our model, which makes the resulting solutions be robust to assumed changes in the problem inputs.
For the optimisation, a classical method called -constraint is applied to the problem in the first step. The solution process is supported by the Gurobi solver [1]. Regarding the very long computation time or shortcomings of the exact approach as the problem size grows, two multi-objective metaheuristics which work based on the Simulated Annealing (SA) are applied. They are the Multi-Objective Simulated Annelaing (MOSA) and the Pareto Simulated Annealing (PSA). For their application some novel solution encoding and operators are defined.
The organisation of this paper is as follows: Section 2 provides a review of a number of related works. The model is covered in Section 3. Section 4 focuses on the solution methodologies. The computational experiments and the related results are presented and compared in Section 5. Finally, Section 6 gives the conclusions of our research and some directions for the future research.

Related works
A literature review of the two fields of berth allocation and quay crane assignment is provided by [6]. In that survey, the BAP is classified based on some spatial and temporal distinctions which are whether the berth space is discrete or continuous, whether the arrival times of vessels are deterministic or stochastic, and whether the handling times of vessels are deterministic or stochastic.
A Particle Swarm Optimisation (PSO) algorithm is presented in [46] to solve the discrete and dynamic BAP of a realistic size within a reasonable time. Xu et al. [52] considers the limitation by water depth and tidal conditions in the BAP. They model the problem as a parallel-machine scheduling with processing set restrictions and divide the time horizon into two periods, in which the processing sets are different. Both the static and dynamic cases of the problem have been taken into account and efficient heuristics are developed for them.
Zhen and Chang [53] studies the development of a robust schedule for berth allocation with a degree of uncertainty in vessels' arrival and operation time. A bi-objective optimisation model is presented that minimises the cost and maximises the robustness of schedules, which is calculated with regard to the buffer times. A heuristic is developed to solve the large-scale cases. Karafa et al. [20] formulates the problem as a two-objective problem considering stochastic vessel handling times. An evolutionary heuristic and a simulation-based Pareto front pruning algorithm are proposed to deal with it.
Lalla-Ruiz and Voß [27] presents two Partial Optimisation Metaheuristic Under Special Intensification Conditions (POPMUSIC) approaches to solve the BAP. These methods work based on the interoperation between metaheuristics and mathematical programming techniques. It is proved that the approaches are able to provide promising results for this problem. The methods proposed in [25] provide excellent results at that time and the presented modelling improvements allow to stretch the limits of problem sizes to be solved efficiently. The work of [20] investigates the berth allocation problem with stochastic vessel handling times as a bi-objective problem. To solve the resulting problem, an evolutionary algorithm-based heuristic and a simulation-based Pareto front pruning algorithm is proposed.
Legato et al. [28] develops a mathematical programming model at the tactical level and a simulation model at the operational level. Their framework uses a beam search heuristic to obtain a weekly plan at the tactical level, and then a simulated annealing based search process is proposed to adjust allocation decisions at the operational level. At this level, randomness in discharge/loading operations is taken into account and modelled by an eventbased Monte Carlo simulator. The authors perform extensive numerical tests and compare all models from a computational perspective.
Buhrkal et al. [8] reviews and describes three main models for the discrete BAP and enhances the performance of one. A mixed integer linear programming formulation and a heuristic are presented for the BAP in [12]. Tavakkoli-Moghaddam et al. [43] models the quay crane scheduling and assignment problem as mixed-integer programming and proposes a genetic algorithm to cope with some real-sized instances.
The multi-quays berth allocation and crane assignment problem under availability restrictions which may arise due to weather conditions or planned maintenance, is investigated in [26]. They apply a mixed-integer programming model and a set of heuristics based on general variable neighborhood search to solve their problem. The research is inspired by a real case of a bulk port in Morocco. Hu [18] considers vessel's fuel consumption and emissions in the Berth Allocation and Quay Crane Assignment Problem (BAQCAP). They apply a non-linear multi-objective mixed-integer programming model which is, afterwards, converted to a second-order mixed-integer cone programming model to deal with the problem's computational intractability. Additionally, the impact of the number of allocated QCs on port operational cost, vessel's fuel consumption and emission are analysed.
Elwany et al. [15] proposes an integrated heuristics-based solution methodology that tackles a continuous BAQCAP. The proposed approach is claimed to produce high quality solutions to such a problem in a relatively short time suggesting its suitability for a practical use. [38] applies a GRASP-based metaheuristic to the BAQCAP aiming at minimisation of the total waiting time elapsed to serve all vessels. They prove that this metaheuristic reduces the waiting time and increases both the berth utilisation and the throughput of QCs.
Türkogulları et al. [48] formulates the problem as a binary integer linear program that is later extended by incorporating the quay crane scheduling problem as well. Zheng [54] proposes an online model considering a hybrid berth which consists of three adjacent small berths together with five quay cranes. The objective is to minimise the makespan or the maximum completion time of the vessels. They consider two sizes for vessels in their problem. [19] applies mixed integer linear programming including multiple objectives. Small-sized test problems have been solved by the LINGO solver to validate the integrated model. [33] considers the problem with the objective of minimising a weighted sum of the waiting time, deviation from desired location and departing delay, for all the vessels.
The berth allocation problem with the consideration of uncertainty factors, including the arrival and operation times of the calling vessels, is focused in [51]. They formulate a bi-objective robust berth allocation model. The solution methodology is a swarm-based metaheuristic to tackle the robustness of the berth allocation policy and solve the proposed model. Idris and Zainuddin [9] surveys a BAQCAP in general seaport container terminals and proposes a model predictive allocation algorithm and preconditioning methods for solving it.
The paper [17] develops a model of fully fuzzy linear programming (FFLP) for the continuous and dynamic BAQCAP which are solved in a separate way. The arrival time of vessels is assumed to be imprecise, so that vessels can be late or early as much as a threshold allows. A case study consisting of ten vessels has been implemented in a mixed integer programming (MIP) solver. A new mathematical formulation for the BQCSP is introduced in [22] which encompasses all associated operations and constraints. A significant assumption in their problem is the existence of preemption in QC tasks.
Nourmohammadzadeh and Voß [34] is a basis for this current research, which provides a similar model for the BQCSP but uses stochastic optimisation to deal with the uncertain arrival times of vessels. The solution approaches in that work are the -constraints and the PSA, too. Li et al. [29] presents an integrated berth allocation and quay crane assignment model considering preventive quay crane maintenance activities. The two objectives are minimising the total turnaround time of vessels and the total penalty cost of quay crane maintenance earliness and tardiness. An -constraint based two-phase iterative heuristic is used to solve the problem. An integrated BAQCP under uncertainties is focused in [42]. A stochastic optimisation approach is applied. Wang et al. [49] considers a BAQCAP regarding carbon emission taxation by giving a bi-objective integer programming model. Their objective is to minimise the total completion delay of all tasks and the total operating costs for all QCs.
An integrated model consisting of waterway scheduling, berth allocation, and quay crane assignment is proposed in [16]. Nguyen [32] develops a new priority-based schedule construction procedure to generate quay crane schedules. For this goal, two new hybrid evolutionary computation methods based on a genetic algorithm (GA) and genetic programming (GP) are devised. [30] investigates the unidirectional quay crane scheduling problem for a stochastic processing time, which requires that all the quay cranes move in the same direction either from bow to stern, or vice versa, throughout the planning horizon. The paper [50] selects algorithms for the BAP under algorithm runtime limits.
A summary of some related works in the field of berth and quay crane assignment are presented in Table 1. There have been several works with regard to the application of simulated annealing (SA)-based metaheuristics for multiobjective optimisation. A review of the state of the art of multiobjective SA algorithms is presented by [3]. Tekinalp and Karsli [44] develops a multiobjective SA for continuous optimisation problems. This algorithm has an adaptive cooling schedule and uses a population of fitness (objective) functions to accurately generate the Pareto front. Saha and Bandyopadhyay [39] presents a multiobjective clustering technique to optimise two objectives, which uses an optimiation method based on the SA as the underlying optimisation criterion.
An SA algorithm for noisy multiobjective optimisation with continuous decision variables is presented in [31]. A remarkable characteristic of this method is determining the performance of a candidate solution by the estimation of the probabilities that the candidate is dominated by the current non-dominated solutions. Cunha and Marques [13] proposes a multiobjective SA with a new generation scheme and reannealing procedures. This is a trajectory-based algorithm, in which diversified generation strategies are used to define candidate solutions at different stages of the search procedure.

Mathematical model
Our mathematical model is a robust version of the one given in [34] with some modifications. This model is presented in the following three subsections explaining its assumptions, notations and formulations, respectively.

Assumptions
The assumptions used in our robust modelling of the problem consisting of the berth capacity scheduling integrated with the scheduling of available quay cranes to serve vessels, abbreviated as the BQCSP, are explained in this part. Two goals are followed in building our model. One is adaptation with the reality and the other is simplification. The berth, where vessels moor at, is assumed to be continuous and of a determined length. Each vessel has a desirable (target) time which is its arrival time and a desirable (target) berthing location according to the storage locations of its unloading and loading containers. The vessels are of three sizes (small, medium and large). Since the vessels can be of a different importance in the planning, we assign each an importance coefficient which is for simplification only based on its size. To conform with the fact that the arrival time of vessels can be unfixed, ten different scenarios are considered for the target berthing times of vessels. They are according to the strategy that ten values are distributed in equal distances within the interval We randomly choose 25% of the quay cranes as those which are unavailable during some periods. For each chosen quay crane, ten different failure times are determined in equal distance within the whole planning horizon. Each crane can have only one breakdown during the planning horizon. An average is considered for the failure duration of crane k as E(F D k ) and three different values are scattered in equal distance within That is, we have 30 scenarios for each crane which is vulnerable to failure in our robust optimisation. The set of scenarios for the availability of crane k Regarding the scenarios, our robust optimisation model is defined as follows: where n is the number of objectives and S is the set of all scenarios. In other words, the objective function is sought to be minimised while the constraints must be satisfied with respect to all scenarios. So, even the worst scenario is considered in the optimisation and the solutions will be robust with regard to all scenarios. This is in accordance with the definition of robustness presented in [4] as "a robust feasible (r-feasible for short) solution to the robust counterpart of (P ) should satisfy all realisations of the constraint from the uncertainty set U , and a robust optimal (r-optimal for short) solution to (P U ) is an r-feasible solution with the best possible value of the objective".

Notations
Our robust mathematical model contains the following notations: The set of quay cranes available at time t k, l, m ∈ QC

Formulations
Objectives ΔT Three objectives are considered for our BQCSP. Equations (1) and (2) formulate the sum of weighted deviations of all vessels from the target berthing locations and times, respectively. The third objective, (3), is the sum of weighted delays of all vessels from their ideal departure times. The elements used in the objective formulations, i.e. ΔB i , ΔT S i and ΔT F i , are calculated in the related constraints according to the consideration of the whole possible occurrences existing in the set of scenarios.
Constraint (4) implies satisfying the required crane-hours of each vessel. The maximum failure duration of each crane between the starting and finish times regarding all the possible scenarios are deducted from the total time. Constraint (5) sets qq k,i,t = 1 for all t values between sc ki and f c ki if the corresponding assignment and schedule are decided. Constraint (6) is responsible for respecting the maximum number of QCs which are allowed to work on a vessel at the same time.
The deviation from the target berthing location is calculated by relations (7) and (8). Likewise, (9) and (10) calculate the deviation from the target berthing times. Equation (11) states that the departure time of each vessel is equal to the latest finish time of all quay cranes which work on it. Constraint (12) calculates the departure delay according to a desirable process time which is based on the assignment of the maximum possible number of cranes to the vessel (Qmax i ). Constraint (13) is for respecting the vessel's length at the berth and assures that no more than one vessel can moor in each point at the same time. Constraint (14) ensures that if vessel j chronologically berths after i (z ij =1), its berthing time is after the departure time of i. Regarding any pair of vessels i and j , they must not collide with each other. It means that they berth at different locations by respecting their lengths or they berth at different times considering the process time of the former. This is guaranteed by (15). Constraint (16) prevents cranes from crossing each other. Each crane cannot serve more than one vessel at a time which is enforced by (17)-(20). This is a model which considers a set of defined scenarios or possibilities for some inputs of the problem while trying to find compromises for the values of its three objectives. Therefore, the resulting solutions are robust even against the worst scenario.

Solution methodologies
Regarding the fact that in multiobjective optimisation each solution corresponds to more than one objective value, solutions can not be easily sorted based on the objective function value. Therefore, a trade-off between the objectives should be considered and there is a concept called dominance. It is said that in the presence of n objectives, a solution x dominates another solution y, x ≺ d y, if the following conditions hold: 1) z i (x) ≤ z i (y), ∀i ∈ 1, 2, ..., n 2) ∃j ∈ 1, ..., n : z j (x) < z j (y) In multiobjective optimisation, we aim at finding a set of globally non-dominated solutions as the best ones instead of only one optimal solution. These solutions are known as Pareto optimal and the set is called Pareto front.
There are various solution methodologies for multiobjective problems. Some of them are classical methods, whereas some are evolutionary or metaheuristic algorithms. Due to the drawbacks of classical methods such as long computational times or the necessity of numerous runs to obtain a set of non-dominated solutions, metaheuristic multi-objective methods are used as alternative approaches for large-sized instances.
In this research, a practical metaheuristic algorithm called Pareto Simulated Annealing is adapted to be applicable to our BQCSP with regard to the high computational complexity of this problem. The performance of this approach is tested by comparing its results with those of an -constraint method which uses an exact solver.

-Constraint
A promising classical multi-objective optimisation approach is the -constraint method, which is proposed by [10]. Unlike some other classical methods, it does not have any problem in finding solutions on non-convex parts of the Pareto curve. It is important because some remarkable solutions may be located on these parts and our final decision can be within them. In this method, one objective is chosen out of all (n) objectives to be minimised. The other objectives are constrained to be less than or equal to some given values. Mathematically expressed, if f i (x) is chosen for minimisation, we have the following problem P ( i ): Each time that the problem is solved by different values of i , one Pareto optimal solution may be found. Therefore, we need multiple solving attempts to search for a set of non-dominated Pareto solutions. In the case of our BQCSP, each of the three objectives can be selected to be minimised and the other two are restricted by two values.

Multiobjective simulated annealing
In the basic single objective SA, we start with an initial solution (s 0 ) and temperature (T 0 ). Then in each iteration a specific number of neighbourhood searches are done. Each time, the fitness of a neighbouring solution is compared with the fitness of the current solution. In case that the neighbouring solution has a better fitness, the current solution is replaced with it, otherwise this happens according to a probability, which is calculated based on the fitness difference and the current temperature. The formula for the calculation of this probability is P = e −(Znew −Z) T , where Z new and Z are the objective function values of the new and the current solution; T is the current temperature. The smaller the difference and the higher the temperature is, the bigger is the replacement probability. If a specific number of neighbours (NS) have been investigated, the SA starts a new iteration by reducing the temperature based on a plan, for example, dividing the current temperature by a constant value CT . By going forward iteration by iteration, the chance of moving to a new solution reduces. Finally, after meeting a termination criterion, the algorithm stops. This criterion can be reaching a maximum number of iterations or passing a number of consecutive iterations (MQ) without any (or with a very low) improvement in the fitness or exceeding an execution time limit (MET ).
Multiobjective Simulated Annealing (MOSA) is originally presented in [40]. This method uses a probabilistic acceptance rule to increase the probability of accepting nondominated solutions. For a problem with n objectives, a random weight is assigned to each objective considering n i=1 λ i = 1. The acceptance probability of a neighbouring solution s of solution s is calculated as follows: where T is the current temperature. Unlike classical SA approaches, which work with a single solution, the MOSA used in this research is population-based. So, neighbourhood searches are done for each solution of the population. The algorithm terminates whenever one of the two conditions, namely passing a number (MQ) of consecutive iterations without finding any new nondominated solution called unsuccessful iterations or reaching a time limit (MET ), arise. The pseudocode of the MOSA is given in Algorithm 1.

Pareto simulated annealing
This multiple criteria metaheuristic algorithm uses the general concept of the classical single objective Simulated Annealing (SA), proposed by [23], and several specific concepts related to its multiobjective nature.
In the PSA algorithm some new ideas related to the multiobjective aspect are added or replace some concepts of the single objective SA. The specific concepts used in PSA are as follows: -Generating solutions or agents: Initially, a set G of generating solutions with a fixed cardinality |G| = φ is randomly produced to provide a sufficiently large search space for PSA. With each of the generating solutions, a separate weight vector λ s = (λ s 1 , ..., λ s n ) is associated such that λ s i ∈ [0, 1] and n i=1 λ s i = 1. Manipulating these weights provides a good dispersion of the resulting Pareto front and makes the solutions of PSA more representative. The generating solutions are treated as "spy agents" which can work almost independently while exchanging information about their position.
-Aggregation function-based acceptance probabilities: In the PSA, when we move from solution s to its neighbour s due to having multiple objectives to be taken into account and considering the concept of domination, one of these mutually exclusive cases can happen: 1) s dominates or is equivalent to s 2) s is dominated by s 3) s is non-dominated with respect to s. In the first case, the new solution is not worse than the current one; therefore, it should be accepted with the probability one. In the second case, the new solution is worse than the current one (s is considered as potentially Pareto-optimal) and should be accepted with a probability less than one in order to avoid getting trapped in local optima. In the third case, s and s are incomparable (and initially non-dominated). The probability of accepting the new solution s for a problem with n objectives is (21). -Management of generating solutions or repulsion: A degree of repulsion α is determined as a very small positive value close to zero. The weight vector λ s associated with a given agent s is modified in order to increase the probability of moving it away from its closest neighbour agents which is nondominated with respect to s. This is achieved by increasing the weights of the objectives for which s is better thans and decreasing the weights of those objectives for which s is worse thans. This can be stated as follows: In the space of normalised objectives, the metric distance between solutions which is used to determine the closest neighbours from s is: Moreover, we need to normalise λ s i in each iteration to satisfy: n i=1 λ s i = 1 -Updating the set of potentially Pareto optimal solutions: The set of potentially Pareto optimal solutions ρ is empty in the beginning of the algorithm. It is updated each time after a new non-dominated solution is generated. Updating the set of potentially Pareto optimal solutions with solution s requires: (1) adding s to ρ if no solution in ρ dominates s and (2) removing all solutions dominated by s from ρ.
The pseudocode of the PSA used in this work is given as Algorithm 2.

Solution encoding
For encoding a solution of our BQCSP, a matrix is defined including |V | (number of vessels) columns each containing the data of one vessel. The number of rows is 2 + 2|QC|. -If no new non-dominated solution has been found during the whole iteration, set q = q + 1; otherwise set q = 0

19
-Decrease the temperature and calculate T = T CT . 20 end 21 -Output the set ρ the corresponding QC is not assigned to the corresponding vessel and values from 0 to 1 are decoded to times in the planning horizon. This encoding structure is depicted in Fig. 1.

Neighbourhood
We define a neighbour solution (s ) of a solution (s) as the one, which is in nbn cells different from s. Therefore, nbn cells are randomly chosen and their values are changed to other allowable random values.

Constraint-handling
The method for handling the constraints of this problem is to respect them during the solution generation. By this way, firstly, a feasible solution is found among a population of randomly generated solutions. We begin alphabetically from the first variable, consider its value to be unknown and find a feasible interval for it according to the given values of the other variables. Then the corresponding value coming from the encoded structure (matrix) explained above is decoded to a real amount within the found interval. This process is repeated for the next variables until all variables are reassigned with new feasible values.

Normalisation of the objectives
In order to better compare the goal achievement in the three objectives, we should make them of the same scale. Hence, an objective value z j is normalised to znorm j as: where Zmax j and Zmin j are the maximum and minimum value of objective j . Zmax j for the three objectives are calculated as follows: So all the objectives are between 0 and 1.

Computational experiments
In this section, firstly, the generation of the test instances is explained. Then, the metrics used to evaluate our multiobjective methods are introduced, parameter tuning of our PSA is described, and finally, the results and comparisons are presented. A Core(TM) i7 computer with 3.10GHz CPU and 16GB of RAM is used for the experiments of this work. PYTHON is used for programming the models and algorithms. In addition, we apply Gurobi [1] for the exact solution processes required in the -constraint.

Generating test instances
The patterns used for the random generation of the required test instances are as follows: The data of the vessels are based on their size. They are shown in Table 2.

Evaluation metrics
In H V is used the most in the multiobjective literature [37]. A very significant advantage of this metric is that it indicates both the accuracy and diversity. It is a unary metric that measures the size of the objective space covered by the set of non-dominated solutions found by the method under evaluation based on a reference point. Regarding the normalised objectives, [1,1,1] is considered as our reference point for the calculation of H V Table 3.

Parameter tuning
The values of initial parameters have a major effect on the performance of metaheuristics.
Setting the parameters at their best values is one of the main challenges in the application of any metaheuristic. The parameters of our MOSA and PSA are set by the response surface method (RSM) [7], which is a design of experiments approach. The instances that we solve include 10 to 1000 vessels. The tuning is done for each problem size (number of vessels) separately based on a half design due to so many factors. The response factor is H V . An interval is given for each parameter value in the beginning, which is determined according to the experience. These intervals and the parameters' values set by the RSM for the two metaheuristics are shown in Table 3 and Table 4, respectively. The parameter as used in the last column is the number of cells or elements in the solution matrix. Independent from the RSM, the solution time limit (MET) is set to 30 min or 1800s to limit the computational effort.

Results
Using the explained patterns, ten random test instances for each size including 10, 20, 50, 100, 200, 500 and 1000 vessels are built. To implement the -constraint method, the third objective (z 3 ) is considered for minimisation, and z 1 and z 2 are constrained. The combination of eleven 1 and 2 values which are scattered in equal distances within [0,1], which means (0, 0.1, 0.2, 0.3,..., 1), are used to constrain z 1 and z 2 . These necessitate altogether 11 × 11 = 121 runs for the -constraint to find distinct non-dominated solutions. A time limit of 10 minutes is considered for each run. The explained -constraint, MOSA and PSA are applied to all of the instances of any size. Figure 2 depicts the number of non-dominated solutions obtained by the methods. As it is evident, the number of non-dominated solutions that the -constraint method can find reduces by increasing the number of involved vessels. The reason is that the problem becomes more complicated for the exact solver and the optimal solution cannot be found within the time limit that we set for each run. So the number of runs without any new nondominated solution increases. On the other hand, the ability of the metaheuristics to locate Fig. 2 The number of non-dominated solutions found by the methods more non-dominated solutions increases as it deals with larger instances. In comparison of the metaheuristic methods, it is observed that the PSA can provide more non-dominated solutions.
In the following, we aim at comparing the methods based on the multiobjective metrics introduced in Section 5.2. As already mentioned, one alternative approach can be solving the berth scheduling and quay crane scheduling part of the problem separately and adding the results together. For this sake, the PSA is applied to each part and each non-dominated solution of one part in added to all the non-dominated solutions of the other part. Finally, the solutions are accumulated and their non-dominated solutions are found as the results of this approach. Figures 3 and 4 provide a comparative view of the results given by the  -constraint, the MOSA, and the PSA with the separated models based on the two metrics of QM and H V . It is observed except for the smallest size including 10 vessels, whereconstraint is better, the PSA method is superior to the other approaches and it has the largest number of representatives in the set of overall non-dominated solutions. It is worth noting that the methods can have some common solutions in this pool. By analysing the H V values, again with the exception of the smallest cases, the PSA is significantly better in terms of this metric as well. It is interesting that from the instance with 50 vessels, the results of the separated models provide higher H V values in comparison with the -constraint. It is again due to the ability of the P SA, which is used in this approach. Figure 5 illustrates the average execution time of only the PSA and the separated models by PSA because except for the smallest case with 10 vessels, the exact solver cannot find  any optimal solution in any of the -constraint runs for larger instances. The time of the separated models by the PSA consists of the times required for each part added together. The total execution time of all the -constraint runs for the smallest instance is over 21500s, while the average elapsed time of the MOSA, the PSA and the separated models are under 468s, 462s and 536s, respectively, even for the largest instance including 1000 vessels. It is observed that the difference between the integrated and the separated optimisation approach in terms of the execution time increases as the problem size grows. It can be seen that the average execution time of the MOSA and PSA are similar.
In the last part, it is aimed at comparing the three methods by a non-parametric method called the Friedman test with the Bergmann-Hommel post hoc procedure [5]. The choice of this method is because of its non-parameteric nature which does not require the normality of results within the groups. Moreover, it is according to [14] one of the best approaches for multiple comparisons. The multiobjective metrics H V , QM and the execution time are statistically compared based on the whole samples of any size to have a more accurate test. These tests are implemented in R [35]. The p-values of the comparisons are shown in Table 5. As the -constraint method reaches its time limit for instances including 20 vessels and more and enough samples cannot be provided, the corresponding tests are not executed and we do not have the related p-values.
As the near zero p-values show, except for the time comparison of MOSA vs. PSA, there are significant differences in the performance of any pair of the methods in terms of the two evaluation metrics and the execution time.
At the end of this part, we can summarise the obtained results as follows: The application of the PSA for the integrated BQCSP results in clearly better solutions in shorter times compared to using it for the separated BSP and QCSP. In addition, regarding the two multiobjective metrics used, the PSA outperforms the MOSA, while the required execution times of them are similar.

Conclusions
In this research, we presented a robust model for the berth and quay crane scheduling problem taking some possible scenarios for the arrival times of the vessels and the availability periods of quay cranes into account. Solution approaches including an MOSA and a PSA with some novel schemes were proposed and examined and their superiority over a classical method with the -constraint with an exact solver is verified. Besides, it is observed that the PSA outperforms the MOSA in terms of two important multiobjective metrics. Moreover, the advantage of the integration of the berth and quay crane scheduling is shown.
For future research, the extension of the model by the inclusion of some new real elements and circumstances can be a direction. One idea in this respect is the investigation in delay management and whether and to which extent vessels and trucks have to wait for each other in case of delays; see, e.g., [24] for the case of public transportation. Other metaheuristic solution methodologies can be tested, too. The uncertainties in other inputs of the problem as well as application of methods such as fuzzy optimisation can also be investigated.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.