遗传算法在离散型工厂选址问题中的应用

一、背景介绍
最近论坛上有一些同学私信要遗传算法的代码，这两天整理了一下，希望对大家有帮助。
1.1 工厂选址问题（转自去年写的退火算法）
工厂选址问题是运筹学中的经典问题之一，它描述的是在综合考虑工厂的建造成本、生产量，产品的运输成本，各地的需求量等因素后，如何制定适当的选址方案和物流运输方案来完成企业的生产经营活动。该问题的研究模型具有相当的普适应，其不仅仅在物流领域，还在生产生活，甚至军事中都有着非常广泛的应用。
以离散工厂选址问题为例，假设有n个工厂为m个配送中心送货，每个配送中心的需求量为bj，j=1,2,…，m。在满足配送中心需求的前提下，可以选取任意工厂组合来生产运输产品。假设如果选定某个工厂，这需要固定成本di（可以理解为建筑成本或者运营成本），i=1,2,…n,且每个工厂的生产能力为ai。如果从工厂i到配送中心j的单位运输成本为cij，那么在满足配送中心的需求量的情况下，如何制定选址方案（选哪些工厂）和运输方案，才能使得总费用最小？
1.2 遗传算法简介
遗传算法（Genetic Algorithm，GA）是一种群体优化算法，其主要思想来自达尔文的生物进化理论：通过生成初始种群（种群可以理解为解的结合）、自然选择、交叉、变异等模拟生物进化的概念来不断优化种群的解，因此它是一种以模拟自然进化过程为导向搜索最优解的方法。
二、问题设置
工厂选址问题的模型这里不赘述，本文的主要目的是如何编写遗传算法代码来解决该问题。问题参数设置如下：
假设有10个工厂（A-J）向20个需求点（2-21）供货，工厂到需求点之间的距离如下，保存在了“locationdata.csv”文件中：

需求地的需求量均设为6，每个工厂的供应量（capality）不等，数值如下：

demand = [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
capality = [30, 40, 50, 55, 45, 48, 36, 50, 60, 35]

假设选定某个工厂供货，这需要一个固定费用（可理解为建筑费用或运行费用），各工厂费用如下：

fix_cost_lst = [65, 50, 40, 30, 20, 50, 45, 30, 35, 25]

三、遗传算法求解
1.导入需要的包，设置参数：

import pandas as pd
import numpy as np
import pulp

data = pd.read_csv('locationdata.csv')
demand = [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
capality = [30, 40, 50, 55, 45, 48, 36, 50, 60, 35]
fix_cost_lst = [65, 50, 40, 30, 20, 50, 45, 30, 35, 25]
dic = {0:'A',1:'B',2:'C',3:'D',4:'E',5:'F',6:'G',7:'H',8:'I',9:'J'}

2.定义Population类，输入种群规模，随机生成初始种群：

class Population:
    def __init__(self,population_size):
        self.population_size = population_size

    def Initial_population(self):
        self.population = [np.random.randint(0,2,10) for i in range(self.population_size)]

3.迭代过程中选址方案的变化是一个部分随机的过程，这个过程中可能会产生需求量得不到满足的选址方案，这时就需要对不满足条件的方案进行调整：如果所有的需求量大于选中工厂的总供应量，那么随机添加一个未被选中的工厂，直到需求量能够得到满足。定义Fix_population类修正整体种群：

class Fix_population:
    def __init__(self,population):
        self.population = population
        self.fixed_population = []

    def Fix_population(self):
        for individual in self.population:
            while sum(demand) > self.Caculate_capality(individual):
                ind = np.random.randint(0,len(individual))
                if individual[ind] == 0:
                    individual[ind] = 1
            self.fixed_population.append(individual)

    def Caculate_capality(self,individual):
        total_capality = 0
        for i in range(len(individual)):
            if individual[i] == 1:
                total_capality += capality[i]
        return total_capality

3.定义费用计算类，其描述和上一篇退火算法一样：

class Caculate_cost:
    def __init__(self,individual):
        self.individual = individual

    def Transportation_problem(self,costs, x_max, y_max):
        row = len(costs)
        col = len(costs[0])
        prob = pulp.LpProblem('Transportation Problem', sense=pulp.LpMinimize)
        var = [[pulp.LpVariable(f'x{i}{j}', lowBound=0, cat=pulp.LpInteger) for j in range(col)] for i in range(row)]
        flatten = lambda x: [y for l in x for y in flatten(l)] if type(x) is list else [x]
        prob += pulp.lpDot(flatten(var), costs.flatten())

        for i in range(row):
            prob += (pulp.lpSum(var[i]) <= x_max[i])
        for j in range(col):
            prob += (pulp.lpSum([var[i][j] for i in range(row)]) == y_max[j])

        prob.solve()
        return {'objective':pulp.value(prob.objective), 'var': [[pulp.value(var[i][j]) for j in range(col)] for i in range(row)]}

    def Translist(self):
        transcost = [list(data[dic[i]]) for i in range(len(self.individual)) if self.individual[i] == 1]
        return transcost

    def Cal_cost(self):
        self.fix_cost = 0
        capality_new = []
        for i in range(len(self.individual)):
            if self.individual[i] == 1:
                self.fix_cost += fix_cost_lst[i]
                capality_new.append(capality[i])
        translist = self.Translist()
        costs = np.array(translist)

        max_plant = capality_new
        max_cultivation = demand
        res = self.Transportation_problem(costs, max_plant, max_cultivation)
        self.transpotation_schedule, self.total_cost = res['var'],res["objective"] + self.fix_cost

4.定义遗传算法类，在其中定义适应度计算、自然选择、交叉、变异、生成子代等方法模拟群体进化过程。其输入：父代种群，输出：自然进化后的子代种群。
这里需要对这些方法做简单说明：
（1）适应度直接采用目标函数值；
（2）自然选择的方法是轮盘赌；
（3）交叉方法采用最简单的单点交叉；
（4）变异方法为以一定概率改变个体的单个基因；
（5）生成子代时已经改变了那些不适应环境的个体，使其适应环境。

class GA:
    def __init__(self,population,cross_probility,mutation_problity):
        self.population = population
        self.cross_probility = cross_probility
        self.mutation_problity = mutation_problity
#计算种群每个个体的适应度，并累加，便于轮盘赌选择
    def Fitness(self):
        self.fitness = []
        self.transpotation_schedule = []
        for individual in self.population:
            demo = Caculate_cost(individual)
            demo.Cal_cost()
            self.fitness.append(demo.total_cost)
            self.transpotation_schedule.append(demo.transpotation_schedule)

        self.probility = [self.fitness[i] / sum(self.fitness) for i in range(len(self.fitness))]
        self.cumsum_probility = [sum(self.probility[0:i + 1]) for i in range(len(self.probility))]
        self.cumsum_probility.insert(0,0)
#寻找每一代个体中的最优解
    def Find_best_result(self):
        best_result = list(zip(self.population,self.transpotation_schedule,self.fitness))
        return [best_result[i] for i in range(len(best_result)) if best_result[i][2] == min(self.fitness)][0]
#轮盘赌选择
    def Select(self):
        self.sub_population = []
        for i in range(len(self.population)):
            probility = np.random.rand()
            for j in range(len(self.cumsum_probility)):
                if self.cumsum_probility[j] < probility <self.cumsum_probility[j+1]:
                    self.sub_population.append(self.population[j])
#单点交叉
    def Crossover(self):
        self.crossover_result = []
        sub_population = np.random.permutation(self.sub_population)
        for i in range(0,len(sub_population),2):
            if np.random.rand() <= self.cross_probility:
                cpoint = np.random.randint(0,len(sub_population[0]) + 1)
                temp1=[]
                temp2=[]
                temp1.extend(sub_population[i][0:cpoint])
                temp1.extend(sub_population[i + 1][cpoint:len(sub_population[i])])
                temp2.extend(sub_population[i + 1][0:cpoint])
                temp2.extend(sub_population[i][cpoint:len(sub_population[i])])
                self.crossover_result.append(temp1)
                self.crossover_result.append(temp2)
            else:
                self.crossover_result.append(sub_population[i])
                self.crossover_result.append(sub_population[i + 1])
#变异
    def Mutation(self):
        self.mutation_result = []
        for individual in self.crossover_result:
            if np.random.rand() <= self.mutation_problity:
                mutation_point = np.random.randint(0,len(individual))
                if individual[mutation_point] == 0:
                    individual[mutation_point] = 1
                else:
                    individual[mutation_point] = 0
            self.mutation_result.append(individual)
#生成子代群体
    def Next_population(self):
        fixed_population = Fix_population(self.mutation_result)
        fixed_population.Fix_population()
        self.next_population = fixed_population.fixed_population

5、主函数部分：
规定种群包含30个个体，进化10代，交叉概率0.5，变异概率0.2

init_population = Population(30)
init_population.Initial_population()
init_population = init_population.population
max_generation = 10

for i in range(max_generation):
    demo = Fix_population(init_population)
    demo.Fix_population()

    GA_demo = GA(demo.fixed_population,0.5,0.2)
    GA_demo.Fitness()
    print('='*80)
    print('第{}代中最优选址方案为:{}'.format(i + 1,GA_demo.Find_best_result()[0]),end = '\n')
    print('其最小费用为:{}'.format(GA_demo.Find_best_result()[2]),end = '\n')
    print('其最优运输方案为:{}'.format(GA_demo.Find_best_result()[1]))
    GA_demo.Select()
    GA_demo.Crossover()
    GA_demo.Mutation()
    GA_demo.Next_population()
    init_population = GA_demo.next_population

四、运行结果
截取部分结果：

第6代中最优选址方案为:[0 1 1 0 0 0 0 1 1 1]
其最小费用为:422.0
其最优运输方案为:[[0.0, 0.0, 6.0, 6.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 6.0, 0.0, 6.0, 6.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 6.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 6.0, 0.0, 2.0, 0.0, 0.0, 0.0], 
[0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 6.0, 0.0],
[6.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 6.0]]
================================================================================
第7代中最优选址方案为:[0 1 0 0 0 0 1 1 1 1]
其最小费用为:409.0
其最优运输方案为:[[0.0, 0.0, 6.0, 6.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 6.0, 6.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0], 
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 6.0, 0.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 6.0, 0.0], 
[6.0, 0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 6.0]]

如上所示，该问题最优个体出现在第7代，最小费用为409。