CN102033491B

CN102033491B - Method for controlling flexible satellite based on feature model

Info

Publication number: CN102033491B
Application number: CN2010102979615A
Authority: CN
Inventors: 孟斌; 吴宏鑫; 杨孟飞
Original assignee: Beijing Institute of Control Engineering
Current assignee: Beijing Institute of Control Engineering
Priority date: 2010-09-29
Filing date: 2010-09-29
Publication date: 2012-08-22
Anticipated expiration: 2030-09-29
Also published as: CN102033491A

Abstract

A characteristic model-based flexible satellite control method, according to the dynamic equation of the flexible satellite, determine its time scale, sampling time, and parameters M and m; determine the coefficient range of the characteristic model according to the obtained variables; use the gradient method Identify the parameters of the characteristic model; design the control law according to the coefficients of the characteristic model obtained from the identification, and feed back the control law to the dynamic equation of the flexible satellite to control the attitude angle of the flexible satellite. The present invention introduces the time scale and sampling period of the flexible satellite, describes the change rate of the flexible satellite, and solves the bottleneck problem of the flexible satellite characteristic modeling; the present invention provides the expression of the parameter range of the flexible satellite characteristic model , qualitatively studied the parameter properties of the characteristic model. From the given parameter range, it can be seen that the boundary of the characteristic model parameters is related to the sampling period, modeling error, system order, and system change rate. It is a flexible satellite based on the characteristic model. The adaptive control of the aircraft has laid a theoretical foundation; the invention is suitable for the characteristic modeling of the dynamics of the attitude of the aircraft, thus laying the foundation for the attitude control of the aircraft based on the characteristic model.

Description

A Characteristic Model-Based Control Method for Flexible Satellites

技术领域 technical field

本发明涉及一种卫星控制方法，特别是涉及一种基于特征模型的挠性卫星控制方法，属于卫星控制技术领域。The invention relates to a satellite control method, in particular to a characteristic model-based flexible satellite control method, which belongs to the technical field of satellite control.

背景技术 Background technique

基于特征模型的全系数自适应控制方法是吴宏鑫院士提出的，经过20多年的研究，在理论和应用上均取得了重要进展，形成了一套实用性很强的自适应控制理论和方法。该方法需要辨识参数少，可以保证闭环系统的暂态性能和稳态性能，具有强鲁棒性和自适应性。特别是该方法的理论思想和工程要点被创造性地应用于飞船返回再入控制，其开伞精度达到世界先进水平。The full-coefficient adaptive control method based on the characteristic model was proposed by Academician Wu Hongxin. After more than 20 years of research, important progress has been made in theory and application, and a set of highly practical adaptive control theory and method has been formed. This method requires fewer identification parameters, can guarantee the transient performance and steady-state performance of the closed-loop system, and has strong robustness and adaptability. In particular, the theoretical ideas and engineering points of this method have been creatively applied to the control of spacecraft return and reentry, and its parachute opening accuracy has reached the world's advanced level.

基于特征模型的全系数自适应控制方法的基本思想是，首先建立系统的特征模型及其参数范围，然后按照特征模型参数设计全系数自适应控制律。一般来说，系统的特征模型用系数有界的二阶时变差分方程描述。在基于特征模型的自适应控制设计中，首先要确定二阶时变差分方程的系数范围，然后在该范围内选取辨识初值，并把每一步辨识的结果投影到该范围内，因此系数范围的确定是基于特征模型控制方法的关键问题之一。线性定常系统特征模型参数范围确定问题已经解决，而对于非线性系统，该问题一直是制约全系数自适应控制设计的瓶颈。The basic idea of the full-coefficient adaptive control method based on the characteristic model is to first establish the characteristic model of the system and its parameter range, and then design the full-coefficient adaptive control law according to the characteristic model parameters. In general, the characteristic model of the system is described by second-order time-varying difference equations with bounded coefficients. In the adaptive control design based on the characteristic model, the coefficient range of the second-order time-varying difference equation must first be determined, and then the identification initial value is selected within this range, and the identification results of each step are projected into this range, so the coefficient range The determination of is one of the key issues in the characteristic model-based control method. The problem of determining the parameter range of the characteristic model of the linear steady system has been solved, but for the nonlinear system, this problem has always been the bottleneck restricting the design of the full coefficient adaptive control.

非线性黄金分割自适应控制、吴宏鑫、王颖、解永春著，宇航学报出版(2002，23(6)：1-8.2)和多变量线性时变系统的特征模型及自适应模糊控制方法，孙多青、吴宏鑫著，宇航学报出版(2005，26(6)：677-681.)中公开了几类特殊形式的非线性系统转化成二阶时变差分方程组的方法，不足之处是所给出的方法只适用于所考虑的特殊形式的非线性系统，并且未考虑建模误差。对于本发明所给出的任意建模误差下挠性卫星特征模型参数界的确定方法，没有公开的具有完整实用意义的方法。Nonlinear Golden Section Adaptive Control, Wu Hongxin, Wang Ying, Xie Yongchun, Acta Astronautics Publishing (2002, 23(6): 1-8.2) and the characteristic model and adaptive fuzzy control method of multivariable linear time-varying systems, Sun Duo Qing, Wu Hongxin, Acta Astronautics Publishing (2005, 26(6): 677-681.) discloses the method of transforming several types of nonlinear systems into second-order time-varying difference equations. The presented method is only applicable to the particular form of nonlinear system considered and does not take into account modeling errors. As for the method for determining the parameter boundary of the flexible satellite characteristic model under any modeling error given by the present invention, there is no disclosed method with complete practical significance.

发明内容 Contents of the invention

本发明的技术解决问题是：克服现有技术的不足，提供一种基于特征模型的挠性卫星控制方法。The technical problem of the present invention is: to overcome the deficiencies of the prior art, and provide a flexible satellite control method based on a feature model.

本发明的技术解决方案是：一种基于特征模型的挠性卫星控制方法，通过以下步骤实现：The technical solution of the present invention is: a kind of flexible satellite control method based on characteristic model, realizes by following steps:

第一步，利用公式(2)确定挠性卫星动力学方程的时间尺度p，In the first step, use formula (2) to determine the time scale p of the dynamic equation of the flexible satellite,

$p p = = min min {{\frac{11}{\sqrt{{M m}_{{\overset{\cdot \cdot}{f f}}_{11}}}},, \frac{11}{\sqrt{{M m}_{\overset{\cdot &Center Dot;}{f f}}}},, \frac{11}{\sqrt{{M m}_{{f f}_{11}}}},, \frac{11}{\sqrt{{M m}_{{f f}_{22}}}},, \frac{11}{\sqrt{{M m}_{u u}}}}} - - - - - - ((22))$

其中挠性卫星动力学方程为公式组(1)，Among them, the dynamic equation of the flexible satellite is the formula group (1),

${[\begin{matrix} \overset{\cdot &Center Dot;}{φ φ} & \overset{\cdot &Center Dot;}{θ θ} & \overset{\cdot &Center Dot;}{ψ ψ} \end{matrix}]}^{T T} = = C C (({x x}_{11})) {[\begin{matrix} {w w}_{x x} & {w w}_{y the y} & {w w}_{z z} \end{matrix}]}^{T T}$

$\{\begin{matrix} {I I}_{s the s} {\overset{\cdot &Center Dot;}{w w}}_{s the s} + + {\overset{~ ~}{w w}}_{s the s} {I I}_{s the s} {w w}_{s the s} + + {F f}_{sl sl} {\overset{\cdot &Center Dot; \cdot &Center Dot;}{η η}}_{l l} + + {F f}_{sr sr} {\overset{\cdot &Center Dot; \cdot &Center Dot;}{η η}}_{r r} = = {T T}_{s the s} \\ {\overset{\cdot &Center Dot; \cdot &Center Dot;}{η η}}_{l l} + + 22 {ξ ξ}_{l l} {w w}_{l l} {\overset{\cdot &Center Dot;}{η η}}_{l l} + + {w w}_{l l}^{22} {\overset{\cdot &Center Dot;}{η η}}_{l l} + + {F f}_{sl sl}^{T T} {\overset{\cdot \cdot}{w w}}_{s the s} = = 00 \\ {\overset{\cdot \cdot \cdot \cdot}{η η}}_{r r} + + 22 {ξ ξ}_{r r} {w w}_{r r} {\overset{\cdot &Center Dot;}{η η}}_{r r} + + {w w}_{r r}^{22} {\overset{\cdot \cdot}{η η}}_{r r} + + {F f}_{sr sr}^{T T} {\overset{\cdot &Center Dot;}{w w}}_{s the s} = = 00 \\ y the y = = [\begin{matrix} φ φ & θ θ & ψ ψ \end{matrix}] \end{matrix} - - - - - - ((11))$

$C C (({x x}_{11})) = = [\begin{matrix} \frac{cos cos θ θ}{cos cos ψ ψ} & 00 & \frac{sin sin θ θ}{cos cos ψ ψ} \\ tan the tan ψ ψ cos cos θ θ & 11 & tan the tan ψ ψ sin sin θ θ \\ - - sin sin θ θ & 00 & cos cos θ θ \end{matrix}],,$ ${\overset{~ ~}{w w}}_{s the s} = = [\begin{matrix} 00 & - - {w w}_{z z} & {w w}_{y the y} \\ {w w}_{z z} & 00 & - - {w w}_{x x} \\ - - {w w}_{y the y} & {w w}_{x x} & 00 \end{matrix}]$

φ，θ，ψ表示挠性卫星的俯仰、偏航和滚动姿态角，[w_x w_y w_z]^T表示卫星相对轨道坐标系的角速度在体坐标系中的坐标，w_s、分别表示挠性卫星中心体的角速度列阵和反对称阵，η_l、η_r分别为挠性卫星左、右太阳翼的模态坐标阵，ξ_l、ξ_r分别为挠性卫星左、右太阳翼的模态阻尼系数，F_sl、F_sr分别为挠性卫星左、右太阳翼与中心体的耦合系数，T_s表示作用在挠性卫星上的外力矩列阵，I_s表示挠性卫星惯量阵，x₁＝[φθψ]^T，y表示挠性卫星输出，w_l、w_r分别为挠性卫星左、右太阳翼的角速度，φ, θ, ψ represent the pitch, yaw and roll attitude angles of the flexible satellite, [w _x w _y w _z ] ^T represents the coordinates of the satellite’s angular velocity relative to the orbital coordinate system in the body coordinate system, w _s , represent the angular velocity array and antisymmetric array of the flexible satellite center body respectively, η _l and η _r are the modal coordinate arrays of the left and right solar wings of the flexible satellite respectively, and ξ _l and ξ _r are the left and right The modal damping coefficient of the solar wing, F _sl and F _sr are the coupling coefficients between the left and right solar wings of the flexible satellite and the central body, respectively, T _s represents the external force matrix acting on the flexible satellite, and I _s represents the flexible Satellite inertia array, x ₁ =[φθψ] ^T , y represents the output of the flexible satellite, w _l and w _r are the angular velocities of the left and right solar wings of the flexible satellite, respectively,

f₁＝C(x₁)w_s， $f_{2} = - I_{s}^{- 1} ({\tilde{w}}_{s} I_{s} w_{s} + F_{sl} {\overset{\cdot \cdot}{η}}_{l} + F_{sr} {\overset{\cdot \cdot}{η}}_{r}),$ $g = I_{s}^{- 1},$ u＝T_s，f＝f₂+gu，f ₁ =C(x ₁ )w _s , $f_{2} = - I_{the s}^{- 1} ({\tilde{w}}_{the s} I_{the s} w_{the s} + f_{sl} {\overset{&Center Dot; &Center Dot;}{η}}_{l} + f_{sr} {\overset{&Center Dot; &Center Dot;}{η}}_{r}),$ $g = I_{the s}^{- 1},$ u=T _s , f=f ₂ +gu,

$M_{{\overset{\cdot}{f}}_{1}} = \max | {\overset{\cdot}{f}}_{1} |,$ $M_{\overset{\cdot}{f}} = \max | \overset{\cdot}{f} |,$ $M_{f_{1}} = \max | f_{1} |,$ $M_{f_{2}} = \max | f_{2} |,$ M_u＝max|gu|； $m_{{\overset{&Center Dot;}{f}}_{1}} = \max | {\overset{\cdot}{f}}_{1} |,$ $m_{\overset{&Center Dot;}{f}} = \max | \overset{&Center Dot;}{f} |,$ $m_{f_{1}} = \max | f_{1} |,$ $m_{f_{2}} = \max | f_{2} |,$ M _u =max|gu|;

第二步，利用公式(3)和第一步得到的时间尺度p确定采样时间尺度h，In the second step, use the formula (3) and the time scale p obtained in the first step to determine the sampling time scale h,

$h h = = \frac{p p}{d d},, d d > > 55 - - - - - - ((33));;$

第三步，利用公式组(4)和第一步中确定的参数f₁、x₁、x₂、g，得到f_1i关于x_1j、x_2j的偏导数和g的上界M，The third step is to use the formula group (4) and the parameters f ₁ , x ₁ , x ₂ , g determined in the first step to obtain the partial derivative of f _1i with respect to x _1j , x _2j and the upper bound M of g,

$| | \frac{&PartialD; &PartialD; {f f}_{11 i i}}{&PartialD; &PartialD; {x x}_{11 j j}} | | \leq \leq M m$

$| | \frac{&PartialD; &PartialD; {f f}_{11 i i}}{&PartialD; &PartialD; {x x}_{22 j j}} | | \leq \leq M m$

||g(k)||≤M (4)||g(k)||≤M (4)

其中x₂＝w_s，i，j＝1，2，３，k＝1，2，…，f_1i、x_1i、x_2i表示f₁、x₁、x₂的第i行；Where x ₂ =w _s , i, j=1, 2, 3, k=1, 2, ..., f _1i , x _1i , x _2i represent the ith row of f ₁ , x ₁ , x ₂ ;

第四步，利用公式(5)和第三步得到的M确定参数m，In the fourth step, use formula (5) and M obtained in the third step to determine the parameter m,

$m m > > \frac{ln ln \frac{{d d}^{22} ϵ ϵ}{66 M m}}{ln ln \frac{{N N}_{x x}}{{C C}_{x x}^{22} + + {N N}_{x x}}} - - 11 - - - - - - ((55))$

其中，N_x＞0表示挠性卫星输出|y_i|即y的上界，i＝1，2，3，ε为建模误差，C_x＞0；Among them, N _x >0 means the flexible satellite output |y _i | is the upper bound of y, i=1, 2, 3, ε is the modeling error, C _x >0;

第五步，根据第二步确定的采样时间尺度h、第三步确定的M和第四步确定的参数m，利用公式组(6)得到挠性卫星特征模型的系数范围，In the fifth step, according to the sampling time scale h determined in the second step, the M determined in the third step and the parameter m determined in the fourth step, the coefficient range of the flexible satellite characteristic model is obtained by using the formula group (6),

$| | {a a}_{i i 11} ((k k)) - - 22 | | \leq \leq Mh mh + + \frac{66 M m ((m m + + 11))}{{d d}^{22}}$

$| | {a a}_{i i 22} ((k k)) + + 11 | | \leq \leq Mh mh + + \frac{66 M m ((m m + + 11))}{{d d}^{22}} - - - - - - ((66))$

|b_ij(k)|≤M²h² |b _ij (k)|≤M ² h ²

其中a_i1、a_i2、b_ij表示特征模型的系数，b_ij∈R，j＝1，2，3，b_i＝[b_i1 b_i2 b_i3]；第六步，利用梯度法对由第五步得到的特征模型的系数a_i1、a_i2、b_ij进行辨识得到辨识后的特征模型的系数

Where a _i1 , a _i2 , and b _ij represent the coefficients of the feature model, b _ij ∈ R, j=1, 2, 3, b _i =[b _i1 b _i2 b _i3 ]; the sixth step, use the gradient method to The coefficients a _i1 , a _i2 , and b _ij of the characteristic model obtained in five steps are identified to obtain the coefficients of the characteristic model after identification

第七步，利用第六步得到的辨识后的特征模型的系数

组成公式(7)的控制率，The seventh step is to use the coefficients of the identified feature model obtained in the sixth step

Constituting the control rate of equation (7),

$u u = = {[\begin{matrix} {\overset{^^}{b b}}_{11}^{T T} & {\overset{^^}{b b}}_{22}^{T T} & {\overset{^^}{b b}}_{33}^{T T} \end{matrix}]}^{- - T T} \times \times diag diag [[{u u}_{11},, {u u}_{22},, {u u}_{33}]] - - - - - - ((77))$

其中，u_i＝u_0i+u_Gi+u_Ii+u_Di，u_0i为维持/跟踪控制律，u_Gi为黄金分割控制率，u_Ii为逻辑积分控制率，u_Di为逻辑微分控制率；Among them, u _i =u _0i +u _Gi +u _Ii +u _Di , u _0i is the maintenance/tracking control law, u _Gi is the golden section control rate, u _Ii is the logic integral control rate, u _Di is the logic differential control rate;

第八步，将第七步确定的控制率代入公式组(1)的挠性卫星动力学方程中，控制挠性卫星的俯仰、偏航和滚动姿态角。In the eighth step, the control rate determined in the seventh step is substituted into the dynamic equation of the flexible satellite in formula group (1) to control the pitch, yaw and roll attitude angles of the flexible satellite.

所述第四步中C_x的确定通过以下步骤完成，The determination of _Cx in the fourth step is accomplished through the following steps,

A4.1、h²|h_i|为连续函数，且h²|h_i|≤ε，得到|x_1i|的范围，A4.1, h ² |h _i | is a continuous function, and h ² |h _i |≤ε, the range of |x _1i | is obtained,

其中

h为采样时间尺度，ε为建模误差；in

h is the sampling time scale, ε is the modeling error;

A4.2、根据A4.1得到的|x_1i|的范围，根据公式|x_1i|≤C_x确定C_x的值。A4.2. According to the range of |x _1i | obtained in A4.1, determine the value of C _x according to the formula |x _1i |≤C _x .

本发明与现有技术相比有益效果为：Compared with the prior art, the present invention has beneficial effects as follows:

(1)本发明引入了挠性卫星的时间尺度和采样周期，刻画了挠性卫星的变化率，解决了挠性卫星特征建模的瓶颈问题；(1) The present invention introduces the time scale and sampling period of the flexible satellite, describes the rate of change of the flexible satellite, and solves the bottleneck problem of flexible satellite feature modeling;

(2)本发明给出了挠性卫星特征模型参数范围的表达式，定性研究了特征模型的参数性质，从所给出的参数范围可以看出，特征模型参数的界与采样周期、建模误差、系统阶数、系统变化率有关，为挠性卫星基于特征模型的自适应控制奠定了理论基础；(2) The present invention has provided the expression of flexible satellite characteristic model parameter scope, has studied the parameter property of characteristic model qualitatively, can find out from the parameter scope provided, the boundary of characteristic model parameter and sampling period, modeling Error, system order, and system change rate are related, laying a theoretical foundation for the adaptive control of flexible satellites based on characteristic models;

(3)本发明适用于飞行器姿态动力学的特征建模，从而为飞行器基于特征模型的姿态控制奠定了基础；(3) The present invention is applicable to the characteristic modeling of aircraft attitude dynamics, thereby laid the foundation for the attitude control of aircraft based on the characteristic model;

(4)本发明的方法简单、明确，适于工程设计。(4) The method of the present invention is simple and clear, and is suitable for engineering design.

附图说明 Description of drawings

图1为本发明流程图。Fig. 1 is the flow chart of the present invention.

具体实施方式 Detailed ways

本发明如图1所示，通过以下步骤实现：As shown in Figure 1, the present invention is realized through the following steps:

1、按照挠性卫星的动力学方程，确定其时间尺度、采样时间，以及参数M和m。1. According to the dynamic equation of the flexible satellite, determine its time scale, sampling time, and parameters M and m.

挠性卫星动力学为：The flexible satellite dynamics are:

$\{\begin{matrix} {I I}_{s the s} {\overset{\cdot &Center Dot;}{w w}}_{s the s} + + {\overset{~ ~}{w w}}_{s the s} {I I}_{s the s} {w w}_{s the s} + + {F f}_{sl sl} {\overset{\cdot &Center Dot; \cdot &Center Dot;}{η η}}_{l l} + + {F f}_{sr sr} {\overset{\cdot &Center Dot; \cdot &Center Dot;}{η η}}_{r r} = = {T T}_{s the s} \\ {\overset{\cdot &Center Dot; \cdot &Center Dot;}{η η}}_{l l} + + 22 {ξ ξ}_{l l} {w w}_{l l} {\overset{\cdot &Center Dot;}{η η}}_{l l} + + {w w}_{l l}^{22} {\overset{\cdot &Center Dot;}{η η}}_{l l} + + {F f}_{sl sl}^{T T} {\overset{\cdot &Center Dot;}{w w}}_{s the s} = = 00 \\ {\overset{\cdot &Center Dot; \cdot &Center Dot;}{η η}}_{r r} + + 22 {ξ ξ}_{r r} {w w}_{r r} {\overset{\cdot \cdot}{η η}}_{r r} + + {w w}_{r r}^{22} {\overset{\cdot &Center Dot;}{η η}}_{r r} + + {F f}_{sr sr}^{T T} {\overset{\cdot \cdot}{w w}}_{s the s} = = 00 \\ y the y = = [\begin{matrix} φ φ & θ θ & ψ ψ \end{matrix}] \end{matrix}$

其中，in,

φ，θ，ψ表示挠性卫星的俯仰、偏航、滚动姿态角，[w_x w_y w_z]^T表示卫星相对轨道坐标系的角速度在体坐标系中的坐标，w_s和

分别表示卫星中心体的角速度列阵和反对称阵，η_l，η_r分别为左右太阳翼的模态坐标阵，ξ_l，ξ_r分别为左右太阳翼的模态阻尼系数，F_sl，F_sr分别为左右太阳翼与中心体的耦合系数，T_s表示作用在卫星上的外力矩列阵，I_s表示卫星惯量阵，x₁＝[φθψ]^T，y表示卫星输出，w_l、w_r分别为挠性卫星左、右太阳翼的角速度。φ, θ, ψ represent the pitch, yaw and roll attitude angles of the flexible satellite, [w _x w _y w _z ] ^T represents the coordinates of the satellite’s angular velocity relative to the orbital coordinate system in the body coordinate system, w _s and

denote the angular velocity array and the antisymmetric array of the satellite center body respectively, η _l , η _r are the modal coordinate arrays of the left and right solar wings respectively, ξ _l , ξ _r are the modal damping coefficients of the left and right solar wings respectively, F _sl , F _sr is the coupling coefficient between the left and right solar wings and the central body, T _s represents the external force matrix array acting on the satellite, I _s represents the satellite inertia array, x ₁ =[φθψ] ^T , y represents the satellite output, w _l , w _r are the angular velocities of the left and right solar wings of the flexible satellite, respectively.

时间尺度p按照下式计算：The time scale p is calculated according to the following formula:

$p p = = min min {{\frac{11}{\sqrt{{M m}_{{\overset{\cdot &Center Dot;}{f f}}_{11}}}},, \frac{11}{\sqrt{{M m}_{\overset{\cdot &Center Dot;}{f f}}}},, \frac{11}{\sqrt{{M m}_{{f f}_{11}}}},, \frac{11}{\sqrt{{M m}_{{f f}_{22}}}},, \frac{11}{\sqrt{{M m}_{u u}}}}}$

其中，in,

f₁＝C(x₁)w_s， $f_{2} = - I_{s}^{- 1} ({\tilde{w}}_{s} I_{s} w_{s} + F_{sl} {\overset{\cdot \cdot}{η}}_{l} + F_{sr} {\overset{\cdot \cdot}{η}}_{r}),$ $g = I_{s}^{- 1},$ u＝T_s，f＝f₂+guf ₁ =C(x ₁ )w _s , $f_{2} = - I_{the s}^{- 1} ({\tilde{w}}_{the s} I_{the s} w_{the s} + f_{sl} {\overset{&Center Dot; &Center Dot;}{η}}_{l} + f_{sr} {\overset{&Center Dot; &Center Dot;}{η}}_{r}),$ $g = I_{the s}^{- 1},$ u=T _s , f=f ₂ +gu

$M_{{\overset{\cdot}{f}}_{1}} = \max | {\overset{\cdot}{f}}_{1} |,$ $M_{\overset{\cdot}{f}} = \max | \overset{\cdot}{f} |,$ $M_{f_{1}} = \max | f_{1} |,$ $M_{f_{2}} = \max | f_{2} |,$ M_u＝max|gu| $m_{{\overset{&Center Dot;}{f}}_{1}} = \max | {\overset{&Center Dot;}{f}}_{1} |,$ $m_{\overset{\cdot}{f}} = \max | \overset{&Center Dot;}{f} |,$ $m_{f_{1}} = \max | f_{1} |,$ $m_{f_{2}} = \max | f_{2} |,$ M _u ＝max|gu|

采样时间尺度h按照下式计算：The sampling time scale h is calculated according to the following formula:

$h h = = \frac{p p}{d d},, d d > > 55$

M为f_1i关于x_1j、x_2j(i，j＝1，2，3)的偏导数和g的上界，即M is the partial derivative of f _1i with respect to x _1j , x _2j (i, j=1, 2, 3) and the upper bound of g, namely

$| \frac{&PartialD; f_{1 i}}{{&PartialD; x}_{1 j}} | \leq M,$ $| \frac{&PartialD; f_{1 i}}{&PartialD; x_{2 j}} | \leq M,$ ||g(k)||≤M，i，j＝1，2，3，k＝1，2，… $| \frac{&PartialD; f_{1 i}}{{&PartialD; x}_{1 j}} | \leq m,$ $| \frac{&PartialD; f_{1 i}}{&PartialD; x_{2 j}} | \leq m,$ ||g(k)||≤M, i, j=1, 2, 3, k=1, 2, ...

其中，x₂＝w_s，f_1i，x_1i，x_2i表示f₁，x₁，x₂，的第i行，i＝1，2，3。Wherein, x ₂ =w _s , f _1i , x _1i , x _2i represent the ith row of f ₁ , x ₁ , x ₂ , i=1, 2, 3.

m按照下式确定：m is determined according to the following formula:

$m m > > \frac{ln ln \frac{{d d}^{22} ϵ ϵ}{66 M m}}{ln ln \frac{{N N}_{x x}}{{C C}_{x x}^{22} + + {N N}_{x x}}} - - 11$

其中，N_x＞0表示系统输出|y_i|的上界，i＝1，2，3；Among them, N _x >0 means the upper bound of the system output |y _i |, i=1, 2, 3;

C_x＞0按照下式确定：记C _x > 0 is determined according to the following formula:

${h h}_{i i} (({x x}_{11},, {x x}_{22})) = = {Σ Σ}_{j j = = 11,, j j &NotEqual; &NotEqual; i i}^{33} \frac{&PartialD; &PartialD; {f f}_{11 i i}}{&PartialD; &PartialD; {x x}_{11 j j}} {\overset{\cdot \cdot}{x x}}_{11 j j} + + \frac{&PartialD; &PartialD; {f f}_{11 i i}}{&PartialD; &PartialD; {x x}_{22}} {f f}_{22} (({x x}_{11},, {x x}_{22})),, i i = = 1,2,3 1,2,3$

h²|h_i|为连续函数，则对于给定的建模误差ε，存在C_x＞0，对于i＝1，2，3，当|x_1i|≤C_x时，h²|h_i|≤ε。h ² |h _i | is a continuous function, then for a given modeling error ε, there exists C _x >0, for i=1, 2, 3, when |x _1i |≤C _x , h ² |h _i |≤ε.

2、根据步骤1中得到的各个变量确定特征模型的系数范围。2. Determine the coefficient range of the characteristic model according to each variable obtained in step 1.

设计一组非线性函数和压缩函数：Design a set of nonlinear and compressive functions:

${s the s}_{i i 11} ((k k)) = = \frac{{y the y}_{i i} ((k k))}{{y the y}_{i i} {((k k))}^{22} + + {y the y}_{i i} {((k k - - 11))}^{22} + + {N N}_{x x}}$

${s the s}_{i i 22} ((k k)) = = \frac{{x x}_{11 i i} ((k k - - 11))}{{y the y}_{i i} {((k k))}^{22} + + {y the y}_{i i} {((k k - - 11))}^{22} + + {N N}_{x x}}$

${s the s}_{i i 33} ((k k)) = = \frac{{N N}_{x x}}{{y the y}_{i i} {((k k))}^{22} + + {y the y}_{i i} {((k k - - 11))}^{22} + + {N N}_{x x}}$

和and

${f f}_{c c,, i i} = = {S S}_{i i} ((k k)) {s the s}_{i i 11} {y the y}_{i i} ((k k)) + + {S S}_{i i} ((k k)) {s the s}_{i i 22} {y the y}_{i i} ((k k - - 11)) + + {s the s}_{i i 33}^{m m + + 11},, i i = = 1,2,3 1,2,3$

其中， $S_{i} (k) = Σ_{j = 0}^{m} s_{i 3} {(k)}^{j}, i = 1,2,3 .$ in, $S_{i} (k) = Σ_{j = 0}^{m} {the s}_{i 3} {(k)}^{j}, i = 1,2,3 .$

特征模型为：The feature model is:

${\overset{~ ~}{y the y}}_{i i} ((k k + + 11)) = = {a a}_{i i 11} ((k k)) {\overset{~ ~}{y the y}}_{i i} ((k k)) + + {a a}_{i i 22} ((k k)) {\overset{~ ~}{y the y}}_{i i} ((k k - - 11)) + + {b b}_{i i} ((k k)) u u ((k k)),, i i = = 1,2,3 1,2,3$

其中，a_i1，a_i2∈R，b_i＝[b_i1 b_i2 b_i3]，b_ij∈R，j＝1，2，3，表示特征模型的系数，

表示特征模型的输出，R表示实数集。Among them, a _i1 , a _i2 ∈ R, b _i = [b _i1 b _i2 b _i3 ], b _ij ∈ R, j=1, 2, 3, representing the coefficients of the feature model,

Represents the output of the feature model, and R represents the set of real numbers.

下面给出y₁＝x₁₁的特征模型，其余状态类同。x₁₁满足如下微分方程：The characteristic model of y ₁ =x ₁₁ is given below, and the rest of the states are similar. x ₁₁ satisfies the following differential equation:

${\overset{\cdot \cdot}{x x}}_{1111} = = {f f}_{1111} (({x x}_{11},, {x x}_{22}))$

对上式求导，可得Deriving the above formula, we can get

${\overset{\cdot &Center Dot; \cdot &Center Dot;}{x x}}_{1111} = = \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{1111}} {\overset{\cdot &Center Dot;}{x x}}_{1111} + + {h h}_{11} (({x x}_{11},, {x x}_{22})) + + \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{22}} gu gu$

其中，in,

${h h}_{11} (({x x}_{11},, {x x}_{22})) = = {Σ Σ}_{j j = = 22}^{33} \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{11 j j}} {\overset{\cdot \cdot}{x x}}_{11 j j} + + \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{22}} {f f}_{22} (({x x}_{11},, {x x}_{22}))$

对上式进行离散化，易得Discretizing the above formula, it is easy to get

${x x}_{1111} ((k k + + 11)) = = 22 {x x}_{1111} ((k k)) - - {x x}_{1111} ((k k - - 11)) + + h h \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{1111}} (({x x}_{1111} ((k k)) - - {x x}_{1111} ((k k - - 11)))) + + {h h}^{22} \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{22}} gu gu ((k k)) + + {h h}^{22} {h h}_{11} ((k k))$

将压缩函数f_c，1乘以上式的最后一项，经整理可得Multiply the compression function f _{c, 1} by the last term of the above formula, and get

x₁₁(k+1)＝a₁₁(k)x₁₁(k)+a₁₂(k)x₁₁(k-1)+b₁(k)u(k)+e₁(k)x ₁₁ (k+1)＝a ₁₁ (k)x ₁₁ (k)+a ₁₂ (k)x ₁₁ (k-1)+b ₁ (k)u(k)+e ₁ (k)

其中，in,

${a a}_{1111} ((k k)) = = 22 + + h h \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{1111}} ((k k)) + + {h h}^{22} {S S}_{11} ((k k)) {h h}_{11} ((k k)) {s the s}_{1111} ((k k))$

${a a}_{1212} ((k k)) = = - - 11 - - h h \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{1111}} ((k k)) + + {h h}^{22} {S S}_{11} ((k k)) {h h}_{11} ((k k)) {s the s}_{1212} ((k k))$

${b b}_{11} ((k k)) = = {h h}^{22} \frac{&PartialD; &PartialD; {f f}_{1111}}{&PartialD; &PartialD; {x x}_{22}} ((k k)) g g ((k k))$

${e e}_{11} ((k k)) = = {h h}^{22} {h h}_{11} {s the s}_{1313}^{m m + + 11}$

其中e₁(k)表示特征模型的建模误差。可知建模误差|e₁(k)|＜ε。由于where e ₁ (k) represents the modeling error of the feature model. It can be seen that the modeling error |e ₁ (k)|<ε. because

|S_i(k)|≤m+1，i＝1，2，3，k＝1，2，…|S _i (k)|≤m+1, i=1, 2, 3, k=1, 2, ...

由h₁的定义，可知

并且易知s₁₁，s₁₂均小于1，因此From the definition of h ₁ , we know

And it is easy to know that s ₁₁ and s ₁₂ are both less than 1, so

$| | {a a}_{i i 22} ((k k)) + + 11 | | \leq \leq Mh mh + + \frac{66 M m ((m m + + 11))}{{d d}^{22}}$

|b_ij(k)|≤M²h²，j＝1，2，3|b _ij (k)|≤M ² h ² , j=1, 2, 3

3、根据步骤2中得到的特征模型系数范围，利用梯度法辨识特征模型的参数。(具体辨识方法见方崇智、萧德云著《过程辨识》第201-229页内容)3. According to the characteristic model coefficient range obtained in step 2, use the gradient method to identify the parameters of the characteristic model. (For specific identification methods, see pages 201-229 of "Process Identification" by Fang Chongzhi and Xiao Deyun)

θ_i＝[a_i1，a_i2，b_i]^T，i＝1，2，3θ _i =[a _i1 , a _i2 , b _i ] ^T , i=1, 2, 3

${\overset{^^}{θ θ}}_{i i} = = {[[{\overset{^^}{a a}}_{i i 11},, {\overset{^^}{a a}}_{i i 22},, {\overset{^^}{b b}}_{i i}]]}^{T T},, i i = = 1,2,3 1,2,3$

其中，λ_i1，λ_i2为可调参数，在辨识过程中，利用步骤2给定的参数范围对辨识结果进行限制。Among them, λ _i1 and λ _i2 are adjustable parameters. During the identification process, the identification result is limited by the parameter range given in step 2.

4根据步骤3中辨识得到的特征模型的系数设计控制律。4 Design the control law according to the coefficients of the characteristic model identified in step 3.

控制律为：The control law is:

$u u = = {[\begin{matrix} {\overset{^^}{b b}}_{11}^{T T} & {\overset{^^}{b b}}_{22}^{T T} & {\overset{^^}{b b}}_{33}^{T T} \end{matrix}]}^{- - T T} \times \times diag diag [[{u u}_{11},, {u u}_{22},, {u u}_{33}]]$

其中，in,

u_i＝u_0i+u_Gi+u_Ii+u_Di，i＝1，2，3u _i =u _0i +u _Gi +u _Ii +u _Di , i=1, 2, 3

${u u}_{00 i i} ((k k)) = = \frac{{y the y}_{ri the ri} ((k k)) - - {\overset{^^}{a a}}_{i i 11} ((k k)) {y the y}_{ri the ri} ((k k)) - - {\overset{^^}{a a}}_{i i 22} ((k k)) {y the y}_{ri the ri} ((k k - - 11))}{{λ λ}_{00 i i}}$

${u u}_{Gi Gi} ((k k)) = = - - \frac{{l l}_{11} {\overset{^^}{a a}}_{i i 11} e e ((k k)) + + {l l}_{22} {\overset{^^}{a a}}_{i i 22} ((k k)) e e ((k k - - 11))}{{λ λ}_{Gi Gi}}$

u_Ii(k)＝u_Ii(k-1)-k_Iie_i(k)u _Ii (k)=u _Ii (k-1)-k _Ii e _i (k)

u_Di(k)＝-k_Die_i(k)u _Di (k)＝-k _Di e _i (k)

l₁＝0.382，l₂＝0.618，e_i(k)＝y_i(k)-y_ri(k)，

k_Ii1＞＞k_Ii2＞0，

或者，k_Ii1，k_Ii2，c_Di，l_Di，λ_0i，λ_Gi为所需调整参数，y_ri(k)为跟踪目标函数。u_0i，u_Gi，u_Ii，和u_Di分别称为维持/跟踪控制律，黄金分割控制率，逻辑积分控制率和逻辑微分控制率(具体的解算方法见科学技术出版社2009出版的吴宏鑫、胡军、解永春著的《基于特征模型的智能自适应控制》，第五章)。l ₁ =0.382, l ₂ =0.618, e _i (k)=y _i (k)-y _ri (k),

k _Ii1 >> k _Ii2 > 0,

or, k _Ii1 , k _Ii2 , c _Di , l _Di , λ _0i , λ _Gi are the required adjustment parameters, and y _ri (k) is the tracking objective function. u _0i , u _Gi , u _Ii , and u _Di are respectively called the maintenance/tracking control law, the golden section control rate, the logic integral control rate and the logic differential control rate (for specific solution methods, see Wu Hongxin published by Science and Technology Press in 2009 , Hu Jun, Xie Yongchun "Intelligent Adaptive Control Based on Feature Model", Chapter V).

5、将步骤4确定的控制率反馈到公式组(1)记载的挠性卫星动力学方程中，用来控制挠性卫星的俯仰、偏航和滚动姿态角。5. Feedback the control rate determined in step 4 to the dynamic equation of the flexible satellite recorded in formula group (1) to control the pitch, yaw and roll attitude angles of the flexible satellite.

本发明未详细说明部分属本领域技术人员公知常识。Parts not described in detail in the present invention belong to the common knowledge of those skilled in the art.

Claims

1. A flexible satellite control method based on a feature model is characterized by comprising the following steps:

in the first step, the time scale p of the flexible satellite dynamic equation is determined by using the formula (2),

wherein the flexible satellite kinetic equation is a formula set (1),

{\tilde{w}}_{s} = [\begin{matrix} 0 & - w_{z} & w_{y} \\ w_{z} & 0 & - w_{x} \\ - w_{y} & w_{x} & 0 \end{matrix}]

phi, theta, psi denotes the pitch, yaw and roll attitude angles of the flexure satellite, [ wx wy wz]^TRepresenting the coordinates of the angular velocity of the satellite relative to the orbital coordinate system in a body coordinate system, w_s、

Respectively representing the angular velocity array and the anti-symmetric array, eta, of the flexible satellite central body_l、η_rModal coordinate arrays, ξ, of the left and right sun wings of a flexible satellite, respectively_l、ξ_rThe modal damping coefficients of the left and right solar wings of the flexible satellite respectively, F_sl、F_srThe coupling coefficients of the left and right solar wings and the central body of the flexible satellite, T_sRepresenting an array of external moments acting on a flexible satellite, I_sRepresenting the flexible satellite inertia matrix, x_l＝[φθψ]^TY denotes the flexible satellite output, w_l、w_rThe angular velocities of the left and right solar wings of the flexible satellite respectively,

f_l＝C(x₁)w_s，

g = I_{s}^{- 1},

u＝T_s，f＝f₂+gu，

M_u＝max|gu|；

secondly, determining a sampling time scale h by using a formula (3) and the time scale p obtained in the first step,

h = \frac{p}{d}, d > 5 - - - (3);

third step, using the formula set (4) and the parameter f determined in the first step₁、x₁、x₂G, obtaining f_1iWith respect to x_1j、x_2jAnd the upper bound M of g and the partial derivative of (c),

<math> <mrow> <mo>|</mo> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>f</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <msub> <mi>x</mi> <mrow> <mn>1</mn> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> <mo>|</mo> <mo>≤</mo> <mi>M</mi> </mrow> </math>

<math> <mrow> <mo>|</mo> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>f</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <msub> <mi>x</mi> <mrow> <mn>2</mn> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> <mo>|</mo> <mo>≤</mo> <mi>M</mi> </mrow> </math>

||g(k)||≤M (4)

wherein x₂＝w_s，i，j＝1，2，3，k＝1，2，…，f_1i、x_1i、x_2iDenotes f₁、x₁、x₂Row i of (1);

step four, determining a parameter M by using the formula (5) and the M obtained in the step three,

wherein N is_x> 0 denotes the Flexible satellite output y_iI is the upper bound of y, i is 1, 2, 3, epsilon is the modeling error, C_x＞0；

Fifthly, obtaining the coefficient range of the flexible satellite characteristic model by using a formula group (6) according to the sampling time scale h determined in the second step, the M determined in the third step and the parameter M determined in the fourth step,

|b_ij(k)|≤M²h²

wherein a is_i1、a_i2、b_ijCoefficients representing a characteristic model, b_ij∈R，j＝1，2，3，b_i＝[b_i1 b_i2 b_i3]，d＞5；

Sixthly, utilizing gradient method to make coefficient a of characteristic model obtained from the fifth step_i1、a_i2、b_ijIdentifying to obtain the coefficient of the identified characteristic model

Seventhly, using the identified coefficient of the characteristic model obtained in the sixth stepThe control rate of the formula (7) is composed,

wherein u is_i＝u_0i+u_Gi+u_Ii+u_Di，u_0iFor maintaining/tracking the control law u_GiFor golden section control rate, u_IiFor the logical integral control rate, u_DiIs a logical differential control rate;

and step eight, substituting the control rate determined in the step seven into the flexible satellite dynamics equation of the formula set (1) to control the pitch, yaw and roll attitude angles of the flexible satellite.

2. The flexible satellite control method based on the feature model according to claim 1, wherein: in the fourth step C_xThe determination of (a) is accomplished by the following steps,

A4.1、h²|h_il is a continuous function, and h²|h_iLess than or equal to epsilon to obtain | x_1iThe range of l is such that,

wherein

<math> <mrow> <msub> <mi>h</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>Σ</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>&NotEqual;</mo> <mi>i</mi> </mrow> <mn>3</mn> </munderover> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>f</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <msub> <mi>x</mi> <mrow> <mn>1</mn> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> <msub> <mover> <mi>x</mi> <mo>·</mo> </mover> <mrow> <mn>1</mn> <mi>j</mi> </mrow> </msub> <mo>+</mo> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>f</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> </mrow> </mfrac> <msub> <mi>f</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mn>1,2,3</mn> <mo>,</mo> </mrow> </math>

h is a sampling time scale, and epsilon is a modeling error;

a4.2, | x obtained according to A4.1_1iRange of | according to the formula | x_1i|≤C_xDetermination of C_xThe value of (c).