Frank-wolfe算法多OD对matlab实现

Frank-wolfe算法原理

在无约束最优化问题的基础上，我们可以进一步来求解约束最优化问题。约束最优化问题的一般形式为：

minf(x)

s.t.gi(x)≥0,i=1,…,m

先考虑gi(x)均为线性函数的情况，此时问题与线性规划的约束条件相同，仅仅是目标函数变成了非线性的。我们可以用泰勒展开对目标函数进行近似，将它线性化。将f(x)在xk处展开，有

f(x)≈f(xk)+∇f(xk)T(x−xk)

故原问题近似于

minf(x)≈f(xk)+∇f(xk)T(x−xk)

s.t.x∈S

其中S为线性约束构成的可行域。去掉常量后，问题可以写为

minf(x)≈∇f(xk)Tx

s.t.x∈S

设此问题的最优解为yk，则直观上dk=yk−xk 应当为原问题的可行下降方向。沿着此方向做一维搜索则可进行一次迭代。为了防止一维搜索的结果超出可行域，我们限制步长0≤λ≤1。注意到线性规划的可行域为凸集，由于xk和yk均为可行点，它们确定的连线均在可行域中。限制步长0≤λ≤1保证了一维搜索的结果在可行域中。

更多FW算法原理相关内容，参考

• 线形约束最优化问题的Frank-Worlfe算法
• frankwolfe算法

Frank-wolfe算法流程

算例

将道路网络抽象为图

对于以下网格，给出网格上各点间的OD需求，计算网格上的交通平衡分布

给定OD对

起始点O 1 1 1 4 4 4 11 11 11 8 8 目的地D 4 11 8 1 11 8 1 4 8 1 4 交通量 1200 1400 1600 1400 1600 1200 1400 1200 1400 1200 1000

搜索每个OD对在网络上的可行径

% 子程序：求解一个OD对间的可行径function possiablePaths = findPath(Graph, partialPath, destination, partialWeight)% findPath按深度优先搜索所有可能的从partialPath出发到destination的路径，这些路径中不包含环路% Graph: 路网图，非无穷或0表示两节点之间直接连通，矩阵值就为路网权值% partialPath: 出发的路径，如果partialPath就一个数，表示这个就是起始点% destination: 目标节点% partialWeight: partialPath的权值，当partialPath为一个数时，partialWeight为0pathLength = length(partialPath);lastNode = partialPath(pathLength); %得到最后一个节点nextNodes = find(0<Graph(lastNode,:) & Graph(lastNode,:)<inf); %根据Graph图得到最后一个节点的下一个节点GLength = length(Graph);possiablePaths = [];if lastNode == destination % 如果lastNode与目标节点相等，则说明partialPath就是从其出发到目标节点的路径，结果只有这一个，直接返回 possiablePaths = partialPath; possiablePaths(GLength + 1) = partialWeight; return;elseif length( find( partialPath == destination ) ) ~= 0 return;end%nextNodes中的数一定大于0,所以为了让nextNodes(i)去掉，先将其赋值为0for i=1:length(nextNodes) if destination == nextNodes(i)  %输出路径  tmpPath = cat(2, partialPath, destination);      %串接成一条完整的路径  tmpPath(GLength + 1) = partialWeight + Graph(lastNode, destination); %延长数组长度至GLength+1, 最后一个元素用于存放该路径的总路阻  possiablePaths( length(possiablePaths) + 1 , : ) = tmpPath;  nextNodes(i) = 0; elseif length( find( partialPath == nextNodes(i) ) ) ~= 0  nextNodes(i) = 0; endendnextNodes = nextNodes(nextNodes ~= 0); %将nextNodes中为0的值去掉，因为下一个节点可能已经遍历过或者它就是目标节点for i=1:length(nextNodes) tmpPath = cat(2, partialPath, nextNodes(i)); tmpPsbPaths = findPath(Graph, tmpPath, destination, partialWeight + Graph(lastNode, nextNodes(i))); possiablePaths = cat(1, possiablePaths, tmpPsbPaths);end

Frank-worlfe算法构造

function [e,Xn,td] = Frank_wolfe(Q,W,Cmax,Mxf)%% 1 给定路网数据，OD需求，路段能力% 计算网络上已知OD集，初始路阻，道路容量，路径路段关系%==========================================================================% Q为OD需求，第一行为O，第二行为D，第三行为OD需求% Cmax为路段能力，是一个节点数乘节点数的矩阵% mxf为路径路段0-1关系，是一个元胞行向量，元素数量为OD数，每一个成员是一个OD对应的路径路段关系ODnum = size(Q,2);W(W == inf) = 1000000;%==========================================================================%% 2 自动求出路径和路段数量，根据路段数量定义路段名，给定初始数据%==========================================================================numf = zeros(1,ODnum);for i = 1:ODnum    numf(i) = size(Mxf{1,i},1); % 计算路径数endnumx = size(Mxf{1,1},2); % 计算路段数n = sqrt(numx);syms lambda realx = sym('x',[1,numx]); % 根据路段数定义路段名cont=0;e=inf;x=x(1:numx); % 路段上的交通量X0=zeros(1,numx); % 路段上的交通量 数值解t=zeros(1,numx); % 路段走行函数%==========================================================================%% 3 构造阻抗函数并求出初始阻抗，此处用BPR函数%=======================================================t=W.*(1+0.15*(x./Cmax).^4);            %路段走行时间函数tt=t;t=W.*(1+0.15*(X0./Cmax).^4);t(isnan(t)) = 1000000;for i = 1:n    t((i-1)*n + i) = 0;endCkrs = cell(1,ODnum);for i = 1:ODnum    Ckrs{1,i} = (Mxf{1,i} * t');    Ckrs{1,i} = Ckrs{1,i}';end%=========================================================%% 4 全有全无配流%=========================================================Min = zeros(ODnum);index = zeros(ODnum);for i = 1:ODnum    [Min(i),index(i)]=min(Ckrs{1,i});endX1 = zeros(1,numx);for i = 1:ODnum    tempmatrix = Mxf{1,i};    X1=tempmatrix(index(i),:).*Q(3,i) + X1;                    %全有全无法为最短路径上的路段分配流量end%=========================================================%% 5 数据更新%=========================================================while e>0.001                          %精度判断    cont=cont+1;                       %迭代次数更新    t=(W).*(1+0.15*(X1./Cmax).^4);     %路段时间跟新    td = t;    t(isnan(t)) = 1000000;    for i = 1:n        t((i-1)*n + i) = 0;    end    for i = 1:ODnum        Ckrs{1,i} = (Mxf{1,i} * t');  %路径时间更新        Ckrs{1,i} = Ckrs{1,i}';    end    Min = zeros(ODnum);    index = zeros(ODnum);    for i = 1:ODnum        [Min(i),index(i)]=min(Ckrs{1,i});    end    %全有全无法求辅助流量    Y1 = zeros(1,numx);    for i = 1:ODnum        tempmatrix = Mxf{1,i};        Y1=tempmatrix(index(i),:).*Q(3,i) + Y1;                    %全有全无法为最短路径上的路段分配流量    end    %Y1=Mxf(index,:).*Q;                    S=Y1-X1;                           %搜索方向    if sum(S) < 0        break;    end    X2=X1+lambda*S;                    %先将X2用X1和lambda进行表达    t=(W).*(1+0.15*(X2./Cmax).^4);     %含lambda的阻抗表达    t(isnan(t)) = 1000000;    f=sum(S.*t,2);                     %2表示按行求和    lambda1 = 0;    lambda1=double(solve(f));          %求解方程，确定步长。    k=length(lambda1);                 %如步长lambda1的解不唯一，取实数，且大于0小于1；    for m=1:k        if lambda1(m,1)>=0&&lam        bda1(m,1)<=1            lambda2=lambda1(m,1);        end    end    X2=X1+lambda2*S;                   %得到下一步的流量值，且进行下一次迭代    e=sqrt(sum((X2-X1).^2))/sum(X1);   %精度计算    X1=X2;                             %流量更新    disp(['迭代次数',num2str(cont),'精度',num2str(e)]);end%==========================================================================Xn = X1;

存在的问题

当网络中的交通量不大时，在迭代计算中利用sovle函数求解λ时，会产生不规律复现求解结果为两个复数的情况，目前认为应该是算法的构建中还有数学思想不够完善的地方导致。
进一步完善，敬请期待。

参考博文：配流07—基于BPR函数的Frank Wolfe算法