00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <robot/control/kinematics/SVDFullSpaceSolver>
00026
00027 using robot::control::kinematics::SVDFullSpaceSolver;
00028
00029 #include <base/SVD>
00030 #include <robot/control/kinematics/solution_error>
00031
00032
00033
00034
00035
00036
00037 using base::IVector;
00038 using base::zeroIVector;
00039 using base::equals;
00040 using base::SVD;
00041
00042 const Real small = 1.0e-04;
00043
00044
00045 SVDFullSpaceSolver::SVDFullSpaceSolver()
00046 : stopOnIllCondition(false)
00047 {
00048 }
00049
00050
00051
00052 const static Real maxCondition = 10000;
00053
00054 base::Matrix SVDFullSpaceSolver::solve(const Matrix& A_in, const Vector& b_in,
00055 array<Int>& dependentRowsEliminated)
00056 {
00057 const Int N = A_in.rows();
00058 const Int M = A_in.cols();
00059
00060 const Vector& b(b_in);
00061 Assert( b.size() == N);
00062
00063 bool nullSpace = b.equals(zeroVector(N));
00064
00065
00066 SVD svd(A_in);
00067
00068 if (!nullSpace && (svd.condition() >= maxCondition)) {
00069 if (stopOnIllCondition)
00070 throw solution_error(Exception("FSP unable to complete solution: A is ill-conditioned (has large SVD condition number="
00071 +base::realToString(svd.condition())+")"));
00072 else
00073 Logln("Warning: A is ill-conditioned; singular values are:" << svd.diag() << " condition=" << svd.condition());
00074 }
00075
00076 const Matrix& U( svd.U() );
00077 const Matrix S( svd.S() );
00078 const Matrix& V( svd.V() );
00079 Matrix Uinv( transpose(U) );
00080
00081
00082
00083 Matrix Sinv( zeroMatrix(M,N) );
00084 for(Int i=0; i< Math::minimum(M, N); i++) {
00085 if (Math::abs(S(i,i)) > small)
00086 Sinv(i,i) = 1/S(i,i);
00087 else {
00088 if (stopOnIllCondition)
00089 throw solution_error(Exception(String("FSP unable to complete solution: a singular value of A is near 0 ")
00090 +"(S("+base::intToString(i)+","+base::intToString(i)+")="+base::realToString(S(i,i))+")"));
00091 else {
00092
00093 Sinv(i,i) = 0;
00094 }
00095 }
00096 }
00097
00098
00099 Matrix g;
00100
00101
00102 if (nullSpace) {
00103
00104
00105 g.resize(M, Math::maximum(SInt(M-N),0) );
00106 if (g.cols() > 0)
00107 for(Int i=0; i<M-N; i++)
00108 matrixColumn(g,i) = matrixColumn(V,N+i);
00109
00110 }
00111 else {
00112
00113 g.resize(M, Math::maximum(SInt(M-N+1),0) );
00114
00115 Vector Utb(Uinv * b);
00116
00117 if (g.cols() > 0) {
00118 matrixColumn(g,0) = V * Sinv * Utb;
00119
00120
00121 for(Int i=0; i<M-N; i++)
00122 matrixColumn(g,i+1)=matrixColumn(V,N+i) + matrixColumn(g,0);
00123 }
00124 }
00125
00126
00127 dependentRowsEliminated.clear();
00128
00129 return g;
00130 }
00131
00132
00133
00134
00135
00136
00137
00138
00139