# Clustering a large ride-matching problem

Background: We are solving a large scale vehicle to person ride-matching problem. The problem is essentially simple (match every person with a vehicle, if possible), yet the problem size is quite large, e.g., 200,000 vehicles and 200,000 people. The problem network has already been divided into 2,000 zones based on geographic attributes, but the zone level assignment is too restrictive as zones are so small and enough number of vehicles may not always be found.

The question: How can we combine zones into clusters such that the deviation between the number of vehicles and the number of people within each cluster is minimized and adjacent zones are clustered (i.e., zones far apart from each other would not fall into the same cluster)? We may assume we would like n-clusters and are not interested in optimizing n.

My failed attempt: Let $$I$$ be a set of zones and $$C=\{0,1,2,...,n\}$$ be a set of clusters. The binary variable, $$y_{ic}=1$$, if zone $$i$$ is assigned to cluster $$c$$. Let $$T_{ij}$$ denote the zonal (euclidean) travel time between zone $$i\in I$$ and zone $$j\in I$$ (considering their centroids). Let $$P_i$$ and $$V_i$$, respectively, denote the number of people and the number of vehicles in zone $$i\in I$$. Finally, let $$\epsilon\in\mathbb{R}^+$$ and $$u\in [0,1]$$ be the maximum allowed deviation between the number of vehicles and people, and the uniformity ratio (i.e., how uniform we want the cluster sizes to be), respectively. I have the following model that I know so far away (indeed, wrong) from what I want.

\begin{alignat}2\min &\quad \sum_{\substack{i,j\in I,\\c\in C}}T_{ij}y_{ic}\tag1\\\text{s.t.}&\quad \sum_{c\in C}y_{ic}=1 \qquad\forall i\in I\tag2\\&\quad\sum_{i\in I}y_{ic} \leq (1+u)\frac{|I|}{|C|} \qquad \forall c\in C\tag3\\&\quad \sum_{i\in I}\left(P_i-V_i\right)y_{ic} \leq \epsilon \qquad \forall c\in C\tag4\\&\quad y_{ic}\in \{0,1\}.\end{alignat}

I expected the above model to collect the zones that are close to each other, but surely it didn't. I have the implementation of the above (ignoring (4)) on a dummy sample ($$|I|=20$$ and $$|C|=2$$ with $$u=0$$) and produces the below output. Different colors indicate which cluster the zone is assigned. In this case (where constraint (4) is ignored), I would at least expect to see a nice partition where no zones on the upper side are red and no zones below are blue. Can the correct objective be $${\displaystyle\min\sum_{\substack{i,j\in I,\\c\in C}}|T_{ij}y_{ic}-T_{ji}y_{jc}|^2}$$ ? If yes, any way to linearize this?

Alternative solutions: I also searched for algorithmic solutions and ended up finding Ward's method of hierarchical clustering and implemented the code following Example 1 in here. It partitions the zones perfectly from a visual standpoint (see below plot) but does not allow factoring in constraints (3)–(4). I have also seen this and that which does not exactly help me. Any suggestion?

For $$i, let $$z_{i,j} \ge 0$$ indicate whether zones $$i$$ and $$j$$ are assigned to the same cluster. You want to minimize $$\sum_{i, with additional constraints $$y_{i,c} + y_{j,c} - 1 \le z_{i,j} \quad \text{for i Note that $$z$$ will automatically be integer-valued without explicitly declaring it to be binary. Because there are $$\binom{|I|}{2}|C|$$ of these constraints, you might want to introduce them dynamically only if they are violated.
Regarding constraint $$(4)$$, do you also want $$\ge -\epsilon$$?
Alternatively, you might think of this in terms of political districting, as in the recent paper Imposing contiguity constraints in political districting models. For $$i,j\in I$$, binary decision variable $$x_{i,j}$$ indicates whether zone $$i$$ is assigned to zone $$j$$, and $$x_{j,j}=1$$ means that zone $$j$$ is the "capital" of the resulting district.
• Rob, first thanks for responding. A couple of things to understand: 1) Why $i<j$? Symmetry-breaking? 2) By saying dynamically, what type of constraints these would be corresponding in Gurobi layout? I have seen lazy-constraints and other different things in Gurobi examples but have never tried understand or implement them. 3) Can I also treat $y_{ic}$ as a continuous variable to possibly speed up? 4) Is $\geq -\epsilon$ for $|P_i-V_i|$? If yes, sure! May 22, 2020 at 0:30
• If you use both $z_{i,j}$ and $z_{j,i}$, these will take the same values, so you would unnecessarily double the number of these variables. I don't use Gurobi, but yes, I mean lazy constraints. You can relax $y_{ic}$, but you will likely get a fractional solution, which you will then need to somehow repair or find a way to interpret. May 22, 2020 at 0:34
• Forgetting about Gurobi, does dynamically introduce mean something like: 1) Solve the problem without these constraints. 2) Evaluate the solution and see if any $(i,j)$ pair violates the constraint. 3) If yes, introduce the constraint for that pair and resolve. 4) Repeat 1,2,3 until answering three as a "No." Wouldn't this be computationally even worse as we check the violation after every run? Maybe, depends on the case... May 22, 2020 at 0:40