A Novel Approach for Optimal Trajectory Design with Multiple Operation Modes of Propulsion System, Part 1

• The specific problem is generating low-thrust fuel-optimal trajectories for spacecraft whose propulsion system has multiple modes of operation.Source: §1 · p. 3
• The problem is especially challenging because low-thrust heliocentric trajectories have weak local-control sensitivity at the final state, creating potential uniqueness and convergence issues.These difficulties are compounded for many-revolution problems, particularly when single shooting is used within homotopy procedures.Source: §1 · pp. 3–5
• Homotopy methods address difficult TPBVPs by solving an easier neighboring problem first and recursively continuing converged solutions toward the original problem, although global convergence is rarely guaranteed.Source: §1 · pp. 3–4
• Indirect methods can provide high-resolution solutions and superior speed performance when they converge, but finding the missing costate boundary conditions remains numerically difficult.Necessary conditions can reveal control structure, yet they do not directly provide the costate conditions that drive the system to the desired extremal.Source: §1 · p. 5
• Composite Smooth Control (CSC) is introduced as a new framework within indirect optimization, motivated partly by a prior framework for continuous multi-phase systems.The paper series targets challenges in characterizing optimal trajectories for systems with multiple propulsion-operation modes.Source: §1 · pp. 3, 5

Scope

• The proposed framework identifies optimal transitions between operation modes for a VIVT engine and the optimal number of engines for an engine cluster in the two-paper set.The paper’s first part addresses variable-Isp, variable-thrust engines, while part 2 addresses engine clusters.Source: §1 · pp. 5–6

Framework

• CSC uses Hyperbolic Tangent Smoothing to connect propulsion controls into a continuous composite control spanning the entire admissible spectrum.The smoothing parameter controls switch sharpness homotopically while retaining transitions between distinct operation modes.Source: §1 · pp. 5–6

Bang-bang controls

• For bang-bang controls, CSC applies HTS to switching functions, making switch times implicit and autonomously determined while incorporating switching-function-triggered constraints.The switching function S is imposed as the constraint S = 0 and serves as the HTS argument.Source: §1 · p. 6

Numerical benefits

• The mathematically continuous composite input significantly enhances convergence for adaptive differential solvers and Newton-type schemes used for TPBVPs.The framework supports unified single-shooting solutions for low-thrust interplanetary heliocentric trajectory optimization without separate algorithms or iterative switch-time determination during integration.Source: §1 · pp. 6–7

Costate modeling

• A second contribution models planetary gravitational perturbations in costate equations through complex-based numerical Jacobian derivatives when modified equinoctial elements represent the trajectory.The Jacobian of the Hamiltonian supplies the costate dynamics.Source: §1 · p. 7

2.1. Power Sub-system Modeling

• The modeled solar contribution changes with spacecraft distance from the Sun and solar-cell efficiency degradation over time, approximated by distance- and time-dependent multiplicative functions.The aging term uses a decay rate σ measured as a percentage.Source: §2.1 · p. 8
• For trajectories entering inside 1 AU, simplified inverse-square power models can produce errors larger than 25%, whereas differences between approximations are negligible beyond 1 AU.More accurate power models are therefore required for inner-planet trajectories and multi-revolution transfers passing inside 1 AU.Source: §2.1 · p. 10

2.2. Variable-Isp, Variable-Thrust Engine Model

• The propulsion model uses input power P and exhaust velocity c = Ispg0 as bounded real-valued controls, allowing simultaneous variation of specific impulse and thrust.Thrust depends on power input, exhaust velocity, and constant engine efficiency, while the resulting thrust acceleration is parameterized by magnitude and a steering unit vector.Source: §2.2 · pp. 10–11

2.3. Secondary-Body Perturbation Modeling

• The heliocentric flight model uses two-body dynamics for preliminary analysis while deriving secondary-body accelerations and corresponding costate dynamics for higher-fidelity perturbation studies.The study specifically investigates the importance of secondary-body perturbations in its test cases, noting prior evidence that power-system modeling can affect solutions more strongly than omitted higher-fidelity gravitational perturbations.Source: §2.3 · pp. 11–12
• Secondary-body gravitational accelerations are represented using gravitational parameters and transformed from the inertial frame into LVLH components defined by osculating unit vectors.The LVLH basis uses radial, transverse, and angular-momentum directions, with a direction-cosine matrix formed from these orthogonal unit vectors.Source: §2.3 · p. 12

2.4. Equations of Motion

• The equations of motion use modified equinoctial elements and spacecraft mass, separating propulsion acceleration from secondary-body perturbations.The formulation separates MEE dynamics from mass-rate dynamics to facilitate derivation of optimality conditions.Source: §2.4 · p. 13
• The complete acceleration includes engine-induced and secondary-body terms, whereas smaller disturbances such as solar radiation pressure are not considered.Secondary-body perturbations represent the gravitational effects of the other planets in the Solar System.Source: §2.4 · p. 13

3. Fuel-Optimal Indirect Formulation for VIVT Engines

• The fuel-optimal indirect formulation derives the Hamiltonian, costate dynamics, boundary conditions, and Jacobian requirements underlying the associated TPBVP.A complex-based derivative method is also presented as an alternative for deriving costate dynamics.Source: §3 · p. 14

3.1. Formulation of the Fuel-Optimal Problem

• The Hamiltonian minimizes thrust direction and power through PMP, with power and thrust selected at their maxima when the switching function S is positive and zero when S is negative.The formulation assumes singular arcs with S = 0 do not occur over finite intervals.Source: §3.1 · pp. 14–15
• The optimal specific impulse follows ∂H/∂c = 0 when admissible, while the nearest bound is selected when the unconstrained value lies outside [cmin, cmax].This max-min implementation is inherently nonsmooth because it requires conditional selection among exhaust-velocity cases.Source: §3.1 · pp. 14, 16
• The framework seeks smooth transitions across the three specific-impulse cases and possible power inputs to avoid nonsmooth on-off controls and enable smoother multi-mode thrust embedding.The stated motivation is to make numerical solutions more convenient, efficient, and accurate.Source: §3.1 · pp. 16–17
• Secondary-body perturbations do not change the optimal-control expressions because they depend only on position and velocity, but they alter costate histories and thereby the resulting control profiles.Costate dynamics are first derived with solar-body perturbations set to zero, with alternatives reserved for perturbed cases.Source: §3.1 · p. 17

3.2. Two-Point Boundary-Value Problem

• For time-fixed rendezvous problems, the spacecraft’s final mass is free, its position and velocity match the target body, seven final equality constraints apply, and the final mass costate is -1.Source: §3.2 · pp. 17–18
• For multi-revolution transfers, the unknown revolution count Nrev is enforced through the true-longitude boundary condition alongside the target state constraints.The target longitude is defined relative to the initial target-body longitude for Nrev = 0.Source: §3.2 · p. 18
• Considering extremals for every feasible revolution count can identify the global extremal for time-fixed rendezvous-type fuel-optimal trajectories.The formulation assumes one local extremal for each feasible integer number of revolutions.Source: §3.2 · p. 18
• The indirect formulation solves a TPBVP by determining seven initial costate components to satisfy seven prescribed boundary constraints.Initial state and mass are specified, while final state and costates depend on the unknown initial costates.Source: §3.2 · p. 19

3.3. Derivation of Costate Dynamics

• In indirect optimization methods, costate derivation and coding can become tedious and error-prone for complicated systems; coordinate choice significantly affects both derivation and solution in low-thrust trajectory optimization.Source: §3.3 · p. 19

3.3.1. Symbolic Approach

• The symbolic method’s derivational complexity does not generally predict numerical convergence efficiency or reliability.For restricted dynamics, secondary bodies do not contribute to the mass-associated costate because spacecraft mass is neglected relative to planetary masses.Source: §3.3 · p. 19, §3.3.1 · p. 21
• Symbolic toolboxes avoid error-prone costate-equation derivation, construct necessary expressions, and automate code generation, while perturbing accelerations affect computational efficiency.Source: §3.3.1 · p. 20

3.3.2. Numerical Approach Using Complex-Based Derivative Approach

• The CX method generates costates and state-transition derivatives numerically, requiring more computation typically but reducing human implementation effort and code size.It can treat interpolated data tables as black boxes, unlike the symbolic approach, although convergence near discontinuous interpolation derivatives can produce anomalies.Source: §3.3.2 · pp. 21–22
• In some problems, CX evaluation becomes faster than lengthy symbolic costate expressions, especially in MATLAB; interpolation discontinuities can still cause occasional anomalies.Local Lagrange interpolation around nearby table entries is proposed to avoid derivative discontinuities.Source: §3.3.2 · p. 22
• Near-machine-precision derivative accuracy is obtained across a wide range of complex step sizes, unlike real-valued finite differences requiring a “Goldilocks” region.The accuracy is reported relative to analytical partial derivatives.Source: §3.3.2 · pp. 22–23
• The CX method numerically evaluates the Hamiltonian rather than constructing a symbolic Jacobian, enabling treatment of complicated dynamics.The Hamiltonian is evaluated using complex perturbations of individual state components.Source: §3.3.2 · p. 23
• Costate reconstruction follows the Euler–Lagrange equations while holding optimal controls fixed during complex perturbations, preventing control-induced changes in the gradient.Controls are computed from the original real state vector, while state-dependent noncontrol terms receive the imaginary perturbation.Source: §3.3.2 · p. 23
• For Cartesian fuel-optimal dynamics, the thrust-direction unit vector is smooth and depends only on the velocity costate, whereas throttle depends on mass through the switching function.This distinction creates the need to preserve the control while evaluating costate derivatives.Source: §3.3.2 · pp. 23–24
• MEE dynamics are more challenging because both optimal direction and switching function depend on the state, requiring careful treatment of state-dependent control terms.The method keeps controls unchanged where they enter the Hamiltonian and perturbs the remaining state-dependent terms.Source: §3.3.2 · p. 24
• The approach is independently parallelizable across the seven state perturbations, although the authors did not use their parallel implementation because the example was computationally undemanding.Each state can be perturbed separately in the complex component.Source: §3.3.2 · pp. 24–25
• Secondary-body gravitational perturbations converge in a small, single-digit number of iterations and take only a few seconds for the considered test case.Convergence is facilitated by starting from an excellent near-optimal solution.Source: §3.3.2 · p. 25
• Applying the CX-based method with secondary-body gravity and MEE dynamics is identified as another contribution of the work.The passage specifically frames these extensions as contributions of the method.Source: §3.3.2 · p. 25

4. Control Smoothing for VIVT Engines

• HTS smoothing enlarges the TPBVP convergence domain, making convergence usually achievable from a moderate number of random initial costate guesses.The method replaces instantaneous switches with controllably sharp, rigorously continuous approximations.Source: §4 · pp. 25–26
• Continuous switch approximations eliminate high-precision switch-time root solving when modern adaptive step-size integrators are used.Without smoothing, conventional integrators must locate state- and costate-dependent force discontinuities and adjust integration steps accordingly.Source: §4 · p. 26
• Composite controls replace discontinuous “if →then” conditions by weighting distinct controls with constraint-based activation functions and summing them into a continuous control.Activation functions approach unity when constraints are satisfied and zero when violated; multiplicative weighting allows one control to dominate smoothly.Source: §4 · p. 26
• For VIVT engines, separate smoothing factors approximate the specific-impulse control and bang-off-bang power control, while singular-type power arcs are excluded because they rarely occur in optimal space flights.The smoothing parameters ρb and ρc can differ by problem and are reduced to approximate Pontryagin-optimal instantaneous switches.Source: §4 · p. 29

4.1. Generalization of Composite Smooth Controls

• The composite smooth-control formulation is generalized by defining each of m scalar control inputs from nb,k smoothly merged building-block controls.The subsection introduces the notation for the i-th building-block control and the composite representation of the k-th control time history.Source: §4.1 · pp. 29–30
• The generalized Composite Smooth Control formulation replaces activation functions with smooth approximations for multiple building-block controls and their associated constraints.For the example with c* ∈ {cmin, cop, cmax}, the three controls have one, two, and one constraints, respectively.Source: §4.1 · p. 30

4.2. Consideration of No-Power Condition

• A user-defined minimum available-power limit Pava,min enables smooth handling of no-power conditions that can occur farther from the Sun in solar electric propulsion.The same smoothing parameter ρb used for bang-bang power behavior is applied to avoid introducing another parameter, and the activation function multiplies the power approximation.Source: §4.2 · pp. 30–31

4.3. Imposing No-Thrust Time Interval and Thrust Vector Direction Constraints

• Smooth activation functions can impose no-thrust intervals, state-dependent thrust patterns, shadow constraints, and specified time-varying or fixed thrust-vector directions.No-thrust intervals support missed-thrust analyses and science-measurement windows; during a constrained interval, the thrust direction is reoriented toward the desired direction from the primer-vector direction.Source: §4.3 · pp. 31–32

4.4. Connection Between Smooth and Non-smooth Necessary Conditions

• Smooth homotopy solutions can approach discontinuous multi-mode switches, recover the original MPBVP optimality conditions, and remain computationally more attractive than formal discontinuous solutions.The smoothing parameters ρb and ρc approach limiting values, while nonzero values can often achieve engineering optimality with negligible sacrifice on optimality.Source: §4.4 · pp. 32–33

Case comparisons

• Solar-array degradation costs about 38 kg of payload, while planetary perturbations have secondary impact on final mass but produce trajectory deviations requiring inclusion in mission design.The comparison considers three cases differing in solar-array degradation and planetary perturbation modeling.Source: §5 · pp. 35–36

Trajectory solution

• The optimal solution makes five revolutions around the Sun, unchanged from the optimal solution reported in; changing thruster type does not alter the count, and VIVT force is approximately on the same scale.Source: §5 · p. 36

Propulsion controls

• The optimal thrust profile stays within the specific-impulse-limited thrust envelope, using maximum thrust with minimum Isp on high-thrust arcs and increasingly longer intermediate-Isp intervals farther from the Sun.The intermediate thrust profile also represents the optimal transition between the envelope boundaries.Source: §5 · p. 37

Power dynamics

• As the trajectory evolves, solar-array degradation becomes more pronounced and power drops with increasing solar distance; the switching function remains entirely positive, reaching its least value at t = 2500 days, almost touching S = 0.Source: §5 · pp. 37–38

Modeling implications

• Power-subsystem dynamics affect final delivered mass more significantly than planetary perturbations in this test case.The CSC implementation can also accommodate a maximum solar-array power limit, although that limit is not included in the considered cases.Source: §5 · p. 38

• 2 degrees and 0.5 degrees are the largest in-plane and out-of-plane steering-angle deviations caused solely by secondary-body perturbations.These deviations are reported for the differences between cases 2 and 3.Source: §5 · p. 40, Figure 8 · p. 41
• The largest planets produce the greatest disturbing accelerations, with Earth generating the largest perturbation among the inner planets.Perturbations can change drastically during close encounters with planets.Source: §5 · pp. 40–41, Figure 9 · p. 42
• Final mass alone does not determine the advantage of a variable-Isp, variable-thrust engine, which requires multidisciplinary design optimization including power, propulsion, solar-array, and propellant masses.The comparison must account for each subsystem’s contribution to total spacecraft mass.Source: §5 · p. 41
• CSC makes indirect optimization numerically tractable for optimal-control problems with complex state dynamics and realistic power-system, solar-array degradation, and planetary-perturbation models.The framework enables single-shooting solutions using MATLAB’s fsolve solver.Source: §5 · p. 41
• More sophisticated solvers, integrators, and compiled C++ or Fortran implementations could accelerate the reported computations.The current numerical treatment uses MATLAB’s fsolve solver.Source: §5 · p. 41

6. Conclusion

• CSC modifies the optimal-control formulation to handle multiple propulsion modes and enables indirect methods to address a broader class of trajectory-optimization problems.The framework is proposed for systems whose mode transitions are not known a priori, including VIVT engines.Source: §6 · pp. 42–43
• The VIVT example demonstrates that specific impulse can be modulated for efficiency while optimal transitions between modes must be determined rather than assumed beforehand.The framework is applied to determine transition times for the VIVT engine.Source: §6 · p. 42
• Treating the resulting multi-point boundary-value problems as two-point problems with continuous control inputs makes them simpler to solve.The approach produces continuous controls while retaining the structure needed to represent hybrid operation.Source: §6 · pp. 42–43
• For fuel-optimal heliocentric flight, accounting for power-subsystem dynamics has a more significant impact on final delivered mass than planetary perturbations in this case.This comparison is stated for the considered heliocentric phase and fuel-optimal problem.Source: §6 · p. 43
• For a complex multi-year, multi-revolution VIVT trajectory, the thruster uses maximum, minimum, and intermediate specific-impulse modes with 16 switches.The example includes all-planet gravity perturbations, modified equinoctial orbital elements, and complex-based derivatives.Source: §6 · p. 43
• The framework incorporates shadow- and time-triggered constraints smoothly, including constraints that produce no-thrust arcs and specified time-varying thrust directions.These capabilities extend the method beyond mode switching alone.Source: §6 · p. 43