Capacitated Facility Location Problem容量有限设施选址问题

Capacitated Facility Location Problem

问题如下

给出所有工厂的容量和开工厂的成本，所有客户的需求，以及客户分配给某个工厂的的分配成本，要求解的问题是：找出一个分配方案，使得总成本最低。 实例数据下载地址：Download: p1-p71

算法思路

思路一：贪心

容易想到的是，从每个用户的角度出发，要使得总成本最小的话，则可以每个用户选择设施的时候都贪心选择一个成本最低的设施。一开始我想的是所有设施都先不开，然后由每个用户去贪心选择，计算选择一个设施需要的成本，如果设施还没开的话，则要加上开设施的成本，最后选择成本最低的那个。但后来仔细一想，其实这样是不太符合贪心策略的，因为用户肯定会倾向于选择已经开过的设施，这样用户就不能真正选择到分配成本更低的设施，所以不能充分利用到贪心的特点。所以不应该让设施的开与不开对选择造成太大的影响，一个简单的做法就是事先决定要开哪些设施，然后在开的设施里面贪心选择。我把所有工厂都开了，然后跑出来的效果确实比一开始的想法要好。 下面是python代码实现： 首先先建立好需要用到的类。

# 贪心算法
def greed(self, status):
    facilities = copy.deepcopy(self.facilities)
    customers = copy.deepcopy(self.customers)
    total_cost = 0
    self.assignment = [-1] * self.cnum
    for i in range(self.fnum):  # 先计算开工厂的费用
        if status[i] == 1:
            total_cost += facilities[i].opening_cost
    for i in range(self.cnum):  # 为每个用户分配一个工厂
        min_cost = 999999999
        for j in range(self.fnum):  # 从已开的工厂中选择有足够剩余量、且费用最小的
            if status[j] == 1 and facilities[j].rest >= customers[i].demand:
                cost = customers[i].assign_cost[j]
                if cost < min_cost:
                    min_cost = cost
                    self.assignment[i] = j
        if self.assignment[i] != -1:  # 加上分配费用
            total_cost += min_cost
            facilities[self.assignment[i]].rest -= customers[i].demand
        else:  # 若没有工厂能容纳下，则无效
            return -1
    self.status = status
    self.total_cost = total_cost
    return self.total_cost

复制

实验结果：（在读数据时把分派成本读错了，本来应该是从设备出发，设备对每个用户分配的成本，但读成了从用户出发，用户对应每个设备的成本，即读反了，所以导致算出来的最小值与参考值有一定差距，但由于数据结构的原因，要改起来比较麻烦，算法才是重点，所以没有改） result_table_greeddetail_greed

思路二：模拟退火法

模拟退火法是对物理退火过程的模拟，所谓退火是冶金专家为了达到某些特种晶体结构重复将金属加热或者或冷却的过程，该过程的控制参数为T。模拟退火法的基本思想是，在系统朝着能量减小的趋势这样一个变化过程中，偶尔允许系统跳到能量较高的状态，以避开局部极小点，最终稳定到全局最小点。 模拟退火法其实就是对爬山法（一种局部搜索算法）的一个改进。爬山法过程如下：

（1）生成第一个可能的解，若是目标则停止；否则转下一步。 （2）从当前可能的解出发，生成新的可能解集。

• 用测试函数测试新的可能解集中的元素，若是解，则停止；若不是解，转下一步。
• 将它与至今测试过的“解”相比较。若它最接近解，则保留作为最佳元素；若它不最接近解则舍弃。

（3）以当前最佳元素为起点，转第（2）步。 爬山法在生成的元素中寻找最优解，很容易陷入局部最优，而模拟退火的一个重要改进是以一定的概率接受差解，从而扩大了搜索空间，更可能进入全局最优。 模拟退火以一个温度作为参数，并当出现差解时，以一个与温度大小有关的概率接受差解，使得温度越低时，这个概率越小，最终达到稳定。并且在每个温度下，以一定的准则判断是否达到平衡状态，比如进行一定次数的内循环、连续200次没有出现更优解等。 对于设施选址问题，将用户的分配情况列为一个数组，作为解空间，则要在这个解空间内进行搜索。 代码如下：

# 模拟退火法
def SA(self):
    # 设置初温
    T = 1000
    for i in range(self.cnum):
        self.assignment[i] = -1
    current = self.assignment
    _next = None
    # good = 0
    # 降温
    while T > 0.01:
        no = 0  # 记录连续多少次没有接收优解了
        # print(T)
        T *= 0.95
        # 内循环直至在该温度下稳定
        for i in range(1000):
            # 生成新解
            _next = self.generate(current)
            while not self.is_valid(_next):
                _next = self.generate(current)
            dE = self.total(_next) - self.total(current)
            if dE < 0:  # 接收更优解
                no = 0
                # good += 1
                current = _next
                if self.total(self.assignment) > self.total(current):  # 记录到目前为止的最优解
                    self.assignment = current
            else:  # 以一定概率接受差解
                no += 1
                rd = random.random()
                if math.exp(-dE / T) >= rd:
                    current = _next
            if no > 200:
                break
    self.total_cost = self.total(self.assignment)
    for i in range(self.cnum):
        self.status[self.assignment[i]] = 1
    self.show()
    # print(good)

复制

生成邻接解的函数如下：

# 通过领域操作生成新解
def generate(self, current):
    customer_1 = random.randint(0, self.cnum - 1)
    customer_2 = random.randint(0, self.cnum - 1)
    _next = copy.deepcopy(current)
    while customer_2 == customer_1:
        customer_2 = random.randint(0, self.cnum - 1)
    method = random.randint(1, 4)
    if method == 1:  # 随机交换两个值
        _next[customer_1] = current[customer_2]
        _next[customer_2] = current[customer_1]
    elif method == 2:  # 随机改变两个值
        _next[customer_1] = random.randint(0, self.fnum - 1)
        _next[customer_2] = random.randint(0, self.fnum - 1)
    elif method == 3:  # 随机将某个值插到另一个位置
        if customer_1 < customer_2:
            for i in range(customer_1, customer_2):
                _next[i] = current[i + 1]
            _next[customer_2] = current[customer_1]
        else:
            a = range(customer_2 + 1, customer_1 + 1)
            for i in reversed(a):
                _next[i] = current[i - 1]
            _next[customer_2] = current[customer_1]
    elif method == 4:  # 倒置
        if customer_1 < customer_2:
            for i in range(customer_2 - customer_1):
                _next[customer_1 + i] = current[customer_2 - i];
        else:
            for i in range(customer_1 - customer_2):
                _next[customer_2 + i] = _next[customer_1 - i]
    return _next

复制

计算总成本的函数：

# 计算某分配的总成本
def total(self, assignment):
    status = [0] * self.fnum
    for i in range(self.cnum):
        status[assignment[i]] = 1
    cost = 0
    for i in range(self.cnum):  # 计算总分配费用
        cost += self.customers[i].assign_cost[assignment[i]]
    for i in range(self.fnum):  # 计算开工厂的总费用
        if status[i] == 1:
            cost += self.facilities[i].opening_cost
    return cost

复制

因为每次生成的邻接解可能是无效解，所以以上用到了is_valid函数，其定义如下：

# 判断一个分配序列是否有效
def is_valid(self, assignment):
    if not isinstance(assignment, list):
        return False
    for i in range(self.fnum):  # 对所有工厂，计算对某工厂总需求量
        total_demand = 0
        for j in range(self.cnum):
            if assignment[j] == i:
                total_demand += self.customers[j].demand
                if total_demand > self.facilities[i].capacity:  # 若对某工厂的总需求超过总容量，则无效
                    return False
    return True

复制

这种算法明显要比上面的贪心耗时很多，但基本上能比贪心找到更优的解。

实验结果： result_table_SAdetail_SA

思路三：贪心+模拟退火法

方法三其实就是前面两种方法的结合，但并不是说把贪心的解作为模拟退火的初始解。在用贪心算法的时候，是把设施的开闭状态固定下来的，是为了避免打开设施的成本影响用户的贪心选择（前面已经讲过），但固定的开闭状态不一定是最佳的，比如前面的贪心是把所有的设施都打开了，但可能关闭几个设施，再贪心，费用可能更小。所以第三个思路就是，求分配数组的时候采用贪心算法，把设施的开或不开作为解空间，然后用模拟退火法求解。

代码如下：

# 模拟退火法2
def SA2(self):
    min_cost = 99999999
    T = 10
    self.status = [1]*self.fnum
    current = self.status
    greed_current = self.greed(current)
    _next = None
    while T > 0.1:
        no = 0  # 记录连续多少次没有接收优解了
        T *= 0.95
        for i in range(12):
            # 生成新解
            _next = self.generate2(current)
            greed_next = self.greed(_next)
            while greed_next == -1:
                _next = self.generate2(current)
                greed_next = self.greed(_next)
            dE = greed_next - greed_current
            if dE < 0:  # 接收更优解
                no = 0
                current = _next
                greed_current = greed_next
                if min_cost > greed_current:  # 记录到目前为止的最优解
                    min_cost = greed_current
                    self.best_assignment = self.assignment
            else:
                no += 1
                rd = random.random()
                if math.exp(-dE / T) >= rd:
                    current = _next
                    greed_current = greed_next
                    if min_cost > greed_current:  # 记录到目前为止的最优解
                        min_cost = greed_current
                        self.best_assignment = self.assignment
            if no > 5:
                break
    self.total_cost = min_cost
    self.assignment = self.best_assignment

    for i in range(self.cnum):
        j = self.assignment[i]
        self.status[j] = 1
    self.show()

# 通过领域操作生成新解
def generate2(self, current):
    facility_1 = random.randint(0, self.fnum - 1)
    facility_2 = random.randint(0, self.fnum - 1)
    _next = copy.deepcopy(current)
    while facility_1 == facility_2:
        facility_2 = random.randint(0, self.fnum - 1)
    method = random.randint(1, 4)
    if method == 1:  # 随机交换两个值
        _next[facility_1] = current[facility_2]
        _next[facility_2] = current[facility_1]
    elif method == 2:  # 随机翻转两个值
        _next[facility_1] = 1 - current[facility_1]
        _next[facility_2] = 1 - current[facility_2]
    elif method == 3:  # 随机将某个值插到另一个位置
        if facility_1 < facility_2:
            for i in range(facility_1, facility_2):
                _next[i] = current[i + 1]
            _next[facility_2] = current[facility_1]
        else:
            a = range(facility_2 + 1, facility_1 + 1)
            for i in reversed(a):
                _next[i] = current[i - 1]
            _next[facility_2] = current[facility_1]
    elif method == 4:  # 倒置
        if facility_1 < facility_2:
            for i in range(facility_2 - facility_1):
                _next[facility_1 + i] = current[facility_2 - i];
        else:
            for i in range(facility_1 - facility_2):
                _next[facility_2 + i] = _next[facility_1 - i]
    return _next

复制

这个思路的解空间比思路二要小很多，所以循环次数比较少，效率更高一些，解的质量甚至比思路二还要好，但由于用了贪心，所以无法找到最优解，有些较小的解用思路二能够到达，而思路三不可能到达。 实验结果： result_table_SA2detail_SA2

实例数据下载地址：Download: p1-p71 所有源代码及实验结果github地址：CFLP