В этом примере показано, как включить информацию о производных в нелинейную оптимизацию на основе задач, когда автоматические производные не применяются, или когда требуется включить гессен, который не предоставляется с помощью автоматического дифференцирования. Включение градиентов или гессена в оптимизацию может дать следующие преимущества:
Более надежные результаты. Шаги конечных разностей иногда достигают точек, где целевая или нелинейная функция ограничения не определена, не является конечной или сложной.
Повышенная точность. Аналитические градиенты могут быть более точными, чем оценки конечных разностей.
Более быстрая сходимость. Включение гессена может привести к более быстрой сходимости, что означает меньшее количество итераций.
Улучшенная производительность. Аналитические градиенты могут быть быстрее для вычисления, чем оценки конечных разностей, особенно для проблем с разреженной структурой. Однако для сложных выражений расчет аналитических градиентов может быть медленнее.
Начиная с R2020b, solve функция может использовать автоматическое дифференцирование для поддерживаемых функций, чтобы предоставлять градиенты решателям. Эти автоматические производные применяются только к градиентам (первые производные), а не гессенам (вторые производные).
Автоматическое дифференцирование применяется, когда вы не используете fcn2optimexpr для создания целевой функции или функции ограничения. Если вам нужно использовать fcn2optimexpr, в этом примере показано, как включить производную информацию.
Способ использования производных в оптимизации на основе задач без автоматической дифференциации состоит в преобразовании проблемы с помощью prob2struct, а затем отредактируйте результирующие функции цели и ограничения. В этом примере показан гибридный подход, в котором автоматическая дифференциация обеспечивает производные для части целевой функции.
С управляющими переменными x и y, используйте целевую функцию
цель = fun1 + fun2.
Включить ограничение, что сумма квадратов x и y не более 4.
fun2 не основан на поддерживаемых функциях для выражений оптимизации; см. раздел Поддерживаемые операции с переменными и выражениями оптимизации. Следовательно, включить fun2 в задаче оптимизации необходимо преобразовать ее в выражение оптимизации с помощью fcn2optimexpr.
Чтобы использовать AD для поддерживаемых функций, настройте проблему без неподдерживаемой функции fun2, и включить fun2 позже.
prob = optimproblem; x = optimvar('x'); y = optimvar('y'); fun1 = 100*(y - x^2)^2 + (1 - x)^2; prob.Objective = fun1; prob.Constraints.cons = x^2 + y^2 <= 4;
Включить производные от fun2, сначала преобразовать проблему без fun2 в структуру с использованием prob2struct.
problem = prob2struct(prob,... 'ObjectiveFunctionName','generatedObjectiveBefore');
Во время преобразования, prob2struct создает файлы функций, представляющие целевую и нелинейную функции ограничения. По умолчанию эти функции имеют имена 'generatedObjective.m' и 'generatedConstraints.m'. Файл целевой функции без fun2 является 'generatedObjectiveBefore.m'.
Созданные функции цели и ограничения включают градиенты.
Вычисление производных, связанных с fun2. При наличии лицензии на символьные математические Toolbox™ можно использовать gradient (инструментарий символьной математики) или jacobian (Символьная математическая панель инструментов), помогающая вычислить производные. См. раздел Расчет градиентов и гессенов с помощью символьных математических Toolbox™.
Подход на основе решателя имеет одну управляющую переменную. Каждая переменная оптимизации (x и y, в этом примере) является частью управляющей переменной. Отображение между переменными оптимизации и одной управляющей переменной можно найти в созданном файле целевой функции. 'generatedObjectiveBefore.m'. Начало файла содержит эти строки кода или подобные строки.
%% Variable indices. xidx = 1; yidx = 2; %% Map solver-based variables to problem-based. x = inputVariables(xidx); y = inputVariables(yidx);
В этом коде управляющая переменная имеет имя inputVariables.
Можно также найти сопоставление с помощью varindex.
idx = varindex(prob); disp(idx.x)
1
disp(idx.y)
2
Полная целевая функция включает fun2:
fun2 = besselj(1,x^2 + y^2);
Используя стандартное исчисление, рассчитайте gradfun2, градиент fun2.
бессель (1, x2 + y2 )/( x2 + y2))].
Править 'generatedObjectiveBefore.m' включать fun2.
%% Compute objective function. arg1 = (y - x.^2); arg2 = 100; arg3 = arg1.^2; arg4 = (1 - x); obj = ((arg2 .* arg3) + arg4.^2); ssq = x^2 + y^2; fun2 = besselj(1,ssq); obj = obj + fun2;
Включение расчетных градиентов в файл целевой функции путем редактирования 'generatedObjectiveBefore.m'. Если имеется версия программного обеспечения, которая не выполняет расчет градиента, включите все эти строки. Если программа выполняет расчет градиента, включите только жирные строки, которые вычисляют градиент fun2.
%% Compute objective gradient.
if nargout > 1
arg5 = 1;
arg6 = zeros([2, 1]);
arg6(xidx,:) = (-(arg5.*2.*(arg4(:)))) + ((-((arg5.*arg2(:)).*2.*(arg1(:)))).*2.*(x(:)));
arg6(yidx,:) = ((arg5.*arg2(:)).*2.*(arg1(:)));
grad = arg6(:);
arg7 = besselj(0,ssq);
arg8 = arg7 - fun2/ssq;
gfun = [2*x*arg8;...
2*y*arg8];
grad = grad + gfun;
endНапомним, что нелинейным ограничением является x^2 + y^2 <= 4. Градиент этой функции ограничения равен 2*[x;y]. Если программа вычисляет градиент зависимости и включает его в созданный файл зависимости, то больше ничего делать не нужно. Если программа не вычисляет градиент ограничения, включите градиент нелинейного ограничения путем редактирования 'generatedConstraints.m' для включения этих строк.
%% Insert gradient calculation here.
% If you include a gradient, notify the solver by setting the
% SpecifyConstraintGradient option to true.
if nargout > 2
cineqGrad = 2*[x;y];
ceqGrad = [];
endЗапустите задачу, используя методы на основе решателя и на основе задачи (без градиента), чтобы увидеть различия. Чтобы запустить задачу на основе решателя с использованием информации о производной, создайте соответствующие опции в структуре задачи.
options = optimoptions('fmincon','SpecifyObjectiveGradient',true,... 'SpecifyConstraintGradient',true); problem.options = options;
Нелинейные проблемы требуют непустого x0 поле в структуре проблемы.
x0 = [1;1]; problem.x0 = x0;
Звонить fmincon о структуре проблемы.
[xsolver,fvalsolver,exitflagsolver,outputsolver] = fmincon(problem)
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.
<stopping criteria details>
xsolver =
1.2494
1.5617
fvalsolver =
-0.0038
exitflagsolver =
1
outputsolver =
struct with fields:
iterations: 15
funcCount: 32
constrviolation: 0
stepsize: 1.5569e-06
algorithm: 'interior-point'
firstorderopt: 2.2058e-08
cgiterations: 7
message: '↵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.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 2.125635e-08,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
bestfeasible: [1×1 struct]Сравните это решение с решением, полученным из solve без производной информации.
init.x = x0(1); init.y = x0(2); f2 = @(x,y)besselj(1,x^2 + y^2); fun2 = fcn2optimexpr(f2,x,y); prob.Objective = prob.Objective + fun2; [xproblem,fvalproblem,exitflagproblem,outputproblem] = solve(prob,init,... "ConstraintDerivative","finite-differences",... "ObjectiveDerivative","finite-differences")
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.
<stopping criteria details>
xproblem =
struct with fields:
x: 1.2494
y: 1.5617
fvalproblem =
-0.0038
exitflagproblem =
OptimalSolution
outputproblem =
struct with fields:
iterations: 15
funcCount: 64
constrviolation: 0
stepsize: 1.5571e-06
algorithm: 'interior-point'
firstorderopt: 6.0139e-07
cgiterations: 7
message: '↵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.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 5.795368e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
bestfeasible: [1×1 struct]
objectivederivative: "finite-differences"
constraintderivative: "closed-form"
solver: 'fmincon'Оба решения сообщают одно и то же значение функции для отображения точности, и оба требуют одинакового количества итераций. Однако решение с градиентами требует только 32 оценки функций по сравнению с 64 для решения без градиентов.
Чтобы включить гессен, необходимо использовать prob2struct, даже если все функции поддерживаются для выражений оптимизации. В этом примере показано, как использовать гессен для fmincon
interior-point алгоритм. fminunc
trust-region алгоритм и fmincon
trust-region-reflective алгоритм использует другой синтаксис; см. Hessian алгоритмы fminunc trust-region или fmincon trust-region-reflective.
Как описано в Гессене для алгоритма fmincon interior-point, гессен - гессен лагранжиана.
| +∑λh,i∇2hi (x). | (1) |
Включите этот гессен, создав файл функции для вычисления гессена, и создав HessianFcn опция для fmincon чтобы использовать гессен. Чтобы создать гессен в этом случае, создайте по отдельности вторые производные целевого и нелинейного ограничений.
Вторая производная матрицы целевой функции для этого примера несколько сложна. Перечисление его целевой функции hessianfun(x) была создана инструментарием «Символьная математика» с использованием того же подхода, который описан в разделе «Расчет градиентов и гессенов с использованием символьных математических Toolbox™».
function hf = hessfun(in1) %HESSFUN % HF = HESSFUN(IN1) % This function was generated by the Symbolic Math Toolbox version 8.6. % 10-Aug-2020 10:50:44 x = in1(1,:); y = in1(2,:); t2 = x.^2; t4 = y.^2; t6 = x.*4.0e+2; t3 = t2.^2; t5 = t4.^2; t7 = -t4; t8 = -t6; t9 = t2+t4; t10 = t2.*t4.*2.0; t11 = besselj(0,t9); t12 = besselj(1,t9); t13 = t2+t7; t14 = 1.0./t9; t16 = t3+t5+t10-2.0; t15 = t14.^2; t17 = t11.*t14.*x.*y.*4.0; t19 = t11.*t13.*t14.*2.0; t18 = -t17; t20 = t12.*t15.*t16.*x.*y.*4.0; t21 = -t20; t22 = t8+t18+t21; hf = reshape([t2.*1.2e+3-t19-y.*4.0e+2-t12.*t15.*... (t2.*-3.0+t4+t2.*t5.*2.0+t3.*t4.*4.0+t2.^3.*2.0).*2.0+2.0,... t22,t22,... t19-t12.*t15.*(t2-t4.*3.0+t2.*t5.*4.0+... t3.*t4.*2.0+t4.^3.*2.0).*2.0+2.0e+2],[2,2]);
Напротив, гессен ограничения нелинейного неравенства прост; это в два раза больше единичной матрицы.
hessianc = 2*eye(2);
Создайте гессен лагранжиана как дескриптор функции.
H = @(x,lam)(hessianfun(x) + hessianc*lam.ineqnonlin);
Создайте параметры для использования этого гессенского языка.
problem.options.HessianFcn = H;
Решите проблему и просмотрите количество итераций и число оценок функций. Раствор примерно такой же, как и ранее.
[xhess,fvalhess,exitflaghess,outputhess] = fmincon(problem); disp(outputhess.iterations) disp(outputhess.funcCount)
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.
8
10На этот раз, fmincon принимает только 8 итераций по сравнению с 15 и только 10 оценок функций по сравнению с 32. Таким образом, обеспечение аналитического расчета Гессена может повысить эффективность процесса решения, но разработка функции для расчета Гессена может быть затруднена.
fcn2optimexpr | prob2struct | varindex