jointsmethod.py 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278
  1. from sympy.physics.mechanics import (Body, Lagrangian, KanesMethod, LagrangesMethod,
  2. RigidBody, Particle)
  3. from sympy.physics.mechanics.method import _Methods
  4. __all__ = ['JointsMethod']
  5. class JointsMethod(_Methods):
  6. """Method for formulating the equations of motion using a set of interconnected bodies with joints.
  7. Parameters
  8. ==========
  9. newtonion : Body or ReferenceFrame
  10. The newtonion(inertial) frame.
  11. *joints : Joint
  12. The joints in the system
  13. Attributes
  14. ==========
  15. q, u : iterable
  16. Iterable of the generalized coordinates and speeds
  17. bodies : iterable
  18. Iterable of Body objects in the system.
  19. loads : iterable
  20. Iterable of (Point, vector) or (ReferenceFrame, vector) tuples
  21. describing the forces on the system.
  22. mass_matrix : Matrix, shape(n, n)
  23. The system's mass matrix
  24. forcing : Matrix, shape(n, 1)
  25. The system's forcing vector
  26. mass_matrix_full : Matrix, shape(2*n, 2*n)
  27. The "mass matrix" for the u's and q's
  28. forcing_full : Matrix, shape(2*n, 1)
  29. The "forcing vector" for the u's and q's
  30. method : KanesMethod or Lagrange's method
  31. Method's object.
  32. kdes : iterable
  33. Iterable of kde in they system.
  34. Examples
  35. ========
  36. This is a simple example for a one degree of freedom translational
  37. spring-mass-damper.
  38. >>> from sympy import symbols
  39. >>> from sympy.physics.mechanics import Body, JointsMethod, PrismaticJoint
  40. >>> from sympy.physics.vector import dynamicsymbols
  41. >>> c, k = symbols('c k')
  42. >>> x, v = dynamicsymbols('x v')
  43. >>> wall = Body('W')
  44. >>> body = Body('B')
  45. >>> J = PrismaticJoint('J', wall, body, coordinates=x, speeds=v)
  46. >>> wall.apply_force(c*v*wall.x, reaction_body=body)
  47. >>> wall.apply_force(k*x*wall.x, reaction_body=body)
  48. >>> method = JointsMethod(wall, J)
  49. >>> method.form_eoms()
  50. Matrix([[-B_mass*Derivative(v(t), t) - c*v(t) - k*x(t)]])
  51. >>> M = method.mass_matrix_full
  52. >>> F = method.forcing_full
  53. >>> rhs = M.LUsolve(F)
  54. >>> rhs
  55. Matrix([
  56. [ v(t)],
  57. [(-c*v(t) - k*x(t))/B_mass]])
  58. Notes
  59. =====
  60. ``JointsMethod`` currently only works with systems that do not have any
  61. configuration or motion constraints.
  62. """
  63. def __init__(self, newtonion, *joints):
  64. if isinstance(newtonion, Body):
  65. self.frame = newtonion.frame
  66. else:
  67. self.frame = newtonion
  68. self._joints = joints
  69. self._bodies = self._generate_bodylist()
  70. self._loads = self._generate_loadlist()
  71. self._q = self._generate_q()
  72. self._u = self._generate_u()
  73. self._kdes = self._generate_kdes()
  74. self._method = None
  75. @property
  76. def bodies(self):
  77. """List of bodies in they system."""
  78. return self._bodies
  79. @property
  80. def loads(self):
  81. """List of loads on the system."""
  82. return self._loads
  83. @property
  84. def q(self):
  85. """List of the generalized coordinates."""
  86. return self._q
  87. @property
  88. def u(self):
  89. """List of the generalized speeds."""
  90. return self._u
  91. @property
  92. def kdes(self):
  93. """List of the generalized coordinates."""
  94. return self._kdes
  95. @property
  96. def forcing_full(self):
  97. """The "forcing vector" for the u's and q's."""
  98. return self.method.forcing_full
  99. @property
  100. def mass_matrix_full(self):
  101. """The "mass matrix" for the u's and q's."""
  102. return self.method.mass_matrix_full
  103. @property
  104. def mass_matrix(self):
  105. """The system's mass matrix."""
  106. return self.method.mass_matrix
  107. @property
  108. def forcing(self):
  109. """The system's forcing vector."""
  110. return self.method.forcing
  111. @property
  112. def method(self):
  113. """Object of method used to form equations of systems."""
  114. return self._method
  115. def _generate_bodylist(self):
  116. bodies = []
  117. for joint in self._joints:
  118. if joint.child not in bodies:
  119. bodies.append(joint.child)
  120. if joint.parent not in bodies:
  121. bodies.append(joint.parent)
  122. return bodies
  123. def _generate_loadlist(self):
  124. load_list = []
  125. for body in self.bodies:
  126. load_list.extend(body.loads)
  127. return load_list
  128. def _generate_q(self):
  129. q_ind = []
  130. for joint in self._joints:
  131. for coordinate in joint.coordinates:
  132. if coordinate in q_ind:
  133. raise ValueError('Coordinates of joints should be unique.')
  134. q_ind.append(coordinate)
  135. return q_ind
  136. def _generate_u(self):
  137. u_ind = []
  138. for joint in self._joints:
  139. for speed in joint.speeds:
  140. if speed in u_ind:
  141. raise ValueError('Speeds of joints should be unique.')
  142. u_ind.append(speed)
  143. return u_ind
  144. def _generate_kdes(self):
  145. kd_ind = []
  146. for joint in self._joints:
  147. kd_ind.extend(joint.kdes)
  148. return kd_ind
  149. def _convert_bodies(self):
  150. # Convert `Body` to `Particle` and `RigidBody`
  151. bodylist = []
  152. for body in self.bodies:
  153. if body.is_rigidbody:
  154. rb = RigidBody(body.name, body.masscenter, body.frame, body.mass,
  155. (body.central_inertia, body.masscenter))
  156. rb.potential_energy = body.potential_energy
  157. bodylist.append(rb)
  158. else:
  159. part = Particle(body.name, body.masscenter, body.mass)
  160. part.potential_energy = body.potential_energy
  161. bodylist.append(part)
  162. return bodylist
  163. def form_eoms(self, method=KanesMethod):
  164. """Method to form system's equation of motions.
  165. Parameters
  166. ==========
  167. method : Class
  168. Class name of method.
  169. Returns
  170. ========
  171. Matrix
  172. Vector of equations of motions.
  173. Examples
  174. ========
  175. This is a simple example for a one degree of freedom translational
  176. spring-mass-damper.
  177. >>> from sympy import S, symbols
  178. >>> from sympy.physics.mechanics import LagrangesMethod, dynamicsymbols, Body
  179. >>> from sympy.physics.mechanics import PrismaticJoint, JointsMethod
  180. >>> q = dynamicsymbols('q')
  181. >>> qd = dynamicsymbols('q', 1)
  182. >>> m, k, b = symbols('m k b')
  183. >>> wall = Body('W')
  184. >>> part = Body('P', mass=m)
  185. >>> part.potential_energy = k * q**2 / S(2)
  186. >>> J = PrismaticJoint('J', wall, part, coordinates=q, speeds=qd)
  187. >>> wall.apply_force(b * qd * wall.x, reaction_body=part)
  188. >>> method = JointsMethod(wall, J)
  189. >>> method.form_eoms(LagrangesMethod)
  190. Matrix([[b*Derivative(q(t), t) + k*q(t) + m*Derivative(q(t), (t, 2))]])
  191. We can also solve for the states using the 'rhs' method.
  192. >>> method.rhs()
  193. Matrix([
  194. [ Derivative(q(t), t)],
  195. [(-b*Derivative(q(t), t) - k*q(t))/m]])
  196. """
  197. bodylist = self._convert_bodies()
  198. if issubclass(method, LagrangesMethod): #LagrangesMethod or similar
  199. L = Lagrangian(self.frame, *bodylist)
  200. self._method = method(L, self.q, self.loads, bodylist, self.frame)
  201. else: #KanesMethod or similar
  202. self._method = method(self.frame, q_ind=self.q, u_ind=self.u, kd_eqs=self.kdes,
  203. forcelist=self.loads, bodies=bodylist)
  204. soln = self.method._form_eoms()
  205. return soln
  206. def rhs(self, inv_method=None):
  207. """Returns equations that can be solved numerically.
  208. Parameters
  209. ==========
  210. inv_method : str
  211. The specific sympy inverse matrix calculation method to use. For a
  212. list of valid methods, see
  213. :meth:`~sympy.matrices.matrices.MatrixBase.inv`
  214. Returns
  215. ========
  216. Matrix
  217. Numerically solveable equations.
  218. See Also
  219. ========
  220. sympy.physics.mechanics.KanesMethod.rhs():
  221. KanesMethod's rhs function.
  222. sympy.physics.mechanics.LagrangesMethod.rhs():
  223. LagrangesMethod's rhs function.
  224. """
  225. return self.method.rhs(inv_method=inv_method)