R-优化X，使所有n的minsumX_i最佳，其中X_n + X_n-1+ X_n-2> = Y_n，其中所有n的Y均已知

问题描述

我有一个长度为50的未知向量X和一个常数的已知向量Y（长度为50）。

我希望找到X，使得对于X_i> = 0，sum（X_i）被最小化，约束为：

X_n + X_ {n-1}> = Y_n

我不确定从R开头。

解决方法

我想您可以尝试CVXR来解决优化问题。

1. 首先，让我们如下定义一个矩阵M

M <- matrix(0,nrow = 10,ncol = 11)
for (i in 1:nrow(M)) {
  for (j in 1:ncol(M)) {
    if (j %in% (i+(0:1))) M[i,j] <- 1
  }
}

看起来像

> M
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]    1    1    0    0    0    0    0    0    0     0     0
 [2,]    0    1    1    0    0    0    0    0    0     0     0
 [3,]    0    0    1    1    0    0    0    0    0     0     0
 [4,]    0    0    0    1    1    0    0    0    0     0     0
 [5,]    0    0    0    0    1    1    0    0    0     0     0
 [6,]    0    0    0    0    0    1    1    0    0     0     0
 [7,]    0    0    0    0    0    0    1    1    0     0     0
 [8,]    0    0    0    0    0    0    0    1    1     0     0
 [9,]    0    0    0    0    0    0    0    0    1     1     0
[10,]    0    0    0    0    0    0    0    0    0     1     1

1. 然后，我们构造目标函数和约束

library(CVXR)
X <- Variable(11)
objective <- Minimize(sum(X))
constraints <- list( X>=0,M%*%X >= Y)
problem <- Problem(objective,constraints)
res <- solve(problem)

1. 最后，我们可以通过X看到res$getValue(X)的值

---

示例 如下给出Y

set.seed(1)
Y <- runif(10)

我们可以得到

Xopt <- res$getValue(X)

> Xopt
              [,1]
 [1,] 1.667850e-07
 [2,] 3.072356e-01
 [3,] 6.488860e-02
 [4,] 6.214644e-01
 [5,] 2.867441e-01
 [6,] 1.486883e-02
 [7,] 8.835218e-01
 [8,] 7.476340e-02
 [9,] 5.860353e-01
[10,] 5.264372e-02
[11,] 9.142897e-03

---

另一个可能的选项可能是pracma::fmincon，例如

pracma:: fmincon(rep(0,11),function(x) sum(x),A = -M,b = -Y,lb = 0,)

,

这可以表示为线性程序

min 1'x
s.t. Ax >= y
     x >= 0

我们可以使用lpSolve包来解决此问题。令n为x的长度，k为每个约束中x的数量。问题主体中定义的约束条件对应于k = 2，但问题的主题定义的约束条件不同，对应于k = 3。将有n-k + 1个约束。

embed(1:n,k)创建一个n-k + 1×k矩阵，每个矩阵的行是A对应行中矩阵的列的索引。例如，对于k = 2，{{1的第一行}}为1：2，因为embed(...)第一行的元素1和2为1，而A的其余行为零。 A输出的第二行是2：3，第三行是3：4，依此类推。然后，我们在各行上应用replace，以将n的零向量embed替换为这些职位。套用作品的方式是，它提供所需内容的转置，因此我们将其转回以获得numeric(n)。

最后，我们运行线性程序。我们可以使用A来检查返回的输出组件。特别是，我们将解决方案显示为str(out)。

out$solution