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?