Attitude Stabilization with Magnetic Actuation

12.07.2012 by jph

• units are SI
• vectors are in satellite body coordinates unless otherwise stated
• center of mass of satellite is at (0,0,0)

Torque generated by square coil

Edges of the coils

We derive the torque for a wire loop of square shape in the y-z plane.
4 edges with the following length and direction vectors

edg=(({{0, L, 0}, {0, 0, -L}, {0, -L, 0}, {0, 0, L}}));

Midpoints of the edges of the coils

The distance from the center of the coil to the midpoint of each edge is r

pnt={(Subscript[ofs, 1]),(Subscript[ofs, 2]),(Subscript[ofs, 3])}+#&/@(({{0, 0, r}, {0, r, 0}, {0, 0, -r}, {0, -r, 0}}));

Torque excerted on satellite

The force F of a current I going through a straight wire in direction L in a magnet field B is F=I (B×L).

We assume that the coil has n windings and that a current (Subscript[I, X]) flows through the wires. Then, the torque accumulates as 
τ=n (∑_(i=1)^4)(Subscript[pnt, i])×(I (B×L)) 
τ=n (∑_(i=1)^4)(Subscript[pnt, i])×(I (B×(Subscript[egd, i])))

For a general B, we obtain

Β={(Subscript[B, 1]),(Subscript[B, 2]),(Subscript[B, 3])}
Tor=n (Subscript[Ι, X]) (∑_(i=1)^4)Cross[pnt[[i]],Cross[Β,edg[[i]]]]//Chop

Relevant constants of satellite (in SI units)

Inertia tensor of satellite (at center of mass)

SI unit of inertia entry is kg * (m^2)

INE=(({{(Subscript[i, 11]), (Subscript[i, 12]), (Subscript[i, 13])}, {(Subscript[i, 12]), (Subscript[i, 22]), (Subscript[i, 23])}, {(Subscript[i, 13]), (Subscript[i, 23]), (Subscript[i, 33])}}))/.{(Subscript[i, 11])->0.005,(Subscript[i, 22])->0.004,(Subscript[i, 33])->0.003,(Subscript[i, 12])->-0.000002,(Subscript[i, 13])->-0.000002,(Subscript[i, 23])->-0.000003};
INE//MatrixForm

The Eigenvectors and Eigenvalues of the inertia tensor determine the long term behaviour of the satellite. More specifically, around which axis the satellite will rotate, if the satellite is subject only to small random distrubances.

{val,vec}=Eigensystem[INE];
Print[#1]&/@vec;

Maximum current in the coils X, Y, Z

the currents can flow in either directions, but are each capped at a maximum amplitude:

(Subscript[IMAX, X])=0.05;
(Subscript[IMAX, Y])=0.09;
(Subscript[IMAX, Z])=0.1;

Force of current carrying straight wire in a magnetic field

The amplitude of the Earth magnetic field ranges between approximately 25000 and 65000 nT, where 1 nT = 10^-9 T.
The SI unit for the entries of B is T = N / (A m)

AMP=65000. 10^-9

Torque excerted on satellite by coils X, Y, Z

using the formulas in the first section
we simply add up:

(Subscript[Tor, X])={0,-1.08  (Subscript[B, 3]) (Subscript[Ι, X]),1.08 (Subscript[B, 2]) (Subscript[Ι, X])};
(Subscript[Tor, Y])={1.08 (Subscript[B, 3]) (Subscript[Ι, Y]),0,-1.08 (Subscript[B, 1])(Subscript[Ι, Y])};
(Subscript[Tor, Z])={-1.08 (Subscript[B, 2]) (Subscript[Ι, Z]),1.08 (Subscript[B, 1]) (Subscript[Ι, Z]),0};
TOR=(Subscript[Tor, X])+(Subscript[Tor, Y])+(Subscript[Tor, Z])

Example: if the magnetic field B is known for instance as
B=0.000065 (0.2,0.6,0.8)
and the currents in the coils at their max
then the torque would be of magnitude about

0.000065 TOR/.{(Subscript[B, 1])->-0.2,(Subscript[B, 2])->0.8,(Subscript[B, 3])->0.5,(Subscript[Ι, X])->(Subscript[IMAX, X]),(Subscript[Ι, Y])->(Subscript[IMAX, Y]),(Subscript[Ι, Z])->(Subscript[IMAX, Z])}//Simplify
%//Norm

Equations of rotational rate

Rotation of satellite under conservation of momentum

With I=INE is the inertia tensor, and ω is the rotational rate vector, and τ the torque, the differential equation is
I.(Subscript[d, t])ω=-ω×(I.ω)+τ

vω={(Subscript[ω, 1])[t],(Subscript[ω, 2])[t],(Subscript[ω, 3])[t]};
dω=(∂_t){(Subscript[ω, 1])[t],(Subscript[ω, 2])[t],(Subscript[ω, 3])[t]};
eqs=0==#1&/@(INE.dω+Cross[vω,INE.vω]-{(Subscript[τ, 1]),(Subscript[τ, 2]),(Subscript[τ, 3])})

Simulation of angular rate and satellite attitude without control torque

Demonstration that mathematica solves the differential equations accurately

Mathematica solves the non-linear (!) differential equations of angular rate ω accurately:
• in the absence of torque, the angular momentum is preserved.
• after one cycle of ω, the orientation of the satellite matches the initial orientation
The differential equation that governs the orientation R is
(Subscript[d, t])R=R.exp(Ω  dt)

f[R_,{dω1_,dω2_,dω3_}]:=R.MatrixExp[ (({{0, -dω3, dω2}, {dω3, 0, -dω1}, {-dω2, dω1, 0}}))]
T=65;
degPsec=0.0174532925;
sol={(Subscript[ω, 1]),(Subscript[ω, 2]),(Subscript[ω, 3])}/.NDSolve[Join[eqs/.{(Subscript[τ, 1])->0,(Subscript[τ, 2])->0,(Subscript[τ, 3])->0},
{(Subscript[ω, 1])[0]==10 degPsec,(Subscript[ω, 2])[0]==3 degPsec,(Subscript[ω, 3])[0]==20 degPsec}],
{(Subscript[ω, 1]),(Subscript[ω, 2]),(Subscript[ω, 3])},
{t,0,T}][[1]];
crv[t_]=#[t]&/@sol;
ParametricPlot3D[crv[t],{t,0,T}]
dt=0.1;
ORI=Fold[f,IdentityMatrix[3],Table[dt crv[k],{k,0,T,dt}]];
ORI//MatrixForm
ORL=FoldList[f,IdentityMatrix[3],Table[dt crv[k],{k,0,T,dt}]];
Dimensions[ORL]
ListPlot[{ORL[[All,1,1]],ORL[[All,2,2]],ORL[[All,3,3]]}]

Simulation of control

• every 5 sec measure magnetic field
• in every 0...5 sec block apply magnetic torquers between t=2...5 sec
• during the simulation time, the external magnetic field vector B remains constant in world coordinates, but is therefore not constant in satellite body coordinates

Bom={0,0,0};
STEP[{Β_,R0_,ω0_,τ0_,t0_,dt_,CTr_,func_}]:=Module[
{Bo=Inverse[R0].Β,R1,ω1=ω0,smp=N[dt] Range[0,5-1]/5,
mag,tor,cτr,t1=t0+dt,seg=Mod[t0,5]},
If[seg==1,Bom=Bo];
cτr=If[2≤seg,
func[Bom,ω0,t0],
{(Subscript[Ι, X])->0,(Subscript[Ι, Y])->0,(Subscript[Ι, Z])->0}
];
mag={(Subscript[B, 1])->Bo[[1]],(Subscript[B, 2])->Bo[[2]],(Subscript[B, 3])->Bo[[3]]};
tor=TOR/.cτr/.mag;
sol={(Subscript[ω, 1]),(Subscript[ω, 2]),(Subscript[ω, 3])}/.NDSolve[Join[eqs/.{(Subscript[τ, 1])->tor[[1]],(Subscript[τ, 2])->tor[[2]],(Subscript[τ, 3])->tor[[3]]},
{(Subscript[ω, 1])[0]==ω0[[1]],(Subscript[ω, 2])[0]==ω0[[2]],(Subscript[ω, 3])[0]==ω0[[3]]}],
{(Subscript[ω, 1]),(Subscript[ω, 2]),(Subscript[ω, 3])},{t,0,dt}][[1]];
crv[t_]=#[t]&/@sol;
R1=Fold[f,R0,Table[smp[[2]] crv[smp[[k]]],{k,Length[smp]}]];
{Β,Orthogonalize[R1],crv[dt],cτr,t1,dt,{(Subscript[Tor, X]),(Subscript[Tor, Y]),(Subscript[Tor, Z])}/.mag/.cτr,func}
];
DISP[sim_]:=Module[{},
Print["Evolution of angular rate, in 3d, and norm of angular rate:"];
Print[{ListPlot[Transpose[sim[[All,3]]],PlotRange->All],ListPointPlot3D[sim[[All,3]]],
ListPlot[Norm/@sim[[All,3]],PlotRange->All]}];
Print["Control parameters: (Subscript[Ι, X]) (Subscript[Ι, Y]) (Subscript[I, Z])"];
Print[{ListPlot[(Subscript[Ι, X])/.sim[[All,4]]],ListPlot[(Subscript[Ι, Y])/.sim[[All,4]]],ListPlot[(Subscript[Ι, Z])/.sim[[All,4]]]}];
Print["Orientation (3 vectors of state matrix):"];
Print[{ListPlot[Transpose[sim[[All,2,1]]]],ListPlot[Transpose[sim[[All,2,2]]]],
ListPlot[Transpose[sim[[All,2,3]]]]}];];

Initial conditions

To compare different control strategies, we define a set of initial conditions

INIT={
{AMP {0,0,1},IdentityMatrix[3],{10,-20,16}degPsec,{(Subscript[Ι, X])->0,(Subscript[Ι, Y])->0,(Subscript[Ι, Z])->0},0,1,0 IdentityMatrix[3]},
{AMP {0,0,1},IdentityMatrix[3],{16,10,-20}degPsec,{(Subscript[Ι, X])->0,(Subscript[Ι, Y])->0,(Subscript[Ι, Z])->0},0,1,0 IdentityMatrix[3]}
};

Control using simple discrete optimization

An alternative control strategy that minimizes the scalar product of torque τ and ω, manages to stabilize attitude slightly faster.

box[cx_,cy_,cz_]:={
If[cx==1,(Subscript[IMAX, X]),-(Subscript[IMAX, X])],
If[cy==1,(Subscript[IMAX, Y]),-(Subscript[IMAX, Y])],
If[cz==1,(Subscript[IMAX, Z]),-(Subscript[IMAX, Z])]
};
STEEP[Bom_,ω0_,t0_]:=Module[
{cor=Flatten[Array[box,{2,2,2}],2],
subbf={(Subscript[B, 1])->Bom[[1]],(Subscript[B, 2])->Bom[[2]],(Subscript[B, 3])->Bom[[3]]},
subst,argsbt={(Subscript[Ι, X])->0,(Subscript[Ι, Y])->0,(Subscript[Ι, Z])->0},min=0,val},
For[i=1,i≤Length[cor],i++,
subst={(Subscript[Ι, X])->cor[[i,1]],(Subscript[Ι, Y])->cor[[i,2]],(Subscript[Ι, Z])->cor[[i,3]]};
val=TOR.ω0/.subst/.subbf;
If[val<min,min=val;argsbt=subst]
];argsbt
]
DISP[NestList[STEP,Append[#1,STEEP],900]]&/@INIT;

Binary export for visualization

CppTensorWrite[file_,tensor_,type_:"Real64"]:=Module[
{ndims=Length[Dimensions[tensor]]},
BinaryWrite[file,ndims,"Integer32"];
BinaryWrite[file,#1,"Integer32"]&/@Dimensions[tensor];
BinaryWrite[file,#1,type]&/@Flatten[Transpose[tensor,Reverse[Range[ndims]]]];
Close[file];
Print[file<>" "<>ToString[Dimensions[tensor]]];
];

ECPP[res_]:=Module[{dir="d:\\cpp9\\attitude\\get\\nb\\"},
CppTensorWrite[dir<>"bfl.dat",Transpose[Append[#1,0]&/@DAT[[All,1]],{2,1}]];
CppTensorWrite[dir<>"pos.dat",Transpose[(({{#1[[1,1]], #1[[1,2]], #1[[1,3]], 0}, {#1[[2,1]], #1[[2,2]], #1[[2,3]], 0}, {#1[[3,1]], #1[[3,2]], #1[[3,3]], 0}, {0, 0, 0, 1}}))&/@DAT[[All,2]],{3,1,2}]];
CppTensorWrite[dir<>"ome.dat",Transpose[Append[#1,0]&/@DAT[[All,3]],{2,1}]];
CppTensorWrite[dir<>"ctr.dat",Transpose[{(Subscript[Ι, X]),(Subscript[Ι, Y]),(Subscript[Ι, Z])}/.DAT[[All,4]],{2,1}]];
CppTensorWrite[dir<>"trq.dat",Transpose[{Append[#1,0]&/@DAT[[All,7,1]],Append[#1,0]&/@DAT[[All,7,2]],Append[#1,0]&/@DAT[[All,7,3]]},{2,3,1}]];
]

DAT=NestList[STEP,Append[INIT[[1]],STEEP],600];
ECPP[DAT]