Задержка дифференциальных уравнений в питоне - PullRequest
Новая
       80

Задержка дифференциальных уравнений в питоне

0 голосов
/ 05 апреля 2019

У меня есть эта проблема управления, где моя цель состоит в том, чтобы выбрать значения D, чтобы Hb(t) оставался в определенном диапазоне, см. Уравнения здесь (если это поможет, я пытаюсь воспроизведите симуляцию, обсуждаемую в этой статье ).

Я не уверен в своей попытке решить эту проблему, в частности, в терминах, где необходимы значения прошлого состояния E и P. Я пытался вычислить их вручную и передать в качестве аргументов моей интегрирующей функции. Но этот подход:

  • имеет повторный код и его трудно понять
  • неверно, я думаю, потому что, как и после решения, я получаю Hb(t) значений, которые я не ожидал бы даже после изменения многих значений параметров.

Есть ли лучший способ решить эти дифференциальные уравнения с состояниями, связанными со временем? 

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

def erythropoiesis_func(s, t, d, term1, term2, term3, term4):
    """
    Inputs
    d     = drug dose (ug/kg per week)
    term1 = avg of past Tp days
    term2 = E_{tot}(t - (Tp + Tm))
    term3 = P(t - (Tp + Tm))
    term4 = sum of past Tj days times normal
    """

    """
    States of model s(E,P,R):

    E: exogenous plasma concentration of EPO
    P: bone marrow concentration of progenitor cells and precursor cells
    R: concentration of red blood cells (RDCs)

    """
    E = s[0]                
    P = s[1]                
    R = s[2]

    # Parameters to vary by each trajectory (patient) but constant for now
    Ep = .36  #np.random.normal(.3588, .0753, 1)
    Cr = .14 #np.random.normal(.1372, .0520, 1)
    Cp = .21 #np.random.normal(.2014, .0640, 1)

    # Constants:
    Vd  = 52.4            # Volume distribution
    E50 = 100 / Vd        # half-maximal effective concentration (drug potency)
    Tp  = 9               # P cell life in days
    Tm  = 4               # Delay from P compartment to R compartment in days
    Tr  = 70              # R cell life in days


    Etot = (d / Vd) + Ep

    Spmr = np.random.normal(Tr , 30, size=Tr).sum()


    # Compute sdot:
    dsdt = np.empty(len(s))
    dsdt[0] = (24/25)*np.log(2)*Etot                      # E'(t)
    dsdt[1] = Cp*(Etot/(E50 + Etot))*P - Cp*term1         # P'(t)
    dsdt[2] = Cr*(term2/(E50 + term2))*term3 - (Cr/Spmr)*term4

    return dsdt

# Initial Condition for the Control
d0 = .5

# Initial Conditions for the States (E, P, R)
s0 = np.array([(d0/52.4), 1., 1.])

# Constants
Tp  = 9
Tm  = 4               
Tr  = 70
E50 = 100 / 52.4
MCH = 2.4

# Time Interval (days)
tt = np.array(range(Tp+Tm+Tr))

# Pre fill vectors
E = np.ones(len(tt)) * s0[0]
P = np.ones(len(tt)) * s0[1]
R = np.ones(len(tt)) * s0[2]
H = MCH * R

d = np.random.uniform(low=0.0, high=1.0, size=len(tt))

# Simulate
term1_vals = []
term4_vals = []
for t in range(len(tt)-1):
    # time to feed to solver
    ts = [tt[t],tt[t+1]]

    # generating term1
    temp1 = []
    for tj in range(Tp):
        if tj <= t:
            temp1.append( (E[t - tj] / (E50 + E[t - tj] )) * P[t - tj] )
    term1_vals.append(np.mean(temp1))

    # collecting term2 and term3 if they exist
    if (Tp + Tm) <= t:
        term2 = E[t-(Tp + Tm)]
        term3 = P[t-(Tp + Tm)]
    else:
        term2 = 0
        term3 = 0

    # generating term4
    temp4 = []
    for tj in range((Tp+Tm),(Tp+Tm+Tr)):
        if tj <= t:
            temp4.append( np.random.normal(Tr , 30, size=1) * (E[t - tj] / (E50 + E[t - tj] )) * P[t - tj] )
    term4_vals.append(np.sum(temp4))

    # calling solver
    s = odeint(erythropoiesis_func,s0,ts,args=(d[t+1],term1_vals[t],term2,term3,term4_vals[t]))

    #storing values
    E[t+1] = s[-1][0]
    P[t+1] = s[-1][1]
    R[t+1] = s[-1][2]
    H[t+1] = MCH * s[-1][2]
    s0 = s[-1]
...