创建 `for` 循环进行静态分析

打开实时脚本

考虑在导体中放置 20 个电子的问题（请参阅使用优化变量的约束静电非线性优化）。电子以最小化其总势能的方式排列，但受到位于体内的约束。本例重点讲解了势能，它自然地表示为嵌套的 for 循环，以及如何转换 for 循环来执行静态分析。

无需修改即可创建静态分析问题

对于 N = 20 个电子，创建问题和优化变量。

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;
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;
for ii = 1:(N-1)
    for jj = (ii+1):N
        energy = energy + ((x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2)^(-1/2);
    end
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 求解问题。对解进行计时。

tic
[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: [20×1 double]
    y: [20×1 double]
    z: [20×1 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.↵↵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, 8.090884e-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.'
            bestfeasible: [1×1 struct]
     objectivederivative: "reverse-AD"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

toc

Elapsed time is 14.483005 seconds.

解时间超过 10 秒（您的结果可能会有所不同）。

修改问题以实现高效静态分析

该视频演示了这些步骤。

为了使优化表达式更有效地工作，提取 for 循环并将其转换为局部函数。首先，选择从 energy = optimexpr; 到最后的 end 的代码行。然后，右键单击选定的代码行并选择转换为本地函数。

编辑结果函数（出现在本示例末尾）如下：

将函数名称更改为 myenergy。
编辑第一行代码如下：

energy = function myenergy(N, x, y, z);

在第二行代码中，将初始值从 optimexpr 更改为 0。

为了利用静态分析，通过调用 fcn2optimexpr 将能量作为目标函数放入问题中。

energy = fcn2optimexpr(@myenergy,N,x,y,z);
elecprob.Objective = energy;

通过调用 solve 求解问题。对解进行计时。

tic
[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: [20×1 double]
    y: [20×1 double]
    z: [20×1 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.↵↵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, 8.090884e-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.'
            bestfeasible: [1×1 struct]
     objectivederivative: "reverse-AD"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

toc

Elapsed time is 0.686725 seconds.

解时间大大减少，降至一秒以下。

与无分析进行比较

您还可以通过将目标函数视为黑匣子来求解问题。黑箱目标不能利用自动微分，因此求解器必须采取有限差分步骤来估计目标梯度。然而，当底层函数能够快速评估时，黑盒就可以快速运行。比较设置 fcn2optimexpr 创建黑盒目标时的解时间。

energy = fcn2optimexpr(@myenergy,N,x,y,z,Analysis="off");
elecprob.Objective = energy;
tic
[sol,fval] = solve(elecprob,init)

Solving problem using fmincon.

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+03.

sol = struct with fields:
    x: [20×1 double]
    y: [20×1 double]
    z: [20×1 double]

fval = 163.3769

toc

Elapsed time is 0.392945 seconds.

求解器未完成优化，因为它超出了函数评估限制。此外，得到的 fval 比前两个解更高。增加函数评估限制并重试。

opts = optimoptions("fmincon",MaxFunctionEvaluations=1e4);
tic
[sol,fval,exitflag,output] = solve(elecprob,init,Options=opts)

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

fval = 163.0099

exitflag = 
    OptimalSolution

output = struct with fields:
              iterations: 96
               funcCount: 5978
         constrviolation: 0
                stepsize: 5.5563e-06
               algorithm: 'interior-point'
           firstorderopt: 6.9362e-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.902142e-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.'
            bestfeasible: [1×1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

toc

Elapsed time is 0.603105 seconds.

解时间与静态分析大致相同，目标函数值也是如此。output 结构表明求解器采用有限差分步骤来估计目标梯度。求解器需要进行数千次函数计算，而之前的情况只需要进行少于 200 次。如果您的函数评估很耗时，则这种方法可能效率低下。

总结

有时，您可以通过转换 for 循环进行静态分析来节省解时间。将 for 循环转换为函数相当简单：突出显示循环，右键单击以显示上下文菜单，选择创建局部函数或函数文件的选项，然后对新函数进行微小编辑。

无论如何，使用 fcn2optimexpr 将 for 循环表达式放入您的问题中。这种做法使您能够自动获得静态分析的任何未来改进，而无需重写代码。

辅助函数

以下代码创建 myenergy 辅助函数。

function energy = myenergy(N, x, y, z)
energy = 0;
for ii = 1:(N-1)
    for jj = (ii+1):N
        energy = energy + ((x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2)^(-1/2);
    end
end
end

另请参阅

fcn2optimexpr

创建 for 循环进行静态分析