# 4-Bar Mechanism Generation

For a school project, I will design and prototype a bicycle brake that uses a four-bar linkage to accomplish its goal.

My Python 3 code does this mechanism generation, and implements methods to remove mechanisms I don't think will work well to hopefully give me a good final mechanism to construct. Overall I wrote this up quite quickly, hence its messiness. I'm hoping to dramatically clean up the code and improve efficiency. Any suggestions on where to begin?

## Background

(This section isn't critical to review the code, but I thought would be nice to include)

To create the mechanism, I need to define a set of three "precision positions" on the coupler curve, which are points the coupler point of the four-bar will pass through.

Since the points I define are just approximate guesses what I think the path should be (approximate linear motion), I loop over a small range for each coordinate I define.

I then also choose the driver rotation between PP1 and PP2 ($\beta_2$) and PP1 and PP3 ($\beta_3$). The $\alpha$'s are solved from the given precision positions (the couplers rotation), as well as the $\delta$'s (the translational movement of the coupler point). The following picture illustrates this:

The following system of equations is solved to obtain a 2 dyads (shown above, cut the mechanism in half between the coupler point, the two links on each side is a dyad). The two dyads together form the 4-bar.

$$\vec{W_A}(e^{i\beta_{2A}} - 1) + \vec{Z_A}(e^{i\alpha_{2A}} - 1) = \delta_2$$ $$\vec{W_A}(e^{i\beta_{3A}} - 1) + \vec{Z_A}(e^{i\alpha_{3A}} - 1) = \delta_3$$

$$\vec{W_B}(e^{i\beta_{2B}} - 1) + \vec{Z_B}(e^{i\beta_{2A}} - 1) = \delta_2$$ $$\vec{W_B}(e^{i\beta_{3B}} - 1) + \vec{Z_B}(e^{i\beta_{3A}} - 1) = \delta_3$$

gen.py:

import numpy as np
import matplotlib.pyplot as plt
import math
from plot import plot_four_bar_list

def check_lengths(four_bar_set):
if (not (np.abs(four_bar_set[0][1]) - np.abs(four_bar_set[1][1])) < 1e-10 or
not (np.abs(four_bar_set[0][2]) - np.abs(four_bar_set[1][2])) < 1e-10 or
not (np.abs(four_bar_set[0][3]) - np.abs(four_bar_set[1][3])) < 1e-10 or
not (np.abs(four_bar_set[0][4]) - np.abs(four_bar_set[1][4])) < 1e-10 or
not (np.abs(four_bar_set[0][1]) - np.abs(four_bar_set[2][1])) < 1e-10 or
not (np.abs(four_bar_set[0][2]) - np.abs(four_bar_set[2][2])) < 1e-10 or
not (np.abs(four_bar_set[0][3]) - np.abs(four_bar_set[2][3])) < 1e-10 or
not (np.abs(four_bar_set[0][4]) - np.abs(four_bar_set[2][4])) < 1e-10):
print('ERROR!')
return

def main():
# precision positions for coupler to pass through
pp1 = (0, 0)
pp2 = (-5, 0)
pp3 = (-10, 0)
pp = [pp1, pp2, pp3]
solutions = []
filtered_solutions = []
delta = []

# calculate changes between precesion positions
for i in range(0, 3):
delta.append(complex(pp[i][0] - pp[0][0], pp[i][1] - pp[0][1]))

# generate set of solutions
for y1 in range(-10, 10, 5):
for y2 in range(-10, 10, 5):
for y3 in range(-10, 10, 5):
delta2 = (complex(pp[1][0] - pp[0][0], y2 - y1))
delta3 = (complex(pp[2][0] - pp[0][0], y3 - y1))
for theta1 in range(-10, 10, 5):
for theta2 in range(-10, 10, 5):
for theta3 in range(-10, 10, 5):
alpha2 = theta2 - theta1
alpha3 = theta3 - theta1
for i in range(0, 120, 40): # beta_2a in degrees
for j in range(0, 120, 40): # beta_3a in degrees
for k in range(0, 360, 120): # beta_2b in degrees
for l in range(0, 360, 120): # beta_3b in degrees
dyad_a = np.matrix([[np.expm1(complex(0, i)), np.expm1(complex(0, alpha2))],
[np.expm1(complex(0, j)), np.expm1(complex(0, alpha3))] ])
dyad_b = np.matrix([[np.expm1(complex(0, k)), np.expm1(complex(0, i))],
[np.expm1(complex(0, l)), np.expm1(complex(0, j))] ])
delta_mat = [delta2, delta3]
try:
ground = [complex(pp1[0], pp1[1]) - part_a[0] - part_a[1], complex(pp1[0], pp1[1]) - part_b[0] - part_b[1]]
# stored in list as ground_a, wa, za, wb, zb, ground_b; then changes to those
four_bar_list = [[ground[0], part_a[0], part_a[1], part_b[0], part_b[1], ground[1]],
[ground[0], part_a[0] * np.exp(complex(0, i)), part_a[1] * np.exp(complex(0, alpha2)), part_b[0] * np.exp(complex(0, k)), part_b[1] * np.exp(complex(0, i)), ground[1]],
[ground[0], part_a[0] * np.exp(complex(0, j)), part_a[1] * np.exp(complex(0, alpha3)), part_b[0] * np.exp(complex(0, l)), part_b[1] * np.exp(complex(0, j)), ground[1]]]
check_lengths(four_bar_list)
solutions.append(four_bar_list)
except:
pass

# filter solutions
print('Number of solutions before filtering: {}'.format(len(solutions)))
for four_bar_list in solutions:
filter_length(solutions, filtered_solutions, four_bar_list)
solutions = filtered_solutions
filtered_solutions = []
print('Number of solutions after filtering based on length: {}'.format(len(solutions)))

for four_bar_list in solutions:
filter_orientation(solutions, filtered_solutions, four_bar_list)
solutions = filtered_solutions
filtered_solutions = []
print('Number of solutions after filtering based on orientation: {}'.format(len(solutions)))

for four_bar_list in solutions:
filter_configuration(solutions, filtered_solutions, four_bar_list)
solutions = filtered_solutions
filtered_solutions = []
print('Number of solutions after filtering based on configurations: {}'.format(len(solutions)))

# for four_bar_list in solutions:
#     plot_four_bar_list(four_bar_list)
if(len(solutions) > 0):
plot_four_bar_list(solutions[0])

if __name__ == "__main__":
main()


filter.py:

import numpy as np
import math

ground_a = (four_bar[0].real, four_bar[0].imag)
joint_a = (ground_a[0] + four_bar[1].real, ground_a[1] + four_bar[1].imag)
joint_c = (joint_a[0] + four_bar[2].real, joint_a[1] + four_bar[2].imag)
joint_b = (joint_c[0] - four_bar[4].real, joint_c[1] - four_bar[4].imag)
ground_b = (joint_b[0] - four_bar[3].real, joint_b[1] - four_bar[3].imag)
joints = [ground_a, joint_a, joint_c, joint_b, ground_b]
return joints

def filter_length(solutions, filtered_solutions, four_bar_list):
# remove any solutions with links greater than 50 mm or less than 10mm
return
return
else:
filtered_solutions.append(four_bar_list)

def filter_orientation(solutions, filtered_solutions, four_bar_list):
for four_bar in four_bar_list:

# only be on one side of the bike
for joint in joints:
if joint[0] < joints[0][0]:
return

# grounds should be around same y-axis
if joints[0][0] - joints[4][0] > 5:
return

# grounds should be around same x-axis
elif joints[0][1] - joints[4][1] > 5:
return

# coupler should be below first ground
elif joints[2][1].item() > joints[0][1].item():
return
else:
filtered_solutions.append(four_bar_list)

def filter_configuration(solutions, filtered_solutions, four_bar_list):
# check circuit change
mu = []
for four_bar in four_bar_list:
r7 = (four_bar[0] + four_bar[1]) - (four_bar[0] + four_bar[5])
r4 = (four_bar[5] + four_bar[4])
psi = np.angle(r4, deg=True) - np.angle(r7, deg=True)
if psi < -180:
psi += 180
elif psi > 180:
psi -= 180
mu.append(psi >= 0)

if len(set(mu)) <= 1:
filtered_solutions.append(four_bar_list)


plot.py:

import matplotlib.pyplot as plt

def plot_four_bar_list(four_bar_list):
x1 = []
y1 = []
x2 = []
y2 = []
x3 = []
y3 = []
x1.append(joint[0])
y1.append(joint[1])
x2.append(joint[0])
y2.append(joint[1])
x3.append(joint[0])
y3.append(joint[1])

fig, ax = plt.subplots(3, 1, sharex='all', sharey='all')
ax[0].plot(x1, y1, 'r-o', label='config_1')
ax[1].plot(x2, y2, 'b-o', label='config_2')
ax[2].plot(x3, y3, 'g-o', label='config_3')
fig.suptitle('4-Bar Setup Between Precision Points')
plt.show()

def show_rotation():
return


## check_lengths

### Avoid repetition

When the same line of code is repeated many many times the information gets lost in the noise and it is not obvious at first glance:

not (np.abs(four_bar_set[a][b]) - np.abs(four_bar_set[c][d])) < 1e-10
not (np.abs(four_bar_set[a][b]) - np.abs(four_bar_set[c][d])) < 1e-10
...


or

In this lines of code only a,b,c and d change, all the rest is equal, how can we do better?

The solution is to make a data structure containing the variable data:

values  = [ # <- Give better name
(0, 1, 1, 1),
(0, 2, 1, 3),
...
]


And the using the any built-in that is the same as reduceing or over a list.

if (any( not (np.abs(four_bar_set[a][b]) - np.abs(four_bar_set[c][d])) < 1e-10 ) for a, b, c, d in values)):
# Handle


We could even go a step further and define an helper function:

 def four_bar_set_error_inside_limit(limit, a, b, c, d):
return (np.abs(four_bar_set[a][b]) - np.abs(four_bar_set[c][d])) < 1e-10


and write:

if (all( four_bar_set_error_inside_limit(limit, a, b, c, d) ) for a,b,c,d in values):
# Handle


Please note that this way the code is immediately clear, also because I removed the double negation of any and not and used just all.

### Avoid printing from inside a function

    print('ERROR!')
return


This way of handling things is one-off, non-standard and hard to reuse. The caller function should call this function to know if a certain condition is met (which condition exactly? Without docstring and with such a general name the user of this function is forced to read the code to understand, that should not happen), then you should let him free to handle this result however he likes, maybe logging a message, maybe crashing, maybe trying another way of calculation.

So in my opinion you should just return the value, so:

### Final version

def four_bar_set_error_inside_limit(tolerance, a, b, c, d):
return (np.abs(four_bar_set[a][b]) - np.abs(four_bar_set[c][d])) < tolerance

def all_lengths_deltas_below_error(four_bar_set, tolerance):
values = [...]
return all(four_bar_set_error_inside_limit(tolerance, a, b, c, d)
for a, b, c, d in values)


Also note that this final version has a tolerance parameter that you can easily pass from the outside while the previous version had the value 1e-10 repeated 8 times and changing it would be time consuming and prone to error forgetting one.

A final improvement would be to add a doctest and a small docstring explaining exactly what a four_bar_set is in real life and how it is represented, and why it is important to check that this property holds true.

## filter_length

### Avoid an out parameter

def filter_length(solutions, filtered_solutions, four_bar_list):
# remove any solutions with links greater than 50 mm or less than 10mm
return
return
else:
filtered_solutions.append(four_bar_list)


In this function filtered_solutions is an out parameter, that is a parameter that is passed into a function for the sole purpose of being modified to contain the output of the function.

This is non-standard in modern Python and there is a much simpler way to handle this, the yield keyword:

def filter_length(solutions, filtered_solutions, four_bar_list):
# remove any solutions with links greater than 50 mm or less than 10mm
if not (20 < np.abs(link) < 80):
return
yield four_bar_list


        # remove any solutions with links greater than 50 mm or less than 10mm