## Assignment 8 (posted on March 11)

### Re: Assignment 8 (posted on March 11)

I'm getting the error that my jacobian is singular to working precision. I use this function for my driving DP1 constraint:

constraint DP1

body_i ground

abar_i 0 1 0

body_j pendulum

abar_j 1 0 0

funct cos(pi/4*sin(2*t)-pi/2)

with the initial pendulum in this configuration

pendulum

r 0 0 -2

p 1 0 0 0

if I have p 0.5 0.5 0.5 0.5, then the newtonRaphson solver blows up (correction gets larger)

constraint DP1

body_i ground

abar_i 0 1 0

body_j pendulum

abar_j 1 0 0

funct cos(pi/4*sin(2*t)-pi/2)

with the initial pendulum in this configuration

pendulum

r 0 0 -2

p 1 0 0 0

if I have p 0.5 0.5 0.5 0.5, then the newtonRaphson solver blows up (correction gets larger)

### Re: Assignment 8 (posted on March 11)

@ Tyler

At time t=1 I get 0 1.31 -1.5113, but what did you put in for p?

At time t=1 I get 0 1.31 -1.5113, but what did you put in for p?

### Re: Assignment 8 (posted on March 11)

For my perpendicular constraint (two DP1's) I have the global x and pendulum's x orthogonal to one another. I have the global x and pendulum's y orthogonal to one another. For x CD, i have the constant matrix as [1 0 0] and sbarQ as [-2 0 0]. For y and z CD I just change the constant matrix to [0 1 0] and [0 0 1].

### Re: Assignment 8 (posted on March 11)

Sorry, I was getting some food.

Your constraint looks good. Rob was having a similar problem. The only difference I have is that my f(t) is cos(pi/2-pi/4*sin(2*t)); I don't think that should change much.

For t=1, if you calculate r = [0 1.31 -1.5113], you can form an A matrix then assume e0 NOT = 0. That leads you to some equations for e0, e1, e2, e3.

Lecture 0211 slides 12 (e0) and 13 (e1,e2,e3)

p = [0.662503 0.247162 0.662503 0.247162]

Your constraint looks good. Rob was having a similar problem. The only difference I have is that my f(t) is cos(pi/2-pi/4*sin(2*t)); I don't think that should change much.

For t=1, if you calculate r = [0 1.31 -1.5113], you can form an A matrix then assume e0 NOT = 0. That leads you to some equations for e0, e1, e2, e3.

Lecture 0211 slides 12 (e0) and 13 (e1,e2,e3)

p = [0.662503 0.247162 0.662503 0.247162]

### Re: Assignment 8 (posted on March 11)

So... I eliminated all my inputs except for an initial r and p.

I didn't ask for an initial guess for p_dot.

You only need a rotational velocity for gamma and I just wait and calculate gamma after newton-raphson is done, so I know I have an accurate A matrix.

Is that alright?

I didn't ask for an initial guess for p_dot.

You only need a rotational velocity for gamma and I just wait and calculate gamma after newton-raphson is done, so I know I have an accurate A matrix.

Is that alright?

### Re: Assignment 8 (posted on March 11)

Christopher - let me know if you figure out what was making your matrix singular. For the sake of moving on with my code, I'm using that pinv that Anne mentioned and then hopefully I'll go back and play around with inputs and get it to work without pinv.

### Re: Assignment 8 (posted on March 11)

I had put in a -I matrix for Ai when body i is ground. It's stemming from your Jacobian being singular (det = 0).

### Re: Assignment 8 (posted on March 11)

Can we solve velocity using pi in stead of pibar? That way it is expressed in the G-RF in which we can use the equation velQ = velP + w_pendulum x rPQ, where velP is the velocity of the COM of the pendulum moving in the G-RF, w_pendulum is the omega vector expressed in the G-RF (which only has entries in the first row due to the axis of rotation being about the global x-axis), and rPQ the vector between the COM of the pendulum and point Q in the L-RF (-2, 0, 0). Can someone verify this?

### Re: Assignment 8 (posted on March 11)

ME751Chris wrote:Can we solve velocity using pi in stead of pibar? That way it is expressed in the G-RF in which we can use the equation velQ = velP + w_pendulum x rPQ, where velP is the velocity of the COM of the pendulum moving in the G-RF, w_pendulum is the omega vector expressed in the G-RF (which only has entries in the first row due to the axis of rotation being about the global x-axis), and rPQ the vector between the COM of the pendulum and point Q in the L-RF (-2, 0, 0). Can someone verify this?

Thats the way I did it.

### Re: Assignment 8 (posted on March 11)

My velocities turned out to be zero for Q, but my accelerations were non-zero (specifically in the 'y' and 'z' direction...'x' direction was pretty much zero). Any reason why? I'm using the equation on slide 8 of 2/16.

Last edited by ME751Chris on Wed Mar 17, 2010 12:44 am, edited 1 time in total.

### Re: Assignment 8 (posted on March 11)

ME751Tyler wrote:Your constraint looks good. Rob was having a similar problem. The only difference I have is that my f(t) is cos(pi/2-pi/4*sin(2*t)); I don't think that should change much.

For t=1, if you calculate r = [0 1.31 -1.5113], you can form an A matrix then assume e0 NOT = 0. That leads you to some equations for e0, e1, e2, e3.

Lecture 0211 slides 12 (e0) and 13 (e1,e2,e3)

p = [0.662503 0.247162 0.662503 0.247162]

Actually, what Tyler did makes a difference and I had a posting in which i mentioned that you would have to go that way. If you use that f(t) it means that your a_i and a_j vectors are preprendicular at time t=0. Here is the posting that i'm talking about, it's an earlier post on this thread:

Dan Negrut wrote:ME751Anne wrote:It says it was singular to a working precision. To get around this I just used pinv (the pseudoinverse), and no problems arise now.

Anne - there is a problem if you choose the axes for the motion to coincide at time t=0. It's too long to explain (I'm working on tomorrow's notes) but the idea is to have for the DP1 motion one axis along the global OY axis, and another axis along the O'x' axis of the pendulum. The Jacobian should be healthy and there should be no need for the pseudoinverse.

Jim came over for office hours today and explained to him. He said he was going to post a more elaborate explanation of what's going on.

If you have the time please to check this please confirm here that it solved your problem.

dan

Again, take a_i to be on ground and parallel with OY axis (as such it'll be [0,1,0]'). Then take {\overbar a_j} to be a unit vector along the O'x' axis of the pendulum; i.e., [1,0,0]. Then the motion will look like Tyler's, and things should work out.

There should be no reason to use the "pinv" that Anne was mentioning.

Can one of you implement this and confirm that it removed the singularity?

Dan

Last edited by Anonymous on Wed Mar 17, 2010 6:57 am, edited 1 time in total.

### Re: Assignment 8 (posted on March 11)

ME751Chris wrote:I'm getting the error that my jacobian is singular to working precision. I use this function for my driving DP1 constraint:

constraint DP1

body_i ground

abar_i 0 1 0

body_j pendulum

abar_j 1 0 0

funct cos(pi/4*sin(2*t)-pi/2)

with the initial pendulum in this configuration

pendulum

r 0 0 -2

p 1 0 0 0

if I have p 0.5 0.5 0.5 0.5, then the newtonRaphson solver blows up (correction gets larger)

Chris - recall that for Newton-Raphson to converge it needs a good starting point. Your starting points are not that good.

Take a look at the initial configuration of the pendulum at t=0. It hangs down. Consequently, the A matrix looks like this:

A =

[ 0 0 1

0 1 0

-1 0 0]

As such, e_0 = sqrt(2)/2, e_1=0, e_2=sqrt(2)/2, e_3=0.

Your r is good.

The immediate question is this: so what should be the initial guess for Newton-Raphson at the next step. A good guess is this: always take as initial guess the configuration at the previous time step. In fact, if I remember correctly this is what the Newton-Raphson code the was provided in the slides uses.

Please comment on this once you get to implement it.

Dan

### Re: Assignment 8 (posted on March 11)

ME751Chris wrote:My velocities turned out to be zero for Q, but my accelerations were non-zero (specifically in the 'y' and 'z' direction...'x' direction was pretty much zero). Any reason why? I'm using the equation on slide 8 of 2/16.

I had a similar problem for a while. I was calculating gamma before I knew what omega was. (I was using initial guess values all the time) I pulled my gamma calculation out of my newton-raphson and computed gamma after i had the velocities. That way I could use the correct omega values.

### Re: Assignment 8 (posted on March 11)

ME751Tyler wrote:ME751Chris wrote:My velocities turned out to be zero for Q, but my accelerations were non-zero (specifically in the 'y' and 'z' direction...'x' direction was pretty much zero). Any reason why? I'm using the equation on slide 8 of 2/16.

I had a similar problem for a while. I was calculating gamma before I knew what omega was. (I was using initial guess values all the time) I pulled my gamma calculation out of my newton-raphson and computed gamma after i had the velocities. That way I could use the correct omega values.

Tyler - have you got correct values for the acceleration?

In general, keep this in mind:

- compute the ingredients for the velocity equation after you came up with the new position.

- compute the ingredients for the acceleration equation after you came up with the new position and new velocity.

dan

40 posts
• Page

**2**of**3**• 1,**2**, 3Return to ME751 Spring 2010: Advanced Computational Multibody Dynamics

### Who is online

Users browsing this forum: No registered users and 1 guest