ADMM算法求解一个简单的例子

2019-11-06 09:16:52

字体：大中小

来源：转载

供稿：网友

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

minx,y(x−1)2+(y−2)2s.t.0≤x≤3,1≤y≤4,2x+3y=5 $/begin{equation}/begin{aligned}/min_{x,y} (x-1)^2+(y-2)^2 //s.t. 0/leq x /leq 3,//1 /leq y /leq 4, //2x+3y=5/end{aligned}/end{equation}$

% 求解下面的最小化问题：% 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% solve y min f(y) = (y-2)^2 + lambda*(2x + 3y -5) + rho/2(2x + 3y -5)^2, s.t. 1<= y <=4% sovle lambda :update lambda = lambda + rho(2x + 3y -5)% 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;param.beta = 1.1; % a constantparam.rho = 0.5;param.rhomax = 2000;[Hx,Fx] = getHession_F('f1');[Hy,Fy] = getHession_F('f2');param.Hx = Hx;param.Fx = Fx;param.Hy = Hy;param.Fy = Fy;% solve PRoblem using admm algrithm[x,y] = solve_admm(param);% disp minimumdisp(['[x,y]:' num2str(x) ',' num2str(y)]);function [H,F] = getHession_F(fn)% fn : function name% H :hessian matrix% F :一次项系数syms x y lambda rho;if strcmp(fn,'f1') f = (x-1)^2 + lambda*(2*x + 3*y -5) + rho/2*(2*x + 3*y -5)^2; H = hessian(f,x); F = (2*lambda + (rho*(12*y - 20))/2 - 2); %%%%%%%%%%%%%%%%%%%%%%%%%%%% fcol = collect(f,{'x'}); disp(fcol);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);endfunction [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;xub = 3;ylb = 1;yub = 4;maxIter = param.maxIter;i = 1;% for plotfunval = zeros(maxIter-1,1);iterNum = zeros(maxIter-1,1);while 1 if i == maxIter break; end % solve x Hxx = eval(Hx); Fxx = eval(Fx); x = quadprog(Hxx,Fxx,[],[],[],[],xlb,xub,[],options);% descend % solve y Hyy = eval(Hy); Fyy = eval(Fy); y = quadprog(Hyy,Fyy,[],[],[],[],ylb,yub,[],options);%descend % update lambda lambda = lambda + rho*(2*x + 3*y -5); % ascend % rho = min(rhomax,beta*rho); funval(i) = compute_fval(x,y); iterNum(i) = i; i = i + 1;endplot(iterNum,funval,'-r');endfunction fval = compute_fval(x,y)fval = (x-1)^2 + (y-2)^2;end

结果：

[x,y]:0.53846,1.3077

这里写图片描述