%%%%%%%%%%%%%???凸優化?二次規劃的障礙法?和?原對偶內點發?%%%%%%%%
%%%%%%%%%%%%%%李月標
%%%%%%%%%%%%%%2011210974
%%%%%%%%%%%%%%2011/11/24
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%?初始化?A，X，b
clear;
clc;
m?=?500;????????????????%行
n?=?1000;???????????????%列
A=randn（mn）;???????????%生成隨機矩陣
RankOfA=rank（A）;
while?RankOfA?~=?m???????????%判斷矩陣的秩，直到符合要求為止
????display（‘A不是滿秩矩陣，重新生成矩陣A！‘）;
????A=randn（mn）;
end
display（‘A是滿秩矩陣，符合要求！‘）;
OrignalX?=?randn（n1）;??????%OrignalX為嚴格可行點
Tempb?=?A*OrignalX;
b?=?Tempb+1;????????????????%使Axp?=?eye（n）;??????????????%對稱正定矩陣
while?min（eig（p））<=0
????%p?=?pascal（n）;??
end
display（‘p為對稱正定矩陣‘）;
q?=?randn（n1）;?????????????%隨機生成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?????????????????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:inter1output（fisrtinter:inter1））;%應該用stairs函數
axis（[0inter1+5010e3]）;
set（gca‘yscale‘‘log‘）;
xlabel（‘迭代次數‘）;
ylabel（‘對偶間隙‘）;
title（‘采用障礙法對偶間隙和Newton迭代次數關系‘）;
%%?采用原對偶內點法
x?=?OrignalX;???????????????????????%初始化嚴格可行點
ErrRe?=?1e-8;???????????????????????%殘差誤差
ErrGAP?=?1e-6;??????????????????????%對偶間隙誤差
Lamd?=?ones（m1）;???????????????????%?初始化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;????????%對偶殘差
???%?dis

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????4729??2011-11-26?17:04??Convex11.m

-----------?---------??----------?-----??----

?????????????????4729????????????????????1