简单配送问题描述

基本配送计划模型是从产地往销地配送产品时，在满足产地的供应量和销地的需求量这一约束前提下，制定配送总成本最低的配送计划。即

简单配送问题示例

考虑有4个产地和5个销地的配送计划问题，供应量和需求量如下:

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

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

个体编码

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

设有个产地，个销地，为产地的供应，为销地的需求，为从产地运输单位数量商品到销地的代价，为从产地运往销地的商品数量。

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

• Step 0: 初始化，优先级初始化为，；
• Step 1: 从供应量和需求量中有一个不为0的产-销地组合中，找到运输代价最小的，并为其分配最大运输量(该组合对应的供应量和需求量中较小者)，；
• Step 2: 更新供应量和需求量，，；
• Step 3: 根据更新后的供应量和需求量分配优先级

• Step 4: 重复Step 1-3直到
• Step 5: 为染色体中仍然为0的剩余位置随机填充优先值

解码

解码过程如下：

• Step 0： 初始化各产销地配对的运输量；
• Step 1： 选择优先度最高的节点；
• Step 2： 如果该节点是销地，则从产地选择运输代价最小的节点，；如果该节点是产地，则从销地选择运输代价最小的节点，；
• Step 3：为选择的产-销地更新可用的最大运输量；
• Step 4：更新剩余的需求量和销量，；
• Step 5： 如果完成了需求量和销量的分配，将优先度归零。，
• 重复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