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

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

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

$${}^{\mathrm{x2}}{+}^{\mathrm{}}y2+\mathrm{}{(}^{\mathrm{z}}\mathrm{+}$$1) 2≤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)$$

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

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

$${}^{\mathrm{x2}}{+}^{\mathrm{}}y2+\mathrm{}{(}^{\mathrm{z}}\mathrm{+}$$1) 2≤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.

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

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

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

$$\underset{}{}\frac{\mathrm{}}{{}_{}{}_{}energy=\sum 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