GEO Station Keeping in GMAT

Question

I am currently trying to implement station keeping maneuvers for a geosynchronous satellite in orbit around Earth using the General Mission Analysis Tool (GMAT) but I am not managing.

I know there are 2 main maneuvers: East-West to correct longitude drift, and North-South to correct inclination but I am not managing to implement these in GMAT.

I have reviewed GMAT's documentation and seen followed the code example for LEO station keeping, which is provided in the samples folder when downloading GMAT, but this deals with altitude station keeping as it is a LEO orbit.

I have also followed an STK tutorial: for maintaining GEO orbit, but am having trouble translating the logic into GMAT's mission section.

Any help or insight provided would be greatly appreciated.

Answer 1

I don't have a specific answer for you, but will offer this up as some sort of consolation prize. GEO stationkeeping is one of the trickiest maneuvers to plan. If you are trying to do it "for real" then you have to factor in the following (in no particular order):

EW perturbations due to the Earth's geopotential are dependent on what longitude you are at. If you haven't already, read up on the stable and unstable points. This will affect the direction of drift and drift rate.

NS perturbations due to solar and lunar gravity.

The size of your longitude box relative to your eccentricity limits and area to mass ratio. A larger area to mass ratio means the daily perturbations due to solar radiation pressure can make it difficult to stay in a tight box.

Eccentricity control strategy. A common method is keeping the eccentricity vector pointed to the sun. There are others. This will affect where in the orbit you do your EW maneuver.

Type of thruster (chemical vs. electric). A low thrust electric thruster has to burn longer and may result in multiple burns a day.

Combined maneuvers - there is nothing (other than potentially your bus hardware constriaints) that says the maneuver has to be pure EW or NS. More advanced spacecraft may combine these along with momentum management thrusting.

Maneuver frequency is generally tied to size of box, type of thrusters, and intangibles like the ops team not wanting to do maneuvers on the weekend.

Attitude control constraints that limit where in the orbit a maneuver may be performed. An older spacecraft I used to operate could not do a maneuver within 2 hours of local noon or midnight due to poorer yaw attitude control in those regions.

This is not meant to discourage you, but to inform you of some of the challenges.

Answer 2

In the end, I managed to implement a station-keeping script for a geosynchronous orbit; it may not represent reality, as most satellites use fixed schedules to maintain orbit, but it was sufficient to solve my problem.

It follows this tutorial for STK closely, with some modifications where necessary. Hopefully this helps anybody trying to implement this in GMAT, at least as a starting point.

The main logic is the following, notably, the use of impulsive burns to correct the orbit:

While elapsed days is under 10 days for example:

if longitude is outside tolerance
    implement East-West station keeping
if inclination is outside tolerance
    implement North-South station keeping

Here is the following GMAT script:

%General Mission Analysis Tool(GMAT) Script


%----------------------------------------
%---------- Spacecraft
%----------------------------------------

Create Spacecraft DefaultSC;
GMAT DefaultSC.DateFormat = UTCGregorian;
GMAT DefaultSC.Epoch = '01 Jan 2000 11:59:28.000';
GMAT DefaultSC.CoordinateSystem = EarthMJ2000Eq;
GMAT DefaultSC.DisplayStateType = Keplerian;
GMAT DefaultSC.SMA = 42164.16900000003;
GMAT DefaultSC.ECC = 1.87903981488134e-16;
GMAT DefaultSC.INC = 10.00000000000001;
GMAT DefaultSC.RAAN = 0;
GMAT DefaultSC.AOP = 0;
GMAT DefaultSC.TA = 0;
GMAT DefaultSC.DryMass = 850;
GMAT DefaultSC.Cd = 2.2;
GMAT DefaultSC.Cr = 1.8;
GMAT DefaultSC.DragArea = 15;
GMAT DefaultSC.SRPArea = 1;
GMAT DefaultSC.SPADDragScaleFactor = 1;
GMAT DefaultSC.SPADSRPScaleFactor = 1;
GMAT DefaultSC.AtmosDensityScaleFactor = 1;
GMAT DefaultSC.ExtendedMassPropertiesModel = 'None';
GMAT DefaultSC.NAIFId = -10001001;
GMAT DefaultSC.NAIFIdReferenceFrame = -9001001;
GMAT DefaultSC.OrbitColor = Red;
GMAT DefaultSC.TargetColor = Teal;
GMAT DefaultSC.OrbitErrorCovariance = [ 1e+70 0 0 0 0 0 ; 0 1e+70 0 0 0 0 ; 0 0 1e+70 0 0 0 ; 0 0 0 1e+70 0 0 ; 0 0 0 0 1e+70 0 ; 0 0 0 0 0 1e+70 ];
GMAT DefaultSC.CdSigma = 1e+70;
GMAT DefaultSC.CrSigma = 1e+70;
GMAT DefaultSC.Id = 'SatId';
GMAT DefaultSC.Attitude = CoordinateSystemFixed;
GMAT DefaultSC.SPADSRPInterpolationMethod = Bilinear;
GMAT DefaultSC.SPADSRPScaleFactorSigma = 1e+70;
GMAT DefaultSC.SPADDragInterpolationMethod = Bilinear;
GMAT DefaultSC.SPADDragScaleFactorSigma = 1e+70;
GMAT DefaultSC.AtmosDensityScaleFactorSigma = 1e+70;
GMAT DefaultSC.ModelFile = 'aura.3ds';
GMAT DefaultSC.ModelOffsetX = 0;
GMAT DefaultSC.ModelOffsetY = 0;
GMAT DefaultSC.ModelOffsetZ = 0;
GMAT DefaultSC.ModelRotationX = 0;
GMAT DefaultSC.ModelRotationY = 0;
GMAT DefaultSC.ModelRotationZ = 0;
GMAT DefaultSC.ModelScale = 1;
GMAT DefaultSC.AttitudeDisplayStateType = 'Quaternion';
GMAT DefaultSC.AttitudeRateDisplayStateType = 'AngularVelocity';
GMAT DefaultSC.AttitudeCoordinateSystem = EarthMJ2000Eq;
GMAT DefaultSC.EulerAngleSequence = '321';


%----------------------------------------
%---------- ForceModels
%----------------------------------------

Create ForceModel DefaultProp_ForceModel;
GMAT DefaultProp_ForceModel.CentralBody = Earth;
GMAT DefaultProp_ForceModel.PrimaryBodies = {Earth};
GMAT DefaultProp_ForceModel.PointMasses = {Luna, Sun};
GMAT DefaultProp_ForceModel.SRP = On;
GMAT DefaultProp_ForceModel.RelativisticCorrection = On;
GMAT DefaultProp_ForceModel.ErrorControl = RSSStep;
GMAT DefaultProp_ForceModel.GravityField.Earth.Degree = 4;
GMAT DefaultProp_ForceModel.GravityField.Earth.Order = 4;
GMAT DefaultProp_ForceModel.GravityField.Earth.StmLimit = 100;
GMAT DefaultProp_ForceModel.GravityField.Earth.PotentialFile = 'JGM2.cof';
GMAT DefaultProp_ForceModel.GravityField.Earth.TideModel = 'None';
GMAT DefaultProp_ForceModel.SRP.Flux = 1367;
GMAT DefaultProp_ForceModel.SRP.SRPModel = Spherical;
GMAT DefaultProp_ForceModel.SRP.Nominal_Sun = 149597870.691;
GMAT DefaultProp_ForceModel.Drag.AtmosphereModel = JacchiaRoberts;
GMAT DefaultProp_ForceModel.Drag.HistoricWeatherSource = 'ConstantFluxAndGeoMag';
GMAT DefaultProp_ForceModel.Drag.PredictedWeatherSource = 'ConstantFluxAndGeoMag';
GMAT DefaultProp_ForceModel.Drag.CSSISpaceWeatherFile = 'SpaceWeather-All-v1.2.txt';
GMAT DefaultProp_ForceModel.Drag.SchattenFile = 'SchattenPredict.txt';
GMAT DefaultProp_ForceModel.Drag.F107 = 150;
GMAT DefaultProp_ForceModel.Drag.F107A = 150;
GMAT DefaultProp_ForceModel.Drag.MagneticIndex = 3;
GMAT DefaultProp_ForceModel.Drag.SchattenErrorModel = 'Nominal';
GMAT DefaultProp_ForceModel.Drag.SchattenTimingModel = 'NominalCycle';
GMAT DefaultProp_ForceModel.Drag.DragModel = 'Spherical';

%----------------------------------------
%---------- Propagators
%----------------------------------------

Create Propagator DefaultProp;
GMAT DefaultProp.FM = DefaultProp_ForceModel;
GMAT DefaultProp.Type = RungeKutta89;
GMAT DefaultProp.InitialStepSize = 60;
GMAT DefaultProp.Accuracy = 9.999999999999999e-12;
GMAT DefaultProp.MinStep = 0.001;
GMAT DefaultProp.MaxStep = 2700;
GMAT DefaultProp.MaxStepAttempts = 50;
GMAT DefaultProp.StopIfAccuracyIsViolated = true;

%----------------------------------------
%---------- Burns
%----------------------------------------

Create ImpulsiveBurn EWBurn;
GMAT EWBurn.CoordinateSystem = Local;
GMAT EWBurn.Origin = Earth;
GMAT EWBurn.Axes = VNB;
GMAT EWBurn.Element1 = 0;
GMAT EWBurn.Element2 = 0;
GMAT EWBurn.Element3 = 0;
GMAT EWBurn.DecrementMass = false;
GMAT EWBurn.Isp = 300;
GMAT EWBurn.GravitationalAccel = 9.81;

Create ImpulsiveBurn NSBurn1;
GMAT NSBurn1.CoordinateSystem = Local;
GMAT NSBurn1.Origin = Earth;
GMAT NSBurn1.Axes = VNB;
GMAT NSBurn1.Element1 = 0;
GMAT NSBurn1.Element2 = 0;
GMAT NSBurn1.Element3 = 0;
GMAT NSBurn1.DecrementMass = false;
GMAT NSBurn1.Isp = 300;
GMAT NSBurn1.GravitationalAccel = 9.81;

Create ImpulsiveBurn NSBurn2;
GMAT NSBurn2.CoordinateSystem = Local;
GMAT NSBurn2.Origin = Earth;
GMAT NSBurn2.Axes = VNB;
GMAT NSBurn2.Element1 = 0;
GMAT NSBurn2.Element2 = 0;
GMAT NSBurn2.Element3 = 0;
GMAT NSBurn2.DecrementMass = false;
GMAT NSBurn2.Isp = 300;
GMAT NSBurn2.GravitationalAccel = 9.81;

%----------------------------------------
%---------- Solvers
%----------------------------------------

Create DifferentialCorrector DefaultDC;
GMAT DefaultDC.ShowProgress = true;
GMAT DefaultDC.ReportStyle = Normal;
GMAT DefaultDC.ReportFile = 'DifferentialCorrectorDefaultDC.data';
GMAT DefaultDC.MaximumIterations = 100;
GMAT DefaultDC.DerivativeMethod = ForwardDifference;
GMAT DefaultDC.Algorithm = NewtonRaphson;

%----------------------------------------
%---------- Subscribers
%----------------------------------------

Create OrbitView DefaultOrbitView;
GMAT DefaultOrbitView.SolverIterations = Current;
GMAT DefaultOrbitView.UpperLeft = [ 0.001354554690145615 0 ];
GMAT DefaultOrbitView.Size = [ 0.2163901117507619 0.5598047192839707 ];
GMAT DefaultOrbitView.RelativeZOrder = 381;
GMAT DefaultOrbitView.Maximized = false;
GMAT DefaultOrbitView.Add = {DefaultSC, Earth};
GMAT DefaultOrbitView.CoordinateSystem = EarthMJ2000Eq;
GMAT DefaultOrbitView.DrawObject = [ true true ];
GMAT DefaultOrbitView.DataCollectFrequency = 1;
GMAT DefaultOrbitView.UpdatePlotFrequency = 50;
GMAT DefaultOrbitView.NumPointsToRedraw = 0;
GMAT DefaultOrbitView.ShowPlot = true;
GMAT DefaultOrbitView.MaxPlotPoints = 20000;
GMAT DefaultOrbitView.ShowLabels = true;
GMAT DefaultOrbitView.ViewPointReference = Earth;
GMAT DefaultOrbitView.ViewPointVector = [ 30000 0 0 ];
GMAT DefaultOrbitView.ViewDirection = Earth;
GMAT DefaultOrbitView.ViewScaleFactor = 1;
GMAT DefaultOrbitView.ViewUpCoordinateSystem = EarthMJ2000Eq;
GMAT DefaultOrbitView.ViewUpAxis = Z;
GMAT DefaultOrbitView.EclipticPlane = Off;
GMAT DefaultOrbitView.XYPlane = On;
GMAT DefaultOrbitView.WireFrame = Off;
GMAT DefaultOrbitView.Axes = On;
GMAT DefaultOrbitView.Grid = Off;
GMAT DefaultOrbitView.SunLine = Off;
GMAT DefaultOrbitView.UseInitialView = On;
GMAT DefaultOrbitView.StarCount = 7000;
GMAT DefaultOrbitView.EnableStars = On;
GMAT DefaultOrbitView.EnableConstellations = On;

Create GroundTrackPlot DefaultGroundTrackPlot;
GMAT DefaultGroundTrackPlot.SolverIterations = Current;
GMAT DefaultGroundTrackPlot.UpperLeft = [ 0.001354554690145615 0.563873067534581 ];
GMAT DefaultGroundTrackPlot.Size = [ 0.2163901117507619 0.5598047192839707 ];
GMAT DefaultGroundTrackPlot.RelativeZOrder = 377;
GMAT DefaultGroundTrackPlot.Maximized = false;
GMAT DefaultGroundTrackPlot.Add = {DefaultSC};
GMAT DefaultGroundTrackPlot.DataCollectFrequency = 1;
GMAT DefaultGroundTrackPlot.UpdatePlotFrequency = 50;
GMAT DefaultGroundTrackPlot.NumPointsToRedraw = 0;
GMAT DefaultGroundTrackPlot.ShowPlot = true;
GMAT DefaultGroundTrackPlot.MaxPlotPoints = 20000;
GMAT DefaultGroundTrackPlot.CentralBody = Earth;
GMAT DefaultGroundTrackPlot.TextureMap = 'ModifiedBlueMarble.jpg';

Create XYPlot XYPlot1;
GMAT XYPlot1.SolverIterations = Current;
GMAT XYPlot1.UpperLeft = [ 0.4118602761982129 0.05240912933220625 ];
GMAT XYPlot1.Size = [ 0.5004061738424046 0.4513947590870668 ];
GMAT XYPlot1.RelativeZOrder = 236;
GMAT XYPlot1.Maximized = false;
GMAT XYPlot1.XVariable = DefaultSC.EarthMJ2000Eq.ModEquinoctialF;
GMAT XYPlot1.YVariables = {DefaultSC.EarthMJ2000Eq.ModEquinoctialG};
GMAT XYPlot1.ShowGrid = true;
GMAT XYPlot1.ShowPlot = true;

Create ReportFile ReportFile1;
GMAT ReportFile1.SolverIterations = Current;
GMAT ReportFile1.UpperLeft = [ 0.09989840839823907 0.1253051261187958 ];
GMAT ReportFile1.Size = [ 0.717236708432103 0.8169243287225386 ];
GMAT ReportFile1.RelativeZOrder = 395;
GMAT ReportFile1.Maximized = false;
GMAT ReportFile1.Filename = 'ReportFile1.txt';
GMAT ReportFile1.Precision = 16;
GMAT ReportFile1.Add = {DefaultSC.ElapsedDays, DefaultSC.EarthMJ2000Eq.RAAN, DefaultSC.EarthMJ2000Eq.RA, DefaultSC.EarthFixed.PlanetodeticLON};
GMAT ReportFile1.WriteHeaders = true;
GMAT ReportFile1.LeftJustify = On;
GMAT ReportFile1.ZeroFill = Off;
GMAT ReportFile1.FixedWidth = true;
GMAT ReportFile1.Delimiter = ' ';
GMAT ReportFile1.ColumnWidth = 23;
GMAT ReportFile1.WriteReport = true;

Create XYPlot XYPlot2;
GMAT XYPlot2.SolverIterations = Current;
GMAT XYPlot2.UpperLeft = [ 0.1537419573315273 0.1993490642799024 ];
GMAT XYPlot2.Size = [ 0.2163901117507619 0.4328722538649308 ];
GMAT XYPlot2.RelativeZOrder = 373;
GMAT XYPlot2.Maximized = false;
GMAT XYPlot2.XVariable = DefaultSC.ElapsedDays;
GMAT XYPlot2.YVariables = {DefaultSC.EarthFixed.PlanetodeticLON};
GMAT XYPlot2.ShowGrid = true;
GMAT XYPlot2.ShowPlot = true;

Create XYPlot XYPlot3;
GMAT XYPlot3.SolverIterations = Current;
GMAT XYPlot3.UpperLeft = [ 0.6525567219776498 0.1334418226200163 ];
GMAT XYPlot3.Size = [ 0.2163901117507619 0.4328722538649308 ];
GMAT XYPlot3.RelativeZOrder = 369;
GMAT XYPlot3.Maximized = false;
GMAT XYPlot3.XVariable = DefaultSC.ElapsedDays;
GMAT XYPlot3.YVariables = {EWBurn.Element1};
GMAT XYPlot3.ShowGrid = true;
GMAT XYPlot3.ShowPlot = true;

Create XYPlot XYPlot4;
GMAT XYPlot4.SolverIterations = Current;
GMAT XYPlot4.UpperLeft = [ 0.3863867253640366 0.4938974776240846 ];
GMAT XYPlot4.Size = [ 0.2086014222824246 0.4328722538649308 ];
GMAT XYPlot4.RelativeZOrder = 385;
GMAT XYPlot4.Maximized = false;
GMAT XYPlot4.XVariable = DefaultSC.ElapsedDays;
GMAT XYPlot4.YVariables = {DefaultSC.EarthMJ2000Eq.INC};
GMAT XYPlot4.ShowGrid = true;
GMAT XYPlot4.ShowPlot = true;

Create XYPlot XYPlot5;
GMAT XYPlot5.SolverIterations = Current;
GMAT XYPlot5.UpperLeft = [ 0.1137286758732738 0.5545224006762468 ];
GMAT XYPlot5.Size = [ 0.5004061738424046 0.4497041420118343 ];
GMAT XYPlot5.RelativeZOrder = 280;
GMAT XYPlot5.Maximized = false;
GMAT XYPlot5.XVariable = DefaultSC.ElapsedDays;
GMAT XYPlot5.YVariables = {NSBurn1.Element1, NSBurn2.Element2, EWBurn.Element1};
GMAT XYPlot5.ShowGrid = true;
GMAT XYPlot5.ShowPlot = true;


%----------------------------------------
%---------- Mission Sequence
%----------------------------------------

BeginMissionSequence;
While 'While ElapsedDays < Threshold' DefaultSC.ElapsedDays < 365
   Propagate 'Prop One Step' DefaultProp(DefaultSC);
   If 'If Lon < | > Tol' DefaultSC.EarthFixed.PlanetodeticLON >= 80 | DefaultSC.EarthFixed.PlanetodeticLON <= 79
      Propagate 'Prop To Apogee' DefaultProp(DefaultSC) {DefaultSC.Earth.Apoapsis};
      Target 'EW Station Keeping' DefaultDC {SolveMode = Solve, ExitMode = DiscardAndContinue, ShowProgressWindow = true};
         Vary 'Vary EWBurn.Element1' DefaultDC(EWBurn.Element1 = 0.0001, {Perturbation = 0.00001, Lower = -1000, Upper = 10000, MaxStep = 0.001, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
         Maneuver 'Apply EWBurn' EWBurn(DefaultSC);
         Propagate 'Prop To Asc. Node' DefaultProp(DefaultSC) {DefaultSC.EarthMJ2000Eq.RA = 0};
         Propagate 'Prop To Asc. Node' DefaultProp(DefaultSC) {DefaultSC.EarthMJ2000Eq.RA = 0};
         Achieve 'Acheive Lon' DefaultDC(DefaultSC.EarthFixed.PlanetodeticLON = 79.67, {Tolerance = 0.001});
      EndTarget;  % For targeter DefaultDC
      Propagate DefaultProp(DefaultSC) {DefaultSC.EarthMJ2000Eq.RA = 0};
   EndIf;
   If 'If Inc < | > Tol' DefaultSC.EarthMJ2000Eq.INC >= 10.05 | DefaultSC.EarthMJ2000Eq.INC <= 9.95
      Propagate 'Prop To Asc. Node' DefaultProp(DefaultSC) {DefaultSC.EarthMJ2000Eq.RA = 0};
      Target 'NS Station Keeping' DefaultDC {SolveMode = Solve, ExitMode = DiscardAndContinue, ShowProgressWindow = true};
         Vary 'Vary NSBurn1.Element1' DefaultDC(NSBurn1.Element1 = 0.001, {Perturbation = 0.0001, Lower = -10, Upper = 10, MaxStep = 0.010, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
         Maneuver 'Apply NSBurn1' NSBurn1(DefaultSC);
         Propagate 'Prop To Asc. Node' DefaultProp(DefaultSC) {DefaultSC.EarthMJ2000Eq.RA = 0};
         Achieve 'Acheive SMA' DefaultDC(DefaultSC.Earth.SMA = 42164.169, {Tolerance = 0.1});
         Vary 'Vary NSBurn2.Element2' DefaultDC(NSBurn2.Element2 = -0.2, {Perturbation = 0.0001, Lower = -1, Upper = 1, MaxStep = 0.010, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
         Maneuver 'Apply NSBurn2' NSBurn2(DefaultSC);
         Achieve 'Acheive Inc' DefaultDC(DefaultSC.EarthMJ2000Eq.INC = 10, {Tolerance = 0.1});
      EndTarget;  % For targeter DefaultDC
   EndIf;
EndWhile;