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;
}