Determining Optimal Parameters in Magnetic Spacecraft Stabilization via Attitude Feedback Renato Bruni1,a) and Fabio Celani2,b) 1

Dep. of Computer, Control, and Management Eng., Sapienza University of Rome, Via Ariosto 25, 00185 Roma, Italy. 2 School of Aerospace Engineering, Sapienza University of Rome, Via Salaria 851, 00138 Roma, Italy. a)

Corresponding author: [email protected] b) [email protected]

Abstract. The attitude control of a spacecraft using magnetorquers can be achieved by a feedback control law which has four design parameters. However, the practical determination of appropriate values for these parameters is a critical open issue. We propose here an innovative systematic approach for finding these values: they should be those that minimize the convergence time to the desired attitude. This a particularly difficult optimization problem, for several reasons: 1) such time cannot be expressed in analytical form as a function of parameters and initial conditions; 2) design parameters may range over very wide intervals; 3) convergence time depends also on the initial conditions of the spacecraft, which are not known in advance. To overcome these difficulties, we present a solution approach based on derivative-free optimization. These algorithms do not need to write analytically the objective function: they only need to compute it in a number of points. We also propose a fast probing technique to identify which regions of the search space have to be explored densely. Finally, we formulate a min-max model to find robust parameters, namely design parameters that minimize convergence time under the worst initial conditions. Results are very promising.

INTRODUCTION The attitude control of a spacecraft that uses magnetorquers as torque actuators is a very important task in astronautics. Many control laws have been designed for this task. A survey of various approaches is in [6]; in particular [1] proposes a feedback control law that, besides measures of the geomagnetic field, requires measures of attitude only. This work shows that attitude stabilization is achieved when the design parameters have certain properties (e.g., they are positive). However, the practical determination of appropriate values for the design parameters is a critical open issue. In this work we propose an innovative systematic approach for determining them: they should be those that minimize the time needed for convergence to the desired attitude. Specifically, we measure this time by considering the Integral Time Absolute Error (ITAE) of the attitude . This turns out to be a particularly difficult optimization problem, for several reasons. The relation between the above ITAE and the design parameters cannot be expressed in analytical form, but only sampled by using software simulations. For this reason, we propose a solution approach based on the use of derivative-free optimization algorithms. These algorithms use no first order information on the objective function. In practice, they work without the need for analytically writing an objective function; they only need to compute it in a number of points by using the above simulations. Moreover, design parameters may range over very wide real intervals, hence the search space of the optimization problem is numerically huge. Consequently, we propose the use of a fast probing technique to identify the ‘promising’ region(s) of the search space, and then we explore them thoroughly. Finally, ITAE also depends on initial conditions. This implies that control parameters that are efficient for specific initial conditions are not guaranteed to remain efficient for different initial conditions. Clearly, there exist a large variety of possible initial conditions, and they are not known a priori. Thus, we initially solve the problem for a particularly meaningful value of initial conditions. Subsequently, we compute the values of the parameters that minimize the ITAE obtained under the worst initial conditions. Such worst initial conditions are not defined in general, but they in turn depend on the adopted design parameters. Therefore, we formulate a min-max problem, whose solutions are robust optimal values for the design parameters. This problem is quite difficult, since the solution of the main minimization problem (upper-level)

needs the solution of a maximization problem (lower-level) at every evaluation of its objective function, and a decomposition is not possible. Many real-world problems share this features, in particular numerous engineering design problems. The proposed approach uses no peculiarities of the considered case and can in principle be extended to other problems with the same structure.

CONTROL ALGORITHM AND OPTIMAL DESIGN PARAMETERS To describe the attitude dynamics of an Earth-orbiting rigid spacecraft and represent the geomagnetic field, we use the standard Geocentric Inertial frame Fi and fix a spacecraft body frame Fb . The goal is aligning Fb to Fi ; thus, the focus will be on the relative kinematics and dynamics of the satellite with respect to the inertial frame. Let q = [q1 q2 q3 q4 ]T = [qTv q4 ]T be the unit quaternion representing rotation of Fb with respect to Fi , and let C(q) be the corresponding attitude matrix. Superscript × applied to any a ∈ R3 denotes the skew symmetric matrix that allows to express the cross product a × b as as the matrix multiplication a× b (as standard in this field, see e.g. [1]). The relative attitude kinematics is given by q˙ = W(q)ω where ω ∈ R3 is the angular rate of Fb w.r.t. Fi resolved in Fb and " # 1 q4 I + q×v W(q) = . (1) −qTv 2 The attitude dynamics in body frame can be expressed by J ω ˙ = −ω× Jω + T where J ∈ R3×3 is the spacecraft inertia matrix, and T is the control torque expressed in Fb . The spacecraft is equipped with three magnetic coils aligned with the Fb axes which generate the magnetic attitude control torque T = mcoils × Bb = −Bb× mcoils where mcoils ∈ R3 is the vector of magnetic dipole moments for the three coils, and Bb is the geomagnetic field at spacecraft expressed in body frame Fb . Let Bi be the geomagnetic field at spacecraft expressed in inertial frame Fi . Note that Bi varies with time at least because of the spacecraft’s motion along the orbit. Then Bb (q, t) = C(q)Bi (t) which shows explicitly the dependence of Bb on both q and t. By grouping together the previous equations the following system is obtained q˙ = Jω ˙ =

W(q)ω −ω× Jω − Bb (q, t)× mcoils

(2)

Since C(q) = I for q = [qTv q4 ]T = ±q¯ where q¯ = [0 0 0 1]T , then the goal is designing control strategies for mcoils so that qv → 0 and ω → 0. The following stabilizing control law is proposed in [1] δ˙ = mcoils

=

α(q − λδ) −m?coils sat

1

m?coils

  Bb×  2 k1 qv + k2 αλW(q)T (q − λδ)

(3)

in which δ ∈ R4 , “sat” denotes the standard saturation function, and m?coils is the saturation limit on each magnetic dipole moment. Note that the previous equation describes an attitude feedback since, besides measures of Bb , it requires measure of attitude q only, and not of attitude rate ω, too. As shown in [1], selecting k1 > 0, k2 > 0, 1 α > 0, λ > 0, and choosing  > 0 small enough, local exponentially stability of equilibrium (q, ω, δ) = (q, ¯ 0, λ q) ¯ is achieved for the closed-loop system (2) (3) if the orbit’s inclination is not too low. However, there are no indications for choosing the feedback parameters k1 , k2 , α, λ, and the scaling factor . We only know that k1 , k2 , α, λ have to be > 0, and  has to be > 0 and smaller than an upper bound  ∗ > 0 which is very difficult to compute. In practice they are mostly determined by a trial-and-error approach, which is exceedingly time-consuming and not systematic. In this work we propose the following approach to determine the feedback parameters. Introducing κ1 =  2 k1 , κ2 = k1 , and β = λ, feedback (3) can be expressed in terms of the four design parameters κ1 , κ2 , α, and β . Now, having set the spacecraft’s initial conditions to specific values, one can determine the values of κ1 , κ2 , α, and β that minimize the settling time t s . However, the objective function t s is not continuous with respect to the design parameters, and this introduces numerical difficulties in solving the optimization problem. Thus, we consider, as alternative objective function, the so called Integral Time Absolute Error (ITAE), denoted by Γ: Z Tf Γ= tkqv (t)kdt (4) 0

where T f is a time chosen sufficiently large. The ITAE has the advantage of being continuous with respect to the design parameters; moreover, minimizing the ITAE leads in general to suboptimal solutions towards the settling time

objective (see [2]). To formulate our optimization problem, we define physically reasonable upper bounds κb1 , κb2 , b α, b β for the design parameters, obtaining a feasible set K = {(κ1 , κ2 , α, β) : 0 ≤ κ1 ≤ κb1 , 0 ≤ κ2 ≤ κb2 , 0 ≤ α ≤ b α, 0 ≤ β ≤ b β}. Now, our optimization problem is: min Γ. (5) (κ1 , κ2 , α, β) ∈ K Given specific initial conditions, problem (5) can be solved to optimality by a suitable use of derivative-free technique described in the following section. However, when changing the initial conditions, that solution may be no longer optimal. Since many different initial conditions may occur in practice, a robust approach would be to find the optimal solution to problem (5) under the worst initial condition. Such a worst case optimization is widely used in similar scenarios, because by doing so we can provide an optimized bound on the objective value notwithstanding the uncertainty regarding the initial conditions. However the worst initial condition is not a priori computable, since it depends on the chosen values of κ1 , κ2 , α, β, so the problem cannot be decomposed and should be solved as a whole. The initial conditions are given by q0 = q(0), ω0 = ω(0), 0 ≤ ψ < 2π which is the argument of the spacecraft at time t = 0, and 0 ≤ α0 < 2π which denotes the right ascension of the Earth’s magnetic dipole at time t = 0 (see [1]). By expressing the set kq0v k ≤ 1 in spherical coordinates (ρ, φ, θ), the set of values describing the initial conditions can be expressed as the following hyperrectangle  d d d H = (ρ, φ, θ, ω0 , ψ, α0 ) : 0 ≤ ρ ≤ 1, 0 ≤ φ ≤ 2π, 0 ≤ θ ≤ π, |ω01 | ≤ ω 01 , |ω02 | ≤ ω 02 , |ω03 | ≤ ω 03 , 0 ≤ ψ ≤ 2π, 0 ≤ α0 ≤ 2π Note that H permits any possible initial attitude and any possible value for ψ and α0 ; it only limits the magnitude of the initial angular rate. Now, minimizing ITAE under the worst initial conditions corresponds to the min-max problem: min (κ1 , κ2 , α, β) ∈ K

max Γ. (ρ, φ, θ, ω0 , ψ, α0 ) ∈ H

(6)

PROPOSED DERIVATIVE-FREE SOLUTION APPROACH Problem (5) may in general be tackled by a global derivative-free optimization algorithm of the type of Direct [4]. Those methods work without the need for analytically writing the objective function; they only need to compute it in a number of points by using simulations. Roughly speaking, this approach is a modification of the standard Lipschitzian approach that eliminates the need to specify a Lipschitz constant. Due to the so-called everywhere dense property, such an algorithm converges to a global optimum of the function if the sampling is dense enough. After considerations on the physical problem, the upper bounds κb1 , κb2 , b α, b β can be set, usually at very high values. However, by doing so, a sufficiently dense exploration of the feasible set K requires a number of function evaluations such that the corresponding computational time would be impracticable. Conversely, a sampling that uses a practically affordable number of function evaluations does not reach solutions sensibly better than random solutions. Indeed, our problem has a large number of local minima, which strongly affect the difficulty of the overall optimization framework. Therefore, we develop a fast probing technique to identify the ‘promising’ regions of the search space, and then we can explore such regions densely. The probing technique is based on the use of the local derivative-free optimization algorithm Sdbox. This algorithm was initially proposed in [5] as a globally convergent algorithm for the minimization of a continuously differentiable function, but it can be practically used to optimize different types of functions as a good compromise between efficiency and convergence properties (see also [3]). It needs no information on the first order derivatives, because a good feasible descent direction is determined by investigating the local behavior of the objective. Most interestingly for our case, given a starting point (κ10 , κ20 , α0 , β0 ) = x0 , this algorithm is aimed at quickly finding good solutions in the vicinity of x0 . A (condensed) description of the proposed approach is: 1. Grid partition the whole feasible set K into a collection of subsets {K1 , . . . , K p }. 2. For each K j , compute the value f j of the solution obtained by maxeval iterations of Sdbox in K j starting from its central point. This upper bound on the value of the best solution in K j constitutes our ‘evaluation’ of K j . 3. Take a number maxsubsets of subsets corresponding to the smallest of the above f j values. 4. Take the region given by the union of those subsets, and ‘convexify’ it by including also the additional subsets required to make it an hyperrectangle K ∗ . Note that K ∗ ⊂ K experimentally is orders of magnitude smaller than K. 5. Switch to Direct algorithm to continue the search in K ∗ . This search can now be dense using reasonable time (see next section), and it produces a solution x∗ . 6. Try to improve x∗ by using maxpost iterations of Sdbox, finally obtaining a solution x∗∗ to problem (5).

On the other hand, problem (6) is composed of an upper level minimization problem and a lower-level maximization one. By renaming the set of variables (κ1 , κ2 , α, β) as x and the set of initial conditions (ρ, φ, θ, ω0 , ψ, α0 ) as y, the first problem has the following form. min g(x) x∈K We solve it by an external loop applying the Direct global strategy to a set K ∗ computed as before. Function g is such that its value on the generic point x¯ is given by the solution of the lower-level problem: g( x¯) = max Γ( x¯, y). y∈H Thus, one lower-level problem must be solved at each iteration of the external loop, and reducing its computation time to a few seconds is crucial. By using our knowledge of the physical problem, we are able to suppose that good solutions of the maximization problem are in the vicinity of extremal values of angular velocity. Hence, we solve the lower-level problem in a nested loop by using the local strategy Sdbox with multi-start from different starting points, in particular from all the eight combinations of extreme angular velocity, obtaining a solution xR∗ to problem (6).

COMPUTATIONAL RESULTS We apply our approach to solve the case study presented in [1]. The spacecraft’s inertia matrix is J = diag[27, 17, 25] kg m2 ; the saturation level for each magnetic dipole moment is m∗coils = 10A m2 . The inclination of the orbit is incl = 87◦ , and the orbit’s altitude is 450 km; the value Ω of Right Ascension of the Ascending Node is 0. Upper bounds κb1 , κb2 , b α, b β are set at (109 , 109 , 104 , 10−3 ). To solve problem (5) for initial condition (ρ, φ, θ, ω0 , ψ, α0 ) = (0, 0, 0, 0.02, 0.02, −0.03, 0.9416, 4.5392) we initially used Direct on the feasible set K. After 20000 iterations and 100961 secs. (about 28 hours), it gave the solution: κ1 = 833333333.333, κ2 = 944444444.444, α = 8611.11111111, β = 0.0005 whose value is ITAE=1.188×109 , that is not much different from random solutions, whose vast majority have values around 1.2 ×109 . We therefore used the described probing technique to determine K ∗ = {(κ1 , κ2 , α, β) : 200000 ≤ κ1 ≤ 260000, 200000000 ≤ κ2 ≤ 310000000, 0 ≤ α ≤ 100, 0 ≤ β ≤ 0.001}. The determination of K ∗ requires 400000 iterations of Sdbox and 201500 secs of computations (about 56 hours). Note that the time required by each simulation is extremely variable. Now, Direct used with the feasible set K ∗ after 3007 iterations and 15 secs. gives the solution: κ1 = 200781.893004, κ2 = 226268861.454, α = 39.3004115226, β = 0.0005 whose value is ITAE=9.072×106 . This solution is orders of magnitude better than before and it is also considerably better than the best solution obtained by trial-and-error in days of work (see [1]), that has ITAE=3.7 ×107 . We then used 5000 additional iterations of Sdbox in 11 secs. to improve the objective in ITAE=9.062×106 with the solution: κ1 = 200781.893004, κ2 = 226268812.597, α = 39.3004115226, β = 0.00050000002366 On the other hand, to solve problem (6) using the above nested loop scheme, with 3005 iterations of Direct for the external loop and 512 iterations of Sdbox in each internal loop, in 26128 secs. (about 7 hours) we have the solution: κ1 = 227777.777778, κ2 = 294444444.444, α = 83.3333333333, β = 0.00061095869532 that has value ITAE= 2.176×107 for the above initial conditions but it is a robust solution: by varying y ∈ H the worst value that is obtained is ITAE=1.357×108 , that is still considerably better than average solutions. More details will appear on an extended version of this work.

REFERENCES [1] [2] [3] [4] [5] [6]

F. Celani, Acta Astronautica 107, 87-96 (2015). D. Graham and R.C. Lawthrop, Transactions of the AIEE 72, 273-288 (1953). T.G. Kolda, R.M. Lewis, V. Torczon, SIAM Review 45(3), 385-482 (2003). D.R. Jones, “DIRECT global optimization”, in Encyclopedia of Optimization, edited by C.A. Floudas and P.M. Pardalos (Springer, Berlin, 2009), 725-735. S. Lucidi and M. Sciandrone, Computational Optimization and Applications 21, 119-142 (2002). E. Silani and M. Lovera, Control Engineering Practice 13(3), 357-371 (2005).