Hot Air Balloon Dynamic Control - Optimization for Decision Science

Prepared by: Stephen Cini (scini@nd.edu, @sscini, 2024) and Oliver Harris (oharris2@nd.edu, @OliverHarris47, 2024)

# Install Pyomo and solvers for Google Colab
import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    !pip install okabeito --quiet
    import helper

    helper.easy_install()

else:
    # Make sure helper.py is installed locally in same location as notebook
    sys.path.insert(0, "../")
    %pip install okabeito --quiet
    import helper

helper.set_plotting_style()

import pyomo.environ as pyo
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pyomo.dae as dae
import idaes
import okabeito
from pyomo.environ import units
from pyomo.util.check_units import assert_units_consistent

Note: you may need to restart the kernel to use updated packages.

Learning Objectives¶

Review syntax and components in Pyomo.DAE
See how change in objective function can affect model behavior
Highlight importance of realistic initialization conditions and parameter values
Carry physical units through a dynamic model, and check them at the right moment (see Units in Pyomo.dae Models)

Introduction¶

In this notebook, we adapt the example in Pyomo.DAE Example: Temperature Control Lab to discuss the model formulation for controlling height in a hot air balloon during a planned hot air balloon ride.

Picture adapted from https://www.real-world-physics-problems.com/hot-air-balloon-physics.html.

Instead of controlling the temperature of the heater and sensor utilizing a heater, we control the height of the hot air balloon by utilizing the lift force generated by the balloon being heated.

Hot Air Balloons¶

A hot air balloon operates based on the principle of buoyancy, where hot air inside the balloon is less dense than the cooler air outside, causing the balloon to rise. The principal forces for hot air balloon operation are: lift, weight, and drag. The objective of our project is to design a Pyomo model that controls a hot air balloon to follow the average path taken by a balloon during a tourism trip. A free body diagram describing how the principal forces act can be seen below.

Free body diagram of a hot air balloon: lift up, weight down, drag opposing motion

Picture adapted from https://charm.stanford.edu/ENGR1052016/AnthonyBombikJayantMukhopadhaya.

Model Formulation: No Set Objective¶

\begin{align*} \max ~~ & 0 \\ \mathrm{s.t.}~~ & \frac{dh}{dt} = v(t) \\ & D(t) = b \cdot v(t)^2 \\ & \frac{dv}{dt} \cdot m = L(t) - m \cdot g - D(t) \\ & h(t_0) = 0 \\ & v(t_0) = 0 \\ & L(t) = L_{ramp} = 10000, \quad \forall t \in [t_0, t_f] \end{align*}

(1)

The equation of motion for the hot air balloon is then:

\begin{align*} & m \cdot \ddot{y} = L - m \cdot g - D \\ & D = b \dot{y}^2 \end{align*}

(2)

Where:

\begin{align*} m & \text{ : Mass (kg)} \\ \ddot{y} & \text{ : Acceleration (m/s}^2\text{)}\\ \dot{y} & \text{ : Velocity (m/s)}\\ L & \text{ : Lift Force (N or kg m/s}^2\text{)} \\ g & \text{ : Gravitational Constant (m/s}^2\text{)} \\ D & \text{ : Drag Force (N)} \\ b & \text{ : Drag Constant (kg/m)} \end{align*}

(3)

Degree of Freedom Analysis¶

Versions 0 --> 1

In our model formulations for versions 0 and 1, we are simulating the balloon dynamics using a step-input (Lift). For plotting purposes, lift was defined as a variable with respect to time. However, a constraint was added that fixed the lift force to 10,000 Newtons. This effectively treats lift as a parameter and not a variable.

Variables: h(t) , v(t) , D(t)
Number of Constraints: 3

$DOF = Variables - Constraints = 3 - 3 = 0$

For versions 0 and 1 of our model, the models have 0 DOF. This matches the 'ipopt' solver output:

Total Number of Variables: 184

Total Number of Equality Constraints: 184

Define Model Constants¶

Here, we define the characteristics of the balloon we will be simulating in this notebook.

# Define Fixed Parameters from model inspiration sources.
#
# Attaching Pyomo units here lets `assert_units_consistent` check every equation
# of motion below, and it makes the derived quantities S and b self-documenting:
# their units are computed, not asserted in a comment.
rho_air = 1.225 * units.kg / units.m**3
g = 9.81 * units.m / units.s**2
radius = 10 * units.m
m_balloon = 750 * units.kg
cD = 0.47  # drag coefficient, dimensionless

S = 4 * np.pi * radius**2  # wetted surface area, m^2
b = 1 / 2 * rho_air * cD * S  # drag constant, kg/m

# Print calculated values. pyo.value() strips the units off so that numpy can
# round the number -- np.round does not know what to do with a Pyomo expression.
print(f"S = {np.round(pyo.value(S), 2)} {units.get_units(S)}")
print(f"b = {np.round(pyo.value(b), 2)} {units.get_units(b)}")

S = 1256.64 m**2
b = 361.75 kg/m

Function Definition¶

Here, we set up the Discretization/Initialization, Solve, and Plotting/Visualization Functions for use in each of our models for different scenarios with the hot air balloon. Similar to the TCLab, we use dae.collocation for discretizing the model, and use the continuous solver 'ipopt'.

# Set-Up the Discretization/Initialization Function for the Model
def discretize_initialize(m, nfe):
    """Discretize and initialize model.

    Inputs:
        m: Model to discretize.
        nfe: Number of finite elements

    Returns:
        m: Discretized model.

    """

    # Check dimensional consistency BEFORE discretizing.
    #
    # This is the last moment at which the check is meaningful. The collocation
    # equations added on the next line multiply state variables by numerical
    # reciprocal-delta-t coefficients that carry no units, because a
    # ContinuousSet cannot hold any -- Pyomo issue #1790, still open:
    #   https://github.com/Pyomo/pyomo/issues/1790
    # Asserting after `apply_to` raises InconsistentUnitsError on a model that is
    # perfectly correct. See ./units_and_pyomo_dae.md
    assert_units_consistent(m)

    # Transform the Model using the 'dae.collocation' strategy.
    discretizer = pyo.TransformationFactory("dae.collocation")
    discretizer.apply_to(m, nfe=nfe, wrt=m.t)  # default ncp = 3

    # Specify Initial Conditions
    m.h[0].fix(0)
    m.v[0].fix(0)

    return m


# Set-Up Solve Function
def solve_model(m, tee=True):
    """Solves model with ipopt.

    Inputs:
        m: Model to solve.
        tee: Whether to print full output

    Returns:
        m: Solved model.

    """

    # Specify Solver
    solver = pyo.SolverFactory("ipopt")
    # Solve the Model
    results = solver.solve(m, tee=tee)
    assert pyo.check_optimal_termination(results), (
        f"Solve failed: status={results.solver.status}, "
        f"termination={results.solver.termination_condition}"
    )

    return m


# Set up list for colors using okabeito color scheme
okabeito_list = [
    okabeito.lightblue,
    okabeito.yellow,
    okabeito.orange,
    okabeito.green,
    okabeito.purple,
    okabeito.red,
    okabeito.blue,
    okabeito.black,
]  # list of colors from the okabeito color pallete

# Linestyles cycled alongside the colours, so a sweep of several series still
# reads when the page is printed in black and white (figures/README.md).
linestyle_list = [
    "-",
    "--",
    "-.",
    ":",
    (0, (3, 1, 1, 1)),
    (0, (5, 1)),
    (0, (1, 1)),
    (0, (3, 1, 1, 1, 1, 1)),
]


# Set-Up the Plotting/Visualization Function
def visualize_model(m):
    """Plotting model function with publication quality for main plots
    used in notebook.

    Inputs:
        m: Model to plot

    Returns:
        fig, ax: Plot information from model. Allows for further additions."""
    # Collect the data
    t = [i for i in m.t]
    h = [pyo.value(m.h[i]) for i in m.t]
    v = [pyo.value(m.v[i]) for i in m.t]
    D = [pyo.value(m.D[i]) for i in m.t]
    L = [pyo.value(m.L[i]) / 1000 for i in m.t]  # plotting in kN

    # Plot with subplots
    fig, axs = plt.subplots(3, 1, dpi=200, figsize=(8, 12.8))

    colors = okabeito_list
    # Plot Height vs Time
    axs[0].plot(t, h, label="height_model", linewidth=3, color=colors[0])
    axs[0].set_ylabel("Height (m)", fontweight="bold", fontsize=16)

    # Plot Velocity vs Time
    axs[1].plot(t, v, linewidth=3, color=colors[0])
    axs[1].set_ylabel("Velocity (m/s)", fontweight="bold", fontsize=16)

    # Lift (m/s^2) vs time
    axs[2].plot(t, L, linewidth=3, color=colors[0])
    axs[2].set_xlabel("Time (s)", fontweight="bold", fontsize=16)
    axs[2].set_ylabel("Lift (kN)", fontweight="bold", fontsize=16)

    # Adjust layout to prevent overlap and add spacing
    axs[0].tick_params(direction="in", labelsize=15)
    axs[1].tick_params(direction="in", labelsize=15)
    axs[2].tick_params(direction="in", labelsize=15)
    fig.tight_layout()
    return fig, axs

Version 0: Simulate the Balloon Dynamics¶

We modeled our hot air balloon with an adapted model from CHARM Lab at Stanford University (https://charm.stanford.edu/ENGR1052016/AnthonyBombikJayantMukhopadhaya). Their analysis was initially done in MATLAB with Simulink to define a controller for the hot air balloon. Our first step was to transfer the model to python and model with zero degrees of freedom to see if we get the expected output.

# Set up our model as a function
def create_balloon_model_v0():
    """Model function. All other models are modified versions.
    Inputs: None
    Returns:
        m: Base model to be discretized, initialized, solved, and plotted."""

    # Create a Pyomo Model
    m = pyo.ConcreteModel("Hot Air Balloon")

    # Define Time Domain. NOTE: a ContinuousSet cannot carry units (Pyomo issue
    # #1790), so the seconds live in this comment. Bounds are in SECONDS.
    m.t = dae.ContinuousSet(bounds=(0, 100))

    # Define Fixed Control Input
    m.L = pyo.Var(m.t, initialize=10000, units=units.N)  # Newtons

    # Define State Variables as Functions of Time
    m.h = pyo.Var(m.t, bounds=(0, None), units=units.m)  # height
    m.v = pyo.Var(m.t, units=units.m / units.s)  # vertical velocity
    m.D = pyo.Var(m.t, units=units.N)  # drag force

    # DerivativeVar defaults to DIMENSIONLESS and does NOT inherit units from the
    # state variable. m.t is time in seconds, so state it: dh/dt is m/s and
    # dv/dt is m/s^2.
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t, units=units.m / units.s)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t, units=units.m / units.s**2)

    # Define Differential Constraints
    def dh_con(m, t):
        return m.dhdt[t] == m.v[t]

    m.dh_con = pyo.Constraint(m.t, rule=dh_con)

    def dv_con(m, t):
        return m.dvdt[t] * m_balloon == m.L[t] - m_balloon * g - m.D[t]

    m.dv_con = pyo.Constraint(m.t, rule=dv_con)

    # Define Algebraic Constraints
    def drag_con(m, t):
        return m.D[t] == b * m.v[t] ** 2

    m.drag_con = pyo.Constraint(m.t, rule=drag_con)

    # Define L
    def l_con(m, t):
        return m.L[t] == 10000 * units.N

    m.l_con = pyo.Constraint(m.t, rule=l_con)

    return m

Once each model is defined, we discretize, initialize, solve, and plot the results.

# Create a Model Instance
model = create_balloon_model_v0()

# Specify Number of Finite Elements
NFE = 10

# Call the Discretization/Initialization Function
discretize_initialize(model, NFE)

# Call the Solve Function
model = solve_model(model)

# Plot the Results
visualize_model(model)

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:      540
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       30

Total number of variables............................:      184
                     variables with only lower bounds:       30
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      184
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 2.64e+03 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 2.81e+06 2.87e+10  -1.0 1.76e+04    -  5.62e-07 2.50e-01f  3
   2  0.0000000e+00 7.01e+05 7.18e+09  -1.0 2.74e+03    -  2.87e-01 1.00e+00h  1
   3  0.0000000e+00 1.74e+05 1.79e+09  -1.0 1.08e+03    -  1.00e+00 1.00e+00h  1
   4  0.0000000e+00 4.30e+04 4.43e+08  -1.0 5.22e+02    -  1.00e+00 1.00e+00h  1
   5  0.0000000e+00 1.01e+04 1.03e+08  -1.0 2.36e+02    -  1.00e+00 1.00e+00h  1
   6  0.0000000e+00 2.00e+03 1.73e+07  -1.0 8.86e+01    -  1.00e+00 1.00e+00h  1
   7  0.0000000e+00 2.14e+02 1.23e+06  -1.0 2.11e+01    -  1.00e+00 1.00e+00h  1
   8  0.0000000e+00 3.89e+00 1.19e+04  -1.0 3.12e+00    -  1.00e+00 1.00e+00h  1
   9  0.0000000e+00 1.35e-03 2.26e+00  -1.0 1.17e-01    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  0.0000000e+00 1.52e-10 1.00e-06  -1.0 7.82e-05    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   1.5234036254696548e-10    1.5234036254696548e-10
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.5234036254696548e-10    1.5234036254696548e-10


Number of objective function evaluations             = 12
Number of objective gradient evaluations             = 11
Number of equality constraint evaluations            = 15
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 11
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 10
Total CPU secs in IPOPT (w/o function evaluations)   =      0.002
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

(<Figure size 1600x2560 with 3 Axes>,
 array([<Axes: ylabel='Height (m)'>, <Axes: ylabel='Velocity (m/s)'>,
        <Axes: xlabel='Time (s)', ylabel='Lift (kN)'>], dtype=object))

Version 1: Correcting the Drag Force Calculation¶

After creating the base model, we noticed a potential problem in the Drag Force calculation. According to the drag force calculation using $v^2$ , the value of drag will always be subtracted from the other forces.

After consulting with Dr. Dowling, drag calculation was adjusted to what is shown in the next iteration of the model, using $D = b|v|v$ . Absolute value is not differentiable, having a kink at 0, so to model in Pyomo, we adjusted to make the function smooth by making $|v| = \sqrt{v^2 + \epsilon^2}$ to take away the singularity at zero.

\begin{align*} \max ~~ & 0 \\ \mathrm{s.t.}~~ & \frac{dh}{dt} = v(t) \\ & D(t) = b \cdot v(t) \sqrt{v(t)^2 + \epsilon^2} \\ & \frac{dv}{dt} \cdot m = L(t) - m \cdot g - D(t) \\ & h(t_0) = 0 \\ & v(t_0) = 0 \\ & L(t) = L_{ramp} = 10000, \quad \forall t \in [t_0, t_f] \end{align*}

(4)

# Set-Up our Model as Function
def create_balloon_model_v1():
    """Modified from v0. Fixed drag."""

    # A velocity-scale epsilon that removes the |v| singularity at v = 0.
    # It is compared against v, so it must carry v's units.
    eps = 1e-6 * units.m / units.s

    # Create a Pyomo Model
    m = pyo.ConcreteModel("Hot Air Balloon")

    # Define Time Domain. NOTE: a ContinuousSet cannot carry units (Pyomo issue
    # #1790), so the seconds live in this comment. Bounds are in SECONDS.
    m.t = dae.ContinuousSet(bounds=(0, 100))

    # Define State Variables as Functions of Time
    m.h = pyo.Var(m.t, bounds=(0, None), units=units.m)  # height
    m.v = pyo.Var(m.t, units=units.m / units.s)  # vertical velocity
    m.D = pyo.Var(m.t, units=units.N)  # drag force

    # Define Fixed Control Input
    m.L = pyo.Var(m.t, initialize=500, units=units.N)

    # DerivativeVar defaults to DIMENSIONLESS and does NOT inherit units from the
    # state variable. m.t is time in seconds, so state it: dh/dt is m/s and
    # dv/dt is m/s^2.
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t, units=units.m / units.s)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t, units=units.m / units.s**2)

    # Reformulate m.D to Add Stability
    def drag_con(m, t):
        return m.D[t] == b * (m.v[t] ** 2 + eps**2) ** (1 / 2) * m.v[t]

    m.drag_con = pyo.Constraint(m.t, rule=drag_con)

    # Define Differential Constraints
    def dh_con(m, t):
        return m.dhdt[t] == m.v[t]

    m.dh_con = pyo.Constraint(m.t, rule=dh_con)

    def dv_con(m, t):
        return m.dvdt[t] * m_balloon == m.L[t] - m_balloon * g - m.D[t]

    m.dv_con = pyo.Constraint(m.t, rule=dv_con)

    # Define L

    def l_con(m, t):
        return m.L[t] == 10000 * units.N

    m.l_con = pyo.Constraint(m.t, rule=l_con)

    return m

# Create a Model Instance
model = create_balloon_model_v1()

# Call the Discretization/Initialization Function
discretize_initialize(model, NFE)

# Call the Solve Function
model = solve_model(model)

# Visualize solution
fig, axs = visualize_model(model)

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:      540
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       30

Total number of variables............................:      184
                     variables with only lower bounds:       30
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      184
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 9.50e+03 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 4.49e+07 4.59e+11  -1.0 1.76e+04    -  5.62e-07 1.00e+00f  1
   2  0.0000000e+00 1.12e+07 1.15e+11  -1.0 8.81e+03    -  1.00e+00 1.00e+00h  1
   3  0.0000000e+00 2.81e+06 2.87e+10  -1.0 4.40e+03    -  1.00e+00 1.00e+00h  1
   4  0.0000000e+00 7.01e+05 7.18e+09  -1.0 2.19e+03    -  1.00e+00 1.00e+00h  1
   5  0.0000000e+00 1.75e+05 1.79e+09  -1.0 1.08e+03    -  1.00e+00 1.00e+00h  1
   6  0.0000000e+00 4.30e+04 4.43e+08  -1.0 5.23e+02    -  1.00e+00 1.00e+00h  1
   7  0.0000000e+00 1.01e+04 1.03e+08  -1.0 2.36e+02    -  1.00e+00 1.00e+00h  1
   8  0.0000000e+00 2.00e+03 1.73e+07  -1.0 8.87e+01    -  1.00e+00 1.00e+00h  1
   9  0.0000000e+00 2.14e+02 1.23e+06  -1.0 2.11e+01    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  0.0000000e+00 3.90e+00 1.19e+04  -1.0 3.12e+00    -  1.00e+00 1.00e+00h  1
  11  0.0000000e+00 1.35e-03 2.26e+00  -1.0 1.17e-01    -  1.00e+00 1.00e+00h  1
  12  0.0000000e+00 1.54e-10 1.00e-06  -1.0 7.84e-05    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 12

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   1.5415935195051134e-10    1.5415935195051134e-10
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.5415935195051134e-10    1.5415935195051134e-10


Number of objective function evaluations             = 13
Number of objective gradient evaluations             = 13
Number of equality constraint evaluations            = 13
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 13
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 12
Total CPU secs in IPOPT (w/o function evaluations)   =      0.003
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

# Uncomment to view model printout after solving
# model.pprint()

Analysis of Model Results: Varying Fixed Lift Value¶

How does the set lift affect the velocity?¶

# Set-Up our Model as Function
def create_balloon_model_v1b():
    """Modified from v0. Fixed drag, made m.L a mutable parameter."""

    # A velocity-scale epsilon that removes the |v| singularity at v = 0.
    # It is compared against v, so it must carry v's units.
    eps = 1e-6 * units.m / units.s

    # Create a Pyomo Model
    m = pyo.ConcreteModel("Hot Air Balloon")

    # Define Time Domain. NOTE: a ContinuousSet cannot carry units (Pyomo issue
    # #1790), so the seconds live in this comment. Bounds are in SECONDS.
    m.t = dae.ContinuousSet(bounds=(0, 100))

    # Define State Variables as Functions of Time
    m.h = pyo.Var(m.t, bounds=(0, None), units=units.m)  # height
    m.v = pyo.Var(m.t, units=units.m / units.s)  # vertical velocity
    m.D = pyo.Var(m.t, units=units.N)  # drag force

    # Define Fixed Control Input
    m.L = pyo.Param(initialize=10000, mutable=True, units=units.N)

    # DerivativeVar defaults to DIMENSIONLESS and does NOT inherit units from the
    # state variable. m.t is time in seconds, so state it: dh/dt is m/s and
    # dv/dt is m/s^2.
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t, units=units.m / units.s)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t, units=units.m / units.s**2)

    # Reformulate m.D to Add Stability

    def drag_con(m, t):
        return m.D[t] == b * (m.v[t] ** 2 + eps**2) ** (1 / 2) * m.v[t]

    m.drag_con = pyo.Constraint(m.t, rule=drag_con)

    # Define Differential Constraints
    def dh_con(m, t):
        return m.dhdt[t] == m.v[t]

    m.dh_con = pyo.Constraint(m.t, rule=dh_con)

    def dv_con(m, t):
        return m.dvdt[t] * m_balloon == m.L - m_balloon * g - m.D[t]

    m.dv_con = pyo.Constraint(m.t, rule=dv_con)

    return m

# Create a Model Instance
m = create_balloon_model_v1b()
nfe = 20

# Check units before discretizing -- see discretize_initialize() above and
# ./units_and_pyomo_dae.md for why the order matters.
assert_units_consistent(m)

# Discretize
# Transform the Model using the 'dae.collocation' strategy.
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(m, nfe=nfe, wrt=m.t)

# Set initial conditions
m.h[0].fix(100)  # Setting initial height above 0 to allow negative movement
m.v[0].fix(0)

fig, axs = plt.subplots(2, 1, dpi=300, figsize=(8, 12.8))
colors = okabeito_list

lift_range = np.linspace(7, 14, 8)
count = 0
for i in lift_range:
    # print(i) Uncomment to print lift_range[i] and track progress
    m.L = i * 1000  # newtons; m.L was declared with units=units.N
    m = solve_model(m, tee=False)
    t = [j for j in m.t]
    h = [m.h[j]() for j in m.t]
    v = [m.v[j]() for j in m.t]
    ls = linestyle_list[count % len(linestyle_list)]
    axs[0].plot(
        t, h, linewidth=3, label=f"L = {i} kN", color=colors[count], linestyle=ls
    )
    axs[1].plot(t, v, linewidth=3, color=colors[count], linestyle=ls)
    plt.xlabel("Time (s)", fontweight="bold", fontsize=16)
    axs[0].set_ylabel("Height (m)", fontweight="bold", fontsize=16)
    axs[1].set_ylabel("Velocity (m/s)", fontweight="bold", fontsize=16)
    count += 1

axs[0].legend()
axs[0].tick_params(direction="in", labelsize=15)
axs[1].tick_params(direction="in", labelsize=15)
plt.tight_layout()
plt.show()

Do these results make sense?

Let’s do a quick calculation to check our work. The lift force needs to be at least equal to the mass of the balloon, which is set to 750 kg. Therefore, the value needs to be at least (750 kg) * (9.81 m/s^2), or ~ 7360 N. We decided to use a value of 10,000 N to see the balloon reach a constant velocity.

To estimate the result of the solver, we calculate the velocity achieved by setting the Equation 1 = 0 for the force balance on the balloon. We then derived the velocity.

\begin{align*} & 0 = L - mg - bv^2 \\ & v = (\frac{L - mg}{b})^{0.5} \\ \end{align*}

(5)

So with the known values of L, mg, and b, the velocity should be 2.7 m/s. This matches the resulting value from the solved model (green line).

Model Formulation - Minimize ISE Objective¶

The next phase of model development involves adding a reference height or path and having the model objective be to minimize the integral of squared error between the height at each model step and the reference.

\begin{align*} \min_{L(t)} ~~ & \int_{t_0}^{t_f} \|~h(t) - r(t)~\|^2\,dt \\ \mathrm{s.t.}~~ & \frac{dh}{dt} = v(t) \\ & D(t) = b \cdot v(t) \sqrt{v(t)^2 + \epsilon^2} \\ & \frac{dv}{dt} \cdot m = L(t) - m \cdot g - D(t) \\ & h(t_0) = 0 \\ & v(t_0) = 0 \\ & L(t_0) = 10000 \end{align*}

(6)

Degree of Freedom Analysis¶

Versions 2 --> 4

In model formulations for versions 2, 3, and 4, the lift force is now treated as a variable with respect to time and is no longer fixed. Lift became the control variable and added a degree of freedom to our system.

Variables: h(t) , v(t) , D(t) , L(t) Number of Constraints: 3

$DOF = 4-3 = 1$

Versions 2 through 4 of our model have 1 DOF allowing the optimizer to freely manipulate our decision variable.

Version 2: Reach a Desired Final Height¶

\begin{align*} \min_{L(t)} ~~ & \|~h(t_f) - h_{desired}~\|^2 \\ \mathrm{s.t.}~~ & \frac{dh}{dt} = v(t) \\ & D(t) = b \cdot v(t) \sqrt{v(t)^2 + \epsilon^2} \\ & \frac{dv}{dt} \cdot m = L(t) - m \cdot g - D(t) \\ & h(t_0) = 0 \\ & v(t_0) = 0 \\ & L(t_0) = 10000 \end{align*}

(7)

In the next iteration of the model, we add an objective to minimize the sum of squared error (SSE) between a target height and the height reached at the end of the model.

This is the first iteration containing the control variable being utilized as a variable in the model. Therefore, we needed to define reasonable bounds for the lift of our balloon.

The lift force, as stated in Engineering Toolbox, is defined by the following equation: $L[t] = V(ρ_{out} - ρ_{in}) \times g$

Volume: We defined parameters regarding the balloon earlier in the notebook. In the referenced CHARM Lab balloon model, they chose to make a hemispherical balloon with a radius of 10 m. For this model, we made a fully spherical balloon of the same radius. So the volume of our balloon is $4/3*π*R^3$ , or roughly 4190 $m^3$ .

Density: The density of air at 100°C (373.15 K) is approximately 0.946 $kg/m^3$ , and at 20 °C is 1.205 $kg/m^3$ .

Gravity: We use the approximation of gravity $9.81 m/s^2$ in all relevant notebook calculations.

For the maximum value of lift, we use the lift force calculation with our size balloon.

For the minimum value, we chose a limiting velocity of -2 $m/s$ for the descent of the balloon, and calculated the lift force required to maintain that velocity using the model equations.

$v = (\frac{L - mg}{b})^{0.5}$

$L = bv^2 + mg$

# Calculating bounds for L
V = 4 / 3 * np.pi * radius**3

# Density of air inside the envelope at ~100 C. The literal needs the same units
# as rho_air, otherwise the subtraction is meaningless (and Pyomo will say so).
rho_hot = 0.946 * units.kg / units.m**3
L_ub = V * (rho_air - rho_hot) * g
print(f"V = {np.round(pyo.value(V), 2)} m^3")
print(f"L_ub = {np.round(pyo.value(L_ub), 2)} N")

# Lift needed to hold a 2 m/s descent: L = b*v*|v| + m*g with v = -2 m/s.
# The 2 here is a SPEED, so it carries m/s.
v_descent = 2 * units.m / units.s
L_lb = b * -1 * v_descent**2 + m_balloon * g
print(f"L_lb = {np.round(pyo.value(L_lb), 2)} N")

V = 4188.79 m^3
L_ub = 11464.68 N
L_lb = 5910.48 N

def create_balloon_model_v2(height_ref):
    """Fixed drag, added objective for difference of final height from height ref, squared.

    Input:
        height_ref: Reference desired final height (m).

        Returns:
        m: Initial model v2.
    """

    # Set a Target Height
    desired_height = height_ref * units.m
    endtime = 200  # seconds

    # A velocity-scale epsilon that removes the |v| singularity at v = 0.
    # It is compared against v, so it must carry v's units.
    eps = 1e-6 * units.m / units.s

    # Create a Pyomo Model
    m = pyo.ConcreteModel("Hot Air Balloon")

    # Define Time Domain. NOTE: a ContinuousSet cannot carry units (Pyomo issue
    # #1790), so the seconds live in this comment. Bounds are in SECONDS.
    m.t = dae.ContinuousSet(bounds=(0, endtime))

    # Define State Variables as Function of Time
    m.h = pyo.Var(m.t, bounds=(0, None), units=units.m)  # height
    # NOTE: Var bounds are NOT units-checked by assert_units_consistent.
    # These are m/s, and only this comment says so.
    m.v = pyo.Var(m.t, bounds=(-2, 5), units=units.m / units.s)
    m.D = pyo.Var(m.t, units=units.N)  # drag force

    # Define Control Variable (Lift) as a Function of Time
    m.L = pyo.Var(m.t, bounds=(L_lb, L_ub), units=units.N)

    # DerivativeVar defaults to DIMENSIONLESS and does NOT inherit units from the
    # state variable. m.t is time in seconds, so state it: dh/dt is m/s and
    # dv/dt is m/s^2.
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t, units=units.m / units.s)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t, units=units.m / units.s**2)

    # Define Derivative Constraints
    # Height derivative constraint (dh/dt = v)
    def dh_con(m, t):
        return m.dhdt[t] == m.v[t]

    m.dh_con = pyo.Constraint(m.t, rule=dh_con)

    # Velocity derivative, acceleration constraint (m * dv/dt = L - mg - D)
    def dv_con(m, t):
        return m.dvdt[t] * m_balloon == m.L[t] - m_balloon * g - m.D[t]

    m.dv_con = pyo.Constraint(m.t, rule=dv_con)

    # Drag force calculation (D = b * v^2)
    def drag_con(m, t):
        return m.D[t] == b * (m.v[t] ** 2 + eps**2) ** (1 / 2) * m.v[t]

    m.drag_con = pyo.Constraint(m.t, rule=drag_con)

    # Set Initial conditions
    m.h[0].fix(0)
    m.v[0].fix(0)
    m.L[0].fix(10000)  # newtons -- .fix() values are not units-checked either

    # Define the objective: Minimize the difference between final height and the desired height
    def objective_rule(m):
        return (
            m.h[endtime] - desired_height
        ) ** 2  # Minimize squared difference at final time

    m.objective = pyo.Objective(rule=objective_rule, sense=pyo.minimize)

    return m

# Create a Model Instance
ref_height = 600
NFE = 20
model = create_balloon_model_v2(ref_height)

# Call the Discretization/Initialization Function
discretize_initialize(model, NFE)

# Call the Solve Function
# To view Ipopt output, set tee = True
model = solve_model(model, tee=False)

# Visualize solution
visualize_model(model)

(<Figure size 1600x2560 with 3 Axes>,
 array([<Axes: ylabel='Height (m)'>, <Axes: ylabel='Velocity (m/s)'>,
        <Axes: xlabel='Time (s)', ylabel='Lift (kN)'>], dtype=object))

Analysis of Model Results: Reach a Desired Final Height¶

What is maximum achievable height in 200 seconds?¶

To answer this question, let’s test a few reference heights in the version 2 model.

NFE = 20
count = 0

# Set figure to hold plots
fig, axs = plt.subplots(2, 1, dpi=300, figsize=(8, 12.8))
colors = okabeito_list

# Assign reference heights
ref_heights = [100, 300, 500, 700]
for i in ref_heights:
    m = create_balloon_model_v2(i)
    discretize_initialize(m, NFE)
    solve_model(m, tee=False)
    t = [j for j in m.t]
    h = [m.h[j]() for j in m.t]
    v = [m.v[j]() for j in m.t]
    ls = linestyle_list[count % len(linestyle_list)]
    axs[0].plot(
        t,
        h,
        linewidth=3,
        label=f"ref_height = {i} m",
        color=colors[count],
        linestyle=ls,
    )
    axs[1].plot(t, v, linewidth=3, color=colors[count], linestyle=ls)
    plt.xlabel("Time (s)", fontweight="bold", fontsize=16)
    axs[0].set_ylabel("Height (m)", fontweight="bold", fontsize=16)
    axs[1].set_ylabel("Velocity (m/s)", fontweight="bold", fontsize=16)
    count += 1

# Format for publication quality
axs[0].legend()
axs[0].tick_params(direction="in", labelsize=15)
axs[1].tick_params(direction="in", labelsize=15)
plt.tight_layout()
plt.show()

From the graph, it looks like the maximum height is just under 700 m.

Check our work

With a quick calculation, let’s check the model solution with a hand calculation. The maximum height would use the maximum allowable lift.

Using our upper bound for lift, and the derived velocity from earlier in the notebook, the velocity achievable would be:

# Maximum velocity uses maximum lift.
#
# Note the units work out on their own: (N / (kg/m))**0.5 = m/s. That is a free
# check on the algebra of the hand calculation.
max_v = ((L_ub - m_balloon * g) / b) ** (0.5)
print(f"Maximum velocity = {np.round(pyo.value(max_v), 2)} {units.get_units(max_v)}")

Maximum velocity = 3.37 m/s

The height would then be $\int_{t_0}^{t_f} v(t) dt$ ,

or $t_f \times v(t)$ when $t_0 = 0$ .

Therefore, the maximum calculated height should be approximately 674 m, matching the plotted result.

Version 3: Follow a Linear Path¶

Now let’s solve the optimal control problem to find a control policy $u(t)$ for the interval $t_0 \leq t \leq t_f$ which causes the output $h(t)$ to track a desired setpoint or reference trajectory $r(t)$ .

\begin{align*} \min_{L(t)} ~~ & \int_{t_0}^{t_f} \|~h(t) - r(t)~\|^2\,dt \\ \mathrm{s.t.}~~ & \frac{dh}{dt} = v(t) \\ & D(t) = b \cdot v(t) \sqrt{v(t)^2 + \epsilon^2} \\ & \frac{dv}{dt} \cdot m = L(t) - m \cdot g - D(t) \\ & h(t_0) = 0 \\ & v(t_0) = 0 \\ & L(t_0) = 10000 \end{align*}

(8)

For the first path, we have the model follow a linear trajectory up defined by a slope.

# make a controlled path up


def linear(t):
    """Linear reference height: a 2.5 m/s climb.

    Inputs:
        t: Time (s), as a bare float -- this is the raw ContinuousSet index,
           which carries no units (Pyomo issue #1790).

    Returns:
        Reference height as a bare float, in METRES. The caller multiplies by
        units.m, because we cannot get units out of `t` itself."""

    return 2.5 * t

def create_balloon_model_v3():
    """Fixed drag, added integrated squared error between model and reference function, objective to minimize ise."""

    # A velocity-scale epsilon that removes the |v| singularity at v = 0.
    # It is compared against v, so it must carry v's units.
    eps = 1e-6 * units.m / units.s

    # Create a Pyomo Model
    m = pyo.ConcreteModel("Hot Air Balloon")

    # Define Time Domain. NOTE: a ContinuousSet cannot carry units (Pyomo issue
    # #1790), so the seconds live in this comment. Bounds are in SECONDS.
    m.t = dae.ContinuousSet(bounds=(0, tf))

    # Define State Variables as Function of Time
    m.h = pyo.Var(m.t, bounds=(0, None), units=units.m)  # height
    m.v = pyo.Var(m.t, units=units.m / units.s)  # vertical velocity
    m.D = pyo.Var(m.t, units=units.N)  # drag force

    # Define Control Variable (Lift) as a Function of Time
    m.L = pyo.Var(m.t, bounds=(L_lb, L_ub), units=units.N)

    # DerivativeVar defaults to DIMENSIONLESS and does NOT inherit units from the
    # state variable. m.t is time in seconds, so state it: dh/dt is m/s and
    # dv/dt is m/s^2.
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t, units=units.m / units.s)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t, units=units.m / units.s**2)

    # Define Derivative Constraints
    # Height derivative constraint (dh/dt = v)
    def dh_con(m, t):
        return m.dhdt[t] == m.v[t]

    m.dh_con = pyo.Constraint(m.t, rule=dh_con)

    # Velocity derivative, acceleration constraint (m * dv/dt = L - mg - D)
    def dv_con(m, t):
        return m.dvdt[t] * m_balloon == m.L[t] - m_balloon * g - m.D[t]

    m.dv_con = pyo.Constraint(m.t, rule=dv_con)

    # Drag force calculation (D = b * v^2)
    def drag_con(m, t):
        return m.D[t] == b * (m.v[t] ** 2 + eps**2) ** (1 / 2) * m.v[t]

    m.drag_con = pyo.Constraint(m.t, rule=drag_con)

    # Set Initial conditions
    m.L[0].fix(10000)  # newtons (.fix() values are not units-checked)

    # Define the objective: Minimize the difference between final height and the desired height
    @m.Integral(m.t)
    def ise(m, t):
        # linear(t) returns a bare float in metres: `t` is the raw ContinuousSet
        # index and carries no units, so we attach them here.
        return (m.h[t] - linear(t) * units.m) ** 2

    # Define the objective function
    @m.Objective(sense=pyo.minimize)
    def objective(m):
        return m.ise

    return m

# Create a Model Instance
tf = 90
nfe = 20
model = create_balloon_model_v3()

# Call the Discretization/Initialization Function
discretize_initialize(model, NFE)

# Call the Solve Function
model = solve_model(model)

# Display the results
# for t in model.t:
#     # if t % 2 == 0:
#         print(f"Time {t}: Height = {pyo.value(model.h[t])}, Velocity = {pyo.value(model.v[t])}, Lift = {pyo.value(model.L[t])}")

# Visualize solution
visualize_model(model)

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     1018
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      120

Total number of variables............................:      363
                     variables with only lower bounds:       60
                variables with lower and upper bounds:       60
                     variables with only upper bounds:        0
Total number of equality constraints.................:      303
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.5188628e+06 2.64e+03 1.01e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.5188502e+06 2.64e+03 1.01e+02  -1.0 1.77e+03    -  2.44e-03 8.71e-04f  1
   2  1.5186759e+06 2.63e+03 1.04e+02  -1.0 1.91e+03    -  1.06e-03 2.95e-03f  1
   3  1.5168410e+06 2.60e+03 1.03e+02  -1.0 2.04e+03    -  1.57e-03 1.34e-02f  1
   4  1.5167698e+06 2.60e+03 1.03e+02  -1.0 2.22e+03    -  1.87e-03 2.65e-04f  1
   5  1.4937278e+06 2.44e+03 6.37e+02  -1.0 2.30e+03    -  1.28e-03 5.87e-02f  1
   6  1.4876568e+06 2.41e+03 6.55e+02  -1.0 2.43e+03    -  2.04e-03 1.34e-02f  1
   7  1.4811581e+06 2.39e+03 6.57e+02  -1.0 2.44e+03    -  4.05e-03 8.86e-03f  1
   8  1.3458283e+06 2.01e+03 9.87e+02  -1.0 2.54e+03    -  6.84e-03 1.57e-01f  1
   9  1.3420699e+06 1.99e+03 9.76e+02  -1.0 1.86e+03  -2.0 1.05e-02 1.14e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.1623797e+06 2.50e+03 4.48e+03  -1.0 2.85e+03  -2.5 7.54e-03 3.38e-01f  1
  11  1.1512002e+06 1.80e+03 3.14e+03  -1.0 1.11e+03  -2.1 8.71e-02 2.88e-01f  1
  12  1.1587460e+06 1.74e+03 3.02e+03  -1.0 6.69e+02  -1.6 3.75e-02 3.52e-02h  1
  13  1.0209893e+06 1.88e+02 3.13e+03  -1.0 7.67e+02  -2.1 5.86e-04 1.00e+00f  1
  14  8.4096869e+05 5.04e+01 4.04e+03  -1.0 2.18e+02  -1.7 6.77e-02 1.00e+00f  1
  15  7.9700441e+05 1.27e+01 8.70e+02  -1.0 1.04e+02  -1.3 9.05e-01 1.00e+00f  1
  16  6.7846696e+05 4.08e+01 5.07e+02  -1.0 1.74e+02  -1.7 3.38e-01 1.00e+00f  1
  17  5.0445588e+05 8.79e+01 3.55e+02  -1.0 2.63e+02  -2.2 6.32e-01 1.00e+00f  1
  18  4.8052369e+05 5.90e+01 4.29e+02  -1.0 1.97e+02  -1.8 8.69e-01 5.49e-01f  1
  19  4.6855432e+05 5.21e+01 3.93e+02  -1.0 3.11e+02  -2.3 9.59e-01 1.07e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  2.6538853e+05 2.41e+03 5.68e+03  -1.0 1.21e+03  -2.7 8.39e-01 1.00e+00f  1
  21  2.3818489e+05 3.49e+02 7.61e+02  -1.0 7.50e+02  -2.3 2.30e-03 1.00e+00f  1
  22  1.5467211e+05 1.04e+02 6.44e+02  -1.0 3.45e+02  -2.8 9.73e-01 1.00e+00f  1
  23  1.3415801e+05 3.77e+01 3.94e+02  -1.0 1.88e+02  -2.4 1.00e+00 1.00e+00f  1
  24  9.0305340e+04 8.74e+01 3.32e+01  -1.0 3.13e+02  -2.8 1.00e+00 1.00e+00f  1
  25  7.8475811e+04 4.11e+01 1.43e+02  -1.0 1.60e+02  -2.4 1.00e+00 1.00e+00f  1
  26  5.3203594e+04 1.90e+02 3.71e+01  -1.0 3.48e+02  -2.9 1.00e+00 1.00e+00f  1
  27  2.4627608e+04 1.16e+02 1.63e+02  -1.0 5.63e+02  -3.4 3.66e-01 1.00e+00f  1
  28  1.9009784e+04 1.90e+01 2.65e+01  -1.0 1.53e+02  -2.9 1.00e+00 1.00e+00f  1
  29  9.7811777e+03 2.54e+02 5.64e+01  -1.0 4.78e+02  -3.4 7.29e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  7.4903736e+03 1.21e+02 9.45e+01  -1.0 2.05e+02  -3.0 1.00e+00 1.00e+00f  1
  31  4.7416733e+03 1.03e+02 6.61e+01  -1.0 2.45e+02  -3.5 1.00e+00 1.00e+00f  1
  32  2.2754878e+03 5.07e+02 1.78e+02  -1.0 8.22e+02  -3.9 1.00e+00 1.00e+00f  1
  33  6.4402017e+02 3.56e+03 7.39e+02  -1.0 4.18e+03  -4.4 1.75e-01 8.92e-01f  1
  34  6.2752648e+02 3.46e+03 7.17e+02  -1.0 1.57e+03  -4.9 1.00e+00 2.89e-02h  1
  35  2.3136066e+02 1.30e+03 2.92e+02  -1.0 1.98e+03  -5.4 9.92e-01 6.36e-01f  1
  36  5.7556626e+01 1.06e+02 2.17e+01  -1.0 8.67e+02  -5.9 1.00e+00 1.00e+00h  1
  37  3.1487527e+01 7.06e+01 8.01e-01  -1.0 9.42e+02  -6.3 5.75e-01 1.00e+00h  1
  38  1.3241605e+01 4.78e+02 2.24e+00  -1.0 3.27e+03    -  7.43e-01 8.63e-01h  1
  39  5.4301539e+00 2.10e+02 1.43e+00  -1.0 1.54e+03    -  5.15e-01 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  1.2734319e+00 1.19e+03 4.43e-01  -1.0 1.76e+03    -  1.00e+00 1.00e+00h  1
  41  4.0712458e-01 1.06e+02 2.98e-02  -1.0 1.61e+03    -  1.00e+00 1.00e+00h  1
  42  4.0696177e-01 2.26e+00 3.08e-04  -1.0 2.41e+02    -  1.00e+00 1.00e+00h  1
  43  4.0811904e-01 3.39e-02 1.42e-05  -1.0 2.55e+01    -  1.00e+00 1.00e+00h  1
  44  2.2549099e-01 1.96e+00 5.37e-03  -2.5 2.07e+02    -  9.28e-01 1.00e+00h  1
  45  3.6572299e-02 6.23e+00 1.02e-03  -2.5 5.62e+02    -  1.00e+00 1.00e+00h  1
  46  1.2564837e-02 1.30e+00 2.18e-04  -2.5 2.43e+02    -  1.00e+00 1.00e+00h  1
  47  1.0187450e-02 4.61e-02 8.98e-06  -2.5 4.65e+01    -  1.00e+00 1.00e+00h  1
  48  1.0102786e-02 4.30e-05 2.83e-08  -2.5 1.59e+00    -  1.00e+00 1.00e+00h  1
  49  1.3082483e-03 9.93e-01 1.78e-04  -3.8 2.09e+02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  1.6383481e-04 2.94e-01 2.43e-05  -3.8 1.06e+02    -  1.00e+00 1.00e+00h  1
  51  8.6281572e-05 6.29e-03 6.93e-07  -3.8 1.61e+01    -  1.00e+00 1.00e+00h  1
  52  8.4483514e-05 3.92e-06 1.50e-09  -3.8 4.17e-01    -  1.00e+00 1.00e+00h  1
  53  6.1634104e-07 2.97e-02 3.42e-06  -5.7 3.46e+01    -  1.00e+00 1.00e+00h  1
  54  1.5812176e-08 2.18e-04 2.24e-08  -5.7 2.95e+00    -  1.00e+00 1.00e+00h  1
  55  1.4889146e-08 5.77e-09 1.85e-11  -5.7 1.60e-02    -  1.00e+00 1.00e+00h  1
  56  9.2271233e-14 6.42e-06 7.29e-10  -8.6 5.08e-01    -  1.00e+00 1.00e+00h  1
  57  2.7530416e-14 9.55e-12 3.12e-14  -8.6 6.20e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 57

                                   (scaled)                 (unscaled)
Objective...............:   3.2761697935018366e-15    2.7530416414891780e-14
Dual infeasibility......:   3.1191097252860303e-14    2.6210604148534641e-13
Constraint violation....:   9.5496943686157465e-12    9.5496943686157465e-12
Complementarity.........:   2.5059065181173505e-09    2.1057715041936681e-08
Overall NLP error.......:   2.5059065181173505e-09    2.1057715041936681e-08


Number of objective function evaluations             = 58
Number of objective gradient evaluations             = 58
Number of equality constraint evaluations            = 58
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 58
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 57
Total CPU secs in IPOPT (w/o function evaluations)   =      0.024
Total CPU secs in NLP function evaluations           =      0.001

EXIT: Optimal Solution Found.

(<Figure size 1600x2560 with 3 Axes>,
 array([<Axes: ylabel='Height (m)'>, <Axes: ylabel='Velocity (m/s)'>,
        <Axes: xlabel='Time (s)', ylabel='Lift (kN)'>], dtype=object))

Version 4: Follow a Dynamic Flight Path¶

Lastly, we have the model follow a flight path containing multiple set heights and velocities.

First, using np.interp, we define the reference path.

# time grid
tf = 6000  # seconds; 100 minutes
dt = 0.5
n = round(tf / dt)
t_grid = np.linspace(0, tf, n + 1)


# setpoint/reference
def r(t):
    return np.interp(
        t,
        [0, 25, 550, 700, 1200, 1350, 1850, 2350, 2850, 4400, 5000, 5900, tf],
        [0, 100, 100, 400, 400, 900, 900, 600, 600, 300, 100, 100, 0],
    )


# plot the setpoint function
fig = plt.figure(dpi=300, figsize=(6.4, 4))
plt.plot(t_grid, r(t_grid), label="setpoint", color=okabeito_list[0], linewidth=3)
plt.xlabel("Time (s)", fontweight="bold", fontsize=16)
plt.ylabel("Height (m)", fontweight="bold", fontsize=16)
plt.tight_layout()
plt.tick_params(direction="in", labelsize=15)
plt.show()

r(t) is then used as the reference path in the previously set model.

def create_balloon_model_v4():

    # A velocity-scale epsilon that removes the |v| singularity at v = 0.
    # It is compared against v, so it must carry v's units.
    eps = 1e-6 * units.m / units.s

    # Create a Pyomo Model
    m = pyo.ConcreteModel("Hot Air Balloon")

    # Define Time Domain. NOTE: a ContinuousSet cannot carry units (Pyomo issue
    # #1790), so the seconds live in this comment. Bounds are in SECONDS.
    m.t = dae.ContinuousSet(bounds=(0, tf))

    # Define State Variables as Function of Time
    m.h = pyo.Var(m.t, bounds=(0, None), units=units.m)  # height
    m.v = pyo.Var(m.t, units=units.m / units.s)  # vertical velocity
    m.D = pyo.Var(m.t, units=units.N)  # drag force

    # Define Control Variable (Lift) as a Function of Time
    m.L = pyo.Var(m.t, bounds=(L_lb, L_ub), units=units.N)

    # DerivativeVar defaults to DIMENSIONLESS and does NOT inherit units from the
    # state variable. m.t is time in seconds, so state it: dh/dt is m/s and
    # dv/dt is m/s^2.
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t, units=units.m / units.s)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t, units=units.m / units.s**2)

    # Define Derivative Constraints
    # Height derivative constraint (dh/dt = v)
    def dh_con(m, t):
        return m.dhdt[t] == m.v[t]

    m.dh_con = pyo.Constraint(m.t, rule=dh_con)

    # Velocity derivative, acceleration constraint (m * dv/dt = L - mg - D)
    def dv_con(m, t):
        return m.dvdt[t] * m_balloon == m.L[t] - m_balloon * g - m.D[t]

    m.dv_con = pyo.Constraint(m.t, rule=dv_con)

    # Define Algebraic Constraints
    # Drag force calculation (D = b * v^2)
    def drag_con(m, t):
        return m.D[t] == b * (m.v[t] ** 2 + eps**2) ** (1 / 2) * m.v[t]

    m.drag_con = pyo.Constraint(m.t, rule=drag_con)

    # Define the objective: Minimize the difference between final height and the desired height
    @m.Integral(m.t)
    def ise(m, t):
        # r(t) interpolates the reference path and returns a bare float in metres;
        # `t` is the raw ContinuousSet index, so attach the units here.
        return (m.h[t] - r(t) * units.m) ** 2

    # Define the objective function
    @m.Objective(sense=pyo.minimize)
    def objective(m):
        return m.ise

    # Set Initial conditions
    m.L[0].fix(10000)  # newtons (.fix() values are not units-checked)

    return m

# Create model instance
tf = 6000
NFE = 60

model = create_balloon_model_v4()

# Call the Discretization/Initialization Function
discretize_initialize(model, NFE)

# Call the Solve Function
model = solve_model(model, tee=False)

# # Display the results
# for t in model.t:
#     if t % 300 == 0:
#         print(f"Time {t}: Height = {pyo.value(model.h[t])}, Velocity = {pyo.value(model.v[t])}, Lift = {pyo.value(model.L[t])}")

# Visualize solution
fig, axs = visualize_model(model)
# Add in plot for set height
h_set = [r(i) for i in model.t]
axs[0].plot(
    model.t,
    h_set,
    label="reference",
    linestyle=":",
    linewidth=3,
    color=okabeito_list[2],
)
axs[0].legend()
plt.tick_params(direction="in", labelsize=15)
plt.xlabel("Time (s)", fontweight="bold", fontsize=16)
plt.ylabel("Lift (kN)", fontweight="bold", fontsize=16)
plt.show()

While solution peaks look very sharp, keep in mind this is over a period of 6000s to simulate a full hot air balloon ride, which are usually around 90 minutes in length. Zooming in to the first quarter of the plot (0 - 1500s).

# Let us zoom in on key areas of the plot

t = [i for i in model.t]
L = [pyo.value(model.L[i]) / 1000 for i in model.t]  # lift in kN

figure = plt.figure(dpi=300, figsize=(6.4, 4))
plt.plot(t, L, linewidth=3, color=colors[0])
plt.xlim(-25, 1500)
plt.tick_params(direction="in", labelsize=15)
plt.xlabel("Time (s)", fontweight="bold", fontsize=16)
plt.ylabel("Lift (kN)", fontweight="bold", fontsize=16)

plt.tight_layout()

plt.show()

Conclusions¶

Take Away Messages¶

It is important to check that initial values and conditions make sense for the physical system.
Understanding the model on paper saves a ton of time coding into pyomo.
Examples like TCLab have a large potential to be modified for a variety of applications including a hot air balloon.

Future Work¶

Potential next steps for this model include:

Making a fuel estimation calculator. The lift would need to be related to the propane flowrate from the burner, and this would allow for a cost minimization tool in trip planning for hot air balloons.
Having the hot air balloon drag and lift calculations depend on height and temperature. This would add further complexity to the model and make the resulting curves more realistic.
Adding in horizontal movement from wind would allow for a multidimensional hot air balloon model.
Taking a specific design for a hot air balloon, a digital twin could be made to match the specifications and conduct uncertainty quantification against historic data.

Extensions of Referenced Material¶

Instead of defining a PID Controller, we implemented the equations of motion for a hot-air balloon in Pyomo and used an objective function to control the decision variable, lift.
The Charm Lab at Stanford only investigated the ascent of a balloon as it traveled to a set height. They did not account for the balloon moving downward; thus, their formulation could not be extended to path planning. We extended their formulation by modifying the Drag equation and defining a flight path for a balloon to reach a set of heights during the balloon’s ascent and descent.

References¶

(Hot Air Balloon Model) Bombik, A., & Mukhopadhaya, J. (2016). Controlling Team Rocket’s Hot Air Balloon. Stanford University. https://charm.stanford.edu/ENGR1052016/AnthonyBombikJayantMukhopadhaya
(Pyomo.DAE Example) Temperature Control Lab: Dowling, A., et al. (2023). Optimal control with Pyomo.DAE. Notre Dame CBE 30324/40324: Process Optimization. Retrieved October 7, 2024, from https://ndcbe.github.io/optimization/notebooks/3/PyomoDAE_TCLab.html
(Hot Air Balloon Physics: Real-World Physics Problems) (n.d.). Hot air balloon physics. Retrieved October 7, 2024, from https://www.real-world-physics-problems.com/hot-air-balloon-physics.html
(Engineering Toolbox: Calculating the Lift Force) The Engineering ToolBox. (n.d.). Hot air balloon lifting force. Retrieved October 7, 2024, from https://www.engineeringtoolbox.com/hot-air-balloon-lifting-force-d_562.html