from gekko import GEKKO

# Initialize Gekko model
m = GEKKO()

#Maximize force of a spring at its preload height h_0 of 1 inches
#The stress at Hs (solid height) must be less than Sy to protect from damage

# Constants
from numpy import pi

# Model Parameters
delta_0 = 0.4 # inches (spring deflection)
h_0 = 1.0 # inches (preload height)
Q = 15e4 # psi
G = 12e6 # psi
S_e = 45e3 # psi
S_f = 1.5 
w = 0.18


# Variables
# inches (wire diameter)
d_wire = m.Var(value=0.07247, lb = 0.01, ub = 0.2) 
# inches (coil diameter)
d_coil = m.Var(value=0.6775, lb = 0.0) 
# number of coils in the spring
n_coils = m.Var(value=7.58898, lb = 0.0) 
# inches (free height spring exerting no force)
h_f = m.Var(value=1.368117, lb = 1.0) 
F = m.Var() # Spring force


# Intermediates
S_y = m.Intermediate((0.44 * Q) / (d_wire**w))
h_s = m.Intermediate(n_coils * d_wire)
k = m.Intermediate((G * d_wire**4)/(8 * d_coil**3 * n_coils))
kW = m.Intermediate((4 * d_coil - d_wire)/(4 * (d_coil-d_wire)) \
                    + 0.62 * d_wire/d_coil)
n = m.Intermediate((8 * kW * d_coil)/(pi * d_wire**3))
tau_max = m.Intermediate((h_f - h_0 + delta_0) * k * n)
tau_min = m.Intermediate((h_f - h_0) * k * n)
tau_mean = m.Intermediate((tau_max + tau_min) / 2)
tau_alt = m.Intermediate((tau_max - tau_min) / 2)
h_def = m.Intermediate(h_0 - delta_0)
tau_hs = m.Intermediate((h_f - h_s) * k * n)

# Equations
m.Equations([
    F == k * (h_f - h_0),
    d_coil / d_wire >= 4,
    d_coil / d_wire <= 16,
    d_coil + d_wire < 0.75,
    h_def - h_s > 0.05,
    tau_alt < S_e / S_f,
    tau_alt + tau_mean < S_y / S_f,
    tau_hs < S_y / S_f
    ])

# Objective function
m.Obj(-F)

# Send to solver
m.solve()

# Print solution
print('Maximum force: ' + str(F[0]))
print('Optimal values: ')
print('d_wire : ' + str(d_wire[0]))
print('d_coil : ' + str(d_coil[0]))
print('n_coils : ' + str(n_coils[0]))
print('h_f : ' + str(h_f[0]))