Results of Solve min-max with maxiter_ext = 10,000 and maxiter_int = 2200.
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
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 , 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 . In particular, Celani  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 optimal values.
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:
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.
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 , and spacecraft not equipped with those sensors, studied in . 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 , as proposed in . 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 ).
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 . The origin of the frame is the center of the Earth, the xi axis coincides with the vernal equinox direction, the zi axis is the axis of rotation of the Earth and points northward, and the yi axis completes the frame (see , Chapter 2.6.1).
Spacecraft body frame . 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 aligned with .
Since the objective is having aligned to , consider the relative kinematics and dynamics of the satellite with respect to the inertial frame. Let the attitude of with respect to be represented by quaternion . The corresponding direction cosine matrix is equal to
where I is the identity matrix (see , Section 5.4). Moreover, given symbol denotes the skew symmetric matrix
so that the cross product can be expressed as the matrix multiplication . The attitude kinematics equation is equal to (see , Section 5.5.3), where is the spacecraft angular velocity resolved in and
The attitude dynamics equation resolved in body frame is equal to , where is the inertia matrix of the spacecraft and T is the control torque resolved in (see ). Three magnetic coils aligned with the axes equip the spacecraft. Thus, the following magnetic attitude control torque is created
In the previous expression, is the column matrix of the magnetic dipole moments for the three coils, and is the Earth magnetic field at spacecraft resolved in body frame (see , Section 7.4.1). Let be the Earth magnetic field at spacecraft expressed in inertial frame . Observe that changes in time, at least because of the motion of the spacecraft along the orbit. Hence,
showing the explicit dependence of on both q and t. The previous equations grouped together lead to the following system
where is the control variable. Let us characterize the dependence on time of , which is equivalent to characterize the time dependence of . Assume a circular orbit with radius R. Through the adoption of the so-called dipole model of the Earth magnetic field (see , Appendix H), obtain:
In the previous equation, is the total dipole strength, is the position of the spacecraft resolved in , and is the column matrix of the direction cosines of . The coordinates of vector are the direction cosines of the Earth magnetic dipole expressed in which are written as follows:
where is the coelevation of the dipole, deg/day is the average angular velocity of the Earth, and is the dipole right ascension at . Set Wb m and as reported in .
Eq. (7) shows that to characterize the time dependence of , one needs to determine an expression for which is the spacecraft position vector resolved in . Define the orbital plane coordinate system ap, bp whose origin is at the Earth center and with ap axis coinciding with the line of nodes. Then, the center of mass of the satellite is positioned in
where is the orbital rate and is spacecraft argument at time . It is possible to determine the mass center coordinates in the inertial frame from (9) by means of a rotation matrix which depends on the inclination of the orbit and on the right ascension of the ascending node (RAAN) (see , Section 2.6.2). By inserting the expressions of the latter coordinates into (7), we obtain an explicit expression for .
Since for where (see (1)), then the goal is designing control strategies for so that and . Reference  proposes the following stabilizing proportional derivative (PD)-like control law, which is a modification of those in [12, 13]:
In the previous equation, is the saturation limit on each magnetic dipole moment, “sat” denotes the standard saturation function, and and are parameters. As shown in , if the orbit’s inclination is not too low, there are large ranges for the values of and which lead to local exponentially stability of equilibrium for the closed-loop systems (6) and (10).
In (11), is an internal state of the controller, are all parameters, and was introduced in (3). Note that the previous equation describes an attitude feedback, since it requires only the measure of attitude and 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 , if the orbit’s inclination is not too low, there are large ranges for the values of , , which lead to local exponential stability of equilibrium 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 , we define the settling time for each component , with , as:
In other words, this is the time needed for to permanently remain under . Value can be set depending on the desired final value of qi. Finally, we can define the settling time ts for the whole quaternion q as the time required by the slowest component of , hence
Now, having set the spacecraft initial conditions to specific values, one can determine the values of the parameters that minimize the settling time ts for 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 and for the gains and . 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 ts on the initial conditions , , and (see (9)) is explicitly indicated,
and the feasible set is . Even though ts obviously depends on and and on initial conditions, it is not possible in practice to express this relation in analytical form. Moreover, ts is discontinuous with respect to , and initial conditions See for example in (, 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 , 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 worst initial 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 and . 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 S includes 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 can 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 in spherical coordinates :
The dependence of on can 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 . 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 , and evaluates the objective function in their central points. The computation of the objective value allows to identify potentially optimal hyperrectangles within . 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 , which can be described as follows. After the selection of a starting point, it cyclically determines a feasible and good descent direction for the objective function, and then it performs a sufficiently large step along such direction. To find the good feasible descent direction, at each iteration , the algorithm computes the local behavior of the objective function at the variation of the -th variable. If a move of length along gives a feasible point where the function is reduced enough, the algorithm applies a linesearch technique to find the stepsize . Otherwise, the algorithm tries the same operations in the opposite direction . If both and are 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  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 as belonging to a feasible set (in our case and ) and the set of initial conditions as belonging to a feasible set (in our case and , since has three components). Moreover, we now use to denote the objective function (in our case ). Its analytical expression is not known; however its value can be computed via a software simulation.
The upper-level minimization problem must 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 obtained for a number maxeval_ext of points . Let be one of them. The evaluation of requires however to solve one lower-lever maximization problem . 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_int total 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
Input: A computable by means of a software simulation for with and with .
Output: A robust optimum value .
Solve the upper-level problem , with by using DIRECT for maxeval_ext evaluations of , and return . The -th evaluation works with point
Solve the lower-level problem
4. Computational results for the case of PD-like control
We apply the above approach to solve the numerical example studied in . The inertia matrix of the spacecraft is equal to ; the saturation level for each magnetic dipole moment is given by . The orbit has an inclination of and 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 for . Simulation shows that T for the considered orbit; since , it follows that . Assume that the resolution of the attitude sensor is in terms of quaternion component. In this case, an upper bound for can be computed by enforcing that each component of the term (see (10)) does not exceed the saturation limit when and are orthogonal, and for one index , and for the other ones. This is achieved if is smaller than ; thus, we set . Assuming that the attitude rate sensor has a resolution of rad/s, by a parallel argument it follows that an upper bound for is given by ; thus, we set .
4.1 Fixed initial conditions
Consider an initial state characterized by attitude equal to the target attitude (which corresponds to having and 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 by setting rad.
Gains and have been found by trial-and-error in  as and . 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 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 (i.e., and any and ), the same argument rad, 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 and yield a settling time .
4.2 Variable initial conditions
We now search for and using the robust optimization approach. We allow initial conditions to vary as follows: , , , , , , and . 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.
|DIRECT + DIRECT||DIRECT + SDBOX single start||DIRECT + SDBOX multistart|
When using SDBOX multistart, the starting points are , , , , 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 even 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 and , 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.
|Gains||DIRECT||SDBOX multistart||Grid search|
|eval = 40,000||eval = 40,000 (5000 × 8)||eval = 78,125 (57)|
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 (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 (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 (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 ts is 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) , indicated by :
where represents the Euclidean norm and Tf is 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 are equilibria of the closed-loop system (6) and (11), it is very unlikely that becomes at some finite time. Note that this would make not continuously differentiable with respect to the design parameters because is not continuously differentiable at .
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 ).
By introducing physically reasonable upper bounds for the design parameters, we obtain the feasible set . 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 worst spacecraft 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 . The initial conditions of the spacecraft are given by , , , and . We chose the set of their possible values as:
Note that S includes all possible initial attitudes, all possible right ascensions of the Earth magnetic dipole at time , 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 in spherical coordinates :
The dependence of on can 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 as belonging to a feasible set (in our case and ) and the set of initial conditions as belonging to a feasible set (in our case and , since possesses three components). Indicate by the function giving the objective value (in our case ). When is fixed (the initial conditions are assigned), we simply write in it; when is fixed (the design parameters are assigned), we write in it. Problem (28) can be written as
and may be tackled by a global derivative-free optimization algorithm of the type of DIRECT , 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 dense property, 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 , since design parameters have to be greater than or equal to zero. Values can typically be set to very large values based on some physical considerations. For instance, 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 Fx requires 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 using 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” would 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 . This algorithm was originally presented in  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 , this algorithm should be able to find good solutions in short times in the neighborhood of . 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
Input: A vector and a function computable by means of a software simulation for any . Values for the parameters , , , , .
Output: A solution approximating one vector in .
Normalize and grid partition the whole feasible set into a collection of hyperrectangles similarly to the initial phase of DIRECT.
For each , compute the value of the solution obtained by iterations of SDBOX in starting from its central point. This is an upper bound on the value of the best solution in and constitutes our “evaluation” of .
Take a number of hyperrectangles corresponding to the smallest of the above values.
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 .
Switch to DIRECT algorithm to continue the search in allowing function evaluations. This search can now be dense using reasonable time, and it gives a solution .
Try to improve by using iterations of a local search method, finally obtaining a solution to problem (28).
Depending on the practical case, the number of hyperrectangles in the partitioning phase (step 1) and the number of iterations in 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 depends not only on but also on the span of each , which in turn depends on . However, there is no need that all the hyperrectangles in the initial partition have the same size. Their size could be growing with the absolute values of the coordinates so that is kept smaller. The number in the selection phase (step 3) must be chosen so that remains much smaller than ; otherwise the benefits of the proposed procedure is reduced. The number of function evaluations in the standard DIRECT phase (step 5) is chosen so that the search in is dense enough. This is now possible because of the size reduction in . The post-optimization phase (step 6) can be executed with SDBOX or with CS-DFN , 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 f to 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 , and if we set to be sufficiently large, we obtain a solution which satisfies necessary conditions for a local optimum. Moreover, should approximate one of the global optima, because, given the dense search in , it should provide one of the global minima of , and the regions of , 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 g such that its value at the generic point 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 g corresponding to different points of . Let be one of them; then the evaluation of requires the solution of one lower-level maximization problem . Thus, the lower-level problem must be solved up to maxeval_ext times. 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 Fy does 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 , , and . 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 for problem (33). We report here the whole procedure:
Procedure 3: Solve min-max combining global and local search
Input: A computable by means of a software simulation for any and . Values for the parameters maxeval_ext, maxeval_int, .
Output: A robust solution approximating one vector in .
Solve the upper-level problem , with ,
by using Procedure 2 with maxeval_ext total evaluations of , and return .
Given , the evaluation of is performed by the internal loop.
Take , and solve the lower-level problem
by using multistart local search with starting points,
performing evaluations of for 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 . The spacecraft inertia matrix is , and the saturation level for each dipole moment is . The orbit inclination is , and the orbit altitude is 450 km; the right ascension of the ascending node is equal to 0. Upper bounds are selected as .
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  by trial and error search achieves a value of ITAE=, while the vast majority of the solutions have the ITAE limit value of about . This upper limit is related to the value of Tf in 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 Tf is 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 K with, 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 that 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 . Note that this solution is also considerably better than the best solution obtained by trial and error in days of work.
|DIRECT 50,000 + CS-DFN 1000|| = 913405022.139|
|1,142,470,478.101||262,000 + 5460 s|
|DIRECT 100,000 + CS-DFN 1000|| = 500798437.500|
|1,129,234,873.703||565,200 + 5040 s|
|Procedure 2: Combining global and local strategies|| = 246494.579020|
|8,021,573.4077||201,500 + 15 + 4 s|
The promising region of the feasible set K has been identified with the previously described probing technique. The values of and p are 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 = 10, the evaluation of each hyperrectangle would require roughly 10 s. Thus, to finish within 2 or 3 days of computation, we select , which is obtained by making 64 partitions on the domain of , 64 partitions on the domain of , 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 actually 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 , 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 , as in Step 5 of Procedure 2, after 3007 iterations and only 15 s, we obtain the solution:
whose value is ITAE = 9.072. Note that the volume of the set K* is considerably smaller than that of K: it is only of the volume of K. To better appreciate the difference, consider that an exploration of K with the same degree of density used on K* would require s that roughly corresponds to more than 7 years, if the simulation times on K were 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 , the solution (39) is improved in 4 s to ITAE = with the solution:
We also tested in this post-optimization phase the local search SDBOX, which moves only along the coordinate directions . With 5000 additional iterations of SDBOX that requires 11 s, solution (39) is improved to the following new solution, which has ITAE = 9.062
Solution (40) corresponds to of 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 s. Then, we search for a lower bound 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 . Consider now simple rotations about each single body axis, and for each simple rotation, compute lower bounds , , 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 . Value can be computed using the equation which describes rotation about the body axis, which is given by
where is the roll angle. The amplitude of torque is limited by an upper limit , which can be found using (4). Indeed, numerical simulations show that T. Since (see (5)), then . Moreover, each component of is bounded by and then . Thus, N m. Next, the minimum time to bring the state of system (42) subject to the constraint , from the initial state to the final state , is given by (see , Section 7.2)
Similar considerations hold for rotations about and axes leading to
Then, . 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 , 20 on that of , and 10 on that of . We therefore perform 40,000 evaluations of for the identification of the new and then another 3000 + 1000 evaluations of to solve the problem on this , for a total of .
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 , , and . 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 . Run times are greater in the rest of ; however, we try to keep them under control by allowing to exit the internal loop when the value of is large enough to reach the limit value for ITAE of about . 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 for the initial conditions (37) instead of ITAE of solution (40). However, this solution is a robust solution: by varying the initial conditions in , the worst value that can be obtained is ITAE that is still considerably better than average solutions, whose vast majority has the limit value for ITAE of . As a comparison, the worst value achievable by varying the initial conditions in for solution (40) is ITAE that 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 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 K allowing 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.
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.