В этом примере показано, как решить задачу оптимизации, имеющую линейные или квадратичные целевые и квадратичные ограничения неравенства. В примере генерируются и используются градиент и гессен целевой функции и функции ограничения.
Предположим, что ваша проблема имеет форму
+ c)
подлежит
где 1 ≤ i ≤ m Предположим, что по меньшей мере один Hi является ненулевым; в противном случае можно использовать quadprog или linprog для решения этой проблемы. При ненулевом Hi ограничения нелинейны, что означает fmincon является подходящим решателем в соответствии с таблицей решений оптимизации.
В примере предполагается, что квадратичные матрицы симметричны без потери общности. Несимметричную матрицу H (или Q) можно заменить эквивалентной симметричной версией (H + HT )/2.
Если x имеет N компонентов, то Q и Hi являются N-на-N матрицами, f и ki являются N-by-1 векторами, а c и di являются скалярами.
Сформулировать проблему с помощью fmincon синтаксис. Предположим, что x и f являются векторами столбцов. (x - вектор столбца, когда начальный вектор x0 является вектором столбца.)
function [y,grady] = quadobj(x,Q,f,c) y = 1/2*x'*Q*x + f'*x + c; if nargout > 1 grady = Q*x + f; end
Для обеспечения согласованности и простоты индексирования поместите каждую квадратную матрицу ограничений в один массив ячеек. Аналогично, поместите линейные и постоянные члены в массивы ячеек.
function [y,yeq,grady,gradyeq] = quadconstr(x,H,k,d) jj = length(H); % jj is the number of inequality constraints y = zeros(1,jj); for i = 1:jj y(i) = 1/2*x'*H{i}*x + k{i}'*x + d{i}; end yeq = []; if nargout > 2 grady = zeros(length(x),jj); for i = 1:jj grady(:,i) = H{i}*x + k{i}; end end gradyeq = [];
Предположим, что у вас есть следующая проблема.
Q = [3,2,1;
2,4,0;
1,0,5];
f = [-24;-48;-130];
c = -2;
rng default % For reproducibility
% Two sets of random quadratic constraints:
H{1} = gallery('randcorr',3); % Random positive definite matrix
H{2} = gallery('randcorr',3);
k{1} = randn(3,1);
k{2} = randn(3,1);
d{1} = randn;
d{2} = randn;Создайте функцию Гессена. Гессен лагранжиана задаётся уравнением
+∑λi∇2ceqi (x).
fmincon вычисляет приблизительное множество множителей Лагранжа λ i и упаковывает их в структуру. Для включения гессенского используйте следующую функцию.
function hess = quadhess(x,lambda,Q,H) hess = Q; jj = length(H); % jj is the number of inequality constraints for i = 1:jj hess = hess + lambda.ineqnonlin(i)*H{i}; end
Используйте fmincon interior-point алгоритм наиболее эффективного решения проблемы. Этот алгоритм принимает функцию Гессена, которую вы предоставляете. Задайте эти параметры.
options = optimoptions(@fmincon,'Algorithm','interior-point',... 'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,... 'HessianFcn',@(x,lambda)quadhess(x,lambda,Q,H));
Звонить fmincon для решения проблемы.
fun = @(x)quadobj(x,Q,f,c); nonlconstr = @(x)quadconstr(x,H,k,d); x0 = [0;0;0]; % Column vector [x,fval,eflag,output,lambda] = fmincon(fun,x0,... [],[],[],[],[],[],nonlconstr,options);
Проверьте множители Лагранжа.
lambda.ineqnonlin
ans = 12.8412 39.2337
Оба нелинейных множителя неравенства ненулевые, поэтому оба квадратичных ограничения активны в решении.
Алгоритм внутренней точки с градиентами и гессенским является эффективным. Просмотр количества аналитических отчетов по функциям.
output
output =
iterations: 9
funcCount: 10
constrviolation: 0
stepsize: 5.3547e-04
algorithm: 'interior-point'
firstorderopt: 1.5851e-05
cgiterations: 0
message: 'Local minimum found that satisfies the constraints.
Optimization compl...'fmincon для решения проблемы требуется всего 10 оценок функций.
Сравните этот результат с решением без гессенского.
options.HessianFcn = [];
[x2,fval2,eflag2,output2,lambda2] = fmincon(fun,[0;0;0],...
[],[],[],[],[],[],nonlconstr,options);
output2output2 =
iterations: 17
funcCount: 22
constrviolation: 0
stepsize: 2.8475e-04
algorithm: 'interior-point'
firstorderopt: 1.7680e-05
cgiterations: 0
message: 'Local minimum found that satisfies the constraints.
Optimization compl...'В этом случае fmincon занимает примерно в два раза больше итераций и оценок функций. Решения те же, что и в пределах допусков.
Если также имеются квадратичные ограничения равенства, можно использовать по существу тот же метод. Проблема та же, с дополнительными ограничениями
qi = 0.
Переформулируйте ограничения, чтобы использовать переменные Ji, pi и qi. lambda.eqnonlin(i) структура имеет множители Лагранжа для ограничений равенства.