A new stabilization algorithm for mechanical systems with underactuation degree one based on controlled Lagrangians

A new stabilization control algorithm based on controlled Lagrangian method for a class of mechanical systems with underactuation degree one is presented in this paper. Firstly, a desired controlled system with the Lagrangian structure and desired properties is constructed. By equating the underactuated system with the desired system, the matching condition and controller structure are determined. A sufficient condition for the matching condition to be held is derived, and from the sufficient condition desired kinetic energy, potential energy, gyroscopic forces and dissipative forces of the desired system can be solved explicitly. Compared with the existing matching conditions, to solve proposed sufficient condition at most two partial differential equations needs to be solved, and the rests are all algebraic equations, which is easier to solve. An algorithm to solve this sufficient condition is given in detail, and a nonlinear smooth feedback control law can be obtained to stabilize the underactuated systems. Finally, the novel control algorithm is applied to achieve almost global stability for a vertical takeoff and landing aircraft and to locally stabilize a Pendubot with two degrees of freedom at the highest point. Simulation results demonstrate effectiveness of the proposed method.


Introduction
The controlled Lagrangian (CL) method is a nonlinear controller design technique from an energy perspective. Because the CL method has great advantages and theoretical value in the control problem of underactuated mechanical systems such as inverted pendulum, acrobot, aircraft, satellite and so on, the method has been widely studied and developed by many scholars and researchers in theory and application since it was first proposed in [1]. For instance, Bloch, Leonard, and Marsden in [2] designed a controller through kinetic energy shaping to stabilize a class of underactuated system whose driving variables are symmetric group variables, and they extended the method in [2] by introducing potential energy shaping in [3], which enables the full state of mechanical system asymptotically stable. Further, they presented a constructive approach to the determination of stabilizing control laws for a class of symmetry Lagrangian systems described by the Euler-Poincaré equations in [4] and realized the stabilization of a class of discrete systems by extending the CL method in [33]. The author in [5]- [7] analyzed the feasibility of the CL method, and the author in [8]- [9] managed to increase degrees of freedom of the system to achieve the purpose of broadening application fields of the CL method and improving control performances. In order to reduce the complexity of solving matching conditions, Hamberg, Auckly, and Kapitanski in [10] [11] proposed the λ -method that makes the matching conditions simplified and innovates the construction of controlled kinetic energy through Riemann measurement. Therefore, after that, the research object was no longer limited to the system with kinetic energy symmetry. Since the λ -method was improved, a simpler matching condition vequation was obtained in [12]. Dong introduced gyroscopic forces in his doctoral thesis [13], which further enhanced λ -method and obtained more general matching condition-s, and Craig Woolsey described the design steps of the CL method after the introduction of artificial gyroscopic forces and pointed out in [14] that the artificial gyroscopic forces can increase the degree of freedom for expanding stable area of the balance point and adjusting performances of the closedloop system. In [15], the criteria for matching and stabilizing of several types of mechanical systems were given and the constraints of matching conditions were relaxed by adding force shaping.
For the asymptotic stabilization problem of general n degrees of freedom underactuated mechanical systems with input coupling, many scholars have done special research on the energy shaping method based on Lagrangian framework and Hamilton framework. Shaft in [19] used the symplectic geometry theory to study the energy shaping of mechanical systems based on the Hamilton framework and then demonstrated the role of geometric mechanics in control. The authors in [20] proposed passivity-based interconnection and damping allocation control. A set of explicit solutions of partial differential equations in matching conditions under some conditions were obtained. In the same year, the author in [21] assumed that the controlled kinetic energy has the same cyclic variable as the original kinetic energy, which makes the kinetic energy equation become a set of ordinary differential equations. Similar techniques have been described for general controlled Hamiltonian and Lagrangian systems with underactuation degree one in [22], [23] and [24], respectively. Chang D. in [25] and Blankenstein G et al. in [26] explored equivalence between CL method and IDA-PBC method and found that they are equivalent for simple mechanical systems. In order to simplify the design of controller, in [18], Li and Huo obtained a sufficient condition satisfying the matching condition by extending the results of [17] which comprised the new K-equation and one Pequation cascaded, the regular condition, and the explicit gyroscopic forces. Further, Viola, Ortega, and Banavar proved that the inhomogeneous terms in K-equations can be eliminated or simplified by modifying the target dynamics and introducing a change of coordinates in the original system in [16]. In summary, so far the main existing matching conditions consist of two sets of cascaded partial differential equations. The controller design process based on the CL method is still complicated. Therefore, it is quite necessary to develop a simpler stabilization control algorithm for underactuated mechanical systems.
In this paper, we are encouraged by [17] [18] [31] [34] to further develop the CL method for a class of mechanical systems with underactuation degree one. Main contributions of the paper are summarized as follows: 1. A feasible sufficient condition satisfying the matching condition is derived. By introducing a vector w(q) into the sufficient condition, the processes of solving the desired kinetic energy matrix and gyroscopic matrix are completely decoupled. After w(q) is first solved from a partial differential equation, the desired kinetic energy and the gyroscopic force can be obtained sequentially. In addition, more free parameters provided by introducing w(q) are helpful for solving the desired potential energy with desired properties. 2. According to whether the inertia matrix of the underactuated system is a constant matrix, two different approaches to solve the sufficient conditions are provided respectively. 3. The gyroscopic matrix is configured as a specific structure, which has fewer unknown gyroscopic force parameters so that the computational complexity of solving the gyroscopic matrix is greatly reduced. Compared with the result in [18], the number of gyroscopic force parameters to be solved is reduced from nC 2 n = n 2 (n−1) Last but not least, Compared with the methods in [2] [3] [4] [7] needed using the knowledge of differential geometry such as group theory, the method proposed in this paper only use matrix theory in the whole derivation process, which is more conducive to understanding and application by engineers and technicians.
In order to verify the effectiveness of the proposed approach, the control algorithm is applied to an almost globally stabilizing scheme for a vertical takeoff and landing aircraft with strong input coupling, and obtain a controller for a two degrees of freedom Pendubot that can be asymptotically stable at the highest point within a given range, respectively. In both cases above, we all obtain a nonlinear feedback control law to give the system desired performance and rigorously prove the asymptotic stability of the closed-loop system. The simulation experiments are carried out on these two examples to confirm the correctness of the control algorithm. The remainder of this paper is organized as follows. In Section 2, the underactuated system model and control objective are presented. A desired controlled system is designed in Section 3. In Section 4, the controller structure and matching condition are determined and a feasible sufficient condition satisfying the matching condition is derived. We solve this sufficient condition to obtain the control law in Section 5. In Section 6, the algorithm proposed in this paper is applied to a vertical take-off and landing aircraft and a two degrees of freedom Pendubot, and the effectiveness of the controller is verified by simulation experiments. Section 7 concludes this paper and proposes possible research directions in the future.

System model and control objective
Consider a n degrees of freedom (DOF) mechanical system with underactuation degree one, where n 2. Let q = [q 1 , · · · , q n ] T andq = [q 1 , · · · ,q n ] T represent its generalized coordinates and generalized velocities, respectively. Denote u = B(q)v as the non-conservative generalized force applied to the system, where v is independent control input and B(q) is input coupling matrix with rank B(q) = n − 1. The kinetic energy and potential energy of the system are expressed as E k (q,q) = 1 2q T M(q)q, where M(q) is inertia matrix satisfying M(q) = M T (q) > 0, and E p (q) : R n → R, respectively. Then, the Lagrangian of the system can be written as By the use of Euler-Lagrange equation, the system dynamical equation can be obtained as follows: Objective of this paper is to design control v such that the underactuated system (2) can be stabilized to the desired equilibrium point (q,q) = (0, 0) from any initial state (q(0),q(0)).

Desired controlled system design
In order to design a controller based on the CL method, we first construct a desired controlled system with the desired controlled kinetic energy E k (q,q) = 1 2q T M(q)q, where M(q) = M T (q) > 0 is a positive definite matrix, and the desired controlled potential energy E p (q) ≥ 0, which possesses the minimum value zero only at q = 0. Then, the Lagrangian of the desired controlled system can be written as Denotingū as the non-conservative generalized force applied to the desired controlled system and using the Euler-Lagrange equation, we can obtain the desired controlled system as follows: Design desired generalized forcē where gyroscopic matrix G(q,q) satisfies G(q,q) = −G T (q,q) (4) gives the desired closed-loop system equation Its a equilibrium point is (q,q) = (0, 0). For the desired closed-loop system (6), its mechanical is a positive definite function of q andq, so it can be used as a Lyapunov function of system (6). The derivative of E(q,q) along trajectory the system (6) can be calculated as follows: From the Lyapunov stability theorem, we know that the desired closed-loop system is stable at the equilibrium point (q,q) = (0, 0).

Determination of controller structure and matching condition
When the underactuated system (2) is equivalent to the desired system (6), the accelerationq in the two systems must be equal. From (2), we can get the acceleration of the underactuated system as The accelerationq of the desired closed-loop system (6) can be obtained as follows: Equating (8) with (9) gives Since rank B(q) = n − 1, there exists an invertible matrix The first n − 1 lines of (11) show controller structure, which can be written as and the last line of (11) gives the matching condition to be satisfied by the desired system (6): where e T n = [0, · · · , 0, 1] ∈ R n . Define vectors The matching condition (13) can be written as If M(q), E p (q), G(q,q) and D(q) can be solved from the matching condition (16), the desired controller can be obtained by substituting them into (12).

A sufficient condition for matching condition
In order to solve the matching condition (16) easily, a feasible sufficient condition satisfying the matching condition is derived firstly. It is easy to verify following results: where e i ∈ R n denotes the i-th column of identity matrix I n ∈ R n×n . Considering the definition of w T (q) and using the same method to derive (17) we have In addition, we can prove that, for any x = [x 1 , · · · , x n ] T ∈ R n , wherẽ Substituting (17)- (19) into the matching condition (16) and using (20), we get Define the matching condition (21) can be converted to A sufficient condition to satisfy the matching condition (22) is that the following three equations hold simultaneously: Furthermore, in order to realize decoupling of M(q), G(q,q), w(q) in (23), we can prove that if w n = 0, (23) is equivalent to following equations: where (26) is a partial differential equation about w(q) only, and A(q) + A T (q) [−nn] indicates that the element at the n-th row and the n-th column of matrix A(q) + A T (q) is replaced by zero. In fact, since (27) and (23) are the same except for the nth row and n-th column element, it is only necessary to prove that (26) is equivalent to a nn = 0 under the condition of w n = 0, where A(q) [a i j ]. Using the definition of w(q), A(q), the property of skew symmetric matrix G(q,q) and (27), from (26) we obtain Obviously, a nn = 0, if w n = 0.
To sum up, we conclude that a sufficient condition for the matching condition (16) to be true is that w n = 0 and (24)- (27) are satisfied simultaneously. The advantage of this sufficient condition is that it realizes decoupling of the matching condition, such that the desired E p (q), M(q), G(q,q) and D(q) can be solved separately.

Solution of matching condition
Regarding the sufficient conditions (24)- (27), we find that (26) is a partial differential equation only with respect to w(q), so we first solve (26) to obtain w(q). Then, solving the potential energy equation (24) gives E p (q). Finally, using w(q), we can design M(q) and obtain G(q,q) and D(q) by solving (27) and (25), respectively. The solution process is described as following algorithm in detail.
Step 1: According to whether M(q) is a constant matrix, we can solve w(q) by two different approaches from (26). Design where the constant matrix Γ ∈ R n×n is to be solved. Substituting (29) into (28) yields A sufficient condition to satisfy (30) is Design Γ as a symmetric matrix, (31) is equivalent to a linear equation as follows: where Γ = e T 1 Γ , · · · , e T n Γ T ∈ R n 2 . In (32), there are 1 2 n(n− 1) independent equations and 1 2 n(n + 1) unknowns. Since 1 2 n(n + 1) − 1 2 n(n − 1) = n > 0, there must exist a non-zero solution Γ for (32), and it contains n free parameters at least. In order to ensure that w n = 0 and M(q) to be designed is a positive definite matrix, the free parameters in Γ must satisfy Case 2: when M(q) is not a constant matrix, w(q) is solved from (26). Since (26) contains n unknowns, among them n − 1 unknowns can be regarded as free parameters when solving w(q). Similar to Case 1, if we want to guarantee that w n = 0 and M(q) to be designed is a positive definite matrix, the n − 1 free parameters in w(q) must satisfy Remark 1: In fact, using the definition of w(q), we know that if the designed M(q) is a positive definite matrix, it follows that This implies that p T (q)M(q)w(q) > 0 is a necessary condition for designed M(q) to be a positive definite matrix. Roughly speaking, if the free parameters in Γ or w(q) do not satisfy p T (q)M(q)w(q) > 0, we can never get a desired positive definite matrix M(q) at all. Remark 2: In order to facilitate solving the desired potential energy E p (q) with an minimum value zero only at q = 0, the free parameters in Γ or w(q) should be reserved as much as possible.
Step 2: Solve the desired potential energy E p (q). Substituting the obtained w(q) into (24) gives This is a quasilinear partial differential equation. Solve it to get the expression of E p (q). By choosing the integral constants and the free parameters in w(q) such that E p (q) satisfies which ensure that E p (q) has a unique minimum value zero at q = 0. It is worth noting that so far, w(q) has been completely determined after E p (q) is solved.
Step 3: Design the desired positive definite matrix M(q). Since p T (q)M(q)w(q) > 0, we must be able to design parameters m kk (q) > 0, k = 1, · · · , n − 1, and satisfy According to w T (q) = p T (q)M(q)M −1 (q), we design where Step 4: Solve gyroscopic matrix G(q,q) from (27). Using the definition of A(q), we have It is a linear algebraic equation with respect to gyroscopic force functions and can be solved as follows. 1) When M(q) is a constant matrix and p(q) is a constant vector, we have H(q) ≡ 0 n×n . According to (39), we can chooseG(q) = 0 n×(n−1) , which implies that it is unnecessary to introduce gyroscopic forces into the desired closed-loop system (6).
Denoting matrixD(q) = [α 1 (q), · · · , α n−1 (q)], the desired dissipative matrix D can be designed as where d 1 , · · · , d n−1 are positive numbers. It is easy to verify that the obtained D(q) shown in (44) is a symmetric semipositive definite matrix and satisfies (25). So far, the desired E p (q), M(q), G(q,q) and D(q) have been completely solved from the sufficient conditions (24)- (27). Substituting them into (12), we can get the controller expression. So above process to solve the matching condition is also an algorithm to design controller for the underactuated system (2).
Besides, we have proved that the closed-loop system is stable in Section 3. For a specific system, we can prove that if there is no non-zero trajectory in the set ofĖ(q,q) ≡ 0, then the closed-loop system is asymptotically stable based on the LaSalle invariant principle.

Vertical Takeoff and Landing Aircraft
Our first example is a vertical takeoff and landing (VTOL) aircraft moving in the vertical x − y plane as shown in Fig.1, its dynamical model is given by [27][28] [34] as follows: where (x, y) indicates its position, θ is roll angle, m is mass, J is moment of inertia, g is gravitational acceleration, ε > 0 is a parameter describing the relationship between rolling torque and lateral acceleration, B is input coupling matrix and v is control input, s (·) sin(·) and c (·) cos(·). Our objective is to design control law v with the method proposed in this paper so that the VTOL aircraft is asymptotically stabilized from any initial state to the desired equilibrium point Next, follow the steps described in Section 5 to design the controller v such that q → [0, 0, 0] T .

2-DOF Pendubot
The second example is a 2-DOF Pendubot shown in Figure 3. Its first joint (rod 1) is actuated by a driving force u 1 ∈ R and the second joint (rod 2) is unactuated (u 2 ≡ 0). Both the rod 1 and rod 2 are regarded as rigid bodies. Define q = [θ 1 , θ 2 ] T , where θ 1 is the angle between the positive direction of y-axis and rod 1, and θ 2 is the relative angle between rod 1 and rod 2, both of them are positive clockwise. m i and l i represent mass and length of the i-th (i = 1, 2) rod, respectively. Let l Ci denote the distance from centroid C i to O i and I O i denote the moment of inertia of rod i with respect to O i . As mentioned in [29]- [32], by using the Euler-Lagrange equation, we can obtain the dynamical equation of the Pendubot as where the constants The control objective is to design control law u 1 with the method described in this paper so that the 2-DOF Pendubot can be locally asymptotically stabilized from an initial state near the equilibrium point (q,q) = (0, 0) to the equilibrium point.
Next, follow the steps described in Section 5 to design the controller.
We carried out simulation studies on the 2-DOF Pendubot with MATLAB. The model parameters are introduced from [31] and shown in Table 1. The initial state is set to (q(0),q(0)) = (π/6, −π/6, 0.1, 0.1) and the controller parameters are chosen as a 0 = 2.5, d 1 = 6, p 1 = 16,w 2 = −1. The change curves of the position θ 1 , θ 2 and angular velocitẏ θ 1 ,θ 2 of the 2-DOF Pendubot, control input u 1 , and desired potential E p (q) are depicted in Fig. 4, respectively. It can be seen from Fig.4 that the state of the 2-DOF Pendubot quickly converge to the desired equilibrium point (q,q) = (0, 0) from initial values, and the closed-loop system achieves local asymptotic stability, which verifies the validity of the nonlinear smooth control law designed in this paper.

Conclusion and Future Research
In this paper, a novel stabilization control algorithm based on the CL method is proposed for a class of mechanical systems with underactuation degree one. Our main research results include that one simpler sufficient condition satisfying the matching condition is derived, and the specific solution method for this sufficient condition, the complete control law expression, and the stability proof of closed-loop system are given in detail. Compared with the existing methods of the same type, our advantage is that the obtained sufficient condition at most need to solve two quasi-linear partial differential equations, and the rest are algebraic equations, which greatly reduces the difficulty and computational complexity of the matching condition. In addition, we only use matrix theory knowledge rather than group theory in the process of deriving matching conditions and characterizing controllers, which is more conducive to the direct application and understanding of engineers and technicians. Finally, using the method given in this paper, we design a global stabilization control strategy for a VTOL aircraft and a local stabilization control law for a 2-DOF Pendubot. Simulation results demonstrate the effectiveness of our control algorithm. In the future research, the work of this paper will be extended in following three aspects. First, explore how to extend it to mechanical systems with more under-actuation degrees. Second, apply the CL method to trajectory tracking and path following problems for underactuated mechanical systems. Third, deal with under-actuated mechanical systems with model uncertainty and external disturbance using the CL method.