STK与matlab交互 Astrogator模块（5）

最新推荐文章于 2024-05-11 23:42:30 发布

矢志航天的阿洪

最新推荐文章于 2024-05-11 23:42:30 发布

阅读量721

点赞数 23

文章标签： matlab 开发语言

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

版权

一、背景介绍

在（4）中，我们介绍了在确定初始时刻即为机动时刻的情况下，通过改变机动的时间从而确定最优转移时间的方法，在本节中，考虑实际的应用的情况，通常实现交会对接任务的时候，是不会直接达到目标星的位置，而是到达其相对距离较近的地方，进行下一步的操作。在本节中，利用相对位置坐标系来确定最终目标的位置。

二、仿真步骤

（1）在本节中，给定到达的时刻为05 Feb 2024 04:00:00.000 UTCG，到达的位置在其VVLH坐标系下面的坐标为 $(0,0,10km,1.4589m/s,0,0)$ ，根据CW方程的定义，该卫星之后将以绕飞的形式环绕在目标星的周围。在这里涉及到VVLH坐标系和地球惯性坐标系的转换，在本文中，直接利用STK构造一个在该时刻满足该位置的卫星，得到其在该时刻地球惯性系的报告数据，即可得到该时刻的惯性坐标系的位置，即整个机动过程的末状态。

下面代码是构建该绕飞卫星，并获取轨道机动末端位置的整个过程。

VMC_SatName='raofei';
satellite3= root.CurrentScenario.Children.New('eSatellite', VMC_SatName);
satellite3.SetPropagatorType('ePropagatorAstrogator'); 
satellite3.Propagator;
Mission_time='5 Feb 2024 04:00:00';% 在给定任务的第十天要求到达该位置
StopTime = '10 Feb 2024 04:00:00.000';     % 场景结束时间
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList Initial_State Propagate']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.CoordinateSystem "Satellite/target VVLH"'])
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.CoordinateType Cartesian'])
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.X 0 km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Y 0 km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Z 10 km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vx ',num2str(dx),' km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vy 0 km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vz 0 km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Epoch ',Mission_time,' UTCG']);
root.ExecuteCommand(['Astrogator */Satellite/raofei RunMCS']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Propagate.StoppingConditions Epoch']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Propagate.StoppingConditions.Epoch.TripValue ',StopTime,' UTCG']);%
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Propagate.Propagator Earth_J2']);
data_raofei=root.ExecuteCommand(['Report_RM */Satellite/raofei Style "Inertial Position Velocity" TimePeriod "',Mission_time,'" "',StopTime,'" TimeStep 3600']);
struct=regexp(data_raofei.Item(1),',','split');
End_rx=str2double(cell2mat(struct(2)));
End_ry=str2double(cell2mat(struct(3)));
End_rz=str2double(cell2mat(struct(4)));
End_vx=str2double(cell2mat(struct(5)));
End_vy=str2double(cell2mat(struct(6)));
End_vz=str2double(cell2mat(struct(7)));

（2）参考（4）的方法，构造时间和脉冲大小的函数，这里的时间指卫星在初始状态到第一次施加脉冲之间的时间，本文中均以一分钟作为时间步长。其设计思路即将每一时刻的速度和位置记录下来，通过STK自带的LAMBERT问题求解器，得到该时刻对应的脉冲大小。其代码如下：

function J=lambert2(t_min)
% t 即设定的转移时间 单位minute
global End_rx End_ry End_rz End_vx End_vy End_vz root
time_step=60;%为后面的固定步长遍历来寻找最优来做准备
t_min=floor(t_min);
StartTime='26 Jan 2024 04:00:00.000';
MissionTime = '5 Feb 2024 04:00:00.000';  
t_all=10*24*60;
sec=(14400-t_min)*60;
data_blue=root.ExecuteCommand(['Report_RM */Satellite/blue Style "Inertial Position Velocity" TimePeriod "',StartTime,'" "',MissionTime,'" TimeStep ',num2str(time_step)]);
data_Line=data_blue.count;
for i=1:data_Line-2
    struct=regexp(data_blue.Item(i),',','split');
    Start_time(i)=struct(1);
end
Start_time_str=string(Start_time(t_min+1));% 找到对应时间对应的报告
for i=1:data_Line-2
    struct=regexp(data_blue.Item(i),',','split');
    time_start=string(struct(1));
    if time_start==Start_time_str
        rx=cell2mat(struct(2));
        ry=cell2mat(struct(3));
        rz=cell2mat(struct(4));
        vx=cell2mat(struct(5));
        vy=cell2mat(struct(6));
        vz=cell2mat(struct(7));
        break
    else
        continue
    end
end

root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert InitEpoch "',char(time_start),'" 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' ]);

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));
J=delta_all;
fprintf('%.4f\n',J)
end

（3）利用matlab自带的fminbnd函数对其进行优化求解，设置时间的上界为9天，下界为0。主函数的代码如下：

clear;clc
%% 本文将实现兰伯特制导，末位置为对目标绕飞所需要的时间和地点
uiApplication = actxGetRunningServer('STK11.application');
% Get our IAgStkObjectRoot interface
global root End_rx End_ry End_rz End_vx End_vy End_vz
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']);

%% 给出到达目标星的时间 确定绕飞的位置 
Mission_time='5 Feb 2024 04:00:00';% 在给定任务的第十天要求到达该位置
%% 利用STK报表的方式生成一个该时刻绕飞的卫星
VMC_SatName='raofei';
R_norm=sqrt(End_rx^2+End_ry^2+End_rz^2);
V_norm=sqrt(End_vx^2+End_vy^2+End_vz^2);
R_vec=[End_rx,End_ry,End_rz]/R_norm;
V_vec=[End_vx,End_vy,End_vz]/V_norm;
Z_vec=-R_vec;
Y_vec=cross(Z_vec,V_vec);
X_vec=cross(Y_vec,Z_vec);
miu=3.986e5;
omega=sqrt(miu/R_norm^3);
z=10;    %绕飞短半轴
dx=2*omega*z; 
satellite3= root.CurrentScenario.Children.New('eSatellite', VMC_SatName);
satellite3.SetPropagatorType('ePropagatorAstrogator'); 
satellite3.Propagator;
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList Initial_State Propagate']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.CoordinateSystem "Satellite/target VVLH"'])
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.CoordinateType Cartesian'])
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.X 0 km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Y 0 km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Z 10 km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vx ',num2str(dx),' km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vy 0 km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vz 0 km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Epoch ',Mission_time,' UTCG']);
root.ExecuteCommand(['Astrogator */Satellite/raofei RunMCS']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Propagate.StoppingConditions Epoch']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Propagate.StoppingConditions.Epoch.TripValue ',StopTime,' UTCG']);%
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Propagate.Propagator Earth_J2']);
%% 该绕飞卫星的初始时刻位置，即为LAMBERT航天器需要到达的末端位置
data_raofei=root.ExecuteCommand(['Report_RM */Satellite/raofei Style "Inertial Position Velocity" TimePeriod "',Mission_time,'" "',StopTime,'" TimeStep 3600']);
struct=regexp(data_raofei.Item(1),',','split');
End_rx=str2double(cell2mat(struct(2)));
End_ry=str2double(cell2mat(struct(3)));
End_rz=str2double(cell2mat(struct(4)));
End_vx=str2double(cell2mat(struct(5)));
End_vy=str2double(cell2mat(struct(6)));
End_vz=str2double(cell2mat(struct(7)));
root.ExecuteCommand(['ComponentBrowser */ Duplicate "Design Tools" "Lambert Solver" myLambert']);

%% 利用函数优化的方法找到最优机动时刻，已知初始出发时间为26 Jan 2024 04:00:00.000 
lb=0;
ub=9*24*60;
options = optimset('Display','iter','MaxIter',200,'PlotFcns',@optimplotfval)
[x,fval]=fminbnd(@lambert2,lb,ub,options)
lambert2(x)

（4）最终对结果进行验证,上图是绕飞卫星的位置和速度，下图是通过lambert求解器后，按照lambert轨道机动得到的最终位置和速度，其VVLH坐标系下的形式如下图所示。由于LAMBERT求解器的误差，所以最终到的蓝方绕飞星并不能实现真正的绕飞，因此在这里还需要进行修正，求出的最终结果是在29 Jan 2024 12:27:00.000 UTCG 时机动，需要施加的脉冲大小为5237.6007m/s

三、总结与展望

通过（5）的操作，最终实现了绕飞的整个过程，但在本文中开始绕飞的入轨点是事先给定的，理论上来说，只要在VVLH坐标系的绕飞椭圆里的任一状态都能够保持绕飞状态，因此，最终的状态点也能据此进行优化，就意味着初状态和末状态都是变化的，初状态与实施机动的时间有关，而末状态与其进入绕飞的轨迹点有关。下一节将会围绕这个问题进行进一步的讨论。创作不易，感兴趣的朋友欢迎咨询闲鱼账号：爱stk的鱼酱，将提供1对1答疑服务

矢志航天的阿洪

关注

23
点赞
踩
11

收藏

觉得还不错? 一键收藏
2
评论
STK与matlab交互 Astrogator模块（5）

通过（5）的操作，最终实现了绕飞的整个过程，但在本文中开始绕飞的入轨点是事先给定的，理论上来说，只要在VVLH坐标系的绕飞椭圆里的任一状态都能够保持绕飞状态，因此，最终的状态点也能据此进行优化，就意味着初状态和末状态都是变化的，初状态与实施机动的时间有关，而末状态与其进入绕飞的轨迹点有关。下一节将会围绕这个问题进行进一步的讨论。创作不易，感兴趣的朋友欢迎咨询闲鱼账号：爱stk的鱼酱，将提供1对1答疑服务。
复制链接

扫一扫