Runtime Optimization of sympy code using numpy or scipy

Question

I'm new to programming and it's my first question here. I spend a lot of time on this problem. It's an abstract from a programme I'm writing for my thesis. Maybe it's trivial to some of you experts but I might even lack mathematical knowledge.

I want to solve following set of equations for b1, b2 and c3:

0 = cos(b1)+ cos(b2)- 0.0166
0 = sin(b1)+ sin(b2)+ 0.3077*c3 - 0.6278
0 = cos(b1)- cos(b2)+ 5.4155*c3 - 4.3547

b1 and b2 are angles and should therefore between [0, 2*pi] and c3 should be between [0,1], but this is not necessary. I can filter results later.

I found a solution using sympy:

import sympy as sy
b1 = sy.Symbol('b1',real=True)
b2 = sy.Symbol('b2',real=True)
c3 = sy.Symbol('c3',real=True)
a = sy.cos(b1)+sy.cos(b2)-0.0166
b = sy.sin(b1)+sy.sin(b2)+0.3077*c3-0.6278
c = sy.cos(b1)-sy.cos(b2)+5.4155*c3-4.3547
Erg = sy.solve([a,b,c,],[b1,b2,c3],dict=True)

[{c3: -0.4634, b1: 2.7245, b2: 0.3739}]

The results are reasonable, but unfortunately this leads to a calculation time of 6s which would result in a total calculation time of over 7hrs for my programme. Please help me find a faster solution. I tried scipy

from scipy.optimize import fsolve
import numpy as np

def equations(p):
    b1, b2, c3 = p
    return (np.cos(b1)+np.cos(b2)-0.0166, np.sin(b1)+np.sin(b2)+0.3077*c3-0.6278,np.cos(b1)-np.cos(b2)+5.4155*c3-4.3547)

b1, b2, c3 =  fsolve(equations, (-0.4634, 2.7245, 0.3739))

This results in under 1s in (-9.8418e-14, 5.6621e-15,-7.5495e-14). I don't know which number belongs to which variable, but they don't make any sense anyway. Another option would be optimizing the equation for minimal c3. Please don't hesitate to ask if I haven't been clear about anything.

Answer 1

If I fix the missing parenthesis in your numpy example, I get a different result:

In [40]: b1,b2,c3
Out[40]: (0.21407625679722384, 2.8598524043642226, 0.4463029985826017)

And anyway your input points are not really good, since your variables are mixed up. Your function is defined in the order b1,b2,c3, yet the starting point corresponds to the symbolic solution of c3,b1,b2. Changing the order gives me

In [42]: b1,b2,c3
Out[42]: (2.9450733464142882, 0.07278767695320594, 1.1693849677611501)

Note that for a given b1,b2 you'll also have a solution with b1+2*k*pi,b2*2*l*pi where k,l are integers.

I didn't put these numbers into your equations, but I'm pretty sure they solve them. The problem is, there are a bunch of solutions. As I said, your first two variables are periodic, and even then there can be a lot of other solutions with c3. If you use a bunch of different starting points for fsolve, you'll get a bunch of different solutions.

What's worse: your solution doesn't enforce the 0<c3<1 bound.

---

On the mathematics: due to the diversity of possible solutions, you might want to simplify the numerical problem on paper first.

For instance, the first equation is special, as it doesn't contain c3. You can relate b1 to b2, albeit in a nonlinear way. For each b1 you'll know that

b2=+-acos(0.0166-cos(b1))+2*k*pi

Which is already a constraint. Furthermore, adding Eq. 1 to Eq. 3 gives you

0 = 2*cos(b1)+ 5.4155*c3 - 4.3547-0.0166

in other words, c3 is linearly dependent on cos(b1). It might help if you introduce this into Eq. 2, getting a connection between functions of b1 and b2. Essentially you have variables as b1,b2(b1),c3(b1), with the single Eq. 2 to solve.

Another note: I'm guessing that the large number of equations comes from a lot of combinations for numerical factors in your equations. Isn't that right? If so, you might want to introduce some symbolic constants in the equations instead of the floating-point literals, and try to solve that. If you succeed, you don't have to solve a huge number of sets of nonlinear equations: you only have to solve one set, then substitute your various parameters.