Big error message.

[Peter230655]

I am trying to apply opty to my woodpecker simulation. (The normal simulation can be found on pst-notebooks) but I am getting the error below:

When setting up the forces, I use sm.piecewise((...), (..)) which in some earlier simulation did not affect opty.

I attached the program below the error. (I do not know, why it looks so ugly. I used ``` as a first and last line, as usual)

Thanks for any help!!

--------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
File c:\Users\Peter\anaconda3\envs\math\Lib\site-packages\opty\utils.py:701, in ufuncify_matrix(args, expr, const, tmp_dir, parallel, show_compile_output)
    700 try:
--> 701     cython_module = importlib.import_module(d['file_prefix'])
    702 except ImportError as error:

File c:\Users\Peter\anaconda3\envs\math\Lib\importlib\__init__.py:126, in import_module(name, package)
    125         level += 1
--> 126 return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1204, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1176, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1140, in _find_and_load_unlocked(name, import_)

ModuleNotFoundError: No module named 'ufuncify_matrix_0'

The above exception was the direct cause of the following exception:

ImportError                               Traceback (most recent call last)
Cell In[2], line 149
    142 final_state_constraints    = {yR: -10}
    144 instance_constraints = (
    145     *tuple(xi.subs({t: t0}) - xi_val for xi, xi_val in initial_state_constraints.items()), 
    146     *tuple(xi.subs({t: tf}) - xi_val for xi, xi_val in final_state_constraints.items()),
    147     )
--> 149 prob = Problem(
    150     obj, 
    151     obj_grad,
    152     EOM, 
    153     state_symbols, 
    154     num_nodes, 
    155     interval_value,
    156     known_parameter_map=par_map,
    157     instance_constraints=instance_constraints,
    158     integration_method=methode)
    160 prob.add_option('max_iter', 3000)  # default is 3000
    162 # Find the optimal solution.

File c:\Users\Peter\anaconda3\envs\math\Lib\site-packages\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\anaconda3\envs\math\Lib\site-packages\opty\direct_collocation.py:161, 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)
    159 self.obj = obj
    160 self.obj_grad = obj_grad
--> 161 self.con = self.collocator.generate_constraint_function()
    162 logging.info('Constraint function generated.')
    163 self.con_jac = self.collocator.generate_jacobian_function()

File c:\Users\Peter\anaconda3\envs\math\Lib\site-packages\opty\direct_collocation.py:1582, in ConstraintCollocator.generate_constraint_function(self)
   1579 """Returns a function which evaluates the constraints given the
   1580 array of free optimization variables."""
   1581 logging.info('Generating constraint function.')
-> 1582 self._gen_multi_arg_con_func()
   1583 return self._wrap_constraint_funcs(self._multi_arg_con_func, 'con')

File c:\Users\Peter\anaconda3\envs\math\Lib\site-packages\opty\direct_collocation.py:1031, in ConstraintCollocator._gen_multi_arg_con_func(self)
   1028     adjacent_stop = None
   1030 logging.info('Compiling the constraint function.')
-> 1031 f = ufuncify_matrix(args, self.discrete_eom,
   1032                     const=constant_syms + (h_sym,),
   1033                     tmp_dir=self.tmp_dir, parallel=self.parallel,
   1034                     show_compile_output=self.show_compile_output)
   1036 def constraints(state_values, specified_values, constant_values,
   1037                 interval_value):
   1038     """Returns a vector of constraint values given all of the
   1039     unknowns in the equations of motion over the 2, ..., N time
   1040     steps.
   (...)
   1059 
   1060     """

File c:\Users\Peter\anaconda3\envs\math\Lib\site-packages\opty\utils.py:706, in ufuncify_matrix(args, expr, const, tmp_dir, parallel, show_compile_output)
    702     except ImportError as error:
    703         msg = ('Unable to import the compiled Cython module {}, '
    704                'compilation likely failed. STDERR output from '
    705                'compilation:\n{}')
--> 706         raise ImportError(msg.format(d['file_prefix'], stderr)) from error
    707 finally:
    708     module_counter += 1

ImportError: Unable to import the compiled Cython module ufuncify_matrix_0, compilation likely failed. STDERR output from compilation:
c:\Users\Peter\anaconda3\envs\math\Lib\site-packages\Cython\Compiler\Main.py:381: FutureWarning: Cython directive 'language_level' not set, using '3str' for now (Py3). This has changed from earlier releases! File: C:\Users\Peter\AppData\Local\Temp\tmpyw2qqq9u.opty_ufuncify_compile\ufuncify_matrix_0.pyx
  tree = Parsing.p_module(s, pxd, full_module_name)
error: command 'C:\\Program Files\\Microsoft Visual Studio\\2022\\Community\\VC\\Tools\\MSVC\\14.38.33130\\bin\\HostX86\\x64\\cl.exe' failed with exit code 2

this is the program:

-------------------------------

import sympy.physics.mechanics as me from collections import OrderedDict import time import numpy as np import sympy as sm from scipy.interpolate import CubicSpline

import opty from opty.direct_collocation import Problem from opty.utils import create_objective_function import matplotlib.pyplot as plt import matplotlib from IPython.display import HTML matplotlib.rcParams['animation.embed_limit'] = 2**128 from matplotlib.animation import FuncAnimation print('opty version: ', opty.version) print('sympy version: ', sm.version)

N, AR, AW = sm.symbols('N, AR, AW', cls = me.ReferenceFrame) O, DmcR, DmcW, P = sm.symbols('O, DmcR, DmcW, P', cls=me.Point)

t = me.dynamicsymbols._t

qR, uR, qW, uW = me.dynamicsymbols('qR, uR, qW, uW') yR, uyR = me.dynamicsymbols('yR, uyR')

mR, mW, g, iZZR, iZZW, ro, d, max_kipp, kR, kW, reibung, korrekt = sm.symbols('mR, mW, g, iZZR, iZZW, ro, d, max_kipp, kR, kW, reibung, korrekt')

set up the geometry

O.set_vel(N, 0)

AR.orient_axis(N, qR, N.z) AR.set_ang_vel(N, uR N.z) AW.orient_axis(AR, qW, AR.z) rot = AW.ang_vel_in(N) AW.set_ang_vel(AR, uW AR.z) rot1 = AW.ang_vel_in(N)

DmcR.set_pos(O, yRN.y) DmcR.set_vel(N, uyRN.y)

P.set_pos(DmcR, roAR.x - ro/5.AR.y) P.v2pt_theory(DmcR, N, AR)

DmcW.set_pos(P, d*AW.x) DmcW.v2pt_theory(P, N, AW)

iR = me.inertia(AR, 0., 0., iZZR) iW = me.inertia(AW, 0., 0., iZZW) Ring = me.RigidBody('Ring', DmcR, AR, mR, (iR, DmcR)) Woodpecker = me.RigidBody('Woodpecker', DmcW, AW, mW, (iW, DmcW)) BODY = [Ring, Woodpecker]

FG = [(DmcR, -mRgN.y), (DmcW, -mWgN.y)] # gravity FT1 = [(AW, -kW (qW-korrekt) N.z), (AR, kW (qW-korrekt) N.z)] # torque if woodpecker is not 'in line' with the ring

When - max_kipp > qR or qR > max_kipp, the ring touches the peg

faktor1 = sm.Piecewise((1., -max_kipp > qR), (1., qR >= max_kipp), (0., True))

the 'penetration' is qR - max_kipp if qR > max_kipp, and -(qR + max_kipp) if qR < -max_kipp

eindring = sm.Piecewise((qR - max_kipp, qR > max_kipp), (-qR - max_kipp, qR < -max_kipp), (0., True))

This penetration will create a torque on AR. I must determine the direction of this torque

faktor2 = sm.Piecewise((-1., max_kipp < qR), (1., True)) kraft = kR eindring faktor1 faktor2 FT2 = [(AR, kraft N.z)]

friction acting on DmcR to stop the sliding. I made it speed dependent to avoid that the ring might be driven upward.

kraft2 = sm.tanh(5 kraft) FR = [(DmcR, -uyR reibung kraft2 N.y)]
FL = FG + FT1 + FT2 + FR hilfs = [kraft, -uyR reibungkraft2]

set up the equations of motion

q_ind = [qR, qW, yR] u_ind = [uR, uW, uyR] kd = sm.Matrix([i - j.diff(t) for i, j in zip(u_ind, q_ind)])

KM = me.KanesMethod(N, q_ind=q_ind, u_ind=u_ind, kd_eqs=kd) (fr, frstar) = KM.kanes_equations(BODY, FL) MM = KM.mass_matrix_full force = KM.forcing_full

EOM = kd.col_join(fr + frstar)

print('EOM DS', me.find_dynamicsymbols(EOM)) print('EOM free symbols', EOM.free_symbols) print(f'EOM has {sm.count_ops(EOM)} operations, {sm.count_ops(sm.cse(EOM))} operations after cse', '\n')

state_symbols = tuple(q_ind + u_ind) laenge = len(state_symbols) constant_symbols = tuple((mR, mW, g, iZZR, iZZW, ro, d, max_kipp, korrekt)) specified_symbols = [] unknown_symbols = [kR, kW, reibung]

methode = 'backward euler'

duration = 10. # The duration of the simulation. num_nodes = 100 # The number of nodes to use in the direct collocation problem. t0, tf = 0., duration # The initial and final times interval_value = duration / (num_nodes - 1) # The interval value between nodes.

WERT = [0., 0., 0., 0., 0., 0.] initialguess = np.array([wert for wert in WERT for in range(num_nodes)] + [1.e6, 2.e3, 0.75])

Specify the known system parameters.

par_map = OrderedDict() par_map[mR] = 1.0 # Mass of the ring par_map[mW] = 15 # Mass of the woodpecker par_map[g] = 9.81 # Acceleration due to gravity par_map[iZZR] = 1. # Moment of inertia of the ring par_map[iZZW] = 10. # Moment of inertia of the woodpecker par_map[ro] = 1.0 # Distance of the mass center of the ring from DmcR par_map[d] = 1.0 # Distance of the mass center of the woodpecker from P par_map[max_kipp] = np.deg2rad(2.) # The maximum angle the ring can be tilted par_map[korrekt] = np.deg2rad(10.) # The angle the ring is tilted in the beginning

objective = sm.Integral(kR2 + kW2 + reibung**2, t) # minimal power

Specify the objective function and it's gradient.

obj, obj_grad = create_objective_function( objective, state_symbols, specified_symbols, unknown_symbols, num_nodes, interval_value, methode)

====================================================================================

initial_state_constraints = { qR: np.deg2rad(0), qW: np.deg2rad(-30), yR: 0., uR: 0., uW: 0., uyR: 0.}

final_state_constraints = {yR: -10}

instance_constraints = ( tuple(xi.subs({t: t0}) - xi_val for xi, xi_val in initial_state_constraints.items()), tuple(xi.subs({t: tf}) - xi_val for xi, xi_val in final_state_constraints.items()), )

prob = Problem( obj, obj_grad, EOM, state_symbols, num_nodes, interval_value, known_parameter_map=par_map, instance_constraints=instance_constraints, integration_method=methode)

prob.add_option('max_iter', 3000) # default is 3000

Find the optimal solution.

solution, info = prob.solve(initial_guess) print('message from optimizer:', info['status_msg']) print('Iterations needed',len(prob.obj_value)) print(f"objective value {info['obj_val']:.3e}")

print(f'kR: {solution[-3]:.3f}') print(f'kW: {solution[-2]:.3f}') print(f'reibung: {solution[-1]:.3f}')

prob.plot_constraint_violations(solution) prob.plot_trajectories(solution)

[Peter230655]

It seems to be related to piecewise, since when I remove it (the simulation makes no sense then) the error disappears. I will try to replace piecewise with combinations of 'differentiable step functions' to get around this issue - but maybe still interesting on its own?

[moorepants]

Piecewise is not differentiable.

[Peter230655]

I replaced the piecewise 'logic' using differentialble approximations of step functions $h(x,b) = \dfrac{1}{1 + e^{a(x - b)}}$ , and the error went away. I think, this should be possible in general, not 100% sure.