В этом примере показано, как включать производную информацию в нелинейную основанную на проблеме оптимизацию, когда автоматические производные не применяются, или, когда это необходимо, включать Гессиан, который не обеспечивается с помощью автоматического дифференцирования. Включая градиенты или Гессиан в оптимизации может принести следующую пользу:
Больше устойчивых результатов. Конечные шаги дифференцирования иногда достигают точек, где цель или нелинейная ограничительная функция являются неопределенными, не конечными, или комплексными.
Увеличенная точность. Аналитические градиенты могут быть более точными, чем оценки конечной разности.
Более быстрая сходимость. Включая Гессиан может привести к более быстрой сходимости, означая меньше итераций.
Улучшенная производительность. Аналитические градиенты могут быть быстрее, чтобы вычислить, чем оценки конечной разности, специально для проблем с разреженной структурой. Для сложных выражений, однако, аналитические градиенты могут быть медленнее, чтобы вычислить.
Начиная в R2020b, solve функция может использовать автоматическое дифференцирование для поддерживаемых функций для того, чтобы предоставить градиенты к решателям. Эти автоматические производные применяются только к градиентам (первые производные), не Гессианы (вторые производные).
Автоматическое дифференцирование применяется, когда вы не используете fcn2optimexpr создать ограничительная функция или цель. Если необходимо использовать fcn2optimexpr, этот пример показывает, как включать производную информацию.
Способ использовать производные в основанной на проблеме оптимизации без автоматического дифференцирования состоит в том, чтобы преобразовать ваше использование задач prob2struct, и затем отредактируйте получившуюся цель и ограничительные функции. Этот пример показывает гибридный подход, где автоматическое дифференцирование предоставляет производные для части целевой функции.
С 2D контрольными переменными x и y, используйте целевую функцию
Включайте ограничение что сумма квадратов x и y не больше, чем 4.
fun2 не основан на поддерживаемых функциях для выражений оптимизации; смотрите Поддерживаемые Операции на Переменных и выражениях Оптимизации. Поэтому включать fun2 в задаче оптимизации необходимо преобразовать его в использование выражения оптимизации fcn2optimexprfun2 не поддерживается, потому что это использует неподдерживаемый sech функция. Однако sech(x) = 1/cosh(x), и cosh поддерживаемая функция. Использовать 1/cosh и автоматическое дифференцирование, смотрите Эффект Автоматического Дифференцирования в Основанной на проблеме Оптимизации.)
Чтобы использовать AD на поддерживаемых функциях, настройте проблему без неподдерживаемого функционального fun2, и включайте fun2 позже.
prob = optimproblem; x = optimvar('x',2); y = optimvar('y',2); fun1 = 100*(x(2) - x(1)^2)^2 + (1 - x(1))^2; fun2 = @(x,y)exp(-sum((x - y).^2))*exp(-exp(-y(1)))*sech(y(2)); prob.Objective = fun1; prob.Constraints.cons = sum(x.^2 + y.^2) <= 4;
Включать производные fun2, сначала преобразуйте проблему без fun2 к использованию структуры prob2struct.
problem = prob2struct(prob,... 'ObjectiveFunctionName','generatedObjectiveBefore');
Во время преобразования, prob2struct создает файлы функции, которые представляют объективные и нелинейные ограничительные функции. По умолчанию эти функции имеют имена 'generatedObjective.m' и 'generatedConstraints.m'. Файл целевой функции без fun2 'generatedObjectiveBefore.m'.
Сгенерированная цель и ограничительные функции включают градиенты.
Вычислите производные, сопоставленные с fun2. Если у вас есть лицензия Symbolic Math Toolbox™, можно использовать gradient (Symbolic Math Toolbox) или jacobian (Symbolic Math Toolbox) функция, чтобы помочь вычислить производные. Смотрите Вычисляют Градиенты и Гессианы Используя Symbolic Math Toolbox™.
Основанный на решателе подход имеет одну контрольную переменную. Каждая переменная оптимизации (x и y, в этом примере), фрагмент контрольной переменной. Можно найти отображение между переменными оптимизации и одной контрольной переменной в сгенерированном файле целевой функции 'generatedObjectiveBefore.m'. Начало файла содержит эти строки кода или подобные линии.
%% Variable indices. xidx = 1:2; yidx = 3:4; %% Map solver-based variables to problem-based. x = inputVariables(xidx); x = x(:); y = inputVariables(yidx); y = y(:);
В этом коде контрольная переменная имеет имя inputVariables.
В качестве альтернативы можно найти отображение при помощи varindex.
idx = varindex(prob); disp(idx.x)
1 2
disp(idx.y)
3 4
Полная целевая функция включает отрицание fun2:
fun2 = @(x,y)exp(-sum((x - y).^2))*exp(-exp(-y(1)))*sech(y(2));
Используя стандартное исчисление, вычислите gradfun2, градиент fun2.
где
Отредактируйте 'generatedObjectiveBefore.m' включать fun2.
%% Compute objective function. arg2 = x(1); arg4 = (x(2) - arg2.^2); arg5 = 100; arg6 = arg4.^2; arg8 = (1 - x(1)); obj = ((arg5 .* arg6) + arg8.^2); fun2 = exp(-sum((x - y).^2))*exp(-exp(-y(1)))*sech(y(2)); obj = obj - fun2;
Включайте расчетные градиенты в файл целевой функции путем редактирования 'generatedObjectiveBefore.m'. Если у вас есть версия программного обеспечения, которая не выполняет вычисление градиента, включает все эти линии. Если ваше программное обеспечение выполняет вычисление градиента, включайте только полужирные линии, которые вычисляют градиент fun2.
%% Compute objective gradient.
if nargout > 1
arg9 = 1;
arg10 = (-(arg9.*2.*(arg8(:))));
arg11 = zeros([2,size(arg10,2)]);
arg11(1,:) = arg10;
arg12 = ((arg9.*arg5(:)).*2.*(arg4(:)));
arg13 = zeros([2,size(arg12,2)]);
arg13(2,:) = arg12;
arg14 = ((-((arg9.*arg5(:)).*2.*(arg4(:)))).*2.*(arg2(:)));
arg15 = zeros([2,size(arg14,2)]);
arg15(1,:) = arg14;
arg16 = zeros([4, 1]);
arg16(xidx,:) = arg11 + arg13 + arg15;
grad = arg16(:);
% Calculate gradient of fun2
d = fun2;
sigma1 = -2*(x(1) - y(1))*d;
sigma2 = -2*(x(2) - y(2))*d;
gradfun2 = zeros(4, 1);
gradfun2(1) = sigma1;
gradfun2(2) = sigma2;
gradfun2(3) = -sigma1 + d*exp(-y(1));
gradfun2(4) = -sigma2 - d*tanh(y(2));
% Subtract the gradient of fun2
grad = grad - gradfun2;
endВспомните, что нелинейным ограничением является x(1)^2 + x(2)^2 + y(1)^2 + y(2)^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;2;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 =
0.8671
0.7505
1.0433
0.5140
fvalsolver =
-0.5500
exitflagsolver =
1
outputsolver =
struct with fields:
iterations: 46
funcCount: 77
constrviolation: 0
stepsize: 7.4091e-06
algorithm: 'interior-point'
firstorderopt: 7.5203e-07
cgiterations: 9
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, 7.520304e-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.↵↵'Сравните это решение с тем, полученным из solve без производной информации.
init.x = x0(1:2); init.y = x0(3:4); prob.Objective = prob.Objective - fcn2optimexpr(fun2,x,y); [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: [2×1 double]
y: [2×1 double]
fvalproblem =
-0.5500
exitflagproblem =
OptimalSolution
outputproblem =
struct with fields:
iterations: 47
funcCount: 269
constrviolation: 0
stepsize: 8.9896e-08
algorithm: 'interior-point'
firstorderopt: 1.3377e-08
cgiterations: 9
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, 1.337675e-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]
solver: 'fmincon'Оба решения сообщают о том же значении функции точности отображения, и оба требуют примерно того же количества итераций (46 информации о градиенте использования, 47 без). Однако решение с градиентами требует только 77 вычислений функции, по сравнению с 269 для решения без градиентов.
Чтобы включать Гессиан, необходимо использовать prob2struct, даже если все ваши функции поддерживаются для выражений оптимизации. В этом примере показано, как использовать Гессиан для fmincon
interior-point алгоритм. fminunc
trust-region алгоритм и fmincon
trust-region-reflective использование алгоритма различный синтаксис; смотрите Гессиан для fminunc доверительной области или области доверия fmincon отражающие алгоритмы.
Как описано в Гессиане для fmincon алгоритма внутренней точки Гессианом является Гессиан функции Лагранжа.
| (1) |
Включайте этот Гессиан путем создания файла функции, чтобы вычислить Гессиан и создания HessianFcn опция для fmincon использовать Гессиан. Чтобы создать Гессиан в этом случае, создайте вторые производные объективных и нелинейных ограничений отдельно.
Вторая производная матрица целевой функции для этого примера является сложной. Его целевая функция, перечисляющая hessianfun(x) был создан Symbolic Math Toolbox с помощью того же подхода как описано в, Вычисляют Градиенты и Гессианы Используя Symbolic Math Toolbox™.
function hessfun = hessianfun(in1) %HESSIANFUN % HESSFUN = HESSIANFUN(IN1) % This function was generated by Symbolic Math Toolbox version 8.5. % 15-Mar-2020 10:12:23 x1 = in1(1,:); x2 = in1(2,:); y1 = in1(3,:); y2 = in1(4,:); t2 = cosh(y2); t3 = sinh(y2); t4 = x1.*2.0; t5 = x2.*2.0; t6 = y1.*2.0; t7 = y2.*2.0; t8 = -y1; t10 = -y2; t15 = x1.*4.0e+2; t9 = -t6; t11 = -t7; t12 = 1.0./t2; t14 = exp(t8); t16 = t8+x1; t17 = t10+x2; t18 = -t15; t13 = t12.^2; t19 = t4+t9; t20 = t5+t11; t21 = t16.^2; t22 = t17.^2; t23 = -t14; t24 = exp(t23); t25 = t19.^2; t26 = t20.^2; t27 = -t21; t28 = -t22; t29 = t27+t28; t30 = exp(t29); t31 = t12.*t24.*t30.*2.0; t33 = t3.*t13.*t14.*t24.*t30; t34 = t3.*t13.*t19.*t24.*t30; t35 = t3.*t13.*t20.*t24.*t30; t36 = t12.*t24.*t25.*t30; t37 = t12.*t24.*t26.*t30; t42 = t12.*t14.*t19.*t24.*t30; t43 = t12.*t14.*t20.*t24.*t30; t44 = t12.*t20.*t23.*t24.*t30; t45 = t12.*t19.*t20.*t24.*t30; t32 = -t31; t38 = -t34; t39 = -t35; t40 = -t36; t41 = -t37; t46 = -t45; t49 = t43+t45; t47 = t18+t46; t48 = t38+t45; t50 = t32+t37+t39; t51 = t32+t36+t42; t52 = t33+t34+t44+t46; hessfun = reshape([t31+t40-x2.*4.0e+2+x1.^2.*1.2e+3+2.0,t47,... t51,t48,t47,t31+t41+2.0e+2,t49,t50,t51,... t49,t31+t40-t42.*2.0+t12.*t14.*t24.*t30-t12.*t24.*t30.*exp(t9),... t52,t48,t50,t52,t35.*2.0+t41+t12.*t24.*t30.*3.0-t3.^2.*t12.^3.*t24.*t30.*2.0],... [4,4]);
В отличие от этого Гессиан нелинейного ограничения неравенства прост; это - дважды единичная матрица.
hessianc = 2*eye(4);
Создайте Гессиан функции Лагранжа как указатель на функцию.
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.
31
59На этот раз, fmincon берет только 31 итерацию по сравнению с более чем 45 и только 59 вычислений функции по сравнению с более чем 75. Таким образом, обеспечение аналитического вычисления Гессиана может повысить эффективность процесса решения, но разработка функции, чтобы вычислить Гессиан может затруднить.
fcn2optimexpr | prob2struct | varindex