Start with a binary variable x_{i, k} that takes value 1 if city i is in group k and 0 otherwise.
Each city must be in a group, so we add the constraint \sum_{k} x_{i, k} = 1 for every i.
The total population of group k is then Q_{k} = \sum_{i} x_{i, k} q_{i}.
Adding bounds Q_{min} \leq Q_{k} \leq Q_{max} (as you suggested initially) will ensure that groups have similar population levels.
Then, there’s the question of the distance. We want to minimize the sum of distances between pairs of cities in a same group.
Let z_{i, j} be a binary variable that takes value 1 if cities i and j are in the same group, and 0 otherwise. Then the total distance is just \sum_{i, j} d_{i, j} z_{i, j}. (note that I’m counting each pair twice here).
To ensure that z_{i, j} = 1 if and only if cities i and j are in the same group, one way is to add the constraints z_{i, j} \geq x_{i, k} + x_{j, k} - 1 for every pair i, j and every k.
It is easy to verify that, if cities i, j are in group k, then z_{i, j} >= 1+ 1 - 1 = 1.
Otherwise, x_{i, k} + x_{j, k} cannot exceed 1, thus z_{i, j} \geq 0 and, since we are minimizing and z_{i, j} has positive objective coefficient, we get z_{i, j} = 0 in any optimal solution.
The size of this model is thus:
- N \times k + k + N \times N variables
- N + N \times N \times k constraints
I tried this with 10 cities, CPLEX solved the model in a fraction of second.
For 30 cities, it took CPLEX 10s.
I let you imagine what happens with 3,000 cities…