import numpy as np
from gekko import GEKKO

class model():
    def __init__(self, ab, cd, ef):
        self.ab = ab
        self.cd = cd
        self.ef = ef

    """Objective Function"""
    'Entry level'
    def fun(self,x):
        paramA = self.paramA()
        return paramA + ((np.pi * x[0] * x[1] * x[2]) * 1e-4)**0.81

    def paramA(self):
        if self.ab >20:
            results = self.ab + self.cd ** self.ef
            return results
        else:
            results = self.ab*2 - self.cd ** self.ef
            return results

'Process inputs'
## Here I just define the initial process values. It will be used for pDrop and Q calculations
obj = model(360, 1.20, 300)

m = GEKKO()
x = []
x1 = m.Var(value=20,lb=20, ub=6555)  #integer=True
x2 = m.Var(value=1,lb=1,ub=10000)  #integer=True
x3 = m.sos1([30, 42, 45, 55])
x3.value = 1.0
x = [x1, x2, x3]


m.Equation((x1 * x2* x3) * 1e-6 >= 50)


m.Obj(obj.fun(x))

# Change to True to initialize with IPOPT
init = False
if init:
    m.options.SOLVER=3
    m.solve(disp=False) # Solve

m.options.SOLVER=1
m.solve(disp=True) # Solve

print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('Objective: ' + str(m.options.objfcnval))

x1 = m.Var(value=20,lb=20, ub=6555)  #integer=True
x2 = m.Var(value=1,lb=1,ub=10000)  #integer=True
x3 = m.sos1([30, 42, 45, 55])