function PHR()
global lumda sigma  l m;
x=[0 2];            %初始點
l=0;                %等式
m=1;                %不等式
lumda=1;            %乘子向量
sigma=5;            %初始罰因子
g=2;                %放大系數(shù)
eps=1e-16;          %控制誤差
sita=0.5; 
temp=terminate(x);
while 1
    c=confun(x);
    x=fminunc(@M, x);
    fai=terminate(x);
    if fai<eps
        opt=x;
        break;
    end
    if fai/temp>=sita
        sigma=g*sigma;
    end
    temp=fai;
    for i=1:m
        for j=1:l
            lumda(i)=lumda(i)-sigma*c(i);
            for j=(l+1):m
                if lumda(i)-sigma*c(i)>=0
                    lumda(i)=lumda(i)-sigma*c(i);
                else
                    lumda(i)=0;
                end
            end
        end
    end
end