python - 使用python GEKKO模拟矢量ODE

我正在尝试学习如何使用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？

最佳答案

模拟模型的一种方法：

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

关于python - 使用python GEKKO模拟矢量ODE，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/57979580/