Main Content

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

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

f(x)=i=1n-1((xi2)(xi+12+1)+(xi+12)(xi2+1)),

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

创建问题

browneq.mat 文件包含矩阵 Aeqbeq,它们代表线性约束 Aeq*x = beqAeq 矩阵有 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)
     3
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

相关主题