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

Этот пример включает проводящее тело, заданное следующими неравенствами. Для каждого электрона с координатами $$(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 года.