Notebook source code: examples/goddard/goddard.ipynb
Run the notebook yourself on binder

Goddard problem (python)

Thumbnail

This well-known problem (see for instance [1],[2]) models the ascent of a rocket through the atmosphere, and we restrict here ourselves to vertical (monodimensional) trajectories. The state variables are the altitude, speed and mass of the rocket during the flight, for a total dimension of 3. The rocket is subject to gravity, thrust and drag forces. The final time is fixed, and the objective is to reach a maximal altitude with a prescribed fuel consumption, ie a fixed final mass.

\[\begin{split}\left \lbrace \begin{array}{l} Max\ h(T)\\ \dot h = v\\ \dot v = \frac{c u - D(h,v)}{m} - g_0\\ \dot m = - u\\ u(\cdot) \in [0,9.52551]\\ h(0)=0,\ v(0)=0,\ m(0)=214.839\\ m(T) = 67.983310\\ \end{array} \right .\end{split}\]

The drag term is a function of speed and altitude defined as \(D(h,v)=\alpha v^2\rho(h)\), with the volumic mass given by the approximate model \(\rho(h)=e^{-\beta h}\).

In the following we use the parameters \(T=206.661,\alpha=0.01227,\beta=0.000145, g_0=9.81, c=2060\).

The Hamiltonian is an affine function of the control, so singular arcs may occur.

References

1. R.H. Goddard. A Method of Reaching Extreme Altitudes, volume 71(2) of Smithsonian Miscellaneous Collections. Smithsonian institution, City of Washington, 1919.

2. 8. Seywald and E.M. Cliff. Goddard problem in presence of a dynamic pressure limit. Journal of Guidance, Control, and Dynamics, 16(4):776–781, 1993.

!pygmentize problem.cpp

// +++DRAFT+++ This class implements the OCP functions
// It derives from the generic class bocop3OCPBase
// OCP functions are defined with templates since they will be called
// from both the NLP solver (double arguments) and AD tool (ad_double arguments)
//#pragma once

#include <OCP.h>
// ///////////////////////////////////////////////////////////////////


// ///////////////////////////////////////////////////////////////////


template <typename Variable>
void OCP::finalCost(double initial_time, double final_time, const Variable *initial_state, const Variable *final_state, const Variable *parameters, const double *constants, Variable &final_cost)
{
  // CRITERION FOR GODDARD PROBLEM
  // MAXIMIZE FINAL HEIGHT
  final_cost = -final_state[0];
}

template <typename Variable>
void OCP::dynamics(double time, const Variable *state, const Variable *control, const Variable *parameters, const double *constants, Variable *state_dynamics)
{
  // DYNAMICS FOR GODDARD PROBLEM
  // dh/dt = v
  // dv/dt = (c u - Drag(h,v)) / m - g0
  // dm/dt =-u

  double tf = constants[0];
  double alpha = constants[1];
  double beta = constants[2];
  double g0 = constants[3];
  double c = constants[4];

  Variable h = state[0];
  Variable v = state[1];
  Variable m = state[2];
  Variable u = control[0];

  // drag term
  Variable D = alpha*v*v * exp(-beta*h);

  state_dynamics[0] = v;
  state_dynamics[1] = (c*u - D) / m - g0;
  state_dynamics[2] = - u;

  // rescale dynamics from [0,tf] to [0,1]
  for (size_t i=0; i<3; i++)
    state_dynamics[i] *= tf;
}

template <typename Variable>
void OCP::boundaryConditions(double initial_time, double final_time, const Variable *initial_state, const Variable *final_state, const Variable *parameters, const double *constants, Variable *boundary_conditions)
{
  // INITIAL/FINAL CONDITIONS FOR GODDARD PROBLEM
  // r0, v0, m0, and mf
  boundary_conditions[0] = initial_state[0];
  boundary_conditions[1] = initial_state[1];
  boundary_conditions[2] = initial_state[2];
  boundary_conditions[3] = final_state[2];
}

template <typename Variable>
void OCP::pathConstraints(double time, const Variable *state, const Variable *control, const Variable *parameters, const double *constants, Variable *path_constraints)
{}

void OCP::preProcessing()
{}

// ///////////////////////////////////////////////////////////////////
// explicit template instanciation for template functions, with double and double_ad 
// +++ could be in an included separate file ? 
// but needs to be done for aux functions too ? APPARENTLY NOT !
template void OCP::finalCost<double>(double initial_time, double final_time, const double *initial_state, const double *final_state, const double *parameters, const double *constants, double &final_cost);
template void OCP::dynamics<double>(double time, const double *state, const double *control, const double *parameters, const double *constants, double *state_dynamics);
template void OCP::boundaryConditions<double>(double initial_time, double final_time, const double *initial_state, const double *final_state, const double *parameters, const double *constants, double *boundary_conditions);
template void OCP::pathConstraints<double>(double time, const double *state, const double *control, const double *parameters, const double *constants, double *path_constraints);

template void OCP::finalCost<double_ad>(double initial_time, double final_time, const double_ad *initial_state, const double_ad *final_state, const double_ad *parameters, const double *constants, double_ad &final_cost);
template void OCP::dynamics<double_ad>(double time, const double_ad *state, const double_ad *control, const double_ad *parameters, const double *constants, double_ad *state_dynamics);
template void OCP::boundaryConditions<double_ad>(double initial_time, double final_time, const double_ad *initial_state, const double_ad *final_state, const double_ad *parameters, const double *constants, double_ad *boundary_conditions);
template void OCP::pathConstraints<double_ad>(double time, const double_ad *state, const double_ad *control, const double_ad *parameters, const double *constants, double_ad *path_constraints);

%matplotlib inline
import bocop
print(bocop.__version__)

3.1.0

problem_path = "." # using local problem definition
bocop.build(problem_path)#, cmake_options = '-DCMAKE_CXX_COMPILER=g++')

[EXEC] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/Users/aherasim/ct-gallery/examples/goddard /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop']
>    -- The C compiler identification is AppleClang 12.0.0.12000032
>    -- The CXX compiler identification is AppleClang 12.0.0.12000032
>    -- Detecting C compiler ABI info
>    -- Detecting C compiler ABI info - done
>    -- Check for working C compiler: /Library/Developer/CommandLineTools/usr/bin/cc - skipped
>    -- Detecting C compile features
>    -- Detecting C compile features - done
>    -- Detecting CXX compiler ABI info
>    -- Detecting CXX compiler ABI info - done
>    -- Check for working CXX compiler: /Library/Developer/CommandLineTools/usr/bin/c++ - skipped
>    -- Detecting CXX compile features
>    -- Detecting CXX compile features - done
>    -- Problem path: /Users/aherasim/ct-gallery/examples/goddard
>    -- Using CPPAD found at /Users/aherasim/anaconda3/envs/ct/include/cppad/..
>    -- Using IPOPT found at /Users/aherasim/anaconda3/envs/ct/lib/libipopt.dylib
>    -- Found Python3: /Users/aherasim/anaconda3/envs/ct/bin/python3.7 (found version "3.7.10") found components: Interpreter Development Development.Module Development.Embed
>    -- Found SWIG: /Users/aherasim/anaconda3/envs/ct/bin/swig (found suitable version "4.0.2", minimum required is "4")
>    -- Build type: Debug
>    -- Configuring done
>    -- Generating done
>    -- Build files have been written to: /Users/aherasim/ct-gallery/examples/goddard/build
[DONE] > ['cmake -DCMAKE_BUILD_TYPE=Debug -DPROBLEM_DIR=/Users/aherasim/ct-gallery/examples/goddard /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop']
[EXEC] > make
>    [  3%] Building CXX object src/CMakeFiles/bocop.dir/AD/dOCPCppAD.cpp.o
>    [  7%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dOCP.cpp.o
>    [ 10%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dODE.cpp.o
>    [ 14%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dControl.cpp.o
>    [ 17%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/dState.cpp.o
>    [ 21%] Building CXX object src/CMakeFiles/bocop.dir/DOCP/solution.cpp.o
>    [ 25%] Building CXX object src/CMakeFiles/bocop.dir/NLP/NLPSolverIpopt.cpp.o
>    /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/NLP/NLPSolverIpopt.cpp:38:8: warning: 'finalize_solution' overrides a member function but is not marked 'override' [-Winconsistent-missing-override]
>      void finalize_solution(Ipopt::SolverReturn status, Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U, Ipopt::Index m, const Ipopt::Number* g,
>           ^
>    /Users/aherasim/anaconda3/envs/ct/include/coin/IpTNLP.hpp:534:17: note: overridden virtual function is here
>       virtual void finalize_solution(
>                    ^
>    1 warning generated.
>    [ 28%] Building CXX object src/CMakeFiles/bocop.dir/OCP/OCP.cpp.o
>    [ 32%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools.cpp.o
>    [ 35%] Building CXX object src/CMakeFiles/bocop.dir/tools/tools_interpolation.cpp.o
>    [ 39%] Building CXX object src/CMakeFiles/bocop.dir/Users/aherasim/ct-gallery/examples/goddard/problem.cpp.o
>    [ 42%] Linking CXX shared library /Users/aherasim/ct-gallery/examples/goddard/libbocop.dylib
>    ld: warning: -pie being ignored. It is only used when linking a main executable
>    [ 42%] Built target bocop
>    Scanning dependencies of target bocopwrapper_swig_compilation
>    [ 46%] Swig compile /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/bocopwrapper.i for python
>    Language subdirectory: python
>    Search paths:
>       ./
>       /Users/aherasim/anaconda3/envs/ct/include/python3.7m/
>       AD/
>       DOCP/
>       NLP/
>       OCP/
>       tools/
>       /Users/aherasim/anaconda3/envs/ct/include/cppad/../
>       /Users/aherasim/anaconda3/envs/ct/include/coin/
>       /Users/aherasim/anaconda3/envs/ct/include/python3.7m/
>       /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/AD/
>       /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/DOCP/
>       /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/NLP/
>       /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/OCP/
>       /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/tools/
>       ./swig_lib/python/
>       /Users/aherasim/anaconda3/envs/ct/share/swig/4.0.2/python/
>       ./swig_lib/
>       /Users/aherasim/anaconda3/envs/ct/share/swig/4.0.2/
>    Preprocessing...
>    Starting language-specific parse...
>    Processing types...
>    C++ analysis...
>    Processing nested classes...
>    Generating wrappers...
>    [ 46%] Built target bocopwrapper_swig_compilation
>    [ 50%] Building CXX object src/CMakeFiles/bocopwrapper.dir/CMakeFiles/bocopwrapper.dir/bocopwrapperPYTHON_wrap.cxx.o
>    [ 53%] Linking CXX shared module /Users/aherasim/ct-gallery/examples/goddard/_bocopwrapper.so
>    ld: warning: -pie being ignored. It is only used when linking a main executable
>    -- Moving python modules to /Users/aherasim/ct-gallery/examples/goddard
>    [ 53%] Built target bocopwrapper
>    [ 57%] Building CXX object src/CMakeFiles/bocopApp.dir/AD/dOCPCppAD.cpp.o
>    [ 60%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dOCP.cpp.o
>    [ 64%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dODE.cpp.o
>    [ 67%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dControl.cpp.o
>    [ 71%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/dState.cpp.o
>    [ 75%] Building CXX object src/CMakeFiles/bocopApp.dir/DOCP/solution.cpp.o
>    [ 78%] Building CXX object src/CMakeFiles/bocopApp.dir/NLP/NLPSolverIpopt.cpp.o
>    /Users/aherasim/anaconda3/envs/ct/lib/python3.7/site-packages/bocop/src/NLP/NLPSolverIpopt.cpp:38:8: warning: 'finalize_solution' overrides a member function but is not marked 'override' [-Winconsistent-missing-override]
>      void finalize_solution(Ipopt::SolverReturn status, Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U, Ipopt::Index m, const Ipopt::Number* g,
>           ^
>    /Users/aherasim/anaconda3/envs/ct/include/coin/IpTNLP.hpp:534:17: note: overridden virtual function is here
>       virtual void finalize_solution(
>                    ^
>    1 warning generated.
>    [ 82%] Building CXX object src/CMakeFiles/bocopApp.dir/OCP/OCP.cpp.o
>    [ 85%] Building CXX object src/CMakeFiles/bocopApp.dir/tools/tools.cpp.o
>    [ 89%] Building CXX object src/CMakeFiles/bocopApp.dir/tools/tools_interpolation.cpp.o
>    [ 92%] Building CXX object src/CMakeFiles/bocopApp.dir/Users/aherasim/ct-gallery/examples/goddard/problem.cpp.o
>    [ 96%] Building CXX object src/CMakeFiles/bocopApp.dir/main.cpp.o
>    [100%] Linking CXX executable /Users/aherasim/ct-gallery/examples/goddard/bocopApp
>    [100%] Built target bocopApp
[DONE] > make
Done

0

bocop.run(problem_path, graph=1)

Done

%matplotlib widget
import matplotlib.pyplot as plt

solution = bocop.readSolution(problem_path + "/problem.sol")
bocop.low_diagnose(solution)

print("Bocop returns status {} with objective {:2.4g} and constraint violation {:2.4g}".format(solution.status,solution.objective,solution.constraints))
p0 = []
for i in range(solution.dim_state):
    p0.append(solution.costate[i][0])
print("Costate at first time stage (t0+h/2): ",p0)
print("Multipliers for initial conditions: ",solution.boundarycond_multipliers[0:solution.dim_state])

previous

Biomass maximization with substrate depletion

next

Goddard problem (julia)