Matlab最优化实现方法之一维非线性优化：fminbnd

最新推荐文章于 2024-07-17 13:00:09 发布

原创最新推荐文章于 2024-07-17 13:00:09 发布 · 5k 阅读

18 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #线性规划

本文介绍了在Matlab中解决一维非线性优化问题的方法，包括进退法、黄金分割法、斐波那契法、牛顿法和割线法，并详细讲解了fminbnd函数的使用，通过实例展示了如何应用这些方法寻找目标函数的极值点。

前言

在实际建模过程中由于诸多因素的影响，其目标函数或约束条件很难用线性关系式表达，针对此类非线性问题的优化称为非线性规划问题，其基本定义为一个n元函数在一组等式或不等式的约束条件下的极值问题，且目标函数和约束条件至少有一个是未知量的非线性函数。而其中寻求一元函数在某区间上的极点的优化方法，称为一维优化问题。

一维最优化方法

在将实际的物理问题抽象为数学模型时，一般按以下步骤来分析：

1）确定自变量；

2）确定优化目标；

3）确定约束条件；

4）确定自变量范围；

而非线性规划问题的一般形式可以表示为：

$minf\left ( x \right ),x \in E^{n} \\ s.t.\left \{ \begin{matrix}h_{i} \left (x \right )=0 \left ( i=1,2,3,... \right ) \\ g_{j}\left ( x \right )\leqslant 0\left ( j=1,2,3,... \right ) \end{matrix} \right \}$

其中， $x=[x_{1},x_{2},x_{3},...,x_{n}]^{T}$ 称为模型的决策变量， $f$ 称为目标函数， $h_{i}(i=1,2,3,...,m)$ 和 $g_{j}(j=1,2,3,...,l)$ 称为约束函数，前者为等式约束，后者为不等式约束。

常见的一维最优化方法主要有黄金分割法、切线法、插值法等，其基本遵循一个思想：

从某一个出水点 $x^{(0)}$ 出发，沿某个适当选择的方向 $p^{(0)}$ （通常为目标函数的下降方法）进行搜索，得到目标值较小的点 $x^{(1)}$ ；再从 $x^{(1)}$ 出发，沿选择的方向 $p^{(1)}$ 进行搜索，得到比目标函数更小的点 $x^{(2)}$ 。

而如何确定目标函数下一次的寻优方向，我们以 $p^{(k+1)}=\lambda _{k}p^{(k)} ;x^{(k+1)}=x^{(k)}+p^{(k+1)}$ 来表示。而针对一维最优化问题，实质上转化成求解步长因子 $\lambda _{k}$ 的极值问题。

而针对一维搜索问题的基本步骤，常见如下：

1）确定搜索区间；

2）采用逐步缩小区间的方法或函数逼近，确定函数的极值点。

常见的方法主要有进退法、黄金分割法、斐波那契法、牛顿法、割线法等。

进退法

进退法的基本思想是指：

单谷函数 $f\left ( x\right )$ ，其极小值位于搜索区间[a,b]内，对于任意的 $x_{1},x_{2}\in [a,b]$ ，如果 $f\left ( x_{1} \right )< f(x_{2})$ ，则 $[a,x_{2}]$ 为极小点的搜索区间；如果 $f(x_{1})>f(x_{2})$ ，则 $[x_{1},b]$ 为极大点的搜索区间。

进退法的代码示例如下，仅供参考：

function [min_x,max_x]=minJT(f,x0,h0,eps)
% the opotimum method of advance and retreat
% 进退法是用来确定搜索区间的算法
% 本次自定义函数 目的是为了寻求极小值的区间
% 目标函数f
% 初始点x0
% 精度eps
if nargin==3
    eps=1E-6;
end
x1=x0;
k=0;
h=h0;
while 1
    x4=x1+h;
    k=k+1;
    f4=subs(f,symvar(f),x4);
    f1=subs(f,symvar(f),x1);
    if f4<f1
        x2=x1;
        x1=x4;
        f2=f1;
        f1=f4;
        h=2*h;
    else
        if k==1
            h=-h;
            x2=x4;
            f2=f4;
        else
            x3=x2;
            x2=x1;
            x1=x4;
            break
        end
    end
end
min_x=min(x1,x3);
max_x=x1+x3-min_x;
end

黄金分割法

黄金分割法的基本思想与进退法基本一致，其对 $x_{1},x_{2}$ 的取值做出了如下要求：

$\left \{ \begin{matrix}x_{1}=a+\lambda (b-a) \\ x_{2}=b-\lambda(b-a) \end{matrix} \right.$

其中， $\lambda$ 取0.618。

黄金分割法的代码示例如下，仅供参考：

function [x,minf]=minHJ(f,a,b,N,eps)
% golden section method - 黄金分割法
% 该算法的做法是选择x1和x2，使得两点在区间[a b]上的位置是对称的；
% 同时满足，新区间内的点x3要求与区间内已有的点相对于搜索区间是对称的。
% 目标函数f
% 极值区间左端点a
% 极值区间右端点b
% 精度eps
% 目标函数取极小值时的自变量x
% 目标函数的极小值minf
if nargin==3
N=10000;
eps=1E-6;
elseif nargin==4
eps=1E-6;
end
l=a+0.382*(b-a);
u=a+0.618*(b-a);
k=1;
tol=b-a;
while tol>eps && k<N
fl=subs(f,symvar(f),l);
fu=subs(f,symvar(f),u);
if fl>fu
a=l;
l=u;
u=a+0.618*(b-a);
else
b=u;
u=l;
l=a+0.382*(b-a);
end
k=k+1;
tol=abs(b-a);
end
if k==N
disp('找不到最小值')
x=NaN;
minf=NaN;
return
end
x=(a+b)/2;
minf=double(subs(f,symvar(f),x));
end

斐波那契法

斐波那契法与黄金分割法基本思想一致，但其针对 $\lambda$ 的取值不同：

$\lambda =\frac{F_{n-1}}{F_{n}}$ ，其中 $F_{n}$ 为第n个斐波那契数。

斐波那契法的代码示例如下，仅供参考：

function [min_x,min_s,a,b]=fibonacci(s,a0,b0,eps)
% 斐波那契法：1,1,2,3,5,8,13,21,34,55,89,……
% 当用斐波那契法以n个探索点缩短某一区间时，区间长度第一次缩短为Fn-1/Fn，
% 其后各次分别为Fn-2/Fn-1、Fn-2/Fn-3、……、F1/F2。
% 即：区间[a b]内取两点t1，t2，满足t1=a+Fn-1/Fn*（b-a），t2=a+Fn-2/Fn*（b-a）
% 目标函数f（符号表达式）
% 极值区间左端点a0
% 极值区间右端点b0
% 精度eps
% Output：分割后的极值区间左端点a
% Output：分割后的极值区间右端点b
% Output：极小值对应的自变量min_x
% Output：极小值min_s
fn=floor((b0-a0)/eps)+1;
f0=1;
f1=1;
n=1;
while f1<fn 
    f2=f1;
    f1=f1+f0;
    f0=f2;
    n=n+1;
end
f=ones(n,1);
for i=2:n
    f(i+1)=f(i)+f(i-1);
end
k=1;
t1=a0+f(n-2)/f(n)*(b0-a0);
t2=a0+f(n-1)/f(n)*(b0-a0);
while k<n-1
    % func4为自定义函数，即目标函数
    f1=double(subs(s,symvar(s),t1));
    f2=double(subs(s,symvar(s),t2));
    if f1<f2
        b0=t2;
        t2=t1;
        t1=b0+f(n-1-k)/f(n-k)*(a0-b0);
    else
        a0=t1;
        t1=t2;
        t2=a0+f(n-1-k)/f(n-k)*(b0-a0);
    end
    k=k+1;
end
t2=(a0+b0)/2;
t1=a0+(1/2+1E-4)*(b0-a0);
f1=double(subs(s,symvar(s),t1));
f2=double(subs(s,symvar(s),t2));
if f1<f2
    x=t1;
    a=a0;
    b=t1;
else
    x=t2;
    a=t2;
    b=b0;
end
min_x=(a+b)/2;
min_s=double(subs(s,symvar(s),min_x));
end

牛顿法

将目标函数 $f(x)$ 在点 $x^{(k)}$ 处泰勒展开：

$f(x)\approx Q(x)=f(x^{(k)})+\bigtriangledown f(x^{(k)})^{T}(x-x^{(k)})+\frac{1}{2}(x-x^{(k)})^{T}\bigtriangledown^{2}f(x^{(k)})(x-x^{(k)})$

假设目标函数 $f(x)$ 在点 $x^{(k+1)}$ 处取极大值或极小值，则其一阶导数为0：

$\bigtriangledown Q(x^{(k+1)})=\bigtriangledown f(x^{(k)})+\bigtriangledown^{2}f(x^{(k)})(x^{(k+1)}-x^{(k)})=0$

推导可得：

$x^{k+1}=x^{k}-[\bigtriangledown^{2}f(x^{k})]^{-1}\bigtriangledown f(x^{k})$

牛顿法的代码示例如下，仅供参考：

function [x,minf]=nwfun(f,x0,accuracy)
% 牛顿法：求初始点附近的极值点
% 将目标函数二阶泰勒展开，并取其在某点处微分=0，依次求极值
% 得：x1=x0-inv(Hesse Matrix)*一阶偏导
% 目标函数f（符号表达式）
% 初始点x0（列向量）
% 精度accuracy
% 目标函数取最小值时的自变量值x
% 目标函数的最小值min_f
if nargin<3
    accuracy=1E-6;
end
InVar=symvar(f);% 获取自变量Independent Variable
% D1f=zeros(length(InVar),1);
for i=1:length(InVar)
   D1f(i,1)=diff(f,InVar(i));
end
% D2f=zeros(length(InVar),length(InVar));
for i=1:length(InVar)
    for j=1:length(InVar)
        D2f(i,j)=diff(D1f(i),InVar(j));
    end
end
x=x0;
n=0;
while n==0 || norm(g1)>accuracy
g1=zeros(length(D1f),1);
for i=1:length(D1f)
    k=find(InVar==symvar(D1f(i)));
    g1(i,1)=double(subs(D1f(i),symvar(D1f(i)),x(k)'));
end
g2=zeros(size(D2f));
for i=1:length(D1f)
    for j=1:length(D1f)
        k=find(InVar==symvar(D2f(i,j)));
        g2=double(subs(D2f(i,j),symvar(D2f(i,j)),x(k)'));
    end
end
p=-inv(g2)*g1;
x=x+p;
n=n+1;
if n>10000
    disp('已循环10000次，未找到结果，请降低精度')
    break
end
end
x=x-p;
minf=double(subs(f,symvar(f),x'));
end

割线法

割线法的基本思想与牛顿法一致，但避免了牛顿法求导求逆等过程，加快了寻优效率。其迭代公式如下：

$x^{k+1}=x^{k}-\frac{(x^{k}-x^{k-1})f^{'}(x^{k})}{f^{'}(x^{(k)})-f^{''}(x^{k-1})}$

割线法的代码示例如下，仅供参考：

function [x,min_f]=minGX(f,x0,x1,accuracy)
% 割线法：适用于一元多阶函数;求初始点附近的极值点；
% 目标函数f（符号表达式）
% 初始点x0，x1
% 精度accuracy
% 目标函数取最小值时的自变量x
% 目标函数的最小值min_f
if nargin==3
    accuracy=1E-6;
end
d1f=diff(f);
k=0;
tol=1;
while tol>accuracy
    d1fx1=double(subs(d1f,symvar(d1f),x1));
    d1fx0=double(subs(d1f,symvar(d1f),x0));
    x2=x1-(x1-x0)*d1fx1/(d1fx1-d1fx0);
    k=k+1;
    tol=abs(d1fx1);
    x0=x1;
    x1=x2;
end
x=x2;
min_f=double(subs(f,symvar(f),x));
end

一维最优化之Matlab实现方法：fminbnd

有关Matlab命令fminbnd函数的调用格式如下：（详情见mathwork官网的介绍）

x=fminbnd（fun,x1,x2）：返回目标函数fun在区间[x1,x2]上的极小值；

x=fminbnd（fun,x1,x2,options）：options为优化参数选项，可通过optimset设置；

options	说明
Display	off：不显示输出； iter：显示每一次迭代信息； final：显示最终结果；
MaxFunEvals	函数评价所允许的最大次数
MaxIter	函数所允许的最大迭代次数
TolX	x的容忍度

[x,fval]=fminbnd(……)：x为返回的最小值，fval为目标函数最小值；

[x,fval,exitflag]=fminbnd(……)：exitflag为终止迭代条件，其取值说明如下：

exitflag	说明
1	表示函数收敛到解x
0	表示达到了函数最大评价次数或迭代的最大次数
-1	表示函数不收敛解x
-2	表示输入的区间有误，即x1>x2

[x,fval,evitflag,output]=fminbnd(……)：output为最优化输出信息，格式为结构体，其取值及说明如下表所示：

output	说明
iterations	表示算法的跌代次数
funccount	表示函数赋值的次数
algorithm	表示求解线性规划问题所用算法
message	表示算法终止的信息

Example 01：

求解如下最优化问题：

$minf(x)=(x-5)^{2}+7\\ s.t. x\in(0,10)$

根据需要建立Matlab命令流，如下：

[x,fval,exitflag,output]=fminbnd(@(x)(x-5)^2+7,0,10)

结果显示如下：