如何在R中构建带二进制变量的LP（线性规划）聚类问题

核心问题处理逻辑 #

你当前遇到的核心难点是目标函数中的簇质心x̄ⱼ属于连续决策变量，直接代入原始模型会产生非线性项，无法通过仅支持混合整数线性规划（MILP）的lpSolve直接求解，可按两种方案处理：

• 方案1：将问题转化为带容量约束的K中心点聚类（Capacitated K-Medoids），强制簇质心为现有样本点中的一个，消除x̄ⱼ这类连续变量，直接适配lpSolve的求解能力
• 方案2：使用支持混合整数非线性规划（MINLP）的R工具包直接对原始K均值模型建模

样本数据生成 #

你提供的样本数据生成代码如下：

set.seed(123)
id<- seq(1:50)
lon <- rnorm(50, 88.5, 0.125)
lat <- rnorm(50, 22.4, 0.15)
demand <- round(runif(50, min=20, max=40))
df<- data.frame(id, lon, lat, demand)

方案1：K中心点线性建模（适配lpSolve） #

转化逻辑：引入新的二进制变量z_kj，当第k个样本点被选为第j个簇的中心时取值为1，否则为0。此时簇j的质心坐标等于第k个点的已知坐标，目标函数转化为最小化所有点到所属簇中心点的距离之和，约束条件调整为：

1. 每个点仅属于1个簇：∑j y_ij = 1 对所有i成立
2. 每个簇仅有1个中心：∑k z_kj =1 对所有j成立
3. 点仅可分配到存在中心的簇：y_ij ≤ ∑k z_kj 对所有i,j成立
4. 容量约束：∑i y_ij * q_i ≤ Q_j 对所有j成立
   该方案所有约束和目标函数均为线性，可直接用lpSolve构建求解，适合中小规模样本。

方案2：原始K均值MINLP实现 #

如果需要严格贴合原始问题模型，可使用ompr包搭配ROI求解器直接建模，无需额外转化变量，示例代码如下：

# 安装依赖包（首次运行执行）
# install.packages(c("ompr", "ompr.roi", "ROI", "ROI.plugin.ipop"))
library(ompr)
library(ompr.roi)
library(ROI)

n <- nrow(df)
k <- 3 # 自定义簇数
max_load <- 700 # 自定义单簇最大负载Q_j

# 构建模型
cap_kmeans_model <- MIPModel() |>
  # 二进制变量y[i,j]：点i属于簇j时取1，否则取0
  add_variable(y[i,j], i = 1:n, j = 1:k, type = "binary") |>
  # 连续变量x[j,d]：簇j在d维度的质心坐标，d=1为经度，d=2为纬度
  add_variable(x[j,d], j = 1:k, d = 1:2, type = "continuous") |>
  # 目标函数：最小化所有点到所属簇质心的距离平方和
  set_objective(
    sum_over(
      (df$lon[i] - x[j,1])^2 + (df$lat[i] - x[j,2])^2 * y[i,j], 
      i = 1:n, j = 1:k
    ), 
    "min"
  ) |>
  # 约束1：每个点仅属于1个簇
  add_constraint(sum_over(y[i,j], j = 1:k) == 1, i = 1:n) |>
  # 约束2：质心为簇内点的坐标均值
  add_constraint(
    x[j,1] * sum_over(y[i,j], i = 1:n) == sum_over(df$lon[i] * y[i,j], i = 1:n), 
    j = 1:k
  ) |>
  add_constraint(
    x[j,2] * sum_over(y[i,j], i = 1:n) == sum_over(df$lat[i] * y[i,j], i = 1:n), 
    j = 1:k
  ) |>
  # 约束3：簇容量不超过上限
  add_constraint(sum_over(df$demand[i] * y[i,j], i = 1:n) <= max_load, j = 1:k)

# 求解模型
solved_result <- solve_model(cap_kmeans_model, with_ROI(solver = "ipop"))

内容的提问来源于stack exchange，提问作者Shibaprasad