# 用 python 求拥塞控制模型欧拉数值解

d x 1 d t = y 1 + I z z + R − x 1 \dfrac{dx_1}{dt}=\dfrac{\frac{y_1+I}{z}}{z+R}-x_1 【E = bw / delay 表达式】
d w 1 d t = C ⋅ y 1 + I y 1 + y 2 + y 3 + 3 ⋅ I − w 1 \dfrac{dw_1}{dt}=C\cdot\dfrac{y_1+I}{y_1+y_2+y_3+3\cdot I}-w_1 【w = C * buffer 占率】
d y 1 d t = w 1 ⋅ z − y 1 \dfrac{dy_1}{dt}=w_1\cdot z-y_1 【bdp 公式】

d x 2 d t = y 2 + I z z + R − x 2 \dfrac{dx_2}{dt}=\dfrac{\frac{y_2+I}{z}}{z+R}-x_2
d w 2 d t = C ⋅ y 2 + I y 1 + y 2 + y 3 + 3 ⋅ I − w 2 \dfrac{dw_2}{dt}=C\cdot\dfrac{y_2+I}{y_1+y_2+y_3+3\cdot I}-w_2
d y 2 d t = w 2 ⋅ z − y 2 \dfrac{dy_2}{dt}=w_2\cdot z-y_2

d x 3 d t = y 3 + I z z + R − x 3 \dfrac{dx_3}{dt}=\dfrac{\frac{y_3+I}{z}}{z+R}-x_3
d w 3 d t = C ⋅ y 3 + I y 1 + y 2 + y 3 + 3 ⋅ I − w 3 \dfrac{dw_3}{dt}=C\cdot\dfrac{y_3+I}{y_1+y_2+y_3+3\cdot I}-w_3
d y 3 d t = w 3 ⋅ z − y 3 \dfrac{dy_3}{dt}=w_3\cdot z-y_3

d z d t = β ⋅ ( y 1 + y 2 + y 3 ) \dfrac{dz}{dt}=\beta\cdot(y_1+y_2+y_3) 【排队时延与 buffer 占率成正比】

d y d t = f ( t , y ) , y ( t 0 ) = y 0 \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0

• 离散：选一个时间步长 h，定义离散时间点 t n = t 0 + n h ( n = 0 , 1 , 2 , … ) t_n = t_0 + nh \quad (n = 0, 1, 2, \ldots)
• 更新：使用欧拉法的更新公式： y n + 1 = y n + h ⋅ f ( t n , y n ) y_{n+1} = y_n + h \cdot f(t_n, y_n)
• 迭代：从初始条件开始，迭代计算 y n + 1 y_{n+1} 直到达到所需的时间范围。

#!/opt/homebrew/bin/python3

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator

class Model:
def __init__(self, initial_conditions, T, dt, a, b, c, d, R, C, I):
self.T = T
self.dt = dt
self.a = a
self.b = b
self.c = c
self.d = d
self.R = R
self.C = C
self.I = I
self.len = (len(initial_conditions) - 1) / 3
self.times = np.arange(0, T + dt, dt)

# 初始化状态
self.states = {name: np.zeros_like(self.times) for name in initial_conditions.keys()}
for name, value in initial_conditions.items():
self.states[name][0] = value

def update(self):
for n in range(1, len(self.times)):
for i in range(1, int(self.len) + 1):  # 假设有三个 x, y, w
x_name = f'x{i}'
y_name = f'y{i}'
w_name = f'w{i}'

self.states[x_name][n] = self.states[x_name][n - 1] + self.dt * (
(self.states[y_name][n - 1] + self.I) / self.states['z'][n - 1] /
(self.states['z'][n - 1] + self.R) - self.states[x_name][n - 1]
)

self.states[w_name][n] = self.states[w_name][n - 1] + self.dt * (
self.a * self.C * (self.states[y_name][n - 1] + self.I) /
(self.states['y1'][n - 1] + self.states['y2'][n - 1] + self.states['y3'][n - 1] + 3 * self.I) - self.states[w_name][n - 1]
)

self.states[y_name][n] = self.states[y_name][n - 1] + self.dt * (
self.states[w_name][n] * self.states['z'][n] - self.states[y_name][n - 1]
)
print("y1", self.states[y_name][n]);

self.states['z'][n] = self.states['z'][n - 1] + self.dt * (
self.states['y1'][n - 1] + self.states['y2'][n - 1] + self.states['y3'][n - 1]
)

def plot(self):
plt.figure(figsize = (12, 6))
plt.subplot(2, 1, 1)
plt.plot(self.times, self.states['z'], label = 'z(t)', linestyle = 'dotted')
for i in range(1, int(self.len) + 1):
plt.plot(self.times, self.states[f'y{i}'], label = f'y{i}(t)', linestyle = 'dashed')
plt.plot(self.times, self.states[f'w{i}'], label = f'w{i}(t)', linestyle = 'dashed')
plt.legend()

plt.subplot(2, 1, 2)
ax = plt.gca()
y_major = MultipleLocator(0.2)
ax.yaxis.set_major_locator(y_major)
plt.ylim(top = 0.01)
for i in range(1, int(self.len) + 1):
plt.plot(self.times, self.states[f'x{i}'], label = f'x{i}(t)')

plt.xlabel('time (t)')
plt.ylabel('')
plt.legend()
plt.grid(linestyle = 'dotted', linewidth = 1)
plt.tight_layout()
plt.show()

initial_conditions = {
'x1': 0.00002, 'y1': 6, 'x2': 0.00003, 'y2': 2,
'x3': 0.000001, 'y3': 14, 'w1': 1, 'w2': 1,
'w3': 3, 'z': 0.1
}

T, dt = 50, 0.1
a, b, c, d, R, C, I = 1, 1.000, 0.2, 0.03, 2, 18, 1

model = Model(initial_conditions, T, dt, a, b, c, d, R, C, I)
model.update()
model.plot()


self.states['y1'][n - 1] + self.states['y2'][n - 1] + self.states['y3'][n - 1] + 3 * self.I


sum(self.states['y'][j][n-1] for j in range(self.len)) + len * self.I


x1 = np.zeros_like(times)
y1 = np.zeros_like(times)
x2 = np.zeros_like(times)
y2 = np.zeros_like(times)
x3 = np.zeros_like(times)
y3 = np.zeros_like(times)
w1 = np.zeros_like(times)
w2 = np.zeros_like(times)
w3 = np.zeros_like(times)
z = np.zeros_like(times)

x1[0], y1[0], x2[0], y2[0], x3[0], y3[0], w1[0], w2[0], w3[0], z[0] = x10, y10, x20, y20, x30, y30, w10, w20, w30, z0
for n in range(1, len(times)):
x1[n] = x1[n-1] + dt * (...)
w1[n] = w1[n-1] + dt * (...)
y1[n] = y1[n-1] + dt * (...)
x2[n] = x2[n-1] + dt * (...)
w2[n] = w2[n-1] + dt * (...)
y2[n] = y2[n-1] + dt * (...)
x3[n] = x3[n-1] + dt * (...)
w3[n] = w3[n-1] + dt * (...)
y3[n] = y3[n-1] + dt * (...)
z[n] = z[n-1] + dt * (...)
...


• 14
点赞
• 10
收藏
觉得还不错? 一键收藏
• 1
评论
08-12 1312
11-15
09-05 964
09-10 708
09-06 1431
09-10 194

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。