具有边界约束和带状预条件子的最小化
此示例显示如何使用 fmincon
trust-region-reflective
算法求解边界的非线性问题。当问题稀疏且具有解析梯度和已知结构(例如其 Hessian 模式)时,该算法提供了额外的效率。
带梯度的目标函数
对于给定的 (4 的正倍数),目标函数为
其中 、 且 。本示例末尾的 tbroyfg
辅助函数实现了目标函数,包括其梯度。
该问题对所有 都有边界 。使用 = 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);
3
disp(fval)
270.4790
disp(output.firstorderopt)
0.0163
disp(output.iterations)
7
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);
3
disp(fval2)
270.4790
disp(output2.firstorderopt)
7.5340e-05
disp(output2.iterations)
9
退出标志和目标函数值似乎没有改变。然而,随着迭代次数的增加,一阶最优性测度会大幅下降。计算目标函数值的差异。
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