system.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445
  1. from sympy.core.backend import eye, Matrix, zeros
  2. from sympy.physics.mechanics import dynamicsymbols
  3. from sympy.physics.mechanics.functions import find_dynamicsymbols
  4. __all__ = ['SymbolicSystem']
  5. class SymbolicSystem:
  6. """SymbolicSystem is a class that contains all the information about a
  7. system in a symbolic format such as the equations of motions and the bodies
  8. and loads in the system.
  9. There are three ways that the equations of motion can be described for
  10. Symbolic System:
  11. [1] Explicit form where the kinematics and dynamics are combined
  12. x' = F_1(x, t, r, p)
  13. [2] Implicit form where the kinematics and dynamics are combined
  14. M_2(x, p) x' = F_2(x, t, r, p)
  15. [3] Implicit form where the kinematics and dynamics are separate
  16. M_3(q, p) u' = F_3(q, u, t, r, p)
  17. q' = G(q, u, t, r, p)
  18. where
  19. x : states, e.g. [q, u]
  20. t : time
  21. r : specified (exogenous) inputs
  22. p : constants
  23. q : generalized coordinates
  24. u : generalized speeds
  25. F_1 : right hand side of the combined equations in explicit form
  26. F_2 : right hand side of the combined equations in implicit form
  27. F_3 : right hand side of the dynamical equations in implicit form
  28. M_2 : mass matrix of the combined equations in implicit form
  29. M_3 : mass matrix of the dynamical equations in implicit form
  30. G : right hand side of the kinematical differential equations
  31. Parameters
  32. ==========
  33. coord_states : ordered iterable of functions of time
  34. This input will either be a collection of the coordinates or states
  35. of the system depending on whether or not the speeds are also
  36. given. If speeds are specified this input will be assumed to
  37. be the coordinates otherwise this input will be assumed to
  38. be the states.
  39. right_hand_side : Matrix
  40. This variable is the right hand side of the equations of motion in
  41. any of the forms. The specific form will be assumed depending on
  42. whether a mass matrix or coordinate derivatives are given.
  43. speeds : ordered iterable of functions of time, optional
  44. This is a collection of the generalized speeds of the system. If
  45. given it will be assumed that the first argument (coord_states)
  46. will represent the generalized coordinates of the system.
  47. mass_matrix : Matrix, optional
  48. The matrix of the implicit forms of the equations of motion (forms
  49. [2] and [3]). The distinction between the forms is determined by
  50. whether or not the coordinate derivatives are passed in. If
  51. they are given form [3] will be assumed otherwise form [2] is
  52. assumed.
  53. coordinate_derivatives : Matrix, optional
  54. The right hand side of the kinematical equations in explicit form.
  55. If given it will be assumed that the equations of motion are being
  56. entered in form [3].
  57. alg_con : Iterable, optional
  58. The indexes of the rows in the equations of motion that contain
  59. algebraic constraints instead of differential equations. If the
  60. equations are input in form [3], it will be assumed the indexes are
  61. referencing the mass_matrix/right_hand_side combination and not the
  62. coordinate_derivatives.
  63. output_eqns : Dictionary, optional
  64. Any output equations that are desired to be tracked are stored in a
  65. dictionary where the key corresponds to the name given for the
  66. specific equation and the value is the equation itself in symbolic
  67. form
  68. coord_idxs : Iterable, optional
  69. If coord_states corresponds to the states rather than the
  70. coordinates this variable will tell SymbolicSystem which indexes of
  71. the states correspond to generalized coordinates.
  72. speed_idxs : Iterable, optional
  73. If coord_states corresponds to the states rather than the
  74. coordinates this variable will tell SymbolicSystem which indexes of
  75. the states correspond to generalized speeds.
  76. bodies : iterable of Body/Rigidbody objects, optional
  77. Iterable containing the bodies of the system
  78. loads : iterable of load instances (described below), optional
  79. Iterable containing the loads of the system where forces are given
  80. by (point of application, force vector) and torques are given by
  81. (reference frame acting upon, torque vector). Ex [(point, force),
  82. (ref_frame, torque)]
  83. Attributes
  84. ==========
  85. coordinates : Matrix, shape(n, 1)
  86. This is a matrix containing the generalized coordinates of the system
  87. speeds : Matrix, shape(m, 1)
  88. This is a matrix containing the generalized speeds of the system
  89. states : Matrix, shape(o, 1)
  90. This is a matrix containing the state variables of the system
  91. alg_con : List
  92. This list contains the indices of the algebraic constraints in the
  93. combined equations of motion. The presence of these constraints
  94. requires that a DAE solver be used instead of an ODE solver.
  95. If the system is given in form [3] the alg_con variable will be
  96. adjusted such that it is a representation of the combined kinematics
  97. and dynamics thus make sure it always matches the mass matrix
  98. entered.
  99. dyn_implicit_mat : Matrix, shape(m, m)
  100. This is the M matrix in form [3] of the equations of motion (the mass
  101. matrix or generalized inertia matrix of the dynamical equations of
  102. motion in implicit form).
  103. dyn_implicit_rhs : Matrix, shape(m, 1)
  104. This is the F vector in form [3] of the equations of motion (the right
  105. hand side of the dynamical equations of motion in implicit form).
  106. comb_implicit_mat : Matrix, shape(o, o)
  107. This is the M matrix in form [2] of the equations of motion.
  108. This matrix contains a block diagonal structure where the top
  109. left block (the first rows) represent the matrix in the
  110. implicit form of the kinematical equations and the bottom right
  111. block (the last rows) represent the matrix in the implicit form
  112. of the dynamical equations.
  113. comb_implicit_rhs : Matrix, shape(o, 1)
  114. This is the F vector in form [2] of the equations of motion. The top
  115. part of the vector represents the right hand side of the implicit form
  116. of the kinemaical equations and the bottom of the vector represents the
  117. right hand side of the implicit form of the dynamical equations of
  118. motion.
  119. comb_explicit_rhs : Matrix, shape(o, 1)
  120. This vector represents the right hand side of the combined equations of
  121. motion in explicit form (form [1] from above).
  122. kin_explicit_rhs : Matrix, shape(m, 1)
  123. This is the right hand side of the explicit form of the kinematical
  124. equations of motion as can be seen in form [3] (the G matrix).
  125. output_eqns : Dictionary
  126. If output equations were given they are stored in a dictionary where
  127. the key corresponds to the name given for the specific equation and
  128. the value is the equation itself in symbolic form
  129. bodies : Tuple
  130. If the bodies in the system were given they are stored in a tuple for
  131. future access
  132. loads : Tuple
  133. If the loads in the system were given they are stored in a tuple for
  134. future access. This includes forces and torques where forces are given
  135. by (point of application, force vector) and torques are given by
  136. (reference frame acted upon, torque vector).
  137. Example
  138. =======
  139. As a simple example, the dynamics of a simple pendulum will be input into a
  140. SymbolicSystem object manually. First some imports will be needed and then
  141. symbols will be set up for the length of the pendulum (l), mass at the end
  142. of the pendulum (m), and a constant for gravity (g). ::
  143. >>> from sympy import Matrix, sin, symbols
  144. >>> from sympy.physics.mechanics import dynamicsymbols, SymbolicSystem
  145. >>> l, m, g = symbols('l m g')
  146. The system will be defined by an angle of theta from the vertical and a
  147. generalized speed of omega will be used where omega = theta_dot. ::
  148. >>> theta, omega = dynamicsymbols('theta omega')
  149. Now the equations of motion are ready to be formed and passed to the
  150. SymbolicSystem object. ::
  151. >>> kin_explicit_rhs = Matrix([omega])
  152. >>> dyn_implicit_mat = Matrix([l**2 * m])
  153. >>> dyn_implicit_rhs = Matrix([-g * l * m * sin(theta)])
  154. >>> symsystem = SymbolicSystem([theta], dyn_implicit_rhs, [omega],
  155. ... dyn_implicit_mat)
  156. Notes
  157. =====
  158. m : number of generalized speeds
  159. n : number of generalized coordinates
  160. o : number of states
  161. """
  162. def __init__(self, coord_states, right_hand_side, speeds=None,
  163. mass_matrix=None, coordinate_derivatives=None, alg_con=None,
  164. output_eqns={}, coord_idxs=None, speed_idxs=None, bodies=None,
  165. loads=None):
  166. """Initializes a SymbolicSystem object"""
  167. # Extract information on speeds, coordinates and states
  168. if speeds is None:
  169. self._states = Matrix(coord_states)
  170. if coord_idxs is None:
  171. self._coordinates = None
  172. else:
  173. coords = [coord_states[i] for i in coord_idxs]
  174. self._coordinates = Matrix(coords)
  175. if speed_idxs is None:
  176. self._speeds = None
  177. else:
  178. speeds_inter = [coord_states[i] for i in speed_idxs]
  179. self._speeds = Matrix(speeds_inter)
  180. else:
  181. self._coordinates = Matrix(coord_states)
  182. self._speeds = Matrix(speeds)
  183. self._states = self._coordinates.col_join(self._speeds)
  184. # Extract equations of motion form
  185. if coordinate_derivatives is not None:
  186. self._kin_explicit_rhs = coordinate_derivatives
  187. self._dyn_implicit_rhs = right_hand_side
  188. self._dyn_implicit_mat = mass_matrix
  189. self._comb_implicit_rhs = None
  190. self._comb_implicit_mat = None
  191. self._comb_explicit_rhs = None
  192. elif mass_matrix is not None:
  193. self._kin_explicit_rhs = None
  194. self._dyn_implicit_rhs = None
  195. self._dyn_implicit_mat = None
  196. self._comb_implicit_rhs = right_hand_side
  197. self._comb_implicit_mat = mass_matrix
  198. self._comb_explicit_rhs = None
  199. else:
  200. self._kin_explicit_rhs = None
  201. self._dyn_implicit_rhs = None
  202. self._dyn_implicit_mat = None
  203. self._comb_implicit_rhs = None
  204. self._comb_implicit_mat = None
  205. self._comb_explicit_rhs = right_hand_side
  206. # Set the remainder of the inputs as instance attributes
  207. if alg_con is not None and coordinate_derivatives is not None:
  208. alg_con = [i + len(coordinate_derivatives) for i in alg_con]
  209. self._alg_con = alg_con
  210. self.output_eqns = output_eqns
  211. # Change the body and loads iterables to tuples if they are not tuples
  212. # already
  213. if not isinstance(bodies, tuple) and bodies is not None:
  214. bodies = tuple(bodies)
  215. if not isinstance(loads, tuple) and loads is not None:
  216. loads = tuple(loads)
  217. self._bodies = bodies
  218. self._loads = loads
  219. @property
  220. def coordinates(self):
  221. """Returns the column matrix of the generalized coordinates"""
  222. if self._coordinates is None:
  223. raise AttributeError("The coordinates were not specified.")
  224. else:
  225. return self._coordinates
  226. @property
  227. def speeds(self):
  228. """Returns the column matrix of generalized speeds"""
  229. if self._speeds is None:
  230. raise AttributeError("The speeds were not specified.")
  231. else:
  232. return self._speeds
  233. @property
  234. def states(self):
  235. """Returns the column matrix of the state variables"""
  236. return self._states
  237. @property
  238. def alg_con(self):
  239. """Returns a list with the indices of the rows containing algebraic
  240. constraints in the combined form of the equations of motion"""
  241. return self._alg_con
  242. @property
  243. def dyn_implicit_mat(self):
  244. """Returns the matrix, M, corresponding to the dynamic equations in
  245. implicit form, M x' = F, where the kinematical equations are not
  246. included"""
  247. if self._dyn_implicit_mat is None:
  248. raise AttributeError("dyn_implicit_mat is not specified for "
  249. "equations of motion form [1] or [2].")
  250. else:
  251. return self._dyn_implicit_mat
  252. @property
  253. def dyn_implicit_rhs(self):
  254. """Returns the column matrix, F, corresponding to the dynamic equations
  255. in implicit form, M x' = F, where the kinematical equations are not
  256. included"""
  257. if self._dyn_implicit_rhs is None:
  258. raise AttributeError("dyn_implicit_rhs is not specified for "
  259. "equations of motion form [1] or [2].")
  260. else:
  261. return self._dyn_implicit_rhs
  262. @property
  263. def comb_implicit_mat(self):
  264. """Returns the matrix, M, corresponding to the equations of motion in
  265. implicit form (form [2]), M x' = F, where the kinematical equations are
  266. included"""
  267. if self._comb_implicit_mat is None:
  268. if self._dyn_implicit_mat is not None:
  269. num_kin_eqns = len(self._kin_explicit_rhs)
  270. num_dyn_eqns = len(self._dyn_implicit_rhs)
  271. zeros1 = zeros(num_kin_eqns, num_dyn_eqns)
  272. zeros2 = zeros(num_dyn_eqns, num_kin_eqns)
  273. inter1 = eye(num_kin_eqns).row_join(zeros1)
  274. inter2 = zeros2.row_join(self._dyn_implicit_mat)
  275. self._comb_implicit_mat = inter1.col_join(inter2)
  276. return self._comb_implicit_mat
  277. else:
  278. raise AttributeError("comb_implicit_mat is not specified for "
  279. "equations of motion form [1].")
  280. else:
  281. return self._comb_implicit_mat
  282. @property
  283. def comb_implicit_rhs(self):
  284. """Returns the column matrix, F, corresponding to the equations of
  285. motion in implicit form (form [2]), M x' = F, where the kinematical
  286. equations are included"""
  287. if self._comb_implicit_rhs is None:
  288. if self._dyn_implicit_rhs is not None:
  289. kin_inter = self._kin_explicit_rhs
  290. dyn_inter = self._dyn_implicit_rhs
  291. self._comb_implicit_rhs = kin_inter.col_join(dyn_inter)
  292. return self._comb_implicit_rhs
  293. else:
  294. raise AttributeError("comb_implicit_mat is not specified for "
  295. "equations of motion in form [1].")
  296. else:
  297. return self._comb_implicit_rhs
  298. def compute_explicit_form(self):
  299. """If the explicit right hand side of the combined equations of motion
  300. is to provided upon initialization, this method will calculate it. This
  301. calculation can potentially take awhile to compute."""
  302. if self._comb_explicit_rhs is not None:
  303. raise AttributeError("comb_explicit_rhs is already formed.")
  304. inter1 = getattr(self, 'kin_explicit_rhs', None)
  305. if inter1 is not None:
  306. inter2 = self._dyn_implicit_mat.LUsolve(self._dyn_implicit_rhs)
  307. out = inter1.col_join(inter2)
  308. else:
  309. out = self._comb_implicit_mat.LUsolve(self._comb_implicit_rhs)
  310. self._comb_explicit_rhs = out
  311. @property
  312. def comb_explicit_rhs(self):
  313. """Returns the right hand side of the equations of motion in explicit
  314. form, x' = F, where the kinematical equations are included"""
  315. if self._comb_explicit_rhs is None:
  316. raise AttributeError("Please run .combute_explicit_form before "
  317. "attempting to access comb_explicit_rhs.")
  318. else:
  319. return self._comb_explicit_rhs
  320. @property
  321. def kin_explicit_rhs(self):
  322. """Returns the right hand side of the kinematical equations in explicit
  323. form, q' = G"""
  324. if self._kin_explicit_rhs is None:
  325. raise AttributeError("kin_explicit_rhs is not specified for "
  326. "equations of motion form [1] or [2].")
  327. else:
  328. return self._kin_explicit_rhs
  329. def dynamic_symbols(self):
  330. """Returns a column matrix containing all of the symbols in the system
  331. that depend on time"""
  332. # Create a list of all of the expressions in the equations of motion
  333. if self._comb_explicit_rhs is None:
  334. eom_expressions = (self.comb_implicit_mat[:] +
  335. self.comb_implicit_rhs[:])
  336. else:
  337. eom_expressions = (self._comb_explicit_rhs[:])
  338. functions_of_time = set()
  339. for expr in eom_expressions:
  340. functions_of_time = functions_of_time.union(
  341. find_dynamicsymbols(expr))
  342. functions_of_time = functions_of_time.union(self._states)
  343. return tuple(functions_of_time)
  344. def constant_symbols(self):
  345. """Returns a column matrix containing all of the symbols in the system
  346. that do not depend on time"""
  347. # Create a list of all of the expressions in the equations of motion
  348. if self._comb_explicit_rhs is None:
  349. eom_expressions = (self.comb_implicit_mat[:] +
  350. self.comb_implicit_rhs[:])
  351. else:
  352. eom_expressions = (self._comb_explicit_rhs[:])
  353. constants = set()
  354. for expr in eom_expressions:
  355. constants = constants.union(expr.free_symbols)
  356. constants.remove(dynamicsymbols._t)
  357. return tuple(constants)
  358. @property
  359. def bodies(self):
  360. """Returns the bodies in the system"""
  361. if self._bodies is None:
  362. raise AttributeError("bodies were not specified for the system.")
  363. else:
  364. return self._bodies
  365. @property
  366. def loads(self):
  367. """Returns the loads in the system"""
  368. if self._loads is None:
  369. raise AttributeError("loads were not specified for the system.")
  370. else:
  371. return self._loads