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

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