{'fraction_non_alphanumeric': 0.033785330129305416, 'fraction_numerical': 0.008439180684288821, 'mean_word_length': 4.468011888002502, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Space Objects Maneuvering Prediction via Maximum Causal Entropy Inverse Reinforcement Learning

Inverse Reinforcement Learning (RL) can be used to determine the behavior of Space Objects (SOs) by estimating the reward function that an SO is using for control. The approach discussed in this work can be used to analyze maneuvering of SOs from observational data. The inverse RL problem is solved using maximum causal entropy. This approach determines the optimal reward function that a SO is using while maneuvering with random disturbances by assuming that the observed trajectories are optimal with respect to the SO's own reward function. Lastly, this paper develops results for scenarios involving Low Earth Orbit (LEO) station-keeping and Geostationary Orbit (GEO) station-keeping.


I. Introduction
Space Situational Awareness (SSA) has many definitions depending on the goal at hand, but in general it involves collecting and maintaining knowledge of all space objects (SOs) orbiting the Earth and the space environment. This task is becoming more difficult as the number of objects currently tracked by the U.S. increases due to breakup events and improving tracking capabilities [1]. The Space Surveillance Network (SSN) is tasked with maintaining information on over 22,000 objects, 1,100 of which are active, with a collection of optical and radar sensors. Determining physically significant characteristics, i.e. attributes, that go beyond simple orbital states is a key objective which is required for protecting space capabilities and achieving SSA. For example, the SSN catalog currently includes radar cross-section and a non-conservative force parameter, analogous to a ballistic coefficient, which provides additional SO characterization information beyond position and velocity. Future SSA systems will have to be capable of building a much more detailed picture of SO attributes in order to maintain better knowledge of their characteristics, which ultimately may lead to better tracking capabilities.
This work discusses the use of inverse Reinforcement Learning (RL) to learn the behavior of Space Objects (SOs) from observed orbital motion. The behavior of SOs is estimated using inverse RL to determine the reward function that each SO is using to control. Since SOs having the capability of maneuver are controlled to achieve a particular mission-driven goal, maneuvering can be very subjective and only a data-driven learning approach can reveal the true goal. It is also important to determine what type of behavior a SO is using and if this behavior changes. Inverse RL approaches use optimal control principles to learn what reward function is being used by an agent given observations. The simplest inverse RL approach, discussed in Ref. [2], solves for the reward function using a weighted sum of features. This is the Feature Matching Approach (FMA). The weights determined from the inverse RL calculation are the representation for the expert reward function. The estimated reward function weights can be used to determine the type of behavior mode the SO is following and to classify the model based on libraries of behavior models. These weight vectors can be added to the state of SOs as a way to represent the policy that the SO is currently following and allow for the change of this policy over time and as the behavior changes. Using the inverse RL approach, the optimal control problem is formed as a Markov Decision Process (MDP) where the reward function is not explicitly given. Rather, observations of expert demonstrations for a given task and the goal is given to estimate the reward function that the expert used to derive the demonstration trajectories. It is common to assume that the expert's actions are optimal with respect to the expert's reward function. Unfortunately, several problems exist with this early approach. First, the formulation is limited to discrete MDP system models which contain discrete state-action pairs at some time. Many problems exist with continuous state-actions (i.e. a low-thrust transfer orbit problem), thus, this formulation is limiting

II. Reinforcement Learning
This section provides a brief introduction to reinforcement learning. Given a discrete-time system model, we can denote the state of the system at time step k by x k . The system dynamics provide the transition from x k to x k+1 given u k , where u k ∈ R denotes the current control action, and this transition may be stochastic. Therefore, it is meaningful to represent this transition with a probability distribution x k+1 ∼ p(x k+1 |x k , u k ) and x k , x k+1 ∈ R n denotes the current and next state, respectively. The actions are modeled probabilistically and are generated by a policy u k ∼ π(u k |x k ) where the randomness in the policy can enable exploration of the policy space while also providing optimality for certain classes of control problems.
An agent (the maneuvering SO) has a current state x k ∈ S (orbital elements) at each discrete time step k and chooses an action u k according to a policy π. For the policy π, a reward signal r k is given for a transition to a new state x k+1 . The general objective of RL is to maximize an expectation over the discounted return, J(θ), given as where γ ∈ [0, 1) is a discount factor and θ are policy parameters. Q-learning is a popular RL method which defines a Q-function that represents the total reward or the total "cost-to-go" for a policy π [11]. Once the Q-function is determined, the action with the highest value or estimated total reward is taken at each time step. Therefore, the optimal policy can be determined using the optimal Q-function. The Q-function of a policy π is where the function estimates the total discounted reward for policy π from state x k assuming that action u k is taken and then all following actions are sampled from policy π. Q-network based methods use neural networks parameterized by θ to represent Q π (x k , u k ; θ), but the dependency notation is dropped for simplicity [11]. Q-networks are optimized by minimizing the following loss function, This equation uses the Bellman optimality condition [11] to relate Q π (x k , u k ) to Q π (x k+1 , u k+1 ), and this equation can be optimized using stochastic gradient descent. Given a reward function, RL can be used to determine a Q-function which is optimal with respect to this reward function.

III. Inverse Reinforcement Learning
The basic principle behind inverse RL is to find the expert's reward function, r k (x k , u k ), that explains the expert's behavior given the observations. Reference [5] introduces a maximum causal entropy method to predicting a reward function using the principle of maximum entropy which will be discussed in the following section.

A. Causal Entropy
As discussed in the introduction, FMA and MMP are formulated by matching feature counts. This is fundamentally ill-posed because no true optimal policy will match the feature counts of a stochastic problem. Many different policies will satisfy the constraint. The principle of maximum entropy determines the least committed probability distribution to stay within the stochastic problem constraints [12]. The causally conditioned probability is an extension that is used when side (state) information, S, are determined after each time-step [13]. The causally conditioned probability of an action, A, on S is where N is the length of time, and the probability of A k is only conditioned on S 1:k and A 1:k−1 . The causal entropy is determined from the causally conditioned probability given by where E P(A,S) is the expectation of the joint probability P(A, S) and the causal entropy, H(A N ||S N ), is the sum of all the entropy for the causally conditioned distribution. Thus, the overall causal entropy measures the uncertainty in the causal conditioned probability distribution. If future information of the distribution is used (e.g. S 1:N ), the additional information decreases the causal entropy because there is less uncertainty present. The maximum causal entropy approach predicts a policy (actions) from P(A k |S 1:k , A 1:k−1 ) based on the state information provided, thus yielding an option to predict a reward function. Note that the joint distribution, P(A, S), can be decomposed to P(A, S) = P(A N ||S N )P(S N ||A N −1 ).

B. Optimization for Maximum Causal Entropy
In order to predict an expected reward function from expert trajectories, an optimization problem is defined by maximizing the causal entropy in Eq. (5). The causal distribution is constrained to match the expert's feature function defined by which is just the sum of features through each time-step k. The empirical expected feature function is defined as TheẼ S,A [F (S, A)] statistics are usually due to limited sampling of the expert. Then, an optimization problem can be formed as given: P(S k ||A k−1 ).
From the constraints in Eq. (8), it is assumed that the state dynamics are provided explicitly with a form Therefore, a prediction on the reward function can be found through optimization of Eq. (8).

C. Application to MDPs and LQR
Both MDPs and LQRs provide methods to represent stochastic systems within an optimal control structure. In inverse RL (or inverse optimal control), the goal is to predict an unknown reward function that a policy or controller is optimized about. Predicting this unknown reward function may give control policies that are close to optimal [14]. As discussed previously, the inverse RL problem contains side information and decisions which are the states and actions respectively. This information is both dependent and stochastic. The relation is provided by the dynamics in Eq. (9).

MDPs
For MDPs, a general reward (or objective) function is given by Eq. (1) which gives a discounted reward for future steps in time. Each reward, r k (x k , u k ), is which is a linear function between the weights, θ, and the feature vector (or counts) f x k ,u k . For an MDP, the optimal policy is determined by finding an action for each state with the highest expected cumulative reward. This can occur with either finite or infinite time horizons [15]. The optimal policy is found by solving the the Bellman equations given by which can be alternatively formulated with a value function, V * (S k ), and a Q-function Q * (S k , A k ) given by Equation (12) is recursively solved using dynamic programming and value iteration to obtain the policy π(S k ) which is similar to the previous Q-learning discussion [15]. For the inverse RL problem, the Q-function and the value function use a softmax function rather than a maximum function as given by Otherwise both equations are equivalent and can be solved similarily. The softmax function allows a smooth differentiable interpolation of the maximum of functions.

LQR
A special case exists using LQR by assuming the dynamics are linear and the stochasticity is due to Gaussian noise.
For the classic LQR problem, the reward function, instead, is defined as a cost function, therefore, they are negatives of each other. For the finite horizon N, the total cost is calculated from initial state x 0 and the control sequence U = [u k , u k+1 , · · · , u N −1 ] applied to the dynamics given by where l(x k , u k ) is the running cost and l f (x N ) is the terminal cost. The LQR costs are quadratic given by where q k , r k , Q k , R k , and P k are the running weights (coefficients) and Q N and q N are the terminal weights. The weight matrices, Q k and R k , are positive definite and the block matrix Q k P k P k R k is positive-semidefinite [16]. The running and terminal costs are substituted into Eq. (14), and due to the symmetry in the weight matrices, the total cost is simplified to The optimal control solution is based on minimizing the cost function in terms of the control sequence which is given by To solve for the optimal control solution given by Eq. (17), a value iteration method is used. Value iteration is a method that determines the optimal cost-to-go (value) starting at the final time-step and moving backwards in time minimizing the control sequence. The cost-to-go is defined as where . This is very similar to Eq. (14), but the only difference is that the cost starts from time-step k instead of k = 0. The optimal cost-to-go is calculated similar to Eq. (17) which is At a time-step k, the optimal cost-to-go function is a quadratic function given by where S k , s k , and c k are computed backwards in time using the value iteration method. First, the final conditions S N = Q N , s N = q N , and c N = c are set. This reduces the minimization of the entire control sequence to just a minimization over a control input at a time-step which is the principle of optimality [17]. To find the optimal cost-to-go, the Riccati equations are used to propagate the final conditions backwards in time given by The Q-function is defined as for LQR [11,18]. By substituting Eq. (15) and (20) into Eq. (22), the Q-function has the form From the state update equation, Using the Ricatti solution from Eq. (21), the optimal control policy is in the affine form where K k is the controller given by and l k is the controller offset given by

D. Maximum Causal Entropy with LQR
The causal conditioned probability (or policy) in Eq. (4) is necessary to determine the causal entropy in Eq. (5) which is to be maximized. For the quadratic cost function weights given in Eq. (15), the value function (Eq. (20)) and Q-function (Eq. (24)) can be found recursively through a value iteration method to determine the causal conditioned probability and causal entropy using which is Gibbs distribution [19]. To predict the LQR weights in Eq. (15) from demonstrated behavior, the optimization problem in Eq. (8) must be solved, but the weights can be predicted through gradient-based techniques using the feature expectation in Eq. (7). The feature expectation gives the likelihood of trajectories using some action distribution. By optimizing the likelihood of predicted trajectories with the demonstrated trajectories under an action distribution, weights will be found that converges to the behavior of the demonstrated trajectories. The gradient is the difference between the expected empirical feature function and the expected predicted feature function defined by which is just the difference between the demonstrated and predicted feature vector for the LQR quadratic cost function.
M is the block weight matrix in Eq. (15). The expectation is easily available from the feature vectors since the stochasticity in the system is Gaussian noise. A general expectation for the quadratic feature vectors given in Eq. (29) is where µ m and Σ m are the mean and covariance of m respectively. Therefore, the block weight matrix, M, can be iteratively optimized by a simple gradient update given by but more effective gradient methods can be used. By recursively updating the Q-function and value function for the estimated trajectories and updating the block weight matrix, M, iteratively through gradient decsent, a cost function that represents the expert's demonstrated trajectory can be predicted.

IV. Dynamical Models
To show the viability of inverse RL for SOs, a orbital model is used to describe dynamics. Perturbations due to aerodynamic drag, solar radiation pressure, and J 2 gravitational disturbances are included in the orbtial model formulation. The nonlinear dynamic equations of motion are linearized about a reference orbit operating point for station-keeping control and are used for the cost function prediction via maximum causal entropy.

A. Orbital Dynamics Model
The dynamics for a SO orbiting Earth is described by the two-body problem is given by where r(t) = [x(t), y(t), z(t)] T is the Cartesian coordinate position vector with respect to the geocentric equatorial frame, r(t) is the position magnitude, µ E is the geocentric gravitational constant, and a d (t) are perturbing accelerations. The perturbing accelerations, a d (t), can be expanded to contain external disturbances given by where a c (t) is control by thrusters, a ad (t) is aerodynamic drag, a J 2 (t) is gravitational perturbations, and a sr p (t) is solar radiation pressure.

Control by Thrust
To control a SO, delta-v maneuvers take place over time. An external force from the thrusters is produced which allows the SO to actively control against any perturbations for station-keeping. The acceleration due to a thrust, T, is where m SO (t) is the mass of the SO, r(t) is the velocity vector, and || · || is the L 2 norm. This shows that the accelerations obtained from a thrust, T(t), must coincide with the direction of the velocity vector. For chemical thrusters, the SO mass decreases due to expenditure of the propellant in the rocket. Therefore, the SO mass decreases in time and can be described by a first order ordinary differential equation, where I sp is the propulsion system's specific impulse, g 0 is the standard gravity at sea-level, and m f (t) is the mass of the propellant. Thus, the total SO mass is m SO (t) = m f (t) + m p , where m p is the payload mass of the SO.

General Dynamical System
The two-body problem and the propellant mass loss rate are dependent, so to propagate the orbital dynamics and mass through time, the state vector and the nonlinear dynamical system are defined as The time-step t is continuous and both x(t) and a d (t) vary with time. For simplicity of notation, the dependency of the time-step, (t), is removed. Solutions of the state through time can be solved with numerical solvers (e.g. a Runge-Kutta 4th order method).

B. Aerodynamic Drag
Aerodynamic drag is a disturbance that affects the life-time of the spacecraft. This drag can disturb the SO from its nominal orbit and possibly deorbit the SO without any correctional maneuver. Since fuel used to correct the trajectory from the disturbance is finite, aerodynamic drag is a disturbance that limits the life-span of the spacecraft. The disturbance becomes more apparent for orbits that are smaller in altitude. Therefore, a SO that is situated in a LEO orbit has more aerodynamic drag than a SO in a GEO orbit. Many models exist that relate the atmospheric density with the altitude on Earth [20]. For this paper, the US Standard Atmosphere 1976 (USSA76) [21] is used to model the aerodynamic drag that occurs on the SO while in orbit. The USSA76 is assumed to be a steady-state 1000km, spherically symmetric gaseous shell surrounding Earth. The variations in the density with altitude is approximated by averaging the year-round conditions at mid-latitudes over many years. Therefore, the model gives realistic densities generally that may not coincide with actual values at a specific point.
To find the disturbance, a ad , the aerodynamic drag force, D ad , acts opposite of the SO velocity relative to the atmosphere,ˆ r r el , given by where r r el is The rotational rate of Earth is ω E , and the direction of the relative velocity is a unit vector The drag is computed as where ρ is the density of the atmosphere, C D is the coefficient of drag, A is the area of the spacecraft. Therefore, the acceleration due to aerodynamic drag is

C. Gravitational Perturbations
It is typically assumed that the Earth has a spherically symmetric mass distribution in which the gravitational field is spherically symmetric about the center of the sphere. In reality, the Earth and other planets are not perfect spheres, but they resemble oblate spheroids instead. Thus, the gravitational field varies with the latitude and radius about Earth [22]. Still, the gravity potential is mostly contributed by the spherically symmetric assumption given by where V ss is the gravity potential due to a spherically symmetric mass distribution. The gravitational potential perturbations due to Earth's oblateness is where R is the equatorial radius of Earth, J k are the Earth zonal harmonics, φ is the polar angle defined by and P k are the Legendre polynomials. The J k are constants that are derived due to orbital motion of SO around Earth. Note that the zonal harmonic for J 1 = 0, thus, Earth has a spherical symmetric mass distribution. The dominant gravitational perturbation from Earth's oblateness is due to the J 2 = 0.00108263 zonal harmonic. Higher order zonal harmonic perturbations can be found in Reference [22]. The gravitational potential includes a Legendre polynomial term obtained by the Rodrigues' formula about an arbitrary variable, w, given by which can be expanded to higher order terms with additional zonal harmonics [22]. By using the most dominant perturbation term, J 2 , only the P 2 term is necessary for the gravitational disturbance formulation, but additional zonal harmonic terms can be added, if necessary. The P 2 derived from the Rodriques' formula is By substituting Eq. (46) into Eq. (43) for k = 2, the gravitation potential perturbation is The acceleration due to Earth's oblateness is defined by which is the negative gradient of Φ(r, φ) in terms of r in the geocentric frame. After taking the gradient, a J 2 is simplified to

D. Solar Radiation Pressure
Solar radiation pressure is a electronmagnetic radiation disturbance from the Sun that affects the SO in orbit. The photosphere (visible surface) of the Sun emits radiation intensity given by where S 0 is the Sun power intensity, R 0 is the radius of the photosphere, and R SE is the distance between the Sun and Earth. Therefore, the solar radiation intensity at Earth is S = 1367W/m 2 for the Sun power intensity, 63.15 × 10 6 W/m 2 . The solar radiation pressure for an Earth-orbiting SO is simply where c is the speed of light. To simplify the geometry of the SO as it is disturbed by solar radiation pressure, it is assumed that the cross-sectional geometry of the SO facing the Sun is a cannonball model (i.e. A sr p = πR 2 SO where R SO is the radius of the SO). The acceleration disturbance is given by where the unit vectorû points to the Sun from the SO. Therefore, a sr p is negative because the solar radiation pressure points away from the Sun. The shadow function, ν, is set to 0 or 1 depending on whether the SO is in Earth's shadow or not, respectively. Reference [23] gives an algorithm to determine whether the SO is in Earth's shadow or not. The radiation pressure coefficient, C R , depends on the color of the spacecraft and lies between 1 and 2. If the cross-sectional area is a black body, all the radiation momentum is absorbed which gives C R = 1. If the color of the cross-sectional area is fully reflective, C R = 2 because all the incoming radiation momentum is reflected which doubles the force on the SO. Equation (52) specifically depends on the A sr p /m SO ratio. Solar sails that have large areas and very little mass are affected more by solar radiation pressure disturbances. To simplify theû calculation, it is assumed that the unit vector points from Earth to the sun instead of the SO to the Sun. This assumption is valid because the angle between the Earth-to-Sun line and the SO-to-Sun line is less than 0.02 • due to the large distance the Sun is away from Earth [23]. The,û defined asû where λ is the solar ecliptic longitude, and the unit vectors,Î andĴ , are part of the geocentric ecliptic frame. A simple rotation matrix involving the angle between Earth's equatorial plan and the ecliptic plane, , iŝ which converts the geocentric ecliptic frame to the geocentric equatorial frame. Thus, the solar radiation pressure in the geocentric equatorial frame is To find the solar ecliptic longitude λ, the obliquity, , and the Julian day since the year 2000 are necessary [24]. This equation is where J D is the Julian day to be simulated. The obliquity, , is determined directly given by To determine the solar ecliptic longitude, both the mean longitude and mean anomaly of the Sun must be found which are from Reference [24]. Thus, the solar ecliptic longitude is

Linearization
The dynamics for a SO orbiting Earth is described by the nonlinear equations of motion given by Eq. (32). The disturbing accelerations (e.g. aerodynamic drag, gravitational perturbations, and solar radiation pressure) can be modeled to high levels of accuracy for time into the future, which allows the nonlinear dynamics to be described as a linear time-varying (LTV) system. The disturbing accelerations are driven through a time-dependent functions discussed previously. By taking a first-order Taylor series expansion of Eq. (32) around an operating point for linearization,x and a c , wherex is the predefined nominal geostationary slot for the SO (i.e.x = x geo ), a linear time-varying system is found given by This equation is simplified to a linear state-space equations using perturbations from the operating point given by where , and ∂a c = (a c −ā c ). The linearized A is time-dependent while the B matrix is time-independent. They are both given by This provides an approximate linearization of the nonlinear dynamics about the geostationary slot operating point.

Discretization
The A and B matrices are discretized using a Runge-Kutta fourth-order with a zero-order hold on the control input a c [25,26]. A general Runge-Kutta fourth-order integration is given by Note that the variable k is the discretized time-step and h is the step size for the system. With the zero-order hold on a control input u, the k 1 , k 2 , k 3 , and k 4 are given by By substituting Eq. (65) into Eq. (64), the discretized system from a Runge-Kutta fourth-order with a zero-order hold on control is The A k is the discrete time-varying version of a continuous time-varying A = A(t). Otherwise, for a continuous time-invariant A(t + 1) = A(t), the discretized version will also be time-invariant (i.e. A k+1 = A k ). Therefore, the linear system approximation given by Eq. (62) can be discretized at each time-step, k, using Eq. (66). Although the continuous time linearized system for the orbital equations of motion is time-dependent in A and time-independent in B in Eq. (63), the discrete time version of the same system contains a time-dependent A k and B k .

F. Conversion of Cartesian Coordinates to Orbital Elements
To describe a SO in orbital motion, either the state vector representation in Eq. (36) or an orbital element representation can be used. The state vector representation describes the position and velocity of SO in orbital motion which is equivalent to describing the shape and orientation of the SO in orbit using the classical orbital elements. The six orbital elements that describe the SO in orbital motion are h, i, Ω, e, ω, and θ which are the specific angular momentum, inclination, right ascension of the ascending node, eccentricity, argument of perigee, and true anomaly, respectively. Although h and θ are used, they can be replaced by the semimajor axis, a, and the mean anomaly, M. The angular momentum, h, and the eccentricity, e, describe the shape of the orbit while the inclination, i, right ascension of the ascending node, Ω, and argument of perigee, ω, describe the orientation of the orbit. The true anomaly, θ, describes the position of the SO on its orbit. More specific descriptions of orbital elements can be found in [23,27]. The advantages of using orbital elements is that it shows the explicit shape and orientation of the orbit compared to knowing the position and velocity using Cartesian coordinates at a time-step. Thus, it is easier to determine whether the SO orbit is on an escape trajectory or orbiting at the equatorial plane using orbital elements.
The transformation of Cartesian coordinates in a geocentric equatorial frame to orbital elements is found in [23]. The angular momentum vector is given by The angular momentum orbital element, h, is just the L 2 norm of Eq. (67). The inclination angle lies in the between 0 • and 180 • and is given by where h Z is the z coordinate of the angular momentum vector h. To obtain the right ascension of the ascending node, the node line must be calculated. The node line is calculated as whereK is the unit vector in the z direction of a Cartesian system. The node magnitude, N, is simply the L 2 norm of Eq.
(69). The right ascension of the ascending node is then which resolves the quadrant ambiguity of the arccosine function from 0 • to 360 • . Next, the eccentricity vector defined as where v r is the radial velocity magnitude given by Again, the eccentricity magnitude e is the L 2 norm. The argument of perigee is considering the quadrant ambiguity of the arccosine function from 0 • to 360 • . Lastly, the true anomaly can be solved using taking into account the arccosine ambiguity from 0 • to 360 • . Two singular cases can arise using orbital elements. When the orbit is circular, e = 0 and ω is undefined. To alleviate this problem, ω is set to 0 when the node magnitude, N = 0. When the orbit lies on the equatorial plane, i = 0 and Ω is undefined. Similar to the other case, Ω is set to 0 when the node magnitude, N = 0. Thus, the six orbital elements are able to alternatively describe the SO orbit based on the position and velocity in Cartesian coordinates based on the geocentric equatorial frame.

V. Simulation Results
This section discusses the initial proof-of-concept results for the proposed maximum causal entropy approach to learning the SO's behavior in orbit. Three cases are considered in this section. The first case discusses the inverse RL approach for GEO station-keeping. The next case applies this same approach to a LEO orbit. The inverse reinforcement learning approach is used to learn the reward function the produces the SO's maneuver. The learned reward function can then be used to estimate the behavior of SOs.

A. GEO Results
The first case involves the simulation of a SO in a GEO orbit. In this scenario, the objective of the SO is to stay in GEO orbit under the influence of solar radiation pressure, gravitational perturbations, and aerodynamic drag. The parameters for the SO are given in Table 1. The ephemeris for a SO in an initial GEO orbit is given in Table 2. To maintain orbital station-keeping, the thrusters are fired every two hours for up to 30 minutes if the SO lies outside a 75km box of the nominal GEO orbit. This is to provide a conventional thruster burn schedule for station-keeping maneuvers [28]. To obtain the expert trajectories, a hundred trajectories were simulated with Gaussian white noise. The cost function for these trajectories were determined using the quadratic cost function in Eq. (15) but with constant weight matrices through time. Figure 1(a) shows a single expert trajectory in GEO orbit with an initial perturbation from the nominal GEO orbit. From this figure, the eccentricity, inclination, and right ascension of the ascending node stay approximately zero through time due to stabilizing LQR control to the nominal GEO orbit. This figure explicitly shows the delta-v burns that occur every two hours when the SO is outside the 75km box around the nominal GEO orbit. A very small transient response in each element can be observed due to the burn itself. This control schedule is able to mitigate the variances due to the disturbances in the system. Figure 2(a) shows the true cost that is used for every expert trajectory. This LQR cost, which is what is considered to be unknown, is estimated using maximum causal entropy inverse RL.  Figure 1(b) shows a trajectory which uses the estimated LQR weights determined by maximum causal entropy. The trajectory is initialized with the same initial conditions as Figure 1(a). Similar to the expert trajectory, the SO is able to mitigate disturbances and maintain a GEO orbit with the same delta-v control schedule. The time-history is similar, but Figure 1(b) has a less distinct transient response from the delta-v burns. The estimated cost function is shown in Figure  2(b). The estimated cost function is able to match the behavior of the expert trajectory. Specifically, the maximum error between the true and estimated cost function is 0.04 as shown in Figure 2(c). Even with these small differences, the estimated cost function captures the general trend of the station-keeping maneuver.

B. LEO Results
A LEO simulation scenario is also considered for learning the cost function used for SO station-keeping using the parameters in Table 1. The SO ephemeris for the initial LEO orbit is provided in Table 3. For the initial circular LEO orbit, the objective is for the SO to stay in LEO under the external disturbances. Similar to the GEO case, the thrusters are fired every two hours for up to 30 minutes to maintain a nominal LEO orbit. One hundred expert trajectories were simulated with Gaussian white noise under constant LQR weight matrices. Figure 3(a) shows a single expert trajectory from an initial perturbation from a LEO orbit. From the figure, the eccentricity, right ascension of the ascending node, and inclination stay approximately constant at 0, 150 • , and 50 • , respectively. The oscillations are due to external disturbances on the SO which are not significant enough to perturb the SO outside a 75km box from the nominal LEO orbit in a 24 hour period. Thus, no thruster burns for SO station-keeping are observed. Figure 4(a) shows the true (and unknown) cost used for every expert trajectory. Since no control thrust is used, this contour does not depend on R.  Figure 3(b) shows a trajectory using the estimated LQR weights found using maximum causal entropy. The estimated trajectory is initialized with the same initial conditions as Figure 3(a). The trajectory that is determined using the estimated weights has an almost identical response to the expert trajectory. Also, the SO trajectory obtained from the estimated weights did not use any thruster burns to maintain its nominal orbit. Thus, it follows the behavior of the expert. The estimated cost function is shown in Figure 4(b). The estimated almost matches the behavior of the expert trajectory. In Figure 4(c), the maximum error is 0.015 between the true and estimated cost functions. Since the trajectories obtain from the expert and estimated weights uses no control within the 24 hour period, the behavior from the estimated weights captures the general trend of the LEO station-keeping maneuver.

VI. Conclusion
The objective of the paper is to formulate inverse RL to learn the behavior of SOs from observed orbital motion using maximum causal entropy and provide examples for a SO in GEO and LEO. By setting up expert SO trajectories using an LQR-based cost, the SO behavior is learned using maximum causal entropy inverse RL. It is shown that the learned cost function in LEO or GEO is comparable to the observed expert SO's trajectory in the presence of orbital disturbances. Thus, maximum causual entropy is attractive in learning the reward function that produces the expert SO's maneuver and can be used to estimate the behavior of SOs.