使用线性等式约束的信赖域反射算法进行最小化

fmincon trust-region-reflective 算法可以最小化仅受线性等式约束（无边界或任何其他约束）的非线性目标函数。例如，最小化

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

受到一些线性等式约束。此示例采用 $n = 1000$ 。

创建问题

browneq.mat 文件包含矩阵 Aeq 和 beq，它们代表线性约束 Aeq*x = beq。Aeq 矩阵有 100 行，代表 100 个线性约束（因此 Aeq 是一个 100×1000 的矩阵）。加载 browneq.mat 数据。

load browneq.mat

本示例末尾的 brownfgh 辅助函数实现了目标函数，包括其梯度和 Hessian。

设置选项

trust-region-reflective 算法要求目标函数包含梯度。该算法在目标函数中接受 Hessian。设置选项以包含所有衍生信息。

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

求解问题

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

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

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

A = [];
b = [];
lb = [];
ub = [];
nonlcon = [];

调用 fmincon 以求解问题。

[x,fval,exitflag,output] = ...
   fmincon(@brownfgh,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)

  205.9313

disp(output.constrviolation)

   2.2338e-13

exitflag 的值为 3 表示 fmincon 停止，因为目标函数值的变化小于容差 FunctionTolerance。最终的目标函数值由 fval 给出。约束得到满足，如 output.constrviolation 所示，显示一个非常小的数字。

要自己计算约束违反值，请执行以下代码。

norm(Aeq*x-beq,Inf)

ans = 2.2338e-13

辅助函数

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

function [f,g,H] = brownfgh(x)
%BROWNFGH  Nonlinear minimization problem (function, its gradients
% and Hessian).
% Documentation example.         

%   Copyright 1990-2019 The MathWorks, Inc.

% Evaluate the function.
  n = length(x);
  y = zeros(n,1);
  i = 1:(n-1);
  y(i) = (x(i).^2).^(x(i+1).^2+1)+(x(i+1).^2).^(x(i).^2+1);
  f = sum(y);

% Evaluate the gradient.
  if nargout > 1
     i=1:(n-1);
     g = zeros(n,1);
     g(i) = 2*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
             2*x(i).*((x(i+1).^2).^(x(i).^2+1)).*log(x(i+1).^2);
     g(i+1) = g(i+1)+...
              2*x(i+1).*((x(i).^2).^(x(i+1).^2+1)).*log(x(i).^2)+...
              2*(x(i).^2+1).*x(i+1).*((x(i+1).^2).^(x(i).^2));
  end

% Evaluate the (sparse, symmetric) Hessian matrix
  if nargout > 2
     v = zeros(n,1);
     i = 1:(n-1);
     v(i) = 2*(x(i+1).^2+1).*((x(i).^2).^(x(i+1).^2))+...
            4*(x(i+1).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i).^2).^((x(i+1).^2)-1))+...
            2*((x(i+1).^2).^(x(i).^2+1)).*(log(x(i+1).^2));
     v(i) = v(i)+4*(x(i).^2).*((x(i+1).^2).^(x(i).^2+1)).*((log(x(i+1).^2)).^2);
     v(i+1) = v(i+1)+...
              2*(x(i).^2).^(x(i+1).^2+1).*(log(x(i).^2))+...
              4*(x(i+1).^2).*((x(i).^2).^(x(i+1).^2+1)).*((log(x(i).^2)).^2)+...
              2*(x(i).^2+1).*((x(i+1).^2).^(x(i).^2));
     v(i+1) = v(i+1)+4*(x(i).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i+1).^2).^(x(i).^2-1));
     v0 = v;
     v = zeros(n-1,1);
     v(i) = 4*x(i+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
            4*x(i+1).*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2)).*log(x(i).^2);
     v(i) = v(i)+ 4*x(i+1).*x(i).*((x(i+1).^2).^(x(i).^2)).*log(x(i+1).^2);
     v(i) = v(i)+4*x(i).*((x(i+1).^2).^(x(i).^2)).*x(i+1);
     v1 = v;
     i = [(1:n)';(1:(n-1))'];
     j = [(1:n)';(2:n)'];
     s = [v0;2*v1];
     H = sparse(i,j,s,n,n);
     H = (H+H')/2;
  end
end