bowfinger.de

This is a short piece about what I think the two entities in the title have in common. And the way I found remarkably telling is through signing up to the European Rail Passengers Union, which by the way you should absolutely consider if you are using railways in Europe at all.

First things first: I wanted to sign up to the ERPU because I think they do crucial work to improve the passenger experience on rails in Europe – which sucks right now.

The ERPU is an association under Belgian law, you can sign up online for only 12€/year, and hey, in the payment process, amongst all the usual payment options, I get “iDEAL Wero”.

Wero you say? I can try this new European payment method? Nice! I’m hooked! So, I click.

And there’s a QR code saying “scan me with your banking app”.

Well. My banking app says, this is not a Wero link – and apparently it isn’t.

So I click on the other button that says: Or use internet banking, select your bank.

All right, so I click this.

And it lists … only Belgian banks.

So, what do I do? I take the plane. I pay using Visa/Mastercard/Amexco.

This is exactly what the European Rail experience is today:

You try. You try again. The connection is shoddy, you have no rights in case ticket 1 fails and you miss your train on ticket 2. So, you use Paypal. So, you take the plane.

Here it’s the same. There’s this European payment service Wero, that in principle could be used Europe wide. But implementation on all sides – bank, service provider – is crappy. It’s limited to national borders for yet another completely in-transparent reason. And so, you pay using credit card.

And it’s probably for the exact same reasons too. Like train operators, many banks have their turf, are safe because they’re the default and somewhat state backed (yes, I’m looking at you Sparkassen) and haven’t adapted to demand in decades. Like train regulation, regulators only seem to care what’s going on when shit is on fire (cf. sub prime loan crisis et al., or railway tracks crumbling under twice the load they were constructed for et al.). And like in train land, there’s the other agile, indirectly but heavily subsidised means of payment with a shitload of hidden downsides that you default to. You don’t think cars and planes have massive downsides? Maybe you also think, Mastercard, Visa, PayPal and Apple/Google Pay offer their services for free. Heck no, you simply don’t ever see the bill. And: Yes, et al indeed.

And this is not a dunk on the folks at ERPU. On the contrary: they’ve taken what is available in terms of payment services, integrated it nicely into their web page. Jon is doing actual field work travelling around all the crazy cross border connections that barely work in Europe. No. The horror starts once you leave their web page and enter into the payment service realm.

But it illustrates so very clearly why it’s so bad. And why work like the one of ERPU is needed. Hell, maybe we need more than the European Rail Passenger Union. We need a European Payment Users Union. How’s that?

A nifty feature of sympy is its code generator. This allows generating code in many languages from expressions derived in sympy. I’ve been using this for a while and came across the Common Subexpression Elimination (CSE) function (again), recently.

In earlier* versions of sympy I did not find it of great use. But now I was pleasantly surprised how much it matured. (* I cannot say for sure when I last tried using it; it must’ve been a long while ago.)

What is it and what is it good for?

Let’s say we want to generate code for a PI controller in C. The literature says its transfer function (in DIN notation) reads:

Then furthermore the representation for a discrete time system with time step becomes

In order to efficiently and flexibly calculate the response we want to implement it in a second order IIR (infinite impulse response) filter structure. This structure is basically a representation of as a rational function

so we need to take the expression apart such that we can determine the coefficients a_1, a_2, b_0, … (by convention, the problem is scaled such that a_0 = 1.

This is all very cumbersome to do manually (and I suck at doing this, especially as I insert quite embarrassing mistakes all the time). So, sympy to the rescue:

import sympy as sy

Kp, Ki, Ts = sy.symbols("K_p K_i T_s", real=True)
z = sy.symbols("z", complex=True)
G_z = Kp + Ki*Ts/2 + (Ki*Ts/2 - Kp)/z

display(G_z)

n, d = G_z.ratsimp().as_numer_denom()
n, d = n.as_poly(z), d.as_poly(z)

display(n, d)

We first define the symbols needed. Then set up G_z according to our equation above. After displaying it for good measure, we do a rational simplification ratsimp() on it and split it into numerator n and denominator d.

Now, we represent n and d as polynomials of z. And this is how we obtain our coefficients a1, a2, and so on. We continue:

scale = d.LC()

n, d = (n.as_expr()/scale).as_poly(z), (d.as_expr()/scale).as_poly(z)
N, D = n.degree(), d.degree()

We find the scaling factor as the leading coefficient (LC()) of the denominator d, scale both accordingly and determine the degree of each (this is mostly to determine for how many terms we have to iterate when printing etc.).

Bs = n.all_coeffs()
As = d.all_coeffs()

display(Bs)
display(As)

Now, basically Bs contains all coefficients [b0 b1] and likewise for As. Here’s a shortcut, see if you can spot it.

And now comes the cse magic:

commons, exprs = sy.cse(Bs)

This returns a tuple. The first element is a list of tuples of arbitrarily assigned place holder variables x0, x1, and so on. The second element are the original expressions, but the placeholders are substituted. In this example:

(
  [
    (x0, K_i*T_s/2)
  ],
  [K_p + x0, -K_p + x0]
)

Which means, instead of calculating K_i*T_s/2 twice, we can do it once, assign it to x0 and use that on the expressions.

You can even concatenate As and Bs for good measure (not really necessary here, since the As only contain 1 and 0, but in other controller topologies this may change).

commons, exprs = sy.cse(As + Bs)

# Printing C code for instance
for com in commons:
    print( sy.ccode(com[1], assign_to=repr(com[0])) )

for expr, name in zip(exprs, "a0 a1 b0 b1".split(" ")):
    print( sy.ccode(expr, assign_to=f"coeffs.{name:s}") )

We obtain the common terms commons and the reduced expressions exprs, then print the common terms assigning them to their respective placeholder variables, and finally print the reduced expressions. This produces this five liner in C:

x0 = (1.0/2.0)*K_i*T_s;
coeffs.a0 = 1;
coeffs.a1 = 0;
coeffs.b0 = K_p + x0;
coeffs.b1 = -K_p + x0;

Now, how is this useful you ask? Well, for this toy example it certainly is a lot of work for five lines I could’ve just typed out.

But if you go into a PID controller with damped differential term, and you want your code to be efficient (i.e. have as few operations as possible), this is very handy. Just by changing the symbols and the expression for G_z we obtain this:

x0 = 2*T_v;
x1 = 1.0/(T_s + x0);
x2 = 4*T_v;
x3 = T_i*T_s;
x4 = 2*x3;
x5 = T_i*x2;
x6 = 1.0/(x4 + x5);
x7 = K_p*x4;
x8 = K_p*T_s*x0;
x9 = pow(T_s, 2);
x10 = 4*K_p*T_d*T_i + K_p*x5;
x11 = K_p*x9 + x10;
c.a0 = 1;
c.a1 = -x1*x2;
c.a2 = x1*(-T_s + 2*T_v);
c.b0 = x6*(x11 + x7 + x8);
c.b1 = (K_p*x9 - x10)/(T_i*x0 + x3);
c.b2 = x6*(x11 - x7 - x8);

I did this manually once, and only obtained five or six of the placeholders. That’s the power of cse().

Oh, you need the code in Python for testing? No problem, just use pycode instead (and work around that assign_to parameter).

So, there you have it.