acado: integrator_bdf.hpp Source File

Go to the documentation of this file.
00001 /*
00002  *    This file is part of ACADO Toolkit.
00003  *
00004  *    ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
00005  *    Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
00006  *    Milan Vukov, Rien Quirynen, KU Leuven.
00007  *    Developed within the Optimization in Engineering Center (OPTEC)
00008  *    under supervision of Moritz Diehl. All rights reserved.
00009  *
00010  *    ACADO Toolkit is free software; you can redistribute it and/or
00011  *    modify it under the terms of the GNU Lesser General Public
00012  *    License as published by the Free Software Foundation; either
00013  *    version 3 of the License, or (at your option) any later version.
00014  *
00015  *    ACADO Toolkit is distributed in the hope that it will be useful,
00016  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018  *    Lesser General Public License for more details.
00019  *
00020  *    You should have received a copy of the GNU Lesser General Public
00021  *    License along with ACADO Toolkit; if not, write to the Free Software
00022  *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00023  *
00024  */
00025 
00026 
00027 
00034 #ifndef ACADO_TOOLKIT_INTEGRATOR_BDF_HPP
00035 #define ACADO_TOOLKIT_INTEGRATOR_BDF_HPP
00036 
00037 #include <acado/integrator/integrator_fwd.hpp>
00038 
00039 BEGIN_NAMESPACE_ACADO
00040 
00041 
00052 class IntegratorBDF : public Integrator{
00053 
00054 
00055 friend class SimulationByIntegration;
00056 friend class ShootingMethod;
00057 
00058 //
00059 // PUBLIC MEMBER FUNCTIONS:
00060 //
00061 
00062 public:
00063 
00065     IntegratorBDF( );
00066 
00068     IntegratorBDF( const DifferentialEquation &rhs_ );
00069 
00071     IntegratorBDF( const IntegratorBDF& arg );
00072 
00074     virtual ~IntegratorBDF( );
00075 
00077     virtual IntegratorBDF& operator=( const IntegratorBDF& arg );
00078 
00079 
00081     virtual Integrator* clone() const;
00082 
00083 
00084    // ================================================================================
00085 
00094     virtual returnValue init( const DifferentialEquation &rhs_ );
00095 
00096 
00108     inline returnValue init( const DifferentialEquation &rhs_,
00109                              const Transition           &trs_ );
00110 
00111 
00112     // ================================================================================
00113 
00123     virtual returnValue freezeMesh();
00124 
00125 
00134     virtual returnValue freezeAll();
00135 
00136 
00144     virtual returnValue unfreeze();
00145 
00146 
00147     // ================================================================================
00148 
00149 
00172     virtual returnValue step( int number   );
00173 
00174 
00183     virtual returnValue stop();
00184 
00185 
00190     virtual returnValue setDxInitialization( double *dx0   );
00194 
00195 
00196 
00197     // ================================================================================
00198 
00199 
00200 
00204     virtual int getNumberOfSteps() const;
00205 
00206 
00210     virtual int getNumberOfRejectedSteps() const;
00211 
00212 
00213 
00215     virtual double getStepSize() const;
00216 
00217 
00218 
00219 //
00220 // PROTECTED MEMBER FUNCTIONS:
00221 //
00222 protected:
00223 
00224 
00225 
00229     virtual returnValue evaluate( const DVector &x0    ,
00230                                   const DVector &xa    ,
00231                                   const DVector &p     ,
00232                                   const DVector &u     ,
00233                                   const DVector &w     ,
00234                                   const Grid   &t_      );
00235 
00236     // ================================================================================
00237 
00238 
00244     virtual returnValue evaluateSensitivities();
00245 
00246 
00247     // ================================================================================
00248 
00249 
00254     virtual returnValue setProtectedForwardSeed( const DVector &xSeed     ,
00256                                                  const DVector &pSeed     ,
00258                                                  const DVector &uSeed     ,
00260                                                  const DVector &wSeed     ,
00262                                                  const int    &order      );
00264 
00265 
00266     // ================================================================================
00267 
00268 
00273     virtual returnValue setProtectedBackwardSeed(  const DVector &seed   ,
00274                                                    const int    &order    );
00276 
00277 
00278     // ================================================================================
00279 
00280 
00281 
00285     virtual returnValue getProtectedX(           DVector *xEnd  ) const;
00288 
00289 
00294     virtual returnValue getProtectedForwardSensitivities( DMatrix *Dx  ,
00297                                                           int order    ) const;
00298 
00299 
00300 
00312     virtual returnValue getProtectedBackwardSensitivities( DVector &Dx_x0,
00313                                                            DVector &Dx_p ,
00314                                                            DVector &Dx_u ,
00315                                                            DVector &Dx_w ,
00316                                                            int order      ) const;
00317 
00318 
00319     // ================================================================================
00320 
00322     virtual int getDim() const;
00323 
00324 
00326     virtual int getDimX() const;
00327 
00328 
00330     void initializeButcherTableau();
00331 
00336     returnValue rk_start();
00337 
00338 
00342     returnValue rk_start_solve( int stepnumber );
00343 
00344 
00349     double determinePredictor( int number, BooleanType ini );
00350 
00351 
00356     returnValue determineCorrector( int number, BooleanType ini );
00357 
00358 
00361     void printBDFIntermediateResults();
00362 
00363 
00366     void printBDFfinalResults();
00367 
00368 
00371     void printBDFfinalResults2(DMatrix &div);
00372 
00373 
00376     void printRKIntermediateResults();
00377 
00378 
00383     returnValue decomposeJacobian(int index, DMatrix &J );
00384 
00385 
00389     double applyNewtonStep( int index, double *etakplus1, const double *etak, const DMatrix &J, const double *FFF );
00390 
00391 
00396     void applyMTranspose( int index, double *seed1, DMatrix &J, double *seed2 );
00397 
00398 
00401     returnValue setForwardSeed2( const DVector &xSeed           ,
00403                                  const DVector &pSeed           ,
00405                                  const DVector &uSeed           ,
00407                                  const DVector &wSeed              );
00409 
00410 
00413     virtual returnValue setBackwardSeed2( const DVector &seed      );
00415 
00416 
00419     void determineRKEtaGForward();
00420 
00421 
00424     void determineRKEtaHBackward();
00425 
00426 
00429     void determineRKEtaGForward2();
00430 
00431 
00434     void determineRKEtaHBackward2();
00435 
00436 
00439     void determineBDFEtaGForward( int number );
00440 
00443     void determineBDFEtaHBackward( int number );
00444 
00447     void determineBDFEtaGForward2( int number );
00448 
00451     void determineBDFEtaHBackward2( int number );
00452 
00453 
00456     void deleteAll();
00457 
00458 
00461     void constructAll( const IntegratorBDF& arg );
00462 
00463 
00464 
00471     void allocateMemory( );
00472 
00473 
00479     void initializeVariables();
00480 
00481 
00483     void relaxAlgebraic( double *residuum, double timePoint );
00484 
00485 
00486     void prepareDividedDifferences( DMatrix &div );
00487 
00488     void interpolate( int number_, DMatrix &div, VariablesGrid &poly );
00489 
00490         void logCurrentIntegratorStep(  const DVector& currentX  = emptyConstVector,
00491                                                                         const DVector& currentXA = emptyConstVector
00492                                                                         );
00493 
00494 
00495     void copyBackward( DVector       &Dx_x0,
00496                        DVector       &Dx_p ,
00497                        DVector       &Dx_u ,
00498                        DVector       &Dx_w ,
00499                        const DMatrix &div    ) const;
00500 
00501 
00502 // DATA MEMBERS:
00503 //
00504 protected:
00505 
00506 
00507     // STORAGE OF THE ALGEBRAIC RESIDUUM :
00508     // -----------------------------------
00509     double *initialAlgebraicResiduum;
00510     double        relaxationConstant;
00511 
00512 
00513     // BUTCHER-
00514     // TABLEAU:
00515     // ----------
00516     int      dim               ;  
00517     double **A                 ;  
00518     double  *b4                ;  
00519     double  *b5                ;  
00520     double  *c                 ;  
00523     // RK-STARTER:
00524     // -----------
00525     double   *eta4             ;  
00526     double   *eta5             ;  
00527     double   *eta4_            ;  
00528     double   *eta5_            ;  
00529     double ***k                ;  
00530     double ***k2               ;  
00531     double ***l                ;  
00532     double ***l2               ;  
00533     double   *iseed            ;  
00536     // BDF-METHOD:
00537     // -----------
00538     int      nstep             ; 
00539     double **psi               ; 
00540     double  *psi_              ; 
00541     double **gamma             ; 
00542     DMatrix   nablaY            ; 
00543     DMatrix   nablaY_           ; 
00544     DMatrix   phi               ; 
00545     DMatrix   delta             ; 
00546     DVector   c2                ; 
00549     DMatrix **M                 ; 
00550 //    std::vector< Eigen::HouseholderQR< DMatrix::Base > > qr;
00551     int     *M_index           ; 
00552     int      nOfM              ; 
00553     int      maxNM             ; 
00555     int     *nOfNewtonSteps    ; 
00556     double **eta               ; 
00557     double **eta2              ; 
00558     double   t                 ; 
00559     double  *x                 ; 
00561     double  *F                 ;
00562     double  *F2                ;
00563 
00564 
00565     // OTHERS:
00566     // -------
00567     double *initial_guess      ;  
00568     double  rel_time_scale     ;  
00571     // SENSITIVITIES:
00572     // --------------
00573     int        ndir            ;  
00575     DVector     fseed           ;  
00576     DVector     bseed           ;  
00578     DVector     fseed2          ;  
00579     DVector     bseed2          ;  
00581     double    *G               ;  
00582     double    *etaG            ;  
00583     DMatrix     nablaG          ;  
00584     DMatrix     phiG            ;  
00585     DMatrix     deltaG          ;  
00586     DVector     c2G             ;  
00588     double    *G2              ;  
00589     double    *G3              ;  
00590     double    *etaG2           ;  
00591     double    *etaG3           ;  
00592     DMatrix     nablaG2         ;  
00593     DMatrix     nablaG3         ;  
00594     DMatrix     phiG2           ;  
00595     DMatrix     phiG3           ;  
00596     DMatrix     deltaG2         ;  
00597     DMatrix     deltaG3         ;  
00598     DVector     c2G2            ;  
00599     DVector     c2G3            ;  
00601     double    *H               ;  
00602     double   **etaH            ;  
00603     DMatrix     nablaH          ;  
00604     DMatrix     nablaH_         ;  
00605     DMatrix     deltaH          ;  
00606     DVector     c2H             ;  
00607     double  ***kH              ;  
00608     double   **zH              ;  
00610     double    *H2              ;  
00611     double    *H3              ;  
00612     double   **etaH2           ;  
00613     double   **etaH3           ;  
00614     DMatrix     nablaH2         ;  
00615     DMatrix     nablaH3         ;  
00616     DMatrix     nablaH2_        ;  
00617     DMatrix     nablaH3_        ;  
00618     DMatrix     deltaH2         ;  
00619     DMatrix     deltaH3         ;  
00620     DVector     c2H2            ;  
00621     DVector     c2H3            ;  
00622     double  ***kH2             ;  
00623     double   **zH2             ;  
00624     double  ***kH3             ;  
00625     double   **zH3             ;  
00628     // STORAGE:
00629     // --------
00630     int maxAlloc               ;  
00634     // STATISTICS:
00635     // -----------
00636     RealClock jacComputation    ;
00637     RealClock jacDecomposition  ;
00638     RealClock correctorTime     ;
00639     int       nJacEvaluations   ;
00640 
00641 };
00642 
00643 
00644 CLOSE_NAMESPACE_ACADO
00645 
00646 
00647 
00648 #include <acado/integrator/integrator_bdf.ipp>
00649 
00650 
00651 #endif  // ACADO_TOOLKIT_INTEGRATOR_BDF_HPP
00652 
00653 // end of file.