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

yuyang911220

• %%%%%%%%%%%%%   凸优化 二次规划的障碍法 和 原对偶内点发 %%%%%%%%
• %%%%%%%%%%%%%%李月标
• %%%%%%%%%%%%%%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('对偶残差范数与迭代次数关系');