{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# K-Means 聚类问题 (MSSC)\n", "\n", "**类别：** 选址\n", "\n", "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n", "\n", "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/k-means-clustering-problem-mssc)。\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 问题描述\n", "\n", "**在 K-Means 聚类问题**（最小平方和聚类，MSSC）中，我们希望将一组多维观测点划分为 k 个聚类，每个聚类由其重心定义。每个观测点属于具有最近重心的聚类。更多细节，请参阅 [Wikipedia](http://en.wikipedia.org/wiki/K-means_clustering)。\n", "\n", "### 学习要点\n", "\n", "- 添加 [集合决策变量](https://optagent.pages.dev/guide/modeling/) 来建模聚类\n", "- 使用 [符号 Lambda 函数](https://optagent.pages.dev/guide/modeling/) 来计算每个聚类的重心和方差\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 数据\n", "\n", "数据格式如下：\n", "\n", "- 第 1 行：观测点数和每个观测点的维度数\n", "- 对每个观测点：每个维度上的坐标以及其在最优解中所属的聚类\n" ] }, { "cell_type": "markdown", "id": "f972b430", "metadata": {}, "source": [ "## 建模思路\n", "\n", "K-Means 聚类问题 (MSSC) 的 OptAgent 模型使用 [集合决策变量](https://optagent.pages.dev/guide/modeling/)。每个 set variable 表示一个聚类，其内部的元素表示属于该聚类的观测点。我们使用一个 [**partition** operator](https://optagent.pages.dev/guide/modeling/) 来确保每个观测点恰好属于一个聚类。\n", "\n", "我们使用 [符号 Lambda 函数](https://optagent.pages.dev/guide/modeling/) 对每个聚类内所有观测点在所有维度上的坐标应用 **sum** 算子来计算每个聚类的重心。注意，此求和中项的数量在搜索过程中会随着集合的大小变化而变化。\n", "\n", "然后我们可以计算总方差。一个聚类的方差是该聚类重心与每个观测点之间的欧几里得距离平方之和。与重心类似，我们使用 符号 Lambda 函数 计算每个聚类的方差。目标是最小化这些方差之和。\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python 实现\n" ] }, { "cell_type": "code", "execution_count": null, "id": "720c093d", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "from optagent import OptModel, solve\n", "\n", "def read_elem(filename):\n", " return Path(filename).read_text(encoding=\"utf-8\").split()\n", "\n", "#\n", "# Read instance data\n", "#\n", "def read_instance(filename):\n", " file_it = iter(read_elem(filename))\n", "\n", " # Data properties\n", " nb_observations = int(next(file_it))\n", " nb_dimensions = int(next(file_it))\n", "\n", " coordinates_data = [None] * nb_observations\n", " for o in range(nb_observations):\n", " coordinates_data[o] = [None] * (nb_dimensions)\n", " for d in range(nb_dimensions):\n", " coordinates_data[o][d] = float(next(file_it))\n", " next(file_it) # skip initial clusters\n", "\n", " return nb_observations, nb_dimensions, coordinates_data\n", "\n", "def main(instance_file, output_file=None, time_limit=60, k=2):\n", " nb_observations, nb_dimensions, coordinates_data = read_instance(instance_file)\n", "\n", " model = OptModel()\n", "\n", " # clusters[c] represents the observations assigned to cluster c.\n", " clusters = [model.set(nb_observations) for c in range(k)]\n", "\n", " # Each observation must belong to exactly one cluster.\n", " model.constraint(model.partition(clusters))\n", "\n", " coordinates = model.array(coordinates_data)\n", " variances = []\n", " for cluster in clusters:\n", " size = model.count(cluster)\n", "\n", " # Keep the empty-cluster guard from the Hexaly model.\n", " centroid = []\n", " for d in range(nb_dimensions):\n", " coordinate_lambda = model.lambda_function(\n", " lambda i: model.at(coordinates, i // 1, d)\n", " )\n", " centroid.append(\n", " model.iif(size == 0, 0.0, model.sum(cluster, coordinate_lambda) / size)\n", " )\n", "\n", " dimension_variances = []\n", " for d in range(nb_dimensions):\n", " centroid_d = centroid[d]\n", " dimension_variance_lambda = model.lambda_function(\n", " lambda i: model.pow(\n", " model.at(coordinates, i // 1, d) - centroid_d, 2\n", " )\n", " )\n", " dimension_variances.append(model.sum(cluster, dimension_variance_lambda))\n", " variances.append(model.sum(dimension_variances))\n", "\n", " # Minimize the total within-cluster sum of squares.\n", " total_variance = model.sum(variances)\n", " model.minimize(total_variance)\n", "\n", " solution = solve(model, time_limit_s=float(time_limit))\n", " print(f\"Observations = {nb_observations}; Dimensions = {nb_dimensions}; \"\n", " f\"Clusters = {k}; Total variance = {total_variance.value}; \"\n", " f\"Status = {solution.feasible}\")\n", "\n", " if solution.feasible:\n", " for c, cluster in enumerate(clusters):\n", " print(f\"Cluster {c}: {sorted(cluster.value)}\")\n", " if output_file is not None:\n", " lines = [str(total_variance.value), str(k)]\n", " lines.extend(\" \".join(str(o) for o in sorted(cluster.value)) for cluster in clusters)\n", " Path(output_file).write_text(\"\\n\".join(lines) + \"\\n\", encoding=\"utf-8\")\n", " return solution\n" ] }, { "cell_type": "markdown", "id": "4acdd230", "metadata": {}, "source": [ "## 本地运行\n", "\n", "以下代码格演示如何调用 OptAgent 模型。`k` 使用数据集中的类别数；如需更长搜索，可增大 `time_limit`。\n" ] }, { "cell_type": "code", "execution_count": null, "id": "92e0a449", "metadata": {}, "outputs": [], "source": [ "INSTANCE_DIR = Path.cwd() / \"instances\"\n", "print(\"Instances:\", INSTANCE_DIR)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5e692c97", "metadata": {}, "outputs": [], "source": [ "solution_ruspini = main(INSTANCE_DIR / \"ruspini.dat\", k=4, time_limit=1)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3bfec7f7", "metadata": {}, "outputs": [], "source": [ "solution_iris = main(INSTANCE_DIR / \"iris.dat\", k=3, time_limit=1)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7c6cc9a2", "metadata": {}, "outputs": [], "source": [ "solution_glass = main(INSTANCE_DIR / \"glass.dat\", k=6, time_limit=1)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1ea36a02", "metadata": {}, "outputs": [], "source": [ "solution_segment = main(INSTANCE_DIR / \"segment.dat\", k=7, time_limit=1)\n" ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.6" } }, "nbformat": 4, "nbformat_minor": 5 }