python 拟牛顿法求非线性方程

本章内容：

介绍了无约束和有约束两种类型的非线性优化

一、无约束优化

无约束优化的一般形式为:

如果要求最大值，则

1.1 `fminunc()`

介绍：是MATLAB求解无约束优化的主要函数，算法有：信赖域(trust region)算法和拟牛顿法(quasi-Newton)，详解如下（对算法有要求的可以看看）:

Unconstrained Nonlinear Optimization Algorithmsww2.mathworks.cn

例1、求Banana function的最小值（因其函数图像形似香蕉得名）

一般优化算法中会用到函数的导数信息，故可以将目标函数的导数以函数的形式输入到fminunc中。若不提供导数信息，fminunc内部会用差分代替导数，由于差分方法计算出的导数值与导数函数计算出的值存在误差，所以，若能为fminunc提供目标函数的导数函数信息，更便于计算。

（1）拟牛顿法

函数

function f = BanaFun(x)
f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;

主程序

options = optimoptions('fminunc','display','iter');     %显示迭代过程，由于matlab推荐采用新的设置optimoptions替代optimset，所以这里跟以前有点不同
x = [-1.9,2];                                           %初始迭代点
[x,fval,exitflag,output] = fminunc(@BanaFun,x,options)

结果

> In fminunc (line 395)
%迭代次数  目标函数计算次数 函数值      步长            一阶导数值First-order Iteration  Func-count       f(x)        Step-size       optimality0           3           267.62                      1.23e+031           6          214.416    0.000813405            519  2           9          54.2992              1            331  3          15          5.90157       0.482557           1.46  4          21          5.89006             10           2.58  5          24          5.84193              1           6.56  6          36          4.10288         3.3617           15.2  7          42          4.08488       0.115159           19.9  8          48          3.39007             10           11.3  9          51          3.13363              1           20.7  10          54          2.46989              1           3.18  11          60          2.12375       0.680215           8.62  12          63          1.92978              1           8.85  13          66          1.43583              1            2.9  14          72          1.29696       0.424794            6.5  15          75          1.11139              1           5.84  16          78         0.846834              1           2.07  17          84         0.694787       0.355489           5.72  18          87         0.610241              1           3.12  19          93         0.478524       0.655536           5.39  First-order Iteration  Func-count       f(x)        Step-size       optimality20          96         0.343127              1           1.82  21         102          0.28558            0.5            3.9  22         105         0.194608              1           3.02  23         108         0.135004              1           2.02  24         114        0.0854278       0.431816           2.78  25         117        0.0672604              1           2.68  26         120        0.0377243              1           2.42  27         126        0.0157395       0.424018           1.72  28         129        0.0113891              1           1.66  29         132       0.00400353              1           1.11  30         138      0.000778789       0.478349          0.549  31         141      0.000349468              1          0.452  32         144      3.21971e-05              1          0.113  33         147      6.28136e-06              1         0.0865  34         150      4.74855e-08              1       0.000383  Local minimum found.                                %找到了局部最优点Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.<stopping criteria details>x =                                                 %最优解0.9998    0.9996 fval =                                              %最优函数值4.7485e-08exitflag =1output = 包含以下字段的 struct:iterations: 34funcCount: 150stepsize: 0.0034lssteplength: 1firstorderopt: 3.8341e-04algorithm: 'quasi-newton'               %拟牛顿法        message: 'Local minimum found.…'

（2）信赖域法

函数

function [f,g] = BanaFun(x)
f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
g = [100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2;100*(2*x(2)-2*x(1)^2]); %目标函数的偏导数矩阵

主程序

options = optimoptions('fminunc','Algorithm','trust-region','SpecifyObjectiveGradient',true);
%注意时optimotions不是optimset
%'Algorithm','trust-region'算法选择信赖域法
%调用目标函数导数
x = [-1.9,2];
[x,fval,exitflag,output] = fminunc(@BanaFun,x,options)

1.2`fminsearch()`函数

算法：可变多面体算法

例2、在例1的基础上用fminsearch()函数

函数

function f = BanaFun(x)
f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;

主程序

options = optimset('display','iter');
x = [-1.9,2];
[x,fval,exitflag,output] = fminsearch(@BanaFun,x,options)

结果

优化已终止:当前的 x 满足使用 1.000000e-04 的 OPTIONS.TolX 的终止条件，
F(X) 满足使用 1.000000e-04 的 OPTIONS.TolFun 的收敛条件x =1.0000    1.0000fval =4.0686e-10exitflag =1output = 包含以下字段的 struct:iterations: 114funcCount: 210algorithm: 'Nelder-Mead simplex direct search'message: '优化已终止:…'

其实无约束非线性优化还有许多算法，但都讲篇幅过大，因此这里只讲几种主要的算法。

二、约束最优化

fmincon()是最主要的求解约束最优化的函数

形式：

算法：大规模内点法，SQP算法，基于内点反射信赖域算法等。

约束非线性最优化算法ww2.mathworks.cn

例3、求解如下优化问题

函数

function f = confun(x)
f = -x(1)*x(2)*x(3);

主程序

options = optimoptions('fmincon','Display','iter','Algorithm','sqp');  %选用sqp算法
A = [-1,-2,-2;1,2,2];
b = [0;72];
x0 = [10,10,10];                                                        %初始迭代点
Aeq = [];
beq = [];
lb = [];
ub = [];
nonlcon = [];                                                           %对应标准形式中的ceq
[x,fval,exitflag,output] = fmincon(@confun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)

结果

 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  step    optimality0           4   -1.000000e+03     0.000e+00     1.000e+00     0.000e+00     1.000e+02  1           9   -1.364305e+03     0.000e+00     7.000e-01     3.340e+01     1.674e+02  2          13   -3.291805e+03     0.000e+00     1.000e+00     1.420e+01     1.673e+02  3          17   -3.437592e+03     0.000e+00     1.000e+00     8.824e+00     1.247e+02  4          21   -3.455626e+03     0.000e+00     1.000e+00     2.423e+00     3.448e+00  5          25   -3.455999e+03     0.000e+00     1.000e+00     2.942e-01     1.630e-01  6          35   -3.455999e+03     1.421e-14     1.176e-01     4.405e-03     1.687e-01  7          39   -3.456000e+03     0.000e+00     1.000e+00     2.144e-02     7.989e-02  8          43   -3.456000e+03     0.000e+00     1.000e+00     1.037e-02     1.709e-04  Local minimum found that satisfies the constraints.Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.<stopping criteria details>x =24             12             12       fval =-3456       exitflag =1       output = 包含以下字段的 struct:iterations: 8funcCount: 43algorithm: 'sqp'message: 'Local minimum found that satisfies the con…'constrviolation: 0stepsize: 50/4823lssteplength: 1firstorderopt: 121/707851

三、大规模优化问题举例

例4、求解如下优化问题（含200个变量）

函数

function f = Fun(x)
f = 0;
n = 200'
for ii = 1:n
f = f+(x(ii)-1/ii)^2;
end

主程序

n = 200;
x0 = 10*ones(1,n);     %初始迭代点
options = optimoptions('fminunc','Display','iter','Algorithm','quasi-newton');
[x,fval,exitflag,output] = fminunc(@Fun,x0,options)

结果

fval =1.1140e-14

总结

本文介绍了非线性优化算法，主要分为无约束和有约束两种情况。

优化算法中最重要的是算法的选择，本文并没有对这些算法进行对比，这将在以后的章节中介绍。