Alternative mixed-integer linear programming models of a maritime inventory routing problem

: A single product maritime inventory routing problem is addressed in this paper by exploring the use of continuous and discrete time models. We first present a continuous time model based on time slots for single docks, which is enhanced by reformulating the time assignment constraints. Next, we present a model based on event points to handle parallel docks. A discrete time is also presented based on a single commodity fixed-charge network flow problem (FCNF). All the models are solved for multiple randomly generated instances of different problems to compare their computational efficiency, and to illustrate the solutions obtained. .


Introduction
Maritime transportation is a major mode of transportation covering more than 80% of the world trade by volume (Agra et al, 2013;UNCTAD, 2008). When one actor or cooperating actors in the maritime supply chain have the responsibility of both the transportation of goods and the inventories at the ports, the underlying planning problem is a maritime inventory routing problem (MIRP). Such problems are complex. However, improvements in the fleet utilization and charge/discharge amounts can translate into large cost reductions. This means that there is great potential and need for research in the area of MIRPs (Agra et al, 2013).
Many works related to MIRP have been reported in the literatures. Ronen (1983Ronen ( , 1993 published two reviews on ship scheduling and related areas, where different levels and modes were discussed. Christiansen (2004) reviewed the status and perspectives of ship routing and scheduling in 2004. Andersson et al (2010) surveyed the combined inventory management and routing problem from industrial and modeling aspects, and suggested future research with regard to both further development of the research area and industrial needs. Agra et al (2013) reviewed the advances in MIRP, and studied valid inequalities for a single product MIRP to tighten a discrete time model. Surveys by Andersson et al (2010) and Christiansen and Fagerholt (2009) showed that MIRP has received increasing attention in the last decade. As pointed out by Agra et al (2013), the real MIRPs are of high complexity.
There are many works on modeling of MIRPs. Sherali et al (1999) formulated a mixed-integer programming model based on a discrete time representation for the Kuwait Petroleum Corporation (KPC) problem, and presented an alternative aggregate model that retains the main features of the KPC problem to improve the computational efficiency. Savelsbergh and Song (2008) developed a discrete time model for the inventory routing problem with continuous moves based on the integer multi-commodity flow formulation. Furman et al (2011) developed a mixed-integer programming formulation based on a discrete time representation for vacuum gas oil routing and inventory management. Song and Furman (2013) introduced a flexible modeling framework for inventory routing problem, which can accommodate various practical features using a discrete time representation. Al-Khayyal and Hwang (2007) studied inventory constrained maritime routing and scheduling for multi-commodity liquid bulk. They defined a position ( , ), where denotes a port, and is the arrival number of the vessels at that port within the planning horizon. Then, they formulated the continuous time constraints of different arrivals within one port and different positions related to the same vessel. Li, Karimi and Srinivasan (2010) addressed an inventory service problem in which a chemical MNC uses a fleet of multi-parcel ships with dedicated compartments to move multiple chemicals continuously among its internal and external production and consumption sites.
They presented a MILP model similar to that of Al-Khayyal and Hwang (2007) at an operational level with finer granularity. Siswanto (2011) presented a variation of Al-Khayyal and Hwang's model in which he relaxed the problem to consider an assignment of multiundedicated compartments to products. Bilgen and Ozkarahan (2007) presented a MILP model based on a discrete time representation that integrates blending, loading and transportation decisions simultaneously into one model in order to obtain an optimal solution.
As pointed out by , it is possible to solve small practical instances using commercial solvers. However, it is difficult to solve large scale problems due to the large number of binary variables. A computationally efficient model is critical for solving such problems. Furthermore, a discrete time model is generally much larger than continuous time models. Since most of the works have focused on discrete time models, we explore in this paper the potential of continuous time models.
A single product MIRP problem is studied in this paper. There are two categories of ports: production ports and consumption ports. Each port has an inventory with a certain capacity.
Vessels transport products from production ports to consumption ports to ensure that the inventories of ports neither exceed given capacities, nor lie below given safety levels. The ports have given production/consumption rates and berth capacities. The vessels have given capacities, speeds, preparation times and charge/discharge rates. The objective is to determine the times that vessels visit the ports, and the charged/discharged amounts of vessels at each port in order to minimize total cost while satisfying the inventory constraints.
From the view point of modeling and solution, MIRP is a complex MILP problem. Model efficiency is a key issue, since different models usually exhibit very different computational requirements. We develop different models for continuous and discrete time and compare them by the sizes, solutions and CPU times in 6 problems. The main difference between the problem in this paper and that in Agra et al's paper (2013) is that in our problem, operations (charges/discharges) require preparation times, and the sailing costs involve a fixed part and a variable part that depends on the loads.
The remainder of the paper is organized as follows. In section 2, we describe the maritime inventory routing problem. In section 3, different time presentations are discussed. Then, four models, for continuous time and discrete time, are presented in section 4. The computational results for multiple randomly generated instances of 6 problems are presented in section 5.
Conclusions are finally given in section 6.

Problem statement
(The pictures of the vessels and ports come from http://www.zhongsou.net and http://www.nipic.com respectively) Figure 1 Illustrative structure of maritime route As an illustrative example of the MIRP problem addressed in this paper, we consider as shown in Figure 1, one production port and two consumption ports, which produce or consume products at a certain rate. Vessel V1 departs from the maintenance port (for convenience, we call it here the departure port, denoted as O), visits production port P1 and consumption ports C1, C2 sequentially to transport the product, and finally arrives at the maintenance port (for convenience, we call it here the destination port, denoted as D). Vessel V2 departs from O, visits production port P1 and consumption ports C2 sequentially to transport the product, and finally arrives at D. The task is to schedule the visiting times and charge/discharge amounts of every vessel to reduce transportation and operation cost with certain constraints. For more general problems, there are more ports and vessels, and different vessels sail by different routes.
The MIRP problem can be stated as follows.

Given:
The set of production ports ∈ and consumption ports ∈ with production and consumption rates P and Q . The stock of each port must lie between specified lower and upper bounds, st lo , st . The transportation times tsa between each pair of ports ( , ) are also given.
The set of vessels ∈ with upper bound load la up are given. For each vessel, the route is given by the 0-1 parameter x . Fixed durations for charge/discharge, tfi , and charge/discharge rate rv , are given. Vessels cannot wait at a port for longer than a time limit, TimeWUp . The time horizon, TH is given.

Determine:
The charge/discharge times and amounts of the vessels at ports are to be determined.

Goal:
Minimize the transportation, operation and waiting cost.
Consider the following illustrative problem and its solution with two vessels, one production port and two consumption ports. The profiles of the loads of the two vessels are depicted in Figures 2 and 3. As seen in Figure 2, V1 departs from the departure port O and arrives at the production port P1 at the beginning of the schedule period, and starts to charge. When a certain amount of product is charged, V1 departs from P1 for the first consumption port C1.
After a certain time of sailing, V1 arrives at C1 and waits for operation, then discharges. After that, V1 sails to C2, and discharges at C2 without waiting. Finally, V1 departs from C2 and sails to the destination port D. In Figure 3, a similar scenario for V2 is depicted. V2 visits P1 and C2 only.
In Figures 4 to 6, the profiles of stocks of the three ports are depicted. For P1, the stock increases constantly at a certain rate, except when there is vessel charging. The stock will decrease rapidly when a vessel is charging at P1. For C1 and C2, the stocks decrease with certain rates except when there is vessel discharging.

Time representation
For scheduling problems, there are two kinds of time representations, discrete time and continuous time representations. As pointed out by Floudas et al (2004), the continuous time representation can be classified into two groups, time slots and event points. We first describe the time representations before developing the models.

a) Time slot representation
Assuming first the case of single docks, we define time slots for each port. At every port, the operation of each vessel, including preparation to operate, is assigned a time slot as seen in  For port , time slots are defined with the start and end times as continuous variables. The vessels' continuous times for starting to prepare at and departing from the port (finishing operation) must match with the start and end times of a time slot, respectively. If there are vessels visiting port , then time slots are defined for port . In the illustrative problem, there are two vessels, V1 and V2, visiting ports P1 and C2, so two time slots, k1 and k2, are defined for P1 and C2, respectively, involving start times and end times . For C1, only one time slot is defined. According to Figure 7, the start and end points of the two time slots are assigned to the two vessels respectively, with the start points assigned to the start to prepare times, and the end times to the departure times.
Consider now the time slots for the vessels. Recall the example illustrated in Figures 2-6. In Figure 2, vessel V1 visits all the three ports, so it is assigned three time slots (each slot corresponds to one port), and the three oblique line sections correspond to the slots. In Figure   3, vessel V2 visits two ports, so it is assigned two time slots, and there are two oblique line sections. The horizontal sections indicate that the vessels are sailing in the sea, or staying in some port without operating. In Figures 4 and 6, there are two time slots. In Figure 6, the four turning points correspond to the two couples of start and end points. The actual start points are a bit earlier than the corresponding turning points because there are preparation times for charge/discharge. But the difference cannot be seen clearly since the preparation time (0.5 day) is much smaller. In Figure 4, the end time of the first slot and start time of second slot coincides with each other. This means that the middle turning point corresponds to both the end time of first slot and start time of second slot. In Figure 5, since only one vessel visits Port C1, only one time slot is defined.

b) Event point representation
In order to be able to handle parallel docks, we define event points for each port. Starting to prepare for operation, starting to operate and departing from the port (finishing operation) matches with different event points as seen in Figure 8. The number of the points is three times that of the vessels (e.g. six event points for two vessels). There are similar descriptions on the profiles of the ports stocks and the vessels loads as in the case of the time slots.

Figure 8 Event point representation
In Figure 8, the operation times (including the preparation time) of vessels V1 and V2 at port i overlap each other, that is, the two vessels operate simultaneously during time period [to v2,i , te v1,i ], since the time of starting to prepare for operation of vessel V2 is later than that of V1, and earlier than the time of finishing operation of V1. That is to v1,i < to v2,i and te v1,i > to v1,i . Hence, berth capacity of port i can be greater than one.  The basic idea of the time slot based model is to assign the operations to the corresponding ports. Figure 10 illustrates the basic idea, where x vij is a 0-1 parameter indicating whether vessel visits port and in succession. The 0-1 binary variables are defined to denote whether vessel visits port at time slot . If a vessel arrives at a port from another port, a time slot is assigned to it, and each slot can be assigned to it, but a slot can only be assigned to no more than one vessel, and no more than one slot can be assigned to one vessel. Then, we have the following constraints.

Time slot constraints
The subsequent slot cannot start before the preceding one ends, and the end time of a slot of port is later than the start time.
where indicates the last time slot of port .

Assignment constraints
If vessel visits port , a time slot of port is assigned to vessel . Except for the departure port, vessel visits port means that it must come from a preceding port.
No more than one vessel can operate at any port at the same time.
If a time slot is assigned to a vessel, it means that its start time to prepare for operation at and departing time from the port are forced to the start and end time of the slot, respectively.
Otherwise, if the time slot is not assigned to a vessel, the constraints are ignored.
We usually reformulate the disjunction in (5) into the following big-M inequalities:

Operational time constraints
If a vessel departs from port to port , its arrival time at equals to its departure time at plus the sailing time.
A vessel can only start to prepare for charge/discharge operation after it arrives at a port.
If x vij = 0, equality holds for an optimal solution because waiting incurs in a cost.
The interval between the arrival time and the start time to prepare for operation is the waiting time.
If a vessel visits a port, it departs after a certain operation time, involving the preparation and variable operation time, which depends on the amount of charged/discharged product.
At the departure and destination ports, any vessel departs immediately after it operates.
A vessel cannot wait longer than for a specified time limit.
Any vessel must depart from the destination port before the time horizon ends ≤ , ∀ ∈ (14)

Stock balance constraints
A port's stock at the end of a slot equals to that at the beginning plus/minus the amount of production/consumption and that discharged/charged by the vessel.
A port's stock at the beginning of the subsequent slot equals to that at the end of the preceding slot plus/minus the amount of production/consumption.
A port's stock at the beginning of the first time slot equals to the initial stock plus/minus the production/consumption amount.
A port's stock at the end of time horizon equals to that at the end of the last slot plus/minus the amount of production/consumption The ports' stock must lie between lower and upper bounds. Figure 11 Sketch of vessel's load calculation

Load balance constraints
Generally, a vessel visits more than one port. When a vessel visits a port, its load changes according to the operations. Calculation of a vessel's load involves different ports as seen in Figure 11. The load of V1 increases after visiting P1, and decreases after visiting C1 and C2.
The amounts of increase/decrease depend on the operation times.
If vessel visits ports and successively, its load before operation at port equals to that at port plus/minus the charge/discharge amounts at port .
Since , , and are variables, constraints (26) and (27) The equations in (26) and (27) Loads must lie between lower and upper bounds (the lower bound is generally zero).

Objective function
The objective is to minimize the total cost , including the sailing cost (fixed sailing cost and variable sailing cost dependent on the load), waiting cost, and operation cost.
where , , and are the corresponding cost coefficients. In this way, the continuous time model based on time slots is as follows:

Reformulation of the continuous time model based on time slot
where is an auxiliary continuous variable.
Define new variables = and = , = and = , to reformulate the disjunction in (5) in nonlinear form.

Number of variables and constraints of and
Adding all the above, we have that the number of the special constraints in RS is 2 * |V||I|+2 * ∑ | | +2 * |V| * ∑ | | + 2 * ∑ | |. Compared with (4 * |V| + 1) * ∑ | |, it depends on the parameters which is lager between the two numbers of the equations of S and RS.
S and RS are both time slot based models, for which, if some of the slots overlap with each other, it means that there are parallel docks. However, for instance, the stock constraints (15) and (16)

Continuous time model based on event points
To formulate the problem for parallel docks, we develop a continuous time model based on event points.
First, we define the following 0-1 binary variables as follows.
, vessel starts to prepare at port at time point , vessel finishes preparation and starts to operate at port at time point , vessel finishes operation at port at time point

Time point constraints
Time points keep an increasing order.
where indicates the last time point for port .

Assignment constraints
Event points are assigned to different operation times as follows.
Equations (59)-(61) state that if an event point is assigned to an operation, it equals to the corresponding variable.

Visiting constraints
At any port and at any time, the number of vessels, which have started to prepare, and have not finished operation, cannot exceed its upper bound.

� ��
where Cap indicates the berth capacity of port .
Vessel at port i can only start and finish operation once if it visits that port.
Vessel at departure port can only start to prepare and finish operation once.

Operation time constraints
If a vessel departs from port to port , its arrival time at equals to its departure time at plus the sailing time.
The interval between the arrival time and the starting to prepare time is the waiting time.
If vessel visits port , there is a fixed preparation time If vessel visits port , the variable operation time depends on the charged/discharged amount.
At the departure and destination ports, vessel departs when it starts the operation of charge/discharge.
The waiting time cannot exceed its upper bound.
Vessels must depart from the destination port before the time horizon ends.

Stock balance constraints
A port's stock at the first event point equals to the initial stock plus/minus the production/ consumption amount.
The difference between the stocks of the successive event points is the summation of the production/consumption amount and the charged/discharged amount. If vessel finishes preparation and does not finish operation at port i, the product is charged/discharged at a certain rate. Figure 12 An illustrative profile of a production port The stock at the end of the time horizon equals to that at the last point, plus/minus the production/consumption amount.
Stocks must lie between lower and upper bounds.

Load balance constraints
If vessel visits port and successively, its load before operation at port equals to that at port plus/minus the charge/discharge amount at port .
The charge/discharge time depends on the amount.
The vessel's load cannot exceed its lower and upper bounds

Objective function
The objective function is the same as equation (34) in the time slot based models. Hence, the model is as follows:
The event point based model can handle parallel docks. However, the size of the model becomes much larger due to a large number of 0-1 assignment binary variables. Agra et al. (2013) developed a discrete time model for the maritime inventory routing problem.

Discrete time model
We present in this section a similar discrete time model, but in contrast to Agra et al. (2013) who consider only variable costs for operation and fixed costs for sailing, we consider fixed and variable costs for operation as well as fixed and variable costs for sailing. As shown below, the model corresponds to a fixed-charge network flow (FCNF) problem. Figure 13 Operation at consumption port for vessel in the fixed charge flow network In Figure 13, the squares indicate connections between the stocks at successive time points, and the discharged and consumed amounts in current time interval. The circles indicate the connections between the binary operation variables , , , and . We then have the following constraints.

Routing constraints
Vessel can only sail between ports and if there is a route planned from port to for vessel .
where sa vijt indicates that vessel is sailing between ports i and j at time t.
The number of with value 1 for any ( , , ) equals to the total time vessel sailing between ports and .
where is a parameter which indicates the transportation time between ports and for vessel .
Any vessel must start and finish a schedule.

Fixed charge network flow constraints
A vessel cannot charge/discharge and prepare simultaneously.
Flow conservation I. Vessel departs or operates at time if it operates or prepares to operate at time − 1.
Vessel cannot wait after operation.
Vessel can only depart after operation.

Stock balance constraints
Stock of port at time + 1 equals to that at time minus/plus the charge/discharge amount plus/minus the production/consumption amount.
where is the last time point of the discrete time point.
Any stock must lie between specified lower and upper bounds.

Load balance constraints
Load of vessel at time + 1 equals to that at time plus the charge amount minus the discharge amount.
Load of any vessel at any time must lie between lower and upper bounds.

Upper bounds for waiting time
The waiting time cannot exceed its upper bound

Preparation time for charging/discharging
The number of of vit with value 1 equals to the total preparation time

Objective function
The objective is to minimizing the total cost ′ , including the sailing cost (fixed sailing cost and variable sailing cost dependent on the load), waiting cost, and operation cost.
The discrete time model is as follows Even though the discrete time model is relatively straightforward to formulate, the size can be much larger because of a large number of potential event time interval assignments due to the time discretization. Furthermore, the quality of the solution will be worse if the time is discretized with longer intervals.

Numerical results
To test the four models, we considered multiple randomly generated instances of the six test problems indicated in Table 1. Specifically, Tten feasible instances for each problem were selected among dozens of instances that were generated with random original stocks of the ports. Since the random generation does not always lead to feasible instances, only feasible ones are selected for each problem. We solve the different models on all the instances using CPLEX 12.5.0.0 within 10 hours limit of CPU time. In the discrete time models, the time is discretized uniformly with time intervals of 0.5 days. For the discrete time model as well as for the continuous time models we do not consider the addition of valid inequalities. While these have been developed for the discrete case (see Agra et al., 2013), at present none have been developed for the continuous case. Therefore, we have conducted the comparison without valid inequalities since our major goal is to establish the scope and potential of the continuous time models. The Aall models are solved in GAMS 24.01 using a CPU Inter Xron E3110 @3.0GHz with RAM 8.0Gb.

Problem 1
All the four models for problem 1 are solved for 10 random instances. The statistics and solutions are shown in Table 2. From Table 2, we can see that the size of model E is much larger than those of models S and RS, and that of model D is much larger than that of model E. As a result, the CPU times required by models S and RS are quite small, and those of models E and D are much longer. The optimal objective of model D is 7.1% higher than that of model E, because the time domain is discretized. And as seen in Table 2, models S and RS are tighter than models E and D. The schedule obtained by model S for one of the instances of Problem 1 is illustrative in Figure 12.

Problem 2
All the four models for problem 2 are solved for 10 random instances. The statistics and solutions are shown in Table 3. There are similar trends as in problem 1. The number of nodes for model RS is one third as that of model S. This means RS performs better than S. The required CPU time of model D is ten times as long as that of model E. The optimal objective of D is 15.1% higher than the ones of the other models.

Problem 3
All the four models for problem 3 are solved for 10 random instances. The statistics and solutions are shown in Table 4.

Problem 4
All the 4 models on problem 4 are solved For 10 random instances. The statistics and solutions are shown in Table 5. Problem 4 is a larger scale problem. Note that models E and D cannot be solved within the time limit of 360000 seconds (100 hours), and no feasible solutions can be found with that time limit. In Table 5, we can see that models S and RS are very efficient, with RS being faster than S.

Problems 5 and 6
Problems 5 and 6 involve parallel docks. Therefore, models S and RS do not apply. Tables 6   and 7 show the statistics and solutions of models E and D for 10 random instances. The results show that model D takes longer time than model E for problem 5, while model E takes longer time than model D for problem 6. In both cases model E obtains a lower cost objective than model D. The schedule of one of the instances of problem 5 is shown in Figure   13.

Conclusions
A single product maritime inventory routing problem has been addressed in this paper. Three In summary, for the maritime inventory model, the computational results have shown that continuous time models have the potential of being more efficient than discrete time models.
The time slot based models can be solved effectively for the problem with single docks.
However, for large sized problems with parallel docks, even the continuous time event point model can require long CPU times. Hence, valid inequalities and/or decomposition algorithms may be required to effectively solve these problems.