Appunti Space missions and systems

Space Missions and Systems Appunti

Martina Santoianni

May 2023

1 Orbit Determination

Each body constitutes a guided by certain laws that define its state and predict its future

dynamical system

evolution. For spacecrafts the minimum variables requires are position and velocity at a certain time. Usually

we use angles, range and range-rate to have more information about the since it is not

dynamical model,

directly measurable. As we imagine, the predicted state differs from the real state if there are

- errors of inadequacy

- measurement noise

- errors in numerical procedure

As a consequence the state determination must be repeated, it is impossible to propagate the state of a

spacecraft for very long time without incurring in serious discrepancies. The time interval between updates

depends on the accuracy of the dynamical model, the quality of observations and the user needs.

These methods can be applied to every dynamical system. They are pillars of every control system:

where the flight dynamics are the dynamical model used, the sensors are those that provide the range

of observable quantities, the observer that is dynamical model + observable quantities. These three give an

estimate of the state.

There may be external torques that move the system away from the desired state then we need to

x (t),

d

compare it with hence the Since they’re usually in disagreement, they produce an error

x̂(t), computed state.

and we want it to be zero. In this situation, the tells us what is necessary to do

−

e(t) = x̂(t) x (t) regulator

d

to put the s/c back on the desired trajectory.

The regulator controls the that make the maneuver effective (like thrusters).

actuators

After this correction, the state is measured again to verify that the s/c is in the desired state (corrections

usually made to enter the sphere of influence of a planet or escape from it - "stocastic" maneuver, loaded on

board).

The is a preloaded controller (thrust) that acts with the regulator that leaves small ∆v,

feed-forward

coming from the preprogrammed thuster, in the feed forward.

Other basic concepts are:

• a system is controllable if it can be brought from a state to a state in a finite

X X

Controllability: 1 2

time by applying appropriate control inputs (forces).

u 1

• a system is observable if there exists a set of observations that uniquely determines the

y

Observability:

state.

Of course not all systems are controllable or observable.

Figure 1: example: 3 trajectories

Here we consider three trajectories:

Example: 3 trajectories

- the true

- estimated by observer (ground station)

- nominal (best guess from linearisation we have before next observation)

The equations of motion and their solution:

Figure 2: equations

Since the initial state is not directly accessible, we can relate the state at time by writing the equations

t j

in matrix form:

The previous equations can be rewritten in compact form as: (1)

−

J(X ) = Y G(X , t) = 0

0 0

2

Figure 3: eqs in matrix form

Figure 4: observables relation with state

where are the observed observables and are the computed observables that depend on (non-linear

Y G X

0

functions of the state).

The equations can be solved iteratively using Newton-Raphson method, the process is stopped when the

vector is sufficiently small n+1

||X −

J X n|| < ϵ.

0

0

• In real orbit determination problems, the dependence of the state at time t on the initial condition is

almost non-linear.

• The equations leading to the solution of the OD problem are non-linear and must be solved by means of

an iterative procedure, initialised at a best guess initial state .

X

0

• The non-linearity of the equations implies that multiple solutions are possible. The closer is to the

X

0

true state, the smaller the risk of reaching a wrong state and the faster is the convergence.

1.1 Linearisation procedure

We start by writing the equations in terms of deviations from an initial state. State vector:

3 



X

Y 

 (2)





Ẋ

X = 

 

 Ẏ 

 g

Figure 5: linearisation

where is the and is called We obtain a linearised

∗

X X

reference trajectory nominal trajectory.

equation of motion around the starred trajectory.

We also need to linearise the expressions of the observables with respect to the reference trajectory, they

change in this way: Figure 6: linearisation of observables

The linearised equations for observation deviation can be compactly written like: (3)

y = H̃x + ϵ

Figure 7: Linearised eqs compact

For the state deviation we have: (4)

ẋ(t) = A(t)x(t)

4

Figure 8: Linear estimation problem

We need to correct the initial guess in order to have observable quantities that can be confronted with

X

0

our trajectory. So we vary the initial state around the reference trajectory to apply the iterative procedure and

to have a better match between computed and observed observables.

The corrected guess for the initial state is X + δx

0

Figure 9: Procedure

is a matrix of all components of the force derived with respect to all the components of the state vector.

∂F

∂x(t)

is called Dynamical matrix.

A

1.2 State transition matrix

This matrix is solution to the linearised equations of motion. (5)

x(t) = Φ(t, x )x

0 0

The matrix is called "mapping matrix" and maps the state deviation into a deviation of the observables

H(t)

(partial derivatives of observables with respect to initial time). It’s a long and thin matrix since it is with

nxm

(position and velocity of s/c) and number of observables.

n = 6 m

So by changing the IC we obtain the change of computed observables at any later time.

Properties of STM:

1.3 OD problem

given an epoch the state propagation equations and the observation-state relation is:

t

Statement: k (6)

x(t) = Φ(t, t )x

k k (7)

y = Hx + ϵ

k

5

Figure 10: state transition matrix

Figure 11: Properties and equations we use

find the best estimate of .

x k

In general the number of observations largely exceeds the number of state variables, so that the system is

always undetermined (m knowns, the observations and unknowns, the errors and state variables).

m + n

The best solution to this problem is the one that minimises the errors, the cost function;

X (8)

2

ϵ i

we take the square of the residuals and the product must be minimum, close to zero.

T

ϵ ϵ = J(x )

k

This problem is called since we have to find the minimum of the sum of the square

least square problem

of the residuals.

We can always compute (8) and check if after the minimisation, the error is similar to the previous iterations,

this may mean that the quality of the observables may be compromised.

We obtain the residuals like: (9)

−

ϵ = y Hx k

6

that is "observed observables - computed observables".

we compute: 1 1 (10)

T T

− −

J(x) = ϵ ϵ = (y Hx) (y Hx)

2 2

∂J(x)

and minimise. So we compute the derivative .

∂x

Figure 12: Demonstration of least square solution

This is the best estimate of in the least squares sense. The solution exists only if hence

T ̸

x det(H H) = 0,

k

if the matrix is non-singular and invertible.

The solution is a combination of the linear observations and the mapping matrix. Once we find the solution,

we have a new trajectory which becomes the new reference trajectory (restart the process to minimise the error).

1.4 Estimation

This is the process of determining the best estimate of the parameters of the models:

• State (position and velocity) of s/c at a certain time

• State, mass and gravity field

• Magnitude and direction of a non-gravitational force (ex: SRP), useful to know if it is necessary to change

orbit

• Accuracy of a measurement type (errors models)

Generally, we may need to iterate more that once to have convergence, especially if the problem is strongly

non-linear.

Fundamentals of orbit determination:

• For s/c orbit estimation, a linear least square approach is used.

- Linearisation: problems are non-linear and require linearisation to determine a closed form solution

- Iteration: required to reach acceptable solution

- Minimisation of cost function, usually quadratic

- Use Gaussian statistics when possible (normal distribution)

- Monte Carlo simulations when statistics is more complicated

For Gaussian statistics we use the normal distribution:

1 (x−x)2

− (11)

√

p(x) = e 2

2σ

2πσ

The ground antenna provides range of measurements of the satellite.

Example: Motion in a central field 7

(a) (b)

b b

Figure 13: Motion in a central field

1.5 STM and Lambert problem

Pseudo-Lambert problem: given a s/c at at time , find the orbit that places it in at time . (Problem

r t r t

1 1 2 2

often encountered for rendez-vous). Figure 14: Lambert problem

Our goal is to determine the velocity variation at that takes us as close as possible to at .

δv r r t

1 1 2 2

This problem has a closed form solution for a two-body system. To deal with more complicated situations,

we can consider a linearised dynamics about the classical solution (that gives a that doesn’t bring the s/c

v 1

in the desired position due to many perturbing actions). We start by writing the propagation equation of the

linearised dynamics, with state deviation:

x

δr (12)

x = δv ∂r ∂r

δr (13)

1

∂r ∂v

x(t) = Φ(t, t )x(t ) = 1 1

1 1 ∂v ∂v δv

1

∂r ∂v

1 1

where last vector is the vector of initial conditions, by changing we make the s/c closer to the desired

x(t )

1

trajectory.

We compute the state using these equations and the solution leads to an error This error can be

x(t ) e.

2

minimised iteratively by solving the equations: ∂r

∂r

−e 0 (14)

∂r ∂v

= 1 1

∂v ∂v

δv δv

1

2 ∂r ∂v

1 1

∂r ∂r (15)

−1

−e −(

= δv => δv = ) e

1 1

∂v ∂v

1 1

We null the vector by imposing the equal to We stop the iteration when the norm of the error is

δv

e -e.

1

sufficiently small for our needs. 8

1.6 Lunar return trajectory

Orbit the guarantees return to Earth should the propulsion system fail (free return trajectory).

Figure 15: Trajectory sketch

In order to land safely the s/c must reach Earth at an altitude of about 100 km with an angle −7.5.

γ =

Since we control velocity we can set additional condition: the landing occurs near recovery fleet, located at

v 1

angle from the entry point.

θ Figure 16: Conditions

Figure 17: Calculation

1.7 Least squares solution

This solution is a random vector that depends on the noise realisation, we want to investigate the statistics of

this solution calculated with the ensemble average that gives contribution of multiple experiments performed

x̂,

with different instruments. 9

time and ensemble average are equal.

Ergodic hypothesis: Figure 18: Statistical interpretation

Taking the average leads to −1

T T

< x̂ >= x + (H H) H < ϵ >

that, supposing it is the true state, the state is this happens if the noise is (flattened

ϵ

unbiased, zero mean

residuals): (the estimate is in average equal to the true state).

< x̂ >= x

The measurement most prone to bias is the station can easily show a bias of 5 m circa, that is a

range,

systematic error, not random. If there’s a systematic error: (16)

< x̂ ≯ = x (17)

̸

< ϵ >= ϵ̄ = 0 (18)

−1

T T

< x̂ >= x + (H H) H ϵ̄

doesn’t coincide with the true state anymore. All solutions correspond to different realisation of will be

ϵ,

around the point and not around the true state (it is useful to estimate the state bias). So errors can be:

< x̂ >

- Random: normal distribution

- Systematic: all measurements show the same error

It is important to understand what is the uncertainty in the LS estimate:

assume where is the of the measurement errors.

T T

< ϵ >= 0, < ϵϵ >= I, ϵϵ Covariance Matrix

−1

T T

−

x̂ x = (H H) H ϵ

T T

−1 −1 −1 −1

T T T T T T T T T

− −

(x̂ x)(x̂ x) = (H H) H ϵ (H H) H ϵ = (H H) H ϵ ϵ H (H H)

T

−1 −1

T T T T T

− −

< (x̂ x)(x̂ x) >= (H H) H < ϵϵ > H (H H) =

T

−1 −1

T T

(H H) = (H H)

is a symmetric matrix called , the inverse, is the

−1

T T

H H P = (H H)

normal matrix. covariance matrix:

its elements provide the variance-covariance of the estimated parameters: the diagonal elements are the variances

of the state vector elements, the off-diagonal terms are their covariances.

The variance is the probability density of variables (probability of finding a random variable in an infinites-

imal interval Using the Gauss formulation, the variance defines the width of the "bell" function.

x + dx.

x

Having as state vector , random vector, the variances of position and velocity are:

x = v 2 2

−

< (x̂ x) >= σ x

2 2

−

< (v̂ v) >= σ v

10

We define as the average of product of variables:

covariance − −

< (x̂ x)(v̂ v) >= σ

xv

From the second assumption we notice that e terms on diagonal while the

t 21 2

< ϵϵ >= I < ϵ >= σ = 1,

1

terms off diagonal show no correlation between a measurement and another done at a later time. The covariance

matrix of the solution is :

−1

T

(H H) 2

σ

σ xv

x

P = 2

σ σ

xv v

To sum up:

• Cost function, sum of squared residuals: 1 1

1 X

T 2 T

− −

ϵϵ = ϵ = (y Hx) (y Hx)

J(x) = i

2 2 2 2

∂J ∂ J

T > 0

=0 2

∂x ∂x

• minimise cost function, since it’s non-negative it must have a minimum

∂J T T

−(y − −H −

= 0 = HX) H = (y HX)

∂x T T

(H H)x̂ = H y

• Least squares solution −1

T T

x̂ = (H H) H y

• Residuals: −

ϵ = y H x̂

Now we have computed the average, but we want to know also the of the root mean square of

variance

each component of the state vector, that says how much confidence we can have in the solution.

We have defined the covariance matrix as: ¯ˆ ¯ˆ ¯ˆ

2 2

σ σ < x > < x v >

xv

x

P = = ¯ˆ ¯ˆ ¯ˆ

2 2

σ σ < v x > < v

xv v ¯ˆ

the ensemble average of position and velocity errors, found by (estimate - true) with estimator

−

x = x̂ x

variance .

2 2

−

< (x̂ x) >= σ x

1.8 Weighted least-squares solution

Usually not all data are equally accurate, since tracking may occur from different stations with different char-

acteristics or some data may be noisier due to weather conditions, so we need to give more weight to good data.

In this solutions, each residual is given a weight inversely proportional to its standard deviation.

w i −1

T T T T

(H W H)x = H W y x

ˆ = (H W H) H W y

k k

is called we assume that it’s diagonal, but this assumption may be dropped if observa-

W weight matrix,

tions at different times are correlated.

So we write: 2

ϵ

X i

2 2

−

< (x̂ x) >= σ =⇒ J(x) =

x 2

σ i

11

Figure 19: Definition of weights in the LS

where the elements in the sum are pure numbers, since they’re divided by the variance (that gives more

weight to good data, as → ∞ ⇒ →

σ J 0).

i , and we can define also a

Hence the weight is given by 1 = w covariance weighted matrix:

i

2

σ i  

1 .. 0

2

σ 1 (19)

T 1 .. 0

J(x) = ϵ wϵ =⇒ W =  

2

σ

 

2

0 ... ...

by deriving: ∂J T

−H −

= W (y Hx) = 0

∂x

k

T T

W y + H W Hx = 0

H −1

T T T T

(H W H)x = H W y =⇒ x̂ = (H W H) H W y

which is the weighted least squares solution.

In this case WLS the matrix to be inverted to find the covariance is instead of in the

T T

(H W H) (H H)

LS. If is not full-rank, the WLS solution doesn’t exist, this means that some parameters may not be

T

H W H

observable (case of OD problems where the normal matrix is close to singularity).

1.9 Weighted LS with a priori information

In most OD problems a priori information is available on the state vector. The knowledge of the initial estimate

and its covariance may be used together with the data to obtain a better estimate of the state vector. We

minimise the augmented functional:

1 1

T T

− − − −

J(x ) = (y Hx) (y Hx) + ( x

¯ x ) W̄ (

x

¯ x )

k k k k k

k 2 2

where is the and the last term comes from additional observations, is

W̄ x

a priori covariance matrix

k k

the a priori unbiased value of state vector, is the covariance matrix of residuals (diagonal).

−1

W

We can show that the new WLS solution with a priori information is (20)

−1

T T

x

ˆ = (H W H + W̄ ) (H W H + W̄ x