Method and apparatus for determining the angular momentum vector of a satellite

Abstract

The invention relates to a method for determining the angular momentum vector of a satellite, located in an external magnetic field, relative to a satellite-integral orthogonal system of coordinates. The satellite has a torque generating system to generate pulse-shaped torques around all three axes. A magnetometer measuring on three axes is used, whose measured values are fed to a computer. In the latter, at successive points in time and using the current measured values, an observer equation is integrated, with new estimated values being determined for the angular momentum vector. Using these as well as the current measured values, an inequality is continuously checked; when it is not fulfilled, the torque generating system operates briefly, so that finally the fulfillment of the inequality is forced and identity of the last estimated value that has been determined with the rotary pulse vector to be determined is achieved.

Description

BACKGROUND AND SUMMARY OF THE INVENTION

The invention relates to a method and apparatus for determining the angular momentum vector of a satellite located in an external magnetic field relative to a satellite integral orthogonal system of coordinates.

In satellites, especially earth satellites, it is frequently necessary to determine the rotational speed (rotational speed vector .omega.), or the angular momentum vector H which is related to rotational speed by the inertia matrix of the satellite. This need arises not only when the satellite has already reached its final circumferential orbit, but also when it is still on the transfer orbit. Knowing the angular momentum vector or the rotational speed is necessary, for example, to establish the initial state of the satellite before commencing position control maneuvers, and to control its movement during the performance of such a maneuver. It is also necessary in order to check the resultant final state.

Usually special sensors such as gyroscope packages and solar sensor arrays are used to determine the rotational speed. Gyroscope packages have the disadvantage that they are expensive and suffer from a certain functional risk, so that redundant design is necessary. Using solar sensors for the same purpose is disadvantageous because they cannot supply any measured values when the spacecraft is in the earth's shadow, and are also subject to limitations of the visual field.

The goal of the present invention is to provide a highly reliable arrangement for determining satellite angular momentum which minimizes equipment cost, and which in particular avoids the above-mentioned disadvantages of the previous methods.

This goal is achieved by the method and apparatus according to the invention, which uses only the measured values of a magnetometer along three axes to determine the angular momentum vector H. (It is presumed that the satellite is in an external magnetic field B. This assumption is always fulfilled in earth satellites, especially those that are in low circular orbits that are inclined relative to the plane of the Equator, and are used preferentially.) Measured values from other measuring instruments are therefore not required. In particular, the use of expensive gyroscope packages and solar sensors are eliminated.

According to the invention, the measured values from the magnetometer are fed to a computer, where they are processed to derive an estimated value H of the angular momentum vector. The process is an iterative one in which H is estimated repeatedly, with the current measured values of the magnetometer being entered each time the integration is performed, and with new estimated values for H resulting each time. The estimated values converge with the desired angular momentum vector H which is thus obtained directly as a result of repeated integration- Besides the measured values from the magnetometer, the computer also receives as parameters the known elements of the inertia matrix of the satellite, it being assumed that a .satellite-integral orthogonal system of coordinates serves as the basis, whose axes coincide with the principal inertia axes of the satellite.

The repeated integrations are initially performed without additional measures until the check of a condition (in the form of an inequality) indicates that a small change in the angular momentum of the satellite must be performed with the aid of a torque generating system of the satellite, for example by orientation control jets or other momentum generators that may be provided. In addition to the new estimated value H just determined, the current measured values of the magnetometer and the inertia matrix of the satellite enter into this condition, and the required change in angular momentum .DELTA.H is dimensioned so that it results in the inequality being satisfied.

The invention takes advantage of the fact that an external magnetic field can be used as a reference not only for orienting the satellite, but also for its rotation. However, while an exact knowledge of the magnetic field must be presumed to determine its orientation, the same is not absolutely necessary for determining the angular momentum vector of the satellite, as will be shown below. Successful implementation of the method according to the invention is likewise not prevented by the fact that the external magnetic field is constantly changing in magnitude and direction, as is the case for example in a steeply inclined circular near-earth orbit. Thus, although the direct measurements of the magnetometer can only pick up those components of the angular momentum vector H that are oriented perpendicularly to the external magnetic field, this difficulty may be overcome by a suitable choice of the observer equation to be integrated repeatedly, and of the values occurring therein.

Thus, according to the invention, in addition to the computer itself, a magnetometer and a torque generating system are used to measure the components of the external magnetic field and to generate torque pulses, which deliberately produce a small change in the angular momentum of the satellite. Constantly updated measured values are required for a certain period of time, and the torque generating system may have to be used several times.

The invention is based on the fact that the following relationship exists between an external magnetic field (magnetic field vector B) and the rotational speed vector .omega. of the satellite:

B+.omega..times.B=QB.sub.i ( 1) in which the vectors .omega. and B are expressed in relation to the satellite-integral system of coordinates, while vector B.sub.i represents the magnetic field vector in an inertial system of coordinates and Q represents the transformation matrix from an inertial to a satellite-integral system of coordinates. This equation also takes into account the fact that the external magnetic field can change locally. This is expressed in the equation as the rate of change with respect to time B.sub.i (that is, the first time derivation of B.sub.i) of the inertial magnetic field vector, and includes the transformation matrix Q that relates the two systems of coordinates.

If the satellite orbit (and hence the external magnetic field associated with it) as well as the position of the satellite are not known sufficiently accurately, the following system of equations can be prepared, disregarding the righthand side of Equation (1):

B.sub.m +.omega..times.B.sub.m =0 (2) .omega..sub.m .multidot.B.sub.m =0. (3)

Equation (2) contains as measured values the magnetic field vector B.sub.m (which is defined in the satellite-integral system of coordinates, and unavoidably includes measurement errors) as well as the rotational speed vector .omega..sub.m (which can be derived from B.sub.m). In this equation, errors in the magnetic field measurements and (as noted above) the right-hand side of Equation (1) are disregarded. Moreover, this equation does not include any component of the real rotational speed vector .omega. that may be oriented parallel to magnetic field vector B.sub.m. The latter is expressed in Equation (3). The solution to the system of equations is given by: ##EQU1## where the normalized magnetic field vector b is as follows: ##EQU2## wherein: ##EQU3##

A quadratic matrix A of rank 3, with three mutually orthogonal line vectors, is then introduced, one of which (for example the one in the third line) is identical to the normalized magnetic field vector b: ##EQU4## with the following expressions expressing orthogonality:

a.sub.1 .multidot.a.sub.2 =a.sub.1 .multidot.b=a.sub.2 .multidot.b=0(7) as well as the additional properties a.sub.1.sup.2 =a.sub.2.sup.2 =1 detA=1,A.sup.-1 =A.sup.T. (8) If .vertline.b.sub.z .vertline.=1, a.sub.1 and a.sub.2 can be chosen as follows: a.sub.1 =(1,0,0).sup.T, a.sub.2 =(0,1,0).sup.T signb.sub.z. Otherwise the following can be used for example for the two line vectors a.sub.1,2 : ##EQU5## with the above conditions again being observed.

The rotational speed vector .omega., as mentioned above, can be resolved into parts that are oriented parallel to and vertical to the magnetic field vectors B and b:

.omega.=.omega..sub.p +.omega..sub.s. (12)

Because b is a unit vector, and .omega..sub.s is identical to the only measurable component .omega..sub.m that has already been mentioned, the above can be expressed as:

.omega.=(.omega..multidot.b)b+.omega..sub.m. (13)

Using matrix A, the following is obtained:

A.omega..sub.m =A.omega.-A(.omega..multidot.b)b. (14)

Finally, with .omega.=I.sup.-1 H and observing conditions (7) while introducing the new vector parameter Y:

Y.ident.A.omega..sub.m =AI.sup.-1 H-(0,0,b.multidot..omega.).sup.T.(15)

Here I.sup.-1 is the inverse matrix of inertia matrix I of the satellite.

The equation of motion of the satellite is:

H+.omega..times.H=M , H=I.omega. (16) where H is the time derivative of the angular momentum vector H and M is the vector of the acting external torque. In theory, if M and I are both known, the angular momentum vector H can be obtained by integration of this equation. However, since these values are not completely susceptible to measurement, the following observer equation is prepared by introducing estimated values H, .omega., and M: H+.omega..times.H-I.sup.-1 A.sup.T K(Y-AI.sup.-1 H)=M. (17)

Here, in addition to the above-mentioned matrices A and the transposed matrix A.sup.T as well as I.sup.-1 and the vector value Y defined in Equation (15), there is also a starting matrix K. The latter is determined so that, after repeated iteration of the observer equation (17) using repeatedly determined estimated values H and the respective current measured values b, the difference vector

.epsilon..ident.H-H converges toward zero. Thus, the estimated value H determined with each of the repeated integration, finally matches the angular momentum vector H to be determined.

The following differential equation is obtained for vector value .epsilon. taking into account .omega.=I.sup.-1 H and Equation (15) by subtraction of (16) and (17):

.epsilon.+(I.sup.-1 H).times.H-(I.sup.-1 H).times.H+I.sup.-1 A.sup.T K�AI.sup.-1 .epsilon.-(0,0,b.multidot..omega.).sup.T !=0. (18)

The vector term (0,0,b.multidot..omega.).sup.T that contains the parallel component of .omega. in direction b can be eliminated by setting all of the coefficients in the third column in the starting matrix K equal to zero: ##EQU6##

There is a correlation between matrices A and K such that in K the coefficients are made equal to zero in precisely that column whose number coincides with the number of the line in A whose line vector is identical with the normalized magnetic field vector b. In the special case under discussion here, it is the third column of K and the third line of A that are involved. Observer equation (17), which contains as the vector equation three scalar equations for components H.sub.i (i=x,y,z) must now be numerically integrated repeatedly at short time intervals, (taking into account a starting value for H, for example H.sub.i =0), each time taking into account the latest measured value b by which matrix A is determined. This results in new estimated values H each time, which are used as starting values in the following integration. In order to assure that the difference .epsilon. converges toward zero in the foreseeable future, it is sufficient that: ##EQU7##

With a proper choice of the satellite integral system of coordinates (namely in accordance with the main inertia axes of the satellite), the inertia matrix I is a diagonal matrix, in which consequently only the diagonal coefficients I.sub.i (i=x,y,z) occur, according to ##EQU8##

Thus, in view of condition (19), it follows from equation (18) and condition (20) that:

.epsilon..sup.T .multidot..epsilon.=-u.sup.T Pu, u.ident.AI.sup.-1 .epsilon.(22) with ##EQU9## and ##EQU10##

Here, P, N, and R are quadratic matrices and the right-hand side of the first equation (22) has a quadratic form, for which condition (20) is known to be fulfilled when matrix P is positively definite or positively semi-definite. This is equivalent to saying that when matrix P is symmetrical, the following applies for its coefficients P.sub.lm (l,m=1,2,3) and for its determinants:

P.sub.11 .gtoreq.0 P.sub.11 P.sub.22 -P.sub.12.sup.2 .gtoreq.0 (26) detP.gtoreq.0 where the > sign corresponds to the requirement that matrix P must be positively definite and the equality sign corresponds to the requirement that P must be positively semi-definite. In particular, all the diagonal elements of P must be larger than (positively definite) or equal to (positively semi-definite) zero.

As indicated by Equations (23) to (25), matrix P is determined by the coefficients of starting matrix K, measured values b.sub.i (i=x,y,z) (matrices A and A.sup.T), and the coefficients I.sub.i of inertia matrix I and estimated values H.sub.i of the components of the estimated angular momentum vector H (matrix R). When all the matrices are multiplied out, it develops that in the resultant matrix P a coefficient with the form (Ib.times.b).sup.T H occurs that is independent of the coefficient k.sub.lm of matrix K and occurs at one place on the diagonal. From the requirements, it then follows that the following must be true:

(Ib.times.b).sup.T .multidot.H.gtoreq.0. (27)

A particularly skillful choice of the coefficients of matrix K without affecting limitation of generality, consists in the fact that when condition (19) is satisfied, a relationship is established between the remaining coefficients of K and coefficient n.sub.lm of matrix N such that the following is true: ##EQU11## Here, P has the particularly simple form ##EQU12##

For coefficients c.sub.1 and c.sub.2, on the basis of equation (26), c.sub.1,2 .gtoreq.0, which can always be achieved easily by appropriate setting.

The coefficients of matrix P all have the same dimension and consequently C.sub.1,2 have the dimension of IH (b has dimension 1) or I.sup.2 .omega., i.e. (Nm sec.sup.2).sup.2 /sec. Accordingly, the following is suggested as one possibility for determining coefficients c.sub.1,2 : ##EQU13##

Hence, the requirement that c.sub.1,2 be greater than 0 is met. The corresponding requirement is also met for the remaining coefficients P.sub.33 of matrix P. For this purpose, after each integration that is performed that supplies a new value for H, using the corresponding latest measured value for b, it is necessary to check whether inequality (27) is fulfilled or not. If so, no intervention with the torque generating system is required. Otherwise, however, H must be modified by an intervention such that the inequality is once again fulfilled. For this purpose, a suitably determined small angular momentum change .DELTA.H can be effected according to:

.DELTA.H=(M.sub.x .DELTA.t.sub.x, M.sub.y .DELTA.t.sub.y M.sub.z .DELTA.t.sub.z).sup.T ( 31) with torques M.sub.i (i=x,y,z) being exerted in the form of pulses over each time interval .DELTA.t.sub.i (for example by the position control jets that act in all three axes). The torques M.sub.i can be made positive or negative in each case; it is merely necessary to fulfill the following inequality as a result: (Ib.times.b).sup.T .multidot.(H+.DELTA.H).gtoreq.0. (32) After introducing vector Z Z.ident.(z.sub.x,z.sub.y,z.sub.z).sup.T =(Ib).times.b (33) the further requirement is obtained: z.sub.x M.sub.x .DELTA.t.sub.x +z.sub.y M.sub.y .DELTA.t.sub.y +z.sub.x M.sub.z .DELTA.t.sub.z.gtoreq.-( IB.times.b).sup.T .multidot.H=-z.sub.x H.sub.x -z.sub.y H.sub.y Hy-z.sub.z H.sub.z. (34)

This inequality can be fulfilled component wise, for example. If the undesirable case Z.multidot.H<0 occurs according to inequality (27), the following must apply for at least one component:

z.sub.i .multidot.H.sub.i <0. In this case, the choice can be made as follows: z.sub.i .DELTA.H.sub.i .ident.z.sub.i M.sub.i .DELTA.t.sub.i .gtoreq..vertline.z.sub.i H.sub.i .vertline. (36) or ##EQU14## Other objects, advantages and novel features of the present invention will become apparent from the following detailed description of the invention when considered in conjunction with the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWING

The single FIGURE of the drawing shows a conceptual block diagram of the arrangement for angular momentum determination according to the invention.

DETAILED DESCRIPTION OF THE DRAWING

The FIGURE shows a magnetometer 1, a computer 10, memories 2 and 20, and a system of jets 3 acting in three axes. Magnetometer 1 measures the current values of magnetic field vector B(t.sub.k) at successive time points t.sub.k (k.gtoreq.1) which are separated from one another by time or scanning intervals .DELTA..sub.t =t.sub.k-1 :

B(t.sub.k)=(B.sub.x (t.sub.k),B.sub.y (t.sub.k),B.sub.z (t.sub.k)).sup.T(38) and feeds the measured values to a memory block 22. Next the time derivative B(t.sub.k) of the magnetic field vector B(t.sub.k) is determined in a computer block 11, for example with the aid of the statement ##EQU15## where s is the Laplacian or differential operator and T is a time constant, and also by means of the definition ##EQU16## From this it follows that: C+TC=B.sub.m (41) ##EQU17## since, because of (39) and (40): C=sC=B. (43)

From Equation (42) we have the following for the numerical processing: ##EQU18##

Computer block 11 determines B(t.sub.k) using these equations, with parameters T and .DELTA.t being taken from memory 2 and the values relating to the previous point in time t.sub.k-1 being taken from memory block 22. The current value B(t.sub.k) is obtained from magnetometer 1.

The auxiliary values C(t.sub.k) are each fed to a memory block 21 and the values of C(t.sub.k-1) already stored there are removed from the memory. (The values B(t.sub.0) and C(t.sub.0)=B(t.sub.0) are entered as initial values into memory blocks 22 and 21 at the beginning, and are called up during the first computing step at time point t.sub.1 for computer 11.)

The derivative B(t.sub.k) is input to another computer block 12 where the rotational speed vector .omega..sub.m is determined in accordance with Formula (4), using the current value of B(t.sub.k).

Then, in a computer block 13, vector b(t.sub.k), matrix A(t.sub.k) , and vector Y(t.sub.k) are determined in accordance with equations (5), (6), and (10) as well as (15), and these values are entered in a memory block 23. Matrix A(t.sub.k) is fed to another computer block 14 where, after calculation of matrix R(t.sub.k) according to (25) by means of inertia matrix I and the estimated value H(t.sub.k-1) taken from a memory 25, matrix N(t.sub.k) is determined by multiplication according to (24). Then the starting matrix K(t.sub.k) is determined in such fashion that the matrix P(t.sub.k) formed from K, K.sup.T, and N according to (23) is positively definite or positively semi-definite. This can take place for example in accordance with equation (28), so that expression (29) results for P. At the beginning, initial value H(t.sub.0) must be fed to memory 25, so that a corresponding value is available to computer block 14 for the first computer step (k=1). The current value for matrix K(t.sub.k) is entered in a memory 24.

The observer equation (17) is then integrated in computer 15. For this purpose, the following form is used as a basis for numerical processing: ##EQU19##

The parameters for .DELTA.t and I are taken from memory 2, and the matrices and vectors related to the previous point in time t.sub.k-1 are taken from memories 23, 24, 25, and 29. This integration is performed for each current point in time t.sub.k and the result H(t.sub.k) is fed first to memory block 25 and then to computer block 16, which determines component P.sub.33 of matrix P in accordance with specification (29), using the current value of b(t.sub.k) and the previously given parameter value I. In a comparator 17, the scalar value P.sub.33 is then checked to determine whether P.sub.33 .gtoreq.0, as required by specification (27). If so, the value zero is fed to memory 29 as the current value of .DELTA.H(t.sub.k). Otherwise, the critical values for triggering torque pulses M.sub.i .DELTA.t.sub.i (or a corresponding change in the angular momentum .DELTA.H according to (31)), are determined in a computer 19, equation (32). The length of time intervals (pulse durations) .DELTA.t.sub.i (i=x,y,z) can then be determined in accordance with (37), with M.sub.i representing the torques that can be generated in the three axes, whose preset values can be obtained from memory block 2.

The change in angular momentum .DELTA.H(t.sub.k) to be formed in accordance with (31) is then entered in memory 29 and the values obtained for the time intervals .DELTA.t.sub.i (t.sub.k) are fed to jet system 3 to produce the corresponding torque pulses M.sub.i .DELTA.t.sub.i.

The computation cycle described above is performed at the beginning of each new scanning interval .DELTA.t.sub.k during a time that is very short relative to the duration of .DELTA.t.sub.k, (For example such computation requires about 2 msec, compared to the scanning interval of approximately 250 msec.) The computation process converges more rapidly (i.e. the number of required integration steps in computer 15 becomes smaller) the greater the extent to which requirements (32) and (27) are met. When the torque generating system is used (jet system) this therefore depends on the degree of "overfulfillment" of requirement (37), with which however an increase in energy consumption, for example jet fuel, is associated.

Although the invention has been described and illustrated in detail, it is to be clearly understood that the same is by way of illustration and example, and is not to be taken by way of limitation. The spirit and scope of the present invention are to be limited only by the terms of the appended claims.

Claims

1. Method for determining an angular momentum vector H of an object located in an externally generated magnetic field B, relative to an orthogonal system of coordinates whose axes x,y,z coincide with main inertia axes of the object, said object having a magnetometer for measuring magnetic field strength in the direction of said x, y and z axes, and a torque generating system for producing torque pulses about all three axes, said method comprising the steps of:

measuring magnetic field strength B of said externally generated magnetic field along said x, y and z axes by means of said magnetometer, wherein B=(B.sub.i).sup.T (i=x.sub.1 y,z);

entering measured values of said magnetic field strength B into a computer;

said computer calculating estimated values H of angular momentum of said object at successive points in time, for each calculation using most recently measured values of said magnetic field strength and integrating the observer equation

H=H.times.I.sup.-1 H+I.sup.-1 A.sup.T K(Y-AI.sup.-1 H)+M

with ##EQU20## wherein M is an estimated value for an angular torque vector; I is an inertia matrix of the satellite, and I.sup.-1 is its inverse matrix; A is a quadratic matrix of rank 3 with three mutually orthogonal line vectors, one of which coincides with a normalized magnetic field vector b.ident.B/.vertline.B.vertline., and A.sup.T is its transposed matrix; K is a quadratic starting matrix K=(k.sub.lm) with rank 3 (l,m=1,2,3), with coefficients k.sub.lm equal to zero in a column whose number coincides with a number of the line in matrix A in which the normalized magnetic field vector b is located, and whose other coefficients must be determined in such fashion that the matrix ##EQU21## is positively definite or semi-definite; with ##EQU22## wherein I.sub.x,y,z represent diagonal components of inertia matrix I; continuously checking calculated values H and current measured values of said normalized magnetic field vector b to determine whether an inequality

(Ib.times.b).sup.T H.gtoreq.0

is satisfied;

if said inequality is not satisfied, said torque generating system applying a torque M=(M.sub.x,M.sub.y,M.sub.z).sup.T to said object for time intervals .DELTA.t.sub.x, .DELTA.t.sub.y, .DELTA.t.sub.z, to cause a change in angular momentum .DELTA.H=(M.sub.i .DELTA.t.sub.i).sup.T whereby said inequality is satisfied; and

performing successive integrations until such integrations no longer result in a change in .vertline.H.vertline., whereby the calculated value of H is equal to an actual angular momentum vector H of said object.

2. Method according to claim 1 wherein the following relationship is created between the elements of matrix K.ident.(k.sub.lm) and those of matrix N.ident.(n.sub.lm): ##EQU23## where c.sub.1,2 >0, so that: ##EQU24##

3. Method according to claim 1 wherein time intervals .DELTA.t.sub.i according to the requirement ##EQU25## are determined, with the following still being valid:

sign M.sub.i =sign�(Ib).times.b!.sub.i =-signH.sub.i.

4. Method according to claim 1 wherein two line vectors of matrix A that do not coincide with the normalized magnetic field vector b in the case .vertline.b.sub.z .vertline.=1 become

a.sub.1 =(1,0,0).sup.T, a.sub.2 =(0,1,0).sup.T sign b.sub.z

and in the case .vertline.b.sub.z .vertline.<1 become ##EQU26## with the following for A: ##EQU27##

5. Method for determining an angular momentum vector H of an object located in an externally generated magnetic field B, relative to an orthogonal system of coordinates whose axes x,y,z coincide with main inertia axes of the object, said object having a magnetometer for measuring magnetic field strength on said x, y and z axes, and a torque generating system for producing torque pulses about all three axes, said method comprising the steps of:

measuring magnetic field strength B of said externally generated magnetic field along said x, y and z axes by means of said magnetometers, wherein B=(B.sub.i).sup.T (i=x.sub.1 y,z);

calculating estimated values H of angular momentum of said object at successive points in time using measured values of magnetic field strength according to a predetermined formula using said measured values of magnetic field strength;

after each calculation of H, using a most recently calculated value of H, current measured values of said magnetic field strength and inertia parameters of said object to determine whether a predetermined convergence criterion is satisfied;

if said predetermined convergence criterion is not satisfied, said torque generating system applying incremental torque impulses to said object on said x, y and z axes, said impulses having a magnitude and duration such that said convergence criterion is satisfied; and

repeating said process until successive calculations result in no change in said estimated value of angular momentum, whereby said estimated value of angular momentum equals actual angular momentum of said object.

6. Apparatus for determining an angular momentum vector H of a satellite located in an externally generated magnetic field, relative to an orthogonal system of coordinates whose axes x, y and z coincide with main inertia axes of the object, said apparatus comprising:

a magnetometer for measuring magnetic field strength of said externally generated magnetic field along said x, y and z axes;

a computer coupled to receive measured values of magnetic field strength from said magnetometer, for calculating estimated values H of angular momentum of said satellite at successive points in time according to a predetermined formula using said measured values of magnetic field strength;

means responsive to each calculation of an estimated value of angular momentum, for determining whether a predetermined convergence criterion is satisfied, based on a most recently calculated value of angular momentum, current measured values of said magnetic field strength and inertia parameters of said satellite;

means operative when said convergence criterion is not satisfied for applying incremental torque impulses to said satellite on said x, y and z axes, said impulses having a magnitude and duration such that said convergence criterion is satisfied; and

means for reading said estimated value of angular momentum as an actual value of angular momentum of said satellite when successive calculations of said estimated angular momentum H produce no change therein.