Propagation Design Issues

[DJCinSB]

It does look like the SPICE SPK file supports multiple objects for propagation at once. Given that feature, we need to be able to identify which object in the SPK file is used for each SPK based object in GMAT. I think that means we need an identifier for each object we pull from the SPICE file -- looks like that is either a character string or a NAIF ID. That would be a new (and often unused) parameter for the SpacePoint class, since we want to have the ability, eventually, to propagate any of our SpacePoints -- celestial bodies (so we can do comets and asteroids and moons from SPICE files), spacecraft, ground stations, and so forth.

I thought about and discarded the idea of using the object's name as the identifier in the ephemeris file. Using a separate field from the object's name gives us the flexibility to use a script specific name rather than being constrained by the identifier in the SPK file -- or other ephemeris source, if we add support for other types of precalculated ephemerides.

Another alternative is to associate the file with the SpacePoint, rather than with the propagator. We'd still need the identifier in this case, though. One other wrinkle -- if I'm reading the SPICE specs correctly, we also need to be able to identify a second ephem source for the location of the reference body (typically the Earth) for the SPICE based propagations at least, unless we can guarantee that that body is already in the SPICE file. (I'm deriving this from statements I'm reading that say that SPICE files give the location of one body with respect to another, so they need to be able to find both bodies in a SPICE compatible format at the desired epoch.) I'll call this second reference the SSEphem in the scripting below. If we take the "file on SpacePoint" approach, the scripting would be something like this:

Code: Select all

Create Spacecraft Clementine
Clementine.EphemID = -40   % Clementine's NAIF ID
Clementine.Ephemeris = DSPSE.SPK
Clementine.SSEphem = SolarBodies.SPK

Create Asteroid Geographos
Geographos.EphemID = 2006513   % No clue about the real NAIF ID for Geographos...
Geographos.Ephemeris = DSPSE.SPK
Geographos.SSEphem = SolarBodies.SPK

Create Propagator prop
prop.type = SPICE
prop.StepSize = 300    % seconds

Create Variable i
For i = 0 : 2000
   Propagate prop(Clementine, Geographos)
EndFor


I think I prefer associating the file with the propagator, though, like this:

Code: Select all

Create Spacecraft Clementine
Clementine.EphemID = -40   % Clementine's NAIF ID

Create Asteroid Geographos
Geographos.EphemID = 2006513   % No clue about the real NAIF ID for Geographos...

Create Propagator prop
prop.type = SPICE
prop.StepSize = 300    % seconds
prop.Ephemeris = DSPSE.SPK
prop.SSEphem = SolarBodies.SPK

Create Variable i
For i = 0 : 2000
   Propagate prop(Clementine, Geographos)
EndFor


One downside to this approach: if we set the ephem sources on the propagator, that forces a different propagator when propagating multiple spacecraft where they are all in different ephemeris files. For SPICE, though, there are tools to merge the SPK files together.

[shughes]

Conrad is very busy these days, and we discussed some items that I think must be resolved soon, so I'll summarize our discussion here. Sorry if I misrepresent your thoughts Conrad!!!

By numerically integrating the equations of motion, you end up with a set of discrete times at which you know the states. Different components of an analysis system often need the states at times that are not the same as the set of times naturally chosen by the integrator: OD measurement models, product generator, and discrete event tracking among others.

We need to determine in general how GMAT will provide the state at times not contained in the set determined during numerical integration. It's not clear if this should always be handled by the propagator, but I (Steve, not Conrad) think that there are some times when the propagator will need to handle it. Here are three approaches used currently:

• Dense output integrator (suggested by Montenbruck and Gill)
• Interpolation (Used in GTDS, GEODYN, and suggested by Montenbruck and Gill
• Propagate to necessary times (I believe FreeFlyer does this for eclipse entry/exit times. )

At a high level, I (Steve, not Conrad) think the propagator would propagate to some criteria, perhaps just one step, and publish the state to different components that require it. These components would determine if they need to do anything, i.e, calculate a measurement, determine root crossing, etc. Then the component would send back a time, or a guess for the time that it needs the state. The propagator, or some other component would use this information and one of the three methods above to determine the state at the requested time. This could be iterative in the case of root finding.

Well there's lot's more to discuss here, but that's enough for a first post.... Thoughts?

[DJCinSB]

UPDATE: This message and the preceding one have been moved into a separate thread: Pause Points in Propagation (see http://gmat.ed-pages.com/forum/viewtopic.php?f=3&t=104&start=0&st=0&sk=t&sd=a).


I think we may want to break this into its own discussion -- if you agree, I'll move the messages around so we don't lose anything.

My first thought here is that we should start by implementing an interpolator approach and see how well it performs. I think that is likely the simplest first step. If we do that, there are refinements we'll need to figure out as well -- for instance, suppose we have stepped a formation past a desired point, and then interpolate the triggering epoch. Do we then prop to that epoch? I think that depends on what the point is used for -- if it's just a reporting event (shadow entry, entry into visibility at a station, etc) then we might not need to do propagation, depending on the accuracy needed for the triggering event. If it's a point where some action occurs -- start a maneuver, for example -- then we should do the prop step.

I do think that propagating to the necessary times is probably overkill, particularly if it means multiple prop steps to find the desired epoch. We should be able to do better than that via interpolation, possibly adding a single prop to be close for the final interpolation step.

I would like to add a dense output integrator. But I'm being selfish there -- I'm not convinced it is needed for GMAT's use, I just like implementing integrators!

[shughes]

ok, if you split out this thread, I'll make comments to your latest post.

[DJCinSB]

We need to be able to prop the state vector to each measurement epoch. What if we add a estimator-specific stopping condition for measurement data, something like this:

Code: Select all

Estimate batchExt
...
   Propagate prop(sat1, sat2) {ToMeasurement}
EndEstimate

that lets the estimator drive propagation to a measurement epoch. We would also, of course, need the ability to interpolate, integrate to, or dense output generate state information at the light time delayed epochs. But we do also need some way to stop at measurement epochs, and this seems like one way to do that part -- basically, the Propagate command (or other prop enabled command if we open that part up) would ask the estimator for the next measurement epoch and prop to it. Then we might script a full estimation process something like this:

Code: Select all

Create Spacecraft ODSat;
Create Propagator prop;
Create GroundStation WhiteSands

%----------------------------------------
%---------- Input Data Specification
%----------------------------------------
Create MeasurementFile ODSatTrackData
ODSatTrackData.Filename = /home/mwilkins/data/moredata.txt;
ODSatTrackData.FileFormat = SLR;
ODSatTrackData.MeasurementTypes = {'WhiteSands.TDRSS2WayRange', 'WhiteSands.TDRSS2WayDoppler'};

%----------------------------------------
%---------- Estimator Properties
%----------------------------------------
Create BatchEstimator BatchLS
BatchLS.MaximumIterations = 25;
BatchLS.AbsoluteConvergenceTolerance = 1e-5;
BatchLS.SolutionEpoch = 01 Jan 2008 12:00:00.000;
BatchLS.Measurements  = {'ODSatTDRSS'}

%---- The estimation loop
Estimate BatchLS

   % Define the measurements
   Measurement BatchLS('ODSatTDRSS');

   %  Define solvefor, consider, and neglect parameters for parameters
   % These data are initialized to the guess for the spacecraft state at the estimation epoch:
   SolveFor BatchLS(ODSat.X = 7200, RelTol = 1)
   SolveFor BatchLS(ODSat.Y = 0, RelTol = 1)
   SolveFor BatchLS(ODSat.Z = 1000, RelTol = 1)
   SolveFor BatchLS(ODSat.VX = 7.3, RelTol = 0.1)
   SolveFor BatchLS(ODSat.VY = 1, RelTol = 0.1)
   SolveFor BatchLS(ODSat.VZ = 2, RelTol = 0.1)
           
   Propagate LowEarthProp(ODSat) {ToMeasurement};

EndEstimate

[shughes]

We need to be able to prop the state vector to each measurement epoch. What if we add a estimator-specific stopping condition for measurement data, something like this:

Code: Select all

   Propagate LowEarthProp(ODSat) {ToMeasurement};


If I'm understanding correctly, the above line would essentially function like a loop or sequence of propagates? If there were 100 measurements, then this line would propagate to all one hundred times. I'm not sure how to extend this to more complicated sequences that have maneuvers etc. Here's a first cut

Code: Select all


   Vary BatchLS( dV.V = 1 );

   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes,ODSat.UTCModJulian = ManeuverEpoch};
   Maneuver dV(ODSat)
   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes};



One concern I have is that this treats some items in the list as pause conditions and others as stop conditions, which may be confusing

[DJCinSB]

We need to be able to prop the state vector to each measurement epoch. What if we add a estimator-specific stopping condition for measurement data, something like this:

Code: Select all

   Propagate LowEarthProp(ODSat) {ToMeasurement};


If I'm understanding correctly, the above line would essentially function like a loop or sequence of propagates? If there were 100 measurements, then this line would propagate to all one hundred times.

Yes, pausing for other processing at each measurement.

I'm not sure how to extend this to more complicated sequences that have maneuvers etc. Here's a first cut

Code: Select all


   Vary BatchLS( dV.V = 1 );

   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes,ODSat.UTCModJulian = ManeuverEpoch};
   Maneuver dV(ODSat)
   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes};



One concern I have is that this treats some items in the list as pause conditions and others as stop conditions, which may be confusing

I was thinking we'd do this type of thing as a multiple stopping condition construction like you have in the first Propagate command. So this line:

Code: Select all

   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes, ODSat.UTCModJulian = ManeuverEpoch};


is pair of stopping conditions, with intermediate pauses at each measurement epoch. The propagation will stop when either (1) the last measurement epoch is reached or (2) the maneuver epoch is reached, whichever one comes first.

There is some coordination that needs to be managed to be sure that the paired Propagates in your example work correctly. What I mean is that with this scripting:

Code: Select all

   Vary BatchLS( dV.V = 1 );

   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes,ODSat.UTCModJulian = ManeuverEpoch};
   Maneuver dV(ODSat)
   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes};


we want to be sure that if the last measurement occurred before the maneuver epoch, then we don't do any propagation afterwards. We also need to do something to propagate to the maneuver epoch, so I'd probably script it like this if I'm not sure about where the epochs fall:

Code: Select all

   Vary BatchLS( dV.V = 1 );

   Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes, ODSat.UTCModJulian = ManeuverEpoch};
   If ODSat.UTCModJulian >= ManeuverEpoch  % >= in case we are microsecs past the maneuver epoch
      Maneuver dV(ODSat)
      Propagate LowEarthProp(ODSat) {BatchLS.MeasurementTimes};
   EndIf


[shughes]

I like the way Matlab handles what they call Event Functions, and think we should consider using the standard they have developed for the stopping condition issues above. If we went this way, I think the issue of stopping on measurement times, eclipses, or whatever would be straight forward from a user perspective. Below is a brief explanation of what Matlab does... see the URL at the bottom for Matlab's Documentation.

In Matlab, the user can provide as part of the odeset structure a function of the following form:

[VALUE, ISTERMINAL, DIRECTION] = MyEventFunction(T, Y)

T is the integration time,Y is the integration state vector. VALUE is the EventFunction value, ISTERMINAL is a flag that tells matlab whether the integration should stop when a root crossing is found, and DIRECTION describes the root crossing direction that the ode solvers should search for. DIRECTION = -1 means track decreasing function root values, 0 means find all roots, and 1 means find all increasing function root value.

During propagation, matlab tracks the root crossings as requested and records their times and the state values at the event times. If we treat measurements, eclipses, rise-set, as Matlab does, we remove the complexity of the stopping condition definition ( for complex problems ) from propagator and place it on EventFunction. This is contingent upon having something like EventManager that I mentioned in an offline email to you Darrel.

I'll provide an example in another post.


http://www.mathworks.com/access/helpdes ... mit=Search

[shughes]

Here is a discrete event example for eclipse

Code: Select all

Create DiscreteEvent EarthSatEclipse
EarthSatEclipse.Type = 'UmbralEclipse';
EarthSatEclipse.File = Shadows.txt;
EarthSatEclipse.ShadowBody = 'Earth';
EarthSatEclipse.Secondary  = 'Sat';
EarthSatEclipse.Action     = {Report}
EarthSatEclipse.Direction  = 0;
EarthSatEclipse.Active      = true

For this discrete event function, UmbralEclipse would be an internal function based on distance from the umbral shadow cone or some other quantity. GMAT would track Earth shadows on Spacecraft named Sat. When an eclipse was encountered data for the event would be written to Shadows.txt, and we would track Umbra entry and exit. GMAT would continue propagation for this event. If the line said

Code: Select all

EarthSatEclipse.Action     = {Report,Stop}

then the report line would be written and propagation would stop.

[shughes]

Here is a proposed measument model example

Code: Select all

Create Measurement SatGSRange
SatGSRange.Type         = TwoWayRange
SatGSRange.FileFormat   = B3;
SatGSRange.ObservationDataFile =  .\Meas\SatCanberraRange.txt
SatGSRange.StopTime     = 01 Jan 2000 12:30:00.000;
SatGSRange.StartTime    = 02 Jan 2000 12:30:00.000;
SatGSRange.Corrections  = {LightTime, Troposphere, Ionosphere};
SatGSRange.Participants = {Sat, Canberra};
SatGSRange.Sensors      = {Transponder1, Antenna1};
SatGSRange.SimulationTimeStep  = 180;
SatGSRange.SimluationDataFile  = .\Meas\SimulSatCanberraRange.txt
SatGSRange.Action              = Continue; %  Options are: {Continue, Stop, Report....) ;

Type is the measurement type.. there would be many of these!!
FileFormat the format of the file containing the observed values
ObservationDataFile the file containing the observed quantities
StartTime Optional field to specify the start time of the measurements. When in data simulation mode this flag would tell when GMAT should start simulating measurements. Or it could tell GMAT to only use a range of measurements from the input file.
StopTime Optional field to specify the start time of the measurements. When in data simulation mode this flag would tell when GMAT should stopsimulating measurements. Or it could tell GMAT to only use a range of measurements from the input file.
Corrections Tells GMAT what measurement corrections to apply and what models to use.
Participants Tell GMAT what system resources are used in the measurement. These would be the inputs the to specific measurement function describe the Type field. For different measurement types, this list would be different.
Sensors Tells GMAT what sensors to use if they are sensors on the participants. If not, there would be a default for the measurement type.
SimulationTimeStep If in data simulation mode, This would be the frequency of the measurement simulation.
SimluationDataFile If in data simulation mode, this would be the data file that the simulated measurements are written to.
Action Telss GMAT what to do when a measurement is reached. This would likely be "Continue" most of the time.