CN110705150B - 一类工程结构超高维大规模多约束非线性优化方法 - Google Patents

Abstract

本发明公开了一类工程结构超高维大规模多约束非线性优化方法,利用非线性优化迭代过程中各子问题的纯量目标函数二阶导具有对称矩阵的特性，在Newton‑Krylov迭代算法的基础上，提出了一种新的正交化算法以及对其负定元素的校正处理算法，该方法的线性方程组系数矩阵仅为一对称三对角矩阵，且由于改进了迭代方向的稳健性，不仅节省了系数矩阵及其后续计算的工作量，且存储量、计算效率及其稳健性得以明显改进。

Description

一类工程结构超高维大规模多约束非线性优化方法

技术领域

本发明涉及轻型工程结构技术领域，尤其涉及一类工程结构超高维大规模多约束非线性优化方法。

背景技术

在轻型工程结构设计领域，如航空航天工程装备的壳体及内部骨架结构，结构减重或者说构件的截面尺寸优化是具有重要意义的设计工作任务。减轻重量意味着材料用量减少，这又意味着装备的能效比提高以及成本的降低。结构重量减轻时，结构构件中应力、位移等响应量呈增大趋势，同时影响结构整体的力学行为特性。因此，控制结构构件中的应力、位移等的响应量水平以及结构整体行为的特性，并极大减轻结构重量是保证结构安全使用且提高装备能效比的一项核心关键技术。

基于数值模型驱动的设计方法是现代结构工程中广泛应用的技术手段，也是提高工程结构设计品质及其精细化水平的基本技术途径。随着不断提高的计算机硬件能力，结构有限元模型的精细化程度亦随之提高，这使得设置结构构件截面的可设计变量具有了更多的灵活性和操作性，现行工程结构的设计变量数目已可达10³以上的量级。为获得安全有效的高品质轻量化结构，工程结构设计中的性能约束控制亦向精细化发展，约束规模数也可达10³量级或更多。以航空航天工程装备的壳体结构为例，结构优化模型中的设计变量通常为蒙皮或骨架构件的截面厚度或截面形状的几何变量，通常结构按总体内力水平分成数十或数百个区域，一个局部区域设置一个变量或多个变量；同时在局部区域上控制其应力或应变的强度安全水平。除结构局域的内力约束控制之外，还包括受压或剪切内力区域的结构稳定性控制水平；同时，结构尚需满足整体力学性能，这包括结构体的变形挠度、区段截面的扭转角以及避免操纵控制失效的其他刚度条件和结构的动力学固有品质特性等。

针对超高维变量与大规模多性态性能约束的轻重量结构数值优化设计问题，基于模型数值梯度的优化算法是现代公认的一类高效率工程实用算法技术。然而，超高维变量的大规模多性态力学性能与其结构重量函数的数值梯度计算，以及在此基础之上的非线性数值优化迭代仍然是这类算法中一项极具挑战的研究课题，解决优化迭代效率、存储量、计算量以及算法稳健性之间的矛盾，一直是学术界极其活跃的研究方向。文献“Iterativemethods for linear and nonlinear equations”中公开了一个基于伽辽金数值计算原理的Arnoldi方法及其Newton-Krylov迭代算法，该算法在大规模非线性矢值函数梯度计算以及非线性迭代求解方面是学术界近年来的一类标志性算法，其优点在于利用二阶导矩阵与向量的乘积形式，采用高维方向导数的近似解法解决了各Newton步中的二阶导计算问题，并利用子空间寻优方法初步解决了计算中的存储量问题。但缺点是在生成Krylov子空间的迭代过程中需存储Gram-Schmidt正交化算法的各相关系数以及存在正交化基向量的重复运算等问题。对数值优化问题而言，随设计变量规模的增加，这类算法在求解序列子目标优化问题中的计算量呈几何级数增长，从而未能有效解决好算法的计算量以及存储问题。

发明内容

针对上述问题，本发明提供了一种工程结构超高维大规模多约束非线性优化方法。

为实现上述目的，本发明采用以下技术方案：

一种工程结构超高维大规模多约束非线性优化方法，包括如下步骤：

1)构建结构力学数值分析的有限元模型，置优化迭代计数k＝0，设置设计变量初值x₀及优化计算收敛精度ε，将设计变量初值x₀赋给需优化设计的结构件，其中，设计变量包括：工程结构单元数量以及几何属性参数组；

2)计算分析结构的性能约束函数、结构的重量目标函数以及所述性能约束函数和结构的重量目标函数各自关于设计变量的一阶梯度向量；

3)利用拉格朗日乘子法的增广形式,将所述结构的重量目标函数与性能约束函数组合成无约束化的子问题目标函数Φ_k；

4)对子问题的无约束化目标函数在给定设计变量点处进行线化处理，得到子问题无约束化目标函数的线化方程组系数矩阵，选择子空间的维数m，采用Krylov子空间迭代方法计算该线化方程组系数矩阵在子空间中的最优解；

5)根据所述线化方程组系数矩阵的对称性，对无约束化目标函数的线化方程组系数矩阵实施非Gram-Schmidt算法的正交化改造，获得非完全正交化的标准基向量组以及Krylov子空间最优解计算所需的矩阵分解式；

6)对矩阵分解式中的三对角对称矩阵进行特征值谱分解，进而代入子空间最优解方程，计算得子空间标准正交基向量组合系数的近似解基向量组合系数；

7)对所述的三对角对称矩阵的特征值及组合系数的近似解实施稳健性校正处理，即逐一检查其所有特征值，当特征值小于零时，将相应的组合系数解反号；当特征值等于零时，将相应的组合系数为零，并计算线性搜索方向及线性搜索设计变量值的充分解以及残量；

8)判断所述残量是否小于或等于预设的迭代误差精度值：是，则转至步骤9)；否，则将迭代过程计算得到的一组设计变量值x_k替代所述设计变量初始值，优化迭代计数k＝k+1,并转至步骤2)开始新一轮优化迭代计算；

9)判断所得设计变量值是否达到优化收敛精度ε：是，则结束优化；否，则迭代计数k加1，将迭代过程计算得到的一组设计变量值替代前次设计变量值，修改需优化的结构件参数，转至步骤2)继续迭代优化。

进一步，所述步骤5)中的非Gram-Schmidt算法的正交化改造，包括如下步骤:

2-1)取q₁＝r⁰/||r⁰||₂，计算α₁＝(q₁,A q₁)，r¹＝A q₁-α₁q₁，取β₁＝||r¹||₂，q₂＝r¹/||r¹||₂；

2-2)取j≧2，计算α_j＝(q_j,A q_j)，r^j＝A q_j-α_j q_j-β_j-1q_j-1，β_j＝||r^j||₂，q_j+1＝r^j/||r^j||₂；

2-3)当j＝m+1时，结束计算；否则，j＝j+1，返回步骤2-2)；

上式中，r⁰为第k次迭代求解中计算线化方程组的初始残余向量，||·||₂表示模运算，q₁为单位残余向量，A为线化方程组系数构成的矩阵，Aq₁表示矩阵与向量的乘积运算，(·,·)表示两向量的点积运算，α₁即为向量点积运算数值，β₁为向量的模；m为选定的基向量数，上标j表示第j次循环计算，r¹,...,r^j,...,r^m为所提算法循环计算得到的第1至第m个残余向量；

上述3个计算步骤即得非完全正交化的标准基向量组

并由理论可知所述计算步骤可写作以下矩阵分解式：

其中，U_n×m＝(q₁,q₂,…,q_m)，

称为m阶三对角阵，

U_n×(m+1)＝(q₁,q₂,…,q_m,q_m+1)。

进一步，所述步骤6)中的三对角对称矩阵特征值谱分解以及最优近似解基向量组合系数，包括如下步骤:

3-1)利用三对角矩阵的对称性，对三对角矩阵T_m×m进行特征值谱分解为：

T_m×m＝Q·S·Q^T,S＝Diag(λ₁,…,λ_m)

其中，Q为单位上三角矩阵，S为T_m×m的特征值对角矩阵，λ_j，j＝1,…，m

为T_m×m阵的特征值，上标T表示矩阵的转置，diag表示对角阵；

3-2)将T_m×m代入前述矩阵分解式，对子空间组合系数最优解予以近似表达，得到如下解耦方程：

式中，y′_m＝Q^Ty_m，β＝||r⁰||₂，

为

去掉最后一个元素，

3-3)计算得子空间组合系数最优解的近似值：

进一步，所述步骤7)中的稳健性校正算法处理及充分解，包括如下步骤：

4-1)当特征值λ_j<0，将相应的组合系数解反号，即修改y′_j＝-y′_j；

4-2)当特征值λ_j＝0，相应的组合系数y′_j＝0；

4-3)计算线性搜索方向向量s_k＝U_n×mQ_m×my′_m，并置线性搜索步长α＝1，

4-4)置设计变量解x_k+1＝x_k+αs_k，x_k为本轮子问题设计变量的初始值；

4-5)计算子问题本次搜索的函数值Φ(x_k+1)，若Φ(x_k+1)<Φ(x_k)则结束；否则，置α＝tα，t＝0.85，返回前步，直至满足，并得充分解x_k。

本发明有益效果为：

本发明利用非线性优化迭代过程中各子问题的纯量目标函数二阶导具有对称矩阵的特性，在Newton-Krylov迭代算法的基础上，提出了一种新的正交化算法以及对其负定元素的校正处理算法，该方法的线性方程组系数矩阵仅为一对称三对角矩阵，且由于改进了迭代方向的稳健性，不仅极大节省了系数矩阵及其后续计算的工作量，且存储量、计算效率及其稳健性得以明显改进。

附图说明

图1是本发明一类工程结构超高维大规模多约束非线性优化方法的流程图；

图2是工程均质阶梯梁在满应力约束下的重量优化模型图。

具体实施方式

本发明使用的数值优化模型在方法总体上属于一类无约束化模型，即将一个多约束条件下的单目标函数数值优化问题利用增广的Lagrange方法，将原约束方程与原目标函数在每个非线性迭代计算过程中组合成一个称之为无约束化的子目标非线性纯量函数的极值问题求解。在每个子问题的无约束化非线性函数求解中，本发明采用Newton步的子空间迭代方法。本发明一类工程结构超高维大规模多约束非线性优化方法的基本特征在于：改变了每个Newton步方程迭代过程中形成正交基向量的计算策略，从而得到工作效率上的极大提升，明显降低了计算工作量和过程中的存储量，又因调整了迭代方向的稳健性，对数值优化结果有明显改观，为轻型结构超高维大规模约束的非线性优化设计问题提供了新的算法技术途径。具体包括以下步骤：

步骤1:工程结构超高维大规模多约束非线性优化问题原型。如前所述，轻量化结构的目标函数是获取其极小化重量值以及使其极小化的设计变量取值，重量涉及到结构中各类构件截面尺寸参数的非线性运算及线性累加。因此，工程结构优化的结构重量函数可记作：

式中，x＝[x₁,…,x_n]^T为设计变量列向量，指结构需优化的各构件截面几何参数，对超高维问题，通常n>100；f为结构体需优化计算的总重量，w_i为第i个优化设计的构件重量，M为需优化的结构构件数；Arg为提取函数变量值的算子符号。

结构优化设计中的约束总体上分为两类：一类称为工艺约束，或尺寸限约束，是指每个设计变量在优化过程中的可变动范围，显然，其下限必须大于零，上限原则上可以任意，工程中可明确判断结构构件截面上所允许的最大尺寸，太大会失去工程使用意义，通常记作：

另一类约束即为结构各类力学性能约束，如前所述，包括结构整体的静动力学特性要求，或结构局部乃至关键点的强度安全控制限等。性能约束是指性能函数在设计变量上的取值与限制界限间的非严格不等式关系，工程上对此有清晰的量化要求。性能函数是通过结构体力学状态方程组联立计算的，除极其简单情况可给出解析形式外，通常只能在一组设计变量取值条件下计算其离散值，这表明力学性能函数是数学意义上的隐函数。性能约束的非严格不等式在数学优化模型以及算法中的意义不大，故用严格不等式表达工程约束要求即可。另外，工程中个别特例会要求等式约束，也是数学模型中常见的约束类型。由此，结构力学上的各类性能约束条件可表示为：

大规模多约束通常指m_e+m_i1+m_i2是一个较大的数。数学优化模型通常采用严格小于和等于零的单边规范形式来表达前述约束，即：

上式中的三个分母项用于消除不同力学约束间不必要的量纲效应，δ用于保证分母非零，一般取0.001即可。

于是，工程结构超高维大规模多约束非线性优化问题的数学符号原型可记作：

式中，m_i＝m_i1+m_i2，I_me和I_mi为等式和不等式的指标集。

步骤2：非线性优化子问题的无约束化目标函数。将初始给定或迭代过程计算得到的一组设计变量值代入结构体状态方程求解重量函数、性能约束函数及其一阶导列向量函数的数值，即意味着一个优化子问题的开始。对于复杂的非线性数学优化问题，通常需要对式(5)的原型做简化处理，以更有效的方式求解一个子问题的最优解。本发明采用了一类无约束化子问题的处理方法，即利用Lagrange乘子的增广形式将原目标与约束综合成一个纯量的非线性函数，称作子问题的无约束化目标函数，或简称子问题目标函数。该类简化方法的最大优点在于不改变原目标函数及性能约束函数的非线性性质，且不增加优化变量，其形式写作：

式中，下标k表示第k个子问题，x_k即为该子问题的起始设计变量值，鉴于各子问题的目标函数形式均一致，故后文中略去子问题的下标标记；

为Lagrange乘子列向量；

为罚参数列向量；z＝{z₁,…z_mi}为人工松弛变量列。

为降低子问题寻优计算的难度，对式(6)可实施以下预处理：其一，每个子问题在起始点迭代计算前预设{λ,σ}的取值，并在本次计算中保持不变，直至下一个子问题开始前更新之；另外，σ及式(6)中第二个和式的λ需限定在大于零的实数范围内。其二，运用迭代寻优的基本原理将松弛变量列z用不等式约束的相应关系替代，且由此建立对有效约束的部分筛选原则。寻优计算的基本原理是获取目标函数关于各变量梯度为零的解。由此，令

可得：

对上式讨论可知，仅当发生λ^j+σ^jg_j(x)<0时，z_j≠0，将其情况代入式(6)可得该式的第二个和式结果为常数。优化函数中包含常数对优化迭代过程及其结果均无意义，故在子问题目标函数迭代中用

予以有效约束的部分筛选。

步骤3：大规模多约束非线性目标函数的Newton迭代运算框架。如前所述，对于一个大规模复杂多约束非线性目标函数的寻优求解，将其转化为一个个子问题的无约束化目标函数求解。求解一个子问题的目标函数是指求解非线性方程组

在给定起始点x_k的线化形式，或称线化方程。Newton法给出在点x_k处的线化形式为

式中，[Φ″(x_k)]为Φ(x)关于设计变量在点x_k处的二阶导对称方阵，又称为Hessian阵，随设计变量维数的增大，其规模很大。对工程结构问题而言，严重的问题是该阵难以从结构体状态方程中求解，即使可能，其工作量也无法承受；s_k为线化形式的解变量，也是优化问题设计变量解过程的第k次增量；

为子问题目标函数在点x_k处的梯度向量，不难由结构体状态方程直接求解。

获得子问题目标函数线化形式的解s_k则意味着求解了一个Newton步，并将

x_k+1＝x_k+s_k (9)

作为新的起始点，进入下一个子问题目标函数的构造与求解。该过程称为优化问题的外迭代过程。

步骤4：子问题目标函数线化方程的Krylov子空间迭代解过程。鉴于前述Hessian阵的计算难度，工程中常采用迭代方法求解式(8)。Krylov子空间迭代方法是近年来发展的一类高效求解线性方程组的优秀迭代算法，该算法基于伽辽金广义正交原理。

设x＊为一线性方程组Ax＝b的真解，任意给定一个起始点x₀，称z＝x＊-x₀为距离真解的误差向量，或简称误差。在误差式两端同乘A可得：

Az＝Ax＊-Ax₀＝b-Ax₀＝r₀ (10)

式中，r₀称为残余向量或残量。上式(10)与原方程式(8)等价，并表明了右端项及未知量的意义。若A阵已知且数值特性良好，则理论上可解，但工程中难以直接计算。另外，Ax＝b是式(8)的一般记法，后文中必要时用A来替代[Φ″(x_k)]，即A意指对称阵。

Krylov子空间迭代算法是在一m维子空间K_m上寻找与r₀相差最小的z_m，其中，0<m<<n，n为矩阵A的阶次。该原理可表示为：

式中，

指向量2-范数的平方。Krylov子空间迭代算法的主要过程如下：

1.设置m的整数取值，通常m<100；设置迭代误差精度ε,通常取ε＝1.0×10^-6；

2.构造子空间K_m的一组基

记K_m＝Span{r⁰,A r⁰,A²r⁰,…,A^m-1r⁰}；

3.按Gram-Schmidt步骤将

转换为K_m的一组标准正交基

即K_m＝Span{v₁,v₂,…,v_m}，且(v_i，v_j)＝δ_ij，其中，(·,·)表示向量点积；

4.将K_m上的任一向量表达成该标准正交基向量的线性组合，即z＝V_n×my_m，

其中，y_m＝[y₁,…,y_m]^T为线性组合系数，是需求解的待定系数向量；

5.将z的表达式代入式(11)计算得线性组合系数y_m及x_m＝x₀+V_n×m y_m；

6.计算||r_m||＝||b-Ax_m||，若||r_m||<ε则结束；否则，置x₀＝x_m返回第2步重新开始。此过程称为优化求解的内迭代过程。

上述Krylov子空间迭代算法中有两个关键:

的计算和第3步的计算效率。

计算

中的各项可写作更一般的形式，即计算A·ρ，其中||ρ||＝1。本发明仅对式(8)给出A·ρ的计算方法，即A＝[Φ″(x_k)]。由高维方向导数的定义可知：

式中，ε应取足够小实数，则解的近似性就足够好。上式的近似计算即解决了工程中A矩阵的计算难题，且存储量也极大减少。

运用上述第3步的算法结果，推导可得如下矩阵分解式为：

式中，

其中，h_ij＝(Av_j,v_i)，且(v_j,v_i)＝δ_ij，为循环计算所得。

需要说明的是，Gram-Schmidt正交化的计算量约为(m+1)(m+2)m·n/3，当m＝10²和n＝10³，计算量约为3亿次，计算量惊人。

步骤5：本发明的不完全三对角正交化过程。借助于式(8)系数矩阵[Φ″(x_k)]的对称性，并暂假设该矩阵正定。据此，本发明对

的Gram-Schmidt正交化过程进行了算法改造，算法描述如下：

1.取q₁＝r⁰/||r⁰||₂，计算α₁＝(q₁,A q₁)，r₁＝A q₁-α₁q₁，取β₁＝||r¹||₂，q₂＝r¹/||r¹||₂；

2.取j≧2，计算α_j＝(q_j,A q_j)，r^j＝A q_j-α_j q_j-β_j-1q_j-1，β_j＝||r^j||₂，q_j+1＝r^j/||r^j||₂；

3.当j＝m+1时，结束计算；否则，j＝j+1，返回第2步。

上述算法可得标准正交基向量组

同理，上述算法可表示成以下矩阵分解式：

其中，U_n×m＝(q₁,q₂,…,q_m)，

称为m阶三对角阵，

β_m＝||r^m||₂q_m+₁＝r^m/||r^m||₂,

U_n×(m+1)＝(q₁,q₂,…,q_m,q_m+1)；

注意，式(14)的计算量约为(5m-1)·n，当m＝10²和n＝10³，其计算量约为50万次，较现有Gram-Schmidt的正交化步骤计算量减少3个数量级。

步骤6：式(11)的快速近似解算法。式(14)的正交化算法使本发明后续求解式(11)的算法更为简洁。取K_m＝Span{q₁,q₂,…,q_m}＝Span{U_n×m},该子空间上的任一向量表示为z＝U_n×my_m，y_m为待定系数向量。代入式(11)可得：

式中，β＝||r⁰||₂，

找上式最小值即等于求解冗余线性方程组：

即：

求解式(16)的算法第1步是采用追赶法将下对角行元素消为零，即得：

其中，α′₁＝α₁，q′₀＝β,

i＝2,…,m

第2步是消除上式的冗余性并进行高效高精度近似求解，即去掉系数矩阵的最后一行元素以及右端项的最后一个元素，使其为一适定线性方程组。对消除冗余后的线性方程组实施由下向上的回代求解即可得：

并可得残量

以及x_m＝x₀+U_n×my_m。

若||r_m||₂<ε则结束计算；否则，置x₀＝x_m返回步骤5，开始新一轮计算。

步骤7：子问题目标函数优化近似解的射线步处理算法。实际上，本发明的前述步骤5和6所得解对式(8)而言是近似的，原因之一在于残量||r_m||₂经历了大量运算，由于计算机精度有限，舍入误差污染严重；另外，Aq_i也是近似计算，且式(18)所得解s_k＝x₀-x_m是子问题目标函数Φ_k下降的间接近似解。为此，需要对解s_k再做处理，以保证本次子问题优化迭代的充分性。本发明采用常规的简单射线步做法：s_k＝U_n×mQ_m×my′_m

1.置s_k＝U_n×m y_m，α＝1，其中，下标k为第k个子问题，m为近似解空间维数；

2.x_k+1＝x_k+αs_k；

3.计算子问题本次搜索的函数值Φ(x_k+1)，若Φ(x_k+1)<Φ(x_k)则结束；否则，置α＝tα，t＝0.85，返回前步。

步骤8：子问题目标函数迭代稳健性校正处理算法。实际上，前文算法的一个前提是要求子问题目标函数在x_k点处沿任意方向必须是凸的，才能保证[Φ″(x_k)]的半正定性，继而算法得以顺利执行，否则可能导致算法的中断或无解。为解决此类求解的稳健性，本发明提出了以下校正算法，以保证有效执行并得其解。算法原理在于关注到式(14)中三对角矩阵T_m×m的对称性，这可利用对称矩阵的特征值分解算法将其变换为：

T_m×m＝Q·S·Q^T,S＝Diag(λ₁,…,λ_m) (19)

其中，S为T_m×m的特征值，但对角元中可能存在零和负值；Q为单位上三角阵，即：

于是，去掉式(16)最后一行，并将其中T_m×m改写为式(19)，即可得其解耦形式：

式中，y′_m＝Q^Ty_m，β＝||r⁰||₂，

为

去掉最后一个元素。

求解过程中，当λ_j<0，修改y′_j＝-y′_j；当λ_j＝0，y′_j＝0。

置s_k＝U_n×mQ_m×my′_m，α＝1，返回步骤7中的循环计算，从第2步开始，直至得到本次子问题目标函数的充分极值解。

本发明利用非线性优化迭代过程中各子问题的纯量目标函数二阶导具有对称矩阵的特性，在Newton-Krylov迭代算法的基础上，提出了一种新的正交化算法以及对其负定元素的校正处理算法，该方法的线性方程组系数矩阵仅为一对称三对角矩阵，且由于改进了迭代方向的稳健性，不仅极大节省了系数矩阵及其后续计算的工作量，且存储量、计算效率及其稳健性得以明显改进。

本发明的完整迭代优化流程如图1所示。关于本发明的特点及效果，下面将结合算例、附图和实施例做进一步说明。

实施例

如图1所示，本实施例1为10个无约束优化问题的国际标准测试算例，相当于本发明无约束化优化子问题的复杂显函数，设计变量数为100。具体过程包括以下步骤：

步骤1：按本发明的算法框图编制软件程序。软件设计中，需编制10个显形式函数及其梯度计算的子例程，用于非线性数值优化迭代过程中调用。这包括：

1.函数名称：Extended Rosenbrock

F_2i(x)＝1-x_2i-1

i＝1，2……50；初值点：-1.2，1，1.2，1，……。

2.函数名称：Augmented Rosenbrock

F_4i-2(x)＝1-x_4i-3

F_4i(x)＝x_4i

i＝1，2……25；初值点：-1.2，1，-1，20，-1.2，1，-1，20，……。

3.函数名称：Modified Rosenbrock

i＝1,2……50；初值点：-1.8，-1，-1.8，-1，……。

4.函数名称：Augmented Powell badly scaled

F_3i-2(x)＝10⁴x_3i-2x_3i-1-1

F_3i-1(x)＝exp(-x_3i-2)+exp(-x_3i-1)-1.0001

i＝1,2……34；初值点：0，1，-4，0，1，-4，……。

5.函数名称：Tridimensional vally

F_3i-1(x)＝10(sin(x_3i-2)-x_3i-1)

F_3i(x)＝10(cos(x_3i-2)-x_3i)

c₁＝1.003344481605351,c₂＝-3.344481605351171×10^-3

i＝1，2……34；初值点：-4，1，2，-4，1，2，……。

6.函数名称：Shifted and augmented trigonometric function with anEuclidean sphere

i＝1，2……99；初值点：0，0，……。

7.函数名称：Trigonometric-exponential system,Trig exp 1

F₁₀₀＝4x₁₀₀-x₉₉ exp(x₉₉-x₁₀₀)-3

i＝1，2……50；初值点：0，0，……。

8.函数名称：Singular Broyden problem

F_2i-1(x)＝[(3-2x_2i-1)x_2i-1-2x_2i+1]²

F_2i(x)＝[(3-2x_2i)x_2i-x_2i-1-2x_2i+1]²

F₁₀₀(x)＝[(3-2x₁₀₀)x₁₀₀-x₉₉+1]²

i＝1，2……50；初值点：-1，-1，……。

9.函数名称：Broyden tridiagonal function

F_2i-1(x)＝(0.5x_2i-1-3)x_2i-1+2x_2i-1

F_2i(x)＝(0.5x_2i-3)x_2i+x_2i-1+2x_2i-1

F₁₀₀(x)＝(0.5x₁₀₀-3)x₁₀₀+x₉₉-1

i＝1，2……50；初值点：-1，-1，……。

10.函数名称：Tridiagonal system

i＝1，2……50；初值点：12，12，……。

上述国际标准测试算例均以矢值函数形式给出，为将其转换为标准纯量函数，即本发明算法处理的子问题显形式，采用以下公式转换：

F(x)为由前述各矢值函数组成的列向量。

软件程序中，再增加非线性优化迭代收敛的判据如下：

1.正常收敛条件：

x*为理论解；

2.非正常条件：|Φ_k-Φ_k+1|/|Φ_k|≤10^-6；或迭代次数k>300。

步骤2：因前述函数均为超越函数或多项式的组合形式，故其二阶导数有显形式解，为对比本发明算法的优势，本实施例除采用本发明算法外，还利用二阶导的显形式直接构造Newton精确迭代算法，即本发明中的步骤3解算过程。

步骤3：由前述10个显函数的纯量形式经计算可知，本实施例的最小纯量函数值全部为零。表1列出了本实施例的全部数值优化迭代计算结果，其中，NN为优化迭代计算次数，FE为优化结果。注：表中横杠表述计算不满足收敛条件而中断。

表1本发明算法与Newton算法的结果对比一览表

由本实施例步骤完成的结果可以得出以下基本结论：

1.本发明的新算法较常规Newton算法的稳健性明显提高。本发明算法只有函数5在规定的迭代次数内未找到有效解；函数7陷入了局部极值。而常规Newton法，仅得到了4个理论解，函数1、3和5在规定的迭代次数内未能找到有效解，函数6、7和8陷入了局部极值。

2.本发明的新算法较常规Newton法收敛效率更高。虽函数6和7的本发明算法迭代次数超过常规Newton法，但本发明算法得到了函数6的理论解；而函数7的极值也远小于常规Newton法的解。

实施例2

如图1、2所示，本实施例2为一工程均质材料的矩形截面阶梯梁，阶梯梁长度L＝1000mm，弹性模量E＝70GPa，泊松比ν＝0.3，密度ρ＝1，梁的左端面为固支撑，右端面施加集中剪力P＝625N，如图2所示。为检验本发明算法结构轻重量优化的迭代稳健性及其效率，阶梯梁沿长度分成不同段数n，以变动设计变量总数，阶梯梁横截面的宽和高设为设计变量，{W_i，H_i}，i＝1,,…，n；即设计变量总数为2n。性能约束函数控制梁各段的最大弯曲应力，即σ_maxi≦400MPa,i＝1,,…，n。重量函数

Δl_i为阶梯梁各段的长度。具体过程包括以下步骤：

步骤1：在工程软件Nastran中建立直梁有限元数值模型，按优化设计变量的分段数设置梁单元数量、材料属性组以及几何属性参数组{a_i，b_i}i＝1,,…，n；本实施例的梁分段截面宽高初参数均取50mm，设计变量的上下限取为[10，50]。由于本实施例为悬臂梁结构，载荷形式较为简单，可按宽度的下限值计算理论解，给出不同分段数的阶梯梁在最大弯曲应力约束控制下的理论重量结果，详见表2中最后一列数据，以用于校验和对比不同优化算法的有效性及其稳健性。

步骤2：运用Nastran软件建立在应力约束控制下的重量优化模型，并分别调用该软件推荐用户使用的MFD和SQP数值优化方法进行数值优化迭代计算。提取优化计算结果，表2第二和第三行块分别列入了Nastran软件输出的不同分段阶梯梁的数值优化结果。

步骤3：编制本发明数值优化算法软件程序。本实施例可通过两种途径实现与本发明优化算法程序的数据传递：其一，编制本实施例函数及梯度计算的子例程，这包括：阶梯梁的重量函数及其关于设计变量的一阶导数：

和阶梯梁各段最大弯曲应力及其关于设计变量的一阶导数：

式中，x_i为阶梯梁各段左端的长度坐标值。

其二，编制调用Nastran软件外部执行的驱动程序以及读取Nastran软件输出数据的接口子例程，并完成与本发明优化算法程序链接。优化算法在迭代过程中直接调用Nastran软件的重量函数、性能约束函数及其一阶导函数的数值计算数据。

步骤4：运行本发明数值优化算法迭代程序，所得阶梯梁不同分段数的数值优化结果见表2第一行块。对比结果可以得出以下基本结论：

1.本发明的非线性数值优化算法具有充分的计算稳健性以及高效率，与理论解误差一般不大于1％，在200个变量范围内优化计算时间与变量数呈粗略线性关系；

2.Nastran软件的MFD方法不仅较本发明算法的计算效率低，而随设计变量数的增多，优化的重量结果严重偏离理论结果，表明严重背离了应力约束控制；

3.Nastran软件的SQP方法较本发明算法的计算效率偏低，不同变量数下的优化计算结果虽有效但误差较大；当设计变量数增至1000时，该算法中断，可能与内存需求量大等原因有关。

表2本发明算法与Nastran软件算法的结果对比一览表

注：表中的横杠表示软件不能输出计算结果。

Claims (4)

1.一种工程结构超高维大规模多约束非线性优化方法，其特征在于，包括如下步骤：

1)构建结构力学数值分析的有限元模型，置优化迭代计数k＝0，设置设计变量初值x₀及优化计算收敛精度ε，将设计变量初值x₀赋给需优化设计的结构件，其中，设计变量包括：工程结构单元数量以及几何属性参数组；

2)计算分析结构的性能约束函数、结构的重量目标函数以及所述性能约束函数和结构的重量目标函数各自关于设计变量的一阶梯度向量；

3)利用拉格朗日乘子法的增广形式,将所述结构的重量目标函数与性能约束函数组合成无约束化的子问题目标函数Φ_k；

4)对子问题的无约束化目标函数在给定设计变量点处进行线化处理，得到子问题无约束化目标函数的线化方程组系数矩阵，选择子空间的维数m，采用Krylov子空间迭代方法计算该线化方程组系数矩阵在子空间中的最优解；

5)根据所述线化方程组系数矩阵的对称性，对无约束化目标函数的线化方程组系数矩阵实施非Gram-Schmidt算法的正交化改造，获得非完全正交化的标准基向量组以及Krylov子空间最优解计算所需的矩阵分解式；

6)对矩阵分解式中的三对角对称矩阵进行特征值谱分解，进而代入子空间最优解方程，计算得子空间标准正交基向量组合系数的近似解基向量组合系数；

7)对所述的三对角对称矩阵的特征值及组合系数的近似解实施稳健性校正处理，即逐一检查其所有特征值，当特征值小于零时，将相应的组合系数解反号；当特征值等于零时，将相应的组合系数为零，并计算线性搜索方向及线性搜索设计变量值的充分解以及残量；

8)判断所述残量是否小于或等于预设的迭代误差精度值：是，则转至步骤9)；否，则将迭代过程计算得到的一组设计变量值x_k替代所述设计变量初始值，优化迭代计数k＝k+1,并转至步骤2)开始新一轮优化迭代计算；

9)判断所得设计变量值是否达到优化收敛精度ε：是，则结束优化；否，则迭代计数k加1，将迭代过程计算得到的一组设计变量值替代前次设计变量值，修改需优化的结构件参数，转至步骤2)继续迭代优化。

2.根据权利要求1所述的一种工程结构超高维大规模多约束非线性优化方法，其特征在于，所述步骤5)中的非Gram-Schmidt算法的正交化改造，包括如下步骤：

2-1)取q₁＝r⁰/||r⁰||₂，计算α₁＝(q₁，Aq₁)，r¹＝Aq₁-α₁q₁，取β₁＝||r¹||₂，q₂＝r¹/||r¹||₂；

2-2)取j≥2，计算α_j＝(q_j，Aq_j)，r^j＝Aq_j-α_jq_j-β_j-1q_j-1，β_j＝||r^j||₂，q_j+1＝r^j/||r^j||₂；

2-3)当j＝m+1时，结束计算；否则，j＝j+1，返回步骤2-2)；

上式中，r⁰为第k次迭代求解中计算线化方程组的初始残余向量，||·||₂表示模运算，q₁为单位残余向量，A为线化方程组系数构成的矩阵，Aq₁表示矩阵与向量的乘积运算，(·，·)表示两向量的点积运算，α₁即为向量点积运算数值，β₁为向量的模；m为选定的基向量数，上标j表示第j次循环计算，r¹，...，r^j，...，r^m为所提算法循环计算得到的第1至第m个残余向量；

上述3个计算步骤即得非完全正交化的标准基向量组

并由理论可知所述计算步骤可写作以下矩阵分解式：

其中，U_n×m＝(q₁，q₂，…，q_m)，

称为m阶三对角阵，β_m＝||r_m||₂q_m+1＝r^m/||r^m||₂，

U_n×(m+1)＝(q₁，q₂，…，q_m，q_m+1)。

3.根据权利要求1所述的一种工程结构超高维大规模多约束非线性优化方法，其特征在于，所述步骤6)中的三对角对称矩阵特征值谱分解以及最优近似解基向量组合系数，包括如下步骤：

3-1)利用三对角矩阵的对称性，对三对角矩阵T_m×m进行特征值谱分解为：

T_m×m＝Q·S·Q^T，S＝Diag(λ₁，…，λ_m)

其中，Q为单位上三角矩阵，S为T_m×m的特征值对角矩阵，λ_j，j＝1，…，m为T_m×m阵的特征值，上标T表示矩阵的转置，diag表示对角阵；

3-2)将T_m×m代入前述矩阵分解式，对子空间组合系数最优解予以近似表达，得到如下解耦方程：

式中，y′_m＝Q^Ty_m，β＝||r⁰||₂，

为

去掉最后一个元素，

3-3)计算得子空间组合系数最优解的近似值：

4.根据权利要求1所述的一种工程结构超高维大规模多约束非线性优化方法，其特征在于，所述步骤7)中的稳健性校正算法处理及充分解，包括如下步骤：

4-1)当特征值λ_j＜0，将相应的组合系数解反号，即修改y′_j＝-y′_j；

4-2)当特征值λ_j＝0，相应的组合系数y′_j＝0；

4-3)计算线性搜索方向向量s_k＝U_n×mQ_m×my′_m，并置线性搜索步长α＝1；

4-4)置设计变量解x_k+1＝x_k+αs_k，x_k为本轮子问题设计变量的初始值；

4-5)计算子问题本次搜索的函数值Φ(x_k+1)，若Φ(x_k+1)＜Φ(x_k)则结束；否则，置α＝tα，t＝0.85，返回前步，直至满足，并得充分解x_k。

Images