STK与matlab交互 Astrogator模块（14）

一、背景介绍

高轨卫星的轨道保持。与任何其它轨道状态一样，地球同步轨道也会受到各种扰动力的影响，这些摄动力会影响GEO卫星在位置方面的稳定性。摄动的主要来源是地球的非地球位势、太阳辐射压力和第三体效应（主要是月球和太阳）。总体影响是经度漂移（其方向取决于经度本身）和倾斜度随时间的变化。

GEO卫星运营商将每个航天器放在一个定义明确的盒子里，以免干扰周围空间中的其它卫星。这项活动被称为“站岗”。在这里，我们将了解如何使用自动序列设计和执行台位保持。我们将创建一个东西向测站保持序列（基于经度控制）和一个南北向测站维持序列。（纯倾斜度控制）

二、场景初始化

首先根据OrbitWizard创建一个参考星Reference，将目标星控制在参考星的Box框架内，将其定位在13°的位置，代码如下所示

%%  本程序用于构建新场景
close all
clear
clc
%% 分析要求1 
%%  *******************  以下为联动 STK 代码
%  ------------  利用Com库方式调用STK引擎
% Get reference to running STK instance
uiApplication = actxGetRunningServer('STK12.application');
% Get our IAgStkObjectRoot interface
root = uiApplication.Personality2;

checkempty = root.Children.Count;
if checkempty ~= 0
    root.CurrentScenario.Unload
    root.CloseScenario;
    
end

%  ---------- 创建新的场景
root.NewScenario('GEOStationKeeping');
StartTime='1 Apr 2016 00:00:00.000';
StopTime='18 Jul 2017 00:00:00.000';
root.ExecuteCommand(['SetAnalysisTimePeriod * "',StartTime,'" "',StopTime,'"']);
root.ExecuteCommand('Animate * Reset');
% root.LoadScenario('D:\STK 11 (x64)\SatelliteDamage\AreaCoverage.sc'); 
% 
Sat_Name='Reference';
satellite=root.CurrentScenario.Children.New('eSatellite',Sat_Name);
root.ExecuteCommand(['OrbitWizard */Satellite/',Sat_Name,' Geosynchronous Inclination 0 SubsatellitePoint 13.0 Color green']);
% 使地固系显现出来
satellite.VO.OrbitSystems.FixedByWindow.IsVisible=1;
satellite.VO.OrbitSystems.InertialByWindow.IsVisible=0;
satellite.VO.Proximity.GeoBox.IsVisible=1;
satellite.VO.Proximity.GeoBox.Longitude=13;
satellite.VO.Proximity.GeoBox.NorthSouth=0.5;
satellite.VO.Proximity.GeoBox.EastWest=0.5;
satellite.VO.Proximity.GeoBox.Radius=42166.3;

satellite.Graphics.Attributes.Inherit=0;
satellite.Graphics.Attributes.LabelVisible=0;
satellite.Graphics.Attributes.IsOrbitVisible=0;

添加一颗GEO_Sat，使用Astrogator模块进行轨道计算，代码如下

%% 插入一颗机动星
Sat_Name2='GEO_Sat';
satellite2=root.CurrentScenario.Children.New('eSatellite',Sat_Name2);
satellite2.SetPropagatorType('ePropagatorAstrogator'); 
satellite2.Propagator;
sma=42166.3;
Ecc=0;
Inc=0;
w=0;
LonAsc=13;
TA=352.96;
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList Initial_State Propagate']);
InitialState=satellite2.Propagator.MainSequence.Item(0);
%% 初始化卫星参数
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.CoordinateType Modified Keplerian']);
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.InitialState.Epoch ',StartTime,' UTCG']);
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.sma ',num2str(sma),' km']);
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.ecc ',num2str(Ecc)]);
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.inc ',num2str(Inc),' deg']);
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.w ',num2str(w),' deg']);

root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.LAN ',num2str(LonAsc),' deg']);
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.TA ',num2str(TA),' deg']);

%% 
Propagate=satellite2.Propagator.MainSequence.Item(1);
Propagate.StoppingConditions.Item(0).Properties.Trip=86400*365.25;
Propagate.StoppingConditions.Item(0).Properties.Inherited=0;

satellite2.VO.OrbitSystems.Add('Satellite/Reference VVLH System')
satellite2.VO.OrbitSystems.InertialByWindow.IsVisible=0;
satellite2.VO.Pass.TrackData.PassData.Orbit.SetLeadDataType('eDataNone');
satellite2.VO.Pass.TrackData.PassData.Orbit.SetTrailDataType('eDataTime');
satellite2.VO.Pass.TrackData.PassData.Orbit.TrailData.Time=720*3600;
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' RunMCS']);

观察可以发现，由于精确轨道预报模型，卫星一直处于东漂状态，我们需要将卫星保持在站点保持的边界内，由于东漂离开边界，我们将首先实施东西方保持程序。

2.1 添加东西位置保持部分

在本文中需要利用Utilities的菜单，选择Component Browser的选项，选择Calculation Objects的目录，选择Keplerian Elems的目录。

选择Inclination，点击Duplicate，重新命名为TOD Inclination。选择坐标系为Earth TOD

其操作对应的MATLAB代码为

root.ExecuteCommand(['ComponentBrowser */ Duplicate "Calculation Objects" Inclination "TOD"']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Calculation Objects" "TOD" CoordSystem "CentralBody/Earth TOD"']);

添加一个East-West位置保持的程序自动序列

然后定义这个AutoSequence，首先插入一个Sequence，把它重命名为EW_Station_Keeping。在Sequence下面添加一个Propagate，把它命名为Apogee，把Propagate颜色设置为红色，把其停止情况设置为Apoapsis，取消Duration

satellite2.Propagator.AutoSequence.Add('EW Station Keeping');

satellite2.Propagator.AutoSequence.Item(1).Sequence.Insert('eVASegmentTypeSequence','EW_Station_Keeping','-');
EW=satellite2.Propagator.AutoSequence.Item(1).Sequence.Item(0).Segment;
EW.Insert('eVASegmentTypePropagate','Stop_on_Apogee','-');
EW.Item(0).Properties.Color=255;
stopcondition=EW.Item(0).StoppingConditions;
stopcondition.Add('Apoapsis');
stopcondition.Remove('Duration');

添加一个目标序列，命名为Target Turn Around，在该目标序列下面，添加一个机动序列，命名为EW_Burn，一个预报序列命名为Prop to Second Node，一个Return序列，种类为Enable(except Profiles bypass)，一个预报序列，命名为Prop to Ascending Node

当在近地点时，EW Burn部分施加Delta-V整体的大小是一个独立的变量。在这个之和，我们将一直运行直到第二个点到来。这是因为我们将在检查结果前等一会儿（记住，在这个时候，卫星正在东漂移）。拥有了Return部分，把它设置为bypass，序列Prop to Ascending Node只有在微分修正器收敛的时候运行。

EW.Insert('eVASegmentTypeTargetSequence','Target_Turn_Around','-');
% 
Target=EW.Item(1).Segment;
Target.Insert('eVASegmentTypeManeuver','EW_Burn','-');
Target.Insert('eVASegmentTypePropagate','Prop_to_Second_Node','-');
Prop_to_Second_Node=Target.Item(1);
Prop_to_Second_Node.Properties.Color=42495;
Target.Insert('eVASegmentTypeReturn','return','-');
Return=Target.Item(2);
Return.ReturnControlToParentSequence='eVAReturnControlEnableExceptProfilesBypass';
Target.Insert('eVASegmentTypePropagate','Prop_to_Asending_Node','-');
Prop_to_Asending_Node=Target.Item(3);
Prop_to_Asending_Node.Properties.Color=65535;

修改EW Burn 将其转换为Attitude Control to Thrust Vector，将施加坐标转换为VNC(Earth)，选择X（Velocity）作为控制变量。

Maneuver=Target.Item(0);
Maneuver.Maneuver.SetAttitudeControlType('eVAAttitudeControlThrustVector');
Maneuver.Maneuver.AttitudeControl.ThrustAxesName='Satellite VNC(Earth)';
Maneuver.EnableControlParameter('eVAControlManeuverImpulsiveCartesianX');

定义Prop to Second Node部分。添加一个Ascending Node停止情况。将重复次数设置为2。它将运行直到第二个点到达。选择参考坐标系为Earth TOD。移除Duration情况

Prop_to_Second_Node.StoppingConditions.Add('AscendingNode');
Prop_to_Second_Node.StoppingConditions.Remove('Duration');
Prop_to_Second_Node.StoppingConditions.Item(0).Properties.RepeatCount=2;
Prop_to_Second_Node.StoppingConditions.Item(0).Properties.CoordSystem='CentralBody/Earth TOD';

定义Prop to Ascending Node部分，添加一个Ascending Node停止情况，将坐标设置为TOD，添加约束条件，选择UserDefined选项，将其重新命名为Min_Longitude，将Criteria的值设置为Greater Than Minimum。将CalcObject设置为Geodetc>Longitude，将Tolerance设置为0.05；

Prop_to_Asending_Node.StoppingConditions.Add('AscendingNode');
Prop_to_Asending_Node.StoppingConditions.Item(1).Properties.CoordSystem='CentralBody/Earth TOD';
Prop_to_Asending_Node.StoppingConditions.Item(1).Properties.Constraints.Add('UserDefined');
Prop_to_Asending_Node.StoppingConditions.Item(1).Properties.Constraints.Item(0).Name='Min_Longitude';
Prop_to_Asending_Node.StoppingConditions.Item(1).Properties.Constraints.Item(0).Criteria='eVACriteriaGreaterThanMinimum';
Prop_to_Asending_Node.StoppingConditions.Item(1).Properties.Constraints.Item(0).CalcObjectName='Longitude';
Prop_to_Asending_Node.StoppingConditions.Item(1).Properties.Constraints.Item(0).Tolerance=0.05;
Prop_to_Asending_Node.StoppingConditions.Remove('Duration');

选择结果，选择Prop to Ascending Node部分。添加Math作为独立的变量，改变ComponentName名字为Minimum_Longitude，将计算对象设置为Longitude。

Prop_to_Asending_Node.Results.Add('Minimum_Value');
Prop_to_Asending_Node.Results.Item(0).Name='Minimum_longitude';
Prop_to_Asending_Node.Results.Item(0).CalcObjectName='Longitude';

选择目标序列，打开Differential Corrector。将ImpulsiveMnver.Cartesian.X的组件作为控制参数。选择Pertubation为0.01m/s，MaxStep为1m/s。这个卫星将会进行小的调整来保持轨道，我们想要我们最大步长和摄动在一个小的规模下。将Minimum_Longitude作为结果，将期望值调整到12.55°。这颗卫星在它机动期间会保持在盒子里，我们允许保护边界来设置我们期望的值到12.55°，将Tolerance设置为0.001deg，将Action field设置为Run active profiles

Differential=EW.Item(1).Profiles.Item(0);
ControlParameters=Differential.ControlParameters;
Results=Differential.Results;
Vx=ControlParameters.Item(0);
Vx.Enable=1;
Vx.Perturbation=0.01;
Vx.MaxStep=1;
Min_lon=Results.Item(0);
Min_lon.Enable=1;
Min_lon.DesiredValue=12.55;
Min_lon.Tolerance=0.001;
EW.Item(1).Action='eVATargetSeqActionRunActiveProfiles';

2.2 添加南北位置保持部分

参照Automatic Sequence Browse命名为NS Station Keeping，双击NS Station Keeping来编辑，添加的序列如下图所示

satellite2.Propagator.AutoSequence.Add('NS Station Keeping');
satellite2.Propagator.AutoSequence.Item(2).Sequence.Insert('eVASegmentTypeTargetSequence','NS_SK_Drift','-');
NS=satellite2.Propagator.AutoSequence.Item(2).Sequence.Item(0).Segment;
NS.Insert('eVASegmentTypeManeuver','NS_Burn','-');

将姿态控制设置为Thrust Vector，设置控制系为VNC(Earth)，选择X和Y作为的控制变量，设置Y为-200m/s。定义结果为NS Burn部分，选择结果为Math-Difference作为变量，选择ComponentName为SMA_Diff。改变CalcObject为Semimajor Axis。设置DifferenceOrder为CurrentMinusInitial。将Elems>TOD Inclination为second component variable.

Maneuver2=NS.Item(0);
Maneuver2.Maneuver.SetAttitudeControlType('eVAAttitudeControlThrustVector');
Maneuver2.EnableControlParameter('eVAControlManeuverImpulsiveCartesianX');
Maneuver2.EnableControlParameter('eVAControlManeuverImpulsiveCartesianY');
Maneuver2.Maneuver.AttitudeControl.Y=-200;
% 定义边界条件
Maneuver2.Results.Add('Difference');
Maneuver2.Results.Item(0).Name='SMA_Diff';
Maneuver2.Results.Item(0).DifferenceOrder='eVADifferenceOrderCurrentMinusInitial';
Maneuver2.Results.Item(0).CalcObjectName='Semimajor Axis';
Maneuver2.Results.Add('TOD');

定义Thrust Component

选择NS SK Drift。双击Differential Corrector，将ImpulsiveMnvr.Cartesian.X,Y。将X，Y component，摄动设置为0.1m/s，最大步数为10m/s。同时设置好SMA_Diff的期望值为0，允许成都为0.01km。将生成模式改为By tolerance，将模式改为Run Active Profiles

Differential2=satellite2.Propagator.AutoSequence.Item(2).Sequence.Item(0).Profiles.Item(0);
ControlParameters2=Differential2.ControlParameters;
Results2=Differential2.Results;
Vx=ControlParameters2.Item(0);
Vx.Enable=1;
Vx.Perturbation=0.1;
Vx.MaxStep=10;
Vy=ControlParameters2.Item(1);
Vy.Enable=1;
Vy.Perturbation=0.1;
Vy.MaxStep=10;
SMA_Diff=Results2.Item(0);
SMA_Diff.Enable=1;
SMA_Diff.DesiredValue=0;
SMA_Diff.Tolerance=0.01;
SMA_Diff.ScalingMethod='eVADCScalingMethodTolerance';
TOD=Results2.Item(1);
TOD.Enable=1;
TOD.DesiredValue=0;
TOD.Tolerance=0.001;
TOD.ScalingMethod='eVADCScalingMethodTolerance';

satellite2.Propagator.AutoSequence.Item(2).Sequence.Item(0).Action='eVATargetSeqActionRunActiveProfiles';

2.3创造额外的停止时间

在定义两个站位保持的自动序列后，我们现在专注于主序列部分。我们将添加新的停止条件来限制卫星如何工作。在同时两件事情发生。卫星接近了GEO盒子到达东边界，EW StationKeeping机动被需求。卫星接近了GEO盒子到达南边界。一个NS StationKeeping机动被需求

创建一个额外的停止条件。把GEO_Sat's卫星的Orbit，选择Propagate部分，添加Apoapsis停止条件，把它命名为Apoapsis EW。在Sequence部分EW Station Keeping，将重复次数设置为3次，选择最大预报时间10yr。选择OK关闭。

satellite2.Propagator.MainSequence.Remove('Propagate');
satellite2.Propagator.MainSequence.Insert('eVASegmentTypePropagate','Apoapsis_EW','-');
Apoapsis_EW=satellite2.Propagator.MainSequence.Item(1);

Apoapsis_EW.StoppingConditions.Add('Apoapsis');
Apoapsis_EW.StoppingConditions.Remove('Duration');
Apoapsis_EW.StoppingConditions.Item(0).Properties.RepeatCount=3;
Apoapsis_EW.StoppingConditions.Item(0).Properties.Sequence='EW Station Keeping';
Apoapsis_EW.MaxPropagationTime=10*365.25*86400;

建立一个前提条件在东西站点保持的情况，我们将添加一个情况来限制航天器沿着经度的距离，点击该条件的before条件，选择longitude作为停止情况，设置值为13.5deg。

Before=Apoapsis_EW.StoppingConditions.Item(0).Properties.BeforeConditions;
Before.Add('Longitude');
Before.Item(0).Properties.Trip=13.5;

添加Ascending Node NS Stopping Condition。接下来添加一个停止情况触发North South StationKeeping 自由序列。添加一个AscendingNode停止情况，命名为AscendingNode NS，选择序列为NS Station Keeping。将坐标系设置为Earth TOD。设置最大时间为10yr。

Apoapsis_EW.StoppingConditions.Add('AscendingNode');
Apoapsis_EW.StoppingConditions.Item(1).Properties.Sequence='NS Station Keeping';
Apoapsis_EW.StoppingConditions.Item(1).Properties.CoordSystem='CentralBody/Earth TOD';

添加一个约束条件为NS Station Keeping。我们将添加一个约束条件来限制航天器的轨道倾角。选择该约束条件为UserDefined，双击ComponentName，改名为MaxInclination。选择Criteria为Great than。将CalcObject换成Keplerian TOD Inclination，将坐标系设置为Earth TOD坐标系，最后设置Values为0.45°

Constraint1=Apoapsis_EW.StoppingConditions.Item(1).Properties.Constraint;
Constraint1.Add('UserDefined');
Constraint1.Item(0).Name='MaxInclination';
Constraint1.Item(0).Criteria='eVACriteriaGreaterThan';
MaxInclination=Constraint1.Item(0);
MaxInclination.CalcObjectName='TOD';
MaxInclination.Value=0.45;

这样确保了运行会在倾角到达0.45°之前停止。

创建一个最大的停止情况。因为最大条件不会运行450天。代码如下

satellite2.Propagator.MainSequence.Item(1).StoppingConditions.Add('AlwaysTripped');
MaxDuration=satellite2.Propagator.MainSequence.Item(1).StoppingConditions.Item(2);
MaxDuration.Name='MaxDuration';
MaxDuration.Properties.Constraints.Add('UserDefined');
MaxDuration.Properties.Constraints.Item(0).Name='MaxDuration';
MaxDuration.Properties.Constraints.Item(0).Criteria='eVACriteriaGreaterThan';
MaxDuration.Properties.Constraints.Item(0).CalcObjectName='Duration';
MaxDuration.Properties.Constraints.Item(0).Value=450*86400;

最后 Clear Graphic 点击运行MCS，再点击Clear Graphic来移除画好的轨迹

root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' RunMCS']);
root.ExecuteCommand(['Astrogator */Satellite/',Sat_Name2,' ClearDWCGraphics');