Using Inverse Optimization to Learn Cost Functions in Generalized Nash Games
UUsing Inverse Optimization to Learn Cost Functions inGeneralized Nash Games
STEPHANIE ALLEN, UNIVERSITY OF MARYLAND, COLLEGE PARKJOHN P. DICKERSON, UNIVERSITY OF MARYLAND, COLLEGE PARKSTEVEN A. GABRIEL, UNIVERSITY OF MARYLAND, COLLEGE PARK & NORWEGIANUNIVERSITY OF SCIENCE AND TECHNOLOGY
As demonstrated by Ratliff et al. [65], inverse optimization can be used to recover the objective functionparameters of players in multi-player Nash games. These games involve the optimization problems of multipleplayers in which the players can affect each other in their objective functions. In generalized Nash equilibriumproblems (GNEPs), a player’s set of feasible actions is also impacted by the actions taken by other players inthe game; see Facchinei and Kanzow [25] for more background on this problem. One example of such impactcomes in the form of joint/“coupled” constraints as referenced by Rosen [67], Harker [36], and Facchineiet al. [24] which involve other players’ variables in the constraints of the feasible region. We extend theframework of Ratliff et al. [65] to find inverse optimization solutions for the class of GNEPs with jointconstraints. The resulting formulation is then applied to a simulated multi-player transportation problem ona road network. Also, we provide some theoretical results related to this transportation problem regardingruntime of the extended framework as well as uniqueness and nonuniqueness of solutions to our simulationexperiments. We see that our model recovers parameterizations that produce the same flow patterns as theoriginal parameterizations and that this holds true across multiple networks, different assumptions regardingplayers’ perceived costs, and the majority of restrictive capacity settings and the associated numbers of players.Code for the project can be found at: https://github.com/sallen7/IO_GNEP
Working paper. Contact Stephanie Allen [email protected] with questions. a r X i v : . [ m a t h . O C ] F e b llen, Dickerson, & Gabriel In traditional optimization problems, a model with exogenous data is given, and the goal is to finda feasible solution that maximizes or minimizes a given objective function. In inverse optimization (IO) problems, a particular solution or set of solutions is instead given as input, and the goal is tofind the parameters of the optimization problem that resulted in that observed solution or set ofsolutions [1, 6, 11, 12, 32, 72, 87]. Applications of inverse optimization are myriad, including thosein medical [10], energy [69], and transportation areas [86].As a more specific example of inverse optimization, we may take as input observations of strategicbehavior from rational, utility-maximizing players acting in equilibrium in a non-cooperative gameand then infer the utility functions being optimized by those players’ behaviors, as was done bythose such as Ratliff et al. [65], Waugh et al. [79, 80], and Risanger et al. [66]. More specifically,Ratliff et al. [65] observe the outputs of a multi-player Nash game in which the players can affecteach other’s objective functions in their individualized optimization problems and, then, utilizethose outputs to parameterize the players’ objective functions. In the present paper, we focus onthe context of a multi-player transportation generalized Nash game with a joint constraint in whichplayers travel across a road network with the goal of minimizing travel time, and we utilize theirobserved behavior to parameterize their underlying cost functions in the game (see Figure 1). Inparameterizing this model, we extend Ratliff et al. [65]’s framework to problems with “coupled”/jointconstraints as referenced by Rosen [67], Harker [36], and Facchinei et al. [24] which involve otherplayers’ variables in the constraints of the feasible region. This means we are parameterizing theobjective functions of players in a game in which the players can affect not only each other’sobjective functions but also each other’s feasible regions. These games with joint constraints are aparticular instance of generalized Nash equilibrium problems (GNEP) [24, 25, 29, 36, 67]. In the caseof the transportation problem, the joint constraint is a capacity constraint that requires the flowacross all players on the various arcs to stay below a certain threshold.
Fig. 1. Framework for the Transportation Game
Although we focus on this transportationproblem, we propose that our extension of theRatliff et al. [65] framework works beyond thisproblem, and it will be presented more gener-ally in Section 3. We validate its use on small butillustrative transportation networks, showingvia simulated experiments with multiple play-ers that we can recover objective function pa-rameters for the players using observed playerbehavior in GNEPs with joint constraints suchthat the same flow patterns hold under the re-covered parameterizations as under the originalcost parameterizations. We also provide somesupporting theoretical results regarding a spe-cial case of our transportation problem and the existence of unique solutions to this transportationproblem.
Work due to Keshavarz et al. [47] on parameterizing the objective function of a convex program hasbeen influential in the inverse optimization community and has been discussed and referenced inmultiple other inverse optimization papers such as [2, 6, 23, 65]. Ratliff et al. [65] apply the ideas ofthat paper to the context of game theory by making the argument that the KKT conditions [9] and llen, Dickerson, & Gabriel llen, Dickerson, & Gabriel
We contribute to the literature by extending the parameterization framework due to Ratliff et al. [65]to the context of GNEP games in which there are joint constraints among the players’ actions [24,25, 29, 36, 67]. The approach is motivated in particular by a transportation game in which there isa joint flow constraint for arcs in a road network, due to the fact that roads naturally have capacitylimitations. We differ from the transportation papers that combine equilibrium problems withinverse optimization, such as Bertsimas et al. [6], Zhang et al. [86], and the other papers listedabove involving estimating parameters, including the non-atomic routing game literature of Thaiand Bayen [73] and Bhaskar et al. [7], because our method models traffic at the level of the playerand includes this joint constraint. Our approach also differs from Kuleshov and Schrijvers [50]who model traffic at the level of the player because of this joint constraint. What’s more, whileHemmecke et al. [41] discuss theoretical results with regard to inverse optimization and GNEPs,we pursue an applied framework to recover parameterizations for cost functions. Finally, althoughone of Xu et al. [83]’s models includes a joint constraint, we focus on extending the frameworkby Ratliff et al. [65] and Keshavarz et al. [47] as opposed to their use of Ahuja and Orlin [1]’sframework. Our proposed approach can take data from simulations (as in this work) or real-worldobservations (in the case of a real-world policy-maker) involving equilibrium solutions to theseproblems, and can parameterize the objective functions of these players. The contribution to theliterature is significant because we extend a framework that has been influential in the literature(as seen in [2, 6, 23, 65]) to this new class of problems.
We examine the problem of combining 𝑁 -players’ shortest path problems to create a traffic gamein which the 𝑁 players affect each other’s costs on the network via interaction terms in theirobjective functions as well as via a joint capacity constraint in the feasible region. The problemparallels what is known as an atomic routing game [68] which features multiple players each withthe ability to affect the outcome of the traffic game and which can be represented as a multi-userNash equilibrium game [16, 61]. However, atomic routing games as presented in Roughgarden[68] do not have joint constraints. We do note that this model was also inspired by Xu et al. [83]who formulate their problem as a set of inverse shortest path problems and who did include a jointcapacity constraint in one of their formulations.For a network with 𝑛 arcs and 𝑚 nodes, let 𝑥 𝑖 ∈ R 𝑛 to be the flow decision variables for thenetwork arcs for each player 𝑖 ∈ [ 𝑁 ] , which can be fractional. The positive diagonal matrix C 𝑖 ∈ R 𝑛 × 𝑛 represents the coefficients for the interactions between player 𝑖 ’s flow and the aggregateflow, capturing the dynamic costs associated with increasing flow from the players. The vector¯ 𝑐 𝑖 ∈ R 𝑛 represents the base cost for traveling along each link in the network; this is also called thefree flow travel time by various papers [6, 15, 70, 75, 85, 86]. We also take 𝐷 ∈ {− , , } 𝑚 × 𝑛 as theconservation-of-flow matrix which keeps track of the net flow at the nodes [54] and 𝛼 ∈ R 𝑛 asthe collective capacity for the 𝑁 players’ flow on each arc. The vector 𝑓 𝑖 ∈ R 𝑚 for each player 𝑖 indicates the starting (in negative units of flow) and ending (in positive units of flow) points in thenetwork of the flow for each player 𝑖 ∈ [ 𝑁 ] . As an example for 𝑓 𝑖 , if player 𝑖 = llen, Dickerson, & Gabriel
42 and ends at node 12 in a network with 16 nodes, then their 𝑓 vector would have a dimensionof 16 𝑥
1, and we would place a − 𝑖 is:min 𝑥 𝑖 𝑥 𝑇𝑖 C 𝑖 (cid:32) 𝑁 ∑︁ 𝑗 = 𝑥 𝑗 (cid:33) + ¯ 𝑐 𝑖𝑇 𝑥 𝑖 (1a) 𝐷𝑥 𝑖 = 𝑓 𝑖 (dual variable: 𝑣 𝑖 ) (1b) 𝑥 𝑖 ≥ 𝑢 𝑖 ) (1c) 𝑁 ∑︁ 𝑗 = 𝑥 𝑗 ≤ 𝛼 (dual variable: ¯ 𝑢 ) (1d)This is a 𝑁 -player traffic game in which the 𝑁 players affect each other’s costs on the network viathe objective function, which has an interaction term between player 𝑖 ’s flow and all other players’flow and a term representing the base cost for traveling along each arc. Also, there is a sharedconstraint in the flows (1d). This constraint makes this problem a generalized Nash equilibriumproblem [24, 25, 29, 36, 67]. We also have the conservation of flow constraints in (1b) for eachplayer 𝑖 and the non-negativity constraints in (1c) for each player 𝑖 . With this notation defined, weare able to state the goal of the paper.Research Goal 2.1. The goal in this paper is to find the C 𝑖 diagonal matrix and the ¯ 𝑐 𝑖 vectorassociated with this problem for each player 𝑖 . We will examine the case in which these matrices and vectors are equal across all 𝑁 players andin which they are different across all 𝑁 players. The first case represents the situation in whicheach player faces the same interaction costs and free flow travel times as all other players which isa common assumption in aggregate models such as [6, 74, 75]. The second case represents a viewin which the players do not have a consensus view of the interaction costs or the free flow traveltimes, as proposed by [13, 14, 83].In order to obtain these matrix and vector parameters, we extend the parameterization formu-lation proposed by Ratliff et al. [65] to GNEPs with joint constraints. Ratliff et al. [65] call theirformulation a residual model after the similar Keshavarz et al. [47] model. Our extended residualmodel is defined as follows for the 𝑁 players [see, e.g., 47, 49, 65]: min C 𝑖 , ¯ 𝑐 𝑖 ,𝑣 𝑘𝑖 ,𝑢 𝑘𝑖 , ¯ 𝑢 𝑘 ∑︁ 𝑘 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) stationarity (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) 𝑁 ∑︁ 𝑖 = ||C 𝑖 (cid:169)(cid:173)(cid:171) 𝑥 𝑘𝑖 + ∑︁ 𝑗 ≠ 𝑖 𝑥 𝑘𝑗 (cid:170)(cid:174)(cid:172) + ¯ 𝑐 𝑖 + 𝐷 𝑇 𝑣 𝑘𝑖 − 𝑢 𝑘𝑖 + ¯ 𝑢 𝑘 || + complementarity (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) 𝑁 ∑︁ 𝑖 = ||( 𝑥 𝑘𝑖 ) 𝑇 𝑢 𝑘𝑖 || + complementarity (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) ||( 𝛼 − ∑︁ 𝑗 𝑥 𝑘𝑗 ) 𝑇 ( ¯ 𝑢 𝑘 )|| (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) 𝐿 ≤ 𝑑𝑖𝑎𝑔 (C 𝑖 ) ≤ 𝑈 ∀ 𝑖 (2a) 𝐿 ≤ ¯ 𝑐 𝑖 ≤ 𝑈 ∀ 𝑖 (2b)¯ 𝑢 𝑘 ≥ ∀ 𝑘 𝑢 𝑘𝑖 ≥ ∀ 𝑖, 𝑘 (2c)At a high level, our objective function naturally decomposes into three parts (labeled above as“stationarity,” “complementarity llen, Dickerson, & Gabriel • The stationarity component is a norm corresponding to the stationarity conditions of theoptimization problem for each player 𝑖 and each data piece 𝑘 , with the data pieces representingthe flow patterns we observe between the origin-destination pairs on the network. For thisnorm, the 𝑥 𝑘𝑖 are data of the flow patterns for each player 𝑖 and each origin-destination pair 𝑘 , and our variables are the diagonal of the C 𝑖 ∀ 𝑖 , the vector ¯ 𝑐 𝑖 ∀ 𝑖 , 𝑢 𝑘𝑖 as the dual variables forconstraint 1c for all 𝑖, 𝑘 , 𝑣 𝑘𝑖 as the dual variables for constraint 1b for all 𝑖, 𝑘 , and ¯ 𝑢 𝑘 as the dualvariables for constraint 1d for each data piece/origin-destination pair 𝑘 . An origin-destination(OD) pair indicates the beginning (origin) and ending (destination) points of a player’s path.We have one set of ¯ 𝑢 𝑘 dual variables for each 𝑘 that are shared across the players, which willbe further explained in Section 3. • The complementarity norm refers to the complementarity conditions for the nonnega-tivity constraint across the data and the players • The complementarity norm refers to the complementarity conditions that pair with thejoint constraint for each data piece.In conjunction, these three components of the objective allow us to minimize the error acrossthe relevant KKT conditions for all players and all data. The constraints on this objective includerequiring that the diagonal of the C 𝑖 matrices remain between two bounds 𝐿 ∈ R 𝑛 and 𝑈 ∈ R 𝑛 as captured in (2a) and that the vectors ¯ 𝑐 𝑖 remain between two bounds 𝐿 ∈ R 𝑛 and 𝑈 ∈ R 𝑛 ascaptured in (2b). The use of upper and lower bounds for the C 𝑖 and ¯ 𝑐 𝑖 was inspired by papers suchas [6, 12, 47, 49, 65] which discuss utilizing “prior information” [47] pertaining to the parameters tobe estimated to aid the mathematical program in finding viable parameterizations. We rememberthat the diagonals of C 𝑖 represent the interaction costs for each player 𝑖 and the vectors ¯ 𝑐 𝑖 representthe free flow travel cost for each player 𝑖 so, by placing bounds on these variables, their estimatedvalues will be on the same order of magnitude as the original values that generated the simulationdata. The final constraint (2c) on the objective requires that the dual variables ¯ 𝑢 𝑘 and 𝑢 𝑘𝑖 arenon-negative because they correspond with inequality constraints.One important point about this model is that it can be solved in polynomial time. From Boydand Vandenberghe [9], we know we can convert the L-1 norm components into linear constraints,with positive variables that bound the constraints from above and below and that replace the L-1terms in the objective function. This leads to a linear program, which we also know from [9] can besolved in polynomial time via the interior point method. Furthermore, the growth in the number ofvariables as a function of the inputs also grows as a polynomial function, the details of which canbe seen in Appendix A.1. Overall, we summarize all of this in the following lemma:Lemma 2.1. Problem (2) can be solved in polynomial time, and the number of variables grows as apolynomial function of the inputs.
Having established the complexity of our inverse optimization model, we can now discuss anexample of our framework.Example 2.1.
We randomly generate a parameterization for the arcs on a 4x4 grid which consistsof 48 total arcs. There are 48 arcs on this grid because there are 4 sets of 3 edges moving down thegrid and then the same pattern moving across the grid; therefore, we multiply ( )( ) by 2 and, then,multiply it by 2 again because we assume two arcs going in opposite directions per edge. We generatethis parameterization by uniformly choosing 48 numbers between 2 and 10 for the diagonal of C 𝑖 such We do want to ensure that a set of parameterizations results in a convex Nash game since we are dealing with a setof minimization problems across players, and convex games result in Nash equilibrium solutions [61, 65, 67]. Due to thefact that our minimization problems for all of the players are linear in each player’s variables, we do not require any extraconstraints [61, 65, 67], but this might be required for other applications as demonstrated by [49].llen, Dickerson, & Gabriel that C = ... = C 𝑁 = C and uniformly choosing 48 numbers between 2 and 10 for the vector ¯ 𝑐 𝑖 suchthat ¯ 𝑐 = ... = ¯ 𝑐 𝑁 = ¯ 𝑐 . We set the number of players as 𝑁 = and then generate the flows for theseparameterizations using all 240 origin-destination (OD) pairs - a number obtained by calculating (cid:0) (cid:1) ,since there are 16 nodes in a 4x4 grid, and then multiplying by 2 to account for our assumption thatthe arcs from origin 𝑎 to destination 𝑏 might be different from the arc from origin 𝑏 to destination 𝑎 .For each OD pair, all of the players each start this one unit of flow at the origin and flow this one unitof flow to the destination, meaning that for 𝑁 = players there are 10 units of flow on the network.We then input the resulting flows under these origin-destination pairs into the residual model above (2)to obtain the parameterizations for C and ¯ 𝑐 . With regard to the results of example 2.1, the original parameterizations and the parameteriza-tions obtained by the inverse optimization residual model do not align. Indeed, the Frobenius normdifference between the two C matrices is 2.4777, and the 2-norm difference between the two param-eterizations for the ¯ 𝑐 vector is 22.4281. However, if we instead use an error metric that comparesthe flows under the original parameterizations versus the flows under the IO parameterizations, wesee promising results. Indeed, when we compare the flows among the 10 players for all 240 origindestination pairs across all 48 arcs on the network, we find that the difference between the flowpatterns using a Frobenius norm metric is 8.1369e-06. We utilize the Frobenius norm because itallows us to calculate essentially a 2-norm between the two flow matrices, one resulting from theoriginal costs and one resulting from the IO costs. Consequently, as indicated by the low Frobeniusnorm error, the IO solution leads to the same type of flow pattern observed under the originalparameterization for the set of OD pair configurations considered, thus making the parametersobtained by the IO residual model a reasonable solution to the inverse optimization problem.We have now explained the framework at a basic level. Next, Section 3 discusses the connectionbetween the residual model [47, 65] and the generalized Nash equilibrium problem (GNEP) alongwith some other theoretical desiderata and preliminaries. Section 4 presents the explicit modelsused and the experimental framework. Section 5 presents our experimental results, and Section 6concludes with a brief discussion and open problems. In this section, we describe our simulation approach, form the residual model more abstractly, andthen make the connection between the simulation conditions and the residual model. In Section 3.1,we discuss the connection between the subclass of Generalized Nash Equilibrium Problems andstacked KKT conditions; we use these stacked KKT conditions as our simulation model and thusour data generation procedure. In Section 3.2, we explain the residual model utilized by Ratliff et al.[65]. Finally, in Section 3.3, we propose and prove a collorary stating that we can use the stackedKKT conditions in the residual model to recover the parameterizations for the objective functionsof the players in a GNEP.
Using the notation from Facchinei and Kanzow [25], for each player 𝑣 in a generalized Nash problem,the goal is to solve the following optimization problem, with 𝑁 as the number of players, 𝑛 𝑣 as thenumber of variables for player 𝑣 , 𝑥 𝑣 ∈ R 𝑛 𝑣 as the variables for player 𝑣 , x − 𝑣 as containing all otherplayers’ variables, 𝑋 𝑣 ( x − 𝑣 ) ⊆ R 𝑛 𝑣 as the “feasible set” parameterized by the other players’ variables, 𝜃 : R ¯ 𝑛 → R 𝑁 with ¯ 𝑛 = (cid:205) 𝑣 𝑛 𝑣 , and 𝜃 𝑣 : R ¯ 𝑛 → R (with the recognition that the last function has x − 𝑣 as parameters in its functional form once we move to the GNEP problem)min 𝑥 𝑣 𝜃 𝑣 ( 𝑥 𝑣 , x − 𝑣 ) (3a) llen, Dickerson, & Gabriel 𝑥 𝑣 ∈ 𝑋 𝑣 ( x − 𝑣 ) (3b)Facchinei and Kanzow [25] and Gabriel et al. [29] convey that a solution to this problem is an x such that no player can minimize their objective further, fixing the other players’ variables. We canform a problem to find this solution (known as a generalized Nash equilibrium) as long as we havethe Convexity Assumption stated in [25]: Convexity Assumption [25].
For every player 𝑣 and every x − 𝑣 , the objective function 𝜃 (· , x − 𝑣 ) is convex and the set 𝑋 𝑣 ( x − 𝑣 ) is closed and convex.With this Convexity Assumption, the following theorem holds in which the solution of a quasi-variational inequality QVI ( X ( x ) , F ( x )) problem is defined as x ∗ ∈ X ( x ∗ ) ⊆ R ¯ 𝑛 such that F ( x ∗ ) 𝑇 ( y − x ∗ ) ≥ ∀ y ∈ X ( x ∗ ) (4)Theorem 3.1 (Theorem 3.3 of [25]). Let a GNEP be given, satisfying the Convexity Assumption,and suppose further that the 𝜃 𝑣 functions are 𝐶 for all 𝑣 . Then a point x is a generalized Nashequilibrium if and only if it is a solution of the quasi-variational inequality QVI ( X ( x ) , F ( x )) , where X ( x ) = 𝑁 (cid:214) 𝑣 = 𝑋 𝑣 ( x − 𝑣 ) , F ( x ) = (∇ 𝑥 𝑣 𝜃 𝑣 ( x )) 𝑁𝑣 = (5)This QVI representation of the GNEP can be further simplified to a related variational inequality (VI)under certain conditions on the 𝑋 𝑣 ( x − 𝑣 ) . For a variational inequality, the goal is to find x ∗ ∈ X ⊆ R ¯ 𝑛 ,with X not parameterized by the variables, such that F ( x ∗ ) 𝑇 ( y − x ∗ ) ≥ , ∀ y ∈ X (6)One way to transform the QVI into a VI is via joint convexity of the GNEP, a property that appliesif the following definition holds pertaining to the 𝑋 𝑣 ( x − 𝑣 ) sets: Definition 3.2 (Definition 3.6 of [25]).
Let a GNEP be given, satisfying the convexity assumption.We say that this GNEP is jointly convex if for some closed, convex X ⊆ R ¯ 𝑛 and all 𝑣 = , ..., 𝑁 , wehave 𝑋 𝑣 ( x − 𝑣 ) = { 𝑥 𝑣 ∈ R 𝑛 𝑣 : ( 𝑥 𝑣 , x − 𝑣 ) ∈ X } (7)This definition means that the 𝑥 𝑣 and x − 𝑣 are all part of the same set X . This means that we cancombine the constraints from each set 𝑋 𝑣 ( x − 𝑣 ) into one set X , which also implies that we willhave one copy of any constraints known as joint constraints that depend upon all the players’variables in the set X . Furthermore, by having one copy of the joint constraints, this equates thedual variables for this constraint across the players [24, 25, 29, 36, 67]. As a result, we can utilizethe following theorem:Theorem 3.3 (Theorem 3.9 of [25]). Let a jointly convex GNEP be given with 𝐶 functions 𝜃 𝑣 .Then every solution of the 𝑉 𝐼 ( X , F ) (where X is the set in the definition of joint convexity and, as usual, F ( x ) = (∇ 𝑥 𝑣 𝜃 𝑣 ( x )) 𝑁𝑣 = ), is also a solution of the GNEP. This theorem states that we can utilize the VI with the combined X set instead of the QVI with X ( x ) to obtain a solution to the original generalized Nash problem [24–26, 29, 36, 67]. For ourtransportation problem, the X set consists of the constraints of all of the players plus one copy ofthe joint constraints: X = (cid:40) x : 𝐷𝑥 𝑣 = 𝑓 𝑣 ∀ 𝑣, 𝑥 𝑣 ≥ ∀ 𝑣, 𝑁 ∑︁ 𝑣 = 𝑥 𝑣 ≤ 𝛼 (cid:41) (8) llen, Dickerson, & Gabriel Let X be a nonempty, compact and convex subset of R ¯ 𝑛 andlet F be a continuous mapping from X into R ¯ 𝑛 . Then there exists a solution to the problem VI ( X , F ) . Due to the fact that our X is nonempty, compact, and convex as a result of the linear constraints,the lower bound provided by the non-negativity constraints, and the upper bound provided by thecapacity constraint and due to the fact that our function F is continuous, we know that a solutionexists for our VI. Furthermore, according to Harker [36] and Gabriel et al. [29], we can obtain aunique solution to this VI if the set X is compact and convex and if F defined by (6) is strictlymonotone, which is defined in the following definition: Definition 3.5 (Strict Monotonicity ([29, 37, 62])).
A function F is strictly monotone if ( F ( x ) − F ( y )) 𝑇 ( x − y ) > , ∀ x , y ∈ X , x ≠ y (9)We achieve both criteria with our traffic game, as defined and discussed in Section 2, because theequality constraints are linear, the non-joint inequality constraints are convex, the joint constraintsare convex [25, 36], the constraints create bounds on the feasible region, and the F defined in (5)composed of the gradients of the objective functions of the players is strictly monotone [29, 36].Strict monotonicity of our F defined in (5) is implied by showing it is strongly monotone [29, 37, 62],a property captured by the following definition: Definition 3.6 (Strong Monotonicity ([29, 37, 62])).
A function F is strongly monotone if for some 𝛼 > ( F ( x ) − F ( y )) 𝑇 ( x − y ) ≥ 𝛼 || x − y || , ∀ x , y ∈ X , x ≠ y (10)We recognize that a similar proposition regarding players that all have the same cost functionsand same origin-destination pairs appeared in Orda et al. [61] as described by Cominetti et al. [16],but our proof is specialized to this problem.Lemma 3.7. The F defined by (5) and created by the gradients of the objective functions defined in(1) is strongly monotone if we assume C = C = ... = C 𝑁 . See Appendix A.2 for the proof. Note though that, just because the VI form for the GNEP has onesolution (due to the strong monotonicity property), this does not mean that the original generalizedNash problem without the assumption of equal multipliers for the joint constraints would not haveother solutions (see Facchinei and Kanzow [25] for a simple example in Section 1). Furthermore, ifthe C 𝑖 matrices are not all equal, then we are not guaranteed to have a strongly monotone F .Lemma 3.8. The F defined by (5) and created by the gradients of the objective functions defined in(1) is not in general strongly monotone if we assume the C 𝑖 are not all equal. See Appendix A.3 for the counterexample. Finally, with the VI formulation, we can then use awell-known theorem from Facchinei and Pang [26] to convert the VI into a mixed complementarityproblem. Facchinei and Pang [26] define the X = { x ∈ R ˆ 𝑛 : 𝐴 x ≤ 𝑏, 𝐺 x = 𝑑 } such that 𝐴 ∈ R ˆ 𝑚 × ˆ 𝑛 , 𝐺 ∈ R ˆ 𝑙 × ˆ 𝑛 , 𝑏 ∈ R ˆ 𝑚 , 𝑑 ∈ R ˆ 𝑙 . They have the following proposition:Proposition 3.9 (Proposition 1.2.1 from [26]). Let X be defined as above. A vector x solves theVI ( X , F ) if and only if these exists 𝜆 ∈ R ˆ 𝑚 and 𝜇 ∈ R ˆ 𝑙 such that = F ( x ) + 𝐺 𝑇 𝜇 + 𝐴 𝑇 𝜆 (11a)0 = 𝑑 − 𝐺 x (11b)0 ≤ 𝜆 ⊥ 𝑏 − 𝐴 x ≥ llen, Dickerson, & Gabriel X and F , so long as the X is a set of linear inequality and linear equality constraints. Therefore,since we have linear equality and linear inequality constraints in our X , we can take the X definedin (8) and write the following system for all players 𝑖 and for a given origin destination pair encodedin 𝑓 𝑖 ∀ 𝑖 : C 𝑖 (cid:32) 𝑥 𝑖 + ∑︁ 𝑗 ≠ 𝑖 𝑥 𝑗 (cid:33) + ¯ 𝑐 𝑖 + 𝐷 𝑇 𝑣 𝑖 − 𝑢 𝑖 + ¯ 𝑢 = ∀ 𝑖 = , ..., 𝑁 (12a)0 ≤ 𝑥 𝑖 ⊥ 𝑢 𝑖 ≥ ∀ 𝑖 = , ..., 𝑁 (12b)0 ≤ 𝛼 − (cid:32)∑︁ 𝑗 𝑥 𝑗 (cid:33) ⊥ ¯ 𝑢 ≥ 𝐷𝑥 𝑖 = 𝑓 𝑖 , 𝑣 𝑖 free ∀ 𝑖 = , ..., 𝑁 (12d)with dual variables of 𝑢 𝑖 ∀ 𝑖, ¯ 𝑢, and 𝑣 𝑖 ∀ 𝑖 . We note that the system established in (12) will be solvedrepeatedly for different right-hand side 𝑓 𝑖 vectors to generate our simulation data. This will bediscussed further in Section 4. Having explained that we will use a mixed complementarity form of our generalized Nash game(1) to produce our simulation data by operating under the assumption of equal multipliers for thejoint constraint [24, 25, 29, 36, 67], we can explain the inverse optimization model we use to extractparameters from this simulation data. We utilize the residual model showcased in Ratliff et al. [65],which is quite similar to the one presented in Keshavarz et al. [47], with the difference being thatRatliff et al. [65] sum over players in a game. As a note to the reader, we will explain this modelutilizing the notation from Ratliff et al. [65]’s more recent paper Konstantakopoulos et al. [49].The model starts with the assumption that each player solves an optimization problem such asthe following: min { 𝑓 𝑖 ( 𝑥 𝑖 , 𝑥 − 𝑖 )| 𝑥 𝑖 ∈ B 𝑖 } with 𝑓 𝑖 ( 𝑥 ; 𝜃 𝑖 ) = ⟨ 𝜙 𝑖 ( 𝑥 𝑖 , 𝑥 − 𝑖 ) , 𝜁 𝑖 ⟩ + ¯ 𝑓 𝑖 ( 𝑥 ) and B 𝑖 = { 𝑥 𝑖 | ℎ 𝑖,𝑗 ( 𝑥 𝑖 ) ≤ , 𝑗 = , ..., 𝑙 𝑖 , 𝑔 𝑖,𝑞 ( 𝑥 ) = , 𝑞 = , ..., 𝑙𝑙 𝑖 } , in which 𝑖 represents the player, 𝑙 𝑖 represents the number of inequality constraints forplayer 𝑖 , 𝑙𝑙 𝑖 represents the number of equality constraints for player 𝑖 , 𝜁 𝑖 represents the linearparameterization of the objective function for player 𝑖 , and ¯ 𝑓 𝑖 is the “known” term in the objectivefunction [49]. In our game, each player also has one joint constraint inequality in the set B 𝑖 throughwhich the players can influence each others’ feasible regions. Extending this model to accommodatethis joint constraint will be covered in Section 3.3.If we take the KKT conditions of this problem [9], we obtain the following relevant “residuals”[47, 49, 65], in which 𝐷 𝑖 represents the derivative with respect to 𝑥 𝑖 : 𝑟 ( 𝑘 ) 𝑠,𝑖 ( 𝜁 𝑖 , 𝜇 𝑖 , 𝜈 𝑖 ) = 𝐷 𝑖 𝑓 𝑖 (cid:16) 𝑥 ( 𝑘 ) 𝑖 , 𝑥 ( 𝑘 )− 𝑖 (cid:17) + 𝑙 𝑖 ∑︁ 𝑗 = 𝜇 𝑗𝑖 𝐷 𝑖 ℎ 𝑖,𝑗 (cid:16) 𝑥 ( 𝑘 ) 𝑖 (cid:17) + 𝑙𝑙 𝑖 ∑︁ 𝑞 = 𝜈 𝑞𝑖 𝐷 𝑖 𝑔 𝑖,𝑞 (cid:16) 𝑥 ( 𝑘 ) 𝑖 (cid:17) (13a) 𝑟 𝑗, ( 𝑘 ) 𝑐,𝑖 ( 𝜇 𝑖 ) = 𝜇 𝑗𝑖 ℎ 𝑖,𝑗 (cid:16) 𝑥 ( 𝑘 ) 𝑖 (cid:17) , 𝑗 ∈ { , ..., 𝑙 𝑖 } (13b)with 𝑘 representing the “kth observation” in the game, 𝑠 representing the label for the stationarityresidual, and 𝑐 representing the label for the complementarity residual [49, 65]. We then have x ( 𝑘 ) data for each instance 𝑘 of the game as a composite vector containing all the players’ variablevalues, and we attempt to choose 𝜁 , 𝜇, 𝜈 such that we minimize the residuals presented above. This llen, Dickerson, & Gabriel 𝜇,𝜁,𝜈 𝑁 ∑︁ 𝑖 = 𝜂 𝑖 ∑︁ 𝑘 = 𝜒 𝑖 (cid:16) 𝑟 ( 𝑘 ) 𝑠,𝑖 ( 𝜁 , 𝜇, 𝜈 ) , 𝑟 ( 𝑘 ) 𝑐,𝑖 ( 𝜇 ) (cid:17) (14a) 𝜁 𝑖 ∈ 𝑍 𝑖 , 𝜇 𝑖 ≥ , ∀ 𝑖 ∈ { , ..., 𝑁 } (14b)with 𝑁 representing the total number of players and 𝜂 𝑖 representing the number of times a player 𝑖 engages in the game which, in our application, means the number of origin-destination runs theplayer undertakes. We assume in our implementation that 𝜂 𝑖 is equal across all players, meaningthat each player engages in each iteration of the game. It is also important to note that there areseparate parameterizations for each player 𝑖 in the form of 𝜁 𝑖 . Furthermore, these parameterizationseach belong to their own convex set 𝑍 𝑖 . Finally, we specify that 𝜒 𝑖 is a “non-negative convex penaltyfunction” [49, 65] which, according to [9], is applied to each set of residuals such that a sum ofresiduals is formed for each set of residuals. With background on the conversion of the generalized Nash problem to a mixed complementarityform and the use of the residual model for inverse optimization, we presently discuss extending theresidual inverse optimization model to accommodate the joint constraint in our generalized Nashequilibrium problem. Before we present our corollary regarding incorporating joint constraints intothe inverse optimization model, we state a theorem from Facchinei et al. [24] in modified notationfrom Facchinei and Kanzow [25].This theorem requires the following specification of X for the purposes of the theorem as: X = { x ∈ R ¯ 𝑛 : ℎ 𝑞 ( x ) ≤ 𝑞 = , ..., 𝑚, 𝑔 𝑗 ( x ) = 𝑗 = , ..., 𝑝 } (15)where 𝑚 represents in this case the total set of all inequality constraints across all players and 𝑝 represents the total set of equality constraints across all players. In the set X , we incorporate allof the constraints for all of the players, including one copy of the joint constraints. We could alsowrite the subset of constraints for each player 𝑣 as: 𝑋 𝑣 ( x − 𝑣 ) = { 𝑥 𝑣 : ℎ 𝑞 ( 𝑥 𝑣 , x − 𝑣 ) ≤ 𝑞 = , ..., 𝑚, 𝑔 𝑗 ( 𝑥, x − 𝑣 ) = 𝑗 = , ..., 𝑝 } (16)For each of these subsets 𝑋 𝑣 ( x − 𝑣 ) , all of the constraints for all of the players still exist, but we aredifferentiating between player 𝑣 ’s variables and the other players’ variables more explicitly; eachsubset 𝑋 𝑣 ( x − 𝑣 ) also has its own copy of the joint constraints. The set X and the subsets 𝑋 𝑣 ( x − 𝑣 ) ∀ 𝑣 lead to two statements of the KKT conditions. First for the 𝑋 𝑣 ( x − 𝑣 ) subsets, according to [24–26], ifwe assume that a constraint qualification holds for each player, a solution for each player 𝑣 willcorrespond with a KKT point that satisfies the following KKT conditions for each individual player 𝑣 : ∇ 𝑥 𝑣 𝜃 𝑣 ( 𝑥 𝑣 , x − 𝑣 ) + 𝑚 ∑︁ 𝑞 = 𝜆 𝑣𝑞 ∇ 𝑥 𝑣 ℎ 𝑞 ( 𝑥 𝑣 , x − 𝑣 ) + 𝑝 ∑︁ 𝑗 = 𝜈 𝑣𝑗 ∇ 𝑥 𝑣 𝑔 𝑗 ( 𝑥 𝑣 , x − 𝑣 ) = ≤ 𝜆 𝑣𝑞 ⊥ ℎ 𝑞 ( 𝑥 𝑣 , x − 𝑣 ) ≤ , ∀ 𝑞 (17b) 𝑔 𝑗 ( 𝑥 𝑣 , x − 𝑣 ) = , 𝜈 𝑣𝑗 free , ∀ 𝑗 (17c)Second for the X set, we know from [5, 24, 26] that, if a constraint qualification is satisfied, then theKKT conditions are satisfied by a solution to the VI (assuming the presence of relevant multipliers)with F defined as (5) and X defined as (15). We also know from [26] that, if ℎ 𝑞 are convex and 𝑔 𝑗 llen, Dickerson, & Gabriel F defined as (5) and X defined as (15). Therefore, we have the resulting KKT conditions from the VI as: F ( x ) + 𝑚 ∑︁ 𝑞 = 𝜆 𝑞 ∇ x ℎ 𝑞 ( x ) + 𝑝 ∑︁ 𝑗 = 𝜈 𝑗 ∇ x 𝑔 𝑗 ( x ) = ≤ 𝜆 𝑞 ⊥ ℎ 𝑞 ( x ) ≤ , ∀ 𝑞 (18b) 𝑔 𝑗 ( x ) = , 𝜈 𝑗 free , ∀ 𝑗 (18c)We can now convey Theorem 4.8 from [25] who cite [24, 36] which states that we can move betweenthe two forms of the KKT conditions (17) and (18) so long as we assume the multipliers 𝜆 𝑣 , 𝜈 𝑣 for (17)are equal across all players 𝑣 . This statement of the theorem is slightly different from either [24]’sstatement or [25]’s statement because we have added equality constraints and their associatedmultipliers into our formulation of these KKT conditions.Theorem 3.10 (Theorem 4.8 [25]). Consider the jointly convex GNEP with ℎ 𝑞 ∀ 𝑞, 𝑔 𝑗 ∀ 𝑗, 𝜃 𝑣 being 𝐶 . Then the following statements hold: • Let ¯ x be a solution of the VI ( X , F ) (with X defined as (15), F defined as (5), and VI defined as (6))such that the KKT conditions (18) hold with some multipliers ¯ 𝜆 and ¯ 𝜈 . Then ¯ x is a solution of theGNEP, and the corresponding KKT conditions (17) are satisfied with 𝜆 = ... = 𝜆 𝑁 = ¯ 𝜆 and with 𝜈 = 𝜈 = ... = 𝜈 𝑁 = ¯ 𝜈 . • Conversely, assume that ¯ x is a solution of the GNEP such that the KKT conditions (17) aresatisfied with 𝜆 = ... = 𝜆 𝑁 and with 𝜈 = 𝜈 = ... = 𝜈 𝑁 . Then ( ¯ x , ¯ 𝜆, ¯ 𝜈 ) with ¯ 𝜆 = 𝜆 and with ¯ 𝜈 = 𝜈 is a KKT point of VI ( X , F ) , and ¯ x itself is a solution of VI ( X , F ) . Proof. As a brief proof of this theorem (which is based off of, and is a relatively direct extensionof, a related proof by Facchinei et al. [24]), for the first statement, we notice a correspondencebetween (17a) concatenated for all players 𝑣 and (18a) as long as the multipliers are equal acrossplayers and are equal to the ¯ 𝜆, ¯ 𝜈 multipliers. This correspondence continues for (17b) & (18b) and(17c) & (18c) so long as the multipliers are equal. Therefore, the ( ¯ x , ¯ 𝜆, ¯ 𝜈 ) from (18) is a KKT point forthe conditions (17) for each player. According to [24], as long as the conditions (17) are sufficientfor optimality for each player, which is true if we assume F is convex, ℎ 𝑞 are convex, and 𝑔 𝑗 areaffine [9], ¯ x with multipliers ¯ 𝜆, ¯ 𝜈 is a solution for each player’s optimization problem, thus makingit a GNEP solution.For the second statement, we point to the fact that, now that we are given that the multipliersare equal across all of the players, then there is a correspondence between (17) (concatenated for allplayers 𝑣 for the stationarity constraint) and (18) as long as ¯ 𝜆 = 𝜆 = ... = 𝜆 𝑁 and ¯ 𝜈 = 𝜈 = ... = 𝜈 𝑁 .Consequently, we have a KKT point for the (18) conditions that comes from the given solution to(17). We know from [26] that if there is a KKT point for the conditions (18) and we assume theconstraints ℎ 𝑞 are convex and 𝑔 𝑗 are affine, then the ¯ x associated with the KKT point is a solutionfor the VI. □ Our technique, as discussed in Section 2, can be looked at as a corollary to this theorem, where westate that we can use the stacked KKT conditions with one copy of the joint/shared constraints toparameterize the cost functions of the players.Corollary 3.11.
With the assumption of equal multipliers for the joint constraint(s) in a GNEP,the stationarity conditions and complementarity conditions of the VI KKT conditions (18) may be usedto form a residual model of the form seen in [47, 49, 65], with 𝜒 as a “non-negative convex penality llen, Dickerson, & Gabriel function” as in those papers and in [9] and with 𝑘 representing the multiple data points x 𝑘 : min 𝜁 𝑖 ,𝜆 𝑞 ,𝜈 𝑗 ∑︁ 𝑘 𝜒 (cid:32) F ( x 𝑘 ) + 𝑚 ∑︁ 𝑞 = 𝜆 𝑞 ∇ x ℎ 𝑞 ( x 𝑘 ) + 𝑝 ∑︁ 𝑗 = 𝜈 𝑗 ∇ x 𝑔 𝑗 ( x 𝑘 ) (cid:33) + ∑︁ 𝑘,𝑞 𝜒 (cid:16) 𝜆 𝑞 ℎ 𝑞 ( x 𝑘 ) (cid:17) (19a) 𝜆 𝑞 ≥ ∀ 𝑞, 𝜁 𝑖 ∈ 𝑍 𝑖 ∀ 𝑖 (19b)Proof. We know from Theorem 3.10 that there is a correspondence between the solutions to(17) for all players 𝜈 and the solution to (18) so long as the multipliers are equal across these twosets of KKT conditions. Therefore, we can take the stationarity condition and the complementarityconditions of (18) as the residuals for the GNEP problem and form the same kind of residual problemas (14). □ As a result of corollary 3.11, we can find parameterizations for players’ cost functions in generalizedNash games with joint constraints. Specifically, in our transportation problem, we can find C 𝑖 and ¯ 𝑐 𝑖 for players 𝑖 ∈ [ 𝑁 ] . Next, in Section 4, we explain our experimental framework, walking the readerthrough our simulation set-up and the specifications for our residual model for the transportationproblem. To test and explore our inverse optimization framework for generalized Nash Equilibrium games,we utilize various grid networks, which can be expressed as a certain number of nodes verticallyand a certain number of nodes horizontally (2 nodes by 2 nodes, 3 nodes by 3 nodes, 4 nodesby 4 nodes, and 5 nodes by 5 nodes), and the Sioux Falls network [51, 76] for the traffic problempresented in the introduction; examples of these networks are visualized in Figures 2a and 2b,respectively. (a) 4 nodes x 4 nodes Grid (b) Sioux Falls Network [51, 76]Fig. 2. Networks utilized for testing inverse optimization framework
To generate the solutions to use in our inverse optimization framework, we utilize a slightlyvaried form of the KKT conditions listed above in (12) in which we substitute the stationarityconstraint (12a) into the nonnegative complementarity constraints (12b) thus producing (20a)[9, 25, 29]: llen, Dickerson, & Gabriel ≤ 𝑥 𝑖 ⊥ C 𝑖 (cid:32) 𝑥 𝑖 + ∑︁ 𝑗 ≠ 𝑖 𝑥 𝑗 (cid:33) + ¯ 𝑐 𝑖 + 𝐷 𝑇 𝑣 𝑖 + ¯ 𝑢 ≥ ∀ 𝑖 = , ..., 𝑁 (20a)0 ≤ 𝛼 − (cid:32)∑︁ 𝑗 𝑥 𝑗 (cid:33) ⊥ ¯ 𝑢 ≥ 𝐷𝑥 𝑖 = 𝑓 𝑖 , 𝑣 𝑖 free ∀ 𝑖 = , ..., 𝑁 (20c)We utilize the implementation of PATH [18, 27] from the software GAMS [30, 31] for this problem.As stated in the introduction, our goal is to solve the inverse optimization residual model forpoints generated from the GNE problem, specifically with the purpose of finding the parameters C 𝑖 and ¯ 𝑐 𝑖 for each player 𝑖 . In some of the experiments, these C 𝑖 and ¯ 𝑐 𝑖 will be assumed the sameacross all players and, in others, they will be assumed different across all players. As a reminder ofthe residual model original presented in Section 2, we present it again here [47, 49, 65]: min C 𝑖 , ¯ 𝑐 𝑖 ,𝑣 𝑘𝑖 ,𝑢 𝑘𝑖 , ¯ 𝑢 𝑘 ∑︁ 𝑘 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) stationarity (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) 𝑁 ∑︁ 𝑖 = ||C 𝑖 (cid:169)(cid:173)(cid:171) 𝑥 𝑘𝑖 + ∑︁ 𝑗 ≠ 𝑖 𝑥 𝑘𝑗 (cid:170)(cid:174)(cid:172) + ¯ 𝑐 𝑖 + 𝐷 𝑇 𝑣 𝑘𝑖 − 𝑢 𝑘𝑖 + ¯ 𝑢 𝑘 || + complementarity (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) 𝑁 ∑︁ 𝑖 = ||( 𝑥 𝑘𝑖 ) 𝑇 𝑢 𝑘𝑖 || + complementarity (cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123) ||( 𝛼 − ∑︁ 𝑗 𝑥 𝑘𝑗 ) 𝑇 ( ¯ 𝑢 𝑘 )|| (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) 𝐿 ≤ 𝑑𝑖𝑎𝑔 (C 𝑖 ) ≤ 𝑈 ∀ 𝑖 (21a) 𝐿 ≤ ¯ 𝑐 𝑖 ≤ 𝑈 ∀ 𝑖 (21b)¯ 𝑢 𝑘 ≥ ∀ 𝑘 𝑢 𝑘𝑖 ≥ ∀ 𝑖, 𝑘 (21c)In this mathematical program, we minimize the error across the relevant KKT conditions for allplayers and all data such that the diagonal of the C 𝑖 matrices and the ¯ 𝑐 𝑖 vectors remain betweenthe two bounds. The 𝐿 , 𝐿 ∈ R 𝑛 and 𝑈 , 𝑈 ∈ R 𝑛 bounds are determined by the experimentationbounds for those variables; for example, if we randomly choose the diagonal of C 𝑖 to uniformlycome from between 1 and 5, then the upper and lower bounds for the diagonal will be 1 and 5. Weuse the 1-norm to preserve the convexity of the problem and to create linear terms as discussed inSection 2.For our experiments, we follow a general framework in which various elements can be modified.The general framework is as follows: ALGORITHM 1:
Basic Experimental Framework
Input: 𝐷 conservation-of-flow matrix, 𝑛 nodes, 𝑚 arcs, 𝛼 capacity level, 𝑁 number of players, bounds 𝐿 , 𝐿 , 𝑈 , 𝑈 Exhaustively generate P , the set of all unique origin-destination pairs for a given network for each pair in P do Solve the mixed complementarity problem (20) form of GNEP (1) for RHS 𝑓 𝑖 ∀ 𝑖 constructed from thecurrent pair, with one unit of flow for each playerStore the 𝑁 player flow vectors 𝑥 𝑖 𝑖 = , ..., 𝑁 (which are the data points) end Input the sets of data points into the IO mathematical program (21) to obtain estimates for C 𝑖 and ¯ 𝑐 𝑖 ∀ 𝑖 We again note that, in some experiments, we assume C 𝑖 and ¯ 𝑐 𝑖 are the same across all 𝑖 and, inothers, we assume C 𝑖 and ¯ 𝑐 𝑖 are different across all 𝑖 . We also utilize all P origin-destination pairsfor a given network because we want to gain a detailed picture of flow on that network. However,we recognize that there are many more configurations of possible flow from which we could have llen, Dickerson, & Gabriel In this section, we present the numerical results of experiments with the proposed inverse opti-mization framework to parameterize the GNEP (1) using the two types of networks showcased inSection 4. To implement the experimental setup described in Section 4, we used a data generationand optimization pipeline including state of the art solvers from PATH [18, 27] in GAMS and Gurobi[34]. This included randomly generating from a uniform distribution the original costs for the C 𝑖 and ¯ 𝑐 𝑖 parameters in the simulation model outlined in (20). See Appendix A.4 for more informationabout the experimental setup.The boxplots in this section denote two different measures of error. First, there is the value of theobjective function of the residual model, denoting how closely the chosen parameterizations causethe input data to satisfy the KKT conditions [47, 49, 65]. Second, there is the total flow error metric,which measures the difference between the flows under both the original parameterization and theIO residual model parameterization using the same set of OD pairs P . The Frobenius-norm metricis utilized, which acts as a vector 2-norm for matrices. This is useful compared to other matrixnorms because it calculates a total difference between the two matrices, not a maximal differencealong rows or columns or an abstract eigenvalue metric as in the case of the more traditional matrixnorms [81]. This leads to the following metric. Definition 5.1 (Flow error).
The flow error is calculated by first taking the squared differencebetween the flow under the original parameterization and the flow under the IO parameterizationfor all origin-destination pairings, all players, and all arcs. Then, a sum is taken across thesesquared differences and, finally, the square root of this summation is taken. This error representsthe Frobenius norm between the two sets of flows.This error is also normalized by the total number of arc flows, which is calculated by multiplyingthe number of origin-destination pairings, the number of players, and the number of arcs. It formsthe following normalized metric.
Definition 5.2 (Normalized flow error).
The normalized flow error is calculated by dividing theflow error by a factor created by multiplying the number of origin-destination pairings, the numberof players, and the number of arcs. It represents a per unit level of error.Note that data are collected regarding the differences between the original parameterizationsand the IO parameterizations, and it is found that there are differences between the two. However,the important metrics are the objective function value and the flow error metrics because the firstmeasures how well the IO parameterization fits the model and the model data and because thesecond showcases whether or not the IO parameterization leads to the same flow patterns observedfor the OD pairs under the original parameterization. Note also that the objective function valuesfor the experiments can be found in Appendix A.6 and timing information for the experiments canbe found in Appendix A.7.
Overall, the results from the four experimental groups (with two in the same cost randomized costscategory and two in the different cost randomized costs category) are encouraging. The maximumobjective function values for all the experiments across the groups are on the order of 1e-6, and For all of the box plots, the “whiskers” are placed at quartile 1 - 1.5 (quartile 3 - quartile 1) and at quartile 3 + 1.5(quartile 3 - quartile 1), and the “dots” are outliers [43].llen, Dickerson, & Gabriel F functions as defined by (5) across the number of arcs present in the grids (2x2-5x5) and in SiouxFalls as well as the number of players, with the exception of one trial in the 5x5 grid and 10 playerscase, are strongly monotone functions. As indicated in Lemma 3.8, we cannot make this guaranteein general for players with different interaction costs. This matters because the different startingpoint scheme for solving the simulation model is bolstered by the positive definiteness assumptionand, as can be seen in Appendix A.5, there is less of a chance that the matrix will be positive definiteas the number of players increases. Therefore, this suggests further computational work that couldbe done to see if non-positive definite matrices of our type work in general in our simulationframework. It should also be noted that, for different costs, the experiments with 𝛼 = 𝑁 = 𝛼 values of ( . ) 𝑁 validatesthe extension of the framework, because some of the OD pairs involved start, end, or both nodeswhere only two arcs were coming out of or into the node, which meant the capacity constraintswere guaranteed to be tight. Table 1 details each of the experiment groups for this subsection and provides an experiment groupnumber to which we will refer in this section. In this subsection, we assume that C 𝑖 and ¯ 𝑐 𝑖 are equalacross all players 𝑖 such that C = C = ... = C 𝑁 = C and ¯ 𝑐 = ¯ 𝑐 = ... = ¯ 𝑐 𝑁 = ¯ 𝑐 . 𝑁 ) C Bounds ¯c Bounds 𝛼 . ( 𝑁 ) , 𝑁 . ( 𝑁 ) , 𝑁 Table 1. Experimental details for the “Same Costs” setting. Here, C Bounds indicating the range for random-ization and the bounds used in the IO mathematical program for the diagonal of the C , ¯ 𝑐 Bounds indicatesthe same for the ¯ 𝑐 parameters, and 𝛼 refers to that parameter value Experiment group 1 (Grid) iterates over generated 2x2, 3x3, 4x4, and 5x5 grids, each with 10trials of randomly chosen parameters, which the inverse optimization residual model attempts toestimate. For each grid, three different player numbers ( 𝑁 = , ,
10) and two different 𝛼 values,one set to half the number of players for each 𝑁 and one set to exactly the number of players foreach 𝑁 , are considered. The graphs for this experiment group display many box plots, and thesystem utilized for labeling the box plots is Grid Size/Number of Players/Alpha Value. Each setof eight moving from left to right indicate the same number of players. With regard to flow error, Note that the one non-positive definite matrix in the 5x5 grid with 10 players and 𝛼 =
10 example appears to work inthe simulation framework because the flow error and objective function values are both low.llen, Dickerson, & Gabriel 𝛼 values, as grid size increases, thetotal flow error also increases. It is important to note that everything on this graph is still on theorder of 1e-7, but this trend is also evident. It likely results from the fact that, as grid size increases,the number of arcs and number of OD pairs also increases and, since this measure is calculatedacross all OD pairs, all players, and all arcs, then even if consistent error were assumed across allarcs, the norm would have to increase. Indeed, upon examining the accompanying Figure 3b, it isapparent that the normalized flow error, calculated as described in Definition 5.2, decreases as gridsize increase across all of the subsets of players and 𝛼 values. (a) Flow Error for Experiment Group 1(b) Normalized Flow Error for Experiment Group 1Fig. 3. Flow Error Metrics for Experiment Group 1: The labeling at the bottom of the graphs indicates attributesof the boxplot, specifically the Grid Size/Number of Players/Alpha Value Experiment group 2 (Sioux Falls) again consists of 10 trials for each of 𝑁 = , ,
10 players andthe two 𝛼 = ( . ) 𝑁 , 𝑁 values. The labeling for the graphics is set to Number of Players/AlphaValue. In examining flow errors, Figure 4a showcases that the flow errors do not appear to follow aset pattern when moving between the different numbers of players, yet they are all still small andon the order of 1e-7 or lower. However, in Figure 4b, the median flow error decreases as number ofplayers increases, even on this very small scale (1e-12). llen, Dickerson, & Gabriel (a) Flow Error for Experiment Group 2 (b) Normalized Flow Error for Experiment Group 2Fig. 4. Flow Error Metrics for Experiment Group 2: The labeling at the bottom of the graphs indicates attributesof the boxplot, specifically the Number of Players/Alpha Value Next, we assume that the costs are not the same across all of the players, meaning that a new 𝐶 𝑖 anda new ¯ 𝑐 𝑖 are drawn for each player 𝑖 . Table 2 describes each experiment group in this subsection; asin the previous subsection, an experiment group number is assigned to the two sets of experimentgroups. 𝑁 ) C i Bounds ¯c i Bounds 𝛼 . ( 𝑁 ) , 𝑁 . ( 𝑁 ) , 𝑁 Table 2. Tables with Experiment Details for Different Costs. Note that C 𝑖 Bounds indicating the rangefor randomization and the bounds used in the IO mathematical program for the diagonal of the C 𝑖 , ¯ 𝑐 𝑖 Boundsindicates the same for the ¯ 𝑐 𝑖 parameters, and 𝛼 refers to that parameter value Experiment group 3 is a mirror image of experiment group 1, except that all of the players do nothave the same costs. One issue with this experiment group was that not all of the trials finishedfor the case of the 5x5 grid, 10 players, and 𝛼 =
5; only 6 of the trials finished in under 24 hours, so they are the ones included in the box plots. Similar to the flow error for the grids under thesame costs (Figure 3a), Figure 5a demonstrates that the median flow error increases as the grid sizeincreases for all the subsets of player number and 𝛼 value. Figure 5b illustrates much of the samedecrease in median normalized flow as grid size increases (among the subsets) as in Figure 3b forsame costs across all players.Experiment group 4 (Sioux Falls) again consists of 10 trials for each of 𝑁 = , ,
10 players andthe two 𝛼 = ( . ) 𝑁 , 𝑁 values. The labeling for the graphics is Number of Players/Alpha Value.It should be noted that the trials did not finish at the time of this posting for the 𝑁 =
10 and 𝛼 = 𝑁 =
10 and 𝛼 = This means about 24 hours were given for each trial before the trial was stopped, with the exception of one that wasstopped by the computer before convergence.llen, Dickerson, & Gabriel (a) Flow Error for Experiment Group 3 (b) Normalized Flow Error for Experiment Group 3Fig. 5. Flow Error Metrics for Experiment Group 3: The labeling at the bottom of the graphs indicates attributesof the boxplot, specifically the Number of Players/Alpha Value. Note: Only 6 trials are included for 5/10/5.0,see note in the text.(a) Flow Error for Experiment Group 4 (b) Normalized Flow Error for Experiment Group 4Fig. 6. Flow Error Metrics for Experiment Group 4: The labeling at the bottom of the graphs indicates attributesof the boxplot, specifically the Number of Players/Alpha Value. Note: The trials did not finish in time for 𝑁 = and 𝛼 = , hence that boxplot is not included. Ratliff et al. [65] applied the framework due to Keshavarz et al. [47] to multi-player Nash games.We extended this framework to the case of generalized Nash games with joint/shared constraintsby proving that we can use the VI KKT conditions in conjunction with a residual model to recover aparameterization for the objective functions of the players in the multi-player game. We have seenthat, although our model may not recover the original parameterization, it recovers a parameteriza-tion that produces the same flow patterns as the original parameterization. This holds true acrossmultiple grid sizes, the Sioux Falls network, different assumptions regarding players’ perceivedcosts, and the majority of restrictive 𝛼 capacity settings and the associated numbers of players.Further research could extend our model to the setting where the multipliers for the sharedconstraint are not assumed to be equal, especially in the case when the cost functions differ betweenthe players [29]. For more on this, we suggest to readers the work of Nabetani et al. [57]. In addition,we could extend the research to a real-world traffic situation with real data, in which we wouldlikely need to incorporate more players into our framework and different origin-destination trafficpatterns, including ones in which not all the players leave from the same origin node and head tothe same destination node. Furthermore, we could attempt to address some of the issues with theNash approach, including looking at a low regret condition on the players as in [58, 79, 80]. Wecould also explore making our model more robust to the issues indicated by [63] such as rationalityassumptions and incorrect models. llen, Dickerson, & Gabriel ACKNOWLEDGMENTS
Allen was partially funded by a Graduate Fellowship in STEM Diversity while completing thisresearch (formerly known as a National Physical Science Consortium Fellowship). Allen was alsosupported by a Flagship Fellowship and a Dean’s Fellowship from the University of Maryland-College Park and has worked for Johns Hopkins University Applied Physics Lab in the summers.Dickerson was supported in part by NSF CAREER Award IIS-1846237, NSF D-ISN Award
REFERENCES [1] Ravindra K Ahuja and James B Orlin. 2001. Inverse optimization.
Operations Research
49, 5 (2001), 771–783.[2] Anil Aswani, Zuo-Jun Shen, and Auyon Siddiq. 2018. Inverse optimization with noisy data.
Operations Research
66, 3(2018), 870–892.[3] Andreas Bärmann, Alexander Martin, Sebastian Pokutta, and Oskar Schneider. 2018. An online-learning approach toinverse optimization. arXiv preprint arXiv:1810.12997 (2018).[4] Andreas Bärmann, Sebastian Pokutta, and Oskar Schneider. 2017. Emulating the expert: inverse optimization throughonline learning. In
International Conference on Machine Learning . 400–410.[5] Mokhtar S Bazaraa, Hanif D Sherali, and Chitharanjan M Shetty. 2013.
Nonlinear programming: theory and algorithms .John Wiley & Sons.[6] Dimitris Bertsimas, Vishal Gupta, and Ioannis Ch Paschalidis. 2015. Data-driven estimation in equilibrium usinginverse optimization.
Mathematical Programming
Games and Economic Behavior
118 (2019), 533–569.[8] John R Birge and Derek F Holmes. 1992. Efficient solution of two-stage stochastic linear programs using interior pointmethods.
Computational Optimization and Applications
1, 3 (1992), 245–276.[9] Stephen Boyd and Lieven Vandenberghe. 2004.
Convex optimization . Cambridge university press.[10] Timothy CY Chan, Tim Craig, Taewoo Lee, and Michael B Sharpe. 2014. Generalized inverse multiobjective optimizationwith application to cancer therapy.
Operations Research
62, 3 (2014), 680–695.[11] Timothy CY Chan and Neal Kaw. 2020. Inverse optimization for the recovery of constraint parameters.
EuropeanJournal of Operational Research
Management Science
65, 3 (2019), 1115–1135.[13] Joseph YJ Chow and Shadi Djavadian. 2015. Activity-based market equilibrium for capacitated multimodal transportsystems.
Transportation Research Procedia
Transportation Research Part B: Methodological
46, 3 (2012), 463–479.[15] Joseph YJ Chow, Stephen G Ritchie, and Kyungsoo Jeong. 2014. Nonlinear inverse optimization for parameter estimationof commodity-vehicle-decoupled freight assignment.
Transportation Research Part E: Logistics and TransportationReview
67 (2014), 71–91.[16] Roberto Cominetti, José R Correa, and Nicolás E Stier-Moses. 2006. Network games with atomic players. In
InternationalColloquium on Automata, Languages, and Programming . Springer, 525–536.[17] Richard W. Cottle, Jong-Shi Pang, and Richard E. Stone. 1992.
The Linear Complementarity Problem . Academic Press,Inc.[18] Steven P Dirkse and Michael C Ferris. 1995. The path solver: a nommonotone stabilization scheme for mixedcomplementarity problems.
Optimization Methods and Software
5, 2 (1995), 123–156.[19] Chaosheng Dong, Yiran Chen, and Bo Zeng. 2018. Generalized Inverse Optimization through Online Learning. In
Advances in Neural Information Processing Systems . 86–95.[20] CW Duin and A Volgenant. 2006. Some inverse optimization problems under the Hamming distance.
European Journalof operational research [21] B Curtis Eaves. 1971. On the basic theorem of complementarity. Mathematical Programming
1, 1 (1971), 68–75.[22] Rune Elvik. 2014. A review of game-theoretic models of road user behaviour.
Accident Analysis & Prevention
62 (2014),388–396.[23] Peyman Mohajerin Esfahani, Soroosh Shafieezadeh-Abadeh, Grani A Hanasusanto, and Daniel Kuhn. 2018. Data-driveninverse optimization with imperfect information.
Mathematical Programming
Operations Research Letters
35, 2 (2007), 159–164.[25] Francisco Facchinei and Christian Kanzow. 2010. Generalized Nash equilibrium problems.
Annals of OperationsResearch
Finite-dimensional variational inequalities and complementarity problems .Springer Science & Business Media.[27] Michael C Ferris and Todd S Munson. 2000. Complementarity problems in GAMS and the PATH solver.
Journal ofEconomic Dynamics and Control
Complementaritymodeling in energy markets
EuropeanJournal of Operational Research
Proceedings of the 7th Python in Science Conference , Gaël Varoquaux, Travis Vaught, and JarrodMillman (Eds.). Pasadena, CA USA, 11 – 15.[36] Patrick T Harker. 1991. Generalized Nash games and quasi-variational inequalities.
European journal of Operationalresearch
54, 1 (1991), 81–94.[37] Patrick T Harker and Jong-Shi Pang. 1990. Finite-dimensional variational inequality and nonlinear complementarityproblems: a survey of theory, algorithms and applications.
Mathematical programming
48, 1-3 (1990), 161–220.[38] William E Hart, Carl D Laird, Jean-Paul Watson, David L Woodruff, Gabriel A Hackebeil, Bethany L Nicholson, andJohn D Siirola. 2017.
Pyomo—Optimization Modeling in Python: Second Edition . Vol. 67. Springer Optimization and ItsApplications.[39] William E Hart, Jean-Paul Watson, and David L Woodruff. 2011. Pyomo: modeling and solving mathematical programsin Python.
Mathematical Programming Computation
3, 3 (2011), 219.[40] Philip Hartman and Guido Stampacchia. 1966. On some non-linear elliptic differential-functional equations.
Actamathematica arXiv preprint arXiv:0903.4577 (2009).[42] R.A. Horn and C.R. Johnson. 1985.
Matrix Analysis . Cambridge University Press.[43] John Hunter, Darren Dale, Eric Firing, Michael Droettboom, and Matplotlib development team. April 2020. mat-plotlib.pyplot.boxplot. https://matplotlib.org/3.2.1/api/_as_gen/matplotlib.pyplot.boxplot.html.[44] J. D. Hunter. 2007. Matplotlib: A 2D graphics environment.
Computing in Science & Engineering
9, 3 (2007), 90–95.https://doi.org/10.1109/MCSE.2007.55[45] Sanjay Jain and Nitin Arya. 2013. An inverse capacitated transportation problem.
IOSR Journal of Mathematics
5, 4(2013), 24–27.[46] CHARLES R Johnson. 1970. Positive definite matrices.
The American Mathematical Monthly
77, 3 (1970), 259–264.[47] Arezou Keshavarz, Yang Wang, and Stephen Boyd. 2011. Imputing a convex objective function. In . IEEE, 613–619.[48] Ioannis C Konstantakopoulos, Andrew R Barkan, Shiying He, Tanya Veeravalli, Huihan Liu, and Costas Spanos. 2019.A deep learning and gamification approach to improving human-building interaction and energy efficiency in smartinfrastructure.
Applied energy
237 (2019), 810–821.[49] Ioannis C Konstantakopoulos, Lillian J Ratliff, Ming Jin, S Shankar Sastry, and Costas J Spanos. 2017. A robust utilitylearning framework via inverse optimization.
IEEE Transactions on Control Systems Technology
26, 3 (2017), 954–970.[50] Volodymyr Kuleshov and Okke Schrijvers. 2015. Inverse game theory: Learning utilities in succinct games. In
International Conference on Web and Internet Economics . Springer, 413–427.llen, Dickerson, & Gabriel [51] Larry J LeBlanc, Edward K Morlok, and William P Pierskalla. 1975. An efficient approach to solving the road networkequilibrium traffic assignment problem. Transportation research
Operations Research
39, 5 (1991), 757–770.[54] Patrice Marcotte and Michael Patriksson. 2007. Traffic equilibrium.
Handbooks in Operations Research and ManagementScience
14 (2007), 623–713.[55] MATLAB. 2020.
Version 9.8.0.1323502 (R2020a) . The MathWorks, Inc., Natick, Massachusetts.[56] Wes McKinney et al. 2010. Data structures for statistical computing in python. In
Proceedings of the 9th Python inScience Conference , Vol. 445. Austin, TX, 51–56.[57] Koichi Nabetani, Paul Tseng, and Masao Fukushima. 2011. Parametrized variational inequality approaches to generalizedNash equilibrium problems with shared constraints.
Computational Optimization and Applications
48, 3 (2011), 423–452.[58] Denis Nekipelov, Vasilis Syrgkanis, and Eva Tardos. 2015. Econometrics for learning agents. In
Proceedings of theSixteenth ACM Conference on Economics and Computation . 1–18.[59] Thai Dung Nguyen. 2010.
Application of robust and inverse optimization in transportation . Ph.D. Dissertation. Mas-sachusetts Institute of Technology.[60] Travis E. Oliphant. 2006.
A guide to NumPy . USA: Trelgol Publishing.[61] Ariel Orda, Raphael Rom, and Nahum Shimkin. 1993. Competitive routing in multiuser communication networks.
IEEE/ACM Transactions on networking
1, 5 (1993), 510–521.[62] James M Ortega and Werner C Rheinboldt. 2000.
Iterative solution of nonlinear equations in several variables . SIAM.[63] Alexander Peysakhovich, Christian Kroer, and Adam Lerer. 2019. Robust multi-agent counterfactual prediction. In
Advances in Neural Information Processing Systems . 3077–3087.[64] Lillian J Ratliff, Samuel A Burden, and S Shankar Sastry. 2013. Characterization and computation of local Nashequilibria in continuous games. In . IEEE, 917–924.[65] Lillian J Ratliff, Ming Jin, Ioannis C Konstantakopoulos, Costas Spanos, and S Shankar Sastry. 2014. Social game forbuilding energy efficiency: Incentive design. In . IEEE, 1011–1018.[66] Simon Risanger, Stein-Erik Fleten, and Steven A Gabriel. 2020. Inverse Equilibrium Analysis of Oligopolistic ElectricityMarkets.
IEEE Transactions on Power Systems
35, 6 (2020), 4159–4166.[67] J Ben Rosen. 1965. Existence and uniqueness of equilibrium points for concave n-person games.
Econometrica: Journalof the Econometric Society (1965), 520–534.[68] Tim Roughgarden. [n.d.]. Routing Games. In
Algorithmic Game Theory . Cambridge University Press.[69] Javier Saez-Gallego and Juan M Morales. 2017. Short-term forecasting of price-responsive loads using inverseoptimization.
IEEE Transactions on Smart Grid
9, 5 (2017), 4805–4814.[70] Enrico Siri, Silvia Siri, and Simona Sacone. 2020. A progressive traffic assignment procedure on networks affected bydisruptive events. In . IEEE, 130–135.[71] Oliver Stein and Nathan Sudermann-Merx. 2018. The noncooperative transportation problem and linear generalizedNash games.
European Journal of Operational Research arXivpreprint arXiv:2006.08923 (2020).[73] Jérôme Thai and Alexandre M Bayen. 2017. Learnability of edge cost functions in routing games. In . IEEE, 6422–6429.[74] Jérôme Thai and Alexandre M Bayen. 2018. Imputing a variational inequality function or a convex objective function:A robust approach.
J. Math. Anal. Appl. . IEEE, 689–695.[76] Transportation Networks for Research Core Team. June 2020. Transportation Networks for Research. https://github.com/bstabler/TransportationNetworks.[77] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski,Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. JarrodMillman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, İlhan Polat, Yu Feng,Eric W. Moore, Jake Vand erPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero,Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1. 0Contributors. 2020. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.
Nature Methods
17 (2020),261–272. https://doi.org/10.1038/s41592-019-0686-2llen, Dickerson, & Gabriel [78] Stéfan van der Walt, S Chris Colbert, and Gael Varoquaux. 2011. The NumPy array: a structure for efficient numericalcomputation. Computing in Science & Engineering
13, 2 (2011), 22–30.[79] Kevin Waugh, Brian D Ziebart, and J Andrew Bagnell. 2011. Computational Rationalization: The Inverse EquilibriumProblem. (2011).[80] Kevin Waugh, Brian D Ziebart, and J Andrew Bagnell. 2013. Computational rationalization: The inverse equilibriumproblem. arXiv preprint arXiv:1308.3506 (2013).[81] Eric W. Weisstein. From MathWorld–A Wolfram Web Resource. Matrix Norm.https://mathworld.wolfram.com/MatrixNorm.html.[82] Eric W. Weisstein. From MathWorld–A Wolfram Web Resource. Positive Definite Matrix.https://mathworld.wolfram.com/PositiveDefiniteMatrix.html.[83] Susan Jia Xu, Mehdi Nourinejad, Xuebo Lai, and Joseph YJ Chow. 2018. Network learning via multiagent inversetransportation problems.
Transportation Science
52, 6 (2018), 1347–1364.[84] He Zhang, Yuelong Su, Lihui Peng, and Danya Yao. 2010. A review of game theory applications in transportationanalysis. In . IEEE, 152–157.[85] Jing Zhang and Ioannis Ch Paschalidis. 2017. Data-driven estimation of travel latency cost functions via inverseoptimization in multi-class transportation networks. In . IEEE, 6295–6300.[86] Jing Zhang, Sepideh Pourazarm, Christos G Cassandras, and Ioannis Ch Paschalidis. 2018. The price of anarchy intransportation networks: Data-driven evaluation and reduction strategies.
Proc. IEEE
Computational Optimization and Applications
49, 2 (2011), 271–297.
A SUPPLEMENTARY MATERIALSA.1 Variable Number Calculations for Lemma 2.1
In this subsection, we will provide the exact variable counts for problem (2) as the grid size increasesfrom 2x2 to 3x3 to 4x4, and so on, as well as a more general variable count for graphs involving anarbitrary number of arcs. These variable counts are polynomial in the problem input size and, whencompared with the polynomial-time complexity of solving linear programs, supports Lemma 2.1 inthe main paper. We will define 𝑁 as the number of players and 𝑚 as the total number of nodes. Grid:
We first focus on the Grid family of networks, as used in Section 5 in the main paper. Wenote the following: Number of Origin-Destination Pairs: ( ) (cid:18) 𝑚 (cid:19) (22a)Number of Arcs in a √ 𝑚𝑥 √ 𝑚 Grid: (√ 𝑚 − )(√ 𝑚 )( )( ) (22b)For the same costs across all players, we have to consider the variables diag (C) , ¯ 𝑐, 𝑣 𝑘𝑖 , 𝑢 𝑘𝑖 , ¯ 𝑢 𝑘 , ∀ 𝑖, 𝑘 along with the variables utilized in creating the 1-norms for the residual model. This variable countcomes out to be for our code:13 𝑁𝑚 − 𝑁𝑚 / − 𝑁𝑚 + 𝑁𝑚 / + 𝑚 − 𝑚 / − 𝑚 + 𝑚 / (23)For different costs across all players, we consider the variables diag (C 𝑖 ) , ¯ 𝑐 𝑖 , 𝑣 𝑘𝑖 , 𝑢 𝑘𝑖 , ¯ 𝑢 𝑘 , ∀ 𝑖, 𝑘 alongwith the variables utilized in creating the 1 norms for the residual model. This variable count comesout to be for our code:21 𝑁𝑚 − 𝑁𝑚 / − 𝑁𝑚 + 𝑁𝑚 / + 𝑚 − 𝑚 / − 𝑚 + 𝑚 / (24)Therefore, in big O notation, the order of both variable counts is O ( 𝑁𝑚 ) . General graphs:
For any graph, we now define 𝑎 as the number of arcs, with 𝑁 as the number ofplayers and 𝑚 as the total number of nodes. Again, the number of OD pairs will be ( ) (cid:0) 𝑚 (cid:1) . The llen, Dickerson, & Gabriel 𝑎𝑁𝑚 − 𝑎𝑁𝑚 + 𝑚 𝑁 − 𝑚 𝑁 + 𝑎𝑚 − 𝑎𝑚 (25)The variable count for not the same costs across all players is for our code:5 𝑎𝑁𝑚 − 𝑎𝑁𝑚 + 𝑚 𝑁 − 𝑚 𝑁 + 𝑎𝑚 − 𝑎𝑚 (26)Therefore, in big O notation, the order of both variable counts is O ( 𝑎𝑁𝑚 ) . A.2 Proof of Lemma 3.7
Proof. F is defined for (1) as follows: F ( x ) = C 𝑥 + C 𝑥 + ... + C 𝑥 𝑁 + ¯ 𝑐 C 𝑥 + C 𝑥 + ... + C 𝑥 𝑁 + ¯ 𝑐 ... C 𝑁 𝑥 + C 𝑁 𝑥 + ... + C 𝑁 𝑥 𝑁 + ¯ 𝑐 𝑁 (27)In order to prove strong monotonicity, we must show that there exists a scalar 𝛼 > ( F ( x ) − F ( y )) 𝑇 ( x − y ) ≥ 𝛼 || x − y || ∀ x , y ∈ X , x ≠ y . (28)Through some algebra, it is not difficult to show that: ( F ( x ) − F ( y )) 𝑇 ( x − y ) = ( x − y ) 𝑇 C C . . . C C C . . . C ... ... . . . ... C 𝑁 C 𝑁 . . . C 𝑁 𝑇 ( x − y ) (29)If we assume that C = C = ... = C 𝑁 = C , we have = ( x − y ) 𝑇 C C . . . CC C . . . C ... ... . . . ... C C . . . C ( x − y ) (30)From Proposition 2.2.10 in Cottle et al. [17], we know that, for symmetric, real valued matrices 𝑀 and the smallest eigenvalue of 𝑀 referenced as 𝜆 , we can write 𝜆 || 𝑥 || ≤ 𝑥 𝑇 𝑀𝑥, ∀ 𝑥 . Therefore,to prove the matrix in (30) is positive definite, we simply need to prove that its smallest eigenvalueis positive. We begin this process by splitting the matrix in (30) into two pieces 𝐴 = C . . . C . . . ... ... . . . ... . . . C , 𝐵 = C C . . .
CC C . . . C ... ... . . . ... C C . . . C (31)We state that 𝐴 + 𝐵 = 𝑀 , with M as the original matrix, and we note that 𝐴, 𝐵, 𝑀 ∈ R 𝑛𝑁 × 𝑛𝑁 , with 𝑛 as the number of arcs and 𝑁 as the number of players. First, we prove that the eigenvalues of 𝐵 are { , ..., , 𝑁 ( 𝑒𝑖𝑔 ( 𝐶 ))} . For this proof, we notice that the eigenvalue-eigenvector equation for 𝐵 is 𝐵𝑧 = 𝜆𝑧 (32) C 𝑧 + ... + C 𝑧 𝑁 = 𝜆𝑧 𝑖 , 𝑖 = , ..., 𝑁 (33) llen, Dickerson, & Gabriel 𝑧 𝑖 ∈ R 𝑛 , ∀ 𝑖 (not all 𝑧 𝑖 =
0) and 𝜆 ∈ R . Because ∀ 𝑖 the left-hand sides of (33) are the same,we can equate the right hand sides to say: 𝜆𝑧 = 𝜆𝑧 = ... = 𝜆𝑧 𝑁 (34)This provides two cases for the eigenvalues. Either 𝜆 =
0, or the 𝑧 𝑖 are such that 𝑧 = 𝑧 = ... = 𝑧 𝑁 = 𝑤 ≠ 𝑤 ∈ R 𝑛 . This 𝑤 can be substituted back into the equation in (33) to say 𝑁 C 𝑤 = 𝜆𝑤 (36)which means the 𝜆 are the eigenvalues of 𝑁 C and, since we already know that C is a positivediagonal matrix, the eigenvalues are just those entries on the diagonal multiplied by 𝑁 . Becausethese diagonal values are positive and 𝑁 is positive, then the resultant eigenvalues are positive.Next, we use this information to prove that 𝑀 is a positive definite matrix. Since we know that 𝐴 and 𝐵 are real, symmetric matrices, Weyl’s Inequality [42, Theorem 4.3.1] states: 𝜆 ( 𝐴 ) + 𝜆 ( 𝐵 ) ≤ 𝜆 ( 𝐴 + 𝐵 ) = 𝜆 ( 𝑀 ) (37)We know from above that 𝜆 ( 𝐵 ) =
0. We also know that 𝜆 ( 𝐴 ) > C matrices on its diagonal and, since C has all positive values, then we know all of theeigenvalues of 𝐴 are positive. Consequently, 𝜆 ( 𝐴 ) + 𝜆 ( 𝐵 ) > 𝜆 ( 𝑀 ) >
0. As a result,we can write ( x − y ) 𝑇 C C . . . CC C . . . C ... ... . . . ... C C . . . C ( x − y ) ≥ 𝜆 ( 𝑀 )|| x − y || (38)and we know 𝜆 ( 𝑀 ) is a positive number. Thus, we have proven strong monotonicity of F . □ A.3 Counterexample and MATLAB Code for Lemma 3.8
We will provide a counter example to demonstrate this lemma. In MATLAB, if we set the rng seed to1 with the ’twister’ option, we draw four arcs using the rand function and multiply these numbersby 1000. We then form a C matrix by putting these four numbers along the diagonal. Next, we formthe following matrix: 𝐴 = C C C C( . )C ( . )C ( . )C ( . )C( . )C ( . )C ( . )C ( . )C( . )C ( . )C ( . )C ( . )C (39)Therefore, we see that we have formed a four player matrix game in which all of the players do havedifferent costs due to the different factors multiplying the C matrices. According to [46, 82], a matrixwith real values 𝐴 is positive definite if and only if ( 𝐴 + 𝐴 𝑇 ) is positive definite. We can check thesmallest eigenvalue of this resultant matrix to find out if it is positive definite. Using MATLAB’seig function, we discover that the smallest eigenvalue of the matrix ( 𝐴 + 𝐴 𝑇 ) is -7.7307e+04.Therefore, 𝐴 is not positive definite, thus demonstrating that, when the costs are not the sameacross all players, we are not guaranteed to have a strongly monotone F as defined by (5). See belowfor MATLAB code.This is the MATLAB code utilized for the counterexample for Lemma 3.8: llen, Dickerson, & Gabriel rng(1,'twister');v=1000*rand(1,4);C=diag(v);M=[2*C C C C500.1*C 2*500.1*C 500.1*C 500.1*C600.7*C 600.7*C 2*600.7*C 600.7*C700.8*C 700.8*C 700.8*C 2*700.8*C];w=eig(0.5*(M+M'));minval=min(w) A.4 Additional Experimental Set-Up Details
The PATH solver [18, 27] in GAMS generates the data used in the inverse optimization residualmodel [47, 65]. The only default altered for the PATH solver was the tolerance threshold for thesolver to finish, which changed from 1e-6 to 1e-8. In the scripts, the PATH solver can begin at afew different starting points in case one of the start points fails due to solver error. For the casein which the costs are the same across all players, F defined by (5) is strictly monotone, so theproblem has a unique solution, which means the starting point does not matter. For the case inwhich the costs are not the same across all players, multiple starting points are still provided toensure that a solution is obtained to the problem. Interestingly, all of the different cost matrices forour experiments are positive definite, except for one in grid size 5x5 and player number 10. Then,the pyomo [38, 39] package in Python is utilized to construct the inverse optimization residualmodel to find the parameterizations. The Gurobi solver [34] is utilized to solve the pyomo models,and the BarQCPConvTol and the BarConvTol parameters are set to 1e-14. The documentation forBarConvTol states that smaller values for this parameter can lead to “a more accurate solution”[33]. The experiments are run on machines with 4 cores and 32 GB of RAM.The original costs under which the simulation data is generated have to be selected, and thiswas done randomly. MATLAB [55] was utilized to generate the random sets of numbers needed,specifically the unifrnd random number function with 5 as the rng seed. This function draws from auniform distribution the original costs for the C 𝑖 and ¯ 𝑐 𝑖 parameters in the simulation model outlinedin (20). For the same costs across all players set-up, one C and one ¯ 𝑐 are drawn but, when the costsvary across all the players, a C 𝑖 and a ¯ 𝑐 𝑖 are drawn for each player 𝑖 . For each experimental set upwith a specified graph, number of players, and 𝛼 level (the capacity level for the joint constraint),10 trials are run. A.5 Positive Definiteness of Different Costs Matrices with MATLAB rng(5)
In Figure 7, we plot a sampling of the minimum eigenvalues for the “symmetric part” [82] of randommatrices of the form: 𝐴 = C C . . . C C C . . . C ... ... . . . ... C 𝑁 C 𝑁 . . . C 𝑁 (40) In different costs, we use slightly different starting point options between when we run the simulation for the originalcosts and when we run the simulation for the IO costs. This doesn’t appear to affect the results, as will be seen in the nextfew subsections. We also did a check for 𝛼 =
10, grid 5x5, and 𝑁 =
10 with the start points set to the same values for theoriginal costs and IO costs, and the resulting errors were the same as before.llen, Dickerson, & Gabriel C 𝑖 are diagonal matrices whose diagonals are chosen as random numbers from the uniformdistribution from 1 to 5 using the MATLAB unifrnd function. We choose to set the arc number to76 (same number as the Sioux Falls arc number), and we gradually increase the number of players 𝑁 from 2 to 15. We do use the seed of 5 for rng but, since we draw all of the numbers in onecontinuous session, the minimum eigenvalues here do not necessarily match the eigenvalues of thematrices that correspond with our experiments. As discussed in Appendix A.3 and according to[46, 82], we can check the positive definiteness of the matrix 𝐴 by finding the minimum eigenvalueof 0 . ( 𝐴 + 𝐴 𝑇 ) , which is what we do for the 10 random matrices corresponding with each playersetting. In return, we obtain Figure 7, and we see that, as player number increases, there is moreand more of a chance that the minimum eigenvalue is negative, thus making the matrix not positivedefinite.This is important because we can only guarantee unique solutions to our problem when the 𝐴 matrix is positive definite, as explained in Section 3.1. This allows us to have multiple startingpoints in our code for the simulation experiments. Fig. 7. Minimum Eigenvalues for 10 Trials, 𝑁 = − , rng(5), and Arc Number=76, which is the same numberof arcs as Sioux Falls A.6 Objective Function Values for the Experiment Groups
Overall, it is important to remember that, for all the experiments, the objective function valuesare small, on the order of 1e-6 or lower. However, it is also relevant to identify patterns wherepresent. For experiment group 1, in Figure 8, increasing grid size sometimes leads to a highermedian objective function value in certain subsets of the experiments, including when there are 2players and 𝛼 =
1, when there are 5 players and 𝛼 = .
5, and when there are 10 players and 𝛼 = 𝛼 = For experiment group 2, in Figure 9, the objective valuesare overall quite low, with one outlier for 10 players and 𝛼 = Note that there are a few negative values, but this isn’t concerning because their corresponding flow errors are stilllow.llen, Dickerson, & Gabriel For experiment group 3, in Figure 10, although the objective function values of the inverseoptimization residual model are on the order of 1e-6, the median values under the 𝛼 = ( . ) 𝑁 capacity are bigger than the 𝛼 = 𝑁 median values. One theory for this trend is that the moreconstrained problems are more difficult to solve and hence lead to larger objective function values.The ranges for these objective function values is about 1.5e-6 higher in the positive direction thanthe corresponding range for the values in Experiment group 1 (see Figure 8). For experiment group4, Figure 11 demonstrates that most of the values fall below 1e-6, except for the 𝑁 = , 𝛼 = . Fig. 8. Objective Function Values for Experiment Group 1: The labeling at the bottom of the graph indicatesattributes of the boxplot, specifically the Grid Size/Number of Players/Alpha Value There do again appear to be some negative values but, again, this isn’t concerning because their corresponding flowerrors are still low.llen, Dickerson, & Gabriel Fig. 9. Objective Function Values for Experiment Group 2: The labeling at the bottom of the graph indicatesattributes of the boxplot, specifically the Number of Players/Alpha ValueFig. 10. Objective Function Values for Experiment Group 3: The labeling at the bottom of the graph indicatesattributes of the boxplot, specifically the Number of Players/Alpha Value. Note: Only 6 trials are included for5/10/5.0, see note in the text. llen, Dickerson, & Gabriel Fig. 11. Objective Function Values for Experiment Group 4: The labeling at the bottom of the graph indicatesattributes of the boxplot, specifically the Number of Players/Alpha Value. Note: The trials did not finish intime for 𝑁 = and 𝛼 = , hence that boxplot is not included. A.7 Timing for the Experiment Groups
In this subsection of the Appendix, we present the timing graphs for the experiment groups. InFigure 12, timing increases with grid size for all of the subsets of 𝛼 and number of players. In Figure13, the timing increases as the player number increases. Both of these trends are understandablebecause of the effect on variable count that both grid size and player number have, as demonstratedby Appendix A.1. In Figure 14, we see that all of the timing data is obscured by the large outlier ofthe experiment with grid 5x5, 10 players, and 𝛼 =
5. This experiment took quite a while and, aspreviously noted, only 6 of the trials for the experiment finished in under 24 hours (the 1st, 2nd,3rd, 4th, 5th, and 7th). We hypothesize that having the largest grid, the most number of players,the more constrained 𝛼 value of the two, and the more difficult task of determining different sets ofcost parameters for each player caused the large time amounts. In Figure 15, again noting that thetrials did not finish in time for 𝑁 =
10 and 𝛼 =
5, it appears that there are not trends to discuss,although note that all the experiments but 𝑁 = , 𝛼 = . llen, Dickerson, & Gabriel Fig. 12. Timing for Experiment Group 1: The labeling at the bottom of the graph indicates attributes of theboxplot, specifically the Grid Size/Number of Players/Alpha ValueFig. 13. Timing for Experiment Group 2: The labeling at the bottom of the graph indicates attributes of theboxplot, specifically the Number of Players/Alpha Value llen, Dickerson, & Gabriel Fig. 14. Timing for Experiment Group 3: The labeling at the bottom of the graph indicates attributes of theboxplot, specifically the Grid Size/Number of Players/Alpha Value. Note: Only 6 trials are included for 5/10/5.0,see note in the text.Fig. 15. Timing for Experiment Group 4: The labeling at the bottom of the graph indicates attributes of theboxplot, specifically the Number of Players/Alpha Value. Note: The trials did not finish in time for 𝑁 = and 𝛼 = , hence that boxplot is not included. B CODE ATTRIBUTION
In this section, we would like to cite the various packages and software that we used in the project.We would also like to note that the code for this project was built off of another code base that isfor another on-going project (not published yet). • Python Packages: – pyomo [38, 39] – networkx [35] – matplotlib [44] llen, Dickerson, & Gabriel – numpy [60, 77, 78] – scipy [77] – pandas [56] • MATLAB 9.8.0.1323502 (R2020a) [55] • GAMS [30, 31] with PATH solver [18, 27] (PATH website [28]) • Solvers: – gurobi Version 9.1.1 [34]; Gurobi Documentation [52] • Data: Sioux Falls data [51, 76]These are some of the helpful resources we utilized for writing our code: • These articles [8, 53] were helpful in determining some of the constraints that equalizednecessary variables in my code. ••