Hi, this question may be somewhat answered already in the forums, but the topic is not completely clear to me so I am asking nevertheless. I am working in an Orbit Determination project and initially was using the BatchLeast squares estimator of orekit directly. Due to some convergence problems and the need to use very specific step selection algorithms I ended up implementing my own estimator. It is quite complex (I actually need to simplify it in some places), and it even has its own observation model implementations that account for light time travel with derivatives computation. Even though it seems to work fine, at the end of the iterations, when it is already close to meeting the criteria for convergence, there is strong numerical instability. My guess is that in that region the problem behaves very noisily, due to the nature of the numerical propagator precision limits and maybe even problems with my implementation of the observation models. I am trying to make a modification so that all the physical part of the problem computation (measurement residuals + derivatives) is handled by Orekit, and I want to do it using the tools Orekit provides as much as possible.
I want to use the âgetParameterDerivativesâ method of the EstimatedMeasurement object I compute from my list of ObservedMeasurement (which in synthetic scenarios can be generated, but in general would be obtained from parsing files with measurements that come from actual sensors). This EstimatedMeasurement is computed at the current state estimate (at the instant of the measurement) of a given iteration, which in turns comes from the initial orbit state being estimated. I am using a propagator builder to activate the Orbital parameters being estimated (and extra parameters if necesary) in the hopes that derivative information is carried on up to the measurement estimation and that âgetParameterDerivativesâ comes with the information I need for each measurement in the batch.
That was my hope, but it is throwing errors that indicate those derivatives wrt the orbital parameters are not available at the EstimatedMeasurement. I am surelly doing something wrong, or this is not the way I am suppose to get derivatives of measurements.
partials = estimated.getParameterDerivatives(driver)
Traceback (most recent call last):
File "<string>", line 1, in <module>
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
Java stacktrace:
org.orekit.errors.OrekitIllegalArgumentException: nombre de parĂĄmetro no conocido Px = -618947.0541776249, parĂĄmetro conocido <none>
at org.orekit.estimation.measurements.EstimatedMeasurement.getParameterDerivatives(EstimatedMeasurement.java:147)
This code snippet is a test, which gives the previous error in the call to âgetParameterDerivativesâ:
def solve_least_squares_orekit_way(
self,
max_iterations=10
) -> bool:
"""
Least squares using Orekit's measurement models and parameter drivers.
"""
if self._builder is None:
self._initialize_builder()
# Get parameter drivers
orbital_drivers = self._builder.getOrbitalParametersDrivers()
prop_drivers = self._builder.getPropagationParametersDrivers()
# Select which parameters to estimate
for driver in orbital_drivers.getDrivers():
driver.setSelected(True) # Always estimate orbital state
# Select propagation parameters if configured
if self._config.solver_options.doCD:
cd_driver = prop_drivers.findByName("drag coefficient")
if cd_driver:
cd_driver.setSelected(True)
if self._config.solver_options.doCR:
cr_driver = prop_drivers.findByName("reflection coefficient")
if cr_driver:
cr_driver.setSelected(True)
# Get all selected drivers in order
all_drivers = []
for driver in orbital_drivers.getDrivers():
if driver.isSelected():
all_drivers.append(driver)
for driver in prop_drivers.getDrivers():
if driver.isSelected():
all_drivers.append(driver)
n_params = len(all_drivers)
measurements, sensor_configs = zip(*self._measurements)
for iteration in range(max_iterations):
# Build propagator with current normalized parameter values
propagator = self._builder.buildPropagator(
self._builder.getSelectedNormalizedParameters()
)
Y_list = [] # Residuals
H_list = [] # Jacobians
W_list = [] # Weights (1/sigma^2)
for meas, sensor in zip(measurements, sensor_configs):
# Propagate to measurement epoch
state = propagator.propagate(meas.getDate())
# Use Orekit's measurement model to get prediction and derivatives
estimated = meas.estimate(0, 0, [state])
# Get predicted and observed values
y_pred = np.array(estimated.getEstimatedValue())
y_obs = np.array(meas.getObservedValue())
n_components = len(y_obs)
# Convert to km for numerical stability (optional)
if meas.measurementType in [Range.MEASUREMENT_TYPE, RangeRate.MEASUREMENT_TYPE]:
y_pred = y_pred / 1000.0
y_obs = y_obs / 1000.0
# Compute residual (observed - predicted)
residual = y_obs - y_pred
Y_list.append(residual)
# Build Jacobian matrix for this measurement
# Shape: [n_components, n_params]
H_meas = np.zeros((n_components, n_params))
for param_idx, driver in enumerate(all_drivers):
# Get partial derivatives w.r.t. this parameter
# This returns an array with one element per measurement component
partials = estimated.getParameterDerivatives(driver)
# Fill the corresponding column in the Jacobian
for comp_idx in range(n_components):
H_meas[comp_idx, param_idx] = partials[comp_idx]
# Scale Jacobian if you scaled the residuals
if meas.measurementType in [Range.MEASUREMENT_TYPE, RangeRate.MEASUREMENT_TYPE]:
H_meas = H_meas / 1000.0
H_list.append(H_meas)
# Get measurement uncertainty (weight)
sigma = np.array(meas.getTheoreticalStandardDeviation())
if meas.measurementType in [Range.MEASUREMENT_TYPE, RangeRate.MEASUREMENT_TYPE]:
sigma = sigma / 1000.0
# Store weights (one per component)
for comp_idx in range(n_components):
W_list.append(1.0 / sigma[comp_idx]**2)
# Stack all residuals, Jacobians, and weights
Y = np.concatenate(Y_list)
H = np.vstack(H_list)
W = np.array(W_list)
# Solve weighted least squares: H^T W H dx = H^T W Y
sqrtW = np.sqrt(W)
Hw = H * sqrtW[:, None]
Yw = Y * sqrtW
# Solve for parameter corrections in normalized space
try:
dx_norm, _, _, _ = np.linalg.lstsq(Hw, Yw, rcond=None)
except np.linalg.LinAlgError as e:
self.logger.error(f"Linear solve failed: {e}")
return False
# Update normalized parameters
for param_idx, driver in enumerate(all_drivers):
current = driver.getNormalizedValue()
driver.setNormalizedValue(current + dx_norm[param_idx])
# Check convergence (in normalized space)
norm_dx = np.linalg.norm(dx_norm)
self.logger.info(
f"Iteration {iteration+1}: ||dx||={norm_dx:.3e}, "
f"cost={0.5 * np.sum(Yw**2):.6e}"
)
if norm_dx < 1e-3: # Normalized threshold
self.logger.info("Converged!")
# Extract final state from propagator
final_prop = self._builder.buildPropagator(
self._builder.getSelectedNormalizedParameters()
)
final_state = final_prop.getInitialState()
pv = final_state.getPVCoordinates(self._inertial_frame)
pos = pv.getPosition().toArray()
vel = pv.getVelocity().toArray()
self._initial_state[:6] = np.concatenate([pos, vel])
# Update extended state with propagation parameters
if len(self._initial_state) > 6:
state_idx = 6
for driver in prop_drivers.getDrivers():
if driver.isSelected():
self._initial_state[state_idx] = driver.getValue()
state_idx += 1
return True
self.logger.warning(f"Did not converge in {max_iterations} iterations")
return False
This is just a simplfied version of the BLS I am trying to perform, but the important part is the computation of the jacobian. If I can do that, then the rest should be easy.
As a side question, is the way I am using âsetNormalizedValueâ correct? I mean, after the step is selected, if I want to work in scaled/normalized space within the solver (and asuming the system is properly scaled before computing the step) is that how I have to âcorrectâ the state for proper use of the builder? I ask this because I am also worried about numerical instability due to working in physical units (meters and meters/s for the cartesian state, which at some point may lead to corrections below the numerical representation limit if I do it in physical units, ignore the conversion to Km, I may drop that later).
Thanks in advance, I hope this isnt too much of an obvious question about how to use Orekit.