使用线性等式约束的信赖域反射算法进行最小化
fmincon
trust-region-reflective
算法可以最小化仅受线性等式约束(无边界或任何其他约束)的非线性目标函数。例如,最小化
受到一些线性等式约束。此示例采用 。
创建问题
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)
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