GSoC Week 3: Generalized Linearizer Class

Posted on June 06, 2014

This week I refactored the Linearizer class to fit a general system of equations. This means that any of the following can be absent, and the linearization still is valid:

$$ \begin{array}[rl] \\f_{c}(q, t) &= 0_{lx1} \\ f_{v}(q, u, t) &= 0_{mx1} \\ f_{a}(q, \dot{q}, u, \dot{u}, t) &= 0_{mx1} \\ f_{0}(q, \dot{q}, t) + f_{1}(q, u, t) &= 0_{nx1} \\ f_{2}(q, \dot{u}, t) + f_{3}(q, \dot{q}, u, r, t) &= 0_{(o-m)x1} \\ \end{array} $$

with

$$ \begin{array}[rl] \\q, \dot{q} & \in \mathbb{R}^n \\ u, \dot{u} & \in \mathbb{R}^o \\ r & \in \mathbb{R}^s \end{array} $$

Note that vectors can be absent too. This means the system can be entirely $q$ or $u$ components. More on this later.

Much of the rest of the week was spent finishing up all methods related to the KanesMethod class. This was mostly cleaning up the code base, and fixing all the tests in place for the previous linearize method. As of this point, the linearization routines accomplish everything that was in place before, but in a more general and extensible way. This is all included in the LinearizerClass branch on my GitHub repo. There is a pull request on my own master branch open right now for code review here. I'll leave this up for another day or two before closing it and making all pull request to Sympy/master.

I also started work on linearizing a system of equations presented in matrix form. This still needs lots of code improvements (it's inefficient), but it currently works. Below is a demo.

The system is the same system I use for my research - a dual solenoid actuator. There are two dependent variables: lam1 and lam2. The dynamics are thus expressed by 3 dynamic differential equations, and two constraint equations.

In [1]:

from sympy import symbols, Matrix
from sympy.physics.mechanics import dynamicsymbols, mprint
from sympy.physics.mechanics.linearize import Linearizer, matrix_to_linearizer

# Create constant symbols
B1, B2 = symbols('B1, B2')
d1, d2, D = symbols('d1, d2, D')
R1, R2 = symbols('R1, R2')
m = symbols('m')
t = symbols('t')

# Create dynamicsymbols
i1, i2 = dynamicsymbols('i1, i2')
V1, V2 = dynamicsymbols('V1, V2')
x = dynamicsymbols('x')
xdot = x.diff()
lam1, lam2 = dynamicsymbols('lam1, lam2')
lam1dot, lam2dot = dynamicsymbols('lam1, lam2', 1)

# Define vectors. Because all speeds are derivatives of the coordinates,
# the system is best represented as just a vector of u. The dependent
# states are lam1 and lam2
u = Matrix([i1, i2, xdot, x, lam1, lam2])
udep = u[4:]

# Define system of equations
eom = Matrix([B1/(d1 + x)*i1 - lam1,
              B2/(D + d2 - x)*i2 - lam2,
              -lam1dot - R1*i1 + V1,
              -lam2dot - R2*i2 + V2,
              -xdot.diff() + (B2/(d2 + D - x)**2 * i2**2 - B1/(d1 + x)**2 * i1**2)/m])

# Perform the linearization
linearizer = matrix_to_linearizer(eom, u=u, udep=udep)
A, B = linearizer.linearize(A_and_B=True)

In [2]:

mprint(A)

Matrix([
[(B1*x' - R1*(d1 + x)**2)/(B1*(d1 + x)),                                               0,                        ((d1 + x)*i1' - 2*i1*x')/(d1 + x)**2,      2*i1/(d1 + x)],
[                                     0, -(B2*x' + R2*(D + d2 - x)**2)/(B2*(D + d2 - x)),               -((D + d2 - x)*i2' + 2*i2*x')/(D + d2 - x)**2, -2*i2/(D + d2 - x)],
[                                     0,                                               0,                                                           0,                  1],
[              -2*B1*i1/(m*(d1 + x)**2),                     2*B2*i2/(m*(D + d2 - x)**2), 2*B1*i1**2/(m*(d1 + x)**3) + 2*B2*i2**2/(m*(D + d2 - x)**3),                  0]])

In [3]:

mprint(B)

Matrix([
[(d1 + x)/B1,               0],
[          0, (D + d2 - x)/B2],
[          0,               0],
[          0,               0]])

You'll have to trust me that the above representations are the same ones that I derived for my research project. This is really exciting progress! If everything works out (and it looks like it will), it will be possible to derive a system of dynamic equations anywhere in sympy, and then linearize it using the functionality I implement in the linearize function.

Currently the system of equations is must formed as a matrix, and then passed to matrix_to_linearizer to return a Linearizer object. Once this code is cleaned up, it will be used internally by the general linearize function to convert the system to a Linearizer, perform the linearization, and return the results. I plan on implementing this functionality next week.

The code used above is found in my matrix2linearizer branch on GitHub here. As before, I made a pull request on my own rep, so people may review the code before i submit it to Sympy proper. The PR is here.