拉格朗日松弛基本概念及在约束最短路径问题中的应用

拉格朗日松弛简介

拉格朗日松弛（Lagrangian Relaxation）是一种约束优化问题里处理约束的思想。其将约束分为简单约束和困难约束，通过一个拉格朗日乘子将困难约束罚至目标函数上。拉格朗日松弛是1971年Held和Krap在研究旅行商问题（travelling salesman problem）时提出的。

假设我们考虑问题(不一定是凸问题)

[egin{array}{l l} minlimits_{m{x} }& f(m{x}) \ ext{s.t.} & g_i(m{x}) leq 0;i = 1,cdots,m\ & m{x} in X\ end{array} ]

其中(X)是(mathbb{R}^n)中的一个有限集。有时候，(g_i(m{x}) leq 0;i = 1,cdots,m)使求解非常困难，称为困难约束，而(m{x} in X)是容易求解的简单约束。拉格朗日松弛就是利用引入的拉格朗日乘子(lambda_i ge 0)将困难约束作为惩罚项写到目标函数里:

[egin{array}{l l} minlimits_{m{x},lambda}& f(m{x})+Sigma_{j=1}^m lambda_i g_i(m{x}) \ ext{s.t.} & m{x} in X\ end{array} ]

并将其称作拉格朗日函数，即(L(lambda) = minlimits_{m{x}} {f(m{x})+Sigma_{j=1}^m lambda_ig_i(m{x})|m{x} in X})。注意拉格朗日函数是关于(lambda)的分片线性函数，因为它是一系列仿射函数的逐点最小函数，所以是凹函数。我们称

[egin{array}{l l} maxlimits_{lambda} & L(lambda)\ end{array} ]

为拉格朗日对偶(Lagrangian dual)，有些地方也称为拉格朗日乘子问题(Lagrangian multiplier problem)。困难约束可以是等式约束也可以是不等式约束。上述做法降低了求解的复杂程度，并使我们得到一个原问题最优解的下界。有时这个下界可以非常接近甚至等于最优解。

写在中间

拉格朗日松弛（Lagrangian Relaxation）是我们凸优化课程读书报告的主题，读书报告全文可见 LagrangianRelaxation.pdf。

读书报告的内容包括：

• 拉格朗日松弛的概念
  • 困难约束为线性等式时
  • 困难约束为线性不等式时
  • 求解拉格朗日对偶：次梯度法（Subgradient）
• 例子：使用拉格朗日松弛求解约束最短路径问题
  • 约束最短路径问题
  • 一个具体的例子
    • (L(lambda))中(lambda)取值的影响
    • 求解(L^*)

报告里的内容就不在博客里再写一遍了（太麻烦了），大家可以直接去看全文 LagrangianRelaxation.pdf。这里只说一下关于程序实现的部分。

使用拉格朗日松弛求解约束优化问题的程序

强烈建议大家通读过全文之后再来看程序。

编程使用Matlab。这里求解凸优化问题使用了CVX工具包，安装和教程可见：CVX: Matlab Software for Disciplined Convex Programming。

程序1：指定乘子(lambda)，求解松弛后的问题

% Constrainted Shortest Path
% 在指定拉格朗日乘子lambda的前提下，
% 使用Lagrangian Relaxation 对 Constrainted Shortest Path 问题求解

clc
clear all

nInf = 10000;
% parametters
Cost = [nInf, 1, 10, nInf, nInf, nInf;
        nInf, nInf, nInf, 1, 2, nInf;
        nInf, 1, nInf, 5, 12, nInf;
        nInf, nInf, nInf, nInf, 10, 1;
        nInf, nInf, nInf, nInf, nInf, 2;
        nInf, nInf, nInf, nInf, nInf, nInf];
    
Time = [0, 10, 3, nInf, nInf, nInf;
              nInf, 0, nInf, 1, 3, nInf;
              nInf, 2, 0, 7, 3, nInf;
              nInf, nInf, nInf, 0,1, 7;
              nInf, nInf, nInf, nInf, 0, 2;
              nInf, nInf, nInf, nInf, nInf, 0];
          
n = 6; % number of nodes
T = 10; % max total time

lambda = 1; % lagrangian multiplier

cvx_begin % quiet
    variable X(n, n)  
    minimize( sum(sum( Cost .* X) )  +   lambda * (sum(sum(Time .* X)) - T) )
    subject to

        sum(X(1, :)) - sum(X(:, 1)) == 1;
        for i = 2 : n-1
            sum(X(i, :)) - sum(X(:, i)) == 0;
        end
        sum(X(n, :)) - sum(X(:, n)) == -1;

        0 <= X <= 1;  % linear relaxation for the integer constraints
cvx_end

程序2：使用次梯度法找到最佳的(lambda)值（求解对偶问题）

注意cvx_solver Mosek。为了获得更高的求解精度，这里使用了一个求解优化问题的商业求解器：Mosek。其和CVX配合使用的方法可见 Using MOSEK with CVX。

% Constrainted Shortest Path
% 使用 subgradient 方法找到最佳的 lambda 值

clc
clear all

nInf = 10000;

Cost = [nInf, 1, 10, nInf, nInf, nInf;
        nInf, nInf, nInf, 1, 2, nInf;
        nInf, 1, nInf, 5, 12, nInf;
        nInf, nInf, nInf, nInf, 10, 1;
        nInf, nInf, nInf, nInf, nInf, 2;
        nInf, nInf, nInf, nInf, nInf, nInf];
    
Time = [0, 10, 3, nInf, nInf, nInf;
              nInf, 0, nInf, 1, 3, nInf;
              nInf, 2, 0, 7, 3, nInf;
              nInf, nInf, nInf, 0,1, 7;
              nInf, nInf, nInf, nInf, 0, 2;
              nInf, nInf, nInf, nInf, nInf, 0];   
          
n = 6;
T = 10;

% formulate the cost constraints as Ax = b
A = zeros(n, n * n);
for i = 1 : n
    for j = 0 : n - 1
        if j == i - 1
            A(i, j*n + 1 : j*n + n) = - ones(1, n);
            A(i, j*n + i) = 0;
        else
            A(i, j * n + i) = 1;
        end
    end
end

b = zeros(n, 1);
b(1, 1) = 1;
b(n, 1) = -1;


% subgradient iteration
k = 0;
lambda = 1;
while (1)
    [x, L] = cvx_optimize(Cost, Time, n, T, A, b, lambda);
    
    g = sum(sum(Time(:) .* x)) - T;
    
    fprintf('k: %i, lambda: %f, L: %f, g: %f 
', k, lambda, L, g);
    
    if g <= 0 && lambda * g == 0
        break;
    end
    
    theta =  0.1 * (20 - L) / (norm(g) ^ 2);
    
    lambda = max(0, lambda + theta' * g);
    
    k = k + 1;
    
end


function [x, fval] = cvx_optimize(Cost, Time, n, T, A, b, lambda)
    cvx_solver Mosek
    cvx_begin quiet
        variable x(n * n)
        minimize( sum(sum( Cost(:) .* x) )  +   lambda * (sum(sum(Time(:) .* x)) - T) )
        subject to
            A * x <= b;
            0 <= x <= 1
    cvx_end
    
    fval = sum(sum( Cost(:) .* x) )  +   lambda * (sum(sum(Time(:) .* x)) - T);
end