使用 GEKKO 求解约束 x(0) + x(2) =0 的最优控制问题

Question

我开始学习 Gekko，我正在测试最优控制问题。我正在尝试用 Gekko 解决以下最优控制问题

这个问题的解是(x_1(t) = (t-2)^2 - 2) 如何建立约束 x(0) + x(2) = 0？ 我的代码给了我一个错误的解决方案。

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

使用m.integral 或m.vsum() 沿时间方向创建时间加权求和或垂直求和。这是一个复制精确解决方案的解决方案。

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()

一个问题是x1(1-e8)+x1(2-1e-8)=0 被用作约束而不是x1(0)+x1(2)=0。数值解应该几乎相等，1e-8 可以进一步减少。