算法分析期末作业：Capacitated Facility Location Problem

求解算法

• 贪心算法求解
  贪心策略是不从整体最优上加以考虑，做出的是在某种意义上的局部最优解。这种方法能在比较短时间内得到问题的解，但只能对一些满足无后效性的问题求得全局最优解。而这个问题就具有后效性，而且暴力求解复杂度为n^m,n为顾客数量，m为工厂数量。
  本题的总费用由开厂费用和配送费用组成，基于这点贪心策略考虑的方向有两个，一个是每次开厂都优先用完容量满足尽可能多的顾客从而节省开厂费用，另一个是顾客都选择最近的厂从而节省配送费用。
  工厂费用优先算法步骤：
  (1).选择还没开启的工厂中开启费用最小的厂S加入开启的队列Open，如过没有可开启的工厂跳转(4)。
  (2).选取S附近未得到分配，且配送费用最小的顾客C，如果工厂内剩余容量足够则跳转(3),否则跳转(1)。
  (3).将C加入S的配送顾客队列Q（S)，如果全部顾客得到分配则结束，否则跳转(2)。
  (4).分配结束，计算总费用
  顾客距离优先算法步骤：
  (1).遍历顾客列表，选择一个顾客C执行(2)
  (2).寻找顾客C最近的工厂S，如果工厂未开启则开启工厂，如果工厂的容量容量足够则将顾客C加入S的配送顾客队列Q(S)并跳转(1)，否则重复(2)
  考虑了样例的特点是大部分样例的配送费用的比重是大于开厂的费用的，基于这点选择第二种效果更好，下面提供的也是第二种贪心策略的代码。（第一种贪心策略经过编码测试，只有很小部分样例结果优于第二种，所以不上传实现代码）

2. 模拟退火迭代法求解
   模拟退火算法是基于迭代求解策略的一种随机寻优算法。算法从某一较高初温出发，然后温度参数不断下调，系统能量朝着减少的趋势变化。在解空间中随机寻找目标函数的解，在局部最优解能概率性地跳出并最终趋于全局最优。
   算法步骤：
   (1).将顾客随机分配到各工厂当中，生成一个初始解s0，设定初始温度temperature，降温速率ratio,迭代次数t,迭代次数变化率vt。
   (2).进行t次如下操作，给当前解一个随机的位移（即进行局部搜索寻找一个新解，本例中随机选取一个顾客换到一个离他较近的工厂），计算总费用变化offset。 如果总费用的变化offset<0则直接接受该解为当前解。如果offset>0，则以一个随温度变化而变化的概率exp(-offset/temperature)接受该解，否则不接受该解。
   (3).更新温度temperature = temperature * ratio 和迭代次数t = t * vt,如果温度下降到预定值LIMIT，算法结束，否则继续(2)

C++源代码

//贪心算法
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <cfloat>
#include <ctime>
#include <string>

using namespace std;

//将结果信息输出
void debug(int k, float totaltime, 
    vector<bool> isOpen, vector<int> toFacility, vector<float> usedCap, int minCost) {
    ofstream ofile;
    ofile.open("result1.txt", ios::app);

    ofile << "p" << k << "," << minCost << "," << totaltime << ","; 
    for(int i = 0; i < isOpen.size(); i++){
        ofile << isOpen[i];
        if (i == isOpen.size()-1) {
            ofile << ",";
        }
        else{
            ofile << " ";
        }
    }
    for(int i = 0; i < toFacility.size(); i++){
        ofile << toFacility[i];
        if (i == toFacility.size()-1) {
            ofile << "\n";
        }
        else{
            ofile << " ";
        }
    }

    ofile.close();    
}

//读取数据文件
void readdata(string filename, vector<float> &capacity, vector<float> &fixedCost, vector<float> &demand, vector<vector<float> > &disCost) {
    ifstream ifile;
    ifile.open(filename, ios::in);

    int fNum, cNum;
    ifile >> fNum >> cNum;

    capacity.assign(fNum,0);
    fixedCost.assign(fNum,0);
    for(int i = 0; i < fNum; i++){
        ifile >> capacity[i] >> fixedCost[i];
    }
    
    demand.assign(cNum,0);
    for(int i = 0; i < cNum; i++){
        ifile >> demand[i];
    }

    disCost.assign(fNum,vector<float>(cNum,0));
    for(int i = 0; i < fNum; i++){
        for(int j = 0; j < cNum; j++){
            ifile >> disCost[i][j];
        }
    }

    ifile.close();
}

int main(int argc, char const *argv[]) {
    for(int k = 1; k <= 71; k++) {
        //程序开始计时
        clock_t start, end;
        start = clock();
        //读取数据文件
        string filename = "Instances/p" + to_string(k);
        vector<float> capacity, fixedCost, demand;
        vector<vector<float> > disCost;
        readdata(filename, capacity, fixedCost, demand, disCost);  
        //存储结果的数据结构
        int minCost = 0;
        vector<bool> isOpen(capacity.size(), false);
        vector<int> toFacility(demand.size(), -1);
        vector<float> usedCap(capacity.size(), 0);
        //贪心策略
        for(int i = 0; i < demand.size(); i++){
            vector<float> cost(capacity.size(), 0);

            for(int j = 0; j < capacity.size(); j++){
                cost[j] = disCost[j][i];
            }
            //每次取最近的能满足顾客要求的工厂
            vector<float>::iterator smallest;
            int index;
            while(true) {
                smallest = min_element(cost.begin(), cost.end());
                index = distance(cost.begin(), smallest);
                if(capacity[index] - usedCap[index] >= demand[i]) {
                    break;
                }
                else{
                    cost[index] = FLT_MAX;
                }
            }
            //记录分配结果
            if(!isOpen[index]) {
                minCost += fixedCost[index];
                isOpen[index] = true;
            }
            minCost += disCost[index][i];
            usedCap[index] += demand[i];
            toFacility[i] = index;
            
        }
        //结束程序计时
        end = clock();
        double totaltime = (double)(end - start) / CLOCKS_PER_SEC;
        debug(k,totaltime,isOpen, toFacility, usedCap, minCost);
    }
    
    return 0;
}

//模拟退火算法
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <cfloat>
#include <string>
#include <cmath>
#include <ctime>

using namespace std;

//模拟退火使用的参数
float TEMPERATURE = 100.0;  //初温
float T_RATIO = 0.98;       //降温速率
float LIMIT = 1.0;          //迭代停止温度条件
int ITERATION_TIME = 500;   //每个降温阶段的迭代次数
float ITERATION_RATIO = 1.01;//迭代次数随降温的提升比例

//将结果信息输出
void debug(int k, float totaltime, 
    vector<bool> isOpen, vector<int> toFacility, vector<float> usedCap, int minCost) {
    ofstream ofile;
    ofile.open("result2.txt", ios::app);

    ofile << "p" << k << "," << minCost << "," << totaltime << ","; 
    for(int i = 0; i < isOpen.size(); i++){
        ofile << isOpen[i];
        if (i == isOpen.size()-1) {
            ofile << ",";
        }
        else{
            ofile << " ";
        }
    }
    for(int i = 0; i < toFacility.size(); i++){
        ofile << toFacility[i];
        if (i == toFacility.size()-1) {
            ofile << "\n";
        }
        else{
            ofile << " ";
        }
    }

    ofile.close();    
}
//生成初始解
void init(vector<float> capacity, vector<float> fixedCost, vector<float> demand, vector<vector<float> > disCost, 
    vector<bool> &isOpen, vector<int> &toFacility, vector<float> &usedCap, int &minCost) {
    for(int i = 0; i < isOpen.size(); i++) {
        isOpen[i] = true;
        minCost += fixedCost[i];
    }

    for(int i = 0; i < demand.size(); i++) {
        int pos = rand() % capacity.size();
        for(int j = 0; j < capacity.size(); j++) {
            int index = (pos+j)%capacity.size();
            if(capacity[index] - usedCap[index] >= demand[index]) {
                usedCap[index] += demand[i];
                toFacility[i] = index;
                minCost += disCost[index][i];
                break;
            }
        }
    }
}
//模拟退火迭代函数
void iteration(vector<float> capacity, vector<float> fixedCost, vector<float> demand, vector<vector<float> > disCost, 
    vector<bool> &isOpen, vector<int> &toFacility, vector<float> &usedCap, int &minCost) {
    float temperature = TEMPERATURE;
    int iteration_time = ITERATION_TIME;
    while(temperature > LIMIT) {
        for(int i = 0; i < iteration_time; i++) {
            int pos_c = rand() % demand.size();
            float costbefore = disCost[toFacility[pos_c]][pos_c];
            if(usedCap[toFacility[pos_c]] == demand[pos_c]) {
                costbefore += fixedCost[toFacility[pos_c]];
            }
            
            int pos_f = rand() % capacity.size();
            float costafter = disCost[pos_f][pos_c];
            if(!isOpen[pos_f]) {
                costafter += fixedCost[pos_f];
            }

            if(pos_f != toFacility[pos_c] && capacity[pos_f] - usedCap[pos_f] >= demand[pos_c] && 
            (costafter < costbefore || 1.0*rand()/RAND_MAX < exp(-(costafter-costbefore) / temperature))) {
                minCost += costafter - costbefore;
                if(usedCap[toFacility[pos_c]] == demand[pos_c]) {
                    isOpen[toFacility[pos_c]] = false;
                }
                if(!isOpen[pos_f]) {
                    isOpen[pos_f] = true;
                }
                usedCap[toFacility[pos_c]] -= demand[pos_c];
                usedCap[pos_f] += demand[pos_c];
                toFacility[pos_c] = pos_f;
            }
        }
        iteration_time *= ITERATION_RATIO;
        temperature *= T_RATIO;
    }

}
 //读取数据文件
void readdata(string filename, vector<float> &capacity, vector<float> &fixedCost, vector<float> &demand, vector<vector<float> > &disCost) {
    ifstream ifile;
    ifile.open(filename, ios::in);

    int fNum, cNum;
    ifile >> fNum >> cNum;

    capacity.assign(fNum,0);
    fixedCost.assign(fNum,0);
    for(int i = 0; i < fNum; i++){
        ifile >> capacity[i] >> fixedCost[i];
    }
    
    demand.assign(cNum,0);
    for(int i = 0; i < cNum; i++){
        ifile >> demand[i];
    }

    disCost.assign(fNum,vector<float>(cNum,0));
    for(int i = 0; i < fNum; i++){
        for(int j = 0; j < cNum; j++){
            ifile >> disCost[i][j];
        }
    }

    ifile.close();
}

int main(int argc, char const *argv[]) {
    srand((unsigned int)time(NULL));
    for(int k = 1; k <= 71; k++) {
        //程序开始计时
        clock_t start, end;
        start = clock(); 
        //读取数据文件
        string filename = "Instances/p" + to_string(k);
        vector<float> capacity, fixedCost, demand;
        vector<vector<float> > disCost;
        readdata(filename, capacity, fixedCost, demand, disCost);
        //存储结果的数据结构
        int minCost = 0;
        vector<bool> isOpen(capacity.size(), false);
        vector<int> toFacility(demand.size(), -1);
        vector<float> usedCap(capacity.size(), 0);
        //生成初始解
        init(capacity, fixedCost, demand, disCost, isOpen, toFacility, usedCap, minCost);
        //进入模拟退火迭代
        iteration(capacity, fixedCost, demand, disCost, isOpen, toFacility, usedCap, minCost);
        //结束程序计时
        end = clock();
        double totaltime = (double)(end - start) / CLOCKS_PER_SEC;
        debug(k,totaltime,isOpen, toFacility, usedCap, minCost);
    }
    return 0;
}

结果表格

表格格式为

Problems(样例)	Results(结果)	Time(用时)	Open(工厂状态)	Facility(顾客分配情况)
p1	8848	0.030	0 1 0 1 0 0 1 1 0 1…	1 3 6 7 9…

贪心策略结果
https://blog.youkuaiyun.com/m0_37779608/article/details/85226459
模拟退火算法结果
https://blog.youkuaiyun.com/m0_37779608/article/details/85226480

可以看到两种算法用时都在几十毫秒，模拟退火算法用时平均较长，但求解结果相比贪心算法要好很多。