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
老师提出，是否可以将此算法与精确算法结合
Phase 2-1
在M*中，几乎所有instances所选的目标仓库数目都是在5左右，因此老师认为可以从1个仓库开始往上遍历，并且每遍历完开设n个仓库时，使用以下判断：
如果 n+1个最小开设费用 - 当前best n 开设费用，也就是最小需要增加的仓库开设费用 < 最多能节省的运费 ，那这个时候就无法证明开设n个仓库一定最好，可能开设n+1个更好，因此需要遍历n+1个仓库
小例子上解的较快，但例子稍微大一点（例如500x500的MR里的instances）就需要很长时间，此前跑了24h也没跑出来。Code在最后附上。
Phase 2-2
老师后来提出，可以使用整数规划（CPLEX）代替挨个遍历，并且加上上述的约束条件，且在CPLEX中加上条件：解必须优于当前最优解，这样可能可以再提高一点速度。并且可以将启发式算法结合进该精确算法，即基于上述精确算法，在CPLEX中所加的条件为优于启发式算法得到的精确解。这样得到的效果确实好了很多，但是速度仍然差于当前对于UFLP的最优精确算法，所以该方法被抛弃了。
利用整数规划的精确算法的代码也附在最后了。
Phase 3
当前解决UFLP的latest算法有两个BinCSA与BinSSA，BinCSA已知代码，也向BinSSA要到了代码不过是Objective-C的，在matlab上运行速度极慢，为了公平起见，我尽力去复现了BinSSA的算法，在论文给出的instance上，都达到了最优解，但均值没有达到。代码放到结尾处了。
（如果是将现有Objective-c直接对照改为c++/c本不是是一件难事，但是作者给的代码与论文上内容不符…，所以还是自己看着论文重新写的）
Phase 3 experimental results
后续在我的算法上具体调整了一些细节和一些参数，跑了M*上的所有例子，并且对于1000x1000和2000x2000的例子，根据生成方法分别随机生成了9个和4个。
M*上的所有例子都跑了30个种子
Phase 4
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 4 process
1）随机解策略的调整：
生成10000个随机解策略：从性价比最好的n个仓库中随机挑选20个仓库开设。且设定阈值：n若小于100则n=总仓库数目，反之n=100
2）顾客分配策略调整：
若顾客与该轮固定的仓库的距离大于上一次最优解中分配得到的距离，则该顾客被选入cplex中计算。
3）非cplex部分仓库固定策略
之前的方法是cplex范围右边全部固定为开设，左边全部固定为关闭，现在左边仍然全部固定为关闭，但右边为随机开设一部分仓库。
4）cplex范围调整
cplex范围为k*当前最优解数目且设定大于10的阈值。且cplex范围每迭代10次随当前最优解变化一次。
5）更新群体的参数调整
更新群体的分数由两部分影响：该解的海明距离和最优值。参数是用来调整两部分的权重。之前的策略会让该权重随迭代次数增加变化，现直接固定为定值。
Phase 4 conclusion
该阶段一直在调整各种小细节，挺折磨人的，看到糟糕的效果很难受，但这个阶段有点麻木了，遇到问题总是喜欢各种调一些小参数，但可能在当前的种子下效果好，换一个种子效果又变差了。
Phase 4仍在持续更新中…后续会将我的代码放上来。
Code
BinSSA
复现 A binary social spider algorithm for uncapacitated facility location problem 的代码如下
//
// main.cpp
// BinSSA
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string.h>
#include <time.h>
#include <vector>
#include <math.h>
#include <algorithm>
#include <assert.h>
#include<stdbool.h>
#include<float.h>
using 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 cost
for(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;i<N;i++) {
distances[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; i<whNum; i++)
if (warehouse[i] == 1)
cost = cost + f[i];
for(i=0; i<cusNum; i++)
{
min_cost = DBL_MAX;
for(j=0; j<whNum; j++)
{
if (warehouse[j] == 1 && dis[i][j]<min_cost)
min_cost = dis[i][j];
}
cost = cost + min_cost;
}
return cost;
}
void computeVibration (int l, double sigma_bar) {
int j, k;
double V;
double max = target_max;
int target = -1;
int flag = 1;
for(j=0;j<N;j++) {
flag = 1;
for(int i=0;i<whNum;i++)
if(x_spiders[l][i] != target_vibr[i]) flag = 0;
if(flag == 1) continue;
V = vibration[j] * exp((-distances[l][j]/(sigma_bar * ra))+1e-5);
if(V > max) {
target = j;
max = V;
}
}
if(target != -1) {
target_max = max;
inactive_deg = 0;
for(j=0;j<whNum;j++) {
target_vibr[j] = x_spiders[target][j];
}
} else{
inactive_deg++;
}
}
void SpiderWalk (int l) {
int i , j, r;
double temp;
//dimension_mask update
p_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;j<whNum;j++) {
dimension_mask[j] = rand() / double(RAND_MAX) < p_m;
}
}
double p_following;
for(j=0;j<whNum;j++) {
x_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 S2
if(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;i<whNum;i++) {
r = rand() % N;
// cout<<x_spiders[l][i] <<" "<<x_spiders[r][i]<<endl;
if(x_spiders[l][i] == 1) nn_1++;
else nn_0++;
if(x_spiders[l][i] == 1 && x_spiders[r][i] == 1) m_11++;
else if(x_spiders[l][i] == 1 && x_spiders[r][i] == 0) m_10++;
else if(x_spiders[l][i] == 0 && x_spiders[r][i] == 1) m_01++;
}
// cout<<m_11<<" "<<m_10<<" "<<m_01<<" "<<nn_0<<" "<<nn_1<<" "<<phi<<endl;
// exit(0);
double value, value_min = 999999;
for(M_10 = 0; M_10 <= nn_0; M_10++) {
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<<M_10<<" "<<M_01<<" "<<M_11<<" "<< value <<endl;
if(value < value_min) {
value_min = value;
M01 = M_01;
M10 = M_10;
M11 = M_11;
}
}
}
for(i=0;i<whNum;i++) {
if(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--;
}
else
candidate_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;j<whNum;j++) {
r_1 = (double)(rand()) /RAND_MAX;
if(r_1 < ST) {
if(dimension_mask[j] == x_spiders[r_2][j]) {//same
if(x_spiders[l][j] == 1)//1
candidate_spider[l][j] = 1;
else//0
candidate_spider[l][j] = 0;
}
else {//diff
if(x_spiders[l][j] == 1)//1
candidate_spider[l][j] = 0;
else//0
candidate_spider[l][j] = 1;
}
}
else
candidate_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"<<endl;
}
for(i=0;i<record;i++) {
outfile << history_obj[i]<<" ";
}
outfile<<endl;
for(i=0;i<record;i++) {
outfile << history_time[i]<<" ";
}
outfile<<endl;
outfile.close();
}
bool cmp (double X[], double Y[]) {
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;i<N;i++) x_memory[i] = new int[whNum];
double sigma_bar = 0;
double *sigma, *x_spider_bar;
sigma = new double[whNum];
x_spider_bar = new double [whNum];
int distance;
//Initialization
for(i=0;i<N;i++) {
for(j=0;j<whNum;j++) {
x_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;i<N;i++) {
obj_memory[i] = UFLP(x_spiders[i]);
}
for(j=0;j<whNum;j++) {
dimension_mask[j] = rand() % 2;
target_vibr[j] = rand() % 2;
}
//iteraion begin
while(time_cost < time_max) {
// cout<<"iter: "<<iter<<endl;
//fitness evaluation
for(i=0;i<N;i++) {
obj_spiders[i] = UFLP(x_spiders[i]);
if(obj_spiders[i] < obj_memory[i])
{
obj_memory[i] = obj_spiders[i];
for(j=0; j<whNum; j++)
x_memory[i][j] = x_spiders[i][j];
}
}
if (iter>1) {
for(i=0;i<N;i++) {
if (obj_spiders[i] > candidate_memory[i]){
for(j=0;j<whNum;j++){
x_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; j<whNum; j++)
x_memory[i][j] = x_spiders[i][j];
}
if (obj_memory[i] < global_best) {
global_best = obj_memory[i];
// cout<<global_best<<" "<<iter << endl;
best_time = (double) (1.0 * (clock() - start_time) / CLOCKS_PER_SEC);
history_obj[record] = obj_memory[i];
history_time[record++] = best_time;
}
}
}
}
for(i=0;i<N;i++) {
for(j=i; j<N;j++) {
distances[i][j] = 0;
distances[j][i] = distances[i][j];
}
}
//mahattan distance
for(i=0;i<N;i++) {
for(j=i+1; j<N;j++) {
distances[i][j] = 0;
for(int k =0;k<whNum;k++) {
distances[i][j] += abs(x_spiders[i][k] - x_spiders[j][k]);
}
distances[j][i] = distances[i][j];
}
}
for(i=0;i<N;i++) {
vibration[i] = log(1.0/(obj_spiders[i] - C) + 1 );
}
//compute sigma_bar
for(i=0;i<whNum;i++) {
x_spider_bar[i] = 0;
for(j=0;j<N;j++)
x_spider_bar[i] += x_spiders[j][i];
x_spider_bar[i] /= N;
}
sigma_bar = 0;
for(i=0;i<whNum;i++) {
sigma[i] = 0;
for(j=0;j<N;j++) {
sigma[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; i<N; i++)
computeVibration(i, sigma_bar);
for(i=0; i<N; i++)
SpiderWalk(i);
//candidate
phi = phi_max - (phi_max - phi_min) * iter / MAX_ITER ;
for(i=0; i<N; i++) {
r = 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<<endl;
// cout<<"BEST: " << global_best <<endl;
// cout<<"time is " << best_time << endl;
// cout<<"cost time: "<< (double) (1.0 * (clock() - start_time) / CLOCKS_PER_SEC)<<endl;
DeleteMemory();
}
exact algorithm 1
Phase 2-1 version 1的精确算法
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#include<string.h>
#include<stdbool.h>
#include<float.h>
#include <iostream>
#include <algorithm>
#include <fstream>
using namespace std;
string File_Name;
int whNum;
int cusNum;
double *f;
double **dis;
double postage_min = 0;
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 cost
for(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();
}
double ComputeSolution (int *warehouse) {
int i, j;
double distance_min;
double value = 0;
for(i=0;i<cusNum;i++) {
distance_min = 999999;
for(j=0;j<whNum;j++) {
if(warehouse[j] == 1) {
if(dis[i][j] < distance_min) {
distance_min = dis[i][j];
}
}
}
value += distance_min;
}
return value;
}
bool Cmp (double X[], double Y[]) {
return X[0] < Y[0];
}
int main() {
File_Name = "/Users/jiangtongjun/Desktop/UFLP/M/R/MR1";
DataInput();
//compute postage_min
int i, j;
double distance;
double distance_min;
double postage;
double obj;
double obj_best;
double postage_best;
int *best_solution;
best_solution = new int[whNum];
for(i=0;i<cusNum;i++) {
distance_min = 999999;
for(j=0;j<whNum;j++) {
if(dis[i][j] < distance_min) {
distance_min = dis[i][j];
}
}
postage_min += distance_min;
}
//C(1000,5)
int flag;
int *warehouse;
warehouse = new int [whNum];
for(i=0;i<whNum;i++) warehouse[i] = 0;
int *index;
index = new int [10];
int last_index;
int last_num;
int k = 5;
int offset;
double fix_best;
double fix;
double fix_min;
double *wh_rank_f;
wh_rank_f = new double [whNum];
for(i=0;i<whNum;i++) {
wh_rank_f[i] = f[i];
}
sort(wh_rank_f,wh_rank_f + whNum);
while(1) {
flag = 0;
fix = 0; //固定k个仓库时的仓库开设费用
obj_best = 999999;
offset = 0;
for(i=0;i<k;i++) {
index[i] = offset;
offset++;
}//初始化固定k个仓库时的仓库序号
for(i=0;i<pow(whNum, k);i++) {
for(i=0;i<whNum;i++) warehouse[i] = 0;
fix = 0;
for(j=0;j<k;j++) {
warehouse[index[j]] = 1;
}
//compute Solution;
postage = ComputeSolution(warehouse);//计算当前开设情况的运费
for(j=0;j<whNum;j++) {
if(warehouse[j] == 1)
fix += f[j];
}
obj = postage + fix;
cout<<obj<<": ";
for(j=0;j<whNum;j++) {
if(warehouse[j] == 1)
cout<<j<<" ";
}
cout<<endl;
if(obj < obj_best) {
postage_best = postage;
fix_best = fix;
obj_best = obj;
for(j=0;j<whNum;j++) {
best_solution[j] = warehouse[j];
}
}
//[996,997,998,999]
//[1,2,3,4]
if(index[k-1] == whNum - 1) {
last_num = 1;
last_index = k - 1;
for(j=k-2;j>=0;j--) {
if(index[j] == whNum - 1 - last_num) {
last_index = j;
}
last_num++;
}
if(last_index == 0) {
flag= 1;
}
if(flag == 0) {
index[last_index-1]++;
for(j=last_index-1;j<k-1;j++) {
index[j+1] = index[j] + 1;
}
}
}
else {
index[k-1]++;
}
if(flag == 1) break;
}
//judge
fix = 0;
for(i=0;i<k+1;i++)
fix += wh_rank_f[i];
fix_min = fix - fix_best;
if(fix_min > postage_min - postage_min) {
cout<<k+1<<"is feasible"<<endl;
k++;
}
else {
cout<<obj_best<<endl;
for(i=0;i<whNum;i++) {
if(best_solution[i]==1)
cout<<i<<" ";
}
break;
}
}
}
exact algorithm 2
Phase 2-3 version2的精确算法
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#include<string.h>
#include<stdbool.h>
#include<float.h>
#include <iostream>
#include <algorithm>
#include <fstream>
#include<ilcplex/ilocplex.h>
using namespace std;
string File_Name;
int whNum, cusNum;
double *f;
double **dis;
float global_best = 999999;
int *best_solution;
typedef IloArray<IloNumVarArray> NumVarMatrix;
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 cost
for(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 cplex(int k) {
int i, j;
IloEnv env;
IloModel model(env);
IloNumArray y( env, whNum);
IloNumVarArray open(env, whNum, 0, 1, ILOINT);
IloNumArray open2( env, whNum);
NumVarMatrix x(env, cusNum);
for (i = 0; i < cusNum; i++) { //x[i][j] 顾客i是否选择仓库j
x[i] = IloNumVarArray(env, whNum, 0, 1, ILOINT);
}
IloExpr obj_value(env);
IloExpr wh_sum(env, ILOINT);
for(i = 0; i < whNum; i++) {//仓库i 与 顾客j
obj_value += f[i] * open[i];
wh_sum += open[i];
for(j = 0; j < cusNum; j++) {
obj_value += dis[j][i] * x[j][i];
model.add(x[j][i] <= open[i]);
}
}
model.add(wh_sum == k);//仓库个数限制
//cout << wh_sum;
for(i = 0; i < cusNum; i++) {
IloExpr x_num(env);
for(j = 0; j < whNum; j++) {
x_num += x[i][j];
}
model.add(x_num == 1);
}
// model.add(obj_value < global_best);
model.add(IloMinimize(env, obj_value));
IloCplex cplex(model);
cplex.setOut(env.getNullStream());
cplex.setParam(IloCplex::Param::Threads, 1);
cplex.solve();
cplex.getValues(y, open);
cout<<"wh: ";
for(i = 0; i < whNum; i++) {
best_solution[i] = y[i];
if(y[i] == 1)
cout<<i<<" ";
}
cout<<endl;
if(cplex.getObjValue() < global_best)
global_best = cplex.getObjValue();
// cout << cplex.getObjValue();
cplex.clear();
cplex.end();
y.clear();
y.end();
open.clear();
open.end();
x.clear();
x.end();
obj_value.clear();
obj_value.end();
model.end();
env.end();
}
double ComputeSolution (int *warehouse) {
int i, j;
double distance_min;
double value = 0;
for(i=0;i<cusNum;i++) {
distance_min = 999999;
for(j=0;j<whNum;j++) {
if(warehouse[j] == 1) {
if(dis[i][j] < distance_min) {
distance_min = dis[i][j];
}
}
}
value += distance_min;
}
return value;
}
int main() {
int i, j;
File_Name = "../M/R/MR1";
DataInput();
double fix_min, fix;
double postage_min = 0, postage_fix;//最低运费
double distance_min;
double *wh_rank_f;
wh_rank_f = new double [whNum];
for(i=0;i<whNum;i++) {
wh_rank_f[i] = f[i];
}
// for(i=0;i<whNum;i++)
// cout<<i<<": "<<wh_rank_f[i]<<endl;
sort(wh_rank_f,wh_rank_f + whNum);
// cout<<"++++++++++++++"<<endl;
// for(i=0;i<whNum;i++)
// cout<<i<<": "<<wh_rank_f[i]<<endl;
//initiliazation
best_solution = new int [whNum];
for(i=0;i<cusNum;i++) {
distance_min = 999999;
for(j=0;j<whNum;j++) {
if(dis[i][j] < distance_min) {
distance_min = dis[i][j];
}
}
postage_min += distance_min;
}
//仓库的个数固定 限制
//最优解的值 > 上一次cplex的值
//
int warehouse = 1;
while(1) {
// Cplex(warehouse);
cplex(warehouse + 1);
fix = 0;
for(i=0;i<whNum;i++)
if(best_solution[i] == 1)
fix += f[i];
//judge
//1->2、2->3、、、warehouse->warehouse+1 至少增加仓库开设费用
fix_min = 0;
for(i=0;i<warehouse+1;i++) {
fix_min += wh_rank_f[i];
}
fix_min = fix_min - fix;
postage_fix = ComputeSolution(best_solution);
cout<<"fix_min: "<< fix_min<<endl;
cout<<"postage_fix: "<< postage_fix <<endl;
cout<<"postage_min: "<< postage_min <<endl;
//至少增加的开设费用 > 最多能节省的运费
if(fix_min > postage_fix - postage_min) {
cout<<global_best<<": "<<endl;
for(j=0;j<whNum;j++) {
if(best_solution[j]==1)
cout<<j<<" ";
}
break;
}
else { //至少增加的开设费用 < 最多能节省的运费
cout<<global_best<<endl;
cout<<warehouse + 1<<"is feasible"<<endl;
cout<<endl;
warehouse++;
}
}
}
最后，感谢Dr.He，也感谢Nenu李宏博老师的指导。
扫一扫
1193
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
