STK与matlab交互 Astrogator模块（3）

矢志航天的阿洪

已于 2024-03-15 10:37:42 修改

阅读量1.2k

点赞数 26

文章标签： matlab 交互

于 2024-03-15 10:34:59 首次发布

本文链接：https://blog.csdn.net/weixin_57997461/article/details/136731841

版权

一、问题介绍

结合（2）介绍的STK中自带的LAMBERT求解器，本文主要介绍matlab与STK交互，自动化完成上述任务需求。

二、操作步骤

（1）建立场景和初始化卫星。（文章1详细介绍，这里只给出相应的代码）得到的场景如下图所示

%% 利用STK自带的LAMBERT库来实现目标拦截
uiApplication = actxGetRunningServer('STK11.application');
% Get our IAgStkObjectRoot interface
root = uiApplication.Personality2;
checkempty = root.Children.Count;
if checkempty ~= 0
    root.CurrentScenario.Unload
    root.CloseScenario;
end
%% 根据你的需要设定场景的名称
root.NewScenario('optimal_control');
StartTime = '26 Jan 2024 04:00:00.000';    % 场景开始时间
StopTime = '10 Feb 2024 04:00:00.000';     % 场景结束时间
root.ExecuteCommand(['SetAnalysisTimePeriod * "',StartTime,'" "',StopTime,'"']);
root.ExecuteCommand(' Animate * Reset');
Sat_name1='blue';
satellite1= root.CurrentScenario.Children.New('eSatellite', Sat_name1);
satellite1.SetPropagatorType('ePropagatorAstrogator'); 
satellite1.Propagator;
%% 生成蓝方星和拦截红方星的初始轨道数据
e=0;Sat_TA=0;Sat_LAN=165;
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList Initial_State Propagate']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.CoordinateType Modified Keplerian']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Epoch ',StartTime,' UTCG']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.ElementType "Kozai-Izsak Mean"']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.Period 86169.6 sec']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.ecc ',num2str(e)]);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.inc 0 deg']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.w 0 deg']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.TA ',num2str(Sat_TA),' deg']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.LAN ',num2str(Sat_LAN),' deg']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Propagate.StoppingConditions Epoch']);
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Propagate.StoppingConditions.Epoch.TripValue ',StopTime,' UTCG']);%
root.ExecuteCommand(['Astrogator */Satellite/blue SetValue MainSequence.SegmentList.Propagate.Propagator Earth_J2']);
Sat_name2='target';
satellite2= root.CurrentScenario.Children.New('eSatellite', Sat_name2);
satellite2.SetPropagatorType('ePropagatorAstrogator'); 
satellite2.Propagator
e=0;Sat_TA=0;Sat_LAN=105;
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList Initial_State Propagate']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.CoordinateType Modified Keplerian']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Epoch ',StartTime,' UTCG']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.ElementType "Kozai-Izsak Mean"']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.Period 86169.6 sec']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.ecc ',num2str(e)]);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.inc 0 deg']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.w 0 deg']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.TA ',num2str(Sat_TA),' deg']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Initial_State.InitialState.Keplerian.LAN ',num2str(Sat_LAN),' deg']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Propagate.StoppingConditions Epoch']);
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Propagate.StoppingConditions.Epoch.TripValue ',StopTime,' UTCG']);%
root.ExecuteCommand(['Astrogator */Satellite/target SetValue MainSequence.SegmentList.Propagate.Propagator Earth_J2']);
root.ExecuteCommand(['Astrogator */Satellite/blue RunMCS']);
root.ExecuteCommand(['Astrogator */Satellite/target RunMCS']);

（2）通过生成报告的方式，来得到蓝方和目标星的位置，根据lambert固定时间拦截的要求，需要对被拦截星末端时刻的位置进行预报，从而确定到达的末位置。

（2.1）初始位置是根据blue星的初始位置生成的，代码如下

data_blue=root.ExecuteCommand(['Report_RM */Satellite/blue Style "Inertial Position Velocity" TimePeriod "',StartTime,'" "',StopTime,'" TimeStep 3600']);
%% 根据蓝方的初始状态得到设置的轨道机动的起点
struct=regexp(data_blue.Item(1),',','split');
rx=str2double(cell2mat(struct(2)));
ry=str2double(cell2mat(struct(3)));
rz=str2double(cell2mat(struct(4)));
vx=str2double(cell2mat(struct(5)));
vy=str2double(cell2mat(struct(6)));
vz=str2double(cell2mat(struct(7)));

（2.2）末端位置则是需要根据任务时间的要求，来预报目标星到达的位置，在这里，设置总的任务时间为3天，即预报目标星在3天后的位置，作为其末端位置。

number_time=4320;%经历4320分钟，即29 Jan 04:00:00.000到达
sec=4320*60;
% 整个任务的时间从26 Jan 2024 04:00:00.000 到 10 Feb 2024 04:00:00.000

mission_day=15;
mission_minute=15*24*60;
% 根据给定的时间，获得其最终位置
data_red=root.ExecuteCommand(['Report_RM */Satellite/target Style "Inertial Position Velocity" TimePeriod "',StartTime,'" "',StopTime,'" TimeStep ',num2str(time_step)]);
starttime='2024-1-26 04:00:00';%开始时间转化为datetime格式,datetime读取不了StartTime的写法
time=datetime(starttime);
minute=minutes(1:mission_minute);
option_fixtime=time+minute;% 这是为了后面变换number_time的数值做准备
% 根据特定的值来获取其末端位置
time_end=datestr(option_fixtime(number_time));
time_end=strrep(time_end,'-',' ');% 与STK的时间格式匹配
time_end=string(time_end)+'.000';% 与STK的时间格式匹配
data_Line=data_red.count;
for i=1:data_Line-2
    struct=regexp(data_red.Item(i),',','split');
    End_time=string(struct(1))
    if End_time==time_end 
        End_rx=cell2mat(struct(2));
        End_ry=cell2mat(struct(3));
        End_rz=cell2mat(struct(4));
        End_vx=cell2mat(struct(5));
        End_vy=cell2mat(struct(6));
        End_vz=cell2mat(struct(7));
        break
    else
        continue
    end
end

（3）调用STK的Design tools，来完成整个问题的求解，代码如下，在STK中打开后的场景如下图所示

root.ExecuteCommand(['ComponentBrowser */ Duplicate "Design Tools" "Lambert Solver" myLambert']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitEpoch "',StartTime,'" UTCG']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitPosVel.Rx ',num2str(rx),' km']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitPosVel.Ry ',num2str(ry),' km']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitPosVel.Rz ',num2str(rz),' km']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitPosVel.Vx ',num2str(vx),' km*sec^-1']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitPosVel.Vy ',num2str(vy),' km*sec^-1']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitPosVel.Vz ',num2str(vz),' km*sec^-1']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert FinlPosVel.Rx ',num2str(End_rx),' km']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert FinlPosVel.Ry ',num2str(End_ry),' km']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert FinlPosVel.Rz ',num2str(End_rz),' km']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert FinlPosVel.Vx ',num2str(End_vx),' km*sec^-1']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert FinlPosVel.Vy ',num2str(End_vy),' km*sec^-1']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert FinlPosVel.Vz ',num2str(End_vz),' km*sec^-1']);
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert MinimumTOF ',num2str(sec), ' sec']);
root.ExecuteCommand(['ComponentBrowser */ LambertCompute "Design Tools" myLambert' ]);

（4）数据处理和分析，为什么要使用STK与matlab交互的最大原因就是这一部分，本文仅讨论如何获取过程的总脉冲大小。代码如下，得到的结果如下图所示

result=root.ExecuteCommand(['ComponentBrowser_RM */ GetValue "Design Tools" myLambert LambertResult']);
%% 获取整个机动过程的总脉冲变量
total_delta=result.Item(23);
V_begin=strfind(total_delta,'=');
V_end=strfind(total_delta,'m*sec^-1');
delta_all=str2num(total_delta(V_begin+1:V_end-1));
fprintf('转移时间为%.4f sec，消耗的总脉冲为%.4f m/s\n',sec,delta_all);

通过这一节，发现在3天机动的情况下，其消耗的脉冲大小很大，很难在实际应用中有所应用，所以，需要修改机动的时间，然而不同的机动时间，消耗的脉冲也是不一样的。为了实现一个较优的控制，下一节将利用matlab中优化算法的方式，综合时间和燃耗这两个因素，实现一个相对较优的机动策略。（本人闲鱼账号：爱stk的鱼酱）如有疑问，欢迎咨询。