Time optimal control for the deployment of a tethered satellite allowing for a massive tether

The time optimal deployment of a satellite from a space ship is studied. In order to take the mass and lateral oscillations of the tether into account, a discretized model of the tether is build. Applying Pontryagin’s Maximum Principle a time-optimal deployment from a trivial downhanging configuration close to the space ship to another one farther away is computed. It is found, that the obtained solution displays a difficult switching pattern and during the variation of the initial length different kinds of bifurcations occur, leading to discontinuous variations of the optimal solution candidates.


Introduction
The deployment and retrieval of tethered satellites is an important and difficult operation in space missions.
Starting from a position close to the space ship the satellite should be steered to a stable stationary position farther away. The control input should optimally be applied by tension control, that is by a tension force at the outlet of the tether.
Several important goals for the control design have been listed in the review [9], one of which is the reduction of lateral oscillations during the process. Also in the book [1] several aspects in the dynamics and control of tethered systems are investigated.
For safety reasons a PD controller, called Kissel's law, is used commonly during space flights. This controller leads to an exponential decay of the lateral oscillations, but takes very long to complete the mission. In order to find out, how much time could be saved theoretically, a time optimal solution was investigated in [8]. It turned out, that this solution was by an order of magnitude faster than the conventional strategy, but it used a bang-bang control, which could cause unwanted oscillations in the tether. Even, if such control strategies would never be performed, the obtained control patterns should prove useful for smoother control schemes.
Most studies on deployment control, e.g. [11], assume a massless tether, which leads to a very simple set of differential equations and still yields reliable results. In [8] a comparison between the control of the simple model with massless tether and a finite element calculation of a massive tether with the tension force obtained from the simple model is performed and it is This article is dedicated to Prof. Francesco Benedettini. demonstrated, that the control is able to steer the satellite close to the target position in short time.
The focus of the present work is an investigation of possible time optimal controls for a massive tether. For the tether a simple discrete model is created and the governing equations of motion are derived. Unfortunately these equations of motion become extremely lengthy even for quite coarse discretizations and can be treated only numerically and by symbolic algebra. For a reasonable set of parameters a valid optimal control solution could be obtained. This investigation is certainly not intended to provide a control strategy, which should be used in real world applications, but to learn about the typical shape of optimally controlled deployment trajectories. The obtained information should prove useful in designing efficient trajectories in commonly used control methods.
In Sect. 2 the mechanical model is presented and the ingredients for the Lagrangian equations of motion, the kinetic and potential energy, and the control and nonconservative forces acting on the systems are introduced. Also the setup of the optimal control equations using Pontryagin's Maximum Principle is presented. In Sect. 3 the numerically obtained optimal solution is shown and discussed. In order to demonstrate the strong dependence of the derived solution on the initial condition, also the variation of a one-parameter family of optimal solutions with different initial lengths of the tether is presented. In the studied parameter range we observe large variations of the switching points, where the bang-bang control jumps between its extremal values. Further the solution family, which is calculated by a continuation algorithm, displays a non-monotonic behaviour, generating multiple candidates for the optimal control strategy.

Model description
As displayed in Fig. 1, we consider a satellite, which is connected to the main station by an inextensible tether. The heavy main station is assumed to move with constant velocity along a circular Keplerian orbit. The connecting tether is regarded as an inextensible string of constant cross-section and without bending stiffness. Before the satellite is released from the spaceship, the tether of length ' L is reeled on a drum, which is fixed at the space ship; during the deployment process the tether unwinds from the drum. By applying a torque on the drum, a tension force is exerted at the tether, which is used to control the motion of the satellite. For our calculations we neglect the rotational motion of the drum and assume, that the stored part of the tether moves with the main station. Since the tension and velocity of the tether change discontinuously at the outlet, so-called Carnot Energy Loss terms have to be added to the tether equations.
The satellite is considered as point mass moving in the orbital plane in a perfectly circular gravitational field. Since the tether length ' is much smaller than the orbital radius r M , we may simplify the equations of motion by assuming near-field dynamics: The resultant of the centrifugal and gravitational force is proportional to the distance from the Keplerian orbit and points into the vertical direction.
At the start of the deployment process the system rests in the vertically downhanging relative equilibrium with released tether length ' 0 . By varying the tension force at the outlet, the satellite should be steered to another relative equilibrium with released length ' T .

Lumped mass modelling of the massive tether
The partial differential equations for a massive tether are derived in [3]. Since we investigate a control problem, we first discretize the continuous string. Motivated by the sinusoidal shape of the transversal tether oscillations in numerical simulations we first tried to apply a Ritz-Galerkin-type discretization using low order sine waves about the connecting line between main station and the satellite. Due to the complicated dependence of the time varying length of the released tether and the assumed shape of the configuration that approach failed and we then tried a  Fig. 1 Tethered satellite system, consisting of the space ship on a circular orbit, the satellite in the orbital plane, and the connecting tether. The local x-coordinate points into the flight direction, and y points toward the center of the earth lumped-mass approach, as already proposed in [4] for elastic beams and sketched in Fig. 2: The current shape of the deployed tether is regarded as a chain of N straight pendula of length hðtÞ ¼ 'ðtÞ=N and mass lhðtÞ, where l denotes the mass of the tether per unit length. The first pendulum is attached to the outlet at the main station; the angle between the ith pendulum and the local vertical direction y is denoted by u i . Of course, the individual pendula in this model should not be regarded as rigid rods, but just as quite coarse parametrization of the currently deployed tether configuration. When the tether becomes longer, the pendula grow accordingly. Since the tether is released from the main station, the material points of the tether move along the pendula.
The position of the tether at the relative length n along the ith pendulum in the local frame is given by for the position of the satellite we get At the tether outlet the control force F is acting against the direction of the string The virtual work of this control force is therefore given by The gravity potential of a mass m at r ¼ ðx; yÞ T in the near-field dynamics is obtained as the second order approximation of the gravity field with X denoting the orbital angular velocity. From this expression we find the augmented potential With (4) the gravity potential of the system is given by In order to calculate the kinetic energy, the velocities of a point r ¼ ðx; yÞ T in the local frame have to be expressed in the inertial frame using the relation The kinetic energy of the satellite is given by For the kinetic energy of the tether we first have to calculate the velocity of the particles, which are located at the position s along the tether, expressed in the local frame. Due to the variation of the tether length that velocity differs from dr i =dt, because the particles move along the pendulum chain with relative speed _ 'ð1 À s='Þ. The velocity in the local frame of the particles, which are located at the positions r i ðnÞ of the pendula, is therefore given by Now the kinetic energy of the released tether is given by The complete kinetic energy is then obtained as Finally we have to consider the Carnot energy loss terms at the drums. According to [3] the instant acceleration of the tether leads to a nonconservative force which can be easily derived by comparing the Lagrangian equations of motion for an inextensible string, which is initially folded to a negligible volume and pulled by a horizontal force F, and the momentum balance equation for this system: Assume that the unfolded length of the string is denoted by x(t), then its kinetic energy is given by and with the virtual work dA ¼ Fdx the generalized force is obtained as Q ¼ F. Due to the absence of further forces the Lagrangian equations of motion are given by d dt The momentum of the unfolded string is given by therefore the momentum balance _ I ¼ F yields the equation which differs from (12) by the term l _ x 2 =2. Together with (3) we find the generalized work for the length variable There occur no non-conservative forces for the angular variables Before calculating the Lagrangian differential equations, the variables are scaled according to It should be noted, that the rescaled time t now corresponds to the ''orbital time'' for the main station: A revolution around the earth takes 2p time units. The parameter l denotes the ratio of the deployed tether mass to the mass of the subsatellite. The dynamics of the tethered satellite is then given by the second order system d dt for the degrees of freedom q ¼ ð'; u 1 ; . . .; u N Þ T . Due to the variation of the tether length the symbolic calculation of these equation is very elaborate even for N ¼ 2 and was carried out using symbolic algebra [10].
With the symmetric mass-matrix M with entries with the abbreviations

Optimal control problem
We search for the tension force F(t), which steers the tethered satellite from the straight downhanging configuration with initial length 'ð0Þ ¼ ' 0 to the downhanging configuration with final length 'ðTÞ ¼ ' T [ ' 0 in shortest time. The tension force F(t) has to satisfy the inequality constraints In order to avoid a slack tether and too large deviations from the straight configuration, we choose F min [ 0.
The upper boundary F max should be larger than the required force F s to hold the system in the final equilibrium, but not too large for safety reasons. In our numerical experiments we choose the scaled values The static equilibrium force is given by In order to derive the differential equations for the Optimal Control problem [5], the system (18) is rewritten as explicit first order system Since we are looking for the fastest solution, the utility function is given by Introducing the adjoint variables p i , we build the Hamilton function The optimal tension force F H is determined by the Maximum Principle Since the tension force F occurs linearly in f ðx; FÞ, the switching function S is given by The singular case, that the switching function S vanishes in an interval, was never observed in our calculations.
The adjoint variables satisfy the differential equations The required boundary conditions read Since the shortest time T isn't specified, the additional boundary condition HðxðTÞ; pðTÞ; FðTÞÞ ¼ 0 ð32Þ has to be satisfied.
Remark The conversion of the Lagrange equations (18) to the explicit system involves the symbolic inversion of the mass matrix, which is already very cumbersome for N ¼ 2. For the adjoint differential equations (29) the partial derivatives of that inverse w.r.t. the state variables would be required. In order to avoid these tedious calculations, one can exploit the identity If the state equations (24) are written in the implicit form where f ðx; FÞ ¼ M À1 ðxÞgðx; FÞ, the corresponding entries in the adjoint equations become where P ¼ M À1 ðxÞp. Since already the calculation of _ x ¼ M À1 ðxÞgðx; FÞ requires a factorization of M, the additional step to solve the equation can be carried out with the same factorization. The derivatives of M w.r.t. x i are much easier to calculate than those of its inverse.

Numerical solution of the boundary value problem
Due to their complicated form the differential equations can only be solved numerically. We use the Fortran solver Boundsco [6], which solves two-point boundary value problems (BVPs) with switching conditions using a multiple shooting method. Because the endpoints of the integration interval have to be specified in advance, we map the unknown integration interval to the unit interval by introducing a scaled time t ¼ Ts; withs 2 ½0; 1: The differential equations (24) For the unknown variable T we introduce another trivial differential equation dT=ds ¼ 0. Self-evidently the boundary conditions (31) and (32) apply at the right endpoint s ¼ 1.
The 4N þ 5 boundary conditions (30), (31), (32) and the differential equations for the 4N þ 5 variables x i , p i and T form the boundary value problem. Since the Hamiltonian depends linearly on the tension force F, we expect a bang-bang control, with the tension force F(t) jumping between its limiting values. At the switching points s i the switching function S(t) vanishes In Boundsco the switching points are treated as internal additional variables, which are determined by (37).
In order to find a solution with bang-bang control, the user has to provide the number of switching points and a proper initial guess for their values. Furthermore the switching sequence has to be fixed at the entry to Boundsco.
Since the differential equations (24) are independent of the co-state variables and the boundary conditions (30) and (31) involve only the state variables x i , we expect to need at least 2N þ 1 switching points to satisfy all prescribed boundary conditions and obtain a unique solution for the co-state variables: For an arbitrarily prescribed time interval T the solution of the fully specified initial value problem (30) for the x i will generically satisfy not any boundary condition (31) at t ¼ T. By varying T we may expect to satisfy one terminal condition. In order to satisfy the remaining boundary conditions, 2N þ 1 further adjustable parameters are needed.
On the other hand more boundary conditions for the 2N þ 2 co-state variables are required: The optimality condition (32) provides one linear constraint for pðTÞ, the remaining 2N þ 1 equations are provided by the switching conditions, which involve the costate variables.

Numerical pathfollowing method
Although Boundsco is a very robust and efficient solver for BVPs, we still have to provide sufficiently good initial guesses for the state and co-state variables and especially for the switching points. If the initial guesses for the switching points are too inaccurate, the solver will almost certainly fail to solve the problem.
A very efficient method for providing good initial guesses is the so-called homotopy or continuation method [7]: If the BVP depends on some parameter k and a solution y 0 for a certain value k 0 is known, that solution will usually be a good guess for a neighbouring parameter value k ¼ k 0 þ dk. With this method one-parameter families of trajectories can be obtained and ''better'' solutions can be found.
That method was successfully applied in [8] to obtain a valid solution for the deployment problem with a massless tether: In the first run the end conditions, the three switching points, and the time interval T were prescribed and a solution was obtained, which satisfied none of the remaining boundary conditions. Starting with this trajectory and varying the time interval T, a solution could be found, which satisfied one initial condition.
In the next loop that initial condition was enforced and T was allowed to vary freely. By varying one of the switching points, another initial condition could be satisfied. After another loop all homogeneous initial conditions could be satisfied and some 'ð0Þ quite close to the final value ' T was obtained. Applying the homotopy method the initial value 'ð0Þ could easily be decreased down to ' 0 ¼ 0:05.

Numerical results for the massive tether
with N ¼ 2 In order to study the behaviour of the tethered satellite with a lightweight tether we started with the simple case N ¼ 2 and l ¼ 0:2. Depending on the tether material and deployed tether length the value of l is usually much smaller (about 0.01), but the relatively large value causes slower transversal oscillations. After a series of homotopy loops we arrived at a solution with ' 0 ¼ 0:1 and 7 switching points, where the control F changed between the F min and F max , the trajectory in the co-rotating orbital frame is displayed in Fig. 3. The parts of the path, where the tension force F ¼ F max , is displayed by a heavy line, the switching points are indicated by black circles. The trajectory shows a quite strong deviation from the local vertical direction.
The trajectory of the subsatellite is quite similar to the obtained solution for the massless tether (see [8]), but there occur three further short intervals of maximum tension, which straighten the tether configuration before the final braking period. One of these pulling periods occurs right at the beginning of the manoeuvre.
For the same boundary conditions the control with a massless tether requires only 3 switching points: Initially the satellite is released (almost) freely. Since it moves approximately on a Keplerian ellipse with smaller major axis than the radius of the main station, it moves ahead in the orbital plane. In order to make it move back to the local vertical, maximum tension is applied in the second phase. After another free flight period it is steered to the target position using again maximum tension. If the same control were applied to the light-weight tether, quite heavy lateral oscillations would persist beyond the control interval.
The temporal behaviour of the tension force F(t) and the switching function S(t) is displayed in Fig. 4. It shows a quite strong oscillation of the switching function. Three pairs of switching points close to t ¼ 0:2, t ¼ 2:2, and t ¼ 2:55 are quite narrow, such that the influence of the short change in F is hardly visible in the following figures.
The switching behaviour in the time interval [1.7, 2.8] is displayed in Fig. 5.
In Figs. 6 and 7 the evolution of d'=dt and of the angular variables u i ðtÞ is displayed. During maximum tension the length change rate d'=dt decreases. From Fig. 7 it is visible, that the transversal oscillations are quite heavy during the first long pulling period. When the final braking period starts, the tether is almost

Variation of the calculated solution for varying initial conditions
The previously displayed solutions correspond to a fixed initial value ' 0 ¼ 0:1. (According to (17) the final length is scaled to ' T ¼ 1.) By applying a continuation strategy solutions with different values for ' 0 can be computed. As displayed in Fig. 8, ' 0 need not vary monotonically during that continuation, but may also turn around and lead to different solutions for the same boundary conditions: In Fig. 8 the variation of a selection of interior switching points s i over the initial length ' 0 is displayed. Starting with a solution with ' 0 % 0:14 and 7 switches in the control force and decreasing ' 0 , a pair of switching points s 4 and s 5 has to be inserted at P 1 : ' 0 ¼ 0:133. Decreasing the new solution (''A'' in Fig. 8) with 9 switches further down to P 1 : ' 0 ¼ 0:0815, again a pair s 6 and s 7 has to be inserted. Along the solution branch ''B'' with 11 switches the initial length ' 0 increases (dashed line) until P 3 : ' 0 ¼ 0:127, where the switching points s 4 and s 5 coalesce and disappear. Along the new solution branch ''C'' with 9 switches the value ' 0 decreases again and extends below ' 0 ¼ 0:08. After traversing the loop a pair of switching points has changed its position significantly: The lower pair s 4;5 , which is created at P 1 , vanishes at P 3 along the branch ''B''. The higher pair s 6;7 is born at the start of branch ''B'' at P 2 and persists along branch ''C''. The remaining switching points, e.g. s 3 and s 8 , vary continuously along the branches. All three branches of solutions exist for ' 0 ¼ 0:1; in Fig. 9 the switching function S(t) is displayed for the three solutions with initial length ' 0 ¼ 0:1. It looks rather similar, but solution ''A'' has two zeroes close to t ¼ 1:8, while C has two zeroes at about t ¼ 2:2. The intermediate solution B has all four zeroes.
Since the difference between the required times T is very small, the quantity T À T M is displayed in Fig. 10 for the three solution branches, where T M ¼ ðT A þ T C Þ=2. The swallowtail-shaped cost function intersects itself at ' 0 % 0:12 For ' [ 0:12 solution ''A'' is more efficient, whereas for '\0:12 ''C'' performs better. The intermediate solution ''B'' with 11 switching points takes always longer than the other candidates and serves only to connect the two branches.
Throughout the common domain of existence of the two branches ''A'' and ''C'' both solutions are local minima. At the intersection of the cost functions T A and T C at ' 0 % 0:12 a so-called ''Maxwell Catastrophe'' [2] occurs, where the optimal solution candidate changes discontinuously. Of course, the difference between T A and T B is extremely small (less than a second for a manoeuvre lasting almost an hour) and plays no practical role. Nevertheless, the birth and death of the switching points needs to be taken into account during the variation of initial conditions and parameters in the system.

Another bifurcation scenario
An almost similar bifurcation occurs for slightly larger values of ' 0 : Close to ' 0 ¼ 0:18 the homotopy branch shows a non-monotonic variation of ' 0 , as depicted in Fig. 11. The distorted cost function T þ 4' 0 is displayed in Fig. 12. It shows again the shape of a swallow-tail, with two different solutions ''A'' and ''C'' with the same cost function at the intersection point. The intermediate solution ''B'' at the upper branch again takes longer and just connects both efficient branches.
The switching function S(t) is displayed in Fig. 13: While the left switching points vary rapidly during the loop, the remaining ones are almost not affected.
The behaviour of _ 'ðtÞ and u 1 ðtÞ À u 2 ðtÞ for the three solution candidates ''A'', ''B'', and ''C'' is displayed in Figs. 14 and 15: Again mainly the phase after the first pulling period is affected strongly. Since the influence of the variation of the first switching points on the variation of u 1 þ u 2 and on the trajectory of the sub-satellite is very small, these quantities are not drawn. Contrary to the first bifurcation scenario, where the non-monotonic variation of ' 0 along the continuation path was caused by the birth and deletion of switching points, in this case the switching points vary smoothly. Both cases demonstrate, that it is important not to calculate only one solution for the intended boundary conditions, but to explore also the vicinity of that solution in order to maybe detect more efficient solutions.

Conclusions
For a tethered satellite system with a massive tether a simple discretized model has been derived. For this  model a time-optimal solution has been obtained numerically. The overall solution shape looks quite similar to that of the corresponding model with zero tether mass. But there are some significant differences: • The system with a massive tether requires significantly more switching points for its bang-bang solution than the massless one. Some of these switching points are quite likely not relevant for the practical applications, because they are extremely close and show almost no influence on the performance. • As expected, the strongest oscillation of the tether occurs during the pulling face, when the curved tether becomes stretched. Since there are no damping forces present, these fast oscillations can only be extinguished by releasing the tension force at the right phase. Before starting the final braking stage the lateral oscillations should already be extinguished. • The numerical detection of the obtained solution was quite involved. According to the author's estimate the calculation of optimal trajectories for more realistic parameters (finer discretization, lighter tether, smaller ratio ' 0 =' T ) would require a huge amount of work and it is quite unlikely that a robust and efficient feedback law could be designed to take care of perturbations to the computed solution. • The occurrence of turning points in the continuation procedure-which also have to be expected, if some parameters in the model are variedindicates, that it is important to explore the behaviour of the system for larger ranges of parameters and initial conditions. Otherwise one might end up at some sub-optimal candidate. In the author's experience such loops occur frequently, when the overall structure of the solution changes. Looking closely at the change in the solution after passing the turning points might provide valuable insight in the control strategy.