Passive Compliance Control of Aerial Manipulators
Abstract
This paper presents a passive compliance control for aerial manipulators to achieve stable environmental interactions. The main challenge is the absence of actuation along bodyplanar directions of the aerial vehicle which might be required during the interaction to preserve passivity. The controller proposed in this paper guarantees passivity of the manipulator through a proper choice of endeffector coordinates, and that of vehicle fuselage is guaranteed by exploiting time domain passivity technique. Simulation studies validate the proposed approach.
I Introduction
After successful achievements in unmanned aerial vehicle (UAV) studies for nonactive missions such as surveillance and remote sensing, aerial manipulation is an emerging research topic. By mounting a manipulator (or multiple manipulators) to a UAV, we can exploit further capabilities of aerial platforms. UAV equipped with a manipulator (UAVM hereinafter), for instance, allows us to extend the missions to active ones such as grasping and manipulation which may include interaction with environments.
A number of studies have been performed to achieve successful aerial manipulation in various perspectives; e.g., UAVM design [1, 2, 3, 4], intelligence [5, 6, 7], and modeling methodologies [8, 9]. In addition to these, control of UAVM is also extensively studied. For example, [10, 11, 12, 8, 13, 14] studied stability of UAVM systems without considering manipulation tasks explicitly. [15, 16, 17, 18] tackled actual aerial manipulation tasks, but the compliance interaction control was not the main scope of these studies mainly because the UAVs were not equipped with multi degrees of freedom (DoF) robotic manipulators. In contrast, the system considered in this work is equipped with a multiDoF robotic manipulator as shown in Fig. 1 (see [11] for more details). To achieve further capabilities of the multiDOF manipulators, [19] and [20] formulated Cartesian impedance control approaches for UAVM systems.
In this paper, we extend our previous work [21] in which stable flight control for UAVM was studied, however the interaction with the environment was not considered. In particular, we extend the formulation from using joint level control to using endeffector control in order to achieve compliance behavior in the task space. Realizing stable environmental interaction, however, is not trivial when the UAV is underactuated because the UAV lacks actuation along body linear , directions (see Fig. 1b).^{1}^{1}1In fact, fully actuated platforms are recently being studied [22, 23]. During the interaction, the interaction force (or torque) propagates through the manipulator to the fuselage. As a result, forces along the linear ,directions are generated. These forces can not be directly handled due to the lack of actuation and a stable interaction with the environment is not guaranteed.
The contribution of this paper is to design a passive compliance control for underactuated UAVM systems to ensure stable interactions with the passive environments. Passivity of UAV’s directional translation dynamics and UAV’s rotational dynamics is rather straightforward by virtue of the collocated actuations. This paper exploits time domain passivity approach [24, 25] in which passivity observer/passivity controller (PO/PC) is used to render the system passive. However, passivity of the compliance controller of the manipulator, which may result in forces along nonactuated directions, is not trivial. In this paper, it will be shown that the passivity is guaranteed if the compliance controller is formulated in the UAV fuselage frame (i.e., using in Fig. 1b). As a result, stability is guaranteed when the controlled UAVM interacts with passive environments.
The resulting controller has three additional benefits. First, although the UAV dynamics and the manipulator dynamics are highly coupled, they are controlled independently except for the gravity compensation. This may significantly reduce the effort and time for implementation. Second, the resulting control law is almost modelfree; the total mass and gravity vector are the only required modeling parameters. Third, the endeffector compliance control is formulated using which is expressed in the fuselage frame. For example, this can be used for visual servoing control, where the desired values are usually expressed in the same frame .
Ii Modeling of UAVManipulator Systems
Considering the following assumption, this section introduces mathematical modeling of the UAVM.
Assumption 1
This paper considers nonredundant robot manipulators and assumes that the manipulator is not in singular configurations. In particular, this paper will consider 6DOF manipulator to exploit full task space.
In Fig. 1b, , , and represent global, UAV fuselage and endeffector frames respectively. The frame is located at the center of mass (CoM) of the UAV, not UAVM. In addition, UAVM velocity vectors are defined by
(1) 
where denotes the body twist of the frame . is the generalized coordinates of the manipulator. After introducing a commonly used UAVM dynamics expressed in , the coordinate will be transformed into .
Equation of motion of the UAVM in the coordinate is given by
(2) 
with inertia matrix , Coriolis/centrifugal matrix , and gravity vector . is the RPY angle associated with the rotation matrix which represents the rotation from to . represents the control command in the body frame, and is given by
(3) 
where is UAV thrust, is the torque around the CoM of UAV, and is the joint torque of the manipulator. Note that, by expressing the dynamic model in the body frame, the first two elements of are zero. UAVM external force/moment vector is given by
(4) 
where represents the body wrench due to the environmental interaction.
The Jacobian matrix defines the relation between the velocity vector and by
(5) 
Using this, the UAVM dynamics (2) can be rewritten in coordinates as
(6) 
where and are the inertia and Coriolis/centrifugal matrices in the new coordinate system. Moreover, using
(7) 
which is valid due to Assumption 1, (6) can be rewritten as
(8) 
where are defined by
(9)  
(10) 
Note that influences body direction dynamics, and influences UAV rotational dynamics. However, influences the whole UAVM dynamics.
Later, the endeffector compliance control will be formulated using defined by
(11) 
where is the RPY angle associated with , and is the displacement vector of with respect to , represented in .
For future convenience, notations used throughout the paper are summarized as follows.

, , : RPY (rollpitchyaw) angles.

.

.


: vector pointing with respect to , represented in .

: Simplification for position in .

and


Frames are omitted for body velocities.

and .

Recall body twists and in (1).


.

: Maps Euler angle rate to the angular velocity.

.

.


Vector components: . This rule also applies to the other vectors.

Inertia matrix can be partitioned as follows.
(12) Here, the subscripts , , and represent ‘translational’, ‘rotational’, and ‘manipulator’, respectively. Similarly, represents the first three elements of .
Finally, the following assumptions on passivity of environment and UAVM are made.
Assumption 2
Manipulator may interact with the environment of which the inputoutput (I/O) pair is strictly passive.
Assumption 3
The dynamics of UAVM is subjected to drag although it is omitted in (8) for simplicity. Therefore, the UAVM dynamics is output strictly passive.
Iii Passive Compliance Control of UAVM
Iiia Control goal
This paper tackles a compliance control of aerial manipulator with the following particular scenario. (i) The UAVM system approaches a target position (freeflight), (ii) and the manipulator endeffector interacts with a passive environment while the UAV is keeping its desired position. As stated earlier, this problem is not trivial as the typical UAV lacks actuations along body  and directions. As shown in Fig. 2a, the forces acting on the endeffector will propagate through the UAVM body, and will eventually result in body  and directional forces in the UAV fuselage. Since there is no actuation to counteract these forces, stability of the resulting closedloop dynamics cannot be guaranteed.
At this point, it is interesting to note that, in (8)(10), the UAV control inputs (, ) do not directly influence the manipulator dynamics. In contrast, the manipulator control input () influences the UAV dynamics directly. Keeping these in mind, this paper takes the following control strategy.

UAV control inputs and only take care of stabilization of UAV dynamics. Regardless of the manipulator dynamics (which may have environmental interaction), the UAV tries to maintain its desired position. Time domain passivity approach (in particular, PO/PC) will be applied to passivate the UAV controller.

Manipulator control is designed to render PD compliance behavior at the endeffector. However, because this control input applies forces along UAV’s body , directions which cannot be handled, will be designed to be intrinsically passive (i.e., passive without passivation technique such as PO/PC). Later, it will be shown that the compliance controller is intrinsically passive if it is formulated using .
More specifically, UAVM controller should fulfill the following control goals.

Asymptotic stability for the freeflight.^{2}^{2}2Based on our scenario of interest, UAV positioning is a global mission to approach the target (using, e.g., GPS), and manipulation is a local mission (using, e.g., vision system on the UAV fuselage).

The desired values of UAV are given in the global frame: , , with zero derivatives. Note that desired roll and pitch angles are zero because, otherwise, the UAV will move.^{3}^{3}3 Unlike most of UAV controllers that update desired orientation in the control loop to achieve zero translational error [3, 8, 14], the desired roll and pitch angles of the proposed controller are always zero. Nonzero translational error perturbs the stable rotational dynamics, which tries to achieve zero roll and pitch angles, in such a way that the error is decreased.

The desired value of the endeffector is given in the fuselage frame: with zero derivative.


Stable interaction with passive environment.

The UAV control should satisfy passivity of I/O pairs and .

The compliance control should satisfy passivity of the I/O pair .

It should be remarked that stable interaction in this paper indicates asymptotic stability of the entire control loop including the passive environment; see also Fig. 3. Therefore, the UAVM will converge to an equilibrium point at which it can balance the interaction forces, as shown in Fig. 2b. We would like to emphasize that the UAVM automatically converges to this equilibrium point without calculating it.
IiiB Control design
This section presents a passive and stable compliance controller for UAVM, as an extension of our previous study [21]. The proposed controller in this paper differs from the one in [21] in the following aspects:

Formulation is extended from joint space tracking control to endeffector compliance control.

UAV and manipulator are controlled independently.

The main interest of [21] was stability in freeflight. This paper additionally tackles stable environmental interaction.

The position controller of the UAV (fuselage) includes an integralterm in order to maintain its position against interaction forces.
In the following, the UAV controller and manipulator controller will be presented.
IiiB1 UAV control
The UAV control law is given by:
(13)  
(14) 
where , are control gains for UAV control. and , which will be defined shortly, are activated to passivate the UAV dynamics only when the passivity condition is violated [24].
In (14), is obtained by integrating the reference acceleration defined by
(15) 
with
(16)  
(17) 
Here, are diagonal gain matrices. Integral action is included in (17) to hold UAV’s desired position against interaction forces. Since the integrator is a nonpassive element, it may result in energy generation. In this paper, PO/PC techniques will be exploited to ensure passivity, and therefore, nonpassive actions (e.g., integral) will be corrected if needed.
The UAV position error perturbs the rotational dynamics via in (15).^{4}^{4}4Note that is given by . Here, is the total mass, is the skewsymmetric operator, and is the position of the CoM of the overall UAVM system from the origin of . Because the UAV , positions are indirectly controlled by orientation of the UAV, gains for rotational dynamics (, ) should be selected larger than those for translational dynamics (, ). Otherwise, the perturbation will be too large, and may result in insufficient performance. For more details on the gain selection strategy, one may refer to Section IVA in [21].
Note that, in the UAV controller (13)(14), manipulator variable does not appear explicitly, except for the gravity compensation (i.e., controlled independently). Conceptually, leads to and to . is indirectly controlled by the reference acceleration tracking of which is designed to achieve and that imply stable error dynamics:
(18)  
(19) 
because and .
On the other hand, from the energy point of view, passivity of (13)(14) is not guaranteed. To overcome this, the following PO is applied to check if the passivity is violated:
(20)  
(21) 
If the passivity condition is violated, the time varying damping terms and are applied to guarantee the passivity of UAV dynamics:
(22)  
(23) 
where is the sampling time and , are initially stored energy.
IiiB2 Manipulator control (compliance control)
The manipulator control input to realize compliance behavior of the endeffector is given by
(24) 
where are stiffness and damping gains respectively, and . Similar to the UAV control, the UAV variables do not appear in the manipulator control except for the gravity compensation.
To show the intrinsic passivity of the compliance controller, we begin with the following lemma.
Lemma 1
Consider an arbitrary system with passive I/O pair . Premultiplication of the input by a matrix and postmultiplication of the output by preserve the passivity. Namely, the new I/O pair is passive, where and .
By assumption, . The passivity of the new I/O pair is trivial because .
The following lemma shows that can be obtained by a coordinate transformation of .
Lemma 2
can be expressed as
(25) 
Using
(26)  
(27) 
the following relation holds:
(28) 
Furthermore, noting that
(29) 
Finally, the following theorem shows the intrinsic passivity of compliance controller (24) using Lemma 1 and 2.
Theorem 1
The compliance controller (24) satisfies passivity of the I/O pair .
We begin with the fact that setpoint velocity PI control (which is equivalent to the setpoint position PD control) is passive. Using Lemma 2, the block diagram from to can be described by Fig. 4. Hence the I/O pair is passive by Lemma 1.
Please note that, to preserve passivity, has to be postmultiplied to , and this leads to as a control variable. The following section presents stability and passivity analysis of the controlled UAVM system.
IiiC Stability of the controlled UAVM system
The following theorem states stability during the freeflight.
Theorem 2 (Asymptotic stability for free flight)
See Appendix.
Remark 1
In [21], stability is analyzed based on (i) perfect regulation and (ii) perfect tracking. To achieved these, feedback linearization was used in [21]. However, feedback linearizing action includes coupling between UAV and manipulator, and consequently, the resulting control law becomes highly model dependent. In principle, modelbased controllers will probably outperform the modelfree ones, but at the cost of increased implementation complexity. In contrast, (13)(14) and (24) do not include feedback linearizing action. In the stability proof, this paper uses twotime scale analysis. This is a reasonable choice because the unactuated ,translational dynamics is slower than the actuated dynamics because the former is indirectly controlled by the rotational dynamics of UAV.
Summarizing the discussion in Section IIIB, the closedloop dynamics can be described by Fig. 3. Using skewsymmetric property of , the controlled UAVM (blue dashed box in Fig. 3) can be represented as feedback interconnections of passive subsystems. The following theorem states the stable environmental interaction.
Theorem 3 (Stable environmental interaction)
Assume that the environment can be modeled as a massspringdamper system so that Assumption 2 is valid. The overall closedloop dynamics including the environment is asymptotically stable to a certain equilibrium point.
The I/O pair of the controlled UAVM is output strictly passive due to Assumption 3. By small extension of Lemma 1, it can be easily shown that pre/postmultiplication of the I/O by preserves the output strict passivity. Therefore, the I/O pair is output strictly passive. Noting that Theorem 2 implies zerostate observability of the controlled UAVM, and that the I/O pair of the environment is strictly passive, asymptotic stability of the entire closedloop shown in Fig. 3 can be concluded by applying Theorem 6.3 in [26].
This theorem requires to model the environment using spring and damper elements. Because this argument must hold even with negligibly small damping value, the environment will be modeled as a pure spring in the simulation validations as an extreme case.
Iv Simulation Validation
To validate the proposed approach, the UAVM in Fig. 1 was simulated using rigid body dynamics (2). The mass and inertia of UAV were and . Due to Assumption 1, the 7th joint of the DLR light weight robot (LWR) manipulator was fixed. For every simulation, initial manipulator position was (recall that is the position/orientation of the endeffector from the UAV fuselage). For compliance controller, the stiffness and damping gains of translation were , , and those of rotation were , , respectively.
Iva Asymptotic stability during free flight
To validate stability during free flight, the following task was performed.

The desired UAV position: .

The desired UAV orientation: .

The desired endeffector position: At , .
Here, means the step command.
The results are shown in Fig. 5. As the UAV position error occurred at the beginning (because of the step command), the reference acceleration defined in (15) was excited. Consequently, UAV roll and pitch angles were perturbed to reduce the UAV position error, and asymptotic stability could be achieved. Notice that the manipulator was also perturbed due to the UAV motion, because of the dynamic coupling between UAV and manipulator. Recall that the controller (13)(14) is (almost) modelfree and does not cancel out the dynamic coupling. Conversely, at , UAV dynamics was perturbed due to the manipulator’s motion (see the magnified view in the first two rows), and was stabilized in a short instant.
IvB Stable interaction with passive environment
To validate stable environmental interaction, the following two tasks were performed.
IvB1 Task 1
In this task, physical wall located at in global direction was simulated using spring, as shown in Fig. 6a. Desired positions were set as follows.

The UAV position: (hovering).

The UAV orientation: .

The endeffector position: .
The simulation results are shown in Fig. 7. After contact occurred (around ), the external force pushed the UAV away, and consequently, external force became zero because the contact was lost. Due to the interaction, the UAV position was disturbed by about . The contact occurred again as the UAV recovered its desired position; see the third and fourth rows of Fig. 7. During the contact, UAVM converged to a certain equilibrium point at which it can balance the external force; converged to a certain value (recall Fig. 2b). Notice that the UAV converged to this equilibrium point automatically by virtue of passivity, without calculating the orientation that balances the interaction force. Note also that oscillation occurred on the endeffector (and hence on the interaction force as well) because the environment was modeled as a pure spring, but the was dissipated eventually. In this simulation, PC for UAV dynamics was not activated because the and were always positive.
IvB2 Task 2
In this task, we consider the forces acting on the manipulator along every direction. Therefore, the environment is modeled as springs in every translational direction with stiffness of , as shown in Fig. 6b. Note that this task considers only environmental interaction, whereas the previous task includes both freeflight phase and interaction phase. Desired positions were set as follows.

The UAV position: (hovering).

The UAV orientation: .

The endeffector position: .
The simulation results are shown in Fig. 8. Similar to task 1, the UAV orientation converged to a certain equilibrium point that balances the interaction force, while the UAV position converged to the desired position. However, in this task, mainly because of the integral action, the passivity condition was broken in body linear direction, and the PC was activated. By virtue of the PC, the energy was maintained to be positive (fourth row of Fig. 8). In addition, the energy of I/O port was observed to validate the passivity of the compliance controller. Passivity of this I/O port was maintained during the interaction because the energy was always less than the initially stored energy (see the fifth row of Fig. 8; 12.125 is the initially stored energy^{5}^{5}5 Initial spring displacement is and the stiffness was .). Therefore, we can conclude that the controlled UAVM was passive because every subblock in Fig. 3 was passive. As a result, stable interaction with the environment could be achieved.
V Conclusion
This paper presents a passive compliance control for UAVM to ensure stable interaction with passive environments. The key finding of this study is that the compliance controller satisfies passivity if the position of the endeffector is represented in the UAV fuselage frame. In addition to the passive manipulator controller, PO/PC technique is applied to ensure passivity of UAV controller. As a result, the controlled UAVM can interact with the passive environment stably. Simulation studies validate stability of free flight and stable environmental interaction.
Appendix
Since the control laws (13)(14) and (24) are the extension of our previous work [21], this section presents only the sketch of the proof. To begin with, let us express a new coordinates which has instead of .
(30) 
Recall introduced in (29). Using this coordinates, the UAVM dynamics is
(31) 
with properly defined and .
In the following, we apply twotime scale (also known as singular perturbation) analysis in which the linear , directional dynamics becomes slow dynamics and rotational dynamics become fast dynamics. This analysis is reasonable because the linear , directional motions are the consequences of orientation of the UAV.
For analysis purpose, let , , , and . Also, let us define the fast time scale by
(32) 
and the derivative of with respect to is defined as
(33) 
Note that and are the frozen variables in the fast time scale, because first two rows of (31) which represent linear  and directional dynamics can be expressed as
(34) 
Then, as , (31) can be written as
(35) 
Here, denotes the last block matrix (the subscript ‘’ stands for reduced). Therefore, it is trivial that the states converge to , , and in the fast time scale .
References
 [1] M. Orsag, C. Korpela, and P. Oh, “Modeling and control of mmuav: Mobile manipulating unmanned aerial vehicle,” Journal of Intelligent & Robotic Systems, pp. 1–14, 2013.
 [2] A. Suarez, G. Heredia, and A. Ollero, “Lightweight compliant arm with compliant finger for aerial manipulation and inspection,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2016, pp. 4449–4454.
 [3] M. Ryll, G. Muscio, F. Pierri, E. Cataldi, G. Antonelli, F. Caccavale, and A. Franchi, “6d physical interaction with a fully actuated aerial robot,” in IEEE International Conference on Robotics and Automation (ICRA), 2017.
 [4] M. J. Kim, J. Lin, K. Kondak, D. Lee, and C. Ott, “Oscillation damping control of pendulumlike manipulation platform using moving masses,” in Symposium on Robot Control (SYROCO), 2018.
 [5] M. Karrer, M. Kamel, R. Siegwart, and M. Chli, “Realtime dense surface reconstruction for aerial manipulation,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2016, pp. 1601–1608.
 [6] R. Rossi, A. SantamariaNavarro, J. AndradeCetto, and P. Rocco, “Trajectory generation for unmanned aerial manipulators through quadratic programming,” IEEE Robotics and Automation Letters, vol. 2, no. 2, pp. 389–396, 2017.
 [7] A. Pumarola, A. Vakhitov, A. Agudo, A. Sanfeliu, and F. MorenoNoguer, “Plslam: Realtime monocular visual slam with points and lines,” in IEEE International Conference on Robotics and Automation (ICRA), 2017.
 [8] H. Yang and D. Lee, “Dynamics and control of quadrotor with robotic manipulator,” in IEEE International Conference on Robotics and Automation (ICRA), 2014, pp. 5544–5549.
 [9] G. Garofalo, B. Henze, J. Englsberger, and C. Ott, “On the inertially decoupled structure of the floating base robot dynamics,” IFACPapersOnLine, vol. 48, no. 1, pp. 322–327, 2015.
 [10] S. Kim, S. Choi, and H. J. Kim, “Aerial manipulation using a quadrotor with a two dof robotic arm,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2013, pp. 4990–4995.
 [11] F. Huber, K. Kondak, K. Krieger, D. Sommer, M. Schwarzbach, M. Laiacker, I. Kossyk, S. Parusel, S. Haddadin, and A. AlbuSchäffer, “First analysis and experiments in aerial manipulation using fully actuated redundant robot arm,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2013, pp. 3452–3457.
 [12] K. Kondak, F. Huber, M. Schwarzbach, M. Laiacker, D. Sommer, M. Bejar, and A. Ollero, “Aerial manipulation robot composed of an autonomous helicopter and a 7 degrees of freedom industrial manipulator,” in IEEE International Conference on Robotics and Automation (ICRA), 2014, pp. 2107–2112.
 [13] M. Tognon, B. Yüksel, G. Buondonno, and A. Franchi, “Dynamic decentralized control for protocentric aerial manipulators,” in IEEE International Conference onRobotics and Automation (ICRA), 2017, pp. 6375–6380.
 [14] G. Garofalo, F. Beck, and C. Ott, “Taskspace tracking control for underactuated aerial manipulators,” in European Control Conference (ECC), 2018.
 [15] D. Mellinger, Q. Lindsey, M. Shomin, and V. Kumar, “Design, modeling, estimation and control for aerial grasping and manipulation,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2011, pp. 2668–2673.
 [16] J. Fink, N. Michael, S. Kim, and V. Kumar, “Planning and control for cooperative manipulation and transportation with aerial robots,” The International Journal of Robotics Research, vol. 30, no. 3, pp. 324–334, 2011.
 [17] A. JimenezCano, J. Martin, G. Heredia, A. Ollero, and R. Cano, “Control of an aerial robot with multilink arm for assembly tasks,” in IEEE International Conference on Robotics and Automation (ICRA), 2013, pp. 4916–4921.
 [18] M. Orsag, C. Korpela, M. Pekala, and P. Oh, “Stability control in aerial manipulation,” in American Control Conference (ACC), 2013, pp. 5581–5586.
 [19] V. Lippiello and F. Ruggiero, “Cartesian impedance control of a uav with a robotic arm,” IFAC Proceedings Volumes, vol. 45, no. 22, pp. 704–709, 2012.
 [20] ——, “Exploiting redundancy in cartesian impedance control of uavs equipped with a robotic arm,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2012, pp. 3768–3773.
 [21] M. J. Kim, K. Kondak, and C. Ott, “A stabilizing controller for regulation of uav with manipulator,” IEEE Robotics and Automation Letters, vol. PP, no. 99, pp. 1–1, 2018.
 [22] S. Rajappa, M. Ryll, H. H. Bülthoff, and A. Franchi, “Modeling, control and design optimization for a fullyactuated hexarotor aerial vehicle with tilted propellers,” in IEEE International Conference on Robotics and Automation (ICRA), 2015, pp. 4006–4013.
 [23] M. Tognon and A. Franchi, “Omnidirectional aerial vehicles with unidirectional thrusters: Theory, optimal design, and control,” IEEE Robotics and Automation Letters, vol. PP, no. 99, pp. 1–1, 2018.
 [24] B. Hannaford and J.H. Ryu, “Timedomain passivity control of haptic interfaces,” Robotics and Automation, IEEE Transactions on, vol. 18, no. 1, pp. 1–10, 2002.
 [25] J.H. Ryu, D.S. Kwon, and B. Hannaford, “Stable teleoperation with timedomain passivity control,” Robotics and Automation, IEEE Transactions on, vol. 20, no. 2, pp. 365–373, 2004.
 [26] H. K. Khalil and J. Grizzle, Nonlinear systems. Prentice hall Upper Saddle River, 2002, vol. 3.