[资料]

Matlab采用障碍法及原对偶内点法解决不等式约束凸优化程序

2012-3-6 15:26:50 14176 matlab

0 matlab采用障碍法及原对偶内点法解决不等式约束凸优化问题 %%%%%%%%%%%%% 凸优化二次规划的障碍法和原对偶内点发 %%%%%%%% %%%%%%%%%%%%%%李月标 %%%%%%%%%%%%%%2011210974 %%%%%%%%%%%%%%2011/11/24 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 初始化 A，X，b clear; clc; m = 500; %行 n = 1000; %列 A=randn(m,n); %生成随机矩阵 RankOfA=rank(A); while RankOfA ~= m %判断矩阵的秩，直到符合要求为止 display('A不是满秩矩阵，重新生成矩阵A！'); A=randn(m,n); end display('A是满秩矩阵，符合要求！'); OrignalX = randn(n,1); %OrignalX为严格可行点 Tempb = AOrignalX; b = Tempb+1; %使Ax p = eye(n); %对称正定矩阵 while min(eig(p))<=0 %p = pascal(n); end display('p为对称正定矩阵'); q = randn(n,1); %随机生成q %% 两种方法的公共变量 MaxCaltime = 100; %设置最大计算次数 x = OrignalX; %初始化严格可行点 u=20; %初始化20 alpha = 0.01; %回溯搜索参数 beta = 0.5; %% 采用障碍法 t = 2; ErrNT = 1e-6; %Newton减量误差 ErrGap = 1e-3; %对偶间隙 flag = 1; %用来判断第一次完成外部迭代的变量 %因为第一次迭代时还没确定对偶间隙 %所以完成第一外部迭代后才能获得对偶间隙 for inter1 = 1:MaxCalTime display(num2str(inter1)); % 第一步，中心点步骤（采用newton 回溯直线搜索方法） grad = t(px + q) + A'(1./(b-Ax)); %求梯度 hess = tp + A'diag(1./(b-Ax).^2)A; %求hess矩阵 Xnt = -hessgrad; %求牛顿方向 Lamd2 = -grad'Xnt; %牛顿减量的平方 output(inter1) = mu/t; %输出对偶间隙 % 判断停止准则 if Lamd2 <=2ErrNT %如果找到最优解则变化t flag = flag + 1; if flag == 2 %记录第一次外部迭代的次数，因为这个次数之前的对偶间隙是无效的 fisrtinter = inter1; end display('外部迭代结束一次'); output(inter1) = m/t; %输出对偶间隙 if(m/t < ErrGap) %判断是否达到最优值 break; end t = t * u; %增加t值 continue; end %回溯直线搜索 step = 1; while(min(b-A(x+stepXnt))<=0) step = betastep; end while(t(0.5(x+stepXnt)'p(x+stepXnt)+q'(x+stepXnt))-sum(log(b-A(x+stepXnt)))>=t(0.5(x)'p(x)+q'(x))-sum(log(b-A(x)))+alphastepgrad'Xnt) step = betastep; end %更新x的值 x = x + stepXnt; end figure(1); stairs(fisrtinter:inter1,output(fisrtinter:inter1));%应该用stairs函数 axis([0,inter1+5,0,10e3]); set(gca,'yscale','log'); xlabel('迭代次数'); ylabel('对偶间隙'); title('采用障碍法对偶间隙和Newton迭代次数关系'); %% 采用原对偶内点法 x = OrignalX; %初始化严格可行点 ErrRe = 1e-8; %残差误差 ErrGAP = 1e-6; %对偶间隙误差 Lamd = ones(m,1); % 初始化lamd MaxCalTime = 100; for inter2 = 1:MaxCalTime display(num2str(inter2)); fx = Ax - b; %约束函数 gap = -fx'Lamd; %计算代理对偶间隙 output2(inter2) = gap; %输出代理对偶间隙 t = um/gap; %确定t Rdual = px+q + A'Lamd; %对偶残差 % display(num2str(norm(Rdual))); Rcent = -diag(Lamd)fx-(1/t)ones(m,1);%中心残差 output3(inter2) = norm(Rdual); %输出对偶残差范数 % 判断是否满足停止条件 if((gap break; end % 求delety,即DeltX，DeltLamd sol = -[p,A';-diag(Lamd)A,-diag(fx)][Rdual;Rcent]; %书上P546公式 DeltX = sol(1:n); DeltLamd = sol(n+1:n+m); % 获得初始步长 [row,clo] = find(DeltLamd<0); step = min(0.99,0.99min(-Lamd(row)./DeltLamd(row))); %回溯直线搜索 while(min(A(x + stepDeltX)-b)>=0) %使f(x+)<0 step = stepbeta; end Xnew = x + stepDeltX; LamdNew = Lamd + stepDeltLamd; while(norm([p(x + stepDeltX)+q + A'(Lamd + stepDeltLamd);-diag( Lamd + stepDeltLamd)(A(x + stepDeltX)-b)-1/tones(m,1)])... >(1-alphastep)norm([Rdual;Rcent])) display('step:'); step = stepbeta; end %更新搜索方向 x = x + stepDeltX; Lamd = Lamd + stepDeltLamd; end figure(2); figure2 = semilogy(output2,'-'); xlabel('迭代次数'); ylabel('代理对偶间隙'); title('代理对偶间隙与迭代次数关系'); figure(3); figure3 = semilogy(output3,'-'); xlabel('迭代次数'); ylabel('对偶残差范数'); title('对偶残差范数与迭代次数关系'); 复制代码 0
举报淘帖0 只看该作者相关推荐 • 怎样去实现基于matlab约束优化的惩罚函数法呢 1932 • 怎样去编写一种有约束复合型法MATLAB主程序呢 1490 • MATLAB编程求解优化设计 2218 • 约束优化问题大致分为哪几类 1074 • 凸优化第二章之26对偶锥与广义不等式 905 • 如何用matlab求不等式方程 9175 • 机器学习基础\|深入理解拉格朗日乘子法 4912 • SVM经典入门系列（三） 2745 • Matlab优化工具箱学习 43347 • 优化入门绝顶精品！！！ 4627 2条评论发表评论只看该作者