STK与matlab交互 Astrogator模块（6）

最新推荐文章于 2024-06-14 00:39:34 发布

矢志航天的阿洪

最新推荐文章于 2024-06-14 00:39:34 发布

阅读量546

点赞数 11

文章标签： matlab 开发语言

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

版权

一、背景介绍

在（5）的结尾中，给出了一个关于入轨点优化的情况，即只要在绕飞椭圆上的某一点的位置和速度不受外力的情况下，将会一直保持绕飞状态，一个自然绕飞周期通常对应的是目标星的周期，因此在本文中，将在matlab中使用CW方程递推，VVLH坐标系下的相对位置和速度，然后利用matlab与STK交互，将该数据放入STK中实现从相对坐标系到绝对坐标系的转化，即该任务时刻，5 Feb 2024 04:00:00.000 UTCG能够实现绕飞的到达位置有很多种。现在将该位置的选择和转移的时间作为同时考虑的两个变量进行优化。

二、CW方程

CW方程有很多参考的地方，其解析形式也在很多参考书上明确列出，作者使用matlab对该公式解析化的整个过程进行了一个编程，感兴趣的同学可以作为一个参考。

%% 常微分方程的解析推导
syms n t
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;0 0 0 0 0 2*n;0 -n^2 0 0 0 0;0 0 3*n^2 -2*n 0 0];
[V,J]=jordan(A);
V_inv=inv(V);
%% exp(J*t)的矩阵形式
exp_J=[1 t 0 0 0 0;0 1 0 0 0 0;0 0 cos(n*t)-sin(n*t)*i 0 0 0;0 0 0 cos(n*t)+sin(n*t)*i 0 0;0 0 0 0 cos(n*t)-sin(n*t)*i 0;0 0 0 0 0 cos(n*t)+sin(n*t)*i];
exp_A=V*exp_J*V_inv;
simplify(exp_A);

因此，该转移方程代码如下所示

function X_t=CW(X_0,n,t)
% X_0 为初始状态
% n   为角速度
% t   为时间
Matrix=[1,           0,  6*n*t - 6*sin(n*t), (4*sin(n*t) - 3*n*t)/n,          0, -(2*cos(n*t) - 2)/n;
        0,    cos(n*t),                   0,                      0, sin(n*t)/n,                   0;
        0,           0,      4 - 3*cos(n*t),     (2*cos(n*t) - 2)/n,          0,          sin(n*t)/n;
        0,           0, -6*n*(cos(n*t) - 1),         4*cos(n*t) - 3,          0,          2*sin(n*t);
        0, -n*sin(n*t),                   0,                      0,   cos(n*t),                   0;
        0,           0,        3*n*sin(n*t),            -2*sin(n*t),          0,            cos(n*t)];
X_t=Matrix*X_0;
end

由此根据（5）中给出能绕飞的VVLH位置和速度使用VVLH方程递推，所有的入轨点均能实现椭圆绕飞，将这些点存储为一个集合。

miu=3.986e5;
omega=2*pi/86169.6
z=10;    %绕飞短半轴
dx=2*omega*z; 
%% 利用CW方程选择能绕飞的点，根据理论知识，VVLH下绕飞椭圆的每个状态再不施加外部脉冲的时候，其绕飞椭圆均保持不变。
%% 根据初始条件 目标星的周期为86169.6sec 为了和前面报表的时间保持一致，本文取的时间步长为60sec
t=0:60:86169.6;
X=[0,0,z,dx,0,0]';xx=zeros(length(t),6);
%% xx的长度 即包含了n个入轨点，也意味着可选择入轨的点包括这n个选择，现在从这n个入轨点进行优化和选择
for i=1:length(t)
    t=(i-1)*60;
    xx(i,:)=CW(X,omega,t);
end

三、优化函数

根据（5）以及第二部分的介绍，本文的优化量包括末位置的入轨点选择和开始实施机动的时间，由于涉及到两个变量，本文使用的优化函数使用fminsearch。优化函数的代码如下所示

function J=lambert3(x)
% x(1) 为绕飞选择的时间
% x(2) 为选择入轨点

global  root xx
%% 通过改变入轨点，修改末端进入的位置和速度
indexdot=floor(x(2));
insert_dot=xx(indexdot,:);
rrxx=insert_dot(1);
rryy=insert_dot(2);
rrzz=insert_dot(3);
vvxx=insert_dot(4);
vvyy=insert_dot(5);
vvzz=insert_dot(6);
StartTime='26 Jan 2024 04:00:00.000';
MissionTime = '5 Feb 2024 04:00:00.000';  
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 ',num2str(rrxx),' km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Y ',num2str(rryy),' km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Z ',num2str(rrzz),' km']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vx ',num2str(vvxx),' km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vy ',num2str(vvyy),' km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Cartesian.Vz ',num2str(vvzz),' km/sec']);
root.ExecuteCommand(['Astrogator */Satellite/raofei SetValue MainSequence.SegmentList.Initial_State.InitialState.Epoch ',MissionTime,' 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 "',MissionTime,'" "',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)));
%% 通过改变时间推移，修改初始的位置和速度

time_step=60;%为后面的固定步长遍历来寻找最优来做准备
t_min=floor(x(1));
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

运行程序的主函数如下所示

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

miu=3.986e5;
omega=2*pi/86169.6
z=10;    %绕飞短半轴
dx=2*omega*z; 
%% 利用CW方程选择能绕飞的点，根据理论知识，VVLH下绕飞椭圆的每个状态再不施加外部脉冲的时候，其绕飞椭圆均保持不变。
%% 根据初始条件 目标星的周期为86169.6sec 为了和前面报表的时间保持一致，本文取的时间步长为60sec
t=0:60:86169.6;
X=[0,0,z,dx,0,0]';xx=zeros(length(t),6);
%% xx的长度 即包含了n个入轨点，也意味着可选择入轨的点包括这n个选择，现在从这n个入轨点进行优化和选择
for i=1:length(t)
    t=(i-1)*60;
    xx(i,:)=CW(X,omega,t);
end

satellite3= root.CurrentScenario.Children.New('eSatellite', VMC_SatName);
satellite3.SetPropagatorType('ePropagatorAstrogator'); 
satellite3.Propagator;root.ExecuteCommand(['ComponentBrowser */ Duplicate "Design Tools" "Lambert Solver" myLambert']);

%% 利用函数优化的方法找到最优机动时刻以及最优入轨点，已知初始出发时间为26 Jan 2024 04:00:00.000 
options = optimset('Display','iter','MaxIter',200,'PlotFcns',@optimplotfval)
[x,fval]=fminsearch(@lambert3,[9000,4],options)
J=lambert3(x)
%% 将最优结果呈现在STK中
root.ExecuteCommand(['ComponentBrowser */ SetValue "Design Tools" myLambert SequenceName ones']);
root.ExecuteCommand(['ComponentBrowser */ LambertConstructSequence "Design Tools" myLambert']);
root.ExecuteCommand(['ComponentBrowser */ LambertAddToCB "Design Tools" myLambert']);
satellite1.Unload
Sat_best='blue_best';
satellite4= root.CurrentScenario.Children.New('eSatellite', Sat_best);
satellite4.SetPropagatorType('ePropagatorAstrogator'); 
satellite4.Propagator;
time_maneuver=root.ExecuteCommand(['ComponentBrowser_RM */ GetValue "Design Tools" myLambert InitEpoch']);
str_start=strfind(time_maneuver.Item(0),'=')
str_end=strfind(time_maneuver.Item(0),'UTCG')
time_ma=time_maneuver.Item(0)
time=time_ma(str_start+2:str_end-2)
root.ExecuteCommand(['Astrogator */Satellite/blue_best SetValue MainSequence.SegmentList ones Propagate']);
root.ExecuteCommand(['Astrogator */Satellite/blue_best SetValue MainSequence.SegmentList.Propagate.StoppingConditions Epoch']);
root.ExecuteCommand(['Astrogator */Satellite/blue_best SetValue MainSequence.SegmentList.Propagate.StoppingConditions.Epoch.TripValue ',StopTime,' UTCG']);%
root.ExecuteCommand(['Astrogator */Satellite/blue_best SetValue MainSequence.SegmentList.Propagate.Propagator Earth_J2']);
root.ExecuteCommand(['Astrogator */Satellite/blue_best RunMCS']);

%% 找出机动时间与消耗脉冲的关系
%% 不显示红方和蓝方的固定坐标系下的轨迹
root.ExecuteCommand(['VO */Satellite/raofei OrbitSystem Modify System "InertialByWindow" Show Off']);
% root.ExecuteCommand(['VO */Satellite/blue OrbitSystem Modify System "InertialByWindow" Show Off']);
root.ExecuteCommand(['VO */Satellite/raofei OrbitSystem Add System VVLH Satellite/target Color red']);
root.ExecuteCommand(['VO */Satellite/target OrbitSystem Modify System "InertialByWindow" Show Off']);
root.ExecuteCommand(['VO */Satellite/blue_best OrbitSystem Add System VVLH Satellite/target Color red']);
%%

四、结果分析

经过最优计算，得到的最优机动时间为1 Feb 2024 12:15:00.000 UTCG，到达的VVLH点为

矢志航天的阿洪

关注

11
点赞
踩
9

收藏

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

在（5）的结尾中，给出了一个关于入轨点优化的情况，即只要在绕飞椭圆上的某一点的位置和速度不受外力的情况下，将会一直保持绕飞状态，一个自然绕飞周期通常对应的是目标星的周期，因此在本文中，将在matlab中使用CW方程递推，VVLH坐标系下的相对位置和速度，然后利用matlab与STK交互，将该数据放入STK中实现从相对坐标系到绝对坐标系的转化，即该任务时刻，5 Feb 2024 04:00:00.000 UTCG能够实现绕飞的到达位置有很多种。现在将该位置的选择和转移的时间作为同时考虑的两个变量进行优化。
复制链接

扫一扫