[发明专利]基于变质心滚控模式的再入飞行器轨迹优化方法

权利要求书

1.一种基于变质心滚控模式的再入飞行器轨迹优化方法，其特征在于采用以下步骤：

步骤一、建立一维变质心滚控再入飞行器数学模型；

根据牛顿第二定理，建立在半速度坐标系的壳体质心平动的动力学矢量方程为

Mt(∂V∂t+ω'×V+2ωe×V+ωe×(ωe×(re+r)))+mz(δ2ρzδt2+2ωT×δρzδt+dωTdt×ρz+ωT×(ωT×ρz))=F---(1)]]>

式中，M_t为包含可移动质量块的弹体的总质量；m_z为可沿弹体系z轴移动的质量块质量；表示壳体质心相对再入系的速度V沿半速度坐标系的投影对时间的一阶导数；ω是半速度坐标系相对于再入系的旋转角速度；ω_e是地球自转角速度；r_e是地心惯性系原点指向再入系原点的矢径；r是再入系原点指向壳体质心的矢径；ω_T是弹体系相对于地心惯性系的旋转角速度；ρ_z表示可移动质量块与壳体质心的相对位移；δρ_zδt表示ρ_z沿弹体系的投影对时间的一阶导数；δ²_ρzδt²表示ρ_z沿再入系的投影对时间的二阶导数；F为弹体受到的合外力矢量；

根据动量矩定理，建立在地心惯性系的绕壳体转动的动力学矢量方程为

(J+ΔJ)δωTδt+ωT×((J+ΔJ)ωT)+Mc+Mgz+Mkz=MR]]>

式中，J为壳体在弹体系下的转动惯量矩阵；△J为可移动质量块所引起的附加转动惯量，ρ_zx、ρ_zz为ρ_z沿弹体系ox₁、oz₁轴的投影；ω_T为弹体系相对于地心惯性系的旋转角速度；M_c为控制力矩；M_gz为质量块的附加相对惯性力所引起的附加相对惯性力矩；M_kz为质量块的附加哥式力所引起的附加哥式力矩；M_R为弹体受到的气动力R对壳体质心O的力矩；

J=Jxx000Jyy000Jzz,ΔJ=mz(1-mzMt)ρzz20-ρzzρzx0ρzz2+ρzx20-ρzzρzx0ρzx2]]>

Mc=mzMtρz×R,Mgz=mz(1-mzMt)ρz×δ2ρzδt2,Mkz=2mz(1-mzMt)ρz×(ωT×δρzδt)]]>

基于下述假设条件：

①“瞬时平衡”假设；

②α＝α_b，β＝0；

③忽略可动质量块位移、速度以及加速度对弹头运动的影响；

则获得简化后的变质心再入飞行器壳体质心运动方程组

MtdVdt=Mtgxh+Rxcosαb-Rysinαb-Mt[ωexhωeyhryh-((ωeyh)2+(ωezh)2)rxh+ωexhωexhrzh]---(3)]]>

MtVcosσdθdt=Mtgyh+Rxcosvsinαb+Rycosvcosαb-Rzsinv-2MtVωezh-Mt[ωeyhωezhrzh-((ωexh)2+(ωezh)2)ryh+ωexhωeyhrxh]---(4)]]>

-MtVdσdt=Mtgzh+Rxsinvsinαb+Rysinvcosαb+Rzcosv+2MtVωeyh-Mt[ωexhωezhrxh-((ωexh)2+(ωeyh)2)rzh+ωeyhωezhryh]---(5)]]>

dxdt=Vcosσcosθ---(6)]]>

dydt=Vcosσsinθ---(7)]]>

dzdt=-Vsinσ---(8)]]>

式中，V为壳体质心在再入系的相对速度大小；[x,y,z]^T为再入系原点指向壳体质心的矢径沿再入系的投影；θ,σ,ν为速度坐标系与再入系之间的速度倾角、航迹偏角和倾侧角；α_b为配平攻角，是关于马赫数和高度的函数；[R_x,R_y,R_z]^T为弹头受到的气动力R沿弹体坐标系的投影；为重力角速度沿半速度坐标系的投影；为地球自转角速度ω_e沿半速度坐标系的投影；为地心指向壳体质心的矢径r_e+r)沿半速度坐标系的投影；

ωexhωeyhωezh=ωecosθcosθ-sinσ-sinθcosθ0sinσcosθsinσsinθcosσsinφ0sinA0cosA0-cosφ0sinA0]]>

rxhryhrzh=cosσcosθcosσsinθ-sinσ-sinθcosθ0sinσcosθsinσsinθcosσxy+Rez]]>

gxhgyhgzh=μ[(rxh)2+(ryh)2+(rzh)2]32rxhryhrzh]]>

其中，φ₀,A₀为再入坐标系与地心坐标系之间的地心纬度和方位角；R_e为地球平均半径R_e＝6371km；μ为地球引力系数，μ＝3.986005×10¹⁴；

步骤二、建立一维变质心滚控再入飞行器轨迹优化模型；

对于给定的受控系统，寻找控制变量U(t)∈R^m，使系统由指定状态转移到期望状态，且使性能指标J取最优值；

再入飞行器的倾侧角ν(t)为控制变量，倾侧角速度作为控制变量；u表示倾侧角的变化速率；

以时间作为离散自变量，对其进行离散，得到离散时间序列T；

离散时间序列T

t₀＝t₁<t₂<…<t_N-1<t_N＝t_f

其中，t₀为初始时刻；t_f为终端时刻(自由或固定)；

以控制变量为设计变量，在离散时间点上，将连续空间的控制变量函数离散为控制序列U；

离散控制序列U

[u₁,u₂,…,u_N-1,u_N]

时间节点之间的控制变量u(t)的值，通过二次样条插值函数S₂(t)获取

u(t)=S2(t)=-12Mi-1(t-ti)2ti-ti-1+12Mi(t-ti-1)2ti-ti-1+12Mi-1(ti-ti-1)+ui-1t∈[ti-1,ti],i={2,3,...,N}---(9)]]>

其中，在给定M₁＝u′(t₁)的值之后，根据下面的迭代公式依次计算M_i；

Mi=2(ui-ui-1)ti-ti-1-Mi-1,i={2,3,...,N}---(10)]]>

选取中间变量的作为控制变量时，变质心再入飞行器的质心运动方程组为

MtdVdt=Mtgxh+Rxcosαb-Rysinαb-Mt[ωexhωeyhryh-((ωeyh)2+(ωezh)2)rxh+rxh+ωexhωexhrzh]MtVcosσdθdt=Mtgyh+Rxcosvsinαb+Rycosvcosαb+cosvcosαb-Rzsinv-2MtVωezh-Mt[ωeyhωezhrzh-((ωexh)2+(ωezh)2)ryh+ωexhωeyhrxh]---(11)Mtgzh+Rxsinvsinαb+Rysinvcosαb+Rzcosv+2MtVωeyh-Mt[ωexhωezhrxh-((ωexh)2+(ωeyh)2)rzh+ωeyhωezhryh]dxdt=Vcosσcosθdydt=Vcosσsinθdzdt=-Vsinσdvdt=u]]>

令状态向量X＝[V,θ,σ,x,y,z,ν]^T，则再入飞行器的质心运动方程组简记为

X.(t)=f(X(t),u(t),t)---(12)]]>

结合再入飞行器的作战任务需求，本发明的性能指标表示为

J=(xf-xm)2(yf-ym)2+(zf-zm)2---(13)]]>

式中，x_f,y_f,z_f为终端时刻弹头的位置；x_m,y_m,z_m为所要攻击的目标位置；

完整的边界条件和过程约束包括：

①飞行时间约束

t_f<t_max           (14)

②终端速度大小、速度倾角约束

V_f>V_min,θ_f<θ_min            (15)

③轴向、横向过载约束

nx=|RxMtg|<nxmaxny=|RyMtg|<nymax---(16)]]>

④一维变质心滚控能力约束

MxJxx-RylzmaxmzMtJxx<du(t)dt<MxJxx+RylzmaxmzMtJxx---(17)]]>

式中，M_x为弹头受到的气动力矩沿弹体纵轴的投影；l_zmax为可移动质量块的最大移动行程；

⑤攻角约束

α_min<α_b<α_max         (18)

⑥控制约束

u_min<u_k<u_max,k＝1,2,…,N       (19)

综合式(9)—式(19)，一维变质心滚控再入飞行器轨迹优化模型为

minJ=(xf-xm)2+(yf-ym)2(zf-zm)2]]>

tf<tmaxVf>Vminθf<θmin|RxMtg|<nxmax|RyMtg|<nymaxMxJxx-RylzmaxmzMtJxx<du(t)dt<MxJxx+RylzmaxmzMtJxxαmin<αb<αmaxumin<uk<umax,k=1,2,...,N---(20)]]>

步骤三、一维变质心滚控再入飞行器轨迹优化模型的求解；

采用模拟退火算法求解一维变质心滚控再入飞行器轨迹优化模型，其基本步骤为：

①给定冷却进度表参数、离散时间序列T及其迭代初始离散控制序列U₀；

②根据初始控制序列U₀，对运动方程组进行数值积分，获取描述运动轨迹的状态变量X₀(t),t∈[t₀,t_f]，计算性能指标函数值J(X₀(t))；

③温度参数TE＝TE_k时作L_k次试探搜索；

根据当前控制序列U_k产生一个随机向量Z_k得到U_k邻域的新的试探点

Uk'=Uk+Zk]]>

对运动方程组进行数值积分，获取描述运动轨迹的状态变量计算性能指标函数值

④产生一个在(0,1)上均匀分布的随机数η，计算出在给定当前迭代点U_k和温度TE_k下与Metropolis接受准则相对应的转移概率P，有：

p=1,J(Xk'(t))≤J(Xk(t))J(Xk(t))-J(Xk'(t))Tke,J(Xk'(t))>J(Xk(t))]]>

若η≤P，则接受新解Uk=Uk',]]>J(Xk(t))=J(Xk'(t)),]]>否则U_k不变；

⑤试探点搜索小于L_k次，转步骤②，否则转步骤⑤；

⑥满足迭代终止条件？

⑦是，则停止，当前解为近似全局最优解；否，则转步骤⑥；

根据给定的温度衰减函数，产生新的温度控制参数TE_k+1及Markov链的长度L_k+1；

⑧重复步骤③-⑦，直到找到最优解U^*。