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得到的最优解附近利用禁忌算法搜索
- 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;
ILOSTLBEGIN
string 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:11
for(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];
else
IW_new[h] -= R[h][i];
if(x[j] == 1)
IW_new[h] += R[h][j];
else
IW_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
1624
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
