ScipyOptimizeDriver中的SLSQP只执行一次迭代，耗时很长，然后退出
SLSQP in ScipyOptimizeDriver only executes one iteration, takes a very long time, then exits
我正在尝试使用 SLSQP 来优化机翼的攻角，以将停滞点放置在所需位置。这纯粹是一个测试用例，用于检查我计算停滞位置分音的方法是否有效。
当 运行 使用 COBYLA 时，优化在 47 次迭代后收敛到正确的 alpha (6.04144912)。当 运行 使用 SLSQP 时，它完成一次迭代，然后挂起 非常 很长时间（10、20 分钟或更长时间，我没有准确计时），然后退出值不正确。输出是：
Driver debug print for iter coord: rank0:ScipyOptimize_SLSQP|0
--------------------------------------------------------------
Design Vars
{'alpha': array([0.5])}
Nonlinear constraints
None
Linear constraints
None
Objectives
{'obj_cmp.obj': array([0.00023868])}
Driver debug print for iter coord: rank0:ScipyOptimize_SLSQP|1
--------------------------------------------------------------
Design Vars
{'alpha': array([0.5])}
Nonlinear constraints
None
Linear constraints
None
Objectives
{'obj_cmp.obj': array([0.00023868])}
Optimization terminated successfully. (Exit mode 0)
Current function value: 0.0002386835700364719
Iterations: 1
Function evaluations: 1
Gradient evaluations: 1
Optimization Complete
-----------------------------------
Finished optimisation
为什么 SLSQP 会这样不正常？据我所知，当我查看 check_partials().
时，没有不正确的分析导数
代码很长，所以我把它放在Pastebin这里：
- 核心：https://pastebin.com/fKJpnWHp
- 无粘性：https://pastebin.com/7Cmac5GF
- 机翼坐标（NACA64-012）：https://pastebin.com/UZHXEsr6
你问了两个问题，但答案最终彼此无关：
- 为什么使用 SLSQP 时模型很慢，而使用 COBYLA 时模型很快
- 为什么 SLSQP 在一次迭代后停止？
1) 为什么SLSQP这么慢？
COBYLA 是一种无梯度方法。 SLSQP 使用梯度。因此，可以肯定的是，当 SLSQP 要求提供衍生品时（COBYLA 从未这样做过），速度就会放缓。
那是我首先去看的地方。计算导数分两步进行：a) 计算每个分量的分项和 b) 使用这些分项求解线性系统以计算总和。减速必须在这两个步骤之一中。
既然您可以 运行 check_partials 而不会有太多麻烦，步骤 (a) 不太可能是罪魁祸首。所以这意味着步骤 (b) 可能是我们需要加快速度的地方。
我在你的模型上 运行 摘要实用程序 (openmdao summary core.py) 并看到了这个：
============== Problem Summary ============
Groups: 9
Components: 36
Max tree depth: 4
Design variables: 1 Total size: 1
Nonlinear Constraints: 0 Total size: 0
equality: 0 0
inequality: 0 0
Linear Constraints: 0 Total size: 0
equality: 0 0
inequality: 0 0
Objectives: 1 Total size: 1
Input variables: 87 Total size: 1661820
Output variables: 44 Total size: 1169614
Total connections: 87 Total transfer data size: 1661820
然后我生成了你的模型的 N2 并看到了这个：
所以我们有一个 1169614 个元素长的输出向量，这意味着你的线性系统是一个大约 1e6x1e6 的矩阵。这相当大，您正在使用 DirectSolver 尝试 compute/store 分解它。这就是减速的根源。使用 DirectSolvers 非常适合较小的模型（根据经验，输出向量应少于 10000 个元素）。对于较大的，您需要更加小心并使用更高级的线性求解器。
在你的情况下，我们可以从 N2 中看到你的模型中的任何地方都没有耦合（N2 的下三角形中没有耦合）。像这样的纯前馈模型可以使用更简单和更快的 LinearRunOnce 求解器（如果您没有设置任何其他内容，这是默认设置）。所以我关闭了你模型中的所有 DirectSolvers，并且导数变得有效即时。让你的 N2 看起来像这样：
最佳线性求解器的选择与模型密切相关。要考虑的一个因素是计算成本，另一个是数值鲁棒性。这个问题在 Section 5.3 of the OpenMDAO paper 中有一些详细的介绍，我不会在这里涵盖所有内容。但这里非常简要地总结了关键考虑因素。
刚开始使用 OpenMDAO 时，使用 DirectSolver 是最简单而且通常是最快的选择。它很简单，因为它不需要考虑您的模型结构，而且速度很快，因为对于小型模型，OpenMDAO 可以 assemble 将雅可比矩阵转换为密集或稀疏矩阵，并提供直接分解。然而，对于较大的模型（或具有非常大的输出向量的模型），计算因式分解的成本高得令人望而却步。在这种情况下，您需要更有意识地分解求解器结构，并使用其他线性求解器（有时与直接求解器结合使用---see Section 5.3 of OpenMDAO paper, and this OpenMDAO doc）。
你说你想使用 DirectSolver 来利用稀疏雅可比存储。这是一个很好的直觉，但 OpenMDAO 的结构方式无论如何都不是问题。我们现在已经深入了解了，但既然你问了，我会给出一个简短的总结解释。从 OpenMDAO 3.7 开始，只有 DirectSolver 完全需要 assembled 雅可比行列式（事实上，是线性求解器本身决定了它所连接的任何系统）。所有其他 LinearSolvers 都使用 DictionaryJacobian（它存储每个子 jac，键控到 [of-var, wrt-var] 对）。每个 sub-jac 可以存储为密集或稀疏（取决于您如何声明该特定偏导数）。字典雅可比矩阵实际上是一种稀疏矩阵形式，尽管不是传统矩阵。这里的关键要点是，如果您使用 LinearRunOnce（或任何其他求解器），那么无论如何您都会获得内存高效的数据存储。只有 DirectSolver 转换为实际矩阵对象的更传统的程序集。
关于内存分配的问题。我从 openmdao 文档中借用了这张图片
2) 为什么SLSQP迭代一次就停止了？
基于梯度的优化对缩放非常敏感。我在你允许的设计 space 中绘制了你的 objective 函数并得到了这个：
所以我们可以看到最小值大约是 6 度，但是 objective 值很小（大约 1e-4）。
作为一般的经验法则，让 objective 达到大约 1 的数量级是个好主意（我们有一个 scaling report 功能可以帮助实现这一点）。我添加了一个与您的 objective:
数量级有关的参考
p.model.add_objective('obj', ref=1e-4)
然后我得到了一个不错的结果：
Optimization terminated successfully (Exit mode 0)
Current function value: [3.02197589e-11]
Iterations: 7
Function evaluations: 9
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Finished optimization
alpha = [6.04143334]
time: 2.1188600063323975 seconds
不幸的是，基于梯度的优化很难缩放。从将 objective/constraints 扩展到 order-1 开始是一个不错的经验法则，但对于更复杂的问题，您通常需要调整超出此范围的东西。
我正在尝试使用 SLSQP 来优化机翼的攻角，以将停滞点放置在所需位置。这纯粹是一个测试用例，用于检查我计算停滞位置分音的方法是否有效。
当 运行 使用 COBYLA 时，优化在 47 次迭代后收敛到正确的 alpha (6.04144912)。当 运行 使用 SLSQP 时，它完成一次迭代，然后挂起 非常 很长时间（10、20 分钟或更长时间，我没有准确计时），然后退出值不正确。输出是：
Driver debug print for iter coord: rank0:ScipyOptimize_SLSQP|0
--------------------------------------------------------------
Design Vars
{'alpha': array([0.5])}
Nonlinear constraints
None
Linear constraints
None
Objectives
{'obj_cmp.obj': array([0.00023868])}
Driver debug print for iter coord: rank0:ScipyOptimize_SLSQP|1
--------------------------------------------------------------
Design Vars
{'alpha': array([0.5])}
Nonlinear constraints
None
Linear constraints
None
Objectives
{'obj_cmp.obj': array([0.00023868])}
Optimization terminated successfully. (Exit mode 0)
Current function value: 0.0002386835700364719
Iterations: 1
Function evaluations: 1
Gradient evaluations: 1
Optimization Complete
-----------------------------------
Finished optimisation
为什么 SLSQP 会这样不正常？据我所知，当我查看 check_partials().
时，没有不正确的分析导数代码很长，所以我把它放在Pastebin这里：
- 核心：https://pastebin.com/fKJpnWHp
- 无粘性：https://pastebin.com/7Cmac5GF
- 机翼坐标（NACA64-012）：https://pastebin.com/UZHXEsr6
你问了两个问题，但答案最终彼此无关：
- 为什么使用 SLSQP 时模型很慢，而使用 COBYLA 时模型很快
- 为什么 SLSQP 在一次迭代后停止？
1) 为什么SLSQP这么慢？
COBYLA 是一种无梯度方法。 SLSQP 使用梯度。因此，可以肯定的是，当 SLSQP 要求提供衍生品时（COBYLA 从未这样做过），速度就会放缓。
那是我首先去看的地方。计算导数分两步进行：a) 计算每个分量的分项和 b) 使用这些分项求解线性系统以计算总和。减速必须在这两个步骤之一中。
既然您可以 运行 check_partials 而不会有太多麻烦，步骤 (a) 不太可能是罪魁祸首。所以这意味着步骤 (b) 可能是我们需要加快速度的地方。
我在你的模型上 运行 摘要实用程序 (openmdao summary core.py) 并看到了这个：
============== Problem Summary ============
Groups: 9
Components: 36
Max tree depth: 4
Design variables: 1 Total size: 1
Nonlinear Constraints: 0 Total size: 0
equality: 0 0
inequality: 0 0
Linear Constraints: 0 Total size: 0
equality: 0 0
inequality: 0 0
Objectives: 1 Total size: 1
Input variables: 87 Total size: 1661820
Output variables: 44 Total size: 1169614
Total connections: 87 Total transfer data size: 1661820
然后我生成了你的模型的 N2 并看到了这个：
所以我们有一个 1169614 个元素长的输出向量，这意味着你的线性系统是一个大约 1e6x1e6 的矩阵。这相当大，您正在使用 DirectSolver 尝试 compute/store 分解它。这就是减速的根源。使用 DirectSolvers 非常适合较小的模型（根据经验，输出向量应少于 10000 个元素）。对于较大的，您需要更加小心并使用更高级的线性求解器。
在你的情况下，我们可以从 N2 中看到你的模型中的任何地方都没有耦合（N2 的下三角形中没有耦合）。像这样的纯前馈模型可以使用更简单和更快的 LinearRunOnce 求解器（如果您没有设置任何其他内容，这是默认设置）。所以我关闭了你模型中的所有 DirectSolvers，并且导数变得有效即时。让你的 N2 看起来像这样：
最佳线性求解器的选择与模型密切相关。要考虑的一个因素是计算成本，另一个是数值鲁棒性。这个问题在 Section 5.3 of the OpenMDAO paper 中有一些详细的介绍，我不会在这里涵盖所有内容。但这里非常简要地总结了关键考虑因素。
刚开始使用 OpenMDAO 时，使用 DirectSolver 是最简单而且通常是最快的选择。它很简单，因为它不需要考虑您的模型结构，而且速度很快，因为对于小型模型，OpenMDAO 可以 assemble 将雅可比矩阵转换为密集或稀疏矩阵，并提供直接分解。然而，对于较大的模型（或具有非常大的输出向量的模型），计算因式分解的成本高得令人望而却步。在这种情况下，您需要更有意识地分解求解器结构，并使用其他线性求解器（有时与直接求解器结合使用---see Section 5.3 of OpenMDAO paper, and this OpenMDAO doc）。
你说你想使用 DirectSolver 来利用稀疏雅可比存储。这是一个很好的直觉，但 OpenMDAO 的结构方式无论如何都不是问题。我们现在已经深入了解了，但既然你问了，我会给出一个简短的总结解释。从 OpenMDAO 3.7 开始，只有 DirectSolver 完全需要 assembled 雅可比行列式（事实上，是线性求解器本身决定了它所连接的任何系统）。所有其他 LinearSolvers 都使用 DictionaryJacobian（它存储每个子 jac，键控到 [of-var, wrt-var] 对）。每个 sub-jac 可以存储为密集或稀疏（取决于您如何声明该特定偏导数）。字典雅可比矩阵实际上是一种稀疏矩阵形式，尽管不是传统矩阵。这里的关键要点是，如果您使用 LinearRunOnce（或任何其他求解器），那么无论如何您都会获得内存高效的数据存储。只有 DirectSolver 转换为实际矩阵对象的更传统的程序集。
关于内存分配的问题。我从 openmdao 文档中借用了这张图片
2) 为什么SLSQP迭代一次就停止了？
基于梯度的优化对缩放非常敏感。我在你允许的设计 space 中绘制了你的 objective 函数并得到了这个：
p.model.add_objective('obj', ref=1e-4)
然后我得到了一个不错的结果：
Optimization terminated successfully (Exit mode 0)
Current function value: [3.02197589e-11]
Iterations: 7
Function evaluations: 9
Gradient evaluations: 7
Optimization Complete
-----------------------------------
Finished optimization
alpha = [6.04143334]
time: 2.1188600063323975 seconds
不幸的是，基于梯度的优化很难缩放。从将 objective/constraints 扩展到 order-1 开始是一个不错的经验法则，但对于更复杂的问题，您通常需要调整超出此范围的东西。