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

分别采用障碍法和原对偶内点法对含有等式和不等式约束的凸优化问题用matlab进行求解

%%%%%%%%%%%%%   凸优化 二次规划的障碍法 和 原对偶内点发 %%%%%%%%
%%%%%%%%%%%%%%李月标
%%%%%%%%%%%%%%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 < 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: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