外点罚函数法的matlab编程与使用

一 函数使用说明

外点罚函数penalty(fun,g,h,x0)，可以用来求解含有非线性不等式和等式约束的优化问题。在调用函数时首先需要用户新建三个函数即目标函数，不等式约束和等式约束，然后输入penalty(@fun,@g,@h,x)即可，fun,g,h为用户自己定义的三个函数的名称，x为初始点坐标。

二 主函数说明

penalty(x)函数主要利用外点罚函数法来对目标函数进行优化，首先根据用户给出的三个函数（fun(x)，g(x)，h(x)）构造一个新的函数TextF(x)，即
TextF(x)=fun(x)+r ( ∑_(j=1)^m▒〖max[0,g_j (x)]^2 〗 + ∑_(k=1)^l▒〖[h_k (x)]^2 〗 )。
根据经验值确定罚因子r0=1，权数c=5，精度0.001。每次迭代都由本轮的r值得出一个最优点（使用鲍威尔法），当步长足够小时结束迭代，当r大于障碍数R时，令r=Inf,再进行一次迭代后结束迭代，结束迭代后返回当前x和fx。

三 子函数说明

fun函数为用户输入的目标函数。
g函数为用户输入的不等式约束函数。
h函数为用户输入的等式约束函数。
PowellMethod函数主要利用鲍威尔法来求得每个r值对应目标函数的最优点。鲍威尔法每轮沿其搜索组进行一次搜索，并获得一个新的方向，为了保证每个方向线性不相关，根据判别式来确定新方向是否替换原来的一个方向。满足终止条件或者迭代次数过多，结束迭代，迭代结束后返回最优点x和fx。
OneDimensionSearch函数是对给定的函数进行一次一维搜索，需要给定初始点和搜索方向。先进行一次外推法得到搜索区间，然后利用黄金分割法来确定本次一维搜索的最优点和最优值。
waitui函数是利用初始点和初始步长来进行搜索和判断，当函数值出现“等等”，“高低高”或者“高等等”时即可得到函数最优点的区间范围。
golden函数根据最优点的区间，对区间进行黄金分割，第一轮进行两次黄金分割，后面每轮进行一次黄金分割和替换，逐步逼近最优点的位置。当时搜索步长足够小或者迭代次数太多时结束迭代，迭代结束后返回最优点x和fx。
四 实例

外点罚函数法主程序
function [X,FX]=penalty(fun,g,h,x0)
%外点罚函数法，需要给出fun函数，g函数，和h函数以及初始点
%设置全局变量r，给出罚因子r0，权数为5，设置精度，定义全局TextF函数
global TextF;
TextF=@(x)fun(x)+r*(sum(max(g(x),0).^2)+sum(h(x).^2));
r=1;
c=5;
x=x0;
sign=0;
err=0.001;
Lastx=x;
Lastfx=TextF(x);
%开始迭代
while(1)
%设置上一次迭代的x值为Lastx，改变罚因子的值，并得到本次迭代最优值x
    Lastx=x;
r=c*r;    %每次循环重新定义TextF函数
TextF=@(x)fun(x)+r*(sum(max(g(x),0).^2)+sum(h(x).^2));
    sign=sign+1;
    [x,fx]=PowellMethod(Lastx);
    %当迭代次数大于100次或者满足精度时返回此时的x和x对应的函数值
    if sign >100
        r=Inf;
        [x,fx]=PowellMethod(Lastx);
        break;
    end
    if norm(x-Lastx)<err
        break;
    end 
end
X=x;
FX=fx;
end

目标函数
function y= fun(x)
%目标函数
y=fun(x);
end


不等式约束
function y= g(x)
%新建函数为m个不等式约束gi<=0,如果没有不等式约束，则令y=0
y(1)=g1(x);
y(2)=g2(x);
...
y(m)=gm(x);
end


等式约束
function y= h(x)
%新建函数为L个等式约束hi=0,如果没有等式约束，则令y=0
y(1)=h1(x);
y(2)=h2(x);
...
y(l)=hl(x);
end

鲍威尔法
function [Matrix,Value]=PowellMethod(X0)
global TextF;     %声明全局函数TextF
%设置初始迭代轮次k=1；迭代精度：Err；初始点；初始搜索方向组；
n=length(X0);
S=eye(n);
Err=0.001;
flag=0;
JZ=zeros(1,3); %预定义F0\F2\F3目标函数值
X=zeros(3,2);  %预定义F0\F2\F3对应的极值点坐标（x1，x2）

%迭代循环计算
while(1)
    flag=flag+1;
%从X0出发，依次沿Si(i=1,2,..,n)进行n次一维搜索，得到n个一维极小点
    [X(1,:),JZ(1)] = OneDimensionSearch(X0,S(1,:),Err);
    for i=2:n
        [X(i,:),JZ(i)] = OneDimensionSearch(X(i-1,:),S(i,:),Err);
    end

    %连接X0,Xn，构成新的共轭方向SK
    SK=X(n,:)-X0;
    %沿共轭方向SK计算X0的新映射点
    Xmap=2*X(n,:)-X0;

    %计算K轮初始点、终点和映射点的函数值
    Val_0 = TextF(X0);
    Val_1 = TextF(X(n,:));
    Val_2 = TextF(Xmap);
    %计算k轮中各相邻极小点目标函数的差值，并求出其中的最大差值及其相应的方向
    Diff=zeros(1,2);
    Diff(1) = Val_0-JZ(1);
    for k = 2:n
        Diff(k) =JZ(k-1)-JZ(k);
    end
    [maxDiff,m] = max(Diff);  

    %用判别条件式检验原方向组是否需要替换
    if(Val_2<Val_0 && (Val_0-2*Val_1+Val_2)*(Val_0-Val_1-maxDiff)<0.5*maxDiff*(Val_0-Val_2)^2 )
        %满足判别式,即为非线性相关 
        %由Xn出发沿方向SK进行一维搜索，求出该方向的极小点X(k),并将X(k)作为K+1轮迭代的初始点
        [X(n+1,:),JZ(n+1)] = OneDimensionSearch(X(n,:),SK(1,:),Err);
        KX0=X0;   %K轮初始点，即是本轮初始点
        X0 = X(n+1,:);  %K+1轮初始点，即是下轮初始点
        %去掉方向S(mk)，而将方向SK作为K+1轮迭代的最末一个方向。
        S(m,:)=[];
        S=[S;SK];

        %判断迭代终止条件,若满足判别条件，则跳出循环，输出最优值，否则进行K+1轮迭代
        delta = X(n,:)-KX0;
        if( sqrt(sum(delta.*delta))<=Err )
            Matrix = X(n+1,:);
            Value = TextF(Matrix); %输出最优解和最优变换参数
            break;
        end
    %判别条件式不成立，重新赋初始值，使用原来的方向组，进行下轮循环    
    else
        if Val_2>Val_1
            X0 = X(n,:);
        else
            X0 = Xmap;
         
        end
        %判断迭代终止条件
        delta = X(n+1,:)-X0;
        if( sqrt(sum(delta.*delta))<=Err )
            Matrix = X(n+1,:);
            Value = TextF(Matrix); %输出最优解和最优变换参数
            break;
        end
    end
    %搜索次数大于300次时，停止迭代
    if flag>300
        Matrix = X(n,:);
        Value = TextF(Matrix); 
        break;
    end
end
Matrix;
Value;
end

一维搜索
function [x,fx]= OneDimensionSearch(X0,S,Err)
%一维搜索，给出初始点，搜索方向和精度
%先进行一次外推初始步长为0.01S，再进行黄金分割，精度为,001
[a,b]=waitui(X0,0.01*S); 
[x,fx]=golden(a,b,0.001);
end 

外推法
function [a,b]=waitui(x0,h)
global TextF;   %声明全局函数TextF
%外推法的基本形式，需要给出初始点和搜索方向
b=x0+h;a=x0;
if  TextF(b)~=TextF(a)      %当第一次搜索两函数值不相等开始搜索，否则返回
%第一次的搜索方向错误时，改变方向
if TextF(b)>TextF(x0)
    h=-h;
    x1=x0-h;
    x2=x0;
    x3=x0+h;
else
    x1=x0;
    x2=x0+h;
    x3=x2+h;
end
%开始搜索，步长依次加倍，高低高或者高等等时停止
while TextF(x3)<TextF(x2)
    h=2*h;
    x1=x2;
    x2=x3;
    x3=x2+h;
    if(TextF(x3)==TextF(x2))
        break;
    end
end
a=x1;
b=x3;
end

黄金分割法
function [min,fmin]=golden(low,up,length)
global TextF;   %声明全局函数TextF
%黄金分割法初始化，需要给出最下限和最上限以及精度
%第0轮
flag=0;
a1=low;a3=up;
f1=TextF(a1);f3=TextF(a3);
a11=a1+0.382*(a3-a1);f11=TextF(a11);
a12=a1+0.618*(a3-a1);f12=TextF(a12);
lowbound=a1;
upbound=a3;
%第flag轮，结束循环标志为步长足够小，或者进行了300轮
while norm(a3-a1)>=length
    flag=flag+1;
    lowbound=a1;
    upbound=a3;
    if f12>f11
      a3=a12;f3=f12;
      a12=a11;f12=f11;
      a11=a1+0.382*(a3-a1);
      f11=TextF(a11);
    else
      a1=a11;f1=f11;
      a11=a12;f11=f12;
      a12=a1+0.618*(a3-a1);
      f12=TextF(a12);
    end
    if flag>300
        break;
    end
end
min=(lowbound+upbound)/2;
fmin=TextF(min);