# Crane Hook optimization
from gekko import GEKKO
from numpy import pi

m = GEKKO(remote=False)

# Constants
F = 100                 # load (N)
S_y = 430               # yield strength of steel (N/mm^2)
                        #  Grade D fine-carbon steel (ASTM A255)         

# Variables
r_o  = m.Var()          # outer radius (mm)
r_i  = m.Var(lb=1.5)    # inner radius (mm)
b    = m.Var(lb=0.2)    # hook width (mm)
V    = m.Var()          # hook volume (mm^3)

# Intermediates
h    = r_o - r_i        # hook height (mm)
R    = (r_o + r_i)/2    # radius of the centroid (mm) 
r_n  = h/m.log(r_o/r_i) # radius of the neutral axis (mm)
e    = R - r_n          # R - r_n (mm)
M    = F * R            # bending moment due to the load
c_o  = r_o - r_n        # distance from outer to neutral(mm)
c_i  = r_n - r_i        # distance from inner to neutral(mm)
Area = b*h              # cross-sectional area (mm^2)
o_A = (M*c_o/(Area*e*r_o)) # outer stress
o_B = (M*c_i/(Area*e*r_i)) # inner stress

# Equations
m.Equations([
        V == pi*(r_o**2-r_i**2)*b, # volume calculation
        o_A < S_y,      # yield stress @ A < yield strength
        o_B < S_y,      # yield stress @ B < yield strength
        r_i < r_o       # constraint for feasibility
        ])

# Objective
m.Minimize(V)

# Solve
m.options.SOLVER = 3
m.solve()

print('Optimal Volume: ' + str(V[0]))
print('Optimal outer radius: ' + str(r_o[0]))
print('Optimal inner radius: ' + str(r_i[0]))
print('Optimal hook width: ' + str(b[0]))