本篇博客涉及贪心算法以及模拟退火算法解决单源设施选址问题

问题描述

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
The objective is to minimize the sum of the opening cost and the assignment cost.
The total demand assigned to a facility must not exceed its capacity.

数据集格式

此处AssignmentCost 有No.Customer行，No.Facility列。即每个用户对于不同的设备有不同的Assignment cost

思路分析

贪心算法

贪心算法是比较简单但不能够获得最优解的方法。使用贪心算法可以考虑较少的情况。我在贪心算法的编写中，忽略了对于用户分配顺序的考虑，直接从前到后进行分配。对于每个用户，首先通过用户的需求DEMAND来筛选出可用的设备，接着通过比较这些设备的opencost以及用户在这些设备上的assignment cost，选出其中的最小值。
显然，贪心算法得到的解并不是最好，并且可能导致后续用户没有可用分配

模拟退火算法

在本题中，我选择将所有用户的设备分配情况作为一个design，而总的opencost + assignment cost则作为目标方程。在生成新的Design时，一般情况下是选择当前解的领域解，但需要考虑这样的情况

因此，一直从后往前选择的话可能会变成下山法，所以我采用直接从前往后进行选择。即产生的解为随机解

代码实现

1. 初始化

import random
import math
import matplotlib.pyplot as plt

CAPACITY = []
OPENCOST = []
ASSIGNCOST = []
OPENSTATUS = []
DEMAND = []
FACNUM = 0
CUSNUM = 0
Num = 0

#读取文件，获取各项变量
def init(filepath):
    with open(filepath, 'r') as f:
        global FACNUM
        global CUSNUM
        global CAPACITY
        global OPENCOST
        global ASSIGNCOST
        global OPENSTATUS
        FACNUM, CUSNUM = f.readline().strip("\n").split()
        FACNUM = int(FACNUM)
        CUSNUM = int(CUSNUM)
        OPENSTATUS = [0] * FACNUM
        for i in range(FACNUM):
            line = f.readline().strip("\n").split()
            CAPACITY.append(int(line[0]))
            OPENCOST.append(int(line[1]))

        for i in range(int(CUSNUM / 10)):
            line = f.readline().strip("\n").replace(' ', '').split(".")
            for j in range(10):
                DEMAND.append(int(line[j]))

        for i in range(CUSNUM):
            linetoread = int(FACNUM / 10)
            temp = []
            for j in range(linetoread):
                line = f.readline().strip("\n").replace(' ', '').split(".")
                for k in range(10):
                    temp.append(int(line[k]))

            ASSIGNCOST.append(temp)

    print(ASSIGNCOST)

2. 贪心算法

def greedy():
    global CAPACITY

    capacity = CAPACITY.copy()
    openstatus = [0] * FACNUM
    totalcost = 0
    assign = [-1] * CUSNUM
    assigncost = 0
    opencost = 0

    for cus in range(CUSNUM):
    
        facIndex = []       # the facilities may be chosen
        for i in range(FACNUM):
            if capacity[i] >= DEMAND[cus]:
                facIndex.append(i)

        # ASSIGNCOST for customer i
        asforeach = ASSIGNCOST[cus]

        # ASSIGNCOST + OPENCOST for customer i
        temp = [sum(x) for x in zip(asforeach, OPENCOST)]

        # select the min cost index
        sindex = facIndex[0]
        for i in facIndex:
            if temp[i] < temp[sindex]:
                sindex = i

        openstatus[sindex] = 1
        capacity[sindex] = capacity[sindex] - DEMAND[cus]
        assign[cus] = sindex

        assigncost += asforeach[sindex]

    for i in range(FACNUM):
        opencost += openstatus[i] * OPENCOST[i]

    totalcost = assigncost + opencost

    print(assign)
    print("assignment cost: %d"%assigncost)
    print("opencost: %d"%opencost)
    print(totalcost)
    print("openstatus: {}".format(openstatus))

    return assign

3. 模拟退火

class Anneal:
    assign = []  # list of Customer's choice , len = NumofCus
    capacity = []
    openstatus = []
    totalcost = 0
    assigncost = 0
    opencost = 0

    def __init__(self):
        self.openstatus = [0] * FACNUM

        self.assign = [-1] * CUSNUM
        # 注意！！list的赋值是传引用，不是真正的copy。应该换为list.copy
        self.capacity = CAPACITY.copy()

        for i in range(CUSNUM):
            available = []
            for j in range(FACNUM):
                if self.capacity[j] >= DEMAND[i]:
                    available.append(j)


            choose = random.choice(available)
            self.capacity[choose] = self.capacity[choose] - DEMAND[i]
            self.assign[i] = choose

            self.assigncost += ASSIGNCOST[i][self.assign[i]]
            self.openstatus[self.assign[i]] = 1

        for j in range(FACNUM):
            self.opencost += self.openstatus[j] * OPENCOST[j]

        self.totalcost = self.assigncost + self.opencost



    def f(self):
        return self.totalcost


def simulated_Annealing():
    # Number of cycles
    n = 50
    # Number of trials per cycle
    m = 50
    # Number of accepted solutions
    na = 0.0
    # Probability of accepting worse solution at the start
    p1 = 0.7
    # Probability of accepting worse solution at the end
    p50 = 0.001
    # Initial temperature
    t1 = -1.0 / math.log(p1)
    # Final temperature
    t50 = -1.0 / math.log(p50)
    # Fractional reduction every cycle
    frac = (t50 / t1) ** (1.0 / (n - 1.0))

    # Initialize x
    x = []
    newassign = Anneal()
    x.append(newassign)

    xi = x[0]
    na = na + 1.0

    # Current best results so far
    xc = x[0]
    fc = xc.f()
    
    fs = []
    fs.append(fc)

    # Current temperature
    t = t1

    # delta Average
    delta_avg = 0.0

    for i in range(n):
        print('Cycle: ' + str(i) + ' with Temperature: ' + str(t))
        for j in range(m):
            # Generate new trial points
            xi = Anneal()
            delta = abs(xi.f() - fc)

            if xi.f() > fc:
                # Initialize delta_avg if a worse solution was found
                #   on the first iteration
                if i == 0 and j == 0:
                    delta_avg = delta

                # objective function is worse
                # generate probability of acceptance
                p = math.exp(-delta / (delta_avg * t))

                # determine whether to accept worse point
                if random.random() < p:
                    # accept the worse solution
                    accept = True
                else:
                    # don't accept the worse solution
                    accept = False
            else:
                accept = True

            if accept is True:
                xc = xi
                fc = xi.f()
                # increment number of accepted solutions
                na = na + 1.0
                # update delta_avg
                delta_avg = (delta_avg * (na - 1.0) + delta) / na

        # Record the best x values at the end of every cycle
        x.append(xc)
        fs.append(fc)

        # Lower the temperature for next cycle
        t = frac * t

    # print solution
    print('Best solution: ' + str(xc.assign))
    print('Best objective: ' + str(fc))
    print('assign cost: %d' % xc.assigncost)
    print('opencost: %d' % xc.opencost)
    print('openstatus:{}'.format(xc.openstatus))

    fig = plt.figure()
    ax1 = fig.add_subplot(211)
    ax1.plot(fs, 'r.-')
    ax1.legend(['Objective'])
    ax2 = fig.add_subplot(212)
    ax2.plot(xc.assign, 'b.')
    ax2.legend(['assignment'])

    # Save the figure as a PNG
    plt.savefig('iterations.png')

    plt.show()

运行结果

1. 贪心
   

2. 模拟退火
   
   

3. 所有测例结果

结果分析

从结果看到，模拟退火过程最后得到的结果并不是很好，原因可能在于模拟退火选择的初始态以及在寻找新解的过程都是随机的。可能的改进为：

1. 使用贪心算法的结果作为初始值
2. 产生新解的过程可以通过：
   • 从中间往前改变
   • 仅仅进行点交换
3. 对于产生的较差的结果，可以通过爬山法的结果进行替换