Blender  V3.3
WSDLSSolver.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 "WSDLSSolver.hpp"
10 #include <cstdio>
11 
12 namespace iTaSC {
13 
15  m_ns(0), m_nc(0), m_nq(0)
16 
17 {
18  // default maximum speed: 50 rad/s
19  m_qmax = 50.0;
20 }
21 
23 }
24 
25 bool WSDLSSolver::init(unsigned int _nq, unsigned int _nc, const std::vector<bool>& gc)
26 {
27  if (_nc == 0 || _nq == 0 || gc.size() != _nc)
28  return false;
29  m_nc = _nc;
30  m_nq = _nq;
31  m_ns = std::min(m_nc,m_nq);
32  m_AWq = e_zero_matrix(m_nc,m_nq);
33  m_WyAWq = e_zero_matrix(m_nc,m_nq);
34  m_WyAWqt = e_zero_matrix(m_nq,m_nc);
35  m_S = e_zero_vector(std::max(m_nc,m_nq));
36  m_Wy_ydot = e_zero_vector(m_nc);
37  m_ytask = gc;
38  if (m_nq > m_nc) {
39  m_transpose = true;
40  m_temp = e_zero_vector(m_nc);
41  m_U = e_zero_matrix(m_nc,m_nc);
42  m_V = e_zero_matrix(m_nq,m_nc);
43  m_WqV = e_zero_matrix(m_nq,m_nc);
44  } else {
45  m_transpose = false;
46  m_temp = e_zero_vector(m_nq);
47  m_U = e_zero_matrix(m_nc,m_nq);
48  m_V = e_zero_matrix(m_nq,m_nq);
49  m_WqV = e_zero_matrix(m_nq,m_nq);
50  }
51  return true;
52 }
53 
54 bool WSDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef)
55 {
56  unsigned int i, j, l;
57  e_scalar N, M;
58 
59  // Create the Weighted jacobian
60  m_AWq.noalias() = A*Wq;
61  for (i=0; i<m_nc; i++)
62  m_WyAWq.row(i) = Wy(i)*m_AWq.row(i);
63 
64  // Compute the SVD of the weighted jacobian
65  int ret;
66  if (m_transpose) {
67  m_WyAWqt = m_WyAWq.transpose();
68  ret = KDL::svd_eigen_HH(m_WyAWqt,m_V,m_S,m_U,m_temp);
69  } else {
70  ret = KDL::svd_eigen_HH(m_WyAWq,m_U,m_S,m_V,m_temp);
71  }
72  if(ret<0)
73  return false;
74 
75  m_Wy_ydot = Wy.array() * ydot.array();
76  m_WqV.noalias() = Wq*m_V;
77  qdot.setZero();
78  e_scalar maxDeltaS = e_scalar(0.0);
79  e_scalar prevS = e_scalar(0.0);
80  e_scalar maxS = e_scalar(1.0);
81  for(i=0;i<m_ns;++i) {
82  e_scalar norm, mag, alpha, _qmax, Sinv, vmax, damp;
83  e_scalar S = m_S(i);
84  bool prev;
85  if (S < KDL::epsilon)
86  break;
87  Sinv = e_scalar(1.)/S;
88  if (i > 0) {
89  if ((prevS-S) > maxDeltaS) {
90  maxDeltaS = (prevS-S);
91  maxS = prevS;
92  }
93  }
94  N = M = e_scalar(0.);
95  for (l=0, prev=m_ytask[0], norm=e_scalar(0.); l<m_nc; l++) {
96  if (prev == m_ytask[l]) {
97  norm += m_U(l,i)*m_U(l,i);
98  } else {
99  N += std::sqrt(norm);
100  norm = m_U(l,i)*m_U(l,i);
101  }
102  prev = m_ytask[l];
103  }
104  N += std::sqrt(norm);
105  for (j=0; j<m_nq; j++) {
106  for (l=0, prev=m_ytask[0], norm=e_scalar(0.), mag=e_scalar(0.); l<m_nc; l++) {
107  if (prev == m_ytask[l]) {
108  norm += m_WyAWq(l,j)*m_WyAWq(l,j);
109  } else {
110  mag += std::sqrt(norm);
111  norm = m_WyAWq(l,j)*m_WyAWq(l,j);
112  }
113  prev = m_ytask[l];
114  }
115  mag += std::sqrt(norm);
116  M += fabs(m_V(j,i))*mag;
117  }
118  M *= Sinv;
119  alpha = m_U.col(i).dot(m_Wy_ydot);
120  _qmax = (N < M) ? m_qmax*N/M : m_qmax;
121  vmax = m_WqV.col(i).array().abs().maxCoeff();
122  norm = fabs(Sinv*alpha*vmax);
123  if (norm > _qmax) {
124  damp = Sinv*alpha*_qmax/norm;
125  } else {
126  damp = Sinv*alpha;
127  }
128  qdot += m_WqV.col(i)*damp;
129  prevS = S;
130  }
131  if (maxDeltaS == e_scalar(0.0))
132  nlcoef = e_scalar(KDL::epsilon);
133  else
134  nlcoef = (maxS-maxDeltaS)/maxS;
135  return true;
136 }
137 
138 }
sqrt(x)+1/max(0
ATTR_WARN_UNUSED_RESULT const BMLoop * l
#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: WSDLSSolver.cpp:25
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: WSDLSSolver.cpp:54
virtual ~WSDLSSolver()
Definition: WSDLSSolver.cpp:22
#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
#define M
#define N
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)
SymEdge< T > * prev(const SymEdge< T > *se)
Definition: delaunay_2d.cc:105
return ret
#define min(a, b)
Definition: sort.c:35
float max