Omega 127 (2024) 103082

Contents lists available at ScienceDirect

Omega
journal homepage: www.elsevier.com/locate/omega

Robust inventory routing problem under uncertain demand and risk-averse

criterion✩
Yuqiang Feng, Ada Che ∗, Na Tian
School of Management, Northwestern Polytechnical University, Xi’an, China

ARTICLE INFO ABSTRACT

Keywords: Inventory routing problem (IRP), which plays an important role in implementing vendor-managed-inventory
Inventory routing problem strategy, is to determine the optimal timing and quantity of products to be delivered as well as the
Uncertain demand optimal vehicle delivery routes from the vendor to retailers. Retailers’ demand is usually uncertain and their
Ambiguous distribution
distributions are ambiguous. We utilize limited historical demand data of retailers to construct a series of
Risk-averse criterion
scenarios to characterize demand uncertainties. Moreover, we develop a framework to construct ambiguity
Distributionally robust optimization
sets and describe the ambiguities in demand distributions. To cater to the new feature of IRP as being exposed
to downside risk and balance the mean and risk level of cost, we incorporate the worst-case mean-conditional
value-at-risk (M-CVaR) criterion into the objective function and present a new distributionally robust IRP
model. We transform the proposed model into a tractable formulation based on duality theory. A case study
demonstrates the feasibility of the proposed method. The results indicate that our model can provide a robust
delivery solution for immunizing against the influence of ambiguous demand distributions. In addition, an
out-of-sample performance analysis shows that the delivery solution provided by our model can effectively
reduce the product shortage rate.

1. Introduction costs of vendors and retailers’ inventory costs and improves the supply
chain’s overall efficiency.
The vendor managed inventory (VMI) strategy is of great impor- The inventory routing problem (IRP) proposed by Bell et al. [7]
tance in the current business environment. In a traditional business en- has received considerable attention since the implementation of the
vironment, information in supply chains is usually fragmented, which VMI strategy. IRP integrates the vehicle routing problem and inven-
means there is a lack of effective communication and coordination tory management problem and optimizes the efficiency of the supply
between different links. This often results in vendors being unable chain on the basis of vendors (upstream) and retailers’ (downstream)
to adjust their policies in a timely manner to satisfy customer re- information-sharing and establishing cooperation. VMI is an agreement
quirements [1]. An obvious characteristic of modern supply chains is that promotes cooperation among different links in the supply chain,
their high complexity. Collaboration and information-sharing among and IRP is the practice that executes this agreement. In numerous
various links in the supply chain are emphasized [2]. In this context, cases, Archetti & Speranza [8] identified that the delivery solution
the VMI management strategy is proposed. This strategy advocates from vendors to retailers developed by the IRP reduced the system cost
cooperation and information-sharing between vendors and retailers in by an average of 23.97% and a maximum of 33.63% compared with
real-time [3,4]. In this strategy, vendors behave as central decision- the delivery solution developed by the traditional approach, that is,
makers in determining timing, quantity, and routing of deliveries to
optimizing the inventory and routing problems separately. We present
retailers to minimize the total cost [5]. Using the VMI strategy, ven-
two real cases where a business firm benefited from IRP. The first case
dors can usually obtain more uniform utilization of production and
is about Midea, which is a famous home appliance company in China.
transportation resources. At the same time, retailers could focus less on
Before implementing VMI and IRP, Midea often held high inventory
inventory management activities and more on their core businesses [6].
levels of components and needed to hold at least 5–7 days of inventory.
In general, the VMI strategy reduces the production and transportation
After implementing these two policies, it was reduced to nearly 3

✩ Area: Production Management, Scheduling and Logistics. This manuscript was processed by Associate Editor Li.
∗ Corresponding author.
E-mail addresses: fengyuqiang2021@126.com (Y. Feng), ache@nwpu.edu.cn (A. Che), tianna@nwpu.edu.cn (N. Tian).
1
https://zhuanlan.zhihu.com/p/402936642
2
https://bbs.kuguanyi.com/thread-1631-1-1.html

https://doi.org/10.1016/j.omega.2024.103082
Received 7 September 2023; Accepted 23 March 2024
Available online 26 March 2024
0305-0483/© 2024 Elsevier Ltd. All rights reserved.
Y. Feng et al. Omega 127 (2024) 103082

days.1 Another case focuses on the collaboration between Nestle and value based on a given probability level. More importantly, compared
Carrefour. Before implementing VMI and IRP, the product arrival rates to VaR, CVaR is now widely accepted by decision-makers because it
of Nestle to Carrefour and Carrefour to retail supermarkets are 80% possesses the properties of a coherent risk measure, that is, subaddi-
and 70%, respectively. After implementing these two policies, these tivity, positive homogeneity, monotonicity, and positive homogeneity.
two rates were raised to 95% and 90%.2 In addition, other prominent Moreover, the CVaR part can be transformed into linear form under
companies have implemented the IRP management mode, such as P&G a discrete probability distribution. The traditional CVaR, based on an
and Walmart [9], Amazon and IKEA [10]. Moreover, some IRP variants, exact probability distribution, has been applied to the risk measure-
such as Location IRP [11], maritime IRP [12], IRP with perishable ment of various operations management problems. However, for the
products [13], have recently been generated. These cases and related IRP problem in this study, the probability distribution information
studies demonstrate the pivotal functions of IRP in the supply chain. of retailers’ demand is not precisely known, which results in CVaR
In the current business environment, demand is characterized by un- not being applied directly. In this context, under a given downside
certainty because it is often influenced by various factors, e.g., dynamic risk, we incorporate the worst-case CVaR measure based on a set of
product prices, product diversity, and product promotions. Uncertain discrete probability distributions into the decision-making process of
demand exacerbates the bullwhip effect of the supply chain, which IRP. Moreover, considering that a risk-averse IRP decision-maker might
further results in uneven capacity utilization, more inventory build-up, prefer a delivery solution that balances expectation cost and risk, we
shortages, delayed deliveries, and confusion in production planning. In introduce the balance parameter and develop the worst-case M-CVaR
this context, neglecting demand uncertainty can lead to large deviations criterion in the objective function.
in the costs and planned costs (i.e., budget) of the enterprise. Tradition- Stimulated by the aforementioned observations, this study con-
ally, uncertain demands are assumed to be random or fuzzy variables, structs two different ambiguity sets to characterize the ambiguities
and their exact distributions are typically assumed to be known. Under of demands distributions of retailers, utilizes the worst-case M-CVaR
this assumption, stochastic optimization (SO) and fuzzy optimization criterion to balance the mean and risk level of cost, and proposes a new
(FO) methods can address uncertain demand. However, in most cases, distributionally robust M-CVaR IRP model to seek an optimal robust
only partial information on the distribution of uncertain demand can be delivery solution. The contributions of this study are summarized as
obtained, and the precise distribution is likely to be challenging to es- follows.
timate [14]. Therefore, it is strict to assume the distribution of demand
(1) This study develops a framework to describe the uncertain de-
is precise for the above two optimization methods to some extent. In
mands of retailers with inexact distributional information. We do
this context, how to optimize the IRP under uncertain demand with
not generate demand scenarios based on the given distribution
ambiguous distribution has become a problem that must be addressed
or assume that the distributions of demands are exactly known
urgently.
as in the existing literature. Based on limited historical demand
This is a promising strategy for applying the distributionally robust
data of retailers, we construct a sequence of demand scenarios to
optimization (DRO) method to address uncertain demand with ambigu-
depict uncertain demands and two ambiguity sets to characterize
ous distribution. The DRO method originates from the SO method but
the ambiguous distributions.
does not require the exact distribution of uncertain parameter to be
(2) This study proposes a novel distributionally robust IRP model.
known, and provides an appropriate optimization framework for solv-
We employ the worst-case M-CVaR criterion over the constructed
ing the problem with partial distribution information of the uncertain
ambiguity sets to balance the mean and risk level of retailers’ in-
parameter [15,16]. The DRO method uses partially known distribu-
ventory and shortage penalty costs for the risk-averse manager.
tional information to construct an ambiguity set and requires that the
To the best of our knowledge, it is the first time that such a
distribution of the uncertain parameter belongs to the constructed set.
hybrid modeling technique has been used that does not require
Based on the worst-case distribution of uncertain parameter over the
exactly known distributional information in the literature related
constructed ambiguity set, the optimization problem with uncertain
to IRP.
parameter is reformulated into a deterministic optimization problem
(3) This study reformulates the proposed distributionally robust IRP
via a series of techniques, such as the duality theory. This study uses
model into tractable forms using duality theory. We validate the
limited historical demand data to construct ambiguity set, and the DRO
following two critical points through numerical experiments. If
method is used to address the ambiguity in demand distribution in IRP.
ambiguity sets contain only nominal distribution, the optimal so-
One new feature of IRP in the current business environment is its
lution generated by the presented distributionally robust model
exposure to downside risk. Herein, the downside risk is defined as the
is the same as that generated by the stochastic model. If ambigu-
probability of the cost exceeding a threshold value [17], has become
ity sets contain more possible distributions, the performance of
an essential guide for making various decisions, including IRP deci-
the delivery solution provided by the distributionally robust IRP
sions [18]. In terms of IRP, the decision-maker usually needs to design a
model is superior to that provided by the stochastic IRP model.
delivery solution for the planning period. There is no doubt that the cost
(i.e., budget) of the delivery solution must be prepared in advance. This The rest of this study is organized as follows. In Section 2, we review
budget is tied up in the planning period so that the enterprise might and analyze the relevant literature from the perspective of different
miss other investment opportunities. Therefore, the IRP decision-maker types of distributional information of uncertain demands. Section 3 first
urgently needs to know the budget and some information that cost describes the details of the studied problem and then establishes a new
over the budget under a given downside risk. Consequently, the vendor distributionally robust M-CVaR IRP model with uncertain demands. In
needs incorporate an appropriate downside risk measure into the IRP Section 4, we use limited historical demand data to construct ambiguity
decision-making process. Classical semi-deviation and square-deviation sets. In Section 5, we conduct model analysis and derive the tractable
measure the deviation between cost and budget, neglecting the down- form of the developed model. We conduct a case study to demonstrate
side risk. Additionally, a model with a square-deviation is usually a the feasibility of the proposed model in Section 6. We summarize some
quadratic programming model, which is not easily solved. Value-at- management insights in Section 7. Conclusion is presented in Section 8.
risk (VaR) and CVaR are well-known downside risk measures [19]. VaR
only estimates the maximum possible cost at a given probability level 2. Literature review
(downside risk parameter) for a planning horizon. However, it fails
to capture the situation where the cost exceeds the VaR value. CVaR, This section reviews the literature that is related to our work and
which was proposed by Rockafellar & Uryasev [20], not only provides shows the existing research gaps on IRP to highlight our contributions.
the VaR value but also estimates the expected cost that exceeds VaR We review the related literature from the following three aspects.

2
Y. Feng et al. Omega 127 (2024) 103082

2.1. IRP with known distributional information of demands 2.2. IRP with unknown distributional information of demands

Under the condition that the distributional information is precisely The robust optimization (RO) method can be utilized to deal with
known, some researchers have used the SO method to address the demands uncertainties in the IRP under the condition that the dis-
uncertain demands in IRP. Chen & Lin [21], Huang & Lin [22], Yu tributional information is unknown. Solyalı et al. [47] studied an
et al. [23], Abdul Rahim et al. [24], Soysal et al. [25], Micheli & IRP in which a vendor distributes a single product to retailers facing
Mantella [26], and Kumar et al. [27] described the uncertain demands dynamic uncertain demands within a finite discrete time horizon. They
using gaussian distribution and proposed corresponding stochastic IRP constructed a polyhedral uncertainty set to characterize the uncertain
models, where the mean and standard deviation are given. Based on
demands, proposed robust and adjustable robust mixed-integer linear
risk aversion, Chen & Lin [21] investigated a hedge-based stochastic
programming (MILP) models for IRP, and designed a branch-and-cut
IRP system to hedge the demand volatility with maximal net present
algorithm to solve the model. Jafarkhan & Yaghoubi [48] and Liu
value (NPV). Kleywegt et al. [28] developed a stochastic IRP model
et al. [49] set a scenario set for the uncertain demands and applied the
based on a Markov decision process and introduced an approximation
scenario-based RO method proposed by Mulvey & Ruszczyński [50] to
method to solve the proposed model. In studies by Hvattum et al. [29],
Bertazzi et al. [30], and Bertazzi et al. [31], demands were specified get an IRP solution that is progressively less sensitive to realizations of
through general discrete distributions. Hvattum et al. [29] developed demands from the given scenario set. Gholami-Zanjani et al. [51] stud-
heuristic methods based on finite-scenario trees to approximate the ied a robust IRP with uncertain demands and transportation costs and
underlying stochastic IRP. Juan et al. [32] employed a multi-start transformed the proposed robust model into three different MILP forms
meta-heuristic with monte carlo simulation and proposed a hybrid based on polyhedral, box, and polyhedral-box uncertainty sets. Fardi
simulation–optimization method to address the single-period IRP with et al. [52] presented a multi-depot, multi-vehicle robust IRP model
stochastic demands. Gruler et al. [33] simulated customer demands based on a budget uncertainty set, where products from different depots
from a skewed geometric distribution to incorporate demand uncer- were delivered to different retailers. Golsefidi & Jokar [53] explored
tainty and designed a meta-heuristic with a simulation to solve the the production–inventory-routing problem with a reverse product flow.
single-period stochastic IRP model. Based on given discrete and normal They simultaneously considered multiple uncertain conditions: retail-
demand distributions of medical drug, Nikzad et al. [34] formulated a ers’ demands, the quantity of defective products in each period, and
two-stage stochastic IRP model, a stochastic IRP model with continuous reproduction cost. They used the RO method to address demands uncer-
probabilistic constraints, and a stochastic IRP model with discrete tainties based on a budget uncertainty set. Shang et al. [11] constructed
probabilistic constraints. They designed a two-phase meta-heuristic a budget uncertainty set to characterize the demands uncertainties
algorithm to solve the presented models. In terms of the IRP with and proposed a robust location-IRP model. By applying the conversion
perishable products, in the study by Onggo et al. [35], retailers’ de- approach proposed by Bertsimas & Sim [54], the developed robust
mand in each period was assumed to be uncertain, and log-normal model was transformed into a MILP formulation.
probability distributions were used to characterize these uncertainties. The RO method assumes that demand falls into a predefined un-
Furthermore, Soysal [36], Alvarez et al. [37], and Mousavi, Bashiri,
certainty set, does not consider distributional information and usually
& Nikzad [38] assumed that the demands of retailers change under
yields over-conservative solutions [55]. This motivates us to fully use
different scenarios and that the probabilities of occurrences for these
the distributional information of retailers in the optimization process
scenarios are known in advance.
of IRP.
Based on the condition that the distributional information of de-
mands in IRP is exactly known, except for the SO method, many
scholars also used the FO method to address the demands uncertainties. 2.3. IRP with partially known distributional information of demands
Niakan & Rahimi [39] studied a health-care IRP and assumed the
demands of hospitals for drug in each period are triangular fuzzy Under the condition that distributional information is partially
variables. Rahimi et al. [40] considered that the uncertain demands are known, few scholars have used the DRO method to deal with the
triangular fuzzy variables and developed a multi-objective IRP model demands uncertainties in the IRP. Following the idea proposed by
from perspectives of profit, service level, and green criteria. Rahimi Delage & Ye [56], Liu et al. [10] used the estimation errors of the
et al. [41] assumed that uncertain demands follow trapezoidal possi- mean and covariance matrix of demands to develop an ambiguity set,
bility distributions and proposed a possibilistic chance constrained pro- and proposed a distributionally robust IRP model to maximize the joint
gramming model. They converted the proposed model into a tractable probability of no stockout at the end of each period. Li et al. [57]
form and designed a branch-and-cut algorithm to solve it. Hu et al. [42] assumed that demands vary under different scenarios, used support
considered that the uncertain demands are triangular fuzzy variables, and mean information to develop an ambiguity set scenario-based, and
proposed a credibilistic goal programming model for the IRP with proposed a distributionally robust IRP model to minimize the expec-
hazardous products. Based on fuzzy demands, Ghasemkhani et al. [43] tation of cost. Cui et al. [58] utilized scenario-based support, mean,
investigated an IRP with multiple perishable products, multi-period, mean absolute deviation, and correlation information to characterize
time windows, and heterogeneous fleets to maximize the total profit. the ambiguity of demand distribution. They first proposed a service
Avila-Torres & Arratia-Martinez [44] used triangular fuzzy numbers
violation index (SVI) to estimate the risk of service violations at the
to portray uncertain demands, constructed a fuzzy IRP model, and
inventory level and then developed a distributionally robust IRP model
applied the 𝛼-cut method proposed by Adamo [45] to transform the
to minimize the SVI.
constraints associated with uncertain demands into deterministic con-
It is worth noting that these three studies do not balance the mean
straints. Ghasemkhani et al. [46] assumed the uncertain demands are
and risk level based on the ambiguous distributions of demands. In this
trapezoidal and triangular fuzzy variables, respectively, and trans-
formed the developed credibility-based chance-constrained IRP model study, we develop a worst-case M-CVaR criterion to balance the mean
into a deterministic model based on the credibility theory. and risk level of retailers’ inventory and shortage penalty costs although
Suppose the precise distributional information of demands is easily the DRO method is also used to address ambiguities.
accessible. In that case, the decision-maker can use the SO or FO To pinpoint the research gaps within the existing studies and show-
method to address demands uncertainties in IRP at hand. Nevertheless, case the innovations of this work on IRP, we outline the research
it is often challenging for the IRP decision-maker to obtain the precise features of the relevant literature in Table 1. We identify three research
distributional information of demands of retailers. This encourages us gaps by observing Table 1. First, there is limited literature on optimiz-
not to make strict assumptions about the exact distribution and to adopt ing the IRP by employing the DRO method under the condition that
the DRO method to handle the demands uncertainties. the distributions of uncertain demands are ambiguous. Second, no prior

3
Y. Feng et al. Omega 127 (2024) 103082

Table 1
The comparison of the relevant literature and this study.
Reference DID OM Risk aversion Developing framework Objective
to handle uncertainty
Known Unknown Partially SO FO RO DRO Downside risk
known
√ √
Kleywegt et al. [28] – – – – – – – Minimize total cost
√ √
Chen & Lin [21] – – – – – – – Maximize NPV
√ √
Hvattum et al. [29] – – – – – – – Minimize total cost
√ √
Huang & Lin [22] – – – – – – – Minimize total cost
√ √
Yu et al. [23] – – – – – – – Minimize total cost
√ √
Solyalı et al. [47] – – – – – – – Minimize total cost
√ √
Abdul Rahim et al. [24] – – – – – – – Minimize total cost
√ √
Juan et al. [32] – – – – – – – Minimize expected cost
√ √
Niakan & Rahimi [39] – – – – – – – Minimize total cost
√ √
Bertazzi et al. [31] – – – – – – – Minimize expected cost
√ √
Soysal et al. [25] – – – – – – – Minimize total cost
√ √
Soysal [36] – – – – – – – Minimize total cost
√ √
Rahimi et al. [40] – – – – – – – Multi-objective
√ √
Kazemi et al. [41] – – – – – – – Minimize total cost
√ √
Jafarkhan & Yaghoubi [48] – – – – – – – Minimize total cost
√ √
Micheli & Mantella [26] – – – – – – – Minimize total cost
√ √
Gruler et al. [33] – – – – – – – Minimize expected cost
√ √
Hu et al. [42] – – – – – – – Minimize total cost
√ √
Fardi et al. [52] – – – – – – – Minimize total cost
√ √
Gholami-Zanjani et al. [51] – – – – – – – Minimize total cost
√ √
Liu et al. [10] – – – – – – – Maximize service level
√ √
Nikzad et al. [34] – – – – – – – Minimize total cost
√ √
Ghasemkhani et al. [43] – – – – – – – Maximize total profit
√ √
Onggo et al. [35] – – – – – – – Minimize expected cost
√ √
Golsefidi & Jokar [53] – – – – – – – Minimize total cost
√ √
Liu et al. [49] – – – – – – – Minimize total cost
√ √
Alvarez et al. [37] – – – – – – – Minimize total cost
√ √
Avila-Torres & Arratia-Martinez [44] – – – – – – – Minimize total cost
√ √
Ghasemkhani et al. [46] – – – – – – – Minimize total cost
√ √
Kumar et al. [27] – – – – – – – Minimize total cost
√ √
Mousavi, Bashiri, & Nikzad [38] – – – – – – – Minimize total cost
√ √
Shang et al. [11] – – – – – – – Minimize total cost
√ √
Cui et al. [58] – – – – – – – Minimize service violation index
√ √
Li et al. [57] – – – – – – – Minimize total cost
√ √ √ √
This research Minimize worst-case M-CVaR

DID and OM denote the distribution information of demands and optimization method, respectively.

studies have considered the downside risk in the decision-making pro- types of decisions and minimize total cost. The first determines which
cess of IRP under the risk-averse criterion. Third, no relevant literature retailers should be visited during each period. The second determines
has developed a reasonable framework to make it more convenient for the optimal delivery quantity for retailers in each period. The third
the IRP decision-maker to address uncertain demands. involves designing the optimal vehicle routing decisions in each period,
To address these gaps, this study designs a framework to deal i.e., determining the optimal number of vehicles to complete delivery
with uncertain demands based on limited historical demand data and and their optimal delivery routes. There are 𝐾 homogeneous vehicles
constructs a distributionally robust M-CVaR IRP model with a downside in the distribution center. Delivery vehicles start at the distribution
risk measure to balance the mean and risk level of cost. The constructed center and deliver products to retailers. After delivery, the vehicles
model can flexibly degenerate to one in different decision-making are returned to the distribution center. Retailers are not involved in
perspective by appropriately setting the inputs. In particular, when the decision-making. In each period, the delivery quantity cannot be split
balance parameters are equal to 0, the upper bound of the CVaR value if the vendor delivers the products to the retailer. If a retailer has a
of the total cost can be obtained for a given downside risk. shortage in the 𝑡th period and if the vendor delivers products to it in
the 𝑡 + 1th period, then the products will take priority in filling the
3. Distributionally robust M-CVaR IRP model shortage. Retailers’ demand is uncertain and belongs to the scenario
set where the probabilities of demand scenarios are ambiguous.
In this section, we first describe the studied IRP and then develop We can define a complete graph 𝒢 = {𝒱 , 𝒜 } to represent the
a new distributionally robust M-CVaR IRP model with ambiguous above IRP, where 𝒱 = {0, … , 𝑁, 𝑁 + 1} is the set of nodes, and
distributions of uncertain demands. 𝒜 = {(𝑖, 𝑗) ∶ 𝑖, 𝑗 ∈ 𝑉 , 𝑖 < 𝑗} is the set of feasible edges. The nodes
involve a vendor and 𝑁 retailers. Node 0 and node 𝑁 + 1 represent the
3.1. Problem description vendor. The edges refer to all possible paths between any two nodes.
In Fig. 1, an example of a vendor supplying 7 retailers is utilized to
This study investigates the IRP, which is a combination of inventory display the structure of the studied IRP. From Fig. 1, we observe that in
management and vehicle routing problems. The system consists of a the (𝑡 − 1)th period, vendor delivers the products with 3 vehicles, while
vendor (distribution center) and 𝑁 geographically dispersed retailers. in the 𝑡th and (𝑡 + 1)th periods, vendor just uses 2 vehicles for delivery.
The vendor is responsible for delivering the products to the retailers. In addition, the vehicle does not visit the 2nd and 6th retailers in the
The demands of these retailers are independent of each other. However, 𝑡th period, and the 5th and 7th retailers in the (𝑡 + 1)th period.
the vendor does not know the quantity of products that should be In addition to the above statements, this study considers the fol-
delivered to each retailer. The planning horizon length is 𝑇 . The lowing assumptions: (i) the vendor has enough products to satisfy the
vendor (the central decision-maker) must simultaneously make three demands of the retailers; (ii) the inventory and shortage penalty costs

4
Y. Feng et al. Omega 127 (2024) 103082

Fig. 1. Illustration of the studied IRP.

of the IRP system only arise from the warehouse at the retailers (node 𝑖, 3.2.1. Constructing objective function
𝑖 ∈ 𝒩 ); (iii) the residual value of the remaining products after the end In our model, the objective function comprises two components. The
of the sales period is not considered; and (iv) during each period, the first component is built as follows:
number of deliveries made by the vendor to each retailer is less than ∑ ∑ ∑
or equal to 1. These assumptions were also made by Cui et al. [58], TC1 (𝒙, 𝒒, 𝒛, 𝒚) = 𝑐𝑖,𝑗 𝑥𝑡,𝑖,𝑗 + 𝑐𝑖ℎ I0,𝑖 ,
𝑡∈𝒯 (𝑖,𝑗)∈𝒜 𝑖∈𝒩
Gruler et al. [33], Micheli & Mantella [26], and Soysal et al. [25].
The required sets, parameters, and variables used are listed as where 𝒙 = (𝑥𝑡,𝑖,𝑗 )𝑡∈𝒯 ,(𝑖,𝑗)∈𝒜 , 𝒒 = (𝑞𝑡,𝑖 )𝑡∈𝒯 ,𝑖∈𝒩 , 𝒛 = (𝑧𝑡,𝑖,𝑗 )𝑡∈𝒯 ,(𝑖,𝑗)∈𝒜 , 𝒚 =
follows: (𝑦𝑡,𝑖 )𝑡∈𝒯 ,𝑖∈𝒩 . The first part of TC1 (𝒙, 𝒒, 𝒛, 𝒚) is the transportation cost of
Sets: the IRP system. The second part is the initial inventory cost of the IRP
𝒩 : Set of all retailers, indice 𝑖, 𝑗 ∈ 𝒩 = {1, 2, … , 𝑁}; system.
𝒱 : Set of all nodes, indice 𝑖, 𝑗 ∈ 𝒱 ; The second component is the worst-case M-CVaR value of the
𝒯 : Set of periods, indice 𝑡, 𝑡′ ∈ 𝒯 = {1, 2, … , 𝑇 }; inventory and shortage penalty costs in the IRP system. For node 𝑖,
𝒮 : Set of demand scenarios, indice 𝑠 ∈ 𝒮 = {1, 2, … , 𝑆}; 𝑖 ∈ 𝒩 , in each period and in the 𝑠th demand scenario, the inventory
Parameters: cost and shortage penalty cost correspond to the cost associated with
excess inventory (i.e., 𝛼𝑡,𝑖 𝑠 > 0) and unfilled demand (i.e., 𝛿 𝑠 > 0),
𝑐𝑖,𝑗 : The transportation cost between node 𝑖 and node 𝑗; 𝑡,𝑖
d𝑠𝑡,𝑖 : The demand of the 𝑖th retailer in the 𝑡th period of the 𝑠th respectively. The sum of the inventory and shortage penalty costs of
node 𝑖 in the 𝑡th period can be expressed as 𝑐𝑖ℎ 𝛼𝑡,𝑖 𝑠 + 𝑐 𝑝 𝛿 𝑠 . In the 𝑠th
demand scenario; 𝑖 𝑡,𝑖
P𝑠𝑖 : The probability of occurrence for the 𝑠th demand scenario of demand scenario, the total inventory and shortage penalty costs can be
the 𝑖th retailer; expressed as:
I0,𝑖 : The initial inventory level of the 𝑖th retailer; ∑ ∑( )
TC𝑠2 (𝜶 𝑠 , 𝜹𝑠 ) = 𝑐𝑖ℎ 𝛼𝑡,𝑖
𝑠
+ 𝑐𝑖𝑝 𝛿𝑡,𝑖
𝑠
,
𝐶𝑖 : The maximum inventory level of the 𝑖th retailer; 𝑖∈𝒩 𝑡∈𝒯
𝑄 : The maximum load of vehicles; 𝑠 ) 𝑠 𝑠
where 𝜶 𝑠 = (𝛼𝑡,𝑖 𝑡∈𝒯 ,𝑖∈𝒩 and 𝜹 = (𝛿𝑡,𝑖 )𝑡∈𝒯 ,𝑖∈𝒩 .
𝑐𝑖ℎ : Holding cost unit product of the 𝑖th retailer;
We introduce the auxiliary variable F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) to denote the sum of
𝑐𝑖𝑝 : Penalty cost unit product shortage of the 𝑖th retailer.
the inventory and shortage( penalty costs) of node 𝑖 in demand scenario
Decision variables: ∑ 𝑠 + 𝑐 𝑝 𝛿 𝑠 , where 𝜶 𝑠 = (𝛼 𝑠 )
𝑠, i.e., F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) = 𝑡∈𝒯 𝑐𝑖ℎ 𝛼𝑡,𝑖 𝑖 𝑡,𝑖 𝑖
𝑠
𝑡,𝑖 𝑡∈𝒯 , 𝜹𝑖 =
𝑥𝑡,𝑖,𝑗 : Binary variables. In the 𝑡th period, if a vehicle travels (𝑖, 𝑗) in 𝑠 )
(𝛿𝑡,𝑖 𝑡∈𝒯 . The total inventory and shortage penalty costs in the 𝑠th
any direction, it is 1; otherwise, it is 0; ∑
demand scenario TC𝑠2 (𝜶 𝑠 , 𝜹𝑠 ) can be expressed as 𝑖∈𝒩 F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ).
𝑦𝑡,𝑖 : Binary variables. In the 𝑡th period, if the 𝑖th retailer is visited
In reality, it is a difficult task to obtain the precise probability
by any vehicle, it is 1; otherwise, it is 0;
distribution of demand for each retailer. In other words, the distribution
𝑧𝑡,𝑖,𝑗 : Non-negative continuous variables. In the 𝑡th period, 𝑧𝑡,𝑖,𝑗 and
is itself a source of uncertainty. Therefore, it is reasonable to assume
𝑧𝑡,𝑗,𝑖 represent the load and the residual capacity of the vehicle traveling
that the ambiguous discrete probability distribution 𝐏𝑖 = (P𝑠𝑖 )𝑠∈𝒮 of
from the 𝑖th retailer to the 𝑗th retailer, respectively.
demand for the 𝑖th retailer resides in an ambiguity set 𝒫𝑖 which we
𝑞𝑡,𝑖 : Non-negative continuous variables. The delivery quantity of
construct in the next Section, that is, 𝐏𝑖 ∈ 𝒫𝑖 . In this case, to balance
products from the vendor to the 𝑖th retailer in the 𝑡th period;
the mean and risk level of the sum of inventory and shortage penalty
I𝑠𝑡,𝑖 : The inventory level of the 𝑖th retailer in the end of the 𝑡th period
costs for node 𝑖, 𝑖 ∈ 𝒩 , we develop the following worst-case M-CVaR
in demand scenario 𝑠;
𝑠 : Non-negative continuous variables. The quantity of surplus
criterion:
𝛼𝑡,𝑖 { [ ] [ ]}
products of the 𝑖th retailer in the end of the 𝑡th period in demand max 𝛽𝑖 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 )𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) ,
𝐏𝑖 ∈𝒫𝑖
scenario 𝑠;
𝑠 : Non-negative continuous variables. The quantity of products
𝛿𝑡,𝑖 where 𝛽𝑖 denotes the balance parameter at node 𝑖. 1−𝜅𝑖 is the downside
shortage of the 𝑖th retailer in the end of the 𝑡th period in demand risk and reflects the probability of F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) exceeds the threshold value
[ ]
scenario 𝑠. (i.e., 𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) ) of node 𝑖. In the demand scenario 𝑠, the
realization value of the random variable F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) is F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ).
3.2. Model formulation In summary, we formulate the following objective function:
[ ]
∑ { [ ] [ ]}
TC1 (𝒙, 𝒒, 𝒛, 𝒚) + max 𝛽𝑖 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 )𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) .
In this subsection, we formulate the objective function and con- 𝑖∈𝒩
𝐏𝑖 ∈𝒫𝑖

straints of the proposed distributionally robust model.

5
Y. Feng et al. Omega 127 (2024) 103082

Remark 1. Herein, we explain why the CVaR is suitable. Two signifi- The following constraint (4) means the cumulative sum of the flows
cant risks are associated with IRP in the current business environment: originating from a node equals the maximum load of vehicle minus the
downside risk and shortage risk. The latter can be reduced by increasing delivery quantity absorbed by this node.
shortage penalty cost 𝑐𝑖𝑝 , ∀𝑖 ∈ 𝒩 . However, regardless of how the value ∑
𝑧𝑡,𝑖,𝑗 = 𝑄𝑦𝑡,𝑖 − 𝑞𝑡,𝑖 , ∀ 𝑡 ∈ 𝒯 , 𝑖 ∈ 𝒩 . (4)
of 𝑐𝑖𝑝 is adjusted, estimating the possible cost under a given downside
𝑗∈𝒱 ,𝑖≠𝑗
risk is impossible. VaR and CVaR are well-known downside risk mea-
sures and can be utilized to accomplish this task. Compared with VaR, The total quantity starting from the distribution center equals the
CVaR has the following advantages: ① CVaR satisfies subadditivity. total delivery quantities for retailers in the 𝑡th period. There are no
The CVaR value of the inventory and shortage penalty costs for all products arrive at the distribution center when routes are terminated
nodes is lower than the sum of the CVaR value of each node. This in the 𝑡th period. These two constraint are expressed as
finding is consistent with the theory of finance risk diversification. If we ∑ ∑
𝑧𝑡,0,𝑖 = 𝑞𝑡,𝑖 , ∀𝑡 ∈ 𝒯 , (5)
use VaR, we cannot obtain the corresponding conclusion because VaR 𝑖∈𝒩 𝑖∈𝒩
does not satisfy subadditivity. ② The CVaR part can be converted to a ∑
𝑧𝑡,𝑖,𝑁+1 = 0, ∀𝑡 ∈ 𝒯 . (6)
linear form under discrete demand probability distribution. However, 𝑖∈𝒩
VaR may include multiple local extremes with a discrete probability
Constraint (7) limits the range of decision variables 𝒛.
distribution [59]. ③ CVaR can provide the mean value of the cost
over the VaR value besides the VaR value of the IRP system. VaR only 0 ≤ 𝑧𝑡,𝑖,𝑗 ≤ 𝑄, ∀ 𝑡 ∈ 𝒯 , 𝑖, 𝑗 ∈ 𝒱 , 𝑖 ≠ 𝑗. (7)
provides the maximum cost that the IRP system may suffer for a given
downside risk and does not capture situations in which the cost exceeds Constraint (8) represents the sum of inventory and shortage penalty
the VaR value. This situation (albeit a small probability event) is also an costs of the 𝑖th node in the 𝑠th demand scenario.
important reference for the IRP decision-maker when making delivery ∑( )
F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) = 𝑐𝑖ℎ 𝛼𝑡,𝑖
𝑠
+ 𝑐𝑖𝑝 𝛿𝑡,𝑖
𝑠
, ∀𝑠 ∈ 𝒮 , 𝑖 ∈ 𝒩 . (8)
solutions. Therefore, it is suitable to incorporate the CVaR downside 𝑡∈𝒯
risk measure into the IRP decision-making process. 𝑠 and 𝛿 𝑠 .
Constraint (9) is established to compute 𝛼𝑡,𝑖 𝑡,𝑖
{
Remark 2. Here, we explain the motivations for incorporating the d𝑠1,𝑖 + 𝛼1,𝑖
𝑠 = I + 𝑞 + 𝛿𝑠 , ∀𝑠 ∈ 𝒮 , 𝑖 ∈ 𝒩 , (a)
0,𝑖 1,𝑖 1,𝑖
worst-case M-CVaR criterion into the objective function. First, the IRP (9)
d𝑠𝑡,𝑖 𝑠
+ 𝛼𝑡,𝑖 = I𝑠𝑡−1,𝑖 𝑠 ,
+ 𝑞𝑡,𝑖 + 𝛿𝑡,𝑖 ∀𝑠 ∈ 𝒮 , 𝑡 ∈ {2, … , 𝑇 }, 𝑖 ∈ 𝒩 . (b)
decision-maker usually can only obtain partial distributional informa-
tion of uncertain demands and is unable to compute the precise mean The inventory balance constraint in the 𝑠th demand scenario and
value and CVaR value of the cost. If the worst-case criterion over the 𝑡th period for retailer 𝑖 is developed as follows:
the developed ambiguity sets is employed, the IRP decision-maker can
∑
𝑡 ∑
𝑡
obtain a robust solution immune to the influence of the ambiguities I𝑠𝑡,𝑖 = I0,𝑖 + 𝑞𝑡′ ,𝑖 − d𝑠𝑡′ ,𝑖 , ∀𝑠 ∈ 𝒮 , 𝑡 ∈ 𝒯 , 𝑖 ∈ 𝒩 . (10)
in distributions. This prompts us to consider the worst-case criterion 𝑡′ =1 𝑡′ =1
for the inventory and shortage penalty costs. Second, for a risk-averse The retailer’s inventory level cannot exceed its maximum inventory
IRP decision-maker, even though he or she hates risk, the delivery level, which is regulated by the following constraint:
solution might not be made solely based on risk aversion in most cases.
It is advisable to balance the mean and risk level of the cost [60]. I𝑠𝑡,𝑖 ≤ 𝐶𝑖 , ∀𝑠 ∈ 𝒮 , 𝑡 ∈ 𝒯 , 𝑖 ∈ 𝒩 . (11)
Third, the IRP decision-maker can alternatively put the sum of the
The types of variables 𝒚, 𝒙, 𝒒, 𝐅, 𝐈, 𝜶, and 𝜷 are specified by the
worst-case mean value of inventory and shortage penalty costs in the
following constraint:
objective function and formulate the following distributionally robust
[ ]
CVaR constraints: 𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) ≤ ℜ𝑖 , 𝐏𝑖 ∈ 𝒫𝑖 , 𝑖 ∈ 𝒩 , where 𝒚 ∈ {0, 1}𝑇 ×𝑁 , 𝒙 ∈ {0, 1}𝑇 ×|𝒜 | , 𝒒 ∈ R𝑇+×𝑁 , (12)
parameter ℜ𝑖 is the upper bound estimated by the IRP decision-maker 𝑆×𝑇 ×𝑁 𝑆×𝑇 ×𝑁 ×𝑁
𝐅∈ R𝑆×𝑁
+ , 𝐈∈R , 𝜶∈ R+ , 𝜹∈ R𝑆×𝑇
+ , (13)
with respect to the CVaR value of inventory and shortage penalty costs
of node 𝑖 over the developed ambiguity set. However, it is difficult where 𝐅 = (F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ))𝑠∈𝒮 ,𝑖∈𝒩 , 𝐈 = (I𝑠𝑡,𝑖 )𝑠∈𝒮 ,𝑡∈𝒯 ,𝑖∈𝒩 , 𝜶 = (𝛼𝑡,𝑖
𝑠 )
𝑠∈𝒮 ,𝑡∈𝒯 ,𝑖∈𝒩 ,
for the IRP decision-maker to accurately estimate ℜ = (ℜ𝑖 )𝑖∈𝒩 . The 𝑠 )
and 𝜷 = (𝛽𝑡,𝑖 .
𝑠∈𝒮 ,𝑡∈𝒯 ,𝑖∈𝒩
above reasons motivate us to introduce the balance parameter 𝛽𝑖 and
incorporate the worst-case M-CVaR criterion into the objective function
3.2.3. Formulating distributionally robust M-CVaR IRP model
to seek a delivery solution that trade-offs the mean and CVaR.
Given the above problem description, we integrate the formulated
objective function and constraints to develop the following distribution-
ally robust M-CVaR IRP model:
3.2.2. Constructing constraints
The following constraint (1) specifies the visit degree for each ⎧min TC1 (𝒙, 𝒒, 𝒛, 𝒚)
retailer in the each period. ⎪ 𝜼,𝐅 [ ]
⎪ ∑ { [ ] [ ]}
∑ ∑ ⎨ + max 𝛽𝑖 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 )𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 )
𝑥𝑡,𝑗,𝑖 + 𝑥𝑡,𝑖,𝑗 = 2𝑦𝑡,𝑖 , ∀𝑡 ∈ 𝒯 , 𝑖 ∈ 𝒩 . (1) ⎪ 𝐏𝑖 ∈𝒫𝑖
𝑖∈𝒩
𝑗∈𝒱 ,𝑗<𝑖 𝑗∈𝒱 ,𝑗>𝑖 ⎪
⎩ s.t. Constraints ∶ (1)–(13),
The number of vehicles back to the distribution center equals the
number of vehicles leaving the distribution center and this cannot where 𝜼 = (𝒙, 𝒒, 𝒚, 𝒛, 𝐈, 𝜶, 𝜷).
exceed the total number of vehicles. We reformulate this constraint as Because 𝐏𝑖 ∈ 𝒫𝑖 , the model constructed above is a semi-infinite
∑ ∑ programming model that is usually difficult to solve directly. In the fol-
𝑥𝑡,0,𝑗 = 𝑥𝑡,𝑗,𝑁+1 ≤ 𝐾, ∀𝑡 ∈ 𝒯 . (2) lowing sections, we will construct appropriate ambiguity sets, analyze
𝑗∈𝒩 𝑗∈𝒩 the characteristics of the above model, and derive its tractable form.
Constraint (3) represents that the total flow equals the maximum
load of vehicle if and only if the edge (𝑖, 𝑗) is traversed by a vehicle in 4. Using limited historical data to construct ambiguity set
the 𝑡th period.
Most existing studies used the following two approaches to charac-
𝑧𝑡,𝑖,𝑗 + 𝑧𝑡,𝑗,𝑖 = 𝑄𝑥𝑡,𝑖,𝑗 , ∀ 𝑡 ∈ 𝒯 , (𝑖, 𝑗) ∈ 𝒜 . (3) terize demands uncertainties for retailers. First, the demand scenarios

6
Y. Feng et al. Omega 127 (2024) 103082

were generated based on the distribution specified in advance to char- constructed ambiguity sets also allow the IRP decision-maker to con-
acterize the demands uncertainties. In this approach, it is assumed trol the robustness of delivery solution according to their preferences.
that the probability of occurrence of each demand scenario is known Third, the distributionally robust counterpart can be transformed into
and remains constant. The second method uses a specified stochastic computationally tractable forms under these two ambiguity sets. The
or fuzzy distribution to characterize the demands uncertainties. The following contents will confirm these three points.
implementation of these two approaches relies on the strict assumption
that the distributions of demands for retailers are known. 4.2. Framework of constructing ambiguity sets
Considering this, for each retailer, we do not apply the above strict
assumption and employ adding a disturbance to the nominal demand We describe our idea by designing the following framework. To
distribution to construct ambiguity sets to characterize the ambiguity construct the set of demand scenarios and calculate the probabili-
of distribution. Specifically, we first use limited historical demand ties of occurrences of demand scenarios for each retailer, we first
data of each retailer to estimate its nominal demand distribution. The need to capture some representative sampling demand points based
disturbance vector is then used to characterize the perturbation. on historical demand data, and then calculate the probability of each
sampling demand point. Subsequently, the demand scenarios faced by
each retailer are constructed based on the sampling demand points. The
4.1. The forms of constructed ambiguity sets
probability of occurrence of each demand scenario is calculated based
on the probabilities of the sampling demand points. To refine the set
Two different types of ambiguity sets are constructed to characterize of demand scenarios while reducing model complexity, we select the
the ambiguity in distribution. For retailer 𝑖, 𝑖 ∈ 𝒩 , the first ambiguity more concerned demand scenarios from the constructed set to form a
set is denoted by 𝒫𝑖1 and is expressed as follows: new set of demand scenarios. Finally, we construct ambiguity sets based
{ } on the regularized probability vector of occurrence of the selected
𝒫𝑖1 = 𝐏𝑖 ∣ 𝐏𝑖 = 𝐏0𝑖 + 𝐀𝑖 𝜻 𝑖 , 𝒆𝖳 𝐀𝑖 𝜻 𝑖 = 0, 𝐏0𝑖 + 𝐀𝑖 𝜻 𝑖 ≥ 𝟎, ‖𝜻 𝑖 ‖ ≤ 𝛤𝑖 , (14)
demand scenarios and set the scale parameters for the ambiguity sets.
where 𝐏0𝑖 denotes the nominal (i.e., most likely) probability vector of We outline the following five steps for constructing the aforementioned
the occurrence of demand scenarios for retailer 𝑖, 𝐀𝑖 ∈ R𝑆×𝑆 represents ambiguity sets.
scaling matrix, 𝜻 𝑖 denotes the disturbance vector, and 𝒆 represents the First, the representative sampling demand points are captured. We
unit column vector. use low demand (dL𝑖 ), medium demand (dM 𝑖 ), or high demand (d𝑖 )
H

to represent the representative sampling demand points for retailer 𝑖

Remark 3. In 𝒫𝑖1 , the specific forms of ‖𝜻 𝑖 ‖ can be set to ‖𝜻 𝑖 ‖1 in each period. The IRP decision-maker can set more representative
and ‖𝜻 𝑖 ‖2 . The tractable reformulation of the worst-case M-CVaR is a sampling demand points. However, the more representative sampling
linear programming form and second-order cone programming (SOCP) demand points are captured, the higher the complexity of constructing
form based on ‖𝜻 𝑖 ‖1 and ‖𝜻 𝑖 ‖2 , respectively. The decision-maker can the set of scenarios. This can be observed from the following content.
regulate the robustness of a decision by adjusting 𝛤𝑖 . Here, we describe We now demonstrate how to capture representative sampling demand
this procedure. We assume that the dimension of 𝜻 is 2 and draw the points.
diagram (Fig. 2(a)) of ‖𝜻‖1 ≤ 𝛤 with 𝛤 = 0.8, 1, 1.2. The range of Based on the historical demand data set (denote D𝑖 ) of retailer 𝑖, the
‖𝜻‖1 ≤ 𝛤 increases as 𝛤 increases. The scale of 𝒫𝑖1 also increases. Under IRP decision-maker sorts the data in D𝑖 in increasing order, dividing the
interval [dmin max
larger-scale 𝒫𝑖1 , 𝑖 ∈ 𝒩 , the robustness of the decision is higher. Hence, 𝑖 , d𝑖 ] into three equal-length and disjoint subintervals,
the IRP decision-maker can modulate the robustness of the delivery the midpoint of the first subinterval being dL𝑖 , the midpoint of the
solution by modulating the value of 𝛤𝑖 , 𝑖 ∈ 𝒩 . second subinterval being dM 𝑖 and the midpoint of the third subinterval
being dH𝑖 . Therefore, the following formulas can be established:
Alternatively, we consider not using the norm of 𝜻 𝑖 to limit its
dmax
𝑖 − dmin
𝑖
disturbance range but instead setting upper and lower bounds. In this dL𝑖 = dmin + , (16)
𝑖 6
context, for retailer 𝑖, 𝑖 ∈ 𝒩 , we construct the second ambiguity set,
which is denoted by 𝒫𝑖2 and expressed as follows: dmax
𝑖 − dmin
𝑖
dM min
𝑖 = d𝑖 + , (17)
2
2
𝒫𝑖 = {𝐏𝑖 ∣ 𝐏𝑖 = 𝐏0𝑖 𝖳
+ 𝜻 𝑖 , 𝒆 𝜻 𝑖 = 0, 𝜻𝐿
𝑖 ≤ 𝜻𝑖 ≤ 𝜻𝑈
𝑖 }, (15) 5(dmax
𝑖 − dmin
𝑖 )
dH min
𝑖 = d𝑖 + , (18)
where 𝜻𝐿
and 𝜻𝑈 denote the upper and lower bounds of the dis- 6
𝑖 𝑖
turbance vector 𝜻 𝑖 , respectively. Furthermore, 𝜻 𝐿 0 where dmax
and dmin
are the maximum and minimum values in set D𝑖 ,
𝑖 ≥ −𝐏𝑖 guarantees 𝑖 𝑖
nonnegativity of 𝐏𝑖 . respectively.
Second, the probabilities of representative sampling demand points
Remark 4. We can introduce the scale parameters 𝜗𝑖 into the upper are calculated. For retailer 𝑖, letting P(dL𝑖 ), P(dM H
𝑖 ), and P(d𝑖 ) denote the
and lower bounds of 𝜻 𝑖 such that the IRP decision-maker can adjust probabilities of occurrences for dL𝑖 , dM
𝑖 , and dH , respectively. Naturally,
𝑖
the scale of 𝒫𝑖2 to modulate the robustness of the delivery solution the IRP decision-maker can obtain P(dL𝑖 ), P(dM H
𝑖 ), and P(d𝑖 ) by calcu-
by modulating 𝜗𝑖 . Below, we demonstrate the realization of this point lating the ratio of the data volume belonging to the corresponding
using the following example: we assume that the dimension of 𝜻 is 2 and subintervals to the total data volume in D𝑖 . Therefore, the following
𝐏0 = (0.6, 0.4), and set −𝜗𝐏0 ≤ 𝜻 ≤ 𝜗𝐏0 . We draw the diagram (Fig. 2(b)) formulas can be established:
of these three different upper and lower bounds under 𝜗 = 0.4, 0.5, and dmax −dmin
Number[dmin min +
𝑖 , d𝑖
𝑖
3
𝑖
]
0.6. As 𝜗 increases, the range of −𝜗𝐏0 ≤ 𝜻 ≤ 𝜗𝐏0 increases. This further P(dL𝑖 ) = , (19)
|D𝑖 |
increases the scale of ambiguity set 𝒫𝑖2 . Based on the larger-scale 𝒫𝑖2 ,
𝑖 ∈ 𝒩 , the robustness of the obtained delivery solution is higher. dmax −dmin 2(dmax −dmin
𝑖 )
Number[dmin
𝑖 + 𝑖
3
𝑖
, dmin
𝑖 + 𝑖
3
]
P(dM
𝑖 ) = , (20)
The developed ambiguity sets are relevant and useful for the prob- |D𝑖 |
lem we studied. The reasons for this are as follows. First, the probability 2(dmax −dmin
𝑖 )
Number[dmin
𝑖 + 𝑖
3
, dmax
𝑖 ]
distribution of demand is discrete for each retailer. Ambiguity sets P(dH (21)
𝑖 )= |D𝑖 |
,
are constructed by adding a disturbance to the discrete nominal de-
mand distribution. Therefore, the aforementioned ambiguity sets can where |D𝑖 | is the cardinality of D𝑖 , Number[𝑎, 𝑏] denotes the number of
characterize the distributions ambiguities for retailers. Second, the elements in the interval [𝑎, 𝑏] and D𝑖 .

7
Y. Feng et al. Omega 127 (2024) 103082

Fig. 2. The diagrams of ‖𝜻‖1 ≤ 𝛤 and 𝜻 𝐿 ≤ 𝜻 ≤ 𝜻 𝑈 .

Table 2 the disturbance range information of the disturbance vector to some

Demand scenarios faced by retailer 𝑖 with 𝑇 = 2.
extent. Herein, to provide more flexibility for the IRP decision-maker to
0 0
Demand scenario P𝑖,𝑠 Demand scenario P𝑖,𝑠 regulate the scales of the ambiguity sets 𝒫𝑖1 and 𝒫𝑖2 , we introduce scale
DS1,𝑖 = {dL𝑖 , dL𝑖 } P(dL𝑖 ) ⋅ P(dL𝑖 ) DS6,𝑖 = {dM , dH } P(dM ) ⋅ P(dH ) parameters 𝜒𝑖 and 𝜗𝑖 into the scaling matrix and bounds, respectively.
𝑖 𝑖 𝑖 𝑖
DS2,𝑖 = {dL𝑖 , dM P(dL𝑖 ) ⋅ P(dM DS7,𝑖 = {dH L
P(dH L The specific operations are as follows: we let 𝐀𝑖 = 𝜒𝑖 ⋅ M (M represents
𝑖 } 𝑖 ) 𝑖 , d𝑖 } 𝑖 ) ⋅ P(d𝑖 )
DS3,𝑖 = {dL𝑖 , dH } P(dL𝑖 ) ⋅ P(dH ) DS8,𝑖 = {dH , dM } P(dH ) ⋅ P(dM ) identify matrix) in the 𝒫𝑖1 , and let 𝜻 𝑈 𝐿 0
𝑖 = −𝜻 𝑖 = 𝜗𝑖 ⋅ 𝐏𝑖 in the 𝒫𝑖 .
2
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
DS4,𝑖 = {dM L
P(dM L
DS9,𝑖 = {dH H
P(dH H In general, 𝜒𝑖 and 𝜗𝑖 take values in the interval [0, 1]. In this context,
𝑖 , d𝑖 } 𝑖 ) ⋅ P(d𝑖 ) 𝑖 , d𝑖 } 𝑖 ) ⋅ P(d𝑖 )
DS5,𝑖 = {dM , dM } P(dM ) ⋅ P(dM ) the constructed ambiguity sets possess the following properties: ① if
𝑖 𝑖 𝑖 𝑖
parameters 𝜒𝑖 = 0 and 𝜗𝑖 = 0, the ambiguity sets only contain
the nominal distribution and the developed distributionally robust M-
CVaR IRP model degenerates to the stochastic M-CVaR IRP model;
Third, demand scenarios are constructed, and the initial probabili- ② if parameters 𝜒𝑖 , 𝜗𝑖 → 1, the scales of 𝒫𝑖1 and 𝒫𝑖2 become larger.
ties of occurrences of these scenarios are calculated. When the number Additionally, for the ambiguity set 𝒫𝑖1 , the value of the parameter 𝛤𝑖
of periods is 𝑇 , the IRP decision-maker can construct 3𝑇 demand can also modulate its scale. By performing these operations, the scale
scenarios (DS𝑠,𝑖 , 𝑠 ∈ {1, … , 3𝑇 }) for retailer 𝑖, 𝑖 ∈ 𝒩 . The demand of the ambiguity sets could be regulated by changing the values of the
scenario DS𝑠,𝑖 comprises 𝑇 elements, each of these 𝑇 elements is dL𝑖 , dM
𝑖 , scale parameters.
0
or dH
𝑖 . For retailer 𝑖, denoting P𝑖,𝑠 as the nominal value of the probability
We employ the matrix 𝐃 to represent the entire historical demand
0 data for all retailers. In this matrix, the number of rows corresponds
of the 𝑠th demand scenario, P𝑖,𝑠 can be calculated by the following
formula: to the number of retailers, and the number of columns represents the
∏ amount of demand data. In other words, each row of 𝐃 represents the
0
P𝑖,𝑠 = P(𝓁𝑖,𝑡 ), (22) historical demand data of a retailer. 𝐃𝐒𝑖 = (DS𝑖,𝑠 )𝑠∈𝒮 represents the set
𝑡∈𝒯 of final demand scenarios for retailer 𝑖, ∀𝑖 ∈ 𝒩 . We develop Algorithm 1
where 𝓁𝑖,𝑡 is the 𝑡th element of scenario DS𝑠,𝑖 , P(𝓁𝑖,𝑡 ) is the probability (incorporating the 5 steps mentioned above) to obtain {𝒫𝑖1 , 𝒫𝑖2 , 𝐃𝐒𝑖 }𝑖∈𝒩
of 𝓁𝑖,𝑡 . based on 𝐃.
Here, we consider the case in which retailer 𝑖 with 𝑇 = 2 to In this section, we have developed two different ambiguity sets. If
demonstrate its demand scenarios. In this context, there are nine de- the decision-maker can gather more information about the scale matrix
mand scenarios: DS𝑠,𝑖 , 𝑠 ∈ {1, … , 9}. We denote the probability of in detail, the ambiguity set 𝒫𝑖1 can be used to characterize distributions
0 ambiguities. If the decision-maker can gather the information about the
occurrence of scenario DS𝑠,𝑖 as P𝑖,𝑠 . The nine demand scenarios and
their probabilities are listed in Table 2. upper and lower bounds of the disturbance vector, the ambiguity set 𝒫𝑖2
can be exploited to characterize distributions ambiguities. Of course,
Fourth, scenarios selection and probability vector regularization.
the decision-maker can also directly use the scale parameters to set
A larger 𝑇 results in the generation of a larger number of demand
the scale matrix or the bounds and use the 𝒫𝑖1 or 𝒫𝑖2 to characterize
scenarios. For each retailer, we rank the demand scenarios according
distributions ambiguities. In the next section, we derive the tractable
to their probability of occurrence and select the first 𝑆 scenarios from
form of the presented distributionally robust M-CVaR IRP model under
3𝑇 scenarios. After selecting the scenarios, the probability vector of the
the constructed ambiguity sets.
selected demand scenarios does not satisfy the regularization. There-
fore, we construct the following normalization equation to normalize 5. Model analysis and processing
this set of probabilities:
0 We present some important analyses of the developed distribution-
P𝑖,𝑠
P0𝑖,𝑠 = , (23) ally robust M-CVaR IRP model and reformulate it.
∑𝑆 0
𝑠=1 P𝑖,𝑠
5.1. Model analysis
where P0𝑖,𝑠 represents the nominal probability of occurrence of demand
scenario 𝑠 for retailer 𝑖 after normalization. In this part, we analyze the characteristics of the constructed model.
Fifth, the scale parameters of the ambiguity sets are set. The scaling First, the objective value of the distributionally robust M-CVaR IRP
matrix 𝐀𝑖 in 𝒫𝑖1 , the upper bound 𝜻 𝑈 𝐿 2
𝑖 and lower bound 𝜻 𝑖 in 𝒫𝑖 reflect model provides an upper bound for the CVaR value of the total cost

8
Y. Feng et al. Omega 127 (2024) 103082

Algorithm 1: Constructing ambiguity sets and returning demand Second, the constructed distributionally robust M-CVaR IRP model
scenarios. is flexible. The decision-maker only needs to vary the inputs appro-
Input: 𝐃, 𝑆, 𝑇 , 𝑁; priately so that the model can be degraded to one based on different
Output: {𝒫𝑖1 , 𝒫𝑖2 , 𝐃𝐒𝑖 }𝑖∈𝒩 ; decision perspectives. By varying the values of the balance parameters,
for 𝑖 = 1 to 𝑁 do 𝛽𝑖 , 𝑖 ∈ 𝒩 , our distributionally robust M-CVaR IRP model can degenerate
Generating ten variables: dmax , dmin , dL , dM , dH , P(dL ), P(dM ), into the distributionally robust IRP model with different risk types. By
P(dH ), 𝜒, 𝜗; setting the scale parameters of the ambiguity sets (14)–(15) to 0, the
Generating two empty cells: 𝐃𝐒𝟏 and 𝐃𝐒 (the dimensions are distributionally robust IRP model can degenerate into the stochastic IRP
3𝑇 and 𝑆); model. Detailed cases are listed in Table 3.
Generating two empty arrays: 𝐏𝟏 and 𝐏0 (the dimensions are
3𝑇 and 𝑆); 5.2. Model processing
dmax ← the maximum value of 𝐃(𝑖, ∶);
dmin ← the minimum value of 𝐃(𝑖, ∶); In this part, we first re-express the objective function and then give
dL , dM , dH ← applying eqs. (16)-(18) based on dmax and dmin ; the tractable reformulation of the proposed model under the above two
P(dL ), P(dM ), P(dH ) ← applying eqs.(19)-(21); different ambiguity sets constructed.
for 𝑠 = 1 to 3𝑇 do
𝐃𝐒𝟏(𝑠) ← using the idea in Table 2 based on dL , dM , and dH 5.2.1. Re-expressing objective function
to construct demand scenario; By referring to the definition of 𝐂𝐕𝐚𝐑𝜅,𝐏 given by Rockafellar &
𝐏𝟏(𝑠) ← using eq.(22) based on 𝐃𝐒𝟏(𝑠), P(dL ), P(dM ), and Uryasev [20], we can reformulate the CVaR value of the 𝑖th retailer,
[ ]
P(dH ) to calculate probability; that is, 𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) as
Ranking the elements of 𝐏𝟏 in decreasing order. { }
1 [ ]
Ranking the elements of 𝐃𝐒𝟏 refer to the sequence of elements min 𝜙𝑖 + 𝐄𝐏𝑖 max{F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) − 𝜙𝑖 , 0} ,
𝜙𝑖 ∈R 1 − 𝜅𝑖
ranking in 𝐏𝟏.
𝐃𝐒 ← selecting the first 𝑆 elements of 𝐃𝐒𝟏; where the confidence level parameter 𝜅𝑖 belongs to (0, 1) and reflects
𝐏0 ← selecting the first 𝑆 elements of 𝐏𝟏 and regularizing them the probability that F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) is lower than 𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 .
{ [ ]
using eq.(23); In light of the above definition, max𝐏𝑖 ∈𝒫𝑖 𝛽𝑖 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 )
[ ] }
𝜒,𝜗 ← taking two suitable values from [0, 1];
𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) is equivalent to
𝛤 ← taking an appropriately sized value;
{ { }}
Constructing 𝒫𝑖1 and 𝒫𝑖2 based on {𝐏0 , 𝜒, 𝛤 } and {𝐏0 , 𝜗}, [ ] 1 [ ]
max 𝛽𝑖 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 ) min 𝜙𝑖 + 𝐄 max{F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) − 𝜙𝑖 , 0} .
respectively; 𝐏𝑖 ∈𝒫𝑖 𝜙𝑖 ∈R 1 − 𝜅𝑖 𝐏𝑖
Based on the Corollary 37.3.2 given by Rockafellar [61], we can
exchange the order of the operators max𝐏𝑖 ∈𝒫𝑖 and min𝜙𝑖 ∈R . Hence, the
worst-case M-CVaR value of the sum of inventory and shortage penalty
based on 𝛽𝑖 = 0, 𝑖 ∈ 𝒩 . We explain this point in detail. We use TC to de- costs for retailer 𝑖 is reformulated as
∑ { { }}
note the total cost of the IRP, i.e., TC = TC1 (𝒙, 𝒒, 𝒛, 𝒚) + 𝑖∈𝒩 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ), [ ] 1 [ ]
min max 𝛽𝑖 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 ) 𝜙𝑖 + 𝐄 max{F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) − 𝜙𝑖 , 0} .
where F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) is a discrete random variable. The worst-case CVaR 𝜙𝑖 ∈R 𝐏𝑖 ∈𝒫𝑖 1 − 𝜅𝑖 𝐏𝑖
value of TC is max𝐏∈𝒫 𝐂𝐕𝐚𝐑𝜅,𝐏 [TC], where 𝜅 = 𝜅𝑖 for 𝑖 ∈ 𝒩 , and 𝒫 =
{ ⨂ } Given the retailer 𝑖 and demand scenario 𝑠, we let auxiliary variable
𝐏 ∣ 𝐏 = 𝑖∈𝒩 𝐏𝑖 , 𝐏𝑖 ∈ 𝒫𝑖 , 𝑖 ∈ 𝒩 . It is difficult for the IRP decision-
𝜔𝑠𝑖 represent max{F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) − 𝜙𝑖 , 0}. Then, the above formula can be
maker to compute max𝐏∈𝒫 𝐂𝐕𝐚𝐑𝜅,𝐏 [TC] because of the inability to
reformulated as
derive its tractable form. However, when the balance parameter 𝛽𝑖 = 0, { { }}
𝑖 ∈ 𝒩 , the developed distributionally robust M-CVaR IRP model degen- ∑ 1 ∑ 𝑠 𝑠
𝑠 𝑠 𝑠 𝑠
min max 𝛽𝑖 P𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 ) 𝜙𝑖 + P𝜔
erates into the distributionally robust CVaR IRP model, whose optimal 𝜙𝑖 ∈R 𝐏𝑖 ∈𝒫𝑖
𝑠∈𝒮
1 − 𝜅𝑖 𝑠∈𝒮 𝑖 𝑖
objective value provides an upper bound on max𝐏∈𝒫 𝐂𝐕𝐚𝐑𝜅,𝐏 [TC]. This { { }}
1
can be verified by the following relationship: = min max 𝛽𝑖 𝐅𝖳𝑖 𝐏𝑖 + (1 − 𝛽𝑖 ) 𝜙𝑖 + 𝝎𝖳 𝐏 ,
𝜙𝑖 ∈R 𝐏𝑖 ∈𝒫𝑖 1 − 𝜅𝑖 𝑖 𝑖
[ ]
∑ where 𝐅𝑖 = (F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ))𝑠∈𝒮 and 𝝎𝑖 = (𝜔𝑠𝑖 )𝑠∈𝒮 .
max 𝐂𝐕𝐚𝐑𝜅,𝐏 [TC] = TC1 (𝒙, 𝒒, 𝒛, 𝒚) + max 𝐂𝐕𝐚𝐑𝜅,𝐏 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) After performing the above transformation, the proposed model is
𝐏∈𝒫 𝐏∈𝒫
𝑖∈𝒩
∑ [ ] reformulated as follows:
≤ TC1 (𝒙, 𝒒, 𝒛, 𝒚) + max 𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) . ∑
𝐏𝑖 ∈𝒫𝑖 ⎧min TC1 (𝒙, 𝒒, 𝒛, 𝒚) + (1 − 𝛽𝑖 ) 𝑖∈𝒩 𝜙𝑖
𝑖∈𝒩
⎪ 𝜼,𝐅,𝝓,𝝎
[ { }]
(24) ⎪ ∑ 𝖳 1−𝛽𝑖 𝖳
⎪ + 𝑖∈𝒩 max𝐏𝑖 ∈𝒫𝑖 𝛽𝑖 𝐅𝑖 𝐏𝑖 + 1−𝜅𝑖 𝝎𝑖 𝐏𝑖 (a)
We prove that inequality (24) holds in the following. For the ambiguity ⎪
⎨ s.t. 𝜔𝑠𝑖 ≥ F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) − 𝜙𝑖 , 𝜙𝑖 ∈ R, ∀𝑖 ∈ 𝒩 , 𝑠 ∈ 𝒮 , (b) (25)
set 𝒫𝑖 of node 𝑖, we denote the worst-case distribution as 𝐏∗𝑖 . Then, ⎪
[ ] [ ] ⎪
we obtain max𝐏𝑖 ∈𝒫𝑖 𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) = 𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏∗ F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) . Letting 𝜔𝑠𝑖 ≥ 0, ∀𝑖 ∈ 𝒩 , 𝑠 ∈ 𝒮 , (c)
⨂ 𝑖 ⎪
𝐏∗ = 𝑖∈𝒩 𝐏∗𝑖 , based on the condition that the demands of all retailers ⎪ Constraints ∶ (1)–(13), (d)
⎩
are independent of each other, we can obtain 𝐏∗ is the worst-case
[∑ ]
distribution in 𝒫 . The equation max𝐏∈𝒫 𝐂𝐕𝐚𝐑𝜅,𝐏 𝑖∈𝒩 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) = where 𝝓 = (𝜙𝑖 )𝑖∈𝒩 and 𝝎 = (𝝎𝑖 )𝑖∈𝒩 .
[∑ ]
𝐂𝐕𝐚𝐑𝜅,𝐏∗ 𝑖∈𝒩 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) holds. According to subadditivity, we obtain
[ ] 5.2.2. Deriving robust counterparts
∑ ∑ [ ] This part provides the computable form of objective (25a). For given
𝐂𝐕𝐚𝐑𝜅,𝐏∗ F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) ≤ 𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏∗ F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) .
𝑖 retailer 𝑖, 𝑖 ∈ 𝒩 , we must derive the equivalent form of the following
𝑖∈𝒩 𝑖∈𝒩
problem based on the ambiguity sets (14)–(15):
[∑ ] ∑ { }
1 − 𝛽𝑖 𝖳
Therefore, we prove that max𝐏∈𝒫 𝐂𝐕𝐚𝐑𝜅,𝐏 𝑖∈𝒩 F𝑖 (𝜶 , 𝜹 ) ≤ 𝑖∈𝒩 max 𝛽𝑖 𝐅𝖳𝑖 𝐏𝑖 + 𝝎𝑖 𝐏𝑖 , (26)
[ ] ∑𝑖 𝑖 𝐏𝑖 ∈𝒫𝑖 1 − 𝜅𝑖
max𝐏𝑖 ∈𝒫𝑖 𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) holds. TC1 (𝒙, 𝒒, 𝒛, 𝒚) + 𝑖∈𝒩 max𝐏𝑖 ∈𝒫𝑖
[ ]
𝐂𝐕𝐚𝐑𝜅𝑖 ,𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) is an upper bound of max𝐏∈𝒫 𝐂𝐕𝐚𝐑𝜅,𝐏 [TC]. where the decision variable is the distribution 𝐏𝑖 .

9
Y. Feng et al. Omega 127 (2024) 103082

Table 3
Degenerate version of the developed distributionally robust M-CVaR IRP model under different inputs.
Risk-neutral Risk-averse
𝜷=𝟏 𝟏>𝜷>𝟎 𝜷=𝟎
𝜞,𝝑 = 𝟎 Stochastic Mean IRP model Stochastic M-CVaR IRP model Stochastic CVaR IRP model
𝝌, 𝜞 , 𝝑 ≥ 𝟎 Distributionally robust Mean IRP model Distributionally robust M-CVaR IRP model Distributionally robust CVaR IRP model

Note: 𝝌 = (𝜒𝑖 )𝑖∈𝒩 , 𝜞 = (𝛤𝑖 )𝑖∈𝒩 , 𝝑 = (𝜗𝑖 )𝑖∈𝒩 , and 𝜷 = (𝜗𝑖 )𝑖∈𝒩 .

∗ ∗ ∗
∗
First, the tractable form of (26) with ambiguity set (14) is given by Proposition 1. If the solution (𝜼∗ , 𝝓∗ , 𝐅 , 𝝎 , 𝝂 ∗ , 𝜰 , 𝝅 ∗ , 𝝂̂ ∗ , 𝜰̂ , 𝝅̂ ∗ ) is
Theorem 1. ∗ ∗ ∗ ∗
optimal for model (27), then (𝜼 , 𝝓 , 𝐅 , 𝝎 ) optimizes model (25). Con-
versely, if (𝜼1∗ , 𝝓1∗ , 𝐅1∗ , 𝝎1∗ ) is the optimal solution of model (25), then
Theorem 1. The maximization problem (26) with ambiguity set (14) can ̂ ∗ , 𝜰̂ 1∗ , 𝝅1
(𝜼1∗ , 𝝓1∗ , 𝐅1∗ , 𝝎1∗ , 𝝂1∗ , 𝜰 1∗ , 𝝅1∗ , 𝝂1 ̂ ∗ ) optimizes the model (27),
be transformed into the following minimization problem: and the solution (𝐅1∗ , 𝝂1∗ , 𝜰 1∗ , 𝝅1∗ ) is optimal for
⎧ { } 1 − 𝛽𝑖 { 𝖳 0 𝖳
}
⎧min ∑{ }
⎫
⎪min 𝛽𝑖 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈𝑖 + 𝝎𝑖 𝐏𝑖 + 𝜰̂ 𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈̂𝑖 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈𝑖
⎪ 1 − 𝜅𝑖 ⎪ ⎪
∑ ⎪ 𝑖∈𝒩 ⎪
⎪ s.t. ‖𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 ‖∗ ≤ 𝜈𝑖 , max 𝐅𝖳𝑖 𝐏𝑖 = ⎨ s.t.
⎨ 1 ‖𝐀𝑖 𝐅𝑖 + 𝐀𝑖 𝜰 𝑖 − 𝐀𝑖 𝒆𝜋𝑖 ‖∗ ≤ 𝜈𝑖 , ∀𝑖 ∈ 𝒩 ,⎬
𝖳 𝖳 𝖳
⎪ 𝑖∈𝒩 𝐏𝑖 ∈𝒫𝑖 ⎪ ⎪
‖𝐀𝖳 𝝎𝑖 + 𝐀𝖳 𝜰̂ 𝑖 − 𝐀𝖳 𝒆𝜋̂ 𝑖 ‖∗ ≤ 𝜈̂𝑖 , ⎪ 𝜰 𝑖 ∈ R𝑆+ , 𝜋𝑖 ∈ R, 𝜈𝑖 ∈ R+ , ∀𝑖 ∈ 𝒩 , ⎪
⎪ 𝑖 𝑖 𝑖
⎩ ⎭
⎪ 𝜰 𝑖 ∈ R𝑆+ , 𝜋𝑖 ∈ R, 𝜈𝑖 ∈ R+ , 𝜰̂ 𝑖 ∈ R𝑆+ , 𝜋̂ 𝑖 ∈ R, 𝜈̂𝑖 ∈ R+ ,
⎩ (28)
where decision variables are 𝜰 𝑖 , 𝜋𝑖 , 𝜈𝑖 , 𝜰̂ 𝑖 , 𝜋̂ 𝑖 , 𝜈̂𝑖 . ‖ ⋅ ‖∗ is the dual norm
and the solution ̂ ∗ , 𝜰̂ 1∗ , 𝝅1
(𝝎1∗ , 𝝂1 ̂ ∗) is optimal for
of ‖ ⋅ ‖. This problem is a linear programming form and SOCP form under
‖𝜻 𝑖 ‖1 and ‖𝜻 𝑖 ‖2 , respectively. ⎧ ⎧ 𝖳 ⎫ ⎫
⎪ ∑ ⎪ 𝝎𝖳𝑖 𝐏0𝑖 + 𝜰̂ 𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈̂𝑖 ⎪ ⎪
⎪min ⎨ ⎬ ⎪
Proof. See Appendix A for this proof process. □ ∑ 𝝎𝖳𝑖 𝐏𝑖 ⎪ 𝑖∈𝒩 ⎪
1 − 𝜅𝑖 ⎪ ⎪
max
1 1−𝜅
=⎨ ⎩ ⎭ ⎬
Second, the tractable form of (26) with ambiguity set (15) is given 𝑖∈𝒩 𝐏𝑖 ∈𝒫𝑖 𝑖 ⎪ ⎪
⎪ s.t. ‖𝐀𝖳 𝝎𝑖 + 𝐀𝖳 𝜰̂ 𝑖 − 𝐀𝖳 𝒆𝜋̂ 𝑖 ‖∗ ≤ 𝜈̂𝑖 ,
𝑖 𝑖 𝑖 ∀𝑖 ∈ 𝒩 ,⎪
by Theorem 2. ⎪ ⎪
⎩ 𝜰̂ 𝑖 ∈ R𝑆+ , 𝜋̂ 𝑖 ∈ R, 𝜈̂𝑖 ∈ R+ , ∀𝑖 ∈ 𝒩 . ⎭
Theorem 2. The maximization problem (26) with ambiguity set (15) can (29)
be transformed into the following linear programming problem:
⎧min { }
⎪ 𝛽𝑖 𝐅𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 𝐿 𝖳
𝑖 ) 𝝅 𝑖 − (𝜻 𝑖 ) 𝜸 𝑖
⎪ Proof. See Appendix C for this proof process. □
1 − 𝛽𝑖 { 𝖳 0 𝖳 𝖳
}
⎪ + 𝝎𝑖 𝐏𝑖 + (𝜻 𝑈 𝑖 ) 𝝅̂ 𝑖 − (𝜻 𝐿
𝑖 ) 𝜸̂ 𝑖
⎪ 1 − 𝜅𝑖 Based on the ambiguity set (15) and using Theorem 2, we can
⎨ s.t. 𝒆𝜏𝑖 + 𝝆𝑖 − 𝜸 𝑖 = 𝐅𝑖 , reformulate the model (25) as follows:
⎪
⎪ 𝒆𝜏̂𝑖 + 𝝆̂ 𝑖 − 𝜸̂ 𝑖 = 𝝎𝑖 , ⎧ ∑
⎪ ⎪min𝜼′′ TC1 (𝒙, 𝒒, 𝒛, 𝒚) + 𝑖∈𝒩 (1 − 𝛽𝑖 )𝜙𝑖
⎪ 𝜏𝑖 ∈ R, 𝝆𝑖 ∈ R𝑆+ , 𝜸 𝑖 ∈ R𝑆+ , 𝜏̂𝑖 ∈ R, 𝝆̂ 𝑖 ∈ R𝑆+ , 𝜸̂ 𝑖 ∈ R𝑆+ , ⎪ ∑ { }
⎩ ⎪ + 𝑖∈𝒩 𝛽𝑖 𝐅𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈𝑖 )𝖳 𝝆𝑖 − (𝜻 𝐿𝑖 )𝖳 𝜸 𝑖
⎪ { 𝖳 0 𝑈𝖳 }
where decision variables 𝜏𝑖 , 𝝆𝑖 , 𝜸 𝑖 , 𝜏̂𝑖 , 𝝆̂ 𝑖 , and 𝜸̂ 𝑖 . ⎪ ∑ 𝝎𝑖 𝐏𝑖 +(𝜻 𝑖 ) 𝝆̂ 𝑖 −(𝜻 𝐿 𝖳̂
𝑖 ) 𝜸 𝑖
⎪ + 𝑖∈𝒩 (1 − 𝛽𝑖 ) 1−𝜅
(a)
𝑖
⎪
Proof. See Appendix B for this proof process. □ ⎨s.t. 𝒆𝜏𝑖 + 𝝆𝑖 − 𝜸 𝑖 = 𝐅𝑖 , ∀𝑖 ∈ 𝒩 , (b)
⎪
⎪ 𝒆𝜏̂𝑖 + 𝝆̂ 𝑖 − 𝜸̂ 𝑖 = 𝝎𝑖 , ∀𝑖 ∈ 𝒩 , (c)
⎪
5.3. Computationally tractable reformulation ⎪
⎪ 𝜏𝑖 ∈ R, 𝝆𝑖 ∈ R𝑆+ , 𝜸 𝑖 ∈ R𝑆+ , 𝜏̂𝑖 ∈ R, 𝝆̂ 𝑖 ∈ R𝑆+ , 𝜸̂ 𝑖 ∈ R𝑆+ , ∀𝑖 ∈ 𝒩 , (d)
⎪
We now present the final computationally tractable reformulation of ⎪ Constraints ∶ (1)–(13), (25b)–(25c), (e)
⎩
the proposed model. Based on the ambiguity set (14), by Theorem 1,
we can reformulate model (25) as follows: (30)

⎧min ∑ ∑ { } where 𝜼′′ = (𝜼, 𝝓, 𝐅, 𝝎, 𝜏, 𝝆, 𝜸, 𝜏, ̂ 𝝆, ̂ 𝜏 = (𝜏𝑖 )𝑖∈𝒩 , 𝝆 = (𝝆𝑖 )𝑖∈𝒩 , 𝜸 =

̂ 𝜸),
⎪ 𝜼′ TC1 (𝒙, 𝒒, 𝒛, 𝒚) + 𝑖∈𝒩 (1 − 𝛽𝑖 )𝜙𝑖 + 𝑖∈𝒩 𝛽𝑖 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈𝑖
{ } (𝜸 𝑖 )𝑖∈[𝑁] , 𝜏̂ = (𝜏̂𝑖 )𝑖∈𝒩 , 𝝆̂ = (𝝆̂ 𝑖 )𝑖∈𝒩 , 𝜸̂ = (𝜸̂ 𝑖 )𝑖∈𝒩 .
⎪ ∑ 𝖳
𝝎𝖳𝑖 𝐏0𝑖 +𝜰̂ 𝑖 𝐏0𝑖 +𝛤𝑖 𝜈̂𝑖
⎪ + 𝑖∈𝒩 (1 − 𝛽𝑖 ) (a) In addition, we introduce the following valid inequalities given by
⎪ 1−𝜅 𝑖
Manousakis et al. [62], to tighten the proposed model:
⎪
⎪s.t. ‖𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 ‖∗ ≤ 𝜈𝑖 , ∀𝑖 ∈ 𝒩 , (b)
⎨ 𝑥𝑡,0,𝑖 ≤ 𝑦𝑡,𝑖 , ∀𝑡 ∈ 𝒯 , 𝑖 ∈ 𝒩 , (31)
⎪ ‖𝐀𝖳𝑖 𝝎𝑖 + 𝐀𝖳𝑖 𝜰̂ 𝑖 − 𝐀𝖳𝑖 𝒆𝜋̂ 𝑖 ‖∗ ≤ 𝜈̂𝑖 , ∀𝑖 ∈ 𝒩 , (c)
⎪ 𝑥𝑡,𝑖,𝑗 ≤ 𝑦𝑡,𝑖 , ∀𝑡 ∈ 𝒯 , 𝑖, 𝑗 ∈ 𝒩 , 𝑖 < 𝑗, (32)
⎪
⎪ 𝜰 𝑖 ∈ R𝑆+ , 𝜋𝑖 ∈ R, 𝜈𝑖 ∈ R+ , 𝜰̂ 𝑖 ∈ R𝑆+ , 𝜋̂ 𝑖 ∈ R, 𝜈̂𝑖 ∈ R+ , ∀𝑖 ∈ 𝒩 , (d)
⎪ 𝑥𝑡,𝑖,𝑗 ≤ 𝑦𝑡,𝑗 , ∀𝑡 ∈ 𝒯 , 𝑖, 𝑗 ∈ 𝒩 , 𝑖 < 𝑗, (33)
⎪ Constraints ∶ (1)–(13), (25b)–(25c), (e) ∑
⎩ 𝑧𝑡,𝑖,𝑗 ≥ 2𝑞𝑡,𝑗 , ∀𝑡 ∈ 𝒯 , 𝑗 ∈ 𝒩 . (34)
𝑖∈𝒱 ,𝑖≠𝑗
(27)
Constraint (31) means that the edge (0, 𝑖) is visited if and only if
where 𝜼′ = (𝜼, 𝝓, 𝐅, 𝝎, 𝝂, 𝜰 , 𝝅, 𝝂, ̂ 𝜰̂ , 𝝅),
̂ 𝝂 = (𝜈𝑖 )𝑖∈𝒩 , 𝜰 = (𝜰 𝑖 )𝑖∈𝒩 , 𝝅 = the 𝑖th retailer is visited in the 𝑡th period. Constraints (32)–(33) denote
(𝜋𝑖 )𝑖∈𝒩 , 𝝂̂ = (𝜈̂𝑖 )𝑖∈𝒩 , 𝜰̂ = (𝜰̂ 𝑖 )𝑖∈𝒩 , 𝝅̂ = (𝜋̂ 𝑖 )𝑖∈𝒩 . that the edge (𝑖, 𝑗) be traveled in the 𝑡th period if and only if both node
Additionally, Proposition 1 ensures that the optimal solution of 𝑖 and node 𝑗 are visited. Constraint (34) sets a lower bound on the sum
model (27) is that of model (25) with the ambiguity set (14). of the flows from node 𝑗.

10
Y. Feng et al. Omega 127 (2024) 103082

Model (27) is an MILP form when ‖𝜻 𝑖 ‖1 appears in the ambiguity

set (14). In contrast, if ‖𝜻 𝑖 ‖2 is used in the ambiguity set (14), then the
model (27) is an MISOCP form. Model (30) is a MILP form. In the next
section, we demonstrate the validity of these models using a case study.

6. Case study

This section conducts a case study to validate the feasibility of the

developed model. In Section 6.1, we give the description for the case
and data. In Section 6.2, we construct the ambiguity sets based on the
developed framework. We present the computational results of models
in Section 6.3. In Sections 6.4–6.5, we assess the optimal solutions
provided by models in out-of-sample performance and robustness. In
Section 6.6, we compare the proposed distributionally robust M-CVaR
IRP model with distributionally robust mean-semideviation (M-SD)
IRP model. Sensitivity analyses of key parameters are performed in
Sections 6.7–6.9. In Section 6.10, we further validate the feasibility of
the proposed method by applying it to instances of different scales. We
perform all experiments on an Inter(R) Core(TM) i5-12500H 2.50 GHz
computer with 16 GB RAM operating under Windows 10 (64bit). All
models are solved by the solver CPLEX 12.10.0.

6.1. Descriptions for problem and data

We consider a practical IRP case based on the one used by Soysal

et al. [25], which consists of one vendor and 20 retailers. The locations
of these nodes come from the Google map and are shown in Fig. 3. The
vendor (located in Torbali) is a tomato supplier. The retailers consist
of 11 physically existing supermarkets (Nodes 1-11) in different cities
in Turkey (Uria, Kusadasi, Didim, Aydin, Odemis, Salihli, Turgutlu,
Manisa, Soma, Dikili, Izmir) and 9 virtual supermarkets (Nodes 12-20).
The vendor delivers tomatoes to supermarkets. The demand of each
supermarket is uncertain. The vendor needs to make three decisions: (1)
the delivery timing and quantity from the vendor to each supermarket;
(2) the optimal number of vehicles to complete delivery; (3) the optimal Fig. 3. Locations of vendor and retailers.
delivery routes of vehicles.
Limited historical demands and the maximum inventory levels of
supermarkets, and the transportation costs between nodes are listed of historical demand data of supermarkets, which do not facilitate the
in supplementary material. The transportation costs are calculated by construction of demand scenarios and the calculation of the nominal
multiplying the transportation distance between nodes by the vehicle’s probability distribution. Considering this, we randomly generate 500
fuel consumption per unit distance and then multiplying by the unit historical demands in the interval [dmin max
𝑖 , d𝑖 ] to form a new historical
oil price. Without loss of generality, the initial inventory of each demand data set for supermarket 𝑖, 𝑖 ∈ 𝒩 . Supermarket 𝑖 faces dL𝑖 ,
supermarket is set to 20% of the demand in the first period. By referring dM H
𝑖 , or d𝑖 in each of the three periods. Based on the newly generated
to Soysal et al. [25], holding cost (𝑐𝑖ℎ , 𝑖 ∈ 𝒩 ) at supermarkets is set to demand data set, by conducting the first and second steps to construct
10% of the average marketplace selling price of tomatoes in that region the ambiguity set, we obtain the values of dL𝑖 , P(dL𝑖 ), dM M H
𝑖 , P(d𝑖 ), d𝑖 ,
and equals 0.06 ¤/kg. We first set the unit penalty cost (𝑐𝑖𝑝 , 𝑖 ∈ 𝒩 ) of H
and P(d𝑖 ) for supermarket 𝑖, 𝑖 ∈ 𝒩 . The results are summarized
product shortage equals to 0.08 ¤/kg to solve the model and present in Table 4. In the third step, we construct 27 demand scenarios for
the results. Then, in the comparison section of models, we will compare each supermarket. In the fourth step, we set the number of selected
the results provided by the models in three cases, 𝑐𝑖𝑝 = 0.08 ¤/kg, demand scenarios to 18. The final nominal value of probability vector
𝑐𝑖𝑝 = 0.09 ¤/kg, and 𝑐𝑖𝑝 = 0.12 ¤/kg. The vendor has 6 homogeneous for the selected demand scenarios after normalizing is presented in
vehicles for delivering products. The maximum load of any vehicle is supplementary material. We set 𝜒𝑖 = 0.2, 𝛤𝑖 = 1, and 𝜗𝑖 = 0.2 in the
5000kg. Moreover, we set 𝑇 = 3 in this study, i.e., considering the fifth step.
vendor to design order and routing solutions for three periods.
Considering that MILP is easier to solve than MISOCP under the 6.3. Computational results of models
same conditions, we select ‖𝜻 𝑖 ‖1 in the ambiguity set (14). We validate
that using ‖𝜻 𝑖 ‖2 is also feasible in (14) in Appendix D. Additionally, Without loss of generality, for all 𝑖 ∈ 𝒩 , 𝜒𝑖 is set to the same value.
owing to the fact that the geometry of ‖𝜻 𝑖 ‖1 is a polyhedron form, we 𝜗𝑖 , 𝛽𝑖 , and 𝜅𝑖 are similarly set up. For convenience, we denote 𝜒𝑖 , 𝜗𝑖 , 𝛽𝑖 ,
denote the model consisting of the model (27) and constraints (31)– and 𝜅𝑖 uniformly by 𝜒, 𝜗, 𝛽, and 𝜅, respectively.
(34) as Poly model. Because the geometry of the upper and lower bounds Based on 𝜒 = 0.2, 𝛽 = 0.5, and 𝜅 = 0.9, the results of the Poly
of 𝜻 𝑖 is a box form, we denote the model consisting of the model (30) model are as follows. The objective value provided by the Poly model
∑
and constraints (31)–(34) as Box model. is 5239.63, where TC1 = 1032.3. We denote the sum (i.e., 𝑖∈𝒩 𝐄[F𝑖 (𝜶 𝑖 ,
∑
𝜹𝑖 )]) of mean value, the sum (i.e., 𝑖∈𝒩 𝜙𝑖 ) of VaR value, and the sum
∑
6.2. Constructing ambiguity sets (i.e., 𝑖∈𝒩 𝐂𝐕𝐚𝐑[F𝑖 (𝜶 𝑖 , 𝜹𝑖 )]) of CVaR value of inventory and shortage
penalty costs as SOM, SOV, and SOCV, respectively. By calculating,
This part constructs the ambiguity sets (14)–(15) according to the SOM = 2946.9, and SOV = SOCV = 5467.8. The optimal delivery
proposed Algorithm 1 in Section 4. We have access to only four periods quantity provided by the Poly model in each period is shown in Table 5.

11
Y. Feng et al. Omega 127 (2024) 103082

Table 4
Low, medium, and high demands and their probabilities for each supermarket.
dL𝑖 P(dL𝑖 ) dM
𝑖
P(dM
𝑖
) dH
𝑖
P(dH
𝑖
) dL𝑖 P(dL𝑖 ) dM
𝑖
P(dM
𝑖
) dH
𝑖
P(dH
𝑖
)
R1 500 0.334 700 0.352 900 0.314 R11 2617 0.344 2849 0.332 3081 0.324
R2 1284 0.358 1450 0.32 1617 0.322 R12 467 0.364 600 0.354 734 0.282
R3 628 0.276 876 0.372 1125 0.352 R13 984 0.35 1350 0.34 1716 0.31
R4 750 0.338 1450 0.33 2150 0.332 R14 627 0.314 876 0.346 1125 0.34
R5 1000 0.338 1200 0.322 1400 0.34 R15 767 0.37 1499 0.314 2231 0.316
R6 568 0.292 901 0.356 1234 0.352 R16 817 0.336 1050 0.322 1284 0.342
R7 550 0.318 650 0.342 750 0.34 R17 566 0.358 897 0.31 1229 0.332
R8 569 0.316 1101 0.34 1634 0.344 R18 550 0.372 650 0.328 750 0.3
R9 550 0.334 850 0.348 1150 0.318 R19 500 0.32 900 0.362 1300 0.318
R10 516 0.332 948 0.316 1380 0.352 R20 517 0.322 749 0.322 981 0.356

Note: R𝑖,𝑖∈{1,…,20} denotes the 𝑖th supermarket.

The optimal number of vehicles to complete the product delivery and ambiguities in demands distributions. Under 𝛽 = 0.5 𝜅 = 0.9, 𝜒 = 0.2
the optimal delivery routes of vehicles in each period are shown in and 𝜗 = 0.2, the PDRs of the Poly and Box models are 739.36 and
Fig. 5. From Fig. 5, the vendor uses five vehicles in the first period 193.22, respectively. We find that the above Box model is superior to
and three vehicles in the second and third periods to complete product the Poly model with respect to PDR.
delivery. The optimal delivery route for each vehicle is shown in Fig. 5.
Based on 𝜗 = 0.2, 𝛽 = 0.5, and 𝜅 = 0.9, the Box model results are as
6.4. Out-of-sample performance
follows. The objective value provided by the Box model is 4693.49,
where TC1 = 991.27, SOM = 2529.3, SOV = 3418.14, and SOCV =
4875.1. Regardless TC1 , SOM or SOCV, the Box model is lower than This part assesses the optimal delivery quantities provided by the
the Poly model with 𝜒 = 0.2, 𝛽 = 0.5, 𝜅 = 0.9. This is also observed in Poly model, Box model, and nominal model in terms of out-of-sample
Fig. 4. The optimal delivery quantity, the optimal number of vehicles performance. First, for 𝛽 = 0.5, 𝜅 = 0.9, 𝜒 = 0.2, and 𝜗 = 0.2, we
to complete delivery of products and the optimal delivery routes of solve the above three models based on 𝑐𝑖𝑝 = 0.08, 𝑐𝑖𝑝 = 0.09, and
vehicles in each period are shown in Table 5 and Fig. 5. We observe 𝑐𝑖𝑝 = 0.12, respectively. Subsequently, we randomly generate an out-
of-sample demand set based on [dmin max
that the optimal delivery quantities determined by two models differed 𝑖 (1 − 𝛥), d𝑖 (1 + 𝛥)], 𝑖 ∈ 𝒩 , where
significantly. For example, in the second period, the Poly model chooses 𝛥 regulates the range of generated out-of-sample demand. We use 𝐿
not to deliver to the first supermarket, wherea the Box model chooses to to denote the size of this set and let 𝐿 = 3000. A larger 𝛥 implies
deliver. For more details on the differences between these two models greater perturbation in the set of generated out-sample demands. Fi-
in terms of delivery quantity, refer to Table 5. The optimal number nally, for the optimal delivery quantity 𝒒 ∗ obtained from each model,
of vehicles required to complete the product delivery and the optimal we compare the performance of the following three indicators: (i) the
delivery routes of vehicles in each period provided by these two models average value (denoted as Aver) of the sum of inventory and penalty
are also different. In the Box model, the vendor uses four (instead of costs for each sample, (ii) the worst value (denoted as Wor) of sum of
five) vehicles to complete the product delivery in the first period. The inventory and penalty costs about each sample, and (iii) the shortage
optimal delivery route for each vehicle is shown in Fig. 5. rate (denoted as SR). These three indicators are defined as follows:
Based on 𝛽 = 0.5 and 𝜅 = 0.9, the results of the nominal model ∑ ∑ ∑ { }
ℎ 𝑠 𝑝 𝑠
(the stochastic M-CVaR IRP model) are as follows. The nominal model, 𝑠∈ℒ 𝑡∈𝒯 𝑖∈𝒩 max 𝑐𝑖 min{I𝑡,𝑖 , 𝐶𝑖 }, −𝑐𝑖 I𝑡,𝑖
Aver = ,
which uses the nominal probability values of demand scenarios and { 𝐿 }
neglects the ambiguities of these probabilities, usually generates weak ∑∑ { }
Wor = max max 𝑐𝑖ℎ min{I𝑠𝑡,𝑖 , 𝐶𝑖 }, −𝑐𝑖𝑝 I𝑠𝑡,𝑖 ,
robust delivery solution. The objective value determined by the nom- 𝑠∈ℒ
𝑡∈𝒯 𝑖∈𝒩
inal model is 4500.27, where TC1 = 994.02, SOM = 2366.8, SOV = ∑ ∑ ∑
𝑠∈ℒ 𝑡∈𝒯 𝑖∈𝒩 1(I𝑠𝑡,𝑖 < 0)
3352.55, and SOCV = 4645.7. From Fig. 4, we observe that SOM, SOV, SR = ,
SOCV, and the optimal objective values determined by the above Poly 𝐿×𝑇 ×𝑁
∗ − d𝑠 , I𝑠 = I , ∀𝑠 ∈ ℒ , 𝑡 ∈ 𝒯 , 𝑖 ∈ 𝒩 . ℒ denotes
where I𝑠𝑡,𝑖 = I𝑠𝑡−1,𝑖 + 𝑞𝑡,𝑖
and Box models are all higher than those provided by the nominal 𝑡,𝑖 0,𝑖 0,𝑖
model. Table 5 and Fig. 5 present the detailed order and routing the index set of out-of-sample demand (i.e., ℒ = {1, … , 𝐿}). 1(event) is
solutions. Through Table 5 and Fig. 5, we can see that the optimal order defined as the indicator function of the event; 1(event) = 1 if the event
and routing solutions in each period determined by the nominal model is true and = 0 otherwise.
are not the same as those of corresponding Poly and Box models. We For convenience, we develop Algorithm 2 to perform the above out-
conclude that the ambiguities in probabilities significantly affect the of-sample performance tests. The test results of the optimal delivery
optimal delivery solution. quantities provided by these three models are listed in Table 6.
We explain the difference between the distributionally robust M- From Table 6, we find that the performance of the optimal delivery
CVaR IRP model (i.e., the Poly and Box models) and its nominal model quantity provided by our distributionally robust M-CVaR IRP model
in terms of the optimal objective. Based on the same 𝛽 and 𝜅, it is (Poly model or Box model) outperforms that provided by the nominal
reasonable that the optimal objective value determined by the nominal model either Aver or Wor in most cases of (𝑐𝑖𝑝 , 𝛥). The Aver or Wor
model is lower than those of the Poly and Box models. The nominal performance of the delivery quantity offered by the nominal model
model makes decisions based on deterministic distributions of demands performs better in a few other cases of (𝑐𝑖𝑝 , 𝛥), but not much better
and ignores ambiguities. However, the Poly and Box models consider than that of our model. Notably, the Poly model provides the delivery
ambiguities in demands distributions and make decisions based on lim- quantity with the smallest SR in all cases. Additionally, the performance
ited distributionally information. The price of distributional robustness of the optimal delivery quantity given by the Box model is also better
(denoted as PDR) is introduced to further elaborate on the difference of than that given by the nominal model in the indicator SR. With fixed
the optimal objective. Here, the PDR equals the objective value of the 𝛥, the indicator SR for the optimal delivery quantities provided by
Poly or Box model minus the objective value of the nominal model. The the three models decreases as the unit penalty cost 𝑐𝑖𝑝 increases. This
indicator PDR implies the price paid by the optimal delivery solution suggests that the decision-maker can reduce the shortage risk of the IRP
determined by the Poly or Box model to immunize the effect of the system by tuning the value of 𝑐𝑖𝑝 .

12
Y. Feng et al. Omega 127 (2024) 103082

Fig. 4. TC1 , SOM, SOCV, Objective determined by the models with 𝛽 = 0.5 and 𝜅 = 0.9.

Table 5
Delivery quantity for each supermarket given by the Poly, Box, and nominal models.
Supermarket
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
∙ Poly model with 𝜒 = 0.2, 𝛽 = 0.5, and 𝜅 = 0.9. Total
𝑡=1 920 1252 1028 1230 1145 869 702 1170 848 996 1114 630 1186 1207 1307 978 969 1010 700 640 19 901
𝑡=2 414 1121 876 750 1147 891 840 1009 792 516 467 643 1126 876 767 1000 703 550 500 0 14 988
𝑡=3 586 876 876 750 901 750 748 991 550 516 600 567 876 876 767 939 566 407 900 0 14 042
∙ Box model with 𝜗 = 0.2, 𝛽 = 0.5, and 𝜅 = 0.9. Total
𝑡=1 920 1170 988 1230 1145 869 590 1170 892 996 1114 460 1110 1005 1307 875 737 1010 700 640 18 928
𝑡=2 0 1203 913 750 1147 891 750 1009 705 516 467 601 1202 876 767 1182 798 0 500 0 14 277
𝑡=3 1000 876 879 750 901 750 599 991 594 516 600 524 876 983 767 884 703 800 900 0 14 899
∙ Nominal model with 𝛽 = 0.5 and 𝜅 = 0.9. Total
𝑡=1 920 1170 813 1230 1190 869 590 1170 690 996 1114 460 1110 1005 1307 1019 737 1010 400 640 18 440
𝑡=2 0 1203 1088 750 1057 891 650 1009 851 516 467 600 1202 876 767 927 798 0 800 0 14 452
𝑡=3 1000 876 876 750 946 750 650 991 649 516 600 468 876 876 767 995 703 800 900 0 14 989

Algorithm 2: Calculating out-of-sample performance of 𝒒 ∗ . the advantage of the optimal delivery quantity provided by our dis-
max ℎ 𝑝 tributionally robust M-CVaR IRP model. From Fig. 6, we observe that
Input: {dmin
𝑖 , d𝑖 , 𝑐𝑖 , 𝑐𝑖 , 𝐶𝑖 }𝑖∈𝒩 , 𝒒∗ ,
𝐿, 𝛥;
the frequency distribution histograms and fitted profiles of the perfor-
Output: Aver, Wor, SR;
mances of the optimal delivery quantities provided by the Poly and
Generating two empty (𝐿 × 𝑇 × 𝑁)-dimensional matrices: D and I;
Box models are to the left of those given by the nominal model. This
Generating a empty 𝐿-dimensional array: TCP;
directly demonstrates the superiority of our distributionally robust M-
Generating a int variable count = 0;
CVaR IRP model, which provides a promising delivery quantity with
for 𝑠 = 1 to 𝐿 do
excellent out-of-sample performance.
for 𝑡 = 1 to 𝑇 do
for 𝑖 = 1 to 𝑁 do
D𝑠𝑡,𝑖 ← randomly generate value from interval
max
6.5. Feasibility (Robustness) analysis of delivery solution
[dmin
𝑖 (1 − 𝛥), d𝑖 (1 + 𝛥)];
𝑠 𝑠 ∗
I𝑡,𝑖 = I𝑡−1,𝑖 + 𝑞𝑡,𝑖 − D𝑠𝑡,𝑖 ;
if I𝑠𝑡,𝑖 < 0 then We now evaluate the feasibility (robustness) of the optimal delivery
count = count +1; solution provided by the Poly and Box models. Given a set of out-of-
{ } sample demands, a delivery solution is feasible if the cost incurred
∑ ∑
TCP(𝑠)= 𝑡∈𝒯 𝑖∈𝒩 max 𝑐𝑖ℎ min{I𝑠𝑡,𝑖 , 𝐶𝑖 }, −𝑐𝑖𝑝 I𝑠𝑡,𝑖 ; by the delivery solution is less than the expected cost provided by
the model. Given 𝐿 demand samples, we define the metric feasibility
Aver ← Calculating the average value of TCP;
probability (FP) as the ratio of the number of delivery solution is feasible
Wor ← Calculating the maximum value of TCP;
count to 𝐿, which is calculated using the following formula:
SR ← Calculating the rate 𝐿×𝑇 ;
×𝑁
∑𝐿 ( ∑ ∑ { } )
𝑝
𝑠=1 1 TC∗1 + 𝑡∈𝒯 𝑖∈𝒩 max 𝑐𝑖ℎ min{I𝑠𝑡,𝑖 , 𝐶𝑖 }, −𝑐𝑖 I𝑠𝑡,𝑖 ≤ TC∗1 + SOM
FP =
(∑ { 𝐿 } )
∑𝐿 ∑
1 max 𝑐𝑖ℎ min{I𝑠𝑡,𝑖 , 𝐶𝑖 }, −𝑐𝑖𝑝 I𝑠𝑡,𝑖 ≤ SOM
Based on the cases (𝑐𝑖𝑝 , 𝛥)
= (0.08, 0.4) and (0.09, 0.4), we draw a =
𝑠=1 𝑡∈𝒯 𝑖∈𝒩
,
frequency distribution histogram (see Fig. 6) to further demonstrate 𝐿

13
Y. Feng et al. Omega 127 (2024) 103082

Fig. 5. Delivery routes provided by models under 𝛽 = 0.5 and 𝜅 = 0.9.

Table 6
Test results of out-of-sample performance of the optimal delivery quantity.
(𝑐𝑖𝑝 , 𝛥) Model (𝑐𝑖𝑝 , 𝛥) Model (𝑐𝑖𝑝 , 𝛥) Model
Poly Box Nominal Poly Box Nominal Poly Box Nominal
(0.08, 0) Aver 2694.17 2625.51 2632.31 (0.09, 0) 2978.92 2887.64 2887.37 (0.12, 0) 3745.37 3610.36 3598.06
Worst 3671.72 3627.1 3652.52 3896.91 3819.06 3857.25 5141.3 5034.2 5055.5
SR 45.16% 52.53% 54.41% 43.41% 51.14% 53.08% 42.65% 47.59% 48.35%
(0.08, 0.1) Aver 2999.77 2973.19 2990.93 (0.09, 0.1) 3337.32 3294.12 3307.13 (0.12, 0.1) 4228.11 4147.03 4144.16
Worst 4122.62 4075.62 4130.46 4680.39 4797.63 4876.83 6062.06 6025.62 6034.14
SR 51.05% 58.05% 59.70% 49.61% 57.03% 58.77% 48.47% 54.59% 55.49%
(0.08, 0.2) Aver 3360.83 3370.04 3395.45 (0.09, 0.2) 3726.2 3723.98 3748.44 (0.12, 0.2) 4755.62 4721.42 4726.50
Worst 4863.46 4912.28 4944.14 5221.35 5343.87 5411.91 6872.66 6933 6989.22
SR 56.24% 62.42% 63.85% 54.37% 61.31% 62.79% 53.52% 59.52% 60.37%
(0.08, 0.3) Aver 3751.29 3788.46 3818.99 (0.09, 0.3) 4160.38 4192.05 4225.34 (0.12, 0.3) 5320.56 5325.01 5336.77
Worst 5266.68 5305.94 5332.22 5919.54 5968.68 6017.73 7716.54 7918.08 7944.42
SR 59.75% 65.19% 66.40% 58.16% 64.58% 65.94% 56.90% 62.64% 63.43%
(0.08, 0.4) Aver 4178.33 4235.99 4271.26 (0.09, 0.4) 4632.11 4690.9 4729.15 (0.12, 0.4) 5963.13 6007.15 6024.05
Worst 6077.24 6214.62 6296.9 6599.19 6687.72 6789.87 8814.98 8987.34 9050.22
SR 62.90% 67.67% 68.75% 61% 66.65% 67.86% 60.33% 65.69% 66.39%
(0.08, 0.5) Aver 4611.82 4689.37 4727.55 (0.09, 0.5) 5065.51 5145.08 5189.43 (0.12, 0.5) 9890.48 9965.18 10029.68
Worst 6723.4 6845.96 6931 7888.23 8075.34 8146.26 6603.41 6672.13 6693.07
SR 65.05% 69.22% 70.15% 63.35% 68.34% 69.39% 62.71% 67.51% 68.11%

where TC∗1 denotes the transportation cost provided by the Poly or Box Poly model increases with 𝛤 . The same phenomenon also occurs for
model. The second equation holds because transportation cost do not parameter 𝜎. This is due to the fact that as 𝛤 and 𝜎 become larger, the
change. scale of the ambiguity set (14) also becomes larger. In addition, the
The value of the indicator FP reflects the robustness of the delivery Poly model is superior to the nominal model for all 𝛤 . With Fig. 7, we
solution. A larger FP implies that the delivery solution is highly feasible observe that the metric FP of the optimal delivery solution provided
in the out-of-sample scenario, indicating that the delivery solution is by the Box model increases with 𝜗 (all of these cases are better than
highly robust. For the Poly model with 𝛽 = 0.5 and 𝜅 = 0.9, we set 𝛤𝑖 , the nominal model). This is because the scale of the ambiguity set (15)
𝑖 ∈ 𝒩 , to be equal and denote 𝛤𝑖 by 𝛤 . We set 𝛤 to take values from 0.4 increases with 𝜗. The experimental results verify that the IRP decision-
to 1.6 with an interval of 0.2. Additionally, we set 𝜎 = 𝜒 to 0, 0.2, 0.4, maker can adjust the robustness of the delivery decision by adjusting
and 0.6. The distributionally robust M-CVaR IRP model degenerates the corresponding scale parameters in Remarks 3–4.
into the stochastic M-CVaR IRP model if 𝜎 = 0. In this setup, we obtain
28 different (𝜎, 𝛤 ). Regarding the Box model with 𝛽 = 0.5 and 𝜅 = 0.9, 6.6. Comparing with distributionally robust mean-semideviation IRP model
and we set 𝜗 to take values from 0.2 to 1.8 with an interval of 0.05.
Based on the out-of-sample set with 𝐿 = 3000 and 𝛥 = 0, the test results In this part, we compare the distributionally robust M-CVaR IRP
of the indicator FP of the optimal delivery solutions provided by these model with the distributionally robust mean-semideviation (M-SD) IRP
two models are listed in Table 7 and Fig. 7. model to provide guidelines for the decision-maker to manage the
From Table 7, with the parameter 𝜎 fixed, FP increases with 𝛤 , risks of the IRP system. Semideviation is also a well-known risk-averse
indicating that the robustness of the delivery solution provided by the measure that only quantifies the positive deviation between cost and

14
Y. Feng et al. Omega 127 (2024) 103082

Fig. 6. Frequency distribution histogram of performance of the delivery quantities given by the above models.

Table 7
Feasibility probability of the delivery solution provided by the Poly model.
𝛤 = 0.4 𝛤 = 0.6 𝛤 = 0.8 𝛤 =1 𝛤 = 1.2 𝛤 = 1.4 𝛤 = 1.6
𝜎 =0 14.63% 14.63% 14.63% 14.63% 14.63% 14.63% 14.63%
𝜎 = 0.2 66.27% 66.3% 75.67% 84.73% 90.57% 93.7% 96.57%
𝜎 = 0.4 75.7% 90.7% 96.57% 98.93% 99.77% 99.97% 100%
𝜎 = 0.6 90.37% 98.03% 99.77% 100% 100% 100% 100%

expected cost and does not consider downside risk. The distributionally re-expressed as:
robust M-SD IRP model is established as follows: ∑ { }
⎧ min TC1 (𝒙, 𝒒, 𝒛, 𝒚) + max 𝛽𝑖 𝐅𝖳𝑖 𝐏𝑖 + (1 − 𝛽𝑖 )𝓵 𝖳𝑖 𝐏𝑖
[ ] ⎪𝜼,𝐅,𝓵 𝑖∈𝒩
𝐏𝑖 ∈𝒫𝑖
⎧ ∑ { [ ] [ ]} ⎪
⎪min TC1 (𝒙, 𝒒, 𝒛, 𝒚) + max 𝛽𝑖 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) + (1 − 𝛽𝑖 )𝐒𝐃𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 )
𝐏𝑖 ∈𝒫𝑖
⎪ s.t. 𝓁𝑖𝑠 ≥ F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) − min 𝐅𝖳𝑖 𝐏𝑖 , ∀𝑠 ∈ 𝒮 , 𝑖 ∈ 𝒩 ,
⎨ 𝜼,𝐅 𝑖∈𝒩 ⎨ 𝐏𝑖 ∈𝒫𝑖
⎪ ⎪
⎩ s.t. Constraints ∶ (1)–(13), ⎪ 𝓁𝑖𝑠 ≥ 0, ∀𝑠 ∈ 𝒮 , 𝑖 ∈ 𝒩 ,
(35) ⎪ Constraints ∶ (1)–(13),
⎩
[ ] [ { [ ] }] where 𝓵 = (𝓵 𝑖 )𝑖∈𝒩 , and 𝓵 𝑖 = (𝓁𝑖𝑠 )𝑠∈𝒮 .
where 𝐒𝐃𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) = 𝐄𝐏𝑖 max F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) − 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) , 0 .
Based on the developed ambiguity sets (14)–(15), we derive the
For given 𝑖 and 𝑠, we introduce auxiliary variable 𝓁𝑖𝑠 to represent tractable reformulation form of model (35) in Appendix E. We calculate
[ ]
max{F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) − 𝐄𝐏𝑖 F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) , 0}. Then, the aforementioned model is these two tractable reformulations based on 𝛽 = 0.5, 𝑆 = 18, and 𝑐𝑖𝑝 =

15
Y. Feng et al. Omega 127 (2024) 103082

types of models is close based on the same ambiguity set and

parameter 𝛥. With respect to the Worst and SR performances,
the distributionally robust M-CVaR IRP model outperforms the
distributionally robust M-SD IRP model. This advantage is more
pronounced, especially when parameter 𝛥 is large. Through
Fig. 9, we see that the frequency distribution histogram and
fitted profile of the performance of the optimal delivery quantity
provided by the distributionally robust M-CVaR IRP model is to
the left of that given by the distributionally robust M-SD IRP
model. This directly substantiates the superiority of our model.

6.7. Sensitivity analysis of 𝜒 and 𝜗

The decision-maker can regulate the scales of ambiguity sets via

adjusting the scale parameters 𝜒 and 𝜗. To investigate the influence
of the changes in the above scale parameters on the objective value
offered by models, we execute the sensitivity analyses on 𝜒 and 𝜗. We
set 𝛽 = 0.5, 𝜅 = 0.9, and 𝛩 = 𝜒 = 𝜗 to take values from 0 to 0.7 with
Fig. 7. Feasibility probability of the delivery solution provided by the Box model. an interval of 0.05. If 𝛩 gradually increases, then the ambiguity sets
(14)–(15) will also increase in scale. Fig. 10 shows the results.
In Fig. 10, we observe that the objective value increases with 𝛩 for
0.08. The optimal delivery quantity and route provided by the distribu- both the Poly and Box models. This is mainly because as 𝛩 gradually
tionally robust M-SD IRP model are presented in Table 8 and Fig. 8. The increases, the Poly and Box models will pay more price to resist the
test results for the out-of-sample performance of the optimal delivery influence of the ambiguities in demands distributions. In addition, the
quantity is presented in Table 9. We draw the frequency distribution following key finding is captured: for either the Poly or Box model,
histogram (see Fig. 9) of the performance of the delivery quantities their optimal objective values are consistent with the objective value
given by distributionally robust M-CVaR and M-SD IRP models under determined by the nominal model when 𝛩 = 0. This verifies that the
the out-of-sample set with 𝛥 = 0.4. Detailed comparisons of these two ambiguity sets contain only the nominal distribution and the distri-
distributionally robust models based on different risk-aversion criteria butionally robust M-CVaR IRP model degenerates into the stochastic
are presented below. M-CVaR IRP model when 𝛩 = 0.

(1) The comparison of budget information. We denote the sum

∑ 6.8. Sensitivity analysis of 𝑆
(i.e., 𝑖∈𝒩 𝐒𝐃[F𝑖 (𝜶 𝑖 , 𝜹𝑖 )]) of SD of inventory and shortage penalty
costs as SOSD. The objective value provided by the distribution-
In this part, to investigate the influence of the changes in the
ally robust M-SD IRP model with ambiguity set (14) is 2823.07,
number (i.e., 𝑆) of selected demand scenarios on the objective value
where TC1 = 991.27, SOM = 2929.64, and SOSD = 733.94.
and PDR, we execute the sensitivity analyses on 𝑆. We set 𝛽 = 0.5,
The objective value provided by this model with ambiguity
𝜅 = 0.9, 𝜒 = 0.2, 𝜗 = 0.2, and 𝑆 to take values from 12 to 27 with an
set (15) is 2492.08, where TC1 = 968.8, SOM = 2555.32,
interval of 3. In particular, 𝑆 = 27 means that all demand scenarios are
and SOSD = 491.25. Based on the same ambiguity set, the
selected. Fig. 11 reports the results.
values of these two different risk-averse types of models with
In Fig. 11, we observe that the optimal objective values increase
respect to TC1 and SOM are close. The distributionally robust
with 𝑆 for both the Poly and Box models. The same phenomenon occurs
M-SD IRP model provides the expected cost (i.e., SOM) and
with PDRs. We explain the phenomenon that the optimal objective
the information (i.e., SOCD) that the cost exceeds SOM. The
values increase as 𝑆. A larger 𝑆 implies that the optimal delivery
constructed distributionally robust M-CVaR IRP model provides
solution generated by the models can adapt more demand scenarios.
SOM, the maximum possible cost (i.e., SOV) at a given probabil-
In addition, the dimension of probability vector becomes larger when
ity level (downside risk parameter), and also provides informa-
𝑆 is larger. In this case, the Poly and Box models will pay more price
tion (i.e., SOCV) that the cost exceeds SOV. The distributionally
to immunize the influence of ambiguous distributions. Therefore, it is
robust M-CVaR IRP model is more comprehensive.
not surprising that the optimal objective values and PDRs increase as
(2) The comparison of the delivery solution. There is a clear dif-
𝑆 for both the Poly and Box models.
ference between the optimal delivery solutions proposed by
these two different risk-averse types of models. In terms of the
delivery quantity, by comparing Tables 5 and 8, we observe 6.9. Sensitivity analyses of 𝛽 and 𝜅
that the distributionally robust M-CVaR IRP model prefers to
deliver more products (19901 > 18834 and 18928 > 15000) to In this part, we conduct sensitivity analyses of 𝛽 and 𝜅 to explore
the retailers in the initial stage to build up a relatively high the effects of changes in these two parameters on the objective and PDR
inventory. However, the distributionally robust M-SD IRP model of the Poly and Box models. We set 𝜒 = 0.2 for the Poly model and
is more willing to take more replenishment actions (14042 < 𝜗 = 0.2 for the Box model. We set 𝜅 to take values in [0.65, 0.9] with
14888 and 14277 < 17635) to retailers later in the sale. There is an interval of 0.05 and set 𝛽 to take values [0.3, 0.8] with an interval of
also a more pronounced difference between these two different 0.1. Figs. 12–13 present the results.
risk-averse types of models in both the number and delivery We first elaborate on the changes in the objective value and PDR
routes of vehicles. This can be captured by comparing Figs. 5 of the Poly model. In Fig. 12(a), we can observe that the objective
and 8. decreases as 𝛽 increases if 𝜅 is fixed. The objective value increases as 𝜅
(3) The comparison of out-of-sample-performance. By comparing and finally becomes stable if 𝛽 is fixed. According to Fig. 12(b), the PDR
Tables 6 and 9, we observe that the performance Aver of the decreases as 𝛽 increases when 𝜅 is constant. In addition, by fixing 𝛽, the
delivery quantities provided by these two different risk-averse PDR first increases and then decreases with 𝜅. In addition, a similar

16
Y. Feng et al. Omega 127 (2024) 103082

Table 8
Delivery quantity provided by the distributionally robust M-SD IRP model under different ambiguity sets.
Supermarket
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
∙ Distributionally robust M-SD IRP model with ambiguity set (14). Total
𝑡=1 920 1252 944 1230 1145 869 590 1170 848 996 1114 522 1110 893 1307 874 700 1010 700 640 18834
𝑡=2 0 1121 960 750 1147 891 650 1009 792 516 467 514 1202 1005 767 1184 835 0 500 0 14310
𝑡=3 1000 876 876 750 901 750 664 991 550 516 600 567 876 918 767 883 703 800 900 0 14888
∙ Distributionally robust M-SD IRP model with ambiguity set (15). Total
𝑡=1 0 1099 823 1230 1019 661 490 1170 690 996 1114 460 1044 756 960 913 524 411 0 640 15000
𝑡=2 1420 1274 898 750 1171 1086 702 740 764 516 467 520 1259 927 1114 1140 938 749 1200 0 17635
𝑡=3 0 876 1059 750 1003 763 717 1260 736 516 600 608 885 1130 767 889 775 579 900 0 14813

Table 9
Out-of-sample performance of delivery quantity provided by the distributionally robust M-SD IRP model.
Ambiguity set Ambiguity set Ambiguity set
(14) (15) (14) (15) (14) (15)
𝛥=0 Aver 2610.67 2669.98 𝛥 = 0.2 3369.84 3485.47 𝛥 = 0.4 4218.79 4373.17
Worst 3699.84 3835.32 4976.76 5195.52 6341.3 6636.9
SR 52.95% 61.48% 62.61% 69.05% 67.78% 73%
𝛥 = 0.1 Aver 2971.5 3061.7 𝛥 = 0.3 3784.26 3924.57 𝛥 = 0.5 4656.3 4821.36
Worst 4286.14 4493.82 5594.8 5898.14 7031.16 7243.6
SR 58.62% 65.92% 65.82% 71.58% 69.2% 73.97%

Fig. 8. Delivery routes provided by distributionally robust M-SD IRP model with different ambiguity sets.

Fig. 9. Frequency distribution histogram of performance of the delivery quantities given by different distributionally robust risk-averse IRP model under out-of-sample set with
𝛥 = 0.4.

17
Y. Feng et al. Omega 127 (2024) 103082

the specified solution time limit for small- and medium-scale instances,
e.g., the instances with (𝑁, 𝑇 ) = (20, 3), (𝑁, 𝑇 ) = (25, 2), (𝑁, 𝑇 ) =
(30, 2). For large-scale instances, e.g., instances with (𝑁, 𝑇 ) = (20, 4),
(𝑁, 𝑇 ) = (30, 3), (𝑁, 𝑇 ) = (35, 2), (𝑁, 𝑇 ) = (35, 3), the optimal solutions
of these two models cannot be found in most cases within the specified
solving time limit. However, the solution gaps are significantly less than
1% when the optimal solution is not found. We consider these gaps
acceptable to IRP decision-maker. The settings of 𝑆 and 𝑐𝑖𝑝 affect the
solution time and gap. In most cases, the solution time and gap increase
with 𝑆. Regarding the effect of 𝑐𝑖𝑝 on solving the model, a larger 𝑐𝑖𝑝
helps to solve the model in some cases. However, there are also cases
in which a smaller 𝑐𝑖𝑝 helps to solve the model.
For large-scale instances related to 𝑁 and 𝑇 , the main reason the op-
timal solution cannot be found within a reasonable time is the increase
in binary variables in the models. These binary variables are self-
contained in the IRP model and are not introduced when transforming
the proposed distributionally robust M-CVaR IRP model into a tractable
reformulation. Meanwhile, we are convinced that the following two
different research focuses are worth investigating: first, exploring how
to balance the expectation and risk level of the cost of IRP system based
on the ambiguous distributions of uncertain demands under downside
risk, and second, exploring how to design efficient algorithms to solve
large-scale instances.

7. Management insights
Fig. 10. The objective values under different 𝛩.

In light of the above experimental results and analyses, we summa-

rize the following management insights as a guide for IRP managers.
phenomenon for 𝜅 and 𝛽 also occurs with the Box model. This can be
First, demands for retailers are usually uncertain. In this context, it
verified by observing Fig. 13.
may be inappropriate for decision-maker to use a model with determin-
We also explore the values of SOM and SOCV determined by the
istic demands to make the IRP decision. In addition, the experimental
Poly and Box models based on different 𝛽 and 𝜅. By fixing 𝜒 = 0.2
results suggest a remarkable difference between the optimal solutions
and 𝜗 = 0.2, we set 𝛽 to take values in [0.1, 0.9] with an interval of
provided by our distributionally robust M-CVaR IRP model and the
0.05 and set 𝜅 to take values [0.7, 0.9] with an interval of 0.1. Fig. 14
nominal stochastic model. The decision-maker should not neglect the
depicts the variations of SOM and SOCV values. In both the Poly and
effect of distributions ambiguities on the optimal delivery solution and
Box models, the SOCV value increases with 𝛽, whereas the SOM value
can apply our distributionally robust M-CVaR IRP model to obtain the
decreases with 𝛽. Moreover, we also observe that the SOM value for a
optimal robust solution. The proposed distributionally robust M-CVaR
larger 𝜅 is always overtops that of a smaller 𝜅. The SOCV value follows
IRP model significantly upgrades the existing literature.
the same trend.
Second, sensitivity analyses of the scale parameters (i.e., 𝜒, 𝛤 , and
6.10. Test experiment in instances of different scales 𝜗) of the ambiguity sets show that these three parameters significantly
affect the optimal solution provided by the corresponding model. Based
The above experiments verify the feasibility of our method for the on larger-scale parameters, although the optimal solutions provided by
case with 20 retailers and 3 periods. In this part, other instances of the Poly or Box model can resist the influence of a larger perturbation
different scales are used to test the feasibility of the proposed method. of distribution, the corresponding cost increases. If the IRP decision-
First, based on the aforementioned case in Section 6.1, we test the maker is conservative, then he/she should set larger-scale parameters.
feasibility for the case with 20 retailers and 4 periods under different Otherwise, they can choose smaller-scale parameters.
𝑆. Then, under different 𝑆, we test the feasibility in cases where the Third, it is critical to set the number (i.e., 𝑆) of selected scenarios.
number of retailers is set to 25, 30, 35, and the number of periods is The sensitivity analysis shows that 𝑆 affects the optimal solution. At a
set to 2 and 3. The data for the first case are the same as those used larger 𝑆, the optimal solution has higher robustness (i.e., the optimal
in Section 6.1. The original data for the remaining cases are generated solution is feasible for many demand scenarios), but the corresponding
with reference to Coelho & Laporte [63], who assume that the inventory cost increases. The decision-maker should set 𝑆 according to the risk
at each retailer must be greater than or equal to 0 at the end of each preference. If the IRP decision-maker is conservative, he/she should set
period and do not consider product shortages. Product shortages are a larger 𝑆. Otherwise, they could select a smaller 𝑆.
permitted in this study. The following 𝑐𝑖𝑝 values are used in all cases. Fourth, parameters 𝛽 and 𝜅 affect the optimal solution. A smaller 𝛽
First, 𝑐𝑖𝑝 is set to be the same as 𝑐𝑖ℎ , following the idea of Li et al. [57]. value indicates that the risk level of the sum of inventory and shortage
Additionally, we set 𝑐𝑖𝑝 = 2𝑐𝑖ℎ and 𝑐𝑖𝑝 = 5𝑐𝑖ℎ to consider the case where penalty costs for each retailer is emphasized. A higher 𝜅 value indicates
𝑐𝑖𝑝 > 𝑐𝑖ℎ . Further explanations of the parameter value settings for the that a higher degree of risk aversion of the IRP decision-maker. When
latter cases are provided in Appendix F. The ambiguity sets and demand the decision-maker focuses on risk, a smaller 𝛽 and larger 𝜅 should be
scenarios of for cases are constructed using Algorithm 1. We set 𝜒 = 0.2, set. Otherwise, a larger 𝛽 and smaller 𝜅 could be set.
𝜗 = 0.2, 𝛽 = 0.5, and 𝜅 = 0.9. The solution time for each model is limited Fifth, although this study is devoted to classical IRP under uncertain
to 14400 seconds. The computational results of the Poly and Box models demands with ambiguous distributions and risk-averse, the proposed
for instances with different scales are listed in Table 10. methods can also be applied to other IRP variants with appropriate
From Table 10, we observe the following results. The optimal adjustments. If the decision-maker possesses limited historical demand
solutions of the Poly and Box models can be found in most cases within data of retailers, the proposed framework for constructing demand

18
Y. Feng et al. Omega 127 (2024) 103082

Fig. 11. The objective values and PDRs provided by the Poly and Box models under different 𝑆.

Fig. 12. The objective value and PDR provided by the Poly model under different 𝛽 and 𝜅.

Fig. 13. The objective value and PDR provided by the Box model under different 𝛽 and 𝜅.

19
Y. Feng et al. Omega 127 (2024) 103082

Fig. 14. SOM and SOCV determined by the Poly and Box models with 𝜅 = 0.9 under different 𝛽.

Table 10
The computational results of the Poly and Box models under instances with different scales.
𝑁-𝑇 𝑆 𝑐𝑖𝑝 = 𝑐𝑖ℎ 𝑐𝑖𝑝 = 2𝑐𝑖ℎ 𝑐𝑖𝑝 = 5𝑐𝑖ℎ
Ploy model Box model Ploy model Box model Ploy model Box model
UB Gap Time(s) UB Gap Time(s) UB Gap Time(s) UB Gap Time(s) UB Gap Time(s) UB Gap Time(s)
20-3 12 2554.78 0.01% 169 2508.11 0.01% 226 3433.31 0.01% 5634 3357.17 0.01% 1312 5849.69 0.05% 14400 5665.47 0.01% 513
18 4308.87 0.01% 75 3912.15 0.01% 1005 7080.18 0.01% 49 6252.67 0.01% 20 15512.25 0.01% 73 13247.16 0.01% 11
24 4537.75 0.01% 131 4010.61 0.01% 192 7627.63 0.01% 176 6539.31 0.01% 23 16959.6 0.01% 71 14150.34 0.01% 14
20-4 12 3198.43 0.01% 12902 3138.18 0.01% 1670 4359.4 0.54% 14400 4256.79 0.24% 14400 7482.73 0.59% 14400 7248.77 0.56% 14400
18 3476.31 0.26% 14400 3377.12 0.17% 14400 4828.18 0.62% 14400 4624.71 0.39% 14400 8876.87 0.45% 14400 8291.19 0.39% 14400
24 3500.35 0.52% 14400 3376.42 0.17% 14400 4894.4 0.81% 14400 4627.29 0.29% 14400 9029.12 0.56% 14400 8286.92 0.42% 14400
30 3644.65 0.22% 14400 3467.32 0.01% 10383 5255.13 0.71% 14400 4850.8 0.45% 14400 9999.52 0.25% 14400 8904.55 0.09% 14400
40 3981.68 0.75% 14400 3703.89 0.01% 8044 5914.89 0.45% 14400 5355.64 0.57% 14400 11754.9 0.07% 14400 10270.96 0.09% 14400
25-2 5 2529.9 0.01% 177 2513.24 0.01% 128 3061.3 0.01% 141 3035.77 0.01% 34 3594.28 0.01% 184 3558.91 0.01% 151
9 2857.2 0.01% 787 2804.59 0.01% 192 3539.88 0.01% 411 3485.14 0.01% 144 4148.77 0.01% 845 4068.28 0.01% 4313
25-3 12 4132.16 0.14% 14400 4064.71 0.05% 14400 5175.19 0.17% 14400 5071.45 0.13% 14400 6221.13 0.15% 14400 6086.02 0.17% 14400
18 4259.23 0.21% 14400 4176.24 0.18% 14400 5357.68 0.23% 14400 5217.53 0.20% 14400 6488.17 0.18% 14400 6278.51 0.14% 14400
24 4527.43 0.22% 14400 4417.51 0.05% 14400 5750.85 0.23% 14400 5562.51 0.15% 14400 6961.73 0.12% 14400 6693.82 0.14% 14400
30-2 5 3058.15 0.07% 14400 3035.55 0.01% 1597 3682.86 0.01% 466 3654.36 0.01% 150 4320.25 0.01% 7775 4277 0.01% 4367
9 3516.77 0.12% 14400 3451.25 0.01% 1943 4358.49 0.07% 14400 4291.04 0.01% 1978 5101.58 0.01% 115 5001.95 0.01% 380
30-3 12 4960.68 0.19% 14400 4877.86 0.12% 14400 6195.55 0.13% 14400 6070.41 0.12% 14400 7447.53 0.14% 14400 7288.22 0.17% 14400
18 5188.53 0.22% 14400 5083.16 0.19% 14400 6514.01 0.2% 14400 6343.84 0.19% 14400 7894.99 0.16% 14400 7630.76 0.12% 14400
24 5525.31 0.29% 14400 5385.95 0.17% 14400 7009.99 0.2% 14400 6781.45 0.15% 14400 8486.73 0.14% 14400 8153.77 0.18% 14400
35-2 5 3424.63 0.29% 14400 3401.37 0.28% 14400 4112.78 0.18% 14400 4071.43 0.13% 14400 4825.02 0.11% 14400 4774.94 0.09% 14400
9 4088.18 0.25% 14400 3401.37 0.25% 14400 5063.59 0.07% 14400 4076.72 0.12% 14400 6090.13 0.21% 14400 4774.94 0.07% 14400
35-3 12 5526.33 0.52% 14400 5431.28 0.36% 14400 6903.64 0.29% 14400 6774.27 0.30% 14400 8371.51 0.33% 14400 8202.31 0.27% 14400
18 5858.63 0.75% 14400 5720.63 0.33% 14400 7308.93 0.32% 14400 7136.05 0.31% 14400 8933.94 0.56% 14400 8658.61 0.29% 14400
24 6372.67 0.68% 14400 6192.61 0.29% 14400 8101.66 0.76% 14400 7829.29 0.30% 14400 9835.33 0.64% 14400 9444.79 0.46% 14400

Note: UB denotes upper bound . When CPLEX finds the optimal solution, the solver outputs Gap=0.01%.

20
Y. Feng et al. Omega 127 (2024) 103082

scenarios and ambiguity sets also applies to other IRP variants. The pro- 2022JM-427], and the Innovation Foundation for Doctor Dissertation
posed method for transforming the worst-case M-CVaR is also suitable of Northwestern Polytechnical University.
for other distributionally robust M-CVaR IRP variant models.
Appendix A. Proof of Theorem 1
8. Conclusions
In this part, we prove Theorem 1.
This study investigates the IRP for which the distributional infor- Proof. For retailer 𝑖, 𝑖 ∈ 𝒩 , denoting 𝛽𝑖 as 𝛽, we can easily obtain that
mation of uncertain demands is partially known. In describing the max𝐏 ∈𝒫 P 𝛽𝐅𝖳𝑖 𝐏𝑖 in the maximization problem (26) is equivalent to
𝑖 𝑖
demands uncertainties of retailers, we do not use demand scenarios { }
𝛽𝐅𝑖 𝐏𝑖 + 𝛽 max 𝐅𝖳𝑖 𝐀𝑖 𝜻 𝑖 ∣ 𝒆𝖳 𝐀𝑖 𝜻 𝑖 = 0, 𝐏0𝑖 + 𝐀𝑖 𝜻 𝑖 ≥ 𝟎, ‖𝜻 𝑖 ‖ ≤ 𝛤𝑖 .
𝖳 0
with precise probabilities of occurrence or specify fixed stochastic 𝜻𝑖
or fuzzy distributions as in the existing literature. Based on limited The optimal value of the above internal maximization problem is
historical demand data, we construct a series of demand scenarios to equivalent to
describe demands uncertainties and develop a framework to construct { }
the ambiguity set. According to the developed framework, two different − min −𝐅𝖳𝑖 𝐀𝑖 𝜻 𝑖 ∣ 𝒆𝖳 𝐀𝑖 𝜻 𝑖 = 0, 𝐏0𝑖 + 𝐀𝑖 𝜻 𝑖 ≥ 𝟎, ‖𝜻 𝑖 ‖ ≤ 𝛤𝑖 .
𝜻𝑖
ambiguity sets are constructed to characterize the ambiguities in the
For the above inner minimization problem, the Lagrange function
probabilities of scenarios occurring. We propose a new distributionally
is shown below
robust M-CVaR IRP model to balance the mean and risk level of
inventory and shortage penalty costs. Based on a given downside risk, 𝑖 (𝜈𝑖 , 𝜰 𝑖 , 𝜋𝑖 , 𝜻 𝑖 ) = −𝐅𝖳𝑖 𝐀𝑖 𝜻 𝑖 + 𝜰 𝖳𝑖 (−𝐏0𝑖 − 𝐀𝑖 𝜻) + 𝜋𝑖 𝒆𝖳 𝐀𝑖 𝜻 𝑖 + 𝜈𝑖 (‖𝜻 𝑖 ‖ − 𝛤𝑖 ),
the proposed model offers an upper bound for the CVaR value of the where 𝜈𝑖 ∈ R+ , 𝜰 𝑖 ∈ R𝑆+ , and 𝜋𝑖 ∈ R are introduced Lagrange
total cost by setting the balance parameters to zero. The constructed multipliers.
model can degenerate to one in different decision-making perspective We set the Lagrangian function 𝑖 (𝜈𝑖 , 𝜰 𝑖 , 𝜋𝑖 , 𝜻 𝑖 ) takes its minimum
by appropriately setting the inputs. value with respect to 𝜻 𝑖 , we obtain the corresponding Lagrange dual
We adopt different duality theories based on different ambiguity function, which is shown as follows:
sets to reformulate the proposed distributionally robust M-CVaR IRP
model into a tractable form. A case study is conducted to demon- 𝑔𝑖 (𝜈𝑖 , 𝜰 𝑖 , 𝜋𝑖 ) = min (𝜈𝑖 , 𝜰 𝑖 , 𝜋𝑖 , 𝜻 𝑖 )
𝜻𝑖
strate the feasibility of the distributionally robust M-CVaR IRP model. { [ ]}
= min (−𝜰 𝖳𝑖 𝐏0𝑖 − 𝛤𝑖 𝜈𝑖 ) − 𝐅𝖳𝑖 𝐀𝑖 𝜻 𝑖 + 𝜰 𝖳𝑖 𝐀𝑖 𝜻 𝑖 − 𝜋𝑖 𝒆𝖳 𝐀𝑖 𝜻 𝑖 − 𝜈𝑖 ‖𝜻 𝑖 ‖
Through experiments and comparative analyses, we observe that am- 𝜻𝑖
[ ]
biguities in demands distributions affect the optimal IRP solution. In = (−𝜰 𝖳𝑖 𝐏0𝑖 − 𝛤𝑖 𝜈𝑖 ) − max (𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 )𝖳 𝜻 𝑖 − 𝜈𝑖 ‖𝜻 𝑖 ‖
𝜻𝑖
addition, the out-of-sample performance analysis shows that the opti-
mal delivery quantity provided by our distributionally robust M-CVaR = (−𝜰 𝖳𝑖 𝐏0𝑖 − 𝛤𝑖 𝜈𝑖 ) − 𝑓 ∗ (𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 ),
IRP model has an obvious advantage over the traditional stochastic where
model in reducing the product shortage rate and can better cope with {
0, ‖𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 ‖∗ ≤ 𝜈𝑖 ,
significant uncertainties in demand. 𝑓𝑖∗ (𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 ) =
We summarize the following proposals for future research on IRP. ∞, other cases,
It is worthwhile for researchers to explore how to design efficient algo- is the conjugate function about = 𝜈𝑖 ‖𝜻‖ [64]. As 𝑔𝑖 (𝜈𝑖 , 𝜰 𝑖 , 𝜋𝑖 )
𝑓 ∗ (𝜻)
rithms (exact decomposition or heuristics) to solve the distributionally generates lower bounds with any 𝜰 𝑖 ∈ R𝑆+ and 𝜈𝑖 ∈ R+ , it follows
robust M-CVaR IRP model and obtain the optimal delivery solution that max𝜈𝑖 ,𝜰 𝑖 ,𝜋𝑖 𝑔(𝜈𝑖 , 𝜰 𝑖 , 𝜋𝑖 ) can be equivalently formulated as the below
for large-scale instances. This study is conducted based on an ambigu- maximization programming model:
ous discrete probability distribution of demand. Another meaningful
⎧
research direction is to explore how to construct an ambiguity set ⎪ 𝜈 max − 𝜰 𝖳𝑖 𝐏0𝑖 − 𝛤𝑖 𝜈𝑖
,𝜰 ,𝜋
based on the continuous probability distribution of demand by utilizing ⎪ 𝑖 𝑖 𝑖
⎨ s.t. ‖𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 ‖∗ ≤ 𝜈𝑖 ,
limited historical data and how to transform the distributionally robust ⎪
M-CVaR IRP model. ⎪ 𝜰 𝑖 ∈ R𝑆+ , 𝜈𝑖 ∈ R+ .
⎩

CRediT authorship contribution statement Additionally, in the transformation of the primal and dual problems
described above, the strong duality holds because the Slater’s condition
Yuqiang Feng: Data curation, Methodology, Software, Writing – holds. Finally, we obtain that max𝐏 ∈𝒫 P 𝛽𝐅𝖳𝑖 𝐏𝑖 in the maximization
𝑖 𝑖
original draft. Ada Che: Conceptualization, Methodology, Validation, problem (26) is equivalent to
Writing – review & editing. Na Tian: Supervision, Validation, Visual- ⎧
⎪𝜈 max − 𝜰 𝖳𝑖 𝐏0𝑖 − 𝛤𝑖 𝜈𝑖
ization. 𝑖 ,𝜰 𝑖 ,𝜋𝑖
⎪
𝛽𝐅𝖳𝑖 𝐏0𝑖 − 𝛽 ⎨ s.t. ‖𝐀𝖳 𝐅 + 𝐀𝖳 𝜰 − 𝐀𝖳 𝒆𝜋 ‖ ≤ 𝜈 ,
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 ∗ 𝑖
Declaration of competing interest ⎪
⎪ 𝜰 𝑖 ∈ R𝑆+ , 𝜈𝑖 ∈ R+
⎩
The authors declare that they have no conflicts of interest. ⎧
⎪𝜈 min 𝛽(𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈𝑖 )
𝑖 ,𝜰 𝑖 ,𝜋𝑖
Data availability ⎪
= ⎨ s.t. ‖𝐀𝖳 𝐅 + 𝐀𝖳 𝜰 − 𝐀𝖳 𝒆𝜋 ‖ ≤ 𝜈 ,
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 ∗ 𝑖
⎪
Data will be made available on request. ⎪ 𝜰 𝑖 ∈ R𝑆+ , 𝜈𝑖 ∈ R+ .
⎩
By applying the above proof techniques again, it follows that
Acknowledgments 1−𝛽 𝖳
max𝐏 ∈𝒫 P 1−𝜅 𝝎𝑖 𝐏𝑖 in the maximization problem (26) is equivalent to
𝑖 𝑖 𝑖
The authors are especially thankful to editors and anonymous re- ⎧ 𝖳
viewers for their valuable comments, which help us to improve the pa- ⎪ max − 𝜰̂ 𝑖 𝐏0𝑖 − 𝛤𝑖 𝜈̂𝑖
𝜈̂𝑖 ,𝜰̂ 𝑖 ,𝜋̂ 𝑖
per a lot. This study is supported by the National Natural Science Foun- 1−𝛽 𝖳 0 1−𝛽 ⎪
1 − 𝜅𝑖 𝑖 𝑖 1 − 𝜅𝑖 ⎨
𝝎 𝐏 − s.t. ‖𝐀𝖳𝑖 𝝎𝑖 + 𝐀𝖳𝑖 𝜰̂ 𝑖 − 𝐀𝖳𝑖 𝒆𝜋̂ 𝑖 ‖∗ ≤ 𝜈̂𝑖 ,
dation of China [grant numbers 72310107003 and 72271201], the Nat- ⎪
ural Science Foundation of Shaanxi Province of China [grant numbers ⎪ 𝜰̂ 𝑖 ≥ 𝟎, 𝜈̂𝑖 ≥ 0
⎩

21
Y. Feng et al. Omega 127 (2024) 103082

⎧ ⎧ 𝖳 0 𝖳 0 ⎫ Appendix C. Proof of Proposition 1

⎪ ⎪ 𝝎 𝐏𝑖 + 𝜰̂ 𝑖 𝐏𝑖 + 𝛤𝑖 𝜈̂𝑖 ⎪
⎪ min (1 − 𝛽) ⎨ ⎬
⎪𝜈̂𝑖 ,𝜰̂ 𝑖 ,𝜋̂𝑖 ⎪ 1 − 𝜅𝑖 ⎪
= ⎨ ⎩ ⎭ In this part, we prove Proposition 1 using proof by contradiction.
⎪
⎪ s.t. ‖𝐀𝖳 𝝎𝑖 + 𝐀𝖳 𝜰̂ 𝑖 − 𝐀𝖳 𝒆𝜋̂ 𝑖 ‖∗ ≤ 𝜈̂𝑖 ,
𝑖 𝑖 𝑖 ∗ ∗
⎪ 𝜰̂ 𝑖 ∈ R𝑆+ , 𝜈̂𝑖 ∈ R+ , Proof. Under the ambiguity set (14), we assume that (𝜼∗ , 𝝓∗ , 𝐅 , 𝝎 , 𝝂 ∗ ,
⎩ ∗ ∗
𝜰 , 𝝅 ∗ , 𝝂̂ ∗ , 𝜰̂ , 𝝅̂ ∗ ) is the optimal solution to the model (27). We recall
∗ ∗
where 𝜈̂𝑖 ∈ R+ , 𝜰̂ 𝑖 ∈ R𝑆+ , and 𝜋̂ 𝑖 ∈ R+ are Lagrange multipliers. that 𝐅 = (𝐅𝑖 )𝑖∈𝒩 and 𝝎 = (𝝎𝑖 )𝑖∈𝒩 . The following relationships can be
In summary, based on the ambiguity set (14), the maximization obtained:
problem (26) can be re-expressed as { }
∑
∗𝖳
max 𝐅𝑖 𝐏𝑖
⎧ { } 1−𝛽 { 𝖳 0 𝖳
}
𝐏𝑖 ∈𝒫𝑖P
⎪min 𝛽 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈𝑖 + 𝝎𝑖 𝐏𝑖 + 𝜰̂ 𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈̂𝑖 𝑖∈𝒩
{ }
⎪ 1 − 𝜅𝑖 ∑ { ∗𝖳 }
⎪ s.t. = 𝐅∗𝖳 0 𝖳 0
𝑖 𝐏𝑖 + max 𝐅𝑖 𝐀𝑖 𝜻 𝑖 ∣ 𝒆 𝐀𝑖 𝜻 𝑖 = 0, 𝐏𝑖 + 𝐀𝑖 𝜻 𝑖 ≥ 𝟎, ‖𝜻 𝑖 ‖ ≤ 𝛤𝑖
‖𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 𝑖 − 𝐀𝖳𝑖 𝒆𝜋𝑖 ‖∗ ≤ 𝜈𝑖 ,
⎨ □ 𝑖∈𝒩
𝜻𝑖
⎪ ‖𝐀𝖳𝑖 𝝎𝑖 + 𝐀𝖳𝑖 𝜰̂ 𝑖 − 𝐀𝖳𝑖 𝒆𝜋̂ 𝑖 ‖∗ ≤ 𝜈̂𝑖 , ∑{ }
⎪ ≤ 𝐅∗𝖳 0
𝑖 𝐏𝑖 + 𝜰 ∗𝖳 0
𝑖 𝐏𝑖 + 𝜈𝑖∗ ,
⎪ 𝜰 𝑖 ∈ R𝑆+ , 𝜈𝑖 ∈ R+ , 𝜰̂ 𝑖 ∈ R𝑆+ , 𝜈̂𝑖 ∈ R+ . 𝑖∈𝒩
⎩
where the inequality holds due to Lagrangian weak duality.
In the same way, the following relationships also can be obtained:
Appendix B. Proof of Theorem 2 { }
∑ 𝝎∗𝖳
𝑖 𝐏𝑖
max
𝑖∈𝒩 𝐏𝑖 ∈𝒫𝑖P 1 − 𝜅𝑖
In this part, we prove Theorem 2. { { }}
∑ 𝝎∗𝖳 0
𝑖 𝐏𝑖 𝑖 𝐀𝑖 𝜻 𝑖 ∣ 𝒆 𝐀𝑖 𝜻 𝑖 = 0, 𝐏𝑖 + 𝐀𝑖 𝜻 𝑖 ≥ 𝟎, ‖𝜻 𝑖 ‖ ≤ 𝛤𝑖
max𝜻 𝑖 𝝎∗𝖳 𝖳 0
= +
Proof. For retailer 𝑖, 𝑖 ∈ 𝒩 , we can easily obtain that max𝐏 ∈𝒫 B 𝛽𝐅𝖳𝑖 𝐏𝑖 𝑖∈𝒩
1 − 𝜅𝑖 1 − 𝜅𝑖
𝑖 𝑖
in the maximization problem (26) is equivalent to { ∗𝖳 }
∑ 𝝎∗𝖳 0
𝑖 𝐏𝑖 + 𝜰̂ 𝑖 𝐏0𝑖 + 𝜈̂𝑖∗
{ } ≤ .
𝛽𝐅𝖳𝑖 𝐏0𝑖 + 𝛽 max 𝐅𝖳𝑖 𝜻 𝑖 ∣ 𝐞𝖳 𝜻 𝑖 = 0, 𝜻 𝐿 𝑈
𝑖 ≤ 𝜻𝑖 ≤ 𝜻𝑖 . 𝑖∈𝒩
1 − 𝜅𝑖
𝜻𝑖
At the same time, 𝜼∗ and 𝝓∗ are feasible to constraints (1)–(13),
For the above maximization problem, in light of the duality theory ∗ ∗
(25b)–(25c). Therefore, the solution (𝜼∗ , 𝝓∗ , 𝐅 , 𝝎 ) is feasible for model
of linear optimization, its dual problem is expressed as follows:
(25).
𝖳 𝐿 𝖳
⎧𝜏 min (𝜻 𝑈
𝑖 ) 𝝆𝑖 − (𝜻 𝑖 ) 𝜸 𝑖 Letting the solution (𝜼2∗ , 𝝓2∗ , 𝐅2∗ , 𝝎2∗ ) is optimal for model (25),
,𝜸 ,𝝆
⎪𝑖 𝑖 𝑖 we obtain
⎪ s.t. 𝒆𝜏𝑖 + 𝝆𝑖 − 𝜸 𝑖 = 𝐅𝑖 , [
⎨ { }]
∑ ∑ 1−𝛽
⎪ 𝝆𝑖 ∈ R𝑆+ , TC1 (𝜼2∗ ) + (1 − 𝛽) 𝜙2∗𝑖 + 𝛽𝐅2∗𝖳𝑖 max
𝐏𝑖 + 𝝎2∗𝖳
𝑖 𝐏𝑖
⎪ 𝑖∈𝒩 𝑖∈[𝑁] 𝐏𝑖 ∈𝒫𝑖
P 1 − 𝜅𝑖
⎩ 𝜸 𝑖 ∈ R𝑆+ , [ { }]
∑ ∑ 1 − 𝛽 ∗𝖳
where 𝜏𝑖 ∈ R, 𝝆𝑖 ∈ R𝑆+ , and 𝜸 𝑖 ∈ R𝑆+ are the introduced dual variables ≤ TC1 (𝜼∗ ) + (1 − 𝛽) 𝜙∗𝑖 + max 𝛽𝐅∗𝖳 𝐏𝑖 + 𝝎 𝐏 .
𝑖∈𝒩 𝑖∈[𝑁] 𝐏𝑖 ∈𝒫𝑖
P 𝑖 1 − 𝜅𝑖 𝑖 𝑖
which respond to the above constraints. The strong duality holds as the
{ }
set 𝐞𝖳 𝜻 𝑖 = 0, 𝜻 𝐿 𝑈
𝑖 ≤ 𝜻𝑖 ≤ 𝜻𝑖 is nonempty and bounded. Letting (𝐅2∗ , 𝝂2∗ , 𝜰 2∗ , 𝝅2∗ ) be the optimal solution to model (28).
It follows that max𝐏 ∈𝒫 B 𝛽𝐅𝖳𝑖 𝐏𝑖 in the maximization problem (26) By strong duality, we obtain
𝑖 𝑖
is equivalent to { }
∑ ∑{ }
𝛽(𝐅𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 𝐿 𝖳 max 𝐅∗𝖳 𝐏 = 𝐅2∗𝖳 0 ∗𝖳 0 ∗
𝑖 𝐏𝑖 + 𝜰 2𝑖 𝐏𝑖 + 𝜈2𝑖 .
⎧𝜏 min
,𝜸 ,𝝆 𝑖 ) 𝝆𝑖 − (𝜻 𝑖 ) 𝜸 𝑖 )
𝑖
𝐏𝑖 ∈𝒫𝑖P
𝑖
⎪𝑖 𝑖 𝑖 𝑖∈𝒩 𝑖∈𝒩
⎪ s.t. 𝒆𝜏𝑖 + 𝝆𝑖 − 𝜸 𝑖 = 𝐅𝑖 ,
⎨ ̂ ∗ , 𝜰̂ 2∗ , 𝝅2
Similarly, we let (𝝎2∗ , 𝝂2 ̂ ∗ ) be the optimal solution to
⎪ 𝝆𝑖 ∈ R𝑆+ ,
⎪ model (29). By strong duality, we obtain
⎩ 𝜸 𝑖 ∈ R𝑆+ . { } { }
∑ 𝝎2∗𝖳
𝑖 𝐏𝑖
∑ 𝝎2∗𝖳 0 ̂ ∗𝖳 0 ̂ ∗
𝑖 𝐏𝑖 + 𝜰 2𝑖 𝐏𝑖 + 𝜈2 𝑖
By applying the above proof techniques again, we obtain that max = .
1−𝛽 𝖳 𝑖∈𝒩 𝐏𝑖 ∈𝒫𝑖P 1 − 𝜅𝑖 𝑖∈𝒩
1 − 𝜅𝑖
max𝐏 ∈𝒫 B 1−𝜅 𝝎𝑖 𝐏𝑖 in the maximization problem (26) is equivalent to
𝑖 𝑖 𝑖
{ } From the other constraints in model (25), (28), and (29), we
⎧ 𝝎𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 ̂ − (𝜻 𝐿 )𝖳 𝜸̂
𝑖 ) 𝝆 𝑖 𝑖 𝑖 find that (𝜼2∗ , 𝝓2∗ , 𝐅2∗ , 𝝎2∗ , 𝝂2∗ , 𝜰 2∗ , 𝝅2∗ , 𝝂2 ̂ ∗ , 𝜰̂ 2∗ , 𝝅2
̂ ∗ ) is optimal for
⎪ min (1 − 𝛽)
⎪𝜏̂𝑖 ,𝜸̂ 𝑖 ,𝝆̂ 𝑖 1 − 𝜅𝑖 model (27), which conflicts with the above assumption that the so-
∗
⎪ ∗ ∗ ∗
lution (𝜼∗ , 𝝓∗ , 𝐅 , 𝝎 , 𝝂 ∗ , 𝜰 , 𝝅 ∗ , 𝝂̂ ∗ , 𝜰̂ , 𝝅̂ ∗ ) is optimal for model (27).
⎨ s.t. 𝒆𝜏̂𝑖 + 𝝆̂ 𝑖 − 𝜸̂ 𝑖 = 𝝎𝑖 , ∗ ∗
⎪ Therefore, the solution (𝜼∗ , 𝝓∗ , 𝐅 , 𝝎 ) is optimal for model (25).
⎪ 𝝆̂ 𝑖 ∈ R𝑆+ ,
On the contrary, we assume that the solution (𝜼1∗ , 𝝓1∗ , 𝐅1∗ , 𝝎1∗ )
⎪ 𝜸̂ 𝑖 ∈ R𝑆+ ,
⎩ is optimal for model (25), the solution (𝐅1∗ , 𝝂1∗ , 𝜰 1∗ , 𝝅1∗ ) is optimal
for model (28), and the solution (𝝎1∗ , 𝝂1 ̂ ∗ , 𝜰̂ 1∗ , 𝝅1
̂ ∗ ) is optimal for
where 𝜏̂𝑖 ∈ R, 𝝆̂ 𝑖 ∈ R𝑆+ , and 𝜸̂ 𝑖 ∈ R𝑆+ are the introduced dual variables.
In summary, under ambiguity set (15), the maximization problem model (29). If (𝜼1∗ , 𝝓1∗ , 𝐅1∗ , 𝝎1∗ , 𝝂1∗ , 𝜰 1∗ , 𝝅1∗ , 𝝂1 ̂ ∗ , 𝜰̂ 1∗ , 𝝅1
̂ ∗ ) is not the
(26) can be re-expressed as optimal solution about the model (27), then there exists a solution
̂ ∗ , 𝜰̂ 3∗ , 𝝅3
(𝜼3∗ , 𝝓3∗ , 𝐅3∗ , 𝝎3∗ , 𝝂3∗ , 𝜰 3∗ , 𝝅3∗ , 𝝂3 ̂ ∗ ) such that
⎧ { } 1−𝛽 { 𝖳 0 }
𝛽 𝐅𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈𝑖 )𝖳 𝝅 𝑖 − (𝜻 𝐿𝑖 )𝖳 𝜸 𝑖 + 𝝎𝑖 𝐏𝑖 + (𝜻 𝑈𝑖 )𝖳 𝝅̂ 𝑖 − (𝜻 𝐿𝑖 )𝖳 𝜸̂ 𝑖 [ { }]
⎪ min 1 − 𝜅𝑖 ∑ ∑ 1−𝛽
⎪ TC1 (𝜼3∗ ) + (1 − 𝛽) 𝜙3∗𝑖 + maxP 𝛽𝐅3∗𝖳𝑖 𝐏𝑖 + 1 − 𝜅 𝝎3𝑖 𝐏𝑖
∗𝖳
⎪s.t. 𝒆𝜏𝑖 + 𝝆𝑖 − 𝜸 𝑖 = 𝐅𝑖 , 𝑖∈𝒩 𝑖∈𝒩
𝐏𝑖 ∈𝒫𝑖 𝑖
⎨ □ [ { }]
⎪ ∑ ∑ 1−𝛽
∗ ∗ ∗𝖳
⎪
𝒆𝜏̂𝑖 + 𝝆̂ 𝑖 − 𝜸̂ 𝑖 = 𝝎𝑖 , ≤ TC1 (𝜼1 ) + (1 − 𝛽) 𝜙1𝑖 + maxP 𝛽𝐅1𝑖 𝐏𝑖 + 𝝎1∗𝖳
𝑖 𝐏𝑖
𝐏𝑖 ∈𝒫𝑖 1 − 𝜅𝑖
⎪ 𝝆𝑖 ∈ R𝑆+ , 𝜸 𝑖 ∈ R𝑆+ , 𝝆̂ 𝑖 ∈ R𝑆+ , 𝜸̂ 𝑖 ∈ R𝑆+ . 𝑖∈𝒩 𝑖∈𝒩
⎩

22
Y. Feng et al. Omega 127 (2024) 103082

Table A.1
Delivery quantity offered by the Ball model.
Supermarket Total
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
𝑡=1 920 1321 1028 1230 1145 870 773 1099 844 996 1114 712 1110 1207 1307 1074 969 1010 698 640 20067
𝑡=2 500 984 876 750 1147 890 852 1152 796 516 467 682 1202 876 767 784 703 403 502 0 14849
𝑡=3 500 944 876 750 901 750 508 919 550 516 600 467 876 876 767 1083 566 697 900 0 14046

Fig. A1. Delivery route provided by the Ball model under 𝛽 = 0.5 and 𝜅 = 0.9.

According to the first half of this proof, we can see that the solution
(𝜼3∗ , 𝝓3∗ , 𝐅3∗ , 𝝎3∗ ) is optimal for model (25), which conflicts with the
assumption that the solution (𝜼1∗ , 𝝓1∗ , 𝐅1∗ , 𝝎1∗ ) is optimal for model
̂ ∗ , 𝜰̂ 1∗ , 𝝅1
(25). Hence, the solution (𝜼1∗ , 𝝓1∗ , 𝐅1∗ , 𝝎1∗ , 𝝂1∗ , 𝜰 1∗ , 𝝅1∗ , 𝝂1 ̂ ∗)
is optimal for model (27). □

Appendix D. Verifying the feasibility of using ‖𝜻 𝒊 ‖𝟐 in (14)

First, we verify that using ‖𝜻 𝑖 ‖2 is also feasible in the ambiguity set

(14). Because the geometry of ‖𝜻 𝑖 ‖2 is a ball form, we denote the model
consisting of model (27) and constraints (31)–(34) as Ball model. The
settings of the parameters of Ball model are consistent with those of
the parameters of Poly model in Section 6.3. The experimental results
of the Ball model are as follows.
The objective value is 5457.62, where TC1 = 1034.79, SOM =
3397.78, SOV = 5440.96, and SOCV = 5447.87. The objective value of
the Ball model is higher than that of the Poly model. This is because
the scale of the ‖𝜻 𝑖 ‖2 -based ambiguity set is larger than the scale of
‖𝜻 𝑖 ‖1 -based ambiguity set for the same 𝛤𝑖 . Meanwhile, the difference Fig. A2. Frequency distribution histogram of performance of the delivery quantities
in the objective values of these two models is mainly caused by the given by the Ball and Nominal models.
large difference in SOM value of these two models. Regarding TC1 ,
SOV and SOCV, these two models provide closer values. Table A.1 and
Fig. A1 present the optimal delivery quantity and route offered by the found within the specified solving time limit under 𝑆 = 27. Although
Ball model. We find that there is a significant difference between the the optimal solution of the Ball model is also found in some cases,
Ball model and Poly model in terms of optimal delivery quantity, but the solution time used is much more than that with the Poly model.
in terms of optimal delivery route, the difference is smaller. However, we declare again that both models are theoretically feasible.
Moreover, we use Algorithm 2 to perform the performance tests. The
test results are listed in Table A.2. We draw the frequency distribution
histogram (see Fig. A2) to further demonstrate the advantage of the Appendix E. Tractable reformulation of model (35)
optimal delivery quantity provided by the Ball model in out-of-sample
performance. By comparing Table A.2 and Table 6, we see that the
Aver and Wor performances of the Ball model prevail over those of Based on the ambiguity set (14), we derive the tractable reformula-
the nominal model but slightly worse than those of the Poly model in tion of the distributionally robust M-SD IRP model (35) in Section 6.6.
most cases. The Ball model has the best SR performance for all 𝛥. From We need to deduce the tractable forms of max𝐏𝑖 ∈𝒫𝑖 𝐅𝖳𝑖 𝐏𝑖 , max𝐏𝑖 ∈𝒫𝑖 𝓵 𝖳𝑖 𝐏𝑖 ,
Fig. A2, we observe that the frequency distribution histogram and fitted and − min𝐏𝑖 ∈𝒫𝑖 𝐅𝖳𝑖 𝐏𝑖 . According to the proof process of Theorem 1, we
profile of the performance of the optimal delivery quantity given by can express the aforementioned three parts respectively as:
the Ball model are to the left of that given by the nominal model. This
⎧
shows that the Ball model can also offer a promising delivery solution min
⎪𝜈1 ,𝜰 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 1𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈1𝑖
𝑖 1𝑖 ,𝜋1𝑖
with excellent out-of-sample performance. ⎪
Second, we demonstrate the difference between the Poly model and ⎨ s.t. ‖𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 1𝑖 − 𝐀𝖳𝑖 𝒆𝜋1𝑖 ‖∗ ≤ 𝜈1𝑖 ,
⎪
Ball model in terms of computation. We set 𝑆 to take values from 3 to 27 ⎪ 𝜰 1𝑖 ∈ R𝑆+ , 𝜋1𝑖 ∈ R, 𝜈1𝑖 ∈ R+ ,
with an interval of 3. The solving time is limited to 14400 seconds under ⎩
each 𝑆. The computational results of these two models is organized in ⎧
Table A.3. min
⎪𝜈2 ,𝜰 𝓵 𝖳𝑖 𝐏0𝑖 + 𝜰 2𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈2𝑖
2 ,𝜋2
With Table A.3, we observe that the optimal solution of the Poly ⎪ 𝑖 𝑖 𝑖
⎨ s.t. ‖𝐀𝖳𝑖 𝓵 𝑖 + 𝐀𝖳𝑖 𝜰 2𝑖 − 𝐀𝖳𝑖 𝒆𝜋2𝑖 ‖∗ ≤ 𝜈2𝑖 ,
model can be found in a shorter time, which is significantly lower than ⎪
the specified time limit. The optimal solution of the Ball model is not ⎪ 𝜰 2𝑖 ∈ R𝑆+ , 𝜋2𝑖 ∈ R, 𝜈2𝑖 ∈ R+ ,
⎩

23
Y. Feng et al. Omega 127 (2024) 103082

Table A.2
Out-of-sample performance of the delivery quantity provided by the Ball model.
Indicator Value Indicator Value Indicator Value
𝛥=0 Aver 2713.35 𝛥 = 0.2 Aver 3368.91 𝛥 = 0.4 Aver 4162.3
Worst 3593.16 Worst 4641.02 Worst 6144.5
SR 44.5% SR 55.36% SR 62.2%
𝛥 = 0.1 Aver 3023.2 𝛥 = 0.3 Aver 3762.8 𝛥 = 0.5 Aver 4604.16
Worst 4076.82 Worst 5539.7 Worst 6941.78
SR 51.2% SR 59.44% SR 64.74%

Table A.3
The computational results of the Poly and Ball models under different 𝑆.
Poly model Ball model Poly model Ball model
UB Gap Time(s) UB Gap Time(s) UB Gap Time(s) UB Gap Time(s)
𝑆 =6 2150.59 0.01% 354 2186.66 0.01% 1920 𝑆 = 18 5239.63 0.01% 1001 5457.62 0.01% 5163
𝑆 =9 2404.38 0.01% 100 2482.73 0.01% 220 𝑆 = 21 5507.84 0.01% 550 5767.88 0.01% 2040
𝑆 = 12 2872.73 0.01% 308 2971.57 0.01% 1589 𝑆 = 24 5573.27 0.01% 128 5871.85 0.01% 1705
𝑆 = 15 2893.99 0.01% 621 3009.32 0.01% 1523 𝑆 = 27 5733.64 0.01% 454 6048.39 0.10% 14400

Note: UB denotes upper bound. When CPLEX finds the optimal solution, the solver outputs Gap=0.01%.

The tractable forms of the distributionally robust M-SD IRP model

⎧
min
⎪𝜈3 ,𝜰 − 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 3𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈3𝑖 (35) with the ambiguity set (15) is shown as:
3 ,𝜋3
⎪ 𝑖 𝑖 𝑖 ∑{ ( )
⎨ s.t. ‖ − 𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 3𝑖 − 𝐀𝖳𝑖 𝒆𝜋3𝑖 ‖∗ ≤ 𝜈3𝑖 , ⎧min TC1 (𝒙, 𝒒, 𝒛, 𝒚) + 𝛽𝑖 𝐅𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 𝐿 𝖳
𝑖 ) 𝝆1𝑖 − (𝜻 𝑖 ) 𝜸1𝑖
⎪ ⎪ 𝑖∈𝒩
⎪ 𝜰 3𝑖 ∈ R𝑆+ , 𝜋3𝑖 ∈ R, 𝜈3𝑖 ∈ R+ , ⎪ ( )}
⎩ ⎪ + (1 − 𝛽𝑖 ) 𝓵 𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 𝐿 𝖳
𝑖 ) 𝝆2𝑖 − (𝜻 𝑖 ) 𝜸2𝑖
⎪ 𝑠 𝑠 𝑠 𝑠 𝖳 0 𝑈 𝖳 𝐿 𝖳
where 𝜈1𝑖 , 𝜰 1𝑖 , 𝜋1𝑖 , 𝜈2𝑖 , 𝜰 2𝑖 , 𝜋2𝑖 , 𝜈3𝑖 , 𝜰 3𝑖 , and 𝜋3𝑖 are introduced dual ⎪s.t. 𝓁𝑖 ≥ F𝑖 (𝜶 𝑖 , 𝜹𝑖 ) − 𝐅𝑖 𝐏𝑖 + (𝜻 𝑖 ) 𝝆3𝑖 − (𝜻 𝑖 ) 𝜸3𝑖 ,
variables. ⎪ 𝓁𝑖𝑠 ≥ 0, ∀𝑠 ∈ 𝒮 , 𝑖 ∈ 𝒩 ,
⎪
The tractable forms of the distributionally robust M-SD IRP model ⎨ 𝒆𝜏1𝑖 + 𝝆1𝑖 − 𝜸1𝑖 = 𝐅𝑖 , 𝒆𝜏2𝑖 + 𝝆2𝑖 − 𝜸2𝑖 = 𝓵 𝑖 ,
(35) with the ambiguity set (14) is shown as: ⎪
⎪ 𝒆𝜏3𝑖 + 𝝆3𝑖 − 𝜸3𝑖 = −𝐅𝑖 , ∀𝑖 ∈ 𝒩 ,
⎧min ∑{ ( ) ⎪
TC1 (𝒙, 𝒒, 𝒛, 𝒚) + 𝛽𝑖 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 1𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈1𝑖 ⎪ 𝝆1𝑖 , 𝜸1𝑖 , 𝝆2𝑖 , 𝜸2𝑖 , 𝝆3𝑖 , 𝜸3𝑖 ∈ R𝑆+ , ∀𝑖 ∈ 𝒩 ,
⎪
⎪ 𝑖∈𝒩 ⎪
( )} ⎪ 𝜏1𝑖 , 𝜏2𝑖 , 𝜏3𝑖 ∈ R, ∀𝑖 ∈ 𝒩 ,
⎪ + (1 − 𝛽𝑖 ) 𝓵 𝖳𝑖 𝐏0𝑖 + 𝜰 2𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈2𝑖
⎪ ⎪
⎪s.t. 𝓁𝑖𝑠 ≥ F𝑠𝑖 (𝜶 𝑠𝑖 , 𝜹𝑠𝑖 ) − 𝐅𝖳𝑖 𝐏0𝑖 + 𝜰 3𝖳𝑖 𝐏0𝑖 + 𝛤𝑖 𝜈3𝑖 , 𝓁𝑖𝑠 ≥ 0, ∀𝑠 ∈ 𝒮 , 𝑖 ∈ 𝒩 , ⎩ Constraints ∶ (1)–(13).
⎪
⎪ ‖𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 1𝑖 − 𝐀𝖳𝑖 𝒆𝜋1𝑖 ‖∗ ≤ 𝜈1𝑖 ,
⎪
⎨ ‖𝐀𝖳𝑖 𝓵 𝑖 + 𝐀𝖳𝑖 𝜰 2𝑖 − 𝐀𝖳𝑖 𝒆𝜋2𝑖 ‖∗ ≤ 𝜈2𝑖 , ∀𝑖 ∈ 𝒩 , Appendix F. Note on parameters setting for large-scale cases
⎪
⎪ ‖ − 𝐀𝖳𝑖 𝐅𝑖 + 𝐀𝖳𝑖 𝜰 3𝑖 − 𝐀𝖳𝑖 𝒆𝜋3𝑖 ‖∗ ≤ 𝜈3𝑖 , ∀𝑖 ∈ 𝒩 ,
⎪
⎪ 𝜰 1𝑖 , 𝜰 2𝑖 , 𝜰 3𝑖 ∈ R𝑆+ , ∀𝑖 ∈ 𝒩 , We give explanations about the setting of parameter values in the in-
⎪ stances with 25, 30, and 35 retailers, respectively. Coelho & Laporte [63]
⎪ 𝜋1𝑖 , 𝜈1𝑖 , 𝜋2𝑖 , 𝜈2𝑖 , 𝜋3𝑖 , 𝜈3𝑖 ∈ R+ , ∀𝑖 ∈ 𝒩 ,
⎪ does not allow product shortages to occur. The transportation cost
⎪ Constraints ∶ (1)–(13). between nodes generated by Coelho & Laporte [63] is much larger than
⎩
the inventory costs at nodes. When we use these original data directly,
Next, we deduce the tractable forms of max𝐏𝑖 ∈𝒫𝑖 𝐅𝖳𝑖 𝐏𝑖 , max𝐏𝑖 ∈𝒫𝑖 𝓵 𝖳𝑖 𝐏𝑖 , the vendor does not deliver products to the majority of retailers because
and − min𝐏𝑖 ∈𝒫𝑖 𝐅𝖳𝑖 𝐏𝑖 under the ambiguity set (15). In the light of the this study allows for product shortages. In this context, we modify the
proof process of Theorem 2, we can reformulate the aforementioned transportation cost to one-fiftieth of the original value, the inventory
three parts respectively as: cost to ten times the original value, and also the initial inventory to one-
fourth of the original value. Other parameter values (e.g., the maximum
⎧
⎪𝜏1 min 𝐅𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 𝐿 𝖳
𝑖 ) 𝝆1𝑖 − (𝜻 𝑖 ) 𝜸1𝑖
load of vehicles, the maximum inventory capacity of retailers, the
𝑖 ,𝜸1𝑖 ,𝝆1𝑖
⎪ historical demand data of retailers) are set in accordance with Coelho
⎨ s.t. 𝒆𝜏1𝑖 + 𝝆1𝑖 − 𝜸1𝑖 = 𝐅𝑖 , & Laporte [63].
⎪
⎪ 𝝆1𝑖 ∈ R𝑆+ , 𝜸1𝑖 ∈ R𝑆+ , 𝜏1𝑖 ∈ R,
⎩
Appendix G. Supplementary material
⎧
⎪𝜏2 min 𝓵 𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 𝐿 𝖳
𝑖 ) 𝝆2𝑖 − (𝜻 𝑖 ) 𝜸2𝑖
𝑖 ,𝜸2𝑖 ,𝝆2𝑖
⎪ Supplementary material mainly lists the following information: lim-
⎨ s.t. 𝒆𝜏2𝑖 + 𝝆2𝑖 − 𝜸2𝑖 = 𝓵 𝑖 , ited historical demand data, the maximum inventory level, the trans-
⎪
⎪ 𝝆2𝑖 ∈ R𝑆+ , 𝜸2𝑖 ∈ R𝑆+ , 𝜏2𝑖 ∈ R, portation costs between nodes, the demand scenarios, and the final
⎩ nominal value of probability vector after normalizing.
⎧
⎪𝜏3 min − 𝐅𝖳𝑖 𝐏0𝑖 + (𝜻 𝑈 𝖳 𝐿 𝖳
𝑖 ) 𝝆3𝑖 − (𝜻 𝑖 ) 𝜸3𝑖
𝑖 ,𝜸3𝑖 ,𝝆3𝑖 References
⎪
⎨ s.t. 𝒆𝜏3𝑖 + 𝝆3𝑖 − 𝜸3𝑖 = −𝐅𝑖 ,
⎪ [1] Ha AY, Luo H, Shang W. Supplier encroachment, information sharing,
⎪ 𝝆3𝑖 ∈ R𝑆+ , 𝜸3𝑖 ∈ R𝑆+ , 𝜏3𝑖 ∈ R,
⎩ and channel structure in online retail platforms. Prod Oper Manage
2022;31(3):1235–51.
where 𝝆1𝑖 , 𝜸1𝑖 , 𝜏1𝑖 , 𝝆2𝑖 , 𝜸2𝑖 , 𝜏2𝑖 , 𝝆3𝑖 , 𝜸3𝑖 , and 𝜏3𝑖 are introduced dual [2] Raweewan M, Ferrell Jr WG. Information sharing in supply chain collaboration.
variables. Comput Ind Eng 2018;126:269–81.

24
Y. Feng et al. Omega 127 (2024) 103082

[3] Dong Y, Dresner M, Yao Y. Beyond information sharing: An empirical analysis [34] Nikzad E, Bashiri M, Oliveira F. Two-stage stochastic programming approach for
of vendor-managed-inventory. Prod Oper Manage 2014;23(5):817–28. the medical drug inventory routing problem under uncertainty. Comput Ind Eng
[4] Ru J, Shi R, Zhang J. When does a supply chain member benefit from 2019;128:358–70.
vendor-managed-inventory? Prod Oper Manage 2018;27(5):807–21. [35] Onggo BS, Panadero J, Corlu CG, Juan AA. Agri-food supply chains with
[5] Mahmutoğulları Ö, Yaman H. A branch-and-cut algorithm for the inventory stochastic demands: A multi-period inventory routing problem with perishable
routing problem with product substitution. Omega 2023;115:102752. products. Simul Model Pract Theory 2019;97:101970.
[6] Shao S, Lai KK, Ge B. A multi-period inventory routing problem with procure- [36] Soysal M. Closed-loop inventory routing problem for returnable transport items.
ment decisions: a case in China. Ann Oper Res 2022;1–29. http://dx.doi.org/10. Transp Res D 2016;48:31–45.
1007/s10479-021-04345-0. [37] Alvarez A, Cordeau JF, Jans R, Munari P, Morabito R. Inventory routing under
[7] Bell WJ, Dalberto LM, Fisher ML, Greenfield AJ, Jaikumar R, Kedia P, et al. stochastic supply and demand. Omega 2021;102:102304.
Improving the distribution of industrial gases with an on-line computerized [38] Mousavi R, Bashiri M, Nikzad E. Stochastic production routing problem for
routing and scheduling optimizer. Interfaces 1983;13(6):4–23. perishable products: Modeling and a solution algorithm. Comput Oper Res
[8] Archetti C, Speranza MG. The inventory routing problem: the value of 2022;142:105725.
integration. Int Trans Oper Res 2016;23(3):393–407. [39] Niakan F, Rahimi M. A multi-objective healthcare inventory routing problem: a
[9] Yu Y, Chu F, Chen H. A stackelberg game and its improvement in a VMI system fuzzy possibilistic approach. Transp Res E 2015;80:74–94.
with a manufacturing vendor. European J Oper Res 2009;192(3):929–48. [40] Rahimi M, Baboli A, Rekik Y. Multi-objective inventory routing problem: A
[10] Liu M, Liu X, Chu F, Zheng F, Chu C. Distributionally robust inventory routing stochastic model to consider profit, service level and green criteria. Transp Res
problem to maximize the service level under limited budget. Transp Res E E 2017;101:59–83.
2019;126:190–211. [41] Kazemi SM, Rabbani M, Tavakkoli-Moghaddam R, Shahreza FA. Blood inventory-
[11] Shang X, Zhang G, Jia B, Almanaseer M. The healthcare supply location- routing problem under uncertainty. J Intell Fuzzy Systems 2017;32(1):467–81.
inventory-routing problem: A robust approach. Transp Res E 2022;158:102588. [42] Hu H, Li J, Li X. A credibilistic goal programming model for inventory routing
[12] Shaabani H, Hoff A, Hvattum LM, Laporte G. A matheuristic for the multi-product problem with hazardous materials. Soft Comput 2018;22(17):5803–16.
maritime inventory routing problem. Comput Oper Res 2023;154:106214. [43] Ghasemkhani A, Tavakkoli-Moghaddam R, Shahnejat-Bushehri S, Momen S,
[13] Zhang Y, Chu F, Che A, Li Y. Closed-loop inventory routing problem for Tavakkoli-Moghaddam H. An integrated production inventory routing prob-
perishable food with returnable transport items selection. Int J Prod Res lem for multi perishable products with fuzzy demands and time windows.
2023;1–21. IFAC-PapersOnLine 2019;52(13):523–8.
[14] Zhang P, Liu Y, Yang G, Zhang G. A multi-objective distributionally robust model [44] Avila-Torres PA, Arratia-Martinez NM. Inventory routing problem with fuzzy
for sustainable last mile relief network design problem. Ann Oper Res 2020;1–42. demand and deliveries with priority. In: Artificial intelligence in industry 4.0
[15] Feng Y, Chen Y, Liu Y. Optimising two-stage robust supplier selection and and 5G technology. 2022, p. 267–82.
order allocation problem under risk-averse criterion. Int J Prod Res 2022;1–25. [45] Adamo JM. Fuzzy decision trees. Fuzzy Sets and Systems 1980;4(3):207–19.
http://dx.doi.org/10.1080/00207543.2022.2127963. [46] Ghasemkhani A, Tavakkoli-Moghaddam R, Rahimi Y, Shahnejat-Bushehri S,
[16] Wang W, Yang K, Yang L, Gao Z. Distributionally robust chance-constrained Tavakkoli-Moghaddam H. Integrated production–inventory-routing problem for
programming for multi-period emergency resource allocation and vehicle routing multi-perishable products under uncertainty by meta-heuristic algorithms. Int J
in disaster response operations. Omega 2023;120:102915. Prod Res 2022;60(9):2766–86.
[17] Li B, Arreola-Risa A. Minimizing conditional value-at-risk under a modified [47] Solyalı O, Cordeau JF, Laporte G. Robust inventory routing under demand
basestock policy. Prod Oper Manage 2022;31(4):1822–38. uncertainty. Transp Sci 2012;46(3):327–40.
[18] Wang L, Yao DD. Production planning with risk hedging under a conditional [48] Jafarkhan F, Yaghoubi S. An efficient solution method for the flexible and robust
value at risk objective. Oper Res 2023. http://dx.doi.org/10.1287/opre.2022. inventory-routing of red blood cells. Comput Ind Eng 2018;117:191–206.
2423. [49] Liu P, Hendalianpour A, Razmi J, Sangari MS. A solution algorithm for inte-
[19] Sehgal R, Sharma A, Mansini R. Worst-case analysis of Omega-VaR ratio grated production–inventory-routing of perishable goods with transshipment and
optimization model. Omega 2023;114:102730. uncertain demand. Complex Intell Syst 2021;7(3):1349–65.
[20] Rockafellar RT, Uryasev S. Optimization of conditional value-at-risk. J Risk [50] Mulvey JM, Ruszczyński A. A new scenario decomposition method for large-scale
2000;2:21–42. stochastic optimization. Oper Res 1995;43(3):477–90.
[21] Chen YM, Lin CT. A coordinated approach to hedge the risks in stochastic [51] Gholami-Zanjani SM, Gholipour S, Rabbani M. An integrated approach for robust
inventory-routing problem. Comput Ind Eng 2009;56(3):1095–112. inventory routing problem in a three-echelon distribution system. Int J Logist
[22] Huang SH, Lin PC. A modified ant colony optimization algorithm for multi- Syst Manag 2019;32(3–4):414–36.
item inventory routing problems with demand uncertainty. Transp Res E [52] Fardi K, Jafarzadeh Ghoushchi S, Hafezalkotob A. An extended robust approach
2010;46(5):598–611. for a cooperative inventory routing problem. Expert Syst Appl 2019;116:310–27.
[23] Yu Y, Chu C, Chen H, Chu F. Large scale stochastic inventory routing [53] Golsefidi AH, Jokar MRA. A robust optimization approach for the production–
problems with split delivery and service level constraints. Ann Oper Res inventory-routing problem with simultaneous pickup and delivery. Comput Ind
2012;197(1):135–58. Eng 2020;143:106388.
[24] Abdul Rahim MKI, Zhong Y, Aghezzaf EH, Aouam T. Modelling and solving the [54] Bertsimas D, Sim M. The price of robustness. Oper Res 2004;52(1):35–53.
multiperiod inventory-routing problem with stochastic stationary demand rates. [55] Feng Y, Liu Y, Chen Y. Distributionally robust location–allocation models of
Int J Prod Res 2014;52(14):4351–63. distribution centers for fresh products with uncertain demands. Expert Syst Appl
[25] Soysal M, Bloemhof-Ruwaard JM, Haijema R, Van Der Vorst JG. Modeling 2022;209:118180.
an inventory routing problem for perishable products with environmental [56] Delage E, Ye Y. Distributionally robust optimization under moment uncertainty
considerations and demand uncertainty. Int J Prod Econ 2015;164:118–33. with application to data-driven problems. Oper Res 2010;58(3):595–612.
[26] Micheli GJ, Mantella F. Modelling an environmentally-extended inventory rout- [57] Li R, Cui Z, Kuo YH, et al. Scenario-based distributionally robust optimization
ing problem with demand uncertainty and a heterogeneous fleet under carbon for the stochastic inventory routing problem. Transp Res E 2023;176:103193.
control policies. Int J Prod Econ 2018;204:316–27. [58] Cui Z, Long DZ, Qi J, Zhang L. The inventory routing problem under uncertainty.
[27] Kumar M, Kumar D, Saini P, Pratap S. Inventory routing model for perishable Oper Res 2023;71(1):378–95.
products toward circular economy. Comput Ind Eng 2022;169:108220. [59] Zhu S, Fukushima M. Worst-case conditional value-at-risk with application to
[28] Kleywegt AJ, Nori VS, Savelsbergh MW. Dynamic programming approximations robust portfolio management. Oper Res 2009;57(5):1155–68.
for a stochastic inventory routing problem. Transp Sci 2004;38(1):42–70. [60] Weskamp C, Koberstein A, Schwartz F, Suhl L, Voß S. A two-stage stochastic
[29] Hvattum LM, Løkketangen A, Laporte G. Scenario tree-based heuristics for programming approach for identifying optimal postponement strategies in supply
stochastic inventory-routing problems. INFORMS J Comput 2009;21(2):268–85. chains with uncertain demand. Omega 2019;83:123–38.
[30] Bertazzi L, Bosco A, Guerriero F, Lagana D. A stochastic inventory routing [61] Rockafellar RT. Convex analysis. Princeton: Princeton University Press; 1970.
problem with stock-out. Transp Res C 2013;27:89–107. [62] Manousakis E, Repoussis P, Zachariadis E, Tarantilis C. Improved branch-and-cut
[31] Bertazzi L, Bosco A, Laganà D. Managing stochastic demand in an inventory for the inventory routing problem based on a two-commodity flow formulation.
routing problem with transportation procurement. Omega 2015;56:112–21. European J Oper Res 2021;290(3):870–85.
[32] Juan AA, Grasman SE, Caceres-Cruz J, Bektaş T. A simheuristic algorithm for the [63] Coelho LC, Laporte G. A branch-and-cut algorithm for the multi-product
single-period stochastic inventory-routing problem with stock-outs. Simul Model multi-vehicle inventory-routing problem. Int J Prod Res 2023;1(23–24):7156–69.
Pract Theory 2014;46:40–52. [64] Boyd S, Boyd SP, Vandenberghe L. Convex optimization. Cambridge University
[33] Gruler A, Panadero J, de Armas J, Pérez JAM, Juan AA. Combining variable Press; 2004.
neighborhood search with simulation for the inventory routing problem with
stochastic demands and stock-outs. Comput Ind Eng 2018;123:278–88.

25