Allocation under a general substitution structure

European Journal of Operational Research 277 (2019) 492–506

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Production, Manufacturing, Transportation and Logistics

Allocation under a general substitution structure Yulia Tseytlin a, H. Sebastian Heese b,∗ a b

Institute for Supply Chain Management, EBS University, Burgstrasse 5, Oestrich-Winkel, 65375 Germany Poole College of Management, NC State University, Raleigh, NC 27695-7229, USA

a r t i c l e

i n f o

Article history: Received 8 January 2018 Accepted 27 February 2019 Available online 5 March 2019 Keywords: Inventory Substitution Allocation Transportation problem Monge sequences

a b s t r a c t We develop a novel solution approach for the single-period multi-product inventory allocation problem with general substitution structure. First, we construct an initial allocation sequence and present conditions for existence of a Monge sequence, implying optimality of a greedy solution. We prove that the greedy allocation along the initial sequence is optimal even if no Monge sequence exists, when in the solution all demand is satisfied and all inventory is used. For all other cases, we develop a correction algorithm and prove optimality of the resulting allocation. The worst-case computational complexity of our solution is superior to existing algorithms for the structurally related transportation problem.

1. Introduction The problem of managing inventory of substitutable products is important in many industries and, accordingly, has received wide attention in the academic literature (e.g., Shumsky & Zhang, 2009). The inventory problem with substitution generally contains two stages. First, before demand is known, a decision is made for each product regarding the quantity to produce or to order. Then, after demand is realized, the problem is to optimally allocate the inventory to that demand. In this paper, we focus on the second-stage allocation problem, that is, determining an optimal matching of a given inventory to a known demand. Specifically, we study a single-period multi-product allocation problem under a general substitution structure. The allocation problem, as a part of the inventory management problem under substitution, has been widely studied; Lang (2010) and Shin, Park, Lee, and Benton (2015) provide recent surveys of this literature. Out of this broad stream, the contributions by Hsu and Bassok (1999), Bassok, Anupindi, and Akella (1999), Chen, Yao, and Zheng (1999), Rao, Swaminathan, and Zhang (2004) and Bansal (2010) are most related to our work; these studies consider different variations of the single-period multi-product inventory problem and solve the allocation problem as a part of the general problem. In the following discussion, we focus on their results regarding the allocation problem.

∗

Corresponding author. E-mail addresses: [email protected] (Y. Tseytlin), [email protected] (H.S. Heese). https://doi.org/10.1016/j.ejor.2019.02.049 0377-2217/© 2019 Elsevier B.V. All rights reserved.

© 2019 Elsevier B.V. All rights reserved.

The above-mentioned papers assume a full downward substitution structure, that is, a quality-based ranking of products, where a higher-quality (higher-ranked) product can serve as substitute for all lower-quality (lower-ranked) products. They also assume constant or additive substitution costs. Bassok et al. (1999), Chen et al. (1999), Rao et al. (2004) and Bansal (2010) make the additional assumption that higher quality (better performance) is associated with both higher effective revenue (selling price + shortage costs) and higher salvage value. These papers, in different variations, propose greedy allocation algorithms which are optimal under these assumptions. Bassok et al. (1999) and Bansal (2010) use the concept of Monge sequences (cf. Burkard, 2007) in their studies. Hsu and Bassok (1999) relax the assumption of product ordering according to effective revenue and salvage value, and develop an efficient greedy algorithm to solve the allocation problem in a substitution problem with random yield. The full downward substitution structure assumed in these studies clearly does not capture many practical settings. Our research is motivated by an allocation problem observed during our collaboration with a leading manufacturer of carbon-based products, for a product portfolio consisting of cylindrical graphite electrodes, which differ in terms of grade (quality), length, and diameter. When allocating existing inventory to standing customer orders, if needed, products of superior grade and/or larger dimensions can be used to satisfy demands for lower-quality and/or smaller products. Since products differ along multiple dimensions, they cannot be ordered along a single dimension, as under the full downward substitution structure. In addition, a downward substitution structure implicitly assumes vertically differentiated products and thus does not apply when products differ along a

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

horizontal dimension. In this paper, we consider a general substitution structure and develop a new solution procedure for the resulting allocation problem. There exists research on substitution structures that are more general than full downward substitution (e.g., Bitran & Leong, 1992), but these studies provide heuristics for the overall substitution problem and not optimal solutions for the allocation problem. Consistent with our observations at the manufacturer or carbon-based products and with the literature described above, we assume that substitution costs are additive, which is structurally equivalent to considering constant substitution costs (e.g., Bassok et al., 1999; Rao et al., 2004). As noted in the literature (e.g., Bassok et al., 1999), the allocation problem can be viewed as a transportation problem, where the sources are the available product inventories, the destinations are demands for the products, and the arcs between sources and destinations represent feasible stock-to-demand allocations. The transportation problem (as a linear programming formulation) was first considered by Hitchcock (1941) and Kantorovich (1942), and it has been studied extensively since then, with particular focus on finding efficient solution algorithms (Alon, Cosares, Hochbaum, & Shamir, 1989). The transportation problem is a special case of the minimum cost problem (cf. Ahuja, Magnanti, & Orlin, 1989); to the best of our knowledge, the fastest general algorithm for these two problems is due to Orlin (1993). Our work closely relates to the notion of Monge sequences. Hoffman (1963) has shown that the transportation problem (and other optimization problems) can be solved efficiently by a greedy algorithm, if the problem cost coefficients fulfill a certain Monge property, or form a Monge sequence, named so due to an analogy with Monge’s observations from 1781 (Burkard, 2007). Since the study of Hoffman (1963), the Monge properties have been applied in a wide spectrum of discrete optimization problems and beyond; comprehensive reviews can be found in Burkard, Klinz, and Rudolf (1996) and Burkard (2007). In particular, Alon et al. (1989) and Shamir and Dietrich (1990) (see also Dietrich, 1990; Shamir, 1993) provide algorithms to construct Monge sequences or determine that none exists for general and restricted (where cost matrices have infinite entries) transportation problems, respectively. Queyranne, Spieksma, and Tardella (1998) consider transportation problems with infinite entries under the assumption that non-forbidden entries form a sublattice, implying that the infinite entries respect the Monge property; this assumption does not necessarily hold in our context. A recent study by Weiß, Knust, Shakhlevich, and Waldherr (2016) considers assignment problems with nearly Monge matrices, that is, matrices where finite entries satisfy the Monge property. The authors develop a linear-time algorithm to solve the problem, noting that their approach is not suitable for transportation problems. Our study is positioned between the two above-mentioned research streams: we consider a more general case of the allocation problem and a more special case of the transportation problem. Specifically, we develop a new solution procedure for the allocation problem with general substitution structure and additive substitution costs. We approach the problem as follows. First, we present an algorithm to construct a sequence of arcs (initial allocation sequence) and define conditions for it to be a Monge sequence; we prove that these conditions are both necessary and sufficient for existence of a Monge sequence in this problem. Existence of a Monge sequence facilitates the derivation of an optimal allocation for every inventory and demand instance in linear time using a simple greedy algorithm. We prove that the solution obtained by a greedy allocation along the initial sequence is optimal even in the absence of a Monge sequence, if it does not contain underages (unsatisfied demand) or overages (unused inventory). The main contribution of our paper is the development

493

of algorithms to correct the allocation for all other cases, and to prove optimality of the obtained solution. The correction is based on recognizing violation cases where underages and/or overages should be moved to other columns/rows; the modification is performed on each case separately, along the order of the associated marginal improvements in the objective function. Our solution approach thus is comprised of two parts. The first part – constructing an allocation sequence and determining the existence of a Monge sequence – is executed only once for any given cost structure. The second part – finding an initial allocation and correcting it when necessary – is executed for each supply and demand instance. We show that the computational complexity of our solution is superior to existing algorithms for more general problems. The remainder of the paper is organized as follows. In Section 2, we present notation, assumptions, the allocation problem formulation, and relevant definitions. In Section 3, we construct an initial sequence of arcs and define conditions for existence of a Monge sequence. We prove that a greedy allocation along this sequence is optimal, if it does not involve any underages or overages. For the case with underages or overages, in Section 4 we present correction algorithms and prove that they lead to an optimal solution to the allocation problem. An analysis of computational complexity and performance follows in Section 5. Finally, we discuss our results and ideas for future research in the concluding Section 6. Auxiliary results are provided in Appendix A, proofs of the results contained in the main part of the paper are given in Appendix B, some secondary algorithms are presented in Appendix C and details on the numerical study are given in Appendix D. 2. The model In this section we present notation and assumptions and introduce the graph and matrix representations of the problem. We then formulate the allocation problem and provide the relevant definitions. We consider a single-period model with N product types. Unsatisfied demand (underage) is lost and unused stock (overage) is salvaged. The selling price, the shortage penalty cost, and the salvage value of product i are denoted by pi , π i and si , respectively. All costs and revenues are linear, proportional to quantities. Satisfying demand of product j with inventory of product i incurs an additional substitution cost bij ; bi j = ∞, when substitution of j with i is not possible. As mentioned above, our problem can be viewed as the transportation problem; thus it can be represented as bipartite directed graph with the source (inventory) nodes on the left side and the destination (demand) nodes on the right side. Equivalently, the problem can be represented with a weight (cost) matrix W where rows are sources and columns are destinations. The arc (matrix entry) (i, j) between node i on the left and node j on the right, implies satisfying demand of product j with inventory of product i. As substitution is not always possible for all product pairs, this is a restricted transportation problem (cf., Shamir & Dietrich, 1990). We introduce two additional nodes, a super-source U (matrix row) and a super-sink O (matrix column), to account for possible underage and overage, respectively, thus ensuring feasibility. We use E and F to denote the set of admissible and forbidden arcs, respectively, and E as the set of admissible arcs in the N × N matrix, without the nodes U and O. We also define e = |E |, f = |F |, and e = |E  |. To simplify exposition, we formulate the problem as a minimization problem, omitting the (constant) revenue associated with total demand and accounting for lost revenue as part of the underage costs. Starting with a given inventory position, we also omit production or order costs without any loss of generality. We define underage cost of product i as ui = pi + πi and overage cost of product i as oi = −si . The costs on the arcs (matrix entries) then

494

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

are as follows: for substitution ∀(i, j) ∈ E , wi j = bi j , for forbidden arcs ∀(i, j) ∈ F, wi j = ∞, for underage ∀ j ∈ {1, . . . , N}, wU j = u j , for overage ∀i ∈ {1, . . . , N}, wiO = oi , and wUO = 0. We make two assumptions regarding the substitution costs. First, we assume that substitution costs are additive, that is, the cost associated with using product i to substitute for product j can be written as a sum of two independent costs, where the first depends only on i and the second only on j, namely ∀(i, j ) ∈ E  : bi j = αi + β j (cf., Bassok et al., 1999). Second, consistent with the related literature (e.g., Assumption (9) in Hsu and Bassok, 1999), we assume that, whenever feasible, substitution is more profitable than simultaneous overage and underage. (We assume this inequality to hold strictly in order to avoid unnecessary expositional complications from the otherwise arising ambiguity.) Assumption 1. u j + oi > αi + β j

∀i, j ∈ {1, . . . , N}

We note that the assumption is required to hold for all couples i, j and not only for the admissible arcs; that is, if (i, j) ∈ F and there is a way to use the stock of i and satisfy the demand on j, we assume this to be more beneficial than having overage on i and underage on j. We note that α i < ∞, ∀i and β j < ∞, ∀j, as the substitution costs on the admissible arcs are finite (for (i, j) ∈ F: bi j = αi + β j ). Continuing with the notation, we use q to denote the vector of supplies, where qi is the quantity of product i in stock. Similarly, we use d to denote the demand vector, where dj is the quantity of product j demanded. Finally, let Xij denote the allocation on arc (i, j), ∀i ∈ {1, . . . , N} ∪ {O} and ∀ j ∈ {1, . . . , N} ∪ {U }. Note that Xi j = 0 for all forbidden arcs (i, j) ∈ F. Using C to denote total costs, the optimization problem we consider is

min C = Xi j



wi j · Xi j ,

(1)

(i, j )∈E

subject to:



Xi j = qi

∀i = 1, . . . , N,

(2)

Xi j = d j

∀ j = 1, . . . , N,

(3)

j∈{1,...,N}∪{O}

 i∈{1,...,N}∪{U }

Xi j ≥ 0

Xi j = 0

∀(i, j ) ∈ E, ∀(i, j ) ∈ F.

3. Sequence construction, greedy allocation, and sufficient conditions for optimality In this section, we first construct an allocation sequence (in Algorithm 1) and then derive a necessary and sufficient condition for Monge sequence existence. If the allocation sequence is a Monge sequence, greedy allocation (as in Algorithm 2) is optimal for all possible demand and supply vectors. Otherwise, we show that the greedy allocation is nevertheless optimal if the resulting allocation does not contain underages or overages. Algorithm 1 rearranges the weight matrix W so that the columns and the rows are ordered according to u j − β j and oi − αi , respectively, both in non-ascending order. For simplicity of exposition, in the following we assume that all products differ in terms of overage and underage costs. (The case where some products have the same underage or overage costs is addressed with the alternative Algorithm 7, in Appendix C.) In the resulting matrix, the upper left entry (1,1) has maximal underage and overage costs, and should have the highest priority for allocation; entries (2,1) and (1,2) should receive allocation next; and so on. For the remainder of the paper, we refer to arcs according to their position in the newly arranged matrix; for example, we use (1,1) to denote the arc with maximal underage and overage costs. The sequence of arcs S constructed in Algorithm 1 will form the basis for the subsequent inventory allocation and the correction algorithms. In the process of constructing S, we also create two sub-sequences of S: sequence Se of admissible arcs (used for the allocation) and sequence Sf of forbidden arcs (used for the correction). Algorithm 1 Constructing sequence S. 1. Initialization and weight matrix rearrangement: (a) S, Se and S f are empty ordered sets. (b) Sort columns j ∈ {1, . . . , N} by u j − β j , in non-ascending order. (c) Sort rows i ∈ {1, . . . , N} by oi − αi , in non-ascending order. 2. Filling S, Se , and S f (a) For k = 2 to 2N: For i = min(k − 1, N ) down to max(1, k − N ): Append (i, k − i ) to S. If (i, k − i ) ∈ E, append it to Se ; else append it to S f . (b) For i = 1 to N: add (i, O ) to S. (c) For j = 1 to N: add (U, j ) to S.

(4)

(5)

The objective (1) is to minimize total costs. Constraint (2) ensures inventory balance for each product type in stock and constraint (3) ensures demand balance for each product type demanded. We next present the definition of Monge sequences: for an m × n matrix W, a permutation S of the indices (i, j) is called Monge sequence if, for every 1 ≤ i, r ≤ m, 1 ≤ j, s ≤ n, whenever (i, j) precedes both (r, j) and (i, s), then wi j + wrs ≤ wis + wr j (cf., Hoffman, 1963). In our paper, we use the alternative definition of Shamir and Dietrich (1990), which is adjusted specifically for restricted transportation problems: a permutation S on the set of admissible arcs E is called a Monge sequence if, for every (i, j) ∈ E and for every r = i, s = j, the following two conditions are both satisfied: (MS1) If (r, j), (i, s), (r, s) ∈ E and wi j + wrs > wis + wr j then either (r, j) or (i, s) precede (i, j) in S. (MS2) If (r, j), (i, s) ∈ E and (r, s) ∈ F then either (r, j) or (i, s) precede (i, j) in S.

We next analyze conditions under which Se is a Monge sequence. As defined in Section 2, the Monge sequence conditions need to hold for every four arcs (quadruple), at least three of which are admissible. Such quadruples can be divided into two groups; either all four arcs are admissible (condition (MS1)) or one arc is forbidden (condition (MS2)). We can see that, due to the additive substitution costs, for every four arcs (i, j), (r, j), (i, s), (r, s) ∈ E it holds that

wi j + wrs = bi j + brs = αi + β j + αr + βs = bis + br j = wis + wr j , (6) that is the cost of satisfying demand of product j with product i and demand of product s with product r is equal to the cost of satisfying demand of product j with product r and demand of product s with product i. Hence, condition (MS1) is satisfied in equality for sequence Se , for every quadruple of admissible arcs that do not include overage or underage nodes. As Se is ordered according to underage and overage costs, condition (MS1) is satisfied for the overage and underage arcs as well. These two statements are formally proved in Lemma 1 in Appendix A. For our problem, we reformulate condition (MS2) as condition (C1) below (the

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

495

equivalence of these two conditions is proven in Lemma 2 in Appendix A). Condition (C1): For any quadruple of arcs with (r, s) ∈ F and (r, j), (i, s), (i, j) ∈ E, the following two conditions do not hold simultaneously: (C1.1) u j − β j > us − βs (or s > j), (C1.2) oi − αi > or − αr (or r > i).

Fig. 1. Example (U) – underage correction.

Proposition 1. Se created in Algorithm 1 is a Monge sequence if and only if condition (C1) holds. Proposition 1 stems from Lemmas 1 and 2. Theorem 1 provides a more general structural property, establishing condition (C1) as a necessary and sufficient condition for existence of a Monge sequence in the setting considered in this paper.

Fig. 2. Example (O) – overage correction.

Theorem 1. In the allocation problem with additive substitution costs, a Monge sequence exists if and only if condition (C1) holds. The following Algorithm 2 allocates supplies q to demands d along sequence Se in a greedy fashion, that is, each arc (i, j) receives the minimum of supply at node i and demand at node j. Algorithm 2 is a variation of a well-known algorithm used to obtain initial solutions for transportation problems (cf. Alon et al., 1989; Shamir & Dietrich, 1990). Algorithm 2 Greedy allocation along sequence Se . For every arc (i, j ) in Se : – – – – –

ai j = min(qi , d j ) Xi j = ai j qi = qi − ai j d j = d j − ai j Move to next (i, j )

If Se is a Monge sequence, the allocation obtained through Algorithm 2 is optimal for all supply and demand vectors (Shamir and Dietrich, 1990, Theorem 2.2). We next consider an alternative case where Se is not a Monge sequence. In Theorem 2, we prove that the greedy allocation obtained through Algorithm 2 is optimal, if it contains no underage or overage. For solutions with underage and/or overage, in the following Section 4, we provide correction algorithms and show that this correction leads to an optimal solution. Theorem 2. If Algorithm 2 results in a solution X with XiO = 0, ∀i ∈ {1, . . . , N} and XU j = 0, ∀ j ∈ {1, . . . , N}, then X is optimal. 4. Correction of overages and underages If Se is not a Monge sequence and the solution from Algorithm 2 contains underages or overages, the resulting allocation might not be optimal. In Section 4.1, we characterize different sub-optimalities or violation cases. In Section 4.2, we present Algorithm 3 that corrects a single violation case. In Section 4.3, we determine the order in which the violation cases should be addressed and present the overall correction Algorithms 4 and 5. With Theorem 3, we prove that the corrected allocation after Algorithms 4 and 5 is optimal. 4.1. Violation cases We distinguish among three ways in which a violation of (C1) can lead to suboptimal overage or underage allocation in the solution obtained by Algorithm 2. In the following, we present one example for each of the three violation types.

Fig. 3. Example (UO) – underage and overage correction.

Example (U): Consider the following four arcs (r, s) ∈ F, (r, j), (i, s), (i, j) ∈ E that violate (C1), namely r > i, s > j, and (i, j) precedes (r, j) and (i, s) in Se . There is an additional arc (r, k) ∈ E such that us − βs > uk − βk . Assume d j = ds = dk = qi = qr = 1 and all other demands and supplies are 0. The allocation received from Algorithm 2, Xi j = XUs = Xrk = 1, is suboptimal, as underage on k would be preferable to underage on s. Changing the allocation from (i, j) to (i, s) and (r, j) and moving underage to column k (Xis = Xr j = XUk = 1) leads to an optimal solution. The suboptimal and the corrected allocations are given in Fig. 1 (here and in the following figures the ∞-entries imply forbidden arcs and the empty cells imply admissible arcs with zero allocation). Example (O): Given the same base example as above, now assume there is an additional arc (l, s) such that or − αr > ol − αl . Assume that d j = ds = qi = qr = ql = 1 and all other demands and supplies are 0. The allocation received from Algorithm 2, Xi j = XrO = Xls = 1, is suboptimal, as overage on l would be preferable to overage on r. Changing the allocation from (i, j) to (i, s) and (r, j) and moving overage to row l (Xis = Xr j = XlO = 1) leads to an optimal solution (see Fig. 2). Example (UO): The third violation type contains both overage and underage, that is, there exists an arc (r, s) ∈ F where XrO > 0 and XUs > 0. Continuing our example from above, such a case can happen when d j = ds = qi = qr = 1 and all other demands and supplies are 0. The allocation received from Algorithm 2, Xi j = XUs = XrO = 1, is suboptimal due to Assumption 1. Changing the allocation from (i, j) to (i, s) and (r, j) (Xis = Xr j = 1) leads to an optimal solution (see Fig. 3). Thus, violations of (C1) can lead to three violation types, where it might be beneficial to change the allocation: (U) move underage from column j to column s, where u j − β j > us − βs (i.e., it is more beneficial to have underage on s than on j); (O) move overage from row i to row r, where oi − αi > or − αr (i.e., it is more beneficial to have overage on r than on i); or (UO) simultaneously move underage from column j and overage from row i (which is beneficial due to Assumption 1). There might be multiple violation cases of each type; we analyze the order of the cases correction in Section 4.3. Each violation case can be attributed to a forbidden arc – (r, s) in our examples (each forbidden arc can have more than one case). We thus characterize each violation case by the associated forbidden arc and the violation type. In addition, cases of type (U) are characterized by the target column for moving the underage (s), and cases of type (O) by the target row for moving the overage (r). This way we denote each case by ((i, j), T, x), where (i, j) is the

496

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

associated forbidden arc, T the type of the case (U/O/UO) and x the target column for type (U) or the target row for type (O). When the whole specification is not required, for expositional simplicity, we might refer to a case through a single index denoting the order of the case in the correction sequence (specified in Section 4.3). Before elaborating the solution process in Sections 4.2 and 4.3, we provide an intuitive explanation. In the examples above (consider in particular Example (O)), the suboptimal allocation (on arc (r, O)) is moved to the preferred arc (arc (l, O)), while ensuring that the allocation sums on each row and each column remain unchanged (to satisfy constraints (2) and (3)). In this simple example, the move takes several steps: from (r, O) to (r, j), from (r, j) to (i, j), from (i, j) to (i, s), from (i, s) to (l, s), and finally from (l, s) to (l, O). If we represent the matrix as a graph, this can be seen as a path that connects (r, O) with (l, O) (or (l, s)), where the allocation is moved along the path from each odd node to the succeeding node, and each pair of succeeding nodes is either in the same row or in the same column, thus preserving the sums in the rows/columns. In more complex examples, the allocation can be moved along multiple paths, each comprised of multiple nodes. This allocation movement can be seen as a flow, from a source to a terminal, and we search the paths until the maximal possible flow (the minimum of the allocation on the source and on the terminal) is reached or all paths are exhausted, such that it is not possible to correct the case further. We measure the profitability of a case by the per-unit marginal cost improvement gained from its potential correction, that is, the difference in the objective function obtained by reallocating one unit. As there might be multiple cases, and correction of one case can affect the other, we need to proceed in such an order that the correction of a case will never impede the correction of a more beneficial case. We show that this approach leads to an optimal solution.

4.2. Correction of a single violation case In this section, we describe the correction process of a single violation case l = ((i, j ), T , x ). For our next algorithms, we represent the weight matrix as a graph G, whose nodes correspond to the admissible entries (as in Alon et al., 1989). Two nodes, including the overage and the underage entries, are connected by an edge if they are in the same row or in the same column. For each case l, following the definition in Itai and Shiloach (1979), we build a directed flow network Nl = (Gl , Sl , Tl , cl ), where Nl consists of the directed graph Gl (a subgraph of G), the source Sl , the terminal Tl and the capacity function cl on the graph edges. The source and the terminal represent the entries the allocation on which we aim to modify, and the capacities are determined according to the current allocation. On this network we find a maximal flow function from Sl to Tl , according to which the allocation is modified. The formal procedure for correcting a single violation case is presented in Algorithm 3. We start constructing the graph from node (i, s), where s = x, if T = U, and s = O otherwise. The final node is (r, j), where r = x, if T = O, and r = U otherwise. The choice to start building the graph from a row rather than from a column is arbitrary and without loss of generality. We find paths (sequences of nodes) from the source to the terminal, along which the allocation can be moved. For each node in the graph we maintain the depth d, which represents the distance from the source; the exact depth is not important for the algorithm, only its parity, which we denote by p(d). When constructing a path, each node with odd depth can be connected to the nodes in the same row, and each node with even depth to the nodes in the same column with positive allocation. For each odd node, the capacity on incoming and outgoing arcs equals the allocation on this entry.

Our algorithm is a variation of the augmenting path maximal flow algorithm of Ford and Fulkerson (1956), which finds a path from the source to the sink at each iteration and augments a maximal possible flow along this path. This is done by recognizing a bottleneck arc (arc with minimal capacity) of the path and increasing the flow according to the bottleneck capacity. Following a similar logic, we push a maximal possible flow through each path while constructing it, until the total maximal possible flow is reached or all paths are exhausted, implying that the violation case cannot be corrected further. At each iteration k, we construct a path Pk from the source to the terminal, find bottleneck edges and alter the allocation accordingly; the allocation from each odd node is moved to the subsequent even node along the path. Next, the bottlenecks are deleted from the graph, thus dividing the path Pk into several partial paths. We save these partial paths, in order not to search for them again, while maintaining for each node v a variable pv that indicates the partial path that v belongs to (at each iteration each node can belong to at most one partial path). At the next iteration, we aim to connect the node before the first bottleneck with the node after the last bottleneck, this way creating a new path from the source to the terminal. The path construction procedure is inspired by the maximal flow algorithm of Itai and Shiloach (1979). Algorithm 3 does not change any arcs involving underage or overage other than those associated with the targeted violation case. In Lemma 4 in Appendix A, we show that correcting a violation case more than done by Algorithm 3, without changing allocations on other overage or underage arcs, is not possible. 4.3. Correction algorithm and optimality proof Next we determine in what order to address the different violation cases. As the correction of one case can affect the correction of the following cases, we address the cases in order of their profitability. Denote by (i,j),T,x the per-unit marginal cost improvement from the correction of case ((i, j), T, x), that is, the difference in the objective function obtained by reallocating one unit with Algorithm 3. This marginal cost improvements for the three violation types then are: (U) (i, j ),U,s = (u j − β j ) − (us − βs ), (O) (i, j ),O,r = (oi − αi ) − (or − αr ), (UO) (i, j ),UO = (u j + oi ) − (αi + β j ). We correct the allocation in Algorithm 4, along the following order: first, cases of type (UO) are considered, according to their order in sequence Sf (created in Algorithm 1). Then cases of type (U) are addressed, in such an order that case ((i, j), U, l) is chosen first if j is the smallest-indexed column with underage and l is the largest-indexed column such that (i, l) has a positive allocation. Finally, cases of type (O) are corrected, in such an order that case ((i, j), O, l) is chosen first if i is the smallest-indexed row with overage and l is the largest-indexed row such that (l, j) has a positive allocation. (Lemma 14 in Appendix A shows that, after correcting all (UO) cases, the corrections of (U) and (O) cases are independent, hence it does not matter whether we start with type (U) or (O).) Note that correction of each case may affect the existence of subsequent cases. For example, consider a forbidden arc (i, j) with three violation cases: (UO), (U) and (O). After the correction of case (UO) with Algorithm 3, it is possible that (1) no overage on i and no underage on j are left, (2) only overage on i is left, (3) only underage on j is left, or (4) both overage on i and underage on j are left, and it is not possible to transfer more flow from row i to column j (the case becomes not-correctable). In (2) only case (O) is left to correct, in (3) only case (U), and in (1) and (4) neither. This way, after addressing a case, the subsequent cases may not exist anymore or may be not-correctable (in Algorithm 4, set NC

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

497

Algorithm 3 Correcting a single violation case l = ((i, j ), T , x ). 1. Initialization: (a) Initialize Gl – an empty directed graph. Mark all nodes in graph G as not-visited, with empty p. Set counter k = 1. (b) Initialize the starting node (i, s ) and the final node (r, j ): i. If T = U, x = s: start with node (i, s ) and finish with node (U, j ). ii. If T = O, x = r: start with node (i, O ) and finish with node (r, j ). iii. If T = UO: start with node (i, O ) and finish with node (U, j ). (c) Connect the source Sl to (i, s ). Set dis = 1 and mark (i, s ) as visited. Connect (r, j ) to the terminal Tl . (d) Set the capacity on the edge (Sl , (i, s )) to Xis and on the edge ((r, j ), Tl ) to Xr j . (e) Initialize two paths P s = (Sl , (i, s )) and Pt = ((r, j ), Tl ), update p(i,s ) and p(r, j ) accordingly, and set two nodes v1 = (i, s ) and v2 = (r, j ). (f) Define the maximal flow to be transferred through the network as m f = min(Xis , Xr j ). 2. Find the next path Pk (if none exists, exit the algorithm): (a) Pk = P s . Assume v1 = (z, y ). (b) If dv1 is even, connect v1 to a not-visited node in the same column (not including (U, y )) with positive allocation, i.e., to w = (q, y ) ∈ E  q = z : Xqy > 0, if (v1, w ) is not deleted. (c) If dv1 is odd, check whether (z, j ) ∈ E. If so, connect v1 to w = (z, j ) and go to Step 3. Otherwise, connect v1 to a not-visited node in the same row (not including (z, O )), on which the allocation could be increased, i.e., to w = (z, q ) ∈ E  q = y : dq − Xzq − XUq > 0, if (v1, w ) is not deleted. (d) If there is no such w: if v1 = (i, s ), stop: no further path exists. Otherwise, set v1 to be the previous node u in Pk , after deleting v1 and the edge (u, v1 ) from Pk and from G. (e) If pw is empty (w does not belong to any partial path), add edge (v1, w ) to Pk and to Gl , set dw = dv1 + 1, mark w as visited, set v1 = w and go to Step 2b. (f) Otherwise (w belongs to some partial path, say P Pι , which begins with v1ι and ends with v2ι ): i. If p(dv1 ) = p(dv1ι ) (different depth parity),and w = v2ι , set w to be the next node in P Pι . ii. If P Pι = Pt , go to Step 3 (the new path is found). iii. Connect Pk to P Pι with edge (v1, w ). Add edge (v1, w ) to Gl , set dw = dv1 + 1, mark nodes in P Pι from w to v2ι as visited and set v1 = v2ι . Shorten P Pι by setting v2ι to be the node before w. Go to Step 2b. 3. Save the edges from v2 to w as a partial path, remove them from Pt , and add the remaining edges of Pt to Pk . Determine capacities of the new edges (that were added to the graph during this iteration): for every node v = (z, y ) with odd depth (where Xzy > 0), the capacity on the incoming and outgoing edges of v is Xzy . Mark the residual capacity of each such edge a, res(a ), to be equal to the capacity. Delete the opposite directions of the new edges in the path: for each (v, w ) ∈ Pk , delete (w, v ) from G. 4. Find the bottlenecks, i.e., edges where the residual capacity equals mina∈Pk res(a ). Set fk = mina∈Pk res(a ). Increase the flow along Pk by fk units. If the total flow through the network has reached m f , stop. Otherwise, reduce the residual capacity of the Pk edges by fk , update v1 and v2 to be the starting node of the first bottleneck arc and the ending node of the last bottleneck arc respectively. 5. Modify the allocation: for each node v = (z, y ) in Pk : (a) Denote by out (v ) and in(v ) the sets of all outgoing and incoming arcs.  (b) If dv is odd: Xzy = Xzy − a∈out (v ) fk (a ).  (c) If dv is even: Xzy = Xzy + a∈in(v ) fk (a ). 6. Delete the bottlenecks and the related (saturated) nodes from the graph. This divides Pk into several paths. Update P s and Pt : P s from the source to v1 and Pt from v2 to the terminal, and save all other partial paths. Update p’s of the nodes along the path. Mark all nodes on the partial paths and on Pt as not-visited. 7. k = k + 1, go to Step 2.

ensures that not-correctable cases are not addressed). To illustrate not-correctable cases, consider example (U) on page 7. If (r, j) ∈ F, the algorithm would attempt to correct the allocation but would not find a path. The case thus is not-correctable, and the current allocation is optimal. However, it might also happen that new violation cases are created when the allocation changes. As allocations are always moved inside the same column or row (i.e., the sum in each column/row remains the same), the existing underage cases may be shifted to other rows and the existing overage cases to other columns. For example, a correction of some violation case might involve moving allocation to (r, j). If there exists an s < j with XUs > 0 and (r, s) ∈ F, a new violation case ((r, s), U, j) is created. We note that no new correctable (UO) cases are created throughout Algorithm 4 (cf. Lemma 8 in Appendix A), hence these cases can be sorted in advance, along the order of Sf . For the (U) and (O) cases, due to their dynamicity, at each iteration the algorithm searches for the next violation case in a way that the resulting order adheres to the criteria specified above (rather than creating a correction sequence and update it after each case correction).

The next result shows that, in Algorithm 4, a correction of a case early on in the sequence will never impede the correction of a more beneficial case. Proposition 2. Consider any two cases a and b that are corrected in Algorithm 4. If case a is corrected before case b, then either a ≥ b or the order of these two cases in the sequence does not matter. During the correction process in Algorithm 4, it could happen that the solution contains arcs (i, j) ∈ E, such that either XUj > 0 and Xis > 0, s > j, or XiO > 0 and Xrj > 0, r > i. This is clearly suboptimal and is corrected in Algorithm 5 below. The following theorem proves that the given correction algorithms lead to an optimal allocation. Theorem 3. The solution obtained after Algorithms 4 and 5 is optimal. The proof of optimality of the obtained solution utilizes several auxiliary Lemmas, which are all provided in Appendix A. In Lemma 5, we show that it is not possible to improve a case

498

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

Algorithm 4 Correcting the allocation. 1. Initialize NC – an empty arc set. 2. While there are unaddressed violation cases: (a) Identify the next violation case ((i, j ), T , x ): i. While there are unaddressed (UO) cases, consider the next ((i, j ), UO ) case along the order of S f . Go to Step 2b. ii. While there are unaddressed (U) cases, consider the next ((i, j ), U, s ) case such that j is minimal among all columns with XU j > 0, s > j is maximal with Xis > 0 and (i, j ) ∈ / NC. Go to Step 2b. iii. While there are unaddressed (O) cases, consider the next ((i, j ), O, r ) case such that i is minimal among all rows with XiO > 0, r > i is maximal with Xr j > 0 and (i, j ) ∈ / NC. (b) Attempt to correct case ((i, j ), T , x ) with Algorithm 3. (c) Update NC: i. For convenience, mark r = U and s = O, if T = UO; r = U and s = x, if T = U; r = x and s = O, if T = O. ii. If, after the correction, Xr j > 0 and Xis > 0, then add (i, j ) to NC.

Algorithm 5 Correcting overage and underage vs. admissible arcs. 1. Underage: For j = 1 up to N − 2, such that XU j > 0: For s = N down to j + 1, such that u j − β j > us − βs and ∃i : (i, j ) ∈ E, Xis > 0: i. δ = min(XU j , Xis ) ii. Xi j = Xi j + δ iii. Xis = Xis − δ iv. XUs = XUs + δ v. XU j = XU j − δ 2. Overage: For i = 1 up to N − 2, such that XiO > 0: For r = N down to i + 1, such that oi − αi > or − αr , ∃ j : (i, j ) ∈ E, Xr j > 0: i. δ = min(XiO , Xr j ) ii. Xi j = Xi j + δ iii. Xr j = Xr j − δ iv. XrO = XrO + δ v. XiO = XiO − δ

further after having addressed other cases later in the sequence. Lemma 6 shows that throughout the solution process, there exists no simultaneous underage and overage for any admissible arc. In Lemma 7, we show that the results from Lemmas 5 and 6 continue to hold after Algorithm 5, and Lemma 8 proves that Algorithm 5 creates no additional correctable cases. In Algorithm 3 we do not change O/U arcs besides the ones related to the considered case; in Lemma 9 we show that the final solution cannot be improved by changing these arcs. 5. Analysis of complexity and computational performance In this section, first, in Section 5.1, we analyze the complexity of our solution approach. Then, in Section 5.2, we conduct a brief numerical study to benchmark the computational performance of our algorithm against a well-established algorithm for the more general minimum cost flow problem. 5.1. Complexity analysis As discussed previously, the solution is comprised of two parts. In the first part, we construct an allocation sequence and deter-

mine whether a Monge sequence exists. This part is executed only once for any given cost structure. The second part, that is, finding an initial allocation and correcting it (if required), is executed for each supply and demand instance. We first analyze the worst-case complexity of operations performed once for each cost structure: creating an allocation sequence with Algorithm 1 and checking whether a Monge sequence exists through Condition (C1). For every forbidden arc (r, s) (with r > 1, s > 1) and for every admissible arcs (i, j) with i < r, j < s, such that u j − β j > us − βs and oi − αi > or − αr , we check whether (i, s) and (r, j) are admissible. If so, then (C1) does not hold and a Monge sequence does not exist. If (i, s) ∈ F, we do not need to check other forbidden arcs in column s versus (i, j) and other admissible arcs in row i versus (r, s); if (r, j) ∈ F, we do not need to check other forbidden arcs in row r versus (i, j) and other admissible arcs in column j versus (r, s). This procedure is formalized in Algorithm 6 in Appendix C. If (C1) holds and there are different products with the same underage and/or overage costs, an arbitrary ordering of these products might cause a violation of (MS2). The procedure of constructing the allocation sequence to avoid this violation, as an alternative to Algorithm 1, is presented in Algorithm 7 in Appendix C. Proposition 3. The total complexity of single-time operations (Algorithms 1, 6 and 7) is O (N 2 ). Next we evaluate the complexity of determining an optimal allocation for given supply and demand vectors. This consists of two parts - finding an initial allocation with Algorithm 2, which is optimal if a Monge sequence exists or if no underage/overage is left, and modifying the allocation otherwise (as presented in Section 4). Before formalizing the results of the complexity evaluation in Proposition 4, we discuss some related issues, regarding Algorithm 3. We note that several steps can be added to the algorithm in order to enhance its efficiency (for clarity these additions are not stated in the algorithm itself). First of all, we could add a pre-processing step that orders columns and rows by entries with/without allocation; columns/rows with zero demand/supply should be omitted from the graph. Throughout the algorithm, sets of columns to be connected to an odd node at each row and similarly, sets of rows to be connected to an even node at each column could be maintained, to be used at Steps 2b and 2c. Additional efficiency improvements can be achieved if, before running Algorithm 3 for each violation case, we dispose of quadruples where all four arcs have positive allocations. Formally, if there exist four arcs (i, j), (r, j), (i, s), (r, s) ∈ E (r, i, j, s ∈ {1, . . . , N}) such that Xij > 0, Xis > 0, Xrj > 0 and Xrs > 0, we modify the allocation in a way that at least one of the arcs in the quadruple has zero allocation. In the solution obtained after Algorithm 2, there are no such quadruples, as the algorithm assigns to each arc the maximum possible before continuing to the next arc. However, modifying the allocation during case corrections can create quadruples where all four arcs have positive allocations. We modify the allocation on such arcs by moving the smallest allocation to the arc in the same row/column, while keeping the sum in each row and in each column the same. Due to the additivity assumption in (6), such a change does not affect the objective function. This procedure is formalized in Algorithm 8, which is contained in Appendix C. Proposition 4. The total complexity of finding an optimal solution after constructing the allocation sequence is: •

•

O (N ), if a Monge sequence exists or if no underage/overage is left (Algorithm 2); O (N 3 ), in all other cases (Algorithms 4 and 5, including Algorithms 3 and 8).

To provide some intuition as to why the complexity of our solution approach is O (N 3 ), we note that the maximal number of

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

steps to correct a single case is O (N 2 ) and the maximal number of cases to address is O (N 2 ). However, if each case takes the maximal number of steps to correct, then at most O (N ) cases are addressed. Alternatively, if each of the O (N 2 ) cases are addressed, then the number of steps for each case cannot exceed O (N ). In conclusion, the total complexity of single-time operations is O (N 2 ). If a Monge sequence exists or if the initial allocation contains no overages or underages, the complexity of solving the problem for a given demand and supply instance then is O (N ). If neither of these cases applies, the complexity of each run is O (N 3 ). 5.2. Numerical analysis of computational performance In this section, we summarize some insights from a brief numerical study conducted to benchmark the computational performance of our algorithm against a well-established algorithm for the minimum cost flow problem. Specifically, we compare the running time of our algorithm with the performance of Goldberg’s CS2 algorithm, which is one of the leading algorithms for the minimum cost flow problem. CS2 is the second version of the well-known cost-scaling algorithm of Goldberg (1997), which is a practical implementation of the cost-scaling push-relabel method. CS2 is highly robust and reliable, and it has served as a benchmark for solving minimum cost flow problems for a long time; it has also been found efficient and robust in a recent experimental evaluation of minimum cost algorithms (Kovács, 2015). We use the C++ code of CS2 from Ababei (2009), adjusted to our data structure. We also code our algorithm in C++. Both algorithms are compiled in Microsoft Visual Studio Express (version 12.0) and run on an Intel Core i7-3520M CPU with 2.9 gigahertz. We consider different scenarios in terms of problem size, cost structure and supply/demand instances. We vary the number of products N from 10 0 0 to 70 0 0, with a step size of 10 0 0. For each scenario, we independently sample overage costs oi and underage costs ui for all products from a uniform distribution with either low (20,40) or high (50,150) support. The values for α and β are drawn from a uniform distribution over (1,20). Each arc (i, j ), i, j ∈ {1, . . . , N}, i = j, is randomly assigned a 0 (forbidden arc) or a 1 (admissible arc), each with equal probability; arcs (i, i) are always set as admissible. Supply and demand of each product are drawn from a Normal distribution, with either a low (200) or a high (300) mean, and with either a low (10%) or high (50%) coefficient of variation; we truncate the normal distribution, re-sampling every realization smaller than zero or larger than twice the mean. Varying these parameters in a full factorial fashion and creating 10 scenarios for each configuration, we run a total of 4,480 scenarios. Overall, our algorithm performs the computations in about half the time needed by the benchmark CS2 algorithm (51.1% reduction in computational time). While the relative improvement in computational time varies across the individual scenarios (from a minimum of 4.2% to a maximum of 81.3%), the improvement is fairly consistent for the vast majority of the scenarios, with the 10th percentile at 36.8% and the 90th percentile at 64.0%. We observe that changes in the parameters do not significantly affect the relative improvement in computational time, besides slightly higher average improvements for scenarios with a larger number of products and for scenarios with a mismatch between demands and supplies in terms of average volume. (Appendix D provides some additional details on the results of our study.) We note that our algorithm would be running faster in settings with a given cost structure, where allocations are repeatedly determined for different supply and demand instances (e.g., a manufacturer deciding on product allocation every week). In such a case, the first part of our solution only needs to be executed once for all scenarios, and not for every scenario as assumed in the computational study.

499

6. Discussion and conclusion We study a single-period multi-product allocation problem with additive substitution costs. We extend prior research to a setting with a general substitution structure and suggest a novel solution approach. First, we construct an allocation sequence and present conditions for existence of a Monge sequence. If a Monge sequence exists, a greedy allocation along the allocation sequence is optimal. Even if no Monge sequence exists, we prove that the solution obtained with a greedy algorithm is optimal when it does not contain underages and overages. For all other cases we provide a correction algorithm and prove optimality of the resulting allocation. The correction procedure is based on recognizing violation cases where underage and/or overage allocations should be changed. The modification is performed on each case separately, in the order of the associated marginal improvements in the objective function. Our solution approach thus is comprised of two parts; while the first (sequence construction) is executed only once for each problem cost structure, the second (allocation and correction) needs to be executed for any given supply and demand instance. We show that the complexity of the first part is O (N 2 ), which is better than the complexity of the algorithm in Shamir (1993), that constructs a Monge sequence (or determines that none exists) for a restricted transportation problem in O (N 3 · log(N )) steps. The complexity of the second part is O (N ), if a Monge sequence exists or if the initial (greedy) allocation does not contain underages or overages; if a correction is required, the complexity increases to O (N 3 ), which still outperforms the O (N 3 · log(N )) complexity of the fastest known algorithm for the general transportation problem, presented in Orlin (1993). We show that our solution outperforms the existing solutions not only theoretically but also practically by conducting a numerical study and comparing our running times to those of the well-established CS2 algorithm of Goldberg (1997). The full downward substitution structure, as considered in Bassok et al. (1999), Chen et al. (1999), Rao et al. (2004), Bansal (2010) and Hsu and Bassok (1999), is a special case of the general structure analyzed in this paper. As discussed in Section 1, all these papers but the last assume that the products are ordered according to both underage and overage costs. Under this assumption, a Monge sequence always exists in the full downward substitution problem; Hsu and Bassok (1999) state that without this assumption a Monge sequence does not necessarily exist and develop a solution of complexity of O (N · log(N )), without checking for existence of a Monge sequence. However, one still could use condition (C1) to check for existence of a Monge sequence: if a Monge sequence exists, it only needs to be constructed once; the problem then can be solved for every supply and demand vector in O (N ) steps. For problem structures where no Monge sequence exists, our solution is less efficient than the algorithm of Hsu and Bassok (1999), as we address a more general case. An interesting substitution cost structure, which is more general than full downward substitution but more specific than the general structure, is that of a multi-dimensional ordering, where products differ along several dimensions. In the motivating example from the carbon industry, products differ by quality, length and diameter, and a product can serve as a substitute if and only if it is equal or superior on all three dimensions. Under such a structure, existence of a Monge sequence can be ensured for certain conditions on underage and overage costs. Another interesting problem arises when demand and supply vectors are unity, a special case of the linear sum assignment problem. Our solution for such problems is much easier than for general problems, as there can be at most one violation case on each forbidden arc and each case is corrected with a single path.

500

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

The worst-case computational complexity then reduces to O (N 2 ), being more efficient than the best algorithm for the general assignment problem, which, to the best of our knowledge, is O (N 3 ) (Burkard & Cela, 1999). To conclude, we note that some practical settings might not fit the assumption of additive substitution costs. For example, consider two products i and j, for which both (i, j) and (j, i) are admissible. If bii = b j j = 0 and bi j + b ji > 0, then the costs for the quadruple (i, i), (i, j), (j, i), (j, j) are not additive. For such a setting our solution would need to be adjusted. In the allocation sequences S and Se , arcs of type (i, i) should precede all other arcs. Further, not all violation cases might need to be corrected, and possible cost savings should be checked for each case. Developing an exact solution procedure for the setting with general substitution costs could be an interesting idea for future research. Appendix A. Additional results Lemma 1. For sequence Se constructed in Algorithm 1, Monge sequence condition (MS1) always holds. Proof of Lemma 1. We show that (MS1) holds in Se , i.e., for every (i, j) ∈ E and for every r = i, s = j: if (r, j), (i, s), (r, s) ∈ E and wi j + wrs > wis + wr j , then either (r, j) or (i, s) precedes (i, j) in Se . Assume to the contrary that there exists (i, j) ∈ E and r = i, s = j such that (r, j), (i, s), (r, s) ∈ E and wi j + wrs > wis + wr j , but (i, j) precedes both (r, j) and (i, s) in Se . We check all possible cases of the nodes i, j, r, s: for each case, we show how the rule wi j + wrs > wis + wr j implies a contradiction: •

•

i, j, r, s ∈ {1, . . . , N}: bi j + brs > bis + br j . However, this holds at equality, implying a contradiction. One or more nodes are U or O: – i = U: then (i, j) cannot precede (r, j) as arcs from U node were added last to Se . Similarly, if j = O then (i, j) cannot precede (i, s) as arcs to O node were added last to Se . – r = U: bi j + us > bis + u j ⇒ us − βs > u j − β j . Then (i, s) precedes (i, j) in sequence Se , due to the order of construction in Algorithm 1. – s = O: bi j + or > br j + oi ⇒ or − αr > oi − αi . Then (r, j) precedes (i, j) in sequence Se , due to the order of construction in Algorithm 1. – r = U and s = O: bi j > u j + oi . This is a contradiction to Assumption 1. 

Lemma 2. For sequence Se constructed in Algorithm 1, condition (C1) is equivalent to Monge sequence condition (MS2). Proof of Lemma 2. First we show that if (C1) is true then (MS2) holds, namely for every (i, j) ∈ E and for every r = i, s = j, if (r, j), (i, s) ∈ E and (r, s) ∈ F, either (r, j) or (i, s) precedes (i, j) in Se . Assume to the contrary that there exists (i, j) ∈ E and r = i, s = j such that (r, j), (i, s) ∈ E, (r, s) ∈ F and (i, j) precedes both (r, j) and (i, s). According to the order of adding arcs to Se in Algorithm 1, column j precedes column s, i.e., u j − β j > us − βs , and row i precedes row r, i.e., oi − αi > or − αr , which contradicts (C1). Next we show that if (C1) is not true, (MS2) does not hold. If for some (r, s) ∈ F, (r, j), (i, s), (i, j) ∈ E, oi − αi > or − αr and u j − β j > us − βs , then (i, j) precedes (r, j) and (i, s) in Se . Then (MS2) does not hold in sequence Se for these four arcs.  For two solutions to the allocation problem, X and Y, we create a matrix δ , that has the same structure as the weight matrix W, where for each entry (i, j), δ ij is the allocation difference between the two solutions, i.e., δi j = Yi j − Xi j . We add all arcs (i, j) where δ ij = 0 (i.e., arcs with changed allocation) to a set P. We use PU

and PO , respectively, to denote the subsets of P with underage and overage, and P for the subset of all other arcs in P. Lemma 3. For any two solutions to the allocation problem, X and Y, the total cost change for all arcs (i, j), such that there are no underage/overage arcs in the same column and row in P ((i, O), (U, j) ∈ P), is 0. Proof of Lemma 3. For each arc (i, j) ∈ P, the value of the objective function changes by δi j · bi j = δi j · (αi + β j ) (follows from the additive cost assumption). Due to constraints (2) and (3), the sum of changed quantities on the arcs in each row and each column is 0. Hence, in each considered column j the terms β j cancel out and in each considered row i the terms α i cancel out.  Lemma 4. For any case ((i, j), T, x), after Algorithm 3 it is not possible to add flow from row i to column j without using additional overage or underage arcs. Proof of Lemma 4. In Algorithm 3, we consider all possible paths on the arcs inside the matrix until the maximal flow is found or no more paths are left. Each node (weight matrix entry) is considered until it is saturated or leads to a dead-end. Hence, it is not possible to add more flow, without using additional overage or underage arcs.  Lemma 5. After any case a, that creased, without Algorithm 4 there

correcting case b with Algorithm 4, the flow for was considered before case b, could not be inusing additional O/U arcs. In other words, after are no more correctable violation cases.

Proof of Lemma 5. Assume the following sequence of events: t1 – correcting case a (by Lemma 4, it is corrected maximally); after the correction there is still some overage and/or underage left (depending on the type of the case); t2 – correcting case b; t3 – addressing case a again. Aiming for a contradiction, assume that in t3 we are able to add more flow, i.e., transfer more overage/underage, thus improving case a further. Assume b is the first case after a that allows adding more flow. Then the flow graph in t3 includes at least one node that has been changed in t2 - assume n2 (this arc did not have allocation at t1 but has allocation after t2 ). That is, in t2 there has been some flow path from n1 to n2 (in the same row) and then from n2 to n3 (in the same column). Hence we could use the path n1 − n2 − n3 at t1 , implying the flow there was not maximal, which is a contradiction to Lemma 4.  Lemma 6. At each correction step of Algorithm 4, there exists no arc (i, j) ∈ E such that XiO > 0 and XUj > 0. Proof of Lemma 6. We prove this result by induction. In the solution obtained from Algorithm 2 the statement holds, due to the order of sequence construction in Algorithm 1: arcs of type (i, O) and (U, j) are the last in sequence Se , hence first maximum possible would be allocated on (i, j) and only then on overage or underage. We assume the statement is true for case n and show that it is then is also true for case n + 1. As the statement holds before we start correcting case n + 1, the correction can violate it only if overage or underage is moved to another row/column. When correcting cases of type (UO), we are getting rid of underage and overage simultaneously and do not move them to other rows/ columns. Hence we need to show that when moving overage/underage for cases (O)/(U), we do not create (i, j) ∈ E such that XiO > 0 and XUj > 0. Consider case n + 1 = ((r, s ), U, j ) (the proof for (O) cases is symmetric), where underage is moved to column j, i.e., there exists a path from row r to column s. Assume XiO > 0 for some i. Due to the induction assumption, (i, s) ∈ F. Hence there was a case ((i, s), UO) that has been addressed before the current case (as cases (UO) are addressed before cases (U)) and it has not been possible to create more flow from row i to

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

501

column s; due to Lemma 5 it is not possible now, either. If (i, j) ∈ E, the path (i, j ) − (r, j ), connected to the path from row r to column s, creates a path from i to s, implying a contradiction. 

Lemma 9. In the solution obtained after Algorithms 4 and 5, it is not possible to add flow to a case using other O/U arcs, without affecting more beneficial case/s.

Lemma 7. Lemmas 5 and 6 hold after Algorithm 5.

Proof of Lemma 9. Aiming for contradiction, assume that there exists a case associated with arc (i, j) which we can further improve by changing other O/U arcs, without affecting more beneficial cases. That is, more flow can be created by adding underage/overage arcs to the network of this case. We consider case ((i, j), UO) or ((i, j), U, k); the proof for overage cases is symmetric. Due to constraints (2) and (3), if we change underage/overage arcs, the allocation on at least one O/U arc will increase and on at least one O/U arc decrease. Assume we find a flow path where overage on some row r decreases – this means that XrO > 0 and XUj > 0. By Lemmas 5–7, (r, j) ∈ F and one cannot correct case ((r, j), UO) more. However now we find a path from row r to column j, which is a contradiction. If we find a flow path where in addition underage on some column s decreases, the argument is similar for case ((r, s), UO). If we correct case ((i, j), U, k) and find a flow path where only underage on some column s decreases (no overage change), this means that after Algorithms 4 and 5, Xik > 0 and XUs > 0. If s < k, then (i, s) ∈ F due to Lemma 8 and we cannot correct case ((i, s), U, k) more. However now we find a path from row i to column s, which is a contradiction. If s > k: there exists column l where the underage allocation increases and thus the allocation on some arc (p, l) decreases. If l > j, then (p, j) ∈ F due to Lemma 8 and a case ((p, j), U, l) cannot be corrected more. However now we find a path from row p to column j, which is a contradiction. If l < j: there exists row t such that Xts increases. Clearly (t, l) ∈ F (otherwise the change should be undone as in Algorithm 5); then a case ((t, l), U, s) is created, which is more beneficial than ((i, j), U, k). 

Proof of Lemma 7. In Algorithm 5, we are moving allocation inside the same row when correcting underage or inside the same column when correcting overage. Lemma 5 holds because, if we can correct some case after moving allocations at the same row/column, we could have corrected the case before the change. In Lemma 6, we show that after Algorithm 4 there exists no arc (i, j) ∈ E, such that XiO > 0 and XUj > 0. Assume that after Algorithm 5 such arc exists. Assume in Algorithm 5 underage is moved from column s to column j versus some row r (the proof for overage is symmetric), and assume there exists a row i: XiO > 0 such that (i, j) ∈ E. Due to Lemma 6, (i, s) ∈ F and due to Lemma 5 the case ((i, s), UO) could not be corrected more. However, there exists an additional correction path (i, O ) − (i, j ) − (r, j ) − (r, s ) − (U, s ), which is a contradiction. If XiO became positive after another change in Algorithm 5, i.e., overage was moved from row k to row i, the argument is similar for case ((k, s), UO).  Lemma 8. Consider the solution obtained after Algorithm 5. If there exist some i and j such that: • • •

XUj > 0 and Xik > 0, k > j, or XiO > 0 and Xlj > 0, l > i, or XUj > 0 and XiO > 0,

then (i, j) ∈ F and the case ((i, j), U, k) / ((i, j), O, l) / ((i, j), UO) is not correctable. Proof of Lemma 8. Aiming for contradiction, assume that after Algorithm 5 there exist i and j such that XUj > 0 and Xik > 0, k > j (the proof for overage is symmetric), and the arc (i, j) is either (I) admissible or (II) forbidden when the case ((i, j), U, k) could be corrected. Occurrences of (I) are considered during Algorithm 5. Due to Lemma 5, after Algorithm 4 there are no more correctable violation cases. Then either (a) underage was moved during Algorithm 5 to column j from some column s < j, or (b) allocation was moved during Algorithm 5 to (i, k): (b1) from some arc (i, l), l > k when correcting underage or (b2) from some arc (t, k), t > i when correcting overage; or (c) both (a) and (b) (denoted with (c1) and (c2), respectively). If (a), assume underage was moved versus allocation change in row r. If (i, s) ∈ E, underage should have been moved in Algorithm 5 from s to k before the move from s to j, as k > j. If (i, s) ∈ F: if (I) ((i, j) ∈ E) then the case ((i, s), U, k) could have been corrected on path (i, k ) − (i, j ) − (r, j ) − (r, s ) − (U, s ), and if (II) ((i, j) ∈ F) – on a path from row i to column j connected with the path (r, j ) − (r, s ); this is a contradiction to Lemma 5. If (b1), the argument is similar for arc (i, j) and case ((i, j), U, l), and if (c1) for arc (i, s) and case ((i, s), U, l). If (b2), (i, j) ∈ F due to Lemma 6, and as the case ((i, j), UO) could not be corrected then the case cannot be corrected now, either. If (c2), similarly for arc (i, s) and case ((i, s), UO). Regarding the (UO) cases, due to Lemmas 6 and 7, (i, j) ∈ F. When (UO) cases are corrected, underage and overage are reallocated to the entries inside the matrix. Hence, a new (UO) case could be created if underage (or overage, or both) is moved from some column s, s < j, when correcting some case ((r, s), U, j) (or correcting underage in Algorithm 5). Due to Lemma 6, there is a case ((i, s), UO) that has been addressed before ((r, s), U, j) and is no more correctable. Hence, ((i, j), UO) is also not correctable, as otherwise we could have connected the path from row i to column j with the path from (r, j) to column s (used to correct ((r, s), U, j)), this way correcting ((i, s), UO). 

Lemma 10. For each arc (i, j) ∈ F, the cost improvement from correcting a case of type (UO) associated with this arc is higher than the cost improvement from correcting cases of types (U) or (O) associated with this arc. Proof of Lemma 10. Assume there exist ated with arc (i, j): ((i, j), UO), ((i, j), U, The cost improvement from correcting the than the cost improvement from correcting Assumption 1 for arc (i, s):

three cases associs) and ((i, j), O, r). (UO) case is higher the (U) case, due to

(i, j ),UO = u j − β j + oi − αi > u j − β j − (us − βs ) = (i, j ),U,s . (A.1) Similarly, the cost improvement from correcting the (UO) case is higher than the cost improvement from correcting the (O) case, due to Assumption 1 for arc (r, j):

(i, j ),UO = oi − αi + u j − β j > oi − αi − (or − αr ) = (i, j ),O,r . (A.2)  Lemma 11. For every two cases ((i, j), UO) and ((r, s), UO): if (i, j) precedes (r, s) in sequence Sf then either (i,j),UO ≥ (r,s),UO or the correction order of the two cases does not matter. Proof of Lemma 11. If r = i, i.e., the two forbidden arcs are in the same row, the following holds because j < s:

(i, j ),UO = u j + oi − αi − β j ≥ us + oi − βs − αi = (i,s),UO

(A.3)

If j = s, i.e. the two forbidden arcs are in the same column, the following holds because i < r:

(i, j ),UO = u j + oi − αi − β j ≥ u j + or − β j − αr = (r, j ),UO

(A.4)

The two equations above imply that if i < r and j < s, then

(i,j),UO ≥ (r,s),UO . However, if (a) j > s, i < r or (b) j < s, i > r, it might

502

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

be unclear in which case the cost improvement is higher. Next, we show that it does not matter which case to address first. We note that due to Lemma 6, (i, s) ∈ F and (r, j) ∈ F, hence there exist cases ((i, s), UO) and ((r, j), UO). If (a) then ((i, s), UO) would be the first case to address. After its correction either XUs = 0 and the case ((r, s), UO) does not exist anymore; and/or XiO = 0 and the case ((i, j), UO) does not exist anymore; or XUs > 0 and XiO > 0 and there does not exist any path from row i to column s. Then there does not exist a common path for the cases ((i, j), UO) and ((r, s), UO), because if there was any, the case ((i, s), UO) could have been further corrected. Then the correction paths of these two cases are independent and it does not matter which case to address first. If (b), then ((r, j), UO) would be the first case to address and the argument is similar.  Lemma 12. For every two cases ((i, j), U, l) and ((r, s), U, k): if ((i, j), U, l) precedes ((r, s), U, k) in Algorithm 4, then (i,j),U,l ≥ (r,s),U,k or the correction of ((i, j), U, l) does not impede the correction of ((r, s), U, k). Similarly, for every two cases ((i, j), O, l) and ((r, s), O, k): if ((i, j), O, l) precedes ((r, s), O, k) in Algorithm 4, then (i,j),O,l ≥ (r,s),O,k or the correction of ((i, j), O, l) does not impede the correction of ((r, s), O, k). Proof of Lemma 12. The proof is for underage cases; the argument for overage cases is similar. As the columns are ordered according to the underage costs, it holds: (1) If j = s: ((i, j), U, l) precedes ((r, j), U, k) if l > k. Then (i,j),U,l ≥ (r,j),U,k . (2) If l = k: ((i, j), U, l) precedes ((r, s), U, l) if j < s. Then (i,j),U,l ≥ (r,s),U,l . The argument is similar when j < s and k < l. For the cases ((i, j), U, l) and ((r, s), U, k) where j < s and l < k (or j > s and l > k), it is not immediately clear which case is more profitable. We prove that correction of ((i, j), U, l) does not impede the correction of ((r, s), U, k). Assume ((i, j), U, l) is addressed at time t1 . We show that if ((r, s), U, k) is correctable before t1 , then it is correctable at t2 , after ((i, j), U, l) has been corrected. Aiming for contradiction, assume that ((r, s), U, k) is not correctable at t2 . This means that allocation on some node n1 that could have been transferred to column s has been moved to column j at t1 . Consider (r, j): if (r, j) ∈ E, a path from row r to column s can be created by connecting (r, j) to the path from n1 to column j, in the opposite direction, and then from n1 to column s. If (r, j) ∈ F, then a case ((r, j), U, k) should have been addressed before t1 (as l < k). As after its correction both XUj > 0 and Xrk > 0, the case is not correctable. However it could have been corrected further by taking the path from row r to n1 connected to the path from n1 to column j – a contradiction.  Lemma 13. For every two cases ((i, j), UO) and ((r, s), U, k), either

(i,j),UO ≥ (r,s),U,k or the order of the two cases does not matter. Similarly, for every two cases ((i, j), UO) and ((r, s), O, k), either (i,j),UO ≥ (r,s),O,k or the order of the two cases does not matter. Proof of Lemma 13. Consider two cases ((i, j), UO) and ((r, s), U, k) (the proof for overage cases is similar). For j < s, we see in (A.5) that the cost improvement from correcting the (UO) case is higher than the cost improvement of correcting the (U) case; the first inequality is due to j < s and the second due to Assumption 1 for arc (i, k).

(i, j ),UO = u j − β j + oi − αi ≥ us − βs + oi − αi > us − βs − (uk − βk ) = (r,s ),U,k (A.5) If j > s then (i, s) ∈ F due to Lemma 6, hence there exists a case ((i, s), UO) which is addressed before ((i, j), UO). After its correction either XUs = 0 and the case ((r, s), U, k) does not exist anymore; and/or XiO = 0 and the case ((i, j), UO) does not exist anymore; or XUs > 0 and XiO > 0 and there does not exist any path from row

i to column s. Then there does not exist a common path for the cases ((i, j), UO) and ((r, s), U, k), because if there was any, the case ((i, s), UO) could have been further corrected. Then the correction paths of these two cases are independent and it does not matter which case to address first.  Lemma 14. After addressing all (UO) cases, correction of (O) and (U) cases is independent. Proof of Lemma 14. Consider two cases ((i, j), U, l) and ((r, s), O, k) and assume, aiming for a contradiction, that there exists a common path to correct both of them. Due to Lemma 6, (r, j) ∈ F and there exists a case ((r, j), UO) which could not be corrected any more. However, if we can find paths from row r to column s and from row i to column j and these paths are connected, then we can find a path from row r to column j and correct the case ((r, j), UO) further, which is a contradiction to Lemma 5.  Lemma 15. Assume that when correcting a case ((i, j), UO) (similarly, ((i, j), U, x) or ((i, j), O, x)), in Step 2c of Algorithm 3, node (z, y) cannot be connected to column j, that is, (z, j) ∈ F. If both underage and overage are left after correcting this case, the cases on (z, j) cannot be corrected. Proof of Lemma 15. Assume some case ((z, j), U, x) (or ((z, j), O, x)), that is addressed after ((i, j), UO), can be corrected, that is, there exists at least one path from row z to column j. Then this path could have been used to correct the case ((i, j), UO) further –  contradiction. Appendix B. Proofs of the results given in the main part of the paper The technical lemmas referred to in the following are contained in Appendix A. Proof of Proposition 1. In Lemma 1 we show that (MS1) always holds in Se . In Lemma 2 we show that condition (C1) is equivalent to (MS2) for sequence Se . Following the definition in Shamir and Dietrich (1990), Se thus is a Monge sequence if and only if condition (C1) holds.  Proof of Theorem 1. According to Proposition 1, if (C1) holds, then Se is a Monge sequence. We prove by contradiction that this condition is necessary for a Monge sequence existence. Assume that (C1) is violated, that is, there exist four arcs (r, s) ∈ F, (r, j), (i, s), (i, j) ∈ E such that (C1.1) u j − β j > us − βs , (C1.2) oi − αi > or − αr , and there exists a Monge sequence. According to (MS2), either (r, j) or (i, s) precede (i, j) in the Monge sequence. •

•

If (i, s) precedes (i, j), consider nodes i, j, U, s. – If (i, s) precedes both (i, j) and (U, s), then wis + wU j ≤ wi j + wUs , according to (MS1). us − βs ≤ u j − β j contradicts (C1.1). – If (U, s) precedes (i, s), then consider nodes i, O, U, s. us + oi ≤ αi + βs contradicts Assumption 1. If (r, j) precedes (i, j), consider nodes i, j, r, O. – If (r, j) precedes both (i, j) and (r, O), then wr j + wiO ≤ wi j + wrO , according to (MS1). oi − αi ≤ or − αr contradicts (C1.2). – If (r, O) precedes (r, j), then consider nodes r, O, U, j. u j + or ≤ αr + β j contradicts Assumption 1. 

Proof of Theorem 2. Aiming for contradiction, assume that there exists a solution Y that is better than X, i.e., CY < CX . For each arc (i, j), denote δi j = Yi j − Xi j . Add all arcs (i, j) where δ ij = 0 to a set P. We use PU and PO , respectively, to denote the subsets of P with underage and overage, and P for the subset of all other arcs in P. The sum of changes on the arcs in P in each row and each column is 0, since the total on each column and each row is the same in X and Y, due to constraints (2) and (3). If the allocation from i to j

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

in Y is smaller than in X, then the allocations from i to some other demand node s and to j from some inventory node r have to be larger. In each column and in each row there are either zero or at least two changes: at least one with δ ij > 0 and at least one with δ ij < 0. If P does not include arcs of type (i, O) or (U, j), the objective function values are equal in X and Y, due to Lemma 3 – a contradiction. Assume P includes arcs of type (i, O) / (U, j). For every arc (i, O) ∈ P, δ iO > 0 as XiO = 0 for all i. As the sum of changes in the column O is 0, there exists an arc (U, j) in P such that δ Uj > 0; then there exists an arc (i, j) where in Y there is overage and underage simultaneously while in X there is no overage or underage. Similarly, for every arc (U, j), δ Uj > 0, and there exists a corresponding arc (i, O) that changed. We analyze the change in the objective function between X and Y, that is, CY − CX . The total change for all arcs (i, j) such that there are no underage/overage arcs in the same column and row in P ((i, O), (U, j) ∈ P) is 0, due to Lemma 3. The total change for all other arcs equals   which i:(i,O )∈PO j:( j,U )∈PU min (|δiO |, |δU j | ) (u j + oi − (αi + β j )), is positive, due to Assumption 1. This implies CY > CX , and thus a contradiction.  Proof of Proposition 2. Jointly Lemmas 10–14 imply the proof of the proposition. According to Lemma 10, for each forbidden arc it is more beneficial to address cases of type (UO) associated with this arc before cases of types (U) and (O) associated with this arc. Lemma 11 addresses the order of correction for (UO) type cases and Lemma 12 the order of correction for (U) and (O) type cases. Lemma 13 shows that when comparing any (UO) case with any (U) or (O) case, either the cost improvement from correcting the (UO) case is higher or the order of the two cases does not matter. Finally, Lemma 14 shows that after addressing all (UO) cases, the correction of (O) and (U) cases is independent, and it does not matter which type to address first.  Proof of Theorem 3. Assume that there exists a solution Y that is better than X, i.e., CY < CX . The changes in the allocation between X and Y are captured in matrix δ , where δi j = Yi j − Xi j . We consider all arcs (i, j) where δ ij = 0, that is, where the allocation differs between X and Y. Such arcs are added to a set P which consists of three sets: PU (underage arcs), PO (overage arcs) and P (all other arcs). The total change for all arcs (i, j), such that there are no underage/overage arcs in the same column and row in P ((i, O), (U, j) ∈ P), is 0, due to Lemma 3. Below we show that the total change for all other arcs is non-negative, i.e., CY ≥ CX , thus reaching a contradiction. If PU and PO are not empty, due to constraints (2) and (3), for every arc (i, O) or (U, j) in P there exists at least one arc (r, O) or (U, s) in P with δ of the opposite sign. The cost change from X to Y can be presented similar to the three violation types identified previously: (U) For two arcs (U, j) and (U, s), such that δ Uj < 0 and δ Us > 0: min(|δU j |, |δUs | ) · (us − βs − (u j − β j )) (O) For two arcs (i, O) and (r, O), such that δ iO < 0 and δ rO > 0: min(|δiO |, |δrO | ) · (or − αr − (oi − αi )). (UO) For two arcs (U, j) and (i, O), such that sign(δU j ) = sign(δiO ): sign(δU j ) · min(|δU j |, |δiO | ) · (u j + oi − (αi + β j )) When the cost change for some two arcs is positive, Y is less beneficial than X for this case. Next we show that if for some case the cost change is negative then it is necessarily positive for some other case and the total sum for these cases is non-negative, which leads to contradiction. Assume Y is better than X for some case b. Due to Lemma 8, after Algorithms 4 and 5 there are no correctable violation cases. By Lemma 4, the case has been corrected by Algorithm 3 in a

503

maximal way and by Lemmas 5 and 7, it could not further be improved when correcting subsequent cases. By Lemma 9, case b also cannot be improved by altering other O/U arcs, unless if Y made worse another, more beneficial case, i.e., flow was added to case b by reducing flow for case a < b, which is not profitable due to Proposition 2. By Lemmas 6 and 7, after Algorithms 4 and 5 there does not exist any arc (i, j) ∈ E such that XiO > 0 and XUj > 0. If in Y there are arcs like this, this is not beneficial due to Assumption 1. Due to Lemma 8, there does not exist any arc (i, j) ∈ E such that XUj > 0 and this underage could be swapped beneficially with another column, or XiO > 0 or and this overage could be swapped beneficially with another row.  Proof of Proposition 3. Algorithm 1: Reordering the weight matrix in Step 1 (sorting according to underage and overage costs) requires O (N log N ) operations for rows and for columns. In Step 2 each arc is considered once; hence, the complexity of the algorithm is O (N 2 ). Algorithm 6: In the computationally worst case scenario, no violation of (C1) can be found for any checked quadruple. Then, each row/column is considered at most O (N ) times, as if in some iteration (i, s) and/or (r, j) are forbidden, certain other quadruples are not required to be checked (Steps 2 and 3). Hence the complexity of the algorithm is O (N 2 ). Algorithm 7: The division into sets (Step 1) can be done during the sorting procedures in Step 1 of Algorithm 1, without affecting the complexity. In the worst case, all admissible arcs have the same underage/overage costs with some other arc/s, and all columns/rows need to be addressed in Step 2. When ordering the columns inside each set in Step 2(b)i, we consider each arc once. Sorting the rows in Step 2(b)ii is done in O (N · log(N )) time. When constructing the sequence in Step 3, each admissible arc is considered once. Hence, the worst-case complexity of the algorithm is O (N 2 ).  Proof of Proposition 4. Algorithm 2: After allocating on an arc, either a demand or supply (or both) is satisfied, hence this column or row (or both) would not be considered again. Hence, the algorithm complexity is O (N ). If a Monge sequence exists or if no underage/overage is left, this is the complexity of finding an optimal solution. Algorithm 3: The pre-processing step of ordering columns and rows by entries with/without allocation can be done in O (N · log(N )) steps. There are two types of computationally worst cases: In the first case, we might need to go over all admissible entries to construct a single path. At each iteration, entries that lead to a dead end are deleted from the graph not to be considered again (Step 2d). Hence, if finding a path takes O (e ) steps, there would be no more paths, and the complexity of the algorithm is O (N 2 ). In the second case, multiple paths are considered. We show that the maximum number of addressing each row/column is O (N ). Consider node (z, y) with even depth. If it is connected to all entries in column y (i.e., all entries in column y have positive allocations), then, at each other column, there exists at most one positive allocation for these rows, as there are no quadruples where all four arcs have positive allocations (due to Algorithm 8). The maximal length of each path in Algorithm 3 is O (N ), as there are two entries in the same row and two entries in the same column. After the path construction, we go over it again in Steps 3–6 to find bottlenecks, delete the second direction and save partial paths; no complexity is added. Hence, the total worst case complexity of Algorithm 3 is O (N 2 ). Algorithm 8: After modifying the allocation on a quadruple, at least one entry has zero allocation. Then, this entry should not be considered in other quadruples, and the other quadruples in the same row/column should not be checked. Hence, the number of addressed quadruples is O (N 2 ), and this is the worst-case

504

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

complexity of the algorithm. We note that all quadruples need to be checked only if Algorithm 3 affected all admissible arcs; in other words, the number of steps in Algorithm 8 is of the same order of magnitude as that in Algorithm 3. Algorithm 4: We next evaluate the total complexity of obtaining an optimal solution. In the worst scenario, each forbidden arc has one or more violation cases associated with it. As discussed earlier, after correcting a case, we might not need to address some of the subsequent cases. For example, if after correcting a case ((i, j), UO) both underage and overage are left (in the worst scenario), other cases on (i, j) are eliminated. This way, maximal number of cases to address is O (N 2 ). In the following, we show that if each case takes the maximal number of steps to correct (O (N 2 )), then at most O (N ) cases should be addressed. Alternatively, if all O (N 2 ) cases are addressed, then number of steps for each case would not exceed O (N ). In both options the total complexity is O (N 3 ). If after correcting a case ((i, j), UO) both underage and overage are left, other cases on (r, j), r = i and on (i, s), s = j still might need to be corrected. We show that not all those cases should be addressed. Consider Step 2c of Algorithm 3 if a node cannot be connected to column j, i.e., (z, j) ∈ F. If underage and overage are left after correcting ((i, j), UO), the cases on (z, j) also cannot be corrected (see Lemma 15 in Appendix A), so all cases on arc (z, j) can be eliminated. Assuming the longest running time of Algorithm 3, we would need to check only one case per column, which leads to the total complexity of O (N 3 ). If all cases are checked, then each time in Step 2c of Algorithm 3, (z, j) ∈ E. Then number of steps for correcting a case would not exceed O (N ). A similar logic applies to violation cases of types (U) and (O). If after correcting a case ((i, j), U, s) (or ((i, j), UO)), Xis = 0 (XrO = 0) and XUj > 0, then all other (U) cases on (i, j) still should be addressed. However, if the ((i, j), U, s) case correction requires O (N 2 ) steps, this means that the flow on all paths (but the last) is maximum possible, hence further cases on this arc do not need to consider these paths again. Hence again the complexity is O (N 3 ). The argument is symmetrical for cases of type (O). To finish the proof, we show that a search for the next (U)/(O) case in Step 2a (the (UO) cases are found in advance) does not increase the complexity. Consider the underage cases (the argument for the overage cases is similar): as the underage is moved only to the right, we do not need to search each time for the smallest-indexed column with underage – it can be remembered from the previous case. The search for the arc with positive allocation on the largest-indexed column takes O (log(N )) steps. Thus the complexity does not increase. Algorithm 5: In each step, after modifying the allocation on (i, j), either underage/overage entry or (i, s)/(r, j) entry has zero allocation and should not be checked in further steps. Hence the algorithm worst-case complexity is O (N 2 ).  Appendix C. Algorithms 6–8

Algorithm 6 Checking condition (C1). For every arc (r, s ) ∈ S f , r > 1, s > 1, starting from the last arc in S f : For every arc (i, j ) ∈ Se , i < r, j < s, starting from the first arc in Se : 1. If (i, s ), (r, j ) ∈ E, u j − β j > us − βs and oi − αi > or − αr – stop: (C1) is violated and no Monge sequence exists. 2. If (i, s ) ∈ F , do not check (i, j ) versus all (k, s ) ∈ F , i < k < r, and all (i, l ) ∈ E, j < l < s, versus (r, s ). 3. If (r, j ) ∈ F , do not check (i, j ) versus all (r, k ) ∈ F , j < k < s, and all (l, j ) ∈ E, i < l < r, versus (r, s ). 4. Go to the next quadruple.

Algorithm 7 Constructing sequence Se for products with equal underage/overage cost. 1. Divide all admissible arcs E  into sets with the same underage and/or overage costs. 2. For each set A, |A| > 1: (a) Denote the columns in set A by CA and the rows by RA . (b) Arrange the arcs inside each set: i. For every row i ∈ RA , sort the arcs (i, j ) : j ∈ CA , first admissible arcs and then forbidden arcs; disregarding rows with no admissible arcs. Def note the column of first arc in row i with ji (A ) and the column of last admissible arc in row i with ji∗ (A ). ii. If |RA | > 1, sort the rows i ∈ RA according to ji∗ (A ), in non-ascending order. Denote the last row by i∗ (A ). 3. Construct sequence Se : (a) Choose the next set A with maximal oi and u j . While A is not empty: i. Insert (i, j ), where i = i∗ (A ) and j = ji∗ (A ), into Se . f ii. A = A \ (i, j ). If j = ji (A ), i∗ (A ) = i∗ (A ) − 1, oth∗ ∗ erwise ji (A ) = ji (A ) − 1. (b) Add the underage and the overage arcs.

Remark. As the arcs are inserted according to their underage and overage costs (as in Algorithm 1) and as the order of the arcs inside each set ensures no violation to (MS2) for the arcs with equal underage/overage costs, the algorithm constructs a sequence that is equivalent to that constructed by Algorithm 1 in the sense that Proposition 1, Theorem 1 and all related lemmas continue to hold. Algorithm 8 Quadruple correction (addition to Algorithm 3). For any four arcs (i, j ), (r, j ), (i, s ), (r, s ) ∈ E (r, i, j, s ∈ {1, . . . , N} ) with Xi j > 0, Xis > 0, Xr j > 0 and Xrs > 0, change the allocation to reduce at least one entry to zero: 1. δ = min(Xi j , Xr j , Xis , Xrs ) 2. If Xi j = δ or Xrs = δ : (a) Xi j = Xi j − δ (b) Xrs = Xrs − δ (c) Xr j = Xr j + δ (d) Xis = Xis + δ 3. Otherwise: (a) Xi j = Xi j + δ (b) Xrs = Xrs + δ (c) Xr j = Xr j − δ (d) Xis = Xis − δ

Appendix D. Details on the numerical study presented in Section 5.2 In this section, we present additional details on the results of the numerical study we have conducted and summarized in Section 5.2. Table 1 provides an overview of the parameter values that we varied in a full factorial fashion to create our 4480 scenarios (drawing 10 replications for each of the 448 configurations). In Tables 2-5, the columns “CS2 algorithm” and “Our algorithm” provide information on the running times of the CS2 algorithm and of our algorithm, respectively, both in seconds. The last column then provides information on the improvement rate, which is

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506 Table 1 Parameter values used for the numeric study. Parameter

Description

Values

N oi

Number of products Overage costs

{10 0 0, 20 0 0, . . . , 70 0 0}

ui

Underage costs

α, β

Components of additive substitution cost (Any substitution is feasible with 50% probability) Supply: mean (Normally distributed, truncated) Supply: coefficient of variation (Normally distributed, truncated) Demand: mean (Normally distributed, truncated) Demand: coefficient of variation (Normally distributed, truncated)

μsupply covsupply

μdemand covdemand

L: Uniform [20,40] H: Uniform [50,150] L: Uniform [20,40] H: Uniform [50,150] Uniform [1,20]

Table 5 Average computational times and improvement rates by volume and variability in demand and supply.

μsupply

covsupply

μdemand

covdemand

CS2 Our Improvement algorithm algorithm (%)

H

H

H

H L H L H L H L H L H L H L H L

10.66 10.27 11.75 11.51 10.23 9.76 11.73 10.91 11.38 11.32 10.26 10.00 11.22 9.74 10.07 9.69

L L

L: 200 H: 300 L: 10% H: 50% L: 200 H: 300 L: 10% H: 50%

505

H L

L

H

H L

L

H L

4.67 4.70 4.60 4.66 4.68 4.58 4.63 4.57 4.71 4.60 4.76 4.67 4.72 4.59 4.80 4.55

50.8 50.3 54.7 54.2 49.8 48.0 55.6 53.0 52.6 53.5 49.9 48.5 51.5 48.9 48.3 47.6

Table 2 Statistics on computational times and improvement rates.

Mean Min 10th percentile 20th percentile 30th percentile 40th percentile 50th percentile 60th percentile 70th percentile 80th percentile 90th percentile Max

CS2 algorithm

Our algorithm

Improvement (%)

10.66 0.22 0.34 1.44 3.10 4.35 7.31 11.03 14.82 19.99 26.71 46.72

4.65 0.12 0.20 0.78 1.59 2.03 3.43 5.29 6.58 8.56 11.15 16.85

51.1 4.2 36.8 42.7 46.5 49.5 52.1 54.7 57.4 60.3 64.0 81.3

Table 3 Average computational times and improvement rates by number of products. N

CS2 algorithm

Our algorithm

Improvement (%)

10 0 0 20 0 0 30 0 0 40 0 0 50 0 0 60 0 0 70 0 0

0.33 1.52 3.82 7.39 12.62 19.54 29.36

0.19 0.80 1.88 3.44 5.75 8.53 11.99

42.0 46.7 49.5 52.4 53.5 55.3 58.1

Table 4 Average computational times and improvement rates by overage and underage costs. oi

ui

CS2 algorithm

Our algorithm

Improvement (%)

H

H L H L

10.84 10.53 10.80 10.46

4.66 4.66 4.63 4.67

51.8 51.2 51.0 50.3

L

the percentage reduction in computational time achieved by our algorithm, that is, [(computational time of CS2 algorithm) minus (computational time of our algorithm)]/(computational time of CS2 algorithm). We first provide some summary statistics of computational times and performance improvement in Table 2. Tables 3–5 then contain some insights with regard to the impact of the different parameters on our results. Specifically, we provide average computational times and average improvement rates for different numbers of products (Table 3), for different underage and overage costs (Table 4), and for different supply and demand configurations (Table 5).

Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ejor.2019.02.049. References Ababei, C. (2009). C++ implementation of Goldberg’s CS2 scaling minimum-cost flow algorithm. Online. https://github.com/eigenpi/CS2-CPP. Ahuja, R. K., Magnanti, T. L., & Orlin, J. B. (1989). Chapter IV: Network flows. In Handbooks in operations research and management science (pp. 211–369). Elsevier. Alon, N., Cosares, S., Hochbaum, D. S., & Shamir, R. (1989). An algorithm for the detection and construction of Monge sequences. Linear Algebra and its Applications, 114, 669–680. Bansal, S. (2010). Optimal assortments of vertically differentiated products: analytical solution and properties. Dissertation. Austin, TX: The University of Texas at Austin. Bassok, Y., Anupindi, R., & Akella, R. (1999). Single-period multiproduct inventory models with substitution. Operations Research, 47(4), 632–642. Bitran, G. R., & Leong, T. Y. (1992). Deterministic approximations to co-production problems with service constraints and random yields. Management Science, 38(5), 724–742. Burkard, R. E. (2007). Monge properties, discrete convexity and applications. European Journal of Operational Research, 176(1), 1–14. Burkard, R. E., & Cela, E. (1999). Linear assignment problems and extensions. Handbook of combinatorial optimization (pp. 75–149). US: Springer. Burkard, R. E., Klinz, B., & Rudolf, R. (1996). Perspectives of Monge properties in optimization. Discrete Applied Mathematics, 70(2), 95–161. Chen, J., Yao, D. D., & Zheng, S. (1999). Optimal control of a multi-product inventory system with substitution. In Proceedings of the thirty-eighth IEEE conference on decision and control: 1 (pp. 468–473). Dietrich, B. L. (1990). Monge sequences, antimatroids, and the transportation problem with forbidden arcs. Linear Algebra and its Applications, 139, 133–145. Ford, L. R., Jr., & Fulkerson, D. R. (1956). Maximal flow through a network. Canadian Journal of Mathematics, 8(3), 399–404. Goldberg, A. V. (1997). An efficient implementation of a scaling minimum-cost flow algorithm. Journal of Algorithms, 22, 1–29. Hitchcock, F. L. (1941). The distribution of a product from several sources to numerous localities. Journal of Mathematics and Physics, 20(1), 224–230. Hoffman, A. J. (1963). On simple linear programming problems. Proceedings of Symposia in Pure Mathematics, 7, 317–327. Hsu, A., & Bassok, Y. (1999). Random yield and random demand in a production system with downward substitution. Operations Research, 47(2), 277–290. Itai, A., & Shiloach, Y. (1979). Maximum flow in planar networks. SIAM Journal on Computing, 8(2), 135–150. Kantorovich, L. V. (1942). On the translocation of masses. USSR National Academy of Sciences, 37(7–8), 199–201. Kovács, P. (2015). Minimum-cost flow algorithms: An experimental evaluation. Optimization Methods and Software, 30(1), 94–127. Lang, J. C. (2010). Production and inventory management with substitution. Springer. Orlin, J. B. (1993). A faster strongly polynomial minimum cost flow algorithm. Operations Research, 41(2), 338–350. Queyranne, M., Spieksma, F., & Tardella, F. (1998). A general class of greedily solvable linear programs. Mathematics of Operations Research, 23(4), 892–908. Rao, U. S., Swaminathan, J. M., & Zhang, J. (2004). Multi-product inventory planning with downward substitution, stochastic demand and setup costs. IIE Transactions, 36, 59–71.

506

Y. Tseytlin and H.S. Heese / European Journal of Operational Research 277 (2019) 492–506

Shamir, R. (1993). A fast algorithm for constructing Monge sequences in transportation problems with forbidden arcs. Discrete Mathematics, 114(1–3), 435–444. Shamir, R., & Dietrich, B. L. (1990). Characterization and algorithms for greedily solvable transportation problems. In Proceedings of the first annual ACM-SIAM symposium on discrete algorithms. Philadelphia, PA. Shin, H., Park, S., Lee, E., & Benton, W. C. (2015). A classification of the literature on the planning of substitutable products. European Journal of Operational Research, 246(3), 686–699.

Shumsky, R. A., & Zhang, F. (2009). Dynamic capacity management with substitution. Operations Research, 57(3), 671–684. Weiß, C., Knust, S., Shakhlevich, N. V., & Waldherr, S. (2016). The assignment problem with nearly Monge arrays and incompatible partner indices. Discrete Applied Mathematics, 211, 183–203.