基于R的加权矩阵损伤和最小化约束优化问题求助

看起来你的问题核心出在约束函数的错误定义，以及没有意识到这个问题其实可以转化为更易求解的线性规划问题。让我们一步步拆解并解决：

1. 约束函数的错误分析 #

你原本的约束2代码逻辑有误：

damage_per_country <- rowSums(D_weighted)/sum(D_weighted)
c2 = damage_per_country/sum(D_weighted)

这里的c2计算的是每个国家的损伤占比再除以全局总损伤，这完全偏离了“各国接收损伤均等分布”的需求。正确的约束应该是：每个国家的损伤占全局总损伤的比例等于1/N（N是国家数量），也就是直接让damage_per_country等于targ_dist，而不需要再除以总损伤。

此外，你的约束集合存在冗余：sum(damage_per_country)恒等于1，而targ_dist的和也是1，这会导致solnp在处理线性化问题时出现矛盾，进而提示无可行解。

2. 修正非线性优化（solnp）的代码 #

如果坚持使用Rsolnp，我们需要修正约束函数，同时去掉冗余的约束（只保留N-1个均等约束，因为所有比例相等的话，最后一个比例可以由前N-1个推导出来）：

library(Rsolnp)

# 假设D是你的对称方阵，先定义一个示例D方便测试
N <- 3
D <- matrix(c(1,0.2,0.1,
              0.2,1,0.3,
              0.1,0.3,1), nrow=N, byrow=TRUE)

# 目标函数：最小化全局损伤（逻辑正确，无需修改）
damage <- function(weights) {
  D_weighted <- t(t(D)*weights)
  return(sum(D_weighted))
}

# 修正后的约束函数
constr <- function(weights) {
  D_weighted <- t(t(D)*weights)
  damage_per_country <- rowSums(D_weighted)
  
  # 约束1：权重和为1
  c1 <- sum(weights)
  # 约束2：第2到N国的损伤 = 第1国的损伤（避免冗余）
  c2 <- damage_per_country[2:N] - damage_per_country[1]
  
  return(c(c1, c2))
}

# 目标约束值：权重和=1，其余约束差值=0
eqB <- c(1, rep(0, N-1))

# 初始权重
startweights <- rep(1/N, N)

# 执行优化
opt_weights <- solnp(
  pars = startweights,
  fun = damage,
  eqfun = constr,
  eqB = eqB,
  LB = rep(0, N),
  UB = rep(1, N),
  control=list(outer.iter=1000,trace=1, tol= 1e-6)
)

# 查看结果
opt_weights$pars

这里把“损伤占比均等”转化为“所有国家的绝对损伤值相等”，这样约束变成了线性的差值为0，避免了冗余，也简化了非线性程度。

3. 更高效的线性规划方法 #

其实你的问题本质是线性规划问题：

• 目标函数sum(D_weighted)是权重向量的线性组合（等于t(weights) %*% colSums(D)，因为sum(D_weighted) = sum_j weights[j] * sum_i D[i,j]）
• 所有约束都是线性的（权重和为1，各国损伤相等，权重非负）

我们可以用lpSolve包来求解，速度更快且更稳定：

library(lpSolve)

# 目标函数系数：每个权重对应的系数是D的列和（因为sum(D_weighted) = sum_j w_j * sum_i D[i,j]）
obj_coef <- colSums(D)

# 约束矩阵：
# 第一行：sum(w) = 1 → 所有元素为1
# 第2到N行：sum_j (D[i,j] - D[1,j])w_j = 0 → 对应位置为D[i,j]-D[1,j]
const_mat <- rbind(
  rep(1, N),
  t(D[2:N,] - D[1,])
)

# 约束方向：所有约束都是等于（"="）
const_dir <- rep("=", nrow(const_mat))

# 约束值：第一行是1，其余是0
const_val <- c(1, rep(0, N-1))

# 求解线性规划
lp_result <- lp(
  direction = "min",
  objective.in = obj_coef,
  const.mat = const_mat,
  const.dir = const_dir,
  const.rhs = const_val,
  all.int = FALSE,
  lower.tail = rep(0, N),
  upper.tail = rep(1, N)
)

# 查看最优权重
lp_result$solution

这种方法完全避开了非线性优化的收敛问题，因为线性规划的求解器处理这类问题非常成熟。

为什么你的原始代码不收敛？ #

1. 约束函数的计算逻辑错误，导致约束条件完全不符合需求；
2. 冗余的约束导致solnp在构建线性化问题时无法找到可行解；
3. 不必要地使用非线性约束处理本可以线性化的问题，增加了求解难度。

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