使用优化变量的约束静电非线性优化

打开实时脚本

考虑将 20 个电子放置在导体中的静电问题。电子会以最小化其总势能的方式排列，但受到位于体内的约束。所有电子至少都在物体的边界上。电子是无法区分的，因此问题没有唯一的最小值（在一个解中排列电子会给出另一个有效的解）。此示例受到 Dolan、Moré 和 Munson [1] 的启发。

当您使用优化变量（基于问题的方法）时，此示例的目标和非线性约束函数都是优化变量和表达式支持的运算。因此，solve 使用自动微分来计算梯度。请参阅Optimization Toolbox 中的自动微分。由于没有自动微分，此示例会因达到 MaxFunctionEvaluations 容差而提前停止。有关使用 Symbolic Math Toolbox ™ 的等效基于求解器的示例，请参阅使用 Symbolic Math Toolbox 计算梯度和 Hessian。

问题几何

此例涉及由下列不等式定义的导电体。对于每个坐标为 $(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, 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;

上述约束命令创建了一个 N 约束向量。使用 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 将其转换为表达式。请参阅将非线性函数转换为优化表达式。对于仅使用可微函数的首选方法，将绝对值约束写为四个线性不等式。每个约束命令返回一个 N 约束向量。

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)
    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;

运行优化

从以 [0,0,-1] 为中心、半径为 1/2 的球体上随机分布的电子开始优化。

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: 'Local minimum found that satisfies the constraints....'
            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] Dolan, Elizabeth D., Jorge J. Moré, and Todd S. Munson.“Benchmarking Optimization Software with COPS 3.0.”Argonne National Laboratory Technical Report ANL/MCS-TM-273, February 2004.