Описание задачи

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

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

$$x^2+y^2+(z+1{)}^2\le 1$$.

Это тело выглядит как пирамида на сфере.

[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 = 
$$\left(\begin{array}{c}{x}_{11}\\ x_{12}\\ x_{13}\\ x_{21}\\ x_{22}\\ x_{23}\\ x_{31}\\ x_{32}\\ x_{33}\\ x_{41}\\ x_{42}\\ x_{43}\\ x_{51}\\ x_{52}\\ x_{53}\\ x_{61}\\ x_{62}\\ x_{63}\\ x_{71}\\ x_{72}\\ x_{73}\\ x_{81}\\ x_{82}\\ x_{83}\\ x_{91}\\ x_{92}\\ x_{93}\\ x_{101}\\ x_{102}\\ x_{103}\end{array}\right)$$

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

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

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

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)

$$\left(\begin{array}{c}{x}_{11}+x_{12}+x_{13}\\ x_{12}-x_{11}+x_{13}\\ x_{11}-x_{12}+x_{13}\\ x_{13}-x_{12}-x_{11}\\ x_{21}+x_{22}+x_{23}\\ x_{22}-x_{21}+x_{23}\\ x_{21}-x_{22}+x_{23}\\ x_{23}-x_{22}-x_{21}\\ x_{31}+x_{32}+x_{33}\\ x_{32}-x_{31}+x_{33}\\ x_{31}-x_{32}+x_{33}\\ x_{33}-x_{32}-x_{31}\\ x_{41}+x_{42}+x_{43}\\ x_{42}-x_{41}+x_{43}\\ x_{41}-x_{42}+x_{43}\\ x_{43}-x_{42}-x_{41}\\ x_{51}+x_{52}+x_{53}\\ x_{52}-x_{51}+x_{53}\\ x_{51}-x_{52}+x_{53}\\ x_{53}-x_{52}-x_{51}\\ x_{61}+x_{62}+x_{63}\\ x_{62}-x_{61}+x_{63}\\ x_{61}-x_{62}+x_{63}\\ x_{63}-x_{62}-x_{61}\\ x_{71}+x_{72}+x_{73}\\ x_{72}-x_{71}+x_{73}\\ x_{71}-x_{72}+x_{73}\\ x_{73}-x_{72}-x_{71}\\ x_{81}+x_{82}+x_{83}\\ x_{82}-x_{81}+x_{83}\\ x_{81}-x_{82}+x_{83}\\ x_{83}-x_{82}-x_{81}\\ x_{91}+x_{92}+x_{93}\\ x_{92}-x_{91}+x_{93}\\ x_{91}-x_{92}+x_{93}\\ x_{93}-x_{92}-x_{91}\\ x_{101}+x_{102}+x_{103}\\ x_{102}-x_{101}+x_{103}\\ x_{101}-x_{102}+x_{103}\\ x_{103}-x_{102}-x_{101}\end{array}\right)$$

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

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

$$x^2+y^2+(z+1{)}^2\le 1$$.

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

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первый столбец матрицы gradc. Эта форма верна, как описано в «Нелинейных ограничениях».

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

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

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

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

$$energy=\sum _{i<j}\frac{1}{\left|x_i-x_j\right|}$$.

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

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

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

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

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

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

Гессиан целевой функции, 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