精华内容
下载资源
问答
  • e=0.005; X0=[0.1;0.2;0.7]; 符号 x1 x2 x3; f=0.5*(x1^2+x2^2+0.1*x3^2)+0.55*x3; [X,fval]=frank_wolfe(f,X0,e,[],[],[1 1 1],[1],[0;0;0],[1;1;1])
  • Frank-wolfe算法多OD对matlab实现

    万次阅读 多人点赞 2017-08-25 14:59:56
    Frank-wolfe算法多OD对matlab实现 Frank-wolfe算法多OD对matlab实现 Frank-wolfe算法原理 Frank-wolfe算法流程 算例 将道路网络抽象为图 给定OD对 关键函数及完整流程 1. 搜索每个OD对在网络上的可行径 2. ...

    Frank-wolfe算法多OD对matlab实现

    Frank-wolfe算法原理

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

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

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

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

    故原问题近似于

    minf(x)f(xk)+f(xk)T(xxk)minf(x)≈f(xk)+∇f(xk)T(x−xk)
    s.t.xSs.t.x∈S

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

    minf(x)f(xk)Txminf(x)≈∇f(xk)Tx
    s.t.xSs.t.x∈S

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

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

    Frank-wolfe算法流程

    FW算法

    算例

    将道路网络抽象为图

    在武汉地区选择一个大型小区,其路网经过梳理后如下图:(其中粗实线表示主干道,而次干道和支路并没有直接分开,可参考道路通行能力上的差异)
    这里写图片描述
    对于该网络,建立邻接矩阵转化为图,给出图上各点间的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

    关键函数及完整流程

    1. 搜索每个OD对在网络上的可行径
    % 子程序:求解一个OD对间的可行径
    function possiablePaths = findPath(Graph, partialPath, destination, partialWeight)
    % findPath按深度优先搜索所有可能的从partialPath出发到destination的路径,这些路径中不包含环路
    % Graph: 路网图,非无穷或0表示两节点之间直接连通,矩阵值就为路网权值
    % partialPath: 出发的路径,如果partialPath就一个数,表示这个就是起始点
    % destination: 目标节点
    % partialWeight: partialPath的权值,当partialPath为一个数时,partialWeight为0
    pathLength = 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)去掉,先将其赋值为0
    for 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;
     end
    end
    nextNodes = 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
    2. 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); % 计算路径数
    end
    numx = size(Mxf{1,1},2); % 计算路段数
    n = sqrt(numx);
    syms lambda real
    x = 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;
    end
    Ckrs = 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});
    end
    X1 = 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;
    3. 主函数
    clc;clear;
    %% 构建通行时间费用阻抗矩阵
    % 网络中的距离设置,inf表示两点之间无直接相连
    n = 11;distance = zeros(n);
    distance(1,2) = 145;distance(1,5) = 445;
    distance(2,3) = 262;distance(2,6) = 192;
    distance(3,4) = 375;distance(3,7) = 174;
    distance(4,11) = 468;
    distance(5,8) = 142;
    distance(6,[7,9]) = [264 278];distance(7,10) = 289;
    distance(8,9) = 334;distance(9,10) = 242;distance(10,11) = 369;
    
    
    distance = distance + distance';
    distance(distance == 0) = inf;
    distance(1:n+1:n^2) = 0;
    distance = distance/1000;% 对角线元素替换成0
    
    % 由道路等级决定的道路设计速度
    v0 = zeros(n);
    v0(1,2) = 50;v0(1,5) = 40;
    v0(2,3) = 50;v0(2,6) = 30;
    v0(3,4) = 50;v0(3,7) = 30;
    v0(4,11) = 50;
    v0(5,8) = 50;
    v0(6,[7,9]) = [30 30];v0(7,10) = 30;
    v0(8,9) =50;v0(9,10) = 50;v0(10,11) = 50;
    
    
    v0 = v0 + v0';
    v0(v0 == 0) = 0.1; 
    v0(1:n+1:n^2) = 0; % 对角线元素替换成0
    
    t0 = distance./v0 * 3600;
    t0(isnan(t0)) = 0;
    %t0 = t0(:);
    
    %% 设置OD矩阵
    OD = [1 1   1   4   4   4   11  11  11  8   8 8;
    4   11  8   1   11  8   1   4   8   1   4 11;
    1200    1400    1600    1400    1600    1200    1400    1200    1400    1200    1000 1300];
    % OD矩阵的第一行为O,第二行为D,第三行为该OD对上的OD值
    
    %% 设置道路容量
    Cmax = zeros(n);
    Cmax(1,2) = 50;Cmax(1,5) = 40;
    Cmax(2,3) = 50;Cmax(2,6) = 30;
    Cmax(3,4) = 50;Cmax(3,7) = 30;
    Cmax(4,11) = 50;
    Cmax(5,8) = 50;
    Cmax(6,[7,9]) = [30 30];Cmax(7,10) = 30;
    Cmax(8,9) =50;Cmax(9,10) = 50;Cmax(10,11) = 50;
    
    
    Cmax(Cmax == 50) = 1400;
    Cmax(Cmax == 40) = 1000;
    Cmax(Cmax == 30) = 700;
    Cmax = Cmax + Cmax';
    Cmax1 = Cmax(:)';
    Cmax1(Cmax1 == 0) = 0.01;
    % 构造关联矩阵的矩阵向量
    
    odnum = size(OD,2);
    mxf = cell(1,odnum);
    for i = 1:odnum
        O = OD(1,i);D = OD(2,i);
        a = findPath(t0,O,D,0);
        a = a(:,1:size(a,2)-1);
        waynum = size(a,1);
        way = zeros(waynum,n*n);
        waytemp = zeros(n);
        for j = 1:waynum % 对每条可行径
            for x = 1:size(a,2)-1
                if a(j,x+1) ~= 0 
                    waytemp(a(j,x),a(j,x+1)) = 1;
                else
                    break;
                end
            end
            way(j,:) = waytemp(:);
            waytemp = zeros(n);
        end
        mxf{1,i} = way;
    end
    
    % 计算交通在道路网络的分布
    tc = t0;
    t0 = t0(:)';
    [e,Xn,td] = Frank_wolfe(OD,t0,Cmax1,mxf);
    X = zeros(n);
    for i = 1:n
        X(:,i) = Xn(((i-1)*n+1):i*n)';
    end
    for i = 1:n
        tdd(:,i) = td(((i-1)*n+1):i*n)';
    end
    tdd(tdd == 1000000) = 0;
    tz = sum(sum(tdd));
    
    tz0 = sum(sum(tc));
    
    %% 数据导出
    % 定义道路等级
    zhugandao = [1,2;2,3;3,4;4,11;11,10;10,9;9,8;8,5]';
    cigandao = [1,5]';
    zhilu = [2,6;6,9;3,7;7,10;6,7]';
    
    
    for i = 1:size(zhugandao,2)
        x= zhugandao(1,i);
        y = zhugandao(2,i);
    
        dis_zhu(1,i) = X(x,y); % 正向交通量
        dis_zhu(2,i) = X(y,x); % 反向交通量
        dis_zhu(3,i) = Cmax(x,y); % 路段容量
        dis_zhu(4,i) = dis_zhu(1,i)/dis_zhu(3,i); % 正向v/c比
        dis_zhu(5,i) = dis_zhu(2,i)/dis_zhu(3,i); % 反向v/c比
        dis_zhu(6,i) = 0;
        dis_zhu(7,i) = tdd(x,y); % 路段时间
        dis_zhu(8,i) = tdd(y,x); % 反向路段时间
        dis_zhu(9,i) = (dis_zhu(7,i)+dis_zhu(8,i))/2; % 平均路段时间
        dis_zhu(10,i) = distance(x,y); % 路段长度
        dis_zhu(11,i) = dis_zhu(9,i)/3600;
        dis_zhu(12,i) = dis_zhu(10,i)/dis_zhu(11,i); % 路段速度
    end
    zhugandao = [zhugandao;dis_zhu];
    zhugandao = sortrows(zhugandao',1);
    zhugandao = zhugandao';
    
    for i = 1:size(cigandao,2)
        x= cigandao(1,i);
        y = cigandao(2,i);
    
        dis_ci(1,i) = X(x,y); % 正向交通量
        dis_ci(2,i) = X(y,x); % 反向交通量
        dis_ci(3,i) = Cmax(x,y); % 路段容量
        dis_ci(4,i) = dis_ci(1,i)/dis_ci(3,i); % 正向负载量
        dis_ci(5,i) = dis_ci(2,i)/dis_ci(3,i); % 反向负载量
        dis_ci(6,i) = 0;
        dis_ci(7,i) = tdd(x,y);
        dis_ci(8,i) = tdd(y,x);
        dis_ci(9,i) = (dis_ci(7,i)+dis_ci(8,i))/2;
        dis_ci(10,i) = distance(x,y);
        dis_ci(11,i) = dis_ci(9,i)/3600;
        dis_ci(12,i) = dis_ci(10,i)/dis_ci(11,i);
    end
    cigandao = [cigandao;dis_ci];
    
    for i = 1:size(zhilu,2)
        x= zhilu(1,i);
        y = zhilu(2,i);
        dis_zhi(1,i) = X(x,y); % 正向交通量
        dis_zhi(2,i) = X(y,x); % 反向交通量
        dis_zhi(3,i) = Cmax(x,y); % 路段容量
        dis_zhi(4,i) = dis_zhi(1,i)/dis_zhi(3,i); % 正向负载量
        dis_zhi(5,i) = dis_zhi(2,i)/dis_zhi(3,i); % 反向负载量
        dis_zhi(6,i) = 0;
        dis_zhi(7,i) = tdd(x,y);
        dis_zhi(8,i) = tdd(y,x);
        dis_zhi(9,i) = (dis_zhi(7,i)+dis_zhi(8,i))/2;
        dis_zhi(10,i) = distance(x,y);
        dis_zhi(11,i) = dis_zhi(9,i)/3600;
        dis_zhi(12,i) = dis_zhi(10,i)/dis_zhi(11,i);
    end
    zhilu = [zhilu;dis_zhi];
    zhilu = sortrows(zhilu',1);
    zhilu = zhilu';
    result = [zhugandao,cigandao,zhilu] % 路段信息按列排布,每行的含义参考上文注释

    结果如下:
    这里写图片描述

    存在的问题

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

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

    展开全文
  • Frank-Wolfe算法是用于求解交通流量分配问题的经典算法,但该算法是基于路段(Link-Based)的交通流量分配算法,无法用于求解路径交通流量。针对此问题,提出一种用于求解路径交通量的改进Frank-Wolfe算法。通过在...
  • 不精确一维搜索wolfe算法

    热门讨论 2010-06-10 08:43:40
    不精确一维搜索wolfe算法。C++环境编译
  • 本文介绍了约束问题中常用的两种算法——乘子法和 Frank-Wolfe 算法, 并通过matlab编程实现了以上两种算法,并对实际问题进行求解。

    前言
    本文介绍了约束最优化问题中常用的两种算法,乘子法和 Frank-Wolfe 算法。
    最后通过matlab编程实现了以上两种算法,并对实际问题进行求解。

    1. 符号定义
      假设约束最优化问题
      minf(x),(xRn)\min f(x), 其中(x \in R^n)
      s.t. gi(x)0,(iI={1,2,,m1})s.t.\ g_i(x) \geq 0, (i \in I=\{1,2,\cdots, m_1\})
        hj(x)=0,(jE={m1+1,m1+2,,m}))h_j(x) = 0, (j \in E = \{m_1+1, m_1+2, \cdots, m\}))
      记函数f(x)的一阶导数为f(x)\nabla f(x),二阶导数为2f(x)\nabla^2 f(x).

    一、乘子法

    1. 简介
      乘子法(multiplier method)是一种约束极小化的算法,通过引入Lagrange函数来求解得到最优解的。

    2. 适用范围:无约束最优化问题

    3. 算法步骤
      步0:取初始点x(0)Rnx^{(0)} \in R^n,初始乘子向量λ(0)\lambda^{(0)}. 给定罚参数序列{μk}\{\mu_k \},精度ϵ>0\epsilon>0. 令k=0.
      步1:下式构造增广Lagrange函数Lμ(x,λ)L_{\mu} (x, \lambda).
      Lμ(x,λ)=f(x)+12μ1iI[min2{μgi(x)λi,0}λi2]jEλihi(x)+12μjEhi2(x)L_{\mu}(x, \lambda) = f(x) + \frac{1}{2}\mu^{-1}\sum_{i \in I}[min^2 \{\mu g_i(x)-\lambda_i, 0\}-\lambda_i^2] - \sum_{j \in E}\lambda_ih_i(x) + \frac{1}{2}\mu\sum_{j \in E} h_i^2(x)
      步2:以x(k1)x^{(k-1)}作为初始点(k=0时,初始点任意),求解无约束问题minLμk(x,λ(k))\min L_{\mu_k}(x, \lambda^{(k)})得解x(k)x^{(k)}.
      步3:若hE(x(k))+min{gI(x(k)),μk1λI(k)}ϵ||h_E(x^{(k)})|| + ||\min \{g_I(x^{(k)}), \mu^{-1}_k\lambda_I^{(k)}\}|| \leq \epsilon则得解x(k)x^{(k)}.
      步4:在下式中取x=x(k),λ=λ(k)x = x^{(k)}, \lambda = \lambda^{(k)}λ(k+1)\lambda^{(k+1)},令k=k+1, 转步1.
      λj+={λjμhj(x),(jE)max{λiμigi(x),0},(iI)\lambda_j^+ = \left\{ \begin{array}{ll} %该矩阵一共3列,每一列都居中放置 \lambda_j - \mu h_j(x), & (j \in E)\\ \max \{\lambda_i - \mu_ig_i(x), 0\}, & (i \in I)\\ \end{array} \right.
      其中,x,λ,μx, \lambda, \mu表示当前迭代点的值,λj+\lambda^+_j表示下一次迭代的乘子向量.

    4. matlab实现
      说明:函数中使用的initpt函数(用于初始化x0x_0)是为了自己编写的辅助函数,这几个函数的具体实现可见最下方的链接处。

    function [iter_num, time, epsilonk] = method_of_multipliers(nprob, n, x, epsilon, max_iter)
    % 乘子法   2020.06.25
    % @author: 豆奶
    % 函数功能:实现乘子法     % 比较两个字符串 strcmp(a, 'De_Jong')
    % 输入:
    % nprob: 问题编号, 
    % n: 维度
    % x: 自变量x
    % epsilon: int, 
    % max_iter: 最大迭代次数
    % 输出:
    % iter_num: 迭代次数
    % time: 计算时间
    % epsilonk: KKT系统模的最终值
    
    t1=clock;  % 计时开始
    % 步0:取初始点x,初始乘子向量. 给定罚参数序列,精度
    if nargin==2
        x = initpt(n, nprob);
        epsilon = 1e-5;
        max_iter = 100;  
    end
    mu = 0.5;
    xk = x;
    % 设置目标函数,约束函数
    if strcmp(nprob, 'De_Jong')
        f = @(x) sum(x.^2);
        g1 =@(x) x+5.12;
        g2 =@(x) 5.12-x;
        g =@(x) [g1(x); g2(x)];
        h =@(x) 0;
        lambda_g = 10*ones(2*n, 1);  % (2n, 1)   和约束条件的个数有关
        lambda_h = 0;
    elseif strcmp(nprob, 'Rastrigin')
        f = @(x)10*n + sum(x.^2 - 10*cos(2*pi*x));  
        g1 =@(x) x+5.12;
        g2 =@(x) 5.12-x;
        g =@(x) [g1(x); g2(x)];
        h =@(x) 0;
        lambda_g = ones(2*n, 1);  % (2n, 1)   和约束条件的个数有关
        lambda_h = 0;
    end
    for k=1:max_iter 
        % 步1:由(12.27)构造增广Lagrange函数L.
        L =@(x) f(x) + 0.5.*(1/mu)*sum(min(mu*g(x)-lambda_g, 0).^2-lambda_g.^2) - sum(lambda_h'*h(x)) + 0.5*mu*sum(h(x).^2);
        % 步2:以x_{k-1}作为初始点(k=0时,初始点任意),求解无约束问题L得解x_k.
        xk = fmincon(L, xk);
        % 步3:判断是否小于epsilon,以跳出循环.
        eps_k = norm(h(xk)) + norm(min(g(xk), 1/mu*lambda_g));
        if eps_k<=epsilon
            break;
        end 
        % 步4:在(12.28)中取x=xk, lambda=lambda_k
        lambda_g = max(lambda_g-mu*g(xk), 0);  % update lambda_g, i \in I
        lambda_h = lambda_h - mu*h(xk);        % update lambda_h, j \in E
    end
    % 获取输出值
    iter_num = k;
    syms x;                     %通过符号变量将匿名函数转换为符号函数
    gxk = mu*g(xk) - lambda_g;  % 由于min无法直接求导,因此将这部分割裂开来,独自求导计算
    %dg = diff(g(x));  % 由于本次测试采用的两个测试问题的g的导数均为实数,即(1,-1),失去了自变量x,因此改为手动求导
    dg = ones(2*n, 1);
    dg(n+1:2*n) = -1;
    gxk(gxk>0) = 0;
    dgxk = gxk'*dg;
    y = f(x) - sum(lambda_h'*h(x)) + 0.5*mu*sum(h(x).^2);   % y=L(x);
    grad=matlabFunction(diff(y));   %通过matlabFunction将符号函数转换为匿名函数,并求导
    epsilonk = [grad(xk)+dgxk; h(xk); min(g(xk),lambda_g)];   
    epsilonk = norm(epsilonk);      
    t2=clock;
    time = etime(t2,t1);
    fprintf('\t问题编号\t\t\t算法\t\t运行次数\t\t运行时间\t\t\tKKT残量\n');
    fprintf('\t%s\t\t乘子法\t%5d\t\t%2f\t\t%5d\n', nprob, iter_num, time, epsilonk);
    return;
    end
    

    二、 Frank-Wolfe 算法

    1. 简介
      解约束问题可行方向法与求解无约束问题的下降算法类似,但对约束问题而言,我们关心的是问题的可行点。可行方向法从问题的可行点出发,在该点的可行方向中,寻找使目标函数下降的方向。然后沿该方向进行线性搜索,得到一个可行点。如此得到一个点列,并希望该点列能终止于问题的解,或其极限点是问题的解。
      Frank-Wolfe 算法便是利用求解约束问题可行方向法的思想,该算法通过解下面线性规划问题来确定可行下降方向。 其搜索方向d(k)d^{(k)}可理解为可行点x(k)x^{(k)}处的可行方向中的最速下降法,是最速下降法的一种推广。
      minyDLf(x(k))+f(x(k))T(yx(k)) \min_{y \in D_L} f(x^{(k)}) + \nabla f(x^{(k)})^T(y - x^{(k)})
      或等价的
      minyDLf(x(k))T(yx(k)) \min_{y \in D_L} \nabla f(x^{(k)})^T(y - x^{(k)})
      其中DLD_L是问题的可行域.

    2. 适用范围:适用于求解约束函数为线性函数的最优化问题

    3. 算法步骤
      步0:取初始点x(0)DLx^{(0)} \in D_L, 精度ϵ>0\epsilon > 0. 令k=0.
      步1: 求解下面线性规划问题得y(k)y^{(k)}. 令d(k)=y(k)x(k)d^{(k)} = y^{(k)} - x^{(k)}.
      minyDLf(x(k))T(yx(k)) \min_{y \in D_L} \nabla f(x^{(k)})^T(y - x^{(k)})
      步2:若f(x(k))Td(k)ϵ|\nabla f(x^{(k)})^Td^{(k)}| \leq \epsilon,则得解x(k)x^{(k)}. 否则,转步3.
      步3:由线性搜索
      min0t1f(x(k)+td(k))=f(x(k)+tkd(k))\min_{0\leq t \leq 1} f(x^{(k)} + td^{(k)}) = f(x^{(k)} + t_kd^{(k)})确定步长tkt_k. 令x(k+1)=x(k)+tkd(k),k=k+1x^{(k+1)} = x^{(k)} + t_kd^{(k)}, k=k+1. 转步1.

    4. matlab实现
      说明:函数中使用的initpt函数(用于初始化x0x_0)是为了自己编写的辅助函数,这几个函数的具体实现可见最下方的链接处。

    function [iter_num, time, epsilonk] = Frank_Wolfe(nprob, n, x, epsilon, max_iter)
    % Frank Wolfe算法   2020.06.25
    % @author: 豆奶
    % 函数功能:实现Frank Wolfe法     
    % 输入:
    % nprob: 问题编号, 
    % n: 维度
    % x: 自变量x
    % epsilon: int, 
    % max_iter: 最大迭代次数
    % 输出:
    % iter_num: 迭代次数
    % time: 计算时间
    % epsilonk: KKT系统模的最终值
    
    t1=clock;   % 计时开始
    if nargin==2
        % 步0:取初始点x和精度epsilon
        x = initpt(n, nprob);  % 选取初始点要注意在可行域内 
        epsilon = 0.0001;
        max_iter = 100;
    end
    xk = x;
    % 设置目标函数,约束函数,这一部分f和df也可以改成参数的形式传入.
    if strcmp(nprob, 'De_Jong')
        f = @(x) sum(x.^2);
        df = @(x) 2.*x;
        g1 =@(x) x+5.12;
        g2 =@(x) 5.12-x;
        g =@(x) [g1(x); g2(x)];
        h =@(x) 0;
    elseif strcmp(nprob, 'Rastrigin')
        f = @(x)10*n + sum(x.^2 - 10*cos(2*pi*x));   % n为常量
        df = @(x) 2*x + 20*pi*sin(2*pi*x);   % f 的导数
        g1 =@(x) 5.12+x;
        g2 =@(x) 5.12-x;
        g =@(x) [g1(x); g2(x)];
        h =@(x) 0;
    end
    syms t;
    for k=0:max_iter
        dfx = df(xk);
        lb = -5.12*ones(n, 1);
        ub = 5.12*ones(n, 1);
        % 步1:利用matlab的linprog函数求解线性规划问题(13.8)得y
        % 参数说明: linprog(f,A,b,Aeq,beq,LB,UB, X0)  
        y = linprog(dfx,[],[],[],[],lb,ub);% 这里要改(6.22), 改好了(6.22)
        d = y - xk;
        % 步2:判断是否满足终止条件
        if norm(dfx'*d)/(1+norm(f(xk)))<epsilon
            break;
        end
        % 步3:利用matlab的fmincon函数求解线性搜索确定步长t
        ft = matlabFunction(f(xk+t*d));  %通过matlabFunction将符号函数转换为匿名函数,因为fmincon的第一个输入是匿名函数
        % 参数说明: fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) %最后两个输入这里不需要
        [t_ans] = fmincon(ft,0.5,[],[],[],[],0,1);  % 这里要改(6.21)  改好了(6.22)
        xk = xk + t_ans*d;   
    end
    % 获取输出值
    iter_num = k;
    lambda_g = zeros(2*n, 1);  % 由于Frank_Wolfe算法没有用到lagrange函数,故lambda设为0
    epsilonk = [df(xk); h(xk); min(g(xk), lambda_g)]; 
    epsilonk = norm(epsilonk); 
    t2=clock;
    time = etime(t2,t1);
    fprintf('\t问题名称\t\t\t算法\t\t运行次数\t\t运行时间\t\t\tKKT系统残量\n');
    fprintf('\t%s\t\tFrank_Wolfe\t%5d\t\t%2f\t\t%5d\n', nprob, iter_num, time, epsilonk);
    return;
    end
    

    参考文献:

    • 数值最优化算法与理论(第二版)

    完整代码与函数问题测试github(点我)

    如有不对之处,还望指出,我会进行修改的,谢谢!

    展开全文
  • 针对用户均衡交通分配问题,提出一种可以避免穷举网络中的所有路径的基于Frank Wolfe算法的路径交通量求解方法。它在已知一组满足用户均衡规则的基于终点的路段交通量和交通网络中各个 OD(origin destination)对间的...
  • 用于快速且可扩展的套索回归的随机Frank-Wolfe算法(C ++代码) 免责声明: 仅在GNU通用公共许可证下,该软件可用于非商业研究目的。 然而,尽管GNU通用公共许可证的任何规定,该软件可能无法接触后用于商业用途的...
  • TAC接受的论文“用于凸和非凸问题的分散式Frank-Wolfe算法”的同伴代码将于2017年发表。[预印本:] 内容 矩阵完成示例 DeFW_MC_SyntheticData.m-图1-4的仿真代码 DeFW_MC_MLData.m-图5的仿真代码 LASSO范例 DeFW_...
  • Deep Frank-Wolfe用于神经网络优化 该存储库包含pytorch中的论文的实现》。 如果您将这项工作用于...注意:您可能对我们的后续算法感兴趣,该算法具有明显的收敛保证,并且在我们的实验中优于DFW。 要求 此代码适用于
  • 在对偶理论作用下,将约束正项几何...利用Frank-wolte算法以及几何规划和约束条件的特点,为有多个变量的几何规划构造出了一种有效的间接算法,而且此方法更适用于困难度大于零的几何规划问题,实验表明此方法是可行的。
  • 针对在飞机路线和机组路线恢复方案中忽略了旅客行程衔接的问题,在飞机路线恢复方案的基础上,通过构建时空网络,建立了不正常航班旅客流恢复的多商品流模型,设计了Danzig-Wolfe分解算法进行求解,给出旅客流恢复方案。...
  • 配流07—基于BPR函数的Frank Wolfe算法

    千次阅读 热门讨论 2016-12-21 16:07:27
    disp(' 《基于BPR函数的Frank Wolfe算法》'); disp('运行环境:MATLAB 8.3.0.532 '); disp('制 作 人:兰州交通大学 刘志祥'); disp('Q Q:531548824'); disp('====================================================...

    一、问题描述

    在道路网中,已知OD需求,路段走行时间,路段能力和路径路段关系,求流量的均衡分配结果。


    二、算法描述

    此处只给出大的步骤,精确的算法描述见第三节——算法程序。

    step1:给定路网数据,OD需求,路段能力

    step2:自动求出路径和路段数量,根据路段数量定义路段名,给定初始数据

    step3:构造阻抗函数并求出初始阻抗,此处用BPR函数

    step4:全有全无配流

    step5:数据更新

    step6:求目标函数值

    step7:输出计算结果


    三、算法程序


    clear
    clc
    disp('========================================================================');
    disp('                      《基于BPR函数的Frank Wolfe算法》');
    disp('运行环境:MATLAB 8.3.0.532 ');
    disp('制 作 人:兰州交通大学   刘志祥');
    disp('Q      Q:531548824');
    disp('=========================================================================');
    
    %% 1 给定路网数据,OD需求,路段能力
    %算例1
    %==========================================================================
    % Q=1000;
    % W=[5 6 6 8 3];
    % Cmax=[600 500 600 500 700];
    % Mxf=[1 0 0 1 0;1 0 1 0 1;0 1 0 0 1];
    %==========================================================================
    
    %算例2
    %==========================================================================
    Q=3000;                                 %OD需求
    W=[5 9 6 7 1];                          %路段初始阻抗
    Cmax=700*ones(1,5);                     %路段能力
    Mxf =[1 0 0 1 0;0 1 0 0 1;1 0 1 0 1];   %路径路段0-1关系
    %==========================================================================
    
    %% 2 自动求出路径和路段数量,根据路段数量定义路段名,给定初始数据
    %==========================================================================
    numf=size(Mxf,1);
    numx=size(Mxf,2);
    syms lambda real
    for i=1:numx
        syms x(i) real;
    end
    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);
    Ckrs=(Mxf*t')';                        %路径的走行时间初值
    %==========================================================================
    
    %% 4 全有全无配流
    %==========================================================================
    [Min,index]=min(Ckrs);
    X1=Mxf(index,:).*Q;                    %全有全无法为最短路径上的路段分配流量
    %==========================================================================
    
    %% 5 数据更新
    %==========================================================================
    while e>1e-3                           %精度判断
        cont=cont+1;                       %迭代次数更新
        t=(W).*(1+0.15*(X1./Cmax).^4);     %路段时间跟新
        Ckrs=(Mxf*t')';                    %路径时间更新
        [Min,index]=min(Ckrs);
        Y1=Mxf(index,:).*Q;                %全有全无法求辅助流量
        S=Y1-X1;                           %搜索方向
        X2=X1+lambda*S;                    %先将X2用X1和lambda进行表达
        t=(W).*(1+0.15*(X2./Cmax).^4);     %含lambda的阻抗表达
        f=sum(S.*t,2);                     %2表示按行求和
        lambda1=double(solve(f));          %求解方程,确定步长。
        k=length(lambda1);                 %如步长lambda1的解不唯一,取实数,且大于0小于1;
        for m=1:k
            if lambda1(m,1)>=0&&lambda1(m,1)<=1
                lambda2=lambda1(m,1);
            end
        end
        X2=X1+lambda2*S;                   %得到下一步的流量值,且进行下一次迭代
        e=sqrt(sum((X2-X1).^2))/sum(X1);   %精度计算
        X1=X2;                             %流量更新
    end
    %==========================================================================
    
    %% 6 求目标函数值
    %==========================================================================
    Xx=zeros(numx,1);                      %积分下界
    Xn=X1;                                 %积分上界
    Zf=zeros(numx,1);                      %目标值元素初始化
    for i=1:numx
        Zf(i)=int(tt(i),Xx(i),Xn(i));      %对每一个路径积分
    end
    Z=sum(Zf);                             %总目标=各路径阻抗求和
    %==========================================================================
    
    %% 7 输出计算结果
    %==========================================================================
    disp('*************************************************************************')
    disp(['     迭代次数:',num2str(cont)]);
    disp(['     误 差 值:',num2str(e)]);
    disp(['     配流结果:',num2str(Xn)]);
    disp(['     路径阻抗:',num2str(Ckrs)]);
    disp(['     目 标 值:',num2str(Z)]);
    disp('*************************************************************************')
    %==========================================================================

    四、算例及运行结果

    1.算例

    路网如图1,已知起点v1,终点v4,路段编号如图1所示,路段初始阻抗为W=[5 9 6 7 1],对应的路段能力为Cmax=[500 700 800 600 700],求流量分配结果。


    2.运行结果

    >> Frank_wolfe_BPR_jingdian

    ========================================================================
                          《基于BPR函数的Frank Wolfe算法》
    运行环境:MATLAB 8.3.0.532 
    制 作 人:兰州交通大学   刘志祥
    Q      Q:531548824
    =========================================================================
    *************************************************************************
         迭代次数:13
         误 差 值:0.00083827
         配流结果:1347.6659      1652.3341      304.76636      1042.8996      1957.1004
         路径阻抗:61.5966      60.7377      61.0709
         目 标 值:62802.3986
    *************************************************************************
    3.数据整理

    路段流量表
    路段 a1 a2 a3 a4 a5
    流量 1347.7 1652.3 304.8 1042.9 1957.1


    五、结论及展望

    基于BPR函数的FW算法是科学有效的配流算法,可以看到最终3条路径的阻抗基本相等,这符合均衡配流的预期结果。


    六、代码下载

    代码下载地址:待更新...



    展开全文
  • 有关交通分配的关键分配算法, 值得一看
  • wolfe 算法_最优化课课后作业笔记

    千次阅读 2016-09-11 10:48:37
    from numpy import * import sys MAXINT=sys.maxint def gradient_(X): x1=X[0];x2=X[1] return array([(-400*(x2-x1**2)*x1-2*(1-x1)),200*(x2-x1**2)]) def f(X): x1=X[0];x2=X[1] return 100*(x2-
    from numpy import *
    import sys
    MAXINT=sys.maxint
    def gradient_(X):
        x1=X[0];x2=X[1]
        return array([(-400*(x2-x1**2)*x1-2*(1-x1)),200*(x2-x1**2)])
    def f(X):
        x1=X[0];x2=X[1]
        return 100*(x2-x1**2)**2+(1-x1)**2
    if "__main__":
        '''
        Xk:point
        Pk:direction
        '''
        Xk=array([-1,1])
        Pk=array([1,1])
    
        mu_=0.1
        sigma=0.5
        Fk=f(Xk)
        Gk=gradient_(Xk)
        print 'Gk=%s ^T,Fk=%s'%(Gk,Fk)
        print sum(gradient_(Xk)*Pk)
    
        alpha=1;b=MAXINT;a=0;j=0
        flag=1
        i=0
        while(flag):
            Xk1=Xk+alpha*Pk
            Fk1=f(Xk1)
            print i,alpha,Xk1,Fk1
            if Fk-Fk1<-(mu_*alpha*sum(gradient_(Xk)*Pk)):
                b=alpha
                alpha=(alpha+a)*0.5
                i=i+1
            elif sum(gradient_(Xk1)*Pk)<sigma*sum(gradient_(Xk)*Pk):
                a=alpha
                alpha=min(2*alpha,(a+b)*0.5)
                i=1+1
            else:
                flag=0
        Xk1 = Xk + alpha * Pk
        Fk1 = f(Xk1)
        print alpha,Fk1,Xk1
    
    
    展开全文
  • F r a nk - Wo lfe ( F W) 算法是一类广泛应用于求解交通分配问题的算法. 它具有容 易编程实现, 所需内存少的特点. 但是该算法收敛速度较慢, 不能得到路径信息. 为了 提高算法的效率, 本文研究三种流量更新策略( al ...
  • 最优化算法wolfe搜索

    2017-06-04 15:21:58
    本程序为matlab编写的最优化一维搜索方法wolfe
  • Frank-Wolfe算法的步骤总结如下: 选择初值。初始点x0∈Sx0∈S,给定允许误差ϵ>0ϵ>0,令k=0k=0。 求解近似线性规划: min∇f(xk)Txs.t.x∈S min∇f(xk)Txs.t.x∈S 得最优解ykyk。 构造可行下降方向。令dk=yk−...
  • 配流06—frank_wolfe配流算法

    千次阅读 2015-10-17 02:21:15
    %说明:本程序为frankwolfe算法,通过输入总流量、在路网情况已知时配流并求出费用值。
  • Frank-Wolfe近邻牛顿算法 介绍 该算法可以解决以下约束凸优化问题: 其中,是自一致的,是紧的凸集,其线性优化预言很容易找到。 先决条件 该代码在Matlab R2018b下进行了测试,不需要其他MATLAB工具箱。 在运行示例...
  • 将结构正割法应用到拟牛顿算法中,利用目标函数的梯度信息和函数值信息,引入拟牛顿方程,采用Wolfe线搜索准则,给出了求解无约束优化问题的一个新算法,并在一定条件下证明了新算法的收敛性和超线性收敛性。
  • 提出了一种记忆梯度法的主要参数以的新形式,分析了该算法Wolfe-Power搜索下的全局收敛性,适合解决大型优化问题。
  • 《DFP算法+wolfe性非线性搜索解决无约束问题的matlab程序》由会员分享,可在线阅读,更多相关《DFP算法+wolfe性非线性搜索解决无约束问题的matlab程序(2页珍藏版)》请在人人文库网上搜索。1、function x,val,k=dfp...
  • 本文提出一类广义Wolfe线搜索模型,并且把它与著名的BFGS方法相结合,对于所得到的算法证明了:对于凸函数算法具有全局收敛性和超线性收敛速度。这推广了参考文献的结果。
  • 翻译自numerical optimization Chapter3 Line Search MethodsWolfe条件(或强Wolfe条件)是广泛使用并且较为有效的终止条件。下面描述一个一维的线搜索过程,该过程确保对于给定的任意参数...算法分为两个步骤,第一步

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 854
精华内容 341
关键字:

wolfe算法