多学科设计优化笔记

1主要研究内容

（一）解耦方法

解耦方法是多学科设计优化的重要组成部分。解耦的目的是为了使整个复杂的多学科系统分解为便为处理的学科子系统。虽然分解后的这些学科子系统总的计算工作量甚至将大于系统在没有分解之前的工作量，但分解带来的好处是分解后的各个学科子系统可以在本学科内独立地进行本学科相关的分析与计算工作，学科内的计算不需要考虑与其它学科的耦合作用。这种分解优化策略还与当前的多处理器、并行计算的特点相一致，大大提高了优化的效率。

多学科设计优化的分解方法大多是借鉴运筹学中的分解原理。现已存在多种分解方法，它们都具有如下的特点：即通过系统分析、学科分析、系统灵敏度分析、学科灵敏度分析、系统层优化、学科层优化及学科的协同等使系统构成一个整体。其中，协同使分解后的系统保持了原系统内的耦合关系。

多学科设计优化的分解策略主要分为层次型、非层次型及混合型。在层次型系统中，整个系统形成一个金字塔形，数据从顶端的父层开始，因此下端的子层可接受上端父层的数据，而同层之间不能直接进行数据的交流。在非层次型系统

中，并没有父与子层的关系，故系统中数据可自由交换没有任何限制。由层次型系统与非层次型系统则可构成混合型系统。

多学科系统的分解可以采用原来的经验或例子，也可根据对系统中各个模块之间数据的流向而进行系统的分解或使用一些正规的方法如N-平方法可对多学科系统进行分解[21]。

（二）灵敏度分析

灵敏度分析的基本思想就是为了求得输入变量改变一个单位时输出变量的变化大小。在传统的单学科设计优化问题中，灵敏度分析主要用来计算最优点处优化点对目标函数或约束函数的灵敏度，常通过计算最优点处目标函数或约束函数的导数来实现灵敏度分析。采用的方法有有限微分法、分析法、自动微分法等。从原理上来说，多学科设计优化中采用的灵敏度分析方法与单学科设计优化相同，但由于多学科设计优化中常常具有多个学科分析模块，各个学科分析模块之间存在着耦合关系，即一个学科分析模块的输入变量可能又是另一个学科的输出变量，因此在进行灵敏度分析，计算学科输入变量对输出变量（本学科或其它学科的学科分析模块输出变量）的影响时通常要涉及到其它学科分析模块，传统的灵敏度计算方法已不能胜任多学科设计优化的要求。另一方面，多学科设计优化中灵敏度分析的变量维数也大大多于单学科设计优化问题，所需的数据量也大为的增加。传统的灵敏度分析方法不能简单地推广到多学科设计优化环境中，因此，在多学科设计优化过程中，使用一种新的灵敏度分析手段成为必然。

九十年代初，Sobieski导出了用于耦合系统灵敏度分析的全局灵敏度方程。GSE 是一种能有效计算相互耦合多学科灵敏度的方法，该方法直接从隐函数原理推导而来，精确性较高。Sobieski导出的算法有两种，一种是基于每个学科分析模块残数的导数；一种是基于各个学科分析模块输出对其输入的导数。其中，以第二种方法应用的最为广泛，首先，各个学科独立的完成各自的灵敏度计算，然后集成分析计算全局灵敏度。多学科设计优化的灵敏度分析的作用主要在于跟踪学科之间相互影响的功能，判断耦合性的大小。同时，在一些多学科设计优化方法中，通过全局灵敏度方程得到的全导数被用来对状态变量进行近似处理。

除了以上介绍的灵敏度分析方法外，另外一种灵敏度分析方法是采用神经网

络技术。通过神经网络模拟学科分析模块而实现灵敏度分析计算。神经网络的训练过程通常是在感兴趣的样本点上输出对输入连接权值的不断调整。因此，对于一个训练好的神经网络而言，其连接权值就可以认为是对敏感度的度量。换而言之，它可以被用于求在感兴趣的范围内系统输出对输入的敏感度，也是一种行之有效的方法。

（三）优化算法

在多学科设计优化环境中，优化算法参与优化的方式有两种。如果多学科设计优化的规划是采用单学科优化策略，则只需一个优化算法程序。如果多学科设计优化的规划是采用系统分解与协同的策略，则需要多种或多个优化算法程序。如系统分解后，在系统层和各学科层都需要配置一个优化算法程序，由于系统层和各学科层优化任务的不同，所采用的优化算法也可能不同。因此，在一个多学科设计优化问题中，可能需要综合应用几种算法，才能取得较好的效果。

多学科设计优化所采用的是已有的在单学科优化中所采用的成熟算法，根据问题的不同特点或性质可采用不同算法。如线性规划；以微分算法为基础的非线性规划：序列二次规划、广义既约梯度法、最速下降法、Newton 法、共扼梯度法、变尺度法、可行方向法、Powell 法、罚函数法与变步长下山法；遗传算法：自适应进化法，差分进化法，模拟退火法；随机搜索算法：Latin 超立方体法；整数规划等。到目前为止，在多学科设计优化的应用中，还没有发现一种对任何情况都适用的优化算法。但是对于特定的数学模型，总有一种或者某几种优化算法能取得较好的效果。对于同一个优化任务，也可以采用不同的算法进行计算及比较(因为算法本身有其局限性)来获得任务的最优解。

在这里需要指出的是，序列二次规划算法（SQP）是基于导数的算法中收敛最快的算法，并在优化之前的初始点不必可行的，因此，在多学科设计优化环境中，序列二次规划算法是采用得较为广泛的一种算法。

多学科设计优化算法进一步研究的重点在于开发出一些能解决设计全过程中出现的难解、完全NP、不可微非光滑等问题的高效的、对数学性态没有特殊要求的、具有并行处理特点的优化算法，以适应MDO 问题发展的需要。

（四）不确定性

在多学科设计优化过程中，如果不考虑不确定性因素的影响，则这时得到的优化解也可称之为确定性的优化解。通常，这些确定性的优化解处于设计空间的边界上。但考虑到实际工作环境中的复杂性，特别是各种确定性因素的影响，都将使确定性的优化解产生某些偏差，而更为重要的是将使某些约束条件失效。因此，如果要得到一个具有实际应用价值的优化解，则需要在进行多学科设计优化的同时，考虑如何减小不确定性因素的影响或者将不确定性的大小控制在一个合理的范围内。

由于多学科系统的耦合性，不确定因素在多学科系统中交叉传播，使不确定性的定量计算费用较大。如果伴随着多学科设计优化过程，计算费用将成倍增长。因此，需要根据问题的性质选择恰当的数学方法处理多学科系统中的不确定量大小。已有的数学方法有概率方法、模糊分析法与区间数学法等。但无论那一种数学方法，由于考虑不确定性因素的影响，多学科设计优化将变得复杂。为了减少计算费用，有时可以利用一些简化的手段或者近似处理方法。

（五）近似方法

由于多学科设计优化的复杂性，在某些情况下，直接采用学科模块中的学科分析参与多学科优化将带来许多问题，如将使优化的难度增大，效率低等。其主要表现在以下三个方面，1、对于设计变量较多的情况，优化算法调用目标函数与约束函数的次数将增多，相应地对学科分析的调用次数也将增多，导致的结果是计算量增大，优化时间变长。2、有时各个学科的分析模块驻留在分布于不同计算机或处于不同的地点，不利于将各个学科模块放置于一起进行多学科设计优化。3、某些学科模块分析结果可能有噪声或具有不连续性的特点，如果直接参与优化不仅引入了噪声因素，而且不利于采用梯度优化算法。

由于以上原因，在多学科设计优化中一般不将优化方法中的设计空间搜索程序与学科分析模块直接耦合，而采用了对目标函数与约束函数近似处理的方法。在优化过程中，目标函数或约束函数的近似模型直接参与优化，并经过反复迭代，使优化结果趋近并最终得到原目标函数的优化解。方法避免了在优化的搜索过程

中频繁地进行极其耗时的学科分析，从而大大降低了设计空间搜索工作的计算量。该方法类似于单学科优化的序列优化法。

在应用于工程系统的传统优化近似技术中，常采用基于导数的近似方法，如线性近似和二次近似，这些技术通常为局部近似。在多学科设计优化环境中，更多的是采用全局近似。目前，全局近似的一种形式是利用“可变复杂性模型"，也就是在优化过程中同时使用较简单与较复杂的模型，复杂模型为简单模型提供修正尺度，这里的简单模型就可以看作是全局近似。另一种特别适用于多学科设计优化问题的全局近似方法是响应面法。其中一种是用简单的函数(如多项式)来近似代替原函数，另一种实施途径是采用人工神经网络法。

MDO 在应用全局近似时，还需要应用试验设计技术进行数据的采样工作，如响应面法常用到的正交矩阵法。在这里，试验设计的目的是为了使用较少的计算量获得较高的近似精度。

（六）人机交互

多学科设计优化大大提高了解决复杂系统的优化效率，但整个设计过程并不能一步到位，自动完成。在设计过程中，还需要通过人机交互的作用进行人工干预。工程师将控制设计的流程、将设计者的判断与创造力融入到整个设计过程之中。因此，这种人机交互功能在多学科设计优化环境中是至关重要的。

人机交互功能的特点是能灵活地选择自变量与因变量，为了方便设计人员迅速获取这些信息，并能将这些信息以平面图、三维图、不同颜色的图形、正交投影图及动画形式图形等展示出来。其中，动画图形不仅可以显示振动系统中的动态变化行为，而且可以显示设计优化的迭代过程中设计的变化过程。当前，动画展示手段是采用虚拟现实的技术。以上对信息的展示主要局限于变量，而多学科设计优化环境中还需展示的是优化过程中的数据流。将数据流形象地展示给工程设计人员，将有利于设计人员设计思想的连续性，从而激发产生出许多的设计创新点。如Parkinson[22]就开发出一个在多学科设计优化环境中可显示优化中变量、约束、目标函数大小的交互窗口。

人机交互的另一方面是设计人员对优化进程的控制。简单的方法是菜单形式，而更高层次的控制手段则是以二次开发的形式实现对优化进程的控制。它将设计

人员从设计好的菜单式控制中解放出来，更有效地设置优化设计中的计算顺序，减少计算的费用。如Lee[23]在多学科优化环境Engineous中，采用二次开发与人工智能AI 相结合的技术操作与控制优化多学科优化进程。

（七）建模方法

物理系统通常是由多个子系统构成的一个整体，其各个子系统分别代表了系统的某个零件或局部的物理现象。多学科设计优化则是对构成系统的多个子系统进行协同优化以实现对整个系统的优化。因此，建立合适的子系统模型是多学科设计优化的前提。在多学科系统中，各个学科子系统模型是以学科分析模块形式来表现，例如：结构学科的有限元分析模块，气动学科的流体力学模块等。这些模块之间进行的数据交换代表了整个多学科系统的耦合情况。由各个学科分析模块及其模块之间的数据传输如实的反映了现实中的多学科系统的实际情况。因此，多学科优化模型的建立需要考虑的是学科分析模块的建立及如何处理模块之间的数据传输。

多学科系统中的数据传输在多学科优化中需要消耗一定的计算费用。它首先体现在模块之间传输的数据作为输入量进行计算之前都需要经过一次预处理转换。例如，气动学科的流体力学模块的输出量为位移，作为输入量传输给结构学科的有限元模块之前，需要将位移转换为负载量的形式参与有限元模块分析计算。数据传输的计算费用的另一方面是灵敏度方法的计算增大了模块间数据的传输量。针对数据传输所导致计算费用的增加，其处理方法有两种。一是对传输的数据进行近似处理，原来大量的数据传输变为近似模型的模型系数的传输。方法二是将模块间耦合强度高、数据传输量大的模块整合为一个模块来处理，从而部分减少了数据的传输量。

在建立多学科优化设计的数学模型时，需要强调的不能追求数学模型的精度而造成计算费用的大幅度上升，也不能为了减少计算费用而不顾一切地牺牲精确性。为此，采用如下方法来协调模块的计算精度与计算费用问题。首先，多学科优化模型中的学科分析模型与原学科分析模块的原理相同，但忽略了其中许多不重要的细节部分，得到精度较低的近似结果。其二，相对单学科优化中的学科分析模块，多学科优化不直接采用原来复杂的学科分析模块，而根据其它理论构造

计算效率高的分析模块。第三，在同一个学科中，还可建立学科分析模块的各种计算精度、计算费用不同的版本。例如，复杂的分析模块用来处理学科输出计算，简单的分析模块用来计算学科状态变量等。

近年来，随着全生命周期设计概念的提出，多学科优化系统的建立也增加了一些新的学科及其分析模块，这些学科及其模块分别代表了整个系统的经济性、可维护性、制造性与可操作性等，因此研究产品全生命周期建模技术，以适应MDO 的不同需要成为多学科设计建模技术研究的一项新任务。

（八）面向对象的分析

工程设计中的每一步都是以提出问题与解决问题的方式即所谓“What if”方式进行，为了在每一步中快速地解决问题，设计人员需要许多辅助的手段与工具，这些辅助的手段与工具可统称为面向对象分析。其主要包括有不同层次即从简单、近似到复杂、精确的分析方法的选择，设计变量改变时对相应部分进行的快速再分析，灵敏度分析，数据管理与可视化。其中，不同层次的分析将在数学模型中讨论，灵敏度分析将也将在后面作专门讨论，而数据管理与可视化内容较广泛，可参考文献[24]。

快速再分析的目的是使分析代码以一种较灵活的方式组织起来，当设计变量发生变化时，只需要计算尽可能少的代码就可得到分析结果。通常代码是以严格结构化的代码形式组成，每一组设计变量对应的问题将有不同的代码执行选项，如若直接采用结构化的代码组成形式进行快速预测再分析将很困难，而且因为如果对所有问题进行快速预测再分析数据的储存将是巨大的，因此采用严格结构化的代码形式进行问题结果预测是不可行的。

改进的方法是建立系统的模块化形式，构造模块之间、模块与系统外部输入的信息依存关系，这些信息依存关系将被保存与记录，当设计发生变化时，通过这些信息依存关系将找出受到影响的模块，给出为计算预测结果需要再次执行的模块。在进行快速再分析时，对于给出的设计变量的改变，首先需要检查是否保存有设计变量原来执行过的结果，如果有，则可直接得到结果而不需再次分析计算。否则，将利用信息依存关系得到需要进行再次分析计算的最少的一组模块及它们的执行循序，最后分析计算得到所需的结果。这种快速再分析方法将在设计

发生改变时，利用最小的计算量获得分析计算的结果。Kroo[25]与Gage[26]开发出了一种称为PASS 的快速再分析工具，它使用文档的形式保存信息依存关系。

2最优灵敏度分析方法

3滞后耦合伴随方法

4可变复杂度方法

最初的可变复杂度方法（Variable-Complexity Method，简称VCM）是用来发展基于梯度的优化技术，整个设计过程由一组优化循环组成。在各优化循环开始的时候采用精确分析，而近似分析则运用于各循环内的后续计算。精确分析获得一个比例因子，用这个比例因子修正近似分析，而且这个比例因子可以在设计过程中不断更新(Dudley et al., 1995; Hutchison et al., 1994)。VCM 的基本流程为(Alexandrov and Lewis, 2002b)：

(1) 分别应用高精度模型和低精度模型计算初始点X0处的目标函数值和约束值；

(2) 求解比例因子；

(3) 使用比例修正后的低精度模型进行优化搜索；

(4) 对优化结果应用高精度模型进行评估；

的结合。在高速民用运输机（High-Speed Civil Transport，简称HSCT）的设计中，成功地应用了该方法。在研究中，首先对由全因子试验设计获得1296 个方案进行初步可行性分析，去除明显不合理的设计，余下157 个设计方案。然后，对这些方案进行低精度模型的分析，建立二次多项式响应面，并根据方差分析（Analysis of Variance, 简称ANOVA）去掉对响应值影响小的项；再根据D-最优设计（D-optimal point selection）选择50 个方案进行高精度分析，并以其结果建立响应面模型（采用方差分析后的响应面形式），以此作为优化对象，该方法取得了良好的结果(Giunta et al., 1995)。

5响应面更新策略

5.1加实验点

5.2Interval Reduction

Once a system optimization is conducted, it is important to check the accuracy ofthe response surface fit. A simple calculation was used to compare values taken from theresponse surfaces to those produced directly from the BB analyses. First, the values {Z |Y | w}* are found by the system optimization. Then these values were used as BB inputparameters, and the BB’s were re-optimized. The difference YBB*-YRS* was taken andnormalized by YBB* to give a measure of the “goodness of fit” obtained in the quadraticrespons e surface model.

The measure of error due to lack of fit was used as a termination criteria. If, at theend of one iteration, the values YRS* and YBB* do not agree by a predetermined tolerance,the response surfaces are reformulated by reducing the intervals around

system optimalvalues.

Reducing the interval sizes has the general property of improving the fit of theresponse surface. Indeed, as the intervals are reduced, the fit becomes more and morelinear. It is desired that the new interval be centered on the anticipated optimum value,and that the new interval makes physical sense. Since the Y and w variables are notsubject to any clear physical limitations, their interval bounds are allowed to “drift”across the original interval bounds. However, the Z variab les are subject to such physicallimitations so their intervals must always lie within the original interval bounds. Thesedifferences lead to the use of two different interval reduction methods.

Intervals on system level variables Z are reduced in such a way as to ensure thesubsequent interval would fall within the physical bounds described by the originalinterval. For example, let XLi and XHi denote the respective lower and upper limits onthe interval used to generate a response surface. If the termination tolerance is notsatisfied, then the interval is reduced by 20% and centered on the system optimal value.

If any of the new XLi+1is less than XLi , the value is returned to XLi. Similarly, if any ofthe XHi+1 is larger than XHi the interval limit is returned to XHi. This operation shrinksthe interval over which points are placed to generate the response surfaces. It should benoted, however, that reducing the intervals in this manner does not ensure that the newinterval is centered on the optimum found in the previous iteration.

The factor by which the interval is reduced each iteration was chosen somewhatarbitrarily. Since the interval is reduced each iteration (without the possibility of beingincreased), it is important to reduce the interval size gradually to prevent prematureconvergence to an unreliable optimum.

The intervals for coupling variables Y and weight factors w were not subject tothe limitation of falling within the original interval. Therefore, the intervals are allowedto shift solely on the basis of where the optimum was found on the previous iteration.

This carries the added benefit of centering the RS’s in the middle of the design

space withrespect to Y and w, so satisfaction of the coupling constraints is always possible.

A simple algorithm was adapted from [16] for re-sizing the intervals on Z and w.

First, a factor K is defined for reducing the interval size. This factor is not directlyapplied to the current interval. Instead, an adjustment factor A i is determined by linearinterpolation as:

whereC i is the center of the current interval, Opti is the optimum found in the currentinterval, and Li is the lower value of the current interval. The actual reduction factor canthen be found as:

The result of this process is that the interval size will not be reduced if the optimum pointis located on a boundary, while it will be reduced by a factor of K when the optimumpoint is located at the center.

Additionally, a history of optimum location with respect to bounds is kept. If thebest design point hits the same bound twice, the interval is expanded by a fixed amount.If the optimum oscillates from upper to lower bounds between iterations, the interval isreduced by a fixed amount.

5.3Modified Interval Reduction（自创）

此方法与original interval reduction类似，只是when we follow the routine of original one to reduce the interval, number of existed experiment point may be reduced to the one not enough to perform RSM. 因而在reducing the interval时，bring in an out-screening mechanism, namely, 在确定在reduced interval中存在足够的实验点时按照original interval reduction进行操作，否则在以上一cycle的优化点为中心的情况下，选择能保证存在足够实验点的最小interval，遇到边界时

将interval的另一边界向相反的方向移动以保证存在足够的点。

5.4实验点选取方法

1）选择上一cycle的解点即最优点

2）类似SKO或MSKO的取点方法，取得全局最优点利用max( expented-improvement)

5.5Move limit（Move limits definition in structural

optimizationwith sequential linear

programming.Part I: Optimization algorithm）

Since Pope introduced move limits in the early 70s,researchers proposed a variety of techniques to definethe move limits. Haftka and Gurdal [1] suggests tochoose as move limit the 10–30% of the value that adesign variable has at the beginning of the iterationcycle. The move limits are gradually shrunk as the designapproaches the optimum and they should be reducedif the cost function does not improve or theconstraints are more violated than at the previous iteration.Vanderplaats and Kodyalam [6] adopt the followingcriterion: the move limits are very large in thefirst design cycles and they are then reduced by 50% ifthe design violates the constraints more than at theprevious iteration; anyhow the move limits are neverreduced to less than 25%. John et al. [7] solve the linearizedsub-problems with the Simplex algorithm. Movelimits are reduced by means of the a parameter each timethe design does not improve in the current iterate. Theoptimal step size a is determined with reanalysis techniquesbased on quadratic approximations of the costfunction.

原则：They concluded that a proper choice of themove limits should ensure that the objective functionimproves monotonously, that each intermediate solutionis feasible, that the design variable movement is controlledin order to maintain the approximation error at areasonable level.

It can be seen from Fig. 1 that the value of the trustregion parameter q gets smaller than one when the designvariables are perturbed so that the cost functionimproves. In such a case, the linear model computes thesmallest value for the trust region parameter while highorder approximation models result in smoothing thevariations of the q parameter. Therefore, the linearmodel may be considered the most conservative schemefor the purpose of shrinking move limits.

6MDO基本理论（飞行器ＭＤＯ过程）6.1定义

6.2基本优化组织模式

6.3单级优化过程