具有边界约束和带状预条件子的最小化

此示例显示如何使用 fmincon trust-region-reflective 算法求解边界的非线性问题。当问题稀疏且具有解析梯度和已知结构（例如其 Hessian 模式）时，该算法提供了额外的效率。

带梯度的目标函数

对于给定的 $n$ （4 的正倍数），目标函数为

$f (x) = 1 + \sum_{i = 1}^{n} {| (3 - 2 x_{i}) x_{i} - x_{i - 1} - x_{i + 1} + 1 |}^{p} + \sum_{i = 1}^{n / 2} {| x_{i} + x_{i + n / 2} |}^{p},$

其中 $p = 7 / 3$ 、 $x_{0} = 0$ 且 $x_{n + 1} = 0$ 。本示例末尾的 tbroyfg 辅助函数实现了目标函数，包括其梯度。

该问题对所有 $i$ 都有边界 $- 10 \leq x_{i} \leq 10$ 。使用 $n$ = 800。

n = 800;
lb = -10*ones(n,1);
ub = -lb;

Hessian 模式

Hessian 矩阵的稀疏模式是预先确定的，并存储在文件 tbroyhstr.mat 中。正如您在下面的间谍图中所看到的，该问题的 Hessian 的稀疏结构是带状的。

load tbroyhstr
spy(Hstr)

在此图中，中心条纹本身是一个五带矩阵。下图更清晰地显示了矩阵。

spy(Hstr(1:20,1:20))

问题选项

设置选项以使用 trust-region-reflective 算法。该算法要求您将 SpecifyObjectiveGradient 选项设置为 true。

另外，使用 optimoptions 将 HessPattern 选项设置为 Hstr。如果您没有为具有明显稀疏结构的大型问题设置此选项，则该问题将使用大量内存和计算，因为 fmincon 尝试对 640,000 个非零条目的完整 Hessian 矩阵使用有限差分。

options = optimoptions('fmincon','SpecifyObjectiveGradient',true,'HessPattern',Hstr,...
    'Algorithm','trust-region-reflective');

求解问题

对于奇数索引，将初始点设置为 -1；对于偶数索引，将初始点设置为 +1。

x0 = -ones(n,1);
x0(2:2:n) = 1;

该问题没有线性或非线性约束，因此将这些参数设置为 []。

A = [];
b = [];
Aeq = [];
beq = [];
nonlcon = [];

调用 fmincon 以求解问题。

[x,fval,exitflag,output] = ... 
   fmincon(@tbroyfg,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

Local minimum possible.

fmincon stopped because the final change in function value relative to 
its initial value is less than the value of the function tolerance.

检查解和求解过程

检查退出标志、目标函数值、一阶最优性测度和求解器迭代次数。

disp(exitflag);

disp(fval)

  270.4790

disp(output.firstorderopt)

    0.0163

disp(output.iterations)

fmincon 不需要经过多次迭代就能得出解。但是该解具有相对较高的一阶最优性测度，这就是退出标志没有首选值 1 的原因。

改进解

尝试使用五带预条件子而不是默认的对角预条件子。使用 optimoptions，将 PrecondBandWidth 选项设置为 2 并再次求解问题。（带宽是上对角线或下对角线的数量，不计算主对角线。）

options.PrecondBandWidth = 2;
[x2,fval2,exitflag2,output2] = ... 
   fmincon(@tbroyfg,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

Local minimum possible.

fmincon stopped because the final change in function value relative to 
its initial value is less than the value of the function tolerance.

disp(exitflag2);

disp(fval2)

  270.4790

disp(output2.firstorderopt)

   7.5340e-05

disp(output2.iterations)

退出标志和目标函数值似乎没有改变。然而，随着迭代次数的增加，一阶最优性测度会大幅下降。计算目标函数值的差异。

disp(fval2 - fval)

  -2.9005e-07

目标函数值微小地减少。该解主要改进了一阶最优性测度，而不是目标函数。

辅助函数

以下代码创建 tbroyfg 辅助函数。

function [f,grad] = tbroyfg(x,dummy)
%TBROYFG Objective function and its gradients for nonlinear minimization
% with bound constraints and banded preconditioner.
% Documentation example.

%   Copyright 1990-2008 The MathWorks, Inc.

n = length(x);  % n should be a multiple of 4
p = 7/3;
y=zeros(n,1);
i = 2:(n-1);
y(i) = abs((3-2*x(i)).*x(i) - x(i-1) - x(i+1)+1).^p;
y(n) = abs((3-2*x(n)).*x(n) - x(n-1)+1).^p;
y(1) = abs((3-2*x(1)).*x(1) - x(2)+1).^p;
j = 1:(n/2);
z = zeros(length(j),1);
z(j) = abs(x(j) + x(j+n/2)).^p;
f = 1 + sum(y) + sum(z);
%
% Evaluate the gradient.
if nargout > 1
   g = zeros(n,1);
   t = zeros(n,1);
   i = 2:(n-1);
   t(i) = (3-2*x(i)).*x(i) - x(i-1) - x(i+1) + 1;
   g(i) = p*abs(t(i)).^(p-1).*sign(t(i)).*(3-4*x(i));
   g(i-1) = g(i-1) - p*abs(t(i)).^(p-1).*sign(t(i));
   g(i+1) = g(i+1) - p*abs(t(i)).^(p-1).*sign(t(i));
   tt = (3-2*x(n)).*x(n) - x(n-1) + 1;
   g(n) = g(n) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(n));
   g(n-1) = g(n-1) - p*abs(tt).^(p-1).*sign(tt);
   tt = (3-2*x(1)).*x(1)-x(2)+1;
   g(1) = g(1) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(1));
   g(2) = g(2) - p*abs(tt).^(p-1).*sign(tt);
   j = 1:(n/2);
   t(j) = x(j) + x(j+n/2);
   g(j) = g(j) + p*abs(t(j)).^(p-1).*sign(t(j));
   jj = j + (n/2);
   g(jj) = g(jj) + p*abs(t(j)).^(p-1).*sign(t(j));
   grad = g;
end
end