Capacitated Facility Location Problem

Capacitated Facility Location Problem

标签（空格分隔）： 课堂项目

姓名：李**
学号：16340114
题目：Capacitated Facility Location Problem

---

题目概要

有n个工厂与m个顾客。每个顾客有自己的需求(demands)，每个工厂也有自己能满足的需求上限(capacity)。开启一个工厂需要一定花费，将一个顾客分配给一个工厂也需要一定花费，每个工厂开启的花费可能不同，一个顾客分配给不同工厂的花费也可能不同，不同顾客分配给同一工厂的花费也可能不同。问如何在不超过工厂capacity的情况下，满足所有顾客的需求，并使总花费最小。

思路

使用贪心算法或者动态规划是可以解决这一题的，但是问题的解空间很大，使用这种精确求解方法开销巨大，而且策略也十分难构思出来。近几节课老师讲了几种近似求解的方法，我觉得很适合把这些方法用在这个项目中。在这个项目中，我使用了爬山法与模拟退火两种方法。
　　在爬山法与模拟退火中，需要明确解要以怎样的形式表示，如何搜索其邻域，如何计算一个解的花费。下面，逐一讲解。

解的形式

记工厂数量为n，顾客数量为m。那么解的形式为一个长度为m的序列：$[a_0,a_1,...,a_{m-1}]$，其中$a_i$的取值范围为$[0,n)，其含义为顾客i分配给$a_i$工厂。其中被分配了顾客的工厂为开启状态，没有被分配顾客的工厂为关闭状态。（这个解的形式是和小伙伴们讨论出来的~，应该是一个比较简洁的表示形式了）

合法性检验

考虑到工厂有容量上限，上述解的形式不一定是正确，因此需要验证解的合法性。
　　给定一个解，计算其中每个工厂接受的顾客数量，如果存在工厂的顾客接受数量超过了该工厂的容量，即这一个解不合法。反之，这一个解合法。

领域操作

搜索neighbour的策略，在解的领域中通过以下随机一个方法产生neighbour：
　　1. 随机挑选解中的两个位置，交换值。即随机挑选两个顾客，把他们被分配的工厂调换。
　　2. 随机挑选解中的三个位置，交换值。即随机挑选三个顾客，把他们被分配的工厂调换。
　　3. 随机挑选一个顾客，重新赋值。即随机挑选一个顾客，重新将他分配给新的工厂。

解的花费计算

给定一个解，先计算每个顾客分配给对应工厂的花费（一开始读数据还读错了，要好好看readme呀~），再判断各个工厂是否开启，计算开启的工厂所需的花费，上述数值的总和即为一个解的花费。

爬山法

确定了解的形式与领域操作后。实现爬山算法就非常容易了。
　　爬山法的核心思想就是随机产生一个解，在解的领域内进行搜索，在邻居中接受更好的解（花费更少的解），再以新的解为中心搜索其领域，直到解的领域中没有更好的解为止。
　　在这个项目中，先随机生成一个解，注意，这个解可能是不合法的，要验证其合法性。具体做法就是重复生成随机解，直到解是合法的为止。
　　然后进入迭代过程：对当前解进行邻域操作，发现其邻居，当然，邻居也可能是不合法的。同样，重复进行邻域操作，直到邻居合法为止。将合法的邻居与当前的解的花费作比较，取花费更少的为接受解。然后进入下一轮迭代。
　　在本项目中，终止条件为执行了了一定的次数的迭代。经过测试，执行了m*n次操作后，解趋于稳定。使用动态迭代次数，可以使小规模的测例更快得出结果，也能使大规模的测例有足够次数的搜索去得到爬山法的“山顶”。

结果

写程序输出到文件里再粘贴过来，一开始还手填，蠢死了

instances	result	time
p1	12230	0.02
p2	10707	0.02
p3	11833	0.02
p4	14157	0.02
p5	11953	0.02
p6	10658	0.02
p7	12931	0.02
p8	14699	0.02
p9	11500	0.02
p10	11142	0.02
p11	12524	0.02
p12	14829	0.02
p13	13068	0.04
p14	11360	0.04
p15	13725	0.04
p16	15956	0.04
p17	12659	0.04
p18	12093	0.04
p19	13970	0.04
p20	17962	0.04
p21	12679	0.04
p22	10287	0.04
p23	14297	0.04
p24	18004	0.04
p25	20303	0.30
p26	17215	0.32
p27	24295	0.30
p28	30215	0.30
p29	21672	0.30
p30	19276	0.30
p31	23704	0.31
p32	29554	0.31
p33	21437	0.30
p34	19346	0.32
p35	20317	0.31
p36	26582	0.31
p37	21445	0.30
p38	17939	0.30
p39	20845	0.30
p40	26350	0.31
p41	9021	0.05
p42	10910	0.07
p43	9864	0.09
p44	10361	0.13
p45	12034	0.07
p46	10946	0.10
p47	8937	0.08
p48	9725	0.07
p49	10665	0.09
p50	11738	0.08
p51	12940	0.10
p52	13785	0.33
p53	15423	0.11
p54	11771	0.59
p55	13256	0.11
p56	31433	0.49
p57	40322	0.48
p58	62294	0.49
p59	47995	0.48
p60	31643	0.49
p61	41138	0.50
p62	62816	0.49
p63	48179	0.50
p64	31449	0.50
p65	40746	0.50
p66	60972	0.47
p67	48093	0.49
p68	31412	0.48
p69	38786	0.50
p70	60302	0.48
p71	47916	0.49

详细结果

https://blog.csdn.net/Ray0758/article/details/85227460

模拟退火

模拟退火与爬山法的思想很接近，不同之处为模拟退火有一个温度变量，在温度较高时有较高的概率接受稍差的解。
　　与爬山法相似，先生成一个合法的初始解。
　　设置初始温度$T=100$
　　然后进入迭代过程：对当前解进行邻域操作，计算概率
　　$$p=e^{\frac{-diff\times k}{T}},diff=cost(s_i^{{}^{\prime }})-cost(s_i)$$
　　这意味着有p的概率接受邻居。注意到，当邻居比当前解花费更少时，p为大于1，即100%接受新解。当邻居比当前解花费更多时，温度高的情况下概率p还是比较大的，而温度低时该概率p趋向0。而为了使邻居比当前解差时的p落在一个合理的区间内，添加常数k进行修正。经过几次测试，diff的范围大部分落在在400_{500区间，取k=0.1进行修正。这样在10}90时p的概率能落在比较大的范围内，使模拟退火取得更好的结果。
　　上述迭代过程进行$m\ast n/10$次（测试得数据，理由同爬山法）。同一温度下，$m\ast n/10$次迭代完成后。通过$T=0.95\times T$修改温度，直到温度降至0.1以下，结束模拟退火。

结果

instances	result	time
p1	9356	0.21
p2	8038	0.23
p3	9828	0.22
p4	11892	0.23
p5	9323	0.24
p6	8036	0.24
p7	10297	0.23
p8	11903	0.23
p9	8828	0.20
p10	7860	0.20
p11	10235	0.22
p12	11214	0.27
p13	9293	0.51
p14	7907	0.45
p15	11822	0.46
p16	11442	0.43
p17	9541	0.43
p18	7646	0.41
p19	10890	0.42
p20	12762	0.41
p21	9264	0.42
p22	8562	0.46
p23	11140	0.48
p24	13067	0.53
p25	14624	3.98
p26	12534	3.85
p27	14567	3.88
p28	19242	3.84
p29	14497	4.13
p30	12545	3.85
p31	15039	3.84
p32	18035	3.94
p33	14418	4.01
p34	12493	3.68
p35	15045	3.96
p36	17805	4.05
p37	14744	3.89
p38	12480	4.02
p39	14547	3.86
p40	17251	3.87
p41	7224	0.49
p42	7605	0.84
p43	7038	1.27
p44	7336	0.57
p45	7716	0.88
p46	7933	1.15
p47	6740	0.54
p48	7878	0.88
p49	6802	1.30
p50	9447	0.67
p51	9136	1.25
p52	9616	0.80
p53	10645	1.48
p54	9513	0.67
p55	9665	1.23
p56	23961	6.30
p57	32723	6.29
p58	48905	6.34
p59	37336	6.41
p60	24153	6.14
p61	31307	6.11
p62	49810	6.08
p63	39088	6.13
p64	23666	6.09
p65	31966	6.15
p66	48052	6.28
p67	38597	6.16
p68	24192	6.12
p69	31634	6.12
p70	48803	6.23
p71	37993	6.62

详细结果

https://blog.csdn.net/Ray0758/article/details/85227481

两种方法对比

计较两个算法得出的结果，从直观上看爬山法的计算速度特别快，但是解的质量却难以令人满意，原因为爬山法很容易陷在局部最优解上，难以找到全局最优。而模拟退火的速度比爬山法稍慢，规模小的测例运算得算是很快了，因为模拟退火接受了稍差的解，有机会跳出局部最优，从而找到更好的解。将两个算法的解与网上能找到最优解比较，爬山法的解约为最优解的1.2_{1.3倍，而模拟退火的解大多数为最优解的1.0}1.1倍。

具体实现

主要使用了python的list数据结构进行解的存储与邻域操作，其他的按照思路写就好，没什么特别的。实现的时候使用了不少global变量，虽然不太好，应该塞进一个类里的，但是简单起见嘛，哈哈一些变量名和函数名都写得清楚的地方就没有多少注释，简洁至上。

心得

本来以为这个项目会很难的，但动手做起来还是挺轻松的，没有想像中那么难。在人工智能课上也做过遗传算法的实现，对做这个项目也有很大帮助，有不少地方是共同的。做这个项目收获还是很大的，毕竟写出了一个还算不错的SA算法~~，虽然和大佬做的比起来还差了不少~~。

源码：

运行环境：python3
模拟退火:

import queue
import random
import numpy
import time

# global variable
facilitiesNum = 0
customerNum = 0
facilitiesCost = []
facilitiesCapacity = []
customersDemands = []
customersCost = []

# randomly assigns customers to facilities
def initial():
    global facilitiesNum
    global customerNum
    plan = [0 for i in range(customerNum)]
    for i in range(0, customerNum):
        plan[i] = random.randint(0, facilitiesNum-1)

    return plan

# check whether customers' demands exceed capacity
def isValid(plan):
    global facilitiesNum
    global customerNum
    global facilitiesCapacity
    global customersDemands

    allocation = [0 for i in range(facilitiesNum)]
    for i in range(0, customerNum):
        facility = plan[i]
        allocation[facility] += customersDemands[i]

    for i in range(0, facilitiesNum):
        if allocation[i] > facilitiesCapacity[i]:
            return False

    return True


def evaluate(plan):
    global facilitiesNum
    global customerNum
    global facilitiesCost
    global customersCost

    isOpen = [False for i in range(facilitiesNum)]
    totalCost = 0

    for i in range(0, customerNum):
        facility = plan[i]
        isOpen[facility] = True
        totalCost += customersCost[facility][i]

    for j in range(0, facilitiesNum):
        if isOpen[j]:
            totalCost += facilitiesCost[j]

    return totalCost



def selectNeighbour(plan):
    global facilitiesNum
    global customerNum

    p = random.random()
    neighbour = plan[:]

    # select two customers
    cut1 = random.randint(0, customerNum-1)
    cut2 = random.randint(0, customerNum-1)

    if (cut1 > cut2):
        cut1, cut2 = cut2, cut1


    if p < 0.3:
        # swap their facilities
        temp = neighbour[cut1]
        neighbour[cut1] = neighbour[cut2]
        neighbour[cut2] = temp
    elif p < 0.6:
        # assign one of them to new facility
        neighbour[cut1] = random.randint(0, facilitiesNum-1)
    else:
        # select one more, and swap theirs
        cut3 = random.randint(0, customerNum-1)
        temp = neighbour[cut1]
        neighbour[cut1] = neighbour[cut2]
        neighbour[cut2] = neighbour[cut3]
        neighbour[cut3] = temp

    return neighbour

# function to solve a instance of file
def solve(filename):
    global facilitiesNum
    global customerNum
    global facilitiesCost
    global customersCost
    global facilitiesCapacity
    global customersDemands

    # reset the global variable
    facilitiesNum = 0
    customerNum = 0
    facilitiesCost = []
    facilitiesCapacity = []
    customersDemands = []
    customersCost = []

    # read file
    file = open(filename, "r")
    data = queue.Queue()

    # split all the word
    for line in file:
        nums = line.split()
        for num in nums:
            data.put(num)

    # get important number: m and n
    facilitiesNum = int(float(data.get()))
    customerNum = int(float(data.get()))

    # get facilities' infomation
    for i in range(0, facilitiesNum):
        facilitiesCapacity.append(int(float(data.get())))
        facilitiesCost.append(int(float(data.get())))

    # get customers' demands
    for i in range(0, customerNum):
        customersDemands.append(int(float(data.get())))

    # get customers' assignments cost
    customersCost = [[0 for i in range(customerNum)] for j in range(facilitiesNum)]
    for i in range(0, facilitiesNum):
        for j in range(0, customerNum):
            customersCost[i][j] = int(float(data.get()))

    # SA
    # initial a solution utill it is valid
    plan = initial()
    while not isValid(plan):
        plan = initial()

    iteration = facilitiesNum * customerNum // 10
    currentCost = evaluate(plan)

    T = 100

    while T > 0.1:
        # inner iteration
        for i in range(0, iteration):
            # find neighbour utill valid
            neighbour = selectNeighbour(plan)
            while not isValid(neighbour):
                neighbour = selectNeighbour(plan)

            neighbourCost = evaluate(neighbour)
            diff = neighbourCost - currentCost

            p = numpy.exp(-diff*0.1/(T))

            if random.random() < p:
                plan = neighbour
                currentCost = neighbourCost
        # change T
        T = 0.95 * T

    return currentCost, plan


# main
lastTime = time.time()
for i in range(1, 72):
    # get cost and solution
    cost, plan = solve("./instances/p" + str(i))
    print("p%d: %d" % (i, cost))

    # calculate runtime
    now = time.time()
    print("Time: ", now - lastTime)
    lastTime = now

    # get facilities state
    isOpen = [0 for i in range(facilitiesNum)]
    for i in range(0, customerNum):
        facility = plan[i]
        isOpen[facility] = 1

    print("Facilities open or not:")
    print(isOpen)
    print("The assignment of customers to facilities:")
    print(plan)
    print()

爬山法：（不同的部分，其余部分相同）

iteration = facilitiesNum * customerNum
currentCost = evaluate(plan)

for i in range(0, iteration):
    neighbour = selectNeighbour(plan)
    while not isValid(neighbour):
        neighbour = selectNeighbour(plan)

    neighbourCost = evaluate(neighbour)
    if neighbourCost < currentCost:
        plan = neighbour
        currentCost = neighbourCost

输出成老师要求的格式，只要改改输出，改成输出到文件里就好。代码就不贴了。