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 <mao@gnu.org>
* Tested : GLPK v4.45
*
* Waiver: To the extent possible under law, Andrew Makhorin
* <mao@gnu.org> 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 轉換為合適的格式,呼叫相應的求解器,然後嘗試識別可行解,最後將解轉換回原始問題例項,以便進行後續分析和報告。