使用 python GEKKO 模拟矢量 ODE

Question

我正在尝试学习如何使用 GEKKO python package. As a first step, I would like to simulate a simple vector ODE: dx/dt = A*x, where A is a matrix and x is a vector. All the ODE examples I have seen for GEKKO have been for scalar ODEs and the array example from the online 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()

使用 python GEKKO 模拟矢量 ODE

Simulating a vector ODE using python GEKKO

python

gekko