4-Bar Mechanism Generation

Question

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 filter import links_to_joints, filter_length, filter_orientation, filter_configuration, filter_mechanical_advantage
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):
                    theta1 = math.radians(theta1)
                    for theta2 in range(-10, 10, 5):
                        theta2 = math.radians(theta2)
                        for theta3 in range(-10, 10, 5):
                            theta3 = math.radians(theta1)
                            alpha2 = theta2 - theta1
                            alpha3 = theta3 - theta1
                            for i in range(0, 120, 40): # beta_2a in degrees
                                i = math.radians(i)
                                for j in range(0, 120, 40): # beta_3a in degrees
                                    j = math.radians(j)
                                    for k in range(0, 360, 120): # beta_2b in degrees
                                        k = math.radians(k)
                                        for l in range(0, 360, 120): # beta_3b in degrees
                                            l = math.radians(l)
                                            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:
                                                part_a = np.linalg.solve(dyad_a, delta_mat)
                                                part_b = np.linalg.solve(dyad_b, delta_mat)
                                                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

def links_to_joints(four_bar):
    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):
    for link in four_bar_list[0]:
        # remove any solutions with links greater than 50 mm or less than 10mm
        if np.abs(link) < 20:
            return
        elif np.abs(link) > 80:
            return
        else:
            filtered_solutions.append(four_bar_list)

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

        # 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
from filter import links_to_joints

def plot_four_bar_list(four_bar_list):
    x1 = []
    y1 = []
    x2 = []
    y2 = []
    x3 = []
    y3 = []
    for joint in links_to_joints(four_bar_list[0]):
        x1.append(joint[0])
        y1.append(joint[1])
    for joint in links_to_joints(four_bar_list[1]):
        x2.append(joint[0])
        y2.append(joint[1])
    for joint in links_to_joints(four_bar_list[2]):
        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

Answer 1

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):
    for link in four_bar_list[0]:
        # remove any solutions with links greater than 50 mm or less than 10mm
        if np.abs(link) < 20:
            return
        elif np.abs(link) > 80:
            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):
    for link in four_bar_list[0]:
        # 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

Wrong comments are worse than no comments

        # remove any solutions with links greater than 50 mm or less than 10mm
        if np.abs(link) < 20:
            return
        elif np.abs(link) > 80:
            return

The comment states 50 and 10 while the code states 20 and 80 so something is clearly wrong.

Maybe you are using different units of measure between the comment and the code, or maybe you changed your mind halfway between writing your comment and the code, or maybe you changed the code and forgotten the comment.

Anyhow this is confusing, also because no link to an original source or material is given in order to discern which between the comment and the code is correct.

This is another example of why comments stating the obvious and repeating the code should not be written in the first place (https://en.wikipedia.org/wiki/Single_source_of_truth and https://en.wikipedia.org/wiki/Dont_repeat_yourself principles).