本文介绍了如何在航天工程中使用fminsearch函数优化无约束的多变量函数，以找到目标星到达位置的最佳转移时间。通过Matlab与STK的交互，设计了函数并讨论了遗传算法的局限性，提出了将Matlab自编的兰伯特问题求解器集成到STK中的改进方向。
一、优化算法
由于整个操作过程涉及到与STK交互，使用常规的遍历找最优和遗传算法产生大量的数据通常是不合适的，本文主要介绍了fminsearch函数，使用无导数法计算无约束的多变量函数的最小值。fminsearch有几个输入参数，结构体通常包含的参数：PlotFcns：绘制算法执行过程种的各个进度测量值。
二、函数设计
（1）根据分析知道，时间带来的变化是目标星到达位置的不同，STK能够预报出每分钟的位置和速度，利用STK和matlab交互，将报告里的数据导入matlab，然后根据通过修改转移的时间，得到的末端位置和速度均不同，通过（3）介绍的互联功能，将得到的数据写入LAMBERT求解器，得到耗费的总脉冲大小。由于整个分析时间很大，我们取每一步长为一分钟，因此在matlab中产生报表的命令如下：
data_red=root.ExecuteCommand(['Report_RM */Satellite/target Style "Inertial Position Velocity" TimePeriod "',StartTime,'" "',StopTime,'" TimeStep ',num2str(time_step)]);
% time_step这里取值为60（2）由于STK报告的格式，固定时间转移通常需要将时间转化为时刻，即根据初始时刻26 Jan 04:00:00.000，开始计算固定的时间对应的到达时刻，例如：固定时间为1天，那么其对应的时刻即为27 Jan 04:00:00.000。从而得到该时刻的位置和速度，在函数设计中，本文直接将所有时间与时刻的字符串序列一一对应起来
data_Line=data_red.count;
for i=1:data_Line-2
struct=regexp(data_red.Item(i),',','split');
End_time(i)=struct(1);
end
End_time_str=string(End_time(t_min+1));% 找到对应时间对应的报告
for i=1:data_Line-2
struct=regexp(data_red.Item(i),',','split');
time_end=string(struct(1));
if time_end==End_time_str
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）因此整个问题的唯一变量是时间，评价的指标设为总脉冲数量，因此本文直接给出的函数为
function J=lambert(t_min)
% t 即设定的转移时间 单位minute
global rx ry rz vx vy vz root
time_step=60;%为后面的固定步长遍历来寻找最优来做准备
% 整个任务的时间从26 Jan 2024 04:00:00.000 到 10 Feb 2024 04:00:00.000
mission_day=15;
mission_minute=15*24*60;
% 有一个时间的计算
t_min=floor(t_min);
StartTime='26 Jan 2024 04:00:00.000';
StopTime = '10 Feb 2024 04:00:00.000';
sec=t_min*60;
data_red=root.ExecuteCommand(['Report_RM */Satellite/target Style "Inertial Position Velocity" TimePeriod "',StartTime,'" "',StopTime,'" TimeStep ',num2str(time_step)]);
data_Line=data_red.count;
for i=1:data_Line-2
struct=regexp(data_red.Item(i),',','split');
End_time(i)=struct(1);
end
End_time_str=string(End_time(t_min+1));% 找到对应时间对应的报告
for i=1:data_Line-2
struct=regexp(data_red.Item(i),',','split');
time_end=string(struct(1));
if time_end==End_time_str
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
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' ]);
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（4）在得到了需要优化的函数后，使用优化算法对其进行优化，由于整个问题只涉及到一个变量，可以考虑使用（一）中介绍的fminsearch函数，设置初始想定的最优转移时间为9000分钟，上界和下界分别设为1天和15天，下面给出主函数及运行结果图。
clear;clc
%% 利用STK自带的LAMBERT库来实现目标拦截
uiApplication = actxGetRunningServer('STK11.application');
% Get our IAgStkObjectRoot interface
global root
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']);
%% 生成报告，得到轨道数据
data_blue=root.ExecuteCommand(['Report_RM */Satellite/blue Style "Inertial Position Velocity" TimePeriod "',StartTime,'" "',StopTime,'" TimeStep 3600']);
%% 根据蓝方的初始状态得到设置的轨道机动的起点
struct=regexp(data_blue.Item(1),',','split');
global rx ry rz vx vy vz
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)));
root.ExecuteCommand(['ComponentBrowser */ Duplicate "Design Tools" "Lambert Solver" myLambert']);
%% 遗传算法计算最优机动时刻
lb=24*60;
ub=15*24*60;
options = optimset('Display','iter','MaxIter',200,'PlotFcns',@optimplotfval)
[x,fval]=fminsearch(@lambert,9000,options)
三、不足与改进
本节设定的情况是在蓝星初始时刻机动的情况下，通过改变转移时间来确定的最优转移方案，步长也是选择60s为一步进行报告生成，同时由于matlab与STK交互的速度限制，该函数几乎很难用遗传算法等进化算法。之后将考虑将matlab自编算法实现兰伯特问题求解，再将其代入STK中仿真实现。下一节将会考虑在确定任务时刻和最终点的情况下，通过修改初始机动时间来优化实现最优转移。
1871
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
