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 [email protected]: > 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 <(530)%20601-9791> > > > On Sun, Nov 8, 2020 at 10:31 AM Bas <[email protected]> wrote: > >> 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. >> [image: sympy_linearise.png] >> 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: >> >>> 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/8c30ce07-b76b-4aab-9614-09db09181034n%40googlegroups.com >> >> <https://groups.google.com/d/msgid/sympy/8c30ce07-b76b-4aab-9614-09db09181034n%40googlegroups.com?utm_medium=email&utm_source=footer> >> . >> > -- 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/4ba67dac-d72f-49b1-801f-cadfc1759c52n%40googlegroups.com.
