Рассмотрите задачу электростатики размещения 10 электронов в органе по проведению. Электроны расположат себя способом, который минимизирует их общую потенциальную энергию согласно ограничению лжи в теле. Все электроны находятся на контуре тела как минимум. Электроны неразличимы, таким образом, проблема не имеет никакого уникального минимума (перестановка электронов в одном решении дает другое допустимое решение). Этот пример был вдохновлен Доланом, Море и Мансоном [1].
Для эквивалентного основанного на решателе примера с помощью Symbolic Math Toolbox™ смотрите, что Symbolic Math Toolbox™ Вычисляет Градиенты и Гессианы.
Этот пример вовлекает орган по проведению, заданный следующими неравенствами. Для каждого электрона с координатами ,
Эти ограничения формируют тело, которое похоже на пирамиду на сфере. Чтобы просмотреть тело, введите следующий код.
[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)
Небольшой разрыв существует между верхними и более низкими поверхностями фигуры. Этот разрыв является артефактом общей стандартной программы графического вывода, используемой, чтобы создать фигуру. Стандартная программа стирает любую прямоугольную закрашенную фигуру на одной поверхности, которая касается другой поверхности.
Проблема имеет десять электронов. Ограничения дают границы на каждом и значение от –1 до 1, и значение от –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;
Предыдущая ограничительная команда создает вектор десяти ограничений. Просмотрите ограничительный вектор использование show
.
show(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;
Целевая функция является потенциальной энергией системы, которая является суммой по каждой электронной паре инверсии их расстояний:
Задайте целевую функцию как выражение оптимизации.
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)
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. <stopping criteria details>
sol = struct with fields:
x: [10×1 double]
y: [10×1 double]
z: [10×1 double]
fval = 34.1365
exitflag = OptimalSolution
output = struct with fields:
iterations: 94
funcCount: 2959
constrviolation: 0
stepsize: 7.0129e-06
algorithm: 'interior-point'
firstorderopt: 2.0944e-06
cgiterations: 0
message: '↵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>↵↵Optimization completed: The relative first-order optimality measure, 6.083643e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
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.