python - 使用 python GEKKO 模拟向量 ODE

Question

我正在尝试学习如何使用GEKKO python 包。第一步，我想模拟一个简单的向量 ODE:dx/dt = A*x，其中 A 是矩阵，x 是向量。我见过的 GEKKO 的所有 ODE 示例都是标量 ODE，数组示例来自在线 documentation没有显示在声明方程时如何合并 .dt() 方法。我知道在声明方程时可以使用列表，所以我认为这样的方法是可行的:

import numpy as np
from gekko import GEKKO
m=GEKKO()
m.time=np.linspace(0.,1.,10)
N=5
A=np.ones([N,N])
x=np.ones(N)
x=m.Var(value=x)
A=m.Param(value=A)
for i in range(N):
    for j in range(N):
        m.Equation(x[i].dt() += A[i][j] * x[j])

m.options.IMODE=4
m.solve()

但是这段代码会因为两个原因而失败:1) += 不是 Equation 方法的有效比较，2) python 提示 x[i].dt() 不是 x[i 的有效属性](这是一个 np.float64)。那么如果可能的话，我该如何在 GEKKO 中模拟这个 ODE？

Answer 1

模拟模型的一种方法:

dx/dt = A x

是将x声明为数组，并使用np.dot()与A的每一行进行矩阵乘法。

import numpy as np
from gekko import GEKKO
m=GEKKO()
m.time=np.linspace(0.,1.,10)
N=5
A=np.ones((N,N))
ic = array([1., 1., 1., 1., 1.])
x=m.Array(m.Var,N,value=0.) #initialize to zero
for i in range(N):
    x[i].value = ic[i] #set to some initial condition
m.Equations([x[i].dt()==np.dot(A[i,:],x) for i in range(N)])
m.options.IMODE=4
m.solve()

import matplotlib.pyplot as plt
for i in range(N):
    plt.plot(m.time,x[i].value)
plt.show()

另一种方法是使用 state space object in Gekko

dx/dt = A x + B u

y = C x + D u

其中B=0、C=0 和D=0。两种方法给出相同的结果。

import numpy as np
from gekko import GEKKO
m=GEKKO()
m.time=np.linspace(0.,1.,10)
N=5
A=np.ones((N,N))
B=np.zeros((N,1))
C=np.zeros((1,N))
x,y,u = m.state_space(A,B,C,D=None)
for i in range(N):
    x[i].value=1
m.options.IMODE=4
m.solve()

import matplotlib.pyplot as plt
for i in range(N):
    plt.plot(m.time,x[i].value)
plt.show()