Blender  V3.3
WDLSSolver.cpp
Go to the documentation of this file.
1 /* SPDX-License-Identifier: LGPL-2.1-or-later
2  * Copyright 2009 Ruben Smits. */
3 
8 #include "WDLSSolver.hpp"
10 
11 namespace iTaSC {
12 
13 WDLSSolver::WDLSSolver() : m_lambda(0.5), m_epsilon(0.1)
14 {
15  // maximum joint velocity
16  m_qmax = 50.0;
17 }
18 
20 }
21 
22 bool WDLSSolver::init(unsigned int nq, unsigned int nc, const std::vector<bool>& gc)
23 {
24  m_ns = std::min(nc,nq);
25  m_AWq = e_zero_matrix(nc,nq);
26  m_WyAWq = e_zero_matrix(nc,nq);
27  m_WyAWqt = e_zero_matrix(nq,nc);
28  m_S = e_zero_vector(std::max(nc,nq));
29  m_Wy_ydot = e_zero_vector(nc);
30  if (nq > nc) {
31  m_transpose = true;
32  m_temp = e_zero_vector(nc);
33  m_U = e_zero_matrix(nc,nc);
34  m_V = e_zero_matrix(nq,nc);
35  m_WqV = e_zero_matrix(nq,nc);
36  } else {
37  m_transpose = false;
38  m_temp = e_zero_vector(nq);
39  m_U = e_zero_matrix(nc,nq);
40  m_V = e_zero_matrix(nq,nq);
41  m_WqV = e_zero_matrix(nq,nq);
42  }
43  return true;
44 }
45 
46 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)
47 {
48  double alpha, vmax, norm;
49  // Create the Weighted jacobian
50  m_AWq = A*Wq;
51  for (int i=0; i<Wy.size(); i++)
52  m_WyAWq.row(i) = Wy(i)*m_AWq.row(i);
53 
54  // Compute the SVD of the weighted jacobian
55  int ret;
56  if (m_transpose) {
57  m_WyAWqt = m_WyAWq.transpose();
58  ret = KDL::svd_eigen_HH(m_WyAWqt,m_V,m_S,m_U,m_temp);
59  } else {
60  ret = KDL::svd_eigen_HH(m_WyAWq,m_U,m_S,m_V,m_temp);
61  }
62  if(ret<0)
63  return false;
64 
65  m_WqV.noalias() = Wq*m_V;
66 
67  //Wy*ydot
68  m_Wy_ydot = Wy.array() * ydot.array();
69  //S^-1*U'*Wy*ydot
70  e_scalar maxDeltaS = e_scalar(0.0);
71  e_scalar prevS = e_scalar(0.0);
72  e_scalar maxS = e_scalar(1.0);
73  e_scalar S, lambda;
74  qdot.setZero();
75  for(int i=0;i<m_ns;++i) {
76  S = m_S(i);
77  if (S <= KDL::epsilon)
78  break;
79  if (i > 0 && (prevS-S) > maxDeltaS) {
80  maxDeltaS = (prevS-S);
81  maxS = prevS;
82  }
83  lambda = (S < m_epsilon) ? (e_scalar(1.0)-KDL::sqr(S/m_epsilon))*m_lambda*m_lambda : e_scalar(0.0);
84  alpha = m_U.col(i).dot(m_Wy_ydot)*S/(S*S+lambda);
85  vmax = m_WqV.col(i).array().abs().maxCoeff();
86  norm = fabs(alpha*vmax);
87  if (norm > m_qmax) {
88  qdot += m_WqV.col(i)*(alpha*m_qmax/norm);
89  } else {
90  qdot += m_WqV.col(i)*alpha;
91  }
92  prevS = S;
93  }
94  if (maxDeltaS == e_scalar(0.0))
95  nlcoef = e_scalar(KDL::epsilon);
96  else
97  nlcoef = (maxS-maxDeltaS)/maxS;
98  return true;
99 }
100 
101 }
#define A
SIMD_FORCE_INLINE btScalar norm() const
Return the norm (length) of the vector.
Definition: btVector3.h:263
virtual bool init(unsigned int nq, unsigned int nc, const std::vector< bool > &gc)
Definition: WDLSSolver.cpp:22
virtual ~WDLSSolver()
Definition: WDLSSolver.cpp:19
virtual bool solve(const e_matrix &A, const e_vector &Wy, const e_vector &ydot, const e_matrix &Wq, e_vector &qdot, e_scalar &nlcoef)
Definition: WDLSSolver.cpp:46
#define e_vector
Definition: eigen_types.hpp:38
#define e_scalar
Definition: eigen_types.hpp:37
#define e_zero_vector
Definition: eigen_types.hpp:39
#define e_zero_matrix
Definition: eigen_types.hpp:44
#define e_matrix
Definition: eigen_types.hpp:40
ccl_device_inline float2 fabs(const float2 &a)
Definition: math_float2.h:222
double epsilon
default precision while comparing with Equal(..,..) functions. Initialized at 0.0000001.
Definition: utility.cpp:22
int svd_eigen_HH(const Eigen::MatrixBase< MatrixA > &A, Eigen::MatrixBase< MatrixUV > &U, Eigen::MatrixBase< VectorS > &S, Eigen::MatrixBase< MatrixUV > &V, Eigen::MatrixBase< VectorS > &tmp, int maxiter=150)
INLINE Rall1d< T, V, S > sqr(const Rall1d< T, V, S > &arg)
Definition: rall1d.h:351
return ret
#define min(a, b)
Definition: sort.c:35
float max