Моделирование векторного ODE с использованием Python GEKKO

Я пытаюсь научиться использовать пакет python GEKKO. В качестве первого шага я хотел бы смоделировать простое векторное ОДУ: dx / dt = A * x, где A - матрица, а x - вектор. Все примеры ODE, которые я видел для GEKKO, были для скалярных ODE и пример массива из онлайн документация не показывает, как включить метод .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). Итак, как мне смоделировать это ODE в GEKKO, если это вообще возможно?


python gekko
person Mitch    schedule 17.09.2019    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Один из способов смоделировать модель:

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

Другой способ - использовать объект пространства состояний в 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()
person John Hedengren    schedule 18.09.2019