В этом примере показано, как использовать функции Symbolic Math Toolbox™ jacobian и matlabFunction предоставить аналитические производные решателям оптимизации. Optimization Toolbox™ решатели обычно более точны и эффективны, когда вы задаете градиенты и Гессианы целевых и ограничительных функций.
Основанная на проблеме оптимизация может вычислять и использовать градиенты автоматически; см. Раздел «Автоматическая дифференциация» в Optimization Toolbox. Для основанного на проблеме примера, использующего автоматическую дифференциацию, см. Ограниченная электростатическая нелинейная оптимизация, основанная на проблеме.
Существуют несколько факторов при использовании символьных вычислений с оптимизационными функциями:
Цель оптимизации и ограничительные функции должны быть определены в терминах вектора, скажем x. Однако символьные переменные являются скалярными или комплексно-значимыми, а не векторно-значимыми. Это требует, чтобы вы перемещались между векторами и скалярами.
Оптимизационные градиенты, а иногда и Гессианы, предполагается вычислить в теле целевой функции или функции ограничений. Это означает, что символьный градиент или Гессиан должен быть помещен в соответствующее место в файле целевой или ограничительной функции или указателе на функцию.
Вычисление градиентов и Гессианов символически может быть длительным. Поэтому вы должны выполнить это вычисление только один раз, и сгенерировать код, через matlabFunction, для вызова во время выполнения решателя.
Оценка символьных выражений с помощью subs функция занимает много времени. Гораздо эффективнее использовать matlabFunction.
matlabFunction генерирует код, который зависит от ориентации входа векторов. Начиная с fmincon вызывает целевую функцию с векторами-столбцами, вы должны быть осторожны, чтобы вызвать matlabFunction с векторами-столбцами символьных переменных.
Целевой функцией для минимизации является:
Эта функция положительная, с уникальным минимальным значением нуля, достигнутым в x1 = 4/3, x2 =(4/3)^3 - 4/3 = 1.0370...
Мы записываем независимые переменные как x1 и x2 потому что в этой форме они могут использоваться как символьные переменные. Как компоненты вектора x они будут написаны x(1) и x(2). Функция имеет твистый овраг, как изображено на графике ниже.
syms x1 x2 real x = [x1;x2]; % column vector of symbolic variables f = log(1 + 3*(x2 - (x1^3 - x1))^2 + (x1 - 4/3)^2)
f =
fsurf(f,[-2 2],'ShowContours','on') view(127,38)
Вычислите градиент и Гессиан f:
gradf = jacobian(f,x).' % column gradfgradf =
hessf = jacobian(gradf,x)
hessf =
The fminunc решатель рассчитывает пройти в векторе x, и, с SpecifyObjectiveGradient значение опции установлено в true и HessianFcn значение опции установлено в 'objective', ожидает список трех выходов: [f(x),gradf(x),hessf(x)].
matlabFunction генерирует точно этот список трех выходов из списка трех входов. Кроме того, использование vars опция, matlabFunction принимает векторные входы.
fh = matlabFunction(f,gradf,hessf,'vars',{x});Теперь решите задачу минимизации, начиная с точки [-1,2]:
options = optimoptions('fminunc', ... 'SpecifyObjectiveGradient', true, ... 'HessianFcn', 'objective', ... 'Algorithm','trust-region', ... 'Display','final'); [xfinal,fval,exitflag,output] = fminunc(fh,[-1;2],options)
Local minimum possible. fminunc stopped because the final change in function value relative to its initial value is less than the value of the function tolerance.
xfinal = 2×1
1.3333
1.0370
fval = 7.6623e-12
exitflag = 3
output = struct with fields:
iterations: 14
funcCount: 15
stepsize: 0.0027
cgiterations: 11
firstorderopt: 3.4391e-05
algorithm: 'trust-region'
message: '...'
constrviolation: []
Сравните это с количеством итераций, не используя никакой градиентной или гессианской информации. Для этого требуется 'quasi-newton' алгоритм.
options = optimoptions('fminunc','Display','final','Algorithm','quasi-newton'); fh2 = matlabFunction(f,'vars',{x}); % fh2 = objective with no gradient or Hessian [xfinal,fval,exitflag,output2] = fminunc(fh2,[-1;2],options)
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
xfinal = 2×1
1.3333
1.0371
fval = 2.1985e-11
exitflag = 1
output2 = struct with fields:
iterations: 18
funcCount: 81
stepsize: 2.1164e-04
lssteplength: 1
firstorderopt: 2.4587e-06
algorithm: 'quasi-newton'
message: '...'
Количество итераций ниже при использовании градиентов и гессианов, и вычислений функции значительно меньше:
sprintf(['There were %d iterations using gradient' ... ' and Hessian, but %d without them.'], ... output.iterations,output2.iterations)
ans = 'There were 14 iterations using gradient and Hessian, but 18 without them.'
sprintf(['There were %d function evaluations using gradient' ... ' and Hessian, but %d without them.'], ... output.funcCount,output2.funcCount)
ans = 'There were 15 function evaluations using gradient and Hessian, but 81 without them.'
Мы рассматриваем одну и ту же целевую функцию и начальную точку, но теперь имеем два нелинейных ограничения:
Ограничения держат оптимизацию подальше от глобальной минимальной точки [1.333,1.037]. Визуализируйте два ограничения:
[X,Y] = meshgrid(-2:.01:3); Z = (5*sinh(Y./5) >= X.^4); % Z=1 where the first constraint is satisfied, Z=0 otherwise Z = Z+ 2*(5*tanh(X./5) >= Y.^2 - 1); % Z=2 where the second is satisfied, Z=3 where both are surf(X,Y,Z,'LineStyle','none'); fig = gcf; fig.Color = 'w'; % white background view(0,90) hold on plot3(.4396, .0373, 4,'o','MarkerEdgeColor','r','MarkerSize',8); % best point xlabel('x') ylabel('y') hold off
Мы построили небольшой красный круг вокруг оптимальной точки.
Вот график объективной функции по выполнимой области, область, которая удовлетворяет оба ограничения, изображенные выше в темно-красном, наряду с маленьким красным кругом вокруг оптимальной точки:
W = log(1 + 3*(Y - (X.^3 - X)).^2 + (X - 4/3).^2); % W = the objective function W(Z < 3) = nan; % plot only where the constraints are satisfied surf(X,Y,W,'LineStyle','none'); view(68,20) hold on plot3(.4396, .0373, .8152,'o','MarkerEdgeColor','r', ... 'MarkerSize',8); % best point xlabel('x') ylabel('y') zlabel('z') hold off
Нелинейные ограничения должны быть записаны в форме c(x) <= 0. Мы вычисляем все символические ограничения и их производные и помещаем их в указатель на функцию, используя matlabFunction.
градиенты ограничений должны быть векторами-столбцами; они должны быть помещены в целевую функцию как матрица, причем каждый столбец матрицы представляет градиент одной ограничительной функции. Это транспонирование формы, сгенерированной jacobian, так что мы возьмем транспонирование ниже.
Мы помещаем нелинейные ограничения в указатель на функцию. fmincon ожидает, что нелинейные ограничения и градиенты будут выведены в порядке [c ceq gradc gradceq]. Поскольку нет нелинейных ограничений равенства, мы выводим [] для ceq и gradceq.
c1 = x1^4 - 5*sinh(x2/5); c2 = x2^2 - 5*tanh(x1/5) - 1; c = [c1 c2]; gradc = jacobian(c,x).'; % transpose to put in correct form constraint = matlabFunction(c,[],gradc,[],'vars',{x});
Алгоритм внутренней точки требует, чтобы его функция Гессия была записана как отдельная функция, вместо того, чтобы быть частью целевой функции. Это потому, что нелинейно ограниченная функция должна включать эти ограничения в свой Гессиан. Его Гессиан - Гессиан Лагрангиан; для получения дополнительной информации см. Руководство пользователя.
Функция Гессиана принимает два входных параметров: вектор положения x, и лямбда структуры множителя Лагранжа. Части лямбда-структуры, которые вы используете для нелинейных ограничений, lambda.ineqnonlin и lambda.eqnonlin. Для ограничения тока нет линейных равенств, поэтому мы используем два множителя lambda.ineqnonlin(1) и lambda.ineqnonlin(2).
Мы вычислили Гессиана целевой функции в первом примере. Теперь мы вычисляем Гессианы двух ограничительных функций и делаем указатель на функцию версии с matlabFunction.
hessc1 = jacobian(gradc(:,1),x); % constraint = first c column hessc2 = jacobian(gradc(:,2),x); hessfh = matlabFunction(hessf,'vars',{x}); hessc1h = matlabFunction(hessc1,'vars',{x}); hessc2h = matlabFunction(hessc2,'vars',{x});
Чтобы сделать окончательный Гессиан, мы сложим три Гессиана вместе, добавив соответствующие множители Лагранжа к ограничительным функциям.
myhess = @(x,lambda)(hessfh(x) + ... lambda.ineqnonlin(1)*hessc1h(x) + ... lambda.ineqnonlin(2)*hessc2h(x));
Установите опции, чтобы использовать алгоритм внутренней точки, градиент и Гессиан, имеют целевую функцию, возвращающую и цель, и градиент, и запустите решатель:
options = optimoptions('fmincon', ... 'Algorithm','interior-point', ... 'SpecifyObjectiveGradient',true, ... 'SpecifyConstraintGradient',true, ... 'HessianFcn',myhess, ... 'Display','final'); % fh2 = objective without Hessian fh2 = matlabFunction(f,gradf,'vars',{x}); [xfinal,fval,exitflag,output] = fmincon(fh2,[-1;2],... [],[],[],[],[],[],constraint,options)
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.
xfinal = 2×1
0.4396
0.0373
fval = 0.8152
exitflag = 1
output = struct with fields:
iterations: 10
funcCount: 13
constrviolation: 0
stepsize: 1.9160e-06
algorithm: 'interior-point'
firstorderopt: 1.9217e-08
cgiterations: 0
message: '...'
bestfeasible: [1x1 struct]
Снова, решатель делает много меньше итераций и вычислений функции с предоставленными градиентом и Гессианом, чем когда они не:
options = optimoptions('fmincon','Algorithm','interior-point',... 'Display','final'); % fh3 = objective without gradient or Hessian fh3 = matlabFunction(f,'vars',{x}); % constraint without gradient: constraint = matlabFunction(c,[],'vars',{x}); [xfinal,fval,exitflag,output2] = fmincon(fh3,[-1;2],... [],[],[],[],[],[],constraint,options)
Feasible point with lower objective function value found. 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.
xfinal = 2×1
0.4396
0.0373
fval = 0.8152
exitflag = 1
output2 = struct with fields:
iterations: 18
funcCount: 57
constrviolation: 0
stepsize: 9.6632e-07
algorithm: 'interior-point'
firstorderopt: 3.8435e-07
cgiterations: 0
message: '...'
bestfeasible: [1x1 struct]
sprintf(['There were %d iterations using gradient' ... ' and Hessian, but %d without them.'],... output.iterations,output2.iterations)
ans = 'There were 10 iterations using gradient and Hessian, but 18 without them.'
sprintf(['There were %d function evaluations using gradient' ... ' and Hessian, but %d without them.'], ... output.funcCount,output2.funcCount)
ans = 'There were 13 function evaluations using gradient and Hessian, but 57 without them.'
Символические переменные, используемые в этом примере, были приняты реальными. Чтобы удалить это предположение из рабочей области символьного двигателя, недостаточно удалить переменные. Вы должны очистить допущения переменных с помощью синтаксиса
assume([x1,x2],'clear')Все допущения удаляются, когда выход следующей команды пуст
assumptions([x1,x2])
ans = Empty sym: 1-by-0