from gekko import GEKKO
m = GEKKO()
Y = {0: 10, 1: 20, 2: 15, 3: 40, 4: 8}
Z = {0: 20, 1: 30, 2: 10, 3: 25, 4: 18}

x = m.Array(m.Var,(5,5),value=0,lb=0,ub=1,integer=True)
for i in range(5):
    # 5 time slots, 5 order options for machine Y and Z
    m.Equation(m.sum([x[i,j] for j in range(5)])==1)
    m.Equation(m.sum([x[j,i] for j in range(5)])==1)

# time to process on Y and Z
ty = m.Array(m.Var,5)
tz = m.Array(m.Var,5)
for i in range(5):
    # time to process on Y
    m.Equation(ty[i] == m.sum([x[i,j]*Y[j] for j in range(5)]))
    # time to process on Z
    m.Equation(tz[i] == m.sum([x[i,j]*Z[j] for j in range(5)]))

# delay time on Z
d = m.Array(m.Var,5,lb=0)
s = m.Array(m.Var,5,lb=0)
m.Equation(d[0]==ty[0])
m.Equation(s[0]==0)
for i in range(1,5):
    m.Equation(d[i]>=tz[i-1]-ty[i]+s[i])
m.Minimize(m.sum(d))
m.Minimize(1e-3*m.sum(s))

m.options.SOLVER=1
m.solve()

print('Order of processing, rows=order, columns=product')
print(x)
print('time on machine Y')
print(ty)
print('time on machine Z')
print(tz)
print('delay')
print(d)