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、1
    c_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矩阵非零元素个数为4
    c_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;

    // Exitflag
    c_int exitflag = 0;
    // Workspace structures
    OSQPWorkspace *work;
    OSQPSettings  *settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings));
    OSQPData      *data     = (OSQPData *)c_malloc(sizeof(OSQPData));

    // Populate data
    if (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 default
    if (settings) {
        osqp_set_default_settings(settings);
        settings->alpha = 1.0; // Change alpha parameter
    }

    // Setup workspace
    exitflag = osqp_setup(&work, data, settings);

    // Solve Problem
    osqp_solve(work);

    // Cleanup
    if (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;
}