from gekko import GEKKO
m = GEKKO()

# revenue per dv
r_dv = 550  
# wet mass  (kg)
m_0 = m.Var(value=150650)                   
# dry mass. (kg)
m_f = m.Var(lb=20000)                       
# effective exhaust velocity. m/s
v_0 = m.Var(value=3500,lb=2500,ub=4500)     
# velocity of the vehicle to orbital velocity at 242 km
dv = m.Var(lb=9400, ub=20200)               
profit = m.Var() # revenue - cost

# Intermediates
c_fuel = m.Intermediate(4.154*(m_0 - m_f)) # Cost of fuel
c_dry = m.Intermediate(154.36*m_f) # Cost of rocket
c_exhaust = m.Intermediate(75*v_0) # Cost of exhaust velocity
cost = m.Intermediate(c_fuel + c_dry + c_exhaust) 
revenue = m.Intermediate(dv*r_dv) # Revenue

# Equations
m.Equations([
        dv == v_0*m.log((m_0/m_f)),
        m_0 >= 2*m_f,
        profit == revenue - cost
        ])

# Objective
m.Maximize(profit)
m.options.SOLVER = 3
m.solve()

print('wet mass: ' , str(m_0[0]))
print('dry mass: ' , str(m_f[0]))
print('dv: ' + str(dv[0]))
print('v_0: ' + str(v_0[0]))