bicycleIftommBenchmark.py
You can view and download this file on Github: bicycleIftommBenchmark.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: bicycle Iftomm benchmark example
5# https://www.iftomm-multibody.org/benchmark/problem/Uncontrolled_bicycle/
6#
7# Author: Johannes Gerstmayr
8# Date: 2021-06-22
9#
10# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
11#
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14
15import exudyn as exu
16from exudyn.itemInterface import *
17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
18import exudyn.graphics as graphics
19
20from math import sin, cos, pi
21import numpy as np
22
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25
26
27#%%++++++++++++++++++++++++++++++++++++++++++++++++
28#coordinate system according to IFToMM:
29#note: here, wheels are rotated as local wheel axis=x, z points upwards in EXUDYN model
30# ox P2
31# oooo o
32# oooo o
33# +++ oooo o +++
34# + + ooo o +
35# + oooo + o +
36# + xP1+ + xP3+
37# + + + +
38# + + + +
39# +++ +++
40# x---------------------------x----------------------> x
41# | <- w ->
42# v z
43
44
45#parameters
46sZ = -1 #switch z coordinate compared to IFToMM description
47w = 1.02 #wheel base (distance of wheel centers)
48c = 0.08 #trail
49lam = pi/10 #steer axis tilt (rad)
50g = [0,0,9.81*sZ] #gravity in m/s^2
51
52#rear wheel R:
53rR = 0.3 #rear wheel radius
54mR = 2 #rear wheel mass
55IRxx = 0.0603 #rear wheel inertia xx = zz ; but in EXUDYN, x=rotation axis!
56IRyy = 0.12 #rear wheel inertia yy (around wheel axis)
57inertiaR = RigidBodyInertia(mass=mR, inertiaTensor=np.array([[IRyy,0,0],[0,IRxx,0],[0,0,IRxx]]))
58
59#front wheel F:
60rF = 0.35 #rear wheel radius
61mF = 3 #rear wheel mass
62IFxx = 0.1405 #rear wheel inertia xx = zz ; but in EXUDYN, x=rotation axis!
63IFyy = 0.28 #rear wheel inertia yy (around wheel axis)
64inertiaF = RigidBodyInertia(mass=mF, inertiaTensor=np.array([[IFyy,0,0],[0,IFxx,0],[0,0,IFxx]]))
65
66#rear body B:
67xB = 0.3 #COM x
68zB = -0.9*sZ #COM z
69bCOM = np.array([xB, 0, zB])
70mB = 85 #rear body (and driver) mass
71zz=-1
72inertiaB = RigidBodyInertia(mass=mB,
73 inertiaTensor=np.array([[9.2,0,2.4*zz],[0,11,0],[2.4*zz,0,2.8]]),
74 # inertiaTensor=np.diag([1,1,1]),
75 com=np.zeros(3)) #reference position = COM for this body
76 # com=bCOM)
77
78#front Handlebar H:
79xH = 0.9 #COM x
80zH = -0.7*sZ #COM z
81hCOM = np.array([xH, 0, zH])
82mH = 4 #handle bar mass
83inertiaH = RigidBodyInertia(mass=mH,
84 inertiaTensor=np.array([[0.05892, 0, -0.00756*zz],[0,0.06,0],[-0.00756*zz, 0, 0.00708]]),
85 # inertiaTensor=np.diag([1,1,1]),
86 com=np.zeros(3)) #reference position = COM for this body
87 # com=hCOM)
88
89#geometrical parameters for joints
90P1 = np.array([0,0,-0.3*sZ])
91P2 = np.array([0.82188470506, 0, -0.85595086466*sZ])
92P3 = np.array([w, 0, -0.35*sZ])
93
94
95#stable velocity limits according to linear theory:
96vMin = 4.29238253634111
97vMax = 6.02426201538837
98
99maneuver = 'M1'
100if maneuver == 'M1':
101 vX0 = 4. #initial forward velocity in x-direction
102 omegaX0 = 0.05 #initial roll velocity around x axis
103elif maneuver == 'M2':
104 vX0 = 4.6 #initial forward velocity in x-direction
105 omegaX0 = 0.5 #initial roll velocity around x axis
106elif maneuver == 'M3':
107 vX0 = 8 #initial forward velocity in x-direction
108 omegaX0 = 0.05 #initial roll velocity around x axis
109
110omegaRy0 = -vX0/rR*sZ #initial angular velocity of rear wheel
111omegaFy0 = -vX0/rF*sZ #initial angular velocity of front wheel
112
113#%%++++++++++++++++++++++++++++++++++++++++++++++++
114#visualization:
115dY = 0.02
116#graphicsFrame = graphics.Brick(centerPoint=[0,0,0],size=[dFoot*1.1,0.7*rFoot,0.7*rFoot], color=graphics.color.lightred)
117graphicsR = graphics.Cylinder(pAxis=[-1*dY,0,0], vAxis=[dY*2,0,0], radius=rR, color=graphics.color.steelblue, nTiles=4)
118graphicsF = graphics.Cylinder(pAxis=[-1*dY,0,0], vAxis=[dY*2,0,0], radius=rF, color=graphics.color.steelblue, nTiles=4)
119graphicsB = graphics.Cylinder(pAxis=P1-bCOM, vAxis=P2-P1, radius=dY*1.5, color=graphics.color.lightred)
120graphicsB2 = graphics.Sphere(point=[0,0,0], radius=3*dY, color=graphics.color.lightgrey)
121graphicsH = graphics.Cylinder(pAxis=P3-hCOM, vAxis=P2-P3, radius=dY*1.3, color=graphics.color.lightgreen)
122
123#option to track motion of bicycle
124#old version; does not work any more to set state during graphics user function => use trackMarker
125if False:
126 #add user function to track bicycle frame
127 def UFgraphics(mbs, objectNum):
128 n = mbs.variables['nTrackNode']
129 p = mbs.GetNodeOutput(n,exu.OutputVariableType.Position,
130 configuration=exu.ConfigurationType.Visualization)
131 rs=SC.renderer.GetState()
132 A = np.array(rs['modelRotation'])
133 p = A.T @ p
134 rs['centerPoint']=[p[0],p[1],p[2]]
135 SC.renderer.SetState(rs)
136 return []
137
138 #add object with graphics user function
139 oGround2 = mbs.AddObject(ObjectGround(visualization=
140 VObjectGround(graphicsDataUserFunction=UFgraphics)))
141#add rigid bodies
142#rear wheel
143resultR = mbs.CreateRigidBody(
144 referencePosition = P1,
145 referenceRotationMatrix = RotationMatrixZ(np.pi*0.5),
146 initialVelocity = [vX0, omegaX0*P1[2]*sZ, 0],
147 initialAngularVelocity = [omegaX0, omegaRy0, 0],
148 inertia = inertiaR,
149 gravity = g,
150 graphicsDataList = [graphicsR],
151 returnDict = True)
152
153nR, bR = resultR['nodeNumber'], resultR['bodyNumber']
154
155mbs.variables['nTrackNode'] = nR #node to be tracked
156
157#front wheel
158resultF = mbs.CreateRigidBody(
159 referencePosition = P3,
160 referenceRotationMatrix = RotationMatrixZ(pi*0.5),
161 initialVelocity = [vX0, omegaX0*P3[2]*sZ, 0],
162 initialAngularVelocity = [omegaX0, omegaFy0, 0],
163 inertia = inertiaF,
164 gravity = g,
165 nodeType = exu.NodeType.RotationEulerParameters,
166 graphicsDataList = [graphicsF],
167 returnDict = True
168)
169nF, bF = resultF['nodeNumber'], resultF['bodyNumber']
170
171#read body
172resultB = mbs.CreateRigidBody(
173 referencePosition = bCOM,
174 initialVelocity = [vX0, omegaX0*bCOM[2]*sZ, 0],
175 initialAngularVelocity = [omegaX0, 0, 0],
176 inertia = inertiaB,
177 gravity = g,
178 nodeType = exu.NodeType.RotationEulerParameters,
179 graphicsDataList = [graphicsB, graphicsB2],
180 returnDict = True
181)
182nB, bB = resultB['nodeNumber'], resultB['bodyNumber']
183
184#handle
185resultH = mbs.CreateRigidBody(
186 referencePosition = hCOM,
187 initialVelocity = [vX0, omegaX0*hCOM[2]*sZ, 0],
188 initialAngularVelocity = [omegaX0, 0, 0],
189 inertia = inertiaH,
190 gravity = g,
191 nodeType = exu.NodeType.RotationEulerParameters,
192 graphicsDataList = [graphicsH],
193 returnDict = True
194)
195nH, bH = resultH['nodeNumber'], resultH['bodyNumber']
196
197
198#%%++++++++++++++++++++++++++++++++++++++++++++++++
199#ground body and marker
200gGround = graphics.CheckerBoard(point=[0,0,0], size=50, nTiles=64)
201oGround = mbs.CreateGround(graphicsDataList=[gGround])
202markerGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
203
204markerR = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bR, localPosition=[0,0,0]))
205markerF = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bF, localPosition=[0,0,0]))
206markerB1 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bB, localPosition=P1-bCOM))
207
208sMarkerR = mbs.AddSensor(SensorMarker(markerNumber=markerR, outputVariableType=exu.OutputVariableType.Position))
209sMarkerB1= mbs.AddSensor(SensorMarker(markerNumber=markerB1,outputVariableType=exu.OutputVariableType.Position))
210
211#%%++++++++++++++++++++++++++++++++++++++++++++++++
212#add joints:
213useJoints = True
214if useJoints:
215 oJointRW = mbs.CreateRevoluteJoint(itemNumbers=[bR, bB], position=P1, axis=[0,1,0],
216 axisRadius=0.5*dY, axisLength=5*dY)
217 oJointFW = mbs.CreateRevoluteJoint(itemNumbers=[bF, bH], position=P3, axis=[0,1,0],
218 axisRadius=0.5*dY, axisLength=5*dY)
219 oJointSteer = mbs.CreateRevoluteJoint(itemNumbers=[bB, bH],
220 position=P2-bCOM, useGlobalFrame=False,
221 axis=RotationMatrixY(-lam) @ [0,0,1],
222 axisRadius=0.5*dY, axisLength=5*dY)
223#%%++++++++++++++++++++++++++++++++++++++++++++++++
224#add 'rolling disc' for wheels:
225cStiffness = 5e4*10 #spring stiffness: 50N==>F/k = u = 0.001m (penetration)
226cDamping = cStiffness*0.05*0.1 #think on a one-mass spring damper
227nGenericR = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
228if False:
229 oRollingR=mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[markerGround, markerR],
230 nodeNumber = nGenericR,
231 discRadius=rR,
232 planeNormal=[0,0,1],
233 dryFriction=[0.8,0.8],
234 dryFrictionProportionalZone=1e-2,
235 rollingViscousFriction=0.,
236 contactStiffness=cStiffness,
237 contactDamping=cDamping,
238 #activeConnector = False, #set to false to deactivated
239 visualization=VObjectConnectorRollingDiscPenalty(show=True,
240 discWidth=dY, color=graphics.color.blue)))
241
242 nGenericF = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
243 oRollingF=mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[markerGround, markerF],
244 nodeNumber = nGenericF,
245 discRadius=rF,
246 planeNormal=[0,0,1],
247 dryFriction=[0.8,0.8],
248 dryFrictionProportionalZone=1e-2,
249 rollingViscousFriction=0.,
250 contactStiffness=cStiffness,
251 contactDamping=cDamping,
252 #activeConnector = False, #set to false to deactivated
253 visualization=VObjectConnectorRollingDiscPenalty(show=True, discWidth=dY, color=graphics.color.blue)))
254else:
255 if True:
256 oRollingR=mbs.AddObject(ObjectJointRollingDisc(markerNumbers=[markerGround, markerR],
257 discRadius=rR,
258 visualization=VObjectJointRollingDisc(show=True, discWidth=dY, color=graphics.color.blue)))
259
260 oRollingF=mbs.AddObject(ObjectJointRollingDisc(markerNumbers=[markerGround, markerF],
261 discRadius=rF,
262 visualization=VObjectJointRollingDisc(show=True, discWidth=dY, color=graphics.color.blue)))
263
264
265
266#%%++++++++++++++++++++++++++++++++++++++++++++++++
267#add sensors
268addSensors = True
269if addSensors:
270 sForwardVel = mbs.AddSensor(SensorBody(bodyNumber=bB, fileName='solution/bicycleBvelLocal.txt',
271 localPosition=P1-bCOM,
272 outputVariableType=exu.OutputVariableType.VelocityLocal))
273
274 sBAngVelLocal = mbs.AddSensor(SensorBody(bodyNumber=bB, fileName='solution/bicycleBangVelLocal.txt',
275 outputVariableType=exu.OutputVariableType.AngularVelocityLocal))
276 sBrot = mbs.AddSensor(SensorBody(bodyNumber=bB, fileName='solution/bicycleBrot.txt',
277 outputVariableType=exu.OutputVariableType.Rotation))
278
279
280 bodies = [bB, bH, bR, bF]
281 massBodies = [mB, mH, mR, mF]
282 inertiaBodies = [inertiaB.inertiaTensor,
283 inertiaH.inertiaTensor,
284 inertiaR.inertiaTensor,
285 inertiaF.inertiaTensor]
286
287 nBodies = len(bodies)
288 sList = []
289 for b in bodies:
290 sPosCOM = mbs.AddSensor(SensorBody(bodyNumber=b, writeToFile=False,
291 outputVariableType=exu.OutputVariableType.Position))
292 sVelCOM = mbs.AddSensor(SensorBody(bodyNumber=b, writeToFile=False,
293 outputVariableType=exu.OutputVariableType.Velocity))
294 sAngVelLocal = mbs.AddSensor(SensorBody(bodyNumber=b, writeToFile=False,
295 outputVariableType=exu.OutputVariableType.AngularVelocityLocal))
296
297 sList += [sPosCOM,sVelCOM,sAngVelLocal]
298
299 if useJoints:
300 sSteerAngle = mbs.AddSensor(SensorObject(objectNumber=oJointSteer, fileName='solution/bicycleSteerAngle.txt',
301 outputVariableType=exu.OutputVariableType.Rotation))
302 sSteerVel = mbs.AddSensor(SensorObject(objectNumber=oJointSteer, fileName='solution/bicycleSteerVelocity.txt',
303 outputVariableType=exu.OutputVariableType.AngularVelocityLocal))
304
305
306 #create user joint for kinetic and potential energy
307 def UFsensorEnergy(mbs, t, sensorNumbers, factors, configuration):
308 E = 0
309 P = 0
310 for i in range(nBodies):
311 pos = mbs.GetSensorValues(sensorNumbers[i*3+0])
312 vel = mbs.GetSensorValues(sensorNumbers[i*3+1]) #vel
313 omega = mbs.GetSensorValues(sensorNumbers[i*3+2]) #ang vel local
314 E += 0.5 * np.linalg.norm(vel)**2 * massBodies[i]
315 E += 0.5 * np.array(omega) @ inertiaBodies[i] @ omega
316
317 P -= np.dot(g,pos)*massBodies[i]
318 return [P, E, E+P] #return potential, kinetic and total mechanical energy
319
320 sEnergy = mbs.AddSensor(SensorUserFunction(sensorNumbers=sList, #factors=[180/pi],
321 fileName='solution/sensorKineticPotentialEnergy.txt',
322 sensorUserFunction=UFsensorEnergy))
323
324 def UFsensorResults(mbs, t, sensorNumbers, factors, configuration):
325 energy = mbs.GetSensorValues(sensorNumbers[0])
326 forwardVel = mbs.GetSensorValues(sensorNumbers[1])
327 angVelLocalB = mbs.GetSensorValues(sensorNumbers[2])
328 rotB = mbs.GetSensorValues(sensorNumbers[3])
329 steerAngle = mbs.GetSensorValues(sensorNumbers[4])
330 steerVel = mbs.GetSensorValues(sensorNumbers[5])
331 return [rotB[0], angVelLocalB[0], forwardVel[0], energy[0], energy[1], energy[2], -steerAngle[2], -steerVel[2]] #return kinetic, potential and total mechanical energy
332
333 # 1=roll angle, 2=roll angular velocity, 3=forward speed, 4=potential energy, 5=kinetic energy, 6=mechanical energy, 7=steer angle, and 8=steer velocity
334 sResults = mbs.AddSensor(SensorUserFunction(sensorNumbers=[sEnergy,sForwardVel,sBAngVelLocal,sBrot, sSteerAngle, sSteerVel],
335 fileName='solution/sensorResults'+maneuver+'.txt',
336 sensorUserFunction=UFsensorResults))
337
338#%%++++++++++++++++++++++++++++++++++++++++++++++++
339#simulate:
340mbs.Assemble()
341
342pR = mbs.GetSensorValues(sMarkerR)
343pB1 = mbs.GetSensorValues(sMarkerB1)
344print("pR=",pR)
345print("pB1=",pB1)
346simulationSettings = exu.SimulationSettings() #takes currently set values or default values
347
348tEnd = 5
349h=0.001 #use small step size to detext contact switching
350
351simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
352simulationSettings.timeIntegration.endTime = tEnd
353simulationSettings.solution.file.write= False #set False for CPU performance measurement
354simulationSettings.solution.sensors.writePeriod = 0.01
355simulationSettings.solution.precision = 16
356
357simulationSettings.timeIntegration.verboseMode = 1
358#simulationSettings.linearSolver.ignoreSingularJacobian = True
359
360# simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
361# simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
362simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.7
363simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
364simulationSettings.timeIntegration.newton.useModifiedNewton = True
365
366#SC.visualizationSettings.general.useMultiThreadedRendering = False
367SC.visualizationSettings.nodes.show = True
368SC.visualizationSettings.nodes.drawNodesAsPoint = False
369SC.visualizationSettings.nodes.showBasis = True
370SC.visualizationSettings.nodes.basisSize = 0.015
371
372SC.visualizationSettings.view0.camera.trackMarker = markerR #rear wheel
373
374if False: #record animation frames:
375 SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
376 SC.visualizationSettings.view0.window.renderWindowSize=[1600,1024]
377 SC.visualizationSettings.openGL.multiSampling = 4
378 simulationSettings.solution.recordImagesInterval = 0.02
379
380SC.visualizationSettings.general.autoFitScene = False #use loaded render state
381useGraphics = True
382if useGraphics:
383 SC.renderer.Start()
384 SC.renderer.RestoreSavedState()
385 SC.renderer.DoIdleTasks()
386
387mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2)
388#mbs.SolveDynamic(simulationSettings, showHints=True)
389
390
391#%%+++++++++++++++++++++++++++++
392if useGraphics:
393 SC.renderer.DoIdleTasks()
394 SC.renderer.Stop() #safely close rendering window!
395
396#%%++++++++++++++++++++++++++++++++++++++++++++++q+++++++
397if addSensors:
398 #plot results
399
400
401
402 import matplotlib.pyplot as plt
403 import matplotlib.ticker as ticker
404 plt.close('all')
405
406
407 # mbs.PlotSensor(sensorNumbers=[sBpos,sBpos,sBpos], components=[0,1,2])
408 #plt.figure('lateral position')
409 #mbs.PlotSensor(sensorNumbers=[sBpos], components=[1])
410
411 plt.figure('forward velocity')
412 mbs.PlotSensor(sensorNumbers=[sForwardVel], components=[0])
413 # mbs.PlotSensor(sensorNumbers=[sBvelLocal,sBvelLocal,sBvelLocal], components=[0,1,2])
414
415 # plt.figure('local ang velocities')
416 # mbs.PlotSensor(sensorNumbers=[sBAngVelLocal,sBAngVelLocal,sBAngVelLocal], components=[0,1,2])
417 # if False:
418 # import matplotlib.pyplot as plt
419 # import matplotlib.ticker as ticker
420
421 # 1=roll angle, 2=roll angular velocity, 3=forward speed, 4=potential energy, 5=kinetic energy, 6=mechanical energy, 7=steer angle, and 8=steer velocity
422 data = np.loadtxt(OutputFilePath('solution/uncontrolledBicycleGonzalez.txt'))#, comments='#', delimiter='')
423 plt.plot(data[:,0], data[:,9], 'b:',label='')
424
425 data2 = np.loadtxt(OutputFilePath('solution/uncontrolledBicycleSanjurjo.txt'))#, comments='#', delimiter='')
426 plt.plot(data2[:,0], data2[:,3+8], 'g:',label='')
427
428 plt.figure('steer vel')
429 mbs.PlotSensor(sensorNumbers=[sSteerVel], components=[2])
430 plt.plot(data2[:,0], -data2[:,8+8], 'g:',label='')
431
432 plt.figure('steer ang')
433 mbs.PlotSensor(sensorNumbers=[sSteerAngle], components=[2])
434 plt.plot(data2[:,0], -data2[:,7+8], 'g:',label='')
435
436 plt.figure('roll')
437 mbs.PlotSensor(sensorNumbers=[sBrot], components=[0])
438 plt.plot(data2[:,0], data2[:,1+8], 'g:',label='')
439
440 plt.figure('roll ang vel')
441 mbs.PlotSensor(sensorNumbers=[sBAngVelLocal], components=[0])
442 plt.plot(data2[:,0], data2[:,2+8], 'g:',label='')
443
444
445 plt.figure('potential energy')
446 mbs.PlotSensor(sensorNumbers=[sEnergy], components=[0])
447 plt.plot(data2[:,0], data2[:,4+8], 'g:',label='')
448
449 plt.figure('kinetic energy')
450 mbs.PlotSensor(sensorNumbers=[sEnergy], components=[1])
451 plt.plot(data2[:,0], data2[:,5+8], 'g:',label='')
452
453 plt.figure('total energy')
454 mbs.PlotSensor(sensorNumbers=[sEnergy], components=[2])
455
456 dataE = np.loadtxt(OutputFilePath('solution/sensorKineticPotentialEnergy.txt'), comments='#', delimiter=',')
457 performance = 100*(max(dataE[:,3]) - min(dataE[:,3])) / dataE[0,3]
458 print("performance = ", performance, "(must by < 1e-3)")
459
460 #CPU performance with 20 seconds simulation time, maneuver 2
461 #performance = 0.000915 < 0.001: max h=0.014; CPU time = 0.596 seconds on Intel Core i9
462 #reference solution computed with:
463 # performance = 6.423e-06: max h=0.001; CPU time = 5.5935 seconds on Intel Core i9
464
465
466#%%+++++++++++++++++
467#merge result files for IFToMM
468if True:
469 dataM1 = np.loadtxt(OutputFilePath('solution/sensorResultsM1.txt'), comments='#', delimiter=',')
470 dataM2 = np.loadtxt(OutputFilePath('solution/sensorResultsM2.txt'), comments='#', delimiter=',')
471 dataM3 = np.loadtxt(OutputFilePath('solution/sensorResultsM3.txt'), comments='#', delimiter=',')
472
473 data = np.hstack((dataM1,dataM2[:,1:],dataM3[:,1:]))
474 np.savetxt('solution/bicycleResultsIFToMM.txt', data, fmt='%1.15e')
475
476# benchmark results:
477# 6.423e-06
478# 5.5935
479# Intel(R) Core(TM) i9-7940X CPU @ 3.10GHz
480# Simulated using C++/Python library EXUDYN, see https://github.com/jgerstmayr/EXUDYN.
481# Solved using implicit trapezoidal rule (energy conserving) with Index 2 constraint reduction, using redundant coordinate formulation. Rigid bodies are modeled with Euler parameters.
482# With a larger step size of 0.014 we still obtain a performance <0.001, but have CPU time of 0.596 seconds.