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
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