Optimization Toolbox 教程
2024-04-07假设要求下面这个函数的最小值:
绘制函数以查看它具有最小值的位置。
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];
设置优化选项以使用 默认 算法。此步骤可确保教程在每个 MATLAB 版本中的工作方式是相同的。
options=optimoptions('fminunc','Algorithm','quasi-newton');
在求解器执行计算时查看迭代。
options.Display='iter';
调用 ,它是一个无约束非线性最小化求解函数。
[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 文件,然后使用嵌套函数。
以上一示例中的目标函数为例。
用 (a,b,c) 对该函数进行参数化,如下所示:
此函数是原始目标函数的移位和缩放版本。
MATLAB 文件函数
假设有一个名为 的 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)
调用 以求最小值。
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
嵌套函数
假设有 函数,它将目标实现为一个嵌套函数。
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) 对于嵌套目标函数 可见。外部函数 调用 并传递目标函数 。
定义参数、初始估计值和选项:
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
这两种方法产生相同的答案,因此您可以使用您认为最方便的一种。
还是上面那个问题,现在我们给它添加一个约束条件:
约束集是一个倾斜椭圆的内部。查看与倾斜椭圆一起绘制的目标函数的轮廓。
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');
求解器要求非线性约束函数给出两个输出:一个用于非线性不等式,一个用于非线性等式。为了给出这两个输出,我们使用 函数来编写约束。
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
通过提供梯度,可以更高效、更准确地求解优化问题。与前一个示例一样,此示例求解不等式约束问题
为了给 提供 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 文件 中。
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);
由于 不需要使用有限差分来估计梯度,因此求解器应具有更少的函数计数。设置选项以在每次迭代时显示结果。
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.
在前面的示例中对梯度进行了很好的估计,因此在此示例中的迭代与之相似。
显示求解器找到的解。
xold=x
xold=2×1
-0.9727
0.4686
查看在该解处的函数值。
minfval=fval
minfval=-0.2449
查看函数计算的总次数。
Fgradevals=output.funcCount
Fgradevals=14
将此数字与无梯度的函数计算次数进行比较。
Fevals
Fevals=38
此示例继续使用梯度,并求解同一约束问题
.
在本例中,您可以通过覆盖默认终止条件( 和 )来获得更准确的解。 内点算法的默认值为 和 。
覆盖这两个默认终止条件。
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
如果除了梯度之外,您还提供黑塞矩阵,求解器会更加精确和高效。
的内点算法将黑塞矩阵作为一个单独的函数(不是目标函数的一部分)。黑塞函数 H(x,lambda) 计算拉格朗日的黑塞矩阵;请参阅 适用于 fmincon 内点算法的黑塞函数。
求解器计算值 和 ;您的黑塞函数会告知求解器如何使用这些值。
此示例有一个不等式约束,因此黑塞矩阵在 函数中给出。
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