cashocs：基于伴随法的形状优化和最优控制软件

软件X 13（2021）100646原始软件出版物cashocs：一个基于伴随的计算形状优化和最优控制软件塞巴斯蒂安·布劳特Fraunhofer ITWM，Kaiserlestern，GermanyTU Kaiserlestern，Kaiserlestern，Germanyar t i cl e i nf o文章历史记录：收到2020年2020年12月9日收到修订版2020年12月10日接受保留字：PDE约束优化伴随法形状优化最优控制a b st ra ct偏微分方程约束优化问题的求解在科学和工业的许多领域中起着重要的作用。在这项工作中，我们提出了cashocs，一个新的软件包编写的Python，自动解决这些问题的最优控制和形状优化的背景下。该软件cashocs实现了连续伴随方法的离散化，该方法以自动化的方式导出必要的伴随系统和（形状）导数。由于cashocs基于有限元软件FEniCS，因此继承了其简单、高级的用户界面。这使得使用我们的软件定义和解决PDE约束优化问题变得简单。在本文中，我们讨论了卡肖克的设计和功能，并展示了其简单的可用性和适用性。©2020作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0.3用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-20-00059Code Ocean compute capsule NAGNU GPL v3.0（或更高版本）使用git的代码版本控制系统使用Python、FeNiCS、NumPy、PETSc、meshio、Gmsh的软件代码语言、工具和服务编译要求、操作环境依赖性FeNiCS、meshio、Gmsh、matplotlib如果可用，链接到开发人员文档/手册https://cashocs.readthedocs.io/问题支持电子邮件sebastian. itwm.fraunhofer.de1. 动机和意义由偏微分方程（PDE）及其数值解约束的形状优化和最优控制问题在科学和工业的许多领域都很重要：例如，它们用于化学反应器[1]、玻璃冷却过程[2]和半导体[3]的优化以及冷却系统[4]、飞机[5]和电机[6]的优化设计。为了解决这些问题，通常采用所谓的伴随方法，该方法便于计算问题的（形状）梯度，可以用于数值求解然而，对于复杂的，耦合的，或高度非线性的问题，如从工业应用中产生的问题，即使是推导必要的方程的伴随方法是一个非常复杂和容易出错的任务。电子邮件地址：sebastian. itwm.fraunhofer.de。https://doi.org/10.1016/j.softx.2020.100646因此，不再手动执行伴随方法是不可行的（参见，例如，[7]）。由于这些原因，最近已经有很多努力来自动化求解PDE约束优化问题的任务，产生了dolfin-adjoint [8]和Fireshape [9]等软件，有限 元 软 件 NGSolve [10] 的 形 状 优 化 功 能 ， 以 及 我 们 的 软 件cashocs。CASHOC与其他软件包的区别在于其新颖的方法，即仅使用自动微分来推导伴随系统和（形状）导数，同时在所有其余方面实现和自动化连续伴随方法的离散化。这意味着优化算法与所有所需的操作一起被实现为底层无限维操作的离散化。上述操作包括，例如，根据计算的（形状）导数确定（形状）梯度，2352-7110/©2020作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx塞巴斯蒂安·布劳特软件X 13（2021）1006462∈∈×→J=∈J×→⊆∈∈y，=Fig. 1. 卡肖克的建筑状态方程和伴随方程，标量积的计算，以及投影算子的使用。因此，计算的（形状）导数仅用作特别是，cashoc实现了与底层优化问题密切相关的连续、无限维优化算法的离散化，而其他软件包使用外部优化库[8，9]，或者要求用户自己实现这些算法[10]。我们的方法导致了独特的功能，如离散化和解决状态和伴随系统的不同的可能性，以及选择标量积的（形状）梯度的计算，也产生了网格独立的行为，如第3节所示。此外，cashocs是这些软件包中唯一一个实现了形状优化问题的重新网格化功能。1.1. 数学背景让我们从陈述我们的软件可以解决的优化问题最优控制问题具有以下形式：min（y，u）S.T.e（y，u）0和uUad，（1）y， u其中uU和yY是控制和状态变量，U和Y是适当的Banach空间，并且容许控制的集合U和U用于对控制变量的附加约束进行建模。此外，委员会认为，：是UR是成本泛函，e：Y UZ是PDE约束，我们在以下弱意义找到y∈Y，使得对所有p ∈ Z，有<$e（y，u），p<$Z<$，Z = 0。给你，Z 表示Z的拓扑对偶空间，且x ∈N，x∈N特别地，这意味着PDE约束是在域上给出的，并且该域的形状受到优化。问题（1）和（2）是卡肖克可以解决的问题的原型，我们请读者参阅第3节的说明性例子。如前所述，这类问题通常用伴随方法解决，其推导超出了本文的范围。因此，我们建议读者参考[11，12]和[132. 软件描述2.1. 软件构架为了用现金流解决优化问题，用户必须做以下事情。首先，他们必须在Python脚本中实现问题，包括计算网格、状态系统和成本泛函的定义。做因此，它们可以使用与定义问题相同的语法在FeNiCS [16，17]中，只有很小的修改，导致一个简单的，高层次的用户界面，支持许多重要类型的优化问题。其次，用户必须定义一个配置文件，该文件指定状态系统和优化算法的解决方案的参数，该文件被加载到用户脚本中。然后，可以设置优化问题使用 cashocs.OptimalControlProblem或cashocs。形状优化问题，分别用相应类的solve方法解决它。在内部，我们的软件利用符号自动微分法，统一形式语言[16，18]的能力来计算所需的（形状）导数和变分公式，表示Z和xZ的对偶配对。形状优化问题具有以下形式：Z，Z伴随系统此外，cashocs使用FEniCS生成和编译C++代码，用于问题的有限元组装，PETSc [19]用于解决出现的线性最小J（y，λ）标准差e（y，n）= 0且n ∈ A，（2）其中，y再次是状态变量，并且容许域A的集合用于结合附加的几何约束。我们在下面解释PDE约束e（y，y）0： 意识淡薄求y∈Y（n），使得e（y，对所有p∈Z（n）.系统，这使得问题的解决非常有效。卡肖克体系结构的示意图见图11。1.一、2.2. 软件功能我们的软件可以处理稳态和瞬态条件下PDE约束的线性和非线性系统，只要它们可以在FEniCS中实现为变分公式（序列）。此外，cashoc使用投影技术处理附加的控制约束，并且可以用于解决状态约束问题，例如，通过塞巴斯蒂安·布劳特软件X 13（2021）100646322{001×2Ω表1图二、 问题（3）的数值解，用 戴 元 NC G 方 法 计算。3. 说明性实例需要迭代以达到问题（3）的停止准则。nGDNCGL-BFGS牛顿163210613233106164331061128331061正则化（参见，例如，[11]）。我们提出了两个受泊松方程约束的模型问题https://cashocs.readthedocs.io/en/latest/tutorial_index。为 了 证 明 我 们 的 软 件 请 注 意 ， 可 以 在https://cashocs.readthedocs.io/en/latest/tutorial_index.html 的 教程中找到使用cashoc的各种其他示例。3.1. 最优控制作为模型最优控制问题，我们考虑[11]html来详细描述我们的软件最小J（y，u）=1（y−yd）2dx+αu2dx（三）以下算法可用于cashocs中的形状优化和最优控制问题• 梯度下降法，• 非线性共轭梯度法（NCG），• 有限记忆BFGS（L-BFGS）方法。注意，特别是对于形状优化，这些算法对应于现有技术，其中L-BFGS方法在[14]中引入，NCG方法在[15]中引入。此外，本发明还服从−y=u在美国，y= 0on r= 0这个最优控制问题有一个跟踪型成本功能与Tikhonov正则化的控制变量。PDE约束由具有齐次Dirichlet边界条件的Poisson方程给出，并且控制变量u作为右手侧进入PDE。该PDE约束的弱公式由下式给出：找到y∈H1（n），使得下列优化算法仅适用于最优控制问题• 截断牛顿法，• 一种注意，对于最优控制问题，所有方法也可以使用投影对于所有p∈ H1（n），dx = 0。对于这个例子，让我们使用α=（0，1）2，α=1e−4，y d（x）= x2（1 − x1）x2（1 − x2）。（四）技术.用户可以使用配置文件调整这些算法的行为，例如，可以修改相对和绝对停止容限、最大迭代次数以及其它算法特定的参数。CASHOCS的其他特征包括，对状态和伴随系统使用不同离散化的可能性，用于求解耦合系统的Picard迭代的实现，指定用于计算（形状）梯度的标量积的可能性，以及形状优化问题的重新网格化，其利用网格生成软件Gmsh [20]。对于域的离散化，我们使用一个统一的三角形网格划分成n n个正方形，减半创建三角形。为了解决这个cashocs问题，我们可以使用清单1所示的代码，我们将在下面简要讨论。请注意，由于我们的软件是基于FeNiCS的，我们建议读者参考[16，第1章]，其中使用几个描述性示例解释了FeNiCS的语法在清单1中，我们首先在第1行和第2行导入FEniCS和cashoc。接下来，我们用UnitSquareMesh函数定义网格，并设置体积测量以进行积分，见第5-7行。随后，我们在第10行中定义线性拉格朗日元的函数空间，并定义函数y，p和u。这些用于定义ΩΩ塞巴斯蒂安·布劳特软件X 13（2021）1006464=J={212123456789101112131415161718192021222324252627282930清单1：使用cashocs解决问题（3）1234567清单2：问题（3）的最小配置文件js.ini。第17行中的PDE约束的弱形式，其中函数p扮演测试函数的角色。在第19和20行中，定义了泊松问题的Dirichlet边界条件请注意，到目前为止，我们只使用来自FeNiCS的命令，并进行了以下微小的更改。而不是将y定义为TrialFunction，将p定义为TestFunction，现在两者都是并且对于用于最优控制问题的配置文件的详细描述，我们参考www.example.com/en/latest/demos/optimal_control/doc_soft.html。https://cashocs.readthedocs.io/注意，用于计算成本函数梯度的标量积可以由用户确定，如教程中所述。默认配置使用L2（k）标量积，这适用于我们的模型问题.用戴元非线性CG法计算的最优控制和状态图如图所示。二、此外，表1显示了优化算法需要使用n16、32、64和128细分来解决一系列更精细网格的问题所需的迭代量。我们观察到，所有的算法都表现出与网格无关的行为，因为无论离散化如何，它们基本上都需要相同数量的迭代来收敛。3.2. 形状优化作为形状优化的模型问题，我们考虑[15，21]函数对象。此外，我们不使用（线性）PDE约束的左侧和右侧定义PDE约束，而是定义min（y，y）y，ydxΩ（五）它就像我们在FEniCS中的非线性变分问题一样，类似于（1）和（4）中的形式。在第23和24行中，我们定义了期望状态，其在第25行中用于定义成本函数。同样，我们只使用FEniCS命令进行这些操作。要调用cashocs来解决这个问题，我们所要做的就是在第28行将配置文件加载到脚本中，在第29行初始化OptimalControlProblem，然后调用它的solve方法。总之，我们只需要添加三行额外的代码来解决这个问题。注意，清单2中显示了代码的最小配置文件，受−y=fin，y= 0on r= 0对于这个问题，PDE约束再次由具有齐次Dirichlet边界条件的Poisson问题给出，因此其弱形式由（4）给出，其中u由f代替。我们继续类似于[15，21]，并使用R2中的单位圆作为域R0的初始猜测，并使用右手边的单位圆作为域R0的初始猜测。如果我们使用f（x）= 2. 5（x1+ 0. 4− x2）2 + x2+ x2− 1。fromfenicsimportImportcashocs#definemeshanddvolumemeasuren = 64mesh=UnitSquareMesh（n，n）dx=Measure（#functionspaceforlinearLagragelementsV=Function Space（mesh，#定义状态、伴随和控制变量y=Function（V）p=Function（V）u=Function（V）#definehewekfhePDee=inner（grad（y），grad（p））<$dx−u<$p <$dx#定义的边界条件bdry=CompiledSubDomain（#定义期望的状态和成本函数x=SpatialCordinate（mesh）J=0。5pow（y− y_d，2）<$dx+1e−4/2<$pow（u，2）<$dx#solvetheomizy_d=pow（x[0]，2）<$（1−x[0]）<$ pow（x[1]，2）<$（1−x[1]）ocp=cashocs。optimalControlProblem（e，bcs，J，y，u，p，cfg）ocp.solve（）[操作时间限制]algorithm=ncgrtol= 1e−3maximum_iterati ns=50#在一张纸上写着#...塞巴斯蒂安·布劳特软件X 13（2021）1006465==-V=V+V1234567891011121314151617181920212223242526272829303132清单3：使用cashocs解决问题（5）1234567891011清单4：问题（5）的最小配置文件js.ini。我们通过将圆划分为n个较小的条带，然后均匀地网格化，使用均匀的三角形网格来离散N0这个问题可以通过使用清单3中提供的代码使用cashocs来解决，我们将在下面简要讨论如前所述，我们参考[16，第1章]详细介绍FEniCS的语法，我们也将其用于cashocs中的问题定义。该代码与清单1中的代码非常相似， 有一个泊松方程作为PDE约束。我们通过导入FEniCS和cashoc来启动脚本。然后，我们定义网格和体积测量，现在使用函数UnitDiscMesh，在5-7行。对于泊松方程的离散化，我们再次使用线性拉格朗日元，其对应的函数空间在第10行中定义，并且函数y和p被定义为在 12 和 13 行 。 此 后 ， 我 们 在 第 16 行 和 第 17 行 中 使 用SpatialCoordinate至于最优控制问题，与传统FEniCS语法的唯一主要区别是y和p是Function对象，并且PDE约束是在（2）和 （ 4 ） 的 意 义 上 编 写 的 。 随 后 ， 我 们 为 边 界 建 立 了FEniCSMeshFunction此外，这用于定义通过配置文件固定哪些边界清单4的第7行和第8行）。最后，我们在第27行中定义成本泛函。为了用cashoc解决这个问题，我们类似于清单1，首先加载配置文件，然后设置ShapeOptimizationProblem，最后在第30-32行调用它的solve方法注意，最小的清单4显示了这个问题的配置文件。 一个德-用 于 形 状 优 化 的 配 置 文 件 的 详 细 讨 论 可 以 在www.example.comshape_optimization/doc_optimization.html上找到。https://cashocs.readthedocs.io/en/latest/demos/注意，用于计算形状梯度的标量积基于线性弹性方程（参见，例如，[15、21、22]）。相应的双线性形式a由下式给出：a （ V， W）=ε2μ ε （ V）： ε （ W）+λdiv（V ）div（ W）+δV·Wdx，Ω其中ε（）1/2（DD）是雅可比矩阵的对称部分。弹性参数的默认值为µ1、λ 0和δ0，可以通过配置文件更改。 图1给出了用戴元NCG方法计算的最优域上的最优状态图。3 .第三章。fromfenicsimportImportcashocs#definemeshanddvolumemeasuren = 64mesh=UnitDiscMesh。create（MPI. comm_world，n，1，2）dx=Measure（#functionspaceofflinearLagrangeelementsV=Function Space（mesh，#状态和伴随变量y=Function（V）p=Function（V）#右手边x=Spat ialCordinate（mesh）F=2. 5个pow（x[0]+0. 4−pow（x[1]，2），2）+pow（x[0]，2）+pow（x[1]，2）−1#定义PDE约束e=inner（grad（y），grad（p））<$dx− f<$p <$dx#定义的边界条件bdry=CompiledSubDomain（' on_bound a r y '）m f_b d r y = M e s h F un c ti o n（' s i ze _ t '，m es h，d i m = 1）bd r y.mark（mf_bdry，1）bcs=DiricletBC（V，Constant（0），mf_bdry，1）#成本泛函J = y<$dx#solveepro blemwithcasocscfg=cashocs。create_config（' c on f i g. inisop=cashocs. ShapeOpt imin Problem（e，bcs，J，y，p，mf_bdry，cfg）sop.solve（）[操作时间限制]algorithm=ncgrtol= 5e−3maximum_iterati ns=50[ShapeGradient]shape_bdry_def = [1]shape_bdry_fix =[]#在一张纸上写着#...塞巴斯蒂安·布劳特软件X 13（2021）1006466=图3.第三章。 问题（5）最优域上最优状态y的数值解用Dai-Yuan NCG方法计算的表2需要迭代以达到问题（5）的停止准则。nGDNCGL-BFGS164620123247191164471911128471911此外，我们还显示了在表2中的单位圆的n16，32，64和128条连续精细离散化的算法所需的迭代次数。和前面一样，我们看到迭代次数基本上保持不变，与离散化无关，这表明我们也有形状优化问题的网格4. 影响我们的软件使用户能够以自动化的方式处理复杂、耦合和高度非线性的PDE约束优化问题。用户只需要使用与在FEniCS中定义这些对象基本相同的语法来定义PDE约束和成本函数。借助高级用户界面，只需添加三行额外的代码即可解决相应的优化问题。如第3节所示，我们实现连续伴随方法的离散化的方法导致优化算法的网格独立行为，使我们的软件对科学和工业具有吸引力。实际上，在文献[1]中，cashocs已经被用于处理化学微反应器中的参数识别和最优控制的高度非线性优化问题。在[15]中，它也被用于形状优化的NCG方法的数值基准。此外，在Fraunhofer ITWM中使用cashocs来解决工业应用中的PDE约束优化问题。由于我们的软件的通用性，它可以处理许多重要的类的成本泛函和PDE约束，它可以应用于许多相关的问题，在科学和工业，自动化的解决方案，在一个有效的和用户友好的方式。5. 结论本文介绍了一个求解偏微分方程约束形状优化和最优控制问题的软件--卡肖克该软件自动导出所需的伴随系统和（形状）导数，并实现连续伴随方法的离散化我们的软件继承了FEniCS高级用户界面，允许直接定义和求解PDE约束优化问题。此外，用户仍然保留对优化的许多重要参数的控制，从PDE的解决方案到优化算法，这允许他们对问题的数值解进行精确调整。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作确认作者感谢德国弗劳恩霍夫工业数学研究所ITWM的财政支持。引用[1]Blauth S，Leithäuser C，Pinnau R.微通道反应器中Sabatier过程的优化控制2020，arXiv：2007.12457。[2] Pinnau R，Thömmes G.玻璃冷却过程的最优边界控制数学方法应用科学2004;27（11）：1261-81. http://dx.doi.org/10.1002/mma.500.[3]Hinze M，Pinnau R.半导体设计的最佳控制方法。数学模型与方法应用科学2002;12（1）：89-107. http://dx.doi.org/10的网站。1142/S0218202502001568.[4]Blauth S，Leithäuser C，Pinnau R.微通道冷却系统的形状敏感性分析。数学分析应用杂志2020;492（2）：124476.网址：//dx.doi.org/10.1016/j.jmaa.2020.124476网站。[5]Schmidt S，Ilic C，Schulz VH，Gauger NR. 基于形状演算的三维大AIAA J2013;51（11）：2615-27。http://dx.doi.org/10.2514/1.J052245网站。[6][10]杨文，杨文.考虑非线性静磁效应的电机形状优化。SIAM J Sci Comput2015;37（6）：B1002-25. http://dx.doi.org/10.1137/15100477X，3436555.[7] Mitusch SK，Funke SW，Dokken JS. Dolfin-adjoint 2018.1：FeNiCS和Firedrake 的 自 动 伴 随 。 J Open Source Softw 2019;4 （ 38 ） ： 1292.http://dx.doi.org/10.21105/joss.01292网站。[8]Dokken JS，Mitusch SK，Funke SW。FEniCS和Firedrake中瞬态PDE2020，arXiv：2001.10058.[9]Paganini A，Wechsung F. Fireform：Firedrake的形状优化工具箱。2020年，arXiv：2005.07264。[10]Gangl P，Sturm K，Neunteufel M，Schöberl J.在NGSolve中完全和半自动化的形状区分Struct Multiplayer Optim 2020.网址：//dx.doi.org/10.1007/s00158-020-02742-w网站。[11]Hinze M，Pinnau R，Ulbrich M，Ulbrich S. PDE约束优化。数学建模：理论与 应 用 （ Mathematical Modeling ： Theory and Applications ） Springer ，New York; 2009，p. xii+270. http://dx.doi.org/10.1007/978-1-4020-8839-1.[12] 特勒茨湾偏微分方程的最优控制在：数学研究生，112，美国数学学会，普罗维登斯，RI;2010年，第112页。xvi+399.http://dx.doi.org/10.1090/gsm/112网站。[13] Delfour MC，Zolésio J-P.《几何学》，第二版。Advances in Design andControl ， 22 ， SocietyforIndustrialandAppliedMathematics（SIAM），Philadelphia，PA;2011，p.xxiv+622.http://dx.doi.org/10.1137/1的网站。9780898719826。[14]放大图片作者：Schulz VH，Siebenborn M，Welker K.基于Steklov-Poincaré型 优 化 的 PDE 约 束 形 状 优 化 。 SIAM J Optim 2016;26 （ 4 ） ： 2800-19.http://dx.doi.org/10.1137/15M1029369网站。[15]Blauth S.基于Steklov-Poincaré-Type的PDE约束形状优化的非线性共轭梯度法。2020，arXiv：2007.12891.[16]李志华，李志华，等.微分方程的有限元法求解.北京：清华大学出版社，1998。Springer; 2012，http://dx.doi.org/10.1007/978-3-642-23099-8。[17]Alnæs MS ， Blechta J ， Hake J ， Johansson A ， Kehlet B ， Logg A ，Richardson C，Ring J，Rognes ME，Wells GN. FeNiCS项目版本1.5。ArchNumer Softwa 2015;3 （100）. http://dx.doi.org/10.11588/ans.2015.100.20553网站。[18]刘伟，刘伟，刘伟.统一形式语言中的自动形状区分。Struct Multiplexp Optim2019;60（5）：1813-20. http://dx.doi.org/10.1007/s00158-019-02281-z网站。塞巴斯蒂安·布劳特软件X 13（2021）1006467[19]Balay S，Abhyankar S，Adams MF ，Brown J，Brune P ，Buschelman K，Dalcin L ，Dener A， Eijkhout V， Gropp WD ，Karpeyev D， Kaushik D，Knepley MG，May DA，McInnes LC，Mills RT，Munson T，Rupp K，SananP，Smith BF，Zampini S，Zhang H，Zhang H. PETSc用户手册。Tech.众议员ANL-95/11- 修订版3.13，阿贡国家实验室; 2020年，https://www.mcs.anl。gov/petsc.[20]Geuzaine C，Remacle J-F. Gmsh：一个内置预处理和后处理设施的三维有限元 网 格 生 成 器 。 Int J Numer Methods Eng 2009;79 （ 11 ） ： 1309-31.http://dx.doi.org/10.1002/nme.2579网站。[21][10]杨文，李文，李文.基于约束网格变形的一阶和二阶形状优化。SIAM J SciComput 2020;42 （ 2 ） ： A1200-25. http://dx.doi.org/10.1137/19M1241465 网站。[22] 作者 声明： John W.偏微 分方程 约束形 状优化的 曲面造 型计算 比较。ComputMethodsApplMath2016;16（3）：485-96.http://dx.doi.org/10.1515/cmam-2016-0009网站。