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

, S. Rubbert

& A. Rizvanolli

Fraunhofer CML, Hamburg, Germany

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

were visualized by [38] as

graphs, with the nodes N

∗

and edges A

. 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

is visualized in orange and

divided into regular arcs and source and sink arcs. As

an example of the following- and preceding nodes sets,

(1,5) and PN

(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 ,

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:

( )

Revenue , , :

v v v v

n n n n a a

v V n N

C ta b l PP l BF b TC ta







=  −  − 





  

(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:

( ) ( )

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 , ,

p t p t p t p t

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:

( )

v v v

t t p

vi vi PT l t T v V

−



= +     



(4)

This constraint ensures that the vessel inventory vi

is correctly updated.

Berth condition:

( )

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:

pi PC n N  

(8)

The inventory pi of a port can never exceed its

capacity PC.

Vessel inventory bounds:

vi VC t T v V   

(9)

The inventory vi of a vessel can never exceed its

capacity VC.

Berth limit:

( )

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:

 

Set of consmption ports CP p P PT=  = −

While not formally a variable of the model, we will

frequently refer to the expression

( )

max

, , ,

p CP

q T sink

q T p T TT

vi PP

ta TC BF















−  − 













(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

((q,|T|),(p,|T|+TTq,p)):=

((q,|T|-TTq,p),(p,|T|).

Above we have denoted the upper and lower

bounds of the expected revenue er

( )

max

p CP

U PP







=−







(12)

( )

, , ,

: max

q T p T TT

p CP p CP

L BF TC







   







= − +









   





(13)

Finally, we can specify the objective function and

constraints of our general model.

Extended model objective function:

( )

, , , :

v v v v v

n n n n a a

v V n N

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

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

v v v

epr ta U v V







    







(15)

( )

, ,sink

v v v v

mer epr ta U v V







 −  −   







(16)

These constraints are equivalent to

( )

, ,sink

v v v

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:

v v v

mer erp U v V    

(17)

( )( )

v v v v v v

er mer er U L erp v V  + − −  

(18)

These constraints are equivalent to mer

:=max{e

,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 ,

PR pr PR n N=  

where

max

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

−

0.737

0.662

0.84

0.691

0.745

0.71

T16_V3_P4_tt2

1188

−

0.6

0.721

0.766

0.717

0.689

0.716

0.71

T16_V4_P4_tt4

1350

−

0.72

0.772

0.661

0.706

0.728

0.718

T15_V5_P4_tt5

1453

−

0.89

0.816

0.735

0.763

0.452

0.792

0.821

0.826

T20_V3_P4_tt3

1515

−

0.684

0.664

0.668

0.718

0.684

0.593

T45_V7_P4_tt9

7688

−

0.649

−

0.139

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 ,

tt tt tt

matrix of travel time 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

15s

30s

opt

15s

30s

T15_V5_P4_tt5_prmin0

1453

1713

7204.8

5709.3

5918.0

5953.8

T15_V5_P4_tt5vprmin20

1453

1713

7184.2

7204.8

5283.5

5483.6

3257.9

T15_V5_P4_tt5_prmin99

1453

1713

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

4273.9

4609.1

4570.7

T16_V3_P4_tt2_prmin20

1188

1277

6205.6

6433.7

4473.8

4930.3

4610.7

T16_V3_P4_tt2_prmin99

1188

1277

6418.0

6433.7

−

3860.2

T16_V3_P4_tt3vprmin0

1135

1270

5899.3

4077.3

4393.2

4187.4

T16_V3_P4_tt3_prmin20

1135

1270

5859.6

5864.2

4320.0

3879.5

4925.0

T16_V3_P4_tt3vprmin80

1135

1270

5859.6

−

T16_V3_P4_tt4_prmin29

1082

1263

5269.5

4014.8

4192.4

4538.4

T16_V4_P4_tt4vprmin0

1350

1546

6877.3

4855.6

5004.2

4936.5

T16_V4_P4_tt4_prmin20

1350

1546

6779.9

6868.2

4881.5

5305.6

4538.8

T16_V4_P4_tt4_prmin99

1350

1546

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

6696.0

6966.1

4417.1

4290.1

4469.9

T20_V3_P4_tt3_prmin85

1515

1594

6388.3

6942.1

−

T25_V2_P4_tt5_prmin15

1356

1530

4263.4

4693.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

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

5372.2

5390.0

−

T45_V1_P2_prmin23

630

1033

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.4

−

T45_V4_P2_tt9_prmin39

1912

2500

7502.0

7559.0

−

T45_V6_P2_tt9_prmin39

2772

3478

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):924–935, 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):509–525, 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):164–179, 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):1515–1536, 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):2995–3030,

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):19–26,

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):86–94, 2011.

[14] Leandro C. Coelho, Jean-Francois Cordeau, and Gilbert

Laporte. Thirty years of inventory routing.

Transportation Science, 48(1):1–19, 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):39–51, 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:201–215, 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):3013–3056, 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:1–22, 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:1–17, 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:120106–120121, 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

145–156. 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):499–524, 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:21–34, 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 59–94, 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):870–885, 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):825–847,

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):657–665, 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.