[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Orekit Developers] recursive orbital estimation

brianh4321@hotmail.com a écrit :


Hello/Bonjour


Hi Brian and everyone,



I've been testing a recursive orbital estimator using by importing the orekit
libraries directly into matlab. Being somewhat new to the field I've had many
questions and Luc and his team have been very gracious with their help. To
further these discussions in a larger domain I wish to post the following
questions.

Not sure if anyone works much with matlab, but I've been testing an recursive
orbital estimator using the new NumericalPropagatorWithJacobians (which imports 
directly into matlab).

1. The orbital state in numerical propagator is equinoctial expressed in
"true latitude", however most books seem to express equinoctial in terms of
longitude. I seem to find that while the equations match most reference as
longitude, the description of your functions seems to replace latitude for
longitude. Is there something I'm missing?
I'd rather think I missed something myself. We have already seen some other cases where the names historically used and thought to be general ones were in fact local names that did not match global consensus.
I'll check for this one but am pretty sure I used a local name here also. If so, I think we should change the name. The proper way to do this is to first add getters with the new names while simply deprecating the other ones, and when the next major release is published, we remove the deprecated names. This simplifies smooth transitions for existing code. 



2. Also, most reference expresses equinoctial in terms of mean
latitude/longitude as oppose to true latitude/longitude. I was hoping to
find a reference that expresses the measurement Jacobian as a function of
your state. Do anyone know of such a reference? (A Jacobian relating your state 
and a cartesian state would also work). I'm currently computing the
measurement Jacobian with finite differencing, but I'd love to have a closed
form.
Concerning only the mean/true values, you can convert derivatives using total differentials. Here are some equations (using LaTeX notation) you can use: 

dv = \frac{\xi^2}{\varepsilon^3} dM + \frac{(1+\xi)\nu}{e\varepsilon^2} de

where

\xi = 1+e\cos v
\varepsilon = \sqrt{1-e^2}
Since \lambda_v = \Omega + \omega + v and \lambda_M = \Omega + \omega + M one can easily use the first equation to compute d\lambda_v from d\lambda_M.
We are going to add state derivatives between all parameters types and Cartesian parameters in the code base very soon, so converting Jacobians should be a matter of performing simple 6x6 linear algebra.
The current code in the Git repository has been vastly improved with respect to jacobians recently. Mainly the experimental NumericalPropagatorWithJacobians class introduced with 5.0 has been completely merged within the regular NumericalPropagator class using a new "additional equations" feature that allow people to add their own equations to be integrated together with the trajectory. For now only the Cartesian Jacobians are in the repository, the equinoctial Jacobians will follow and the remaining ones afterwards.
Véronique will present this in the upcoming ISSFD conference in Brazil at the end of the month. We will probably put the paper online after the conference, and we will certainly explain in more details how to use it in some tutorials in the wiki.
3. Here is a sub-section of the Kalman estimator code. It works but the answers 
don't seem to converge with the accuracy I'd expect. My initial state
covariance is somewhat of a guess, I don't have much experience with
equinoctial and they are somewhat of a guess. The problem is derived from an
example in Satellite Orbits, pg 291 (by Oliver MontenBruck)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pest = diag([10e5,1e-2,1e-2,1e-4,1e-3,1e-3].^2);
Xest= Xref0;
Xref= Xref0;
b=zeros(3,N_obs);
sigma_angle= .01*pi/180;
sigma_range= 10;
Q= zeros(6);
linear= 1;
for ii=1:N_obs
    date= propagator.getInitialState().getDate();
    dt= Obs{ii,1}.durationFrom(date);

    % Propagate
    predictedState = propagator.propagate(date.shiftedBy(dt), PhiEq, dFdP);
The line above will need some change with the current code. As we merged NumericalPropagatorWithJacobians and NumericalPropagator, we preserved the general interface for the propagate method so it does not accept the PhiEq and dFdP extra parameters anymore. The new way to handle this is to register an instance of the PartialDerivativesEquations class to the numerical propagator as follows: 

  NumericalPropagator propagator = new NumericalPropagator(integrator);
  PartialDerivativesEquations partialDerivatives =
    new PartialDerivativesEquations(propagator);
As of today, you should also call propagator.addAdditionalEquations(additionalEquations); but this need will be removed in the next few days, as the method will be called automatically at the end of PartialDerivativesEquations constructor. I cannot commit the change yet as I am struggling with another problem at the moment.
After that, you will call propagate with only the target date as argument and will retrieve the Jacobians using 

partialDerivatives.additionalEquations(PhiEq, dFdP);



    predictedPos= predictedState.getPVCoordinates().getPosition();
    % Predicted Measurements
    Elev= staF.getElevation(predictedPos,inertialFrame,Obs{ii,1});
    Azim= staF.getAzimuth(predictedPos,inertialFrame,Obs{ii,1});
    Range= staF.getRange(predictedPos,inertialFrame,Obs{ii,1});
    % Store predicted reference and compute Phi, Inverse Measurement Cov(R)
    % and Measurement Jacobian H (finite differencing)
    oXref= Xref;
    Xref= [ predictedState.getA();
	predictedState.getEquinoctialEx();
	predictedState.getEquinoctialEy();
	predictedState.getHx();
	predictedState.getHy();
	predictedState.getLv()];
    mPhiEq= ja2mat(PhiEq);
    mPhiEq= mPhiEq(1:6,1:6);
    Rinv= diag([sigma_angle*cos(Elev),sigma_angle,sigma_range].^2);
    G= computeH(staF,predictedState);

    % Time Update
    Pbar= mPhiEq*Pest*mPhiEq.'+Q;
    if(linear)
	% Linearized Kalman Update
	Xbar= Xref+ mPhiEq*(Xest-oXref);
	% Measurement Update
	b(:,ii)= [Obs{ii,2};Obs{ii,3};Obs{ii,4}]-[Azim;Elev;Range];
	K= Pbar*G'*inv(G*Pbar*G'+Rinv);
	Xest= Xbar+K*(b(:,ii)-G*(Xbar-Xref));
	Pest= Pbar-K*G*Pbar;
	propagator.setInitialState(predictedState);
    else
	% Extended Kalman Update
	Xbar= Xref;
	b(:,ii)= [Obs{ii,2};Obs{ii,3};Obs{ii,4}]-[Azim;Elev;Range];
	K= Pbar*G'*inv(G*Pbar*G'+Rinv);
	Pest= Pbar-K*G*Pbar;
	Xest= Xbar+K*b(:,ii);
	stateEst= EquinoctialOrbit(Xest(1),Xest(2), Xest(3), Xest(4), ...
	    Xest(5), Xest(6),EquinoctialOrbit.TRUE_LATITUDE_ARGUMENT,...
	    FramesFactory.getEME2000(), Obs{ii,1}, mu);
	propagator.setInitialState(SpacecraftState(stateEst));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I am not sure, but it may be related to the intialization of the Jacobians when you reset the initial state. In the new implementation, there is a dedicated setInitialJacobians to reset the start state of the jacobians before calling propagate. Depending on what you do with your Jacobians, you may either need to reset them to Identity and Zero each time you reset the state (so they are computed with respect to the new initial state at each sub-propagation) or you may force them to the value they have at the end of the previous propagation so they are computed with respect to the original initial state). 

Luc



Any help or comment are appreciated.

Brian






----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.