Topic: Newbie needs help

Hi all,

I'm new to both Python and Openopt, so I think calling me a newbie would be an understatement tongue
I've used AMPL (with IPOPT solver) in the past, but I need more power so I'm looking at OpenOpt.

I'm trying to use FuncDesigner and OpenOpt to automatically generate (and then optimize) a system of equations that describe the motion of some rigid-bodies in 2D over N discreet timesteps. The solution that ralg is returning (all zeros) is not what I expect. I think I'm making a mistake when I define the constraints - possibly in the way I'm using oovars in arrays? See my code:

from FuncDesigner import *
from openopt import NLP

from RigidBody import RigidBody

pTimeStepSize = 0.1 #0.1 seconds per timestep
nTimesteps = int(8) #8 timesteps
nDOF = 3    #2 translational, 1 rotational
pGlobalGravity = [0.0,-9.81,0.0] #9.81m/s**2 is gravity at Earth sea-level

#add some bodies
bodyList = []
bodyList.append(RigidBody(10.5,2.0)) #long thin rod; 10.5kg, 2.0m length
bodyList.append(RigidBody(5.5,2.0)) #long thin rod; 5.5kg, 2.0m length

#make constraints
constraints = []
q = [ [ [ None for t in range(nTimesteps)] for d in range(nDOF)] for b in range(len(bodyList))]

for b in xrange(len(bodyList)):
    for d in xrange(nDOF):
        for t in xrange(nTimesteps):
            q[b][d][t] = oovar('q['+str(b)+']['+str(d)+']['+str(t)+']')
        for t in xrange(1,nTimesteps-1):
            f = bodyList[b].mass[d]*pGlobalGravity[d] #force of gravity (f=ma)
            m = bodyList[b].mass[d]
            a = ((q[b][d][t+1] - 2 * q[b][d][t] + q[b][d][t-1])/pTimeStepSize**2) #2nd time derivative of q, by finite difference formula
            constraints.append(f==m*a) #Newton's second law

#start point
startPoint = {}
for b in xrange(len(bodyList)):
    for d in xrange(nDOF):
        for t in xrange(nTimesteps):
            startPoint[q[b][d][t]] = 0

#objective
objective = q[0][0][0]**2 #dummy objective, just for testing purposes

#solve it
p = NLP(objective, startPoint, constraints = constraints)
r = p.solve('ralg')
print(r.xf)

The result is:

new rigid body: length = 2.0 mass=[10.5, 10.5, 3.5]
new rigid body: length = 2.0 mass=[5.5, 5.5, 1.8333333333333333]
--------------------------------------------------
solver: ralg   problem: unnamed    type: NLP   goal: minimum
 iter    objFunVal    log10(maxResidual)   
    0  0.000e+00               2.01 
    1  0.000e+00            -100.00 
istop: 14(move direction has all-zero coords)
Solver:   Time Elapsed = 0.01         CPU Time Elapsed = 0.00716551751364
objFunValue: 0 (feasible, MaxResidual = 0)
{q[0][0][0]: array([ 0.]), q[0][0][1]: array([ 0.]), q[0][0][2]: array([ 0.]), q[0][0][3]: array([ 0.]), q[0][0][4]: array([ 0.]), q[0][0][5]: array([ 0.]), q[0][0][6]: array([ 0.]), q[0][0][7]: array([ 0.]), q[1][2][0]: array([ 0.]), q[1][2][1]: array([ 0.]), q[1][2][2]: array([ 0.]), q[1][2][3]: array([ 0.]), q[1][2][4]: array([ 0.]), q[1][2][5]: array([ 0.]), q[1][2][6]: array([ 0.]), q[1][2][7]: array([ 0.]), q[1][0][0]: array([ 0.]), q[1][0][1]: array([ 0.]), q[1][0][2]: array([ 0.]), q[1][0][3]: array([ 0.]), q[1][0][4]: array([ 0.]), q[1][0][5]: array([ 0.]), q[1][0][6]: array([ 0.]), q[1][0][7]: array([ 0.]), q[0][1][0]: array([ 0.]), q[0][1][1]: array([ 0.]), q[0][1][2]: array([ 0.]), q[0][1][3]: array([ 0.]), q[0][1][4]: array([ 0.]), q[0][1][5]: array([ 0.]), q[0][1][6]: array([ 0.]), q[0][1][7]: array([ 0.]), q[1][1][0]: array([ 0.]), q[1][1][1]: array([ 0.]), q[1][1][2]: array([ 0.]), q[1][1][3]: array([ 0.]), q[1][1][4]: array([ 0.]), q[1][1][5]: array([ 0.]), q[1][1][6]: array([ 0.]), q[1][1][7]: array([ 0.]), q[0][2][0]: array([ 0.]), q[0][2][1]: array([ 0.]), q[0][2][2]: array([ 0.]), q[0][2][3]: array([ 0.]), q[0][2][4]: array([ 0.]), q[0][2][5]: array([ 0.]), q[0][2][6]: array([ 0.]), q[0][2][7]: array([ 0.])}

Also, is there any way to print out the constraints that I'm generating (so I can do a sanity-check, and make sure they are what I *think* they are)?

Any help appreciated!

Thanks,
Spacecookies

Re: Newbie needs help

hi Spacecookies,

For the non-negaative objective involved (q[0][0][0]**2) ralg the zero result is correct (since all constraints are satisfied).

Also, I have no Python module RigidBody to perform any code tests.

If you'll have any other questions, please provide information in an appropriate format for an openopt.org/Applications entry.

Regards, D.

Re: Newbie needs help

Ah, my apologies... RigidBody is a pretty simple class that just stores mass and length (and calculates a moment of inertia):

class RigidBody:
    '''Represents a rigid body.'''

    def __init__(self, mass, length):
    '''Constructor'''
    self.length = length
    self.mass = [mass,mass,1.0/12.0 * mass * length**2]  #moment of inertia for long thin rod (from http://en.wikipedia.org/wiki/List_of_moments_of_inertia)
    print 'new rigid body: length = '+str(length)+' mass='+str(self.mass)

I'm pretty sure that the all zero solution should violate the constraints (it's impossible for a body to hover at the origin because gravity will accelerate it downwards). The same constraints in AMPL + IPOPT result in a parabolic trajectory for each body.

Is there anyway to debug the generated constraints? To print them perhaps, so I can make sure they are what I *think* they are?

please provide information in an appropriate format for an openopt.org/Applications entry

"Spacecookes, MSc. Candidate in the School of Interactive Arts + Technology (www.siat.sfu.ca) at Simon Fraser University, Canada, using OpenOpt to simulate the motion of rigid bodies"

Thanks!

Last edited by Spacecookies (2010-11-25 07:23:34)

Re: Newbie needs help

Ah I think I solved it...
If I change the objective from:

#objective
objective = q[0][0][0]**2

to:

#objective
objective = q[0][0][0]**2.0

then it works as I expect. Apparently Python was treating the objective as an integer type?

Re: Newbie needs help

No, using q[0][0][0]**2.0 instead of q[0][0][0]**2 just adds the constraint q[0][0][0] >= 0, and thus the results are different.
I think ralg works incorrectly, because all you have are some linear equality constraints, that are handled in a very special way (via modification of space dilation matrix), and your objective in start point has all-zero gradient, thus dilated gradient also has all-zero coords. I'll take a look what's wrong in current ralg implementation, currently I would recommend you either using another solver (e.g. ipopt works good on it, yielding not all-zero coords), or using another start point, where q[0][0][0] is non-zero. If you have windows or other problems with ipopt installation, you could try using our sage-server (see link from left side of our website).

Regards, D.

Re: Newbie needs help

I have committed the fix for ralg, now it works ok with the situation.

Re: Newbie needs help

Dmitrey wrote:

I have committed the fix for ralg, now it works ok with the situation.

Thanks Dmitrey! I think that fixed it smile


By the way, is it possible to avoid specifying a startpoint to the NLP() method? I don't particularly care (or even know!) where it should start, so automatically setting it to all zeros would be convenient for me.

Re: Newbie needs help

Spacecookies wrote:

By the way, is it possible to avoid specifying a startpoint to the NLP() method? I don't particularly care (or even know!) where it should start, so automatically setting it to all zeros would be convenient for me.

Maybe for your situation you don't care, but generally time elapsed for solution depends on start point very much. Currently you have to provide it, maybe in future it will be reworked.