Main Content

Optimization Toolbox 教程

本教程包括多个示例,说明如何使用两个非线性优化求解器 fminuncfmincon,以及如何设置选项。此示例中概述的原则也适用于其他非线性求解器,如 fgoalattainfminimaxlsqnonlinlsqcurvefitfsolve

教程示例涵盖下列任务:

  • 最小化一个目标函数

  • 最小化带附加参数的同一函数

  • 最小化带约束的目标函数

  • 通过提供梯度或黑塞矩阵,或通过更改选项来获得更高效或更准确的解。

无约束优化示例

假设要求下面这个函数的最小值:

xexp(-(x2+y2))+(x2+y2)/20.

绘制函数以查看它具有最小值的位置。

f = @(x,y) x.*exp(-x.^2-y.^2)+(x.^2+y.^2)/20;
fsurf(f,[-2,2],'ShowContours','on')

绘图显示最小值在点 (-1/2,0) 附近。

通常将目标函数定义为一个 MATLAB® 文件。在这种情况下,函数足够简单,可以定义为匿名函数。

fun = @(x) f(x(1),x(2));

为找到解设置一个初始点。

x0 = [-.5; 0];

设置优化选项以使用 fminunc 默认 'quasi-newton' 算法。此步骤可确保教程在每个 MATLAB 版本中的工作方式是相同的。

options = optimoptions('fminunc','Algorithm','quasi-newton');

在求解器执行计算时查看迭代。

options.Display = 'iter';

调用 fminunc,它是一个无约束非线性最小化求解函数。

[x, fval, exitflag, output] = fminunc(fun,x0,options);
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           3          -0.3769                         0.339
     1           6        -0.379694              1          0.286  
     2           9        -0.405023              1         0.0284  
     3          12        -0.405233              1        0.00386  
     4          15        -0.405237              1       3.17e-05  
     5          18        -0.405237              1       3.35e-08  

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

显示求解器找到的解。

uncx = x
uncx = 2×1

   -0.6691
    0.0000

查看在该解处的函数值。

uncf = fval
uncf = -0.4052

这些示例使用函数计算的次数来衡量求解效率。查看函数计算的总次数。

output.funcCount
ans = 18

带附加参数的无约束优化示例

接下来,将额外参数作为附加参量传递给目标函数,首先使用 MATLAB 文件,然后使用嵌套函数。

以上一示例中的目标函数为例。

f(x,y)=xexp(-(x2+y2))+(x2+y2)/20.

用 (a,b,c) 对该函数进行参数化,如下所示:

f(x,y,a,b,c)=(x-a)exp(-((x-a)2+(y-b)2))+((x-a)2+(y-b)2)/c.

此函数是原始目标函数的移位和缩放版本。

MATLAB 文件函数

假设有一个名为 bowlpeakfun 的 MATLAB 文件目标函数,定义如下。

type bowlpeakfun
function y = bowlpeakfun(x, a, b, c)
%BOWLPEAKFUN Objective function for parameter passing in TUTDEMO.

%   Copyright 2008 The MathWorks, Inc.

y = (x(1)-a).*exp(-((x(1)-a).^2+(x(2)-b).^2))+((x(1)-a).^2+(x(2)-b).^2)/c;

定义参数。

a = 2;
b = 3;
c = 10;

创建 MATLAB 文件的匿名函数句柄。

f = @(x)bowlpeakfun(x,a,b,c)
f = function_handle with value:
    @(x)bowlpeakfun(x,a,b,c)

调用 fminunc 以求最小值。

x0 = [-.5; 0];
options = optimoptions('fminunc','Algorithm','quasi-newton');
[x, fval] = fminunc(f,x0,options)
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
x = 2×1

    1.3639
    3.0000

fval = -0.3840

嵌套函数

假设有 nestedbowlpeak 函数,它将目标实现为一个嵌套函数。

type nestedbowlpeak
function [x,fval] =  nestedbowlpeak(a,b,c,x0,options)
%NESTEDBOWLPEAK Nested function for parameter passing in TUTDEMO.

%   Copyright 2008 The MathWorks, Inc.

[x,fval] = fminunc(@nestedfun,x0,options);  
    function y = nestedfun(x)
      y = (x(1)-a).*exp(-((x(1)-a).^2+(x(2)-b).^2))+((x(1)-a).^2+(x(2)-b).^2)/c;    
    end
end

参数 (a,b,c) 对于嵌套目标函数 nestedfun 可见。外部函数 nestedbowlpeak 调用 fminunc 并传递目标函数 nestedfun

定义参数、初始估计值和选项:

a = 2;
b = 3;
c = 10;
x0 = [-.5; 0];
options = optimoptions('fminunc','Algorithm','quasi-newton');

运行优化:

[x,fval] =  nestedbowlpeak(a,b,c,x0,options)
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
x = 2×1

    1.3639
    3.0000

fval = -0.3840

这两种方法产生相同的答案,因此您可以使用您认为最方便的一种。

约束优化示例:不相等性

还是上面那个问题,现在我们给它添加一个约束条件:

minimize xexp(-(x2+y2))+(x2+y2)/20,

subject to xy/2+(x+2)2+(y-2)2/22.

约束集是一个倾斜椭圆的内部。查看与倾斜椭圆一起绘制的目标函数的轮廓。

f = @(x,y) x.*exp(-x.^2-y.^2)+(x.^2+y.^2)/20;
g = @(x,y) x.*y/2+(x+2).^2+(y-2).^2/2-2;
fimplicit(g)
axis([-6 0 -1 7])
hold on
fcontour(f)
plot(-.9727,.4685,'ro');
legend('constraint','f contours','minimum');
hold off

绘图显示椭圆内目标函数的最小值出现在椭圆的右下角附近。在计算绘制的最小值之前,猜测解出现在以下位置。

x0 = [-2 1];

设置优化选项以使用内点算法,并显示每次迭代的结果。

options = optimoptions('fmincon','Algorithm','interior-point','Display','iter');

求解器要求非线性约束函数给出两个输出:一个用于非线性不等式,一个用于非线性等式。为了给出这两个输出,我们使用 deal 函数来编写约束。

gfun = @(x) deal(g(x(1),x(2)),[]);

调用非线性约束求解器。问题没有线性等式、不等式或边界,因此对这些参量传递 [ ]。

[x,fval,exitflag,output] = fmincon(fun,x0,[],[],[],[],[],[],gfun,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       3    2.365241e-01    0.000e+00    1.972e-01
    1       6    1.748504e-01    0.000e+00    1.734e-01    2.260e-01
    2      10   -1.570560e-01    0.000e+00    2.608e-01    9.347e-01
    3      14   -6.629160e-02    0.000e+00    1.241e-01    3.103e-01
    4      17   -1.584082e-01    0.000e+00    7.934e-02    1.826e-01
    5      20   -2.349124e-01    0.000e+00    1.912e-02    1.571e-01
    6      23   -2.255299e-01    0.000e+00    1.955e-02    1.993e-02
    7      26   -2.444225e-01    0.000e+00    4.293e-03    3.821e-02
    8      29   -2.446931e-01    0.000e+00    8.100e-04    4.035e-03
    9      32   -2.446933e-01    0.000e+00    1.999e-04    8.126e-04
   10      35   -2.448531e-01    0.000e+00    4.004e-05    3.289e-04
   11      38   -2.448927e-01    0.000e+00    4.036e-07    8.156e-05

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

显示求解器找到的解。

x
x = 1×2

   -0.9727    0.4686

查看在该解处的函数值。

fval
fval = -0.2449

查看函数计算的总次数。

Fevals = output.funcCount
Fevals = 38

在解处满足不等式约束。

[c, ceq] = gfun(x)
c = -2.4608e-06
ceq =

     []

由于 c(x) 接近 0,约束是活动的,这意味着它影响解。回想一下无约束解。

uncx
uncx = 2×1

   -0.6691
    0.0000

回想一下无约束目标函数。

uncf
uncf = -0.4052

查看约束使解发生了多大的移动并使目标函数的值增加了多少。

fval-uncf
ans = 0.1603

约束优化示例:用户提供的梯度

通过提供梯度,可以更高效、更准确地求解优化问题。与前一个示例一样,此示例求解不等式约束问题

minimize xexp(-(x2+y2))+(x2+y2)/20,

subject to xy/2+(x+2)2+(y-2)2/22.

为了给 fmincon 提供 f(x) 的梯度,以 MATLAB 文件的形式编写目标函数。

type onehump
function [f,gf] = onehump(x)
% ONEHUMP Helper function for Tutorial for the Optimization Toolbox demo

%   Copyright 2008-2009 The MathWorks, Inc.

r = x(1)^2 + x(2)^2;
s = exp(-r);
f = x(1)*s+r/20;

if nargout > 1
   gf = [(1-2*x(1)^2)*s+x(1)/10;
       -2*x(1)*x(2)*s+x(2)/10];
end

约束及其梯度包含在 MATLAB 文件 tiltellipse 中。

type tiltellipse
function [c,ceq,gc,gceq] = tiltellipse(x)
% TILTELLIPSE Helper function for Tutorial for the Optimization Toolbox demo

%   Copyright 2008-2009 The MathWorks, Inc.

c = x(1)*x(2)/2 + (x(1)+2)^2 + (x(2)-2)^2/2 - 2;
ceq = [];

if nargout > 2
   gc = [x(2)/2+2*(x(1)+2);
       x(1)/2+x(2)-2];
   gceq = [];
end

为找到解设置一个初始点。

x0 = [-2; 1];

设置优化选项,以使用与上例相同的算法进行比较。

options = optimoptions('fmincon','Algorithm','interior-point');

设置在目标函数和约束函数中使用梯度信息的选项。注意:必须启用这些选项,否则梯度信息将被忽略。

options = optimoptions(options,...
    'SpecifyObjectiveGradient',true,...
    'SpecifyConstraintGradient',true);

由于 fmincon 不需要使用有限差分来估计梯度,因此求解器应具有更少的函数计数。设置选项以在每次迭代时显示结果。

options.Display = 'iter';

调用求解器。

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    2.365241e-01    0.000e+00    1.972e-01
    1       2    1.748504e-01    0.000e+00    1.734e-01    2.260e-01
    2       4   -1.570560e-01    0.000e+00    2.608e-01    9.347e-01
    3       6   -6.629161e-02    0.000e+00    1.241e-01    3.103e-01
    4       7   -1.584082e-01    0.000e+00    7.934e-02    1.826e-01
    5       8   -2.349124e-01    0.000e+00    1.912e-02    1.571e-01
    6       9   -2.255299e-01    0.000e+00    1.955e-02    1.993e-02
    7      10   -2.444225e-01    0.000e+00    4.293e-03    3.821e-02
    8      11   -2.446931e-01    0.000e+00    8.100e-04    4.035e-03
    9      12   -2.446933e-01    0.000e+00    1.999e-04    8.126e-04
   10      13   -2.448531e-01    0.000e+00    4.004e-05    3.289e-04
   11      14   -2.448927e-01    0.000e+00    4.036e-07    8.156e-05

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

fmincon 在前面的示例中对梯度进行了很好的估计,因此在此示例中的迭代与之相似。

显示求解器找到的解。

xold = x
xold = 2×1

   -0.9727
    0.4686

查看在该解处的函数值。

minfval = fval
minfval = -0.2449

查看函数计算的总次数。

Fgradevals = output.funcCount
Fgradevals = 14

将此数字与无梯度的函数计算次数进行比较。

Fevals
Fevals = 38

约束优化示例:更改默认终止容差

此示例继续使用梯度,并求解同一约束问题

minimize xexp(-(x2+y2))+(x2+y2)/20,

subject to xy/2+(x+2)2+(y-2)2/22.

在本例中,您可以通过覆盖默认终止条件(options.StepToleranceoptions.OptimalityTolerance)来获得更准确的解。fmincon 内点算法的默认值为 options.StepTolerance = 1e-10options.OptimalityTolerance = 1e-6

覆盖这两个默认终止条件。

options = optimoptions(options,...
    'StepTolerance',1e-15,...
    'OptimalityTolerance',1e-8);

调用求解器。

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    2.365241e-01    0.000e+00    1.972e-01
    1       2    1.748504e-01    0.000e+00    1.734e-01    2.260e-01
    2       4   -1.570560e-01    0.000e+00    2.608e-01    9.347e-01
    3       6   -6.629161e-02    0.000e+00    1.241e-01    3.103e-01
    4       7   -1.584082e-01    0.000e+00    7.934e-02    1.826e-01
    5       8   -2.349124e-01    0.000e+00    1.912e-02    1.571e-01
    6       9   -2.255299e-01    0.000e+00    1.955e-02    1.993e-02
    7      10   -2.444225e-01    0.000e+00    4.293e-03    3.821e-02
    8      11   -2.446931e-01    0.000e+00    8.100e-04    4.035e-03
    9      12   -2.446933e-01    0.000e+00    1.999e-04    8.126e-04
   10      13   -2.448531e-01    0.000e+00    4.004e-05    3.289e-04
   11      14   -2.448927e-01    0.000e+00    4.036e-07    8.156e-05
   12      15   -2.448931e-01    0.000e+00    4.000e-09    8.230e-07

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

为了更准确地查看新容差造成的差异,在解中显示更多小数位数。

format long

显示求解器找到的解。

x
x = 2×1

  -0.972742227363546
   0.468569289098342

将这些值与前面示例中的值进行比较。

xold
xold = 2×1

  -0.972742694488360
   0.468569966693330

确定值的变化。

x - xold
ans = 2×1
10-6 ×

   0.467124813385844
  -0.677594988729435

查看在该解处的函数值。

fval
fval = 
  -0.244893137879894

查看解改进了多少。

fval - minfval
ans = 
    -3.996450220755676e-07

答案为负值,因为新解更小。

查看函数计算的总次数。

output.funcCount
ans = 
    15

将此数字与使用用户提供的梯度和默认容差求解的示例中的函数计算数进行比较。

Fgradevals
Fgradevals = 
    14

约束优化示例:用户提供的黑塞矩阵

如果除了梯度之外,您还提供黑塞矩阵,求解器会更加精确和高效。

fmincon 的内点算法将黑塞矩阵作为一个单独的函数(不是目标函数的一部分)。黑塞函数 H(x,lambda) 计算拉格朗日的黑塞矩阵;请参阅 适用于 fmincon 内点算法的黑塞函数

求解器计算值 lambda.ineqnonlinlambda.eqlin;您的黑塞函数会告知求解器如何使用这些值。

此示例有一个不等式约束,因此黑塞矩阵在 hessfordemo 函数中给出。

type hessfordemo
function H = hessfordemo(x,lambda)
% HESSFORDEMO Helper function for Tutorial for the Optimization Toolbox demo

%   Copyright 2008-2009 The MathWorks, Inc.

s = exp(-(x(1)^2+x(2)^2));
H = [2*x(1)*(2*x(1)^2-3)*s+1/10, 2*x(2)*(2*x(1)^2-1)*s;
        2*x(2)*(2*x(1)^2-1)*s, 2*x(1)*(2*x(2)^2-1)*s+1/10];
hessc = [2,1/2;1/2,1];
H = H + lambda.ineqnonlin(1)*hessc;

为了使用黑塞矩阵,您需要适当地设置选项。

options = optimoptions('fmincon',...
    'Algorithm','interior-point',...
    'SpecifyConstraintGradient',true,...
    'SpecifyObjectiveGradient',true,...
    'HessianFcn',@hessfordemo);

容差设置为默认值,这样应得到更少的函数计数。设置选项以在每次迭代时显示结果。

options.Display = 'iter';

调用求解器。

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    2.365241e-01    0.000e+00    1.972e-01
    1       3    5.821325e-02    0.000e+00    1.443e-01    8.728e-01
    2       5   -1.218829e-01    0.000e+00    1.007e-01    4.927e-01
    3       6   -1.421167e-01    0.000e+00    8.486e-02    5.165e-02
    4       7   -2.261916e-01    0.000e+00    1.989e-02    1.667e-01
    5       8   -2.433609e-01    0.000e+00    1.537e-03    3.486e-02
    6       9   -2.446875e-01    0.000e+00    2.057e-04    2.727e-03
    7      10   -2.448911e-01    0.000e+00    2.068e-06    4.191e-04
    8      11   -2.448931e-01    0.000e+00    2.001e-08    4.218e-06

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

结果显示迭代有变化且次数更少。

显示求解器找到的解。

x
x = 2×1

  -0.972742246093537
   0.468569316215571

查看在该解处的函数值。

fval
fval = 
  -0.244893121872758

查看函数计算的总次数。

output.funcCount
ans = 
    11

将此数字与仅使用梯度计算和相同默认容差求解的示例中的函数计算次数进行比较。

Fgradevals
Fgradevals = 
    14

相关主题