ADMM算法求解一个简单的例子
2017-02-27 15:39
232 查看
求解下面的带有等式约束和简单的边框约束(box constraints)的优化问题:
minx,y(x−1)2+(y−2)2s.t.0≤x≤3,1≤y≤4,2x+3y=5
结果:
从上图我们可以看到,算法在第5次迭代后就几乎收敛。
minx,y(x−1)2+(y−2)2s.t.0≤x≤3,1≤y≤4,2x+3y=5
% 求解下面的最小化问题: % 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 constant param.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 minimum disp(['[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); end 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; xub = 3; ylb = 1; yub = 4; maxIter = param.maxIter; i = 1; % for plot funval = 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; end plot(iterNum,funval,'-r'); end function fval = compute_fval(x,y) fval = (x-1)^2 + (y-2)^2; end
结果:
[x,y]:0.53846,1.3077
从上图我们可以看到,算法在第5次迭代后就几乎收敛。
相关文章推荐
- 【转】一个用Visual C#做组件的简单例子
- 今天的问题:一个简单的例子,请帮我解开“接口实现Java‘隐藏实现细目’”的迷惑。
- 一个简单的tcp filter的例子
- Struts 教程I:一个用jbuilder X 做的最简单的helloworld的struts例子
- 将本地CSV格式文件内容上传到服务器的一个简单例子
- 一个简单的sturts-menu例子
- [导入]一个简单的用JS调用WebService的例子
- 一个简单的IoC例子(抄袭)
- VB Api简单入门(2)-一个简单的例子
- SYBASE ASE 12.0 上一个横表转纵表的简单例子
- web.config文件自定义配置节的使用方法的一个简单例子
- 一个简单的XML Schema的例子
- 一个简单的用JS调用WebService的例子
- [Struts]学习日记1 - 一个简单的例子
- 蛙蛙推荐: 用web服务传递Dataset的一个简单例子
- 一个简单例子表示fixed functional VS/Assemble VS/HLSI VS的例子
- 一个简单的Eclipse插件开发的例子——HelloWorld【转载】
- 在MFC下使用OpenGL的一个简单的例子
- 一个简单的以太网广播收发例子
- Spring 入门(一个简单的例子)