# Application of Sum of Squares Method in NonlinearH∞ Control for Satellite Attitude Maneuvers

Authors: Fanwei Meng; Dini Wang; Penghui Yang; Guanzhou Xie
Journal: Complexity
 (2019)

---

## Abstract

The Hamilton–Jacobi–Issacs (HJI) inequality is the most basic relation in nonlinearH∞ design, to which no effective analytical solution is currently available. The sum of squares (SOS) method can numerically solve nonlinear problems that are not easy to solve analytically, but it still cannot solve HJI inequalities directly. In this paper, an HJI inequality suitable for SOS is firstly derived to solve the problem of nonconvex optimization. Then, the problems of SOS in nonlinear H∞ design are analyzed in detail. Finally, a two-step iterative design method for solving nonlinear H∞ control is presented. The first step is to design an adjustable nonlinear state feedback of the gain array of the system using SOS. The second step is to solve the L2 gain of the system; the optimization problem is solved by a graphical analytical method. In the iterative design, a diagonally dominant design idea is proposed to reduce the numerical error of SOS. The nonlinear H∞ control design of a polynomial system for large satellite attitude maneuvers is taken as our example. Simulation results show that the SOS method is comparable to the LMI method used for linear systems, and it is expected to find a broad range of applications in the analysis and design of nonlinear systems.


---

## Body

## 1. Introduction

Space vehicles, underwater vehicles, and mechanical arms are often required to make large-angle and fast maneuvers during use. Due to the cross-coupling between velocities, the system exhibits nonlinear rigid body dynamics when maneuvering at large angles. Generally, the analysis of nonlinear systems is conducted by using a Lyapunov function. However, such a Lyapunov function is often difficult to construct due to lack of a systematic and convenient design method.L
2 gain control is also known as nonlinear H∞ control [1, 2]. This is because when we go to the time domain, the H∞ norm, although defined in terms of the transfer function of a linear system, it would be a L2-induced norm and thus is called L2 gain in a nonlinear system. L2 gain control refers to using L2 gain as the performance index of the system design and reducing L2 gain to a minimum. This paper studies the design of the nonlinear state feedback law with L2 gain as the performance index. Although the HJI inequality is theoretically capable of solving nonlinear H∞ control [1, 2], no effective analytical solution is currently available to solve this inequality. The relatively recent emergence of SOS [3–5] provides a new promising approach to solving nonlinear H∞ control problems. The SOS method works by using SOS polynomials to study nonlinear systems. To go beyond the nonlinear properties of the object itself and use a higher-than-quadratic Lyapunov function or design a higher-order nonlinear control law, we must study polynomials in general form. If the polynomial of the corresponding system can be arranged in SOS form, it must be nonnegative. Although it is a newly developed method, SOS has shown several unique advantages in some important applications, such as estimation of regions of attraction in nonlinear systems [6, 7], satellite attitude control under large maneuvers [8, 9], aircraft attitude control [10], nonlinear model predictive control [11], and stability analysis of time-delay systems [12, 13]. Recent years have also seen the emergence of several SOS methods for nonlinear H∞ control [14–16], but problems exist with the application of these methods. The existing SOS methods need some formulas in solving problems, so they will be reviewed later in Section 4, based on which we propose a new method for solving nonlinear H∞ control. In literature [17], SOS is used to iterate between the control and the Lyapunov function to minimize the H∞ norm. In solving the controller, the optimization problem is decomposed into a series of feasible semidefinite programming problems using a relaxation method, and then the optimal L2 gain is obtained. In literature [18], a novel SOS-based adaptive sliding mode control for polynomial systems consisting of uncertainties and input nonlinearities is proposed. Vafamand et al. [19] propose an approach that uniquely considers the stability of input saturated polynomial systems together with the nonlinear control law. Another advantage of this approach over the recently guaranteed cost controller design methods is that it can solve the proposed conditions without relying on any iterative techniques. Comparable to the LMI method used in linear systems, the SOS method can solve nonlinear problems that are difficult to solve analytically and will find a wide variety of applications in the analysis and design of nonlinear systems. For example, Vafamand and Khorshidi [20] propose a new polynomial observer and controller based on SOS. The proposed approach employs the polynomial representation and numerical SOS convex optimization technique to design a novel polynomial synchronizer for hyper (chaotic) systems. The focus of our paper is how to use the SOS method to define a nonlinear H∞ control law and then design a nonlinear H∞ control for large satellite attitude maneuvers as an example to illustrate the use of this method.


## 2. System Equations and Basic Formulas

We assume that the nonlinear system can be arranged into the following state-dependent linear-like differential equation [14]:(1)x˙=Axx+B1xw,(2)z=C1xx,where x is the state variable, w is the external input, z is the performance output, and Ax, B1x, and C1x are set as the polynomial matrices for x. Such systems can also be called polynomial nonlinear systems.For a given scalarγ≥0, if there exists(3)∫0Tzt2dt≤γ2∫0Twt2dt,∀w∈L20T,for any T≥0 and x0=0, the L2 gain of the system is said to be less than or equal to γ [1].The stability andL2 gain of nonlinear polynomial systems (1) and (2) can be determined by the following Hamilton–Jacobi inequality [1, 2]:(4)∂V∂xAxx+14γ2∂V∂xB1xB1Tx∂VT∂x+xTC1TxC1xx≤0,where Vx is the storage function, Vx≥0, and Vx0=0. If the system is observed in the zero state [1, 2], then Vx is the Lyapunov function.According to equation (4), the Schur complementary lemma can be used to easily get the bounded real lemma for the following polynomial nonlinear systems parallel to the linear system.Lemma 1 (see [15]).
If a positive definite Lyapunov functionVx can satisfy the following condition, then the system of (1) and (2) is stable and the L2 gain is less than or equal to γ:(5)∂V∂xAxx+xTC1TxC1xx∂V∂xB1xB1Tx∂VT∂x−4γ2I≤0.Equation (4) is a basic inequality corresponding to general systems (1) and (2). When considering the control problem, a control input u is added to system (1). For the design problem of H∞, the performance output z also needs to have the weight of the control input u, so as to make a proper compromise between x and u in the design, that is, between the size of the error and the size of the control quantity. Therefore, the system equation in the nonlinear H∞ control problem is(6)x˙=Axx+B1xw+B2xu,(7)z=C1xxD12xu,where B2x and D12x are also polynomial matrices.When equation (7) is considered, the third term of Hamilton–Jacobi inequality (4) is(8)C1xxD12xuTC1xxD12xu=xTC1TC1x+uTDu,D=D12TD12.Letu be state feedback. The HJI inequality for solving this state feedback can be deduced according to equations (6), (7), and (3) [1, 2]:(9)∂V∂xAxx+14∂V∂x1γ2B1xB1Tx−B2xD−1xB2Tx∂VT∂x+xTC1TxC1xx≤0,and the corresponding control law is(10)u=−12D−1xB2Tx∂VT∂x.TheB2D−1B2T term in equation (9) is obtained by substituting equation (10) into uTDu. Equation (9) is to solve the HJI inequality controlled by nonlinear H∞, and the state feedback law u can be obtained by substituting the obtained storage function Vx into equation (10). But this HJI inequality is not easy to solve. The emergence of the SOS method has made us think about the possibility of using SOS to solve this HJI inequality.


## 3. SOS Method

SOS is the polynomial of sum of squares. A polynomial is a linear combination of a finite number of monomials. For instance, a polynomial(11)px1,x2=x12+2x14+2x13x2−x12x22+5x24,is made up of 5 monomials with two variables. For a polynomial px1,…,xn≜px, if there is a polynomial f1x,…,fmx, then px can be written as a sum of squares, namely,(12)px=∑i=1mfi2x.It is clear that each SOS polynomial is nonnegative or represented aspx≥0. A set of these SOS polynomials are represented by ∑x. If a polynomial is an SOS, it is written as follows:(13)px∈∑x.SOS polynomials can also be expressed in the following special quadratic form:(14)px=zTxQzx,where Q is a semidefinite symmetric matrix, zx is a column vector composed of monomials of order less than or equal to d, and polynomial px is of order less than or equal to 2d. For instance, for the polynomial in equation (10),(15)zx≔x1x12x1x2x22,Q≔1000021−0.501000−0.505.Equation (14) shows that the solution of SOS can be reduced to a linear matrix inequality (LMI) problem. Now this problem can be addressed by using a software product named SOSTOOLS, which is downloadable off the Internet [3].


## 4. NonlinearH∞ Control

SOS is a numerical method to solve nonlinear problems that are difficult to solve analytically. Therefore, attempts have been made to use SOS to solve HJI inequalities in nonlinearH∞ control [14], but things have not been going well. The following is an analysis of the solution of the state feedback problem for the difficulties encountered in using SOS to solve HJI inequalities. This state feedback will also be an integral part of the approach proposed in this paper.For the state feedback problem, the equation of the system is set as follows:(16)x˙=Axx+Bxu,where Ax and Bx are polynomial matrices of x, x∈Rn,andu∈Rm.The Lyapunov function of the system is given as follows:(17)Vx=xTPx,where the P matrix is a constant matrix.For the system in equation (16), if the following state feedback law is adopted,(18)u=Fxx,then in solving for dV/dt,(19)∂V∂xAxx+BxFxx,there will exist the multiplication of P and F, which does not constitute a convex problem. Here, we can use a common approach in the LMI method to take the inverse matrix of P matrix in the solution [18]. This approach involves the following steps. First, the nonlinear state feedback law is given as follows:(20)u=KxPx.According to equations (16)–(20), we can obtain(21)dVdtxt=xTPAx+BxKxP+Ax+BxKxPTPx.Asymptotic stability requiresdV/dt<0, so this term in the above curly braces is required to be negative definite, that is,(22)PAx+ATxP+PBxKx+KTxBTxP<0.There should be an affine relation among the variables to solve in the inequality processed by SOS, but no such relation exists between theP and Kx matrices to solve in equation (22). This makes it difficult to obtain controller (20) by using SOS to solve the P and Kx matrices in equation (22). Based on the above analysis and the findings in literature [14], we can obtain the following Theorem 1 to solve this problem.Theorem 1.
For system (16), if there exists a symmetric positive definite matrix Q of n×n and a polynomial matrix Kx of n×n, and they satisfy(23)AxQ+QATx+BxKx+KTxBTx<0,then for all x the state feedback stability problem is solvable and the controller of the stable system is given as u=KxPx, as shown in (20), where P=Q−1.This theorem is easy to prove. Noticing thatP of equation (17) is positive definite, we set Q=P−1, so Q is also positive definite. We can get formula (23) by multiplying the left and right sides of the terms in formula (22) by a Q matrix, respectively. Theorem 1 shows that, after taking the inverse of P matrix, both the Q and Kx matrices to solve in the inequality of equation (23) are affine; they therefore become solvable by using the SOS method. After solving (23), P can be obtained from P=Q−1, and then substituted into (20) to obtain the control law.Now, let us analyze the problem of solving HJI inequality (9). For the linear-like system in equation (6), its Lyapunov function is generally taken as follows [14]:(24)Vx=xTPxx,where Px is a polynomial matrix. In order to illustrate the problems in the calculation, matrix P is tentatively defined as a constant matrix, and at the same time(25)∂V∂x=2xTP.By substituting equation (25) into equation (9), it can be seen that the HJI inequality is a Riccati inequality, not an affine one. In the LMI method, a common approach to solving a Riccati inequality is by inverting its solution to obtain a linear inequality [17]. This approach is used in literature [14] to solve the HJI inequality, which features the use of Q=P−1 to form an affine relation for equation (9) and thereby makes related problems solvable by using SOS. But for nonlinear H infinity control, the solution to HJI requires the control law u from equation (10), which, as you can see from equations (10) and (25), requires the inverse P−1 to be turned over. This leads to a problem: the inverse of a polynomial matrix can by no means still be a polynomial matrix. In other words, the approach in literature [14] is incapable of solving problems on the SOSTOOLS software. In fact, literature [14] and some subsequent literatures, for instance [16], only give examples of P matrix being a constant matrix. However, if P is a constant matrix and the input matrix B2 of the system is also set as a constant matrix, it can be seen from equations (10) and (25) that the resulting control law u is a linear control law. This brings us back to the early work of nonlinear H∞ control [2].Zheng and Wu [15] propose another iterative solution to nonlinear H∞ control by SOS. This method uses inequality (5) in the bounded real lemma (see Theorem 2).Theorem 2.
For system (26),(26)x˙=Axx+B2xu+B1xw.If Lyapunov functionV=1/2xTPx makes −S an SOS, the closed-loop system formed by (10) and (26) is stable, and the L2 gain of the closed-loop system is less than γ, where(27)S≔∂V∂xAxx+B2xux+xTC1TC1x∂V∂xB1xuTxB1Tx∂VT∂x−4γ2I0ux0−D−1≤0.Taking Lemma1 into account we can easily prove Theorem 2 by simply substituting equations (26) and (7) into equation (5).Zheng and Wu [15] rewrite inequality (27) into an iterative form as follows:(28)∂Vi∂xAxx+B2xui−1x+xTC1TC1x∂Vi∂xB1xui−1TxB1Tx∂ViT∂x−4γi2I0ui−1x0−D−1≤0.For each fixed state feedback lawui−1, inequality (28) is a convex problem, which can be solved by SOS to get Vix and γi. Here, the first (initial) Lyapunov function is taken as V0=1/2xTP0x, and the first u0 is obtained by borrowing equation (10). The Lyapunov function of the subsequent steps is taken as follows:(29)Vix=12MTxPiMx,where Mx is the column vector of a monomial with a set order. This Mx is obtained when SOS is used to solve inequality (28). Equation (10) is used to calculate ui after obtaining Vix. If the observed γi value after each iteration is decreasing, the iteration can be continued. Otherwise, it will be stopped.The primary problem with this method lies in the first step, which requires a solution forV0=1/2xTP0x. Noticing that P0 is a constant matrix, if there is a solution for V0, it means that there exists a linear state feedback solution, which is also a globally stable one. This puts a limit on the nonlinear system; otherwise, other constraints should be added to the SOS algorithm. The second problem with this method is that there is no rigorous theoretical basis for the use of equation (10) for calculating ui in each step, which is only based on a “borrowing” from the HJI inequality. Therefore, this SOS-based method of solving nonlinear H∞ control is limited and imperfect. However, at least it can obtain the nonlinear control law through the higher-order Lyapunov function of equation (29), which is a step forward from the first method that solves the HJI inequality directly. Nonlinear H∞ control can be solved by HJI inequality (9) or nonlinear bounded real inequality (27). However, as mentioned above, it is difficult to solve these two inequalities directly with the SOS method. By contrast, the affine inequality is obtained by inverting the matrix P, but the solution (equation (20)) of the state feedback problem is different from that (equation (10)) of HJI. In equation (20), the nonlinear control law can be obtained by solving Kx when P is a constant matrix; in the HJI problem, however, the SOS method is restricted by the inverse of polynomial matrix. This demonstrates that the state feedback problem is simple and can give full play to the characteristics of SOS.Based on the above understanding, a new iterative method for solving nonlinearH∞ control is proposed in this paper; each iteration is calculated in two steps. The specific design steps are as follows:Step 1. By solving the state feedback problem through equation (23), the nonlinear control law uix as shown in equation (20) is obtainedStep 2. By solving Hamilton–Jacobi inequality (4) with a graphical analytical method, the L2 gain γi of the system can be obtained (see the following example for an illustration)If theγi value in Step 2 does not meet the requirements, modify the state feedback design in Step 1 and iterate until the γi value reaches its minimum.Based on this understanding, a new iterative method for solving nonlinearH∞ control is proposed in this paper; each iteration is calculated in two steps. The first step is to solve the state feedback problem and obtain the nonlinear control law uix. The second step is to use a graphical analytic method to solve Hamilton–Jacobi inequality (4) and obtain the L2 gain γi of the system (see the following example for an illustration). If this γi value does not meet the requirements, modify the state feedback design and repeat the above calculation. The modified state feedback design here refers to the feedback gain Kx in the modified equation (20). Generally speaking, a large L2 value is usually caused by the fact that the feedback gain is not large enough. This requires an increase in the feedback gain. Another cause of a large L2 value is an excessively large feedback gain. This would make the system prone to vibration and thus lead to a large L2 value. A decrease in the feedback gain is therefore required. The decision to increase or decrease can be made based on the response curve (see Figure 1). Note that this Kx is a polynomial matrix, and each element in Kx is a (polynomial) decision variable in SOS design. Hence, as long as the value range of each element in Kx is limited, the whole state feedback can be redesigned. However, the value range of decision variables can be limited by adding an inequality constraint to algorithm [3], which is a simple task in SOS. For instance, if we want a coefficient in Kx to have a range of K6<−2, which is −2−K6>0, we can add −2−K6 to the SOS program by writing a statement as follows:(30)sosineqprog,−2−coeff_6.SOS is to solve a set of polynomial inequalities, to which some inequality constraints on a certain coefficient are added. Solving the state feedback with the gain constraint is the first step, followed by solving theL2 gain with a graphical method. If the first turn requires an increase (or decrease) in the feedback gain, we can go back and modify the constraint range in the first step and take the second turn to gradually reach the optimal performance value γ. This new method not only gives full play to the characteristics of SOS, and it is also simple and easy to use.Figure 1
Time response of angular velocityw when Φ=200∘.


## 5. NonlinearH∞ Control Example

### 5.1. Large-Angle Satellite Attitude Maneuver Model

A satellite attitude maneuver control is taken as an example to illustrate the two-step iterative method proposed in this paper. This is a nonlinear system with 6 state variables. There have not been other effective methods available to solve the nonlinearH∞ control problem of such high-dimensional systems. The first step here is to use state feedback to find the control law u. The motion equations for the satellite include a dynamic equation and a kinematic differential equation. The satellite is defined as a rigid body, and its dynamic equation is given as follows:(31)Ixw˙x+Iz−Iywywz=Tx,Iyw˙y+Ix−Izwzwx=Ty,Izw˙z+Iy−Ixwxwy=Tz,where Ix,Iy,andIz are the moment of inertia of the corresponding axis, wx,wy,andwz are the angular velocity components about the corresponding axis, and Tx,Ty,andTz are the control torques of the corresponding axis.According to Euler’s theorem, any displacement of a rigid body around a fixed point can be achieved by rotating by an angle about an axis that runs through at that point. This axis is called the instantaneous axis of rotation. In attitude control, it is called eigenaxis, which is denoted byk, k=n1n2n3T; the angel of rotation is represented by Φ. In this paper, modified Rodrigue parameters (MRPs) are used to represent the attitude of the satellite [19]. These parameters can be applied to 360° rotations (about the characteristic axis). The relationship between the MRP parameter (vector) and k,Φ is(32)σ=ktanΦ4.When MRP is adopted, the kinematic differential equation for the satellite is [19](33)σ˙=Ωσw,where(34)Ωσ=141−σ2+2σ122σ1σ2−σ32σ1σ2+σ22σ2σ1+σ31−σ2+2σ222σ2σ3−σ12σ3σ1−σ22σ3σ2+σ11−σ2+2σ32,w=wxwywzT,σ2=σ12+σ22+σ32.Equations (31) and (33) constitute the nonlinear differential equation for satellite attitude motion. If the state vector is set as x=wxwywzσ1σ2σ3T, equations (31) and (33) can be arranged into a state-dependent linear-like equation as follows:(35)x˙=Axx+Bxu,where(36)Ax=0a12000000a23000a3100000a41a52a53000a51a52a53000a61a62a63000,Bx=1/Ix0001/Iy0001/Iz000000000.Specific parameters are as follows:(37)a12=Iy−IzIxx3,a23=Iz−IxIyx1,a31=Ix−IyIzx2,a41=141+x42−x52−x62,a42=142x4x5−x6,a43=142x4x6+x5,a51=142x5x4+x6,a52=141−x42+x52−x62,a53=142x5x6−x4,a61=142x6x4−x5,a62=142x6x5+x4,a63=141−x42−x52+x62.


### 5.2. SOS Solving Control Law

The inertial matrix of the rigid body satellite is set asI=diag15,16,12.5kg⋅m2. After substituting I into equation (35), SOS can be used for state feedback design according to the idea of equation (23).Here, the state feedback law obtained in literature [20] is taken as the solution for the first step as follows:(38)u1=−2.5317x1−1.14912x2x3−3.0285x4−3.0282x4x42+x52+x62,u2=−0.61048x1x3−2.7422x2−3.2208x5−3.2204x5x42+x52+x62,u3=1.4787x1x2−2.1356x3−2.5208x6−2.5207x6x42+x52+x62.What should be clear here is that inequality (23) of the state feedback problem in this example is a polynomial matrix of 6×6. The so-called polynomial matrix means that every single variable in the matrix is a polynomial or a higher-order polynomial with 6 variables. Therefore, the solution of each coefficient, i.e., the number of decision variables in the solution process, requires a large amount of calculation, and we can only rely on SOS to solve this inequality constraint.Figures1–3 are the response curves of attitude (σ), angular velocity (w), and control input (u) under the action of this control law. The initial angular velocity is zero, and the attitude at the initial moment corresponds to a rotation by an angle of Φ=200∘ about the characteristic axis k, with k=0.5028−0.66930.5469T. The corresponding attitude parameter MRP is σ0=0.5992−0.79760.6517T.Figure 2
Time response of attitudeσ when Φ=200∘.Figure 3
Time response of control inputu when Φ=200∘.


#### 5.2.1. About Saturation Design

As shown in Figures1–3, the designed control law has good regulating performance at a large angle (200∘). However, it can be seen from equation (32) that when the required rotation angle is larger, σ will be large and the required control input will be very large. Therefore, a saturation link can be added at the output end of u, which accords with the actual situation of satellite control; the actuator has a limited torque. In this case, the saturation value is ±5 N·m. This is performed on the theoretical basis that the object equations (31) and (33) in attitude control are passive. When saturation limiting occurs, the negative feedback connection of the system is broken, which is equivalent to an open loop, while the open loop connection between two passive links is always stable. Therefore, σt is going to converge in spite of a large σ0 value. When reaching the limit value, the system will return to the continuous negative feedback working state. Figures 4–6 show the adjustment waveform when Φ=320∘ after addition of the saturation link in the above calculation example. In the case of large deviation saturation, the system still converges and all the parameters are within the normal working range.Figure 4
Time response of attitudeσ when Φ=320∘ after saturation.Figure 5
Time response of angular velocityw when Φ=320∘ after saturation.Figure 6
Time response of control inputu when Φ=320∘ after saturation.


#### 5.2.2. About Numerical Error

As shown in the expressions ofAx and Bx, the Ax and Bx arrays in this example are sparse matrices, most of which are elements equal to zero, and the matrix operation in numerical calculation introduces room for numerical errors. For B matrix, in particular, its upper half is a diagonal matrix, indicating that the three input channels are independent. This is because in the system, B and K are directly multiplied (see equations (16) and (20)). If B is a diagonal matrix, its coefficients can be put into K first. B is set as an I matrix. After calculating the solution, B can be separated from BK to get U (see equation (20)). But in the case of a large moment of inertia, such as Ix=300kg⋅m2, the coefficient 1/Ix of B is small, and the number 0.0033 is not much different from the zero value of the off-diagonal elements in B matrix. This makes the diagonal dominance of B matrix insignificant. In other words, a degree of coupling exists between the channels equivalent to B matrix in numerical operation. Although the numerical solution can be obtained, the response of the system will be very slow due to the need to eliminate the effect of such coupling on stability. Tong et al.[21] point out that SOS has a drawback of slow convergence, which is actually a problem of numerical error rather than the method itself. In this paper, we deal with this problem by classifying the coefficients of B matrix into K and taking the upper half of B matrix as I matrix, i.e., diag1,1,1. This is because B and K are directly multiplied in the system (see equations (16) and (20)). If B is a diagonal matrix, its coefficients can be put into K first, and then B is set as an I matrix. After the solution is calculated, B can be separated from BK to get u (see equation (20)). When B is taken as I matrix, the diagonal dominance of B becomes noticeable, reducing the possibility of numerical errors.


### 5.3.L2 Gain of the System Solved by the Graphical Method

After the control law is obtained, theL2 gain of the system is further solved to see whether γ meets the performance requirements. The L2 gain can be obtained by solving the bounded real lemma (equation (27)). Based on the analysis we have made earlier on solving the bounded real lemma directly, a direct solution to Hamilton–Jacobi inequality (4) is proposed. The left-hand side of inequality (4) is defined as Hamilton function H [22, 23], i.e.,(39)H≔∂V∂xAxx+14γ2∂V∂xBxBTx∂VT∂x+xTCTxCxx≤0.This Hamilton function corresponds to the bounded real Lemma1. When the feedback is added to form a closed-loop system (see equation (26)), this function becomes(40)H≔∂V∂xAxx+B2xux+14γ2∂V∂xB1xB1Tx∂VT∂x+xTC1TxC1xx≤0.Here,L2 gain is weighted only to the output variable of the system and not to the control variable u, i.e., we set D12=0 in equation (7). In the attitude control of this example, the performance output is the weight of the attitude variable σ=σ1σ2σ3T (see equation (32)), so(41)C1=000q2000000q2000000q2,where q2 is the weight coefficient of each attitude variable, which is a decision variable in the optimization problem here.The idea of this paper is to solve theL2 gain of the system after the state feedback design is complete. So we have a positive definite Lypunov function in this step. As defined in literature [20], Vx=xTPx, and the solution is its inverse matrix Q=P−1, i.e.,(42)Q=diag1430.9378,1426.6854,1429.3421,1771.9941,1771.9938,1771.9940.Since a positive definite matrix remains positive definite after being multiplied by a positive number, we take this Lyapunov functionVx as the basic function and multiply it by a coefficient kv to form a new positive definite function in formula (40), which is another decision variable in the optimization process.Now the variables in Hamilton formula (40) are x,γ,q2,andkv, which can be written as Hx,γ,q2,kv. Thus, L2 gain solution translates into an optimization problem as follows:(43)minimizeγsubject toHx,γ,q2,kv≤0.Note thatq2 in the formula is a weighting of the output. If q2=1, then the meaning of L2 gain is made clearer, so q2=1 can be taken to start optimization. In this way, the only decision variables are γ and kv. To solve this optimization problem, we need to see whether H is less than or equal to zero in equation (43). If this is a second-order system, then as long as the state variables x1 and x2 are divided into grids, the H values at each grid point will form a surface. If the highest point of the surface is less than or equal to zero, the optimal solution is obtained. However, the state variables in this example are six-dimensional, and a large amount of data are involved in numerical optimization, which is difficult to express graphically. From the perspective of a specific Hamilton function, the hyperplane expressed by equation (43) should be continuous and without any mutation point. So, we can take a trajectory that goes diagonally through the entire state space to find the optimal value. Specifically, in equation (43), the state variable of each calculation point is set to the same value, that is, x1=x2=x3=x4=x5=x6=c, so that a total of 2×104 points from −10 to 10 are selected for optimization. When kv=104, we get γmin=0.42. Figure 7 shows the Hamilton function corresponding to this γmin (solid line). As a way of checking, if this Hamilton function starts going upward when γ=0.41, then H≥0. Figure 7 shows that for this type of optimization problem, since the hyperplane of the Hamilton function is continuous, a slanting trajectory is sufficient to reflect whether H is always less than zero. It is convenient and practical to select such a trajectory for optimization.Figure 7
Illustration of the Hamilton function.To examine theL2 gain of the designed system, the x,y,andz axes of the satellite are, respectively, added with segment constant disturbances wx,wy, and wz:(44)wxt=2N⋅m,0≤t<5s,0N⋅m,t≥5s,wyt=0N⋅m,0≤t<5s,1N⋅m,t≥5s,wzt=0.Figure8 shows the response curve of attitude output (σ) under the action of these disturbance torques.Figure 8
Output response under disturbance.According towt and output response curve σt, the truncation of 2-norm z2,T/w2,T from input to output in the finite range 0,T can be calculated, as shown in Figure 9. The formulas are as follows:(45)z2,T=∫0Tzt2dt1/2=∫0Tσ12+σ22+σ32dt1/2,w2,T=∫0Tw12+w22+w32dt1/2.Figure 9
Truncation of 2-norm.Theγ value of 0.42 obtained by the above design is indicated by a dotted line in Figure 9. According to Figure 9 and equation (3), the L2 gain of the designed system is less than or equal to γ=0.42. Of course, this only applies to a particular perturbation signal because the gain will be less than or equal to 0.42 for any truncated L2 function. The L2 gain in this case is a disturbance attenuation property of the system. If the disturbance attenuation property does not meet the design requirements, a modification is required to the state feedback design, i.e., to the feedback gain Kx in equation (20), and then the L2 gain γi+1 should be calculated again.


## 5.1. Large-Angle Satellite Attitude Maneuver Model

A satellite attitude maneuver control is taken as an example to illustrate the two-step iterative method proposed in this paper. This is a nonlinear system with 6 state variables. There have not been other effective methods available to solve the nonlinearH∞ control problem of such high-dimensional systems. The first step here is to use state feedback to find the control law u. The motion equations for the satellite include a dynamic equation and a kinematic differential equation. The satellite is defined as a rigid body, and its dynamic equation is given as follows:(31)Ixw˙x+Iz−Iywywz=Tx,Iyw˙y+Ix−Izwzwx=Ty,Izw˙z+Iy−Ixwxwy=Tz,where Ix,Iy,andIz are the moment of inertia of the corresponding axis, wx,wy,andwz are the angular velocity components about the corresponding axis, and Tx,Ty,andTz are the control torques of the corresponding axis.According to Euler’s theorem, any displacement of a rigid body around a fixed point can be achieved by rotating by an angle about an axis that runs through at that point. This axis is called the instantaneous axis of rotation. In attitude control, it is called eigenaxis, which is denoted byk, k=n1n2n3T; the angel of rotation is represented by Φ. In this paper, modified Rodrigue parameters (MRPs) are used to represent the attitude of the satellite [19]. These parameters can be applied to 360° rotations (about the characteristic axis). The relationship between the MRP parameter (vector) and k,Φ is(32)σ=ktanΦ4.When MRP is adopted, the kinematic differential equation for the satellite is [19](33)σ˙=Ωσw,where(34)Ωσ=141−σ2+2σ122σ1σ2−σ32σ1σ2+σ22σ2σ1+σ31−σ2+2σ222σ2σ3−σ12σ3σ1−σ22σ3σ2+σ11−σ2+2σ32,w=wxwywzT,σ2=σ12+σ22+σ32.Equations (31) and (33) constitute the nonlinear differential equation for satellite attitude motion. If the state vector is set as x=wxwywzσ1σ2σ3T, equations (31) and (33) can be arranged into a state-dependent linear-like equation as follows:(35)x˙=Axx+Bxu,where(36)Ax=0a12000000a23000a3100000a41a52a53000a51a52a53000a61a62a63000,Bx=1/Ix0001/Iy0001/Iz000000000.Specific parameters are as follows:(37)a12=Iy−IzIxx3,a23=Iz−IxIyx1,a31=Ix−IyIzx2,a41=141+x42−x52−x62,a42=142x4x5−x6,a43=142x4x6+x5,a51=142x5x4+x6,a52=141−x42+x52−x62,a53=142x5x6−x4,a61=142x6x4−x5,a62=142x6x5+x4,a63=141−x42−x52+x62.


## 5.2. SOS Solving Control Law

The inertial matrix of the rigid body satellite is set asI=diag15,16,12.5kg⋅m2. After substituting I into equation (35), SOS can be used for state feedback design according to the idea of equation (23).Here, the state feedback law obtained in literature [20] is taken as the solution for the first step as follows:(38)u1=−2.5317x1−1.14912x2x3−3.0285x4−3.0282x4x42+x52+x62,u2=−0.61048x1x3−2.7422x2−3.2208x5−3.2204x5x42+x52+x62,u3=1.4787x1x2−2.1356x3−2.5208x6−2.5207x6x42+x52+x62.What should be clear here is that inequality (23) of the state feedback problem in this example is a polynomial matrix of 6×6. The so-called polynomial matrix means that every single variable in the matrix is a polynomial or a higher-order polynomial with 6 variables. Therefore, the solution of each coefficient, i.e., the number of decision variables in the solution process, requires a large amount of calculation, and we can only rely on SOS to solve this inequality constraint.Figures1–3 are the response curves of attitude (σ), angular velocity (w), and control input (u) under the action of this control law. The initial angular velocity is zero, and the attitude at the initial moment corresponds to a rotation by an angle of Φ=200∘ about the characteristic axis k, with k=0.5028−0.66930.5469T. The corresponding attitude parameter MRP is σ0=0.5992−0.79760.6517T.Figure 2
Time response of attitudeσ when Φ=200∘.Figure 3
Time response of control inputu when Φ=200∘.


### 5.2.1. About Saturation Design

As shown in Figures1–3, the designed control law has good regulating performance at a large angle (200∘). However, it can be seen from equation (32) that when the required rotation angle is larger, σ will be large and the required control input will be very large. Therefore, a saturation link can be added at the output end of u, which accords with the actual situation of satellite control; the actuator has a limited torque. In this case, the saturation value is ±5 N·m. This is performed on the theoretical basis that the object equations (31) and (33) in attitude control are passive. When saturation limiting occurs, the negative feedback connection of the system is broken, which is equivalent to an open loop, while the open loop connection between two passive links is always stable. Therefore, σt is going to converge in spite of a large σ0 value. When reaching the limit value, the system will return to the continuous negative feedback working state. Figures 4–6 show the adjustment waveform when Φ=320∘ after addition of the saturation link in the above calculation example. In the case of large deviation saturation, the system still converges and all the parameters are within the normal working range.Figure 4
Time response of attitudeσ when Φ=320∘ after saturation.Figure 5
Time response of angular velocityw when Φ=320∘ after saturation.Figure 6
Time response of control inputu when Φ=320∘ after saturation.


### 5.2.2. About Numerical Error

As shown in the expressions ofAx and Bx, the Ax and Bx arrays in this example are sparse matrices, most of which are elements equal to zero, and the matrix operation in numerical calculation introduces room for numerical errors. For B matrix, in particular, its upper half is a diagonal matrix, indicating that the three input channels are independent. This is because in the system, B and K are directly multiplied (see equations (16) and (20)). If B is a diagonal matrix, its coefficients can be put into K first. B is set as an I matrix. After calculating the solution, B can be separated from BK to get U (see equation (20)). But in the case of a large moment of inertia, such as Ix=300kg⋅m2, the coefficient 1/Ix of B is small, and the number 0.0033 is not much different from the zero value of the off-diagonal elements in B matrix. This makes the diagonal dominance of B matrix insignificant. In other words, a degree of coupling exists between the channels equivalent to B matrix in numerical operation. Although the numerical solution can be obtained, the response of the system will be very slow due to the need to eliminate the effect of such coupling on stability. Tong et al.[21] point out that SOS has a drawback of slow convergence, which is actually a problem of numerical error rather than the method itself. In this paper, we deal with this problem by classifying the coefficients of B matrix into K and taking the upper half of B matrix as I matrix, i.e., diag1,1,1. This is because B and K are directly multiplied in the system (see equations (16) and (20)). If B is a diagonal matrix, its coefficients can be put into K first, and then B is set as an I matrix. After the solution is calculated, B can be separated from BK to get u (see equation (20)). When B is taken as I matrix, the diagonal dominance of B becomes noticeable, reducing the possibility of numerical errors.


## 5.2.1. About Saturation Design

As shown in Figures1–3, the designed control law has good regulating performance at a large angle (200∘). However, it can be seen from equation (32) that when the required rotation angle is larger, σ will be large and the required control input will be very large. Therefore, a saturation link can be added at the output end of u, which accords with the actual situation of satellite control; the actuator has a limited torque. In this case, the saturation value is ±5 N·m. This is performed on the theoretical basis that the object equations (31) and (33) in attitude control are passive. When saturation limiting occurs, the negative feedback connection of the system is broken, which is equivalent to an open loop, while the open loop connection between two passive links is always stable. Therefore, σt is going to converge in spite of a large σ0 value. When reaching the limit value, the system will return to the continuous negative feedback working state. Figures 4–6 show the adjustment waveform when Φ=320∘ after addition of the saturation link in the above calculation example. In the case of large deviation saturation, the system still converges and all the parameters are within the normal working range.Figure 4
Time response of attitudeσ when Φ=320∘ after saturation.Figure 5
Time response of angular velocityw when Φ=320∘ after saturation.Figure 6
Time response of control inputu when Φ=320∘ after saturation.


## 5.2.2. About Numerical Error

As shown in the expressions ofAx and Bx, the Ax and Bx arrays in this example are sparse matrices, most of which are elements equal to zero, and the matrix operation in numerical calculation introduces room for numerical errors. For B matrix, in particular, its upper half is a diagonal matrix, indicating that the three input channels are independent. This is because in the system, B and K are directly multiplied (see equations (16) and (20)). If B is a diagonal matrix, its coefficients can be put into K first. B is set as an I matrix. After calculating the solution, B can be separated from BK to get U (see equation (20)). But in the case of a large moment of inertia, such as Ix=300kg⋅m2, the coefficient 1/Ix of B is small, and the number 0.0033 is not much different from the zero value of the off-diagonal elements in B matrix. This makes the diagonal dominance of B matrix insignificant. In other words, a degree of coupling exists between the channels equivalent to B matrix in numerical operation. Although the numerical solution can be obtained, the response of the system will be very slow due to the need to eliminate the effect of such coupling on stability. Tong et al.[21] point out that SOS has a drawback of slow convergence, which is actually a problem of numerical error rather than the method itself. In this paper, we deal with this problem by classifying the coefficients of B matrix into K and taking the upper half of B matrix as I matrix, i.e., diag1,1,1. This is because B and K are directly multiplied in the system (see equations (16) and (20)). If B is a diagonal matrix, its coefficients can be put into K first, and then B is set as an I matrix. After the solution is calculated, B can be separated from BK to get u (see equation (20)). When B is taken as I matrix, the diagonal dominance of B becomes noticeable, reducing the possibility of numerical errors.


## 5.3.L2 Gain of the System Solved by the Graphical Method

After the control law is obtained, theL2 gain of the system is further solved to see whether γ meets the performance requirements. The L2 gain can be obtained by solving the bounded real lemma (equation (27)). Based on the analysis we have made earlier on solving the bounded real lemma directly, a direct solution to Hamilton–Jacobi inequality (4) is proposed. The left-hand side of inequality (4) is defined as Hamilton function H [22, 23], i.e.,(39)H≔∂V∂xAxx+14γ2∂V∂xBxBTx∂VT∂x+xTCTxCxx≤0.This Hamilton function corresponds to the bounded real Lemma1. When the feedback is added to form a closed-loop system (see equation (26)), this function becomes(40)H≔∂V∂xAxx+B2xux+14γ2∂V∂xB1xB1Tx∂VT∂x+xTC1TxC1xx≤0.Here,L2 gain is weighted only to the output variable of the system and not to the control variable u, i.e., we set D12=0 in equation (7). In the attitude control of this example, the performance output is the weight of the attitude variable σ=σ1σ2σ3T (see equation (32)), so(41)C1=000q2000000q2000000q2,where q2 is the weight coefficient of each attitude variable, which is a decision variable in the optimization problem here.The idea of this paper is to solve theL2 gain of the system after the state feedback design is complete. So we have a positive definite Lypunov function in this step. As defined in literature [20], Vx=xTPx, and the solution is its inverse matrix Q=P−1, i.e.,(42)Q=diag1430.9378,1426.6854,1429.3421,1771.9941,1771.9938,1771.9940.Since a positive definite matrix remains positive definite after being multiplied by a positive number, we take this Lyapunov functionVx as the basic function and multiply it by a coefficient kv to form a new positive definite function in formula (40), which is another decision variable in the optimization process.Now the variables in Hamilton formula (40) are x,γ,q2,andkv, which can be written as Hx,γ,q2,kv. Thus, L2 gain solution translates into an optimization problem as follows:(43)minimizeγsubject toHx,γ,q2,kv≤0.Note thatq2 in the formula is a weighting of the output. If q2=1, then the meaning of L2 gain is made clearer, so q2=1 can be taken to start optimization. In this way, the only decision variables are γ and kv. To solve this optimization problem, we need to see whether H is less than or equal to zero in equation (43). If this is a second-order system, then as long as the state variables x1 and x2 are divided into grids, the H values at each grid point will form a surface. If the highest point of the surface is less than or equal to zero, the optimal solution is obtained. However, the state variables in this example are six-dimensional, and a large amount of data are involved in numerical optimization, which is difficult to express graphically. From the perspective of a specific Hamilton function, the hyperplane expressed by equation (43) should be continuous and without any mutation point. So, we can take a trajectory that goes diagonally through the entire state space to find the optimal value. Specifically, in equation (43), the state variable of each calculation point is set to the same value, that is, x1=x2=x3=x4=x5=x6=c, so that a total of 2×104 points from −10 to 10 are selected for optimization. When kv=104, we get γmin=0.42. Figure 7 shows the Hamilton function corresponding to this γmin (solid line). As a way of checking, if this Hamilton function starts going upward when γ=0.41, then H≥0. Figure 7 shows that for this type of optimization problem, since the hyperplane of the Hamilton function is continuous, a slanting trajectory is sufficient to reflect whether H is always less than zero. It is convenient and practical to select such a trajectory for optimization.Figure 7
Illustration of the Hamilton function.To examine theL2 gain of the designed system, the x,y,andz axes of the satellite are, respectively, added with segment constant disturbances wx,wy, and wz:(44)wxt=2N⋅m,0≤t<5s,0N⋅m,t≥5s,wyt=0N⋅m,0≤t<5s,1N⋅m,t≥5s,wzt=0.Figure8 shows the response curve of attitude output (σ) under the action of these disturbance torques.Figure 8
Output response under disturbance.According towt and output response curve σt, the truncation of 2-norm z2,T/w2,T from input to output in the finite range 0,T can be calculated, as shown in Figure 9. The formulas are as follows:(45)z2,T=∫0Tzt2dt1/2=∫0Tσ12+σ22+σ32dt1/2,w2,T=∫0Tw12+w22+w32dt1/2.Figure 9
Truncation of 2-norm.Theγ value of 0.42 obtained by the above design is indicated by a dotted line in Figure 9. According to Figure 9 and equation (3), the L2 gain of the designed system is less than or equal to γ=0.42. Of course, this only applies to a particular perturbation signal because the gain will be less than or equal to 0.42 for any truncated L2 function. The L2 gain in this case is a disturbance attenuation property of the system. If the disturbance attenuation property does not meet the design requirements, a modification is required to the state feedback design, i.e., to the feedback gain Kx in equation (20), and then the L2 gain γi+1 should be calculated again.


## 6. Conclusions

As a numerical method, SOS can solve nonlinear problems that are difficult to solve analytically. However, this method still suffers from some limitations for its extensive use, including the inability to solve HJI inequalities or bounded real inequalities directly and at once. This paper proposes a two-step method to solve the nonlinearH∞ control problem. This method gives full play to the advantages SOS when it comes to solving the state feedback and modifying the state feedback solution, while allowing the use of a graphical analytical method to find the optimal solution for the L2 gain. Our proposed method proves to be a simple and practical solution to the nonlinear H∞ control problem of polynomial systems.



---

*Source: 5124108-2019-11-30.xml*