在 C++ 中使用 GLPK 求解线性规划

Zhenbo Li February 12, 2018

最近，参加了一个提交答案类的编程比赛，有一道题可用线性规划解决。搜索发现，GLPK (GNU Linear Programming Kit) 是一个免费的线性规划计算库，可以方便地被 C/C++ 代码调用。现将基本使用方法整理如下：

准备工作

安装 GLPK

多数的发行版都应该提供了 GLPK 的包。在 Fedora 下，只需运行
sudo dnf install glpk glpk-devel

编译 GLPK

Fedora 将 glpk.h 存在了 /usr/include/glpk.h，因此，不需添加指令即可找到头文件。如果你是手动编译安装的 GLPK，那可能需要使用 -I 指定头文件目录。链接的指令为 -lglpk。如果使用的是 C++ 的话，在调用头文件时，记得声明 extern。你可以试着使用 g++ -lglpk test.cpp 编译如下代码

{% highlight C %} #include extern "C"{ #include "glpk.h" } int main() { std::cout << glp_version() << std::endl; return 0; } {% endhighlight %}

在我撰写本文时（2018年2月初），GLPK 的版本应当为 4.61。

实战演练

GLPK 官方手册上的例子有点让人混淆。在这里，我将采用更精简的例子。
Maximize \[z=10x_1+6x_2\] Subject to \[x_1+x_2\leqslant200\] \[x_1+2x_2\geqslant10\] \[3x_1+x_2\leqslant275.5\] where all variables are non-negative \[x_1\geqslant 0, x_2\geqslant 0\]

对于三个约束条件，我们可以创建三个辅助变量(auxiliary variables)，将问题转化为如下形式：
Maximize \[z=10x_1+6x_2\] Subject to \[p=x_1+x_2\] \[q=x_1+2x_2\] \[r=3x_1+x_2\] where all variables are non-negative \[x_1\geqslant 0, x_2\geqslant 0\] \[0\leqslant p\leqslant200,q\geqslant10,0\leqslant r\leqslant275.5\]

现在，可以将我们的问题输入程序了。GLPK 将各辅助变量看作行(row)，将原有的变量看作列(column)，用一个矩阵表示辅助变量和原有变量的关系。

{% highlight c %} #include extern "C"{ #include "glpk.h" }

int main() { initialize: glp_prob *lp; lp = glp_create_prob(); glp_set_obj_dir(lp, GLP_MAX); auxiliary_variables_rows: glp_add_rows(lp, 3); glp_set_row_name(lp, 1, "p"); glp_set_row_bnds(lp, 1, GLP_DB, 0.0, 200.0); glp_set_row_name(lp, 2, "q"); glp_set_row_bnds(lp, 2, GLP_LO, 10.0, 0.0); glp_set_row_name(lp, 3, "r"); glp_set_row_bnds(lp, 3, GLP_DB, 0.0, 275.5);

variables_columns: glp_add_cols(lp, 2); glp_set_col_name(lp, 1, "x1"); glp_set_col_bnds(lp, 1, GLP_LO, 0.0, 0.0); glp_set_col_name(lp, 2, "x2"); glp_set_col_bnds(lp, 2, GLP_LO, 0.0, 0.0); to_maximize: glp_set_obj_coef(lp, 1, 10.0); glp_set_obj_coef(lp, 2, 6.0);

constrant_matrix: int ia[7], ja[7]; double ar[7]; ia[1] = 1, ja[1] = 1, ar[1] = 1; ia[2] = 1, ja[2] = 2, ar[2] = 1; // p = x1 + x2 ia[3] = 2, ja[3] = 1, ar[3] = 1; ia[4] = 2, ja[4] = 2, ar[4] = 2; // q = x1 + 2x2 ia[5] = 3, ja[5] = 1, ar[5] = 3; ia[6] = 3, ja[6] = 2, ar[6] = 1; // r = 3x1 + x2 glp_load_matrix(lp, 6, ia, ja, ar);

calculate: glp_simplex(lp, NULL);

output: double z, x1, x2; z = glp_get_obj_val(lp); x1 = glp_get_col_prim(lp, 1); x2 = glp_get_col_prim(lp, 2); printf("z = %lf, x1 = %lf, x2 = %lf\n", z, x1, x2);

cleanup: glp_delete_prob(lp); return 0; } /* GLPK Simplex Optimizer, v4.61 3 rows, 2 columns, 6 non-zeros 0: obj = -0.000000000e+00 inf = 1.000e+01 (1) 1: obj = 3.000000000e+01 inf = 0.000e+00 (0)

OPTIMAL LP SOLUTION FOUND z = 1351.000000, x1 = 37.750000, x2 = 162.250000 */ {% endhighlight %}

第一次看到示例代码，对于 constrant_matrix 部分可能会比较费解。在这里，ia 表示行号(第几个辅助变量)，ja 表示列号(第几个变量)，而 ar 的类型是 double，表示 constrant matrix 中的系数。将这些条件成功导入后，使用 glp_simplex 就可以求解了。

一个常见的需求是，需要求出对应的整数解。使用 glpk，这一问题也很好解决。
{% highlight c %} set_variables_to_integer: glp_set_col_kind(lp, 1, GLP_IV); glp_set_col_kind(lp, 2, GLP_IV); calculate: glp_simplex(lp, NULL); glp_intopt(lp, NULL); output: double z, x1, x2; z = glp_mip_obj_val(lp); x1 = glp_mip_col_val(lp, 1); x2 = glp_mip_col_val(lp, 2); printf("z = %lf, x1 = %lf, x2 = %lf\n", z, x1, x2); /* GLPK Simplex Optimizer, v4.61 3 rows, 2 columns, 6 non-zeros 0: obj = -0.000000000e+00 inf = 1.000e+01 (1) 1: obj = 3.000000000e+01 inf = 0.000e+00 (0)

OPTIMAL LP SOLUTION FOUND GLPK Integer Optimizer, v4.61 3 rows, 2 columns, 6 non-zeros 2 integer variables, none of which are binary Integer optimization begins...

Solution found by heuristic: 1346

INTEGER OPTIMAL SOLUTION FOUND z = 1348.000000, x1 = 37.000000, x2 = 163.000000 */ {% endhighlight %}

在运行过glp_simplex后，再运行glp_intopt即可得到整数解。但要注意的是，提取整数解的命令是 glp_mip_obj_val / glp_mip_col_val。

PS 写这篇 blog 时，我第一次使用了 mathjax。只能说，写作体验超乎想象地好。

讲一个逸闻吧。之前我曾把线性规划和高斯消元的时间复杂度弄混了，以至于我成功将大量NP问题转化为了线性规划问题后得到了“多项式”的解。

{% include mathjax.html %}