Using Maneuver Estimation tutorial, receiving error "cannot normalize a zero norm vector"

Orekit MATLAB usage
1

Hello,
I am recreating the Maneuver Estimation tutorial (ManeuverEstimation (Orekit Tutorials 11.4-SNAPSHOT API)) in Matlab, with all the same values from the tutorial except that I am not including any drag in the force modeling or as an estimated parameter.
I am receiving the following error when I run my code:

org.orekit.errors.OrekitException: cannot normalize a zero norm vector
at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:511)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:425)
at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:140)
at org.orekit.estimation.leastsquares.AbstractBatchLSModel.value(AbstractBatchLSModel.java:319)
at
org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:615)
at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:163)
at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:435)
Caused by: org.hipparchus.exception.MathRuntimeException: cannot normalize a zero norm vector
at org.hipparchus.geometry.euclidean.threed.Vector3D.normalize(Vector3D.java:318)
at org.orekit.forces.maneuvers.propulsion.ThrustPropulsionModel.getAcceleration(ThrustPropulsionModel.java:113)
at org.orekit.forces.maneuvers.Maneuver.acceleration(Maneuver.java:203)
at org.orekit.forces.maneuvers.jacobians.TriggerDate.resetState(TriggerDate.java:259)
at org.orekit.forces.maneuvers.trigger.AbstractManeuverTriggers.applyResetters(AbstractManeuverTriggers.java:185)
at org.orekit.forces.maneuvers.trigger.IntervalEventTrigger$Handler.resetState(IntervalEventTrigger.java:210)
at org.orekit.forces.maneuvers.trigger.IntervalEventTrigger$Handler.resetState(IntervalEventTrigger.java:183)
at org.orekit.propagation.events.AbstractDetector.resetState(AbstractDetector.java:185)
at
org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedEventDetector.resetState(AbstractIntegratedPropagator.java:922)
at org.hipparchus.ode.events.EventState.doEvent(EventState.java:500)
at org.hipparchus.ode.AbstractIntegrator.acceptStep(AbstractIntegrator.java:363)
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:285)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:477)
… 7 more

A sampling of my measurements is below:

2022-09-13T12:00:00.000Z PV ----- 463272.147749054 6992267.905691480 1555273.655311108 -4346.548757579 -1038.371604749 5963.139652441
2022-09-13T12:00:10.000Z PV ----- 419803.216684626 6981917.516779606 1614895.745307146 -4346.619793056 -1039.283266380 5962.907011863
2022-09-13T12:00:20.000Z PV ----- 376331.080492988 6971512.613902257 1674523.135533920 -4346.680607224 -1040.236676222 5962.714843002
2022-09-13T12:00:30.000Z PV ----- 332869.864556574 6961120.872731240 1734157.911362629 -4346.722057009 -1041.153670566 5962.493058494
2022-09-13T12:00:40.000Z PV ----- 289405.405203702 6950693.396136270 1793770.079844198 -4346.759444377 -1042.073515034 5962.243837872
2022-09-13T12:00:50.000Z PV ----- 245944.880270599 6940258.646056903 1853401.581105252 -4346.781807465 -1042.996803519 5962.011614959
2022-09-13T12:01:00.000Z PV ----- 202445.485850803 6929835.266586544 1913010.184411759 -4346.814418186 -1043.891207417 5961.760260108
2022-09-13T12:01:10.000Z PV ----- 159010.628113158 6919385.771570295 1972631.955872234 -4346.854427079 -1044.823114368 5961.507242000
2022-09-13T12:01:20.000Z PV ----- 115503.375871160 6908929.359183694 2032236.820163853 -4346.849584498 -1045.710765215 5961.237855919
2022-09-13T12:01:30.000Z PV ----- 72043.900903971 6898471.926309989 2091858.757810761 -4346.865639508 -1046.614597581 5960.972417016
2022-09-13T12:01:40.000Z PV ----- 28591.614317047 6888002.456938398 2151457.352380632 -4346.900638022 -1047.533601477 5960.696394999
2022-09-13T12:01:50.000Z PV ----- -14872.606814529 6877507.550788619 2211042.809691511 -4346.884852296 -1048.445860409 5960.409473730
2022-09-13T12:02:00.000Z PV ----- -58341.177017025 6867036.139279542 2270685.791957838 -4346.892933885 -1049.321146864 5960.118634893
2022-09-13T12:02:10.000Z PV ----- -101821.285431113 6856512.299247375 2330264.668991085 -4346.862336892 -1050.196873780 5959.832420731
2022-09-13T12:02:20.000Z PV ----- -145288.880207306 6846024.063038516 2389869.498986119 -4346.841640675 -1051.091195312 5959.541109460
2022-09-13T12:02:30.000Z PV ----- -188754.298127625 6835526.967262852 2449435.083000462 -4346.841224131 -1051.963978911 5959.205711617
2022-09-13T12:02:40.000Z PV ----- -232225.992985060 6824976.967238658 2509041.913707445 -4346.807048580 -1052.836264942 5958.912286610
2022-09-13T12:02:50.000Z PV ----- -275713.795363755 6814452.730576607 2568634.881663586 -4346.773421931 -1053.705052351 5958.577883703

And here is some information from my BatchLSEstimator:
measurementsParametersDrivers =
[]

orbitalParametersDrivers =
[Px = 463262.733938, Py = 6992295.407876, Pz = 1555261.269564, Vx = -4346.558024, Vy = -1038.3768, Vz = 5963.136782]

propagatorParametersDrivers =
[MAN_0_START = 19439.211, MAN_0_STOP = 23760.789, MAN_0thrust = 0.6, MAN_1_DURATION = 4320.074, MAN_1_MEDIAN = 24629.369, MAN_1thrust = 0.6, MAN_2_DURATION = 1800.0, MAN_2_MEDIAN = 36033.351, MAN_2thrust = 0.6, reflection coefficient = 2.0]

Any helpful hints about what might be causing the “cannot normalize a zero norm vector” error would be welcome.
Thank you,
Erin

2

Hi @erinf,

Looks like one of your thrust vectors is 0.
Can you show us how you define your maneuvers ?

Maxime

3

I have been looking further into this, and I am not sure that I understand what kinds of values are expected in each maneuver driver field. For instance, I receive errors when I try to set an AbsoluteDate as the value for a “_START” driver, even though I think AbsoluteDate is the correct type of value to use.
My code is below.
Thank you for your help,
Erin

I’m defining my maneuvers as follows (from a yaml file that I read as “inputData”, in Matlab):

%Maneuvers
if isfield(inputData.propagator.forceModels,‘maneuvers’) && ~isempty(inputData.propagator.forceModels.maneuvers)
maneuvers = inputData.propagator.forceModels.maneuvers;
for ii = 1:1:size(maneuvers,2)
maneuver = maneuvers{1,ii};
name = maneuver.name;
thrust = maneuver.thrust.values{1};
isp = maneuver.isp;
direction = org.hipparchus.geometry.euclidean.threed.Vector3D(maneuver.direction{1},…
maneuver.direction{2},maneuver.direction{3})

    %Compute values of firingDate and duration from parameters to estimate
    duration = 0.0;
    if isfield(maneuver,'startDate') && isfield(maneuver,'stopDate')
        firingDate = org.orekit.time.AbsoluteDate(java.lang.String(maneuver.startDate.values{1}),org.orekit.time.TimeScalesFactory.getUTC())
        thisStopDate = org.orekit.time.AbsoluteDate(java.lang.String(maneuver.stopDate.values{1}),org.orekit.time.TimeScalesFactory.getUTC());
        duration = thisStopDate.durationFrom(firingDate)
    elseif isfield(maneuver, 'medianDate') && isfield(maneuver,'duration')
        duration = maneuver.duration.values{1}
        thisMedianDateValues1 = org.orekit.time.AbsoluteDate(java.lang.String(maneuver.medianDate.values{1}),org.orekit.time.TimeScalesFactory.getUTC());
        firingDate = thisMedianDateValues1.shiftedBy(-duration/2)
        %         else
        %             throw org.orekit.errors.OrekitException(OrekitMessages.INCONSISTENT_SELECTION, maneuver.startDate, maneuver.stopDate,...
        %                 maneuver.medianDate, maneuver.duration);
    end
            
    constantThrustPropulsionModel = org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel(thrust, isp, direction, name);
    constThrustMnv = org.orekit.forces.maneuvers.ConstantThrustManeuver(builder.getAttitudeProvider(),...
            org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers(name, firingDate, duration),...
            constantThrustPropulsionModel);
    builder.addForceModel(constThrustMnv);
    org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers(name, firingDate, duration).getStartDate()
    drivers = constThrustMnv.getParametersDrivers()
    
    if isfield(maneuver,'startDate')
        for jj = 1:1:drivers.size()
            if drivers.get(jj-1).getName().equals(string([name, '_START']))
                drivers.get(jj-1).setValue(firingDate); %org.orekit.time.AbsoluteDate(java.lang.String(maneuver.startDate.values{1}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                if maneuver.startDate.isEstimated
                    drivers.get(jj-1).setSelected(true);
                    drivers.get(jj-1).setMaxValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.startDate.values{3}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                    drivers.get(jj-1).setMinValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.startDate.values{2}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                end
            end
        end
    end

    if isfield(maneuver,'stopDate')
        for jj = 1:1:drivers.size()
            if drivers.get(jj-1).getName().equals(string([name, '_STOP']))
                drivers.get(jj-1).setValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.stopDate.values{1}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                if maneuver.stopDate.isEstimated
                    drivers.get(jj-1).setSelected(true);
                    drivers.get(jj-1).setMaxValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.stopDate.values{3}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                    drivers.get(jj-1).setMinValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.stopDate.values{2}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                end
            end
        end
    end

    if isfield(maneuver,'medianDate')
        for jj = 1:1:drivers.size()
            if drivers.get(jj-1).getName().equals(string([name, '_MEDIAN']))
                drivers.get(jj-1).setValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.medianDate.values{1}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                if maneuver.medianDate.isEstimated
                    drivers.get(jj-1).setSelected(true);
                    drivers.get(jj-1).setMaxValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.medianDate.values{3}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                    drivers.get(jj-1).setMinValue(org.orekit.time.AbsoluteDate(java.lang.String(maneuver.medianDate.values{2}),org.orekit.time.TimeScalesFactory.getUTC()).durationFrom(builder.getInitialOrbitDate()));
                end
            end
        end
    end

    if isfield(maneuver,'duration')
        for jj = 1:1:drivers.size()
            if drivers.get(jj-1).getName().equals(string([name, '_DURATION']))
                drivers.get(jj-1).setValue(duration);
                if maneuver.duration.isEstimated
                    drivers.get(jj-1).setSelected(true);
                    drivers.get(jj-1).setMinValue(maneuver.duration.values{2});
                    drivers.get(jj-1).setMaxValue(maneuver.duration.values{3});
                end
            end
        end
    end
    
    if isfield(maneuver,'thrust')
        for jj = 1:1:drivers.size()
            if drivers.get(jj-1).getName().equals(string([name, 'thrust']))
                drivers.get(jj-1).setValue(thrust);
                if maneuver.thrust.isEstimated
                    drivers.get(jj-1).setSelected(true);
                    drivers.get(jj-1).setMinValue(maneuver.thrust.values{2});
                    drivers.get(jj-1).setMaxValue(maneuver.thrust.values{3});
                end
            end
        end
    end
     
    if isfield(maneuver,'flowRate')
        for jj = 1:1:drivers.size()
            if drivers.get(jj-1).getName().equals(string([name, 'flow rate']))
                drivers.get(jj-1).setValue(maneuver.flowRate.values{1});
                if maneuver.flowRate.isEstimated
                    drivers.get(jj-1).setSelected(true);
                    drivers.get(jj-1).setMinValue(maneuver.flowRate.values{2});
                    drivers.get(jj-1).setMaxValue(maneuver.flowRate.values{3});
                end
            end
        end
    end
end

end

4

Here is the yaml information for the measurement generation (based on maneuver-estimation-generate-measurements.yaml from the orekit tutorial):
- name: “MAN_1”
startDate:
values: [“2022-09-13T18:16:29.000”, “2000-01-01T00:00:00.000”, “2100-01-01T00:00:00.000”]
isEstimated: false
stopDate:
values: [“2022-09-13T19:24:20.000”, “2000-01-01T00:00:00.000”, “2100-01-01T00:00:00.000”]
isEstimated: false
thrust:
values: [0.7, 0.0, 1.0]
isEstimated: false
isp: 2900
direction: [1.0, 0.0, 0.0]
- name: “MAN_2”
medianDate:
values: [“2022-09-13T22:03:26.000”, “2000-01-01T00:00:00.000”, “2100-01-01T00:00:00.000”]
isEstimated: false
duration:
values: [1821.0, 0.0, 4000.0]
isEstimated: false
thrust:
values: [0.5, 0.0, 1.0]
isEstimated: false
isp: 2900
direction: [0.0, 0.0, 1.0]

Here is the yaml information for the estimation (based on maneuver-estimation.yaml from the orekit tutorial):
- name: “MAN_1”
medianDate:
values: [“2022-09-13T18:50:29.369”, “2000-01-01T00:00:00.000”, “2100-01-01T00:00:00.000”]
isEstimated: true
duration:
values: [4320.074, 4000, 5000]
isEstimated: true
thrust:
values: [0.6, 0.0, 1.0]
isEstimated: true
isp: 2900
direction: [1.0, 0.0, 0.0]
- name: “MAN_2”
medianDate:
values: [“2022-09-13T22:00:33.351”, “2000-01-01T00:00:00.000”, “2100-01-01T00:00:00.000”]
isEstimated: true
duration:
values: [1800.0, 0.0, 4000.0]
isEstimated: true
thrust:
values: [0.6, 0.0, 1.0]
isEstimated: true
isp: 2900
direction: [0.0, 0.0, 1.0]

5

No, it’s a DateDriver object so it’s actually an offset (in seconds) from a base date (the one you gave it at creation).

Thats a lot of code :sweat_smile:
Can you try constantThrustPropulsionModel.getDirection() after the creation of the object constantThrustPropulsionModel ?

6

Okay, thank you for the clarification on the DateDriver.
Just to be sure, DateBasedManeuverTriggers takes AbsoluteDate for the start date/firing date, right?

To answer your question:
constantThrustPropulsionModel = org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel(thrust, isp, direction, name);
constantThrustPropulsionModel.getDirection()

The code above produces {1; 0; 0} and then {0; 0; 1} (since I’m iterating over 2 maneuvers).

Then I have this block of code:
constThrustMnv = org.orekit.forces.maneuvers.ConstantThrustManeuver(builder.getAttitudeProvider(),…
org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers(name, firingDate, duration),…
constantThrustPropulsionModel);
builder.addForceModel(constThrustMnv);
org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers(name, firingDate, duration).getStartDate()
drivers = constThrustMnv.getParametersDrivers()

And that code produces the following:
2022-09-13T18:16:29.000Z

drivers = [MAN_1thrust = 0.7, MAN_1flow rate = -2.4613839623605163E-5, MAN_1_START = 0.0, MAN_1_STOP = 0.0, MAN_1_MEDIAN = 0.0, MAN_1_DURATION = 4071.0]

I wonder why “drivers” has START, STOP, and MEDIAN values all 0.0? Could that be the issue?

7

Yes

Because START, STOP and MEDIAN are deviations in seconds from the initial dates in input.
So, at initialisation, the deviations are 0s. And then they are updated when performing the estimation.

8

Hello @MaximeJ,
Thank you very much for your help so far. I am making good progress. I am no longer receiving the “cannot normalize a zero norm vector” error. But when I try to estimate my maneuvers, all of the final estimated values are reasonable except the dates (start, stop, and median dates).

I have one remaining question which I think will solve my issue if you can help me understand.
First I build a constantThrustPropulsion Model like this:
constantThrustPropulsionModel = org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel(thrust, isp, direction, name);

Then I build a constantThrustManeuver like this:
constThrustMnv = org.orekit.forces.maneuvers.ConstantThrustManeuver(builder.getAttitudeProvider(),…
org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers(name, firingDate, duration),…
constantThrustPropulsionModel);

Then I add the constantThrustManeuver to my propagator builder like this:
builder.addForceModel(constThrustMnv);

Then I find that if I check the drivers as:
drivers = constThrustMnv.getParametersDrivers()
The result is, for instance:
drivers = [MAN_1thrust = 0.6, MAN_1flow rate = -2.1097576820232998E-5, MAN_1_START = 0.0, MAN_1_STOP = 0.0, MAN_1_MEDIAN = 0.0, MAN_1_DURATION = 4320.074]

This result for drivers is confusing to me based on what you said above because every driver value which is not a date is the actual initial guess that I provided. (For instance, MAN_1thrust = 0.6 matches the input I provided for thrust, and MAN_1_DURATION = 4320.074 matches the input I provided for duration.) But the drivers which are dates (MAN_1_START = 0.0, MAN_1_STOP = 0.0, MAN_1_MEDIAN = 0.0) are all zero.

Additionally, the “-log.out” file shows the following:
Estimated propagation parameters changes:
1 MAN_0_START -253.9566 (final value: -253.9566)
2 MAN_0_STOP 157.2162 (final value: 157.2162)
3 MAN_0thrust -0.099872 (final value: 0.50013)
4 MAN_1_DURATION -249.2154 (final value: 4070.8586)
5 MAN_1_MEDIAN -4.227 (final value: -4.227)
6 MAN_1thrust 0.099652 (final value: 0.69965)
7 MAN_2_DURATION 26.3 (final value: 1826.3)
8 MAN_2_MEDIAN 172.6666 (final value: 172.6666)
9 MAN_2thrust -0.10097 (final value: 0.49903)

Here, MAN_1thrust final value 0.69965 is very close to the actual value that I used to generate the measurements (actual value 0.7), and the same is true for the other thrust values and the duration values. But the START, STOP, and MEDIAN final values still just seem to be showing some change from 0 rather than a change from an initial guess.

Am I missing somewhere that I need to provide an initial guess for the START, STOP, and MEDIAN dates? I do not think I am successfully setting these initial guesses for dates when I run the ConstantThrustManeuver constructor as I showed above:
constThrustMnv = org.orekit.forces.maneuvers.ConstantThrustManeuver(builder.getAttitudeProvider(),…
org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers(name, firingDate, duration),…
constantThrustPropulsionModel);

As always, thank you for any help!
Erin

9

Update:
I have now figured out that, in the case of dates, if I add the final estimated values (e.g. MAN_0_START -253.9566 (final value: -253.9566)) as seconds to the initial guesses that I input, then I get the final date estimate that I am looking for.
(Yay!)
Final question: In the case of dates parameters (start, stop, and median), is there any way to check what initial guess has been loaded into the parameters drivers, after it has been loaded?

In other words, the following code doesn’t tell me what initial guess has been loaded for START, STOP, or MEDIAN maneuver dates:
drivers = constThrustMnv.getParametersDrivers()
Because the result is:
drivers = [MAN_1thrust = 0.6, MAN_1flow rate = -2.1097576820232998E-5, MAN_1_START = 0.0, MAN_1_STOP = 0.0, MAN_1_MEDIAN = 0.0, MAN_1_DURATION = 4320.074]

Is there a way to read out the initial guesses for dates?
Thank you again,
Erin

10

Hi @erinf ,

Yes you simply have to cast the driver to a DateDriver and use method getBaseDate(). But I don’t know how to do the cast in Matlab.
In Java, supposing that you have isolated the driver for the MAN_1_START parameter.
Cast it: DateDriver dateDriver = (DateDriver) man1StartDriver;
Get the base date: AbsoluteDate man1StartBaseDate = dateDriver.getBaseDate();

Hope this helps,
Maxime