使用CMake和VS2019编译OSQP库文件并在C/C++中调用求解

前言

OSQP(Operator Splitting Quadratic Program)是一个求解凸二次规划问题的开源求解器，官网地址为https://osqp.org/。

OSQP的内核优化算法是交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)，具体的算法原理可以阅读运筹OR帷幄的专栏文章。

OSQP使用时无需调用其他第三方库，其内置的线性方程组求解算法是QDLDL，在求解中小规模的问题时效果较好；如果想要求解大规模的问题，官方建议使用Intel的MKL Pardiso，也可以使用其他求解器，具体可以参考官方帮助文档。这里只介绍OSQP的基本安装及调用方法。

---

一、库文件编译

下载OSQP源码后，使用CMake和VS2019可以非常简单地实现库文件编译。也可以直接下载我之前已经编译好的0.6.2版本库文件。

1. 源码下载

首先准备OSQP和QDLDL的源码，下载地址分别为OSQP源码和QDLDL源码。下载成功后分别解压，并将QDLDL的源码放到OSQP的lin_sys\direct\qdldl\qdldl_sources文件夹中。

2. 使用CMake生成VS工程

安装CMake后打开，在"Where is the source code"栏选择刚刚下载的OSQP的源码位置，在"Where to build the binaries"栏选择工程的生成位置，然后点击左下的"Configure"按钮进行配置。

在弹出来的窗口中选择编译环境及编译平台，像我的常用开发环境是64为Windows下的VS2019，选择好generator和platform后点击"Finish"进入下一步。

等待自动配置完成后，下方信息框将打印出"Configuring done"提示，中间将出现所有可供配置的选项，可以根据自己的需求进行选择。比如常用的VS工程配置为Debug和Release，则可在"CMAKE_CONFIGURATION_TYPES"栏选择对应值；如果想使用MKL Pardiso求解器，则应选择"ENABLE_MKL_PARDISO"项；如果希望看到求解器打印的迭代过程及结果等信息，则应选择"PRINTING"项。所有项按需选择完成后，再次点击左下的"Configure"按钮完成配置。

 等待再次配置完成后，中间的所有配置项将变白，点击"Generate"按钮开始生成。

等待信息框出现"Generating done"提示，说明VS工程生成成功。此时可以在刚刚选择的工程生成位置看到"osqp.sln"工程文件。

3. 使用VS编译得到库文件

使用VS2019打开刚刚生成的"osqp.sln"工程文件，选择想要的Debug或者Release模式后，生成解决方案即可。生成成功后，可在工程中的out文件夹中看到对应的库文件。

使用时将源码位置中include文件夹和刚刚生成的库文件包含在工程中即可调用OSQP求解凸二次规划问题了。

二、调用求解

凸二次规划问题的标准数学形式为：

求解时，包含"osqp.h"头文件后，只需要将问题规模、矩阵P、A和向量q、l、u输入即可调用OSQP，但是需要注意的是OSQP需要使用CSC matrix进行稀疏矩阵的输入。下面使用官方帮助文档中的Demo进行说明，其对应的C代码可以在examples文件夹中找到。

1. 凸二次规划问题的数据输入

// Load problem data
c_float P_x[3] = { 4.0, 1.0, 2.0, };
c_int   P_nnz  = 3;
c_int   P_i[3] = { 0, 0, 1, };
c_int   P_p[3] = { 0, 1, 3, };
c_float q[2]   = { 1.0, 1.0, };
c_float A_x[4] = { 1.0, 1.0, 1.0, 1.0, };
c_int   A_nnz  = 4;
c_int   A_i[4] = { 0, 1, 0, 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;

// Populate data
OSQPData *data = (OSQPData *)c_malloc(sizeof(OSQPData));
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;
}

// Clean workspace
if (data) {
    if (data->A) c_free(data->A);
    if (data->P) c_free(data->P);
    c_free(data);
}

待求解的凸二次规划问题数据存储在data变量中，其数据类型为OSQPData结构体，是已经定义好的数据类型，其成员分别表示：n-待求解QP问题的决策变量个数；m-线性等式及不等式约束个数之和；P-目标函数的Hessian矩阵的上三角矩阵（由于Hessian矩阵为对称矩阵，只需输入上三角部分即可）；q-目标函数的一次项系数向量；A-线性约束的系数矩阵；l-线性约束的下界；u-线性约束的上界。如果某个约束为线性等式约束，那么其对应的上下界相等。

OSQP使用按列压缩的CSC稀疏矩阵，即将矩阵中每一列的非零元素统计好后依次存储起来。以矩阵P为例：其行数为n，列数为n；P_nnz表示所有非零元素的个数，因为只需输入上三角部分，所以等于3；P_x的维数等于非零元素的个数P_nnz，其中存储的是所有非零元素的值，分别为第一列的4、第二列的1,2；P_i的维数与P_x相同，表示的是非零元素对应的行标，比如4和1是第0行，2是第1行，那么P_i中存储的就是{0,0,1}；P_p的维数等于P的列数n加1，所以也是3，可以理解为其中存储的第1个元素减去第0个元素（第0个元素总是0）是P中第一列的非零元素的个数1，第2个元素减去第1个元素是P中第二列的非零元素的个数2。

需要注意的是data声明时使用了malloc，所以求解完成之后需要使用free将其空间释放掉。

2. 求解器设置

// Define solver settings as default
OSQPSettings *settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings));
if (settings) osqp_set_default_settings(settings);

// Clean workspace
if (settings)  c_free(settings);

求解器设置项存储在settings变量中，其数据类型为OSQPSettings结构体，同样是已经定义好的数据类型，其数据成员为可以设置的选项，具体说明可以参考官方说明文档的Solver settings。也可以直接调用osqp_set_default_settings()函数将所有的选项设置为默认值。

同样settings声明时使用了malloc，求解完成之后需要使用free将其空间释放掉。

3. 调用求解及结果输出

// Workspace structures
OSQPWorkspace *work;

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

// Solve Problem
exitflag = osqp_solve(work);

// Print solution
if (work->info[0].status_val == 1) {
    std::cout << "Solution =" << std::endl;
    for (int idx = 0; idx < n; idx++) {
        std::cout << work->solution[0].x[idx] << std::endl;
    }
}

// Clean workspace
osqp_cleanup(work);

OSQP的所有输入和输出都在work中，其数据类型为OSQPWorkspace结构体，同样是已经定义好的数据类型，其数据成员包括问题数据、求解器选项、以及求解结果等。使用osqp_setup()函数将工作空间配置好后，调用osqp_solve()函数即可求解凸二次规划问题，如果编译库文件时选择了打印求解过程及结果"PRINTING"项，那么信息会自动打印，也可以通过work的成员solution将最优决策变量输出。求解完成之后需要使用osqp_cleanup()函数将工作空间清理掉。

调用osqp_setup()和osqp_solve()函数可以返回一些错误，OSQP的求解状态存储在work->info的status_val成员中，它们的值代表的含义可以参考官方说明文档的Status values and errors。

前面Demo的最终求解结果如下：