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