GLPK全称为GNU Linear Programming Kit(GNU线性规划工具包)，可从 https://sourceforge.net/projects/winglpk/ 下载源码及二进制库，最新版本为4.65。也可从 https://ftp.gnu.org/gnu/glpk/ 下载，仅包含源码，最新版本为5.0。
GLPK是用C实现的，旨在解决线性规划(LP, Linear Programming)、混合整数线性规划(MIP, Mixed Integer linear Programming)和其他相关问题。MIP问题是一种LP问题，其中某些变量需要额外满足整数条件。
之前在 https://blog.csdn.net/fengbingchun/article/details/146341096 中介绍过通过Pyomo调用glpk解决饮食问题，这里通过调用GLPK API接口实现，测试代码如下：
int test_glpk()
{
constexpr int total_foods{ 9 };
const std::vector<std::string> foods{ "food1", "food2", "food3", "food4", "food5", "food6", "food7", "food8", "food9" };
const std::vector<float> unit_prices{ 0.886, 0.863, 0.855, 0.917, 0.237, 0.856, 0.833, 0.904, 0.876 };
const std::vector<std::array<float, 2>> bounds{ {1.0, 55.0}, {1.0, 55.0}, {1.0, 55.0}, {1.0, 55.0}, {7.0, 9.0}, {1.0, 55.0}, {1.0, 55.0}, {1.0, 55.0}, {1.0, 55.0} };
constexpr int num_nutrients{ 4 };
const std::vector<std::string> nutrients{ "nutrient1", "nutrient2", "nutrient3", "nutrient4" };
const std::vector<std::array<float, 4>> nutrient_values{ {0.21, 65.10, 0.72, 88.5}, {0.08, 64.58, 0.63, 76.9}, {0.22, 64.81, 0.53, 86.1}, {0.58, 65.84, 1.09, 57.8}, {0.033, 46.07, 14.15, 0},
{0.059, 65.25, 0.39, 87.2}, {0.14, 64.25, 0.57, 93.7}, {0.033, 63.06, 1.36, 79.0}, {0.076, 65.20, 0.59, 99.0} };
const std::vector<std::array<float, 2>> nutrient_limit{ {0., 49.}, {6200., 6230.}, {9.9, 782.}, {6500., 7000.} };
constexpr float total_quantity{ 99. };
const std::string mandatory_food{ "food5" };
constexpr int num_select_foods{ 5 };
if (foods.size() != total_foods) {
std::cerr << "Error: number of foods mismatch" << std::endl;
return -1;
}
if (unit_prices.size() != total_foods) {
std::cerr << "Error: number of unit_prices mismatch" << std::endl;
return -1;
}
if (bounds.size() != total_foods) {
std::cerr << "Error: number of bounds mismatch" << std::endl;
return -1;
}
if (nutrient_values.size() != total_foods) {
std::cerr << "Error: number of nutrient_values mismatch" << std::endl;
return -1;
}
if (nutrients.size() != num_nutrients) {
std::cerr << "Error: number of nutrients mismatch" << std::endl;
return -1;
}
if (nutrient_limit.size() != num_nutrients) {
std::cerr << "Error: number of nutrient_limit mismatch" << std::endl;
return -1;
}
int mandatroy_food_index{ 0 };
for (auto i = 0; i < total_foods; ++i) {
if (mandatory_food == foods[i]) {
mandatroy_food_index = i;
break;
}
}
auto lp = glp_create_prob(); // create problem object
glp_set_prob_name(lp, "diet problem");
glp_set_obj_dir(lp, GLP_MIN); // minimize objective function
// add col variables
glp_add_cols(lp, 2 * total_foods); // first total_foods: number of each food; last total_foods: binary variables, whether the food is selected
for (auto j = 0; j < total_foods; ++j) {
// set food number variable
glp_set_col_name(lp, j + 1, foods[j].c_str());
glp_set_col_kind(lp, j + 1, GLP_CV); // GLP_IV: number of foods can only be an integer, use MIP; default is GLP_CV
glp_set_col_bnds(lp, j + 1, GLP_LO, 0.0, 0.0);
glp_set_obj_coef(lp, j + 1, unit_prices[j]); // set cost coefficient
// set binary select variables
std::string name = "S_" + foods[j];
glp_set_col_name(lp, total_foods + j + 1, name.c_str());
glp_set_col_kind(lp, total_foods + j + 1, GLP_BV);
glp_set_obj_coef(lp, total_foods + j + 1, 0.);
}
// add row constraints
glp_add_rows(lp, num_nutrients);
for (auto i = 0; i < num_nutrients; ++i) {
glp_set_row_name(lp, i + 1, nutrients[i].c_str());
glp_set_row_bnds(lp, i + 1, GLP_DB, nutrient_limit[i][0], nutrient_limit[i][1]);
}
// add select constraints
glp_add_rows(lp, 2);
glp_set_row_name(lp, num_nutrients + 1, "select_total");
glp_set_row_bnds(lp, num_nutrients + 1, GLP_FX, num_select_foods, num_select_foods);
// add mandatory food constraint
glp_set_row_name(lp, num_nutrients + 2, "must_have_food");
glp_set_row_bnds(lp, num_nutrients + 2, GLP_FX, 1., 1.);
// add total_quantity constraint
glp_add_rows(lp, 1);
glp_set_row_name(lp, num_nutrients + 3, "total_quantity");
glp_set_row_bnds(lp, num_nutrients + 3, GLP_FX, total_quantity, total_quantity);
// constraint matrix:
// glp_load_matrix(lp, ne, ia, ja, ar): for k = 1,..., ne, where ia[k] is the row index, ja[k] is the column index, and ar[k] is a numeric value of corresponding constraint coefficient
// parameter ne specifies the total number of (non-zero) elements in the matrix to be loaded
int ia[1 + 1000]{}, ja[1 + 1000]{}, k{ 0 };
double ar[1 + 1000]{};
// nutrients constraint
for (auto i = 0; i < num_nutrients; ++i) {
for (auto j = 0; j < total_foods; ++j) {
k++;
ia[k] = i + 1, ja[k] = j + 1;
ar[k] = nutrient_values[j][i];
}
}
// select total constraint
for (auto j = 0; j < total_foods; ++j) {
k++;
ia[k] = num_nutrients + 1, ja[k] = total_foods + j + 1;
ar[k] = 1.;
}
// mandatory food constraint
k++;
ia[k] = num_nutrients + 2, ja[k] = total_foods + mandatroy_food_index + 1;
ar[k] = 1.;
// total_quantity constraint
for (auto j = 0; j < total_foods; ++j) {
k++;
ia[k] = num_nutrients + 3, ja[k] = j + 1;
ar[k] = 1.;
}
// add relationship constraint between food quantity and select variable
for (auto j = 0; j < total_foods; ++j) {
// lower constraint: quantity_j >= min_j * select_j
glp_add_rows(lp, 1);
glp_set_row_name(lp, num_nutrients + 4 + 2 * j, "quantity_min");
glp_set_row_bnds(lp, num_nutrients + 4 + 2 * j, GLP_LO, 0., 0.);
// quantity_j coefficient: 1.0
k++;
ia[k] = num_nutrients + 4 + 2 * j, ja[k] = j + 1;
ar[k] = 1.;
// select_j cofficient: -min_j
k++;
ia[k] = num_nutrients + 4 + 2 * j, ja[k] = total_foods + j + 1;
ar[k] = -bounds[j][0];
// upper constraint: quantity_j <= max_j * select_j
glp_add_rows(lp, 1);
glp_set_row_name(lp, num_nutrients + 4 + 2 * j + 1, "quantity_max");
glp_set_row_bnds(lp, num_nutrients + 4 + 2 * j + 1, GLP_UP, 0., 0.);
k++;
ia[k] = num_nutrients + 4 + 2 * j + 1, ja[k] = j + 1;
ar[k] = 1.;
k++;
ia[k] = num_nutrients + 4 + 2 * j + 1, ja[k] = total_foods + j + 1;
ar[k] = -bounds[j][1];
}
glp_load_matrix(lp, k, ia, ja, ar);
// verify constraint matrix and variables
std::cout << "number of constraints: " << glp_get_num_rows(lp) << std::endl;
for (auto i = 1; i <= glp_get_num_rows(lp); ++i)
std::cout << "i: " << i << "; type: " << glp_get_row_type(lp, i) << "; lower: " << glp_get_row_lb(lp, i) << "; upper: " << glp_get_row_ub(lp, i) << std::endl;
std::cout << "number of variables: " << glp_get_num_cols(lp) << std::endl;
for (auto j = 1; j <= glp_get_num_cols(lp); ++j) // Note: the value of glp_get_col_prim is different in different positions
std::cout << "j: " << j << "; kind: " << glp_get_col_kind(lp, j) << "; primal value: " << glp_get_col_prim(lp, j) << "; lower: " << glp_get_col_lb(lp, j) << "; upper: " << glp_get_col_ub(lp, j) << "; coef: " << glp_get_obj_coef(lp, j) << std::endl;
auto ret = glp_write_lp(lp, nullptr, "../../../testdata/model.lp");
if (ret != 0) {
std::cerr << "Error: failed to write problem data, error code: " << ret << std::endl;
glp_delete_prob(lp);
return -1;
}
double quantity_sum{ 0. };
if (0) { // LP(Linear Programming)
glp_smcp parm{};
glp_init_smcp(&parm);
parm.msg_lev = GLP_MSG_ERR; // warning and error messages only
parm.presolve = GLP_ON;
ret = glp_simplex(lp, &parm); // solve LP problem with the primal or dual simplex method
if (ret != 0) {
std::cerr << "Error: failed to solve: error code: " << ret << std::endl;
glp_delete_prob(lp);
return -1;
}
ret = glp_get_status(lp);
if (ret != GLP_OPT) {
std::cerr << "Error: there is no optimal solution, status: " << ret << std::endl;
glp_delete_prob(lp);
return -1;
}
std::cout << "minimum cost: " << glp_get_obj_val(lp) << std::endl;
for (auto j = 0; j < total_foods; ++j) {
auto qty = glp_get_col_prim(lp, j + 1);
auto selected = glp_get_col_prim(lp, total_foods + j + 1);
if (selected > 0.5) {
std::cout << foods[j] << ": quantity: " << qty << ", limit: [" << bounds[j][0] << "," << bounds[j][1] << "]" << std::endl;
quantity_sum += qty;
}
}
}
else { // MIP(Mixed Integer linear Programming)
glp_iocp parm;
glp_init_iocp(&parm);
parm.presolve = GLP_ON;
parm.msg_lev = GLP_MSG_ERR; // close information output prompt
ret = glp_intopt(lp, &parm);
if (ret != 0) {
std::cerr << "Error: failed to solve: error code: " << ret << std::endl;
glp_delete_prob(lp);
return -1;
}
ret = glp_mip_status(lp);
if (ret != GLP_OPT) {
std::cerr << "Error: there is no optimal solution, status: " << ret << std::endl;
glp_delete_prob(lp);
return -1;
}
std::cout << "minimum cost: " << glp_mip_obj_val(lp) << std::endl;
for (auto j = 0; j < total_foods; ++j) {
auto qty = glp_mip_col_val(lp, j + 1);
auto selected = glp_mip_col_val(lp, total_foods + j + 1);
if (selected > 0.5) {
std::cout << foods[j] << ": quantity: " << qty << ", limit: [" << bounds[j][0] << "," << bounds[j][1] << "]" << std::endl;
quantity_sum += qty;
}
}
}
std::cout << "result quantity sum: " << quantity_sum << "; require quantity sum: " << total_quantity << std::endl;
glp_delete_prob(lp);
return 0;
}测试数据说明：饮食问题，目标是选择以最低成本满足每日营养需求的食物，测试数据无实际意义
total_foods：食物种类总数。
foods：每种食物的名称。
unit_prices：每种食物的价格。
bounds：每种食物可取的食物份数限制，变量的界限。
num_nutrients：每种食物含有的营养成分个数，每种食物含有相同的营养成分个数。
nutrients：每种食物含有的营养成分名称。
nutrient_values：每种食物含有的营养成分含量。
nutrient_limit：每种营养成分消耗总量约束。
total_quantity：被选中食物总份数约束。
mandatory_food：被选中食物中必选要包含的食物约束。
num_select_foods：被选中食物种类数约束。
使用到的GLPK函数：接口详细描述参考源码中的doc/glpk.pdf
glp_create_prob：创建一个新的问题对象(problem object)，该对象最初是"empty"，即没有行和列。该函数返回指向所创建对象的指针，该指针应在对该对象的任何后续操作中使用。
glp_delete_prob：删除问题对象，释放分配给该对象的所有内存。
glp_set_prob_name：给问题对象分配问题名，长度要求1至255个字符，如果参数名称为nullptr或空字符串，则将删除问题对象的现有问题名。
glp_set_obj_dir：设置(更改)目标函数优化方向标志，GLP_MIN表示最小化目标函数；GLP_MAX表示最大化目标函数。默认情况，最小化目标函数。
glp_add_cols：将新列添加到问题对象。新列始终添加到列表的末尾，因此现有列的序数不会改变。每个新添加的列初始值固定为零，并且约束系数列表为空。
glp_set_col_name：分配(更改)第j列列名，长度要求1至255个字符，如果参数名称为nullptr或空字符串，则将删除第j列的现有名称。
glp_set_col_kind：按照参数类型指定的方式设置(更改)第j列的类型(kind)。GLP_CV表示连续变量；GLP_IV表示整数变量；GLP_BV表示二元变量。
glp_get_col_kind：返回第j列的类型(kind)。
glp_set_col_bnds(lp,j,type,lb,ub)：设置(更改)第j列类型和界限。支持的类型如下图所示：如果列没有下限(lower bound)，则忽略参数lb。如果列没有上限(upper bound)，则忽略参数ub。如果列是固定(fixed)类型，则仅使用参数lb，而忽略参数ub。添加到问题对象后，每列的初始值都固定为零，即其类型为GLP_FX，且两个界限均为0。
glp_get_num_cols：返回指定问题对象中的列数。
glp_set_obj_coef(lp,j,coef)：设置(更改)第j列目标系数(objective coefficient)。新的目标系数值由参数coef指定。
glp_get_obj_coef：返回第j列的目标系数。
glp_get_col_prim：返回与第j列关联变量的原始值(primal value)。注：此函数在调用glp_simplex前后得到值不同
glp_get_col_lb：返回第j列的下限，即相应变量的下限。如果该列没有下限，则返回-DBL_MAX。
glp_get_col_ub：返回第j列的上限，即相应变量的上限。如果该列没有上限，则返回+DBL_MAX。
glp_add_rows：将新行添加到问题对象。新行始终添加到行列表的末尾，因此现有行的序号不会改变。添加的每个新行最初都是free(unbounded)，并且约束系数列表为空。
glp_set_row_name：分配(更改)第i行行名，长度要求1至255个字符，如果参数名称为nullptr或空字符串，则将删除第i行的现有名称。
glp_set_row_bnds(lp,i,type,lb,ub)：设置(更改)第i行类型和界限。支持的类型与glp_set_col_bnds相同。如果行没有下限(lower bound)，则忽略参数lb。如果行没有上限(upper bound)，则忽略参数ub。如果行是固定(fixed)类型，则仅使用参数lb，而忽略参数ub。添加到问题对象中的每一行最初都是free，即其类型为GLP_FR。
glp_get_num_rows：返回指定问题对象中的行数。
glp_get_row_type：返回第i行的类型，支持的类型与glp_set_col_bnds相同。
glp_get_row_lb：返回第i行的下限，如果该行没有下限，则返回-DBL_MAX。
glp_get_row_ub：返回第i行的上限，如果该行没有上限，则返回+DBL_MAX。
glp_load_matrix(lp,ne,ia,ja,ar)：加载(替换)整个约束矩阵(constraint matrix)，将传入数组ia,ja和ar的约束矩阵加载到指定的问题对象中。加载前，约束矩阵的当前内容将被销毁。
约束系数(约束矩阵的元素)应指定为三元组(triplets)(ia[k], ja[k], ar[k])，其中k=1, ..., ne，其中ia[k]是行索引，ja[k]是列索引，ar[k]是相应约束系数(constraint coefficient)的数值。参数ne指定要加载的矩阵中(非零)元素的总数。不允许使用索引相同的系数(同一行和列的索引组合只能出现一次)。允许使用零系数，但是，它们不会存储在约束矩阵中。如果参数ne为0，则参数ia、ja和/或ar可以指定为nullptr。
glp_write_lp：将LP/MIP问题数据以CPLEX LP格式写入文本文件。可用于检查约束是否正确。
glp_init_smcp：使用默认值初始化simplex求解器(solver)使用的控制参数。
glp_simplex：使用simplex方法求解线性规划问题，并将计算结果存储回问题对象。
glp_get_status：报告指定问题对象当前基本解的通用状态，返回可能的状态如下图所示：
glp_get_obj_val：返回目标函数的当前值。求目标函数的最小值或最大值，由glp_set_obj_dir指定。
glp_init_iocp：使用默认值初始化branch-and-cut求解器使用的控制参数。
glp_intopt：使用branch-and-cut方法求解混合整数线性规划问题。
glp_mip_status：报告MIP求解器找到的MIP解决方案的状态，返回可能的状态如下图所示：
glp_mip_obj_val：返回MIP解决方案的目标函数值。求目标函数的最小值或最大值，由glp_set_obj_dir指定。
glp_mip_col_val：返回与MIP解决方案的第j列关联变量的值。
使用到的MIP接口包括：glp_set_col_kind/glp_get_col_kind、glp_init_iocp、glp_intopt、glp_mip_status、glp_mip_obj_val、glp_mip_col_val。
按照以上给的测试数据，共需18个变量(cols)；共需25个约束(rows)。执行结果如下图所示：
注：在GLPK中，C API要求矩阵元素从索引1开始，row[i]对应第i个约束条件，col[j]对应第j个决策变量。通过约束矩阵即ia(行)、ja(列)、ar(系数)关联二者。
GitHub：https://github.com/fengbingchun/Messy_Test
