On seeking efficient Pareto optimal points in multi-player minimum cost flow problems with application to transportation systems
NNoname manuscript No. (will be inserted by the editor)
On seeking efficient Pareto optimal points inmulti-player minimum cost flow problems withapplication to transportation systems
Shuvomoy Das Gupta & Lacra PavelAbstract
In this paper, we propose a multi-player extension of the minimum cost flow probleminspired by a transportation problem that arises in modern transportation industry. We associateone player with each arc of a directed network, each trying to minimize its cost function subjectto the network flow constraints. In our model, the cost function can be any general nonlinearfunction, and the flow through each arc is an integer. We present algorithms to compute efficientPareto optimal point(s ), where the maximum possible number of players (but not all) minimizetheir cost functions simultaneously. The computed Pareto optimal points are Nash equilibriumsif the problem is transformed into a finite static game in normal form.
In recent years, product transportation systems are increasingly being dominated by retailerssuch as Amazon, Alibaba, and Walmart, who utilize e-commerce solutions to fulfill customersupply chain expectations [2, 3, 4, 5]. In the supply chain strategy of these retailers, productslocated at different warehouses are shipped to geographically dispersed retail centers by differ-ent transportation organizations. These transportation organizations ( carriers ) compete amongthemselves and transport goods between warehouses and retail centers over multiple transporta-tion links. For example, Amazon uses FedEx, UPS (United Parcel Service), AAR (Associationof American Railroads), and other competing organizations to provide transportation services[6, Chapter 11].Product shipment from a warehouse to a retail center requires contracting multiple competingcarriers, e.g. , a common shipment may comprise of (i) a trucking shipment from the warehouseto a railhead provided by FedEx, then (ii) a rail shipment provided by AAR, and finally (iii)a trucking shipment from the rail yard to the retail center provided by UPS. It is common
A preliminary version of this work without proofs appeared in [1].Shuvomoy Das GuptaResearch & Technology DepartmentThales Canada, Transportation Solutions105 Moatfield Drive, Toronto, Ontario, CanadaE-mail: [email protected] PavelSystems Control GroupDepartment of Electrical & Computer EngineeringUniversity of Toronto10 King’s College Road, Toronto, Ontario, CanadaE-mail: [email protected] a r X i v : . [ m a t h . O C ] M a y hat different competing carriers operate over the same transportation link, e.g. , both FedExand UPS provide trucking shipment services for Amazon over the same transportation link [4,Section 9.1]. The goal of each carrier is to maximize its profit (minimize its cost). So a relevantquestion in this regard is how to determine a good socially-optimal solution for the competingcarriers.We can formulate the multi-carrier transportation setup described above as a multi-playerextension of the well-known minimum cost flow problems. The transportation setup can bemodeled by a directed network, where a warehouse is a supply node and a retail center isa demand node. A transportation link is an arc, and a carrier ( e.g. , FedEx, UPS) in chargeof transporting products over that transportation link is a player. The products transportedover the directed network are modeled as the flow through the network, and customer supplychain expectations can be modeled as mass-balance constraints. Each of the carriers is trying tomaximize its profit by maximizing the total number of products that it transports. Carriers arecompeting for a limited resource, namely the total number of products to be transported. Notethat one carrier making the maximum profit may impact the other carriers in a negative mannerand even violate customer supply chain expectations. Our goal is to define a socially-optimalsolution concept in this setup and to provide algorithms to calculate such a solution.The problem above can be generalized as a multi-player minimum cost flow problem [1]. Weassociate one player with one arc of the directed and connected network graph. Parallel arcsbetween two nodes denote two competing carriers over the same transportation link. Each ofthe players is trying to minimize its cost function, subject to the network flow constraints. Theflow through each arc is taken to be an integer. This assumption incurs no significant loss ofgenerality because by suitable scaling we can use the integer model to obtain a real-valued flowvector to any desired degree of accuracy. Naturally, defining an efficient solution concept is ofinterest in such a problem setup.The unambiguously best choice would be a utopian vector optimal solution, which minimizesall the objectives simultaneously. However, this is unlikely to exist in practice [7, page 176]. A generic Pareto optimal point, where none of the objective functions can be improved withoutworsening some of the other objective values, is a better solution concept. However, there canbe numerous such generic Pareto optimal points, many of which would be poor in quality orefficiency [7, Section 4.7.5]. In this paper, we investigate an efficient Pareto optimal point [8,Section 2.2] as a good compromise solution that finds a balance between the utopian vectoroptimality and the generic Pareto optimality. We present algorithms to compute an efficientPareto optimal point, i.e. , a Pareto optimal point where the maximum possible (but not all)number of players minimize their cost functions simultaneously. Related work.
The classic version of the minimum cost flow problem has a linear cost functionfor which polynomial time algorithms exist [9, Chapter 10], even when the network structure( e.g. , nodes and arcs) is subject to uncertainty [10]. The polynomial runtime can be improvedto a linear runtime when the minimum cost flow problem has a special structure [11]. However,for nonlinear cost functions, results exist for very specific cases, and no work seems to existfor arbitrary nonlinear functions. The most commonly used nonlinear cost function is the fixedcharge function, where the cost on the arc is if the flow is zero and affine otherwise. Networkflow problems with fixed charge costs are studied in [12, 13, 14], where the integrality conditionon the flow is not considered. Minimum cost flow problems with concave cost functions arestudied in [15, 16, 17, 18, 19]. Minimum cost flow problems with piece-wise linear cost functionsis investigated in [20, 21, 22]. A dynamic domain contraction algorithm for nonconvex piece-8wise linear network flow problems is proposed in [23]. A particle swarm optimization basedhybrid algorithm to solve the minimum concave cost network flow problem is investigated in224]. Finding Pareto optimal points in multi-objective network flow problem with integer flowshas been limited so far to linear cost functions and two objectives [25, 26, 27, 28]. A multi-player minimum cost flow problem with nonconvex cost functions is explored in [1]. Integermulti-commodity flow problems are investigated in [29, 30, 31]; the underlying problems areoptimization problems in these papers. Contributions.
In this paper, we propose an extension of the minimum cost flow problem to amulti-player setup and construct algorithms to compute efficient Pareto optimal solutions. Ourproblem can be interpreted as a multi-objective optimization problem [7, Section 4.7.5] with theobjective vector consisting of a number of univariate general nonlinear cost functions subjectto the network flow constraints and integer flows. In comparison with existing literature, we donot require the cost function to be of any specific structure. The only assumption on the costfunctions is that they are proper. In contrast to relevant works in multi-objective network flowproblems, our objective vector has an arbitrary number of components; however, each componentis a function of a decoupled single variable. We extend this setup to a problem class that is strictly larger than, but that contains the network flow problems. We develop our algorithms for thislarger class.We show that, although in its original form the problem has coupled constraints bindingevery player, there exists an equivalent variable transformation that decouples the optimizationproblems for a maximal number of players. Solving these decoupled optimization problems canpotentially lead to a significant reduction in the number of candidate points to be searched.We use the solutions of these decoupled optimization problems to reduce the size of the setof candidate efficient Pareto optimal solutions even further using algebraic geometry. Then wepresent algorithms to compute efficient Pareto optimal points that depend on a certain semi-algebraic set being nonempty. We also present a penalty based approach applicable when thatsemi-algebraic set is empty; such an approach can be of value to network administrators andpolicy makers. To the best of our knowledge, our methodology is novel. The computed efficientPareto optimal point has some desirable properties: (i) it is a good compromise solution betweenthe utopian vector optimality and the generic Pareto optimality, and (ii) it is a Nash equilibriumif we convert our setup into a finite static game in normal form.The rest of the paper is organized as follows. Section 1 presents notation and notions usedin the paper. Section 2 describes the problem for directed networks and the extension to astrictly larger class. In Section 3, we transform the problem under consideration into decoupledoptimization problems for a number of players and reformulate the optimization problems forthe rest of the players using consensus constraints. Section 4 presents algorithms for computingefficient Pareto optimal points for our problem if a certain semi-algebraic set is nonempty. Section5 discusses a penalty based approach if the semi-algebraic set is empty. Section 6 presents anillustrative numerical example of our methodology in a transportation setup. Finally, in Section7 we present some concluding remarks regarding our methodology.
Proofs are provided inthe Appendix.
Notation and notions.
We denote the sets of real numbers, integers, and natural numbersby R , Z , and N , respectively. The i th column, j th row, and ( i, j ) th component of a matrix A ∈ R m × n is denoted by A i , a Tj , and a ij , respectively. The submatrix of a matrix A ∈ R m × n ,which constitutes of its rows r , r + 1 , . . . , r and columns c , c + 1 , . . . , c , is denoted by A [ r : r ,c : c ] ∈ R ( r − r +1) × ( c − c +1) . If we make two copies of a vector x ∈ R n , then the copies aredenoted by x (1) and x (2) . By I i,j ∈ R n × n , we denote a matrix that has a on its ( i, j ) th positionand everywhere else. If C and D are two nonempty sets, then C + D = { x + y | x ∈ C, y ∈ D } .If A ∈ R m × n is a matrix and C is a nonempty set containing n -dimensional points, then3 C = { Ax | x ∈ C } . The set of consecutive integers from 1 to n is denoted by [ n ] = { , , . . . , n } and m to n is denoted by [ m : n ] = { m, m + 1 , . . . , n } . If we have two vectors x, y ∈ R n , then x (cid:23) y means ( ∀ i ∈ [ n ]) x i ≥ y i , and we write x − y ∈ R n + . Depending on the context, can represent the scalar zero, a columnvector of zeros, or a matrix with all the entries zero, e.g. , if we say x ∈ R n and x = 0 , then represents a column vector of n zeros. On the other hand, if we say A ∈ R m × n and A = 0 , then represents a matrix of m rows and n columns with all entries zero.Consider a standard form polyhedron { x ∈ R n | Ax = b, x (cid:23) } , where A ∈ R m × n is a fullrow rank matrix, and b ∈ R m . A basis matrix of this polyhedron is constructed as follows. Wepick m linearly independent columns of A and construct the m × m invertible square submatrixout of those columns; the resultant matrix is called a basis matrix. The concept of basis matrixis pivotal in simplex algorithm, which is used to solve linear programming problem. Suppose wehave a polyhedron defined by linear equality and inequality constraints in R n . Then ˜ x ∈ R n iscalled a basic solution of the polyhedron, if all the equality constraints are active at ˜ x , and outof all the active constraints (both equality and inequality) that are active at ˜ x , there are n ofthem that are linearly independent. Let G = ( M , A ) be a directed connected graph associated with a network, where M = [ m + 1] isthe set of nodes, and A = [ n ] is the set of (directed) arcs. With each arc j ∈ A , we associate oneplayer, which we call the j th player. The variable controlled by the j th player is the nonnegativeinteger flow on arc j , denoted by x j ∈ Z . Each player is trying to minimize a proper cost function f j : Z → R , subject to the network flow constraints. We assume each of the cost functions isproper, i.e. , for all i ∈ [ n ] we have −∞ (cid:54)(cid:54)∈ f i ( Z ) , and dom f i = { x i ∈ Z n | f i ( x i ) < + ∞} (cid:54) = ∅ . There is an upper bound u j , which limits how much flow the j th player can carry through arc j .Without any loss of generality, we take the lower bound on every arc to be [9, Page 39]. Thesupply or demand of flow at each node i ∈ M is denoted by b i . If b i > , then i is a supply node;if b i < , then i is a demand node with a demand of − b i , and if b i = 0 , then i is a trans-shipmentnode. We allow parallel arcs to exist between two nodes.In any minimum cost flow problem there are three types of constraints. (i) Mass balance constraint. The mass balance constraint states that for any node, the outflowminus inflow must equal the supply/demand of the node. We describe the constraint using the node-arc incidence matrix . Let us fix a particular ordering of the arcs, and let x ∈ Z n be theresultant vector of flows. First, we define the augmented node-arc incidence matrix ˜ A , whereeach row corresponds to a node, and each column corresponds to an arc. The symbol ˜ a ij denotesthe ( i, j ) th entry of ˜ A that corresponds to the i th node and the j th arc; ˜ a ij is if i is the startnode of the j th arc, − if i is the end node of the j th arc, and otherwise. Note that parallelarcs will correspond to different columns with same entries in the augmented node-arc incidencematrix. So every column of ˜ A has exactly two nonzero entries, one equal to and one equal to − indicating the start node and the end node of the associated arc. Denote, ˜ b = ( b , . . . , b m , b m +1 ) .Then in matrix notation, we write the mass balance constraint as: ˜ Ax = ˜ b. The sum of the rowsof ˜ A is equal to zero vector, so one of the constraints associated with the rows of the linearsystem ˜ Ax − ˜ b is redundant, and by removing the last row of the linear system, we can arriveat a system, Ax = b, where A = ˜ A [1: m, n ] is the node-arc incidence matrix , and b = ˜ b [1: m ] . Thevector b is also called the resource vector . Now A is a full row rank matrix under the assumptionof G being connected and (cid:80) i ∈N b i = 0 [32, Corollary 7.1].4 ii) Flow bound constraint. The flow on any arc must satisfy the lower bound and capacityconstraints, i.e. , (cid:22) x (cid:22) u . The flow bound constraint can often be relaxed or omitted inpractice [7, pages 550-551]. In such cases, the flow direction is flexible, and overflow is allowedsubject to a suitable penalty. (iii) Integrality constraint. The flow on any arc is integer-valued, i.e. , x ∈ Z n . This does notincur a significant loss of generality (see Remark 1 below).So the constraint set, which we denote by P can be written as, P = { x ∈ Z n | Ax = b, (cid:22) x (cid:22) u } , (1)and the subset of P containing only the equality constraints is denoted by Q , i.e. , Q = { x ∈ Z n | Ax = b } . (2)Consider a set of players denoted by [ n ] . The decision variable controlled by the i th playeris x i ∈ Z , i.e. , each player has to take an integer-valued action. The vector formed by all thedecision variables is denoted by x = ( x , x , . . . , x n ) ∈ Z n . By x − i ∈ Z n − , we denote the vectorformed by all the players decision variables except i th player’s decision variable. To put emphasison the i th player’s variable we sometimes write x as ( x i , x − i ) . Each player has a cost function f i ( x i ) : Z → R , which depends on its variable x i . The goal of the i th player for i ∈ [ n ] , givenother players’ strategies x − i ∈ Z n − , is to solve the minimization problemminimize x i f i ( x i ) subject to A i x i + n (cid:88) j =1 ,j (cid:54) = i A j x j = b (cid:22) ( x i , x − i ) (cid:22) ux ∈ Z n . (3)Our objective is to calculate efficient Pareto optimal points for the problem. We define vectoroptimal points first, then Pareto optimal point, and finally efficient Pareto optimal point. Definition 1 ( Vector optimal point ) In problem (3), a point x vo ∈ P is vector optimal if itsatisfies the following: ( ∀ ˜ x ∈ P ) ( ∀ i ∈ [ n ]) f i (˜ x i ) ≥ f i ( x vo i ) . Definition 2 (Pareto optimal point)
In problem (3), a point x po ∈ P is Pareto optimal if itsatisfies the following: there does not exist another point ˜ x ∈ P such that ( ∀ i ∈ [ n ]) f i (˜ x i ) ≤ f i ( x po i ) , (4)with at least one index j ∈ [ n ] satisfying f j (˜ x j ) < f j ( x po j ) . Definition 3 (Efficient Pareto optimal point)
In problem (3), a point x ∗ is an efficientPareto optimal solution , if it is Pareto optimal and it achieves partial vector optimality over amaximal subset of [ n ] ; i.e. , x ∗ satisfies (4) and the set S ⊆ [ n ] that satisfies ( ∀ ˜ x ∈ P ) ( ∀ i ∈ S ) f i (˜ x i ) ≥ f i ( x ∗ i ) , is maximal. 5 emark 1 In our model, we have taken the flow through any arc of the network to be an integer.However this assumption does not incur a significant loss of generality because we can use ourinteger model to obtain a real valued Pareto optimal solution to an arbitrary degree of accuracyby using the following scaling technique [9, page 545]. Suppose we want a real valued Paretooptimal solution x ∗ . Such a real valued Pareto optimal solution corresponds to a modified versionof problem (3) with the last constraint being changed to x ∈ R n . In practice, we always havean estimate of how many points after the decimal point we need to consider. So in the modifiedproblem we substitute each x i for i ∈ [ n ] with y i /α , where y i ∈ Z , and α is chosen dependingon the desired degree of accuracy ( e.g. , α = , or , or larger depending on how manypoints after the decimal point we are interested in). Then we proceed with our methodologydescribed in the subsequent sections to compute Pareto optimal solutions over integers. Let y ∗ be one such integer-valued Pareto optimal solution. Then x ∗ = ( x ∗ i ) ni =1 = (cid:0) α y ∗ i (cid:1) ni =1 correspondsto a real-valued Pareto optimal solution to the degree of accuracy of /α . Remark 2
We can formulate our problem as an n person finite static game in normal form [33,pages 88-91]. In problem (3), a player i ∈ [ n ] has a finite but possibly astronomical number ofalternatives to choose from the feasible set. Let m i be the number of feasible alternatives availableto player i . Furthermore, define the index set M i = [ m i ] = { , . . . , m i } with a typical elementof the set designated as n i , which corresponds to some flow x i . If player chooses a strategy n ∈ M , player chooses a strategy n ∈ M , and so on for all the other players, then the costincurred to player i is a single number a i n ··· n n that can be determined from problem (3). Theordered tuple of all these numbers (over i ∈ [ n ] ), i.e. , (cid:0) a n ··· n n , a n ··· n n , . . . , a n n ··· n n (cid:1) , constitutesthe corresponding unique outcome of the game. For a strategy ( n · · · n n ) that violates any of theconstraints in problem (3), the cost is taken as + ∞ . Players make their decisions independently,and each player unilaterally seeks the minimum possible loss, of course by also taking intoaccount the possible rational and feasible choices of the other players. The noncooperative Nashequilibrium solution concept within the context of this n -person game can be described as follows. Definition 4 (Noncooperative Nash equilibrium) [33, page 88] An n -tuple of strategies (cid:0) n Nash1 , . . . , n Nash n (cid:1) with n Nash i ∈ M i for all i ∈ [ n ] , is said to constitute a noncooperative Nashequilibrium solution for the aforementioned n -person nonzero-sum static finite game in normalform if the following n inequalities are satisfied for all n i ∈ M i and all i ∈ [ n ] : a i, Nash = a i n Nash1 n Nash2 ··· n Nash n ≤ a i n Nash1 n Nash2 ··· n i ··· n Nash n . (5)The flow corresponding to (cid:0) n Nash1 , . . . , n Nash n (cid:1) is denoted by x Nash = (cid:0) x Nash1 , . . . , x
Nash n (cid:1) and iscalled the noncooperative Nash equilibrium flow . Here, the n -tuple ( a , Nash , a , Nash , . . . , a n, Nash ) is known as a noncooperative (Nash) equilibrium outcome of the n -person game in normal form.Note that the strategy associated with an efficient Pareto optimal solution x ∗ in Definition 3also satisfies (5) in Definition 4, thus it is a noncooperative Nash equilibrium flow.We now extend the class of problems that we are going to investigate, which is strictly larger thanthe class defined by problem (3) and contains it. We will develop our algorithms for this largerclass. Everything defined in the previous subsection still holds, and we extend the constraint set P and the equality constraint set Q as follows. The structure of P is still that of a standardform integer polyhedron, i.e. , P = { x ∈ Z n | Ax = b, (cid:22) x (cid:22) u } , where A is a full row rankmatrix, but it may not necessarily be a node-arc incidence matrix only. Denote the convex hullof the points in P by conv P . Consider the relaxed polyhedron relaxed P = { x ∈ R n | Ax = b, (cid:22) x (cid:22) u } , where we have relaxed the condition of x being an integer-valued vector. We nowimpose the following assumption. 6 ssumption 1 For any integer-valued b , relaxed P has at least one integer-valued basic solu-tion. As vertices of a polyhedron are also basic solutions [32, page 50], if conv P and relaxed P share at least one vertex, Assumption 1 will be satisfied. We can see immediately that if A is a node-arc incidence matrix, then P will belong to this class as conv P = relaxed P fornetwork flow problems [34, Chapter 19]. In other practical cases of interest, the matrix cansatisfy Assumption 1, e.g. , matrices with upper or lower triangular square submatrices withdiagonal entries 1, sparse matrices with m variables appearing only once with coefficients one, etc .Moreover, at the expense of adding slack variables (thus making a larger dimensional problem),we can turn the problem under consideration into one satisfying Assumption 1, though thecomputational price may be heavy.In the rest of the paper, whenever we mention (1), (2), and (3), they correspond to thislarger class of problems containing the network flow problems, and the full row rank matrix A is associated with this larger class. So the results developed in the subsequent sections will holdfor a network flow setup. Remark 3
Before proceeding any further, we recall that integer programming problems are
N P -hard, and even determining the existence of one feasible point in P is N P hard [35, page 242].So problem (3) is at least as hard.
In this section, we describe how to transform problem (3) into n − m decoupled optimizationproblems for the last n − m players and how to reformulate the optimization problems for the restof the players using consensus constraints. These transformation and reformulation are necessaryfor the development of our algorithms.3.1 Decoupling optimization problems for the last n − m playersFirst, we present the following lemma. Recall that, an integer square matrix is unimodular if itsdeterminant is ± . Lemma 1
Assumption 1 holds if and only if we can extract a unimodular basis matrix from A . Without any loss of generality, we rearrange the columns of the matrix A so that the unimodularbasis matrix constitutes the first m columns, i.e. , if A = [ A | A | . . . | A m | A m +1 | . . . | A n ] ,then det([ A | A | . . . | A m ]) = ± , and we reindex the variables accordingly. Let us denote [ A | A | . . . | A m ] = B , so A = [ B | A m +1 | . . . | A n ] . Next we present the following Lemma. Lemma 2
Let, C = B − A and d = B − b . Then C ∈ Z m × n , d ∈ Z m , and the sets Q and P (defined in (1) and (2) , respectively) have the equivalent representations: Q = { x ∈ Z n | Cx = d } , (6) P = { x ∈ Z n | Cx = d, (cid:22) x (cid:22) u } . (7)Before we present the next result, we recall the following definitions and facts. A full row rankmatrix A ∈ Z m × n is in Hermite normal form , if it is of the structure [ B | , where B ∈ Z m × m isinvertible and lower triangular, and represents { } m × n − m [34, Section 4.1]. Any full row rankinteger matrix can be brought into the Hermite normal form using elementary integer columnoperations in polynomial time in the size of the matrix [35, Page 243].7 emma 3 The matrix C , as defined in Lemma 2, can be brought into the Hermite normal form [ I | by elementary integer column operations, more specifically by adding integer multiple of onecolumn to another column. Lemma 4
There exists a unimodular matrix U such that CU = [ I | . The following theorem is key in transforming the problem into an equivalent form with m decoupled optimization problems for players m + 1 , m + 2 , . . . , n . Theorem 1
The constraint set Q defined in (6) is nonempty. Furthermore, any vector x ∈ Q can be maximally decomposed in terms of a new variable z ∈ Z n − m as follows: x ∈ Q ⇔ ∃ z ∈ Z n − m such that x = d − h T zd − h T z ... d m − h Tm zz ... z n − m , (8) where d i is the i th component of d = B − b , and h Ti ∈ Z × n − m is the i th row of B − A [1: m,m +1: n ] . We can transform our problem using the new variable z . The advantage of this transformation isthat for player m + 1 , m + 2 , . . . , n , we have decoupled optimization problems. By solving thesedecoupled problems, we can reduce the constraint set significantly (especially when the numberof minimizers for the constraint set are small).From Theorem 1, x i = z i − m for i ∈ [ m +1 : n ] . For problem (3), we can write the optimizationproblem for any player m + i for i ∈ [ n − m ] as follows.minimize z i f i ( z i ) subject to ≤ z i ≤ u i z i ∈ Z . (9)Each of these optimization problems is a decoupled univariate optimization problem, whichcan be easily solved graphically. We can optimize over real numbers, find the minimizers ofthe resultant relaxed optimization problem, determine whether the floor or ceiling of such aminimizer results in the minimum cost, and pick that as a minimizer of the original problem.Solving n − m decoupled optimization problems immediately reduces the constraint set intoa much smaller set. Let us denote the set of different optimal solutions for player m + i for i ∈ [ n − m ] as D i = { z i, , z i, , . . . , z i,p i } sorted from smaller to larger, where p i is the totalnumber of minimizers. Define, D = × n − mi =1 D i . Note that D (cid:54) = ∅ .3.2 Consensus reformulation for the first m playersIn this section, we transform the optimization problems for the first m players using (8) andconsensus constraints. Consider the optimization problems for the first m players in variable z ,which have coupled costs due to (8). We deal with the issue by introducing consensus constraints [36, Section 5.2]. We provide each player i ∈ [ m ] with its own local copy of z , denoted by z ( i ) ∈ Z n − m , which acts as its decision variable. This local copy has to satisfy the followingconditions. First , using (8) for any i ∈ [ m ] , x i = d i − h Ti z ( i ) . The copy z ( i ) has to be in consensus8ith the rest of the first m players, i.e. , z ( i ) = z ( j ) for all j ∈ [ m ] \ { i } . Second , the copy z ( i ) hasto satisfy the flow bound constraints, i.e. , ≤ d i − h Ti z ( i ) ≤ u i for all i ∈ [ m ] . Third , for the last n − m players z i ∈ D i , as obtained from the solutions of the decoupled optimization problems(9), so z ( i ) has to be in D , i.e. , for all i ∈ [ m ] we have z ( i ) = z ∈ D ⇔ ( ∀ j ∈ [ n − m ]) z ( i ) j ∈ D j . Combining the aforementioned conditions, for all i ∈ [ m ] , the i th player’s optimization problemin variable z ( i ) can be written as:minimize z ( i ) ¯ f i (cid:16) z ( i ) (cid:17) = f i ( d i − h Ti z ( i ) ) subject to z ( i ) = z ( j ) , j ∈ [ m ] \ { i } ≤ d i − h Ti z ( i ) ≤ u i z ( i ) j ∈ D j , j ∈ [ n − m ] . (10)An integer linear inequality constraint α ≤ v ≤ β, where α, β, v ∈ Z , is equivalent to v ∈{ α, α + 1 , . . . , β } ⇔ ( v − α )( v − α − · · · ( v − β ) = 0 . Using this fact, we write the last twoconstraints in (10) in polynomial forms as follows. q i ( z ( i ) ) = ( d i − h Ti z ( i ) )( d i − h Ti z ( i ) − · · · ( d i − h Ti z ( i ) − u i ) = 0 , (11) r j ( z ( i ) ) = ( z ( i ) j − z j, )( z ( i ) j − z j, ) . . . ( z ( i ) j − z j,p i ) = 0 , j ∈ [ n − m ] . (12)Hence for all i ∈ [ m ] any feasible z ( i ) for problem (10) comes from the following set: F = m (cid:92) k =1 { z ∈ Z n − m | q k ( z ) = 0 , and ( ∀ j ∈ [ n − m ]) r j ( z ) = 0 } = { z ∈ Z n − m | ( ∀ k ∈ [ m ]) q k ( z ) = 0 , and ( ∀ j ∈ [ n − m ]) r j ( z ) = 0 } . (13)In (13), the intersection in the first line ensures that the consensus constraints are satisfied, andthe second line just expands the first. So the optimization problem (10) is equivalent tominimize z ( i ) ¯ f i (cid:16) z ( i ) (cid:17) subject to z ( i ) ∈ F , (14)for i ∈ [ m ] . Thus each of these players is optimizing over a common constraint set F . So findingthe points in F is of interest, which we discuss next. In this section, first, we review some necessary background on algebraic geometry, and then wepresent a theorem to check if F is nonempty and provide an algorithm to compute the points ina nonempty F . Finally, we present our algorithm to compute efficient Pareto optimal points. Indevising our algorithm, we use algebraic geometry rather than integer programming techniquesfor the following reasons. First, in this way, we are able to provide an algebraic geometriccharacterization for the set of efficient Pareto optimal solutions for our problem.
Second, we canshow that this set is nonempty if and only if the reduced Groebner basis (disucussed in Section4.1 below) of a certain set associated with the problem is not equal to { } (Theorem 3). Third, the mentioned result has an algorithmic significance: the reduced Groebner basis can be used toconstruct algorithms to calculate efficient Pareto optimal points.9.1 Background on algebraic geometryA monomial in variables x = ( x , x , . . . , x n ) is a product of the structure x α = x α · · · x α n n , where α = ( α , . . . , α n ) ∈ N n . A polynomial is an expression that is the sum of a finite numberof terms with each term being a monomial times a real or complex coefficient. The set of allreal polynomials in x = ( x , . . . , x n ) with real and complex coefficients are denoted by R [ x ] and C [ x ] , respectively, with the variable ordering x > x > · · · > x n . The ideal generated by f , f , . . . , f m ∈ C [ x ] is the set ideal { f , . . . , f m } = (cid:40) m (cid:88) i =1 h i f i | ( ∀ i ∈ [ m ]) h i ∈ C [ x ] (cid:41) . Consider f , f , . . . , f s which are polynomials in C [ x ] . The affine variety V of f , f , . . . , f m isgiven by V ( f , . . . , f m ) = { x ∈ C n | ( ∀ i ∈ [ m ]) f i ( x ) = 0 } . (15)A monomial order on C [ x , . . . , x n ] is a relation, denoted by (cid:31) , on the set of monomials x α , α ∈ N n satisfying the following. First , it is a total order; second , if x α (cid:31) x β and x γ isany monomial, then x α + γ (cid:31) x β + γ ; third , every nonempty subset of N n has a smallest elementunder (cid:31) . We will use lexicographic order , where we say x α (cid:31) lex x β if and only if the left mostnonzero entry of α − β is positive. Suppose that we are given a monomial order (cid:31) and apolynomial f ( x ) = (cid:80) α ∈ S f α x α . The leading term of the polynomial with respect to (cid:31) , denotedby lt (cid:31) ( f ) , is that monomial f α x α with f α (cid:54) = 0 , such that x α (cid:31) x β for all other monomials x β with f β (cid:54) = 0 . The monomial x α is called the leading monomial of f . Consider a nonzeroideal I . The set of the leading terms for the polynomials in I is denoted by lt (cid:31) ( I ) . Thuslt (cid:31) ( I ) = { cx α | ( ∃ f ∈ I ) lt (cid:31) ( f ) = cx α } . By ideal { lt (cid:31) ( I ) } with respect to (cid:31) , we denote theideal generated by the elements of lt (cid:31) ( I ) .A Groebner basis G (cid:31) of an ideal I with respect to the monomial order (cid:31) is a finite set ofpolynomials g , . . . , g t ∈ I such that ideal { lt (cid:31) ( I ) } = ideal { lt (cid:31) ( g ) , . . . , lt (cid:31) ( g t ) } . A reducedGroebner basis G reduced , (cid:31) for an ideal I is a Groebner basis for I with respect to monomialorder (cid:31) such that, for any f ∈ G reduced , (cid:31) , the coefficient associated with lt (cid:31) ( f ) is , and for all f ∈ G reduced , (cid:31) , no monomial of f lies in ideal { lt (cid:31) ( G \ { f } ) } . For a nonzero ideal I and givenmonomial ordering, the reduced Groebner basis is unique [Proposition 6, 37, Page 92]. Suppose I = ideal { f , . . . , f m } ⊆ C [ x , . . . , x n ] . Then for any l ∈ [ n ] , the l th elimination ideal for I isdefined by I l = I ∩ C [ x l +1 , . . . , x n ] . Let I ⊆ C [ x , . . . , x n ] be an ideal, and let G be a Groebnerbasis of I with respect to lexicographic order with x (cid:31) x (cid:31) . . . (cid:31) x n . Then for every integer l ∈ { , n − } , the set G l = G ∩ C [ x l +1 , . . . , x n ] is a Groebner basis for the l th elimination ideal I l .4.2 Nonemptyness of F We will use the following theorem in proving the results in this section.
Theorem 2 (Weak Nullstellensatz [Theorem 1, 37, page 170]) Consider f , f , . . . , f s as poly-nomials in C [ x , x , . . . , x n ] . If V ( f , f , . . . , f m ) = ∅ (see (15) for definition of V ), then ideal { f , f , . . . , f m } = C [ x , x , . . . , x n ] . First, we present the following result. 10 heorem 3
The set F is nonempty if and only if G reduced , (cid:31) (cid:54) = { } , where G reduced , (cid:31) is the reduced Groebner basis of ideal { q , . . . , q m , r , . . . , r n − m } with respect toany ordering.Remark 4 In the proof, we have shown that feasibility of the system (26) in C n − m is equivalentto its feasibility in Z n − m . So V ( q , . . . , q m , r , . . . , r n − m ) ∩ Z n − m = V ( q , . . . , q m , r , . . . , r n − m ) . (16)If we are interested in just verifying the feasibility of the polynomial system, then calculatinga reduced Groebner basis with respect to any ordering suffices. However, if we are interested inextracting the feasible points, then we choose lexicographic ordering as lexicographic orderingallows us to use algebraic elimination theory. There are many computer algebra packages to com-pute reduced Groebner basis such as Macaulry2, SINGULAR, FGb, Maple, and Mathematica.We now describe how to extract the points in F .Suppose G reduced , (cid:31) (cid:54) = { } . Naturally, the next question is how to compute points in F ? Inthe next section, we will show that the points in F are related to the efficient Pareto optimalpoints that we are seeking. We now briefly discuss systematic methods for extracting F basedon algebraic elimination theory, a branch of computational algebraic geometry. For details onelimination theory, we refer the interested readers to [37, Chapter 3]. First, we present thefollowing lemma. Lemma 5
Suppose G reduced , (cid:31) (cid:54) = { } . Then F = V ( G reduced , (cid:31) lex ) (cid:54) = ∅ . Algorithm 1 below calculates all the points in F , when G reduced , (cid:31) lex (cid:54) = { } . Lemma 6 provesits accuracy.
Lemma 6
Algorithm 1 correctly calculates all the points in F , when it is nonempty. F Suppose G reduced , (cid:31) lex (cid:54) = { } , and using Algorithm 1 we have computed F . We now proposeAlgorithm 2 and show that the resultant points are Pareto optimal.We have the following results for Algorithm 2. Lemma 7
In Algorithm 2, for all i ∈ [ m − , we have F ∗ s i +1 ⊆ F ∗ s i ⊆ F ( F ∗ s i defined in (19) ). Lemma 8
In Algorithm 2, for any i ∈ [ m ] , x s i ∈ X ∗ s i (the set X ∗ i is defined in Step 1, 3 ofAlgorithm 2) if and only if z ∗ ∈ F ∗ s i . Furthermore, z ∗ ∈ F ∗ s i solves the following optimizationproblem minimize z f s i ( d s i − h Ts i z )subject to z ∈ F ∗ s i − , for all i ∈ [2 : m ] . In Algorithm 2, at no stage can F ∗ s i get empty. Lemma 9
Suppose
F (cid:54) = ∅ . Then in Algorithm 2, F ∗ s i is nonempty for any i ∈ [ m ] . Theorem 4
For any z ∗ ∈ F ∗ s m , x ∗ = ( d − h T z ∗ , . . . , d m − h Tm z ∗ , z ∗ , . . . , z ∗ n − m ) (21) is an efficient Pareto optimal point. lgorithm 1 Extracting the points in F Input:
Polynomial system q i ( z ) = 0 for i ∈ [ m ] , and r j ( z ) = 0 for j ∈ [ n − m ] , G reduced , (cid:31) lex (cid:54) = { } . Output:
The set F .Step 1. – Calculate the set G n − m − = G reduced , (cid:31) lex ∩ C [ z n − m ] , which is a Groebner basis of the ( n − m ) thelimination ideal of ideal { q , . . . , q m , r , . . . , r n − m } that consists of univariate polynomials in z n − m as an implication of [37, page 116, Theorem 2]. – Find the variety of G n − m − , denoted by V ( G n − m − ) , which contains the list all possible z n − m coordinates for the points in F .Step 2. – Calculate G n − m − = G reduced , (cid:31) lex ∩ C [ z n − m − , z n − m ] , which is again a Groebner basis of the ( n − m − th elimination ideal of ideal { q , . . . , q m , r , . . . , r n − m } and consists of bivariate polynomialsin z n − m and z n − m − . – From Step 1, we already have the z n − m coordinates for the points in F . So by substituting those | V ( G n − m − ) | values in G n − m − , we arrive at a set of univariate polynomials in z n − m − , whichwe denote by { ¯ G ( i ) n − m − } | V ( G n − m − ) | i =1 . – For all i = 1 , , . . . , | V ( G n − m − ) | , find the variety of ¯ G ( i ) n − m − , denoted by V ( ¯ G ( i ) n − m − ) , whichcontains the list all possible z n − m − coordinates associated with a particular z n − m ∈ V ( G n − m − ) .We now have all the possible ( z n − m − , z n − m ) coordinates of F .Step 3. – We repeat this procedure for G n − m − , G n − m − , . . . , G . In the end, we have set of all points in F . return F . F is empty In the case that F is empty, we design a penalty based approach to solve a penalized version ofproblem (14) that can be of use to network administrators and policy makers. First, for i ∈ [ m ] ,using (13), (12), and (11), problem (14) can be written in the following equivalent form:minimize z ( i ) ∈ D ¯ f i (cid:16) z ( i ) (cid:17) subject to ( ∀ k ∈ [ m ]) q k ( z ( i ) ) = 0 . (22)When F = ∅ , we relate the problem above to a penalized version, which is a standardpractice in operations research literature [38, Page 278]. In this penalized version, we disregardthe equality constraints q k ( z ( i ) ) = 0 for k ∈ [ m ] , rather we we augment the cost with a termthat penalizes the violation of these equality constraints. So a penalized version of the problem(22) is as follows: minimize z ( i ) ¯ f i (cid:16) z ( i ) (cid:17) + m (cid:88) k =1 γ k p (cid:16) q k ( z ( i ) ) (cid:17) , subject to z ( i ) ∈ D. (23)12 lgorithm 2 Computing the set of solutions to problem (14).
Input:
The optimization problem (14) for any i ∈ [ m ] , F (cid:54) = ∅ . Output:
Efficient Pareto optimal solutions for problem (3).Step 1. for i = 1 , . . . , m X i := (cid:110) d i − h Ti z ( i ) | z ( i ) ∈ F (cid:111) = d i − h Ti F , ( ∀ x i ∈ X i ) ( X i ) − ( x i ) := { z ( i ) ∈ F | x i = d i − h Ti z ( i ) . } (17) end for Step 2.Sort the elements of the { X i } mi =1 s with respect to cardinality of the elements in a descending order.Denote the index set of the sorted set by { s , s , . . . , s m } such that | X | s ≥ | X | s ≥ · · · ≥ | X | s m . Step 3. for i ∈ [ m ] Solve the univariate optimization problem minimize x si f s i ( x s i )subject to x s i ∈ X s i , (18)and denote the set of solutions by X ∗ s i . Set F ∗ s i := (cid:91) x si ∈ X ∗ si ( X ∗ s i ) − ( x s i ) ⊆ F , (19) if i ≤ m X s i +1 := (cid:110) d s i +1 − h Ts i +1 z | z ∈ F ∗ s i (cid:111) . (20) end ifend forreturn F ∗ s m . where p : R → R + is a penalty function, and γ k is a positive penalty parameter. Some commonpenalty functions are: – exact penalty: p : x (cid:55)→ x , – power barrier: p : x (cid:55)→ | x | , etc. We have already shown D to be a nonempty set. From a network point of view, the penalizedproblem (23) has the following interpretation. For the first m arcs in the directed network underconsideration, rather than having a strict flow bound constraint (see the discussion in Section13), we have a penalty when x i < or x i > u i for i ∈ [ m ] . The flow bound constraint is stillmaintained for i ∈ [ m + 1 : n ] . In this regard, the original problem defined by (3) has thefollowing penalized version:minimize x i f i ( x i ) + (cid:18) penalty for violation of ≤ x i ≤ u i for player i ∈ [ m ] (cid:19) subject to Ax = A ( x i , x − i ) = A i x i + n (cid:88) j =1 ,j (cid:54) = i A j x j = b ≤ x i ≤ u i , i ∈ [ m + 1 : n ] x ∈ Z n . (24)Problem (24) can be considered as a network flow problem, where flow direction is flexibleor overflow is allowed subject to a suitable penalty, which is a quite realistic scenario in practice[7, pages 550-551]. With this penalized problem, we can proceed as follows. In the developmentsof Section 4.3 set: ¯ f i (cid:16) z ( i ) (cid:17) := ¯ f i (cid:16) z ( i ) (cid:17) + m (cid:88) k =1 γ k p (cid:16) q k ( z ( i ) ) (cid:17) , F := D (cid:54) = ∅ , and then apply Algorithm 2, which will calculate efficient Pareto optimal points for the penalizedproblem (24).The described penalty scheme can be of use to network administrators and policy makers toenforce a decision making architecture. Such an architecture would allow the players to makeindependent decisions while ensuring that (i) total amount of flow is conserved in the networkby maintaining the mass balance constraint , (ii) the flow bound constraint is strictly enforcedfor the last n − m players and is softened for the first m players by imposing penalty, and yet (iii) an efficient Pareto optimal point for the penalized problem can be achieved by the playerswhere none of their objective functions can be improved without worsening some of the otherplayers’ objective values.From both the players’ and the policy maker’s point of views, the penalty based approachmakes sense and can be considered fair for the following reason. In the exact version, each ofthe last n − m players gets to minimize its optimization problem (9) in a decoupled manner,whereas each of the first m players is solving a more restrictive optimization problem (18). Socutting each of the first m players some slack by softening the flow bound constraint, where itcan carry some extra flow or flow in opposite direction by paying a penalty, can be consideredfair. In this section, we present an illustrative numerical example of our methodology in a transporta-tion setup. We have used
Wolfram Mathematica 10 for numerical computation.
Problem setup.
Consider the following directed network with 5 nodes and 16 arcs as shown inFigure 1. Nodes 2 and 4 represent two retail centers with demands for 13 and 11 units of acertain product. The warehouses are denoted by nodes 1 and 3, which supply 9 and 15 units,respectively. Node 5 is a trans-shipment node. Different modes of shipment from one node to14 ig. 1:
Network under consideration other is represented by the arcs in the figure, and these shipments are carried out by differentorganizations (carriers). The cost of a certain shipment depends on the number of productsshipped; it is nonlinear and not necessarily convex. With each arc, we associate one carrier(player).The cost functions for the players are listed in Table 1. The cost functions associated withplayers 5, 7, 10, and 15 are convex, and the cost functions associated with the rest of the playersare nonconvex. Each of the players is trying to minimize its cost. The number of products carriedby each player has to be non-negative, and the capacity of each player is represented by u = (10 , , , , , , , , , , , , , , , , where u i represents the maximum capacity for the number of products carried by player i .We seek efficient Pareto optimal points in this setup. Computing node-arc incidence matrix and resource vector.
First, we compute the node-arc inci-dence matrix associated with the network under consideration by following the procedure men-tioned in the description of the mass balance constraint in Section 2. The resultant augmentednode-arc incidence matrix of the network is: ˜ A = − − − − − − − − − − − − − − − − layer Cost function1 − x − x + x − x + 1 x − x + x − x + x + 1 x − x + x − x + 485 x − x + x − x + x + 5 ( x − − x + x − x + x + 10 | x − | x − x + x − x + x − x + x + 1 − x + x − x + ( x − x − x + x − x + 90 x − x + x − x + 15 x − x + x − x + 285 x − x + x − x + 4165 | x − | x + 1 , if ≤ x ≤ , if ≤ x ≤ x + 1) , if ≤ x ≤ − x + x − x + 330 , else Table 1:
Cost functions for the players in the network considered
The vector ˜ b , where ˜ b i represents the demand or supply at node i , can be computed byrecording the given demand or supply at each node and is given by ˜ b = (9 , − , , − , . Asdiscussed in the description of the mass balance constraint in Section 2, we compute the full rowrank node-arc incidence matrix A by removing the 5th row of ˜ A and the resource vector b byremoving the last component of ˜ b , i.e. , b = (9 , − , , − . So A has 4 rows and 16 columns,and it is: A = − − − − − − − − − − − − , where we associate player i with the i th column of A . Reindexing the variables.
To obtain the notational advantage as discussed in the comment afterLemma 1, we next rearrange the columns of the matrix A so that the unimodular basis matrixcomprising of the last 4 columns of A constitutes the first columns in the new A i.e. , if A = [ A | A | . . . | A | A | . . . | A ] then det([ A | A | . . . | A ]) = ± . We reindex thevariables as follows. The last 4 players are reindexed as the first four players, and players 1 to12 are reindexed as players 5 to 16. Due to this reindexing procedure, we need to reindex theupper bound vector u , where the last 4 indices become the first 4 indices, and the indices 1 to12 become indices 5 to 16, i.e. , u = (5 , , , , , , , , , , , , , , , . b would not change due to the reindexing procedure. We perform ourcomputation on these reindexed variables and then revert the resultant efficient Pareto optimalsolutions to the original indexing. Problem transformation.
Now we are in a position to transform the problem by following themethodology described in Section 3. First, we decouple the optimization problems for the last12 players. In order to achieve that, we follow the constructive proofs of Lemma 4 and Theorem1 to arrive at the following representation of the variable x in terms of the new variable z : x = x x x x x x x x x x x x x x x x = z + z + z − z − z − z + 9 − z + z + z + z − z − − z + z + z + z − z + 15 − z − z + z + z + z − z z z z z z z z z z z z , (25)where the advantage of this transformation is that, for players 5 to 16, we have decoupledunivariate optimization problems of the form (9) in z , and by solving these decoupled problems,we can reduce the constraint set significantly. Solving the decoupled optimization problems.
Next we solve the aforementioned decoupled uni-variate optimization problems for the last 12 players by following the procedure described inthe last paragraph of Subsection 3.1. The solution set is given by Table 2; in this table, the setsof optimal solutions for players , , . . . , are given by D , D , . . . , D , respectively. Solvingthese 12 decoupled optimization problems immediately reduces the constraint set into a muchsmaller set, which we define as D = × i =1 D i . Consensus reformulation for the first 4 players.
Now we are in a position to transform theoptimization problems for the first 4 players using (8) and consensus constraints as describedin Subsection 3.2. By following the straightforward procedure of Subsection 3.2, we arrive at 4optimization problems of the form (14), where each of the players 1 to 4 is optimizing its costfunction over a common constraint set denoted by F . Next we compute the points in F . Computing the points in F . To compute the points in F , we follow the methodology describedin Subsection 4.2. Finding the points in F requires computing the reduced Groebner basis withrespect to lexicographic ordering, denoted by G reduced , (cid:31) lex . We compute G reduced , (cid:31) lex using the GroebnerBasis function in
Wolfram Mathematica 10 , and we find it to be not equal to { } .Hence F is nonempty due to Lemma 5. Next we compute the points in F by following Algorithm1; the list of computed points is given by Table 3.17 { } D { } D { } D { , } D { , } D { } D { } D {1} D { } D { } D { } D { , , , , } Table 2: D for the network under considerationElements z z z z z z z z z z z z Table 3:
Table listing elements of F Computing the efficient Pareto optimal points.
The final step is computing the efficient Paretooptimal points from F by executing the three steps described in Algorithm 2. After applyingAlgorithm 2, we arrive at two efficient Pareto optimal points in variable z as follows: (1 , , , , , , , , , , , , and (1 , , , , , , , , , , , . Using the relationship between x and z in (25), we can express these efficient Pareto optimalpoints in variable x as follows: (5 , , , , , , , , , , , , , , , , and (3 , , , , , , , , , , , , , , , . Note that these efficient Pareto optimal points above are associated with the reindexed x .By reversing this reindexing, where indices 5 to 16 will be indices 1 to 12, and indices 1 to 4 willbe indices 13 to 16, we arrive at the efficient Pareto optimal points in our original variable x asfollows: (1 , , , , , , , , , , , , , , , , and (1 , , , , , , , , , , , , , , , . Conclusions
In this paper, we have proposed a multi-player extension of the minimum cost flow probleminspired by a multi-player transportation problem. We have associated one player with each arcof a directed connected network, each trying to minimize its cost subject to the network flowconstraints. The cost can be any general nonlinear function, and the flow through each arc isinteger-valued. In this multi-player setup, we have presented algorithms to compute an efficientPareto optimal point , which is a good compromise solution between the utopian vector optimalityand the generic Pareto optimality and is a Nash equilibrium if our problem is transformed intoan n -person finite static game in normal form.Some concluding remarks on the limitations of our methodology are as follows. First , at theheart of our methodology is the transformation provided by Theorem 1, which decouples theoptimization problems for the last n − m players. Each of these decoupled optimization problemsis univariate over a discrete interval and is easy to solve. This can potentially allow us to workin a much smaller search space. So if we have a system where n − m > m ⇔ m < n , then it willbe convenient from a numerical point of view. Second , computation of Pareto optimal pointsdepends on determining the points in F using Groebner basis. Calculating Groebner basis canbe numerically challenging for large system [37, pages 111-112], though significant speed-up hasbeen achieved in recent years by computer algebra packages such as Macaulry2, SINGULAR,FGb, and Mathematica. References
1. Shuvomoy Das Gupta and Lacra Pavel. Multi-player minimum cost flow problems with non-convex costs and integer flows. In
Decision and Control (CDC), 2016 IEEE 55th Conferenceon , pages 7617–7622. IEEE, 2016.2. Efraim Turban, Jon Outland, David King, Jae Kyu Lee, Ting-Peng Liang, and Debor-rah C Turban.
Electronic Commerce 2018: A Managerial and Social Networks Perspective .Springer, Berlin, 2017.3. Scott Galloway.
The Four: The Hidden DNA of Amazon, Apple, Facebook, and Google .Portfolio, New York, 2017.4. Ling Li.
Managing Supply Chain and Logistics: Competitive Strategy for a Sustainable Future .World Scientific Publishing Co Inc, 2014.5. Victor M Mendez, Carlos A Monje Jr, and Vinn White. Beyond traffic: Trends and choices2045: A national dialogue about future transportation opportunities and challenges. In
Disrupting Mobility , pages 3–20. Springer, Cham, 2017.6. Timothy M Laseter and Elliot Rabinovich.
Internet retail operations: integrating theory andpractice for managers . CRC Press, 2011.7. Stephen Boyd and Lieven Vandenberghe.
Convex optimization . Cambridge university press,2004.8. Kaisa Miettinen.
Nonlinear multiobjective optimization , volume 12. Springer Science &Business Media, 2012.9. Ravindra K Ahuja, Thomas L Magnanti, and James B Orlin.
Network Flows: Theory,Algorithms, and Applications . Prentice Hall, 1993.190. Dimitris Bertsimas, Ebrahim Nasrabadi, and Sebastian Stiller. Robust and adaptive networkflows.
Operations Research , 61(5):1218–1242, 2013.11. Balachandran Vaidyanathan and Ravindra K Ahuja. Fast algorithms for specially structuredminimum cost flow problems with applications.
Operations research , 58(6):1681–1696, 2010.12. Thomas L Magnanti and Richard T Wong. Network design and transportation planning:Models and algorithms.
Transportation science , 18(1):1–55, 1984.13. Stephen C Graves and James B Orlin. A minimum concave-cost dynamic network flowproblem with an application to lot-sizing.
Networks , 15(1):59–71, 1985.14. Mark S Daskin.
Network and discrete location: models, algorithms, and applications . JohnWiley & Sons, 2011.15. Stefan Feltenmark and Per Olof Lindberg.
Network methods for head-dependent hydro powerscheduling . Springer, 1997.16. Bt Yaged. Minimum cost routing for static network models.
Networks , 1(2):139–172, 1971.17. Qie He, Shabbir Ahmed, and George L Nemhauser. Minimum concave cost flow over a gridnetwork.
Mathematical Programming , 150(1):79–98, 2015.18. Hoang Tuy, Saied Ghannadan, Athanasios Migdalas, and Peter Värbrand. The minimumconcave cost network flow problem with fixed numbers of sources and nonlinear arc costs.
Journal of Global Optimization , 6(2):135–151, 1995.19. Dalila BMM Fontes, Eleni Hadjiconstantinou, and Nicos Christofides. A branch-and-boundalgorithm for concave network flow problems.
Journal of Global Optimization , 34(1):127–155,2006.20. Willard I Zangwill. Minimum concave cost flows in certain networks.
Management Science ,14(7):429–450, 1968.21. S Jorjani and BW Lamar. Cash flow management network models with quantity discounting.
Omega , 22(2):149–155, 1994.22. David Gamarnik, Devavrat Shah, and Yehua Wei. Belief propagation for min-cost networkflow: Convergence and correctness.
Operations research , 60(2):410–428, 2012.23. Dukwon Kim and Panos M Pardalos. A dynamic domain contraction algorithm for noncon-vex piecewise linear network flow problems.
Journal of Global optimization , 17(1-4):225–234,2000.24. Shangyao Yan, Yu-Lin Shih, and Wang-Tsang Lee. A particle swarm optimization-basedhybrid algorithm for minimum concave cost network flow problems.
Journal of Global Op-timization , 49(4):539–559, 2011.25. Andrea Raith and Matthias Ehrgott. A two-phase algorithm for the biobjective integerminimum cost flow problem.
Computers & Operations Research , 36(6):1945–1954, 2009.26. Haijune Lee and P Simin Pulat. Bicriteria network flow problems: Integer case.
EuropeanJournal of Operational Research , 66(1):148–157, 1993.27. Augusto Eusébio and José Rui Figueira. Finding non-dominated solutions in bi-objectiveinteger network flow problems.
Computers & Operations Research , 36(9):2554–2564, 2009.208. Salvador Hernández, Srinivas Peeta, and George Kalafatas. A less-than-truckload carriercollaboration planning problem under dynamic capacities.
Transportation Research Part E:Logistics and Transportation Review , 47(6):933–946, 2011.29. Cynthia Barnhart, Christopher A Hane, and Pamela H Vance. Integer multicommodity flowproblems. In
Integer Programming and Combinatorial Optimization , pages 58–71. Springer,Berlin, 1996.30. Lorenzo Brunetta, Michele Conforti, and Matteo Fischetti. A polyhedral approach to aninteger multicommodity flow problem.
Discrete Applied Mathematics , 101(1):13–36, 2000.31. Asuman E Ozdaglar and Dimitri P Bertsekas. Optimal solution of integer multicommodityflow problems with application in optical networks. In
Frontiers in global optimization , pages411–435. Springer, Boston, 2004.32. Dimitris Bertsimas and John N Tsitsiklis.
Introduction to linear optimization , volume 6.Athena Scientific Belmont, MA, 1997.33. Tamer Basar and Geert Jan Olsder.
Dynamic noncooperative game theory , volume 23. SIAM,1995.34. Alexander Schrijver.
Theory of linear and integer programming . John Wiley & Sons, 1998.35. Dimitris Bertsimas and Robert Weismantel.
Optimization over integers , volume 13. DynamicIdeas Belmont, 2005.36. Neal Parikh and Stephen Boyd. Proximal algorithms.
Foundations and Trends in optimiza-tion , 1(3):123–231, 2013.37. David Cox, John Little, and Donal O’shea.
Ideals, varieties, and algorithms , volume 3.Springer, 2007.38. Dimitri P. Bertsekas, Asuman E. Ozdaglar, and Angelia Nedic.
Convex analysis and opti-mization . Athena scientific optimization and computation series. Athena Scientific, Belmont(Mass.), 2003. ISBN 1-886529-45-0. 21 ppendixProof of Lemma 1
Any basic solution of relaxed P can be constructed as follows. Take m linearly independentcolumns A B (1) , A B (2) , . . . , A B ( m ) , where B (1) , B (2) , . . . , B ( m ) are the indices of those columns. Construct the basis matrix B = (cid:2) A B (1) A B (2) · · · A B ( m ) (cid:3) . Set x B = ( x B (1) , x B (2) , . . . , x B ( m ) ) = B − b , which is called the basic variable in linear pro-gramming theory. Set the rest of the components of x (called the nonbasic variable) to be zero, i.e. , x NB = ( x NB (1) , x NB (2) , . . . , x NB ( n − m ) ) = (0 , , . . . , . The resultant x will be a basic solution of relaxed P [32, Theorem 2.4].Suppose there is an integer basic solution in relaxed P . Denote that integer basic solutionby ¯ x and the associated basis matrix by ¯ B . The nonbasic variables are integers (zeros) in everybasic solution, so the basic variables ¯ x ¯ B = ¯ B − b has to be an integer vector for any integer b .Now from Cramer’s rule [35, Proposition 3.1], ( ∀ b ∈ Z m ) (cid:18) ¯ x ¯ B = ¯ B − b ∈ Z m ⇔ ∀ i ∈ [ m ] ¯ x ¯ B ( i ) = det ¯ B i det ¯ B ∈ Z (cid:19) , where ¯ B i is the same as ¯ B , except the i th column has been replaced with ¯ b . Now det ¯ B i ∈ Z because b is an integer vector. So having integer basic solution is equivalent to det ¯ B = ± , i.e. , ¯ B is unimodular. (cid:4) Proof of Lemma 2
First, note that unimodularity of B is equivalent to unimodularity of B − , which we showeasily as follows. First note that, det (cid:0) B − (cid:1) = B = ± = ± . Each component of B − is asubdeterminant of B divided by det B . As B ∈ Z m × m and det B = ± , each component of B − is an integer. Thus B − is a unimodular matrix. Similarly, unimodularity of B − implies B isunimodular.So C = B − A ∈ Z m × n and d = B − b ∈ Z m as A and b are integer matrices. As multiplyingboth sides of a linear system by an invertible matrix does not change the solution, we have Ax = b ⇔ B − Ax = B − b ⇔ Cx = d . Thus we arrive at the claim. (cid:4) Proof of Lemma 3
First, recall that the following operations on a matrix are called elementary integer columnoperations : (i) adding an integer multiple of one column to another column, (ii) exchanging twocolumns, and (iii) multiplying a column by -1.The matrix C is of the form: C = B − A = B − [ B | A m +1 | A m +2 | · · · | A n ]= · · · C ,m +1 C ,m +2 · · · C ,n · · · C ,m +1 C ,m +2 · · · C ,n ... ... ... ... ... ... ... ... · · · C m,m +1 C m,m +2 · · · C m,n C . For all j = m + 1 , m + 2 , . . . , n , we multiply the first column of C , C = e of C by − C ,j and then add it to C i . Thus C is transformed to · · · · · ·
00 1 · · · C ,m +1 C ,m +2 · · · C ,n ... ... ... ... ... ... ... ... · · · C m,m +1 C m,m +2 · · · C m,n Similarly for column indices, i = 2 , , . . . , m , respectively we do the following. For j = m +1 , m +2 , · · · , n , we multiply the i th column e i with − C i,j and add it to C j . In the end, the Hermitenormal form becomes: · · · · · ·
00 1 · · · · · · ... ... ... ... ... ... ... ... · · · · · · . The steps describing the process in Lemma 3 can be summarized by Algorithm 3, which weare going to use to prove Lemma 4. (cid:4)
Algorithm 3
Converting C to [ I | procedure Converting C to [ I | for i := 1 , , . . . , m do for j := m + 1 , m + 2 , . . . , n do C j := C j − C i,j e i end for end for end procedure Proof of Lemma 4
If we multiply column C i of a matrix C with an integer factor γ and subsequently add itto another column C j , then it is equivalent to right multiplying the matrix C with a matrix I + γI ij (recall that the matrix I ij has a one in ( i, j ) th position and zero everywhere else). Notethat I + γI ij is a triangular matrix with diagonal entries being one, γ being on the ( i, j ) thposition, and zero everywhere else. As the determinant of a triangular matrix is the product ofits diagonal entries, det( I + γI ij ) = 1 . So I + γI ij is a unimodular matrix. Furthermore, step 4of the procedure above to convert C to Hermite normal form, i.e. , C j = C j − C i,j e i is equivalentto left multiplying the current matrix with I − C i,j I ij . So the inner loop of the procedure aboveover j = m + 1 , m + 2 , . . . , n (lines 2-4) can be achieved by left multiplying the current matrixwith U i = n (cid:89) j = m +1 ( I − C ij I ij ) = ( I − C i,m +1 I i,m +1 )( I − C i,m +2 I i,m +2 ) · · · ( I − C i,n I i,n ) As each of the matrices in the product is a unimodular matrix and determinant of multiplicationof square matrices of same size is equal to multiplication of determinant of those matrices, wehave det( U i ) = 1 . So U i is a unimodular matrix. Structurally the i th row of U i , denoted by u Ti ,has a on i th position, has − C i,j on j th position for j = m + 1 , m + 2 , . . . , n , and zero every-where else. Any other k th row ( k (cid:54) = i ) of U i is e Tk . So we can convert C to its Hermite normalform by repeatedly left multiplying C with U , U , . . . , U m , respectively. This is equivalent toleft multiplying C with one single matrix U = (cid:81) mi =1 U i . The final matrix U is unimodular as it23s multiplication of unimodular matrices. (cid:4) Proof of Theorem 1
Let y = U − x . As U is unimodular, U − is also unimodular. Hence x ∈ Z n ⇔ y ∈ Z n . Let y = ( y , y ) , where y ∈ Z m and y ∈ Z n − m . Then Q (cid:54) = ∅⇔∃ x ∈ Z n ( Cx = d ) ⇔∃ y ∈ Z n (cid:18) CU y = (cid:2) I (cid:3) (cid:20) y y (cid:21) = y = d (cid:19) . As d = B − b ∈ Z m (Lemma 2), by taking y = (cid:20) y y (cid:21) = (cid:20) B − bz (cid:21) ∈ Z n , where z ∈ Z n − m , wecan satisfy the condition above. Thus Q is nonempty.Now for any x , we can maximally decompose it in terms of z as: x ∈ Q ⇔ x = U y = (cid:20) I m − B − A [1: m,m +1: n ] n − m × m I n − m . (cid:21) (cid:20) B − bz (cid:21) = (cid:20) B − (cid:0) b − A [1: m,m +1: n ] z (cid:1) z (cid:21) = d − h T zd − h T z ... d m − h Tm zz ... z n − m , where the last n − m variables are completely decoupled in terms of z . Note that the decompo-sition is maximal by construction and cannot be extended any further in general cases. (cid:4) Proof of Theorem 3
The proof sketch is as follows. The elements of F are the solution of the polynomial system: ( ∀ i ∈ [ m ]) q i ( z ) = 0( ∀ j ∈ [ n − m ]) r j ( z ) = 0 . (26)We prove in two steps. In step , we show that the polynomial system (26) is feasible if and onlyif / ∈ ideal { q , . . . , q m , r , . . . , r n − m } . Then in step , we show that, / ∈ ideal { q , . . . , q m , r , . . . , r n − m } is equivalent to G reduced , (cid:31) (cid:54) = { } . Step 1. Polynomial system (26) is feasible if and only if / ∈ ideal { q , . . . , q m , r , . . . , r n − m } . We prove necessity first and then sufficiency.( ⇒ ) 24ssume system (26) is feasible, but ∈ ideal { q , . . . , q m , r , . . . , r n − m } . That means ( ∃ h , . . . , h m , s , . . . , s n − m ∈ C [ z ]) (cid:0) ∀ z ∈ C n − m (cid:1) m (cid:88) i =1 h i ( z ) q i ( z ) + n − m (cid:88) i =1 s i ( z ) r i ( z ) , (27)and (cid:0) ∃ ¯ z ∈ Z n − m (cid:1) (cid:18) ( ∀ i ∈ [ m ]) q i (¯ z ) = 0 , ( ∀ i ∈ [ n − m ]) r i (¯ z ) = 0 (cid:19) . (28)Putting z = ¯ z in (27) and using (28) we get , so we have a contradiction.( ⇐ )We want to show that if (cid:54)∈ ideal { q , · · · , q m , r , · · · , r n − m } , then the polynomial system(26) is feasible. We prove the contrapositive again: if the polynomial system (26) is infeasible inintegers, then ∈ ideal { q , · · · , q m , r , · · · , r n − m } .First, we show that, feasibility of the system in C n − m is equivalent to feasibility in Z n − m . As Z n − m ⊂ C n − m , if the polynomial system is infeasible in C n − m , it is infeasible in Z n − m . Also,if the polynomial system is infeasible in Z n − m , then it will be infeasible in C n − m . We showthis by contradiction. Assume the system system is infeasible in Z n − m but feasible in C n − m i.e. , C n − m \ Z n − m . Let that feasible solution be ˜ z ∈ C n − m \ Z n − m , so there is at least one componentof it (say ˜ i ) such that ˜ z ˜ i ∈ C \ Z . Now r ˜ i (˜ z ) = (˜ z ˜ i − z i, )(˜ z ˜ i − z i, ) . . . (˜ z ˜ i − z i,p i ) , where, by construction, each of the elements of the set D i = { z i, , z i, , . . . , z i,p i } are integers anddifferent from each other, so each component in the product r ˜ i (˜ z ) are nonzero complex numberswith the absence of complex conjugates. So r ˜ i (˜ z ) (cid:54) = 0 , which is a contradiction.If the polynomial system is infeasible in C n − m , then it is equivalent to saying that, V ( q , . . . , q m , r , . . . , r n − m ) = ∅ , where V has been defined in (15). Then using the Weak Nullstellensatz , we have ideal { q , . . . , q m , r , . . . , r n − m } = C [ z , z , . . . , z n − m ] . As ∈ C [ z , z , . . . , z n − m ] , this means ∈ ideal { q , . . . , q m , r , . . . , r n − m } . Step 2. Now we show / ∈ ideal { q , . . . , q m , r , . . . , r n − m } is equivalent to G reduced , (cid:31) (cid:54) = { } . Wehave shown that feasibility of the system in C n − m is equivalent to feasibility in Z n − m . As a result,we can work over complex numbers, which is a algebraically closed field. This allows us to apply consistency algorithm [37, page 172], which states that / ∈ ideal { q , . . . , q m , r , . . . , r n − m } ifand only if G reduced , (cid:31) (cid:54) = { } . (cid:4) Proof of Lemma 5
By Theorem 3, we have
F (cid:54) = ∅ . So from the definition of affine variety in (15) and (13) we canwrite F as: F = V ( q , . . . , q m , r , . . . , r n − m ) ∩ Z n − m . F = V ( q , . . . , q m , r , . . . , r n − m ) (cid:54) = ∅ . By definition of basis, we have ideal { q , . . . , q m , r , . . . , r n − m } = ideal { G reduced , (cid:31) lex } , whichimplies V ( q , . . . , q m , r , . . . , r n − m ) = V ( G reduced , (cid:31) lex ) due to [Proposition 2, 37, page 32]. So F = V ( G reduced , (cid:31) lex ) . (cid:4) Proof of Lemma 6
Using Theorem 3, we have G reduced , (cid:31) lex (cid:54) = ∅ . So by the elimination theorem [Theorem 2, 37,page 116] V ( G n − m − ) is nonempty and will contain the list of all possible z n − m coordinates forthe points in F . As G reduced , (cid:31) lex (cid:54) = ∅ , when moving from one step to the next, not all the affinevarieties associated with the univariate polynomials (after replacing the previous coordinatesinto the elimination ideal) can be empty due to the extension theorem [Theorem 3, 37, page118]. Using this logic repeatedly, the final step will give us F = V ( G reduced , (cid:31) lex ) (cid:54) = ∅ . (cid:4) Proof of Lemma 7
Follows from (19), (17) and (20). (cid:4)
Proof of Lemma 8
For any i ∈ [ m ] , x s i ∈ X ∗ s i ⇔ d s i − h Ts i z ∈ X ∗ s i ⇔ z ∈ F ∗ s i , where the second line follows from (8). So (cid:32) minimize x si f s i ( x s i )subject to x s i ∈ X s i (cid:33) = (cid:32) minimize x si f s i ( d s i − h Ts i z )subject to z ∈ F ∗ s i − (cid:33) , where the second line follows from (20) in Algorithm 2. (cid:4) Proof of Lemma 9 As F (cid:54) = ∅ , X s := (cid:8) d s − h Ts z ( s ) | z ( s ) ∈ F (cid:9) (cid:54) = ∅ . Assume, for i ∈ [ m ] , we have F ∗ s i (cid:54) = ∅ .Then X s i +1 := (cid:110) d s i +1 − h Ts i +1 z | z ∈ F ∗ s i (cid:111) = d s i +1 − h Ts i +1 F ∗ s i (cid:54) = ∅ . The subsequent optimizationproblem is minimize x si +1 f s i +1 ( x s i +1 )subject to x s i +1 ∈ X s i +1 . As we are optimizing over a finite and countable set, a minimizer will exist. So X ∗ s i +1 (cid:54) = ∅ . Hence F ∗ s i = (cid:83) x i ∈ X ∗ si ( X ∗ s i ) − ( x i ) (cid:54) = ∅ . So for any i = 1 , . . . , m , we have F ∗ s i nonempty. (cid:4) Proof of Theorem 4
26e want to show that x ∗ is feasible, and for every x ∈ P, i ∈ [ n ] , if we have f i ( x ∗ i ) ≥ f i ( x i ) , thenfor every j ∈ [ n ] , f j ( x ∗ j ) = f j ( x j ) . Using (21) and (8), we can translate the Pareto optimalitycondition in z as follows. Consider a z ∈ Z n − m such that (0 , . . . , (cid:22) (cid:16)(cid:0) d i − h Ti z (cid:1) mi =1 (cid:17) , z ) (cid:22) ( u , . . . , u n ) , (29)and suppose (cid:0)(cid:0) f i ( d i − h Ti z ∗ ) mi =1 (cid:1) , ( f m + i ( z ∗ i )) n − mi =1 (cid:1) (cid:23) (cid:0)(cid:0) f i ( d i − h Ti z ) mi =1 (cid:1) , ( f m + i ( z i )) n − mi =1 (cid:1) . (30)Then we want to show that: (cid:0)(cid:0) f i ( d i − h Ti z ∗ ) mi =1 (cid:1) , ( f m + i ( z ∗ i )) n − mi =1 (cid:1) = (cid:0)(cid:0) f i ( d i − h Ti z ) mi =1 (cid:1) , ( f m + i ( z i )) n − mi =1 (cid:1) (31)Let us start with the last n − m rows. As, z ∗ ∈ F ∗ s m ⊆ F ⊆ D and by construction, D = × n − mi =1 D i , where any element of D i is a minimizer of (9), so ( f m + i ( z ∗ i )) n − mi =1 (cid:23) ( f m + i ( z i )) n − mi =1 implies ( f m + i ( z ∗ i )) n − mi =1 = ( f m + i ( z i )) n − mi =1 . In the subsequent steps, it suffices to confine z ∈ D , as otherwise last n − m inequalitiesof (30) will be violated. Now let us consider the first m inequalities of (30) As discussed inSection 3.2, any z is in D which obeys the inequalities ≤ d i − h Ti z ≤ u i for i = 1 , . . . , m ; thisis equivalent to z ∈ F ⊆ D . Consider, s ∈ [ m ] . Lemmas 7 and 8 implies that z ∗ solves thefollowing optimization problem min z { f s ( d s − h Ts z ) | z ∈ F } = min x s { f s ( x s ) | x s ∈ X s } ,which has solution x ∗ s ∈ X ∗ s ⇔ z ∗ ∈ F ∗ s ⊇ F ∗ s m . So f s ( d s − h Ts z ) ≤ f s ( d s − h Ts z ∗ ) implies f s ( d s − h Ts z ) = f s ( d s − h Ts z ∗ ) and z ∈ F ∗ s .Consider s ∈ [ m ] \ { s } . First, note that z ∈ F ∗ s , otherwise f s ( d s − h Ts z ) ≤ f s ( d s − h Ts z ∗ ) will not hold. Now the x ∗ s associated with z ∗ solves the following optimization problem min x s { f s ( x s ) | x s ∈ X s } = min z { f s ( d s − h Ts z ) | z ∈ F ∗ s } , where an optimal solution to thefirst line will be in X ∗ s and the optimal solution to the second line will be in F ∗ s (Lemma 8). Socombining f s ( d s − h Ts z ) ≤ f s ( d s − h Ts z ∗ ) and z ∈ F ∗ s implies f s ( d s − h Ts z ) = f s ( d s − h Ts z ∗ ) . Repeating a similar argument for i = s , s , . . . , s m , we can show that ( ∀ i ∈ { s , . . . , s m } ) f s ( d s − h Ts z ) = f s i ( d s i − h Ts i z ∗ ) . Thus we have arrived at (31). (cid:4)(cid:4)