Drent et al.

Dynamic Supply Mode Selection

Efficient Emission Reduction Through Dynamic Supply Mode Selection

Melvin Drent School of Industrial Engineering, Eindhoven University of
Technology, Eindhoven, the Netherlands, PO BOX 513, 5600MB<EMAIL_ADDRESS>Poulad Moradi, Joachim Arts Luxembourg Centre for Logistics and Supply Chain
Management, University of Luxembourg, Luxembourg City, Luxembourg, 6, rue
Richard Coudenhove-Kalergi L-1359, {poulad.moradi<EMAIL_ADDRESS>
Reducing the carbon footprint of global supply chains is a challenge for many
companies. Governmental emission regulations are increasingly stringent, and
consumers are increasingly environmentally conscious. Companies should
therefore integrate carbon emissions in their supply chain decision making. In
this paper, we study the inbound supply mode and inventory management decision
making for a company that sells an assortment of products. Stochastic demand
for each product arrives periodically and unmet demand is backlogged. Each
product has two distinct supply modes that differ in terms of their carbon
emissions, speed, and costs. The company needs to decide when to ship how much
using which supply mode such that total holding, backlog, and procurement
costs are minimized while the emissions associated with different supply modes
across the assortment remains below a target level. We assume that shipment
decisions for each product are governed by a dual-index policy for which we
optimize the parameters. We formulate this decision problem as a mixed integer
linear program that we solve through Dantzig-wolfe decomposition. We benchmark
our decision model against two state-of-the-art approaches in a large test-bed
based on real-life carbon emissions data. Relative to our decision model, the
first benchmark lacks the flexibility to dynamically ship products with two
supply modes while the second benchmark makes supply mode decisions for each
product individually. Our computational experiment shows that our decision
model can outperform the first and second benchmark by up to 15 and 40
percent, respectively, for moderate carbon emission reduction targets.

Inventory; sustainability; policy comparison; multi-item; carbon emissions

##  1 Introduction

The transportation sector has consistently been one of the most polluting
European sectors for more than a decade now, and it is projected to remain so
for the foreseeable future (European Environment Agency 2020). This,
unfortunately, appears to be a trend that stretches beyond Europe. Recent
analysis indicates that the G20 countries, currently responsible for 80% of
the global greenhouse gas (GHG) emissions, will see an increase of 60% in
their transportation sector emissions by 2050 (Vieweg et al. 2018). Prominent
global climate targets, such as the ones outlined in the Paris agreement, will
soon become unattainable (European Environment Agency 2020, United Nations
Environment Programme 2020).

In light of the above, the European Union (EU) recently announced the Green
Deal, a framework containing climate targets and policy initiatives that sets
the EU on a path to reach carbon neutrality by 2050. The Green Deal is legally
enshrined in the European Climate Law, which states that member states are
legally committed to meet the targets, and face penalties in case they do not
meet these targets. Being among the most polluting sectors, a key part of the
Green Deal relates to policy initiatives that impact the transportation
sector. For instance, the EU plans to extend the European emissions trading
scheme (EU ETS) to include both road and maritime transport (Abnett 2020,
European Commission 2020). Under the ETS, which until now includes only air
transport, the EU enforces a cap on the total amount of GHG emissions from
sectors covered by the scheme. The EU also investigates whether to increase
fossil fuel taxation, thereby effectively raising the price of GHG emissions.
More and more companies are reducing their emissions voluntarily as part of
their corporate social responsibility. If not penalized by governments,
companies that excessively pollute might still lose revenues as
environmentally conscious customers take their business elsewhere (Dong et al.
2019).

The developments described above highlight the urgency for companies to
explicitly incorporate GHG emissions in their supply chain decision making. In
this paper, we study the inbound supply mode and inventory decision making of
a company that sells an assortment of products which are sourced from outside
suppliers. The company wishes to keep the total GHG emissions associated with
using different supply modes across the assortment below a certain target
level in the most economically viable manner. As is often the case in
practice, the company may rely on a third party logistics (3PL) provider for
the inbound transport of the products. 3PLs typically offer several transport
modes for the transportation of products – these may differ in terms of
transportation costs, transit times, and GHG emissions. Alternatively, there
may be different suppliers (e.g. a near and offshore supplier) for a product
that naturally have different costs and transportation emissions. In the
remainder of the paper, we will use the terminology of a 3PL provider that
offers multiple transportation modes.

The company can utilize the heterogeneity in the fleet of the 3PL to its
advantage. While some transport modes are low emitting but slow, others may be
fast but result in more emissions. Fast transport modes also typically come at
the expense of a cost premium, and yet they are often relied upon when
responsiveness is required (e.g. in case of imminent stock outs). Thus the
company should rely dynamically on both transport modes. Implementing this
holistically across the entire assortment of products allows the company to
reduce emissions significantly for products for which it is relatively cost-
efficient to do so and less for products for which this is more expensive. It
additionally enables the company to reduce the total inventory and
transportation costs by shipping the majority of products with a relatively
cheap but slow transport modes while simultaneously resorting to faster but
more expensive and often more polluting transport modes whenever expedited
shipments are needed. While the advantages of dynamically selecting different
transport modes are evident, two important and interrelated questions remain:

  1. 1.

When should the company ship how many units of which product with what
transport mode?

  2. 2.

What is the value of dynamically shipping products with different transport
modes?

These questions are interesting but also intricate when one wishes to answer
them for an entire assortment of products where the combined total of GHG
emissions from transportation must not exceed a certain target level.

To tractably answer the questions above, we focus on the setting where the 3PL
offers two distinct transport modes for the transport of each product (or,
equivalently, the setting where the company has already decided on the two
transport modes for each product). These transport modes need not be the same
for every product; they will depend on the characteristics of the suppliers as
well as the 3PL (e.g., some products can be transported using aircraft or
rolling stock while others can be shipped via inland waterways or ocean
shipping). We consider long distance and/or high volume transport lanes where
the impact of the transport mode decisions of any particular individual
shipper on the actual shipping and carbon footprint is negligible. The company
decides periodically how many units it wishes to transport with what transport
mode and incurs mode specific unit transportation costs. Shipments arrive at
the company after a deterministic transit time that depends on the transport
mode that is used. Demand for each product in every period is stochastic and
independent and identically distributed across periods. Any demand in excess
of on-hand inventory is backlogged and satisfied in later periods. The company
incurs per unit holding and backorder penalty costs, and the specific cost
parameters may vary from product to product. The company seeks to minimize the
long-run average holding, backorder, and transportation costs while keeping
the total long-run average GHG emissions from transportation of the entire
assortment below a certain target level.

It is well-known that the optimal policy for the inventory system described
above is complex, even in the simplest case of a single product and absent of
the emission constraint (Whittemore and Saunders 1977, Feng et al. 2006). For
the control of each product, we therefore use a heuristic policy that is
originally due to Veeraraghavan and Scheller-Wolf (2008). They show that their
so-called dual-index policy performs quite well compared to the optimal
policy. The dual-index policy tracks two inventory positions for each product:
The slow inventory position, which equals the on-hand inventory plus all in-
transit products minus backlog, and the fast inventory position, which is
defined similarly but includes only those in-transit products that are due to
arrive within the transit time of the fastest transport mode. Under the dual-
index policy, we place orders with both modes such that these inventory
positions are kept at (or above) certain target levels, also referred to as
base-stock levels. As such, the dual-index policy dynamically prescribes
shipment quantities for both transport modes based on the net inventory level
and the number of products that are still in-transit. To find the optimal
base-stock levels for the entire assortment of products, we formulate the
decision problem as a non-linear non-convex integer programming problem. A
partition reformulation of this problem allows us to use column generation
techniques to solve the decision problem. These techniques enable us to
decompose the complex multi-product decision problem into simpler sub-problems
per product. Leveraging a separability result of Veeraraghavan and Scheller-
Wolf (2008), we show that this sub-problem constitutes a special Newsvendor
problem that can be solved efficiently through a simulation-based optimization
procedure.

The main contributions of this paper are:

  1. 1.

We are the first to study dynamic mode selection for an assortment of products
with stochastic demand where the total average GHG emissions from the inbound
transport of those products must be kept below a certain target level.

  2. 2.

We provide a tractable optimization model that finds a tight lower bound on
the optimal solution as well as near-optimal feasible solutions within
reasonable time. We show that our mathematical formulation of the decision
problem allows us to decompose the non-linear non-convex integer programming
problem into sub-problems per product. We leverage results from Veeraraghavan
and Scheller-Wolf (2008) to show that the sub-problems can be solved
efficiently through a one dimensional search procedure in which each instance
constitutes a Newsvendor type problem that is readily solved through
simulation.

  3. 3.

We perform an extensive computational experiment based on data from different
industries. Through these experiments:

    1. i.

We establish the value of dynamic mode selection by comparing our model with a
model in which only one transport mode per assortment product can be used.
This value can go up to 15 percent in cost savings;

    2. ii.

We show that decomposing an aggregate carbon emission reduction target into
targets for each product in the assortment individually is financially
detrimental. Our holistic approach can lead to cost savings of over 40 percent
relative to the approach with reduction targets per individual product;

    3. iii.

We find that the emission differences between transport modes relative to the
cost difference between modes is the main determinant of emission reduction
potential for a given assortment. In our experiments we find that, 20 percent
of the products for which this ratio is highest contribute between 59 and 94
percent of the emission reduction.

The remainder of this paper is organized as follows. In Section 2, we review
the existing literature and position our work within the literature. Section 3
contains the model description as well as the mathematical formulation of the
decision problem. A column generation procedure to solve the decision problem
is provided in Section 4. We subsequently report on an extensive computational
experiment in Section 5, and we provide concluding remarks in Section 6.

##  2 Related literature

This paper integrates carbon emissions from inbound transportation into an
inventory system with two supply modes. As such, our work contributes to the
large stream of literature that studies multi-mode or multi-supplier inventory
systems. For an excellent overview of such systems, we refer the reader to the
review papers of Thomas and Tyworth (2006), Engebrethsen and Dauzére-Pérés
(2019) and Svoboda et al. (2021), and references therein. We also contribute
to the extensive body of literature that revolves around the integration of
environmental aspects into supply chain decision making; see Dekker et al.
(2012), Brandenburg et al. (2014), Barbosa-Póvoa et al. (2018), and references
therein, for an overview of this field. In what follows, we focus on
contributions that are most relevant to the present paper.

The decision how many products to order from which supplier is considered a
canonical problem in the inventory management literature. It has been studied
extensively since the sixties, mostly under the assumption that lead times are
deterministic, that unmet demand is backlogged, and that only two distinct
suppliers are at the disposal of the decision maker; the fastest being more
expensive than the slowest. Fukuda (1964) and Whittemore and Saunders (1977)
were the first to study this system. Assuming periodic review, they show that
its optimal policy is a simple base-stock rule only under the assumption that
the difference between the lead times of both suppliers is one period. For
general lead time differences, the optimal policy is complex and can only be
computed through dynamic programming for small instances. Since then, most
researchers have focused on developing well-performing heuristic polices for
which the best control policy parameters can be tractably obtained.

In this paper, we rely on the so-called dual-index policy to decide upon the
shipment sizes for both transport modes for each product. Under this policy,
which is originally due to Veeraraghavan and Scheller-Wolf (2008), two
different inventory positions are kept track off: One that includes all
outstanding shipments and one that includes only those outstanding shipments
that are due to arrive within the lead time time of the fastest mode.
Veeraraghavan and Scheller-Wolf (2008) show numerically that the dual-index
policy performs well compared to the optimal policy. In fact, Drent and Arts
(2022) show that the dual-index policy is asymptotically optimal as the cost
of the fastest transport mode and the backorder penalty cost become large
simultaneously. The policy has received quite some attention in recent years
(see, e.g., Sheopuri et al. 2010, Arts et al. 2011, Sun and Van Mieghem 2019).
We employ the dual-index policy because it is intuitive, has good performance,
and can be optimized efficiently. Unlike the present work, the dual-index
policy has so far been studied exclusively in single product settings under
conventional cost criteria absent of any emission considerations.

Within the transportation literature, inventory systems with multiple
transport modes have received considerable attention too. To properly embed
the present work in the existing literature, we group contributions to this
field into two categories depending on the modelling assumptions regarding the
usage of the available transport modes (c.f. Engebrethsen and Dauzére-Pérés
2019). The first category, which we refer to as dynamic mode selection, is
concerned with inventory systems in which multiple transport modes are used
simultaneously over a given (possibly infinite) planning horizon. Since we
study an infinite horizon periodic review inventory model in which products
can be transported with two distinct modes in each period, our work falls into
this category – as do all the inventory papers with two suppliers described so
far. Only few papers exist in this category that explicitly account for carbon
emissions, and the few that do differ substantially from the present work in
terms of modelling choices as well as analysis. They either assume
deterministic demand and a finite horizon (Palak et al. 2014) or study the
closely related yet different problem of splitting an order among several
transport modes (Konur et al. 2017). While not explicitly modeling carbon
emissions, Dong et al. (2018) and Lemmens et al. (2019) also study the benefit
of dynamically switching between multiple transport modes in the context of
multi-modal transport. They show that this can lead to more usage of less
polluting transport modes without compromising on costs or responsiveness.
Different from our work, all papers mentioned above consider the inventory
control and transport mode decisions for a single product only.

The second category concerns inventory systems in which a single transport
mode is selected a priori at the start of a planning horizon; all
replenishment orders until the end of that planning horizon are then shipped
with this mode. We refer to this category as static mode selection. Two papers
belonging to this category are particularly relevant to our work. Hoen et al.
(2014b) study a periodic review inventory system under backlogging where
inbound transport is outsourced to a 3PL that offers multiple transport modes.
Assuming base-stock control for each mode, they are interested in selecting
the transportation mode that leads to the lowest long-run average total cost
consisting of holding, backlogging, ordering, and emission costs. For
calculating transportation emissions, they rely on the well-known NTM
methodology (we discuss this methodology in more detail later in Section 3.1
and Appendix 7). We extend Hoen et al. (2014b) in two important directions.
First, we move from static to dynamic mode selection, thereby incorporating
the flexibility to dynamically switch between different transport modes for
each product. Second, we consider an assortment of products under a single
constraint on the total average transportation emissions from those products.
Hoen et al. (2014a) consider a similar constraint in a multi-product variant
of the setting of Hoen et al. (2014b) under the assumption that demand is
deterministic and inversely related to the price set by the decision maker.
They show that because of the portfolio effect of such an assortment-wide
emission constraint, carbon emissions from transportation can be reduced
substantially at hardly any additional cost.

The dual-index policy studied in this paper has the appealing feature that it
can mimic static mode selection. This is useful in our computational
experiment where we establish the added value of dynamic mode selection over
static mode selection. A closely related paper in that respect is Berling and
Martínez de Albéniz (2016) who study dynamic speed optimization of a single
transport mode in a single-product stochastic inventory problem. They show
that the value of dynamically controlling the speed of outstanding shipments,
as opposed to a static speed policy, can be significant, both financially and
from a carbon emission perspective.

Our review so far has almost exclusively revolved around papers on multi-
period inventory systems. We note that there is also a stream of literature
that integrates carbon emissions into single period multi-supplier models, see
e.g., Rosic̆ and Jammernegg (2013), Arıkan and Jammernegg (2014), and Chen and
Wang (2016). Similar to the majority of the papers discussed so far, these
papers focus on single-product settings.

##  3 Model description

In this section, we first provide a description of the inventory system under
consideration and introduce the notation that we use throughout this paper. We
then describe the policy we propose to dynamically ship products with two
transport modes. We conclude with providing a mathematical formulation of the
decision problem.

###  3.1 Description and notation

We consider a company that sells an assortment of products. The inventories
for these products are replenished from external suppliers through a third
party logistics provider (3PL). A 3PL often offers several transport modes. We
focus on the setting where the company has already decided upon two distinct
transport modes that it would like to use for the transport of each product.
These two transport modes will differ in terms of costs, lead times,
emissions, or a combination thereof. Given these two transport modes for each
product, the operational question that remains is how many units of each
product the company should transport using which transport mode at what time
so that costs –holding, backlog, and ordering– are minimized and an overall
emission constraint is met. Companies will increasingly impose such
constraints, either voluntarily or due to government regulation.

The inventory system runs in discrete time with $t\in\mathbb{N}_{0}$ denoting
the period index. Without loss of generality, we assume that the period is of
unit length and coincides with the review epoch. Let
$\mathnormal{J}=\\{1,2,\ldots,|J|\\}$ denote the nonempty set of products that
the company offers for sale. Demand for product $j\in\mathnormal{J}$ across
periods is a sequence of non-negative independent and identically distributed
(i.i.d.) random variables $\\{D_{j}^{t}\\}$. Any demand in excess of on-hand
inventory is backlogged. Let $I^{t}_{j}$ denote the net inventory level (on-
hand inventory minus backlog) of product $j$ at the beginning of period $t$
after any outstanding orders have arrived. Each unit of product $j$ in on-hand
inventory $(I^{t}_{j}-D^{t}_{j})^{+}$ carried over to the next period and
incurs a holding cost $h_{j}>0$. Similarly, each unit of product $j$ in
backlog $(D^{t}_{j}-I^{t}_{j})^{+}$ incurs a penalty cost $p_{j}>0$. Here we
use the standard notation $x^{+}=\max(0,x)$.

Each product can be shipped using two distinct transport modes from one
supplier (or, equivalently, using one or two distinct transport modes from two
distinct suppliers). Let $\mathnormal{M}=\\{f,s\\}$ denote the set of
available transport modes, where we use $f$ and $s$ to refer to the faster and
slower transport mode, respectively. Associated with the transport of one unit
of product $j\in\mathnormal{J}$ with mode $m\in\mathnormal{M}$ is a cost
$c_{j,m}\geq 0$, a deterministic lead time $l_{j,m}\in\mathbb{N}_{0}$, a
distance traveled from the supplier to the company $d_{j,m}>0$, and a certain
number of units CO2 emission $e_{j,m}\geq 0$. The weight of one unit of
product $j$ is denoted $w_{j}>0$. Recall that the company outsources its
transport to a 3PL provider and hence has no control over the actual shipping.
We therefore consider variable emissions that depend only on product and
transport mode specific characteristics as well as on distance traveled, and
we refrain from incorporating a fixed emission factor per actual shipment.
This is a reasonable assumption for long distance and/or high volume transport
lanes where the impact of the decisions of any particular individual shipper
on the carbon footprint are negligible. In line with previous literature that
models transportation emissions in the context of mode selection (e.g., Hoen
et al. 2014a, b), we endow $e_{j,m}$ with the following structure which is
based on the NTM methodology:

$e_{j,m}=w_{j}(a_{m}+d_{j,m}b_{m}),$ (1)

where $a_{m}\geq 0$ and $b_{m}>0$ are a fixed and variable transport mode
specific emission constant, respectively. Consistent with the NTM methodology,
we assume that each product is shipped with an averagely loaded transport
mode. We define the lead-time difference between the fast and slow mode as
$l_{j}=l_{j,s}-l_{j,f}\geq 0$ for each product $j\in\mathnormal{J}$.
Conventional literature on dual model problems (e.g., Sheopuri et al. 2010)
imposes the assumption that the cost premium of using the fast mode does not
exceed the lead-time difference multiplied by the penalty costs, i.e.,
$(c_{j,f}-c_{j,s})<l_{j}p_{j}$ to ensure that using the fast mode is
attractive. We do not impose this assumption as whether using the fast supply
mode will also depend on the target carbon reduction. When the fast supply
mode is less polluting than the fast mode, the fast mode may become attractive
even when $(c_{j,f}-c_{j,s})<l_{j}p_{j}$. Conversely, the fast supply mode may
become unattractive when the fast supply mode is more polluting, even when
$(c_{j,f}-c_{j,s})<l_{j}p_{j}$. Thus our model allows for situations where,
e.g., the expensive transport mode is either the fastest and most polluting or
the fastest and least polluting. Finally, the amount of items of product
$j\in\mathnormal{J}$ to be shipped with transport mode $m\in\mathnormal{M}$ in
period $t$ is denoted by $Q^{t}_{j,m}$. With this notation, observe that
shipments $Q^{t-{l_{j,f}}}_{j,f}$ and $Q^{t-{l_{j,s}}}_{j,s}$ arrive in period
$t$ so that we can write the following recursion for the inventory level
$I_{j}^{t}$ of each product $j$:

$I_{j}^{t}=I_{j}^{t-1}-D_{j}^{t-1}+Q^{t-{l_{j,f}}}_{j,f}+Q^{t-{l_{j,s}}}_{j,s}.$

All notation introduced so far as well as notation that we will introduce
later is summarized in Table 1.

Table 1: Overview of notation. Notation |  Description
---|---
Sets |
$\mathnormal{J}$ |  Assortment; Set of all products.
$\mathnormal{M}$ |  Set of available transport modes, i.e. $\mathnormal{M}=\\{f,s\\}$.
Input |
$D^{t}_{j}$ |  Random demand for product $j\in\mathnormal{J}$ in period $t\in\mathbb{N}_{0}$.
$p_{j}$ |  Penalty cost for one unit of product $j\in\mathnormal{J}$ in backlog carried over to the next period.
$h_{j}$ |  Holding cost for one unit of product $j\in\mathnormal{J}$ in on-hand inventory carried over to the next period.
$c_{j,m}$ |  Cost of shipping one unit of product $j\in\mathnormal{J}$ with transport mode $m\in\mathnormal{M}$.
$l_{j,m}$ |  Transportation lead time for product $j\in\mathnormal{J}$ by transport mode $m\in\mathnormal{M}$.
$l_{j}$ |  Transportation lead time difference between the fast and slow mode for product $j\in\mathnormal{J}$, i.e. $l_{j}=l_{j,s}-l_{j,f}$.
$d_{j,m}$ |  Distance for the transport of product $j\in\mathnormal{J}$ with transport mode $m\in\mathnormal{M}$.
$a_{m}$ |  Fixed emission constant corresponding with transport mode $m\in\mathnormal{M}$.
$b_{m}$ |  Variable emission constant corresponding with transport mode $m\in\mathnormal{M}$.
$w_{j}$ |  Unit weight of product $j\in\mathnormal{J}.$
$e_{j,m}$ |  Total units CO2 emission associated with shipping one unit of product $j\in\mathnormal{J}$ with transport mode $m\in\mathnormal{M}$.
$\mathcal{E}^{max}$ |  The maximally allowable carbon emissions for the transport of the entire assortment of products.
Decision variables |
$S_{j,m}$ |  Base-stock level for product $j\in\mathnormal{J}$ and transport mode $m\in\mathnormal{M}$.
$\Delta_{j}$ |  Difference between the slow and fast base-stock level for product $j\in\mathnormal{J}$, i.e. $S_{j,s}-S_{j,f}$.
$\mathbf{S}_{f}$ |  The vector $(S_{1,f},S_{2,f},\ldots,S_{|\mathnormal{J}|,f})$.
$\mathbf{\Delta}$ |  The vector $(\Delta_{1},\Delta_{2},\ldots,\Delta_{|\mathnormal{J}|})$.
State variables |
$I^{t}_{j}$ |  Inventory level of product $j\in\mathnormal{J}$ at the beginning of period $t\in\mathbb{N}_{0}$ after orders have arrived.
$IP^{t}_{j,f}$ |  Fast inventory position of product $j\in\mathnormal{J}$ in period $t\in\mathbb{N}_{0}$ before shipping orders.
$IP^{t}_{j,s}$ |  Slow inventory position of product $j\in\mathnormal{J}$ in period $t\in\mathbb{N}_{0}$ after shipping orders with the fast transport mode.
$Q^{t}_{j,m}$ |  Amount of product $j\in\mathnormal{J}$ shipped with transport mode $m\in\mathnormal{M}$ in period $t\in\mathbb{N}_{0}$.
$O^{t}_{j}$ |  The overshoot of product $j\in\mathnormal{J}$ in period $t\in\mathbb{N}_{0}$, i.e. $(IP^{t}_{j,f}-S_{j,f})^{+}.$
Output of model |
$C(\mathbf{S}_{f},\mathbf{\Delta})$ |  Total long-run average holding, backlog, and ordering costs under a given control policy $(\mathbf{S}_{f},\mathbf{\Delta})$.
$\mathnormal{E}(\mathbf{S}_{f},\mathbf{\Delta})$ |  Total emissions from transportation under a given control policy $(\mathbf{S}_{f},\mathbf{\Delta})$.
$C^{UB}_{P}(C^{LB}_{P})$ |  Upper (lower) bound for the optimal solution to Problem $(P)$.

###  3.2 Control policy

It is well-known that even for the simplest case where $|J|=1$ and absent of
the emission constraint, the policy that prescribes the optimal shipment
quantities is complex and can only be computed through dynamic programming for
very small instances that are arguably not representative for practice
(Whittemore and Saunders 1977, Feng et al. 2006). For the control of this
inventory system, we therefore use a heuristic policy that is originally due
to Veeraraghavan and Scheller-Wolf (2008). They show numerically that their
so-called dual-index policy performs quite well compared to the optimal
policy. The dual index policy tracks two indices: One that contains all orders
that are still in-transit and one that contains only those in-transit orders
that are due to arrive within the lead time of the fast mode. Based on these
outstanding orders, the policy dynamically ships orders with both modes to
keep these indices at certain target levels. In line with standard inventory
management nomenclature, we also refer to these target levels as base-stock
levels. More specifically, the policy operates as follows. At the beginning of
every period $t$ after orders $Q^{t-{l_{j,f}}}_{j,f}$ and
$Q^{t-{l_{j,s}}}_{j,s}$ have arrived, we review the fast inventory position,
which includes all in-transit orders – i.e. shipped with both the slow and the
fast transport mode – that will arrive within the lead time of the fast
transport mode:

$IP^{t}_{j,f}=I^{t}_{j}+\sum_{k=t-l_{j,f}+1}^{t-1}Q^{k}_{j,f}+\sum_{k=t-l_{j,s}+1}^{t-l_{j}}Q^{k}_{j,s}.$

Then, if necessary, we place order $Q^{t}_{j,f}$ with the fast transport mode
to raise the fast inventory position to its target level $S_{j,f}$. That is,
the amount of product $j$ shipped in period $t$ with the fast transport mode
equals:

$Q^{t}_{j,f}=(S_{j,f}-IP^{t}_{j,f})^{+}.$

After placing the fast shipment order, we inspect the slow inventory position,
which includes the fast order just placed

$IP^{t}_{j,s}=I^{t}_{j}+\sum\limits_{k=t-l_{j,f}+1}^{t}Q^{k}_{j,f}+\sum\limits_{k=t-l_{j,s}+1}^{t-1}Q^{k}_{j,s},$

and ship an order with the slow transport mode such that this inventory
position is raised to its target level $S_{j,s}$, with $S_{j,s}\geq S_{j,f}$
since the fast inventory position is contained in the slow inventory position.
Thus the amount of product $j$ shipped with the slow transport mode in period
$t$ equals:

$Q^{t}_{j,s}=S_{j,s}-IP^{t}_{j,s}.$

Note that, contrary to the fast inventory position, the slow inventory
position can never be larger than its base stock level $S_{j,s}$. After
shipping both orders, demand $D_{j}^{t}$ is satisfied or backlogged, depending
on whether there is sufficient inventory available or not. The period then
concludes with charging holding or backlog costs.

The order of events in a period $t$ for each product $j$ is thus as follows:

  1. 1.

Orders $Q^{t-{l_{j,f}}}_{j,f}$ and $Q^{t-{l_{j,s}}}_{j,s}$ arrive with the
fast and slow transport mode, respectively, and are added to the on-hand
inventory $I^{t}_{j}$.

  2. 2.

Review the fast inventory position and ship order $Q^{t}_{j,f}$ with the fast
transport mode at unit cost $c_{j,f}$.

  3. 3.

Review the slow inventory position and ship order $Q^{t}_{j,s}$ with the slow
transport mode at unit cost $c_{j,s}$.

  4. 4.

Demand $D^{t}_{j}$ occurs and is satisfied from on-hand inventory if possible,
and otherwise backlogged.

  5. 5.

Incur a cost $h_{j}$ for any unit in on-hand inventory
$(I_{j}^{t}-D^{t}_{j})^{+}$ and a cost $p_{j}$ for any unit in backlog
$(D^{t}_{j}-I^{t}_{j})^{+}$.

Observe that under a dual-index policy, slow orders entering the information
horizon of the fast transport mode may cause the fast inventory position to
exceed its target level. The amount by which the fast inventory position
exceeds its target level is referred to as the overshoot. The fast inventory
position in period $t$ after placing orders with both modes thus equals
$S_{j,f}+O^{t}_{j}$, where $O^{t}_{j}$ denotes the overshoot for product $j$
in period $t$:

$O^{t}_{j}=IP^{t}_{j,f}+Q^{t}_{j,f}-S_{j,f}=(IP^{t}_{j,f}-S_{j,f})^{+}.$

Later, in Section 4.2, we shall see that computing the steady state
distribution of the overshoot is crucial for determining the performance of a
given control policy for a single product.

We furthermore define $\Delta_{j}=S_{j,s}-S_{j,f}$, $j\in\mathnormal{J}$, so
that the control policy for a product can be specified in terms of its base-
stock levels $S_{j,s}$ and $S_{j,f}$ or in terms of its base-stock level for
the fast transport mode $S_{j,f}$ and the difference $\Delta_{j}$. We mostly
use the latter specification in our subsequent analysis. A control policy
$(\mathbf{S}_{f},\Delta)$ for the entire assortment of products consists of
the vectors $\mathbf{S}_{f}=(S_{1,f},S_{2,f},\ldots,S_{|J|,f})$ and
$\mathbf{\Delta}=(\Delta_{1},\Delta_{2},\ldots,\Delta_{|J|})$.

In what follows, for all sequences of random variables $X^{t}$, we define
their stationary expectation as
$\mathbb{E}[X]=\lim_{T\rightarrow\infty}(1/T)\sum\nolimits_{t=0}^{T}X^{t}$ and
their distribution as $\mathbb{P}(X\leq
x)=\lim_{T\rightarrow\infty}(1/T)\sum\nolimits_{t=0}^{T}\mathbbm{1}\\{{X^{t}\leq
x}\\}$, where $\mathbbm{1}\\{A\\}$ is the indicator function for the event
$A$. Whenever we drop the period index $t$ we refer to the generic stationary
random variable $X$ with expectation and distribution as defined above.

###  3.3 Decision problem

For a given control policy $(\mathbf{S}_{f},\mathbf{\Delta})$, we define the
total long-run average holding, backlog, and ordering costs per period for the
entire assortment of products as

$\displaystyle C(\mathbf{S}_{f},\mathbf{\Delta})$
$\displaystyle=\sum\limits_{j\in\mathnormal{J}}C_{j}(S_{j,f},\Delta_{j})$
$\displaystyle=\sum\limits_{j\in\mathnormal{J}}\left(h_{j}\mathbb{E}[(I_{j}-D_{j})^{+}]+p_{j}\mathbb{E}[(D_{j}-I_{j})^{+}]+\sum\limits_{m\in\mathnormal{M}}c_{j,m}\mathbb{E}[Q_{j,m}]\right),$
(2)

and the total emissions as

$\mathnormal{E}(\mathbf{S}_{f},\mathbf{\Delta})=\sum\limits_{j\in\mathnormal{J}}\mathnormal{E}_{j}(S_{j,f},\Delta_{j})=\sum\limits_{j\in\mathnormal{J}}\sum\limits_{m\in\mathnormal{M}}e_{j,m}\mathbb{E}[Q_{j,m}],$

where it is understood that the expectation operators are conditional on the
control policy $(\mathbf{S}_{f},\mathbf{\Delta})$. Veeraraghavan and Scheller-
Wolf (2008) show that $C(\mathbf{S}_{f},\mathbf{\Delta})$ is well-defined for
any control policy $(\mathbf{S}_{f},\mathbf{\Delta})$ as long as
$\mathbb{E}[D_{j}]<\infty$ for all products $j\in\mathnormal{J}$.

The objective of our decision problem is to minimize the total long-run
average costs while keeping the total emissions below a target level
$\mathcal{E}^{max}$. Combining the above-mentioned leads to the following
mathematical formulation of our decision problem which we refer to as problem
$(P)$:

$\displaystyle(P)\qquad\qquad$ $\displaystyle\displaystyle\min$ $\displaystyle
C(\mathbf{S}_{f},\mathbf{\Delta})$ subject to
$\displaystyle\mathit{E}(\mathbf{S}_{f},\mathbf{\Delta})\leq\mathcal{E}^{max},$
$\displaystyle\mathbf{S}_{f}\in\mathbb{R}^{|J|},\quad\mathbf{\Delta}\in\mathbb{R}_{0}^{|J|}.$

Let ($\mathbf{S}_{f}^{*},\mathbf{\Delta}^{*}$) denote an optimal solution to
problem $(P)$ and let $C_{P}$ be the corresponding optimal cost. Note that
Problem $(P)$ is a non-linear non-convex knapsack problem where more than one
copy of each item can be selected. It is well-known that even the simplest
types of such knapsack problems are $\mathcal{NP}$-hard (e.g. Kellerer et al.
2004). Since our knapsack is more complex, we conclude that Problem $(P)$ also
falls in that same complexity class; it is hence likely that also for our
problem no exact polynomial time solution algorithm exists.

We remark that Problem $(P)$ enables companies to reduce carbon emissions from
their inbound logistics by imposing a constraint on the maximally allowable
carbon emissions. This is particularly useful for companies that seek to
reduce carbon emissions proactively. However, companies may also take a
reactive position and make supply mode decisions based only on inventory and
transport costs. This cost will also include a carbon emission price component
in regions where emissions are subject to carbon pricing mechanisms such as
carbon crediting or taxing. We now briefly show that our model and analysis
also apply to that setting. To that end, let $c^{e}$ denote the price of one
unit of carbon emissions. This price can also depend on the transport mode $m$
and/or product $j$, but for ease of exposition we omit those dependencies. The
long-run average costs per period in 2 can now be redefined as follows:

$\displaystyle\tilde{C}(\mathbf{S}_{f},\mathbf{\Delta})$
$\displaystyle=\sum\limits_{j\in\mathnormal{J}}\left(h_{j}\mathbb{E}[(I_{j}-D_{j})^{+}]+p_{j}\mathbb{E}[(D_{j}-I_{j})^{+}]+\sum\limits_{m\in\mathnormal{M}}(c_{j,m}+c^{e}e_{j,m})\mathbb{E}[Q_{j,m}]\right).$

The decision problem is now to minimize the long-run average costs per period:

$\displaystyle(\tilde{P})\qquad\qquad$ $\displaystyle\displaystyle\min$
$\displaystyle\tilde{C}(\mathbf{S}_{f},\mathbf{\Delta})$ subject to
$\displaystyle\mathbf{S}_{f}\in\mathbb{R}^{|J|},\quad\mathbf{\Delta}\in\mathbb{R}_{0}^{|J|}.$

Problem $(\tilde{P})$ is less complex than the original decision problem since
it does not involve a constraint that links the individual products. As such,
Problem $(\tilde{P})$ can be decomposed in $|J|$ product specific problems,
each of which can be solved individually. The column generation sub-problem
that we will discuss in Section 4.2 has a similar structure as the product
specific problems of $(\tilde{P})$, and the solution method we discuss there
thus also applies to $(\tilde{P})$.

##  4 Analysis

This section focuses on finding the optimal control policy for Problem $(P)$.
Our approach relies on the technique of column generation – also named
Dantzig-Wolfe decomposition after its pioneers (Dantzig and Wolfe 1960). This
technique enables a natural decomposition of the original multi-product
decision problem into smaller single-product problems that have more
structure. We refer the interested reader to Lübbecke and Desrosiers (2005)
for a comprehensive survey on column generation. Below, we first explain how
we apply column generation to Problem $(P)$, and we then describe a
simulation-based optimization method for solving the sub-problem of this
column generation procedure.

###  4.1 The column generation procedure

We first reformulate decision problem $(P)$ as an integer linear program in
which each binary decision variable corresponds to a certain combination of
values for the decision variables of our original decision problem. We
subsequently relax the integrality constraint and we call this problem the
master problem $(\mathit{MP})$. Formally, let $K_{j}$ be the set of all
possible dual index policies for product $j\in\mathnormal{J}$. Each policy
$k\in K_{j}$ is determined by its policy parameters $S^{k}_{j,f}$ and
$\Delta^{k}_{j}$. Let $x^{k}_{j}\in\\{0,1\\}$ denote the decision variable
that indicates whether policy $k\in K_{j}$ is selected ($x^{k}_{j}=1$) for
product $j\in\mathnormal{J}$ or not ($x^{k}_{j}=0$). By relaxing the
integrality constraint on this binary decision variable, we arrive at the
mathematical formulation of the master problem $(\mathit{MP})$:

$\displaystyle(\mathit{MP})\quad\quad$ $\displaystyle\displaystyle\min$
$\displaystyle\sum\limits_{j\in\mathnormal{J}}\sum\limits_{k\in\mathnormal{K}_{j}}C_{j}(S^{k}_{j,f},\Delta_{j}^{k})x_{j}^{k}$
(3) subject to
$\displaystyle\sum\limits_{j\in\mathnormal{J}}\sum\limits_{k\in\mathnormal{K}_{j}}\mathit{E}_{j}(S_{j,f}^{k},\Delta^{k}_{j})x_{j}^{k}\leq\mathcal{E}^{max},$
(4) $\displaystyle\sum\limits_{k\in\mathnormal{K}_{j}}x_{j}^{k}=1,$
$\displaystyle\quad\forall j\in\mathnormal{J}$ (5) $\displaystyle\mbox{
}x_{j}^{k}\geq 0,$ $\displaystyle\quad\forall j\in\mathnormal{J},\forall
k\in\mathnormal{K}_{j}.$

Let $C^{\mathit{LB}}_{P}$ denote the optimal cost for master problem
$\mathit{(MP)}$. Due to the linear relaxation of the integrality constraint on
$x^{k}_{j}$, an optimal cost $C^{\mathit{LB}}_{P}$ is also a lower bound on
the optimal cost for Problem $(P)$, $C_{P}$.

Due to its large number of decision variables, master problem $(\mathit{MP})$
is solved using column generation. To this end, we first restrict master
problem $(\mathit{MP})$ to a small subset $\tilde{K}_{j}\subseteq K_{j}$ of
trivial policies per product $j\in J$ (i.e. columns) that are feasible for
Problem $(P)$ (and thus also for Problem $(\mathit{MP})$). Such a trivial
policy prescribes, for instance, to ship orders exclusively with the least
polluting mode. This restricted problem is referred to as the restricted
master problem $(\mathit{RMP})$. We then solve $(\mathit{RMP})$ to optimality,
and we are interested in new policies $k\in K_{j}\setminus\tilde{K}_{j}$,
$j\in\mathnormal{J}$, that will improve the objective value of
$(\mathit{RMP})$ if they are added to $\tilde{K}_{j}$. Such policies $k\in
K_{j}\setminus\tilde{K}_{j}$ are identified through solving a column
generation sub-problem for each product $j\in J$. The objective function of
such a sub-problem is the reduced cost as a function of the policy with
respect to the current dual variables obtained through solving
$(\mathit{RMP})$ to optimality. If a policy $k\in K_{j}\setminus\tilde{K}_{j}$
has a negative reduced cost, then adding that policy as a column to
$\tilde{K}_{j}$ in $(\mathit{RMP})$ will reduce the objective value of
$(\mathit{RMP})$. More formally, the column generation sub-problem for product
$j\in J$ has the following form:

$\displaystyle(\mathit{SP(j)})~{}$
$\displaystyle~{}\displaystyle\min~{}h_{j}\mathbb{E}[(I_{j}-D_{j})^{+}]+p_{j}\mathbb{E}[(D_{j}-I_{j})^{+}]+\sum_{m\in\mathnormal{M}}(c_{j,m}-\eta
e_{j,m})\mathbb{E}[Q_{j,m}]-\upsilon_{j},$ $\displaystyle\text{subject
to}\quad{S}_{j,f}\in\mathbb{R},\quad\Delta_{j}\in\mathbb{R}_{0},$

where $\eta$ denotes the dual variable of $(\mathit{RMP})$ that corresponds
with emission constraint (4) and $\upsilon_{j}$ denotes the dual variable of
$(\mathit{RMP})$ that corresponds with constraint (5) that assures that for
each product $j\in\mathnormal{J}$ a convex combination of policies is chosen.
Note that these dual variables can also be interpreted as the Lagrange
multipliers of relaxing the corresponding constraints (Lübbecke and Desrosiers
2005). If for product $j\in J$ a feasible solution for $(\mathit{SP(j)})$
exists with a negative objective value (i.e. a negative reduced cost), then
this policy is added to $\tilde{K}_{j}$ since the objective value of
$(\mathit{RMP})$ can be improved when solved with the enlarged set
$\tilde{K}_{j}$.

We continue with iterating between optimizing $(\mathit{RMP})$ and finding new
policies through solving $(\mathit{SP(j))}$, $j\in\mathnormal{J}$, until no
product for which there is a policy with a negative reduced cost exists. An
optimal solution for $(\mathit{RMP})$ is then also an optimal solution for
$(\mathit{MP})$. If this optimal solution contains integer values only, then
it is also an optimal solution for $(P)$. If this is not the case, then we
solve $(\mathit{RMP})$ one last time as an integer program to find an integer
solution for $(\mathit{RMP})$, which is then also a feasible solution for
$(P)$. Recent inventory literature has shown that solving the restricted
master program as an integer program to arrive at an integer solution leads to
good performance in terms of optimality gaps (e.g., Drent and Arts 2021,
Haubitz and Thonemann 2021), and often outperforms alternative approaches such
as local searches or rounding procedures (Alvarez et al. 2013). The
corresponding cost of the resulting feasible solution is also an upper bound,
denoted $C_{\mathit{UB}}^{P}$, for $C^{P}$.

In the next section, we provide a simulation-based optimization procedure to
solve the column generation sub-problem $(\mathit{SP}(j))$.

###  4.2 Solving the column generation sub-problem

The column generation sub-problem $(\mathit{SP(j)})$ has the same structure as
the problem studied by Veeraraghavan and Scheller-Wolf (2008). We follow their
simulation-based optimization procedure to solve $(\mathit{SP(j)})$. This
procedure is grounded in the following separability result that allows us to
find the optimal $S_{j,f}$ for a given $\Delta_{j}$ as the solution to a
special Newsvendor problem.

###### Lemma 4.1

(Veeraraghavan and Scheller-Wolf 2008, Proposition 4.1) The distributions of
the overshoot $O_{j}$, the fast transport mode shipment size $Q_{j,f}$, and
the slow transport mode shipment size $Q_{j,s}$ are functions of $\Delta_{j}$
only, independent of $S_{j,f}$.

Let $O^{t}_{j}(\Delta_{j})$ denote the overshoot of product $j$ in period $t$
for a given $\Delta_{j}$. Recall that the fast inventory position of product
$j$ in period $t$ after shipping equals $S_{j,f}+O^{t}_{j}(\Delta_{j})$.
Consequently, for the net inventory level of product $j$ in each period $t$,
we can also write

$I_{j}^{t}=S_{j,f}-\left(\sum\nolimits_{k=t-l_{j,f}}^{t-1}D_{j}^{k}-O^{t-l_{j,f}}_{j}(\Delta_{j})\right).$
(6)

By plugging (6) in the objective function of $(\mathit{SP(j)})$, and using
Lemma 4.1 as well as the fact that in each period the overshoot is independent
of the demand and $S_{j,f}$, we readily recognize that for given $\Delta_{j}$
the objective function is convex in $S_{j,f}$. This implies the following
result.

###### Lemma 4.2

(Veeraraghavan and Scheller-Wolf 2008, Theorem 4.1) The optimal base-stock
level $S_{j,f}^{*}$ for a given $\Delta_{j}$, denoted
$S_{j,f}^{*}(\Delta_{j})$, equals

$S_{j,f}^{*}(\Delta_{j})=\inf\left\\{S_{j,f}\in\mathbb{R}:\mathbb{P}\left(\sum\nolimits_{i=1}^{l_{j,f}+1}D_{i}-O_{j}(\Delta_{j})\leq
S_{j,f}\right)\geq\frac{p_{j}}{p_{j}+h_{j}}\right\\}.$

It now remains to calculate the objective value of $\mathit{SUB}(j)$ for given
$\Delta_{j}$ and corresponding $S_{j,f}^{*}(\Delta_{j})$. Observe that in each
period immediately after shipping orders with both modes, the slow inventory
position equals the fast inventory position plus the overshoot and all
remaining outstanding slow orders. Since these inventory positions are equal
to their respective base-stock levels following order placement, we have for
each product $j$ in each period $t$:

$S_{j,s}=S_{j,f}+O^{t}_{j}+\sum\limits_{k=0}^{l_{j}-1}Q^{t-k}_{j,s}.$ (7)

From (7) it follows that
$\mathbb{E}[Q_{j,s}]=(\Delta_{j}-\mathbb{E}[O_{j}])/l_{j}$. Since under
backlogging the sum of both orders must on average be equal to the period
demand, we finally find
$\mathbb{E}[Q_{j,f}]=\mathbb{E}[D_{j}]-\mathbb{E}[Q_{j,s}]$.

To solve $(\mathit{SP(j)})$, $j\in\mathnormal{J}$, to optimality, it thus
suffices to perform a one-dimensional search over $\Delta_{j}$. For each
$\Delta_{j}$, we compute the stationary distribution of the overshoot. With
this stationary distribution we readily find the optimal base-stock level
$S_{j,f}^{*}(\Delta_{j})$ through Lemma 4.2 and the total reduced cost through
the identities following Equation $\eqref{eq:totalOrders}$. As there is in
general no closed-form expression for the stationary distribution of the
overshoot, we follow Veeraraghavan and Scheller-Wolf (2008) and rely on
simulation to compute this distribution.

Note that our optimization model and analysis readily extends to settings
where the slow transport modes of all (or some) products are operated
according to any other rule that depends only on the current overshoot as well
as all in-transit orders that are not yet included in the fast inventory
position. That is, any other rule that depends only on the information state
$(O_{j}^{t},Q^{t-1}_{j,s},Q^{t-2}_{j,s},\ldots,Q^{t-l_{j}+1}_{j,s})$,
$j\in\mathnormal{J}$, $t\in\mathbb{N}_{0}$. Most well-performing control
policies satisfy this condition, e.g., the Capped Dual-Index policy (Sun and
Van Mieghem 2019), the Tailored Base-Surge policy (Allon and Van Mieghem
2010), and the Projected Expedited Inventory Position policy (Drent and Arts
2022). Sheopuri et al. (2010) show that for such control policies, the
stationary distribution of the overshoot is a function of only the
parameter(s) for operating the slow transport mode, and that consequently a
Newsvendor result similar to Lemma 4.2 holds for all such policies.

##  5 Computational experiment

Our test-bed has three different types of assortments of products, each
representing a different type of industry. The first assortment type consists
solely of products for which emissions from the slowest transport mode are
less than the emissions from the faster transport mode. This assortment is
inspired by apparel goods that are delivered from Vietnam to Europe by sea
transport as the slowest mode and by air transport as the fastest mode. In
this example, the fast supply mode has a higher carbon footprint. The opposite
holds for the second assortment type. Here we are inspired by industrial goods
that are delivered from China to Europe by sea transport as the slowest mode
and from Germany by truck as the fastest mode. In this case, the slow supply
mode is associated with more emissions as goods are transported over a longer
distance. The third assortment has products of both types.

We perform a parametric computational experiment. The base case is set up as
follows. We consider 100 products for each assortment type, i.e. $|J|=100$.
The input parameters in the base case are identical for all three assortments,
except for the carbon emissions from transportation. For each product $j\in
J$, the period demand $D_{j}$ follows a negative binomial distribution. To
create heterogeneous assortments, the parameters of this negative binomial
distribution are randomly drawn from two separate distributions for each
product $j$. The mean $\mu_{D_{j}}:=\mathbb{E}[D_{j}]$ is randomly drawn from
a gamma distribution with mean 100 and coefficient of variation of 0.5. The
coefficient of variation $CV_{D_{j}}:=\sqrt{\mbox{\bf
Var}[D_{j}]}/\mu_{D_{j}}$ is randomly drawn from a shifted beta distribution
with mean 0.9, standard deviation of 0.25, and shifted to the right by 0.3.
Since low demand products typically have higher holding cost, the holding cost
$h_{j}$ is negatively correlated with the mean demand $\mu_{D_{j}}$ of each
product $j$ through a Gaussian copula with a fitting covariance matrix. In
particular, $h_{j}$ is drawn from a gamma distribution with mean 1 and
coefficient of variation equal to $CV_{D_{j}}$, with a Pearson correlation
coefficient of -0.5. Details regarding our approach to generate correlated
random numbers are relegated to Appendix 8.

We set $l_{j,f}$ and $c_{j,s}$ to 0 for all products, and focus on $l_{j,s}$
and $c_{j,f}$, which now coincide with the lead time difference and the cost
premium of product $j\in J$, respectively. We set $l_{j,f}=3$ for all
products. The back-order penalty cost $p_{j}$ for product $j$ is a function of
its holding cost $h_{j}$. The ratio between the holding and the penalty cost
is an important determinant of the service level in an inventory system.
Therefore we set $p_{j}=\psi_{j}\chi_{j}^{p}h_{j}$ where $\psi_{j}$ is a
parameter we use to control the ratio between $p_{j}$ and $h$, and
$\chi_{j}^{p}$ is a random perturbation. That is, $\chi_{j}^{p}$ has a shifted
beta distribution with mean 0.98, standard deviation 0.1, and shifted 0.02 to
the right.

The cost premium $c_{j,f}$ of the fast transport mode of product $j$ equals
$\chi_{j}^{c}p_{j}l_{j,f}$, where $\chi_{j}^{c}$ has a beta distribution with
mean 0.25 and standard deviation 0.1. Table 2 provides a summary of how we
randomly generated the products of the base case; we shortly explain how we
randomly generated the emission units for these products for all the three
assortment types. In Table 2, $NB(\mu,cv)$ denotes a negative binomial random
variable with mean $\mu$ and coefficient of variation $cv$, $\Gamma(\mu,cv)$
denotes a gamma random variable with mean $\mu$ and coefficient of variation
$cv$, and $B(\mu,\sigma,s)$ denotes a beta random variable with mean $\mu$ and
standard deviation $\sigma$ that is shifted to the right by $s$; if we drop
$s$ then this beta random variable is not shifted, i.e.
$B(\mu,\sigma,s)=_{d}B(\mu,\sigma)+s$ where $=_{d}$ denotes equality in
distribution.

Table 2: Generating the base case input parameters. Input parameter | Generation
---|---
$|\mathnormal{J}|$ | 100
$D_{j}$ |  $\mathnormal{NB}(\mu_{D_{j}},\text{CV}_{D_{j}})$, with $\mu_{D_{j}}\sim\Gamma(100,0.5)$ and $\text{CV}_{D_{j}}\sim\mathnormal{B}(0.9,0.25,0.3)$
$l_{j,f}$ | 0
$l_{j,s}$ | 3
$h_{j}$ |  $\Gamma(1,0.5)$, with $\rho_{\mu_{D_{j}},h_{j}}=-0.5$
$p_{j}$ |  $\psi_{p}\chi_{j}^{p}h_{j}$, with $\psi_{p}=9$ and $\chi_{j}^{p}\sim\mathnormal{B}(0.98,0.1,0.02)$
$c_{j,s}$ | 0
$c_{j,f}$ |  $\chi_{j}^{c}p_{j}l_{j,f}$, with $\chi_{j}^{c}\sim\mathnormal{B}(0.25,0.1)$

As explained in Section 3, we rely on the NTM framework (NTM 2015) to set
emissions based on the structure of equation (1). We apply the NTM framework
to data from the UN Comtrade Database (United Nations Comtrade Database 2020)
to obtain sample emission units; details regarding this methodology are
relegated to appendix 7. We then apply maximum likelihood estimation on these
sample unit emissions to obtain three distinct sets of two distribution
functions; two for each assortment type of products. The number of carbon
emission units $e_{j,m}$ for product $j\in J$ with transport mode $m\in M$ are
then randomly drawn from these distributions. These distributions are
presented in Table 3. In this table, $LN(\mu,\sigma)$ denotes a random
variable whose logarithm is normally distributed with mean $\mu$ and standard
deviation $\sigma$, and $WB(\lambda,k)$ denotes a Weibull random variable with
scale $\lambda$ and shape $k$.

Table 3: Generating the emission units for the base case. Assortment type | Emission parameter | Generation
---|---|---
1 | $e_{j,s}$ | $\Gamma(0.35,0.21)$
$e_{j,f}$ |  $e_{j,s}$ \+ $\mathnormal{LN}(1.52,0.21)$
2 | $e_{j,f}$ | $\Gamma(0.19,1.27)$
$e_{j,s}$ |  $e_{j,f}$ \+ $\Gamma(2.19,1.27)$
3 | $e_{j,f}$ | $\mathnormal{WB}(0.87,0.77)$
$e_{j,s}$ | $\Gamma(3.31,1.34)$

The total allowable carbon emissions from transportation
$\mathcal{E}^{\text{max}}$ is set as a percentage of the total reducible
carbon emissions. For each instance of the test-bed, the reducible carbon
emissions is defined as the difference between the total amount of carbon
emissions of the control policy that is optimal for Problem $(P)$ absent of
the emission constraint and the total amount of carbon emissions of the
control policy that leads to the lowest possible total carbon emissions. The
latter implies that each product is only shipped with its least polluting
transport mode. Under the dual-index policy, setting $\Delta_{j}$ to zero
implies that all orders for product $j\in J$ are shipped with its fastest
transport mode. Alternatively, letting $S_{j,f}$ go to $-\infty$ implies that
all orders for product $j$ are shipped with its slowest transport mode
(Veeraraghavan and Scheller-Wolf 2008).

To evaluate the effectiveness of the column generation procedure in solving
Problem $(P)$, we compute for each instance of the test-bed the relative
difference between the total average cost under a feasible solution and the
corresponding lower bound. That is,

$\%\mathit{GAP}=100\cdot\frac{C^{\mathit{UB}}_{P}-C^{\mathit{LB}}_{P}}{C^{\mathit{LB}}_{P}},$

where $C^{\mathit{LB}}_{P}$ and $C^{\mathit{UB}}_{P}$ are obtained using the
methods described in Section 4. In what follows, we also refer to this
feasible solution as dynamic mode selection (DMS). Hence the long run average
cost of dynamic mode selection equals $C^{\mathit{UB}}_{P}$.

To quantify the benefit of using two transport modes dynamically rather than
relying statically on one transport mode, we define for each instance of the
test-bed a benchmark instance in which we can only select one transport mode
for each product. As described above, the dual-index policy can mimic static
mode selection (SMS). Hence, to find a feasible static mode selection solution
to this benchmark instance of Problem $(P)$, we apply our column generation
procedure of Section 4 in which we restrict the solution space for each
product $j\in\mathnormal{J}$ such that all orders are shipped with either the
fastest or the slowest transport mode. The mathematical formulation for the
static mode selection approach as well as the benchmark approach described in
the next paragraph can be found in Appendix 9. (Note that in the Master
Problem $(\mathit{MP})$ for this approach the set of possible policies $K_{j}$
for each product $i\in\mathnormal{J}$ contains only two single transport mode
policies.) The long run average cost of this solution is denoted
$C^{\mathit{SMS}}_{P}$. To quantify the value of dynamic mode selection, we
compare for each instance of the test-bed the long run average cost of static
mode selection with the long run average cost of dynamic mode selection. That
is,

$\%\mathit{SMS}=100\cdot\frac{C^{\mathit{SMS}}_{P}-C^{UB}_{P}}{C^{\mathit{UB}}_{P}},$

where $\%SMS$ indicates the relative increase in the long run average cost
when the company chooses to rely on only one transport mode for each product
in meeting an assortment wide constraint on total emissions.

To quantify the portfolio effect, we define for each instance of the test-bed
an additional benchmark in which we enforce an emission constraint for each
product $j\in J$. This emission constraint is set as a percentage of the total
reducible emissions of the individual product rather than of the entire
assortment of products. We compute a feasible solution for this benchmark
instance with the column generation procedure of Section 4, which can readily
be modified so that it can be applied to settings where we have additional
emission constraints. (An alternate solution procedure is to use a simple
search procedure per product.) This approach, which we refer to as blanket
mode selection (BMS), seems a plausible approach for most practitioners.
Indeed, they can consider each product individually absent of the complicating
linking emission constraint and yet they are guaranteed that the total
emissions of the entire assortment is kept below the target level. The long
run average cost of the blanket mode selection approach is denoted
$C^{BMS}_{P}$. To quantify the portfolio effect, we compare for each instance
of the test-bed the long run average cost of blanket mode selection with the
long run average cost of dynamic mode selection. That is,

$\%\mathit{BMS}=100\cdot\frac{C^{BMS}_{P}-C^{UB}_{P}}{C^{UB}_{P}},$

where $\%\mathit{BMS}$ indicates the relative increase in the long run average
cost when the company enforces an emission constraint on each individual
product rather than one single constraint for the emissions from the entire
assortment.

In solving the column generation sub-problem for each product, we simulated 10
samples of 9500 time periods following a 5000 time periods warm-up . The width
of the 95 percent confidence interval of the long run average cost per period
for each product was no larger than 3 percent of its corresponding point
estimate for each instance of the column generation sub problem that we
solved. The average computational time of our column generation procedure was
23 minutes, 5 minutes for blanket dynamic mode selection, and less than 5 sec
for static mode selection.

###  5.1 Results for the base case

Figure 1 presents the normalized optimal average costs of each approach for
each assortment group under emission targets that range from from 0 to 100
percent of the total reducible carbon emissions. Observe that for each
assortment group, all approaches have the same performance when the emission
target is set at 100 percent of the total reducible emissions. In this case,
all approaches solely utilize the least polluting transport mode.
Alternatively, when we impose no target on the emissions from transportation,
then both dual mode approaches perform equally well while the static mode
selection approach seems to perform the poorest over all possible emission
targets. Indeed, static mode selection is around 15 percent more expensive
than both dual dual mode approaches for all assortment types when
transportation emissions are not constrained. Based on Figure 1, we conclude
that dynamic mode selection, as opposed to static and blanket mode selection,
has great potential to efficiently curb carbon emissions from transportation
at relatively little additional costs across all assortment types.

Figure 1: Optimal normalized absolute costs of each approach for different
targets on the reducible emissions.

We explicitly compare our dynamic mode selection with the benchmark approaches
in Figure 2, which presents the $\%\mathit{SMS}$ and $\%\mathit{BMS}$
percentages for each assortment group under emission targets that range from
from 0 to 100 percent of the total reducible carbon emissions. The figure
indicates that the performance of static mode selection over dynamic mode
selection is consistent across all assortments. The relative increase in its
total cost over dynamic mode selection is the largest when there is no
emission target, and gradually decreases as the emission constraint tightens.
At moderate carbon emission targets, around 40 to 60 percent of the total
reducible emissions, static mode selection still leads to increases in the
total average cost per period of around 10 to 15 percent for all assortment
types.

The performance of the blanket mode selection approach depends on the specific
assortment type. Figure 2(a) illustrates that when the unit emissions from the
fast transport mode are more than those from the slow transport mode, the
performance of the blanket mode selection approach seems to be quite
reasonable. This can be explained as in this setting, the cost of the fastest
transport mode is larger than the cost of the slowest transport mode. The most
polluting transport mode is thus also the most expensive transport mode. The
portfolio effect is limited for this assortment type.

For the other two assortment types, however, the cheapest transport mode is
not necessarily also the least polluting transport mode, and the portfolio
effect is more prevalent. Figures 2(b) and 2(c) show that $\%\mathit{BMS}$ can
be more than 35 and 20 percent in assortment type 2 and 3, respectively, under
carbon emission reduction targets of 50 percent. The static mode selection
approach, which also takes advantage of the portfolio effect, even outperforms
the blanket mode selection approach for quite some emission targets.

Figure 2: Relative surplus of the optimal cost of the alternative approaches
(BMS and SMS) compared to the DMS approach for different targets on the
reducible emissions.

Figure 3 shows the usage of the fast supply mode as a function of the carbon
emission reduction target under dynamic mode selection as measured by

$\%F=100\cdot\frac{1}{|J|}\sum_{j\in
J}\frac{\mathbb{E}[Q_{j,f}]}{\mathbb{E}[Q_{j,f}]+\mathbb{E}[Q_{j,s}]}.$

Dynamic mode selection is economically attractive regardless of emissions
targets. Carbon reduction targets make dynamic mode selection even more
attractive for assortment type 2, and mixed assortments, but not for
assortments of type 1.

Figure 3: Percentage of products shipped with the fast mode for different
carbon emission reduction targets.

To recapitulate, the value of dynamically shipping products with two transport
modes simultaneously rather than statically selecting one transport mode a
priori is quite large. Regardless of the assortment type, $\%\mathit{SMS}$ is
in between 5 and 15 percent for emission reduction targets up to 90 percent.
The portfolio effect depends on the specific assortment type. If the least
polluting transport mode of each product is also its cheapest transport mode,
then the fastest and most polluting transport modes are typically only relied
upon in case of imminent backorders. This behavior remains in case of an
assortment-wide emission target, and the portfolio effect is consequently
rather limited. If the least polluting transport modes are not necessarily the
cheapest transport modes, then there is substantial value to be reaped in
optimizing the assortment of products under a single emission constraint
rather than under separate emission constraints for each individual product.
Indeed, $\%\mathit{BMS}$ can go up to 40 and 20 percent for assortment type 2
and 3, respectively.

Table 4 below presents the average relative slack in the emission constraints
for each assortment over the different emission targets considered in the base
case analysis. The table shows that due to the binary nature of the static
mode selection approach, the total average emissions under this approach are
often substantially lower than the target level. This leads to particularly
poor performance for assortment type 1. We observed in our computational
experiments that for this assortment type, the static mode selection approach
selects the cheapest and thus least polluting transport mode for almost all
products under each emission target.

Table 4: Average slack in emission constraints for the base case analysis. | Assortment type
---|---
Approach | 1 | 2 | 3
DMS | 0.00% | 0.00% | 0.00%
BMS | 0.08% | 0.61% | 0.12%
SMS | 12.16% | 0.17% | 0.97%

The average $\%GAP$ of the base case over all emission targets is less than
0.01 percent, indicating that the column generation procedure finds feasible
solutions that are close to optimal. Such a low average $\%GAP$ occurs because
there can be at most 1 product for which the optimal solution to Problem
$(MP)$ is fractional. Indeed, Problem $(MP)$ has $|J|+1$ constraints and an
optimal solution for this problem has the same number of basic variables.
Constraint (5) assures that for each product $j\in J$ a convex combination of
policies is chosen. As such, there is at least one basic variable for each
product $j$. This implies that there is at most 1 product for which the
optimal solution to Problem $(MP)$ is fractional.

###  5.2 Determinant of emission reduction potential

The base case analysis of our DMS approach indicates that the emission
reductions are not evenly distributed across the products. Products can be
ordered by their contribution to emission reduction following the DMS
optimization. In this manner we can construct the cumulative reduction in
emissions as shown in the Lorenz curves (Figure 4) with the dashed line.
Figure 4 shows that 20% of the products in assortments 1 through 3 account for
61.22%, 94.19%, and 91.88% of the emission reduction, respectively. This
suggests that most of the emission reduction can be achieved by using dynamic
mode selection for a limited subset of a given assortment. Although it is
possible to determine the limited subset of products that account for most of
the emission reduction after performing the DMS optimization, it would be
convenient to know which products to focus on without having to solve a
sophisticated optimization problem. Suppose we order products in increasing
order of $\frac{|e_{j,f}-e_{j,s}|}{c_{j,f}-c_{j,s}}$, i.e., we sort products
according to the how much emission can be saved by using the least polluting
transport mode relative to the additional cost of the faster transport mode.
Figure 4 shows the cumulative emission saving by products ordered this way
with the solid line. Here we see that focusing on the 20% of products in
assortments 1 through 3 for which $\frac{|e_{j,f}-e_{j,s}|}{c_{j,f}-c_{j,s}}$
is highest, already achieves 58.93%, 93.63%, and 82.64% of the potential
emission reduction, respectively. Thus firms seeking to minimize the carbon
footprint of their inbound logistics should focus their attention on products
for which the difference in emission in different transport modes is large
relative to the additional cost of fast transportation modes. That is emission
differences between modes relative to cost difference between modes is the
main determinant of emission reduction potential for a given assortment.

Figure 4: The cumulative emissions reduction share of items arranged based on
different criteria.

###  5.3 Comparative statics

In this section, we study how changes in the input parameters with respect to
the base case affect the performance of the blanket mode selection and the
static mode selection approach. In what follows, we keep the emission target
fixed at a 50 percent reduction of the reducible emissions, and we study the
effects of changing a certain input parameter while generating the other input
parameters as in the base case, i.e. as in Table 2. We also investigate the
effect of scaling the emission differences between the most polluting and
least polluting transport modes when the target on the total emissions is kept
fixed. To achieve this, we first generate emission units as in the base case.
We subsequently change the emission units of the most polluting transport mode
through scaling $\mid e_{j,s}-e_{j,f}\mid$ by a constant $\delta_{e}$ while
keeping the emission units from the least polluting mode fixed at its base
level. The changes in the parameters we investigate are summarized in Table 5.

Table 5: Changes in the base case input parameters. Parameter | Generation | Base case | Changes
---|---|---|---
$|\mathnormal{J}|$ |  | 100 | $|\mathnormal{J}|\in\\{40,60,80\\}$
$\mu_{D_{j}}$ | $\Gamma(100,CV_{\mu_{D_{j}}})$ | 0.5 | $CV_{\mu_{D_{j}}}\in\\{0.3,0.4,0.6,0.7\\}$
$CV_{D_{j}}$ | $B(0.9,0.25,s_{D_{j}})$ | 0.3 | $s_{D_{j}}\in\\{0.2,0.25,0.35,0.4\\}$
$\rho_{\mu_{D_{j}},h_{j}}$ |  | -0.5 | $\rho_{\mu_{D_{j}},h_{j}}\in\\{-0.3,-0.4,-0.6,-0.7\\}$
$\psi_{p}$ |  | 9 | $\psi_{p}\in\\{3,4,5,19,99\\}$
$\chi_{j}^{c}$ | $B(0.25,\sigma_{c})$ | 0.1 | $\sigma_{c}\in\\{0.15,0.2,0.3,0.35\\}$
$l_{j,s}$ |  | 3 | $l_{j,s}\in\\{2,4\\}$
$\delta_{e}$ |  | 1 | $\delta_{e}\in\\{0.8,0.9,1.1,1.2\\}$

Figure 5 shows the effect of changing the coefficient of variation of the
gamma distribution from which we sample the mean demands per period for each
product. The figure indicates that this effect is relatively limited. With
respect to the base case, both $\%\mathit{BMS}$ and $\%\mathit{SMS}$ only
change up to 1 percent point for all three assortment types. We can draw a
similar conclusion for the effect of changing the Pearson correlation
coefficient between the holding cost and the mean demand per period of each
product. Figure 6 shows that both $\%\mathit{BMS}$ and $\%\mathit{SMS}$ change
at most 1 percent point with respect to the base case for all three assortment
types.

Figure 5: Effect of changing $CV_{\mu_{D_{j}}}$ while keeping the rest of the
parameters as in the base case.
Figure 6: Effect of changing $\rho_{\mu_{D_{j}},h_{j}}$ while keeping the rest
of the parameters as in the base case.

Alterations in the shift parameter of the beta distribution from which we
sample the coefficient of variation of the demand per period for each product
has a relatively moderate effect on both $\%\mathit{BMS}$ and
$\%\mathit{SMS}$. Figure 7 illustrates that for all assortment types, the
$\%\mathit{SMS}$ tends to increase in the variability of the demand while the
$\%\mathit{BMS}$ decreases. This indicates that the flexibility to dynamically
ship products with two transport modes has particular merit in highly variable
demand settings.

Figure 7: Effect of changing $s_{D_{j}}$ while keeping the rest of the
parameters as in the base case.

Figure 8 indicates that for all assortment types, an increase (decrease) in
the lead time difference between the fastest and the slowest transport modes
of each product leads to an increase (decrease) in both $\%\mathit{SMS}$ and
$\%\mathit{BMS}$. The blanket mode selection approach seems to be more
susceptible to changes in the lead time difference than the static mode
selection approach. For assortment type 2, for instance, an increase in the
lead time difference to 4 leads to an increase in $\%\mathit{BMS}$ of 7
percent points with respect to the base case. By contrast, $\%\mathit{SMS}$
increases only slightly by 0.5 percent points. This can be attributed to the
fact that the blanket mode selection approach imposes constraints on the
emissions of each individual product while the static mode selection approach
imposes a single constraint on the entire assortment of products.

Figure 8: Effect of changing $l_{j,s}$ while keeping the rest of the
parameters as in the base case.

Figure 9 illustrates the effect of changing the critical ratio for all
products through varying $\psi_{p}$. We conclude that this effect is quite
large. For assortment type 1, for instance, $\%\mathit{SMS}$ varies from 20
percent to 2 percent. While $\%\mathit{SMS}$ seems to decrease in the critical
ratio for all products, $\%\mathit{BMS}$ tends to increase. For assortment
type 2 and 3, for instance, $\%\mathit{BMS}$ increases from 15 percent to over
80 percent. These effects can be explained by the fact that as the critical
ratios of all products approach 1, our dynamic mode selection approach will
mimic the static mode selection approach in which the assortment wide emission
constraint is met by relying on less polluting transport mode for product for
which this is relatively cheap to do so. The blanket mode selection approach,
however, must meet emission targets for each product individually which leads
to poor performance if we increase the critical ratios for all products.

Figure 9: Effect of changing $\psi_{p}$ while keeping the rest of the
parameters as in the base case.

Figure 10 indicates that $\%\mathit{SMS}$ decreases in the cost of the fast
transport mode. This can again be explained by the fact that our dynamic mode
selection approach will also rely more on the cheaper transport mode as the
cost premium for the fast transport mode increases, and that consequently the
gap with the static mode selection approach decreases. By contrast,
$\%\mathit{BMS}$ increases in the cost of the fast transport mode for
assortment type 2 and 3. This can be attributed to the fact that the blanket
mode selection approach, contrary to the other two approaches, imposes
itemized emission constraints and relying on the most expensive but least
polluting transport mode is therefore inevitable. Note that this is not true
for assortment type 1 because there the fast, expensive transport mode is also
the most polluting transport mode. We can draw similar conclusions for the
effects scaling the emission units of the most polluting transport mode, see
Figure 11. Finally, Figure 12 indicates that the impact of the assortment size
on both $\%\mathit{SMS}$ and $\%\mathit{BMS}$ is relatively limited. With
respect to the base case, both $\%\mathit{BMS}$ and $\%\mathit{SMS}$ only
change up to 2 percent point for all three assortment types.

Figure 10: Effect of changing $\chi_{j}^{c}$ while keeping the rest of the
parameters as in the base case.
Figure 11: Effect of changing $\delta_{e}$ while keeping the rest of the
parameters as in the base case.
Figure 12: Effect of changing $|J|$ while keeping the rest of the parameters
as in the base case.

##  6 Concluding remarks

As carbon emissions from the transportation sector are projected to increase
over the next decades, it is important for companies to rethink their supply
chain strategies and explicitly incorporate carbon emissions into their
decision making. In this paper, we have studied the inbound transport and
inventory management decision making for a company that sells an assortment of
products. The company wishes to minimize inventory costs while keeping the
total emissions from the inbound transport of the entire assortment below a
certain target level. Each product can be shipped using two distinct transport
modes. As each mode has its own merits, we have proposed a dynamic mode
selection model that allows the company to ship products with either mode
depending on when one mode is more favorable than the other. Since the optimal
policy for dual transport mode problems are known to be complex, we have
assumed that shipment shipment quantities for each product are governed by a
dual-index policy. We have formulated the resulting decision problem as a
mixed integer linear program that we have solved through a column generation
solution procedure. This column generation procedure decomposes the complex
multi-product problem into smaller sub-problems per product. These sub-
problems are readily solved through a simple bisection search over Newsvendor
type problems.

In an extensive computational experiment, we have compared the performance of
our dynamic mode selection approach with two alternative approaches that are
considered state of the art. The first benchmark, static mode selection, lacks
the flexibility to dynamically ship products with two transport modes; it
rather selects one transport mode for each product a priori. The second
benchmark, blanket mode selection, does have the flexibility to rely on two
transport modes simultaneously but it makes transport decisions for each
product individually rather than holistically for the entire assortment. Our
computational experiments indicate that the value of our dynamic mode
selection approach over the blanket mode selection approach is particularly
high for assortments of products for which the fastest transport modes are not
necessarily the most polluting transport modes. For such settings, our dynamic
mode selection approach can reduce the long run average costs by 40 percent
under moderate carbon emission targets. These huge savings can be attributed
to the portfolio effect inherent to our approach. The computational
experiments further indicate that dynamic mode selection can significantly
outperform static mode selection. Under moderate emission targets, dynamically
relying on two transport modes rather than a single transport mode can lead to
cost savings of up to 15 percent.

Future studies can extend the current model by studying other settings with
multiple transport modes such as multiple echelons in a serial system (e.g.
Lawson and Porteus 2000, Arts and Kiesmüller 2013) or assembly systems (e.g.
Angelus and Özer 2016). Alternatively one may consider more sophisticated dual
mode heuristic policies such as the projected expedited inventory position
policy (Drent and Arts 2022), capped base-stock policy (Sun and Van Mieghem
2019), or vector base-stock policy (Sheopuri et al. 2010).

The research of the first author is supported by the National Research Fund of
Luxembourg through AFR grant 12451704. The authors gratefully acknowledge
Ranit Sinha for his comments on earlier versions of the paper.

## References

  * Abnett (2020) Abnett K (2020) World’s largest carbon market faces revamp under draft eu plan. _Reuters_ .
  * Allon and Van Mieghem (2010) Allon G, Van Mieghem J (2010) Global dual sourcing: Tailored base-surge allocation to near-and offshore production. _Management Science_ 56(1):110–124.
  * Alvarez et al. (2013) Alvarez E, Van der Heijden M, Zijm W (2013) The selective use of emergency shipments for service-contract differentiation. _International Journal of Production Economics_ 143(2):518–526.
  * Angelus and Özer (2016) Angelus A, Özer O (2016) Knowledge you can act on: Optimal policies for assembly systems with expediting and advance demand information. _Operations Research_ 64(6):1338–1371.
  * Arts and Kiesmüller (2013) Arts J, Kiesmüller G (2013) Analysis of a two-echelon inventory system with two supply modes. _European Journal of Operational Research_ 225(2):263–272.
  * Arts et al. (2011) Arts J, Van Vuuren M, Kiesmüller G (2011) Efficient optimization of the dual-index policy using markov chains. _IIE Transactions_ 43(8):604–620.
  * Arıkan and Jammernegg (2014) Arıkan E, Jammernegg W (2014) The single period inventory model under dual sourcing and product carbon footprint constraint. _International Journal of Production Economics_ 157:15–23.
  * Barbosa-Póvoa et al. (2018) Barbosa-Póvoa A, da Silva C, Carvalho A (2018) Opportunities and challenges in sustainable supply chain: An operations research perspective. _European Journal of Operational Research_ 268(2):399–431.
  * Berling and Martínez de Albéniz (2016) Berling P, Martínez de Albéniz V (2016) Dynamic speed optimization in supply chains with stochastic demand. _Transportation Science_ 50(3):1114–1127.
  * Brandenburg et al. (2014) Brandenburg M, Govindan K, Sarkis J, Seuring S (2014) Quantitative models for sustainable supply chain management: Developments and directions. _European Journal of Operational Research_ 233(2):299–312.
  * Chen and Wang (2016) Chen X, Wang X (2016) Effects of carbon emission reduction policies on transportation mode selections with stochastic demand. _Transportation Research Part E: Logistics and Transportation Review_ 90:196–205.
  * Dantzig and Wolfe (1960) Dantzig G, Wolfe P (1960) Decomposition principle for linear programs. _Operations Research_ 8(1):101–111.
  * Dekker et al. (2012) Dekker R, Bloemhof J, Mallidis I (2012) Operations research for green logistics – an overview of aspects, issues, contributions and challenges. _European Journal of Operational Research_ 219(3):671–679.
  * Dong et al. (2018) Dong C, Transchel S, Hoberg K (2018) An inventory control model for modal split transport: A tailored base-surge approach. _European Journal of Operational Research_ 264(1):89–105.
  * Dong et al. (2019) Dong Q Cand Li, Shen B, Tong X (2019) Sustainability in supply chains with behavioral concerns. _Sustainability_ 11(15).
  * Drent and Arts (2021) Drent M, Arts J (2021) Expediting in two-echelon spare parts inventory systems. _Manufacturing & Service Operations Management_ 23(6):1431–1448.
  * Drent and Arts (2022) Drent M, Arts J (2022) Effective dual-sourcing through inventory projection. _arXiv preprint 2207.12182_ .
  * Engebrethsen and Dauzére-Pérés (2019) Engebrethsen E, Dauzére-Pérés S (2019) Transportation mode selection in inventory models: A literature review. _European Journal of Operational Research_ 279(1):1–25.
  * European Commission (2020) European Commission (2020) State of the Union: Q & A on the 2030 Climate Target Plan .
  * European Environment Agency (2020) European Environment Agency (2020) _Trends and projections in Europe 2020: Tracking progress towards Europe’s climate and energy targets_.
  * Feng et al. (2006) Feng Q, Sethi S, Yan H, Zhang H (2006) Are base-stock policies optimal in inventory problems with multiple delivery modes? _Operations Research_ 54(4):801–807.
  * Fukuda (1964) Fukuda Y (1964) Optimal policies for the inventory problem with negotiable leadtime. _Management Science_ 10(4):690–708.
  * Gut (2009) Gut A (2009) _An Intermediate Course in Probability_ (Springer), second edition.
  * Haubitz and Thonemann (2021) Haubitz C, Thonemann U (2021) How to change a running system—controlling the transition to optimized spare parts inventory policies. _Production and Operations Management_ 30(5):1386–1405.
  * Hoen et al. (2014a) Hoen K, Tan T, Fransoo J, Van Houtum G (2014a) Effect of carbon emission regulations on transport mode selection under stochastic demand. _Flexible Service Manufacturing Journal_ 26:170–195.
  * Hoen et al. (2014b) Hoen K, Tan T, Fransoo J, Van Houtum G (2014b) Switching transport modes to meet voluntary carbon emission targets. _Transportation Science_ 48(4):592–608.
  * International Maritime Organization (2011) International Maritime Organization (2011) Resolution MEPC.203(62) amendments to the annex of the protocol of 1997 to amend the international convention for the prevention of pollution from ships, 1973, as modified by the protocol of 1978 relating thereto (inclusion of regulations on energy efficiency for ships in MARPOL annex vi) .
  * Kellerer et al. (2004) Kellerer H, Pferschy U, Pisinger D (2004) _Knapsack problems._ (Springer, Berlin).
  * Konur et al. (2017) Konur D, Campbell J, Monfared S (2017) Economic and environmental considerations in a stochastic inventory control model with order splitting under different delivery schedules among suppliers. _Omega_ 71:46–65.
  * Lawson and Porteus (2000) Lawson D, Porteus E (2000) Multistage inventory management with expediting. _Operations Research_ 48(6):878–893.
  * Lemmens et al. (2019) Lemmens N, Gijsbrechts J, Boute R (2019) Synchromodality in the physical internet – dual sourcing and real-time switching between transport modes. _European Transport Research Review_ 1(10):518–526.
  * Lübbecke and Desrosiers (2005) Lübbecke M, Desrosiers J (2005) Selected topics in column generation. _Operations Research_ 53(6):1007–1023.
  * Nelsen (2006) Nelsen R (2006) _An Introduction to Copulas_ (Springer), second edition.
  * NTM (2015) NTM (2015) Methods and Manuals. URL https://www.transportmeasures.org/en/wiki/manuals/.
  * NTM (n.d.) NTM (n.d.) NTMCalc Basic 4.0. https://www.transportmeasures.org/ntmcalc/v4/basic/index.html#/, Accessed: 2021-7-20.
  * Palak et al. (2014) Palak G, Ekşioğlu S, Geunes J (2014) Analyzing the impacts of carbon regulatory mechanisms on supplier and mode selection decisions: An application to a biofuel supply chain. _International Journal of Production Economics_ 154:198–216.
  * Rosic̆ and Jammernegg (2013) Rosic̆ H, Jammernegg W (2013) The economic and environmental performance of dual sourcing: A newsvendor approach. _International Journal of Production Economics_ 143(1):109–119.
  * Sheopuri et al. (2010) Sheopuri A, Janakiraman G, Seshadri S (2010) New policies for the stochastic inventory control problem with two supply sources. _Operations Research_ 58(3):734–745.
  * Sun and Van Mieghem (2019) Sun J, Van Mieghem J (2019) Robust dual sourcing inventory management: Optimality of capped dual index policies and smoothing. _Manufacturing & Service Operations Management_ 21(4):912–931.
  * Svoboda et al. (2021) Svoboda J, Minner S, Yao M (2021) Typology and literature review on multiple supplier inventory control models. _European Journal of Operational Research_ 293(1):1–23.
  * Thomas and Tyworth (2006) Thomas D, Tyworth J (2006) Pooling lead-time risk by order splitting: A critical review. _Transportation Research Part E: Logistics and Transportation Review_ 42(4):245–257.
  * United Nations Comtrade Database (2020) United Nations Comtrade Database (2020) URL https://comtrade.un.org.
  * United Nations Conference on Trade and Development (2020) United Nations Conference on Trade and Development (2020) Review of maritime transport .
  * United Nations Environment Programme (2020) United Nations Environment Programme (2020) _Emissions Gap Report 2020_ (UNEP DTU Partnership).
  * Veeraraghavan and Scheller-Wolf (2008) Veeraraghavan S, Scheller-Wolf A (2008) Now or later: A simple policy for effective dual sourcing in capacitated systems. _Operations Research_ 56(4):850–864.
  * Vieweg et al. (2018) Vieweg M, Bongardt D, Hochfeld C, Jung A, Scherer E, Adib R, Guerra F (2018) _Towards Decarbonising Transport – A 2018 Stocktake on Sectoral Ambition in the G20_.
  * Whittemore and Saunders (1977) Whittemore A, Saunders S (1977) Optimal inventory under stochastic demand with two supply options. _SIAM Journal on Applied Mathematics_ 32(2):293–305.

##  7 Carbon accounting

In this section, we briefly explain how we determine the distribution
functions that we use for pseudo-random generation of the unit emissions in
our computational experiments. We utilize the United Nations Comtrade Database
(2020) to calculate the average unit weights for 122 groups of products
imported by The Netherlands in 2020. These product groups consist of two
categories: (i) apparel goods that are imported from Vietnam and (ii)
industrial goods that are imported from China and Germany. We consider air
transport (from Tan Son Nhat international airport) and sea transport (from
Haiphong port) as the fast and slow transport mode for the apparel category,
respectively. For the industrial goods category, we assume sea transport from
Shanghai, China, as the slow mode and road transportation from Stuttgart,
Germany, as the fast mode.

We rely on the Network for Transport Measures methodology (NTM 2015) to model
and measure transportation emissions based on equation 1. This model has been
widely used in literature (e.g., Hoen et al. 2014a, b). Following the NTM
methodology, we first compute the overall carbon emissions generated by a
single vehicle and then allocate a proportion of those emissions to each
freight unit carried by the vehicle.

Sea transportation. We assume that all products are shipped via container. The
average age of the container fleet worldwide is around 12 years and the
average vessel size (dwt) of container ships with age 10-14 is 43,993 ton
(United Nations Conference on Trade and Development 2020). Based on section 7
of the NTM framework and resolutions of the Marine Environment Protection
Committee (International Maritime Organization 2011), we approximate sea
transportation emissions in kilograms of $CO_{2}$ of one unit of a certain
product with weight $w$ (in kilogram) for a certain trip with distance $d$ (in
kilometers) using the following relation,

$e_{sea}=w\cdot EI_{ship}\cdot 10^{-3}\cdot d$

where $EI_{ship}$ is kilograms of $CO_{2}$ emissions per kilogram weight per
kilometer. Furthermore $EI_{ship}$ is computed through,

$EI_{ship}=\frac{(a\cdot dwt^{-c})/(PDR_{ship}\cdot LCU)}{1.852}$

where $a$ and $c$ are constants, $dwt$ is the deadweight tonnage of the ship,
$LCU$ is average load capacity utilization, $F(LCU)$ is fuel consumption as a
function of load, and $PDR_{ship}$ is the payload of the ship. $1.852$ is the
nautical mile to km conversion coefficient. For a container ship, NTM
methodology states: $a=0.17422,c=0.201,LCU=0.70,F(LCU)=1$, and $PDR=0.8$.
Succinctly, we have for the total emissions in kilograms of $CO_{2}$ of one
unit of a product with weight $w$ (in kilogram) for a sea trip with distance
$d$ (in kilometers)

$e_{sea}=w(1.996\cdot 10^{-5}\cdot d).$ (8)

Air transportation. Our calculations for the emissions of air transportation
are based on section 8 of the NTM Framework. We consider an Airbus A310-300 F
as the aircraft. Based on the May 2021 Air Cargo Market Analysis of The
International Air Transport Association, we assume an average international
cargo load factor of 65%. Following the NTM Framework, we have the following
relation for air transportation emissions

$e_{air}=\frac{w}{c_{max}}(CEF+VEF\cdot d),$

where $c_{max}$ is the maximum freight load, $CEF$ is the constant emissions
factor, and $VEF$ is the variable emissions factor. $CEF$ and $VEF$ are the
outcomes of applying a linear regression on real data provided by the NTM. We
obtain the $CEF$ and $VEF$ parameters via interpolation over the associated
tables provided by the NTM. We furthermore assume $c_{max}=39,000kg$ as per
section 8.3.1 and perform the interpolation on table 4.1 of section 8.2.1.
Succinctly, we have for the total emissions in kilograms of $CO_{2}$ of one
unit of a product with weight $w$ (in kilogram) for an air trip with distance
$d$ (in kilometers)

$e_{air}=w(1.525\cdot 10^{-1}+4.938\cdot 10^{-4}\cdot d).$ (9)

Road transportation. We rely on Hoen et al. (2014a) to obtain the emission
units of road transportation. They too rely on the NTM framework to estimate
$CO_{2}$ emissions from road transportation in Europe. In particular, they
approximate the total emissions in kilograms of $CO_{2}$ of one unit of a
product with weight $w$ (in kilogram) for a truck trip with distance $d$ (in
kilometers) as

$e_{road}=w(3.214\cdot 10^{-4}+4.836\cdot 10^{-5}\cdot d).$ (10)

Distances. We use NTMCalc Basic 4.0 (NTM n.d.), which is an online tool
provided by the NTM for approximating emissions, to calculate travel distances
between the origin destination pairs as described at the beginning of this
section. Based on this tool we find that the sea distance between Haiphong and
Rotterdam is 9,610 nautical miles (17,798 km) and the sea distance between
Shanghai and Rotterdam is 10,525 nautical miles (19,492 km). The distance
traveled by aircraft between Tan Son Nhat international airport and Rotterdam
The Hague Airport is 10,073 km, and the road distance between Stuttgart and
Rotterdam is 633 km. With these distances, we compute the total kilogram
$CO_{2}$ emissions of one unit of product with weight $w$ for each mode-trip
category using Equations (8)-(10). We call these emission coefficients and
they are presented in Table 6 below.

Table 6: Emissions coefficients for each mode-trip category. Industry | Slow Mode | Fast Mode
---|---|---
Apparel | $3.552\cdot 10^{-1}$ | $5.127$
Industrial | $3.891\cdot 10^{-1}$ | $3.093\cdot 10^{-2}$

Fitting distribution functions. We use the emission coefficients from Table 6
to calculate the unit emissions for the 122 groups of products mentioned at
the beginning of this section. We subsequently use maximum likelihood
estimation on the resulting emission units to find distribution functions from
which we can sample the emission units of the fast and slow transport modes
for all three assortment types in our computational experiment. The emission
unit distributions for assortment type 1 are based on the apparel category,
the emission unit distributions for assortment type 2 are based on the
industrial category, and the emission unit distributions for assortment type 3
are based on both categories. The final distribution functions are provided in
Table 3.

##  8 Generating correlated random numbers

Suppose $X$ and $Y$ are two real random variables with marginal distribution
functions $F$ and $G$, respectively. Suppose their joint distribution is bi-
variate standard normal $\mathcal{N}_{\rho}$ with Pearson’s correlation
coefficient $\rho=\mbox{\bf Cov}(X,Y)/(\sqrt{\mbox{\bf Var}[X]}\sqrt{\mbox{\bf
Var}[Y]})$ (Nelsen 2006). Let $Z$ be a vector of size two with independent
random elements that have standard normal distributions $\Phi$, and let $W=AZ$
be a linear combination of $Z$ with

$A=\begin{bmatrix}1&0\\\ \rho&\sqrt{1-\rho^{2}}\end{bmatrix}.$

It can be shown that $W$ has a bivariate normal distribution
$\mathcal{N}_{\rho}$ with covariance matrix $\Sigma=AA^{T}$ (see, e.g., Gut
2009). We use this result to sample from $X$ and $Y$ as follows:

  1. 1.

Generate the vector $Z=\begin{bmatrix}Z_{1}\\\ Z_{2}\end{bmatrix}$ by
independently sampling from a standard normal distribution function,

  2. 2.

Calculate the bivariate normal sample $W=\begin{bmatrix}W_{1}\\\
W_{2}\end{bmatrix}=AZ$,

  3. 3.

Generate the required samples by inversion $\begin{bmatrix}X\\\
Y\end{bmatrix}=\begin{bmatrix}F^{-1}(\Phi(W_{1}))\\\
G^{-1}(\Phi(W_{2}))\end{bmatrix}$.

##  9 Benchmark approaches

The mathematical formulation for the blanket mode selection approach, which
enforces emission constraints $\mathcal{E}_{j}^{max}$ per item $j\in J$, is
called Problem $(\mathit{BMS})$ and is given as follows:

$\displaystyle(\mathit{BMS})\qquad\qquad$ $\displaystyle\displaystyle\min$
$\displaystyle C(\mathbf{S}_{f},\mathbf{\Delta})$ subject to
$\displaystyle\mathnormal{E}_{j}(S_{j,f},\Delta_{j})\leq\mathcal{E}_{j}^{max},$
$\displaystyle\qquad\forall j\in J,$
$\displaystyle(S_{j,f},\Delta_{j})\in(\mathbb{R}\times\mathbb{R}_{0}),\quad$
$\displaystyle\qquad\forall j\in J.$

The mathematical formulation for the static mode selection approach, which
selects one transportation mode per item, is called Problem $(\mathit{SMS})$
and is given as follows:

$\displaystyle(\mathit{SMS})\qquad\qquad$ $\displaystyle\min$
$\displaystyle\sum_{j\in J}C_{j,s}(S^{*}_{j,s})x_{j,s}+\sum_{j\in
J}C_{j,f}(S^{*}_{j,f})x_{j,f}$ subject to $\displaystyle\sum_{j\in
J}E_{j,s}x_{j,s}+\sum_{j\in J}E_{j,f}x_{j,f}\leq\mathcal{E}^{max},$
$\displaystyle x_{j,s},x_{j,f}\in\\{0,1\\},$ $\displaystyle\forall j\in J,$

where
$C_{j,s}(S_{j,s}):=h_{j}\mathbb{E}\left[\left(S_{j,s}-\textstyle\sum_{t=0}^{l_{j,s}}D^{t}_{j}\right)^{+}\right]+p_{j}\mathbb{E}\left[\left(\textstyle\sum_{t=0}^{l_{j,s}}D^{t}_{j}-S_{j,s}\right)^{+}\right]+c_{j,s}\mathbb{E}[D_{j}]$,
and
$C_{j,f}(S_{j,f}):=h_{j}\mathbb{E}\left[\left(S_{j,f}-\textstyle\sum_{t=0}^{l_{j,f}}D^{t}_{j}\right)^{+}\right]+p_{j}\mathbb{E}\left[\left(\textstyle\sum_{t=0}^{l_{j,f}}D^{t}_{j}-S_{j,f}\right)^{+}\right]+c_{j,f}\mathbb{E}[D_{j}]$
are the average cost-rate for using the slow and fast transport mode
respectively, $S^{*}_{j,s}:=\operatorname{argmin}_{S_{j,s}}C_{j,s}(S_{j,s})$,
and $S^{*}_{j,f}:=\operatorname{argmin}_{S_{j,f}}C_{j,f}(S_{j,f})$ are the
base-stock levels that minimize inventory related costs when exclusively using
the slow and fast transport mode respectively, and
$E_{j,s}:=e_{j,s}\mathbb{E}[D_{j}]$ and $E_{j,f}:=e_{j,f}\mathbb{E}[D_{j}]$
are the emissions per time unit of shipping exclusively with the slow and fast
transport mode respectively.