思路
分别使用模拟退火算法与禁忌搜索算法求解。
因为题目的精确求解算法时间复杂度很高,所以使用元启发式算法来求解。将原题目分为两个问题:1. 开放哪些工厂; 2. 用户分配给哪些工厂。 这样,就可以用一串01串来表示开放工厂列表,用一个整数串表示用户的工厂分配。
因为工厂开放之间彼此没有什么规律,因此原则上在邻域搜索的时候,除非出现了大量的无解情况,否则不应该新加入工厂。可以有两种探索邻域的方法:一种是将工厂A的顾客A移交给工厂B,一种是将工厂A的顾客A与工厂B的顾客B相交换。其中后者一般用于工厂B的空间不足以放下顾客A时。
同时考虑到很难验证是否有解,同时合法的情况(剩余容量大于0)能从非法的情况(剩余容量小于0)中产生,所以允许非法的情况出现在迭代过程中,只是要有一定的惩罚值。
算法架构
贪婪自适应搜索算法
用于两种算法中产生初始解。这种算法将会一直开放工厂,直到工厂的总容量大于用户的总需求,并在此过程中将一些用户指定到工厂中。在工厂开放完毕之后,将需求最大的顾客指定到剩余容量最多的工厂中(无论是否合法),直到所有顾客都指定完毕。
工厂开放的时候,通过贪婪函数对每个未开放的工厂进行求值,其值为开放这个工厂的代价与所有能分配到该工厂的顾客分配到该工厂的总代价除以分配的顾客数量(分配的顾客顺序按照其对该工厂的代价升序排序)。这样求出的值中选取最小值,将140%*最小值作为阈值,在阈值之下的工厂纳入到随机选择列表中,从中随机选出要开放的工厂。
模拟退火算法
通过贪婪自适应算法产生若干个初始解(目前是每个问题500个初始解),然后用模拟退火对这些初始解进行第二阶段的邻域搜索。
参数
初始温度:400,终止温度:1,温度变化公式:K = 0.97*K,在100次没有更新最优解后进行跳出。
模拟退火的邻域选取中,如果不存在不合法解,则随机选取两个工厂,随机选择其中一个工厂的用户,如果这个用户能够移动到第二个工厂(第二个工厂的剩余容量大于第一个工厂的用户需求),则进行移动;不然,再次随机选取第二个工厂的用户,两者进行交换(无论是否合法)。
如果存在不合法解,则从剩余容量小于0的工厂中按照超出容量加权随机,选择一个工厂;再从剩余容量大于0的工厂中,按照剩余容量加权随机,选择一个工厂。用户移动同上。
禁忌搜索算法
同样使用贪婪自适应算法产生若干个初始解,然后搜索第二阶段,随机产生若干个邻域,每次从这些邻域中选择满足禁忌条件的最小值,并将其加入到禁忌条件中。
参数
禁忌表长度:用户数量的四分之一,200次没有更新最优解则跳出,每次随机产生50个邻域
禁忌表为pair<int,int>
的列表,禁止用户再次分配到指定工厂。
禁忌表的邻域搜索与模拟退火略有不同。
对于工厂的选择,在没有不合法解的时候,按照剩余容量加权随机选择第二个工厂,第一个工厂完全随机。如果存在不合法解的话,同上。
对于用户选择,会穷举所有可能的用户情况。考虑到两个工厂的容量有限,因此可能性不会太多。在其中如果满足禁忌列表要求,或者不满足禁忌列表要求但是比最优解还小的解,会被纳入考量中。
源代码
源代码分为主程序main.cpp与头文件Solution.h
主程序
main.cpp
#include <iostream>
#include "Solution.h"
#include <ctime>
#include <string>
using namespace std;
#define OUTPUT 1 //output to file?
//#define DEBUG
double solveUsingSA(string filename)
{
Solution s(filename);
vector<Solve> initSet;
int factoryNum = s.getFactorySize();
double times;
times = 500;
for(int i = 0;i<times;i++)
{
initSet.push_back(s.initializeSolution());
}
Solve ss = s.SA(initSet);
#ifdef DEBUG
cout<<"best solution"<<endl;
ss.print();
#endif
if(!ss.fitness)//如果解不合法,则尝试修复
{
s.repair(ss);
#ifdef DEBUG
cout<<"after repair"<<endl;
ss.print();
#endif
}
while(!ss.fitness)//如果仍然不符合限制,则继续修复
{
#ifdef DEBUG
cout<<"not fit"<<endl;
#endif
//如果不符合限制,则尝试打开工厂
//对于每个给定的初始解,随机打开一个关闭的工厂
for(int i = 0;i<times;i++)
{
Solve& solve = initSet[i];
vector<int> factory;
for(int j = 0;j<solve.openList.size();j++)//因为如果有多于一次开放工厂的可能,所以这部分要放在循环内
{
if(!solve.openList[j])
{
factory.push_back(j);
}
}
int index = rand()%factory.size();
index = factory[index];
solve.openList[index] = true;
solve.restCapacity[index] = s.facCapacity[index];
}
ss = s.SA(initSet);
#ifdef DEBUG
cout<<"\n\n"<<endl;
ss.print();
cout<<"\n\n"<<endl;
#endif
}
ss.value = s.judgeValue(ss);//因为可能修改,所以进行更新
#ifdef OUTPUT
s.outputTofile("sa",ss);
#endif
return ss.value;
}
void outputSASolve()
{
ofstream fout("ans_sa_total.txt",ios::app);
vector<string> filenameSet;
for(int i = 1;i<=71;i++)
{
filenameSet.push_back("p"+to_string(i));
}
for(int i = 0;i<filenameSet.size();i++)
{
time_t begining,ending;
begining = time(NULL);
double ans = solveUsingSA(filenameSet[i]);
ending = time(NULL);
#ifdef OUTPUT
fout<<filenameSet[i]<<"\t"<<ans<<"\t"<<difftime(ending,begining)<<endl;
#endif
cout<<filenameSet[i]<<"\t"<<ans<<"\t"<<difftime(ending,begining)<<endl;
}
fout.close();
}
double solveUsingTS(string filename)
{
Solution s(filename);
vector<Solve> initSet;
for(int i = 0;i<100;i++)
{
initSet.push_back(s.initializeSolution());
}
Solve ss = s.TabuSearch(initSet);
if(!ss.fitness)
{
s.repair(ss);
}
ss.value = s.TabuSearchJudge(ss);
#ifdef DEBUG
cout<<ss.value<<endl;
ss.print();
#endif
#ifdef OUTPUT
s.outputTofile("ts",ss);
#endif
return ss.value;
}
void outputTSSolve()
{
ofstream fout("ans_ts_total.txt",ios::app);
vector<string> filenameSet;
for(int i = 1;i<=71;i++)
{
filenameSet.push_back("p"+to_string(i));
}
for(int i = 0;i<filenameSet.size();i++)
{
time_t begining,ending;
begining = time(NULL);
double ans = solveUsingTS(filenameSet[i]);
ending = time(NULL);
#ifdef OUTPUT
fout<<filenameSet[i]<<"\t"<<ans<<"\t"<<difftime(ending,begining)<<endl;
#endif
cout<<filenameSet[i]<<"\t"<<ans<<"\t"<<difftime(ending,begining)<<endl;
}
fout.close();
}
int main(void)
{
outputSASolve();
outputTSSolve();
return 0;
}
Solution.h
#ifndef SOLUTION_H
#define SOLUTION_H
#include <string>
#include <vector>
#include <map>
#include <algorithm>
#include <iostream>
#include <ctime>
#include <cstdlib>
#include <utility>
#include <fstream>
#include <cmath>
//for Tabu Search
/*
用来记录禁忌搜索的操作(两个工厂之间交换顾客,或者一个工厂给予另外一个工厂顾客)
并记录这样变动后的大小
*/
struct Oper
{
int j1,j2;
int i1,i2;
int weight;
Oper(int j1,int j2,int i1,int i2,int w)
{
this->j1 = j1;
this->j2 = j2;
this->i1 = i1;
this->i2 = i2;
weight = w;
}
};
#define MAX 0xFFFFFFFF
using namespace std;
struct Solve
{
vector<bool> openList; //开放工厂列表,一系列01串,用来表示工厂是否开着
vector<int> assignmentList; // 顾客的分配列表,表示顾客分配给了哪个工厂
//工厂相关
vector<int> restCapacity;
vector<int> useCapacity;
bool fitness;
double value;
Solve(int n,int m)
{
openList = vector<bool>(n,0);
assignmentList = vector<int>(m,-1);
restCapacity = vector<int>(n,0);
useCapacity = vector<int>(n,0);
fitness = true;
}
Solve()
{
fitness = true;
}
Solve(const Solve& other)
{
openList = other.openList;
assignmentList = other.assignmentList;
restCapacity = other.restCapacity;
useCapacity = other.useCapacity;
fitness = other.fitness;
value = other.value;
}
//如果工厂j的容量在操作后小于0,返回false,不然返回true
bool factoryUseCapacity(int j,int value)
{
restCapacity[j] -= value;
useCapacity[j] += value;
if(restCapacity[j]<0)
{
fitness = false;
}
return fitness;
}
//如果工厂j的容量在操作后小于0,返回false,不然返回true
/*
如果某个工厂过去流量为负,突然见为正了,则可能从一个不合法解变为合法解
*/
bool factoryFreeCapacity(int j,int value)
{
bool change = false;
if(restCapacity[j]<0 && restCapacity[j] + value>=0)
change = true;
restCapacity[j] += value;
useCapacity[j] -= value;
if(change)
{
fitness = true;
for(int i = 0 ;i<restCapacity.size();i++)
{
if(restCapacity[i]<0)
{
fitness = false;
break;
}
}
}
return fitness;
}
void print()
{
cout<<"Open List"<<endl;
for(int i = 0;i<openList.size();i++)
{
cout<<openList[i]<<" ";
}
cout<<endl;
cout<<"assignment list"<<endl;
for(int i = 0;i<assignmentList.size();i++)
{
cout<<assignmentList[i]<<" ";
}
cout<<endl;
cout<<"use"<<endl;
for(int i = 0;i<useCapacity.size();i++)
{
cout<<useCapacity[i]<<" ";
}
cout<<endl;
cout<<"rest"<<endl;
for(int i = 0;i<restCapacity.size();i++)
{
cout<<restCapacity[i]<<" ";
}
cout<<endl;
}
};
bool cmp1(const pair<int,int>& a,const pair<int,int>& b)
{
return a.second<b.second;
}
bool cmp2(const pair<int,int>& a,const pair<int,int>& b)
{
return a.second>b.second;
}
class Solution
{
public:
int clientNum,factoryNum; // 客户数目 与 工厂数目
vector<int> facCapacity; //各个工厂的容量
vector<int> facOpenCost; //各个工厂的花费
vector<int> cliDemand; //用户的需求
vector<vector<double>> assignmentCost;//外层为工厂号,内层为用户号
//贪婪估计函数,用来为贪婪自适应随机算法提供参考数据
/*
通过对没有分配且可以分配到该工厂(尽量填满)的顾客进行按照代价升序排序,以尽量小的代价获得更多的流量。
*/
double greedyFunction(int j,Solve& s)
{
if(s.openList[j])
{
return 0;
}
double totalCost = 0;
int num =