期末项目：PI-CFLP问题的求解

github地址

思路

分别使用模拟退火算法与禁忌搜索算法求解。

因为题目的精确求解算法时间复杂度很高，所以使用元启发式算法来求解。将原题目分为两个问题：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 = 0;
        vector<pair<int,int>> cost;//工厂j到各个用户i的标号与代价
        for(int i = 0;i<clientNum;i++)
        {
            if(s.assignmentList[i]==-1)
            {
                cost.push_back(make_pair(i,assignmentCost[j][i]));
            }
        }

        sort(cost.begin(),cost.end(),cmp1);//升序排列代价

        int restCapacity = s.restCapacity[j];
        for(int i = 0;i<cost.size();i++)
        {
            int client = cost[i].first;
            int asCost = cost[i].second;
            if(restCapacity - cliDemand[client]<0)//到达底线，取消
            {
                break;
            }
            else
            {
                restCapacity -= cliDemand[client];
            }
            num++;
            totalCost += asCost;
        }

        return (totalCost + facOpenCost[j])/num;
    }

    //局部搜索估值函数
    double judgeValue(const Solve& s,bool switcher = true)
    {
        double totalCost = 0;
        for(int i = 0;i<factoryNum;i++)
        {
            if(s.openList[i])
            {
                totalCost += facOpenCost[i];
            }

            if(s.restCapacity[i]<0)
            {
                if(switcher)
                    totalCost += (double)abs(s.restCapacity[i])/facCapacity[i]*facOpenCost[i];
                else
                    totalCost += MAX;
            }
        }

        for(int i = 0;i<clientNum;i++)
        {
            totalCost += assignmentCost[s.assignmentList[i]][i];
        }
        return totalCost;
    }
    //根据给定解返回其随机邻域中的一个，并进行赋值
    /*
    邻域操作：
        如果解没有不合法的地方：
            随机选取两个工厂
                随机选取第一个工厂的客户i，如果第二个工厂有空位，则丢过去。不然，交换。

        如果解有不合法的地方：
            按照剩余概率加权与超出概率加权选择两个工厂，优先将超出工厂的流量移除，不然则交换。
    */
    Solve generateNeighbourRandom(const Solve& origin)
    {
        Solve s = origin;
        int j1,j2;
        if(s.fitness)//如果该解是合法解
        {
            j1 = rand()%factoryNum,j2 = rand()%factoryNum;
            while(j1==j2 || !s.openList[j1] || !s.openList[j2])
            {
                j1 = rand()%factoryNum;
                j2 = rand()%factoryNum;
            }
        }
        else
        {
            //如果存在不合法解，则将不合法的工厂的超出流量进行加权后随机选择，合法工厂的剩余流量加权后随机选择，这样子可以更可能选出尽量剩余多的工厂来接受流量。
            int totalRest = 0;
            int totalOut = 0;
            vector<pair<int,int>> restFactory;
            vector<pair<int,int>> outFactory;
            for(int i = 0;i<factoryNum;i++)
            {
                if(s.restCapacity[i]>0 && s.openList[i])
                {
                    totalRest += s.restCapacity[i];
                    restFactory.push_back(make_pair(i,totalRest));
                }
                else if(s.restCapacity[i]<0 && s.openList[i])
                {
                    totalOut += abs(s.restCapacity[i]);
                    outFactory.push_back(make_pair(i,totalOut));
                }
            }

            if(totalOut == 0)//可能出现错漏
            {
                s.fitness = true;
                return generateNeighbourRandom(s);
            }
            //找到大于指针的最小值
            int index1 = rand()%totalOut;
            int index2 = rand()%totalRest;
            j1 = outFactory[outFactory.size()-1].first;
            j2 = restFactory[restFactory.size()-1].first;
            int j1Judge = totalOut,j2Judge = totalRest;
            for(int i = 0;i<outFactory.size();i++)
            {
                if(outFactory[i].second>index1 && outFactory[i].second<j1Judge)
                {
                    j1Judge = outFactory[i].second;
                    j1 = outFactory[i].first;
                }
            }
            for(int i = 0;i<restFactory.size();i++)
            {
                if(restFactory[i].second>index2 && restFactory[i].second<j2Judge)
                {
                    j2Judge = restFactory[i].second;
                    j2 = restFactory[i].first;
                }
            }
        }
        //寻找要交换的顾客
        vector<int> clientI,clientJ;
        for(int i = 0 ;i<clientNum;i++)
        {
            if(s.assignmentList[i] == j1)
            {
                clientI.push_back(i);
            }
            else if(s.assignmentList[i] == j2)
            {
                clientJ.push_back(i);
            }
        }

        if(clientI.size()==0)//有可能出现工厂没有分配人的情况，这时候随机开始下一次
        {
            return generateNeighbourRandom(origin);
        }
        //随机选择一个顾客，进行移动或者交换
        int index1 = rand()%clientI.size();

        int i1 = clientI[index1];
        if(cliDemand[i1] < s.restCapacity[j2])//如果工厂j2能够接受i1，则移动
        {
            s.factoryFreeCapacity(j1,cliDemand[i1]);
            s.factoryUseCapacity(j2,cliDemand[i1]);
            s.assignmentList[i1] = j2;
        }
        else//不然，则进行交换
        {
            int index2 = rand()%clientJ.size();
            int i2 = clientJ[index2];
            s.factoryFreeCapacity(j1,cliDemand[i1]);
            s.factoryUseCapacity(j2,cliDemand[i1]);
            s.factoryFreeCapacity(j2,cliDemand[i2]);
            s.factoryUseCapacity(j1,cliDemand[i2]);
            s.assignmentList[i1] = j2;
            s.assignmentList[i2] = j1;
        }
        s.value = judgeValue(s);

        return s;
    }
    //产生操作，用于禁忌搜索
    Oper generateNeighbourOper(const Solve& origin,vector<pair<int,int>>& tabuList,int minNum)
    {
        Solve s = origin;
        int j1,j2;
        if(s.fitness)//如果该解是合法解，则用容量加权随机求j2（因为j2是接受者，显然，一个容量更大的接收者会更好）
        {
            //求加权平均值
            int capaSum = 0;
            int pointer = 0;
            vector<pair<int,int>> capacity;
            for(int i = 0;i<s.openList.size();i++)
            {
                if(!s.openList[i])
                    continue;

                capaSum += s.restCapacity[i];
                capacity.push_back(make_pair(i,capaSum));
            }

            pointer = rand()%capaSum;
            j2 = capacity[capacity.size()-1].first;
            int j2Judge = capaSum;
            for(int i = 0;i<capacity.size();i++)
            {
                if(capacity[i].second>pointer && capacity[i].second<j2Judge)
                {
                    j2Judge = capacity[i].second;
                    j2 = capacity[i].first;
                }
            }

            j1 = rand()%factoryNum;
            while(!s.openList[j1] || j1==j2)
            {
                j1 = rand()%factoryNum;
            }
        }
        else
        {
            //如果存在不合法解，则将不合法的工厂的超出流量进行加权后随机选择，合法工厂的剩余流量加权后随机选择，这样子可以更可能选出尽量剩余多的工厂来接受流量。
            int totalRest = 0;
            int totalOut = 0;
            vector<pair<int,int>> restFactory;
            vector<pair<int,int>> outFactory;
            for(int i = 0;i<factoryNum;i++)
            {
                if(s.restCapacity[i]>0 && s.openList[i])
                {
                    totalRest += s.restCapacity[i];
                    restFactory.push_back(make_pair(i,totalRest));
                }
                else if(s.restCapacity[i]<0 && s.openList[i])
                {
                    totalOut += abs(s.restCapacity[i]);
                    outFactory.push_back(make_pair(i,totalOut));
                }
            }

            if(totalOut == 0)//可能出现错漏
            {
                s.fitness = true;
                return generateNeighbourOper(s,tabuList,minNum);
            }
            //找到大于指针的最小值
            int index1 = rand()%totalOut;
            int index2 = rand()%totalRest;
            j1 = outFactory[outFactory.size()-1].first;
            j2 = restFactory[restFactory.size()-1].first;
            int j1Judge = totalOut,j2Judge = totalRest;
            for(int i = 0;i<outFactory.size();i++)
            {
                if(outFactory[i].second>index1 && outFactory[i].second<j1Judge)
                {
                    j1Judge = outFactory[i].second;
                    j1 = outFactory[i].first;
                }
            }
            for(int i = 0;i<restFactory.size();i++)
            {
                if(restFactory[i].second>index2 && restFactory[i].second<j2Judge)
                {
                    j2Judge = restFactory[i].second;
                    j2 = restFactory[i].first;
                }
            }
        }
        //寻找要交换的顾客
        vector<int> clientI,clientJ;
        for(int i = 0 ;i<clientNum;i++)
        {
            if(s.assignmentList[i] == j1)
            {
                clientI.push_back(i);
            }
            else if(s.assignmentList[i] == j2)
            {
                clientJ.push_back(i);
            }
        }

        if(clientI.size()==0)//有可能出现工厂没有分配人的情况，这时候随机开始下一次
        {
            return generateNeighbourOper(origin,tabuList,minNum);
        }
        //预先遍历所有客户，找到最优解。注：要考虑到禁忌列表的因素
        int index1 = rand()%clientI.size();
        int minW = MAX;
        Oper opera(-1,-1,-1,-1,-1);
        for(int i = 0;i<clientI.size();i++)
        {
            int i1 = clientI[i];
            if(cliDemand[i1] < s.restCapacity[j2])
            {
                //该操作需要i1 j2不在禁忌列表中
                bool finding = false;
                for(int i = 0;i<tabuList.size();i++)
                {
                    if(tabuList[i].first==i1&&tabuList[i].second==j2)
                    {
                        finding = true;
                        break;
                    }
                }
                if(assignmentCost[j2][i1] - assignmentCost[j1][i1] < minW && !finding)
                {
                    minW = assignmentCost[j2][i1] - assignmentCost[j1][i1];
                    opera = Oper(j1,j2,i1,-1,minW);
                }
                else if(s.value + assignmentCost[j2][i1] - assignmentCost[j1][i1] < minNum && finding && assignmentCost[j2][i1] - assignmentCost[j1][i1] < minW)
                {
                    minW = assignmentCost[j2][i1] - assignmentCost[j1][i1];
                    opera = Oper(j1,j2,i1,-1,minW);
                }
            }
            else
            {
                for(int j = 0;j<clientJ.size();j++)
                {
                    int i2 = clientJ[j];

                    bool finding = false;
                    for(int i = 0;i<tabuList.size();i++)
                    {
                        if(tabuList[i].first==i1&&tabuList[i].second==j2 || tabuList[i].first==i2&&tabuList[i].second==j1)
                        {
                            finding = true;
                            break;
                        }
                    }
                    if(assignmentCost[j2][i1] + assignmentCost[j1][i2] - assignmentCost[j1][i1] - assignmentCost[j2][i2] <minW && !finding)
                    {
                        minW = assignmentCost[j2][i1] + assignmentCost[j1][i2] - assignmentCost[j1][i1] - assignmentCost[j2][i2];
                        opera = Oper(j1,j2,i1,i2,minW);
                    }
                    else if(s.value + assignmentCost[j2][i1] + assignmentCost[j1][i2] - assignmentCost[j1][i1] - assignmentCost[j2][i2] < minNum && finding && assignmentCost[j2][i1] + assignmentCost[j1][i2] - assignmentCost[j1][i1] - assignmentCost[j2][i2] <minW)
                    {
                        minW = assignmentCost[j2][i1] + assignmentCost[j1][i2] - assignmentCost[j1][i1] - assignmentCost[j2][i2];
                        opera = Oper(j1,j2,i1,i2,minW);
                    }
                }
            }
        }

        return opera;
    }
public:
    //使用文件内容初始化
    Solution(string filename)
    {
        srand((unsigned)time(NULL));
        ifstream fin(filename);
        fin>>factoryNum>>clientNum;
        facCapacity = vector<int>(factoryNum,0);
        facOpenCost = vector<int>(factoryNum,0);
        cliDemand = vector<int>(clientNum,0);
        for(int i = 0;i<factoryNum;i++)
        {
            fin>>facCapacity[i]>>facOpenCost[i];
        }
        for(int i = 0;i<clientNum;i++)
        {
            double x;
            fin>>x;
            cliDemand[i] = x;
        }
        for(int i = 0;i<factoryNum;i++)
        {
            vector<double> cost(clientNum,0);
            for(int j = 0;j<clientNum;j++)
            {
                fin>>cost[j];
            }
            assignmentCost.push_back(cost);
        }
        //对于每个用户建立一个排序表，按照各个工厂的费用顺序排序

        fin.close();
    }

    //使用贪婪随机自适应算法产生初始解并返回
    Solve initializeSolution()
    {
        double alpha = 0.4;
        Solve s(factoryNum,clientNum);
        for(int j = 0;j<factoryNum;j++)
        {
            s.restCapacity[j] = facCapacity[j];
        }
        int totalDemand = 0;//尚未满足的用户需求
        for(int i = 0;i<clientNum;i++)
        {
            totalDemand += cliDemand[i];
        }
        int satisFiedDemand = 0;
        //对于所有的工厂，检查其是否能够满足用户需求
        for(int j = 0;j<factoryNum;j++)
        {
            //数组greedyCost保存所有工厂的贪婪判断值
            vector<double> greedyCost(factoryNum,0);
            double greedyMin = MAX;
            for(int i = 0;i<factoryNum;i++)
            {
                if(s.openList[i])
                {
                    greedyCost[i] = MAX;
                }
                else
                {
                    greedyCost[i] = greedyFunction(i,s);
                    if(greedyCost[i]<greedyMin)
                    {
                        greedyMin = greedyCost[i];
                    }
                }
            }
            //产生候选解集
            vector<int> RCL;
            double threshold = (1+alpha)*greedyMin;
            for(int i = 0;i<greedyCost.size();i++)
            {
                if(greedyCost[i]<threshold)
                {
                    RCL.push_back(i);
                }
            }
            int index = rand()%RCL.size();
            int openFactory = RCL[index];
            s.openList[openFactory] = true;

            //求出需要赋值的客户列表（对所有用户到该工厂的代价升序排序，从最小取起）
            vector<pair<int,int>> cost;//工厂j到各个用户i的标号与代价
            for(int i = 0;i<clientNum;i++)
            {
                if(s.assignmentList[i]==-1)
                {
                    cost.push_back(make_pair(i,assignmentCost[openFactory][i]));
                }
            }

            sort(cost.begin(),cost.end(),[](const pair<int,int>& a,const pair<int,int>& b)
            {
                return a.second<b.second;
            });
            //进行赋值，修改解
            for(int i = 0;i<cost.size();i++)
            {
                int client = cost[i].first;
                if(s.restCapacity[openFactory] - cliDemand[client]<0)//如果剩余空间不足，则退出
                {
                    break;
                }

                s.factoryUseCapacity(openFactory,cliDemand[client]);//加入顾客
                s.assignmentList[client] = openFactory;
            }
            //检查是否满足需求
            satisFiedDemand += s.useCapacity[openFactory];
            if(satisFiedDemand>=totalDemand)
            {
                break;
            }

        }
        //对于没有赋值的客户进行赋值。取没有赋值的客户按需求从大到小排列，有剩余容量的工厂从大到小排列，以此装入
        vector<pair<int,int>> restClient;
        vector<pair<int,int>> restFactory;
        for(int i = 0;i<clientNum;i++)
        {
            if(s.assignmentList[i]==-1)
            {
                restClient.push_back(make_pair(i,cliDemand[i]));
            }
        }
        for(int j = 0;j<clientNum;j++)
        {
            if(!s.openList[j])
            {
                restFactory.push_back(make_pair(j,facCapacity[j]));
            }
        }
        sort(restClient.begin(),restClient.end(),cmp1);
        sort(restFactory.begin(),restFactory.end(),cmp1);
        //对于剩余客户进行分配，每次分配后，重新排序
        while(restClient.size()>0)
        {
            int client = restClient[restClient.size()-1].first;
            int factory = restFactory[restFactory.size()-1].first;
            s.assignmentList[client] = factory;
            s.factoryUseCapacity(factory,cliDemand[client]);

            restClient.pop_back();
            restFactory[restFactory.size()-1].second -= cliDemand[client];

            sort(restClient.begin(),restClient.end(),cmp1);
            sort(restFactory.begin(),restFactory.end(),cmp1);
        }
        return s;
    }

    //根据产生的初始解进行模拟退火并返回
    Solve SA(vector<Solve> initSet)
    {
        Solve best;
        double bestAns = MAX;//保存最小的解
        for(int i = 0;i<initSet.size();i++)
        {
            Solve s = initSet[i];
            s.value = judgeValue(s);
            if(s.value<bestAns)
            {
                bestAns = s.value;
                best = s;
            }

            double temperature = 400;
            double endTemperature = 1;
            double ratio = 0.97;
            int endCounting = 0;

            while(temperature>endTemperature)
            {
                //随机产生邻域
                //如果是好解，则无条件接纳
                    //保存全局最优解
                //如果是差解，则有选择接纳
                Solve ss = generateNeighbourRandom(s);
                if(ss.value<s.value)
                {
                    s = ss;
                    if(s.value<bestAns)//最大值更新
                    {
                        bestAns = s.value;
                        best = s;
                    }
                    endCounting = 0;
                }
                else
                {
                    endCounting++;
                    double p = exp((s.value-ss.value)/temperature);
                    double pointer = (double)1.0*rand()/RAND_MAX;
                    if(p>pointer)
                    {
                        s = ss;
                    }
                }
                if(endCounting > 100)
                {
                    break;
                }

                temperature *= ratio;
            }
        }
        return best;
    }

    void repair(Solve& ss)
    {
        //假设有一个工厂因为某些原因引入了一个客户导致出现负数的情况
        int factory = 0;
        for(int i = 0;i<ss.restCapacity.size();i++)
        {
            if(ss.restCapacity[i]<0)
            {
                factory = ss.restCapacity[i];
            }
        }
        vector<int> customers;
        for(int i = 0;i<ss.assignmentList.size();i++)
        {
            if(ss.assignmentList[i]==factory)
            {
                customers.push_back(i);
            }
        }
        //对于每一个用户，检查移走它们的费用与可能性
        for(int i = 0;i<customers.size();i++)
        {
            int c = customers[i];
            int minCost = MAX;
            int orderFactory = -1;
            for(int j = 0;j<ss.openList.size();j++)
            {
                if(ss.openList[j])//如果开放，能放得下，还够小
                {
                    if(ss.restCapacity[j]>=cliDemand[c] && assignmentCost[j][c]<minCost)
                    {
                        minCost = assignmentCost[j][c];
                        orderFactory = j;
                    }
                }
            }
            //如果有工厂能够接纳，则进行移动
            if(orderFactory != -1)
            {
                ss.assignmentList[c] = orderFactory;
                ss.factoryFreeCapacity(factory,cliDemand[c]);
                ss.factoryUseCapacity(orderFactory,cliDemand[c]);
                ss.fitness = true;
                for(int j = 0 ;j < ss.restCapacity.size();j++)
                {
                    if(ss.restCapacity[j]<0)
                    {
                        ss.fitness = false;
                        break;
                    }
                }
                if(ss.fitness)
                    break;
            }

        }
    }

    int getFactorySize()
    {
        return factoryNum;
    }
    void outputTofile(string filename,Solve s)
    {
        string newFilename = \"ans_\"+filename+\".txt\";
        ofstream fout(newFilename,ios::app);
        fout<<s.value<<endl;
        for(int i = 0;i<s.openList.size();i++)
        {
            fout<<s.openList[i]<<\" \";
        }
        fout<<endl;
        for(int i = 0;i<s.assignmentList.size();i++)
        {
            fout<<s.assignmentList[i]<<\" \";
        }
        fout<<endl;
        fout.close();
    }

    void clear()
    {
        facCapacity.clear();
        facOpenCost.clear();
        cliDemand.clear();
        assignmentCost.clear();
    }

/*
    禁忌搜索用的判断函数（主要是之前的判断函数似乎效果不好，试用一下其他的）
*/
   int TabuSearchJudge(const Solve& s,bool debug = false)
   {
        int penalty = 100;//惩罚系数，每超出一单位容量惩罚50
        int totalCost = 0;
        int penaltyc = 0;
        for(int i = 0;i<factoryNum;i++)
        {
            if(s.openList[i])
            {
                totalCost += facOpenCost[i];
            }

            if(s.restCapacity[i]<0)
            {
                totalCost += penalty*abs(s.restCapacity[i]);
            }
        }

        for(int i = 0;i<clientNum;i++)
        {
            totalCost += assignmentCost[s.assignmentList[i]][i];
        }
        return totalCost;
   }


    //根据产生的初始解进行禁忌搜索并返回
    //用一个循环列表来作为禁忌搜索的表格

    /*
    禁忌搜索具体步骤：
        在邻域中找到一个合适的且最好的（符合禁忌列表要求）
        如果能够找到一个比最优解更好的解，可以特赦

        中期与长期暂时不做
        在不考虑到工厂开放费用的情况下，一个客户从工厂A到工厂B，得益于其差值。因为A是不可选的，所以B选择越小的，越能得到更优解。
        如果放不下客户A，可以考虑进行交换。
        上述两种情况，都应该优先从小的考虑。
    */

    Solve TabuSearch(vector<Solve>& initSet)
    {
        Solve best;
        int bestValue = MAX;
        srand((unsigned)time(NULL));
        vector<pair<int,int>> tabuList(clientNum/4,make_pair(-1,-1));//禁忌列表表示在这段范围内，客户i不入工厂j
        int tabuPoint = 0;
        int tabuEnd = 0;
        //tabuList为循环队列，tabuPoint==tabuEnd时为空，其余情况下由tabuPoint -> tabuEnd，到达tabuEnd则停止。如果队列满，则有tabuPoint = tabuEnd + 1

        //initSet为外层构建好的初始解，大部分可行（不排除不可行的情况）
        /*
        当没有达到停止要求的时候
            创建移动的候选集合，运用局部搜索的思想，各个工厂间均运用贪心的思想（即用户会优先选择代价更小的工厂），进行搜索
        如果100次搜索没有结果，会尝试对最好解进行加深操作
        如果1000次搜索没有结果，会尝试扩大搜索范围

        */

        for(int i = 0 ;i<initSet.size();i++)
        {
            cout<<i<<endl;
            Solve s = initSet[i];
            s.value = TabuSearchJudge(s);
            if(bestValue<s.value)
            {
                best = s;
                bestValue = s.value;
            }
            int times = 0;
            //允许不相容情况的出现
            /*如果有500次没有更新，则退出*/

            while(times<200)
            {
                //暴力搜索
                //500次，从开放列表中进行遍历
                vector<Oper> oper;
                //通过多次随机操作产生大量可行解
                for(int j = 0;j<50;j++)
                {
                    oper.push_back(generateNeighbourOper(s,tabuList,bestValue));
                }
                sort(oper.begin(),oper.end(),[](const Oper &a, const Oper& b){
                    return a.weight<b.weight;
                });
                //找到禁忌表中允许的最小值，或者打破禁忌的特赦值
                for(int i = 0;i<oper.size();i++)
                {
                    int j1 = oper[i].j1,j2 = oper[i].j2;
                    int i1 = oper[i].i1,i2 = oper[i].i2;
                    bool pass = true;
                    //说明没有进一步选项
                    if(i1 == -1 && i2 == -1 && j1 == -1 && j2 == -1)
                        break;
                    //检查禁忌表，看是否存在项目
                    for(int j = tabuPoint;j<tabuEnd;j++)
                    {
                        if(j == tabuList.size())
                        {
                            j = 0;
                        }

                        if(tabuList[j].first == i1 && tabuList[j].second == j1 || tabuList[j].first==i2 && tabuList[j].second == j2)
                        {
                            pass = false;
                            break;
                        }
                    }
                    if(!pass && s.value + oper[i].weight > bestValue)//判断是否存在特赦
                    {
                        continue;
                    }
                    else//找到最小值
                    {
                        if(i2==-1)//单向移动
                        {
                            s.assignmentList[i1] = j2;
                            s.factoryFreeCapacity(j1,cliDemand[i1]);
                            s.factoryUseCapacity(j2,cliDemand[i1]);
                            tabuList[tabuEnd] = make_pair(i1,j1);
                            tabuEnd++;
                            if(tabuEnd == tabuList.size())
                                tabuEnd = 0;
                            if(tabuPoint == tabuEnd)
                            {
                                tabuPoint++;
                                if(tabuPoint == tabuList.size())
                                    tabuPoint = 0;
                            }
                        }
                        else
                        {
                            s.assignmentList[i1] = j2;
                            s.assignmentList[i2] = j1;
                            s.factoryFreeCapacity(j1,cliDemand[i1]);
                            s.factoryUseCapacity(j2,cliDemand[i1]);
                            s.factoryFreeCapacity(j2,cliDemand[i2]);
                            s.factoryUseCapacity(j1,cliDemand[i2]);

                            tabuList[tabuEnd] = make_pair(i1,j1);
                            tabuEnd++;
                            if(tabuEnd == tabuList.size())
                                tabuEnd = 0;
                            if(tabuPoint == tabuEnd)
                            {
                                tabuPoint++;