Optimising cargo loading and ship scheduling in tidal areas

This paper describes a framework that combines decision theory and stochastic optimisation techniques to address tide routing (i.e. optimisation of cargo loading and ship scheduling decisions in tidal ports and shallow seas). Unlike weather routing, tidal routing has been little investigated so far, especially from the perspective of risk analysis. Considering the journey of a bulk carrier between N ports, a shipping decision model is designed to compute cargo loading and scheduling decisions, given the time series of the sea level point forecasts in these ports. Two procedures based on particle swarm optimisation and Monte Carlo simulations are used to solve the shipping net benefit constrained optimisation problem. The outputs of probabilistic risk minimisation are compared with those of net benefit maximisation, the latter including the possibility of a ‘rule-of-the-thumb’ safety margin. Distributional robustness is discussed as well, with respect to the modelling of sea level residuals. Our technique is assessed on two realistic case studies in British ports. Results show that the decision taking into account the stochastic dimension of sea levels is not only robust in real port and weather conditions, but also closer to optimality than standard practices using a fixed safety margin. Furthermore, it is shown that the proposed technique remains more interesting when sea level variations are artificially increased beyond the extremes of the current residual


Ship scheduling in tidal areas
A ship's draft is the distance between the waterline and the bottom of the hull. It is a fundamental characteristic of a ship and forms a major constraint in terms of scheduling or cargo loading decisions because a poor choice can lead to grounding in tidal areas or shallow waters. Yet the research on ship loading has mostly focused on operations safety and logistic aspects (see for instance a review in Christiansen et al. (2007)). The question of scheduling with time-varying draft was not tackled until recently, when Kelareva and colleagues developed a deterministic procedure to optimise ship scheduling and cargo loading decisions of multiple vessels at a single port (Kelareva, 2011;Kelareva et al., 2012). Their procedure is based on the increasingly popular concept of dynamic under-keel clearance.
The under-keel clearance is the distance between the deepest underwater point of the ship and the seabed. In the traditional static approach, the under-keel clearance is computed as the difference between the water depth (combining channel depth and tide prediction) and the nominal ship draft. The objective is then to maintain an under-keel clearance at least equal to either a given minimum value or a given percentage of the ship draft, depending on port policy. On the contrary, the dynamic approach deducts from the channel depth and predicted tide, not only the nominal draft, but also a number of allowances accounting for the dynamical responses of the hull to its environment (squat, heeling, wave, water density variation), the tidal prediction error and the variability of bathymetry (Galor, 2008). Kelareva (2011) use short-term predictions of the dynamic under-keel clearance provided by the DUKC software (OMC International, 1993, described in Kelareva et al. (2012; O'Brien et al. (2002)). Specifically, from real-time environmental measurements (water depths, wind, waves, current) and ship information (trim, speed, acceleration), the physical responses to the ship moving in a dynamic environment are computed and the dynamic under-keel clearance is estimated. The optimal cargo loading and short term ship scheduling decisions, given this estimation, are then computed. (O'Brien et al., 2002) report two case studies showing the added value of using such a dynamic under-keel clearance approach in port operations, for both shippers (freight savings and increase in export value) and port operators (reduced dredging costs in the long term, increased ship departure/arrival windows and consequently reduced congestion, contribute scientific knowledge to estimation of the minimal under-keel clearance).
Such a solution is based on real-time measurement of the sea state and provides under keel clearance information for the upcoming tide-window only (Kelareva et al., 2012). Being deterministic, safety margins have to be introduced as the under-keel clearance is only estimated a priori. One can ask whether taking into account the stochastic nature of sea levels (and, consequently, the under-keel clearance) could reduce this safety margin to some theoretical minimum -this is one of the aspects investigated in the current paper. Besides, the planification horizon allowed by the procedure described above is relatively short (one tide) -the current work addresses relatively longer time scales.
The work of Kelareva et al. (2012) was extended to a shipping cost optimisation problem for a fleet considering time-varying draft restrictions at waypoints, variable ship speed and cargo loads as well as flow control through busy waterways (Kelareva, 2014). The specific waterway ship scheduling problem was later formulated by Lalla-Ruiz et al. (2016), who integrated tide as a constraint in their approach to optimally schedule the flow of incoming and outgoing ships through different shipping channels (so that the waiting times were globally minimised).
Similarly, researchers focusing on the berth allocation problem, which aims at scheduling berth and crane allocation to optimise port throughput, introduced tide as a constraint only quite recently. While early works (Xu et al., 2012;Du et al., 2015) were more concerned with the quantification of the economic impact of tides on port operations, recent studies developed practical models and solutions for berth scheduling optimisation (Dadashi et al., 2017;Zhen et al., 2017) or quay crane allocation (Yu et al., 2017) in tidal ports.

Shipping optimisation in stochastic environments
Maritime transportation is an activity particularly subject to risk, i.e. the possibility of a loss, due to the complex dynamics and stochastic nature of the multi-dimensional environment in which it takes place. From the weather at sea to port variables (berth availability, loading/unloading works), including the volatility of bunker fuel prices, a range of uncertain factors condition the outputs of a shipping operation. In spite of its significant impacts on shipping productivity, the issue of uncertainty has remained marginal in the research on maritime transportation until recently. Indeed, as stressed by Song and Furman (2013), due to the complexity and intractability of some shipping problems, authors introduce simplifications (constant speed, single cargo type, basic weather model, etc) that are different from one study to another, making comparison difficult. The introduction of stochasticity is often limited to the modelling of a single or a very limited number of factors (e.g. weather (Azaron and Kianfar, 2003), market demand (Chuang et al., 2010), weather and berth occupation (Agra et al., 2015)). Water depth is also an uncertain factor that should not be neglected. Although tide forecasts used to predict the water depths in shallow seas are traditionally given by harmonic analysis from past observations, a range of causes can modulate the observed water levels. These encompass weather influence, river discharge, the interaction between currents, shallow water seabed and ship traffic (NOREL, 2014) and lead to significant deviations between astronomical tides and actual water level observations (called residuals hereafter: the difference between observations and predictions). Flowerdew et al. (2010) estimate that the root mean square error on the high tide predictions in UK tide stations is typically 10 cm and rises to 29 cm for high tidal range stations. Makarynskyy et al. (2004) note that sea level residuals can amount to 30% of the total measured sea level in Hillarys Boat Harbour, Western Australia.
The uncertainty about future water depths has a significant impact on shipping optimisation. First, as shown in the case study presented in Section 2.1, even for a small-sized carrier of horizontal dimensions 85 m ×15 m, one additional centimetre of under-keel clearance can be turned into an extra freight of 13.05 metric tons (mt) whose value ranges from US$ 2, 556 for a single hold of malting barley (Agriculture & Horticulture Development Board, 2017) to US$ 223, 477 for a single hold of tin (Quandl, 2017) with little increase of operational costs in short journeys 1 . Secondly, when it costs thousands of dollars a day to operate the same vessel, missing a tide-window because of a negative anomaly in the water depth is significantly costly to the shipper, to say nothing about the cost of grounding and its potential environmental consequences. Another economic justification is found in O'Brien et al. (2002), who shows how the use of the DUKC software (deterministic dynamic under-keel clearance estimation) allowed 123 vessels to load an additional 743,246 tonnes of coal (an average of 6,042 tonnes per vessel) in the Port of Hay Point, Australia, in the 1996/1997 financial year. Resulting freight savings amounted to approximately US$ 7, 500, 000 and the increase in export earnings summed up to US$ 30, 000, 000.

Robustness in shipping optimisation
In most of the approaches mentioned in Section 1.1, water depths are considered as perfectly predictable variables. When they are not (see (O'Brien et al., 2002) as well as studies on the probabilistic risk assessment of ship grounding in ports (Gucma, 2004;Gucma and Schoeneich, 2008;Abaei et al., 2018)), an allowance, accounting for tide (and possibly bathymetry) prediction error, is introduced. To the knowledge of the authors, the modelling of this source of uncertainty is not discussed in the literature. It is consequently worth investigating the robustness and optimality of such modelling, as the introduction of safety margins always decreases shipping benefit (Kelareva, 2014).
Although Kelareva et al. (2012) introduced a conservative 15-minute departure window for each departure/arrival decision, the authors justified the slack as a way to take into account the inertia of large ships in port operations rather than to account for sea level uncertainties. The large operational costs of ships tend to prevent the shippers from adding significant slack in their schedule (Christiansen et al., 2007), as a ship is only productive when it is sailing. Again, it is worth investigating the robustness of a dynamic under-keel clearance-based shipping optimisation with respect to real port conditions, namely delays.
An original approach to robustness in ship routing and scheduling is finally found in Christiansen and Fagerholt (2002), who introduced the concept of risky arrival. A penalty cost proportional to the risk of a given schedule is integrated in the optimisation procedure of the transportation cost of a fleet. The work of Brown et al. (1997) should also be mentioned as it questions the applicability of mathematical optimisation in a real port context. The authors especially highlight the situation where a small change in the model inputs leads to a radically different optimal solution. The concept of persistence is introduced as a new feature of the optimisation model so that small changes in the input values do not drastically change the nature of the optimal solution.

Objective and contribution
The present work aims at filling a gap in the field of ship routing by explicitly considering and modelling the uncertainty in tide prediction. A robust analysis of cargo loading and ship scheduling decisions in tidal areas is drawn through a realistic case study. The question at hand is: how can we optimise the cargo loading and ship scheduling decisions given imperfect sea level forecasts without foregoing safety?
To this purpose, a maritime shipping decision model is introduced. The model assumes that an industrial operator has sea level forecasts at N ports, at a given time t 0 , over a prediction horizon T . On this basis, the operator has to decide the total amount of a given commodity to load at the first port, and the fraction of this cargo that will be delivered in each of the remaining ports, as well as the estimated arrival and departure times in each port. The deliveries all have to satisfy the constraints of the inventory routing problem, namely to match a given demand in each port. We assume that all ports have unlimited storage capacities.
Our model computes the 'optimal' solution to such a problem by taking into account the uncertainty in actual sea levels w.r.t. tide predictions. However, the reader must note that our framework does not address the uncertainty associated with the under-keel clearance arising from dynamical responses to the sea state (heeling, heaving, squat effect), nor from the bathymetry. We limit our approach to uncertainty about still water levels resulting from deviations to the tide predictions. These additional dynamical sources of uncertainty could still be integrated to a similar approach in order to address the open water problem (see for instance Briggs et al. (2013) for empirical methods to estimate the squat allowance, Quy et al. (2007) to quantify the ship response to waves and Drwięga et al. (2017) to assess the heel components).
In the following, our model allows us to demonstrate the economic potential of a robust underkeel clearance optimisation. Beyond the application to industrial shipping, for which the bulk cargo load is quite flexible, this work wants to raise awareness of the economic potential for small vessels (mini-bulkers), cheap commodities (grains) and small ports strongly affected by tidal effects (i.e. limited dredging). In the current context of transportation greening (Davarzani et al., 2016), we expect this to be an important area for future applications.
The reader must keep in mind that, to clearly demonstrate the potential of the proposed tide routing approach, we deliberately omit the uncertainties associated with weather, as well as berth availability and cargo handling capacity. With the increase of slow steaming practices (Mallidis et al., 2018), weather does not currently represent a serious limitation since generous journey times are planned and ships are no longer expected to be at the maximum of their performance. The authors aim to further study the limitation represented by berth congestion in future works. For many of the small ports that the current work targets, neglecting the aforementioned uncertainties does not represent a significant issue, as the overall ship and cargo flow is not at its full capacity. For the larger and busier ports, this is indeed a question to ask: Is the added value gained from pure tide routing lost in the variability associated with berthing and cargo handling? Or does tide routing helps to smooth port traffic? This is a subject of future work.
Section 1 has introduced the motivations for the investigation of robust cargo loading and scheduling optimisation in tidal areas and outlined the state of the art around this issue. Section 2 presents the case studies and sets up the shipping decision model. In Section 3, the uncertainty on port sea level forecasts is discussed and a robust alternative to the deterministic decision-making process is presented. Section 4 describes the implementation as well as the modelling of sea level uncertainty. Section 5 discusses the results of our approach, compared to standard techniques based on the results from two realistic case studies. Finally the findings are summarised in Section 6 and perspectives are opened.

Case study
To illustrate the approach in this paper, a case study is presented, which gives the reader context for the model development that is detailed later. We consider a farm cooperative that owns a small-size bulk carrier. The company uses it to carry various farm commodities between ports along the British coast, especially along the route: Liverpool-Portsmouth-Lowestoft.
Given a freight unit value of US$ 195.61 per metric ton (for a malting barley freight (Agriculture & Horticulture Development Board, 2017)), 1 cm of additional draft equals an extra freight of 13.05 mt on the first vessel which conveys an extra profit of US$ 2,556 (case study 1, cf. parameters in Table 1). As described before, although a heavier ship will consume more fuel, for small vessels and short sea voyages, it remains much more profitable for the company to increase the overall cargo loading if possible.
We assume that at time t 0 , given the demand constraints on commodity X in the delivery ports, the cooperative has to decide the total amount of X to load in the departure port, as well as departure scheduling in each port. To this purpose, the company uses the long term harmonic tide forecasts as sea level forecasts. Indeed, the more recent and accurate models are not available in all ports. Besides, traditional tide forecasts remain the main source of water level information for many shippers. In the following, we use 'decision' to refer to this set of loading/scheduling decisions and 'benefit' to refer to the net benefit resulting from the implementation of the said decision in actual conditions.
Hence the problem of interest: given the economic, vessel and port parameters (Table 1), given the information on available quays and maintained depths (Table 2), given the tide predictions in all ports 2 : 1. What is the optimal decision, if the harmonic tide forecasts were considered as perfect (from now on called the 'standard approach')?
2. Is this decision robust to actual port and sea level conditions?
3. What is an optimal and robust decision if the uncertainty on tide forecasts is taken into account?
4. What shipping benefit can be guaranteed, given a predefined level of acceptable error, from the robust solution?
5. How do the robust solution and guaranteed benefit depend on the model of the tide residuals?
6. Is the procedure robust to unseen (i.e. extreme) sea level variations?
To answer these questions, a first case study is implemented for t 0 = 13/01/2017 − 07 : 30 : 00 and N = 2 ports. We compute the optimal solution according to our approach and compare it w.r.t. the standard one. We assess its distributional robustness as well. A second larger-scale analysis is then analysed: 175 different t 0 are considered, from July 2016 to December 2016. We compare the performance of our model's decision with that of a standard approach, in terms of daily benefit and robustness across this range of t 0 values.

Model overview
The model used here is a simplified representation of the maritime inventory routing problem. A material is produced at a given rate in a production site (called the loading port) and consumed at other sites (called unloading ports), at specified rates. Given storage capacities in the production and consumption locations, what is the optimal design of routes and fleet schedule that minimises the shipping costs (sailing and port costs) without interrupting any of the production or the consumption in the aforementioned sites? The optimisation is made on an industrial shipping basis. In other words, the shipper owns the material to be shipped and wants to maximise the net benefit of the shipment (the value of the cargo  loaded minus the shipping costs). The fleet consists of a single bulk carrier or general cargo ship and the study is restricted to the N − 1 legs of length l j , j = {1, ..., N − 1} between the loading (departure, p 1 ) and unloading (arrival, p N ) ports with a constant ship speed, v, provided by the ship specifications. From this, the goal is to optimise the decision vector d consisting of the departure time t j and the cargo m j shipped from each port p j given the overall tide predictions available at time t 0 spanning the horizon T in the entrance channels of all ports, given constraints on the demand a j in each port and given constraints from the ship design (carrying capacity), safety at sea (minimum acceptable water under keel), port management (opening times and price bands for port labour). For now, unlimited storage capacities in all ports are assumed. The question of rate of production in the departure port (i.e. offer) is not taken into account. The ship is assumed to be in the departure port at time t 0 with empty tanks and the most recent predictionsX j (t) at the shipper's disposal for the sea levels in all ports p j , over the horizon T . Time is discretized with the time step ∆t (following the precision in the sea level prediction and observation time series). Here and in the following, in order to simplify the notations, t j will be relative to the origin of our time axis t 0 .

Model description
The model takes time series of sea level point-forecasts in all ports of call as inputs. Given contextual parameters regarding the journey, including ship characteristics, freight and port management, generic constraints about acceptable under-keel clearance, latest arrival time and cargo load, demand constraints in delivery ports and, finally, the net return computation rule for a journey, it computes the optimal cargo loading and departure time by means of a particle swarm optimisation (PSO) solver. Table 1 defines the model's input parameters. A few comments and justifications are provided here.

Journey parameters
The operational speed v is assumed to be fixed and constant over the journey (as it is often the case in maritime shipping models). Operational port costs are subject to price bands. Although most often docks and loading / unloading operations are accessible 24 hours a day 7 days a week, the cost of such operations depends on the local port schedule, e.g. midweek vs weekend periods for Liverpool port (Peel Ports Group, 2017a,b). The simple price band framework allows us to simulate a range of situations: night vs days, week days vs weekends, bank holidays. Finally, the safety margin coefficient α in terms of legally required under-keel clearance to use the confined navigation channel of port p is set to 10% of the laden ship draft as this is usual practice at limited speeds (NOREL, 2014). The open sea version would require adding a 30% margin to the dynamical draft.

Sea level input variables
The sea level point predictions in each port are harmonic tide forecasts available online through the British Oceanographic Data Center portal. The time step, ∆t = 15 minutes, sets a minimum bound on the resolution of our departure time solution. In real port conditions, cargo ships cannot be expected to be exactly on time. Reducing this lower bound would consequently not be realistic. As the tide height can change quickly, and because we are dealing with additional centimetres of under keel clearance, it would also not be judicious to increase ∆t too much. Indeed, the sea level within 1 hour (or even 30 minutes) could change significantly with respect to the small variations we are interested in. Consequently, a time step of 15 minutes seems a suitable trade-off.

Model variables
The ship draft, a key element in shipping planning and realisation, is a function of the cargo load as well as the fuel mass f (t) in tanks at the time t of interest. Considering Archimedes' principle and the equilibrium of forces in a gravitational field, the draft r can be estimated from the equality between ship's weight and water displacement. In a simple approximation (barge ship), this leads to: where r 50 is the half laden ship's draft, S the ship's horizontal area, m max its carrying capacity, ρ the water density. The function f (t) is computed by taking into account the fuel consumption rates at sea f s and at port f p respectively, the time already spent at sea and at port respectively at t, as well as the total fuel load necessary to move the ship from one port to another and (un)load material. Dynamical effects such as the squat effect or the heel due to the wind and the wave responses can reduce the under-keel clearance temporarily. They are not taken into account here beyond the safety margins αr(t) as, again, we consider the still water problem.

Constraints
The ship's cargo and scheduling have to satisfy some constraints. First, at any stage, the cargo load m j cannot exceed the tank capacity m max and must fit with the requirements for safe structural behaviour of the hull (m j ≥ m min ), as well as with the demand constraints in the next ports to visit (m j ≥ N k=j+1 a k ). In the following: m min is taken as the minimum of the structural constraint and the economic constraint.
The fuel load necessary to carry the ship and its cargo m j over the distance l = N −1 k=j l k at speed v and load/unload the freight at rate u p in port p must be subtracted from m max : Second, to enter/leave port p j at time t, the water depth must be greater than the ship draft plus the safety margin: Third, the ship cannot leave port p j before the cargo is (un)loaded and must arrive before the horizon T is reached, so:

Shipping return
The problem is to find the optimal combination of decisions d * = (t * j , m * j ), j = {1, ..., N − 1} that maximises the net benefit B, where: The gross value V is the merchant value of the cargo: with C c the unit value of the freight. From there, we subtract the operational costs of the journey, starting from t 0 (time of decision) with an empty ship and finishing at t a + m N u N after unloading the material in port p N where t a is the arrival time in the last port of call. These charges encompass the propulsion costs: where T s is the total time spent at sea and and T p , T p * the total times spent at port p within and outside normal work hours respectively and C f is the fuel unit price. Operational charges also include usage costs: with C u the hourly usage cost (staff) of the ship. Finally, port costs have to be included: where · is a ceiling operator and C p , C bp * , C bp * , the daily port fee, hourly manutention prices in normal hours and outside normal hours in port p respectively. Z is the cost of not making the delivery in time (i.e within the horizon T ). Depending on the aim of the user, Z can also take into account the negative externalities on the environment and society of a grounding (Z → −∞) or simply the loss for the shipper (Z = −V − (O + P + U )).

A PROBABILISTIC APPROACH TO DECISION MAKING
Using the model described above, one can choose an optimisation technique (e.g. particle swarm optimisation or simulated annealing) to compute the optimal decision to take at time t 0 , according to the sea level forecast time seriesX j (t), for the ports p j , j = {1, ..., N }. Such a calculation does not consider the actual stochastic behaviour of the water depth. Local sea levels are influenced by a range of factors, including weather. A residual e j (t) = X j (t) −X j (t) between the predictions and the observations can lead to either a regret (e j > 0: the shipper could have loaded more or departed earlier) or a loss (e j < 0: in order to adjust to the actual water level the journey is delayed, or a grounding can happen). In other words, the resulting solution is risky. It does not tolerate a negative deviation to prediction nor port delays. In order to take account of the uncertainty on the output of a given decision, we must introduce a risk measure.
Risk is a polysemous notion. This is reflected in the many works that have been published in order to identify and classify the variety of definitions (from Kaplan and Garrick (1981) to Aven (2012) and Goerlandt and Montewka (2015)). In the field of maritime transportation, an analysis of risk-related publications spanning over forty years  led by Goerlandt and Montewka (2015) shows that the majority of the works rely on four definitions: (a) Risk is the expected value of the loss; (b) Risk is a combination of scenarios, their probability and the extent of their consequences, represented as a triplet; (c) Risk is the possibility of a loss; or (d) Risk is the probability of an undesired event.
Although simplistic, definition (a) has the advantage of easing comparisons between two options as the information about the possible scenarios and their consequences is synthesized into a single number. In the present work, 'loss' takes the meaning of the loss in profit due to the fact that decision d ∈ D is taken at time t 0 based on imperfect forecastsX j ∈ X of the environment state X j ∈ X . Let F j be the cumulative distribution function over X j , which is conditional on information on the prior values of X j and possible other information. LetF j be a predictive distribution of X j (that is a distribution overX j ) provided by the forecaster at t 0 . LetX j (t) be a point forecast time series of X j (t) over time [t 0 , t 0 + T ], B(., .) : D × X → the utility function (namely the net benefit of the journey based on decision d) and y(·) : X → D an optimal action function defined by: The loss function L(., .) : D × [0, 1] → is then defined by Granger and Machina (2006) as: for allX j , X j ∈ X . In other words, the utility of the decision made under uncertainty B y(F j ), X j is compared to the utility resulting from the decision made under perfect knowledge of the future B y(X j ), X j . The loss associated with a given decision can only be evaluated a posteriori as it requires the knowledge of the exact future states of the environment, that are not known at the time of the decision. Hence the recourse to the expected loss which only requests the actual knowledge on the possible values of these future states. We consequently define the risk R of taking a shipping decision d as: Looking more closely at the definition of the loss which we aim to minimise (the expectation over the space of sea level residuals), one can notice that minimising E L d, X j F j is equivalent to finding the decision d * that maximises the expected benefitB The decision minimising R would be, from a frequentist viewpoint, the one that, on average, over a large number of journeys, produces the maximal net benefit. The theoretical expectation addresses both the feasibility and the performance (high return) of the candidate solution, since the cost −Z → ∞ of a grounding would prevent any solution with the least probability of grounding to be returned as optimal.

IMPLEMENTATION
The problem of deterministic shipping optimisation was defined in Section 2.3.5. It consists of finding the combination of decisions d * = (t * j , m * j ), j = 1, ..., N − 1 maximising the net benefit of the shipping given sea level forecastsX j , j = {1, ..., N }. Similarly, the risk minimisation problem consists of finding a set of decisions maximising the objective function (or risk function) defined in Section 3.
Both are constrained 2(N-1)-dimensional optimisation tasks whose objective functions are not continuous nor differentiable. As a result, classical analytical optimisation techniques cannot be used. Hence the call to derivative-free algorithms such as particle swarm optimisation to estimate d * . A range of other computational methods could have been implemented as well. However, PSO was chosen because it generally demonstrates good convergence and execution speed properties in addition to its simplicity of implementation. A review and comparison of the derivative-free approaches is provided in Rios and Sahinidis (2013). PSO is an iterative stochastic optimisation technique that imitates the natural swarm behaviour of a bird flock (Eberhart and Kennedy, 1995). At each iteration, the elements (particles) of the flock explore the search space in a semi-random way and evaluate the fitness (value of the function to optimise) of their positions. They share the information so that their next move is influenced by both their own findings and the findings of the other members of the swarm. The algorithm stops when the desired number of iterations is reached and the position with optimum fitness is returned. Algorithm 1 describes the procedure and our implementation choices.
Because the risk function defined in Section 3 cannot be written in closed form due to the definition of the net benefit B, it is natural to turn to Monte Carlo simulations to estimate them, within the PSO Algorithm 1 Particle Swarm Optimisation procedure 1. Initialise randomly the position d i of each particle i in the search space D and set their initial velocity vector to 0.

For each step j:
(a) For each particle i: .., N }}) resulting from the shipping decision d i with actual sea levels X j (t), or the expected benefitB(d i )). This is the fitness of position d i . ii. Update the personal best b i (j) of each particle i.e its position with optimal fitness among the set of previous iterations. Similarly, identify and update the global best g(j), that is the best solution among the positions visited by the whole swarm so far. iii. Move each particle according to the following equation of motion: where the velocity is defined by: The cognitive c 1 and social c 2 coefficients are set up so as to optimise the ratio between individual exploitation and social interaction while the linearly decreasing inertia weight ω(j) limits 'velocity explosion'. The diagonal matrices R 1 and R 2 introduce stochasticity in the walk of the particles.
3. Stop when the maximum number of steps is reached or when there is no change in the global optimum for a given number of steps. Return the position with optimal fitness. procedure. Algorithm 2 shows the general approach, now referred to as R P SO . B P SO refers to the "deterministic" optimisation of the shipping benefit (by means of Algorithm 1), that is without taking into account any uncertainty on the sea level forecasts (although technically PSO is a stochastic technique).
Hereafter we name nominal state the forecasted sea-level state.
Algorithm 2 Procedure R P SO 1. Initialise randomly the position d i of each particle i in the search space D and set their initial velocity vector to 0.

For each step s:
(a) For each particle i: i. Sample time series of tide residuals e j (t) in ports j = {1, ..., N } from a distribution model. In this first study, we assume the residuals to be independent from port to port at a given time, and between any two given times at a given port. The spatial independence is not a strong assumption as long as ports are not too close. On the other side, the time independence can be discussed as on short durations the residuals show correlation. ii. Compute the net benefit B(d i , {X j (t), j = {1, ..., N }} for the simulated sea level conditions. These are given by the nominal state modified by the tide residuals, namely at port j: X j (t) =X j (t) + e j (t). iii. Repeat steps 2(a)i to 2(a)ii until the number N s of simulated environments requested to compute the empirical risk function is reached. iv. Estimate the latter from the N s outputs of step 2(a)iii.
(b) Move particles according to the general PSO procedure described in Algorithm 1, in the search space, step by step. Here the objective function to be maximised are the risk functions computed in step 2(a)iv, e.g. the expected benefit or the expected loss.
(c) Stop when the maximal number of steps is reached or when there is no change in the global optimum for a given number of steps and return the position with optimal fitness.

Sampling
Mathematically R is, within a constant, the expectation of the economic output of a given decision when the sea levels in both departure and arrival ports vary around their nominal state (the predicted one). Such a definition implies that R is model-dependent: its accuracy depends on the quality of the modelling of sea level residual distributions. In this section, we present the results of an analysis of these residuals in both departure and arrival ports of our first case study, namely Portsmouth and Liverpool. The dataset used for the modelling and then the testing consists of sea level residuals sampled every 15 minutes between 00:15-01/01/2006 and 23:45-31/12/2016 UTC, in each port. We split it into two parts: even years (dataset D e ) and uneven years (D u ). The former is used for modelling the sea level residuals by means of best-fit distributions. It is then used for the shipping optimisation procedure per se. Finally, D u is used as validation set, to perform simulations and gather statistics on the performance of the optimisation outputs.
Three distributions were tested: Gaussian, Logistic and Gaussian mixture model (GMM). The number of components in the Gaussian mixture models were chosen so as to minimise the Akaike Information Criterion (McLachlan and Peel, 2004). This criterion assesses the goodness-of-fit of a model to a data set while introducing a penalty that increases with the number of free parameters requiring estimation. The aim is to find the optimal trade-off between model complexity and loss of information.
Kolmogorov-Smirnov statistics were computed to quantify the goodness-of-fit of each model. They reject at the 1% significance level the null hypothesis that the residuals follow a Gaussian or Logistic distribution for both ports. A graphical analysis of the three models ( Figure 1) shows that the Gaussian distribution significantly under-represents the small deviations of sea level observations with respect to tide predictions. Hence the introduction of the Gaussian mixture model, that globally represents the original residual distribution with greater fidelity. Besides, the GMM is able to capture the long tails that the single Gaussian or Logistic cannot. This could be important, as extreme events are usually in the tails. This analysis will first be used to assess the distributional robustness of the optimisation procedures in Section 5, by analysing the effect of the residual modelling on the optimisation results. Besides, on a standard desktop computer running Linux, sampling from a Logistic distribution is about 10 times quicker than from a Gaussian distribution and 15 times quicker than from a 5 component mixture distribution. Since feasibility is at stake in operational research, in the following sections we check whether the difference in risk outputs and its implication in real-world decision making justify the added complexity of the GMM input model.

RESULTS & DISCUSSION
Initially, we present the result of a study between N = 2 ports, allowing us to assess the distributional robustness of our probabilistic approach and justify implementation choices for the second study, a larger performance analysis with N = 3 ports.
All the results in terms of benefit B will be expressed as multiples of B 0 = US$ 363, 550 (resp. B 0 = US$ 190, 530 for the second case study). We also set the cost of not making the delivery in time to Z = −V − (O + P + U ). Negative benefits would thus imply a grounding or the impossibility to reach the arrival port within the specified time horizon.

Deterministic approach
The B P SO procedure recommends the ship to leave Portsmouth Harbour at 11 : 45 UTC on January 13th 2017 with an overall barley freight of 4, 475 mt. The precision on these recommendations is estimated to 3 mt in freight and 15 minutes in time as standard deviations (from 1,000 independent runs). Figure 2 presents a mapping of the final shipping benefit over the decision search space D, given the forecast a priori at hand and given perfect forecasts, i.e. the a posteriori exact observations of the sea level depths. The optimal decisions according to B P SO in each scenario differ by one tide cycle  in time and about 400 mt in cargo load. In other words, the deterministic solution under imperfect harmonic predictions is far away from optimality in the real-world of non-zero residuals. Besides, it is quite straightforward to see on these maps that both solutions are very sensitive to perturbations. A 15mn departure/arrival shift or a negative error in the actual sea levels both shift the expected benefit from maximum to the negative area. One way to get over the second limitation is to improve the accuracy of sea level forecasts. This is currently achieved by means of storm surge models. To take into account the local weather perturbations, these models use atmospheric forecasts as forcing in shallow-water hydrodynamic simulations e.g. the CS3 storm surge model covering the sea of the northwest European continental shelf (Flowerdew et al., 2010). Nevertheless, whatever the accuracy reached, these forecasts cannot prevent the issue of port perturbations and delays. Hence it seems reasonable to develop a robust solution instead of a single deterministic optimisation.

Risk model
We now use R P SO to compute the optimal shipping decision under uncertain sea levels. The risk metric presented in Section 3 is combined with one of the three sea level residuals distribution models under consideration. Table 3 reports the statistical results of each combination as regards the optimal cargo load, departure time and the resulting guaranteed benefit at the error rate of 2 %, that is the 2% percentile B .98 . The latter is estimated from 100,000 Monte Carlo simulations. In order to prevent a methodological bias, these simulations sample the sea level by means of bootstrapping (over dataset D u , c.f. Section 4.0.1).
As the purpose of the R P SO procedure is to support decision-making, it is necessary to analyse the

Reduction in standard deviation (%)
Logistic Gaussian GMM Sea level residual modelling Figure 3: Performance of each optimisation approach (sea level residuals distribution) from the perspective of the reduction of the guaranteed benefit at the error level of 2% and the standard deviation of the actual shipping benefit, with respect to the performances of the "deterministic" solution based on sea level forecasts alone. 100,000 Monte Carlo simulations are used to compute these statistics, with bootstrap sampling.
consequences of the above results as regards their translation in terms of practical shipping decision. The overall majority of the computed departure times are located within a 15 minutes time slot centered on 00 : 15. Taking into account the relative inertia of large vessels and generally slow port dynamics (from decision to subsequent actions), this range of uncertainty can be seen as a buffer to consider in the decision-making schedule. Trying to increase the precision on t would be meaningless considering the real world context of a maritime shipping problem. As regards the distribution impact, Logistic sampling produces more conservative loads than the GMM approach and further again, than the Gaussian one. The difference between the maximal and minimal loads above-mentioned is in the range of 25 mt, that is in our case study less than 2 centimetres of draft. This leads to close guaranteed benefits B .98 . Figure 3 summarises most of the information discussed above: a Logistic sampling will produce more stable (smaller variance) outcomes than the other residual models. It also shows that the approach can be said distributionally robust. Indeed, the ranges of the reduction in standard deviation and in guaranteed benefit when the underlying distribution varies are much smaller (close to 0.06 and 0.006 % respectively). Three observations can be highlighted as well. First, in this particular case study, the stochastic optimisation allows the owner to (in most of the configurations) save money as the guaranteed benefit is above the expected benefit of the deterministic decision in real conditions (B = 0.164). Second, the spatial organisation of the points underlines a general pattern in robust optimisation: the guaranteed benefit increases at the cost of the increase in variance (Gotoh et al., 2015). Finally, as noted by Gotoh et al. (2015), the variation in actual benefit is about one order of magnitude smaller than the reduction in its standard deviation. Figure 4 summarises all the above considerations in a 3-dimensional view of the optimisation problem. A map of the expected benefit (whose maximum corresponds to the minimal risk defined in the previous sections) is estimated with bootstrap sampling for each couple (t, m) of the search space, as well as a map of the actual benefit variance on a smaller area of the search space. On top of both maps, are reported the decision suggested by the net benefit optimisation from sea level forecasts, perfect forecasts  (i.e. perfect knowledge of the future) and our optimisation approach with a set of residual modellings.

Summary of results
Concretely, as the owner of the company, you could use the benefit optimisation decision that is based on the deterministic harmonic forecasts, load 4, 475 mt of barley and cast off at 11 : 45. However the outcome of this decision, given the actual observations of sea levels is −2.49B 0 . This is much less desirable than the benefit 2.12B 0 that you could make if you knew the future perfectly and left Portsmouth port at 11:45 with 4, 705 mt on board. Using the stochastic optimisation method developed in this paper, you could load cargo between 4, 200 and 4, 225 mt, raise anchor at 11:45 and get a net benefit from 2.06B 0 to 2.07B 0 .
If these decisions were reported in Figure 2(b) (mapping based on actual sea level conditions), one could notice that a port re-scheduling of up to 1.5 hours (earlier or delay) would not substantially change the benefit, nor a variation (in standard limits) in sea level conditions. Besides, Figure 4 reminds that the variance in the actual benefit is substantially reduced for our solutions, contrary to the variance of the deterministic proposition. In other words, the approach R P SO proposes a robust solution. This is true for any sampling distribution although a Gaussian generally leads to solutions with slightly less predictable economic outcomes. Recalling the questions raised in the motivation of the problem (Section 1), in this case study, our stochastic approach demonstrated to be economically valuable with respect to the standard (deterministic) approach. Besides, a simple Logistic modelling of the residuals is enough to produce quality results, similar to those gained by means of a GMM.
One can note that the cargo load output m can be turned into a safety margin ∆r to be deducted from the maximum draft that would have been allowed given the sea level tide forecasts at hand at t 0 (procedure B P SO ). For future works, it would be interesting to compare ∆r with what a "non-stochastic" commercial software would suggest on a similar problem, so as to assess the quality and potential added value of our model.

Case study 2: 3-port analysis
The first case study was a relatively simple example, chosen to show the potential of a probabilistic approach of tide routing, especially for tide-sensitive ports (Portsmouth in our illustration).
In the following, to provide a more representative analysis, the approach is applied to 175 different decision times t 0 between July 2016 and December 2016, each spaced by at least 24 hours. One additional port is also added to the analysis, with the chosen route: Liverpool-Portsmouth-Lowestoft. Again, we first compute the perfect decision, given a perfect knowledge of future sea level conditions in the three ports. The deterministic and probabilistic optimal decisions given tide predictions are then computed. Note that the probabilistic decisions are, given the results in the previous section, computed using the Logistic modelling of the sea level residuals. Besides, since this approach is more time consuming than a standard deterministic approach, we restrict the computations to 50 t 0 , randomly sampled Det. SM=1m Figure 6: Variation of the net benefit when the residual at each port of transit is artificially modified. We sample them from their original distribution whose mean is shifted to the value indicated along the x-axis. For the probabilistic approach (green), the variation in benefit is measured w.r.t. the unperturbed actual net benefit. For all deterministic approaches, the variation is measured w.r.t. the net benefit of the probabilistic approach resulting from the same perturbed conditions. Q 0.1 and Q 0.9 represent the quantile 0.1 and 0.9 of the set of results ∆B computed over the 50 tests.
to represent the various trends in the whole set of 175 t 0 . In addition to the deterministic decision, we add a deterministic decision taking into account a rule-of-the-thumb safety margin on the sea depth, as it is common practice in the maritime shipping industry. Static safety margins of 1m and 0.5m were both investigated in our experiments. We compare the performances of the three different approaches in terms of net benefit in actual conditions. Figure 5(a) shows that 17% of the journeys cannot be fulfilled in the given horizon (i.e. cannot reach the final port, with the prescribed cargo load, because of low sea levels) if the deterministic decision is used without a safety margin. This score is lowered to zero by using the probabilistic approach or the deterministic one including a safety margin of 1 or 0.5m, confirming the robustness of our approach. Moreover, it is clear from Figure 5(b) that the net benefit in these 'critical' situations is significantly higher when using the probabilistic approach than the safe deterministic ones. This shows that switching from a traditional rule-of-the-thumb static safety margin, often too conservative, to a flexible safety margin provided by the probabilistic approach and taking into account the port calls to come, facilitates time savings and/or increased loading, improving the company's overall net benefit without foregoing safety. Such a solution can be said both robust and near-optimal.

Robustness to extreme sea level variations
One could argue that, as the probabilistic approach is based on the modeling of sea level residuals (itself fitted with archived observations), the results might be sensitive to extreme residuals. To analyse a possible lack of robustness to unseen sea level variations, we use the setting of case study 2. From the scheduling solutions for each of the 50 tests, we artificially modify the residuals in a time-window of ±45 minutes around the departure time in each port of transit. The perturbation procedure is the following: instead of the observed residuals, we sample them from the actual residual distribution (in each port of interest) whose mean is shifted to the 1) minimum residual ever observed; 2) quantile 0.1 of the residual distribution; 3) median of residuals; 4) quantile 0.9 of their distribution; 5) maximum residual ever observed. The impact on the net benefit of the journey is assessed and results are summarised in Figure 6. For the probabilistic approach, we measure the variation in benefit with respect to the unperturbed actual net benefit. For all deterministic approaches, the variation is measured with respect to the net benefit of the probabilistic approach resulting from the same perturbed conditions. Figure 6 shows that the net benefit resulting from the probabilistic approach is not sensitive to the more extreme residuals, whether negative or positive. The approach remains, in all conditions, more attractive regarding the actual benefit than the deterministic approaches with safety margins. Clearly, as a consequence of its conservative nature, the probabilistic approach cannot profit from the windfall effect generated by extreme positive residuals as much as a 0-safety margin approach.

Limitations
As stated in Section 3, the risk measure was chosen because its definition allowed us to address both feasibility and performance in terms of robust optimisation. However in practice, as detailed in the next section, R(d) is estimated from Monte Carlo simulations of the shipping journey subject to various residual scenarios. These Monte Carlo simulations investigate a smaller uncertainty set than a theoretical expectation. The modelling of residuals is indeed based on historical data and potentially not conservative enough, though Section 5.3 showed that in the context of our case study, the approach was robust to unseen previous conditions. Besides, calling Monte Carlo techniques implies that the number of sampled scenarios is limited, which is even more true if real-world applicability (computation time) of the decision-support tool is at stake.
We would like to conclude this section by further justifying one of the assumptions in our model. We chose not to consider the possible restrictions in terms of actual water depth during the loading or unloading steps. For more operational decision-making support, these additional constraints should be integrated. In our case study and generally speaking for small vessels, results are not affected by this simplification. As long as (un)loading rates are high and the loads small, the loading/unloading stages are very limited in time and the increasing ship draft matches the rising tide (which is the only potentially problematical scenario).

CONCLUSION & FUTURE WORKS
This study introduced a decision model for robust cargo loading and ship scheduling in tidal areas. We associated a risk measure to each possible shipping decision. This measure was defined as the expected economic loss of taking the decision in an uncertain environment (sea levels), that is the loss with respect to the net benefit that could have been achieved if actual sea levels were perfectly known in advance. We developed a stochastic approach based on particle swarm optimisation and Monte Carlo simulations to estimate the decision that minimised risk. Results from a Portsmouth-Liverpool case study showed that this solution was both robust and optimal with regard to real port and sea level conditions. We also addressed the question of residual modelling and the resultant issue of distributionally robust optimisation. Thereon, the impact of a change of residual model on the optimisation outputs was negligible in terms of decision-making. A final application to 3 ports confirmed the added-value of our approach compared to standard practices.
While both the probabilistic and classical, 'rule-of-thumb', approaches can be considered robust (to, for example, port delays, forecasting errors etc.), the probabilistic approach was shown to be closer to optimality. Both case studies show the relevance of our approach for tide sensitive ports, small capacity carriers and cheap commodities. Finally, by analysing artificial extreme sea level variations, the robustness of this approach to unseen residuals and its efficiency over existing 'rule-of-the-thumb' practices was demonstrated.
To address the computation time and underconservative historical modelling of residual issues, it would be interesting to define sounder uncertainty sets on which the risk metric would then be applied.
Another avenue of research is finer modelling of the sea level residuals, taking into account the cyclic character of data as well as results in the relevant literature, like those of Horsburgh and Wilson (2007) who noted patterns in the observation of highest weather-based surges. Because positive and negative deviations in tide prediction do not have the same effect on the end-user (that is the shipper), more attention could be given to their respective modelling as well as to the way of treating them through an appropriate asymetrical utility function, beyond what has already been done by focusing on the net benefit.
In practice, it will likely be necessary to make a more complex analysis. Shipping is a multidimensional activity. Loading / unloading a ship and leaving / entering a port require external support. We have analysed the robustness of shipping decisions under uncertain sea levels. However from congestion in waterways to berth availability, crane and tug allocation, a range of uncertain factors should also be included in the analysis of a robust optimal shipping decision.
Similarly, the uncertainty on the exact local water depth was assumed to come from the sea surface: the possibility of a weather-induced deviation to tides. Yet a range of factors can also locally modify in space and time the water depth: currents, sedimentation, vessel traffic for instance. Including the uncertainty on the lower part of the water column, at the sea floor, would consequently be interesting.
We assumed that the total fuel costs did not change significantly on a given journey when the cargo load is slightly increased. Our study would benefit from an analysis of the increase in fuel costs with the added cargo value as a function of the weather and ship characteristics (fuel consumption increasing with bad weather).
Finally, we intend to analyse our tide-routing problem from the perspective of existing weather routing solutions. The specificities of tide routing could be introduced in the dynamic criteria and constraints of such approaches.