A General Cooperative Optimization Approach for Distributing Service Points in Mobility Applications

: This article presents a cooperative optimization approach (COA) for distributing service 1 points for mobility applications, which generalizes and reﬁnes a previously proposed method. 2 COA is an iterative framework for optimizing service point locations, combining an optimization 3 component with user interaction on a large scale and a machine learning component that learns 4 user needs and provides the objective function for the optimization. The previously proposed 5 COA was designed for mobility applications in which single service points are sufﬁcient for 6 satisfying individual user demand. This framework is generalized here for applications in which 7 the satisfaction of demand relies on the existence of two or more suitably located service stations, 8 such as in the case of bike/car sharing systems. A new matrix factorization model is used as 9 surrogate objective function for the optimization, allowing to learn and exploit similar preferences 10 among users w.r.t. service point locations. Based on this surrogate objective function, a mixed 11 integer linear program is solved to generate an optimized solution to the problem w.r.t. the 12 currently known user information. User interaction, reﬁnement of the matrix factorization, and 13 optimization are iterated. An experimental evaluation analyzes the performance of COA with 14 special consideration of the number of user interactions required to ﬁnd near optimal solutions. 15 The algorithm is tested on artiﬁcial instances as well as instances derived from real-world taxi data 16 from Manhattan. Results show that the approach can effectively solve instances with hundreds 17 of potential service point locations and thousands of users while keeping the user interactions 18 reasonably low. A bound on the number of user interactions required to obtain full knowledge of 19 user preferences is derived, and results show that with 50% of performed user interactions the 20 solutions generated by COA feature optimality gaps of only 1.45% on average. 21


Introduction
A fundamental ingredient for optimizing the locations of service points in mobility 25 applications, such as charging stations for electric vehicles or pickup and drop-off 26 stations for car/bike sharing systems, is the distribution of existing customer demand 27 to be potentially fulfilled in the considered geographical area. While there exists a 28 vast amount of literature regarding setting up service points for mobility applications, worthless if there does not also exist a suitable location near the destination for returning 48 the vehicle. Furthermore, some customers might use multiple modes of transport for 49 a single trip [11]. Consequently, some more distant service station for returning the 50 vehicle might be acceptable if this place is well connected by public transport used for an 51 additional last leg [12]. Thus, there also might be alternatives for fulfilling demand that 52 cannot all be exactly specified by potential users. The example with an additional leg by 53 public transport also illustrates that geographical closeness is not always the deciding 54 factor. 55 To address these issues Jatschka et al. [13] proposed the concept of a Cooperative data acquisition [14], the direct integration of users into the whole location planning 64 process, a possibly stronger emotional link of the users to the product, and ultimately that solutions generated by COA feature optimality gaps of only 1.45% on average with 93 50% of performed user interactions. 94 The next Section discusses related work. In Section 3 we introduce our formal 95 problem setting, the Generalized Service Point Distribution Problem (GSPDP), which reflects 96 the essence of various location problems for mobility services. Afterwards, in Section 97 4, we detail the cooperative optimization framework for solving the GSPDP. Section 98 5 describes how benchmark scenarios for testing were generated, and in Section 6 we 99 present and discuss experimental results. Section 7 concludes this work with an outlook 100 on promising future work. given potential locations as economically as possible in order to fulfill a certain level of 105 customer demand. Similarly to a p-median or maximal covering location model, see 106 [18], we limit the facilities to be opened. However, instead of imposing a direct limit on 107 the number of facilities, in our formulation a maximum budget for setting up facilities 108 is specified and each facility is associated with setup costs. In general our problem can 109 be classified as an uncapacitated fixed charge Facility Location Problem (FLP), see [19], 110 however, without explicitly given demands (see below). Moreover, our problem also 111 has similarities to stochastic or robust FLPs [20] in which problem inputs, such as the 112 demand, may be uncertain. For such problems uncertain inputs are usually modeled via 113 random variables [21] or scenario-based approaches [22]. For a survey on FLPs see [23]. 114 In our problem formulation we have mobility applications in mind, such as the such systems to the best of our knowledge all existing work essentially assumes customer 118 demand to be estimated upfront. For example, [24] and [8] use parking information to 119 identify promising locations for electric vehicle charging stations. [25] locate charging 120 stations for an on-demand bus system using taxi probe data of Tokyo. Moreover, census 121 data are commonly used for estimating demand for car sharing systems [26], bike sharing 122 systems [27], or for setting up electric vehicle charging stations [28]. Data mining and 123 other related techniques are also often employed to detect traffic patterns for bike sharing 124 systems [29,30] as well as car sharing systems [31]. 125 Ciari et al. [32] recognize the difficulties in making demand predictions for new 126 transport options based on estimated data and therefore propose to use an activity-based 127 microsimulation technique for the modeling of car sharing demand. The simulation is 128 done with help of the travel demand simulator MATSim [33]. 129 There also exist some works that take user preferences into account, e.g., in [34] a 130 car sharing system is designed based on different assumptions on the behavior of users. 131 The authors come to the conclusion that providing real-time information to customers 132 can greatly improve the service level of a system if users are willing to visit a different 133 station if their preferred one does not have a vehicle available. 134 In our approach we substantially deviate from this traditional way of acquiring 135 existing demand upfront and instead resort to an interactive approach. Potential future 136 customers are directly incorporated in the optimization process as an integral part by 137 iteratively providing feedback on meaningfully constructed location scenarios. In this 138 way we learn user demands on-the-fly and may avoid errors due to unreliable a priori 139 estimations. For a survey on interactive optimization algorithms in general see [14]. 140 The performance of interactive algorithms is strongly influenced by the quality of the 141 feedback given by the interactors. Too many interactions with a user will eventually 142 result in user exhaustion [35], negatively influencing the reliability of the obtained 143 feedback. Additionally, interacting with users can be quite time consuming, even when 144 dealing with a single user. Hence, in order to keep the interactions with users low one 145 can resort to surrogate-based optimization approaches [36,37]. Surrogate models are 146 typically machine learning models serving as proxy of functions that are difficult to 147 evaluate [38]. In [16] as well as in this contribution we make use of a matrix factorization 148 [39] based surrogate model. Matrix factorization is a collaborative filtering technique 149 which is frequently used in recommender systems, see e.g., [40]. 150 As already pointed out, the basic concept of COA was already presented in [13], 151 where we made use of an adaptive surrogate model [41]. The underlying structure 152 of this surrogate model is formed by a large set of individual smaller models, i.e., 153 one machine learning model for each combination of user and potential service point 154 location, which are trained with the feedback of the respective users. While initially 155 being instantiated by simple linear models, these machine learning models are step-wise 156 upgraded as needed to more complex regressors during the course of the algorithm in 157 order to cope with possibly encountered higher complexity. Specifically in [13], each 158 linear model can be upgraded to a neural network. The number of neurons in the hidden 159 layer is increased whenever the training error exceeds a certain threshold. Overfitting Variable Neighborhood Search (VNS) [42] was used as optimization core to generate 162 new solutions w.r.t. the current surrogate function. In [43] the performance of the VNS 163 optimization core was compared to an optimization core using a population-based 164 iterated greedy approach [44]. Unfortunately, this first realization of COA exhibits severe 165 limitations in the scalability to larger numbers of potential service point locations and/or 166 users, in particular as all users are considered independently of each other.

167
The current work builds upon the observation that in a larger user base there are 168 typically users sharing the same or similar needs or preferences. Identifying these shared 169 demands and exploiting them to improve scalability as well as to reduce the required 170 feedback per user is a main goal here. While the basic principles of COA remain the 171 same, major changes are performed in the way the approach interacts with users, how 172 the feedback of users is processed, as well as how new candidate solutions are generated.

173
Moreover, the adaptive surrogate model from [13] realized by the large set of underlying 174 simpler machine learning models is replaced by a single matrix factorization based 175 model that is able to exploit said similarities between users. Besides our aforementioned 176 preliminary conference paper already sketching the application of a matrix factorization 177 within COA [16], to the best of our knowledge there exists no further work on interactive 178 optimization approaches for location planning in mobility applications. unit of satisfied customer demand a prize q > 0 is earned.

189
A solution to the GSPDP is a subset X ⊆ V of all locations where service points are to be set up. A solution X is feasible if its total fixed costs do not exceed the maximum budget B, i.e., Given the set of users U, we assume that each user u ∈ U has a certain set of use all SPRs over all users. Note that in this notation, different users never share the same 220 SPR labels, although labels may refer to similar SPRs.

221
For indicating how suitable a location is w.r.t. to an SPR we define values w r,v ∈ [0, 1]

222
indicating the suitability of a service point at location v ∈ V to satisfy the needs of user 223 u ∈ U concerning SPR r ∈ R u,c in the use case c ∈ C u . A value of w r,v = 1 represents 224 perfect suitability while a value of zero means that location v is unsuitable; values in 225 between indicate partial suitability.

226
With these suitability values in mind, the objective of the GSPDP is to maximize In the first term of this objective function, the obtained prize for the expected total 227 satisfied demand is determined by considering for each user u, each use case c, and each 228 SPR r a most suitable location v ∈ V at which a service point is to be opened (v ∈ X).

229
Over all SPRs of a use case, the minimum of the obtained suitability values is taken 230 so that the full demand is only fulfilled when for each SPR an ideally suited service 231 station is planned, and no demand is fulfilled as soon as one of the SPRs does not have 232 an appropriate service point. The second term of the objective function represents the 233 total maintenance costs for the service stations.

234
By linearizing the above objective function, the GSPDP can be formulated as a mixed integer linear program (MILP) with the following variables. Binary variables x v indicate whether or not a service point is deployed at location v ∈ V, i.e., the binary vector x = (x v ) v∈V is the incidence vector of a corresponding solution X ⊂ V. Additional variables h r,v are used to indicate the actually used location v ∈ V for each SPR r ∈ R.
The degree to which a use case c ∈ C u of a user u ∈ U can be satisfied is expressed by continuous variables y u,c ∈ [0, 1]. The GSPDP is then stated as follows.
In correspondence to the definition of f , the objective value is calculated in (3)   the GSPDP instance is set to p, and the prize for a unit of covered demand q is set to 260 one. The suitability value w r,v is set to one for v ∈ V and r ∈ R if facility i can satisfy the 261 demand of demand node j, i.e., j ∈ F i , and zero otherwise.  the quality of obtained feedback when too much information is asked from a user is, for 291 example, discussed in [35].

292
Hence, interaction with users needs to be kept to a minimum and should be done 293 wisely to extract as much meaningful information as possible. Moreover, users must be 294 confronted with easy questions whose answers at the same time provide strong guidance 295 for the target system. Based on this philosophy, we present a Cooperative Optimization  In Section 4.3 we formally define how users can interact with COA. Note however 302 that COA does not put a strict limit on the number of allowed interactions with each 303 user. Instead, we will measure the effectivity of our COA framework by the number of 304 user interactions the framework requires for generating a (close to) optimal solution.  First, the FC is called starting an initialization phase by asking each user u ∈ U to 321 specify the user's use cases C u , associated SPRs R u,c , as well as corresponding demands   him in respect to a use case requirement r ∈ R u .

328
The feedback obtained from the users is processed in the EC. The EC maintains and continuously updates a surrogate suitability functionw Θ (r, v) approximating the suitability values w r,v of service point locations v ∈ V w.r.t. SPR r ∈ R without interacting with the respective user. This function is realized by a machine learning model with parameters  (weights) Θ. Based on this surrogate function, the EC also provides the surrogate objective functioñ with which a candidate solution X can be approximately evaluated. This fast approxi-   The whole process is repeated until some termination criterion is reached, e.g., the 343 discrepancy of user feedback and the results of the EC is small enough or a maximum   In the next subsections we describe each component's functionality in more detail. Input : an instance of the GSPDP Output: best solutionX * ⊆ V found 1: Feedback Component: 2: for u ∈ U do 3: obtain use cases C u with associated demands D u,c and service point requirements R u,c , ∀c ∈ C u , from user u; The FC is responsible for extracting as much information as possible from the users 366 with as few interactions as necessary in order to not to fatigue them.

367
In the context of the COA framework, user interactions are understood as letting

374
Let S ⊆ V be such a location scenario provided to a user u in respect to one of 375 the user's SPRs r ∈ R u . It is then assumed that the user returns as evaluation of the 376 scenario S either a best suited location v r,S ∈ S and the corresponding suitability value 377 w(r, v r,S ) > 0 or the information that none of the locations of the scenario S is suitable.

378
The latter case implies that w r,v = 0 for all v ∈ S. In case multiple locations are equally 379 well suited, we assume that the user selects one of them at random. It is assumed here 380 that the suitability of a location w.r.t. an SPR can be specified by the user on a five-valued 381 scale from zero, i.e., completely unsuitable, to one, i.e., perfectly suitable; a more fine 382 grained evaluation would not make much practical sense.

383
Clearly, this definition of user interaction is simplified and idealized, in particular 384 as we assume here that all users always give precise answers. In a real application, the 385 uncertainty of user feedback and the possibility of misbehaving users who intentionally 386 give misleading answers also need to be considered among other aspects. Moreover, it 387 would be meaningful to extend the possibilities of user feedback. For example, users 388 could be allowed to optionally rate more than one suitable locations for an SPR in one 389 scenario or to make suggestions which locations to additionally include in a scenario as 390 we considered it in [13]. thus, w r,v r,S ≥ w r,v , ∀v ∈ S. By w UB r,v , the SMC maintains for each SPR r ∈ R, and each 400 location v ∈ V the so far best obtained upper bound on each w r,v ; initially, w UB r,v = 1. that it can only be guaranteed that V w (r) is completely known once the user returns that 408 none of the locations in the last scenario S V r are suitable for r. Consequently, V w (r) will 409 be completely known after |V w (r)| + 1 user interactions.

410
Hence, an upper bound I UB u on the total number of required interactions with user u for completely identifying all relevant locations for all of his/her use cases is While this value is unknown in a real-world scenario, it allows us to establish a measure 411 of quality on how well our strategy for presenting scenarios to users performs within 412 our testing environments.

413
The same combination of two strategies for generating scenarios w.r.t. an SPR r ∈ R 414 as in [16] is used. The first strategy generates scenarios according to the approach according to S V r and S * r , respectively, with ς V and ς * being strategy parameters.

427
The EC processes the user feedback obtained from the FC and provides the means for evaluating candidate solutions without relying on users in particular within the OC. The exact objective function f from (2), which is based on the mostly unknown suitability values w r,v with r ∈ R, v ∈ V, is approximated by the surrogate objective functionf Θ , cf. (11), making use of the following surrogate suitability functioñ preferences is matrix factorization [39], which we also apply here. Given an incomplete 439 matrix containing ratings R = (w i,j ) i∈U, j∈P for a set of users U over a set of products 440 P, the idea behind matrix factorization is to decompose this matrix into two smaller 441 matrices, a user/feature matrix ξ and a product/feature matrix ν, such that the product 442 of these two matrices approximates the original matrix. An unknown rating, i.e., a rating 443 not contained in the original matrix R, can then be estimated as the dot product of the 444 corresponding feature vectors in matrix ξ and matrix ν, respectively.

445
Moreover, we also want to exploit the fact that only a small fraction of the locations in V is typically relevant for the SPR of a user and that unknown ratings are not missing at random. In our problem users are always asked to rate the most suitable location of a scenario. Therefore, known ratings tend to be biased towards more positive values while unrated locations are likely to have a low suitability for a user w.r.t. an SPR. A matrix factorization approach that takes such considerations into account has been suggested by [17]. Traditionally, the rating matrix R is factorized by solving the optimization problem where E is a loss function for measuring the error between the actual and the predicted ratings and ρ is a regularization term. In [17] this minimization problem is expanded by adding a bias term for unknown ratings towards a certain valueŵ, i.e., values are more likely zero than being greater than zero, we set the bias targetŵ = 0.

461
Having obtained matrices ξ and ν, an unknown value of W is approximated by the dot product of the respective feature vectors rounded to the nearest of the five discrete suitability values we defined, i.e., The trainable parameters of g Θ are therefore Θ = (ξ, ν).

462
Our loss function for the matrix factorization is and randomized block coordinate descent [47] is used to minimizing it. go beyond the scope of this paper, we therefore leave this task for future work.

478
In the following it is described how test instances for evaluating COA have been

495
Test instances from the groups EVC and CSS are generated with the same approach. The budget is assumed to be B = 7.5 · n so that roughly 10% of the stations with 501 average costs can be expected to be opened.

502
The number of use cases for each user u ∈ U is chosen randomly according to Each SPR r ∈ R u,c of a use case c also is associated with a particular geographic location q r ∈ {0, . . . , L − 1} 2 . In order to model similarities in the users' SPRs, these locations are selected in the following correlated way. First ten attraction points A with uniform random coordinates are selected from {0, . . . , L − 1} 2 . Then, each use case location is derived by randomly choosing one of these attraction points (a x , a y ) ∈ A and adding a small individual offset to the coordinates, i.e., where N (·, ·) denotes a random value sampled from a normal distribution with the 508 respectively given mean value and standard deviation σ v . If obtained coordinates are 509 not in {0, . . . , L − 1} 2 a new attraction point is chosen and the deviation is re-sampled.

510
A service point location v ∈ V receives a rating w.r.t. an SPR r according to a sigmoidal decay function applied to the Euclidean distance, and is also perturbed by a Gaussian noise with a standard deviation of σ r : The parameters of the sigmoid function are chosen so that w r,v decreases as the dis-511 tance between v and q r increases and becomes approximately zero at a distance larger 512 than twelve. Additionally, as motivated in Section 3, we discretize the rating w r,v by 513 rounding to the closest value in {0, 0.25, 0.5, 0.75, 1}, obtaining w r,v . Hence, w r,v = 514 4 · min(1, max(0, w r,v )) + 0.5 /4.

515
Six sets of 30 benchmark instances were generated for EVC as well as CSS. As 516 detailed in Table 1, these sets consider n ∈ {100, 200, 300} potential service point loca-

548
The number of use cases for each user u ∈ U is again chosen randomly according to a shifted Poisson distribution with offset one, expected value three, and a maximum value of five. Each of these use cases c ∈ C u is associated with an individual demand D u,c randomly chosen from {5, . . . , 50} and the two SPRs representing the origin and destination of a trip chosen from Q uniformly at random. A rating for an SPR r is calculated for each v ∈ V via the sigmoidal decay function where sp(r, v) refers to the length of the shortest path between location v and the SPR

561
The COA framework including the FC were implemented in Python 3.8.   Recall that the number of attraction points is the same for all EVC and CSS instances.

611
Therefore, for instances with a higher number of users it is generally easier to find better Note that an increased number of user interactions does not only imply a larger trainings set for the surrogate function but also results in better upper bounds w UB r,v for locations v ∈ V w.r.t. to an SPR r ∈ R. Therefore, to gain a better understanding of how much the surrogate function actually contributes to finding an optimized solution, we study what happens when the learning surrogate suitability functionw Θ is replaced by the naive function with no learning capabilities   MAN Figure 6. Average optimality gaps with standard deviations as shaded areas obtained by COA using the surrogate suitability functions with learning capabilities (w Θ ) and without (w bl ) plotted over the interaction level.
Therefore, the plots on the right side of Figure 7   Finally, in Table 3

693
In this paper the previously introduced Cooperative Optimization Algorithm was 694 generalized to be applicable to more application scenarios and to larger instances with