Вычислите градиенты и гессианы Используя Symbolic Math Toolbox™

В этом примере показано, как получить быстрее и больше надежных решений нелинейных задач оптимизации с помощью fmincon наряду с функциями Symbolic Math Toolbox. Если у вас есть лицензия Symbolic Math Toolbox, можно легко вычислить аналитические градиенты и Гессианы для цели и ограничительных функций с помощью этих функций Symbolic Math Toolbox:

  • jacobian (Symbolic Math Toolbox) генерирует градиент скалярной функции и генерирует матрицу частных производных вектор-функции. Так, например, можно получить матрицу Гессиана (вторые производные целевой функции) путем применения jacobian к градиенту. В этом примере показано, как использовать jacobian сгенерировать символьные градиенты и Гессианы ограничительных функций и цели.

  • matlabFunction (Symbolic Math Toolbox) генерирует или анонимную функцию или файл, который вычисляет значения символьного выражения. В этом примере показано, как использовать matlabFunction сгенерировать файлы, которые оценивают цель и ограничительные функции и их производные в произвольных точках.

Функции Symbolic Math Toolbox имеют различные синтаксисы и структуры по сравнению с функциями Optimization Toolbox™. В частности, символьные переменные являются действительными или комплексными скалярами, тогда как функции Optimization Toolbox передают аргументы вектора. Так, необходимо сделать несколько шагов, чтобы символически сгенерировать целевую функцию, ограничения и все их необходимые производные, в форме, подходящей для алгоритма внутренней точки fmincon.

Основанная на проблеме оптимизация может вычислить и использовать градиенты автоматически; смотрите Автоматическое Дифференцирование в Optimization Toolbox. Для подхода, основанного на проблеме к этой проблеме с помощью автоматического дифференцирования смотрите Ограниченную Электростатическую Нелинейную Оптимизацию, Основанную на проблеме.

Описание проблемы

Рассмотрите задачу электростатики размещения 10 электронов в органе по проведению. Электроны располагают себя способом, который минимизирует их общую потенциальную энергию согласно ограничению лжи в теле. Все электроны находятся на контуре тела как минимум. Электроны неразличимы, таким образом, никакой уникальный минимум не существует для этой проблемы (переставляющий электроны в одном решении, дает другое допустимое решение). Этот пример вдохновлен Доланом, Море и Мансоном [58].

Этим примером является орган по проведению, заданный неравенствами

z-|x|-|y|

x2+y2+(z+1)21.

Это тело похоже на пирамиду на сфере.

[X,Y] = meshgrid(-1:.01:1);
Z1 = -abs(X) - abs(Y);
Z2 = -1 - sqrt(1 - X.^2 - Y.^2);
Z2 = real(Z2);
W1 = Z1; W2 = Z2;
W1(Z1 < Z2) = nan; % Only plot points where Z1 > Z2
W2(Z1 < Z2) = nan; % Only plot points where Z1 > Z2
hand = figure; % Handle to the figure, for more plotting later
set(gcf,'Color','w') % White background
surf(X,Y,W1,'LineStyle','none');
hold on
surf(X,Y,W2,'LineStyle','none');
view(-44,18)

Небольшой разрыв существует между верхними и более низкими поверхностями фигуры. Этот разрыв является артефактом общей стандартной программы графического вывода, используемой, чтобы создать фигуру. Стандартная программа стирает любую прямоугольную закрашенную фигуру на одной поверхности, которая касается другой поверхности.

Создание переменных

Сгенерируйте символьный векторный x как 30 1 вектор, состоявший из действительных символьных переменных xij, с i между 1 и 10, и j между 1 и 3. Эти переменные представляют три координаты электронного i: xi1 соответствует x координата, xi2 соответствует y координата и xi3 соответствует z координата.

x = cell(3, 10);
for i = 1:10
    for j = 1:3
        x{j,i} = sprintf('x%d%d',i,j);
    end
end
x = x(:); % Now x is a 30-by-1 vector
x = sym(x, 'real');

Отобразите векторный x.

x
x = 

(x11x12x13x21x22x23x31x32x33x41x42x43x51x52x53x61x62x63x71x72x73x81x82x83x91x92x93x101x102x103)[x11; x12; x13; x21; x22; x23; x31; x32; x33; x41; x42; x43; x51; x52; x53; x61; x62; x63; x71; x72; x73; x81; x82; x83; x91; x92; x93; x101; x102; x103]

Включайте линейные ограничения

Запишите линейные ограничения

z-|x|-|y|

в виде набора четырех линейных неравенств для каждого электрона:

xi3 - xi1 - xi2 ≤ 0
xi3 - xi1 + xi2 ≤ 0
xi3 + xi1 - xi2 ≤ 0
xi3 + xi1 + xi2 ≤ 0

Так, эта проблема имеет в общей сложности 40 линейных неравенств.

Запишите неравенства структурированным способом.

B = [1,1,1;-1,1,1;1,-1,1;-1,-1,1];

A = zeros(40,30);
for i=1:10
    A(4*i-3:4*i,3*i-2:3*i) = B;
end

b = zeros(40,1);

Вы видите тот A*x ≤ b представляет неравенства.

disp(A*x)

(x11+x12+x13x12-x11+x13x11-x12+x13x13-x12-x11x21+x22+x23x22-x21+x23x21-x22+x23x23-x22-x21x31+x32+x33x32-x31+x33x31-x32+x33x33-x32-x31x41+x42+x43x42-x41+x43x41-x42+x43x43-x42-x41x51+x52+x53x52-x51+x53x51-x52+x53x53-x52-x51x61+x62+x63x62-x61+x63x61-x62+x63x63-x62-x61x71+x72+x73x72-x71+x73x71-x72+x73x73-x72-x71x81+x82+x83x82-x81+x83x81-x82+x83x83-x82-x81x91+x92+x93x92-x91+x93x91-x92+x93x93-x92-x91x101+x102+x103x102-x101+x103x101-x102+x103x103-x102-x101)[x11 + x12 + x13; x12 - x11 + x13; x11 - x12 + x13; x13 - x12 - x11; x21 + x22 + x23; x22 - x21 + x23; x21 - x22 + x23; x23 - x22 - x21; x31 + x32 + x33; x32 - x31 + x33; x31 - x32 + x33; x33 - x32 - x31; x41 + x42 + x43; x42 - x41 + x43; x41 - x42 + x43; x43 - x42 - x41; x51 + x52 + x53; x52 - x51 + x53; x51 - x52 + x53; x53 - x52 - x51; x61 + x62 + x63; x62 - x61 + x63; x61 - x62 + x63; x63 - x62 - x61; x71 + x72 + x73; x72 - x71 + x73; x71 - x72 + x73; x73 - x72 - x71; x81 + x82 + x83; x82 - x81 + x83; x81 - x82 + x83; x83 - x82 - x81; x91 + x92 + x93; x92 - x91 + x93; x91 - x92 + x93; x93 - x92 - x91; x101 + x102 + x103; x102 - x101 + x103; x101 - x102 + x103; x103 - x102 - x101]

Создайте нелинейные ограничения и их градиенты и гессианы

Нелинейные ограничения также структурированы.

x2+y2+(z+1)21.

Сгенерируйте ограничения и их градиенты и Гессианы.

c = sym(zeros(1,10));
i = 1:10;
c = (x(3*i-2).^2 + x(3*i-1).^2 + (x(3*i)+1).^2 - 1).';

gradc = jacobian(c,x).'; % .' performs transpose

hessc = cell(1, 10);
for i = 1:10
    hessc{i} = jacobian(gradc(:,i),x);
end

Ограничительный вектор c вектор-строка и градиент c(i) представлен в iстолбец th матричного gradc. Эта форма правильна, как описано в Нелинейных Ограничениях.

Матрицы Гессиана, hessc{1}..., hessc{10}, являются квадратными и симметричными. Этот пример хранит их в массиве ячеек, который лучше, чем хранение их в отдельных переменных, таких как hessc1..., hessc10.

Используйте .' синтаксис, чтобы транспонировать. ' синтаксис означает сопряженное транспонирование, которое имеет различные символьные производные.

Создайте целевую функцию и ее градиент и гессиан

Целевая функция, потенциальная энергия, является суммой инверсий расстояний между каждой электронной парой.

energy=i<j1|xi-xj|.

Расстояние является квадратным корнем из суммы квадратов различий в компонентах векторов.

Вычислите энергию (целевая функция) и ее градиент и Гессиан.

energy = sym(0);
for i = 1:3:25
    for j = i+3:3:28
        dist = x(i:i+2) - x(j:j+2);
        energy = energy + 1/sqrt(dist.'*dist);
    end
end

gradenergy = jacobian(energy,x).';

hessenergy = jacobian(gradenergy,x);

Создайте файл целевой функции

Целевая функция имеет два выходных параметров, energy и gradenergy. Поместите обе функции в один вектор при вызове matlabFunction сокращать количество подвыражений что matlabFunction генерирует, и возвратить градиент только когда функция вызова (fmincon в этом случае), запрашивает оба выходных параметров. Можно поместить получившиеся файлы в любую папку на пути MATLAB. В этом случае разместите их в свою текущую папку.

currdir = [pwd filesep]; % You might need to use currdir = pwd 
filename = [currdir,'demoenergy.m'];
matlabFunction(energy,gradenergy,'file',filename,'vars',{x});

matlabFunction возвращает energy как первый выход и gradenergy как второе. Функция также берет один входной вектор {x} вместо списка входных параметров x11..., x103.

Получившийся файл demoenergy.m содержит следующие линии или подобные единицы:

function [energy,gradenergy] = demoenergy(in1)
%DEMOENERGY
%   [ENERGY,GRADENERGY] = DEMOENERGY(IN1)
...
x101 = in1(28,:);
...
energy = 1./t140.^(1./2) + ...;
if nargout > 1
    ...
    gradenergy = [(t174.*(t185 - 2.*x11))./2 - ...];
end

Эта функция имеет правильную форму для целевой функции с градиентом; смотрите Пишущие Скалярные Целевые функции.

Создайте ограничительный файл функции

Сгенерируйте нелинейную ограничительную функцию и поместите ее в правильный формат.

filename = [currdir,'democonstr.m'];
matlabFunction(c,[],gradc,[],'file',filename,'vars',{x},...
    'outputs',{'c','ceq','gradc','gradceq'});

Получившийся файл democonstr.m содержит следующие линии или подобные единицы:

function [c,ceq,gradc,gradceq] = democonstr(in1)
%DEMOCONSTR
%    [C,CEQ,GRADC,GRADCEQ] = DEMOCONSTR(IN1)
...
x101 = in1(28,:);
...
c = [t417.^2 + ...];
if nargout > 1
    ceq = [];
end
if nargout > 2
    gradc = [2.*x11,...];
end
if nargout > 3
    gradceq = [];
end

Эта функция имеет правильную форму для ограничительной функции с градиентом; смотрите Нелинейные Ограничения.

Сгенерируйте файлы гессиана

Чтобы сгенерировать Гессиан функции Лагранжа для проблемы, сначала сгенерируйте файлы для энергетического Гессиана и ограничительных Гессианов.

Гессиан целевой функции, hessenergy, большое символьное выражение, содержащее более чем 150 000 символов, как показано путем оценки size(char(hessenergy)). Так, выполнение matlabFunction(hessenergy) занимает значительное количество времени.

Сгенерируйте файл hessenergy.m.

filename = [currdir,'hessenergy.m'];
matlabFunction(hessenergy,'file',filename,'vars',{x});

В отличие от этого Гессианы ограничительных функций малы и быстры, чтобы вычислить.

for i = 1:10
    ii = num2str(i);
    thename = ['hessc',ii];
    filename = [currdir,thename,'.m'];
    matlabFunction(hessc{i},'file',filename,'vars',{x});
end

После генерации всех файлов для цели и ограничений, помещает их вместе с соответствующими множителями Лагранжа в функции помощника hessfinal.m, чей код появляется в конце этого примера.

Запустите оптимизацию

Начните оптимизацию с электронов, распределенных случайным образом на сфере радиуса 1/2 сосредоточенный в [0,0, –1].

rng default % For reproducibility
Xinitial = randn(3,10); % Columns are normal 3-D vectors
for j=1:10
    Xinitial(:,j) = Xinitial(:,j)/norm(Xinitial(:,j))/2;
    % This normalizes to a 1/2-sphere
end
Xinitial(3,:) = Xinitial(3,:) - 1; % Center at [0,0,-1]
Xinitial = Xinitial(:); % Convert to a column vector

Установите опции использовать interior-point алгоритм, и использовать градиенты и Гессиан.

options = optimoptions(@fmincon,'Algorithm','interior-point','SpecifyObjectiveGradient',true,...
    'SpecifyConstraintGradient',true,'HessianFcn',@hessfinal,'Display','final');

Вызовите fmincon.

[xfinal,fval,exitflag,output] = fmincon(@demoenergy,Xinitial,...
    A,b,[],[],[],[],@democonstr,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.

<stopping criteria details>

Просмотрите значение целевой функции, выходной флаг, количество итераций и количество вычислений функции.

disp(fval)
   34.1365
disp(exitflag)
     1
disp(output.iterations)
    19
disp(output.funcCount)
    28

Даже при том, что исходные положения электронов случайны, конечные положения почти симметричны.

for i = 1:10
    plot3(xfinal(3*i-2),xfinal(3*i-1),xfinal(3*i),'r.','MarkerSize',25);
end

Чтобы исследовать целую 3-D геометрию, вращайте фигуру.

rotate3d

figure(hand)

Сравните с оптимизацией без градиентов и гессианов

Использование градиентов и Гессианов делает оптимизацию запущенной быстрее и более точно. Чтобы сравнить ту же оптимизацию, не используя градиента или информации о Гессиане, установите опции не использовать градиенты и Гессианы.

options = optimoptions(@fmincon,'Algorithm','interior-point',...
    'Display','final');
[xfinal2,fval2,exitflag2,output2] = fmincon(@demoenergy,Xinitial,...
    A,b,[],[],[],[],@democonstr,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.

<stopping criteria details>

Просмотрите значение целевой функции, выходной флаг, количество итераций и количество вычислений функции.

disp(fval2)
   34.1365
disp(exitflag2)
     1
disp(output2.iterations)
    80
disp(output2.funcCount)
        2525

Сравните количество итераций и количество вычислений функции с и без производной информации.

tbl = table([output.iterations;output2.iterations],[output.funcCount;output2.funcCount],...
    'RowNames',{'With Derivatives','Without Derivatives'},'VariableNames',{'Iterations','Fevals'})
tbl=2×2 table
                           Iterations    Fevals
                           __________    ______

    With Derivatives           19           28 
    Without Derivatives        80         2525 

Очистите символьные переменные предположения

Символьные переменные в этом примере имеют предположение, что они действительны в символьной рабочей области механизма. Удаление переменных не очищает это предположение от символьной рабочей области механизма. Очистите переменные предположения при помощи syms.

syms x

Проверьте, что предположения пусты.

assumptions(x)
 
ans =
 
Empty sym: 1-by-0
 

Функция помощника

Следующий код создает hessfinal функция помощника.

function H = hessfinal(X,lambda)
%
% Call the function hessenergy to start
H = hessenergy(X);

% Add the Lagrange multipliers * the constraint Hessians
H = H + hessc1(X) * lambda.ineqnonlin(1);
H = H + hessc2(X) * lambda.ineqnonlin(2);
H = H + hessc3(X) * lambda.ineqnonlin(3);
H = H + hessc4(X) * lambda.ineqnonlin(4);
H = H + hessc5(X) * lambda.ineqnonlin(5);
H = H + hessc6(X) * lambda.ineqnonlin(6);
H = H + hessc7(X) * lambda.ineqnonlin(7);
H = H + hessc8(X) * lambda.ineqnonlin(8);
H = H + hessc9(X) * lambda.ineqnonlin(9);
H = H + hessc10(X) * lambda.ineqnonlin(10);
end

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте