The study of forces in action, known as kinetics, combined with the properties of motion without reference to forces, or kinematics, forms the study of dynamic motion. In 1775 Leonhard Euler, a Swiss mathematician, began the formal study of kinematics [1]. Following this Baron William Thomson, better known as Lord Kelvin, and a physicist named Guthrie Tait [5] showed that system of dynamic motion maintains energy conservation. A decade later in 1876, Theoretische Kinematik by Franz Reuleaux was republished in English [6] and became a staple in engineering courses from Italy to Russia and the United States by taking practices of engineers and compiling them into a single source document. Prior to the end of the century Thomas Wallace Wright compiled necessary topics for kinematics into a single disquisition [7] that was specifically designed to allow mathematics as simple as basic geometry and trigonometry while still allowing calculus for more compact forms. Five years later in 1903 John Theodore Merz, a German British chemist attempted to unify the kinematic field of thought by publishing a history of European thought on the subject [8].

In a philosophical effort to give credit to the works of Chasle, Hamilton and Rodrigues, Sir Edmund Whittaker published four editions in 1904 [9-12]. In contrast to the approach of Whittaker, Irving Porter Church focused on the engineering approach and published his own guide in 1908 [13]. In an effort to build on Wrights work [14] and reduce the mathematical requirements further Christian Hugo Eduard Study, another German mathematician published reference [15] in 1913. Five years later Andrew Gray argued that the key facet to understand dynamics is kinematics in reference [16]. Gray’s assertion was reemphasized in 1954 when Clifford Truesdell stated, “I cannot too strongly urge that a kinematical result is a result valid forever, no matter how time and fashion may change the ”laws” of physics” [17] and then again in 1967 when G.C. Fox stated “It is my belief that students have difficulty with mechanics because of an inadequate knowledge of kinematics” [18].

In the second half of 1950’s Kane focused on kinetics and expounded on D’Alembert’s version of Lagrange’s equation [20,21]. Computational capabilities drastically approved over the following decade and were utilized by Fang in order to develop a way to mathematically simulate rigid body rotational kinematics [22]. Henderson, in the late 1970’s, further elaborated on the kinematic relationships between rotational sequences as well as the four-dimensional quaternion and direction cosine matrices [23,24]. From this history, the sequential rotation sequences were formed with the two specific sequences being termed the aerospace sequence about repeating axes, sometimes referred to as the “Tait-Bryan angles”, and the orbital sequence using axis-type repetitions, or the “proper Euler angles”.

In the 1980’s, ubiquitous acceptance of two rotation sequences (the so-called “orbital” sequence, and the so–called “aerospace” sequence) appeared in the literature and mature textbooks. [25] In the following decade of the 1990’s renewed interest in improvements was evidenced amidst the ubiquitously accepted paradigms elaborated in a comprehensive survey presented in [26] and evaluation of linearization schemes in reference [27]. Renewed interest in the errors resulting from a chosen kinematic selection was sparked at conference in the United Kingdom at the turn of the century [28] and elaborated for attitude estimation the same year in the Journal of Guidance, Control, and Dynamics. [29] The differences in choices for parameterizing attitude was presented geometrically in 2013. [30] From the available rotation choices, the question may be posed: of the twelve possible sequences, which is the best to convert inertial coordinates to body coordinates? In 2018, Smeresky and Rizzo [31] studied the sequence known as the so–called “aerospace sequence” (also referred to as “3-2-1”) and concluded that the yaw and pitch angles do not accurately depict a rotation about the body’s $z$–axis and $y$–axis respectively. Conversely, only the roll angle, or last angle in the rotation sequence, is accurate in the rotation about the body’s $x$-axis. These facts drive some of the innovations in this manuscript. In 2019, [32] presented the full analytical derivations of the attitude error kinematics equations.

Two figures of merit are used to evaluate the twelve rotation sequences: mean and standard deviation of the representation error to indicate how the tested sequence compares to true roll, pitch and yaw angles and computational burden as shown by total calculation time which shows relative numerical superiority. As an aside, this manuscript separates itself from the work of Smeresky, et. Al, [31] by calculating absolute representation error and precise execution time of the rotation sequence as opposed to tracking errors and end-to-end simulation run-time. Instead, the kinematics portions of this manuscript repeat and validate just–published work. [33]. Meanwhile, the application of Pontryagin’s minimization methods [34] follow the procedures laid out by Sandberg [35] as elaborated in [36] are utilized for the trajectory optimization portions of this manuscript’s proposals, importantly notice time minimization is proposed here, where Sandberg [35] and Raigoza [37] optimized for minimum fuel usage and Raigoza proposed autonomous collision avoidance during reentry.

The results from kinematic analysis show that the so-called aerospace sequence and the reverse (123 sequence) are contrasting in error and computational time with the aerospace sequence being superior between the two. Additionally, the 123 sequence results in a significantly slower computational time than all other rotation sequences. Symmetric rotation sequences, or those that rotate about the same axis first and third in order, average slower computational times than non-symmetric, or rotations involving all three axes, despite requiring the mathematical process and steps to solve for Euler angles. Of all the non-symmetric rotation sequences, the “aerospace” rotation sequence (also called the “321”) was the fastest, while the fastest symmetrical sequence was the “232” demonstrating lower computational requirements.

There are many terms that can be associated with an optimum trajectory; whether the problem is looking for the path of least resistance, the path that provides the best fuel burn, or the path that takes the least amount of time to get from point a to point b. Most of the time, these paths are not direct lines connecting the two points. Knowing the starting and end points turn the selection of the flight path into a boundary value problem solved using ordinary differential equations that can be solved in a multitude of ways. The first, and possibly least complex, method is known as the shooting method where the problem is turned into an initial value problem. A major downfall to the shooting method is initial guesses are required for any unknown initial values. If the provided initial guesses are even minutely off, the entire problem will “implode” and fail to provide a viable solution as the problem is extremely sensitive to initial guess errors [NEW 1]. A second method is known as the collocation method where the boundary value problem is provided a set number of points to form a mesh between the start and end points and this mesh is then continuously analyzed to form an optimal path between the two [NEW2]. As the dynamics and variables used are all known positions, the collocation method eliminates the errors induced by guessing errors but does typically require more computing power than the shooting method. The MATLAB function bvp4c is a built-in program that uses the collocation method. Lastly, multiple companies have developed their own proprietary programs that will solve boundary value problems to a set accuracy. This manuscript uses the program “DIDO” to solve the dynamic optimal problem for a missile set in geosynchronous orbit to a target on the earth’s surface. An advantage to using DIDO is its ability to solve complex dynamic problems without requiring the user to generate the Hamiltonian or any costates or endpoint vectors [NEW3]. By not having to code the Hamiltonian or costate equations, the opportunities for human error in coding is limited to the initial complex dynamic functions.

In the case of a missile delivered from orbit, the amount of fuel used, while a consideration, does not outweigh the necessity to strike the target as fast as possible. A solid rocket motor will burn at a set rate and does not require the oxygen of the atmosphere to burn, and as such can be planned for as a set constant depending on flight time. To validate the results achieved, the Hamiltonian and all costates are developed to concur with the products from DIDO’s program. By manually computing the Hamiltonian and associated equations, natural points are produced that can be analyzed to ensure the validity of the simulation. The results from this manuscript produce a minimum time trajectory from the launch position to chosen target site with a flight time of two hours, thirty-four minutes and forty-six seconds, a fraction of the time required to launch in the opposite direction.

Just last year, [42] proposed and investigated the potential of using suborbital space vehicles for the transportation of cargo to precise delivery locations. Towards achievement of the necessary high–precision of delivery to ensure safety, this manuscript describes selection of kinematics representations to minimize errors and also includes a formal time–minimizing optimization problem seeking rapid delivery capabilities. Three-dimensional (rotational and translational) motion with six degrees of freedom when applied to an object delivered from geostationary orbit. Due to size constraints in modern systems and to minimize necessary computing power, efficient schemes are desired to provide vehicle guidance. Determining the most accurate and computationally efficient kinematic representation is, therefore, essential to producing an effective space delivery capability.

A state variable is one that describes the corresponding state of a dynamic system and can be written as a function of a control variable and the dynamic function of that system. When the state dynamic functions are known, the performance function can be written as Equation 1, where J is the standard performance function, x is a state vector, u is the control variable, E is the endpoint performance (also known as the cost function) and F is the running performance (also known as the running cost).

To minimize the performance function, the state variables are normalized using a common unit of measurement, CU, to ensure that the units make sense. The common unit is used to create a costate identified by λ as shown in Equation 2.

Each state and costate variable do not just consist of a single value but is instead a vector and combines to create the value ω as shown in Equation 3.

Knowing this combination, the units of Equation 1 are the common unit established to form the costate vector since the endpoint performance is composed of the common unit and the running performance is the common unit per time unit integrated.

Pontryagin used a series of equations to determine a candidate optimal solution written as “HAMVET” [36]. These letters represent the sequence: Hamiltonian, Adjoint Equation, Minimization of the Hamiltonian, Value Condition, Evolution Condition, and Transversality. The process is called HAzMAT in [34].

There are a few different methods in which to attempt to determine an optimal solution to the boundary value problem. The first, known as the shooting method, is a simple iterative approach in which the boundary value problem is transformed into an initial value problem using the known initial values and subsequently guessing on the remaining unknown initial values. These values are then iterated forward through the state dynamic equations to see if the final conditions set forth are met. If the required conditions are not met, the initial conditions that were guessed are changed and equations are run again. This process repeats until the conditions do meet at which point the solution becomes a candidate for the optimal solution pending verification and validation. The shooting method is extremely sensitive to small errors in the initial values used for the problem as these errors are then magnified through the dynamic state equations during the iteration.

The second method of solving the boundary value problem is known as the collocation method. The collocation method creates a mesh of a set number of points between the initial condition to the final condition and then adjusts each of these points until all conditions of the BVP are satisfied through the dynamic equations given. The collocation method requires more guesses than the shooting method as any unknown final conditions will need to be guessed in addition to the initial conditions. Additionally, this method requires a great deal more computing power.

Lastly, some companies or corporations create their own software designed to solve the boundary value problem. These programs are typically not either shooting or collocation methods, but instead generate their own algorithm to produce the solution. For this project, a program called “DIDO” [43] is used that takes the problem with the dynamic equations and formats Pontryagin’s Principle into a boundary value problem and then produces a candidate solution. To understand the solution that DIDO gives, it is important that the use works through Pontryagin’s Principle and conduct verification or validation in order to determine if the candidate solution is truly optimal.

For kinetics, this manuscript will discuss the use of Newton’s and Euler’s principles and specifically how they combine to form Chasle’s theorem to form a direct line to explaining the six degrees of freedom in kinetics. Direction cosines and the matrices they form, or quaternions, describe the body of an object in terms of a predefined axis.

In accordance with Euler’s equations of motion, the sums of externally applied torques will angularly accelerate a proportional mass moment of inertia as depicted in Equation 14, where τ is torque, J is the rotation tensor, and ω is the angular velocity.

Similarly, an object of a defined mass will accelerate proportionally to the sum of forces acting on that mass in accordance with Newton’s second law as depicted in Equation 15 where F is the force, m is the defined mass, a is the acceleration, ω is angular velocity, r’ is the position vector, and v’ is velocity vector.

Newton’s law applies in the inertial frame and allows us to express our results in a set of coordinates in another defined frame of reference.

Disturbance forces and torques act on a body in accordance with equations (14) and (15) and are accounted for within the Simulink model for simulation found in Appendix (A).

By utilizing direction cosines, we are able to take a set of coordinates relative to the non-rotating inertial frame and express them in terms of the body frame to account for motion in all three degrees as depicted in Figure 1, where all lines are unit vectors and the thin solid line represents the non-rotating inertial frame and each dash line increasing in weight indicates the first prime, second prime and finally the body frame respectively.

The DCM is labeled in the order that the mathematical calculations are solved, right to left in accordance with algebraic multiplication rules to pre-multiply the total transformation matrix. Equation (16) represents the three mathematical calculations for a “3-2-1” rotation and is then multiplied to generate the total DCM in Equation (17).

A mathematical problem arising from the DCM is the inability to pull values directly from the matrix due to quadrant ambiguity presented by the law of sines and cosines. To solve the quadrant ambiguity problem portions of the DCM have to be individually selected that allow use of the atan2 function to solve for a specific angle by cancelling out other portions, as shown in Equation 18 in which the atan2 function distinguishes between the four standard reference quadrants.

An additional mathematical error from the law of sines and cosines is singularities in the DCM from solutions equaling 1/0. Lastly, the DCM, based on the chosen order of rotation, results in errors in the final result with only the “last” rotation having zero error due to being the only rotation about the body axis. To account for zero error in all three Euler angles, an additional rotation needs to be completed for each angle in the “last” position.

One way to remove the singularity problem presented by DCMs is to use the quaternion where the three-dimensional rotation is executed around a 4th axis; a scaled, stretched and normalized eigenvector, instead of the three standard axes. This requires induction of three imaginary axes of motion as opposed to the standard inertial frame and is shown in Figure (1c) and Equation (19) where q represents the different values of the quaternion and η is the angle that allows a single rotation [1,2]. The quaternion works so long as we enforce the quaternion normalization condition shown in Equation (19). Specifically, Equation 19 shows a direct comparison to the DCM shown in Equation (17).

The guidance of objects is known as the outer loop function and generates commands for the inner loop discussed in section 2.5 by accounting for three degrees of freedom of translational motion governed by Newton’s 2nd Law expressed in Equation 15.

Due to the transport theorem, ω couples the three degrees of freedom of translational and rotational motion. Additionally, equations 14 and 15 are internally coupled within the inertia matrix with ω, causing motion in all three degrees with inputs to only one. A non-linear decoupling controller is placed in a feedback control, or a desired output is placed in a feed forward control, to decouple the nonlinear couple motion as displayed in Equation 20 where u is the controller and ${\omega }_d$ is the desired output angular velocity.

By using Newton’s and Euler’s equations and invoking Chasle’s method, we are able to express six degrees of freedom for a three-dimensional movement of an object. DCMs are then used to describe the position of a body in terms of a defined axes, with quaternions used to eliminate singularities and validate the DCM. From these equations we pull out the true Euler angles.

The dynamic equations of motion describing the vehicle in flight are defined as shown in Equation 21, where r is the distance from the center of the planet in kilometers, θ is the longitude angle, φ is the latitude angle, v is the planet relative velocity, γ is the planet relative flight path angle (or pitch), and ψ is the planet relative heading azimuth (or yaw).