113
1 INTRODUCTION
Maritime transport is the backbone of global trade,
carrying over 80% of its volume [40]. In a cost-
competitive market and additionally tightening
environmental regulations, efficient operations are
more important than ever. One approach to
maximizing efficiency already starts in the planning
stage of operations: Using mathematical optimization.
Formalizing a planning procedure as a mathematical
problem automatizes much of the planning, allows
adaptions to changing circumstances by the click of a
button and offers a systematic approach to finding the
optimal, most resource-efficient plan. Despite its
obvious advantages, such mathematical models are
rarely found in the maritime industry, especially when
it comes to core operations of shipping companies e.g.
ship routing in long and short term, ship loading, crew
scheduling, et cetera. Most of these problems are
combinatorial. Due to harsh scaling behavior, state-of-
the-art classical solvers generally struggle to efficiently
solve practical instances, and thus are usually replaced
by heuristic methods. Such solution approaches tend to
be designed and implemented specifically for each
problem to best fit companies’ individual
requirements.
Quantum computing is an emerging field that
provides a new computation approach compared to
classical computers. It is assumed to have the potential
to improve this scaling behavior [4, 24, 39, 35].
Therefore, it is of high interest to investigate how this
Benchmarking the Maritime Inventory Routing Problem
on a Quantum Annealing-Hybrid System
O. Szal
1
, S. Rubbert
2
& A. Rizvanolli
1
1
Fraunhofer CML, Hamburg, Germany
2
eleQtron GmbH, Siegen, Germany
ABSTRACT: The maritime inventory routing problem (MIRP) is an optimization task aimed to increase the
efficiency of the distribution of bulk products by sea. It combines the routing of a fleet of heterogeneous vessels
between capacitated supplying and demanding ports with the inventory handling at the involved facilities. We
consider a well-studied and general MILP-model variant and introduce modelling adaptations to reduce end-of-
horizon effects. The primary goal is to investigate the capabilities and limitations of current large-scale quantum-
based optimization platforms as a new solution method for MIRPs. We thus benchmark the computational
performance of D-Wave’s quantum-classical hybrid solver on our model by comparing it to results obtained with
CPLEX as a classical state-of-the-art solution method. The test instances cover a range of different parameter
scales, ranging from 2 to 4 ports, fleet size of 2 to 7 vessels and up to 45 discrete time periods. The benchmark
results show that the hybrid system fails to find solutions in the same time as CPLEX for about half the problem
instances. In particular, it struggles to explore tight solution spaces of larger instances. The hybrid solutions that
were found vary in quality, averaging to about 65% to 75% of the classically computed objective values. For
improved results we believe that the problem formulation needs to be changed to a regime better suited for the
hybrid solver, e.g. by incorporating quadratic terms.
http://www.transnav.eu
the International Journal
on Marine Navigation
and Safety of Sea Transportation
Volume 19
Number 1
March 2025
DOI: 10.12716/1001.19.01.14
114
technology can help to solve combinatorial
optimization problems from the shipping industry.
By considering a mathematical model from [34] that
represents a relevant inventory and routing problem
for the shipping industry, this paper aims to:
present a more realistic model by adapting the
model such that end-of-horizon effects are reduced,
provide some first insights into the potentials of
quantum technologies for such problems by setting
up and running calculations with a quantum
annealing-hybrid system and analyzing the results
by comparison with solutions of an exact-solver on
a classical machine.
One relevant planning problem in logistics is the
Maritime Inventoury Routing Problem (MIRP). It
combines vessel routing with inventory management
for efficient sea trade. Due to this combination, such
models facilitate the business of vertically integrated
companies responsible for the transportation of bulk or
liquid goods. There exist various adaptions
characterized by company-specific procedures and
needs which require dedicated formulation and
solution approaches [e.g. 16, 13, 32, 7].
In this paper, we introduce an extension to a
prominent and general variation of MIRP [c.f. 34, 19], a
formulation with discrete-time and arc-flow,
formalized as Mixed-Integer Linear Programming
(MILP). Furthermore, we run calculations on instances
that are inspired by Papageorgiou’s benchmark library
MIRPLIB [34, 3].
At the Fraunhofer Center for Maritime Logistics
and Services we focus our applied research on the
identification of optimization problems that arise in the
maritime industry. Moreover, we formalize the most
relevant problems and provide our partners from the
industry with mathematical models and algorithms
that solve real-world data instances. From the
exchange with some of our partner shipping
companies, we can state that the considered MIRP is in
a similar form relevant to their operations. This
advocates its potential for broad applicability in real-
world scenarios. It features a heterogeneous fleet of
vessels, many-to-many route structure, with both split-
loads and split-deliveries allowed, as well as the
possibility of time-varying parameters. The possibility
for parameters such as supply and demand to change
over time enables mid-to long term tactical planning.
Regardless of the size of the planning horizon, a major
challenge can be its finiteness. Especially when
systematically optimizing, a finite time frame can hide
subsequent consequences of decisions made towards
its end. Those hidden consequences can lead to
opportunity cost. We introduce a preventive measure
against such end-of-horizon effects by adding a reward
based on the final state of the vessels.
Apart from challenges when modelling MIRPs,
solving them is not trivial either. As previously
mentioned, Quantum computing is considered a
technology with the potential to improve
computational scaling. While universal quantum
computers are not yet able to solve problems that
exceed a few bits [4], quantum annealers are
technically more mature. They gain this advantage by
giving up universality. Indeed, as special purpose
machines they are constrained to solving a specific type
of unconstrained optmization problems by a
probabilistic heuristic algorithm. We aim to test such a
device to gain insight into this technology’s potentials
for the MIRP and to identify current challenges. To
achieve this, we run some MIRP instances on the D-
Wave’s quantum annealing based hybrid system.
Previous works have considered the application of
quantum annealing only for simpler problems like the
capacitated vehicle routing problem. In order to
evaluate the quality of these results, we compare them
with the solutions of a state-of-the-art exact solver for
MILPs, CPLEX [15]. Furthermore, we investigate the
impact of the computation time limit for the hybrid
system on the quality of the solution.
The paper is structured as follows. In Section 2 we
present a short literature review on MIRPs, modelling
and solution methods. The problem is introduced in
detail in Section 3. Section 4 contains the mathematical
model of the MIRP problem including the
mathematical formulation which adresses the end-of-
horizon effects. Section 5 presents the numerical
experiments and the analysis of the results.
Conclusions and an outlook can be found in Section 6.
2 LITERATURE REVIEW
Due to the novelty of the topic of quantum computing
in operations research, only a few works have been
conducted which address the solution of related
problems by means of quantum computing. To better
contextualize our work in both areas, we focus on
briefly explaining the relevant notions that divide the
MIRP literature, along with pointing out related work,
including quantum. For those interested in more
comprehensive literature surveys regarding MIRP
studies, we refer to [34, Section 2] for a categorization
of the material published prior to the last decade, as
well as to [19] for a detailed review of the last decade.
Well cited examples of studies on the general Inventory
Routing literature include [9] and [14].
The reader interested in the latest development
related to quantum optimization in general should
consult [4].
2.1 Problem and modelling variations
There are multiple variations of the MIRP established
in the literature. They can be classified along multiple
dimensions, which include:
length of the planning horizon (operational, tactical
or strategic),
single or multiple products (in dedicated or
undedicated compartments),
route structure,
cost model,
transshipment [see e.g. 27],
vessel speed optimization [see e.g. 17, 18],
uncertainty [see e.g. 6, 12, 37].
The route structure and cost model typically
depend on the mode of ship operation considered (e.g.
tramp or industrial). Out of the variations above,
especially MIRP applications with uncertainty are
gaining popularity [26]. The basic idea is to not search
for the optimal solution of a MIRP, but for a good
solution that is stable under pertubations in the
parameters. Such perturbations can range from
115
uncertain vessel- or berth availability and delays in
voyages to varying production and consumption rates.
Given a problem formulation of an MIRP, a variety
of mathematical models can be chosen for
representation. Such models have to balance accuracy
and computational complexity. For instance, the notion
of time can be implemented with discrete or
continuous variables, the latter usually being easier to
solve. However, continuous time is primarily used on
the operational level due to difficult implementation of
parameters like supply and demand changing over
longer time scales. While Mixed-Integer Linear
Programs have become the standard in the MIRP
community, some few studies include non-linearities,
e.g. to incorporate variable vessel speed [17, 18] or
special restrictions, e.g. by tides [36].
[19] pointed out that MIRPs, as a real world
application, should be posed as open ended problems.
But most models assume a finite planning horizon, to
be interpreted as a finite section of that, and thus suffer
undesirable end-of-horizon effects. Their importance
was highlighted and discussed by [10], who compared
the performance of four Inventory Routing Problem
(IRP) formulations with regard to such effects.
However, the considered IRPs are land-based, i.e.
feature several differences and are not directly
applicable to the maritime setting. Documented
attempts to tackle the end-of-horizon effects in MIRPs
include [8], who proposed a single-vehicle continuous-
time model with additional inventory bounding
contraints at the end of the planning horizon. In
particular, such hard bounds are supposed to adress
unrealized supply and unmet demand near the end of
horizon and hinder infeasibilities arising beyond it.
More stringent measures were employed by [41], who
removed the end of horizon completely by restricting
to cyclic solutions of a continuous-time MIRP. While
only applicable under certain (periodic) conditions and
less optimal than non- cyclic solutions on one cycle
period, we consider it a promising method. After all,
these solutions have the potential to be more efficient
when repeated over long time periods.
2.2 Solution methods
As a combinatorial optimization problem, MIRP is
commonly formalized as an MILP. But due to
undesirable scaling behavior, real-world-relevant
instances typically grow out of scope to solve exactly,
and therefore require the use of heuristic methods.
Especially common are heuristics which produce
smaller MILPs by decomposing or restricting the
original problem. For instance, rolling horizon
heuristics rely on consecutively optimizing over small
sections of the planning horizon, and are often favored
on strategic problems like anual delivery plans [e.g.
27]. A variety of heuristics were compared by [5] and
[31], with the former study suggesting that combining
multiple heuristics yields the best results on practical
problem sizes.
Other solution algorithms which could provide
similar benefit are quantum heuristic methods. Despite
the current limitations of quantum hardware, there
already exists some related work. For instance, [22]
investigated the possibility to cast IRPs as Quadratic
Binary Unconstrained Optimization (QUBO)
formulations, and simulated instances running on
different quantum algorithms. Some of the considered
problem formulations can be interpreted as a
homogenoeous fleet MIRP with port inventory bounds
transcribed by fixed time windows for delivery and
pickup of orders of size below vessel capacity. For
numerical tests, full loads and non-split-deliveries
were assumed. Cases with no possible direct QUBO
formulation could still be efficiently reduced to QUBO
by use of the ADMM decomposition heuristic [21]. In a
prior study, [20] proposed a hybrid solution for
another IRP-variant, the capacitated vehicle routing
problem (CVRP), employing an actual quantum
annealer in the process. In particular, the approach is
based on a two-phase heuristic com- prised of
clustering and routing, decomposing the CVRP into a
Knapsack and Traveling Salesman Problem,
respectively. The former was solved classically, the
latter with a D-Wave annealer. In the same year, [25]
proposed an alternative QUBO formulation of the
CVRP with various additional timing constraints. The
model was validated by running small size instances
on the same annealer.
We are not aware of attempts to solve an MIRP on a
quantum computer, yet.
2.3 D-Wave annealing-based hybrid system
Quantum computers have the potential to solve hard
mathematical problems, including combinatorial
optimization, with better scaling behavior than
classical computers [30, 4]. There are various
approaches to building quantum computers. This does
not only include different hardware platforms for
qubits, but also different calculation models. Most
prominently, there are
universal gate-based quantum computers [30],
measurement-based quantum computers [11],
and annealers [23].
In opposition to the first two, annealers are special-
purpose machines which can only run the annealing
algorithm [23], a probabilistic heuristic method to solve
QUBO problems.
D-Wave is a pioneer of quantum annealers and has
commercialized its first quantum annealer in 2011. By
now, they offer access the quantum annealer
”Advantage” with roughly 5000 qubits [28].
Additionally, D-Wave developed their quantum
annealers into quantum-classical hybrid systems, so
each annealer is paired with some high performance
classical resources. The combination of both allows to
run hybrid algorithms that solve more generic
problems than QUBOs with just 5000 binary variables.
It can solve larger problems, it can handle other
variable types than just binaries and it even supports
constraints. The hybrid algorithms that D-Wave uses
are not public, therefore we can analyze their results
but not give detailed explanations for their
performance. Importantly, these algorithms are
samplers, i.e, return samples from the solution space
for a given calculation time. While these samples are
biased towards optimization, there is no guarantee that
the global optimum is found.
116
3 MIRP: PROBLEM DESCRIPTION
In this paper we consider a tactical single-product
MIRP with many-to-many route structure and
inventory tracking at every port. It is based on problem
formulations previously introduced in the literature
[c.f. 38, 34].
Specifically, the problem emerges from the
following initial setup: We are given a number of ports
with an increasing or declining inventory level of a
single fungible product. This can be interpreted as
supply or demand of the product. In order for the
inventories to not run full or empty, a heterogeneous
fleet of vessels is employed for the transportation of the
product. Any vessel is able to travel between any two
ports and load or unload a portion of its freight to
contribute to their inventory and manage the supply
and demand. Figure 1 sketches this process.
Figure 1. Illustration of the Maritime Inventory Routing
Problem.
Now that we have introduced the general concept,
let us take a look at the parameters which define our
MIRP formulation in detail.
Port information:
port type (production or consumption, i.e. supply
or demand)
storage capacity
initial inventory level
range of supply or demand
range of loading rates
number of loading berths (in ports there are
unlimited ships allowed but only few can load at a
time)
product price
berth fees
Vessel information:
travel time between any two ports (assumed fixed)
travel cost between any two ports
initial inventory level
storage capacity
starting point
The aim of the MIRP is to coordinate the vessels in
such a way that the trading profit is maximized within
a finite time frame. This coordination must comply
with a range of constraints, enforcing the travel times,
loading rates, capacities and the number of berths in
the ports.
This generic MIRP formulation sets a highly
adaptable baseline model. Variations of it appear e.g.
with fixed production and consumption rates [19] or
with full loads and unloads and soft inventory bounds
[33, 31, 29]. Particularly due to its flexibility, a variety
of company-specific requirements can be
implemented. In MIRP applications the optimization is
usually conducted over a finite planning horizon. An
arising issue of basic profit maximization models is
that the reward is realized at the end of several steps,
including loading, travelling and unloading of the
vessels. If a trade cannot be transacted during the
optimization time frame, there is missing incentive for
the earlier steps. Not taking these earlier steps can lead
to opportunity cost beyond the planning horizon,
especially when a consumption port is running out of
supply. In our mathematical model of the MIRP, we
add a term to counteract some of these effects, which
are labeled in the literature as end-of-horizon effects
[10, 19]. Effectively we are rewarding trades before
they are finalized.
Note that our approach only works in the case of
profit maximization. If instead costs are minimized,
our addition to the model is trivial, as then the
additional term vanishes.
4 MATHEMATICAL PROBLEM FORMULATION
Our formalization of the maritime inventory routing
problem is, to the best of our understanding, an
extension to the model proposed by [34]. For
readability purposes, we reintroduce the model as a
whole. Readers who are sufficiently familiar with the
base model can skip ahead to Section 4.2. Our extension
introduces an expectation value for the revenue
generated by unloading non-empty vessels outside the
optimization time frame, in order to mitigate end-of-
horizon effects. Before introducing equations in detail,
we first give a qualitative overview of the model with
references to the corresponding equations in Section
4.1 and 4.2.
Objective (to be maximized) (14):
+ product sales
travel cost
berth cost
+ expected vessel revenue beyond the planning
horizon
Bounding constraints:
port inventories are bounded (8)
vessel inventories are bounded (9)
berth places are limited (10)
port production and consumption rates are
bounded (7)
loading rates are bounded (6)
Consistency contraints:
every vessel follows a realizable path (2) port
inventories are updated (3)
vessel inventories are updated (4)
a vessel can only occupy a berth if it has travelled to
the corresponding port (5)
a vessel can only exchange product with a port if it
is at a berth of that port (6)
the expected vessel revenue is only added for
vessels which:
are still in the system at the end of horizon (15,
16)
are expected to make revenue with the next
unload to a consumption port after horizon end
(17, 18)
117
4.1 Base model
As a first step to a formal model definition, let us
introduce some sets, variables, parameters and finally
the objective and constraints.
The sets N
and A
v
were visualized by [38] as
graphs, with the nodes N
and edges A
v
. While the
regular nodes (p, t) N simply represent a port p at the
end of a given time period t, the source and sink nodes
SN only augment the model as well defined start- and
end points. This visualization clarifies what a
consistent path of a vessel is: A subset of arcs and
nodes, connecting source and sink, with each regular
node connected to exactly one following node and one
preceding node, while source and sink only have one
connection each. Figure 2 shows an example of the
described graph for two ports, one vessel, and a time
frame of seven periods.
Figure 2. Time-space graph for seven time periods T = {1, ...,
7}, two ports P = {0, 1} and one vessel V = {0}.
Additional parameters visible in the picture are a
vessel starting node of (0, 2) and a travel time of 2
periods. The set of arcs A
0
is visualized in orange and
divided into regular arcs and source and sink arcs. As
an example of the following- and preceding nodes sets,
FN
0
(1,5) and PN
0
(1,5) are plotted in green and red,
respectively, as well as the node (1, 5) in blue.
At first glance a cost model that only includes
travel-, product- and berth cost might seem
oversimplified, but other costs such as port fees, vessel
operation costs or (un)loading costs can be absorbed
into those parameters.
Now that we have defined the travel time and
starting nodes, we can formalize the definition of the
following nodes set:
( )
( )
,
,sink
: , sink ,
v
vv
n p q m
sn n source
FN q s N s t TT n p t FN for some node m of earlier time
else
=
= = + =
Base model objective function:
(1)
Note that since the load l is defined to be positive,
we allow negative values of the product price PP to
represent buying goods at a production port.
Base model constraints:
Let us now state the constraints:
Arc flow:
( ) ( )
*
,,
1
0,
1
vv
nn
vv
n m m n
m FN m PN
n source
ta ta n N n N v V
n sink

+=
=
−=

(2)
Setting the number of active following nodes and
preceding nodes equal for each regular node prohibits
vessels to be spontaneously created or vanish. Due to
the first and third case in the above equation, the route
begins in source and ends in sink. In short, this
constraint forces consistent paths.
Port inventory update:
( ) ( ) ( ) ( )
( )
, , 1 , ,
,
v
p
p t p t p t p t
vV
pi pi PT pr l p t N

= +



(3)
This constraint ensures that the port inventory pi is
correctly updated, considering both
production/consumption and loading.
Vessel inventory update:
( )
1
,
,
v v v
t t p
pt
pP
vi vi PT l t T v V
= +
(4)
This constraint ensures that the vessel inventory vi
is correctly updated.
Berth condition:
( )
,
,
v
n
vv
n
pt
m PN
b ta n N v V
(5)
This guarantees that a vessel v is only able to occupy
a berth at node n if it is actually located there.
118
Loading condition and loading bounds:
min max
,
v v v
n n n n n
LR b l LR b n N v V
(6)
This constraint limits loading speeds and
completely prohibits loading if the vessel is not at a
berth.
Production/consumption bounds:
min max
n n n
PR pr PR n N
(7)
The amount of product pr produced or consumed at
a port within one period is bounded by the minimum
and maximum production/consumption rates PR
min
and PR
max
.
Port inventory bounds:
nn
pi PC n N
(8)
The inventory pi of a port can never exceed its
capacity PC.
Vessel inventory bounds:
,
vv
t
vi VC t T v V
(9)
The inventory vi of a vessel can never exceed its
capacity VC.
Berth limit:
( )
,
v
np
vV
b B n p t N
(10)
The number of vessels occupying a berth at a port
p P is limited by the berth count Bp.
4.2 Extended model - Damping end-of-horizon effects
End-of-horizon effects are myopic operational
decisions in the solution of optimization problems that
occur due to the finiteness assumption of the planning
horizon. In MIRP applications they manifest as an
aversion to travelling and loading near the end of
horizon, simply because these actions cannot be
rewarded in the remaining time. The resulting increase
in unrealized supply and unmet demand can be
directly adressed by imposing hard or soft bounds on
the final port inventories. But while hard bounds might
induce infeasibilities, soft bounds can be difficult to
motivate practically. Furthermore, our choice to
maximize revenue already helps to mitigate low
running consumption port inventories, although for
some instances (e.g. with economically unfavorable
ports) that problem still may occur.
Our approach is to give an incentive for vessels to
load cargo, even if there is no time left to discharge at a
consumption port. In particular, we choose to reward
filled vessel inventories at the end of the planning
horizon by adding an estimated future revenue to the
objective function, assuming that they will unload all
of their freight with the next visit to a consumption
port.
To describe the extended model structure, we need
to introduce the set of consumption ports and some
new variables:
:1
p
Set of consmption ports CP p P PT= =
While not formally a variable of the model, we will
frequently refer to the expression
( )
( )
( )
( )
( )
( )
( )
( )
,
,
,
,,
max
, , ,
,
1
:
qp
v
T
pT
v
v
vv
p CP
pT
q T sink
q T p T TT
qP
pT
vi PP
er
VC
CP
ta TC BF
LR
+




=







(11)
as the expected revenue. For this to be well defined, the
travel cost tensor has to be extended past the horizon.
We estimate it to be the same cost as travelling just
before the end of horizon, i.e. TC
v
((q,|T|),(p,|T|+TTq,p)):=
TC
v
((q,|T|-TTq,p),(p,|T|).
Above we have denoted the upper and lower
bounds of the expected revenue er
v
by
( )
( )
( )
,
,
max
,
:
pT
v
pT
v
p CP
pT
BF
VC
U PP
CP
LR


=−


(12)
( )
( )
( )
( )
,
,
, , ,
1
: max
qp
vv
qP
pT
q T p T TT
p CP p CP
L BF TC
CP
+





= +






(13)
Finally, we can specify the objective function and
constraints of our general model.
Extended model objective function:
( )
( )
, , , :
v
v v v v v
n n n n a a
v V n N
aA
CE ta b l epr PP l BF b TC ta epr


= +



(14)
The constraints of the extended model are the
linearization of the variable definition
( )
( )
, ,sink
: max ,0
v v v
qT
qP
epr ta er


=


The first factor is a boolean that indicates whether
vessel v has gone into sink before the end of horizon.
Additional constraints:
Vessel end-of-horizon existence:
( )
( )
, ,sink
0
v v v
qT
qP
epr ta U v V




(15)
( )
( )
, ,sink
01
v v v v
qT
qP
mer epr ta U v V




(16)
These constraints are equivalent to
( )
( )
, ,sink
:
v v v
qT
qP
epr ta mer


=


119
The reward is not considered for vessels that have gone
into sink before the end of horizon.
Positiveness of expected vessel revenue:
0
v v v
mer erp U v V
(17)
( )( )
1
v v v v v v
er mer er U L erp v V +
(18)
These constraints are equivalent to mer
v
:=max{e
rv
,0},
ultimately making sure that we only reward vessels
that would make revenue by unloading their cargo
after the end of horizon.
5 BENCHMARKS
Our goal is to investigate the capabilities of quantum
annealing on an operational problem of significant
practical relevance. To this end, we compare the
performance of a D-Wave quantum annealing-hybrid
system vs. a laptop with CPLEX on instances of the
above MIRP formulation. In the following we explain
the initiated steps and present the results.
At this point we only compare the generic solvers
and do not use MIRP-specific algorithms or heuristics.
Thus, we require relatively small problem instances for
benchmarking. The most important parameters
deciding the size of any MIRP instance are the number
of vessels, the number of ports and the time frame.
Additionally, travel times effectively act as a scale for
the time frame. We used Pyomo to create 30 instances
covering different sizes and parameters, and saved
them as LP-files. The LP-format serves as a common
ground that both CPLEX and D-Wave are able to read
(with only minor differences in their conventions).
Note that the D-Wave software stack supports the
”constrained quadratic model”-class [1], which
supports integers, floats and constraints. As the name
suggests, it also supports quadratic terms, which we do
not employ in this work.
We ran all 30 problem instances as 5, 15 and 30
second calculations on both the D-Wave system and on
a laptop with CPLEX. Additionally, we used CPLEX to
calculate the global optimum for those instances where
this was possible in less than 24 hours.
As an overview of the results, Figure 3 displays the
ratios of objective values of the solutions found by the
D-Wave solver and CPLEX. These ratios are exactly the
relative solution quality of the annealing-hybrid
system compared to CPLEX. The instances are sorted
along the x-axis by the number of variables. For a
detailed description of all instances see the Appendix,
including the results data in Table 2.
We find that for the given linear problem instances,
CPLEX performs better than the D-Wave solver on
every instance and calculation time. While the hybrid
solver in the best case reaches up to 89.03% of the
optimal solution, CPLEX finds the optimal solution at
least for 9 instances in the limit of 5 seconds, 16
instances in 15 seconds and 18 instances in 30 seconds.
Moreover, on the instances for which the hybrid solver
finds a feasible solution, the worst performance of
CPLEX, that is 88.72%, is close to the best performance
of the hybrid solver.
We find no trivial correlation between the number
of variables and relative solution quality. If there is any
such correlation, we can mainly observe a clustering of
feasible D-Wave solutions at just above 1000 variables,
but this is also where most instances lie. Interestingly,
we see no improvement in increasing the time limit, as
the hybrid solver found even less feasible solutions in
15 seconds than in 5. And while the number of found
feasible solutions went up to 19 for a limit of 30
seconds, the solution quality decreased.
Figure 3. Ratio of the objective values of solutions found by
the D-Wave solver and CPLEX.
Objective values of inf represent unfeasible
solutions (CPLEX gave feasible solutions for every but
one instance; the corresponding marker is excluded).
The x-axis runs logarithmically in the number of
variables of the instances. Included in the plot are the
mean values and standard deviations of all finite ratios.
The lack of an improvement with longer calculation
times suggests a limit to viable problem sizes of the
underlying linear formulation. Such a limit would
oppose any scaling advantage.
Among the 30 problem instances, there are 6 triplets
of instances, which only differ by a lower bound prmin
on the production/consumption of all ports:
( )
min max
min
min ,
nn
PR pr PR n N=
where
max
n
PR
is set to a
fixed value for each port. We adjust prmin to vary the
number of feasible solutions. The smaller prmin, the
larger the feasible-solution space. In Table 1 we take a
closer look at the benchmarking results for these
instances.
Table 1. Relative solution quality for triplet instances
triplet
vars
prmin
D-Wave obj/CPLEX obj
5 s
15 s
30 s
T16_V3_P4_tt3
1135
80
20
0.737
0.662
0.84
0
0.691
0.745
0.71
T16_V3_P4_tt2
1188
99
0.6
20
0.721
0.766
0.717
0
0.689
0.716
0.71
T16_V4_P4_tt4
1350
99
20
0.72
0.772
0.661
0
0.706
0.728
0.718
T15_V5_P4_tt5
1453
99
0.89
0.816
20
0.735
0.763
0.452
0
0.792
0.821
0.826
T20_V3_P4_tt3
1515
85
20
0.684
0.664
0.668
0
0.718
0.684
0.593
T45_V7_P4_tt9
7688
26
10
0.649
0.139
0
0.457
0.659
0.194
Comparison of the relative solution quality of D-
Wave for different sizes of solution space (growing
with decreasing prmin). The triplet name indicates
parameters |T|,|V |,|P| and travel times for each
120
triplet. The vars column lists the number of variables in
those instances.
We observe that as prmin decreases, the hybrid solver
is able to identify an increasing numbers of feasible
solutions. So while it is hard to tell whether tighter
constraints support or hinder branch-and-bound
solvers like CPLEX, larger feasible-solution spaces do
seem to help the D-Wave hybrid solver. More
precisely, as a probabilistic machine, the annealer is
more likely to find a feasible solution if there are more.
However, we do not observe an improvement in the
relative solution quality.
6 CONCLUSION AND OUTLOOK
In this paper we have reintroduced a broadly
applicable variant of the Maritime Inventory Routing
Problem with the capability of being specialized to a
wide array of potential needs and restrictions. To
propose an adaption that is more realistic in the setting
of open-ended operational procedure, we have
introduced a reward based on the final vessel
inventory. Our approach aims to damp end-of-horizon
effects, while still maintaining generality and flexibility
of the model.
Being interested in the computational potential of
quantum computing for the MIRP, we in addition
benchmarked a CPLEX solver (on a laptop) against a
D-Wave quantum annealing-hybrid system on 30
MIRP problem instances of MILP type. In our setup
CPLEX delivered better solutions on every single
instance. The hybrid solver did not succeed in solving
all problems. Our key observations for the considered
formulation are:
The D-Wave hybrid solver does not, or just barely
profits from longer computation times.
The number of feasible solutions strongly affects the
probability for the D-Wave hybrid solver to find a
solution.
The outcome of the benchmark, i.e. a simple laptop
performing better than the D-Wave hybrid solver,
highlights a major challenge when applying quantum
computing to operational research problems: that it is
crucial to adapt the model formulation of the problem
to the capabilities and strengths of the solver. Indeed,
with Mixed-Integer Linear Programs being the
standard in classical optimization, the benchmark was
probably favorable for the classical system. A quantum
annealer natively solves QUBOs, i.e. quadratic
problems. Specifically, D-Wave’s hybrid system with
its constrained quadratic model is not a purely linear
solver. Unless we actually utilize quadratic terms, we
are benchmarking a more generic algorithm against a
more specialized one on the more specialized problem
class. Of course, that is likely to be advantageous for
the specialized algorithm. Therefore, in order to utilize
(hybrid) quantum systems to their full potential, we
might have to go beyond mixed-integer linear
formalizations. We are planning to do that in future
work.
Other potential research is pursuable on the subject
of end-of-horizon effects, especially due to so far
limited attention in the literature. Our approach could
for instance be modified or combined with other
methods like soft final inventory bounds. However,
while working on these extensions, we noticed that
there is a lack of an objective, or at least formal measure
of end-of-horizon effects. Therefore, benchmarking
different approaches to prevent those effects is
currently not possible. In order to choose our reward-
based method out of other possible approaches, we
played around with a measure based on the following
procedure resembling a rolling horizon framework. By
extending the horizon of an MIRP instance and
comparing the cumulative objective value of the
instance with end-of-horizon considerations and the
extension with the objective of the whole extended
instance, we could assess how well the instance
extends into the future and if infeasibilities occur. A
study of a similar notion or other proposals could be
addressed in future research.
ACKNOWLEDGMENTS
We thank the Hanseatic City of Hamburg for funding this
research within its programme ”Quantum computing for
shipping and maritime logistics in Hamburg” [2].
APPENDIX: BENCHMARK DETAILS
In order to generate different MIRP instances, we do
not randomly sweep all parameters. We change the
time frame |T|, the number of vessels |V |, the number
of ports |P|, as well as the parameters tt and prmin,
which set travel times and production/consumption
rates as follows:
( )
min max
min
1 1 1 1
1 1 2 2
1 2 1 2
1 2 2 1
min / min ,
nn
tt tt tt
tt
matrix of travel time TT
tt
tt
production consumtion rate PR pr PR n N
+++


+

=

+

+

=
As the size of the matrix indicates, we assume a
maximum of 4 ports, and restrict to lower indices if less
are employed. In our test instances the parameter prmin
to control minimum production/consumption rates is
chosen low enough to guarantee stability between the
supply and demand in product. In particular, we make
sure that the minimum supply does not surpass the
maximum demand, and both the minimum supply and
demand does not exceed the maximum turnover
volume. Note however that for some instances prmin has
been reduced significantly to compare the solvers
performance on small and large solution spaces.
The following are the exact parameters chosen to
stay fixed for all our test instances. Many of these
parameters are copied from the MIRP instance
LR1_1_DR1_3_VC1_V7a of the Benchmarking library
MIRPLIB [3]. With only one production port available,
we technically restrict to a one-to-many route
structure. In the cases where less than 4 ports or 7
vessels are used, only the values with lower indices are
considered.
121
Table 2. Benchmark results
name
vars
cons
CPLEX
D-Wave
5s
15s
30s
opt
5s
15s
30s
T15_V5_P4_tt5_prmin0
1453
1713
7204.8
7204.8
7204.8
7204.8
5709.3
5918.0
5953.8
T15_V5_P4_tt5vprmin20
1453
1713
7184.2
7184.2
7204.8
7204.8
5283.5
5483.6
3257.9
T15_V5_P4_tt5_prmin99
1453
1713
7204.8
7204.8
7204.8
7204.8
6414.4
5882.0
T16_V3_P4_tt1_prmin99
1241
1284
6365.8
6428.5
6433.7
6433.7
4412.9
5635.8
5451.5
T16_V3_P4_tt2_prmin0
1188
1277
6200.3
6433.7
6433.7
6433.7
4273.9
4609.1
4570.7
T16_V3_P4_tt2_prmin20
1188
1277
6205.6
6433.7
6433.7
6433.7
4473.8
4930.3
4610.7
T16_V3_P4_tt2_prmin99
1188
1277
6418.0
6433.7
6433.7
6433.7
3860.2
T16_V3_P4_tt3vprmin0
1135
1270
5899.3
5899.3
5899.3
5899.3
4077.3
4393.2
4187.4
T16_V3_P4_tt3_prmin20
1135
1270
5859.6
5864.2
5864.2
5864.2
4320.0
3879.5
4925.0
T16_V3_P4_tt3vprmin80
1135
1270
5859.6
5859.6
5859.6
5859.6
T16_V3_P4_tt4_prmin29
1082
1263
5269.5
5269.5
5269.5
5269.5
4014.8
4192.4
4538.4
T16_V4_P4_tt4vprmin0
1350
1546
6877.3
6877.3
6877.3
6877.3
4855.6
5004.2
4936.5
T16_V4_P4_tt4_prmin20
1350
1546
6779.9
6868.2
6868.2
6868.2
4881.5
5305.6
4538.8
T16_V4_P4_tt4_prmin99
1350
1546
6797.5
6797.5
6797.5
6797.5
T20_V3_P4_tt3_prmin0
1515
1594
6311.5
6508.9
6457.9
6966.1
4534.8
4451.8
3832.5
T20_V3_P4_tt3_prmin20
1515
1594
6457.9
6457.9
6696.0
6966.1
4417.1
4290.1
4469.9
T20_V3_P4_tt3_prmin85
1515
1594
6388.3
6942.1
6942.1
6942.1
T25_V2_P4_tt5_prmin15
1356
1530
4263.4
4693.3
4805.3
4805.3
2846.6
2967.7
3148.0
T30_V2_P4_tt5_prmin15
1686
1840
4753.6
4855.6
4831.4
4860.2
2790.1
T35_V2_P2_tt5_838
1192
4919.8
4919.8
4919.8
4919.8
4253.8
T35_V2_P3_tt5_prmin23
1372
1671
4938.2
4940.0
5530.1
5536.6
T35_V2_P4_tt5vprmin15
2016
2150
5225.8
5225.8
5372.2
5390.0
T45_V1_P2_prmin23
630
1033
3267.2
3267.2
3267.2
3267.2
T45_V1_P4_tt4_prmin9
1519
1924
3711.6
3789.4
3797.7
3813.6
2307.3
2345.0
T45_V3_P2_tt9_prmin39
1488
2013
6316.9
6534.5
6534.5
6534.4
T45_V4_P2_tt9_prmin39
1912
2500
7502.0
7559.0
7559.0
7559.0
T45_V6_P2_tt9_prmin39
2772
3478
9608.1
9608.1
9608.1
9608.1
T45_V7_P4_tt9_prmin0
7688
6873
8095.7
8665.2
13166.6
13939.2
3698.9
5710.2
2558.4
T45_V7_P4_tt9_prmin10
7688
6873
4903.8
9185.0
10382.3
13727.1
3182.1
1447.8
T45_V7_P4_tt9_prmin26
7688
6873
9488.6
11510.8
13668.0
Although most of the above parameters come with
an explicit time dependency, for our benchmarks we
have chosen them constant in time. Finally, we present
the instances and results in Table 2.
The laptop used for the classical benchmark was a
Dell Inc. Latitude 5411 with Intel(R) Core(TM) i7-
10850H CPU @ 2.70GHz, 16GB RAM and NVIDIA
GeForce MX250. On it we ran CPLEX 12.10.0.0 on
standard settings, except for the finite computation
time.
REFERENCES
[1] dwave-system documentation, 2021.
[2] Hamburg soll an die spitze im quantencomputing.
https://www.hamburg.de/bwi/medien/16613674/ 2022-
10-25-bwi-quantencomputing/, 2022. Accessed: 2023-11-
01.
[3] MIRPLIB a library of maritime inventory routing
problems. https://mirplib.scl.gatech.edu/ instances, 2023.
Accessed: 2023-11-01.
[4] Amira Abbas et al. Quantum Optimization: Potential,
Challenges, and the Path Forward. 12 2023. [5] Agostinho
Agra, Marielle Christiansen, Alexandrino Delgado, and
Luidi Simonetti. Hybrid heuristics for a short sea
inventory routing problem. European Journal of
Operational Research, 236(3):924935, 2014. Vehicle
Routing and Distribution Logistics.
[6] Agostinho Agra, Marielle Christiansen, Lars Magnus
Hvattum, and Filipe Rodrigues. Robust optimization for
a maritime inventory routing problem. Transportation
Science, 52(3):509525, 2018.
[7] Agostinho Agra, Marielle Christiansen, Kristine S Ivarsøy,
Ida Elise Solhaug, and Asgeir Tomasgard. Combined ship
routing and inventory management in the salmon
farming industry. Annals of operations research, 253:799
823, 2017.
[8] Agostinho Agra, Marielle Christiansen, and Laurence
Wolsey. Improved models for a single vehicle
continuous-time inventory routing problem with pickups
and deliveries. European Journal of Operational
Research, 297(1):164179, 2022.
[9] Henrik Andersson, Arild Hoff, Marielle Christiansen,
Geir Hasle, and Arne Løkketangen. Industrial aspects
and literature survey: Combined inventory management
and routing. Computers & Operations Research,
37(9):15151536, 2010.
[10] Mohamed Ben Ahmed, Onyemaechi Linda Okoronkwo,
Edwin Chimezie Okoronkwo, and Lars Magnus
Hvattum. Long-term effects of short planning horizons
for inventory routing problems. International
Transactions in Operational Research, 29(5):29953030,
2022.
[11] Hans J Briegel, David E Browne, Wolfgang Dur, Robert
Raussendorf, and Maarten Van den Nest. Measurement-
based quantum computation. Nature Physics, 5(1):1926,
2009.
122
[12] Jaeyoung Cho, Gino J Lim, Seon Jin Kim, and Taofeek
Biobaku. Liquefied natural gas inventory routing
problem under uncertain weather conditions.
International Journal of Production Economics, 204:18
29, 2018.
[13] Marielle Christiansen, Kjetil Fagerholt, Truls Flatberg,
Øyvind Haugen, Oddvar Kloster, and Erik H. Lund.
Maritime inventory routing with multiple products: A
case study from the cement industry. European Journal of
Operational Research, 208(1):8694, 2011.
[14] Leandro C. Coelho, Jean-Francois Cordeau, and Gilbert
Laporte. Thirty years of inventory routing.
Transportation Science, 48(1):119, 2014.
[15] IBM ILOG Cplex. V12. 1: User’s manual for cplex.
International Business Machines Corporation, 46(53):157,
2009.
[16] Stephane Dauzere-Peres, Atle Nordli, Asmund Olstad,
Kjetil Haugen, Ulrich Koester, Myrstad Per Olav, Geir
Teistklub, and Alf Reistad. Omya hustadmarmor
optimizes its supply chain for delivering calcium
carbonate slurry to european paper manufacturers.
Interfaces, 37(1):3951, 2007.
[17] Arijit De, Vamsee Krishna Reddy Mamanduru, Angappa
Gunasekaran, Nachiappan Subramanian, and Manoj
Kumar Tiwari. Composite particle algorithm for
sustainable integrated dynamic ship routing and
scheduling optimization. Computers & Industrial
Engineering, 96:201215, 2016.
[18] Line Eide, Gro Cesilie Hahjem Ardal, Nataliia Evsikova,
Lars Magnus Hvattum, and Sebastian Urrutia. Load-
dependent speed optimization in maritime inventory
routing. Computers & Operations Research, 123:105051,
2020.
[19] Kjetil Fagerholt, Lars Magnus Hvattum, Dimitri J.
Papageorgiou, and Sebastian Urrutia. Maritime
inventory routing: recent trends and future directions.
International Transactions in Operational Re- search,
30(6):30133056, 2023.
[20] Sebastian Feld, Christoph Roch, Thomas Gabor,
Christian Seidel, Florian Neukart, Isabella Galter,
Wolfgang Mauerer, and Claudia Linnhoff-Popien. A
hybrid solution method for the capacitated vehicle
routing problem using a quantum annealer. Frontiers in
ICT, 6:13, 2019.
[21] Claudio Gambella and Andrea Simonetto. Multiblock
admm heuristics for mixed-binary optimization on
classical and quantum computers. IEEE Transactions on
Quantum Engineering, 1:122, 2020.
[22] Stuart Harwood, Claudio Gambella, Dimitar Trenev,
Andrea Simonetto, David Bernal, and Donny Greenberg.
Formulating and solving routing problems on quantum
computers. IEEE Transactions on Quantum Engineering,
2:117, 2021.
[23] Philipp Hauke, Helmut G Katzgraber, Wolfgang
Lechner, Hidetoshi Nishimori, and William D Oliver.
Perspectives of quantum annealing: methods and
implementations. Reports on Progress in Physics,
83(5):054401, May 2020.
[24] Sovanmonynuth Heng, Dongmin Kim, Taekyung Kim,
and Youngsun Han. How to solve combinatorial
optimization problems using real quantum machines: A
recent survey. IEEE Access, 10:120106120121, 2022.
[25] Hirotaka Irie, Goragot Wongpaisarnsin, Masayoshi
Terabe, Akira Miki, and Shinichirou Taguchi. Quan- tum
annealing of vehicle routing problem with time, state and
capacity. In Quantum Technology and Optimization
Problems: First International Workshop, QTOP 2019,
Munich, Germany, March 18, 2019, Proceedings 1, pages
145156. Springer, 2019.
[26] Jana Ksciuk, Stefan Kuhlemann, Kevin Tierney, and
Achim Koberstein. Uncertainty in maritime ship routing
and scheduling: a literature review. European Journal of
Operational Research, 308(2):499524, 2023.
[27] Mingyu Li and Peter Schu¨tz. Planning annual lng
deliveries with transshipment. Energies, 13(6):1490, 2020.
[28] Catherine McGeoch and Pau Farr´e. The d-wave
advantage system: An overview. D-Wave Systems Inc.,
Burnaby, BC, Canada, Tech. Rep, 2020.
[29]Lluis-Miquel Munguia, Shabbir Ahmed, David A. Bader,
George L. Nemhauser, Yufen Shao, and Dimitri J.
Papageorgiou. Tailoring parallel alternating criteria
search for domain specific mips: Application to maritime
inventory routing. Computers & Operations Research,
111:2134, 2019.
[30] M.A. Nielsen and I.L. Chuang. Quantum Computation
and Quantum Information: 10th Anniversary Edition.
Cambridge University Press, 2010.
[31] Dimitri J Papageorgiou, Myun-Seok Cheon, Stuart
Harwood, Francisco Trespalacios, and George L
Nemhauser. Recent progress using matheuristics for
strategic maritime inventory routing. Modeling,
computing and data handling methodologies for
maritime transportation, pages 5994, 2018.
[32] Dimitri J. Papageorgiou, Myun-Seok Cheon, George
Nemhauser, and Joel Sokol. Approximate dynamic
programming for a class of long-horizon maritime
inventory routing problems. Transportation Science,
49(4):870885, 2015.
[33] Dimitri J. Papageorgiou, Ahmet B. Keha, George L.
Nemhauser, and Joel Sokol. Two-stage decom- position
algorithms for single product maritime inventory
routing. INFORMS Journal on Computing, 26(4):825847,
2014.
[34] Dimitri J. Papageorgiou, George L. Nemhauser, Joel
Sokol, Myun-Seok Cheon, and Ahmet B. Keha. Mirplib
a library of maritime inventory routing problem
instances: Survey, core model, and benchmark results.
European Journal of Operational Research, 235(2):350
366, 2014. Maritime Logistics.
[35] Niklas Pirnay, Vincent Ulitzsch, Frederik Wilde, Jens
Eisert, and Jean-Pierre Seifert. An in-principle super-
polynomial quantum advantage for approximating
combinatorial optimization problems. 12 2022.
[36] F Radan, SMT Fatemi Ghomi, SMJ Mirzapour Al-e
hashem, and Moeen Sammak Jalali. Maritime inventory
routing problem considering weather conditions and tide
at ports. Transportation Research Record, 2677(5):934
950, 2023.
[37] Filipe Rodrigues, Agostinho Agra, Marielle
Christiansen, Lars Magnus Hvattum, and Cristina
Requejo. Comparing techniques for modelling
uncertainty in a maritime inventory routing problem.
European Journal of Operational Research, 277(3):831
845, 2019.
[38] Jin-Hwa Song and Kevin C. Furman. A maritime
inventory routing problem: Practical approach.
Computers & Operations Research, 40(3):657665, 2013.
Transport Scheduling.
[39] Byron Tasseff, Tameem Albash, Zachary Morrell, Marc
Vuffray, Andrey Y. Lokhov, Sidhant Misra, and Carleton
Coffrin. On the emerging potential of quantum annealing
hardware for combinatorial optimization. 2022.
[40] United Nations Conference on Trade and Development
(UNCTAD). Review of maritime transport 2022.
[41] Amir Zojaji, Kiarash Soltaniani, Lars Magnus Hvattum,
and Sebastia´n Urrutia. Cyclic solutions to a maritime
inventory routing problem. Maritime Transport
Research, 3:100074, 2022.