Open access peer-reviewed chapter

Parameter Optimization for Spacecraft Attitude Stabilization Using Magnetorquers

Written By

Renato Bruni and Fabio Celani

Submitted: February 19th, 2019 Reviewed: August 16th, 2019 Published: January 15th, 2020

DOI: 10.5772/intechopen.89197

Chapter metrics overview

524 Chapter Downloads

View Full Metrics


The attitude stabilization of a spacecraft that uses magnetorquers as torque actuators is a very important task. Depending on the availability of angular rate sensors on the spacecraft, control laws can be designed either by using measurements of both attitude and attitude rate or by using measurements of attitude only. The parameters of both types of control laws are usually determined by means of a simple trial-and-error approach. Evidently, such an approach has several drawbacks. This chapter describes recently developed systematic approaches for determining the parameters using derivative-free optimization techniques. These approaches allow to find the parameter values that minimize the settling time of the attitude error or an indirect measure of this error. However, such cost indices depend also on initial conditions of the spacecraft, which are not known in advance. Thus, a min-max optimization problem is formulated, whose solution provides values of the parameters minimizing the chosen cost index under the worst initial conditions. The chapter also provides numerical results showing the effectiveness of the described approaches.


  • derivative-free optimization
  • spacecraft attitude control
  • robust optimization
  • min-max problems
  • magnetic actuators

1. Introduction

A magnetorquer(or magnetic torquer) is a torque actuator widely used for attitude control in satellites, especially those flying in low Earth orbits. The magnetorquer generates a magnetic dipole that interacts with the Earth magnetic field, thus generating a torque used to control the spacecraft attitude. The control torque generated by magnetorquers is constrained to belong to a plane orthogonal to the Earth magnetic field; hence, magnetorquers may be supported by additional torque actuators to achieve full three-axis control (see [1], Chapter 7). However, attitude control systems using only magnetorquers represent a feasible option especially for low-cost satellites or for satellites with a fault in the main attitude control system. Therefore, attitude control of spacecraft using magnetorquers is a very important topic in aerospace engineering.

Many control laws have been designed in such a setting, and an overview can be found in [2]. In particular, Celani [3] shows that attitude stabilization using only magnetorquers can be achieved by proportional derivative-like (PD-like) control which requires measurements of both attitude and attitude rate. The same work shows that attitude stabilization can also be obtained by using attitude-only feedback, which has the advantage of not requiring the installation of angular rate sensors (rate gyros), thus saving in cost, volume, and weight. Those control laws contain parameters. Thus, numerical values must be assigned to those parameters to practically implement the laws. A common way to find those values is by using a so-called trial-and-error approach, consisting in trying several values and rejecting those which are not acceptable, and eventually choosing the best among the acceptable ones. This approach is affected by two important drawbacks: (a) it is often very time-consuming, and (b) it is not systematic. Thus, even if an obtained solution is satisfactory, it is not known if protracting the search could lead to better solutions and which would be the magnitude of the possible improvement.

On the other hand, a systematic approach for determining the parameters should aim at finding the values which minimize the settling time of the attitude error. Such an approach has been recently proposed in [4, 5, 6]. However, this is not an easy objective to pursue. The settling time depends not only on the parameters but also on the initial conditions of the spacecraft. To overcome this issue, the above works propose to compute the values of the parameters that minimize the settling time obtained under the worst initial conditions, so as to provide averagely good results and to set as bound the worst-case behavior. Hence, in the above works, the problem is modeled as a min-max problem, and the obtained parameters’ values are called robust optimalvalues.

This min-max problem is considerably challenging, since solving the main minimization problem (upper-level) requires solving a maximization problem (lower-level) at every evaluation of its objective function. A decomposition is not possible because the worst initial conditions are not determined in general: they, in turn, depend on the adopted parameters. Optimizing this min-max problem involves other two mathematical issues:

  1. Settling time cannot be expressed as an analytical function of parameters and initial conditions; therefore the specification of an explicit optimization model is not possible.

  2. Settling time is not even continuous with respect to both parameters and initial conditions, in the sense that small variations of the latter may result in a substantial gap in the variation of the settling time.

To overcome these issues, the approach proposed in [4, 5, 6] relies on the use of derivative-free optimization algorithms as building blocks. These techniques do not need first-order information on the objective function nor do they need its analytical expression. They only need to compute the objective function over a number of points by means of simulations.

For the design of the specific optimization algorithm, a distinction must be done between spacecraft equipped with angular rate sensors, studied in [5], and spacecraft not equipped with those sensors, studied in [6]. This because, in the first case, the control is based on the availability of measurements of both attitude and attitude rates and is realized through a PD-like control, while in the second case, the control can rely only on attitude feedback, and it is called attitude-only feedback control.

The first case, though presenting the mentioned difficulties, turns out to be the easier among the two, since the PD-like control law has only two parameters, and their robust optimal values can be successfully determined by means of a global search optimization procedure of the type of DIRECT algorithm [7], as proposed in [5]. In the second case, the control law contains four parameters, and their range of variation is wider than the first case. Thus, the determination of the robust optimal parameters becomes even harder, and a very complex derivative-free optimization algorithm based on a combination of both local and global search had to be developed. Moreover, numerical experience has shown that dealing with a discontinuous objective function makes the determination of the optimal solution very challenging even when using derivative-free algorithms. Thus, in Ref. [4, 6], the objective function is changed to the so-called integral time absolute error (ITAE) which is continuous with respect to both the parameters and the initial conditions. Such a change is acceptable since it has been shown that minimizing the ITAE is approximately equivalent to minimizing the settling time (see [8]).

This chapter describes in detail the abovementioned approaches to the determination of robust optimal values for the parameters in order to minimize the settling time or the ITAE obtained under the worst initial conditions. The exposition will highlight in particular the following main contributions: (i) the definition of a new systematic approach for the determination of the parameters for the PD-like control algorithm and for the attitude-only feedback, (ii) the formulation of a min-max optimization model to find a robust optimal solution to both problems, and (iii) the development of derivative-free optimization strategies to tackle the min-max problems.

The chapter is organized as follows: Section 2 defines the spacecraft model and the control algorithms for the two cases mentioned above; Section 3 describes the optimization model and the solution algorithm for the case of PD-like control; Section 4 provides some computational experience for this first case; Section 5 describes the optimization model and the solution algorithm for the case of attitude-only feedback control; Section 6 provides again computational experience for this second case.


2. Spacecraft model and control algorithms

The following coordinate frames are introduced to describe the attitude dynamics of an Earth-orbiting spacecraft and the Earth magnetic field:

Earth-centered inertial frame Fi. The origin of the frame is the center of the Earth, the xiaxis coincides with the vernal equinox direction, the ziaxis is the axis of rotation of the Earth and points northward, and the yiaxis completes the frame (see [1], Chapter 2.6.1).

Spacecraft body frame Fb. Its origin is in the spacecraft mass center. The axes are attached to the spacecraft and are selected so that the (inertial) pointing objective is having Fbaligned with Fi.

Since the objective is having Fbaligned to Fi, consider the relative kinematics and dynamics of the satellite with respect to the inertial frame. Let the attitude of Fbwith respect to Fibe represented by quaternion q=q1q2q3q4T=qvTq4T. The corresponding direction cosine matrix is equal to


where Iis the identity matrix (see [9], Section 5.4). Moreover, given aR3symbol a×denotes the skew symmetric matrix


so that the cross product a×bcan be expressed as the matrix multiplication a×b. The attitude kinematics equation is equal to q̇=Wqω(see [9], Section 5.5.3), where ωR3is the spacecraft angular velocity resolved in Fband


The attitude dynamics equation resolved in body frame is equal to Jω̇=ω×+T, where JR3×3is the inertia matrix of the spacecraft and Tis the control torque resolved in Fb(see [9]). Three magnetic coils aligned with the Fbaxes equip the spacecraft. Thus, the following magnetic attitude control torque is created


In the previous expression, mcoilsR3is the column matrix of the magnetic dipole moments for the three coils, and Bbis the Earth magnetic field at spacecraft resolved in body frame Fb(see [1], Section 7.4.1). Let Bibe the Earth magnetic field at spacecraft expressed in inertial frame Fi. Observe that Bichanges in time, at least because of the motion of the spacecraft along the orbit. Hence,


showing the explicit dependence of Bbon both qand t. The previous equations grouped together lead to the following system


where mcoilsis the control variable. Let us characterize the dependence on time of Bbqt, which is equivalent to characterize the time dependence of Bit. Assume a circular orbit with radius R. Through the adoption of the so-called dipole model of the Earth magnetic field (see [10], Appendix H), obtain:


In the previous equation, μmis the total dipole strength, ritis the position of the spacecraft resolved in Fi, and r̂itis the column matrix of the direction cosines of rit. The coordinates of vector m̂itare the direction cosines of the Earth magnetic dipole expressed in Fiwhich are written as follows:


where θmis the coelevation of the dipole, ωe=360.99deg/day is the average angular velocity of the Earth, and α0is the dipole right ascension at t=0. Set μm=7.7461015Wb m and θm=170.0as reported in [11].

Eq. (7) shows that to characterize the time dependence of Bit, one needs to determine an expression for ritwhich is the spacecraft position vector resolved in Fi. Define the orbital plane coordinate system ap, bpwhose origin is at the Earth center and with apaxis coinciding with the line of nodes. Then, the center of mass of the satellite is positioned in


where nis the orbital rate and ψis spacecraft argument at time t=0. It is possible to determine the mass center coordinates in the inertial frame Fifrom (9) by means of a rotation matrix which depends on the inclination inclof the orbit and on the right ascension of the ascending node (RAAN) Ω(see [1], Section 2.6.2). By inserting the expressions of the latter coordinates into (7), we obtain an explicit expression for Bit.

Since Cq=Ifor q=qvTq4T=±q¯where q¯=0001T(see (1)), then the goal is designing control strategies for mcoilsso that qv0and ω0. Reference [3] proposes the following stabilizing proportional derivative (PD)-like control law, which is a modification of those in [12, 13]:


In the previous equation, mcoilsis the saturation limit on each magnetic dipole moment, “sat” denotes the standard saturation function, and κpand κdare parameters. As shown in [3], if the orbit’s inclination inclis not too low, there are large ranges for the values of κp>0and κd>0which lead to local exponentially stability of equilibrium qω=q¯0for the closed-loop systems (6) and (10).

The following attitude-only feedback, obtained as modification of one described in [12], is also proposed in [3]


In (11), δR4is an internal state of the controller, κ1κ2αβare all parameters, and Wqwas introduced in (3). Note that the previous equation describes an attitude feedback, since it requires only the measure of attitude qand not of attitude rate ω. It has the advantage of not requiring the installation of rate gyros obtaining savings in cost, volume, and weight. As shown in [3], if the orbit’s inclination inclis not too low, there are large ranges for the values of κ1>0, κ2>0α>0, β>0which lead to local exponential stability of equilibrium qωδ=q¯01ϵλq¯achieved for the closed-loop system (6) and (11).

For both control laws, there are no precise indications for choosing the parameters. In practice, they are usually computed recurring to a trial-and-error approach. This causes the two main limitations: (1) the computation is quite time-consuming, and (2) it is not systematic. This means that when a satisfactory set of parameter values is finally obtained, it is unknown whether continuing the search could allow to discover new parameter values producing a better performance of the closed-loop system. And, if the search is continued, it is unknown when it should be stopped and what could be the possible improvements obtained with this additional search. In any case, it can easily happen to neglect values providing an overall better performance, unless all possible values are tried. However, such an exhaustive search is generally impracticable, because the search space is way too vast to be completely explored. Therefore, we describe here a different approach to find the values of the feedback parameters. Since the desired attitude is obtained when qv=0, we define the settling time tsifor each component qi, with i1,2,3, as:


In other words, this is the time needed for qito permanently remain under ν. Value 0<ν<1can be set depending on the desired final value of qi. Finally, we can define the settling time tsfor the whole quaternion qas the time required by the slowest component of qv, hence


Now, having set the spacecraft initial conditions to specific values, one can determine the values of the parameters that minimize the settling time tsfor each control law (10) and (11).


3. Determination of optimal parameters for the case of PD-like control

In the case of PD-like control (10), the minimization of the settling time can be formulated as follows


In order to practically solve the problem, the feasible set should be reasonably bounded and not infinite. Hence, we define two upper bounds κp̂and κd̂for the gains κpand κd. These values can usually be determined for the specific problem; we do this for our case study in Section 3.3. Thus, the problem becomes as follows, where the dependence of tson the initial conditions q0q0, ω0ω0, and ψ(see (9)) is explicitly indicated,


and the feasible set is K=κpκd:0κpκp̂0<κdκd̂. Even though tsobviously depends on κpand κdand on initial conditions, it is not possible in practice to express this relation in analytical form. Moreover, tsis discontinuous with respect to κp, κdand initial conditions See for example in ([14], p. 233) the prove that the settling time of the step response is not continuous with respect to the system parameters.

Now, problem (15) has the following features: (i) an analytic expression of the objective function cannot be given, and (ii) the objective function is not continuous with respect to the decision variables. Therefore, we cannot use standard optimization techniques. Instead, we need to use derivative-free optimization.

3.1 From the simple min to the min-max problem

Given the initial conditions q0ω0ψ, problem (15) can be solved to optimality by using the global optimization derivative-free technique described in Section 3.2. However, when the initial conditions change, then that solution may be no longer optimal. An example of this is given in Section 4.1. Therefore, a robust approach consists in searching for the optimal solution to (15) under the worstinitial condition. Such a worst-case approach is commonly used in similar cases. In other words, we accept to pay something, in terms of objective function values, in the easy cases, but we obtain in return advantages in the more difficult situations. However, the worst initial condition is not a priori known, since it depends on the chosen values of κpand κd. Therefore, the problem cannot be decomposed and should be solved as a whole.

The set of values describing the initial conditions is given by:


Set Sincludes any possible initial attitude and any possible initial argument ψfor the spacecraft; only the magnitude of the initial angular rate is limited. The worst-case minimization of tscan be modeled as follows:


To use the derivative-free algorithm explained in Section 3.2, the feasible set of each optimization problem must be now converted into a hyperrectangle. This can be obtained by expressing q0v1in spherical coordinates ρϕθ:


The dependence of tson q0can be expressed as dependence on the variables ρϕθ, as follows. We introduce the hyperrectangle:


Now, the min-max problem (17) can be equivalently reformulated as follows:


3.2 The whole derivative-free optimization approach

This section explains the min-max procedure to solve problem (20). Its building blocks are given by the following derivative-free algorithms:

3.2.1 The global strategy

Lipschitzian methods constitute a main approach in non-differentiable optimization. However, they are limited by their requirement of knowing the value of the Lipschitz constant. On the other hand, the search may be conducted without knowing the Lipschitz constant if we use DIRECT-type algorithms [15]. These techniques are based on the partitioning of the feasible region into hyperrectangles and on their examination in a specific order. The feasible region starts as a single hyperrectangle that is internally normalized to a unit hyperrectangle. At each iteration k, the algorithm partitions the hyperrectangles identified in the previous (k−1)th iteration, obtaining a collection of smaller hyperrectangles Hk=H1Hk, and evaluates the objective function in their central points. The computation of the objective value allows to identify potentially optimal hyperrectangles within Hk. These hyperrectangles are further partitioned and investigated in the next (k+1)th iteration of the algorithm, while the rest of them is simply discarded. The algorithm stops when the hyperrectangle size becomes very small or when the maximum number of iterations is reached. This algorithm guarantees to converge to the global optimum of the function if sampling is dense enough. However, a large number of function evaluations might be needed to obtain dense sampling, and consequently large computational times could be required.

3.2.2 The local strategy

When the function under optimization is sufficiently regular, it can be optimized by using gradient-based methods, which require the computation of the first-order derivatives of the objective function. However, the basic strategy underlying these methods could be replicated even if no information on the derivatives is available. This is the basic idea of the SDBOX algorithm [16], which can be described as follows. After the selection of a starting point, it cyclically determines a feasible and gooddescent direction for the objective function, and then it performs a sufficientlylarge step along such direction. To find the good feasible descent direction, at each iteration k, the algorithm computes the local behavior of the objective function at the variation of the i-th variable. If a move of length αalong digives a feasible point where the function is reduced enough, the algorithm applies a linesearch technique dito find the stepsize αk. Otherwise, the algorithm tries the same operations in the opposite direction di. If both diand diare not able to obtain a sufficient decrease, the stepsize αis decreased. The linesearch technique does not require information on the slope of the objective function. Note that the convergence of the algorithm was proved in [16] for minimization of a continuously differentiable function; however, the same technique is often used to optimize different types of functions. This approach is considerably faster than the global strategy; however it is well known that a poor choice of the initial point leads to poor solutions.

3.2.3 The whole procedure

To simplify the description of the robust optimization approach, we now rename the set of variables κpκdas xbelonging to a feasible set Fx=lbixiubii=1nRn(in our case Fx=Kand n=2) and the set of initial conditions ρϕθω0ψas ybelonging to a feasible set Fy=lbjyjubjj=1mRm(in our case Fy=Hand m=7, since ω0has three components). Moreover, we now use gxyto denote the objective function (in our case ts). Its analytical expression is not known; however its value can be computed via a software simulation.

The upper-level minimization problem minxFxgxymust be solved by means of a global strategy applied one time: no starting points can be considered here; hence, there is no ground for a local strategy. Thus, the procedure will contain an external loop implementing the global strategy. This loop essentially computes the value of gobtained for a number maxeval_extof points x. Let x¯be one of them. The evaluation of x¯requires however to solve one lower-lever maximization problem maxyFygx¯y. Hence, the lower-lever problem should be solved many times, and its corresponding computation time must be reduced to a few seconds. We can solve it in three ways: (1) with the global strategy but using a necessarily small maximum number of function evaluations maxeval_int; (2) with the local strategy, allowing the same number of function evaluations from a single start; and (3) with the local strategy, using a multistart from different starting points, still allowing maxeval_inttotal function evaluations. This last option allows to use information that happens to be available in the choice of the starting points. In our case, from the physics of the problem, we suppose that good solutions are in the vicinity of extremal values of angular velocity, so we take them as starting points. We report the whole procedure:

Procedure 1: Solve min-max

  1. Input:A gxycomputable by means of a software simulation for xilbiubiwith i=1,,nand yjlbjubjwith j=1,,m.

  2. Output:A robust optimum value x=argminxFxmaxyFygxy.

  3. External loop:

  4. Solve the upper-level problem minxFxfx, with fx=maxyFygxyby using DIRECT for maxeval_extevaluations of f, and return x. The k-th evaluation works with point xk

  5.   Internal loop:

  Solve the lower-level problem maxyFygxky

  by using


4. Computational results for the case of PD-like control

We apply the above approach to solve the numerical example studied in [3]. The inertia matrix of the spacecraft is equal to J=diag271725kgm2; the saturation level for each magnetic dipole moment is given by mcoils=10Am2. The orbit has an inclination of incl=87and an altitude of 450 km. The orbital period in these conditions is about 5609 s. The value Ωof RAAN is 0. Given these values, we can compute an upper bound κp̂for κp. Simulation shows that BitBmin=2.4×105T for the considered orbit; since Bbqt=Bit, it follows that BbqtBmin. Assume that the resolution of the attitude sensor is qs=0.04in terms of quaternion component. In this case, an upper bound for κpcan be computed by enforcing that each component of the term Bbqt×κpqv(see (10)) does not exceed the saturation limit mcoils=10Am2when Bbqtand qvare orthogonal, Bbqt=Bminand qi=qrfor one index i1,2,3, and qi=0for the other ones. This is achieved if κpis smaller than mcoils/Bminqr=1.0417×108; thus, we set κp̂=108. Assuming that the attitude rate sensor has a resolution of ωr=4×104rad/s, by a parallel argument it follows that an upper bound for κdis given by mcoils/Bminωr=1.0417×109; thus, we set κd̂=109.

4.1 Fixed initial conditions

Consider an initial state characterized by attitude equal to the target attitude q0=q¯(which corresponds to having ρ=0and any value for ϕand θ) and by the following initial angular rate:


This example corresponds to a spacecraft possessing the desired attitude and no angular momentum, which received an impact from a small object leading to an instantaneous change in the spacecraft angular rate. We assign to ψa random value over the interval [0,2π[by setting ψ=0.332rad.

Gains κpand κdhave been found by trial-and-error in [3] as κp=2×105and κd=3×108. The settling time corresponding to these values of the gains is ts = 4.002 orbital periods, while the vast majority of possible gain values would produce settling times larger than 10 orbital periods. On the contrary, if we use the optimization algorithm described in Section 3.2, we obtain the following fixed initial condition optimal gains after only 32,821 iterations (i.e., 251 s of computation):


The corresponding optimal settling time is ts=1.775 orbital periods. Hence, the optimization algorithm can provide a significant improvement in convergence speed.

However, the above gains might become no longer optimal if we assume different initial conditions. Consider, for example, the case of initial conditions with the same attitude q0=q¯(i.e., ρ=0and any ϕand θ), the same argument ψ=0.332rad, but a different initial angular rate


The above gains yield a settling time ts = 3.562, and they are no longer optimal, as shown in Section 4.2. Indeed, the optimal values of the gains would need to be determined for every possible initial condition, which is clearly impossible in practice. As seen in Section 3.1, we can instead search for the gain values providing the best performance under the worst initial condition. It has been found that the worst initial conditions corresponding to gains (22) are


Under conditions (24), the above values of κpand κdyield a settling time ts=4.338.

4.2 Variable initial conditions

We now search for κpand κdusing the robust optimization approach. We allow initial conditions to vary as follows: 0ρ1, 0ϕ2π, 0θπ, ω010.1, ω020.1, ω030.1, and 0ψ2π. We allow a maximum of 2,200 function evaluations in each internal loop, in order to practically solve the problem. Table 1 reports the results of Procedure Solve min-max with maxeval_ext = 10,000 and maxeval_int = 2200.

ext eval10,64111,12710,011
Time (s)172,242.657,255.189,640.8

Table 1.

Results of Solve min-max with maxiter_ext = 10,000 and maxiter_int = 2200.

When using SDBOX multistart, the starting points are ρ=0.5, ϕ=π, θ=π/2, ψ=π, and all the eight combinations of extreme angular velocities,,,

These results show that (i) SDBOX multistart is able to reach the highest value of the objective tseven with the small number of computations allowed in the internal loop and (ii) the different results in the solution of the lower-level problems cause a different evolution of the upper-level search. To understand which among the three points above is the best choice, we continue the analysis in Table 2. We fix parameters κpand κd, and we solve the maximization problem with increased accuracy by allowing more function evaluations. We do this by means of global strategy, local multistart strategy (which proved to dominate the single start one), and an exhaustive grid search, which is much slower but used here as reference.

GainsDIRECTSDBOX multistartGrid search
eval = 40,000eval = 40,000 (5000 × 8)eval = 78,125 (57)

Table 2.

Accurate evaluation of the previously obtained gains.

This analysis shows that the first point has the smallest maximum settling time. In conclusion, the robust optimal solution is:


To evaluate the performance of the above solution (26), we try it with the three different initial conditions reported above. If the initial state is (21), the settling time becomes ts=1.945(instead of 1.775). Hence, there is a small worsening. Clearly, no improvement was possible because gains (22) are optimal for state (21). If the initial state is given by (23), we obtain ts=2.958(instead of 3.562). Hence, there is an improvement. This could be expected, though not guaranteed, since gains (22) have a higher worst-case result than gains (26).

Finally, if the initial state is given by (24), then we obtain ts=3.370(instead of 4.338). Hence, there is a substantial improvement. Indeed, this was certain, since gains (26) have a worst-case result considerably better than gains (22). In conclusion, the gains computed with the proposed robust optimization approach offer improvements in the difficult situations. Moreover, even if they may be suboptimal in the easier cases, the balance appears beneficial.


5. Determination of optimal parameters for the case of attitude-only feedback control

As observed, the objective function tsis not continuous with respect to the parameters, and this introduces numerical difficulties in solving the optimization problem. Thus, one can consider an alternative objective function named integral time absolute error (ITAE) [8], indicated by Γ:


where represents the Euclidean norm and Tfis a time which is selected large enough. The ITAE carries the benefit of being continuous with respect to the design parameters. Continuous differentiability cannot be guaranteed. However, since qωδ=±q¯01βq¯are equilibria of the closed-loop system (6) and (11), it is very unlikely that qvbecomes 000Tat some finite time. Note that this would make Γnot continuously differentiable with respect to the design parameters because qvis not continuously differentiable at qv=000T.

Minimizing the ITAE is known to lead to nearly optimal solutions with regard to the settling time. In fact, in very simple situations, it is shown analytically that minimizing the ITAE leads to solutions that minimize the settling time. In more complex scenarios, it was numerically shown that minimizing the ITAE gives solutions that are very close to the optimal ones in terms of settling time (see [8]).

By introducing physically reasonable upper bounds κ1̂,κ2̂,α̂,β̂for the design parameters, we obtain the feasible set K={κ1κ2αβ:0κ1κ1̂,0κ2κ2̂,0αα̂,0ββ̂}. Now, our optimization problem is:


Given specific initial conditions of the spacecraft, problem (28) can be solved by a suitable use of derivative-free techniques. However, if the initial conditions of the spacecraft change, that solution might be no longer optimal. Since several different initial conditions for the spacecraft are possible in practical situations, a robust solution is an optimal solution to problem (28) under the worstspacecraft initial conditions. Such a worst-case optimization is widely employed in such scenarios, because by adopting the latter approach, we can give an efficient bound on the objective value in spite of the uncertainty on the spacecraft initial conditions. However, the worst initial conditions for the spacecraft cannot be determined a priori, since they depend on the selected values of κ1,κ2,α,β. The initial conditions of the spacecraft are given by q0=q0, ω0=ω0, 0ψ<2π, and 0α0<2π. We chose the set of their possible values as:


Note that Sincludes all possible initial attitudes, all possible right ascensions of the Earth magnetic dipole at time t=0, and all possible initial arguments ψfor the spacecraft. It only constrains the amplitude of the initial angular velocity. The minimization of Γunder the worst spacecraft initial conditions is equivalent to formulating the following min-max problem:


To apply the optimization techniques, we convert the feasible set of each optimization problem into a hyperrectangle by expressing the set q0v1in spherical coordinates ρϕθ:


The dependence of Γon q0can now be expressed as dependence on the variables ρϕθ. Consequently, after having introduced the hyperrectangle


the min-max problem (30) can be equivalently reformulated as follows:


5.1 Combining global and local search

To simplify the description of the proposed approach, we now rename the set of design parameters κ1κ2αβas xbelonging to a feasible set Fx=lbxixiubxii=1nRn(in our case Fx=Kand n=4) and the set of initial conditions ρϕθω0ψα0as ybelonging to a feasible set Fy=lbyjyjubyjj=1mRm(in our case Fy=Hand m=8, since ω0possesses three components). Indicate by fxythe function giving the objective value (in our case Γ). When yis fixed (the initial conditions are assigned), we simply write y¯in it; when xis fixed (the design parameters are assigned), we write x¯in it. Problem (28) can be written as


and may be tackled by a global derivative-free optimization algorithm of the type of DIRECT [7], already described in Section 3. 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. Due to the so-called everywhere denseproperty, such an algorithm reaches a global optimum if the sampling is dense enough. However, a dense search may require that the function is evaluated many times.

In our case, all lbxi=0, since design parameters have to be greater than or equal to zero. Values ubxican typically be set to very large values based on some physical considerations. For instance, κ1̂can easily become equal to 109, since its maximum feasible value can be determined knowing the saturation level of the coil moments, the minimum amplitude of the geomagnetic field, and the attitude sensor resolution. However, when considering these values, a sufficiently dense exploration of the feasible set Fxrequires a number of function evaluations such that the corresponding run time is impracticable. On the other hand, a sampling that uses a practically sustainable number of function evaluations does not produce solutions substantially better than random solutions. The presence of a large number of local minima makes the optimization task particularly difficult. In such conditions, the evaluation of the generic Hhusing only one point may be very inaccurate at the first iterations of the algorithm, because the dimension of the hyperrectangles is too large. By proceeding with the iterations, the hyperrectangles becomes smaller, but their number, and consequently the run time, increases exceedingly.

On the other hand, if a fast but effective probing technique to early identify the “promising” Hhwould be available, then one could explore densely only such promising regions in reasonable times. Hence, a probing technique based on the use of the local derivative-free optimization algorithm SDBOX was proposed in [6]. This algorithm was originally presented in [16] as a globally convergent algorithm for the minimization of a continuously differentiable function, but in practice it can be employed to optimize different types of functions as a good trade-off between efficiency and convergence properties. Most interestingly for our case, given an initial guess κ1κ2αβ=x, this algorithm should be able to find good solutions in short times in the neighborhood of x. We describe below a solution approach combining these global and local strategies to solve problem (28).

Procedure 2: Solve min combining global and local search

  1. Input:A vector y¯Fyand a function fxy¯computable by means of a software simulation for any xFx. Values for the parameters p, maxeval, maxsubsets, maxiter, maxpost.

  2. Output:A solution xapproximating one vector in argminxFxfxy¯.

    1. Normalize and grid partition the whole feasible set Fxinto a collection of hyperrectangles H1=H1Hpsimilarly to the initial phase of DIRECT.

    2. For each HhH1, compute the value fhof the solution obtained by maxevaliterations of SDBOX in Hhstarting from its central point. This is an upper bound on the value of the best solution in Hhand constitutes our “evaluation” of Hh.

    3. Take a number maxsubsetsof hyperrectangles corresponding to the smallest of the above fhvalues.

    4. Take the region given by the union of those subsets, and “convexify” it by including also the additional subsets required to convert it into an hyperrectangle Fx.

    5. Switch to DIRECT algorithm to continue the search in Fxallowing maxiterfunction evaluations. This search can now be dense using reasonable time, and it gives a solution x.

    6. Try to improve xby using maxpostiterations of a local search method, finally obtaining a solution xto problem (28).

Depending on the practical case, the number of hyperrectangles pin the partitioning phase (step 1) and the number of iterations maxevalin the evaluation phase (step 2) must be selected in order to allow a fixed computational time to the evaluation phase. Indeed, we search for a compromise between speed and effectiveness at this stage. Those parameters must be tuned by also considering that the accuracy of the evaluation of each Hhdepends not only on maxevalbut also on the span of each Hh, which in turn depends on p. However, there is no need that all the hyperrectangles in the initial partition H1have the same size. Their size could be growing with the absolute values of the coordinates so that pis kept smaller. The number maxsubsetsin the selection phase (step 3) must be chosen so that Fxremains much smaller than Fx; otherwise the benefits of the proposed procedure is reduced. The number of function evaluations maxiterin the standard DIRECT phase (step 5) is chosen so that the search in Fxis dense enough. This is now possible because of the size reduction in Fx. The post-optimization phase (step 6) can be executed with SDBOX or with CS-DFN [17], which is another linesearch-based method which uses a dense set of search directions and not only the coordinates ones. CS-DFN does not require fto be continuously differentiable; however, it could be more computationally expensive on smooth problems. For example, in our case, Γis at least continuous, and even if there are no theoretical arguments to ensure its continuous differentiability, in practice this property may often occur. In conclusion, we can choose between SDBOX and CS-DFN depending on the presence or absence of the continuous differentiability of f, and if we set maxpostto be sufficiently large, we obtain a solution xwhich satisfies necessary conditions for a local optimum. Moreover, xshould approximate one of the global optima, because, given the dense search in F, it should provide one of the global minima of F, and the regions of F\F, which were less “promising,” should not contain better solutions. Note that the accuracy and the properties of the evaluation phase can be modified, depending on the computational request of the practical case, either by modifying the number of iterations or even by employing an alternative evaluation algorithm, still keeping the same algorithmic framework.

On the other hand, problem (33) is made of an upper-level minimization problem and a lower-level maximization one. By changing the names of the set of variables as explained at the beginning of the section, the problem becomes


with function gsuch that its value at the generic point x¯is given by the solution of the lower-level problem:


The upper-level problem is solved through an external loop applying the combination of local and global search described as Procedure 2. This loop computes, using a parameter maxeval_ext, to define the overall maximum number of function evaluations and the value of gcorresponding to different points of Fx. Let x¯be one of them; then the evaluation of gx¯requires the solution of one lower-level maximization problem maxyFyfx¯y. Thus, the lower-level problem must be solved up to maxeval_exttimes. Consequently, solving it in a few seconds is crucial.

For this problem, the global strategy either would perform a very poor search or would need excessive time. Therefore, the local strategy appears to be the only feasible choice for the lower-level problem. However, taking as initial guess the center of the feasible set Fydoes not lead to good solutions of the maximization problem within the limited available time.

In this case, based on the problem physics, it can be assumed that good solutions of the maximization problem are in the neighborhood of extreme values of angular velocity. Therefore, we solve the lower-level problem in a nested loop by using a local search with multistart. We take as initial guesses the eight combinations of extreme values for the three components of the angular velocity ±ω01̂, ±ω02̂, and ±ω03̂. This can be done with SDBOX or CS-DFN, depending on the presence or absence of the continuous differentiability of f. In our case, Γis at least continuous also with regard to the initial conditions of the spacecraft. Again, there are no theoretical arguments to ensure its continuous differentiability; however in practice this may often happen. For each solution of the lower-level problem, we allow a necessarily small maximum number of function evaluations maxeval_int. In conclusion, this produces a solution xRfor problem (33). We report here the whole procedure:

Procedure 3: Solve min-max combining global and local search

  1. Input:A fxycomputable by means of a software simulation for any xFxand yFy. Values for the parameters maxeval_ext, maxeval_int, s.

  2. Output:A robust solution xRapproximating one vector in argminxFxmaxyFyfxy.

  3.    External loop:

Solve the upper-level problem minxFxgx, with gx¯=maxyFyfx¯yx¯Fx,

by using Procedure 2 with maxeval_exttotal evaluations of g, and return xR.

Given x¯, the evaluation of gx¯is performed by the internal loop.

  1.       Internal loop:

      Take x¯, and solve the lower-level problem maxyFyfx¯y

      by using multistart local search with sstarting points,

      performing maxeval_intsevaluations of ffor each of them.


6. Computational results for the case of attitude-only feedback control

We apply the above method to solve the numerical example presented in [3]. The spacecraft inertia matrix is J=diag271725kgm2, and the saturation level for each dipole moment is mcoils=10Am2. The orbit inclination is incl=87, and the orbit altitude is 450 km; the right ascension of the ascending node Ωis equal to 0. Upper bounds κ1̂,κ2̂,α̂,β̂are selected as 109109104103.

At the beginning we consider the easier situation of known spacecraft initial conditions, and we present, in Section 6.1, the results of Procedure 2 in solving this problem. We report also a comparison with the classical DIRECT method. Subsequently, in Section 6.2, the more realistic case of a spacecraft having variable initial conditions is considered.

6.1 Fixed initial conditions

We consider here the case of the above described spacecraft with known fixed initial conditions; thus we deal with problem (28) using the following values:


The best solution obtained in [3] by trial and error search achieves a value of ITAE=3.7×107, while the vast majority of the solutions have the ITAE limit value of about 1.2×109. This upper limit is related to the value of Tfin the definition of ITAE (27), which is chosen equal to 56,009 s corresponding to 10 orbital periods. Practically this means that when we reach the limit value of ITAE, the corresponding settling time would be roughly larger than 10 orbital periods. Then, that solution is not an attractive one, and we are not interested in determining it with further precision. Even if this choice causes a flattening in the values of ITAE, the use of such a finite Tfis necessary to run the simulations that compute the ITAE in practice.

Table 3 reports two solution attempts carried on with the standard DIRECT algorithm on the whole feasible set Kwith, respectively, 50,000 and 100,000 iterations, followed by 1000 iterations of local search refinement using CS-DFN. In spite of the significant computational effort (the running times of these experiments, respectively, correspond to about 3 days and 1 week), the obtained solutions have values of ITAE greater than 1.1×109that is not very different from the ITAE value of a random solution. Evidently, the search was not dense enough to explore the feasible set as it would be needed. Other similar attempts with standard DIRECT do not achieve better results. The same table also reports Procedure 2 of Section 5.1. As observable, this procedure is able to reach a much better solution, with a value of ITAE of about 8×106. Note that this solution is also considerably better than the best solution obtained by trial and error in days of work.

AlgorithmSolutionObj. valueTime
DIRECT 50,000 + CS-DFN 1000κ1 = 913405022.139
κ2 = 195426826.870
α = 9794.170752422
β = 0.000000000000
1,142,470,478.101262,000 + 5460 s
DIRECT 100,000 + CS-DFN 1000κ1 = 500798437.500
κ2 = 159788790.177
α = 9996.679486501
β = 0.000000000003
1,129,234,873.703565,200 + 5040 s
Procedure 2: Combining global and local strategiesκ1 = 246494.579020
κ2 = 233333315.349
α = 92.5925925927
β = 0.000129629629
8,021,573.4077201,500 + 15 + 4 s

Table 3.

Comparison of Procedure 2 and standard DIRECT with 50,000 or 100,000 iterations.

The promising region of the feasible set Khas been identified with the previously described probing technique. The values of maxelavand pare selected so that this step is executed in reasonable time, according to the following considerations. One single function evaluation takes a time which is very variable and goes from fractions of seconds to several tenths of seconds; but a very rough average function evaluation time can be assessed as equal to 1 s. By selecting maxeval = 10, the evaluation of each hyperrectangle would require roughly 10 s. Thus, to finish within 2 or 3 days of computation, we select p=40,960, which is obtained by making 64 partitions on the domain of κ1, 64 partitions on the domain of κ2, and 10 partitions on the domain of α. The intervals corresponding to these partitions do not have the same size; they increase with the absolute value of the coordinates. The following hyperrectangle is obtained as convexification of the collection of the most promising regions, as in Step 4 of Procedure 2:


The determination of Kactually needs 409,600 iterations of SDBOX and 201,500 s of computations (about 56 hours). Note that the time necessary for each function evaluation is generally much faster in K, where it can be less than 0.01 s than in the rest of K. Now, by applying standard DIRECT strategy over the feasible set K, as in Step 5 of Procedure 2, after 3007 iterations and only 15 s, we obtain the solution:


whose value is ITAE = 9.072×106. Note that the volume of the set K*is considerably smaller than that of K: it is only 1/15151515.15of the volume of K. To better appreciate the difference, consider that an exploration of Kwith the same degree of density used on K*would require 15×15,151,515.15=227,272,727.25 s that roughly corresponds to more than 7 years, if the simulation times on Kwere the same as on K*. Since they are often much slower, the time needed would be even more.

Solution (39) can be further refined by using the local search strategy, as in Step 6 of Procedure 2. By performing 1000 iterations of the local search CS-DFN, which can move along a dense set of directions [17], the solution (39) is improved in 4 s to ITAE = 8.021×106with the solution:


We also tested in this post-optimization phase the local search SDBOX, which moves only along the coordinate directions [16]. With 5000 additional iterations of SDBOX that requires 11 s, solution (39) is improved to the following new solution, which has ITAE = 9.062×106


Solution (40) corresponds to xof Procedure 2. Note that this is actually an approximation of a theoretical optimal solution. Since there is probably no optimal solutions available for comparison, we evaluate it by using the following considerations.

The settling time (see (12) and (13)) corresponding to (40) is computed as ts=6280 s. Then, we search for a lower bound ts¯on the minimum settling time of the considered case, through physical considerations. The initial conditions (37) correspond to having the spacecraft with the desired attitude but with a nonzero initial angular velocity ω0= Consider now simple rotations about each single body axis, and for each simple rotation, compute lower bounds tsx¯, tsy¯, tsz¯of the times necessary to move to the desired attitude with final zero angular rate. Then, a rough lower bound for the settling time is given by ts¯=maxtsx¯tsy¯tsz¯. Value tsx¯can be computed using the equation which describes rotation about the xbody axis, which is given by


where ϕis the roll angle. The amplitude of torque Txis limited by an upper limit T, which can be found using (4). Indeed, numerical simulations show that BiB=5105T. Since Bb=Bi(see (5)), then BbB. Moreover, each component of mcoilsis bounded by mcoils=10Am2and then mcoils3mcoils. Thus, T=3Bmcoils=53104N m. Next, the minimum time to bring the state of system (42) subject to the constraint TxT, from the initial state ϕ=0ϕ̇=ωx0to the final state ϕ=0ϕ̇=0, is given by (see [18], Section 7.2)


Similar considerations hold for rotations about yand zaxes leading to


Then, ts¯=maxtsx¯tsy¯tsz¯=2091s. Therefore, solution (40) takes just about 70 min more than the minimum time necessary to rotate the spacecraft about a single body axis at the maximum speed allowed by the available magnetorquers so that it reaches the desired rest position. Thus, solution (40) does not appear to be too far from an optimal solution. Note also that the above computed bound is very conservative, in the sense that the minimum spacecraft evolution time surely cannot require less, though it could very easily require more.

6.2 Variable initial conditions

Now, problem (33) is solved by employing the above described Procedure 3. We select 20 partitions on the domain of κ1, 20 on that of κ2, and 10 on that of α. We therefore perform 40,000 evaluations of gxfor the identification of the new Kand then another 3000 + 1000 evaluations of gxto solve the problem on this K, for a total of maxeval_ext=44,000.

For the internal loop, as described in Procedure 3, we employ as initial guesses the eight combinations of extreme values for the three components of the angular velocity ±ω01̂, ±ω02̂, and ±ω03̂. We allow 64 iterations per point, for a total of 512 iterations per single lower-level problem, which requires slightly less than 3 s in average within K. Run times are greater in the rest of K; however, we try to keep them under control by allowing to exit the internal loop when the value of fx¯yis large enough to reach the limit value for ITAE of about 1.2×109. Then, the whole nested loop procedure provides the following solution in about 342,000 s (about 95 h)


The ITAE value of this solution is 2.176×107for the initial conditions (37) instead of ITAE =8.021×106of solution (40). However, this solution is a robust solution: by varying the initial conditions in H, the worst value that can be obtained is ITAE =1.357×108that is still considerably better than average solutions, whose vast majority has the limit value for ITAE of 1.2×109. As a comparison, the worst value achievable by varying the initial conditions in Hfor solution (40) is ITAE =1.103×109that is very close to the limit value for ITAE. Indeed, the above limit value for ITAE is very easily obtainable for almost any generic tuple κ¯1κ¯2α¯β¯by simply searching for difficult initial conditions. Note also that this value of ITAE is obtained in some attempts in solving problem (33) by using standard DIRECT algorithm on the whole feasible set Kallowing no more than 1 day of computation. More time-consuming attempts in solving problem (33) using standard DIRECT algorithm cannot be practically accomplished. This holds for the following motivations. In the case of fixed initial conditions (Section 6.1), 1 week of computation was not enough to explore the search space. In this case, there is an internal loop requiring 512 function evaluations instead of one single function evaluation. Hence, any serious attempt would need to allocate months of computation and would obtain results probably similar to those obtained in Section 6.1.


7. Conclusions

The attitude of a spacecraft can be controlled using only magnetorquers by means of a PD-like control or an attitude feedback. However, for both control laws, design parameters must be assigned. These parameters may be conveniently selected so that they minimize the spacecraft settling time, or an indirect measure of it, either for fixed initial conditions of the spacecraft or under the worst initial conditions. This latter choice gives an upper bound on the minimum value of the objective obtainable by varying the initial conditions. This chapter has described solution approaches based on an innovative use of global and local derivative-free optimization techniques to practically solve these computationally demanding problems. This approach is able to provide robust solutions to the considered application in reasonable times.


  1. 1. Sidi MJ. Spacecraft Dynamics and control. New York: Cambridge University Press; 1997
  2. 2. Silani E, Lovera M. Magnetic spacecraft attitude control: A survey and some new results. Control Engineering Practice. 2005;13(3):357-371
  3. 3. Celani F. Robust three-axis attitude stabilization for inertial pointing spacecraft using magnetorquers. Acta Astronautica. 2015;107:87-96
  4. 4. Bruni R, Celani F. Determining optimal parameters in magnetic spacecraft stabilization via attitude feedback. AIP Proceedings. 2016;1776:090032
  5. 5. Bruni R, Celani F. A robust optimization approach for magnetic spacecraft attitude stabilization. Journal of Optimization Theory and Applications. 2017;173:994-1012
  6. 6. Bruni R, Celani F. Combining global and local strategies to optimize parameters in magnetic spacecraft control via attitude feedback. Journal of Optimization Theory and Applications. 2019;181:997-1014
  7. 7. Jones DR, Perttunen CD, Stuckman BE. Lipschitzian optimization without the Lipschitz constant. Journal of Optimization Theory and Application. 1993;79(1):157-181
  8. 8. Graham D, Lawthrop RC. The synthesis of optimum response: Criteria and standard forms, part 2. Transactions of the AIEE. 1953;72:273-288
  9. 9. Wie B. Space Vehicle Dynamics and Control. Reston, VA: AIAA; 2008
  10. 10. Wertz JR, editor. Spacecraft Attitude Determination and Control. Norwell, MA: Kluwer Academic; 1978
  11. 11. Rodriguez-Vazquez AL, Martin-Prats MA, Bernelli-Zazzera F. Full magnetic satellite attitude control using ASRE method. In: 1st IAA Conference on Dynamics and Control of Space Systems. 2012
  12. 12. Lovera M, Astolfi A. Spacecraft attitude control using magnetic actuators. Automatica. 2004;40(8):1405-1414
  13. 13. Lovera M, Astolfi A. Global magnetic attitude control of inertially pointing spacecraft. Journal of Guidance, Control, and Dynamics. 2005;28(5):1065-1067
  14. 14. Ogata K. Modern Control Engineering. 4th ed. Upper Saddle River, NJ: Prentice Hall; 2002
  15. 15. Jones DR. DIRECT global optimization. In: Floudas CA, Pardalos PM, editors. Encyclopedia of Optimization. Berlin: Springer; 2009. pp. 725-735
  16. 16. Lucidi S, Sciandrone M. A derivative-free algorithm for bound constrained optimization. Computational Optimization and Applications. 2002;21:119-142
  17. 17. Fasano G, Liuzzi G, Lucidi S, Rinaldi F. A linesearch-based derivative-free approach for nonsmooth constrained optimization. SIAM Journal on Optimization. 2014;24(3):959-992
  18. 18. Athans M, Falb PL. Optimal Control: An Introduction to the Theory and Its Applications. New York: Dover; 2007

Written By

Renato Bruni and Fabio Celani

Submitted: February 19th, 2019 Reviewed: August 16th, 2019 Published: January 15th, 2020