Question: I am starting to learn Gekko and I am testing optimal control problems. I am trying to solve the following optimal control problem with Gekko  The solution of this problem is (x_1(t) = (t-2)^2 - 2) How to build the constraint x(0) + x(2) = 0? My code gives me a wrong solution. m = GEKKO(remote=False) # initialize gekko nt = 101 m.time = np.linspace(0,2,nt) #end_loc = nt-1 # Variables x1 = m.CV(fixed_initial=False) x2 = m.CV(fixed_initial=False) x3 = m.Var(value=0) #u = m.Var(value=0,lb=-2,ub=2) u = m.MV(fixed_initial=False,lb=-2,ub=2) u.STATUS = 1 p = np.zeros(nt) # mark final time point p[-1] = 1.0 final = m.Param(value=p) p1 = np.zeros(nt) p1[0] = 1.0 p1[-1] = 1.0 infin = m.Param(value=p1) # Equations m.Equation(x1.dt()==x2) m.Equation(x2.dt()==u) m.Equation(x3.dt()==x1) # Constraints m.Equation(infin*x1==0) m.Equation(final*x2==0) m.Obj(x3*final) # Objective function #m.fix(x2,pos=end_loc,val=0.0) m.options.IMODE = 6 # optimal control mode #m.solve(disp=True) # solve m.solve(disp=False) # solve plt.figure(1) # plot results plt.plot(m.time,x1.value,'k-',label=r'$x_1$') plt.plot(m.time,x2.value,'b-',label=r'$x_2$') plt.plot(m.time,x3.value,'g-',label=r'$x_3$') plt.plot(m.time,u.value,'r--',label=r'$u$') plt.legend(loc='best') plt.xlabel('Time') plt.ylabel('Value') plt.show() plt.figure(1) # plot results plt.plot(m.time,x1.value,'k-',label=r'$x_1$') plt.plot(m.time,(m.time-2)**2-2,'g--',label=r'$\hat x_1$') plt.legend(loc='best') plt.xlabel('Time') plt.ylabel('Value') plt.show() 

Answer: Use m.integral or m.vsum() to create a time weighted summation or vertical summation along the time direction. Here is a solution that replicates the exact solution.  from gekko import GEKKO import numpy as np import matplotlib.pyplot as plt m = GEKKO(remote=True) # initialize gekko nt = 501 m.time = np.linspace(0,2,nt) # insert a small time step at the beginning and end # for final and initial condition equation x(1e-8)+x(2-1e-8)=0 # this doesn't change the numerical solution significantly m.time = np.insert(m.time,1,1e-8) m.time = np.insert(m.time,len(m.time)-1,2.0-1e-8) nt += 2 # Variables x1 = m.Var(fixed_initial=False,lb=-100,ub=100) x2 = m.Var(fixed_initial=False) u = m.MV(fixed_initial=False,lb=-2,ub=2) u.STATUS = 1 p = np.zeros(nt) # mark final time point p[-2] = 1.0 final = m.Param(value=p) q = p.copy() q[1] = 1.0 final_initial = m.Param(value=q) xfi = m.Var() m.Equation(xfi==final_initial*x1) # Equations m.Equation(x1.dt()==x2) m.Equation(x2.dt()==u) m.Equation(final*m.vsum(xfi)==0) m.Equation(final*x2==0) m.Minimize(m.integral(x1)*final) # Objective function m.options.IMODE = 6 # optimal control mode m.options.NODES = 2 m.solve(disp=True) # solve plt.figure(1) # plot results plt.subplot(2,1,1) plt.plot(m.time,x1.value,'k-',label=r'$x_1$') plt.plot(m.time,x2.value,'b-',label=r'$x_2$') plt.plot(m.time,u.value,'--',color='orange',label=r'$u$') plt.legend(loc='best') plt.ylabel('Value') plt.subplot(2,1,2) plt.plot(m.time,x1.value,'k-',label=r'$x_1$') plt.plot(m.time,(m.time-2)**2-2,'r--',label=r'$\hat x_1$') plt.legend(loc='best') plt.xlabel('Time') plt.ylabel('Value') plt.show()  One issue is that x1(1-e8)+x1(2-1e-8)=0 is used as a constraint instead of x1(0)+x1(2)=0. The numerical solutions should be nearly equivalent and the 1e-8 can be further reduced.