Ограниченная электростатическая нелинейная оптимизация, основанная на проблеме

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

Для эквивалентного основанного на решателе примера с помощью Symbolic Math Toolbox™ смотрите, что Symbolic Math Toolbox Вычисляет Градиенты и Гессианы.

Проблемная геометрия

Этот пример вовлекает орган по проведению, заданный следующими неравенствами. Для каждого электрона с координатами $(x, y, z)$ ,

$\begin{array}{l} z \leq - | x | - | y | \\ x^{2} + y^{2} + (z + 1)^{2} \leq 1 . \end{array}$

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

[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$ и $y$ значение от –1 до 1, и $z$ значение от –2 до 0. Задайте переменные для проблемы.

N = 10;
x = optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y = optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z = optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob = optimproblem;

Задайте ограничения

Проблема имеет два типа ограничений. Первым, сферическим ограничением, является простое полиномиальное неравенство для каждого электрона отдельно. Задайте это сферическое ограничение.

elecprob.Constraints.spherec = (x.^2 + y.^2 + (z+1).^2) <= 1;

Предыдущая ограничительная команда создает вектор десяти ограничений. Просмотрите ограничительный вектор использование showconstr.

showconstr(elecprob.Constraints.spherec)

  ((x.^2 + y.^2) + (z + 1).^2) <= arg_RHS

  where:

    arg2 = 1;
    arg1 = arg2([1 1 1 1 1 1 1 1 1 1]);
    arg_RHS = arg1(:);

Второй тип ограничения в проблеме линеен. Можно выразить линейные ограничения по-разному. Например, можно использовать функцию abs, чтобы представлять ограничение абсолютного значения. Чтобы выразить ограничения этот путь, запишите функцию MATLAB и преобразуйте его в выражение с помощью fcn2optimexpr. Для другого подхода запишите ограничение абсолютного значения как четыре линейных неравенства. Каждая ограничительная команда возвращает вектор десяти ограничений.

elecprob.Constraints.plane1 = z <= -x-y;
elecprob.Constraints.plane2 = z <= -x+y;
elecprob.Constraints.plane3 = z <= x-y;
elecprob.Constraints.plane4 = z <= x+y;

Задайте целевую функцию

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

$e n e r g y = \sum_{i < j} \frac{1}{‖ e l e c t r o n (i) - e l e c t r o n (j) ‖} .$

Задайте целевую функцию как выражение оптимизации.

energy = optimexpr(1);
for ii = 1:(N-1)
    for jj = (ii+1):N
        tempe = (x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2;
        energy = energy + tempe^(-1/2);
    end
end
elecprob.Objective = energy;

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

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

rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);

Решите проблему путем вызова solve.

[sol,fval,exitflag,output] = solve(elecprob,init)

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.

sol = struct with fields:
    x: [10x1 double]
    y: [10x1 double]
    z: [10x1 double]

fval = 34.1365

exitflag = 
    OptimalSolution

output = struct with fields:
         iterations: 95
          funcCount: 2990
    constrviolation: 0
           stepsize: 1.8867e-06
          algorithm: 'interior-point'
      firstorderopt: 1.9841e-06
       cgiterations: 0
            message: '...'
             solver: 'fmincon'

Просмотрите решение

Постройте решение как точки на органе по проведению.

figure(hand)
plot3(sol.x,sol.y,sol.z,'r.','MarkerSize',25)
hold off

Электроны распределяются справедливо равномерно на ограничительном контуре. Электроны находятся в основном на ребрах и точке пирамиды.

Ссылка

[1] Долан, Элизабет Д., Хорхе Х. Море и Тодд С. Мансон. “Тестируя программного обеспечения оптимизации в сравнении с эталоном с COPS 3.0”. Технический отчет национальной лаборатории Аргонна ANL/MCS-TM-273, февраль 2004.

Документация