Feb 12, 2018 - 在 C++ 中使用 GLPK 求解线性规划

最近,参加了一个提交答案类的编程比赛,有一道题可用线性规划解决。搜索发现,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 编译如下代码

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

在我撰写本文时(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),用一个矩阵表示辅助变量和原有变量的关系。

#include <cstdio>
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)
*     4: obj =   1.351000000e+03 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
z = 1351.000000, x1 = 37.750000, x2 = 162.250000
*/

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

一个常见的需求是,需要求出对应的整数解。使用 glpk,这一问题也很好解决。

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)
*     4: obj =   1.351000000e+03 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...
+     4: mip =     not found yet <=              +inf        (1; 0)
Solution found by heuristic: 1346
+     6: >>>>>   1.348000000e+03 <=   1.348000000e+03   0.0% (1; 1)
+     6: mip =   1.348000000e+03 <=     tree is empty   0.0% (0; 3)
INTEGER OPTIMAL SOLUTION FOUND
z = 1348.000000, x1 = 37.000000, x2 = 163.000000
*/

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

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

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

Sep 20, 2015 - 用牛刀 LaTeX 完成简单的排版任务

假期里,我们被布置了一篇报告作为作业。排版要求如下:

正文为宋体、小四号、1.25倍行距。标题号第一层用一、二……第二层用1、2……第三层用(1)、(2)……。

同时,我们的报告要求写在特定纸张上。换句话说,是要用特定的页眉页脚。如图: Example

在布置作业时,我们也得到了一个Word 模板。刚刚读完 《LaTex 入门》,就想着用 LaTeX 来完成这项工作。

自定义纸张

自定义纸张可以用 wallpaper 宏包 来处理。这个包的功能很简单,刚好符合本例的需求。将原有的模板导出成 pdf,然后使用该宏包加载即可。

\CenterWallPaper{1}{background.pdf}

自定义章节标题格式

新的 CTeX 套件 已经能完善地支持标题格式。本例中,可以使用

\ctexset {
	section = {
		name = {,、},
		number = \chinese{section}
	},
	subsection = {
		name = {,、},
		number = \arabic{subsection},
	},
	subsubsection = {
		name = {(,)、},
		numbering = true,
		number = \arabic{subsubsection}
	}
}

不过,这一功能比较新。如果使用的是较旧的套件,可能就无法正确编译。因此,建议首先升级至 texlive 2015.

Jun 7, 2015 - 从 OI 到世界

又到了一年的高考季,而我,当年泡在机房里的高中生,已经毕业了两年。在参加 OI 比赛时,侧重点都在算法上,因而忽视了很多 C 语言的特性。在刚参与开源项目时,也多走了不少弯路。所以,我整理了一个小小的列表,希望这些问题能对你起到些许帮助。

预处理

预处理的潜在风险

C 语言中,预处理与编译是独立的阶段。很多时候,这会引入额外的风险。

  1. 举出至少两个可能带来风险的,在 OI 中较为常用的宏,说明其危险性
  2. inline, const 等方式重写,在不增加额外代价的情况下解决该问题

常用的宏

虽然宏很危险,但很多时候,宏有其不可替代的作用。

  1. 举出至少一个例子,说明宏的必要性
  2. 解释 __FILE__, __FUNCTION__, __LINE__ 等宏的作用

链接器

与许多带有 import 的语言不同,C 语言中链接器与编译器是分开的。链接器的机制,并不是三言两语就可以说清的。所以,在这里,我只是提几个常见的基础问题。

  1. C 语言中,什么是函数的声明?什么是函数的定义?缺失了其中的一个会发生什么?
  2. #include 时发生了什么?为什么有时没有包含 stdio.h 也可以使用 printf 函数?
  3. 什么是动态链接?有什么好处?如何使用?
  4. 头文件中的 extern "C" 有什么用?如果没有,会发生什么?

实战演练

此部分没有固定答案,希望大家的思路不要被拘束。

“带模版”的 qsort

stdlib.h 中的 qsort 有四个参数。它们的意义都是什么?为什么要如此定义?
尝试自己实现一个 my_qsort 函数,要求:

  1. 具有普适性
  2. 在平均情况下为 O(n lgn) 的复杂度。

Better String

C 语言中并没有内置的 string,而只有 char *。请举若干例子,说明这导致了哪些问题。
尝试在不使用编译器扩展特性的情况下,实现一个字符串类型。要求如下:

  1. 求字符串长度的时间复杂度为 O(1)
  2. string.hstrlen, strstr, strcat, strcpy 兼容,或是重写对应的函数
  3. 额外的时间和空间代价必须是常数级的




参考信息

我并不会对这些问题给出标准答案。如果你感觉没什么头绪,可以参考如下信息

举出至少一个例子,说明宏的必要性

  1. Windows API 中的 VARIANT (Wine 中的实现)
  2. 用 C 语言实现链表 等数据结构

C 语言中,什么是函数的声明?什么是函数的定义?缺失了其中的一个会发生什么?

参见 C 语言程序设计 现代方法 第九章 函数

什么是动态链接?有什么好处?如何使用?

参见 程序员的自我修养 链接、装载与库 第7章(Linux) 第9章(Windows)

头文件中的 extern "C" 有什么用?如果没有,会发生什么?

参见 程序员的自我修养 链接、装载与库 3.5.4
Wikipedia Name Mangling

“带模版”的 qsort

参见 C 语言程序设计 现代方法 17.7.2 指针的高级应用

Better String

主要有两种思路:

  1. 使用结构体封装,如 bstrlib
  2. 在指针前附加信息,如 SDS, BSTR

SDS 项目给出了两种实现方式的对比

书籍推荐

一站式学习C编程 程序员的自我修养:链接、装载与库 C语言程序设计:现代方法(第2版) C陷阱与缺陷