OSQP二次规划求解库使用说明

OSQP二次规划求解库使用说明

贺志国
2023.5.10

1. 凸二次规划的一般表达式

$$min\frac{1}{2}{x}^{T}Px+q^{T}xs.t.l\le Ax\le u$$
其中，$P$称为内核矩阵，要求是半正定矩阵。正定矩阵在复数域上是共轭对称矩阵（Hermite矩阵），在实数域上是对称矩阵。因为我们在实数域内求解，$P$也是一个对称矩阵。$q$称为偏移向量，$A$称为仿射矩阵。

凸二次规划的求解算法主要有如下几种：

• 有效集法：Active Set Method
• 内点法：Interior Point Method
• 一阶方法：First Order Method
• 交替方向乘子法：Alternating Direction Method of Multipliers（ADMM）

OSQP利用交替方向乘子法（ADMM）求解，qpOASES则利用有效集法，从网上给出的测评报告看，OSQP库的求解效率比qpOASES库要高很多。

2. 使用OSQP求解凸二次规划的简单示例

为了形象地说明使用OSQP求解凸二次规划的过程，下面给出一个简单的示例：
$$min\frac{1}{2}{x}^T\left[\begin{array}{cc}4& 1\\ 1& 2\end{array}\right]x+{\left[\begin{array}{c}1\\ 1\end{array}\right]}^{T}x$$
s.t. （subject to）：
$$\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]\le \left[\begin{array}{cc}1& 1\\ 1& 0\\ 0& 1\end{array}\right]x\le \left[\begin{array}{c}1\\ 0.7\\ 0.7\end{array}\right]$$
本示例中，内核矩阵$P=\left[\begin{array}{cc}4& 1\\ 1& 2\end{array}\right]$的元素全部为实数，因此它一定是对称矩阵。对于对称矩阵，只需要存储上三角矩阵即可，这也是OSQP存储稀疏矩阵的常见做法，当然也是该库提高计算效率的有效措施。

具体而言，OSQP库使用csc_matrix(indices, indptr, data, csc = Compressed Sparse Column)的方式来存储一个矩阵。csc_matrix存储格式定义如下：

• indices is array of row indices
• data is array of corresponding nonzero values
• indptr points to column starts in indices and data

下面仍然以内核矩阵$P=\left[\begin{array}{cc}4& 1\\ 1& 2\end{array}\right]$为例进行说明。首先，将内核矩阵简化为上三角矩阵：
$$P=\left[\begin{array}{cc}4& 1\\ 0& 2\end{array}\right]$$
$P$的非零元素个数为3，记为：P_nnz = 3 (nnz: number of non-zero elements)

$P$中非零元素为：4, 1, 2，记为：P_x = {4, 1, 2}；

$P$中非零元素行索引为：0, 0, 1，这是按照先第一列，再第二列。。。，最后第N列的顺序标记。第一列的非零元素是4，行索引为0；第二列的非零元素为1, 2，行索引分别为：0, 1，于是记为：P_i = {0, 0, 1}；

$P$中第一列非零元素是1个，第二列非零元素是2个，记为：P_p = {0, 1, 3}。这个记法有些反人类，也很难形象理解，容我再详细解释一番。P_p元素个数为矩阵列数加1，P_p的第一个元素为0，P_p中前后两个元素的差值表示矩阵P相应列的非零元素个数。例如：P_p = {0, 1, 3}表示：1 - 0 = 1 代表矩阵P第一列元素非零个数为１个，3 - 1 = 2 代表矩阵P第二列元素非零个数为２个。助记方法：P_i (i代表indices), P_p(p代表indptr)。

根据上述介绍，下面给出该示例的求解代码：

#include "osqp.h"int main(int argc, char **argv) {// P矩阵是对称矩阵，只存储上三角部分,下三角部分视为零值c_float P_x[3] = {4.0, 1.0, 2.0, }; // P矩阵非零元素个数为3// nnz = number of non-zero elements.c_int P_nnz = 3; // P矩阵非零元素行索引，按照先第一列，再第二列。。。// 的顺序排列。下例表示非零元素的行序号为0、0、1c_int P_i[3] = {0, 0, 1, }; // 1-0=1表示第0列非零元素个数// 3-1=2表示第1列非零元素个数c_int P_p[3] = {0, 1, 3, }; c_float q[2] = {1.0, 1.0, };// A矩阵非零元素c_float A_x[4] = {1.0, 1.0, 1.0, 1.0, };// A矩阵非零元素个数为4c_int A_nnz = 4;// A矩阵非零元素的行序号为0、1、0、2(按列排序）c_int A_i[4] = {0, 1, 0, 2, };// 2-0=2表示第0列2个元素非零// 4-2=2表示第1列2个元素非零c_int A_p[3] = {0, 2, 4, };c_float l[3] = {1.0, 0.0, 0.0, };c_float u[3] = {1.0, 0.7, 0.7, };c_int n = 2;c_int m = 3;// Exitflagc_int exitflag = 0;// Workspace structuresOSQPWorkspace *work;OSQPSettings  *settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings));OSQPData      *data     = (OSQPData *)c_malloc(sizeof(OSQPData));// Populate dataif (data) {data->n = n;data->m = m;data->P = csc_matrix(data->n, data->n, P_nnz, P_x, P_i, P_p);data->q = q;data->A = csc_matrix(data->m, data->n, A_nnz, A_x, A_i, A_p);data->l = l;data->u = u;}// Define solver settings as defaultif (settings) {osqp_set_default_settings(settings);settings->alpha = 1.0; // Change alpha parameter}// Setup workspaceexitflag = osqp_setup(&work, data, settings);// Solve Problemosqp_solve(work);// Cleanupif (data) {if (data->A) c_free(data->A);if (data->P) c_free(data->P);c_free(data);}if (settings) c_free(settings);return exitflag;
}