【物流选址】基于matlab佛洛依德算法求解物流选址问题【含Matlab源码 892期】

一、物流选址简介

1 LRP问题本质
LRP问题的本质是为了满足区域类客户的配送需求，需要建立一个多级配送网络，一般由配送枢纽、配送中心、配送末端网点三级构成（其他类型可类比）。要求满足客户的需求量，不同设施有其吞吐容量和服务能力约束，不同级设施之间配送单位费用均不同，客户也有配送的时间约束。在满足多种约束下，求整个配送网络总成本最小的选址和需求分配方案。

2 LRP重难点
LRP问题求解的重难点一般是两个：选址和需求分配。选址的难度是在不同级设施备选点中选择合适数量和位置的设施组合；需求分配的难度是在不同分配路径方案中选择不同路径方案的组合，其中，因为每个客户的需求量从源头流经多级配送设施的流量是未知的，而且这个未知的变量是MGH个的。（M,G,H分别是1，2，3级配送设施的备选点数量），流量取值是连续的变量。因此，通过转化为多变量的目标函数求最值的线性或者非线性函数求解是不可能做到的。

3 LRP思路
因此，如何克服以上难点成为了我们的关键问题。首先，选址的找到一个最有的排列组合，可以通过万能钥匙的遗传算法求解。然后需求分配则需要按照最小费用最大流问题的思路进行分配。最小费用最大流问题指的是，在一个由一个源和汇和中间节点构成的配送网络中，每个路径都有其成本权重和流量约束，在满足流量约束的条件下找到最小配送成本的配送路径，这种方法就称为最小费用最大流方法。
借鉴了这种方法之后，再根据模型的约束进行改进，得到一个新的最小费用最大流方法。

二、弗洛伊德(Floyd)算法简介

1 弗洛伊德(Floyd)算法介绍
1)和Dijkstra算法一样，弗洛伊德(Floyd)算法也是一种用于寻找给定的加权图中顶点间最短路径的算法。该算法名称以创始人之一、1978年图灵奖获得者、斯坦福大学计算机科学系教授罗伯特.弗洛伊德命名
2)弗洛伊德算法(Floyd)计算图中各个顶点之间的最短路径
3)迪杰斯特拉算法用于计算图中某一个顶点到其他顶点的最短路径。
4)弗洛伊德算法VS迪杰斯特拉算法:迪杰斯特拉算法通过选定的被访问项点，求出从出发访问顶点到其他项点的最短路径;弗洛伊德算法中每个顶点都是出发访问点，所以需要将每一个顶点看做被访问顶点，求出从每一个顶点到其他项点的最短路径。

2 弗洛伊德(Floyd)算法最佳应用-最短路径

1)胜利乡有7个村庄(A,B,C,D,E,E, G)
2)各个村庄的距离用边线表示(权)，比如A-B距离5公里
3)问:如何计算出各村庄到其它各村庄的最短距离?

3 弗洛伊德(Floyd)算法图解分析
1)设置顶点vi到顶点vk的最短路径己知为Lik，顶点vk到vj的最短路径已知为Lkj，顶点vi到vj的路径为Lij，则vi到vj的最短路径为: min((Lik+Lkj),Lij)， vk的取值为图中所有顶点，则可获得vi到vj的最短路径
2)至于vi到vk的最短路径Lik或者vk到vj的最短路径Lkj,是以同样的方式获得
3)弗洛伊德(Floyd)算法图解分析举例说明
3.1 弗洛伊德算法的步骤:
初始状态：

第一轮循环中，以A(下标为: 0)作为中间顶点[即把A作为中间顶点的所有情况都进行遍历,
就会得到更新距离表和前驱关系]，距离表和前驱关系更新为:

将A做为中间顶点的情况有：
C->A->G:9,C->G:N
C->A->B:12,C->B:N
G->A->B:7,G->B:3
通过比较，得出最小的值，然后更新距离表和前驱关系表
以此类推。。。

三、案例及部分源代码

1 案例
某公司有6个建筑工地要开工，每个工地的位置（用平面坐标系a，b表示，距离单位：千米）及水泥日用量d（吨）由下表给出。目前急需新建两个料场，位置待定，日储量各有20吨。假设从料场到工地之间均有直线道路相连。试制定两个料场的位置以及每天的供应计划，即从A，B两料场分别向各工地运送多少吨水泥，使得总的运输吨千米数最小。

2 部分源代码

x0 = [3.0000  5.0000  0.0707  7.0000  0  0.9293  0  0  3.9293  0  6.0000 ...10.0707  6.3875  4.3943  5.7511  7.1867]';
A=[1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0];
b=[20;20];
Aeq=[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0;...0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0;...0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0;...0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0;...0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0;...0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0];
beq=[3;5;4;7;6;11];
VLB=[0;0;0;0;0;0;0;0;0;0;0;0];
VUB=[];
[x,fval]=fmincon('f1',x0,A,b,Aeq,beq,VLB,VUB)
function [x,val,k]=bfgs(fun,gfun,x0,varargin)
%功能: 用BFGS算法求解无约束问题:  min f(x)
%输入: x0是初始点, fun, gfun分别是目标函数及其梯度;
% varargin是输入的可变参数变量, 简单调用bfgs时可以忽略它,
% 但若其它程序循环调用该程序时将发挥重要的作用
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
maxk=500;   %给出最大迭代次数
rho=0.55; sigma1=0.4; epsilon1=1e-5;
k=0;   n=length(x0);
Bk=eye(n);   %Bk=feval('Hess',x0);
while(k<maxk)gk=feval(gfun,x0,varargin{:}); %计算梯度if(norm(gk)<epsilon1), break; end  %检验终止准则dk=-Bk\gk;  %解方程组, 计算搜索方向m=0; mk=0;while(m<20)   % 用Armijo搜索求步长 newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma1*rho^m*gk'*dk)mk=m; break;endm=m+1;end%BFGS校正x=x0+rho^mk*dk;  sk=x-x0;  yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;     x0=x;
end
function g=df1(x)
g = [sqrt((x(13)-1.25)^2+(x(14)-1.25)^2), sqrt((x(13)-8.75)^2+(x(14)-0.75)^2),...sqrt((x(13)-0.5)^2+(x(14)-4.75)^2), sqrt((x(13)-5.75)^2+(x(14)-5.0)^2),...sqrt((x(13)-3.0)^2+(x(14)-6.5)^2), sqrt((x(13)-7.25)^2+(x(14)-7.25)^2),...sqrt((x(15)-1.25)^2+(x(16)-1.25)^2), sqrt((x(15)-8.75)^2+(x(16)-0.75)^2),...sqrt((x(15)-0.5)^2+(x(16)-4.75)^2), sqrt((x(15)-5.75)^2+(x(16)-5.0)^2),...sqrt((x(15)-3.0)^2+(x(16)-6.5)^2), sqrt((x(15)-7.25)^2+(x(16)-7.25)^2)]';
g(13)=x(1)*(x(13)-1.25)/g(1)+x(2)*(x(13)-8.75)/g(2)+x(3)*(x(13)-0.5)/g(3)+...x(4)*(x(13)-5.75)/g(4)+x(5)*(x(13)-3.0)/g(5)+x(6)*(x(13)-7.25)/g(6);
g(14)=x(1)*(x(14)-1.25)/g(1)+x(2)*(x(14)-0.75)/g(2)+x(3)*(x(14)-4.75)/g(3)+...x(4)*(x(14)-5.0)/g(4)+x(5)*(x(14)-6.5)/g(5)+x(6)*(x(14)-7.25)/g(6);
g(15)=x(7)*(x(15)-1.25)/g(7)+x(8)*(x(15)-8.75)/g(8)+x(9)*(x(15)-0.5)/g(9)+...x(10)*(x(15)-5.75)/g(10)+x(11)*(x(15)-3.0)/g(11)+x(12)*(x(15)-7.25)/g(12);
g(16)=x(7)*(x(16)-1.25)/g(7)+x(8)*(x(16)-0.75)/g(8)+x(9)*(x(16)-4.75)/g(9)+...x(10)*(x(16)-5.0)/g(10)+x(11)*(x(16)-6.5)/g(11)+x(12)*(x(16)-7.25)/g(12);

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.