2022-06-14至2022-08-11 关于复现MKP算法的总结与反思

Prerequisite

自2022年6月14日至2022年8月11日的时间内，我致力于完成A Hybrid Approach for the 0–1 Multidimensional Knapsack problem 论文的复现工作，此次是我第一次进行组合优化方向的学习工作，下面介绍该工作内容发展过程以及该工作结束后的总结反思。

Introduction of — A Hybrid Approach for the 0–1 Multidimensional Knapsack problem

01 knapsack

众所周知的01背包问题，formulation如下所示

MKP 01

将knapsack 01问题可归约为MKP 01
即现在不仅仅是一个单独的属性容量限制，有多个属性容量限制，在满足所有容量限制的情况下要使得背包中的物品价值最大。

method introduction

A Hybrid Approach for the 0–1 Multidimensional Knapsack problem是结合 linear programming & Tabu Search ，实现一个优化的启发式算法来解决MKP01

该方法分为两个步骤，具体如下：

Phase 1: simplex phase

计算解决松弛的MKP来得到一个带小数的最优解
此阶段需要用到整数规划，此处用的是IBM的商业求解器CPLEX

Phase 2: TS phase

在Phase1得到的最优解附近利用禁忌算法搜索

1. Search Space Reduction to gain the Neighborhood
   

2)Tabu management: Tabu list 通过RCS和running list来实现更新

running list：记录当前所有的迭代
RCS：每一次更新running list后，RCS就会通过回溯running list得到当前需要禁忌的对象，以更新tabu list。

3)Evaluation Function
评价函数：MKP中每一个解的评价函数的值是该解超出约束的值

Complish the Code

之前只是看着题目刷题，第一次做照着方法写代码的工作，下面是在coding时遇到的技术学习：
1）增强了C和部分C++语言的掌握能力，如多维数组排序、复杂函数调用、vector容器使用等，能够自主实现所需求的功能。
2）掌握对CPLEX的使用
3）目前debug仅仅只会通过cout和exit一起来间接查看错误点。

Problems

下图表示某一个instance的效果，其中p1的z*表示论文要求效果，p2的z_min表示目前所达到效果。

以下几点为当时推测问题点：

随机数

1、delta[k] = 2 * (u + q - k)
2、delta_max ≈ delta[k]
论文中未给出delta_max具体为多少，因此猜测：
delta_max = delta[k] + rand
经过不断尝试随机值后，效果有一定的改善，但该随机值并未固定。这是每换一个例子且得到效果差时，我第一个进行不断调整的数。

k值的计算

利用在introduction中提及的k值计算方法，paper中得到的值为130（instance 1），与我们得到的值有[1,10]的差距，因此继续检查k的计算方法，其中涉及到z（如下图所示），在paper中仍未给出具体的计算方法，因此一直在自行调整。

Conclusion

虽然最终算法性能未完全复现，但是在此过程中仍然获得了启发式学习的知识、代码能力增强、自我学习能力增强等等。
1）精确算法和启发式算法的学习。例如在coding中就分别使用了CPLEX和Tabu Search。
2）文献的阅读。一开始，一些术语例如“relaxed”, "benchmark"等都不明白是什么，阅读进度很慢，我也很难理解内容。现采用方法：“先读abstarct（了解论文的大致方法）再分别去对应章节细读“。
3) Coding。第一次按照论文给出算法写的代码中，每一个函数里几乎都存在bug。其中总是在数组的长度上出现越界的情况上出错，最开始还找不到哪里出错，也不太会debug，甚至一行一行肉眼看代码来找问题点，后来学会用输出和assert来debug（当然现在是会用断点了）。
4) 碰到一些例子效果不好的情况，总是在一些小地方一直调，想掩盖代码的问题。但这样几乎是无济于事。

Code

下面放出当时的复现代码

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string.h>
#include <time.h>
#include <vector>
#include <math.h>
#include <algorithm>
#include<ilcplex/ilocplex.h>
using namespace std;ILOSTLBEGINstring File_Name;
int **R;
int *P;
int *B;
int N;
int M;
int k ;
int randZ;
vector<float> x_bar;
int seed = 15;void Initializing() {int i;int j;ifstream FIC;ofstream FIC2;FIC.open(File_Name);if (FIC.fail()) {cout << "### Error open, File_Name " << File_Name << endl;exit(0);}FIC >> N >> M;P = new int [N];B = new int [N];R = new int *[M];for(i = 0; i < M; i++) {R[i] = new int [N];}while (! FIC.eof()){for(i = 0; i < N; i++) {FIC >> P[i];}for(i = 0; i < M; i++) {for(j = 0; j < N; j++) {FIC >> R[i][j];//[属性][物品]}}for(j = 0;j < M;j++) {FIC >> B[j];}}FIC.close();
}void Initializing2()
{int i,j;ifstream FIC;ofstream FIC2;FIC.open(File_Name);if (FIC.fail()){cout << "### Erreur open, File_Name " << File_Name << endl;exit(0);}FIC >> N >> M;P = new int [N];B = new int [N];R = new int *[M];for(i = 0; i < M; i++) {R[i] = new int [N];}while (! FIC.eof()){for(i = 0; i < N; i++) {FIC >> P[i];}for(j = 0; j < M; j++) {FIC >> B[j];}for(i = 0; i < M; i++) {for(j = 0; j < N; j++) {FIC >> R[i][j];}}}FIC.close();
}bool cmp1(float a[], float b[])
{return a[0] > b[0];
}bool cmp2(float a[], float b[])
{return a[0] < b[0];
}typedef IloArray<IloIntVarArray>    IntVarMatrix;
typedef IloArray<IloNumVarArray>    NumVarMatrix;
typedef IloArray<IloIntArray>        IntMatrix;
typedef IloArray<IloNumArray>        NumMatrix;int computeValue (vector<int> &x) {int value = 0;for(int i = 0; i < N; i++)value += x[i] * P[i];return value;
}int computeV_b(vector<int> &x) {vector<int> A(M, 0);int v_b = 0;for(int i = 0; i < M; i++) {for(int j = 0; j < N; j++) {A[i] += R[i][j] * x[j];}if(A[i] > B[i]) {v_b += A[i] - B[i];}}return v_b;
}void computeZ()
{float **M2 = new float*[N];int i;int num = rand() % (N / 10 + 1) + N / 10;//[N/10,N/5] instances:1-4,6-8
//    int num = rand() % (N / 20 + 1) + N / 10;//[N/10,3*N/20] instance:5
//    int num = rand() % (N / 144 + 1) + N / 8;// [8/N,19N/144] instances:9-10
//    int num =rand() % (N / 224 + 1) + 15 * N / 112;//[15N/112,N/7] instance:11for(i = 0; i < N; i++) {M2[i] = new float[2];M2[i][0] = P[i];M2[i][1] = i;}sort(M2, M2 + N, cmp2);vector<int> x_1(N, 0);vector<int> x_2(N, 0);int v_b;int item = 0;do {x_1[M2[item][1]] = 1;v_b = computeV_b(x_1);if(v_b > 0) {x_1 = x_2;break;}else {x_2 = x_1;}item++;}while (item < num);randZ = computeValue(x_2);
}int computeK_max() {IloEnv env;IloModel model( env );IloInt i;IloInt j;IloNumVarArray X(env, N, 0, 1,ILOFLOAT);IloExpr k_max(env);IloNumArray x( env, N );for( i = 0; i < M; i++ ) {IloExpr v1( env );for( j = 0; j < N; j++ )v1 += (R[i][j] * X[j]);model.add( v1 <= B[i]);}IloExpr v2( env );for( i = 0; i < N; i++) {v2 += X[i] * P[i];k_max += X[i];}model.add( v2 >= randZ + 1);model.add(IloMaximize(env, k_max));IloCplex cplex(model);cplex.setOut(env.getNullStream());cplex.setParam(IloCplex::Param::RandomSeed,seed);cplex.setParam(IloCplex::Param::Threads, 1);cplex.solve();cplex.getValues(x, X);double obj = cplex.getObjValue();env.end();return ceil(obj);
}int computeK_min() {IloEnv env;IloModel model(env);IloInt i;IloInt j;IloNumVarArray X(env, N, 0, 1,ILOFLOAT);IloExpr k_min(env);IloNumArray x(env, N);for(i = 0; i < M; i++) {IloExpr v1(env);for( j = 0; j < N; j++ ) {v1 += (R[i][j] * X[j]);}model.add( v1 <= B[i]);}IloExpr v2( env );for(i = 0; i < N; i++) {v2 += X[i] * P[i];k_min += X[i];}model.add(v2 >= randZ + 1);model.add(IloMinimize(env, k_min));IloCplex cplex(model);cplex.setOut(env.getNullStream());cplex.setParam(IloCplex::Param::RandomSeed,seed);cplex.setParam(IloCplex::Param::Threads, 1);cplex.solve();cplex.getValues(x, X);double obj = cplex.getObjValue();env.end();return floor(obj);}float *computeX_bar() {IloEnv env;IloModel model(env);IloInt i;IloInt j;IloNumArray x(env, N);IloExpr value(env);IloNumVarArray X(env, N, 0.0, 1.0,ILOFLOAT);IloExpr sum(env);for (i = 0; i < M; i++) {IloExpr v(env);for( j = 0; j < N; j++ ) {v += (R[i][j] * X[j]);}model.add( v <= B[i] );}for(i = 0; i < N; i++)value += (P[i] * X[i]);for(i = 0; i < N; i++)sum += X[i];model.add(sum == k);model.add(IloMaximize(env, value));IloCplex cplex(model);cplex.setOut(env.getNullStream());cplex.setParam(IloCplex::Param::RandomSeed,seed);cplex.setParam(IloCplex::Param::Threads, 1);cplex.solve();cplex.getValues(x, X);env.end();float *x_bar  = new float[N];for(i = 0; i < N; i++) {x_bar[i] = x[i];}return x_bar;
}void swap(float* &a) {for(int i = 0; i < N; i++)x_bar.push_back(a[i]);
}int computeSum(vector<int> &x) {int sum = 0;for(int i = 0; i < N; i++) {sum += x[i];}return sum;
}vector<int> computeXinit () {float **M = new float*[N];for(int i = 0; i < N; i++) {M[i] = new float[2];M[i][0] = x_bar[i];M[i][1] = i ;}sort(M, M + N, cmp1);int index[N];for(int i = 0; i < N; i++) {index[i] = M[i][1];}vector<int> x_init(N,0);for(int i = 0; i < k; i++) {x_init[index[i]] = 1;}return x_init;
}float computeBias(vector<int> &x) {float bias = 0;for(int l = 0; l < N; l++) {bias += abs(x[l] *1.0 - x_bar[l]);}return bias;
}int main(int argc, char **argv)
{string outfilename;string data_out;File_Name = "/Users/jiangtongjun/Library/Containers/com.tencent.xinWeChat/Data/Library/Application\ Support/com.tencent.xinWeChat/2.0b4.0.9/396ec1d845a72fdc3aaa185bb6f9c7eb/Message/MessageTemp/aef63e547d5e98cce87554f37a88a811/File/benchmark2/MK_GK7.DAT";srand(seed);
//    Initializing();Initializing2();int erl = 0;int i;int j;int v_b;int z_min = 0;int z_max;int z;int iter = 0;int v_min;float bias;int bias_max;float bias_old;int i2 = 0;int j2 = 0;int flag;int h = 0;int *RL = new int[500000];int **tabu = new int*[N];int *x_best = new int[N];vector<int> RCS;vector<int> x_temp2;int *IW = new int[M];int *IW_new= new int[M];float *x = new float[N];float *xBar = new float[N];int z_old;int e = 1000;time_t now_time=time(NULL);tm*  t_tm = localtime(&now_time);time_t mk_time_1 = mktime(t_tm);int l;for(i = 0; i < N; ++i) {tabu[i] = new int[N];}for(i = 0; i < N; ++i) {for(j = 0; j < N; ++j) {tabu[i][j] = -1;}}computeZ();k = computeK_max() - computeK_min() + 1;cout<<"k:"<<k<<endl;float *x_temp = computeX_bar();swap(x_temp);int u = 0;int q = 0;for(i = 0; i < N; ++i) {if((x_bar[i] > 0.0) && (x_bar[i] < 1.0)) {q++;}else if( x_bar[i] == 1.0 ) {u++;}}bias_max = 2 * (q + u - k) - 10;cout<<"bias_max:"<<bias_max<<endl;x_temp2 =computeXinit();for(i = 0; i < N; ++i) {x[i] = x_temp2[i];}for(i = 0; i < N; ++i) {xBar[i] = x_bar[i];}for(i = 0; i < M; ++i) {IW[i] = 0;for(j = 0; j < N; j++) {IW[i] += R[i][j] * x[j];}}z = computeValue(x_temp2);v_b = computeV_b(x_temp2);bias = computeBias(x_temp2);if (v_b == 0) {z_min = z;for (i = 0; i < N; ++i)x_best[i] = x[i];}else {z_min = 0;for( l = 0; l < N; ++l)x_best[l] = 0;}while ((v_min < INT_MAX) && (erl < 500000)) {v_min = INT_MAX;z_max = INT_MIN;for(i = 0; i < N; ++i) {for(j = i + 1; j < N; ++j) {if (x[i]+x[j] != 1) {continue;}if(tabu[i][j] != iter) {bias_old = bias;x[i] = 1 - x[i];x[j] = 1 - x[j];if(bias > bias_max - 2.0)bias = bias - abs(1 - x[i] * 1.0 - xBar[i]) - abs(1 - x[j] * 1.0 - xBar[j]) + abs(x[i] * 1.0 - xBar[i]) + abs(x[j] * 1.0 - xBar[j]);if(bias > bias_max) {x[i] = 1 - x[i];x[j] = 1 - x[j];bias = bias_old;continue;}z_old = z;if(x[i] == 1) {z = z + P[i];}else {z = z - P[i];}if(x[j] == 1) {z = z + P[j];}else {z = z - P[j];}if (z <= z_min) {x[i] = 1 - x[i];x[j] = 1 - x[j];z = z_old;bias = bias_old;continue;}v_b = 0;for(h = 0; h < M; ++h) {IW_new[h] = IW[h];if(x[i] == 1)IW_new[h] +=  R[h][i];elseIW_new[h] -=  R[h][i];if(x[j] == 1)IW_new[h] +=  R[h][j];elseIW_new[h] -=  R[h][j];if(IW_new[h] > B[h])v_b += IW_new[h] - B[h];}if((v_b < v_min) || ((v_b == v_min)&&(z > z_max))) {i2 = i;j2 = j;v_min = v_b;z_max = z;}x[i] = 1 - x[i];x[j] = 1 - x[j];z = z_old;bias = bias_old;}}}if (v_min != INT_MAX) {x[i2] = 1 - x[i2];x[j2] = 1 - x[j2];bias = bias - abs(1 - x[i2] * 1.0 - xBar[i2]) - abs(1 - x[j2] * 1.0 - xBar[j2]) + abs(x[i2] * 1.0 - xBar[i2]) + abs(x[j2] * 1.0 - xBar[j2]);z = z_max;for(h = 0; h < M; ++h) {if(x[i2] == 1) {IW[h] +=  R[h][i2];}else {IW[h] -=  R[h][i2];}if(x[j2] == 1) {IW[h] +=  R[h][j2];}else {IW[h] -=  R[h][j2];}}if (v_min == 0) {erl = 0;z_min = z;cout << "z_min :"<< z_min <<endl;for(i = 0; i < N; ++i) {x_best[i] = x[i];}}else {iter++;RL[erl++] = i2;RL[erl++] = j2;i = erl;while(i != 0) {--i;j = RL[i];flag = 0;if(find(RCS.begin(), RCS.end(), j) != RCS.end()) {RCS.erase(std::remove(RCS.begin(), RCS.end(), j), RCS.end());flag = 1;}if(flag == 0) {RCS.push_back(j);}if(RCS.size() == 2) {tabu[RCS[0]][RCS[1]] = iter;tabu[RCS[1]][RCS[0]] = iter;}}RCS.resize(0);}}time_t now_time=time(NULL);tm*  t_tm = localtime(&now_time);time_t mk_time_2 = mktime(t_tm);if(mk_time_2 - mk_time_1 > 36000)break;if(iter > e) {cout << iter << " "<<erl <<endl;e += 1000;}}cout <<"Z_min:" << z_min << endl;
}

最后，感谢Dr.He