为适应未来高威胁作战环境的迫切需要，对飞行器的航程提出了越来越高的要求^[1]。由于发射平台的约束，飞行器的质量与尺寸往往有着严苛的限制，要实现高性能的战术技术指标，除了先进总体布局及新理论、新材料、新技术的应用外，必须有高性能的推进系统与之相匹配^[2]。

固体火箭发动机(Solid rocket motor, SRM)具有结构简单、快速反应能力强、工作可靠性高等优点，广泛应用于运载火箭、导弹等飞行器系统中^[3]。传统上在飞行器设计初期，发动机设计单位按给定质量与尺寸要求给出初步的性能参数，供总体设计评估，在不满足航程要求时再迭代设计，由于传统的设计方法没有直接引入航程能力等总体性能的影响^[4]，故难以得到系统最优需求的发动机设计方案。因此在飞行器设计初期，除了不断提升SRM自身性能外，基于总体性能驱动的一体化设计优化是提高飞行器航程的有效途径之一。后者需要结合飞行器总体布局合理设计发动机的内弹道特性，并通过飞行器轨迹合理分配与使用发动机的能量。

早期的一体化设计优化方法受限于计算资源和时间成本，往往采用大量近似解析公式、经验数据，适用范围受限，对于新型高性能飞行器，其设计结果往往与后续详细设计的性能差别较大，容易导致设计的反复^[4-5]。为保证设计优化结果的可信度，需要在分析模型中采用大量的高精度数值计算，这导致进行一次方案评估的计算时间成本很高。若直接采用传统参数优化方法调用精确分析模型进行设计优化，则优化总时间难以承受。因此，有必要采用近似优化(Approximate optimization, AO)方法^[6-11]进行求解，即采用计算代价较小的代理模型来代替复杂且耗时的数值分析，从而以较小的时间成本获得较优的设计方案。

AO方法的策略主要分为：直接法(One-shot approach, OA)^[6-8]和序列迭代法(Updating approach, UA)^[7-11]。OA策略仅构建一次代理模型，并通过一定准则校验代理模型的全局精度，而后在该代理模型上进行优化。UA策略则采用较少的初始样本点构建低精度代理模型，并根据一定加点准则，在序列优化过程中不断更新代理模型的信息，直至满足优化停止条件。相比于OA策略，UA策略无需校验代理模型的全局精度，特别是对于高维优化问题，UA策略调用耗时模型的次数更少。近年来，基于UA策略的序列近似优化(Sequential approximate optimization, SAO)方法得到了广泛发展^[7-11]。例如，基于Kriging代理模型的高效全局优化(Efficient global optimization, EGO)算法^[8]，其定义了目标函数值改善期望(Expected improvement, EI)，并选择EI最大值点作为新的样本点，该算法具有较好的全局搜索性，但在最优解附近收敛速度较慢^[9]；基于目标函数预测值最小(Minimizing the predicted objective function, MP)准则的优化算法，其属于贪婪算法，收敛速度较快，但容易陷入局部最优解；结合多种加点准则的并行加点算法^{[7, 9, 11]}，其能够弥补单一加点准则的不足，目前仍是研究的热点。

本文以固体动力助推滑翔式飞行器为背景，针对其规模约束严格、飞行环境严酷及高性能需求的特点，提出面向航程能力的SRM方案设计优化方法，建立复杂三维装药发动机性能分析的参数化模型，以及基于自适应伪谱法的多约束飞行器航程能力评估模型，提出基于Kriging代理模型的改进混合整数序列优化算法，开展发动机与飞行器轨迹一体化设计优化。

1 模型建立 1.1 参数优化问题

不失一般性，参数优化问题可表示为

$ \begin{aligned} & \min f(z) \\ & \text { s.t. }\left\{\begin{array}{l} g_i(z) \leqslant 0\left(i=1, 2, \cdots, n_g\right) \\ h_j(z)=0\left(j=1, 2, \cdots, n_h\right) \end{array}\right. \end{aligned} $

(1)

式中：z为设计变量，f(z)为目标函数，g_i(z)为第i个不等式约束，n_g为不等式约束数量，h_j(z)为第j个等式约束，n_h为等式约束数量。

对于SRM总体方案设计优化问题，本文以单级助推滑翔式飞行器为例，在不改变飞行器气动外形的条件下，通过设计优化发动机的装药参数和喷管参数，即设计变量z，来提升飞行器的航程能力，即目标函数f(z)，其中要求发动机总质量不大于要求值，即不等式约束g_i(z)，发动机的外径、总长度和装药质量为固定值，即等式约束h_j(z)。

$ \left\{\begin{array}{l} f(z)=\left.s_{\max }\right|_z \\ h_1(z)=\left.d_{\mathrm{M}}\right|_z=d_{\mathrm{Md}} \\ h_2(z)=\left.l_{\mathrm{M}}\right|_z=l_{\mathrm{Md}} \\ h_3(z)=\left.m_{\mathrm{p}}\right|_z=m_{\mathrm{pd}} \\ g_1(z)=\left.m_{\mathrm{M}}\right|_z \leqslant m_{\mathrm{M} \max } \end{array}\right. $

(2)

式中：s为飞行器的航程；d_M、d_Md分别为发动机外径和要求值；l_M、l_Md分别为发动机总长度和要求值；m_p、m_pd分别为装药质量和要求值；m_M为发动机总质量；下标max表示最大值或上限值；|_z表示设计方案z所对应的结果。

值得注意的是，设计变量z中除了实数型变量外，往往还包括整数型变量，例如翼柱型装药的翼槽数量等。因此，式(1)所示的参数优化问题为混合整数优化问题。

1.2 发动机性能分析模型 1.2.1 装药模型

考虑到翼柱形装药结构完整性好，体积装填分数高，燃面可调范围大，目前广泛应用于各类SRM中^[12-13]。本文发动机装药采用前、后三角翼柱形构型，如图 1所示。

1.2.2 壳体模型

发动机前、后封头均采用椭球比相同的椭球封头，壳体圆筒段的厚度δ_c1采用最大应变能强度理论确定，椭球封头厚度沿椭球表面从边缘向中心逐渐加厚，顶点处的厚度δ_c2采用最大应力强度理论确定，其计算式^[13]为

$ \left\{\begin{array}{l} \delta_{\mathrm{c} 1}=\frac{k_{\mathrm{c}} \varphi_{\mathrm{c}} p_{\mathrm{cmax}} d_1}{2.3 \xi_{\mathrm{c}}\left[\sigma_{\mathrm{c}}\right]-k_{\mathrm{c}} \varphi_{\mathrm{c}} p_{\mathrm{cmax}}} \\ \delta_{\mathrm{c} 2}=\frac{k_{\mathrm{c}} \varphi_{\mathrm{c}} p_{\mathrm{cmax}} d_1 k_{\mathrm{m}}}{4 \xi_{\mathrm{c}}\left[\sigma_{\mathrm{c}}\right]-k_{\mathrm{c}} \varphi_{\mathrm{c}} p_{\mathrm{cmax}} k_{\mathrm{m}}} \end{array}\right. $

(3)

式中：p_c为燃烧室压强；k_c、φ_c、ξ_c分别为安全系数、压力波动系数和焊缝强度系数，本文分别取为1.25、1.10和0.90；[σ_c]为壳体材料许用应力；d₁为圆筒的内径；k_m为前、后封头的椭球比，取为2。

1.2.3 喷管模型

喷管采用特型喷管，入口段采用椭圆形型面，扩张段采用优化3次多项式型面，入口段与扩张段之间采用圆弧过渡。其中扩张段型面的表达式为^[13]

$ \left\{\begin{array}{l} y_{\mathrm{n}}=k_{\mathrm{n} 0}+k_{\mathrm{n} 1}\left(x_{\mathrm{n}}-l_{\mathrm{nb}}-l_{\mathrm{nt}}\right)+k_{\mathrm{n} 2}\left(x_{\mathrm{n}}-l_{\mathrm{nb}}-l_{\mathrm{nt}}\right)^3 \\ k_{\mathrm{n} 0}=r_{\mathrm{nt}}+r_{\mathrm{nm}}\left(1-\cos \alpha_{\mathrm{nm}}\right) \\ k_{\mathrm{n} 1}=\tan \alpha_{\mathrm{nm}} \\ k_{\mathrm{n} 2}=\left(\tan \alpha_{\mathrm{ne}}-\tan \alpha_{\mathrm{nm}}\right) /\left(3 l_{\mathrm{ne}}^2\right) \\ l_{\mathrm{nb}}=r_{\mathrm{nt}} \sin \frac{\pi}{4} \\ l_{\mathrm{nt}}=r_{\mathrm{nm}} \sin \alpha_{\mathrm{nm}} \\ l_{\mathrm{ne}}=3\left(r_{\mathrm{ne}}-k_{\mathrm{n} 0}\right) /\left(\tan \alpha_{\mathrm{ne}}+2 \tan \alpha_{\mathrm{nm}}\right) \end{array}\right. $

(4)

式中：x_n为扩张段型面上的点到喷管入口段前端面的距离；y_n为扩张段型面上的点到旋成轴线的距离；l_nb、l_nt、l_ne分别为喷管收敛段、过渡段和扩张段的长度；r_nt、r_ne分别为喉部半径和扩张段出口半径；r_nm为过渡段圆弧半径，本文取r_nm = r_nt；α_nm、α_ne分别为扩张段入口和出口的扩张半角，本文分别取值为26°和22°。

1.2.4 质量模型

发动机总质量m_M可以表示为

$ \left\{\begin{array}{l} m_{\mathrm{M}}=m_{\mathrm{p}}+m_{\mathrm{c}}+m_{\mathrm{n}}+m_{\mathrm{a}} \\ m_{\mathrm{p}}=\rho_{\mathrm{p}} V_{\mathrm{p}} \\ m_{\mathrm{c}}=\rho_{\mathrm{c}} V_{\mathrm{c}} \\ m_{\mathrm{n}}=\rho_{\mathrm{nb}} V_{\mathrm{nb}}+\rho_{\mathrm{ne}} V_{\mathrm{ne}} \end{array}\right. $

(5)

式中：ρ_p、V_p分别为装药的密度和体积；ρ_c、V_c分别为燃烧室壳体的材料密度和体积；ρ_q、V_q分别为燃烧室绝热层的材料密度和体积；ρ_nb、V_nb分别为喉衬的材料密度和体积；ρ_ne、V_ne分别为扩张段绝热层的材料密度和体积；m_a为其他装置质量，包括：点火装置、发火装置、连接装置等，按经验估算。其中V_nb和V_ne采用文献[14]方法计算，其余部件体积由几何模型计算获得。不同时刻发动机相对飞行器头部顶点处的质心位置x_GM亦由几何模型计算获得。

1.2.5 推力模型

对于长径比较小的SRM，采用零维内弹道模型进行性能预示可满足工程设计的精度要求^[12]。零维内弹道模型假设燃烧室内的燃速、压强、温度等参数随空间分布均匀，根据质量守恒关系，则可建立燃烧室内的流动控制方程为

$ \left\{\begin{array}{l} \dot{p}_{\mathrm{c}}=\frac{R_{\mathrm{g}} T_{\mathrm{g}}}{V_{\mathrm{ce}}}\left(A_{\mathrm{p}} a_0 p_{\mathrm{c}}^{n_0}\left(\rho_{\mathrm{p}}-\frac{p_{\mathrm{c}}}{R_{\mathrm{g}} T_{\mathrm{g}}}\right)-\frac{\pi r_{\mathrm{nt}}^2 p_{\mathrm{c}}}{C^*}\right) \\ \dot{V}_{\mathrm{ce}}=A_{\mathrm{p}} a_0 p_{\mathrm{c}}^{n_0} \\ \dot{e}_{\mathrm{p}}=a_0 p_{\mathrm{c}}^{n_0} \\ \dot{r}_{\mathrm{nt}}=a_{\mathrm{nt}} \\ \dot{m}_{\mathrm{p}}=-\rho_{\mathrm{p}} A_{\mathrm{p}} a_0 p_{\mathrm{c}}^{n_0} \end{array}\right. $

(6)

式中：V_ce为燃烧室空腔自由容积；R_g、T_g分别为燃气的气体常数和燃烧温度；e_p为燃烧厚度；A_p为装药在燃烧厚度为e_p时的燃面面积；a₀、n₀分别为装药在标准设计状态下的燃速系数和燃速压强指数；C^*为装药的特征速度；a_nt为喉部烧蚀速率，取为0.24 mm/s。

则发动机推力大小F_T为

$ \left\{\begin{array}{l} F_{\mathrm{T}}=F_{\mathrm{T} 0}-\pi r_{\mathrm{ne}}^2 p_{\mathrm{a}} \\ F_{\mathrm{T} 0}=\pi r_{\mathrm{nt}}^2 C_{\mathrm{F} 0} p_{\mathrm{c}} \\ C_{\mathrm{F} 0}=\mathit{\Gamma} \sqrt{\frac{2 \gamma_{\mathrm{g}}}{\gamma_{\mathrm{g}}-1}\left[1-\left(\frac{p_{\mathrm{e}}}{p_{\mathrm{c}}}\right)^{\frac{\gamma_{\mathrm{g}}-1}{\gamma_{\mathrm{g}}}}\right]}+\varepsilon_{\mathrm{e}} \frac{p_{\mathrm{e}}}{p_{\mathrm{c}}} \\ \mathit{\Gamma}=\sqrt{\gamma_{\mathrm{g}}}\left(\frac{2}{\gamma_{\mathrm{g}}+1}\right)^{\frac{\gamma_{\mathrm{g}}+1}{2\left(\mathrm{~g}_{\mathrm{g}}-1\right)}} \end{array}\right. $

(7)

式中：F_T0为动推力，γ_g为燃气比热比，p_e为喷管出口压力，p_a为环境压力，ε_e为喷管膨胀比，其中ε_e与p_e存在关系为

$ {\varepsilon _{\rm{e}}} = \frac{{r_{{\rm{ne}}}^2}}{{r_{{\rm{nt}}}^2}} = {\left( {\frac{2}{{{\gamma _{\rm{g}}} + 1}}} \right)^{\frac{1}{{{\gamma _{\rm{g}}} - 1}}}}\sqrt {\frac{{{\gamma _{\rm{g}}} - 1}}{{{\gamma _{\rm{g}}} + 1}}} /\sqrt {\frac{{{p_{\rm{e}}}^{\frac{2}{{{\gamma _{\rm{g}}}}}}}}{{{p_{\rm{c}}}}} - {{\frac{{{p_{\rm{e}}}}}{{{p_{\rm{c}}}}}}^{\frac{{{\gamma _{\rm{g}}} + 1}}{{{\gamma _{\rm{g}}}}}}}} $

(8)

可采用牛顿迭代法求解式(8)以获得p_e。

对于复杂三维装药的燃面计算问题，最小距离函数(Minimum distance function, MDF)法能够自动处理燃烧过程中装药几何特征消失/增加等拓扑变化，具有较好的稳定性、通用性和计算精度^{[13, 15]}。该方法首先建立装药体网格到初始燃面面网格的最小距离函数，根据平行层燃烧定律，距离值等于燃烧厚度e_p的等值面即为当前燃面A_p，从而为式(6)提供输入数据。详细的计算步骤见文献[15]，本文不再赘述。

1.3 航程能力评估模型 1.3.1 飞行器运动模型

为有效进行航程能力分析，在方案设计优化阶段可忽略地球扁率和自转影响，建立纵向平面有动力无侧滑的质点运动方程。同时，在方程中引入迎角变化率，以便考虑控制系统限制的影响。则运动方程可表示为

$ \left\{\begin{array}{l} \dot{v}=\frac{F_{\mathrm{T}} \cos \alpha-F_{\mathrm{D}}}{\left(m_{\mathrm{T}}+m_{\mathrm{M}}\right)}+\frac{\mu}{\left(R_{\mathrm{E}}+h\right)^2} \sin \theta \\ \dot{\theta}=\frac{F_{\mathrm{T}} \sin \alpha+F_{\mathrm{L}}}{\left(m_{\mathrm{T}}+m_{\mathrm{M}}\right) v}+\left(\frac{\mu}{v\left(R_{\mathrm{E}}+h\right)^2}+\frac{v}{R_{\mathrm{E}}+h}\right) \cos \theta \\ \dot{h}=v \sin \theta \\ \dot{s}=\frac{R_{\mathrm{E}} v \cos \theta}{R_{\mathrm{E}}+h} \\ \dot{\alpha}=u_{\mathrm{d}} \end{array}\right. $

(9)

式中：v、θ、h、s、α分别为飞行器的速度、航迹角、高度、航程和迎角；u_d为迎角变化率控制指令；m_T为飞行器有效载荷质量，其值为500 kg；μ、R_E分别为地球引力常数和半径，其值分别为3.986 005×10¹⁴ m³/s²和6 371 004 m；F_D、F_L分别为飞行器的阻力和升力，其可表示为

$ \left\{\begin{array}{l} F_{\mathrm{D}}=C_{\mathrm{D}} q S_{\mathrm{R}} \\ F_{\mathrm{L}}=C_{\mathrm{L}} q S_{\mathrm{R}} \\ C_{\mathrm{D}}=C_{\mathrm{A}}(M a, \alpha, \delta) \cos \alpha+C_{\mathrm{N}}(M a, \alpha, \delta) \sin \alpha \\ C_{\mathrm{L}}=C_{\mathrm{N}}(M a, \alpha, \delta) \cos \alpha-C_{\mathrm{A}}(M a, \alpha, \delta) \sin \alpha \\ M a=v / a(h) \\ q=0.5 \rho(h) v^2 \end{array}\right. $

(10)

式中：C_A、C_N分别为轴向力系数和法向力系数；C_D、C_L分别为阻力系数和升力系数；S_R为飞行器的参考面积；Ma为马赫数；q为动压；ρ、a分别为当地大气密度和当地声速；δ为俯仰舵偏角，根据“瞬时平衡假设”有

$ C_{\mathrm{MZ}}(M a, \alpha, \delta)+C_{\mathrm{N}}(M a, \alpha, \delta) \frac{m_{\mathrm{T}} x_{\mathrm{GT}}+m_{\mathrm{M}} x_{\mathrm{GM}}}{L_{\mathrm{R}}\left(m_{\mathrm{T}}+m_{\mathrm{M}}\right)}=0 $

(11)

式中：C_MZ、x_GT分别为相对飞行器头部顶点处的俯仰力矩系数和有效载荷质心位置，x_GT = 1.8 m；L_R为飞行器的参考长度，值为6 m。由于式(10)难以解析求解，故本文采用牛顿迭代法来获得当前Ma和α所对应的δ。

1.3.2 飞行约束条件

助推滑翔式飞行器在大气层内高速飞行，其力热环境严酷，因此，在航程评估时需要考虑防热、载荷、控制的限制约束，即

$ \left\{\begin{array}{l} 0 \leqslant \frac{2.37 \times 10^{-7}}{\sqrt{R_{\mathrm{N}}}}\left(\frac{\gamma-1}{\gamma}\right)^{0.25}\left(\frac{\gamma_{\mathrm{h}}+1}{\gamma_{\mathrm{h}}-1}\right)^{0.25} \sqrt{\rho} v^3 \leqslant V_{Q_{\max }} \\ 0 \leqslant \frac{\left|\left(F_{\mathrm{T}}+F_{\mathrm{D}}\right) \sin \alpha+F_{\mathrm{L}} \cos \alpha\right|}{\left(m_{\mathrm{T}}+m_{\mathrm{M}}\right) g_0} \leqslant n_{y \max } \\ q_{\min }<q \leqslant q_{\max } \\ \delta_{\min }<\delta \leqslant \delta_{\max } \\ \alpha_{\min }<\alpha \leqslant \alpha_{\max } \\ u_{\mathrm{d} \min }<u_{\mathrm{d}} \leqslant u_{\mathrm{d} \max } \end{array}\right. $

(12)

式中：R_N为驻点曲率半径，值为0.03；γ为完全气体比热比，其值为1.4；γ_h为高温气体比热比，其值为1.2；V_Q为驻点热流密度；n_y为法向过载；g₀为海平面引力加速度，其值为9.806 6 m/s²；下标min为最小值或下限值。

此外，飞行末端还需要进行能量管理，即

$ \left\{\begin{array}{l} v_{\mathrm{fmin}} \leqslant v\left(t_{\mathrm{f}}\right) \leqslant v_{\mathrm{f} \max } \\ \theta_{\mathrm{fmin}} \leqslant \theta\left(t_{\mathrm{f}}\right) \leqslant \theta_{\mathrm{fmax}} \\ \alpha_{\mathrm{f} \min } \leqslant \alpha\left(t_{\mathrm{f}}\right) \leqslant \alpha_{\mathrm{f} \max } \end{array}\right. $

(13)

式中：t_f、v_f、θ_f分别为末段的时刻、速度和航迹角。

1.3.3 基于最优控制方法的航程能力评估

当给定一组设计方案z时，根据发动机性能分析模型方法可获得所对应的发动机推力F_T、发动机质量m_M、发动机质心位置x_GM随时间t的变化关系。至此，对于式(2)中目标函数的计算，则可转化为对连续时间最优控制问题的求解，即

$ \begin{aligned} & f(z)=\left.s_{\max }\right|_z=-s\left(t_{\mathrm{f}}\right) \\ & \text { s.t. }\left\{\begin{array}{l} \dot{\boldsymbol{x}}(t)=\boldsymbol{f}(\boldsymbol{x}(t), \boldsymbol{u}(t), t) \\ \boldsymbol{C}_{\min } \leqslant \boldsymbol{C}(\boldsymbol{x}(t), \boldsymbol{u}(t), t) \leqslant \boldsymbol{C}_{\max } \\ \boldsymbol{\varphi}_{\min } \leqslant \boldsymbol{\varphi}\left(\boldsymbol{x}\left(t_0\right), t_0, \boldsymbol{x}\left(t_{\mathrm{f}}\right), t_{\mathrm{f}}\right) \leqslant \boldsymbol{\varphi}_{\max } \end{array}\right. \end{aligned} $

(14)

式中：状态变量x=[v θ h s α]^T；控制变量u=u_d；C为路径约束，见式(12)；φ为端点约束，见式(13)。

上述航程能力评估的实质为求解最优控制问题，即为内层优化。为了保证式(1)所示的参数优化问题能够稳健有效进行，在保证式(14)的求解精度前提下，需要对求解方法的快速性和稳定性提出要求。自适应伪谱法将连续时间最优控制问题转化为非线性规划问题，根据一定准则来保证问题转化的精度，这类方法对于初值不敏感，近年来广泛应用于飞行器轨迹优化中^[16-19]。

本文采用改进自适应Legendre-Gauss-Radau (LGR)伪谱法^[19]作为内层优化求解器。文献[19]验证了该方法对于多约束飞行器轨迹优化问题的求解快速性和可靠性，该方法根据当前解的相对误差来调整配点区间数量，在精度不满足的区间内根据Legendre近似理论增加配点或进行分割，在精度满足的区间内根据多项式系数分析减少配点或进行合并^[18-19]，直至计算精度满足要求。

1.4 发动机方案评估模型

综上所述，进行单次发动机总体方案评估的流程如图 3所示。

2 改进混合整数序列近似优化算法 2.1 Kriging代理模型

Kriging代理模型能够提供任意位置处函数的预测值以及该预测值的统计学误差估计。该误差估计信息可用于指导新样本点的选择。

Kriging代理模型利用已计算的n_s个输入变量(样本点)处的输出值进行线性加权，并利用插值来预测未知输入变量(预测点)z处的输出估计值。根据最小均方误差无偏估计原则^[20]，则最优输出估计值$\hat{f}(z)$为

$ \left\{\begin{array}{l} \hat{f}(z)=\boldsymbol{B}(z) \boldsymbol{b}_{\mathrm{s}}^*+\boldsymbol{r}(z)^{\mathrm{T}} \boldsymbol{G}_{\mathrm{s}}^* \\ \boldsymbol{b}_{\mathrm{s}}^*=\left(\boldsymbol{B}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{R}_{\mathrm{s}}^{-1} \boldsymbol{B}_{\mathrm{s}}\right)^{-1} \boldsymbol{B}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{R}_{\mathrm{s}}^{-1} f_{\mathrm{s}} \\ \boldsymbol{G}_{\mathrm{s}}^*=\boldsymbol{R}_{\mathrm{s}}^{-1} \boldsymbol{f}_{\mathrm{s}}-\boldsymbol{R}_{\mathrm{s}}^{-1} \boldsymbol{B}_{\mathrm{s}} \boldsymbol{b}_{\mathrm{s}}^* \\ \boldsymbol{f}_{\mathrm{s}}=\left[f\left(z_1\right), \cdots, f\left(z_{n_{\mathrm{s}}}\right)\right]^{\mathrm{T}} \\ \boldsymbol{B}_{\mathrm{s}}=\left[\boldsymbol{B}\left(z_1\right), \cdots, \boldsymbol{B}\left(z_{n_{\mathrm{s}}}\right)\right]^{\mathrm{T}} \\ \boldsymbol{R}_s=\left[\begin{array}{ccc} R\left(z_1, z_1\right) & \cdots & R\left(z_1, z_{n_{\mathrm{s}}}\right) \\ \vdots & \ddots & \vdots \\ R\left(z_{n_{\mathrm{s}}}, z_1\right) & \cdots & R\left(z_{n_{\mathrm{s}}}, z_{n_{\mathrm{s}}}\right) \end{array}\right] \\ \boldsymbol{r}(z)=\left[R\left(z, z_1\right), \cdots, R\left(z, z_{n_{\mathrm{s}}}\right)\right]^{\mathrm{T}} \\ R\left(z_i, z_j\right)=\prod\limits_{k=1}^{n_z} R^{(k)}\left(z_i^{(k)}-z_j^{(k)}\right) \end{array}\right. $

(15)

式中：B(z)为预测点的回归基函数矩阵；b_s^*为最优回归系数列向量；B_s为样本点回归基函数矩阵；B(z)b_s^*为f(z)的数学期望；r(z)为预测点的相关函数矩阵；G_s^*为最优相关系数列向量；R_s为样本点的相关函数矩阵；r(z)^TG_s^*体现f(z)的高斯静态随机过程；R为相关函数；R^(k)为相关基函数；上标(k)为输入变量第k维变量所对应的参数。

预测点均方根误差的最优估计值$\hat{\sigma}(z)$为

$ \left\{\begin{array}{l} \hat{\sigma}(z)=\sigma_s \sqrt{1-\boldsymbol{r}(z)^{\mathrm{T}} \boldsymbol{R}_{\mathrm{s}}^{-1} \boldsymbol{r}(z)+\boldsymbol{F}(z)^{\mathrm{T}}\left(\boldsymbol{B}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{R}_{\mathrm{s}}^{-1} \boldsymbol{B}_{\mathrm{s}}\right)^{-1} \boldsymbol{F}(z)} \\ \sigma_{\mathrm{s}}=\sqrt{\frac{\left(f_{\mathrm{s}}-\boldsymbol{B}_{\mathrm{s}} b_{\mathrm{s}}^*\right)^{\mathrm{T}} \boldsymbol{R}_{\mathrm{s}}^{-1}\left(f_{\mathrm{s}}-\boldsymbol{B}_{\mathrm{s}} b_{\mathrm{s}}^*\right)}{n_{\mathrm{s}}}} \\ \boldsymbol{F}(z)=\boldsymbol{B}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{R}_{\mathrm{s}}^{-1} \boldsymbol{r}(z)-\boldsymbol{B}(z) \end{array}\right. $

(16)

本文中回归基函数采用一次多项式模型^{[9, 20]}，相关基函数采用3次样条模型^{[9, 20]}，即

$ \left\{ {\begin{array}{*{20}{left}} {{\boldsymbol{B}}(z) = \left[ {\begin{array}{*{20}{l}} 1&{{z^{(1)}}}& \cdots &{{z^{\left( {{n_z}} \right)}}} \end{array}} \right]}\\ {{R^{(k)}}(\Delta ) = \left\{ {\begin{array}{*{20}{l}} {1 - 15{\xi ^{(k)2}} + 30{\xi ^{(k)3}}, }&{\left( {0≤{\xi _k}≤0.2} \right)}\\ {1.25{{\left( {1 - {\xi ^{(k)}}} \right)}^3}, }&{\left( {0.2 < {\xi _k} < 1.0} \right)}\\ {0, }&{\left( {{\xi _k}≥1.0} \right)} \end{array}} \right.}\\ {{\xi ^{(k)}} = {\theta ^{(k)}}|\Delta |} \end{array}} \right. $

(17)

式中：θ^(k)为模型控制参数，可根据极大似然原则^[20]，求解参数优化问题式(18)获得，即

$ \underset{\left\{\theta^{(1)}, \theta^{(2)}, \cdots, \theta^{\left.\left(n_z\right)\right\}}\right.}{\operatorname{argmin}}\left(n_s \ln \left(\sigma_s^2\right)+\ln \left(\operatorname{det} \boldsymbol{R}_s\right)\right) $

(18)

2.2 混合整数LHD初始采样

提高投影均匀性和空间均布性，能够有效丰富设计空间的信息获取^[7]。拉丁超立方设计(Latin hypercube design, LHD)采样自然满足投影均匀性，已广泛应用于各类近似优化算法的初始采样中。LHD采样点可表示为

$ \left\{\begin{array}{c} {\left[\begin{array}{c} z_1 \\ \vdots \\ z_{n_{\mathrm{s}}} \end{array}\right]=\left[\begin{array}{ccc} z_1^{(1)} & \cdots & z_1^{n_z} \\ \vdots & \ddots & \vdots \\ z_{n_{\mathrm{s}}}^{(1)} & \cdots & z_{n_{\mathrm{s}}}^{n_z} \end{array}\right]=\left[\begin{array}{lll} z^{(1)} & \cdots & z^{\left(n_z\right)} \end{array}\right]} \\ z^{(i)}=\left[\begin{array}{c} z_1^{(i)} \\ \vdots \\ z_{n_{\mathrm{s}}}^{(i)} \end{array}\right]=z_{\min }^{(i)}+\frac{\left(z_{\max }^{(i)}-z_{\min }^{(i)}\right)\left(\boldsymbol{\pi}_i+\rho_i\right)}{n_{\mathrm{s}}} \\ \left(i=1, 2, \cdots, n_z\right) \end{array}\right. $

(19)

式中：π_i=[π_i1, π_i2, …, π_ins]^T为由{1, 2, …, n_s}组成的随机序列；ρ_i=[ρ_i1, ρ_i2, …, ρ_ins]^T为n_s个0~1之间随机数组成的列向量；z_max⁽ⁱ⁾、z_min⁽ⁱ⁾分别为设计变量z第i维的搜索上界和下界。

由于LHD采样具有随机性，难以保证空间均布性，故本文将π_i和ρ_i作为优化参数以改善空间均布性，其中优化准则采用Maximin准则^[7]，即

$ \mathop {{\mathop{\rm argmin}\nolimits} }\limits_{\left\{ {{\pi _1},{\pi _2}, \cdots ,{\pi _{{n_z}}},{\rho _1},{\rho _2}, \cdots ,{\rho _{{n_z}}}} \right\}} \left( {\mathop { - \min }\limits_{1≤i,j≤{n_{\rm{s}}},i \ne j} {{\left\| {\frac{{{{\boldsymbol{\pi }}_i} + {{\boldsymbol{\rho }}_i} - {{\boldsymbol{\pi }}_j} - {{\boldsymbol{\rho }}_j}}}{{{n_{\rm{s}}}}}} \right\|}_2}} \right) $

(20)

对于整型变量，可将其视作实数型变量处理，并通过四舍五入取整以确保变量始终为整数^[21]。

2.3 局部增强的EI加点准则

最大目标函数值改善期望(Expected improvement, EI)加点准则同时利用Kriging代理模型提供的输出估计值和均方根误差估计值，具有较好的全局搜索能力^[8]。该准则假设代理模型预测值Y(z)满足正态分布，其均值为$\hat{f}(z)$，标准差为$\hat{\sigma}(z)$。则可定义输出值Y(z)相对于当前样本点集中最优真实输出值f_min的改进函数I(z)为

$ I(z)= \begin{cases}f_{\min }-Y(z), & \left(Y(z)<f_{\min }\right) \\ 0, & \left(Y(z) \geqslant f_{\min }\right)\end{cases} $

(21)

则I(z)的数学期望E[I(z)]为

$ \left\{ {\begin{array}{*{20}{l}} {E[I(z)] = \left\{ {\begin{array}{*{20}{l}} {{E_1}(z) + {E_2}(z), }&{(\hat \sigma (z) > 0)}\\ {0, }&{(\hat \sigma (z) = 0)} \end{array}} \right.}\\ {{E_1}(z) = \left[ {{f_{\min }} - \hat f(z)} \right]\mathit{\Phi} \left[ {\frac{{{f_{\min }} - \hat f(z)}}{{\hat \sigma (z)}}} \right]}\\ {{E_2}(z) = \hat \sigma (z)\varphi \left[ {\frac{{{f_{\min }} - \hat f(z)}}{{\hat \sigma (z)}}} \right]} \end{array}} \right. $

(22)

式中：Φ[ ·]为标准正态累积分布函数，φ[·]为标准正态概率密度函数。

EI加点准则即选择E[I(z)]值(称为EI值)最大的样本点。对于混合整数参数优化问题，可在满足等式约束、不等式约束以及整数约束^[21]的可行集内搜索EI值最大点作为新的采样点，即

$ \begin{aligned} & \min (-E[I(z)]) \\ & \text { s.t. } \begin{cases}\hat{g}_i(z) \leqslant 0, & \left(i=1, \cdots, n_g\right) \\ \hat{h}_j(z)=0, & \left(j=1, \cdots, n_h\right) \\ z^{(k)} \in \mathbb{Z}, & \left(k=1, \cdots, n_{\mathrm{I}}\right)\end{cases} \end{aligned} $

(23)

式中：$\hat{g}_i(z)$为第i个不等式约束函数或其代理模型, $\hat{h}_j(z)$为第j个等式约束函数或其代理模型, n_I为整型设计变量的数量。

由式(22)可知，EI加点准则可能选择$\hat{f}(z)$值较小或$\hat{\sigma}(z)$较大的点，随着加点的进行，最优解附近的EI值可能较小，而样本点稀疏的区域因$\hat{\sigma}(z)$较大导致EI值较大，使得后续加点倾向于样本点更加稀疏的区域，从而导致最优解的收敛速度较慢^[10]。

本文针对上述问题提出一种局部增强的EI加点准则(简称为LEEI加点准则)。其基本思想是：在当前最优解附近限定一个小的局部搜索区域，在该区域内采用EI加点准则，以提升局部搜索性能；同时在全局搜索区域内采用EI加点准则，以保留原算法的全局收敛性；引入比例参数对局部搜索范围进行动态调整，以控制算法收敛。

则引入局部搜索范围的EI加点优化问题(23)可改写为

$ \begin{aligned} & \min (-E[I(z)]) \\ & \text { s.t. }\left\{\begin{array}{l} \hat{g}_i(z) \leqslant 0, \quad\left(i=1, \cdots, n_g\right) \\ \hat{h}_j(z)=0, \quad\left(j=1, \cdots, n_h\right) \\ z^{(k)} \in \mathbb{Z}, \quad\left(k=1, \cdots, n_{\mathrm{I}}\right) \\ \operatorname{round}\left(r_{\min }^{(k)}\right) \leqslant z^{(k)} \leqslant \operatorname{round}\left(r_{\max }^{(k)}\right) \\ r_{\min }^{(l)} \leqslant z^{(l)} \leqslant r_{\max }^{(l)} \quad\left(l=n_{\mathrm{I}}+1, \cdots, n_z\right) \\ r_{\min }^{(n)}=z^{*(n)}-r_D\left(z_{\max }^{(n)}-z_{\min }^{(n)}\right) \\ r_{\max }^{(n)}=z^{*(n)}+r_D\left(z_{\max }^{(n)}-z_{\min }^{(n)}\right) \end{array}\right. \end{aligned} $

(24)

式中：round(·)为四舍五入取整算子，r_D为搜索空间缩放因子，z^*(n)为当前最优解z^*为的第n维值。

2.4 算法流程

基于LEEI加点准则的改进混合整数SAO算法的具体计算步骤为：

1) 确定初始采样点数量(可取n_s=10n_z)，求解参数优化问题式(20)获得初始样本点集{z₁, z₂, …, z_ns}，计算各样本点真实的目标函数值和约束函数值；初始化迭代次数k=0，并设置最大迭代次数k_max(可取k_max=200)和r_D的初始值(可取r_D=0.2)。

2) 若当前样本点集中存在可行解，则获取真实目标函数最小的可行解作为当前最优解z_k^*；否则获取真实约束函数违反程度最小的解作为当前最优解z_k^*，即

$ \mathop {\min }\limits_{z \in \left\{ {{z_1},{z_2}, \cdots ,{z_{\left. {{n_{\rm{s}}}} \right\}}}} \right.} \left( {\sum\limits_{i = 0}^{{n_g}} {\max } \left( {{g_i}(z),0} \right) + \sum\limits_{j = 0}^{{n_h}} {\left| {{h_j}(z)} \right|} } \right) $

(25)

3) 令k=k+1，基于当前样本点集构建目标函数和约束函数的Kriging代理模型，分别求解参数优化问题式(23)和式(24)获得全局EI点z_G和局部EI点z_L。

4) 若z_G在局部搜索区域内(即z_G满足式(24)中的约束3、4)，则将z_L加入样本点集中, 计算z_L真实的目标函数值和约束函数值；否则将z_G和z_L同时加入样本点集中，分别计算z_G和z_L真实的目标函数值和约束函数值。

5) 运行步骤2)，若当前最优解z_k^*与z_L相等，则令r_D=max(0.5r_D, r_Dmin)；否则令r_D=min(2r_D, r_Dmax)；其中r_Dmax、r_Dmin分别为最大和最小的搜索空间缩放因子(可取r_Dmax=0.4，r_Dmin=0.001)。

6) 判断收敛条件式(26)，其中ε_i(i=1, 2, 3, 4)为收敛精度(可根据目标函数和约束函数的值域对应选择ε_i的大小)。若满足则停止；否则返回步骤3)。

$ \left\{\begin{array}{l} \left|\hat{f}\left(z_k^*\right)-f\left(z_k^*\right)\right| \leqslant \varepsilon_1\left(\left|f\left(z_k^*\right)\right|+10^{-3}\right) \\ \left|f\left(z_k^*\right)-f\left(z_{k-1}^*\right)\right| \leqslant \varepsilon_2\left(\left|f\left(z_k^*\right)\right|+10^{-3}\right) \\ g\left(z_k^*\right) \leqslant \varepsilon_3 \\ \left|h\left(z_k^*\right)\right| \leqslant \varepsilon_4 \\ k \geqslant k_{\max } \end{array}\right. $

(26)

3 仿真验证

本文燃烧室壳体材料采用40SiMnCrMoV，绝热层材料采用EPDM，推进剂采用HTPB/AP/Al，喷管喉衬和扩张段材料采用碳/碳复合材料。共优化17个设计变量，包括前、后翼槽数量共2个整型设计变量，以及装药参数、喷管参数、燃速系数共15个实数型设计变量，为保证几何参数化模型的拓扑正确，定义设计变量与几何参数的关系为

$ \left\{\begin{array}{l} n_{\mathrm{w} 1}=z^{(1)} \\ n_{\mathrm{w} 2}=z^{(2)} \\ h_1=z^{(3)} \\ \gamma_1=z^{(4)} \\ \gamma_2=z^{(5)} \\ h_2=z^{(8)} \\ \gamma_3=z^{(9)} \\ \gamma_4=z^{(10)} \\ d_1=z^{(13)} \\ L_1=z^{(14)} \\ r_{\mathrm{nt}}=z^{(15)} \\ \varepsilon_{\mathrm{e}}=z^{(16)} \\ r_1=z^{(6)} \cos \left(\left(\gamma_1+\gamma_2\right) / 2\right) / \cos \left(\left|\gamma_1-\gamma_2\right| / 2\right) \\ r_2=z^{(11)} \cos \left(\left(\gamma_3+\gamma_4\right) / 2\right) / \cos \left(\left|\gamma_3-\gamma_4\right| / 2\right) \\ \delta_1=15+z^{(7)}\left(\delta_{1 \max }-15\right) \\ \delta_2=15+z^{(12)}\left(\delta_{2 \max }-15\right) \\ a_0=z^{(17)} \\ d_{\mathrm{M}}=d_1+2 \delta_{\mathrm{c} 1} \\ l_{\mathrm{M}}=L_1+d_1 / k_{\mathrm{m}}+2 \delta_{\mathrm{c} 2}+l_{\mathrm{nb}}+l_{\mathrm{nt}}+l_{\mathrm{ne}} \end{array}\right. $

(27)

设计变量搜索范围见表 1。设计约束、飞行条件定义见表 2，其中上标P和H分别为发动机工作和不工作时的参数，其余均为全程参数。

本文基于PYTHON脚本语言驱动ABAQUS 6.14实现发动机壳体、绝热层、装药和喷管的几何参数化建模和四面体/三角形网格自动划分；采用文献[15]方法实现燃面计算；采用文献[19]方法实现改进自适应LGR伪谱法；采用经典的DACE工具箱实现Kriging代理模型的构建；基于MATLAB R2019b平台编写混合整数LHD初始采样程序，以及基于LEEI加点准则的改进混合整数SAO算法(简称为LEEI算法)程序，其中混合整数参数优化问题均采用MATLAB自带的混合整数遗传算法(Genetic algorithm, GA)进行求解，LEEI算法中的收敛精度分别设置为：ε₁=ε₂=10^-3，ε₃=0，ε₄=10^-1)；仿真计算环境为Windows 7 2.10 GHz，32.0 GB内存，优化结果如图 4所示。

为说明设计优化的有效性，以初始采样中的最优方案作为优化前的方案(由算法流程步骤2)获得)进行对比，其中优化前的最大航程为554.90 km，装药前、后翼槽数n_w1、n_w2分别为8和6，装药质量m_p为593.82 kg，不满足等式约束要求；采用LEEI算法优化后，最大航程为626.42 km(如图 4(a)所示)，相比优化前提升了12.89%，装药的前、后翼槽数n_w1、n_w2均为8个，发动机外径d_M、总长度l_M和装药质量m_p分别为0.44，2.90，589.97 kg，发动机总质量m_M为737.40 kg，其中等式约束函数违反程度为0.03(≤ε₄=10^-1)，不等式约束函数违反程度为0(≤ε₃=0)，均满足式(2)中的等式和不等式约束要求，验证了LEEI算法对于等式和不等式约束处理的有效性。

由图 4(b)~(e)可知，飞行器的迎角、迎角变化率、驻点热流密度、动压、法向过载和配平舵偏角均全程满足相应的约束要求。其中迎角在初始爬升时达到了约束上界，即最大爬升迎角是航程的限制因素之一；法向过载在发动机工作时达到了约束上界，而驻点热密度全程均未达到约束边界，可知对于本文飞行器，法向过载为主要限制因素，其制约了可用迎角和发动机推力大小；最高点动压达到了约束下限，由于控制系统对动压的限制，导致最大飞行高度受到制约，进而影响了航程；配平舵偏角全程变化平稳，仅在滑翔末段达到了约束下界，即限制了可用迎角。由图 4(d)可知，为提高航程，飞行末端的飞行速度达到了约束下界，而航迹角达到了上界，两者均满足相应约束要求。由此可知，优化结果满足全部飞行约束要求，验证了LEEI算法结果的正确性。

下面为对比分析优化结果，设置2个整型变量为若干固定值，并分别对剩余的15个实数型变量进行优化，优化对比结果如图 5(a)~(d)所示。

由图 5(a)可知，优化方案的最大航程优于其他情况；由图 5(b)~(d)可知，前8后8翼槽方案的推力接近于双推力的形式，即采用短时间大推力辅助飞行器爬升，而后采用长时间小推力持续加速，由于发动机长度、外径以及装药质量被严格限制，因此当翼槽数较少时，则需要翼面和翼厚更大，以满足尺寸质量的约束条件，从而导致发动机的初始推力相对更大，工作时间相对更短；对于翼槽数较多时，结果相反；此外，可以发现优化方案对应的关机点高度和速度也相对更大，从而其最大航程也更大。

为进一步验证LEEI算法的性能，分别采用基于目标函数预测值最小加点准则的SAO算法(简称为MP算法)、基于传统EI加点准则的SAO算法(简称为EI算法)和GA算法对本文问题进行求解与对比，其中MP算法来自文献[14]，EI算法来自MATLAB软件包Surrogate Model Toolbox。为剥离初始采样对SAO算法的性能影响，LEEI、MP、EI算法均采用本文混合整数LHD初始采样算法进行初始采样，采样数量为170个，算法性能对比结果如图 6、7和表 3所示。

由图 6、7可知，随着序列加点与迭代优化，3种SAO算法均得到收敛。MP算法收敛速度最快，满足全部约束要求，但收敛时的目标函数值较大，即优化的最大航程较小。LEEI和EI算法的结果与GA算法结果一致，均满足全部约束要求，且优于MP算法的结果，其中EI算法初始收敛速度较快，但后期收敛速度很慢，耗时函数调用次数约为MP算法的2倍；GA算法虽然优化结果较好，但耗时函数调用次数过大，约为SAO算法的30~60倍；本文LEEI算法的耗时函数调用次数多于MP算法，但相比于EI算法显著减少，由于增加了局部搜索控制策略，从而有效改善了收敛速度，同时保留了原算法的全局搜索性，验证了LEEI算法的有效性。

4 结论

1) LEEI算法能够有效求解SRM方案设计优化的混合整数参数优化问题，优化后飞行器的最大航程提升了12.89%，且设计结果满足全部约束条件。

2) LEEI算法的优化结果与GA算法的优化结果相一致，但耗时函数的调用次数仅为后者的2%。

3) LEEI算法具有一定全局搜索能力，且相比于EI算法有效提高了局部收敛速度。