from gekko import GEKKO
from numpy import pi

m = GEKKO(remote=False)

# Constants
Px = 10              # Transverse load (lb)      
Py = 25              # Axial load (lb)
Ln = 50              # Beam length (in)
rho = 0.3            # Material density (lb/in^3)
S_y = 30000          # Yield stress (psi)
E = 30e6             # Young's modulus

# Variables
b = m.Var(lb=0.5)       # Beam width (in)
d = m.Var(lb=1e-3)      # Beam depth (in)
weight = m.Var(lb=1e-3) # Beam weight (lb)

# Intermediate calculations
I_zz = b*d**3/12         # Moment of Inertia 
x_stress = m.Var(ub=S_y) # Compressive stress due to Px
y_stress = m.Var(ub=S_y) # Compressive stress due to Py
# Intermediate: Axial Buckling Load
Py_ab = m.Intermediate(pi**2*E*I_zz/(4*Ln**2)) 

# Equations
m.Equations([
        weight == b*d*Ln*rho,
        y_stress*(b*d)==Py,
        x_stress*(2*I_zz)==Px*Ln*d,
        Py_ab >= Py,
        b <= 2*d
        ])

# Objective
m.Minimize(weight)

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

print('Weight: ' + str(weight[0]))
print('Depth: ' + str(d[0]))
print('Width: ' + str(b[0]))
print('Vertical stress: ' + str(x_stress[0]))
print('Horizontal stress: ' + str(y_stress[0]))
print('Axial Buckling Load: ' + str(Py_ab[0]))