matlab编写的进退法，黄金分割法，抛物线法（二次插值法），wolfe不精确一维搜索

这是我最优化方法课程的编程作业，贴来和大家分享，后续会继续发来一些最优化的程序。、

以下程序由matlab编写

程序简介

jintuifa.m
进退法，用于确定下单峰区间.根据最优化方法（天津大学出版社）20页算法1.4.3编写。
输出：[left right]  为下单峰区间
输入：y  x  x0 step                 
      y为函数，x为函数y的变量,x0  step（>0）分别为初始点,初始步长

golddiv.m
黄金分割法。根据最优化方法（天津大学出版社）17页算法1.4.2编写。
当保留的区间长度|b-a|<=epsilon时停止迭代
输出：[best_x best_fx]  best_x为最优的x值，best_fx为最优的函数值
输入： y x a b epsilon     
       y为函数，x为函数y的变量,a,b为下单峰区间[a,b],epsilon为精确度

paowuxianfa.m
抛物线法（二次插值法）。根据最优化方法（天津大学出版社）22页算法1.4.5编写。
输出：[best_x best_fx]   best_x为最优的x值，best_fx为最优的函数值
输入：x1 x0 x2 epsilon1 epsilon2   x1,x0,x2为已知的三点且满足f1>f0<f2,epsilon为精确度

wolfe.m
Wolfe不精确一维搜索。根据最优化方法（天津大学出版社）24页算法1.4.6编写。
输出：[alpha xk1]  alpha为要求的步长，xk1为x(k+1)是得出的下一个点
输入：xk pk mu sigma   xk为初始点,pk为方向，mu和sigma为参数. 一般mu属于（0，1/2）,sigma属于(mu,1)
f.m
老师测试的例子函数例4
g.m
老师作业例子4函数求梯度

test.m
一个调用jintuifa.m golddiv.m paowuxianfa.m的计算的例子

程序详细代码

f.m
%最优化编程作业中的例子
 function result = f(x)
 %result = x.^4+2.*x+4;%例1 区间为[-1,0]
 %result = exp(-x)+x.^2;%例2 区间为[0,1]
result =100*(x(2)-x(1)^2)^2+(1-x(1))^2;%例4

g.m

%求在点x处的梯度，对应编程作业中的例四 这里的x是向量，代表x1,x2……
%author liuxi BIT
function result =g(x)

syms x1 x2;%定义符号变量
f=100*(x2-x1^2)^2+(1-x1)^2;
result_1=diff(f,x1);%函数f对符号变量x1求一阶微分
result_2=diff(f,x2);

result(1)=subs(result_1,[x1 x2],[x(1) x(2)]);%将result_1中的符号变量x1和x2替换为普通变量x(1),x(2)
result(2)=subs(result_2,[x1 x2],[x(1) x(2)]);

golddiv.m

%黄金分割法。根据最优化方法（天津大学出版社）17页算法1.4.2编写。
%v1.0 author: liuxi BIT
%best_x为最优的x值，best_fx为最优的函数值，y为函数，x为函数y的变量，a,b为下单峰区间[a,b],epsilon为精确度
%当保留的区间长度|b-a|<=epsilon时停止迭代
function [best_x best_fx]=golddiv(y,x,a,b,epsilon)
if nargin==4
    epsilon=0.0000001;;%设置默认的epsilon
end

x2=a+0.618*(b-a);%step1
f2=subs(y,x,x2);

x1=a+0.382*(b-a);%step2
f1=subs(y,x,x1);

while(abs(b-a)>epsilon)
    if f1<f2 %step4
        b=x2;x2=x1;f2=f1;
        x1=a+0.382*(b-a);%转step2
        f1=subs(y,x,x1);
    elseif f1==f2
        a=x1;b=x2;
        x2=a+0.618*(b-a);f2=subs(y,x,x2);%转step1
        x1=a+0.382*(b-a);f1=subs(y,x,x1);
    else
        a=x1;x1=x2;f1=f2;
        x2=a+0.618*(b-a);%step5
        f2=subs(y,x,x2);
    end
end%while
best_x=(a+b)/2;%最优的x值
best_fx=subs(y,x,best_x);%最优的函数值

jintuifa.m

%进退法，用于确定下单峰区间.根据最优化方法（天津大学出版社）20页算法1.4.3编写。
%v1.0 author: liuxi BIT
%[left right]为下单峰区间，y为函数，x为函数y的变量，x0为初始点（默认为0）,step（>0）为初始步长（默认为0.01）
function [left right] = jintuifa(y,x,x0,step)
if nargin==3%当只有三个参数时，默认设置步长为0.01
    step=0.01;
end
if nargin==2
    x0=0;%当只有两个参数时，默认设置初始点为0
    step=0.01;
end

f0=subs(y,{x},{x0});%step1 求f(x0) 将函数y中变量x替换为x0
x1=x0+step;%step2
f1=subs(y,{x},{x1});

if (f1<=f0)%step3 step4
    step=2*step;
    x2=x1+step;
    f2=subs(y,{x},{x2});
    while(f1>f2)
        x0=x1; x1=x2;
        f0=f1; f1=f2;
        step=2*step;
        x2=x1+step;
        f2=subs(y,{x},{x2});
    end
    left=x0
    right=x2
else%step5 step6
    step=2*step;
    x2=x1-step;
    f2=subs(y,{x},{x2});
    while(f0>f2)
        x1=x0; x0=x2;
        f1=f0; f0=f2;
        step=2*step;
        x2=x1-step;
        f2=subs(y,{x},{x2});
    end
    left=x2;right=x1;
end

 

paowuxianfa.m

%抛物线法（二次插值法）。根据最优化方法（天津大学出版社）22页算法1.4.5编写。
%v1.0 author: liuxi BIT
%best_x为最优的x值，best_fx为最优的函数值，y为函数，x为函数y的变量，x1,x0,x2为已知的三点且满足f1>f0<f2,epsilon为精确度
function[best_x best_fx]=paowuxianfa(y,x,x1,x0,x2,epsilon1,epsilon2)
if nargin==5
    epsilon1=0.0000001;
    epsilon2=0.0000001;
end
%求f(x1),f(x0),f(x2)
f1=subs(y,x,x1);
f0=subs(y,x,x0);
f2=subs(y,x,x2);
k=0
while abs(x1-x2)>epsilon1&&abs((x2-x0)*f1+(x1-x2)*f0+(x0-x1)*f2)>epsilon2%step1 step2
    x3=0.5*((x2^2-x0^2)*f1+(x1^2-x2^2)*f0+(x0^2-x1^2)*f2)/((x2-x0)*f1+(x1-x2)*f0+(x0-x1)*f2);%step3   x3对应算法中的x一拔，就是x上面有一个“—”
    f3=subs(y,x,x3);
    if f0-f3<0%转step6
        if x0<x3
            x2=x3;f2=f3;
        else
            x1=x3;f1=f3;
        end%if x0<x3
    elseif f0-f3==0%转step7
        if x0<x3
            x1=x0;x2=x3;x0=0.5*(x1+x2);
            f1=f0;f2=f3;f0=subs(y,x,x0);
        elseif x0==x3%转step8
            x4=0.5*(x1+x0);f4=subs(y,x,x4);%x4对应算法中的x一冒，就是x上面有一个“^”
            if f4<f0
                x2=x0;x0=x4,f2=f0;f0=f4;
            elseif f4==f0
                x1=x4;x2=x0;x0=0.5*(x1+x2);
                f1=f4;f2=f0,f0=subs(y,x,x0);
            else
                x1=x4;f1=f4;
            end%if f4<f0
        else%转step9
            x1=x3;x2=x0;x0=0.5*(x1+x2);
             f1=f3;f2=f0;f0=subs(y,x,x0);
        end% if x0<x3
    else %转step5
        if x0>x3
            x2=x0;x0=x3;f2=f0;f0=f3;
        else
            x1=x0;x0=x3;f1=f0;f0=f3;
        end
    end%if f0-f3<0
    k=k+1
end%while
best_x=x0;
best_fx=f0;

wolfe.m

%Wolfe不精确一维搜索。根据最优化方法（天津大学出版社）24页算法1.4.6编写。
%v1.0 author: liuxi BIT
%alpha为要求的步长，xk1为x(k+1)是得出的下一个点，xk为初始点,pk为方向，mu和sigma为参数，
%一般mu属于（0，1/2）,sigma属于(mu,1)
function[alpha xk1] =wolfe(xk,pk,mu,sigma)
if nargin<=2
    mu=0.1;%默认设置
    sigma=0.5;
end
if nargin==0
    xk=[0 0];%这个是为了测试的那个例子方便而设的，通用性就不需要这段了
    pk=[1 0];
end

a=0;b=inf;alpha=1;j=0;%step1 j用来标记迭代次数 Inf表示无穷大
xk1=xk+alpha*pk;%step2 xk1代表算法里的x(k+1) , fk1代表f(k+1), gk1代表g(k+1)
fk1=f(xk1);
gk1=g(xk1);
fk=f(xk);
gk=g(xk);

while  fk-fk1<-mu*alpha*gk*pk'
    j=j+1;
    b=alpha;
    alpha=0.5*(alpha+a);
    xk1=xk+alpha*pk;
    fk1=f(xk1);
    %gk1=g(xk1);
end

while gk1*pk'<sigma*gk*pk'
    j=j+1;
    a=alpha;
    alpha=min([2*alpha 0.5*(a+b)]);
    xk1=xk+alpha*pk;
    %fk1=f(xk1);
    gk1=g(xk1);
end
xk1=xk+alpha*pk;

test.m

%format long
syms x;
y = x.^4+2.*x+4;
[left right] = jintuifa(y,x,-1);
[best_x best_fx]=golddiv(y,x,left,right);
[best_x best_fx]=paowuxianfa(y,x,left,(left+right)/2,right);