Goddart's rocket problem (Bocop)

This is a Bocop implementation of Goddart's rocket problem.

Problem.def

# This file defines all dimensions and parameters
# values for your problem :
 
# Initial and final time :
time.free string final
time.initial double 0
time.final double 1
 
# Dimensions :
state.dimension integer 3
control.dimension integer 1
algebraic.dimension integer 0
parameter.dimension integer 1
constant.dimension integer 6
boundarycond.dimension integer 4
constraint.dimension integer 1
 
# Discretization :
discretization.steps integer 100
discretization.method string gauss
 
# Optimization :
optimization.type string batch
batch.type integer 0
batch.index integer 5
batch.nrange integer 3
batch.lowerbound double 0.3
batch.upperbound double 0.9
batch.directory string Batch-C
 
# Initialization :
initialization.type string from_init_file
initialization.file string none
 
# Parameter identification :
paramid.type string false
paramid.separator string ,
paramid.file string no_directory
paramid.dimension integer 0
 
# Names :
state.0 string position
state.1 string speed
state.2 string mass
control.0 string acceleration_u
parameter.0 string finaltime
boundarycond.0 string r(0)
boundarycond.1 string v(0)
boundarycond.2 string m(0)
boundarycond.3 string r(f)
constraint.0 string drag_minus_C
constant.0 string Tmax
constant.1 string A
constant.2 string k
constant.3 string r0
constant.4 string b
constant.5 string C
 
# Solution file :
solution.file string problem.sol
 
# Iteration output frequency :
iteration.output.frequency integer 0

Problem.constants

# This file contains the values of the constants of your problem.
# Number of constants used in your problem : 
6
 
# Values of the constants : 
3.5
310
500
1
7
0.6

Problem.bounds

# This file contains all the bounds of your problem.
# Bounds are stored in standard format : 
# [lower bound]  [upper bound] [type of bound]
 
# Dimensions (i&f conditions, y, u, z, p, path constraints) :
4 3 1 0 1 1
 
# Bounds for the initial and final conditions :
1 1 equal
0 0 equal
1 1 equal
1.01 2e+020 lower
 
# Bounds for the state variables :
>0:1:0 1 2e+020 lower
>1:1:1 0 2e+020 lower
>2:1:2 0 2e+020 lower
 
# Bounds for the control variables :
0 1 both
 
# Bounds for the algebraic variables :
 
# Bounds for the optimization parameters :
0 2e+020 lower
 
# Bounds for the path constraints :
-2e+020 0 upper

criterion.tpp

#include "header_criterion"
{
	// CRITERION FOR GODDARD PROBLEM
 
	// MAXIMIZE FINAL MASS
	criterion = -final_state[2];
}

dynamics.tpp

#include "header_dynamics"
{
	// DYNAMICS FOR GODDARD PROBLEM
	// dr/dt = v
	// dv/dt = (Thrust(u) - Drag(r,v)) / m - grav(r)
	// dm/dt = -b*|u|
 
	double Tmax = constants[0];
	double A = constants[1];
	double k = constants[2];
	double r0 = constants[3];
	double b = constants[4];
 
	Tdouble r = state[0];
	Tdouble v = state[1];
	Tdouble m = state[2];
 
	state_dynamics[0] = v;
	state_dynamics[1] = (thrust(control[0],Tmax) - drag(r,v,A,k,r0)) / m - grav(r);
	state_dynamics[2] = - b * control[0];
}

boundarycond.tpp

#include "header_boundarycond"
{
	// INITIAL CONDITIONS FOR GODDARD PROBLEM
	// r0 = 1    v0 = 0   m0 = 1
	// MODELED AS 1 <= r0 <= 1, etc
	boundary_conditions[0] = initial_state[0];
	boundary_conditions[1] = initial_state[1];
	boundary_conditions[2] = initial_state[2];
 
	// FINAL CONDITIONS FOR GODDARD PROBLEM
  	// rf >= 1.01   MODELED AS   1.01 <= rf
	boundary_conditions[3] = final_state[0];
}

pathcond.tpp

#include "header_pathcond"
{
	// CONSTRAINT ON MAX DRAG FOR GODDARD PROBLEM
	// Drag <= C ie Drag - C <= 0
 
	double A = constants[1];
	double k = constants[2];
	double r0 = constants[3];
	double C = constants[5];
 
	Tdouble r = state[0];
	Tdouble v = state[1];
 
	path_constraints[0] = drag(r,v,A,k,r0) - C;
}

dependencies.tpp

#include "./grav.tpp"
#include "./drag.tpp"
#include "./thrust.tpp"

drag.tpp

#include "adolc/adolc.h"
#include "adolc/adouble.h"
#include <cmath>
 
// FUNCTION FOR GODDARD DRAG
// drag = 310 v^2 exp (-500(r-1))
 
// Arguments:
// r: radius
// v: velocity
 
template<class Tdouble> Tdouble drag(Tdouble r, Tdouble v, double A, double k, double r0)
{
 
	Tdouble drag;
	drag = A * v * v * exp(-k*(fabs(r)-r0));
	return drag;
 
}

grav.tpp

#include "adolc/adolc.h"
#include "adolc/adouble.h"
#include <cmath>
 
// FUNCTION FOR GRAVITY 
// g = 1 / r^2
 
// Arguments:
// r: radius
 
template<class Tdouble> Tdouble grav(Tdouble r)
{
 
	Tdouble grav;
	grav = 1e0 / r / r;
	return grav;
 
}

thrust.tpp

#include "adolc/adolc.h"
#include "adolc/adouble.h"
#include <cmath>
 
// FUNCTION FOR THRUST (GODDARD)
// T = u * Tmax
 
// Arguments:
// r: radius
 
template<class Tdouble> Tdouble thrust(Tdouble u, double Tmax)
{
	Tdouble thrust;
	thrust = u * Tmax;
	return thrust;
}