电子技术论坛_中国专业的电子工程师学习交流社区-中电网技术论坛

标题: Matlab采用障碍法及原对偶内点法解决不等式约束凸优化程序 [打印本页]

作者: yuyang911220 时间: 2017-2-23 17:04 标题: 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 = 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('对偶残差范数与迭代次数关系');

欢迎光临电子技术论坛_中国专业的电子工程师学习交流社区-中电网技术论坛 (http://bbs.eccn.com/)