from gekko import GEKKO

m = GEKKO()

#%% Constants
pi = m.Const(3.14159,'pi')
P = 2300        # compressive load (kg_f)
o_y = 450       # yield stress (kg_f/cm^2)
E = 0.65e6      # elasticity (kg_f/cm^2)
p = 0.0020      # weight density (kg_f/cm^3)
l = 300         # length of the column (cm)

#%% Variables
d = m.Var(value=8.0,lb=2.0,ub=14.0) # mean diameter (cm)
t = m.Var(value=0.3,lb=0.2,ub=0.8)  # thickness (cm)
cost = m.Var()

#%% Intermediates
d_i = m.Intermediate(d - t)
d_o = m.Intermediate(d + t)
W = m.Intermediate(p*l*pi*(d_o**2 - d_i**2)/4) # weight (kgf)
o_i = m.Intermediate(P/(pi*d*t))    # induced stress

# second moment of area of the cross section of the column
I = m.Intermediate((pi/64)*(d_o**4 - d_i**4))

# buckling stress (Euler buckling load/cross-sectional area)
o_b = m.Intermediate((pi**2*E*I/l**2)*(1/(pi*d*t)))

#%% Equations
m.Equations([
        o_i - o_y <= 0,
        o_i - o_b <= 0,
        cost == 5*W + 2*d
        ])

#%% Objective
m.Minimize(cost)

#%% Solve and print solution
m.options.SOLVER = 1
m.solve()

print('Optimal cost: ' + str(cost[0]))
print('Optimal mean diameter: ' + str(d[0]))
print('Optimal thickness: ' + str(t[0]))