Prerequisite

自2022年8月12日至2023年3月10日+的时间内，参加了院内的组合优化方向研究组，我致力于完成优化the uncapacitated facility location problem（UFLP）问题，以下是我的工作内容发展过程以及该工作结束后的总结反思。

Introduction of UFLP

m个顾客，n个仓库，仓库有开设成本fi，每个顾客与仓库之间有运费cij
在n个仓库中选择某几个开设需要使得仓库开设成本+顾客与开设的仓库之间的运费总和最小。

Algorithm introduction

介于组内有相似问题即MKP，故老师授予我用于MKP的新方法，让我挪过来使用。

original algorithm

1）生成随机解作为群体使用
2）根据群体得到分数将所有仓库排序，取中间一部分放入cplex求解，左边部分固定为0，右边部分投票数足够高则被固定为1，得到两部分解拼凑后为新解
3）在cplex中有两个计算部分，分别是仓库和顾客，但若将所有顾客分给cplex则cplex压力过高时间计算会很长，因此分一部分顾客给固定的仓库，其余给cplex。
4）新解与固定的群体随机解对比，若优于群体解则留下，否则淘汰

由于效果不佳，方法不断在进行调整，下面介绍方法修改历程：

Phase 1

Phase 1 process

1）随机解过差，因此对随机解进行调整，优化随机解后再利用其作为群体。

2)在做仓库排序时，涉及到仓库性价比，仓库性价比调整为由两个因素影响：开设费用排名和离顾客最近的仓库排名。利用排名而不是数据本身是为了归一化。

3）在计算大规模instance时，若给cplex数目太多则会造成cplex结果很差所得新解无效，因此需要尽可能降低cplex压力，缩短cplex范围，减少cplex顾客分配。

4）cplex范围采取动态变化，最开始需要搜索分数偏低的仓库，随迭代次数cplex范围右移。

5）仓库固定策略中的参数调整
非cplex部分的仓库投票数大于一个阈值才能被固定为开设。阈值=非cplex部分仓库的平均投票数*zeta 。此处调整的参数即为zeta。

6）顾客分配策略调整

7）上图为群体更新策略，更新分数与解的目标函数值F、海明距离D以及解中所选仓库的性价比相关

8）对于难以辨别的仓库，在达到一定迭代次数后，用cplex搜索投票数目最高的一部分仓库，以达到精确搜索。

Phase 1 experimental results

以下为在M* benchmark上的测试结果，与已知代码的对比算法BinSCA进行了对比

Phase 2

当前解决UFLP的latest算法有两个BinCSA与BinSSA，BinCSA已知代码，也向BinSSA要到了代码不过是Objective-C的，在matlab上运行速度极慢，为了公平起见，我尽力去复现了BinSSA的算法，在论文给出的instance上，都达到了最优解，但均值没有达到。代码放到结尾处了。
（如果是将现有Objective-c直接对照改为c++/c本不是是一件难事，但是作者给的代码与论文上内容不符…，所以还是自己看着论文重新写的）

Phase 2 experimental results

后续在我的算法上具体调整了一些细节和一些参数，跑了M*上的所有例子，并且对于1000x1000和2000x2000的例子，根据生成方法分别随机生成了9个和4个。

M*上的所有例子都跑了30个种子

Phase 3

Phase1和2的操作是为了提升效果而做了很多细节上的修改，给了很多具体的值。

例如1、在2000个仓库中取前30个性价比最好的仓库随机固定5个，重复20000次得到20000个随机解，这是我看着答案做题，显然是不可取的。

这个问题是怎么发现的呢？
现有两个benchmark：M和ORLIB，M上跑的效果很不错，但是在ORLIB上很糟糕，ORLIB还只是小例子。这是因为ORLIB中每个例子最优解中选择的仓库性价比是不一定高的，而M中的都很高，所以M很容易能通过搜索投票数目最高的几个仓库来搜索到最优解。而这里的“5”是因为M*中的最优解总是只选了5个仓库。

2、在1、的基础上目标仓库的投票值一定会在前15左右，那么前几次cplex随便搜索一下，等着第一次的投票最高的仓库集体搜索就能轻易得到最优解，此处的投票最高的仓库集体搜索数目，也是我按照不同规模例子上目标仓库的聚集程度来定的，规模小一点的，那么目标仓库就更容易在投票数最高的前几，规模大一点的，目标仓库就会难聚集一点，因此随着规模增大，cplex仓库集体搜索数目也随之增大，从10增加到15左右。

上述两个操作使得方法本身无效化，ORLIB的存在使得我必须得修改方法。

Phase 3 process

1）随机解策略的调整：
生成10000个随机解策略：从性价比最好的n个仓库中随机挑选20个仓库开设。且设定阈值：n若小于100则n=总仓库数目，反之n=100

2）顾客分配策略调整：
若顾客与该轮固定的仓库的距离大于上一次最优解中分配得到的距离，则该顾客被选入cplex中计算。

3）非cplex部分仓库固定策略
之前的方法是cplex范围右边全部固定为开设，左边全部固定为关闭，现在左边仍然全部固定为关闭，但右边为随机开设一部分仓库。

4）cplex范围调整
cplex范围为k*当前最优解数目且设定大于10的阈值。且cplex范围每迭代10次随当前最优解变化一次。

5）更新群体的参数调整
更新群体的分数由两部分影响：该解的海明距离和最优值。参数是用来调整两部分的权重。之前的策略会让该权重随迭代次数增加变化，现直接固定为定值。

Phase 3 conclusion

该阶段一直在调整各种小细节，挺折磨人的，看到糟糕的效果很难受，但这个阶段有点麻木了，遇到问题总是喜欢各种调一些小参数，但可能在当前的种子下效果好，换一个种子效果又变差了。
Phase 3仍在持续更新中…后续会将我的代码放上来。

Code

BinSSA

复现 A binary social spider algorithm for uncapacitated facility location problem 的代码如下

//
//  main.cpp
//  BinSSA#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include
#includeusing namespace std;string File_Name;
double start_time;
int N = 40;          // Number of spiders (Population Size)
int MAX_ITER = 10000;
double ra = 1, p_c = 0.7, p_m = 0.1, ST = 0.3;
int iter;
int record = 0;
double target_max = -99999;
int inactive_deg;
int *target_vibr;
int whNum, cusNum;
double *f, **dis;
double phi, phi_max = 0.9, phi_min = 0.5;
int **x_spiders;
int **candidate_spider;
double *candidate_memory,  **x_spiders_new, *x_spiders_previous;
double *vibration;
int *dimension_mask;
long long int **distances;
double **spider_rank;
double *history_obj;
int *history_time;void DataInput() {int *temp_1;int *temp_2;ifstream FIC;ofstream FIC2;FIC.open(File_Name);if (FIC.fail()) {cout << "### Error open, File_Name " << File_Name << endl;exit(0);}FIC >> whNum >> cusNum;f = new double [whNum];temp_1 = new int [whNum];temp_2 = new int [cusNum];//读取fixed costfor(int i = 0; i < whNum; i++) {FIC >> temp_1[i] >> f[i];}dis = new double *[cusNum];//[customer][facility]for(int i = 0; i < cusNum; i++) {dis[i] = new  double [whNum];}//读取customer与各个仓库的距离for(int i = 0; i < cusNum; i++) {FIC >> temp_2[i];
//        FIC >> t;for(int j = 0; j < whNum; j++) {FIC >> dis[i][j];}}FIC.close();
}void AssignMemory () {int i;history_time = new int [1000];history_obj = new double[1000];distances = new long long int *[N];target_vibr = new int [whNum];x_spiders = new int *[N];x_spiders_new = new double *[N];x_spiders_previous = new double[whNum];candidate_spider = new int *[N];candidate_memory = new double [N];dimension_mask = new int[whNum];spider_rank = new double *[2*N];vibration = new double [N];for(i=0;idistances[i] = new long long int[N];x_spiders[i] = new int[whNum];x_spiders_new[i] = new double[whNum];candidate_spider[i] = new int [whNum];}for(i=0;i<2*N;i++) {spider_rank[i] = new double[whNum + 1];}}void DeleteMemory () {delete [] target_vibr;delete []x_spiders;delete []x_spiders_new;delete []x_spiders_previous;delete []candidate_spider;delete []candidate_memory;delete []dimension_mask;delete []spider_rank;delete []vibration;delete  [] distances;
}double UFLP(int *warehouse)
{double cost = 0;double min_cost = 0;int i,j;for(i=0; imin_cost = DBL_MAX;for(j=0; jif (warehouse[j] == 1 && dis[i][j]int j, k;double V;double max = target_max;int target = -1;int flag = 1;for(j=0;jflag = 1;for(int i=0;i max) {target = j;max = V;}}if(target != -1) {target_max = max;inactive_deg = 0;for(j=0;jtarget_vibr[j] = x_spiders[target][j];}} else{inactive_deg++;}
}void SpiderWalk (int l) {int i , j, r;double temp;//dimension_mask updatep_m = 0.1;if(rand() / double(RAND_MAX) > pow(p_c, inactive_deg)) {inactive_deg = 0;p_m *= rand() / double(RAND_MAX);for(j=0;jdimension_mask[j] = rand() / double(RAND_MAX) < p_m;}}double p_following;for(j=0;jx_spiders_previous[j] *= rand() / double(RAND_MAX);p_following = dimension_mask[j] ? x_spiders[rand() % N][j] : target_vibr[j];x_spiders_previous[j] += (p_following - x_spiders_new[l][j]) * (rand() / double(RAND_MAX));x_spiders_new[l][j] += x_spiders_previous[j];temp = 1 / (1 + exp(-x_spiders_new[l][j]));// transfer function S2if(temp > 0.51) x_spiders[l][j] = 1;else x_spiders[l][j] = 0;}}void Similarity(int l) {int i, r;int M_11, M_10,M_01;int M11,M10,M01;int nn_1 = 0, nn_0 = 0;int m_01 = 0,m_10 = 0,m_11 = 0;for(i=0;ir = rand() % N;
//        cout<for(M_01 = 0; M_01 <= nn_1; M_01++) {M_11 = nn_1 - M_01;value = abs( 1 - 2.0*M_11 / (2.0*M_11 + M_10 + M_01) - phi*(1 - 2.0*m_11/(2.0*m_11 + m_10 + m_01)) );
//            cout<value_min = value;M01 = M_01;M10 = M_10;M11 = M_11;}}}for(i=0;iif(M11 > 0 && x_spiders[l][i] == 1) {candidate_spider[l][i] = 1;M11--;}else if(M10 > 0 && x_spiders[l][i] == 0) {candidate_spider[l][i] = 1;M10--;}elsecandidate_spider[l][i] = 0;}}void logicGate (int l) {int j, r_2;//[0,N-1];double r_1;r_2 = rand() % N;for(j=0;jr_1 = (double)(rand()) /RAND_MAX;if(r_1 < ST) {if(dimension_mask[j] == x_spiders[r_2][j]) {//sameif(x_spiders[l][j] == 1)//1candidate_spider[l][j] = 1;else//0candidate_spider[l][j] = 0;}else {//diffif(x_spiders[l][j] == 1)//1candidate_spider[l][j] = 0;else//0candidate_spider[l][j] = 1;}}elsecandidate_spider[l][j] = x_spiders[l][j];}}void DataOutput (string Out_file_name) {int i,j;ofstream outfile;outfile.open(Out_file_name, ios::app);if(!outfile) {cout<<"fail"<outfile << history_obj[i]<<" ";}outfile<outfile << history_time[i]<<" ";}outfile<return X[0] < Y[0];
}int main(int argc, char **argv) {double start_time = clock();;double time_max, time_cost=0;string Out_file_name;// int seed = 8;// File_Name = "../M/S/MS1";int seed = atoi(argv[1]);File_Name = argv[2];time_max = atoi(argv[3]);Out_file_name = argv[4];srand(seed);DataInput();AssignMemory();double best_time;double r;int i, j;double C = -1E-100;double global_best = DBL_MAX;double *obj_spiders, *obj_memory;int **x_memory;obj_spiders = new double[N];obj_memory = new double[N];x_memory = new int *[N];for(i=0;ifor(j=0;jx_spiders[i][j] = rand() % 2;//[0,1]x_spiders_new[i][j] = x_spiders[i][j];x_spiders_previous[j] = 0;x_memory[i][j] = x_spiders[i][j];spider_rank[i][j+1]  = x_spiders[i][j];}}for(i=0;iobj_memory[i] = UFLP(x_spiders[i]);}for(j=0;jdimension_mask[j] = rand() % 2;target_vibr[j] = rand() % 2;}//iteraion beginwhile(time_cost < time_max) {//    cout<<"iter: "<obj_spiders[i] = UFLP(x_spiders[i]);if(obj_spiders[i] < obj_memory[i]){obj_memory[i] = obj_spiders[i];for(j=0; j1) {for(i=0;iif (obj_spiders[i] > candidate_memory[i]){for(j=0;jx_spiders[i][j] = candidate_spider[i][j];}obj_spiders[i] = candidate_memory[i];if(obj_spiders[i] < obj_memory[i]) {obj_memory[i] = obj_spiders[i];for(j=0; jglobal_best = obj_memory[i];//    cout<for(j=i; jdistances[i][j] = 0;distances[j][i] = distances[i][j];}}//mahattan distancefor(i=0;ifor(j=i+1; jdistances[i][j] = 0;for(int k =0;kdistances[i][j] += abs(x_spiders[i][k] - x_spiders[j][k]);}distances[j][i] = distances[i][j];}}for(i=0;ivibration[i] = log(1.0/(obj_spiders[i] - C) + 1 );}//compute sigma_barfor(i=0;ix_spider_bar[i] = 0;for(j=0;jsigma[i] = 0;for(j=0;jsigma[i] += (x_spider_bar[i] - x_spiders[j][i]) * (x_spider_bar[i] - x_spiders[j][i]);}sigma[i] = sqrt(sigma[i] / N);sigma_bar += sigma[i];}sigma_bar /= whNum;//    cout << sigma_bar<<" ";//    exit(1);for(i=0; ir = rand()*1.0/RAND_MAX; //[a,b] [0,1]if(r < 0.5) {Similarity(i);}else {logicGate(i);}candidate_memory[i] = UFLP(candidate_spider[i]);}time_cost = (double) (1.0 * (clock() - start_time) / CLOCKS_PER_SEC);iter++;}//DataOutput(Out_file_name);
//    cout<