障碍法和原对偶内点法解决含有不等式约束的凸优化问题

%%%%%%%%%%%%% 凸优化 二次规划的障碍法 和 原对偶内点发 %%%%%%%% %%%%%%%%%%%%%%李月标 %%%%%%%%%%%%%%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 = A*OrignalX; b = Tempb+1; %使Ax<b 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*(p*x + q) + A'*(1./(b-A*x)); %求梯度 hess = t*p + A'*diag(1./(b-A*x).^2)*A; %求hess矩阵 Xnt = -hess\grad; %求牛顿方向 Lamd2 = -grad'*Xnt; %牛顿减量的平方 output(inter1) = m*u/t; %输出对偶间隙 % 判断停止准则 if Lamd2 <=2*ErrNT %如果找到最优解则变化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+step*Xnt))<=0) step = beta*step; end while(t*(0.5*(x+step*Xnt)'*p*(x+step*Xnt)+q'*(x+step*Xnt))-sum(log(b-A*(x+step*Xnt)))>=t*(0.5*(x)'*p*(x)+q'*(x))-sum(log(b-A*(x)))+alpha*step*grad'*Xnt) step = beta*step; end %更新x的值 x = x + step*Xnt; 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 = A*x - b; %约束函数 gap = -fx'*Lamd; %计算代理对偶间隙 output2(inter2) = gap; %输出代理对偶间隙 t = u*m/gap; %确定t Rdual = p*x+q + A'*Lamd; %对偶残差 % display(num2str(norm(Rdual))); Rcent = -diag(Lamd)*fx-(1/t)*ones(m,1);%中心残差 output3(inter2) = norm(Rdual); %输出对偶残差范数 % 判断是否满足停止条件 if((gap<ErrGAP)&&norm(Rdual)<ErrRe) 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.99*min(-Lamd(row)./DeltLamd(row))); %回溯直线搜索 while(min(A*(x + step*DeltX)-b)>=0) %使f(x+)<0 step = step*beta; end Xnew = x + step*DeltX; LamdNew = Lamd + step*DeltLamd; while(norm([p*(x + step*DeltX)+q + A'*(Lamd + step*DeltLamd);-diag( Lamd + step*DeltLamd)*(A*(x + step*DeltX)-b)-1/t*ones(m,1)])... >(1-alpha*step)*norm([Rdual;Rcent])) display('step:'); step = step*beta; end %更新搜索方向 x = x + step*DeltX; Lamd = Lamd + step*DeltLamd; end figure(2); figure2 = semilogy(output2,'*-'); xlabel('迭代次数'); ylabel('代理对偶间隙'); title('代理对偶间隙与迭代次数关系'); figure(3); figure3 = semilogy(output3,'*-'); xlabel('迭代次数'); ylabel('对偶残差范数'); title('对偶残差范数与迭代次数关系');