34 namespace ublas = boost::numeric::ublas;
36 template<
typename T,
int>
52 enum { _E_order = O };
60 template<
typename Data,
80 for ( uint __i = 0; __i < M_grad.size(); ++__i )
82 M_grad[ __i ].resize( (
size_t )
_E_n );
83 M_hess[ __i ].resize( (
size_t )_E_n, (
size_t )_E_n );
86 typedef typename Data::value_type numerical_type;
87 typedef ublas::vector<double> numerical_gradient_type;
94 typename array_fixed<numerical_type,_E_f>::type,
96 typename array_fixed<numerical_type,_E_h>::type,
97 typename array_fixed<numerical_type,_E_ng>::type
103 typename array_fixed<numerical_gradient_type,_E_f>::type,
105 typename array_fixed<numerical_gradient_type,_E_h>::type,
106 typename array_fixed<numerical_gradient_type,_E_g>::type
112 typename array_fixed<numerical_hessian_type,_E_f>::type,
114 typename array_fixed<numerical_hessian_type,_E_h>::type,
115 typename array_fixed<numerical_hessian_type,_E_g>::type
129 numerical_type value(
int __i )
const
133 numerical_type & value(
int __i )
147 numerical_gradient_type & gradient(
int __i )
149 return M_grad[ __i ];
151 numerical_gradient_type gradient(
int __i )
const
153 return M_grad[ __i ];
156 numerical_type & gradient(
int __i,
int __j )
158 return M_grad[ __i]( __j );
160 numerical_type gradient(
int __i,
int __j )
const
162 return M_grad[ __i]( __j );
176 return M_hess[ __i ];
180 return M_hess[ __i ];
183 numerical_type & hessian(
int __i,
int __j,
int __k )
185 return M_hess[__i]( __j, __k );
187 numerical_type hessian(
int __i,
int __j,
int __k )
const
189 return M_hess[__i]( __j, __k );
203 template<
typename Data>
212 template<
int Diff = 0>
215 enum { diff_order = Diff };
216 typedef diff<diff_order> type;
226 template<
typename Data>
233 template<
int Diff = 0>
236 enum { diff_order = Diff };
238 typedef diff<diff_order> type;
250 template<
typename Data>
255 template<
int Diff = 0>
258 enum { diff_order = Diff };
259 typedef diff<diff_order> type;
273 enum { _E_dummy = N };
282 template<
typename Data>
309 typedef Data data_type;
329 typedef ublas::matrix<double> matrix_type;
333 typedef typename IF<
_E_g!=0,
336 typedef typename IF<
_E_h!=0,
367 template<
typename Defs>
385 return M_xlower( i );
391 return M_xupper( i );
405 return __x - M_xlower;
409 return M_xupper - __x;
418 return __x_definitions[__i][0] >=
value_type( .5 );
439 typedef typename mpl::if_<( Order == 0 ),
441 typename mpl::if_<Order == 1,
445 typedef typename mpl::if_<(
_E_g!=0 ),
449 typedef typename mpl::if_<
_E_h!=0,
461 __pt->f( __x, __adt );
463 to ( __adt, __fx, mpl::int_<Order>() );
465 void to ( ad_type
const& __adt,
f_type& __fx, mpl::int_<0> )
467 __fx.value( 0 ) = __adt.value();
469 void to ( ad_type
const& __adt,
f_type& __fx, mpl::int_<1> )
471 __fx.value( 0 ) = __adt.value();
472 __fx.gradient( 0 ) = __adt.grad();
474 void to ( ad_type
const& __adt,
f_type& __fx, mpl::int_<2> )
476 __fx.value( 0 ) = __adt.value();
477 __fx.gradient( 0 ) = __adt.grad();
478 __fx.hessian( 0 ) = __adt.hessian();
486 inequalities_array_type __adt;
488 __pt->g( __x, __adt );
489 to ( __adt, __gx, mpl::int_<Order>() );
492 for (
int __i = 0; __i < _E_n; ++__i )
494 __gx.value(
_E_g + __i ) = __pt->
x_l( __i ) - __x( __i );
495 __gx.value(
_E_g + _E_n + __i ) = __x( __i ) - __pt->
x_u( __i );
498 void to ( inequalities_array_type
const& __adt, g_type& __gx, mpl::int_<0> )
500 for (
int __i = 0; __i <
_E_g; ++__i )
502 __gx.value( __i ) = __adt[ __i ].value();
507 void to ( inequalities_array_type
const& __adt, g_type& __gx, mpl::int_<1> )
509 for (
int __i = 0; __i <
_E_g; ++__i )
511 __gx.value( __i ) = __adt[ __i ].value();
513 for (
int __j = 0; __j < _E_n; ++__j )
516 __gx.gradient( __i,__j ) = __adt[ __i ].grad( __j );
521 void to ( inequalities_array_type
const& __adt, g_type& __gx, mpl::int_<2> )
523 for (
int __i = 0; __i <
_E_g; ++__i )
525 __gx.value( __i ) = __adt[ __i ].value();
527 for (
int __j = 0; __j < _E_n; ++__j )
529 __gx.gradient( __i, __j ) = __adt[ __i ].grad( __j );
531 for (
int __k = 0; __k <= __j; ++__k )
533 __gx.hessian( __i, __k, __j ) =__adt[ __i ].hessian( __k, __j );
543 equalities_array_type __adt;
545 __pt->h( __x, __adt );
546 to ( __adt, __hx, mpl::int_<Order>() );
548 void to ( equalities_array_type
const& __adt, h_type& __hx, mpl::int_<0> )
550 for (
int __i = 0; __i <
_E_h; ++__i )
552 __hx.value( __i ) = __adt[ __i ].value();
555 void to ( equalities_array_type
const& __adt, h_type& __hx, mpl::int_<1> )
557 for (
int __i = 0; __i <
_E_h; ++__i )
559 __hx.value( __i ) = __adt[ __i ].value();
561 for (
int __j = 0; __j < _E_n; ++__j )
563 __hx.gradient( __i,__j ) = __adt[ __i ].grad( __j );
567 void to ( equalities_array_type
const& __adt, h_type& __hx, mpl::int_<2> )
569 for (
int __i = 0; __i <
_E_h; ++__i )
571 __hx.value( __i ) = __adt[ __i ].value();
573 for (
int __j = 0; __j < _E_n; ++__j )
575 __hx.gradient( __i, __j ) = __adt[ __i ].grad( __j );
577 for (
int __k = 0; __k <= __j; ++__k )
579 __hx.hessian( __i, __k, __j ) =__adt[ __i ].hessian( __k, __j );
594 template<
int Order,
typename Type>
614 void print_complete( std::string
const&,
vector_type const& _x )
const;
616 void print_bound_constraints()
const;
638 template<
typename Data>
646 M_xlower = ublas::scalar_vector<value_type>( M_xlower.size(), -inf );
647 M_xupper = ublas::scalar_vector<value_type>( M_xupper.size(), +inf );
651 template<
typename Data>
664 template<
typename Data>
669 template<
typename Data>
670 template<
typename Defs>
688 for (
int i = 0; i < _E_n; i++ )
690 __x_definitions[ i][0 ] = x_definitions[ i][0 ];
691 __x_definitions[ i][1 ] = x_definitions[ i][1 ];
692 __x_definitions[ i][2 ] = x_definitions[ i][2 ];
696 M_xlower( i ) = __x_definitions[ i][1 ];
697 M_xupper( i ) = __x_definitions[ i][2 ];
706 template<
typename Data>
712 for (
int i = 0; i < _E_n; i++ )
715 __x( i-offset ) = M_x0( i );
723 template<
typename Data>
729 for (
int i = 0; i < _E_n; i++ )
732 M_x0( i ) = __x( i-offset );
740 template<
typename Data>
744 for (
int i = 0; i < _E_n; i++ )
748 if ( __x_definitions[ i][1 ] <= -inf && __x_definitions[ i][2 ] >= inf )
751 if ( __x_definitions[ i][1 ] <= -inf && __x_definitions[ i][2 ] < inf )
752 M_x0( i ) = __x_definitions[ i][2 ] -
value_type( 10 );
754 if ( __x_definitions[ i][1 ] > -inf && __x_definitions[ i][2 ] >= inf )
755 M_x0( i ) = __x_definitions[ i][1 ] + value_type( 10 );
757 if ( __x_definitions[ i][1 ] > -inf && __x_definitions[ i][2 ] < inf )
758 M_x0( i ) = 0.5 * ( __x_definitions[ i][1 ] + __x_definitions[ i][2 ] );
762 M_x0( i ) = __x_definitions[ i][1 ];
767 template<
typename Data>
771 for (
int i = 0; i < _E_n; i++ )
773 if ( __custom( i ) < __x_definitions[ i][1 ] )
775 std::ostringstream __ex;
776 __ex <<
"Initialization Error: initial x^0[ " << i <<
" ] < xlower[ " << i <<
" ]\n"
777 <<
" * x^0[ " << i <<
" ] = " << __custom( i ) <<
"\n"
778 <<
" * xlower[ " << i <<
" ] = " << __x_definitions[ i][1 ] <<
"\n";
780 throw std::invalid_argument( __ex.str() );
783 if ( __custom( i ) > __x_definitions[ i][2 ] )
785 std::ostringstream __ex;
786 __ex <<
"Initialization Error: initial x^0[ " << i <<
" ] > xupper[ " << i <<
" ]\n"
787 <<
" * x^0[ " << i <<
" ] = " << __custom( i ) <<
"\n"
788 <<
" * xupper[ " << i <<
" ] = " << __x_definitions[ i][2 ] <<
"\n";
790 throw std::invalid_argument( __ex.str() );
793 M_x0( i ) = __custom( i );
803 template<
typename Data>
810 for (
int i = 0; i < _E_n; i++ )
813 std::cout <<
" x [ " << i <<
" ] = " << _x( i-offset );
817 std::cout <<
" x [ " << i <<
" ] = " << __x_definitions[ i][1 ] <<
" ( variable \"off\" )";
825 template<
typename Data>
830 for (
int i = 0; i < _E_n; i++ )
834 if ( __x_definitions[ i][1 ] == -inf )
835 std::cerr << std::setw( 8 ) << std::right <<
" -Inf" <<
" <= ";
837 else if ( abs( __x_definitions[ i][1 ] ) < 1e4 )
838 std::cerr << std::setw( 8 ) << std::right << __x_definitions[ i][1 ] <<
" <= ";
841 std::cerr << std::setw( 8 ) << std::right << __x_definitions[ i][1 ] <<
" <= ";
843 std::cerr <<
"x[" << i <<
"] <= ";
845 if ( __x_definitions[ i][2 ] == inf )
846 std::cerr << std::left <<
"+Inf";
848 else if ( fabs( __x_definitions[ i][2 ] ) < 1e4 )
849 std::cerr << std::left << __x_definitions[ i][2 ];
852 std::cerr << std::left << __x_definitions[ i][2 ];
857 std::cerr <<
"x[" << i <<
" ] := " << __x_definitions[ i][1 ];
864 template<
typename Data>
866 problem<Data>::print_complete( std::string
const& __mesg, vector_type
const& _x )
const
868 std::cout << __mesg <<
"\n";
871 print_bound_constraints();
874 template<
typename Data>
876 problem<Data>::print_complete( vector_type
const& _x,
877 vector_type
const& _s )
const
879 print_bound_constraints();
884 this->evaluate ( _x, __fx, diff_order<2>() );
886 std::cout <<
"fx.gradient:\n"
887 << __fx.gradient( 0 ) <<
"\n";
889 std::cout <<
"fx.hessian:\n"
890 << __fx.hessian( 0 ) <<
"\n";
894 this->evaluate ( _x, __gx, diff_order<2>() );
896 for (
int m = 0; m < _E_g; m++ )
898 std::cout <<
"gx(" << m <<
"):\n"
899 << __gx.value( m ) <<
"\n";
901 std::cout <<
"gx(" << m <<
").gradient:\n"
902 << __gx.gradient( m ) <<
"\n";
904 std::cout <<
"gx(" << m <<
").hessian:\n"
905 << __gx.hessian( m ) <<
"\n";
909 this->evaluate ( _x, __hx, diff_order<2>() );
911 for (
int m = 0; m < _E_h; m++ )
913 std::cout <<
"hx(" << m <<
"):\n"
914 << __hx.value( m ) <<
"\n";
916 std::cout <<
"hx(" << m <<
").gradient:\n"
917 << __hx.gradient( m ) <<
"\n";
919 std::cout <<
"hx(" << m <<
").hessian:\n"
920 << __hx.hessian( m ) <<
"\n";
925 template<
typename Data>
927 problem<Data>::print_stationary_x( vector_type
const& _x,
928 vector_type
const& _lambda_l,
929 vector_type
const& _lambda_u )
const
932 vector_type __sum( _E_n );
935 this->print_bound_constraints();
937 std::cerr.precision ( 5 );
939 std::cerr <<
"\n\nSolution:\n";
942 for (
int i = 0; i < _E_n; i++ )
946 std::cerr <<
" x*[ " << i <<
" ] = " << std::setw( 10 ) << _x( i-offset );
948 if ( std::abs( _lambda_l( i-offset ) ) > 1e-6 )
949 std::cerr <<
" ( lower bound active )";
951 if ( std::abs( _lambda_u( i-offset ) ) > 1e-6 )
952 std::cerr <<
" ( upper bound active )";
957 std::cerr <<
" x*[ " << i <<
" ] = " << std::setw( 10 ) << __x_definitions[ i][1 ];
958 std::cerr <<
" ( variable \"off\" ) ";
966 this->evaluate( _x, __f_x, diff_order<2>() );
967 std::cerr <<
"f(x*) = " << __f_x.value( 0 ) <<
"\n"
968 <<
"f'(x*) = " << __f_x.gradient( 0 ) <<
"\n"
969 <<
"f''(x*) = " << __f_x.hessian( 0 ) <<
"\n";
971 std::cerr <<
"\nStationarity conditions (LS):\n"
972 <<
" \\nabla f - \\lambda_l + \\lambda_u = residual\n"
973 <<
"----------------------------------------------------------\n";
977 for (
int i = 0; i < _E_n; i++ )
980 __sum( i-offset ) = __f_x.gradient( 0, i-offset ) - _lambda_l( i-offset ) + _lambda_u( i-offset );
982 std::cerr <<
" " << std::setw( 10 ) << __f_x.gradient( 0, i-offset ) <<
" - "
983 <<
" " << std::setw( 10 ) << _lambda_l( i-offset ) <<
" + "
984 <<
" " << std::setw( 10 ) << _lambda_u( i-offset ) <<
" = "
985 <<
" " << std::setw( 10 ) << __sum( i-offset ) <<
"\n";
991 std::cerr <<
" x x x x x x x x x x x x x x x x x x x x x x x x x x x\n";
995 std::cerr <<
"\n =========================================================";
996 std::cerr <<
"\n # Norm of residuals ------------------------> " ;
997 std::cerr <<
" " << std::setprecision( 2 ) << std::setw( 8 ) << norm_2( __sum ) <<
" #";
998 std::cerr <<
"\n =========================================================\n";
1000 std::cerr <<
"\nComplimentarity conditions (LS):\n"
1001 <<
" \\lambda_l * ( l - x ) = 0.000000 | \\lambda_u * ( x - u ) = 0.00000\n"
1002 <<
"----------------------------------------------|----------------------------------------------\n";
1006 for (
int i = 0; i < _E_n; i++ )
1009 __prod = _lambda_l( i-offset ) * ( __x_definitions[ i][1 ] - _x( i-offset ) );
1011 std::cerr <<
" " << std::setw( 9 ) << _lambda_l( i-offset ) <<
" * "
1012 <<
" " << std::setw( 15 ) << __x_definitions[ i][1 ] - _x( i-offset ) <<
" = "
1013 <<
" " << std::setw( 9 ) << __prod <<
" |";
1015 __prod = _lambda_u( i-offset ) * ( _x( i-offset ) - __x_definitions[ i][2 ] );
1017 std::cerr <<
" " << std::setw( 9 ) << _lambda_u( i-offset ) <<
" * "
1018 <<
" " << std::setw( 15 ) << _x( i-offset )-__x_definitions[ i][2 ] <<
" = "
1019 <<
" " << std::setw( 9 ) << __prod <<
"\n";
1025 std::cerr <<
" x x x x x x x x x x x x x x x x x x x x x x |"
1026 <<
" x x x x x x x x x x x x x x x x x x x x x x\n";
1032 template<
typename Data>
1034 problem<Data>::print_stationaryN_x( vector_type
const& _x,
1035 vector_type
const& _s,
1036 vector_type
const& _lambda_h,
1037 vector_type
const& _lambda_g )
const
1039 std::cerr <<
"\n\nVariables and (variable) constraints:\n";
1041 this->print_bound_constraints();
1043 std::cerr <<
"\n\nSolution:\n\n";
1045 std::cerr.precision ( 5 );
1049 for (
int i = 0; i < _E_n; i++ )
1053 std::cerr <<
" x*[ " << i <<
" ] = " << std::setw( 10 ) << _x( i-offset );
1055 if ( _x( i-offset ) < M_xlower( i-offset ) ||
1056 _x( i-offset ) > M_xupper( i-offset ) )
1058 std::cerr <<
" <- ( violated ! )";
1063 value_type lam_l = _lambda_g( _E_g+i-offset );
1064 value_type lam_u = _lambda_g( _E_g+_E_n+i-offset );
1066 if ( std::abs( lam_l ) > 1e-4 )
1067 std::cerr <<
" ( lower bound active, lambda_l " << std::setw( 10 ) << lam_l;
1069 if ( std::abs( lam_u ) > 1e-4 )
1070 std::cerr <<
" ( upper bound active, lambda_u " << std::setw( 10 ) << lam_u;
1076 std::cerr <<
" x*[ " << i <<
" ] = " << std::setw( 10 ) << __x_definitions[ i][1 ];
1077 std::cerr <<
" ( variable \"off\" ) ";
1086 this->evaluate ( _x, __fx, diff_order<1>() );
1088 std::cerr <<
" f(x*) = " << std::setw( 10 ) << __fx.value( 0 ) <<
"\n";
1092 this->evaluate ( _x, __gx, diff_order<0>() );
1094 for (
int m = 0; m < _E_g; m++ )
1096 if ( __gx.value( m ) < 0 )
1098 if ( std::abs( std::fabs( _lambda_g( m ) ) ) > 1e-3 )
1100 std::cerr <<
" g[" << m <<
"](x*) = " << std::setw( 10 ) << __gx.value( m )
1101 <<
" ( satisfied - active, lambda = " << std::setw( 10 ) << _lambda_g( m );
1106 std::cerr <<
" g[" << m <<
"](x*) = " << std::setw( 10 ) << __gx.value( m )
1107 <<
" ( satisfied - inactive )";
1114 std::cerr <<
" g[" << m <<
"](x*) = " << std::setw( 10 ) << __gx.value( m )
1115 <<
" <- ( violated ! )";
1123 this->evaluate ( _x, __hx, diff_order<0>() );
1125 for (
int m = 0; m < _E_h; m++ )
1127 std::cerr <<
" h[" << m <<
"](x*) = " << std::setw( 10 ) << __hx.value( m );
1129 if ( std::abs( __hx.value( m ) < 1e-3 ) )
1130 std::cerr <<
" ( satisfied )";
1133 std::cerr <<
" <- ( violated ! )";
number of objective functions
Definition: problem.hpp:293
super::ad_1_type ad_1_type
automatic differentiation type of order one
Definition: problem.hpp:316
Definition: solverlinear.hpp:33
void custom_initialize_x0(vector_type &__custom)
use custom starting point
Definition: problem.hpp:769
super::ad_0_type::value_type value_type
numerical type
Definition: problem.hpp:322
number of objective functions
Definition: problem.hpp:68
IF< Id==FUNCTIONAL, typename array_fixed< numerical_type, _E_f >::type, typename IF< Id==EQUALITIES, typename array_fixed< numerical_type, _E_h >::type, typename array_fixed< numerical_type, _E_ng >::type >::Result >::Result value_type
value type for the functional
Definition: problem.hpp:99
value_type x_l(unsigned int i) const
Definition: problem.hpp:383
problem()
default constructor: initialize with infinite bounds
Definition: problem.hpp:639
void initialize_x0()
initialize the starting point x0
Definition: problem.hpp:742
ublas::vector< double > vector_type
vector type
Definition: problem.hpp:325
Data super
super class for the problem.
Definition: problem.hpp:308
number of equality constraints
Definition: problem.hpp:295
super::ad_0_type ad_0_type
automatic differentiation type of order one
Definition: problem.hpp:313
functional identifier
Definition: problem.hpp:44
number of multiplers
Definition: problem.hpp:297
super::ad_2_type ad_2_type
automatic differentiation type of order two
Definition: problem.hpp:319
~problem()
destructor
Definition: problem.hpp:665
Optimization problem specifications.
Definition: problem.hpp:283
number of inequalities
Definition: problem.hpp:69
number of entries in the matrix
Definition: problem.hpp:72
define the inequalities type
Definition: problem.hpp:227
IF< Id==FUNCTIONAL, typename array_fixed< numerical_gradient_type, _E_f >::type, typename IF< Id==EQUALITIES, typename array_fixed< numerical_gradient_type, _E_h >::type, typename array_fixed< numerical_gradient_type, _E_g >::type >::Result >::Result gradient_type
value type for the functional gradient
Definition: problem.hpp:108
void evaluate(vector_type const &__x, Type &__fx, diff_order< Order >) const
evaluate components of the problem Type can be #) the objective functions #) the inequality contraint...
Definition: problem.hpp:595
void copy_x0_to_x(vector_type &_x)
copy initialization x0 into x
Definition: problem.hpp:708
compute the value of the functionals, equalities and inequalities associated with the problem ...
Definition: problem.hpp:431
Definition: problem.hpp:62
define the equalities type
Definition: problem.hpp:251
defines the functional type
Definition: problem.hpp:204
number of control variables
Definition: problem.hpp:292
dummy data type
Definition: problem.hpp:271
number of inequality constraints
Definition: problem.hpp:294
equalities identifier
Definition: problem.hpp:45
number of control variables
Definition: problem.hpp:67
void print(vector_type const &_x) const
prints current value of x
Definition: problem.hpp:805
void define_problem(Defs)
define the bounds of the problem
Definition: problem.hpp:672
void copy_x_to_x0(vector_type const &_x)
copy x to initialization x0
Definition: problem.hpp:725
ublas::symmetric_matrix< double, ublas::upper > numerical_hessian_type
matrix type
Definition: problem.hpp:89
int nA() const
Definition: problem.hpp:377
bool isOn(int __i) const
Definition: problem.hpp:416
ublas::symmetric_matrix< double, ublas::upper > symmetric_matrix_type
matrix type
Definition: problem.hpp:328
size of the multipliers matrix
Definition: problem.hpp:298
size of informations stored for inequalities
Definition: problem.hpp:70
number of equalities
Definition: problem.hpp:71
IF< Id==FUNCTIONAL, typename array_fixed< numerical_hessian_type, _E_f >::type, typename IF< Id==EQUALITIES, typename array_fixed< numerical_hessian_type, _E_h >::type, typename array_fixed< numerical_hessian_type, _E_g >::type >::Result >::Result hessian_type
value type for the functional hessian: in packed format
Definition: problem.hpp:117
size of the matrix
Definition: problem.hpp:296
inequalities identifier
Definition: problem.hpp:46
value_type x_u(unsigned int i) const
Definition: problem.hpp:389
int n() const
Definition: problem.hpp:371