本页对应的英文页面已更新,但尚未翻译。 若要查看最新内容,请点击此处访问英文页面。

Optimization Toolbox™ 教程

此示例说明如何使用两个非线性优化求解器以及如何设置选项。我们在此示例中使用的非线性求解器是 fminuncfmincon

此示例中概述的所有原则都适用于其他非线性求解器,如 fgoalattainfminimaxlsqnonlinlsqcurvefitfsolve

该示例从最小化目标函数开始,进而扩展到使用附加参数最小化同一函数。之后,该示例说明如何在约束条件下最小化目标函数,并最终展示如何通过提供梯度或 Hessian 矩阵,或通过更改一些选项来获得更高效和/或更准确的解。

无约束优化示例

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

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 的默认大规模算法,因为该算法要求提供目标函数的梯度:

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.

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

方法 1: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

方法 2:嵌套函数

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

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.StepTolerance and options.OptimalityTolerance) 使求解更准确。我们继续使用梯度。fmincon 的内点算法的默认值为 options.StepTolerance = 1e-10,options.OptimalityTolerance = 1e-6。

覆盖两个默认终止条件:X 和 fval 的终止容差。

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

约束优化示例:使用用户提供的 Hessian 矩阵

如果您不仅给出梯度,还给出 Hessian 矩阵,求解器会更加精确和高效。

fmincon 的内点求解器将 Hessian 矩阵作为一个单独的函数(不是目标函数的一部分)。Hessian 函数 H(x,lambda) 应计算拉格朗日的 Hessian 矩阵;关于其定义,请参阅《用户指南》。

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

在此问题中,我们只有一个不等式约束,因此 Hessian 矩阵为:

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;

为了使用 Hessian 矩阵,您需要适当地设置选项:

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

相关主题