LPV MODELING AND ADAPTIVE MPC OF A TILT-ROTOR AIRCRAFT By Shen Qu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Doctor of Philosophy 2021 ABSTRACT LPV MODELING AND ADAPTIVE MPC OF A TILT-ROTOR AIRCRAFT By Shen Qu This dissertation studies the modeling and control of an urban air mobility (UAM) electric vertical take-off and landing (eVTOL) aircraft with tilt-rotors, covering a brief UAM survey, analytical nonlinear flight modeling, control-oriented linear parameter varying (LPV) modeling, adaptive model predictive control (MPC) design and optimization with simulation validation. The core UAM aircraft control framework adopts the adaptive MPC with dynamic reference compensation (DRC) based on developed LPV models with certain modifications and improvements to adaptive MPC due to specific UAM aircraft characteristics. Both hovering and tilt transition control cases are studied. Furthermore, motor failure cases are also investigated under hovering and tilting flight conditions, where motor failure is introduced by limiting motor output power to the corresponding level of failure and the adaptive MPC design is conducted by modifying associated hard constraints. Note that failure recovery control is used to study safety redundancy of UAM aircrafts. UAM has attracted tremendous attention from both transportation industry and academia com- munity. It is viewed as the promising future of next-generation transportation and is expected to be a quiet, fast, clean, efficient, and safe point-to-point transportation. The most promising and popular UAM concept is eVTOL aircrafts. Equipped with a distributed electric propulsion (DEP) system, eVTOL aircrafts use multiple electric motors and propellers to provide the required lift force for vertical takeoff and landing, and tilt-rotors (or extra push rotors) are used to cruise horizontally. The fixed-wing design allows for improving cruise efficiency during the flight. How- ever, aeromechanics and flight dynamics of tiltrotor eVTOL aircrafts are complicated due to the increased system complexity. Thus, a comprehensive study through dynamic modeling and control of eVTOLs is necessary, which will form a foundation for future dynamic analysis, controller design and optimization. This dissertation firstly introduces an analytical flight dynamic formulation for UAM eVTOL aircrafts that feature tiltrotors for fixed-wing level flight, vertical takeoff and landing. In this ana- lytical formulation, a nonlinear rigid-body dynamic model is developed by incorporating multiple tiltrotor dynamics, their gyroscopic and inertial coupling effects. A quasi-steady aerodynamic formulation is implemented to calculate the aerodynamic loads on all lifting surfaces. In addition to the conventional control surfaces for fixed-wing aircrafts, such as elevator, both tilt angle and individual rotor rotational speed are considered as control inputs in the formulation of nonlinear flight dynamic model. These nonlinear dynamics is then linearized with respect to a set of trimmed flight conditions of interest (or a tilt operational trajectory) to render the corresponding linear time-invariant state-space models used to create an LPV model as a function of flight condition to be used for adaptive MPC control design. Based on the above analytical model, the aircraft control strategies for hovering and tilt transition operations are investigated using adaptive MPC based on the developed LPV models. To model the aircraft dynamics variation, a family of linearized time-invariant models are obtained by trimming the nonlinear aircraft model at multiple equilibrium conditions or linearizing along a predetermined flight trajectory, and the associated LPV model is obtained by linking the linear time-invariant models using scheduling parameter(s), where the system variation is modeled by the measured scheduling parameter(s) under current flight condition. The proposed adaptive MPC strategy is developed to optimize the system output performance, including the rigid-body aircraft velocity and altitude, and solved by using quadratic programming optimization with dynamic reference compensation (DRC), subject to a set of time-varying hard constraints representing the available propeller acceleration calculated based on the available motor power. The adaptive MPC simulations are conducted based on both developed LPV and nonlinear rigid-body models, and the simulation results demonstrate that the designed adaptive MPC controllers are able to achieve desired closed-loop performance for the target UAM aircraft. Copyright by SHEN QU 2021 ACKNOWLEDGEMENTS These years of my Ph.D. study experience is a meaningful and unforgettable time during my life. Since I joined the Michigan State University in 2016, I have been feeling so lucky for having this valuable opportunity to study at this beautiful university, with the supports and help received from many people. Without it, it would be impossible to complete this dissertation and my Ph.D. study. Firstly, I would like to express my most appreciation to my advisor, Professor Guoming (George) Zhu. As a new Ph.D. student, I was taught to complete my course work with a good GPA, to get familiar with the experimental research tools in the lab, and finally to conduct academic research and publish the associated research work. My advisor is a model for me not only in academic aspect for his outstanding research ability in control theory and applications but also as a person how he deals with his family, friends, and students with kindness and selflessness. I also want to express my appreciation to my advisor’s wife, Mrs. Hongli Zhu, who always provides us good career advice, and also organizes the annual lab gatherings, making us feel warm and being cared. Also, I would like to thank Professor Weihua Su from University of Alabama and Professor Sean Swei from Khalifa University, who are both erudite researchers in aerospace engineering. They always provide me a lot of supports on my UAM (urban air mobility) research, which helps me to conduct my research on a correct and efficient path. From these meetings with Professors Su and Swei, along with Professor Zhu, I wanted to work with them forever since their minds and ideas are so valuable and interesting but I can only learn from them with limited project duration. Besides, I would like to thank Professors Ranjan Mukherjee, Vaibhav Srivastava, and Zhaojian Li for being my guidance committee members. They always provide me instructions and help with fast responses under their busy work schedule. Also, their courses, lectures and talks broaden my knowledge in control theories, as they are all specialized in certain aspects with outstanding controls knowledge and achievements. v I experienced a wonderful time working and having fun with all the peers in our ERC lab including Prof. Tianyi HE, Prof. Liang Liu, Prof. Ali Al-jiboory, Dr. Ruitao Song, Dr. Yifan Men, Prof. Ali Alhajjar, Dr. Yingxu Wang, Dr. Ruixue Li, Dr. Wenpeng Wei, Dr. Anuj Pal, Mr. Jian Tang, Dr. Huan Li, Dr. Chengsheng Miao, Ms. Dawei Hu, Mr. Yu He, Mr. Corey Gamache and Mr. Lingyun Hua. I really appreciate their support during my Ph.D. study at MSU. In my daily life, my deep gratitude is given to my companion, Agumon (Gu Meng), a lovely stuffed Digimon I bought in 2016 who stays with me all the time. Also, I would like to express my gratitude to my best friends in China, including Ms. Yan Yan (Cheng Cheng), Mr. Chengye Ma (Ma Dia), and many others. Although we can only meet a few times during these years, our friendships remains strong, which are extremely precious in modern society nowadays. I also had a great time enjoying life and doing hobbies like snowboarding, biking, car racing with my friends including Dr. Di Zhang, Dr. Na Pang, Mr. Zhan Liu, Dr. Jialin Liu, Mr. Jiamu Hu, Mr. Yiheng Zhang and many others. Lastly, I would like to express my most gratitude to my parents, Professors Yirong Ou and Linyan Qu. They always give me the best supports they can, allow me to study in an area I am interested in, and cultivate many hobbies in my childhood. TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x KEY TO SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Introduction of urban air mobility (UAM) . . . . . . . . . . . . . . . . . . . . . . 1 1.2 LPV Modeling and Adaptive MPC . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Linear parameter varying modeling . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Model predictive control . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.3 Adaptive MPC based on an LPV system model . . . . . . . . . . . . . . . 4 1.3 Motivations of adaptive LPV-MPC control for UAMs . . . . . . . . . . . . . . . . 5 1.4 Organization of dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 CHAPTER 2 LPV MODELING AND ADAPTIVE LPV-MPC . . . . . . . . . . . . . . . 7 2.1 LPV model development based on a nonlinear model . . . . . . . . . . . . . . . . 7 2.2 Modeling of UAM aircraft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1 General theoretical formulation of tiltrotor aircraft . . . . . . . . . . . . . . 8 2.2.2 Linearized equation and state-space equation . . . . . . . . . . . . . . . . 10 2.2.3 Dynamic formulation validation . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.4 Specific information of the tilt-rotor aircraft . . . . . . . . . . . . . . . . . 14 2.3 Demonstration example: LPV model of the eVTOL aircraft during transition . . . . 16 2.3.1 Flight state generation based on tilt angle . . . . . . . . . . . . . . . . . . 16 2.3.2 LPV modeling during transition process . . . . . . . . . . . . . . . . . . . 17 2.4 Modes alignment and LPV model error evaluation . . . . . . . . . . . . . . . . . . 18 2.4.1 Coordinate transformation . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4.2 System norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4.3 σ-shifted norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.4 Mode alignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4.5 LPV model accuracy evaluation . . . . . . . . . . . . . . . . . . . . . . . 22 2.5 Adaptive LPV-MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5.2 Model predictive control . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5.3 MPC with feedforward control information . . . . . . . . . . . . . . . . . 26 2.5.4 Discrete affine LPV system with equilibrium condition variation . . . . . . 27 2.5.5 Adaptive LPV-MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 CHAPTER 3 UAM HOVERING FAILURE CONTROL STUDY . . . . . . . . . . . . . . 30 3.1 LPV model of a hovering aircraft with single motor failure scheduled by motor speed 30 3.2 Formulation of adaptive LPV-MPC under hovering motor failure . . . . . . . . . . 34 3.2.1 Motor failure simulation with input constraint variation . . . . . . . . . . . 34 vii 3.2.2 Dynamic reference compensation . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Case study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3.1 Rotor failure simulation architecture . . . . . . . . . . . . . . . . . . . . . 37 3.3.2 Simulation cases study based on LPV aircraft model . . . . . . . . . . . . 38 3.3.2.1 Hovering failure Case A . . . . . . . . . . . . . . . . . . . . . . 38 3.3.2.2 Hovering failure Case B . . . . . . . . . . . . . . . . . . . . . . 43 3.3.2.3 Hovering failure Case C . . . . . . . . . . . . . . . . . . . . . . 46 3.3.2.4 Behaviors comparison of hovering cases under disturbance . . . . 52 3.3.3 Simulation cases study based on analytical aircraft model . . . . . . . . . . 62 CHAPTER 4 UAM TRANSITION/TILTING OPTIMIZATION AND CONTROL . . . . . 73 4.1 Steady-state dynamics analyzation and preliminary control study . . . . . . . . . . 75 4.2 Reference tilting trajectory generation and optimization . . . . . . . . . . . . . . . 82 4.2.1 Simplified aircraft force analysis model . . . . . . . . . . . . . . . . . . . 82 4.2.2 Tilt transition optimization . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3 Control design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.3.1 Linearization along a nominal trajectory . . . . . . . . . . . . . . . . . . . 87 4.3.2 Affine LPV model along the nominal trajectory . . . . . . . . . . . . . . . 88 4.3.3 Adaptive MPC of aircraft transition . . . . . . . . . . . . . . . . . . . . . 90 4.3.4 Adaptive MPC formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.4 Simulation Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 CHAPTER 5 TILT TRANSITION WITH MOTOR FAILURE . . . . . . . . . . . . . . . 101 5.1 2D LPV model of aircraft transition with motor failure . . . . . . . . . . . . . . . 101 5.1.1 Model error evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.1.2 Adaptive MPC formulation for tilt transition with rotor failure . . . . . . . 103 5.2 Simulation results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2.1 Transition with front-right motor failure . . . . . . . . . . . . . . . . . . . 105 5.2.2 Transition with mid-right motor failure . . . . . . . . . . . . . . . . . . . . 107 5.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 CHAPTER 6 CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 110 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 viii LIST OF TABLES Table 2.1: Inertial properties of UAM aircraft . . . . . . . . . . . . . . . . . . . . . . . . . 15 Table 2.2: Aerodynamic properties of wing and v-tail of UAM aircraft . . . . . . . . . . . . 15 Table 3.1: Equilibrium condition of hovering failure . . . . . . . . . . . . . . . . . . . . . 32 Table 3.2: System input and output definition . . . . . . . . . . . . . . . . . . . . . . . . . 34 Table 3.3: Equilibrium condition of hovering failure - Case A . . . . . . . . . . . . . . . . 38 Table 3.4: Adaptive MPC design specification parameters . . . . . . . . . . . . . . . . . . 40 Table 3.5: Equilibrium condition of hovering failure - Case B . . . . . . . . . . . . . . . . 46 Table 3.6: MPC controller spec parameters for hovering - Case C . . . . . . . . . . . . . . 47 Table 3.7: Condition definition of comparison cases . . . . . . . . . . . . . . . . . . . . . 52 Table 3.8: Adaptive MPC controller spec parameters . . . . . . . . . . . . . . . . . . . . . 62 Table 3.9: RMSE of simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Table 4.1: Equilibrium condition of tilting . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Table 4.2: Information of aircraft modes . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Table 4.3: System input and output definition . . . . . . . . . . . . . . . . . . . . . . . . . 79 Table 4.4: Aircraft parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Table 4.5: Adaptive MPC design specification parameters . . . . . . . . . . . . . . . . . . 94 Table 4.6: Performance evaluation of simulation results . . . . . . . . . . . . . . . . . . . 100 ix LIST OF FIGURES Figure 1.1: Roadmap of this dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Figure 2.1: Global and body reference frames of a rigid-body tiltrotor aircraft . . . . . . . . 9 Figure 2.2: A free rigid-body with a tiltrotor . . . . . . . . . . . . . . . . . . . . . . . . . 12 Figure 2.3: Structure of tiltrotor model validation process . . . . . . . . . . . . . . . . . . 12 Figure 2.4: Rotor speed under the actuation torques . . . . . . . . . . . . . . . . . . . . . . 13 Figure 2.5: Body linear velocities of validation model . . . . . . . . . . . . . . . . . . . . 13 Figure 2.6: Body angular velocities of validation model . . . . . . . . . . . . . . . . . . . 14 Figure 2.7: Geometry of sample tiltrotor UAM aircraft (top view) . . . . . . . . . . . . . . 16 Figure 2.8: LPV model accuracy evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 2.9: Structure of adaptive MPC controller . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 2.10: Equilibrium condition variation . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Figure 3.1: Coordinate definition of UAM aircraft . . . . . . . . . . . . . . . . . . . . . . 33 Figure 3.2: Motor failure simulation structure . . . . . . . . . . . . . . . . . . . . . . . . . 37 Figure 3.3: Open-loop linear speed response of hovering failure - Case A . . . . . . . . . . 39 Figure 3.4: Open-loop angular position response of hovering failure - Case A . . . . . . . . 39 Figure 3.5: Linear speed of hovering failure - Case A . . . . . . . . . . . . . . . . . . . . . 41 Figure 3.6: Angular speed of hovering failure - Case A . . . . . . . . . . . . . . . . . . . . 41 Figure 3.7: Angular position of hovering failure - Case A . . . . . . . . . . . . . . . . . . 42 Figure 3.8: Propeller speed of hovering failure - Case A . . . . . . . . . . . . . . . . . . . 42 Figure 3.9: Motor power at hovering failure simulation - Case B . . . . . . . . . . . . . . . 43 Figure 3.10: Linear speed of hovering failure - Case B . . . . . . . . . . . . . . . . . . . . . 44 x Figure 3.11: Angular position of hovering failure - Case B . . . . . . . . . . . . . . . . . . . 44 Figure 3.12: Propeller speed of hovering failure - Case B . . . . . . . . . . . . . . . . . . . 45 Figure 3.13: Propeller power of hovering failure - Case B . . . . . . . . . . . . . . . . . . . 45 Figure 3.14: Linear speed of hovering failure - Case C . . . . . . . . . . . . . . . . . . . . . 48 Figure 3.15: Angular speed of hovering failure - Case C . . . . . . . . . . . . . . . . . . . . 48 Figure 3.16: Linear displacement of hovering failure - Case C . . . . . . . . . . . . . . . . . 49 Figure 3.17: Angular position of hovering failure - Case C . . . . . . . . . . . . . . . . . . . 49 Figure 3.18: Propeller speed of hovering failure - Case C . . . . . . . . . . . . . . . . . . . 50 Figure 3.19: Propeller power of hovering failure - Case C . . . . . . . . . . . . . . . . . . . 50 Figure 3.20: First set of Motor torque of hovering failure - Case C . . . . . . . . . . . . . . 51 Figure 3.21: Second set of Motor torque of hovering failure - Case C . . . . . . . . . . . . . 51 Figure 3.22: Linear speed comparison - Case C1 . . . . . . . . . . . . . . . . . . . . . . . . 53 Figure 3.23: Linear speed comparison - Case C2 . . . . . . . . . . . . . . . . . . . . . . . . 53 Figure 3.24: Linear speed comparison - Case C3 . . . . . . . . . . . . . . . . . . . . . . . . 54 Figure 3.25: Linear speed comparison - Case C4 . . . . . . . . . . . . . . . . . . . . . . . . 54 Figure 3.26: Linear displacement comparison - Case C1 . . . . . . . . . . . . . . . . . . . . 54 Figure 3.27: Linear displacement comparison - Case C2 . . . . . . . . . . . . . . . . . . . . 55 Figure 3.28: Linear displacement comparison - Case C3 . . . . . . . . . . . . . . . . . . . . 55 Figure 3.29: Linear displacement comparison - Case C4 . . . . . . . . . . . . . . . . . . . . 56 Figure 3.30: Angular speed comparison - Case C1 . . . . . . . . . . . . . . . . . . . . . . . 56 Figure 3.31: Angular speed comparison - Case C2 . . . . . . . . . . . . . . . . . . . . . . . 56 Figure 3.32: Angular speed comparison - Case C3 . . . . . . . . . . . . . . . . . . . . . . . 57 Figure 3.33: Angular speed comparison - Case C4 . . . . . . . . . . . . . . . . . . . . . . . 57 xi Figure 3.34: Angular position comparison - Case C1 . . . . . . . . . . . . . . . . . . . . . . 57 Figure 3.35: Angular position comparison - Case C2 . . . . . . . . . . . . . . . . . . . . . . 58 Figure 3.36: Angular position comparison - Case C3 . . . . . . . . . . . . . . . . . . . . . . 58 Figure 3.37: Angular position comparison - Case C4 . . . . . . . . . . . . . . . . . . . . . . 58 Figure 3.38: Propeller speed comparison - Case C1 . . . . . . . . . . . . . . . . . . . . . . 59 Figure 3.39: Propeller speed comparison - Case C2 . . . . . . . . . . . . . . . . . . . . . . 59 Figure 3.40: Propeller speed comparison - Case C3 . . . . . . . . . . . . . . . . . . . . . . 59 Figure 3.41: Propeller speed comparison - Case C4 . . . . . . . . . . . . . . . . . . . . . . 60 Figure 3.42: Motor power comparison - Case C1 . . . . . . . . . . . . . . . . . . . . . . . . 60 Figure 3.43: Motor power comparison - Case C2 . . . . . . . . . . . . . . . . . . . . . . . . 60 Figure 3.44: Motor power comparison - Case C3 . . . . . . . . . . . . . . . . . . . . . . . . 61 Figure 3.45: Motor power comparison - Case C4 . . . . . . . . . . . . . . . . . . . . . . . . 61 Figure 3.46: Vertical speed and acceleration of ideal failure case . . . . . . . . . . . . . . . 63 Figure 3.47: Linear speed of propeller #1 failure situations . . . . . . . . . . . . . . . . . . 64 Figure 3.48: Angular speed of propeller #1 failure situations . . . . . . . . . . . . . . . . . 64 Figure 3.49: Speed of propellers in Case #1 and #3 . . . . . . . . . . . . . . . . . . . . . . 66 Figure 3.50: Linear speed of propeller #4 failure situations . . . . . . . . . . . . . . . . . . 67 Figure 3.51: Angular speed of propeller #4 failure situations . . . . . . . . . . . . . . . . . 67 Figure 3.52: Speed of propellers in Case HF6 . . . . . . . . . . . . . . . . . . . . . . . . . 68 Figure 3.53: Linear speed comparison of AMPC and MPC . . . . . . . . . . . . . . . . . . 69 Figure 3.54: Angular speed comparison of AMPC and MPC . . . . . . . . . . . . . . . . . 69 Figure 3.55: Speed of propellers comparison of AMPC and MPC . . . . . . . . . . . . . . . 70 Figure 4.1: Root loci during transition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 xii Figure 4.2: Level flight speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Figure 4.3: Level flight acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Figure 4.4: Vertical speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Figure 4.5: Vertical displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Figure 4.6: Pitch rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 Figure 4.7: Pitch angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 Figure 4.8: Propeller speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 Figure 4.9: Geometry of UAM aircraft with dynamic force definition . . . . . . . . . . . . 82 Figure 4.10: Comparison of three reference trajectories during transition flight . . . . . . . . 87 Figure 4.11: Variation of aircraft dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Figure 4.12: Adaptive MPC architecture for transition control . . . . . . . . . . . . . . . . . 92 Figure 4.13: Simulation response - Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Figure 4.14: Control efforts - Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Figure 4.15: Simulation response - Cases 2 and 3 . . . . . . . . . . . . . . . . . . . . . . . 97 Figure 4.16: Control efforts - Cases 2 (left) and 3 (right) . . . . . . . . . . . . . . . . . . . . 97 Figure 4.17: Power and energy comparison for Cases 2 and 3 . . . . . . . . . . . . . . . . . 98 Figure 4.18: Simulation response and disturbance of Cases 4 and 5 . . . . . . . . . . . . . . 98 Figure 4.19: Control efforts - Cases 4 (left) and 5 (right) . . . . . . . . . . . . . . . . . . . . 99 Figure 5.1: Model error evaluation with front-right rotor failure. . . . . . . . . . . . . . . . 104 Figure 5.2: Model error evaluation with mid-left rotor failure. . . . . . . . . . . . . . . . . 104 Figure 5.3: Adaptive MPC based on 2D LPV model . . . . . . . . . . . . . . . . . . . . . 105 Figure 5.4: Simulated responses of tilt transition with front-right motor failure . . . . . . . 106 Figure 5.5: Control effort of tilt transition with front-right motor failure . . . . . . . . . . . 107 xiii Figure 5.6: Simulated responses of transition with mid-right motor failure . . . . . . . . . . 107 Figure 5.7: Control efforts of tilt-transition with mid-right motor failure . . . . . . . . . . . 108 Figure 6.1: A free rigid-body with a tilt-rotor . . . . . . . . . . . . . . . . . . . . . . . . . 114 Figure 6.2: Aerodynamic frame and load components . . . . . . . . . . . . . . . . . . . . 120 Figure 6.3: Rotor speed and thrust during nominal vertical takeoff (identical speed and thrust for all rotors) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Figure 6.4: Individual rotor speed and thrust if rotor 1 rpm is zero . . . . . . . . . . . . . . 125 Figure 6.5: Individual rotor speed and thrust if rotor 3 rpm is zero . . . . . . . . . . . . . . 125 Figure 6.6: Trimmed aircraft body pitch and ruddervator deflection during level flight transition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Figure 6.7: Trimmed rotor tilt angle during level flight transition . . . . . . . . . . . . . . . 127 Figure 6.8: Trimmed rotor spin rate during level flight transition . . . . . . . . . . . . . . . 127 xiv KEY TO SYMBOLS A state-space system matrix A rotor disk area a0 aerodynamic frames along lifting surface B state-space control input matrix B body-fixed frame b local frame C rotational transformation matrix CT rotor thrust coefficient d rotor mass radial location vector E aerodynamic influence matrices F rotor thrust e engine frame G global frame g gravitational acceleration vector H angular momentum I mass moment of inertia L rotor pylon length vector M BB Inertia matrix of full aircraft C BB damping matrix of full aircraft m mass N rotor speed p position vector q trim variable R load vector R rotor radius xv r rotor frame S trim cost function T kinetic energy T rotor torque U potential energy u system control input vector V volume, v translational velocity W ext external work x system state vector yÛ airfoil velocity component along a0y , zÛ airfoil velocity component along a0z , αÛ airfoil pitch velocity β rigid-body velocity vector Γ control input of rotor spin angle γ rotor spin angle vector δe , δa , δr deflection of elevator, aileron, and rudder 0 nonlinear equilibrium ζ quaternions θ angular position Ξ control input of rotor tilt angle ξ rotor tilt angle vector ρ air density, Ωζ matrix of quaternions differential equation ω angular velocity xvi CHAPTER 1 INTRODUCTION 1.1 Introduction of urban air mobility (UAM) Urban air mobility (UAM) has attracted tremendous attention from the transportation industry and academic community. It is viewed as the promising future of next-generation transportation and is expected to be a quiet, fast, clean, efficient, and safe point-to-point transportation. UAM aircrafts will, without doubt, reshape the current transportation system and extend it from two- dimension mobility on the ground surface to the three-dimension one in the air space. Not only traffic congestion in the urban area can be largely relieved, but also human’s activities and life will be entirely changed when UAM aircrafts get matured. Unlike traditional transportation systems, like railways, buses, taxis, bicycles, UAM aircrafts do not need infrastructures built up along travel routes. This feature makes it possible that air vehicles (aicrafts) are not limited by built-up routes, and air traffic is much more flexible and efficient. The desire to develop UAM aircrafts is not a new idea. However, recent technological advances in materials, battery technology, autonomous systems, electric propulsion systems, and so on make it feasible for commercial operations. The market of UAM aircrafts is huge to explore and is experiencing rapid boosting in its early stage. In the past five or six years, over 50 companies or research institutes released prototypes or concept products to the public. One of the big news is that, in 2016, DARPA funded Aurora Flight Sciences with 89.4 million to develop a vertical take-off and landing (VTOL) distributed hybrid-electric experimental plane to explore the UAM feasibility. UAM aircrafts could be the next revolutionary industry to drive the economy and technology advances. Since the complete concept of UAM is still under exploration, a consensus has not yet been achieved completely among the community. Different companies are following their design phi- losophy and guidelines in UAM aircraft development. The most promising and popular concept 1 of realizing UAMs is eVTOL (electric Vertical Take-Off and Landing) aircrafts. Equipped with distributed electric propulsion (DEP) system, VTOL (Vertical Take-Off and Landin) aircrafts use multiple electric motors and propellers to provide the lift force for vertical takeoff and landing. Tiltrotors (or extra push rotors) are used for cruising flight horizontally. The fixed-wing design allows for improved aircraft cruise efficiency. This eVTOL concept has been adopted by the most leading VTOL companies, including Boeing Aurora, Airbus A3, Joby Aviation, and Lilium. On the other hand, related areas are widely investigated, including market prediction [1], voyage planning [2, 3, 4, 5], power and energy modeling and optimization [6, 7, 8, 9, 10, 11], dynamics analysis [12] and failure study [13]. Compared to traditional helicopters, tiltrotor eVTOL aircrafts have not only an enhanced diversification of power distribution, including position and tilting freedom, but also a faster-driven response with the advantages of utilizing electric motors. However, aeromechanics and flight dynamics of tiltrotor eVTOL aircrafts are complicated due to the increased system control freedom. Thus, a comprehensive study and dynamic modeling of the eVTOL are necessary and can be used as a foundation of future dynamic analysis, control strategy development, and case study. A tilt-rotor aircraft is a complex system consisting of flexible and rigid structural components connected each other. When coupled with unsteady aerodynamics, the aeromechanics and flight dynamics of tiltrotor eVTOL aircrafts are highly nonlinear, non-periodic, and challenging to simulate [14, 15]. In rotorcraft aircraft studies, a systematic multibody dynamic simulation approach [16, 17, 18, 19] is usually applied to accurately model all the essential structural components. 1.2 LPV Modeling and Adaptive MPC 1.2.1 Linear parameter varying modeling Linear Parameter-Varying (LPV) system modeling method and its associated control theories have drawn great interest in recent years [20, 21, 22, 23, 24, 25, 26, 27]. The advantage of LPV modeling methodology is that it presents system nonlinear dynamics using a linear system architecture, where the nonlinear system dynamics is captured by the variation of linear system parameters dependent 2 on measured scheduling parameters. Consider an LPV system described by Û = A(ρ(t))x(t) + B(ρ(t))u(t)   x(t)   (1.1)   y(k) = C(ρ(t))x(t) + D(ρ(t))u(t)    where x(t), y(t) and u(t) denote the system state, output and input vectors, respectively. Instead of being constant, system matrices A, B, C and D are scheduled by a time-varying measurable vector ρ(t). Thus, the nonlinear system dynamics can be captured by system matrix variation scheduled by the real-time measured signals. 1.2.2 Model predictive control Model predictive control (MPC) is a widely used advanced control method that solves process control problems for the current time step using predicted future information. Based on the system dynamic model, the MPC calculates the control effort u(k) at current time step k by solving a finite horizon optimization problem from k to k + N − 1 with the control solution u(k + m) (m = 0, 1, ..., N − 1), and only the current control solution u(k) (or first a few controls) are used for closed-loop control. The whole optimization process repeats at the next control step (or after next a few control steps). Nowadays, the reference model for linear MPC is normally described by a linear discrete-time state-space model shown in Eq. (1.2) x(k + 1) = Ax(k) + Bu(k) (1.2) where x(k) and u(k) denote the state and control input vectors, respectively. Then, an optimal control problem is formulated based on the following objective function. N−1 Õ N−1 Õ J= min[xT (k + N)Px(k + N) + xT (k + m)Qx(k + m) + uT (k + m)Ru(k + m)] u m=0 m=0 (1.3) subject to system (1.2) and E x + Fu ≤ h where N denotes the length of prediction horizon; P ≥ 0, Q ≥ 0 and R > 0 are weighting matrices for penalizing the final states, states and control inputs, respectively; E and F are constrain 3 coefficient matrices, and h is a constraint vector. Thus, Eqs. (1.2) and (1.3) define a quadratic program problem with hard constraints in a finite horizon that can be solved using many existing MPC design algorithms. 1.2.3 Adaptive MPC based on an LPV system model For nonlinear systems, linear MPC cannot lead to good closed-loop system performsnce since a simple linearized model cannot describe the nonlinear dynamics accurately. Nonlinear MPC theories have been studied in recent years for improving control performance in many aspects. However, the computational cost of nonlinear MPC strategies are always very or extreme high due to the nonlinear optimization process. Note that to make the linear MPC real-time implementable, a linear time-invariant (LTI) system model is often desired for the MPC design. With the development of adaptive MPC [28, 29] and linear parameter-varying (LPV) modeling method [30, 27, 31], adaptive MPC designed based on an LPV system model has been widely used for parameter-varying systems, which is aimed to improve the closed-loop system performance [32, 33]. Note that an LPV model represents a nonlinear system using linear differential equations with varying parameters as functions of a set of measurements called scheduling parameters. A discrete-time LPV system model is shown in Eq. (1.4)  x(k + 1) = A(ρ(k))x(k) + B(ρ(k))u(k)    (1.4)  y(k) = C(ρ(k))x(k) + D(ρ(k))u(k)     where ρ is the scheduling parameter vector measured in real-time and used for control. On the one hand, compared with conventional MPC, at each time step, the dynamic model used for MPC design is adaptively varied based on the current system operational condition, which means that the change of nonlinear system dynamics will be captured in the LPV model and updated at each MPC design step. On the other hand, compared with the nonlinear MPC, the adaptive LPV-MPC reduces the computational cost significantly since the linear quadratic optimization process can still be utilized. 4 1.3 Motivations of adaptive LPV-MPC control for UAMs To control a highly nonlinear system, nonlinear MPC method is naturally adoppted to optimize the closed-loop system performance. In addition, compared to linear MPC, the nonlinear MPC design iteratively calculates the optimal control over the prediction horizon since the nonlinear optimization problem is often nonconvex, leading to fairly high computational cost. Thus, the adaptive LPV-MPC control method becomes a good candidate for nonlinear systems with signifi- cantly reduced optimization complexity. At each control step, the linear reference model is updated with the current scheduling parameter along with hard constraints, design weightings, and dynamic reference correspondingly. For the tilt-rotor aircraft control, the adaptive LPV-MPC has the following advantages. Firstly, its capability of handling a set of constraints [34, 35, 36] is essential since limited motor power is one of most crucial constraints to be considered. Secondly, comparing with traditional PID (proportional-integral-derivative) control requiring a huge amount of tuning effort to achieve sub- optimal performance at best, the adaptive LPV-MPC uses a specific set of design weightings to shape the target output performance, making it easy to tune. Last, comparing with the nonlinear MPC, the computational efficiency of adaptive LPV-MPC shows the potential capability of improv- ing the aircraft performance with significantly reduced computational resources and is suitable for real-time implementation. 1.4 Organization of dissertation In this dissertation, a study of the Urban Air Mobility (UAM) is proposed, including LPV modeling and MPC design for multiple flight operational modes. After introducing urban air mobility and the adaptive LPV-MPC in Chapter 1, the LPV modeling method and adaptive LPV- MPC are discussed in Chapter 2 with detailed analytical nonlinear model of tilt-rotor aircraft in the Appendix A. Also, a systematic procedure of forming LPV model for adaptive MPC design based on a set of trimmed (linearized) LTI models at multiple equilibrium conditions from nonlinear model is proposed and demonstrated. Investigation of UAM hovering failure modeling and control 5 is discussed in Chapters 3. In Chapter 4, a preliminary steady state transition study is introduced at the beginning, followed with the transition optimization for minimal energy utilization during the transition based on the LPV model created using trajectory linearization method, and adaptive MPC-DRC designed for achieving desired transition performance. Chapter 5 studied the tilt transition subject to motor failure based on a two-parameter LPV model proposed. Lastly, the conclusions and future work are summarized in Chapter 6. The structure of this dissertation is shown in Figure 1.1 Figure 1.1: Roadmap of this dissertation 6 CHAPTER 2 LPV MODELING AND ADAPTIVE LPV-MPC 2.1 LPV model development based on a nonlinear model Consider a nonlinear system Û = f (x(t), u(t), ρ(t)) x(t) (2.1) y(t) = g(x(t), u(t), ρ(t)) where f and g are differentiable; x(t) ∈ Rm , y(t) ∈ Rn and u(t) ∈ R p denote the system state, output and input, respectively; and ρ(t) ∈ Rq is a measurable parameter vector called scheduling parameter vector that varies based on nonlinear system characteristics, and also, ρ is assumed to be continuously differentiable in time within its admissible trajectory set P and its variation rate is bounded to a known compact subset P. Û Thus, the magnitude and variation rate bounds of the scheduling parameter ρ can be expressed as ρ ∈ P = {ρimin ≤ ρi ≤ ρimax, i = 1, ..., q} (2.2) i i i ρÛ ∈ P = { ρÛ min ≤ ρÛ ≤ ρÛ max, i = 1, ..., q} Û To simplify the expression, time dependence notation t will be omitted through the remaining of dessertation. Assuming that there is a family of equilibrium points (x0 (ρ), u0 (ρ)) for the nonlinear system (4.11) such that f (x0 (ρ), u0 (ρ), ρ) = 0 (2.3) y0 (ρ) = g(x0 (ρ), u0 (ρ), ρ) Defining difference values of system state, input and output vectors as ∆x = x−x0 (ρ), ∆u = u−u0 (ρ) and ∆y = y − y0 (ρ), the linearization equation can be formulated in (2.4) based on Taylor expansion of f and g. ∆ xÛ = ∇ f |0 ∆x + ∇ f |0 ∆u + Ψ f (∆x, ∆u, ρ) (2.4) y = y0 (ρ) + ∇g|0 ∆x + ∇g|0 ∆u + Ψg (∆x, ∆u, ρ) 7 where (·)|0 means a function evaluated at the given nonlinear equilibrium condition; Ψ f and Ψg denote the high order terms of Taylor expansion. Thus, the linearized nonlinear system (2.1) can be formulated below.  ∆ xÛ = A(ρ)∆x + B(ρ)∆u + N(ρ) ρÛ + Ψ f (∆x, ∆u, ρ)    (2.5)   ∆y = C(ρ)∆x + D(ρ)∆u + Ψg (∆x, ∆u, ρ)    Noting that term N(ρ) ρÛ is based on the existing dependence on ρ for linearization and assumed to be N(ρ) ρÛ ≈ 0. Also, higher order terms of the Tayler expansion Ψ f and Ψg are assumed to be negligible and the LPV form of the original system is defined by  ∆ xÛ = A(ρ)∆x + B(ρ)∆u    (2.6)   ∆y = C(ρ)∆x + D(ρ)∆u    with its discrete form as below with sample period T.  ∆x(k + 1) = Ad (ρ)∆x(k) + Bd (ρ)∆u(k)    (2.7)  ∆y(k) = Cd (ρ)∆x(k) + Dd (ρ)∆u(k)     where Ad, Bd, Cd and Dd are the discrete-time system matrices obtained from A, B, C and D. 2.2 Modeling of UAM aircraft This section (Section 2.2) describes the nonlinear rigid-body model formulation for the target UAM aircraft and its linearization that are from [37] with rigid-body model developed by Dr. Weihua Su at University of Alabama who is a collaborator of this project. The modeling work is simplified and restated in this section and details can be found in Appendix , with the key to symbols stated in the key to symbols section. 2.2.1 General theoretical formulation of tiltrotor aircraft Firstly, the tiltrotor aircraft is considered as a fixed-wing aircraft (modeled as a rigid-body with fixed lifting surfaces) with tiltrotors; see Figure 2.1 for a specific two tiltrotor aircraft case. The 8 contributions of fixed-wing aircraft and tiltrotors are modeled individually and combined at the end. The aircraft’s rigid-body velocity is given by  pÛ B + ω B × p B       vB    β= = (2.8)       ωB  θÛ B           where v B and ω B are the translational and angular velocities of the rigid body, respectively; and pÛ B and θÛ B denote the inertial position of O B . Note that all kinematic quantities are resolved using the B frame unit vectors. G Bz pG/B dm pm/B B ωB Bx OB vB Gz Gy By pB OG Gx Figure 2.1: Global and body reference frames of a rigid-body tiltrotor aircraft By following the Hamilton’s Principle, the governing equation of the complete aircraft is found to be M BB (ξ)βÛ + C BB (β, ξ)β = R B (2.9) where the inertia, damping matrices and load vector are      mB I 3 mB p̃TG/B   mr I3 mr p̃TA/B  M BB =  +     m p̃  B G/B I B  m p̃   r A/B r I + m p̃ p̃T  r A/B A/B      (2.10)      mB ω̃ B m ω̃ p̃T   mr ω̃ B m ω̃ p̃T B B r B  G/B  +  A/B C BB =       mB p̃G/B ω̃ B ω̃ B I B   r A/B ω̃ B ω̃ B Ir + mr p̃ A/B ω̃ B p̃TA/B    m p̃      R B = Rgrav + Riner + Rrate + Rgyro + Rext 9 where I3 denotes the 3-by-3 identity matrix. The gravitational acceleration g is resolved in the B frame. Rgrav , Riner , Rrate , Rgyro and Rext denote the gravity load, inertial load, change of inertia, gyroscopic load, propulsive and aerodynamic loads, respectively, defined as follows. grav       F  I   I 3 3    grav (2.11) = =  mB g +      R  mr g      Mgrav    p̃G/B  p̃ (ξ)  A/B        iner       iner  F   I  3  iner   03×1   = = + (2.12)       R F   Miner    p̃ A/B   Mracc            rate  03×1   = (2.13)    R  Mrate         03×1   Rgyro = (2.14)     Mgyro      Rext = Rprop + Raero (2.15) where 0 0             prop               prop  F   I 3  Be   0   3  Be   = =  C (ξ) F + (2.16)             R   C (ξ) T    Mprop    p̃ A/B (ξ)       I   3       0 0                    E f f aero     aero = (2.17)    R  Em maero      with detailed derivation of these equations stated in Appendix . 2.2.2 Linearized equation and state-space equation Taylor expansion is applied to the governing equation for a given nonlinear equilibrium ( 0 ) and corresponding control input (u0 ) as below.  T  T  T T  T  0 = βT0 βÛ 0 ζT0 ζÛ 0 pG B,0 pÛ G B,0 n o (2.18) T T T T ÜT u0 = δ T δa,0 δr,0 Ξ0 Ξ0 Ξ0 Γ0 Γ0 T T T Û Ü Û e,0 10 After Taylor expansion with high order terms neglected and terms with the same variables grouped together, the linearized equation of motion can be obtained. ! ∂R B ∂C ∂R   BB B M BB | 0 − ∆βÛ + C BB | 0 + β − ∆β ∂ βÛ ∂β 0 0 ∂β 0 0 ∂R B ∂R B ∂R B ∂R B = ∆ζ + ∆δe + ∆δa + ∆δr ∂ζ 0 ∂δe 0 ∂δa 0 ∂δr 0 ∂R B ∂M BB Û ∂C BB   + − β − β ∆Ξ ∂Ξ 0 ∂Ξ 0 0 ∂Ξ 0 0 (2.19) ∂R B Û + ∂R B ∆Ξ Ü + ∂R B ∆ΓÛ + ∂R B ∆ΓÜ + ∆Ξ ∂Ξ 0 Û ∂Ξ 0 Ü ∂ ΓÛ 0 ∂ ΓÜ 0 1 ∂Ωζ (β) 1 ∆ζÛ = − ∆βζ0 − Ωζ (β0 )∆ζ 2 ∂β 0 2 ∂CGB (ζ)   h i G ∆pÛ B = ∆ζ 03 β0 + C (ζ0 ) 03 ∆β GB ∂ζ 0 where (·)| 0 means a function evaluated at the given nonlinear equilibrium. One can omit the ∆ sign in Eq. (2.19), and derive the state-space form of the linearized equation of motion as below. xÛ = Ax + Bu (2.20) where state variable and control input vectors are given by   T T x= β ζ T T pGB h iT (2.21) u = δTe δTa δrT ΞT Ξ ÛT ÜT Ξ ΓÛT ÜT Γ 2.2.3 Dynamic formulation validation In order to validate the developed formulation, a single tiltrotor is added to a rigid body, as shown in Figure 2.2. For simplicity, the origin of body frame B is at the rigid-body center of mass, i.e. pG/B = 0. The pivot point of the tiltrotor, Oe , is chosen to be the same as Ob , O B , and OG , i.e., the pivot point is the same as mass center and B frame origin. The pylon length L is set to be 0.5 m. This model is validated in the Matlab™-Simulink™ environment using Simulink™ Simscape™ Multibody™ toolbox with the procedure described in Figure 2.3. Note that the Simscape™ Multi- 11 bz by ξ ob pe A Bz bx oe γ pb By (a) full system Gz Gy vB pB OB ωB Gx OG Bx ez rz ez L ey A γ γ ex A bz ξ, ξ R d dm rx r pe oe ob by (b) tiltrotor (c) rotor Figure 2.2: A free rigid-body with a tiltrotor Torque Simulink Tilting/Spin Derived Actuation Multibody Motion Dynamic Model Model Body Motion Body Motion Matched? Figure 2.3: Structure of tiltrotor model validation process body™ toolbox uses rotor torque as the control input while the developed analytical formulation utilizes rotor kinematics as the control input. To validate the analytical formulation, predefined torques are applied to the single rotor model in the Simscape™ Multibody™ toolbox to actuate the rotor’s spin and tilt motions. The rotor kinematics is then recorded and applied to the analytical formulation as its actuation. Specifically, from 0 to 1 s, a constant torque of 300 Nm is applied to spin the rotor with an acceleration, followed by a constant torque of 5 Nm applied to the tilt rotor from 1 to 20 s. The system is initially set at rest. From the Simscape™ Multibody™ simulation, the resulting rotor spin and tilt motions are shown in Figure 2.4, which are fed into the analytical formulation as its inputs for the simulation. Finally, the body linear velocity components and its Euler angles are solved using the two models (Simscape™ Multibody™ and proposed analytical 12 Spin rate, Rad/s 150 100 50 0 0 5 10 15 20 Tilting rate, Rad/s 4 2 0 0 5 10 15 20 Time, s Figure 2.4: Rotor speed under the actuation torques Figure 2.5: Body linear velocities of validation model 13 Figure 2.6: Body angular velocities of validation model models) and are compared in Figures 2.5 and 2.6. It can be seen that two responses are matched very well, validating the proposed analytical tiltrotor formulation successfully. 2.2.4 Specific information of the tilt-rotor aircraft Figure 2.7 shows the geometry of a tilt-rotor aircraft model of interest, with 6 rotors at 90 deg orientation and the fuselage not shown. The dashed-dotted line indicates the reference line of lifting surfaces, which is 25% the chord from leading edge. Both main wing and v-tail are modeled as rigid lifting surfaces. Table 2.1 lists the inertial parameters of rigid-body, fixed-wing aircraft and the rotors, and Table 2.2 lists the aerodynamic parameters of the main wing and v-tail. The main wing contains an aileron, as indicated in Figure 2.7, which occupies 25% of the chord, running from 60% to 90% of the wing span. A ruddervator is defined on the v-tail, occupying 25% of the chord, running from 40% to 80% of the tail span. 14 Table 2.1: Inertial properties of UAM aircraft Value Unit Body mass, mB 2240.7276 kg Body moment of inertia, I B,x x 12000 kgm2 Body moment of inertia, I B,yy 9400 kgm2 Body moment of inertia, I B,zz 20000 kgm2 Rotor mass, mr 4.5454 kg Rotor moment of inertia, Ir,x e x 3.5 kgm2 Rotor moment of inertia, Ir,yy e 7.0 kgm2 Rotor moment of inertia, Ir,zz e 3.5 kgm2 Table 2.2: Aerodynamic properties of wing and v-tail of UAM aircraft Wing Tail Unit Airfoil NACA 0012 NACA 0012 — Ref. axis location (from L.E.) 25% 25% — Span 13.720 6.900 m Sweep angle −2.306 0 deg Dihedral angle 0 0 deg Chord (root/tip) 2.075/0.970 1.080/1.080 m Incidence angle (no twist) 3.1598 1.0626 deg 15 11 m 5.3 m 3.51 m 2 0.45125 m 1 5.3 m 4 3 4.71 m 6 5 Figure 2.7: Geometry of sample tiltrotor UAM aircraft (top view) 2.3 Demonstration example: LPV model of the eVTOL aircraft during tran- sition 2.3.1 Flight state generation based on tilt angle To model the aircraft tilt transition process in the LPV form, scheduling parameter vector ρ is a scalar defined as the rotor tilt angle ρ(t) ∈ [0, π2 ] rad, assuming the tilting motion of multiple tiltable rotors are synchronized. The transition process from hovering (helicopter mode) to level flight (airplane mode) can be achieved with both mid-rotors tilted from "tilt-up (ρ(t) = π2 rad)" position to "tilt-forward (ρ(t) = 0 rad)". For generating the reference flight state trajectory, the desired transition motion of aircraft can be formulated as Vre f = f (ρ(t)) (2.22) where Vre f denotes the reference horizontal flight speed. As for rotor tilt angle ρ(t), the scheduling 16 trajectory can be formulated as π ∫ t ρ(t) = + ρ(τ)dτ Û 2 t0 ∫ t (2.23) ρ(t) Û = ρ(τ)dτ Ü t0 where t0 denotes the start point of transition process (the flight is at helicopter mode at t0 ). Assuming that the flight pitch and roll angle, vertical and lateral speed are holding at zero during the transition process, the reference flight state can be generated at the selected tilt angle with solved equilibrium flight condition following steps. 1) Define rotor tilting curve ρ(t) during the transition period as t ∈ [t0, t1 ], where ρ(t0 ) = π2 rad and ρ(t1 ) = 0 rad. 2) Calculated the reference flight speed Vre f based on Eq. (2.22). 3) Define total number N of LTI models to be obtained during the tilt transition process and select the reference flight state (ρi, Vrei ) for i = 1, 2, ..., N. For example, for N = 4, the reference flight f state can be selected as ρ1 = π2 rad, π3 rad, π6 rad and ρ4 = 0 rad with the calculated reference [1,2,3,4] speed Vre f , respectively. 4) Calculate the control surface U0i (ex. propeller speed, elevator deflection, etc.) at each flight transition state (ρi, Vrei ) for i = 1, 2, ..., N. f 5) Linearize the nonlinear aircraft model at each equilibrium flight state (ρi, Vre i ) for i = 1, 2, ..., N f with system matrices Ai , Bi , C i , and Di obtained. 2.3.2 LPV modeling during transition process Consider an LPV system Û = A(ρ(t))∆x(t) + B(ρ(t))∆u(t)   ∆ x(t)   (2.24)   ∆y(t) = C(ρ(t))∆x(t) + D(ρ(t))∆u(t)    17 where system matrices are variated based on scheduling parameter ρ. Interpolation method of LPV system matrices is described as below using A(ρ) as a sample, where other system matrices and equilibrium control surfaces follows the same procedure as well. With the family of N equilibrium flight conditions obtained for ρ(t) = ρi (i = 1, 2, ..., N) at which the rotor tilt angle is exactly equal to one of the obtained equilibrium conditions, system matrices of the LPV model can be interpolated as the same as these of linearized system. Linear interpolation method has been applied to obtain the LPV model when the rotor tilt angle is between two adjacent LTI models as ρ(k) ∈ (ρi, ρi+1 ), i = 1, 2, ..., N − 1. In summary, interpolation of system matrices and control surfaces are stated as i ∀ρ(t) = ρi  A(ρ(t)) = A ,     (2.25)  i + (Ai+1 − Ai ) ρ(t) − ρ , ∀ρ(t) ∈ (ρi, ρi+1 ), i = 1, 2, ..., N − 1 i A(ρ(t)) = A   ρi+1 − ρi    In summary, the system matrices and equilibrium control surfaces are variated based on scheduling parameter ρ(t). 2.4 Modes alignment and LPV model error evaluation Although the LPV modeling method is able to capture the original nonlinear system dynamics, certain modeling errors still exist, which are mainly caused by linear interpolation of adjacent LTI models. In this section, a notion of σ shifted H2 and H∞ norms are introduced [38] to evaluate the LPV modeling accuracy. Also, the mode mismatch problem is discussed here. 2.4.1 Coordinate transformation Consider a collection of N linearized LTI system models obtained from the previous section (assuming Di = 0) below. Û = Ai ∆x(t) + Bi ∆u(t)   ∆ x(t)   i G : i = 1, 2, ..., N (2.26)   ∆y(t) = C i ∆x(t)    18 where Gi denotes the linearized LTI model at scheduling parameter ρ(t) = ρi , and x(t), y(t) and u(t) denote the system state, output and input vectors, respectively. Assuming that there are m modes in each LTI model, coordinate transformation of system matrices are defined below. A i = (Pi )−1 Ai Pi, B i = (Pi )−1 Bi, Ci = C i Pi (2.27) where Pi denotes the coordinate transformation matrix for system Gi . The Jordan canonical form of each LTI model Gi in (2.26) yields Û = A i ∆x(t) + B i ∆u(t)   ∆ x(t)   i G : i = 1, 2, ..., N (2.28)   ∆y(t) = Ci ∆x(t)    and the associated system matrices are partitioned as follows.  1 0 ... 0   i   i A B   1   0 Ai . . . ..   .     Bi  2  , B i =  .2  , Ci = C i C i . . . Cm h i Ai =  i . (2.29)     . . . . . . . . . 0   ..  1 2          0 . . . 0 Aim     i   Bm      Based on the converted Jordan canonical form of the system matrices, the transfer function for each mode in Gi is described by Gij = C ij (sI − Aij )−1 Bij , i = 1, 2, ..., N, j = 1, 2, ..., m, (2.30) and the transfer function of LTI system Gi can be calculated as the summation of the transfer functions of each mode as Õm Õm Gi = Gij (s) = C ij (sI − Aij )−1 Bij , i = 1, 2, ..., N. (2.31) j=1 j=1 2.4.2 System norms For a stable, proper single-input and single-output (SISO) system with a transfer function G(s), the H∞ norm is defined by kG(s)k ∞ = sup |G( jω)| (2.32) ω 19 which is equal to the maximum gain over all frequencies or the peak value in the Bode magnitude plot from the frequency response. k y(t)k kG(s)k ∞ = sup (2.33) u,0 ku(t)k The H∞ of a multi-input and multi-output (MIMO) system is defined as the maximum singular value σ of the transfer function matrix G( jω) over the entire frequency range that can be efficiently calculated following the methods stated in [39]. The H2 norm of a stable SISO system defined by strictly proper transfer function G(s) is defined by s 1 ∫ ∞ kG(s)k 2 = |G( jω)| 2 dω (2.34) 2π −∞ Based on the realization of G(s) as xÛ = Ax + Bu and y = C x, a deterministic interpretation of H2 norm is described in term of its impulse response g(t) below s∫ ∞ kG(s)k 2 = |g(t)| 2 dt = kg(t)k 2 (2.35) 0 where g(t) = Ce At B denotes the system impulse response. Last but not the least, define g(t) = [g1 (t), g2 (t), ..., gnu (t)], where gi (t) is the i-th impulse response with respect to the i-th input and nu is number of inputs, and the H2 of the MIMO system is defined by s ∫ ∞ kG(s)k 2 = trace gT (t)g(t)dt (2.36) 0 Since g(t) = Ce At B in Eq. (2.36), the square of the H2 norm is ∫ ∞ 2 T kG(s)k 2 = trace BT e A t CT Ce AtBdt 0∫ ∞ T (2.37) = traceBT e A t CT Ce AtdtB 0 Letting ∫ ∞ T P= e A t CT Ce Atdt (2.38) 0 leads to p kG(s)k 2 = traceBT PB (2.39) 20 where P is the solution of the following Lyapunov equation AT P + P A + CT C = 0 (2.40) Noting that trace(M N) = trace(N M) for two matrices with compatible dimensions, Eq. (2.37) can be rewritten as ∫ ∞ T kG(s)k 22 = trace Ce At BBT e A t CT dt (2.41) 0 and the H2 norm can also be expressed as q kG(s)k 2 = traceCPc CT (2.42) ∫∞ T where Pc = 0 e At BBT e A t dt is the solution of the following Lyapunov equation APc + Pc AT + BBT = 0 (2.43) 2.4.3 σ-shifted norm Noting that the H2 and H∞ are defined for stable system. However, for each LTI model in the LPV system, there may be some unstable modes. To handle this situation, a relative stability principle is introduced here by shifting the imaginary axis to the right by σ amount, which converts the original unstable modes to stable ones for σ-shifted H2 and H∞ norms denoted by Hσ−2 and Hσ−∞ for future model error evaluation. For a system mode Gij (s), the Hσ−2 is defined by 2 2 Gij (s) = Gij (s + σ) = trace(C ij Pij (σ)C iT j ) (2.44) σ−2 2 where Pij (σ) ≥ 0 is obtained by solving the following Lyapunov equation (Aij − σI)Pij (σ) + Pij (σ)(Aij − σI)T + Bij BiT (2.45) j =0 As for the σ-shifted H∞ of system Gij (s), the Hσ−∞ is defined by 2 Gij (s) = Gij (s + σ) (2.46) σ−∞ ∞ which can be obtained following the same procedure of the unshifted H∞ norm. 21 2.4.4 Mode alignment For an LPV system with a set of N LTI models, consider the adjacent two LTI models between Gi and Gi+1 , where i = 1, 2, ..., N − 1, the sequential order of block diagonal mode may not be aligned due to the variation of system behaviors from one linearized point to another. Then, there could be mode mismatch between Gij and Gi+1 j , leading to large modeling error when linearly interpolating between these two models. This mismatch can be identified using the Hσ−2 mode distance defined by M ij (s) = Gij (s) − Gi+1 k (s) σ−2 , i = 1, 2, ..., N − 1, j, k = 1, 2, ..., m (2.47) where the minimum case of (M ij ) presents the combination of j, k for the most similar mode dynamics of model Gij (s) in Gi+1 k (s). In other words, j , k → (M ij )min denotes the identification of mismatched modes from two LTI models, which needs to be realigned. 2.4.5 LPV model accuracy evaluation Consider the number of N LTI models trimmed from the gridded parameter space within the nonlinear dynamics variation envelop, the LPV model accuracy can be improved by increasing number of models N. In principle, the modeling error will become negligible when N is sufficiently large. For evaluating the model accuracy, one option is to calculate the eigenvalues of realized system matrix A, and analyzing the natural mode variation between every two adjacent models. However, the system matrix B is not considered in this case, which means the sensitivity variation of control efforts is ignored. Recalling the Hσ−2 norm introduced in the previous section and assuming that all the modes are well aligned, the model error between two LTI systems can be evaluated by E i,i+1 = Gi (s) − Gi+1 (s) , i = 1, 2, ..., N − 1, (2.48) σ−2 where the σ is properly selected to make sure that all Gi for i = 1, 2, ..., N are stable. Let  > 0 ( ∈ Rm ) denotes the predefined Hσ−2 bound of allowed LPV modeling error, an additional LTI model will be inserted between models Gi and Gi+1 only if the modeling error E i,i+1 exceeds the 22 giving bound, or E i,i+1 > . In summary, the process of achieving desired LPV model accuracy is shown in Figure 2.8. Figure 2.8: LPV model accuracy evaluation 2.5 Adaptive LPV-MPC 2.5.1 Motivation Compared with linear MPC, nonlinear MPC design is normally solved by an iterative solution process for the optimal control which could make it hard to find a feasible solution. Also, the computational cost could be increased respectively due to the problem complexity which often makes it impossible for real-time applications. On the other hand, to achieve the desired control performance for a nonlinear system, the linear MPC design is based on an LTI model for control optimization, which could limit the achievable performance due to the modeling error between the linear and nonlinear models for the physical system. In this section, the adaptive LPV-MPC method 23 was introduced that not only takes the advantage of linear QP solution efficiency but also considers the nonlinear system dynamics for reducing modeling error. 2.5.2 Model predictive control Note that ∆x(k) in Eq. (2.7) is described by ∆x(k + 1) = x(k) − x0 (k) , Let e(k) = ∆xre f − ∆x(k) denote the tracking error at current time step k, the adaptive MPC design [34, 35, 36] is to find the constrained optimal control law ∆u(k) over a given horizon N that minimizes the quadratic performance defined by " N N−1 # 1 Õ T Õ min e (k + m)Qe(k + m) + ∆uT (k + m)R∆u(k + m) ∆u(k),...,∆u(k+N−1) 2 (2.49) m=0 m=0 subject to G∆u(k + m) ≤ h , m = 0, ..., N − 1 where Q ≥ 0 and R > 0 are the weighting matrices used to penalize tracking error e(k) and control effort ∆u(k), respectively, and G a constraint matrix used to ensure that the control effort ∆u(k) stay within the prescribed constraint vector h. In addition, e(k + m) and ∆u(k + m) represent the predicted error and control input at the m-th time step, respectively. For a finite prediction horizon of N steps, the cost function and the constraint equation in Eq. (2.49) can be reorganized into a compact form below. 1 min [êT (k)Q̂ ê(k) + ∆ûT (k) R̂∆û(k)] û(k) 2 (2.50) subject to Ĝ∆û(k) ≤ ĥ 24 where        e(k)   ∆u(k)   h              e(k + 1)   ∆u(k + 1)   h ê(k) =   ∈ R(N+1)n x , ∆û(k) = , ĥ =  .  ∈ R Nnu       ..  ..  ..    .     .           e(k + N) ∆u(k + N − 1)        h       0 . . . 0  G 0 . . . 0      R        0 R . . . 0  0 G . . . 0 R̂ =  . ..  , Ĝ =  ∈ R Nnu ×Nnu ,     .. .. . . . .  .. .. . . . ...   . . .       0 0 . . . R  0 0 . . . G         0 . . . 0    Q    0 Q . . . 0  Q̂ =  . (N+1)n x ×Nn x . ..  ∈ R   .. .. . .  . . .   0 0 . . . Q     (2.51) Furthermore, the predicted tracking error (over a given horizon N) ê(k) can be described by ê(k) = Âe(k) + B̂∆û(k) , (2.52) where  I  0 0 ... 0           0 ... 0       A(ρ)   B(ρ)    2  =  A (ρ)  ∈ R(N+1)n x ×n x , B̂ = (N+1)n x ×Nnu       A(ρ)B(ρ)  B(ρ) ... 0  ∈ R  ..  .. .. ..      .  . . . 0         ...  N   N−1 (ρ)B(ρ) AN−2 (ρ)B(ρ)   A (ρ) A B(ρ)    As shown in Eq. (4.24), the predicted tracking error ê(k) can be determined based on the cur- rent tracking error e(k) and control vector ∆û(k). Therefore, the optimization problem can be 25 reformulated by utilizing the current information as follows. 1 T 1 min ∆û (k)( R̂ + B̂T Q̂ B̂)∆û(k) + eT (k) ÂT Q̂ B̂∆û(k) + eT (k) ÂT Q̂ Âe(k) ∆û(k) 2 2 (2.53) subject to Ĝ∆û(k) ≤ ĥ(k) Note that the optimal control law ∆û(k) obtained from above is also the optimal solution to the constrained optimization problem defined in (2.49). For real-time control at the current time step k given the measured or estimated tracking error e(k) = ∆xre f − ∆x(k), the minimization problem described in Eq. (4.25) can be solved using the quadratic programming (QP) solver in Matlab [40]. Instead of applying only the first control entry ∆u(k + 0) at current time step k and repeating the optimization process in next time step k + 1, a control horizon of Nc ≤ N is chosen so that the first Nc control entries [u(k + 0), u(k + 1), ..., u(k + Nc − 1)] are used between current time step k and the time step k + Nc − 1. Subsequently, the optimization process is repeated at time step k + Nc − 1 to solve for the next control effort [u(k + Nc ), u(k + Nc + 1), ..., u(k + 2Nc − 1)] [41]. For MPC of LPV systems, certain closed-loop system stability can be achieved based on the theories developed in [33], which could lead to very conservative control with high computational complexity. On the other hand, a reasonable selection of tuning parameters and a well-posed optimization problem usually deliver stable closed-loop performance. Benefiting from the linear interpolation of LTI models, the parameter variation of prediction model between each control design step is very small due to continuous system dynamics and small sample period. Thus, in this study, the MPC stability is not guaranteed in theory but demonstrated in simulations. Closed-loop stability will be part of future work. 2.5.3 MPC with feedforward control information For the system with multiple control inputs, there is always a possible situation that certain control inputs are controlled with predefined feedforward information. Due to the coupled system dynamics, directly isolating the predefined feedforward control from MPC may lead to undesired system 26 performance since the predictive control is not accurate without considering the isolated feedforward terms. Thus, it is necessary to implement the feedforward information into the MPC framework. With the advantage of constraint handling ability of MPC, the feedforward control information can be implemented into the MPC optimization process as u p − ε ≤ ∆u f (k) ≤ u p + ε (2.54) where ∆u f (k) ∈ ∆u(k) denotes the feedforward controls and u p is the precalculated control efforts for achieving certain system performance, where ε is chosen to be small enough so that ∆u f (k) ≈ u p . 2.5.4 Discrete affine LPV system with equilibrium condition variation With the set of LTI models obtained by trimming the nonlinear model under different equilibrium conditions, the obtained discrete-time LTI system models are stated in Eq. (2.55), where i presents the LTI model index i i  ∆x(k + 1) = A ∆x(k) + B ∆u(k)    (2.55)     i i ∆y(k) = C ∆x(k) + D ∆t(k)  LPV modeling method [38] is applied by linearly interpreting the state-space model between operational conditions using a varying parameter. Considering the equilibrium condition variation, Eq. (2.7) is extended to the following discrete affine LPV model stated in Eq. (2.56)  x(k + 1) = x0 (ρ(k + 1)) + A(ρ(k))(x(k) − x0 (ρ(k))) + B(ρ(k))(u(k) − u0 (ρ(k)))    (2.56)  y(k) = y0 (ρ(k)) + C(ρ(k))x(k) + D(ρ(k))u(k)     where ∆x, ∆u and ∆y in Eq. (2.7) are presented in the difference form of system states, inputs, and outputs (x(k), u(k) and y(k)) with its equilibrium conditions x0 (ρ(k)), u0 (ρ(k)) and y0 (ρ(k)), where variables x0 (ρ(k)), u0 (ρ(k)), and y0 (ρ(k)) are linearly interpolated equilibrium conditions of the adjacent LTI models following the same method as interpolating the system matrices (see example for A(ρ)). It is worth to mention that ρ(k + 1) in Eq. (2.56) cannot be directly obtained since it depends on the system response in the next time step. Thus, it is approximately substituted by ρ(k + 1) ≈ ρ(k), assuming that the equilibrium condition from the current step to the next is negligible with high sampling rate. 27 2.5.5 Adaptive LPV-MPC Figure 2.9: Structure of adaptive MPC controller The adaptive MPC architecture is shown in Figure 2.9, where ∆x and ∆u are the deviations of controlled states and inputs. Based on the real-time state reference signal ∆x re f and the system state feedback signal ∆x, the MPC generates a set of optimized control signal from ∆u(k) to ∆u(k + N − 1), where control signals within the control horizon from ∆u(k) to ∆u(k + Nc − 1) are selected as the optimal state-feedback control ∆u that will be combined with u0 (ρ) to form control signal u = ∆u + u0 (ρ) and ρ(t) is the scheduling parameter. In this study, all the system states are assumed to be measurable due to available aircraft sensor set, that is, y = x (or C = I). The nonlinear aircraft model presents the nonlinear flight dynamics with system states/outputs x measured. Then, the required control feedback signal ∆x is obtained in the form of ∆x = x − x0 (ρ) that is fed back to the MPC for closed-loop flight control, where x0 (ρ) is the nominal states. On the other hand, the basic principle of the model variation is presented in Figure 2.10. Noting that the change of x0 (ρ) and u0 (ρ) denote the changing of steady-state system equilibrium condition. Thus, for real-time simulations, actual state x(t) and input u(t) are constructed as x(t) = x0 (ρ) + ∆x and u(t) = u0 (ρ) + ∆u, where x0 (ρ) and u0 (ρ) are calculated based on the scheduling parameter ρ(t) with ∆x = x − x0 (ρ) and ∆u determined by measured state x and MPC control, respectively. 28 Transition from normal to failed (1) scheduling parameter 𝜌 = 𝜌(𝑡) variation 𝜌 = 𝜌(𝑡 + 1) 𝑥 𝑥( (𝜌) Trim 𝑢( (𝜌) 𝑥 𝑥( (𝜌) Trim 𝑢( (𝜌) + + conditions conditions + 𝑢 + 𝑢 Trimmed Trimmed ∆𝑥 linear model ∆𝑥 linear model at 𝜌(𝑡) at 𝜌(𝑡 + 1) MPC ∆𝑢 MPC ∆𝑢 controller controller Figure 2.10: Equilibrium condition variation UA-MSU-SS Page 1 29 CHAPTER 3 UAM HOVERING FAILURE CONTROL STUDY In this chapter, motor failure control of a six-rotor electric vertical take-off and landing (eVTOL) urban air mobility (UAM) aircraft is investigated using adaptive model predictive control (MPC) based on the linear parameter-varying (LPV) model developed using the nonlinear rigid-body aircraft model. For capturing the aircraft dynamics under motor failure conditions, a family of linearized models are obtained by trimming the nonlinear aircraft model at multiple equilibrium conditions and the LPV model is obtained by linking the linear models using the failed rotor speed, where the system transition from healthy to failure is modeled by a scheduling parameter calculated based on failed rotor speed caused by reduced motor peak power after failure. The proposed adaptive MPC is developed to optimize the system output performance, including the rigid-body aircraft velocity and altitude, by using quadratic programming optimization with reference compensation subject to a set of time-varying constraints representing the current available propeller acceleration calculated based on the motor power. Simulation study is conducted based on the developed LPV control design and original nonlinear rigid-body model, and the simulation results demonstrate that the designed adaptive MPC controller is able to recover and maintain the aircraft at desired stable condition after motor failure. 3.1 LPV model of a hovering aircraft with single motor failure scheduled by motor speed In this section, the aircraft dynamics under single motor failure conditions is modeled as an LPV system. The motor power P of a rotor spinning at speed n can be calculated by Eq. (3.1) P = C p ρ a n3 D 5 (3.1) where Cp , ρa and D denote the power coefficient, air density and propeller diameter, respectively. Assuming that the target motor power and speed are Pnh and nh , respectively, under normal aircraft hovering situation, when the motor failed at certain percentage p f , meaning the maximum avalable 30 power Pmax, f of failed motor equal to p f × Pnh , and in this case the rotor speed will reduce from nh to n f , where n f is corresponding to Pmax, f . Assuming that the available power limits of the rest normal rotors are high enough to compensate the loss of lift force due to failed motor and are able to hold the aircraft at hovering flight condition, a new hovering equilibrium control can be found based on the failed maximum motor power Pmax, f or speed n f . Assuming that the failed motor is operated at its maximum available speed due to failed motor power limit, a set of N equilibrium points maintaining the aircraft hovering under a single motor failure for a given set of N available failed-motor power levels, P1 = Pnh > P2 > ... > P N = Pmax, f can be found, the assoiciated maximum failed rotor speed ni, f can be calculated based on Eq. (3.1) and the rest healthy rotor speed vectors ni (i = 1, 2, ..., N) maintaining hovering condition can also be found. As a result, the steady-state state vectors x0i (i = 1, 2, ..., N) can be found, so does the associated equilibrium control ui0 (i = 1, 2, ..., N) by trimming the nonlinear model (A.1.42). With the obtained x0i and ui0 (i = 1, 2, ..., N), LTI model set (Ai, Bi, C i, Di ) (i = 1, 2, ..., N), defined in Eq. (2.55), can be obtained. In summary, the set of LTI hovering aircraft models with single motor failure can be formulated based on the following procedure. 1) Define the location of failed motor. 2) Define total number N of LTI models to be obtained during motor failure transition, along with the associated available failed motor power Pi for i = 1, 2, ..., N. 3) Calculate the failed motor speed nif based on available failed rotor power Pi for i = 1, 2, ..., N using Eq. (3.1). 4) Calculate healthy rotor speed vectors ni = [ni1, ni2, ..., nip−1 ]T (i = 1, 2, ..., N), where p is the total number of rotors (see Table 3.1), assuming that the failed motor is located at p to simplify the notation. 5) Calculate equilibrium condition equilibrium conditions x0i and ui0 (i = 1, 2, ..., N) based on solved rotor speed values ni, f and ni (i = 1, 2, ..., N). 31 Table 3.1: Equilibrium condition of hovering failure Available failed motor power P1 P2 ... PN Failed rotor speed (rad/s) n1f n2f ... n Nf Other rotors speed (rad/s) n1 n2 ... nN 6) Linearize the nonlinear aircraft model based on the equilibrium condition (rotor speed), with system matrices Ai, Bi, C i, Di obtained. Then, the LPV system model can be obtained by linearly interpreting the neighboring state- space model between operational conditions using a varying parameter calculated from the failed rotor speed. Recall the LPV model is a linear state-space model with its coefficient matrices as an affine function of a time-varying parameter vector called scheduling parameter vector. In this study, the scheduling parameter vector is a scalar represented by ρ and chosen to be the dynamic power of the failed rotor as ρ(t) = P(t) that is calculated based on the real-time measured failed rotor speed n(t) (see Eq. (3.2)) ρ(t) = P(t) = Cp ρa n3 (t)D5 (3.2) Then, the LPV system can be generated following the process stated in Section 2.3.2. The LTI state-space models are obtained by linearizing (trimming) the nonlinear rigid-body model at each motor failure levels Pi (i = 1, 2, ..., N). Noting that the trimmed model stated in the previous chapter cannot be directly used for control design since some of the elements in the input channel are coupled each other, which means that directly designing the controller based on the trimmed model may lead to unrealistic solutions. Among the elements of input channels, the propeller speed, tilt speed and tilt angle are not independent. In fact, they are related to their driven torques and associated resistant (aero-dynamic load, damping) torques. As a result, they need to be added to the system state vector. Noting that the trimmed dynamics of propeller spin and tilt acceleration can be regarded as the equivalent driven torques since they are obtained under the equilibrium condition. To unify the expression, they are named controlled acceleration in this 32 dissertation. Letting the propeller acceleration be n, Û the spin dynamics can be represented by Eq. (3.3) below. nÛ = −Kn2 + u (3.3) where K = Cp ρD5 /Ip is based on Eq. (3.1) and u is the controlled propeller spin acceleration. Linearizing Eq. (3.3) at given propeller speed n yields the linearized model stated in Eq. (3.4), where −2nK is the obtained coefficient. δ nÛ = −2nKδn + δu (3.4) Now we are ready to construct the control design model system matrix Ac , Bc below Ac (1 : 12, 1 : 12) = Asys Ac (1 : 12, 13 : 18) = Bsys (1 : 12, 22 : 27) (3.5) Ac (13 : 18, 13 : 18) = −2nK · I6×6 Bc (1 : 12, 1 : 6) = Bsys (1 : 12, 28 : 33) (3.6) Bc (13 : 18, 1 : 6) = I6 as the dimension of system matrices changed from Asys ∈ R12×12 and Bsys ∈ R12×33 to Ac ∈ R18×18 and Bc ∈ R18×6 . Noting that all the unmentioned elements in Ac and Bc are zeros. System matrix Cc is set to identify (that is, Cc = I18 ), assuming that all the system states are measurable, where aircraft body coordinates are defined as y, z, and x; see Figure 3.1. Figure 3.1: Coordinate definition of UAM aircraft Definitions of input and output channels are summarized in Table 3.2. 33 Table 3.2: System input and output definition System Input Order System State Order Front right/left Propeller acceleration u1 /u2 Rigid-body velocity x1 to x6 Mid right/left Propeller acceleration u3 /u4 Euler angle x7 to x9 Rear right/left Propeller acceleration u5 /u6 Body inertial position x10 to x12 Propeller speed x13 to x18 System Output Order System Output Order Lateral speed (Vx ) y1 Roll angle (θ y ) y7 Longitudinal speed (Vy ) y2 Pitch angle (θ x ) y8 Vertical speed (Vz ) y3 Yaw angle (−θ z ) y9 Pitch rate (θÛx ) y4 Lateral displacement (D x ) y10 Roll rate (θÛy ) y5 Longitudinal displacement (D y ) y11 Opposite value of yaw rate (θÛz ) y6 Vertical displacement (Dz ) y12 3.2 Formulation of adaptive LPV-MPC under hovering motor failure 3.2.1 Motor failure simulation with input constraint variation Utilizing the constraint handling ability of the MPC, a motor failure simulation structure is intro- duced in this section. Firstly, the required power P to spin a propeller at speed n can be calculated using Eq. (3.1) Letting the propeller-motor assembly rotational inertia be Ip and the propeller acceleration n, Û the propeller rotational dynamics can be formed in Eq. (3.7) −Cp ρa D5 n2 + T nÛ = (3.7) Ip where T is the propeller driven torque generated from the motor, and term −Cp ρa D5 n2 calculates the resistant torque at current spin speed n. By defining the resistant coefficient K = Cp ρa D5 /Ip and the equivalent motor driven acceleration u = T/Ip , the propeller rotational dynamics can be expressed by Eq. (3.8) due to Eq. (3.7). nÛ = −Kn2 + u (3.8) 34 Consider a failed motor that can only provide r(t)% of power compared to the power (Phn ) used in normal hovering condition, the maximum propeller driven torque Tmax can be calculated based on Eq. (3.9) Ph ∗ r(t)% Tmax = (3.9) n(t) where Phn is the motor power at hovering. After converting Eq. (3.9) to the same unit as that used for control input (dividing the equation by propeller-motor assembly inertia Ip ), Eq. (3.10) presents the real-time MPC control constraint, where umax = Tmax /Ip can be considered as the equivalent propeller acceleration resulted by the motor driven torque. Ph ∗ r(t)% umax = (3.10) Ip n(t) After motor failure occurs under hovering condition, the propeller acceleration n(t) Û becomes negative based on Eq. (3.8), where u < Kn2 (t) leads to the reduction of propeller speed. In summary, the control input constraints umax = f (n, r) used for the MPC design are real-time calculated as a function of the propeller speed n(t) and available power percentage r(t). 3.2.2 Dynamic reference compensation To stabilize the vehicle under hovering condition, the reference signal for each output is set to 0, which means that the control target is to hold the system at its equilibrium condition. It was found, during calibrating the adaptive MPC controller, that even though large penalization was used in the weighting matrix on certain system outputs (including linear speeds and angular positions), these outputs have very slow responses, resulting in large errors; and sometimes the system is not even stable. In theory, a sufficient long prediction horizon will allow the MPC to solve the optimal control considering all the weighted outputs. However, the computational cost will be too high to be practical. To achieve desired system performance with feasible computational cost, a dynamic reference compensation method was developed with details explained below. Note that since only the incre- mental propeller acceleration is controlled directly thorugh ∆u, the propeller speed (or propeller 35 thrust) responses are slow since significant part of control u is contributed by u0 that is desired for maintaining steady-state condition and has nothing to do with the MPC design. To stabilize the hovering operation, the angular speeds (such as pitch, roll, and yaw rates) are directly affected by the propeller thrusts. As for the longitudinal/lateral speed, in the hovering case without any tilt operation, it can be controlled by manipulating the flight pitch/roll angle to change the direction of total thrust vector formed by all the propellers. On the other side, the flight angular position control can be achieved by setting the reference angular speed based on the feedback angular position with calibrated negative gain. Last but not the least, consider the vertical drop which is always crucial for hovering control, it will be recovered and stabilized by setting the reference vertical speed correspondingly. In summary, instead of finding well-calibrated combination overall weighting matrices for all system outputs with relative long prediction horizon, the output weighting matrix used in the adaptive MPC control design for flight hovering control is selected to mainly focus on the angular and vertical speeds that are directly affected by the control inputs. As a result, instead of holding the reference of angular/vertical speeds at zero, the calibrated dynamic compensation matrix is formulated in Eq. (3.11).    y1     re f     y   0 0 0 0 0 a16   y   3    2  re f      y   0 a 0 a 0 0  y7     4   22 24  = (3.11)       yre f  a  5   31 0 a33 0 0 0     y8         re f     y   0 0 0 0 a45 0   y   6     9    y12      re f where yi and yi (i = 3, 4, 5, 6) denotes the system outputs and their reference signals, with the definitions stated in Table 4.3. Observed from the simulation results during the controller development, the system responses are not smooth nor satisfactory after using this structure. Thus, re f constraints have been placed on yi (i = 1, 2, 7, 8, 9, 12) and yi (i = 3, 4, 5, 6) in Eq. (3.11), and re f the calculated reference signal yi (i = 3, 4, 5, 6) are filtered using a low-pass first order transfer 36 re f re f re f re f function such that Ŷi (s) = H(s)Yi (s), where Ŷi (s) and Yi (s) are the Laplace transforms re f re f of the ŷi (filtered reference signal) and Yi . This improves the MPC performance significantly with a stable flight. 3.3 Case study Multiple hovering failure cases are studied in this chapter under different failure conditions and evaluated using both LPV and analytical nonlinear rigid-bod aircraft models, where the open-loop and closed-loop system responses, and required motor power consumption are analyzed. The aircraft body coordinates (x, y, z) are defined as x pointing at the aitcraft right direction, y at the front direction, and z at the up direction; see Figure 3.1 for reference. Adaptive MPC 3.3.1 Rotor failure simulation architecture The simulation structure for propeller motor failure is showed in Figure 3.2. Failure/ Propeller System dynamics Torque drop speed drop Varying Figure 3.2: Motor failure simulation structure Assuming that the flight starts with normal hovering, the motor failure occurs with a sudden motor output torque drop, causing a gradual propeller speed reduction due to the fact that the motor cannot provide the desired torque to maintain the target propeller speed. Meanwhile, as the failed propeller speed is used to calculated the scheduling parameter of system dynamics, the system LPV model parameters will also be changed correspondingly. UA-MSU-SS Page 1 37 3.3.2 Simulation cases study based on LPV aircraft model 3.3.2.1 Hovering failure Case A For hovering failure Case A, the front right propeller (see Figure 2.7 for #1 propeller location) was selected for the failure study based on three trimmed models of r(t) = 100% (no failure), 66% (34% failure), and 33% (67% failure) motor power available, where the equilibrium conditions for each model are described in Table 3.3. Table 3.3: Equilibrium condition of hovering failure - Case A Equilibrium condition 100% power 66% power 33% power Front right speed (rad/s) 100.7 87.7 69.6 Front left speed (rad/s) 100.7 108.6 116.8 Mid right speed (rad/s) 100.7 108.4 116.6 Mid left speed (rad/s) 100.7 101.0 101.6 Rear right speed (rad/s) 100.7 100.8 101.2 Rear left speed (rad/s) 100.7 96.5 90.9 Assuming that the front right motor could only provide 47.6% of normal hovering power, the open-loop vehicle dynamic responses are showed in Figures 3.3 and 3.4. As observed from the angular position responses, pitch and roll angle oscillations can be observed with a magnitude around 0.035 radians, indicating unstable linear speed drifts of Vx and Vy , while Vz remains at zero or no vertical speed change. This behavior can be attributed to the LPV model structure since it uses the failed propeller speed as scheduling parameter generated based on the current failed motor speed and multiple steady-state equilibrium conditions. The open-loop response simulation confirms that it is necessary to design a closed-loop controller to stabilize the vehicle since using the steady-state open-loop compensation cannot stabilize the aircraft after motor failure. The adaptive MPC controller was designed by tuning the output weighting matrix to penalize 38 Figure 3.3: Open-loop linear speed response of hovering failure - Case A Figure 3.4: Open-loop angular position response of hovering failure - Case A 39 both pitch and roll rates (two key unstable states observed in the open-loop simulation study), and the control design parameters are listed in Table 3.4. Table 3.4: Adaptive MPC design specification parameters Parameter Value Parameter Value Q[1,1] 1 Q[2,2] 1 Q[3,3] 10002 Q[4,4] 4002 Q[5,5] 4002 Q[6,6] 102 R[1:6,1:6] 0.12 ∗ I6 Step size 1 ms Prediction Horizon 5 steps Control Horizon 1 steps min ∆u[1:6] 2 (t) −Kn[1:6] max ∆u[1:6] 2 (t)+max u −Kn[1:6] [1:6] a15 −0.3 a22 0.01 a24 −1 a31 −0.002 a33 −1 a45 −1 s.t y[1,2] [−5, 5] s.t y[7,8] [−0.5, 0.5] re f s.t y[4,5,6] [−0.03.0.03] H(s) 1 0.2s+1 The closed-loop system response are shown in Figures 3.5 to 3.8. From these plots, it is evident that the vehicle linear and angular speeds are stable with fairly small oscillation magnitudes at the beginning of motor failure. As for the propeller speed, the failed front right (FR) propeller speed decreased from 961.6 RPM to 750.6 RPM in 3 seconds with the central symmetry rear left propeller speed dropped correspondently. In order to provide enough thrust to hold the flight at the desired hovering condition, the front left (#2 propeller in Figure 2.7) and mid right (#3 in Figure 2.7) propeller speed increased from 961.6 RPM to 1075.3 RPM and 1073.7 RPM, respectively. Finally, the propeller speed of mid-left (#4 in Figure 2.7) , rear-right (#5 in Figure 2.7) and real-left (#6 in Figure 2.7) are near their normal hovering speeds. Based on the method stated in Eq. (3.1), the motor power consumption is showed in Figure 3.9. To hold the aircraft hovering with 47.6% available front right motor (#1) power, the required max motor power is 140% of normal hovering one with a well-tuned/designed adaptive MPC controller. 40 Figure 3.5: Linear speed of hovering failure - Case A Figure 3.6: Angular speed of hovering failure - Case A 41 Figure 3.7: Angular position of hovering failure - Case A Figure 3.8: Propeller speed of hovering failure - Case A 42 Figure 3.9: Motor power at hovering failure simulation - Case B 3.3.2.2 Hovering failure Case B In this case, we still assumed a 47.6% available front right (#1 in Figure 2.7) motor power but added detailed dynamics including propeller air drag in the model trimming process to the system model, which results in symmetric equilibrium conditions shown in Table 3.5, where the rear left (#6 in Figure 2.7) propeller speed is synchronized with the failed front right (#1 in Figure 2.7) one and the other four healthy propellers are operating at the same speed higher than their normal hovering speeds to compensate the thrust lost due to speed reduction of front right (#1 in Figure 2.7) and rear left (#6 in Figure 2.7) so that the hovering altitude can be maintained. Following a similar controller design process to what is stated in hovering Case A, the closed- loop vehicle rigid body speed responses are shown in Figures 3.10 and 3.11 with the related propeller speeds and powers showed in Figures 3.12 and 3.13. It is obvious that the aircraft is stabilized after front left (#1 in Figure 2.7) propeller is failed using the symmetric trimmed model since in this case the steady-state positions of pitch, roll, and yaw are not disturbed under closed-loop control and the control object is mainly focused on holding the aircraft at its desired altitude by regulating the vertical speed to zero. Note that the control designed 43 Figure 3.10: Linear speed of hovering failure - Case B Figure 3.11: Angular position of hovering failure - Case B 44 Figure 3.12: Propeller speed of hovering failure - Case B Figure 3.13: Propeller power of hovering failure - Case B 45 Table 3.5: Equilibrium condition of hovering failure - Case B Equilibrium condition 100% power 66% power 33% power Front right speed (rad/s) 100.7 87.7 69.6 Front left speed (rad/s) 100.7 106.7 113.2 Mid right speed (rad/s) 100.7 106.7 113.2 Mid left speed (rad/s) 100.7 106.7 113.2 Rear right speed (rad/s) 100.7 106.7 113.2 Rear left speed (rad/s) 100.7 87.7 69.6 based on the linear interpolated trimmed model (hovering Case A) could not hold the vehicle at the desired altitude under motor failure transitional operation. By looking into Figures 3.12 and 3.13, the failed front right propeller speed dropped from 961.6 RPM to 750.7 RPM, synchronized with the rear left propeller speed, under 52.4% power loss. The rest of propellers are accelerated to 1049.3 RPM with 129.9% power of normal hovering. 3.3.2.3 Hovering failure Case C In reality, when one motor failed instantly, it is impossible to activate the failure control without any delay showed in the previous ideal situation (Case B). Therefore, it is necessary to further study the control performance under certain non-equilibrium condition since after the motor failure, the system states could shift to non-equilibrium point. In this case study, a -2.9 deg/s initial pitching speed was used as the system initial condition, leading to the aircraft nose dipping. The adaptive MPC controller was tuned with large weighting placed on reference angular and vertical speed. The detailed controller parameters are shown in Table 3.8, where Q and R are calibration weights defined in Eq. (2.49) with all the unmentioned elements to be zeros. Also, U[1:6] denotes the real-time saturated control inputs which are the controlled rotational acceleration for each rotor, where each umax is calculated based on Eq. (3.10). Following the motor failure simulation structure used in Case 2 for validating the stability of 46 Table 3.6: MPC controller spec parameters for hovering - Case C Parameter Value Parameter Value Q[1,1] 1 Q[2,2] 1 Q[3,3] 10002 Q[4,4] 4002 Q[5,5] 4002 Q[6,6] 10002 R[1:6,1:6] 0.12 ∗ I6 Step size 1 ms Prediction Horizon 5 steps Control Horizon 1 steps min ∆u[1:6] 2 (t) −Kn[1:6] max ∆u[1:6] 2 (t)+max u −Kn[1:6] [1:6] a16 0 a22 0.06 a24 −0.45 a31 −0.025 a33 −0.3 a45 −1 s.t y[1,2] [−5, 5] s.t y[7,8] [−0.5, 0.5] re f s.t y[4,5,6] [−0.03.0.03] H(s) 1 0.2s+1 closed-loop MPC controlled system, the aircraft rigid-body speeds are shown in Figures 4.18 and 4.19. With -2.9 deg/s initial pitching rate shown at the right side plot in blue line, the vehicle accelerates toward the front-right direction with a slight opposite motion at the very beginning. The maximum magnitude of flight linear speeds are 2.3 km/h, 4.6km/h, and 0.2 km/h in x, y and z directions, respectively and finally they converge to zero in 40 seconds. The associated displacements of hovering failure recovery operation are summarized in Fig- ures 3.16 and 3.17. From the beginning of failure to stabilizing the aircraft in 40th seconds, the vehicle drifts 10.3 meters towards front, 8.6 meters towards right, and 1.1 meters down (altitude drop) with all angular speeds converged to zero. Propeller speed and power are shown in Figures 3.18 and 3.19. They are similar to these results for Case B, which is mainly due to individual motor control capability. In this case, the maximum motor power output is around 130.6% that is slightly larger than Case 2 value of 129.9%. Finally, motor driven torques are plotted in Figures 3.20 and 3.21, assuming a 4:1 gear ration from motor to propeller. Noting that the propeller spin direction is 47 Figure 3.14: Linear speed of hovering failure - Case C Figure 3.15: Angular speed of hovering failure - Case C 48 Figure 3.16: Linear displacement of hovering failure - Case C Figure 3.17: Angular position of hovering failure - Case C 49 Figure 3.18: Propeller speed of hovering failure - Case C Figure 3.19: Propeller power of hovering failure - Case C 50 opposite to adjacent propeller in the same row. To simplify the expression, driving torques of all propellers are normalized, indicating that positively increased torque leads to accelerating the associated propeller. Figure 3.20: First set of Motor torque of hovering failure - Case C Figure 3.21: Second set of Motor torque of hovering failure - Case C 51 By analyzing the nonzero initial condition such as pitching down, the controlled motor reaction can be explained intuitively. The most significant torque drop happened at rear left (#6 in Figure 2.7) and right (#5 in Figure 2.7) propellers to prevent the vehicle from pitching down with increased driving torque for the front left (#2 in Figure 2.7) propeller. However, due to the motor failure, the front right motor cannot provide enough thrust to balance the vehicle, leading to the utilization of extra power of mid two propellers. Last but not the least, failed front right motors slightly variated between 10 and 25 seconds to balance the flight yaw motion, indicating the significance of coordinated control under failed condition. 3.3.2.4 Behaviors comparison of hovering cases under disturbance Based on the obtained trimmed UAM vehicle models under different failure situations and the adaptive controllers designed, four disturbance tests are studied and discussed for comparing the flight behaviors under different disturbances and motor failure percentages stated in Table 3.7. Table 3.7: Condition definition of comparison cases Case Initial pitch rate (deg/s) Available power of front right motor (% of normal hovering) C1 -2.9 100 C2 -5.8 100 C3 -2.9 33 C4 -5.8 33 The response comparison between Cases C1 and C2 or between Cases C3 and C4 shows how the UAM behavior changes as a function of disturbance magnitude. By comparing Cases C1 and C3 or Cases C2 and C4, it can be seen that the impact to UAM aircraft behaviors from the available motor power for handling the identical disturbance. In Case C4, the controller has been recalibrated by increasing roll-rate weighting to handle this severe situation. Except for the failed front right motor, all the other available motor powers are not saturated to observe the required peak power levels. The comparison results are shown below. 52 The aircraft linear speeds and displacements are shown in Figures 3.22 to 3.29, and the ve- hicle angular speeds and positions are shown in Figures 3.30 to 3.37. Obviously, the oscillation magnitudes of rigid-body velocity increase as the initial pitch rate or motor failure percentage goes up. For investigating the control actuation behaviors, the propeller speed and power are shown in Figures 3.38 to 3.45. It can be seen from the plots that variations of propeller speed and power increase under larger disturbance or failure situations. Last but not the least, peak propeller power occurs at the front left (FL) motor for all four cases which are 104.5%, 109.2%, 144.4% and 148.9% from Cases C1 to C4, respectively. Figure 3.22: Linear speed comparison - Case C1 Figure 3.23: Linear speed comparison - Case C2 53 Figure 3.24: Linear speed comparison - Case C3 Figure 3.25: Linear speed comparison - Case C4 Figure 3.26: Linear displacement comparison - Case C1 54 Figure 3.27: Linear displacement comparison - Case C2 Figure 3.28: Linear displacement comparison - Case C3 55 Figure 3.29: Linear displacement comparison - Case C4 Figure 3.30: Angular speed comparison - Case C1 Figure 3.31: Angular speed comparison - Case C2 56 Figure 3.32: Angular speed comparison - Case C3 Figure 3.33: Angular speed comparison - Case C4 Figure 3.34: Angular position comparison - Case C1 57 Figure 3.35: Angular position comparison - Case C2 Figure 3.36: Angular position comparison - Case C3 Figure 3.37: Angular position comparison - Case C4 58 Figure 3.38: Propeller speed comparison - Case C1 Figure 3.39: Propeller speed comparison - Case C2 Figure 3.40: Propeller speed comparison - Case C3 59 Figure 3.41: Propeller speed comparison - Case C4 Figure 3.42: Motor power comparison - Case C1 Figure 3.43: Motor power comparison - Case C2 60 Figure 3.44: Motor power comparison - Case C3 Figure 3.45: Motor power comparison - Case C4 61 3.3.3 Simulation cases study based on analytical aircraft model In this section, case study based on nonlinear analytical aircraft model has been conducted for practical applications. Control design parameters are carefully tuned based on the rigid-body aircraft model responses with the simulation structure shown in Figure 2.9, and the detailed control design parameters are listed in Table 3.8, where these design parameters are defined in Eqs. (2.49) and (3.11). Table 3.8: Adaptive MPC controller spec parameters Parameter Value Parameter Value Q[1,1] 1 Q[2,2] 1 Q[3,3] 10002 Q[4,4] 200002 Q[5,5] 200002 Q[6,6] 10002 R[1:6,1:6] 0.12 ∗ I6 Step size 1 ms Prediction Horizon 20 steps Control Horizon 4 steps min ∆u[1:6] 2 (t) −Kn[1:6] max ∆u[1:6] 2 (t)+max u −Kn[1:6] [1:6] a16 −0.3 a22 0.1 a24 −1.5 a31 −0.1 a33 −1.5 a45 −1 s.t y[1,2] [−5, 5] s.t y[7,8] [−0.5, 0.5] re f s.t y[4,5,6] [−0.03.0.03] H(s) 1 0.2s+1 Firstly, we assumed that the available front right motor power failed from 100% to 66% to 33.3% (models 1 to 3 in Table 3.3) with stable hovering as the initial aircraft condition. To compare the system responses between LPV and nonlinear rigid-body models, the designed adaptive MPC controller was used to conduct closed-loop simulations based on both LPV and nonlinear rigid- body models, respectively, where the simulation results are denoted as LPV and rigid-body models in Figure 3.46. In this ideal failure situation, the rear left (#6 in Figure 2.7) propeller speed is synchronized with the failed front right (#1 in Figure 2.7) one and the rest four healthy propellers 62 are operated at the same speed higher than their normal hovering speeds to compensate the thrust lost due to speed reduction of front right (#1 in Figure 2.7) and rear left (#6 in Figure 2.7) propellers. Benefiting from these symmetric propeller responses, the oscillation dynamics only appears for the vertical speed; see Figure 3.46 for vertical speed and acceleration. Figure 3.46: Vertical speed and acceleration of ideal failure case In both cases, the aircraft vertical speed Vz and its acceleration VÛz are stablized in 5 seconds with a magnitude of speed oscillation less than 0.13 m/s. Small response differences can be observed at the beginning of simulation as the maximum vertical speed drop for the rigid-body model is slightly larger (-0.127 m/s), comparing with the LPV one (-0.119 m/s). For the rest of studies in this chapter, the simulation studies are based on the closed-loop system performance of nonlinear rigid-body aircraft model to make the simulation study practical. In practice, it is necessary to further study the control performance under certain non-zero condition since the aircraft could operate under non-equilibrium condition at the same time when motor failure happens. The scenario that UAM vehicle experience a sharp gust disturbance is considered in this study. The gust disturbance is assumed to generate a sudden pitch rate θÛx = −2.9 deg/s, leading to the aircraft nose dipping at the beginning (0 sec) of simulation. Three cases are studied in this section, assuming the aircraft is under normally hovering initially, with front-right motor (#1) operating under healthy condition (100%), partially failed condition 66% or 33% of normal hovering power, denoted as Cases HF1, HF2 and HF3, respectively. The vehicle rigid-body velocity are shown in Figures 3.47 and 3.48. 63 Figure 3.47: Linear speed of propeller #1 failure situations Figure 3.48: Angular speed of propeller #1 failure situations 64 With a θÛx = −2.9 deg/s initial pitching rate shown at the top plot of Figure 3.48, the vehicle moves forward at the very beginning with a maximum longitudinal speed Vy of 0.208, 0.213, and 0.22 m/s for Cases HF1, HF2, and HF3, respectively. Also, all three cases have very similar responses for Vx , θÛx , and θÛy under different available power levels for rotor #1. As for the vertical speed Vz , the oscillation magnitude increases as case number goes up from HF1 to HF3 due to rising motor failure severity. Last but not the least, it can be observed that the oscillation magnitude of θÛz also goes up as case number increases with no oscillation for Case HF1. This motion is mainly affected by the unbalanced air drag from propellers, which can be eliminated with symmetric propeller control actuation (see Case HF1). However, under motor failure situations (Cases HF2 and HF3), unsymmetric control could happen since the six motors are operated under different available power levels. To further understand and compare the vehicle behaviors, propeller speeds of Cases HF1 and HF3 are shown in Figure 3.49, where different colors present the propeller locations (e.g., FR: front-right) with the solid and dotted lines denote Cases HF1 and HF3, respectively. Recall the initial pitching down motion caused by simulated gust disturbance, controlled propeller motions can be explained intuitively. In general, the vehicle requires a pitch-up acceleration to recover from the initial pitching-down motion, which can be achieved by increasing these thrusts generated by front propellers and reducing rear ones. This is validated by the propeller speed responses of Case HF1 (solid line). However, due to the limited front-right motor power in Case HF3, the front-right (#1) motor cannot provide enough thrust to balance the vehicle, leading to extra effort of the middle two propellers (MR: #3 and ML: #4). With the failure gets severe, except the failed front-right propeller and its central symmetry one (RL: #6), all other propeller speeds increases correspondingly, with the peak speed for the front-left (FL: #2) one. Last, the failed front-right propeller speed is slightly varying between 3 and 6 sec to balance the vehicle yaw motion, indicating the significance of coordinated control even under failed condition. In summary, the adaptive MPC controller is able to stabilize the vehicle in 10 seconds for all three failure cases. Note that the failure studies for propellers #1, #2, #5, and #6 are the same since their locations 65 Figure 3.49: Speed of propellers in Case #1 and #3 are symmetric; see Figure 2.7. However, to validate the overall control performance with all possible single motor failure, the failure study of propeller #3 or #4 becomes necessary since its failure is different from that of propeller #1. As a result, propeller #4 is selected for failure study with the selected equilibrium conditions same as these for propeller #1 at 100% and 33% available power of normal hovering condition. The same adaptive MPC (see Table 3.8 for control design parameters) as these used for propeller #1 failure study is used for propeller #4 study, where 33% available power for the failed #4 motor with initial conditions θÛx = −2.9 deg/s, θÛx = −4.64 deg/s and θÛ = −2.9θÛx + −2.9θÛy deg/s are denoted by Cases HF4, HF5 and HF6, respectively. The rigid-body velocities are shown in Figures 3.50 and 3.51. In general, the aircraft response trends and characteristics are similar to these for Cases HF1, HF2 and HF3. It is worth to mention that for Cases HF4 and HF5 with initial pitch motion, the oscillations of roll and yaw angles are eliminated, compared with Cases HF2 and HF3, which is mainly benefited from the symmetric responses of left and right propellers. Propeller speeds are also plotted in Figure 3.52, where the highest speed occurs at the front left propeller that generates the most efficient thrust for eliminating 66 Figure 3.50: Linear speed of propeller #4 failure situations Figure 3.51: Angular speed of propeller #4 failure situations 67 Figure 3.52: Speed of propellers in Case HF6 the joint motion caused by initial pitch and roll disturbances. To validate the capability of proposed adaptive MPC, a normal MPC was designed and its design weighting is optimized based on the model with healthy motors (system at ρ = 1) with the same reference dynamic compensation and motor failure constraints. Note that u0 in the comparison simulations is the same and scheduled by ρ presenting the failed motor/propeller speed drop, but ∆u is different since ∆u from the normal MPC is designed based on the LTI system at ρ = 1 and ∆u from the adaptive MPC is based on the LPV model with ρ varying from 1 to 3. Assuming that both failure situations are the same as Case HF3, that is, the initial pitching speed is θÛx = −2.9 deg/s with 33% of normal hovering power available for the front-right (FR #1) motor, the rigid-body velocity are shown in Figures 3.53 and 3.54. It is obvious that flight speed oscillations are smaller in all x, y and z directions under adaptive MPC. The angular speeds are shown in Figure 3.54 with the related reference signal also plotted in dotted- and dashed-lines for adaptive and normal MPC, respectively. The angular speeds approach their reference signals under both control strategies, and the reference signals are also gradually 68 Figure 3.53: Linear speed comparison of AMPC and MPC Figure 3.54: Angular speed comparison of AMPC and MPC 69 converging to zero, with a stable flight. The most significant difference is shown at the pitch rate θÛy , where relatively large oscillation can be observed for the normal MPC. This is due to the fast initial propeller failure response, as the variation of control signal u = u0 + ∆u is dominated by u0 , leading to very similar responses in both cases. However, after 0.9 seconds, the variation of u0 becomes small and ∆u takes over the control gradually. In the normal MPC, the system model referred in the controller is not able to capture the aircraft dynamics variation caused by the failure but the adaptive MPC does, which leads to larger oscillation in θÛy . The propeller speed responses are also shown in Figure 3.55. It can be observed that the speed drop trajectories of the failed front-right (FR: #1) propeller and symmetric rear-left (RL: #6) one are the same for both cases at the beginning since it is a natural process because of the insufficient propeller driven torques. As for the other four propellers, the magnified region shows the start of propeller speeds deviations between the adaptive and normal MPC, since ∆u takes over the control as discussed above, leading to different recovering responses, especially the adaptive MPC strategy leads to a much smaller roll rate oscillation. Figure 3.55: Speed of propellers comparison of AMPC and MPC 70 In summary, the designed adaptive MPC is able to stabilize the aircraft under possible propeller failure cases with different failure percentages, propeller locations, and initial disturbances. Also, the adaptive MPC improves the closed-loop system performance since the LPV model used for MPC design varies based on the real-time flight condition, which reduces the modeling error between the control design model and the actual physical system. Lastly, the RMSE (root mean square error) of all studied cases are summarized in Table 4.6 with zero references. In summary, this section applies the adaptive MPC strategy to hovering control under propeller motor failure for a six-rotor UAM eVTOL (electric Vertical Take-Off and Landing) aircraft. An LPV model is developed based on the set of LTI models and used for designing adaptive MPC controller. The simulation results show that the designed adaptive MPC controller is able to stabilized the vehicle under multiple motor failure cases, and the overall RMSE of the vehicle rigid-body velocity oscillations can be significantly reduced using the adaptive MPC by 8.77% compared with normal MPC. 71 Table 3.9: RMSE of simulation results Case HF1 HF2 HF3(AMPC) HF3(MPC) Failed Propeller FR FR FR FR Available Power 100% 66% 33% 33% Direction of disturbance θÛx θÛx θÛx θÛx Magnitude of disturbance (deg/s) -2.9 -2.9 -2.9 -2.9 Vx 0 0.0010 0.0015 0.0108 Vy 0.0913 0.0933 0.0968 0.1003 Vz 0 0.0124 0.0283 0.0301 θÛx 0.5800 0.5870 0.5962 0.5814 θÛy 0 0.0044 0.0063 0.1308 θÛz 0 0.2031 0.3041 0.2791 Sum 0.6720 0.9013 1.0333 1.1325 Case HF4 HF5 HF6 Failed Propeller ML ML ML Available Power 33% 33% 33% Direction of disturbance θÛx θÛx θÛx & θÛy Magnitude of disturbance (deg/s) -2.9 -4.64 -2.9 & -2.9 Vx 0 0 0.0917 Vy 0.0937 0.3170 0.0952 Vz 0.0228 0.0233 0.0237 θÛx 0.5873 1.0650 0.5902 θÛy 0 0 0.5903 θÛz 0 0 0.1609 Sum 0.7038 1.4054 1.5520 72 CHAPTER 4 UAM TRANSITION/TILTING OPTIMIZATION AND CONTROL Many advanced flight control concepts were proposed aiming to improve the stability, flight qualities and performance of the tilt-rotor aircraft transition. In [42], the transition flight control of XV-15 was investigated using model predictive control (MPC) and time-varying linear control, respectively, and in [43] an adaptive model inversion control was used to improve the closed-loop performance and also to reduce the control development cycle. Contrary to the conventional tilt-rotor aircraft, the DEP (distributed electric propulsion) enabled hybrid eVTOL vehicles offer much more operational stability and agility, allowing for a high degree of fault tolerance and robustness. However, the aerodynamic complexity of transition flight from hovering to level flight and the flight control design for stable transition are still a challenging subject of research. The development of tractable tilt-rotor articulation profile for hybrid eVTOL vehicles is also an open question. Given the vehicle specifications and operational constraints; including the flying qualities and passenger ride quality, the goal of articulation profile design is to develop a tilting trajectory that is best suited for the type of aircraft at its given flight conditions. In [44], the tilting trajectory was optimized for a quadrotor aircraft using the ant colony optimization algorithm, and in [45] the takeoff trajectory of a tilt-rotor aircraft was optimized using the direct allocation method. In [46, 47] the dynamic modeling and control of a hybrid eVTOL vehicle configured with four synchronized tilt-rotors were considered. There are other eVTOL vehicles that utilize a DEP system with various number of lift-rotors and tilt-rotors. For example, Joby S4 and Uber eCRM-001 are both configured with six rotors, but Joby S4 has only two tilt-rotors for synchronous articulation, whereas in eCRM-001 all six rotors can be articulated. Another DEP enabled, high technology readiness level (TRL) hybrid eVTOL vehicle is Solvay VA-1X, it features four tilt-rotors and four lift-rotors, separately. It is worth to mention that not all rotors needs to be tiltable due to the following two reasons. Firstly, the vertical thrust required for maintaining hovering is always much larger than the longitudinal thrust required for maintaining the desired cruising speed, which means 73 that only part of rotors needs to be tilted. Secondly, in order to minimize the weight added by the tilting mechanism leads to using minimal number of tilting rotors. For example, the Uber eCRM-001 and Solvay VA-1X have two and four tiltable rotors for six and eight rotors in total, respectively. Thus, regarding the proposed UAM configuration, the middle two rotors are tiltable, and the two front and two rear ones are fixed vertically. Also, optimization of the trajectory also becomes an open question that is aimed to generate the optimal reference tilt angle and aircraft longitudinal speed for minimal energy usage with peak motor power and passenger comfort level constraints. For example, at the beginning of tilting, the achievable longitudinal acceleration is up-limited by the tilt-rotor angle since when the rotor angle is close to the vertical position, the available longitudinal thrust is very small. As a result, maximum available longitudinal acceleration is a function of tilting angle and it needs to be considered during the optimization process. In [44], the tilting trajectory was optimized for a quadrotor aircraft using the ant colony optimization algorithm, and the take-off trajectory of a tilt-rotor aircraft was optimized in [45] using the direct allocation method. In this chapter, we study three tilt-rotor articulation profiles for transition flight, subject to hardware constraints and flying/ride quality specifications. In particular, by leveraging the dis- tributed rotors for reconfigurable thrust vector, more aggressive and optimal tilting profiles can be developed, where the optimal profile consumes the minimum total propulsion energy. A control- oriented LPV model [30, 27, 31] is developed to characterize the aircraft flight transition process along a prescribed articulation profile. Through linearization of nonlinear aircraft model with respect to a set of equilibrium flight conditions, a set of linear time-invariant models can be at- tained [22, 26]. However, such a modeling approach leads to a large modeling error, because of the presence of nonzero longitudinal vehicle acceleration during the transition flight. In order to reduce the LPV modeling error, a trajectory linearization method is adopted in this study, hence a set of discretized and linearized models can be obtained through trajectory linearization along the tilting profile, furthermore, these discrete-time linearized models are then used to construct the LPV aircraft transition flight model. The adaptive MPC method [28, 29, 32, 33] is used for 74 designing the transition flight controllers based on the developed LPV model [34, 35, 36]. Note that the MPC utilizes a specific set of penalty weightings to shape the target output performance. In this chapter, in order to improve the closed-loop system performance and reduce the real-time computational cost, two critical modifications are introduced to the proposed adaptive MPC law, namely, integration of feedforward control via dynamic reference compensation (DRC) and uti- lization of event-driven/time-varying weightings and hard constraints. Finally, the performance of the proposed adaptive MPC-DRC strategy is evaluated in Matlab™-Simulink™ environment using the nonlinear rigid-body eVTOL model. Additional simulations are conducted to demonstrate the disturbance rejection performance of the proposed adaptive MPC-DRC transition flight controllers. 4.1 Steady-state dynamics analyzation and preliminary control study Different from hovering cases, the propeller tilting motions should be considered in the tilt transition process. The relationship between tilt angle α and desired level flight speed v becomes a trajectory planning problem as v is a function of α. As a result, the tilt control can be formulated as a tracking problem for given target tilt angle and level flight speed trajectories. Based on the predefined tilting transition dynamics, multiple trimmed models can be obtained. In this section, the steady state dynamics of aircraft is investigated at the trimmed points listed in Table 4.1. The eigenvalues of each model are showed in Figure 4.1, where blue, green, red, and black colors represent the poles associated with tilt angles of 89, 60, 30, and 0 degree, respectively. Pole locations of five dominant aircraft modes are plotted and denoted by Modes 1 to 5 using "×," "o," "+," "∗," and "4," respectively, and associated mode physical meaning and values are summarized in Table 4.2. Note that modes 2, 4 and 5 are complex and their natural frequencies ωn and damping ratio ζn (n = 2, 4, 5) can be calculated based on Eq. (4.1) Re(λn ) q ωn = Re(λn )2 + Im(λn )2, ζn = (4.1) ωn where Re(λn ) and Im(λn ) denote the real and imaginary parts of the nth eigenvalue. By looking into natural frequencies and damping ratios, the aircraft oscillation frequencies and decay rates can be obtained. 75 Table 4.1: Equilibrium condition of tilting Equilibrium conditions Model 1 Model 2 Model 3 Model 4 Tilt angle (deg) 89 60 30 0 Forward speed (m/s) 1 22.7 45.3 68 Front right speed (rad/s) 102.3 120.0 96.0 0 Front left speed (rad/s) 102.3 120.0 96.0 0 Mid right speed (rad/s) 102.2 33.0 37.1 40.1 Mid left speed (rad/s) 102.2 33.0 37.1 40.1 Rear right speed (rad/s) 102.1 96.3 26.7 0 Rear left speed (rad/s) 102.1 96.3 26.7 0 Elevator deflection (rad) 0.349 0.349 0.138 0 89 Hovering M4 Tilt angle (deg) 60 M4 M5 30 M3 M2 M1 0 Level flight Figure 4.1: Root loci during transition 76 The rolling and short period modes (#1 and #4) remain stable during the entire tilt transition process with increasing natural frequencies. The Dutch roll mode (#2) varies from unstable to stable with increasing damping ratio, which is benefited from the increased fixed wing aerodynamic lift. As for the phugoid mode (#5), it finally becomes stable when the aircraft achieves level flight (cruising) condition. Last but not the least, note that the spiral mode (#3) is always slow and unstable, which is common characteristics of aircraft. It means that a small roll angle change caused by external disturbance could slowly diverge the vehicle rolling angle, yaw angle, and altitude, leading to a spiral flight trajectory towards the ground. The aircraft dynamics under multiple transition conditions are modeled as an LPV system based on four linearized models defined in Table 4.1. And the obtained LTI system matrices are formulated with system inputs, states and outputs defined in Table 4.3. Following a similar controller design process stated in hovering failure study (last chapter), the preliminary case study is introduced aiming to control the aircraft from hovering to level flight in 45 seconds with a cruising speed of 68 m/s. Firstly, it can be seen from Figures 4.2 and 4.3 that the level flight speed accelerated from 0 to 68 m/s from beginning to end. However, there is a certain overshoot of the aircraft level flight speed when the designed controller was evaluated using both LPV and nonlinear models around the 45-th second. Besides, the acceleration trends are similar for both cases with their certain magnitude difference, which could be caused by the LPV modeling error using the interpolated linear models obtained by linearizing the nonlinear one at equilibrium condition. The vertical speed and displacement are shown in Figures 4.4 and 4.5. Figure 4.2: Level flight speed 77 Table 4.2: Information of aircraft modes Mode ID Rigid body mode 89deg 1m/s 60deg 22.7m/s 1 Rolling -0.3008 -2.5855 Dutch roll 0.0892±0.2293i -0.0587±0.7145i 2 ωn 0.246 0.7169 ζ -0.3625 0.0819 3 Spiral 0.0049 0.0362 Short period -0.1117±0.1043i -1.1696±1.6673i 4 ωn 0.1528 2.0366 ζ 0.7309 0.5743 5 Phugoid 0.0605±0.087i 0.0088±0.18i ωn 0.106 0.1802 ζ -0.5709 -0.0488 Mode ID Rigid body mode 30deg 45.3m/s 0deg 68m/s 1 Rolling -5.0809 -7.6005 Dutch roll -0.1396±1.248i -0.2071 ±1.8156i 2 ωn 1.2558 1.8274 ζ 0.1112 0.1133 3 Spiral 0.0240 0.0147 Short period -2.3239 ± 3.3413i -3.4814±5.0149i 4 ωn 4.07 6.1049 ζ 0.571 0.5703 5 Phugoid 0.0011 ± 0.1803i -0.003±0.1804i ωn 0.1804 0.1804 ζ -0.0061 0.0166 78 Table 4.3: System input and output definition System Input Order System State Order Front R/L rotor acceleration u1 /u2 Rigid-body velocity x1 to x6 Mid R/Lrotor acceleration u3 /u4 Euler angle x7 to x9 Rear R/L rotor acceleration u5 /u6 Body inertial position x10 to x12 Mid R/L rotor tilt driven torque u7 /u8 Rotor speed x13 to x18 Elevator Deflection u9 Mid right/left rotor tilt angle x19 to x20 Mid R/L rotor tilt speed x21 to x22 System Output Order System Output Order Lateral speed (Vx ) y1 Roll angle y7 Longitudinal speed (Vy ) y2 Pitch angle y8 Vertical speed (Vz ) y3 Yaw angle y9 Pitch rate (θÛx ) y4 Lateral displacement y10 Roll rate (θÛy ) y5 Longitudinal displacement y11 Opposite value of yaw rate (θÛz ) y6 Vertical displacement y12 Mid R/L rotor tilt angle y13/14 Mid R/L rotor tilt speed y15/16 Figure 4.3: Level flight acceleration 79 Figure 4.4: Vertical speed Figure 4.5: Vertical displacement Although variations of vertical speed from both models are small at the beginning, difference of vertical displacement increases as the time goes up due to the accumulated error. After the 45-th second, large speed drop is observed due to the control mode switch from normal tilt transition (MPC) to the level flight (linear control), and finally, for both cases, vertical speed and displacement errors are reduced down to zero. Figures 4.2, 4.6, and 4.7 show how vehicle pitch angle and rate vary with respect to the level flight acceleration. In the first 45 seconds, the aircraft pitches down to bring the net total thrust vector of both front and rear four propellers forward to provide additional thrust for aircraft acceleration. Then, to reduce the speed overshoot shown in Figure 4.2, the flight pitches up at the 45th second to decelerate the vehicle to reach the equilibrium speed. The lateral speed, roll and yaw rates are held at zero stably during the tilt transition since there is no excitation to change them due to synchronized left and right propeller speeds. As a result, propeller speed plots are simplified as front (blue), mid (black), and rear (red); see Figure 4.8. All propeller speeds track the trimmed solution accurately and finally converge to the wing-lift level flight equilibrium solution at about 70th second for both cases. However, the response differences from LPV and 80 nonlinear models gradually increase because of the accumulative modeling errors, especially after the 45-th second. Figure 4.6: Pitch rate Figure 4.7: Pitch angle Figure 4.8: Propeller speed In summary, the steady-state aircraft transition can be achieved based on the steady-state LPV transition model. However, the control performance is limited mainly due to two reasons. Firstly, 81 the rotor tilt angle and related speed are predefined which may not be achievable during transition. Secondly, because of the nonzero level flight acceleration during the transition, there exists certain LPV modeling error caused by trimming the nonlinear aircraft model at their static equilibrium conditions. 4.2 Reference tilting trajectory generation and optimization To optimize tilting performance, a tracking control problem was formulated to follow a reference tilting trajectory including both tilting angle and level flight speed. As a result, it is necessary to find the optimal tilting trajectory based on the control design LPV model, along with the level flight speed. 4.2.1 Simplified aircraft force analysis model Figure 4.9: Geometry of UAM aircraft with dynamic force definition Firstly, as shown in the left drawing of Figure 4.9, the flight pitch and roll angles, vertical and lateral speed values are assumed to be zero during the tilting process, and the total drag Fd and lift Fl generated by rotors and wing are with respect to CG (center of gravity) of aircraft. To be more specific, each of front and rear rotors generates the same synchronized thrust Tl and each of mid two tilted rotors produces thrust Tt . Action points of Tl and Tt can be moved to CG since the rotor locations are symmetric. Noting that the pivot point of mid tiltable rotor is 0.5 m above the CG, 82 meaning that thrust vector Tt will slightly deviate from the aircraft CG once the rotor tilts for the vertical position, and this small deviation is omitted to simplify the energy optimization discussed later. mg and ma are the gravity and required force for accelerating the aircraft at a. And the middle rotor tilt angle ρ is defined as ρ(t) ∈ [0, 90] deg, and it is also assumed that the tilting motion of mid two rotors are synchronized. Based on the above assumptions, the force balancing equations in both vertical and longitudinal direction are stated as follows:  Fl + 4Tl + 2sin(ρ(t))Tt = mg    (4.2)  2cos(ρ(t))Tt − Fd = ma(t)     where the calculations of each force element are stated in Eq. (4.3)  F = 1 C ρ Av (t)2, F = 1 C ρ Av (t)2  a l a l 2 l d 2 d  l  (4.3)  2 5 2  Tt = Ct ρa nt (t) D , Tl = Ct ρa nl (t) D   5  where nl and nt are the lift and tilt rotor speed, Cl and Cd are the lift and drag coefficients, and ρa , A and D present the air density, wing area and rotor diameter, respectively. During the aircraft transition from hovering mode to level flight one, the aircraft gravity mg will be balanced mainly by the vertical rotors thrust 4Tl at the beginning, and gradually shifting to the aerodynamic lift force Fl as level flight speed increases. At the end of transition, the aerodynamic lift force Fl will be equal to the aircraft weight mg = Fl when the aircraft cruises at target speed vc . The transition process from hovering mode to level flight one is defined as the middle rotors tilted from "tilt-up (ρ(t) = 90 deg)" position to "tilt-forward (ρ(t) = 0 deg)", where the level flight speed increased gradually from 0 m/s to the target cruising speed of 68 m/s. 4.2.2 Tilt transition optimization The goal of transition flight optimization study is to develop the most efficient tilt-rotor articulation profile that minimizes the total energy cost to accomplish the transition flight. To achieve this, it is critical to carefully modulate the transition duration, subjected to the constraints on allowable peak rotor thrust, aircraft flying qualities, and passenger ride quality. In this study, three cases of transition flight strategy are considered with respect to a scheduled transition timeline. 83 As a baseline study, Case 1 presents an ordinary tilt-rotor articulation profile, featuring an initial ramping-up period, followed by cruising at constant tilt-rate, and finally, winding down to complete the transition flight. More specifically, from 0 to 1s, the tilt-rotor articulation is to linearly accelerate and reach the targeted tilt-rate of 2 deg/s, and from 1s to 45s, the tilt-rate is to coast at 2 deg/s, finally, from 45s to 46s, the tilt-rate is to reduce linearly to zero. A simple calculation over this profile confirms that the total angular distance traveled is 90-deg. At the beginning of transition flight t0 , the tilt-rotor rotational speeds are held at a predetermined maximum speed Np , and as tilt-rotors begin to articulate the aircraft starts to accelerate and pick up forward velocity, hence the longitudinal acceleration reaches the maximum allowable amax at t1 . From t1 to t2 , the tilt-rotor thrusts are regulated so as to maintain the longitudinal acceleration at amax , while the aircraft continues to gain forward velocity. As a result, the aerodynamic lift force induced by fixed-wing emerges, correspondingly, the four lift-rotor thrust levels can be tailored to maintain altitude-hold. From t2 to t3 , the longitudinal acceleration is to reduce from amax to zero by adjusting the tilt-rotor thrust level, and the aircraft is to reach or continue the desired level flight velocity Vc . The targeted coasting tilt-rate of 2 deg/s in Case 1 is chosen based on the early study conducted for a twin tilt-rotor aircraft XV-15. However, with the proposed distributed electric propulsion platform, a higher tilt-rate can be commanded without compromising the vehicle stability during transition flight. The primary goal of rapid tilt-rotor maneuver is to develop vehicle forward velocity as soon as possible, hence the aerodynamic lift can be generated much sooner so as to save the total energy. Therefore, in Case 2 we propose an aggressive tilt-rotor articulation profile aiming for the aircraft to quickly reach the maximum longitudinal acceleration amax at t1 by commanding maximum tilt-rate ρÛ max and maximum thrust level Tt1 . At t1 the corresponding vehicle forward velocity Vy (t1 ) = Vy1 . Specifically, from t1 to t3 , the vehicle reference acceleration and velocity profiles are proposed as follows,  amax , ∀t ∈ [t1, t2 ]    a(t) = (4.4)  t3 − t  amax   , ∀t ∈ [t2, t3 ]  t3 − t2 84 ∫ t 3 Vy (t) = a(t)dt + Vy1, s.t.Vy (t1 ) = Vy1, Vy (t3 ) = Vc (4.5) t1 where amax , Vc , Vy1 , and t1 are known, t2 and t3 can be calculated respectively based on the predefined level flight convergence duration ts , where ts = t3 − t2 . In this case, ts is chosen to be 5 s, t1 and t2 are calculated as 3.215 and 36.58 s, respectively, hence t3 = 41.58 s. Apparently, the proposed aggressive tilt-rotor articulation profile can indeed expedite the transition flight when compared to the baseline profile. In Case 3, we consider the same scheduled timeline as in Case 2, but propose an optimal articulation profile that consumes minimum energy. In particular, from t1 to t3 , total energy E and instant power P(t) of all six rotors can be given by ∫ t 3 E= P(t)dt, P(t) = 4Cp ρa nl (t)3 D5 + 2Cp ρa nt (t)3 D5 (4.6) t1 where Cp denotes the power coefficient, nl (t) and nt (t) are lift- and tilt-rotor rotational speeds derived from Eqs. (4.2) and (4.3) as 1 C ρ AV (t)2 3 #1 +  "   ma(t) 2 d a y nt (t) =   2Ct ρa D5 cos(ρ(t))    (4.7)   #1 3 D5 − 1 C ρ AV (t)2 3 2 sin(ρ(t))C ρ  " mg − n (t)   t a t 2 l a y  nl (t) =   5  4C ρ   t a D where a(t) and Vy (t) are known. By substituting nl (t) and nt (t) into Eq. (4.7), the problem of energy minimization is now converted into finding ρ(t) that minimizes the total energy cost E, i.e. ∫ t ∫ t h 3 3 4Cp ρa nl (t)3 D5 + 2Cp ρa nt (t)3 D5 dt . i min E = min P(t)dt = min (4.8) ρ(t) t1 ρ(t) t1 Note that it is difficult to directly solve the optimization problem as presented in the form of Eq. (4.8). Thus, we propose that the integral in Eq. (4.8) be approximated using trapezoidal quadrature rule by dividing integration horizon [t1, t3 ] into Nd sections with a step size δ = (t3 − t1 )/Nd , i.e. NÕ d −1 δ ∫ t 3 P(t)dt ≈ (P(k) + P(k + 1)) , (4.9) t1 2 k=0 85 and the approximation accuracy can be guaranteed if Nd is chosen to be sufficiently large. We assume that tilt angle ρ ∈ [ρi , ρi+1 ] can be linearly interpolated between the two sample points, ρi and ρi+1 . Note that Pk = P(k) ≥ 0 is a constrained nonlinear function of ρ k = ρ(k) defined by Eqs. (4.8) and (4.9). Hence, the optimal tilt angle ρ k , for k = 0, 1, . . . , Nd − 1, can be obtained by solving the following optimization problem min Pk (ρ k ) for k = 0, 1, . . . , Nd − 1 . (4.10) ρk By substituting relevant aircraft parameters, given in Table 4.4, into Eq. (4.10), the optimal tilt-rotor articulation profile from t1 to t3 can be attained by solving dPk /dρ k = 0, for k = 0, 1, . . . , Nd − 1. Figure 4.10 shows the reference aircraft forward velocity Vy (t), the tilt-rotor articulation profile ρ(t), and the tilt-rotor rotational speed nt (t) for the three cases described above. They all have successfully transitioned to the desired level flight velocity within the specified duration, however, the total energy cost for Case 1 is 3.05 kWh, for Case 2 is 2.64 kWh, and for Case 3 is 2.56 kWh. This further proves the advantage of the proposed minimum energy reference profile. Table 4.4: Aircraft parameters Parameter Value Parameter Value Parameter Value Cl 0.3141 Cd 0.0166 Cp 8.4 × 10−5 Ct 5.6 × 10−4 ρa, kg/m3 1.225 A, m2 25 D, m 3.51 amax, m/s2 1.85 ρÛ max, deg/s 9 t s, s 5 Np, r pm 1146 g, m/s2 9.81 4.3 Control design This section presents an LPV model that captures the tilt transition flight dynamics as the two tilt-rotor thrust vectors are articulated from vertical to horizontal. 86 Figure 4.10: Comparison of three reference trajectories during transition flight 4.3.1 Linearization along a nominal trajectory Consider a nonlinear aircraft model in the state-space form described by Û = f (x(t), u(t)) x(t) (4.11) where f is a differentiable function of state x(t) ∈ Rn x and control input u(t) ∈ Rnu . We assume that x(t) = x0 (t) + ∆x(t) and u(t) = u0 (t) + ∆u(t), where x0 (t) and u0 (t) denote the nominal state and control, respectively, and ∆x(t) and ∆u(t) a small deviation from the nominal conditions. In addition, x0 (t) and u0 (t) satisfy xÛ0 (t) = f (x0 (t), u0 (t)) , (4.12) therefore, Eq. (4.11) can be rewritten as xÛ0 (t) + ∆ x(t) Û = f (x0 (t) + ∆x(t), u0 (t) + ∆u(t)) . (4.13) By applying the Taylor expansion to the above, we obtain xÛ0 (t) + ∆ x(t) Û = f (x0 (t), u0 (t)) + A∆x(t) + B∆u(t) + Ψ f (∆x(t), ∆ x(t),Û ∆u(t)) (4.14) 87 where A and B represent the sensitivity matrices at the given nominal conditions, and Ψ f denotes the collection of higher order terms that are negligible. Substituting Eq. (4.12) into Eq. (4.14) yields Û = A∆x(t) + B∆u(t) . ∆ x(t) (4.15) Since ∆ x(t) Û = x(t)Û − xÛ0 (t), Eq. (4.15) can be rewritten as below, assuming that the system output y(t) is a linear function of system states, where C is a constant matrix. Û = xÛ0 (t) + A∆x(t) + B∆u(t)   x(t)   (4.16)   y(t) = C x(t)    which can be converted into the discrete-time form with a sampling time T as  x(k + 1) − x0 (k + 1) = xÛ0 (k)T + A∆x(k) + B∆u(k)    (4.17)   y(k) = C x(k) , k = 0, 1, ...    where (A, B) are the discrete-time system matrices derived from (A, B), and x0 (k) and xÛ0 (k) represent the nominal state and its derivative at k. 4.3.2 Affine LPV model along the nominal trajectory An LPV model can be developed by linearly interpolating the state-space models at two adjacent operational conditions along the nominal trajectory over a varying parameter (or scheduling param- eter) [38]. We assume that the tilt-rotor angular position ρ(t) is measurable in real-time. Consider a family of M linearized system models obtained along the given nominal tilting trajectory, which covers from the start to the end of tilting sequence with ρ as the scheduling parameter. This family of linear models can be defined as i i i i  x(k + 1) = x0 (k + 1) + xÛ0 (k)T + A ∆x(k) + B ∆u(k)    (4.18)  y(k) = C x(k)     88 where i = 1, 2, ..., M denotes the i-th linearized model at ρi . Therefore, the discrete-time affine LPV model can be formulated as follows,  x(k + 1) = x0 (ρ(k + 1)) + xÛ0 (ρ(k))T + A(ρ(k))∆x(k) + B(ρ(k))∆u(k)    (4.19)  y(k) = C x(k)     where ∆x(k) = x(k) − x0 (ρ(k)) and ∆u(k) = u(k) − u0 (ρ(k)), and the system matrices A(ρ(k)) and B(ρ(k)) are given by ρ(k) − ρi A(ρ(k)) = Ai + (Ai+1 − Ai ) , ρ(k) ∈ [ρi, ρi+1 ], i = 1, 2, ..., M − 1 ρi+1 − ρi (4.20) i i+1 i ρ(k) − ρi B(ρ(k)) = B + (B − B ) , ρ(k) ∈ [ρi, ρi+1 ], i = 1, 2, ..., M − 1 ρ −ρ i+1 i It should be noted that ρ(k + 1) in Eq. (4.19) cannot be measured at time step k, however, if the sampling rate is chosen to be sufficiently high, then the difference between the two time steps is negligible, hence we can assume that ρ(k + 1) ≈ ρ(k) when T is very small. Therefore, the resulting LPV model approximates the nonlinear aircraft dynamics along the desired operational trajectory as a function of ρ(k), in which the sampling period T is chosen to be 1 msec. The following list summarizes steps of LPV modeling process. 1) Define the aircraft transition trajectory in terms of ρ(t), Vy (t), and a(t) over the transition flight period [t0, t3 ], based on the approaches specified in Section 4.2. 2) Define the number (M) of LTI models to be linearized along the tilting trajectory and select the reference points ρi for i = 1, 2, ..., M along the tilting trajectory. For example, when M = 4, the reference tilting angles can be selected as ρ1 = 90 deg, ρ2 = 60 deg, ρ3 = 30 deg and ρ4 = 0 deg. In this study, we choose M = 20. 3) Calculate the nominal control input ui0 , e.g. rotor speeds, corresponding to ρi for i = 1, 2, ..., M, based on the nonlinear aircraft model using Newton-Raphson search method. 4) Linearize the nonlinear aircraft model at each tilting angle ρi and [x0i , xÛ0i , ui0 ] to obtain the system matrices Ai and Bi for i = 1, 2, ..., M. 89 5) Assemble the affine LPV model in the form of Eq. (4.19) with the matrix interpolations defined in Eq. (4.20). 4.3.3 Adaptive MPC of aircraft transition From hovering to level flight mode, the varying aircraft dynamics is mainly due to the gradually increased level flight speed and rotor tilt angle change. As shown in Figure 4.11, for tracking the desired level flight speed at start of tilting, it was found through calibrating the adaptive MPC control law that even though large penalization was used for output weighting Q related to level flight speed, the level flight speed performance does not improve much, leading to large tracking error. The main reason is that both rotors are at near tilt-up position, and the control of level flight speed can only be achieved by pitching down aircraft to change rotor thrust vector forward or backward, resulting in slow response due to relatively slow aircraft rigid body dynamics. On the other hand, when the aircraft is close to the level flight mode as both tiltable rotors are becoming push rotors, the level flight speed can be directly controlled by manipulating the tilting rotor thrust for aircraft acceleration or deceleration. Figure 4.11: Variation of aircraft dynamics To track the level flight speed near the hovering mode, a dynamic reference compensation method was proposed. Since the aircraft pitch rate is directly affected by the rotor thrusts, the level flight re f speed can be compensated by using a real-time feedback pitch rate reference ∆x4 = Gl (ρ)∆x2 , where ∆x2 is the level flight speed error, and Gl (ρ) is the sensitivity gain of dynamic compensation. As a result, the MPC design has a small weighting Q[2,2] (ρ) and a large sensitivity gain Gl (ρ) 90 for level flight speed control at the beginning of tilting since the aircraft is close to the hovering mode; and error weighting Q[2,2] (ρ) will be gradually increased as the tilting transition continuous, while the sensitivity gain Gl (ρ) reduced correspondingly. Using the same principle, the aircraft pitch angle control can be achieved by adjusting the reference angular speed based on the feedback angular position with a calibratable negative gain. Thus, the reference pitch rate can be defined re f by ∆x4 = −K p ∆x8 + Gl (ρ)∆x2 . In summary, the desired level flight tracking performance is achieved by adjusting pitch rate reference as a function of level flight speed error and pitch angle along with increases penalization of level flight speed error. Following the same principle, for the purpose of maintaining the aircraft altitude, the reference vertical speed ∆x3 = −Kv ∆x12 is compensated by the deviation of altitude ∆x12 , where Kv is the selected sensitivity gain. Lastly, since the mid-rotor tilt angle ρ(t) can be tracked directly by regulating its acceleration ρ(t) Ü (or u8,9 = ρ(t)),Ü directly, that is also the nominal control u0 along the tilting trajectory. Then, in the MPC scheme, ∆u8,9 are constrained by −ε ≤ ∆u8,9 ≤ +ε, where ε is chosen to allow the tilting angle to be mainly controlled by the nominal control signal. In this study, the penalized tracking outputs are the longitudinal speed ∆x2 , vertical speed ∆x3 re f re f and pitch rate x4 , where the reference vertical speed ∆x3 and pitch rate ∆x4 will be adjusted based on the related system states (see the discussion above and Table 4.5) and the reference vertical re f speed ∆x2 will be held at zero. 4.3.4 Adaptive MPC formulation The proposed adaptive MPC integrated with a dynamic reference compensation (DRC) architecture is shown in Figure 4.12, where ∆x and ∆u denote the deviations of the controlled states and inputs as described in Eq. (4.19). Based on the real-time state reference signal ∆x re f and the system state feedback signal ∆x, the adaptive MPC law generates a set of optimized control signals ∆u(k + m), m = 0, 1, ..., N − 1, from which the control signals within the control horizon, i.e. [∆u(k) · · · ∆u(k + Nc − 1)] and Nc < N, are selected as the optimal state-feedback control ∆u and to be combined with the nominal control input u0 (ρ) to form a complete control input, 91 u = ∆u + u0 (ρ). In this study, all the system states x are assumed to be measurable, and the required control feedback signal ∆x is defined as ∆x = x − x0 (ρ), where x0 (ρ) denotes the nominal state. Figure 4.12: Adaptive MPC architecture for transition control Letting e(k) = ∆xre f − ∆x(k) to denote the tracking error at current time step k, the adaptive MPC design [34, 35, 36] is to find the constrained optimal control law ∆u(k) over a given horizon N that minimizes the constrained quadratic performance defined by N N−1 1 Õ T Õ min { e (k + m)Qe(k + m) + ∆uT (k + m)R∆u(k + m)} ∆u(k),...,∆u(k+N−1) 2 (4.21) m=0 m=0 subject to G∆u(k + m) ≤ h , m = 0, ..., N − 1 where Q ≥ 0 and R > 0 are the weighting matrices used to penalize the tracking error e(k) and control effort ∆u(k), respectively, and G a constraint matrix used to ensure that the control effort ∆u(k) stay within the prescribed bound h. In addition, e(k + m) and ∆u(k + m) represent the predicted error and control input at m-th time step. For a finite prediction horizon of N steps, the cost function and the constraint equation in Eq. (2.49) can be rewritten in a compact form as follows, 1 min [êT (k)Q̂ ê(k) + ∆ûT (k) R̂∆û(k)] û(k) 2 (4.22) subject to Ĝ∆û(k) ≤ ĥ 92 where        e(k)   ∆u(k)   h              e(k + 1)   ∆u(k + 1)   h ê(k) =   ∈ R(N+1)n x , ∆û(k) = , ĥ =  .  ∈ R Nnu       ..  ..  ..    .     .           e(k + N) ∆u(k + N − 1)        h       0 . . . 0  G 0 . . . 0      R        0 R . . . 0  0 G . . . 0 R̂ =  . ..  , Ĝ =  ∈ R Nnu ×Nnu ,     .. .. . . . .  .. .. . . . ...   . . .       0 0 . . . R  0 0 . . . G         0 . . . 0    Q    0 Q . . . 0  Q̂ =  . (N+1)n x ×Nn x . ..  ∈ R   .. .. . .  . . .   0 0 . . . Q     (4.23) Furthermore, the predicted tracking error (over a given horizon) ê(k) can be described by ê(k) = Âe(k) + B̂∆û(k) , (4.24) where the predicted tracking error ê(k) can be determined based on the current tracking error e(k) and control vector ∆û(k). Therefore, the optimization problem can be reformulated by utilizing the current information as 1 T 1 min ∆û (k)( R̂ + B̂T Q̂ B̂)∆û(k) + eT (k) ÂT Q̂ B̂∆û(k) + eT (k) ÂT Q̂ Âe(k) ∆û(k) 2 2 (4.25) subject to Ĝ∆û(k) ≤ ĥ(k) 4.4 Simulation Results and Discussions Based on the simulation structure showed in Figure 2.9, the MPC law is carefully calibrated to achieve the desired aircraft transition performance, and the specific design parameters are listed 93 90−ρ in in Table 4.5, where G1 (ρ) = 90 is linearly increased from 0 to 1 along the transition process 90−ρ from ρ = 90 deg to ρ = 0 deg, and G2 (ρ) = 1 − 90 is linearly decreased from 1 to 0, respectively. The maximum level flight acceleration is constrained at 1.85m/s2 for flight comfortableness, and the mid-rotor tilting speed is restricted to be below −9deg/s for this simulation study. Five cases are studied, including normal transition (Case 1), optimal transition (Case 2) with its baseline comparison (Case 3), and transition under gust wind disturbance (Cases 4 and 5), where all five simulation studies are performed using the nonlinear aircraft model (see [48] for details). Firstly, Table 4.5: Adaptive MPC design specification parameters Parameter Value Parameter Value Q[2,2] G1 (ρ) ∗ 4002 Q[3,3] 2002 Q[4,4] 40002 R[1:6,1:6] 0.012 ∗ I6 R[9,9] 0.012 Step size 1 ms Prediction Horizon 4 steps Control Horizon 2 steps min ∆u[1:6] 2 −Kn[10:15] (t) max ∆u[1:6] 2 −Kn[10:15] (t)+max u[10:15] min ∆u[7,8] −ε max ∆u[7,8] ε min ∆u[9] −0.2 max ∆u[9] 0.2 re f re f ∆x3 −∆x12 ∆x4 −∆x8 + 0.01G2 (ρ)∆x2 Case 1 is studied as defined in Section 4.2, where the references ρ(t) and vre f (t) are showed in Figure 4.10 and the associated simulation results are plotted in Figures 4.13 and 4.14. It can be seen in Figure 4.13 that the MPC law tracks the desired tilting angle and level flight speed well, but with certain small oscillations of pitch angle α and vertical speed Vz . This is mainly due to the LPV model formed by interpreting a finite number of linearized models along the transition process, leading to both modeling and nominal control (u0 (ρ)) errors. The associated control efforts are shown in Figure 4.14, where propeller speeds are controlled correspondingly to maintain the altitude during the transition process. Also, based on tilting geometry, the mid-rotor thrust is no longer passing the aircraft center of gravity once the rotor tilts from its vertical position 94 Figure 4.13: Simulation response - Case 1 Figure 4.14: Control efforts - Case 1 95 (ρ = 90 deg), causing aircraft pitch motion. As a result, there is an obvious difference between the front and rear rotor speed from 0 to 35 seconds to reduce the aircraft pitch motion. Then, as for Cases 2 and 3 discussed in Section 4.2, the simulation results are showed in Figure 4.15. The aircraft tracks longitudinal speed Vy for both Cases 2 and 3. However, compared with Case 1, oscillations of pitch angle α and vertical speed Vz become more intensive, which is mainly caused by the fast mid-rotor tilting at the beginning. Rotor speeds and altitude deviation are showed in Figure 4.16, where the active manipulations can be observed from 0 to 10 seconds to stabilize the aircraft during the initial fast tilting. To compare Case 2 with the optimal energy solution (Case 3), both power and energy of Cases 2 and 3 are plotted in solid and dashed curves in Figure 4.17. The total mid-rotor power (summation of two) of Case 3 holds at a higher level in the solid black curve than that of Case 2 in the dashed black one. This is due to relatively larger rotor tilt angle ρ(t) of Case 3 at the beginning of tilting, leading to a larger demand of mid-rotor thrust to achieve the same longitudinal pushing force. Meanwhile, in Case 3, the untilted front and rear rotors are operating at lower power conditions since vertical thrust required is reduced. As shown in the lower plot of Figure 4.17, at 50 sec, the total energy cost of Case 3 is 2.735 kw.h, while that of Case 2 is 2.825 kW.h, a 3.3% reduction. Last but not the least, the disturbance rejection performance has been studied in Cases 4 and 5. The gust disturbance is modeled as an external vertical acceleration d(t) added to the aircraft in the shape of d(t) = K(1 − cos(2πt/T))/2, where K and T denote the magnitude and period of the gust disturbance, respectively. The magnitude K is estimated based on the advisory circular (No. 25-7D) of the U.S. Department of Transportation Federal Aviation Administration, where a gust wind at speed of 12.86m/s (25knot) is considered as the level of atmospheric gust disturbance. The maximum vertical acceleration occurs once the gust wind is in the head-on direction of the aircraft, leading to a sudden increment of wing lift force due to equivalent suddenly-increased airspeed. Based on the target aircraft parameters, magnitude K is estimated to be 4m/s2 with a period of 2 or 4 seconds for Cases 4 and 5, respectively. Based on the normal transition process defined in Case 1, one period of gust disturbance is added at 10, 20, 30, and 40 sec of the transition 96 Figure 4.15: Simulation response - Cases 2 and 3 Figure 4.16: Control efforts - Cases 2 (left) and 3 (right) 97 Figure 4.17: Power and energy comparison for Cases 2 and 3 process, and the simulation results of Cases 4 and 5 are compared in Figure 4.18. Figure 4.18: Simulation response and disturbance of Cases 4 and 5 With four of the gust disturbances kicked in at multiple stages, it can be observed in the left plot of Figure 4.18 that there are obvious oscillations occurred in vertical speed Vz and pitch angle α. And longitudinal speed (Vy ) deviates slightly from its reference, especially when the first disturbance is added for both cases and tracks the reference at the end of transition. It is also worth to mention that the ability to reject gust disturbance is gradually improved when transition 98 Figure 4.19: Control efforts - Cases 4 (left) and 5 (right) from hovering to cruising mode, which can be explained from two aspects. Firstly, from aircraft aerodynamic lift and equivalent disturbance forces showed in the right upper plot of Figure 4.18, the magnitude of disturbance force becomes less significant comparing with the aerodynamic lift forces. Secondly, benefiting from the enhanced aerodynamic force with the increased level flight speed, the elevation actuation sensitivity is reduced, which increases the capability of controlling aircraft pitch motion. To further investigate the control responses, rotor speeds and elevation deflections are showed in Figure 4.19. It can be seen from the responses that rotor speed N reduces each time that the gust disturbance is added to reduce the positive vertical acceleration caused by external disturbance d(t). Also, compared with the four untilted rotors, the speed reduction of the tilted mid-rotors is gradually attenuated from 10s to 40s. This is because of the variation of rotor tilt angle (from ρ(t) = 90 degs to ρ(t) = 0 deg), leading to reduced vertical thrust generated by both tilted rotors and elevation actuation sensitivity. As for the manipulation of elevator θ e , the deflection magnitude at each disturbance rejection stage is gradually reduced with the ongoing transition because of the increased elevator actuation sensitivity explained before. Overall, the designed controller is able to stabilize the aircraft under the two gust disturbance conditions studied, with disturbance rejection performance improved gradually from hovering to cruising mode. In summary, the designed adaptive MPC is able to realize the aircraft tilting transition con- 99 trol with different transition flight envelopes. The overall comparison of five case studies are summarized in Table 4.6 for transition tracking performance, energy consumption, etc. Table 4.6: Performance evaluation of simulation results Case #1 #2 #3 #4 #5 RMSE of Vy 0.2040 0.1506 0.2904 0.2957 0.5571 RMSE of Vz 0.0630 0.0934 0.0748 0.4869 0.6509 RMSE of α 0.1802 0.3024 0.3347 0.4698 0.7169 Total Energy cost (kW.h) 3.243 2.8252 2.7352 3.122 3.044 Peak vertical acceleration (m/s2 ) 0.9153 -0.5501 0.4530 -3.0195 -2.5614 4.5 Conclusion In this chapter, a tilting transition optimization and control for the six-rotor eVTOL (electric vertical takeoff and landing) aircraft with two middle rotors tiltable. Firstly, based on a nonlinear aircraft model developed earlier, a simplified force analysis model is introduced for studying the reference flight states during the tilting transition and an optimization problem is formulated to obtain a transition strategy with minimum total energy. Secondly, a set of linearized models are obtained along with the nominal transition trajectory, based on the nonlinear rigid-body aircraft model, and an affine LPV model is constructed by linearly interpreting the neighboring linearized models. Thirdly, the adaptive MPC method is used for designing the transition control strategy with time-varying constraints, penalty weightings, and feedforward compensations. Five simulation case studies were conducted and simulation results show that the optimal tilting strategy uses minimal transition energy and the MPC scheme is robust to the gust disturbance during the transition. 100 CHAPTER 5 TILT TRANSITION WITH MOTOR FAILURE During the aircraft tilt transition process, motors may be failed due to various reasons, which could lead to dangerous situations. In the hovering control under motor failure, a 1-D control-oriented LPV model has been proposed which is scheduled by the failed motor speed, and the control target is to stabilize the aircraft at the hovering condition. However, in the tilt transition control, the situation becomes even more complicated. Regarding the tracking control targets, there are multiple possibilities when motor failure occurred during transition under different operational conditions. The aircraft can be controlled to track the original nominal transition trajectories, hold at the current condition, or tilt back to the hovering condition. In this research, to study the aircraft dynamics over full transition profile under motor failure, the goal of control is to track the original nominal aircraft transition trajectories when one motor is failed at the beginning of transition. 5.1 2D LPV model of aircraft transition with motor failure Noting that when the aircraft is operated under motor failure condition, the failed rotor speed will deviate from the trimmed nominal operational condition, which means that the locally linearized model is no longer accurate. Following the same principle of generating the 1-D LPV model for normal tilt transition process, one option is still linearizing the nonlinear aircraft model under all operational conditions. However, it requires a large number of offline linearized models with the increased parameter dimension from one to two. For example, consider obtain 20 LTI models of a single rotor failure at each trimmed condition along the transition trajectory with 20 trimmed points, it requires 400 LTI models to be obtained and saved in the real-time controller. Thus, in this section, a novel 2-D LPV model is proposed based on the nominal 1-D LPV tilt transition model with online adaptive parameter. 101 Based on the nonlinear aircraft model, a linearization equation can be expressed as ∂R B (Γ) Û (.)∆βÛ + (.)∆β = ∆ΓÛ + (.) (5.1) ∂ ΓÛ 0 where β is the rigid-body velocity of the aircraft, R B is the loads applied to the aircraft and ΓÛ denotes the speed vector of six rotors. The notation of (.) presents the omitted linearization terms which are not effected by the variation of rotor speed. The load vector R B is constituted as R B = Rgrav + Riner + Rrate + Rgyro + Rprop + Raero , where only terms of gyroscopic load Rgyro and proposive load Rprop are correlated to the propeller speed ΓÛ as below ∂R B (Γ) Û ∂Rgyro (Γ) Û ∂Rprop (Γ)Û = + (5.2) ∂ ΓÛ 0 ∂ ΓÛ 0 ∂ ΓÛ 0 where proposive load is formulated as Rprop = (.)F + (.)T (5.3) Noting that F and T are the thrust and torque applied to rotors which are calculated by F = CT ρAΓÛ 2 R2 (5.4) T = k ΓÛ 2 (5.5) where CT , ρ, A and k are constant coefficients. Based on Eqs. (5.3), (5.4), and (5.5), Rprop can 2 2 ∂Rprop (Γ) Û be converted to a proportional form of ΓÛ as Rprop = K ΓÛ , which means = 2ΓK. Û As ∂Γ Û 0 for gyroscopic load Rgyro , based on the calculation from aircraft specification parameters, it has been observed that the gyroscopic load Rgyro is much smaller than external load Rprop . Thus, by ∂R B (ΓÛ 0 ) omitting the effects from the variation of term Rgyro , the rotor speed dependent term in ∂ ΓÛ 0 Eq. (5.1) can be approximated as ∂R B (Γ) Û ∂R B (ΓÛ 0 ) ΓÛ = ( ) (5.6) ∂ ΓÛ 0 ∂ ΓÛ 0 ΓÛ0 ∂R B (ΓÛ 0 ) where and ΓÛ 0 are the offline linearized system parameter and rotor speed along the ∂ ΓÛ 0 nominal trajectory. In this way, the 2-D LPV model can be formulated based on the 1-D nominal 102 LPV model, with the rotor speed dependent system parameters online variated based on rotor speed vector Γ.Û Converted from the nominal transition LPV system stated in Eq. (4.19), the 2-D LPV model can be formulated as  x(k + 1) = x0 (ρ(k + 1)) + xÛ0 (ρ(k))T + A(ρ(k), Γ(k))∆x(k) Û + B(ρ(k), Γ(k))∆u(k) Û    (5.7)  y(k) = C x(k)     In summary, for modeling the aircraft transition process with motor failure, this 2-D LPV model is constructed based on the nominal transition LPV model scheduled by rotor tilt angle ρ, used as the first scheduling parameter, with part of the system model elements porpotional to failed rotor speed Γ,Û as the second scheduling parameter measured online. 5.1.1 Model error evaluation Recalling the Hσ−2 norm introduced previously and assuming that all the modes are well aligned, the model error between the exact trimmed linear model Gtrim (s) under desired rotor failure condition and the approximated 2-D LPV model G2dl pv can be evaluated by E = Gtrim (s) − G2dl pv (s) σ−2 /kGtrim (s)k σ−2 (5.8) where σ is properly selected to make sure both G2dl pv (s + σ) and Gtrim (s + σ) are stable under the entire tilt-transition operation with rotor failure. Consider the rotor failure of either front-right and mid-left rotor with different rotor speed drop percentage at multiple transition stages (e.g., 85, 50 and 25 deg), the model error evaluation results are shown in Figures 5.1 and 5.2, where the nominal LPV models are also evaluated for comparison. It can be observed that the modeling errors of the 2-D LPV model are below 5% and slightly increased due to deviation of the rotor speed from its nominal condition. 5.1.2 Adaptive MPC formulation for tilt transition with rotor failure Modified from the adaptive MPC structure used in the previous study, the 2-D LPV model based adaptive MPC architecture is shown in Figure 5.3 where the system matrices used for MPC 103 Figure 5.1: Model error evaluation with front-right rotor failure. Figure 5.2: Model error evaluation with mid-left rotor failure. 104 Figure 5.3: Adaptive MPC based on 2D LPV model prediction are online calculated based on scheduling parameter ρ and feedback state signal (rotor speed). ∆xre f and ∆x(k) are the reference and state vectors at current time step k. Denoting e(k) = ∆xre f − ∆x(k), the adaptive MPC design [34, 35, 36] is to find the constrained optimal control law ∆u(k) over a given horizon N that minimizes the constrained quadratic performance defined by N N−1 1 Õ T Õ min { e (k + m)Qe(k + m) + ∆uT (k + m)R∆u(k + m)} ∆u(k),...,∆u(k+N−1) 2 (5.9) m=0 m=0 subject to G∆u(k + m) ≤ h , m = 0, ..., N − 1 This is the same as optimization formulation defined in previous normal transition control, where the control effort ∆u can be solved using QP method. 5.2 Simulation results and discussions This section provides simulation results of tilt transition with front-right or mid-right rotor failure, which covers all single-rotor failure conditions. 5.2.1 Transition with front-right motor failure Firstly, we assumed that the available front-right motor power failed to 50% of the hovering power, which is a 60% motor power loss, assuming peak motor power Pmax is Pmax = 1.25Phover . Based 105 on the same control parameters defined in Table 4.5, the closed-loop simulation results are presented below, where the level flight speed Vy , pitch angle α and vertical speed Vz are shown in Figure 5.4. Figure 5.4: Simulated responses of tilt transition with front-right motor failure It can be observed that the aircraft completes its tilting process and reaches the target level speed at the end of 50 sec. However, the vertical speed drops down at the beginning of the tilting process since the lift thrusts are limited due to the rotor failure. Also, the aircraft is holding a positive pitch angle during the transition, which is because of the dynamic reference compensation developed in Section 4.3.3. With the increased aircraft longitudinal speed, the pitch up motion transfers the aircraft to a climbing pose to compensate the altitude drop. On the other hand, the level flight speed tracking performance is sacrificed relatively due to the climbing mode. The rotor speed and elevator deflection are shown in Figure 5.5. Noting that the vehicle roll and yaw motions can still be balanced with left and right symmetric propeller thrusts. With the front rotor speed drops down due to the rotor failure, the rear rotor speed reduced correspondingly to balance the aircraft pitch motion. The elevator control is activated after 10 sec with the effective aerodynamics generated from the increased level flight speed. 106 Figure 5.5: Control effort of tilt transition with front-right motor failure 5.2.2 Transition with mid-right motor failure Following the same procedure of investigating, the front-right motor failure in the above section, the mid-tight motor failure is studied with the results shown in Figure 5.6. Figure 5.6: Simulated responses of transition with mid-right motor failure 107 At the initial stage, the level flight tracking performance is relatively poor because of the insufficient failed tilt-rotor thrust. After a quick pitching up, the aircraft keeps pitching down from 15 sec to 38 sec which brings the thrust vector of the untilted front and rear rotors toward the front direction to achieve the level flight acceleration. The rotor speed and elevator deflection are shown in Figure 5.7, where the front and rear rotors are operated at higher speed than that without failure to hold the aircraft altitude since the lack of lift force is caused not only by the reduced mid rotor thrust but also the slower longitudinal speed than the reference one, leading to a decreased wing lift. Finally, at 50 sec, the aircraft completes its tilt transition process and reaches the targeted cruising speed. Figure 5.7: Control efforts of tilt-transition with mid-right motor failure 5.3 Conclusion This chapter presents a tilt-transition modeling and control method for a six-rotor eVTOL (elec- tric vertical takeoff and landing) aircraft with two tiltable middle rotors under single rotor failure. Based on the nonlinear aircraft model developed earlier, an online identified 2-D LPV modeling method has been developed, considering the aircraft dynamic characteristics with improved the 108 reference model accuracy. The σ shifted H2 norm is used for evaluating the 2-D LPV model error caused by online updated LPV model, and then, adaptive MPC method is used for designing the transition control strategy with time-varying hard constraints and dynamic reference compensation. Simulation case studies were conducted and show that the 2-D LPV model based MPC scheme is robust to the rotor failure conditions during the tilt-transition. As a summary, this chapter demon- strates the effectiveness of adaptive MPC scheme based on a 2-D online updated LPV model to achieve desired tilt-transition performance under single motor failure. The future work is to study the aircraft failure control with emergent trajectory planning. 109 CHAPTER 6 CONCLUSIONS AND FUTURE WORK 6.1 Conclusions In this dissertation, the urban air mobility (UAM) aircraft research conducted includes control- oriented LPV modeling and adaptive LPV-MPC scheme of hovering and tilt-transition cases with or without motor failure. The main contributions are summarized as follows. (1) The Urban Air Mobility aircraft investigated is a fixed-wing planform with several tilt-rotors built on the wing members or fuselage. An analytical nonlinear rigid-body model is developed to model aeromechanics and flight dynamics of the target tiltrotor eVTOL aircraft. Note that the model is highly nonlinear due to rotor kinematics (rotor spin-rate and tilt angle). For control design, the developed nonlinear model was linearized about a given nonlinear equilibrium condition of interest (such as hovering) or along a desired trajectory for tilt transition to obtain a set of linear state-space models and an LPV (linear parameter-varying) model is obtained by linking these linearized model together using one or two scheduling parameter(s). Then, an adaptive LPV-MPC framework with DRC (direct reference compensation) was developed to handle the nonlinear system dynamics and to achieve desired performance. (2) The UAM aircraft hovering control problem was studied with single rotor failure at multiple possible locations at multiple failure power levels. The adaptive LPV-MPC was applied to hovering failure cases to maintain the aircraft at desired stable hovering condition. The available power level of propeller driven motors is considered in control design using the time- varying input hard constraints. Dynamic reference compensation was also used to improve vehicle recovery performance, with reduced computational cost. The simulation results show that the designed adaptive MPC controller is able to stabilized the vehicle under all single motor failure cases, and the overall RMSE of the vehicle rigid-body velocity oscillations can be 110 significantly reduced using the adaptive LPV-MPC by 8.77% compared with the normal MPC. (3) The UAM aircraft tilt-transition control problem was investigated to achieve smooth transition from hovering to level flight, using the same adaptive LPV-MPC framework developed for hovering failure study. Based on the nonlinear aircraft rigid-body model developed earlier, a simplified force balancing model is developed to find optimal tilt-angle and level flight speed profiles that minimize total tilt-transition energy. Also, a set of linearized models are obtained along with the nominal tilt-transition trajectory based on the nonlinear rigid-body aircraft model, and are used to form the affine LPV model. Simulation case studies were conducted and simulation results validated that the developed optimal tilt-transition strategy uses minimal energy and the adaptive LPV-MPC scheme with DRC is robust to the simulated gust disturbances during the tilt-transition. (4) The UAM aircraft tilt-transition control under motor failure was also conducted using a novel 2- D LPV modeling method that uses the same number of linearized models as that without failure by modifying the 1-D LPV model by adding the second scheduling parameter analytically, which significantly reduces required number of linearized models to be stored in the real-time controller. Note that comparing with the 1-D LPV model developed for tilt-transition without rotor failure, the 2-D LPV model incorporates motor-failure aircraft dynamics to improve the control-oriented model accuracy. A σ shifted H2 norm is used to evaluate the 2-D LPV model errors and to determine required number of linearized models. Simulation case studies were conducted and show that the 2-D LPV model based MPC scheme with DRC is robust to all single rotor failure conditions during the tilt-transition. 6.2 Future work The UAM aircraft tilt-transition process can be further studied considering additional aerody- namic actuators including Aileron and Ruddervators. Also, experimental study shall be conducted with a scaled-down aircraft. 111 APPENDIX 112 APPENDIX ANALYTICAL MODELING OF A UAM AIRCRAFT A.1 Theoretical Formulation Firstly, the rigid-body and tilt-rotor kinematics is considered. A tilt-rotor aircraft is considered as a fixed-wing aircraft (modeled as a rigid body, including fixed lifting surfaces) with tiltrotor(s). The fixed-wing aircraft and tiltrotors are modeled independently and combined in the end. Note that rigid-body aircraft velocity is given by  pÛ B + ω B × p B       vB    β= = (A.1.1)       ωB  θÛ B           where v B and ω B are the translational and angular velocities of the rigid body, respectively. Note that all the kinematic quantities are resolved using the B frame unit vectors in Figure 6.1 that illustrates a single tiltrotor built on a rigid body. The local frame b with the origin at ob is established at a point in the aircraft to define the rotor’s tilt motion, with pb denoting the position vector from O B to ob . Let C Bb be the rotational transformation matrix from b frame to B frame. For a rigid aircraft, both pb and C Bb are constant upon the initialization. Vector pe (more accurately peb , as it is resolved and relative to the b frame) denotes the constant offset between the motor/rotor assembly’s tilting axis located at oe and the b frame. The two vectors pb and pe can be combined for a rigid aircraft. However, they are separated herein for the convenience of accommodating flexible aircraft in future studies. As highlighted in Figure 6.1, the 1-D tilt angle of the motor/rotor assembly is ξ, measured as the rotation of the motor frame e (fixed at the root of pylon) about the b x -axis. The rotor rotation (spin) angle is γ, measured as the rotation of the rotor frame r (fixed at the center of the rotor) about the e y -axis, and d is the radial position of an infinitesimal mass dmr of the rotor. Except for the tilt and spin motions, the motor/rotor assembly is rigid. With these reference frames, we can 113 bz by ξ ob pe A Bz bx oe γ pb By (a) full system Gz Gy vB pB OB ωB Gx OG Bx ez rz ez L ey A γ γ ex A bz ξ, ξ R d dm rx r pe oe ob by (b) tiltrotor (c) rotor Figure 6.1: A free rigid-body with a tilt-rotor define the following vectors, ξ Û 0 0         d                                        Ûξb = 0 , e γÛ = γÛ , e L = L , and d = 0 , r (A.1.2)                                             0 0 0 0                     where the superscripts indicate reference frames in which the vectors are defined and resolved. For simplicity, these superscripts are omitted in the following equations. The transformation matrices between the local frames are  cos γ 0 sin γ  1 0 0         er C (γ) =  0 Cbe (ξ) = 0 cos ξ − sin ξ  . 0  , (A.1.3)     1     − sin γ 0 cos γ  0 sin ξ cos ξ          Let ξÛ and γÛ be the tilt and spin rates, respectively. The rates of these transformation matrices are CÛ er (γ, γ) Û = γC Û̃ er (γ) (A.1.4) Û be (ξ, ξ) C Û = ξC Û̃ be (ξ) where operator [·]˜ refers to a skew symmetric matrix to complete the cross product of two vectors, 114 i.e., a × b = [ã] {b} = [b̃]T {a}. In other words, if a = [a1, a2, a3 ]T , then  0 −a3 a2      ã =  a3 0 −a1  . (A.1.5)     −a2 a1 0      Additional rotational transformation matrix between the body and engine frames can be found by C Be = C Bb Cbe (A.1.6) The angular velocity and acceleration of the rotor (or the rotor frame r) are ωr = ω B + C Bb ξÛ + C Bb Cbe (ξ)γÛ (A.1.7) ωÛr = ω ÛB + C Bb ξÜ + C Bb C Û be (ξ)γÛ + C Bb Cbe (ξ)γÜ Both are already transformed to the B frame. The position of A, assumed to be the mass center of the rotor, is p A(ξ) = p B + p A/B (ξ) (A.1.8) p A/B (ξ) = pb + C Bb p e + C Bb Cbe (ξ)L where p A/B (ξ) is its relative position with respect to the B frame. The velocity of point A is obtained by differentiating its position, leading to       v A = v B + ω B × pb + C Bb pe + C Bb Cbe L + C Bb ξÛ × C Bb Cbe L (A.1.9) A.1.1 Energy and external work of complete aircraft The kinetic energy of the complete aircraft is T = TB + Tr (A.1.10) and the kinetic energy of the rigid-body is 1   1 TB = mB (v B · v B ) + mB v B · ω B × pG/B + ω B · H B (A.1.11) 2 2 where mB is the total mass of the rigid-body, I B the moment of inertia about the B frame, and H B the corresponding angular momentum of the body. Again, all these inertia quantities exclude the 115 contributions of tiltrotor. In addition, one can refer to Figure 6.1 to find the kinetic energy of the tiltrotor, which is given by 1 1 Tr = mr (v A · v A) + ωr · Hr , (A.1.12) 2 2 in which the rotor’s translational and rotational kinetic energy is decoupled since A is considered as the mass center of the rotor, and Hr denotes the angular momentum of the rotor about A, resolved in the B frame, given by Hr = Ir · ωr (A.1.13) where the rotor inertia Ir is firstly calculated in the motor frame (e) and transformed to the body frame (B), and it is given by Ir = C Bb Cbe Ire Ceb CbB (A.1.14) The time derivative of the angular momentum Hr is given by dHr Û = Ir · ωr + Ir · ω Û r + ω B × (Ir · ωr ) (A.1.15) dt where the time derivative of rotor inertia is   IÛr = C Bb ξC Û̃ be (ξ)Ie Ceb (ξ) − Cbe (ξ)Ie Ceb (ξ)ξÛ̃ CbB (A.1.16) r r The potential energy of the complete aircraft comes from the gravity effect and is given by     U = −mB g · p B − pG/B × mB g · θ B − mr g · p B − p A/B × mr g · θ B (A.1.17) The work done by external loads is W ext = Fext ext B · pB + MB · θB (A.1.18) where Fext ext B and M B are the resultant external force and moment, respectively, applied about the B frame from both the aircraft body and tiltrotor. They include the propulsive and aerodynamic loads (as well as loads due to the deployment of wing control surfaces). The reaction forces and moments between the rigid-body aircraft and tiltrotor are internal loads and do not contribute to the formulation. 116 A.1.2 Governing equations of motion By following the Hamilton’s Principle, the governing equation of the complete aircraft is found to be M BB (ξ)βÛ + C BB (β, ξ)β = R B (A.1.19) where the inertia and damping matrices and load vector are      mB I 3 mB p̃TG/B   mr I3 mr p̃TA/B  M BB =  +     m p̃  B G/B I B  m p̃   r A/B r I + m p̃ p̃T  r A/B A/B      (A.1.20)      mB ω̃ B m ω̃ p̃T   mr ω̃ B m ω̃ p̃T B B r B  G/B  +  A/B C BB =       mB p̃G/B ω̃ B ω̃ B I B   r A/B ω̃ B ω̃ B Ir + mr p̃ A/B ω̃ B p̃TA/B    m p̃      R B = Rgrav + Riner + Rrate + Rgyro + Rext where I3 denotes a 3-by-3 identity matrix. The gravitational acceleration g is resolved in the B frame. The term Ir + mr p̃ A/B p̃TA/B in M BB is nothing but the rotor’s mass moment of inertia transferred to the B frame, after applying the Parallel Axis Theorem. Equation (A.1.19) is sufficient to solve for the aircraft’s dynamic response. However, for the convenience of tracking the propagation of the aircraft body frame, additional kinematic equations are established. Given the rigid-body angular velocity ω B = [ω B1, ω B2, ω B3 ]T , one may use quaternions (ζ = [ζ0, ζ1, ζ2, ζ3 ]T ) to describe the rigid-body orientation as follows, 1 ζÛ = − Ωζ (β)ζ (A.1.21) 2 where  0 ω ω ω    B1 B2 B3     −ω  B1 0 −ω B3 ω B2  Ωζ (β) =   (A.1.22) −ω  B2 ω B3 0 −ω  B1    −ω B3 −ω B2 ω B1 0      117 The rotational transformation matrix C BG is also given in terms of the quaternions, i.e.,  2 ζ + ζ 2 − ζ 2 − ζ 2 2 (ζ1 ζ2 + ζ0 ζ3 ) 2 (ζ1 ζ3 − ζ0 ζ2 )    0 1 2 3 C BG (ζ) =  2 (ζ1 ζ2 − ζ0 ζ3 ) ζ02 − ζ12 + ζ22 − ζ32 (A.1.23)   2 (ζ2 ζ3 + ζ0 ζ1 )     2 (ζ1 ζ3 + ζ0 ζ2 )  2 (ζ2 ζ3 − ζ0 ζ1 ) 2 2 2 ζ0 − ζ1 − ζ2 + ζ3  2    which leads to the inertial translational velocity of the B frame: h i pÛ G B = CGB v B = CGB 03 β (A.1.24)  T where CGB = C BG . Finally, the full set of nonlinear flight dynamic governing equations is described by M BB (ξ)βÛ + C BB (β, ξ)β = R B 1 ζÛ = − Ωζ (β)ζ (A.1.25) 2 h i G pÛ B = C (ζ) 03 β GB As a post-processing, the Euler angles can be calculated from the quaternions. A.1.3 Inertial and gyroscopic loads Among the components of the load vectors in Eq. (A.1.20), the gravity load is  Fgrav         I   I 3 3   Rgrav = =  mB g + (A.1.26)     mr g       Mgrav    p̃G/B  p̃  A/B  (ξ)        The inertial load (Riner , Rrate ) and gyroscopic load (Rgyro ) are all originated from the kinetic energy. The mass of rotor may have inertial force due to tilt-acceleration, centrifugal force, and Coriolis force, given by Finer = −mr C Bb ξÜ × Cbe L + C Bb ξÛ × ξÛ × Cbe L  h  i n h  i o n h   i o + 2ω B × C Bb ξÛ × Cbe L (A.1.27) More inertial loads, coming from the accelerated rotation of the rotor considered as a rigid body, is given by Mracc = −Ir · C Bb ξÜ + C Bb Cbe γÜ   (A.1.28) 118 Therefore, the equivalent inertial load is  Finer        iner   I  3  iner   03×1   = = + (A.1.29)      R F   Miner    p̃ A/B   Mracc          In addition, the orientation of the moment of inertia varies in time for a tilting rotor, with its rate evaluated in the body frame given by Eq. (A.1.16). Such a change of inertia generates a moment, given by    03×1   Rrate = (A.1.30)     Mrate      where Mrate = −IÛr · ω B + C Bb ξÛ + C Bb Cbe γÛ   (A.1.31) Finally, the gyroscopic load is    03×1   Rgyro = (A.1.32)     Mgyro      where Mgyro = −ω B × Ir · C Bb ξÛ + C Bb Cbe γÛ − Ir · C Bb ξC h  i   Û̃ be γÛ (A.1.33) A.1.4 Propulsion In this study, Rext in Eq. (A.1.20) includes both propulsive and aerodynamic loads. For simplicity, the static thrust produced by a single rotor is πN 2 2   F = CT ρA R = CT ρAγÛ 2 R2 (A.1.34) 30 The aerodynamic torque applied on the rotor is T = (−1)n−1 k γÛ 2 (A.1.35) where k is a constant determined by the air density, geometry of the rotor, etc. The thrust and torque are aligned along the e y -axis, whose orientation changes with the rotor’s tilt angle. Hence, 119 the propulsive load in the body frame B is 0 0             prop               prop  F   I 3  Be   0   3  Be     = =  C (ξ) F + (A.1.36)           R   C (ξ) T    Mprop    p̃ A/B (ξ)       I   3       0 0                   A.1.5 Aerodynamics The aerodynamic load is firstly calculated at each section along the span of lifting surfaces. In doing so, the lifting surfaces are meshed with 1-D elements along their span. Referring to Figure 6.2 for airfoil velocity components, the 2-D quasi-steady aerodynamic loads on each thin airfoil section undergoing arbitrary motions in an incompressible inviscid subsonic flow are described by 1 αÛ     zÛ 2 lac = πρb (−Üz + yÛ αÛ − d α) 2 Ü + 2πρb yÛ − + b − d + ρb yÛ 2 clδ δ yÛ 2 yÛ 1 1 1     mac = πρb3 zÜ − yÛ αÛ − b − d αÜ + 2ρb2 yÛ 2 cmδ δ (A.1.37) 2 8 2 αÛ 2   2 zÛ dac = −2πρb yÛ − − d + ρb yÛ 2 cdδ δ yÛ yÛ where δ is the trailing-edge flap deflection angle, b the semichord, and d the distance of midchord in front of the reference axis. The sectional aerodynamic loads are transformed to the body frame B, given by oT f aero = C Ba1 0 dac lac n  T (A.1.38) aero = C 1 mac + lac 1 b + d 0 0   m Ba 2 Bz lac a1z a0z By Bx mac y a0y α a.c. dac V∞ a1y z zero-lift line d b b Figure 6.2: Aerodynamic frame and load components 120 The distributed aerodynamic force and moments are further integrated to obtain the resultant aerodynamic load about the B frame, i.e.,  E f f aero     Raero = (A.1.39)     Em maero      with E f and Em being the influence matrices determined by the numerical integration scheme. A.1.6 Linearized equation and state-space equation To facilitate the process, one can define the following control inputs that include rotor tilt and spin kinematics, respectively: n o ΞT = ξ1 · · · ξn n o (A.1.40) T Γ = γ1 · · · γ n Note their difference between the vectors defined in Eq. (A.1.2). The Taylor’s expansion is utilized to the governing equation about a given nonlinear equilibrium ( 0 ) and the corresponding control input (u0 ) given by  T T  T  T  T0 = βT0 βÛ 0 ζT0 ζÛ 0 pG pÛ G B,0 B,0 n o (A.1.41) T T T T T u0 = δT δT δT ΞT Ξ Û0 Ξ Ü 0 ΓÛ 0 ΓÜ 0 e,0 a,0 r,0 0 After the expansion, high order terms are neglected, and terms with the same variables are grouped together, leading to the linearized equation of motion as follows, ! ∂R B ∂C ∂R   BB B M BB | 0 − ∆βÛ + C BB | 0 + β − ∆β ∂ βÛ ∂β 0 0 ∂β 0 0 ∂R B ∂R B ∂R B ∂R B = ∆ζ + ∆δe + ∆δa + ∆δr ∂ζ 0 ∂δe 0 ∂δa 0 ∂δr 0 ∂R B ∂M BB Û ∂C BB   + − β − β ∆Ξ ∂Ξ 0 ∂Ξ 0 0 ∂Ξ 0 0 (A.1.42) ∂R B Û + ∂R B ∆Ξ Ü + ∂R B ∆ΓÛ + ∂R B ∆ΓÜ + ∆Ξ ∂Ξ 0 Û ∂Ξ 0 Ü ∂ ΓÛ 0 ∂ ΓÜ 0 1 ∂Ωζ (β) 1 ∆ζÛ = − ∆βζ0 − Ωζ (β0 )∆ζ 2 ∂β 0 2 ∂CGB (ζ)   h i G ∆pÛ B = ∆ζ 03 β0 + C (ζ0 ) 03 ∆β GB ∂ζ 0 121 where (·)| 0 means a quantity evaluated at the given nonlinear equilibrium. One can omit the ∆ sign in Eq. (A.1.42), and derive the state-space form of the linearized equation of motion as xÛ = Ax + Bu (A.1.43) where the state variable and control input are given by   T  T x = β ζ T T pGB h i (A.1.44) uT = δTe δTa δrT ΞT Ξ ÛT ÜT Ξ ΓÛT ÜT Γ A.1.7 Trim solutions Two types of trim problems are considered in this study. The first is to trim for a steady longitudinal flight with a constant speed, including both steady-level flight and constant climbing, in which the equilibrium equation is reduced to R B = 0. In these flights, one only has to satisfy three components in the equilibrium equation. Therefore, the cost function is given by n o ST = R B,2 R B,3 R B,4 (A.1.45) and the trim valuables are n o qT = αB δe N (A.1.46) where the body incidence angle αB is equal to the pitch angle. The solution starts from an initial guess of trim valuables q0 . The Newton-Raphson method is used to search for the solution of q that satisfies S(q) = 0. At each step, the search direction is determined by   −1 ∂S ∆qi = − S ∂q i i (A.1.47) qi+1 = qi + ∆qi ∂S in which the Jacobian is evaluated at each step i using the finite difference method. The search ∂q continues until a user-defined tolerance is satisfied, i.e., ||S(q)||2 ≤ εtrim (A.1.48) 122 Another trim problem is to find the required rotor speed during the vertical take-off. The aerodynamics does not participate in the solution, and the body incidence angle can be reasonably set to zero. Therefore, the trim variable consists of rotor speeds only, i.e., n o qT = N1 N2 · · · Nn (A.1.49) whereas the cost function of trimming for the vertical flight becomes n o ST = R B,3 − (mB + mr ) az R B,4 R B,5 R B,6 (A.1.50) where az is the required vertical acceleration during the take-off. The approach to trim for the vertical flight is the same as level flight. However, the Jacobian is usually not a square matrix as the number of rotors being used for trim might be greater than 3. In this case, the Moore-Penrose pseudo-inverse of the Jacobian is used to search for the feasible trim variables. Usually, the solution is not unique if one starts from different initial conditions. The preferred solution should result in the minimum infinity-norm of q (min ||q||∞ ), which corresponds to the lowest possible peak rotor speed. If any rotor has a malfunction, the corresponding rotor speed can be set to zero or reduce to a specific value. It is then removed from the list of trim valuables from Eq. (A.1.49) during the search. A.2 Trim Studies Figure 2.7 shows the geometry of an urban air mobility aircraft model of interest, with 6 rotors at 90 deg orientation and the fuselage not shown. The dashed-dotted line indicates the reference line of lifting surfaces, which is 25% the chord from the leading edge. Both the main wing and v-tail are modeled as rigid lifting surfaces. Table 2.1 lists the inertial parameters of rigid-body, fixed-wing aircraft and the rotors, and Table 2.2 lists the aerodynamic parameters of the main wing and v-tail. The main wing contains an aileron, as indicated in Figure 2.7, which occupies 25% of the chord, running from 60% to 90% of the wing span. A ruddervator is defined on the v-tail, occupying 25% of the chord, running from 40% to 80% of the tail span. 123 1200 Rotor speed, rpm 1100 1000 900 0 1 2 3 4 5 6000 Rotor thrust, N 5000 4000 3000 0 1 2 3 4 5 Vertical acceleration, m/s 2 Figure 6.3: Rotor speed and thrust during nominal vertical takeoff (identical speed and thrust for all rotors) A.2.1 Trim for vertical takeoff Figure 6.3 plots the trimmed rotor speed and thrust for the vertical takeoff if all 6 rotors work under the nominal conditions, where all rotors have the same speed and thrust. The plot can be used to find the required rotor speed and thrust during the initial vertical takeoff acceleration, as well as the continuing flight at a constant vertical speed. As only the static thrust is considered in the study, the trim solution turns out to be independent of the takeoff speed. Next, rotor 1 is considered malfunctioned. When this happens, the rotor speed is set to zero. Hence, there is no thrust output from that rotor. For the aircraft to continue the takeoff at the required acceleration, other rotors’ speed must be adjusted, so the thrusts are re-distributed to maintain the vertical flight. Figure 6.4 shows the re-distributed rotor speed and thrust. As can be seen, the speed of rotor 6 is also knocked down to zero, while the reminders speed up by 22.47%, resulting in a 50% increase in individual thrust. However, the total thrust remains unchanged, as required by trim. Figure 6.5 shows the trim result when rotor 3 has a malfunction. The trim result indicates the same quantitative change in rotor outputs to balance the vehicle and maintain the vertical flight. However, the distribution of rotor speed and thrust among the specific functioning rotors is different from the previous case. The following studies will only consider rotor 1 being nonfunctional. 124 1600 9000 Rotor 2, 3, 4, & 5 1400 8000 Rotor 2, 3, 4, & 5 7000 1200 6000 Rotor speed, rpm Rotor thrust, N 1000 5000 800 Rotor 1 4000 Rotor 1 600 Rotor 2 Rotor 2 Rotor 3 3000 Rotor 3 400 Rotor 4 Rotor 4 Rotor 5 2000 Rotor 5 200 Rotor 1 & 6 Rotor 6 1000 Rotor 1 & 6 Rotor 6 All normal All normal 0 0 0 1 2 3 4 5 0 1 2 3 4 5 2 2 Vertical acceleration, m/s Vertical acceleration, m/s Figure 6.4: Individual rotor speed and thrust if rotor 1 rpm is zero 1600 9000 Rotor 1, 2, 5, & 6 1400 8000 Rotor 1, 2, 5, & 6 7000 1200 6000 Rotor speed, rpm Rotor thrust, N 1000 5000 800 Rotor 1 4000 Rotor 1 600 Rotor 2 Rotor 2 Rotor 3 3000 Rotor 3 400 Rotor 4 Rotor 4 Rotor 5 2000 Rotor 5 200 Rotor 3 & 4 Rotor 6 1000 Rotor 3 & 4 Rotor 6 All normal All normal 0 0 0 1 2 3 4 5 0 1 2 3 4 5 Vertical acceleration, m/s 2 Vertical acceleration, m/s 2 Figure 6.5: Individual rotor speed and thrust if rotor 3 rpm is zero A.2.2 Trim for level transition flight In this section, the aircraft is already bright up to the full altitude 304.8m after the vertical takeoff. It should then transit from the vertical takeoff to level flight. It is a transient process, with rotors tilting from upright to forward (horizontal) position. During the transition, the aircraft’s level flight speed is increased from zero to 68 m/s. For simplicity, the transition is assumed to be completed as a quasi-steady trimmed process. The body pitching and ruddervator deflection, along with rotor tilt angle and spin rate, are considered trim variables in the process. Only the middle two rotors (3 and 4) are tilted, as two rotors are sufficient to provide the thrust once they are in the horizontal 125 position. In addition, the remaining vertical rotors can provide a significant amount of lift during the transition. Figure 6.6 shows the aircraft body pitch and ruddervator deflection during level flight transition. Figures 6.7 and 6.8 illustrate the rotors’ tilt angle and spin rate, respectively, during the transition. Rotor kinematics is symmetric during the transition. Therefore, the rotors not shown in Figures 6.7 and 6.8 should take negative values. Note that the spin rates of the non-tilting rotors are specified, such that the trim results return small body pitch angles. This is for the consideration of the ride quality of the vehicle. Further optimization study can be carried out to find the rotor kinematics with minimum power consumption during the transition. 0.4 Body pitch, deg 0.3 0.2 0.1 0 0 10 20 30 40 50 60 70 Ruddervator def., deg 4 0 -4 -8 -12 0 10 20 30 40 50 60 70 Level flight speed, m/s Figure 6.6: Trimmed aircraft body pitch and ruddervator deflection during level flight transition 126 100 80 Rotor tilt angle, deg 60 40 20 Rotor 1 Rotor 4 Rotor 5 0 0 10 20 30 40 50 60 70 Level flight speed, m/s Figure 6.7: Trimmed rotor tilt angle during level flight transition 1200 Rotor 1 Rotor 4 1000 Rotor 5 Rotor spin rate, rpm 800 600 400 200 0 0 10 20 30 40 50 60 70 Level flight speed, m/s Figure 6.8: Trimmed rotor spin rate during level flight transition 127 BIBLIOGRAPHY 128 BIBLIOGRAPHY [1] R. Binder, L. Garrow, B. German, P. Mokhtarian, M. Daskilewicz, and T. Douthat, “If you fly it, will commuters come? predicting demand for evtol urban air trips,” in AIAA Conference, Atlanta, Georgia, 2018, pp. 1–41. [2] G. E. Dorrington et al., “Performance of electric vertical take-off and landing (evtol) hovering craft,” in AIAC18: 18th Australian International Aerospace Congress (2019): HUMS-11th Defence Science and Technology (DST) International Conference on Health and Usage Mon- itoring (HUMS 2019): ISSFD-27th International Symposium on Space Flight Dynamics (ISSFD). Engineers Australia, Royal Aeronautical Society., 2019, p. 84. [3] I. C. Kleinbekman, M. A. Mitici, and P. Wei, “evtol arrival sequencing and scheduling for on-demand urban air mobility,” in 2018 IEEE/AIAA 37th Digital Avionics Systems Conference (DASC). IEEE, 2018, pp. 1–7. [4] B. German, M. Daskilewicz, T. K. Hamilton, and M. M. Warren, “Cargo delivery in by passenger evtol aircraft: A case study in the san francisco bay area,” in 2018 AIAA Aerospace Sciences Meeting, 2018, p. 2006. [5] P. Pradeep and P. Wei, “Heuristic approach for arrival sequencing and scheduling for evtol aircraft in on-demand urban air mobility,” in 2018 IEEE/AIAA 37th Digital Avionics Systems Conference (DASC). IEEE, 2018, pp. 1–7. [6] M. Skuhersky, “A first-principle power and energy model for evtol vehicles,” Ph.D. disserta- tion, 2019. [7] P. Pradeep and P. Wei, “Energy efficient arrival with rta constraint for urban evtol operations,” in 2018 AIAA Aerospace Sciences Meeting, 2018, p. 2008. [8] M. Daskilewicz, B. German, M. Warren, L. A. Garrow, S.-S. Boddupalli, and T. H. Douthat, “Progress in vertiport placement and estimating aircraft range requirements for evtol daily commuting,” in 2018 Aviation Technology, Integration, and Operations Conference, 2018, p. 2884. [9] P. Payuhavorakulchai, “Cost analysis of evtol configuration design for air ambulance in japan,” 2019. [10] J. M. Vegh, E. Botero, M. Clarke, J. Smart, and J. Alonso, “Current capabilities and challenges of ndarc and suave for evtol aircraft design and analysis,” in AIAA Propulsion and Energy 2019 Forum, 2019, p. 4505. [11] S.-S. Boddupalli, “Estimating demand for an electric vertical landing and takeoff (evtol) air taxi service using discrete choice modeling,” Ph.D. dissertation, Georgia Institute of Technology, 2019. 129 [12] T. Lombaerts, J. Kaneshige, S. Schuet, G. Hardy, B. L. Aponso, and K. H. Shish, “Nonlinear dynamic inversion based attitude control for a hovering quad tiltrotor evtol vehicle,” in AIAA Scitech 2019 Forum, 2019, p. 0134. [13] P.-M. Basset, B. D. Vu, P. Beaumier, G. Reboul, and B. Ortun, “Models and methods at onera for the presizing of evtol hybrid aircraft including analysis of failure scenarios,” 2018. [14] H. Yeo and W. Johnson, “Assessment of comprehensive analysis calculation of structural loads on rotors,” in American Helicopter Society 60th Annual Forum Proceedings. American Helicopter Society, 2004. [15] F. Kendoul, I. Fantoni, and R. Lozano, “Modeling and control of a small autonomous aircraft having two tilting rotors,” IEEE Transactions on Robotics, vol. 22, no. 6, pp. 1297–1302, 2006. [16] O. A. Bauchau and N. K. Kang, “A multibody formulation for helicopter structural dynamic analysis,” Journal of the American Helicopter Society, vol. 38, no. 2, pp. 3–14, 1993. [17] W. Johnson, “Rotorcraft dynamics models for a comprehensive analysis,” in American Heli- copter Society 54th Annual Forum Proceedings. American Helicopter Society, 1998. [18] A. Abhishek, A. Datta, and I. Chopra, “Prediction of uh-60a structural loads using multibody analysis and swashplate dynamics,” in American Helicopter Society 62nd Annual Forum Proceedings. American Helicopter Society, 2006. [19] D. H. Hodges, H. Saberi, and R. A. Ormiston, “Development of nonlinear beam elements for rotorcraft comprehensive analyses,” Journal of the American Helicopter Society, vol. 52, no. 1, pp. 36–48, 2007. [20] S. Lim, “Analysis and control of linear parameter-varying systems,” Ph.D. dissertation, Cite- seer, 1999. [21] J. Mohammadpour and C. W. Scherer, Control of linear parameter varying systems with applications. Springer Science & Business Media, 2012. [22] T. He, G. G. Zhu, S. S.-M. Swei, and W. Su, “Smooth-switching lpv control for vibration suppression of a flexible airplane wing,” Aerospace Science and Technology, vol. 84, pp. 895–903, 2019. [23] T. He, G. G. Zhu, and S. S.-M. Swei, “Smooth switching lpv dynamic output-feedback control,” International Journal of Control, Automation and Systems, vol. 18, no. 6, pp. 1367–1377, 2020. [24] T. He, G. G. Zhu, S. S.-M. Swei, and W. Su, “Optimal sensor placement for flexible wings using the greedy algorithm,” in 2020 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM). IEEE, 2020, pp. 577–582. [25] T. He, G. G. Zhu, and X. Chen, “A two-step lmi scheme for h2-hinf control design,” in 2020 American Control Conference (ACC). IEEE, 2020, pp. 1545–1550. 130 [26] T. He, A. K. Al-Jiboory, G. G. Zhu, S. S.-M. Swei, and W. Su, “Application of icc lpv control to a blended-wing-body airplane with guaranteed hinf performance,” Aerospace Science and Technology, vol. 81, pp. 88–98, 2018. [27] S. Qu, T. He, and G. G. Zhu, “Engine egr valve modeling and switched lpv control considering nonlinear dry friction,” IEEE/ASME Transactions on Mechatronics, vol. 25, no. 3, pp. 1668– 1678, 2020. [28] H. Fukushima, T.-H. Kim, and T. Sugie, “Adaptive model predictive control for a class of constrained linear systems based on the comparison model,” Automatica, vol. 43, no. 2, pp. 301–308, 2007. [29] J.-S. Kim, “Recent advances in adaptive mpc,” in ICCAS 2010. IEEE, 2010, pp. 218–222. [30] J. S. Shamma, “An overview of lpv systems,” in Control of linear parameter varying systems with applications. Springer, 2012, pp. 3–26. [31] S. Qu, T. He, and G. G. Zhu, “Lpv modeling and switched control for egr valves with dry friction,” in 2019 IEEE Conference on Control Technology and Applications (CCTA). IEEE, 2019, pp. 400–405. [32] B. Ding, “Dynamic output feedback mpc for lpv systems via near-optimal solutions,” in Proceedings of the 30th Chinese control conference. IEEE, 2011, pp. 3340–3345. [33] A. Casavola, D. Famularo, and G. Franze, “A feedback min-max mpc algorithm for lpv systems subject to bounded rates of change of parameters,” IEEE Transactions on Automatic Control, vol. 47, no. 7, pp. 1147–1153, 2002. [34] E. F. Camacho and C. B. Alba, Model predictive control. Springer Science & Business Media, 2013. [35] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000. [36] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control: Theory and practice—a survey,” Automatica, vol. 25, no. 3, pp. 335–348, 1989. [37] W. Su, S. Qu, G. G. Zhu, S. S.-M. Swei, M. Hashimoto, and T. Zeng, “A control-oriented dynamic model of tiltrotor aircraft for urban air mobility,” in 2021 AIAA SciTech Forum, 2021, p. 0091. [38] A. K. Al-Jiboory, G. Zhu, S. S.-M. Swei, W. Su, and N. T. Nguyen, “Lpv modeling of a flexible wing aircraft using modal alignment and adaptive gridding methods,” Aerospace science and technology, vol. 66, pp. 92–102, 2017. [39] S. Gumussoy and W. Michiels, “Computing hi nf norms of time-delay systems,” in Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference. IEEE, 2009, pp. 263–268. 131 [40] C. Schmid and L. T. Biegler, “Quadratic programming methods for reduced hessian sqp,” Computers & chemical engineering, vol. 18, no. 9, pp. 817–832, 1994. [41] T. H. Phuong, M. P. Belov, and D. Van Thuy, “Adaptive model predictive control for nonlinear elastic electrical transmission servo drives,” in 2019 IEEE Conference of Russian Young Researchers in Electrical and Electronic Engineering (EIConRus). IEEE, 2019, pp. 704– 708. [42] R. K. Mehra, R. K. Prasanth, and S. Gopalaswamy, “XV-15 tiltrotor flight control system design using model predictive control,” in 1998 IEEE Aerospace Conference, Snowmass, CO, Mar. 1998. [43] R. T. Rysdyk and A. J. Calise, “Adaptive model inversion flight control for tilt-rotor aircraft,” Journal of guidance, control, and dynamics, vol. 22, no. 3, pp. 402–407, 1999. [44] Z. Lyu, Z. Wang, D. Duan, L. Lin, J. Li, Y. Yang, Y. Chen, and Y. Li, “Tilting path optimization of tilt quad rotor in conversion process based on ant colony optimization algorithm,” IEEE Access, vol. 8, pp. 140 777–140 791, 2020. [45] J. Gao, Q. Zhang, J. Chen, and X. Wang, “Take-off trajectory optimization of tilt-rotor aircraft based on direct allocation method,” in IOP Conference Series: Materials Science and Engineering, vol. 768, no. 4. IOP Publishing, 2020, p. 042004. [46] X. Wang and L. Cai, “Mathematical modeling and control of a tilt-rotor aircraft,” Aerospace Science and Technology, vol. 47, pp. 473–492, 2015. [47] E. Cetinsoy, S. Dikyar, C. Hançer, K. Oner, E. Sirimoglu, M. Unel, and M. Aksit, “Design and construction of a novel quad tilt-wing uav,” Mechatronics, vol. 22, no. 6, pp. 723–745, 2012. [48] W. Su, S. Qu, G. G. Zhu, S. S.-M. Swei, M. Hashimoto, and T. Zeng, “A control-oriented dynamic model of tiltrotor aircraft for urban air mobility,” in AIAA Scitech 2021 Forum, 2021, p. 0091. 132