GLPK/将 GLPK 与其他求解器软件包混合
用户有时喜欢使用 GLPK 生成问题,但随后使用其他求解器解决它们。或者反过来,将由其他求解器软件包生成的问题导入 GLPK。此页面重点介绍涉及高级编程的选项。
首先,术语“外部求解器软件包”指的是像 CPLEX 或 lp_solve 这样的产品。组合两个软件包的动机多种多样,这里不讨论。
在 GLPK 和另一个求解器软件包之间切换的选项可以在三个级别上解决
- 使用 shell 脚本依次调用每个求解器,通过命令行,同时通过文件系统使用通用问题格式(如 MPS)交换数据
- 使用来自 GLPK 和外部软件包的 API 调用进行高级编程,但仍然通过文件系统交换数据
- 使用来自 GLPK 和外部软件包的 API 调用进行高级编程,但使用内存数据结构而不是文件系统
有关标准问题格式的更多信息,请参阅 互操作性页面。
Osi 计划在本书的 其他地方有介绍。
使用 2.7 或更高版本,COIN-OR CBC 求解器可以在 GLPK mathprog 语言中编译为输入建模语言,尽管目前不支持输出功能,例如数据库输出或控制台 printf 语句。COIN-OR CBC 求解器已为 MIP 求解技术实现了多线程操作,并且目前能够解决比 GLPK 大得多、更复杂得多的 MIP 问题。
COIN-OR CBC 求解器 2.7 版本目前使用 GLPK 4.45,尽管计划在 2.8 版本中使用 4.47。
要在 Linux 或 Max OSX 平台上编译和运行 CBC 求解器,请按照以下步骤操作(稳定版 2.7 使用 Ubuntu 12.10 64 位的全新安装为例)
~$ sudo apt-get install build-essential ~$ sudo apt-get install autoconf ~$ sudo apt-get install subversion ~$ svn co https://projects.coin-or.org/svn/Cbc/stable/2.7 coin-Cbc ~$ cd coin-Cbc/ThirdParty/Glpk ~/coin-Cbc/ThirdParty/Glpk$ ./get.Glpk ~/coin-Cbc/ThirdParty/Glpk$ ./configure ~/coin-Cbc/ThirdParty/Glpk$ cd ../.. ~/coin-Cbc$ ./configure CFLAGS='-m64 -O3' CC=gcc-4.7 --enable-gnu-packages --enable-cbc-parallel --enable-debug -C --prefix=/usr/local ~/coin-Cbc$ make -j 8 [if you want to use parallel compilation - just 'make' otherwise] ~/coin-Cbc$ sudo make install
现在将您的 mathprog 模型复制到机器上并运行模型,将所有计算的变量写入 csv 文件
cbc modelname.mod% -randomCbcSeed 0 -randomSeed 0 -proximity on -threads 8 -printi csv -solve -solu modelname.csv
请注意模型名称末尾的“%”,表明它是一个 mathprog 模型。
最简单的高级方法是使用文件系统在 GLPK 和外部求解器之间传输问题,然后传输解决方案。
例如,GLPK 可以以自由 MPS 格式输出问题,调用外部求解器,然后将解决方案加载回一个新的(第二个)glpk 问题对象。在这种情况下,API 例程glp_write_prob, glp_read_mip和glp_write_mip应该有用。在硬盘 I/O 导致性能下降的极少数情况下,可以部署 RAM 驱动器。
如果外部求解器具有合适的 API,则整个练习可以在程序空间中进行。
以下示例使用 CPLEX。CPLEX 是 IBM 的领先商业求解器套件,于 2009 年被 IBM 收购。在这里,GLPK 用于构建 GLPK 问题,将该问题转换为 CPLEX 问题,调用 CPLEX 解决问题,最后恢复解决方案。以下公共领域 C 代码标题为3party.c正是这样做的。用户应该注意3party.c使用 GLPK 的内部(非公共)数据结构 - 尽管代码可以轻松地重写为使用发布的 API 数据结构。偏离公共接口始终是一个不好的做法。
3party.c | ||
---|---|---|
/* Filename : 3party.c
* Date : July 2011
* Author : Andrew Makhorin <[email protected]>
* Tested : GLPK v4.45
*
* Waiver: To the extent possible under law, Andrew Makhorin
* <[email protected]> has waived all copyright and related or
* neighboring rights to this program code.
* http://creativecommons.org/publicdomain/zero/1.0/
*
* Caution: This code uses internal (non-public) data structures
* and should be rewritten to make use of published APIs instead.
*/
#include "cplex.h"
#include "glpapi.h"
static void cplex_error(CPXENVptr env, const char *who, int status)
{ char buffer[511+1];
if (CPXgeterrorstring(env, status, buffer) == NULL)
xprintf("%s failed\nCPLEX Error %5d: Unknown error code.\n",
who, status);
else
xprintf("%s failed\n%s", who, buffer);
return;
}
int solve_mip(glp_prob *P, const glp_iocp *parm)
{ CPXENVptr env;
CPXLPptr prob = NULL;
CPXFILEptr logfile = NULL;
GLPROW *row;
GLPCOL *col;
GLPAIJ *aij;
int status, numrows, numcols, objsen, numnz, m, n;
int *matbeg, *matcnt, *matind, loc, i, j, k;
char *sense, *ctype;
int *ref = NULL, ret, *indices;
double *matval, *obj, *lb, *ub, *rhs, *rngval, *x, sum;
/* initialize CPLEX environment */
env = CPXopenCPLEX(&status);
if (env == NULL)
{ xprintf("CPXopenCPLEX failed; CPLEX Error %5d.\n", status);
ret = GLP_EFAIL;
goto done;
}
xprintf("CPLEX version is %s\n", CPXversion(env));
/* set CPLEX log file */
logfile = CPXfopen("CONOUT$", "w");
xassert(logfile != NULL);
status = CPXsetlogfile(env, logfile);
if (status)
{ cplex_error(env, "CPXsetlogfile", status);
ret = GLP_EFAIL;
goto done;
}
/* create CPLEX problem object */
prob = CPXcreateprob(env, &status, "MIP");
if (prob == NULL)
{ cplex_error(env, "CPXcreateprob", status);
ret = GLP_EFAIL;
goto done;
}
/* determine original number of rows and columns */
m = P->m, n = P->n;
/* build the row reference array;
ref[i], i = 1,...,m, is the number which i-th row would have
in the problem object after removing all free rows;
if i-th row is free, ref[i] is set to 0;
also count the number of rows and constraint coefficients in
CPLEX problem object */
ref = xcalloc(1+m, sizeof(int));
numrows = 0, numnz = 0;
for (i = 1; i <= m; i++)
{ row = P->row[i];
if (row->type == GLP_FR)
ref[i] = 0;
else
{ ref[i] = ++numrows;
for (aij = row->ptr; aij != NULL; aij = aij->r_next)
numnz++;
}
}
/* the set of columns includes one additional column fixed at 1
to account the constant term of the objective function */
numcols = n+1;
/* allocate working arrays */
obj = xcalloc(numcols, sizeof(double));
sense = xcalloc(numrows, sizeof(char));
rhs = xcalloc(numrows, sizeof(double));
rngval = xcalloc(numrows, sizeof(double));
matbeg = xcalloc(numcols, sizeof(int));
matcnt = xcalloc(numcols, sizeof(int));
matind = xcalloc(numnz, sizeof(int));
matval = xcalloc(numnz, sizeof(double));
lb = xcalloc(numcols, sizeof(double));
ub = xcalloc(numcols, sizeof(double));
/* set objective sense */
if (P->dir == GLP_MIN)
objsen = CPX_MIN;
else if (P->dir == GLP_MAX)
objsen = CPX_MAX;
else
xassert(P != P);
/* set row attributes */
for (k = 1; k <= m; k++)
{ row = P->row[k];
i = ref[k];
if (i == 0) continue;
if (row->type == GLP_LO)
{ sense[i-1] = 'G';
rhs[i-1] = row->lb;
rngval[i-1] = 0.0;
}
else if (row->type == GLP_UP)
{ sense[i-1] = 'L';
rhs[i-1] = row->ub;
rngval[i-1] = 0.0;
}
else if (row->type == GLP_DB)
{ sense[i-1] = 'R';
rhs[i-1] = row->lb;
rngval[i-1] = row->ub - row->lb;
}
else if (row->type == GLP_FX)
{ sense[i-1] = 'E';
rhs[i-1] = row->lb;
rngval[i-1] = 0.0;
}
else
xassert(row != row);
}
/* set column attributes */
for (j = 1; j <= n; j++)
{ col = P->col[j];
obj[j-1] = col->coef;
if (col->type == GLP_FR)
{ lb[j-1] = -CPX_INFBOUND;
ub[j-1] = +CPX_INFBOUND;
}
else if (col->type == GLP_LO)
{ lb[j-1] = col->lb;
ub[j-1] = +CPX_INFBOUND;
}
else if (col->type == GLP_UP)
{ lb[j-1] = -CPX_INFBOUND;
ub[j-1] = col->ub;
}
else if (col->type == GLP_DB)
{ lb[j-1] = col->lb;
ub[j-1] = col->ub;
}
else if (col->type == GLP_FX)
lb[j-1] = ub[j-1] = col->lb;
else
xassert(col != col);
}
obj[numcols-1] = P->c0;
lb[numcols-1] = ub[numcols-1] = 1.0;
/* build the constraint matrix in column-wise format */
loc = 0;
for (j = 1; j <= n; j++)
{ col = P->col[j];
matbeg[j-1] = loc;
for (aij = col->ptr; aij != NULL; aij = aij->c_next)
{ i = ref[aij->row->i];
if (i != 0)
{ matind[loc] = i-1;
matval[loc] = aij->val;
loc++;
}
}
matcnt[j-1] = loc - matbeg[j-1];
}
matbeg[numcols-1] = loc;
matcnt[numcols-1] = 0;
xassert(loc == numnz);
/* copy problem data to CPLEX problem object */
status = CPXcopylp(env, prob, numcols, numrows, objsen, obj, rhs,
sense, matbeg, matcnt, matind, matval, lb, ub, rngval);
/* free working arrays */
xfree(obj);
xfree(sense);
xfree(rhs);
xfree(rngval);
xfree(matbeg);
xfree(matcnt);
xfree(matind);
xfree(matval);
xfree(lb);
xfree(ub);
/* check if copying problem data is successful */
if (status)
{ cplex_error(env, "CPXcopylp", status);
ret = GLP_EFAIL;
goto done;
}
/* change problem type to MIP */
status = CPXchgprobtype(env, prob, CPXPROB_MILP);
if (status)
{ cplex_error(env, "CPXchgprobtype", status);
ret = GLP_EFAIL;
goto done;
}
/* set column types */
indices = xcalloc(numcols+1, sizeof(int));
ctype = xcalloc(numcols+1, sizeof(char));
for (j = 1; j <= n; j++)
{ col = P->col[j];
indices[j-1] = j-1;
if (col->kind == GLP_CV)
ctype[j-1] = 'C';
else if (col->kind == GLP_IV)
ctype[j-1] = 'I';
else
xassert(col != col);
}
indices[numcols-1] = numcols-1;
ctype[numcols-1] = 'C';
status = CPXchgctype(env, prob, numcols, indices, ctype);
xfree(indices);
xfree(ctype);
if (status)
{ cplex_error(env, "CPXchgctype", status);
ret = GLP_EFAIL;
goto done;
}
/* set MIP node log display level */
status = CPXsetintparam(env, CPX_PARAM_MIPDISPLAY, 4);
if (status)
{ cplex_error(env, "CPXsetintparam(CPX_PARAM_MIPDISPLAY)",
status);
ret = GLP_EFAIL;
goto done;
}
/* set solution time limit */
if (parm->tm_lim < INT_MAX)
{ status = CPXsetdblparam(env, CPX_PARAM_TILIM,
(double)parm->tm_lim / 1000.0);
if (status)
{ cplex_error(env, "CPXsetdblparam(CPX_PARAM_TILIM)",
status);
ret = GLP_EFAIL;
goto done;
}
}
/* try to solve the problem */
status = CPXmipopt(env, prob);
if (status)
{ cplex_error(env, "CPXmipopt", status);
ret = GLP_EFAIL;
goto done;
}
/* determine solution status */
status = CPXgetstat(env, prob);
switch (status)
{ case CPXMIP_OPTIMAL:
case CPXMIP_OPTIMAL_TOL:
P->mip_stat = GLP_OPT;
break;
case CPXMIP_TIME_LIM_FEAS:
P->mip_stat = GLP_FEAS;
break;
case CPXMIP_TIME_LIM_INFEAS:
P->mip_stat = GLP_UNDEF;
ret = GLP_ETMLIM;
goto done;
default:
xprintf("CPXgetstat returned %d\n", status);
ret = GLP_EFAIL;
goto done;
}
/* obtain column values */
x = xcalloc(numcols, sizeof(double));
status = CPXgetmipx(env, prob, x, 0, numcols-1);
if (status)
{ cplex_error(env, "CPXgetmipx", status);
xfree(x);
P->mip_stat = GLP_UNDEF;
ret = GLP_EFAIL;
goto done;
}
for (j = 1; j <= n; j++)
{ double t;
col = P->col[j];
if (col->kind == GLP_IV)
{ t = floor(x[j-1] + 0.5);
xassert(fabs(x[j-1] - t) <= 1e-3);
x[j-1] = t;
}
col->mipx = x[j-1];
}
xfree(x);
/* calculate objective value */
sum = P->c0;
for (j = 1; j <= n; j++)
{ col = P->col[j];
sum += col->coef * col->mipx;
}
P->mip_obj = sum;
/* calculate row values */
for (i = 1; i <= m; i++)
{ row = P->row[i];
sum = 0.0;
for (aij = row->ptr; aij != NULL; aij = aij->r_next)
sum += aij->val * aij->col->mipx;
row->mipx = sum;
}
ret = 0;
done: if (ref != NULL)
xfree(ref);
if (prob != NULL)
CPXfreeprob(env, &prob);
if (logfile != NULL)
fclose(logfile);
if (env != NULL)
CPXcloseCPLEX(&env);
return ret;
}
/* eof */
|
CPLEX 是 IBM 的专有优化软件包。这个 2011 年中旬的 主题讨论了 CPLEX 和 GLPK 之间解决方案的转换。
从 4.46 版本开始,GLPK 捆绑了由瑞典查尔姆斯理工大学的 Niklas Eén 和 Niklas Sörensson 开发的 MiniSat CNF-SAT 或 合取范式 可满足性问题 求解器。MiniSat 的宽松 MIT 许可证 允许它作为 GLPK 的一部分进行分发。MiniSat 包装器调用名为glp_minisat1.
4.47 版本引入了glp_intfeas1API,它提供了 MiniSat 代码的原生实现 - 请注意,此函数在此版本中标记为“暂定”。GLPSOL 选项--minisat现在重定向到这个新的求解器。
可满足性 (SAT) 问题涉及确定给定的布尔公式是否可以通过其变量的适当赋值计算为 TRUE。合取范式 (CNF) 指示了布尔公式表示方式的限制。以其他方式表示的布尔公式可以转换为 CNF。
GLPK 官方文档 文件doc/cnfsat.pdf包含此功能的完整参考。因此,这些资料不会在此重复。请参阅 获取 GLPK。
GLPK 以 DIMACS 格式读取和写入 CNF 可满足性问题(另请参阅 GLPK 格式)。对问题规范施加某些限制,使其仍然是有效的 CNF-SAT 问题。这些限制可以通过glp_check_cnfsat公共调用进行显式检查,并且在 GLPSOL--minisat选项部署时会自动确认。
在 GLPK 中,CNF-SAT 被认为是 MIP 的特例,其中所有变量都是二进制的,所有约束都是覆盖不等式。也就是说,特定的 CNF-SAT 问题可以存储在类型为glp_prob的 GLPK 问题对象中,就好像它是 MIP 实例一样。特别是,这意味着 CNF-SAT 问题可以使用glp_intopt求解器解决,但调用glp_minisat1或glp_infeas1应该是问题特定的,效率应该高得多。
在后一种情况下,GLPK 将其 MIP 转换为合适的格式,调用相应的求解器,然后尝试识别可行解,最后将解转换回原始问题实例,以便进行后续分析和报告。