Energy-efficient control of nanosatellites during distributed region formation flying

Reza HAGHIGHI,Chee Khiang PANG

1.Department of Electrical and Computer Engineering,National University of Singapore,Singapore 117583

2.Engineering Cluster,Singapore Institute of Techno logy,Singapore 138683

1 Introduction

Advances in nanotechnology have led to the development of nanosatellites capable of performing missions that were previously done by large and costly satellites[1].The evolution of low-cost nanosatellites have helped m ost of the universities around the world to build their own nanosatellites.Many of the nanosatellites were built in a standard format known as a CubeSat,a10 cm cubewith amass of about1kg[2].CubeSatscan be made up from one to several units.The concept of CubeSat was originally developed at the Stanford University’s Space System s Developm ent Laboratory and the first CubeSats were launched in 2003[3].Lately,a project called QB50 is funded by the European Union to launch anetwork of 50 CubeSats[4].Although hundreds of nanosatellites have been launched into space,m ost of them lack propulsion system s.Hence,once nanosatel-lites reach their designated orbits,they passively spin in space.Nanosatellites deployed in the low Earth orbit(LEO)are affected by the atm ospheric drag which could slow down their motion.This will cause the nanosatellites to de-orbit,decrease their altitude and eventually burn up in the atmosphere.

Recently there has been increasing interests on development of propulsion system s for nanosatellites.Generally,propulsion system s utilized for satellites can be classified into cold gas,chemical thrusters,and electric thrusters.Although chemical propulsion system s can be easilymanufactured and they provide higher thrust,they are not suitable choice for nanosatellites,as nanosatellites with explosive pressurized tanks are not allowed to piggyback on rockets that carry them into space.In 2013,NASA’s Game Changing Development(GCD)program awarded to three team s(that include Busek Co.Inc.,Massachusetts Institute of Technology,and the NASA’s Jet Propulsion Laboratory)to develop small,highly efficient propulsion system s specifically designed to enable satellites like CubeSats to orient them selves to the desired attitudes,m aneuver,and even change their own orbits.The team s are developing M icrofluidic Electrospray Propulsion(MEP)thrusters[5]based on three concepts including High Aspect Ratio Porous Surface(HARPS),Scalable ion Electrospray Propulsion System(S-iEPS),and indium m icrofluidic electrospray propulsion system.The above m entioned developm ents on propulsion system s for nanosatellites help to actively control the nanosatellites in space.

Formation flying of nanosatellites has been the major interest of researcherson space studies in the pastyears.The very first attempt on formation flying of nanosatellites goes back to SNAP-1 with a m ass of 6.5 kg launched in 2000[6].In 2014,the CanX-4 and CanX-5(developed by the University of Toronto)were launched together to demonstrate the on-orbit formation flying.The formation flying is controlled by the Canadian Nanosatellite Advanced Propulsion System(CNAPS)which is a liquefied cold-gas thruster system[7].In 2015,the Aerocube 8 A and B(developed by the Aerospace Corporation of El Segundo)were lunched to demonstrate de-orbit maneuvers and formation flying by using the Scalable ion-Electrospray Propulsion System(SiEPro)[8].

To achieve formation flying,various techniques have been developed.In[9-11],the satellite formation flying was addressed by considering the attitude coordination control problem.In[12,13],the finite-time attitude synchronization has been investigated.Another strategy for the formation flying of satellites is using the relative motion control[14-16].Chung et al.[17]considered the synchronization of both attitude and translational motion of a group of satellites.Ahn et al.[18]considered an iterative learning control scheme to ensure trajectory-keeping between the leader and follower satellites.Morgan and Chung[19]presented a real-time optimal control algorithm for a linearized model of a swarm of satellites.The above mentioned methods on the satellite formation flying based on relative motion control rely on the assumption of availability of position of leader in its orbital slot to all followers.However,this assumption may be unrealistic as the leader may encounter various disturbances during flying in the LEO such as drag force,hence,the prep lanned leader’s information may not be accurate.In addition,as the number of nanosatellites increases the methods that rely on the virtual structure approach for the formation control cannot provide satisfactory result as the constraint relationships among the members become more complicated.An efficient method to deal with formation of large number of nanosatellites is using the region formation control[20].In the region formation control,a region is defined based on a desired reference such that satellites are derived into the specified region and remain inside during movement.

In this paper,we present an energy-efficient method for distributed region formation flying of nanosatellites under directed the switching communication topology.Here,it is assumed that only a portion of nanosatellites have know ledge of the leader’s orbital elements.Hence,a distributed estimator is proposed so that other follower nanosatellites estimate the position of the leader in itsorbital slot based on the switching communication topology.Using thees timated position of leadernanosatellite,an optimal region following formation method based on the receding horizon control(RHC)is developed.An algorithm is presented to obtain the optimal input for the distributed region-based formation flying.The simulation results are presented to illustrate the effectiveness of the proposed methodology.The contribution of the paper is two folds:first,a distributed estimator is developed for nanosatellites under directed switching communication topology to estim ate the position of the leader satellite.Existing methods on formation flying are based on the assumption of the availability of the leader’s orbital elements to all followers.The sta-bility analysis is performed to show the convergence of the proposed estimator.Second,an optimal control methodology is proposed based on the RHC approach to achieve the distributed region formation flying of nanosatellites.Note that,a prelim inary version of this paper was presented in ICCA 2016[21].This paper presents the extended version with the following extensions:i)The convergence analysis of the proposed estimator is perform ed.ii)The stability analysis of the proposed energy-efficient formation flying method is demonstrated.iii)An algorithm is presented to obtain the optimal input for the distributed region-based formation flying.

2 Relative motion equations for satellite formation flying

Fig.1 The ECI frame and the LVLH frame.

In order to achieve the satellite formation flying,a reference satellite(which can be virtual)is assigned to follow the desired orbit while the follower satellites maintain the desired formation with respect to the reference satellite.The motion equations of satellites are defined with respect to two coordinate fram es as illustrated in Fig.1.First,the Earth centered inertial(ECI)frame shown by(X,Y,Z)and the local vertical local horizontal(LVLH)frame shown by(x,y,z)which is attached to the reference satellite.The motion of a satellite within an orbit is com p letely described by six orbital elements,which areasem i-major axis,eeccentricity,ioorbit inclination,ωpargument of Perigee,Ω the right ascension of the ascending node,and ν true anomaly.As illustrated in Fig.1,the motion equation of the reference satellite on the desired orbit is defined based on the ECI frame,and the relative motion equations of the follower satellites are defined based on the LVLH fram e.

Let|ro|and|h|denote the geocentric distance and the magnitude of the angular momentum,respectively.The angular velocities and angular accelerations of the LVLH frame are expressed as follow s(for more details refer to[22]):

and

whereiois the orbit inclination of the reference satellite,θ=ωp+ν is the argument of latitude,vxis the radial velocity,andkJ2is a constant defined as

such thatJ2is the second zonal harmonic coefficient of the Earth,μ is the Earth’s gravitational constant,andReis the Earth’s equatorial radius.

2.1 Dynam ic equation of the reference satellite under the J2 perturbation

The motion equation of the reference satellite considering theJ2perturbation can be expressed as follow s(for more details refer to[22]):

whereandf(?)is defined as

Property 1Since parameters|ro|and|h|are lower bounded,the functionf(?)satisfies the Lipschitz condition on a set Π with the Lipschitz constant ?,i.e.,

2.2 Dynam ic equations of follower satellites under the J2 perturbation

Consider a group ofNsatellites such that the position of satelliteiwith respect to the reference satellite in LVLH fram e is specified by ρi=[xi yi zi]T.The motion equations of the follower satellites with respect to the reference satellite in LVLH frame is expressed as follows(for more details refer to[23]):

whereuiare control inputs fori=1,...,N,andS(?)and G(ρi,?)are defined as

and

such that

whererizis the projection of ρion theZaxis of the ECI frame,obtained as follows:

2.3 Dynamics of propulsion system

Here,we introduce the utilized propulsion system.Nanosatellites with an approximate mass of 5 kg require 50μN to 5m N thrust.An applicable propulsion system for nanosatellites is electric propulsion system.Next,w e provide the relationship between the acceleration of nanosatellites and the consumed power by the electrospray propulsion systems.

An important factor that indicates a thruster efficiency is called the specific impulse which is defined as follows:

whereFiis the thrust,andis the mass flow rate of theith follower nanosatellite,andg0is the gravitational acceleration constant.The specific impulse for electrospray propulsion system s is defined as follow s:

wheremqiis the charge to mass ratio andVaiis the applied voltage.Using(11)the trust is obtained as follow s:

The m ass flow rate of ions is related to the ion beam current,Iai,and the charge to mass ratio as

The electric power efficiency of the thruster is defined as the ratio of the beam power,Wai,divided by the input powerW i,

From(13)-(15)the thrust of an electrospray propulsion system can be computed as follow s:

3 Distributed estimation of the reference orbital elements

In satellites formation flying for few number of satellites,we can assume that all of the follower satellites have know ledge about the reference orbital parameters.However,as the number of satellites increases,this assumption is no more realistic.Hence,here we assume that the follower satellites estim ate the reference orbital parameters in a distributed manner.Moreover,during the information exchange among the satellites,the communication topology may changes due to possible m essage dropouts in the communication channels,communication link failures or even failures of some satellites.Hence,it is reasonable to assume that the communication network has a switching topology[24].

Let Gσ(t)=(V,Eσ(t))denotes the dynamic communication graph among satellitessuch that V={1,...,N}is thenode set,Eσ(t)is theedge set,andσ(t):[0,+∞)→ P isa right-continuouspiecewise constantswitchingsignal where P is thesetofallpossible interconnection topologies.The laplacian matrix associated with the communication graph Gσ(t)is represented by Lσ(t)such that

where(j,i)∈ Eσ(t).

Lemm a 1[25]Lettj,j=0,1,2,...be an infinite sequence of nonem pty,bounded and contiguous time intervals such thatt0=0,tj+1-tj≥TwhereT＞0.Suppose that continuous functionV(t):satisfies the following conditions:

i)V(t)is lower bounded;

ii)is differentiable and non-positive on each interval[tj,tj+1);

iii)is bounded over[0,+∞)in the sense that there exists a positive constant such that

Thenast→∞.

Assum ption 1Over each intervalfor some integerp∈P;in addition,the dynamic graph Gσ(t)is directed and every node is reachable from the reference node 0.

Letbe the estimation of ? defined by(4)by theith follower nanosatellite.The proposed distributed estimator is expressed as follow s:

wherekais a positive constant,f(·)is defined by(5),andbi＞0 if satelliteihas an access to the reference satellite andbi=0 otherw ise.The dynamics of the estimation error can be obtained as

whereThe closed-loop error dynamic can be written as

where1N=[1 ···1]T∈ RNandsuch that

Theorem 1Under Assum ption 1,the proposed estimator(20)yields the convergence of

ProofUsing Lemma2 in theappendix,and Assum ption 1,one can concludematrix Hσ(t)ispositivestable in each non-switching interval.Consequently,there exist positive definitematricesPσ(t)andQσ(t)which satisfy the Lyapunov equationwhereandQσ(t)are constants in each non-switching interval and possibly time-varying in switching instances.

Consider the following Lyapunov function candidate:

whereI6∈R6is an identity matrix.The Lyapunov function candidateVeis differentiable at any time except for switching instants.We first show thatat any non-switching instants.DifferentiatingVeand substituting the closed-loop equation(22)into it,yields

Using Property 1,we obtain

w here λmin(·)and λmax(·)denote the minim um and maximum eigenvalue of a matrix,respectively.We can assureis negative by choosing(note that we do not need the exact know ledge of matricesPσ(t)andQσ(t).An approxim ate know ledge of the bounds of the eigenvalues of them suffices.A conservative way is to choose the gainkalarge enough to satisfy the condition).Therefore,we can concludeVeis bounded and consequentlyis bounded.From(22)w e haveis bounded,therefore w e can concludeis bounded on each non-switching interval i.e.there exists a positive constant Υ such that

Invoking Lemma 1,we haveTherefore,which leads to□

Rem ark 1In Theorem 1 we assume the leader nanosatellite is globally reachable at each non-switching instance.This assumption can be further relaxed to the case where the leader is jointly reachable.Note that the leader is jointly reachable if each interval[ti,ti+1)contains a sequence of non-overlapping,contiguous subinterval where the topology graph in invariant in each subinterval;and in the graph obtained from the union of graphs of subintervals,every node is reachable from the reference node 0.We can show that the proposed estimator(20)under jointly reachability condition yields the convergence ofby using the Cauchy’s convergence criteria and Lemma 3.3 in[26].

Rep lace ? with its estimationobtained by theith follower nanosatellite,then equation(17)is written as follows:

4 Region-based formation flying strategy

An efficient method to deal with formation of swarm of nanosatellites is using region formation control[27,28].In region formation control,a region is defined based on a desired reference such that satellites are derived into the specified region and remain inside during the movement.This method has an advantage when we are dealing with a large number of nanosatellites.Here,the desired region is specified by intersection ofmsubregions,as follow s[20]:

For the satelliteithe region function,πi(ρi),is less than or equal to zero when the satelliteiis located inside or on the boundary of the desired region,otherw ise πi(ρi)is greater than zero.

Assum p tion 2are continuous scalar functions with continuous partial derivatives.

An exam p le of the region function in 3D space is

for sphere shape with the radiusr.The potential region function for the satelliteiis defined as follow s:

whereklis a positive constant,κ is a constant which should be chosen greater than 2 such that the formation potential functions at least belong to class C2.

In order to keep the minimum distance among satellites and avoid collisions,interactive potential function is introduced as follow s:

where Δρijrefers to the relative distance between the agentiand the agentj,kijare positive constants,Nirepresents the set of neighbors of theith nanosatellite,ε is defined in such a way that the interactive potential functions at least belong to class C2,and

whereaanddare positive constants.

Next,w e present tw o formation flying methods using the region following technique.At first a formation flying control based on the sliding mode approach is presented.Then,an energy-efficient method is proposed by utilizing the RHC.

5 Region formation flying based on sliding mode approach

Formation flying based on the sliding method approach hasbeen frequentlyutilized[29,30].Theexisting methods on formation flying control are based on the virtual structure.However,as the number of nanosatellites in the group increases,it is difficult to specify the desired position for each nanosatellite.Hence,a region based formation flying provides more efficient solution to the satellite formation flying problem.Here,we enhance the exiting results by introducing a smooth region based formation flying controller.

To achieve the region based formation flying control,a reference vector is defined as follow s:

where

Using the reference vector sliding vectors are defined as follow s:

By differentiating equation(36)with respect to time,is obtained as follows:

Substituting equations(36)and(37)into(27),the following equation is obtained:

The control law to achieve the distributed region formation flying is introduced as follow s:

Before stating the theorem,w e express the following definition on solutions of a system with discontinuous right-hand side.

De finition 1(Filippov solution[31]) Consider the system

where f∈Rn×R→Rnis Lebesgue m easurable and essentially locally bounded.A vector functionx(t)is called a solution of(41)on the interval[0,∞)ifx(t)is absolutely continuous and for almost allt∈ [0,∞)

where

Theorem 2Consider a group of nanosatellites represented by the relative motion equations(27),which are following the reference nanosatellite represented by(5).Suppose the communication topology among the nanosatellites is represented by a directed switching graph.Under assumptions 1 and 2,the distributed region based formation flying control introduced by(39)and the estimator(20),result in the convergence of Δ?ito zero,for alli=1,...,N,ast→ ∞,where Δ?i= Δξi+Δψi.

ProofTo examine the stability of the overall system described by equation(40),the following(smooth,time independent)Lyapunov-like candidate is proposed:

Then,we have

whereSinceVis C∞insi,ρi,and Δρij,

Using the closed-loop equations(40),and considering the skew-symmetricity of the matrixyields:

where

Using thenon-smooth version ofthe LaSalle’s Invariance principle[32,33],we have

Rem ark 2Theorem 2 show s that the sum of the region error Δξiand the interaction force Δψiis zero in steady-state for each nanosatellite.If the desired region is large enough,then by choosing a higher gain for the region following term we can guarantee that both the region error and the interaction force go to zero.

6 Energy-efficient control of nanosatellites formation flying

In this section,we present an optimal method based on the RHC to achieve the energy-efficient formation flying of a group of nanosatellites.We should note that,although the continuous-time representation of the RHC would bemore natural,the development of the RHC law is much more difficult.Moreover,the overall optimization procedure has to be continuously repeated after each extremely small sam p ling time,which is computationally intractable.In fact,the sam p ling time must be larger than the com putation time required to perform the optimization,noting that the RHC is a com putationally expensive approach.Considering the aforementioned issues,it is more reasonable to formulate the RHC in discrete-time even though the proposed estimator is developed in continuous time.

LetThe discretized model of(27)can be expressed as follows:

Next,the formulation of the decentralized RHC for the region based formation control of nanosatellites is presented.

6.1 RHC formulation for region based formation control

In the RHC the current control action is obtained for a system by solving,at each sampling time,a finite horizon optimization problem using the current state of the system as the initial state.Then the optimization yields an optimal input sequence and the first elem ent of the optimal input sequence is applied.

In the proposed formulation of the RHC for the formation control,each nanosatellite formulates an optimization problemover a finite time horizon.In the optimization formulation,each nanosatellite takes its current position and the last known positions of its neighboring nanosatellites as the initial condition.During the optimization process each nanosatellite optimizes independently its predicted trajectory to follow as close as possible the set ofnsamples of its future desired trajectory.Next,the cost function for the RHC is presented.

6.2 Cost function formulation

Given the current state of nanosatellitei,i.e.qi(0),the RHC is to compute the optimal input sequence which minimizes the following cost function:

wherekis the time step,ρi(0)refers to the current state ofnanosatellitei,W i(0)is the currentinput,nis the horizon length,and Φ(·)andL(·)are defined as follow s:such that R is the weight matrix.

To obtain the desired performance,the cost function(52)is solved subject to the following constraint.

6.2.1Kinem atic model equality constraint

Each nanosatellite must satisfy at all sampling times the discrete kinematic model presented in(51).The idea behind such a requirement is to obtain optimized trajectories(over each finite time horizon)that nanosatellites would be able to follow.This requirem ent is incorporated into the finite time horizon optimization problem as equality constraint as follows:

6.3 Augm ented cost function

The constraint can be incorporated into the original cost function(52)using the Lagrange multipliers condition,as follows:

where λiare the Lagrange multiplier vectors fori=1,...,N.We can define the Ham iltonian function as follows:

Therefore,(70)can be rewritten as

In order to find the desired input that minimizes the augmented cost function(57),the derivative of(57)is taken as follows:

Parameters λkfork=0,...,ncan be computed by the following recursive law s:

Therefore,(58)is simplified as follow s:

Using the gradient method,the following iterative form ula can be utilized to obtain the optimal input sequence:

w here Δiis the step length.Now,the following theorem can be expressed.

Theorem 3Consider a group of nanosatellites represented by the relative motion equations(27),which are following the reference nanosatellite represented by(5).Suppose the communication topology among the nanosatellites is represented by a directed switching graph.Under assumptions 1 and 2,the distributed optimal scheme based on the RHC approach obtained by(65)and the estimator(20),result in the convergence of Δ?ito zero,for alli=1,...,N,aswhere

ProofIn order to show the convergence of signals,we utilize the Lyapunov-based contraction analysis[34,35].Hence,w e consider the following Lyapunovlike candidate:

Taking the variation ofVT,and simplifying it,yields

Since the initial condition is fixed,its variation is zero.In addition,substituting the incremental input defined by(65),yields

Invoking Remark 2,we can guarantee that both the region error Δξiand the interaction force Δψigo to zero for each nanosatellite.

6.4 Mod i fied augmented cost function

The RHC computes the augmented cost function(70)fork=1 tok=n-1.It is known that increasingngives better performance for future decision;in other words,the current control input is effected more by future control action.This is advantageous as a smoother input can be obtained.However,a disadvantage of the increasingnis that the com putationalcostincreases too.In real-time applications the com putational cost is the major problem that must be tackled.Hence,we introduce a modified augmented cost function such that for the large values ofnthe com putational cost remains the same.In the modified cost function the cost function is computed for somekwhich is expressed as follows:

where step denotes the step size of the horizon length.

6.5 Algorithm for energy-efficient formation flying of a group of nanosatelites

We propose Algorithm 1 to solve the constraint nonlinear optimization problem presented in the previous section.

Algorithm 1(Algorithm for energy-efficient formation flying of a group of nanosatelites)

Algorithm 1 produces the desired input sequence at each time instance fora future finite time horizon,which the first elem ent of each sequence is utilized as the input.

7 Sim u lations

In this section,we present simulation results to illustrate the performance of the proposed energy-efficient method for the region-based nanosatellite formation flying.For the numerical Simulations,w e used MATLAB software.In the simulation,we consider 20 nanosatellites consisted of a leader and followers.The m ass of each nanosatellite was set toM i=5 kg,and the specifications of their propulsion system s are as follows:power 10W,specific impulseIspi=1500s,and efficiency ηi=70%.The orbital elements of the leader was considered as sem i-major axisa=6945.034km,eccentricitye=0.0002,orbit inclinationio=0.5236 rad,argument of Perigee ωp=0.3491 rad,the right ascension of the ascending node Ω=0.8727 rad,and true anomaly ν=-0.3491 rad.Hence,the orbital angular momentum and the initial orbital radius is obtained as follow s:

where gravitational parameter of earth isμ =398600km3/s2.The number of orbits was set to 1,yields to mean motionn0=0.0011 rad/s,semi-latusrectusp0=6945km,periodT0=5760s,and total timetf=5760s.The second zonal harmonic coefficient of Earth was considered asJ2=1082.626×10-6,yieldskJ2=2.6332×1010.Any two nanosatellites with the probability of 0.7 were able to directly communicate with each other.The desired region was considered as sphere with the radius of 10km and the minim um distance between the nanosatellites was set to 6 km.The gains w ere set aska=5,kij=10,κ=ε=2.2,kl=1,and Δi=0.001.

The estimations of the leader’s orbital elements are depicted in Fig.2.The motion of nanosatellites is depicted in Fig.3.Position of nanosatellites at different time instances is depicted in Fig.4.The powerconsum ed by each satellite is depicted in Fig.5.

Fig.2 The leader’s orbital elements(depicted by a thick dashed line)and the estimations of the leader’s orbital elements(depicted by thin solid lines).(a)The leader’s orbital radius and its estimation by followers.(b)The leader’s radial velocity and its estimation by followers.(c)The leader’s orbital angular momentum and its estimation by followers.(d)The leader’s right ascension of the ascending node and its estimation by followers.(e)The leader’s inclination and its estimation by followers.(f)The leader’s argument of latitude and its estimation by followers.

Fig.3 The motion of nanosatellites

Fig.4 Position of nanosatellites at different time instances.In this scenario 19 nanosatellites were assigned to follow the reference nanosatellite within the desired spherical region.(a)Initial position of nanosatellites.The desired spherical region is shown in transparent surface.(b)Position of nanosatellites in halfof the orbit.It is illustrated that the nanosatellites spread evenly within the desired region while following the reference nanosatellines.(c)Final position of nanosatellites.It can be observed that due to J2 effect on the motion of reference satellite the final position differs with the initial position.

Fig.5 Power consumption for all nanosatellites during formation flying based on the proposed method.

For comparison purposes,we applied the region formation flying controller based on the sliding m ode approach presented in Section 5.The initial conditions and the control gains were considered the same.The power consumption by each nanosatellite is depicted in Fig.6.For better comparison,the average power consumption for both methods are depicted in Fig.7.It show s that the proposed energy-efficient method presented in Section6,improved the power consumption for each nanosatellite by an average of 27.3%with respect to the region formation flying controller based on the sliding m ode approach.

Fig.6 Power consumption for all nanosatellites during formation flying based on the sliding mode controller.

Fig.7 Average power consumption of all nanosatellites for both methods.

8 Conclusions

In this paper,we have developed an energy-efficient control methodology to achieve distributed formation flying of nanosatellites.Com pared to existing methods that assume the availability of the leader’s position to all follower satellites,here a distributed estimator has been presented so that follower nanosatellitesestim ate the position of leader in its orbital slot.Using the estimated position of leader nanosatellite,an optimal region following formation method based on the RHC has been proposed.Subsequently,an algorithm has been presented to solve the proposed energy-efficient formation flying of a group of nanosatellites.The simulation result has been presented to illustrate the improvement in the power consumption of each nanosatellite with respect to the region formation flying controller based on the sliding mode approach.

Appendix

Lemm a 2[36]Consider a digraph G with the corresponding laplacian matrix L,which has a directed spanning tree;and a diagonal matrix B with nonnegative diagonal entries.The matrix H=L+B is positive stable,if and only if,at least one of the diagonal entries of B associated with the root node of the directed spanning tree,is positive.

References

[1]D.Selva,D.Krejci.A survey and assessment of the capabilities of Cubesats for Earth observation.Acta Astronautica,2012,74:50-68.

[2]H.Helvajian,S.W.Janson(eds.).Sm all Satellites:Past,Present,and Future.El Segundo,California:The Aerospace Press,2008.

[3]M.Swartwout.The first one hundred CubeSats:A statistical look.Journal of Sm all Satellites,2013,2(2):213-233.

[4]E.Gill,P.Sundaram oorthy,J.Bouwm eester,et al.formation flying within a constellation of nano-satellites:TheQB50m ission.Acta Astronautica,2013,82(1):110-117.

[5]C.M.Marrese-Reading,J.Mueller,W.C.West.Microfluidic Electrospray Thruster.US Patent#8,850,792,2014.

[6]C.I.Underwood,G.Richardson,J.Savignol.In-orbit results from the SNAP-1 nanosatellite and its future potential.Philosophical Transactions of the Royal Society A,2003,361(1802):199-203.

[7]N.G.O rr,J.K.Eyer,B.P.Larouche,et al.Precision formation flight:The CanX-4 and CanX-5 dual nanosatellite m ission.Proceedings of the 21st AIAA/USU Conference on Sm all Satellites,Noordw ijk,Netherlands:European Space Agency,2008.

[8]NASA.National Space Science Data Center:http://nssdc.gsfc.nasa.gov/nm c/spacecraftDisplay.do?id=2015-025J.

[9]M.C.VanDyke,C.D.Hall.Decentralized coordinated attitude control within a formation of spacecraft.Journal of Guidance,Control,and Dynamics,2006,29(5):1101-1109.

[10]E.Jin,X.Jiang,Z.Sun.Robust decentralized attitude coordination control of spacecraft formation.System s&Control Letters,2008,57(7):567-577.

[11]A.Zou,K.Kum ar.Distributed attitude coordination control for spacecraft formation flying.IEEE Transactions on Aerospace and Electronics System s,2012,48(2):1329-1346.

[12]J.Zhou,Q.Hu,M.I.Friswell.Decentralized finite time attitude synchronization control of satellite formation flying.Journal of Guidance,Contro l,and Dynamics,2013,36(1):185-195.

[13]N.Zhou,Y.Xia,M.Wang,et al.Finite-time attitude control of mu ltiple rigid spacecraft using term inal sliding m ode.International Journal of Robust and Non linear Contro l,2015,25(12):1862-1876.

[14]D.P.Scharf,F.Y.Hadaegh,S.R.Ploen.A survey of spacecraft formation flying guidance and control-Part II:Control.Proceedings of the Am erican Control Conference,Boston:IEEE,2004:2976-2985.

[15]R.Kristiansen,P.J.Nicklasson.Spacecraft formation flying:a review and new results on state feedback control.Acta Astronautica,2009,65(11/12):1537-1552.

[16]R. Haghighi, C. K. Pang. Concurrent attitude-position control of under-actuated nanosatellites for formation flying.IEEE International Conference on Control and Autom ation,Kathm andu,Nepal:IEEE,2016:1014-1019.

[17]S.-J. Chung, U. Ahsun, J.-J.E. Slotine. Application of synchronization to formation flying spacecraft:Lagrangian approach.Journal of Guidance,Control,and Dynamics,2009,32(2):512-526.

[18]H.S.Ahn,K.L.Moore,Y.Q.Chen.Trajectory-keeping in satellite formation flying via robustperiodic learning control.International Journal of Robust&Nonlinear Control,2010,20(14):1655-1666.

[19]D.Morgan,S.-J.Chung,F.Y.Hadaegh.Model predictive control of swarm s of spacecraft using sequential convex p rogramm ing.Journal of Guidance,Control,and Dynamics,2014,37(6):1725-1740.

[20]R.Haghighi,C.C.Cheah.Multi-group coordination control for robot swarms.Automatica,2012,48(10):2526-2534.

[21]R.Haghighi,C.K.Pang.Distributed optimal formation flying control of a group of nanosatellites.IEEE International Conference on Contro land Autom ation,Kathm andu,Nepal:IEEE,2016:1020-1025.

[22]D.J.Morgan.Swarm keeping strategies for spacecraft under J2 and atmospheric drag perturbations.Journal of Guidance Control and Dynamics,2012,35(5):1492-1506.

[23]G.Xu,D.Wang.Nonlinear dynamic equations of satellite relative motion around an ob late earth.Journal of Guidance,Control,and Dynamics,2008,31(5):1521-1524.

[24]S.-J.Chung,S.Bandyopadhyay,I.Chang,et al.Phase synchronization control of com plex networks of Lagrangian system s on adaptive digraphs.Automatica,2013,49(5):1148-1161.

[25]Y.Su,J.Huang.Stability of a class of linear switching system s with applications to two consensus problem s.IEEE Transactions on Automatic Control,2012,57(6):1420-1430.

[26]W.Zhu.Consensus of m ultiagent system s with switching jointly reachable interconnection and time delays.IEEE Transactions on System s,Man,And Cybernetics-Part A:System s And Hum ans,2012,42(2):348-358.

[27]R.Haghighi,C.C.Cheah.Two-layered framework for distributed m ultiagent formation following.IEEE Transactions on Control System s Technology,2016,24(3):1134-1139.

[28]R.Haghighi,C.C.Cheah.Optical manipulation of multiple groups of m icroobjects using robotic tweezers.IEEE Transactions on Robotics,2016,32(2):275-285.

[29]H.Liu,J.F.Li.Term inal sliding m ode control for spacecraft formation flying.IEEE Transactions on Aerospace and Electronic System s,2009,45(3):835-846.

[30]J.Wang,Z.Sun.6-DOF robust adaptive term inal sliding m ode control for spacecraft formation flying.Acta Astronautica,2012,73:76-87.

[31]A.F.Filippov.Differential Equations with Discontinuous Righthand Sides.Norwell:Kluwer,1988.

[32]D.Shevitz,B.Paden.Lyapunov stability theory of nonsmooth system s.IEEE Transactions on Automatic Control,1994,39(9):1910-1914.

[33]N.Fischer,R.Kam alapurkar,W.E.Dixon.LaSalle-Yoshizawa corollaries for nonsmooth system s.IEEE Transactions on Automatic Control,2013,58(9):2333-2338.

[34]W.Lohm iller,J.-J.E.Slotine.On contraction analysis fornonlinear system s.Automatica,1998,34(6):683-696.

[35]F.Forni,R.Sepulchre.A differential Lyapunov framework for contraction analysis.IEEE Transactions on Automatic Control,2014,59(3):614-628.

[36]J.Hu,Y.Hong.Leader-follower coordination of m ulti-agent system s with coupling time delays.Physica A,2007,374(7346):853-863.