使用osqp求解简单二次规划问题

一、问题描述

已知：
$$min(x_1-1{)}^2+(x_2-2{)}^{2}s.t.0⩽x_1⩽1.5,1⩽x_2⩽2.5$$
目标函数为二元二次函数，可行域为线性、凸集，此为二次规划问题，可将其转换成二次规划表达式再进行求解。相关数学概念参考另一篇： 最优化问题基础理论概述。

二、数学推导

1. 目标函数处理

$$f(x_1,x_2)=(x_1-1{)}^2+(x_2-2{)}^2=x_1^2+x_2^2-2x_1-4x_2+C$$

其中，常数项用$C$表示；

令，$X=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$，则
$$\begin{array}{rl}f(x_1,x_2)& =[x_1,x_2][x_1,x_2{]}^T+[-2,-4][x_1,x_2{]}^T\\ & =X^{T}X+[-2,-4]X\\ & =\frac{1}{2}{X}^T\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right]X+{\left[\begin{array}{c}-2\\ -4\end{array}\right]}^{T}X\\ & =\frac{1}{2}{X}^{T}PX+Q^{T}X\end{array}$$

其中，$P=\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right]，Q=\left[\begin{array}{c}-2\\ -4\end{array}\right]$
 
关于为什么要写成 $\frac{1}{2}{X}^{T}PX$ 形式，因为此时 $P$ 为目标函数的海塞矩阵，具体参看 此链接。

2. 约束条件处理

$$\begin{array}{r}\{\begin{array}{c}0⩽x_1⩽1.5\\ 1⩽x_2⩽2.5\end{array}⟺\left[\begin{array}{c}0\\ 1\end{array}\right]⩽\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]⩽\left[\begin{array}{c}1.5\\ 2.5\end{array}\right]⟺\left[\begin{array}{c}0\\ 1\end{array}\right]⩽\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]X⩽\left[\begin{array}{c}1.5\\ 2.5\end{array}\right]\end{array}$$
令$L_B=\left[\begin{array}{c}0\\ 1\end{array}\right],A=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],U_B=\left[\begin{array}{c}1.5\\ 2.5\end{array}\right]$ ，整理得约束条件如下：
$$L_B⩽AX⩽U_B$$

三、代码编写

  由步骤 二、数学推导 得到5个矩阵：

• $P$ : 二次型矩阵(实对称矩阵)；
• $Q$ : 一次项矩阵；
• $U_B$ : 上边界矩阵；
• $L_B$ : 下边界矩阵；
• $A$ : 边界系数矩阵；

  现在根据这5个矩阵进行代码编写，是使用osqp进行二次型规划问题构建及求解。

代码如下：

Eigen::SparseMatrix<double> P(2, 2); // P， 二次型矩阵
Eigen::VectorXd Q(2);                // Q， 一次项向量
Eigen::SparseMatrix<double> A(2, 2); // 单位阵
Eigen::VectorXd lowerBound(2);       // 下边界向量
Eigen::VectorXd upperBound(2);       // 上边界向量

P.insert(0, 0) = 2.0;
P.insert(1, 1) = 2.0;
std::cout << "\033[34m" << "P:" << std::endl
            << P << "\033[0m" << std::endl;

A.insert(0, 0) = 1.0;
A.insert(1, 1) = 1.0;
std::cout << "\033[34m" << "A:" << std::endl
            << A << "\033[0m" << std::endl;

Q << -2, -4;
std::cout << "\033[34m" << "Q:" << std::endl
            << Q << "\033[0m" << std::endl;

lowerBound << 0.0, 1.0;
upperBound << 1.5, 2.5;

// Step 1： 创建求解器
OsqpEigen::Solver solver;
// Step 2： 设置(提升求解速度)
solver.settings()->setVerbosity(false);
solver.settings()->setWarmStart(true);

// Step 3： 初始化(7部分)
solver.data()->setNumberOfVariables(2);   // 变量数
solver.data()->setNumberOfConstraints(2); // 约束数
if (!solver.data()->setHessianMatrix(P))  // 海塞矩阵
{
    return;
}
if (!solver.data()->setGradient(Q)) // Q矩阵
{
    return;
}
if (!solver.data()->setLinearConstraintsMatrix(A)) // 线性约束矩阵A
{
    return;
}
if (!solver.data()->setLowerBound(lowerBound)) // 下边界矩阵
{
    return;
}
if (!solver.data()->setUpperBound(upperBound)) // 上边界矩阵
{
    return;
}

if (!solver.initSolver())
{
    return;
}

// Step 4：求解
Eigen::VectorXd QPSolution;
if (solver.solveProblem() != OsqpEigen::ErrorExitFlag::NoError)
{
    return;
}
QPSolution = solver.getSolution();
std::cout << "\033[1;32m" << "QPSolution:" << std::endl
            << QPSolution << "\033[0m" << std::endl;

运行结果如下:

可见，当$x_1=1,x_2=2$ 时目标函数取得最小。