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

Рассмотрите задачу электростатики размещения 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, since we'll plot more 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 вектор, состоявший из действительных символьных переменных xiji между 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)$$

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

Нелинейные ограничения

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

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..., hesssc10.

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

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

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

$$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 may need to use currdir = pwd 
filename = [currdir,'demoenergy.m'];
matlabFunction(energy,gradenergy,'file',filename,'vars',{x});

Этот синтаксис вызывает matlabFunction возвратить энергию как первый выход и 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

Установите опции использовать алгоритм внутренней точки и использовать градиенты и Гессиан.

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

rotate3d
figure(hand)

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

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

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

   34.1365

disp(exitflag2)

     1

disp(output2.iterations)

    78

disp(output2.funcCount)

        2463

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

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        78         2463 

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

Символьные переменные в этом примере имеют предположение в символьной рабочей области механизма, что они действительны. Чтобы очистить это предположение от символьной рабочей области механизма, не достаточно удалить переменные. Очистите переменные предположения при помощи 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