Solving optimal control problem with constraint x(0) + x(2) =0 with GEKKO

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 1

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.