使用 GEKKO 优化框架形成约束方程系统

Question

我正在尝试使用形式的双积分器动力学来解决一个简单的最小时间最优控制问题，

dx1/dt = x2
dx2/dt = u

配合GEKKO优化框架如下：

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

model = GEKKO(remote=False)

x1_initial = 0.0
x1_final = 10.0

x2_initial = 0.0
x2_final = 0.0

t_initial = 0.0
t_final = 25.0
num_timesteps = 1000
dt = (t_final - t_initial) / num_timesteps

x = model.Array(model.Var, (2, num_timesteps + 1))
u = model.Array(model.Var, num_timesteps + 1)
tf = model.Var()

for k in range(num_timesteps + 1):
    u[k].lower = -0.4
    u[k].upper = 0.4
    u[k].value = 0.0

for k in range(num_timesteps + 1):
    x[0, k].value = 5.0
    x[1, k].value = 0.0

tf.lower = t_initial
tf.upper = t_final
tf.value = t_final
dt = (tf - t_initial) / num_timesteps

def f(x, u, k):
    return np.array([x[1,k], u[k]])

for k in range(num_timesteps):
    model.Equations([x[:, k + 1] == x[:, k] + (dt/2.0)*(f(x, u, k + 1) + f(x, u, k))])
    # model.Equation(x[0, k + 1] == x[0, k] + (dt/2.0)*(x[1, k + 1] + x[1, k]))
    # model.Equation(x[1, k + 1] == x[1, k] + (dt/2.0)*(u[k + 1] + u[k]))

model.Equation(x[0, 0] == x1_initial)
model.Equation(x[0, num_timesteps] == x1_final)

model.Equation(x[1, 0] == x2_initial)
model.Equation(x[1, num_timesteps] == x2_final)

model.Minimize(tf)
model.options.solver = 3
model.solve()

# Plotting results
t = np.linspace(t_initial, tf.value, num_timesteps + 1)

u_optimal = []
for k in range(num_timesteps + 1):
    u_optimal.append(u[k].value)

x1_optimal = []
for k in range(num_timesteps + 1):
    x1_optimal.append(x[0, k].value)

x2_optimal = []
for k in range(num_timesteps + 1):
    x2_optimal.append(x[1, k].value)

plt.figure()
plt.plot(t, u_optimal)
plt.xlabel('time (s)')
plt.ylabel('u(t)')
plt.grid()

plt.figure()
plt.plot(t, x1_optimal)
plt.xlabel('time (s)')
plt.ylabel('x1(t)')
plt.grid()

plt.figure()
plt.plot(t, x2_optimal)
plt.xlabel('time (s)')
plt.ylabel('x2(t)')
plt.grid()

plt.show()

我要做的是使用梯形积分形成一个等式约束系统，然后使用 GEKKO 求解该系统以获得最佳控制输入。但是，使用函数定义，

def f(x, u, k):
    return np.array([x[1,k], u[k]])

结合等式约束系统，

for k in range(num_timesteps):
    model.Equations([x[:, k + 1] == x[:, k] + (dt/2.0)*(f(x, u, k + 1) + f(x, u, k))])

给我以下错误，

Exception: @error: Equation Definition
 Equation without an equality (=) or inequality (>,<)
 false
 STOPPING...

我在上面的代码 sn-p 中添加了两行注释代码，这将允许程序正确运行，但我希望避免将每个方程分开，因为我想扩展它对于处理更复杂系统动力学的问题，也可以使用更复杂的搭配方法而不是梯形方法。

我知道 GEKKO 有一些很好的动态优化功能，但我希望自己尝试实现各种直接搭配方法以更好地理解理论。

Answer 1

有一些orthogonal collocation on finite elements in the Machine Learning and Dynamic Optimization course 的例子。

from __future__ import division
import numpy as np
from scipy.optimize import fsolve
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# final time
tf = 1.0

# solve with ODEINT (for comparison)
def model(x,t):
    u = 4.0
    return (-x**2 + u)/5.0
t = np.linspace(0,tf,20)
y0 = 0
y = odeint(model,y0,t)
plt.figure(1)
plt.plot(t,y,'r-',label='ODEINT')

# ----------------------------------------------------
# Approach #1 - Write the model equations in Python
# ----------------------------------------------------
# define collocation matrices
def colloc(n):
    if (n==2):
        NC = np.array([[1.0]])
    if (n==3):
        NC = np.array([[0.75,-0.25], \
                       [1.00, 0.00]])
    if (n==4):
        NC = np.array([[0.436,-0.281, 0.121], \
                       [0.614, 0.064, 0.0461], \
                       [0.603, 0.230, 0.167]])
    if (n==5):
        NC = np.array([[0.278, -0.202, 0.169, -0.071], \
                       [0.398,  0.069, 0.064, -0.031], \
                       [0.387,  0.234, 0.278, -0.071], \
                       [0.389,  0.222, 0.389,  0.000]])
    if (n==6):
        NC = np.array([[0.191, -0.147, 0.139, -0.113, 0.047],
                       [0.276,  0.059, 0.051, -0.050, 0.022],
                       [0.267,  0.193, 0.252, -0.114, 0.045],
                       [0.269,  0.178, 0.384,  0.032, 0.019],
                       [0.269,  0.181, 0.374,  0.110, 0.067]])
    return NC

# define collocation points from Lobatto quadrature
def tc(n):
    if (n==2):
        time = np.array([0.0,1.0])
    if (n==3):
        time = np.array([0.0,0.5,1.0])
    if (n==4):
        time = np.array([0.0, \
                         0.5-np.sqrt(5)/10.0, \
                         0.5+np.sqrt(5)/10.0, \
                         1.0])
    if (n==5):
        time = np.array([0.0,0.5-np.sqrt(21)/14.0, \
                         0.5,0.5+np.sqrt(21)/14.0, 1])
    if (n==6):
        time = np.array([0.0, \
                         0.5-np.sqrt((7.0+2.0*np.sqrt(7.0))/21.0)/2.0, \
                         0.5-np.sqrt((7.0-2.0*np.sqrt(7.0))/21.0)/2.0, \
                         0.5+np.sqrt((7.0-2.0*np.sqrt(7.0))/21.0)/2.0, \
                         0.5+np.sqrt((7.0+2.0*np.sqrt(7.0))/21.0)/2.0, \
                         1.0])
    return time*tf

# solve with SciPy fsolve
def myFunction(z,*param):
    n = param[0]
    m = param[1]
    # rename z as x and xdot variables
    x = np.empty(n-1)
    xdot = np.empty(n-1)
    x[0:n-1] = z[0:n-1]
    xdot[0:n-1] = z[n-1:m]

    # initial condition (x0)
    x0 = 0.0
    # input parameter (u)
    u = 4.0
    # final time
    tn = tf

    # function evaluation residuals
    F = np.empty(m)
    # nonlinear differential equations at each node
    # 5 dx/dt = -x^2 + u
    F[0:n-1] = 5.0 * xdot[0:n-1] + x[0:n-1]**2 - u
    # collocation equations
    # tn * NC * xdot = x - x0
    NC = colloc(n)
    F[n-1:m] = tn * np.dot(NC,xdot) - x + x0 * np.ones(n-1)
    return F

sol_py = np.empty(5) # store 5 results
for i in range(2,7):
    n = i
    m = (i-1)*2
    zGuess = np.ones(m)
    z = fsolve(myFunction,zGuess,args=(n,m))
    # add to plot
    yc = np.insert(z[0:n-1],0,0)
    plt.plot(tc(n),yc,'o',markersize=10,label='Nodes = '+str(i))
    # store just the last x[n] value
    sol_py[i-2] = z[n-2]
plt.legend(loc='best')

# ----------------------------------------------------
# Approach #2 - Write model in APMonitor and let
#   modeling language create the collocation equations
# ----------------------------------------------------
# load GEKKO
from gekko import GEKKO

sol_apm = np.empty(5) # store 5 results
i = 0
for nodes in range(2,7):
    m = GEKKO(remote=False)

    u = m.Param(value=4)
    x = m.Var(value=0)
    m.Equation(5*x.dt() == -x**2 + u)

    m.time = [0,tf]

    m.options.imode = 4
    m.options.time_shift = 0
    m.options.nodes = nodes

    m.solve() # solve problem
    sol_apm[i] = x.value[-1] # store solution (last point)
    i += 1

# print the solutions
print(sol_py)
print(sol_apm)

# show plot
plt.ylabel('x(t)')
plt.xlabel('time')
plt.show()

您可以定义具有相同名称的变量（例如x）或使用m.Array(m.Var,n) 定义变量。要查看的一件事是通过在发送m.solve() 命令之前使用m.open_folder() 打开运行文件夹来查看模型文件。使用文本编辑器查看该文件夹中的.apm 文件。