全国大学生数学建模竞赛题获奖论文

承诺书

我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）：

我们的报名参赛队号为（8位数字组成的编号）：

所属学校（请填写完整的全名）：

参赛队员(打印并签名) ：1.

2.

3.

指导教师或指导教师组负责人(打印并签名)：

（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。）

日期：年月日

赛区评阅编号（由赛区组委会评阅前进行编2014高教社杯全国大学生数学建模竞赛

编号专用页

赛区评阅编号（由赛区组委会评阅前进行编号）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

嫦娥三号软着陆轨道设计与控制策略

摘要

在进行载人登月或月面勘测时，需要使飞行器实现月面软着陆以保证人员或设备的安全，但关键问题是着陆轨道与控制策略的设计。本文通过物理中的力学知识以及协方差分析等方法，进行了合理的轨道设计及优化。

针对问题一，对嫦娥三号软着陆的轨道以及六个阶段进行分析，通过机械能守恒定律、开普勒三定律等力学知识，建立了动力学模型。因为嫦娥三号绕月球运行的轨道是偏心率很小的椭圆，所以可以近似看作圆周轨道运动，然后迅速减速进入椭圆轨道，由动能改变量等于重力势能改变量及开普勒第二定律，算出着陆器在近月点与远月点的速度大小分别是s和s，方向沿运行轨道切线方向。然后根据质点运动学知识求出近月点与着陆点水平距离，进而利用坐标正反算软件算出近月点的经纬度为,，进而由空间解析几何知识得出了远月点的坐标（，，），并采用Matlab软件画出近月点和远月点在三维空间中的示意图。

针对问题二，嫦娥三号着陆轨道近月点和远月点的位置以及相应速度的大小与方向确定后，需要描述的是嫦娥三号软着陆过程中在不同阶段的运动状态，进而确定出嫦娥三号着陆轨道。由于轨道的设计要以燃料消耗最优为出发点，所以

可以在Matlab的平台上采用SFLA[]1优化方法，建立优化模型。将软着陆的动力

学方程做归一处理，经过将软着陆轨道离散化，从而将轨道优化问题转变为参数优化问题。通过仿真实验，作出嫦娥三号在软着陆过程中径向速度、推力控制角以及月心距的变化曲线，即设计出了最优软着陆轨道。

针对问题三，在一般的发射任务中，软着陆轨道修正都会选取将着陆器送到满足要求的目标轨道上（例如形成满足条件的环月轨道）的方式，而并非送到目标点上，这是因为后者需要选择合适的目标点使得轨道修正的能耗不会太大，且着陆器还需要在目标点进行变轨从而使得实际轨道与标称轨道重合。考虑到轨道参数的误差相对于轨道参数的标称值是小量，因此可以将轨道运动方程进行线性化，从而得到能够反映轨道参数偏差量的传播关系的误差方程。因此该问题采用协方差分析的方法，将着陆器发动机的一些技术指标的误差作为待考察的随机误差源，通过考虑嫦娥三号的运动轨迹进而评估位置误差和速度误差对飞行轨道的影响。最后，通过对变量F的敏感性分析，当F在1500N到6000N时，位移变化较小，运动轨迹影响较小，因此变量F对运动轨迹不敏感；当F在6000N到7500N时，位移变化较大，对运动轨迹影响较大，因此变量F对运动轨迹比较敏感。

通过仿真计算等验证，说明了建立的模型和计算结果都是可靠的。

关键词：动力学模型，轨道优化，混合蛙跳算法，协方差分析法

一、问题重述

嫦娥三号将在北京时间2013年12月14号在月球表面实施软着陆。嫦娥三号如何实现软着陆以及能否成功成为外界关注的焦点。目前，全球仅有美国、前苏联成功实施了13次无人月球表面软着陆。嫦娥三号的轨道和嫦娥二号一样，它从地球到月球的路程需要4~5天，到了月球以后，还不能直接登录月球，先得让月球“捕捉”嫦娥三号，使之成为月球的卫星，然后绕行。一开始沿着绕月球的大椭圆轨道运行，接着需要调整轨道，让它离月球越来越近，一直调整到降落轨道，再根据地面指令，在虹湾地区软着陆。

在月球上软着陆时不能用降落伞，因为月球是真空，降落伞毫无用处，探测器系统原来的初速度很大，加上月球的引力作用，下降的速度会越来越快，这时必须降低它的降落速度。在嫦娥三号的着陆器下方有一些发动机，可以产生向上的推力，减低它的下降速度。当它距月面100米高时，地球上的测控人员看不到现场的情况，因此要交给嫦娥三号自己去判断，从而选择相对平坦的地方降落。嫦娥三号从100米高的地方慢慢下降，落到距离月面4米高的地方关闭发动机，自由落体到月球表面，实现软着陆。其着陆轨道的基本要求是：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道，着陆轨道为近月点至着陆点，其软着陆过程分为6个阶段（主减速段、快速调整段、粗避障段、精避障段、缓速下降段、自由落体段），尽量减少软着陆过程中的燃料消耗。本文尝试解决以下问题：

问题一：确定着陆轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。

问题二：确定嫦娥三号着陆轨道和在6个阶段的最优控制策略。

问题三：对于设计的着陆轨道和控制策略做出相应的误差分析和敏感性分析。

二、问题分析

对嫦娥三号软着陆轨道的设计与控制为一个最优控制问题，要求保证准确地在月球预定区域内实现软着陆，根据软着陆过程中给定的6个阶段在关键点所处的状态设计嫦娥三号软着陆过程的运行轨道，并尽量减少燃料消耗。

(一) 问题一

要使嫦娥三号准确地在月球预定区域内实现软着陆，首先需要解决的是近月点位置的选择，根据附件2中嫦娥三号软着陆过程示意图及着陆过程中的主减速阶段和快速减速阶段计算出了着陆器着陆过程的水平位移，进而确定了近月点的位置，由空间解析几何知识得出远月点的位置，并采用Matlab软件画出示意图。考虑到嫦娥三号绕月球运行的轨道是偏心率很小的椭圆，所以可以近似看作圆周轨道运动，然后迅速减速进入椭圆轨道，又由机械能守恒定律及开普勒第二定律建立动力学模型算出在近月点与远月点的速度和方向。

(二) 问题二

嫦娥三号着陆轨道近月点和远月点的位置以及相应速度的大小与方向确定后，需要描述是嫦娥三号软着陆过程中在不同阶段的运动状态，进而确定出嫦娥三号着陆轨道。由于在整个着陆过程中力的方向和大小在不断的变化，如何确定

每个阶段的运行时间和位移以及相应阶段在关键点力的方向和大小，使每个阶段燃料消耗最少是这个问题的关键。我们将软着陆的动力学方程做归一处理，经过将软着陆轨道离散化，从而将轨道优化问题转变为参数优化问题，最后设计了蛙跳法作为优化方法。

(三) 问题三

本问题需要解决的是对于问题一、问题二轨道的设计进行误差分析和敏感性分析。由于前两个问题中的模型比较理想化，未考虑着陆器的质量随时间的变化，与实际情况中嫦娥三号的速度与位移存在较大的误差，为了减小误差，使它与实际轨道较为相符，采用协方差的方法，通过考虑嫦娥三号的运动轨迹继而评估位置误差和速度误差对飞行轨道的影响。推力为定值时，在不同数值的推力作用下，水平位移受其影响在不断变化，因此判断F对水平位移的敏感性。

三、模型假设

1.月球的偏心率为0；

2.由环月圆轨道进入椭圆轨道推力做的功较小，忽略不计；

3.忽略月球自转和倾斜角；

4.无穷远处万有引力势能为零；

5.月球半径取月球平均半径；

6.地球对嫦娥三号的引力忽略不计；

7.月球表面无大气层；

四、符号说明

五、模型建立与求解

（一）问题一 1问题分析

嫦娥三号在进入椭圆轨道前在距月面100的环月轨道做匀速圆周运动，然后迅速减速进入近月点高度约15公里，远月点高度约为100公里的椭圆轨道，由于绕月球运行的轨道是偏心率很小的椭圆，可以近似看作圆周轨道运动，由能量守恒定律及开普勒第二定律，算出着陆器在近月点与远月点的速度，其方向沿运行轨道切线方向。对嫦娥三号软着陆过程进行分析，将嫦娥三号在软着陆过程中受到的推力分解为竖直方向与水平方向，大小视为恒定的平均值，通过质点运动学的公式算出主减速阶段的时间，代入到水平方向的运动方程中算出水平位移并利用坐标正反算软件算出近月点的经纬度，进而由空间解析几何知识得出了远月点的位置，并采用Matlab 软件画出近月点和远月点在三维空间中的示意图。 2动力学模型的建立与求解 模型一

在环圆轨道进入椭圆轨道时，运用能量守恒定律及开普勒第二定律可知：飞船在向心力的作用下沿轨道运行时，其掠面速度是恒定的，可列出两个方程：

R g m v m E v p ??-??=+??月近远2

22

1m 21 （1） 近远v r ?=?R v （2）

若取无穷远处为万有引力势能的零点，则嫦娥三号在椭圆轨道上运行的势能

为：

r

GMm

dr r GMm E P -

==?∞r 2 （3） 又因

月月

g m R GM ?=2

m

（4） 故有r R g E P 2

m 月

月??-= （5）

由此解得s km v /692.1=近，s km v /633.1=远，其方向分别与运行轨道相切。 在嫦娥三号卫星软着陆过程中，可以将主减速阶段看作类平抛运动，水平方向初速度为近月点速度，即s km v /692.1=近，竖直方向速度为=0竖v 0s m k ，完成主减速阶段后进入快速调整阶段，此时的速度57/m s 可近似看作为竖直方向的速度，故竖直方向的速度为s m v /571=竖，主减速阶段竖直方向距离变化量m h 120001=?，对竖直方向进行分析，因为飞船的推力大小和方向随时间在不断变化，竖直方向的分力大小也在改变，为了方便计算，将嫦娥三号在软着陆过程

中受到的推力分解为竖直方向与水平方向，大小视为恒定的平均值，于是竖直方向的力大小恒定。由机械能守恒定律和牛顿第二定律，列出下式：

)(2

12

021111竖竖竖月v v m h F h mg -=??-? （6）

11-ma F mg =竖月 （7）

月

g h t 1

12?=

（8） 由式（6）、（7）、（8）得出s t s m a N F 421,/135.0,1.35951211===竖。

查阅相关资料知嫦娥三号着陆过程的快速调整阶段时间s t 202=。在水平阶段，推力的水平分力可以看作一平均值且大小恒定。在主减速阶段和快速调整阶段，水平方向可以近似看作匀减速直线运动，末速度s v /m 0=水。由运动学公式，得：

x a v v 22

22=-近水 (10)

求得km x s m a 8.372/84.322==，。

近月点投影在月球表面的点B 的经纬度即为近月点的经纬度。已知着陆点A 的经纬度，将其代入到坐标正反算和经纬度与XY 转化软件，结合余弦定理

2

2

222cos R X R R ?-+=α (11)

得近月点B 点的经纬度为,，故近月点的位置为距月球表面15公里处，经纬度为，。软件部分实现过程见如下图1、图2、图3。

图 1

图 2

图 3

以月球的球心为坐标原点，即球心O 坐标为（0，0，0），以赤道平面为XOY 平面，垂直赤道平面过球心的轴为Z 轴，建立空间直角坐标系，将着陆点所在的象限为第一象限，并由式（11）写出球面上各点对应点在空间直角坐标系下的坐标

θ

?θ?

θsin sin cos cos cos R Z R Y R X === (12)则有B ，，，如图4，设远月点C 的坐标为(x ，y ，z)。

图4 近月点与着陆点位置

则由利用空间解析几何知识得方程组

k z y x =-=-=-90.554018.1076034.11670 (13) 2222)90.554()18.1076()34.1167()013.173715100(-+-+-=++z y x (14) 解方程组求出远月点C 点的坐标为（，，）。

（二） 问题二 1 问题分析

假设不考虑质量随时间的变化，对每个阶段进行分析：在第一阶段和第二阶段嫦娥三号做类平抛运动；第三阶段做匀减速运动，使速度减小为零；第四阶段做自由落体并通过其四周的发动机选择适宜方位；第五阶段打开主发动机，做匀减速运动；第六阶段自由落体到预定着陆点。

由于以上模型都是在假设条件下的理想化模型，不符合实际情况，应对着陆的运行轨道进行优化。在椭圆轨道的近月点制动发动机点火，以抵消登月器的初始动能和势能，从而使得着陆器在水平速度被基本抵消之后相对月面速度降为0，以垂直姿态降落到月面，实现所谓的软着陆。着陆器的大部分燃料消耗在制动发动上，所以月球软着陆轨道的优化控制应以燃耗最优性为出发点，同时要兼顾降落到月面时的安全性。 2 动态优化模型的建立 模型二

在模型一中已得出主减速阶段和快速调整阶段的水平位移x=，时间

s t t t 44121=+=，快速调整后速度为s m v /7.592=。在粗避障阶段，嫦娥三号做

匀减速运动，初速度为s m v /7.592=，末速度为零，高度m h 23002=?，有

232

22h a v ?=，332t a v =，求解得出结果2330.77/,77.53a m s t s ==。

在精避障阶段，在竖直方向受到推力和重力合力作用下的自由落体运动，初速度为0m/s ，降落高度m h 703=?，水平方向周围小发动机调整位置，粗步避开大陨石坑。得

23

44544

12h a t v a t ?

?=??

?=? 求解方程组得s m v s t /12.15,26.954==。

缓速下降阶段，推力向上，做匀减速运动，末速度为零，高度m h 264=?。

255455

2v a h v t t ?=??

?

=??

求解方程组得算出s t s m a 45.3,/38.4525==

自由落体阶段，高度m h 45=?。有652

1

t g h 月=

?，得出s t 90.46=。 模型三、混合蛙跳模型

假设着陆轨道在纵向平面内，建立图5所示的平面极坐标系

图 5

其中，坐标原点在月心O ，y 轴指向初始轨道近月点，x 轴指向登月器开始运动方向。软着陆的动力学方程描述如式（15）：

2

2sin cos 2F v r m r r v F v mr r F m I

????

??βμ?=-+?ω???=?

?

θ=ω?

?βω

?ω=--

?

?=-?

? （15） 月球软着陆轨道是一个服从两点边值约束条件问题，由于起点处于霍曼转移轨道

的近月点，故满足的初始条件为：

00/v m s =，01753r km =，00rad γ=，409.6510/rad s -ω=?，终端约束条件为软

着陆，即需满足：0/f v m s =，f r R =，0/r rad s ω=，其中R=。

待优化变量为制动发动机的推力方向角()t β和终端时刻f t ，优化的性能指标

为在满足上述初始条件和终端约束条件的前提下，着陆过程中消耗的燃料最少，即：因为在发动机推力大小可调变的前提下得出的月球探测器软着陆最优制动方案，实际上是一个发动机推力大小恒定、推力与速度夹角变化的制动过程[]6，因此式（16）描述的性能指标也相当于令终端时刻f t 最短。

混合蛙跳算法模拟青蛙在寻找食物时，按族群分类进行思想传递的过程，将全局信息交换和局部深度搜索相结合，局部搜索使得思想在局部个体间传递，混合策略使得局部间的思想得到交换。通过这种全局信息交换和局部深度探索，算法能够跳出局部极值点，朝着全局最优的方向进行。

混合蛙跳算法的搜索优化过程如下：对于N 维组合优化问题，可将每只青蛙表示为()n x x x ,...,,flog 21=,其中i x 为第i 维的可能解；在进入迭优化之前，随机产生m 只青蛙组成初始种群，以后在每次进化迭代过程中，首先根据青蛙的适应值大小将所有青蛙按降序排列；然后再将重新排列后的青蛙分到n 个群体中，分配的方法为：第一只青蛙进入第一个蛙群，第二只青蛙进入第二个蛙群，第n 只青蛙进入第n 个群体，第n+1只青蛙进入到第一个蛙群，第n+2只青蛙进入到第二个蛙群，依此类推，只到所有青蛙分配完毕；当蛙群分好之后，针对每个种群中的最差青蛙，分别提高它们的适应值；最后将所有青蛙混合，进入下一次迭代循环。 3模型求解

将青蛙代入轨道动力学方程中解算得到终端时刻着陆器的状态，最后将末端状态代入式（16）

-

--

-

-

---+-+-+=f f f f f f f T r t r r r t v t v J J )())(())((22

2

1α??α （16）

其中，21,αα为惩罚因子。

计算性能指标，由于SFLA 的适应度取大值，而式（16）取小值为优，所以最终需对式（16）的计算值取负号。在得到每个青蛙的适应度后就可以按照上述的SFLA []1优化迭代步骤进行迭代优化了。

本文在Matlab 平台上采用SFLA 优化月球软着陆轨道做了仿真实验，以验证和分析此方法的应用可行性和优化效果。为了更方便地看出优化效果，本文采用文献[2]仿真的着陆器参数为：推力F=1350N ，比冲I=?s,初始质量0m =600kg 。轨道离散化段数取N=15,则待优化量为16个控制方向角i β和一个终端时刻f t 共17维参数。

SFLA 的参数值设置为：青蛙总数m=260，蛙群数n=20，外循环次数0N =100，内循环次数10i =N ，每只青蛙由17维待优化参数组成。惩罚因子

()

02i/11max 1,10N -?=，211.8?=?，惩罚因子设计成时变的目的主要是在迭代后期加强对不满足终端速度和终端位置的惩罚力度。拟合函数的阶数取n=4,求解微分方程采用龙格—库塔4阶法，计算步长取为1s 。表1将本文的平均优化结果和

文献[2]的优化结果做了对比，然后给出了某次仿真的优化结果图：

本文优化结果文献[2]优化结果约束要求终端径向速度（m/s）0

终端月心距（km）1738 终端角速度（rad/s）0

消耗燃料（kg）最少

喷气时间（s）无

表1 优化结果

图6 径向速度变化曲线

图7 推力控制角变化曲线

图8 月心距变化曲线

（三）问题三 1 问题三的分析

相比发射地球卫星而言，着陆器在地月空间飞行距离较大、飞行时间较长，由于力学模型和测控误差等的影响，实际转移轨道相对目标轨道会存在一些偏离，而且随着着陆器在轨运行，这些偏差还会被不同程度的放大。实际上，在一般的发射任务中，轨道修正都会选取将着陆器送到满足要求的目标轨道上（例如形成满足条件的环月轨道）的方式，而并非送到目标点上，这是因为后者需要选择合适的目标点使得轨道修正的能耗不会太大，且着陆器还需要在目标点进行变轨从而使得实际轨道在相空间中与标称轨道重合。考虑到轨道参数的误差之相对于轨道参数的标称值是小量，因此可以将轨道运动方程进行线性化，从而得到能够反映轨道参数偏差量的传播关系的误差方程。 2双二体模型建立[]3

误差分析：

反映轨道位置和速度误差的线性化方程如下

T r v g v r

r ?

?

→→→

→→→??=???????=??

??

（17） 其中→

→

??-=r R

g r g 3)(月， 且

2

22z y x r r r R ++= （18）

写成状态方程形式：

???

?? ?????????

?

?=?????

?

?

????→→→→?

?v r O I G r r 0 （19） 其中→

→??=

T

r

g

G 。

3模型求解 令

???

?

?

????=?

??

? ??=→→v r X I G F ,00 （20） 则式（18）变为

X F X ?=?

（21）

下面推导矩阵F 的表达式：

333333

()()()()()T T

x y z g u G r R r r g g g g r I r r r r r r r →

→

→→→?

???==-?????

??????---?-??????????

??月月月月 （22） 式中z y x r r r ,,是嫦娥三号在月心惯性坐标系里的轨道位置坐标。则

???

?

? ?

?

??-

-=→→2

333r r r I r g G T

月 （23）

()????

?? ?

?=?????

??=→

→222z y

z x

z z y y x

y z x y x x z y x z y x T

r r r r r r r r r r r r r r r r r r r r r r r （24） 将式（23），（24）代入（19），得：

003-13300033-13000333-11

000

00

0100

00001000

22

3555

22

355

522

3

）（）（）（月月月月月

月月月月R

r R g R r r g R r r g R r r g R r R g v R r r g R r r g R r r g R r R g F z

z y z x z y y y x y

x y

x x ---=

对式（20）积分，得到：

)0()(X e t X t F ?=?

式中

222

2

()()()()2!3!4!!

()!F t

i N

i i F t F t F t F t e

I F t n t F i ?=????=+?++++

+

?=∑

取前6阶截断，即：

)!(6

i t F e

i

i i

t

F ?=∑=?

得到计算误差方程的迭代方程：

)()(i t F i t X e t t X ?=?+

t F e ?相当于P 阵，由于误差方程是时变方程，因此每一步迭代都需要重新计算P

阵需要利用标称轨道参数数据。进而得到协方差矩阵迭代方程：

T

i i i P C P C i 1=+

求出结果，并画出图像，如表2，图9，图10，然后进行分析。

向月飞行轨道的初始轨道位置和速度误差由嫦娥三号的发射入轨精度决定，若着陆器在飞行中进行轨道修正，则经过轨道修正以后的轨道位置误差将由导航误差决定，速度误差将由姿态误差和制导误差决定。上述误差决定了轨道误差协方差分析的计算初始条件，表2给出了在不进行中途轨道修正情况下，在地心惯性坐标系里，初始轨道位置误差和初始速度误差对轨道终点的位置和速度误差的影响。图9和图10给出了着陆器从近月轨道入轨点开始至进入月球轨道为止的相应轨道位置和速度的总误差的时间历程。

表2 初始轨道位置和速度误差对轨道终点误差的影响

图9 轨道位置总误差时间历程

图10 速度总误差时间历

灵敏度分析：

运动学公式：

t

a v x a a h t m F mg a 水水

水竖

竖

竖2

1

m F 20+===

-=

求出位移，列表画图

由图像可知，推力为定值时，在不同数值的力推力作用下，水平位移在不断变化，F 在1500N 到6000N 时，位移变化较小，运动轨迹影响较小，因此变量F 对运动轨迹不敏感；当F 在6000N 到7500N 时，位移变化较大，对运动轨迹影响较大，因此变量F 对运动轨迹比较敏感。

六、模型评价

模型优缺点

混合蛙跳算法的优点是比较强的全局寻优能力和很高的优化精度，经SFLA 优化的轨道与传统方法优化的结果相差无几，某些指标甚至优于传统方法。同时

另一个优点是实现比较简单，没有其他智能算法的那样具有比较复杂的应用流程，也不需要设置太多的参数，因此它的应用操作上更加简单。双二体模型的具体一个优点是相对于其他智能传统算法，在应用上更加简单和易于操作。

七、参考文献

[1]郭景录，付平，登月软着陆轨道优化算法研究

[2]段佳佳，徐世杰，朱建丰，基于蚁群算法的月球软着陆轨迹优化[J],宇航学报，2008，29（2）：476-481.

[3]孙宝，荣思远，向月飞行轨道误差分析，哈尔滨工业大学，哈尔滨，150001.

[4]王大轶，李铁寿，马兴瑞.月球最优软着陆两点边值问题的数值解法[J].航天控制，2000，（3）：44-49.

[5]郗晓宁，等，月球探测器轨道设计[M].北京：国防工业出版社，2001.

附录

近月点与着陆点三维坐标系下的位置：

hold on;

a=;b=;

theta=0::2*pi;

phi=theta';

x=a*cos(phi)*cos(theta);

y=b*cos(phi)*sin(theta);

z=b*sin(phi)*ones(size(theta));

mesh(x,y,z);

shading interp

x=a*cos*1/180)*cos180);

y=b*cos*1/180)*sin180);

z=b*sin180);

plot3(x,y,z,'*')

x1=(a+15)*cos*1/180)*cos180);

y1=(b+15)*cos*1/180)*sin180);

z1=(b+15)*sin180);

plot3(x1,y1,z1,'b+')

推力与位移关系图：

x=[1500 3000 4500 6000 7500];

y=[ ];

plot(x,y,'b')