时间戳

首页 > 解决方案 > 使用 python GEKKO 模拟向量 ODE

python - 使用 python GEKKO 模拟向量 ODE

问题描述

我正在尝试学习如何使用GEKKO python 包。作为第一步，我想模拟一个简单的向量 ODE：dx/dt = A*x，其中 A 是矩阵，x 是向量。我为 GEKKO 看到的所有 ODE 示例都是针对标量 ODE 的，在线文档中的数组示例没有显示如何在声明 Equation 时合并 .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？

标签： pythongekko

解决方案

模拟模型的一种方法：

dx/dt = A x

是声明x为 Array 并用于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()

另一种方法是使用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()

推荐阅读