自由终端时间、积分目标和微分方程作为约束

我正在尝试解决一个最优控制问题，该问题涉及最小化具有固定状态但自由终端时间的积分目标。这是一个相对简单的问题，可以通过解析来解决。 Gekko 的解决方案与解析不符。

我不确定我做错了什么。我遵循了几个 Gekko 示例来解决这个问题。任何帮助深表感谢。

我遇到的另一个问题是如何让 Gekko 自动计算控制的初始值。最优控制总是从指定的初始控制猜测开始。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# create GEKKO model
m = GEKKO()
# time points
n = 501
tm = np.linspace(0, 1, n)
m.time = tm

# Variables
x1 = m.Var(value=1)  # x1
x2 = m.Var(value=2)  # x2
# u = m.Var(value=-1)  # control variable used as normal var
u = m.MV(value=-1) # manipulative variable
u.STATUS = 1
u.DCOST = 1e-5

p = np.zeros(n)
p[-1] = 1.0
final = m.Param(value=p)

# FV
tf = m.FV(value=10.0, lb=0.0, ub=100.0)
tf.STATUS = 1

# equations
m.Equation(x1.dt()/tf == x2)
m.Equation(x2.dt()/tf == u)


# Final conditions
soft = True
if soft:
    # soft terminal constraint
    m.Minimize(final*1e5*(x1-3)**2)
    # m.Minimize(final*1e5*(x2-2)**2)
else:
    # hard terminal constraint
    x1f = m.Param()
    m.free(x1f)
    m.fix_final(x1f, 3)
    # connect endpoint parameters to x1 and x2
    m.Equations([x1f == x1])

# Objective Function
obj = m.Intermediate(tf*final*m.integral(0.5*u**2))

m.Minimize(final*obj)

m.options.IMODE = 6
m.options.NODES = 2
m.options.SOLVER = 3
m.options.MAX_ITER = 500
# m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
m.solve(disp=False)


# Create a figure
plt.figure(figsize=(10, 4))
plt.subplot(2, 2, 1)
# plt.plot([0,1],[1/9,1/9],'r:',label=r'$x<\frac{1}{9}$')
plt.plot(tm, x1.value, 'k-', lw=2, label=r'$x1$')
plt.ylabel('x1')
plt.legend(loc='best')
plt.subplot(2, 2, 2)
plt.plot(tm, x2.value, 'b--', lw=2, label=r'$x2$')
plt.ylabel('x2')
plt.legend(loc='best')
plt.subplot(2, 2, 3)
plt.plot(tm, u.value, 'r--', lw=2, label=r'$u$')
plt.ylabel('control')
plt.legend(loc='best')
plt.xlabel('Time')
plt.subplot(2, 2, 4)
plt.plot(tm, obj.value, 'g-', lw=2, label=r'$\frac{1}{2} \int u^2$')
plt.text(0.5, 3.0, 'Final Value = '+str(np.round(obj.value[-1], 2)))
plt.ylabel('Objective')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()

以下是一些修改：

# u = m.MV(value=-1)
u = m.MV(value=-1,fixed_initial=False)

#obj = m.Intermediate(tf*final*m.integral(0.5*u**2))
obj = m.Intermediate(m.integral(0.5*u**2))

m.options.NODES = 3 # increase accuracy

如果添加一个约束tf<=3然后它给出与上面相同的解决方案。

然而，如果你放松tf约束到<=100那么就有更好的解决方案了。

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# create GEKKO model
m = GEKKO()
# time points
n = 501
tm = np.linspace(0, 1, n)
m.time = tm

# Variables
x1 = m.Var(value=1)  # x1
x2 = m.Var(value=2)  # x2
u = m.MV(value=-1,fixed_initial=False) # manipulated variable
u.STATUS = 1
u.DCOST = 1e-5

p = np.zeros(n)
p[-1] = 1.0
final = m.Param(value=p)

# FV
tf = m.FV(value=10.0, lb=0.0, ub=100.0)
tf.STATUS = 1

# equations
m.Equation(x1.dt()/tf == x2)
m.Equation(x2.dt()/tf == u)


# Final conditions
soft = True
if soft:
    # soft terminal constraint
    m.Minimize(final*1e5*(x1-3)**2)
    # m.Minimize(final*1e5*(x2-2)**2)
else:
    # hard terminal constraint
    x1f = m.Param()
    m.free(x1f)
    m.fix_final(x1f, 3)
    # connect endpoint parameters to x1 and x2
    m.Equations([x1f == x1])

# Objective Function
obj = m.Intermediate(m.integral(0.5*u**2))

m.Minimize(final*obj)

m.options.IMODE = 6
m.options.NODES = 3
m.options.SOLVER = 3
m.options.MAX_ITER = 500
# m.options.MV_TYPE = 0
m.options.DIAGLEVEL = 0
m.solve(disp=True)

# Create a figure
tm = tm*tf.value[0]
plt.figure(figsize=(10, 4))
plt.subplot(2, 2, 1)
# plt.plot([0,1],[1/9,1/9],'r:',label=r'$x<\frac{1}{9}$')
plt.plot(tm, x1.value, 'k-', lw=2, label=r'$x1$')
plt.ylabel('x1')
plt.legend(loc='best')
plt.subplot(2, 2, 2)
plt.plot(tm, x2.value, 'b--', lw=2, label=r'$x2$')
plt.ylabel('x2')
plt.legend(loc='best')
plt.subplot(2, 2, 3)
plt.plot(tm, u.value, 'r--', lw=2, label=r'$u$')
plt.ylabel('control')
plt.legend(loc='best')
plt.xlabel('Time')
plt.subplot(2, 2, 4)
plt.plot(tm, obj.value, 'g-', lw=2, label=r'$\frac{1}{2} \int u^2$')
plt.text(0.5, 3.0, 'Final Value = '+str(np.round(obj.value[-1], 2)))
plt.ylabel('Objective')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()