diff --git a/include/robust-equilibrium-lib/static_equilibrium.hh b/include/robust-equilibrium-lib/static_equilibrium.hh index cfb26a4e34069bb950aa42e08245e01af62495c7..1eed5f438b9ac38c01186ccdb22b57cb1125e5e2 100644 --- a/include/robust-equilibrium-lib/static_equilibrium.hh +++ b/include/robust-equilibrium-lib/static_equilibrium.hh @@ -53,8 +53,23 @@ private: Matrix62 m_D; Vector6 m_d; + /** Coefficient used for converting the robustness measure in Newtons */ + double m_b0_to_emax_coefficient; + bool computePolytopeProjection(Cref_matrix6X v); + /** + * @brief Given the minimum coefficient of the contact force generators it computes + * the minimum norm of force error necessary to result in a contact force being on + * the friction cone boundaries. + * @param b0 Minimum coefficient of the contact force generators. + * @return Minimum norm of the force error necessary to result in a contact force being + * on the friction cone boundaries. + */ + double convert_b0_to_emax(double b0); + + double convert_emax_to_b0(double emax); + public: /** @@ -78,12 +93,12 @@ public: * In other words the gravity vecotr is (0, 0, -9.81). * @param contactPoints List of N 3d contact points as an Nx3 matrix. * @param contactNormals List of N 3d contact normal directions as an Nx3 matrix. - * @param frictionCoefficients List of N friction coefficients. + * @param frictionCoefficient The contact friction coefficient. * @param alg Algorithm to use for testing equilibrium. * @return True if the operation succeeded, false otherwise. */ bool setNewContacts(Cref_matrixX3 contactPoints, Cref_matrixX3 contactNormals, - Cref_vectorX frictionCoefficients, StaticEquilibriumAlgorithm alg); + double frictionCoefficient, StaticEquilibriumAlgorithm alg); /** * @brief Compute a measure of the robustness of the equilibrium of the specified com position. diff --git a/src/static_equilibrium.cpp b/src/static_equilibrium.cpp index a5ec8244a7fca2e4520a7158a0a89160f4426ba0..f7e4705f41defa1d0bb98d1f3cc69f4a0c88c799 100644 --- a/src/static_equilibrium.cpp +++ b/src/static_equilibrium.cpp @@ -47,10 +47,9 @@ StaticEquilibrium::StaticEquilibrium(string name, double mass, unsigned int gene } bool StaticEquilibrium::setNewContacts(Cref_matrixX3 contactPoints, Cref_matrixX3 contactNormals, - Cref_vectorX frictionCoefficients, StaticEquilibriumAlgorithm alg) + double frictionCoefficient, StaticEquilibriumAlgorithm alg) { assert(contactPoints.rows()==contactNormals.rows()); - assert(contactPoints.rows()==frictionCoefficients.rows()); if(alg==STATIC_EQUILIBRIUM_ALGORITHM_IP) { @@ -65,7 +64,6 @@ bool StaticEquilibrium::setNewContacts(Cref_matrixX3 contactPoints, Cref_matrixX m_algorithm = alg; - // compute tangent directions long int c = contactPoints.rows(); unsigned int &cg = m_generatorsPerContact; double theta, delta_theta=2*M_PI/cg; @@ -101,8 +99,8 @@ bool StaticEquilibrium::setNewContacts(Cref_matrixX3 contactPoints, Cref_matrixX theta = 0.0; for(int j=0; j<cg; j++) { - G.col(j) = frictionCoefficients(i)*sin(theta)*T1 - + frictionCoefficients(i)*cos(theta)*T2 + G.col(j) = frictionCoefficient*sin(theta)*T1 + + frictionCoefficient*cos(theta)*T2 + contactNormals.row(i).transpose(); G.col(j).normalize(); // SEND_DEBUG_MSG("Contact "+toString(i)+" generator "+toString(j)+" = "+toString(G.col(j).transpose())); @@ -113,6 +111,16 @@ bool StaticEquilibrium::setNewContacts(Cref_matrixX3 contactPoints, Cref_matrixX m_G_centr.block(0,cg*i,6,cg) = A * G; } + // Compute the coefficient to convert b0 to e_max + Vector3 f0 = Vector3::Zero(); + for(int j=0; j<cg; j++) + f0 += G.col(j); // sum of the contact generators + // Compute the distance between the friction cone boundaries and + // the sum of the contact generators, which is e_max when b0=1. + // When b0!=1 we just multiply b0 times this value. + // This value depends only on the number of generators and the friction coefficient + m_b0_to_emax_coefficient = (f0.cross(G.col(0))).norm(); + if(m_algorithm==STATIC_EQUILIBRIUM_ALGORITHM_PP) { if(!computePolytopeProjection(m_G_centr)) @@ -160,7 +168,7 @@ LP_status StaticEquilibrium::computeEquilibriumRobustness(Cref_vector2 com, doub LP_status lpStatus = m_solver->solve(c, lb, ub, A, Alb, Aub, b_b0); if(lpStatus==LP_STATUS_OPTIMAL) { - robustness = -1.0*m_solver->getObjectiveValue(); + robustness = convert_b0_to_emax(-1.0*m_solver->getObjectiveValue()); return lpStatus; } @@ -197,7 +205,7 @@ LP_status StaticEquilibrium::computeEquilibriumRobustness(Cref_vector2 com, doub LP_status lpStatus_primal = m_solver->solve(c, lb, ub, A, Alb, Aub, b_b0); if(lpStatus_primal==LP_STATUS_OPTIMAL) { - robustness = -1.0*m_solver->getObjectiveValue(); + robustness = convert_b0_to_emax(-1.0*m_solver->getObjectiveValue()); return lpStatus_primal; } @@ -233,7 +241,7 @@ LP_status StaticEquilibrium::computeEquilibriumRobustness(Cref_vector2 com, doub LP_status lpStatus_dual = m_solver->solve(c, lb, ub, A, Alb, Aub, v); if(lpStatus_dual==LP_STATUS_OPTIMAL) { - robustness = m_solver->getObjectiveValue(); + robustness = convert_b0_to_emax(m_solver->getObjectiveValue()); return lpStatus_dual; } @@ -273,6 +281,7 @@ LP_status StaticEquilibrium::checkRobustEquilibrium(Cref_vector2 com, bool &equi LP_status StaticEquilibrium::findExtremumOverLine(Cref_vector2 a, Cref_vector2 a0, double e_max, Ref_vector2 com) { const long m = m_G_centr.cols(); // number of gravito-inertial wrench generators + double b0 = convert_emax_to_b0(e_max); if(m_algorithm==STATIC_EQUILIBRIUM_ALGORITHM_LP) { @@ -294,7 +303,7 @@ LP_status StaticEquilibrium::findExtremumOverLine(Cref_vector2 a, Cref_vector2 a VectorX lb = -VectorX::Zero(m+1); lb(m) = -1e5; VectorX ub = VectorX::Ones(m+1)*1e10; - Vector6 Alb = m_D*a0 + m_d - m_G_centr*VectorX::Ones(m)*e_max; + Vector6 Alb = m_D*a0 + m_d - m_G_centr*VectorX::Ones(m)*b0; Vector6 Aub = Alb; Matrix6X A = Matrix6X::Zero(6, m+1); A.leftCols(m) = m_G_centr; @@ -338,7 +347,7 @@ LP_status StaticEquilibrium::findExtremumOverLine(Cref_vector2 a, Cref_vector2 a G is the matrix whose columns are the gravito-inertial wrench generators */ Vector6 v; - Vector6 c = m_D*a0 + m_d - m_G_centr*VectorX::Ones(m)*e_max; + Vector6 c = m_D*a0 + m_d - m_G_centr*VectorX::Ones(m)*b0; Vector6 lb = Vector6::Ones()*-1e10; Vector6 ub = Vector6::Ones()*1e10; VectorX Alb = VectorX::Zero(m+1); @@ -441,4 +450,14 @@ bool StaticEquilibrium::computePolytopeProjection(Cref_matrix6X v) return true; } +double StaticEquilibrium::convert_b0_to_emax(double b0) +{ + return (b0*m_b0_to_emax_coefficient); +} + +double StaticEquilibrium::convert_emax_to_b0(double emax) +{ + return (emax/m_b0_to_emax_coefficient); +} + } // end namespace robust_equilibrium diff --git a/test/test_static_equilibrium.cpp b/test/test_static_equilibrium.cpp index feb9cc8d72e83ee1389f23c82313f8d291e83876..101846d1880a7636aec8dc29dbe8435605ce7aad 100644 --- a/test/test_static_equilibrium.cpp +++ b/test/test_static_equilibrium.cpp @@ -303,8 +303,6 @@ int main() MatrixXX contact_rpy = MatrixXX::Zero(N_CONTACTS, 3); MatrixXX p = MatrixXX::Zero(4*N_CONTACTS,3); // contact points MatrixXX N = MatrixXX::Zero(4*N_CONTACTS,3); // contact normals - VectorX frictionCoefficients(4*N_CONTACTS); - frictionCoefficients.fill(mu); // Generate contact positions and orientations bool collision; @@ -398,7 +396,7 @@ int main() } getProfiler().start(PERF_LP_PREPARATION); - if(!solver_LP_oases.setNewContacts(p, N, frictionCoefficients, STATIC_EQUILIBRIUM_ALGORITHM_LP)) + if(!solver_LP_oases.setNewContacts(p, N, mu, STATIC_EQUILIBRIUM_ALGORITHM_LP)) { printf("Error while setting new contacts"); return -1; @@ -406,7 +404,7 @@ int main() getProfiler().stop(PERF_LP_PREPARATION); getProfiler().start(PERF_LP_PREPARATION); - if(!solver_LP2_oases.setNewContacts(p, N, frictionCoefficients, STATIC_EQUILIBRIUM_ALGORITHM_LP2)) + if(!solver_LP2_oases.setNewContacts(p, N, mu, STATIC_EQUILIBRIUM_ALGORITHM_LP2)) { printf("Error while setting new contacts"); return -1; @@ -414,7 +412,7 @@ int main() getProfiler().stop(PERF_LP_PREPARATION); getProfiler().start(PERF_LP_PREPARATION); - if(!solver_DLP_oases.setNewContacts(p, N, frictionCoefficients, STATIC_EQUILIBRIUM_ALGORITHM_DLP)) + if(!solver_DLP_oases.setNewContacts(p, N, mu, STATIC_EQUILIBRIUM_ALGORITHM_DLP)) { printf("Error while setting new contacts"); return -1; @@ -422,7 +420,7 @@ int main() getProfiler().stop(PERF_LP_PREPARATION); getProfiler().start(PERF_PP); - bool res = solver_PP.setNewContacts(p, N, frictionCoefficients, STATIC_EQUILIBRIUM_ALGORITHM_PP); + bool res = solver_PP.setNewContacts(p, N, mu, STATIC_EQUILIBRIUM_ALGORITHM_PP); getProfiler().stop(PERF_PP); if(!res) { @@ -432,7 +430,7 @@ int main() #ifdef CLP_FOUND getProfiler().start(PERF_LP_PREPARATION); - if(!solver_LP_coin.setNewContacts(p, N, frictionCoefficients, STATIC_EQUILIBRIUM_ALGORITHM_LP)) + if(!solver_LP_coin.setNewContacts(p, N, mu, STATIC_EQUILIBRIUM_ALGORITHM_LP)) { printf("Error while setting new contacts"); return -1; @@ -440,7 +438,7 @@ int main() getProfiler().stop(PERF_LP_PREPARATION); getProfiler().start(PERF_LP_PREPARATION); - if(!solver_LP2_coin.setNewContacts(p, N, frictionCoefficients, STATIC_EQUILIBRIUM_ALGORITHM_LP2)) + if(!solver_LP2_coin.setNewContacts(p, N, mu, STATIC_EQUILIBRIUM_ALGORITHM_LP2)) { printf("Error while setting new contacts"); return -1; @@ -448,7 +446,7 @@ int main() getProfiler().stop(PERF_LP_PREPARATION); getProfiler().start(PERF_LP_PREPARATION); - if(!solver_DLP_coin.setNewContacts(p, N, frictionCoefficients, STATIC_EQUILIBRIUM_ALGORITHM_DLP)) + if(!solver_DLP_coin.setNewContacts(p, N, mu, STATIC_EQUILIBRIUM_ALGORITHM_DLP)) { printf("Error while setting new contacts"); return -1;