On a romantic date 12,000 kilometers from home, I sat in a hot air balloon wondering how the pilot managed to steer the aircraft. After all, it seemed like the balloon only had controls to go up and down, yet the flight path seemed so intentional. I wanted to understand the dynamics of these balloons, and build a virtual controller capable of maneuvering them.

Balloon Dynamics

The dynamics of hot air balloons has been well-studied. We use the derivation from Badgwell1 as the reference. Badgwell modeled the vertical dynamics of an AX7-772. The balloon that I rode in was a larger and heavier one, but we will build our simulation and controllers for an AX7-77 because it is useful to have a reference implementation to compare against.

We will summarize the dynamics at a high level and extend it to incorporate horizontal motion. Define the following terms. Most of the derivation is omitted for simplicity.

  • Let \(\textbf{p} = \left(p_\textrm{x}, p_\textrm{y}, p_\textrm{z}\right)\) be the position of the balloon.
  • Let \(\textbf{w} = \left(w_\textrm{x}, w_\textrm{y}, w_\textrm{z}\right)\) be the wind velocity at the current position of the balloon.
  • Let \(m_\textrm{t}\) be the total mass of the balloon, including the air in the envelope.
  • Let \(m_\textrm{p}\) be the payload mass of the balloon, excluding the air in the envelope.
  • Let \(k_\textrm{d}\) be the proportional constant of the drag force on the balloon.
  • Let \(F_\textrm{l}\) be the lift force generated by the hotter air in the envelope.
  • Let \(g\) be the acceleration due to gravity.

The acceleration in each component is given by the following equation. We assume the drag constant is the same in all directions, which treats the hot air balloon as a sphere. This is not a perfectly accurate model, but should be good enough given that the envelope represents most of the balloon’s surface area and is approximately spherical.

\[\frac{d^2\textbf{p}}{dt^2} = \frac{1}{m_\text{t}} \left[ k_\text{d} \left| \textbf{w} - \frac{d\textbf{p}}{dt} \right| \left[ \textbf{w} - \frac{d\textbf{p}}{dt} \right] + \begin{pmatrix} 0 \\ 0 \\ F_\text{l} - m_\text{p} g \end{pmatrix} \right]\]

The drag force is a function of the relative wind velocity. This is a generalization of Badgwell’s model where the wind velocity is zero everywhere. The lift force is a function of the height of the balloon and the temperature of air in the envelope. The height is necessary to determine the outside air temperature, which Badgwell models as linearly decreasing with height.

The air temperature inside the envelope is controlled by the fuel valve and the vent valve. Increasing the fuel valve position heats up the air in the envelope, while increasing the vent valve position lets hot air out from the top of the envelope, effectively cooling the air inside. There is also heat dissipation as a result of the cooler outside air. Therefore, the derivative of envelope air temperature is a function of the fuel and vent valve positions along with the temperatures of the inside and outside air.

Modelling Wind

To move horizontally, a hot air balloon relies on varying wind velocities at different altitudes. A uniform wind field is not very interesting to simulate since the balloon will always move in the same direction on the horizontal plane. However, a detailed model of localized wind patterns is beyond the scope of this post. Instead, we rely on a simplified wind model based on random control points.

Define the wind field as a function that takes in the position of the balloon and returns the wind velocity at that position. Given the boundaries of the simulation, we place evenly spaced control points. For each control point, we generate a random wind velocity within a specified magnitude. The output wind velocity at any point is determined through trilinear interpolation in the regular grid.

control_points = (
    np.linspace(-1000, 1000, 20),
    np.linspace(-1000, 1000, 20),
    np.linspace(0, 2000, 20),
)
control_vectors = (
    np.random.uniform(-5, 5, size=(20, 20, 20)),
    np.random.uniform(-5, 5, size=(20, 20, 20)),
    np.random.uniform(-1, 1, size=(20, 20, 20)),
)
wind_vector = (
    interpn(control_points, control_vectors[i], position)
    for i in range(3)
)

The above snippet shows an example of evaluating a random wind field defined in a 2 km × 2 km × 2 km grid. The horizontal wind magnitude is up to 5 m/s, and the vertical wind magnitude is up to 1 m/s. We show SciPy’s interpn3 function here, but the actual implementation relies on compiled Numba4 code since the wind field is evaluated in the hot path of the simulation.

The diagram above shows a 2D slice of the wind field from the previous example with an arbitrary seed. We can see that interpolation with control points does a fairly good job of making interesting wind fields while ensuring that there are no discontinuities.

Implementation

We are interested in simulating the hot air balloon with the aforementioned dynamics under varying wind fields and with different controllers. This sections goes into some of the implementation details. The full source code is available here.

Balloon

The core simulation loop for the balloon operates in fixed size time steps. In each step, we compute the derivative of the state vector consisting of position, velocity, and temperature. We then use odeint5 to compute the new state vector. Note that the wind field is evaluated in each call of the derivative function for more simulation accuracy, but it is also okay to assume that the wind field is constant for each step if the interval is small enough.

class Balloon:
    # ...

    def step(self, time_step: float):
        time_start = self.time
        time_end = time_start + time_step
        time_span = [time_start, time_end]

        x_start = np.empty(7, dtype=np.float64)
        x_start[0:3] = self.position
        x_start[3:6] = self.velocity
        x_start[6] = self.temperature
        x_end = odeint(self.derivative, x_start, time_span)[-1]

        self.position = x_end[0:3]
        self.velocity = x_end[3:6]
        self.temperature = x_end[6]

        if self.position[-1] <= 0.0:
            self.position = np.array([self.position.x, self.position.y, 0.0])
            self.velocity = np.array([0.0, 0.0, 0.0])

        self.time = time_end

A simplified step function is shown above. Like Badgwell’s implementation, we handle cases where the balloon is on the ground after a step by zeroing the velocity. This does not affect takeoff since the balloon should only ever be on the ground after a step if it is already on the ground with no vertical velocity or it is descending. In the derivative, we also prevent horizontal wind velocity from affecting the balloon if it is on the ground so that there is no horizontal movement unless the balloon is in the air.

Badgwell’s derivation uses a dimensionless model. However, we want to recover the dimensions in all of our simulation outputs. In the full implementation, the Balloon class interacts with the outside world in SI units, but stores its internal state in dimensionless values to simplify derivative calculations. The dimensioned wind field is passed in to a Balloon instance, and its outputs are scaled appropriately in the derivative.

Controller

The balloon is controlled by the fuel and vent valve positions. Each of them can be independently set to a percentage value and is assumed to be fixed in a given time step. We generalize the controller as a function that takes in the current observable state of the balloon and outputs the valve positions that should be applied before the next simulation step.

class FixedController:
    def __init__(self, output: ControllerOutput):
        self.output = output

    def __call__(self, input: ControllerInput) -> ControllerOutput:
        return self.output

class SequenceController:
    def __init__(self, *controllers: Tuple[float, Controller]):
        self.controllers = sorted(controllers, key=lambda t: t[0], reverse=True)
        self.last_controller = None

    def __call__(self, input: ControllerInput) -> ControllerOutput:
        while self.controllers and self.controllers[-1][0] <= input.time:
            self.last_controller = self.controllers.pop()[1]

        if self.last_controller is None:
            return ControllerOutput(fuel=0.0, vent=0.0)

        return self.last_controller(input)

We show two simple controllers that are useful for testing and tuning. The FixedController always returns the same output. The SequenceController takes in a sequence of controllers and switches between them based on the input time. More sophisticated controllers (described in a later section) can leverage other properties of the input state such as position and velocity.

Simulation

Given the Balloon and Controller, we can run the simulation using a simple loop that is shown below. For most of our experiments, we set the time step to one second and the total time to two hours. This is similar to the real flight time of a hot air balloon, allows enough time to observe the behavior of different controllers, and can be fully simulated in just a few seconds.

def simulate(
    balloon: Balloon,
    controller: Controller,
    time_step: float,
    total_time: float,
) -> Monitor:
    monitor = Monitor()
    monitor.update(balloon)

    start_time = balloon.get_time()
    num_steps = int(math.ceil((total_time - start_time) / time_step))

    for _ in range(num_steps):
        controller_input = get_controller_input(balloon)
        controller_output = controller(controller_input)
        apply_controller_output(balloon, controller_output)
        balloon.step(time_step)
        monitor.update(balloon)

    return monitor

The simulation updates and returns a Monitor instance. This is used to track the internal state of the balloon over time and has methods for data interpolation, plotting, and animation. As an example, we show the path of the reference simulation6 under a random wind field with horizontal magnitude up to 10 m/s in the diagram below.

Moving Vertically

Before we can take full advantage of the wind field, we first need to develop a controller that can move the balloon to a target vertical position. To be somewhat realistic, we put a 4 m/s cap7 on the maximum ascent and descent speeds even though the balloon can theoretically reach higher vertical speeds.

Controller

It turns out that PID controllers work surprisingly well for this use case even without any linearization. However, we cannot directly use a single PID controller with the setpoint as the desired vertical position because there would be no way to control the maximum velocity. Instead, we use a cascade of two PID controllers. The first controller targets a constant vertical position, and its output is used as the setpoint for the second controller that targets a constant vertical velocity.

We discretize the output of the vertical position PID controller to 0.1 m/s since more fine-grained control is likely not realistic. The vertical velocity PID controller needs to output the positions for both the fuel and vent valves. Although it is possible for both to be non-zero at the same time, we use a single output where negative values are vent and positive values are fuel (i.e., vent and fuel are mutually exclusive). We also discretize the vent valve positions to 1%.

Tuning

Both the vertical velocity and vertical position PID controllers need to be tuned in that order. Instead of manually trying gain values, we apply bayesian optimization8 to find the approximate best parameters. For both PID controllers, we use the negative mean absolute error as the objective function since it is simple and penalizes time away from the target sequence. More complex objectives can be used if we want to strictly prohibit overshooting at the cost of slower convergence.

controller = SequenceController(
    (0.0, VerticalPositionController(1000.0, k_p, k_i, k_d)),
    (1500.0, VerticalPositionController(500.0, k_p, k_i, k_d)),
    (3000.0, VerticalPositionController(750.0, k_p, k_i, k_d)),
    (4500.0, VerticalPositionController(250.0, k_p, k_i, k_d)),
)

The target sequence used to tune the vertical position PID controller is shown above. The construction is fairly arbitrary, but it allows the balloon to ascend and descend at varying altitudes. When using bayesian optimization, we consider P-only, I-only, D-only, PI, PD, and PID gain configurations since zero values on the boundaries are hard for the optimizer to explore by itself. We get fairly good results after only around 100 iterations.

Results

After tuning, we find that vertical velocity uses a full PID controller with a large derivative gain, while vertical position uses a P-only controller. This makes sense because the balloon is a high-inertia system. Envelope temperature changes slowly in response to fuel and vent valve positions, and the balloon also has a lot of momentum due to its mass.

We can see from the chart above that the tuned vertical position controller is able to effectively track the target sequence while ensuring that velocity does not exceed 4 m/s when ascending or descending. This is not a realistic depiction of how a pilot would control the balloon since we allow the fuel and vent valve positions to change once per second. We could increase the sample time of the PID controllers so that the outputs change slower, but the existing configuration is sufficient for our purposes.

Moving Horizontally

With a controller that can target arbitrary vertical positions, we can use it move the balloon horizontally by taking advantage of differing wind velocities at different heights. The general approach for all of the horizontal controllers is the same. Given the current state of the balloon and the wind field, output a desired height that the balloon should go to. This is then passed into the vertical position controller to adjust the fuel and vent valve positions accordingly.

Fixed Controller

This is a simple controller that is just used for establishing a baseline. It holds the balloon at a fixed height of 500 m. We expect all other controllers to perform better than this.

Greedy Controller

This controller greedily finds the best height that the balloon should go to at every step. The idea is that moving to a height where the wind vector points in the direction of the target will get us closer to the target. To bound computation time, we discretize the input space into grids of 100 m. The wind field is sampled at the center of each grid, and moving to a grid targets the height at the center of the grid. A rough outline of the controller logic is below.

  1. If the balloon is in the same horizontal grid as the target, we have effectively reached the target. Output the target height.
  2. Compute the normalized vector from the balloon’s current position to the target, ignoring the vertical component. Call this \(\mathbf{v_t}\).
  3. For each grid above and below the balloon’s current grid, compute the normalized vector for the wind velocity at that grid, ignoring the vertical component. Call this \(\mathbf{v_w}\).
  4. Compute the cosine similarity between \(\mathbf{v_t}\) and \(\mathbf{v_w}\) at each grid from the previous step. The grid with the highest similarity is the best grid to move to. Output the height at the center of that grid.

Search Controller

This controller performs path finding in 3D space. The idea is that we can output the sequence of heights along the shortest path from the current position to the target. Similar to the greedy controller, we discretize the input space into grids of 100 m. These grids are then used as vertices in the search graph. Edges are determined based on the direction of the wind at the center of each grid.

It is not guaranteed that the balloon will actually follow the exact search path because edges are very approximate. Therefore, we want the controller to be able to efficiently switch to a new path if the balloon deviates. In addition, the target may not be reachable at all given the wind field. However, the search should still be able to find the reachable grid which is closest to the target grid. Given these considerations, we can construct the search graph in the following way.

  1. Add all grids as vertices. Call these vertices \(v_{i,j,k}\) where \(i\), \(j\), and \(k\) are the indices of the grid.
  2. Add an unreachable vertex \(v_u\). This vertex will only show up in the search path if the target grid is not reachable from the current grid.
  3. Define the weight of an edge as \(\left( w_u, w_r \right)\). \(w_u\) is the cost of moving to the unreachable vertex. \(w_r\) is the cost of moving to any other vertex. Comparison on weights should treat any non-zero \(w_u\) as larger than any non-zero \(w_r\). This is already how tuple comparison works in Python.
  4. For each \(v_{i,j,k}\), add an edge to one of its four horizontally adjacent neighbors \(v_{i \pm 1,j \pm 1,k}\) by looking at the dominant horizontal component of the wind vector at grid \(\left( i,j,k \right)\). The weight of each edge is \(\left( 0, c \right)\) where \(c\) is the time it takes for the balloon to move between the two grids taking into account the projected wind velocity.
  5. For each \(v_{i,j,k}\), add edges to \(v_{i,j,k+1}\) and \(v_{i,j,k-1}\). The balloon can always move vertically. The weight of these edges are \(\left( 0, c \right)\) where \(c\) is the time it takes for the balloon to move between the two grids assuming the maximum vertical speed of 4 m/s.
  6. For each \(v_{i,j,k}\), add an edge to \(v_u\). The weight of each edge is \(\left( c, 0 \right)\) where \(c\) is the distance from grid \(\left( i,j,k \right)\) to the target grid.
  7. Add an edge from \(v_u\) to the target grid with weight \(\left( 0, 0 \right)\).

Instead of searching in the forward graph, we run Dijkstra’s algorithm on the reverse graph starting from the target grid. When the search terminates, we will have the shortest path from every grid to the target grid. Even if the balloon deviates from the chosen path, we do not need to rerun the search. This allows us to implement the controller as follows.

  1. Construct the graph once on controller initialization. Run Dijkstra’s algorithm on the reverse graph starting from the target grid. Store the output in a parent map.
  2. On controller invocation, if the balloon is in the same grid as the previous step, do nothing because the shortest path cannot change unless the balloon moves to a different grid.
  3. Look up the current grid in the parent map. Because we reversed the graph, the parent of the current grid will be the next grid to move to. If the next grid exists and is not the unreachable grid, output the height at the center of the next grid. Otherwise, continue holding at the current height because we have gotten as close as we can to the target.

Results

To compare the effectiveness of different controllers, we evaluate how close their trajectories get in the horizontal plane to a random target point 2 km away from the starting point under a random wind field with horizontal magnitude up to 10 m/s. Each controller is simulated 100 times with the same set of seeds. Note that we do not care about vertical distance to the target point since the vertical position is fully controllable. The results are shown below. Note that the last bucket is \([1900, \infty)\).

As expected, the search controller significantly outperforms the greedy controller and the fixed controller in most cases. There are instances where the greedy controller and search controller encounter wind fields that make it impossible to move towards the target. In around half of these simulations, the search controller is able to get within 100 m of the target. This is reinforced in the summary statistics below.

Controller Mean Median Standard Deviation
Fixed 1688 1796 345
Greedy 729 473 707
Search 424 162 542

For a given seed, we show the trajectory of all three controllers below. We can see that the search controller gets to within a few meters of the target through a fairly complex path. The greedy controller does eventually get somewhat close to the target, but nowhere near as close as the search controller. The fixed controller eventually gets stuck in a dead zone with little to no wind velocity.

Extensions

There are a few extensions that might be interesting to explore. First, we assume near-perfect knowledge of the wind field in our controllers, but this is not realistic. Pilots do plot out their courses using some wind data, but they do not know the exact horizontal velocity of the balloon until they are in the air. It would be interesting to develop a controller that can perform path finding with uncertainty in the search space. Rather than getting closest to a chosen point, we may want the controller to find the most likely points that are reachable. An application of this is to plan the best route for the balloon.

When I was on the hot air balloon flight, I heard our pilot communicating wind data with other pilots in the same fleet. This led me to wonder if there was some way for multiple balloons to collaborate and find paths more precisely as a group by sharing wind data. In addition, the pilots had to make sure that the balloons did not collide with each other during landing and takeoff. It would be interesting to further explore the idea of multi-agent path finding when applied to hot air balloons.

References

  1. Badgwell, Thomas (2017). Dynamic Simulation of a Hot Air Balloon. ↩

  2. Head Balloons (2024). AX7-77. ↩

  3. The SciPy Community (2024). interpn - SciPy Manual. ↩

  4. Anaconda (2024). Numba: A High Performance Python Compiler. ↩

  5. The SciPy Community (2024). odeint - SciPy Manual. ↩

  6. Badgwell, Thomas (2017). Hot Air Balloon Simulation and Control. ↩

  7. Kämpf, Peter (2014). How quickly can a hot air balloon climb?. ↩

  8. Nogueira, Fernando (2014). Bayesian Optimization: Open source constrained global optimization tool for Python. ↩