import numpy as np
from gekko import GEKKO
m = GEKKO(remote=False)
# input  (1) = u
# states (5) = x[1] to x[5]
# output (1) = y

# Constants
n = 5
K = m.Const(4)

# Parameter
u = m.Param(value=3,lb=0,ub=10)

# Variable array
x = m.Array(m.Var,5,lb=0)

# Intermediate
y = m.Intermediate(x[4])

# Equations
m.Equation( x[0].dt()  + x[0] == K*u)
m.Equations([x[i].dt() + x[i] == x[i-1]
             for i in range(1,n)])

# Objective
m.Minimize((y-5)**2)

# Define time
m.time = np.linspace(0,20)

# Solve dynamic simulation
m.options.IMODE = 4
m.solve()

# Show solution
import matplotlib.pyplot as plt
for i in range(n):
    plt.plot(m.time,x[i],label=f'$x_{i}$')
plt.plot(m.time,y,'k--',label='y')
plt.xlabel('Time'); plt.grid(); plt.legend()
plt.show()