Skip to content

Commit

Permalink
added LCP solver to thirdParty . working on joint limit example.
Browse files Browse the repository at this point in the history
git-svn-id: https://svn.csail.mit.edu/locomotion/robotlib/trunk@3552 c9849af7-e679-4ec6-a44e-fc146a885bd3
  • Loading branch information
russt committed Aug 17, 2012
1 parent 95f2b73 commit 4af540b
Show file tree
Hide file tree
Showing 13 changed files with 403 additions and 25 deletions.
27 changes: 16 additions & 11 deletions configure.m
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@
addpath([conf.root,'/systems/trajectories']);
addpath([conf.root,'/util']);
addpath([conf.root,'/util/obstacles']);
addpath([conf.root,'/thirdParty/path']);
addpath([conf.root,'/thirdParty/spatial']);
addpath([conf.root,'/thirdParty/graphviz2mat']);

Expand All @@ -74,17 +75,14 @@
% check for all dependencies

v=ver('simulink');
if (isempty(v)) conf.simulink_enabled = false;
else
ind = find(v.Version,'.','first');
major_ver = str2num(v.Version(1:ind));
minor_ver = str2num(v.Version((ind+2):end));
if (major_ver < 7 || (major_ver==7 && minor_ver < 3))
warning('Some features of Drake reguires SIMULINK version 7.3 or above.');
% haven't actually tested with lower versions
conf.simulink_enabled = false;
end
conf.simulink_enabled = true;
if (isempty(v))
conf.simulink_enabled = false;
elseif verLessThan('simulink','7.3')
warning('Some features of Drake reguires SIMULINK version 7.3 or above.');
% haven't actually tested with lower versions
conf.simulink_enabled = false;
else
conf.simulink_enabled = true;
end

conf.lcm_enabled = logical(exist('lcm.lcm.LCM'));
Expand Down Expand Up @@ -146,6 +144,13 @@
disp(' YALMIP not found. YALMIP support will be disabled. To re-enable, install YALMIP and rerun configure.');
end

conf.pathlcp_enabled = ~isempty(getenv('PATH_LICENSE_STRING'));
if (~conf.pathlcp_enabled)
disp('The PATH LCP solver (in the thirdparty directory) needs you to get the setup the license: http://pages.cs.wisc.edu/~ferris/path.html');
disp('I recommend adding a setenv(''PATH_LICENSE_STRING'',...) line to your startup.m');
disp('The LCP solver will be disabled');
end

conf.simulationconstructionset_enabled = ...
logical(exist('com.yobotics.simulationconstructionset.Robot'));
if (~conf.simulationconstructionset_enabled)
Expand Down
16 changes: 8 additions & 8 deletions systems/plants/Manipulator.m
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,16 @@
[H,C,B] = manipulatorDynamics(obj,q,qd);
Hinv = inv(H);

if (obj.num_u>0) tau=B*u; else tau=zeros(obj.num_q,1); end
tau = tau + computeConstraintForce(obj,q,qd,H,C,B,Hinv);
if (obj.num_u>0) tau=B*u - C; else tau=-C; end
tau = tau + computeConstraintForce(obj,q,qd,H,tau,Hinv);

qdd = Hinv*(tau - C);
qdd = Hinv*tau;
% note that I used to do this (instead of calling inv(H)):
% qdd = H\(tau - C)
% qdd = H\tau
% but I already have and use Hinv, so use it again here
end

function constraint_force = computeConstraintForce(obj,q,qd,H,C,B,Hinv)
function constraint_force = computeConstraintForce(obj,q,qd,H,tau,Hinv)
phi=[]; psi=[];
if (obj.num_position_constraints>0 && obj.num_velocity_constraints>0)
[phi,J,dJ] = geval(@obj.positionConstraints,q);
Expand All @@ -48,20 +48,20 @@
dpsidqd = dpsi(:,obj.num_q+1:end);

term1=Hinv*[J;dpsidqd]';
term2=Hinv*(tau-C);
term2=Hinv*tau;

constraint_force = -[J;dpsidqd]'*inv([J*term1;dpsidqd*term1])*[J*term2 + Jdotqd + alpha*J*qd; dpsidqd*term2 + dpsidq*qd + beta*psi];
elseif (obj.num_position_constraints>0) % note: it didn't work to just have dpsidq,etc=[], so it seems like the best solution is to handle each case...
[phi,J,dJ] = geval(@obj.positionConstraints,q);
Jdotqd = dJ*reshape(qd*qd',obj.num_q^2,1);

constraint_force = -J'*(inv(J*Hinv*J')*(J*Hinv*(tau-C) + Jdotqd + alpha*J*qd));
constraint_force = -J'*(inv(J*Hinv*J')*(J*Hinv*tau + Jdotqd + alpha*J*qd));
elseif (obj.num_velocity_constraints>0)
[psi,J] = geval(@obj.velocityConstraints,q,qd);
dpsidq = J(:,1:obj.num_q);
dpsidqd = J(:,obj.num_q+1:end);

constraint_force = -dpsidqd'*inv(dpsidqd*Hinv*dpsidqd')*(dpsidq*qd + dpsidqd*Hinv*(tau-C)+beta*psi);
constraint_force = -dpsidqd'*inv(dpsidqd*Hinv*dpsidqd')*(dpsidq*qd + dpsidqd*Hinv*tau+beta*psi);
else
constraint_force = 0*q;
end
Expand Down
16 changes: 10 additions & 6 deletions systems/plants/PlanarRigidBodyManipulatorWContact.m
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,9 @@

methods
function obj=PlanarRigidBodyManipulatorWContact(model)

checkDependency('pathlcp_enabled');

S = warning('off','Drake:PlanarRigidBodyManipulator:UnsupportedJointLimits');
obj = obj@PlanarRigidBodyManipulator(model);
warning(S);
Expand All @@ -33,18 +36,19 @@
end
end

function constraint_force = computeConstraintForce(obj,q,qd,H,C,B,Hinv)
function constraint_force = computeConstraintForce(obj,q,qd,H,tau,Hinv)
[phi,J,Jdot] = jointLimits(obj,q);
phidot = J*qd;
lambda=zeros(size(J,1),1);

epsinv = 1/.01;
phi
phidot
epsinv = 0; %1/.01;
lcp_ind = (phi<=0 & phidot <= -epsinv*phi);

if ~isempty(lcp_ind)
error('not implemented yet');
if any(lcp_ind)
Jlcp = J(lcp_ind,:);
M = Jlcp*Hinv*Jlcp';
b = Jdot(lcp_ind,:)*qd + Jlcp*Hinv*tau - 2*epsinv*phidot(lcp_ind) - epsinv^2*phi(lcp_ind);
lambda(lcp_ind) = pathlcp(M,b);
end
constraint_force = J'*lambda;
end
Expand Down
Binary file added thirdParty/path/lcppath.mexa64
Binary file not shown.
Binary file added thirdParty/path/lcppath.mexmaci64
Binary file not shown.
Binary file added thirdParty/path/lcppath.mexw32
Binary file not shown.
Binary file added thirdParty/path/lcppath.mexw64
Binary file not shown.
Binary file added thirdParty/path/mcppath.mexa64
Binary file not shown.
Binary file added thirdParty/path/mcppath.mexmaci64
Binary file not shown.
Binary file added thirdParty/path/mcppath.mexw32
Binary file not shown.
Binary file added thirdParty/path/mcppath.mexw64
Binary file not shown.
172 changes: 172 additions & 0 deletions thirdParty/path/pathlcp.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,172 @@
function [z,mu] = pathlcp(M,q,l,u,z,A,b,t,mu)
% pathlcp(M,q,l,u,z,A,b,t,mu)
%
% Solve the standard linear complementarity problem using PATH:
% z >= 0, Mz + q >= 0, z'*(Mz + q) = 0
%
% Required input:
% M(n,n) - matrix
% q(n) - vector
%
% Output:
% z(n) - solution
% mu(m) - multipliers (if polyhedral constraints are present)
%
% Optional input:
% l(n) - lower bounds default: zero
% u(n) - upper bounds default: infinity
% z(n) - starting point default: zero
% A(m,n) - polyhedral constraint matrix default: empty
% b(m) - polyhedral right-hand side default: empty
% t(m) - type of polyhedral constraint default: 1
% < 0: less than or equal
% 0: equation
% > 0: greater than or equal
% mu(m) - starting value for multipliers default: zero
%
% The optional lower and upper bounds are used to define a linear mixed
% complementarity problem (box constrained variational inequality).
% l <= z <= u
% where l_i < z_i < u_i => (Mz + q)_i = 0
% l_i = z => (Mz + q)_i >= 0
% u_i = z => (Mz + q)_i <= 0
%
% The optional constraints are used to define a polyhedrally constrained
% variational inequality. These are transformed internally to a standard
% mixed complementarity problem. The polyhedral constraints are of the
% form
% Ax ? b
% where ? can be <=, =, or >= depending on the type specified for each
% constraint.

Big = 1e20;

if (nargin < 2)
error('two input arguments required for lcp(M, q)');
end

if (~issparse(M))
M = sparse(M); % Make sure M is sparse
end
q = full(q(:)); % Make sure q is a column vector

[mm,mn] = size(M); % Get the size of the inputs
n = length(q);

if (mm ~= mn | mm ~= n)
error('dimensions of M and q must match');
end

if (n == 0)
error('empty model');
end

if (nargin < 3 | isempty(l))
l = zeros(n,1);
end

if (nargin < 4 | isempty(u))
u = Big*ones(n,1);
end

if (nargin < 5 | isempty(z))
z = zeros(n,1);
end

z = full(z(:)); l = full(l(:)); u = full(u(:));
if (length(z) ~= n | length(l) ~= n | length(u) ~= n)
error('Input arguments are of incompatible sizes');
end

l = max(l,-Big*ones(n,1));
u = min(u,Big*ones(n,1));
z = min(max(z,l),u);

m = 0;
if (nargin > 5)
if (nargin < 7)
error('Polyhedral constraints require A and b');
end

if (~issparse(A))
A = sparse(A);
end
b = full(b(:));

m = length(b);

if (m > 0)

[am, an] = size(A);

if (am ~= m | an ~= n)
error('Polyhedral constraints of incompatible sizes');
end

if (nargin < 8 | isempty(t))
t = ones(m,1);
end

if (nargin < 9 | isempty(mu))
mu = zeros(m,1);
end

t = full(t(:)); mu = full(mu(:));
if (length(t) ~= m | length(mu) ~= m)
error('Polyhedral input arguments are of incompatible sizes');
end

l_p = -Big*ones(m,1);
u_p = Big*ones(m,1);

idx = find(t > 0);
if (length(idx) > 0)
l_p(idx) = zeros(length(idx),1);
end

idx = find(t < 0);
if (length(idx) > 0)
u_p(idx) = zeros(length(idx),1);
end

mu = min(max(mu,l_p),u_p);

M = [M -A'; A sparse(m,m)];
q = [q; -b];

z = [z; mu];
l = [l; l_p];
u = [u; u_p];
else
if (nargin >= 9 & ~isempty(mu))
error('No polyhedral constraints -- multipliers set.');
end

if (nargin >= 8 & ~isempty(t))
error('No polyhedral constraints -- equation types set.');
end
end
end

idx = find(l > u);
if length(idx) > 0
error('Bounds infeasible.');
end

nnzJ = nnz(M);

[status, ttime] = lcppath(n+m, nnzJ, z, l, u, M, q);

if (status ~= 1)
status,
error('Path fails to solve problem');
end

mu = [];
if (m > 0)
mu = z(n+1:n+m);
z = z(1:n);
end

return;

Loading

0 comments on commit 4af540b

Please sign in to comment.