integrator_bdf.hpp

Go to the documentation of this file.

/*
 *    This file is part of ACADO Toolkit.
 *
 *    ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
 *    Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
 *    Milan Vukov, Rien Quirynen, KU Leuven.
 *    Developed within the Optimization in Engineering Center (OPTEC)
 *    under supervision of Moritz Diehl. All rights reserved.
 *
 *    ACADO Toolkit is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation; either
 *    version 3 of the License, or (at your option) any later version.
 *
 *    ACADO Toolkit is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *    Lesser General Public License for more details.
 *
 *    You should have received a copy of the GNU Lesser General Public
 *    License along with ACADO Toolkit; if not, write to the Free Software
 *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 *
 */



#ifndef ACADO_TOOLKIT_INTEGRATOR_BDF_HPP
#define ACADO_TOOLKIT_INTEGRATOR_BDF_HPP

#include <acado/integrator/integrator_fwd.hpp>

BEGIN_NAMESPACE_ACADO


class IntegratorBDF : public Integrator{


friend class SimulationByIntegration;
friend class ShootingMethod;

//
// PUBLIC MEMBER FUNCTIONS:
//

public:

    IntegratorBDF( );

    IntegratorBDF( const DifferentialEquation &rhs_ );

    IntegratorBDF( const IntegratorBDF& arg );

    virtual ~IntegratorBDF( );

    virtual IntegratorBDF& operator=( const IntegratorBDF& arg );


    virtual Integrator* clone() const;


   // ================================================================================

    virtual returnValue init( const DifferentialEquation &rhs_ );


    inline returnValue init( const DifferentialEquation &rhs_,
                             const Transition           &trs_ );


    // ================================================================================

    virtual returnValue freezeMesh();


    virtual returnValue freezeAll();


    virtual returnValue unfreeze();


    // ================================================================================


    virtual returnValue step( int number   );


    virtual returnValue stop();


    virtual returnValue setDxInitialization( double *dx0   );



    // ================================================================================



    virtual int getNumberOfSteps() const;


    virtual int getNumberOfRejectedSteps() const;



    virtual double getStepSize() const;



//
// PROTECTED MEMBER FUNCTIONS:
//
protected:



    virtual returnValue evaluate( const DVector &x0    ,
                                  const DVector &xa    ,
                                  const DVector &p     ,
                                  const DVector &u     ,
                                  const DVector &w     ,
                                  const Grid   &t_      );

    // ================================================================================


    virtual returnValue evaluateSensitivities();


    // ================================================================================


    virtual returnValue setProtectedForwardSeed( const DVector &xSeed     ,
                                                 const DVector &pSeed     ,
                                                 const DVector &uSeed     ,
                                                 const DVector &wSeed     ,
                                                 const int    &order      );


    // ================================================================================


    virtual returnValue setProtectedBackwardSeed(  const DVector &seed   ,
                                                   const int    &order    );


    // ================================================================================



    virtual returnValue getProtectedX(           DVector *xEnd  ) const;


    virtual returnValue getProtectedForwardSensitivities( DMatrix *Dx  ,
                                                          int order    ) const;



    virtual returnValue getProtectedBackwardSensitivities( DVector &Dx_x0,
                                                           DVector &Dx_p ,
                                                           DVector &Dx_u ,
                                                           DVector &Dx_w ,
                                                           int order      ) const;


    // ================================================================================

    virtual int getDim() const;


    virtual int getDimX() const;


    void initializeButcherTableau();

    returnValue rk_start();


    returnValue rk_start_solve( int stepnumber );


    double determinePredictor( int number, BooleanType ini );


    returnValue determineCorrector( int number, BooleanType ini );


    void printBDFIntermediateResults();


    void printBDFfinalResults();


    void printBDFfinalResults2(DMatrix &div);


    void printRKIntermediateResults();


    returnValue decomposeJacobian(int index, DMatrix &J );


    double applyNewtonStep( int index, double *etakplus1, const double *etak, const DMatrix &J, const double *FFF );


    void applyMTranspose( int index, double *seed1, DMatrix &J, double *seed2 );


    returnValue setForwardSeed2( const DVector &xSeed           ,
                                 const DVector &pSeed           ,
                                 const DVector &uSeed           ,
                                 const DVector &wSeed              );


    virtual returnValue setBackwardSeed2( const DVector &seed      );


    void determineRKEtaGForward();


    void determineRKEtaHBackward();


    void determineRKEtaGForward2();


    void determineRKEtaHBackward2();


    void determineBDFEtaGForward( int number );

    void determineBDFEtaHBackward( int number );

    void determineBDFEtaGForward2( int number );

    void determineBDFEtaHBackward2( int number );


    void deleteAll();


    void constructAll( const IntegratorBDF& arg );



    void allocateMemory( );


    void initializeVariables();


    void relaxAlgebraic( double *residuum, double timePoint );


    void prepareDividedDifferences( DMatrix &div );

    void interpolate( int number_, DMatrix &div, VariablesGrid &poly );

        void logCurrentIntegratorStep(  const DVector& currentX  = emptyConstVector,
                                                                        const DVector& currentXA = emptyConstVector
                                                                        );


    void copyBackward( DVector       &Dx_x0,
                       DVector       &Dx_p ,
                       DVector       &Dx_u ,
                       DVector       &Dx_w ,
                       const DMatrix &div    ) const;


// DATA MEMBERS:
//
protected:


    // STORAGE OF THE ALGEBRAIC RESIDUUM :
    // -----------------------------------
    double *initialAlgebraicResiduum;
    double        relaxationConstant;


    // BUTCHER-
    // TABLEAU:
    // ----------
    int      dim               ;  
    double **A                 ;  
    double  *b4                ;  
    double  *b5                ;  
    double  *c                 ;  
    // RK-STARTER:
    // -----------
    double   *eta4             ;  
    double   *eta5             ;  
    double   *eta4_            ;  
    double   *eta5_            ;  
    double ***k                ;  
    double ***k2               ;  
    double ***l                ;  
    double ***l2               ;  
    double   *iseed            ;  
    // BDF-METHOD:
    // -----------
    int      nstep             ; 
    double **psi               ; 
    double  *psi_              ; 
    double **gamma             ; 
    DMatrix   nablaY            ; 
    DMatrix   nablaY_           ; 
    DMatrix   phi               ; 
    DMatrix   delta             ; 
    DVector   c2                ; 
    DMatrix **M                 ; 
//    std::vector< Eigen::HouseholderQR< DMatrix::Base > > qr;
    int     *M_index           ; 
    int      nOfM              ; 
    int      maxNM             ; 
    int     *nOfNewtonSteps    ; 
    double **eta               ; 
    double **eta2              ; 
    double   t                 ; 
    double  *x                 ; 
    double  *F                 ;
    double  *F2                ;


    // OTHERS:
    // -------
    double *initial_guess      ;  
    double  rel_time_scale     ;  
    // SENSITIVITIES:
    // --------------
    int        ndir            ;  
    DVector     fseed           ;  
    DVector     bseed           ;  
    DVector     fseed2          ;  
    DVector     bseed2          ;  
    double    *G               ;  
    double    *etaG            ;  
    DMatrix     nablaG          ;  
    DMatrix     phiG            ;  
    DMatrix     deltaG          ;  
    DVector     c2G             ;  
    double    *G2              ;  
    double    *G3              ;  
    double    *etaG2           ;  
    double    *etaG3           ;  
    DMatrix     nablaG2         ;  
    DMatrix     nablaG3         ;  
    DMatrix     phiG2           ;  
    DMatrix     phiG3           ;  
    DMatrix     deltaG2         ;  
    DMatrix     deltaG3         ;  
    DVector     c2G2            ;  
    DVector     c2G3            ;  
    double    *H               ;  
    double   **etaH            ;  
    DMatrix     nablaH          ;  
    DMatrix     nablaH_         ;  
    DMatrix     deltaH          ;  
    DVector     c2H             ;  
    double  ***kH              ;  
    double   **zH              ;  
    double    *H2              ;  
    double    *H3              ;  
    double   **etaH2           ;  
    double   **etaH3           ;  
    DMatrix     nablaH2         ;  
    DMatrix     nablaH3         ;  
    DMatrix     nablaH2_        ;  
    DMatrix     nablaH3_        ;  
    DMatrix     deltaH2         ;  
    DMatrix     deltaH3         ;  
    DVector     c2H2            ;  
    DVector     c2H3            ;  
    double  ***kH2             ;  
    double   **zH2             ;  
    double  ***kH3             ;  
    double   **zH3             ;  
    // STORAGE:
    // --------
    int maxAlloc               ;  
    // STATISTICS:
    // -----------
    RealClock jacComputation    ;
    RealClock jacDecomposition  ;
    RealClock correctorTime     ;
    int       nJacEvaluations   ;

};


CLOSE_NAMESPACE_ACADO



#include <acado/integrator/integrator_bdf.ipp>


#endif  // ACADO_TOOLKIT_INTEGRATOR_BDF_HPP

// end of file.