遗传算法实践(八) 简单配送问题

作者: ChaoesLuol | 来源:发表于2019-09-16 09:12 被阅读0次

遗传算法实践(八) 简单配送问题
遗传算法综述及简单应用实例
遗传算法实践(九) VRP问题
遗传算法实践(三) 背包问题
遗传算法helloworld级别的python实现（结果可视化）
模拟退火算法(附代码)
遗传算法实践(二) 配词问题
遗传算法实践(十) VRPTW问题求解
遗传算法简单介绍与MATLAB实现（二）
旅行商问题的近似最优解（局部搜索、模拟退火、遗传算法）

简单配送问题描述

基本配送计划模型是从产地 $i=1,2,...,I$ 往销地 $j=1,2,...,J$ 配送产品时，在满足产地的供应量和销地的需求量这一约束前提下，制定配送总成本最低的配送计划。即
$min \ z = \sum_{i=1}^{I}\sum_{j=1}^{J}c_{ij}x_{ij} \\ s.t. \sum_{j=1}^J x_{ij} \le a_i,\ i=1,2,3,...,I \\ \sum_{i=1}^Ix_{ij}\ge b_j,\ j=1,2,...,J \\ x_{ij}\ge0,\ \forall i,j \\ \sum_{i=1}^I a_i=\sum_{j=1}^J b_j$

简单配送问题示例

考虑有4个产地和5个销地的配送计划问题，供应量和需求量如下:
$a_1=8, a_2=4, a_3=12, a_4=6 \\ b_1=3, b_2=5, b_3=10, b_4=7, b_5=5$

产地和销地之间每单位数量商品的运输成本为：

	销地1	销地2	销地3	销地4	销地5
产地1	8.6860922	9.35319846	4.04839742	6.64822502	2.21060537
产地2	2.09958268	7.93166898	5.69745574	4.18627465	6.11192203
产地3	7.05011006	9.65028246	5.79551805	2.74632325	5.96366128
产地4	6.66228395	8.73934611	4.89579209	6.13292158	5.65538837

遗传算法求解简单配送问题

个体编码

这里采用基于优先级的编码方式。

设有 $I$ 个产地， $J$ 个销地， $b_i$ 为产地 $i$ 的供应， $d_j$ 为销地 $j$ 的需求， $c_{ij}$ 为从产地 $i$ 运输单位数量商品到销地 $j$ 的代价， $y_{ij}$ 为从产地 $i$ 运往销地 $j$ 的商品数量。

染色体编码 $v$ 长度为 $I+J$ ，前 $I$ 个位置代表产地，后 $J$ 个位置代表销地，其编码过程如下：

Step 0: 初始化，优先级初始化为 $p\leftarrow (I+J)$ ， $v_t\leftarrow 0,\forall t\in(I+J)$ ；
Step 1: 从供应量和需求量中有一个不为0的产-销地组合中，找到运输代价最小的，并为其分配最大运输量(该组合对应的供应量和需求量中较小者)， $(i^\prime,j^\prime)\leftarrow argmin\{(c_{ij}|y_{ij}\ne 0), (b_{i^\prime}=y_{ij}||d_{j^\prime}=y_{ij})\}$ ；
Step 2: 更新供应量和需求量， $b_{i^\prime}=b_{i^\prime}-y_{i^\prime j^\prime}$ ， $d_{j^\prime}=d_{j^\prime}-y_{i^\prime j^\prime}$ ；
Step 3: 根据更新后的供应量和需求量分配优先级

$if\ b_{i^\prime}=0, then\ v_{i^\prime}\leftarrow p, p\leftarrow p-1; \\if\ d_{j^\prime}=0, then\ v_{I+j^\prime}\leftarrow p, p\leftarrow p-1;$

Step 4: 重复Step 1-3直到 $(b_{i^\prime}=0,\ \forall i^\prime \in I)\&(d_{j^\prime}=0,\ \forall j^\prime \in J)$
Step 5: 为染色体中仍然为0的剩余位置随机填充优先值

解码

解码过程如下：

Step 0：初始化各产销地配对的运输量 $y_{ij}\leftarrow 0, \forall i\in I, j\in J$ ；
Step 1：选择优先度最高的节点 $l\leftarrow argmax\{v\}$ ；
Step 2：如果该节点是销地，则从产地选择运输代价最小的节点 $if\ l\in J, then\ j\leftarrow l$ ， $i*\leftarrow argmin\{ c_{il}|v_i\ne 0, i\in I\}$ ；如果该节点是产地，则从销地选择运输代价最小的节点 $Else, \ i\leftarrow l$ ， $j*\leftarrow argmin\{ c_{lj}|v_j\ne 0, j\in J\}$ ；
Step 3：为选择的产-销地更新可用的最大运输量 $y_{i*j*}\leftarrow min\{b_{i*},d_{j*}\}$ ；
Step 4：更新剩余的需求量和销量 $b_i\leftarrow b_i-y_{i*j*}$ ， $d_j\leftarrow d_j - y_{i*j*}$ ；
Step 5：如果完成了需求量和销量的分配，将优先度归零。 $if\ b_{i*}=0,then\ v_{i*}\leftarrow 0$ ， $if\ d_{j*}=0,then\ v_{I+j*}\leftarrow 0$
重复Step 1-5直到优先度编码全部归零。

交叉操作

交叉操作采用PMX交叉

突变操作

使用shuffleIndex突变

代码示例

完整代码如下：

## 环境设置
import numpy as np
import matplotlib.pyplot as plt
from deap import base, tools, creator, algorithms
import random

params = {
    'font.family': 'serif',
    'figure.dpi': 300,
    'savefig.dpi': 300,
    'font.size': 12,
    'legend.fontsize': 'small'
}
plt.rcParams.update(params)

from copy import deepcopy
## 问题描述
creator.create('FitnessMin', base.Fitness, weights=(-1.0,))
creator.create('Individual', list, fitness=creator.FitnessMin)

## 个体编码
sourceFlow = [8, 4, 12, 6]
destinationFlow = [3, 5, 10, 7, 5]

def matEncode(sourceFlow=sourceFlow, destinationFlow=destinationFlow):
    '''生成可行的配送量矩阵'''
    rowNum = len(sourceFlow)
    colNum = len(destinationFlow)
    posSet = list(range(rowNum * colNum))  # 当前未填充的位置集合
    mat = np.zeros((rowNum, colNum))  # 初始化染色体矩阵
    sourceFlowCopy = deepcopy(sourceFlow)
    destinationFlowCopy = deepcopy(destinationFlow)
    while posSet:
        randIdx = random.randint(0, len(posSet)-1)  # 随机选取位置集中的一个元素
        posNum = posSet[randIdx]
        rowPos = (posNum - 1)//colNum + 1  # 第rowPos行
        colPos = posNum % colNum  # 第colPos列
        # 填充矩阵相应位置
        mat[rowPos-1, colPos -
            1] = min(sourceFlowCopy[rowPos-1], destinationFlowCopy[colPos-1])
        # 根据填充值，修改剩余的产量和销量
        sourceFlowCopy[rowPos-1] -= mat[rowPos-1, colPos-1]
        destinationFlowCopy[colPos-1] -= mat[rowPos-1, colPos-1]
        # 从位置集合中删去已经被填充的位置
        posSet = posSet[:randIdx] + posSet[randIdx+1:]
    return mat.tolist()

# Cost Matrix
costMat = [
    [8.6860922, 9.35319846, 4.04839742, 6.64822502, 2.21060537],
    [2.09958268, 7.93166898, 5.69745574, 4.18627465, 6.11192203],
    [7.05011006, 9.65028246, 5.79551805, 2.74632325, 5.96366128],
    [6.66228395, 8.73934611, 4.89579209, 6.13292158, 5.65538837]
]

## 基于优先级的编码
def priorityCoding(matCode, costMat=costMat, sourceFlow=sourceFlow, destinationFlow=destinationFlow):
    '''
    从给定的配送量矩阵和代价矩阵,生成基于优先级的编码
    输入: matCode -- 一个满足sourceFlow和destinationFlow的可行解
    costMat -- 代价矩阵,给出由一个source到一个destination的代价
    sourceFlow, destinationFlow -- 每个source的输出能力和destination的接受能力
    输出: 长度为len(sourceFlow) + len(destinationFlow)的编码
    '''
    priorityCode = [0] * (len(sourceFlow) + len(destinationFlow)) # 初始化编码
    priority = len(priorityCode) # 初始化优先级
    # 复制矩阵,防止改变原值
    sourceFlowCopy = deepcopy(sourceFlow)
    destinationFlowCopy = deepcopy(destinationFlow)
    matCodeCopy = deepcopy(matCode)
    costMatCopy = deepcopy(costMat)
    largeNum = 1e5 # 设定一个足够大的数
    # 当流量没有完全分配时,执行迭代
    while not (np.any(sourceFlowCopy) and np.any(destinationFlowCopy)):
        # 为配送量为0的连接分配较大代价
        costMatCopy = np.where(matCodeCopy==0, largeNum, costMatCopy)
        i,j = np.where(costMatCopy = np.min(costMatCopy)) # 最小运输代价所对应的source和destination
        # 更新剩余source和destination流量
        sourceFlowCopy[i] -= matCodeCopy[i][j]
        destinationFlowCopy[j] -= matCodeCopy[i][j]
        # 当剩余sourceFlow或者destinationFlow为0时,分配优先度编码
        if sourceFlowCopy[i] == 0:
            priorityCode[i] = priority
            priority -= 1
        if destinationFlowCopy[j] == 0:
            priorityCode[j + len(sourceFlowCopy)] = priority
            priority -= 1
    # 为剩余位置填充优先级,因为该边上实质流量为0,因此事实上非有效边,可以随意分配流量
    perm = np.random.permutation(range(1, priority + 1))
    priorityCode = np.asarray(priorityCode)
    priorityCode[priorityCode==0] = perm
    return priorityCode


## 解码
def priorityDecoding(priorityCode, costMat=costMat, sourceFlow=sourceFlow, destinationFlow=destinationFlow):
    '''根据优先度编码解码回配送方案的矩阵编码
    输入:priorityCode -- 基于优先级的编码, 长度为len(sourceFlow) + len(destinationFlow)
    costMat -- 代价矩阵,给出由一个source到一个destination的代价
    sourceFlow, destinationFlow -- 每个source的输出能力和destination的接受能力
    输出:matCode -- 一个满足sourceFlow和destinationFlow的可行解矩阵编码'''
    # 初始化矩阵编码
    matCode = np.zeros((len(sourceFlow), len(destinationFlow)))
    # 复制矩阵,防止改变原值
    sourceFlowCopy = deepcopy(sourceFlow)
    destinationFlowCopy = deepcopy(destinationFlow)
    costMatCopy = np.array(deepcopy(costMat))
    priorityCodeCopy = deepcopy(priorityCode)
    largeNum = 1e5 # 设定一个足够大的数
    # 列出source Node和destination Node
    sourceNodeList = list(range(1, len(sourceFlow)+1))
    destinationNodeList = list(range(len(sourceFlow)+1,
                               len(destinationFlow)+len(sourceFlow)+1))
    nodeList = sourceNodeList + destinationNodeList
    while np.any(priorityCodeCopy):
        # 选择优先度最高的节点    
        nodeSelected = np.asarray(nodeList)[np.argmax(priorityCodeCopy)]
        # 为剩余流量为0的行和列分配一个大数作为运输代价
        rowIdx = [i for i in range(len(sourceFlowCopy)) if sourceFlowCopy[i]==0]
        colIdx = [i for i in range(len(destinationFlowCopy)) if destinationFlowCopy[i]==0]
        for row in rowIdx:
            costMatCopy[row,:] = [largeNum]*len(costMatCopy[row,:])
        for col in colIdx:
            costMatCopy[:,col] = [largeNum]*len(costMatCopy[:,col])
        # 如果选中的节点在供方,则从需方选择对应运输代价最小的节点
        if nodeSelected in sourceNodeList:
            sourceNodeIdx = nodeSelected -1 # 作为index,比节点标号小1
            destinationNodeIdx = np.argmin(costMatCopy[sourceNodeIdx,:])
        # 如果选中的节点在需方,则从供方选择对应运输代价最小的节点
        else:
            destinationNodeIdx = nodeSelected -1 - len(sourceFlow)
            sourceNodeIdx = np.argmin(costMatCopy[:,destinationNodeIdx])
        # 更新选中边上的流量
        matCode[sourceNodeIdx][destinationNodeIdx] = min(sourceFlowCopy[sourceNodeIdx],
                                                        destinationFlowCopy[destinationNodeIdx])
        # 更新剩余流量
        sourceFlowCopy[sourceNodeIdx] -= matCode[sourceNodeIdx][destinationNodeIdx]
        destinationFlowCopy[destinationNodeIdx] -= matCode[sourceNodeIdx][destinationNodeIdx]
        # 更新优先度
        if sourceFlowCopy[sourceNodeIdx] == 0:
            priorityCodeCopy[sourceNodeIdx] = 0
        if destinationFlowCopy[destinationNodeIdx] == 0:
            priorityCodeCopy[destinationNodeIdx + len(sourceFlowCopy)] = 0
    return matCode

## 评价函数
def evaluate(ind):
    matCode = priorityDecoding(ind) # 将个体优先度编码解码为矩阵编码
    return np.sum(np.sum(np.asarray(costMat)*np.asarray(matCode))),

## DEAP原生的PMX需要个体编码从0开始
def PMX(ind1, ind2):
    ind1Aligned = [ind1[i]-1 for i in range(len(ind1))]
    ind2Aligned = [ind2[i]-1 for i in range(len(ind2))]
    tools.cxPartialyMatched(ind1Aligned, ind2Aligned)
    ind1Recovered = [ind1Aligned[i]+1 for i in range(len(ind1Aligned))]
    ind2Recovered = [ind2Aligned[i]+1 for i in range(len(ind2Aligned))]
    ind1[:] = ind1Recovered
    ind2[:] = ind2Recovered
    return ind1,ind2

## 注册工具
toolbox = base.Toolbox()
toolbox.register('genPriority', priorityCoding, matEncode())
toolbox.register('individual', tools.initIterate, creator.Individual, toolbox.genPriority)
toolbox.register('evaluate', evaluate)
toolbox.register('select', tools.selTournament, tournsize=2)
toolbox.register('mate', PMX)
toolbox.register('mutate', tools.mutShuffleIndexes, indpb=0.5)

## 生成初始族群
toolbox.popSize = 100
toolbox.register('population', tools.initRepeat, list, toolbox.individual)
pop = toolbox.population(toolbox.popSize)

## 记录迭代数据
stats=tools.Statistics(key=lambda ind: ind.fitness.values)
stats.register('min', np.min)
stats.register('avg', np.mean)
stats.register('std', np.std)
hallOfFame = tools.HallOfFame(maxsize=1)

## 遗传算法参数
toolbox.ngen = 500
toolbox.cxpb = 0.8
toolbox.mutpb = 0.1


# ## 遗传算法主程序
pop,logbook=algorithms.eaSimple(pop, toolbox, cxpb=toolbox.cxpb, mutpb=toolbox.mutpb,
                   ngen=toolbox.ngen ,stats=stats, halloffame=hallOfFame, verbose=True)

计算结果如下：

# 输出结果
from pprint import pprint
bestInd = hallOfFame.items[0]
bestFitness = bestInd.fitness.values
print('最佳运输组合为：')
pprint(priorityDecoding(bestInd))
print('该运输组合的代价为：'+str(bestFitness))

# 画出迭代图
minFit = logbook.select('min')
avgFit = logbook.select('avg')
plt.plot(minFit, 'b-', label='Minimum Fitness')
plt.plot(avgFit, 'r-', label='Average Fitness')
plt.xlabel('# Gen')
plt.ylabel('Fitness')
plt.legend(loc='best')
# 计算结果：
#最佳运输组合为：
#array([[0., 0., 3., 0., 5.],
#       [3., 1., 0., 0., 0.],
#       [0., 0., 5., 7., 0.],
#       [0., 4., 2., 0., 0.]])
#该运输组合的代价为：(130.37945775,)

迭代过程如图所示：

Tansportation planning_iterations

遗传算法实践(八) 简单配送问题
简单配送问题描述基本配送计划模型是从产地往销地配送产品时，在满足产地的供应量和销地的需求量这一约束前提下，制定配...
遗传算法综述及简单应用实例
遗传算法综述及简单应用实例
遗传算法实践(九) VRP问题
VRP（Vehicle Routing Problem）问题描述在VRP问题中，假设有一个供求关系系统，车辆从仓...
遗传算法实践(三) 背包问题
背包问题描述背包问题(knapsack problem)是指从多种物品中选择几件物品装满背包。在不超过背包承受重...
遗传算法helloworld级别的python实现（结果可视化）
问题描述：用遗传算法求使得F(X)最大的X,问题来源：莫烦的python教程之遗传算法最终效果：
模拟退火算法(附代码)
上一篇我们讲了旅行商问题及遗传算法来解决此类问题：遗传算法(附代码)[https://www.jianshu.co...
遗传算法实践(二) 配词问题
在非线性函数寻优中我们看到了编码方式对于问题求解的重要性，但是遗传算法中的评价函数设定也对求解有至关重要的作用。如...
遗传算法实践(十) VRPTW问题求解
问题描述车辆配送模型(Vehicle routing problem)是指从配送中心用车辆把物资配送给顾客时，规...
遗传算法简单介绍与MATLAB实现（二）
遗传算法简单介绍与MATLAB实现（二）引入题目一上一篇文章中我们简单的介绍了一了一下遗传算法，其中提到了多元...
旅行商问题的近似最优解（局部搜索、模拟退火、遗传算法）
旅行商问题的近似最优解（局部搜索、模拟退火、遗传算法）关键字：旅行商问题，TSP，局部搜索，模拟退火，遗传算法 ...