linearize - q' ' in A matrix

[Bas]

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 

-----------------

#%% 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])

[Bas]

Hello, 

It think that I found what goes wrong: SymPy seems to return the linearization too quickly.  Of course, I might have made an error in the derivation below. 

I use a simpler system: a cart with a pendulum (code below). The position of the cart (x), and the angle of the pendulum (q) are used as the generalised coordinates. If I linearise this system with SymPy around q=q0, I find in the linearised equations

x'' + c1.q.q'' = c2.Fx

q'' + c3.q.x'' = 0

This system is not linear in (x'',q'') due to the multipication of the states. In this equation ci are constants. If I start from the equations of motion, and do the linearisation myself, I find the following.

 

Is this something that ought to be modified in Sympy, or am I making a mistake? 

Hope it helps.

Bas

The code used in Python:

-------------------

from numpy import pi

from sympy import *

from sympy.physics.mechanics import *

#%% declare the parameters and variables

mp, mc, l, g = symbols('mp, mc, l, g')

# position x, and angle used as coordinates

x, q = dynamicsymbols('x, q')

xd, qd = dynamicsymbols('x, q', 1)

Fx = dynamicsymbols('Fx')

#%% setup the reference frames 

O = ReferenceFrame('O')

A = O.orientnew('A', 'Axis', [q, O.z])

#%% setup the points

Op = Point('Op')

Ap = Op.locatenew('Ap', x*O.x)

Bp = Ap.locatenew('Bp', l*A.x)

# velocity

Op.set_vel(O,0)

Ap.set_vel(O, Ap.pos_from(Op).dt(O))

Bp.v2pt_theory(Ap, O, A)

#%% give some body to it

Apa = Particle('Apa', Ap, mc)

Bpa = Particle('Bpa', Bp, mp)

Bpa.potential_energy = mp*g* Bp.pos_from(Op).express(O).dot(O.y)

#%% put it into Lagrangian

force = [(Ap, Fx*O.x)]

L = Lagrangian(O, Apa, Bpa)

LM = LagrangesMethod(L, [x, q], forcelist=force, frame=O)

LM.form_lagranges_equations()

#%% linearise

op_point={x:0,q:0,qd:0,xd:0}

AA, BB, rr = LM.linearize(q_ind=[x, q], qd_ind=[xd, qd], op_point=op_point, A_and_B=True)

Op woensdag 4 november 2020 om 10:07:56 UTC+1 schreef Bas:

[Jason Moore]

Bas,

I haven't had time to look at this carefully. I'll try to find time this week sometime.

A quick comment, would be to use find_dynamicsymbols() on the equations of motion to make sure you are getting equations with the appropriate variables. KanesMethod and LagrangesMethod may not get them quite to the best final form for subsequent steps (like linearization).

If your system does not have motion contractions or configuration constraints, linearization is pretty straight forward without the linearize() function. You can check the results. The notes for lecture 20 here: https://moorepants.github.io/mae223/schedule.html show how to do that.

It may very well be possible that linearize() isnt' working correctly, but I'll have to find some time to look closely.

Hope that helps.

Jason

moorepants.info
+01 530-601-9791

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/8c30ce07-b76b-4aab-9614-09db09181034n%40googlegroups.com.

[Bas]

Hello Jason,

Thanks for you reaction. A quick check on the dynamic variables for the cart-pendulum system gave x, x'', q, q', q'' and Fx. This seems to be correct to me. The second second derivates and Fx are not states, but are dynamic variables. 

For the simple system I linearised it by hand (see previous post), and that gave me the impression that the linearisation of sin(q)q'' is the origin of the problem: The sine is expanded into Taylor components, but the multiplication with q'' is not expanded afterwards. 

Hope it help, and thanks 

Bas

Op zondag 8 november 2020 om 11:33:11 UTC+1 schreef moore...@gmail.com: