Построение переменных из системной функции дифференциальных уравнений

Question

У меня есть 4-4 системы дифференциальных уравнений в функции (subsystem4), и я решил ее с помощью функции odeint. Мне удалось построить результаты системы. Моя проблема в том, что я хочу построить и некоторые другие уравнения (например, x, y, vcxdot ...), которые включены в ту же функцию (subsystem4), но я получаю NameError: имя 'vcxdot' не определено. Кроме того, я хочу использовать некоторые из этих уравнений (а не только результаты системы уравнений) в качестве входных данных в следующей системе дифференциальных уравнений и построить все уравнения за один и тот же период времени (t). Я сделал это с помощью Matlab-Simulink, но это было намного проще из-за блоков Simulink. Как я могу получить доступ и построить все уравнения функции (subsystem4) и использовать их в качестве входных данных в следующей системе? Я новичок в Python и я использую Python 2.7.12. Заранее спасибо!

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def subsystem4(u,t):
    added_mass_x = 0.03 # kg
    added_mass_y = 0.04
    mb = 0.3 # kg
    m1 = mb-added_mass_x
    m2 = mb-added_mass_y
    l1 = 0.07 # m
    l2 = 0.05 # m
    J = 0.00050797 # kgm^2
    Sa = 0.0110 # m^2
    Cd = 2.44
    Cl = 3.41
    Kd = 0.000655 # kgm^2
    r = 1000 # kg/m^3
    f = 2 # Hz
    c1 = 0.5*r*Sa*Cd
    c2 = 0.5*r*Sa*Cl
    c3 = 0.5*mb*(l1**2)
    c4 = Kd/J
    c5 = (1/(2*J))*(l1**2)*mb*l2
    c6 = (1/(3*J))*(l1**3)*mb


    vcx = u[0]
    vcy = u[1]
    psi = u[2]
    wz = u[3]

    x = 3 + 0.3*np.cos(t)
    y = 0.5 + 0.3*np.sin(t)
    xdot = -0.3*np.sin(t)
    ydot = 0.3*np.cos(t)
    xdotdot = -0.3*np.cos(t)
    ydotdot = -0.3*np.sin(t)
    vcx = xdot*np.cos(psi)-ydot*np.sin(psi)
    vcy = ydot*np.cos(psi)+xdot*np.sin(psi)

    psidot = wz
    vcxdot = xdotdot*np.cos(psi)-xdot*np.sin(psi)*psidot-ydotdot*np.sin(psi)-ydot*np.cos(psi)*psidot
    vcydot = ydotdot*np.cos(psi)-ydot*np.sin(psi)*psidot+xdotdot*np.sin(psi)+xdot*np.cos(psi)*psidot
    g1 = -(m1/c3)*vcxdot+(m2/c3)*vcy*wz-(c1/c3)*vcx*np.sqrt((vcx**2)+(vcy**2))+(c2/c3)*vcy*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)
    g2 = (m2/c3)*vcydot+(m1/c3)*vcx*wz+(c1/c3)*vcy*np.sqrt((vcx**2)+(vcy**2))+(c2/c3)*vcx*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)

    A = 12*np.sin(2*np.pi*f*t+np.pi)
    if A>=0.1:
        wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz**2*np.sign(wz)-c5*g2-c6*np.sqrt((g1**2)+(g2**2))
    elif A<-0.1:
        wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz**2*np.sign(wz)-c5*g2+c6*np.sqrt((g1**2)+(g2**2))
    else:
        wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz**2*np.sign(wz)-c5*g2


    return [vcxdot,vcydot,psidot,wzdot]

u0 = [0,0,0,0]
t = np.linspace(0,15,1000)
u = odeint(subsystem4,u0,t)

vcx = u[:,0]
vcy = u[:,1]
psi = u[:,2]
wz = u[:,3]

plt.figure(1)
plt.subplot(211)
plt.plot(t,vcx,'r-',linewidth=2,label='vcx')
plt.plot(t,vcy,'b--',linewidth=2,label='vcy')
plt.plot(t,psi,'g:',linewidth=2,label='psi')
plt.plot(t,wz,'c',linewidth=2,label='wz')
plt.xlabel('time')
plt.legend()
plt.show()

Answer 1

К непосредственному вопросу о построении производных вы можете получить скорости, непосредственно вызвав функцию ODE для решения,

u = odeint(subsystem4,u0,t)
udot = subsystem4(u.T,t)

и получите отдельные массивы скоростей через

vcxdot,vcydot,psidot,wzdot = udot

---

В этом случае функция включает в себя ветвление, которое не очень удобно для ее векторизованных вызовов. Существуют способы векторизации ветвления, но самый простой обходной путь - это обход вручную точек решения, который медленнее, чем работающая векторизованная реализация. Это снова будет прокручивать список кортежей, таких как odeint, поэтому результат должен быть транспонирован как кортеж списков для «легкого» присвоения переменным одного массива.

udot = [ subsystem4(uk, tk) for uk, tk in zip(u,t) ]; 
vcxdot,vcydot,psidot,wzdot = np.asarray(udot).T

---

Может показаться, что это несколько удваивает вычисление, но не совсем, поскольку точки решения обычно интерполируются из точек внутреннего шага решателя. Оценка функции ODE во время интеграции обычно происходит в точках, отличных от точек решения.

Для других переменных извлеките вычисление положения и скоростей в функции, чтобы иметь постоянную и композицию только в одном месте:

def xy_pos(t): return 3 + 0.3*np.cos(t), 0.5 + 0.3*np.sin(t)
def xy_vel(t): return -0.3*np.sin(t), 0.3*np.cos(t) 
def xy_acc(t): return -0.3*np.cos(t), -0.3*np.sin(t)

или аналогичный, который вы затем можете использовать как внутри функции ODE, так и при подготовке графиков.

---

Скорее всего, Simulink соберет все уравнения всех блоков и сформирует их в одну большую систему ODE, которая затем решается для всего состояния сразу. Вам нужно будет реализовать нечто подобное. Один большой вектор состояния, и каждая подсистема знает свою часть состояния, соотв. производный вектор, чтобы получить его конкретные переменные состояния и записать производные в. Затем для вычисления производных можно использовать значения, передаваемые между подсистемами.

То, что вы пытаетесь сделать, решая подсистемы отдельно, работает только для соотв. скорее всего, приведет к методу интеграции 1-го порядка. Все методы более высокого порядка должны иметь возможность одновременно сдвигать состояние в некотором направлении, рассчитанном на предыдущих этапах метода, и оценивать всю систему там.