Этот пример показывает, как включить информацию производной в нелинейную основанную на проблеме оптимизацию, когда автоматические производные не применяются, или когда вы хотите включить Гессиана, который не предоставляется с помощью автоматической дифференциации. Включение градиентов или Гессиана в оптимизацию может дать следующие преимущества:
Более устойчивые результаты. Конечные шаги дифференцирования иногда достигают точек, где цель или нелинейная функция ограничения неопределенна, а не конечна или комплексна.
Повышенная точность. Аналитические градиенты могут быть более точными, чем конечные оценки различий.
Более быстрое сходимость. Включение Гессиана может привести к более быстрому сходимости, что означает меньше итераций.
Улучшенная производительность. Аналитические градиенты могут вычисляться быстрее, чем конечные оценки различий, особенно для задач с разреженной структурой. Однако для сложных выражений аналитические градиенты могут вычисляться медленнее.
Начиная с R2020b, solve функция может использовать автоматическую дифференциацию для поддерживаемых функций в порядок для передачи градиентов решателям. Эти автоматические производные применяются только к градиентам (первые производные), а не к Гессианам (вторые производные).
Автоматическая дифференциация применяется, когда вы не используете fcn2optimexpr для создания целевой функции или функции ограничений. Если вам нужно использовать fcn2optimexpr, этот пример показывает, как включить информацию о производной.
Способ использовать производные в основанной на проблеме оптимизации без автоматической дифференциации - преобразовать вашу задачу используя prob2struct, а затем отредактируйте получившиеся целевые и ограничительные функции. Этот пример показывает гибридный подход, где автоматическая дифференциация предоставляет производные для части целевой функции.
С управляющими переменными x и y, используйте целевую функцию
Включите ограничение, что сумма квадратов 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. Если у вас есть лицензия Symbolic Math Toolbox™, можно использовать gradient (Symbolic Math Toolbox) или jacobian (Symbolic Math Toolbox) функция, которая помогает вычислить производные. См. «Вычисление градиентов и гессианов с использованием символьных математических 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.
Редактирование '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 алгоритм использует другой синтаксис; см. Гессиан для fminunc trust-region или fmincon trust- области reflective алгоритмов.
Как описано в Гессиане для алгоритма внутренней точки fmincon, Гессиан является Гессианом Лагранжа.
| (1) |
Включите этого Гессиана, создав файл функции для вычисления Гессиана и создав HessianFcn опция для fmincon использовать Гессиана. Чтобы создать Гессиана в этом случае, создайте вторые производные от целевых и нелинейных ограничений отдельно.
Вторая производная матрица целевой функции для этого примера несколько сложна. Описание его целевой функции hessianfun(x) был создан Symbolic Math Toolbox с использованием того же подхода, как описано в Calculate Gradients and Hessians Using Symbolic Math 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