如何在R中使用Clarabel求解带二次约束的线性规划

正确转换步骤 #

1. 正定约束等价转换：由于Q是对称正定矩阵（可通过Cholesky分解验证），原约束sqrt((x-y) %*% Q %*% (x-y)) <= 0.1等价于||S(x-y)||₂ ≤ 0.1，其中S = chol(Q)是Q的Cholesky分解矩阵。
2. 匹配SOCP标准形式：Clarabel的二阶锥（q型锥）要求向量第一个分量≥其余分量的欧几里得范数。这里取v=0.1，u=S(x-y)，即[0.1; S(x-y)]需属于维度为11的二阶锥。
3. 构造约束矩阵与向量：
   • 在原有约束矩阵A后添加一行全0向量（对应v的线性项），再添加-S矩阵（对应u的线性项）。
   • 在原有约束向量b后添加0.1和S%*%y（对应约束的常数项）。

完整可运行代码 #

library(clarabel)
library(Matrix)

# 原始线性规划部分
a = c(1.425, 0.27, -0.085, -0.733, -0.534, 1.402, -0.199, -0.875, -1.586, -2.126)
xl = rep(0, 10)
xu = rep(1, 10)
A = rbind(rep(1, 10), Diagonal(10), -Diagonal(10))
b = c(1, xu, -xl)

# 二次约束测试数据
y = c(0.242, 0.058, 0.204, 0.112, 0.123, 0.143, 0, 0.065, 0, 0.053)
Q = structure(c(1.291, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.018, 0.525, 0,
                0, 0, 0, 0, 0, 0, 0, 0.034, 0.02, 0.543, 0, 0, 0, 0, 0, 0, 0,
                0.012, 0.011, 0.015, 0.536, 0, 0, 0, 0, 0, 0, 0.039, 0.021, 0.036,
                0.012, 0.874, 0, 0, 0, 0, 0, 0.064, 0.02, 0.046, 0.009, 0.054,
                0.735, 0, 0, 0, 0, 0.02, 0.007, 0.019, 0.001, 0.025, 0.035, 1.791,
                0, 0, 0, 0.031, 0.018, 0.029, 0.011, 0.032, 0.043, 0.03, 1.242,
                0, 0, 0.008, 0.008, 0.012, 0.045, 0.006, 0.004, 0.008, 0.007,
                0.55, 0, 0.005, 0.007, 0.009, 0.041, 0.003, 0, 0.01, 0.006, 0.062,
                0.479), .Dim = c(10L, 10L))
Q = as(as(Q + lower.tri(Q) * t(Q), 'CsparseMatrix'), 'symmetricMatrix')

# 转换二次约束为SOCP形式
S = chol(Q)
Sy = drop(S %*% y)

# 构造新的约束矩阵和向量
A_new = rbind(A, matrix(0, nrow = 1, ncol = 10), -S)
b_new = c(b, 0.1, Sy)

# 调用Clarabel求解
res = clarabel(A_new, b_new, -a, cones = list(z=1, l=20, q=11))

# 查看结果
res$x
solver_status_descriptions()[res$status]

错误尝试分析 #

你的尝试核心问题是约束构造不符合Clarabel的SOCP标准格式：

• 错误处理了二次约束的常数项与线性项对应关系，导致约束矩阵和向量的组合无法匹配二阶锥要求。
• 引入的gamma和bb计算逻辑偏离了SOCP约束的标准映射方式，导致约束失效或不可行。

原始问题描述 #

我需要将二次约束转换为Clarabel优化包（二阶锥规划求解器）所需的形式。首先给出问题的线性部分代码：

library(clarabel)
library(Matrix)

a = c(1.425, 0.27, -0.085, -0.733, -0.534, 1.402, -0.199, -0.875, -1.586, -2.126)
xl = rep(0, 10)
xu = rep(1, 10)
A = rbind(rep(1, 10), Diagonal(10), -Diagonal(10))
b = c(1, xu, -xl)

x = clarabel(A, b, -a, cones=list(z=1, l=20))$x

以下是二次约束的测试数据：

y = c(0.242, 0.058, 0.204, 0.112, 0.123, 0.143, 0, 0.065, 0, 0.053)
Q = structure(c(1.291, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.018, 0.525, 0,
                0, 0, 0, 0, 0, 0, 0, 0.034, 0.02, 0.543, 0, 0, 0, 0, 0, 0, 0,
                0.012, 0.011, 0.015, 0.536, 0, 0, 0, 0, 0, 0, 0.039, 0.021, 0.036,
                0.012, 0.874, 0, 0, 0, 0, 0, 0.064, 0.02, 0.046, 0.009, 0.054,
                0.735, 0, 0, 0, 0, 0.02, 0.007, 0.019, 0.001, 0.025, 0.035, 1.791,
                0, 0, 0, 0.031, 0.018, 0.029, 0.011, 0.032, 0.043, 0.03, 1.242,
                0, 0, 0.008, 0.008, 0.012, 0.045, 0.006, 0.004, 0.008, 0.007,
                0.55, 0, 0.005, 0.007, 0.009, 0.041, 0.003, 0, 0.01, 0.006, 0.062,
                0.479), .Dim = c(10L, 10L))
Q = as(as(Q + lower.tri(Q) * t(Q), 'CsparseMatrix'), 'symmetricMatrix')

我要施加的二次约束为sqrt((x-y) %*% Q %*% (x-y)) <= 0.1，等价于x %*% Q %*% x - 2 * x %*% Q %*% y - (0.1^2 - y %*% Q %*% y) <= 0。

通过将二次项放入目标函数，我做了简单搜索得到近似解作为参考：

lambda = 13.8703
x = clarabel(A, b, -a + lambda * -2 * drop(Q %*% y), 2 * lambda * Q, cones=list(z=1, l=20))$x

但我尝试参考各类将QCQP转换为SOCP的指南时，要么得到不可行错误，要么约束被完全忽略。以下是我的尝试（均不可行）：

S = chol(Q)
lin = drop(-2 * Q %*% y)
rhs = drop(0.1^2 - y %*% Q %*% y)
bb = -drop(solve(S, lin))
gamma = -sqrt(drop(crossprod(bb)) + rhs)

res = clarabel(rbind(A, 0, S), c(b, rhs, lin), -a, cones=list(z=1, l=20, q=11))
solver_status_descriptions()[res$status]

res = clarabel(rbind(A, 0, S), c(b, gamma, bb), -a, cones=list(z=1, l=20, q=11))
solver_status_descriptions()[res$status]

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