IMPULSIVE CONTROL OF UNDERACTUATED MECHANICAL SYSTEMS By Sayyed Rouhollah Jafari Tafti A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Mechanical Engineering 2012 ABSTRACT IMPULSIVE CONTROL OF UNDERACTUATED MECHANICAL SYSTEMS By Sayyed Rouhollah Jafari Tafti Although there has been a significant amount of research in designing and analyzing impulsive control systems, there are very few applications of impulsive control in mechanical systems. In the first part of this dissertation, we investigate new control strategies for underactuated mechanical systems based on impulsive inputs. The control problem of underactuated systems is more challenging since such systems have fewer actuators than the number of their degrees of freedom. We first address the important concern related to application of impulsive control in mechanical systems, namely, implementation of impulse-like control inputs using standard hardware. This is done through experimental verification of an impulsive control algorithm for swing-up control of the Pendubot; the control algorithm was developed earlier in our research group. Showing the effectiveness of the impulsive control algorithm in experiments, we develop impulsive control algorithms for swing-up control of the Acrobot; the impulsive control algorithms have distinct advantages over existing algorithms in terms of the time required for swing-up and maximum control torque used by the continuous controller. Impulsive inputs cause jumps in velocity states of mechanical systems, and consequently, produce jumps in Lyapunov function candidates used in control design. This attribute is used to enlarge the region of attraction of equilibria using impulsive inputs at discrete instants of time. Several case studies of underactuated mechanical systems have been presented to demonstrate this benefit of using impulsive control. Another advantage of using impulsive inputs is that such inputs can significantly alter the dynamics of the system in a very short period of time. This property is used to design a safe fall algorithm for humanoid robots undergoing a fall, i.e., after the continuous controller has failed to keep the system trajectories confined to a fixed region around the equilibrium. The algorithm uses impulsive inputs to change the fall direction of the robot to minimize the damage to people and objects in the vicinity, as well as to its own self. In many instances, external disturbances or impact from interaction with the environment can have an adverse effect on system performance. In the second part of this dissertation, we develop control algorithms to mitigate these effects in underactuated biped robots. We first develop a disturbance rejection algorithm for the synthetic-wheel biped to mitigate the effects of external disturbances. A continuous controller is then designed for the synthetic-wheel biped to generate an impact-free walking gait. This gait consumes zero energy in the ideal case and the necessary conditions to achieve this gait are shown to be general and applicable to a range of bipedal robotic systems. An underactuated mechanical system can be non-minimum-phase if its zero dynamics is unstable. We finally investigate the output-tracking problem for linear non-minimum-phase systems. Intermittent output tracking is achieved under the condition of finite preview of the reference trajectory; the control algorithm uses switched inputs, which can be approximated by impulsive inputs in the limit. This dissertation is lovingly dedicated to my family, especially: to my wife, Hamideh, without her constant love, endless patience and deep understanding I could not accomplish this to my parents, Mahmoud and Sedigheh, whose support and encouragement have sustained me throughout my life. iv ACKNOWLEDGMENTS I would like to express my sincere gratitude to Dr. Ranjan Mukherjee for all his support, expertise and dedication throughout my research at the Dynamic Systems and Control Laboratory at Michigan State University. His passion in finding solutions to the challenging problems, his willingness in pursuing new research ideas and his attitude towards his students will always lead me in achieving my future career goals. I would like to thank Drs. Hassan Khalil, Steve Shaw and Jongeun Choi for the willingness to be on my Ph.D. committee in spite of their busy schedule and for their insightful suggestions whcih have enhanced this dissertation. I am also grateful to my colleagues in the Dynamic Systems and Control laboratory from whom I have benefited over the years; Louis Flynn, Frank Mathis, Dr. Aren Hellum, Dr. Assaad Awad Abbass, Dr. Nandagopal Methil, Paul Strefling, David Crouse, Mahdi Jadaliha, Dr. Yunfei Xu, Junho Lee. I would like to sincerely acknowledge all financial supports that I have received from the Department of Mechanical Engineering, the College of Engineering, the Graduate School at Michigan State University and the National Science Foundation. I am deeply indebted to my family and above all, my wife, Hamideh, for their full support, love and encouragement. Without their help, I would not be able to accomplish my Ph.D. career. Last but not the least, I wish to thank all my friends who have supported me in many different ways during my stay at Michigan State University. Their friendship has made it much easier to overcome the difficulties of being far from the family and the home country. v TABLE OF CONTENTS List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Chapter 1 Introduction . . . . . . . . . . . . . . . . 1.1 Motivation and Objectives . . . . . . . . . . . . . 1.2 Swing-up Control of Acrobot . . . . . . . . . . . . 1.3 Stabilizing Control for Underactuated Systems . . 1.4 Safe Fall Control for Humanoid Robots . . . . . . 1.5 Disturbance Rejection for Bipeds . . . . . . . . . 1.6 Optimizing the Energy Consumption for Bipeds . 1.7 Output Tracking Control for Non-Minimum-Phase 1.8 Scope of Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 4 5 7 9 11 15 17 Chapter 2 Swing-up Control of Pendubot: Experimental Validation 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Background: Swing-up Algorithm for Pendubot . . . . . . . . . . . . . 2.3 Experimental Setup and Results . . . . . . . . . . . . . . . . . . . . . . 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 19 20 21 23 Chapter 3 Swing-up Control of Acrobot . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Acrobot Dynamics and Impulsive Effects . . . . . . . . . 3.2.1 Equations of Motion . . . . . . . . . . . . . . . . 3.2.2 Holding Torque . . . . . . . . . . . . . . . . . . . 3.2.3 Impulsive Torque for Sudden Change in Velocity . 3.2.4 Change in Velocity due to Impulsive Torque . . . 3.2.5 Change in Energy due to Impulsive Torque . . . . 3.3 Swing-Up Control of the Acrobot . . . . . . . . . . . . . 3.3.1 Rest-to-Rest Maneuvers of the Second Link . . . 3.3.2 Acrobot Swing-up Using Rest-to-Rest Maneuvers 3.3.2.1 Swing-up Algorithm . . . . . . . . . . . 3.3.2.2 Numerical Simulation . . . . . . . . . . 3.3.3 Increasing Energy of the Acrobot Using Impulsive 3.3.3.1 Swing-up Algorithm . . . . . . . . . . . 3.3.3.2 Numerical Simulation . . . . . . . . . . 3.4 Modification of an Existing Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 25 26 26 27 28 29 31 32 32 34 34 37 38 38 41 43 vi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 3.4.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 45 45 47 51 52 52 54 55 Chapter 4 Enlarging Region of Attraction for Underactuated Systems 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Dynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Impulsive Inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Effect of Impulsive Inputs on Velocity . . . . . . . . . . . . . . . . 4.3 Stabilization Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Case Studies: Linear Stabilizing Controllers . . . . . . . . . . . . . . . . 4.4.1 The Pendubot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1.1 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1.2 Continuous Controller . . . . . . . . . . . . . . . . . . . 4.4.1.3 Numerical Simulations . . . . . . . . . . . . . . . . . . . 4.4.2 The Acrobot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2.1 Dynamics and Control . . . . . . . . . . . . . . . . . . . 4.4.2.2 Numerical Simulations . . . . . . . . . . . . . . . . . . . 4.4.3 Rolling Acrobot . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3.1 Dynamics and Control . . . . . . . . . . . . . . . . . . . 4.4.3.2 Numerical Simulations . . . . . . . . . . . . . . . . . . . 4.5 Extension for Nonlinear Controllers . . . . . . . . . . . . . . . . . . . . . 4.5.1 Inverted Pendulum on an Inclined Plane . . . . . . . . . . . . . . 4.5.1.1 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.1.2 Stabilizing Control Design . . . . . . . . . . . . . . . . . 4.5.1.3 Numerical Simulations . . . . . . . . . . . . . . . . . . . 4.5.2 The Ball and Beam System . . . . . . . . . . . . . . . . . . . . . 4.5.2.1 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.2.2 Stabilizing Control Design . . . . . . . . . . . . . . . . . 4.5.2.3 Numerical Simulations . . . . . . . . . . . . . . . . . . . 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 57 58 58 58 59 61 64 64 64 65 66 69 69 70 72 72 74 78 80 80 81 82 85 85 86 88 90 Chapter 5 Safe Fall Control for Humanoid Robots 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . 5.2 Dynamic Model . . . . . . . . . . . . . . . . . . . . 5.2.1 Equations of Motion . . . . . . . . . . . . . 5.2.2 Lean Line . . . . . . . . . . . . . . . . . . . 5.3 Impulsive Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 91 92 92 94 95 3.5 Existing Control Method . . . . . . . . . Modified Method Using Impulsive Inputs 3.4.2.1 Impulsive System . . . . . . . . 3.4.2.2 Stability of Impulsive System . 3.4.2.3 Numerical Simulation . . . . . 3.4.3 Alternate Continuous Control Design . . 3.4.3.1 Control Design . . . . . . . . . 3.4.3.2 Numerical Simulation . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . vii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 5.5 5.6 Safe Fall Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Chapter 6 Disturbance Rejection for the Synthetic-Wheel Biped 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Torso Effect on Biped Dynamics . . . . . . . . . . . . . . . . 6.2.3 Interchange of Stance and Swing Legs . . . . . . . . . . . . . 6.2.4 Control Design . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Impulsive Torques and Effects . . . . . . . . . . . . . . . . . . . . . 6.3.1 Braking Torque for the Swing Leg . . . . . . . . . . . . . . 6.3.2 Braking Torque for the Torso . . . . . . . . . . . . . . . . . 6.3.3 Impulsive Effect on Velocity . . . . . . . . . . . . . . . . . . 6.4 Disturbance Rejection Algorithm . . . . . . . . . . . . . . . . . . . 6.5 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 106 107 107 108 109 110 113 113 114 115 116 119 123 Chapter 7 Energy-Conserving Gaits for Point-Foot Passive-Ankle Bipeds 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Energy-Conserving Gaits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.1 Assumptions and Definitions . . . . . . . . . . . . . . . . . . . . . . . 7.2.2 Sufficient Conditions for Energy-Conserving Gait . . . . . . . . . . . 7.2.3 Characterization of Trajectories of the Actuated Joints . . . . . . . . 7.3 Three-DOF Bipeds: Two Case Studies . . . . . . . . . . . . . . . . . . . . . 7.3.1 Synthetic-Wheel Biped (SWB) . . . . . . . . . . . . . . . . . . . . . 7.3.1.1 System Description . . . . . . . . . . . . . . . . . . . . . . . 7.3.1.2 Sufficient Conditions for Energy-Conserving Gait . . . . . . 7.3.1.3 Determination of Actuated Joint Trajectories . . . . . . . . 7.3.1.4 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . 7.3.2 Point-Foot Three-Link Biped . . . . . . . . . . . . . . . . . . . . . . 7.4 Five-DOF Point-Foot Biped . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4.1 System Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4.2 Sufficient Conditions for Energy-Conserving Gait . . . . . . . . . . . 7.4.3 Determination of Actuated Joint Trajectories . . . . . . . . . . . . . 7.4.4 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Two-Step Periodic Gaits: Five-DOF Bipeds . . . . . . . . . . . . . . . . . . 7.5.1 Generalization to Two-Step Periodic Gait . . . . . . . . . . . . . . . 7.5.2 Characterization of Trajectories . . . . . . . . . . . . . . . . . . . . . 7.5.2.1 Trajectories for the First Step . . . . . . . . . . . . . . . . . 7.5.2.2 Trajectories for the Second Step . . . . . . . . . . . . . . . . 7.5.3 Five-DOF Point-Foot Biped: Similar Leg Links . . . . . . . . . . . . 7.5.4 Five-DOF Point-Foot Biped: Dissimilar Leg Links . . . . . . . . . . . 7.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 124 125 125 127 128 131 131 131 132 133 135 137 139 139 140 141 144 146 146 147 147 148 150 152 155 viii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 8 8.1 8.2 8.3 8.4 8.5 8.6 Intermittent Output Tracking For Linear Non-Minimum-Phase Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discrete Equivalent Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . Zero Shift for Linear NMP Systems . . . . . . . . . . . . . . . . . . . . . . . 8.3.1 Stable NMP Systems with Distinct Poles - General Case . . . . . . . 8.3.2 Stable NMP Systems with Distinct Poles - Special Case . . . . . . . . 8.3.3 Stable NMP Systems with Repeated Poles . . . . . . . . . . . . . . . Intermittent Output Tracking Control . . . . . . . . . . . . . . . . . . . . . . 8.4.1 Tracking Control for MP Systems . . . . . . . . . . . . . . . . . . . . 8.4.2 Intermittent Output Tracking for NMP Systems . . . . . . . . . . . . 8.4.3 Main Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 9 157 157 158 162 162 165 167 168 169 171 176 177 181 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . 183 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 Appendix A Procedure for Verifying the Condition in Theorem 3.4.1 . . . . 187 Appendix B Equations of Motion: Synthetic-Wheel Biped . . . . . . . . . . 189 Appendix C Equations of Motion: Five-DOF Biped . . . . . . . . . . . . . . . 191 Bibliography . . . . . . . . . . . . . . . ix . . . . . . . . . . . . . . . . . . . 194 LIST OF TABLES Table 1.1 Mechanical Cost of Transport of Bipeds . . . . . . . . . . . . . . . . 13 Table 5.1 Degrees of freedom for the humanoid robot in Figure 5.1 93 Table 5.2 Kinematic and dynamic parameters for the six-dof robot with cylindrical links. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Table 6.1 Parameters of the Synthetic-Wheel Biped [1] . . . . . . . . . . . . . 119 Table 7.1 Kinematic and dynamic parameters of Synthetic Wheel Biped Table 7.2 Kinematic and dynamic parameters of Five-DOF Point-Foot biped with similar leg links . . . . . . . . . . . . . . . . . . . . . . . . . . 144 Table 7.3 Kinematic and dynamic parameters of Five-DOF Point-Foot biped with dissimilar leg links . . . . . . . . . . . . . . . . . . . . . . . . . 152 Table 8.1 Linear NMP systems with repeated poles and their corresponding MP DE systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 x . . . . . . . . . 135 LIST OF FIGURES Figure 2.1 The pendubot is shown in an arbitrary configuration. . . . . . . . . 20 Figure 2.2 Experimental validation of the impulsive control algorithm [2] for pendubot swing-up. . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Figure 3.1 The Acrobot in an arbitrary configuration. . . . . . . . . . . . . . . 27 Figure 3.2 ˙ Plot showing θ1 -θ2 space where Π is positive and negative for θ1 = ˙ θ2 = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Simulation results for swing-up of the Acrobot using rest-to-rest maneuvers of the second link. . . . . . . . . . . . . . . . . . . . . . . . 38 Configuration of the Acrobot at the time when the impulsive torque is applied to increase the energy of the system. . . . . . . . . . . . . 41 Simulation results for swing-up of the Acrobot using impulsive inputs to increase the energy. . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Simulation results for swing-up of the Acrobot: Comparison of the modified impulsive controller (solid line) and the continuous controller (dashed line) proposed by Xin and Kaneda [3]. . . . . . . . . 52 Simulation results for swing-up of the Acrobot: Comparison of the controller in Section 3.4.3 (solid line) with the continuous controller (dashed line) proposed by Xin and Kaneda [3]. . . . . . . . . . . . . 55 ˆ ˙ Plots of the curves V = V and V = 0 and the impulse line correspond to fixed values of q = (ˆ1 , q2 ) for a two-dof system. The prior-impulse q ˆ configuration at A : (ˆ, q − ) is jumped to a configuration at B : (ˆ, q + ) q ˙ q ˙ ˙ < 0 and V < V . . . . . . . . . . . . . . . . . . . . . . . for which V 63 The Pendubot/Acrobot is shown in an arbitrary configuration. For the Pendubot we have τ2 = 0 and for the Acrobot we have τ1 = 0. . 65 Figure 3.3 Figure 3.4 Figure 3.5 Figure 3.6 Figure 3.7 Figure 4.1 Figure 4.2 xi Figure 4.3 Simulation results showing joint angles and joint velocities of the Pendubot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - First simulation 67 Figure 4.4 ˙ ˙ ˙ The curves V = V0 , V = 0, and the impulse line are shown in the θ1 -θ2 plane at the initial configuration of the Pendubot for first simulation, where θ10 = 1.484 and θ20 = 0.087. The impulsive input at t = 0 in Figure 4.3 results in a jump in the configuration from A : (q0 , q0 ) to ˙− ˙ B : (q0 , q0 ), where V (q0 , q0 ) < 0 and V (q0 , q0 ) < V0 . . . . . . . . . ˙+ ˙+ ˙+ 68 Figure 4.5 Simulation results showing joint angles and joint velocities of the Pendubot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - Second simulation 69 Figure 4.6 Simulation results showing joint angles and joint velocities of the Acrobot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - First Simulation . . . 71 Simulation results showing joint angles and joint velocities of the Acrobot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - Second Simulation . . 72 A schematic of the rolling acrobot. The link angles θ1 and θ2 are measured counter-clockwise. . . . . . . . . . . . . . . . . . . . . . . 73 Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with unbounded leg angle originally outside RA - First simulation . . . . . . . . . . . . . . . . . . . . . . 75 Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with unbounded leg angle originally outside RA - Second simulation . . . . . . . . . . . . . . . . . . . . 76 Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with bounded leg angle originally outside RA - First simulation . . . . . . . . . . . . . . . . . . . . . . 77 Figure 4.7 Figure 4.8 Figure 4.9 Figure 4.10 Figure 4.11 xii Figure 4.12 Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with bounded leg angle originally outside RA - Second simulation . . . . . . . . . . . . . . . . . . . . 78 Figure 4.13 Cart-pendulum on a slope . . . . . . . . . . . . . . . . . . . . . . . 81 Figure 4.14 Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the inverted pendulum on a slope originally outside the RA - First simulation . . . . . . . . . . . . . . 83 Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the inverted pendulum on a slope originally outside the RA - Second simulation . . . . . . . . . . . . . 84 Figure 4.16 Ball and beam system . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Figure 4.17 Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the ball and beam system originally outside the RA - First simulation . . . . . . . . . . . . . . . . . . . 88 Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the ball and beam system originally outside the RA - Second simulation . . . . . . . . . . . . . . . . . . 89 Figure 5.1 Seven-dof three-dimensional humanoid robot . . . . . . . . . . . . . 93 Figure 5.2 Orientation of the lean line in the inertial frame . . . . . . . . . . . 94 Figure 5.3 A schematic of the safe fall algorithm using two sets of impulsive inputs. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Figure 4.15 Figure 4.18 Figure 5.4 Simulation results showing all the joint angles and angular velocities for changing the fall direction of the robot . . . . . . . . . . . . . . 104 xiii Figure 5.5 Simulation results showing the lean line angles for changing the fall direction of the robot . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Figure 5.6 Simulation results showing the impulsive inputs for changing the fall direction of the robot . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Figure 6.1 A schematic of the Synthetic-Wheel biped, reproduced from [1]. . . . 107 Figure 6.2 Biped configuration at the time of interchange of stance and swing legs.109 Figure 6.3 Simulation results showing angular positions (rad), angular velocities (rad/s), and control inputs (N.m) of the SWB in rejecting the impulsive disturbances applied at two instants of time. The dashed lines ˙ correspond to the plots of ψ and ψ. . . . . . . . . . . . . . . . . . . 121 Figure 6.4 Simulation results showing angular positions (rad), angular velocities (rad/s), and control inputs (N.m) of the SWB in rejecting the impulsive disturbances applied at five instants of time. The dashed lines ˙ correspond to the plots of ψ and ψ. . . . . . . . . . . . . . . . . . . 122 Figure 7.1 A schematic of an n-link point-foot planar biped. . . . . . . . . . . . 126 Figure 7.2 A schematic of the Synthetic-Wheel biped, reproduced from [1]. . . . 132 Figure 7.3 Phase portraits of the Synthetic-Wheel Biped for different amplitudes of torso oscillation A . . . . . . . . . . . . . . . . . . . . . . . . . . 136 Figure 7.4 Torque plots for an energy-conserving gait for the Synthetic-Wheel Biped (SWB). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 Figure 7.5 Generalized torque plots for an energy-conserving gait for the SyntheticWheel Biped (SWB). . . . . . . . . . . . . . . . . . . . . . . . . . . 137 Figure 7.6 A schematic of the Point-Foot Three-Link. . . . . . . . . . . . . . . 138 Figure 7.7 A schematic of the five-dof point-foot biped. . . . . . . . . . . . . . 140 Figure 7.8 Torque plots for an energy-conserving gait for the Five-DOF PointFoot Biped with similar leg links. . . . . . . . . . . . . . . . . . . . 144 Figure 7.9 Generalized torque plot for an energy-conserving gait for the FiveDOF Point-Foot Biped with similar leg links. . . . . . . . . . . . . . 145 xiv Figure 7.10 Configurations of the five-dof point-foot biped with similar leg links at different instants of time over one step. . . . . . . . . . . . . . . . 146 Figure 7.11 Torque plots for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with similar leg links. . . . . . . 151 Figure 7.12 Generalized torque plot for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with similar leg links. 151 Figure 7.13 Configurations of the five-dof point-foot biped with similar leg links at different instants of time for a two-step periodic gait. The joint configurations at the beginning of the first step and end of the second step are the same but they are different from the joint configuration at the end of the first step. . . . . . . . . . . . . . . . . . . . . . . . 152 Figure 7.14 Torque plots for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with dissimilar leg links. . . . . 153 Figure 7.15 Generalized torque plot for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with dissimilar leg links.153 Figure 7.16 Configurations of the five-dof point-foot biped with dissimilar leg links at different instants of time for a two-step periodic gait. The joint configurations at the beginning of the first step and end of the second step are the same but they are different from the joint configuration at the end of the first step. . . . . . . . . . . . . . . . . . . . . . . . 154 Figure 8.1 The periodic switched input for the system in Eq.(8.1) (in solid line) and the constant input for DE system in Eq.(8.3) (in dashed line) . 159 Figure 8.2 Simulation results for example 8.2.1 showing state variables and inputs for switched-input system (solid line) and discrete equivalent system (dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . 161 Figure 8.3 Simulation results for the first case showing the system output (solid line) and reference trajectory (dashed line) in the first plot, internal states in the second plot and control input (Solid line) and tracking control for MP system (dashed line) in the third plot . . . . . . . . 179 Figure 8.4 Simulation results for the second case showing the system output (solid line) and reference trajectory (dashed line) in the first plot, internal states in the second plot and control input (Solid line) and tracking control for MP system (dashed line) in the third plot . . . . 181 xv Chapter 1 Introduction 1.1 Motivation and Objectives Impulsive systems are dynamical systems which can have discontinuous jumps in their state variables. These jumps are caused by either impulsive inputs which are intentionally applied to control the system or external disturbances and impacts which are generally undesirable. The impulsive inputs are theoretically modeled as Dirac-delta functions with unbounded magnitude and infinitesimal time support. Since the state variables are not always continuous for impulsive systems, control theories developed for continuous systems need to be revised for application to impulsive systems. There has been a fair amount of theoretical research on impulsive control of dynamical systems and credit for some of the early works goes to Pavlidis [4], Gilbert and Harasty [5], Menaldi [6] and Lakshmikantham [7]. Researchers have studied the problems of stability, controllability and observability, optimality (see [8, 9, 10], for example). Impulsive control has been designed for dynamical systems with different control objectives [11, 12, 13]. However, impulsive control can only be applied to systems for which the state variables can be changed instantaneously without violating any physical laws. Population growth control systems, financial models or chaotic dynamical systems are some examples of such systems. In mechanical systems, the discontinuous jumps occur in the velocity states which are directly changed through the inputs of the system such as motor 1 torques. Therefore, it has to be clarified how to implement ideal delta-function inputs using standard actuators before applying impulsive control algorithms to mechanical systems. The implementability concern can be the main reason why there are few works in the literature designing impulsive control for mechanical systems [14, 15, 2, 16, 17]. Weibel et al [14] studied the implementation of impulsive inputs in controlling single-degree-of-freedom Hamiltonian mechanical systems. They derive the delta-function impulsive inputs to make desired jumps in the energy levels of the system. The impulsive inputs are approximated by a series of rectangular pulses and were experimentally verified on a pendulum-cart system. In this research, we first confirm that impulsive inputs can be implemented in mechanical systems with an experimental validation of swing-up control of the Pendubot [2] using impulsive inputs. Unlike the approach presented in [14], the impulsive inputs are obtained from Lagrangian dynamics and are approximated by large-gain continuous inputs applied during short intervals of time. After demonstrating the effectiveness of an impulsive control algorithm in experiments, we develop impulsive control algorithms for underactuated mechanical systems. The control problem of underactuated systems is more challenging since such systems have fewer actuators than the number of their degrees of freedom. An impulsive algorithm is developed for swing-up control of the Acrobot, a benchmark problem in underactuated systems, where impulsive inputs are used to make instantaneous changes in energy level of the acrobot. The impulsive control algorithm has distinct advantages over existing algorithms in terms of the time required for swing-up and maximum control torque used by the continuous controller. Impulsive inputs cause jumps in velocity states of mechanical systems, and consequently, produce jumps in Lyapunov function candidates used in control design. This attribute is used to enlarge the region of attraction of equilibria using impulsive inputs at discrete instants 2 of time. Several case studies of underactuated mechanical systems have been presented to demonstrate this benefit of using impulsive control. Another advantage of using impulsive inputs is that such inputs can significantly alter the dynamics of the system in a very short period of time. This property is used to design a safe fall algorithm for humanoid robots undergoing a fall, i.e., after the continuous controller has failed to keep the system trajectories confined to a fixed region around the equilibrium. The algorithm uses impulsive inputs to change the fall direction of the robot to minimize the damage to people and objects in the vicinity, as well as to its own self. In many instances, external disturbances or impact from interaction with the environment can have an adverse effect on system performance. In the second part of this dissertation, we develop control algorithms to mitigate these effects in underactuated biped robots. We first develop a disturbance rejection algorithm for the synthetic-wheel biped [1] to mitigate the effects of external disturbances. The control algorithm uses impulsive inputs to reject the effects of the external disturbances as well as cancel the undesirable effects of ground impacts during leg interchange. Continuous controller is then designed for the synthetic-wheel biped to generate an impact-free walking gait. This gait consumes zero energy in the ideal case and the necessary conditions to achieve this gait are shown to be general and applicable to a range of bipedal robotic systems. An underactuated mechanical system can be non-minimum-phase if its zero dynamics is unstable. We finally investigate the output-tracking problem for linear single-input-singleoutput non-minimum-phase systems. Intermittent output tracking is achieved under the condition of finite preview of the reference trajectory, i.e., the system output tracks a desired trajectory intermittently in some regular intervals of time rather than continuous tracking of the trajectory. The control algorithm uses switched inputs, which can be approximated 3 by impulsive inputs in the limit. In the remainder of this chapter, we give a brief description and review the literature for each control problem which will be studied in this dissertation. 1.2 Swing-up Control of Acrobot The Acrobot is a two-link robot in the vertical plane with an actuator at the elbow joint and a passive shoulder joint. The dynamic model of the Acrobot has four equilibria, all of which are stabilizable, and the standard control problem refers to swing-up from an arbitrary initial configuration to the equilibrium configuration with the highest potential energy. The Acrobot has identical kinematic structure as that of the Pendubot [2], but swing-up of the Acrobot is more challenging due to the adverse location of its actuated joint. Most of the results in the literature provide solutions that apply either to the Pendubot or the Acrobot. The paper by Boone [18], for example, provides a near-optimal solution to the Acrobot swing-up problem using bang-bang control. More interesting are the papers that can be applied to the Pendubot and the Acrobot alike. Included in this category are the papers by Spong [19, 20] based on feedback linearization, and Kolesnichenko and Shiriaev [21] based on partial stabilization using Lyapunov-like functions. The feedback linearization method proposed by Spong [19] designs reference trajectories for the actuated joint that render the zero dynamics unstable and thereby increases the energy of the system. For the Acrobot dynamic model, this method is highly sensitive to the control gains. Furthermore, the control design is approximated for implementation without a rigorous proof. The partial stabilization method proposed by Kolesnichenko and Shiriaev [21] has been implemented in the Acrobot dynamic model by Xin and Kaneda [3] 4 and Mahindrakar and Banavar [22, 23]. The method uses a Lyapunov-like function and results in an invariant set that contains points other than the trivial solution. For swing-up to the desired configuration, constraints have to be imposed on the control gains and certain initial conditions have to be avoided [3, 22]. 1.3 Stabilizing Control for Underactuated Systems Stabilization of underactuated dynamical systems has been widely studied and several control methods have been proposed. Linear regulation methods such as pole placement and LQR can be designed based on the linearized system but the Region of Attraction (RA ) of the desired equilibrium for such controllers is typically small. The main nonlinear control methods in the literature modify the structure of the original system to yield a closed-loop system for which the desired equilibrium is asymptotically stable. The Controlled Lagrangian method, Interconnection and Damping Assignment Passivity-Based Control (IDA-PBC), and the λ-Method are the three main approaches that have been developed. The Controlled Lagrangian method, proposed by Bloch et al. [24, 25], is an energy-based approach that can be applied to underactuated mechanical systems in Euler-Lagrange form. The kinetic and potential energies of the system are modified to find a desired Lagrangian, called Controlled Lagrangian. The control gains and the system input are then chosen such that the closed-loop equations derived from the Controlled Lagrangian match those of the controlled system. Full phase space stabilization is obtained through kinetic and potential energy shaping. The control procedure can be complicated for systems with several degrees of freedom (dof’s). Passivity-based control can be effectively used to stabilize some underactuated mechanical 5 systems through shaping of the potential energy [26]. However, for most underactuated systems, the total energy needs to be shaped. For such systems, Ortega et al. [27, 28] proposed IDA-PBC which considers the Hamiltonian form of the system. The desired energy function is found by assigning interconnection and damping matrices of the system. Similar to the Controlled Lagrangian method, where the closed-loop system is maintained in EulerLagrange form, the Hamiltonian form of the closed-loop system is preserved in IDA-PBC. The key advantage of IDA-PBC is that it provides freedom to change the interconnection matrix as well as the energy function; the Controlled Lagrangian method can therefore be viewed as a special case of IDA-PBC. The class of underactuated systems which can be stabilized by IDA-PBC is limited to those for which two partial differential equations can be solved. Acosta et al. [29] showed that these equations can be explicitly solved for mechanical systems with underactuation degree one. The λ-Method, proposed by Auckly et al. [30, 31], uses differential geometry to assign a special form for the closed-loop system. This method produces an infinite-dimensional family of control laws to stabilize the underactuated system. Using a transformation, the control laws are derived using matching conditions that are recast as linear first-order partial differential equations. The control laws have been used to stabilize linear systems and the nonlinear ball and beam system [32]. Among other stabilizing controllers, White et al. [33, 34, 35] proposed a direct Lyapunov function approach which results in three matching conditions comprised of ordinary and partial differential equations. The matching conditions can be numerically solved in real time; this makes it better suited for systems with several dof’s compared to the Lagrangian and Hamiltonian approaches. Olfati-Saber proposed a method to transform underactuated mechanical systems into a cascade normal form comprised of linear and nonlinear subsystems 6 [36, 37, 38, 39]. The transformation is chosen such that the control problem is reduced to control of the reduced-order nonlinear subsystem. The backstepping method is then used to design a stabilizing controller for the resulting cascade normal form. Sliding mode controllers have also been proposed to stabilize a class of underactuated systems under uncertainty conditions [40, 41, 42]. 1.4 Safe Fall Control for Humanoid Robots Safety is one of the main concerns behind the independent autonomous existence of humanoid robots in physically interactive human environments. Although the loss of balance and fall are rare in isolated or controlled environments, the situation would be different when the physical separation between robots and humans gradually disappear. A fall of a humanoid robot is a particularly worrisome event; a fall from an upright posture can cause damage to the robot, to delicate and expensive objects in the surrounding or can inflict injury to a nearby human being. Regardless of the substantial progress in humanoid robot balance control strategies, the possibility of fall remains real, even unavoidable. Yet, few comprehensive studies of humanoid fall encompassing fall avoidance, prediction, control and recovery exists in the literature. A humanoid fall may be caused due to external factors such as unexpected or excessive forces, unusual or unknown floor slipperiness, slope or profile of the ground, causing the robot to slip, trip or topple. In these cases the disturbances that threaten balance are larger than what the balance controller can handle. Fall can also result from internal factors such as a malfunction or a failure of actuator, power or communication system where the balance controller is partially or fully incapacitated. 7 A controller dealing with an accidental humanoid fall may have two primary, and distinctly different goals: a) self-damage reduction and b) reduction of damage to others. When a fall occurs in an open space, the self-damage reduction strategy can aim to mitigate the damaging effects of the ground impact. If, however, a robot falls in close proximity to others, its primary objective should be to prevent damage to these objects or minimize injury to persons. Recently, Wilken et al. have reported a third possible goal of a fall controller, that of a deliberate fall of a humanoid soccer goalie [43]. This is an example of a strategic fall. A falling humanoid can be easily appreciated as a formidable system to study and control. During fall, a humanoid behaves as a nonlinear underactuated multi-body system with highDOF, changing morphology and unpredictable contact configurations. Additionally, the robot is subjected to the relentless gravitational pull, which it cannot fully resist, and which tries to bring it down to a crash. Consequently, the time interval between the detection of a fall and the actual event of ground impact becomes a critical parameter for fall control. Time is at a premium during the occurrence of a fall; a single rigid body model of a full sized humanoid indicates that a fall from the vertical upright stationary configuration due to a mild push takes about 800-900 ms [44]. In many situations the time to fall can be significantly shorter, and there is no opportunity for elaborate planning or time-consuming control. Yet, through simulation and experiments we are able to demonstrate that meaningful modification to the default fall behavior can be achieved and damage to the environment can be avoided. Some earlier work on humanoid fall direction change has been reported in [45, 46, 47]. A number of recent papers reported on the damage minimization and prediction aspects of humanoid fall. In their exhaustive work, Fujiwara et al. [48, 49, 50, 51] proposed martial arts type motion for damage reduction, computed optimal falling motions using minimum impact 8 and angular momentum, and fabricated special hardware for fall damage study. Ogata et al. proposed two fall prediction methods based on abnormality detection and predicted Zero Moment Point (ZMP) [52, 53]. The robot improves fall prediction through experimental learning. Renner and Behnke [54] use model-based approach to detect external forces on the robot and Karssen and Wisse [55] use principal component analysis to predict fall. Following human movement based search procedure, Ruiz-del-Solar et al. implemented a low damage fall sequence for soccer robots [56]. In [57] and [58], fall prediction and control are treated together using Gaussian mixture models and Hidden Markov model. Ishida et al. employed servo loop gain shift to reduce shock due to fall [59]. Kanoi et al. worked on fall detection from walk patterns [60]. Fall damage minimization is obviously of natural interest in biomechanics [61]. 1.5 Disturbance Rejection for Bipeds A humanoid robot is eventually aimed to operate in real environments with all interactions and disturbances from humans and surrounding obstacles. Therefore, the control design for humanoids must be able to reject the external disturbances which can be applied during stance or walking phases of the robot. The complicated, nonlinear, unstable and hybrid dynamics of the bipeds make it a challenging problem to find effective balancing strategies for humanoid robots. Based on biomechanical studies [62, 63], Hofmann [64] suggested three basic push recovery strategies for bipeds: Ankle strategy(or CoP balancing) , hip strategy (or CMP balancing) and stepping. If the external disturbances are small, the ankle joint is used to adjust the Center of Pressure (CoP) of the robot within the foot support area to maintain the balance. 9 For greater disturbances, other joints such as hip joint and arms are used to create angular momentum about the Center of Mass (CoM) to counteract the external push. This strategy leads to an effective CoP outside the support polygon which is called Centroidal Moment Point (CMP). For large disturbances, the biped has to take steps to maintain balance. Based on these strategies, different controllers have been proposed to solve push recovery problem in humanoids. ZMP-based approaches have been widely used to design balance controllers for humanoid bipeds [65, 66, 67, 68, 69, 70]. Stephens [71] derived decision surfaces as functions of reference points of a humanoid to determine when each balancing strategy is needed and to predict whether or not the biped can recover from the push. He used simplified models for humanoids which result in inaccurate decision surfaces for actual robots. Yi et al. [72] proposed reinforcement learning to train a high-level controller which chooses the proper strategy for disturbance rejection. Their method is more accurate but it requires off-line procedures to train the controller. Linear Inverted Pendulum Model (LIPM)[73, 74] is used for most of balance controllers using ankle and hip strategies (e.g. in [75]). Since LIPM assumes zero angular momentum about CoM, other models such as Linear Inverted Pendulum model plus Flywheel (LIPF) [76] and Angular Momentum inducing inverted Pendulum Model (AMPM) [77] are proposed to incorporate the significant role of angular momentum generated by internal joints of the humanoid. Pratt, et al. [76] introduced the concept of “Capture point”, a point on the ground where the robot can step to in order to stop and found the analytical expressions for LIPF. Rebula, et al. [78] presented a learning method to compute the capture points for a general humanoid. The proposed method uses the estimated capture points predicted by LIPM as the starting points and finds the desired foot placement to recover from the push. Their method was validated on a three-dimensional humanoid robot with twelve actuated dof’s. Stephens and Atkeson 10 extended the humanoid balance strategies for compliant robots with force-controlled joints [79], [80]. Their standing balance controller uses a model-based approach called Dynamic Balance Force Control (DBFC) to calculate the joint torques based on desired COM motion [79]. Their stepping controller for larger disturbances was based on feedforward torque inputs obtained from DBFC to follow desired CoM trajectories generated by a linear model predictive control [80]. They validated their push recovery strategies through experiments on Sarcos, a force-controlled humanoid. Komura et al. [81] and Adiwanhono et al. [82] proposed balance controllers for walking bipeds which uses AMPM and LIPM, respectively. In [81], the amount of angular momentum needed to counteract the external disturbance in both sagital and frontal planes is calculated based on a new criteria called the “difference of inertia”. The push recovery algorithm in [82] employs a different walking gait after the push which is based on a ZMP-shifting method. 1.6 Optimizing the Energy Consumption for Bipeds With biped walking machines becoming quite common, robotics engineers are increasingly using the metric of Cost of Transport (COT) to compare their energy efficiencies. The COT is a dimensionless measure and is defined as [83] COT = energy used body weight × distance traveled It has been used to quantify the energetics of locomotion in animals and is now being used to compare energy efficiencies of walking machines and other mobile robots. For a biped robot, the COT can be computed based on the total energy used by the robot, or the energy 11 used by the actuators driving the locomotion mechanism. To distinguish the efficiency of the mechanical design from the efficiency of the complete system, researchers ([84], for example) have defined the mechanical cost of transport as Cmt = energy used by actuators body weight × distance traveled Although there is no ambiguity in the definition of Cmt , the “energy used by the actuators” has been calculated in different ways by different researchers. Hobbelen and Wisse [84] made the assumption that actuators expend energy irrespective of whether the mechanical energy is positive or negative. On the other hand, Srinivasan and Ruina [85] assumed that the energy used by the actuators is positive when the mechanical energy is positive but is equal to zero when the mechanical energy is negative. Both Hobbelen and Wisse [84] and Srinivasan and Ruina [85] have attempted to take into account the limitations of technology in regenerating energy from motion and storing it back in the energy storage unit. Their assumptions are certainly not consistent with one another and are likely to change further if actuators other than electric motors are used. To avoid an actuator-specific definition of Cmt , and in light of the fact that technological advancements (in hybrid cars, for example) are indeed making it possible to regenerate and store energy, the original definition of Cmt can be reworded as follows: Cmt = mechanical energy required by gait body weight × distance traveled where T mechanical energy required by gait = i 12 0 ˙ τi θi dt (1.1) ˙ and where τi and θi are the generalized force and generalized velocity of the i-th actuated joint, and T is the time period of each step. Many successful biped designs have been reported in the literature. The Cmt of some of these designs are compiled in Table 1.1 and compared with that of human walking. All of these values were derived from experimental data and are therefore based on the energy used by the actuators. Although different methods may have been used for calculation of the energy used by the actuators, the Cmt values of these designs span a wide range. At one end of the spectrum is Dynamite [83], a passive dynamic walker that can only walk down inclines for certain initial conditions; at the other end of the spectrum is Honda Asimo [92, 93, 94], an active biped that can walk, run, and climb stairs, and has many other functionalities apart from locomotion. Dynamite and Asimo are very different platforms and no meaningful comparisons can be made, but it should be emphasized that Dynamite relies on natural limb dynamics whereas Asimo uses significant amount of energy to overcome natural limb dynamics and track prescribed trajectories. The Cornell Ranger [90] has the Table 1.1: Mechanical Cost of Transport of Bipeds Bipeds Robots Humans Honda Asimo Rabbit Huang et al. - Optimal Trajectory Twente Dribbel MSU Synthetic Wheel Mabel TU Delft Denise MIT Spring Flamingo Cornell Collins 3D Cornell Ranger Dynamite Walking 13 Cmt 1.600 0.380 0.360 0.220 0.140 0.100 0.080 0.070 0.055 0.040 0.040 0.050 Ref. [86] [87] [88] [89] [1] [87] [86] [86] [86] [90] [83] [91] same Cmt as that of Dynamite but it is an active biped. It has a stable gait on relatively flat ground and is highly efficient by design; it was recently recorded to walk 65 km on a single charge. The Cornell Collins 3D [86] is a three-dimensional biped which achieves yaw control through mechanical coupling between its legs and opposite arms. Like the Ranger, it stores energy in a spring during the stance phase and releases it during push-off. It is fairly efficient in terms of energy consumption but it lacks gait stability. The MIT Spring Flamingo is a seven link planar biped with six actuated degrees-of-freedom. The feet and legs are slender compared to the upper body and are driven by series elastic actuators. A virtual model control [95] algorithm was used to make the robot walk over flat terrain and inclines. Denise [96], developed at TU Delft, has five degrees-of-freedom: one at the hip, two at the knees, and two at the ankles. Only the hip joint is actuated, it is driven by two antagonistic pairs of McKibben muscle-like actuators. Denise has arms but they do not add degrees of freedom since they are mechanically linked to the opposing thighs with belts. Similar to Denise, Dribble [89] is underactuated and has four degrees-of-freedom with a single active hip joint. The knee joints are passive but they have active locking mechanisms. It is important to note that both the gait and the controller used to control the gait have significant effect on the energy consumption of biped robots. Rabbit [97] and Mabel [87] are controlled using the Hybrid Zero Dynamics method [98] but they differ widely in their construction and implementation. Rabbit is a conventional stiff robot with highly-geared servomotors whereas Mabel is a compliant robot. A significantly lower Cmt for Mabel can be attributed to its ability to store energy during the gait cycle due to its compliance, and the ability of the controller to exploit this compliance. Huang et al. [88] presented the behaviors of a planar three-link biped for two different controllers. The first controller, implemented experimentally, uses an optimization method to determine optimal trajectories 14 for the boundary conditions of a step. The second controller, described as the self-excited controller, imposes constraints on the degrees-of-freedom of the robot and generates an oscillatory behavior of the system for walking. The self-excited controller results in a lower cost of transport but its Cmt value is not listed in Table 1.1 since it was not implemented experimentally. The MSU Synthetic-Wheel biped [1] has a moderate Cmt with its symmetric gait. An optimal gait [99], that eliminates impulsive forces at the time of swing-leg touchdown, offers to provide a significantly lower Cmt but it has not been implemented in experiments. Among other methods to improve energy efficiency, Lohmeier [100] proposed to increase limb stiffness, move the center of mass of the biped higher, and minimize limb inertia; Asano and Luo [101] and Wisse et al. [102] designed efficient mechanisms such as the bisecting hip mechanism; and Gomes and Ruina [103] proposed to eliminate energy losses due to impact through velocity matching. 1.7 Output Tracking Control for Non-Minimum-Phase Systems The output tracking problem for linear systems has been studied for more than four decades. It deals with finding a control input such that the internal dynamics is stable and the system output asymptotically tracks a desired trajectory. For linear Minimum-Phase (MP) systems, the tracking problem can be easily solved, however, it is difficult to achieve both asymptotic tracking and good transient performance in the case that the system is Non-Minimum-Phase (NMP). If the reference signals to be tracked are generated by an unstable external autonomous linear system, i.e., an exosystem, it is referred to as the output regulation problem or the 15 servomechanism problem [104, 105, 106]. The solution to the output regulation problem for linear systems relies on solving a set of two linear matrix equations which are not solvable if an eigenvalue of the exosystem is identical to a zero of the linear plant [107, 108, 109]. This makes the class of exosystem-generated reference signals even more limited for NMP systems with Right-Half-Plane (RHP) zeros. Inversion-based control techniques have been proposed for output tracking with specified desired trajectory. For a given desired output, the inversion problem is solved to find the desired state variables and the desired control input. The tracking control input is comprised of the desired input as a feedforward component and a stabilizing input as a feedback component. For NMP systems, the inversion techniques in [110] results in causal but unbounded solutions. This problem was resolved using stable inversion [111] for which the control input is non-causal and requires preactuation. Later, Weng and Chen [112] proposed a method based on causal stable inversion by requiring the desired trajectory to have a compact support. Stable inversion techniques cannot be applied to systems with non-hyperbolic internal dynamics, which correspond to imaginary-axis zeros for a linear system. For such cases, approximate stable inversion [113, 114] was proposed to achieve approximate output tracking. For nonlinear systems, Grizzle et al. [115] proved that asymptotically stable zero- dynamics is a necessary condition for asymptotic tracking of desired output trajectories. For NMP systems, only approximate output tracking methods have been developed; in these methods the output is redefined or the unstable internal dynamics is approximated by stable dynamics [116, 117, 118]. The unstable zero dynamics of NMP systems also limits the best achievable transient performance for the output tracking problem. The ideal performance for the cheap control problem is zero for MP systems; it is limited to the energy needed to 16 stabilize the internal dynamics for NMP systems [119, 120, 121]. Discrete-time systems with unstable zeros have the same control limitations as those of NMP continuous-time systems [122] and therefore, control methods based on zeroth-order hold can only achieve approximate output tracking for such systems [123, 124, 125]. By replacing the traditional zeroth-order hold, multirate control [126, 127] was proposed to overcome the problem of unstable zeros. Using a multirate feedforward control, Fujimoto et al. [128] developed a control method for perfect tracking in discrete-time systems. 1.8 Scope of Dissertation The work in this dissertation is organized as follows: In Chapter 2, we provide the experimental results for the swing-up control of the Pendubot to verify an impulsive control algorithm developed earlier in [2]. We develop an energy-based algorithm using impulsive inputs for swing-up control of the Acrobot in Chapter 3. An extension to current swing-up algorithms to include impulsive inputs is also discussed. In Chapter 4, impulsive inputs are used to enlarge the region of attraction of a desired equilibrium for underactuated systems. Several case studies of underactuated mechanical systems are studied to show the effectiveness of the algorithm. In Chapter 5, we use impulsive inputs to propose a safe fall algorithm for Humanoid robots. In particular, we design an impulsive control algorithm to change the fall direction of a six-dof robot to minimize damage to the surrounding objects. In Chapter 6, a disturbance rejection algorithm is developed for the Synthetic-Wheel biped [1]. Impulsive inputs are used to cancel the adverse effects of external impulsive disturbances and ground impact during the leg interchange. In Chapter 7, Energy-conserving gaits with ideal zero-energy consumption are generated for a class of underactuated bipeds by avoiding the 17 ground impact at the leg interchange. In Chapter 8, we propose an intermittent output tracking control method for linear non-minimum-phase systems. The switched inputs used in the control method can be approximated by impulsive inputs in the limit. Finally, we provide some concluding remarks and the directions for the future work in Chapter 9. 18 Chapter 2 Swing-up Control of Pendubot: Experimental Validation 2.1 Introduction Although impulsive inputs which are ideally modeled as Dirac-delta functions, can not be implemented, they can be estimated such that the same impulsive effects are resulted in the state variables. Using numerical simulations, Albahkali et al [2] demonstrated that impulsive inputs can be estimated by large-gain continuous inputs in their swing-up control algorithm for the Pendubot. The large-gain inputs are derived from dynamics which result in a fast convergence of the system velocities to their new values. It is shown that as the gains are chosen larger, the effects on the system dynamics get closer to the ideal impulsive case. In this chapter, the impulsive control algorithm in [2] is experimentally verified using a Pendubot set-up. The experimental results confirms that the impulsive inputs cab be closely estimated by continuous inputs and therefore, implemented using standard hardware. This chapter is organized as follows: In Section 2.2, we briefly review the swing-up control algorithm of the Pendubot in [2]. The experimental set-up and results are presented in Section 2.3 and concluding remarks are provided in Section 2.4. 19 2.2 Background: Swing-up Algorithm for Pendubot We briefly summarize the algorithm proposed by Albahkali, et al. [2] for swing-up control of the pendubot, shown in Figure 2.1: 1. Initialization: Linearize the dynamic equations of the pendubot about the desired ˙ ˙ equilibrium (θ1 , θ1 , θ2 , θ2 ) = (π/2, 0, 0, 0) and design a linear controller to render the equilibrium locally asymptotically stable. 2. First link swing-up: Drive the first link from its initial configuration to any configu˙ ration where θ1 satisfies (π/2 − α) ≤ θ1 ≤ (π/2 + α), θ1 = 0 for some small positive angle α. 3. Swing-up control of the second link: Conduct rest-to-rest maneuvers of the first link about the vertically upright configuration satisfying (π/2 − α) ≤ θ1 ≤ (π/2 + α). Use a continuous torque to drive the first link from rest and stop it with an impulsive torque when it reaches the boundary of the interval [(π/2 − α), (π/2 + α)]. It was shown [2] g θ2 d2 l1 l2 c.m . Y d1 θ1 . c.m X τ1 Figure 2.1: The pendubot is shown in an arbitrary configuration. 20 that this maneuver can increase the energy of the second link. Continue rest-to-rest maneuvers till the energy of the second link is approximately equal to Edes . 4. Stabilization: When the second link reaches its highest potential energy configuration, the pendubot configuration will be within the region of attraction of the desired equilibrium. Invoke the linear controller, designed in the first step, to stabilize the desired equilibrium. 2.3 Experimental Setup and Results Experiments were done with the P-2 pendubot system, produced by Mechatronic Systems, Inc. The kinematic and dynamic parameters of the pendubot are as follows: m1 = 1.0367 kg m2 = 0.5549 kg l1 = 0.1508 m l2 = 0.2667 m d1 = 0.1206 m d2 = 0.1135 m I1 = 0.0031 kgm2 I2 = 0.0035 kgm2 (2.1) The shoulder joint was driven by a DC motor, powered by an amplifier operating in the current mode. Both joints have encoders for joint angle sensing. The pendubot was interfaced with a dSpace DS1104 board and the controller was implemented in the Matlab/Simulink environment for real-time control. The impulsive torque in the third step of the algorithm was realized by applying a high-gain braking torque, discussed in [2]. The experiments were performed with α = 0.175 rad ≈ 10◦ ; the results are shown in Figure 2.2. It can be seen from Figure 2.2 that the second step of the algorithm commences at around t = 1.1 sec and the first link reaches its upright configuration at around t = 3.0 21 sec. The swing-up of the second link (third step) is initiated at approximately t = 3.2 sec. The pendubot configuration reaches the region of attraction of the desired equilibrium immediately after t = 7.0 sec; the linear controller is then invoked for stabilization. The rest-to-rest maneuvers of the first link in the third step of the algorithm can be seen from the subplot of θ1 in Figure 2.2. The two horizontal dotted lines in this plot roughly mark the boundary of oscillation of the first joint. The amplitude of this oscillation is approximately 13.5◦ . This is larger than the value of α ≈ 10◦ and this indicates that the braking torque was insufficient to stop the joint immediately at the boundary of the interval [(π/2 − α), (π/2 + α)]. In fact, it can be seen from the plot that θ1 overshot the boundary significantly during one maneuver. A larger gain for the braking torque could 2 a 0 ◦ −2 27 0 −2 b −6 5 0 −5 10 0 −10 5 0 −5 2 0 −2 θ1 (rad) θ2 (rad) c ˙ θ1 (rad/s) d ˙ θ2 (rad/s) e τ1 (N.m) f E(J) 0 step 1 2 step 2 4 time (s) step 3 6 8 10 step 4 Figure 2.2: Experimental validation of the impulsive control algorithm [2] for pendubot swing-up. 22 not be generated because of the limitation of the power amplifier but this did not pose any significant problem in swing-up of the second link. Although it is possible to design electrical circuits that enable the motor to apply torques that closely approximate an impulsive torque, the experimental results clearly indicate that regular power amplifiers are sufficient for simple implementations of impulsive control. The large peaks in the subplot of τ1 (in Figure 2.2) in the third step of the algorithm correspond to the impulsive braking torques. Each peak is followed by an oscillation in the torque profile - this is due to amplification of measurement noise. The impulsive braking ˙ torque is computed based on θ1 , which is obtained through differentiation of the measurement ˙ of θ1 . The process of differentiation amplifies the noise, which is visible in the plot of θ1 , and results in the torque ripple. ˙ ˙ The subplots of θ1 and θ2 indicate that the joint velocities undergo large changes (relatively speaking) at the time of application of the impulsive torques - this is consistent with the theory of impulsive control. It can be seen from the figure that the energy of the system also undergoes large changes at the time of application of the impulsive torques. Overall, the experimental results are consistent with those obtained through numerical simulations by Albahkali, et al. [2]. 2.4 Conclusion The swing-up algorithm of the Pendubot using impulsive inputs was experimentally validated. Impulsive inputs were approximated by large-gain continuous inputs and are derived from a fast dynamics which result in desired jumps in velocities. Although the impulsive inputs were assumed to provide instantaneous braking of the first link, the experimental 23 results indicate that the impulsive torques need not to be Dirac delta functions and the swing-up algorithm is effective even when the magnitude of the impulsive input is bounded and its time support is not infinitesimal. The magnitude of the impulsive inputs are however larger than the continuous torque required by the algorithm. This should not be seen as a problem since actuators such as motors can apply substantially larger torques than the maximum continuous torque over small time intervals. This is referred to as the peak torque [129] and it can be twice to ten times larger than the maximum continuous torque for different motors. The experiment results match the simulation results and this confirms that impulsive inputs are indeed implementable in mechanical systems using regular actuators and amplifier circuits. 24 Chapter 3 Swing-up Control of Acrobot 3.1 Introduction In this chapter, we present an energy-based method for swing-up control of the Acrobot. The method relies on the use of impulsive inputs together with continuous inputs and was used earlier to develop an algorithm for swing-up of the Pendubot [2]. The algorithm developed for the Pendubot can be directly used for swing-up of the Acrobot but it suffers from slow convergence. A modified algorithm, based on an additional set of impulsive inputs, provides fast convergence; also, it is not sensitive to the choice of control gains and does not impose restrictions on the initial conditions. This algorithm is presented in Section 3.3. The algorithm presented in Section 3.3 is based on the energy of the system and can therefore be viewed as a modification of Spong’s algorithm [19, 20] through inclusion of impulsive inputs in the set of admissible controls. This approach to enlarging the set of admissible controls is general and can be used to modify other algorithms in the literature. This is demonstrated in Section 3.4 through modification of the algorithms proposed by Xin and Kaneda [3] and Mahindrakar and Banavar [22], which are based on the work by Kolesnichenko and Shiriaev [21]. It is shown that impulsive inputs allow us to enlarge the set of continuous inputs for swing-up, and thereby additionally enlarge the set of admissible controls. This provides greater flexibility in control design. This chapter is organized as follows. In Section 3.2 we review the dynamic model of the 25 Acrobot and derive expressions for jump in velocities and change in energy for impulsive elbow torque. In Section 3.3 we present an algorithm for swing-up that does not suffer from limitations. In Section 3.4 we demonstrate the flexibility in control design for Acrobot swing-up by including impulsive inputs in the set of admissible controls. Section 3.5 provides the concluding remarks. 3.2 3.2.1 Acrobot Dynamics and Impulsive Effects Equations of Motion Consider the Acrobot in Figure 3.1. Assuming no friction in the joints, the equation of motion can be written as [19] ¨ ˙ ˙ M(θ)θ + N(θ, θ)θ + G(θ) = T (3.1) θ = [θ1 θ2 ]T , (3.2) where T = [0 τ2 ]T ˙ and M(θ), N(θ, θ), and G(θ), given by the expressions   q2 + q3 C2   q1 + q2 + 2q3 C2 M(θ) =   q2 + q3 C2 q2     ˙2 ˙1 + θ2 ) ˙ −(θ  −θ   q4 C1 + q5 C12  ˙ G(θ) = g  N(θ, θ) = q3 S2    , ˙1 q5 C12 θ 0 (3.3) are the inertia matrix, matrix containing the Coriolis and centrifugal forces, and vector of 26 g θ2 d1 θ1 d2 l1 l2 c.m . Y τ2 . c.m X Figure 3.1: The Acrobot in an arbitrary configuration. gravity forces, respectively. In Eq.(3.3), Sk = sin θk , Ck = cos θk , for k = 1, 2, C12 = cos(θ1 + θ2 ), g is the acceleration due to gravity and qj , j = 1, 2, · · · , 5 are some constants with the following expressions 2 q1 =m1 d2 + m2 l1 + I1 , 1 q3 =m2 l1 d2 , q2 = m2 d2 + I2 2 q4 = m1 d1 + m2 l1 , q5 = m2 d2 (3.4) where mk is the mass of the k-th link, dk is the distance between the k-th joint and its center-of mass, lk is the length of the k-th link and Ik is the mass moment of inertia of the k-th link, for k = 1, 2. 3.2.2 Holding Torque ˙ We compute the torque required to hold the second link fixed, i.e., maintain θ2 = 0. By ˙ ¨ substituting θ2 = θ2 = 0 in Eq.(3.1), we get 27          q4 C1 +q5 C12   0   q1 +q2 +2q3 C2  ¨  0   =  +g   θ1 +   ˙2 τ2h q5 C12 q3 S2 θ1 q2 +q3 C2 (3.5) ¨ Eliminating θ1 from the two equations in Eq.(3.5), the holding torque at the second joint can be expressed as follows ˙2 τ2h = q3 S2 θ1 + 3.2.3 g [(q +q C )q C −(q2 +q3 C2 )q4 C1 ] q1 +q2 +2q3 C2 1 3 2 5 12 (3.6) Impulsive Torque for Sudden Change in Velocity ˙ Consider the action that results in exponential convergence of the second link velocity, θ2 , ˙ to a desired value, θ2des , i.e., ¨ ˙ ˙ θ2 = −k1 (θ2 − θ2des ), k1 > 0 (3.7) ˙ where k1 is a constant that determines the rate of convergence of θ2 . To compute the torque required for this action, we multiply Eq.(3.1) with the inverse of the inertia matrix to get     ¨ 1  −(q2 +q3 C2 )τ2 +h1   θ1    = 2 2 q1 q2 −q3 C2 ¨ (q1 +q2 +2q3 C2 )τ2 +h2 θ2 (3.8) where h1 and h2 are given by the expressions 2 ˙2 ˙ ˙ h1 = q2 q3 (θ1 + θ2 )2 S2 + q3 θ1 S2 C2 + g(q3 q5 C2 C12 − q2 q4 C1 ) 28 (3.9) 2 ˙ ˙ ˙2 h2 =−(θ1 + θ2 )2 (q2 q3 +q3 C2 )S2 −(q1 +q3 C2 )q3 θ1 S2 −g {q3 q5 C2 C12 −(q2 +q3 C2 )q4 C1 +q1 q5 C12 } (3.10) Substitution of Eq.(3.7) into the second equation of (3.8) gives τ2 = −1 ˙ ˙ k (θ − θ )(q q −q 2 C 2 )+h2 q1 +q2 +2q3 C2 1 2 2des 1 2 3 2 (3.11) For a large value of gain k1 , the torque in Eq.(3.11) will be impulsive and converge the second link velocity to the desired velocity in a short period of time. As a special case, consider the impulsive torque that results in instantaneous stopping of the second link. The impulsive braking torque for the second link is obtained by substituting ˙ θ2des = 0 in Eq.(3.11), namely, τ2b = −1 2 2 ˙ k1 θ2 (q1 q2 − q3 C2 ) + h2 q1 + q2 + 2q3 C2 (3.12) When the second joint comes to rest, the braking torque is equal to the holding torque. This can be verified from Eq.(3.6) and Eq.(3.12). 3.2.4 Change in Velocity due to Impulsive Torque The application of an impulsive torque in the active coordinate θ2 results in velocity jumps in both active and passive coordinates. The relationship between the jumps in velocities can be derived from Lagrange’s equations using the approach presented in [1]. Consider the integral of the equations of motion in Eq.(3.1) over the short interval of time ∆t during 29 which the impulsive torque acts, i.e. ∆t ∆t ¨ ˙ ˙ M(θ)θ + N(θ, θ)θ + G(θ) dt = 0 T dt (3.13) 0 The above equation can be rewritten as: ∆t ˙ ˙ M(θ)∆θ+N(θ, θ)∆θ + ∆t G(θ) dt = 0 ∆t ˙ ∆θ ¨ θ dt, T dt, 0 ∆t ∆θ 0 ˙ θ dt (3.14) 0 Since ∆t ≈ 0 and the configuration of the Acrobot does not change during this time, i.e. ∆θ = 0, Eq.(3.14) will be simplified to ˙ M(θ)∆θ = Timp (3.15) where Timp = [0 τ2i ]T is the vector of applied impulses. From Eq.(3.3), the above equation can be decomposed into the following equations: ˙ + ˙− ˙ + ˙− [q1 + q2 + 2q3 C2 ](θ1 − θ1 )+[q2 + q3 C2 ](θ2 − θ2 ) = 0 ˙+ ˙− ˙+ ˙− [q2 + q3 C2 ](θ1 − θ1 )+q2 (θ2 − θ2 ) = τ2i where (3.16) (3.17) ˙− ˙− θk and θk denote the amgular velocities of the k-link, immidiately before and after the application of impulse, respectively. Eq.(3.16) can be written as ˙+ ˙ − ˙+ ˙− [q1 + q2 + 2q3 C2 ](θ1 − θ1 ) = −[q2 + q3 C2 ](θ2 − θ2 ) 30 (3.18) If the velocity of the second link after the impulse is equal to some desired velocity, i.e., ˙+ ˙ θ2 = θ2des , Eq.(3.18) can be used to obtain the velocity of the first link after the impulse: q2 + q3 C2 ˙+ ˙− ˙ ˙− θ1 = θ1 − (θ2des − θ2 ) q1 + q2 + 2q3 C2 3.2.5 (3.19) Change in Energy due to Impulsive Torque The configuration, and hence the potential energy, of the Acrobot will not change over the small period of time that the impulse is applied. Therefore, the change in total energy of the system is due to the change in the kinetic energy, which can be expressed as ∆E = ∆K1 + ∆K2 (3.20) where ∆K1 and ∆K2 denote the change in kinetic energy of the first and second link, and are given by the relations 1 ˙+ ˙− ∆K1 = (I1 + m1 d2 )[(θ1 )2 − (θ1 )2 )] 1 2 1 1 ˙+ ˙+ ˙− ˙− ∆K2 = I2 (θ1 + θ2 )2 − I2 (θ1 + θ2 )2 2 2 1 2 ˙+ ˙+ ˙+ ˙+ ˙+ ˙+ + m2 l1 (θ1 )2 + d2 (θ1 + θ2 )2 + 2l1 d2 θ1 (θ1 + θ2 )C2 2 2 1 2 ˙− ˙ − ˙− ˙− ˙− ˙− − m2 l1 (θ1 )2 + d2 (θ1 + θ2 )2 + 2l1 d2 θ1 (θ1 + θ2 )C2 2 2 (3.21) (3.22) If we assume the velocity of the second link after the impulse to be equal to the desired ˙+ ˙ velocity, i.e., θ2 = θ2des , the change in total energy can be obtained by substituting Eq.(3.19) 31 into Eqs.(3.21) and (3.22) as follows: ∆E = 2 2 q1 q2 − q3 C2 ˙ ˙− (θ2des )2 − (θ2 )2 2(q1 + q2 + 2q3 C2 ) (3.23) Equation Eq.(3.23) implies that the change in total energy of the system depends on the ˙ magnitude of θ2des , i.e., ˙ ˙− ∆E > 0 if |θ2des | > |θ2 | ˙ ˙− ∆E = 0 if |θ2des | = |θ2 | (3.24) ˙ ˙− ∆E < 0 if |θ2des | < |θ2 | ˙ If an impulsive braking torque is used, i.e., θ2des = 0, the change in energy is negative and is given by the relation ∆E = − 3.3 3.3.1 2 2 q1 q2 − q3 C2 ˙ (θ − )2 2(q1 + q2 + 2q3 C2 ) 2 (3.25) Swing-Up Control of the Acrobot Rest-to-Rest Maneuvers of the Second Link Consider a maneuver in which the second joint starts from rest under the application of a continuous torque and is brought back to rest through the application of the impulsive braking torque in Eq.(3.12), where gain k1 has a large value. Taking into account the loss of energy due to impulsive braking given by Eq.(3.25), the change in energy during each 32 rest-to-rest maneuver can be computed as follows ∆E = ˙ τ2c θ2 dt − 2 2 q1 q2 − q3 C2 ˙ (θ − )2 2(q1 + q2 + 2q3 C2 ) 2 (3.26) If the torque in Eq.(3.26) is assumed to have the form ¨ τ2c = k2 θ2 , k2 > 0 (3.27) ˙− the work done by the torque to increase the second link velocity from zero to θ2 , which is the velocity prior to application of the impulsive brake, is expressed as follows ˙ τ2c θ2 dt = ˙ ¨ k2 θ2 θ2 dt = ˙− θ2 0 1 ˙− ˙ ˙ k2 θ2 dθ2 = k2 (θ2 )2 2 (3.28) The total change in energy due to the rest-to-rest maneuver can now be obtained by substituting Eq.(3.28) into Eq.(3.26) ∆E = 2 2 q1 q2 − q3 C2 1 ˙− k2 − (θ2 )2 2 (q1 + q2 + 2q3 C2 ) (3.29) From Eq.(3.29), it can be seen that the gain k2 can be chosen to add or remove energy, namely 2 2 q1 q2 − q3 C2 (q1 + q2 + 2q3 C2 ) 2 2 q1 q2 − q3 C2 k2 < (q1 + q2 + 2q3 C2 ) k2 > 33 ⇒ ∆E > 0 ⇒ ∆E < 0 (3.30) Using Eq.(3.8), the torque expression in Eq.(3.27) can be written as: τ2c = 2 2 q1 q2 − q3 C2 k2 h2 − k2 (q1 + q2 + 2q3 C2 ) (3.31) The above expression for the torque will not be singular if k2 is chosen based on Eq.(3.30). 3.3.2 Acrobot Swing-up Using Rest-to-Rest Maneuvers 3.3.2.1 Swing-up Algorithm A swing-up algorithm for the Acrobot is proposed using rest-to-rest maneuvers of the second link. This algorithm is similar to the one proposed for the Pendubot [2] in the sense that energy of the system is increased through rest-to-rest maneuvers, and impulsive control inputs are used as braking torques. The algorithm is comprised of the following three steps: 1. Initialization: • Linearize the dynamic equations of the Acrobot in Eq.(3.1) about the equilibrium ˙ ˙ configuration (θ1 , θ1 , θ2 , θ2 ) = (π/2, 0, 0, 0). • Design a linear controller to render the equilibrium locally asymptotically stable. Let RA define the region of attraction of the equilibrium. • Choose a small positive constant δ, such that the system configuration lies inside ˙ ˙ RA when θ1 ≈ θ2 ≈ 0 and |P − Edes | < δ, where P is the potential energy and Edes is the total energy of the Acrobot at its upright equilibrium. 2. Increasing the energy of the system: 34 • Conduct rest-to-rest maneuvers to increase E to its desired value, Edes . The restto-rest maneuvers will be implemented with θ2 satisfying −α ≤ θ2 ≤ α, where α is a small positive angle. Using Eq.(3.30), the control gain for the torque expression in Eq.(3.31) is chosen such that ∆E > 0, i.e., the energy is increased through the rest-to-rest maneuvers. In particular, the following procedure will be adopted. The holding torque in Eq.(3.6) will be applied to hold the second link fixed. To initiate motion of the second link in the positive (counter-clockwise) direction, the torque expression in Eq.(3.31) will be used when it is greater than τ2h . To initiate motion in the negative (clockwise) direction, the torque expression in Eq.(3.31) will be used when it is less than τ2h . As the second link approaches the boundary of the interval [−α, α], the braking torque τ2b in Eq.(3.12) will be used with a large value of gain k1 to quickly stop the motion of the second link. • Terminate the rest-to-rest maneuvers when the condition |E−Edes | < δ is satisfied and apply the holding torque in Eq.(3.6) to keep the second link fixed. 3. Stabilization: ˙ With −α ≤ θ2 ≤ α, θ2 = 0, and |E − Edes | < δ, the first link will behave like a pendulum and reach its configuration with maximum potential energy in finite time. ˙ ˙ When the Acrobot has maximum potential energy, i.e. |P −Edes | < δ and θ1 = θ2 = 0, its configuration will be inside RA . Invoke the linear controller, designed in the first step of the algorithm, to stabilize the equilibrium. Remark 3.3.1. For the rest-to-rest maneuvers of the second link in step 2, τ2c in Eq.(3.31) has to be greater (or less) than τ2h to initiate the motion of the second link in the counter35 ˙ clockwise (or clockwise) direction when θ2 = −α (or θ2 = +α) and θ2 = 0. We show that ˙ ˙ this condition can be satisfied by choosing proper values for α. For θ1 = θ2 = 0, Figure 3.2 shows the regions in θ1 -θ2 space where Π = (τ2c − τ2h ) is positive and negative. If the second link is held fixed at θ2 = +α or θ2 = −α, the Acrobot will oscillate like a pendulum and θ1 ˙ will oscillate between its maximum and minimum values where θ1 = 0. The amplitude of oscillation of θ1 will depend on the energy of the Acrobot and the maximum and minimum values of θ1 will correspond to Acrobot configurations that lie at the end points of the dashed lines in Figure 3.2. If the end point of the dashed lines are on opposite sides of the curve Π = 0, as shown in Figure 3.2, it will be possible to initiate the rest-to-rest maneuvers in the desired direction. To ensure that the end points are on opposite sides of the curve Π = 0, we will have to properly choose the value of α, i.e., choose a small value when the Acrobot has low energy. Remark 3.3.2. It may be desirable to choose a fixed value of α, say α, for all the rest-to-rest ¯ maneuvers during swing-up. It follows from Remark 3.3.1 that it will be possible to conduct 0.8 Π=0 θ2 (rad) 0.4 Π>0 α α ¯ 0.0 η ¯ Π<0 -0.4 -0.8 −3π/2 −α Π=0 −π/2 −π θ1 (rad) Π=0 0 π/2 ˙ ˙ Figure 3.2: Plot showing θ1 -θ2 space where Π is positive and negative for θ1 = θ2 = 0. 36 ¯ the rest-to-rest maneuvers successfully if the initial energy of the Acrobot is greater than P , ¯ where P is the potential energy of the Acrobot at the configuration (θ1 , θ2 ) = (α, η ) and α ¯ ¯ ¯ and η are shown in Figure 3.2. ¯ 3.3.2.2 Numerical Simulation We present numerical simulation results using kinematic and dynamic parameters of the Acrobot from Xin and Kaneda [3]: m1 = 1 kg, l1 = 1 m, d1 = 0.5 m, I1 = 0.083 N.m2 m2 = 1 kg, l2 = 2 m, d2 = 1 m, I2 = 0.33 N.m2 (3.32) ˙ ˙ The initial conditions were chosen as (θ1 , θ2 , θ1 , θ2 ) = (−π/4, 0, 0, 0) and the algorithm parameters were set to be δ = 0.1 J and α = 0.17 rad. For the parameters in Eq.(3.32), Edes was computed to be 49.05 J. The simulation results are shown in Figure 3.3. By choosing the gain k2 in Eq.(3.31) appropriately, the total energy of the Acrobot is increased in each rest-to-rest maneuver. The peaks in the torque plot represent impulsive braking torques used in the rest-to-rest maneuvers. The main drawback of the algorithm in Section 3.3.2.1 is that it takes a long time for swing-up. For the simulation in Figure 3.3, the algorithm is switched to the linear controller at about t = 115 sec, which is shown by vertical dotted line. This time is much longer than swing-up times reported in the literature. For example, for the same kinematic and dynamic parameters and lower initial energy level, the algorithm in [3] achieved swing-up in about 8 sec. The swing-up time required by our algorithm will be longer if the initial energy of the Acrobot is lower due to the constraint on the choice of α - see Remark 3.3.1. In the next 37 θ1 (rad) ˙ ˙ θ2 (rad/s) θ1 (rad/s) θ2 (rad) 0.17 0.0 -0.17 5 0 -5 1 0 -1 τ2 (N.m) E(J) 3.5 0.0 -3.5 a b c d 50 e 0 10 0 -10 f 0 20 40 60 80 100 120 time (sec) Figure 3.3: Simulation results for swing-up of the Acrobot using rest-to-rest maneuvers of the second link. section, we modify the algorithm to provide faster swing-up. 3.3.3 Increasing Energy of the Acrobot Using Impulsive Inputs 3.3.3.1 Swing-up Algorithm The algorithm in Section 3.3.2 used impulsive inputs for braking, which reduces the energy of the system. For faster swing-up, we propose to use impulsive inputs to increase the energy of the Acrobot. Specifically, we apply the impulsive torque in Eq.(3.11) with large gain k1 ˙ and θ2des chosen such that the energy of the system increases based on the condition in 38 Eq.(3.24). If the velocity of the second link after the impulse is chosen as, ˙ ˙− θ2des = −λ θ2 , ˙− θ2 = 0, λ>1 (3.33) the second link will have a velocity greater in magnitude and opposite in direction to that prior to the impulse, and the energy of the system will be higher. The change in sign of the velocity is needed to keep the second link angle bounded within some predefined range while the energy is raised to the desired level. Most of the energy gained from application of impulsive torques will be in the form of kinetic energy and needs to be converted into potential energy for swing-up. To convert kinetic energy into potential energy, the torque is set to zero and the second joint is made free. This keeps the total energy of the system constant and allows kinetic energy to be converted into potential energy. The rest-to-rest maneuvers are used subsequently to regulate the energy of the Acrobot to the desired level Edes . The modified algorithm is stated next: 1. Initialization: Same as step 1 in Section 3.3.2.1. 2. Increasing the energy of the system: ˙ • At the initial time, if θ2 = θ2 = 0, apply the impulsive torque in Eq.(3.11), where ˙ θ2des is chosen to be nonzero. From Eq.(3.24), we know that this will increase the energy of the Acrobot. • Using Eq.(3.33), apply the impulsive torque in Eq.(3.11) when the second link approaches a bound of the interval [−γ, γ] from inside, where γ is some positive angle - see Figure 3.4. Note that γ can always be chosen such that the second link approaches the bound of the interval [−γ, γ] for any initial conditions or initial 39 ˙ velocity of the second link, θ2des . Applying these impulsive torques will increase the energy of the system, mainly in the form of kinetic energy. Continue this process till E ≥ Edes . • Release the system by setting τ2 = 0. This will cause the links to swing freely resulting in exchange between kinetic and potential energies while the total energy remains constant. • Stop the second link using the impulsive braking torque in Eq.(3.12) when the ˙ ¨ potential energy of the Acrobot reaches its local maxima, i.e, P = 0, P < 0. The braking action will reduce the total energy of the system according to Eq.(3.25) but the potential energy will remain at its local maxima. 3. Regulation of energy to the desired level: • Same as step 2 in Section 3.3.2.1: Conduct rest-to-rest maneuvers with θ2 ∈ [−α, α]. If the second link is initially out of this bound, it will be brought inside by the first rest-to-rest maneuver. Depending on the current level of energy, the gain k2 is chosen from Eq.(3.30) to regulate (increase or decrease) E to Edes . • Terminate the rest-to-rest maneuvers when the condition |E−Edes | < δ is satisfied and apply the holding torque in Eq.(3.6) to keep the second link fixed. 4. Stabilization: Same as step 3 in Section 3.3.2.1. Remark 3.3.3. The above algorithm is not sensitive to the choice of control gains and does not impose restrictions on the initial conditions. It also provides fast swing-up of the Acrobot - this is demonstrated next through simulations. 40 Y X τ2 ˙+ ˙− θ2 = −λ θ2 +γ −γ Figure 3.4: Configuration of the Acrobot at the time when the impulsive torque is applied to increase the energy of the system. 3.3.3.2 Numerical Simulation We present numerical simulation results using kinematic and dynamic parameters of the Acrobot from Xin and Kaneda [3], given in Eq.(3.32). The parameters of our algorithm were chosen as: λ = 2.5, δ = 0.10 J, γ = 60◦ , α = 5◦ (3.34) ˙ ˙ The initial conditions were chosen at (θ1 , θ2 , θ1 , θ2 ) = (−π/2, 0, 0, 0)1 and the initial im˙ pulsive torque was chosen to result in θ2des = 3.0 rad/s. For the parameters in Eq.(3.32), Edes was computed as 49.05 J. The simulation results are shown in Figure 3.5. The total energy (solid line) is discretely increased through the application of impulsive torques. After five discrete applications of impulsive torque, the condition E ≥ Edes is satisfied and the second joint is released. The potential energy (dashed line) subsequently increases, and 1 This initial configuration corresponds to zero kinetic energy and lowest potential energy and cannot be handled by the algorithms in Xin and Kaneda [3] and Mahindrakar and Banavar [22]. 41 0 −4 8 4 0 10 0 −10 a θ1 (rad) b θ2 (rad) c ˙ θ1 (rad/s) ˙ θ2 (rad/s) 20 d 0 −20 40 0 10 0 −10 e E, P (J) f τ2c (N.m) 100 g 0 −100 0 2.42 3.86 2 step 2 4 time (sec) step 3 τ2i (N.m) 5.35 6 8 step 4 Figure 3.5: Simulation results for swing-up of the Acrobot using impulsive inputs to increase the energy. when it reaches its local maxima, the second link is stopped using an impulsive brake. This event, which occurs at approximately t = 2.42 sec, results in a drop in the total energy of the system. After application of the impulsive brake, the condition |E − Edes | < δ is not satisfied and therefore one rest-to-rest maneuver is used to regulate the energy to Edes . The rest-to-rest maneuver is completed at approximately t = 3.86 sec and the linear controller is invoked at t = 5.35 sec. In comparison to the algorithm proposed in [3], our algorithm requires shorter time for swing-up: 5.35 sec compared to 7.33 sec, and lower maximum continuous torque: 12.5 N.m compared to 20.0 N.m, despite starting from a lower initial energy configuration. 42 Numerical simulations results of the algorithm using Acrobots parameters in Mahindrakar and Banavar [22] can be found in [130]. 3.4 3.4.1 Modification of an Existing Method Existing Control Method Kolesnichenko and Shiriaev [21] developed a general methodology for control of underactuated systems; this method was used for Acrobot swing-up by Xin and Kaneda [3] and Mahindrakar and Banavar [22]. We first review this method and then modify it to include impulsive control inputs. Consider the dynamics of the Acrobot in the following form: x = f (x, τ2c ) ˙ (3.35) ˙ ˙ where x = (θ1 , θ2 , θ1 , θ2 )T and f (x, τ2c ) can be determined using Eq.(3.8). A continuous controller is designed to achieve the following objectives: lim θ2 = 0, t→∞ ˙ lim θ2 = 0, t→∞ lim E = Edes t→∞ (3.36) ˙ With θ2 = 0, θ2 = 0, and E = Edes , the Acrobot behaves like a pendulum and its trajectory converges to a heteroclinic orbit that passes through the desired upright equilibrium. To achieve the control objective, the following “positive storage function” is used: 1 2 ˙ ˙ kp θ2 + kd (θ2 )2 + ke (E − Edes )2 V (θ, θ) = 2 43 (3.37) where kp , kd , ke > 0. The following control input is chosen: τ2c = − ˙ ˙ kp θ2 + kd a(θ, θ) + φ(θ2 ) kd b(θ2 ) + ke (E − Edes ) (3.38) ˙ where a(θ, θ) and b(θ2 ) are given by the expressions ˙ a(θ, θ) = b(θ2 ) = 1 h 2 2 2 q1 q2 − q3 C2 1 (q + q2 + 2q3 C2 ) 2 2 1 q1 q1 − q3 C2 (3.39) and the condition kd min{b(θ2 )} > ke Edes θ2 (3.40) is satisfied to avoid singularity in the control torque in Eq.(3.38). Using Eq.(3.38), the derivative of the storage function is shown to be: ˙ ˙ ˙ ˙ V (θ, θ) = −θ2 φ(θ2 ) (3.41) ˙ ˙ where φ(θ2 ) is chosen such that V ≤ 0, and therefore the system is stable. In [3] and [22], ˙ the function φ(θ2 ) is chosen as: ˙ ˙ φ(θ2 ) = kc θ2 , kc > 0 (3.42) The invariant set for the closed-loop system is given by M = V0 ∪ Ω 44 (3.43) where V0 corresponds to the desired equilibrium and Ω is the finite set of undesirable equilibria: V0 = {x ∈ R4 | V (x) = 0} ˙ Ω = {x ∈ R4 | V (x) > 0 , V ≡ 0} (3.44) In [3] and [22], the gains kp , kd , and ke are chosen such that the equilibria in Ω are hyperbolic and V0 is attractive. This choice of gains ensures that the the control objective is achieved for all initial conditions except a set from which system trajectories do not converge to the desired equilibrium. 3.4.2 Modified Method Using Impulsive Inputs 3.4.2.1 Impulsive System We modify the control method in [3] and [22] to include impulsive inputs. Consider the following impulsive dynamical system: x(t) = f (x, τ2c ), t ∈ η(x) ˙ / ∆x(t) = H(x, τ2i ), t ∈ η(x) (3.45) where x, f (x, τ2c ) and τ2c are given in Eqs.(3.35) and (3.38), and correspond to the closedloop system dynamics of the Acrobot designed by Xin and Kaneda [3] and Mahindrakar and Banavar [22]. H(x, τ2i ) represents a jump in the states of the Acrobot due to the application 45 of the impulse τ2i . Using Eqs.(3.16) and (3.17), it can be shown H(x, τ2i ) = (x+ − x− ) q +q +2q C q +q C = 0 0 − 2 32 2 2 1 2 2 3 2 2 q1 q2 −q3 C2 q1 q2 −q3 C2 T τ2i (3.46) In (3.46), x+ and x− represent the right and left limits of x at the instant of the impulse, respectively. By choosing the impulse τ2i to be τ2i = − 2 2 q1 q2 − q3 C2 ˙ θ− q1 + q2 + 2q3 C2 2 (3.47) ˙+ we get impulsive braking which results in θ2 = 0. This can be verified from H(x, τ2i ), which can be rewritten as: H(x, τ2i ) = 0 0 T q2 + q3 C2 ˙− ˙− θ2 − θ2 q1 + q2 + 2q3 C2 (3.48) We choose the set η(x) in Eq.(3.45) to be: η(x) = {t ∈ [t0 , ∞) | x ∈ ρ(x)} (3.49) ˙ ¨ ρ(x) = {x ∈ R4 | θ2 θ2 < 0} (3.50) ˙ ¨ ρ⊥ (x) = {x ∈ R4 | θ2 θ2 ≥ 0} (3.51) where We define the set ρ⊥ (x) as 46 such that ρ(x) ∩ ρ⊥ (x) = ∅ (3.52) ˙+ Since θ2 = 0 from our choice of τ2i in Eq.(3.47), x+ lies on the boundary of ρ⊥ (x). By ˙ ¨ taking the derivative of θ2 θ2 with respect to time immediately following the impulse, we get ˙ ¨ ¨+ ˙+ + (3) = (θ+ )2 ≥ 0 ¨ D(θ2 θ2 ) = (θ2 )2 + θ2 θ2 2 (3.53) ˙ ¨ which means that θ2 θ2 is nondecreasing after the impulse. It follows from continuity that x will be in the interior of ρ⊥ (x) after an impulse and there exists a nonzero time interval ǫ > 0 before x can enter the set ρ(x). Therefore there is no instant with multiple impulses. 3.4.2.2 Stability of Impulsive System The jump in the positive storage function in Eq.(3.37) due to the impulse in Eq.(3.47) can be computed as follows ∆V = V (x+ ) − V (x− ) 1 1 ˙− = − kd (θ2 )2 + ke (E + −Edes )2 −(E − −Edes )2 2 2 1 1 ˙− = − kd (θ2 )2 + ke (E + +E − −2Edes ) ∆E 2 2 − )2 ˙ ˙− (θ2 ke (θ2 )2 =− kd b(θ2 )+ +ke (E + −Edes ) 2b(θ2 ) 4b(θ2 ) ˙− 1 kd b(θ2 )+ke(E + −Edes ) ˙− 2 ke (θ2 )4 (θ2 ) + 2 =− 2 b(θ2 ) 4b (θ2 ) 47 (3.54) where Eqs.(3.25) and (3.39) were used. From the condition in Eq.(3.40) we have kd b(θ2 ) + ke (E + − Edes ) >0 b(θ2 ) ⇒ ∆V ≤ 0 (3.55) The continuous control input given by Eqs.(3.38) and (3.42) and the impulsive input in Eq.(3.47) results in D+V ≤ 0 ∀t ∈ η(x) / ∆V ≤ 0 ∀t ∈ η(x) (3.56) where D + V denotes the upper right Dini derivative of V with respect to time. Using the stability theorems in [9] and [131], we can claim stability of the impulsive dynamical system in (3.45). To prove the asymptotic convergence to the desired equilibrium, we use the Invariance Principle for impulsive dynamical systems in [132] and [133]. To find the largest invariant set, we define the set ˙ Z ={x ∈ R4 | V = 0, x ∈ ρ⊥ (x)} ∪ {x ∈ R4 | ∆V = 0, x ∈ ρ(x)} (3.57) Using Eqs.(3.41), (3.50) and (3.54), it can be shown ˙ Z = {x ∈ R4 | θ2 = 0, x ∈ ρ⊥ (x)} ∪ ∅ 48 (3.58) Therefore, the largest invariant set M is given by ˙ M = {x ∈ R4 | θ2 ≡ 0, x ∈ ρ⊥ (x)} (3.59) which is the same set as in Eq.(3.43) obtained in [3] and [22]. The invariant set is not simply V0 ; it includes Ω which contains a finite set of equilibrium points. In [3] and [22], a set of conditions were imposed such that all continuous trajectories that do not start in Ω are attracted to V0 . For an impulsive dynamical system we have both continuous and discrete trajectories, and therefore, we must extend the results in [3] and [22] to include discrete trajectories. We need to find conditions such that discrete trajectories do not end in Ω. To this end, we define the set ˙ ¨ S(x) = {x ∈ R4 | θ2 θ2 = 0} (3.60) which is the boundary between ρ(x) and ρ⊥ (x). Let S1 (x) be the subset of S(x) from where trajectories enter ρ(x) from ρ⊥ (x), i.e., S1 (x) = {x ∈ S(x) | d ˙ ¨ (θ θ ) < 0} dt 2 2 (3.61) Each element of S1 is arbitrarily close to a point in ρ(x) from where a discrete trajectory originates and satisfies d ˙ ¨ ˙ (3) ¨ ˙ (3) (θ2 θ2 ) = θ2 θ2 + (θ2 )2 < 0 ⇒ θ2 θ2 < 0 dt 49 (3.62) In addition, define the set S2 as S2 (x) = {x ∈ S(x) | x + H(x, τ2i ) ∈ Ω} (3.63) Each element of S2 is arbitrarily close to a point in ρ(x) from where a discrete trajectory ends in Ω. The set S1 ∩ S2 contains all points from where trajectories of the closed-loop system end in Ω. We are now ready to state the additional condition that needs to be satisfied for asymptotic convergence of trajectories to the desired equilibrium. Theorem 3.4.1. Consider the impulsive dynamical system in Eq.(3.45) with continuous control given by Eqs.(3.38) and (3.42); and impulsive control given by Eq.(3.47) and applied at t ∈ η(x), defined in Eq.(3.49). If the positive gains kp , kd , ke , and kc are chosen such that ˙ (3) θ2 θ2 ≥ 0 ∀ x ∈ S2 (3.64) then no discrete trajectories will end in the set Ω of the invariant set M. Proof. From the definition of S1 in (3.61) it can be shown that the condition in (3.64) implies S1 ∩ S2 = ∅. Remark 3.4.1. The algorithms in [3] and [22] impose constraints on kp , kd and ke values to guarantee convergence of Acrobot trajectories to the upright equilibrium provided that certain initial conditions are avoided. The constraint in (3.64) has to be additionally satisfied for implementation of the modified approach that includes impulsive inputs, presented above. Although the modified approach imposes an additional constraint on the gains, the use of impulsive inputs provides us with greater flexibility in choosing the continuous control inputs. This will be illustrated in Section 3.4.3. 50 3.4.2.3 Numerical Simulation We provide simulation results for the modified impulsive controller presented above and compare them with the results of the continuous controller in [3]. The kinematic and dynamic parameters of the Acrobot are taken from [3] and are given here in Eq.(3.32). The initial conditions and gains were chosen as follows: ˙ ˙ (θ1 , θ2 , θ1 , θ2 ) = (−1.4, 0, 0, 0) 100 kp = 61.2, kd = 35.8, ke = 1, kc = 1+0.14 V (x) (3.65) The initial conditions and the gains kp , kd and ke in Eq.(3.65) are the same as those in [3] and they guarantee convergence of the Acrobot trajectory to the upright equilibrium under continuous control. To establish convergence using the impulsive control approach, we used the procedure in the Appendix A to determine the elements of S2 , defined in (3.63). For the initial conditions and gain values in Eq.(3.65), the positive storage function V (x) for all x ∈ S2 , was found to be larger than the value of V at the initial time. This implies that S2 = ∅ for positive-time solution trajectories and therefore the condition in Eq.(3.64) is trivially satisfied. The simulation results are shown in Figure 3.6. For the parameters in Eq.(3.32), Edes was computed to be 49.05 J. The solid and dashed lines in the plots for energy E, positive storage function V , and control input τ2 correspond to the modified impulsive controller and the controller in [3], respectively. Note that the torque plot for the modified controller includes both continuous and impulsive torques, i.e., τ2 = τ2c + τ2i . Both controllers achieve swing-up in approximately 8 sec, after which a linear controller is invoked for stabilization. 51 10 0 -10 5 0 -5 5 0 -5 10 0 -10 40 0 a θ1 (rad) b θ2 (rad) c ˙ θ1 (rad/s) d ˙ θ2 (rad/s) e E (J) 1000 0 V f 0 g -200 0 τ2 (N.m) 2 4 time (sec) 6 8 10 Figure 3.6: Simulation results for swing-up of the Acrobot: Comparison of the modified impulsive controller (solid line) and the continuous controller (dashed line) proposed by Xin and Kaneda [3]. Although the time required for swing-up is similar, the impulsive controller requires larger continuous torques and multiple impulses within a short interval of time. 3.4.3 Alternate Continuous Control Design 3.4.3.1 Control Design Consider the impulsive dynamical system in Eq.(3.45) with positive storage function given by Eq.(3.37). We choose to apply the same impulsive control input in Eq.(3.47) at time 52 instants given by Eq.(3.49) but modify the continuous control input in Eq.(3.38) as follows: τ2c = − ¨ ˙ kp θ2 + kd a(θ, θ) + kc θ2 , kd b(θ2 ) + ke (E − Edes ) kc > 0 (3.66) This results in ˙ ¨ D + V = −kc θ2 θ2 (3.67) From the definitions for η(x) and ρ(x) in Eq.(3.49) and Eq.(3.50), we have ˙ ¨ D + V = −kc θ2 θ2 ≤ 0 ∀t ∈ η(x) / (3.68) ¨ Using Eqs.(3.8) and (3.39), the expression for θ2 can be found and substituted into Eq.(3.66) to get τ2c = − ˙ kp θ2 + (kd + kc ) a(θ, θ) (kd + kc ) b(θ2 ) + ke (E − Edes ) (3.69) Since b(θ2 ) > 0 and kc , kd , ke > 0, the singularity in the torque expression in Eq.(3.69) can be avoided by simply choosing the gains as follows kd min{b(θ)} > ke Edes θ (3.70) The jump in V due to the impulse is the same as in Eq.(3.54) and therefore the stability of the system is guaranteed - this follows from our discussion in Section 3.4.2.2. The invariant set M is found to be the same as the set in Eq.(3.59), which is identical that derived in [3] and [22]. To prove the convergence to the desired upright configuration, we follow the same procedure as in Section 3.4.2.2. The gain values kp , kd , ke and kc need to be chosen such that the condition in Theorem 3.4.1 is satisfied in addition to the conditions 53 in [3] and [22]. 3.4.3.2 Numerical Simulation We compare the performance of the controller presented above with that of the continuous controller in [3]. The kinematic and dynamic parameters of the Acrobot are taken from [3] and are given in Eq.(3.32). The initial conditions and gains were chosen as: ˙ ˙ (θ1 , θ2 , θ1 , θ2 ) = (−1.4, 0, 0, 0) kp = 61.2, kd = 35.8, ke = 1, kc = 75.6 (3.71) The initial conditions and the gains kp , kd and ke are the same as those in [3] and they guarantee convergence of the Acrobot trajectory to the upright equilibrium under continuous control. To establish convergence using the impulsive control approach, we used the procedure in Appendix A to determine the elements of S2 , defined in Eq.(3.63). For the initial conditions and gain values in Eq.(3.71), the set S2 was found to be empty and therefore the condition in Eq.(3.64) is trivially satisfied. The simulation results are shown in Figure 3.7. As in the last simulation, the solid and dashed lines in the plots for energy E, positive storage function V , and control input τ correspond to our controller and the controller in [3], respectively. Our controller achieves swing up in approximately 7.33 sec. This is slightly faster than the controller in [3], which requires 8 sec. For both controllers, a linear controller was invoked for stabilization after swing-up. It can be noted that the number of impulses used by our controller are fewer and farther apart than those in Figure 3.6. Furthermore, the maximum continuous torque is approximately 15 N.m, which is lower than 20 N.m required by the controller in [3]. 54 0 θ1 (rad) -5 a 0 θ2 (rad) -4 b 0 -4 c 5 0 -5 d 40 e 20 0 1000 500 0 f 80 g 40 0 0 ˙ θ1 (rad/s) ˙ θ2 (rad/s) E (J) V τ (N.m) 2 4 time (sec) 6 8 10 Figure 3.7: Simulation results for swing-up of the Acrobot: Comparison of the controller in Section 3.4.3 (solid line) with the continuous controller (dashed line) proposed by Xin and Kaneda [3]. The impulsive control approach allows us to replace the continuous controller in Eq.(3.38) with the controller in Eq.(3.69) and comparison of simulation results in Figures 3.6 and 3.7 indicate improvement in performance. Numerical simulations results of the algorithm using Acrobots parameters in Mahindrakar and Banavar [22] can be found in [134]. 3.5 Conclusion We proposed two swing-up algorithms for the Acrobot based on continuous and impulsive control inputs. The first algorithm is energy-based and is a modification of the algorithm 55 developed earlier for the Pendubot [2]. The algorithm for the Pendubot can be directly used for swing-up of the Acrobot but provides slow convergence; a modification of the algorithm, based on a set of impulsive inputs that increase the energy of the system, provides rapid swing-up. As compared to other general algorithms in the literature, our algorithm is not sensitive to control gains and does not impose constraints on the control gains or initial conditions. Numerical simulations were presented to show swing-up from an initial condition that corresponds to the equilibrium configuration with the lowest potential energy. This initial condition cannot be handled by existing algorithms in the literature. The idea of enlarging the set of admissible controls to include impulsive input is quite general and can be used to modify existing algorithms in the literature. This is demonstrated by the second swing-up algorithm presented in this chapter, which is a modification of the algorithms in [3] and [22]. The modification based on impulsive inputs imposes an additional constraint on the control gains but provides greater flexibility in choosing the continuous control inputs. This results in improved performance, which is demonstrated through numerical simulations. The numerical simulations also indicate that the additional constraint is trivially satisfied. 56 Chapter 4 Enlarging Region of Attraction for Underactuated Systems 4.1 Introduction The main algorithms for stabilization of underactuated systems were discussed in Section 1.3. Regardless of the algorithm, the system trajectories will not approach the desired equilibrium if the configuration of the system lies outside the region of attraction, RA . In this chapter, we propose a control algorithm which can enlarge the RA of a desired equilibrium for underactuated systems. This implies that it may be possible to stabilize an equilibrium from a configuration which lies outside the RA for an existing controller. Our algorithm uses impulsive inputs in conjunction with the continuous control inputs and its success depends on two sufficient conditions that need to be satisfied by the impulsive inputs. This chapter is organized as follows: In Section 4.2, we model the effect of impulsive inputs on the dynamics of underactuated systems. An algorithm to enlarge the RA of an equilibrium of the system is proposed in section 4.3. In Section 4.4, the algorithm is applied to three underactuated systems, namely, the Pendubot, the Acrobot, and the Rolling Acrobot, for all of which linear stabilizing controllers were designed. In Section 4.5, the algorithm is applied to the inverted pendulum on slope and the ball-beam system for which nonlinear stabilizing controllers were designed. Concluding remarks are provided in section 4.6. 57 4.2 4.2.1 Background Dynamic Model For an underactuated dynamical system, Lagrange’s equations have the form: M(q)¨ + N(q, q)q + G(q) = u q ˙ ˙ (4.1) where q ∈ Rn represents the vector of generalized coordinates, and M(q) ∈ Rn×n , N(q, q) ∈ ˙ Rn×n and G(q) ∈ Rn denote the positive definite inertia matrix, the matrix of centrifugal and Coriolis forces, and the vector of gravitational forces, respectively. u ∈ Rn is the vector of generalized forces and has the following form: u = [0, · · · , 0, u1, · · · , um ]T , m0 B : (ˆ, q + ) q ˙ A : (ˆ, q − ) q ˙ se ul p Im ine L ˙ V <0 V 0 for i = 1, 2, · · · , w, apply impulsive inputs uimp = [u1 , · · · imp , um ]T such that: V (q, q + ) < V (q, q − ) , ˙ ˙ ˙ V (q, q + ) < 0 ˙ (4.14) where hi (q) denotes the i-th constraint function imposed on the states of the system. 63 The effectiveness of the stabilization algorithm presented above will be illustrated through several numerical simulations in the following sections. 4.4 Case Studies: Linear Stabilizing Controllers 4.4.1 The Pendubot 4.4.1.1 Dynamics The Pendubot is shown in Figure 4.2. Assuming no friction in the joints, the equations of motion of the Pendubot are in the form of Eq.(4.1), where q = [θ2 θ1 ]T , τ1 ]T u = [0 (4.15) and M(q), N(q, q) and G(q), given by the expressions ˙   M(q) =  α2 α2 + α3 C 2 α2 + α3 C 2    , α1 + α2 + 2 α3 C 2    ˙ α5 C12 0 θ1     , G(q) = g  N(q, q) = α3 S2  ˙   ˙ ˙ ˙ α4 C1 + α5 C12 −(θ1 + θ2 ) −θ2 (4.16) In Eq.(4.16), Cj = cos θj , Sj = sin θj for j = 1, 2, C12 = cos(θ1 + θ2 ), and the positive constants αj , j = 1, 2, · · · , 5, have the following expressions 2 α1 = m1 d2 + m2 l1 + I1 , 1 α3 = m2 l1 d2 , α2 = m2 d2 + I2 2 α4 = m1 d1 + m2 l1 , 64 α5 = m2 d2 (4.17) g θ2 d1 θ1 d2 l1 l2 c.m . Y τ2 . c.m X τ1 Figure 4.2: The Pendubot/Acrobot is shown in an arbitrary configuration. For the Pendubot we have τ2 = 0 and for the Acrobot we have τ1 = 0. where, for k = 1, 2, mk and Ik are the mass and moment of inertia for the k-th link, and lk and dk are link lengths that are shown in Figure 4.2. The jump in the velocity of the passive coordinate θ2 due to the application of impulse in the θ1 coordinate can be obtained from Eqs.(4.9) and (4.16) as follows: ˙ + ˙− (θ2 − θ2 ) = − 4.4.1.2 α2 + α3 C 2 ˙ + ˙ − (θ1 − θ1 ) α2 (4.18) Continuous Controller ˙ ˙ To stabilize the system about its upright equilibrium (θ1 , θ1 , θ2 , θ2 ) = (π/2, 0, 0, 0), we use the pole-placement method for the linearized system. This controller renders the equilibrium locally asymptotically stable, and this can be established using the Lyapunov function: V = XT P X 65 (4.19) ˙ ˙ In Eq.(4.19), X = (θ1 , θ1 , θ2 , θ2 )T and P is the solution of the Lyapunov equation: (A − BK)T P + P (A − BK) + Q = 0 (4.20) where A and B describe the linearized system in state-space form, K is the vector of control gains, and Q is the identity matrix. In the next section, we use the numerical simulations to show that our algorithm can enlarge the RA of the upright equilibrium. 4.4.1.3 Numerical Simulations Consider the Pendubot with the following parameters: m1 = 1, l1 = 2, d1 = 1, I1 = 0.333 m2 = 1, l2 = 2, d2 = 1, I2 = 0.333 (4.21) The poles of the linearized system were placed at λdes = (−0.5, −1, −1.2, −1.5). ˙ ˙ For the first simulation, the initial condition was assumed to be (θ10 , θ10 , θ20 , θ20 ) = (1.484, −0.45, 0.087, 0.15), where the units are in rad and rad/s. This configuration, which ˙ corresponds to V0 = 15.037 and V0 = 2.7746, lies outside RA of the upright equilibrium. This can be verified through numerical simulations over a short interval of time for which the trajectories go unbounded. Figure 4.3 shows simulation results obtained using our algorithm. Impulsive inputs are applied at seven different time instants: at t = 0, and six more times in the interval t ∈ [0.1, 2.25]. It can be observed that the impulsive inputs result in jumps ˙ ˙ in the joint velocities, and jumps in V to maintain V < 0. Since the initial impulse is 66 80 60 50 a b 0 0 -1 θ2 (deg) c ˙ θ1 (rad/s) ˙ θ2 (rad/s) 1 0 d e 20 10 0 5 0 -5 0 -50 θ1 (deg) V f ˙ V g 0 0.5 1 1.5 2 2.5 time (sec) 3 3.5 τ1 (N.m) 4 4.5 5 Figure 4.3: Simulation results showing joint angles and joint velocities of the Pendubot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - First simulation large in magnitude compared to other impulsive inputs, the negative jump in V can be observed for the first impulse only. The stabilization of the upright equilibrium establishes the effectiveness of our algorithm. The effect of the impulsive input applied at t = 0 is illustrated with the help of Figure ˙ ˙ ˙ 4.4. This figure plots the curves V = V0 , V = 0, and the impulse line in θ1 -θ2 plane at the initial configuration. The impulsive input satisfies the conditions in Eq.(4.12) and causes the configuration to jump from A : (q0 , q0 ) to B : (q0 , q0 ) along the impulse line. At B, ˙− ˙+ ˙ V (q0 , q0 ) < 0 and V (q0 , q0 ) < V0 . A numerical search was used to find the impulsive input ˙+ ˙+ ˙ which results in a maximum decrease in V while ensuring V < 0. ˙ ˙ For the second simulation, the initial condition was assumed to be (θ10 , θ10 , θ20 , θ20 ) 67 6 4 B : (q0 , q0 ) ˙+ ˙ V <0 ˙ θ2 (rad/s) 2 ˙ V =0 ˙ V >0 0 A : (q0 , q0 ) ˙− -2 V = V lse pu Im 0 -4 ne Li -6 -3 -2 -1 0 ˙ θ1 (rad/s) 1 2 3 ˙ ˙ ˙ Figure 4.4: The curves V = V0 , V = 0, and the impulse line are shown in the θ1 -θ2 plane at the initial configuration of the Pendubot for first simulation, where θ10 = 1.484 and θ20 = 0.087. The impulsive input at t = 0 in Figure 4.3 results in a jump in the configuration ˙ from A : (q0 , q0 ) to B : (q0 , q0 ), where V (q0 , q0 ) < 0 and V (q0 , q0 ) < V0 . ˙− ˙+ ˙+ ˙+ = (1.745, 0, 0.174, 0), where the units are in rad and rad/s. This configuration, which cor˙ responds to V0 = 32.489 and V0 = 4.433, lies outside RA of the upright equilibrium. This can be verified through numerical simulations over a short interval of time for which the trajectories go unbounded. Figure 4.5 shows simulation results obtained using our algorithm. Impulsive inputs are applied at four different time instants: at t = 0, and three more times in the interval t ∈ [1.6, 2.25]. It can be observed that the impulsive inputs result in jumps ˙ ˙ in the joint velocities, and jumps in V to maintain V < 0. 68 a 0 -40 2 1 0 θ1 (deg) b 120 100 80 θ2 (deg) ˙ θ1 (rad/s) c 0 -2 d 40 20 0 e 0 -10 60 f ˙ θ2 (rad/s) V ˙ V g τ1 (N.m) 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 time (sec) Figure 4.5: Simulation results showing joint angles and joint velocities of the Pendubot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - Second simulation 4.4.2 The Acrobot 4.4.2.1 Dynamics and Control The Acrobot is shown in Figure 4.2. Assuming no friction in the joints, the equations of motion has the form of Eq.(4.1), where q = [θ1 θ2 ]T , u = [0 τ2 ]T and M(q), N(q, q), and G(q), given by the expressions ˙ 69 (4.22)    α1 + α2 + 2 α3 C 2 α2 + α3 C 2  M(q) =  , α2 + α3 C 2 α2     ˙2 −(θ1 + θ2 ) ˙ ˙  −θ   α4 C1 + α5 C12  G(q) = g  N(q, q) = α3 S2  ˙ ,  ˙1 θ 0 α5 C12 (4.23) The parameters αj , j = 1, 2, · · · , 5 in Eq.(4.23), were defined in Eq.(4.17). The jump in the velocity of the passive coordinate θ1 due to the application of impulse in the θ2 coordinate can be obtained from Eq.(4.9) and Eq.(4.23) as follows: ˙+ ˙ − (θ1 − θ1 ) = − α2 + α3 C 2 ˙ + ˙− (θ2 − θ2 ) α1 + α2 + 2 α3 C 2 (4.24) Using the pole-placement method, a stabilizing controller was designed to render the upright equilibrium locally asymptotically stable. The procedure followed is similar to that used for the Pendubot in Section 4.4.1.2. 4.4.2.2 Numerical Simulations Consider the following parameters for the Acrobot: m1 = 1, l1 = 1, d1 = 0.5, I1 = 0.083 m2 = 1, l2 = 2, d2 = 1, I2 = 0.333 (4.25) The poles of the linearized system were placed at λdes = (−2, −3, −4, −5). ˙ ˙ For the first simulation, the initial conditions are assumed to be (θ10 , θ10 , θ20 , θ20 ) = 70 100 θ1 (deg) a 50 100 b 50 0 0 c -4 5 0 -5 d θ2 (deg) ˙ θ1 (rad/s) ˙ θ2 (rad/s) 10 0 e 0 V f -40 50 0 ˙ V g 0 τ2 (N.m) 0.5 1 1.5 time (sec) 2 2.5 3 Figure 4.6: Simulation results showing joint angles and joint velocities of the Acrobot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - First Simulation ˙ (1.658, 0, 0.087, 0). This configuration corresponds to V0 = 24.56 and V0 = 1.04. Using numerical simulations, it was ascertained to be outside the RA of the upright equilibrium. The simulation results obtained using our algorithm is shown in Figure 4.6. Four impulsive inputs were applied which includes one impulse at t = 0 and three other impulses in the ˙ ˙ interval t ∈ [0.4, 0.9]. The jumps in the joint velocities, and jumps in V to maintain V < 0 can be observed in the figure. Also, the first two negative jumps in V which are large in magnitude can be seen from the plot. The stabilization of the equilibrium from outside the RA establishes the efficacy of our algorithm. ˙ ˙ For the second simulation, the initial conditions were assumed to be at (θ10 , θ10 , θ20 , θ20 ) = ˙ (1.658, 0.08, −0.087, 0.1). This configuration corresponds to V0 = 17.7 and V0 = 0.82 which is outside the RA of the upright equilibrium as well. The simulation results are shown in 71 100 a 0 2 0 -4 5 0 θ1 (deg) b 50 150 θ2 (deg) c ˙ θ1 (rad/s) ˙ θ2 (rad/s) d 20 10 e 0 0 f -40 50 0 -50 V ˙ V g 0 τ2 (N.m) 0.5 1 1.5 time (sec) 2 2.5 3 Figure 4.7: Simulation results showing joint angles and joint velocities of the Acrobot and the Lyapunov function, its time derivative, and the control input for an initial configuration outside RA - Second Simulation Figure 4.7. Four impulsive inputs including the initial impulse were required for this case and the system is ultimately stabilized using the control algorithm. 4.4.3 Rolling Acrobot 4.4.3.1 Dynamics and Control A rolling acrobot, a modified version of the Acrobot, is a two-link planar robot with a passive rolling shoulder joint and an active elbow joint. Figure 4.8 shows a schematic of the robot. Assuming no friction in the joints, the equations of motion take the form of Eq.(4.1) where q = [θ1 θ2 ]T , u = [0 τ2 ]T 72 (4.26) d2 l2 X Y τ2 θ2 g l1 θ1 d1 Figure 4.8: A schematic of the rolling acrobot. The link angles θ1 and θ2 are measured counter-clockwise. and M(q), N(q, q), and G(q), given by the expressions ˙    γ1 + γ2 − 2γ3 C1 − 2γ4 C12 γ2 − γ4 C12  M(q) =  , γ2 − γ4 C12 γ2   ˙1 + 2γ4 S12 θ2 γ4 S12 θ2 ˙ ˙   (γ3 S1 + γ4 S12 )θ N(q, q) =  ˙ , 0 0    γ5 S1 + γ6 S12  G(θ) = g   γ6 S12 (4.27) In Eq.(4.27), the positive constants γi , i = 1, 2, · · · , 6 have the following expressions 2 γ1 =I1 + m1 d2 − 2m1 l1 d1 + (2m1 + m2 )l1 , 1 γ3 =m1 l1 (l1 − d1 ), γ4 = m2 l1 (l2 − d2 ), γ2 = I2 + m2 (l2 − d2 )2 γ5 = m1 (l1 − d1 ), 73 γ6 = m2 (l2 − d2 ) (4.28) The velocity jump in θ1 coordinate due to the application of impulse in θ2 coordinate can be obtained from Eq.(4.9) and Eq.(4.27), namely, ˙ + ˙− (θ1 − θ1 ) = − γ2 − γ4 C12 ˙ + ˙− (θ2 − θ2 ) γ1 + γ2 − 2γ3 C1 − 2γ4 C12 (4.29) The same pole placement controller as in Section 4.4.1.2 is designed for the rolling acrobot ˙ ˙ to render the upright equilibrium (θ1 , θ1 , θ2 , θ2 ) = (0, 0, π, 0) locally asymptotically stable. 4.4.3.2 Numerical Simulations The kinematic and dynamic parameters of the rolling acrobot is chosen as: m1 = 1, l1 = 1, d1 = 0.5, I1 = 0.083 m2 = 2, l2 = 1, d2 = .05, I2 = 0.167 (4.30) For simulation results, we consider two cases of the rolling acrobot; In the first case, we assume there is no limit for θ1 and the system can roll on its arc foot without any constraint. In the second case, the arc foot is limited and the angle θ1 is constrained to be inside a prescribed range during the stabilization process. First Case - Unbounded Foot Arc Angle The poles of the linearized system are placed at λdes = (−0.5, −1, −1.5, −2) using pole ˙ ˙ placement control. For the first simulation, the initial conditions are set to (θ10 , θ10 , θ20 , θ20 ) = ˙ (0.436, 0, 3.491, 0). This configuration corresponds to V0 = 12.5 and V0 = 5.33 and is outside the RA of the upright equilibrium which can be verified through a numerical simulation. The simulation results are shown in Figure 4.9. The applied impulsive inputs can be recognized 74 40 20 0 a 200 θ1 (deg) 150 b 1 0 c -1 0 d -2 -4 θ2 (deg) 5 0 5 0 -5 ˙ θ1 (rad/s) ˙ θ2 (rad/s) e V f 0 -40 g 0 ˙ V 1 2 3 4 time (sec) 5 6 τ2 (N.m) 7 8 Figure 4.9: Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with unbounded leg angle originally outside RA - First simulation ˙ by the picks in the control input subplot. All impulsive inputs result in a negative V and a drop in V as can be seen in the figure. Using these impulsive inputs, the upright equilibrium is stabilized and the configuration which was originally outside the RA , is now attracted to the desired equilibrium. For the second simulation, we choose the initial conditions with initial velocities to be at ˙ ˙ (θ10 , θ10 , θ20 , θ20 ) = (0.087, 0.5, 3.229, 1). This configuration corresponds to V0 = 23.36 and ˙ V0 = −1.195 which is outside the RA of the upright equilibrium as well. The simulation results are shown in Figure 4.10. Three impulsive inputs including the initial impulse were required for this case and the system is ultimately stabilized using the control algorithm. 75 60 20 -20 a 200 150 100 b 2 0 c θ1 (deg) θ2 (deg) ˙ θ1 (rad/s) 0 d ˙ θ2 (rad/s) -4 10 e 0 V f 0 -20 20 0 -20 ˙ V g 0 1 2 3 4 time (sec) 5 6 τ2 (N.m) 7 8 Figure 4.10: Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with unbounded leg angle originally outside RA - Second simulation Second Case - Bounded Foot Arc Angle The leg angle θ1 is constrained as follows: h1 (q) = |θ1 | − η1 ≤ 0 (4.31) where η1 = 15◦ . The poles of the closed-loop linearized system are placed at λdes = (−1, −1.5, −2, −2.5). The modified algorithm mentioned in Remark 4.3.3 is followed to stabilize the upright equilibrium while ensuring the constraint in Eq.(4.31) is satisfied. ˙ ˙ For the first simulation, the initial conditions are assumed to be at (θ10 , θ10 , θ20 , θ20 ) = ˙ (0.087, 0, 3.229, 0). This configuration corresponds to V0 = 1.218 and V0 = 0.366. Using 76 20 a 10 0 200 180 160 b 0.4 0 c 0 -1 -2 θ1 (deg) θ2 (deg) ˙ θ1 (rad/s) d ˙ θ2 (rad/s) 1 e 0 0 -2 V f 0 -5 -10 ˙ V τ2 (N.m) g 0 0.5 1 1.5 2 2.5 time (sec) 3 3.5 4 4.5 5 Figure 4.11: Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with bounded leg angle originally outside RA - First simulation the pole placement controller, it can be verified through a numerical simulation that the constraint in Eq.(4.31) is not satisfied. Figure 4.11 shows the simulation results for our algorithm where impulsive inputs were used to stabilize the upright equilibrium while satisfying the state constraints. Two impulsive inputs were applied for this simulation; the first ˙ impulse, applied at the initial time, is due to V0 > 0 while the second impulsive input is applied when θ1 = η1 to keep θ1 inside its bound. Both impulsive inputs result in negative ˙ V as well as a negative jump in V . Using these impulsive inputs, the upright equilibrium is stabilized without exceeding the bound on θ1 . For the second simulation, we choose the initial conditions with initial velocities to be ˙ ˙ at (θ10 , θ10 , θ20 , θ20 ) = (0, −0.1, 2.967, −0.25). This configuration corresponds to V0 = 1.298 77 0 -10 -20 200 180 160 0 -0.5 θ1 (deg) a b θ2 (deg) c ˙ θ1 (rad/s) ˙ θ2 (rad/s) 1 0 d e 1 0 0 V f ˙ V -4 10 5 0 g 0 τ2 (N.m) 0.5 1 1.5 2 2.5 time (sec) 3 3.5 4 4.5 5 Figure 4.12: Simulation results including the angular positions, angular velocities, Lyapunov function and its time derivative and the control input for a configuration of rolling acrobot with bounded leg angle originally outside RA - Second simulation ˙ and V0 = 0.094 and it is verified that the bound on θ1 will be exceeded using the pole placement controller. The simulation results for our algorithm are shown in Figure 4.12. ˙ As the first simulation, the initial impulsive input is applied since V0 > 0 and the second impulsive input is intended to keep θ1 inside its limit. 4.5 Extension for Nonlinear Controllers The case studies studied in Section 4.4 consider the systems with linear stabilizing controllers which is only valid in a small neighborhood of the desired equilibrium. For such systems, ˙ the derivative of the corresponding Lyapunov function V is more likely to get positive when the system configuration lies outside the RA since the controller is designed based on the 78 ˙ linearized system. For system with nonlinear stabilizing controllers, however, the V is designed to be less than or equal to zero1 for the nonlinear system and therefore, the condition ˙ to apply the impulsive inputs which is V ≥ 0 as discussed in Section 4.3, might never be satisfied. In order to use the benefits of impulsive inputs in enlarging the RA for systems with nonlinear controllers, we modify the stabilization algorithm proposed in Section 4.3 as follows: 1. Initialization: Find the Lyapunov function V (q, q) for the underactuated system in ˙ Eq.(4.1) with stabilizing nonlinear controller u = [u1 , · · · , um ]T . ˆ ˙ 2. Impulsive Inputs: If V ≥ 0 or the system configuration lies outside the set Ω, apply imp imp impulsive inputs uimp = [u1 , · · · , um ]T such that: V (q, q + ) < V (q, q − ) , ˙ ˙ ˙ V (q, q + ) < 0 ˙ (4.32) where Ω is an estimate of the RA calculated based on the Lyapunov function of the closed-loop system. Remark 4.5.1. Compared to the algorithm in Section 4.3, we need an off-line calculation of the set Ω for the algorithm above. Although it might be hard to find a set which estimates the RA of the desired equilibrium closely, any conservative estimate of the RA will be adequate for the proposed algorithm to be applied. Remark 4.5.2. When there exist constraints on the states of the system, following the discussion in Remark 4.3.3, the additional impulsive inputs are applied to ensure the constraints ˙ that ensuring V ≤ 0 does not necessarily mean that the desired equilibrium is stable. This is because the states of the system might get unbounded using the designed controller or the system approaches another equilibrium of the system. 1 Note 79 are satisfied. In the rest of this section, we consider two well-known underactuated systems in the literature with nonlinear stabilizing controllers and apply the stabilization algorithm proposed above to enlarge the RA of the desired equilibrium. 4.5.1 Inverted Pendulum on an Inclined Plane 4.5.1.1 Dynamics Consider the inverted pendulum system on an inclined plane shown in Figure 4.13. Assuming no friction in the joints, the equations of motion has the form of Eq.(4.1) where q = [φ s]T , u = [0 F ]T (4.33) and M(q), N(q, q), and G(q), given by the expressions [25] ˙   β cos(φ − ψ)  , β cos(φ − ψ) γ     0 0   −β sin φ   N(q, q) =  ˙   , G(q) = g  ˙ 0 −γ sin ψ −β sin(φ − ψ)φ  M(q) =  α (4.34) where γ = M + m , α = ml2 , β = ml (4.35) and φ and s are the generalized coordinates as shown in Figure 4.13, m and l are the point mass and the length of the pendulum, M is the mass of the cart and ψ is the slope angle. 80 r φ l g F ψ s Figure 4.13: Cart-pendulum on a slope The jump in the velocity of the passive coordinate φ due to the application of impulsive force in active coordinate s, can be obtained from Eqs.(4.9) and (4.34), namely, ˙ ˙ (φ + − φ − ) = − 4.5.1.2 β cos(φ − ψ) ( s+ − s− ) ˙ ˙ α (4.36) Stabilizing Control Design ˙ The nonlinear controller to stabilize the zero equilibrium at (φ, φ, s, s) = (0, 0, 0, 0) is ob˙ tained from [25] using the Lagrangian Controlled method: F = −γg sin ψ + ˙ ˙ κβ[α sin(φ − ψ)φ2 + D cos(φ − ψ) sin φ] − B ∂H + B c γ y ∂s α− β2 2 γ (κ + 1) cos (φ − ψ) (4.37) where αγ − β 2 cos2 (φ − ψ) , D = −mgl ργ ρ−1 β ǫDγ 2 2 )[sin(φ − ψ) + sin ψ] y , y = s + (κ + H = mgl cos(φ) + 2 γ ρ 2β κ = 1/σ, B= (4.38) and σ, ρ < 0, ǫ > 0, and c > 0 are some scalar constants. The asymptotic stability of the 81 equilibrium is established using the following Lyapunov function which is obtained from the controlled energy of the system [25]: κ 1 ˙ ˙ ˙ V =mgl − αφ2 − β cos(φ − ψ)[s + β cos(φ − ψ)φ]φ 2 γ 1 κ κ ˙ ˙ − γ[s + β cos(φ − ψ)φ]2 + β 2 cos2 (φ − ψ)φ2 ˙ 2 γ 2γ 1 β ˙ − (ρ − 1)γ[s + (κ + 1) cos(φ − ψ)φ]2 − H ˙ 2 γ (4.39) An estimate of the RA is obtained as: ˙ Ωc = z = (φ, φ, s, s) ∈ T Q | V (z) ≤ c ˙ ¯ ¯ (4.40) where T Q is a bounded set, and c is a constant scalar. ¯ 4.5.1.3 Numerical Simulations The parameters for the inverted pendulum system shown in Figure 4.13 are chosen as: m = 0.14 Kg, M = 0.44 Kg, L = 0.215 m, ψ = π/9 rad (4.41) and the control parameters are chosen as: κ = 20, ρ = −0.02, ǫ = 0.00001, c = 0.015 (4.42) For the parameters in Eqs.(4.41) and (4.42), the estimate of the RA in Eq.(4.40) is determined by c = 0.26. We also assume that the pendulum angle is constrained to be inside a pre¯ 82 60 0 -60 4 0 a φ (deg) b ˙ φ (rad/s) 5 0 c 4 0 d 0.5 e 0 0 -2 s (m) s (m/s) ˙ V ˙ V f g 60 F (N) 0 0 5 10 15 20 25 30 time (sec) Figure 4.14: Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the inverted pendulum on a slope originally outside the RA - First simulation defined range, i.e., h1 (q) = |φ| − φmax ≤ 0 where φmax = 1.222 rad = 70 deg. ˙ For the first simulation, we choose the initial conditions of the system as (φ0 , φ0 , s0 , s0 ) = ˙ ˙ (π/6, 4, 3, 3). This configuration corresponds to V0 = 0.564 and V0 = −2.669 and lies outside the set Ωc in Eq.(4.40). From numerical simulations for a short period of time, it can be ¯ confirmed that the sates of the system get unbounded using the stabilizing controller in Eq.(4.37) and therefore, this configuration is outside the RA of the desired equilibrium. Figure 4.14 shows the simulation results for the proposed stabilization algorithm in Section 4.5. The peak in the control input subplot shows the initial impulsive input applied at t = 0 sec due to the fact that the initial configuration is outside the set Ωc . This impulsive input ¯ ˙ causes a large negative jump in V with V < 0 as required in Eq.(4.32). Using the proposed algorithm, the desired equilibrium at the origin is stabilized from a configuration which was 83 80 40 0 4 0 a b ˙ φ (rad/s) c s (m) d e 0 -2 4 0 -4 0.3 φ (deg) s (m/s) ˙ V 0 0 ˙ V f g -0.6 40 F (N) 0 0 5 time (sec) 10 15 Figure 4.15: Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the inverted pendulum on a slope originally outside the RA - Second simulation originally outside the RA . For the second simulation, the pendulum is assumed to be farther from upright equilib˙ rium position compared to the first simulation by choosing (φ0 , φ0 , s0 , s0 ) = (π/3, 4, −2, −3.5). ˙ The corresponding values of the Lyapunov function and its time derivative are calculated ˙ as V0 = 0.251 and V0 = −0.528. Although this configuration lies inside Ωc , the numerical ¯ simulation using the stabilizing controller in Eq.(4.37) shows that the angular position of the pendulum exceeds its limit. Simulation results for our algorithm is shown in Figure 4.15. The only impulsive input is applied a little bit after the starting time when the pendulum angle reaches its limit, φ = 70 deg. The impulsive input kicks back the pendulum inside the allowed region and cause a small drop in V which is enough to stabilize the desired equilibrium without exceeding the limit for φ. 84 4.5.2 The Ball and Beam System 4.5.2.1 Dynamics The ball and beam system is shown in Figure 4.16. The equations of motion is in the form of Eq.(4.1) where q = [r θ]T , u = [0 τ ]T (4.43) and M(q), N(q, q), and G(q), given by the expressions ˙   M(q) =    N(q, q) =  ˙ Jb R2 +M 0 Mr 2 + J + Jb 0  ˙ −Mr θ  , ˙ 2Mr θ 0 0   ,    M sin θ  G(q) = g   Mr cos θ (4.44) In Eq.(4.44), r and θ are the generalized coordinates as shown in Figure 4.16, M, Jb , and R are the mass, mass moment of inertial and radius of the ball, respectively, and J is the mass moment of inertia of the beam. The length of the beam is assumed to be equal to 2L. Using Eqs.(4.9) and (4.44), it can be shown that application of an impulsive input in τ r θ Figure 4.16: Ball and beam system 85 active coordinate θ does not cause any velocity jump in the passive coordinate r, i.e., r+ = r− ˙ ˙ 4.5.2.2 (4.45) Stabilizing Control Design A nonlinear stabilizing controller for the ball and beam system is given using the IDA-PBC method in [28] as follows: τ =√ r 2(L2 + r2) −p2 1 + gr cos θ − kp −g L2 + r 2 + √ 2p1 p2 + √ 1 L2 + r2 p2 2 L2 + r 2 1 r θ − √ arcsinh( ) 2 L 2 kv p1 − p2 2(L2 + r 2 ) sin θ + 2 L + r2 L2 2 + r2 (4.46) where p1 = r , ˙ ˙ p2 = (L2 + r 2 )θ (4.47) and kp and kv are the constant control gains. The control input in Eq.(4.46) was obtained assuming that M =1, J = L2 , Jb = 0 (4.48) ˙ Asymptotic stability of the desired equilibrium at (r, r, θ, θ) = (0, 0, 0, 0) is established ˙ using the following Lyapunov function which is the Hamiltonian of the desired closed-loop 86 system [28]: 1 −1 V = P T Md P + Vd 2 (4.49) where    p1  P = , p2   Md =   √ √ 2 L2 +r2 1 1 2(L2 + r 2 ) kp 1 r 2 θ − √ arcsinh( ) Vd = g(1 − cosθ) + 2 L 2     (4.50) The RA of the equilibrium such that the ball remains on the beam for all the time, i.e., h1 (q) = |r(t)| − L ≤ 0 for all t ≥ 0, is estimated as [28]: ¯ Ω= (q, p) | kp 1 1 kp g r 1 T −1 θ − √ arcsinh( ) < P Md P + g(1 − cosθ) + 2 2 L 8 2kp + g 2 (4.51) which gives a conservative estimate based on the desired Hamiltonian of the system. Remark 4.5.3. Although the coordinate r is constrained for the ball and beam system, application of an impulsive input is not helpful in kicking back the constrained state inside the allowed region when it gets to its limit as proposed in Remark 4.3.3. This is because the impulsive input does not cause any velocity jump in r as derived in Eq.(4.45). The only way ˙ to ensure a configuration of the system approaches the origin without exceeding the limit on ¯ r is to check whether the configuration lies inside the set Ω in Eq.(4.51). 87 4.5.2.3 Numerical Simulations For simulations, the following parameters for the system and the controller are chosen: L = 1, kp = 1, kv = 2 (4.52) ˙ For the first simulation, we choose the initial conditions as (r0 , r0 , θ0 , θ0 ) = (0.8, 0, 0.175, −1). ˙ ˙ ¯ This configuration corresponds to V0 = 1.113 and V0 = −2.439 and lies outside the set Ω in Eq.(4.51). Using the stabilizing controller in Eq.(4.46), the ball will exceed its limit and gets off the beam and therefore, this configuration is outside the RA of the desired equilibrium. Figure 4.17 shows the simulation results for the proposed stabilization algorithm. The initial impulsive input shown by the peak in the control input subplot, is due to the fact that the 1 0.5 0 a 0 b -0.4 10 0 -10 c r (m/s) ˙ 0 d -0.4 ˙ θ (rad/s) r (m) θ (deg) 1 e 0 0 -2 f V ˙ V 50 g 0 0 τ (N.m) 5 time (sec) 10 15 Figure 4.17: Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the ball and beam system originally outside the RA - First simulation 88 ¯ initial configuration is outside the set Ω. This impulsive input causes a large negative jump ˙ in V with V < 0 as required in Eq.(4.32). Using the proposed algorithm, the desired equilibrium at the origin is stabilized without exceeding the limit for r. ˙ For the second simulation, the initial conditions were chosen as (r0 , r0 , θ0 , θ0 ) = (−0.6, ˙ −0.5, 0, 0.5). The corresponding values of the Lyapunov function and its time derivative are ˙ ¯ calculated as V0 = 0.689 and V0 = −1.98. This configuration is also outside the set Ω and it can be shown with the help of numerical simulations that the ball exceeds its limit using the controller in Eq.(4.46). Simulation results for our algorithm is shown in Figure 4.18. The only impulsive input is applied at the initial time causing a negative jump in V . The continuous controller is applied thereafter. The initial impulse in enough to enlarge the RA of the desired equilibrium to include the initial configuration of the system. 0 -0.5 -1 0.5 0 -0.5 10 0 -10 0.5 0 -0.5 0.5 0 0 a r (m) r (m/s) ˙ b c θ (deg) d ˙ θ (rad/s) e V ˙ V -2 f 0 -40 g 0 τ (N.m) 4 time (sec) 8 12 Figure 4.18: Simulation results including the positions and velocities of the system coordinates, Lyapunov function and its time derivative and the control input for a configuration of the ball and beam system originally outside the RA - Second simulation 89 4.6 Conclusion An impulsive control algorithm is proposed to enlarge the region of attraction of stabilized equilibria for underactuated systems. The algorithm requires impulsive inputs to be applied in conjunction with the stabilizing controller. The impulsive inputs are applied when the derivative of the Lyapunov function is non-negative and two conditions are satisfied. These conditions require that the impulsive inputs result in a negative jump in the Lyapunov function and a jump to a configuration where the derivative of the Lyapunov function is strictly negative. Although there is no guarantee that these conditions can be satisfied at any given configuration, numerical simulations show that it is possible to stabilize an equilibrium point from certain configurations that lie outside its region of attraction. The effectiveness of the algorithm is demonstrated using several examples of underactuated mechanical systems with linear and nonlinear continuous controllers. 90 Chapter 5 Safe Fall Control for Humanoid Robots 5.1 Introduction Impulsive control inputs are a natural candidate for emergency control when the available time to influence the dynamic system is limited. Furthermore, due to a very short “time of application” of the control, the system configuration does not evolve much during the application of the control. A direct outcome of this property is that the control can be applied even when the system is close to the boundary of its configuration space. For example, a robot can accept control command even when it is close to its joint limits. In this chapter, we propose a control strategy for changing the default fall direction of the robot so that it avoids contact with surrounding objects or people as a mean of reducing damage to others. A six-dof robot with three actuated joints is considered here to better understand the impulsive algorithm before applying to higher dof humanoids. Note that we consider only those situations in which a fall is caused by an external factor; in particular the sensors and the motor power must remain intact such that the robot can execute a prescribed control strategy. We should clarify that a fall controller is not a balance controller. A fall controller complements, and does not replace, a balance controller. Further, a fall controller is not a 91 push-recovery controller. A push-recovery controller is essentially a balance controller, which specifically deals with external disturbances of larger magnitude when the robot must take a step in order to regain balance [135, 136, 137]. The fall controller is activated only when the default balance controller or the push recovery controller has failed to stabilize the robot. The rest of this chapter is organized as follows: In Section 5.2, dynamic model of the six-dof robot is described. we derive the impulsive effects on the system velocities in Section 5.3. The safe fall algorithm is presented in Section 5.4 which is followed by the numerical simulations in Section 5.5. Concluding remarks are provided in Section 5.6. 5.2 5.2.1 Dynamic Model Equations of Motion We considered a three-dimensional seven-dof robot as shown in Figure 5.1. The degrees of freedom of the robot are listed in Table 5.1. Assuming that the knee joint is locked, i.e., θ1 ≡ 0, we have a six-dof robot with three active joints. Using Lagrange’s formulas, the equations of motion can be derived as follows: M(q)¨ + N(q, q)q + G(q) = T q ˙ ˙ (5.1) q = [α, β, γ, θ2 , θ3 , θ4 ]T , T =[0, 0, 0, τ2 , τ3 , τ4 ]T (5.2) where 92 z4 y4 x4 z3 y3 x3 z 2 Z y2 x2 z1 X Inertial Frame Y y1 Front View Side View x1 Figure 5.1: Seven-dof three-dimensional humanoid robot and M(q) ∈ R6×6 , N(q, q) ∈ R6×6 and G(q) ∈ R6 represent the mass matrix, the matrix ˙ of centrifugal and Coriolis forces and the vector of forces due to gravity, respectively. Full expressions of the matrices M(q), N(q, q) and G(q) are lengthy and are not provided for the ˙ sake simplicity. Table 5.1: Degrees of freedom for the humanoid robot in Figure 5.1 Link No. of dof’s Coordinates Actuation 1 3 α, β, γ: z-x-z Euler angles passive 2 1 θ1 : Rotation about x2 active θ2 : Rotation about x2 3 2 active θ3 : Rotation about z3 4 1 θ4 : Rotation about y4 active 93 5.2.2 Lean Line To better understand the orientation of the robot during fall, we consider the orientation of the line connecting the foot contact point with the ground to the center of mass of the whole system which we call the Lean Line. The vector representing the lean line can be calculated from: 4 n=1 mn rn 4 n=1 mn rcm = (5.3) where mn is the mass of the nth link and rn is the position vector of the center of mass of the nth link in the inertial frame. The orientation of the lean line in the inertial frame can be determined through the angles φ and ψ shown in Figure 5.2 where: rcm = |rcm | [sin(φ) cos(ψ) i + sin(φ) sin(ψ) j + cos(φ) k] Z φ rcm Y ψ X Figure 5.2: Orientation of the lean line in the inertial frame 94 (5.4) 5.3 Impulsive Effects Consider a set of impulsive inputs inputs applied in active coordinates θ2 , θ3 and θ4 at time t = ti . We integrate the equations of motion in Eq.(5.1) over the short period of time, ∆t when the impulses are applied: ti +∆t M(q)¨ dt + q ti ti +∆t N(q, q) q dt + ˙ ˙ ti ti +∆t G(q) dt = ti ti +∆t T dt (5.5) ti Since ∆t is very small and the configuration of the robot does not change over this short interval of time, Eq.(5.5) can be simplified as: M(q)(q + − q − ) = I ˙ ˙ (5.6) where q − and q + are the vectors of angular velocities immediately before and after the ˙ ˙ application of impulses, respectively, and I represents the vector of impulses, i.e.: I = [0, 0, 0, I2 , I3 , I4 ]T (5.7) The first three equations in Eq.(5.6) which correspond to the passive coordinates of the first link are written as: M (q)(q + − q − ) = 0 ˙ ˙ (5.8) where M (q) ∈ R3×6 is obtained from partitioning the mass matrix M(q). Eq.(5.8) can be 95 rewritten as:     + − θ− + − α− ˙ ˙ ˙ ˙  θ2 α 2      ˙+ ˙−   ˙+  −  + M2 (q) θ − θ  = 0 ˙ M1 (q) β − β 3  3       ˙+ ˙− θ4 − θ4 γ+ − γ− ˙ ˙ (5.9) where M1 (q) ∈ R3×3 and M2 (q) ∈ R3×3 are obtained from: M (q) = M1 (q) M2 (q) (5.10) Equation (5.9) can be used to find the relationships between the system velocities immediately before and after the application of impulse. If values of three velocities after the application of impulsive inputs are known, we can solve for the other three velocities after the application of impulsive inputs from Eq.(5.9). As the first case, assume that the desired values for the active coordinate velocities after the ˙+ ˙ ˙+ ˙ ˙+ ˙ application of impulsive inputs are specified as θ2 = θ2des , θ3 = θ3des and θ4 = θ4des , then from Eq.(5.9) we have:       − − + ˙ ˙ ˙ ˙ θ2des − θ2  α  α        ˙  ˙ +  ˙ − −1 ˙−  β  = β  − M1 (q) M2 (q) θ3des − θ  3            − + ˙ ˙− θ4des − θ4 γ ˙ γ ˙ (5.11) where M1 is invertible since it is a diagonal block of the positive definite matrix M(q). As the second case, assume that the desired values of the passive coordinate velocities after ˙ ˙ the application of impulsive inputs are specified as α+ = αdes , β + = βdes and γ + = γdes . ˙ ˙ ˙ ˙ 96 Then from Eq.(5.9) we can write:       − ˙− ˙+ ˙ ˙ αdes − α  θ2  θ2           +  − # ˙ ˙ ˙ ˙ θ  = θ  − M2 (q) M1 (q)  βdes − β −     3  3       − ˙− ˙+ θ4 γdes − γ ˙ ˙ θ4 (5.12) # where M2 (q) denotes the Pseudo-inverse of M2 and it is used since the matrix M2 might be singular. ˙+ ˙+ ˙+ Remark 5.3.1. If the matrix M2 is singular, the velocities [θ2 , θ3 , θ4 ] obtained from Eq.(5.12) ˙ do not exactly correspond to the desired values for the passive coordinates [αdes , βdes , γdes ]. ˙ ˙ The resulted velocities, however, will give the closest solution in terms of the least Euclidean ˙ norm to the active coordinate velocities which correspond to [αdes , βdes , γdes ]. The veloci˙ ˙ ties for the active coordinate obtained from Eq.(5.12) can be substituted into Eq.(5.11) to yield the corresponding velocities of the passive coordinates. To find the applied impulses which result in some desired velocity jumps, we consider the last three equations in Eq.(5.6), namely, ¯ ¯ M (q)(q + − q − ) = I ˙ ˙ (5.13) ¯ ¯ where M (q) ∈ R3×6 and I = [I2 , I3 , I4 ]T . Having the velocities after the application of the impulsive inputs from Eqs.(5.11) and (5.12), the vector of the applied impulses can be calculated using Eq.(5.13). Remark 5.3.2. The ideal impulses obtained from Eq.(5.13) can not be implemented in the standard actuators. However, we can use the high-gain continuous control inputs to estimate 97 these impulses as discussed earlier in Chapter 2. 5.4 Safe Fall Algorithm We propose a safe fall algorithm for the system in Eq.(5.1) using two sets of impulsive inputs. The first set of impulsive inputs is applied in order to reach some desired velocities for the passive coordinates and the second set is intended to brake the active coordinates which are held fixed thereafter. The robot undergoes a free fall for the time period between the impulsive sets. A schematic of the proposed algorithm is shown in Figure 5.3 and can be described as follows: ¯ The fall is realized at A where the lean line angle in Figure 5.2 reaches φ = φA . The first set of impulsive inputs are then applied to reach some desired velocities for passive coordinates which takes the robot configuration to B. After the application of the first set of impulsive inputs, the joint torques are set to zero such that the robot undergoes a free fall till it reaches Free fall with no torque inputs First Impulse to make desired velocity jumps in passive coordinates C : (qC , qC ) ˙ B : (qB , qB ) ˙ A : (qA , qA ) ˙ Fall is recognized at A Second Impulse to brake active coordinates D : (qD , qD ) ˙ Free fall with all active joints locked E : (qE , qE ) ˙ Desired Configuration on floor Start point: (Robot at upright configuration) Figure 5.3: A schematic of the safe fall algorithm using two sets of impulsive inputs. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation. 98 ¯ the configuration at C which is determined by φ = φC . The second set of impulsive torques are applied then to brake the active coordinates. This corresponds to a jump in the robot configuration form C to D. The active coordinates are held locked thereafter till the robot ¯ hits the floor at φ = φE . In particular, the following steps are adopted to design the safe fall algorithm: A. Off-line calculations to find a set of solutions at B: 1. Knowing the desired configuration on the floor and having all active coordinates locked, integrate the equations of motion backward in time to find the configuration at D : (qD , qD ). This configuration is realized when the lean line angle φ, ˙ ¯ shown in Figure 5.2, reaches φ = φC . Since the control objective is to change the fall direction which is mainly described by the angle ψ - see Figure 5.2, we have some freedom in choosing the angular positions of the active coordinates and angular velocities of the passive coordinates at the desired configuration on the floor. Each of these configurations will result in a different solution for D. 2. For each of the solutions at D obtained in step A.1, use the impulsive relations to find the robot configurations at C : (qC , qC ). In particular, from Eq.(5.11) we ˙ have:       ˙ ˙ ˙ θ2C  αC  αD              −1 ˙ ˙ ˙ C  =  βD  − M1 (qD ) M2 (qD ) θ3C  β             ˙4C θ γD ˙ γC ˙ (5.14) ˙ ˙ ˙ where θ2D = θ3D = θ4D = 0. Many solutions can be found for C by choosing a ˙ ˙ ˙ range of velocities for the active joint coordinats, θ2C , θ3C , and θ4C . Also, having different solutions for D obtained in step A.1 results in additional solutions for C. 99 Note that the the angular positions of the active and passive coordinates remain constant during the jump from D to C. 3. For each of the solutions at C, integrate the equations of motion backward in time ¯ to reach the configuration at B : (qB , qB ) which is realized when φ = φA . This is ˙ done assuming a free fall for the robot, i.e., control inputs are zero. Remark 5.4.1. After this step, we have a set of solutions for B which are resulted from different choices of desired configuration at E and different choices of velocities at C. Each of these solutions will indeed result in the desired change in fall direction of the robot, i.e., desired ψ angle when the robot hits the floor. B. On-line calculations to change the fall direction: ¯ 1. Once the fall is realized at φ = φA , among the set of solutions for B obtained ˙ in step A, find the closest solution QB = [qB , q B ] to the current configuration QA = [qA , qA ] where the distance si from each solution QBi = [qBi , qBi ] is defined ˙ ˙ as: si = ||QA − QBi ||2 (5.15) 2. Apply the impulsive inputs to reach the desired velocities for the passive coordi˙ ˙ ˙ nates, (αB , βB , γB ) corresponding to QB . The impulsive relations in Eq.(5.12) ˆ ˆ ˆ can be used to find the velocity jumps in the active coordinates as follows:       ˙ ˙ ˙ ˆ ˙ αB − αA  θ2B  θ2A        ˙      ˆ ˙ ˙ ˙3B  = θ3A  − M2# (qA ) M1 (qA )  βB − βA  θ             ˙ B − γA ˙4A ˙4B γ ˆ ˙ θ θ 100 (5.16) 3. Let the robot have a free fall with applying no torque in the active joints till it gets to C. ¯ 4. When φ = φC , apply the impulsive torques to brake the active coordinates and ¯ hold them fixed till the humanoid hit the floor at φ = φE . Remark 5.4.2. The application of the first set of impulsive inputs in step B.2 will not exactly result in the desired solution QB . This is because: first, the application of impulsive inputs does not change the angular positions and this yields a deviation from QB if qA = qB and second, the angular velocities resulted from the application ˙ of impulsive inputs will not match the angular velocities q B if qA = qB or the matrix M2 (qA ) is singular. Remark 5.4.3. Following the algorithm in step B, there is no guarantee that we reach the exact desired configuration on the ground. This is mainly because we can not exactly reach a desired solution for the configuration at B through applying the first set of impulsive inputs as discussed in Remark 5.4.2. However, we can get close to the desired change in the fall direction of the robot as we get closer to a desired solution at B, obtained in off-line calculations in step A. Our chance of reaching the desired change in the fall direction will be higher for a larger set of of-line solutions for B as is discussed in Remark 5.4.1. 5.5 Numerical Simulations The parameters for the robot in Figure 5.1 are listed in Table 5.2. The control parameters used in the safe fall algorithm are chosen as: ¯ φA = 30 deg, ¯ φC = 70 deg, 101 ¯ φE = 90 deg (5.17) A wide range of velocities are chosen for passive coordinates in Step A.1 and the active coordinates in step A.2 which result in a large set of solutions for B. The desired configuration of the robot on the floor which is used in off-line calculations in Step A, is chosen such that the desired change in the fall direction is 20 deg. The initial configuration of the robot are assumed to be: q0 =[α0 , β0 , γ0 , θ20 , θ30 , θ40 ]T = [90, 10, 0, −10, 0, 0]T ˙ ˙ ˙ ˙ ˙ q0 =[α0 , β0 , γ0 , θ20 , θ30 , θ40 ]T = [0, 5, 0, 0, 0, 0]T ˙ ˙ (5.18) where the units are in deg and deg/s. The simulation results are shown in Figures 5.4, 5.5 and 5.6. The joint positions and velocities of all the coordinates are shown in Figure 5.4 where the impulse moments can be seen from the jumps in the velocity subplots. The lean line angles φ and ψ are shown in Figure 5.5 and Figure 5.6 shows the impulsive inputs for the active coordinates. Starting from the initial configuration, the robot goes under a free fall with zero torques till the fall is realized by the algorithm at φ = 30 deg. The first impulsive inputs are then applied in accordance with Step B.2. The joint torques are set to zero thereafter till the robot reaches the second impulsive instant at φ = 70 deg when the active joints are brought to full stop. The active coordinates will remain locked using a locking mechanism till the robot hit the floor at φ = 90 deg. Table 5.2: Kinematic and dynamic parameters for the six-dof robot with cylindrical links. First link Second link Third link Fourth link Mass (Kg) Length (m) 1 0.5 1 0.5 4 0.8 2 0.3 Radius (m) Distance from C.M. (m) 0.1 0.25 0.1 0.25 0.25 0.4 0.05 0 102 The change in the fall direction can be seen from the ψ plot in Figure 5.5 where it is shown the change in the fall direction angle is close to the desired change angle, ψdes = 20 deg. Note that the dashed lines in Figures 5.4 and 5.5 show the simulation results when a single set of impulsive inputs are used to brake the active joints at φ = 30 deg. Remark 5.5.1. Although the fall direction has been changed using the proposed algorithm, the impulsive inputs result in very large velocity jumps in active coordinates which may cause the robot joints to exceed their limits. This is mainly because it is physically hard to change the fall direction of the robot with current active coordinates and very large impulsive inputs for active coordinates are needed to reach desired changes in the passive coordinates. 5.6 Conclusion A safe fall algorithm was proposed to change the default fall direction of a six-dof robot with three active joints. This is done to avoid hitting the surrounding objects to minimize damage to them. The algorithm take over when the regular control has failed to keep the robot balanced and a fall is unavoidable. Since there is a short period of time before the robot hit the floor, only two sets of impulsive inputs are applied at specific configurations of the robot. The first set of impulsive inputs is applied to give desired velocity jumps to the passive coordinates while the second impulsive set is intended to stop the active coordinates. Simulation results show the effectiveness of the algorithm in changing the fall direction of the robot. For the robot considered here, large impulsive inputs are needed to have a significant change in the fall direction which may cause large velocity jumps and exceeding joint limits for the active coordinates. This issue can be addressed for higher-dof robots with more active joints which make it physically easier to change the fall direction. 103 140 120 100 100 150 α (deg) α (deg/s) ˙ a 0 200 β (deg) 50 0 400 γ (deg) c θ2 (deg) i j 0 θ3 (deg) θ4 (deg) -100 0 γ (deg/s) ˙ ˙ 0 θ2 (deg/s) -200 d 0 -600 100 h 0 0 0 -50 ˙ β (deg/s) b 80 40 g 0.3 0.6 time (sec) e -2000 ˙ θ3 (deg/s) k 2000 f 0 0.9 ˙ -2000 θ4 (deg/s) 0 0.3 0.6 time (sec) l 0.9 Figure 5.4: Simulation results showing all the joint angles and angular velocities for changing the fall direction of the robot 104 100 φ (deg) 80 a 60 40 20 0 20 ψ (deg) b 10 0 0 0.3 time (sec) 0.6 0.9 Figure 5.5: Simulation results showing the lean line angles for changing the fall direction of the robot τ2 (N.m) 150 a 0 τ3 (N.m) 60 b τ4 (N.m) -60 5 c -10 0 0.3 time (sec) 0.6 0.9 Figure 5.6: Simulation results showing the impulsive inputs for changing the fall direction of the robot 105 Chapter 6 Disturbance Rejection for the Synthetic-Wheel Biped 6.1 Introduction In this chapter, we propose a push recovery strategy for the Synthetic-Wheel Biped (SWB) [1]. Stepping strategy is used to maintain balance when the biped is pushed due to some external disturbances. The external disturbances are impulsive in nature, i.e., cause some sudden jumps in the system velocities. To reject these disturbances, our control algorithm uses impulsive inputs. The push recovery algorithm is pretty robust and capable of rejecting any pushes during stance and walking phases. Impulsive inputs are also used to impose some virtual constraints on the motion of the legs and the torso. This work is mainly different from the approaches in the literature because first, the SWB is an underactuated biped with curved feet and no ankle joints and second, the set of admissible controls are enlarged to include impulsive inputs. This chapter is organized as follows: The dynamics and control of the Synthetic-Wheel biped is reviewed in Section 6.2. We derive the expressions for impulsive inputs and their effects on system velocities in Section 6.3. Disturbance rejection algorithm is discussed in Section 6.4 and two numerical examples are presented in Section 6.5. We provide the concluding remarks in Section 6.6. 106 6.2 6.2.1 Background Equations of Motion The SWB [1] is shown in Figure 6.1. It has a torso and two legs with arc-shaped feet. The arc-shaped feet have the same radius as the leg length and therefore the functionality of a wheel is maintained by proper placement of the swing leg. The biped has a mechanism that enables the swing leg to be slightly shorter in length than the stance leg and thereby avoid scuffing with the ground. The dynamics of this mechanism is insignificant compared to the overall dynamics of the system. The SWB is characterized by three generalized coordinates: θ, φ, and ψ. The coordinates φ and ψ are active and are controlled by two actuators; the coordinate θ is passive. The equations of motion are obtained as follows, M(q)¨ + B(q, q)q + G(q) = T q ˙ ˙ (6.1) dt lt torso R x β y dsw φ ψ g R R swing foot θ dst stance foot Figure 6.1: A schematic of the Synthetic-Wheel biped, reproduced from [1]. 107 where q = [θ φ ψ]T , T = [0 τ2 τ3 ]T (6.2) The matrices M(q), N(q, q) and G(q) are represented as: ˙ M(q) = Mij 3×3 , N(q, q) = Nij 3×3 , G(q) = [Gi ]3×1 ˙ (6.3) and are defined in Appendix B. 6.2.2 Torso Effect on Biped Dynamics The specific design of the SWB provides the ability to effectively change the direction of the walking through changing the torso angle. To demonstrate this, consider a “symmetric gait” of the biped which is generated by imposing the following constraints on the motion of the torso and the swing leg [1]: C1 : α = αdes (6.4) C2 : ψ = −2θ where α, depicted in Figure 6.2, is the angle of torso with respect to the vertical axis and is defined as α= θ+φ−π (6.5) The constraint C1 ensures that the torso maintains a desired angle αdes with respect to the vertical axis while constraint C2 ensures that the swing leg is symmetric with respect to the stance leg about the vertical axis at all times. The constrained system has one passive dof 108 with the following dynamics, which is derived by substituting Eq.(6.4) into Eq.(6.1): ¨ ˙ ˙ Mc (θ) θ + Nc (θ, θ) θ + Gc = 0 (6.6) ˙ where the expressions for Mc (θ), Nc (θ, θ) and Gc can be found in [1]. For any set of reasonable ¨ parameter values, it can be verified that θ will be positive for positive angle αdes and vice ˙ versa. This implies that we can slow down the biped velocity θ by applying the constraints ˙ in Eq.(6.4) and choosing αdes to have the opposite sign of θ. 6.2.3 Interchange of Stance and Swing Legs Assuming that the biped has a positive velocity, it will roll on its stance leg and the point of contact with the ground will move from the heel to the toe. The stance and swing legs can be interchanged at any time; but to have the maximum step size, the swing leg should torso α φ = −θ+π+α X Y β swing leg toe θ −ψ=2θ swing leg heel stance leg heel stance leg toe Figure 6.2: Biped configuration at the time of interchange of stance and swing legs. 109 touch down when the heel of the swing leg is right in front of the toe of the stance leg, as seen in Figure 6.2. Without assuming maximum step size, an interchange of stance and swing legs will result in a transformation of the generalized coordinates and their velocities given by the following relations [1]: qnew = P qold (6.7) qnew = P qold ˙ ˙ where qold and qnew denote the values of the vector q before and after the interchange of stance and swing legs, respectively, and the entries of the transformation matrix P are   1   1 0     P =  0 1 −1      0 0 −1 (6.8) which can be obtained from Figure 6.2. The dynamic model in Eq.(6.1) involves kinematic and dynamic parameters of the stance and swing legs and their values need to be interchanged during interchange of stance and swing legs as well. 6.2.4 Control Design We develop the controller for imposing the constraints described in section 6.2.2. For the symmetric gait with constraints in Eq.(6.4), the new set of generalized coordinates is defined as T q= ¯ θ v1 v2 110 (6.9) where v1 = α − αdes v2 = ψ + 2 θ (6.10) The original generalized coordinates φ and ψ are related to the new coordinates v1 and v2 according to the following relations, which can be derived from Eqs.(6.5) and (6.10): φ =v1 − θ + π + αd ψ =v2 − 2 θ (6.11) Substituting Eq.(6.11) into Eq.(6.1), we obtain the dynamics of the system in terms of the new generalized coordinates as follows: ¯ q ¨ ¯ q ¯ ¯ ¯ q ˙ ˙ M (¯) q + N (¯, q ) q + G(¯) = T ¯ (6.12) ¯ ¯ ¯ where M , N , and G have the same dimensions as M, N, and G respectively. The generalized ¨ force corresponding to θ in Eq.(6.12) is zero and this allows us to eliminate θ from the two equations corresponding to the generalized coordinates v1 and v2 . The reduced-order equations have the form ˆ q ¨ ˆ q ¯ ˆ ˆ q ˆ ˙ ˙ M (¯) q + N (¯, q ) q + G(¯) = T ˆ (6.13) ˆ ˆ ˆ where M ∈ R2×2 , N ∈ R2×2 and G ∈ R2×1 , and q = [v1 v2 ]T , ˆ ˆ T = [τ2 τ3 ]T 111 (6.14) Equation (6.13) represents a completely actuated system and we use feedback linearization to design our controller as follows: ˙ ˆ ˆ q ¯ˆ ˆ q ˆ q ˙ ˙ T = N (¯, q )q + G(¯) − M (¯)(Kd q + Kp q ) ˆ ˆ (6.15) where Kd and Kp are diagonal, positive-definite matrices of dimension two. Indeed, substitution of Eq.(6.15) into Eq.(6.13) yields ¨ ˙ q + K d q + Kp q = 0 ˆ ˆ ˆ (6.16) which implies q → 0 as t → ∞. This simply follows that v1 , v2 → 0 as t → ∞, i.e., the ˆ constraints in Eq.(6.4) are satisfied. The controller in Eq.(6.15) has been implemented and experimentally verified in our laboratory biped [1]. The experimental results show that v1 and v2 do not converge to zero but oscillate around zero due to impulsive disturbances from the ground at the time of swing-leg touchdown. The algorithm presented in this chapter provides a way to avoid these impulsive disturbances from the ground while taking steps due to an external disturbance. 112 6.3 Impulsive Torques and Effects 6.3.1 Braking Torque for the Swing Leg ˙ Consider the action that results in exponential convergence of the swing leg velocity ψ to ˙ zero while keeping θ unchanged. To this end, the following dynamics are assumed, ¨ θ =0 ˙ ¨ ψ = − k1 ψ (6.17) ˙ where k1 is a positive constant determining the rate of convergence of ψ to zero. To compute the torque required for this action, we multiply Eq.(6.1) with the inverse of the inertia matrix to obtain   ¨1  θ    1  ¨   θ2 =   A33 (A11 A22 − A2 ) − A22 A2 12 13   ¨ θ3    h1 − A12 A33 τ2 − A13 A22 τ3       h2 +(A11 A33 −A2 )τ2 +A12 A13 τ3  13     h3 +A12 A13 τ2 +(A11A22 −A2 )τ3 12 (6.18) where h1 , h2 and h3 are given by the following expressions, ˙ ˙ ˙ h1 = − A22 A33 (B11 θ + B12 φ + B13 ψ) − A22 A33 G1 + A12 A33 G2 + A22 A13 G3 ˙ ˙ ˙ h2 =A12 A13 (B11 θ+B12 φ+B13 ψ) +A12 A33 G1 +(A2 −A11 A33 )G2 −A12 A13 G3 13 ˙ ˙ ˙ h3 =A22 A13 (B11 θ + B12 φ + B13 ψ) +A13 A22 G1 −A12 A13 G2 +(A2 −A11 A22 )G3 12 Substituting Eq.(6.17) into the first and third equations in Eq.(6.18) results in 113 (6.19) τ2 = 1 A2 h1 −A22 (A11 h1 +A13 h3 ) 12 ˙ +k1A13 A22 ψ A12 A2 A22 +(A2 −A11 A22 )A33 13 12 A13 h1 + A33 h3 ˙ τ3 = 2 − k1 A33 ψ A13 A22 + (A2 − A11 A22 )A33 12 (6.20) If the constant k1 has a very large value, the torque expressions in Eq.(6.20) will be impulsive ˙ in nature and will stop the swing leg relative velocity, ψ, in a very short period of time without changing the stance leg velocity. 6.3.2 Braking Torque for the Torso Consider an action which results in exponential convergence of the torso absolute velocity, α, to zero while maintaining equal and opposite velocities for the stance and swing legs, i.e. ˙ ˙ ˙ ψ = −2 θ, in conformity with the symmetric gait. Using Eq.(6.5), we consider the following dynamics to achieve the goal, ¨ ¨ ˙ ˙ θ + φ = − k2 (θ + φ) ¨ ¨ ˙ ˙ θ + 2ψ = − k3 (θ + 2ψ) (6.21) where k2 and k3 are some positive constants determining the rate of convergence of the velocities to their desired values. The control inputs that result in the dynamics in Eq.(6.21) can be obtained by substituting Eq.(6.18) into Eq.(6.21) and solving for τ2 and τ3 . The massive torque expressions are not explicitly mentioned here for the sake of simplicity. If we choose very large values for k2 and k3 in Eq.(6.21), the resulted torque expressions will be impulsive in nature and stop the motion of the torso in a very short period of time without 114 affecting the symmetric velocity condition associated with the symmetric gait. 6.3.3 Impulsive Effect on Velocity The application of impulsive torques for φ and ψ generalized coordinates results in velocity jumps in all three coordinates θ, φ and ψ. The relationship between the jumps in velocities can be derived from Lagrange’s equations. Consider the integral of the equations of motion in Eq.(6.1) over the short interval of time ∆t during which the impulsive forces and moments act, i.e. ∆t ∆t [M(q)¨ + N(q, q)q + G(q)] dt = q ˙ ˙ 0 T dt (6.22) 0 The above equation can be rewritten as: ∆t M(q)∆q + N(q, q)∆q + ˙ ˙ ∆t G(q) dt = 0 T dt (6.23) 0 Since ∆t is very short time interval and the configuration of the biped does not change during this time, i.e. ∆q = 0, Eq.(6.23) will be simplified to M(q)∆q = Tim ˙ (6.24) where Tim = [0 I2 I3 ]T represents the vector of impulses applied to the biped. The above equation can be decomposed into the following three equations: ˙ ˙ ˙ ˙ ˙ ˙ M11 (θ+ − θ− )+M12(φ+ − φ− )+M13 (ψ + − ψ − ) = 0 (6.25) ˙ ˙ ˙ ˙ ˙ ˙ M21 (θ+ − θ− )+M22(φ+ − φ− )+M23 (ψ + − ψ − ) = I2 (6.26) 115 ˙ ˙ ˙ ˙ ˙ ˙ M31 (θ+ − θ− )+M32(φ+ − φ− )+M33 (ψ + − ψ − ) = I3 (6.27) where the superscripts “−” and “+” denote the values immediately before and after the application of the impulsive inputs, respectively. As the first special case, consider the impulsive action described in Section 6.3.1 in which ˙ ˙ ˙ ψ + = 0 and θ+ = θ− . Using Eq.(6.25), the relative velocity of torso after the impulsive action can be obtained as M ˙ ˙ ˙ φ+ = φ− + 13 ψ − M12 (6.28) As the second case, consider the impulsive action in section 6.3.2 which results in α+ = 0 ˙ ˙ ˙ and ψ + = −2 θ+ . Then the velocity of the swing leg after the impulse can be obtained from Eq.(6.25) as follows: ˙ ˙ ˙ M θ− + M12 φ− + M13 ψ − ˙ θ+ = 11 M11 − M12 − 2M13 (6.29) It should be noted that the impulses I2 and I3 can be computed from Eqs.(6.26) and (6.27). These impulses can also be approximated by the time integral of the torque expressions obtained in sections 6.3.1 for large gain k1 and the torque expressions in section 6.3.2 for large gains k2 and k3 . 6.4 Disturbance Rejection Algorithm In this section, we develop an algorithm to reject the external disturbances applied to the Synthetic-Wheel Biped. The external disturbances are assumed to be of the form of impulsive forces and their effects are modeled by jumps in the generalized velocities. The disturbances are assumed to be large enough that precludes the possibility of stabilizing the upright 116 posture without taking a step. The problem at hand is therefore to reject the disturbances and stabilize the upright posture by taking a few steps in the forward or backward direction. The disturbance rejection algorithm is based on slowing down the biped using the symmetric gait while imposing constraints on the legs and torso. These constraints are described as follows: C3 : − β/2 ≤ θ1 ≤ β/2 (6.30) C4 : − γ ≤ α ≤ γ (6.31) where β is the foot arc angle (see Figure 6.2) and γ is a small positive angle. We propose the following three-step algorithm to reject external impulsive disturbances and stabilize the upright posture of the Synthetic-Wheel Biped: 1. Initialization: (a) Linearize the dynamic equations in Eq.(6.1) about the desired equilibrium point ˙ ˙ ˙ (θ, θ, φ, φ, ψ, ψ) = (0, 0, π, 0, 0, 0). (b) Design a LQR controller for the linearized system to render the desired equilibrium locally asymptotically stable 1 . Let RA be the region of attraction of the desired equilibrium. (c) Define a quadratic Lyapunov function for the linearized system as V = X T SX ˙ ˙ ˙ where X = [θ, θ, φ − π, φ, ψ, ψ]T and S is the solution of the Riccati equation associated to the LQR design. T (d) Define a set point for the Lyapunov function as V0 = X0 SX0 where X0 = 1 This is always possible because the linearized system is controllable 117 [η, 0, αdes − η, 0, −2η, 0]T and η < β/2 is a positive angle. η and αdes values are chosen sufficiently small such that if V ≤ V0 , the biped configuration will lie inside RA and will approach the upright equilibrium without exceeding the bound on θ in Eq.(6.30). 2. Disturbance Rejection: (a) Following the external disturbance, apply the torque expression in Eq.(6.15) to implement the symmetric gait. The torque expression in Eq.(6.15) depends on the value of the torso desired angle, αdes , which is chosen to have the opposite ˙ sign of the stance leg velocity, θ, so the biped slows down. This follows from our discussion in Section 6.2.2. (b) If the stance leg angle θ reaches the boundary of the interval [−β/2, β/2], ap˙ ply the torque expressions in Eq.(6.20) to quickly enforce ψ = 0 while keeping the stance leg velocity unchanged. The stance leg and swing leg are immediately interchanged using Eq.(6.7) and the control torques in step 2(a) are applied there˙ after. By enforcing ψ = 0, we avoid impulsive disturbance from the ground to the biped (swing leg) at the time of leg interchange. It can be seen from Eq.(6.7) that ˙ this also eliminates jumps in θ at the time of leg interchange, which is necessary for smooth walking. (c) If the torso angle α reaches the boundary of the interval [−γ, γ], apply the torque expressions obtained in section 6.3.2 to quickly stop the torso while keeping the velocity of the legs symmetric. Then continue to apply the control torques in step 2(a). ˙ (d) If V ≤ V0 , terminate this step and go to step 3. The stance leg velocity, θ, will 118 ˙ go to zero since αdes was chosen to have a sign opposite to that of θ. Imposing ˙ ˙ the symmetric gait constraints in Eq.(6.4), φ and ψ will approach zero as well. Therefore, By choosing small enough value for αdes , the inequality V ≤ V0 is eventually satisfied. 3. Stabilization: With V ≤ V0 and the biped satisfying the constraints in Eq.(6.30) and Eq.(6.31), the biped configuration will be inside the region of attraction RA. Invoke the linear controller designed in step 1 to stabilize the desired equilibrium. 6.5 Numerical Simulations We present two simulations to demonstrate the effectiveness of the proposed algorithm in rejecting the impulsive disturbances on the Synthetic-Wheel biped. For both simulations, the kinematic and dynamic parameters in Table 6.1 are used. For the first simulation, we assume two impulsive disturbances are applied on the biped Table 6.1: Parameters of the Synthetic-Wheel Biped [1] Inner leg Outer leg Torso Inner leg Outer leg Torso Kinematic parameters Length (m) Foot radius, R (m) Foot arc, β (deg) 0.635 0.635 22.5 0.635 0.635 22.5 0.457 Dynamic parameters Mass (kg) Inertia (kgm2 ) d in Fig.2 (m) 1.64 0.094 0.285 3.64 0.128 0.355 11.87 0.198 0.307 119 at two instants of time which can be modeled by the following values of the angular velocities: ˙ + ˙+ ˙ + (θ1 , θ2 , θ3 ) = (0, −5, −1) , at t = 0 sec ˙ + ˙+ ˙ + (θ1 , θ2 , θ3 ) = (−0.1, 3, −1) , at t = 2 sec (6.32) The initial angular positions and the parameters needed in the algorithm are chosen as: (θ10 , θ20 , θ30 ) = (5◦ , 180◦ , −8◦ ) ˙ η = β/3, αdes = −sgn(θ) 5◦ , γ = 10◦ (6.33) Figure 6.3 shows the simulation results including the plots of the joint angles and angular ˙ velocities and the control inputs. The angular position ψ and the angular velocity ψ of the swing leg are shown by dashed lines. The horizontal lines in the angular position subplots represent the bounds on the stance leg and the torso. Due to the initial impulsive disturbance, the torso gets to its lower bound and the impulsive inputs are applied to keep them inside ˙ the allowed region. The biped starts walking in the negative direction with θ < 0 and the continuous controller will slow down the biped motion by imposing the symmetric gait constraints. when the stance leg gets to its bound, the impulsive inputs are applied to make ˙ ψ = 0 right before the interchange of the legs. When the second impulsive disturbance is applied at t = 2 sec., the biped changes its walking direction to positive and consequently, the continuous controller switches the desired angle of torso to negative to slow down the biped using the symmetric gait constraints. Torso reaches its upper bound after the second disturbance as well and is kicked back inside the region using impulsive inputs. The peaks in the torque plots correspond to impulsive torques which are applied when the stance leg 120 and the torso reach their bounds. As can be seen from Figure 6.3, the biped is able to reject the impulsive disturbance after taking few steps in the negative and positive directions. The linear controller is invoked at t = 6.04 sec. (shown by vertical dashed line) to stabilize the upright configuration of the biped. For the second simulation, we applied five identical impulsive disturbances to the biped at time instants t = 0, 2, 4, 6, 8 sec. which are modeled by the following changes in the velocities: ˙ ˙ ˙ τ3 (N.m) τ2 (N.m) α(rad/s) θ, ψ(rad/s) α(rad) θ, ψ(rad) ˙+ ˙ + ˙ + (θ1 , θ2 , θ3 ) = (−0.5, 6, −2) (6.34) 0.4 0 -0.4 0.17 0 -0.17 a b c 2 0 -2 2 0 -4 45 0 -40 120 0 -200 0 d e f 1 2 3 4 time (sec) 5 6 7 8 Figure 6.3: Simulation results showing angular positions (rad), angular velocities (rad/s), and control inputs (N.m) of the SWB in rejecting the impulsive disturbances applied at two ˙ instants of time. The dashed lines correspond to the plots of ψ and ψ. 121 the initial conditions for joint angles and the other algorithm parameters were chosen as: (θ10 , θ20 , θ30 ) = (0, π, 0) ˙ η = β/3, αdes = −sgn(θ) 5◦ , γ = 20◦ (6.35) The simulation results are shown in Figure 6.4. As in the first simulation, dashed lines represent the swing-leg angle and angular velocity and the horizontal lines in angular positions subplots show the bounds on the stance leg and torso motion. The disturbances make the biped to walk in the positive direction with symmetric gait to recover its upright configuration. In this case, the limit on the torso motion is not reached and all the impulsive θ, ψ(rad) 0.35 0 -0.35 ˙ ˙ τ3 (N.m) τ2 (N.m) α(rad/s) θ, ψ(rad/s) α(rad) ˙ inputs shown by peaks in the torque plots, are applied when the stance leg gets to its bound 0.35 0.17 0 2 0 -2 -4 6 2 -2 20 10 0 -10 100 50 0 0 a b c d e f 5 time (sec) 10 15 Figure 6.4: Simulation results showing angular positions (rad), angular velocities (rad/s), and control inputs (N.m) of the SWB in rejecting the impulsive disturbances applied at five ˙ instants of time. The dashed lines correspond to the plots of ψ and ψ. 122 at θ = β/2. Using the proposed algorithm, the biped is able to reject the disturbances. The vertical dashed line in Figure (6.4) represents the time t = 12.3 sec. when the linear controller is invoked to stabilize the upright equilibrium of the biped. 6.6 Conclusion We proposed an algorithm to reject the effects of impulsive disturbances applied to the Synthetic-Wheel Biped in its stance and walking configurations. The impulsive disturbance is modeled by jumps in system velocities. The algorithm uses a combination of continuous and impulsive control inputs to generate a symmetric gait for the SWB while imposing some constraints on the system states. Using the symmetric gait, the biped takes a few steps in order to reject the effects of the external disturbances and recover its balance. The constraints on the system are imposed to avoid additional impulsive disturbances from the ground at the interchange of legs. Simulation results show the effectiveness of the algorithm in rejecting the impulsive disturbances applied at different instants of time in both stance and walking phases of the biped. 123 Chapter 7 Energy-Conserving Gaits for Point-Foot Passive-Ankle Bipeds 7.1 Introduction In the previous chapter, we studied a case study where impulsive inputs in the form of external disturbance are undesirable and the control objective is to cancel their effects on the system. In this chapter, we consider another undesirable form of impulsive inputs which exists in walking bipeds: the interaction impact from ground applied during the legs interchange. In addition to the velocity jumps for the biped links, energy is lost during the impact. To make a bipedal walking gait more energy efficient, controllers need to be designed to avoid these impacts. In this chapter, we present a general method for designing energy-conserving gaits for point-foot passive-ankle planar bipeds. Assuming no energy dissipation due to friction, an energy-conserving gait is one that has a zero Cmt based on the definition in Eq.(1.1). The sufficient conditions for such a gait are equal potential and kinetic energies at the beginning and end of each step and no energy loss due to impact at the time of swing-leg touchdown. These conditions are presented in Section 7.2. The conditions are applied to two three-dof bipeds in Section 7.3; the results indicate that energy-conserving gaits can be designed for the Synthetic-Wheel biped [1] but not for the biped with point feet. Energy-conserving gaits 124 are designed for five-dof bipeds in Sections 7.4 and 7.5. One-step periodic gaits are designed in Section 7.4 for a biped with similar leg links and more general two-step periodic gaits are designed for bipeds with similar and dissimilar leg links in Section 7.5. Concluding remarks are provided in Section 7.6. 7.2 Energy-Conserving Gaits 7.2.1 Assumptions and Definitions Consider the general n-link planar biped comprised of a torso and two legs, shown in Figure 7.1. The equations of motion of the biped are assumed to be of the form M(q) q + N(q, q) q + G(q) = T ¨ ˙ ˙ (7.1) where q = (θ1 , θ2 , · · · , θn )T and T = (u1 , u2, · · · , un )T are the vectors of generalized coordinates and generalized forces1 , M(q) ∈ Rn×n is the symmetric mass matrix, N(q, q) ∈ Rn×n ˙ is the matrix of centrifugal and Coriolis terms, and G(q) ∈ Rn is the vector of gravitational forces. We assume that A1. Only one leg of the biped is in contact with the ground at any given time, i.e., there is no double-support phase. The leg in contact with the ground is commonly referred to as the stance leg; the other leg is referred to as the swing leg. All joints are measured counter-clockwise with respect to the vertical. They are numbered arbitrarily except that θ1 denotes the angle of the stance 1 Since the generalized coordinates are absolute joint angles, the generalized forces are the net external moments acting on each link. 125 torso θj hip joint θn θ2 y θ1 swing foot stance foot x Figure 7.1: A schematic of an n-link point-foot planar biped. leg link in contact with the ground, and θn denotes the angle of the swing leg link that will come in contact with the ground. For any biped gait, we define the term “step” as follows: Definition 7.2.1. A step is the complete motion of the biped on the stance leg which is followed by interchange of the stance and swing legs. We make the following additional assumptions: A2. The stance and swing legs of the biped have point feet and the stance foot does not slide or leave the ground during a step. A3. All joints of the biped are actuated except the ankle joint of the stance leg, which is passive. A4. The biped walks on level ground and there is no energy dissipation due to friction. 126 At the beginning of each step, the joint angles and their velocities are denoted by i i i q i = (θ1 , θ2 , · · · , θn )T , ˙i ˙i ˙i q i = (θ1 , θ2 , · · · , θn )T ˙ Similarly, the joint angles and their velocities at the end of each step are denoted by f ˙f ˙f ˙ f q f = (θ1 , θ2 , · · · , θn )T ˙ f f q f = (θ1 , θ2 , · · · , θn )T , We also assume that: A5. The joint angles of the stance leg ankle at the beginning and end of each step are f i symmetric with respect to the vertical, i.e., θ1 = −θ1 . The velocity of the swing foot at the time of contact with the ground (at the end of each step) is denoted by f f f Vn = Vn,x ˆ + Vn,y ˆ i j f f where ˆ and ˆ are unit vectors along the inertial x and y axes, and Vn,x and Vn,y are the i j corresponding velocity components. 7.2.2 Sufficient Conditions for Energy-Conserving Gait For a biped, the mechanical energy required by its gait will be zero over each step if the following conditions are satisfied: 1. The potential energies at the beginning and end of each step are identical. This is achieved if f i θj = −θj , 127 j = 1, 2, · · · , n (7.2) 2. The kinetic energies at the beginning and end of each step are identical. This is achieved if ˙i ˙f θj = θj , j = 1, 2, · · · , n (7.3) Additionally, during interchange of the stance and swing legs, the energy of the biped will remain conserved if: 3. Swing foot touchdown does not cause loss of energy. This is achieved if there is no impact due to foot-ground interaction, i.e. f f Vn,x = Vn,y = 0 (7.4) The conditions in Eqs.(7.2) and (7.3) can be easily satisfied for θj , j = 2, 3, · · · , n, since all joints of the biped except the stance leg ankle are actuated. The condition in Eq.(7.2) is satisfied for θ1 as well - this is a direct consequence of assumption A5. Since the stance leg ˙ ankle is passive, the main challenge is to satisfy the condition in Eq.(7.3) for θ1 . 7.2.3 Characterization of Trajectories of the Actuated Joints ˙ The condition in Eq.(7.3) can be satisfied by the ankle joint velocity θ1 if it is an even function2 of θ1 , i.e., e ˙ θ1 = g1 (θ1 ) (7.5) This follows from assumption A5. We now state necessary conditions for the actuated joint trajectories with the help of the following Theorem. 2 We use superscripts “e” and “o” to denote even and odd functions, respectively. 128 Theorem 7.2.1. The conditions in Eqs.(7.2), (7.3) and (7.5) will be satisfied only if the velocities of the actuated joints are even functions of θ1 , i.e. e ˙ θj = gj (θ1 ), j = 2, 3, · · · , n (7.6) ˙ Proof. To guarantee θ1 in Eq.(7.5) is a function of θ1 alone, all actuated joints must be functions of θ1 , i.e. θj = fj (θ1 ), j = 2, 3, · · · , n (7.7) ∂fj , ∂θ1 (7.8) Taking the derivative of Eq.(7.7), we get ˙ ˙ θj = hj (θ1 ) θ1 , hj j = 2, 3, · · · , n The hj ’s can be written as a sum of even and odd functions as follows ˙ ˙ θj = he (θ1 ) + ho (θ1 ) θ1 , j j j = 2, 3, · · · , n (7.9) From Eq.(7.3), we have for j = 2, 3, · · · , n f f i i ˙i ˙f he (θ1 ) + ho (θ1 ) θ1 = he (θ1 ) + ho (θ1 ) θ1 j j j j ⇒ f f i i he (θ1 ) + ho (θ1 ) = he (θ1 ) + ho (θ1 ) j j j j 129 (7.10) Rearranging terms and using Eq.(7.2) and the property of even and odd functions, we get f f i i ho (θ1 ) − ho (θ1 ) = he (θ1 ) − he (θ1 ) j j j j ⇒ i i i i ho (−θ1 ) − ho (θ1 ) = he (θ1 ) − he (−θ1 ) j j j j ⇒ ⇒ i i ho (−θ1 ) = ho (θ1 ) j j ho (θ1 ) = 0 j (7.11) It simply follows that hj (θ1 ) = he (θ1 ), j = 2, 3, · · · , N. Using Eqs.(7.5) and (7.8) we can j now write e ˙ ˙ θj = he (θ1 ) θ1 = he (θ1 ) g1 (θ1 ) j j e gj (θ1 ) (7.12) ˙ which implies that θj , j = 2, 3, · · · , n is even function of θ1 . This completes the proof. Since the derivative of an odd function is always an even function, the theorem above implies that the trajectories of the actuated joints consist of an odd function and a constant. We choose the following general form for the trajectories o θj (θ1 ) = fj (θ1 ) + aj θ1 + bj , j = 2, 3, · · · , n (7.13) where aj and bj are constants, and the term aj θ1 is the linear part of the odd function. The constants bj and the coefficients of the odd functions in Eq.(7.13) will be chosen such that the boundary conditions in Eqs.(7.2), (7.3) and (7.4) are satisfied. Of these boundary ˙ conditions, the condition on θ1 in Eq.(7.3) is satisfied indirectly by choosing the constants and ¨ the coefficients such that θ1 is an odd function. This requires the substitution of Eq.(7.13) in ¨ the expression for θ1 , which is obtained from the equation of motion for θ1 . Although there 130 is no guarantee that trajectories of the form given by Eq.(7.13) can be found for energyconserving gaits for arbitrary bipeds, we show that such trajectories do exist for several planar bipeds. 7.3 Three-DOF Bipeds: Two Case Studies In this section, we design energy-conserving gaits for two 3-dof biped platforms: the SyntheticWheel biped3 [1] and a point-foot three-link biped with two legs and a torso. 7.3.1 Synthetic-Wheel Biped (SWB) 7.3.1.1 System Description The SWB [1], which was described in Section 6.2.1, is shown in Figure 7.2. The equations of motion of the SWB have the same form as Eq.(7.1), with q and T defined as follows: q = [θ φ ψ]T , T = [0 τ2 τ3 ]T (7.14) where τ2 and τ3 are motor torques corresponding to the relative coordinates φ and ψ. The equations of motion of the SWB are provided in Appendix B. Remark 7.3.1. The generalized coordinates in Eq.(7.14) are different from those chosen for the general case in the previous section. This is because the dynamics of the SWB, already presented in [1], is based on relative generalized coordinates shown in Figure 7.2. This difference does not change the sufficient conditions for an energy-conserving gait, namely, 3 The Synthetic-Wheel biped has arc-shaped feet but is investigated here along with point-foot bipeds because of its point contact with the ground. 131 dt lt torso R x β y dsw φ ψ g swing foot R R θ dst stance foot Figure 7.2: A schematic of the Synthetic-Wheel biped, reproduced from [1]. the kinetic and potential energies at the beginning and end of each step are the same, and there is no loss of energy at the time of swing-leg touch down. 7.3.1.2 Sufficient Conditions for Energy-Conserving Gait The boundary conditions for the gait are determined as follows: 1. The potential energies at the beginning and end of each step are identical. This is achieved by θf = −θi , φf = −φi , ψ f = −ψ i (7.15) 2. The kinetic energies at the beginning and end of each step are identical. This is achieved by ˙ ˙ θf = θi , ˙ ˙ φf = φi , ˙ ˙ ψf = ψi (7.16) 3. The synthetic-wheel design ensures that there are no vertical impulsive forces, i.e., 132 f Vn,y = 0. At the time of swing foot touchdown, there will be no horizontal impulsive forces if ˙ ψf = 0 7.3.1.3 (7.17) Determination of Actuated Joint Trajectories Using Eq.(7.13), we consider the following functional dependence for φ and ψ o φ(θ) = (a1 θ + b1 ) + f1 (θ) (7.18) o ψ(θ) = (a2 θ + b2 ) + f2 (θ) (7.19) ˙ ˙ ¨ The condition θi = θf in Eq.(7.16) can be satisfied if θ is an odd function of θ. To this end, ¨ we derive the expression for θ as follows ¨ N θ= D (7.20) where N = N(q, q, q ) and D = D(q) are defined as follows ˙ ¨ ¨ N(q, q, q ) =− φ {It +mt (dt −lt )(dt −lt +R cos [θ + φ])} ˙ ¨ ¨ − ψ {Isw +msw (dsw −R)(dsw −R+R cos [θ+ψ])} ˙ ˙ ˙ + mst (dst − R)(g + Rθ2 ) sin θ + mt (dt − lt ) g + R(θ + φ)2 sin [θ + φ] ˙ ˙ +msw (dsw −R) g+R(θ+ ψ)2 sin [θ+ψ] 133 D(q) = It + Ist + Isw + mst (dst − R)2 + R2 − 2R(−dst + R) cos θ 2 + mt d2 +lt +R2 +2R(dt −lt ) cos(θ + φ)−2dt lt t +msw (dsw −R)2 +R2 −2R(−dsw +R) cos [θ + ψ] ¨ Using Eqs.(7.18) and (7.19), a term-by-term analysis of N and D shows that θ is an odd function of θ for any values of a1 and a2 provided that b1 and b2 are chosen as [138] b1 = mπ, m, n = 0, ±1, ±2, · · · b2 = nπ, (7.21) For the torso to be upright, the choice of m for b1 = mπ should be limited to m = ±1, ±3, ±5, · · · . Similarly, for the swing-leg to point downwards, the choice of n for b2 = nπ should be limited to n = 0, ±2, ±4, · · · . By simply setting b1 = π and b2 = 0 and assuming o o f1 and f2 in Eqs.(7.18) and (7.19) to be sinusoidal functions, we have φ(θ) = (a1 θ + π) + A sin (Bθ) (7.22) ψ(θ) = a2 θ + C sin (Dθ) (7.23) If the SWB moves in the positive x direction with maximum step length, we have θf = −θi = β 2 ψ f = −ψ i = −β (7.24) (7.25) where β is the arc angle of the feet - see Figure 7.2. Using Eq.(7.24) and the boundary conditions corresponding to θ and φ in Eqs.(7.15) and (7.16), the coefficients a1 and B can 134 be determined. Using Eqs.(7.24) and (7.25) and the boundary conditions corresponding to θ and ψ in Eqs.(7.15), (7.16) and (7.17), the coefficients a2 , C, and D can be determined. The energy-conserving trajectories of the actuated joints finally take the form φ(θ) = π − θ + A sin(2πθ/β) (7.26) ψ(θ) = −β sin(πθ/β) (7.27) The choice of the constant A is discussed in the next section. 7.3.1.4 Numerical Simulations The kinematic and dynamic parameters of the SWB are listed in Table 7.1. The constant A corresponds to the amplitude of oscillation of the torso and determines if the biped will succeed or fail to take a step. To illustrate this point, we plot phase portraits of the motion of the SWB in Figure 7.3 for different values of A with initial conditions θi = −β/2 = −11.25 ˙ ˙ deg and θi = 45 deg/s. For A equal to 0.10 and 0.08, the SWB fails to take a step because θ changes sign before the step is completed. As A is reduced, for A equal to 0.06 and smaller, Table 7.1: Kinematic and dynamic parameters of Synthetic Wheel Biped Leg 1 Leg 2 Torso Leg 1 Leg 2 Torso Kinematic parameters lst , lsw , lt R 0.635 m 0.635 m 0.635 m 0.635 m 0.457 m Dynamic parameters mst , msw , mt Ist , Isw , It 2.640 kg 0.111 kg m2 2.640 kg 0.111 kg m2 11.87 kg 0.198 kg m2 135 β 22.5 deg 22.5 deg dst , dsw , dt 0.320 m 0.320 m 0.307 m 100 θi = −11.5 ˙ θi = 45 -0.40 -0.20 0.0 0 ˙ θ (deg/s) 50 0.06 0.10 0.02 0.04 0.08 -50 −15 −5 5 15 θ (deg) Figure 7.3: Phase portraits of the Synthetic-Wheel Biped for different amplitudes of torso oscillation A the biped rolls on its stance foot from heel to toe and is able to take a step. The values of A for which the SWB is able to take a step depends on its inertial properties and its initial velocity. For the particular case of A = 0.0, we plot the torques τ2 and τ3 over one step in Figure 7.4; these torques are responsible for generating the trajectories of φ and ψ in Eqs.(7.26) and (7.27). The mechanical work done at these joints over one step can be expressed as an τ2 , τ3 (Nm) 5.0 θi = −11.25 0.0 τ3 τ2 τ2 τ3 -5.0 -15 θf = 11.25 -5 θ (deg) 5 15 Figure 7.4: Torque plots for an energy-conserving gait for the Synthetic-Wheel Biped (SWB). 136 θ θ τ2 , τ3 (Nm) 5.0 θf = 11.25 θ τ3 θ τ2 0.0 θi = −11.25 -5.0 -15 -5 5 15 θ (deg) Figure 7.5: Generalized torque plots for an energy-conserving gait for the Synthetic-Wheel Biped (SWB). integral over θ, as follows dφ ) dθ dθ t θ dψ ˙ W3 = τ3 ψ dt = τ3 ( ) dθ dθ t θ W2 = ˙ τ2 φ dt = τ2 ( θ θ θ τ2 dθ (7.28) θ τ3 dθ (7.29) θ θ The plots of τ2 and τ3 are shown in Figure 7.5. The area under these curves are equal to zero and this confirms that the gait described by Eqs.(7.26) and (7.27) conserves mechanical energy. 7.3.2 Point-Foot Three-Link Biped The point-foot three-link biped is obtained by replacing the arc-shaped feet on the SWB with point feet. A schematic of this biped is shown in Figure 7.6. Using the same set of generalized coordinates and generalized forces as the SWB, we first look at the sufficient condition for energy-conserving gait that corresponds to foot-ground interaction without 137 dt lt torso lsw x y swing foot dsw φ ψ g lst θ stance foot dst Figure 7.6: A schematic of the Point-Foot Three-Link. impact. This condition, given by Eq.(7.4), can be written as ˙ ˙ ˙ lst cos θf θf −lsw cos (θf +ψ f )θf −lsw cos (θf +ψ f )ψ f = 0 ˙ ˙ ˙ lst sin θf θf −lsw sin (θf +ψ f ) θf −lsw sin (θf +ψ f ) ψ f = 0 (7.30) The above equation was obtained by taking the derivative of the swing foot position coordinates relative to the stance foot - see Figure 7.6. Applying condition Eq.(7.2) next and assuming that the stance leg angle and swing leg angle satisfy (similar to Eqs.(7.24) and (7.25) for the SWB) θf = β/2, θ + ψ = −β/2 138 for some angle β, we get from Eq.(7.30) ˙ ˙ ˙ cos(β/2) lst θf −lsw θf −lsw ψ f = 0 ˙ ˙ ˙ sin(β/2) lst θf +lsw θf +lsw ψ f = 0 (7.31) ˙ ˙ Since β = 0, Eq.(7.31) implies θf = 0, which implies θi = 0 from Eq.(7.3). At the beginning and end of each step, the stance and swing legs will be symmetric with respect to the vertical and have zero velocity; and therefore, the biped will roll on its stance foot only if the torso is inclined with respect to the vertical. The sufficient condition in Eq.(7.2) dictates that the torso angle with respect to the vertical at the end of the step be equal and opposite to that at the beginning of the step. This implies that the biped will only roll back and forth on its stance leg. Clearly, it is not possible to design an energy-conserving gait for the point-foot three-link biped using the sufficient conditions in Eqs.(7.2), (7.3) and (7.4). 7.4 7.4.1 Five-DOF Point-Foot Biped System Description Consider the general five-dof point-foot biped shown in Figure 7.7. The generalized coordinates and generalized forces of the biped are as follows T q= θ1 θ2 θ3 θ4 θ5 T T = u1 u2 u3 u4 u5 (7.32) Assuming that the j-th link is driven by a motor mounted on the (j −1)-th link, j = 2, 3, 4, 5, 139 l3 θ3 torso l4 θ4 l2 θ5 l5 θ2 θ1 y l1 swing foot stance foot x Figure 7.7: A schematic of the five-dof point-foot biped. the generalized force vector can be rewritten in the form T T = −τ2 (τ3 − τ2 ) (τ4 − τ3 ) (τ5 − τ4 ) τ5 (7.33) where τj , j = 2, 3, 4, 5, are the corresponding motor torques. Note that τ1 = 0 since the ankle joint of the biped is passive. The equations of motion of the biped and description of its parameters are provided in Appendix C. 7.4.2 Sufficient Conditions for Energy-Conserving Gait The boundary conditions for the gait are determined as follows: 1. The potential energies at the beginning and end of each step are identical. This is achieved by f i θj = −θj , 140 j = 1, · · · , 5 (7.34) 2. The kinetic energies at the beginning and end of each step are identical. This is achieved by ˙i ˙f θj = θj , j = 1, · · · , 5 (7.35) 3. At the time of swing foot touchdown, there will be no impact due to impulsive forces f f if Vn,x = Vn,y = 0, i.e. f ˙f f ˙f f ˙f f ˙f l1 cos θ1 θ1 +l2 cos θ2 θ2 +l4 cos θ4 θ4 +l5 cos θ5 θ5 = 0 f ˙f f ˙f f ˙f f ˙f l1 sin θ1 θ1 +l2 sin θ2 θ2 + l4 sin θ4 θ4 +l5 sin θ5 θ5 = 0 (7.36) The above equation was obtained by taking the derivative of the swing foot position coordinates relative to the stance foot - see Figure 7.7. 7.4.3 Determination of Actuated Joint Trajectories Using Eq.(7.13), we consider the following functional dependence for θj o θj (θ1 ) = (aj θ1 + bj ) + fj (θ1 ), j = 2, 3, 4, 5 (7.37) ˙i ˙f ¨ The condition θ1 = θ1 in Eq.(7.35) can be satisfied if θ1 is an odd function of θ1 . To this ¨ end, we derive the expression for θ1 as follows N ¨ θ1 = D 141 (7.38) where N = N(q, q, q ) and D = D(q) are defined as follows (see Appendix C): ˙ ¨ 5 N(q, q, q ) = − ˙ ¨ 5 j=2 i=1 5 5 ¨ M(j, i)θj − 5 j=1 i=1 5 ˙ N(i, j)θj − G(j) j=1 M(j, 1) D(q) = j=1 ¨ Using Eq.(7.37), a term-by-term analysis of N and D shows that θ1 is an odd function of θ1 for any values of aj ’s provided that bj ’s are chosen as follows: bj = kj π, kj = 0, ±1, ±2, · · · (7.39) o For a realistic walking configuration, we set b2 = b3 = 0 and b4 = b5 = π. Assuming fj in Eq.(7.37) to be sinusoidal functions, we have θ2 (θ1 ) = a2 θ1 + A sin (Bθ1 ) θ3 (θ1 ) = a3 θ1 + C sin (Dθ1 ) θ4 (θ1 ) = a4 θ1 + π + E sin (F θ1 ) θ5 (θ1 ) = a5 θ1 + π + G sin (Hθ1 ) (7.40) f i We assume θ1 = −θ1 = β/2 for some angle β and choose B = D = F = 2π/β; this choice ensures that actuated joint trajectories of θ2 , θ3 and θ4 include one period of sinusoidal oscillation in each step. For a reasonable clearance between the swing foot and the ground, we choose H differently, equal to 2.5π/β. To have identical configuration of the legs at the beginning and end of each step, links 1 and 5 and links 2 and 4 should have the same length and be symmetrical with respect to 142 the vertical at the beginning of each step. From the definition of the joint angles in Figure 7.7, this implies that i i θ5 = −θ1 + π (7.41) i i θ4 = −θ2 + π (7.42) Furthermore, the biped should have the same speed before and after interchange of the stance and swing legs. This implies ˙f ˙i θ5 = θ1 (7.43) We use Eq.(7.42) to obtain a4 = −a2 and use Eqs.(7.41) and (7.43) to solve for a5 and G in terms of β. The actuated joint trajectories in Eq.(7.40) can now be described with the help of five parameters, namely, a2 , a3 , A, C and E, as shown below: θ2 (θ1 ) = a2 θ1 + A sin(2πθ1 /β) θ3 (θ1 ) = a3 θ1 + C sin(2πθ1 /β) θ4 (θ1 ) = −a2 θ1 + π + E sin(2πθ1 /β) θ5 (θ1 ) = a5 θ1 + π + G sin(2.5πθ1 /β) (7.44) Of the five parameters, a2 , A and E are chosen such that the two boundary conditions in Eq.(7.36) are satisfied. This provides some flexibility in choice of these parameters. The remaining two parameters, a3 and C, are related to the motion of the torso and therefore do not show up in the no-impact boundary conditions in Eq.(7.36); they can be chosen arbitrarily and provide additional flexibility in design of energy-conserving gaits. 143 Table 7.2: Kinematic and dynamic parameters of Five-DOF Point-Foot biped with similar leg links Kinematic and dynamic parameters j lj (m) dj (m) mj (kg) Ij (kg.m2 ) 1, 5 0.600 0.300 2.000 0.060 2, 4 0.400 0.200 2.000 0.027 3 0.600 0.300 4.000 0.120 7.4.4 Numerical Simulations The kinematic and dynamic parameters of the biped with similar leg links (l1 = l5 , l2 = l4 ) are listed in Table 7.2. Although there is significant flexibility in choosing the gait parameters, many choices will result in gaits that are not feasible. For example, certain choices of parameters will result in the swing leg passing through the ground. A set of parameters that result in a feasible gait can be obtained by trial and error; one such set of f θ1 = −22.5 τ2 , τ3 , τ4 , τ5 (Nm) 40.0 20.0 τ3 τ2 τ4 0.0 τ5 -20.0 τ2 τ5 τ4 i θ1 = 22.5 τ3 -40.0 -20 -10 0 10 20 θ1 (deg) Figure 7.8: Torque plots for an energy-conserving gait for the Five-DOF Point-Foot Biped with similar leg links. 144 parameters for a step angle of β = 45.0 deg is given below: a5 = −1.683 G = −0.739 (7.45) a2 = 1.200, a3 = 0.000 A = 0.308, C = −0.550, E = 0.008 ˙i For an initial velocity of θ1 = −100 deg/s, the torques τj , j = 2, 3, 4, 5, are shown in Figure 7.8. The total mechanical work done by the actuators over one step can be expressed as an integral over θ1 , as follows 5 W = j=2 t ˙ τj θj dt = f θ1 i θ1   5 j=2  dθj τj ( ) dθ1 dθ1 f θ1 i θ1 τ θ dθ1 (7.46) The plot of τ θ is shown in Figure 7.9. The area under this curve is equal to zero and this confirms that the gait described by Eqs.(7.44) and (7.45) conserves mechanical energy. The configurations of the biped at different instants of time over a single step, moving from left to f θ1 = −22.5 τ θ (Nm) 100.0 0.0 -100.0 i θ1 = 22.5 -20 -10 0 10 20 θ1 (deg) Figure 7.9: Generalized torque plot for an energy-conserving gait for the Five-DOF PointFoot Biped with similar leg links. 145 Figure 7.10: Configurations of the five-dof point-foot biped with similar leg links at different instants of time over one step. right, are shown in Figure 7.10. It was verified that the swing leg does not scuff the ground. Since actuated joint trajectories of the biped, θj (t), j = 2, 3, 4, 5, are odd functions of the passive coordinate θ1 , a right-to-left step of the biped will have the same set of configurations as in Figure 7.10. 7.5 Two-Step Periodic Gaits: Five-DOF Bipeds 7.5.1 Generalization to Two-Step Periodic Gait The energy-conserving gait designed in the last section is one-step periodic. This is because of Eqs.(7.41), (7.42) and (7.43), which impose additional conditions that make the configurations4 of the five-dof point-foot biped identical at the beginning and end of each step. It is not necessary to have identical configurations at the beginning and end of each step to have an energy-conserving gait, and therefore, energy-conserving gaits need not be one-step periodic. Energy-conserving gaits designed using the sufficient conditions in section 7.2.2 will however be two-step periodic, i.e., result in identical configurations of the biped after two steps. This is because two times application of the boundary conditions in Eqs.(7.2) and (7.3) will revert each link back to its original configuration. 4A biped configuration is described by the position and velocity of its links. 146 7.5.2 Characterization of Trajectories For a two-step periodic gait, the configurations of the biped will be different at the beginning and the end of the first step. Consequently, we will have a different initial configuration for the second step after the interchange of the stance and swing legs. This will result in a new set of parameters for the actuated joint trajectories for the second step. In this section, we modify some conditions in Section 7.4.3 to obtain the trajectories of the actuated joints for each step. In spite of different walking gaits for the first and second steps, we still ensure that the energy is conserved for each step by satisfying the sufficient conditions in Eqs.(7.34), (7.35) and (7.36). 7.5.2.1 Trajectories for the First Step Consider the trajectories in Eq.(7.40) with the parameters B, D, F and H chosen the same as those in Section 7.4.3, i.e., θ2 (θ1 ) = a2 θ1 + A sin(2πθ1 /β) θ3 (θ1 ) = a3 θ1 + C sin(2πθ1 /β) θ4 (θ1 ) = a4 θ1 + π + E sin(2πθ1 /β) θ5 (θ1 ) = a5 θ1 + π + G sin(2.5πθ1 /β) (7.47) Since the condition in Eq.(7.42) is not necessary for a two-step periodic gait, the parameters a2 and a4 are now chosen independently. The condition in Eq.(7.41) is replaced with the following condition which ensures that the swing foot is on the ground at the beginning of the step: i i i i l1 cos θ1 + l2 cos θ2 + l4 cos θ4 + l5 cos θ5 = 0 147 (7.48) In addition, for a two-step periodic gait, we can have different initial velocities of the stance leg in the first and second steps. Therefore, the condition in Eq.(7.43) is modified to: ˙f ˙i θ5 = cv θ1 , cv > 0 (7.49) f i Knowing that θ1 = −θ1 = β/2, Eqs.(7.48) and (7.49) can be solved for a5 and G in terms of β. The parameters A and E are chosen such that the two boundary conditions in Eq.(7.36) are satisfied. The remaining two parameters, a3 and C, which are related to the motion of the torso, can be chosen arbitrarily and provide some flexibility in design of energy-conserving gaits. 7.5.2.2 Trajectories for the Second Step The actuated joint trajectories for the second step should be chosen such that the sufficient boundary conditions for energy conservation are satisfied, and the configuration of the biped at the beginning of the second step (after interchange of the legs) matches the configuration at the end of the first step (before interchange of the legs). This is a direct consequence of having no impact at the instant of swing-foot touchdown at the end of the first step. To this end, the actuated joint trajectories in Eq.(7.47) are modified as follows: ˆ ˆ θ2 (θ1 ) = a2 θ1 + A sin(2πθ1 /β) ˆ ˆ ˆ θ3 (θ1 ) = a3 θ1 + C sin(2πθ1 /β) ˆ ˆ ˆ θ4 (θ1 ) = a4 θ1 + π + E sin(2πθ1 /β) ˆ ˆ ˆ θ5 (θ1 ) = a5 θ1 + π + G sin(2.5πθ1 /β) ˆ 148 (7.50) ˆ ˆ ˆ ˆ ˆ ˆ where aj ’s, A, C, E, G and β are new parameters to be determined. The parameter β ˆ represents the total angle traveled by the new stance leg and is obtained as: i(2) f (1) ˆ β = 2θ1 = 2θ5 = −2[a5 (β/2) + G sin(2.5π/2)] (7.51) where the superscripts (1) and (2) denote the corresponding variables for the first and second step, respectively. To ensure that the position and velocity of the torso is the same before and after the interchange of the legs, we have: i(2) f (1) ⇒ ˆ a3 = −(β/β) a3 ˆ ˙f (1) ˙i(2) θ3 = θ3 ⇒ ˆ β 1 2π ˆ C= [ˆ3 − (a3 − a C)] 2π cv β θ3 = θ3 (7.52) ˙i(1) ˙i(2) ˆ where we used the relation θ1 = cv θ1 from Eq.(7.49) in deriving C. To ensure that the positions of the upper links of the stance and swing legs are the same before and after the interchange, we have: i(2) θ2 i(2) θ4 f (1) ⇒ ˆ a2 = −(β/β) a4 ˆ f (1) ⇒ ˆ a4 = −(β/β) a2 ˆ = θ4 = θ2 (7.53) To ensure continuity in the position and velocity of the lower link of the swing leg, we impose the conditions: i(2) θ5 f (1) = θ1 ˙i(2) ˙f (1) θ5 = θ1 149 (7.54) ˆ The above equations can be solved for the parameters a5 and G. The remaining two paramˆ ˆ ˆ eters A and E in Eq.(7.50) can be obtained by solving the no-impact boundary conditions ˆ ˆ in Eq.(7.36) for the second step. Note that these choices of A and E will guarantee identical velocities for the links 2 and 4 before and after the interchange. Remark 7.5.1. For the two-step periodic gait, the actuated joint trajectories for the first step in Eq.(7.47) and the second step in Eq.(7.50) will be sequentially repeated. 7.5.3 Five-DOF Point-Foot Biped: Similar Leg Links Consider the five-dof point-foot biped with similar leg links, presented in section 7.4 with the same kinematic and dynamic parameters listed in Table 7.2. To design an energy-conserving gait with two-step periodicity, we remove the conditions in Eq.(7.41) and Eq.(7.42) to make the initial configuration of the legs asymmetric with respect to the vertical. The procedure in Section 7.5.2 is followed to determine gait parameters for the first and second steps. A feasible gait is obtained by using the following parameters for the first step, described by Eq.(7.47): β = π/4 a3 = −0.300 C = −0.500 a2 = 1.200, a4 = −1.500 a5 = −1.267 G = −0.321, A = 0.261, 150 E = −0.059 (7.55) θ1 (deg) first step τ2 , τ3 , τ4 , τ5 (Nm) 80.0 f θ1 -10 0 f θ1 = −22.5 10 = −15.5 40.0 τ4 0.0 τ5 τ5 τ2 -40.0 -80.0 τ3 τ3 i θ1 = 22.5 -20 -10 0 10 τ2 τ4 i θ1 = 15.5 20 second step θ1 (deg) Figure 7.11: Torque plots for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with similar leg links. and the following parameters for the second step, described by Eq.(7.50): ˆ β = 0.172π a3 = 0.435 ˆ a2 = 2.176, ˆ a4 = −1.741 a5 = −2.288 ˆ ˆ ˆ C = −0.293 ˆ ˆ G = −0.321, A = 0.276, (7.56) ˆ E = −0.073 θ1 (deg) τ θ (Nm) 150.0 first step -10 0 10 0.0 -150.0 i θ1 = 22.5 -20 -10 0 10 i θ1 = 15.5 20 second step θ1 (deg) Figure 7.12: Generalized torque plot for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with similar leg links. 151 first step second step Figure 7.13: Configurations of the five-dof point-foot biped with similar leg links at different instants of time for a two-step periodic gait. The joint configurations at the beginning of the first step and end of the second step are the same but they are different from the joint configuration at the end of the first step. ˆ The value of cv in Eq.(7.49) was set to unity. Note that β < β, which implies that the second ˙ i(2) ˙i(1) step is shorter than the first step. For initial velocities of θ1 = θ1 = −100 deg/s, the torques τj , j = 2, 3, 4, 5, are shown for the first and second steps in Figure 7.11. The total mechanical work done by the actuators over each step can be expressed by Eq.(7.46). The plot of τ θ for both steps is shown in Figure 7.12. The area under this curve for each step is equal to zero and this confirms that the gait conserves mechanical energy in each step. The configurations of the biped at different instants of time over two steps, moving from left to right, is shown in Figure 7.13. It was verified that the swing leg does not scuff the ground. 7.5.4 Five-DOF Point-Foot Biped: Dissimilar Leg Links Consider the biped in Figure 7.7 and assume the upper link of each leg to have the same length as the lower link of the other leg, i.e. l1 = l4 and l2 = l5 . The kinematic and dynamic Table 7.3: Kinematic and dynamic parameters of Five-DOF Point-Foot biped with dissimilar leg links j (1st step) 1, 4 2, 5 3 Kinematic and dynamic parameters j (2nd step) lj (m) dj (m) mj (kg) Ij (kg.m2 ) 2, 5 0.600 0.300 2.000 0.060 1, 4 0.400 0.200 2.000 0.027 3 0.600 0.300 4.000 0.120 152 θ1 (deg) first step τ2 , τ3 , τ4 , τ5 (Nm) 80.0 0 20 f f θ1 = −31.77 τ5 τ5 τ4 θ1 = −22.5 40.0 0.0 -40.0 -20 τ4 τ2 -80.0 -20 τ2 τ3 τ3 i θ1 = 31.77 i θ1 = 22.5 0 20 θ1 (deg) second step Figure 7.14: Torque plots for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with dissimilar leg links. parameters of the biped are listed in Table 7.3. Since the leg links are not similar, the kinematic and dynamic parameters of the links are interchanged for the first and second steps - see Table 7.3. To generate a two-step periodic gait, we choose an asymmetric configuration of the legs at the beginning of the first step. Following the procedure in Section 7.5.2, a feasible gait is obtained by using the following parameters for the first step, described by θ1 (deg) first step 200.0 -20 τ θ (Nm) 20 f θ1 = −31.77 f 100.0 0 θ1 = −22.5 0.0 -100.0 -200.0 i θ1 = 22.5 -20 0 20 θ1 (deg) i θ1 = 31.77 second step Figure 7.15: Generalized torque plot for an energy-conserving gait with two-step periodicity for the Five-DOF Point-Foot Biped with dissimilar leg links. 153 first step second step Figure 7.16: Configurations of the five-dof point-foot biped with dissimilar leg links at different instants of time for a two-step periodic gait. The joint configurations at the beginning of the first step and end of the second step are the same but they are different from the joint configuration at the end of the first step. Eq.(7.47): β = π/4 a3 = −0.300 C = −0.550 a2 = 1.200, a4 = −0.800 a5 = −2.236 G = −0.458, A = 0.364, (7.57) E = −0.013 and the following parameters for the second step, described by Eq.(7.50): ˆ β = 0.353π a3 = 0.213 ˆ a2 = 0.567, ˆ a4 = −0.850 a5 = −1.292 ˆ ˆ ˆ ˆ G = −0.458, A = 0.223, ˆ C = −0.671 (7.58) ˆ E = 0.152 ˆ The value of cv in Eq.(7.49) was set to unity. Note that β > β, which implies that the second ˙i(1) ˙i(2) step is longer than the first step. For an initial velocity of θ1 = θ1 = −130 deg/s, the torques τj , j = 2, 3, 4, 5, are shown for the first and second steps in Figure 7.14. The total mechanical work done by the actuators over one step can be expressed by Eq.(7.46). The plot of τ θ for both steps is shown in Figure 7.15. The area under this curve for each step is equal to zero and this confirms that the gait conserves mechanical energy in each step. The configurations of the biped at different instants of time over two steps, moving from left to right, are shown in Figure 7.16. It was verified that the swing leg does not scuff the ground. 154 7.6 Conclusion We proposed a general method for designing gaits for point-foot passive-ankle planar bipeds that require zero mechanical energy over each individual step. The proposed method provides trajectories for all actuated dof of the biped which can be tracked using classical control techniques. These trajectories ensure that the biped has identical potential energy and kinetic energy at the beginning and end of each step and there is no loss of energy due to impact at the time of swing-foot touchdown. To show its effectiveness, the proposed method has been applied to several three-dof and five-dof bipeds; the designed gaits are one-step periodic in some cases and two-step periodic in others. The method can be easily applied to bipeds with greater number of dof’s. The proposed method provides significant flexibility in choosing joint trajectories for the actuated dof of the biped. This flexibility can be used to design a wide range of gaits that have different step lengths, walking speeds, and range of joint torques, for example, and can therefore be used for optimizing the energy-conserving gaits. For implementation, controllers will have to be designed to track the desired trajectories accurately. Any deviation from the desired trajectories will result in a gait that is not energy-conserving but continuous dependence of the mechanical energy on solution trajectories dictates that small deviations in the trajectories will only require a small amount of mechanical energy. A small deviation in the desired trajectories can be due to inaccurate tracking or may be introduced purposely to provide uniformly larger ground clearance while walking on rough terrain, for example. An energy-conserving gait requires zero mechanical energy but will have zero energy consumption if the negative work done by the actuators can be recaptured and stored. This will require particular robot capabilities, such as backdriveable actuators and electrical circuits 155 capable of energy regeneration and storage. The boundary conditions presented here provide a simple framework for designing energyconserving gaits but since these conditions are sufficient there may exist other gaits that require zero mechanical energy. Of the three conditions presented, the no-impact boundary condition is necessary to avoid energy dissipation but the other two conditions relating to kinetic and potential energies of the links can be relaxed. In particular, it may be possible to find energy-conserving gaits by requiring the total kinetic energy and potential energy of the biped to be the same at the beginning and end of each step or even multiple steps. 156 Chapter 8 Intermittent Output Tracking For Linear Non-Minimum-Phase Systems 8.1 Introduction Underactuated systems with unstable zero dynamics are Non-Minimum-Phase (NMP) systems. The limitations in output tracking for NMP continuous-time systems motivated us to look for a controller that ensures output matching with the reference trajectory at regular intervals of time. The control objective is similar to the output tracking controllers for discrete-time systems. Although perfect output tracking can be achieved for discretetime systems, we propose a continuous-time controller to intermittently match the system output with the reference trajectory for two reasons: First, in the continuous-time control design, we do not need to incorporate the dynamics of holders and samplers associated with discrete-time sampled systems. Second, the discretized systems require strong conditions to be satisfied to preserve the stability of the sampling zeros [139, 140, 141]. For SISO linear systems with relative degree greater than two, Astrom et al. [142] showed that at least one of the the limiting zeros (the sampling zeros for the limiting case of zero sampling period) is outside the unit circle, irrespective of their location for the original continuous system. This implies that even for a MP continuous-time system, we might face control limitations induced by unstable zeros for the discretized system. 157 In this chapter, we propose an output tracking controller for linear Single-Input SingleOutput (SISO) NMP systems. A continuous-time controller is designed to take the output tracking error to zero at regular intervals of time. The controller for the NMP plant is based on switched inputs, which are derived from the continuous control inputs of a “Discrete Equivalent” (DE) [143] MP system. The control design ensures zero output tracking error at regular intervals of time and guarantees the stability of the internal dynamics. A finite preview of the desired trajectory is required for this method. 8.2 Discrete Equivalent Systems Consider a linear SISO system represented in the following state-space form: x(t) = Ax(t) + Bu(t) , x(t0 ) = x0 ˙ y(t) = Cx(t) (8.1) Let the system input u, shown in Figure 8.1, be defined as: u(t) =   M  0 if (j −1)δ < t − t0 < (j −1+p)δ if j = 1, 2, ... (8.2) (j −1+p)δ < t − t0 < jδ The switching parameter δ > 0 represents the period of the switched input and p ∈ (0, 1) defines the fraction of the time interval when the input has a nonzero constant value M. Definition 8.2.1. [143] Consider the SISO linear time-invariant (LTI) system with constant input u = M: ¯ 158 PSfrag u, u ¯ δ δ M pδ z t0 pδ t1 t2 t Figure 8.1: The periodic switched input for the system in Eq.(8.1) (in solid line) and the constant input for DE system in Eq.(8.3) (in dashed line) ¯x ¯¯ ¯ ˙ x(t) = A¯(t) + B u , x(t0 ) = x0 ¯ ¯¯ y (t) = C x(t) ¯ (8.3) This system is a Discrete Equivalent (DE) of the system in Eq.(8.1) with switched input in Eq.(8.2) if the state variables and the output of two systems match exactly at t = tj , where tj = t0 + jδ and j = 1, 2, · · · , after starting from the same initial conditions at t = t0 . To obtain the relations between switched-input system and its DE, we first find the state variables of the two systems at times t = t0 + jδ, j = 1, 2, · · · . For the switched-input system in Eq.(8.1), we can write: x(t1 ) = eA(t1 −t0 ) x0 + t1 t0 eA(t1 −τ ) Bu(τ ) dτ = eApδ x0 + A−1 (eApδ − I)BM x(t2 ) = eA(t2 −t1 ) x(t1 ) = eAδ x0 + A−1 (eAδ − eA(1−p)δ )BM 159 (8.4) where I represents the identity matrix. For the DE system in Eq.(8.3), we have: ¯ x(t2 ) = eA(t2 −t0 ) x0 + ¯ t2 t0 ¯ ¯ ¯ ¯ ¯ ¯ eA(t2 −τ ) BM dτ = eAδ x0 + A−1 (eAδ − I)BM (8.5) From Definition 8.2.1: x(t2 ) = x(t2 ), y(t2 ) = y (t2 ) ¯ ¯ (8.6) Using Eqs.(8.4) and (8.5), we get: ¯ A =A ¯ B =(eAδ − I)−1 (eAδ − eA(1−p)δ )B ¯ C =C (8.7) We demonstrate the concept of DE systems with the example below. Example 8.2.1. Consider a linear system with the following transfer function: G(s) = (s + 1)(s − 2) (s + 1)(s + 2)(s + 3) (8.8) It has the following state-space representation:    0   −1 0  1       A =  1 −5 −2.45  , B =  0       0 2.45 0 0 160      , C = [1 − 6 − 3.27]   (8.9) If we assume: M = 2 , δ = 0.8 , p = 0.2 (8.10) the DE system will have the following state-space representation:    0   −1 0  0.14     ¯  ¯  A =  1 −5 −2.45  , B =  0.028       0 2.45 0 0.015     ¯  , C = [1 − 6 − 3.27]   (8.11) The transfer function of the DE system can be shown to be: (s + 1)(s + 7.137) ¯ ¯ ¯ ¯ G(s) = C(sI − A)−1 B = (s + 1)(s + 2)(s + 3) a x1 , x1 ¯ 0.5 x2 , x2 ¯ 0 0.04 x3 , x3 ¯ (8.12) b 0 0.1 c 0 2 u, u ¯ d 0 0 0.8 1.6 2.4 3.2 4 time (sec) Figure 8.2: Simulation results for example 8.2.1 showing state variables and inputs for switched-input system (solid line) and discrete equivalent system (dashed line). 161 A simulation of the response of the switched-input system and its DE is shown in Figure 8.2. It can be seen that all state variables of the two systems match exactly at regular intervals of time, δ = 0.8 for u = 2. 8.3 Zero Shift for Linear NMP Systems The poles of the DE system are identical to the poles of the switched-input system since ¯ A = A in Eq.(8.7). The location of the zeros of the systems can however be different as shown in Example 8.2.1. In this example, the original system has a RHP zero at s = 2. This zero is moved to the LHP at s = −7.137 for the DE system. This indicates that a NMP system can have a DE system which is MP. 8.3.1 Stable NMP Systems with Distinct Poles - General Case Consider the stable linear system with distinct poles λi , i = 1, · · · , n: b sn−1 + b2 sn−2 + · · · + bn−1 s + bn Y (s) = 1 U(s) (s − λ1 )(s − λ2 ) · · · (s − λn ) = c1 c2 cn + +···+ s − λ1 s − λ2 s − λn (8.13) where bi ’s are the constant coefficients of the zero polynomial and ci ’s are constant residues of partial fraction expansion of the system transfer function. Assume that the system has at least one zero in RHP, and therefore it is NMP. The transfer function in Eq.(8.13) can be 162 represented in Diagonal Canonical form [144] as: x =Ad x + Bd u ˙ y =Cd x (8.14) where: T Ad =diag λ1 λ2 · · · λn , Bd = 1 1 · · · 1 C d = c1 c2 · · · cn (8.15) For the corresponding DE system, we have from Eq.(8.7): ¯ ¯ Ad =Ad , Cd = Cd ¯ Bd =(eAd δ − I)−1 (eAd δ − eAd (1−p)δ )Bd = h1 h2 · · · hn T (8.16) where hi ’s are some positive constants obtained using the properties of diagonal matrix Ad : hi = eλi δ − eλi (1−p)δ eλi δ − 1 , i = 1, . . . , n (8.17) The transfer function of the DE system can be evaluated as: c h cn hn c h ¯ ¯ ¯ ¯ G(s) =Cd (sI − Ad )−1 Bd = 1 1 + 2 2 + · · · + s − λ1 s − λ2 s − λn 163 (8.18) The zero polynomial of the DE system (with distinct poles) is obtained as: ¯ Z(s) = c1 h1 (s − λ2 )(s − λ3 ) · · · (s − λn ) + c2 h2 (s − λ1 )(s − λ3 ) · · · (s − λn ) +··· + cn hn (s − λ1 )(s − λ2 ) · · · (s − λn−1 ) (8.19) which can be expanded into the following form: ¯ Z(s) = ¯1 sn−1 + ¯2 sn−2 + ¯3 sn−3 + · · · + ¯n−1 s + ¯n b b b b b (8.20) where n ¯1 = b ci hi i=1 ¯m =(−1)m−1 b n n λj1 λj2 · · · λjm−1 )ci hi , ( i=1 m = 2, · · · , n (8.21) j1 , j2 , · · · , jm−1 = i j1 < j2 < · · · < jm−1 Remark 8.3.1. For the NMP system in Eq.(8.13) and switching parameters δ and p, the DE system will be MP if the zero polynomial in Eq.(8.20) does not have any root in the closed RHP. This can be verified by applying Routh’s stability criterion [144] without actually solving for the zero locations. Remark 8.3.2. The order of the zero polynomial of the original NMP system and its DE may be different. This implies that the number of zeros for the DE system may be different from 164 that of the original system. The number of zeros for both systems is however less than n. 8.3.2 Stable NMP Systems with Distinct Poles - Special Case In this section, we investigate a special case of the system in Eq.(8.13) with n = 2 and stable real modes, i.e., b1 s + b2 c1 c2 Y (s) = = + U(s) (s − λ1 )(s − λ2 ) (s − λ1 ) (s − λ2 ) (8.22) Using Eqs.(8.20) and (8.21), the zero polynomial for the DE system is obtained as follows: ¯ b Z = ¯1 s + ¯2 b (8.23) where ¯1 = c1 h1 + c2 h2 , b ¯2 = −(c1 h1 λ2 + c2 h2 λ1 ) b (8.24) and h1 and h2 can be obtained from Eq.(8.17) Lemma 8.3.1. If λ2 < λ1 < 0, then h2 < h1 . λδ λ(1−p)δ is strictly increasing, i.e., dh > 0 in Proof. We prove by showing that h(λ) = e −e dλ eλδ −1 λ ∈ (−∞, 0). To this end, we define a new variable q: q = eλδ − 1 From λδ < 0, we know that −1 < q < 0. The function h can be rewritten in terms of the 165 new variable q as: h(q) = q + 1 − (q + 1)(1−n) q Therefore, ∂h ∂q 1 + pq − (q + 1)p dh ) (δeλδ ) =( )( )=( dλ ∂q ∂λ q 2 (q + 1)p Since δeλδ > 0 and q 2 (q + 1)p > 0, the sign of (dh/dλ) is identical to the sign of (1 + pq − (q + 1)p ). Since |q| < 1, (q + 1)p can be expanded using the binomial series [145]. Therefore, 1+pq−(1+q)p =1 + pq − ∞   p(p−1) 2 p(p−1)(p−2) 3  p  j q + q + ···)   q = −( 2 3! j=0 j which is always positive for 0 < p < 1 and q < 0. This completes the proof. The following theorem states the sufficient conditions for the DE system of Eq.(8.22) to be MP. Theorem 8.3.1. Consider the second-order linear system in Eq.(8.22), which has two negative real and distinct poles λ1 and λ2 and a RHP zero. Without loss of generality, let λ1 > λ2 . The DE system will be MP if any of the following conditions is satisfied: (a) c1 c2 > 0 (b) c1 c2 < 0 and |c1 | ≥ |c2 | (c) c1 c2 < 0, |c1 | ≤ |c2 | and c1¯1 > 0 b Proof. Since h1 , h2 > 0 and λ1 , λ2 < 0, condition (a) implies that ¯1 and ¯2 in Eq.(8.24) will b b have the same sign, which implies that the DE system is MP. 166 For condition (b), first consider the case that c1 > 0 and c1 ≥ −c2 . Since λ1 > λ2 , Lemma 8.3.1 implies: h1 > h2 ⇒ c1 h1 > −c2 h2 ⇒ ¯1 > 0 , b c1 h1 > −c2 h2 ⇒ c1 h1 λ2 < −c2 h2 λ2 < −c2 h2 λ1 ⇒ ¯2 > 0 b If c1 < 0 and c1 ≤ −c2 , it can be shown that ¯1 , ¯2 < 0. This proves that (b) is a sufficient b b condition for the DE system to be MP. For condition (c), first consider the case that c1 < 0, c2 > −c1 and ¯1 < 0. Since λ1 > λ2 , b Lemma 8.3.1 gives: h1 > h2 ⇒ λ1 h2 > λ1 h1 > λ2 h1 ⇒ λ1 h2 c2 > λ2 h1 c2 > −λ1 h1 c1 ⇒ ¯2 < 0 b For the case that c1 > 0, c1 < −c2 and ¯1 > 0, we can show that ¯2 > 0. Therefore, b b condition (c) will ensure that the DE system is MP. 8.3.3 Stable NMP Systems with Repeated Poles For a linear system with repeated poles, the state-space matrices can be written in the Jordan Canonical Form [144]. Deriving a general expression for the zero polynomial of the DE system is not as simple as the case of distinct poles. This is because of the complexity ¯ in calculating B in Eq.(8.7) due to non-diagonal matrix A. The DE system can nevertheless be designed to be MP even when the original system has RHP zeros - this is illustrated with the help of numerical examples in Table 8.1. Using a numerical search for proper switching parameters δ and p, it is shown that the DE system can be MP when the original NMP 167 Table 8.1: Linear NMP systems with repeated poles and their corresponding MP DE systems Original Zeros Poles δ p DE zeros 5 [-2,-2,-5,-5] 0.95 0.15 [-33.61,-6.72±1.84i] [-1,5] [-2,-2,-2] 0.85 0.3 [-5.86,-1] [3,4] [-4,-4,-5,10] 1 0.1 [-13.49,-9.46±1.45i ] [-1,3±i] [-3,-3,-5±i] 0.5 0.75 [-0.22,-0.59,-29.39] system has repeated poles. It should be noted that Table 8.1 includes an example of a system with relative degree three. For such a system, sampling will necessarily result in a NMP discrete-time system [142]. Remark 8.3.3. In this section, we only considered stable systems. If the system is unstable, it can be first stabilized using a state-feedback controller provided that the system is stabilizable. The switching parameters can then be found for the stabilized system. This will be discussed further in the next section. 8.4 Intermittent Output Tracking Control It is not possible to design controllers for asymptotic output tracking of arbitrary reference trajectories for NMP systems. However, it may be possible to find a Discrete Equivalent Minimum-Phase (DEMP) system for a NMP system through proper choice of switching parameters δ and p. From Definition 8.2.1, the state variables and the output of the NMP system will match with those of DEMP system at regular intervals of time. Therefore, if the output tracking control problem is solved for the DEMP system, the output of the NMP system will also match with the reference trajectory at regular intervals of time δ. We define Intermittent Output Tracking next: Definition 8.4.1. For a linear SISO NMP system, Intermittent Output Tracking (IOT) 168 problem refers to the task of finding a control input u(t) such that for any desired reference trajectory r(t) and any initial condition x0 : a) The system is internally stable, i.e, all state variables are bounded for any bounded r(t). b) The output y(t) matches with r(t) at regular intervals of time δ, i.e., y(t0 + jδ) → r(t0 + jδ) as j → ∞ for j = 1, 2, · · · . To solve the IOT problem, we first review the tracking control design for MP systems. 8.4.1 Tracking Control for MP Systems Consider a MP system which is represented by the following transfer function: m m−1 + · · · + β m−1 s + βm ¯ N(s) = β1 s + β2 s G= n + a sn−1 + · · · + a D(s) s 1 n−1 s + an (8.25) where m < n and β1 = 0. To obtain a state-space representation of the system in the normal, we follow the procedure outlined in [146]. First, the denominator polynomial D(s) is written as: D(s) = Q(s)N(s) + R(s) (8.26) 1 ρ s + q1 sρ−1 + · · · + qρ−1 s + qρ β1 (8.27) where Q(s) = is the quotient polynomial of order ρ = (n − m), where ρ is the relative degree of the system, and R(s) is the remainder polynomial. In normal form, the state-space representation is 169 found to be: ¯ ˙ η =A0 η + B0 Cc ζ ¯ ¯ ˙ ¯ ¯ ¯ ¯ ¯ ζ =Ac ζ + Bc (γ T ζ − β1 C0 η + β1 u) ¯ y =Cc ζ ¯ (8.28) ¯ ˙ where [¯, ζ]T is the vector of state variables with η ∈ Rn−ρ , and ζ = [¯, y , · · · , y (ρ−1)]T ∈ Rρ η ¯ ¯ y ¯ ¯ and  0 1   0 0   . Ac =  . .   0 · · ·   0 ···    · · · 0 0       0  1 · · · 0      . , B = . , C = .. . . . c c . 1 0 ··· 0 0 .       0 0 1       0 0 1 0 γ T = − β1 qρ qρ−1 · · · q2 q1 (8.29) and (A0 , B0 , C0 ) is a minimal realization of transfer function R(s)/N(s). Note that the eigenvalues of A0 are identical to the zeros of the MP system in Eq.(8.25) and therefore A0 is Hurwitz. We assume that the reference trajectory r(t) and its derivatives up to r (ρ) (t) are available on-line, and r (ρ) (t) is piecewise continuous. The tracking error is defined as: e = ζ − Rd = y − r ¯ ¯ ¯ T ˙ y−r ¯ ˙ ··· y (ρ−1) − r (ρ−1) ¯ (8.30) Using the change of variables e = ζ − Rd , the state equations in Eq.(8.28) can be rewritten ¯ ¯ 170 as: ˙ η =A0 η + B0 Cc (¯ + Rd ) ¯ ¯ e ˙ e =Ac e + Bc [γ T (¯ + Rd ) − β1 C0 η + β1 u − r (ρ) ] ¯ ¯ e ¯ ¯ (8.31) The following state-feedback controller, u= ¯ 1 ¯¯ [−γ T (¯ + Rd ) + β1 C0 η + r (ρ) − K e] e ¯ β1 (8.32) ¯ ¯ with control gain K chosen to make (Ac − Bc K) Hurwitz, yields the closed-loop system: ˙ η =A0 η + B0 Cc (¯ + Rd ) ¯ ¯ e ¯ e ˙ e =(Ac − Bc K)¯ ¯ (8.33) ¯ Since (Ac − Bc K) is Hurwitz, e → 0 as t → ∞. The internal stability of the closed-loop ¯ system is also ensured since A0 is Hurwitz. 8.4.2 Intermittent Output Tracking for NMP Systems Consider the normal form of a MP system in Eq.(8.28) which can be rewritten as: ¯x ¯¯ ˙ x(t) = A¯(t) + B u(t) , x(t0 ) = x0 ¯ ¯ ¯¯ y (t) = C x(t) ¯ (8.34) 171 where x = [¯, ζ]T and ¯ η ¯     0  A0 B0 Cc  ¯ ¯  ¯  A =  , C = 0 Cc , B =  β1 Bc −β1 Bc C0 Ac + Bc γ T (8.35) ¯ ¯ ¯ Now, assume that there are switching parameters δ and p such that (A, B, C) in Eq.(8.35) represents a DE system for the following NMP system: x(t) =Ax(t) + Bu(t) , x(t0 ) = x0 ˙ y(t) =Cx(t) (8.36) where the state-space matrices are obtained using Eq.(8.7) as: ¯ ¯ A =A , C = C ¯ B =(eAδ − eA(1−p)δ )−1 (eAδ − I)B (8.37) From Definition 8.2.1, state variables and outputs of the systems in Eqs.(8.34) and (8.36) will be discretely identical provided that over each interval δ, u(t) is constant and u(t) is ¯ switched based on Eq.(8.2). The tracking controller in Eq.(8.32) will however not result in an input that is constant over every time interval δ. To overcome this problem, we will replace u in Eq.(8.32) by its averaged value over every time interval δ. The error due to averaging ¯ is given by the following lemma: 172 Lemma 8.4.1. Consider the system in Eq.(8.34) with the input u replaced by ua : ¯ ¯ ¯x ¯¯ ˙ xa (t) = A¯a (t) + B ua , xa (t0 ) = x0 ¯ ¯ ¯¯ ya (t) = C xa (t) ¯ (8.38) where ua is the average of the Lipschitz control input u(t) over the time interval [t0 , t0 + δ]. ¯ ¯ Define exa (t) = x(t) − xa (t) and eya (t) = y (t) − ya (t), then: ¯ ¯ ¯ ¯ ¯ ¯ exa (t0 +δ) ≤ ¯ |¯ya (t0 +δ)| ≤ e ¯ ¯ eA(t0 +δ−τ ) B L δ 2 sup τ ∈[t0 ,t0 +δ] sup τ ∈[t0 ,t0 +δ] ¯ ¯ ¯ |CeA(t0 +δ−τ ) B| L δ 2 (8.39) where L is the Lipschitz constant for u(t) in time interval [t0 , t0 + δ]. ¯ Proof. For the linear systems in Eqs.(8.34) and (8.38), we can write: exa (t0 + δ) = ¯ t0 +δ t0 ¯ ¯u eA(t0 +δ−τ ) B[¯(τ ) − ua ] dτ ¯ (8.40) Form real analysis [147], for a bounded linear operator T (ˆ) such as the above integral, we x have: T (ˆ) ≤ T x x ˆ (8.41) where x ∈ A and ˆ ˆ T = sup{ T (ˆ) : x ∈ A, x ≤ 1} x ˆ ˆ ˆ 173 (8.42) Applying Eq.(8.41) to Eq.(8.40), we have: t0 +δ exa (t0 + δ) ≤ ≤ ≤ dτ sup t0 τ ∈[t0 ,t0 +δ] sup τ ∈[t0 ,t0 +δ] ¯ ¯u eA(t0 +δ−τ ) B[¯(τ ) − ua ] ¯ ¯ ¯ { eA(t0 +δ−τ ) B sup u(τ ) − ua } δ ¯ ¯ ¯ ¯ eA(t0 +δ−τ ) B L δ 2 (8.43) τ ∈[t0 ,t0 +δ] where we used u(τ ) − ua ≤ Lδ for the Lipschitz function u(t) in Eq.(8.32). The second ¯ ¯ ¯ ¯ ¯ inequality in Eq.(8.39) is similarly proved from eya = C exa . ¯ Let eo (t) = y(t) − r(t) and eo (t) = y (t) − r(t) be output tracking errors for the NMP and ¯ ¯ MP systems in Eqs.(8.36) and (8.34), respectively. The following theorem states a relation between these errors: Theorem 8.4.1. Consider the NMP system in Eq.(8.36) and MP system in Eq.(8.34) for t ∈ [t0 , t0 + δ] with the same initial conditions x(t0 ) = x(t0 ) = x0 . Assume that both systems ¯ are stable. For any desired reference trajectory r(t): (a) State variables x(t) are bounded for any bounded r(t) (b) eo (t0 + δ) = eo (t0 + δ) + eya (t0 + δ) ¯ ¯ if the system input is chosen as u(t) =    ua if ¯  0 t0 ≤ t ≤ t0 + pδ (8.44) if t0 + pδ < t < t0 + δ where 1 t0 +δ u(t) dt ¯ ua = ¯ δ t0 174 (8.45) and u(t) defined in Eq.(8.32) is the tracking control input designed for the MP system over ¯ the time interval [t0 , t0 + δ]. Proof. Eq.(8.38) is a DEMP system of the NMP system in Eq.(8.36) with the control input in Eq.(8.44). Therefore, their state variables will match at t = t0 + δ. Furthermore, the states of the NMP system will be bounded during the entire interval since linear systems do not have finite escape time. From Lemma 8.4.1, the difference between state variables of DEMP system in Eq.(8.38) and MP system in Eq.(8.34) at t = t0 + δ is a bounded value. Since the tracking control u(t) ensures the internal stability for the MP system in Eq.(8.34), ¯ we conclude that all state variables of the NMP system will be bounded for any bounded r(t). This completes the proof of part (a). The proof of part (b) is very similar and is not provided here. Corollary 8.4.1. For any time domain [ti , tf ], the results in Theorem 8.4.1 is valid for each time interval [ti + (j − 1)δ, ti + jδ] where j = 1, 2, · · · , if the initial time t0 and initial conditions x0 in Theorem 8.4.1 are redefined as below for each interval: t0 = ti + (j − 1)δ x0 = x(ti + (j − 1)δ) (8.46) This implies that the tracking control input u(t) for the MP system is designed with new ¯ initial conditions x(t0 ) = x0 for each time interval. ¯ Remark 8.4.1. The conditions for the IOT problem in Definition 8.4.1 can be satisfied for the NMP system in Eq.(8.36) using the results in Theorem 8.4.1 and Corollary 8.4.1 provided that: 175 (a) eya (t0 + δ) ≈ 0, which is achieved by choosing a sufficiently small switching interval δ. ¯ ¯ (b) eo (t0 + δ) → 0 as t → δ, which is achieved by choosing a sufficiently large gain K in ¯ Eq.(8.32). 8.4.3 Main Result We can summarize the proposed controller for intermittent output tracking in the following theorem: Theorem 8.4.2. Consider the linear SISO NMP system in Eq.(8.36). Assume that it is stable, without loss of generality. Let δ and p be the switching parameters, as shown in Figure 8.1, such that Eq.(8.34) is a DEMP system for Eq.(8.36). For sufficiently small δ and sufficiently large control gain, intermittent output tracking in the interval [ti , ∞], as defined in Definition 8.4.1, is achieved if the control input is chosen as: u(t) =    uaj if (j − 1)δ ≤ t − ti ≤ (j − 1 + p)δ ¯ (8.47) 1 ti +jδ u (t) dt ¯ δ ti +(j−1)δ j (8.48)  0 if for j = 1, 2, · · · where uaj = ¯ (j − 1 + p)δ < t − ti < jδ and uj (t) in Eq.(8.32) is the tracking control input designed for the MP system for each time ¯ interval [ti + (j − 1)δ, ti + jδ] with initial conditions given in Eq.(8.46). Proof. The proof is followed from earlier discussions in sections 8.4.1 and 8.4.2. Remark 8.4.2. If the NMP system is not stable, a state-feedback controller can be used first to stabilize the system provided (A, B) are stabilizable. Then, the discussion to design the 176 IOT controller holds for the stabilized system. In this case, the control input in Eq.(8.44) is changed to: u(t) = −Kx(t) +    ua if ¯  0 t0 ≤ t ≤ t0 + pδ (8.49) if t0 + pδ < t < t0 + δ where K is the state-feedback gain chosen such that A − BK is Hurwitz. It should be noted that the switching parameters δ and p are found for the stabilized system (A − BK, B, C). For the case that the NMP system is stable and (A, B) is controllable, the state-feedback control in Eq.(8.49) can be used to move the poles of the system. A relocation of the poles can be helpful in finding appropriate switching parameters for the NMP system. Remark 8.4.3. The first necessary step in solving the IOT problem for a NMP system is to find a DEMP system. This can be done following the discussion in Section 8.3 through proper choice of parameters δ and p. Although there is no guarantee that a DEMP system exists for the general case, a numerical search for a range of δ values can be carried out. If there are multiple choices for the pair (δ, p) for DEMP systems, a trade off between IOT error, eo (t0 + δ), and the input switching frequency can be used to choose the appropriate DEMP system. A smaller value of δ implies less IOT error but higher switching frequency. 8.5 Numerical Simulations In this section, we present two numerical simulations for solving the IOT problem for NMP systems. The system parameters and reference trajectories are chosen in a way that standard continuous-time output regulation methods are inapplicable. We summarize the steps to design the IOT controller as follows: 1. If the NMP system is unstable and stabilizable, design a state-feedback controller to 177 stabilize the system. 2. Find switching parameters δ and p such that all RHP zeros of the stable NMP system are moved to open LHP. This can be done following the discussion in Section 8.3. 3. For the resulting MP system in step 2, design the tracking control in Eq.(8.32) for each time interval [ti + (j − 1)δ, ti + jδ] where j = 1, 2, · · · , with initial conditions given in Eq.(8.46). 4. Apply the control input in Eq.(8.44) for each time interval. In the case of unstable NMP system, the control input is given by Eq.(8.49). For the first simulation, consider the following stable NMP system: G(s) = 10 s(s − 1) (s + 1)2 (s + 10) (8.50) and let the reference trajectory be defined by: r(t) =    exp(t) − 1 ,   exp(3) − 1 , t≤3 (8.51) t>3 This reference trajectory can be generated by the following exosystem:   1 0 w = Sw , S =  ˙  0 0 (8.52) Since the unstable eigenvalue of the exosystem is identical to the RHP zero of the system, this tracking problem can not be solved by output regulation techniques [108]. The choice 178 y, r(t) 20 a 10 η1 , η2 (t) 0 0.4 b 0.2 0 -0.2 u, u(t) ¯ 0 c -20 -40 0 0.75 1.5 3 2.25 3.75 4.5 5.25 6 t(sec) Figure 8.3: Simulation results for the first case showing the system output (solid line) and reference trajectory (dashed line) in the first plot, internal states in the second plot and control input (Solid line) and tracking control for MP system (dashed line) in the third plot of switching parameters δ = 0.75 and p = 0.75 yields the following DEMP system: (s + 110.22)(s + 0.28) ¯ G(s) = −0.17 (s + 1)2 (s + 10) (8.53) Figure 8.3 shows the simulation results for the system in Eq.(8.50). The system output (solid line) and the reference trajectory (dashed line) are presented in the first plot in Figure 8.3. The second plot shows the system internal states, η1 (t) and η2 (t), which represent the internal dynamics of the system in the normal form. The third plot shows the control input for NMP system (solid line) and the tracking control for MP system which are related by Eq.(8.45) for each time interval δ. These time intervals are illustrated by vertical dashed lines in Figure 8.3. Simulation results ensure that the system output meets the reference trajectory with a sufficiently small error while the internal states of the NMP are bounded. 179 For the second simulation, consider the following unstable NMP system: G(s) = (s2 + 1)(s − 15) (s − 1)(s − 2)(s2 + 2s + 5) (8.54) with complex zeros and poles. The reference trajectory is defined by: r(t) = 5 arctan(t) (8.55) This trajectory can not be generated by an exosystem so the output regulation methods can not be applied to this example. Since the NMP system is controllable, a state-feedback control is first designed to stabilize the system by moving the unstable poles s = 1 and s = 2 to s = −1 and s = −2, respectively. Through a numerical search, the switching parameters are chosen as δ = 0.5 and p = 0.7 for the stabilized system and result in the following DEMP system: 2 (s + 25.02)(s + 0.08s + 1) ¯ G(s) = −0.45 (s + 1)(s + 2)(s2 + 2s + 5) (8.56) The simulation results are shown in Figure 8.4. It includes the system output (solid line) and reference trajectory (dashed line) in the first plot, the internal states in the normal form of the system in the second plot and the control input in the third plot. The control input is comprised of a continuous state-feedback control input and a switched control input using (8.49). As can be seen from the figure, the system output is approximately matched with the reference trajectory at regular intervals of time δ shown by vertical dashed lines. Furthermore, the internal states remain bounded for the NMP system. Remark 8.5.1. The deviation of the system output from the desired trajectory during the 180 y, r(t) 8 a 4 0 0.5 η1 , η2 , η3 (t) 0 u(t) b 0 -0.5 10 c -10 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t(sec) Figure 8.4: Simulation results for the second case showing the system output (solid line) and reference trajectory (dashed line) in the first plot, internal states in the second plot and control input (Solid line) and tracking control for MP system (dashed line) in the third plot time interval t ∈ (t0 , t0 + δ) depends on the system parameters as well as the switching parameters δ and p. In the case that there exist multiple choices for switching parameters, another criterion can be developed to minimize this deviation. This can be investigated more as a future work. 8.6 Conclusion A continuous-time control method for intermittent output tracking in linear SISO NMP systems was proposed. The concept of discrete equivalent system have been used to move the unstable zeros to the LHP. The switched control input is obtained from a tracking control input for the discrete equivalent minimum-phase system. The proposed method ensures that the system output matches with the reference trajectory at regular intervals of time while 181 keeping the internal dynamics stable. Although the controller does not guarantee exact output matching with the reference trajectory, the matching error can be made sufficiently small by adjusting the control gains and switching parameters of the input. The reference trajectories to be tracked are not constrained to any specific set, however, a finite preview of the trajectory is required. The switched inputs used in the control can be estimated by impulsive inputs as their magnitudes get larger and their time supports get smaller. 182 Chapter 9 Conclusion and Future Work In this research, we investigated different aspects of impulsive control in underactuated mechanical systems. As the first step, we ensured that impulsive inputs are indeed implementable in mechanical systems using standard hardware. This was proven through an experimental verification of a swing-up algorithm for the Pendubot. Having addressed the implementability concern, we used the properties of impulsive inputs to develop several algorithms for underactuated mechanical systems with different control objectives. Impulsive inputs are capable of causing instantaneous jumps in the energy of a mechanical system. This can be used to develop new energy-based control algorithms or improve current algorithms in the literature. To illustrate this, we considered the swing-up control of the Acrobot, which is a benchmark problem in underactuated mechanical systems and compared our results with those in the literature. To show another aspect of impulsive inputs, we looked at the stabilization problem for underactuated systems. Impulsive inputs can also be useful in this category of controllers as they can result in instantaneous jumps in the Lyapunov function corresponding to a stabilizing controller. We used this property to develop a stabilization algorithm which results in enlarging the region of attraction of equilibria for underactuated systems. Several examples with linear and nonlinear controllers were studied and it was shown that adding impulsive inputs at some specific instants of time to the continuous control inputs can enlarge the region of attraction for a desired equilibrium. As another property of impulsive inputs, they can affect the dynamics of a system significantly and in a 183 very short period of time. This is beneficial especially when there is not much time available for continuous inputs to control the system behavior. We show this advantage with a control problem of changing the fall direction for humanoid robots in order to avoid damaging the surrounding objects. A six-dof robot was considered and two sets of impulsive inputs were applied to change the direction of the fall to a desired angle. Impulsive inputs were also shown to be useful in rejecting external disturbances with the help of balance maintenance control for the Synthetic-Wheel biped. The last two pieces of work are not using impulsive inputs but still related to the concept of impulsive control for underactuated systems. The first one dealt with designing zero-energy walking gaits for underactuated bipeds by avoiding undesirable impacts from the ground at the leg interchange. The second work proposed an intermittent output tracking control for non-minimum-phase systems. The switched inputs can be approximated by impulsive inputs in the limit. We suggest the following to continue this field of research further in the future. Although the experiments on the Pendubot showed satisfactory results, a new amplifier circuit may be designed to get closer to the ideal impulses in experiments. The new circuit will have a large capacitor which enables the amplifier to provide enough current for the motor at the impulsive times. The experimental verification of all the algorithms developed in this dissertation would add a lot to prove that impulsive control is a reliable method for underactuated mechanical systems. For the work presented in Chapter 4, more analysis is needed to predict if the impulsive algorithm is capable of stabilizing the desired equilibrium from the given configuration. The safe fall algorithm in Chapter 5 is a preliminary work to prove the concept of using impulsive inputs for such platforms. A more realistic configuration of a humanoid with higher dof’s and two legs should be studied. There is some flexibility in the algorithm which may be used to optimize the safe fall algorithm for a specific humanoid. For 184 example, the impulsive times which are based on specific configurations of the robot may be changed to get better results. Energy-conserving gaits which are presented in Chapter 7 are not practical yet and the limitations of the generating such gaits in actual bipeds need to be considered. For example, generated gaits for five-dof bipeds require the knee bent backward or the swing-leg touches the floor in the middle of its swing phase which are not practical assumptions. In addition, finding proper electrical circuits which are capable of regenerating and storing energy from negative work during each step is essential. For the algorithm in Chapter 8, future work may focus on extending the current results to multi-input multi-output linear systems. Also, the effect of switching parameters on output deviations from desired trajectory can be investigated to find the best control performance in between intermittent tracking times. 185 APPENDICES 186 Appendix A Procedure for Verifying the Condition in Theorem 3.4.1 Assume that the gains kp , kd , and ke , and the initial condition x0 are chosen based on the conditions stated in [3], [22]. This guarantees asymptotic convergence of the acrobot to its upright configuration using the continuous controller designed in those works. We now outline the procedure to check the condition in Eq.(3.64) in Theorem 3.4.1: ¯ ¯ 1. Find all points x = (θ1 , θ2 , 0, 0) in the set Ω defined in Eq.(3.44) using the relations in ¯ ¯ [3]. To this end, we first solve the following equation to find θ2 : kp ¯ θ2 = g 2 q4 q5 g ¯ sin(θ2 ) 2 2 ¯ q4 +q5 +2q4 q5 cos(θ2 ) [g 2 2 ¯ q4 +q5 +2q4 q5 cos(θ2 )+0.5Edes ] (A.1) ¯ Then find θ1 from the following equation: kp θ2 + (E − Edes )τh = 0 (A.2) ¯ ¯ ˙− ˙ − 2. Find the points (θ1 , θ2 , θ1 , θ2 ) in the set S2 (x) defined in Eq.(3.63). Using Eq.(3.18) ˙+ ˙+ with θ1 = θ2 = 0, we have: ˙− θ1 = − q2 + q3 C2 ˙ θ− q1 + q2 + 2q3 C2 2 187 (A.3) and from the definition of S2 (x) we know: ¨ ˙ θ2 = a(θ, θ− ) + b(θ2 )τc = 0 (A.4) Using the expression for τ2c in Eqs.(3.38) or (3.69), as appropriate, Eqs.(A.3) and (A.4) ˙− ˙− can be solved to find the values of θ1 and θ2 . 3. For the points in S2 (x) which satisfy V (x) ≤ V (x0 ), check the condition in Eq.(3.64) ˙ (3) (3) ¨ by computing θ2 θ2 . θ2 is computed by taking the derivative of θ2 given in Eq.(3.8). 188 Appendix B Equations of Motion: Synthetic-Wheel Biped The symmetric matrix M, matrix N, and vector G in Eq.(7.1) have the following entries: 2 M(1, 1) = It + Ist + Isw +mt d2 +lt +R2 +2R(dt −lt ) cos(θ + φ)−2dt lt t + mst d2 + 2R(R + (dst − R) cos(θ) − dst ) st +msw d2 +2R(R+(dsw −R) cos(θ+ψ)−dsw ) sw M(1, 2) = It + mt (dt − lt )(dt − lt + R cos(θ + φ)) M(1, 3) = Isw + msw (dsw − R)(dsw − R + R cos(θ + ψ)) M(2, 2) = It + mt (dt − lt )2 M(2, 3) = 0 M(3, 3) = Isw + msw (dsw − R)2 ˙ ˙ ˙ N(1, 1) = mt R(lt − dt ) sin(θ + φ)(θ + φ) + mst R(R − dst ) sin(θ)θ ˙ ˙ + msw R(R − dsw ) sin(θ + ψ)(θ + ψ) ˙ ˙ N(1, 2) = mt R(lt − dt ) sin(θ + φ)(θ + φ) ˙ ˙ N(1, 3) = msw R(R − dsw ) sin(θ + ψ)(θ + ψ) 189 N(2, 1) = N(2, 2) = N(2, 3) = 0 N(3, 1) = N(3, 2) = N(3, 3) = 0 G(1) = [mt (lt − dt ) sin(θ + φ) + mst (R − dst ) sin(θ) + msw (R − dsw ) sin(θ + ψ)] g G(2) = mt (lt − dt ) sin(θ + φ) g G(3) = msw (R − dsw ) sin(θ + ψ) g where mt , mst and msw are the masses of the torso, stance leg and swing leg; It , Ist and Isw are the mass moments of inertia of the torso, stance leg and swing leg about their respective center of mass; R is the radius of curvature of the feet; lt is the length of the torso; dt , dst and dsw are the distances shown in Figure 7.2, and g is the acceleration due to gravity. 190 Appendix C Equations of Motion: Five-DOF Biped The symmetric matrix M, matrix N, and vector G in Eq.(7.1) have the following entries: 2 M(1, 1) = I1 + d2 m1 + l1 (m2 + m3 + m4 + m5 ) 1 M(1, 2) = l1 [d2 m2 + l2 (m3 + m4 + m5 )] cos(θ1 − θ2 ) M(1, 3) = d3 l1 m3 cos(θ1 − θ3 ) M(1, 4) = l1 (d4 m4 + l4 m5 ) cos(θ1 − θ4 ) M(1, 5) = d5 l1 m5 cos (θ1 − θ5 ) 2 M(2, 2) = I2 + d2 m2 + l2 (m3 + m4 + m5 ) 2 M(2, 3) = d3 l2 m3 cos(θ2 − θ3 ) M(2, 4) = l2 (d4 m4 + l4 m5 ) cos(θ2 − θ4 ) M(2, 5) = d5 l2 m5 cos(θ2 − θ5 ) M(3, 3) = I3 d2 m3 3 M(3, 4) = M(3, 5) = 0 2 M(4, 4) = I4 + d2 m4 l4 m5 4 M(4, 5) = d5 l4 m5 cos(θ4 − θ5 ) 191 M(5, 5) = I5 + d2 m5 5 N(1, 1) = N(2, 2) = N(3.3) = N(4, 4) = N(5, 5) = 0 ˙ N(1, 2) = l1 [d2 m2 + l2 (m3 + m4 + m5 )] sin(θ1 − θ2 )θ2 ˙ N(1, 3) = d3 l1 m3 sin(θ1 − θ3 )θ3 ˙ N(1, 4) = l1 (d4 m4 + l4 m5 ) sin(θ1 − θ4 )θ4 ˙ N(1, 5) = d5 l1 m5 sin(θ1 − θ5 )θ5 ˙ N(2, 1) = − l1 [d2 m2 +l2 (m3 +m4 +m5 )] sin(θ1 −θ2 )θ1 ˙ N(2, 3) = d3 l2 m3 sin(θ2 − θ3 )θ3 ˙ N(2, 4) = l2 (d4 m4 + l4 m5 ) sin(θ2 − θ4 )θ4 ˙ N(2, 5) = d5 l1 m5 sin(θ1 − θ5 )θ5 ˙ N(3, 1) = − d3 l1 m3 sin(θ1 − θ3 )θ1 ˙ N(3, 2) = − d3 l2 m3 sin(θ2 − θ3 )θ2 N(3, 4) = N(3, 5) = N(4, 3) = N(5, 3) = 0 ˙ N(4, 1) = − l1 (d4 m4 + l4 m5 ) sin(θ1 − θ4 )θ1 ˙ N(4, 2) = − l2 (d4 m4 + l4 m5 ) sin(θ2 − θ4 )θ2 ˙ N(4, 5) = d5 l4 m5 sin(θ4 − θ5 )θ5 ˙ N(5, 1) = − d5 l1 m5 sin(θ1 − θ5 )θ1 ˙ N(5, 2) = − d5 l2 m5 sin(θ2 − θ5 )θ2 ˙ N(5, 4) = − d5 l4 m5 sin (θ4 − θ5 )θ4 G(1) = − g[d1m1 + l1 (m2 + m3 + m4 + m5 )] sin θ1 G(2) = − g[d2m2 + l2 (m3 + m4 + m5 )] sin θ2 192 G(3) = − gd3m3 sin θ3 G(4) = − g(d4m4 + l4 m5 ) sin θ4 G(5) = − gd5m5 sin θ5 where mi , Ii , and li , i = 1, · · · , 5, are the mass, mass moment of inertia, and length of the i-th link; di is the distance of the center of the mass of the i-th link from the i-th joint; and g is the acceleration due to gravity. 193 BIBLIOGRAPHY 194 BIBLIOGRAPHY [1] L. Flynn, R. Jafari, and R. Mukherjee, “Active synthetic-wheel biped with torso,” Robotics, IEEE Transactions on, vol. 26, no. 5, pp. 816–826, 2010. [2] T. Albahkali, R. Mukherjee, and T. Das, “Swing-up control of the pendubot: an impulse–momentum approach,” Robotics, IEEE Transactions on, vol. 25, no. 4, pp. 975–982, 2009. [3] X. Xin and M. Kaneda, “Analysis of the energy-based swing-up control of the Acrobot,” International Journal of Robust and Nonlinear Control, vol. 17, no. 16, pp. 1503–1524, 2007. [4] T. Pavlidis, “Stability of systems described by differential equations containing impulses,” IEEE Transactions on Automatic Control, vol. 12, no. 1, pp. 43–45, 1967. [5] E. Gilbert and G. Harasty, “A class of fixed-time fuel-optimal impulsive control problems and an efficient algorithm for their solution,” Automatic Control, IEEE Transactions on, vol. 16, no. 1, pp. 1–11, 1971. [6] J. Menaldi, “The separation principle for impulse control problems,” Proceedings of the American Mathematical Society, vol. 82, no. 3, pp. 439–445, 1981. [7] V. Lakshmikantham, D. Bainov, and P. Simeonov, Theory of impulsive differential equations. World Scientific Pub Co Inc, 1989. [8] J. Sun, Y. Zhang, and Q. Wu, “Less conservative conditions for asymptotic stability of impulsive control systems,” IEEE Transactions on automatic control, vol. 48, no. 5, p. 829, 2003. [9] I. Stamova, Stability Analysis of Impulsive Functional Differential Equations. Walter de Gruyter, 2009. [10] G. Xie and L. Wang, “Necessary and sufficient conditions for controllability and observability of switched impulsive control systems,” IEEE Transactions on Automatic Control, vol. 49, no. 6, pp. 960–966, 2004. 195 [11] T. Yang, “Impulsive control,” Automatic Control, IEEE Transactions on, vol. 44, no. 5, pp. 1081–1083, 1999. [12] Y. Orlov, “Nonlinear control systems with impulsive inputs,” in Decision and Control, 1997., Proceedings of the 36th IEEE Conference on, vol. 1, pp. 630–635, IEEE, 1997. [13] Z. Guan, G. Chen, and T. Ueta, “On impulsive control of a periodically forced chaotic pendulum system,” Automatic Control, IEEE Transactions on, vol. 45, no. 9, pp. 1724– 1727, 2000. [14] S. Weibel, G. Howell, and J. Baillieul, “Control of single-degree-of-freedom hamiltonian systems with impulsive inputs,” in Decision and Control, 1996., Proceedings of the 35th IEEE, vol. 4, pp. 4661–4666, IEEE, 1996. [15] Y. Aoustin, D. Romero, C. Chevallereau, and S. Aubin, “Impulsive control for a thirteen-link biped,” in 9th IEEE International Workshop on Advanced Motion Control, pp. 439–444, 2006. [16] Y. Orlov, R. Santiesteban, and L. Aguilar, “Impulsive control of a mechanical oscillator with friction,” in American Control Conference, 2009. ACC’09., pp. 3494–3499, IEEE, 2009. [17] N. van de Wouw and R. Leine, “Impulsive control of mechanical motion systems with uncertain friction,” in Decision and Control and European Control Conference (CDCECC), 2011 50th IEEE Conference on, pp. 4176–4182, IEEE, 2011. [18] G. Boone, “Minimum-time control of the acrobot,” in IEEE International Conference on Robotics and Automation, pp. 3281–3287, Citeseer, 1997. [19] M. Spong, “The swing up control problem for the acrobot,” IEEE Control Systems Magazine, vol. 15, no. 1, pp. 49–55, 1995. [20] I. Fantoni, R. Lozano, and M. Spong, “Energy based control of the pendubot,” Automatic Control, IEEE Transactions on, vol. 45, no. 4, pp. 725–729, 2000. [21] O. Kolesnichenko and A. Shiriaev, “Partial stabilization of underactuated euler– lagrange systems via a class of feedback transformations,” Systems & control letters, vol. 45, no. 2, pp. 121–132, 2002. [22] A. Mahindrakar and R. Banavar, “A swing-up of the acrobot based on a simple pendulum strategy,” International Journal of Control, vol. 78, no. 6, pp. 424–429, 2005. 196 [23] R. Banavar and A. Mahindrakar, “A non-smooth control law and time-optimality notions for the acrobot,” International Journal of Control, vol. 78, no. 15, pp. 1166– 1173, 2005. [24] A. Bloch, N. Leonard, and J. Marsden, “Controlled lagrangians and the stabilization of mechanical systems. i. the first matching theorem,” Automatic Control, IEEE Transactions on, vol. 45, no. 12, pp. 2253–2270, 2000. [25] A. Bloch, D. Chang, N. Leonard, and J. Marsden, “Controlled lagrangians and the stabilization of mechanical systems. ii. potential shaping,” Automatic Control, IEEE Transactions on, vol. 46, no. 10, pp. 1556–1571, 2001. [26] R. Ortega, A. Loria, R. Kelly, and L. Praly, “On passivity-based output feedback global stabilization of euler-lagrange systems,” International Journal of Robust and Nonlinear Control, vol. 5, no. 4, pp. 313–323, 1995. [27] R. Ortega, A. Van Der Schaft, B. Maschke, and G. Escobar, “Interconnection and damping assignment passivity-based control of port-controlled hamiltonian systems,” Automatica, vol. 38, no. 4, pp. 585–596, 2002. [28] R. Ortega, M. Spong, F. G´mez-Estern, and G. Blankenstein, “Stabilization of a class o of underactuated mechanical systems via interconnection and damping assignment,” Automatic Control, IEEE Transactions on, vol. 47, no. 8, pp. 1218–1233, 2002. [29] J. Acosta, R. Ortega, A. Astolfi, and A. Mahindrakar, “Interconnection and damping assignment passivity-based control of mechanical systems with underactuation degree one,” Automatic Control, IEEE Transactions on, vol. 50, no. 12, pp. 1936–1955, 2005. [30] D. Auckly, L. Kapitanski, and W. White, “Control of nonlinear underactuated systems,” Communications on Pure Applied Mathematics, vol. 53, pp. 354–369, 2000. [31] D. Auckly and L. Kapitanski, “On the λ-equations for matching control laws,” Siam Journal on Control and Optimization, vol. 41, no. 5, pp. 1372–1388, 2003. [32] F. Andreev, D. Auckly, S. Gosavi, L. Kapitanski, A. Kelkar, and W. White, “Matching, linear systems, and the ball and beam,” Automatica, vol. 38, no. 12, pp. 2147–2152, 2002. [33] W. White, M. Foss, and X. Guo, “A direct lyapunov approach for a class of underactuated mechanical systems,” in American Control Conference, 2006, pp. 103–110, IEEE, 2006. 197 [34] W. White, M. Foss, and X. Guo, “A direct lyapunov approach for stabilization of underactuated mechanical systems,” in American Control Conference, 2007. ACC’07, pp. 4817–4822, IEEE, 2007. [35] W. White, M. Foss, J. Patenaude, X. Guo, and D. Garc´ “Improvements in direct ıa, lyapunov stabilization of underactuated, mechanical systems,” in American Control Conference, 2008, pp. 2927–2932, IEEE, 2008. [36] R. Olfati-Saber, Nonlinear control of underactuated mechanical systems with application to robotics and aerospace vehicles. PhD thesis, Massachusetts Institute of Technology, 2000. [37] R. Olfati-Saber, “Nonlinear control and reduction of underactuated systems with symmetry. i. actuated shape variables case,” in Decision and Control, 2001. Proceedings of the 40th IEEE Conference on, vol. 5, pp. 4158–4163, IEEE, 2001. [38] R. Olfati-Saber, “Nonlinear control and reduction of underactuated systems with symmetry. ii. unactuated shape variables case,” in Decision and Control, 2001. Proceedings of the 40th IEEE Conference on, vol. 5, pp. 4164–4169, IEEE, 2001. [39] R. Olfati-Saber, “Nonlinear control and reduction of underactuated systems with symmetry. iii. input coupling case,” in Decision and Control, 2001. Proceedings of the 40th IEEE Conference on, vol. 4, pp. 3778–3783, IEEE, 2001. ¨ ¨ u [40] R. Xu and U. Ozg¨ ner, “Sliding mode control of a class of underactuated systems,” Automatica, vol. 44, no. 1, pp. 233–241, 2008. [41] S. Riachy, Y. Orlov, T. Floquet, R. Santiesteban, and J. Richard, “Second-order sliding mode control of underactuated mechanical systems i: Local stabilization with application to an inverted pendulum,” International Journal of Robust and Nonlinear Control, vol. 18, no. 4-5, pp. 529–543, 2008. [42] R. Martinez, J. Alvarez, and Y. Orlov, “Hybrid sliding-mode-based control of underactuated systems with dry friction,” Industrial Electronics, IEEE Transactions on, vol. 55, no. 11, pp. 3998–4003, 2008. [43] T. Wilken, M. Missura, and S. Behnke, “Designing falling motions for a humanoid soccer goalie,” in Proceedings of the 4th Workshop on Humanoid Soccer Robots (Humanoids 2009), pp. 79–84, 2009. 198 [44] J.-S. Tan, J. J. Eng, S. R. Robinovitch, and B. Warnick, “Wrist impact velocities are smaller in forward falls than backward falls from standing,” vol. 39, no. 10, pp. 1804– 1811, 2006. [45] S.-K. Yun, A. Goswami, and Y. Sakagami, “Safe fall: Humanoid robot fall direction change through intelligent stepping and inertia shaping,” pp. 781–787, May 2009. [46] U. Nagarajan and A. Goswami, “Generalized direction changing fall control of humanoid robots among multiple objects,” (Anchorage, Alaska), pp. 3316–3322, May 2010. [47] S. Kalyanakrishnan and A. Goswami, “Learning to predict humanoid fall,” vol. 8, no. 2, pp. 245–273, 2011. [48] K. Fujiwara, K. F., H. Saito, S. Kajita, K. Harada, and H. Hirukawa, “Falling motion control of a humanoid robot trained by virtual supplementary tests,” pp. 1077–1082, 2004, New Orleans, LA, USA. [49] K. Fujiwara, S. Kajita, K. Harada, K. Kaneko, M. Morisawa, F. Kanehiro, S. Nakaoka, S. Harada, and H. Hirukawa, “Towards an optimal falling motion for a humanoid robot,” in IEEE-RAS/RSJ International Conference on Humanoid Robots, pp. 524– 529, 2006. [50] K. Fujiwara, F. Kanehiro, S. Kajita, K. Kaneko, K. Yokoi, and H. Hirukawa, “UKEMI: Falling motion control to minimize damage to biped humanoid robot,” pp. 2521–2526, September 30 – October 4, 2002 Lausanne, Switzerland. [51] K. Fujiwara, F. Kanehiro, S. Kajita, and H. Hirukawa, “Safe knee landing of a humansize humanoid robot while falling forward,” pp. 503–508, September 28– October 2 2004, Sendai, Japan. [52] K. Ogata, K. Terada, and Y. Kuniyoshi, “Real-time selection and generation of fall damagae reduction actions for humanoid robots,” in IEEE-RAS/RSJ International Conference on Humanoid Robots, pp. 233–238, Dec. -3 2008, Daejeon, Korea. [53] K. Ogata, K. Terada, and Y. Kuniyoshi, “Falling motion control for humanoid robots while walking,” in IEEE-RAS/RSJ International Conference on Humanoid Robots, Pittsburgh, 2007. [54] R. Renner and S. Behnke, “Instability detection and fall avoidance for a humanoid using attitude sensors and reflexes,” (Beijing), pp. 2967–2973, October 2006. 199 [55] J. G. D. Karssen and M. Wisse, “Fall detection in walking robots by multi-way principal component analysis,” Robotica, vol. 27, no. 2, pp. 249–257, 2009. [56] J. R. del Solar, R. Palma-Amestoy, R. Marchant, I. Parra-Tsunekawa, and P. Zegers, “Learning to fall: Designing low damage fall sequences for humanoid soccer robots,” vol. 57, no. 8, pp. 796–807, 2009. [57] O. H¨hn, J. Gaˇnik, and W. Gerth, Climbing and Walking Robots, ch. Detection and o c Classification of Posture Instabilities of Bipedal Robots, pp. 409–416. Springer Berlin Heidelberg, 2006. [58] O. H¨hn and W. Gerth, “Probabilistic balance monitoring for bipedal robots,” vol. 28, o pp. 245–256, 2009. [59] T. Ishida, Y. Kuroki, and T. Takahashi, “Analysis of motions of a small biped entertainment robot,” pp. 142–147, September 28 – October 2, 2004 Sendai, Japan. [60] R. Kanoi and C. Hartland, “Fall detections in humanoid walk patterns using reservoir computing based control architecture,” in 5th national conference on Control Architecture of Robots, (Douai, France), may 2010. [61] A. Forner Cordero, Human gait, stumble and . . . fall? Twente, Enschede, The Netherlands, 2003. PhD thesis, University of [62] F. Horak and L. Nashner, “Central programming of postural movements: adaptation to altered support-surface configurations,” Journal of Neurophysiology, vol. 55, no. 6, pp. 1369–1381, 1986. [63] B. Maki and W. McIlroy, “The role of limb movements in maintaining upright stance: the change-in-support strategy,” Physical Therapy, vol. 77, no. 5, pp. 488–507, 1997. [64] A. Hofmann, Robust execution of bipedal walking tasks from biomechanical principles. PhD thesis, Massachusetts Institute of Technology, 2005. [65] K. Hirai, M. Hirose, Y. Haikawa, and T. Takenaka, “The development of Honda humanoid robot,” in IEEE International Conference on Robotics and Automation, pp. 1321–1326, 1998. [66] K. Mitobe, G. Capi, and Y. Nasu, “Control of walking robots based on manipulation of the zero moment point,” Robotica, vol. 18, no. 06, pp. 651–657, 2001. 200 [67] J. Park and H. Chung, “ZMP compensation by online trajectory generation for biped robots,” in 1999 IEEE International Conference on Systems, Man, and Cybernetics, 1999. IEEE SMC’99 Conference Proceedings, vol. 4, 1999. [68] M. Abdallah and A. Goswami, “A biomechanically motivated two-phase strategy for biped upright balance control,” in Robotics and Automation, 2005. ICRA 2005. Proceedings of the 2005 IEEE International Conference on, pp. 1996–2001, 2005. [69] V. Prahlad, G. Dip, and C. Meng-Hwee, “Disturbance rejection by online ZMP compensation,” Robotica, vol. 26, no. 01, pp. 9–17, 2007. [70] A. Goswami and V. Kallem, “Rate of change of angular momentum and balance maintenance of biped robots,” in Robotics and Automation, 2004. Proceedings. ICRA’04. 2004 IEEE International Conference on, vol. 4, pp. 3785–3790, IEEE, 2004. [71] B. Stephens, “Humanoid push recovery,” in Humanoid Robots, 2007 7th IEEE-RAS International Conference on, pp. 589–595, Ieee, 2007. [72] S. Yi, B. Zhang, D. Hong, and D. Lee, “Learning full body push recovery control for small humanoid robots,” in Robotics and Automation (ICRA), 2011 IEEE International Conference on, pp. 2047–2052, IEEE, 2011. [73] K. Shuuji and T. Kazuo, “Study of dynamic biped locomotion on rugged terrain: Derivation and application of the linear inverted pendulum mode,” in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 1405–1411, 1991. [74] S. Kajita, F. Kanehiro, K. Kaneko, K. Yokoi, and H. Hirukawa, “The 3d linear inverted pendulum mode: A simple modeling for a biped walking pattern generation,” in Intelligent Robots and Systems, 2001. Proceedings. 2001 IEEE/RSJ International Conference on, vol. 1, pp. 239–246, Ieee, 2001. [75] H. Hemami and P. Camana, “Nonlinear feedback in simple locomotion systems,” Automatic Control, IEEE Transactions on, vol. 21, no. 6, pp. 855–860, 1976. [76] J. Pratt, J. Carff, S. Drakunov, and A. Goswami, “Capture point: A step toward humanoid push recovery,” in Humanoid Robots, 2006 6th IEEE-RAS International Conference on, pp. 200–207, IEEE, 2006. [77] S. Kudoh and T. Komura, “C2 continuous gait-pattern generation for biped robots,” in Intelligent Robots and Systems, 2003.(IROS 2003). Proceedings. 2003 IEEE/RSJ International Conference on, vol. 2, pp. 1135–1140, IEEE, 2003. 201 [78] J. Rebula, F. Canas, J. Pratt, and A. Goswami, “Learning capture points for humanoid push recovery,” in Humanoid Robots, 2007 7th IEEE-RAS International Conference on, pp. 65–72, IEEE, 2007. [79] B. Stephens and C. Atkeson, “Dynamic balance force control for compliant humanoid robots,” in Intelligent Robots and Systems (IROS), 2010 IEEE/RSJ International Conference on, pp. 1248–1255, IEEE, 2010. [80] B. Stephens and C. Atkeson, “Push recovery by stepping for humanoid robots with force controlled joints,” in Humanoid Robots (Humanoids), 2010 10th IEEE-RAS International Conference on, pp. 52–59, IEEE, 2010. [81] T. Komura, H. Leung, S. Kudoh, and J. Kuffner, “A feedback controller for biped humanoids that can counteract large perturbations during gait,” in Robotics and Automation, 2005. ICRA 2005. Proceedings of the 2005 IEEE International Conference on, pp. 1989–1995, IEEE, 2005. [82] A. Adiwahono, C. Chew, W. Huang, and Y. Zheng, “Push recovery controller for bipedal robot walking,” in Advanced Intelligent Mechatronics, 2009. AIM 2009. IEEE/ASME International Conference on, pp. 162–167, IEEE, 2009. [83] T. McGeer, “Passive dynamic walking,” vol. 9, no. 2, pp. 62–82, 1990. [84] D. G. E. Hobbelen and M. Wisse, “Controlling the walking speed in limit cycle walking,” vol. 27, no. 9, pp. 989–1005, 2008. [85] M. Srinivasan and A. Ruina, “Computer optimization of a minimal biped model discovers walking and running,” Nature, vol. 439, no. 5, pp. 72–75, 2006. [86] S. H. Collins and A. Ruina, “A bipedal walking robot with efficient and human-like gait,” pp. 1983–1988, 2005. [87] K. Sreenath, H. W. Park, I. Poulakakis, and J. W. Grizzle, “A compliant hybrid zero dynamics controller for stable, efficient and fast bipedal walking on mabel,” vol. 30, no. 9, pp. 1170–1193, 2010. [88] Q. Huang, T. Hase, and K. Ono., “Optimal trajectory planning method using inequality state constraint for a biped walking robot with upper body mass,” Journal of System Design and Dynamics, vol. 1, pp. 168–179, 2007. 202 [89] E. Dertien, “Dynamic walking with dribbel - design and construction of a passivitybased walking robot,” IEEE Robotics & Automation Magazine, vol. 13, no. 3, pp. 118– 122, 2006. [90] P. Bhounsule, J. Cortell, B. Hendriksen, J. G. D. Karssen, C. Paul, and A. Ruina, “A robot that walked 65 km on a single charge: Energy-effective and reliable locomotion using trajectory optimization and stabilization from reflexes,” submitted. [91] M. J. Donelan, R. Kram, and A. D. Kuo, “Mechanical work for step-to-step transitions is a major determinant of the metabolic cost of human walking,” The Journal of Experimental Biology, vol. 205, pp. 3717–3727, 2002. [92] Y. Sakagami, R. Watanabe, C. Aoyama, S. Matsunaga, N. Higaki, and K. Fujimura, “The intelligent Asimo: System overview and integration,” vol. 3, pp. 2478–2483, 2002. [93] M. Hirose and K. Ogawa, “Honda humanoid robotics development,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 365, pp. 11–19, 2007. [94] H. Hirukawa, F. Kanehiro, K. Kaneko, S. Kajita, K. Fujiwara, Y. Kawai, F. Tomita, S. Hirai, K. Tanie, T. Isozumi, K. Akachi, T. Kawasaki, S. Ota, K. Yokoyama, H. Handa, Y. Fukase, J. Maeda, Y. Nakamura, S. Tachi, and H. Inoue, “Humanoid robotics platforms developed in HRP,” Robotics and Autonomous Systems, vol. 48, pp. 165–175, 2004. [95] J. Pratt, C. M. Chew, A. Torres, P. Dilworth, and G. Pratt, “Virtual model control: An intuitive approach for bipedal locomotion,” vol. 20, no. 2, pp. 129–143, 2001. [96] M. Wisse, G. Feliksdal, J. V. Frankenhuyzen, and B. Moyer, “Passive-based walking robot - Denise, a simple efficient and lightweight biped,” vol. 14, no. 2, pp. 52–62, 2007. [97] C. Chevallereau, G. Abba, Y. Aoustin, F. Plestan, E. R. Westervelt, C. C. de Wit, and J. W. Grizzle, “Rabbit: A testbed for advanced control theory,” IEEE Control Systems Magazine, vol. 23, pp. 57–79, 2003. [98] E. R. Westervelt, J. W. Grizzle, and D. E. Koditschek, “Hybrid zero dynamics of planar biped walkers,” IEEE Transactions on Automatic Control, vol. 48, pp. 42–56, 2003. 203 [99] L. Flynn, R. Jafari, A. Hellum, and R. Mukherjee, “An energy optimal gait for the MSU synthetic-wheel biped,” in Proc. of ASME Dynamic Systems and Control Conference, 2010. [100] S. Lohmeier, T. Buschmann, and H. Ulbrich, “System design and control of anthropomorphic walking robot LOLA,” IEEE/ASME Transactions on Mechatronics, vol. 14, no. 6, pp. 658 – 666, 2009. [101] F. Asano and Z. W. Luo, “Underactuated virtual passive dynamic walking with an upper body,” pp. 2441–2446, 2008. [102] M. Wisse, D. G. E. Hobbelen, and A. L. Schwab, “Adding an upper body to passive dynamic walking robots by means of a bisecting hip mechanism,” vol. 23, no. 1, pp. 112– 123, 2007. [103] M. Gomes and A. Ruina, “Walking model with no energy cost,” Physical Review E, vol. 83, no. 3, p. 032901, 2011. [104] B. Francis, “The linear multivariable regulator problem,” SIAM Journal on Control and Optimization, vol. 15, pp. 486–505, 1977. [105] B. Francis and W. Wonham, “The internal model principle for linear multivariable regulators,” Applied Mathematics and Optimization, vol. 2, no. 2, pp. 170–194, 1975. [106] E. Davison, “The robust control of a servomechanism problem for linear time-invariant multivariable systems,” Automatic Control, IEEE Transactions on, vol. 21, no. 1, pp. 25–34, 1976. [107] G. Bengtsson, “Output regulation and internal models–A frequency domain approach,” Automatica, vol. 13, no. 4, pp. 333–345, 1977. [108] A. Saberi, A. Stoorvogel, and Z. Lin, “Generalized output regulation for linear systems,” in American Control Conference, 1997. Proceedings of the 1997, vol. 6, pp. 3909–3914, IEEE, 2002. [109] R. Marino and P. Tomei, “Output regulation for linear minimum phase systems with unknown order exosystem,” Automatic Control, IEEE Transactions on, vol. 52, no. 10, pp. 2000–2005, 2007. [110] R. Hirschorn, “Invertibility of nonlinear control systems,” SIAM Journal on Control and Optimization, vol. 17, p. 289, 1979. 204 [111] S. Devasia, D. Chen, and B. Paden, “Nonlinear inversion-based output tracking,” Automatic Control, IEEE Transactions on, vol. 41, no. 7, pp. 930–942, 2002. [112] X. Wang and D. Chen, “Causal inversion of nonminimum phase systems,” in IEEE Conference on Decision and Control, vol. 1, pp. 73–78, IEEE; 1998, 2001. [113] S. Devasia, “Output tracking with nonhyperbolic and near nonhyperbolic internal dynamics: Helicopter hover control,” in American Control Conference, 1997. Proceedings of the 1997, vol. 3, pp. 1439–1446, IEEE, 1997. [114] S. Devasia, “Approximated stable inversion for nonlinear systems with nonhyperbolic internal dynamics,” Automatic Control, IEEE Transactions on, vol. 44, no. 7, pp. 1419– 1425, 2002. [115] J. Grizzle, M. Di Benedetto, and F. Lamnabhi-Lagarrigue, “Necessary conditions for asymptotic tracking in nonlinear systems,” Automatic Control, IEEE Transactions on, vol. 39, no. 9, pp. 1782–1794, 2002. [116] S. Gopalswamy and J. Hedrick, “Tracking nonlinear non-minimum phase systems using sliding control,” International Journal of Control, vol. 57, no. 5, pp. 1141–1158, 1993. [117] J. Hauser, S. Sastry, and P. Kokotovic, “Nonlinear control via approximate inputoutput linearization: The ball and beam example,” Automatic Control, IEEE Transactions on, vol. 37, no. 3, pp. 392–398, 1992. [118] J. Romano and S. Singh, “IO map inversion, zero dynamics and flight control,” Aerospace and Electronic Systems, IEEE Transactions on, vol. 26, no. 6, pp. 1022– 1029, 1990. [119] L. Qiu and E. Davison, “Performance limitations of non-minimum phase systems in the servomechanism problem,” Automatica, vol. 29, no. 2, pp. 337–349, 1993. [120] B. Francis, “The optimal linear-quadratic time-invariant regulator with cheap control,” Automatic Control, IEEE Transactions on, vol. 24, no. 4, pp. 616–621, 2002. [121] M. Seron, J. Braslavsky, P. Kokotovic, and D. Mayne, “Feedback limitations in nonlinear systems: From Bode integrals to cheap control,” Automatic Control, IEEE Transactions on, vol. 44, no. 4, pp. 829–833, 2002. [122] O. Toker, J. Chen, and L. Qiu, “Tracking performance limitations in LTI multivariable discrete-time systems,” Circuits and Systems I: Fundamental Theory and Applications, IEEE Transactions on, vol. 49, no. 5, pp. 657–670, 2002. 205 [123] M. Tomizuka, “Zero phase error tracking algorithm for digital control.,” J. Dyn. Syst. Meas. Control., vol. 109, no. 1, pp. 65–68, 1987. [124] E. Gross, M. Tomizuka, and W. Messner, “Cancellation of discrete time unstable zeros by feedforward control,” Journal of dynamic systems, measurement, and control, vol. 116, no. 1, pp. 33–38, 1994. [125] B. Haack and M. Tomizuka, “The effect of adding zeroes to feedforward controllers,” Journal of dynamic systems, measurement, and control, vol. 113, p. 6, 1991. [126] P. Kabamba, “Control of linear systems using generalized sampled-data hold functions,” Automatic Control, IEEE Transactions on, vol. 32, no. 9, pp. 772–783, 1987. [127] T. Mita, Y. Chida, Y. Kaku, and H. Numasato, “Two-delay robust digital control and its applications-avoiding the problem on unstable limiting zeros,” Automatic Control, IEEE Transactions on, vol. 35, no. 8, pp. 962–970, 1990. [128] H. Fujimoto, Y. Hori, and A. Kawamura, “Perfect tracking control based on multirate feedforward control with generalized sampling periods,” Industrial Electronics, IEEE Transactions on, vol. 48, no. 3, pp. 636–644, 2002. [129] H. Toliyat and G. Kliman, Handbook of electric motors, vol. 120. CRC Press, 2004. [130] R. Jafari, F. Mathis, and R. Mukherjee, “Swing-up control of the acrobot: An impulsemomentum approach,” in American Control Conference (ACC), 2011, pp. 262–267, IEEE, 2011. [131] W. Haddad, V. Chellaboina, and S. Nersesov, Impulsive and hybrid dynamical systems: stability, dissipativity, and control. Princeton Univ Pr, 2006. [132] V. Chellaboina, S. Bhat, and W. Haddad, “An invariance principle for nonlinear hybrid and impulsive dynamical systems,” Nonlinear Analysis, vol. 53, no. 3-4, pp. 527–550, 2003. [133] R. Sanfelice, R. Goebel, and A. Teel, “Invariance principles for hybrid systems with connections to detectability and asymptotic stability,” Automatic Control, IEEE Transactions on, vol. 52, no. 12, pp. 2282–2297, 2007. [134] F. Mathis, R. Jafari, and R. Mukherjee, “Efficient swing-up of the acrobot using continuous torque and impulsive braking,” in American Control Conference (ACC), 2011, pp. 268–273, IEEE, 2011. 206 [135] J. Pratt, J. Carff, S. Drakunov, and A. Goswami, “Capture point: A step toward humanoid push recovery,” in EEE-RAS/RSJ International Conference on Humanoid Robots, December, Genoa, Italy 2006. [136] K. Yin, K. Loken, and M. van de Panne, “Simbicon: Simple biped locomotion control,” ACM Trans. Graph., vol. 26, no. 3, p. 105, 2007. [137] B. Stephens, “Integral control of humanoid balance,” 2007. [138] L. Flynn, “Active synthetic wheel prismatic joint biped,” Master’s thesis, Michigan State University, 2009. [139] T. Hagiwara, T. Yuasa, and M. Araki, “Stability of the limiting zeros of sampled-data systems with zero-and first-order holds,” International Journal of Control, vol. 58, no. 6, pp. 1325–1346, 1993. [140] Y. Hayakawa, S. Hosoe, and M. Ito, “On the limiting zeros of sampled multivariable systems,” Systems & Control Letters, vol. 2, no. 5, pp. 292–300, 1983. [141] M. Ishitobi, “A stability condition of zeros of sampled multivariable systems,” Automatic Control, IEEE Transactions on, vol. 45, no. 2, pp. 295–299, 2000. [142] K. Astroem, P. Hagander, and J. Sternby, “Zeros of sampled systems.,” Automatica, vol. 20, no. 1, pp. 31–38, 1984. [143] T. Das and R. Mukherjee, “Shared-sensing and control using reversible transducers,” Control Systems Technology, IEEE Transactions on, vol. 17, no. 1, pp. 242–248, 2009. [144] K. Ogata, Modern control engineering. Prentice Hall, 2009. [145] W. R. Wade, An introduction to analysis. Prentice Hall, 2000. [146] H. Khalil, Nonlinear systems, vol. 3. Prentice hall Englewood Cliffs, NJ, 2002. [147] G. Klambauer, Real analysis. American Elsevier, 1973. 207