admm代码-部分转载-作为自用笔记

求解下面的带有等式约束和简单的边框约束(box constraints)的优化问题：

↓主程序

% 求解下面的最小化问题：
% min_{x,y} (x-1)^2 + (y-2)^2
% s.t. 0 \leq x \leq 3
%      1 \leq y \leq 4
%      2x + 3y = 5

% augumented lagrangian function:
% L_{rho}(x,y,lambda) = (x-1)^2 + (y-2)^2 + lambda*(2x + 3y -5) + rho/2(2x + 3y -5)^2  
%↑构建拉格朗日函数
% solve x  min f(x) = (x-1)^2 +  lambda*(2x + 3y -5) + rho/2(2x + 3y -5)^2,s.t. 0<= x <=3
%↑分成了求解x的子问题
% solve y  min f(y) = (y-2)^2 + lambda*(2x + 3y -5) + rho/2(2x + 3y -5)^2, s.t. 1<= y <=4
%↑分成了求解y的子问题
% sovle lambda :update lambda = lambda + rho(2x + 3y -5)
%↑对lambda进行更新
% rho ,we set rho = min(rho_max,beta*rho) beta is a constant ,we set to 1.1,rho0 = 0.5;

%↓在设置参数（结构数组）
% x0,y0 都为1是一个可行解。
param.x0 = 1;
param.y0 = 1;
param.lambda = 1;
param.maxIter = 30;%迭代次数最多为30
param.beta = 1.1; % a constant
param.rho = 0.5;
param.rhomax = 2000;

[Hx,Fx] = getHession_F('f1');%调用了子程序getHession_F，求解出了x子问题里的H阵和F，并完成了赋值（不过现在得到的F里面还包含着字母y和lambda，要等solve_admm子程序里用eval函数将它编程一个数值
[Hy,Fy] = getHession_F('f2');

param.Hx = Hx;
param.Fx = Fx;
param.Hy = Hy;
param.Fy = Fy;

% solve problem using admm algrithm，通过下面这个子程序用admm算法求解x和y的最优解
[x,y] = solve_admm(param);

% disp minimum
disp(['[x,y]:' num2str(x) ',' num2str(y)]);

求解H，F（因为用quadprog要知道H阵和一次项系数F

function [H,F] = getHession_F(fn)
% fn : function name
% H :hessian matrix
% F :一次项系数，具体可以看goodnotes里笔记
syms x y lambda rho;%syms意为定义多个变量
%↓strcmp函数在比较两个字符串，如果两个字符串相等则输出0（但是我觉得strcmp的用法有点问题
if strcmp(fn,'f1')
    f = (x-1)^2 + lambda*(2*x + 3*y -5) + rho/2*(2*x + 3*y -5)^2;%这里是在构建子问题1里的增广拉格朗日函数
    H = hessian(f,x);%这里是在求H矩阵
    F = (2*lambda + (rho*(12*y - 20))/2 - 2);%这里F的求法有没有更智能的选项？（难道说其实是先运行了下面的collect函数再写了F？
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fcol = collect(f,{'x'});%固定y，默认x为符号变量，按照x相同次幂的合并同类项
    disp(fcol);%直接将内容输出在Matlab命令窗口中
elseif strcmp(fn,'f2')
    f = (y-2)^2 + lambda*(2*x + 3*y -5) + rho/2*(2*x + 3*y -5)^2;
    H = hessian(f,y);
    F = (3*lambda + (rho*(12*x - 30))/2 - 4);
    fcol = collect(f,{'y'});
    disp(fcol);
end
% fcol = collect(f,{'x'});
% fcol = collect(f,{'y'});
% disp(fcol);
end

↓使用quadprog进行交替求解x,y
分解为了x相关的二次规划子问题和y相关的二次规划子问题，子问题中二次规划的H阵和一次项系数都包含了增广拉格朗日函数的影响

%param包含了多个参数，使用param可以将很多个参数封装在一起，这样函数调用的时候比较简洁，也利于编程和维护
%param是个结构数组，把一组彼此相关但类型不同的数据组织在一起，具体可参考帮助文档中的struct命令。
function [x,y] = solve_admm(param)
x = param.x0;
y = param.y0;
lambda = param.lambda;
beta = param.beta;
rho = param.rho;
rhomax = param.rhomax;
Hx = param.Hx;
Fx = param.Fx;
Hy = param.Hy;
Fy = param.Fy;
%%
options = optimoptions('quadprog','Algorithm','interior-point-convex');%在设置优化选项，但没看太懂
xlb = 0;%xlb≤x≤xub
xub = 3;
ylb = 1;%ylb≤y≤yub
yub = 4;
maxIter = param.maxIter;
i = 1;
% for plot
funval = zeros(maxIter-1,1);%运行到第maxIter会直接break，根本没有任何数值可以记录，所以是maxIter-1
iterNum = zeros(maxIter-1,1);
while 1
    if i == maxIter
        break;
    end


    % solve x，求解x子问题，这是一个二次规划问题，同时二次规划的H阵和一次项系数都包含了增广拉格朗日函数的影响
    Hxx = eval(Hx);%eval函数的功能是将字符串转换为matlab可执行语句→用循环末尾更新的拉格朗日乘子lambda、上一次求解出的x和y来更新H阵和一次项系数F，使之成为一个数/纯数矩阵，不包含字母
    Fxx = eval(Fx);
    x = quadprog(Hxx,Fxx,[],[],[],[],xlb,xub,[],options);% descend；调用格式见goodnotes里二次规划笔记


    % solve y，求解y子问题
    Hyy = eval(Hy);
    Fyy = eval(Fy);
    y = quadprog(Hyy,Fyy,[],[],[],[],ylb,yub,[],options);%descend


    % update lambda，在每次迭代后更新拉格朗日乘子lambda（用子问题里找到的x、y的最优解来更新，lambda是个数值
    lambda = lambda + rho*(2*x + 3*y -5); % ascend
    
    
    % rho = min(rhomax,beta*rho);这里是注释了自适应补偿？并没有更新惩罚系数rho
    funval(i) = compute_fval(x,y);
    iterNum(i) = i;
    i = i + 1;
end

plot(iterNum,funval,'-r');

end

function  fval  = compute_fval(x,y)
fval = (x-1)^2 + (y-2)^2;
end