Forgot the first part of the output _________________________________ test_bicycle _________________________________
@slow def test_bicycle(): # Code to get equations of motion for a bicycle modeled as in: # J.P Meijaard, Jim M Papadopoulos, Andy Ruina and A.L Schwab. Linearized # dynamics equations for the balance and steer of a bicycle: a benchmark # and review. Proceedings of The Royal Society (2007) 463, 1955-1982 # doi: 10.1098/rspa.2007.1857 # Note that this code has been crudely ported from Autolev, which is the # reason for some of the unusual naming conventions. It was purposefully as # similar as possible in order to aide debugging. # Declare Coordinates & Speeds # Simple definitions for qdots - qd = u # Speeds are: yaw frame ang. rate, roll frame ang. rate, rear wheel frame # ang. rate (spinning motion), frame ang. rate (pitching motion), steering # frame ang. rate, and front wheel ang. rate (spinning motion). # Wheel positions are ignorable coordinates, so they are not introduced. q1, q2, q4, q5 = dynamicsymbols('q1 q2 q4 q5') q1d, q2d, q4d, q5d = dynamicsymbols('q1 q2 q4 q5', 1) u1, u2, u3, u4, u5, u6 = dynamicsymbols('u1 u2 u3 u4 u5 u6') u1d, u2d, u3d, u4d, u5d, u6d = dynamicsymbols('u1 u2 u3 u4 u5 u6', 1) # Declare System's Parameters WFrad, WRrad, htangle, forkoffset = symbols('WFrad WRrad htangle forkoffset') forklength, framelength, forkcg1 = symbols('forklength framelength forkcg1') forkcg3, framecg1, framecg3, Iwr11 = symbols('forkcg3 framecg1 framecg3 Iwr11') Iwr22, Iwf11, Iwf22, Iframe11 = symbols('Iwr22 Iwf11 Iwf22 Iframe11') Iframe22, Iframe33, Iframe31, Ifork11 = symbols('Iframe22 Iframe33 Iframe31 Ifork11') Ifork22, Ifork33, Ifork31, g = symbols('Ifork22 Ifork33 Ifork31 g') mframe, mfork, mwf, mwr = symbols('mframe mfork mwf mwr') # Set up reference frames for the system # N - inertial # Y - yaw # R - roll # WR - rear wheel, rotation angle is ignorable coordinate so not oriented # Frame - bicycle frame # TempFrame - statically rotated frame for easier reference inertia definition # Fork - bicycle fork # TempFork - statically rotated frame for easier reference inertia definition # WF - front wheel, again posses a ignorable coordinate N = ReferenceFrame('N') Y = N.orientnew('Y', 'Axis', [q1, N.z]) R = Y.orientnew('R', 'Axis', [q2, Y.x]) Frame = R.orientnew('Frame', 'Axis', [q4 + htangle, R.y]) WR = ReferenceFrame('WR') TempFrame = Frame.orientnew('TempFrame', 'Axis', [-htangle, Frame.y]) Fork = Frame.orientnew('Fork', 'Axis', [q5, Frame.x]) TempFork = Fork.orientnew('TempFork', 'Axis', [-htangle, Fork.y]) WF = ReferenceFrame('WF') # Kinematics of the Bicycle First block of code is forming the positions of # the relevant points # rear wheel contact -> rear wheel mass center -> frame mass center + # frame/fork connection -> fork mass center + front wheel mass center -> # front wheel contact point WR_cont = Point('WR_cont') WR_mc = WR_cont.locatenew('WR_mc', WRrad * R.z) Steer = WR_mc.locatenew('Steer', framelength * Frame.z) Frame_mc = WR_mc.locatenew('Frame_mc', - framecg1 * Frame.x + framecg3 * Frame.z) Fork_mc = Steer.locatenew('Fork_mc', - forkcg1 * Fork.x + forkcg3 * Fork.z) WF_mc = Steer.locatenew('WF_mc', forklength * Fork.x + forkoffset * Fork.z) WF_cont = WF_mc.locatenew('WF_cont', WFrad * (dot(Fork.y, Y.z) * Fork.y - Y.z).normalize()) # Set the angular velocity of each frame. # Angular accelerations end up being calculated automatically by # differentiating the angular velocities when first needed. # u1 is yaw rate # u2 is roll rate # u3 is rear wheel rate # u4 is frame pitch rate # u5 is fork steer rate # u6 is front wheel rate Y.set_ang_vel(N, u1 * Y.z) R.set_ang_vel(Y, u2 * R.x) WR.set_ang_vel(Frame, u3 * Frame.y) Frame.set_ang_vel(R, u4 * Frame.y) Fork.set_ang_vel(Frame, u5 * Fork.x) WF.set_ang_vel(Fork, u6 * Fork.y) # Form the velocities of the previously defined points, using the 2 - point # theorem (written out by hand here). Accelerations again are calculated # automatically when first needed. WR_cont.set_vel(N, 0) WR_mc.v2pt_theory(WR_cont, N, WR) Steer.v2pt_theory(WR_mc, N, Frame) Frame_mc.v2pt_theory(WR_mc, N, Frame) Fork_mc.v2pt_theory(Steer, N, Fork) WF_mc.v2pt_theory(Steer, N, Fork) WF_cont.v2pt_theory(WF_mc, N, WF) # Sets the inertias of each body. Uses the inertia frame to construct the # inertia dyadics. Wheel inertias are only defined by principle moments of # inertia, and are in fact constant in the frame and fork reference frames; # it is for this reason that the orientations of the wheels does not need # to be defined. The frame and fork inertias are defined in the 'Temp' # frames which are fixed to the appropriate body frames; this is to allow # easier input of the reference values of the benchmark paper. Note that # due to slightly different orientations, the products of inertia need to # have their signs flipped; this is done later when entering the numerical # value. Frame_I = (inertia(TempFrame, Iframe11, Iframe22, Iframe33, 0, 0, Iframe31), Frame_mc) Fork_I = (inertia(TempFork, Ifork11, Ifork22, Ifork33, 0, 0, Ifork31), Fork_mc) WR_I = (inertia(Frame, Iwr11, Iwr22, Iwr11), WR_mc) WF_I = (inertia(Fork, Iwf11, Iwf22, Iwf11), WF_mc) # Declaration of the RigidBody containers. :: BodyFrame = RigidBody('BodyFrame', Frame_mc, Frame, mframe, Frame_I) BodyFork = RigidBody('BodyFork', Fork_mc, Fork, mfork, Fork_I) BodyWR = RigidBody('BodyWR', WR_mc, WR, mwr, WR_I) BodyWF = RigidBody('BodyWF', WF_mc, WF, mwf, WF_I) # The kinematic differential equations; they are defined quite simply. Each # entry in this list is equal to zero. kd = [q1d - u1, q2d - u2, q4d - u4, q5d - u5] # The nonholonomic constraints are the velocity of the front wheel contact # point dotted into the X, Y, and Z directions; the yaw frame is used as it # is "closer" to the front wheel (1 less DCM connecting them). These # constraints force the velocity of the front wheel contact point to be 0 # in the inertial frame; the X and Y direction constraints enforce a # "no-slip" condition, and the Z direction constraint forces the front # wheel contact point to not move away from the ground frame, essentially # replicating the holonomic constraint which does not allow the frame pitch # to change in an invalid fashion. conlist_speed = [WF_cont.vel(N) & Y.x, WF_cont.vel(N) & Y.y, WF_cont.vel(N) & Y.z] # The holonomic constraint is that the position from the rear wheel contact # point to the front wheel contact point when dotted into the # normal-to-ground plane direction must be zero; effectively that the front # and rear wheel contact points are always touching the ground plane. This # is actually not part of the dynamic equations, but instead is necessary # for the lineraization process. conlist_coord = [WF_cont.pos_from(WR_cont) & Y.z] # The force list; each body has the appropriate gravitational force applied # at its mass center. FL = [(Frame_mc, -mframe * g * Y.z), (Fork_mc, -mfork * g * Y.z), (WF_mc, -mwf * g * Y.z), (WR_mc, -mwr * g * Y.z)] BL = [BodyFrame, BodyFork, BodyWR, BodyWF] # The N frame is the inertial frame, coordinates are supplied in the order # of independent, dependent coordinates, as are the speeds. The kinematic # differential equation are also entered here. Here the dependent speeds # are specified, in the same order they were provided in earlier, along # with the non-holonomic constraints. The dependent coordinate is also # provided, with the holonomic constraint. Again, this is only provided # for the linearization process. KM = KanesMethod(N, q_ind=[q1, q2, q5], q_dependent=[q4], configuration_constraints=conlist_coord, u_ind=[u2, u3, u5], u_dependent=[u1, u4, u6], velocity_constraints=conlist_speed, kd_eqs=kd) (fr, frstar) = KM.kanes_equations(FL, BL) # This is the start of entering in the numerical values from the benchmark # paper to validate the eigen values of the linearized equations from this # model to the reference eigen values. Look at the aforementioned paper for # more information. Some of these are intermediate values, used to # transform values from the paper into the coordinate systems used in this # model. PaperRadRear = 0.3 PaperRadFront = 0.35 HTA = evalf.N(pi / 2 - pi / 10) TrailPaper = 0.08 rake = evalf.N(-(TrailPaper*sin(HTA)-(PaperRadFront*cos(HTA)))) PaperWb = 1.02 PaperFrameCgX = 0.3 PaperFrameCgZ = 0.9 PaperForkCgX = 0.9 PaperForkCgZ = 0.7 FrameLength = evalf.N(PaperWb*sin(HTA)-(rake-(PaperRadFront-PaperRadRear)*cos(HTA))) FrameCGNorm = evalf.N((PaperFrameCgZ - PaperRadRear-(PaperFrameCgX/sin(HTA))*cos(HTA))*sin(HTA)) FrameCGPar = evalf.N((PaperFrameCgX / sin(HTA) + (PaperFrameCgZ - PaperRadRear - PaperFrameCgX / sin(HTA) * cos(HTA)) * cos(HTA))) tempa = evalf.N((PaperForkCgZ - PaperRadFront)) tempb = evalf.N((PaperWb-PaperForkCgX)) tempc = evalf.N(sqrt(tempa**2+tempb**2)) PaperForkL = evalf.N((PaperWb*cos(HTA)-(PaperRadFront-PaperRadRear)*sin(HTA))) ForkCGNorm = evalf.N(rake+(tempc * sin(pi/2-HTA-acos(tempa/tempc)))) ForkCGPar = evalf.N(tempc * cos((pi/2-HTA)-acos(tempa/tempc))-PaperForkL) # Here is the final assembly of the numerical values. The symbol 'v' is the # forward speed of the bicycle (a concept which only makes sense in the # upright, static equilibrium case?). These are in a dictionary which will # later be substituted in. Again the sign on the *product* of inertia # values is flipped here, due to different orientations of coordinate # systems. v = symbols('v') val_dict = {WFrad: PaperRadFront, WRrad: PaperRadRear, htangle: HTA, forkoffset: rake, forklength: PaperForkL, framelength: FrameLength, forkcg1: ForkCGPar, forkcg3: ForkCGNorm, framecg1: FrameCGNorm, framecg3: FrameCGPar, Iwr11: 0.0603, Iwr22: 0.12, Iwf11: 0.1405, Iwf22: 0.28, Ifork11: 0.05892, Ifork22: 0.06, Ifork33: 0.00708, Ifork31: 0.00756, Iframe11: 9.2, Iframe22: 11, Iframe33: 2.8, Iframe31: -2.4, mfork: 4, mframe: 85, mwf: 3, mwr: 2, g: 9.81, q1: 0, q2: 0, q4: 0, q5: 0, u1: 0, u2: 0, u3: v / PaperRadRear, u4: 0, u5: 0, u6: v / PaperRadFront} # Linearizes the forcing vector; the equations are set up as MM udot = # forcing, where MM is the mass matrix, udot is the vector representing the # time derivatives of the generalized speeds, and forcing is a vector which # contains both external forcing terms and internal forcing terms, such as # centripital or coriolis forces. This actually returns a matrix with as # many rows as *total* coordinates and speeds, but only as many columns as # independent coordinates and speeds. > forcing_lin = KM.linearize()[0] On Fri, Sep 26, 2014 at 9:34 AM, Peter Brady <petertbr...@gmail.com> wrote: > No. But there are actually a few tests which pass under our testing > framework (bin/test) but not under py.test. Here's the output > > sympy/physics/mechanics/tests/test_kane3.py:253: > _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ > _ _ _ > > self = <sympy.physics.mechanics.kane.KanesMethod object at 0x7f921ae68208> > kwargs = {} > > def linearize(self, **kwargs): > """ Linearize the equations of motion about a symbolic operating > point. > > If kwarg A_and_B is False (default), returns M, A, B, r for the > linearized form, M*[q', u']^T = A*[q_ind, u_ind]^T + B*r. > > If kwarg A_and_B is True, returns A, B, r for the linearized > form > dx = A*x + B*r, where x = [q_ind, u_ind]^T. Note that this is > computationally intensive if there are many symbolic > parameters. For > this reason, it may be more desirable to use the default > A_and_B=False, > returning M, A, and B. Values may then be substituted in to > these > matrices, and the state space form found as > A = P.T*M.inv()*A, B = P.T*M.inv()*B, where P = > Linearizer.perm_mat. > > In both cases, r is found as all dynamicsymbols in the > equations of > motion that are not part of q, u, q', or u'. They are sorted in > canonical form. > > The operating points may be also entered using the > ``op_point`` kwarg. > This takes a dictionary of {symbol: value}, or a an iterable > of such > dictionaries. The values may be numberic or symbolic. The more > values > you can specify beforehand, the faster this computation will > run. > > As part of the deprecation cycle, the new method will not be > used unless > the kwarg ``new_method`` is set to True. If the kwarg is > missing, or set > to false, the old linearization method will be used. After > next release > the need for this kwarg will be removed. > > For more documentation, please see the ``Linearizer`` class.""" > > if 'new_method' not in kwargs or not kwargs['new_method']: > # User is still using old code. > SymPyDeprecationWarning('The linearize class method has > changed ' > 'to a new interface, the old method is deprecated. To ' > 'use the new method, set the kwarg `new_method=True`. ' > 'For more information, read the docstring ' > 'of `linearize`.').warn() > > return self._old_linearize() > > sympy/physics/mechanics/kane.py:494: > _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ > _ _ _ > > self = <sympy.physics.mechanics.kane.KanesMethod object at 0x7f921ae68208> > > def _old_linearize(self): > """Old method to linearize the equations of motion. Returns a > tuple of > (f_lin_A, f_lin_B, y) for forming [M]qudot = [f_lin_A]qu + > [f_lin_B]y. > > Deprecated in favor of new method using Linearizer class. > Please change > your code to use the new `linearize` method.""" > > if (self._fr is None) or (self._frstar is None): > raise ValueError('Need to compute Fr, Fr* first.') > > # Note that this is now unneccessary, and it should never be > # encountered; I still think it should be in here in case the user > # manually sets these matrices incorrectly. > for i in self.q: > if self._k_kqdot.diff(i) != 0 * self._k_kqdot: > raise ValueError('Matrix K_kqdot must not depend on any > q.') > > t = dynamicsymbols._t > uaux = self._uaux > uauxdot = [diff(i, t) for i in uaux] > # dictionary of auxiliary speeds & derivatives which are equal to > zero > > subdict = dict(list(zip(uaux + uauxdot, [0] * (len(uaux) + > len(uauxdot))))) > > sympy/physics/mechanics/kane.py:522: > _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ > _ _ _ > > self = Matrix(0, 0, []), other = [] > > @wraps(func) > def binary_op_wrapper(self, other): > if hasattr(other, '_op_priority'): > if other._op_priority > self._op_priority: > try: > f = getattr(other, method_name) > except AttributeError: > pass > else: > return f(self) > > return func(self, other) > > sympy/core/decorators.py:118: > _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ > _ _ _ > > self = Matrix(0, 0, []), other = [] > > @call_highest_priority('__radd__') > def __add__(self, other): > > return super(DenseMatrix, self).__add__(_force_mutable(other)) > > sympy/matrices/dense.py:531: > _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ > _ _ _ > > self = Matrix(0, 0, []), other = [] > > def __add__(self, other): > """Return self + other, raising ShapeError if shapes don't > match.""" > if getattr(other, 'is_Matrix', False): > A = self > B = other > if A.shape != B.shape: > raise ShapeError("Matrix size mismatch.") > alst = A.tolist() > blst = B.tolist() > ret = [S.Zero]*A.rows > for i in range(A.shape[0]): > ret[i] = list(map(lambda j, k: j + k, alst[i], blst[i])) > rv = classof(A, B)._new(ret) > if 0 in A.shape: > rv = rv.reshape(*A.shape) > return rv > > raise TypeError('cannot add matrix and %s' % type(other)) > E TypeError: cannot add matrix and <class 'list'> > > sympy/matrices/matrices.py:579: TypeError > ------------------------------- Captured stderr > -------------------------------- > /home/ptb/gitrepos/sympy/sympy/physics/mechanics/kane.py:493: > SymPyDeprecationWarning: > > The linearize class method has changed to a new interface, the old > method is deprecated. To use the new method, set the kwarg > `new_method=True`. For more information, read the docstring of > `linearize`. > > 'of `linearize`.').warn() > > > On Fri, Sep 26, 2014 at 9:28 AM, Jason Moore <moorepa...@gmail.com> wrote: > >> BTW, Peter, did that kane3 test pass on your machine? >> >> >> Jason >> moorepants.info >> +01 530-601-9791 >> >> On Fri, Sep 26, 2014 at 11:27 AM, Jason Moore <moorepa...@gmail.com> >> wrote: >> >>> It is certainly too slow to be a useful test. But that example should be >>> available to run at least before any SymPy release. We have many unit tests >>> in the mechanics package, but that example provides a minimal example that >>> exercises almost all of the functionality in the Kane class. I'd personally >>> like a suite of benchmark multibody dynamics problems with known results >>> that we could run to exercise the mechanics package, with this example as >>> one. >>> >>> But regardless, it is fine to skip it now because it obviously is not >>> working and we haven't fixed that. >>> >>> >>> Jason >>> moorepants.info >>> +01 530-601-9791 >>> >>> On Fri, Sep 26, 2014 at 11:21 AM, James Crist <crist...@umn.edu> wrote: >>> >>>> The `test_kane3` test I think is not a good test in general. It's going >>>> to run slow, and personally I've never had it complete (although my laptop >>>> is slow). It's a fine example of the capabilities of our library, but as a >>>> test - not so good. I doubt the rest of the mechanics team would agree with >>>> me on removing it though, so perhaps a veryslow tag should be used? >>>> >>>> - Jim >>>> >>>> >>>> On Friday, September 26, 2014 10:13:33 AM UTC-5, Peter Brady wrote: >>>>> >>>>> Thanks to some recent PR's I was able to run the slow tests under >>>>> py.test which supports a --durations flags. Here are the slowest of the >>>>> slow tests: >>>>> >>>>> ========================== slowest 15 test durations >>>>> =========================== >>>>> 6745.05s call sympy/physics/mechanics/tests/test_kane3.py:: >>>>> test_bicycle >>>>> 2083.29s call sympy/solvers/tests/test_ode.py::test_1st_exact2 >>>>> 1522.66s call sympy/integrals/tests/test_failing_integrals.py:: >>>>> test_issue_4540 >>>>> 1221.33s call sympy/core/tests/test_wester.py::test_W25 >>>>> 1090.84s call sympy/integrals/tests/test_failing_integrals.py:: >>>>> test_issue_4941 >>>>> 963.31s call sympy/solvers/tests/test_solvers.py::test_issue_6528 >>>>> 437.72s call sympy/integrals/tests/test_heurisch.py::test_pmint_ >>>>> WrightOmega >>>>> 411.27s call sympy/core/tests/test_wester.py::test_W24 >>>>> 138.77s call sympy/polys/tests/test_groebnertools.py::test_ >>>>> benchmark_czichowski_buchberger >>>>> 111.14s call sympy/integrals/tests/test_failing_integrals.py:: >>>>> test_issue_4891 >>>>> 108.96s call sympy/functions/elementary/tests/test_trigonometric. >>>>> py::test_tancot_rewrite_sqrt >>>>> 90.53s call sympy/stats/tests/test_finite_rv.py:: >>>>> test_binomial_symbolic >>>>> 67.98s call sympy/core/tests/test_wester.py::test_W6 >>>>> 62.59s call sympy/simplify/tests/test_hyperexpand.py::test_ >>>>> prudnikov_3 >>>>> 48.68s call sympy/simplify/tests/test_hyperexpand.py::test_ >>>>> prudnikov_8 >>>>> >>>>> >>>>> Currently the timeout limit is set to 180s on the slow tests and due >>>>> to issues with travis, the signal does not appear to always be respected >>>>> leading to travis timeout errors. I think there are two good solutions to >>>>> this problem: >>>>> >>>>> 1. Mark all the tests which report runtimes over 180s as skip >>>>> >>>>> 2. Introduce a new mark 'veryslow` or something like that so that >>>>> these tests can still be run if someone chooses by setting the appropriate >>>>> flag. >>>>> >>>>> Any thoughts? Votes? >>>>> >>>> -- >>>> You received this message because you are subscribed to the Google >>>> Groups "sympy" group. >>>> To unsubscribe from this group and stop receiving emails from it, send >>>> an email to sympy+unsubscr...@googlegroups.com. >>>> To post to this group, send email to sympy@googlegroups.com. >>>> Visit this group at http://groups.google.com/group/sympy. >>>> To view this discussion on the web visit >>>> https://groups.google.com/d/msgid/sympy/97a9bf63-6db9-4b94-870d-ebbeeb59c7d9%40googlegroups.com >>>> <https://groups.google.com/d/msgid/sympy/97a9bf63-6db9-4b94-870d-ebbeeb59c7d9%40googlegroups.com?utm_medium=email&utm_source=footer> >>>> . >>>> >>>> For more options, visit https://groups.google.com/d/optout. >>>> >>> >>> >> -- >> You received this message because you are subscribed to a topic in the >> Google Groups "sympy" group. >> To unsubscribe from this topic, visit >> https://groups.google.com/d/topic/sympy/xDM4NHWP9cA/unsubscribe. >> To unsubscribe from this group and all its topics, send an email to >> sympy+unsubscr...@googlegroups.com. >> To post to this group, send email to sympy@googlegroups.com. >> Visit this group at http://groups.google.com/group/sympy. >> To view this discussion on the web visit >> https://groups.google.com/d/msgid/sympy/CAP7f1AjrdBMsW8WMTT_qbLhUEDgG7MrXQO6gL1jy3Pj6eOVJaQ%40mail.gmail.com >> <https://groups.google.com/d/msgid/sympy/CAP7f1AjrdBMsW8WMTT_qbLhUEDgG7MrXQO6gL1jy3Pj6eOVJaQ%40mail.gmail.com?utm_medium=email&utm_source=footer> >> . >> >> For more options, visit https://groups.google.com/d/optout. >> > > -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to sympy+unsubscr...@googlegroups.com. To post to this group, send email to sympy@googlegroups.com. Visit this group at http://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CALoNiQf-jXbVNGSE6Z6ZM%3D1uo3RE%3Dev3uTLrZHShNiAB187UPQ%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.