Hello all,
First of all, thanks for all the great work that has already been done by
the developers and the community. I recently discovered SymPy, and I use it
frequently now for my job as control engineer. The discussions in this
group and the documentation helped me get a quick start.
However, I do get some behaviour that I cannot explain, and didn't expect.
I try to model a crane on a ship by means of the Lagrangian method (see
below for the code and a schematic). When I linearize the equations of
motion, the A matrix contains accelerations of the independent coordinates.
And these are exactly the variables I try to solve: they are also in the
left-hand-side of the linearized equation M x' = Ax + Br. I expected to
get only numerical values in the matrices M, A, and B.
I hope somebody can tell me why these elements are in there, or/and how I
can get them out of the A matrix?
The system that I want to model, as well as the code, is shown below. The
'ship' is a rigid body, with coordinates xs, ys, and qs. Connected to it is
a point Pp. This point is rigidly connected with the ship. It has no mass,
and is not included in the Lagrangian.
Then there is a pendulum connected to this point. The angle it makes with
respect to the ship angle is qp. Along it's x-axis there is a point Qp,
with some mass. I include the rigid body and the mass at the end of the
pendulum in the Lagrangian, and use [xs, ys, qs, qp] as the generalized
coordinates. There are some forces acting on the ship to keep it in place
(spring/dampers)
For the linearization, I insert the parametric values and the operating
point to minimize calculation time This should also result in fully numeric
matrices. The same approach works great for a single/double pendulum to get
numerical matrices, but here not. For example, element AA[4,2]
= 35*sqrt(2)*qs' '
With the help of pydy I can extract the set of first order differential
equations and simulate them. (this is not in the code below in order to
keep it minimal). These simulations seem correct.
Can anyone point me in the right direction to get the acceleration
variables out of the A matrix??
Thanks for any help
Bas
[image: rect4519.png]
-----------------
#%% import the packages I rely on
import numpy as np
from numpy import pi
from sympy import *
from sympy.physics.mechanics import *
#%% declare the parameters and variables
lc, lp1, ms, m1, Is, K, Ky, bs, bx, by, g = \
symbols('lc, lp1, ms, m1, Is, K, Ky, bs, bx, by, g'
,positive=True, real=True)
# crane angle, considered a parameter
qc = symbols('qc')
# the dynamic variables
qs, qp = dynamicsymbols('qs, qp')
qsd, qpd = dynamicsymbols('qs, qp', 1)
xs, ys = dynamicsymbols('xs, ys')
xsd, ysd = dynamicsymbols('xs, ys', 1)
# numerical values
par = {lc:7, lp1:2, ms:50, m1:10, Is:100, K:10000, Ky:1000, bs:500, bx:10,
by:100, g:9.81, qc:pi/4}
#%% setup the reference frames (as it is 2D, we rotate all around O.z)
# note that they are connected to each other
O = ReferenceFrame('O')
S = O.orientnew('S', 'Axis', [qs, O.z])
P = S.orientnew('P', 'Axis', [qp, O.z])
#%% setup the points
# position
Op = Point('Op')
Sp = Op.locatenew('Sp', xs*O.x + ys*O.y)
Pp = Sp.locatenew('Pp', lc*cos(qc)*S.x + lc*sin(qc)*S.y)
Qp = Pp.locatenew('Qp', lp1*P.x)
# velocity
Op.set_vel(O,0)
Sp.set_vel(O, Sp.pos_from(Op).dt(O))
Pp.v2pt_theory(Sp, O, S)
Qp.v2pt_theory(Pp, O, P)
#%% give some body to it
Izz = Is*outer(S.z, S.z)
Sbody = RigidBody('Sb', Sp, S, ms, (Izz, Sp))
Qparticle = Particle('Qparticle', Qp, m1)
Qparticle.potential_energy =
trigsimp(m1*g*Qp.pos_from(Op).express(O).dot(O.y))
#%% put it into Lagrangian
force = [(S, -K*qs*S.z -bs*qsd*S.z), (Sp, -Ky*ys*O.y - by*ysd*O.y), (Sp,
-bx*xsd*O.x)]
L = Lagrangian(O, Sbody, Qparticle)
LM = LagrangesMethod(L, [xs, ys, qs, qp], forcelist=force, frame=O)
LM.form_lagranges_equations()
#%% linearise
# first set the operating point and the numerical values
op_point = par.copy()
op_point[xs]= 0.0
op_point[ys]= 0.0
op_point[qs]= 0.0
op_point[qp]= -pi/2
op_point[xsd]= 0.0
op_point[ysd]= 0.0
op_point[qsd]= 0.0
op_point[qpd]= 0.0
# go linearise
MM, AA, BB, rr = LM.linearize(q_ind=LM.q, qd_ind=LM.u, op_point=op_point,
A_and_B=False)
# print unexpected result
vprint(AA[4:,:4])
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/1f22de5f-8d89-4375-8d31-97d67ab2624en%40googlegroups.com.
