基于问题的约束约束非线性最小化

此示例说明如何使用基于问题的方法最小化受线性等式约束的非线性函数，其中您可以根据优化变量来构造约束。此示例还展示了如何使用 fcn2optimexpr 将目标函数文件转换为优化表达式。

示例使用线性等式约束的信赖域反射算法进行最小化使用涉及梯度和 Hessian 的基于求解器的方法。使用基于问题的方法求解同样的问题很简单，但需要更多的解时间，因为当问题有不受支持的运算符时，基于问题的方法目前不使用梯度信息。此外，基于问题的方法不使用 Hessian 信息。

创建问题和目标

问题可以表示为最小化下式

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

受到一组线性等式约束 Aeq*x = beq。首先创建优化问题和变量。

prob = optimproblem;
N = 1000;
x = optimvar('x',N);

目标函数位于 brownfgh.m 文件中，运行此示例时该文件可用。您必须使用 fcn2optimexpr 将函数转换为优化表达式，因为优化变量不会出现在指数中。请参阅优化变量和表达式支持的运算和将非线性函数转换为优化表达式。

prob.Objective = fcn2optimexpr(@brownfgh,x,'OutputSize',[1,1]);

包括约束

要获取工作区中的 Aeq 和 beq 矩阵，请执行此命令。

load browneq

在问题中包括线性约束。

prob.Constraints = Aeq*x == beq;

回顾并求解问题

回顾问题目标。

type brownfgh

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

%   Copyright 1990-2008 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

该问题有一百个线性等式约束，因此得到的约束表达式太长，无法包含在示例中。为了显示约束，请取消注释并运行以下行。

% show(prob.Constraints)

将初始点设置为具有字段 x 的结构。

x0.x = -ones(N,1);
x0.x(2:2:N) = 1;

通过调用 solve 求解问题。

[sol,fval,exitflag,output] = solve(prob,x0)

Solving problem using fmincon.

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+03.

sol = struct with fields:
    x: [1000x1 double]

fval = 207.5463

exitflag = 
    SolverLimitExceeded

output = struct with fields:
              iterations: 2
               funcCount: 3007
         constrviolation: 3.5194e-13
                stepsize: 1.9303
               algorithm: 'interior-point'
           firstorderopt: 2.6432
            cgiterations: 0
                 message: 'Solver stopped prematurely....'
            bestfeasible: [1x1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

由于超出了函数评估限制，求解器过早停止。要继续优化，请从终点重新开始优化，并允许进行更多的函数计算。

options = optimoptions(prob,'MaxFunctionEvaluations',1e5);
[sol,fval,exitflag,output] = solve(prob,sol,'Options',options)

Solving problem using fmincon.

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.

sol = struct with fields:
    x: [1000x1 double]

fval = 205.9313

exitflag = 
    OptimalSolution

output = struct with fields:
              iterations: 40
               funcCount: 41086
         constrviolation: 2.8422e-14
                stepsize: 5.7844e-06
               algorithm: 'interior-point'
           firstorderopt: 4.1038e-06
            cgiterations: 4
                 message: 'Local minimum found that satisfies the constraints....'
            bestfeasible: [1x1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

与基于求解器的解进行比较

要使用基于求解器的方法求解问题（如使用线性等式约束的信赖域反射算法进行最小化所示），请将初始点转换为向量。然后设置选项以使用 brownfgh 中提供的梯度和 Hessian 信息。

xstart = x0.x;
fun = @brownfgh;
opts = optimoptions('fmincon','SpecifyObjectiveGradient',true,'HessianFcn','objective',...
    'Algorithm','trust-region-reflective');
[x,fval,exitflag,output] = ...
   fmincon(fun,xstart,[],[],Aeq,beq,[],[],[],opts);

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.

fprintf("Fval = %g\nNumber of iterations = %g\nNumber of function evals = %g.\n",...
    fval,output.iterations,output.funcCount)

Fval = 205.931
Number of iterations = 22
Number of function evals = 23.

使用线性等式约束的信赖域反射算法进行最小化中的基于求解器的解使用目标函数中提供的梯度和 Hessian。通过使用该导数信息，求解器 fmincon 仅使用 23 个函数计算，便可在 22 次迭代中收敛到解。基于求解器的解与基于问题的解具有相同的最终目标函数值。

然而，如果不使用符号数学来构建梯度和 Hessian 函数很困难并且容易出错。有关如何使用符号数学计算导数的示例，请参阅使用 Symbolic Math Toolbox 计算梯度和 Hessian。

另请参阅

fcn2optimexpr