Error message from Problem

[Peter230655]

I am trying to optimize a simulation with variable h, but I am getting the long error message below.
My environment, it uses opty 1.3.0 dev0, runs other simulations, e.g. plot_pendulum_swing_up_variable_duration.py fine.
It does not look to me like the error I get sometimes, realted to the 'C translation'. of course not sure.
I have no idea, where I should look for the error.
Any help is greatly appreciated!

I have added the simulation below this error.

PrintMethodNotImplementedError            Traceback (most recent call last)
File c:\Users\Peter\github_repos\opty\examples-gallery\plot_Mississippi_steamboat.py:50
     43 limit_torque = 10.
     44 bounds = {
     45     tLW: (-limit_torque, limit_torque),
     46     tRW: (-limit_torque, limit_torque),
     47     h: (0.0, 1.0),
     48     }
---> 50 prob = Problem(
     51     obj,
     52     obj_grad,
     53     eom,
     54     state_symbols,
     55     num_nodes,
     56     interval_value,
     57     time_symbol=t,
     58     known_parameter_map=par_map,
     59     instance_constraints=instance_constraints,
     60     bounds=bounds,
     61     )

File c:\users\peter\github_repos\opty\opty\direct_collocation.py:53, in _DocInherit.get_with_inst.<locals>.f(*args, **kwargs)
     51 @wraps(self.mthd, assigned=('__name__', '__module__'))
     52 def f(*args, **kwargs):
---> 53     return self.mthd(obj, *args, **kwargs)

File c:\users\peter\github_repos\opty\opty\direct_collocation.py:163, in Problem.__init__(self, obj, obj_grad, equations_of_motion, state_symbols, num_collocation_nodes, node_time_interval, known_parameter_map, known_trajectory_map, instance_constraints, time_symbol, tmp_dir, integration_method, parallel, bounds, show_compile_output)
    161 self.con = self.collocator.generate_constraint_function()
    162 logging.info('Constraint function generated.')
--> 163 self.con_jac = self.collocator.generate_jacobian_function()
    164 logging.info('Jacobian function generated.')
    166 self.con_jac_rows, self.con_jac_cols = \
    167     self.collocator.jacobian_indices()

File c:\users\peter\github_repos\opty\opty\direct_collocation.py:1600, in ConstraintCollocator.generate_jacobian_function(self)
   1597 """Returns a function which evaluates the Jacobian of the
   1598 constraints given the array of free optimization variables."""
   1599 logging.info('Generating jacobian function.')
-> 1600 self._gen_multi_arg_con_jac_func()
   1601 return self._wrap_constraint_funcs(self._multi_arg_con_jac_func, 'jac')

File c:\users\peter\github_repos\opty\opty\direct_collocation.py:1402, in ConstraintCollocator._gen_multi_arg_con_jac_func(self)
   1398 # This generates a numerical function that evaluates the matrix of
   1399 # partial derivatives. This function returns the non-zero elements
   1400 # needed to build the sparse constraint Jacobian.
   1401 logging.info('Compiling the Jacobian function.')
-> 1402 eval_partials = ufuncify_matrix(args, symbolic_partials,
   1403                                 const=constant_syms + (h_sym,),
   1404                                 tmp_dir=self.tmp_dir,
   1405                                 parallel=self.parallel)
   1407 if isinstance(symbolic_partials, tuple) and len(symbolic_partials) == 2:
   1408     num_rows = symbolic_partials[1][0].shape[0]

File c:\users\peter\github_repos\opty\opty\utils.py:647, in ufuncify_matrix(args, expr, const, tmp_dir, parallel, show_compile_output)
    644 else:
    645     sub_exprs, simple_mat = sm.cse(expr, sm.numbered_symbols('z_'))
--> 647 sub_expr_code = '\n'.join(['double ' + ccode(sub_expr[1], sub_expr[0])
    648                            for sub_expr in sub_exprs])
    650 matrix_code = ccode(simple_mat[0], matrix_sym)
    652 d['eval_code'] = '    ' + '\n    '.join((sub_expr_code + '\n' +
    653                                          matrix_code).split('\n'))

File c:\users\peter\github_repos\opty\opty\utils.py:57, in ccode(expr, assign_to, **settings)
     55 def ccode(expr, assign_to=None, **settings):
     56     """Mimics SymPy's ccode, but uses our printer."""
---> 57     return OptyC99CodePrinter(settings).doprint(expr, assign_to)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\codeprinter.py:172, in CodePrinter.doprint(self, expr, assign_to)
    169 self._not_supported = set()
    170 self._number_symbols = set()
--> 172 lines = self._print(expr).splitlines()
    174 # format the output
    175 if self._settings["human"]:

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\codeprinter.py:400, in CodePrinter._print_Assignment(self, expr)
    398 else:
    399     lhs_code = self._print(lhs)
--> 400     rhs_code = self._print(rhs)
    401     return self._get_statement("%s = %s" % (lhs_code, rhs_code))

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\str.py:57, in StrPrinter._print_Add(self, expr, order)
     55 l = []
     56 for term in terms:
---> 57     t = self._print(term)
     58     if t.startswith('-') and not term.is_Add:
     59         sign = "-"

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\c.py:142, in _as_macro_if_defined.<locals>._meth_wrapper(self, expr, **kwargs)
    140     return '%s%s' % (self.math_macros[expr], self._get_math_macro_suffix(real))
    141 else:
--> 142     return meth(self, expr, **kwargs)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\c.py:272, in C89CodePrinter._print_Mul(self, expr, **kwargs)
    270 @_as_macro_if_defined
    271 def _print_Mul(self, expr, **kwargs):
--> 272     return super()._print_Mul(expr, **kwargs)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\codeprinter.py:565, in CodePrinter._print_Mul(self, expr)
    563     a_str = [self.parenthesize(a[0], 0.5*(PRECEDENCE["Pow"]+PRECEDENCE["Mul"]))]
    564 else:
--> 565     a_str = [self.parenthesize(x, prec) for x in a]
    566 b_str = [self.parenthesize(x, prec) for x in b]
    568 # To parenthesize Pow with exp = -1 and having more than one Symbol

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\str.py:36, in StrPrinter.parenthesize(self, item, level, strict)
     34 def parenthesize(self, item, level, strict=False):
     35     if (precedence(item) < level) or ((not strict) and precedence(item) <= level):
---> 36         return "(%s)" % self._print(item)
     37     else:
     38         return self._print(item)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\str.py:57, in StrPrinter._print_Add(self, expr, order)
     55 l = []
     56 for term in terms:
---> 57     t = self._print(term)
     58     if t.startswith('-') and not term.is_Add:
     59         sign = "-"

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\c.py:142, in _as_macro_if_defined.<locals>._meth_wrapper(self, expr, **kwargs)
    140     return '%s%s' % (self.math_macros[expr], self._get_math_macro_suffix(real))
    141 else:
--> 142     return meth(self, expr, **kwargs)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\c.py:272, in C89CodePrinter._print_Mul(self, expr, **kwargs)
    270 @_as_macro_if_defined
    271 def _print_Mul(self, expr, **kwargs):
--> 272     return super()._print_Mul(expr, **kwargs)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\codeprinter.py:565, in CodePrinter._print_Mul(self, expr)
    563     a_str = [self.parenthesize(a[0], 0.5*(PRECEDENCE["Pow"]+PRECEDENCE["Mul"]))]
    564 else:
--> 565     a_str = [self.parenthesize(x, prec) for x in a]
    566 b_str = [self.parenthesize(x, prec) for x in b]
    568 # To parenthesize Pow with exp = -1 and having more than one Symbol

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\str.py:36, in StrPrinter.parenthesize(self, item, level, strict)
     34 def parenthesize(self, item, level, strict=False):
     35     if (precedence(item) < level) or ((not strict) and precedence(item) <= level):
---> 36         return "(%s)" % self._print(item)
     37     else:
     38         return self._print(item)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\str.py:57, in StrPrinter._print_Add(self, expr, order)
     55 l = []
     56 for term in terms:
---> 57     t = self._print(term)
     58     if t.startswith('-') and not term.is_Add:
     59         sign = "-"

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\c.py:142, in _as_macro_if_defined.<locals>._meth_wrapper(self, expr, **kwargs)
    140     return '%s%s' % (self.math_macros[expr], self._get_math_macro_suffix(real))
    141 else:
--> 142     return meth(self, expr, **kwargs)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\c.py:272, in C89CodePrinter._print_Mul(self, expr, **kwargs)
    270 @_as_macro_if_defined
    271 def _print_Mul(self, expr, **kwargs):
--> 272     return super()._print_Mul(expr, **kwargs)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\codeprinter.py:565, in CodePrinter._print_Mul(self, expr)
    563     a_str = [self.parenthesize(a[0], 0.5*(PRECEDENCE["Pow"]+PRECEDENCE["Mul"]))]
    564 else:
--> 565     a_str = [self.parenthesize(x, prec) for x in a]
    566 b_str = [self.parenthesize(x, prec) for x in b]
    568 # To parenthesize Pow with exp = -1 and having more than one Symbol

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\str.py:36, in StrPrinter.parenthesize(self, item, level, strict)
     34 def parenthesize(self, item, level, strict=False):
     35     if (precedence(item) < level) or ((not strict) and precedence(item) <= level):
---> 36         return "(%s)" % self._print(item)
     37     else:
     38         return self._print(item)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\str.py:57, in StrPrinter._print_Add(self, expr, order)
     55 l = []
     56 for term in terms:
---> 57     t = self._print(term)
     58     if t.startswith('-') and not term.is_Add:
     59         sign = "-"

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\printer.py:331, in Printer._print(self, expr, **kwargs)
    329     printmethod = getattr(self, printmethodname, None)
    330     if printmethod is not None:
--> 331         return printmethod(expr, **kwargs)
    332 # Unknown object, fall back to the emptyPrinter.
    333 return self.emptyPrinter(expr)

File c:\users\peter\github_repos\opty\opty\utils.py:51, in OptyC99CodePrinter._print_Function(self, expr)
     50 def _print_Function(self, expr):
---> 51     name = super()._print_Function(expr)
     52     return name + '_'

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\codeprinter.py:461, in CodePrinter._print_Function(self, expr)
    459     return '%s(%s)' % (self._print(expr.func), ', '.join(map(self._print, expr.args)))
    460 else:
--> 461     return self._print_not_supported(expr)

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\sympy\printing\codeprinter.py:582, in CodePrinter._print_not_supported(self, expr)
    580 def _print_not_supported(self, expr):
    581     if self._settings.get('strict', False):
--> 582         raise PrintMethodNotImplementedError("Unsupported by %s: %s" % (str(type(self)), str(type(expr))) + \
    583                          "\nSet the printer option 'strict' to False in order to generate partially printed code.")
    584     try:
    585         self._not_supported.add(expr)

PrintMethodNotImplementedError: Unsupported by <class 'opty.utils.OptyC99CodePrinter'>: re
Set the printer option 'strict' to False in order to generate partially printed code.

This is the program up to the error:

"""
Mississippi Steamboat
=====================

A boat is modeled as a rectangular with length :math:`a_S` and width :math:`b_S`.
it has a mass :math:`m_S` and is modeled as a rigid body.
Water wheels are attached to the boat on the left and right side.
The wheels have radius :math:`r_W` and mass :math:`m_W` and are modeled as
rigid bodies.
By running the wheels at different speeds, the boat can be steered.
The water speed is assumed to be zero, gravity is disregarded.


**Constants**

- :math:`m_S`: mass of the steamboat [kg]
- :math:`m_W`: mass of the wheels [kg]
- :math:`r_W`: radius of the wheel [m]
- :math:`a_S`: length of the steamboat [m]
- :math:`b_S`: width of the steamboat [m]
- :math:`c_W`: drag coefficient at wheels [kg/m^2]
- :math:`c_S`: drag coefficient at steamboat [kg/m^2]


**States**

- :math:`x`: X - position of the center of the steamboat [m]
- :math:`y`: Y - position of the center of the steamboat [m]
- :math:`q_S`: angle of the steamboat [rad]
- :math:`q_{LW}`: angle of the left wheel [rad]
- :math:`q_{RW}`: angle of the right wheel [rad]
- :math:`u_x`: speed of the steamboat in X direction [m/s]
- :math:`u_y`: speed of the steamboat in Y direction [m/s]
- :math:`u_S`: angular speed of the steamboat [rad/s]
- :math:`u_{LW}`: angular speed of the left wheel [rad/s]
- :math:`u_{RW}`: angular speed of the right wheel [rad/s]

**Specifieds**

- :math:`t_{LW}`: torque applied to the left wheel [Nm]
- :math:`t_{RW}`: torque applied to the right wheel [Nm]

"""
import numpy as np
import sympy as sm
from opty.utils import create_objective_function
import sympy.physics.mechanics as me
from opty.direct_collocation import Problem
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


# %%
# Set up the geometry of the system.
#
# - :math:`N`: inertial frame of reference
# - :math:`O`: origin of the inertial frame of reference
# - :math:`A_S`: body fixed frame of the steamboat
# - :math:`A_{LW}`: body fixed frame of the left wheel
# - :math:`A_{RW}`: body fixed frame of the right wheel
# - :math:`Dmc_S`: center of the steamboat
# - :math:`Dmc_{LW}`: center of the left wheel
# - :math:`Dmc_{RW}`: center of the right wheel
# - :math:`FP_{LW}`: point where the force due to rotation of the left
#   wheel is applied
# - :math:`FP_{RW}`: point where the force due to rotation of the right
#   wheel is applied

N, AS, ALW, ARW = sm.symbols('N, AS, ALW, ARW', cls = me.ReferenceFrame)
DmcS, DmcLW, DmcRW = sm.symbols('DmcS, DmcLW, DmcRW', cls = me.Point)
O, FPLW, FPRW = sm.symbols('O, FPLW, FPRW', cls = me.Point)

q, x, y, qLW, qRW = me.dynamicsymbols('q, x, y, qLW, qRW')
u, ux, uy, uLW, uRW = me.dynamicsymbols('u, ux, uy, uLW, uRW')
tLW, tRW = me.dynamicsymbols('tLW, tRW')

mS, mW, rW, aS, bS, cS, cW = sm.symbols('mS, mW, rW, aS, bS, cS, cW', real=True)

t = me.dynamicsymbols._t
O.set_vel(N, 0)

AS.orient_axis(N, q, N.z)
AS.set_ang_vel(N, u*N.z)
ALW.orient_axis(AS, qLW, AS.x)
ALW.set_ang_vel(AS, uLW*AS.x)
ARW.orient_axis(AS, qRW, AS.x)
ARW.set_ang_vel(AS, uRW*AS.x)

DmcS.set_pos(O, x*N.x + y*N.y)
DmcS.set_vel(N, ux*N.x + uy*N.y)
DmcLW.set_pos(DmcS, -1.1*bS*AS.x)
DmcLW.v2pt_theory(DmcS, N, AS)
DmcRW.set_pos(DmcS, 1.1*bS*AS.x)
DmcRW.v2pt_theory(DmcS, N, AS)

FPLW.set_pos(DmcLW, -rW*N.z)
FPLW.set_vel(N, DmcLW.vel(N) + (-uLW*AS.x.cross(-rW*N.z)))
FPRW.set_pos(DmcRW, -rW*N.z)
FPRW.set_vel(N, DmcRW.vel(N) + uRW*AS.x.cross(-rW*N.z))


# %%
# Set up the drag forces acting on the boat.
#
# The drag force acting on a body moving in a fluid is given by
# :math:`F_D = -c \cdot A \cdot | \bar v|^2 \cdot \dfrac{\bar v}{| \bar v|}`,
# where :math:`c` is the drag coefficient, :math:`\bar v` is the velocity
# vector of the body and :math:`A` is the cross section area of the body facing
# the flow. This may be found here:
#
# https://courses.lumenlearning.com/suny-physics/chapter/5-2-drag-forces/
#
# In order to avoid the issues with .normalize(), I will use the following:
#
# ( :math:`\circ` denotes the inner product )
#
# :math:`F_{D_x} = -c \cdot A \cdot \sqrt{(A.x \circ \bar v)^2 )} \cdot A.x \circ \bar v \cdot A.x`
# :math:`F_{D_y} = -c \cdot A \cdot \sqrt{(A.y \circ \bar v)^2 )} \cdot A.y \circ \bar v \cdot A.y`

helpx = DmcS.vel(N).dot(AS.x)
helpy = DmcS.vel(N).dot(AS.y)

FDx = -cS*aS*sm.sqrt(helpx**2)*helpx*AS.x
FDy = -cS*bS*sm.sqrt(helpy**2)*helpy*AS.y
forces = [(DmcS, FDx + FDy)]

# %%
# Set up the forces acting on the wheels.
# The drag forces are similar to above, except that these forces act only in the AS.y
# direction.

helpy = FPLW.vel(N).dot(AS.y)
FLW = -cW*rW*sm.sqrt(helpy**2)*helpy*AS.y

helpy = FPRW.vel(N).dot(AS.y)
FRW = -cW*rW*sm.sqrt(helpy**2)*helpy*AS.y
forces.append((DmcLW, FLW))
forces.append((DmcRW, FRW))

forces.append((ALW, tLW*AS.x - FLW.cross(rW*N.z)))
forces.append((ARW, tRW*AS.x - FRW.cross(rW*N.z)))

# %%
# If :math:`u \neq 0`, the boat will rotate and a drag torque will act on it.

Tdrag = -cS*aS * u*AS.z.cross(AS.y)
forces.append((AS, Tdrag))

# %%
# Set up the rigid bodies

iXX = 0.5*mW*rW**2
iYY = 0.25*mW*rW**2
iZZ = iYY
J1 = me.inertia(ALW, iXX, iYY, iZZ)
J2 = me.inertia(ARW, iXX, iYY, iZZ)
left_wheel = me.RigidBody('left_wheel', DmcLW, ALW, mW, (J1, DmcLW))
right_wheel = me.RigidBody('right_wheel', DmcRW, ARW, mW, (J2, DmcRW))

iZZ = 1/12 * mS*(aS**2 + bS**2)
J3 = me.inertia(AS, 0, 0, iZZ)
boat = me.RigidBody('boat', DmcS, AS, mS, (J3, DmcS))

bodies = [boat, left_wheel, right_wheel]

# %%
# Set up Kame's equations of motion.

q_ind = [q, x, y, qLW, qRW]
u_ind = [u, ux, uy, uLW, uRW]
kd = sm.Matrix([i - j.diff(t) for j, i in zip(q_ind, u_ind)])

KM = me.KanesMethod(N,
                    q_ind=q_ind,
                    u_ind=u_ind,
                    kd_eqs=kd,
                    )

fr, frstar = KM.kanes_equations(bodies, forces)
eom = kd.col_join(fr + frstar)

print(f'the eoms contain {sm.count_ops(eom):,} operations \n')
print(me.find_dynamicsymbols(eom), '\n')
print(eom.free_symbols, '\n')

# %%
# Set up the optimization problem and solve it.
# --------------------------------------------
#

state_symbols = [q, x, y, qLW, qRW, u, ux, uy, uLW, uRW]
specified_symbols = [tLW, tRW]
constant_symbols = [mS, mW, rW, aS, bS, cS, cW]
num_nodes = 250
h = sm.symbols('h')

# %%
# Specify the known symbols.

par_map = {}
par_map[mS] = 10
par_map[mW] = 1.0
par_map[rW] = 1.0
par_map[aS] = 10
par_map[bS] = 1.0
par_map[cS] = 1.0
par_map[cW] = 1.0

# %%
# Set up the objective function and its gradient.

def obj(free):
    return free[-1]

def obj_grad(free):
    grad = np.zeros_like(free)
    grad[-1] = 1.0
    return grad

duration = (num_nodes - 1)*h
t0, tf = 0.0, duration
interval_value = h


# %%
# Set up the instance constraints the bounds and Problem.

initial_state_constraints = {
    q: -np.pi/2.,
    x: 0.0,
    y: 0.0,
    qLW: 0.0,
    qRW: 0.0,
    u: 0.0,
    ux: 0.0,
    uy: 0.0,
    uLW: 0.0,
    uRW: 0.0,
    }

final_state_constraints = {
    q: -np.pi/2.,
    x: 10,
    y: 10,
    ux: 0.0,
    uy: 0.0,
    }

instance_constraints = (
    q.subs({t: t0}) - initial_state_constraints[q],
    x.subs({t: t0}) - initial_state_constraints[x],
    y.subs({t: t0}) - initial_state_constraints[y],
    qLW.subs({t: t0}) - initial_state_constraints[qLW],
    qRW.subs({t: t0}) - initial_state_constraints[qRW],
    u.subs({t: t0}) - initial_state_constraints[u],
    ux.subs({t: t0}) - initial_state_constraints[ux],
    uy.subs({t: t0}) - initial_state_constraints[uy],
    uLW.subs({t: t0}) - initial_state_constraints[uLW],
    uRW.subs({t: t0}) - initial_state_constraints[uRW],
    q.subs({t: tf}) - final_state_constraints[q],
    x.subs({t: tf}) - final_state_constraints[x],
    y.subs({t: tf}) - final_state_constraints[y],
    ux.subs({t: tf}) - final_state_constraints[ux],
    uy.subs({t: tf}) - final_state_constraints[uy],
    )

limit_torque = 10.
bounds = {
    tLW: (-limit_torque, limit_torque),
    tRW: (-limit_torque, limit_torque),
    h: (0.0, 1.0),
    }

prob = Problem(
    obj,
    obj_grad,
    eom,
    state_symbols,
    num_nodes,
    interval_value,
    time_symbol=t,
    known_parameter_map=par_map,
    instance_constraints=instance_constraints,
    bounds=bounds,
    )
'''

[Peter230655]

The problem was due to my using $\sqrt{a}$ somewhere, with $a \geq 0$ always.
I must admit that the error message sent me on a wild goose chase for a few hours - but no issue, I have time!