WO2016111317A1 - 軌道制御装置および衛星 - Google Patents

軌道制御装置および衛星 Download PDF

Info

Publication number
WO2016111317A1
WO2016111317A1 PCT/JP2016/050262 JP2016050262W WO2016111317A1 WO 2016111317 A1 WO2016111317 A1 WO 2016111317A1 JP 2016050262 W JP2016050262 W JP 2016050262W WO 2016111317 A1 WO2016111317 A1 WO 2016111317A1
Authority
WO
WIPO (PCT)
Prior art keywords
average
thruster
control
injection
satellite
Prior art date
Application number
PCT/JP2016/050262
Other languages
English (en)
French (fr)
Inventor
憲司 北村
克彦 山田
岳也 島
博 末延
Original Assignee
三菱電機株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 三菱電機株式会社 filed Critical 三菱電機株式会社
Priority to EP16735056.0A priority Critical patent/EP3243756B1/en
Priority to US15/540,726 priority patent/US10752384B2/en
Priority to JP2016568739A priority patent/JP6271043B2/ja
Publication of WO2016111317A1 publication Critical patent/WO2016111317A1/ja

Links

Images

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242—Orbits and trajectories
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/244—Spacecraft control systems
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/26—Guiding or controlling apparatus, e.g. for attitude control using jets
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
    • G05D1/08—Control of attitude, i.e. control of roll, pitch, or yaw
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/244—Spacecraft control systems
    • B64G1/245—Attitude control algorithms for spacecraft attitude control

Definitions

  • the present invention relates to an orbit control device for maintaining the orbit of a satellite and a satellite equipped with the orbit control device.
  • a satellite or spacecraft flying in a circular orbit over an altitude of about 36000km above the equator appears to be roughly stationary when viewed from the ground observers because the period of orbital motion and the rotation period of the Earth are almost the same. .
  • Such satellites are called geostationary satellites.
  • geostationary satellites are subject to various perturbations such as the tidal force of the earth and the sun, the perturbation due to the non-uniformity of the earth's gravitational potential, and the solar radiation pressure. Go. For this reason, it is necessary for the satellite to perform thruster injection for correcting variations in latitude and longitude.
  • Thruster injection that corrects longitude fluctuations is generally called east-west control, and control that corrects latitude fluctuations is generally called north-south control.
  • east-west control and north-south control to keep the satellite latitude and longitude within a desired range is generally called orbit maintenance.
  • Non-patent document 1 is an example of a technique for realizing orbit maintenance of a geostationary satellite.
  • orbit maintenance control is performed using a nonlinear optimal control method for a geostationary satellite in which electric propulsion devices are arranged in a petal shape.
  • four thrusters are jetted simultaneously, and the trajectory holding accuracy is about 0.005 °.
  • Non-Patent Document 1 the injection amount of the electric propulsion device is calculated by a non-linear programming method in order to realize highly accurate track maintenance with low propellant consumption. Therefore, it requires a large calculation cost and is not suitable for calculating the injection law.
  • An object of the present invention is to calculate a thruster injection law for performing highly accurate orbit maintenance while suppressing propellant consumption at a low calculation cost.
  • An orbit control device is a satellite orbit control device including four thrusters arranged on a satellite in directions different from each other in directions away from the center of mass of the satellite. From the orbit determination unit that determines the average orbital element and the time change rate of the average orbital element, the target value setting unit that sets the target value of the average orbital element, the average orbital element, the time change rate and the target value of the average orbital element, The control amount calculation unit that calculates the control amount of the average orbital element and the movement of the satellite are expressed by the orbital element, taking into account the combination of the thruster placement angle and the out-of-orbital motion and in-orbital motion due to multiple thruster injection amounts.
  • Control amount calculation unit by combining the thruster injection that mainly controls the out-of-track surface direction and the thruster injection that mainly controls the in-plane direction.
  • a distribution unit that calculates the injection timing and injection amount of the thruster for realizing the calculated control amount of the average orbital element, and a thruster control unit that controls the thruster according to the injection timing and injection amount calculated by the distribution unit; Prepare.
  • FIG. 1 It is a block diagram which shows the structure of the orbit control apparatus of the satellite which concerns on Embodiment 1 of this invention. It is a figure which shows arrangement
  • 3 is a flowchart showing an operation of a control amount calculation unit according to the first embodiment.
  • 3 is a flowchart showing an operation of a distribution unit according to the first embodiment.
  • 10 is a flowchart illustrating an operation of a distribution unit according to the second embodiment.
  • 10 is a flowchart illustrating an operation of a distribution unit according to the third embodiment. It is a figure which shows the example of the hardware constitutions of this invention. It is a figure which shows the example of the hardware constitutions of this invention.
  • the injection amount of the electric propulsion device is determined by nonlinear optimization calculation in order to realize highly accurate orbit maintenance, which requires a large calculation cost, and the injection amount is calculated by a satellite-borne computer. Not suitable.
  • the present invention provides an orbit control device and a satellite that suppress power consumption or the amount of propellant and that maintain the orbit of the satellite with high accuracy.
  • FIG. 1 is a block diagram showing the configuration of a satellite orbit control apparatus according to Embodiment 1 of the present invention.
  • the orbit control device 9 is mounted on the satellite 7.
  • a geostationary satellite is assumed as the satellite 7.
  • the satellite 7 includes four thrusters 91, 92, 93, and 94 and an orbit control device 9 that controls the thrusters.
  • the trajectory control device 9 includes a target value setting unit 1, a control amount calculation unit 2, a distribution unit 5, a trajectory determination unit 8, and a thruster control unit 6.
  • the control amount calculation unit 2 includes a feedforward control amount calculation unit 3 and a feedback control amount calculation unit 4.
  • the orbit control device 9 determines the injection timing and injection amount of each of the thrusters 91, 92, 93, 94 in order to maintain the orbit of the satellite 7, and the thrusters 91, 92, 93, 94 with the determined injection timing and injection amount. Control to spray.
  • the orbit determination unit 8 includes, for example, GPS information, satellite range and range rate information with respect to the ground station, satellite azimuth / elevation angle information obtained by observation with a ground optical camera, and the like.
  • the average orbital element of the satellite 7 and the time change rate of the average orbital element are determined.
  • the average orbit element is an average orbit element obtained by removing a periodically varying component from the orbit element of the satellite.
  • the target value setting unit 1 sets a target value of an average orbit element suitable for realizing the orbit maintenance of the satellite 7.
  • the target value of the average orbital element is a target value of the average longitude longitude, average eccentricity vector, average inclination angle vector, or the like.
  • the control amount calculation unit 2 calculates the average eccentricity of the satellite 7 based on the average orbit element calculated by the orbit determination unit 8, the time change rate of the average orbit element, and the target value of the average orbit element set by the target value setting unit 1.
  • the feedforward control amount and the feedback control amount of the rate vector and the average inclination angle vector are calculated.
  • the distribution unit 5 includes the control amount of the average inclination angle vector calculated by the control amount calculation unit 2, the target value of the average direct point longitude set by the target value setting unit 1, the average trajectory element calculated by the trajectory determination unit 8, and The time change rate of the average orbital element is received, and the control amount of the average eccentricity vector and the average inclination angle vector and the control amount of the average direct point longitude for maintaining the average direct point longitude in the vicinity of the target value are calculated. Then, the injection timing and the injection amount of the thruster for realizing the calculated control amount are calculated.
  • the thruster control unit 6 causes the thrusters 91, 92, 93, and 94 to perform injection according to the injection timing and injection amount of the thruster calculated by the distribution unit 5.
  • thrusters 91, 92, 93, 94 are arranged in a petal shape on the anti-earth surface of the satellite 7.
  • the thrusters 91, 92, 93, 94 are ejected in different directions on a straight line passing through the center of mass of the satellite 7 in a direction away from the central celestial body (Earth) as seen from the center of mass of the satellite 7.
  • FIG. 2 shows the arrangement of the thrusters.
  • M.M. Center of Mass
  • central celestial typically a satellite center of mass C. from the center of the earth M.M.
  • the Z B axis coincides with the velocity direction of the satellite, and the Y B axis forms a right-handed coordinate system with respect to the X B axis and the Z B axis.
  • the thrusters 91, 92, 93, and 94 are directed in the northwest direction, the northeast direction, the southwest direction, and the southeast direction, respectively.
  • the thrusters 91, 92, 93, and 94 are also referred to as NW thrusters, NE thrusters, SW thrusters, and SE thrusters, respectively.
  • the angle ⁇ 2 is an angle thrusters 91, 92, 93 and 94 makes with the Y B axis.
  • the angle phi the line obtained by projecting the line connecting the center of mass of the thruster and satellite 7 thrusters 91, 92, 93 and 94 in the X B Z B plane, an angle formed between the X B axis.
  • FIG. 7 is a diagram illustrating an example of a hardware configuration.
  • the satellite 7 includes a receiving unit 111 that receives a GPS signal, a range and a range rate signal, a trajectory control device 9 that receives a signal from the receiving unit 111 and determines a trajectory to control a thruster, and a control target. Thrusters 91, 92, 93 and 94 are provided.
  • the trajectory control device 9 controls the thruster based on the calculation result of the distribution circuit 5 and the processing circuit 101 that realizes the functions of the target value setting unit 1, the control amount calculation unit 2, the distribution unit 5, and the trajectory determination unit 8.
  • a thruster control unit 6 is provided. Even if the processing circuit is dedicated hardware, a CPU that executes a program stored in a memory (Central Processing Unit, central processing unit, processing unit, arithmetic unit, microprocessor, microcomputer, processor, DSP) It may be.
  • a memory Central Processing Unit, central processing unit, processing unit, arithmetic unit, microprocessor, micro
  • FIG. 8 is a diagram illustrating another example of the hardware configuration of the present embodiment.
  • the satellite 7 includes a receiving unit 111 that receives a GPS signal, a range and a range rate signal, a trajectory control device 9 that receives a signal from the receiving unit 111 and determines a trajectory to control a thruster, and a control target.
  • the thrusters 91, 92, 93, and 94 may be provided.
  • the processing circuit 101 in FIG. 7 corresponds to a processor (CPU) 102 and a memory 103.
  • the processing circuit 101 corresponds to, for example, a single circuit, a composite circuit, a programmed processor, a parallel programmed processor, an ASIC, an FPGA, or a combination of these. To do.
  • the functions of the target value setting unit 1, the control amount calculation unit 2, the distribution unit 5, and the trajectory determination unit 8 may be realized by a processing circuit, or the functions of the units may be realized by a processing circuit. .
  • the processing circuit is the processor (CPU) 102
  • the functions of the target value setting unit 1, the control amount calculation unit 2, the distribution unit 5, and the trajectory determination unit 8 are realized by software, firmware, or a combination of software and firmware.
  • Software and firmware are described as programs and stored in a memory.
  • the processing circuit reads out and executes the program stored in the memory, thereby realizing the function of each unit.
  • the trajectory control device 9 includes a memory for storing a program that, when executed by the processing circuit 101, results in the target setting step, the control amount calculation step, and the distribution step being executed as a result. These programs can be said to cause the computer to execute the procedures and methods of the target value setting unit 1, the control amount calculation unit 2, the distribution unit 5, and the trajectory determination unit 8.
  • the memory corresponds to, for example, a nonvolatile or volatile semiconductor memory such as RAM, ROM, flash memory, EPROM, or EEPROM, a magnetic disk, a flexible disk, an optical disk, a compact disk, a mini disk, a DVD, or the like. To do.
  • a nonvolatile or volatile semiconductor memory such as RAM, ROM, flash memory, EPROM, or EEPROM, a magnetic disk, a flexible disk, an optical disk, a compact disk, a mini disk, a DVD, or the like.
  • the memory includes an average orbital element determined by the orbit determination unit 8, a time change rate of the average orbital element, a target value of the average orbital element set by the target value setting unit 1, and a control amount calculation unit 2
  • the feedforward control amount and the feedback control amount of the average trajectory element calculated in step (1), the injection amount and injection timing of each thruster calculated by the distribution unit 5 are stored.
  • the information stored in the above-mentioned memory is appropriately read out, used for processing, and corrected.
  • a part is implement
  • a part is implement
  • the function of the trajectory determination unit 8 is realized by a processing circuit as dedicated hardware
  • the function of the bag distribution unit 5 is realized by the processing circuit reading and executing a program stored in the memory. Is possible.
  • the processing circuit 101 corresponds to the processor 102 in FIG. 8
  • the memory corresponds to the memory 103 in FIG.
  • the processing circuit can realize the above functions by hardware, software, firmware, or a combination thereof.
  • the operation of each part of the trajectory control device 9 will be described.
  • the target value setting unit 1 sets a target value for the average orbital element of the satellite 7.
  • the average average point longitude ⁇ M defined by ⁇ g is a time average of the average point longitude ⁇ M.
  • the average trajectory element may be used as the average trajectory element.
  • eccentricity vector draws a circular trajectory of radius e n in one year due to the influence of solar radiation pressure.
  • the size of the radius e n is the mass of the satellite 7, determined by the optical constants of the effective area, and a satellite surface.
  • the control amount of the average eccentricity vector in the control amount calculation unit 2 can be suppressed.
  • the control amount calculation unit 2 calculates the control amount of the average orbital element from the average orbital element, the time change rate of the average orbital element, and the target value.
  • the control amount calculation unit 2 includes a feedforward control amount calculation unit 3 and a feedback control amount calculation unit 4.
  • the feedforward control amount calculation unit 3 calculates the feedforward control amount of the average eccentricity vector and the average inclination angle vector of the satellite 7 based on the time change rate of the average orbit element calculated by the orbit determination unit 8.
  • the average eccentricity vector and average inclination angle vector of the satellite 7 are calculated based on the average orbit element calculated by the orbit determination unit 8 and the target value of the average orbit element set by the target value setting unit 1. Calculate the feedback control amount.
  • operations of the feedforward control amount calculation unit 3 and the feedback control amount calculation unit 4 will be described.
  • the feedforward control amount calculation unit 3 calculates the feedforward control amount of the average eccentricity vector and the average inclination angle vector of the satellite 7 based on the time change rate of the average orbit element calculated by the orbit determination unit 8.
  • the time rate of change at time t 0 of the average eccentricity vector e as the 'x0, e' y0 define the time rate of change at time t 0 of the average inclination angle vector i and 'x0, i' y0.
  • changes ⁇ e and ⁇ i per cycle of the average eccentricity vector and average inclination angle vector are expressed by the following equations (1) and (2).
  • the subscript “ T ” on the right shoulder represents transposition.
  • the period T of the orbital movement of the satellite 7 is equal to the rotation period of the earth.
  • a constant multiple of ⁇ e and ⁇ i is used as the feedforward control amount of the average eccentricity vector and the average inclination angle vector.
  • ⁇ e FF ⁇ (3/4) ⁇ e (3)
  • ⁇ i FF ⁇ (3/4) ⁇ i (4)
  • the average eccentricity vector and average inclination angle vector of the satellite 7 are calculated based on the average orbit element calculated by the orbit determination unit 8 and the target value of the average orbit element set by the target value setting unit 1. Calculate the feedback control amount.
  • De e 0 -e ref
  • Di i 0 -i ref.
  • PID control PID control as in the following formulas (5) and (6).
  • k ep , k ed, and k ed are a proportional gain, a differential gain, and an integral gain, respectively, for calculating the feedback control amount of the average eccentricity vector.
  • k ip , k id, and k ii are respectively a proportional gain, a differential gain, and an integral gain for calculating the feedback control amount of the average inclination angle vector.
  • the control amount output by the control amount calculation unit 2 is the sum of the feedforward control amount calculated by the feedforward control amount calculation unit 3 and the feedback control amount calculated by the feedback control amount calculation unit 4. If the control amounts of the average eccentricity vector and the average inclination angle vector output by the control amount calculation unit 2 are defined as ⁇ e and ⁇ i, respectively, ⁇ e and ⁇ i are expressed as the following equations (7) and (8).
  • Control of the average eccentricity vector is control in the in-orbit plane direction of the satellite 7, and control of the average inclination angle vector is control in the out-of-orbit plane direction of the satellite 7.
  • the control amount calculation unit 2 calculates the control amounts of the average eccentricity vector and the average inclination angle vector according to the calculations of the equations (1) to (8).
  • the control amount calculation unit 2 realizes highly accurate trajectory control by the feedforward control amount, and enhances the stability and robustness of the control system by feedback control. The operation of the control amount calculation unit 2 has been described above, but will be described below with reference to a flowchart.
  • FIG. 3 is a flowchart showing the operation of the control amount calculation unit 2 according to the first embodiment.
  • the control amount calculation unit 2 includes a feedforward control amount calculation unit 3 and a feedback control amount calculation unit 4. This will be described below with reference to a flowchart.
  • the feedforward control amount calculation unit 3 receives the time change rate of the average trajectory element from the trajectory determination unit 8 (step S250).
  • step S251 calculates the amount of change ⁇ e, ⁇ i per orbit of the average eccentricity vector and average inclination angle vector (Sstep 251).
  • step S252 values obtained by multiplying ⁇ e and ⁇ i by a constant are set as feedforward control amounts ⁇ eFF and ⁇ iFF (step S252).
  • the feedforward control amount calculation unit 3 executes the above steps S250 to S252.
  • the feedforward control amount calculation unit 3 may be processed by the processing circuit 101 or the processor 102 or may be realized by an independent processing circuit. Further, the feedforward control amount calculation unit 3 reads the time change rate of the average orbital element from the memory, and stores the calculated changes ⁇ e and ⁇ i and the feedforward control amounts ⁇ eFF and ⁇ iFF in the memory.
  • step S253 the feedback control amount calculation unit 4 receives the average trajectory element from the trajectory determination unit 8, and receives the target value of the average trajectory element from the target value setting unit 1 (step S253).
  • step S254 calculates differences De and Di between the trajectory determination value of the average trajectory element and the target value (step S254).
  • step S255 feedback control amounts ⁇ eFB and ⁇ iFB for the average eccentricity vector and the average inclination angle vector are calculated (step S255).
  • the feedback control amount calculation unit 4 executes the above steps S253 to S255.
  • the feedback control amount calculation unit 4 may be processed by the processing circuit 101 or the processor 102, or may be realized by an independent processing circuit. Further, the feedback control amount calculation unit 4 reads out the average trajectory element, the target value of the average trajectory element, and the like from the memory, and stores the calculated differences De and Di and the feedback control amounts ⁇ eFB and ⁇ iFB in the memory.
  • step S256 the feedforward control amount and the feedback control amount are finally added to calculate the control amount of the average eccentricity vector and the average inclination angle vector (step S256). Further, in step S257, the calculated control amount of the average trajectory element is transferred to the distribution unit 5 (step S257).
  • the above-described steps S256 and S257 are executed by the control amount calculation unit 2.
  • the control amount calculation unit 2 stores the control amounts of the average eccentricity vector and the average inclination angle vector in the memory.
  • the control amount of the average eccentricity vector and the average inclination angle vector calculated by the control amount calculation unit 2 the target value of the average direct point longitude set by the target value setting unit 1, and the trajectory determination
  • the average trajectory element calculated by the unit 8 and the time change rate of the average trajectory element are received.
  • the thruster injection timing and the injection amount for realizing the control amount are calculated.
  • of the control amount ⁇ i of the average inclination angle vector at the time t 0 calculated by the control amount calculation unit 2 and the deflection angle ⁇ . That is, ⁇ i is expressed as the following formula (9). ⁇ i [ ⁇ i cos ⁇ ⁇ i sin ⁇ ] T (9) However, ⁇ i > 0 is satisfied.
  • the distribution unit 5 determines the timing of the first north-south control (out-orbit control).
  • control for correcting fluctuations in longitude is generally called east-west control
  • control for correcting fluctuations in latitude is generally called north-south control.
  • An equation related to the tilt angle vector indicates that the optimum timing for the north-south control is when the average longitude of the satellite 7 coincides with ⁇ or ⁇ + ⁇ .
  • ⁇ is the circumference ratio. Therefore, when the average longitude lambda 0 of the satellite 7 at time t 0 coincides with beta or beta + [pi, performs first north-south control at time t 0.
  • mean longitude lambda mean longitude lambda 0 and if 0 matches the beta is divided into two types in the case of matching the beta + [pi, the operation of the distributor 5.
  • the distribution unit 5 calculates the timing and the injection amount of the first and second north-south control.
  • ⁇ 1 ⁇ (10)
  • ⁇ 2 ⁇ + ⁇ (11)
  • SW thrusters and SE thrusters are used in the first north-south control
  • NW thrusters and NE thrusters are used in the second north-south control.
  • the total injection amount of SW thrusters and SE thrusters and f S during the first north-south control, the total injection amount of NW thrusters and NE thruster in the second north-south control is defined as f N.
  • F S and f N are as follows: It is determined as follows.
  • A na ⁇ i ( ⁇ e x cos ⁇ + ⁇ e y sin ⁇ ) / (sin ⁇ cos ⁇ ) (12)
  • B na ⁇ i / cos ⁇ (13)
  • the angle ⁇ is the angle of thruster 91, 92, 93 and 94 makes with the Y B axis, the angle phi, the line segment connecting the center of mass of the thruster and satellite 7 thrusters 91, 92, 93 and 94 line projected onto the X B Z B plane, an angle formed between the X B axis ( Figure 2).
  • the variables A and B have a speed dimension.
  • variable B is obtained by substituting the component of impulse in the Y B axis direction generated by the SW thruster and SE thruster injection into the Gaussian planetary equation relating to the tilt angle vector.
  • variable A is obtained by substituting the components of X B-axis direction impulse caused by the injection of SW thrusters and SE thrusters Gaussian planet equations for eccentricity vector.
  • the distribution unit 5 gives f S and f N according to the following expression (14) according to the sizes of A and B.
  • the distribution unit 5 calculates a control amount for controlling the average eccentricity vector (in-plane control).
  • the average eccentricity vector is controlled by applying a force in the orbital radial direction to the satellite 7 by injecting the NW thruster and the SE thruster equally, or by injecting the NE thruster and the SW thruster equally. Achieve.
  • Average longitude when controlling the mean eccentricity vector and lambda r defines a control amount of the track radially applied to the satellite at the time and f r.
  • ⁇ r is given by the following equation (15).
  • the distribution unit 5 restricts the distribution rate ⁇ S of the SW thruster and the SE thruster in the first north-south control and the distribution rate ⁇ N of the NW thruster and the NE thruster in the second north-south control.
  • the evaluation function J in the optimization problem can be set as the following expression (17), for example.
  • the constraint condition can be set as in the following formula (18), for example.
  • Expression (18) is a constraint condition for controlling the average directly below longitude to the target value ⁇ Mref .
  • .DELTA..LAMBDA Mtotal an average mean directly below the satellite 7 caused by the mean eccentricity vector control mean longitude lambda 1 during the first north-south control and the average longitude lambda 2 2 nd north-south control in the in the average longitude lambda r This is the sum of the changes in point longitude, and can be expressed as the following equation (19).
  • Equation (19) is obtained by adding the first north-south control at the average longitude ⁇ 1, the second north-south control at the average longitude ⁇ 2 , and the average eccentricity vector control at the average longitude ⁇ r to the Gaussian planetary equation for the average nadir point longitude. obtained by substituting X B-axis direction impulse occurring.
  • ⁇ Mref is a target holding value of the longitude point directly below the satellite 7 and is set by the target value setting unit 1.
  • Equation (18) ⁇ MT is the average average point longitude of satellite 7 one cycle after the first north-south control, and can be expressed as Equation (20) below.
  • ⁇ MT ⁇ M0 + ⁇ ' M0 T + (1/2) ⁇ " M0 T 2 + ⁇ Mtotal -(3 / a) sin ⁇ sin ⁇ [(2 ⁇ S -1) f S + (1/2) (2 ⁇ N -1) f N ] T (20)
  • the first three terms are obtained by time-integrating the angular velocity and angular acceleration of the average average direct point longitude immediately before executing the first north-south control.
  • ⁇ ′ M0 , ⁇ ′ M0 , and ⁇ ′′ M0 are the average average nadir longitude, the angular velocity of the average average nadir longitude, and the angular acceleration of the average average nadir longitude of the satellite 7 at time t 0 .
  • the distribution unit 5 is given by the equation (17) while satisfying the constraint condition for controlling the average direct point longitude given by the equation (18) to the target value ⁇ Mref.
  • the distribution ratios ⁇ S and ⁇ N are determined so that the evaluation function J is minimized.
  • the distribution unit 5 determines the injection amounts of the first and second north-south control from the calculated distribution rates ⁇ S and ⁇ N as shown in the following equations (21) and (22).
  • the subscripts NW, NE, SW, and SE of the injection amount f represent an NW thruster, an NE thruster, an SW thruster, and an SE thruster, respectively.
  • the subscript 1 indicates the first time, and the subscript 2 indicates the second time.
  • Formula (21) shows the injection amount of each thruster in the first north-south control.
  • Formula (22) shows the injection amount of each thruster in the second north-south control.
  • f NW2 ⁇ N f N
  • f NE2 (1 ⁇ N ) f N
  • the distribution unit 5 passes the calculated injection timing and injection amount of the north-south control and the average eccentricity vector control to the thruster control unit 6.
  • the thruster control unit 6 controls the thrusters 91, 92, 93, and 94 according to the injection timing and the injection amount calculated by the distribution unit 5.
  • the distribution unit 5 calculates the timing of the first and second north-south control and the injection amount.
  • lambda 1 Average longitude when performing the second north-south control, when referred to as lambda 2 is the following formula (23) given by (24).
  • ⁇ 1 ⁇ + ⁇ (23)
  • ⁇ 2 ⁇ + 2 ⁇ (24)
  • NW thrusters and NE thrusters are used, and in the second north-south control, SW thrusters and SE thrusters are used.
  • the total injection amount of NW thrusters and NE thruster during the first north-south control and f N, a total injection quantity of SW thrusters and SE thrusters in the second north-south control is defined as f S.
  • f S a total injection quantity of SW thrusters and SE thrusters in the second north-south control
  • the orbital length radius and the average orbital angular velocity of the satellite 7 are represented as a and n, respectively.
  • variables A and B depending on ⁇ e and ⁇ i are defined as in the following equations (25) and (26).
  • A na ⁇ i ( ⁇ e x cos ⁇ + ⁇ e y sin ⁇ ) / (sin ⁇ cos ⁇ ) (25)
  • B na ⁇ i / cos ⁇ (26)
  • Angle ⁇ is the angle of thruster 91, 92, 93 and 94 makes with the Y B axis, the angle ⁇ for the line segment connecting the center of mass of the thruster and satellite 7 thrusters 91, 92, 93 and 94 X B Z line segment projected on B plane, an angle formed between the X B axis ( Figure 2).
  • the distribution unit 5 gives f S and f N according to the following expression (27) according to the sizes of A and B.
  • the distribution unit 5 calculates a control amount for controlling the average eccentricity vector.
  • the average eccentricity vector is controlled by applying a force in the orbital radial direction to the satellite 7 by injecting the NW star and SE thrusters equally, or by injecting the NE thruster and SW thrusters equally. Achieve.
  • Average longitude when controlling the mean eccentricity vector and lambda r defines a control amount of the track radially applied to the satellite at the time and f r.
  • ⁇ e average eccentricity vector and the unit vector [rho beta provides a mean longitude lambda r when controlling the mean eccentricity vector by the following equation (28).
  • the distribution unit 5 constrains the distribution rate ⁇ N of the NW thruster and the NE thruster in the first north-south control and the distribution rate ⁇ S of the SW thruster and the SE thruster in the second north-south control.
  • the evaluation function J in the optimization problem can be set as the following expression (30), for example.
  • the constraint condition can be set as in the following equation (31), for example.
  • ⁇ MT ⁇ Mref ⁇ (1/2) ⁇ Mtotal (31)
  • .DELTA..LAMBDA Mtotal is first north-south control in the average longitude lambda 1, 2 nd north-south control in mean longitude lambda 2, and the average mean of satellite 7 caused by the mean eccentricity vector control in the average longitude lambda r This is the sum of the change amount of the direct point longitude and can be expressed as the following equation (32).
  • ⁇ Mtotal 2 tan ⁇ cos ⁇ i +2
  • ⁇ Mref is a target holding value of the longitude directly below the satellite 7 and is set by the target value setting unit 1.
  • ⁇ MT is the average average point longitude of satellite 7 after one orbit cycle has elapsed since the first north-south control, and can be expressed as Equation (33) below.
  • ⁇ MT ⁇ M0 + ⁇ ' M0 T + (1/2) ⁇ " M0 T 2 + ⁇ Mtotal -(3 / a) sin ⁇ sin ⁇ [(2 ⁇ N -1) f N + (1/2) (2 ⁇ S -1) f S ] T (33)
  • ⁇ M0 , ⁇ ′ M0 and ⁇ ′′ M0 are the average average nadir longitude, the average average nadir velocity, and the average average nadir acceleration of the satellite 7 at time t 0, respectively.
  • the distribution unit 5 satisfies the constraint condition given by the equation (31), and the distribution rate ⁇ S so that the evaluation function J given by the equation (30) is minimized.
  • ⁇ N is determined.
  • the distribution unit 5 determines the first and second north-south injection amounts from the calculated distribution rates ⁇ S and ⁇ N as in the following equations (34) and (35).
  • the subscript of the injection amount f is the same as in the equations (21) and (22).
  • Formula (34) shows the injection amount of each thruster in the first north-south control.
  • Formula (35) shows the injection amount of each thruster in the second north-south control.
  • Equations (12) to (14), (16) and (18) to (20) are equations that take into account the coupling of the out-of-orbital motion and the in-orbital motion.
  • the distribution unit 5 solves an equation that takes into account the coupling of out-of-plane motion and in-plane motion by a plurality of thruster injection amounts, and performs thruster injection that mainly controls the out-of-plane direction and thruster injection that mainly controls the in-plane direction.
  • the injection timing and the injection amount of the thruster for realizing the control amount of the average trajectory element calculated by the control amount calculation unit 2 are calculated.
  • the distribution unit 5 passes to the thruster control unit 6 the injection timing and the injection amount of each thruster of the north-south control and the average eccentricity vector control calculated above.
  • the thruster control unit 6 controls the thrusters 91, 92, 93, and 94 according to the injection timing and the injection amount calculated by the distribution unit 5.
  • FIG. 4 is a flowchart showing the operation of the distribution unit 5 according to the first embodiment.
  • the distribution unit 5 first calculates the magnitude and declination of the control amount of the average inclination angle vector represented by the equation (9) (step S01). And the timing of north-south control is determined (step S02).
  • the injection timings ⁇ 1 and ⁇ 2 of the north-south control are calculated by the equations (10) and (11), and the total injections of each.
  • the quantities f S and f N are calculated from the equations (12) to (14) (step S03).
  • the distribution unit 5 stores in the memory the magnitude and declination of the control amount of the average inclination angle vector, the north and south control injection timings ⁇ 1 and ⁇ 2 , and the total injection amounts f S and f N.
  • the injection timing lambda r of the mean eccentricity vector control of injection quantity f r is calculated by equation (15) and (16) (step S04). Then, the distribution ratios ⁇ S and ⁇ N that minimize the evaluation function J of Expression (17) are determined under the constraint condition of Expression (18) (Step S05).
  • the distribution unit 5 calculates the injection amount of each thruster for the first time and the second time by the equations (21) and (22) from the distribution rates ⁇ S and ⁇ N (step S06). Then, the calculated injection timing and injection amount are sent to the thruster controller 6 (step S07).
  • the distribution unit 5 stores the injection timing ⁇ r and the injection amount fr of the average eccentricity vector control , the distribution rates ⁇ S and ⁇ N , and the injection timing and the injection amount in a memory.
  • the distribution unit 5 uses the equations (23) and (24). ), The injection timings ⁇ 1 and ⁇ 2 of the north-south control are calculated, and the respective total injection amounts f S and f N are calculated from the equations (25) to (27) (step S09). The distribution unit 5 stores the north-south control injection timings ⁇ 1 and ⁇ 2 and the total injection amounts f S and f N in a memory.
  • Distributor 5 is then injection timing lambda r and the injection quantity f r of the mean eccentricity vector control, is calculated by equation (28) and (29) (step S10). Then, the distribution ratios ⁇ S and ⁇ N that minimize the evaluation function J of Equation (30) are determined under the constraint condition of Equation (31) (Step S11). The distribution unit 5 calculates the injection amount of each thruster for the first time and the second time by the equations (34) and (35) from the distribution rates ⁇ S and ⁇ N (step S12). Then, the calculated injection timing and injection amount are sent to the thruster controller 6 (step S07). The distribution unit 5 stores in memory the injection timing ⁇ r and injection amount f r of the average eccentricity vector control, the distribution rates ⁇ S and ⁇ N, and the injection amounts of the first and second thrusters.
  • step S02 and step S08 when the average longitude ⁇ 0 does not coincide with the declination ⁇ or ⁇ + ⁇ (step S02; N, step S08; N), neither injection amount calculation nor injection is performed (step S13). End the process.
  • the distribution unit 5 analytically analyzes the Gaussian planetary equation including the coupling between the in-orbit surface and the out-of-orbit surface depending on the thruster arrangement angle and multiple thruster injections. By solving, it is possible to calculate the thruster injection law that suppresses the consumption propellant at a low calculation cost.
  • the satellite can be highly accurately controlled while suppressing power consumption. Can keep track.
  • the optimum north-south control injection timing can be calculated based on the control amount of the average inclination angle vector, and high-precision north-south control can be realized. Further, it becomes possible to simultaneously control the average average direct point longitude, the average inclination angle vector, and the average eccentricity vector in the north-south control.
  • the average eccentricity vector is controlled once per orbit, so that the average eccentricity vector can be held near the target value. It becomes possible to maintain longitude in the vicinity of the target value with high accuracy.
  • the interval between the north-south control and the average eccentricity vector control can be set at a quarter of the orbit, and it is possible to avoid the timing of thruster injection of the north-south control and the average eccentricity vector control from overlapping.
  • Embodiment 2 the timing at which thruster injection is possible in one round of the orbit is fixed in advance. Then, the injection amount of the thruster 2 to 4 times and the propriety of the injection are determined according to the control amount of the average eccentricity vector and the average inclination angle vector calculated by the control amount calculation unit for each timing.
  • the configuration of the second embodiment is the same as that of the first embodiment, only the operation of the distribution unit 5 in the trajectory control device 9 is different.
  • Other target value setting unit 1, control amount calculation unit 2, thruster control unit 6 and trajectory determination unit 8 are the same as those in the first embodiment.
  • the distribution unit 5 controls the average eccentricity vector and the average inclination angle vector control amount calculated by the control amount calculation unit 2, the target value of the average direct point longitude set by the target value setting unit 1,
  • the average trajectory element calculated by the trajectory determination unit 8 and the time change rate of the average trajectory element are received.
  • the injection timing and the injection amount of the thruster for realizing the control amount of the average direct point longitude for maintaining the control amount of the average eccentricity vector and the average inclination angle vector and the average direct point longitude in the vicinity of the target value are calculated.
  • the thruster injection timing is set to a maximum of four times in one cycle.
  • Each timing is expressed as ⁇ 1 , ⁇ 2 , ⁇ 3 , ⁇ 4 using the average longitude ⁇ , and is given by the following equations (48) to (51).
  • ⁇ 1 ⁇ C (48)
  • ⁇ 2 ⁇ 1 + ⁇ / 2 (49)
  • ⁇ 3 ⁇ 1 + ⁇ (50)
  • ⁇ 4 ⁇ 1 + 3 ⁇ / 2 (51)
  • ⁇ C is the first thruster injection timing given in advance.
  • thruster injection is performed with the phase being separated by ⁇ / 2 at equal intervals from the first thruster injection.
  • ⁇ C, ⁇ 1 , ⁇ 2 , ⁇ 3 , and ⁇ 4 can be selected from values other than the above-described method.
  • the number of thruster injections is four here, the number of thruster injections may be suppressed to two from the actually calculated control amount. The following describes the control force required to maintain the trajectory and how to calculate the thruster injection amount.
  • .delta.e i determines the f hi by case analysis on the magnitude relationship .delta.i i. In this case, priority is given to realizing the control amount of the average eccentricity vector. If f hi is determined from the arrangement of the thrusters, f ri is automatically determined. When the case became .delta.e i> .delta.i i would fit the control amount of the average tilt angle vector, the control amount of the average eccentricity vector becomes insufficient. Therefore, in such a case, f hi is determined based on the control amount of the average eccentricity vector.
  • a is an orbital radius
  • n is an average orbital angular velocity
  • ⁇ and ⁇ are angles given in FIG.
  • the thruster injection for giving this f ⁇ i is the distribution of the injection quantity to the SW thruster and SE thruster or the NW thruster and NE thruster when giving the out-of-plane control force f hi (in the case of electric propulsion, the injection) This is achieved by allocating time).
  • This distribution amount can be uniquely determined based on the thruster arrangement angle when f hi and f ⁇ i are given.
  • f hi varies for each control cycle, and in some cases, a sufficient total f ⁇ i cannot be obtained by four times of control. In such a case, the injection of two or one thruster is performed in the first thruster injection, and the necessary f ⁇ is additionally given.
  • the distribution unit 5 passes the above calculation result to the thruster control unit 6.
  • the basic operation in the second embodiment is as described above. However, when aiming at further saving of the propellant, the following control is performed when calculating the out-of-plane control amount f hi .
  • the average inclination angle vector error ⁇ is a sufficiently small value, for example,
  • the number of thruster injections is reduced to three or less by adjusting the control amount only to the control amount necessary for controlling the average eccentricity vector.
  • FIG. 5 is a flowchart showing the operation of the distribution unit 5 according to the second embodiment.
  • Distributor unit 5 controls the amount of average eccentricity vector and the average tilt angle vector, and stores the target value and .delta.e i of the average mean nadir longitude, reads .delta.i i from the memory, the out-of-plane control amount f hi the memory.
  • the distribution unit 5 determines the distribution of the injection amount of each thruster at each thruster injection timing (step S26), and determines the fixed injection timing and the injection amount of each thruster as a thruster control unit. 6 (step S27).
  • the distribution unit 5 stores the distribution of the injection amount of each thruster at each thruster injection timing in the memory.
  • the control amount of the average eccentricity vector and the average inclination angle vector and the control amount of the average direct point longitude for maintaining the average direct point longitude in the vicinity of the target value are realized in the second embodiment.
  • the injection timing and injection amount of the thruster to be calculated are calculated.
  • the injection interval between the thruster injections can be fixed, and the orbit maintenance error in the longitude direction caused by the fluctuation of the control time interval can be reduced.
  • Embodiment 3 In the third embodiment, two out of four thrusters are combined to perform north-south control once per orbit, and two east thrusters are controlled twice per orbit.
  • the configuration of the third embodiment is the same as that of the first embodiment, only the operation of the distribution unit 5 in the trajectory control device 9 is different.
  • Other target value setting unit 1, control amount calculation unit 2, thruster control unit 6, satellite 7 and orbit determination unit 8 are the same as those in the first embodiment. Since only the distribution unit 5 is different from the first embodiment, the distribution unit 5 will be described.
  • the injection amounts of the thrusters 91 (NW thruster), thruster 92 (NE thruster), thruster 93 (SW thruster), and thruster 94 (SE thruster) shown in FIG. 2 are respectively expressed as f NW , f NE , f SW , and f SE. .
  • the north direction injection amount f N f NW + f NE
  • f S f SW + f SE
  • f E f NE + f SE
  • f W f NW + f SW
  • control amount of the distribution unit 5 controls the control amount of the average inclination angle by north injection quantity f N or south direction injection quantity f S. Since the east direction injection amount f E and west direction injection amount f W does not affect the average angle of inclination, is to define the north direction injection amount f N or south direction injection amount f S from only control the amount of the average inclination angle it can.
  • the injection timing can be determined from the deviation angle ⁇ of the control amount ⁇ i of the average inclination angle vector.
  • This north injection quantity f N or south direction injection quantity f S is able to control the average inclination angle be used either, in the case of both calculate the total value of the final injection amount is reduced Adopt one.
  • the average eccentricity and mean average nadir longitude by north injection quantity f N or south direction injection quantity f S is affected by adding a correction amount of that amount, a new average eccentricity vector
  • This correction amount is calculated by substituting the X B axis and Z B axis direction components of the north direction injection amount f N or the south direction injection amount f S into the eccentricity vector and the Gaussian planetary equation for the mean nadir longitude. Can do.
  • eastward injection amount f E and west directions injection quantity f W is a total of two minutes is required. That is, given the eastward injection quantity f E or west injection quantity f W In certain longitude, obtaining a solution which gives an east injection quantity f E or west injection quantity f W at different longitudes.
  • partitioned section 5 determines whether to adopt either north injection quantity f N or south direction injection amount f S as follows.
  • Distributor 5 compares the total value of the East-West injection amount combined injection quantity when using the f N, a sum of the east and west injection amount combined injection quantity when using the f S, the sum The smaller value is adopted as the final injection.
  • the lesser of the way the injection amount in one injection and two injections of east-west south of f N or f S, since longitude will be away when their implementation, these The effect that the timing of injection does not overlap is also obtained.
  • the average tilt angle, the average eccentricity, and the average average directly below longitude in the series of control operations with one north-south and two east-west thrusters. can be controlled to a desired value.
  • the average inclination angle can be affected.
  • the average tilt angle is affected by the east and west thruster injections, it is not possible to easily determine the north and south thruster injection amounts.
  • both or one of the two east-west thruster injections can be planned with only one thruster. If planned in this way, a part of the average inclination angle is controlled by the thruster injection in the east-west direction, so that there is an effect of keeping the total value of the injection amount low.
  • FIG. 6 is a flowchart showing the operation of the distribution unit 5 according to the third embodiment.
  • the distribution unit 5 first receives the control amount of the average eccentricity vector and the average inclination angle vector, the target value of the average average direct point longitude, and the like (step S31). Then, the north direction injection amount f N or the south direction injection amount f S (the north-south injection amount f N or f S ) is determined from the control amount of the average inclination angle vector (step S32).
  • the distribution unit 5 reads out the control amount of the average eccentricity vector and the average inclination angle vector, the target value of the average average direct point longitude, and the like from the memory, and reads the north direction injection amount f N or the south direction injection amount f S (the north-south injection amount f N or f S ) is stored in memory.
  • step S33 the control amount of the determined north-south injection quantity f N or new average eccentricity than f S and corrects the target value of the average mean nadir longitude. Then, determined from the target value of the average mean nadir longitude control of new average eccentricity, eastward injection amount f E and west directions injection quantity f W (east-west injection quantity f E, f W) (steps S34) .
  • Distributor 5 reads the south injection quantity f N or f S from the memory, the control amount of the new average eccentricity, a target value of the average mean nadir longitude east injection amount f E and west directions injection quantity f W ( storing east-west injection quantity f E, the f W) in the memory.
  • Distributor 5 in north-south injection quantity f N or f S, employing a lesser of the total injection quantity obtained by adding the east-west injection amount (step S35). Furthermore, based on the calculation of the thruster injection amount, in the two thruster injections in the east and west, instead of the simultaneous injection of two thrusters that exerts force in the east and west directions, some or all Is recalculated so as to be replaced by the injection of only one thruster (step S36). The calculated thruster injection amount is sent to the thruster controller 6 (step S37). The distribution unit 5 stores the calculated thruster injection amount in a memory.
  • the average inclination angle is obtained by one injection in the north-south direction and two injections in the east-west direction by simple calculation on the premise of the injection combining the two thrusters. It is possible to control the average eccentricity and the average average direct point longitude to desirable values. Further, according to the configuration of the third embodiment, at least one of the two east and west thruster injections is made to be an injection of only one thruster based on a solution based on an injection that combines two thrusters. By obtaining a numerical search solution, it is possible to obtain an effect of controlling the average inclination angle, the average eccentricity, and the average average direct point longitude to desirable values and keeping the total value of the injection amount low.
  • the orbit control device 9 of the embodiment can be applied not only to geostationary satellites.
  • the orbit control device 9 can be applied to, for example, a quasi-zenith satellite.
  • the quasi-zenith satellite looks like an eight-shaped orbit long from north to south when observed from one point on the ground, but the period and altitude are the same as those of geostationary satellites, and the orbital plane is tilted with respect to the equator plane. Therefore, if the above equation is appropriately coordinate-transformed, it can be applied as it is.

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

　衛星（７）は、それぞれが衛星（７）の質量中心から遠ざかる向きに互いに異なる方向に噴射方向を向けて配置されたスラスタ（９１、９２、９３、９４）を備える。軌道決定部（８）は、衛星（７）の平均軌道要素および平均軌道要素の時間変化率を決定する。目標値設定部（１）は、平均軌道要素の目標値を設定する。制御量計算部（２）は、平均軌道要素、平均軌道要素の時間変化率および目標値から、平均軌道要素の制御量を計算する。分配部（５）は、衛星（７）の運動を軌道要素で表現し、スラスタ配置角度と複数回のスラスタ噴射量による面外運動および面内運動の連成を考慮した方程式を解き、主として軌道面外方向を制御するスラスタ噴射および主として軌道面内方向を制御するスラスタ噴射を複数回組み合わせることによって、平均軌道要素の制御量を実現させるためのスラスタ噴射タイミングおよび噴射量を計算する。

Description

軌道制御装置および衛星
　この発明は、衛星の軌道を保持するための軌道制御装置、および軌道制御装置を搭載する衛星に関する。
　赤道上の高度約３６０００km上空を円軌道で飛行する衛星もしくは宇宙飛翔体は、軌道運動の周期と地球の自転周期がほぼ一致するため、地上の観測者から見ておおよそ静止しているように見える。このような衛星は静止衛星と呼ばれている。しかしながら実際には、静止衛星には、地球・太陽の潮汐力、地球の重力ポテンシャルの非一様性による摂動力、太陽輻射圧といった種々の摂動力が働くため、徐々に緯度と経度が変動してゆく。そのため、緯度と経度の変動を補正するためのスラスタ噴射を衛星が実施する必要がある。経度の変動を補正するスラスタ噴射は一般に東西制御と呼ばれ、緯度の変動を補正する制御は一般に南北制御と呼ばれる。東西制御および南北制御を行って衛星の緯度と経度を所望の範囲内に留めることは、一般に軌道保持と呼ばれる。
　静止衛星の軌道保持を実現するための技術として、例えば非特許文献１が挙げられる。非特許文献１の技術では、電気推進装置を花弁状に配置した静止衛星を対象として、非線形最適制御の手法を用いて軌道保持制御を行う。非特許文献１では４本のスラスタを同時に噴射し、軌道保持精度は０．００５°程度である。
　しかしながら、非特許文献１に示される先行技術例では、少ない推薬消費で高精度な軌道保持を実現するために、非線形計画法によって電気推進装置の噴射量を算出している。そのため、大きな計算コストを必要としており、噴射則を算出するのに適していない。
　本発明は、推薬消費量を抑えながら、高精度な軌道保持を行うためのスラスタ噴射則を少ない計算コストで算出することを目的とする。
　この発明に係る軌道制御装置は、それぞれが衛星の質量中心から遠ざかる向きに互いに異なる方向に噴射方向を向けて衛星に配置された４本のスラスタを備える衛星の軌道制御装置であって、衛星の平均軌道要素および平均軌道要素の時間変化率を決定する軌道決定部と、平均軌道要素の目標値を設定する目標値設定部と、平均軌道要素、平均軌道要素の時間変化率および目標値から、平均軌道要素の制御量を計算する制御量計算部と、衛星の運動を軌道要素で表現し、スラスタ配置角度と複数回のスラスタ噴射量による軌道面外運動および軌道面内運動の連成を考慮した方程式を解き、主として軌道面外方向を制御するスラスタ噴射および主として軌道面内方向を制御するスラスタ噴射を複数回組み合わせることによって、制御量計算部で計算した平均軌道要素の制御量を実現させるためのスラスタの噴射タイミングおよび噴射量を計算する分配部と、分配部で計算された噴射タイミングおよび噴射量に従って、スラスタを制御するスラスタ制御部と、を備える。
　この発明によれば、消費推薬を抑えたスラスタ噴射則を少ない計算コストで算出することができる。
この発明の実施の形態１に係る衛星の軌道制御装置の構成を示すブロック図である。 スラスタの配置を示す図である。 実施の形態１に係る制御量計算部の動作を示すフローチャートである。 実施の形態１に係る分配部の動作を示すフローチャートである。 実施の形態２に係る分配部の動作を示すフローチャートである。 実施の形態３に係る分配部の動作を示すフローチャートである。 本発明のハードウェア構成の例を示す図である。 本発明のハードウェア構成の例を示す図である。
　従来は、高精度な軌道保持を実現するために複数の電気推進スラスタを同時に噴射している。電気推進スラスタは大きな電力を消費するため、複数のスラスタを同時に噴射することは実運用上の大きな制約となる。また、先行技術例では高精度な軌道保持を実現するために非線形最適化計算によって電気推進装置の噴射量を決定するため、大きな計算コストを必要とし、衛星搭載計算機で噴射量を計算するのに適していない。本発明は、消費電力または推薬量を抑制するとともに、高精度に衛星の軌道を保持する軌道制御装置および衛星を提供するものである。
　実施の形態１．
　図１はこの発明の実施の形態１に係る衛星の軌道制御装置の構成を示すブロック図である。軌道制御装置９は衛星７に搭載される。実施の形態では、衛星７として静止衛星を想定する。衛星７は、４つのスラスタ９１、９２、９３、９４とそれを制御する軌道制御装置９を備える。軌道制御装置９は、目標値設定部１、制御量計算部２、分配部５、軌道決定部８、およびスラスタ制御部６から構成される。制御量計算部２は、フィードフォワード制御量計算部３とフィードバック制御量計算部４を含む。軌道制御装置９は、衛星７の軌道を保持するためにスラスタ９１、９２、９３、９４それぞれの噴射タイミングと噴射量を決定し、決定した噴射タイミングと噴射量でスラスタ９１、９２、９３、９４に噴射させるよう制御する。
　実施の形態１の軌道制御装置９では、軌道決定部８は、例えばＧＰＳ情報、地上局に対する衛星のレンジ及びレンジレート情報、地上の光学カメラでの観測で得られる衛星の方位角・仰角情報などから、衛星７の平均軌道要素および平均軌道要素の時間変化率を決定する。ここで、平均軌道要素とは、衛星の軌道要素から周期的に変動する成分を除去して得られる、平均的な軌道要素である。目標値設定部１は、衛星７の軌道保持を実現する上で適切な平均軌道要素の目標値を設定する。平均軌道要素の目標値とは、平均直下経度や平均離心率ベクトル，平均傾斜角ベクトルの目標値などである。制御量計算部２は、軌道決定部８で計算した平均軌道要素、平均軌道要素の時間変化率、および目標値設定部１で設定した平均軌道要素の目標値を元に、衛星７の平均離心率ベクトルおよび平均傾斜角ベクトルのフィードフォワード制御量およびフィードバック制御量を計算する。分配部５は、制御量計算部２で計算した平均傾斜角ベクトルの制御量、および目標値設定部１で設定した平均直下点経度の目標値、軌道決定部８で計算した平均軌道要素、および平均軌道要素の時間変化率を受け取り、平均離心率ベクトルと平均傾斜角ベクトルの制御量および平均直下点経度を目標値近傍に維持するための平均直下点経度の制御量を計算する。そして、計算した制御量を実現するためのスラスタの噴射タイミングおよび噴射量を計算する。スラスタ制御部６は、分配部５で計算したスラスタの噴射タイミングと噴射量に従ってスラスタ９１、９２、９３、９４に噴射を実行させる。
　スラスタ９１、９２、９３、９４は、衛星７の反地球面に花弁状に４台配置される。言い換えると、スラスタ９１、９２、９３、９４は、衛星７の質量中心から見て中心天体（地球）から遠ざかる向きに、衛星７の質量中心を通る直線上で、各スラスタが互いに異なる方向に噴射方向を持つように配置される。図２は、スラスタの配置を示す。図２において、Ｃ．Ｍ．（Center of Mass）は、衛星７の質量中心を表す。図２においてＸＢ軸は、中心天体、典型的には地球の中心から衛星質量中心Ｃ．Ｍ．を指す方向に一致し、ＺＢ軸は衛星の速度方向に一致し、ＹＢ軸はＸＢ軸とＺＢ軸に対して右手座標系を構成する。図２において、スラスタ９１、９２、９３および９４は、噴射方向がそれぞれ、北西方向、北東方向、南西方向および南東方向に向いている。スラスタ９１、９２、９３および９４は、それぞれＮＷスラスタ、ＮＥスラスタ、ＳＷスラスタおよびＳＥスラスタとも呼ぶこととする。図２において角度θは、スラスタ９１、９２、９３および９４がＹＢ軸と成す角度である。角度φは、スラスタ９１、９２、９３および９４の各スラスタと衛星７の質量中心を結ぶ線分をＸＢＺＢ平面に投影した線分が、ＸＢ軸と成す角度である。
　次に、本実施の形態のハードウェア構成について説明する。図７は、ハードウェア構成の例を示す図である。図のように、衛星７は、ＧＰＳ信号やレンジ及びレンジレート信号を受信する受信部１１１、受信部１１１からの信号を受けて軌道決定してスラスタの制御を行う軌道制御装置９、並びに制御対象のスラスタ９１、９２、９３および９４を備える。また、軌道制御装置９は、目標値設定部１、制御量計算部２、分配部５および軌道決定部８の各機能を実現する処理回路１０１と、分配部５の計算結果に基づきスラスタを制御するスラスタ制御部６を備える。処理回路は、専用のハードウェアであっても、メモリに格納されるプログラムを実行するＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ、中央処理装置、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、プロセッサ、ＤＳＰともいう）であってもよい。
　また、図８は、本実施の形態のハードウェア構成の別の例を示す図である。図のように、衛星７は、GＰＳ信号やレンジ及びレンジレート信号を受信する受信部１１１、受信部１１１からの信号を受けて軌道決定してスラスタの制御を行う軌道制御装置９、並びに制御対象のスラスタ９１、９２、９３および９４を備えるように構成してもよい。図７の処理回路１０１は、プロセッサ（ＣＰＵ）１０２とメモリ１０３に当たる。
　処理回路１０１が、専用のハードウェアである場合、処理回路１０１は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ，ＦＰＧＡ、またはこれらを組み合わせたものが該当する。目標値設定部１、制御量計算部２、分配部５および軌道決定部８の各部の機能それぞれを処理回路で実現してもよいし、各部の機能をまとめて処理回路で実現してもよい。
　処理回路がプロセッサ（ＣＰＵ）１０２の場合、目標値設定部１、制御量計算部２、分配部５および軌道決定部８の機能は、ソフトウェア、ファームウェア、またはソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェアやファームウェアはプログラムとして記述され、メモリに格納される。処理回路は、メモリに記憶されたプログラムを読み出して実行することにより、各部の機能を実現する。すなわち、軌道制御装置９は、処理回路１０１により実行されるときに、目標設定ステップ、制御量計算ステップ、分配ステップが結果的に実行されることになるプログラムを格納するためのメモリを備える。また、これらのプログラムは、目標値設定部１、制御量計算部２、分配部５および軌道決定部８の手順や方法をコンピュータに実行させるものであるともいえる。ここで、メモリとは、例えば、ＲＡＭ、ＲＯＭ、フラッシュメモリー、ＥＰＲＯＭ、ＥＥＰＲＯＭ等の、不揮発性または揮発性の半導体メモリや、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤ等が該当する。
　メモリは、上述のプログラムの他、軌道決定部８で決定される平均軌道要素および平均軌道要素の時間変化率、目標値設定部１で設定される平均軌道要素の目標値、制御量計算部２で計算される平均軌道要素のフィードフォワード制御量およびフィードバック制御量、分配部５で計算されるおのおののスラスタの噴射量および噴射タイミングを記憶する。上述のメモリに記憶された情報は適宜読みだされ処理に用いられ、修正される。
　なお、目標値設定部１、制御量計算部２、分配部５および軌道決定部８の各機能について、一部を専用のハードウェアで実現し、一部をソフトウェアまたはファームウェアで実現するようにしてもよい。例えば、軌道決定部８については専用のハードウェアとしての処理回路でその機能を実現し、 分配部５については処理回路がメモリに格納されたプログラムを読み出して実行することによってその機能を実現することが可能である。なお、上記において、図８では、処理回路１０１は、図８のプロセッサ１０２に当たり、メモリは、図８のメモリ１０３に当たる。
　このように、処理回路は、ハードウェア、ソフトウェア、ファームウェア、またはこれらの組み合わせによって、上述の各機能を実現することができる。以下、軌道制御装置９の各部の動作を説明する。
　目標値設定部１では、衛星７の平均軌道要素の目標値を設定する。平均軌道要素としては、例えば平均平均直下点経度ΛＭ、平均平均直下点経度角速度Λ'Ｍ、平均平均直下点経度角加速度Λ"Ｍ、平均離心率ベクトルｅｘ、ｅｙ、平均傾斜角ベクトルｉｘ、ｉｙが挙げられる。ただし、平均直下点経度λＭは、衛星の昇交点赤経Ω、近地点引数ω、平均近点離角Ｍ、グリニッジ緯度引数θｇとして、λＭ＝ω＋Ω＋Ｍ－θｇで定義され、平均平均直下点経度ΛＭとは、平均直下点経度λＭの時間平均をとったものである。なお平均軌道要素として他のパラメータを用いてもよいことは言うまでもない。衛星７の平均離心率ベクトルの目標値をｅｘｒｅｆ、ｅｙｒｅｆ、および平均傾斜角ベクトルの目標値をｉｘｒｅｆ、ｉｙｒｅｆ、および平均平均直下点経度の目標値をΛＭｒｅｆと記述する。ここで、角度の単位は、特に明示しない場合には〔rad〕とする。
　平均離心率ベクトルの目標値ｅｘｒｅｆとｅｙｒｅｆの設定方法は、例えばｅｘｒｅｆ＝０、ｅｙｒｅｆ＝０と設定することが挙げられる。その他の設定方法として、離心率ベクトルが太陽輻射圧の影響により１年間で半径ｅｎの円状の軌跡を描くという特性を利用することも可能である。ここで、半径ｅｎの大きさは、衛星７の質量、有効断面積、および衛星表面の光学定数によって定まる。半径ｅｎの大きさが離心率として許容できる場合には、原点を中心とした半径ｅｎの円に沿って、平均離心率ベクトルの目標値ｅｘｒｅｆとｅｙｒｅｆを設定することが可能であり、これにより制御量計算部２における平均離心率ベクトルの制御量を抑えることができる。平均傾斜角ベクトルの目標値ｉｘｒｅｆとｉｙｒｅｆの設定方法は、例えばｉｘｒｅｆ＝０、ｉｙｒｅｆ＝０と設定することが挙げられる。
また、平均平均直下点経度の目標値ΛＭｒｅｆの設定方法は、例えば衛星７を日本の上空に保持する場合にはΛＭｒｅｆ＝１４０degと設定することが挙げられる。
　制御量計算部２は、平均軌道要素、平均軌道要素の時間変化率および目標値から、平均軌道要素の制御量を計算する。制御量計算部２は、フィードフォワード制御量計算部３とフィードバック制御量計算部４から構成される。フィードフォワード制御量計算部３では、軌道決定部８で計算した平均軌道要素の時間変化率をもとに、衛星７の平均離心率ベクトルと平均傾斜角ベクトルのフィードフォワード制御量を計算する。フィードバック制御量計算部４では、軌道決定部８で計算した平均軌道要素および目標値設定部１で設定した平均軌道要素の目標値を元に、衛星７の平均離心率ベクトルと平均傾斜角ベクトルのフィードバック制御量を計算する。以下、フィードフォワード制御量計算部３とフィードバック制御量計算部４の動作について説明する。
　フィードフォワード制御量計算部３では、軌道決定部８で計算した平均軌道要素の時間変化率をもとに、衛星７の平均離心率ベクトルと平均傾斜角ベクトルのフィードフォワード制御量を計算する。平均離心率ベクトルの時刻ｔ０における時間変化率をｅ'ｘ０、ｅ'ｙ０とし、平均傾斜角ベクトルの時刻ｔ０における時間変化率をｉ'ｘ０、ｉ'ｙ０と定義する。このとき、平均離心率ベクトルと平均傾斜角ベクトルの軌道１周期あたりの変化量Δｅ、Δｉは以下の式（１）、（２）で表される。右肩の添え字「Ｔ」は転置を表す。
　　Δｅ＝［ｅ'ｘ０Ｔ　ｅ'ｙ０Ｔ］Ｔ　　　　（１）
　　Δｉ＝［ｉ'ｘ０Ｔ　ｉ'ｙ０Ｔ］Ｔ　　　　（２）
ここで、静止衛星の場合、衛星７の軌道運動の周期Ｔは、地球の自転周期に等しい。平均離心率ベクトルと平均傾斜角ベクトルが軌道１周期あたりにΔｅ、Δｉだけ変化することを見越して、Δｅ、Δｉを定数倍したものを平均離心率ベクトルと平均傾斜角ベクトルのフィードフォワード制御量に設定する。平均離心率ベクトルと平均傾斜角ベクトルそれぞれのフィードフォワード制御量δｅＦＦ、δｉＦＦを、例えば以下の式（３）、（４）のように設定することができる。
　　δｅＦＦ＝－（３／４）Δｅ　　　　（３）
　　δｉＦＦ＝－（３／４）Δｉ　　　　（４）
　以上が、フィードフォワード制御量計算部３の動作の説明である。続いて、フィードバック制御量計算部４の動作について説明する。
　フィードバック制御量計算部４では、軌道決定部８で計算した平均軌道要素および目標値設定部１で設定した平均軌道要素の目標値を元に、衛星７の平均離心率ベクトルと平均傾斜角ベクトルのフィードバック制御量を計算する。時刻ｔ０における平均離心率ベクトルをｅ０＝［ｅｘ０　ｅｙ０］Ｔとし、時刻ｔ０における平均傾斜角ベクトルをｉ０＝［ｉｘ０　ｉｙ０］Ｔと表す。また、目標値設定部１で設定した平均離心率ベクトルと平均傾斜角ベクトルの目標値をｅｒｅｆ＝［ｅｘｒｅｆ　ｅｙｒｅｆ］Ｔ、ｉｒｅｆ＝［ｉｘｒｅｆ　ｉｙｒｅｆ］Ｔと表す。そして、ｅ０とｅｒｅｆの差分をＤｅ＝ｅ０－ｅｒｅｆと定義し、ｉ０とｉｒｅｆの差分をＤｉ＝ｉ０－ｉｒｅｆと定義する。平均離心率ベクトルと平均傾斜角ベクトルのフィードバック制御量をそれぞれδｅＦＢ、δｉＦＢと定義すると、δｅＦＢとδｉＦＢは、ＤｅとＤｉを基に決定される。例えば以下の式（５）、（６）のようにＰＩＤ制御として決定することが可能である。
Figure JPOXMLDOC01-appb-M000001
 　式（５）、（６）において、ｋｅｐ、ｋｅｄおよびｋｅｄはそれぞれ、平均離心率ベクトルのフィードバック制御量を計算するための比例ゲイン、微分ゲインおよび積分ゲインである。ｋｉｐ、ｋｉｄおよびｋｉｉはそれぞれ、平均傾斜角ベクトルのフィードバック制御量を計算するための比例ゲイン、微分ゲインおよび積分ゲインである。
　以上が、フィードバック制御量計算部４の動作の説明である。制御量計算部２が出力する制御量は、フィードフォワード制御量計算部３で計算したフィードフォワード制御量と、フィードバック制御量計算部４で計算したフィードバック制御量を足したものである。
制御量計算部２が出力する平均離心率ベクトルと平均傾斜角ベクトルの制御量をそれぞれδｅ、δｉと定義すると、δｅとδｉは以下の式（７）、（８）のように表される。
　　δｅ＝δｅＦＦ＋δｅＦＢ　　　　（７）
　　δｉ＝δｉＦＦ＋δｉＦＢ　　　　（８）
　平均離心率ベクトルの制御は、衛星７の軌道面内方向の制御であり、平均傾斜角ベクトルの制御は、衛星７の軌道面外方向の制御である。
　制御量計算部２では、式（１）から式（８）の演算に従って、平均離心率ベクトルと平均傾斜角ベクトルの制御量を計算する。制御量計算部２では、フィードフォワード制御量によって高精度な軌道制御を実現し、フィードバック制御によって制御系の安定性・ロバスト性を高めている。以上が制御量計算部２の動作の説明であるが、以下、フローチャートにて説明する。
　図３は、実施の形態１に係る制御量計算部２の動作を表すフローチャートである。制御量計算部２は、フィードフォワード制御量計算部３とフィードバック制御量計算部４から構成される。以下フローチャートに沿って説明する。ステップＳ２５０は、フィードフォワード制御量計算部３が、軌道決定部８から平均軌道要素の時間変化率を受け取る（ステップＳ２５０）。次に、ステップＳ２５１は、平均離心率ベクトルと平均傾斜角ベクトルの軌道１周回あたりの変化量Δｅ、Δｉを算出する（Ｓステップ２５１）。次に、ステップＳ２５２は、Δｅ、Δｉを定数倍したものをフィードフォワード制御量δｅＦＦ，δｉＦＦに設定する（ステップＳ２５２）。以上のステップＳ２５０からＳ２５２までは、フィードフォワード制御量計算部３が実行する。フィードフォワード制御量計算部３は、処理回路１０１、またはプロセッサ１０２で処理してもよいし、独立した処理回路で実現してもよい。また、フィードフォワード制御量計算部３は、平均軌道要素の時間変化率等をメモリから読出し、算出した変化量Δｅ、Δｉ、フィードフォワード制御量δｅＦＦ，δｉＦＦをメモリに記憶する。
　次に、フィードバック制御量計算部４の動作に移る。まず、ステップＳ２５３は、フィードバック制御量計算部４が、軌道決定部８から平均軌道要素を受け取り、目標値設定部１から平均軌道要素の目標値を受け取る（ステップＳ２５３）。次に、ステップＳ２５４は、平均軌道要素の軌道決定値と目標値の差分Ｄｅ，Ｄｉを算出する（ステップＳ２５４）。次に、ステップＳ２５５は、平均離心率ベクトルと平均傾斜角ベクトルのフィードバック制御量δｅＦＢ、δｉＦＢを算出する（ステップＳ２５５）。以上のステップＳ２５３からＳ２５５までは、フィードバック制御量計算部４が実行する。フィードバック制御量計算部４は、処理回路１０１、またはプロセッサ１０２で処理してもよいし、独立した処理回路で実現してもよい。また、フィードバック制御量計算部４は、平均軌道要素、平均軌道要素の目標値等をメモリから読出し、算出した差分Ｄｅ，Ｄｉ、フィードバック制御量δｅＦＢ、δｉＦＢをメモリに記憶する。
　次に、ステップＳ２５６は、フィードフォワード制御量とフィードバック制御量は最終的に加算され、平均離心率ベクトルと平均傾斜角ベクトルの制御量が算出される（ステップＳ２５６）。さらに、ステップＳ２５７では、算出された平均軌道要素の制御量が、分配部５へ渡される（ステップＳ２５７）。以上のステップＳ２５６とＳ２５７とは、制御量計算部２によって実行される。制御量計算部２は、平均離心率ベクトルと平均傾斜角ベクトルの制御量をメモリに記憶する。
　軌道制御装置９の分配部５では、制御量計算部２で計算した平均離心率ベクトルおよび平均傾斜角ベクトルの制御量、目標値設定部１で設定した平均直下点経度の目標値、ならびに軌道決定部８で計算した平均軌道要素および平均軌道要素の時間変化率を受け取る。そして、平均離心率ベクトル（軌道面内方向）と平均傾斜角ベクトル（軌道面外方向）の制御量および平均直下点経度を目標値近傍に維持するための平均直下点経度（軌道面内方向）の制御量を実現するスラスタの噴射タイミングおよび噴射量を計算する。
　分配部５はまず、制御量計算部２で計算した時刻t０における平均傾斜角ベクトルの制御量δｉの大きさδｉ＝|δｉ|と偏角βを計算する。すなわち、δｉを以下の式（９）のように表現する。
　　δｉ＝［δｉcosβ　δｉsinβ］Ｔ　　　　（９）
ただしδｉ＞０を満たす。
　次に分配部５は、１回目の南北制御（軌道面外制御）のタイミングの判定を行う。ここで、経度の変動を補正する制御は一般に東西制御と呼ばれ、緯度の変動を補正する制御は一般に南北制御と呼ばれる。傾斜角ベクトルに関する方程式（例えばガウスの惑星方程式）より、南北制御を行う最適なタイミングは、衛星７の平均経度がβもしくはβ＋πと一致する場合であることが示される。πは円周率である。それ故に、時刻ｔ０における衛星７の平均経度λ０がβもしくはβ＋πと一致する場合に、時刻ｔ０にて１回目の南北制御を行う。平均経度λ０がβにもβ＋πにも一致しない場合には、時刻ｔ０にては噴射量の計算を行わず、南北制御を行わない。以下では、平均経度λ０がβに一致する場合と平均経度λ０がβ＋πに一致する場合の２通りに分けて、分配部５の動作について説明する。
　平均経度λ０がβに一致する場合、分配部５は１回目と２回目の南北制御のタイミングと噴射量を計算する。１回目と２回目の南北制御を行うときの平均経度をそれぞれλ１、λ２と記すと、λ１とλ２は以下の式（１０）、（１１）で与えられる。
　　λ１＝β　　　　　　（１０）
　　λ２＝β＋π　　　　（１１）
　１回目の南北制御ではＳＷスラスタとＳＥスラスタを使用し、２回目の南北制御ではＮＷスラスタとＮＥスラスタを使用する。１回目の南北制御におけるＳＷスラスタとＳＥスラスタの噴射量の合計をｆＳとし、２回目の南北制御におけるＮＷスラスタとＮＥスラスタの噴射量の合計をｆＮと定義する。このとき、ガウスの惑星方程式（Walker M J H, Ireland B, Owens J., A set of modified equinoctial orbit elements. Celestial Mechanics 1985, 36(4): 409-419.）より、ｆＳとｆＮは以下のように定まる。
　まず制御量計算部２において計算した平均離心率ベクトルの制御量δｅの成分をδｅ＝［δｅｘ　δｅｙ］Ｔと表す。また衛星７の軌道長半径ａと平均軌道角速度ｎを用いて、δｅとδｉに依存する変数Ａ、Ｂを以下の式（１２）、（１３）のように定義する。
　　Ａ＝ｎａδｉ（－δｅｘcosβ＋δｅｙsinβ）／（sinθcosφ）　（１２）
　　Ｂ＝ｎａδｉ／cosθ　　　　　　　　　　　　　　　　　　　（１３）
　ここで、角度θはスラスタ９１、９２、９３および９４がＹＢ軸と成す角度であり、角度φは、スラスタ９１、９２、９３および９４の各スラスタと衛星７の質量中心を結ぶ線分をＸＢＺＢ平面に投影した線分が、ＸＢ軸と成す角度である（図２）。
　ここで、上記変数ＡとＢとは、速度の次元を有する。変数Ｂは、傾斜角ベクトルに関するガウス惑星方程式にＳＷスラスタとＳＥスラスタの噴射で生じるYB軸方向の力積の成分を代入することで得られる。また、上記変数Ａは、離心率ベクトルに関するガウス惑星方程式にＳＷスラスタとＳＥスラスタの噴射で生じるXB軸方向の力積の成分を代入することで得られる。
　分配部５は、ＡとＢの大きさに応じて、ｆＳとｆＮを以下の式（１４）のように与える。
Figure JPOXMLDOC01-appb-M000002
　次に分配部５は、平均離心率ベクトルを制御（軌道面内制御）するための制御量の計算を行う。平均離心率ベクトルの制御は、ＮＷスラスタとＳＥスラスタを組み合わせてそれぞれ等しく噴射すること、またはＮＥスラスタとＳＷスラスタを組み合わせてそれぞれ等しく噴射することで、衛星７に軌道半径方向の力を加えることで達成する。平均離心率ベクトルの制御を行うときの平均経度をλｒとし、そのときの衛星に加わる軌道半径方向の制御量をｆｒと定義する。また、偏角β方向の単位ベクトルをρβ＝［cosβ　sinβ］Ｔと定義する。単位ベクトルρβと平均離心率ベクトルの制御量δｅのなす角に応じて、λｒを以下の式（１５）で与える。
Figure JPOXMLDOC01-appb-M000003
 　ここで、δｅ・ρβはδｅとρβの内積を表す。また、ｆｒを以下の式（１６）で与える。
　　ｆｒ＝－ｎａ|δｅ・ρβ|　　　　（１６）
　次に分配部５は、１回目の南北制御におけるＳＷスラスタとＳＥスラスタの噴射量の配分率αＳ、および２回目の南北制御におけるＮＷスラスタとＮＥスラスタの噴射量の配分率αＮを、拘束条件付き最適化問題を解くことで決定する。最適化問題における評価関数Ｊを、例えば以下の式（１７）のように設定することができる。
　　Ｊ＝|（２αＳ－１）ｆＳ－（２αＮ－１）ｆＮ|　　　　（１７）
　式（１７）で表される評価関数は、1回目の南北制御で発生するZB軸方向の力積と，2回目の南北制御で発生するZB軸方向の力積の差が小さくなることを要請している。
　また、拘束条件を例えば以下の式（１８）のように設定できる。
　　ΛＭＴ＝ΛＭｒｅｆ－（１／２）ΔΛＭｔｏｔａｌ　　　　（１８）
　式（１８）は、平均直下点経度を目標値ΛＭｒｅｆに制御するための拘束条件である。式（１８）において、ΔΛＭｔｏｔａｌは、平均経度λ１における１回目の南北制御と平均経度λ２における２回目の南北制御と平均経度λｒにおける平均離心率ベクトル制御によって生じる衛星７の平均平均直下点経度の変化量の総和であり、以下の式（１９）のように表せる。
　　ΔΛＭｔｏｔａｌ＝２tanθcosφδｉ＋２|δｅ・ρβ|　　　　（１９）
　式（１９）は、平均直下点経度に関するガウス惑星方程式に、平均経度λ１における１回目の南北制御、平均経度λ２における２回目の南北制御、および平均経度λｒにおける平均離心率ベクトル制御によって生じるXB軸方向の力積を代入することで得られる。式（１８）において、ΛＭｒｅｆは衛星７の直下点経度の目標保持値であり、目標値設定部１にて設定される。また、式（１８）において、ΛＭＴは１回目の南北制御を行ってから軌道１周期後の衛星７の平均平均直下点経度であり、以下の式（２０）のように表せる。
ΛＭＴ＝ΛＭ０＋Λ'Ｍ０Ｔ＋（１／２）Λ"Ｍ０Ｔ２＋ΔΛＭｔｏｔａｌ
　　－（３／ａ）sinθsinφ［（２αＳ－１）ｆＳ＋（１／２）（２αＮ－１）ｆＮ］Ｔ
　　　　　（２０）
　式（２０）において、最初の３項は，１回目の南北制御を実行する直前の平均平均直下点経度の角速度と角加速度を時間積分することで得られる。また、最後の項は、１回目と２回目の南北制御で生じる平均平均直下点経度の角速度を時間積分することで得られる。ここで、Λ'Ｍ０、Λ'Ｍ０、Λ"Ｍ０はそれぞれ時刻ｔ０における衛星７の平均平均直下点経度、平均平均直下点経度の角速度、および平均平均直下点経度の角加速度である。これらは軌道決定部８にて計算される。分配部５は、式（１８）で与えられる平均直下点経度を目標値ΛＭｒｅｆに制御するための拘束条件を満たしつつ、式（１７）で与えられる評価関数Ｊが最小となるように、配分率αＳとαＮを決定する。
　次に分配部５は、計算した配分率αＳとαＮから、１回目と２回目の南北制御の噴射量を以下の式（２１）、（２２）のように決定する。式（２１）、（２２）で、噴射量ｆの添え字ＮＷ、ＮＥ、ＳＷ、ＳＥはそれぞれ、ＮＷスラスタ、ＮＥスラスタ、ＳＷスラスタ、ＳＥスラスタを表す。また、添え字の１は１回目を、添え字の２は２回目を示す。式（２１）は、１回目の南北制御における各スラスタの噴射量を示す。式（２２）は、２回目の南北制御における各スラスタの噴射量を示す。
　　ｆＮＷ１＝０、ｆＮＥ１＝０、ｆＳＷ１＝αＳｆＳ、ｆＳＥ１＝（１－αＳ）ｆＳ　　　　（２１）
　　ｆＮＷ２＝αＮｆＮ、ｆＮＥ２＝（１－αＮ）ｆＮ、ｆＳＷ２＝０、ｆＳＥ２＝０　　　　（２２）
　分配部５は、計算した南北制御および平均離心率ベクトル制御の噴射タイミングと噴射量をスラスタ制御部６へ渡す。スラスタ制御部６は、分配部５で計算された噴射タイミングおよび噴射量に従って、スラスタ９１、９２、９３および９４を制御する。
　南北制御のタイミング判定に戻って、平均経度λ０がβ＋πに一致する場合の分配部５の動作を説明する。平均経度λ０がβ＋πに一致する場合、分配部５は、１回目と２回目の南北制御のタイミングと噴射量を計算する。１回目と２回目の南北制御を行うときの平均経度をそれぞれλ１、λ２と記すと、λ１とλ２は以下の式（２３）、（２４）で与えられる。
　　λ１＝β＋π　　　　　（２３）
　　λ２＝β＋２π　　　　（２４）
　１回目の南北制御ではＮＷスラスタとＮＥスラスタを使用し、２回目の南北制御ではＳＷスラスタとＳＥスラスタを使用する。１回目の南北制御におけるＮＷスラスタとＮＥスラスタの噴射量の合計をｆＮとし、２回目の南北制御におけるＳＷスラスタとＳＥスラスタの噴射量の合計をｆＳと定義する。以下、ｆＳとｆＮの計算方法について述べる。
　計算した平均離心率ベクトルの制御量δｅの成分をδｅ＝［δｅｘ　δｅｙ］Ｔと表す。また衛星７の軌道長半径と平均軌道角速度をそれぞれａ、ｎと表す。そして、δｅとδｉに依存する変数Ａ、Ｂを以下の式（２５）、（２６）ように定義する。
　　Ａ＝ｎａδｉ（－δｅｘcosβ＋δｅｙsinβ）／（sinθcosφ）　（２５）
　　Ｂ＝ｎａδｉ／cosθ　　　　　　　　　　　　　　　　　　　（２６）
　角度θはスラスタ９１、９２、９３および９４がＹＢ軸と成す角度であり、角度φはスラスタ９１、９２、９３および９４の各スラスタと衛星７の質量中心を結ぶ線分を向けＸＢＺＢ平面に投影した線分が、ＸＢ軸と成す角度である（図２）。
　分配部５は、ＡとＢの大きさに応じて、ｆＳとｆＮを以下の式（２７）のように与える。
Figure JPOXMLDOC01-appb-M000004
　次に分配部５は、平均離心率ベクトルを制御するための制御量の計算を行う。平均離心率ベクトルの制御は、ＮＷスタスタとＳＥスラスタを組み合わせてそれぞれ等しく噴射すること、もしくはＮＥスラスタとＳＷスラスタを組み合わせてそれぞれ等しく噴射することで、衛星７に軌道半径方向の力を加えることで達成する。平均離心率ベクトルの制御を行うときの平均経度をλｒとし、そのときの衛星に加わる軌道半径方向の制御量をｆｒと定義する。また、偏角β方向の単位ベクトルをρβ＝［cosβ　sinβ］Ｔと定義する。単位ベクトルρβと平均離心率ベクトルの制御量δｅのなす角に応じて、平均離心率ベクトルの制御を行うときの平均経度λｒを以下の式（２８）で与える。
Figure JPOXMLDOC01-appb-M000005
 　ここで、δｅ・ρβはδｅとρβの内積を表す。また、衛星に加わる軌道半径方向（XB軸方向）の制御量ｆｒを以下の式（２９）で与える。
　　ｆｒ＝－ｎａ|δｅ・ρβ|　　　　（２９）
　次に分配部５は、１回目の南北制御におけるＮＷスラスタとＮＥスラスタの噴射量の配分率αＮ、および２回目の南北制御におけるＳＷスラスタとＳＥスラスタの噴射量の配分率αＳを、拘束条件付き最適化問題を解くことで決定する。最適化問題における評価関数Ｊを、例えば以下の式（３０）のように設定することができる。
　　Ｊ＝|（２αＳ－１）ｆＳ－（２αＮ－１）ｆＮ|　　　　（３０）
　また、拘束条件を例えば以下の式（３１）のように設定できる。
　　ΛＭＴ＝ΛＭｒｅｆ－（１／２）ΔΛＭｔｏｔａｌ　　　　（３１）
　式（３１）において、ΔΛＭｔｏｔａｌは、平均経度λ１における１回目の南北制御、平均経度λ２における２回目の南北制御、および平均経度λｒにおける平均離心率ベクトル制御によって生じる衛星７の平均平均直下点経度の変化量の総和であり、以下の式（３２）のように表せる。
　　ΔΛＭｔｏｔａｌ＝２tanθcosφδｉ＋２|δｅ・ρβ|　　　　（３２）
　ΛＭｒｅｆは衛星７の直下点経度の目標保持値であり、目標値設定部１にて設定される。また、式（３１）において、ΛＭＴは１回目の南北制御を行ってから軌道１周期経過した後の衛星７の平均平均直下点経度であり、以下の式（３３）のように表せる。
ΛＭＴ＝ΛＭ０＋Λ'Ｍ０Ｔ＋（１／２）Λ"Ｍ０Ｔ２＋ΔΛＭｔｏｔａｌ
　　　－（３／ａ）sinθsinφ［（２αＮ－１）ｆＮ＋（１／２）（２αＳ－１）ｆＳ］Ｔ
　　　　（３３）
　ここで、ΛＭ０、Λ'Ｍ０およびΛ"Ｍ０はそれぞれ時刻ｔ０における衛星７の平均平均直下点経度、平均平均直下点経度の角速度、および平均平均直下点経度の角加速度であり、これらは軌道計算部８にて計算される。分配部５は、式（３１）で与えられる拘束条件を満たしつつ、式（３０）で与えられる評価関数Ｊが最小となるように、配分率αＳとαＮを決定する。
　分配部５は、計算した配分率αＳとαＮから、１回目と２回目の南北制御の噴射量を以下の式（３４）、（３５）のように決定する。噴射量ｆの添え字は、式（２１）および（２２）と同じである。式（３４）は、１回目の南北制御における各スラスタの噴射量を示す。式（３５）は、２回目の南北制御における各スラスタの噴射量を示す。
　　ｆＮＷ１＝αＮｆＮ、ｆＮＥ１＝（１－αＮ）ｆＮ、ｆＳＷ１＝０、ｆＳＥ１＝０　　　　（３４）
　　ｆＮＷ２＝０、ｆＮＥ２＝０、ｆＳＷ２＝αＳｆＳ、ｆＳＥ２＝（１－αＳ）ｆＳ　　　　（３５）
　式（１２）～（１４）、（１６）および（１８）～（２０）は、軌道面外運動および軌道面内運動の連成を考慮した方程式である。分配部５は複数回のスラスタ噴射量による面外運動および面内運動の連成を考慮した方程式を解き、主として軌道面外方向を制御するスラスタ噴射および主として軌道面内方向を制御するスラスタ噴射を複数回組み合わせることによって、制御量計算部２で計算した平均軌道要素の制御量を実現させるためのスラスタの噴射タイミングおよび噴射量を計算する。
　分配部５は、以上で計算した南北制御および平均離心率ベクトル制御の各スラスタの噴射タイミングと噴射量をスラスタ制御部６に渡す。スラスタ制御部６は、分配部５で計算された噴射タイミングおよび噴射量に従って、スラスタ９１、９２、９３および９４を制御する。
　図４は実施の形態１に係る分配部５の動作を示すフローチャートである。分配部５はまず、式（９）で表される平均傾斜角ベクトルの制御量の大きさと偏角を計算する（ステップＳ０１）。そして、南北制御のタイミングを判定する（ステップＳ０２）。衛星７の平均経度λ０が偏角βに一致する場合（ステップＳ０２；Ｙ）、式（１０）および式（１１）で南北制御の噴射タイミングλ１、λ２を計算し、それぞれの合計噴射量ｆＳ、ｆＮを式（１２）から式（１４）で計算する（ステップＳ０３）。分配部５は、平均傾斜角ベクトルの制御量の大きさと偏角、南北制御の噴射タイミングλ１、λ２、合計噴射量ｆＳ、ｆＮをメモリに記憶する。
　分配部５は次に、平均離心率ベクトル制御の噴射タイミングλｒと噴射量ｆｒを、式（１５）および式（１６）で計算する（ステップＳ０４）。そして、式（１８）の拘束条件の下で、式（１７）の評価関数Ｊを最小にする配分率αＳとαＮを決定する（ステップＳ０５）。分配部５は、配分率αＳとαＮから、式（２１）および式（２２）で１回目と２回目の各スラスタの噴射量を計算する（ステップＳ０６）。そして、計算した噴射タイミングと噴射量をスラスタ制御部６に送る（ステップＳ０７）。分配部５は、平均離心率ベクトル制御の噴射タイミングλｒと噴射量ｆｒ、配分率αＳとαＮ、および噴射タイミングと噴射量をメモリに記憶する。
　ステップＳ０２で平均経度λ０が偏角βに一致せず（ステップＳ０２；Ｎ）、λ０がβ＋πに一致する場合（ステップＳ０８；Ｙ）、分配部５は、式（２３）および式（２４）で南北制御の噴射タイミングλ１、λ２を計算し、それぞれの合計噴射量ｆＳ、ｆＮを式（２５）から式（２７）で計算する（ステップＳ０９）。分配部５は、南北制御の噴射タイミングλ１、λ２、および合計噴射量ｆＳ、ｆＮをメモリに記憶する。
　分配部５は次に、平均離心率ベクトル制御の噴射タイミングλｒと噴射量ｆｒを、式（２８）および式（２９）で計算する（ステップＳ１０）。そして、式（３１）の拘束条件の下で、式（３０）の評価関数Ｊを最小にする配分率αＳとαＮを決定する（ステップＳ１１）。分配部５は、配分率αＳとαＮから、式（３４）および式（３５）で１回目と２回目の各スラスタの噴射量を計算する（ステップＳ１２）。そして、計算した噴射タイミングと噴射量をスラスタ制御部６に送る（ステップＳ０７）。分配部５は、平均離心率ベクトル制御の噴射タイミングλｒと噴射量ｆｒ、配分率αＳとαＮおよび１回目と２回目の各スラスタの噴射量をメモリに記憶する。
　ステップＳ０２およびステップＳ０８で、平均経度λ０が偏角βにもβ＋πにも一致しない場合（ステップＳ０２；Ｎ，ステップＳ０８；Ｎ）、噴射量の計算も噴射も行わず（ステップＳ１３）、その回の処理を終了する。
　以上の実施の形態１の構成によれば、分配部５において、スラスタ配置角度と複数回のスラスタ噴射に依存した軌道面内と軌道面外の連成を含んだガウスの惑星方程式を解析的に解くことによって、消費推薬を抑えたスラスタ噴射則を少ない計算コストで算出することができる。
　また、スラスタの同時噴射を１本もしくは２本までに制限した噴射則を計算することができ、消費電力を抑制することができる。具体的には、分配部５によって電気推進装置の同時噴射を１本もしくは２本までに制限した噴射則を少ない計算量で得ることができるので、消費電力を抑制しつつ、かつ高精度に衛星の軌道を保持することができる。そして、平均傾斜角ベクトルの制御量をもとに最適な南北制御の噴射タイミングを計算することができ、高精度な南北制御を実現することができる。また、平均平均直下点経度、平均傾斜角ベクトル、平均離心率ベクトルの制御を南北制御において同時に行うことが可能となる。さらに、２回の南北制御に加えて、平均離心率ベクトルの制御を軌道１周回あたりに１回行うことで、平均離心率ベクトルを目標値近傍に保持することが可能となり、衛星７の直下点経度を目標値近傍に高精度に保持することが可能となる。実施の形態１によれば、南北制御と平均離心率ベクトル制御の間隔を軌道１／４周期あけることが可能となり、南北制御と平均離心率ベクトル制御のスラスタ噴射のタイミングが重なることを回避できる。
　実施の形態２．
　実施の形態２では、軌道一周回においてスラスタ噴射が可能なタイミングを予め固定しておく。そして、その各タイミングについて制御量計算部で計算された平均離心率ベクトルおよび平均傾斜角ベクトルの制御量に応じて２回～４回のスラスタの噴射量および噴射の可否を決定する。実施の形態２の構成は実施の形態１と同様であるが、軌道制御装置９における分配部５の動作のみが異なる。その他の、目標値設定部１、制御量計算部２、スラスタ制御部６および軌道決定部８は実施の形態１と共通である。
　分配部５は、実施の形態１と同様に、制御量計算部２で計算した平均離心率ベクトルおよび平均傾斜角ベクトルの制御量、目標値設定部１で設定した平均直下点経度の目標値、ならびに軌道決定部８で計算した平均軌道要素および平均軌道要素の時間変化率を受け取る。そして、平均離心率ベクトルと平均傾斜角ベクトルの制御量および平均直下点経度を目標値近傍に維持するための平均直下点経度の制御量を実現するスラスタの噴射タイミングおよび噴射量を計算する。
　このとき実施の形態２では、スラスタの噴射タイミングを１周期において最大４回とする。それぞれのタイミングは平均経度λを用いて、λ１、λ２、λ３、λ４と表し次の式（４８）～（５１）で与える。
　　λ１＝λＣ　　　　　　　　　（４８）
　　λ２＝λ１＋π／２　　　　　（４９）
　　λ３＝λ１＋π　　　　　　　（５０）
　　λ４＝λ１＋３π／２　　　　（５１）
　λＣは予め与える１回目のスラスタ噴射タイミングである。実施の形態２においてはこのλＣを昇交点Ω付近の適当な値に固定する。例えばλＣ＝８０［deg］と設定する。そして１回目のスラスタ噴射から等間隔に位相をπ／２ずつ離してスラスタ噴射を行う。但しこのλＣおよびλ１、λ２、λ３、λ４は上記の与え方以外の値を選ぶことも可能である。ここではスラスタの噴射回数を４回としているが実際に計算された制御量から、スラスタ噴射の回数が２回に抑えられる場合もありうる。以下に　軌道の保持に必要な制御力とスラスタの噴射量の計算方法について説明する。
　はじめに分配部５は面外方向の制御力を計算する。まず、制御量計算部２において計算された平均離心率ベクトルの制御量δｅ＝［δｅｘ　δｅｙ］Ｔと平均傾斜角ベクトルの制御量δｉ＝［δｉｘ　δｉｙ］Ｔを受け取る。また目標値設定部１より平均平均直下点経度の目標値ΛＭＴｒｅｆを受け取る。これらを元に分配部５は各スラスタ噴射タイミングにおける面外方向の制御力ｆｈｉ（ｉ＝１～４）を計算する。制御力の計算に当たって次の単位ベクトルρｅｉ＝［sinλｉ　－cosλｉ］Ｔ 、ρｉｉ＝［cosλｉ　sinλｉ］Ｔを用いる。すると各タイミングにおける面外方向の制御力ｆｈｉは次の二つの式を満たすように与えられる。
　　|ｆｈｉ|－|ｆｈ（ｉ＋２）|＝ａｎδｅ・ρｅｉ／（tanθcosφ）　（５２）
　　ｆｈｉ－ｆｈ（ｉ＋２）＝ａｎδｉ・ρｅｉ　　　　　　　　　　（５３）
　θおよびφは、図２で与えられる角度である。
　式（５２）の右辺をδｅｉ、式（５３）の右辺をδｉｉとおいて、δｅｉ、δｉｉの大小関係によって場合分けしてｆｈｉを決定する。この場合分けを行う際には、平均離心率ベクトルの制御量を実現することを優先する。スラスタの配置よりｆｈｉが定まると自動的にｆｒｉは決まってしまう。δｅｉ＞δｉｉとなった場合は平均傾斜角ベクトルの制御量に合わせてしまうと、平均離心率ベクトルの制御量が足りなくなってしまう。そのためこのような場合は平均離心率ベクトルの制御量を基準にｆｈｉを決めることとする。
（１）　|δｅｉ|＞|δｉｉ|のとき
　（１ａ）　δｅｉ＜０　 のとき
　　ｆｈｉ＝sign（δｉｉ）|δｅｉ|、ｆｈ（ｉ＋２）＝０　　　　（５４）
　（１ｂ）　δｅｉ＞０のとき
　　ｆｈｉ＝０、ｆｈ（ｉ＋２）＝－sign（δｉｉ）|δｅｉ|　　　　（５５）
（２）　|δｅｉ|＜|δｉｉ|のとき
　（２ａ）　δｅｉ＋δｉｉ＞０かつδｅｉ－δｉｉ＜０ のとき
　　ｆｈｉ＝（１／２）（δｅｉ＋δｉｉ）、ｆｈ（ｉ＋２）＝（１／２）（δｅｉ－δｉｉ）
　　　　（５６）
　（２ｂ）　δｅｉ＋δｉｉ＜０かつδｅｉ－δｉｉ＞０のとき
　　ｆｈｉ＝－（１／２）（δｅｉ＋δｉｉ）、ｆｈ（ｉ＋２）＝－（１／２）（δｅｉ－δｉｉ）
　　　　（５７）
　次に東西制御（直下点経度の保持）に必要なｆθｉ（ｉ＝１～４） を決定する。分配部５は、次の条件式（５８）を満たすように、ｆθｉ（ｉ＝１～４）を決定する。
Figure JPOXMLDOC01-appb-M000006
 　ａは軌道長半径、ｎは平均軌道角速度、θおよびφは図２で与えられる角度である。
　このｆθｉを与えるためのスラスタ噴射は原則として、面外方向の制御力ｆｈｉを与える際にＳＷスラスタとＳＥスラスタ、または、ＮＷスラスタとＮＥスラスタに対する噴射量の配分（電気推進の場合は噴射時間の配分）を行うことで実現する。この配分量は、ｆｈｉとｆθｉが与えられた時、スラスタ配置角にもとづいて一意に定めることができる。ただし実施の形態２ではｆθｉを決定する際に、ｆｈｉは制御サイクルごとに異なり、場合によっては４回の制御では十分な合計のｆθｉを得られない場合も考えうる。そのような場合は一回目のスラスタ噴射においてスラスタ２本もしくは１本による噴射を行い、必要なｆθを追加で与えることとする。分配部５は、以上の計算結果をスラスタ制御部６へ渡す。
　実施の形態２における基本の動作は以上の通りであるが、推薬のさらなる節約を目指す場合は面外制御量ｆｈｉを計算する際に次のような制御を行う。制御後の平均傾斜角ベクトルの値を予測した際に、その値が予め設定しておいた許容範囲内に収まる場合、すなわち平均傾斜角ベクトルの誤差ιが十分小さい値、例えば|ι|＜０．００３［deg］に収まる場は、平均的な離心率ベクトルの制御に必要な制御量にのみ制御量を合わせることとしてスラスタ噴射の回数を３回以下に削減する。
　図５は、実施の形態２に係る分配部５の動作を示すフローチャートである。分配部５は、まず、平均離心率ベクトルおよび平均傾斜角ベクトルの制御量、平均平均直下点経度の目標値などを受け取る（ステップＳ２１）。そして、平均経度λが設定したλ１（＝λｃ）になったとき、式（５２）および（５３）を満たす面外制御量ｆｈｉを、δｅｉ、δｉｉの大小関係によって場合分けして式（５４）～（５７）に従って計算する（ステップＳ２２）。分配部５は、平均離心率ベクトルおよび平均傾斜角ベクトルの制御量、平均平均直下点経度の目標値およびδｅｉ、δｉｉをメモリから読出し、面外制御量ｆｈｉをメモリに記憶する。
　次に、分配部５は、制御量ｆｈｉを元に、東西制御（直下点経度の保持）に必要なｆθｉ（ｉ＝１～４） を式（５８）を満たすように計算する（ステップＳ２３）。Σｆθｉが東西制御で十分かどうか判定し（ステップＳ２４）、十分でなければ（ステップＳ２４；Ｎ）、Σｆθｉが十分となるようにｆθｉを再計算して追加する（ステップＳ２５）。Σｆθｉが東西制御で十分なら（ステップＳ２４；Ｙ）、再計算を行わない。分配部５は、制御量ｆｈｉをメモリから読出し、東西制御（直下点経度の保持）に必要なｆθｉ（ｉ＝１～４）をメモリに記憶する。
　分配部５は、実施の形態１と同様にして、各スラスタ噴射タイミングにおける各スラスタの噴射量の配分を決定し（ステップＳ２６）、固定されている噴射タイミングと各スラスタの噴射量をスラスタ制御部６に送る（ステップＳ２７）。分配部５は、各スラスタ噴射タイミングにおける各スラスタの噴射量の配分をメモリに記憶する。
　以上のステップＳ２１からステップＳ２７にしたがって、実施の形態２では平均離心率ベクトルと平均傾斜角ベクトルの制御量および平均直下点経度を目標値近傍に維持するための平均直下点経度の制御量を実現するスラスタの噴射タイミングおよび噴射量を計算する。
　以上の実施の形態２の構成と動作によれば、各スラスタ噴射間の噴射間隔を固定することができ、制御時間間隔の変動によって生じる経度方向の軌道保持誤差を低減することができる。この実施の形態２は、最も高精度の東西制御（例えば直下点経度誤差Δλ＝０．００５［deg］以下など）が要求される場合に有用である。
　実施の形態３．
　実施の形態３では、４本のスラスタのうち２本のスラスタを組み合わせて軌道１周回あたり１回南北制御を行い、２本のスラスタを組み合わせて軌道１周回あたり２回東西制御を行う。実施の形態３の構成は実施の形態１と同様であるが、軌道制御装置９における分配部５の動作のみが異なる。その他の、目標値設定部１、制御量計算部２、スラスタ制御部６、衛星７および軌道決定部８は実施の形態１と共通である。実施の形態１と異なる部分は分配部５のみであるので、分配部５について述べる。
　図２に示すスラスタ９１（ＮＷスラスタ）、スラスタ９２（ＮＥスラスタ）、スラスタ９３（ＳＷスラスタ）およびスラスタ９４（ＳＥスラスタ）の各スラスタの噴射量をそれぞれｆＮＷ、ｆＮＥ、ｆＳＷ、ｆＳＥと記す。スラスタを２本同時に噴射する場合を前提にすると、北方向噴射量ｆＮ、南方向噴射量ｆＳ、東方向噴射量ｆＥ、西方向噴射量ｆＷをそれぞれ次のように表すことができる。
　　ｆＮ＝ｆＮＷ＋ｆＮＥ、ｆＳ＝ｆＳＷ＋ｆＳＥ、
　　ｆＥ＝ｆＮＥ＋ｆＳＥ、ｆＷ＝ｆＮＷ＋ｆＳＷ
　分配部５は受け取った軌道要素の制御量のうち、平均傾斜角の制御量を北方向噴射量ｆＮまたは南方向噴射量ｆＳによって制御する。東方向噴射量ｆＥや西方向噴射量ｆＷは平均傾斜角には影響を与えないので、北方向噴射量ｆＮあるいは南方向噴射量ｆＳを平均傾斜角の制御量だけから定めることができる。噴射タイミングは、平均傾斜角ベクトルの制御量δｉの偏角βから決めることができる。この北方向噴射量ｆＮあるいは南方向噴射量ｆＳはどちらを用いても平均傾斜角を制御することができるが、この両方の場合を計算し、最終的に噴射量の合計値が小さくなる方を採用する。
　つぎに分配部５は、北方向噴射量ｆＮあるいは南方向噴射量ｆＳによって平均離心率や平均平均直下点経度は影響を受けるのでその分の補正量を加えて、新たな平均離心率ベクトルの制御量と平均平均直下点経度の目標値を定める。この補正量は、北方向噴射量ｆＮまたは南方向噴射量ｆＳのXB軸およびZB軸方向成分を離心率ベクトルおよび平均直下点経度に関するガウスの惑星方程式に代入することで算出することができる。
　分配部５は、新たな平均離心率ベクトルの制御量と平均平均直下点経度の目標値から、東方向噴射量ｆＥと西方向噴射量ｆＷを定める。平均離心率ベクトルの制御量と平均平均直下点経度の目標値に対応する噴射量を定めるためには、東方向噴射量ｆＥと西方向噴射量ｆＷは合計２回分が必要となる。すなわち、ある経度において東方向噴射量ｆＥあるいは西方向噴射量ｆＷを与え、異なる経度において東方向噴射量ｆＥあるいは西方向噴射量ｆＷを与えるような解を求める。その結果、東西方向の噴射においてはｆＥとｆＷが各１回ずつか、ｆＥが２回か、ｆＷが２回かのいずれかの解が得られる。この合計２回の噴射とｆＮあるいはｆＳの１回の噴射によって一連の制御動作は終了する。したがって一連の制御動作において、合計３回の噴射を行うことになる。
　次に分配部５は、北方向噴射量ｆＮあるいは南方向噴射量ｆＳのいずれを採用するかを次のように決定する。分配部５は、ｆＮを用いたときの東西噴射量もあわせた噴射量の合計値と、ｆＳを用いたときの東西噴射量もあわせた噴射量の合計値を比較して、その合計値の小さな方を最終的な噴射として採用する。このように噴射量の少ない方を用いることにより、南北のｆＮあるいはｆＳの１回の噴射と東西の２回の噴射において、それらを実施するときの経度は離れることになるので、これらの噴射のタイミングが重なることはないという効果も得られる。
　このように南北および東西の噴射において２本のスラスタを組み合わせて用いることにより、一連の制御動作において南北１回、東西２回のスラスタ噴射で、平均傾斜角、平均離心率、平均平均直下点経度を望ましい値に制御することができる。
　また、上記のスラスタ噴射量の計算をもとに、東西における２回のスラスタ噴射において、これを東方向や西方向に力を出すような２本のスラスタの同時噴射に代えて、一部あるいは全部を１本のスラスタだけの噴射にすることで、平均傾斜角にも影響を与えるようにすることができる。この場合、東西のスラスタ噴射によって平均傾斜角が影響を受けるため、南北のスラスタ噴射量を簡単に決定することはできないが、スラスタ２本の同時噴射の解をもとに、たとえば全部を１本のスラスタだけの噴射に代えるように数値探索を行うことによって、２回の東西方向のスラスタ噴射において、その両方あるいは一方を１本のスラスタだけの噴射で計画することができる。このように計画すると東西方向のスラスタ噴射で平均傾斜角の一部を制御するため、噴射量の合計値を低く抑える効果がある。
　図６は、実施の形態３に係る分配部５の動作を示すフローチャートである。分配部５は、まず、平均離心率ベクトルおよび平均傾斜角ベクトルの制御量、平均平均直下点経度の目標値などを受け取る（ステップＳ３１）。そして、このうち平均傾斜角ベクトルの制御量から北方向噴射量ｆＮまたは南方向噴射量ｆＳ（南北噴射量ｆＮまたはｆＳ）を決定する（ステップＳ３２）。分配部５は、平均離心率ベクトルおよび平均傾斜角ベクトルの制御量、平均平均直下点経度の目標値などをメモリから読出し、北方向噴射量ｆＮまたは南方向噴射量ｆＳ（南北噴射量ｆＮまたはｆＳ）をメモリに記憶する。
　つぎに分配部５は、決定した南北噴射量ｆＮまたはｆＳより新たな平均離心率の制御量と平均平均直下点経度の目標値を修正する（ステップＳ３３）。そして、新たな平均離心率の制御量と平均平均直下点経度の目標値から、東方向噴射量ｆＥと西方向噴射量ｆＷ（東西噴射量ｆＥ、ｆＷ）を定める（ステップＳ３４）。分配部５は、南北噴射量ｆＮまたはｆＳをメモリから読出し、新たな平均離心率の制御量、平均平均直下点経度の目標値、東方向噴射量ｆＥと西方向噴射量ｆＷ（東西噴射量ｆＥ、ｆＷ）をメモリに記憶する。
　分配部５は、南北噴射量ｆＮまたはｆＳにおいて、東西噴射量を加えた合計噴射量の少ない方を採用する（ステップＳ３５）。さらに、スラスタの噴射量の計算をもとに、東西における２回のスラスタ噴射において、これを東方向や西方向に力を出すような２本のスラスタの同時噴射に代えて、一部あるいは全部を１本のスラスタだけの噴射で代替するように再計算する（ステップＳ３６）。そして、計算したスラスタの噴射量をスラスタ制御部６に送る（ステップＳ３７）。分配部５は、計算したスラスタの噴射量をメモリに記憶する。
　以上の実施の形態３の構成によれば、２本のスラスタを組み合わせた噴射を前提とすることで、簡単な計算により、南北に１回の噴射、東西に２回の噴射によって、平均傾斜角、平均離心率、平均平均直下点経度を望ましい値に制御することが可能になる。また、実施の形態３の構成によれば、２本のスラスタを組み合わせた噴射を前提とした解に基づいて、東西２回のスラスタ噴射の少なくとも１回をスラスタ１本だけの噴射とするような数値探索の解を求めることにより、平均傾斜角、平均離心率、平均平均直下点経度を望ましい値に制御するとともに、噴射量の合計値を低く抑える効果が得られる。
　実施の形態では、静止衛星を想定して説明した。実施の形態の軌道制御装置９が適用できるのは、静止衛星には限らない。軌道制御装置９は、例えば、準天頂衛星にも適用できる。準天頂衛星は、地上の１点から観測すると南北に長い８の字の軌道に見えるが、周期と高度は静止衛星と同じで軌道面が赤道面に対して傾いている。従って、上述の式を適当に座標変換すれば、そのまま適用できる。
　本発明は、本発明の広義の精神と範囲を逸脱することなく、様々な実施の形態および変形が可能とされるものである。また、上述した実施の形態は、本発明を説明するためのものであり、本発明の範囲を限定するものではない。本発明の範囲は、実施の形態ではなく、特許請求の範囲によって示される。そして、特許請求の範囲内およびそれと同等の発明の意義の範囲内で施される様々な変形が、本発明の範囲内とみなされる。
　本出願は、２０１５年１月９日に出願された、明細書、特許請求の範囲、図、および要約書を含む、日本国特許出願２０１５－００２８６４号に基づく優先権を主張するものである。日本国特許出願２０１５－００２８６４号の開示内容は参照により全体として本出願に含まれる。
　１　目標値設定部、２　制御量計算部、３　フィードフォワード制御量計算部、４　フィードバック制御量計算部、５　分配部、６　スラスタ制御部、７　衛星、８　軌道決定部、９　軌道制御装置、９１，９２，９３，９４　スラスタ。

Claims (9)

  1. 　それぞれが衛星の質量中心から遠ざかる向きに互いに異なる方向に噴射方向を向けて前記衛星に配置された４本のスラスタを備える衛星の軌道制御装置であって、
    　前記衛星の平均軌道要素および平均軌道要素の時間変化率を決定する軌道決定部と、
    　前記平均軌道要素の目標値を設定する目標値設定部と、
    　前記平均軌道要素、前記平均軌道要素の時間変化率および前記目標値から、前記平均軌道要素の制御量を計算する制御量計算部と、
    　前記衛星の運動を軌道要素で表現し、スラスタ配置角度と複数回のスラスタ噴射量による軌道面外運動および軌道面内運動の連成を考慮した方程式を解き、主として軌道面外方向を制御するスラスタ噴射および主として軌道面内方向を制御するスラスタ噴射を複数回組み合わせることによって、前記制御量計算部で計算した前記平均軌道要素の制御量を実現させるための前記スラスタの噴射タイミングおよび噴射量を計算する分配部と、
    　前記分配部で計算された前記噴射タイミングおよび前記噴射量に従って、前記スラスタを制御するスラスタ制御部と、
    　を備える軌道制御装置。
  2. 　前記制御量計算部は、
    　前記軌道決定部にて計算された前記平均軌道要素および前記目標値設定部にて設定された前記平均軌道要素の目標値をもとに、前記平均軌道要素のフィードバック制御量を計算するフィードバック制御量計算部と、
    　前記平均軌道要素の時間変化率をもとに、前記平均軌道要素のフィードフォワード制御量を計算するフィードフォワード制御量計算部と、
    　を含む請求項１に記載の軌道制御装置。
  3. 　前記分配部は、同時に使用するスラスタ本数が２本以下になるように前記スラスタの前記噴射タイミングおよび前記噴射量を計算する請求項１に記載の軌道制御装置。
  4. 　前記制御量計算部は、平均傾斜角ベクトルの制御量を計算し、
    　前記分配部は、前記平均傾斜角ベクトルの制御量の偏角を計算し、前記平均傾斜角ベクトルの制御量の偏角から南北制御の噴射タイミングを決定し、１回目の前記南北制御の噴射タイミングから軌道１／２周期後に２回目の南北制御を実施し、前記１回目の前記南北制御の噴射タイミングから軌道１／４周期後もしくは３／４周期後に平均離心率ベクトルの制御を行うように前記スラスタの前記噴射タイミングおよび前記噴射量を計算する、請求項１に記載の軌道制御装置。
  5. 　前記分配部は、４本のうち２本の前記スラスタを組み合わせて軌道１周回あたり１回南北制御を行い、２本の前記スラスタを組み合わせて軌道１周回あたり２回東西制御を行うように前記スラスタの前記噴射タイミングおよび前記噴射量を計算する、請求項１に記載の軌道制御装置。
  6. 　前記分配部は、前記分配部で計算した前記スラスタの前記噴射タイミングと前記噴射量を初期解として数値探索を実行し、２回の前記東西制御のスラスタ噴射のうち、少なくとも１回は前記スラスタを１本使用しないような前記スラスタの前記噴射タイミングおよび前記噴射量を計算する、請求項５に記載の軌道制御装置。
  7. 　前記制御量計算部は、平均離心率ベクトルおよび平均傾斜角ベクトルの制御量を計算し、
    　前記分配部においては、軌道一周回においてスラスタ噴射が可能なタイミングを予め固定し、その各タイミングについて前記制御量計算部で計算された前記平均離心率ベクトルおよび前記平均傾斜角ベクトルの制御量に応じて、２回から４回の前記スラスタの前記噴射量および噴射の可否を決定する、請求項１に記載の軌道制御装置。
  8. 　それぞれが衛星の質量中心から遠ざかる向きに互いに異なる方向に噴射方向を向けて前記衛星に配置された４本のスラスタと、
    　請求項１から７のいずれか１項に記載の軌道制御装置と、
    　を備える衛星。
  9. 　それぞれが衛星の質量中心から遠ざかる向きに互いに異なる方向に噴射方向を向けて前記衛星に配置された４本のスラスタを備える衛星の軌道制御方法であって、
    　前記衛星の平均軌道要素および平均軌道要素の時間変化率を決定する軌道決定ステップと、
    　前記平均軌道要素の目標値を設定する目標値設定ステップと、
    　前記平均軌道要素、前記平均軌道要素の時間変化率および前記目標値から、前記平均軌道要素の制御量を計算する制御量計算ステップと、
    　前記衛星の運動を軌道要素で表現し、スラスタ配置角度と複数回のスラスタ噴射量による軌道面外運動および軌道面内運動の連成を考慮した方程式を解き、主として軌道面外方向を制御するスラスタ噴射および主として軌道面内方向を制御するスラスタ噴射を複数回組み合わせることによって、前記制御量計算ステップで計算した前記平均軌道要素の制御量を実現させるための前記スラスタの噴射タイミングおよび噴射量を計算する分配ステップと、
    　前記分配ステップで計算された前記噴射タイミングおよび前記噴射量に従って、前記スラスタを制御するスラスタ制御ステップと、
    　を備える軌道制御方法。
PCT/JP2016/050262 2015-01-09 2016-01-06 軌道制御装置および衛星 WO2016111317A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
EP16735056.0A EP3243756B1 (en) 2015-01-09 2016-01-06 Orbit control device and satellite
US15/540,726 US10752384B2 (en) 2015-01-09 2016-01-06 Orbit control device and satellite
JP2016568739A JP6271043B2 (ja) 2015-01-09 2016-01-06 軌道制御装置および衛星

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2015-002864 2015-01-09
JP2015002864 2015-01-09

Publications (1)

Publication Number Publication Date
WO2016111317A1 true WO2016111317A1 (ja) 2016-07-14

Family

ID=56356003

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2016/050262 WO2016111317A1 (ja) 2015-01-09 2016-01-06 軌道制御装置および衛星

Country Status (4)

Country Link
US (1) US10752384B2 (ja)
EP (1) EP3243756B1 (ja)
JP (1) JP6271043B2 (ja)
WO (1) WO2016111317A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2557011A (en) * 2016-11-28 2018-06-13 National Univ Of Defense Technology Spark-based imaging satellite task preprocessing parallelization method
WO2019044735A1 (ja) * 2017-09-01 2019-03-07 三菱電機株式会社 宇宙機制御装置、宇宙機制御方法、およびプログラム
JPWO2021171409A1 (ja) * 2020-02-26 2021-09-02
US11167867B2 (en) * 2016-09-09 2021-11-09 Mitsubishi Electric Corporation Artificial satellite, attitude control system, and attitude control method

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR3043985B1 (fr) * 2015-11-20 2018-04-20 Thales Procede de transfert orbital d'un vaisseau spatial utilisant une poussee continue ou quasi-continue et systeme embarque de pilotage pour la mise en oeuvre d'un tel procede
CN108490973B (zh) * 2018-04-19 2021-04-13 哈尔滨工业大学 航天器编队相对轨道确定方法及装置
CN108873953B (zh) * 2018-08-28 2021-09-07 北京控制工程研究所 一种基于电磁比例阀的高精度压力控制方法及***
CN109343554B (zh) * 2018-11-02 2020-08-21 北京理工大学 一种基于状态转换代价值的启发式航天器任务规划方法
CN112124626B (zh) * 2020-08-27 2022-02-15 中国人民解放军战略支援部队航天工程大学 一种Walker星座构型维持方法和终端设备
CN113184220B (zh) * 2021-04-21 2021-11-19 中国人民解放军63923部队 一种地球同步轨道通信卫星的轨道控制方法及装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5443231A (en) * 1993-11-17 1995-08-22 Hughes Aircraft Company Method and apparatus for a satellite station keeping
US20090020650A1 (en) * 2007-07-17 2009-01-22 Ho Yiu-Hung M System and methods for simultaneous momentum dumping and orbit control

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US3907226A (en) * 1970-07-06 1975-09-23 Hughes Aircraft Co Redundant position and attitude control for spin stabilized devices
JP2518212B2 (ja) * 1986-06-26 1996-07-24 日本電気株式会社 人工衛星の軌道制御方式
JPH03159898A (ja) 1989-11-17 1991-07-09 Toshiba Corp 宇宙航行体の静止軌道保持装置
US5124925A (en) * 1990-01-16 1992-06-23 Space Systems/Loral, Inc. Method for controlling east/west motion of a geostationary satellite
US5813633A (en) * 1995-12-22 1998-09-29 Hughes Electronics Corporation Method and apparatus for stationkeeping a satellite offset by pitch rotation
US5810295A (en) * 1996-07-10 1998-09-22 Hughes Electronics Corporation Method and apparatus for a satellite station keeping
US6135394A (en) * 1998-12-08 2000-10-24 Space Systems/Loral, Inc. Practical method and apparatus for satellite stationkeeping
US6439507B1 (en) * 2000-05-05 2002-08-27 Space Systems/Loral, Inc. Closed-loop spacecraft orbit control

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5443231A (en) * 1993-11-17 1995-08-22 Hughes Aircraft Company Method and apparatus for a satellite station keeping
US20090020650A1 (en) * 2007-07-17 2009-01-22 Ho Yiu-Hung M System and methods for simultaneous momentum dumping and orbit control

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
LOSA, DAMIANA ET AL.: "Electric Station Keeping of Geostationary Satellites: a Differential Inclusion Approach", PROCEEDINGS OF THE 44TH IEEE CONFERENCE ON DECISION AND CONTROL, AND THE EUROPEAN CONTROL CONFERENCE 2005, December 2005 (2005-12-01), pages 7484 - 7489, XP010884902 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11167867B2 (en) * 2016-09-09 2021-11-09 Mitsubishi Electric Corporation Artificial satellite, attitude control system, and attitude control method
GB2557011A (en) * 2016-11-28 2018-06-13 National Univ Of Defense Technology Spark-based imaging satellite task preprocessing parallelization method
GB2557011B (en) * 2016-11-28 2020-12-09 National Univ Of Defense Technology Spark-based imaging satellite task preprocessing parallelization method
WO2019044735A1 (ja) * 2017-09-01 2019-03-07 三菱電機株式会社 宇宙機制御装置、宇宙機制御方法、およびプログラム
JPWO2019044735A1 (ja) * 2017-09-01 2020-04-23 三菱電機株式会社 宇宙機制御装置、宇宙機制御方法、およびプログラム
JPWO2021171409A1 (ja) * 2020-02-26 2021-09-02
WO2021171409A1 (ja) * 2020-02-26 2021-09-02 三菱電機株式会社 軌道姿勢制御装置、人工衛星、軌道姿勢制御方法及びプログラム
JP7204987B2 (ja) 2020-02-26 2023-01-16 三菱電機株式会社 軌道姿勢制御装置、人工衛星、軌道姿勢制御方法及びプログラム

Also Published As

Publication number Publication date
EP3243756A1 (en) 2017-11-15
JPWO2016111317A1 (ja) 2017-06-22
EP3243756A4 (en) 2018-10-03
JP6271043B2 (ja) 2018-01-31
US20170369192A1 (en) 2017-12-28
US10752384B2 (en) 2020-08-25
EP3243756B1 (en) 2020-12-30

Similar Documents

Publication Publication Date Title
JP6271043B2 (ja) 軌道制御装置および衛星
EP3089267B1 (en) Method and system for controlling antenna of mobile communication application system based on double quaternions in mems inertial navigation
Xie et al. Highly constrained entry trajectory generation
CN104015938B (zh) 一种电推进静止轨道卫星的位置保持方法
CN105511490B (zh) 一种静止轨道卫星位置保持-角动量卸载联合控制方法
US7142981B2 (en) Laser range finder closed-loop pointing technology of relative navigation, attitude determination, pointing and tracking for spacecraft rendezvous
CN108490963B (zh) 全电推进卫星电推力器故障模式下的位置保持方法及***
CA2948860C (en) Orbit transfer method for a spacecraft using a continuous or quasi-continuous thrust and embedded driving system for implementing such a method
US9309010B2 (en) Methods and apparatus for controlling a plurality of satellites using node-synchronous eccentricity control
CN105373133B (zh) 一种同步轨道电推进位置保持与角动量卸载联合控制方法
CN110632935B (zh) 一种编队卫星绕飞自主控制方法
CN109539903A (zh) 一种固体运载火箭椭圆转移轨道迭代制导控制方法
JP7292132B2 (ja) 衛星制御装置、観測システム、観測方法、および観測プログラム
US20130105632A1 (en) Method and System for Controlling a Set of at Least Two Satellites Adapted to Provide a Service
CN112486196B (zh) 一种满足严格时间位置约束的飞行器快速轨迹优化方法
CN106094529A (zh) 编队任务多脉冲控制条件下的推力器在轨自主标定方法
CN111552312B (zh) 一种同步轨道卫星共位策略的生成方法及装置
Brady et al. ALHAT system architecture and operational concept
RU2004126645A (ru) Способ управления кластером находящихся на геостационарной орбите спутников (варианты)
Dell’Elce et al. Comparison between analytical and optimal control techniques in the differential drag based rendez-vous
WO2016125145A1 (en) Method and system for station keeping of geo satellites
Wood The evolution of deep space navigation: 1962-1989
CN106005482B (zh) 一种适用于导航倾斜轨道卫星的零偏持续天数确定方法
Yan et al. Formation maintenance and fuel balancing for satellites with implusive control
JP6707206B2 (ja) 宇宙機制御装置、宇宙機制御方法、およびプログラム

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 16735056

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2016568739

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 15540726

Country of ref document: US

REEP Request for entry into the european phase

Ref document number: 2016735056

Country of ref document: EP

NENP Non-entry into the national phase

Ref country code: DE