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

文章目录

一、问题描述
二、数学推导
- 1. 目标函数处理
- 2. 约束条件处理
三、代码编写

一、问题描述

已知：
$min(x_1-1)^2+(x_2-2)^2 \qquad s.t. \ \ 0 \leqslant x_1 \leqslant 1.5,\ \ 1 \leqslant x_2 \leqslant 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{matrix} x_1 \\ x_2 \end{matrix} \right]$ ，则
$\begin{aligned} f(x_1, x_2) &= [x_1, x_2][x_1, x_2]^T+[-2, -4][x_1, x_2]^T \\[2ex] &= X^TX+[-2, -4]X \\[2ex] &=\frac{1}{2} X^T \left[\begin{matrix} 2 &0 \\ 0&2 \end{matrix} \right] X+\left[\begin{matrix} -2 \\ -4 \end{matrix} \right]^TX \\[2ex] &= \frac{1}{2} X^TPX+Q^TX \end{aligned}$

其中， $P=\left[\begin{matrix} 2 &0 \\ 0&2 \end{matrix} \right]，\ Q=\left[\begin{matrix} -2 \\ -4 \end{matrix} \right]$

关于为什么要写成 $\frac{1}{2} X^TPX$ 形式，因为此时 $P$ 为目标函数的海塞矩阵，具体参看此链接。

2. 约束条件处理

$\begin{aligned} \left\{ \begin{array}{} 0 \leqslant x_1 \leqslant 1.5 \\ 1 \leqslant x_2 \leqslant 2.5 \\ \end{array} \right . \quad \Longleftrightarrow \quad \left[\begin{matrix}0 \\1\end{matrix} \right] \leqslant \left[\begin{matrix}x_1 \\x_2\end{matrix} \right] \leqslant \left[\begin{matrix}1.5 \\2.5\end{matrix} \right] \quad \Longleftrightarrow \quad \left[\begin{matrix}0 \\1\end{matrix} \right] \leqslant \left[\begin{matrix}1&0 \\0&1\end{matrix} \right]X \leqslant \left[\begin{matrix}1.5 \\2.5\end{matrix} \right] \end{aligned}$
令 $L_B=\left[\begin{matrix}0 \\1\end{matrix} \right],\ A=\left[\begin{matrix}1&0 \\0&1\end{matrix} \right],\ U_B=\left[\begin{matrix}1.5 \\2.5\end{matrix} \right]$ ，整理得约束条件如下：
$L_B \leqslant AX \leqslant 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;