WO2019044735A1 - 宇宙機制御装置、宇宙機制御方法、およびプログラム - Google Patents

宇宙機制御装置、宇宙機制御方法、およびプログラム Download PDF

Info

Publication number
WO2019044735A1
WO2019044735A1 PCT/JP2018/031486 JP2018031486W WO2019044735A1 WO 2019044735 A1 WO2019044735 A1 WO 2019044735A1 JP 2018031486 W JP2018031486 W JP 2018031486W WO 2019044735 A1 WO2019044735 A1 WO 2019044735A1
Authority
WO
WIPO (PCT)
Prior art keywords
command value
thruster
spacecraft
angular momentum
injection
Prior art date
Application number
PCT/JP2018/031486
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 JP2019539474A priority Critical patent/JP6707206B2/ja
Publication of WO2019044735A1 publication Critical patent/WO2019044735A1/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/26—Guiding or controlling apparatus, e.g. for attitude control using jets

Definitions

  • the present invention relates to a spacecraft control device, a spacecraft control method, and a program for performing orbit holding control and angular momentum unloading of a spacecraft orbiting an orbit around a celestial body.
  • the spacecraft is equipped with an attitude control actuator that controls the attitude of the spacecraft, and a thruster that maintains the longitude and latitude of the spacecraft in a desired range.
  • the spacecraft is, for example, a geostationary satellite, an earth observation satellite, a spacecraft, and the like.
  • a reaction wheel, CMG (Control Moment Gyro) or the like is mounted on a spacecraft.
  • CMG Control Moment Gyro
  • Various disturbance torques such as solar radiation torque and magnetic torque work on the spacecraft. Therefore, it is necessary to rotate the reaction wheel to absorb disturbance torque and maintain the attitude of the spacecraft in a desired direction.
  • unloading of the angular momentum stored in the attitude control actuator is performed by generating unloading torque using a thruster.
  • Patent Documents 1 and 2 disclose techniques for performing orbit holding control for maintaining the longitude and latitude of a spacecraft within a desired range, and unloading of angular momentum.
  • the thruster mounted on the satellite disclosed in Patent Document 1 is fired near the descending node of the orbit and near the ascending node.
  • the control system disclosed in U.S. Pat. No. 5,677,698 uses the optimization of the cost function over a receding horizon subject to constraints on the spacecraft's attitude and input to the thruster so as to unload excess momentum. Generate control input commands that simultaneously control the thruster and spacecraft momentum exchange devices.
  • a propellant is necessary to perform the above-mentioned orbit holding control and angular momentum unloading.
  • the spacecraft can not perform orbit holding control and angular momentum unloading. Therefore, there is a need to reduce the amount of propellant used in track holding control and unloading of angular momentum.
  • the solar cell paddle of the spacecraft tends to be large.
  • the influence of the solar radiation pressure acting on the spacecraft increases.
  • the control amount of east-west control for maintaining the longitude of the satellite within the desired range is increased.
  • the solar radiation pressure torque acting on the spacecraft is also increased, the angular momentum of the attitude control actuator mounted on the spacecraft tends to be saturated.
  • trajectory holding control and angular momentum unloading are performed near the ascending node and the descending node, or in the vicinity of the descending node and the ascending node of the orbit, as disclosed in Patent Document 1, for example. And doing may not be optimal in terms of reducing propellant usage.
  • the control system disclosed in Patent Document 2 optimizes the operation of a spacecraft by using Model Predictive Control (MPC) over reverse horizon.
  • the control system determines the next state of the spacecraft based on the current state of the spacecraft, the current model of the spacecraft, and the current control inputs to the spacecraft. Therefore, the dimension of the optimization variable of MPC becomes a huge number from several tens to several hundreds, and processing becomes complicated, and there is a problem that it is not suitable for processing by a computer mounted on a spacecraft.
  • MPC Model Predictive Control
  • the present invention has been made in view of the above-described circumstances, and the amount of propellant used is small, and processing for performing track holding control and unloading of angular momentum is performed by a computer mounted on a space machine. It is an object of the present invention to provide a spacecraft control device, a spacecraft control method, and a program that are simplified as possible.
  • a spacecraft control device includes an angular momentum calculation unit, a command value calculation unit, a thruster control unit, and a gimbal control unit.
  • the angular momentum calculation unit calculates the angular momentum of the spacecraft orbiting the orbit around the celestial body.
  • the command value calculation unit is a thruster command value instructing to jet at least one of a plurality of thrusters attached to the assembly of the spacecraft via the gimbal mechanism, and the gimbal mechanism at the time of jetting the thruster. Calculate an angle command value that indicates the angle of.
  • the thruster control unit controls the thruster based on the thruster command value.
  • the gimbal control unit controls the gimbal mechanism based on the angle command value.
  • the command value calculation unit obtains the solution of the non-linear programming problem in which the sum of the injection amounts of the thrusters in the plurality of injection sections is an objective function, so that each thruster command of the plurality of injection sections Calculate the value and the angle command value.
  • the first constraint of the nonlinear programming problem is based on the first function and the second function. In the first function, thrusters are controlled based on the thruster command values using thruster command values and angle command values in a plurality of injection sections in the trajectory as arguments, and the gimbal mechanism is controlled based on the corner command values.
  • the change amount of the orbital element representing the spacecraft's orbit and the change amount of the angular momentum are output.
  • the second function takes, as arguments, an orbital element, a time change rate of the orbital element, and an angular momentum, and outputs a control quantity of the orbital element and a control quantity of the angular momentum while the spacecraft goes around the orbit.
  • the first constraint is that the change amount of the track element matches the control amount of the track element, and the change amount of the angular momentum matches the control amount of the angular momentum.
  • a solution of the non-linear programming problem whose objective function is the sum of the injection quantities of the thrusters in the injection section is determined.
  • the amount of propellant used is small, and the processing for performing track holding control and angular momentum unloading is simplified.
  • Block diagram showing a configuration example of the spacecraft control device according to the first embodiment Block diagram showing a configuration example of a command value calculation unit according to Embodiment 1
  • the flowchart which shows an example of operation of command value calculation which the command value calculation part concerning Embodiment 1 performs
  • the flowchart which shows one example of operation of calculation of the 2nd function which the 2nd function calculation section which relates to the execution form 1 does
  • the flowchart which shows an example of operation of command value calculation which the command value calculation part concerning Embodiment 2 performs
  • Configuration diagram showing an example of the hardware configuration of the spacecraft control device according to the embodiment The block diagram which shows another example of the hardware constitutions of the spacecraft control apparatus which concerns on embodiment
  • Embodiment 1 The spacecraft control device according to the first embodiment will be described by way of an example of a spacecraft control device mounted on a geostationary satellite orbiting the earth, which is an example of a spacecraft.
  • the spacecraft control device performs orbit holding control for maintaining the geostationary satellite on a target orbit, and unloading of the angular momentum accumulated in the attitude control actuator for controlling the attitude of the geostationary satellite.
  • the geostationary satellite 10 shown in FIG. 1 includes a structure 5, gimbal mechanisms 31, 32, 33, 34 attached to the structure 5, and thrusters 21, 22, 23, 23 attached to gimbal mechanisms 31, 32, 33, 34, respectively. And 24.
  • the gimbal mechanisms 31, 32, 33, and 34 have the same structure, and are attached to the surface 51 of the outer surface of the assembly 5 opposite to the surface 50 facing the earth.
  • the thrusters 21, 22, 23, 24 are mounted on the surface 51 in different injection directions.
  • the thrusters 21, 22, 23, 24 are also referred to as NW thruster 21, NE thruster 22, SW thruster 23, and SE thruster 24, respectively.
  • the thrusters 21, 22, 23, 24 and the gimbal mechanisms 31, 32, 33, 34 are controlled by a spacecraft controller to be described later.
  • the spacecraft control device controls the angles of the gimbal mechanisms 31, 32, 33, 34 shown in FIG. 1 and is mounted on the assembly 5 through the gimbal mechanisms 31, 32, 33, 34.
  • 23 and 24 perform orbit holding control of the geostationary satellite 10 and unloading of angular momentum.
  • the x H axis, the y H axis, and the z H axis are bases, and the center of mass C.I of the geostationary satellite 10 is used.
  • M. Set an orbital coordinate system whose origin is (Center of Mass) and refer to it as appropriate.
  • x H axis is the center of the earth 4 which is an astronomical object located at the center of the orbit and the center of mass of the geosynchronous satellite 10 C.
  • the z H axis is the normal direction of the orbital plane of the geostationary satellite 10.
  • the y H axis constitutes the coordinates of the right hand system with the x H axis and the z H axis.
  • the gimbal mechanisms 31, 32, 33, 34 have the same structure, and therefore the gimbal mechanism 31 will be described with reference to FIG.
  • z T axis is located parallel to the thrust axis of the thruster 21 attached to the gimbal mechanism 31.
  • the above-described thrusters 21, 22, 23, 24 and the gimbal mechanisms 31, 32, 33, 34 are controlled by the spacecraft controller 1 shown in FIG.
  • the spacecraft control device 1 controls the angles of the gimbal mechanisms 31, 32, 33, 34, and jets at least one of the thrusters 21, 22, 23, 24 to control the orbit holding control of the geostationary satellite 10 Do unloading of angular momentum.
  • the spacecraft control device 1 includes a trajectory determination unit 11 that calculates an orbital element of the geostationary satellite 10, an angular momentum calculation unit 12 that computes the angular momentum of the geostationary satellite 10, and thruster command values for the thrusters 21, 22, 23, 24.
  • the thruster control unit 14 which controls the thrusters 21, 22, 23, 24 based on the thruster command value, and the command value calculation unit 13 which calculates the angle command value for the gimbal mechanism 31, 32, 33, 34; And a gimbal control unit 15 for controlling the gimbal mechanisms 31, 32, 33, 34 based on the above.
  • the orbit determination unit 11 is, for example, GPS information acquired from a GPS (Global Positioning System) receiver mounted on the geostationary satellite 10, a range received from the ground station, a range rate, an azimuth angle of the geostationary satellite 10 Orbital elements of the geostationary satellite 10 are calculated from the information of elevation angle and the like.
  • the range indicates the distance from the ground station to the geostationary satellite 10.
  • the range rate indicates the rate of change of the distance from the ground station to the geostationary satellite 10.
  • the ground station can obtain the azimuth angle and elevation angle of the geostationary satellite 10 from observation using an optical camera on the ground.
  • the orbit determination unit 11 removes periodic fluctuation components from the instantaneous values of the orbit elements of the geostationary satellite 10 to calculate an average orbit element.
  • the orbital elements include, for example, eccentricity, inclination angle and the like.
  • an average orbit element is used as an orbit element of the geostationary satellite 10.
  • the periodic fluctuation component is a periodic fluctuation component that is generated by a perturbation of the earth's gravity, such as tidal force of a perturbation object such as the sun or the moon, and the period is equal to or less than a threshold.
  • the trajectory determination unit 11 sends the calculated average trajectory element to the command value calculation unit 13.
  • the orbit determination unit 11 further calculates the time change rate of the average orbit element, and sends the time change rate of the average orbit element to the command value calculation unit 13.
  • the time rate of change of the mean orbital element indicates the amount of change of the mean orbital element in a sufficiently short time.
  • the angular momentum calculation unit 12 uses the attitude information of the geostationary satellite 10 obtained from sensors such as a star sensor mounted on the geostationary satellite 10, a gyro sensor, a magnetic sensor, an earth sensor, a sun sensor, etc. The attitude angular velocity is determined, and the angular momentum of the geostationary satellite 10 is calculated. The angular momentum calculation unit 12 sends the calculated angular momentum to the command value calculation unit 13.
  • the thruster command value is a command value for at least one of the thrusters 21, 22, 23, 24.
  • the thruster command value includes, for example, the injection timing, the injection phase, the sum of the injection amounts of the thrusters 21, 22, 23, 24 and the ratio of the injection amount to the sum of the injection amounts.
  • the angle command value indicates the gimbal angle ⁇ , ⁇ of at least one of the gimbal mechanisms 31, 32, 33, 34 at the time of injection of the thrusters 21, 22, 23, 24, respectively.
  • the thruster control unit 14 controls the thrusters 21, 22, 23, 24 based on the thruster command value calculated by the command value calculation unit 13.
  • the gimbal control unit 15 controls the gimbal mechanisms 31, 32, 33, 34 based on the angle command value calculated by the command value calculation unit 13.
  • the command value calculation unit 13 calculates a first function calculation unit 131 for obtaining a first function described later, a second function calculation unit 132 for obtaining a second function described later, and a first function
  • a first constraint condition setting unit 133 for determining a first constraint condition in which the output of the second function and the output of the second function coincide, and the second one in which the angles of the gimbal mechanisms 31, 32, 33, 34 are within a defined range
  • the sum of the injection amounts of the thrusters 21, 22, 23, 24 in the plurality of injection sections is an objective function.
  • the command value calculation unit 13 holds in advance the latitude and longitude tolerance of the geostationary satellite 10, the angular momentum tolerance, and the gimbal angles ⁇ and ⁇ tolerances.
  • the first function calculation unit 131 obtains a first function including the following equations (13) to (16) described later, using thruster command values and angle command values in a plurality of injection sections as arguments.
  • thrusters 21, 22, 23, 24 are controlled based on thruster command values which are arguments
  • gimbal mechanisms 31, 32, 33, 34 are controlled based on angular command values which are arguments.
  • the variation of the orbital element and the variation of the angular momentum are output.
  • the second function calculation unit 132 acquires the mean trajectory element and the time change rate of the mean trajectory element from the trajectory determination unit 11, and acquires the angular momentum from the angular momentum calculation unit 12. Then, the second function calculation unit 132 sets the average orbit element, the time change rate of the average orbit element, and the angular momentum of the geostationary satellite 10 as arguments to the following (17), (20), (22), 24) Find a second function including the equation. The second function outputs the control amount of the average orbit element of the geostationary satellite 10 and the control amount of angular momentum of the geostationary satellite 10.
  • the first constraint condition setting unit 133 obtains the first function from the first function calculation unit 131 and obtains the second function from the second function calculation unit 132. Then, the first constraint condition setting unit 133 is a conditional expression indicating that the change amount of the track element coincides with the control amount of the average track element and the change amount of the angular momentum coincides with the control amount of the angular momentum. Find the first constraint that is
  • the second constraint condition setting unit 134 obtains a second constraint condition that is a conditional expression indicating that the angles of the gimbal mechanisms 31, 32, 33, and 34 are within a defined range.
  • the solution solution unit 135 uses an objective function for a nonlinear programming problem in which the sum of the injection amounts of the thrusters 21, 22, 23, 24 in the plurality of injection sections is an objective function. Find a solution that minimizes As a result, thruster command values and angle command values of the plurality of injection sections are calculated.
  • the command value calculation unit 13 holds in advance an allowable range of deviation from the target position of the geostationary satellite 10. Specifically, the target value of the latitude of the geostationary satellite 10? Ref and the holding accuracy?? Max , and the target value of the longitude of the geostationary satellite 10? Ref and the holding accuracy?? Max are held in advance.
  • the latitude target value ⁇ ref and the longitude target value ⁇ ref indicate the target trajectory described above.
  • the latitude holding accuracy ⁇ max and the longitude holding accuracy ⁇ max indicate a defined range including the target orbit in which the geostationary satellite 10 is held.
  • the range of the latitude ⁇ of the geostationary satellite 10 is expressed by the following equation (1).
  • the range of the longitude ⁇ of the geostationary satellite 10 is expressed by the following equation (2). ⁇ ref ⁇ max ⁇ ⁇ ⁇ ⁇ ref + ⁇ max (1) ⁇ ref ⁇ max ⁇ ⁇ ⁇ ⁇ ref + ⁇ max (2)
  • ⁇ ref is 0, and when the observation target of the geostationary satellite 10 is Japan, ⁇ ref is, for example, 140 °.
  • the latitude holding accuracy ⁇ max and the longitude holding accuracy ⁇ max are, eg, 0.1 °.
  • the latitude range of the geostationary satellite 10 is -0.1 ° or more and 0.1 ° or less, and the longitude range of the geostationary satellite 10 is 139.9 ° or more and 140.1 ° or less It is.
  • the command value calculation unit 13 holds in advance an allowable range of deviation of the angular momentum of the geostationary satellite 10 from the target value. Specifically, the command value calculation unit 13 holds in advance the target values h xref , h yref and h zref of the angular momentum of the stationary satellite 10 in the inertial system, and the holding accuracies ⁇ h xmax , ⁇ h ymax and ⁇ h zmax of the angular momentum. ing.
  • the range of angular momentum of the geostationary satellite 10 is expressed by the following equations (3) to (5).
  • h xref , h yref and h zref are all zero.
  • h xref , h yref and h zref are all not zero.
  • ⁇ h xmax , ⁇ h ymax , and ⁇ h zmax are determined according to the magnitude of angular momentum that can be compensated by an attitude control actuator mounted on a stationary satellite 10 such as a reaction wheel or CMG (Control Moment Gyro).
  • ⁇ h xmax , ⁇ h ymax and ⁇ h zmax be, for example, 30 Nms.
  • the command value calculation unit 13 holds in advance the allowable range of the drive angle of the gimbal mechanism. Specifically, the command value calculation unit 13 holds in advance the upper limits ⁇ max and ⁇ max and the lower limits ⁇ max and ⁇ max of the drive angles of the gimbal mechanisms 31, 32, 33 and 34.
  • the upper limits ⁇ max and ⁇ max and the lower limits - ⁇ max and - ⁇ max of the drive angles of the gimbal mechanisms 31, 32, 33 and 34 are the above-mentioned sensors for obtaining attitude information of the geostationary satellite 10, attitude control actuators, attitude It is determined in accordance with the magnitude of the disturbance torque that can be compensated by the attitude control system including a control device that controls the control actuator.
  • the magnitude T of the torque is proportional to the square norm of the gimbal angle, the magnitude F of the thrust of the thrusters 21, 22, 23, 24, and the mounting distance r, as expressed by the following equation (6).
  • the thrust magnitude F of the thrusters 21, 22, 23, 24 is considered to be constant on the orbit. Therefore, the adjustment of the magnitude T of the torque is performed by the gimbal angles ⁇ and ⁇ . By increasing the gimbal angles ⁇ and ⁇ , it is possible to give the geostationary satellite 10 a larger torque. By increasing the gimbal angles ⁇ and ⁇ , the ability to unload the angular momentum stored in the geostationary satellite 10 is improved.
  • the torque acting on the geostationary satellite 10 by the injection of the thrusters 21, 22, 23, 24 is a disturbance torque for the attitude control system of the geostationary satellite 10.
  • the gimbal angles ⁇ and ⁇ are determined according to the upper limit value T max of the disturbance torque that can be compensated by the posture control system.
  • the maximum value ⁇ max of the gimbal angle ⁇ is expressed by the following equation (7).
  • the maximum value ⁇ max of the gimbal angle ⁇ is expressed by the following equation (8). In the following formulas (7) and (8), ⁇ is a positive number less than one.
  • the command value calculation unit 13 holds in advance the ranges of latitude and longitude, the range of angular momentum, and the range of the gimbal angles ⁇ and ⁇ of the geostationary satellite 10 described above.
  • the command value calculation unit 13 calculates a thruster command value and an angle command value by obtaining a solution of the non-linear programming problem based on these tolerances.
  • the command value calculation unit 13 obtains a solution that minimizes the objective function under a first constraint condition described later.
  • the first function calculation unit 131 obtains a first function using thruster command values and corner command values in a plurality of injection sections as arguments (step S11).
  • the second function calculating unit 132 acquires an average trajectory element from the trajectory determining unit 11, and acquires an angular momentum from the angular momentum calculating unit 12 (step S12).
  • the second function calculation unit 132 obtains a second function using the average orbit element, the time change rate of the average orbit element, and the angular momentum of the geostationary satellite 10 as arguments (step S13).
  • the first constraint condition setting unit 133 is a conditional expression that indicates that the change amount of the track element matches the control amount of the average track element, and the change amount of the angular momentum matches the control amount of the angular momentum.
  • a first constraint condition is obtained (step S14).
  • the second constraint condition setting unit 134 obtains a second constraint condition that is a conditional expression indicating that the angles of the gimbal mechanisms 31, 32, 33, 34 are within the defined range (step S15).
  • the solution unit 135 obtains an objective function (step S16).
  • the solution unit 135 generates an objective function for a nonlinear programming problem based on an objective function that is the sum of the first constraint, the second constraint, and the injection amounts of the thrusters 21, 22, 23, 24 in a plurality of injection sections.
  • the thruster command value and the corner command value are calculated by obtaining a solution to be minimized (step S17).
  • the command value calculation unit 13 outputs the thruster command value and the angle command value to the thruster control unit 14 and the gimbal control unit 15, respectively (step S18).
  • the first function calculation unit 131 obtains a first function using thruster command values and corner command values in a plurality of injection sections as arguments.
  • the case where N injection sections are provided per circumference of the track will be described as an example.
  • An example will be described in which the thruster command value is an instruction to fire two of the thrusters 21, 22, 23, 24 in each of the N injection intervals.
  • the combination of the two is either a set of NW thruster 21 and NE thruster 22, a set of SW thruster 23 and SE thruster 24, a set of SW thruster 23 and NW thruster 21, and a set of NE thruster 22 and SE thruster 24. It is.
  • the amount of propellant used can be reduced by injecting only two of the thrusters 21, 22, 23, 24 in one injection section.
  • a j and b j Two of the thrusters 21, 22, 23, 24 to be injected in the j (1 ⁇ j ⁇ N) injection interval are denoted by a j and b j .
  • a j and b j are any of a pair of NW thruster 21 and NE thruster 22, a pair of SW thruster 23 and SE thruster 24, a pair of SW thruster 23 and NW thruster 21, and a pair of NE thruster 22 and SE thruster 24 It is.
  • the thruster a j for injecting the j-th injection interval, when the total injection amount of b j and F j, using the distribution coefficient xi] j to determine the respective injection quantity a j, b j, thrusters a j,
  • the injection amounts f aj and f bj of b j are expressed by the following equations (11) and (12).
  • the distribution coefficient ⁇ j is a coefficient that satisfies the condition of ⁇ 1 ⁇ ⁇ j ⁇ 1.
  • the parameters of the first function, the j-th sum F j of the injection amount in the injection interval, the mean longitude lambda Avgj in the j-th allocation in the injection interval coefficients xi] j, j-th injection interval, and j th injection section The angles ⁇ aj , ⁇ aj , ⁇ bj , ⁇ bj are given.
  • the angles ⁇ aj and ⁇ aj are angles in the j-th injection section of the gimbal mechanisms 31, 32, 33 and 34 for attaching the thrusters a j to the assembly 5.
  • the angles ⁇ bj and ⁇ bj are angles in the j-th injection section of the gimbal mechanisms 31, 32, 33 and 34 for attaching the thrusters b j to the assembly 5.
  • Variation of orbital elements is the output of the first function is based on the linearized Gaussian planet equations, the sum F j of the injection quantity, distribution coefficient xi] j, the mean longitude lambda Avgj, and angle alpha aj, beta It is represented by aj , ⁇ bj , ⁇ bj .
  • the variation of the orbital element includes the variation of the mean eccentricity vector, the mean inclination angle vector, and the mean nadir longitude.
  • the amount of change in angular momentum which is the output of the first function, is the sum of injection amount F j , distribution coefficient ⁇ j , average longitude ⁇ avgj , and angle ⁇ aj , ⁇ based on the linearized Euler equation It is represented by aj , ⁇ bj , ⁇ bj .
  • G i (F j , ⁇ j , ⁇ avgj , ⁇ aj , ⁇ aj , ⁇ bj , ⁇ bj ) is a Gauss's planetary equation for the inclination angle vector.
  • G ⁇ (F j , ⁇ j , ⁇ avgj , ⁇ aj , ⁇ aj , ⁇ bj , ⁇ bj ) is a Gauss's planetary equation for the mean nadir longitude.
  • E h (F j , ⁇ j , ⁇ avgj , ⁇ aj , ⁇ aj , ⁇ bj , ⁇ bj ) is Euler's equation relating to the angular momentum of the geostationary satellite 10.
  • the first function calculator 131 determines the above parameters, obtains a first function including the following equations (13) to (16), and sends the first function to the first constraint setting unit 133.
  • step S ⁇ b> 12 the second function calculator 132 acquires the mean trajectory element and the time change rate of the mean trajectory element from the trajectory determination unit 11, and acquires the angular momentum from the angular momentum calculator 12. Then, in step S13, the second function calculating unit 132 determines parameters using the average orbit element, the time change rate of the average orbit element, and the angular momentum of the geostationary satellite 10 as arguments, and the average orbit element of the geostationary satellite 10 A second function is obtained which outputs the control amount of x and the control amount of the angular momentum of the geostationary satellite 10.
  • the control amount of the average orbit element is the control amount of the average orbit element to be achieved by injecting the thrusters 21, 22, 23, 24.
  • the mean trajectory element to be controlled includes a mean eccentricity vector, a mean tilt angle vector, and a mean nadir longitude.
  • the second function calculator 132 determines the radius e R of a circle centered on the origin of the orbital coordinate system that determines the target value e target of the average eccentricity vector (step S21).
  • the second function calculation unit 132 sets the target value e target of the average eccentricity vector set at the intersection of a circle connecting a radius e R centered on the origin of the orbital coordinate system and the line connecting the earth 4 and the sun It calculates a control amount .DELTA.e c to approximate eccentricity vector (step S22).
  • the second function calculating unit 132 is configured to control one of the annual average inclination angle vector, the monthly average inclination angle vector, and the daily average inclination angle vector according to the magnitude of the latitude holding accuracy ⁇ max as an average inclination to be controlled. It is determined as an angle vector (step S23).
  • the second function calculating unit 132 calculates a control amount .DELTA.i c to the average tilt angle vector determined in step S23 (step S24).
  • the second function calculating unit 132 calculates the control amount ⁇ c of the average straight point longitude (step S25).
  • the second function calculator 132 calculates a control amount ⁇ h c of angular momentum (step S26). As shown in FIG. 7, the second function calculator 132 can perform the processes of steps S21, S23, S25, and S26 in parallel.
  • the second function calculation unit 132 sets the target value e target of the average eccentricity vector at the intersection of the line connecting the circle of radius e R centered on the origin of the orbital coordinate system and the earth 4 and the sun (17) calculates the control amount .DELTA.e c average eccentricity vector formula.
  • e OD is an average eccentricity vector calculated by the trajectory determination unit 11.
  • K e is a positive coefficient.
  • the second term -K e (e OD -e target ) of the following equation (17) corresponds to the feedback control amount.
  • the first term ⁇ e FF in the following equation (17) corresponds to the feedforward control amount.
  • the feedforward control amount compensates for the fluctuation of the average eccentricity vector due to a perturbation force acting on the geostationary satellite 10, for example, solar radiation pressure, moon tide force, sun tide force, and the like.
  • the ⁇ e FF found in the following equation (17) is expressed by the following equation (18).
  • e ′ OD is a time change rate of e OD
  • T c is a cycle in which the geostationary satellite 10 orbits.
  • ⁇ e c ⁇ e FF ⁇ K e (e OD ⁇ e target ) (17)
  • ⁇ e FF e ' OD T c (18)
  • the target value e target of the average eccentricity vector in the above equation (17) is set at the intersection of a circle of radius e R centered on the origin of the orbital coordinate system and a line connecting the earth 4 and the sun.
  • e R is expressed by the following equation (19).
  • ⁇ M is an error of the average trajectory element calculated by the trajectory determination unit 11, an error of the longitude caused by an injection error of the thrusters 21, 22, 23, 24, and the like.
  • ⁇ M is a value that satisfies 0 ⁇ M ⁇ max and is, for example, 0.02 °.
  • the second function calculator 132 controls an average inclination angle vector, for example, an annual average inclination angle vector, a monthly average inclination angle vector, and a day average inclination angle vector.
  • the second function calculator 132 controls an average tilt angle vector suitable for the magnitude of the latitude holding accuracy ⁇ max . For example, for 0.040 ° ⁇ (1 + ⁇ ⁇ ) ⁇ ⁇ max, the second function calculating unit 132 determines the annual average tilt angle vector as a control target.
  • the second function calculating unit 132 determines a control target of the monthly average tilt angle vector. In the case of ⁇ max ⁇ 0.004 ° ⁇ (1 + ⁇ ⁇ ), the second function calculator 132 determines a day average inclination angle vector as a control target.
  • the kappa phi is a positive coefficient, for example, is set to 1.1.
  • Control amount .DELTA.i c of average inclination angle vector is expressed by the following equation (20).
  • i OD is an average inclination angle vector calculated by the trajectory determination unit 11.
  • i target is a target value of a defined average tilt angle vector.
  • K i is a positive constant.
  • the second term -K i (i OD -i target ) in the following equation (20) corresponds to a feedback control amount.
  • the first term ⁇ i FF in the following equation (20) corresponds to the feedforward control amount.
  • the feedforward control amount compensates for the fluctuation of the average tilt angle vector due to the perturbation force acting on the geostationary satellite 10.
  • ⁇ i FF is expressed by the following equation (21).
  • di FF / dt is the time change rate of the average tilt angle vector.
  • ⁇ i c ⁇ i FF ⁇ K i (i OD ⁇ i target ) (20)
  • the second function calculation unit 132 calculates the control amount ⁇ c of the average straight point longitude represented by the following equation (22).
  • ⁇ OD is the average nadir longitude calculated by the trajectory determination unit 11.
  • ⁇ target is a target value of the average direct longitude point.
  • K ⁇ is a positive coefficient.
  • the second term ⁇ K ⁇ ( ⁇ OD ⁇ target ) in the following equation (22) corresponds to a feedback control amount.
  • the first term ⁇ EF in the following equation (22) corresponds to the feedforward control amount.
  • ⁇ EF is expressed by the following equation (23).
  • ⁇ ′ OD is the angular velocity of the average direct longitude point.
  • ⁇ ′ ′ OD is the angular acceleration of the average straight point longitude.
  • ⁇ c ⁇ FF ⁇ K ⁇ ( ⁇ OD ⁇ target ) (22)
  • the second function calculator 132 calculates a control amount ⁇ h c of angular momentum represented by the following equation (24).
  • h AD is the angular momentum of the geostationary satellite 10 in the inertial system calculated by the angular momentum calculation unit 12.
  • the second term -K h (h AD -h ref ) of the following equation (24) corresponds to the feedback control amount.
  • Kh is a positive coefficient.
  • the first term ⁇ h FF in the following equation (24) corresponds to the feedforward control amount.
  • the feedforward control compensates for the amount of change in angular momentum in the inertial system caused by the perturbation torque.
  • the most dominant torque among the perturbation torques is the torque due to the solar radiation pressure.
  • ⁇ SRP the solar radiation pressure torque acting on the geostationary satellite 10
  • ⁇ h FF is approximated by the following equation (25).
  • h ref is a target value of angular momentum in the inertial system, and is expressed by the following equation (26).
  • Each component of the target value h ref of the angular momentum is held in advance by the command value calculation unit 13.
  • the second function calculation unit 132 sends the second function including the expressions (17), (20), (22), and (24) to the first constraint condition setting unit 133.
  • the first constraint condition setting unit 133 is a conditional expression that indicates that the change amount of the track element matches the control amount of the average track element, and the change amount of the angular momentum matches the control amount of the angular momentum. Determine the first constraint.
  • the second constraint condition setting unit 134 obtains a second constraint condition that is a conditional expression indicating that the angles of the gimbal mechanisms 31, 32, 33, and 34 are within a defined range.
  • the second constraint condition setting unit 134 sends the second constraint condition including the following equations (31) to (34) to the solution unit 135.
  • ⁇ max ⁇ ⁇ aj ⁇ ⁇ max (31) ⁇ max ⁇ ⁇ aj ⁇ ⁇ max (32) ⁇ max ⁇ ⁇ bj ⁇ ⁇ max (33) ⁇ max ⁇ ⁇ bj ⁇ ⁇ max (34)
  • the solution unit 135 generates an objective function for a nonlinear programming problem based on an objective function that is the sum of the first constraint, the second constraint, and the injection amounts of the thrusters 21, 22, 23, 24 in a plurality of injection sections. By obtaining a solution to be minimized, thruster command values and angle command values of each of the plurality of injection sections are calculated.
  • the objective function is expressed by the following equation (37).
  • the solution solving unit 135 obtains a solution of the total injection amount F j , the distribution coefficient ⁇ j , the average longitude ⁇ avgj , and the angles ⁇ aj , ⁇ aj , ⁇ bj , and ⁇ bj in each of the injection sections.
  • methods such as sequential quadratic programming and interior point method are used.
  • the solution solution unit 135 obtains the injection timing including the start position and the end position of the injection section of the thrusters 21, 22, 23, 24 from the solution of ⁇ avgj obtained by solving the non-linear programming problem.
  • the solution unit 135 obtains the injection amount of the thrusters 21, 22, 23, 24 in the injection section from the solution of F j and ⁇ j obtained by solving the non-linear programming problem. Further, the solution unit 135 obtains gimbal angles at the time of jetting the thrusters 21, 22, 23, 24 from the solutions of ⁇ aj , ⁇ aj , ⁇ bj , ⁇ bj obtained by solving the non-linear programming problem.
  • the solution solution unit 135 sends thruster command values indicating the injection timing and injection amount of the thrusters 21, 22, 23, 24 to the thruster control unit 14, and sends an angle command value indicating the gimbal angle to the gimbal control unit 15.
  • the thruster control unit 14 controls and jets the thrusters 21, 22, 23, 24 based on the thruster command value.
  • the gimbal control unit 15 controls the gimbal angles of the gimbal mechanisms 31, 32, 33, 34 based on the angle command values. Solve the non-linear programming problem as mentioned above and calculate thruster command value and angle command value, and simplify the process for carrying out track maintenance control and unloading of angular momentum while reducing the amount of propellant used. It is possible to
  • the spacecraft control device 1 finds a solution that minimizes the objective function which is the sum of the injection amounts of the thrusters in the non-linear programming problem, and obtains thruster command values and angle commands. By calculating the value, it is possible to simplify the process for performing the track holding control and the unloading of the angular momentum while suppressing the amount of propellant used.
  • the thruster 21 when it is necessary to fire the thrusters 21, 22, 23, 24 at a position separated from the ascending intersection and the descending intersection, the thruster 21 is also provided during one orbit. , 22, 23, 24 even when it is necessary to inject a plurality of times, it is possible to carry out a process for carrying out track holding control and unloading of angular momentum while suppressing the amount of propellant used.
  • the constraints of the nonlinear programming problem in the first embodiment are the first constraint and the second constraint
  • the nonlinear programming problem may further include other constraints.
  • the configuration of the spacecraft control device 1 according to the second embodiment is the same as the configuration of the spacecraft control device 1 according to the first embodiment shown in FIG.
  • FIG. 8 is a diagram showing an example of a configuration of a command value calculation unit according to Embodiment 2 of the present invention.
  • the command value calculation unit 16 according to the second embodiment further includes a third constraint condition setting unit 136 in addition to the configuration of the command value calculation unit 13 according to the first embodiment shown in FIG.
  • the third constraint condition setting unit 136 sets a third constraint condition, which is a conditional expression indicating that the angular momentum is within a defined range, at least either immediately before and / or immediately after each of the plurality of injection sections. Then, the solution is sent to the solution unit 135.
  • the command value calculation unit 16 calculates the thruster command value and the angle command value.
  • the processes of steps S11 to S16 and S18 are the same as the processes shown in FIG.
  • the third constraint condition setting unit 136 determines that the angular momentum is determined within at least one of immediately before and / or immediately after injection of each of the plurality of injection sections.
  • a third constraint condition which is a conditional expression indicating that there is, is obtained (step S19).
  • the solution unit 135 obtains an objective function (step S16).
  • the solution solving unit 135 is a nonlinear program based on an objective function that is the sum of the first constraint, the second constraint, the third constraint, and the injection amounts of the thrusters 21, 22, 23, 24 in a plurality of injection sections.
  • a thruster command value and an angle command value are calculated by obtaining a solution that minimizes the objective function (step S20).
  • the subsequent processing is the same as the processing shown in FIG.
  • the third constraint condition setting unit 136 is a third constraint that is a conditional expression indicating that the angular momentum is within a defined range at least either immediately before or after each of the plurality of injection sections. Determine the condition.
  • the third constraint condition setting unit 136 obtains a third constraint condition that is a conditional expression indicating that at least one of the following equations (38) and (39) holds.
  • H j ( ⁇ ) in the following equation (38) is the angular momentum immediately before the injection of the j-th injection section.
  • the angular momentum immediately before the injection of the j-th injection interval is the angular momentum of the geostationary satellite 10 immediately before the start of the injection of the thrusters 21, 22, 23, 24 in the j-th injection interval.
  • H j (+) in the following equation (39) is the angular momentum immediately after the injection of the j-th injection section.
  • the angular momentum immediately after the j-th injection interval is the angular momentum of the geostationary satellite 10 immediately after the end of the injection of the thrusters 21, 22, 23, 24 in the j-th injection interval.
  • h j ( ⁇ ) and h j (+) are functions of variables F j , ⁇ j , ⁇ avgj , ⁇ aj , ⁇ aj , ⁇ bj , ⁇ bj .
  • the angular momentum is within the defined range immediately before and / or immediately after each of the plurality of injection sections.
  • the angular momentum of the geostationary satellite 10 can be maintained within the holding range by making the third constraint condition of the nonlinear programming problem. Therefore, it is possible to carry out track holding control and unloading of the angular momentum without saturating the angular momentum of the posture control actuator.
  • FIG. 10 is a diagram illustrating an example of a hardware configuration.
  • FIG. 10 shows an example in which the spacecraft controller 1 is realized by the processing circuit 101.
  • the geostationary satellite 10 includes a sensor 110 for observing the position, velocity, attitude angle, and angular velocity of the satellite, the spacecraft control device 1, thrusters 21, 22, 23, 24, and gimbal mechanisms 31, 32, 33. , 34.
  • the processing circuit 101 may be, for example, a single circuit, a composite circuit, a programmed processor, a parallel programmed processor, an application specific integrated circuit (ASIC), an FPGA (field programmable) Gate Array) or a combination thereof is applicable.
  • the processing circuit 101 may realize each of the functions of the trajectory determination unit 11, the angular momentum calculation unit 12, the command value calculation unit 13, the thruster control unit 14, and the gimbal control unit 15, or It may be realized by the processing circuit 101.
  • FIG. 11 is a diagram showing another example of the hardware configuration according to the present embodiment.
  • FIG. 11 shows an example in which the spacecraft control device 1 is realized by a processor (CPU: Central Processing Unit) 102 and a memory 103.
  • the processor 102 and the memory 103 in FIG. 11 correspond to the processing circuit 101 in FIG.
  • the functions of the trajectory determination unit 11, the angular momentum calculation unit 12, the command value calculation unit 13, the thruster control unit 14, and the gimbal control unit 15 are software, It is realized by firmware or a combination of software and firmware.
  • Software and firmware are described as a program and stored in the memory 103.
  • the processor 102 implements the functions of the respective units by reading and executing the program stored in the memory 103. That is, the memory 103 stores programs for executing the processing of the trajectory determination unit 11, the angular momentum calculation unit 12, the command value calculation unit 13, the thruster control unit 14, and the gimbal control unit 15.
  • the memory 103 is, for example, nonvolatile or random access memory (RAM), read only memory (ROM), flash memory, erasable programmable read only memory (EPROM), electrically erasable and programmable read only memory (EEPROM), or the like. It includes volatile semiconductor memory, magnetic disk, flexible disk, optical disk, compact disk, mini disk, DVD (Digital Versatile Disc) and the like.
  • RAM nonvolatile or random access memory
  • ROM read only memory
  • EPROM erasable programmable read only memory
  • EEPROM electrically erasable and programmable read only memory
  • the memory 103 also includes a trajectory determination value determined by the trajectory determination unit 11, a satellite angular momentum value set by the angular momentum calculation unit 12, a thruster command value calculated by the command value calculation unit 13, and a gimbal Store the corner command value.
  • the information stored in the memory 103 described above is appropriately read and used for processing of the trajectory determination unit 11, angular momentum calculation unit 12, command value calculation unit 13, thruster control unit 14, and gimbal control unit 15, and correction Be done.
  • the functions of the trajectory determination unit 11, the angular momentum calculation unit 12, the command value calculation unit 13, the thruster control unit 14, and the gimbal control unit 15 are partially realized by dedicated hardware, and partially performed by software or It may be realized by firmware.
  • the function of the trajectory determination unit 11 is realized by the processing circuit 101 as dedicated hardware, and the command value calculation unit 13 is implemented by reading and executing the program stored in the memory 103 by the processor 102. It is possible to realize the function.
  • Embodiments of the present invention are not limited to the above-described embodiments.
  • the spacecraft is not limited to the geostationary satellite 10, but is any flying object orbiting around the celestial body.
  • An example of a flying object is an earth observation satellite, a spacecraft, and the like.
  • the number of thrusters 21, 22, 23, 24 is arbitrary, and the number and combination of thrusters 21, 22, 23, 24 to be injected in each injection section are also arbitrary.
  • the degree of freedom of the gimbal mechanisms 31, 32, 33, 34 may be one.
  • the first function is the variable F j , ⁇ j , ⁇ avgj , ⁇ aj , ⁇ bj It is a function.
  • h j ( ⁇ ) and h j (+) are functions of variables F j , ⁇ j , ⁇ avgj , ⁇ aj , ⁇ bj .
  • the thrusters 21, 22, 23, 24 are attached to the spacecraft structure 5, but the number of thrusters 21, 22, 23, 24 attached to the spacecraft structure 5 is arbitrary. .
  • the angle ⁇ and the angle ⁇ of each of the thrusters 21, 22, 23, 24 may be the same as or different from each other.
  • the latitude and longitude range, and the latitude and longitude retention accuracy are not limited to the above-described example, and can be arbitrarily determined.
  • the second function calculation unit 132 may calculate the time change ratio of the average trajectory element.
  • the trajectory determination unit 11 may acquire the activation component from the external device instead of calculating the trajectory component.
  • the trajectory determination unit 11 may send the instantaneous value of the trajectory element to the command value calculation unit 13.
  • the solution solution unit 135 may calculate the thruster command value and the gimbal command value by finding a solution of the nonlinear programming problem under the first constraint. Alternatively, the solution unit 135 may obtain the solution of the nonlinear programming problem under the first and third constraints. Further, the solution solution unit 135 may calculate the thrust command value and the gimbal command value by obtaining a semi-optimal solution that reduces the objective function under the first constraint condition.

Landscapes

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

Abstract

角運動量算出部（１２）は、宇宙機の角運動量を算出する。指令値算出部（１３）は、宇宙機の軌道を定める軌道要素の変化量と軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致する第１の制約条件の下で、複数の噴射区間におけるスラスタ（２１，２２，２３，２４）の噴射量の合計が目的関数である非線形計画問題の解を求めることで、スラスタ指令値および角指令値を算出する。スラスタ制御部（１４）は、スラスタ指令値に基づき、スラスタ（２１，２２，２３，２４）の噴射を実行する。ジンバル制御部（１５）は、角指令値に基づき、ジンバル機構（３１，３２，３３，３４）の角度を調節する。

Description

宇宙機制御装置、宇宙機制御方法、およびプログラム
　本発明は、天体の回りの軌道を周回する宇宙機の軌道保持制御および角運動量のアンローディングを行う宇宙機制御装置、宇宙機制御方法、およびプログラムに関する。
　宇宙機には、宇宙機の姿勢を制御する姿勢制御アクチュエータ、および宇宙機の経度および緯度を所望の範囲に維持するスラスタが搭載される。宇宙機は、例えば、静止衛星、地球観測衛星、宇宙船等である。姿勢制御アクチュエータの一例として、リアクションホイール、ＣＭＧ（Control Moment Gyro）等が宇宙機に搭載される。宇宙機には、太陽輻射トルク、磁気トルク等の種々の外乱トルクが働く。そこで、リアクションホイールを回転させて、外乱トルクを吸収して、宇宙機の姿勢を所望の向きに維持する必要がある。リアクションホイールの回転数には限度があるため、リアクションホイールの回転数の増加を抑え、リアクションホイールの回転数を正常範囲内に維持することが行われる。これを姿勢制御アクチュエータに蓄積された角運動量のアンローディングという。詳細には、スラスタを用いてアンローディングトルクを発生させることで、アンローディングが行われる。
　宇宙機の経度および緯度を所望の範囲に維持する軌道保持制御と、角運動量のアンローディングとを行う技術が特許文献１，２に開示されている。特許文献１に開示される衛星に搭載されたスラスタは、軌道の下降ノードの近傍、および上昇ノードの近傍で噴射される。特許文献２に開示される制御システムは、超過運動量をアンロードするように、宇宙船の姿勢およびスラスタへの入力に対する制約を条件とした後退ホライズンにわたるコスト関数の最適化を用いて、宇宙船のスラスタおよび宇宙船の運動量交換デバイスを同時に制御する制御入力コマンドを生成する。
特開平９－２７７９９７号公報 特開２０１６－１２４５３８号公報
　上述の軌道保持制御と角運動量のアンローディングとを行うためには、推薬が必要である。推薬が枯渇すると、宇宙機は軌道保持制御および角運動量のアンローディングを行うことができない。そのため、軌道保持制御と角運動量のアンローディングとにおける推薬の使用量の低減が求められている。
　一方で、宇宙機の通信能力を増大させるために宇宙機の発生電力の大電力化が求められ、宇宙機の太陽電池パドルが大型化する傾向にある。太陽電池パドルが大型化すると、宇宙機に働く太陽輻射圧の影響が増大する。その結果、太陽輻射圧の影響によって増大する衛星の軌道の離心率を補償するために、衛星の経度を所望の範囲内に維持する東西制御の制御量が増大する。また宇宙機に働く太陽輻射圧トルクも増大するため、宇宙機に搭載された姿勢制御アクチュエータの角運動量が飽和しやすくなる。この場合、例えば、昇交点近傍および降交点近傍の二箇所、または、特許文献１に開示されるように、軌道の下降ノードの近傍および上昇ノードの近傍において、軌道保持制御と角運動量のアンローディングとを行うことは、推薬の使用量の削減の観点から、最適ではないことがある。
　特許文献２に開示される制御システムは、後退ホライズンにわたるモデル予測制御（Model Predictive Control：ＭＰＣ）を用いて、宇宙船の動作の最適化を行う。この制御システムは、宇宙船の現在の状態、宇宙船の現在のモデル、および宇宙船への現在の制御入力に基づいて、宇宙船の次の状態を求める。そのため、ＭＰＣの最適化変数の次元が数十から数百と膨大な数となって、処理が複雑化するため、宇宙機に搭載される計算機での処理に適さないという課題がある。
　本発明は上述の事情に鑑みてなされたものであり、推薬の使用量が少なく、軌道保持制御と角運動量のアンローディングとを行うための処理が、宇宙機に搭載される計算機で処理が可能となるように簡易化された宇宙機制御装置、宇宙機制御方法、およびプログラムを提供することを目的とする。
　上記目的を達成するため、本発明に係る宇宙機制御装置は、角運動量算出部、指令値算出部、スラスタ制御部、およびジンバル制御部を備える。角運動量算出部は、天体の周りの軌道を周回する宇宙機が有する角運動量を算出する。指令値算出部は、それぞれが宇宙機の構体に、ジンバル機構を介して取り付けられた複数のスラスタの内、少なくともいずれかのスラスタの噴射を指示するスラスタ指令値、およびスラスタの噴射時におけるジンバル機構の角度を指示する角指令値を算出する。スラスタ制御部は、スラスタ指令値に基づき、スラスタを制御する。ジンバル制御部は、角指令値に基づき、ジンバル機構を制御する。指令値算出部は、第１の制約条件の下で、複数の噴射区間におけるスラスタの噴射量の合計が目的関数である非線形計画問題の解を求めることで、複数の噴射区間のそれぞれのスラスタ指令値および角指令値を算出する。非線形計画問題の第１の制約条件は、第１の関数および第２の関数に基づく。第１の関数は、軌道における複数の噴射区間でのスラスタ指令値および角指令値を引数として、該スラスタ指令値に基づいてスラスタが制御され、該角指令値に基づいてジンバル機構が制御された場合に、宇宙機が軌道を一周する間の、宇宙機の軌道を表す軌道要素の変化量と角運動量の変化量とを出力する。第２の関数は、軌道要素、軌道要素の時間変化率、および角運動量を引数として、宇宙機が軌道を一周する間の、軌道要素の制御量と角運動量の制御量とを出力する。第１の制約条件は、軌道要素の変化量と軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することである。
　本発明によれば、軌道要素の変化量と軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致する第１の制約条件の下で、複数の噴射区間におけるスラスタの噴射量の合計が目的関数である非線形計画問題の解が求められる。非線形計画問題の解を求めて、スラスタ指令値および角指令値を算出するため、推薬の使用量が少なく、軌道保持制御および角運動量のアンローディングを行うための処理が簡易化される。
本発明の実施の形態１における静止衛星の構成の例を示す図 実施の形態１におけるスラスタの宇宙機への取付例を示す図 実施の形態１におけるスラスタおよびジンバル機構の例を示す図 実施の形態１に係る宇宙機制御装置の構成例を示すブロック図 実施の形態１に係る指令値算出部の構成例を示すブロック図 実施の形態１に係る指令値算出部が行う指令値算出の動作の一例を示すフローチャート 実施の形態１に係る第２関数算出部が行う第２の関数の算出の動作の一例を示すフローチャート 本発明の実施の形態２に係る指令値算出部の構成例を示す図 実施の形態２に係る指令値算出部が行う指令値算出の動作の一例を示すフローチャート 実施の形態に係る宇宙機制御装置のハードウェア構成の一例を示す構成図 実施の形態に係る宇宙機制御装置のハードウェア構成の他の一例を示す構成図
　以下、本発明の実施の形態に係る宇宙機制御装置について図面を参照して詳細に説明する。なお図中、同一または同等の部分には同一の符号を付す。
　（実施の形態１）
　実施の形態１に係る宇宙機制御装置を、宇宙機の一例である地球の周りを周回する静止衛星に搭載される宇宙機制御装置を例に説明する。この宇宙機制御装置は、静止衛星を目標とする軌道上に維持する軌道保持制御と、静止衛星の姿勢を制御する姿勢制御アクチュエータに蓄積された角運動量のアンローディングとを行う。
　上述したように実施の形態１に係る宇宙機制御装置の制御対象である静止衛星を図１に示す。図１に示す静止衛星１０は、構体５と、構体５に取り付けられるジンバル機構３１，３２，３３，３４と、ジンバル機構３１，３２，３３，３４のそれぞれに取り付けられるスラスタ２１，２２，２３，２４と、を備える。ジンバル機構３１，３２，３３，３４は、いずれも同じ構造であり、構体５の外面の内、地球に向いている面５０の反対側に位置する面５１に取り付けられる。スラスタ２１，２２，２３，２４は、噴射方向が互いに異なる向きで面５１に取り付けられる。スラスタ２１，２２，２３，２４をそれぞれ、ＮＷスラスタ２１、ＮＥスラスタ２２、ＳＷスラスタ２３、およびＳＥスラスタ２４とも呼ぶこととする。スラスタ２１，２２，２３，２４およびジンバル機構３１，３２，３３，３４は、後述する宇宙機制御装置によって制御される。詳細には、宇宙機制御装置は、図１に示すジンバル機構３１，３２，３３，３４の角度を制御し、ジンバル機構３１，３２，３３，３４を介して構体５に取り付けられるスラスタ２１，２２，２３，２４を噴射させることで、静止衛星１０の軌道保持制御と角運動量のアンローディングとを行う。
　図２に示すように、理解を容易にするため、ｘＨ軸、ｙＨ軸、ｚＨ軸を基底とし、静止衛星１０の質量中心Ｃ．Ｍ．（Center of Mass）を原点とする軌道座標系を設定し、適宜参照する。ｘＨ軸は、軌道の中心に位置する天体である地球４の中心と静止衛星１０の質量中心Ｃ．Ｍ．とを通る。ｚＨ軸は、静止衛星１０の軌道面の法線方向である。ｙＨ軸は、ｘＨ軸およびｚＨ軸と右手系の座標を構成する。ジンバル機構３１，３２，３３，３４のそれぞれから質量中心Ｃ．Ｍ．までの距離を示す取付距離ｒは互いに同じである。ｘＨｙＨ平面に投影された質量中心Ｃ．Ｍ．からスラスタ２１，２２，２３，２４のそれぞれを指すベクトルと、ｘＨ軸とがなす鋭角である角度ψは互いに同じである。ｘＨｚＨ平面に投影された質量中心Ｃ．Ｍ．からスラスタ２１，２２，２３，２４のそれぞれを指すベクトルと、ｚＨ軸とがなす鋭角である角度γは互いに同じである。
　ジンバル機構３１，３２，３３，３４の構造はいずれも同じであるから、ジンバル機構３１について図３を用いて説明する。ジンバル機構３１の回転軸を、互いに直交するｘＴ軸およびｙＴ軸とする。ｘＴ軸およびｙＴ軸に直交する軸をｚＴ軸とする。ｚＴ軸は、ジンバル機構３１に取り付けられるスラスタ２１の推力軸に平行に位置する。ｘＴ軸周りのジンバル角をαとし、ｙＴ軸周りのジンバル角をβとする。ジンバル角α，βが０である場合、ｚＴ軸は静止衛星１０の質量中心Ｃ．Ｍ．を貫き、スラスタ２１を噴射しても静止衛星１０にトルクが働かない。ジンバル角α，βの少なくとも１つが０でない場合は、スラスタ２１の噴射時に静止衛星１０にトルクが働く。
　上述のスラスタ２１，２２，２３，２４およびジンバル機構３１，３２，３３，３４は、図４に示す宇宙機制御装置１によって制御される。宇宙機制御装置１は、ジンバル機構３１，３２，３３，３４の角度を制御し、スラスタ２１，２２，２３，２４の内、少なくとも１つを噴射させることで、静止衛星１０の軌道保持制御と角運動量のアンローディングとを行う。
　宇宙機制御装置１は、静止衛星１０の軌道要素を算出する軌道決定部１１と、静止衛星１０の角運動量を算出する角運動量算出部１２と、スラスタ２１，２２，２３，２４に対するスラスタ指令値およびジンバル機構３１，３２，３３，３４に対する角指令値を算出する指令値算出部１３と、スラスタ指令値に基づきスラスタ２１，２２，２３，２４を制御するスラスタ制御部１４と、角指令値に基づきジンバル機構３１，３２，３３，３４を制御するジンバル制御部１５とを備える。
　軌道決定部１１は、例えば、静止衛星１０に搭載されたＧＰＳ（Global Positioning System：全地球測位システム）受信機から取得したＧＰＳ情報、地上局から受信したレンジ、レンジレート、静止衛星１０の方位角および仰角の情報等から、静止衛星１０の軌道要素を算出する。なおレンジは、地上局から静止衛星１０までの距離を示す。レンジレートは、地上局から静止衛星１０までの距離の変化の割合を示す。地上局は、地上の光学カメラを用いた観測から静止衛星１０の方位角および仰角を得ることができる。軌道決定部１１は、静止衛星１０の軌道要素の瞬時値から周期的な変動成分を取り除いて平均軌道要素を算出する。軌道要素は、例えば、離心率、傾斜角等を含む。以下の説明において、静止衛星１０の軌道要素として平均軌道要素を用いる。周期的な変動成分とは、地球重力の摂動、例えば太陽、月等の摂動天体の潮汐力等によって発生し、周期が閾値以下の周期的な変動成分である。軌道決定部１１は、算出した平均軌道要素を指令値算出部１３に送る。軌道決定部１１はさらに、平均軌道要素の時間変化率を算出し、平均軌道要素の時間変化率を指令値算出部１３に送る。平均軌道要素の時間変化率は、十分に短い時間における平均軌道要素の変化量を示す。
　角運動量算出部１２は、静止衛星１０に搭載された恒星センサ、ジャイロセンサ、磁気センサ、地球センサ、太陽センサ等のセンサから得られる静止衛星１０の姿勢情報を用いて、静止衛星１０の姿勢および姿勢角速度を求め、静止衛星１０の角運動量を算出する。角運動量算出部１２は、算出した角運動量を指令値算出部１３に送る。
　指令値算出部１３は、詳細については後述するが、非線形計画問題の解を求めることで、静止衛星１０を目標とする軌道から定められた範囲内に保持するための、スラスタ指令値および角指令値を算出する。なおスラスタ指令値は、スラスタ２１，２２，２３，２４の内、少なくともいずれかのスラスタに対する指令値である。詳細にはスラスタ指令値は、例えば、噴射タイミング、噴射位相、スラスタ２１，２２，２３，２４の噴射量の合計、噴射量の合計に対する噴射量の割合等を含む。角指令値は、スラスタ２１，２２，２３，２４のそれぞれの噴射時におけるジンバル機構３１，３２，３３，３４の少なくともいずれかのジンバル角α，βを示す。
　スラスタ制御部１４は、指令値算出部１３が算出したスラスタ指令値に基づき、スラスタ２１，２２，２３，２４を制御する。ジンバル制御部１５は、指令値算出部１３が算出した角指令値に基づきジンバル機構３１，３２，３３，３４を制御する。
　指令値算出部１３は、図５に示すように、後述する第１の関数を求める第１関数算出部１３１と、後述する第２の関数を求める第２関数算出部１３２と、第１の関数の出力と第２の関数の出力とが一致する第１の制約条件を決める第１制約条件設定部１３３と、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にある第２の制約条件を決める第２制約条件設定部１３４と、第１の制約条件と第２の制約条件の下で、複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計が目的関数である非線形計画問題の解を求めることで、スラスタ指令値および角指令値を算出する求解部１３５と、を備える。なお指令値算出部１３は、静止衛星１０の緯度および経度の許容範囲、角運動量の許容範囲、ならびにジンバル角α，βの許容範囲を予め保持している。
　第１関数算出部１３１は、複数の噴射区間でのスラスタ指令値および角指令値を引数として、後述する下記（１３）－（１６）式を含む第１の関数を求める。第１の関数は、引数であるスラスタ指令値に基づいてスラスタ２１，２２，２３，２４が制御され、引数である角指令値に基づいてジンバル機構３１，３２，３３，３４が制御された場合に、静止衛星１０が軌道を一周する間の、軌道要素の変化量および角運動量の変化量を出力する。
　第２関数算出部１３２は、軌道決定部１１から平均軌道要素および平均軌道要素の時間変化率を取得し、角運動量算出部１２から角運動量を取得する。そして、第２関数算出部１３２は、平均軌道要素、平均軌道要素の時間変化率、および静止衛星１０が有する角運動量を引数として、後述する下記（１７），（２０），（２２），（２４）式を含む第２の関数を求める。第２の関数は、静止衛星１０の平均軌道要素の制御量および静止衛星１０が有する角運動量の制御量を出力する。
　第１制約条件設定部１３３は、第１の関数を第１関数算出部１３１から取得し、第２の関数を第２関数算出部１３２から取得する。そして、第１制約条件設定部１３３は、軌道要素の変化量と平均軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することを示す条件式である第１の制約条件を求める。
　第２制約条件設定部１３４は、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にあることを示す条件式である第２の制約条件を求める。
　求解部１３５は、第１の制約条件と第２の制約条件の下で、複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計が目的関数である非線形計画問題について、目的関数を最小化する解を求める。その結果、複数の噴射区間のそれぞれのスラスタ指令値および角指令値が算出される。
　指令値算出部１３が予め保持している許容範囲の詳細について説明する。指令値算出部１３は、静止衛星１０の目標位置からのずれの許容範囲を予め保持している。詳細には、静止衛星１０の緯度の目標値φｒｅｆおよび保持精度Δφｍａｘ、ならびに、静止衛星１０の経度の目標値λｒｅｆおよび保持精度Δλｍａｘを予め保持している。緯度の目標値φｒｅｆおよび経度の目標値λｒｅｆは、上述の目標とする軌道を示す。緯度の保持精度Δφｍａｘおよび経度の保持精度Δλｍａｘは、静止衛星１０が保持される、目標とする軌道を含む定められた範囲を示す。
　静止衛星１０の緯度φの範囲は下記（１）式で表される。静止衛星１０の経度λの範囲は、下記（２）式で表される。
　　φｒｅｆ－Δφｍａｘ≦φ≦φｒｅｆ＋Δφｍａｘ　　　（１）
　　λｒｅｆ－Δλｍａｘ≦λ≦λｒｅｆ＋Δλｍａｘ　　　（２）
　静止衛星１０の場合、φｒｅｆは０であり、静止衛星１０の観測対象が日本であれば、λｒｅｆは例えば１４０°である。緯度の保持精度Δφｍａｘおよび経度の保持精度Δλｍａｘを、例えば０．１°とする。この場合、静止衛星１０の緯度の範囲は－０．１°以上、かつ、０．１°以下であり、静止衛星１０の経度の範囲は、１３９．９°以上、かつ、１４０．１°以下である。
　さらに、指令値算出部１３は、静止衛星１０の角運動量の目標値からのずれの許容範囲を予め保持している。詳細には、指令値算出部１３は、慣性系における静止衛星１０の角運動量の目標値ｈｘｒｅｆ，ｈｙｒｅｆ，ｈｚｒｅｆ、および角運動量の保持精度Δｈｘｍａｘ，Δｈｙｍａｘ，Δｈｚｍａｘを予め保持している。
　静止衛星１０の角運動量の範囲は、下記（３）－（５）式で表される。
　　ｈｘｒｅｆ－Δｈｘｍａｘ≦ｈｘ≦ｈｘｒｅｆ＋Δｈｘｍａｘ　　　（３）
　　ｈｙｒｅｆ－Δｈｙｍａｘ≦ｈｙ≦ｈｙｒｅｆ＋Δｈｙｍａｘ　　　（４）
　　ｈｚｒｅｆ－Δｈｚｍａｘ≦ｈｚ≦ｈｚｒｅｆ＋Δｈｚｍａｘ　　　（５）
　例えば、ゼロモーメンタム衛星の場合は、ｈｘｒｅｆ、ｈｙｒｅｆ、および、ｈｚｒｅｆはいずれも０である。バイアスモーメンタム衛星の場合は、ｈｘｒｅｆ、ｈｙｒｅｆ、ｈｚｒｅｆはいずれも０でない。Δｈｘｍａｘ、Δｈｙｍａｘ、および、Δｈｚｍａｘは、リアクションホイール、ＣＭＧ（Control Moment Gyro）等の静止衛星１０に搭載された姿勢制御アクチュエータで補償可能な角運動量の大きさに応じて定められる。Δｈｘｍａｘ、Δｈｙｍａｘ、および、Δｈｚｍａｘを、例えば３０Ｎｍｓとする。
　さらに、指令値算出部１３は、ジンバル機構の駆動角度の許容範囲を予め保持している。詳細には、指令値算出部１３は、ジンバル機構３１，３２，３３，３４の駆動角度の上限αｍａｘ，βｍａｘおよび下限－αｍａｘ，－βｍａｘを予め保持している。
　ジンバル機構３１，３２，３３，３４の駆動角度の上限αｍａｘ、βｍａｘおよび下限－αｍａｘ，－βｍａｘは、静止衛星１０の姿勢情報を取得するための上述のセンサ、姿勢制御アクチュエータ、姿勢制御アクチュエータを制御する制御装置等を含む姿勢制御系が補償可能な外乱トルクの大きさに応じて定められる。ジンバル機構３１，３２，３３，３４を回転させて、スラスタ２１，２２，２３，２４を噴射することで、静止衛星１０にトルクが働く。トルクの大きさＴは、下記（６）式で表されるように、ジンバル角の二乗ノルム、スラスタ２１，２２，２３，２４の推力の大きさＦ、および取付距離ｒに比例する。
Figure JPOXMLDOC01-appb-M000001
　静止衛星１０の打ち上げ後は、取付距離ｒを変更することはできない。スラスタ２１，２２，２３，２４の推力の大きさＦは、軌道上で一定とみなされる。そのため、トルクの大きさＴの調節は、ジンバル角α，βによって行われる。ジンバル角α，βを大きくすることで、より大きなトルクを静止衛星１０に与えることが可能である。ジンバル角α，βを大きくすることで、静止衛星１０に蓄積された角運動量をアンローディングする能力は向上する。しかしながら、スラスタ２１，２２，２３，２４の噴射によって静止衛星１０に作用するトルクは、静止衛星１０の姿勢制御系にとっては、外乱トルクである。そのため、トルクの大きさＴが大きすぎると、姿勢制御系が不安定になり、静止衛星１０の姿勢を保持できなくなる。そのため、ジンバル角α，βは、姿勢制御系で補償可能な外乱トルクの上限値Ｔｍａｘに応じて定められる。ジンバル角αの最大値αｍａｘは、下記（７）式で表される。ジンバル角βの最大値βｍａｘは、下記（８）式で表される。下記（７），（８）式におけるκは、１未満の正数である。
Figure JPOXMLDOC01-appb-M000002
Figure JPOXMLDOC01-appb-M000003
　ジンバル角α，βの範囲は、下記（９），（１０）式で表される。
　　－αｍａｘ≦α≦αｍａｘ　　　（９）
　　－βｍａｘ≦β≦βｍａｘ　　　（１０）
　指令値算出部１３は、上述の静止衛星１０の緯度および経度の範囲、角運動量の範囲、ならびにジンバル角α，βの範囲を予め保持している。指令値算出部１３は、これらの許容範囲に基づいて、非線形計画問題の解を求めることで、スラスタ指令値および角指令値を算出する。詳細には、指令値算出部１３は、後述する第１の制約条件の下で、目的関数を最小化する解を求める。図６を用いて、指令値算出部１３の全体の流れを説明した後、個々のステップについて順に詳細に説明する。第１関数算出部１３１は、複数の噴射区間でのスラスタ指令値および角指令値を引数として、第１の関数を求める（ステップＳ１１）。第２関数算出部１３２は、軌道決定部１１から平均軌道要素を取得し、角運動量算出部１２から角運動量を取得する（ステップＳ１２）。第２関数算出部１３２は、平均軌道要素、平均軌道要素の時間変化率、および静止衛星１０が有する角運動量を引数として、第２の関数を求める（ステップＳ１３）。第１制約条件設定部１３３は、軌道要素の変化量と平均軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することを示す条件式である第１の制約条件を求める（ステップＳ１４）。第２制約条件設定部１３４は、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にあることを示す条件式である第２の制約条件を求める（ステップＳ１５）。求解部１３５は、目的関数を求める（ステップＳ１６）。求解部１３５は、第１の制約条件、第２の制約条件、および複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計である目的関数に基づく非線形計画問題について、目的関数を最小化する解を求めることで、スラスタ指令値および角指令値を算出する（ステップＳ１７）。指令値算出部１３は、スラスタ指令値および角指令値をそれぞれ、スラスタ制御部１４およびジンバル制御部１５に出力する（ステップＳ１８）。
　図６の各ステップの詳細について順に説明する。上記ステップＳ１１の処理の詳細について説明する。第１関数算出部１３１は、複数の噴射区間でのスラスタ指令値および角指令値を引数として、第１の関数を求める。軌道の一周あたりＮ個の噴射区間が設けられる場合を例にして説明する。スラスタ指令値が、Ｎ個の噴射区間のそれぞれにおいて、スラスタ２１，２２，２３，２４の内、２本に対する噴射の指示である場合を例にして説明する。２本の組み合わせは、ＮＷスラスタ２１とＮＥスラスタ２２の組、ＳＷスラスタ２３とＳＥスラスタ２４の組、ＳＷスラスタ２３とＮＷスラスタ２１の組、および、ＮＥスラスタ２２とＳＥスラスタ２４の組のいずれかである。１つの噴射区間でスラスタ２１，２２，２３，２４の内、２本のみ噴射することで、推薬の使用量を少なくすることが可能である。
　ｊ（１≦ｊ≦Ｎ）番目の噴射区間において噴射するスラスタ２１，２２，２３，２４の内の２本を、ａｊ，ｂｊとする。ａｊとｂｊは、ＮＷスラスタ２１とＮＥスラスタ２２の組、ＳＷスラスタ２３とＳＥスラスタ２４の組、ＳＷスラスタ２３とＮＷスラスタ２１の組、および、ＮＥスラスタ２２とＳＥスラスタ２４の組のいずれかである。ｊ番目の噴射区間において噴射するスラスタａｊ，ｂｊの噴射量の合計をＦｊとすると、ａｊ，ｂｊのそれぞれの噴射量を決定する配分係数ξｊを用いて、スラスタａｊ，ｂｊのそれぞれの噴射量ｆａｊ，ｆｂｊは、下記（１１），（１２）式で表される。配分係数ξｊは、－１≦ξｊ≦１の条件を満たす係数である。ξｊ＝１の場合には、スラスタａｊのみが噴射し、ξｊ＝－１の場合は、スラスタｂｊのみが噴射し、ξｊ＝０の場合は、スラスタａｊ，ｂｊが等しく噴射する。
Figure JPOXMLDOC01-appb-M000004
Figure JPOXMLDOC01-appb-M000005
　第１の関数のパラメータは、ｊ番目の噴射区間における噴射量の合計Ｆｊ，ｊ番目の噴射区間における配分係数ξｊ，ｊ番目の噴射区間における平均経度λａｖｇｊ、およびｊ番目の噴射区間における角度αａｊ，βａｊ，αｂｊ，βｂｊである。角度αａｊ，βａｊは、スラスタａｊを構体５に取り付けるジンバル機構３１，３２，３３，３４のｊ番目の噴射区間における角度である。角度αｂｊ，βｂｊは、スラスタｂｊを構体５に取り付けるジンバル機構３１，３２，３３，３４のｊ番目の噴射区間における角度である。
　第１の関数の出力である軌道要素の変化量は、線形化されたガウスの惑星方程式に基づいて、噴射量の合計Ｆｊ、配分係数ξｊ、平均経度λａｖｇｊ、および角度αａｊ，βａｊ，αｂｊ，βｂｊで表される。軌道要素の変化量は、平均離心率ベクトル、平均傾斜角ベクトル、および平均直下点経度のそれぞれの変化量を含む。さらに第１の関数の出力である角運動量の変化量は、線形化されたオイラーの方程式に基づいて、噴射量の合計Ｆｊ、配分係数ξｊ、平均経度λａｖｇｊ、および角度αａｊ，βａｊ，αｂｊ，βｂｊで表される。
　Ｎ個の噴射区間におけるスラスタ２１，２２，２３，２４の噴射によって生じる軌道要素の変化量の合計および角運動量の変化量の合計は、下記（１３）－（１６）式で表される。下記（１３）式において、Ｇｅ（Ｆｊ，ξｊ，λａｖｇｊ，αａｊ，βａｊ，αｂｊ，βｂｊ）は、離心率ベクトルについてのガウスの惑星方程式である。下記（１４）式において、Ｇｉ（Ｆｊ，ξｊ，λａｖｇｊ，αａｊ，βａｊ，αｂｊ，βｂｊ）は、傾斜角ベクトルについてのガウスの惑星方程式である。下記（１５）式において、Ｇλ（Ｆｊ，ξｊ，λａｖｇｊ，αａｊ，βａｊ，αｂｊ，βｂｊ）は、平均直下点経度についてのガウスの惑星方程式である。下記（１６）式において、Ｅｈ（Ｆｊ，ξｊ，λａｖｇｊ，αａｊ，βａｊ，αｂｊ，βｂｊ）は、静止衛星１０が有する角運動量に関するオイラーの方程式である。第１関数算出部１３１は、上記パラメータを確定して、下記（１３）－（１６）式を含む第１の関数を求め、第１制約条件設定部１３３に送る。
Figure JPOXMLDOC01-appb-M000006
Figure JPOXMLDOC01-appb-M000007
Figure JPOXMLDOC01-appb-M000008
Figure JPOXMLDOC01-appb-M000009
　図６のステップＳ１１の処理の説明は以上であるから、図６のステップＳ１１と並行して行われるステップＳ１２，Ｓ１３の処理について以下に説明する。ステップＳ１２において、第２関数算出部１３２は、軌道決定部１１から平均軌道要素および平均軌道要素の時間変化率を取得し、角運動量算出部１２から角運動量を取得する。そして、ステップＳ１３において、第２関数算出部１３２は、平均軌道要素、平均軌道要素の時間変化率、および静止衛星１０が有する角運動量を引数として、パラメータを確定し、静止衛星１０の平均軌道要素の制御量および静止衛星１０が有する角運動量の制御量を出力する第２の関数を求める。平均軌道要素の制御量とは、スラスタ２１，２２，２３，２４を噴射することで達成すべき平均軌道要素の制御量である。制御対象となる平均軌道要素は、平均離心率ベクトル、平均傾斜角ベクトル、および平均直下点経度を含む。
　第２関数算出部１３２が第２の関数を算出するステップＳ１３の処理の概要について図７を用いて説明した後に、各ステップの詳細について説明する。第２関数算出部１３２は、平均離心率ベクトルの目標値ｅｔａｒｇｅｔを定める軌道座標系の原点を中心とする円の半径ｅＲを決定する（ステップＳ２１）。第２関数算出部１３２は、軌道座標系の原点を中心とする半径ｅＲの円と、地球４と太陽とを結ぶ線分の交点に設定される平均離心率ベクトルの目標値ｅｔａｒｇｅｔに平均離心率ベクトルを近づけるための制御量Δｅｃを算出する（ステップＳ２２）。第２関数算出部１３２は、緯度の保持精度Δφｍａｘの大きさに応じて、年平均傾斜角ベクトル、月平均傾斜角ベクトル、および日平均傾斜角ベクトルのいずれかを、制御対象である平均傾斜角ベクトルとして決定する（ステップＳ２３）。第２関数算出部１３２は、ステップＳ２３で決定した平均傾斜角ベクトルに対する制御量Δｉｃを算出する（ステップＳ２４）。第２関数算出部１３２は、平均直下点経度の制御量Δλｃを算出する（ステップＳ２５）。第２関数算出部１３２は、角運動量の制御量Δｈｃを算出する（ステップＳ２６）。図７に示すように、第２関数算出部１３２は、ステップＳ２１，Ｓ２３，Ｓ２５，Ｓ２６の処理を並行に行うことができる。
　上記ステップＳ２１，Ｓ２２の処理について説明する。第２関数算出部１３２は、軌道座標系の原点を中心とする半径ｅＲの円と、地球４と太陽とを結ぶ線分の交点に平均離心率ベクトルの目標値ｅｔａｒｇｅｔを設定し、下記（１７）式で表される平均離心率ベクトルの制御量Δｅｃを算出する。下記（１７）式において、ｅＯＤは、軌道決定部１１が算出した平均離心率ベクトルである。Ｋｅは正の係数である。下記（１７）式の第２項－Ｋｅ（ｅＯＤ－ｅｔａｒｇｅｔ）は、フィードバック制御量に相当する。下記（１７）式における第１項－δｅＦＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、静止衛星１０に働く摂動力、例えば、太陽輻射圧、月の潮汐力、太陽の潮汐力等による、平均離心率ベクトルの変動を補償する。下記（１７）式におえるδｅＦＦは、下記（１８）式で表される。下記（１８）式におけるｅ’ＯＤは、ｅＯＤの時間変化率であり、Ｔｃは、静止衛星１０が軌道を周回する周期である。
　　Δｅｃ＝－δｅＦＦ－Ｋｅ（ｅＯＤ－ｅｔａｒｇｅｔ）　　　（１７）
　　δｅＦＦ＝ｅ’ＯＤＴｃ　　　（１８）
　上記（１７）式における平均離心率ベクトルの目標値ｅｔａｒｇｅｔは、軌道座標系の原点を中心とする半径ｅＲの円と、地球４と太陽とを結ぶ線分の交点に設定される。ｅＲは、下記（１９）式で表される。ΔλＭは、軌道決定部１１が算出する平均軌道要素の誤差、スラスタ２１，２２，２３，２４の噴射誤差等によって生じる経度の誤差である。ΔλＭは、０＜ΔλＭ＜Δλｍａｘを満たす値であり、例えば、０．０２°である。
Figure JPOXMLDOC01-appb-M000010
　図７のステップＳ２１，Ｓ２２の処理については以上であるから、図７においてステップＳ２１，Ｓ２２と並行して行われるステップＳ２３，Ｓ２４の処理について説明する。第２関数算出部１３２は、平均傾斜角ベクトル、例えば、年平均傾斜角ベクトル、月平均傾斜角ベクトル、および日平均傾斜角ベクトル等を制御する。第２関数算出部１３２は、緯度の保持精度Δφｍａｘの大きさに適した平均傾斜角ベクトルを制御する。例えば、０．０４０°×（１＋κφ）≦Δφｍａｘの場合、第２関数算出部１３２は、年平均傾斜角ベクトルを制御対象として決定する。０．００４０°×（１＋κφ）＜Δφｍａｘ＜０．０４０°×（１＋κφ）の場合、第２関数算出部１３２は、月平均傾斜角ベクトルを制御対象として決定する。Δφｍａｘ＜０．００４°×（１＋κφ）の場合、第２関数算出部１３２は、日平均傾斜角ベクトルを制御対象として決定する。κφは、正の係数であり、例えば、１．１に設定される。
　平均傾斜角ベクトルの制御量Δｉｃは、下記（２０）式で表される。下記（２０）式において、ｉＯＤは、軌道決定部１１が算出した平均傾斜角ベクトルである。ｉｔａｒｇｅｔは、定められた平均傾斜角ベクトルの目標値である。Ｋｉは、正の定数である。下記（２０）式における第２項－Ｋｉ（ｉＯＤ－ｉｔａｒｇｅｔ）は、フィードバック制御量に相当する。下記（２０）式における第１項－δｉＦＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、静止衛星１０に働く摂動力による、平均傾斜角ベクトルの変動を補償する。δｉＦＦは、下記（２１）式で表される。下記（２１）式において、ｄｉＦＦ／ｄｔは、平均傾斜角ベクトルの時間変化率である。
　　Δｉｃ＝－δｉＦＦ－Ｋｉ（ｉＯＤ－ｉｔａｒｇｅｔ）　　　（２０）
Figure JPOXMLDOC01-appb-M000011
　図７のステップＳ２３，Ｓ２４の説明は以上であるから、図７においてステップＳ２３，Ｓ２４と並行して行われるステップＳ２５の処理について説明する。第２関数算出部１３２は、下記（２２）式で表される平均直下点経度の制御量Δλｃを算出する。下記（２２）式において、λＯＤは、軌道決定部１１で算出された平均直下点経度である。下記（２２）式において、λｔａｒｇｅｔは、平均直下点経度の目標値である。Ｋλは正の係数である。下記（２２）式における第２項－Ｋλ（λＯＤ－λｔａｒｇｅｔ）は、フィードバック制御量に相当する。下記（２２）式における第１項－δλＥＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、静止衛星１０に働く摂動力による、平均直下点経度の変動を補償する。δλＥＦは、下記（２３）式で表される。下記（２３）式において、λ’ＯＤは、平均直下点経度の角速度である。下記（２３）式においてλ”ＯＤは、平均直下点経度の角加速度である。
　　Δλｃ＝－δλＦＦ－Ｋλ（λＯＤ－λｔａｒｇｅｔ）　　　（２２）
Figure JPOXMLDOC01-appb-M000012
　図７のステップＳ２５の説明は以上であるから、図７においてステップＳ２５と並行して行われるステップＳ２６の処理について説明する。第２関数算出部１３２は、下記（２４）式で表される角運動量の制御量Δｈｃを算出する。下記（２４）式において、ｈＡＤは、角運動量算出部１２で算出された慣性系における静止衛星１０の角運動量である。下記（２４）式の第２項－Ｋｈ（ｈＡＤ－ｈｒｅｆ）は、フィードバック制御量に相当する。Ｋｈは、正の係数である。下記（２４）式の第１項－δｈＦＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、摂動トルクによって生じる慣性系での角運動量の変動量を補償する。摂動トルクの内、最も支配的なトルクは、太陽輻射圧によるトルクである。静止衛星１０に働く太陽輻射圧トルクをτＳＲＰとし、δｈＦＦを下記（２５）式で近似する。ｈｒｅｆは、慣性系における角運動量の目標値であり、下記（２６）式で表される。角運動量の目標値ｈｒｅｆの各成分は、予め指令値算出部１３が保持している。
　　Δｈｃ＝－δｈＦＦ－Ｋｈ（ｈＡＤ－ｈｒｅｆ）　　　（２４）
　　δｈＦＦ＝ＴｃτＳＲＰ　　　（２５）
　　ｈｒｅｆ＝［ｈｘｒｅｆ，ｈｙｒｅｆ，ｈｚｒｅｆ］Ｔ　　　（２６）
　第２関数算出部１３２は、上記（１７），（２０），（２２），（２４）式を含む第２の関数を第１制約条件設定部１３３に送る。
　図６のステップＳ１２，１３の処理の説明は以上であるから、後続のステップＳ１４の処理について説明する。第１制約条件設定部１３３は、軌道要素の変化量と平均軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することを示す条件式である第１の制約条件を求める。第１制約条件設定部１３３は、下記（２７）－（３０）式を含む第１の制約条件を求解部１３５に送る。
　　Δｅｃ＝δｅ　　　（２７）
　　Δｉｃ＝δｉ　　　（２８）
　　Δλｃ＝δλ　　　（２９）
　　Δｈｃ＝δｈ　　　（３０）
　図６のステップＳ１４の処理の説明は以上であるから、後続のステップＳ１５の処理について説明する。第２制約条件設定部１３４は、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にあることを示す条件式である第２の制約条件を求める。第２制約条件設定部１３４は、下記（３１）－（３４）式を含む第２の制約条件を求解部１３５に送る。
　　－αｍａｘ≦αａｊ≦αｍａｘ　　　（３１）
　　－βｍａｘ≦βａｊ≦βｍａｘ　　　（３２）
　　－αｍａｘ≦αｂｊ≦αｍａｘ　　　（３３）
　　－βｍａｘ≦βｂｊ≦βｍａｘ　　　（３４）
　スラスタ２１，２２，２３，２４の噴射量は０より大きいことは明らかであるから、下記（３５）式が成り立つ。また上述のように、噴射区間ｊにおける配分係数については、下記（３６）式が成り立つ。
　　０≦Ｆｊ　　　（３５）
　　－１≦ξｊ≦１　　　（３６）
　図６のステップＳ１５の説明は以上であるから、後続のステップＳ１６，Ｓ１７の処理について説明する。求解部１３５は、第１の制約条件、第２の制約条件、および複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計である目的関数に基づく非線形計画問題について、目的関数を最小化する解を求めることで、複数の噴射区間のそれぞれのスラスタ指令値および角指令値を算出する。目的関数は、下記（３７）式で表される。
Figure JPOXMLDOC01-appb-M000013
　求解部１３５は、噴射区間のそれぞれにおける、噴射量の合計Ｆｊ、配分係数ξｊ、平均経度λａｖｇｊ、および角度αａｊ，βａｊ，αｂｊ，βｂｊの解を求める。非線形計画問題を解くにあたっては、逐次二次計画法、内点法等の手法が用いられる。求解部１３５は、非線形計画問題を解くことで得られるλａｖｇｊの解から、スラスタ２１，２２，２３，２４の噴射区間の開始位置および終了位置を含む噴射タイミングを得る。また求解部１３５は、非線形計画問題を解くことで得られるＦｊ，ξｊの解から、噴射区間におけるスラスタ２１，２２，２３，２４の噴射量を得る。さらに求解部１３５は、非線形計画問題を解くことで得られるαａｊ，βａｊ，αｂｊ，βｂｊの解から、スラスタ２１，２２，２３，２４を噴射する際のジンバル角を得る。
　図６のステップＳ１６，Ｓ１７の説明は以上であるから、後続のステップＳ１８の処理について説明する。求解部１３５は、スラスタ２１，２２，２３，２４の噴射タイミングおよび噴射量を示すスラスタ指令値を、スラスタ制御部１４に送り、ジンバル角を示す角指令値を、ジンバル制御部１５に送る。
　スラスタ制御部１４は、スラスタ指令値に基づいてスラスタ２１，２２，２３，２４を制御して噴射させる。ジンバル制御部１５は、角指令値に基づいてジンバル機構３１，３２，３３，３４のジンバル角を制御する。上述のように非線形計画問題を解いて、スラスタ指令値および角指令値を算出することで、推薬の使用量を少なくしながら、軌道保持制御および角運動量のアンローディングを行うための処理を簡易化することが可能である。
　以上説明したとおり、本実施の形態１に係る宇宙機制御装置１によれば、非線形計画問題について、スラスタの噴射量の合計である目的関数を最小化する解を求め、スラスタ指令値および角指令値を算出することで、推薬の使用量を抑制しながら、軌道保持制御および角運動量のアンローディングを行うための処理を簡易化することが可能である。
　実施の形態１に係る宇宙機制御装置１によれば、昇交点および降交点から離隔した位置でスラスタ２１，２２，２３，２４を噴射する必要がある場合、また軌道を一周する間にスラスタ２１，２２，２３，２４を複数回噴射する必要がある場合でも、推薬の使用量を抑制しながら、軌道保持制御および角運動量のアンローディングを行うための処理を行うことができる。
　（実施の形態２）
　実施の形態１における非線形計画問題の制約条件は、第１の制約条件と第２の制約条件であったが、非線形計画問題は、さらに他の制約条件を含んでもよい。実施の形態２に係る宇宙機制御装置１の構成は、図４に示す実施の形態１に係る宇宙機制御装置１の構成と同様である。図８は、本発明の実施の形態２に係る指令値算出部の構成例を示す図である。実施の形態２に係る指令値算出部１６は、図５に示す実施の形態１に係る指令値算出部１３の構成に加えて、第３制約条件設定部１３６をさらに備える。第３制約条件設定部１３６は、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを示す条件式である第３の制約条件を求め、求解部１３５に送る。
　指令値算出部１６がスラスタ指令値および角指令値を算出する処理の概略について図９を用いて説明する。ステップＳ１１－Ｓ１６，Ｓ１８の処理は、図６に示す処理と同様である。実施の形態２に係る指令値算出部１６においては、第３制約条件設定部１３６が、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを示す条件式である第３の制約条件を求める（ステップＳ１９）。求解部１３５は、目的関数を求める（ステップＳ１６）。求解部１３５は、第１の制約条件、第２の制約条件、第３の制約条件、および複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計である目的関数に基づく非線形計画問題について、目的関数を最小化する解を求めることで、スラスタ指令値および角指令値を算出する（ステップＳ２０）。後続の処理は、図６に示す処理と同様である。
　上記（３１）－（３４）式における、αｍａｘ，βｍａｘを小さくすると、スラスタ２１，２２，２３，２４の噴射中に静止衛星１０に作用するトルクは小さくなる。すなわち、スラスタ２１，２２，２３，２４の噴射によって生じる角運動量の変動が抑制される。αｍａｘ，βｍａｘを小さくすることで、角運動量の変動が抑制されても、噴射区間の噴射直前または噴射直後において、角運動量が、上記（３）－（５）式で規定される範囲内にないことがある。そこで、第３制約条件設定部１３６は、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを示す条件式である第３の制約条件を求める。すなわち、第３制約条件設定部１３６は、下記（３８），（３９）式の少なくともいずれかが成立することを示す条件式である第３の制約条件を求める。下記（３８）式におけるｈｊ （－）は、ｊ番目の噴射区間の噴射直前における角運動量である。ｊ番目の噴射区間の噴射直前における角運動量とは、ｊ番目の噴射区間におけるスラスタ２１，２２，２３，２４の噴射開始直前の静止衛星１０の角運動量である。下記（３９）式におけるｈｊ （＋）は、ｊ番目の噴射区間の噴射直後における角運動量である。ｊ番目の噴射区間の噴射直後における角運動量とは、ｊ番目の噴射区間におけるスラスタ２１，２２，２３，２４の噴射終了直後の静止衛星１０の角運動量である。ｈｊ （－）およびｈｊ （＋）は、変数Ｆｊ，ξｊ，λａｖｇｊ，αａｊ，βａｊ，αｂｊ，βｂｊの関数である。
　　ｈｒｅｆ－Δｈｍａｘ≦ｈｊ （－）≦ｈｒｅｆ＋ｈｍａｘ　　　（３８）
　　ｈｒｅｆ－Δｈｍａｘ≦ｈｊ （＋）≦ｈｒｅｆ＋ｈｍａｘ　　　（３９）
　ｈｊ （－）の各成分について、上記（３）－（５）式と同様に、下記（４０）－（４２）が成り立つ。またｈｊ （＋）の各成分について、上記（３）－（５）式と同様に、下記（４３）－（４５）が成り立つ。
　　ｈｘｒｅｆ－Δｈｘｍａｘ≦ｈｘｊ （－）≦ｈｘｒｅｆ＋Δｈｘｍａｘ　　　（４０）
　　ｈｙｒｅｆ－Δｈｙｍａｘ≦ｈｙｊ （－）≦ｈｙｒｅｆ＋Δｈｙｍａｘ　　　（４１）
　　ｈｚｒｅｆ－Δｈｚｍａｘ≦ｈｚｊ （－）≦ｈｚｒｅｆ＋Δｈｚｍａｘ　　　（４２）
　　ｈｘｒｅｆ－Δｈｘｍａｘ≦ｈｘｊ （＋）≦ｈｘｒｅｆ＋Δｈｘｍａｘ　　　（４３）
　　ｈｙｒｅｆ－Δｈｙｍａｘ≦ｈｙｊ （＋）≦ｈｙｒｅｆ＋Δｈｙｍａｘ　　　（４４）
　　ｈｚｒｅｆ－Δｈｚｍａｘ≦ｈｚｊ （＋）≦ｈｚｒｅｆ＋Δｈｚｍａｘ　　　（４５）
　以上説明したとおり、本実施の形態２に係る宇宙機制御装置１によれば、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを、非線形計画問題の第３の制約条件とすることで、静止衛星１０の角運動量を保持範囲内に維持することができる。そのため、姿勢制御アクチュエータの角運動量を飽和させることなく、軌道保持制御および角運動量のアンローディングを行うことが可能である。
　次に、本実施の形態のハードウェア構成について説明する。図１０は、ハードウェア構成の例を示す図である。図１０は、宇宙機制御装置１が処理回路１０１によって実現される例を示している。図１０のように、静止衛星１０は、衛星の位置、速度、姿勢角、角速度を観測するセンサ１１０、宇宙機制御装置１、スラスタ２１，２２，２３，２４、およびジンバル機構３１，３２，３３，３４を有する。
　処理回路１０１が、専用のハードウェアである場合、処理回路１０１は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Application Specific Integrated Circuit）、ＦＰＧＡ（Field Programmable Gate Array）、またはこれらを組み合わせたものが該当する。軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の各部の機能それぞれを処理回路１０１で実現してもよいし、各部の機能をまとめて処理回路１０１で実現してもよい。
　図１１は、本実施の形態のハードウェア構成の別の例を示す図である。図１１は、宇宙機制御装置１がプロセッサ（ＣＰＵ：Central Processing Unit）１０２とメモリ１０３によって実現される例を示している。図１１のプロセッサ１０２とメモリ１０３は、図１０の処理回路１０１に相当する。
　宇宙機制御装置１がプロセッサ１０２とメモリ１０３によって実現される場合、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の機能は、ソフトウェア、ファームウェア、またはソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェアやファームウェアはプログラムとして記述され、メモリ１０３に格納される。プロセッサ１０２は、メモリ１０３に記憶されたプログラムを読み出して実行することにより、各部の機能を実現する。すなわち、メモリ１０３には、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の処理を実行するためのプログラムが格納される。メモリ１０３は、例えば、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read-Only Memory）、フラッシュメモリー、ＥＰＲＯＭ（Erasable Programmable Read Only Memory）、ＥＥＰＲＯＭ（Electrically Erasable and Programmable Read-Only Memory）等の、不揮発性または揮発性の半導体メモリ、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤ（Digital Versatile Disc）等を含む。
　メモリ１０３は、上述のプログラムの他、軌道決定部１１で決定される軌道決定値、角運動量算出部１２で設定される衛星角運動量値、指令値算出部１３で算出されるスラスタ指令値およびジンバル角指令値を記憶する。上述のメモリ１０３に記憶された情報は適宜読みだされて、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の処理に用いられ、修正される。
　なお、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の各機能について、一部を専用のハードウェアで実現し、一部をソフトウェアまたはファームウェアで実現するようにしてもよい。例えば、軌道決定部１１については専用のハードウェアとしての処理回路１０１でその機能を実現し、指令値算出部１３については、プロセッサ１０２がメモリ１０３に格納されたプログラムを読み出して実行することによってその機能を実現することが可能である。
　本発明の実施の形態は上述の実施の形態に限られない。宇宙機は、静止衛星１０に限られず、天体の周りの軌道を周回する任意の飛行体である。飛行体の一例は、地球観測衛星、宇宙船等である。スラスタ２１，２２，２３，２４の本数は任意であり、それぞれの噴射区間において噴射するスラスタ２１，２２，２３，２４の本数および組み合わせも任意である。またジンバル機構３１，３２，３３，３４の自由度は１でもよい。ジンバル機構３１，３２，３３，３４の自由度が１である場合、上記数式のβ＝０となるため、第１の関数は、変数Ｆｊ，ξｊ，λａｖｇｊ，αａｊ，αｂｊの関数である。β＝０の場合、第１の関数と同様に、ｈｊ （－）およびｈｊ （＋）は、変数Ｆｊ，ξｊ，λａｖｇｊ，αａｊ，αｂｊの関数である。図１の例では、宇宙機の構体５に、スラスタ２１，２２，２３，２４が取り付けられているが、宇宙機の構体５に取り付けられるスラスタ２１，２２，２３，２４の数は任意である。スラスタ２１，２２，２３，２４のそれぞれの角度ψおよび角度γは互いに同じでもよいし、異なってもよい。
　緯度および経度の範囲、ならびに、緯度および経度の保持精度は、上述の例に限られず、任意に定めることができる。上述のように、軌道決定部１１が平均軌道要素の時間変化率を算出する代わりに、第２関数算出部１３２が平均軌道要素の時間変化率を算出してもよい。軌道決定部１１は、軌道要素を算出する代わりに、外部装置から起動要素を取得してもよい。軌道決定部１１は、軌道要素の瞬時値を指令値算出部１３に送ってもよい。第２関数算出部１３２は、図７に示す、平均離心率ベクトルの制御量Δｅｃの算出、平均傾斜角ベクトルΔｉｃの算出、平均直下点経度の制御量Δλｃの算出、および、角運動量の制御量Δｈｃの算出を順次行ってもよい。
　求解部１３５は、第１の制約条件の下で、非線形計画問題の解を求めることで、スラスタ指令値およびジンバル指令値を算出してもよい。あるいは、求解部１３５は、第１の制約条件と第３の制約条件の下で、非線形計画問題の解を求めてもよい。また、求解部１３５は、第１の制約条件の下で、目的関数を小さくする準最適解を求めることで、スラスト指令値およびジンバル指令値を算出してもよい。
　本発明は、本発明の広義の精神と範囲を逸脱することなく、様々な実施の形態及び変形が可能とされるものである。また、上述した実施の形態は、この発明を説明するためのものであり、本発明の範囲を限定するものではない。すなわち、本発明の範囲は、実施の形態ではなく、特許請求の範囲によって示される。そして、特許請求の範囲内及びそれと同等の発明の意義の範囲内で施される様々な変形が、この発明の範囲内とみなされる。
　本出願は、２０１７年９月１日に出願された、日本国特許出願特願２０１７－１６８６４９号に基づく。本明細書中に日本国特許出願特願２０１７－１６８６４９号の明細書、特許請求の範囲、図面全体を参照として取り込むものとする。
　１　宇宙機制御装置、４　地球、５　構体、１０　静止衛星、１１　軌道決定部、１２　角運動量算出部、１３，１６　指令値算出部、１４　スラスタ制御部、１５　ジンバル制御部、２１，２２，２３，２４　スラスタ、３１，３２，３３，３４　ジンバル機構、５０，５１　面、１０１　処理回路、１０２　プロセッサ、１０３　メモリ、１１０　センサ、１３１　第１関数算出部、１３２　第２関数算出部、１３３　第１制約条件設定部、１３４　第２制約条件設定部、１３５　求解部、１３６　第３制約条件設定部。

Claims (14)

  1. 　天体の周りの軌道を周回する宇宙機が有する角運動量を算出する角運動量算出部と、
    　それぞれが前記宇宙機の構体に、ジンバル機構を介して取り付けられた複数のスラスタの内、少なくともいずれかのスラスタの噴射を指示するスラスタ指令値、および前記スラスタの噴射時における前記ジンバル機構の角度を指示する角指令値を算出する指令値算出部と、
    　前記スラスタ指令値に基づき、前記スラスタを制御するスラスタ制御部と、
    　前記角指令値に基づき、前記ジンバル機構を制御するジンバル制御部と、
    　を備え、
    　前記指令値算出部は、前記軌道における複数の噴射区間での前記スラスタ指令値および前記角指令値を引数として、該スラスタ指令値に基づいて前記スラスタが制御され、該角指令値に基づいて前記ジンバル機構が制御された場合に、前記宇宙機が前記軌道を一周する間の、前記宇宙機の軌道を表す軌道要素の変化量と前記角運動量の変化量とを出力する第１の関数、ならびに、前記軌道要素、前記軌道要素の時間変化率、および前記角運動量を引数として、前記宇宙機が前記軌道を一周する間の、前記軌道要素の制御量および前記角運動量の制御量を出力する第２の関数に基づき、前記軌道要素の変化量と前記軌道要素の制御量とが一致し、かつ、前記角運動量の変化量と前記角運動量の制御量とが一致する第１の制約条件の下で、前記複数の噴射区間における前記スラスタの噴射量の合計が目的関数である非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    　宇宙機制御装置。
  2. 　前記指令値算出部は、前記目的関数を最小化する前記非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    　請求項１に記載の宇宙機制御装置。
  3. 　前記指令値算出部は、前記ジンバル機構の角度が定められた範囲内にある第２の制約条件の下で、前記非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    　請求項１または２に記載の宇宙機制御装置。
  4. 　前記指令値算出部は、前記ジンバル機構の角度が定められた範囲内にある第２の制約条件と、前記複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、前記角運動量が定められた範囲内にある第３の制約条件との下で、前記非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    　請求項１または２に記載の宇宙機制御装置。
  5. 　前記第２の関数が算出する前記軌道要素の制御量は、前記宇宙機が前記軌道を一周する間に、前記宇宙機に作用する摂動力によって生じる前記軌道要素の変動を補正するフィードフォワード制御量、および、前記軌道要素と前記軌道要素の目標値とに基づくフィードバック制御量を含む、
    　請求項１から４のいずれか１項に記載の宇宙機制御装置。
  6. 　前記第２の関数が算出する前記角運動量の制御量は、前記宇宙機が前記軌道を一周する間に、前記宇宙機に作用する摂動トルクによって生じる前記角運動量の変動を補正するフィードフォワード制御量、および、前記角運動量と前記角運動量の目標値とに基づくフィードバック制御量を含む、
    　請求項１から５のいずれか１項に記載の宇宙機制御装置。
  7. 　前記軌道要素は、離心率ベクトル、傾斜角ベクトル、および直下点経度を含む、
    　請求項１から６のいずれか１項に記載の宇宙機制御装置。
  8. 　前記第２の関数は、前記宇宙機の軌道保持制御において求められる緯度の保持精度に応じた、年平均傾斜角ベクトル、月平均傾斜角ベクトル、および日平均傾斜角ベクトルのいずれかに対する制御量を含む前記軌道要素の制御量を出力する、
    　請求項７に記載の宇宙機制御装置。
  9. 　前記構体の外面の内、いずれかの面は、常に前記天体に向いており、
    　前記複数のスラスタは、噴射方向が互いに異なる向きで、前記構体の外面の内、常に前記天体に向いている面と反対側に位置する面に取り付けられる、
    　請求項１から８のいずれか１項に記載の宇宙機制御装置。
  10. 　前記ジンバル機構の自由度は２であり、
    　前記複数のスラスタの個数は４つであり、
    　前記スラスタは、前記宇宙機の質量中心から離れる方向に伸び、
    　前記噴射区間において４つの前記スラスタの内、２つの前記スラスタが噴射される、
    　請求項１から９のいずれか１項に記載の宇宙機制御装置。
  11. 　前記ジンバル機構の角度は、前記スラスタの噴射によって前記宇宙機に作用するトルクを、前記宇宙機が有する姿勢制御系で補償可能な外乱トルクの上限値以下にする範囲内の値である、
    　請求項１から１０のいずれか１項に記載の宇宙機制御装置。
  12. 　前記スラスタ指令値は、前記スラスタの噴射位相、前記複数のスラスタの噴射量の合計、前記合計において前記スラスタ指令値によって噴射が指示される前記スラスタの噴射量の割合である、
    　請求項１から１１のいずれか１項に記載の宇宙機制御装置。
  13. 　軌道における複数の噴射区間でのスラスタ指令値および角指令値を引数として、該スラスタ指令値に基づいてスラスタが制御され、該角指令値に基づいてジンバル機構が制御された場合に、軌道要素の変化量および角運動量の変化量を出力する第１の関数、ならびに、軌道要素、前記軌道要素の時間変化率、および前記角運動量を引数として、前記軌道要素の制御量および前記角運動量の制御量を出力する第２の関数に基づき、予め定めた制約条件を満たしつつ、前記複数の噴射区間における前記スラスタの噴射量の合計が目的関数である非線形計画問題について、前記目的関数の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    　宇宙機制御方法。
  14. 　コンピュータを、
    　天体の周りの軌道を周回する宇宙機が有する角運動量を算出する角運動量算出部、
    　それぞれが前記宇宙機の構体に、ジンバル機構を介して取り付けられた複数のスラスタの内、少なくともいずれかのスラスタの噴射を指示するスラスタ指令値、および前記スラスタの噴射時における前記ジンバル機構の角度を指示する角指令値を算出する指令値算出部、
    　前記スラスタ指令値に基づき、前記スラスタを制御するスラスタ制御部、および、
    　前記角指令値に基づき、前記ジンバル機構を制御するジンバル制御部、
    　として機能させるためのプログラムであって、
    　前記指令値算出部は、前記軌道における複数の噴射区間での前記スラスタ指令値および前記角指令値を引数として、該スラスタ指令値に基づいて前記スラスタが制御され、該角指令値に基づいて前記ジンバル機構が制御された場合に、前記宇宙機が前記軌道を一周する間の、前記宇宙機の軌道を表す軌道要素の変化量と前記角運動量の変化量とを出力する第１の関数、ならびに、前記軌道要素、前記軌道要素の時間変化率、および前記角運動量を引数として、前記宇宙機が前記軌道を一周する間の、前記軌道要素の制御量および前記角運動量の制御量を出力する第２の関数に基づき、前記軌道要素の変化量と前記軌道要素の制御量とが一致し、かつ、前記角運動量の変化量と前記角運動量の制御量とが一致する第１の制約条件の下で、前記複数の噴射区間における前記スラスタの噴射量の合計が目的関数である非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    　プログラム。
PCT/JP2018/031486 2017-09-01 2018-08-27 宇宙機制御装置、宇宙機制御方法、およびプログラム WO2019044735A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2019539474A JP6707206B2 (ja) 2017-09-01 2018-08-27 宇宙機制御装置、宇宙機制御方法、およびプログラム

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2017-168649 2017-09-01
JP2017168649 2017-09-01

Publications (1)

Publication Number Publication Date
WO2019044735A1 true WO2019044735A1 (ja) 2019-03-07

Family

ID=65527282

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2018/031486 WO2019044735A1 (ja) 2017-09-01 2018-08-27 宇宙機制御装置、宇宙機制御方法、およびプログラム

Country Status (2)

Country Link
JP (1) JP6707206B2 (ja)
WO (1) WO2019044735A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021011257A (ja) * 2019-07-03 2021-02-04 三菱電機株式会社 結合された天文座標系の非線形モデル予測制御

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH06144398A (ja) * 1992-11-04 1994-05-24 Mitsubishi Electric Corp 人工衛星の姿勢制御装置
JPH09277997A (ja) * 1995-12-22 1997-10-28 He Holdings Inc Dba Hughes Electron ジンバル制御されたスラスタを使用したモーメンタム制動
JP2006240375A (ja) * 2005-03-01 2006-09-14 Mitsubishi Electric Corp 人工衛星の姿勢制御装置
JP2012051387A (ja) * 2010-08-31 2012-03-15 Mitsubishi Electric Corp 宇宙機の姿勢制御装置
JP2016124538A (ja) * 2015-01-07 2016-07-11 三菱電機株式会社 宇宙船の動作を制御する方法および制御システム、並びに宇宙船
WO2016111317A1 (ja) * 2015-01-09 2016-07-14 三菱電機株式会社 軌道制御装置および衛星
US9522746B1 (en) * 2015-08-27 2016-12-20 The Boeing Company Attitude slew methodology for space vehicles using gimbaled low-thrust propulsion subsystem

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH06144398A (ja) * 1992-11-04 1994-05-24 Mitsubishi Electric Corp 人工衛星の姿勢制御装置
JPH09277997A (ja) * 1995-12-22 1997-10-28 He Holdings Inc Dba Hughes Electron ジンバル制御されたスラスタを使用したモーメンタム制動
JP2006240375A (ja) * 2005-03-01 2006-09-14 Mitsubishi Electric Corp 人工衛星の姿勢制御装置
JP2012051387A (ja) * 2010-08-31 2012-03-15 Mitsubishi Electric Corp 宇宙機の姿勢制御装置
JP2016124538A (ja) * 2015-01-07 2016-07-11 三菱電機株式会社 宇宙船の動作を制御する方法および制御システム、並びに宇宙船
WO2016111317A1 (ja) * 2015-01-09 2016-07-14 三菱電機株式会社 軌道制御装置および衛星
US9522746B1 (en) * 2015-08-27 2016-12-20 The Boeing Company Attitude slew methodology for space vehicles using gimbaled low-thrust propulsion subsystem

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021011257A (ja) * 2019-07-03 2021-02-04 三菱電機株式会社 結合された天文座標系の非線形モデル予測制御
JP7275078B2 (ja) 2019-07-03 2023-05-17 三菱電機株式会社 結合された天文座標系の非線形モデル予測制御

Also Published As

Publication number Publication date
JP6707206B2 (ja) 2020-06-10
JPWO2019044735A1 (ja) 2020-04-23

Similar Documents

Publication Publication Date Title
US8113468B2 (en) Precision attitude control system for gimbaled thruster
EP0769736B1 (en) Method for inclined orbit attitude control for momentum bias spacecraft
JP6656380B2 (ja) 宇宙機の動作を制御する方法及びシステム、宇宙機
US5984236A (en) Momentum unloading using gimbaled thrusters
US5098041A (en) Attitude control system for momentum-biased spacecraft
Bedrossian et al. Zero-propellant maneuver guidance
US6237876B1 (en) Methods for using satellite state vector prediction to provide three-axis satellite attitude control
US5528502A (en) Satellite orbit maintenance system
JPH03189297A (ja) 衛星のロール及びヨー姿勢の制御方法
US20220097872A1 (en) Attitude control system and method
CN107850900A (zh) 用于小型卫星的快速转动和安置系统
EP3854698B1 (en) Orientation control device, satellite, orientation control method, and program
Sugimura et al. Attitude determination and control system for nadir pointing using magnetorquer and magnetometer
US7464898B1 (en) Precision thrust/sun tracking attitude control system for gimbaled thruster
JP6707206B2 (ja) 宇宙機制御装置、宇宙機制御方法、およびプログラム
JP2021011257A (ja) 結合された天文座標系の非線形モデル予測制御
WO2016125145A1 (en) Method and system for station keeping of geo satellites
US5996942A (en) Autonomous solar torque management
Wong et al. An attitude control design for the Cassini spacecraft
Koh et al. In-orbit performance of attitude control system in DubaiSat-2
US20240182185A1 (en) Phased control of multiple spacecraft during a low-thrust orbit transfer maneuver
Ashok et al. Cmg configuration and steering approach for spacecraft rapid maneuvers
Ali et al. In-flight Correction of the Satellite Orientation Parameter during Target Mode.
Daniel et al. MESSENGER SPACECRAFT POINTING OPTIONS
Brás et al. Evaluating strategies for ground track acquisition and orbital phasing

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: 18850852

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2019539474

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 18850852

Country of ref document: EP

Kind code of ref document: A1