免疫优化算法在物流配送中心选址中的应用(matlab实现)

12.1 理论基础

12.1.1 物流中心选址问题

        随着世界经济的快速发展以及现代科学技术的进步，物流业作为国民经济的一个新兴服务部门，正在全球范围内迅速发展。物流业的发展给社会的生产和管理、人们的生活和就业乃至政府的职能以及社会的法律制度等都带来巨大的影响，因此物流业被认为是国民经济发展的动脉和基础产业，被形象地喻为促进经济发展的“加速器”。

        在物流系统的运作中，配送中心的任务就是根据各个用户的需求及时、准确和经济地配送商品货物。配送中心是连接供应商与客户的中间桥梁，其选址方式往往决定着物流的配送距离和配送模式，进而影响着物流系统的运作效率。另外，物流中心的位置一旦被确定，其位置难以再改变。因此，研究物流配送中心的选址具有重要的理论意义和现实应用意义。一般说来，物流中心选址模型是非凸和非光滑的带有复杂约束的非线性规划模型，属NP-hard问题。

12.1.2 免疫算法的基本思想

        生物免疫系统是一个高度进化的生物系统，它旨在区分外部有害抗原和自身组织，从而保持有机体的稳定。从计算角度看，生物免疫系统是一个高度并行、分布、自适应和自组织的系统，具有很强的学习、识别和记忆能力。

        免疫系统具有如下特征：

        (1)产生多样抗体的能力。通过细胞的分裂和分化作用，免疫系统可产生大量的抗体来抵御各种抗原。

        (2)自我调节机构。免疫系统具有维持免疫平衡的机制，通过对抗体的抑制和促进作用，能自我调节产生适当数量的必要抗体。

        (3)免疫记忆功能。产生抗体的部分细胞会作为记忆细胞被保存下来，对于今后侵入的同类抗原，相应的记忆细胞会迅速激发而产生大量的抗体。

        免疫算法(immune algorithm)正是受生物免疫系统启发，在免疫学理论基础上发展起来的一种新兴的智能计算方法。它利用免疫系统的多样性产生和维持机制来保持群体的多样性，克服了一般寻优过程尤其是多峰函数寻优过程中难处理的“早熟”问题，最终求得全局最优解。与其他智能算法相比，免疫算法的研究起步较晚，其发展历史只有短短二十几年。Farmer等于1986年率先基于免疫网络学说构造了免疫系统的动态模型，并探讨了免疫系统与其他人工智能方法的联系，从而开创了免疫系统的研究。

        免疫算法和遗传算法都是采用群体搜索策略，并且强调群体中个体间的信息交换，因此有许多相似之处，比如两者具有大致相同的算法结构，都要经过“初始种群产生→评价标准计算 →种群间个体信息交换→新种群产生”这一循环过程，最终以较大的概率获得问题的最优解； 另外，两者本质上具有并行性，具有与其他智能计算方法结合的固有优势等。 免疫算法和遗传算法之间也存在一些区别，主要表现为对个体的评价、选择及产生的方式不同。遗传算法中个体的评价是通过计算个体适应度得到的，算法选择父代个体的唯一标准是个体适应度；而免疫算法对个体的评价则是通过计算亲和度(affinity)得到的，个体的选择也是以亲和度为基础进行的。个体的亲和度包括抗体-抗原之间的亲和度(匹配程度)和抗体- 抗体之间的亲和度(相似程度),它反映了真实的免疫系统的多样性，因此免疫算法对个体的评价更加全面，其个体选择方式也更为合理。此外，遗传算法通过交叉、变异等遗传操作产生新个体，而在免疫算法中，虽然交叉、变异等固有的遗传操作也被广泛应用，但是新抗体的产生还可以借助克隆选择、免疫记忆、疫苗接种等遗传算法中所欠缺的机理，同时免疫算法中还对抗体的产生进行促进或者抑制，体现了免疫反应的自我调节功能，保证了个体的多样性。

        本案例把免疫优化算法用于物流配送中心选址问题中。在考虑该问题的约束条件和优化目标的基础上，建立了物流配送中心选址问题的数学模型，并采用免疫优化算法求解最佳物流配送中心选址模型

12.2 案例背景

12.2.1 问题描述

        在物流配送中心选址模型中做如下假设：

        (1)配送中心的规模容量总可以满足需求点需求，并由其配送辐射范围内的需求量确定；

        (2)一个需求点仅由一个配送中心供应；

        (3)不考虑工厂到配送中心的运输费用。

        基于以上假设，建立如下模型。该模型是一个选址/分配模型，在满足距离上限的情况下，需要从n个需求点中找出配送中心并向各需求点配送物品。目标函数是各配送中心到需求点的需求量和距离值的乘积之和最小，目标函数为

        式(12-2)保证每个需求点只能由一个配送中心服务；式(12-3)确保需求点的需求量只能被设为配送中心的点供应，即没有配送中心的地点不会有客户；式(12-4)规定了被选为配送中心的数量为p;式(12-5)表示变量Zij和hj是0-1变量；式(12-6)保证了需求点在配送中心可配送到的范围内。

12.2.2 解题思路及步骤

        1.算法流程

        免疫算法流程如图12-1所示。

        免疫算法具体实现步骤如下：

        (1)分析问题。对问题及其解的特性进行分析，设计解的合适表达形式。

        (2)产生初始抗体群。随机产生N个个体并从记忆库中提取m个个体构成初始群体，其中m为记忆库中个体的数量。

        (3)对上述群体中各个抗体进行评价。在本算法中对个体的评价是以个体的期望繁殖率P为标准的。

        (4)形成父代群体。将初始群体按期望繁殖率P进行降序排列，并取前N个个体构成父代群体；同时取前m个个体存入记忆库中。

        (5)判断是否满足结束条件，是则结束；反之，则继续下一步操作。

        (6)新群体的产生。基于步骤(4)的计算结果对抗体群体进行选择、交叉、变异操作得到新群体，再从记忆库中取出记忆的个体，共同构成新一代群体。

        (7)转去执行步骤(3)。

        2.初始抗体群的产生

        如果记忆库非空，则初始抗体群从记忆库中选择生成。否则，在可行解空间随机产生初始抗体群。此处采用简单编码方式。每个选址方案可形成一个长度为p的抗体(p表示配送中心数量),每个抗体代表被选为配送中心的需求点的序列。例如，考虑包含31个需求点的问题。1,2,…,31代表需求点的序号。从中选出6个作为配送中心。抗体[2 7 15 21 29 11]代表一个可行解，它表示2,7,15,21,29,11被选为配送中心。这种编码方式能够满足约束条件式(12-4)和式(12-5)。

        3.解的多样性评价

        (1)抗体与抗原间亲和力

        抗体与抗原之间的亲和力用于表示抗体对抗原的识别程度，此处针对上述配送中心选址模型设计亲和力函数Av

        (2)抗体与抗体间亲和力

        抗体与抗体之间的亲和力反映了抗体之间的相似程度。此处借鉴由Forrest等提出的R位连续方法计算抗体与抗体间的亲和力。R位连续方法实际是一种部分匹配规则。该方法的关键是确定一个R值，代表亲和度判定的阈值。两种个体编码有超过R位或者连续R位的编码相同，则表示这两种抗体近似“相同”,否则表示两种个体不同。此处抗原的编码方法，各位之间不需考虑排序，可参考变形的R位连续方法计算抗体间亲和度，即

        免疫算法在抑制高浓度个体时，与抗原亲和度最高的个体也可能因其浓度高而受到抑制，从而导致已求得的最优解丢失，因此采取精英保留策略，在每次更新记忆库时，先将与抗原亲和度最高的若干个个体存入记忆库，再按照期望繁殖概率将剩余群体中优秀个体存入记忆库。

4.免疫操作

        (1)选择：按照轮盘赌选择机制进行选择操作，个体被选择的概率即为式(12-10)计算出的期望繁殖概率。

        (2)交叉：本文采用单点交叉法进行交叉操作。

        (3)变异：采用常用的变异方法，即随机选择变异位进行变异。

12.3 MATLAB程序实现

        在MATLAB软件中编程实现免疫优化算法求解物流配送中心选址问题，代码如下。

12.3.1 免疫算法主函数

%% 免疫优化算法在物流配送中心选址中的应用
%% 清空环境
clc
clear

%% 算法基本参数           
sizepop=50;           % 种群规模
overbest=10;          % 记忆库容量
MAXGEN=100;            % 迭代次数
pcross=0.5;           % 交叉概率
pmutation=0.4;        % 变异概率
ps=0.95;              % 多样性评价参数
length=6;             % 配送中心数
M=sizepop+overbest;

%% step1 识别抗原,将种群信息定义为一个结构体
individuals = struct('fitness',zeros(1,M), 'concentration',zeros(1,M),'excellence',zeros(1,M),'chrom',[]);
%% step2 产生初始抗体群
individuals.chrom = popinit(M,length);
trace=[]; %记录每代最个体优适应度和平均适应度

%% 迭代寻优
for iii=1:MAXGEN

     %% step3 抗体群多样性评价
     for i=1:M
         individuals.fitness(i) = fitness(individuals.chrom(i,:));      % 抗体与抗原亲和度(适应度值）计算
         individuals.concentration(i) = concentration(i,M,individuals); % 抗体浓度计算
     end
     % 综合亲和度和浓度评价抗体优秀程度，得出繁殖概率
     individuals.excellence = excellence(individuals,M,ps);
          
     % 记录当代最佳个体和种群平均适应度
     [best,index] = min(individuals.fitness);   % 找出最优适应度 
     bestchrom = individuals.chrom(index,:);    % 找出最优个体
     average = mean(individuals.fitness);       % 计算平均适应度
     trace = [trace;best,average];              % 记录
     
     %% step4 根据excellence，形成父代群，更新记忆库（加入精英保留策略，可由s控制）
     bestindividuals = bestselect(individuals,M,overbest);   % 更新记忆库
     individuals = bestselect(individuals,M,sizepop);        % 形成父代群

     %% step5 选择，交叉，变异操作，再加入记忆库中抗体，产生新种群
     individuals = Select(individuals,sizepop);                                                             % 选择
     individuals.chrom = Cross(pcross,individuals.chrom,sizepop,length);                                    % 交叉
     individuals.chrom = Mutation(pmutation,individuals.chrom,sizepop,length);   % 变异
     individuals = incorporate(individuals,sizepop,bestindividuals,overbest);                               % 加入记忆库中抗体      

end

%% 画出免疫算法收敛曲线
figure(1)
plot(trace(:,1));
hold on
plot(trace(:,2),'--');
legend('最优适应度值','平均适应度值')
title('免疫算法收敛曲线','fontsize',12)
xlabel('迭代次数','fontsize',12)
ylabel('适应度值','fontsize',12)

%% 画出配送中心选址图
%城市坐标
city_coordinate=[1304,2312;3639,1315;4177,2244;3712,1399;3488,1535;3326,1556;3238,1229;4196,1044;4312,790;4386,570;
                 3007,1970;2562,1756;2788,1491;2381,1676;1332,695;3715,1678;3918,2179;4061,2370;3780,2212;3676,2578;
                 4029,2838;4263,2931;3429,1908;3507,2376;3394,2643;3439,3201;2935,3240;3140,3550;2545,2357;2778,2826;2370,2975];
carge=[20,90,90,60,70,70,40,90,90,70,60,40,40,40,20,80,90,70,100,50,50,50,80,70,80,40,40,60,70,50,30];
%找出最近配送点
for i=1:31
    distance(i,:)=dist(city_coordinate(i,:),city_coordinate(bestchrom,:)');
end
[a,b]=min(distance');

index=cell(1,length);

for i=1:length
%计算各个派送点的地址
index{i}=find(b==i);
end
figure(2)
title('最优规划派送路线')
cargox=city_coordinate(bestchrom,1);
cargoy=city_coordinate(bestchrom,2);
plot(cargox,cargoy,'rs','LineWidth',2,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor','b',...
    'MarkerSize',20)
hold on

plot(city_coordinate(:,1),city_coordinate(:,2),'o','LineWidth',2,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor','g',...
    'MarkerSize',10)

for i=1:31
    x=[city_coordinate(i,1),city_coordinate(bestchrom(b(i)),1)];
    y=[city_coordinate(i,2),city_coordinate(bestchrom(b(i)),2)];
    plot(x,y,'c');hold on
end

12.3.2 多样性评价

        1.适应度计算

        函数Fitness用于计算个体的适应度值：

function fit=fitness(individual)
%计算个体适应度值
%individual    input      个体
%fit           output     适应度值
%城市坐标
city_coordinate=[1304,2312;3639,1315;4177,2244;3712,1399;3488,1535;3326,1556;3238,1229;4196,1044;4312,790;4386,570;
                 3007,1970;2562,1756;2788,1491;2381,1676;1332,695;3715,1678;3918,2179;4061,2370;3780,2212;3676,2578;
                 4029,2838;4263,2931;3429,1908;3507,2376;3394,2643;3439,3201;2935,3240;3140,3550;2545,2357;2778,2826;2370,2975];
%货物量
carge=[20,90,90,60,70,70,40,90,90,70,60,40,40,40,20,80,90,70,100,50,50,50,80,70,80,40,40,60,70,50,30];

%找出最近配送点
for i=1:31
    distance(i,:)=dist(city_coordinate(i,:),city_coordinate(individual,:)');
end

[a,b]=min(distance');

%计算费用
for i=1:31  
    expense(i)=carge(i)*a(i);
end

%距离大于3000取一个惩罚值

fit=sum(expense) + 4.0e+4*length(find(a>3000));


end

        2.相似度计算

        函数similar计算个体之间的相似程度：

function resemble = similar(individual1,individual2)
% 计算个体individual1和individual2的相似度
% individual1,individual2    input     两个个体
% resemble                   output     相似度

k=zeros(1,length(individual1));
for i=1:length(individual1)
    if find(individual1(i)==individual2)
        k(i)=1;
    end
end

resemble=sum(k)/length(individual1);

end

        3.浓度计算

        函数concentration用于计算个体之间的浓度：

function concentration = concentration(i,M,individuals)
% 计算个体浓度值
% i              input      第i个抗体
% M              input      种群规模
% individuals    input     个体
% concentration  output     浓度值

concentration=0;
for j=1:M
    xsd=similar(individuals.chrom(i,:),individuals.chrom(j,:));  % 第i个体与种群个体间的相似度
    % 相似度大于阀值
    if xsd>0.7
        concentration=concentration+1;
    end
end

concentration=concentration/M;

end

        4.期望繁殖概率计算

        函数excellence计算个体繁殖概率：

function exc=excellence(individuals,M,ps)
% 计算个体繁殖概率
% individuals    input      种群
% M              input      种群规模
% ps             input      多样性评价参数
% exc            output     繁殖概率

fit = 1./individuals.fitness;
sumfit = sum(fit);
con = individuals.concentration;
sumcon = sum(con);
for i=1:M
    exc(i) = fit(i)/sumfit*ps +con(i)/sumcon*(1-ps); 
end

end

12.3.3 免疫操作

        1.选择函数根据个体适应度值采用轮盘赌法选择个体

function ret=Select(individuals,sizepop)
% 轮盘赌选择
% individuals input  : 种群信息
% sizepop     input  : 种群规模
% ret         output : 选择后得到的种群

excellence=individuals.excellence;
pselect=excellence./sum(excellence);
% 事实上 pselect = excellence；

index=[]; 
for i=1:sizepop   % 转sizepop次轮盘
    pick=rand;
    while pick==0    
        pick=rand;        
    end
    for j=1:sizepop   
        pick=pick-pselect(j);        
        if pick<0        
            index=[index j];
            break;  % 寻找落入的区间，此次转轮盘选中了染色体j
        end
    end
end
% 注意：在转sizepop次轮盘的过程中，有可能会重复选择某些染色体

individuals.chrom=individuals.chrom(index,:);
individuals.fitness=individuals.fitness(index);
individuals.concentration=individuals.concentration(index);
individuals.excellence=individuals.excellence(index);
ret=individuals;

end
 

        2.交叉操作采用实数交叉法进行交叉

function ret=Cross(pcross,chrom,sizepop,length)
% 交叉操作
% pcorss                input  : 交叉概率
% chrom                 input  : 抗体群
% sizepop               input  : 种群规模
% length                input  : 抗体长度
% ret                   output : 交叉得到的抗体群

% 每一轮for循环中，可能会进行一次交叉操作，随机选择染色体是和交叉位置，是否进行交叉操作则由交叉概率（continue）控制
for i=1:sizepop  
    
    % 随机选择两个染色体进行交叉
    pick=rand;
    while prod(pick)==0
        pick=rand(1);
    end
    
    if pick>pcross
        continue;
    end
    
    % 找出交叉个体
    index(1)=unidrnd(sizepop);
    index(2)=unidrnd(sizepop);
    while index(2)==index(1)
        index(2)=unidrnd(sizepop);
    end
    
    % 选择交叉位置
    pos=ceil(length*rand);
    while pos==1
        pos=ceil(length*rand);
    end

    % 个体交叉
    chrom1=chrom(index(1),:);
    chrom2=chrom(index(2),:);
    
    k=chrom1(pos:length);
    chrom1(pos:length)=chrom2(pos:length);
    chrom2(pos:length)=k; 
    
    % 满足约束条件赋予新种群
    flag1=test(chrom(index(1),:));
    flag2=test(chrom(index(2),:));
    
    if flag1*flag2==1
        chrom(index(1),:)=chrom1;
        chrom(index(2),:)=chrom2;
    end
    
end

ret=chrom;
end

        3.变异操作采用实数变异法进行变异

function ret=Mutation(pmutation,chrom,sizepop,length1)
% 变异操作
% pmutation        input  : 变异概率
% chrom            input  : 抗体群
% sizepop          input  : 种群规模
% iii              input  : 进化代数
% MAXGEN           input  : 最大进化代数
% length1          input  : 抗体长度
% ret              output : 变异得到的抗体群
% 每一轮for循环中，可能会进行一次变异操作，染色体是随机选择的，变异位置也是随机选择的
for i=1:sizepop   
    
    % 变异概率
    pick=rand;
    while pick==0
        pick=rand;
    end
    index=unidrnd(sizepop);

   % 判断是否变异
    if pick>pmutation
        continue;
    end
    
    pos=unidrnd(length1);
    while pos==1
        pos=unidrnd(length1);
    end
    
    nchrom=chrom(index,:);
    nchrom(pos)=unidrnd(31);
    while length(unique(nchrom))==(length1-1)
        nchrom(pos)=unidrnd(31);
    end
    
    flag=test(nchrom);
    if flag==1
        chrom(index,:)=nchrom;
    end
    
end

ret=chrom;
end
 

        4.产生新种群

function newindividuals = incorporate(individuals,sizepop,bestindividuals,overbest)
% 将记忆库中抗体加入，形成新种群
% individuals         input          抗体群
% sizepop             input          抗体数
% bestindividuals     input          记忆库
% overbest            input          记忆库容量

m = sizepop+overbest;
newindividuals = struct('fitness',zeros(1,m), 'concentration',zeros(1,m),'excellence',zeros(1,m),'chrom',[]);

% 遗传操作得到的抗体
for i=1:sizepop
    newindividuals.fitness(i) = individuals.fitness(i);   
    newindividuals.concentration(i) = individuals.concentration(i);   
    newindividuals.excellence(i) = individuals.excellence(i);   
    newindividuals.chrom(i,:) = individuals.chrom(i,:);   
end
% 记忆库中抗体
for i=sizepop+1:m
    newindividuals.fitness(i) = bestindividuals.fitness(i-sizepop);   
    newindividuals.concentration(i) = bestindividuals.concentration(i-sizepop);   
    newindividuals.excellence(i) = bestindividuals.excellence(i-sizepop);   
    newindividuals.chrom(i,:) = bestindividuals.chrom(i-sizepop,:);   
end

end

12.3.4 仿真实验

        为证明算法的可行性和有效性，采集了全国31个城市的坐标，每个用户的位置及其物资 需求量由表12-1中给出，这里的物资需求量是经过规范化处理后的数值，并不代表实际值。从中选择6个作为物流配送中心。根据配送中心选址模型，按照免疫算法步骤对算例进行求解，算法的参数分别为：种群规模为50,记忆库容量为10,迭代次数为100,交叉概率为0.5,变异概率为0.4,多样性评价参数设为0.95,求得配送中心的选址方案为[18 25 5 27 9 14],此方案下以各需求点需求量为权重的距离和为5.68×10^5。

        免疫算法收敛曲线如图12-2所示。得到的物流配送中心选址方案如图12-3所示。

        图12-3中方框表示配送中心，圆点表示城市点，若点间有连线表示该城市点的物资由连接的配送中心配送。

12.4 案例扩展

        人工免疫算法的研究还处于起步阶段，并且由于免疫机理复杂，系统庞大，人工免疫系统可借鉴的成果不多，免疫算法在系统建模，算法等方面存在诸多问题。目前，免疫算法的研究热点包括如下几个方面：

        (1)免疫算法性能研究。免疫算法性能包括收敛性、动态性能以及有效性，目前研究较少，研究的有效途径可以参考进化算法和人工神经网络的相关研究成果。

        (2)免疫算法与进化算法等其他智能算法的对比研究。该方法不仅能够深入认识人工免疫系统特点，还可以促进多种智能策略的互补融合。

        (3)免疫算法在网络、智能系统和鲁棒系统中的应用。神经网络、内分泌及免疫这三大调节系统相互联系、相互补充和配合、相互制约的机理为基于人工免疫系统的智能综合集成提供了生物学基础，网络和智能成为免疫算法发展的不可缺少的特征，也是其重要应用领域。免疫算法能增强系统的鲁棒性，而且免疫性与鲁棒性之间存在的必然联系使得免疫算法将在鲁棒系统中得到较好的应用。