Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add benchmark example for dircon #107

Merged
merged 34 commits into from
Nov 6, 2019
Merged

Add benchmark example for dircon #107

merged 34 commits into from
Nov 6, 2019

Conversation

yminchen
Copy link
Contributor

@yminchen yminchen commented Oct 21, 2019

Fix bugs in dircon and add a benchmark example for dircon's speed.

Detail changes:
Dircon:

  • fix a bug. Use the jacobian that doesn't skip any constraints in the equations of motion
  • fix bugs in the index of decision variable when adding constraint
  • pass in surface normal into constructor of dircon_position_data
  • add surface (e.g. ground) incline feature to dircon_position_data
  • use linear constraint instead of LorentzConeConstraint for friction cone (to increase speed)
  • add a utility function to detect if a MBP is quaternion floating based
  • add unit norm quaternion constraint and also add a slack variable in dynamic constraint

Example:

  • Trajectory optimization of Cassie squating

@yminchen yminchen requested a review from mposa October 21, 2019 19:57
@@ -109,7 +110,7 @@ void DirconKinematicDataSet<T>::updateData(const Context<T>& context,
right_hand_side_ = -right_hand_side_ +
plant_.MakeActuationMatrix() * input +
plant_.CalcGravityGeneralizedForces(context) +
getJ().transpose() * forces;
J_.transpose() * forces;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be getJWithoutSkipping? That might be the best way to avoid a similar error down the road.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, getJWithoutSkipping() is probably better here

@@ -26,6 +26,8 @@ using drake::VectorX;
using drake::MatrixX;
using Eigen::VectorXd;
using Eigen::MatrixXd;
using std::cout;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably better to not bloat files with using's for debugging code that shouldn't be merged to master anyway.

@@ -69,15 +70,17 @@ class DirconDynamicConstraint : public DirconAbstractConstraint<T> {
DirconDynamicConstraint(const drake::multibody::MultibodyPlant<T>& plant,
DirconKinematicDataSet<T>& constraints,
int num_positions, int num_velocities, int num_inputs,
int num_kinematic_constraints);
int num_kinematic_constraints_wo_skipping,
int num_quat_slack);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is num_quat_slack for? I know this file needs documentation badly, and that's on me, but it would be good to document any new additions.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

num_quat_slack is the dimension of the slack variable for quaternion. It's always 1 if the MBP is quaternion floating based. Otherwise, it's 0. I'll add documentation about this.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't this be documented in the header (here)?

@@ -13,7 +13,8 @@ class DirconPositionData : public DirconKinematicData<T> {
public:
DirconPositionData(const drake::multibody::MultibodyPlant<T>& plant,
const drake::multibody::Body<T>& body,
Eigen::Vector3d pt, bool isXZ = false);
Eigen::Vector3d pt, bool isXZ = false,
Eigen::Vector2d surface_slope_roll_pitch = Eigen::Vector2d::Zero());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't surface normal, as a Vector3d be easier? Roll-pitch is a little ambiguous.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree. I replaced roll and pitch with normal.

/// With the default initial guess, the solving time is about 100 seconds.

// Constraint to fix the position of a point on a body (for initial guess)
class BodyPonitPositionConstraint : public DirconAbstractConstraint<double> {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo in class name BodyPonitPositionConstraint

ContactInfo contact_info(xa, idxa);

VectorXd q_desired = VectorXd::Zero(n_q);
q_desired << -0.0872062,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where did these magic numbers come from (comment would be useful). Can these be rounded?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It came from inverse kinematics where cassie was standing on the ground.
I also added documentation and rounded some of the numbers.

examples/Cassie/run_dircon_squatting.cc Show resolved Hide resolved

// For N-1 timesteps, add a constraint which depends on the knot
// value along with the state and input vectors at that knot and the
// next.

//Adding dynamic constraints
auto constraint = std::make_shared<DirconDynamicConstraint<T>>(plant_, *constraints_[i], num_quat_slack);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't the quaternion constraint be created and added here, rather than in the main file?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. I made this change in the new commits

@yminchen
Copy link
Contributor Author

yminchen commented Oct 26, 2019

@mposa I updated the code and it's ready for another review

Copy link
Contributor

@mposa mposa left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some comments. It looks like you auto-formatted some of the white space, and a lot of it ended up looking a little off. Are you sure it's correct?

template <typename T>
QuaternionNormConstraint<T>::QuaternionNormConstraint() :
DirconAbstractConstraint<T>(1, 4,
VectorXd::Zero(1),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Spacing here looks off?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I forgot to update the code spacing of this file

@@ -69,15 +70,17 @@ class DirconDynamicConstraint : public DirconAbstractConstraint<T> {
DirconDynamicConstraint(const drake::multibody::MultibodyPlant<T>& plant,
DirconKinematicDataSet<T>& constraints,
int num_positions, int num_velocities, int num_inputs,
int num_kinematic_constraints);
int num_kinematic_constraints_wo_skipping,
int num_quat_slack);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't this be documented in the header (here)?

isXZ_(isXZ) {
const Body<T>& body, Vector3d pt, bool isXZ,
Vector3d surface_normal) :
DirconKinematicData<T>(plant, isXZ ? 2 : 3),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Spacing here looks strange?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I'll update the code style of this file

TXZ_ << 1, 0, 0,
0, 0, 1;

Vector3d z_hat(0,0,1);
surface_normal.normalize();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does it make sense to supply a surface_normal value if isXZ = true? If these two arguments cannot be supplied together, then this should either be prevented via the API or, at minimum, checked at runtime.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I use both of them in the example of planar five-link robot example.

normal_xz << normal(0)/L, normal(2)/L;
d_xz << -normal_xz(1), normal_xz(0);
// builds the basis
Vector2d normal_xz(0,1);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Related to earlier comment re: normal and isXZ. This looks like a change where, previously, we supported non-flat normal vectors for planar contact. Is this no longer the case? If not, why are we removing this functionality?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Currently there are c, cdot and cddot in dircon_kinemaitcs_data. We want c=0 in the normal direction of the surface, so we rotate the contact position into the surface frame. (I guess we can generalize this more by adding a translation, if there is a specific position where we want the robot to touch the surface? Currently the surface has to contain the global origin.)

Since we rotate c, cdot and cddot because of the above reason, we don't need to rotate it again in the addFixedNormalFrictionConstraints().

vector<double> maximum_timestep,
vector<DirconKinematicDataSet<T>*> constraints,
vector<DirconOptions> options)
: MultipleShooting(plant.num_actuators(),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Spacing looks off here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I couldn't tell how it is off. I used a code style formatter to format hybrid_dircon.cc.
The style should be the same (if I didn't miss any thing) as Drake's. https://github.com/RobotLocomotion/drake/blob/master/systems/trajectory_optimization/multiple_shooting.cc#L63-L92

num_modes_(num_time_samples.size()),
mode_lengths_(num_time_samples),
v_post_impact_vars_(NewContinuousVariables(
plant.num_velocities() *
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment. I couldn't tell how it is off. I think this is the same as drake's style?

constraints_[i]->countConstraintsWithoutSkipping()
* (num_time_samples[i] - 1), "v_c[" + std::to_string(i) + "]"));
// slack variables used to scale quaternion norm to 1 in the dynamic constraints.
quaternion_slack_vars_.push_back(NewContinuousVariables(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't very readable. Why not

if (is_quaternion) {
// Add variables

}

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel they are equally readable, but I will change to

if (is_quaternion) {
// Add variables

}


// For N-1 timesteps, add a constraint which depends on the knot
// value along with the state and input vectors at that knot and the
// next.


//Adding quaternion norm constraint
if (multibody::isQuaternion(plant)) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should use the boolean you've created, is_quaternion here.

if (multibody::isQuaternion(plant)) {
auto quat_norm_constraint =
std::make_shared<QuaternionNormConstraint<T>>();
for (int j = (i == 0) ? 0 : 1; j < mode_lengths_[i]; j++) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment here to explain the use of x ? y : z.

@yminchen
Copy link
Contributor Author

yminchen commented Nov 4, 2019

@mposa I've added some documentation and updated code spacing, but I'm not sure if this is the code style you want.

@mposa
Copy link
Contributor

mposa commented Nov 4, 2019

Per face to face with Yu-Ming, add an option to fix the seed for the RNG so runtime comparisons are consistent.

Copy link
Contributor

@mposa mposa left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lgtm. Thank you!

@yminchen yminchen merged commit ea86f21 into master Nov 6, 2019
yminchen added a commit that referenced this pull request Aug 4, 2020
**Detail changes**:
Dircon:
- fix a bug. Use the jacobian that doesn't skip any constraints in the equations of motion
- fix bugs in the index of decision variable when adding constraint 
- pass in surface normal into constructor of dircon_position_data
- add surface (e.g. ground) incline feature to dircon_position_data
- use linear constraint instead of LorentzConeConstraint for friction cone (to increase speed)
- add a utility function to detect if a MBP is quaternion floating based
- add unit norm quaternion constraint and also add a slack variable in dynamic constraint 

Example:
- Trajectory optimization of Cassie squating
MaDaO009 pushed a commit to KodlabPenn/dairlib that referenced this pull request May 23, 2023
**Detail changes**:
Dircon:
- fix a bug. Use the jacobian that doesn't skip any constraints in the equations of motion
- fix bugs in the index of decision variable when adding constraint 
- pass in surface normal into constructor of dircon_position_data
- add surface (e.g. ground) incline feature to dircon_position_data
- use linear constraint instead of LorentzConeConstraint for friction cone (to increase speed)
- add a utility function to detect if a MBP is quaternion floating based
- add unit norm quaternion constraint and also add a slack variable in dynamic constraint 

Example:
- Trajectory optimization of Cassie squating
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants