算法project——Capacitated Facility Location Problem

Capacitated Facility Location Problem是一个NP难的问题，为了在有限的时间内解决这个问题，我采用了两种启发式算法来解决它。第一种是模拟退火算法，第二种是禁忌搜索算法。两种算法应该说各有千秋。

下面是Capacitated Facility Location Problem的大致内容：

 Suppose there are n facilities and m customers. We wish to choose:

(1) which of the n facilities to open

(2) the assignment of customers to facilities

(3) The total demand assigned to a facility must not exceed its capacity

The objective is to minimize the sum of the opening cost and the assignment cost. 

简单来说就是每个工厂会有一个开启费用（没有顾客分配到这个工厂，则这个工厂就不会开启），然后把一个客户分配到一个工厂会有一个费用。问如何分配所有顾客使费用和最小。

除了上面所说的以外，每个工厂（设施）会有一个容量，然后每个客户会有一个容量的需求。然后分配到这个工厂的客户的容量需求之和不能超过这个工厂的容量。

数据集格式：

There are currently 71 data files. The format of these data files is:

 

|J| |I|

s_1 f_1

s_2 f_2

  ...

s_|J| f_|J|

d_1 d_2 d_3 … d_|I|

c_{11} c_{12} c_{13} … c_{1|I|}

c_{21} c_{22} c_{23} … c_{2|I|}

   ...    ...    ...   ....

c_{|J|1} c_{|J|2} c_{|J|3} … c_{|J||I|}

 

where:

 

|J| is the number of potential facility locations;

|I| is the number of customers;

s_j (j=1,...,|J|) is the capacity of facility j;

f_j (j=1,...,|J|) is the fixed cost of opening facility j;

d_i (i=1,...,|I|) is the demand of customer i;

c_{ji} (j=1,...,|J|), (i=1,...,|I|) is the cost of allocating all the demand of customer i to facility j.

 

首先我采用的是模拟退火算法：

首先随机分配初始解，并且判断其合法性。

我的模拟退火算法主要采用了两种邻域操作：

（1）将一个顾客随机分到另一种工厂去。

（2）将两个顾客所在的工厂交换。

当然，这两种操作必须保证容量。其实只要工厂容量允许，那么两次第一种操作就可以达到第二种操作的效果，但是我也想不到更多简单的邻域操作了。。。。。其他的邻域操作写起来都不容易。

我设计的初温t=1500，然后每次*0.97，外循环迭代3000次，这是为了让开始的动荡大一些。然后内循环次数将随着数据集的大小来改变，定为工厂数*顾客数/10。

为了防止陷入局部最优，我采用了重升温的操作，如果每30次外循环得到相同的总费用cost，那么将进行重升温操作，这里我定义了一个重升温参数ncnt, 重升温将把当前的温度乘ncnt，知道当前得到的cost改变。实际上，通过观察总费用cost的迭代过程，我观察到如果定义一个常数ncnt，那么每次的重升温仍有可能跳不出当前的局部最优而再次收敛。所以ncnt将随重升温次数的增多而变大。如果历史最优值成功更新（即得到了更好的结果），那么ncnt参数将重置为3。

因为需要重升温，所以必须保存历史最优值。

可以看到随着迭代次数增加，费用陷入局部最优9256，但是通过重升温，它找到了更优的局部最小值9135。 

 

坑点：在做第一种局部操作的时候，一定要注意如果选到的客户刚好在选到的工厂中，要如何操作。

下面给出代码：

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <sstream>
#include <cstdlib>
#include <cmath>
#include <ctime>


using namespace std;

struct factory{
    int cost;
    int capacity;
    int cnt;
    factory(int c1, int c2, int c3):cost(c1), capacity(c2), cnt(c3){
    }
};

struct customer
{
    int fac;
    int demand;
    vector<int> costs;
};

vector<factory> factorys;
vector<customer> customers;
vector<factory> bfactorys;
vector<customer> bcustomers;

int SA_search()
{
    int i, j, cnt = 0, cost = 0, pastcost = 0, bestcost = 1000000, past = 0;
    int lcost = 0;
    float ncnt = 0;
    float t;
    for(i=0;i<customers.size();i++)
    {
        int ran = rand()%factorys.size();
        if(factorys[ran].capacity - factorys[ran].cnt < customers[i].demand)
        {
            i--;
            continue;
        }
        factorys[ran].cnt += customers[i].demand;
        customers[i].fac = ran;
        cost += (customers[i].costs)[ran];
    }
    for(i=0;i<factorys.size();i++)
    {
        if(factorys[i].cnt!=0)
            cost += factorys[i].cost;
    }
    cout<<cost<<endl;
    t = 1500;
    cnt = 0;
    for(i=0;i<3000;i++)
    {
        for(j=0;j< customers.size() * factorys.size() / 10;j++)
        {
            int ran = rand()%10;
            if(ran<5)
            {
                int changedCus = rand()%customers.size();
                int changedFac = rand()%factorys.size();
                if(changedFac == customers[changedCus].fac)
                    continue;
                if(factorys[changedFac].capacity - factorys[changedFac].cnt < customers[changedCus].demand)
                    continue;
                customer cc = customers[changedCus];
                //past = cc.fac;
                int changedCost = - (cc.costs)[cc.fac] + (cc.costs)[changedFac];
                if(factorys[changedFac].cnt==0)
                    changedCost += factorys[changedFac].cost;
                if(factorys[cc.fac].cnt == cc.demand)
                    changedCost -= factorys[cc.fac].cost;
                int r1 = rand() % 32767;
                float r2 = r1 / 32767.0;

                if(exp(-changedCost/t)> r2)
                {
                    factorys[cc.fac].cnt -=cc.demand;
                    customers[changedCus].fac = changedFac;
                    factorys[changedFac].cnt +=cc.demand;
                    cost += changedCost;
                }
            }
            else
            {
                int changedCus1 = rand()% customers.size();
                int changedCus2 = rand()% customers.size();
                customer cc1 = customers[changedCus1];
                customer cc2 = customers[changedCus2];
                if(cc1.demand>cc2.demand)
                {
                    if(factorys[cc2.fac].capacity-factorys[cc2.fac].cnt<cc1.demand-cc2.demand)
                        continue;
                }
                else
                {
                    if(factorys[cc1.fac].capacity-factorys[cc1.fac].cnt<cc2.demand-cc1.demand)
                        continue;
                }
                int changedCost = (cc1.costs)[cc2.fac] + (cc2.costs)[cc1.fac];
                changedCost -= ((cc1.costs)[cc1.fac] + (cc2.costs)[cc2.fac]);
                int r1 = rand() % 32767;
                float r2 = r1  / 32767.0;
                if(exp(-changedCost/t)> r2)
                {
                    factorys[cc1.fac].cnt -= cc1.demand;
                    factorys[cc1.fac].cnt += cc2.demand;
                    factorys[cc2.fac].cnt -= cc2.demand;
                    factorys[cc2.fac].cnt += cc1.demand;
                    int tmp = cc1.fac;
                    customers[changedCus1].fac = cc2.fac;
                    customers[changedCus2].fac = tmp;
                    cost += changedCost;
                }
            }
        }

        t*=0.97;
        if(i%5==0)
        {
            if(cost==pastcost)
                cnt+=1;
            else
                cnt=0;
            if(cnt>=6)
            {
                t *= ncnt;
                ncnt++;
            }
            pastcost = cost;
            cout<<i<<" cost: "<<cost<<" t: "<<t<<" ncnt:"<<ncnt<<endl;
        }
        if(cost<bestcost)
        {
            bfactorys = factorys;
            bcustomers = customers;
            bestcost = cost;
            ncnt = 3;
        }
    }
    lcost = 0;
    for(i=0;i<bcustomers.size();i++)
    {
        lcost += (bcustomers[i].costs)[bcustomers[i].fac];
    }

    for(i=0;i<bfactorys.size();i++)
    {
        if(bfactorys[i].cnt!=0)
            lcost += bfactorys[i].cost;
    }
    cout<<lcost<<endl;
    return bestcost;
}

接下来是禁忌搜索：

禁忌搜索其实就是在爬山法上加入禁忌表，减少陷入局部最优的可能性。

我只采用了一种领域操作，就是上面模拟退火的第一种领域操作。通过1000次全局搜索，选择不在禁忌表中的最优值。禁忌表中存的是选择的客户和工厂。

禁忌表长度为10+工厂个数。

当搜索道德最优值小于历史最优值，则有可能出现破除禁忌的操作。即无视禁忌表取最优值。禁忌表以滑动数组的方式存放（因为stl中的queue不提供遍历操作）。

下面给出代码：

int Tabu_search()
{
    int tabulength = 10 + factorys.size();
    vector<int> tabuList;
    vector<int> tabuList2;
    int i, j, k, q, cnt = 0, cost = 0, bestcost = 1000000, mincost = 0,tfront = 0, minj = 0, mink = 0;
    for(i=0;i<customers.size();i++)
    {
        int ran = rand()%factorys.size();
        if(factorys[ran].capacity - factorys[ran].cnt < customers[i].demand)
        {
            i--;
            continue;
        }
        factorys[ran].cnt += customers[i].demand;
        customers[i].fac = ran;
        cost += (customers[i].costs)[ran];
    }
    for(i=0;i<factorys.size();i++)
    {
        if(factorys[i].cnt!=0)
            cost += factorys[i].cost;
    }
    cout<<cost<<endl;
    for(i=0;i<1000;i++)
    {
        mincost = 1000000;
        for(j=0;j<customers.size();j++)
        {
            for(k=0;k<factorys.size();k++)
            {
                int changedCus = j;
                int changedFac = k;
                if(changedFac == customers[changedCus].fac)
                    continue;
                if(factorys[changedFac].capacity - factorys[changedFac].cnt < customers[changedCus].demand)
                    continue;
                customer cc = customers[changedCus];
                int changedCost = - (cc.costs)[cc.fac] + (cc.costs)[changedFac];
                if(factorys[changedFac].cnt==0)
                    changedCost += factorys[changedFac].cost;
                if(factorys[cc.fac].cnt == cc.demand)
                    changedCost -= factorys[cc.fac].cost;
                if(cost+changedCost<bestcost)
                {
                    minj = j;
                    mink = k;
                    mincost = changedCost;
                }
                for(q=tfront;q<tabuList.size();q++)
                {
                    if(j!=tabuList[q])
                        continue;
                    if(k!=tabuList2[q])
                        continue;
                    break;
                }
                if(q!=tabuList.size())
                {
                    //cout<<i<<" "<<q<<" "<<tabuList.size();
                    continue;
                }
                if(changedCost < mincost)
                {
                    minj = j;
                    mink = k;
                    mincost = changedCost;
                }
            }
        }

        factorys[customers[minj].fac].cnt -=customers[minj].demand;
        customers[minj].fac = mink;
        factorys[mink].cnt +=customers[minj].demand;
        cost += mincost;
        tabuList.push_back(minj);
        tabuList2.push_back(mink);
        if(tabuList.size()-tfront>tabulength)
            tfront++;
        if(cost<bestcost)
        {
            bfactorys = factorys;
            bcustomers = customers;
            bestcost = cost;

        }
    }
    return bestcost;
}

 最后给出样例结果表格：

随着数据集增大，花费的时间增多。

	SA	time	Tabu	time
p1	8876	0.249	9038	0.461
p2	7920	0.235	8010	0.456
p3	9314	0.253	10010	0.42
p4	11121	0.242	11975	0.365
p5	8838	0.226	9503	0.385
p6	7855	0.233	8143	0.309
p7	9699	0.232	10099	0.412
p8	11436	0.236	12189	0.381
p9	8561	0.266	9040	0.537
p10	7704	0.259	7726	0.566
p11	9027	0.293	9726	0.575
p12	10826	0.285	10905	0.546
p13	8252	0.631	8762	1.326
p14	7219	0.644	7215	1.313
p15	9317	0.645	9635	1.328
p16	11482	0.655	12106	1.454
p17	8930	0.648	9129	1.433
p18	7125	0.66	7125	1.277
p19	9295	0.639	9790	1.354
p20	11573	0.65	11060	1.259
p21	8169	0.678	8919	1.491
p22	7389	0.669	7933	1.479
p23	9502	0.68	10129	1.514
p24	10756	0.671	11836	1.505
p25	13004	5.021	13259	12.4
p26	11306	3.601	13455	7.986
p27	13506	3.644	15037	7.832
p28	15706	3.632	18308	7.662
p29	13575	3.514	15723	7.518
p30	11889	3.562	12717	7.443
p31	14102	3.535	16413	7.442
p32	16395	3.536	21400	7.456
p33	12796	3.647	15126	7.846
p34	11285	3.617	11822	7.843
p35	13650	3.833	16420	9.08
p36	15004	4.011	16439	8.198
p37	12200	3.709	15126	8.94
p38	11149	4.006	12767	8.751
p39	12984	4.135	13931	8.478
p40	14984	3.75	19983	8.32
p41	7098	0.511	7112	0.939
p42	7080	1.138	8092	2.245
p43	6656	1.628	9123	3.935
p44	7133	0.504	7234	0.771
p45	7127	1.59	8251	3.959
p46	6864	1.626	8347	3.63
p47	6249	0.448	6555	0.727
p48	6189	1.52	8149	3.348
p49	7082	1.75	8058	3.653
p50	9239	0.515	9663	0.867
p51	8178	1.37	9686	2.884
p52	9369	0.507	10477	0.741
p53	9711	1.293	10848	2.655
p54	9273	0.501	9562	0.694
p55	8143	1.321	9540	2.833
p56	22746	4.66	23560	10.443
p57	29993	4.624	32882	10.414
p58	45117	4.696	53882	10.579
p59	35557	4.703	39121	10.239
p60	23532	4.84	23882	10.784
p61	30068	5.024	32260	11.114
p62	45490	4.854	52560	11.059
p63	32775	4.726	37423	10.955
p64	23114	4.857	23458	10.602
p65	30210	4.796	31799	10.851
p66	43632	4.815	52560	10.835
p67	34688	5.125	39671	10.509
p68	22810	5.108	23630	10.614
p69	30396	4.825	32260	10.701
p70	44716	4.963	52560	11.388
p71	34913	4.79	39232	10.586

下面列出前几个样例的输出:

p1

p2

p3

 

github地址：

https://github.com/smiletomisery/Capacitated-Facility-Location-Problem