Meaningless results when implementing second order derivatives using Tuple notation in pyomo.DAE

[samphillips38]

#When implementing second order derivatives using the (model.t, model.t) notation, I do not get reasonable results. Here's an example with the cart-pole swing-up problem. Implementing the second order derivatives as first order with respect to the first order derivative produces a reasonable result:

# Define constants
m1 = 1
m2 = 0.3
g = 9.81

# Problem global parameters
d = 1.0
d_max = 1.0
u_max = 20.0

# Construct model
model = pyo.ConcreteModel("Cart-Pole Swing-up")

# Parameters
model.T = pyo.Param(default=2)
model.l = pyo.Param(default=0.5)

# Create time
model.t = pyoDAE.ContinuousSet(bounds=(0, model.T))

# Create State variables
model.q1 = pyo.Var(model.t, bounds=(-d_max, d_max))
model.q2 = pyo.Var(model.t)

# Create control variables
model.u = pyo.Var(model.t, bounds=(-u_max, u_max))

# Create derivatives
model.q1_dot = pyoDAE.DerivativeVar(model.q1, wrt=model.t)
model.q2_dot = pyoDAE.DerivativeVar(model.q2, wrt=model.t)

model.q1_ddot = pyoDAE.DerivativeVar(model.q1, wrt=(model.t, model.t)) # This doesn't work
model.q2_ddot = pyoDAE.DerivativeVar(model.q2, wrt=(model.t, model.t)) # This doesn't work
# model.q1_ddot = pyoDAE.DerivativeVar(model.q1_dot, wrt=model.t) # This works
# model.q2_ddot = pyoDAE.DerivativeVar(model.q2_dot, wrt=model.t) # This works

# Initial and final constraints
def _init(model):
    yield model.q1[0] == 0
    yield model.q2[0] == 0
    yield model.q1_dot[0] == 0
    yield model.q2_dot[0] == 0

def _final(model):
    yield model.q1[model.T] == d
    yield model.q2[model.T] == np.pi
    yield model.q1_dot[model.T] == 0
    yield model.q2_dot[model.T] == 0

model.initial_constraints = pyo.ConstraintList(rule=_init)
model.final_constraints = pyo.ConstraintList(rule=_final)

def dynamic_constraint_1(model, t):
    q1_ddot = (model.l*m2*pyo.sin(model.q2[t])*model.q2_dot[t]**2 + model.u[t] + m2*g*pyo.cos(model.q2[t])*pyo.sin(model.q2[t])) / (m1 + m2*(1 - pyo.cos(model.q2[t])**2))
    return q1_ddot == model.q1_ddot[t]

def dynamic_constraint_2(model, t):
    q2_ddot = - (model.l*m2*pyo.cos(model.q2[t])*pyo.sin(model.q2[t])*model.q2_dot[t]**2 + model.u[t]*pyo.cos(model.q2[t]) + (m1 + m2)*g*pyo.sin(model.q2[t])) / (model.l*m1 + model.l*m2*(1 - pyo.cos(model.q2[t])**2))
    return q2_ddot == model.q2_ddot[t]

model.dynamic_constraint_1 = pyo.Constraint(model.t, rule=dynamic_constraint_1)
model.dynamic_constraint_2 = pyo.Constraint(model.t, rule=dynamic_constraint_2)

# New Objective function method - Differential equation in J
model.J = pyo.Var(model.t)
model.dJdt = pyoDAE.DerivativeVar(model.J, wrt=model.t)
model.dJdt_constraint = pyo.Constraint(model.t, rule=lambda model, t: model.dJdt[t] == model.u[t]**2)

model.objective = pyo.Objective(expr=model.J[model.T] - model.J[0])

discretizer = pyo.TransformationFactory('dae.collocation')
discretizer.apply_to(model, wrt=model.t, nfe=30, ncp=4, scheme='LAGRANGE-RADAU')

opt = pyo.SolverFactory("ipopt")
opt.solve(model, tee=True)
model.display()

Here's some plots to show a comparison of the outputs: