k-Means Clustering

Partition observations into k clusters by minimizing the within-cluster squared distance, a nonconvex model that, unlike the classical algorithm, also accepts must-link, cannot-link and cluster-size constraints.

Notebook
Python / Knitro API Python / Pyomo Julia / JuMP

Introduction

\(k\)-means clustering is a classical problem in data science. It consists in finding a partition of a given set of observations into a given number of clusters.

\(k\)-means clustering being a classical problem in data science, algorithms to solve it are already available, and they are fast. The advantage of nonlinear programming shows when additional constraints need to be taken into account. These extensions are usually not supported by the algorithms available in the data science libraries. Here we illustrate three of them and how to model them with nonlinear programming:

  • Force some pairs of observations to be in the same cluster
  • Force some pairs of observations to be in different clusters
  • Impose a minimum and maximum size for clusters

Problem description

Inputs

  • \(N\) observations with their coordinates \(X_i\), \(i = 1, \dots, N\) in the considered space
  • A number of clusters to build \(K\)

Problem: partition the \(N\) observations into \(K\) clusters such that each observation belongs to exactly one cluster.

Objective: minimize the total sum of Euclidean distances between each observation and the center of its cluster

\[ \min \sum_{k=1}^K \sum_{i \in \mathrm{cluster}_k} \| X_i - c_k \|^2_2 \]

with \(c_k\) the coordinates of the center of cluster \(k\).

Mathematical model

Now, we show how to model the \(k\)-means clustering problem using nonlinear programming.

Variables

  • \(a^k_i \in \{0, 1\}\), \(i = 1, \dots, N\), \(k = 1, \dots, K\): a boolean variable which is 1 if the observation \(i\) is in cluster \(k\), else 0
  • \(d^k_i \in \mathbb{R}^+\), \(i = 1, \dots, N\), \(k = 1, \dots, K\): a positive variable which represents the Euclidean squared distance between observation \(i\) and cluster \(k\)
  • \(c^k \in \mathbb{R}^2\), \(k = 1, \dots, K\): position of the barycenter of the cluster \(k\)

Objective

\[ \min \sum_{i=1}^N \sum_{k=1}^K a^k_i d^k_i \]

Constraints

  • Each observation has to be in one cluster

\[ \forall i = 1, \dots, N \qquad \sum_{k=1}^K a^k_i = 1 \]

  • Coordinates of the barycenter of the clusters

\[ \forall k = 1, \dots, K \qquad c_k \sum_{i=1}^N a^k_i = \sum_{i=1}^N a^k_i X_i \]

  • Square distance formula

\[ \forall i = 1, \dots, N, \quad \forall k = 1, \dots, K \qquad d^k_i = \| X_i - c_k \|^2_2 \]

The products \(a^k_i d^k_i\) and \(c_k a^k_i\) make this a nonconvex mixed-integer nonlinear program.

Model extensions

In some cases, additional constraints must be taken into account. We add the following sets of constraints to the model above.

Cannot-link and must-link observations. In some cases, some pairs of observations might be known to belong to the same cluster or to different clusters. Let \(P^\text{must}\) be the set of sets of observations that must belong to the same cluster, and \(P^\text{cannot}\) the set of pairs of observations that must belong to different clusters:

\[ \forall k = 1, \dots, K, \quad \forall S \in P^\text{must}, \quad \forall i \in S \qquad a^k_i = a^k_{S[0]} \]

\[ \forall k = 1, \dots, K, \quad \forall i_1, i_2 \in P^\text{cannot} \qquad a^k_{i_1} + a^k_{i_2} \leq 1 \]

Constrained clustering. Another common case in clustering is to have lower and/or upper bounds on the cluster sizes. Let \(u_\text{limit}\) be the maximum number of observations in a cluster and \(l_\text{limit}\) the minimum:

\[ \forall k = 1, \dots, K \qquad l_\text{limit} \leq \sum_{i=1}^N a^k_i \leq u_\text{limit} \]

Input data

An instance is the coordinates of the observations and the number of clusters to build. Clusters carries total_distance, the total squared distance the solve reached, along with the cluster label of every observation and the cluster centers.

from dataclasses import dataclass


@dataclass
class Clusters:
    total_distance: float
    labels: list
    centers: list

    @property
    def num_clusters(self):
        return len(self.centers)

    def cluster(self, k):
        return [i for i, label in enumerate(self.labels) if label == k]

The observations are a CSV of X and Y columns under shared/data. link_groups picks the must-link and cannot-link sets out of the positions themselves rather than naming observations by index.

import csv
from pathlib import Path

from plot import draw_clusters, draw_link_groups, draw_observations
from report import report_clusters

DATA_DIR = Path("shared") / "data"


def load_observations(name):
    with (DATA_DIR / name).open(encoding="utf-8") as handle:
        return [(float(row["X"]), float(row["Y"])) for row in csv.DictReader(handle)]


def link_groups(positions, band_y=20, band_x=(0, 20), split=10):
    left, right = [], []
    for i, (x, y) in enumerate(positions):
        if y >= band_y and band_x[0] < x < band_x[1]:
            (left if x < split else right).append(i)

    def gap(pair):
        i, j = pair
        return sum((a - b) ** 2 for a, b in zip(positions[i], positions[j]))

    closest = min(((i, j) for i in left for j in right), key=gap)
    return [left, right], list(closest)


positions = load_observations("kmeans.csv")
num_clusters = 3

Both axes take one scale, never one each. The model measures a Euclidean distance to a center, and an axis stretched to fill the canvas would draw clusters that do not match the arithmetic underneath them.

draw_observations(positions, "Observations")

A baseline with scikit-learn

Scikit-learn ships a \(k\)-means implementation, so we have something to compare against.

def baseline_clusters(positions, num_clusters):
    from sklearn.cluster import KMeans

    kmeans = KMeans(n_clusters=num_clusters, random_state=0).fit(positions)
    centers = [tuple(center) for center in kmeans.cluster_centers_]
    total_distance = kmeans.inertia_
    return Clusters(total_distance, list(kmeans.labels_), centers)


baseline = baseline_clusters(positions, num_clusters)
print(f"Total squared distance (scikit-learn): {baseline.total_distance:.2f}")
Total squared distance (scikit-learn): 26243.63
draw_clusters(positions, baseline, "Scikit-learn solution")

Model implementation

import knitro

minimize_distance(positions, num_clusters, ...) builds the model, solves it with Knitro, and returns Clusters. The keyword arguments are the side constraints, all optional and all left empty for the first solve.

Both the objective \(\sum_i \sum_k a^k_i d^k_i\) and the barycenter constraint multiply two variables, and add_constraint takes them as written.

def minimize_distance(
    positions,
    num_clusters,
    *,
    must_link=(),
    cannot_link=(),
    min_size=None,
    max_size=None,
):
    n = len(positions)
    prob = knitro.Problem()
    observations, clusters = range(n), range(num_clusters)
    pairs = [(i, k) for i in observations for k in clusters]

    a = {pair: prob.add_variable(vtype=knitro.KN_VARTYPE_BINARY) for pair in pairs}
    d = {pair: prob.add_variable(lb=0) for pair in pairs}
    cx = {k: prob.add_variable() for k in clusters}
    cy = {k: prob.add_variable() for k in clusters}

    prob.add_objective(prob.nsum(a[pair] * d[pair] for pair in pairs))

    for i in observations:
        prob.add_constraint(prob.nsum(a[i, k] for k in clusters) == 1)

    for k in clusters:
        size = prob.nsum(a[i, k] for i in observations)
        total_x = prob.nsum(a[i, k] * positions[i][0] for i in observations)
        total_y = prob.nsum(a[i, k] * positions[i][1] for i in observations)
        prob.add_constraint(cx[k] * size == total_x)
        prob.add_constraint(cy[k] * size == total_y)

    for i in observations:
        x, y = positions[i]
        for k in clusters:
            prob.add_constraint(d[i, k] == (x - cx[k]) ** 2 + (y - cy[k]) ** 2)

    for group in must_link:
        for i in group[1:]:
            for k in clusters:
                prob.add_constraint(a[i, k] == a[group[0], k])
    for x in range(len(cannot_link)):
        for y in range(x + 1, len(cannot_link)):
            for k in clusters:
                prob.add_constraint(a[cannot_link[x], k] + a[cannot_link[y], k] <= 1)
    for k in clusters:
        if max_size is not None:
            prob.add_constraint(prob.nsum(a[i, k] for i in observations) <= max_size)
        if min_size is not None:
            prob.add_constraint(prob.nsum(a[i, k] for i in observations) >= min_size)

    prob.solve()

    labels = [max(clusters, key=lambda k: a[i, k].value) for i in observations]
    centers = [(cx[k].value, cy[k].value) for k in clusters]
    total_distance = prob.get_attr(knitro.KN_ATTR_OBJ_VALUE)
    return Clusters(total_distance, labels, centers)

Comparison with scikit-learn

Let’s solve the same instance with Knitro.

clusters = minimize_distance(positions, num_clusters)
=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

concurrent_evals         0
feastol                  1e-06
feastol_abs              1e-06
opttol                   1e-06
opttol_abs               0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
         The Knitro mixed integer solver is designed for convex problems.
         For non-convex problems it is only a heuristic, and the reported
         bounds and optimality claims cannot be verified.


Problem Characteristics                     |           Presolved
-----------------------
Problem type: MIQCQP
Objective: minimize / quadratic
Number of variables:                   1806 |                          1806
  bounds:         lower     upper     range |     lower     upper     range
                    900         0       900 |       900         0       900
                             free     fixed |                free     fixed
                                6         0 |                   6         0
                  cont.    binary   integer |     cont.    binary   integer
                    906       900         0 |       906       900         0
Number of expressions:                10829 |                             0
Number of constraints:                 1206 |                          1206
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:           300         0         0 |       300         0         0
  quadratic:        906         0         0 |       906         0         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0      5400           |         0      5400          
  quadratic:       1800      3606      2706 |      1800      3606      2706
  total:           1800      5406      2706 |      1800      5406      2706

Knitro using Branch and Bound method with 8 threads.

Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.

Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [9e-02, 4e+01] |                [9e-02, 8e+01]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 2e+03]

Root node relaxation
--------------------

 Iter      Objective      Feasibility        Optimality       Time 
                             error              error        (secs)
 ----      ---------      -----------        ----------      ------
    0        3695.50          1109.15          0.280016       0.113
    1        61427.1          1.46486           110.560       0.121
    2        61427.8      4.43006e-02           5.14429       0.129
    3        61420.7          2.21936           4.06169       0.135
    4        61418.2          1.98961           7.23102       0.143
    5        61410.1         0.502792       3.09268e-02       0.152
    6        61402.1         0.788797           4.11273       0.159
    7        61392.4          1.17870           4.30275       0.166
    8        61109.0          2.80582           65.0264       0.172
    9        61035.3          1.92984           57.0354       0.179
   10        60749.6         0.281936           34.1836       0.186
   11        60560.3         0.204399           86.0664       0.196
   12        60243.9          5.12656           250.824       0.205
   13        59967.6          2.79784          0.286995       0.214
   14        59319.6          16.8188           554.049       0.222
   15        58880.8          10.9504           553.930       0.228
   16        59068.7          3.97989           626.526       0.237
   17        59159.4         0.971826           811.995       0.245
   18        59176.5         0.590024           852.490       0.251
   19        59201.4         0.178497           914.084       0.260
   20        59173.5         0.709994           991.480       0.267
   30        57938.6          2.52692           1464.59       0.333
   40        55428.6          5.78292           1132.77       0.443
   50        48255.2         0.235207           756.542       0.522
   60        40811.6          1.76421           2048.20       0.574
   70        40591.5      2.91997e-03           41.2780       0.649
   80        39836.4      8.32420e-03           18.9218       0.740
   90        39241.5      2.42676e-03           30.5086       0.826
  100        38607.6      2.60814e-03           38.9491       0.910
  110        38071.8      1.27730e-03           8.86080       0.997
  120        37526.3      7.76244e-03           19.9145       1.085
  130        36959.6         0.145227           33.9095       1.172
  140        36408.4      3.97384e-02           33.7348       1.248
  150        35763.1      3.03368e-02           45.1469       1.340
  160        35390.9      4.28555e-04           51.7920       1.435
  170        34894.7         0.132845           47.3950       1.524
  180        34565.5      9.12914e-03           56.4190       1.623
  190        34043.2      2.20938e-02           57.4717       1.704
  200        33528.7      6.11846e-03           23.6403       1.782
  300        30904.2      7.29125e-02           28.4540       2.499
  400        29248.2      1.21879e-03           4.05594       3.210
  500        27987.3      1.74104e-02           98.3979       3.968
  600        26273.3      1.76612e-05           4.71021       4.579

Tree search
-----------

       Nodes        Best solution   Best bound      Gap       Time 
   Expl  |  Unexpl      value         value                  (secs)
   ---------------  -------------   ----------      ---      ------
      1       0      26240.7 LEAF      26240.7      0.00%     4.714

EXIT: Optimal solution found (assuming convexity).

HINT: The problem may be a non-convex mixed-integer problem.  Set
      mip_multistart=1 to enable a mixed-integer multistart heuristic,
      which may improve the chances of finding the global solution.

Final Statistics for MIP
------------------------
Final objective value               =  2.62406575903795310e+04
Final bound value                   =  2.62406575903795310e+04
Final optimality gap (abs / rel)    =  0.00000e+00 / 0.00000e+00 (0.00%)
# of root cutting plane rounds      =  0
# of restarts                       =  0
# of nodes processed                =  1 (4.702s)
# of strong branching evaluations   =  0 (0.000s)
# of function evaluations           =  0 (0.000s)
# of gradient evaluations           =  0 (0.000s)
# of hessian evaluations            =  0 (0.000s)
# of hessian-vector evaluations     =  0
# of subproblems processed          =  1 (4.702s)
Total program time (secs)           =  4.71431 (6.011 CPU time)
Time spent in evaluations (secs)    =  0.00000

Cuts statistics (gen / add)
---------------------------
Knapsack cuts                       =  0 / 0
Mixed-integer rounding cuts         =  0 / 0
Flow-cover cuts                     =  0 / 0
Probing cuts                        =  0 / 0

Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump                    =  0 / 0 / 0.000s
Rounding heuristic                  =  0 / 0 / 0.000s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  3 / 0 / 0.648s

===========================================================================
Show the full outputHide the full output
report_clusters(positions, clusters)
cluster  observations   centre x   centre y    distance  share
--------------------------------------------------------------
      1            82      20.10       4.89     6981.86  26.6%
      2           118       1.48       4.56     9569.17  36.5%
      3           100      11.15      21.20     9689.63  36.9%

Total squared distance 26240.66, cluster sizes 82 to 118
draw_clusters(positions, clusters, "Knitro solution")

The two totals agree to four significant figures, and Knitro takes a few times longer to get there. Neither is proved optimal. The model is nonconvex, so Knitro returns a local optimum, and scikit-learn restarts from several random initializations for the same reason. On the plain problem there is nothing to gain by switching. The rest of this page is about what the dedicated algorithm cannot do.

Cluster sizes

A second dataset, split into 10 clusters.

positions_2 = load_observations("kmeans2.csv")
draw_observations(positions_2, "Observations")

Left alone, the clusters come out very uneven:

unbalanced = minimize_distance(positions_2, 10)
=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

concurrent_evals         0
feastol                  1e-06
feastol_abs              1e-06
opttol                   1e-06
opttol_abs               0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
         The Knitro mixed integer solver is designed for convex problems.
         For non-convex problems it is only a heuristic, and the reported
         bounds and optimality claims cannot be verified.


Problem Characteristics                     |           Presolved
-----------------------
Problem type: MIQCQP
Objective: minimize / quadratic
Number of variables:                   1820 |                          1820
  bounds:         lower     upper     range |     lower     upper     range
                    900         0       900 |       900         0       900
                             free     fixed |                free     fixed
                               20         0 |                  20         0
                  cont.    binary   integer |     cont.    binary   integer
                    920       900         0 |       920       900         0
Number of expressions:                10262 |                             0
Number of constraints:                 1010 |                          1010
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            90         0         0 |        90         0         0
  quadratic:        920         0         0 |       920         0         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0      5400           |         0      5400          
  quadratic:       1800      3620      2720 |      1800      3620      2720
  total:           1800      5420      2720 |      1800      5420      2720

Knitro using Branch and Bound method with 8 threads.

Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.

Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [9e-02, 4e+01] |                [9e-02, 8e+01]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+00] |                [1e+00, 2e+03]

Root node relaxation
--------------------

 Iter      Objective      Feasibility        Optimality       Time 
                             error              error        (secs)
 ----      ---------      -----------        ----------      ------
    0        1353.86          1109.12          0.255056       0.022
    1        1408.83          1105.95           2.77395       0.029
    2        2652.19          1034.35           7.47763       0.034
    3        6845.10          776.437           11.0477       0.040
    4        7481.17          735.078           11.9023       0.045
    5        16440.7          113.768           26.4435       0.056
    6        17832.6          3.02890           19.4946       0.063
    7        17720.8          1.83506           117.490       0.072
    8        17380.8          8.73919           109.398       0.080
    9        17260.1          7.84965           97.2386       0.089
   10        17117.0         0.216981          0.881562       0.101
   11        17008.1          2.85782           6.93578       0.111
   12        16790.6          1.85313          0.693047       0.122
   13        16184.9          13.7071           23.8737       0.132
   14        15650.1          12.1496           25.2712       0.139
   15        15181.2         0.494600           28.1361       0.150
   16        14728.0          1.05001          0.783843       0.161
   17        14411.0         0.362034          0.924467       0.172
   18        14181.4         0.121771           1.06941       0.183
   19        13390.7         0.556190           66.0772       0.196
   20        13180.7         0.556865           24.4949       0.210
   30        12653.3         0.766207           16.0930       0.316
   40        12338.4         0.201530           4.75078       0.426
   50        11919.5         0.677553           36.7197       0.606
   60        11253.2          1.54308           70.0946       0.730
   70        10735.9         0.103027           6.92797       0.889
   80        10003.9      1.92900e-02           21.3788       1.009
   90        9136.95         0.372813           31.5789       1.159
  100        8425.64         0.267365           23.1098       1.345
  110        7825.28      3.01980e-03           30.4255       1.509
  120        7337.86      6.46777e-02           33.8307       1.615
  130        6880.39      7.13774e-02           62.5398       1.717
  140        6714.58      2.61975e-02           72.6677       1.864
  150        6389.23      4.47745e-02           18.5179       2.034
  160        6229.16      7.31402e-03           17.7292       2.215
  170        5962.68      1.14265e-02           24.1955       2.367
  180        5775.81      3.12063e-04           26.7308       2.467
  190        5495.03      2.34993e-02           5.66067       2.620
  200        5307.96      4.99778e-02           5.22593       2.798
  300        4252.44      1.59170e-04           3.38690       4.080
  400        3224.25         0.532993           32.7134       5.217
  500        3097.29         0.393481           7.80063       6.249

Tree search
-----------

       Nodes        Best solution   Best bound      Gap       Time 
   Expl  |  Unexpl      value         value                  (secs)
   ---------------  -------------   ----------      ---      ------
      1       0      2216.00 LEAF      2216.00      0.00%     6.516

EXIT: Optimal solution found (assuming convexity).

HINT: The problem may be a non-convex mixed-integer problem.  Set
      mip_multistart=1 to enable a mixed-integer multistart heuristic,
      which may improve the chances of finding the global solution.

Final Statistics for MIP
------------------------
Final objective value               =  2.21599548133198414e+03
Final bound value                   =  2.21599548133198414e+03
Final optimality gap (abs / rel)    =  0.00000e+00 / 0.00000e+00 (0.00%)
# of root cutting plane rounds      =  0
# of restarts                       =  0
# of nodes processed                =  1 (6.509s)
# of strong branching evaluations   =  0 (0.000s)
# of function evaluations           =  0 (0.000s)
# of gradient evaluations           =  0 (0.000s)
# of hessian evaluations            =  0 (0.000s)
# of hessian-vector evaluations     =  0
# of subproblems processed          =  1 (6.509s)
Total program time (secs)           =  6.51689 (7.148 CPU time)
Time spent in evaluations (secs)    =  0.00000

Cuts statistics (gen / add)
---------------------------
Knapsack cuts                       =  0 / 0
Mixed-integer rounding cuts         =  0 / 0
Flow-cover cuts                     =  0 / 0
Probing cuts                        =  0 / 0

Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump                    =  0 / 0 / 0.000s
Rounding heuristic                  =  0 / 0 / 0.000s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  4 / 0 / 0.598s

===========================================================================
Show the full outputHide the full output
report_clusters(positions_2, unbalanced)
cluster  observations   centre x   centre y    distance  share
--------------------------------------------------------------
      1            12       4.28      -3.10      363.83  16.4%
      2             9      21.96      -0.36      346.84  15.7%
      3            12      16.32      14.02      233.73  10.5%
      4            10       3.47      14.18      107.94   4.9%
      5             4      21.25      30.16      218.62   9.9%
      6             5      -9.11       3.37      206.05   9.3%
      7            11       2.29       6.02      184.16   8.3%
      8             8      11.95       7.85       84.15   3.8%
      9            11       7.57      23.11      150.44   6.8%
     10             8      28.42      12.19      320.23  14.5%

Total squared distance 2216.00, cluster sizes 4 to 12
draw_clusters(positions_2, unbalanced, "Clusters of various sizes")

Bounding every cluster between 8 and 12 observations evens them out:

balanced = minimize_distance(positions_2, 10, min_size=8, max_size=12)
=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

concurrent_evals         0
feastol                  1e-06
feastol_abs              1e-06
opttol                   1e-06
opttol_abs               0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
         The Knitro mixed integer solver is designed for convex problems.
         For non-convex problems it is only a heuristic, and the reported
         bounds and optimality claims cannot be verified.


Problem Characteristics                     |           Presolved
-----------------------
Problem type: MIQCQP
Objective: minimize / quadratic
Number of variables:                   1820 |                          1820
  bounds:         lower     upper     range |     lower     upper     range
                    900         0       900 |       900         0       900
                             free     fixed |                free     fixed
                               20         0 |                  20         0
                  cont.    binary   integer |     cont.    binary   integer
                    920       900         0 |       920       900         0
Number of expressions:                10262 |                             0
Number of constraints:                 1030 |                          1030
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            90        20         0 |        90        20         0
  quadratic:        920         0         0 |       920         0         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0      7200           |         0      7200          
  quadratic:       1800      3620      2720 |      1800      3620      2720
  total:           1800      7220      2720 |      1800      7220      2720

Knitro using Branch and Bound method with 8 threads.

Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.

Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [9e-02, 4e+01] |                [9e-02, 8e+01]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 1e+01] |                [1e+00, 2e+03]

Root node relaxation
--------------------

 Iter      Objective      Feasibility        Optimality       Time 
                             error              error        (secs)
 ----      ---------      -----------        ----------      ------
    0        1290.90          1109.12          0.259250       0.034
    1        1317.08          1107.51           1.34656       0.048
    2        4567.51          905.105           12.4905       0.060
    3        8146.36          669.513           26.3210       0.071
    4        9577.93          570.390           32.6265       0.082
    5        12411.4          384.825           15.7704       0.101
    6        17787.7         0.186804           7.79838       0.111
    7        17669.2         0.808774           3.70947       0.121
    8        17517.9         0.724213           5.80388       0.134
    9        17392.2         0.335118          0.305357       0.150
   10        17290.1         0.230633           3.28107       0.163
   11        17074.1         0.935114           4.89089       0.172
   12        16928.1         0.710185           5.50899       0.185
   13        16781.0         0.761296           10.8371       0.198
   14        16247.8          2.47114           5.63185       0.207
   15        15609.7          3.57355           8.46401       0.221
   16        14968.3          4.25448           19.1055       0.233
   17        13993.0          5.08790           33.5284       0.247
   18        13117.5          5.26773           27.2724       0.260
   19        12375.7          5.00643           38.8465       0.273
   20        12165.2          4.83416           38.8589       0.286
   30        9968.46          4.43857           28.0790       0.419
   40        8717.25          3.12368           17.4706       0.547
   50        7777.59          1.79943           14.5063       0.691
   60        6765.33         0.869958           7.20423       0.853
   70        6059.06      3.61800e-03           8.77838       1.026
   80        5515.16      2.08992e-04           14.2629       1.165
   90        5123.48      1.00147e-02           15.8843       1.295
  100        4537.94      6.40318e-03           18.6595       1.430
  110        4198.10         0.182432           9.07008       1.567
  120        3838.68      1.55369e-02           7.01654       1.700
  130        3639.15      6.25700e-02           12.3614       1.836
  140        3367.10         0.375746           18.5609       1.971
  150        3144.98         0.390591           9.83157       2.108
  160        3040.63         0.464382           14.9895       2.261
  170        2952.00         0.243887           26.3686       2.383
  180        2856.68         0.185019           25.5800       2.526
  190        2679.96      7.19675e-04           19.9252       2.658
  200        2502.63         0.249926           9.08155       2.794

Tree search
-----------

       Nodes        Best solution   Best bound      Gap       Time 
   Expl  |  Unexpl      value         value                  (secs)
   ---------------  -------------   ----------      ---      ------
      1       2      2440.50 FCRD         -inf                2.943
      3       0      2440.50           2440.50     -0.00%     3.802

EXIT: Optimal solution found (assuming convexity).

HINT: The problem may be a non-convex mixed-integer problem.  Set
      mip_multistart=1 to enable a mixed-integer multistart heuristic,
      which may improve the chances of finding the global solution.

Final Statistics for MIP
------------------------
Final objective value               =  2.44050077279231937e+03
Final bound value                   =  2.44050086853016546e+03
Final optimality gap (abs / rel)    =  -9.57378e-05 / -3.92288e-08 (-0.00%)
# of root cutting plane rounds      =  1
# of restarts                       =  0
# of nodes processed                =  3 (2.925s)
# of strong branching evaluations   =  0 (0.000s)
# of function evaluations           =  0 (0.000s)
# of gradient evaluations           =  0 (0.000s)
# of hessian evaluations            =  0 (0.000s)
# of hessian-vector evaluations     =  0
# of subproblems processed          =  2 (2.925s)
Total program time (secs)           =  3.80384 (5.225 CPU time)
Time spent in evaluations (secs)    =  0.00000

Cuts statistics (gen / add)
---------------------------
Knapsack cuts                       =  0 / 0
Mixed-integer rounding cuts         =  0 / 0
Flow-cover cuts                     =  0 / 0
Probing cuts                        =  0 / 0

Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump                    =  0 / 0 / 0.000s
Rounding heuristic                  =  1 / 1 / 0.003s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  6 / 0 / 1.842s

===========================================================================
Show the full outputHide the full output
report_clusters(positions_2, balanced)
cluster  observations   centre x   centre y    distance  share
--------------------------------------------------------------
      1             8      27.57      10.22      240.86   9.9%
      2             8      20.86      -2.12      332.24  13.6%
      3             8      18.04      10.06      108.93   4.5%
      4             8      11.20      14.07       62.54   2.6%
      5            10       7.08      23.17      123.30   5.1%
      6             9       2.12      13.74       91.66   3.8%
      7            11       1.94      -2.79      253.39  10.4%
      8             8      -5.82       5.25      403.75  16.5%
      9             8      20.62      25.25      608.19  24.9%
     10            12       8.29       5.11      215.64   8.8%

Total squared distance 2440.50, cluster sizes 8 to 12
draw_clusters(positions_2, balanced, "Clusters with restricted sizes")

Balancing is not free. The bounds pull observations away from their nearest center, so the best balanced clustering costs more than the best unbalanced one. The two solves above need not show that, and may even come out the other way round, since each stops at a local optimum and neither is the best of its problem. What matters is that the requirement can be stated at all, and that is the reason to write the problem out rather than call a library.