J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 Journal of Mechatronics, Electrical Power, and Vehicular Technology e-ISSN:2088-6985 p-ISSN: 2087-3379 www.mevjournal.com © 2015 RCEPM - LIPI All rights reserved. Open access under CC BY-NC-SA license. Accreditation Number: 633/AU/P2MI-LIPI/03/2015. doi: 10.14203/j.mev.2015.v6.9-18 A MODIFIED GAIN SCHEDULING CONTROLLER BY CONSIDERING THE SPARSENESS PROPERTY OF UAV QUADROTORS M. Qodar Abdurrohman a, *, Reka Inovan a , Ahmad Ataka a , Hilton Tnunay a , Ardhimas Wimbo a , Iswanto a , Adha Cahyadi a , Yoshio Yamamoto b a Department of Electrical Engineering and Information Technology, Gadjah Mada University Jalan Grafika 2, Yogyakarta 55281 b Department of Precision Engineering, School of Engineering, Tokai University 1117 Hiratsuka-Shi, Kanagawa-Ken, Japan Received 12 November 2014; received in revised form 17 March 2015; accepted 17 March 2015 Published online 30 July 2015 Abstract This work presented the gain scheduling based LQR for Quadrotor systems. From the original nonlinear model, the system is always controllable and observable in various equilibrium points. Moreover, the linearized systems have a unique property that is known as sparse system. Hence, in order to implement the most efficient state feedback controller, post-filter and pre-filter were introduced to transform the state coordinate to decrease coupling between states. Finally, the gain scheduling systems using these facts was proposed. The system behavior was tested using the proposed controller. The numerical studies showed the effectiveness of the controller to achieve desired altitude, attitude, and its ability during the disturbance. Keywords: quadcopters; sparse system; linearization; gain scheduling; pole-placement. I. INTRODUCTION Quadcopter is one of Unmanned Aerial Vehicles that become popular and having much attention recently, especially from the researchers and hobbyist in aeromodelling. Several factors that contribute to its popularity are its reliability in maneuvering, its ability to be flown indoors, and easier to model and control [1], [2]. One of the most important problem on the quadcopter stems comes from the fact that quadcopter is essentially not a stable system, both in stabilization and trajectory following. Therefore, special considerations are needed in designing the control system for stabilizing or maneuvering. Existing control theories in controlling quadcopter are widely varied. The most commonly used is the conventional PID control [3], mainly due to its simple structure that is easy enough to be designed and implemented in varied systems, including quadcopter [1]. The drawback of PID controller is the gain that set for optimum in some specific conditions. In order to get the better results, the controller has to be adaptive so that it can adjust the controller gain to adapt to the position and attitude change of the quadrotor. Many people have tried to design this adaptive control such as Gaikwad [4] with auto-tuning PID Loop Shaping and Liu [5] who design self-adaptive PID based on the least-square method. Another approach is proposed to control a quadcopter using PD controller equipped with active force control to reject uncertainty disturbance by estimating disturbance torque value [6]. One of the challenges in designing controller for quadrotor is that the system is non-linear, a very common to linearize first. The basic limitation of the controller design via standard linearization is the fact that the control is guaranteed to work only in the neighborhood of a single equilibrium point. Gain scheduling is a technique to design a controller of non-linear system by linearization the system at several equilibrium points, designing the controller at each point, and implementing the family of linear controllers as a single controller with varying gain or parameter [7]. This paper also present the linearization of the simplified model of quadrotor based on [1] using * Corresponding Author.Tel: +6281804092541 E-mail: m.qodar.a@mail.ugm.ac.id M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 10 the gain scheduling linearization at some equilibrium points. Then, the controllability and observability of the system for various equilibrium points is proposed. After that, using the linearized state equations, the state feedback controllers is obtained and applied to control the altitude and attitude of the quadrotor non-linear model. For implementation purpose, a controller that focused on improvement of pole placement method is proposed by restructuring the state variable of linearized quadcopter dynamics to its Jordan form to emphasize the sparseness of the dynamics. II. SYSTEM OVERVIEW A. Notation Explanation This section consists of explanation on the derivation of state-space equations of quadcopter system. The linearization is performed to formulate the transfer function of the quadcopter plant. The model is shown in Figure 1. 12 states are used for this state-space model. The position in world frame is denoted as = [x y z], while the roll, pitch, yaw angles denoted as  = []. The velocity due to x-axis, y-axis, z- axis denoted as  = [ẋ ẏ ż], and the angular velocity due to x-axis, y-axis, z-axis denoted as  = [ṗ q̇ ṙ]. The state variable and its system input is set as x = [        ], and u = [        ]. Hence it was set that x = [ x1, . . . , x12 ] T and u = [ u1, . . . , u4 ] T B. Translational and Rotational Analysis Based on Newton’s second law of translational motion, this equation is obtained: 𝐹 = 𝑚𝑣 + (𝜔 × 𝑚𝑣)  where 4 and v = 3. From Figure 1, the forces which is worked on the quadcopter is obtained as: 𝐹 = 𝐹𝑔 − 𝐹thrust  𝐹 = 0 0 𝑚𝑔 𝑇 +𝑤 𝑅𝐵 0 0 𝑇 𝑇  Therefore, equation 3 can be expressed as: 𝑣 = 1 𝑚 0 0 1 − 𝑅 0 0 𝑇 𝐵 𝑊 − 𝑝 𝑞 𝑟 × 𝑥 𝑦 𝑧  where m is the mass of quadcopter, T is vertical thrust of quadcopter against gravity and W RB is the rotation matrix from body-frame to world- frame or inertial-frame, where: 𝑊 𝑅𝐵 = 𝑐𝜃 𝑐𝜓 𝑠𝜙 𝑠𝜃 𝑐𝜓 − 𝑐𝜙 𝑠𝜓 𝑐𝜙 𝑠𝜃 𝑐𝜓 + 𝑠𝜙 𝑠𝜓 𝑐𝜃 𝑠𝜓 𝑠𝜙 𝑠𝜃 𝑠𝜓 + 𝑐𝜙 𝑐𝜓 𝑐𝜙 𝑠𝜃 𝑠𝜓 − 𝑠𝜙 𝑐𝜓 𝑠𝜃 𝑠𝜙 𝑐𝜃 𝑐𝜙 𝑐𝜃  𝑣 = 1 𝑚 −𝑇(𝑐𝜙 𝑠𝜃 𝑐𝜓 + 𝑠𝜙 𝑠𝜓 ) −𝑇(𝑐𝜙 𝑠𝜃 𝑠𝜓 − 𝑠𝜙 𝑐𝜓 ) 𝑚𝑔 − 𝑇𝑐𝜙 𝑐𝜃 − 𝑅 𝑞 𝑧 − 𝑟 𝑦 𝑟 𝑥 − 𝑝 𝑧 𝑝 𝑦 − 𝑞 𝑥 𝐵 𝑊  Assumed that ṗ, q̇, ṙ, ẋ, ẏ, ż equal to zero, then: 𝑥 = − 1 𝑚 𝑇(𝑐𝜙 𝑠𝜃 𝑐𝜓 + 𝑠𝜙 𝑠𝜓 )  𝑦 = − 1 𝑚 𝑇(𝑐𝜙 𝑠𝜃 𝑠𝜓 − 𝑠𝜙 𝑐𝜓 )  𝑧 = 𝑔 − 1 𝑚 𝑇𝑐𝜙 𝑐𝜃  Using rigid body rotational law, 𝛤 is: 𝛤 = 𝐼𝜔 + (𝜔 × 𝐼𝜔)   where, I is moment of inertia of quadcopter as: 𝐼 = 𝐼𝑥 0 0 0 𝐼𝑦 0 0 0 𝐼𝑧   Then, equation 10 can be written as: 𝐼𝜔 = 𝜏𝑥 𝜏𝑦 𝜏𝑧 − 𝑝 𝑞 𝑟 × 𝐼𝑥 𝑝 𝐼𝑦 𝑞 𝐼𝑧𝑟  and 𝜏𝑥 = 𝑑𝑏 𝜔4 2 − 𝜔2 2  𝜏𝑦 = 𝑑𝑏 𝜔1 2 − 𝜔3 2  𝜏𝑧 = 𝑘 𝜔1 2 − 𝜔2 2 + 𝜔3 2−𝜔4 2  𝑇 = 𝑏(𝜔1 2 + 𝜔2 2 + 𝜔3 2 +𝜔4 2 )  Thus, these equations for angular acceleration of quadcopter is: 𝑝 = 𝑑𝑏 𝐼𝑥 (𝜔4 2 − 𝜔2 2 ) − 𝐼𝑧 −𝐼𝑦 𝐼𝑥 𝑞 𝑟   𝑞 = 𝑑𝑏 𝐼𝑦 (𝜔1 2 − 𝜔3 2 ) − 𝐼𝑥 −𝐼𝑧 𝐼𝑦 𝑝 𝑟   𝑟 = 𝑘 𝐼𝑧 (𝜔1 2 − 𝜔2 2 + 𝜔3 2−𝜔4 2 ) − 𝐼𝑦 −𝐼𝑥 𝐼𝑧 𝑝 𝑞   Figure 1. Quadcopter axis M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 11 The roll, pitch and yaw (RPY) rates which is a function of angular velocity were derived using inverted Jacobian matrix, denoted as: 𝑊−1 = 1 𝑐𝜃 𝑐𝜃 𝑠𝜙 𝑠𝜃 𝑐𝜙 𝑠𝜃 0 𝑐𝜙 𝑐𝜃 −𝑠𝜙 𝑐𝜃 0 𝑠𝜙 𝑐𝜙  and the relation between RPY rates and angular velocity is expressed by matrix: 𝜙 𝜃 𝜓 = 𝑊−1 𝑝 𝑞 𝑟  So, the roll, pitch, yaw rates: 𝜙 = 𝑝 + 𝑠𝜙 𝑡𝜃 𝑞 + 𝑐𝜙 𝑡𝜃 𝑟 = 0  𝜃 = 𝑐𝜙 𝑞 + 𝑠𝜙 𝑟 = 0  ψ = sϕ cθ q + cϕ cθ r = 0   C. Linearization The non-linear model of quadrotor will be linearized at the equilibrium points to make the system more amenable. In order to do the linearization, it is required to find the equilibrium point of the system, hence 0 = 𝑓(𝜒 𝛼, 𝛽, 𝛾, 𝛿 ). It is trivial that part of equibrium points are as follow: 𝑥1 = 𝑥 = 𝑥7 = 0  𝑥2 = 𝑥 = 𝑥8 = 0  𝑥3 = 𝑥 = 𝑥9 = 0  𝑥 4 = 𝜙 = 𝑥10 + 𝑠𝑥4 𝑡𝑥5 𝑥11 + 𝑐𝑥4 𝑡𝑥5 𝑥12 = 0  𝑥 5 = 𝜃 = 𝑐𝑥4 𝑥11 + 𝑠𝑥4 𝑥12 = 0  𝑥 6 = 𝜓 = 𝑠𝑥 4 𝑐𝑥 5 𝑥11 + 𝑐𝑥 4 𝑐𝑥 5 𝑥12 = 0  𝑥 7 = 𝑥 = − 1 𝑚 𝑇(𝑐𝑥4 𝑠𝑥5 𝑐𝑥6 + 𝑠𝑥4 𝑠𝑥6 ) = 0  𝑥 8 = 𝑦 = − 1 𝑚 𝑇(𝑐𝑥4 𝑠𝑥5 𝑠𝑥6 − 𝑠𝑥4 𝑐𝑥6 ) = 0  𝑥 9 = 𝑧 = 𝑔 − 1 𝑚 𝑇(𝑐𝑥4 𝑐𝑥5 ) = 0  𝑥 10 = 𝑝 = 𝑑𝑏 𝐼𝑥 𝑢4 − 𝑢2 − 𝐼𝑧 −𝐼𝑦 𝐼𝑥 𝑥11 𝑥12 = 0  𝑥 11 = 𝑞 = 𝑑𝑏 𝐼𝑦 (𝑢1 − 𝑢3 ) − 𝐼𝑥 −𝐼𝑧 𝐼𝑦 𝑥10 𝑥12  𝑥 12 = 𝑟 = 𝑘 𝐼𝑧 (𝑢1 − 𝑢2 + 𝑢3 − 𝑢4 ) − 𝐼𝑦 −𝐼𝑥 𝐼𝑧 𝑥10 𝑥11  From equation (31) and (32), it was obtained: 𝑠2𝑥5 + 𝑡 2𝑥4 = 0  The only solution for this equation is x4 = 0 and x5 = 0. By combining the results with equation (28), (29), and (30) following state can be obtained as x10 = 0, x11 = 0 and x12 = 0. Using assumption that the equilibrium point is located at certain positions in Cartesian coordinate (x, y, z) and at some yaw angle positions is defined by x = y = , z = and = . So, the complete list of the state variables value in this equilibrium, X(), can be written as x1 = α, x2 = β, x3 = γ, x4 =0, x5 = 0, x6 = δ, x7= 0, x8 = 0, x9 = 0, x10 = 0, x11= 0, and x12= 0. Remark 1. It can be noted that in the equilibrium, unless the arbitrary position and the yaw angle, all of the state are zeros. The state equation for the linearized state space model is represented by: 𝑥 = 𝐴𝑥 + 𝐵𝑢 𝑦 = 𝐶𝑥 + 𝐷𝑢.  where the matrix A and B can be found by using following equation. 𝐴 = 𝜕 𝑓1 𝜕𝑥1 |𝑋(𝛼, 𝛽, 𝛾, 𝛿) ⋯ 𝜕 𝑓1 𝜕𝑥12 |𝑋(𝛼, 𝛽, 𝛾, 𝛿) ⋮ ⋱ ⋮ 𝜕𝑓12 𝜕𝑥1 |𝑋(𝛼, 𝛽, 𝛾, 𝛿) ⋯ 𝜕𝑓12 𝜕𝑥12 |𝑋(𝛼, 𝛽, 𝛾, 𝛿)  𝐵 = 𝜕𝑓1 𝜕𝑢1 |𝑋(𝛼, 𝛽, 𝛾, 𝛿) ⋯ 𝜕𝑓1 𝜕𝑢12 |𝑋(𝛼, 𝛽, 𝛾, 𝛿) ⋮ ⋱ ⋮ 𝜕𝑓12 𝜕 𝑢1 |𝑋(𝛼, 𝛽, 𝛾, 𝛿) ⋯ 𝜕𝑓12 𝜕𝑢12 |𝑋(𝛼, 𝛽, 𝛾, 𝛿)  By careful calculation it was found that: 𝐴(12×12) = 𝑂 6×6 𝐼 6×6 𝑂 2×2 𝑁 2×2 𝑂 2×1 𝐼 6×6 𝑂 4×6   and 𝐵 12 ×4 = 𝑂 8×4 𝑀 4×4 (42) where O is zero matrix and I is identity matrix. While N and M can be defined as: 𝑁 2×2 = −𝑔𝑠𝛿 −𝑔𝑐𝛿 𝑔𝑐𝛿 𝑔𝑠𝛿  𝑀 2×2 = − 𝑏 𝑚 − 𝑏 𝑚 0 − 𝑑𝑏 𝐼𝑥 − 𝑏 𝑚 − 𝑏 𝑚 0 𝑑𝑏 𝐼𝑥 𝑑𝑏 𝐼𝑦 0 𝑘 𝐼𝑧 − 𝑘 𝐼𝑧 − 𝑑𝑏 𝐼𝑦 0 𝑘 𝐼𝑧 − 𝑘 𝐼𝑧  For simulation purpose, the constants which are going to be used is set as g = 9.81 m/s 2 , ix = 0.0820 kg.m 2 , iy = 0.0845 kg.m 2 , iz = 0.1377 kg.m 2 , b = 1.2953 x 10 -5 kg.m, d = 0.165 m, k = 1.0368 x 10 -7 kg.m 2 , m = 4.34 kg. The output of this quadrotor model is defined by the vector y = [x y z ] T , so the matrix C can be written as: M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 12 𝐶 4×12 = 𝐼 3×3 𝑂 3×9 𝑂 1×3 𝐿 1×3  where 𝐿 4×12 = [0 0 1 0 0 0 0 0 0]  Finally, the linearization model of this quadcopter is: 𝑥 = 𝐴12 ×12 𝑥 + 𝐵12 ×4𝑢  𝑦 = 𝐶12 ×1𝑥  Remark 2. It can be seen from equation (41) to (45) that most of the components of the linearized state space model are zero components. Remark 3. Here matrices A, B, and C as sparse matrices and the corresponding state equation is called by the systems with sparseness property. D. Controlability and Observability Before designing the controller or the observer, it has been checked the controllability and observability of the system. Define the controllability matrix as: 𝑃𝑐 = [𝐵 𝐴𝐵 𝐴2𝐵 ⋯ 𝐴𝑛−1𝐵]  where n defines the order of the system. In this case, the matrix can be written as: 𝑃𝑐 (12 ×48 ) = [𝐵 𝐴𝐵 𝐴 2𝐵 ⋯ 𝐴11 𝐵]  𝑃𝑐 4×12 = 𝑂 8×4 𝑀 4×4 𝑂 2×4 𝑀 4×4 𝑂 6×4 𝑂 6×4 𝑅 2×4 𝑂 4×4 𝑅 2×4 𝑂 10 ×4 𝑂 8×4 𝑂 12 ×36 (51) where 𝑅 2×4 = − 𝑔𝑐𝛿 𝑑𝑏 𝐼𝑦 𝑔𝑐𝛿 𝑑𝑏 𝐼𝑥 − 𝑔𝑠𝛿 𝑑𝑏 𝐼𝑦 − 𝑔𝑐𝛿 𝑑𝑏 𝐼𝑥 𝑔𝑐𝛿 𝑑𝑏 𝐼𝑦 − 𝑔𝑐𝛿 𝑑𝑏 𝐼𝑥 𝑔𝑠𝛿 𝑑𝑏 𝐼𝑦 𝑔𝑐𝛿 𝑑𝑏 𝐼𝑥  From the above equation is concluded the following proposition. Proposition 1. The linearized sytems with components in equation (41) to (45) are controllable regardless of the value of and  Proof. The controllability matrix equation (51) has unique configuration of the matrix element. Only prove that the rank of the controllability matrix will always be full rank is needed, i.e., rank of 12, regardless of the value of and . Using standard reduced row escelon form, the systems should be converted into a perfect triangular matrix thus the prove that the system is controllable can be concluded. By the definition, the matrix observability can be written as: 𝑃𝑜 (48 ×12 ) = 𝐶 𝐶𝐴 𝐶𝐴 2 ⋯ 𝐶𝐴11 𝑇  𝑃𝑜(48×12) = 𝜒1 4×12 𝜒2 4×12 𝜒3 4×12 𝜒4 4×12 𝑂 32×12 𝑇  where 𝜒1 4×12 = 𝐼 3×3 𝑂 3×9 𝑂 1×3 𝐿 1×9  𝜒2 4×12 = 𝑂 4×6 𝐼 3×3 𝑂 1×3 𝑂 4×2 𝑍 4×1  𝜒3 4×12 = 𝑂 4×3 𝑁 2×2 𝑂 4×7 𝑂 2×2  𝜒4 4×12 = 𝑂 2×9 𝑁 2×2 𝑂 4×1 𝑂 2×2  and 𝑍 4×1 = 0 0 0 1 𝑇  then the observability of the systems as summarized in the following preposition can also be concluded. Preposition 2. The above linearized system is also always observable regardless of the variation of the value of and . Proof. The proof is very similar with the proof of the first preposition, thus to save space, it is omited. Remark 4. What makes this result interesting is the fact that the controllability and observability depend on the equilibrium point in the beginning, which is at x = y = , z = and = Eventually, the property of matrix A and B depended only on the yaw angle value. From the above remark, the system can be controlled by more advanced adaptive controller technique such as gain scheduling without concerning the controllability and observability of the system. For gain scheduling, the controller matrix K will depend on the value of yaw angle and will be both controllable and observable for all . III. CONTROLLER IMPLEMENTATION Before discussing the gain scheduling that is going to explain in the next section, the novel method for implementing the controller is proposed. It can be observed from the derived model in the previous section that the linearized dynamics of the quadcopter are dominated by chain of integrator. The controller design strategy is based on exploiting this fact by restructuring M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 13 the state variable to its Jordan form to emphasize the sparseness of the dynamics [8]. From the restructured states, a pole placement controller is designed. Both of the designed controllers are then re-transformed to the original state space. The aim of controller design in this article is to formulate a gain matrix K for state-feedback controller, such that the eigen value of a matrix (A+BK) coincides with the desired dynamics pole. The method employed in this article is based on a notion of matrix similarity. A pair of matrix A and B are called similar if there is a similarity transformation P such that: 𝐵 = 𝑃𝐴𝑃−1  One of the interesting characteristics of similarity transforms is that the eigen value of the matrices are preserved under the similarity transformation. In this article, the transformed matrix of a matrix A is denoted as A. The similarity transform used in this article are the transformation matrix from the original state notation to a new state notation x composed of: 𝒙′ = [𝜖𝑥 𝑥 𝑥 𝛼1𝜑𝑟 + 𝛼2𝜑𝑝 𝛼1𝜑 𝑟 + 𝛼2𝜑 𝑝 𝜖𝑦 𝑦 𝑦 𝛼3𝜑𝑟 + 𝛼4𝜑𝑝 𝛼3𝜑 𝑟 + 𝛼4𝜑 𝑝 𝜖𝑧 𝑧 𝑧 𝜖𝑦 𝜑𝑦 𝜑 𝑦 ] 𝑇  The matrix used for this similarity transform can be found on the appendix. The resulting matrix from this similarity transform is: 𝐴′ = 𝐻5×5 0 0 0 0 0 𝐻5×5 0 0 0 𝐻5×5 0 0 0 0 𝐻5×5   with Hn × n is an n × n matrix such that: 𝐻𝑛 ×𝑛 = 0 𝐼𝑛 −1×𝑛−1 0 0   By changing the last row of H matrix with row vector [a1 . . . an] the H matrix become a controller canonical matrix with characteristics polynomial as: 𝐻 𝑠 = 𝑎1 + 𝑎2𝑠 2 + ⋯ + 𝑎𝑛 𝑠 𝑛 + 𝑠𝑛 +1  The other interesting property of A is that by separation principle, the poles of each H can designed without regarding the poles of other blocks. This property permits to do pole placement of each H matrix separately. For a Hn × n matrix, n number of poles (p1 . . . pn) are selected. The desired characteristics polynomial is given by: 𝐻𝑖𝑑 𝑠 = (𝑠 − 𝑝𝑖 ) 𝑛 𝑖=1 = 𝑎1 + 𝑎2𝑠 2 + ⋯ + 𝑎𝑛 𝑠 𝑛 + 𝑠𝑛 +1  The coefficients of the desired polynomial (a1 . . . an) are then inserted as [a1 . . . an] to the last row of H. Such the desired dynamics matrix of the systems is given by: 𝐴𝑖𝑑 ′ = 𝐴′ + (𝐵𝐾)′ 𝐵𝑖𝑑 ′   The Bid is the collection of H blocks last rows that control the dynamics of the systems. By simple algebraic manipulation, the desired B and K can be derived from a Bid. The K are calculated as: 𝐾 = 𝐵† (𝑃−1𝐵𝑖𝑑 ′ 𝑃)  The B † matrix is defined such that for every vector v this relation hold: 𝑣 = 𝐵†𝐵𝑣   Matrix B has a special structure as: 𝐵† = [𝑂4×8 𝐵𝑠𝑢𝑏 𝑇 𝑂4×4 ]  then matrix B † can be defined as: 𝐵† = [𝑂4×8 𝐵𝑠𝑢𝑏 −1 𝑂4×4 ]  The Bsub matrix has same structure with the Haar Wavelet Analysis matrix [9], these structures ensures that Bsub matrix are invertible as long the coefficients are not zero. A. Block Diagram of Controller In this work, the gain will be computed using LQR via gain scheduling. However, pole placement is used in this section for clarity. The naive implementation of the controller will need 48 multiplication and addition. This stems from the fact that the designed gain matrix is not sparse under the original state base. The approach in implementing the designed controller is by using a pre-filter and post-filter before the controller to transform the signal between bases. The formula of gain as given in equation (67) can be divided as follows. 𝐾 = 𝐵†𝑃−1 𝑃𝑜𝑠𝑡 −𝑓𝑖𝑙𝑡𝑒𝑟 𝐵𝑖𝑑 ′ 𝐶𝑜𝑛𝑡𝑟𝑜𝑙𝑙𝑒𝑟 𝑃 𝑃𝑟𝑒 −𝑓𝑖𝑙𝑡𝑒𝑟   The value of B † are fully parameterized by the physical construction of the quadcopter. Thus the gains of this part are static for each quadcopter type. This fact is reflected in designing the controller by only using static gain and adder for this part of post-filter. The P -1 part of the post- filter would only permute the position of control signal, thus does not need any mathematical operation to be implemented. The block diagram of the implementation is shown in Figure 2. M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 14 The pre-filter part of the implementation is composed of permutation of state structure, eight multiplication blocks and four addition blocks. In the original gain-scheduling method in the next section, the entire gain has to be calculated for each yaw value. However, by using this implementation the equivalent process is achieved by changing the coefficient { a1, a2, a3, a4 }. The block diagram of the implementation is shown in Figure 3. The implementation for the controller with desired poles in -1 is shown in the block diagram in Figure 4. B. Numerical Experiment The numerical experiment is conducted by using Quadrotor Model from Peter Corke’s Robotic Toolbox [1]. The controller is implemented using simulink block. There is two numerical experiment performed. The first experiment shows the ability of the controller to achieve desired height. The second experiment shows the ability of the controller to stabilize its height given a force impulse. Time constant for controller in Figure 5 is 0.998, it closely resembles a linear system with corresponding poles in -1, while the time constant for controller in Figure 6 is 0.644. The controlled system also differs with the linear system by the existence of overshoot in both controllers. The second experiment is conducted by giving force impulse during the period of 7s to 8s. This experiment is aimed to shows the ability of the controller to correct disturbance due to external forces. The experiment is conducted using a controller with poles in -2. The result of this experiment is shown in Figure 7. It is shown Figure 2. Block Diagram Implementation of post – filter Figure 3. Block Diagram Implementation of pre-filter Figure 4. Block Diagram Implementation of controller Figure 5. Step Response for controller with -1 poles Figure 6. Step Response for controller with -2 poles Figure 7. Controller Response to Disturbance M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 15 that the controller is able to correct the error due to force impulse. IV. CONTROL SYSTEM DESIGN Knowing the controllability and observability property that is independent of the position and yaw movement of the systems, the usage gain scheduling controller is opted. For simplicity, the gain in each of system’s structural changing due to the change of eign structure is going to be computed using LQR formalism. Moreover, in order to give more insight, hereby the controller directly was implemented into the original nonlinear system equation (7) to (23) using the series of numerical studies. A. Linear Quadratic Regulator Linear quadratic regulator is a control method using state feedback law u = Kx to minimize the cost function, defined as: 𝐽(𝑢) = 𝑥𝑇𝑄𝑥 + 𝑢𝑇𝑅𝑢 𝑑𝑡 ∞ 0  where Q is weight matrix for state energy and R is weight matrix for input energy. The matrix K can be derived from equation (73). 𝐾 = 𝑅−1 𝐵𝑇 𝑆 + 𝑁𝑇  while the matrix S is solution for the Riccati equation: 𝐴𝑇𝑆 + 𝑆𝐴 − 𝑆𝐵 + 𝑁 𝑅−1 𝐵𝑇𝑆 + 𝑁𝑇 + 𝑄 = 0 In order to solve matrix K, one have to give numerical value for  so that all numerical values of matrix A and B can be obtained. For simplicity, choose the weight matrix Q and R as Q = qI12×12, R = rI4×4, where q = 10,000,000,000 and r = 0.0000000001. A small value for weight matrix R is chosen because of the minimal energy of input signal is desired. In order to make the quadrotor able to maintain its altitude, a great amount of energy for input signal is used. However, in order to include the effect of motor saturation, a maximum and minimum boundary of the motor’s rotational speed is also used (which is square root of the input signal) according to the [1], where i;max = 1,000 rpm and i;min = 700 rpm for all i = {i : 1 ≤ i ≤ 4|i  Z }. Here, it was choosen the value for  = 0.5. The matrix K obtained from the LQR is: 𝐾 4×12 = 𝜉 × [𝑈 4×4 𝑉 4×4 𝑊 4×4 ]  where  is a constant, its value is 1.0 × 10 10 and 𝑈 4×4 = −0.6205 −0.3390 0.3390 −0.6250 −0.5000 0.0000 −0.5000 −3.8392 0.6205 0.3390 −0.3390 0.6250 −0.5000 0.0000 −0.5000 3.8392 ,  𝑉 4×4 = 3.8392 0.5000 0.0000 −0.5000 −0.9007 −0.4921 0.4921 −0.9007 −3.8392 0.5000 0.0000 −0.5000 0.9007 0.4921 −0.4921 0.9007 ,  𝑊 4×4 = −0.5000 0.0000 −0.5000 −0.7071 0.7071 0.5000 0.0000 −0.5000 −0.5000 0.0000 −0.5000 0.7071 −0.7071 0.5000 0.0000 −0.5000  Using the matrix, which has been mentioned above to produce input signal based on the state feedback law u = Kx and fed the input signal to then on-linear model of quadrotor. In this part, the result for non-linear model quadrotor attitude and altitude control using the state feedback controller that is obtained by using the linearized model property of quadrotor is presented. Figure 8 shows the quadrotor altitude with initial condition z = 0 to steady state condition z = 5. While, Figure 9, shows the quadrotor yaw angle with initial condition = 0 to steady state condition  = 0.5. As the picture shown, the state feedback controller obtained from the linearized model of the quadrotor can be concluded to works well. There is no overshoot seen and the system does not need much time to reach the stability. When changing the value of variable = 1, the result was got as seen in Figure 10 and Figure 9. Altitude yaw, for δ = 0.5 Figure 8. Altitude Z, for δ = 0.5 M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 16 Figure 11. Through this simulation, it can be concluded that proposed controller is effective. B. Observer To verify the observability of the system, linear feedback observer also known as Luenberger Observer is applied. The state space equation of the estimated state can be written as follow. 𝑥 = 𝑓 𝑥 + 𝐿(𝑦 − 𝐶𝑥 )  where L is the Luenberger matrix gain and f(x) is nonlinear function describing the quadrotor model as explained in equation (25) until (36). It can be determined that the matrix L using the pole placement strategy. The eigen value of matrix (A + LC) is designed to have value 1.5 times the poles of the controller, which is the eigen value of matrix (A + KB). To make sure that the error of the observer, the difference between the real state and estimated state, is close to zero fast enough before the states value is forced to close to the desired value by the controller. The value of matrix L is: 𝐿 12 ×4 = 𝜁 × 𝑋 4×4 𝑌 4×4 𝑍 4×4 𝑇  where 𝜁 is a constant, its value is 1.0 × 10 6 and 𝑋 4×4 = 0.5319 −0.0475 −0.0475 0.0438 0.0081 0.0005 0.0319 −0.0091 0.0080 0.0319 −0.9690 0.1999 0.1124 −0.1307 0.0788 −0.0269 , () 𝑌 4×4 = −1.4502 0.0675 0.0005 −0.0091 −0.0730 0.0128 −0.1307 0.5135 4.3321 −0.3868 −0.3867 −0.3568 0.0656 0.0042 0.2596 −0.0741 , () 𝑍 4×4 = 0.0121 0.0478 −1.0013 0.2066 0.1686 −0.1960 0.0814 −0.0278 −1.4986 0.0698 0.0008 −0.0136 −0.0754 0.0132 −0.1960 0.7702 () In Figure 12 until 15, the real states drew can be seen in blue line and the estimated states in red line are convergent. For all the states, the observer seems to works well. Through this simulation, the observability of the linearized model of quadrotor has been approved. Figure 10. Altitude Z, for δ = 1 Figure 11. Altitude yaw, for δ = 1 Figure 12. State 1 to state 3 M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 17 Figure 13. State 4 to state 6 Figure 14. State 7 to state 9 Figure 15. State 10 to state 12 M.Q. Abdurrohman et al. / J. Mechatron. Electr. Power Veh. Technol 06 (2015) 9–18 18 V. CONCLUSION In this paper, the linearized model of quadrotor simplified model is obtained and show that the system is controllable and observable regardless of the value of reference position in Cartesian coordinate as well as the reference yaw angle. It is found from the linearization that the system has mostly zero components thus can be considered as a sparse system. By rearranging into its respective Jordan form finally via similarity transformation the number of components can be reduced from 48 to 8 plus one permutation and one addition blocks. A numerical study of the implementation said that the scenario worked well. Finally, the gain scheduling controller whose gains are designed via LQR approach is proposed. The simulation results said that the proposed controller that is implemented in the original system using the sparseness property was effective. APPENDIX The similarity matrix for transforming the states is given as: 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 𝑎1 0 0 0 0 0 0 0 0 0 0 𝑎2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 𝑎1 0 0 𝑎2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 𝑎3 0 0 0 0 0 0 0 0 1 0 𝑎4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 𝑎3 0 0 𝑎4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ACKNOWLEDGEMENT The authors would like to thanks to Japan International for Cooperation Agency (JICA) for their kind support in providing the Collaborative Grant for Researches Alumnae. REFERENCES [1] P. I. Corke, Robotics, vision and control: fundamental algorithms in MATLAB, ser. Springer tracts in advanced robotics. Berlin: Springer, 2011, no. v. 73. [2] P. Pounds, et al., “Modelling and control of a quad-rotor robot,” in Australasian Conference on Robotics and Automation, Auckland, New Zealand: Australian Robotics and Automation Association Inc., 2006. [3] S. Bouabdallah, et al., “PID vs LQ control techniques applied to an indoor micro quadrotor” in Intelligent Robots and Systems (IROS 2004), vol 3, 2004, pp. 2451-2456. [4] S. Gaikwad, et al., “Auto-tuning PID using loopshaping ideas,” Proceedings of the 1999 IEEE International Conference on Control Applications, vol. 1, 1999, pp. 589-593. [5] X. Liu, et al., “Design of self-adaptive PID controller based on least square method.” 3 rd International Conference on Genetic and Evolutionary Computing, IEEE., Oct. 2009, pp. 527-529. [6] N. Tamami, et al., “Proportional Derivative Active Force Control for “x” Configuration Quadcopter”, Journal of Mechatronics, Electrical Power, and Vehicular Technology, Vol.5, pp. 67-74, 2014. DOI: 10.14203/j.mev.2014.v5.67-74. [Online]. Available: www.mevjournal.com [7] H. K. Khalil, Nonlinear systems, 2nd ed. Upper Saddle River, NJ: Prentice Hall, 1996. [8] K. Busawon, “Control design using jordan controllable canonical form,” Proceedings of the 39th IEEE Conference on Decision and Control, vol. 4, 2000, pp. 3386–3391. [9] S. G. Mallat, A wavelet tour of signal processing the Sparse way. Amsterdam; Boston: Elsevier /Academic Press, 2009.