Python中的系统动力学模型：利用SciPy实现微分方程组的数值求解

Python中的系统动力学模型：利用SciPy实现微分方程组的数值求解

大家好，今天我们将深入探讨如何使用Python和SciPy库构建和求解系统动力学模型。系统动力学是一种用于理解复杂系统行为的建模方法，它通过使用微分方程组来描述系统中各个变量之间的相互作用。 SciPy库提供了强大的数值积分工具，使得我们能够对这些微分方程组进行求解，从而模拟系统的演化过程。

系统动力学建模基础

在开始编写代码之前，我们需要了解系统动力学的基本概念。系统动力学模型的核心是存量 (Stocks) 和 流量 (Flows)。

• 存量 (Stocks)：代表系统中累积的量，例如人口数量、资金量、库存量等。存量的变化率由流量决定。
• 流量 (Flows)：代表进入或离开存量的速率，例如出生率、死亡率、投资额、消耗率等。流量通常是存量和其他辅助变量的函数。
• 辅助变量 (Auxiliary Variables)：用于简化模型，将复杂的计算过程分解成更小的步骤，提高模型的可读性和可维护性。
• 参数 (Parameters)：模型的常数，例如初始人口、利率等。

一个简单的系统动力学模型可以用以下公式表示：

d(Stock)/dt = Inflow - Outflow

这个公式说明，存量随时间的变化率等于流入量减去流出量。

使用Python和SciPy构建系统动力学模型

现在，我们以一个经典的SIR模型 (Susceptible-Infected-Recovered) 为例，演示如何使用Python和SciPy构建系统动力学模型。 SIR模型用于描述传染病在人群中的传播过程。

SIR模型描述:

• S (Susceptible)：易感人群，即尚未感染疾病的人群。
• I (Infected)：感染人群，即正在感染疾病的人群。
• R (Recovered)：康复人群，即已经康复并获得免疫力的人群。

模型假设:

• 人群总数保持不变。
• 感染率与易感人群和感染人群的数量成正比。
• 康复率与感染人群的数量成正比。

微分方程组:

dS/dt = -beta * S * I
dI/dt = beta * S * I - gamma * I
dR/dt = gamma * I

其中：

• beta 是感染率常数。
• gamma 是康复率常数。

Python代码实现:

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

# 定义SIR模型的微分方程组
def sir_model(t, y, beta, gamma):
  """
  SIR模型的微分方程组。

  Args:
    t: 时间 (无意义参数，solve_ivp要求必须有).
    y: 包含S, I, R的数组.
    beta: 感染率常数.
    gamma: 康复率常数.

  Returns:
    包含dS/dt, dI/dt, dR/dt的数组.
  """
  S, I, R = y
  dSdt = -beta * S * I
  dIdt = beta * S * I - gamma * I
  dRdt = gamma * I
  return [dSdt, dIdt, dRdt]

# 设置模型参数
N = 1000  # 总人口
I0 = 1    # 初始感染人数
R0 = 0    # 初始康复人数
S0 = N - I0 - R0  # 初始易感人数
beta = 0.3  # 感染率
gamma = 0.1  # 康复率
t_start = 0  # 模拟起始时间
t_end = 100  # 模拟结束时间
t_eval = np.linspace(t_start, t_end, 100)  # 模拟时间点

# 初始条件
y0 = [S0, I0, R0]

# 使用solve_ivp求解微分方程组
sol = solve_ivp(sir_model, [t_start, t_end], y0, args=(beta, gamma), dense_output=True, t_eval=t_eval)

# 提取结果
S = sol.y[0]
I = sol.y[1]
R = sol.y[2]
t = sol.t

# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(t, S, label='Susceptible')
plt.plot(t, I, label='Infected')
plt.plot(t, R, label='Recovered')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('SIR Model Simulation')
plt.legend()
plt.grid(True)
plt.show()

代码解释:

1. 导入库: 导入numpy用于数值计算，scipy.integrate中的solve_ivp用于求解微分方程组，matplotlib.pyplot用于绘制结果。
2. 定义SIR模型: 定义sir_model函数，该函数接受时间t，状态变量y (包含S, I, R)，以及参数beta和gamma作为输入。函数返回一个包含dS/dt, dI/dt, dR/dt的数组，表示各个状态变量的变化率。
3. 设置模型参数: 设置模型的初始条件和参数，包括总人口N，初始感染人数I0，初始康复人数R0，感染率beta，康复率gamma，模拟起始时间t_start，模拟结束时间t_end，以及模拟时间点t_eval。
4. 初始条件: 将S0, I0, R0封装成列表y0，作为solve_ivp的初始条件。
5. 使用solve_ivp求解微分方程组: 使用solve_ivp函数求解微分方程组。solve_ivp函数接受以下参数：
   • fun: 微分方程组的函数 (sir_model).
   • t_span: 时间区间 [t_start, t_end].
   • y0: 初始条件.
   • args: 传递给微分方程组函数的额外参数 (beta, gamma).
   • dense_output: True表示返回一个连续解，允许在任意时间点计算解.
   • t_eval: 指定需要计算解的时间点.
6. 提取结果: 从sol对象中提取S, I, R的值和时间t。
7. 绘制结果: 使用matplotlib.pyplot绘制S, I, R随时间变化的曲线。

运行这段代码，你将会看到一个SIR模型的模拟结果图，显示了易感人群、感染人群和康复人群随时间的变化趋势。

模型的改进与扩展

SIR模型只是一个简单的示例，我们可以根据实际情况对模型进行改进和扩展。

1. 添加出生率和死亡率:

现实世界中，人口数量会受到出生率和死亡率的影响。我们可以将出生率和死亡率添加到SIR模型中：

dS/dt = mu * N - beta * S * I - mu * S
dI/dt = beta * S * I - gamma * I - mu * I
dR/dt = gamma * I - mu * R

其中：

• mu 是出生率和死亡率 (假设相等)。
• N 是总人口。

2. 添加潜伏期 (SEIR模型):

在某些传染病中，感染者在感染后会经历一个潜伏期，在此期间他们不会传播疾病。我们可以引入一个新的状态变量E (Exposed)来表示潜伏期的人群，从而构建SEIR模型：

dS/dt = -beta * S * I
dE/dt = beta * S * I - sigma * E
dI/dt = sigma * E - gamma * I
dR/dt = gamma * I

其中：

• sigma 是潜伏期的转化率。

3. 添加疫苗接种:

疫苗接种可以降低易感人群的数量，从而减缓疾病的传播。我们可以将疫苗接种的影响添加到SIR模型中：

dS/dt = -beta * S * I - v * S
dI/dt = beta * S * I - gamma * I
dR/dt = gamma * I + v * S

其中：

• v 是疫苗接种率。

你可以根据实际需求，将这些改进和扩展添加到你的模型中。

复杂系统建模：供应链案例

现在，让我们考虑一个更复杂的系统动力学模型：一个简单的供应链模型。这个模型包含三个主要的存量：

• 库存 (Inventory)：零售商的库存量。
• 在途订单 (Orders in Transit)：零售商已下单但尚未收到的订单量。
• 积压订单 (Backlog)：零售商未满足的客户订单量。

流量包括:

• 销售 (Sales)：零售商卖出的产品数量。
• 订单 (Orders)：零售商向供应商发出的订单数量。
• 交货 (Deliveries)：供应商交付给零售商的产品数量。

模型假设:

• 零售商根据库存目标和在途订单来确定订单数量。
• 供应商需要一定的交货时间才能完成订单。
• 如果库存不足，客户订单会积压。

模型方程:

为了简化模型，我们引入一些辅助变量：

• 需求 (Demand)：客户对产品的需求量 (假设为常数)。
• 库存目标 (Inventory Target)：零售商希望维持的库存水平。
• 订单填充时间 (Order Fulfillment Time)：供应商完成订单所需的时间。
• 库存调整时间 (Inventory Adjustment Time)：零售商调整库存所需的时间。
• 积压订单填充时间 (Backlog Fulfillment Time)：零售商填充积压订单所需的时间。

使用这些辅助变量，我们可以定义以下方程：

Inventory Target = Demand * Inventory Adjustment Time
Orders = Demand + (Inventory Target - Inventory) / Inventory Adjustment Time + Backlog / Backlog Fulfillment Time
Deliveries = Orders in Transit / Order Fulfillment Time
Sales = min(Demand, Inventory + Backlog)
d(Inventory)/dt = Deliveries - Sales
d(Orders in Transit)/dt = Orders - Deliveries
d(Backlog)/dt = Demand - Sales

Python代码实现:

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

def supply_chain_model(t, y, demand, inventory_adjustment_time, backlog_fulfillment_time, order_fulfillment_time):
  """
  简单的供应链模型。

  Args:
    t: 时间.
    y: 包含Inventory, Orders in Transit, Backlog的数组.
    demand: 客户需求量.
    inventory_adjustment_time: 库存调整时间.
    backlog_fulfillment_time: 积压订单填充时间.
    order_fulfillment_time: 订单填充时间.

  Returns:
    包含d(Inventory)/dt, d(Orders in Transit)/dt, d(Backlog)/dt的数组.
  """
  inventory, orders_in_transit, backlog = y

  # 辅助变量
  inventory_target = demand * inventory_adjustment_time
  orders = demand + (inventory_target - inventory) / inventory_adjustment_time + backlog / backlog_fulfillment_time
  deliveries = orders_in_transit / order_fulfillment_time
  sales = min(demand, inventory + backlog)

  # 流量
  d_inventory_dt = deliveries - sales
  d_orders_in_transit_dt = orders - deliveries
  d_backlog_dt = demand - sales

  return [d_inventory_dt, d_orders_in_transit_dt, d_backlog_dt]

# 设置模型参数
demand = 100  # 客户需求量
inventory_adjustment_time = 5  # 库存调整时间
backlog_fulfillment_time = 2  # 积压订单填充时间
order_fulfillment_time = 10  # 订单填充时间
t_start = 0
t_end = 100
t_eval = np.linspace(t_start, t_end, 100)

# 初始条件
initial_inventory = 500
initial_orders_in_transit = 0
initial_backlog = 0
y0 = [initial_inventory, initial_orders_in_transit, initial_backlog]

# 使用solve_ivp求解微分方程组
sol = solve_ivp(supply_chain_model, [t_start, t_end], y0, args=(demand, inventory_adjustment_time, backlog_fulfillment_time, order_fulfillment_time), dense_output=True, t_eval=t_eval)

# 提取结果
inventory = sol.y[0]
orders_in_transit = sol.y[1]
backlog = sol.y[2]
t = sol.t

# 绘制结果
plt.figure(figsize=(12, 8))
plt.plot(t, inventory, label='Inventory')
plt.plot(t, orders_in_transit, label='Orders in Transit')
plt.plot(t, backlog, label='Backlog')
plt.xlabel('Time')
plt.ylabel('Quantity')
plt.title('Supply Chain Model Simulation')
plt.legend()
plt.grid(True)
plt.show()

代码解释:

1. 定义supply_chain_model函数: 该函数接受时间t，状态变量y (包含Inventory, Orders in Transit, Backlog)，以及模型参数作为输入。函数计算并返回各个状态变量的变化率。
2. 设置模型参数: 设置模型的初始条件和参数，包括需求demand，库存调整时间inventory_adjustment_time，积压订单填充时间backlog_fulfillment_time，订单填充时间order_fulfillment_time，模拟起始时间t_start，模拟结束时间t_end，以及模拟时间点t_eval。
3. 使用solve_ivp求解微分方程组: 使用solve_ivp函数求解微分方程组，将supply_chain_model函数作为参数传递给solve_ivp。
4. 提取结果: 从sol对象中提取Inventory, Orders in Transit, Backlog的值和时间t。
5. 绘制结果: 使用matplotlib.pyplot绘制Inventory, Orders in Transit, Backlog随时间变化的曲线。

运行这段代码，你将会看到一个简单的供应链模型的模拟结果图，显示了库存、在途订单和积压订单随时间的变化趋势。 通过调整模型的参数，你可以观察不同参数对系统行为的影响，例如订单填充时间对库存波动的影响。

SciPy solve_ivp 的高级应用

solve_ivp 函数提供了许多高级功能，可以帮助我们更精确、更有效地求解微分方程组。

1. 设置容差 (Tolerance):

solve_ivp 允许我们设置绝对容差 (atol) 和相对容差 (rtol)，以控制解的精度。 默认情况下，atol 为 1e-6，rtol 为 1e-3。 减小容差可以提高解的精度，但也会增加计算时间。

sol = solve_ivp(sir_model, [t_start, t_end], y0, args=(beta, gamma), atol=1e-8, rtol=1e-6, t_eval=t_eval)

2. 选择求解器 (Solver):

solve_ivp 提供了多种求解器，例如 RK45 (默认), RK23, DOP853, Radau, BDF, LSODA。 不同的求解器适用于不同类型的微分方程组。 例如，对于刚性 (Stiff) 微分方程组，Radau 或 BDF 求解器可能更有效。

sol = solve_ivp(sir_model, [t_start, t_end], y0, args=(beta, gamma), method='Radau', t_eval=t_eval)

3. 事件处理 (Event Handling):

solve_ivp 允许我们定义事件函数，当满足特定条件时，事件函数会被触发。 这对于模拟系统中的突变或开关行为非常有用。

例如，我们可以定义一个事件函数，当感染人数超过总人口的10%时触发：

def infection_threshold(t, y, beta, gamma):
  return y[1] - 0.1 * N

infection_threshold.terminal = True  # 停止求解
infection_threshold.direction = 1  # 仅在函数值从负变为正时触发

sol = solve_ivp(sir_model, [t_start, t_end], y0, args=(beta, gamma), events=infection_threshold, t_eval=t_eval)

在这个例子中，infection_threshold 函数返回感染人数与总人口的10%之差。 terminal = True 表示当事件触发时，求解过程将会停止。 direction = 1 表示仅在函数值从负变为正时触发事件。

总结：系统动力学建模与数值求解

我们学习了如何使用Python和SciPy库构建和求解系统动力学模型。 从简单的SIR模型到复杂的供应链模型，我们展示了如何使用微分方程组来描述系统行为，并使用solve_ivp函数进行数值求解。 了解了solve_ivp函数的高级功能，例如设置容差、选择求解器和事件处理，这些功能可以帮助我们更精确、更有效地模拟复杂系统的行为。

更多IT精英技术系列讲座，到智猿学院