Simulation and Control of Deformable Autonomous Airships in Turbulent Wind

Abstract. Fixed wing and multirotor UAVs are common in the field of robotics. Solutions for simulation and control of these vehicles are ubiquitous. This is not the case for airships, a simulation of which needs to address unique properties, i) dynamic deformation in response to aerodynamic and control forces, ii) high susceptibility to wind and turbulence at low airspeed, iii) high variability in airship designs regarding placement, direction and vectoring of thrusters and control surfaces. We present a flexible framework for modeling, simulation and control of airships, based on the Robot operating system (ROS), simulation environment (Gazebo) and commercial off the shelf (COTS) electronics, both of which are open source. Based on simulated wind and deformation, we predict substantial effects on controllability, verified in real world flight experiments. All our code is shared as open source, for the benefit of the community and to facilitate lighter-than-air vehicle (LTAV) research. https://github.com/robot-perception-group/airship_simulation


Introduction
In the last decade, the term Unmanned Aerial Vehicles (UAV) has become near-synonymous with so called "drones"small lightweight multirotor craft, kept stable by active control using cheap lightweight MEMS sensors that became ubiquitous as a byproduct of the smartphone revolution [1]. Their popularity is mirrored in the field of aerial robotic research. The vast majority of works cover multirotor craft, and to a lesser degree, fixed wing aircraft and helicopters. Thus, flight electronics, control software and simulation tools are readily available for these types of vehicles [2,3]. This allows researchers to build up a suitable autonomous flying platform for any research task with ease, using commercial off the shelf (COTS) components and open source software.
In contrast, lighter-than-air vehicles (LTAV) such as rigid and nonrigid dirigibles ( Fig. 1) have fallen out of fashion, despite their superiority for some applications [4]. LTAV are uniquely suited as mobile aerial communication relays, for monitoring and wildlife conservation tasks. They produce little noise, have high energy efficiency and long flight times, display benign collision and crash characteristics, pose low danger and cause little environmental impact. However, their comparably high handling complexity, size, lifting gas requirements and cost create an entry barrier for researchers. Unlike for heavier-than-air "drones", there have been no COTS flight controllers that support autonomous dirigible flight. Therefore, guidance and control algorithms have to be implemented for each vehicle -even though various suitable control strategies can be found in the literature [4,5,6,7,8,9,10]. Similar to both rotor craft and fixed wing UAVs, dirigibles come in many types of actuator arrangements: Fixed or vectoring main thrusters, differential thrust, different tail fin arrangements and auxiliary thrusters, single or double hull, etc. Thus, a control algorithm for a specific vehicle might not always be applicable to others.
Development of control algorithms is further complicated by the lack of realistic simulation. Existing robotic simulation environments [2] are often ill-suited for the characteristics of LTAV. As LTAV are highly susceptible to wind and turbulence, it is crucial to accurately model both aerodynamic forces and the effects these have on an inherently non-rigid vehicle. To our knowledge, no existing realtime robotic simulation environment addresses this, although analytic studies of these effects have been conducted in the past [11].
In this work, we attempt to solve the aforementioned issues in the context of LTAV and provide researchers in the field of aerial robotics with the tools they are accustomed to when working with fixed wing and rotor craft. Our contributions, for which all source code is provided 3 , are: 1. A robotic simulation framework, using ROS/Gazebo [2,12] for real time hardware in the loop (HITL) and software in the loop (SITL) simulation of easily customizable airship models in a realistic virtual environment.
(a) Realistic simulation of wind, turbulence [13] and aerodynamic forces on a custom shaped, modular, non-rigid deformable air frame. (b) Simulated deformation of the air frame in response to aerodynamic, control, thrust and collision forces. (c) Simulation of buoyancy variation, changes in rigidity and shape, depending on hull pressure and structural parameters. 2. A baseline guidance and control algorithm for autonomous navigation of airships using hierarchical PI controllers. We employ a generic design, to operate on a large subset of possible airship shape and actuator configurations.
(a) Implemented based on the Librepilot [2,14] open source flight control tool chain, to integrate with existing SITL & HITL frameworks and available COTS flight control hardware as well as the above mentioned simulation framework. (b) Full integration with ROS [12]. This facilitates easy integration of lighter-than-air vehicles into existing and future robotic research projects for higher level tasks. 3. Through simulation experiments, we show the effect of changes in rigidity and wind turbulence on the vehicle's controllability, including the introduction of oscillation modes not present in rigid models. 4. Through real flight experiments, we demonstrate the flight behavior of a 5m autonomous blimp UAV, with control coefficients previously determined in our simulation, thus validating our simulation results.
In Sec. 2, we compare our work to the state of the art. We explain our methodology in Sec. 3. Experiments and results are shown in Sec. 4, followed by a remark on limitations in Sec. 5 and our summarizing conclusions in Sec. 6.

State of the art
The seminal work involving the AURORA airship [4,10,15] describes the development of a complete system, including novel on board electronics, communication, control algorithms and a flight simulator, both for control evaluation and human operator training. The hierarchical PI control algorithm described there is similar to our implementation, but lacks utilization of dynamic thrust vector control. It also makes many simplifications, including flight at constant airspeed. Although a separate hover control scheme is mentioned, which does use thrust vector control, this is not described in detail. Their simulation assumes rigidity of the air frame. Aerodynamics are modeled based on existing wind tunnel data from a different airship with similar proportions. In contrast, we use modern COTS UAV components and freely available open source software, which greatly enhances reproducibility and the ability to integrate with common robotic components [12]. We also model nonrigid deformation effects in our simulation. Our aerodynamic approach is more generic and allows modeling arbitrary dirigibles, even when no detailed aerodynamic measurements are available -based on geometric shape.
The literature review paper [11] covers modeling of airships, including analysis of deformation and its effect on aerodynamics, spanning 100 years of research, from wind tunnel tests in the early 20th century [16] to computational fluid dynamics (CFD) analysis in the 2000s [17]. We take inspiration from some of the presented methods, especially [18], which we apply to real time simulation of aerodynamically relevant components. We consider these rigid individually, but allow non-rigid interactions on joints between them. Our simulation does not take aerodynamic cross-interactions between individual components into account. However, we argue that these can be approximated by artificially adjusting lift and drag coefficients to minimize the overall modeling error. This could be done, for example, by comparison to CFD analysis of the rigid shape. However, this is not covered by our current work.
More recent works on airships often involve coupling control with perception and computer vision, such as [5,9,19,20]. All of these works make simplifying assumptions, such as frame rigidity and disregard of wind turbulence.
A novel airship guidance and control scheme was presented in [8], using a backstepping control strategy. Its performance is mathematically proven utilizing a highly simplified analytical plant model under the assumption of rigidity. No real world experiments have been made. This has inspired a number of more recent theoretical works on various airship guidance schemes [21,22]. However, none of them seem to have been validated in real-world flight. In contrast, we base our model on observed real world flight behavior and validate our results in additional real world flights.
The airship control problem can also be addressed with a model free or learned-model approach [6,7]. These techniques are well suited to deal with the dynamics of a deforming vehicle under changing conditions from a control perspective, and can be directly applied to the real world. However, a simulation, as presented in our work, is still required to validate and compare the performance of learned algorithms under controlled, reproducible conditions.

Simulation
Gazebo is a physics simulator often used in robotics, which seamlessly integrates with the Robot Operating System (ROS) [12]. Gazebo models rigid physical objects (primitives) based on shape, mass and moments of inertia. Using an XML/URDF description file, these component objects can be described and linked together into hierarchical compound objects. Rigid joints are simulated, using very stiff springs with high damping. These joint parameters can be altered for both linear and rotational freedom. Motion limits, stiffness and dampening properties can be set to simulate non-rigid compound objects. The simulator employs a solver engine to calculate object motions under Newtonian physics, taking all joint forces, collisions and gravity into account. This can be extended through plugins that add additional forces to component objects or joints to model actuators or thrusters. The simulator models time in discrete time steps with configurable temporal resolution and can handle most scenarios in real time with sufficient accuracy. Simulated sensors, which provide simulation data measurements for IMU and positioning sensors as ROS messages are also provided by plugins attached to their respective component objects.
We model a blimp (Fig. 2) as a flexible compound object, consisting of rigid primitive component objects for all actuated and non-actuated components. These are fins, rudders, the gondola, thrusters, ballast weights, etc. We model the main hull with at least two separate components, to allow flexibility of the shape as well as higher granularity of forces during rotations. Plugins are added to these primitives to model buoyancy, aerodynamic forces, control actuation and thrust. This setup allows the physical properties of the blimp to be modeled as the sum of its physical components. Mass, dimensions and inertial moments of homogeneous basic shapes can be easily determined by first principles or by measuring the corresponding part on real hardware.
Buoyancy Buoyancy is modeled as a single upward force on sections of the hull based on their volume and coefficients. Lift coefficients can be altered at run-time during the simulation, to experiment with the effects of changing buoyancy.
Wind We model wind as a turbulent flow field, using the Dryden turbulence model (MIL-F-8785C) [13], the implementation of which is from the FlightGear open source flight simulator [3]. It must be noted, that this frequency model simulates both the spacial and temporal variability of turbulence in a single, moving point only. As such, we can only calculate a single, time variable global flow vector, which we assume constant in space; the same assumption that was made by FlightGear. For small airships we deem this is an acceptable limitation. This flow vector is published as a ROS message, so it is available to other plugins at any point in time.
Lift and drag For all aerodynamically relevant primitive objects, including hull sections, we calculate the local lift and drag forces based on the orientation and motion of the object through the flow field, defined by the current wind  flow vector. These are then taken into account by the physics solver to model the effect on the whole vehicle as well as its deformations. We model two types of basic lifting primitives, "quasi-planar" and "quasi-cylindrical". Let f be the flow vector relative to the primitive, consisting of orthogonal components, f x from the front, f y from the left and f z from above. For quasi-planar types f y is ignored. Quasi-cylindrical objects are rotated around their x (forward) axis. This yield a rotated flow r f such that From that point on, they are treated identical to quasi-planar types. We calculate the angle of attack α, the normalized drag vectorf d , normalized lift vectorf l and dynamic pressure q d as Let A be the relevant area of the object, given coefficients c l0 , c d0 , c d1 and α stall . We calculate forces for lift F l and drag F d as Both functions are highly simplified linear approximations (Fig. 3) of the true lift and drag curves. For each primitive object, only 4 aerodynamic coefficients are needed.
The approximate drag coefficients for basic primitive shapes, both for frontal and sideways flow are commonly known, as are lift coefficients of thin rectangular plates or cylinders. Care has to be taken to adjust the coefficients when and set the coefficients for the hull sections accordingly, weighted by their expected drag contribution.
As another example, fins attached to the hull might encounter increased lift due to increased dynamic pressure of the airflow around the hull, which can be modeled by artificially inflated lift coefficients for these surfaces.
If an aerodynamic model exists for a specific vehicle, it should be possible to solve for a "best fit" of the primitive coefficients. However, our approach has the advantage that a roughly approximated aerodynamic model can be created based on geometric shape only, which we argue is a valuable property for evaluating new vehicle designs.
Non-rigidity Gazebo joint properties in compound objects can be manipulated through a ROS service while a simulation is running. We wrote a script that can simultaneously adjust both the rotational and linear spring stiffness, as well as the range of free rotation for all joint connections between the blimp hull and other components such as fins and gondola around all axis of rotation. The same program also adjusts the buoyancy of the blimp. This way, the effects of a drop or increase in lifting gas pressure and volume can be evaluated at any time, both regarding flight behavior and controllability (Fig. 4).

Control
We employ a cascaded PI controller design similar to the controller presented in [4,10], with extensions. Librepilot [14, 2] offers a vehicle-agnostic control hierarchy based on virtual control axes, which are mapped to physical actuators through a configurable matrix. This maps a "Pitch" virtual axis on the elevator servos, the "Yaw" virtual axis on both rudder servos and a yaw thruster, as well as a "Thrust" virtual axis on both main thrusters. We decided not to employ main thruster differential thrust for control, due to the low torsional moment and high strain inflicted between gondola and hull. A simple PI control loop controls the yaw rate, while a hierarchical controller with an outer loop proportional term controls the pitch angle based on an PI inner loop controlling pitching rate. Implementations for these already exist and are utilized both for fixed wing and multirotor control.
In extension, we implemented a novel path follower, which determines the setpoints for the lower level control: Pitch P , Yaw rateẎ , Thrust T and Thrust vector angle 0 ≤ γ ≤ π 2 . Its inputs are a setpoint velocity vector in 3d space v S , the airspeed v I , the velocity in 3d space v = [v x , v y , v z ] and the current attitude of the vehicle in Euler angles q = [q r , q p , q y ]. All the state estimates come from the flight controller's state estimation using an extended Kalman filter (EKF). If no airspeed sensor is available, v I is assumed to be the forward component of v. We compensated the wind to calculate motion relative to the flow field. Let q (v I ) be a function that rotates v I into the world-frame. We then calculate the estimated wind flow vector f l and the corrected velocity in airṽ as For each vehicle a minimum controllable airspeed v min and a maximum achievable airspeed v max are defined. We solve for a wind-corrected setpointṽ S such that If the wind is stronger than the vehicle's maximum speed, there might be no solution with b > 0. In this case, the solution with the highest b is chosen, which minimizes the positional error accumulated under these conditions. Separate controllers are employed for direction C d , airspeed C v , climb rate Cḣand thrust vector C γ . Let a, b be the signed angular difference between a and b while ±li shall denote an integral accumulator limited by ±l i . We then calculateẎ The main addition to [4,10] is the thrust vector term C γ which, in combination with the proportional climb speed cross term P ×T k p P in C v , employs the main thrusters for altitude control. For a vehicle without thrust vector control, the corresponding coefficients γ k p , P ×T k p would remain 0.

Middleware
The flight control algorithm is implemented in C/C++ and executed either on a physical embedded flight controller (Openpilot Revolution) [2] on the LTAV, both in flight or for HITL simulation, or compiled as a computer program for SITL simulation. The Librepilot Ground Control Software (GCS) connects to either via telemetry/networking and offers a HITL/SITL interface, in which the flight controller's sensor measurements are overridden with simulated data and actuator commands are sent to the simulator. We extended this interface, in order to connect with ROS, and as such our Gazebo simulation. We also wrote additional software that interfaces with the flight controller directly (via USB or Networking) and allows HITL/SITL without the GCS. Its main purpose however, is to send sensor and state estimate data to ROS components running on board and in flight. In turn ROS can send the flight controller guidance setpoints. This facilitates high level autonomous control, such as vision or perception based algorithms, through ROS. Of course, these components are useful beyond their application for airships.

Robotic Hardware
We equipped a commercial 5m blimp (Fig. 2) (CloudMedia Sopot, Sopot, Poland), designed for advertising and aerial photography, with an Openpilot Revolution flight controller (including an IMU, 3 axis magnetometer, barometer and processor), a GPS receiver, a digital airspeed sensor, and an onboard computer for data logging and wireless networking.

Experiments and Results
We conducted several experiments, both in simulation as well as real-world validation.

Simulation experiments
Simulation experiment 1 -manual simulated flight -Based on data from physical measurements and observations of a tethered real-world test flight 4 , the simulation model was configured and tested in manual simulated flight. Unknown coefficients were adjusted until simulation and reality were in qualitative agreement 5 .
Simulation experiment 2 -loitering -The control algorithm was implemented and control parameters were tuned in simulation, starting with the lowest level control loop. Given manual setpoint inputs, P and I rate control coefficients were determined for pitch and yaw rate, followed by P coefficient for the pitch control loop. When these were satisfactory, the C d control loop was tuned, followed by Cḣ, C v and C γ . We then followed an iterative tuning procedure, while the simulated blimp was constantly loitering around a position P 0 . Given current position estimate P , the velocity setpoint is calculated as With the simulated blimp loitering, the performance in each controlled domain was monitored and control coefficients for whichever parameter (speed, altitude, course) was deemed worst were adjusted until satisfactory results were achieved (Fig. 5).
Simulation experiment 3 -trajectory following -We employ the path planner in Librepilot to follow a waypoint sequence with 4 waypoints. We denotep = p |p| as the notation for a normalized vector. For any pair of consecutive waypoints P 0 , P 1 and desired velocity v 0 , we determine v S via In equation (19) the term v 0pb moves the vehicle along the desired path vector, while the second term corrects any orthogonal deviance from the path, as p b (p a ·p b ) projects p a on p b . This guidance scheme allows observations of the control behavior in steady state, while the vehicle follow a reproducible trajectory repeatedly (Fig. 6).
Simulation experiment 4 -wind and turbulence -We added wind and turbulence to situations from both simulation experiment 2 and experiment 3 and observe the effects. Wind velocities between v min and v max allowed the simulated blimp to hover on the spot with good accuracy. We noted altitude excursions of several meters due to vertical wind gust components. During the trajectory following, wind from the rear poses a challenge if v min is set too low, as the the blimp is almost stationary in relation to the air-stream. Consequently, the control surfaces have very limited authority. This led to a noticeable overshoot of waypoints due to the time it took for the blimp to reorient. This can be alleviated by commanding a higher v min (Fig. 7).
Simulation experiment 5 -deflation -We re-conducted both simulation experiment 2 and simulation experiment 3. During the experiment we adjusted the rigidity of the simulated blimp to simulate helium loss and observe the effects. Simulated loss of hull pressure and the resulting decrease in rigidity reduced the controllability of the blimp. A reduction in altitude accuracy and oscillations in both pitch, altitude and also in course were observed. The ability of the fins to flex and twist in response to control actuation not only adds additional oscillation modes, especially for the top fin which houses the IMU (both in simulation and on the real-world blimp), it also reduces the control effects, since the fins rotate in response to the control forces until a new force equilibrium is reached. This significantly reduces the overall control force enacted on the hull, since the angle of attack of the fin and attached rudder negate each other. Nevertheless, the blimp remains controllable, as long as sufficient buoyancy is present to maintain altitude (Fig. 8).

Real-world experiment
Loitering outdoor with wind and deflation -We configured the real-world blimp with the control coefficients determined in simulation and engaged position hold mode. Due to a previous ground test incident, the hull had developed significant leakage. As such we were able to reproduce the effects in simulation experiment 5 in real-world, with the blimp's internal pressure far below optimum and dropping. Although we were able to maintain buoyancy by removing ballast weights, the blimp's flight performance was below optimal. As predicted in simulation experiment 5, a high frequency oscillation developed in the yaw axis. The rudder actuation imparted momentum and flexion on the top fin which was picked up by the IMU gyroscope. The real-world blimp was flown with identical control coefficients to those used in simulation experiments. However, due to changes on the electronic speed controller (ESC), the blimp produced higher thrust than anticipated. This might have contributed to oscillations in pitch and altitude, which were similar, but more persistent and of higher frequency than those observed in simulation experiment 5. Our airspeed sensor failed to work and had to be deactivated, leading to v I estimated based on ground speed, but we expect minor impact due to the benign winds that day (< 1m/s average). The altitude and position were controlled reliably for more than 5 minutes of continuous flight time. One observed wind gust caused an altitude excursion, as predicted by simulation experiment 4. No intervention by the pilot was required at any time ( Fig. 9) 6 .

Discussion
Our experiments have confirmed, that changes in rigidity, as simulated in our framework, have significant effects on behavior and control of airships, which is well documented by the literature [11]. Modeling of wind gusts has also exposed notable control behavior, not seen in non-turbulent air. We have shown, that a simulation framework as presented here can be used to develop, tune and evaluate a control algorithm with results comparable to and representative of the real world, with limitations discussed below.

Limitations
Very large airships would benefit from wind and turbulence simulation, that accurately models local variations of the wind flow at the same point in time. This would require the wind model to be modified, to estimate a whole wind field and calculate the wind vector as a function of both space and time for multiple places. Due to the modular approach taken with aerodynamics, it would be relatively straightforward to incorporate this data into a simulation with a sufficiently segmented hull, although this is beyond the scope of our current work.
Our simulation, as conducted for our experiments, is not quantitatively accurate for our vehicle. To achieve quantitative accuracy, or to analyze the simulation error made by the model, it would be necessary to compare the calculated forces to those determined in wind tunnel tests of CFD analysis, which could be done with acceptable effort for the rigid case. This data could be used to solve for "best fitting" coefficients on all component objects, to minimize the simulation error. This would also show, if additional correction coefficients are needed to model aerodynamic interaction between different components, or if the current set of coefficients in combination with higher granularity is sufficient. We would defer this analysis to future work.  turning around after an overshoot leads to significant drift downwind. However, with the wind speed between vmin = 1m/s and vmax = 2m/s the blimp manages to hover in close proximity of the reference point. The "random walk" is an effect of turbulent fluctuations in wind speed and direction. When navigating a waypoint sequence, the ground speed is severely increased on the downwind leg, even at an airspeed of vmin, resulting in large turn radii and overshoot. On the upwind leg the ground speed is greatly reduced, coming to a near halt during strong gusts. The gusty wind poses challenges to the blimp's ability to maintain altitude, vertical gust components cause severe altitude excursions of several meters which the blimp then compensates to the best of its abilities. These are reproducible regardless of the guidance mode, leading to similar excursions at the same time in each experiment, whenever the same wind gusts are encountered. Simulated wind is coming from south-east (135°) with an average velocity of 1.5m/s, simulated turbulence for the Dryden model was set to magnitude 3. In the loiter experiment, this exposes a severe oscillation instability in altitude and pitch that was not present when fully inflated, although the blimp is still able to control itself. In both situations, the altitude accuracy is noticeably decreased. In this plot, the yaw gyroscope signal is plotted, revealing a high frequency oscillation in yaw, as the actuation of the rudder now moves the fin that the simulated gyroscope is mounted to, relative to the blimp. Our current aerodynamic model can not model "prop wash" and the effects it has on impinged control surfaces and fins. This should be addressed by future work to improve the blimp simulation, although solving this would also be beneficial for other vehicles such as multicopters and fixed wing.
Although our model can simulate deformations of the hull as a whole, it can not model surface deformations of the skin such as wrinkling or denting. Temporary changes of the geometric shape of a component could have significant effects on the drag coefficient as mentioned in [11]. These have been observed to occur in our real-world experiment. Although it would be straightforward to adjust the aerodynamic coefficients of component objects in the simulation, this would require knowledge of the relationship between aerodynamic forces, internal pressure, shape and effects on aerodynamics. This is a hard problem and has not been well researched yet, as noted in [11]. It justifies future work.  Figure 9: Real-world experiment -The blimp loiters around a given position in the physical world, while losing helium. The winds were benign, but a single gust at the 160s mark results in a noticeable drop in altitude which the control subsequently corrected, similar to those seen in simulation experiment 4. Due to lack of functional airspeed sensor and higher than tuned-for thrust, the ground speed and turn radius are larger than in the simulation. As predicted in simulation experiment 5, the blimp displays pitch-altitude oscillations throughout the flight. A high frequency oscillation in the yaw axis was also present (drawn down-scaled by factor 100), as predicted.