diff options
Diffstat (limited to 'intern/itasc/WSDLSSolver.cpp')
| -rw-r--r-- | intern/itasc/WSDLSSolver.cpp | 138 |
1 files changed, 138 insertions, 0 deletions
diff --git a/intern/itasc/WSDLSSolver.cpp b/intern/itasc/WSDLSSolver.cpp new file mode 100644 index 00000000000..9f7ebed960a --- /dev/null +++ b/intern/itasc/WSDLSSolver.cpp @@ -0,0 +1,138 @@ +/* $Id$ + * WDLSSolver.hpp.cpp + * + * Created on: Jan 8, 2009 + * Author: rubensmits + */ + +#include "WSDLSSolver.hpp" +#include "kdl/utilities/svd_eigen_HH.hpp" +#include <cstdio> + +namespace iTaSC { + +WSDLSSolver::WSDLSSolver() : + m_ns(0), m_nc(0), m_nq(0) + +{ + // default maximum speed: 50 rad/s + m_qmax = 50.0; +} + +WSDLSSolver::~WSDLSSolver() { +} + +bool WSDLSSolver::init(unsigned int _nq, unsigned int _nc, const std::vector<bool>& gc) +{ + if (_nc == 0 || _nq == 0 || gc.size() != _nc) + return false; + m_nc = _nc; + m_nq = _nq; + m_ns = std::min(m_nc,m_nq); + m_AWq = e_zero_matrix(m_nc,m_nq); + m_WyAWq = e_zero_matrix(m_nc,m_nq); + m_WyAWqt = e_zero_matrix(m_nq,m_nc); + m_S = e_zero_vector(std::max(m_nc,m_nq)); + m_Wy_ydot = e_zero_vector(m_nc); + m_ytask = gc; + if (m_nq > m_nc) { + m_transpose = true; + m_temp = e_zero_vector(m_nc); + m_U = e_zero_matrix(m_nc,m_nc); + m_V = e_zero_matrix(m_nq,m_nc); + m_WqV = e_zero_matrix(m_nq,m_nc); + } else { + m_transpose = false; + m_temp = e_zero_vector(m_nq); + m_U = e_zero_matrix(m_nc,m_nq); + m_V = e_zero_matrix(m_nq,m_nq); + m_WqV = e_zero_matrix(m_nq,m_nq); + } + return true; +} + +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) +{ + unsigned int i, j, l; + e_scalar N, M; + + // Create the Weighted jacobian + m_AWq = (A*Wq).lazy(); + for (i=0; i<m_nc; 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_Wy_ydot = Wy.cwise() * ydot; + m_WqV = (Wq*m_V).lazy(); + qdot.setZero(); + e_scalar maxDeltaS = e_scalar(0.0); + e_scalar prevS = e_scalar(0.0); + e_scalar maxS = e_scalar(1.0); + for(i=0;i<m_ns;++i) { + e_scalar norm, mag, alpha, _qmax, Sinv, vmax, damp; + e_scalar S = m_S(i); + bool prev; + if (S < KDL::epsilon) + break; + Sinv = e_scalar(1.)/S; + if (i > 0) { + if ((prevS-S) > maxDeltaS) { + maxDeltaS = (prevS-S); + maxS = prevS; + } + } + N = M = e_scalar(0.); + for (l=0, prev=m_ytask[0], norm=e_scalar(0.); l<m_nc; l++) { + if (prev == m_ytask[l]) { + norm += m_U(l,i)*m_U(l,i); + } else { + N += std::sqrt(norm); + norm = m_U(l,i)*m_U(l,i); + } + prev = m_ytask[l]; + } + N += std::sqrt(norm); + for (j=0; j<m_nq; j++) { + for (l=0, prev=m_ytask[0], norm=e_scalar(0.), mag=e_scalar(0.); l<m_nc; l++) { + if (prev == m_ytask[l]) { + norm += m_WyAWq(l,j)*m_WyAWq(l,j); + } else { + mag += std::sqrt(norm); + norm = m_WyAWq(l,j)*m_WyAWq(l,j); + } + prev = m_ytask[l]; + } + mag += std::sqrt(norm); + M += fabs(m_V(j,i))*mag; + } + M *= Sinv; + alpha = m_U.col(i).dot(m_Wy_ydot); + _qmax = (N < M) ? m_qmax*N/M : m_qmax; + vmax = m_WqV.col(i).cwise().abs().maxCoeff(); + norm = fabs(Sinv*alpha*vmax); + if (norm > _qmax) { + damp = Sinv*alpha*_qmax/norm; + } else { + damp = Sinv*alpha; + } + qdot += m_WqV.col(i)*damp; + prevS = S; + } + if (maxDeltaS == e_scalar(0.0)) + nlcoef = e_scalar(KDL::epsilon); + else + nlcoef = (maxS-maxDeltaS)/maxS; + return true; +} + +} |