Hot Air Balloon Dynamic Control

Hot Air Balloon Dynamic Control#

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

--2024-11-01 13:51:34--  https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.111.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 6493 (6.3K) [text/plain]
Saving to: ‘helper.py.9’


helper.py.9           0%[                    ]       0  --.-KB/s               
helper.py.9         100%[===================>]   6.34K  --.-KB/s    in 0s      

2024-11-01 13:51:34 (45.4 MB/s) - ‘helper.py.9’ saved [6493/6493]

idaes was found! No need to install.

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

Introduction#

In this notebook, we adapt the example in ../notebooks/3/PyomoDAE_TCLab.ipynb 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 principle 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 principle forces act can be seen below.

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*}\]

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*}\]

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*}\]

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: 364

Total Number of Equality Constraints: 364

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
rho_air = 1.225       # kg/m^3
g = 9.81              # m/s^2
radius = 10          # m
m_balloon = 750       # kg
cD = 0.47

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

# Print calculated values
print(f'S = {np.round(S, 2)} m^2')
print(f'b = {np.round(b, 2)} kg/m')

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 continous 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.

    '''

    # 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
    solver.solve(m, tee=tee)

    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

# 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
    m.t = dae.ContinuousSet(bounds=(0, 100))

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

    # Define State Variables as Functions of Time
    m.h = pyo.Var(m.t, bounds=(0, None))
    m.v = pyo.Var(m.t)
    m.D = pyo.Var(m.t)

    # Define Derivatives of the State Variables
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t)

    # 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
    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.53e-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.5279510989785194e-10    1.5279510989785194e-10
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.5279510989785194e-10    1.5279510989785194e-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.009
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))

../../_images/360654dd3d61d1587c300456fcaff6d5716ee0ec766b2b93f74bfb5bd7adce64.png

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 \(y^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 continuous, having a discontinuity at 0, so to model in Pyomo, we adjusted to make the function continuous by making \(|v| = \sqrt{v^2 + eps}\) 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*}\]

# Set-Up our Model as Function
def create_balloon_model_v1():

    ''' Modified from v0. Fixed drag.'''

    # Define a value ~0 to Help Stabilize the Drag Force Differential Equation
    eps = 1e-6

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

    # Define time Domain
    m.t = dae.ContinuousSet(bounds=(0, 100))

    # Define State Variables as Functions of Time
    m.h = pyo.Var(m.t, bounds=(0, None))
    m.v = pyo.Var(m.t)
    m.D = pyo.Var(m.t)

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

    # Define derivatives
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t)


    # 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
    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.5370460459962487e-10    1.5370460459962487e-10
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.5370460459962487e-10    1.5370460459962487e-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.010
Total CPU secs in NLP function evaluations           =      0.001

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.'''

    # Define a Samll Constant to Help Stablize the Drag Force Differential Equation
    eps = 1e-6

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

    # Define time Domain
    m.t = dae.ContinuousSet(bounds=(0, 100))

    #Define State Variables as Functions of Time
    m.h = pyo.Var(m.t, bounds=(0, None))
    m.v = pyo.Var(m.t)
    m.D = pyo.Var(m.t)

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

    # Define derivatives
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t)


    # 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

# 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
    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]
    axs[0].plot(t, h, linewidth = 3, label = f'L = {i} kN', color = colors[count])
    axs[1].plot(t, v, linewidth = 3, color = colors[count])
    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()

../../_images/1ad62aecf554be12f30d4a33223d3e85858baf3fcba00fc413c9890a0a24ecf3.png

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*}\]

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 involving 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_{u(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*}\]

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 = 0 \)

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_{u(t)} ~~ & \int_{t_0}^{t_f} \|~h(t_f) - h_{desired}~\|^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*}\]

In the next iteration of the model, we add an objective to minimize the sum of squared error (SSE) between a target height 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
L_ub = V*(rho_air-0.946)*g
print(f'V = {np.round(V, 2)} m^3')
print(f'L_ub = {np.round(L_ub, 2)} N')

L_lb = b*-1*(2)**2 + m_balloon*g
print(f'L_lb = {np.round(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, sqaured.

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

        Returns:
        m: Initial model v2.
        '''

    # Set a Target Height
    desired_height = height_ref
    endtime = 200

    # Define a Small Constant to Help Stabilize the Drag Force Constraint
    eps = 1e-6

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

    # Define Time Domain
    m.t = dae.ContinuousSet(bounds=(0, endtime))

    # Define State Variables as Function of Time
    m.h = pyo.Var(m.t, bounds=(0, None))
    m.v = pyo.Var(m.t, bounds=(-2, 5))
    m.D = pyo.Var(m.t)

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

    # Define Derivative Variables
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t)

    # 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)

    # 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))

../../_images/bf8129168e8082b9b1682bc444af0d3341c79b639c9066ea43f31e9680de6f3c.png

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 model.t]
    h = [m.h[j]() for j in m.t]
    v = [m.v[j]() for j in m.t]
    axs[0].plot(t, h, linewidth = 3, label = f'ref_height = {i} m', color = colors[count])
    axs[1].plot(t, v, linewidth = 3, color = colors[count])
    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()

../../_images/0233761c5f189837ece0053272287371f99ef897a2672f8cc0e94fbc7fe53ed0.png

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
max_v = ((L_ub - m_balloon*g)/b)**(0.5)
print(f'Maximum velocity = {np.round(max_v, 2)} m/s')

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 \(T_H(t)\) to track a desired setpoint or reference tracjectory \(SP(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*}\]

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 function with slope 2.5.
    Inputs:
        t: Time (s)

    Returns
        2.5*t:'''

    return 2.5*t

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

    # Define a Samll Constant to Help Stabalize the Drag Force Differential Equation
    eps = 1e-6

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

    # Define Time Domain
    m.t = dae.ContinuousSet(bounds=(0, tf))

    # Define State Variables as Function of Time
    m.h = pyo.Var(m.t, bounds=(0, None))
    m.v = pyo.Var(m.t)
    m.D = pyo.Var(m.t)

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

    # Define Derivative Variables
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t)

    # 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)

    # Define the objective: Minimize the difference between final height and the desired height
    @m.Integral(m.t)
    def ise(m, t):
        return (m.h[t] - linear(t))**2  # Minimize squared difference at final time


    # 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.09e-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.1191096832652552e-14    2.6210603795424299e-13
Constraint violation....:   9.0949470177292824e-12    9.0949470177292824e-12
Complementarity.........:   2.5059065181173699e-09    2.1057715041936844e-08
Overall NLP error.......:   2.5059065181173699e-09    2.1057715041936844e-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.093
Total CPU secs in NLP function evaluations           =      0.006

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))

../../_images/f8586fc80c964f799b2bb8cc1b54adbf2a03fc48d999817910c9397c2dec9e3b.png

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 #sec in 100 minuts 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()

../../_images/a20d36294ba14886ea00ccbfac774f2b010225342fddf1264db86becb1760daa.png

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

def create_balloon_model_v4():


    # Define a Samll Constant to Help Stabalize the Drag Force Differential Equation
    eps = 1e-6

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

    # Define Time Domain
    m.t = dae.ContinuousSet(bounds=(0, tf))

    # Define State Variables as Function of Time
    m.h = pyo.Var(m.t, bounds=(0, None))
    m.v = pyo.Var(m.t)
    m.D = pyo.Var(m.t)

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

    # Define Derivative Variables
    m.dhdt = dae.DerivativeVar(m.h, wrt=m.t)
    m.dvdt = dae.DerivativeVar(m.v, wrt=m.t)

    # 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 Contraints
    # 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):
        return (m.h[t] - r(t))**2  # Minimize squared difference at final time


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

    # Set Initial conditions
    m.L[0].fix(10000)

    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()

../../_images/01f8e7df5d81399d2cc4b21750e64e85b4fb980d415ec76cb95e8098a1874986.png

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()

../../_images/53e2792a86926d8d7faf58f3c1045748b348537de68b8baff3832dca246412b9.png

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 ballon 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

Hot Air Balloon Dynamic Control

Contents

Hot Air Balloon Dynamic Control#

Learning Objectives#

Introduction#

Hot Air Balloons#

Model Formulation: No Set Objective#

Degree of Freedom Analysis#

Define Model Constants#

Function Definition#

Version 0: Simulate the Balloon Dynamics#

Version 1: Correcting the Drag Force Calculation#

Analysis of Model Results: Varying Fixed Lift Value#

How does the set lift affect the velocity?#

Model Formulation - Minimize ISE Objective#

Degree of Freedom Analysis#

Version 2: Reach a Desired Final Height#

Analysis of Model Results: Reach a Desired Final Height#

What is maximum achievable height in 200 seconds?#

Version 3: Follow a Linear Path#

Version 4: Follow a Dynamic Flight Path#

Conclusions#

Take Away Messages#

Future Work#

Extensions of Referenced Material#

References#