1   /* Copyright 2002-2022 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.utils;
18  
19  import java.util.ArrayList;
20  import java.util.Arrays;
21  import java.util.HashMap;
22  import java.util.List;
23  import java.util.Map;
24  
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.linear.MatrixUtils;
27  import org.hipparchus.linear.RealMatrix;
28  import org.hipparchus.linear.RealVector;
29  import org.orekit.attitudes.Attitude;
30  import org.orekit.attitudes.AttitudeProvider;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.errors.OrekitMessages;
33  import org.orekit.propagation.SpacecraftState;
34  import org.orekit.propagation.numerical.NumericalPropagator;
35  import org.orekit.time.AbsoluteDate;
36  
37  /**
38   * Multiple shooting method using only constraints on state vectors of patch points (and possibly on epoch and integration time).
39   * @see "TRAJECTORY DESIGN AND ORBIT MAINTENANCE STRATEGIES IN MULTI-BODY DYNAMICAL REGIMES by Thomas A. Pavlak, Purdue University"
40   * @author William Desprats
41   * @since 10.2
42   */
43  public abstract class AbstractMultipleShooting implements MultipleShooting {
44  
45      /** Patch points along the trajectory. */
46      private List<SpacecraftState> patchedSpacecraftStates;
47  
48      /** Derivatives linked to the Propagators.
49       * @deprecated as of 11.1 not used anymore
50       */
51      @Deprecated
52      private List<org.orekit.propagation.integration.AdditionalEquations> additionalEquations;
53  
54      /** List of Propagators. */
55      private final List<NumericalPropagator> propagatorList;
56  
57      /** Duration of propagation along each arcs. */
58      private final double[] propagationTime;
59  
60      /** Free components of patch points. */
61      private final boolean[] freePatchPointMap;
62  
63      /** Free epoch of patch points. */
64      private final boolean[] freeEpochMap;
65  
66      /** Number of free variables. */
67      private int nFree;
68  
69      /** Number of free epoch. */
70      private int nEpoch;
71  
72      /** Number of constraints. */
73      private int nConstraints;
74  
75      /** Patch points components which are constrained. */
76      private final Map<Integer, Double> mapConstraints;
77  
78      /** True if orbit is closed. */
79      private boolean isClosedOrbit;
80  
81      /** Tolerance on the constraint vector. */
82      private final double tolerance;
83  
84      /** Maximum number of iterations. */
85      private final int maxIter;
86  
87      /** Expected name of the additional equations. */
88      private final String additionalName;
89  
90      /** Simple Constructor.
91       * <p> Standard constructor for multiple shooting </p>
92       * @param initialGuessList initial patch points to be corrected.
93       * @param propagatorList list of propagators associated to each patch point.
94       * @param additionalEquations list of additional equations linked to propagatorList.
95       * @param arcDuration initial guess of the duration of each arc.
96       * @param tolerance convergence tolerance on the constraint vector.
97       * @param additionalName name of the additional equations
98       * @deprecated as of 11.1, replaced by {@link #AbstractMultipleShooting(List, List, double, double, int, String)}
99       */
100     @Deprecated
101     protected AbstractMultipleShooting(final List<SpacecraftState> initialGuessList, final List<NumericalPropagator> propagatorList,
102                                        final List<org.orekit.propagation.integration.AdditionalEquations> additionalEquations,
103                                        final double arcDuration, final double tolerance, final String additionalName) {
104         this(initialGuessList, propagatorList, arcDuration, tolerance, 1, additionalName);
105         this.additionalEquations = additionalEquations;
106     }
107 
108     /** Simple Constructor.
109      * <p> Standard constructor for multiple shooting </p>
110      * @param initialGuessList initial patch points to be corrected.
111      * @param propagatorList list of propagators associated to each patch point.
112      * @param arcDuration initial guess of the duration of each arc.
113      * @param tolerance convergence tolerance on the constraint vector.
114      * @param maxIter maximum number of iterations
115      * @param additionalName name of the additional equations
116      * @since 11.1
117      */
118     protected AbstractMultipleShooting(final List<SpacecraftState> initialGuessList, final List<NumericalPropagator> propagatorList,
119                                        final double arcDuration, final double tolerance, final int maxIter, final String additionalName) {
120         this.patchedSpacecraftStates = initialGuessList;
121         this.propagatorList = propagatorList;
122         this.additionalName = additionalName;
123         // Should check if propagatorList.size() = initialGuessList.size() - 1
124         final int propagationNumber = initialGuessList.size() - 1;
125         propagationTime = new double[propagationNumber];
126         Arrays.fill(propagationTime, arcDuration);
127 
128         // All the patch points are set initially as free variables
129         this.freePatchPointMap = new boolean[6 * initialGuessList.size()]; // epoch
130         Arrays.fill(freePatchPointMap, true);
131 
132         //Except the first one, the epochs of the patch points are set free.
133         this.freeEpochMap = new boolean[initialGuessList.size()];
134         Arrays.fill(freeEpochMap, true);
135         freeEpochMap[0] = false;
136         this.nEpoch = initialGuessList.size() - 1;
137 
138         this.nConstraints = 6 * propagationNumber;
139         this.nFree = 6 * initialGuessList.size() + 1;
140 
141         this.tolerance = tolerance;
142         this.maxIter   = maxIter;
143 
144         // All the additional constraints must be set afterward
145         this.mapConstraints = new HashMap<>();
146     }
147 
148     /** Get a patch point.
149      * @param i index of the patch point
150      * @return state of the patch point
151      * @since 11.1
152      */
153     protected SpacecraftState getPatchPoint(final int i) {
154         return patchedSpacecraftStates.get(i);
155     }
156 
157     /** Set a component of a patch point to free or not.
158      * @param patchNumber Patch point with constraint
159      * @param componentIndex Component of the patch points which are constrained.
160      * @param isFree constraint value
161      */
162     public void setPatchPointComponentFreedom(final int patchNumber, final int componentIndex, final boolean isFree) {
163         if (freePatchPointMap[6 * (patchNumber - 1) +  componentIndex] != isFree) {
164             nFree += isFree ? 1 : -1;
165             freePatchPointMap[6 * (patchNumber - 1) +  componentIndex] = isFree;
166         }
167     }
168 
169     /** Add a constraint on one component of one patch point.
170      * @param patchNumber Patch point with constraint
171      * @param componentIndex Component of the patch points which are constrained.
172      * @param constraintValue constraint value
173      */
174     public void addConstraint(final int patchNumber, final int componentIndex, final double constraintValue) {
175         final int contraintIndex = (patchNumber - 1) * 6 + componentIndex;
176         if (!mapConstraints.containsKey(contraintIndex)) {
177             nConstraints++;
178         }
179         mapConstraints.put(contraintIndex, constraintValue);
180     }
181 
182     /** Set the epoch a patch point to free or not.
183      * @param patchNumber Patch point
184      * @param isFree constraint value
185      */
186     public void setEpochFreedom(final int patchNumber, final boolean isFree) {
187         if (freeEpochMap[patchNumber - 1] != isFree) {
188             freeEpochMap[patchNumber - 1] = isFree;
189             final int eps = isFree ? 1 : -1;
190             nEpoch = nEpoch + eps;
191         }
192     }
193 
194     /** {@inheritDoc} */
195     @Override
196     public List<SpacecraftState> compute() {
197 
198         if (nFree > nConstraints) {
199             throw new OrekitException(OrekitMessages.MULTIPLE_SHOOTING_UNDERCONSTRAINED, nFree, nConstraints);
200         }
201 
202         int iter = 0; // number of iteration
203 
204         double fxNorm = 0;
205 
206         do {
207 
208             final List<SpacecraftState> propagatedSP = propagatePatchedSpacecraftState(); // multi threading see PropagatorsParallelizer
209             final RealMatrix M = computeJacobianMatrix(propagatedSP);
210             final RealVector fx = MatrixUtils.createRealVector(computeConstraint(propagatedSP));
211 
212             // Solve linear system using minimum norm approach
213             // (i.e. minimize difference between solutions from successive iterations,
214             //  in other word try to stay close to initial guess; this is *not* a least squares solution)
215             // see equation 3.12 in Pavlak's thesis
216             final RealMatrix MMt = M.multiplyTransposed(M);
217             final RealVector dx  = M.transposeMultiply(MatrixUtils.inverse(MMt)).operate(fx);
218 
219             // Apply correction from the free variable vector to all the variables (propagation time, pacthSpaceraftState)
220             updateTrajectory(dx);
221 
222             fxNorm = fx.getNorm() / fx.getDimension();
223 
224             iter++;
225 
226         } while (fxNorm > tolerance && iter < maxIter); // Converge within tolerance and under max iterations
227 
228         return patchedSpacecraftStates;
229     }
230 
231     /** Compute the Jacobian matrix of the problem.
232      *  @param propagatedSP List of propagated SpacecraftStates (patch points)
233      *  @return jacobianMatrix Jacobian matrix
234      */
235     private RealMatrix computeJacobianMatrix(final List<SpacecraftState> propagatedSP) {
236 
237         final int npoints         = patchedSpacecraftStates.size();
238         final int epochConstraint = nEpoch == 0 ? 0 : npoints - 1;
239         final int nrows    = getNumberOfConstraints() + epochConstraint;
240         final int ncolumns = getNumberOfFreeVariables() + nEpoch;
241 
242         final RealMatrix M = MatrixUtils.createRealMatrix(nrows, ncolumns);
243 
244         int index = 0;
245         int indexEpoch = nFree;
246         for (int i = 0; i < npoints - 1; i++) {
247 
248             final SpacecraftState finalState = propagatedSP.get(i);
249             // The Jacobian matrix has the following form:
250 
251             //          [     |     |     ]
252             //          [  A  |  B  |  C  ]
253             // DF(X) =  [     |     |     ]
254             //          [-----------------]
255             //          [  0  |  D  |  E  ]
256             //          [-----------------]
257             //          [  F  |  0  |  0  ]
258 
259             //          [     |     ]
260             //          [  A  |  B  ]
261             // DF(X) =  [     |     ]
262             //          [-----------]
263             //          [  C  |  0  ]
264             //
265             // For a problem with all the components of each patch points is free, A is detailed below :
266             //      [ phi1     -I                                   ]
267             //      [         phi2     -I                           ]
268             // A =  [                 ....    ....                  ]   6(n-1)x6n
269             //      [                         ....     ....         ]
270             //      [                                 phin-1    -I  ]
271 
272             // D is computing according to additional constraints
273             // (for now, closed orbit, or a component of a patch point equals to a specified value)
274             //
275 
276             // If the duration of the propagation of each arc is the same :
277             //      [ xdot1f ]
278             //      [ xdot2f ]                      [  -1 ]
279             // B =  [  ....  ]   6(n-1)x1   and D = [ ... ]
280             //      [  ....  ]                      [  -1 ]
281             //      [xdotn-1f]
282 
283             // Otherwise :
284             //      [ xdot1f                                ]
285             //      [        xdot2f                         ]
286             // B =  [                 ....                  ]   6(n-1)x(n-1) and D = -I
287             //      [                        ....           ]
288             //      [                             xdotn-1f  ]
289             //
290             // If the problem is not dependant of the epoch (e.g. CR3BP), the C and E matrices are not computed.
291             // Otherwise :
292             //      [ -dx1f/dtau1                                           0 ]
293             //      [          -dx2f/dtau2                                  0 ]
294             // C =  [                       ....                            0 ]   6(n-1)xn
295             //      [                               ....                    0 ]
296             //      [                                     -dxn-1f/dtaun-1   0 ]
297             //
298             //      [ -1   1  0                 ]
299             //      [     -1   1  0             ]
300             // E =  [          ..   ..          ] n-1xn
301             //      [               ..   ..     ]
302             //      [                    -1   1 ]
303 
304             final PVCoordinates pvf = finalState.getPVCoordinates();
305 
306             // Get State Transition Matrix phi
307             final double[][] phi = getStateTransitionMatrix(finalState);
308 
309             // Matrix A
310             for (int j = 0; j < 6; j++) { // Loop on 6 component of the patch point
311                 if (freePatchPointMap[6 * i + j]) { // If this component is free
312                     for (int k = 0; k < 6; k++) { // Loop on the 6 component of the patch point constraint
313                         M.setEntry(6 * i + k, index, phi[k][j]);
314                     }
315                     if (i > 0) {
316                         M.setEntry(6 * (i - 1) + j, index, -1);
317                     }
318                     index++;
319                 }
320             }
321 
322             // Matrix B
323             final double[][] pvfArray = new double[][] {
324                 {pvf.getVelocity().getX()},
325                 {pvf.getVelocity().getY()},
326                 {pvf.getVelocity().getZ()},
327                 {pvf.getAcceleration().getX()},
328                 {pvf.getAcceleration().getY()},
329                 {pvf.getAcceleration().getZ()}};
330 
331             M.setSubMatrix(pvfArray, 6 * i, nFree - 1);
332 
333             // Matrix C
334             if (freeEpochMap[i]) { // If this component is free
335                 final double[] derivatives = finalState.getAdditionalState("derivatives");
336                 final double[][] subC = new double[6][1];
337                 for (int j = 0; j < 6; j++) { // Loop on 6 component of the patch point
338                     subC[j][0] = -derivatives[derivatives.length - 6 + j];
339                 }
340                 M.setSubMatrix(subC, 6 * i, indexEpoch);
341                 indexEpoch++;
342             }
343         }
344 
345         for (int j = 0; j < 6; j++) { // Loop on 6 component of the patch point
346             if (freePatchPointMap[6 * (npoints - 1) + j]) { // If this component is free
347                 M.setEntry(6 * (npoints - 2) + j, index, -1);
348                 index++;
349             }
350         }
351 
352 
353         // Matrices D and E.
354         if (nEpoch > 0) {
355             final double[][] subDE = computeEpochJacobianMatrix(propagatedSP);
356             M.setSubMatrix(subDE, 6 * (npoints - 1), nFree - 1);
357         }
358 
359         // Matrices F.
360         final double[][] subF = computeAdditionalJacobianMatrix(propagatedSP);
361         if (subF.length > 0) {
362             M.setSubMatrix(subF, 6 * (npoints - 1) + epochConstraint, 0);
363         }
364 
365         return M;
366     }
367 
368     /** Compute the constraint vector of the problem.
369      *  @param propagatedSP propagated SpacecraftState
370      *  @return constraint vector
371      */
372     private double[] computeConstraint(final List<SpacecraftState> propagatedSP) {
373 
374         // The Constraint vector has the following form :
375 
376         //         [ x1f - x2i ]---
377         //         [ x2f - x3i ]   |
378         // F(X) =  [    ...    ] vectors' equality for a continuous trajectory
379         //         [    ...    ]   |
380         //         [xn-1f - xni]---
381         //         [ d2-(d1+T) ]   continuity between epoch
382         //         [    ...    ]   and integration time
383         //         [dn-(dn-1+T)]---
384         //         [    ...    ] additional
385         //         [    ...    ] constraints
386 
387         final int npoints = patchedSpacecraftStates.size();
388 
389         final double[] additionalConstraints = computeAdditionalConstraints(propagatedSP);
390         final boolean epoch = getNumberOfFreeEpoch() > 0;
391 
392         final int nrows = epoch ? getNumberOfConstraints() + npoints - 1 : getNumberOfConstraints();
393 
394         final double[] fx = new double[nrows];
395         for (int i = 0; i < npoints - 1; i++) {
396             final AbsolutePVCoordinates absPvi = patchedSpacecraftStates.get(i + 1).getAbsPVA();
397             final AbsolutePVCoordinates absPvf = propagatedSP.get(i).getAbsPVA();
398             final double[] ecartPos = absPvf.getPosition().subtract(absPvi.getPosition()).toArray();
399             final double[] ecartVel = absPvf.getVelocity().subtract(absPvi.getVelocity()).toArray();
400             for (int j = 0; j < 3; j++) {
401                 fx[6 * i + j] = ecartPos[j];
402                 fx[6 * i + 3 + j] = ecartVel[j];
403             }
404         }
405 
406         int index = 6 * (npoints - 1);
407 
408         if (epoch) {
409             for (int i = 0; i < npoints - 1; i++) {
410                 final double deltaEpoch = patchedSpacecraftStates.get(i + 1).getDate().durationFrom(patchedSpacecraftStates.get(i).getDate());
411                 fx[index] = deltaEpoch - propagationTime[i];
412                 index++;
413             }
414         }
415 
416         for (int i = 0; i < additionalConstraints.length; i++) {
417             fx[index] = additionalConstraints[i];
418             index++;
419         }
420 
421         return fx;
422     }
423 
424 
425     /** Update the trajectory, and the propagation time.
426      *  @param dx correction on the initial vector
427      */
428     private void updateTrajectory(final RealVector dx) {
429         // X = [x1, ..., xn, T1, ..., Tn, d1, ..., dn]
430         // X = [x1, ..., xn, T, d2, ..., dn]
431 
432         final int n = getNumberOfFreeVariables();
433 
434         final boolean epochFree = getNumberOfFreeEpoch() > 0;
435 
436         // Update propagation time
437         //------------------------------------------------------
438         final double deltaT = dx.getEntry(n - 1);
439         for (int i = 0; i < propagationTime.length; i++) {
440             propagationTime[i] = propagationTime[i] - deltaT;
441         }
442 
443         // Update patch points through SpacecrafStates
444         //--------------------------------------------------------------------------------
445 
446         int index = 0;
447         int indexEpoch = 0;
448 
449         for (int i = 0; i < patchedSpacecraftStates.size(); i++) {
450             // Get delta in position and velocity
451             final double[] deltaPV = new double[6];
452             for (int j = 0; j < 6; j++) { // Loop on 6 component of the patch point
453                 if (freePatchPointMap[6 * i + j]) { // If this component is free (is to be updated)
454                     deltaPV[j] = dx.getEntry(index);
455                     index++;
456                 }
457             }
458             final Vector3D deltaP = new Vector3D(deltaPV[0], deltaPV[1], deltaPV[2]);
459             final Vector3D deltaV = new Vector3D(deltaPV[3], deltaPV[4], deltaPV[5]);
460 
461             // Update the PVCoordinates of the patch point
462             final AbsolutePVCoordinates currentAPV = patchedSpacecraftStates.get(i).getAbsPVA();
463             final Vector3D position = currentAPV.getPosition().subtract(deltaP);
464             final Vector3D velocity = currentAPV.getVelocity().subtract(deltaV);
465             final PVCoordinates pv = new PVCoordinates(position, velocity);
466 
467             //Update epoch in the AbsolutePVCoordinates
468             AbsoluteDate epoch = currentAPV.getDate();
469             if (epochFree) {
470                 if (freeEpochMap[i]) {
471                     final double deltaEpoch = dx.getEntry(n + indexEpoch);
472                     epoch = epoch.shiftedBy(-deltaEpoch);
473                     indexEpoch++;
474                 }
475             } else {
476                 if (i > 0) {
477                     epoch = patchedSpacecraftStates.get(i - 1).getDate().shiftedBy(propagationTime[i - 1]);
478                 }
479             }
480             final AbsolutePVCoordinates updatedAPV = new AbsolutePVCoordinates(currentAPV.getFrame(), epoch, pv);
481 
482             //Update attitude epoch
483             // Last point does not have an associated propagator. The previous one is then selected.
484             final int iAttitude = i < getPropagatorList().size() ? i : getPropagatorList().size() - 1;
485             final AttitudeProvider attitudeProvider = getPropagatorList().get(iAttitude).getAttitudeProvider();
486             final Attitude attitude = attitudeProvider.getAttitude(updatedAPV, epoch, currentAPV.getFrame());
487 
488             //Update the SpacecraftState using previously updated attitude and AbsolutePVCoordinates
489             patchedSpacecraftStates.set(i, new SpacecraftState(updatedAPV, attitude));
490         }
491 
492     }
493 
494     /** Compute the Jacobian matrix of the problem.
495      *  @return propagatedSP propagated SpacecraftStates
496      */
497     private List<SpacecraftState> propagatePatchedSpacecraftState() {
498 
499         final int n = patchedSpacecraftStates.size() - 1;
500 
501         final ArrayList<SpacecraftState> propagatedSP = new ArrayList<>(n);
502 
503         for (int i = 0; i < n; i++) {
504 
505             // SpacecraftState initialization
506             final SpacecraftState augmentedInitialState = getAugmentedInitialState(i);
507 
508             // Propagator initialization
509             propagatorList.get(i).setInitialState(augmentedInitialState);
510 
511             final double integrationTime = propagationTime[i];
512 
513             // Propagate trajectory
514             final SpacecraftState finalState = propagatorList.get(i).propagate(augmentedInitialState.getDate().shiftedBy(integrationTime));
515 
516             propagatedSP.add(finalState);
517         }
518 
519         return propagatedSP;
520     }
521 
522     /** Get the state transition matrix.
523      * @param s current spacecraft state
524      * @return the state transition matrix
525      */
526     private double[][] getStateTransitionMatrix(final SpacecraftState s) {
527         // Additional states
528         final DoubleArrayDictionary dictionary = s.getAdditionalStatesValues();
529         // Initialize state transition matrix
530         final int        dim  = 6;
531         final double[][] phiM = new double[dim][dim];
532 
533         // Loop on entry set
534         for (final DoubleArrayDictionary.Entry entry : dictionary.getData()) {
535             // Extract entry name
536             final String name = entry.getKey();
537             if (additionalName.equals(name)) {
538                 final double[] stm = entry.getValue();
539                 for (int i = 0; i < dim; i++) {
540                     for (int j = 0; j < 6; j++) {
541                         phiM[i][j] = stm[dim * i + j];
542                     }
543                 }
544             }
545         }
546 
547         // Return state transition matrix
548         return phiM;
549     }
550 
551     /** Compute a part of the Jacobian matrix with derivatives from epoch.
552      * The CR3BP is a time invariant problem. The derivatives w.r.t. epoch are zero.
553      *  @param propagatedSP propagatedSP
554      *  @return jacobianMatrix Jacobian sub-matrix
555      */
556     protected double[][] computeEpochJacobianMatrix(final List<SpacecraftState> propagatedSP) {
557 
558         final boolean[] map = getFreeEpochMap();
559 
560         final int nFreeEpoch = getNumberOfFreeEpoch();
561         final int ncolumns = 1 + nFreeEpoch;
562         final int nrows = patchedSpacecraftStates.size() - 1;
563 
564         final double[][] M = new double[nrows][ncolumns];
565 
566         // The Jacobian matrix has the following form:
567 
568         //      [-1 -1   1  0                 ]
569         //      [-1     -1   1  0             ]
570         // F =  [..          ..   ..          ]
571         //      [..               ..   ..   0 ]
572         //      [-1                    -1   1 ]
573 
574         int index = 1;
575         for (int i = 0; i < nrows; i++) {
576             M[i][0] = -1;
577             if (map[i]) {
578                 M[i][index] = -1;
579                 index++;
580             }
581             if (map[i + 1]) {
582                 M[i][index] =  1;
583             }
584         }
585 
586         return M;
587     }
588 
589     /** Update the array of additional constraints.
590      * @param startIndex start index
591      * @param fxAdditional array of additional constraints
592      */
593     protected void updateAdditionalConstraints(final int startIndex, final double[] fxAdditional) {
594         int iter = startIndex;
595         for (final Map.Entry<Integer, Double> entry : getConstraintsMap().entrySet()) {
596             // Extract entry values
597             final int    key   = entry.getKey();
598             final double value = entry.getValue();
599             final int np = key / 6;
600             final int nc = key % 6;
601             final AbsolutePVCoordinates absPv = getPatchedSpacecraftState().get(np).getAbsPVA();
602             if (nc < 3) {
603                 fxAdditional[iter] = absPv.getPosition().toArray()[nc] - value;
604             } else {
605                 fxAdditional[iter] = absPv.getVelocity().toArray()[nc - 3] -  value;
606             }
607             iter++;
608         }
609     }
610 
611     /** Compute the additional constraints.
612      *  @param propagatedSP propagated SpacecraftState
613      *  @return fxAdditionnal additional constraints
614      */
615     protected abstract double[] computeAdditionalConstraints(List<SpacecraftState> propagatedSP);
616 
617     /** Compute a part of the Jacobian matrix from additional constraints.
618      *  @param propagatedSP propagatedSP
619      *  @return jacobianMatrix Jacobian sub-matrix
620      */
621     protected abstract double[][] computeAdditionalJacobianMatrix(List<SpacecraftState> propagatedSP);
622 
623 
624     /** Compute the additional state from the additionalEquations.
625      *  @param initialState SpacecraftState without the additional state
626      *  @param additionalEquations2 Additional Equations.
627      *  @return augmentedSP SpacecraftState with the additional state within.
628      *  @deprecated as of 11.1, replaced by {@link #getAugmentedInitialState(int)}
629      */
630     @Deprecated
631     protected SpacecraftState getAugmentedInitialState(final SpacecraftState initialState,
632                                                        final org.orekit.propagation.integration.AdditionalEquations additionalEquations2) {
633         // should never be called, only implementations by derived classes should be called
634         throw new UnsupportedOperationException();
635     }
636 
637     /** Compute the additional state from the additionalEquations.
638      *  @param i index of the state
639      *  @return augmentedSP SpacecraftState with the additional state within.
640      *  @since 11.1
641      */
642     protected SpacecraftState getAugmentedInitialState(final int i) {
643         // FIXME: this base implementation is only intended for version 11.1 to delegate to a deprecated method
644         // it should be removed in 12.0 when getAugmentedInitialState(SpacecraftState, AdditionalDerivativesProvider) is removed
645         // and the method should remain abstract in this class and be implemented by derived classes only
646         return getAugmentedInitialState(patchedSpacecraftStates.get(i), additionalEquations.get(i));
647     }
648 
649     /** Set the constraint of a closed orbit or not.
650      *  @param isClosed true if orbit should be closed
651      */
652     public void setClosedOrbitConstraint(final boolean isClosed) {
653         if (this.isClosedOrbit != isClosed) {
654             nConstraints = nConstraints + (isClosed ? 6 : -6);
655             this.isClosedOrbit = isClosed;
656         }
657     }
658 
659     /** Get the number of free variables.
660      * @return the number of free variables
661      */
662     protected int getNumberOfFreeVariables() {
663         return nFree;
664     }
665 
666     /** Get the number of free epoch.
667      * @return the number of free epoch
668      */
669     protected int getNumberOfFreeEpoch() {
670         return nEpoch;
671     }
672 
673     /** Get the number of constraints.
674      * @return the number of constraints
675      */
676     protected int getNumberOfConstraints() {
677         return nConstraints;
678     }
679 
680     /** Get the flags representing the free components of patch points.
681      * @return an array of flags representing the free components of patch points
682      */
683     protected boolean[] getFreePatchPointMap() {
684         return freePatchPointMap;
685     }
686 
687     /** Get the flags representing the free epoch of patch points.
688      * @return an array of flags representing the free epoch of patch points
689      */
690     protected boolean[] getFreeEpochMap() {
691         return freeEpochMap;
692     }
693 
694     /** Get the map of patch points components which are constrained.
695      * @return a map of patch points components which are constrained
696      */
697     protected Map<Integer, Double> getConstraintsMap() {
698         return mapConstraints;
699     }
700 
701     /** Get the list of patched spacecraft states.
702      * @return a list of patched spacecraft states
703      */
704     protected List<SpacecraftState> getPatchedSpacecraftState() {
705         return patchedSpacecraftStates;
706     }
707 
708     /** Get the list of propagators.
709      * @return a list of propagators
710      */
711     protected List<NumericalPropagator> getPropagatorList() {
712         return propagatorList;
713     }
714 
715     /** Get he flag representing if the orbit is closed.
716      * @return true if orbit is closed
717      */
718     protected boolean isClosedOrbit() {
719         return isClosedOrbit;
720     }
721 
722 }