21 #include <yarp/math/Math.h>
26 #include <IpIpoptApplication.hpp>
29 using namespace yarp::os;
30 using namespace yarp::sig;
31 using namespace yarp::math;
39 namespace optimization
45 Matrix
H=euler2dcm(
x.subVector(3,5));
46 H(0,3)=
x[0];
H(1,3)=
x[1];
H(2,3)=
x[2];
56 for (
size_t i=0; i<_x.length(); i++)
67 const deque<Vector> &
p0;
68 const deque<Vector> &
p1;
78 const deque<Vector> &_p1,
79 const Vector &_min,
const Vector &_max) :
88 virtual void set_x0(
const Vector &x0)
90 size_t len=
std::min(this->x0.length(),x0.length());
91 for (
size_t i=0; i<len; i++)
102 bool get_nlp_info(Ipopt::Index &
n, Ipopt::Index &m, Ipopt::Index &nnz_jac_g,
103 Ipopt::Index &nnz_h_lag, IndexStyleEnum &index_style)
106 m=nnz_jac_g=nnz_h_lag=0;
107 index_style=TNLP::C_STYLE;
114 Ipopt::Index m, Ipopt::Number *g_l, Ipopt::Number *g_u)
116 for (Ipopt::Index i=0; i<
n; i++)
127 bool init_z, Ipopt::Number *z_L, Ipopt::Number *z_U,
128 Ipopt::Index m,
bool init_lambda, Ipopt::Number *lambda)
130 for (Ipopt::Index i=0; i<
n; i++)
137 bool eval_f(Ipopt::Index
n,
const Ipopt::Number *
x,
bool new_x,
138 Ipopt::Number &obj_value)
145 for (
size_t i=0; i<p0.size(); i++)
146 obj_value+=
norm2(p1[i]-
H*p0[i]);
148 obj_value/=p0.size();
156 Ipopt::Number *grad_f)
158 double ca=cos(
x[3]);
double sa=sin(
x[3]);
159 double cb=cos(
x[4]);
double sb=sin(
x[4]);
160 double cg=cos(
x[5]);
double sg=sin(
x[5]);
163 Rza(0,0)=ca; Rza(1,1)=ca; Rza(1,0)=sa; Rza(0,1)=-sa;
164 Matrix dRza=
zeros(4,4);
165 dRza(0,0)=-sa; dRza(1,1)=-sa; dRza(1,0)=ca; dRza(0,1)=-ca;
168 Rzg(0,0)=cg; Rzg(1,1)=cg; Rzg(1,0)=sg; Rzg(0,1)=-sg;
169 Matrix dRzg=
zeros(4,4);
170 dRzg(0,0)=-sg; dRzg(1,1)=-sg; dRzg(1,0)=cg; dRzg(0,1)=-cg;
173 Ryb(0,0)=cb; Ryb(2,2)=cb; Ryb(2,0)=-sb; Ryb(0,2)=sb;
174 Matrix dRyb=
zeros(4,4);
175 dRyb(0,0)=-sb; dRyb(2,2)=-sb; dRyb(2,0)=-cb; dRyb(0,2)=cb;
177 Matrix dHdx0=
zeros(4,4); dHdx0(0,3)=1.0;
178 Matrix dHdx1=
zeros(4,4); dHdx1(1,3)=1.0;
179 Matrix dHdx2=
zeros(4,4); dHdx2(2,3)=1.0;
180 Matrix dHdx3=dRza*Ryb*Rzg;
181 Matrix dHdx4=Rza*dRyb*Rzg;
182 Matrix dHdx5=Rza*Ryb*dRzg;
185 grad_f[0]=grad_f[1]=grad_f[2]=0.0;
186 grad_f[3]=grad_f[4]=grad_f[5]=0.0;
189 for (
size_t i=0; i<p0.size(); i++)
191 Vector d=p1[i]-
H*p0[i];
192 grad_f[0]-=2.0*
dot(d,(dHdx0*p0[i]));
193 grad_f[1]-=2.0*
dot(d,(dHdx1*p0[i]));
194 grad_f[2]-=2.0*
dot(d,(dHdx2*p0[i]));
195 grad_f[3]-=2.0*
dot(d,(dHdx3*p0[i]));
196 grad_f[4]-=2.0*
dot(d,(dHdx4*p0[i]));
197 grad_f[5]-=2.0*
dot(d,(dHdx5*p0[i]));
200 for (Ipopt::Index i=0; i<
n; i++)
201 grad_f[i]/=p0.size();
208 bool eval_g(Ipopt::Index
n,
const Ipopt::Number *
x,
bool new_x,
209 Ipopt::Index m, Ipopt::Number *g)
216 Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index *iRow,
217 Ipopt::Index *jCol, Ipopt::Number *values)
223 bool eval_h(Ipopt::Index
n,
const Ipopt::Number *
x,
bool new_x,
224 Ipopt::Number obj_factor, Ipopt::Index m,
const Ipopt::Number *lambda,
225 bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index *iRow,
226 Ipopt::Index *jCol, Ipopt::Number *values)
233 const Ipopt::Number *
x,
const Ipopt::Number *z_L,
234 const Ipopt::Number *z_U, Ipopt::Index m,
235 const Ipopt::Number *g,
const Ipopt::Number *lambda,
236 Ipopt::Number obj_value,
const Ipopt::IpoptData *ip_data,
237 Ipopt::IpoptCalculatedQuantities *ip_cq)
240 for (Ipopt::Index i=0; i<
n; i++)
252 const deque<Vector> &_p1,
253 const Vector &_min,
const Vector &_max) :
257 bool get_nlp_info(Ipopt::Index &
n, Ipopt::Index &m, Ipopt::Index &nnz_jac_g,
258 Ipopt::Index &nnz_h_lag, IndexStyleEnum &index_style)
260 CalibReferenceWithMatchedPointsNLP::get_nlp_info(
n,m,nnz_jac_g,nnz_h_lag,index_style);
267 bool eval_f(Ipopt::Index
n,
const Ipopt::Number *
x,
bool new_x,
268 Ipopt::Number &obj_value)
280 for (
size_t i=0; i<p0.size(); i++)
281 obj_value+=
norm2(p1[i]-s*(
H*p0[i]));
283 obj_value/=p0.size();
291 Ipopt::Number *grad_f)
293 double ca=cos(
x[3]);
double sa=sin(
x[3]);
294 double cb=cos(
x[4]);
double sb=sin(
x[4]);
295 double cg=cos(
x[5]);
double sg=sin(
x[5]);
298 Rza(0,0)=ca; Rza(1,1)=ca; Rza(1,0)=sa; Rza(0,1)=-sa;
299 Matrix dRza=
zeros(4,4);
300 dRza(0,0)=-sa; dRza(1,1)=-sa; dRza(1,0)=ca; dRza(0,1)=-ca;
303 Rzg(0,0)=cg; Rzg(1,1)=cg; Rzg(1,0)=sg; Rzg(0,1)=-sg;
304 Matrix dRzg=
zeros(4,4);
305 dRzg(0,0)=-sg; dRzg(1,1)=-sg; dRzg(1,0)=cg; dRzg(0,1)=-cg;
308 Ryb(0,0)=cb; Ryb(2,2)=cb; Ryb(2,0)=-sb; Ryb(0,2)=sb;
309 Matrix dRyb=
zeros(4,4);
310 dRyb(0,0)=-sb; dRyb(2,2)=-sb; dRyb(2,0)=-cb; dRyb(0,2)=cb;
312 Matrix dHdx0=
zeros(4,4); dHdx0(0,3)=1.0;
313 Matrix dHdx1=
zeros(4,4); dHdx1(1,3)=1.0;
314 Matrix dHdx2=
zeros(4,4); dHdx2(2,3)=1.0;
315 Matrix dHdx3=dRza*Ryb*Rzg;
316 Matrix dHdx4=Rza*dRyb*Rzg;
317 Matrix dHdx5=Rza*Ryb*dRzg;
320 Matrix dHdx6=
zeros(4,4); dHdx6(0,0)=1.0; dHdx6*=
H;
321 Matrix dHdx7=
zeros(4,4); dHdx7(1,1)=1.0; dHdx7*=
H;
322 Matrix dHdx8=
zeros(4,4); dHdx8(2,2)=1.0; dHdx8*=
H;
330 grad_f[0]=grad_f[1]=grad_f[2]=0.0;
331 grad_f[3]=grad_f[4]=grad_f[5]=0.0;
332 grad_f[6]=grad_f[7]=grad_f[8]=0.0;
335 for (
size_t i=0; i<p0.size(); i++)
337 Vector d=p1[i]-s*(
H*p0[i]);
338 grad_f[0]-=2.0*
dot(d,(dHdx0*p0[i]));
339 grad_f[1]-=2.0*
dot(d,(dHdx1*p0[i]));
340 grad_f[2]-=2.0*
dot(d,(dHdx2*p0[i]));
341 grad_f[3]-=2.0*
dot(d,(dHdx3*p0[i]));
342 grad_f[4]-=2.0*
dot(d,(dHdx4*p0[i]));
343 grad_f[5]-=2.0*
dot(d,(dHdx5*p0[i]));
344 grad_f[6]-=2.0*
dot(d,(dHdx6*p0[i]));
345 grad_f[7]-=2.0*
dot(d,(dHdx7*p0[i]));
346 grad_f[8]-=2.0*
dot(d,(dHdx8*p0[i]));
349 for (Ipopt::Index i=0; i<
n; i++)
350 grad_f[i]/=p0.size();
364 const deque<Vector> &_p1,
365 const Vector &_min,
const Vector &_max) :
369 bool get_nlp_info(Ipopt::Index &
n, Ipopt::Index &m, Ipopt::Index &nnz_jac_g,
370 Ipopt::Index &nnz_h_lag, IndexStyleEnum &index_style)
372 CalibReferenceWithMatchedPointsNLP::get_nlp_info(
n,m,nnz_jac_g,nnz_h_lag,index_style);
379 bool eval_f(Ipopt::Index
n,
const Ipopt::Number *
x,
bool new_x,
380 Ipopt::Number &obj_value)
390 for (
size_t i=0; i<p0.size(); i++)
391 obj_value+=
norm2(p1[i]-s*(
H*p0[i]));
393 obj_value/=p0.size();
401 Ipopt::Number *grad_f)
403 double ca=cos(
x[3]);
double sa=sin(
x[3]);
404 double cb=cos(
x[4]);
double sb=sin(
x[4]);
405 double cg=cos(
x[5]);
double sg=sin(
x[5]);
408 Rza(0,0)=ca; Rza(1,1)=ca; Rza(1,0)=sa; Rza(0,1)=-sa;
409 Matrix dRza=
zeros(4,4);
410 dRza(0,0)=-sa; dRza(1,1)=-sa; dRza(1,0)=ca; dRza(0,1)=-ca;
413 Rzg(0,0)=cg; Rzg(1,1)=cg; Rzg(1,0)=sg; Rzg(0,1)=-sg;
414 Matrix dRzg=
zeros(4,4);
415 dRzg(0,0)=-sg; dRzg(1,1)=-sg; dRzg(1,0)=cg; dRzg(0,1)=-cg;
418 Ryb(0,0)=cb; Ryb(2,2)=cb; Ryb(2,0)=-sb; Ryb(0,2)=sb;
419 Matrix dRyb=
zeros(4,4);
420 dRyb(0,0)=-sb; dRyb(2,2)=-sb; dRyb(2,0)=-cb; dRyb(0,2)=cb;
422 Matrix dHdx0=
zeros(4,4); dHdx0(0,3)=1.0;
423 Matrix dHdx1=
zeros(4,4); dHdx1(1,3)=1.0;
424 Matrix dHdx2=
zeros(4,4); dHdx2(2,3)=1.0;
425 Matrix dHdx3=dRza*Ryb*Rzg;
426 Matrix dHdx4=Rza*dRyb*Rzg;
427 Matrix dHdx5=Rza*Ryb*dRzg;
430 Matrix dHdx6=
zeros(4,4); dHdx6(0,0)=1.0; dHdx6*=
H;
436 grad_f[0]=grad_f[1]=grad_f[2]=0.0;
437 grad_f[3]=grad_f[4]=grad_f[5]=0.0;
441 for (
size_t i=0; i<p0.size(); i++)
443 Vector d=p1[i]-s*(
H*p0[i]);
444 grad_f[0]-=2.0*
dot(d,(dHdx0*p0[i]));
445 grad_f[1]-=2.0*
dot(d,(dHdx1*p0[i]));
446 grad_f[2]-=2.0*
dot(d,(dHdx2*p0[i]));
447 grad_f[3]-=2.0*
dot(d,(dHdx3*p0[i]));
448 grad_f[4]-=2.0*
dot(d,(dHdx4*p0[i]));
449 grad_f[5]-=2.0*
dot(d,(dHdx5*p0[i]));
450 grad_f[6]-=2.0*
dot(d,(dHdx6*p0[i]));
453 for (Ipopt::Index i=0; i<
n; i++)
454 grad_f[i]/=p0.size();
467 CalibReferenceWithMatchedPoints::CalibReferenceWithMatchedPoints()
472 min.resize(6);
max.resize(6);
480 min_s.resize(3); max_s.resize(3);
481 min_s[0]=0.1; max_s[0]=10.0;
482 min_s[1]=0.1; max_s[1]=10.0;
483 min_s[2]=0.1; max_s[2]=10.0;
485 min_s_scalar=0.1; max_s_scalar=10.0;
494 void CalibReferenceWithMatchedPoints::setBounds(
const Matrix &
min,
497 if ((
min.rows()<3) || (
min.cols()<3) ||
498 (
max.rows()<3) || (
max.cols()<3))
501 this->min[0]=
min(0,3); this->max[0]=
max(0,3);
502 this->min[1]=
min(1,3); this->max[1]=
max(1,3);
503 this->min[2]=
min(2,3); this->max[2]=
max(2,3);
505 this->min[3]=
min(0,0); this->max[3]=
max(0,0);
506 this->min[4]=
min(0,1); this->max[4]=
max(0,1);
507 this->min[5]=
min(0,2); this->max[5]=
max(0,2);
512 void CalibReferenceWithMatchedPoints::setBounds(
const Vector &
min,
515 size_t len_min=
std::min(this->min.length(),
min.length());
516 size_t len_max=
std::min(this->max.length(),
max.length());
518 for (
size_t i=0; i<len_min; i++)
521 for (
size_t i=0; i<len_max; i++)
527 void CalibReferenceWithMatchedPoints::setScalingBounds(
const Vector &
min,
530 size_t len_min=
std::min(this->min_s.length(),
min.length());
531 size_t len_max=
std::min(this->max_s.length(),
max.length());
533 for (
size_t i=0; i<len_min; i++)
534 this->min_s[i]=
min[i];
536 for (
size_t i=0; i<len_max; i++)
537 this->max_s[i]=
max[i];
542 void CalibReferenceWithMatchedPoints::setScalingBounds(
const double min,
551 double CalibReferenceWithMatchedPoints::evalError(
const Matrix &
H)
556 for (
size_t i=0; i<p0.size(); i++)
567 bool CalibReferenceWithMatchedPoints::addPoints(
const Vector &p0,
570 if ((p0.length()>=3) && (p1.length()>=3))
572 Vector _p0=p0.subVector(0,2); _p0.push_back(1.0);
573 Vector _p1=p1.subVector(0,2); _p1.push_back(1.0);
575 this->p0.push_back(_p0);
576 this->p1.push_back(_p1);
586 void CalibReferenceWithMatchedPoints::getPoints(deque<Vector> &p0,
587 deque<Vector> &p1)
const
595 void CalibReferenceWithMatchedPoints::clearPoints()
603 bool CalibReferenceWithMatchedPoints::setInitialGuess(
const Matrix &
H)
605 if ((
H.rows()>=4) && (
H.cols()>=4))
607 Vector euler=dcm2euler(
H);
608 x0[0]=
H(0,3); x0[1]=
H(1,3); x0[2]=
H(2,3);
609 x0[3]=euler[0]; x0[4]=euler[1]; x0[5]=euler[2];
619 bool CalibReferenceWithMatchedPoints::setScalingInitialGuess(
const Vector &s)
621 if (s.length()>=s0.length())
623 s0=s.subVector(0,(
unsigned int)s0.length()-1);
632 bool CalibReferenceWithMatchedPoints::setScalingInitialGuess(
const double s)
640 bool CalibReferenceWithMatchedPoints::setCalibrationOptions(
const Property &options)
642 if (options.check(
"max_iter"))
643 max_iter=options.find(
"max_iter").asInt32();
645 if (options.check(
"tol"))
646 tol=options.find(
"tol").asFloat64();
657 Ipopt::SmartPtr<Ipopt::IpoptApplication>
app=
new Ipopt::IpoptApplication;
658 app->Options()->SetNumericValue(
"tol",tol);
659 app->Options()->SetIntegerValue(
"acceptable_iter",0);
660 app->Options()->SetStringValue(
"mu_strategy",
"adaptive");
661 app->Options()->SetIntegerValue(
"max_iter",max_iter);
662 app->Options()->SetStringValue(
"nlp_scaling_method",
"gradient-based");
663 app->Options()->SetStringValue(
"hessian_approximation",
"limited-memory");
664 app->Options()->SetIntegerValue(
"print_level",0);
665 app->Options()->SetStringValue(
"derivative_test",
"none");
671 Ipopt::ApplicationReturnStatus status=
app->OptimizeTNLP(GetRawPtr(nlp));
673 Vector
x=nlp->get_result();
677 return (status==Ipopt::Solve_Succeeded);
690 Ipopt::SmartPtr<Ipopt::IpoptApplication>
app=
new Ipopt::IpoptApplication;
691 app->Options()->SetNumericValue(
"tol",tol);
692 app->Options()->SetIntegerValue(
"acceptable_iter",0);
693 app->Options()->SetStringValue(
"mu_strategy",
"adaptive");
694 app->Options()->SetIntegerValue(
"max_iter",max_iter);
695 app->Options()->SetStringValue(
"nlp_scaling_method",
"gradient-based");
696 app->Options()->SetStringValue(
"hessian_approximation",
"limited-memory");
697 app->Options()->SetIntegerValue(
"print_level",0);
698 app->Options()->SetStringValue(
"derivative_test",
"none");
703 nlp->set_x0(cat(x0,s0));
704 Ipopt::ApplicationReturnStatus status=
app->OptimizeTNLP(GetRawPtr(nlp));
706 Vector
x=nlp->get_result();
710 S(0,0)=s[0]; S(1,1)=s[1]; S(2,2)=s[2];
713 return (status==Ipopt::Solve_Succeeded);
726 Ipopt::SmartPtr<Ipopt::IpoptApplication>
app=
new Ipopt::IpoptApplication;
727 app->Options()->SetNumericValue(
"tol",tol);
728 app->Options()->SetIntegerValue(
"acceptable_iter",0);
729 app->Options()->SetStringValue(
"mu_strategy",
"adaptive");
730 app->Options()->SetIntegerValue(
"max_iter",max_iter);
731 app->Options()->SetStringValue(
"nlp_scaling_method",
"gradient-based");
732 app->Options()->SetStringValue(
"hessian_approximation",
"limited-memory");
733 app->Options()->SetIntegerValue(
"print_level",0);
734 app->Options()->SetStringValue(
"derivative_test",
"none");
739 nlp->set_x0(cat(x0,s0_scalar));
740 Ipopt::ApplicationReturnStatus status=
app->OptimizeTNLP(GetRawPtr(nlp));
742 Vector
x=nlp->get_result();
746 S(0,0)=S(1,1)=S(2,2)=s;
749 return (status==Ipopt::Solve_Succeeded);
bool get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number *x, bool init_z, Ipopt::Number *z_L, Ipopt::Number *z_U, Ipopt::Index m, bool init_lambda, Ipopt::Number *lambda)
bool get_nlp_info(Ipopt::Index &n, Ipopt::Index &m, Ipopt::Index &nnz_jac_g, Ipopt::Index &nnz_h_lag, IndexStyleEnum &index_style)
CalibReferenceWithMatchedPointsNLP(const deque< Vector > &_p0, const deque< Vector > &_p1, const Vector &_min, const Vector &_max)
bool eval_g(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Index m, Ipopt::Number *g)
bool eval_h(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number *lambda, bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index *iRow, Ipopt::Index *jCol, Ipopt::Number *values)
virtual void set_x0(const Vector &x0)
void finalize_solution(Ipopt::SolverReturn status, Ipopt::Index n, const Ipopt::Number *x, const Ipopt::Number *z_L, const Ipopt::Number *z_U, Ipopt::Index m, const Ipopt::Number *g, const Ipopt::Number *lambda, Ipopt::Number obj_value, const Ipopt::IpoptData *ip_data, Ipopt::IpoptCalculatedQuantities *ip_cq)
const deque< Vector > & p1
bool eval_grad_f(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Number *grad_f)
virtual Vector get_result() const
bool eval_jac_g(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index *iRow, Ipopt::Index *jCol, Ipopt::Number *values)
const deque< Vector > & p0
bool eval_f(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Number &obj_value)
bool get_bounds_info(Ipopt::Index n, Ipopt::Number *x_l, Ipopt::Number *x_u, Ipopt::Index m, Ipopt::Number *g_l, Ipopt::Number *g_u)
bool eval_grad_f(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Number *grad_f)
CalibReferenceWithScalarScaledMatchedPointsNLP(const deque< Vector > &_p0, const deque< Vector > &_p1, const Vector &_min, const Vector &_max)
bool get_nlp_info(Ipopt::Index &n, Ipopt::Index &m, Ipopt::Index &nnz_jac_g, Ipopt::Index &nnz_h_lag, IndexStyleEnum &index_style)
bool eval_f(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Number &obj_value)
bool get_nlp_info(Ipopt::Index &n, Ipopt::Index &m, Ipopt::Index &nnz_jac_g, Ipopt::Index &nnz_h_lag, IndexStyleEnum &index_style)
CalibReferenceWithScaledMatchedPointsNLP(const deque< Vector > &_p0, const deque< Vector > &_p1, const Vector &_min, const Vector &_max)
bool eval_grad_f(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Number *grad_f)
bool eval_f(Ipopt::Index n, const Ipopt::Number *x, bool new_x, Ipopt::Number &obj_value)
double dot(const yarp::sig::Matrix &A, int colA, const yarp::sig::Matrix &B, int colB)
Returns the dot product between two vectors given in the form: matrix(:,col).
double norm2(const yarp::sig::Matrix &M, int col)
Returns the squared norm of the vector given in the form: matrix(:,col).
double norm(const yarp::sig::Matrix &M, int col)
Returns the norm of the vector given in the form: matrix(:,col).
Matrix computeH(const Ipopt::Number *x)
This file contains the definition of unique IDs for the body parts and the skin parts of the robot.