WDLSSolver.cpp 2.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105
  1. /** \file itasc/WDLSSolver.cpp
  2. * \ingroup itasc
  3. */
  4. /*
  5. * WDLSSolver.hpp.cpp
  6. *
  7. * Created on: Jan 8, 2009
  8. * Author: rubensmits
  9. */
  10. #include "WDLSSolver.hpp"
  11. #include "kdl/utilities/svd_eigen_HH.hpp"
  12. namespace iTaSC {
  13. WDLSSolver::WDLSSolver() : m_lambda(0.5), m_epsilon(0.1)
  14. {
  15. // maximum joint velocity
  16. m_qmax = 50.0;
  17. }
  18. WDLSSolver::~WDLSSolver() {
  19. }
  20. bool WDLSSolver::init(unsigned int nq, unsigned int nc, const std::vector<bool>& gc)
  21. {
  22. m_ns = std::min(nc,nq);
  23. m_AWq = e_zero_matrix(nc,nq);
  24. m_WyAWq = e_zero_matrix(nc,nq);
  25. m_WyAWqt = e_zero_matrix(nq,nc);
  26. m_S = e_zero_vector(std::max(nc,nq));
  27. m_Wy_ydot = e_zero_vector(nc);
  28. if (nq > nc) {
  29. m_transpose = true;
  30. m_temp = e_zero_vector(nc);
  31. m_U = e_zero_matrix(nc,nc);
  32. m_V = e_zero_matrix(nq,nc);
  33. m_WqV = e_zero_matrix(nq,nc);
  34. } else {
  35. m_transpose = false;
  36. m_temp = e_zero_vector(nq);
  37. m_U = e_zero_matrix(nc,nq);
  38. m_V = e_zero_matrix(nq,nq);
  39. m_WqV = e_zero_matrix(nq,nq);
  40. }
  41. return true;
  42. }
  43. bool WDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef)
  44. {
  45. double alpha, vmax, norm;
  46. // Create the Weighted jacobian
  47. m_AWq = A*Wq;
  48. for (int i=0; i<Wy.size(); i++)
  49. m_WyAWq.row(i) = Wy(i)*m_AWq.row(i);
  50. // Compute the SVD of the weighted jacobian
  51. int ret;
  52. if (m_transpose) {
  53. m_WyAWqt = m_WyAWq.transpose();
  54. ret = KDL::svd_eigen_HH(m_WyAWqt,m_V,m_S,m_U,m_temp);
  55. } else {
  56. ret = KDL::svd_eigen_HH(m_WyAWq,m_U,m_S,m_V,m_temp);
  57. }
  58. if(ret<0)
  59. return false;
  60. m_WqV.noalias() = Wq*m_V;
  61. //Wy*ydot
  62. m_Wy_ydot = Wy.array() * ydot.array();
  63. //S^-1*U'*Wy*ydot
  64. e_scalar maxDeltaS = e_scalar(0.0);
  65. e_scalar prevS = e_scalar(0.0);
  66. e_scalar maxS = e_scalar(1.0);
  67. e_scalar S, lambda;
  68. qdot.setZero();
  69. for(int i=0;i<m_ns;++i) {
  70. S = m_S(i);
  71. if (S <= KDL::epsilon)
  72. break;
  73. if (i > 0 && (prevS-S) > maxDeltaS) {
  74. maxDeltaS = (prevS-S);
  75. maxS = prevS;
  76. }
  77. lambda = (S < m_epsilon) ? (e_scalar(1.0)-KDL::sqr(S/m_epsilon))*m_lambda*m_lambda : e_scalar(0.0);
  78. alpha = m_U.col(i).dot(m_Wy_ydot)*S/(S*S+lambda);
  79. vmax = m_WqV.col(i).array().abs().maxCoeff();
  80. norm = fabs(alpha*vmax);
  81. if (norm > m_qmax) {
  82. qdot += m_WqV.col(i)*(alpha*m_qmax/norm);
  83. } else {
  84. qdot += m_WqV.col(i)*alpha;
  85. }
  86. prevS = S;
  87. }
  88. if (maxDeltaS == e_scalar(0.0))
  89. nlcoef = e_scalar(KDL::epsilon);
  90. else
  91. nlcoef = (maxS-maxDeltaS)/maxS;
  92. return true;
  93. }
  94. }