diff --git a/developerGuide/extendingOptimization/MissionLinker.svg b/developerGuide/extendingOptimization/MissionLinker.svg
new file mode 100644
index 0000000000..637d870daf
--- /dev/null
+++ b/developerGuide/extendingOptimization/MissionLinker.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/developerGuide/extendingOptimization/Missionsegmentblock.svg b/developerGuide/extendingOptimization/Missionsegmentblock.svg
new file mode 100644
index 0000000000..c4ac3230c7
--- /dev/null
+++ b/developerGuide/extendingOptimization/Missionsegmentblock.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/developerGuide/extendingOptimization/basics.rst b/developerGuide/extendingOptimization/basics.rst
new file mode 100644
index 0000000000..01decb8c4d
--- /dev/null
+++ b/developerGuide/extendingOptimization/basics.rst
@@ -0,0 +1,83 @@
+.. _optimization_basics:
+
+Basics
+======
+
+.. _optimization-linearInterpolation:
+
+Linear interpolation
+~~~~~~~~~~~~~~~~~~~~
+
+The interpolation used by the linking methods is linear. Assume we have two available epochs
+:math:`t_0` and :math:`t_1` and the respective function values :math:`f(t_0)` and :math:`f(t_1)`.
+If we want to retrieve the value of :math:`f` at :math:`t_e`, with
+
+ .. math::
+
+ t_0 < t_e < t_1
+
+then, assuming that :math:`e = t_1 - t_0` is small enough, the value can be retrieved as following:
+
+ .. math::
+
+ f(t_e) = f(t_0) + (f(t_1) - f(t_0))(t_e - t_0)/(t_1 - t_0)
+
+.. _optimization-dynamicPointerCast:
+
+Dynamic pointer cast
+~~~~~~~~~~~~~~~~~~~~
+
+Let us consider the class ``DecisionVariableSettings``:
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/decisionVariableSettings.h`
+ :name: decisionVariableSettings-h
+
+ struct DecisionVariableSettings{
+
+ DecisionVariableSettings( std::vector< boost::shared_ptr< SingleDecisionVariableSettings > >
+ multiDecisionVariableSettings ) : decisionVariableSettings_( multiDecisionVariableSettings )
+ { }
+
+ DecisionVariableSettings( boost::shared_ptr< SingleDecisionVariableSettings >
+ singleDecisionVariableSettings )
+ {
+ decisionVariableSettings_.push_back( singleDecisionVariableSettings );
+ }
+
+ std::vector< boost::shared_ptr< SingleDecisionVariableSettings > > decisionVariableSettings_;
+
+ };
+
+The member ``decisionVariableSettings_`` is a vector of ``boost::shared_ptr< SingleDecisionVariableSettings >``,
+but any index can be as well initialized as ``boost::shared_ptr< >`` of any of the
+``SingleDecisionVariableSettings`` derived classes.
+
+For example, if the first member of the constructor parameter ``multiDecisionVariableSettings`` is of the type
+``boost::shared_ptr< SingleKeplerElementDecisionVariableSettings >`` the value is accepted and stored as
+``boost::shared_ptr< SingleDecisionVariableSettings >``.
+
+``SingleKeplerElementDecisionVariableSettings`` contains methods and members that do not belong to
+``SingleDecisionVariableSettings``. In order to access these methods one has to perform a **dynamic pointer cast**.
+
+Let us stick to the example above. That means that ``decisionVariableSettings_[0]`` is stored as
+``boost::shared_ptr< SingleDecisionVariableSettings >`` but contains information of ``SingleKeplerElementDecisionVariableSettings``.
+
+This can be accertained with the flag ``SingleDecisionVariableSettings::decisionVariable_``, which in this case is set
+automatically by the constructor of ``SingleKeplerElementDecisionVariableSettings`` to be:
+
+ .. code-block:: cpp
+
+ decisionVariableSettings_[0]->decisionVariable_ == single_kepler_element_decision_variable;
+
+The dynamic cast to retrieve the information of the derived class is then performed as following:
+
+
+ .. code-block:: cpp
+
+ boost::shared_ptr< SingleKeplerElementDecisionVariableSettings > dynamicCastDecisionVariable =
+ boost::dynamic_pointer_cast< SingleKeplerElementDecisionVariableSettings >( decisionVariableSettings_[0] );
+
+
+Now all the ``SingleKeplerElementDecisionVariableSettings`` are available inside ``dynamicCastDecisionVariable``.
+
diff --git a/developerGuide/extendingOptimization/fluxDiagram.svg b/developerGuide/extendingOptimization/fluxDiagram.svg
new file mode 100644
index 0000000000..2cc3f02adb
--- /dev/null
+++ b/developerGuide/extendingOptimization/fluxDiagram.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/developerGuide/extendingOptimization/index.rst b/developerGuide/extendingOptimization/index.rst
new file mode 100644
index 0000000000..2ad6f3b759
--- /dev/null
+++ b/developerGuide/extendingOptimization/index.rst
@@ -0,0 +1,15 @@
+.. _extendingOptimization:
+
+Extending the optimization namespace
+====================================
+
+The goal of this page is to document how the :ref:`tudatOptimizationIndex` can be updated to support new/additional features of Tudat.
+
+.. toctree::
+ :maxdepth: 2
+
+ basics
+ objectiveFunction
+ missionSegment
+ missionLinker
+ necessaryFeatures
diff --git a/developerGuide/extendingOptimization/missionLinker.rst b/developerGuide/extendingOptimization/missionLinker.rst
new file mode 100644
index 0000000000..bf3e367260
--- /dev/null
+++ b/developerGuide/extendingOptimization/missionLinker.rst
@@ -0,0 +1,382 @@
+.. _extendMissionLinker:
+
+Mission Linker
+==============
+
+What is the Mission Linker?
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The mission linker, here represented by the class ``MissionLinker`` is a tool by whih a number of
+of mission segments (class ``MissionSegmentSettings``) are linked together in a chronological fashion.
+
+The structure to define such an object is the following:
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/missionLinker.h`
+ :name: missionLinker-h
+
+ struct MissionLinker
+ {
+
+ MissionLinker( std::vector< boost::shared_ptr< MissionSegmentSettings > > missionSegmentSettings,
+ boost::shared_ptr< OptimizationSettings > optimizationSettings = boost::make_shared< OptimizationSettings >(),
+ const bool optimizeConfiguration = true )
+ : missionSegmentSettings_( missionSegmentSettings ), optimizationSettings_(optimizationSettings)
+ {
+
+ if( optimizeConfiguration )
+ optimize();
+ }
+
+ ~MissionLinker( ){ }
+
+ std::vector< boost::shared_ptr< MissionSegmentSettings > > missionSegmentSettings_;
+
+ boost::shared_ptr< OptimizationSettings > optimizationSettings_;
+
+ void optimize( void );
+
+ };
+
+Briefly, the ``missionSegmentSettings_`` object is the vector of mission segment settings to be otpimized. The vector represents a series of **blocks**,
+defined in the following section. ``optimizationSettings_`` contains some useful information for the optimization process, which is performed by
+the ``optimize()`` method, explained in details in the section :ref:`optimizationProcess`.
+
+.. _optimizationStyle:
+
+Optimization style
+~~~~~~~~~~~~~~~~~~
+
+Optimization is performed in blocks. Several ``MissionSegmentSettings`` objects are loaded inside
+``MissionLinker`` in a vector. Because of the current state of the art of the tool, multi-objective optimization
+has not been introduced yet, this means that several objective functions cannot be optimized at the same
+time. The software prevents this issue by optimizing a set of mission segments at a time, as shown in the following
+picture:
+
+ .. _figure-missionSegmentBlock:
+
+ .. figure:: Missionsegmentblock.svg
+ :align: center
+
+ Division of optimization blocks in a chain of mission segments
+
+As you can see, the first mission segment does not contain either a decision variable or an objective function, therefore
+it is not incuded in any block. The software will propagate this mission segment and just keep it out of the optimization process.
+Each of the two blocks starts at the first mission segment containing a decision variable and
+ends at the first mission segment with an objective function.
+Block 1 is optimized independently from Block 2, but before performing the optimization on each block, the starting mission
+segment is linked to the previous in the chain.
+Linking is also performed inside the block, between the contiguous mission segments.
+See the section on the :ref:`optimizationProcess` for information on the code.
+
+
+.. _introPagmoProblem:
+
+Introduction of PaGMO problem
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Of course, there is a PaGMO problem involved: the ``LinkingProblem`` class. You may find it in the file ``missionLinker.cpp``.
+See the page :ref:`tudatFeaturesPagmo2` for more information about PaGMO problems and procedures.
+In a PaGMO fashion, the problem is defined as following:
+
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/missionLinker.cpp`
+ :name: missionLinker-cpp
+
+ struct LinkingProblem{
+
+ LinkingProblem( ){ }
+
+ LinkingProblem( std::vector< boost::shared_ptr< MissionSegmentSettings > >& missionSegments ) :
+ missionSegments_( missionSegments )
+ {
+ nOfDecisionVariables_ = getNumberOfDecisionVariables( missionSegments_ );
+ boundaries_ = getBoundariesPair( missionSegments_ );
+ }
+
+ ~LinkingProblem( ){ }
+
+ std::vector< double > fitness( const std::vector< double > &decisionVariables ) const
+ {
+
+ ...
+
+ }
+
+ std::pair< std::vector< double >, std::vector< double > > get_bounds( ) const
+ {
+ return boundaries_;
+ }
+
+ private:
+
+ std::vector< boost::shared_ptr< MissionSegmentSettings > > missionSegments_;
+ std::pair< std::vector< double >, std::vector< double > > boundaries_;
+ std::vector< int > nOfDecisionVariables_;
+
+ };
+
+
+The ``missionSegments_`` vector, in this case, can contain **one and only one block** of mission segments, as defined in the previous
+section :ref:`optimizationStyle`.
+``nOfDecisionVariables_`` is a vector of integer containing the total number of scalar decision variables in each mission segment in the block,
+while ``boundaries_`` are, of course, the pair of boundary vectors in PaGMO style retrieved from the ``decisionVariableSettings_`` objects in each
+mission segment containing a non-``NULL`` ``decisionVariableSettings_`` member.
+The functions ``getNumberOfDecisionVariables()`` and ``getBoundariesPair()`` have been written and employed to extrapolate the information
+from the block of mission segments. You may find them in the file ``missionLinker.cpp``.
+
+The ``fitness()`` method has, obviously, the purpose of retrieving the value of the objective function at the end of the block.
+The code is separately shown below:
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/missionLinker.cpp`
+ :name: missionLinker-cpp
+
+ ...
+
+ std::vector< double > fitness( const std::vector< double > &decisionVariables ) const
+ {
+ int accumulator = 0;
+ Eigen::VectorXd newValues;
+
+ for( unsigned int i = 0; i < missionSegments_.size(); i++ ){
+
+ // Set the (interpolated) final states of the previous mission segment
+ if( i > 0 ){
+
+ // get the interpolated end of simulation's epoch
+ double prevFinalTime = missionSegments_[i-1]->getFinalSimulationEpoch();
+
+ // get the interpolated final state of the simulation's epoch
+ std::map< std::string, Eigen::Vector6d > prevFinalStates = missionSegments_[i-1]->
+ getFinalOrbitalStates( prevFinalTime );
+
+ std::map< std::string, double > prevFinalMasses = std::map< std::string, double >();
+
+ if( missionSegments_[i-1]->hasPropagatorSettings( propagators::body_mass_state ) )
+ prevFinalMasses = missionSegments_[i-1]->getFinalBodyMasses( prevFinalTime );
+
+ // reset the initial time and states of the dynamics simulator
+ missionSegments_[i]->setPreliminaryConditions( prevFinalTime, prevFinalStates,
+ prevFinalMasses );
+
+ }
+
+ if( nOfDecisionVariables_[i] > 0 )
+ {
+ newValues.resize( nOfDecisionVariables_[i] );
+
+ for( int j = 0; j < nOfDecisionVariables_[i]; j++)
+ {
+ newValues[j] = decisionVariables[accumulator];
+ accumulator++;
+ }
+
+ missionSegments_[i]->setPropagationConditions( newValues );
+ }
+
+ missionSegments_[i]->dynamicsSimulator_->integrateEquationsOfMotion(
+ boost::dynamic_pointer_cast< propagators::TranslationalStatePropagatorSettings< > >(
+ missionSegments_[i]->dynamicsSimulator_->getPropagatorSettings() )->getInitialStates() );
+
+ }
+
+ ...
+
+The following figure shows an example on how the problem works:
+
+.. _figure-fitnessExample:
+
+ .. figure:: MissionLinker.svg
+ :align: center
+
+ Visual example of fitness() function in the case of a block of 3 mission segments with
+ different numbers of decision variables.
+
+
+Basically, PaGMO needs to optimize the whole block in each iteration. This means that all the decision variables of all the mission segments in the
+block are to be passed to the ``fitness()`` function. The partition of these values to the different mission segments must then be performed inside ``fitness()``.
+In this example a block with three mission segments, in which the first contains 3 decision variables, the second none and the third one.
+
+Using the function ``getBoundariesPair()`` all the 4 boundaries are packed together into the ``boundaries_`` member. PaGMO uses these boundaries
+to generate values for the decision variables. In this case 4 decision variable values are generated by PaGMO and passed as vector
+``decisionVariables`` into ``fitness()``. The fitness uses the information stored in ``nOfDecisionVariables_`` to decide to which
+mission segments these values belong.
+
+After the propagation of the first mission segment, **linkage** to the next mission segment is performed: this involves retrieving the final
+simulation epoch, states and masses of the propagated objects of the current mission segment and passing them to the next as initial conditions.
+``MissionSegmentSettings::`` built-in methods ``getFinalSimulationEpoch()``, ``getFinalOrbitalstates()``, ``getFinalBodyMasses()`` and
+``setPreliminaryConditions()`` come in handy here.
+
+The PaGMO generated values of the decision variables are then passed to the next mission segment. The code scouts all the items in the ``decisionVariables`` vector and
+groups them in ``Eigen::VectorXd`` objects called ``newValues``, sized accordingly to the number of scalar decision variables in the segment.
+The propagation is then set by using the method ``MissionSegmentSettings::setPropagationConditions``, which uses the information provided
+by the ``decisionVariableSettings_`` member to modify the mission segment's ``dynamicsSimulator_`` accordingly.
+
+
+.. _optimizationProcess:
+
+Optimization process
+~~~~~~~~~~~~~~~~~~~~
+
+Optimization is performed by the method ``optimize()``:
+
+ .. _figure-optimizeFluxDiagram:
+
+ .. figure:: fluxDiagram.svg
+ :align: center
+
+ Flux diagram of the optimize() method
+
+Most of the above alghorithm works to ensure the :ref:`figure-missionSegmentBlock` defined in the section on the :ref:`optimizationStyle`.
+
+The focus of a potential developer, as well as the heart of the process, is the **Optimize block** block on the rightmost diagram branch.
+In here the call for PaGMO procedures is performed.
+Let us dive into the code, found in the file ``missionLinker.cpp``, in the ``MissionLinker::optimize()`` method declaration:
+
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/missionLinker.cpp`
+ :name: missionLinker-cpp
+
+ ...
+
+ if( start_series == 2 )
+ {
+
+ // PAGMO TAKEOVER
+
+ // Create LinkingProblem
+ pagmo::problem optimizationProblem = pagmo::problem{
+ LinkingProblem( missionSegmentsBlock ) };
+ pagmo::algorithm algo;
+
+ // Set the right algorithm according to the opimization settings
+ if( optimizationSettings_->optimizationType_ == global_optimization )
+ algo = pagmo::algorithm{ pagmo::de1220() };
+ else
+ algo = pagmo::algorithm{ pagmo::compass_search() };
+
+ // Set population according to the optimization settings
+ pagmo::population::size_type populationSize = optimizationSettings_->populationSize_;
+ pagmo::island isl = pagmo::island{ algo, optimizationProblem,
+ populationSize };
+
+ unsigned int counter = 0;
+
+ //Start evolving
+ while( true )
+ {
+ counter ++;
+ isl.evolve();
+ for( ; isl.status()!=pagmo::evolve_status::idle; )
+ isl.wait();
+
+ if( optimizationSettings_->verbosity_ > 0 )
+ {
+ // Show results according to verbosity
+ if( counter % optimizationSettings_->verbosity_ == 0)
+ {
+ std::cout << "Evolution n: " << counter << " Objective Function: " << isl.get_population().champion_f()[0] << "\n";
+ fflush(stdout);
+ }
+ }
+ // Stop if tolerance has been reached
+ if( missionSegmentsBlock.back()->objectiveFunctionSettings_->objectiveValueIsSet_ )
+ {
+ if( isl.get_population().champion_f()[0] <= missionSegmentsBlock.back()
+ ->objectiveFunctionSettings_->tolerance_ )
+ break;
+ }
+ // Stop if max number of evolutions has been reached
+ if( optimizationSettings_->stopAtMaxNumberOfEvolutions_ && ( counter == missionSegmentsBlock.back()->objectiveFunctionSettings_->
+ maxNumberOfEvolutions_ ) )
+ break;
+ }
+
+ // END OF PAGMO TAKEOVER
+
+
+ // Perform the propagation with the optimum values
+ unsigned int accumulator = 0;
+ std::vector< double > optimizedVariables = isl.get_population().champion_x();
+ std::vector< int > nOfDecisionVariables = getNumberOfDecisionVariables( missionSegmentsBlock );
+
+ for( unsigned int j = 0; j < missionSegmentsBlock.size(); j++ )
+ {
+ if( nOfDecisionVariables[j] > 0 )
+ {
+
+ Eigen::VectorXd newValues;
+ newValues.resize( nOfDecisionVariables[j] );
+
+ // Partition the optimum values for each missionSegment
+ for( int k = 0; k < nOfDecisionVariables[j]; k++ )
+ {
+ newValues[k] = optimizedVariables[accumulator];
+ accumulator++;
+ }
+
+ // Set the optimum values as propagation conditions
+ missionSegmentsBlock[j]->decisionVariableValues_ = newValues;
+ missionSegmentsBlock[j]->setPropagationConditions(newValues);
+
+ // Integrate the mission segment
+ missionSegmentsBlock[j]->dynamicsSimulator_->integrateEquationsOfMotion(
+ boost::dynamic_pointer_cast<
+ propagators::TranslationalStatePropagatorSettings< > >(
+ missionSegmentsBlock[j]->selectPropagatorSettings(
+ propagators::transational_state ) )->getInitialStates() );
+
+ }
+
+ // Perform the linkage with the next mission segment
+ if( j < missionSegmentsBlock.size() - 1)
+ {
+ // Get the final epoch
+ double prevFinalTime = missionSegmentsBlock[j]->getFinalSimulationEpoch();
+
+ // Get the final states
+ std::map< std::string, Eigen::Vector6d > prevFinalStates = missionSegmentsBlock[j]->
+ getFinalOrbitalStates( prevFinalTime );
+
+ std::map< std::string, double > prevFinalMasses = std::map< std::string, double >();
+
+ // Get the final masses
+ if( missionSegmentsBlock[j]->hasPropagatorSettings( propagators::body_mass_state ) )
+ prevFinalMasses = missionSegmentsBlock[j]->getFinalBodyMasses( prevFinalTime );
+
+ // Reset the initial time and states of the dynamics simulator
+ missionSegmentsBlock[j]->setPreliminaryConditions( prevFinalTime, prevFinalStates,
+ prevFinalMasses );
+ }
+ }
+
+ // End of block
+ start_series = 0;
+ }
+
+ ...
+
+Although it seems a lot is happening here, there is no need to fear. The code starts right after the non-``NULL`` ``objectiveFunction_`` object has
+been found inside a mission segment, which triggers the flag ``start_series`` to be set to 2. You can find the rest of the code in
+``missionLinker.cpp``.
+
+The code inside the **PAGMO TAKEOVER** tags is where the optimization occurs. The ``LinkingProblem`` described in the section on the :ref:`introPagmoProblem`
+is here constructed, passing the block of mission segments, the ``std::vector< boost::share_ptr< MissionSegmentSettings > >`` type object
+``missionSegmentsBlock``, which has been created in the process described in the :ref:`figure-optimizeFluxDiagram`.
+The ``optimizationSettings->optimizationType_`` object contains the type of optimization chosen by the user:
+
+ * ``global_optimization`` calls the PaGMO alghorithm ``de1220()``, which is a kind of differential evolutionary alghorithm;
+ * ``local_optimization`` (implicitly checked in the ``else`` statement) calls the PaGMO alghorithm ``compass_search()``.
+
+At the moment these are the only two types of optimization integrated in Tudat. As you can see, adding more optimization options can be as
+simple as adding an if-case and an enumeration to the ``enum OptimizationTypes`` in the file ``missionLinker.h``.
+
+Despite the fact that the propagation would automatically update the dynamics simulator, this does not happen in this case:
+**PaGMO** creates copies of all the objects involved and performs the propagation in separate threads. As a result, the mission
+segments in our ``missionSegmentBlock`` do not get modified! That is why a final propagation and linking must be performed, which is the purpose
+of the code outside the **PAGMO TAKEOVER** tags.
+
+
diff --git a/developerGuide/extendingOptimization/missionSegment.rst b/developerGuide/extendingOptimization/missionSegment.rst
new file mode 100644
index 0000000000..bb000f3a56
--- /dev/null
+++ b/developerGuide/extendingOptimization/missionSegment.rst
@@ -0,0 +1,179 @@
+.. _extendMissionSegment:
+
+Mission Segment Settings
+========================
+
+What is it
+~~~~~~~~~~
+
+A **mission segment** in this framework is defined by a ``SingleArcDynamicsSimulator``,
+a ``DecisionVariableSettings`` object which establishes what variables in the dynamics simulator
+have to be tuned in order to optimize the mission and an ``ObjectiveFunctionSettings`` (or derived classes)
+object used to define the function to optimize.
+
+The class is ``MissionSegmentSettings``. See the file ``missionSegmentSettings.h`` to visualize the code.
+
+Linkage aiding methods
+~~~~~~~~~~~~~~~~~~~~~~
+
+Most of the methods used to link together two consecutive mission segments can be found in the
+``MissionSegmentSettings`` class. Their codes can be found in the file ``missionSegmentSettings.cpp``.
+
+They all use :ref:`optimization-linearInterpolation` to retrieve the final states.
+
+Their implementation is pretty straightforward. I will leave the developer to explore the code by their own.
+
+Decision variable settings
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Most of the methods to modify decision variables are also found in ``missionSegmentSettings.cpp``.
+Their implementation is also pretty straightforward. The developer must have a previous knowledge of Tudat and
+:ref:`optimization-dynamicPointerCast`.
+
+
+.. _getObjectiveFunctionValue:
+
+The method getObjectiveFunctionValue()
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The retrieval of the value is performed in the class ``MissionSegmentSettings`` using the function
+
+ .. code-block:: cpp
+
+ double getObjectiveFunctionValue( void )
+
+ .. warning:: The objective function based on final integration values must always be interpolated to the final conditions imposed by the user in the termination settings.
+
+Let us take into consideration the case where the objective function is of the type ``FinalSphericalOrbitalElementObjectiveFunctionSettings``.
+In this case the code is going to retrieve, after the propagation, the final spherical orbital element defined by the user in the objective function settings
+in the state map of the desired body.
+Here's the code:
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/missionSegmentSettings.cpp`
+ :name: missionSegmentSettings-cpp
+
+ double MissionSegmentSettings::getObjectiveFunctionValue( void )
+ {
+
+ double finalTime = getFinalSimulationEpoch();
+
+ double objectiveFunctionValue;
+
+ ...
+
+ else if( objectiveFunctionSettings_->objectiveVariable_ == final_spherical_orbital_element_objective_function )
+ {
+ boost::shared_ptr< FinalSphericalOrbitalElementObjectiveFunctionSettings > objectiveFunctionSettings =
+ boost::dynamic_pointer_cast< FinalSphericalOrbitalElementObjectiveFunctionSettings >(
+ objectiveFunctionSettings_);
+ std::map< std::string, Eigen::Vector6d > finalStates = getFinalOrbitalStates( finalTime );
+ Eigen::Vector6d sphericalState = orbital_element_conversions::convertCartesianToSphericalOrbitalState(
+ finalStates[objectiveFunctionSettings->associatedBody_] );
+ objectiveFunctionValue = sphericalState[objectiveFunctionSettings->sphericalOrbitalElement_];
+ }
+
+ ...
+
+ }
+
+First of all notice that the method ``getFinalSimulationEpoch()`` has been used to retrieve the variable finalTime. This is a method of the ``MissionSegmentSettings`` class,
+which calculates the final epoch according to the user defined termination settings. Please visit the relative documentation for more information.
+The member ``objectiveFunctionSettings_`` in ``MissionSegmentSettings`` is stored as a ``boost::shared_ptr< ObjectiveFunctionSettings >`` object.
+In order to retrieve it as a ``FinalSphericalOrbitalElementObjectiveFunctionSettings``, which is a derived class, a dynamic pointer cast of the
+member is performed, after constatating that
+
+ .. code-block:: cpp
+
+ objectiveFunctionSettings_->objectiveVariable_ == final_spherical_orbital_element_objective_function
+
+
+The rest of the code is pretty straightforward. The final orbital state is recovered with the ``MissionSegmentSettings``method
+``getFinalOrbitalStates( double )`` which accepts the epoch to which the states are interpolated.
+
+ .. code-block:: cpp
+
+ std::map< std::string, Eigen::Vector6d > finalStates = getFinalOrbitalStates( finalTime );
+
+conveniently the method maps the orbital states in a map with the name of the body as key.
+
+All the methods shown above to retrieve the final states are based on linear interpolation. In some cases there is not a method directly
+defined, in any case the final result must always be interpolated to the user defined final conditions. For example, for the case of
+``ObjectiveFunctionFromFinalDependentVariableSettings``, the code is the following:
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/missionSegmentSettings.cpp`
+ :name: missionSegmentSettings-cpp
+
+
+ double MissionSegmentSettings::getObjectiveFunctionValue( void )
+ {
+
+ double finalTime = getFinalSimulationEpoch();
+
+ double objectiveFunctionValue;
+ ...
+
+ else if( objectiveFunctionSettings_->objectiveVariable_ == final_dependent_variable_objective_function )
+ {
+ boost::shared_ptr< ObjectiveFunctionFromFinalDependentVariableSettings > objectiveFunctionSettings =
+ boost::dynamic_pointer_cast< ObjectiveFunctionFromFinalDependentVariableSettings >(
+ objectiveFunctionSettings_ );
+ int positionInVector = objectiveFunctionSettings->getPositionInDependentVariablesSaveVector(
+ dynamicsSimulator_ );
+ std::map< double, Eigen::VectorXd > dependentVariableHistory
+ = dynamicsSimulator_->getDependentVariableHistory();
+ std::map< double, Eigen::VectorXd >::reverse_iterator it = dependentVariableHistory.rbegin();
+
+ double time1 = it->first;
+ double value1 = it->second[positionInVector];
+
+ it++;
+
+ double time2 = it->first;
+ double value2 = it->second[positionInVector];
+
+ // If the last two epochs do not contain the final time in their boundaries
+ // scout all the epochs
+ while( !( time2 <= finalTime && time1 >= finalTime) )
+ {
+ time1 = time2;
+ value1 = value2;
+
+ it++;
+
+ time2 = it->first;
+ value2 = it->second[positionInVector];
+
+ }
+
+ objectiveFunctionValue = value2 + ( value1 - value2 )*(finalTime - time2)/(time1 - time2);
+ }
+
+ ...
+ }
+
+The method ``ObjectiveFunctionFromFinalDependentVariableSettings::getPositionInDependentVariablesSaveVector()`` is used in this case
+to retrieve the position of the dependent variable settings defined in the objective function object from the dynamics simulator.
+
+The method ends with the lines:
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/missionSegmentSettings.cpp`
+ :name: missionSegmentSettings-cpp
+
+ ...
+
+ if( objectiveFunctionSettings_->objectiveValueIsSet_ == true )
+ return fabs( objectiveFunctionSettings_->objectiveValue_ - objectiveFunctionValue );
+ else if( objectiveFunctionSettings_->minimizeOrMaximize_ )
+ return objectiveFunctionValue;
+ else
+ return -objectiveFunctionValue;
+
+ ...
+
+The first if-case is for an objective function settings defined with the first constructor, in which case the return value is the
+absolute difference of the retrieved value and the objective value. The other two cases are for an objective function settings defined
+with the second constructor: for the second if-case (minimize) simply the retrieved value is returned, otherwise (maximize) the
+negative retrieved values is retrieved. That's because PaGMO is set to minimize the function.
diff --git a/developerGuide/extendingOptimization/necessaryFeatures.rst b/developerGuide/extendingOptimization/necessaryFeatures.rst
new file mode 100644
index 0000000000..5c0159109f
--- /dev/null
+++ b/developerGuide/extendingOptimization/necessaryFeatures.rst
@@ -0,0 +1,55 @@
+.. _optimizationNecessaryFeatures:
+
+Ideas for future development
+============================
+
+At the moment several features are missing from the optimization package, which should be included to
+provide the best Tudat experience for the users.
+
+Linkage between different reference frames
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The linking of two mission segments in the :ref:`extendMissionLinker` does not take care to recognize the two reference frames are
+centered at the same body. The program should recognize whether the two reference frames are the same, and, if not,
+perform the due translations, maybe aided by the ``ReferenceFrameManager`` class. This part should be implemented in the
+``optimize()`` method of the ``MissionLinker`` class. See file ``Tudat/Optimization/missionLinker.cpp``.
+Beware that at the moment **Tudat do not support different frame orientations**, therefore only the centre of propagation can be
+modified.
+
+Centre of frame for orbital element decision variables
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The classes ``SingleCartesianComponentDecisionVariableSettings``, ``SingleKeplerElementDecisionVariableSettings`` and
+``SingleSphericalOrbitalElementDecisionVariableSettings``, as well as ``SingleDecisionVariableSettings`` using the flags
+``initial_cartesian_state_decision_variable``, ``initial_cartesian_velocity_decision_variable`` and
+``initial_cartesian_position_decision_variable`` allow to use the initial orbital states of a body as decision variables.
+
+At the moment the centre of the reference frame around which these orbital states are defined is the current ``centralBodies`` vector
+associated index in the propagator settings, although the user should be able to define its own reference frame centre.
+
+These modifications should be implemented both in the structures of the above mentioned classes, in the files
+``Tudat/Optimization/decisionVariableSettings.h`` and ``Tudat/Optimization/decisionVariableSettings.cpp`` and in the methods
+``MissionLinker::getInitialOrbitalState( )``, ``MissionLinker::resetOrbitalState( )``, ``MissionLinker::modifyCartesianComponent( )``,
+``MissionLinker::modifyKeplerElement( )`` and ``MissionLinker::modifySphericalComponent( )``, found in the file ``Tudat/Optimization/missionSegmentSettings.cpp``.
+Also here the challenges are mainly about the frames in which each body is centered. The developer can use the ``ReferenceFrameManager`` class for this purpose.
+
+Pareto front multi-objective optimization
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Despite PaGMO allows some methods to create a multi-objective optimization problem, at the moment this option is not yet introduced in
+the optimization namespace. The two options allowed in the ``OptimizationSettings`` class are ``local_optimization`` and ``global_optimization``,
+each referencing to only one optimization alghorithm provided by PaGMO.
+
+
+Other optimization settings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+At the moment the ``OptimizationSettings`` class, defined in the file ``Tudat/Optimization/missionLinker.h`` allows for a sparse number of optimization
+settings. Some more options should be introduced to allow the user a more flexible optimization capability.
+
+
+Miscellaneous
+~~~~~~~~~~~~~
+
+ * An exact Lambert targeter, starting at a certain orbit around a departure Planet and ending at a certain orbit around a target Planet. This can only be implemented with an optimization process.
+ * Other?
diff --git a/developerGuide/extendingOptimization/necessaryFeatures.srt b/developerGuide/extendingOptimization/necessaryFeatures.srt
new file mode 100644
index 0000000000..5c0159109f
--- /dev/null
+++ b/developerGuide/extendingOptimization/necessaryFeatures.srt
@@ -0,0 +1,55 @@
+.. _optimizationNecessaryFeatures:
+
+Ideas for future development
+============================
+
+At the moment several features are missing from the optimization package, which should be included to
+provide the best Tudat experience for the users.
+
+Linkage between different reference frames
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The linking of two mission segments in the :ref:`extendMissionLinker` does not take care to recognize the two reference frames are
+centered at the same body. The program should recognize whether the two reference frames are the same, and, if not,
+perform the due translations, maybe aided by the ``ReferenceFrameManager`` class. This part should be implemented in the
+``optimize()`` method of the ``MissionLinker`` class. See file ``Tudat/Optimization/missionLinker.cpp``.
+Beware that at the moment **Tudat do not support different frame orientations**, therefore only the centre of propagation can be
+modified.
+
+Centre of frame for orbital element decision variables
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The classes ``SingleCartesianComponentDecisionVariableSettings``, ``SingleKeplerElementDecisionVariableSettings`` and
+``SingleSphericalOrbitalElementDecisionVariableSettings``, as well as ``SingleDecisionVariableSettings`` using the flags
+``initial_cartesian_state_decision_variable``, ``initial_cartesian_velocity_decision_variable`` and
+``initial_cartesian_position_decision_variable`` allow to use the initial orbital states of a body as decision variables.
+
+At the moment the centre of the reference frame around which these orbital states are defined is the current ``centralBodies`` vector
+associated index in the propagator settings, although the user should be able to define its own reference frame centre.
+
+These modifications should be implemented both in the structures of the above mentioned classes, in the files
+``Tudat/Optimization/decisionVariableSettings.h`` and ``Tudat/Optimization/decisionVariableSettings.cpp`` and in the methods
+``MissionLinker::getInitialOrbitalState( )``, ``MissionLinker::resetOrbitalState( )``, ``MissionLinker::modifyCartesianComponent( )``,
+``MissionLinker::modifyKeplerElement( )`` and ``MissionLinker::modifySphericalComponent( )``, found in the file ``Tudat/Optimization/missionSegmentSettings.cpp``.
+Also here the challenges are mainly about the frames in which each body is centered. The developer can use the ``ReferenceFrameManager`` class for this purpose.
+
+Pareto front multi-objective optimization
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Despite PaGMO allows some methods to create a multi-objective optimization problem, at the moment this option is not yet introduced in
+the optimization namespace. The two options allowed in the ``OptimizationSettings`` class are ``local_optimization`` and ``global_optimization``,
+each referencing to only one optimization alghorithm provided by PaGMO.
+
+
+Other optimization settings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+At the moment the ``OptimizationSettings`` class, defined in the file ``Tudat/Optimization/missionLinker.h`` allows for a sparse number of optimization
+settings. Some more options should be introduced to allow the user a more flexible optimization capability.
+
+
+Miscellaneous
+~~~~~~~~~~~~~
+
+ * An exact Lambert targeter, starting at a certain orbit around a departure Planet and ending at a certain orbit around a target Planet. This can only be implemented with an optimization process.
+ * Other?
diff --git a/developerGuide/extendingOptimization/objectiveFunction.rst b/developerGuide/extendingOptimization/objectiveFunction.rst
new file mode 100644
index 0000000000..cfa682c303
--- /dev/null
+++ b/developerGuide/extendingOptimization/objectiveFunction.rst
@@ -0,0 +1,88 @@
+.. _extendObjectiveFunction:
+
+Objective Function Settings
+===========================
+
+How does it work
+~~~~~~~~~~~~~~~~
+
+The class ``ObjectiveFunctionSettings`` and its derived are used to retrieve some value from the the last
+mission segment after propagation.
+
+The structure of the base class for the objective function segments is the following:
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/objectiveFunction.h`
+ :name: objectiveFunction-h
+
+ struct ObjectiveFunctionSettings{
+
+ ObjectiveFunctionSettings( ObjectiveFunctionType objectiveVariable,
+ const double objectiveValue, const double tolerance = 1e-3,
+ const unsigned int maxNumberOfEvolutions = 100 ) :
+ objectiveVariable_( objectiveVariable ), objectiveValue_( objectiveValue ), tolerance_( tolerance ),
+ maxNumberOfEvolutions_( maxNumberOfEvolutions ), objectiveValueIsSet_( true ){ }
+
+ ObjectiveFunctionSettings( ObjectiveFunctionType objectiveVariable, const bool minimizeOrMaximize = true,
+ const unsigned int maxNumberOfEvolutions = 100 ) :
+ objectiveVariable_(objectiveVariable), minimizeOrMaximize_( minimizeOrMaximize ),
+ maxNumberOfEvolutions_( maxNumberOfEvolutions ), objectiveValueIsSet_( false ){ }
+
+ virtual ~ObjectiveFunctionSettings( ){ }
+
+ ObjectiveFunctionType objectiveVariable_;
+ double objectiveValue_;
+ double tolerance_;
+ bool minimizeOrMaximize_;
+ unsigned int maxNumberOfEvolutions_;
+ bool objectiveValueIsSet_;
+
+ };
+
+The two constructor are used in two different cases: either an objective value is set and in that case the objective function will
+converge to that value, or the objective value will just be minimized or maximized accordingly to the ``minimizeOrMaximize`` variable
+(``true`` for minimize and ``false`` for maximize).
+
+Types of objective functions currently available
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The current types of objective function are reflected in the enumerator ``ObjectiveFunctionType``.
+
+ .. code-block:: cpp
+ :caption: :class:`Tudat/Optimization/objectiveFunction.h`
+ :name: objectiveFunction-h
+
+ enum ObjectiveFunctionType
+ {
+ final_cartesian_component_objective_function = 0,
+ final_kepler_orbital_element_objective_function = 1,
+ final_spherical_orbital_element_objective_function = 2,
+ final_dependent_variable_objective_function = 3,
+ min_max_dependent_variable_objective_function = 4,
+ user_defined_objective_function = 5
+ };
+
+These enumerations are used in the class ``ObjectiveFunctionSettings`` and its derived to define the member
+``ObjectiveFunctionSettings::objectiveVariable_``.
+This member is not directly defined by the Tudat user, as each enumerator belong to different derived classes with different settings, instead, when creating a
+the derived class the member ``objectiveVariable_`` is set automatically. Here's a
+definition of the above defined enumerations:
+
+ * ``final_cartesian_component_objective_function`` is the value of ``objectiveVariable_`` when an object of class ``FinalCartesianComponentObjectiveFunctionSettings`` is created;
+ * ``final_kepler_orbital_element_objective_function`` is the value of ``objectiveVariable_`` when an object of class ``FinalKeplerOrbitalElementObjectiveFunction`` is created;
+ * ``final_spherical_orbital_element_objective_function`` is the value of ``objectiveVariable_`` when an object of class ``FinalSphericalObitalElementObjectiveFunction`` is created;
+ * ``final_dependent_variable_objective_function`` is the value of ``objectiveVariable_`` when an object of class ``FinalDependentVariableObjectiveFunction``;
+ * ``min_max_dependent_variable_objective_function`` is the value of ``objectiveVariable_`` when an object of class ``MinMaxDependentVariableObjectiveFunction``;
+ * ``user_defined_objective_function`` is the value of ``objectiveVariable_`` when an object of class ``UserDefinedObjectiveFunction`` is created.
+
+You are invited to explore the file ``objectiveFunction.h`` to understand the differences of the various derived classes
+
+Retrieving the objective function value
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+See the section :ref:`getObjectiveFunctionValue` for an understading on how the objective function value is retrieved from the ``dynamicsSimulator_`` member of the
+``MissionSegmentSettings`` class.
+
+
+
+
diff --git a/developerGuide/index.rst b/developerGuide/index.rst
index ce28d7acbb..30d5b4bb8c 100644
--- a/developerGuide/index.rst
+++ b/developerGuide/index.rst
@@ -11,6 +11,7 @@ In these pages of the wiki you will find a number of tutorials that will help yo
githubBasics
extendingJSON/index
extendingMATLAB/index
+ extendingOptimization/index
guideToDebugging
cppTutorials
howToWriteTheWiki
diff --git a/tutorials/applicationWalkthroughs/images/RendezVousXY.svg b/tutorials/applicationWalkthroughs/images/RendezVousXY.svg
new file mode 100644
index 0000000000..9af3c728a8
--- /dev/null
+++ b/tutorials/applicationWalkthroughs/images/RendezVousXY.svg
@@ -0,0 +1,268 @@
+
+
+
diff --git a/tutorials/applicationWalkthroughs/images/RendezVousXZ.svg b/tutorials/applicationWalkthroughs/images/RendezVousXZ.svg
new file mode 100644
index 0000000000..c5f2f77271
--- /dev/null
+++ b/tutorials/applicationWalkthroughs/images/RendezVousXZ.svg
@@ -0,0 +1,200 @@
+
+
+
diff --git a/tutorials/applicationWalkthroughs/index.rst b/tutorials/applicationWalkthroughs/index.rst
index ec9787d288..7a74b82fe3 100644
--- a/tutorials/applicationWalkthroughs/index.rst
+++ b/tutorials/applicationWalkthroughs/index.rst
@@ -29,4 +29,6 @@ Each example is discussed in detail on a separate page. Note that in the first e
variationalEquationsPropagation
earthOrbiterBasicStateEstimation
earthOrbiterStateEstimation
+ satelliteRendezVous
+
diff --git a/tutorials/applicationWalkthroughs/satelliteRendezVous.rst b/tutorials/applicationWalkthroughs/satelliteRendezVous.rst
new file mode 100644
index 0000000000..3f4cfb0d3f
--- /dev/null
+++ b/tutorials/applicationWalkthroughs/satelliteRendezVous.rst
@@ -0,0 +1,435 @@
+.. _satelliteRendezVous:
+
+Mission Linking and Optimization
+================================
+
+In this tutorial we introduce the features of the Optimization package, including the MissionLinker class and how to set up an optimization problem in Tudat.
+
+ The code for this tutorial is given on GitHub and is also located in the tudatBundle at:
+
+ .. code-block:: cpp
+
+ tudatBundle/tudatExampleApplications/satellitePropagatorExamples/SatellitePropagatorExamples/satelliteRendezVousExample.cpp
+
+The code follows the following example: two satellites are approximately on the same orbit around Earth. The two are separated by a small angle, with the Chaser
+falling behind the Target. We need to optimize the initial burn time for the chaser in order to rendez-vous with the target.
+
+These are the initial orbital values for the two satellites:
+
+ .. tabularColumns:: l R R
+
+ +-------------------------------------------+------------+------------+
+ | **Parameter** | **Target** | **Chaser** |
+ +-------------------------------------------+------------+------------+
+ | Semi-major axis :math:`a` [m] | 6782650.0 | 6782650.0 |
+ +-------------------------------------------+------------+------------+
+ | Eccentricity :math:`e` | 0.00050 | 0.00025 |
+ +-------------------------------------------+------------+------------+
+ | Inclination :math:`i` [°] | 51.640 | 51.635 |
+ +-------------------------------------------+------------+------------+
+ | Argument of perigee :math:`\omega` [°] | 235.700 | 84.100 |
+ +-------------------------------------------+------------+------------+
+ | RAAN :math:`\Omega` [°] | 23.400 | 23.400 |
+ +-------------------------------------------+------------+------------+
+ | Initial true anomaly :math:`\theta_0` [°] | 139.800 | 290.500 |
+ +-------------------------------------------+------------+------------+
+
+
+From the table one can discern that, since the two spacecraft have approximately the same inclination and exactly the same RAAN, the angular separation between the two is
+
+ .. math::
+
+ \Delta u_0 = \omega_T + (\theta_0)_T - (\omega_C + (\theta_0)_C) = 0.9 [^\circ]
+
+which makes up for an arc distance of
+
+ .. math::
+
+ \Delta y_0 = a \cdot \Delta \theta_0 \approx 106.604 [\text{km}]
+
+In order to rendez-vous, the Chaser needs to gain velocity w.r.t. the Target. Because of the geometry of the problem, the needed impulse and velocity direction can be approximated
+using the Clohessy Whiltshire equations. The approximation sees the Chaser on the line of the velocity vector, in delay w.r.t. the Target:
+
+ .. math::
+
+ y_0 = -\Delta y_0
+
+We decide to perform a burn parallel to the direction of motion, i.e., all the parameters except for :math:`y_0` and :math:`\dot y_0` are set to 0:
+
+ .. math::
+
+ \dot y_0 = \dfrac{y_0 n}{6 \pi}
+
+where :math:`n = \sqrt{\mu / a^3 }` with :math:`\mu` being the orbital parameter of Earth. The fact that :math:`y_0 < 0` means that **the required burn is retrograde**.
+
+The total impulse required is
+
+ .. math::
+
+ \Delta V = | \dot y_0 | \approx 6.39207 [\text{m/s}]
+
+The propulsion system of the Chaser is modelled on the Orbital Maneouvering System of the STS. The values are given in the following table:
+
+ .. tabularColumns:: l
+
+ +---------------------------------------------+
+ | Initial mass :math:`M_0 = 90116.4` [kg] |
+ +---------------------------------------------+
+ | Thrust magnitude :math:`T = 26700.0` [N] |
+ +---------------------------------------------+
+ | Specific impulse :math:`I_{sp} = 316.0` [s] |
+ +---------------------------------------------+
+
+The initial mass represents an empty space-shuttle with enough fuel for orbital maneouver, de-orbiting and a small margin.
+
+Using the Tsiolkowsky equation:
+
+ .. math::
+
+ \Delta V = I_{sp} g_0 ln( \Lambda )
+
+and the relation with a constant thrust:
+
+ .. math::
+
+ t_b = \dfrac {I_{sp} g_0 M_0}{T} \left( 1 - \dfrac{1}{\Lambda} \right)
+
+The total burn time can be estimated to be :math:`t_b \approx 21.552` [s].
+
+The objetive of this example is to optimize this value to minimize the final separation between the spacecraft. In order to do so we have to set-up two simulations in Tudat: one with a thrust acceleration
+model for the chaser and one with only the gravitational accelerations due to Earth, to represent the coasting phase. We will use the MissionLinker class
+to link the two simulations and define the decision variable and the objective function.
+
+Creating the first mission segment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+For the first simulation we need to define the acceleration model in order to include a thrust in the opposite direction of motion for the Chaser.
+After defining the accelerations of the target, which only include the spherical armonics of Earth up to degree and order 4, the accelerations of
+the Chaser are defined in the following piece of code:
+
+ .. code-block:: cpp
+
+ // Space shuttle OMS thrust (AJ10-190 engine)
+ double thrustMagnitude = 26700.0; //N
+ double specificImpulse = 316.0; //s
+
+ // Thrust opposite to direction of motion
+ boost::shared_ptr< ThrustDirectionGuidanceSettings > thrustDirectionGuidanceSettings =
+ boost::make_shared< ThrustDirectionFromStateGuidanceSettings >(
+ "Earth", true, true);
+ boost::shared_ptr< ThrustEngineSettings > thrustMagnitudeSettings =
+ boost::make_shared< ConstantThrustEngineSettings >(
+ thrustMagnitude, specificImpulse );
+
+ // Define chaser acceleration model settings.
+ std::map< std::string, std::vector< boost::shared_ptr< AccelerationSettings > > > accelerationsOfChaser_1;
+ accelerationsOfChaser_1[ "Earth" ].push_back( boost::make_shared< SphericalHarmonicAccelerationSettings >( 4, 0 ) );
+ accelerationsOfChaser_1[ "Chaser" ].push_back(
+ boost::make_shared< ThrustAccelerationSettings >( thrustDirectionGuidanceSettings, thrustMagnitudeSettings) );
+ accelerationMap_1[ "Chaser" ] = accelerationsOfChaser_1;
+ bodiesToPropagate_1.push_back( "Chaser" );
+ centralBodies_1.push_back( "Earth" );
+
+ // Create acceleration models and propagation settings.
+ basic_astrodynamics::AccelerationMap accelerationModelMap_1 = createAccelerationModelsMap(
+ bodyMap_1, accelerationMap_1, bodiesToPropagate_1, centralBodies_1 );
+
+
+The initial Cartesian states are translated from the initial orbital parameters in the table.
+
+A mass propagator for the Chaser is included in the propagator settings:
+
+ .. code-block:: cpp
+
+ // Preliminary time termination settings for first simulation
+ boost::shared_ptr< PropagationTimeTerminationSettings > timeTerminationSettings_1 =
+ boost::make_shared< propagators::PropagationTimeTerminationSettings >( estimatedBurnTime );
+
+ // State propagator for first simulation
+ boost::shared_ptr< TranslationalStatePropagatorSettings< double > > translationalPropagatorSettings =
+ boost::make_shared< TranslationalStatePropagatorSettings< double > >
+ ( centralBodies_1, accelerationModelMap_1, bodiesToPropagate_1, systemInitialState, timeTerminationSettings_1 );
+
+ // Create mass rate model for the chaser
+ boost::shared_ptr< MassRateModelSettings > massRateModelSettings =
+ boost::make_shared< FromThrustMassModelSettings >( 1 );
+ std::map< std::string, boost::shared_ptr< basic_astrodynamics::MassRateModel > > massRateModels;
+ massRateModels[ "Chaser" ] = createMassRateModel(
+ "Chaser", massRateModelSettings, bodyMap_1, accelerationModelMap_1 );
+
+ // Create mass propagator settings for the chaser
+ std::vector< std::string > bodiesWithMassToPropagate;
+ bodiesWithMassToPropagate.push_back( "Chaser" );
+
+ Eigen::VectorXd initialBodyMasses = Eigen::VectorXd( 1 );
+ initialBodyMasses( 0 ) = chaserMass;
+
+ boost::shared_ptr< MassPropagatorSettings< double > > massPropagatorSettings =
+ boost::make_shared< MassPropagatorSettings< double > >(
+ bodiesWithMassToPropagate, massRateModels, initialBodyMasses, timeTerminationSettings_1 );
+
+ // Create multi-type propagator settings (state + chaser mass propagator)
+ std::vector< boost::shared_ptr< PropagatorSettings< double > > > propagatorSettingsVector;
+ propagatorSettingsVector.push_back( translationalPropagatorSettings);
+ propagatorSettingsVector.push_back( massPropagatorSettings );
+
+ boost::shared_ptr< MultiTypePropagatorSettings< double > > propagatorSettings_1 =
+ boost::make_shared< MultiTypePropagatorSettings< double > >(
+ propagatorSettingsVector, timeTerminationSettings_1 );
+
+The time termination settings are just preliminary, as we want to define the time of this part of the simulation as the decision variable. Nevertheless,
+they need to be included, as the optimizer does not create one by itself.
+
+Afterwards, the integrator settings and the dynamics simulator are created. The chosen fixed step size is very small, so to augment the focus of the time decision variable:
+
+ .. code-block:: cpp
+
+ const double fixedStepSize = 0.001; //s
+
+ // Create integrator settings for first simulation
+ boost::shared_ptr< IntegratorSettings< > > integratorSettings_1 =
+ boost::make_shared< IntegratorSettings< > >
+ ( rungeKutta4, simulationStartEpoch, fixedStepSize );
+
+ // Create dynamics simulator for first simulation
+ boost::shared_ptr< SingleArcDynamicsSimulator< > > dynamicsSimulator_1
+ = boost::make_shared< SingleArcDynamicsSimulator< > >(
+ bodyMap_1, integratorSettings_1, propagatorSettings_1, false, false, false );
+
+
+We now proceed to create the decision variable. In this case we want to optimize the burn time, i.e. the simulation time.
+The boundaries are set to search for the optimum time in a 10 seconds interval :math:`t_b - 5 \text{ s} < t < t_b + 5 \text{ s}` :
+
+ .. code-block::
+
+ // Set the simulation time as decision variable, with boundaries +-5 seconds from
+ // the estimated burn time.
+ boost::shared_ptr< optimization::SingleDecisionVariableSettings > decisionVariable_1 =
+ boost::make_shared< optimization::SingleDecisionVariableSettings >(
+ optimization::simulation_time_decision_variable,
+ estimatedBurnTime - 5.0, estimatedBurnTime + 5.0 );
+
+Since the start epoch is set to 0, we coud also use the option ``optimization::SingleDecisionVariableFromTerminationSettings``, which includes the time
+termination settings. It would be set as following:
+
+ .. code-block:: cpp
+
+ boost::shared_ptr< optimization::SingleDecisionVariableFromTerminationSettings > decisionVariable_1 =
+ boost::make_shared< optimization::SingleDecisionVariableFromTerminationSettings >(
+ dynamicsSimulator_1, estimatedBurnTime - 5.0, estimatedBurnTime + 5.0 );
+
+In order to link the dynamics simulator and the decision variable, we use the class MissionSegment:
+
+ .. code-block:: cpp
+
+ // Create first mission segment with the dynamics simulator and the decision variable
+ boost::shared_ptr< optimization::MissionSegmentSettings > missionSegment_1 =
+ boost::make_shared< optimization::MissionSegmentSettings >( dynamicsSimulator_1,
+ decisionVariable_1 );
+
+Creating the second mission segment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Now we proceed to create the second simulation.
+
+Both the Target and the Chaser have as acceleration the spherical armonics gravitational field of Earth up to order and degree 4.
+
+Since we want to minimize the separation between the two spacecraft we are going to need to save the distance between the two bodies as
+a dependent variable:
+
+
+ .. code-block:: cpp
+
+ // Create dependent variable save settings to retrieve
+ // the distance between the two spacecraft
+ std::vector< boost::shared_ptr< SingleDependentVariableSaveSettings > > dependentVariableSaveSettingsVector;
+
+ dependentVariableSaveSettingsVector.push_back( boost::make_shared< SingleDependentVariableSaveSettings > (
+ relative_distance_dependent_variable, "Chaser", "Target" ) );
+
+ boost::shared_ptr< DependentVariableSaveSettings > dependentVariableSaveSettings =
+ boost::make_shared< DependentVariableSaveSettings >( dependentVariableSaveSettingsVector, false );
+
+
+We also want to propagate the simulation until the separation of the two spacecraft is at least 1000 m. We create, therefore, hibrid termination settings, in order to
+include a stop at the minimum required separation or at an adeguate amount of time, here defined as the estimated rendez-vous time plus a small margin:
+
+ .. code-block:: cpp
+
+ // Create hybrid termination settings to simulate until estimated
+ // rendez vous + margin or until minimum separation
+
+ const double minimumSeparation = 1000; //m
+
+ std::vector< boost::shared_ptr< propagators::PropagationTerminationSettings > > multiTerminationSettings;
+
+ multiTerminationSettings.push_back( boost::make_shared< propagators::PropagationTimeTerminationSettings >(
+ estimatedBurnTime + estimatedRendezVousTime + 10*fixedStepSize) );
+
+ multiTerminationSettings.push_back( boost::make_shared<
+ propagators::PropagationDependentVariableTerminationSettings >(
+ dependentVariableSaveSettingsVector[0], minimumSeparation, true ) );
+
+ boost::shared_ptr< propagators::PropagationHybridTerminationSettings > hybridTerminationSettings =
+ boost::make_shared< propagators::PropagationHybridTerminationSettings >( multiTerminationSettings, true );
+
+The estimated rendez-vous time is calculated using the Clohessy-Wiltshire equations as: :math:`t_r = 2 \pi / n`, i.e. the orbital period of the target.
+
+We then proceed to finalize the simulation with the propagator settings, integrator settings and the dynamics simulator:
+
+ .. code-block:: cpp
+
+ // Create propagator settings for second simulation
+ boost::shared_ptr< TranslationalStatePropagatorSettings< double > > propagatorSettings_2 =
+ boost::make_shared< TranslationalStatePropagatorSettings< double > >
+ ( centralBodies_2, accelerationModelMap_2, bodiesToPropagate_2, Eigen::VectorXd::Zero(12),
+ hybridTerminationSettings, cowell, dependentVariableSaveSettings );
+
+ fixedStepSize = 5.0; //s
+
+ // Create integrator settings for second simulation
+ boost::shared_ptr< IntegratorSettings< > > integratorSettings_2 =
+ boost::make_shared< IntegratorSettings< > >( rungeKutta4, 0.0, fixedStepSize );
+
+ // Create dynamics simulator for second simulation
+ boost::shared_ptr< SingleArcDynamicsSimulator< > > dynamicsSimulator_2 =
+ boost::make_shared< SingleArcDynamicsSimulator< > >( bodyMap_2, integratorSettings_2,
+ propagatorSettings_2, false );
+
+The initial state is defined as a vector of zeros (``Eigen::VectorXd::Zero(12)``) and the initial time is set to 0, since, during the optimization process, they will be automatically set to the (interpolated) final values of the
+previous mission segment. For the second simulation we are allowed to use a larger step size of 5 s, since the simulation is less sensitive (there is no thrust this time).
+
+It is convenient to set the **objective function** as the minimum separation between the two spacecraft, since the defined simulation time is arbitrary it is unlikely that the distance between spacecraft
+will reach its minimum value at the final entry. We use the :ref:`ObjectiveFunction` class as following:
+
+ .. code-block:: cpp
+
+ // Use the minimum separation between spacecraft as objective function
+
+ const int maxNumberOfEvolutions = 30;
+ const double objectiveValue = 0.0;
+
+ boost::shared_ptr< optimization::ObjectiveFunctionFromMinOrMaxDependentVariableSettings > objectiveFunction =
+ boost::make_shared< optimization::ObjectiveFunctionFromMinOrMaxDependentVariableSettings >(
+ dependentVariableSaveSettingsVector[0], objectiveValue, 0, true,
+ minimumSeparation, maxNumberOfEvolutions );
+
+We used the previously defined ``dependentVariableSaveSettingsVector[0]`` to get the dependent variable save settings of the distance between spacecraft.
+The objective value is 0, so that the objective function will trend to converge to that value, but we set the tolerance to ``minimumSeparation``, so that
+the optimization process stops when the minimum distance between spacecraft is inferior to 1000 m. Finally, we want to have a maximum of 30 evolutions.
+
+We then define the :ref:`MissionSegment` object of this second simulation to link the dynamics simulator with the objective function:
+
+ .. code-block:: cpp
+
+ // Create mission segment for second simulation containing the second
+ // dynamics simulator and the objective function
+ boost::shared_ptr< optimization::MissionSegmentSettings > missionSegment_2 =
+ boost::make_shared< optimization::MissionSegmentSettings >( dynamicsSimulator_2,
+ objectiveFunction
+
+
+Linkage and optimization
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Finally, we link together the two mission segments in a chronological order with the :ref:`MissionLinker` class:
+
+
+ .. code-block:: cpp
+
+ // Create object to link the two mission segments and optimize the mission constraint
+
+ std::vector< boost::shared_ptr< optimization::MissionSegmentSettings > > missionSegments;
+
+ missionSegments.push_back( missionSegment_1 );
+ missionSegments.push_back( missionSegment_2 );
+
+ optimization::MissionLinker missionLinker(missionSegments,
+ boost::make_shared< optimization::OptimizationSettings >( optimization::global_optimization,
+ 32, true, 1 ), false );
+
+As :ref:`OptimizationSettings`, we are going to use a global optimization with a population of 32 individuals, set the flag that
+stops the simulation at the maximum number of evolutions to ``true`` and the verbosity to 1, as to visualize the results of each
+iteration.
+
+You could either set the last flag to ``true`` or use the statement
+
+ .. code-block:: cpp
+
+ // Start optimization
+ missionLinker.optimize();
+
+to optimize the configuration.
+
+Since the decision variable is the simulation time of the first mission segment, in order to retrieve the optimum,
+which is set automatically at the end of each simulation, you can directly scout the termination settings of the first dynamics simulator:
+
+ .. code-block:: cpp
+
+ std::cout << "\n Calculated burn time = " << timeTerminationSettings_1->terminationTime_ << "s\n";
+
+The same is done with the optimized results of the two simulations:
+
+ .. code-block:: cpp
+
+ std::map< double, Eigen::VectorXd > integrationResult_1 = dynamicsSimulator_1->getEquationsOfMotionNumericalSolution();
+
+ std::map< double, Eigen::VectorXd > integrationResult_2 = dynamicsSimulator_2->getEquationsOfMotionNumericalSolution();
+
+
+Results
+~~~~~~~
+
+The output for this configuration is:
+
+ .. code-block:: cpp
+
+ Initial spacecraft seaparation = 106.604km
+
+ Required impulse = 6.39207m/s
+
+ Estimated burn time = 21.5519s
+
+ Evolution n: 1 Objective Function: 2806.44
+ Evolution n: 2 Objective Function: 2323.05
+ [...]
+ Evolution n: 6 Objective Function: 2316.05
+ Evolution n: 7 Objective Function: 2315.48
+ [...]
+ Evolution n: 11 Objective Function: 2315.48
+ Evolution n: 12 Objective Function: 2315.3
+ [...]
+ Evolution n: 30 Objective Function: 2315.3
+
+ Calculated burn time = 21.3993s
+
+The calculated burn time is the optimized value. As you can see it is very much similar to the estimated one. Moreover, due to the geometry of the
+problem the optimizer fails to reach a value of 1000 m of separation, stopping at about 2315.3 m.
+
+An image of the simulation shows the problem with our assumptions:
+
+ .. figure:: images/RendezVousXY.svg
+ :scale: 100 %
+ :alt: x-y visualization of rendez-vous in target frame.
+
+ Path of chaser in the target-frame after retrograde burn, planar visualization
+
+
+ .. figure:: images/RendezVousXZ.svg
+ :scale: 100 %
+ :alt: x-z visualization of rendez-vous in target frame.
+
+ Path of chaser in the target-frame after retrograde burn, x-z visualization.
+
+The picures above show the path of the Chaser in the Target frame, where the x-direction is the radial direction from the centre of Earth,
+the z-direction is the direction of the angular momentum of the target w.r.t. the centre of Earth and the y-direction completes the
+right-hand rule of thumb frame.
+
+One can see that even if the two orbits are similar, the two spacecraft are separated initially in the radial direction
+by some 60 km. The burn, opposite to the direction of motion of the Chaser, has also a radial component. Nevertheless, the calculated burn
+brings the orbit at just 2315.3 meters from the target at its minimum separation.
+
+See also: `"Clohessy-Wiltshire equations" (PDF). University of Texas at Austin. Retrieved 23 October 2017. `_
+
+
+
diff --git a/tutorials/tudatFeatures/index.rst b/tutorials/tudatFeatures/index.rst
index 6c16b9589f..c2876100e6 100644
--- a/tutorials/tudatFeatures/index.rst
+++ b/tutorials/tudatFeatures/index.rst
@@ -119,3 +119,4 @@ These pages of the wiki provide further details about critical libraries necessa
propagationSetup/index
estimationSetup/index
otherLibraries/index
+ optimizationTools/index
diff --git a/tutorials/tudatFeatures/optimizationTools/decisionVariableSettings.rst b/tutorials/tudatFeatures/optimizationTools/decisionVariableSettings.rst
new file mode 100644
index 0000000000..b3fe3856b5
--- /dev/null
+++ b/tutorials/tudatFeatures/optimizationTools/decisionVariableSettings.rst
@@ -0,0 +1,152 @@
+.. _decisionVariableSettings:
+
+Decision variable settings
+==========================
+
+The decision variable settings classes allow the user to define some of the variables of a
+``SingleArcDynamicsSimulator`` and its members as decision variables for an optimization process.
+
+SingleDecisionVariableSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The class ``SingleDecisionVariableSettings`` is the base class which allows to define a single variable as
+decision variable.
+
+**Constructor:**
+
+ .. code-block:: cpp
+
+ SingleDecisionVariableSettings( DecisionVariables decisionVariable,
+ double lowerBoundary, double upperBoundary,
+ std::string associatedBody = "" );
+
+ SingleDecisionVariableSettings( DecisionVariables decisionVariable,
+ Eigen::VectorXd& lowerBoundary, Eigen::VectorXd& upperBoundary,
+ std::string associatedBody = "" );
+
+ * The ``decisionVariable`` can be either of the following enumerations:
+
+ * ``simulation_time_decision_variable``: the simulation time (from start epoch to end epoch) is tuned to obtain the optimal solution. It is a scalar quantity, therefore the first constructor can be used;
+ * ``initial_cartesian_state_decision_variable``: the initial Cartesian state is tuned to obtain the optimal value. This is a 6x1 matrix, therefore only the second constructor, with ``Eigen::Vector6d`` boundaries can be used;
+ * ``initial_cartesian_velocity_decision_variable``: the initial Cartesian velocity is tuned to obtain the optimal value. This is a 3x1 matrix, therefore only the second constructor, with ``Eigen::Vector3d`` boundaries can be used;
+ * ``initial_cartesian_position_decision_variable``: the initial Cartesian position is tuned to obtain the optimal value. This is a 3x1 matrix, therefore only the second constructor, with ``Eigen::Vector3d`` boundaries can be used.
+
+ * The ``lowerBoundary`` and ``upperBoundary`` represent the interval in which the optimum is localized. They can be set either as scalar, with the first constructor, or as vectors,
+ with the second constructor, depending on the dimension of the decision variable.
+
+ * The ``associatedBody`` is the name of the body, as determined in the ``bodyMap`` of the dynamics simulator, to which the decision variable applies.
+
+
+SingleDecisionVariableFromTerminationSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to use one of the propagation termination settings as decision variables. It can also be used to modulate the time of simulation,
+if a ``PropagationTimeTerminationSettings`` object is present in the propagator settings.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ SingleDecisionVariableFromTerminationSettings(
+ boost::shared_ptr< propagators::SingleArcDynamicsSimulator< double > > dynamicsSimulator,
+ const double lowerBoundary,
+ const double upperBoundary,
+ int positionInVectorOfTerminationSettings = 0 )
+
+
+ * The ``dynamicsSimulator`` is the dynamics simulator object in which the propagator settings object with the termination
+ settings is found
+ * The ``positionInVectorOfTerminationSettings`` is the index of the termination settings if a ``PropagatorHybridTerminationSettings`` is used.
+
+
+SingleCartesianComponentDecisionVariableSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to use one of the initial Cartesian components as decision variable.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ SingleCartesianComponentDecisionVariableSettings(
+ orbital_elements::CartesianElements cartesianComponent,
+ double lowerBoundary, double upperBoundary,
+ std::string associatedBody )
+
+
+ * The ``cartesianComponent`` can be chosen as either of the following enumerations:
+
+ * ``xCartesianPosition``
+ * ``yCartesianPosition``
+ * ``zCartesianPosition``
+ * ``xCartesianVelocity``
+ * ``yCartesianVelocity``
+ * ``zCartesianVelocity``
+
+
+SingleKeplerElementDecisionVariableSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to define one of the initial Kepler elements as the decision variable.
+
+
+**Constructor:**
+ .. code-block:: cpp
+
+ SingleKeplerElementDecisionVariableSettings(
+ orbital_elements::KeplerianElements keplerElement,
+ std::string centralBody, double lowerBoundary,
+ double upperBoundary, std::string associatedBody ) :
+
+ * The ``keplerElement`` can be chosen as either of the following enumerations:
+
+ * ``semiMajorAxis``
+ * ``eccentricity``
+ * ``inclination``
+ * ``argumentOfPeriapsis``
+ * ``longitudeOfAscendingNode``
+ * ``trueAnomaly``
+
+ * The ``centralBody`` is the name of the massive body containing the orbital parameter used to calculated the
+ Kepler elements. Beware that at the moment there is no option to use a different reference frame than the one used in
+ the simulation, therefore the ``centralBody`` should be the name of the associated central body of the
+ ``associatedBody`` in the propagator settings.
+
+
+
+SingleSphericalOrbitalElementDecisionVariableSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to define one of the initial spherical orbital elements as the decision variable.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ SingleSphericalOrbitalElementDecisionVariableSettings(
+ orbital_elements::SphericalOrbitalStateElements sphericalOrbitalElement,
+ double lowerBoundary, double upperBoundary, std::string associatedBody )
+
+ * The ``sphericalOrbitalElement`` can be chosen as either of the following enumerations:
+
+ * ``radius``
+ * ``latitude``
+ * ``longitude``
+ * ``speed``
+ * ``flightPath``
+ * ``headingAngle``
+
+
+DecisionVariableSettings
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+You can use this class if you want to set multiple decision variables for one mission segment.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ DecisionVariableSettings(
+ std::vector< boost::shared_ptr< SingleDecisionVariableSettings > > multiDecisionVariableSettings )
+
+ * ``multiDecisionVariableSettings`` is a vector of ``boost::shared_ptr< SingleDecisionVariableSettings >`` objects or any
+ of the above mentioned classes.
+
+
+
diff --git a/tutorials/tudatFeatures/optimizationTools/index.rst b/tutorials/tudatFeatures/optimizationTools/index.rst
new file mode 100644
index 0000000000..986d2df5b3
--- /dev/null
+++ b/tutorials/tudatFeatures/optimizationTools/index.rst
@@ -0,0 +1,13 @@
+.. _tudatOptimizationIndex:
+
+Optimization tools
+==================
+
+These pages of the wiki will explain the optimization namespace.
+
+.. toctree::
+
+ decisionVariableSettings
+ objectiveFunctionSettings
+ missionSegmentSettings
+ missionLinker
diff --git a/tutorials/tudatFeatures/optimizationTools/missionLinker.rst b/tutorials/tudatFeatures/optimizationTools/missionLinker.rst
new file mode 100644
index 0000000000..fd9392c5dc
--- /dev/null
+++ b/tutorials/tudatFeatures/optimizationTools/missionLinker.rst
@@ -0,0 +1,63 @@
+ .. _missionLinker:
+
+Mission linker
+==============
+
+The **mission linker** is the tool that allows to link together and optimize a chain of **mission segments**.
+Here it is represented by the class ``MissionLinker``
+
+MissionLinker
+~~~~~~~~~~~~~
+
+**Constructor**
+ .. code-block:: cpp
+
+ MissionLinker(
+ std::vector< boost::shared_ptr< MissionSegmentSettings > > missionSegmentSettings,
+ boost::shared_ptr< OptimizationSettings > optimizationSettings = boost::make_shared< OptimizationSettings >(),
+ const bool optimizeConfiguration = true )
+
+ * The ``missionSegmentSettings`` is the vector of ``MissionSegmentSettings`` objects defining the chain of mission segments to optimize.
+ * The ``optimizationSettings`` is a user defined object which allows some degrees of freedom over the optimization process.
+ * The ``optimizeConfiguration`` flag can be set to ``true`` to start the optimization right after creating the
+ MissionLinker object.
+
+
+**Methods**
+ .. code-block:: cpp
+
+ void optimize( void );
+
+This is the method to call if you want to optimize the configuration. All the optimum values as well as the results
+of the optimal propagation are subsequently stored in the dynamics simulators of each mission segment.
+
+OptimizationSettings
+~~~~~~~~~~~~~~~~~~~~
+
+This class allows the user some degrees of freedom over the optimization process.
+
+**Constructor**
+ .. code-block:: cpp
+
+ OptimizationSettings( OptimizationTypes optimizationType = global_optimization,
+ const int populationSize = 32,
+ const bool stopAtMaxNumberOfEvolutions = true,
+ const int verbosity = 0 )
+
+ * The ``optimizationType`` can be set to either of the following enumerations:
+
+ * ``global_optimization`` will use an evolutionary algorithm to find the optimum value (generally slower, but reliable),
+ * ``local_optimization`` will use a local derivative algorithm to find the optimum value (faster, less accurate);
+
+ * The ``populationSize`` determines the number of individual in a population. The larger the number, the larger the possibilities;
+ to converge in less iterations to the target value, although the optimization takes a proportional amount of time;
+ * The ``stopAtMaxNumberOfEvolutions`` flag can be set to ``false`` if you want your optimization to only stop if convergence is reached (not advisable);
+ * The ``verbosity`` can be set to:
+
+ * 0: no message will be shown during the optimization,
+ * n > 0: it will show the value of the objective function in a message every n iterations.
+
+Example
+~~~~~~~
+
+See an example that uses the ``MissionLinker`` at the page :ref:`satelliteRendezVous`.
diff --git a/tutorials/tudatFeatures/optimizationTools/missionSegmentSettings.rst b/tutorials/tudatFeatures/optimizationTools/missionSegmentSettings.rst
new file mode 100644
index 0000000000..4d81cef532
--- /dev/null
+++ b/tutorials/tudatFeatures/optimizationTools/missionSegmentSettings.rst
@@ -0,0 +1,67 @@
+ .. _missionSegmentSettings:
+
+Mission segment settings
+========================
+
+A **mission segment** in this framework is a piece of simulation placed in a chain with other pieces of simulation.
+Each simulation from the secon on in the chain starts at the end point of the previous simulation.
+
+The simulation is defined by a ``SingleArcDynamicsSimulator< double >`` object, and the mission segment can also have a **decision variable**
+and an **objective function**, which allow to optimize the chain using a :ref:`missionLinker`.
+
+The mission segment is here defined by the object ``MissionSegmentSettings``.
+
+MissionSegmentSettings
+~~~~~~~~~~~~~~~~~~~~~~
+
+The ``MissionSegmentSettings`` class allows to link a decision variable or an objective function to a dynamics simulator.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ MissionSegmentSettings(
+ boost::shared_ptr< propagators::SingleArcDynamicsSimulator< > > dynamicsSimulator,
+ boost::shared_ptr< DecisionVariableSettings > decisionVariableSettings = NULL,
+ boost::shared_ptr< ObjectiveFunctionSettings > objectiveFunctionSettings = NULL )
+
+
+ MissionSegmentSettings(
+ boost::shared_ptr< propagators::SingleArcDynamicsSimulator< > > dynamicsSimulator,
+ boost::shared_ptr< SingleDecisionVariableSettings > singleDecisionVariableSettings,
+ boost::shared_ptr< ObjectiveFunctionSettings > objectiveFunctionSettings = NULL )
+
+
+ MissionSegmentSettings(
+ boost::shared_ptr< propagators::SingleArcDynamicsSimulator< > > dynamicsSimulator,
+ boost::shared_ptr< ObjectiveFunctionSettings > objectiveFunctionSettings )
+
+ * The ``dynamicsSimulator`` is a ``SingleArcDynamicsSimulator< double >`` object containing the pre-assembled features
+ of the simulation;
+ * The ``decisionVariableSettings`` is a ``DecisionVariableSettings`` object containing multiple decision variables retrievable from
+ the ``dynamicsSimulator``;
+ * The ``singleDecisionVariableSettings`` can be used alternatively is only one decision variable is needed;
+ * The ``objectiveFunctionSettings`` is the objective function. A mission segment containing this object must
+ be placed at the end of the chain.
+
+Following the provided constructors one can create mission segments with the following combinations:
+
+ * Mission segment with **only a dynamics simulator**: if placed at the beginning or at the end of the chain this
+ mission segment is not involved in the optimization. It will anyway be linked to the other mission segments and
+ propagated in the optimizaton process if placed between a mission segment with a decision variable and a mission segment
+ with an objective function;
+ * Mission segment with a ``singleDecisionVariableSettings``: this mission segment cannot be placed at the end of the chain.
+ It contains only one decision variable;
+ * Mission segment with a ``singleDecisionVariableSettings`` and a ``objectiveFunctionSettings``: this mission segment can be
+ placed anywhere in the chain and it is sufficient for the optimization process. It contains only one decision variable;
+ * Mission segment with a ``decisionVariableSettings``: this mission segment cannot be placed at the end of the chain.
+ It can have multiple decision variables;
+ * Mission segment with a ``decisionVariableSettings`` and a ``objectiveFunctionSettings``: this mission segment can be
+ placed anywhere in the chain and it is sufficient for the optimization process. It can have multiple decision variables;
+ * Mission segment with an ``objectiveFunctionSettings``: this mission segment must be preceded by at least one
+ mission segment with a ``singleDecisionVariableSettings`` or a ``decisionVariableSettings`` and without
+ an ``objectiveFunctionSettings``.
+
+
+
+
+
diff --git a/tutorials/tudatFeatures/optimizationTools/objectiveFunctionSettings.rst b/tutorials/tudatFeatures/optimizationTools/objectiveFunctionSettings.rst
new file mode 100644
index 0000000000..9192372520
--- /dev/null
+++ b/tutorials/tudatFeatures/optimizationTools/objectiveFunctionSettings.rst
@@ -0,0 +1,327 @@
+.. _objectiveFunctionSettings:
+
+Objective function settings
+===========================
+
+The objective function settings allows to define what is the objective function of the optimization process. Several classes can be
+used, each focused on a certain type of objective function.
+
+ .. warning:: At the moment only single objective function optimization is allowed.
+
+FinalCartesianComponentObjectiveFunctionSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to use one of the final cartesian components of a propagated body as objective function.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ FinalCartesianComponentObjectiveFunctionSettings(
+ orbital_elements::CartesianElements cartesianComponent,
+ const std::string associatedBody,
+ const double objectiveValue,
+ const double tolerance = 1e-3,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+
+ FinalCartesianComponentObjectiveFunctionSettings(
+ orbital_elements::CartesianElements cartesianComponent,
+ const std::string associatedBody,
+ const bool minimizeOrMaximize = true,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+
+The first constructor allows to define an objective value towards which the objective function converges. The second constructor will
+continue to optimize the configuration until the maximum number of evolutions is reached.
+
+ * The ``cartesianComponent`` can be chosed as either of the following enumerations:
+
+ * ``xCartesianPosition``
+ * ``yCartesianPosition``
+ * ``zCartesianPosition``
+ * ``xCartesianVelocity``
+ * ``yCartesianVelocity``
+ * ``zCartesianVelocity``
+
+ * The ``associatedBody`` is the name of the propagated body, as determined in the ``bodyMap`` of the dynamics simulator to which the
+ final cartesian component belongs.
+ * The ``objectiveValue`` is the value towards which the optimization converges.
+ * The ``tolerance`` is the error between the final value and the objective function for which the optimization is allowed to stop.
+ * The ``maxNumberOfEvolutions`` is the maximum number of iterations allowed in the optimization process if the convergence is not reached.
+ * The ``minimizeOrMaximize`` flag is set to:
+
+ * ``true`` if you want to minimize the objective function,
+ * ``false`` if you want to maximize the objective function.
+
+FinalKeplerElementObjectiveFunctionSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to use one of the final Kepler orbital elements of a propagated body as objective function.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ FinalKeplerElementObjectiveFunctionSettings(
+ orbital_elements::KeplerianElements keplerElement,
+ const std::string associatedBody,
+ std::string centralBody,
+ const double objectiveValue,
+ const double tolerance = 1e-3,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+ FinalKeplerElementObjectiveFunctionSettings(
+ orbital_elements::KeplerianElements keplerElement,
+ const std::string associatedBody,
+ std::string centralBody,
+ const bool minimizeOrMaximize = true,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+The first constructor allows to define an objective value towards which the objective function converges. The second constructor will
+continue to optimize the configuration until the maximum number of evolutions is reached.
+
+ * The ``keplerElement`` can be chosen as either of the following enumerations:
+
+ * ``semiMajorAxis``
+ * ``eccentricity``
+ * ``inclination``
+ * ``argumentOfPeriapsis``
+ * ``longitudeOfAscendingNode``
+ * ``trueAnomaly``
+
+ * The ``associatedBody`` is the name of the propagated body, as determined in the ``bodyMap`` of the dynamics simulator to which the
+ final cartesian component belongs.
+ * The ``centralBody`` is the name of the massive body containing the orbital parameter used to calculated the
+ Kepler elements. Beware that at the moment there is no option to use a different reference frame than the one used in
+ the simulation, therefore the ``centralBody`` should be the name of the associated central body of the
+ ``associatedBody`` in the propagator settings.
+ * The ``objectiveValue`` is the value towards which the optimization converges.
+ * The ``tolerance`` is the error between the final value and the objective function for which the optimization is allowed to stop.
+ * The ``maxNumberOfEvolutions`` is the maximum number of iterations allowed in the optimization process if the convergence is not reached.
+ * The ``minimizeOrMaximize`` flag is set to:
+
+ * ``true`` if you want to minimize the objective function,
+ * ``false`` if you want to maximize the objective function.
+
+FinalSphericalOrbitalElementObjectiveFunctionSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to use one of the final spherical orbital elements of a propagated body as objective function.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ FinalSphericalOrbitalElementObjectiveFunctionSettings(
+ orbital_elements::SphericalOrbitalStateElements spericalOrbitalElement,
+ const std::string associatedBody,
+ const double objectiveValue,
+ const double tolerance = 1e-3,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+ FinalSphericalOrbitalElementObjectiveFunctionSettings(
+ orbital_elements::SphericalOrbitalStateElements spericalOrbitalElement,
+ const std::string associatedBody,
+ const bool minimizeOrMaximize = true,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+The first constructor allows to define an objective value towards which the objective function converges. The second constructor will
+continue to optimize the configuration until the maximum number of evolutions is reached.
+
+ * The ``sphericalOrbitalElement`` can be chosen as either of the following enumerations:
+
+ * ``radius``
+ * ``latitude``
+ * ``longitude``
+ * ``speed``
+ * ``flightPath``
+ * ``headingAngle``
+
+ * The ``associatedBody`` is the name of the propagated body, as determined in the ``bodyMap`` of the dynamics simulator to which the
+ final cartesian component belongs.
+ * The ``objectiveValue`` is the value towards which the optimization converges.
+ * The ``tolerance`` is the error between the final value and the objective function for which the optimization is allowed to stop.
+ * The ``maxNumberOfEvolutions`` is the maximum number of iterations allowed in the optimization process if the convergence is not reached.
+ * The ``minimizeOrMaximize`` flag is set to:
+
+ * ``true`` if you want to minimize the objective function,
+ * ``false`` if you want to maximize the objective function.
+
+ObjectiveFunctionFromFinalDependentVariableSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to use the final value of one of the dependent variables saved by the propagator as objective function.
+
+ .. warning:: You need to have set the DependentVariableSaveSettings inside the PropagatorSettings in order to use this class!
+
+**Constructor:**
+ .. code-block:: cpp
+
+ ObjectiveFunctionFromFinalDependentVariableSettings(
+ boost::shared_ptr< propagators::SingleDependentVariableSaveSettings > dependentVariableSaveSettings,
+ const double objectiveValue,
+ const int indexInVectorialVariable = 0,
+ const double tolerance = 1e-3,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+ ObjectiveFunctionFromFinalDependentVariableSettings(
+ boost::shared_ptr< propagators::SingleDependentVariableSaveSettings > dependentVariableSaveSettings,
+ const bool minimizeOrMaximize,
+ const int indexInVectorialVariable = 0,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+The first constructor allows to define an objective value towards which the objective function converges. The second constructor will
+continue to optimize the configuration until the maximum number of evolutions is reached.
+
+ * The ``dependentVariableSettings`` is the ``SingleDependentVariableSaveSettings`` object belonging to the saved variable whose final value
+ you want to use as objective function.
+ * The ``objectiveValue`` is the value towards which the optimization converges.
+ * The ``indexInVectorialVariable`` is the index pointing to the position that you need in the dependent variable in case the latter
+ is not scalar.
+ * The ``tolerance`` is the error between the final value and the objective function for which the optimization is allowed to stop.
+ * The ``maxNumberOfEvolutions`` is the maximum number of iterations allowed in the optimization process if the convergence is not reached.
+ * The ``minimizeOrMaximize`` flag is set to:
+
+ * ``true`` if you want to minimize the objective function,
+ * ``false`` if you want to maximize the objective function.
+
+ObjectiveFunctionFromMinOrMaxDependentVariableSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to use the minimum or maximum value in the propagation history of one of the dependent variables saved by the propagator as objective function.
+
+ .. warning:: You need to have set the DependentVariableSaveSettings inside the PropagatorSettings in order to use this class!
+
+**Constructor:**
+ .. code-block:: cpp
+
+ ObjectiveFunctionFromMinOrMaxDependentVariableSettings(
+ boost::shared_ptr< propagators::SingleDependentVariableSaveSettings > dependentVariableSaveSettings,
+ const double objectiveValue,
+ const int indexInVectorialVariable = 0,
+ const bool minimumOrMaximum = true,
+ const double tolerance = 1e-3,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+ ObjectiveFunctionFromMinOrMaxDependentVariableSettings(
+ boost::shared_ptr< propagators::SingleDependentVariableSaveSettings > dependentVariableSaveSettings,
+ const int indexInVectorialVariable = 0,
+ const bool minimizeOrMaximize = true,
+ const bool minimumOrMaximum = true,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+The first constructor allows to define an objective value towards which the objective function converges. The second constructor will
+continue to optimize the configuration until the maximum number of evolutions is reached.
+
+ * The ``dependentVariableSettings`` is the ``SingleDependentVariableSaveSettings`` object belonging to the saved variable whose final value
+ you want to use as objective function.
+ * The ``objectiveValue`` is the value towards which the optimization converges.
+ * The ``indexInVectorialVariable`` is the index pointing to the position that you need in the dependent variable in case the latter
+ is not scalar.
+ * The ``minimumOrMaximum`` flag is set to:
+
+ * ``true`` if you want to use the **minimum** in the propagation history of the variable,
+ * ``false`` if you want to use the **maximum** in the propagation history of the variable.
+
+ * The ``tolerance`` is the error between the final value and the objective function for which the optimization is allowed to stop.
+ * The ``maxNumberOfEvolutions`` is the maximum number of iterations allowed in the optimization process if the convergence is not reached.
+ * The ``minimizeOrMaximize`` flag is set to:
+
+ * ``true`` if you want to minimize the objective function,
+ * ``false`` if you want to maximize the objective function.
+
+
+UserDefinedObjectiveFunctionSettings
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This class allows to setup an objective function which is not available among the above. It could also be used to define
+a weighted multi-objective optimization problem.
+
+See the paragraph on :ref:`howToSetupUserDefinedObjectiveFunction` for information on how to use this class.
+
+**Constructor:**
+ .. code-block:: cpp
+
+ UserDefinedObjectiveFunctionSettings(
+ boost::function< double () > userDefinedFunction,
+ const double objectiveValue = 0.0,
+ const double tolerance = 1e-3,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+ UserDefinedObjectiveFunctionSettings(
+ boost::function< double () > userDefinedFunction,
+ const bool minimizeOrMaximize,
+ const unsigned int maxNumberOfEvolutions = 100 )
+
+The first constructor allows to define an objective value towards which the objective function converges. The second constructor will
+continue to optimize the configuration until the maximum number of evolutions is reached.
+
+ * The ``maxNumberOfEvolutions`` is the maximum number of iterations allowed in the optimization process if the convergence is not reached.
+ * The ``minimizeOrMaximize`` flag is set to:
+
+ * ``true`` if you want to minimize the objective function,
+ * ``false`` if you want to maximize the objective function.
+
+.. _howToSetupUserDefinedObjectiveFunction:
+
+How to set up a user defined objective function
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+You might have noticed the use of ``boost::function< double () >`` and the fact that the ``UserDefinedObjectiveFunction`` class does not allow
+an objective function with parameters. How do you setup an objective function without parameters?
+
+First of all you will need a support class that takes all the parameters that you need, possibly as pointers, and a ``double`` method that does
+not accept parameters. You can store in this object pointers to e.g. your dynamics simulator or propagator settings or whatever.
+
+ .. code-block:: cpp
+
+ struct ClassToDefineFunction{
+
+ //Constructor. Make sure to input all the needed values
+ //as pointers or boost::shared_ptr
+ ClassToDefineFunction( argsType args ) :
+ members_( args )
+ {
+ ...
+ }
+
+ ~ClassToDefineFunction( ){ }
+
+ //This is the function that provides the calculated
+ //value of the objective function.
+ double myFunction( void )
+ {
+ double value;
+
+ ...
+
+ return value;
+ }
+
+ memberType members_
+
+ };
+
+
+Afterwards, when creating your objective function from ``UserDefinedObjectiveFunction`` you will need to use the ``boost::bind()`` statement to bind the method to a
+``boost::function< double () >`` object. All you need to know is the following sintax, which is the only sintax that you are allowed to use in such a case:
+
+ .. code-block:: cpp
+
+ //Create object from class to define function and input all
+ //the needed arguments
+ boost::shared_ptr< ClassToDefineFunction > objectToDefineFunction =
+ boost::make_shared< ClassToDefineFunction >( args );
+
+ //Create boost::function object from the above object
+ boost::function< double () > myObjectiveFunction =
+ boost::bind( &ClassToDefineFunction::myFunction,
+ objectToDefineFunction )
+
+ //Create object to retrieve objective function from simulation
+ boost::shared_ptr< UserDefinedObjectiveFunction > myObjectiveFunctionObject =
+ boost::make_shared< UserDefinedObjectiveFunction >( myObjectiveFunction )
+
+
+Et voilà. The user defined objective function is set.
+
+
+
diff --git a/tutorials/tudatFeatures/otherLibraries/index.rst b/tutorials/tudatFeatures/otherLibraries/index.rst
index 60aa8d5521..078d65fa9b 100644
--- a/tutorials/tudatFeatures/otherLibraries/index.rst
+++ b/tutorials/tudatFeatures/otherLibraries/index.rst
@@ -10,3 +10,4 @@ These pages of the wiki will help you build a strong knowledge basis to get star
inputOutput
spice
jsoncpp
+ pagmo2
diff --git a/tutorials/tudatFeatures/otherLibraries/pagmo2.rst b/tutorials/tudatFeatures/otherLibraries/pagmo2.rst
new file mode 100644
index 0000000000..0ad66c51de
--- /dev/null
+++ b/tutorials/tudatFeatures/otherLibraries/pagmo2.rst
@@ -0,0 +1,231 @@
+.. _tudatFeaturesPagmo2:
+
+External Libraries: PaGMO 2
+===========================
+
+
+Download and build PaGMO 2
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+Clone or download PaGMO 2 from::
+
+ https://github.com/esa/pagmo2
+
+Make sure to name the source folder :literal:`pagmo2/`.
+
+.. warning:: You need Boost version 1.55 or higher in order to build PaGMO 2.
+
+The PaGMO project contains its own CMakeLists.txt. In order to build PaGMO as a standalone release version simply run CMake in the folder.
+In order to include PaGMO in your project you will need to specify whether you want to build PaGMO or PyGMO (Python interface for PaGMO).
+You can either set either of the two options :literal:`PAGMO_BUILD_PAGMO` or :literal:`PAGMO_BUILD_PYGMO` according to your needs in the CMakeLists.txt file.
+
+.. warning:: Only one of the two flags can be set to ON at anytime, so make sure to switch the other to OFF (by default PAGMO_BUILD_PAGMO is set to ON).
+
+If you are building PaGMO inside your project, you can override the two options using the command :literal:`-option()`
+
+ .. code-block:: cmake
+
+ option(PAGMO_BUILD_PAGMO ON)
+ option(PAGMO_BUILD_PYGMO OFF)
+
+You can simply include PaGMO 2 by adding the source folder as a subdirectory.
+
+ .. code-block:: cmake
+
+ unset(pagmo_LIB_DEPENDS CACHE)
+ add_subdirectory( "my_path/pagmo2/" )
+
+The command :literal:`unset(pagmo_LIB_DEPENDS CACHE)` makes sure to eliminate residual dependencies in the cache in case of re-build of the package.
+
+Create an optimization problem with PaGMO
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Let us define the Himmelblau's function as a PaGMO problem.
+The Himmelblau's function is defined as: :math:`f(x)=(x^2+y-11)^2 + (x+y^2-7)^2` and has
+
+ * one local maximum :math:`f(− 0.270845, − 0.923039) = 181.617`;
+ * four local minima:
+
+ * :math:`f( 3.0 , 2.0 ) = 0.0`;
+ * :math:`f ( − 2.805118 , 3.131312 ) = 0.0`;
+ * :math:`f ( − 3.779310 , − 3.283186 ) = 0.0`;
+ * :math:`f ( 3.584428 , − 1.848126 ) = 0.0`.
+
+A PaGMO problem in C++ is defined as a ``struct`` with the following mandatory methods:
+
+ * ``std::vector< double > fitness( const std::vector< double > &x ) const``
+ * ``std::pair< std::vector< double >, std::vector< double > > get_bounds( ) const``
+
+In a file ``himmelblaus.h`` let us combine the himmelblau's function and the problem and allow the user to define the **box-boundaries** for the optimization search.
+
+ .. code-block:: cpp
+
+ //File: himmelblaus.h
+
+ struct my_problem {
+
+ //Empty constructor
+ my_problem( ){ }
+
+ //Actual constructor allowing the user to define the boundaries
+ my_problem( const double x_min, const double x_max, const double y_min,
+ const double y_max ) :
+ x_min_( x_min ), x_max_( x_max ), y_min_( y_min ), y_max( y_max_ )
+ { }
+
+ // Mandatory, computes the fitness, i.e. the Himmelblau's function
+ std::vector< double > fitness( const std::vector< double > &x ) const
+ {
+ std::vetor< double > return_value;
+
+ return_value.push_back( pow( x[0]*x[0] + x[1] - 11, 2) + pow( x[0] + x[1]*x[1] - 7, 2 ) );
+
+ return return_value;
+
+ }
+
+ // Mandatory, returns the box-bounds
+ std::pair< std::vector< double >, std::vector< double > > get_bounds( ) const
+ {
+ std::pair< std::vector< double >, std::vector< double > > box_bounds;
+
+ box_bounds.first.push_back( x_min_ );
+ box_bounds.first.push_back( y_min_ );
+
+ box_bounds.second.push_back( x_max_ );
+ box_bounds.second.push_back( y_max_ );
+
+ return box_bounds;
+ }
+
+ private:
+
+ //Storage members
+ double x_min_;
+ double x_max_;
+ double y_min_;
+ double y_max_;
+
+ };
+
+The **empty constructor** used in the example above will avoid multi-threading identity related problems associated with the ``pthread`` library.
+
+The use of PaGMO enbedded alghorithms and the user-defined problem is demonstrated in the following code block. Let us define a file called ``main.cpp``.
+
+ .. code-block:: cpp
+
+ //File: main.cpp
+
+ #include
+ #include"pagmo/algorithms/sade.hpp"
+ #include"pagmo/island.hpp"
+ #include"pagmo/problem.hpp"
+ #include"himmelblaus.h"
+
+ int main( )
+ {
+
+ // Create a PaGMO problem from the user-defined Himmelblau's problem
+ // Search within 0:${PAGMO_CXX_FLAGS_DEBUG}>" "$<$:${PAGMO_CXX_FLAGS_RELEASE}>")
+
+ # Set C++ standard to C++11
+ set_property(TARGET application_PagmoHimmelblausExample PROPERTY CXX_STANDARD 11)
+ set_property(TARGET application_PagmoHimmelblausExample PROPERTY CXX_STANDARD_REQUIRED YES)
+ set_property(TARGET application_PagmoHimmelblausExample PROPERTY CXX_EXTENSIONS NO)
+
+ # Link Boost libraries and pthread
+ target_link_libraries(application_PagmoHimmelblausExample ${Boost_LIBRARIES} pthread)
+
+
+
+
+
+
+
+
+
+
+
+