Closed-Loop Endoatmospheric Ascent Guidance

作者信息：
Ping Lu，爱荷华州立大学，爱荷华州艾姆斯市
Hongsheng Sun
Bruce Tsai

总结人： 鲁鹏，北京理工大学宇航学院
2019.05.09

文章结构

1. 引言
2. 大气层内上升制导问题的数学模型
   1. 上升段动力学方程
   2. 最优控制问题
   3. 确定Φ\PhiΦ和α\alphaα的符号
   4. 添加路径约束
3. 数值方法
   1. 有限差分方法
   2. 更新算法
   3. Continuation on Atmospheric Density和初值猜测
4. 仿真
   1. 终端条件和上升时间调整
   2. 开环解
   3. 闭环仿真
5. 结论

• 附录A：三维最优上升问题的协态微分方程
• 附录B：真空最优上升问题解析解

收获

1. 运载器的在惯性坐标系中的动力学方程中的r\boldsymbol{r}r是运载器的地心位置矢量
2. 质量变化率公式m˙=−ηTvac/(g0Isp)\dot{m} = -\eta T_{vac}/(g_{0}I_{sp})m˙=−ηTvac​/(g0​Isp​)中的g0g_{0}g0​是地球表面的重力加速度。探测器月球表面软着陆也使用这个质量变化公式，g0g_{0}g0​ 依然是地球表面的重力加速度
3. 地球自传给垂直发射运载器带来水平方向（相对于发射惯性坐标系）的初始速度
4. 运载器上升段制导中，常用J=1/rf−Vf2/2J = 1/r_{f}-V^{2}_{f}/2J=1/rf​−Vf2​/2 作为目标函数，使JJJ最小，则末端能量−1/rf+Vf2/2-1/r_{f}+V^{2}_{f}/2−1/rf​+Vf2​/2最大

优缺点

to be added …

难理解的部分

约束S1=qα−Qα≤0S_{1} = q\alpha - Q_{\alpha} \leq 0S1​=qα−Qα​≤0零阶约束，因为控制量1b\boldsymbol{1}_{b}1b​直接出现在约束中（通过α\alphaα），约束S3=q−qmax≤0S_{3} = q - q_{max} \leq 0S3​=q−qmax​≤0是一阶约束，因为控制量1b\boldsymbol{1}_{b}1b​出现在S3S_{3}S3​对时间的一阶导数中。

公式推导

文章中的公式有很多值得注意的细节，文章中的错误也一并记录下来。

1. 真空最优上升解

真空主动段动力学方程

归一化后的方程[1]
r˙=VV˙=−r+T(τ)1b\begin{aligned} \dot{\boldsymbol{r}} =& \boldsymbol{V} \\ \dot{\boldsymbol{V}} = -\boldsymbol{r} +& T(\tau)\boldsymbol{1}_{b} \\ \end{aligned} r˙=V˙=−r+​VT(τ)1b​​

动力学方程可以简化为这种形式的前提是引力加速度做如下近似
g=−(μE/r03)r=−ω2r\boldsymbol{g} = -(\mu_{E}/r^{3}_{0})\boldsymbol{r} = -\omega^{2}\boldsymbol{r} g=−(μE​/r03​)r=−ω2r

为了减小引力简化带来的误差，在每一次制导环的开始更新r0r_{0}r0​.

该方程的解析解如下所示[1-3]
[pv(τ)−pr(τ)]=[cos⁡τI3sin⁡τI3−sin⁡τI3cos⁡τI3][pv0−pr0]=Ω(τ)[pv0−pr0]\begin{bmatrix} \boldsymbol{p}_{v}(\tau)\\ -\boldsymbol{p}_{r}(\tau)\\ \end{bmatrix} = \begin{bmatrix} \cos{\tau}{I}_{3} & \sin{\tau}{I}_{3}\\ -\sin{\tau}{I}_{3} & \cos{\tau}{I}_{3}\\ \end{bmatrix} \begin{bmatrix} \boldsymbol{p}_{v0}\\ -\boldsymbol{p}_{r0}\\ \end{bmatrix} = \Omega(\tau) \begin{bmatrix} \boldsymbol{p}_{v0}\\ -\boldsymbol{p}_{r0}\\ \end{bmatrix} [pv​(τ)−pr​(τ)​]=[cosτI3​−sinτI3​​sinτI3​cosτI3​​][pv0​−pr0​​]=Ω(τ)[pv0​−pr0​​]

[r(τ)V(τ)]=Ω(τ)[r0V0]+Γ(τ)W\begin{bmatrix} \boldsymbol{r}(\tau)\\ \boldsymbol{V}(\tau)\\ \end{bmatrix} = \Omega(\tau) \begin{bmatrix} \boldsymbol{r}_{0}\\ \boldsymbol{V}_{0}\\ \end{bmatrix} + \Gamma(\tau) \mathrm{W} [r(τ)V(τ)​]=Ω(τ)[r0​V0​​]+Γ(τ)W

Γ(τ)=[sin⁡τI3−cos⁡τI3cos⁡τI3sin⁡τI3]\Gamma(\tau) = \begin{bmatrix} \sin{\tau}{I}_{3} & -\cos{\tau}{I}_{3}\\ \cos{\tau}{I}_{3} & \sin{\tau}{I}_{3}\\ \end{bmatrix} Γ(τ)=[sinτI3​cosτI3​​−cosτI3​sinτI3​​]

W=[IcIs]=[∫0τ1pv(ζ)cos⁡(ζ)T(ζ)dζ∫0τ1pv(ζ)sin⁡(ζ)T(ζ)dζ]\mathrm{W} = \begin{bmatrix} I_{c} \\ I_{s} \end{bmatrix} = \begin{bmatrix} \int^{\tau}_{0}\boldsymbol{1}_{pv}(\zeta)\cos(\zeta) T(\zeta) d\zeta \\ \int^{\tau}_{0}\boldsymbol{1}_{pv}(\zeta)\sin(\zeta) T(\zeta) d\zeta \\ \end{bmatrix} W=[Ic​Is​​]=[∫0τ​1pv​(ζ)cos(ζ)T(ζ)dζ∫0τ​1pv​(ζ)sin(ζ)T(ζ)dζ​]
其中，I3I_{3}I3​是3×33\times 33×3的单位矩阵。通过公式可以看出初始协态变量的模(pv0Tpv0+pr0Tpr0)(\boldsymbol{p}_{v0}^{T}\boldsymbol{p}_{v0}+\boldsymbol{p}_{r0}^{T}\boldsymbol{p}_{r0})(pv0T​pv0​+pr0T​pr0​)对协态变量之后的方向变化没有任何影响，只影响协态变量的大小。考虑到状态的解析解只和1pv\boldsymbol{1}_{pv}1pv​有关，所以初始协态变量的模对状态也没有影响。

质量推力和发动机阀门

首先声明，程序中我在求解质量推力和发动机阀门大小时未对质量和时间归一化，只是在每个时间区间用对应的引力加速度对推力加速度T(t)T(t)T(t)进行了归一化。

1. 当T<TmaxT<T_{max}T<Tmax​时

m(t)=m0−TvacctT(t)=Tvacm(t)g0η(t)=1\begin{aligned} m(t) =& m_{0} - \frac{T_{vac}}{c}t \\ T(t) =& \frac{T_{vac}}{m(t)g_{0}} \\ \eta(t) =& 1 \\ \end{aligned} m(t)=T(t)=η(t)=​m0​−cTvac​​tm(t)g0​Tvac​​1​

其中，TvacT_{vac}Tvac​是发动机真空推力，g0g_{0}g0​是运载器所在某时间区间[ti,ti+τ]\lbrack t_{i},t_{i}+\tau \rbrack[ti​,ti​+τ]内的引力加速度值，喷嘴等效气流速度（nozzle exit velocity）c=gEIspc = g_{E}I_{sp}c=gE​Isp​，gEg_{E}gE​是地球表面重力加速度，m(t)m(t)m(t)是运载器质量，m0m_{0}m0​是运载器质量初始值，T(t)T(t)T(t)是被g0g_{0}g0​归一化后的推力加速度，η\etaη是发动机阀门。
随着TTT不断增大，假设在tct_{c}tc​时刻，推力加速度达到最大（取Tmax=4geT_{max}=4g_{e}Tmax​=4ge​），此后推力加速度TTT保持最大，并且tct_{c}tc​时刻在[ti,ti+τ]\lbrack t_{i},t_{i}+\tau\rbrack[ti​,ti​+τ]内，那么
Tvac[m0−(Tvac/c)tc]g0=Tmaxg0\frac{T_{vac}}{\lbrack m_{0}-(T_{vac}/c)t_{c}\rbrack g_{0}} = \frac{T_{max}}{g_{0}} [m0​−(Tvac​/c)tc​]g0​Tvac​​=g0​Tmax​​

所以
tc=c(m0Tvac−1Tmax)t_c = c(\frac{m_{0}}{T_{vac}}-\frac{1}{T_{max}}) tc​=c(Tvac​m0​​−Tmax​1​)

1. 当T=TmaxT=T_{max}T=Tmax​后，即t⩾tct \geqslant t_{c}t⩾tc​时

m(t)=mtce−(Tmax/c)tT(t)=Tmaxg0η(t)=Tmaxm(t)Tvac\begin{aligned} m(t) =& m_{t_{c}}e^{-(T_{max}/c)t} \\ T(t) =& \frac{T_{max}}{g_{0}} \\ \eta(t) =& \frac{T_{max}m(t)}{T_{vac}} \\ \end{aligned} m(t)=T(t)=η(t)=​mtc​​e−(Tmax​/c)tg0​Tmax​​Tvac​Tmax​m(t)​​

此时的g0g_{0}g0​是某时间区间[tj,tj+τ]\lbrack t_{j},t_{j}+\tau \rbrack[tj​,tj​+τ]（注意tj+τ⩾tct_{j}+\tau \geqslant t_{c}tj​+τ⩾tc​）内引力加速度，T(t)T(t)T(t)是被g0g_{0}g0​归一化后的推力加速度，mtcm_{t_{c}}mtc​​是tct_{c}tc​时刻运载器的质量。

垂直上升结束时运载器在发射点惯性坐标系中的位置和速度

运载器在垂直上升结束后制导率才开始工作，不同运载器有不同的垂直上升时间t0t_{0}t0​。文献[5]中垂直上升时间t0t_{0}t0​取9s，文献[1]中垂直上升时间t0t_{0}t0​取5s。
垂直上升阶段运载器推力加速度满足最大推力加速度约束，此时运载器所受加速度
avehicle(t)=Tvac−m(t)gEm(t)a_{vehicle}(t) = \frac{T_{vac} - m(t)g_{E}}{m(t)} avehicle​(t)=m(t)Tvac​−m(t)gE​​

所以，对加速度积分可以得到竖直上升的速度v0x(t)v_{0x}(t)v0x​(t)
v0=[v0x(t0),v0y,v0z]Tv0x(t)=−cln⁡(m0−Tvacct)−gEt+cln⁡(m0)\begin{aligned} v_{0} = &\lbrack v_{0x}(t_{0}), v_{0y}, v_{0z}\rbrack ^{T}\\ v_{0x}(t) = -c\ln(&m_{0}-\frac{T_{vac}}{c}t) - g_{E}t + c\ln(m_{0}) \end{aligned} v0​=v0x​(t)=−cln(​[v0x​(t0​),v0y​,v0z​]Tm0​−cTvac​​t)−gE​t+cln(m0​)​

其中v0y,v0zv_{0y},v_{0z}v0y​,v0z​是由地球自转引起的水平面内的速度在y和z轴上的分量。对速度v0x(t)v_{0x}(t)v0x​(t)积分可得r0x(t)r_{0x}(t)r0x​(t)
r0=[r0x(t0),v0yt0,v0zt0]Tr0x(t)=c2Tvac(m0−Tvacct)[ln⁡(m0−Tvacct)−1]−12gEt2+cln⁡(m0)t+c2Tvacm0[1−ln⁡(m0)]\begin{aligned} r_{0} = \lbrack r_{0x}(t_{0}), v_{0y}&t_{0}, v_{0z}t_{0}\rbrack ^{T} \\ r_{0x}(t) = \frac{c^2}{T_{vac}}(m_{0} - \frac{T_{vac}}{c}t) \lbrack \ln(&m_{0} - \frac{T_{vac}}{c}t)-1 \rbrack - \frac{1}{2}g_{E}t^{2} + \\ c\ln(m_{0})t + \frac{c^2}{T_{vac}}m_{0}\lbrack &1 - \ln(m_{0}) \rbrack \end{aligned} r0​=[r0x​(t0​),v0y​r0x​(t)=Tvac​c2​(m0​−cTvac​​t)[ln(cln(m0​)t+Tvac​c2​m0​[​t0​,v0z​t0​]Tm0​−cTvac​​t)−1]−21​gE​t2+1−ln(m0​)]​

终端约束和横截条件得到的约束方程Ψ=0{\Psi} = \mathbf{0}Ψ=0

方程如下[1]
12rfTrf−12rf∗2=0\frac{1}{2}\boldsymbol{r}^{T}_{f}\boldsymbol{r}_{f} - \frac{1}{2}r^{*2}_{f} = 0 21​rfT​rf​−21​rf∗2​=0

1NT(rf×Vf)−∥rf×Vf∥cos⁡i∗=0\boldsymbol{1}^{T}_{N}(\boldsymbol{r}_{f}\times \boldsymbol{V}_{f})-\lVert \boldsymbol{r}_{f}\times \boldsymbol{V}_{f}\rVert \cos{i^{*}} = 0 1NT​(rf​×Vf​)−∥rf​×Vf​∥cosi∗=0

rfTVf−rfVfsin⁡γf∗=0\boldsymbol{r}^{T}_{f}\boldsymbol{V}_{f} - r_{f}V_{f}\sin{\gamma^{*}_{f}} = 0 rfT​Vf​−rf​Vf​sinγf∗​=0

(VfTprf)rf2−(rfTpVf)Vf2+(rfTVf)(Vf2−rfTprf)=0(\boldsymbol{V}^{T}_{f}\boldsymbol{p}_{rf})r^{2}_{f} - (\boldsymbol{r}^{T}_{f}\boldsymbol{p}_{Vf})V^{2}_{f} + (\boldsymbol{r}^{T}_{f}\boldsymbol{V}_{f})(V^{2}_{f} - \boldsymbol{r}^{T}_{f}\boldsymbol{p}_{rf}) = 0 (VfT​prf​)rf2​−(rfT​pVf​)Vf2​+(rfT​Vf​)(Vf2​−rfT​prf​)=0

VfTpVf−Vf2=0\boldsymbol{V}^{T}_{f}\boldsymbol{p}_{Vf} - V^{2}_{f} = 0 VfT​pVf​−Vf2​=0

(hfTprf)[hfT(rf×1N)]+(hfTpVf)[hfT(Vf×1N)]=0(\boldsymbol{h}^{T}_{f}\boldsymbol{p}_{rf})\lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{r}_{f}\times \boldsymbol{1}_{N})\rbrack + (\boldsymbol{h}^{T}_{f}\boldsymbol{p}_{Vf})\lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{V}_{f}\times \boldsymbol{1}_{N})\rbrack = 0 (hfT​prf​)[hfT​(rf​×1N​)]+(hfT​pVf​)[hfT​(Vf​×1N​)]=0

求Ψ{\Psi}Ψ对状态和协态变量的导数

f1f_1f1​对状态变量和协态变量求导
∂f1∂rf=rf,∂f1∂Vf=0,∂f1∂pVf=0,∂f1∂prf=0\frac{\partial f_{1}}{\partial \boldsymbol{r}_{f}} = \boldsymbol{r}_{f}, \frac{\partial f_{1}}{\partial \boldsymbol{V}_{f}} = \boldsymbol{0}, \frac{\partial f_{1}}{\partial \boldsymbol{p}_{Vf}} = \boldsymbol{0}, \frac{\partial f_{1}}{\partial \boldsymbol{p}_{rf}} = \boldsymbol{0} ∂rf​∂f1​​=rf​,∂Vf​∂f1​​=0,∂pVf​∂f1​​=0,∂prf​∂f1​​=0

f2f_2f2​对状态变量和协态变量求导
∂f2∂rf=Vf×1N+(Vf×)T(rf×Vf)cos⁡i∗∥rf×Vf∥\frac{\partial f_{2}}{\partial \boldsymbol{r}_{f}} = \boldsymbol{V}^{\times}_{f}\boldsymbol{1}_{N} + (\boldsymbol{V}^{\times}_{f})^{T}(\boldsymbol{r}_{f}\times \boldsymbol{V}_{f})\frac{\cos{i^{*}}}{\lVert \boldsymbol{r}_{f}\times \boldsymbol{V}_{f}\rVert} ∂rf​∂f2​​=Vf×​1N​+(Vf×​)T(rf​×Vf​)∥rf​×Vf​∥cosi∗​

∂f2∂Vf=(rf×)T1N+rf×(rf×Vf)cos⁡i∗∥rf×Vf∥\frac{\partial f_{2}}{\partial \boldsymbol{V}_{f}} = (\boldsymbol{r}^{\times}_{f})^{T}\boldsymbol{1}_{N} + \boldsymbol{r}^{\times}_{f}(\boldsymbol{r}_{f}\times \boldsymbol{V}_{f})\frac{\cos{i^{*}}}{\lVert \boldsymbol{r}_{f}\times \boldsymbol{V}_{f}\rVert} ∂Vf​∂f2​​=(rf×​)T1N​+rf×​(rf​×Vf​)∥rf​×Vf​∥cosi∗​

∂f2∂pVf=0,∂f2∂prf=0\frac{\partial f_{2}}{\partial \boldsymbol{p}_{Vf}} = \boldsymbol{0}, \frac{\partial f_{2}}{\partial \boldsymbol{p}_{rf}} = \boldsymbol{0} ∂pVf​∂f2​​=0,∂prf​∂f2​​=0

f3f_3f3​对状态变量和协态变量求导
∂f3∂rf=Vf−(Vfsin⁡γ∗)rf/rf\frac{\partial f_{3}}{\partial \boldsymbol{r}_{f}} = \boldsymbol{V}_{f} - (V_{f}\sin{\gamma^{*})}\boldsymbol{r}_{f}/r_{f} ∂rf​∂f3​​=Vf​−(Vf​sinγ∗)rf​/rf​

∂f3∂Vf=rf−(rfsin⁡γ∗)Vf/Vf\frac{\partial f_{3}}{\partial \boldsymbol{V}_{f}} = \boldsymbol{r}_{f} - (r_{f}\sin{\gamma^{*})}\boldsymbol{V}_{f}/V_{f} ∂Vf​∂f3​​=rf​−(rf​sinγ∗)Vf​/Vf​

∂f3∂pVf=0,∂f3∂prf=0\frac{\partial f_{3}}{\partial \boldsymbol{p}_{Vf}} = \boldsymbol{0}, \frac{\partial f_{3}}{\partial \boldsymbol{p}_{rf}} = \boldsymbol{0} ∂pVf​∂f3​​=0,∂prf​∂f3​​=0

f4f_{4}f4​对状态变量和协态变量求导
∂f4∂rf=2(VfTprf)rf−Vf2pvf+Vf2Vf−(rfTprf)Vf−(rfTVf)prf\frac{\partial f_{4}}{\partial \boldsymbol{r}_{f}} = 2(\boldsymbol{V}^{T}_{f}\boldsymbol{p}_{rf})\boldsymbol{r}_{f} - V^{2}_{f}\boldsymbol{p}_{vf} + V^{2}_{f}\boldsymbol{V}_{f} - (\boldsymbol{r}^{T}_{f}\boldsymbol{p}_{rf})\boldsymbol{V}_{f} - (\boldsymbol{r}^{T}_{f}\boldsymbol{V}_{f})\boldsymbol{p}_{rf} ∂rf​∂f4​​=2(VfT​prf​)rf​−Vf2​pvf​+Vf2​Vf​−(rfT​prf​)Vf​−(rfT​Vf​)prf​

∂f4∂Vf=rf2prf−2(rfTpVf)Vf+Vf2rf+2(rfTVf)Vf−(rfTprf)rf\frac{\partial f_{4}}{\partial \boldsymbol{V}_{f}} = r^{2}_{f}\boldsymbol{p}_{rf} - 2(\boldsymbol{r}^{T}_{f}\boldsymbol{p}_{Vf})\boldsymbol{V}_{f} + V^{2}_{f}\boldsymbol{r}_{f} + 2(\boldsymbol{r}^{T}_{f}\boldsymbol{V}_{f})\boldsymbol{V}_{f} - (\boldsymbol{r}^{T}_{f}\boldsymbol{p}_{rf})\boldsymbol{r}_{f} ∂Vf​∂f4​​=rf2​prf​−2(rfT​pVf​)Vf​+Vf2​rf​+2(rfT​Vf​)Vf​−(rfT​prf​)rf​

∂f4∂pVf=−Vf2rf,∂f4∂prf=rf2Vf−(rfTVf)rf\frac{\partial f_{4}}{\partial \boldsymbol{p}_{Vf}} = -V^{2}_{f}\boldsymbol{r}_{f}, \frac{\partial f_{4}}{\partial \boldsymbol{p}_{rf}} = r^{2}_{f}\boldsymbol{V}_{f} - (\boldsymbol{r}^{T}_{f}\boldsymbol{V}_{f})\boldsymbol{r}_{f} ∂pVf​∂f4​​=−Vf2​rf​,∂prf​∂f4​​=rf2​Vf​−(rfT​Vf​)rf​

f5f_{5}f5​对状态变量和协态变量求导
∂f5∂rf=0,∂f5∂Vf=pVf−2Vf\frac{\partial f_{5}}{\partial \boldsymbol{r}_{f}} = \boldsymbol{0}, \frac{\partial f_{5}}{\partial \boldsymbol{V}_{f}} = \boldsymbol{p}_{Vf}-2\boldsymbol{V}_{f} ∂rf​∂f5​​=0,∂Vf​∂f5​​=pVf​−2Vf​

∂f5∂pVf=Vf,∂f5∂prf=0\frac{\partial f_{5}}{\partial \boldsymbol{p}_{Vf}} = \boldsymbol{V}_{f}, \frac{\partial f_{5}}{\partial \boldsymbol{p}_{rf}} = \boldsymbol{0} ∂pVf​∂f5​​=Vf​,∂prf​∂f5​​=0

f6f_{6}f6​对状态变量和协态变量求导
∂f6∂rf=[hfT(rf×1N)]Vf×prf+(hfTprf)[(Vf×)T1N×+(1N×)TVf×]rf+[hfT(Vf×1N)]Vf×pVf+(hfTpVf)[(Vf×)21N]\begin{aligned} \frac{\partial f_{6}}{\partial \boldsymbol{r}_{f}} = &\lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{r}_{f}\times \boldsymbol{1}_{N})\rbrack \boldsymbol{V}^{\times}_{f}\boldsymbol{p}_{rf} + (\boldsymbol{h}^{T}_{f}\boldsymbol{p}_{rf})\lbrack (\boldsymbol{V}^{\times}_{f})^{T}\boldsymbol{1}^{\times}_{N} + (\boldsymbol{1}^{\times}_{N})^{T}\boldsymbol{V}^{\times}_{f}\rbrack \boldsymbol{r}_{f}\\ &+\lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{V}_{f}\times \boldsymbol{1}_{N})\rbrack \boldsymbol{V}^{\times}_{f}\boldsymbol{p}_{Vf} + (\boldsymbol{h}^{T}_{f}\boldsymbol{p}_{Vf}) \lbrack (\boldsymbol{V}^{\times}_{f})^{2}\boldsymbol{1}_{N}\rbrack \end{aligned} ∂rf​∂f6​​=​[hfT​(rf​×1N​)]Vf×​prf​+(hfT​prf​)[(Vf×​)T1N×​+(1N×​)TVf×​]rf​+[hfT​(Vf​×1N​)]Vf×​pVf​+(hfT​pVf​)[(Vf×​)21N​]​

∂f6∂Vf=[hfT(rf×1N)](rf×)Tprf+(hfTprf)[−(rf×)21N]+[hfT(Vf×1N)](rf×)TpVf+(hfTpVf)(rf×1N×+1N×rf×)Vf\begin{aligned} \frac{\partial f_{6}}{\partial \boldsymbol{V}_{f}} = &\lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{r}_{f}\times \boldsymbol{1}_{N})\rbrack (\boldsymbol{r}^{\times}_{f})^{T}\boldsymbol{p}_{rf} + (\boldsymbol{h}^{T}_{f}\boldsymbol{p}_{rf})\lbrack -(\boldsymbol{r}^{\times}_{f})^{2}\boldsymbol{1}_{N}\rbrack \\ &+\lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{V}_{f}\times \boldsymbol{1}_{N})\rbrack (\boldsymbol{r}^{\times}_{f})^{T}\boldsymbol{p}_{Vf} + (\boldsymbol{h}^{T}_{f}\boldsymbol{p}_{Vf})(\boldsymbol{r}^{\times}_{f}\boldsymbol{1}^{\times}_{N} + \boldsymbol{1}^{\times}_{N}\boldsymbol{r}^{\times}_{f}) \boldsymbol{V}_{f} \end{aligned} ∂Vf​∂f6​​=​[hfT​(rf​×1N​)](rf×​)Tprf​+(hfT​prf​)[−(rf×​)21N​]+[hfT​(Vf​×1N​)](rf×​)TpVf​+(hfT​pVf​)(rf×​1N×​+1N×​rf×​)Vf​​

∂f6∂pVf=[hfT(Vf×1N)]hf\frac{\partial f_{6}}{\partial \boldsymbol{p}_{Vf}} = \lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{V}_{f}\times \boldsymbol{1}_{N})\rbrack \boldsymbol{h}_{f} ∂pVf​∂f6​​=[hfT​(Vf​×1N​)]hf​

∂f6∂prf=[hfT(rf×1N)]hf\frac{\partial f_{6}}{\partial \boldsymbol{p}_{rf}} = \lbrack \boldsymbol{h}^{T}_{f}(\boldsymbol{r}_{f}\times \boldsymbol{1}_{N})\rbrack \boldsymbol{h}_{f} ∂prf​∂f6​​=[hfT​(rf​×1N​)]hf​

计算雅克比矩阵

计算雅克比矩阵，用于牛顿迭代法求解协态变量初值，过程跟参考文献[2]方法相同。考虑到完整性，在下面列出推导过程

定义UT=[STλT]{U}^{T}=\left[S^{T}\quad \lambda^{T}\right]UT=[STλT]，JT=[IcT(τ)IsT(τ)]{J}^{T}=\left[I_{c}^{T}(\tau)\quad I_{s}^{T}(\tau)\right]JT=[IcT​(τ)IsT​(τ)]，其中τ=t−ti\tau=t-t_{i}τ=t−ti​，则有下式成立
dU(ti+τ)=B(τ)dU(ti),dUT(t0)=[0dλ0T]d {U}\left(t_{i}+\tau\right)=B(\tau) d{U}\left(t_{i}\right) ,\quad d{U}^{T}\left(t_{0}\right)=\left[ \begin{array}{ll}{0} & {d \lambda_{0}^{T}}\end{array}\right] dU(ti​+τ)=B(τ)dU(ti​),dUT(t0​)=[0​dλ0T​​]

其中
B(τ)=[ΩΓdJ(τ)dλ(ti)0Ω]=[ΩΓ(τ)[dIc/dλ(ti)dIs/dλ(ti)]0Ω]B(\tau)=\left[ \begin{array}{cc}{\Omega} & {\Gamma \frac{d \boldsymbol{J}(\tau)}{d \boldsymbol{\lambda}\left(t_{i}\right)}} \\ \boldsymbol{0} & {\Omega}\end{array}\right]= \begin{bmatrix} \Omega & \Gamma(\tau)\begin{bmatrix} dI_{c}/{d}\lambda(t_{i})\\ {dI_{s}}/{d}\lambda(t_{i})\end{bmatrix}\\ \boldsymbol{0} & \Omega \end{bmatrix} B(τ)=[Ω0​Γdλ(ti​)dJ(τ)​Ω​]=⎣⎡​Ω0​Γ(τ)[dIc​/dλ(ti​)dIs​/dλ(ti​)​]Ω​⎦⎤​

K(τ)=cos⁡(τ)T(τ)∥pV(τ)∥[I3−pV(τ)(pV(τ))T∥pV(τ)∥2][cos⁡(τ)I3sin⁡(τ)I3]K(\tau)=\frac{\cos(\tau)T(\tau)}{\lVert \boldsymbol{p}_{V}(\tau)\rVert}\lbrack {I}_{3} - \frac{\boldsymbol{p}_{V}(\tau)({\boldsymbol{p}_{V}(\tau)})^{T}}{{\lVert \boldsymbol{p}_{V}(\tau)\rVert}^{2}}\rbrack \lbrack \cos(\tau){I}_{3} \enspace \sin(\tau){I}_{3}\rbrack K(τ)=∥pV​(τ)∥cos(τ)T(τ)​[I3​−∥pV​(τ)∥2pV​(τ)(pV​(τ))T​][cos(τ)I3​sin(τ)I3​]

dIc(τ)dλ(ti)=τ[7K(0)+32K(δ)+12K(2δ)+32K(3δ)+7K(4δ)]90\frac{d{I}_{c}(\tau)}{d{\lambda}(t_{i})} = \frac{\tau \lbrack 7{K}(0) + 32{K}(\delta) + 12{K}(2\delta) + 32{K}(3\delta) + 7{K}(4\delta)\rbrack}{90} dλ(ti​)dIc​(τ)​=90τ[7K(0)+32K(δ)+12K(2δ)+32K(3δ)+7K(4δ)]​

dIs(τ)dλ(ti)=τ[32K(δ)tan⁡(δ)+12K(2δ)tan⁡(2δ)+32K(3δ)tan⁡(3δ)+7K(4δ)tan⁡(4δ)]90\frac{d{I}_{s}(\tau)}{d{\lambda}(t_{i})} = \frac{\tau \lbrack 32{K}(\delta)\tan(\delta) + 12{K}(2\delta)\tan(2\delta) + 32{K}(3\delta)\tan(3\delta) + 7{K}(4\delta)\tan(4\delta)\rbrack}{90} dλ(ti​)dIs​(τ)​=90τ[32K(δ)tan(δ)+12K(2δ)tan(2δ)+32K(3δ)tan(3δ)+7K(4δ)tan(4δ)]​

Θ=∏i=1M−1B(ti+1−ti)\Theta = \prod^{M-1}_{i=1} B(t_{i+1}-t_{i}) Θ=i=1∏M−1​B(ti+1​−ti​)

dΨ=ΨSdSf+Ψλdλf=(ΨSΘ12+ΨλΘ22)dλ0d{\Psi} = \Psi_{S}d{S}_{f} + \Psi_{\lambda}d{\lambda}_{f} = (\Psi_{S}\Theta_{12} + \Psi_{\lambda}\Theta_{22} )d{\lambda}_{0} dΨ=ΨS​dSf​+Ψλ​dλf​=(ΨS​Θ12​+Ψλ​Θ22​)dλ0​

其中Θ12\Theta_{12}Θ12​和Θ22\Theta_{22}Θ22​是矩阵Θ\ThetaΘ的分块矩阵，因此雅克比矩阵Jacob\mathrm{Jacob}Jacob
Jacob=dΨdλ0=ΨSΘ12+ΨλΘ22\mathrm{Jacob} = \frac{d{\Psi}}{d{\lambda}_{0}} = \Psi_{S}\Theta_{12} + \Psi_{\lambda}\Theta_{22} Jacob=dλ0​dΨ​=ΨS​Θ12​+Ψλ​Θ22​

2. 大气层内上升解

作者将大气层内上升问题转化为了两点边值问题(TPBVPs)，并使用有限差分法求解两点边值问题。解两点边值问题之前先要推导状态方程和协态方程中涉及的相关方程。

过程约束的处理

处理不等式约束时，当约束S<0S < 0S<0时，系统的协态变量微分方程取−∂H/∂x-\partial{H}/\partial{\boldsymbol{x}}−∂H/∂x当约束S=0S = 0S=0时，才考虑将约束相关的项加入协态变量微分方程。
约束S1=qα−Qα≤0S_{1}=q\alpha - Q_{\alpha} \leq 0S1​=qα−Qα​≤0
当S1=0S_{1} = 0S1​=0时
∂q∂r=∂(12ρVr2)∂r=12Vr2∂ρ∂r+ρ∂Vr∂r=12Vr2ρrrr+ρωˉE×Vr\begin{aligned} \frac{\partial q}{\partial\boldsymbol{r}} =& \frac{\partial(\frac{1}{2}\rho V^{2}_{r})} {\partial\boldsymbol{r}} = \frac{1}{2} V^{2}_{r} \frac{\partial\rho}{\partial{\boldsymbol{r}}} + \rho \frac{\partial V_{r}}{\partial{\boldsymbol{r}}} \\ =& \frac{1}{2} V^{2}_{r} \rho_{r} \frac{\boldsymbol{r}}{r} + \rho \bar{\omega}^{\times}_{E}\boldsymbol{V}_{r} \end{aligned} ∂r∂q​==​∂r∂(21​ρVr2​)​=21​Vr2​∂r∂ρ​+ρ∂r∂Vr​​21​Vr2​ρr​rr​+ρωˉE×​Vr​​

cos⁡α=1bT1Vr⇒−sin⁡α∂α∂V=∂(1bT1Vr)∂V⇒−sin⁡α∂α∂V=1Vr∂(1bTVr)∂V+1bTVr∂1Vr∂V⇒∂α∂V=1Vrsin⁡α(cos⁡α1Vr−1b)\begin{aligned} \cos{\alpha}=\boldsymbol{1}^{T}_{b}\boldsymbol{1}_{Vr} &\Rightarrow -\sin{\alpha}\frac{\partial{\alpha}}{\partial{\boldsymbol{V}}} = \frac{\partial{(\boldsymbol{1}^{T}_{b}\boldsymbol{1}_{Vr})}}{\partial{\boldsymbol{V}}} \\ &\Rightarrow -\sin{\alpha} \frac{\partial{\alpha}}{\partial{\boldsymbol{V}}} = \frac{1}{V_{r}} \frac{\partial{(\boldsymbol{1}^{T}_{b}\boldsymbol{V}_{r})}}{\partial{\boldsymbol{V}}} + \boldsymbol{1}^{T}_{b}\boldsymbol{V}_{r} \frac{\partial \frac{1}{V_{r}}}{\partial\boldsymbol{V}} \\ &\Rightarrow \frac{\partial{\alpha}}{\partial\boldsymbol{V}} = \frac{1}{V_{r}\sin{\alpha}}(\cos{\alpha} \boldsymbol{1}_{Vr} - \boldsymbol{1}_{b}) \end{aligned} cosα=1bT​1Vr​​⇒−sinα∂V∂α​=∂V∂(1bT​1Vr​)​⇒−sinα∂V∂α​=Vr​1​∂V∂(1bT​Vr​)​+1bT​Vr​∂V∂Vr​1​​⇒∂V∂α​=Vr​sinα1​(cosα1Vr​−1b​)​

∂α∂r=(∂V∂r)T∂α∂V=(ωˉE×)TqVrsin⁡α(cos⁡α1Vr−1b)\begin{aligned} \frac{\partial\alpha}{\partial\boldsymbol{r}} =& \left(\frac{\partial\boldsymbol{V}}{\partial\boldsymbol{r}}\right)^{T} \frac{\partial\alpha}{\partial\boldsymbol{V}} \\ =& (\bar{\omega}^{\times}_{E})^{T} \frac{q}{V_{r}\sin{\alpha}}(\cos{\alpha} \boldsymbol{1}_{Vr} - \boldsymbol{1}_{b}) \end{aligned} ∂r∂α​==​(∂r∂V​)T∂V∂α​(ωˉE×​)TVr​sinαq​(cosα1Vr​−1b​)​

∂S1∂r=∂(qα)∂r=α∂q∂r+q∂α∂r=12αρrVr2rr+αρωˉE×Vr+(ωˉE×)TqVrsin⁡α(cos⁡α1Vr−1b)\begin{aligned} \frac{\partial S_{1}}{\partial \boldsymbol{r}} =& \frac{\partial (q\alpha)}{\partial \boldsymbol{r}} = \alpha \frac{\partial{q}}{\partial{\boldsymbol{r}}} + q \frac{\partial{\alpha}}{\partial{\boldsymbol{r}}} \\ =&\frac{1}{2} \alpha \rho_{r} V^{2}_{r}\frac{\boldsymbol{r}}{r} + \alpha \rho \bar{\omega}^{\times}_{E} \boldsymbol{V}_{r} + (\bar{\omega}^{\times}_{E})^{T} \frac{q}{V_{r}\sin{\alpha}}(\cos{\alpha} \boldsymbol{1}_{Vr} - \boldsymbol{1}_{b}) \end{aligned} ∂r∂S1​​==​∂r∂(qα)​=α∂r∂q​+q∂r∂α​21​αρr​Vr2​rr​+αρωˉE×​Vr​+(ωˉE×​)TVr​sinαq​(cosα1Vr​−1b​)​

∂S1∂V=∂(qα)∂V=α∂q∂V+q∂α∂V=αρVr+qVrsin⁡(α)(cos⁡α1Vr−1b)\begin{aligned} \frac{\partial S_{1}}{\partial \boldsymbol{V}} =& \frac{\partial (q\alpha)}{\partial \boldsymbol{V}} = \alpha \frac{\partial{q}}{\partial{\boldsymbol{V}}} + q \frac{\partial{\alpha}}{\partial{\boldsymbol{V}}}\\ =& \alpha \rho \boldsymbol{V}_{r} + \frac{q}{V_{r}\sin(\alpha)} (\cos{\alpha} \boldsymbol{1}_{Vr} - \boldsymbol{1}_{b}) \end{aligned} ∂V∂S1​​==​∂V∂(qα)​=α∂V∂q​+q∂V∂α​αρVr​+Vr​sin(α)q​(cosα1Vr​−1b​)​

其中，ρr=∂ρ/∂r\rho_{r} = \partial{\rho}/\partial{r}ρr​=∂ρ/∂r。

约束S2=T−Tmax≤0S_{2}=T - T_{max} \leq 0S2​=T−Tmax​≤0
当S2=0S_{2} = 0S2​=0时，发动机的阀门根据以下公式调节
η=Tmaxm(t)g0/Tvac\eta = T_{max}m(t)g_{0}/T_{vac} η=Tmax​m(t)g0​/Tvac​

约束S3=q−qmax≤0S_{3}=q - q_{max} \leq 0S3​=q−qmax​≤0
文章中对约束S3S_{3}S3​的处理那块儿，公式有点没弄懂！！当S3=0S_{3} = 0S3​=0时，作者将动压约束转化为了对发动机阀门大小η\etaη的约束，联立以下四个方程[1]
V′=−1r3r+A+T1b+NT=η[Tvac+ΔT(r)]/(m(t)gE)Vr′=V′−ωˉE×V′−VwS3′=(1/(2r))ρrVr2rTV+ρVrTVr′\begin{aligned} \boldsymbol{V}^{\prime} =& -\frac{1}{r^{3}}\boldsymbol{r} + \boldsymbol{A} + T\boldsymbol{1}_{b} + \boldsymbol{N} \\ T = &\eta \left[ T_{vac} + \Delta T(r)\right]/(m(t)g_{E}) \\ \boldsymbol{V}^{\prime}_{r} =& \boldsymbol{V}^{\prime} - \bar{\omega}_{E} \times \boldsymbol{V}^{\prime} - \boldsymbol{V}_{w} \\ S^{\prime}_{3} =& (1/(2r))\rho_{r}V^{2}_{r} \boldsymbol{r}^{T}\boldsymbol{V} + \rho \boldsymbol{V}^{T}_{r}\boldsymbol{V}^{\prime}_{r} \end{aligned} V′=T=Vr′​=S3′​=​−r31​r+A+T1b​+Nη[Tvac​+ΔT(r)]/(m(t)gE​)V′−ωˉE​×V′−Vw​(1/(2r))ρr​Vr2​rTV+ρVrT​Vr′​​

可知
q′(t)=aq(x,1b)η(t)+bq(x,1b)q^{\prime}(t) = a_{q}(\boldsymbol{x},\boldsymbol{1}_{b})\eta(t) + b_{q}(\boldsymbol{x},\boldsymbol{1}_{b}) q′(t)=aq​(x,1b​)η(t)+bq​(x,1b​)

其中
aq(x,1b)=...a_{q}(\boldsymbol{x},\boldsymbol{1}_{b}) = ... aq​(x,1b​)=...

bq(x,1b)=...b_{q}(\boldsymbol{x},\boldsymbol{1}_{b}) = ... bq​(x,1b​)=...

另δ>0\delta > 0δ>0则有
q(t+δ)≈q(t)+q′(t)δ\begin{aligned} q(t + \delta) \approx q(t) + q^{\prime}(t)\delta \end{aligned} q(t+δ)≈q(t)+q′(t)δ​

根据动压约束q(t+δ)−qmax≤0q(t + \delta) - q_{max} \leq 0q(t+δ)−qmax​≤0
η(t)≤qmax⁡−q(t)−bqδaqδ≜ηq\begin{aligned} \eta(t) \leq \frac{q_{\max }-q(t)-b_{q} \delta}{a_{q} \delta} \triangleq \eta_{q} \end{aligned} η(t)≤aq​δqmax​−q(t)−bq​δ​≜ηq​​

将由约束S2S_{2}S2​求得的η\etaη记为ηprb\eta_{\mathrm{prb}}ηprb​，η\etaη的最小值是ηmin⁡\eta_{\min }ηmin​，那么同时考虑约束2和3，η\etaη的取值方法如下[1]
η={ηprb,if ηq>ηprbηq,if ηmin⁡≤ηq≤ηprbηmin,if ηq<ηmin⁡\begin{aligned} \eta = \left \{ \begin{array}{lll} {\eta_{\mathrm{prb}},} & {\text { if }} & {\eta_{q}>\eta_{\mathrm{prb}}} \\ {\eta_{q},} & {\text { if }} & {\eta_{\min } \leq \eta_{q} \leq \eta_{\mathrm{prb}}}\\ {\eta_{\mathrm{min}},} & {\text { if }} & {\eta_{q}<\eta_{\min }}\end{array}\right. \end{aligned} η=⎩⎨⎧​ηprb​,ηq​,ηmin​,​ if  if  if ​ηq​>ηprb​ηmin​≤ηq​≤ηprb​ηq​<ηmin​​​

密度、声速、推力对地心距求导

密度
已知是空气密度是关于运载器所在高度的函数，即ρ1=f(h1)\rho_{1} = f(h_{1})ρ1​=f(h1​)
ρr=∂ρ∂r=∂ρ∂(1+h)=∂ρ∂h(归一化方程)\begin{aligned} \rho_r = \frac{\partial \rho}{\partial r} =& \frac{\partial \rho}{\partial(1 + h)} = \frac{\partial \rho}{\partial h} \end{aligned}\quad \text{(归一化方程)} ρr​=∂r∂ρ​=​∂(1+h)∂ρ​=∂h∂ρ​​(归一化方程)

未归一化密度ρ1=ρ0ρ\rho_{1} = \rho_{0}\rhoρ1​=ρ0​ρ，单位kg/m3kg/m^{3}kg/m3，未归一化运载器高度h1=R0h/1000h_{1} = R_{0}h / 1000h1​=R0​h/1000，单位kmkmkm
ρr=R01000ρ0∂ρ1∂h1\begin{aligned} \rho_{r} = \frac{R_{0}}{1000\rho_{0}} \frac{\partial \rho_{1}}{\partial h_{1}} \end{aligned} ρr​=1000ρ0​R0​​∂h1​∂ρ1​​​

声速
声速是关于运载器所在高度的函数，即Vs1=f(h1)V_{s1} = f(h_{1})Vs1​=f(h1​)，此f(h)f(h)f(h)非计算密度的f(h)f(h)f(h)。未归一化声速Vs1=R0gEVsV_{s1} = \sqrt{R_{0}g_{E}} V_{s}Vs1​=R0​gE​​Vs​，单位m/sm/sm/s，未归一化运载器高度h1=R0h/1000h_{1} = R_{0}h / 1000h1​=R0​h/1000，单位kmkmkm
∂Vs∂r=R01000R0gE∂Vs1∂h1\begin{aligned} \frac{\partial V_{s}}{\partial r} = \frac{R_{0}}{1000\sqrt{R_{0}g_{E}}}\frac{\partial V_{s1}}{\partial h_{1}} \end{aligned} ∂r∂Vs​​=1000R0​gE​​R0​​∂h1​∂Vs1​​​

推力
文献中计算推力的公式T=η[Tvac+ΔT(r)]/(m(t)g0)T = \eta[T_{vac} + \Delta T(r)] / (m(t)g_{0})T=η[Tvac​+ΔT(r)]/(m(t)g0​)，其中ΔT(r)\Delta T(r)ΔT(r)怎么求呢？由文献[4]可知，在弹道计算中，通常采用下式计算推力
P=P0+Se(p0−pH)P = P_{0} + S_{e}(p_{0}-p_{H}) P=P0​+Se​(p0​−pH​)

其中，PPP是发动机推力，P0P_{0}P0​是发动机地面推力，SeS_{e}Se​是喷口载面积，p0p_{0}p0​是地面大气压，pHp_{H}pH​是发动机所在高度大气压。因此ΔT(r)=−SepH\Delta T(r) = -S_{e}p_{H}ΔT(r)=−Se​pH​，进而可得
Tr=∂T∂r=−ηSem(t)gE∂pH∂rT_{r} = \frac{\partial{T}}{\partial{r}} = -\frac{\eta S_{e}}{m(t)g_{E}} \frac{\partial p_{H}}{\partial{r}} Tr​=∂r∂T​=−m(t)gE​ηSe​​∂r∂pH​​

求解最优体轴矢量1b∗\boldsymbol{1}^{*}_{b}1b∗​

在求Y0Y_{0}Y0​过程中，我们已经求出一条满足终端约束的上升轨迹，根据该过程中求得的各量，求解1b∗\boldsymbol{1}^{*}_{b}1b∗​
计算Φ\PhiΦ
根据以下公式可以求解Φ\PhiΦ的绝对值[1]
Φ=1pVT1Vr\Phi=\mathbf{1}_{p_{V}}^{T} \mathbf{1}_{V_{r}} Φ=1pV​T​1Vr​​

再根据文章中的方法确定Φ\PhiΦ的正负。

用牛顿迭代法求解α\alphaα
∂H/∂α=0\partial{H}/\partial{\alpha} = 0∂H/∂α=0可求得下式[1]
tan⁡(Φ−α)(T−A+Nα)−(Aα+N)=0\tan (\Phi-\alpha)\left(T-A+N_{\alpha}\right)-\left(A_{\alpha}+N\right)=0 tan(Φ−α)(T−A+Nα​)−(Aα​+N)=0

要使用牛顿迭代法求解α\alphaα还需求解∂2H/∂α2{\partial^{2}H}/\partial{\alpha^{2}}∂2H/∂α2
∂2H∂α2=A−T−Nαcos⁡2(Φ−α)+tan⁡(Φ−α)(−Aα+Nαα)−(Aαα+Nα)\frac{\partial^{2}H}{\partial{\alpha^{2}}} = \frac{A-T-N_{\alpha}}{\cos^{2}{(\Phi - \alpha)}} + \tan{(\Phi - \alpha)}(-A_{\alpha}+N_{\alpha\alpha}) - (A_{\alpha\alpha} + N_{\alpha}) ∂α2∂2H​=cos2(Φ−α)A−T−Nα​​+tan(Φ−α)(−Aα​+Nαα​)−(Aαα​+Nα​)

其中，Aαα=∂2A/∂α2,Nαα=∂2N/∂α2A_{\alpha\alpha}={\partial^{2}A}/\partial{\alpha^{2}},N_{\alpha\alpha}={\partial^{2}N}/\partial{\alpha^{2}}Aαα​=∂2A/∂α2,Nαα​=∂2N/∂α2。仿真时，气动力AAA和NNN被拟合成攻角α\alphaα的二次多项式形式，很容易求得气动力对攻角的一阶和二阶偏导数。求得α\alphaα后，根据下式求解最优体轴矢量[1]
1b∗=(sin⁡αsin⁡Φ)1pV+[cos⁡α−cos⁡Φcos⁡(Φ−α)sin⁡2Φ]1Vr\mathbf{1}_{b}^{*}=\left(\frac{\sin \alpha}{\sin \Phi}\right) \mathbf{1}_{p_{V}}+\left[\frac{\cos \alpha-\cos \Phi \cos (\Phi-\alpha)}{\sin ^{2} \Phi}\right] \mathbf{1}_{V_{r}} 1b∗​=(sinΦsinα​)1pV​​+[sin2Φcosα−cosΦcos(Φ−α)​]1Vr​​

有限差分法解TPBVPs

两点边值问题中的动力学方程由运载器的动力学方程和协态变量微分方程组成[1]
r′=VV′=−1r3r+A+T1b+N\begin{aligned} \boldsymbol{r}^{\prime} =& \boldsymbol{V} \\ \boldsymbol{V}^{\prime} = -\frac{1}{r^{3}}\boldsymbol{r} + \boldsymbol{A} &+ T\boldsymbol{1}_{b} + \boldsymbol{N} \\ \end{aligned} r′=V′=−r31​r+A​V+T1b​+N​

pr′=1r3pV−[3apvbr4−apvn(Tr−Aρr+12VrCρVs2CAMach∂Vs∂r)+apvn(Nρr−12VrCρVs2CNMach∂Vs∂r)]rr+CρωˉE×{apvb[(CA+12VrVsCAMach)Vr+12CAαVr2∂α∂V]−apvn[(CN+12VrVsCNMach)Vr+12CNαVr2∂α∂V]}pV′=−pr+Cρ[apvb(CA+12VrVsCAMach)−apvn(CN+12VrVsCNMach)+12(apvbCAα−apvnCNα)cos⁡αsin⁡α]×Vr−CρVr2sin⁡α(apvbCAα−apvnCNα)1b\begin{aligned} \boldsymbol{p}^{\prime}_{r} =& \frac{1}{r^{3}}\boldsymbol{p}_{V} - \lbrack \frac{3a_{pvb}}{r^{4}} -a_{pvn} ( T_{r}-A_{\rho r} + \frac{1}{2V_{r}}C_{\rho}V^{2}_{s}C_{A_{Mach}}\frac{\partial V_{s}}{\partial r} ) \\ &+ a_{pvn}(N_{\rho r} - \frac{1}{2V_{r}}C_{\rho}V^{2}_{s}C_{N_{Mach}}\frac{\partial V_{s}}{\partial r})]\frac{\boldsymbol{r}}{r} \\ &+ C_{\rho}\bar{\omega}_{E} \times \{ a_{pvb} [(C_{A} + \frac{1}{2V_{r}}V_{s}C_{A_{Mach}}) \boldsymbol{V}_{r} + \frac{1}{2}C_{A_{\alpha}}V^{2}_{r} \frac{\partial \alpha}{\partial \boldsymbol{V}} ] \\ &- a_{pvn}[(C_{N} + \frac{1}{2V_{r}}V_{s}C_{N_{Mach}}) \boldsymbol{V}_{r} + \frac{1}{2}C_{N_{\alpha}}V^{2}_{r} \frac{\partial \alpha}{\partial \boldsymbol{V}} ] \} \\ \boldsymbol{p}^{\prime}_{V} =& -\boldsymbol{p}_{r} + C_{\rho}[a_{pvb}(C_{A} + \frac{1}{2V_{r}}V_{s}C_{A_{Mach}}) \\ &- a_{pvn}(C_{N} + \frac{1}{2V_{r}}V_{s}C_{N_{Mach}}) + \frac{1}{2}(a_{pvb}C_{A_{\alpha}} - a_{pvn}C_{N_{\alpha}})\frac{\cos{\alpha}}{\sin{\alpha}}] \\ & \times \boldsymbol{V}_{r} - \frac{C_{\rho}V_{r}}{2\sin{\alpha}}(a_{pvb}C_{A_{\alpha}} - a_{pvn}C_{N_{\alpha}})\boldsymbol{1}_{b} \end{aligned} pr′​=pV′​=​r31​pV​−[r43apvb​​−apvn​(Tr​−Aρr​+2Vr​1​Cρ​Vs2​CAMach​​∂r∂Vs​​)+apvn​(Nρr​−2Vr​1​Cρ​Vs2​CNMach​​∂r∂Vs​​)]rr​+Cρ​ωˉE​×{apvb​[(CA​+2Vr​1​Vs​CAMach​​)Vr​+21​CAα​​Vr2​∂V∂α​]−apvn​[(CN​+2Vr​1​Vs​CNMach​​)Vr​+21​CNα​​Vr2​∂V∂α​]}−pr​+Cρ​[apvb​(CA​+2Vr​1​Vs​CAMach​​)−apvn​(CN​+2Vr​1​Vs​CNMach​​)+21​(apvb​CAα​​−apvn​CNα​​)sinαcosα​]×Vr​−2sinαCρ​Vr​​(apvb​CAα​​−apvn​CNα​​)1b​​

其中，ρr=∂ρ/∂r,Tr=∂T/∂r,Cρ=ρ0ρSrefR0/m(t)\rho_{r}=\partial\rho/\partial r, T_{r}=\partial T/\partial r, C_{\rho}=\rho_{0}\rho S_{ref}R_{0}/m(t)ρr​=∂ρ/∂r,Tr​=∂T/∂r,Cρ​=ρ0​ρSref​R0​/m(t)，并且有
Aρr=Vr2SrefR0ρ0CAρr2m(t),Nρr=Vr2SrefR0ρ0CNρr2m(t)A_{\rho r} = \frac{V^2_{r}S_{ref}R_{0}\rho_{0}C_{A}\rho_{r}}{2m(t)}, N_{\rho r} = \frac{V^2_{r}S_{ref}R_{0}\rho_{0}C_{N}\rho_{r}}{2m(t)} Aρr​=2m(t)Vr2​Sref​R0​ρ0​CA​ρr​​,Nρr​=2m(t)Vr2​Sref​R0​ρ0​CN​ρr​​

CAα=∂CA∂α,CNα=∂CN∂α,CAMach=∂CA∂MachC_{A_{\alpha}} = \frac{\partial C_{A}}{\partial \alpha}, C_{N_{\alpha}} = \frac{\partial C_{N}}{\partial \alpha}, C_{A_{Mach}} = \frac{\partial C_{A}}{\partial Mach} CAα​​=∂α∂CA​​,CNα​​=∂α∂CN​​,CAMach​​=∂Mach∂CA​​

CNMach=∂CN∂Mach,apvb=pVT1b=pVcos⁡(Φ−α)C_{N_{Mach}} = \frac{\partial C_{N}}{\partial Mach}, a_{pvb} = \boldsymbol{p}^{T}_{V}\boldsymbol{1}_{b} = p_{V}\cos{(\Phi - \alpha)} CNMach​​=∂Mach∂CN​​,apvb​=pVT​1b​=pV​cos(Φ−α)

apvn=pVT1n=pVsin⁡(Φ−α)a_{pvn} = \boldsymbol{p}^{T}_{V}\boldsymbol{1}_{n} = p_{V}\sin{(\Phi - \alpha)} apvn​=pVT​1n​=pV​sin(Φ−α)

∂α∂V=1Vrsin⁡α(cos⁡α1Vr−1b)\frac{\partial \alpha}{\partial \boldsymbol{V}} = \frac{1}{V_{r}\sin{\alpha}}(\cos{\alpha \boldsymbol{1}_{Vr}} - \boldsymbol{1}_{b}) ∂V∂α​=Vr​sinα1​(cosα1Vr​−1b​)

当加入约束S1=qα−Qα=0S_{1}=q\alpha - Q_{\alpha} = 0S1​=qα−Qα​=0后，协态方程
p′=−∂H∂x−λqα∂S1∂x\begin{aligned} \boldsymbol{p}^{\prime} = -\frac{\partial H}{\partial \boldsymbol{x}} - \lambda_{q\alpha}\frac{\partial S_{1}}{\partial \boldsymbol{x}} \end{aligned} p′=−∂x∂H​−λqα​∂x∂S1​​​

初始条件(x0,p0)(\boldsymbol{x}_{0},\boldsymbol{p}_{0})(x0​,p0​)已知，6个终端条件由入轨条件和最优控制理论求得的横截条件组成。
将时间分成相等的MMM段，一共M+1M+1M+1个节点，则每一段的长度h=tf/Mh = tf/Mh=tf/M。我们取M=100M = 100M=100，那么两点边值问题就转换为如下方程组求解问题

E0=[y0(1:3)y0(4:6)]−[r0V0]E1=y1−y0−hf(t1−12,(y1+y0)2)E2=y2−y1−hf(t2−12,(y2+y1)2)...E100=y100−y99−hf(t100−12,(y100+y99)2)E101=Ψ(y100)\begin{aligned} E_{0} =& \begin{bmatrix} y_{0}(1:3) \\ y_{0}(4:6) \end{bmatrix} - \begin{bmatrix} \boldsymbol{r}_{0} \\ \boldsymbol{V}_{0} \end{bmatrix}\\ E_{1} =& y_{1} - y_{0} - h f(t_{1-\frac{1}{2}},\frac{(y_1+y_0)}{2}) \\ E_{2} =& y_{2} - y_{1} - h f(t_{2-\frac{1}{2}},\frac{(y_2+y_1)}{2}) \\ &... \\ E_{100} =& y_{100} - y_{99} - h f(t_{100-\frac{1}{2}},\frac{(y_{100}+y_{99})}{2}) \\ E_{101} =& {\Psi}(y_{100}) \end{aligned} \\ E0​=E1​=E2​=E100​=E101​=​[y0​(1:3)y0​(4:6)​]−[r0​V0​​]y1​−y0​−hf(t1−21​​,2(y1​+y0​)​)y2​−y1​−hf(t2−21​​,2(y2​+y1​)​)...y100​−y99​−hf(t100−21​​,2(y100​+y99​)​)Ψ(y100​)​

其中y0(1:3)y_{0}(1:3)y0​(1:3)表示y0y_{0}y0​向量的前三个元素组成的向量。记方程的解为Y=[y0T,y1T,y2T,...,y100T]TY = [y_{0}^{T},y_{1}^{T},y_{2}^{T},...,y_{100}^{T}]^{T}Y=[y0T​,y1T​,y2T​,...,y100T​]T，其中y=[rT,VT,prT,pvT]Ty = [\boldsymbol{r}^{T}, \boldsymbol{V}^{T}, \boldsymbol{p}^{T}_{r},\boldsymbol{p}^{T}_{v}]^{T}y=[rT,VT,prT​,pvT​]T。这里需要注意的是前面的推导中我们用过[rT,VT,pvT,−prT]T[\boldsymbol{r}^{T}, \boldsymbol{V}^{T},\boldsymbol{p}^{T}_{v}, -\boldsymbol{p}^{T}_{r}]^{T}[rT,VT,pvT​,−prT​]T这种顺序。

数值方法求∂E/∂Y\partial{E} / \partial{Y}∂E/∂Y
∂E∂Y=[...]1212×1212=[∂E0y0∂E0y1...∂E0y100∂E1y0∂E1y1...∂E1y100...∂E101y0∂E101y1...∂E101y100]\begin{aligned} \frac{\partial{E}}{\partial{Y}} = [...]_{1212\times1212} = \begin{bmatrix} \frac{\partial{E_{0}}}{y_{0}} & \frac{\partial{E_{0}}}{y_{1}} & ... & \frac{\partial{E_{0}}}{y_{100}} \\ \frac{\partial{E_{1}}}{y_{0}} & \frac{\partial{E_{1}}}{y_{1}} & ... & \frac{\partial{E_{1}}}{y_{100}} \\ &... \\ \frac{\partial{E_{101}}}{y_{0}} & \frac{\partial{E_{101}}}{y_{1}} & ... & \frac{\partial{E_{101}}}{y_{100}} \end{bmatrix} \end{aligned} ∂Y∂E​=[...]1212×1212​=⎣⎢⎢⎢⎡​y0​∂E0​​y0​∂E1​​y0​∂E101​​​y1​∂E0​​y1​∂E1​​...y1​∂E101​​​.........​y100​∂E0​​y100​∂E1​​y100​∂E101​​​⎦⎥⎥⎥⎤​​

=[[I6,0]6×1206×12......06×12−I12×12−h∂f(t1−12,(y1+y0)2)∂y0I12×12−hf(,)y1012×12...012×12012×12−I12×12−h∂f(,)∂y1I12×12−hf(,)y2...012×12...............012×12...012×12−I12×12−h∂f(,)∂y99I12×12−hf(,)y10006×12......06×12∂Ψ∂y100]\begin{aligned} = \begin{bmatrix} [I_{6},\mathbf{0}]_{6\times12} & \mathbf{0}_{6\times12} & ... & ... & \mathbf{0}_{6\times12} \\ -I_{12\times12}-h\frac{\partial{f(t_{1-\frac{1}{2}},\frac{(y_1+y_0)}{2})}}{\partial{y_{0}}} & I_{12\times12}-h\frac{f(,)}{y_{1}} & \boldsymbol{0}_{12\times12} & ... & \boldsymbol{0}_{12\times12} \\ \boldsymbol{0}_{12\times12} & -I_{12\times12}-h\frac{\partial{f(,)}}{\partial{y_{1}}} & I_{12\times12}-h\frac{f(,)}{y_{2}} & ... & \mathbf{0}_{12\times12} \\ ... & ... & ... & ... & ... & \\ \mathbf{0}_{12\times12} & ... & \mathbf{0}_{12\times12} & -I_{12\times12}-h\frac{\partial{f(,)}}{\partial{y_{99}}} & I_{12\times12}-h\frac{f(,)}{y_{100}}\\ \mathbf{0}_{6\times12} & ... & ... & \mathbf{0}_{6\times12} & \frac{\partial{\mathbf{\Psi}}}{\partial{y_{100}}} \end{bmatrix} \end{aligned} =⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​[I6​,0]6×12​−I12×12​−h∂y0​∂f(t1−21​​,2(y1​+y0​)​)​012×12​...012×12​06×12​​06×12​I12×12​−hy1​f(,)​−I12×12​−h∂y1​∂f(,)​.........​...012×12​I12×12​−hy2​f(,)​...012×12​...​............−I12×12​−h∂y99​∂f(,)​06×12​​06×12​012×12​012×12​...I12×12​−hy100​f(,)​∂y100​∂Ψ​​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​​

其中，f(,)f(,)f(,)是对应时刻f(tk−12,(yk+yk−1)2)f(t_{k-\frac{1}{2}},\frac{(y_{k}+y_{k-1})}{2})f(tk−21​​,2(yk​+yk−1​)​)的简写，∂Ψ∂y100\frac{\partial{\mathbf{\Psi}}}{\partial{y_{100}}}∂y100​∂Ψ​之前有个小节推导过，是有解析形式的，除此以外其它的偏导数由数值的方法计算。

∂f(tk−12,(yk+yk−1)2)∂yk=∂f(tk−12,(yk+yk−1)2)∂((yk+yk−1)2)∂((yk+yk−1)2)∂yk\frac{\partial{f(t_{k-\frac{1}{2}},\frac{(y_k+y_{k-1})}{2})}}{\partial{y_{k}}} = \frac{\partial{f(t_{k-\frac{1}{2}},\frac{(y_k+y_{k-1})}{2})}}{\partial{(\frac{(y_k+y_{k-1})}{2}})} \frac{\partial{(\frac{(y_k+y_{k-1})}{2}})}{\partial{y_{k}}} ∂yk​∂f(tk−21​​,2(yk​+yk−1​)​)​=∂(2(yk​+yk−1​)​)∂f(tk−21​​,2(yk​+yk−1​)​)​∂yk​∂(2(yk​+yk−1​)​)​

高斯消元法求解搜索方向dj\boldsymbol{d}_{j}dj​
[∂E(Yj−1)∂Y]dj=−E(Yj−1)\begin{aligned} \left[ \frac{\partial{E}(Y_{j-1})}{\partial{Y}} \right] \boldsymbol{d}_{j} = -E(Y_{j-1}) \end{aligned} [∂Y∂E(Yj−1​)​]dj​=−E(Yj−1​)​

由于∂E∂Y\frac{\partial{E}}{\partial{Y}}∂Y∂E​的稀疏特性，我们可以用高斯消元法求解dj\boldsymbol{d}_{j}dj​

参考文献

[1] Lu P , Sun H , Tsai B . Closed-Loop Endoatmospheric Ascent Guidance[J]. Journal of Guidance, Control, and Dynamics, 2003, 26(2):283-294.
[2] Calise A J , Melamed N , Lee S . Design and Evaluation of a Three-Dimensional Optimal Ascent Guidance Algorithm[J]. Journal of Guidance, Control, and Dynamics, 1998, 21(6):867-875.
[3] Mcadoo S F J , Jezewski D J , Dawkins G S . Development of a method for optimal maneuver analysis of complex space missions[J]. Nasa Sti/recon Technical Report N, 1975, 75.
[4] 王志刚, 施志佳. 远程火箭与卫星轨道力学基础[M].西安：西北工业大学出版社，2006.
[5] 傅瑜. 升力式天地往返飞行器自主制导方法研究[D]. 哈尔滨工业大学, 2012.

陆平老师论文Closed-Loop Endoatmospheric Ascent Guidance读后总结相关推荐

1. 计算机网络教学方式探讨论文,学生老师论文,关于关于高中计算机网络教学效率提升相关参考文献资料-免费论文范文...

   导读:本论文是一篇免费优秀的关于学生老师论文范文资料,可用于相关论文写作参考. (山东省新泰市第一中学 山东新泰 271200) 摘 要:随着计算机网络的发展和快速普及,计算机网络教学已经逐步进入高中 ...

2. 闵帆老师论文写作课心得体会——如何写好一篇论文

   目录 前言 一.学术论文 二.Latex排版及数学公式 三.关于英文 四.题目--论文的第一印象 五.摘要及关键词 六.引言--完整的故事 七.文献综述要重视 八.伪代码--论文的核心之一 九.实验- ...

3. 【读论文】Loop Closure Detection for Visual SLAM Systems Using Convolutional Neural Network

   [读论文]Loop Closure Detection for Visual SLAM Systems Using Convolutional Neural Network 发表于2017年,作者是南 ...

4. 解决ajax异步请求数据后swiper不能循环轮播(loop失效)问题、滑动后不能轮播的问题。

   解决ajax异步请求数据后swiper不能循环轮播(loop失效)问题.滑动后不能轮播的问题. 参考文章: (1)解决ajax异步请求数据后swiper不能循环轮播(loop失效)问题.滑动后不能轮播 ...

5. 缺陷检测相关论文阅读总结（记录自己读过的论文主要内容/Ideas）

   缺陷检测相关论文阅读总结(记录自己读过的论文主要内容) Attention!!! 点击论文题目即可访问原文or下载原文PDF文件: 每篇文章的内容包含:内容总结.文章Ideas: 更多关于缺陷检测以及 ...

6. word中给论文标参考文献和文献序号变化后的自动更新

   @word中给论文标参考文献和文献序号变化后的自动更新TOC 论文中参考文献的标注问题,整理出方法后面的童鞋以后避坑 1.文献采用word自动编号,菜单栏单击插入,选择交叉引用,选择相应的文献插入,如 ...

7. 普通人应该怎样做学术、写论文？——对于大多数本硕在读学生来说，是没有能力写一篇突破人类知识界限的论文的

   普通人应该怎样做学术.写论文? --对于大多数本硕在读学生来说,是没有能力写一篇突破人类知识界限的论文的. 然而,评奖.保研.升学都push着我们需要发论文.用论文去证明自己的工作和贡献.所以,作为普 ...

8. 【论文写作】闵帆老师论文写作课程心得体会30篇

   我的博客一直记录代码,还是第一次记录心得体会.谢谢我闵帆老师.此次博客记录了我在这学期上闵帆老师<论文写作>后的一些心得体会.在这节课上,我学习了写论文的注意事项.写论文所用的工具还有论文 ...

9. 【Exploring and Thinking】——闵帆老师论文写作课程学习心得

   开始之前先给自己打个鸡血吧: A journey of a thousand miles begins with a single step. 千里之行始于足下. 文章标题 前言 写学术论文前的准备 ...

最新文章

1. FPGA篇（三）基于FPGA的几种排序算法
2. U3D开发中关于脚本方面的限制-有关IOS反射和JIT的支持问题
3. MySQL中索引的分类和基本操作
4. 计算机专业运动会口号,运动会口号押韵有气势 计算机系霸气口号
5. mongodb mysql json数据_使用MongoDB与MySQL有很多JSON字段？
6. php在html中生成option,使用PHP可以将HTML SELECT/OPTION值设为NULL吗？
7. 电子科大计算机2014级,电子科大-计算机-操作系统实验报告-2014级.docx
8. 宁波大学2020计算机技术复试线,宁波大学2020年考研复试分数线
9. bzoj3620 似乎在梦中见过的样子
10. 参考资料：图片效果展示
11. jeesite实战(三十六)——非status的其他属性In条件查询
12. VsCode之在vue中HTML代码使用自动补全插件
13. 查看编译war包的jdk版本
14. XJOI恺撒加密术1级19段
15. linux蜂鸣器驱动
16. 实战接入腾讯云日志服务
17. matlab 怎么求直线斜率,matlab中如何求近似（不平滑）直线的斜率
18. 该内存不能为 read/written解决办法
19. ARM TrustZone ----ARM信任区
20. 你必须知道的最好的开源WEB 资源

热门文章

1. linux如何重新分区
2. linux网卡slave状态,生产环境中linux bonding 主备模式slave网卡切换的方法
3. J Magn Reson Imaging：磁共振指纹(MRF)动脉自旋标记(ASL)的灌注特性估计
4. 《Activiti工作流框架》专题(一)-Activiti工作流框架基础入门
5. 软件测试工程实训综合管理平台
6. str(n)cpy的注意事项以及memset的简单使用
7. 《Python深度学习》第一章笔记
8. 加壳工具WinLicense使用教程，以v2.3.9.0为例
9. 如何用2SC5200晶体管制作音频放大器
10. ubuntu常用命令大全（转）