|
#!/usr/bin/env python |
|
|
|
import numpy as np |
|
from collections import OrderedDict |
|
from scipy.optimize import LinearConstraint, minimize, BFGS |
|
|
|
|
|
def normlen(v): |
|
''' length of norm of a vector ''' |
|
return np.sqrt(np.sum(np.square(v))) |
|
|
|
|
|
# FREE PARAMETERS |
|
L = 20 # Length of upper triangle side |
|
H = 9 # Height of structure |
|
G = 12 # Length of lower triangle side |
|
T = 0 # Angle (degrees) of top triangle (should probably be 0) |
|
U = 200 # Angle (degrees) of bottom triangle (should probably be 200) |
|
|
|
# GEOMETRY CALCULATION |
|
X, Y, Z = np.eye(3) # Unit vectors |
|
J = L * np.sqrt(3) / 3 # Radius of upper triangle |
|
K = G * np.sqrt(3) / 3 # Radius of lower triangle |
|
# Top triangle points |
|
a, b, c = [np.deg2rad(T + x) for x in range(0, 360, 120)] # Angles in radians |
|
A, B, C = [np.array([np.sin(x) * J, np.cos(x) * J, H]) for x in (a, b, c)] |
|
# Bottom triangle points |
|
d, e, f = [np.deg2rad(U + x) for x in range(0, 360, 120)] # Angles in radians |
|
D, E, F = [np.array([np.sin(x) * K, np.cos(x) * K, 0]) for x in (d, e, f)] |
|
|
|
|
|
def normvec(v): |
|
''' return a normalized vector ''' |
|
return np.array(v) / normlen(v) |
|
|
|
|
|
# structural force vectors |
|
structure = OrderedDict([ |
|
('rod AD', normvec(A - D)), |
|
('rope AB', normvec(B - A)), |
|
('rope AC', normvec(C - A)), |
|
('rope AE', normvec(E - A)), |
|
]) |
|
SV = np.array(list(structure.values())).transpose() |
|
|
|
|
|
# hammock force vectors |
|
hammocks = OrderedDict([ |
|
('hammock AB', normvec(B - A - Z * L / 2)), |
|
('hammock AC', normvec(C - A - Z * L / 2)), |
|
]) |
|
HV = np.array(list(hammocks.values())).transpose() |
|
|
|
|
|
def error(hf, sf): |
|
''' Calculate error in total forces ''' |
|
return normlen(np.dot(SV, sf) + np.dot(HV, hf)) |
|
|
|
|
|
def solve_sf(hf): |
|
''' Solve for structural forces (sf) given hammock forces (hf) ''' |
|
# Convert to array |
|
hf = np.array(hf) |
|
# Total forces at corner must equal zero |
|
# SV * sf + HV * hf = 0, therefore: -HV * hf <= SV * sf <= -HV * hf |
|
equality_constraint = LinearConstraint(SV, -np.dot(HV, hf), -np.dot(HV, hf)) |
|
# Also forces must all be nonnegative |
|
n = len(structure) |
|
positive_constraint = LinearConstraint(np.eye(n), np.zeros(n), np.ones(n) * np.inf) |
|
# We will use both constraints for optimization |
|
constraints = [equality_constraint, positive_constraint] |
|
|
|
# Calculate structure forces as minimum total forces which satisfy constraints |
|
sf = minimize(np.sum, # function to minimize -- want the smallest sum of total forces |
|
np.zeros(n), # initial guess |
|
method='trust-constr', # For when we have LinearConstraint's |
|
constraints=constraints, # Constraints on our solution space |
|
# These two are defaults, but need to be included because of a bug. |
|
jac='2-point', hess=BFGS()) # https://github.com/scipy/scipy/issues/8867 |
|
return sf.x |
|
|
|
|
|
def main(): |
|
# Check geometry |
|
print('Top triangle sides', normlen(A - B), normlen(B - C), normlen(C - A)) |
|
print('Bottom triangle sides', normlen(D - E), normlen(E - F), normlen(F - D)) |
|
print('Rod lengths', normlen(A - D), normlen(B - E), normlen(C - F)) |
|
# Print out structure forces for different hammock loads |
|
for hf in [[200, 200], [200, 0], [0, 200]]: |
|
sf = solve_sf(hf) |
|
print('hammock forces:', list(zip(hammocks.keys(), hf))) |
|
print('error', error(hf, sf)) |
|
print('structure forces:') |
|
for name, force in zip(structure.keys(), sf): |
|
print(f' {name} {force:.2f}') |
|
|
|
|
|
if __name__ == '__main__': |
|
main() |