对偶上升实例-MATLAB代码

正仪

已于 2023-10-21 21:26:50 修改

阅读量3.3k

点赞数 11

分类专栏：基础数学类文章标签： matlab 线性代数矩阵

于 2022-01-28 22:15:54 首次发布

本文链接：https://blog.csdn.net/weixin_43509834/article/details/122737009

版权

基础数学类专栏收录该内容

40 篇文章 10 订阅

订阅专栏

一、本文概述：

本文给出对偶上升法（dual ascent）求解凸优化问题最优解的代码实例。如果您觉得对您有帮助，请点个赞，加个收藏，谢谢！

二、简单实例

本文以下述实例为例，撰写对偶上升法的迭代步骤，并给出最终可运行的MATLAB代码，以便大家上手学习。
1）优化问题为：
$\begin{aligned} &\mathop{min}\limits_{x_1,x_2} \ f(x_1,x_2)=2(x_1-1)^2+(x_2+2)^2\\ &s.t. \begin{cases} x_1\geq2 \\ x_2\geq0 \end{cases} \end{aligned}$

2）拉格朗日函数为：
$\begin{aligned} &L(x_1,x_2,\lambda_1,\lambda_2) \\ &=2(x_1-1)^2+(x_2+2)^2+\lambda_1(2-x_1)+\lambda_2(-x_2)\\ \end{aligned}$

3）决策变量 $x_1,x_2$ 解析解为：
$\begin{aligned} \frac{\partial L(x_1,x_2,\lambda_1,\lambda_2) }{\partial x_1}=0 \\ \frac{\partial L(x_1,x_2,\lambda_1,\lambda_2) }{\partial x_2}=0 \end{aligned}$
求解上述方程，解得：
$\begin{aligned} x_1^*=\frac{\lambda_1}{4}+1 \\ x_2^*=\frac{\lambda_2}{2}-2 \end{aligned}$

4）拉格朗日函数更新公式：
第k+1轮迭代更新公式：
$\begin{aligned} \lambda_1^{k+1}&=\lambda_1^{k}+t_1(2-x_1) \\ \lambda_2^{k+1}&=\lambda_2^{k}+t_1(-x_2) \end{aligned}$

至此，对偶上升方法完毕求解完毕，下面为MATLAB代码！

三、MATLAB代码

clear all;
clc;

x = [100 -120]';   %决策变量初始化   %注：本函数的决策变量为向量！ 

size_x_row  = size(x,1);
size_x_line = size(x,2);
lambda = zeros( size_x_row , size_x_line );   %拉格朗日乘子初始化 
k=1;
G_v = [];  lambda_v = [];  x_v = [];  L_v = [];  f_v = [];
tol = 10^(-5);
t_lambda=0.1;  % 拉格朗日乘子的更新步长
t_x=0.1;      % 决策变量 x 的更新步长
max_iteration = 10000;

for i = 1 : max_iteration
    %% 迭代更新决策变量 x
    x(1) = x(1) - t_x * (  4 * ( x(1) - 1 ) - lambda(1)  );  %% 注意这里可以写解析解，也可以写梯度迭代更新解，也可混合写！
    x(2) = x(2) - t_x * (  2 * ( x(2) + 2 ) - lambda(2)  ); 
    %x(1) = lambda(1)/4+1; 
    %x(2) = lambda(2)/2-2; 
    
    %% 迭代更新拉格朗日乘子 lambda 
    Grad_f = [ (2-x(1)) ; (-x(2)) ] ; 
    lambda = max(   lambda + t_lambda * Grad_f  ,  zeros( size_x_row , size_x_line )    );  
    f =  2*( x(1)-1 )^2 + ( x(2)+2 )^2 ; 
    L =  2*( x(1)-1 )^2 + ( x(2)+2 )^2 + lambda(1)*(2-x(1)) + lambda(2)*(-x(2));
    
    %% 记录历史更新结果
    G_v = [G_v Grad_f];
    x_v = [x_v x];
    lambda_v = [lambda_v lambda];
    f_v = [f_v f]; 
    L_v = [L_v L];
    
    %% 终止条件判断
    % if (k>2) &&( abs(L - L_v(k-1)) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
    if (k>2) &&( abs( f - L ) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
        break;
    end
    k=k+1;
end
if i == max_iteration
    disp(['超过最大迭代次数！']);
end
disp(['迭代次数为：',num2str(k),'次']);
disp(['最优值点为：']);
x 
figure(1)
plot( 2:length(L_v),L_v(2:length(L_v)), '-ro', 'linewidth',2)
disp(['最优值为：',num2str(L),'']);

figure(2)
plot( 2:length(L_v), f_v(2:length(f_v)) - L_v(2:length(L_v)) , '-b*', 'linewidth',2)

四、注意事项及扩展

注意：上述第3)步中，若无法得到解析解，则求梯度亦可。换言之，上述第3)步可用下述代码替代：（计算方法为：拉格朗日函数分别对 $x_1$ 与 $x_2$ 求偏导即可）
$\begin{aligned} x_1^{k+1}&=x_1^{k} - t_2[4(x_1-1)-\lambda_1] \\ x_2^{k+1}&=x_2^{k} - t_2[2(x_2+2)-\lambda_2] \end{aligned}$
MATLAB代码中将解析解的情况注释掉了，用数值解求解的。换言之代码的“第18行第19行”与代码中的“第20行第21行”是等价替换的，运行哪一组都可以！！！当然，也可以混合运行！！！

五、运行结果

六、基础回顾

这里需要注意，不要把梯度下降前的更新正负号写反了，对拉格朗日乘子更新时，其梯度前是加号（即：加梯度），而对决策变量更新时，其梯度前为减号（即：减梯度）。其原因是，前者为梯度上升，后者为梯度下降。我们回顾一下，拉更朗日函数的目标是求解下述优化问题：

$\begin{aligned} \mathop{max}\limits_{\lambda_1,\lambda_2} \left[ \mathop{min}\limits_{x_1,x_2} \ L(x_1,x_2,\lambda_1,\lambda_2)\right] \end{aligned}$
因此，选取拉格朗日乘子 $\lambda_1,\lambda_2$ 以最大化 $L(x_1,x_2,\lambda_1,\lambda_2)$ ，因此对拉格朗日乘子做梯度上升，所以是加号，而选取决策变量 $x_1,x_2$ 以最小化 $L(x_1,x_2,\lambda_1,\lambda_2)$ ，因此对拉格朗日乘子做梯度下降，所以是减号。

另：若有时间，会陆续更新 “精确直线搜索代码” 及 “ADMM的consensus代码”

七、后期更正

终止条件应是原函数f 与拉格朗日函数L 只差足够小。而并非拉格朗日函数L 自身收敛。

八、变换函数

甚至更换函数，并在0附近，该方法依然可行

$\begin{aligned} &\mathop{min}\limits_{x_1,x_2} \ f(x_1,x_2) = \frac{1}{x_1}+x_1^2 + \frac{1}{x_2}+x_2^2\\ &s.t. \begin{cases} 0<x_1\leq1 \\ 0<x_2\leq1 \end{cases} \end{aligned}$

clear all;
clc;

x = [-10 10]';   %决策变量初始化   %注：本函数的决策变量为向量！ 

lambda = ones( 4 , 1 );   %拉格朗日乘子初始化 
k=1;
G_v = [];  lambda_v = [];  x_v = [];  L_v = [];  f_v = [];
tol = 10^(-5);
t_lambda=0.01;  % 拉格朗日乘子的更新步长
epsilon = 0.1;
t_x_1=0.1;      % 决策变量 x 的更新步长
t_x_2=0.1;      % 决策变量 x 的更新步长
max_iteration = 10000;

for i = 1 : max_iteration
    %% 迭代更新决策变量 x
    x(1) = x(1) - t_x_1 * (  - x(1)^(-2) + 2 * x(1) - lambda(1) + lambda(2)   );  %% 注意这里可以写解析解，也可以写梯度迭代更新解，也可混合写！
    x(2) = x(2) - t_x_2 * (  - x(2)^(-2) + 2 * x(2) - lambda(3) + lambda(4) ); 
        
    %% 迭代更新拉格朗日乘子 lambda 
    Grad_f = [ ( -x(1) ) ; ( x(1)-1 ) ;  ( -x(2) ) ; ( x(2)-1 )  ] ; 
    lambda = max(   lambda + t_lambda * Grad_f  ,  zeros( 4 , 1 )    );  
    f =  ( x(1)^(-1) ) + ( x(1)^(2) ) + ( x(2)^(-1) ) + ( x(2)^(2) ) ; 
    L =  ( x(1)^(-1) ) + ( x(1)^(2) ) + ( x(2)^(-1) ) + ( x(2)^(2) ) + lambda(1)*( -x(1) ) + lambda(2)*( x(1)-1 ) + lambda(3)*( -x(2) ) + lambda(4)*( x(2)-1 ) ;
    
    %% 记录历史更新结果
    G_v = [G_v Grad_f];
    x_v = [x_v x];
    lambda_v = [lambda_v lambda];
    f_v = [f_v f]; 
    L_v = [L_v L];
    
    %% 终止条件判断
    % if (k>2) &&( abs(L - L_v(k-1)) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
    if ( k > 2 ) && ( abs( (f-L)  ) < tol ) &&  ( Grad_f(1)<=epsilon )  &&  ( Grad_f(2)<=epsilon )  &&  ( Grad_f(3)<=epsilon )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
        break;
    end
    k=k+1;
end
if i == max_iteration
    disp(['超过最大迭代次数！']);
end
disp(['迭代次数为：',num2str(k),'次']);
disp(['最优值点为：']);
x
figure(1)
plot( 2:length(L_v),L_v(2:length(L_v)), '-ro', 'linewidth',2)
disp(['最优值为：',num2str(L),'']);
figure(2)
plot( 2:length(L_v), f_v(2:length(f_v)) - L_v(2:length(L_v)) , '-ro', 'linewidth',2)
figure(3)
plot( 2:length(x_v(1,:)),x_v(1,2:length(x_v)), '-bo', 'linewidth',2)
grid on
xlabel('Iterate')
ylabel('x_1 ')
title(['x_1' ]) 
figure(4)
plot( 2:length(x_v(2,:)),x_v(2,2:length(x_v)), '-go', 'linewidth',2)
grid on
xlabel('Iterate')
ylabel('x_2 ')
title(['x_2' ])

九、直接添加框约束

1）优化问题为：
$\begin{aligned} &\mathop{min}\limits_{x_1,x_2} \ f(x_1,x_2)=2(x_1-1)^2+(x_2+2)^2\\ &s.t. \begin{cases} x_1 + x_2 \leq 1 \\ 2 \leq x_1 \leq 4 \\ -4 \leq x_2 \leq 0 \end{cases} \end{aligned}$

2）拉格朗日函数为：
$\begin{aligned} &L(x_1,x_2,\lambda_1,\lambda_2) \\ &=2(x_1-1)^2+(x_2+2)^2+\lambda(x_1+x_2-1) \end{aligned}$

2）每轮迭代，对 x 直接限值

$\begin{aligned} x_1^{k+1}&= \min\{ \max\{ x_1^{k} - t_2[4(x_1-1)+\lambda] , 2 \} , 4 \} \\ x_2^{k+1}&= \min\{ \max\{ x_2^{k} - t_2[2(x_2+2)+\lambda] , -4 \} , 0 \} \end{aligned}$

代码

clear all;
clc;

x = [100 -120]';   %决策变量初始化   %注：本函数的决策变量为向量！ 

lambda = 1 ;   %拉格朗日乘子初始化 
k=1;
lambda_v = [];  x_v = [];  L_v = [];  f_v = [];
tol = 10^(-5);
t_lambda=0.1;  % 拉格朗日乘子的更新步长
t_x=0.1;      % 决策变量 x 的更新步长
max_iteration = 10000;

for i = 1 : max_iteration
    %% 迭代更新决策变量 x
    x_old = x ; 
    
    x(1) = min( max( x(1) - t_x * (  4 * ( x(1) - 1 ) + lambda  )  ,   2 )  ,  4  ) ;  %% 注意这里可以写解析解，也可以写梯度迭代更新解，也可混合写！
    x(2) = min( max( x(2) - t_x * (  2 * ( x(2) + 2 ) + lambda  )  ,  -4 )  ,  0  ) ; 
%     x(1)  =  x(1) - t_x * (  4 * ( x(1) - 1 ) + lambda  ) ;  %% 注意这里可以写解析解，也可以写梯度迭代更新解，也可混合写！
%     x(2)  =  x(2) - t_x * (  2 * ( x(2) + 2 ) + lambda  ) ; 
    
    %% 迭代更新拉格朗日乘子 lambda 
    lambda = max(   lambda + t_lambda * ( x(1) + x(2) - 1 )  ,  0  );  
    f =  2*( x(1)-1 )^2 + ( x(2)+2 )^2 ; 
    L =  2*( x(1)-1 )^2 + ( x(2)+2 )^2 + lambda*( x(1) + x(2) - 1 );
    
    %% 记录历史更新结果
    x_v = [x_v x];
    lambda_v = [lambda_v lambda];
    f_v = [f_v f]; 
    L_v = [L_v L];
    
    %% 终止条件判断
    % if (k>2) &&( abs(L - L_v(k-1)) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
    if (k>2) &&( abs( f - L ) < tol )&&( abs( x_old(1) - x(1) ) < tol )&&( abs( x_old(2) - x(2) ) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
        break;
    end
    k=k+1;
end
if i == max_iteration
    disp(['超过最大迭代次数！']);
end
disp(['迭代次数为：',num2str(k),'次']);
disp(['最优值点为：']);
x 
figure(1)
plot( 2:length(L_v),L_v(2:length(L_v)), '-ro', 'linewidth',2)
disp(['最优值为：',num2str(L),'']);

figure(2)
plot( 2:length(L_v), f_v(2:length(f_v)) - L_v(2:length(L_v)) , '-b*', 'linewidth',2)