# Averaging Methods for Design of Spacecraft Hysteresis Damper

Authors: Ricardo Gama; Anna D. Guerman; Ana Seabra; Georgi V. Smirnov
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This work deals with averaging methods for dynamics of attitude stabilization systems. The operation of passive gravity-gradient attitude stabilization systems involving hysteresis rods is described by discontinuous differential equations. We apply recently developed averaging techniques for discontinuous system in order to simplify its analysis and to perform parameter optimization. The results obtained using this analytic method are compared with those of numerical optimization.


---

## Body

## 1. Introduction

Dampers that use magnetic hysteresis rods to dissipate the energy of undesired angular motions occurred during deployment or caused by perturbations are used in attitude control systems of small satellites since 1960s [1]. Mathematical modeling of such systems is quite a difficult task since the majority of existent hysteresis models result in differential equations with discontinuous right-hand side.Analysis of dynamics for attitude control systems with magnetic hysteresis dampers and optimization of their parameters have been done in [2, 3], and the results of these studies have been implemented in real missions [4, 5]. However, these studies lack an accurate theoretical basis for application of averaging methods to such problems.Recently, an adequate mathematical approach has been developed by the authors in [6]. Now we can address a complete mathematical theory for attitude stabilization systems with hysteresis.Consider a differential equation(1)x˙=ϵf(t,x,u),x∈Rn,t≥0
describing a mechanical system with stabilizer. Here u∈U⊂Rk is a parameter. It is assumed that 0≈f(t,0,u) for all t≥0 and u∈U; that is, the velocity of the system near the origin is small. Here we do not assume that zero is an equilibrium position of system (1). The parameter u should be chosen to optimize, in some sense, the behavior of the trajectories. The choice of this parameter can be based on various criteria. Obviously, it is impossible to construct a stabilizer optimal in all aspects. Consider, for example, a linear controllable system. The pole assignment theorem guarantees the existence of a linear feedback yielding a linear differential equation with any given set of eigenvalues, so one can choose a stabilizer with a very high damping speed. However, such a stabilizer is practically useless because of the so-called peak effect (see [7, 8]). Namely, there exists a large deviation of the solutions from the equilibrium position at the beginning of the stabilization process, whenever the module of the eigenvalues is big.The aim of this paper is to develop effective analytical and numerical tools oriented to optimization of stabilizer parameters for passive attitude stabilization system with hysteresis rods.Throughout this paper, we denote the set of real numbers byR and the usual n-dimensional space of vectors with components in R by Rn. We denote by 〈a,b〉 the usual scalar product in Rn and by |·| a norm. By B we denote the closed unit ball, that is, the set of vectors x∈Rn satisfying |x|≤1. The transpose of a matrix A is denoted by A*. The set of positively definite symmetric n×n-matrices is denoted by M(n). If P and Q are two subsets in Rn and λ∈R, we use the following notations: λP={λp∣p∈P}, P+Q={p+q∣p∈P,q∈Q}. The convex hull and the closure of a subset S⊂ℝn are denoted by coS and clS, respectively. The Hausdorff distance between two sets A1,A2⊂Rn is defined as
(2)h(A1,A2)=min{h≥0∣A1⊂A2+hB,A2⊂A1+hB}.
We denote by 𝒮[0,T](F,x0) the set of solutions to the Cauchy problem x˙∈F(t,x),  t∈[0,T],  x(0)=x0, and by ℛ[0,T](F,x0)={x(T)∣x(·)∈𝒮[0,T](F,x0)} the reachability set. We use also the notations 𝒮[0,T](F,C)=∪x0∈C𝒮[0,T](F,x0) and 𝒮[0,T](F)=𝒮[0,T](F,Rn). The closed unit ball in the space of continuous functions f:[0,T]→Rn with the uniform norm, C([0,T],Rn), is denoted by ℬ. The set of locally integrable functions f:[0,∞[→Rn is denoted by L1loc([0,∞[,Rn). The upper limit of a set-valued map F:Rn→Rm is given by
(3)limsupx′→xF(x′)={v=limn→∞vk∣(xk,vk)∈grF,xk⟶x},
where grF stands for the graph of the set valued map F.


## 2. Statement of the Problem

Consider dynamics of a satellite in a circular geocentric orbit. The satellite is equipped with a gravity-gradient attitude control system that includes a number of magnetic hysteresis rods as a damper. The spacecraft's equations of motion [2] can be represented in the normalized form. Denote by x(t,x0,u) the solution to the respective Cauchy problem
(4)x˙=ϵf(t,x,u),x∈Rn,t∈[0,T],x(0)=x0,
where u is a parameter from a compact set U⊂Rk. Define the functions
(5)φi(u)=maxt∈Δimaxx0∈Bi|x(t,x0,u)|i,i=0,m¯.HereΔi⊆[0,T] are compact sets, |·|i are norms in Rn, and Bi={x∈Rn∣|x|i≤bi}. Consider the following mathematical programming problem:
(6)φ0(u)⟶min,φi(u)≤φ-i,i=1,m¯,u∈U.Many problems of stabilization systems' parameters optimization can be written in this form (see [9]). For example, the minimization of the final deviation can be formalized as follows:
(7)maxx0∈B|x(T,x0,u)|⟶min,u∈U,
and the minimization of the maximal deviation of trajectories satisfying certain restrictions at the final moment of time has the form:
(8)maxt∈[0,T]max|x0|=1|x(t,x0,u)|⟶min,max|x0|=1|x(T,x0,u)|≤δ,u∈U.


## 3. Averaging for Discontinuous Systems

The averaging method is one of the most used methods to analyze differential equations of the form(9)x˙=ϵf(t,x),
appearing in the study of nonlinear problems. The idea behind the averaging method is to replace the original equation by the averaged one:
(10)x˙=ϵf-(x)=ϵlimT→∞1T∫0T‍f(t,x)dt.This equation is simpler and has solutions close to the solutions of the original equation. A rigorous justification of the method is given by Bogolyubov's first theorem containing an estimate for the distance between the solutions of the exact and averaged systems on large time intervals [10]. The Samoilenko-Stanzhitskii theorem [11, Theorem  2] which is a generalization of Bogolyubov's second theorem shows that asymptotic stability of the zero equilibrium position of averaged (10) implies that the solutions to original (9) are close to zero on the infinite time interval.For several models of systems with hysteresis, including the passive attitude stabilization systems, the functionf(t,·) appearing in (9) is discontinuous (see, e.g., [3]), and the classical notion of solution and the classical averaging method cannot be used. For such systems Filippov proposed a generalized concept of solution, rewriting problem (9) as a differential inclusion
(11)x˙∈ϵF(t,x),x(0)=x0,
where x→F(·,x) is an upper semicontinuous set-valued map obtained from f(t,·) by Filippov regularization [12, 13]. The use of this concept of solution makes it necessary to generalize the averaging method to differential inclusions. Many results extending Bogolyubov's first theorem to differential inclusions have been obtained (see, e.g., [14, 15]). In the case of Lipschitzian differential inclusions, the problem has been completely solved by Plotnikov [14]. Averaging results for inclusions with upper semicontinuous right-hand side have been obtained by Plotnikov [15] under conditions of Lipschitz continuity of the averaged inclusion and for inclusions with a piecewise Lipschitzian right-hand side. Recently [6], an averaged differential inclusion has been introduced allowing one to prove extensions of Bogolyubov's first theorem and of the Samoilenko-Stanzhitskii theorem for upper semi-continuous differential inclusions and, as a consequence, for discontinuous dynamical systems. Here we outline the main results from [6].LetF:R×Rn→Rn be a set-valued map. Set
(12)I(t1,t2,x,δ)={∫t1t2v(t)dt∣v(·)∈L1loc([0,∞[,Rn),v(t)∈F(t,x+δB)∫}.We denote byF-δ(x) the convex hull of the map
(13)Φ-δ(x)=limsupθ↑1limsupT→∞1(1-θ)TI(θT,T,x,δ)
and define the averaged differential inclusion as
(14)x˙∈F-(x)=⋂δ>0F-δ(x).Note that under Lipschitz condition this map coincides with(15)F̿=limT→∞1T∫0TF(t,x)dt,
if the limit exists in the sense of Hausdorff distance (see [6]).Assume that the following conditions are satisfied:(C1)
clcoF(t,x)=F(t,x), for all (t,x)∈R×Rn;(C2)
the set-valued mapF(t,·) is upper semi-continuous;(C3)
for anyx there exists measurable selection of F(t,x), that is, there exists f(t,x)∈F(t,x) such that t→f(t,x) is measurable for all x;(C4)
there exists a nonnegativeb(·)∈L1loc([0,∞[,R) such that F(t,x)⊂b(t)B for all (t,x)∈[0,+∞[×Rn;(C5)
there exists the limit(16)b=limT→∞1T∫0Tb(t)dt.Under these conditions, the following version of Bogolyubov's first theorem is true.Theorem 1.
LetT>0 and let F:R×Rn→Rn be a set-valued map satisfying conditions (C1)–(C5). Let C∈Rn be a compact set. Then for any η>0 there exists ϵ0>0 such that for any ϵ∈]0,ϵ0[ and any solution x(·)∈𝒮[0,T/ϵ](ϵF,C), there exists a solution x-(·)∈𝒮[0,T/ϵ](ϵF-,C) satisfying
(17)|x(t)-x-(t)|<η,t∈[0,Tϵ].Set(18)Gϵ(τ,y)=F(τϵ,y),G0(y)=F-(y).Next theorem is an extension of the Samoilenko-Stanzhitskii theorem.Theorem 2.
LetF:R×Rn→Rn be a set-valued map satisfying conditions (C1)–(C5). Assume that y=0 is an asymptotically stable equilibrium position of the differential inclusion y˙∈G0(y). Then for any η>0, there exist ϵ0>0 and δ>0 such that 𝒮[0,∞[(Gϵ,δB)⊂ηℬ, whenever ϵ∈]0,ϵ0[.The last theorem shows that if the averaged inclusion has zero as its asymptotically stable equilibrium position, the trajectories of the original inclusion stay in the vicinity of the origin providedϵ>0 and |x0| are sufficiently small.If the averaged inclusion has a special form, we can go further and make some conclusion on the detailed behaviour of the trajectories of the original system. Assume that the averaged inclusion has the form(19)x¯˙∈ϵ(A(u)x-+P-(x-,u)),
where P-(x-,u)⊂c|x-|2B, c>0, the real parts of the matrix A(u) eigenvalues are negative for all u∈U, and the function u→A(u) is continuous for all u∈U. If γ0>0 is sufficiently small, then the set of solutions to the Lyapunov inequality (see [16]) for the matrix A(u),
(20)ℒ(u)={(γ,V)γ≥γ0,V∈M(n),AV+A*V≤-2γV},
is nonempty and compact for all u∈U. Let (γ,V)∈ℒ(U). Denote by |x-|V the Euclidean norm defined by |x-|V=〈x-,Vx-〉. There exist positive constants c1 and c2 satisfying
(21)c1|x-|≤|x¯|V≤c2|x-|,
whenever (γ,V)∈ℒ(U) for some γ.Theorem 3.
Letδ>0,  u∈U, and (γ,V)∈ℒ(U). There exists ϵ0(δ) such that for all ϵ∈]0,ϵ0(δ)[ the condition |x0|V<δ<c12γ/c implies the inequality |x(t,x0,u)|V<3δ/2.This theorem shows that the behavior of the trajectoryx(t,x0,u) can be characterized in terms of the pair (γ,V). The parameter γ is responsible for the damping speed of the process, while the form of the ellipsoid {x∣〈x,Vx〉≤1} describes the amplitude of the deviation of the trajectory from the origin. The aim of parameter choosing can be formulated as follows: maximal value of γ and maximal sphericity of the ellipsoid {x∣〈x,Vx〉≤1}. The latter property guarantees minimal overshooting of the damping process and, as a consequence, the largest region of applicability of the approximation obtained via averaging.


## 4. Choosing Passive Magnetic Stabilizer Parameters

The in-plane oscillations of a satellite moving along a polar circular orbit and equipped with a passive gravity-gradient attitude stabilization system with one hysteresis rod are described by the equation(22)α¨+ω2α=ϵf,
where α is the pitch angle of spacecraft, ϵ is small parameter proportional to the rod's volume, and the force f is given by
(23)f=W(Hτ)(H1e3-H3e1-α(H1e1+H3e3)).Here(24)Hτ=H1e1+H3e3+α(H1e3-H3e1)
describes the projection of the geomagnetic field on the rod axis,
(25)W(Hτ)=Hτ-κ2signH˙τ
is the hysteresis function, κ corresponds to the coercive force,
(26)H1=cost,H3=-2sint,
and t is the argument of latitude of the current point of the orbit [2]. The vector
(27)(e1,e3)=(cosθ,sinθ),θ∈[0,π],
describes the orientation of the hysteresis rod in the satellite body. Equation (22) is equivalent to the system
(28)α˙=β,β˙=-ω2α+ϵf.After the change of variables(29)α=acosωt+bsinωt,(30)β=-aωsinωt+bωcosωt,
one arrives at the system
(31)a˙=-ϵωfsinωt,b˙=ϵωfcosωt.Theorem 4.
Assume thatω is an irrational number. Then the averaged system for (31) is
(32)a¯˙=-ϵ2ω(pa-+qb-)+ra-(a-,b-,θ),b-˙=ϵ2ω(qa--pb-)+rb-(a-,b-,θ),
where
(33)p=9κωe12e32π(1+3e32)3/2,q=32(e12-e32)+6κe1e3π(1+3e32)3/2,|
r
a
-
|
=
O
(
a
-
2
+
b
-
2
), and |rb-|=O(a-2+b-2).Obviously we have(34)(1001)(-p-qq-p)+(-pq-q-p)(1001)=-2p(1001).
Therefore, we see that the linearization of the averaged system always has a Lyapunov function V=a-2+b-2, and the damping speed is determined by the value of p. This means that the peak effect does not take place for the linearization of the averaged system.To maximize the damping, one has to increase the total volume of the hysteresis material on board. However, it is wellknown that the efficiency of a damping rod is increased with the increase of the ratio between the rod's length and its cross-section dimension. Therefore, instead of one massive bar, the attitude control system should use several rather thin rods of the maximum length allowed by the spacecraft geometrical and system restrictions. On the other hand, to minimize the perturbation of the spacecraft angular motion caused by the damping system itself, the direction of total magnetic field in the rods should deviate as little as possible from the direction of the geomagnetic field at the current point of the orbit. Thus, in general case, one should use a system of three equal orthogonal hysteresis rods or a number of such systems. Here we consider in-plane satellite dynamics on a polar orbit, and for such purpose it suffices to analyze a pair of equal orthogonal rods.Orientation of this pair of equal orthogonal rods can be defined as(e1,e3), (-e3,e1), where e1=cosθ and e3=sinθ, and due to the system symmetry it is enough to study the interval θ∈[0,π/2]. If the satellite is equipped with several identical hysteresis rods, the corresponding nonlinear system is
(35)a˙=-ϵω(f1+f2+⋯)sinωt,b˙=ϵω(f1+f2+⋯)cosωt.Here the termsf1,f2,…, describe the interaction of the respective rod with the geomagnetic field. For a couple of equal orthogonal rods and for small deviation from the origin, the forces f1 and f2 are given by
(36)f1=W(Hτ,1)(H1e3-H3e1-α(H1e1+H3e3)),f2=W(Hτ,2)(H1e1+H3e3-α(-H1e3+H3e1)),
respectively. Here
(37)Hτ,1=H1e1+H3e3+α(H1e3-H3e1),Hτ,2=-H1e3+H3e1+α(H1e1+H3e3),W(Hτ,j)=Hτ,j-κ2signH˙τ,j,j=1,2.For this case, the first approximation of the averaged system takes the form(38)a-˙=-ϵ2ωPa-,b-˙=-ϵ2ωPb-,
where
(39)P=9κωe12e32π(1(1+3e32)3/2+1(1+3e12)3/2).An easy calculation shows that the optimal value ofθ is π/4, so e1=e3=2/2. In the next section we numerically analyze the validity of the previous analytical study.


## 5. Numerical Simulations

We approximate problem (6) by the following problem:
(40)φ-0⟶min,|x~(tki,xji,u)|i≤φ-i+ε,i=0,m¯,u∈U,
where t0i=0,  tki∈Δi,  xji∈Bi,  j=1,J¯, and
(41)x~(tk+1i,xji,u)=x~(tki,xji,u)+τf(tk,x~(tki,xji,u),u),τ=tk+1i-tki,k=0,N¯,
is the Euler approximation for the solution x(·,xji,u). Note that this is a hard problem because of the discontinuity of the system. It is necessary to consider very fine partition of the time interval in order to get a good approximation of the solutions to the discontinuous differential equation. For smooth right-hand sides the number of points in the mesh can be significantly reduced (see [9]).Letϵκ>0 be small enough. We consider a satellite with two equal orthogonal hysteresis rods. A typical trajectory of the system is shown in Figure 1.A typical trajectory forϵ=0.25, κ=0.1, ω=0.949.
(a)
Phase plot ofb versus a
(b)
Time evolution of amplitudea2+b2 of pitch oscillations (29)Note that the oscillations of the trajectory do not allow one to characterize its damping speed using the norm at the final moment of time. For this reason, we numerically solve the following problem:(42)maxx0∈δ0Bmaxt∈[T-Tp,T]|x(t,x0,u)|⟶min,u∈U,
where Tp≪T is an interval corresponding to the period of oscillation of solutions. The minimization is done using multistart Nelder-Mead method. For ϵ=0.25, κ=0.1, ω=0.949, δ0=1, T=300π, and Tp=20π, the fastest damping speed is observed for values of θ very close to π/4, that is, to the expected value. Next we consider the problem
(43)maxx0∈δ0Bmaxt∈[0,T]|x(t,x0,u)|⟶min,maxx0∈δ0Bmaxt∈[T-Tp,T]|x(t,x0,u)|≤δ1,u∈U.The results of minimization forϵ=0.25, κ=0.1, ω=0.949, δ0=1, δ1=0.1, T=300π, and Tp=20π show that the problem has many local minima but with the values of the functional very close to 1; that is, the nonlinear system with two orthogonal hysteresis rods has no overshooting for all values of θ. This allows one to conclude that the value θ=π/4 is the best one. It guarantees high damping speed and does not cause peaking.


## 6. Conclusions

This paper is dedicated to the problem of parameter optimization for a gravitationally stabilized satellite with magnetic hysteresis damper. Its motion is described by differential equations with discontinuous right-hand side. The discontinuity is the principal obstacle in the application of the averaging method. Our recently obtained results on averaging of discontinuous systems are applied now to rigorously justify the use of this method for a satellite with hysteresis rods. We consider here the simplest case of in-plane oscillations on a polar circular orbit. Theorem3 shows that the behavior of the system can be characterized in terms of Lyapunov function which can be chosen in order to guarantee the best properties of damping process. Further study is performed numerically, and the simulations are in excellent agreement with the analytical results, confirming also previous studies on hysteresis damping of satellite pitch oscillation in gravity-gradient mode.Our results can also be applied to rigorously justify the use of averaging techniques in analysis of other engineering problems involving differential equations with discontinuous right-hand side.



---

*Source: 483457-2013-06-13.xml*