STK 与 Matlab 的接口及其 在复杂航天任务仿真分析中的应用 STK实验室 二零一一年 八月.

Presentation on theme: "STK 与 Matlab 的接口及其 在复杂航天任务仿真分析中的应用 STK实验室 二零一一年 八月."— Presentation transcript:

1 STK 与 Matlab 的接口及其 在复杂航天任务仿真分析中的应用 STK实验室 二零一一年 八月

2 主要内容 1 STK与Matlab接口简介 2 在复杂航天任务仿真分析中的应用 使用STK与Matlab接口应用场合与必要性
常用接口函数的使用 2 在复杂航天任务仿真分析中的应用 循环——星座对地观测特性分析 迭代——交会对接目标飞行器调相分析 循环+嵌套迭代——交会对接对远距离导引精度分析

3 STK与Matlab互联——两个强大成熟商业软件的强强联合

4 1 STK与Matlab接口简介—必要性 有效的工程应用例子——交会对接远距离导引精度论证

5 1 STK与Matlab接口简介—必要性 远距离导引精度仿真分析流程图

6 1 STK与Matlab接口简介—必要性 高效（快） 可靠（好） 高效+可靠->技术优势->话语权->地位 应用效果
建模与修改快捷 可靠（好） 辅助实现关键技术环节突破 高效+可靠->技术优势->话语权->地位

7 1 STK与Matlab接口简介—互联与接口函数
- 复制文件“agiCleanPath.m”、“agiInit.m”和“stkInit.m”到指定目录 - 在Matlab中运行指令“stkinit” - 选择M文件的位置 - 选择Mex文件的位置 - 选择默认的其他设置

8 Generic Object Commands
1 STK与Matlab接口简介—互联与接口函数 接口函数 分类 接口函数 说明 STK Server Commands （STK服务器命令） stkDefaultConID 返回缺省的连接ID号 stkDefaultHost 返回缺省的主服务器 stkSetDefaultHost 设置与STK的远程连接 Application Commands （应用命令） stkAutoSave 自动保存 stkHomeDir 返回STK安装目录 stkCurrentDir 返回当前目录 stkSetCurrentDir 设置当前目录 stkObjNames 返回STK对象名称 stkTimerStart 启动计时器 stkTimer 计时器时间 Generic Object Commands （普通对象命令） stkCentralBody 返回坐标系中心体 stkLoadObj 导入对象 stkSaveObj 保存对象 stkNewObj 新建对象 stkUnload 卸载对象 stkUnloadChildren 卸载所有子对象 stkRename 重命名对象 stkTimePeriod 返回场景起止时间 stkTimePeriodInSec 返回场景起止时间，以秒为单位 stkValidObj 返回是否存在有效对象 stkDescrip 获取对象的描述 stkSetDescrip 设置对象的描述

9 1 STK与Matlab接口简介—互联与接口函数
Scenario Commands （场景命令） stkAnimTime 返回当前STK动画时间 stkEpoch 返回STK场景历元 stkSetEpoch 设置 STK 场景 历元 stkSetEpochInSec 改变STK 场景 历元, 以秒为单位 stkSyncEpoch 同步aeroToolbox和STK场景历元 stkConvertDate 转换日期格式 stkSetTimePeriod 设置场景起止时间 stkSetTimePeriodInSec 设置场景起止时间, 以秒为单位 stkScenFilePath 返回STK场景文件目录 stkValidScen 返回是否存在有效场景 Basic Vehicle Data （基本飞行器数据） stkTimePeriod 返回航天器或场景起止时间 stkTimePeriodInSec 返回航天器或场景起止时间, 以秒为单位 设置航天器或场景起止时间 Satellite Propagators （卫星属性） stkSetPropCart 设置笛卡尔坐标形式的卫星运动参数 stkSetPropClassical 设置经典轨道根数形式的卫星运动参数 stkSetPropEqui 设置第二类无奇点根数形式的卫星运动参数 stkSetPropSGP4 设置SGP4形式的卫星运动参数 stkPropagate 卫星轨道积分计算 Vehicle Position and Velocity （飞行器的位置速度） stkPosVelCBF 获取中心体固连坐标系下的位置速度 stkPosVelCBI 获取惯性坐标系下的位置速度 stkEphemerisCBF 获取中心体固连坐标系下的星历 stkEphemerisCBI 获取惯性坐标系下的星历 stkSetEphemerisCBF 设置中心体固连坐标系下的星历 stkSetEphemerisCBI 设置惯性坐标系下的星历 stkLoadEphemeris 从文件中载入星历 stkSaveEphemerisCBF 保存中心体固连坐标系下的星历 stkSaveEphemerisCBI 保存惯性坐标系下的星历 stkSaveDynState 将工作空间中所有数据保存成.mat文件 stkDynState 产生动力学状态数据

10 1 STK与Matlab接口简介—互联与接口函数
Vehicle Attitude （飞行器姿态） stkAttitudeCBF 获取中心体固连坐标系中的姿态四元数 stkAttitudeCBI 获取惯性坐标系中姿态四元数 stkSetAttitudeCBF 设置中心体固连坐标系中的姿态四元数 stkSetAttitudeCBI 设置惯性坐标系中姿态四元数 stkLoadAttitude 从文件中载入姿态数据 stkSaveAttitudeCBF 保存中心体固连坐标系中的姿态 stkSaveAttitudeCBI 保存惯性坐标系中姿态 Basic Vehicle Attitude Types （飞行器基本姿态类型） stkSetAttBasic 设置基本姿态类型 stkSetAttSpinSun 设置自旋轴指向太阳 stkSetAttSpinNadir 设置自旋轴指向天底 stkSetAttYawNadir 设置偏航轴指向天底 stkSetAttSpinning 设置自旋轴为给定方向 stkSetAttFixed 设置惯性系中姿态 stkAttOffset 设置姿态偏移量 Aircraft, GroundVehicle and Ship Propagators （飞机，车辆与船舰属性） stkSetWaypoints 设置路径点 stkLoadWaypoints 载入路径点文件 stkSetGreatArcStart 设置起始时间 Missile Propagators （导弹运动积分） stkSetPropBallistic 设置弹道导弹积分参数 Facilities and Targets （设施与目标） stkPosVelCBF 返回中心体固连坐标系中的位置 stkSetFacPosCBF 设置中心体固连坐标系中设施的位置 stkSetFacPosLLA 设置设施的经度、纬度、高度 Area Targets （面目标） stkGetBoundary 设置面目标边界 STK Tools （STK工具） stkAccess 返回过境时间段 stkReport 生成目标报告 stkAccReport 生成过境报告 stkFindData 从报告中提取数据 Connect Commands （连接指令） stkOpen 打开STK连接 stkClose 关闭STK连接 stkExec 向STK发送指令 stkConnect 连接指令发送函数

11 1 STK与Matlab接口简介—常用接口函数
1）初始化需要的接口函数 stkinit—— 完成STK和Matlab的互联和初始化 conid=stkOpen(stkDefaultHost); —— 返回互联成功的主机端口的连接句柄 2）初始窗口管理 if stkValidScen == 1 stkUnload('/*') end —— 如果已经有打开的场景，则关闭场景。

12 1 STK与Matlab接口简介—常用接口函数
也可以通过交互提示进行窗口管理，如： scen_open = stkValidScen; if scen_open == 1 rtn = questdlg('Close the current scenario?'); if ~strcmp(rtn,'Yes') stkClose(conid) return else stkUnload('/*') end —— 如果已经有打开的场景，则弹出提示对话框，询问是否要关闭现有场景。

13 1 STK与Matlab接口简介—常用接口函数
3）建立场景 stkNewObj('/','Scenario','场景名称'); —— 建立给定名称的场景。 stkSetTimePeriod('10 Apr :00:00.0','12 Apr :00:00.0','GREGUTC'); —— 设置场景的起止时间和采用的时间系统。 stkSetEpoch('10 Apr :00:00.0','GREGUTC'); —— 设置场景的历元。 注：在STK输出的各类报告中，时间起点将以该历元作为时间零点。

14 1 STK与Matlab接口简介—常用接口函数
stkSyncEpoch; —— 同步aeroToolbox和STK场景历元。 注：aeroToolbox函数库包含“安装盘：\Program Files\AGI\STK\Matlab”目录下各种atb开头的内部函数，可以与STK场景历元保持独立。STK场景时间发生改变时，通过调用该语句可保持二者时间同步。 rtn = stkConnect(conid,'Animate','Scenario/场景名称','SetValues "10 Apr :00:00.0" '); rtn = stkConnect(conid,'Animate','Scenario/场景名称','Reset'); —— 设置STK场景动画历元。

15 1 STK与Matlab接口简介—常用接口函数
4）建立航天器 stkNewObj('*/','Satellite','航天器名称'); —— 建立卫星。 最基础最常用的接口函数：stkSetPropClassical 语法： stkSetPropClassical('objPath', 'propagator', 'coordSystem', tStart, tStop, dt, orbitEpoch, semimajorAxis, eccentricity,inclination, argOfPerigee, RAAN, meanAnomaly, coordEpoch)

16 1 STK与Matlab接口简介—常用接口函数
objPath – 有效路径。 Propagator – 选用的动力学模型，可选择为 ‘TwoBody’, ‘J2Perturbation’, ‘J4Perturbation’, ‘HPOP’ or ‘PODS’，分别表示“二体”、“J2摄动”、“J4摄动”、“HPOP模型”或“PODS模型”，区分大小写。 coordSystem – 选用的坐标系，可选择为 ‘Fixed’, ‘J2000’, ‘MeanOfDate’, ‘MeanOfEpoch’, ‘TrueOfDate’, ‘TrueOfEpoch’, ‘B1950’, ‘TEMEOfDate’, ‘TEMEOfEpoch’, ‘AlignmentAtEpoch’，区分大小写。 tStart, tStop – 起止时刻，以相对于历元时刻的秒为单位。 dt – 积分步长，以秒为单位。 orbitEpoch – 轨道历元。

17 1 STK与Matlab接口简介—常用接口函数
semimajorAxis, eccentricity,inclination, argOfPerigee, RAAN, meanAnomaly – 6个轨道根数。 其中，若坐标系选为‘AlignmentAtEpoch’，RAAN为升交点地理经度，若坐标系选为‘J2000’，RAAN为升交点赤经。 coordEpoch – 坐标历元, 所有 '...OfEpoch' 选项均需提供对应的坐标历元。 同样重要的接口函数：stkSetPropCart 语法： stkSetPropCart('objPath', 'propagator', 'coordSystem', tStart, tStop, dt, orbitEpoch, pos, vel, coordEpoch)

18 1 STK与Matlab接口简介—常用接口函数
若动力学模型选为“HPOP模型”，则还需要通过stkConnect函数设置各种力模型参数。 第一类：“HPOP force model” 设置重力、第三体引力、太阳光压、固体潮、海潮、辐射压参数。 第二类：“HPOP drag model” 设置大气密度模型及参数。

19 1 STK与Matlab接口简介—常用接口函数
第一类：“HPOP force model” 语法：HPOP <ObjectPath> Force {ForceModel} <ModelParameters> 其中， < >中的参数表示需要用户根据具体目录安装或参数使用情况设置的， {}中的参数表示STK固有的关键字， 例如： stkConnect(conid,'HPOP','Scenario/场景名/Satellite/航天器名','Force SolarRad Off') <ObjectPath>为 'Scenario/场景名/Satellite/航天器名'， {ForceModel}为 SolarRad， <ModelParameters>为 Off。

20 1 STK与Matlab接口简介—常用接口函数
{ForceModel} 和 <ModelParameters>的可选参数如表所示。 {ForceModel} <ModelParameters> 说明 Gravity "<FilePath>" 或<MaxDegreeValue> 或<MaxOrderValue> <MaxDegreeValue>和<MaxOrderValue>取值为0.0~1000.0，<MaxOrderValue> 必须小于等于 <MaxDegreeValue>，并且都小于所选择的重力场文件的最高阶数。 ThirdBodyGravity <CentralBodyName> {On | Off} {GravSource} [<UserGravValue>] {GravSource}有效值可通过CB, JPL的DE文件或用户自定义设置，用户自定义值需为正值。 SolarRad {SolarRadState} [<Coeff> <AreaMassRatio>] [ShadowModel {ModelType}] [BoundaryMitigation {On | Off}] {SolRadState}可设置为On或Off。 如果设置为 On, 还需要设置-100.0~100.0的光压系数<Coeff>，和 0.0~ 的面质比 参数<AreaMassRatioValue>。阴影模型参数{ModelType}]可取为空、圆柱或双圆锥。 SolidTides {Full | Permanent | Off} [IncTimeDep {On | Off}] [MinAmplitude <Amplitude>] 如果设置为Off则不需要再输入其他参数。 <Amplitude> 在0.0~1.0范围内，单位是米。如果输入Permanent则只考虑与时间不相关的永久固体潮部分，如果输入Full则考虑全部固体潮。 OceanTides {On | Off} [MaxDegree <MaxDegreeValue>] [MaxOrder <MaxOrderValue>] [MinAmplitude <Amplitude>] 如果设置为Off不需要再输入其他参数。 <MaxDegreeValue>为2~30之间的整数， <MaxOrderValue> 为0~30之间的整数， <Amplitude>在0.0~1.0范围内，单位是米。 RadiationPressure  Albedo {On | Off}  Thermal {On | Off}  Coefficient <CoeffValue>  AreaMassRatio <AMRValue>  File "<FilePath>" <CoeffValue> 为 -100~100.0之间的实数。 <AMRValue> 为0.0~9999.0之间的实数，单位是m^2/kg。 "<FilePath>" 地面反射模型文件。

21 1 STK与Matlab接口简介—常用接口函数
stkConnect(conid,'HPOP','Scenario/场景名称/Satellite/航天器名称', 'Force Gravity "C:\Program Files\AGI\STK\STKData\CentralBodies\Earth\JGM3.grv" 8 8'); 'Force SolarRad Off'); 'Force ThirdBodyGravity Sun Off'); 'Force ThirdBodyGravity Moon Off');

22 1 STK与Matlab接口简介—常用接口函数
第二类：“HPOP drag model ” 语法：HPOP <ObjectPath> Drag {DragState} <Coeff> <AreaMassRatio> {"AtmDensModel"} 其中，< >与{}中的参数的含义不变。 {DragState} 设置为On或Off。 设置为Off则不需要输入其他参数。 设置为On时， <Coeff>为-10.0~10.0之间的实数， <AreaMassRatio> 为0.0~9999.0之间的实数， {"AtmDensModel"}的有效输入值和相应的<Parameters>输入值见表所示。

23 1 STK与Matlab接口简介—常用接口函数
{"AtmDensModel"} <Parameters> "Harris-Priester" 附加参数可通过手工或flux/geomag 文件输入： Manual <AverageF10.7> File "<FilePath>" "1976 Standard" N/A "Jacchia 1960" "Jacchia 1971" Manual <DailyF10.7> <AverageF10.7> <GeoMagIndex> "Jacchia 1970" "Jacchia-Roberts" "CIRA 1972" "NRLMSISE 2000" "MSISE 1990" "MSIS 1986" "<UserModelName>" "<FilePath>"

24 1 STK与Matlab接口简介—常用接口函数
例如： stkConnect(conid,'HPOP','Scenario/场景名称/Satellite/航天器名称', 'Drag On "NRLMSISE 2000" Manual '); 'Drag On "NRLMSISE 2000" File "C:\Program Files\AGI\STK\DynamicEarthData\stkFluxGeoMag.fxm"');

25 1 STK与Matlab接口简介—常用接口函数
5）建立地面站 stkNewObj(' */ ','Facility','测站名称') —— 建立地面站。 stkSetFacPosLLA('Scenario/场景名称/Facility/测站名称', [phi; namda; h]) —— 设置地面站的经度、纬度和高度。 stkConnect(conid,'SetConstraint','Scenario/场景名称/Facility/测站名称', 'ElevationAngle Min 最低仰角值') —— 设置地面站的最低仰角。 stkConnect(conid,‘SetConstraint’,‘Scenario/场景名称/Facility/测站名称’, ‘Range Max 最大探测距离值'); —— 设置地面站的最大探测距离。

26 1 STK与Matlab接口简介—常用接口函数
6）输出计算结果 stkPropagate('*/Satellite/航天器名称', t_start, t_stop) [T,r,v, cb] = stkEphemerisCBI('Scenario/场景名称/Satellite/航天器名称', dt, t_start, t_stop) —— 运行航天器并输出给定起止时间内的惯性坐标系位置和速度。 注：输出航天器位置、速度之前必须运行stkPropagate函数，否则力模型参数设置不生效。 interval=stkAccess('Scenario/场景名称/Satellite/航天器名称', 'Scenario/场景名称/Facility/测站名称') —— 输出航天器相对于测站的可见时间段。

27 1 STK与Matlab接口简介—常用接口函数
stkClose(conid) 注：STK和Matlab最多同时可创建2个连接，因此每次程序运行结束通常需要关闭连接，以免下次运行重新建立连接时个数超过限制。 采用接口函数输出的数据上虽然不如手工设置时形式多样，但因为可以在Matlab环境下通过编程进行数据处理，依然可以获得多种需要的计算结果，并且可实现手工设置无法完成的循环计算和多重迭代功能。 上述1）~6）过程为利用STK/Matlab接口进行计算分析的常用或基本函数。

28 2 在复杂航天任务分析中的应用—循环 循环——星座对地观测特性分析
设某一对地观测卫星星座含多颗星，各星之间通过保持一定的几何构型，实现对某一特定区域的重点观测。在系统设计阶段，需要分析对某种给定的构型，星座对某一特定区域进行观测时，各星观测该区域的 时间间隔特性，如该特定区域重访时间间隔平均值、最大值等。 对某一给定点，通过手工设置STK中各类参数可计算和输出星座中各星访问该地点的时间间隔，但对于某一给定区域，需要在该区域内密集采样，采用手工输入输出方式已不再可能，因此需要通过Matlab作为计算引擎调用STK循环计算。 以下给出计算过程和结果。

29 2 在复杂航天任务分析中的应用—循环 1）初始化并关闭所有打开的窗口 % 初始化 stkinit
remMachine = stkDefaultHost; % 关闭所有打开的窗口 delete(get(0,'chiSatren')); conid=stkOpen(remMachine); scen_open = stkValidScen; if scen_open == 1 stkUnload('/*') end

30 2 在复杂航天任务分析中的应用—循环 2) 建立场景、航天器和测站 % 建立场景并设置场景属性
stkNewObj('/','Scenario','RemoteSensor'); stkSetTimePeriod('1 Jan :00:00.0','5 Jan :00:0.0','GREGUTC'); stkSetEpoch('1 Jan :00:00.0','GREGUTC'); stkSyncEpoch; rtn = stkConnect(conid,'Animate','Scenario/RemoteSensor', 'SetValues "1 Jan :00:00.0" '); rtn = stkConnect(conid,'Animate','Scenario/RemoteSensor','Reset'); % 建立航天器和测站 stkNewObj('*/','Satellite','Sat1'); stkNewObj('*/','Satellite','Sat2'); stkNewObj('*/','Satellite','Sat3'); stkNewObj('*/','Satellite','Sat4'); stkNewObj('*/','Facility','Station');

31 2 在复杂航天任务分析中的应用—循环 3) 设置星座中各卫星的属性
t_start=0; t_stop=4*(24*3600); dt=60; orbitEpoch=t_start; a_sat=7000*1000; e_sat= 0; i_sat= 60; w_sat= 0; Raan_sat=0; M_sat=0 ; stkSetPropClassical('*/Satellite/Sat1','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_sat,e_sat,i_sat,w_sat,Raan_sat,M_sat); stkPropagate('*/Satellite/Sat1', t_start, t_stop) stkSetPropClassical('*/Satellite/Sat2','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_sat,e_sat,i_sat,w_sat,Raan_sat+pi/4,M_sat+pi); stkPropagate('*/Satellite/Sat2', t_start, t_stop) stkSetPropClassical('*/Satellite/Sat3','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_sat,e_sat,i_sat,w_sat,Raan_sat+pi/2,M_sat+2*pi); stkPropagate('*/Satellite/Sat3', t_start, t_stop) stkSetPropClassical('*/Satellite/Sat4','J2Perturbation','J2000',t_start,t_stop,dt,orbitEpoch,a_sat,e_sat,i_sat,w_sat,Raan_sat+3*pi/4,M_sat+pi); stkPropagate('*/Satellite/Sat4', t_start, t_stop)

32 2 在复杂航天任务分析中的应用—循环 4）循环计算区域内各点的重访时间间隔及最大值 namda_min=0; namda_max=2*pi;
deta_namda=(namda_max-namda_min)/n_namda; phi_min=-60*pi/180; phi_max=60*pi/180; n_phi=30; deta_phi=(phi_max-phi_min)/n_phi; for i=1:n_namda+1 for j=1:n_phi+1 remain=(n_namda+1)*(n_phi+1)-(i-1)*(n_phi+1)-j namda=namda_min+(i-1)*deta_namda; X(i,j)=namda*180/pi; phi=phi_min+(j-1)*deta_phi; Y(i,j)=phi*180/pi; stkSetFacPosLLA('Scenario/RemoteSensor/Facility/Station', [phi; namda; 0]);

33 2 在复杂航天任务分析中的应用—循环 stkConnect(conid,'SetConstraint','Scenario/RemoteSensor/Facility/Station','ElevationAngle Min '); interval1=stkAccess('Scenario/RemoteSensor/Satellite/Sat1','Scenario/RemoteSensor/Facility/Station'); interval2=stkAccess('Scenario/RemoteSensor/Satellite/Sat2','Scenario/RemoteSensor/Facility/Station'); interval3=stkAccess('Scenario/RemoteSensor/Satellite/Sat3','Scenario/RemoteSensor/Facility/Station'); interval4=stkAccess('Scenario/RemoteSensor/Satellite/Sat4','Scenario/RemoteSensor/Facility/Station');

34 2 在复杂航天任务分析中的应用—循环 if size(interval1)==[0 0] & size(interval2)==[0 0] & size(interval3)==[0 0] & size(interval4)==[0 0] Z(i,j)=20; else temp1=[interval1.start,interval2.start,interval3.start,interval4.start,interval1.stop,interval2.stop,interval3.stop,interval4.stop]; temp2=sort(temp1); Z(i,j)=max(diff(temp2))/3600; end

35 2 在复杂航天任务分析中的应用—循环 5）绘制计算结果并关闭连接 save('mat') set(gcf,'color',[1,1,1])
figure surf(X,Y,Z) xlabel('namda/deg') ylabel('phi/deg') stkClose(conid)

36 2 在复杂航天任务分析中的应用—迭代 迭代——交会对接目标飞行器调相分析
交会对接任务中，目标飞行器在飞船入轨前2~3个月先行入轨，目标飞行器入轨后在一天时间内通过两次变轨进入某一高度的近圆轨道，之后通过大气衰减和主动变轨，使2~3个月后飞船入轨时刻目标飞行器轨道的相位、高度和偏心率满足要求。在任务分析阶段，需要根据给定的空间环境设计目标飞行器的圆化轨道高度，在后续2~3个月的时间内尽可能通过大气衰减作用满足各种终端约束条件，减少主动变轨的能量消耗。 第13圈变轨 第4圈变轨

37 2 在复杂航天任务分析中的应用—迭代 为了重点说明STK/Matlab接口的使用，对上述设计进行简化，假定不进行调整相位、高度和偏心率的变轨，设计过程简化为迭代搜索满足飞船入轨时刻相位约束条件的目标飞行器圆化轨道的高度，为单层迭代过程。 需要补充说明的是，在STK中尽管Astrogator模式可实现一定功能的迭代，但在Astrogator模式下空间环境参数无法根据需要进行设置，并且难以设定某些特殊的约束条件，这种情况下可以利用Matlab的编程能力通过STK/Matlab接口实现迭代求解。 以下给出计算过程。

38 2 在复杂航天任务分析中的应用—迭代 1）初始化和并关闭所有打开的窗口 与前述过程相同，不再重复。 2）建立场景和航天器
目标飞行器入轨后，在第4圈近地点变轨抬高远地点高度，在第13圈远地点变轨进入圆轨道，之后进入大气衰减自由飞行阶段。以下建立场景后，再建立3个航天器分别用于不同阶段的轨道预报。 stkNewObj('/','Scenario','GM'); stkSetTimePeriod(‘6 Sep :11:53', 6 Dec :11:53 ','GREGUTC'); stkSetEpoch(6 Sep :11:53 ','GREGUTC'); stkSyncEpoch; rtn = stkConnect(conid,'Animate','Scenario/GM','SetValues " 6 Sep :11:53 " '); rtn = stkConnect(conid,'Animate','Scenario/GM','Reset'); stkNewObj('*/','Satellite','MB_1'); stkNewObj('*/','Satellite','MB_2'); stkNewObj('*/','Satellite','MB_3');

39 2 在复杂航天任务分析中的应用—迭代 3）外推至第4圈近地点 % 设置J2000坐标系中初始轨道根数
t_start_1=0; t_stop_1=0.187*(24*3600); dt=60; orbitEpoch=t_start_1; coordEpoch=t_start_1; a_MB= ; e_MB= ; i_MB=42.747; w_MB= ; Raan_MB= 7.309; M_MB= 0 ;

40 2 在复杂航天任务分析中的应用—迭代 % 选用HPOP模型，设置大气密度模型、重力场模型、光压和第三体摄动模型。
stkSetPropClassical('*/Satellite/MB_1','HPOP','AlignmentAtEpoch',t_start_1,t_stop_1,dt,orbitEpoch,a_MB,e_MB,i_MB,w_MB,Raan_MB,M_MB,coordEpoch); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_1','Drag On "NRLMSISE 2000" Manual '); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_1','Force Gravity "C:\Program Files\AGI\STK\STKData\CentralBodies\Earth\JGM3.grv" 8 8'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_1','Force SolarRad Off'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_1','Force ThirdBodyGravity Sun Off'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_1','Force ThirdBodyGravity Moon Off'); stkPropagate('*/Satellite/MB_1', t_start_1, t_stop_1) [T_MB_1, r_MB_1, v_MB_1, cb] = stkEphemerisCBI('Scenario/GM/Satellite/MB_1',dt, t_start_1, t_stop_1);

41 2 在复杂航天任务分析中的应用—迭代 4）迭代求解第4圈近地点变轨速度增量（或求解圆化轨道高度）
% 采用牛顿迭代法求解，约束条件为第89天目标飞行器轨道面与飞船入轨点共面时刻的相位在114°±10 °范围内。 deta_v_2=9; error=100; while abs(error)>1 Cal_remain; index_RG=find(T_interp>89); error=beta_interp(index_RG(1))-114 deta_v_2=deta_v_2+error/3600; End 其中，Cal_ remain为需要循环计算的用户自定义子函数，用于计算给定一个第4圈近地点变轨速度增量后，从第4圈近地点变轨外推至第13圈远地点并进行轨道圆化变轨后的轨道，轨道外推一直到第89天目标飞行器轨道面与飞船入轨点共面时刻。输出结果外共面时刻的相位。

42 2 在复杂航天任务分析中的应用—迭代 Cal_ remain子函数为： % 第4圈变轨外推至第13圈变轨
pos_2=r_MB_1(:,length(r_MB_1)); vel_2=v_MB_1(:,length(v_MB_1)); ix=vel_2/norm(vel_2); ir=pos_2/norm(pos_2); iz=CrossProduct(ir,ix)/norm(CrossProduct(ir,ix)); iy=CrossProduct(ix,iz); thita_2=0; vel_2=vel_2+deta_v_2*(cos(thita_2)*ix+sin(thita_2)*iy); t_start_2=t_stop_1; t_stop_2=0.782*(24*3600); dt=60; orbitEpoch=t_start_2; stkSetPropCart('*/Satellite/MB_2','HPOP','J2000',t_start_2, t_stop_2, dt, orbitEpoch, pos_2, vel_2) stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_2','Drag On "NRLMSISE 2000" Manual '); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_2','Force Gravity "C:\Program Files\AGI\STK\STKData\CentralBodies\Earth\JGM3.grv" 8 8');

43 2 在复杂航天任务分析中的应用—迭代 stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_2','Force SolarRad Off'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_2','Force ThirdBodyGravity Sun Off'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_2','Force ThirdBodyGravity Moon Off'); stkPropagate('*/Satellite/MB_2', t_start_2, t_stop_2) [T_MB_2, r_MB_2, v_MB_2, cb] = stkEphemerisCBI('Scenario/GM/Satellite/MB_2',dt, t_start_2, t_stop_2); %第13圈变轨外推至飞船入轨时刻 pos_3=r_MB_2(:,length(r_MB_2)); vel_3=v_MB_2(:,length(v_MB_2)); u_e= e+014; ix=vel_3/norm(vel_3); ir=pos_3/norm(pos_3); iz=CrossProduct(ir,ix)/norm(CrossProduct(ir,ix)); iy=CrossProduct(ix,iz); it=CrossProduct(iz,ir); vel_3=sqrt(u_e/norm(pos_3))*it;

44 2 在复杂航天任务分析中的应用—迭代 t_start_3=T_MB_2(length(T_MB_2));
t_stop_3=92*(24*3600); dt=60; orbitEpoch=t_start_3; stkSetPropCart('*/Satellite/MB_3','HPOP','J2000',t_start_3, t_stop_3, dt, orbitEpoch, pos_3, vel_3) stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_3','Drag On "NRLMSISE 2000" Manual '); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_3','Force Gravity "C:\Program Files\AGI\STK\STKData\CentralBodies\Earth\JGM3.grv" 8 8'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_3','Force SolarRad Off'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_3','Force ThirdBodyGravity Sun Off'); stkConnect(conid,'HPOP','Scenario/GM/Satellite/MB_3','Force ThirdBodyGravity Moon Off'); stkPropagate('*/Satellite/MB_3', t_start_3, t_stop_3) [T_MB_3, r_MB_3, v_MB_3, cb] = stkEphemerisCBI('Scenario/GM/Satellite/MB_3',dt, t_start_3, t_stop_3);

45 2 在复杂航天任务分析中的应用—迭代 %计算共面时刻相位
T_MB_all=[T_MB_1,T_MB_2(2:length(T_MB_2)),T_MB_3(2:length(T_MB_3))]'/(24*3600); R_MB_all=[r_MB_1,r_MB_2(:,2:length(r_MB_2)),r_MB_3(:,2:length(r_MB_3))]'; V_MB_all=[v_MB_1,v_MB_2(:,2:length(v_MB_2)),v_MB_3(:,2:length(v_MB_3))]'; data_RGD=load('data_pos_RGD.m'); % 读入入轨点在惯性系中的位置。 T_RGD_all=data_RGD(:,1); r_RGD_all=data_RGD(:,2:4); step_window=1/60/24; t_start=0; t_end=92; T=[t_start:step_window:t_end]'; index1=find(T_MB_all>t_start & T_MB_all<t_end); R_MB=[interp1(T_MB_all(index1),R_MB_all(index1,1),T),interp1(T_MB_all(index1),R_MB_all(index1,2),T),interp1(T_MB_all(index1),R_MB_all(index1,3),T)]; V_MB=[interp1(T_MB_all(index1),V_MB_all(index1,1),T),interp1(T_MB_all(index1),V_MB_all(index1,2),T),interp1(T_MB_all(index1),V_MB_all(index1,3),T)]; index2=find(T_RGD_all>t_start & T_RGD_all<t_end); R_RGD=[interp1(T_RGD_all(index2),r_RGD_all(index2,1),T),interp1(T_RGD_all(index2),r_RGD_all(index2,2),T),interp1(T_RGD_all(index2),r_RGD_all(index2,3),T)];

46 2 在复杂航天任务分析中的应用—迭代 j=1; for i=1:length(T)
AngMomentum_MB=CrossProduct(R_MB(i,:),V_MB(i,:))/norm(CrossProduct(R_MB(i,:),V_MB(i,:))); alpha(i)=acos(Dotproduct(R_RGD(i,:)/norm(R_RGD(i,:)),AngMomentum_MB'))*180/pi; beta(i)=acos(Dotproduct(R_RGD(i,:)/norm(R_RGD(i,:)),R_MB(i,:)/norm(R_MB(i,:))))*180/pi; if DotProduct(CrossProduct(R_RGD(i,:),R_MB(i,:)),AngMomentum_MB)<=0 beta(i)=-beta(i); end if i>1 & alpha(i)<90 & alpha(i-1)>90 index(j)=i; T_interp(j)=interp1([alpha(i),alpha(i-1)],[T(i),T(i-1)],90); beta_interp(j)=interp1([T(i),T(i-1)],[beta(i),beta(i-1)],T_interp(j)); j=j+1;

47 2 在复杂航天任务分析中的应用—迭代 5）绘制计算结果并关闭连接 %迭代收敛后，将迭代结果对应的轨道主要特征绘制出来。
set(gcf,'color',[1,1,1]) subplot(2,1,1) hold on; plot(T_interp,beta_interp,'k.',T_interp,beta_interp,'k'); plot([0 t_end],[104,104],'r',[0 t_end],[114,114],'g',[0 t_end],[124,124],'r') axis([0,t_end,-180,180]) subplot(2,1,2) h=(sqrt(R_MB_all(:,1).^2+R_MB_all(:,2).^2+R_MB_all(:,3).^2) )/1000; plot(T_MB_all,h,'k'); plot([0 t_end],[332,332],'g','LineWidth',2) axis([0,t_end,300,400]) stkClose(conid)

48 2 在复杂航天任务分析中的应用—迭代 飞船入轨时刻 相位取值范围 目标飞行器相位 需要的终点轨道高度 初始圆轨道高度

49 2 在复杂航天任务分析中的应用—循环+嵌套迭代
循环+嵌套迭代——交会对接对远距离导引精度分析 入轨 第1次变轨 （抬高近地点） 第2次变轨 （修正轨道面） 第3次变轨 （抬高远地点） 第4次变轨 第5次变轨 （组合修正） 导引终点 目标 飞船

50 2 在复杂航天任务分析中的应用—循环+嵌套迭代
远距离导引的终端约束为：在给定时刻，飞船与目标飞行器共面，处于比目标飞行器轨道略低的近圆轨道上，且相对目标飞行器相位满足要求。上述约束条件的数学描述即为给定时刻飞船相对于目标飞行器的相对位置和相对速度满足给定要求（标称值及误差范围）。 可用于规划的变量即为各次变轨的时刻及变轨速度增量。为减少可变因素并简化工程实施，变轨时刻通常选择在近地点或远地点，变轨速度增量限定为水平方向（轨道面修正除外）。由于各次变轨均会对终点状态产生影响，因此每次变轨规划都是后续剩余变轨规划变量的统一综合规划。

51 2 在复杂航天任务分析中的应用—循环+嵌套迭代
远距离导引精度分析具体计算过程可分为两个层次的循环或迭代： 最外层循环为误差样本空间中循环多次采样进行仿真计算，对于每一次采样，在得到一组各误差源的具体采样误差值后，加到相应环节进行仿真计算，多次采样计算结束后进行统计分析； 向内第二层循环为各次变轨控制量的迭代计算，属于多变量规划问题，规划的目标量为远距离导引终点飞船相对于目标飞行器的相对位置和速度状态，规划变量的数量和内容在各次变轨规划时各不相同，第一次变轨前规划时，可用于规划的变量有6个，规划目标6个，可得到唯一解；第二次规划时可用的规划变量有7个，规划目标6个，存在两种规划方案；第三次和第四次规划变量和规划目标均为4个，同样是唯一解；第五次规划变量只有2个，规划目标4个，需要通过加权对规划目标进行取舍。通常情况下为了避免嵌套迭代减少循环层次，采用各规划变量同时规划的方式，但对于第二次规划，由于轨道面内控制与轨道面法向控制不相关，只能先规划轨道面内的控制，再规划轨道面法向控制。

52 2 在复杂航天任务分析中的应用—循环+嵌套迭代

53 2 在复杂航天任务分析中的应用—循环+嵌套迭代
1）初始化并关闭所有打开的窗口 同前。 2）建立场景和航天器 stkNewObj('/','Scenario','LongRangeGuidance'); stkSetTimePeriod('10 Apr :00:00.0','12 Apr :00:00.0','GREGUTC'); stkSetEpoch('10 Apr :00:00.0','GREGUTC'); stkSyncEpoch; rtn = stkConnect(conid,'Animate','Scenario/LongRangeGuidance','SetValues "10 Apr :00:00.0" '); rtn = stkConnect(conid,'Animate','Scenario/LongRangeGuidance','Reset'); stkNewObj('*/','Satellite','sat'); 注：因远距离导引过程分段过多，不宜每段定义一个航天器，而是采用定义一个航天器进行轨道计算，再将计算结果针对不同分段用不同变量名保存的方式。

54 2 在复杂航天任务分析中的应用—循环+嵌套迭代
3）设置航天器基本参数 u_e= e+014; a_temp=xx; e_temp= xx; i_temp= xx; w_temp= xx; Raan_temp= xx; M_temp= xx; t_start_temp=0; t_stop_temp=60; orbitEpoch_temp=t_start_temp; dt=60; stkSetPropClassical('*/Satellite/sat','HPOP','J2000',t_start_temp,t_start_temp+dt,dt,orbitEpoch_temp,a_temp,e_temp,i_temp,w_temp,Raan_temp,M_temp); stkConnect(conid,'HPOP','Scenario/LongRangeGuidance/Satellite/sat','Drag On "NRLMSISE 2000" Manual '); stkConnect(conid,'HPOP','Scenario/LongRangeGuidance/Satellite/sat','Force Gravity "C:\Program Files\AGI\STK\STKData\CentralBodies\Earth\JGM3.grv" 8 8'); stkConnect(conid,'HPOP','Scenario/LongRangeGuidance/Satellite/sat','Force SolarRad Off'); stkConnect(conid,'HPOP','Scenario/LongRangeGuidance/Satellite/sat','Force ThirdBodyGravity Sun Off'); stkConnect(conid,'HPOP','Scenario/LongRangeGuidance/Satellite/sat','Force ThirdBodyGravity Moon Off'); stkPropagate('*/Satellite/sat', t_start_temp, t_stop_temp);

55 2 在复杂航天任务分析中的应用—循环+嵌套迭代
4）最外层循环：误差样本空间大样本数采样和结果统计 m_max=3000; for mm=1:m_max Trajectory_Strategy_1 Trajectory_Strategy_2 Trajectory_Strategy_3 Trajectory_Strategy_4 Trajectory_Strategy_5 Trajectory_real Cal_output Draw save mat_USB; End 其中，Trajectory_Strategy_1、Trajectory_Strategy_2… Trajectory_Strategy_5为远距离导引段5次轨控策略计算的自定义子函数，Trajectory_real为实际轨道计算的自定义子函数。每次循环各子程序自动生成各误差源随机误差。

56 谢 谢