GLPK(GNU线性规划工具包)介绍

      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 objectglp_set_prob_name(lp, "diet problem");glp_set_obj_dir(lp, GLP_MIN); // minimize objective function// add col variablesglp_add_cols(lp, 2 * total_foods); // first total_foods: number of each food; last total_foods: binary variables, whether the food is selectedfor (auto j = 0; j < total_foods; ++j) {// set food number variableglp_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_CVglp_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 variablesstd::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 constraintsglp_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 constraintsglp_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 constraintglp_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 constraintglp_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 loadedint ia[1 + 1000]{}, ja[1 + 1000]{}, k{ 0 };double ar[1 + 1000]{};// nutrients constraintfor (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 constraintfor (auto j = 0; j < total_foods; ++j) {k++;ia[k] = num_nutrients + 1, ja[k] = total_foods + j + 1;ar[k] = 1.;}// mandatory food constraintk++;ia[k] = num_nutrients + 2, ja[k] = total_foods + mandatroy_food_index + 1;ar[k] = 1.;// total_quantity constraintfor (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 variablefor (auto j = 0; j < total_foods; ++j) {// lower constraint: quantity_j >= min_j * select_jglp_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.0k++;ia[k] = num_nutrients + 4 + 2 * j, ja[k] = j + 1;ar[k] = 1.;// select_j cofficient: -min_jk++;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_jglp_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 variablesstd::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 positionsstd::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 onlyparm.presolve = GLP_ON;ret = glp_simplex(lp, &parm); // solve LP problem with the primal or dual simplex methodif (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 promptret = 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