使用遗传算法解决城镇设施建造问题

问题建模/配方

约束管理

图1 原始数据

        首先看一下原始数据(图1)，这是一个58 * 58的矩阵。数据显示了每两个城镇的驾驶时间。问题是我们需要在一些城镇建设设施，而所有城镇都可以被设施覆盖。

        其次，有5种类型的设施:3分钟，5分钟，10分钟，15分钟和20分钟。因此，我们需要对这些需求的数据进行预处理。然后根据不同的要求，将实值矩阵转换为二进制矩阵。

        然后，假设要求是20，我们可以看到图2:

• 值1表示这两个城镇可以通过建筑设施到达对方。
• 值为0表示这两个城镇不能同时到达，我也将主对角线的值设置为0。它们是没有意义的，而且会干扰我的算法操作

图2 使用要求为20的二值化数据

        第三，个体的基因由城镇的索引构成，每个数字代表一个城镇。

        第四，根据上面的描述，我们可以个体的是一个数组，数组内的数字由0-57构成，代表这58个城镇。因此，问题可以描述为样本中数字的交叉和突变等。

Fitness函数设计

        fitness函数可以换个词来说可以是cost函数或者loss函数。这个函数就是用来评价设计的算法。在这个问题中，评价的标准是建设设施的数量。fitness值越小越好。

        我的fitness值分为两个部分，一个是连接的城镇的数量，一个是孤立城镇的数量。这很容易理解，有一些城镇距离其他城镇都很远，超出了要求的大小，不与其他城镇连接，这些就是孤立城镇。

孤立的城镇

图3 孤立城镇的计算

        孤独小镇不连接其他小镇，比意味着这个小镇的列和行都被填满0。就意义而言，这个城镇是任何其他城镇都无法到达的。假设2号镇是孤独的，我们可以看到第2列和第2行都是0。

        在此基础上，我们可以得到解。如果有一行表示该数字全部为0。这个小镇一定很孤独。这个值是不变的。但不同的要求有不同数量的孤独小镇

        这个方法的返回值是孤立小镇的数量，以及孤立小镇的编号。

连接的城镇

图4 计算可连接矩阵需要建造设施的数量

         前面提到过，种群中个体的染色体的构成是有0-57的数组，但是数组内的数字顺序是随机的，染色体表示建造顺序，即从染色体的第一个数字开始建造，一直到全部覆盖为止。因此这个函数方法就是用来计算连接小镇需要建造设施的数量。

        这个方法是一个递归方法，染色体也是进行处理过的。染色体需要把孤立小镇的编号进行移除。下面是我设计计算的简单介绍：

• 输入一个全为1的模板
• 根据染色体的小镇编号选择建造的小镇
• 根据这个小镇生成新的模板，并与之前的模板进行融合
• 判断这个模板是不是符合只有孤立的小镇为1，可连接的小镇为0
• 不是的话循环2-4，是的话直接输出

Example

这里用一个4*4的矩阵进行展示，同样是适用于58*58的情况。

小镇二值化矩阵 = [[0, 1, 0, 1]    样本1= [1, 2, 0, 3]  样本2= [0, 1, 3, 2]  
                              [1, 0, 0, 0]     
                              [0, 0, 0, 0] 
                              [1, 0, 0, 0]] 

孤立小镇： 2号小镇， 孤立小镇的数量：1。

将这个孤立小镇的编号从染色体中移除：

样本1= [1, 0, 3]  样本2= [0, 1, 3]

可连接小镇需要的设施计算：

• 样本1：

1. 模板= [1, 1, 1, 1], 设施数量 = 0
2. 样本1[0] = 1，二值化矩阵点乘样本的转置不等于0向量
3. 新的模板=[0, 0, 1, 1]，将新的模板之前的模板进行融合成为模板=[0,0,1,1]，设施数量 = 1
4. 样本1[1] = 0，二值化矩阵点乘样本的转置不等于0向量
5. 新的模板=[0, 0, 1, 0]，将新的模板之前的模板进行融合成为模板=[0,0,1,0]，设施数量 = 2
6. 样本1[1] = 0，二值化矩阵点乘样本的转置等于0向量
7. 结束

对于样本1来说：需要建设的设施数量为1+2=3

• 样本2：

1. 模板= [1, 1, 1, 1], 设施数量 = 0
2. 样本2[0] = 0，二值化矩阵点乘样本的转置不等于0向量
3. 新的模板=[0, 0, 1, 0]，将新的模板之前的模板进行融合成为模板=[0,0,1,0]，设施数量 = 1
4. 样本1[1] = 3，二值化矩阵点乘样本的转置等于0向量
5. 结束

对于样本2来说：需要建造的设施数量为1+1=2

应用

种群初始化

        所有样本由0-57的字符串组成，代表1 - 58个城镇。

        然后，我们需要从样本中去除孤独的城镇。我们需要随机打乱字符串数据。并进行了大量的样本生成。

图5 种群初始化

        如图5所示，如果需求是20。首先，我们将得到4个孤独的城镇，将这些小镇从染色体中去除。其次，我们建立了一个20 * 54的矩阵。第三，将所有样本的数据设置为相同。最后，打乱所有样本中的数据。

        通过以上源代码，生成形状为20*54。染色体中的数据也是随机打乱的。

新种群的生成

杂交

        我们得到第一代。生成形状为20*(58-孤立小镇的数量)。接下来，总体中的每两个样本需要随机重组。两个父母会生两个孩子，并把这些孩子放入种群中。因此，生成规模为40 *(58 -孤立小镇数量)，如图6所示。

图6 杂交​​​​

        杂交其实就是两个父母，将两个排列组合成两个新的排列。

1. 随机在染色体上选择一个杂交断点。
2. 将父母的一部分给第一个孩子，另外一部分给第二个孩子。需要对出生的孩子进行染色体合理性检查。如果出现冲突的数据，进行处理。图7为子代的染色体组合，图8为染色体的合理性检查处理。染色体不合理的情况是个体染色体中出现重复的数字，即重复的在一个小镇上建造。这意味着在另外一个孩子会缺失这个城镇。因此需要找到他们并进行交换。

图7 子代染色体生成

图8 染色体合理性检查处理

 幸存者选择

        在进行一轮杂交之后，种群中的个体的数量翻倍。因此需要对种群中的个体进行选择，只保留一半数量作为新的族群。

图9 新种群的选择

        选择的依据就是fitness值。通过对每一个个体计算fitness值，然后进行排序。选择出最小的20个作为新的种群，种群变异比例由超参数GENERATION_MUTATION决定。具体函数处理步骤如下：

1. 建立一个空的种群
2. 计算现有种群中40个个体的fitness值
3. 通过fitness值从小到大将这些个体进行排序
4. 选择前20个个体放入新的种群。同时记录最好和最差的fitness值。
5. 如果最好和最差的fitness值大小相等，那就进行变异。

变异

        突变是随机打乱样本中的数据。置乱比由超参数MUTION决定。源代码如图10所示。突变步骤如下：

1. 对每个位置数据都进行判断
2. 如果需要交换此数据，将此数据与染色体中的随机位置数据交换

图10 变异

参数的设计以及实验结果

参数设计

固定参数

• 父母数量：20
• 后代数量：20
• 幸存者选择机制：将父母和孩子都放在种群中，通过fitness值排列进行筛选。前20个个体将成为幸存者
• 初始化：随机初始化20个个体作为第一代。
• 停止条件：种群繁殖1000次

超参数

• GENERA TION_MUTATION:用于确定种群中突变的比例。取值范围为0.1 ~ 0.9。他们代表人口的10% - 90%的样本

1. 示例:如果GENERA TION_MUTA TION = 0.3，则表示种群的最后30%需要突变

• MUTATION:用于确定样本中数据变异的比例。取值范围为1 ~ 9，代表样本数据的10% ~ 90%。

1. 例:MUTATION = 0.3，表示样本中染色体有30%的突变可能性。

参数调优

        对于两个超参数的排列和组合。我们可以得到81组数据。在此基础上，为了减少随机性的影响，每组超参数分别进行了5个实验。我画了不同要求下fitness值的变化曲线。因此我总共得到了405条变化曲线(图11)。在超参数初始化中，我将两个参数都设置为0。

图11 实验结果数据

 参数控制

        这两个超级参数需要调整。为了获得算法在不同超参数下求解问题的效果。

        这里，每调整一次超参数，超参数都会增加10%(图12)

图12 参数控制

 实验结果

        现在有5个要求，分别是3分钟、5分钟、10分钟、15分钟、20分钟。因此，这五个需求总共生成了405张图片。

        每个变化曲线如图13所示，显示了第一个运行时10% GENERA TION_MUTA TION和10% MUTATION的结果

图13 实验结果数据案例

        通过对405组数据的处理。我们可以得到以下结果，这是五次运行的平均分。纵坐标为GENERA TION_MUTA TION，横坐标为MUTA TION。

要求: 3分钟

	10%	20%	30%	40%	50%	60%	70%	80%	90%
10%	55	55	55	55	55	55	55	55	55
20%	55	55	55	55	55	55	55	55	55
30%	55	55	55	55	55	55	55	55	55
40%	55	55	55	55	55	55	55	55	55
50%	55	55	55	55	55	55	55	55	55
60%	55	55	55	55	55	55	55	55	55
70%	55	55	55	55	55	55	55	55	55
80%	55	55	55	55	55	55	55	55	55
90%	55	55	55	55	55	55	55	55	55

要求：5分钟

	10%	20%	30%	40%	50%	60%	70%	80%	90%
10%	45	45	45	45	45	45	45	45	45
20%	45	45	45	45	45	45	45	45.2	45
30%	45	45	45	45	45	45	45	45	45
40%	45	45	45	45	45	45	45	45	45
50%	45	45	45	45	45	45	45	45	45
60%	45	45	45	45	45	45	45	45	45
70%	45	45	45	45	45	45	45	45	45
80%	45	45	45	45	45	45	45	45	45
90%	45	45	45	45	45	45	45	45	45

要求：10分钟

	10%	20%	30%	40%	50%	60%	70%	80%	90%
10%	27.8	27.8	27.6	27.8	28	28	27.8	28.2	28
20%	27.6	27.2	27.4	27.6	28.2	27.6	28	28	27.8
30%	27	27.2	27.8	27.2	27.4	27.8	27.4	27.6	27.6
40%	27.8	27.2	27.2	27.4	27	27.8	27.2	27.6	27.8
50%	27.2	27.8	27.4	27.2	27.4	27.6	27.6	27.6	27.2
60%	27	27	27.4	27.6	27.8	27.2	27.8	28.2	27.6
70%	27	27.4	27.2	27.4	27.6	27.4	27.2	27.4	27.4
80%	27.6	27.6	27.4	27.6	27.8	27	27.4	27.4	27.4
90%	27.6	27.4	27.4	27.2	27.2	27.2	27.2	27.4	27.6

要求：15分钟

	10%	20%	30%	40%	50%	60%	70%	80%	90%
10%	21.4	20.8	21	21	21	21.2	21.4	22	21
20%	20.8	20.6	21	20.8	21.2	21.2	20.8	20.8	21
30%	20.4	21	21.2	21.4	21	21	21.2	21.4	21
40%	21	21	20.6	21	20.8	21.2	20.8	20.8	21
50%	20.8	20.8	21.2	20.8	20.8	21.4	20.8	20.8	21
60%	20.8	20.8	21.4	20.8	21.4	21.6	21.2	21.2	21.2
70%	20.8	20.8	20.8	21.2	20.6	21	20.8	20.8	21
80%	20.6	21	20.6	20.4	21	21	20.2	20.6	20.4
90%	20.8	20.4	20.6	20.4	20.8	20.8	21	20.8	21

要求：20分钟

	10%	20%	30%	40%	50%	60%	70%	80%	90%
10%	14.4	14.8	14.2	14	14.2	14.8	15.4	14.2	15
20%	14.8	13.6	13.4	14.4	14.8	14.8	14.2	14.4	15
30%	13.2	14	14.4	13.6	14.8	14.4	14.6	15	14.2
40%	13.4	13.4	13.8	13.4	15	14.4	13.4	14.6	14.2
50%	13.6	13.8	14	13.8	14	14.4	14.4	14.6	14.6
60%	13.6	13.6	14.2	14	14.2	14.8	14.4	14.2	14.6
70%	13.8	13	13.6	13.8	14.2	13.6	14.4	14.2	14.2
80%	13.4	13.6	13.4	14.2	14.4	14.2	14.4	13.8	14.2
90%	13.4	14	13.4	13.4	14	13.8	14	13.8	14.2

完整代码

# 开发作者  : Tian.Z.L
# 开发时间  : 2022/12/2 16:26
# 文件名称  : IntelligentComputing1_3.py
# 开发工具  : PyCharm
# @Function:
import openpyxl
import pandas as pd
import numpy as np
import random
import os

try:
    import cv2.cv2 as cv2
except:
    import cv2
import matplotlib.pyplot as plt

TENSOR_CALCULATE_SUM_COL = np.ones([58, 1])
TENSOR_CALCULATE_SUM_ROW = np.ones([1, 58])
NUMBER_GENERATION = 20
Facility_choice = [3, 5, 10, 15, 20]
FACILITY_NEED = 0
TIMES = 5
random.seed(42)
LONELY_TOWN = 0


def calculate_min(template, array, list):
    '''
    calculate the lower number of facilities needed
    :param template:
    :param array: Binary matrix
    :param list: Sample
    :return: end flag
    '''
    global FACILITY_NEED
    tensor_rows = np.dot(array, template)
    if np.count_nonzero(tensor_rows) == 0 or FACILITY_NEED == 58-LONELY_TOWN:
        return 0
    template_this_iter = np.where(array[int(list[FACILITY_NEED]), :] == 1, 0, 1).reshape(-1, 1)
    template_this_iter[int(list[FACILITY_NEED])] = 0
    template = np.minimum(template, template_this_iter)
    FACILITY_NEED += 1
    calculate_min(template, array, list)


def lonely_town(array):
    """
    calculate the number of lonely towns
    :param array: Binary matrix
    :return: number of lonely towns ,connectable town array
    """
    tensor_rows = np.dot(array, TENSOR_CALCULATE_SUM_COL)
    tensor_clos = np.dot(TENSOR_CALCULATE_SUM_ROW, array)
    index_rows = np.nonzero(tensor_rows)
    index_cols = np.nonzero(tensor_clos)
    lonely_town_array = np.unique(np.hstack((index_rows[0], index_cols[1])))
    lonely_town_num = 58 - len(lonely_town_array)
    connectable_town_array = index_rows[0]
    return lonely_town_num, connectable_town_array


def random_disruption(array_r):
    """
    Random disruption the string
    :param array_r: string input
    :return: string shuffled
    """
    np.random.shuffle(array_r)
    return array_r


def offspring_rationality_check(offspring1, offspring2):
    offspring1_index_array = []
    offspring2_index_array = []
    for i in np.unique(offspring1):
        idx = np.argwhere(offspring1 == i)
        if len(idx.reshape(1, -1).squeeze(0)) == 2:
            offspring1_index_array.append(idx.reshape(1, -1).squeeze(0))
    for i in np.unique(offspring2):
        idx = np.argwhere(offspring2 == i)
        if len(idx.reshape(1, -1).squeeze(0)) == 2:
            offspring2_index_array.append(idx.reshape(1, -1).squeeze(0))
    for i in range(len(offspring1_index_array)):
        x = random.randint(0, 3)
        if x == 0:
            offspring1_cache1 = offspring1[offspring1_index_array[i][0]]
            offspring1[offspring1_index_array[i][0]] = offspring2[offspring2_index_array[i][0]]
            offspring2[offspring2_index_array[i][0]] = offspring1_cache1
        elif x == 1:
            offspring1_cache2 = offspring1[offspring1_index_array[i][0]]
            offspring1[offspring1_index_array[i][0]] = offspring2[offspring2_index_array[i][1]]
            offspring2[offspring2_index_array[i][1]] = offspring1_cache2
        elif x == 2:
            offspring1_cache2 = offspring1[offspring1_index_array[i][1]]
            offspring1[offspring1_index_array[i][1]] = offspring2[offspring2_index_array[i][0]]
            offspring2[offspring2_index_array[i][0]] = offspring1_cache2
        else:
            offspring1_cache2 = offspring1[offspring1_index_array[i][1]]
            offspring1[offspring1_index_array[i][1]] = offspring2[offspring2_index_array[i][1]]
            offspring2[offspring2_index_array[i][1]] = offspring1_cache2
    return offspring1, offspring2


def mutation(array):
    '''
    Mutate the sample
    :param array:sample
    :return:sample mutated
    '''
    global MUTATION
    for i in range(len(array)):
        random_number = random.randint(1, 10)
        if random_number <= MUTATION:
            choise = random.randint(0, 57-LONELY_TOWN)
            cahche = array[i]
            array[i] = array[choise]
            array[choise] = cahche
    return array


def crossover(array1, array2):
    """
    generate two children through crossover
    :param array1: parent1
    :param array2: parent2
    :return: child1, child2
    """
    proportion = random.randint(2, 10)
    offspring_1 = np.hstack((array1[:array1.shape[0] // proportion], array2[array2.shape[0] // proportion:]))
    offspring_2 = np.hstack((array2[:array2.shape[0] // proportion], array1[array1.shape[0] // proportion:]))
    offspring_1, offspring_2 = offspring_rationality_check(offspring_1, offspring_2)
    return offspring_1, offspring_2


def random_generation(ori_array, number_generation):
    """
    Random initialization of the first generation
    :param ori_array: input origin string
    :param number_generation: number of samples in generation
    :return:generation
    """
    array_a = np.zeros([number_generation, 58-LONELY_TOWN])
    for i in range(number_generation):
        array_a[i, :] = ori_array
    for i in range(number_generation):
        np.random.shuffle(array_a[i, :])
        array_a[i, :] = random_disruption(array_a[i, :])
    return array_a


def selection(array_group, binary_matrix):
    """
    select survivors as new generation
    :param array_group:generation.size:40* (58 - number of lonely towns)
    :param binary_matrix: binary matrix.size:58*(58 - number of lonely towns)
    :return:new generation.size:20*(58 - number of lonely towns)
    """
    global FACILITY_NEED, GENERATION_MUTATION, LONELY_TOWN
    loss_list = []
    array_survival = np.zeros([array_group.shape[0] // 2, array_group.shape[1]])
    for i in range(array_group.shape[0]):
        calculate_min(TENSOR_CALCULATE_SUM_COL, binary_matrix, array_group[i, :])
        loss_list.append(FACILITY_NEED + LONELY_TOWN)
        FACILITY_NEED = 0
    index_from_small2big = np.argsort(np.array(loss_list))
    index_survival = index_from_small2big[:20]
    for i in range(len(index_survival)):
        array_survival[i, :] = array_group[index_survival[i], :]
    best_loss = loss_list[index_survival[0]]
    worst_loss = loss_list[index_survival[19]]
    if best_loss == worst_loss:
        for i in range(array_survival.shape[0] - int(array_survival.shape[0] * GENERATION_MUTATION),
                       array_survival.shape[0]):
            array_survival[i, :] = mutation(array_survival[i, :])
    return best_loss, worst_loss, array_survival


def new_born(generation, number):
    """
    Parents cross at random to produce the children generation
    :param generation: parents generation
    :param number: number of samples in the generation
    :return: generation include children generation and parents generation
    """
    arr = np.arange(number)
    np.random.shuffle(arr)
    array_group = np.zeros([number * 2, 58-LONELY_TOWN])
    array_group[:number, :] = generation
    generation = generation[arr]
    for i in range(number // 2):
        offspring_1, offspring_2 = crossover(generation[i * 2, :], generation[i * 2 + 1, :])
        array_group[number + 2 * i, :] = offspring_1
        array_group[number + 2 * i + 1, :] = offspring_2
    return array_group


if __name__ == '__main__':
    df = pd.read_csv(
        r"C:\Users\User\Downloads\Assignment 01 Facility-Location Set-Covering Problem-20221115\drivingtime.csv",
        header=None)
    df.head(58)
    number_generation = 20
    train_epoch = 1000
    time_array = np.array(df)
    time_array = np.where(time_array == 2., 30, time_array)
    GENERATION_MUTATION = 0
    MUTATION = 0
    for gen in range(9):
        GEN_num = (gen+1) * 10
        GENERATION_MUTATION += 0.1
        MUTATION = 0
        for mu in range(9):
            MUTATION += 1
            MU_num = (mu+1) * 10
            os.mkdir(r'C:\Users\User\Desktop\assignment\ResultImage2\G{}_M{}'.format(GEN_num, MU_num))
            print("Mutation scale: {}    Population Mutation scale: {}".format(MUTATION / 10, GENERATION_MUTATION))
            sort_array = np.arange(58)
            for time in range(TIMES):
                print("\ntime: {}".format(time))
                print("---------------------------------")
                for index in range(len(Facility_choice)):
                    # for index in range(4, 5):
                    best_loss_array = []
                    worst_loss_array = []
                    binary_matrix = np.where(time_array <= Facility_choice[index], 1, 0)
                    need, connectable_town_array = lonely_town(binary_matrix)
                    LONELY_TOWN = need
                    # sort_array = np.setdiff1d(sort_array, lonely_town)
                    sort_array_1 = random_generation(connectable_town_array, number_generation)
                    for i in range(train_epoch):
                        # print("epoch:{}/{}".format(i + 1, train_epoch))
                        new_born_group = new_born(sort_array_1, number_generation)
                        best_loss, worst_loss, select_group = selection(new_born_group, binary_matrix)
                        sort_array_1 = select_group
                        best_loss_array.append(best_loss)
                        worst_loss_array.append(worst_loss)
                    worst_loss_array = np.array(worst_loss_array)
                    best_loss_array = np.array(best_loss_array)
                    plt.plot(worst_loss_array, label='worst data')
                    plt.plot(best_loss_array, label='best data')
                    plt.xlabel("generation")
                    plt.ylabel("loss")
                    plt.title('Requirement_{}  Best:{}'.format(Facility_choice[index], np.min(best_loss_array)))
                    plt.legend()
                    plt.savefig(
                        r'C:\Users\User\Desktop\assignment\Facility_{}.jpg'.format(Facility_choice[index]))
                    # plt.show()
                    plt.close()
                    # array = sort_array_1[0, :]
                    print("Facility_" + str(Facility_choice[index]) + ": ", np.min(best_loss_array))
                Facility_3 = cv2.imread('Facility_3.jpg')
                Facility_5 = cv2.imread('Facility_5.jpg')
                Facility_10 = cv2.imread('Facility_10.jpg')
                Facility_15 = cv2.imread('Facility_15.jpg')
                Facility_20 = cv2.imread('Facility_20.jpg')
                result = np.ones([Facility_3.shape[0], Facility_3.shape[1], 3])
                # result = cv2.resize(result, (Facility_3.shape[1], Facility_3.shape[0]))
                resylt = np.vstack(
                    (np.hstack((Facility_3, Facility_5, Facility_10)), np.hstack((Facility_15, Facility_20, result))))
                cv2.imwrite(
                    r'C:\Users\User\Desktop\assignment\ResultImage2\G{}_M{}\generation_{}_mutation_{}_times_{}.bmp'.format(GEN_num, MU_num,
                        GEN_num, MU_num, time), resylt)