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

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

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

[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 later use
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 = 20;
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;

Предыдущая команда зависимости создает вектор из десяти зависимостей. Просмотр вектора зависимости с помощью show.

show(elecprob.Constraints.spherec)

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

  where:

    arg2 = 1;
    arg1 = arg2(ones(1,20));
    arg_RHS = arg1(:);

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

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

Определение целевой функции

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

Определите целевую функцию как выражение оптимизации. Для хорошей производительности запишите целевую функцию векторизированным образом. См. раздел Создание проблем эффективной оптимизации.

energy = optimexpr(1);
for ii = 1:(N-1)
    jj = (ii+1):N;
    tempe = (x(ii) - x(jj)).^2 + (y(ii) - y(jj)).^2 + (z(ii) - z(jj)).^2;
    energy = energy + sum(tempe.^(-1/2));
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)

Solving problem using fmincon.

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: [20x1 double]
    y: [20x1 double]
    z: [20x1 double]

fval = 163.0099

exitflag = 
    OptimalSolution

output = struct with fields:
              iterations: 94
               funcCount: 150
         constrviolation: 0
                stepsize: 2.8395e-05
               algorithm: 'interior-point'
           firstorderopt: 8.1308e-06
            cgiterations: 0
                 message: '...'
            bestfeasible: [1x1 struct]
     objectivederivative: "reverse-AD"
    constraintderivative: "closed-form"
                  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 года.