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

昨天使用 scipy 的 odeint 模拟了 E_best 的微分方程组模型(参见 用 python scipy 库模拟拥塞控制模型)，但我觉得那个模型中处理 z 时不够优雅，只是一个负反馈，并未体现 “排队时延与 buffer 占用率成比例增长” 的事实，所以今天我直接用 inflight-in-buffer 重新建模，就优雅多了。

设 x 为 E = bw / delay，y 为 inflight-in-buffer 但不包括余量，w 为到达瓶颈的速率(等价于 sender 的 pacing rate 但又不是，具体为何，看我前面的文章)，z 为排队时延，并且有以下参数，R 为传播时延，C 为瓶颈带宽，I 为 inflight 余量。新模型如下：

$\frac{d{x}_1}{dt}=\frac{\frac{y_1+I}{z}}{z+R}-x_1$ 【E = bw / delay 表达式】
$\frac{d{w}_1}{dt}=C\cdot \frac{y_1+I}{y_1+y_2+y_3+3\cdot I}-w_1$ 【w = C * buffer 占率】
$\frac{d{y}_1}{dt}=w_1\cdot z-y_1$ 【bdp 公式】

$\frac{d{x}_2}{dt}=\frac{\frac{y_2+I}{z}}{z+R}-x_2$
$\frac{d{w}_2}{dt}=C\cdot \frac{y_2+I}{y_1+y_2+y_3+3\cdot I}-w_2$
$\frac{d{y}_2}{dt}=w_2\cdot z-y_2$

$\frac{d{x}_3}{dt}=\frac{\frac{y_3+I}{z}}{z+R}-x_3$
$\frac{d{w}_3}{dt}=C\cdot \frac{y_3+I}{y_1+y_2+y_3+3\cdot I}-w_3$
$\frac{d{y}_3}{dt}=w_3\cdot z-y_3$

$\frac{dz}{dt}=\beta \cdot (y_1+y_2+y_3)$ 【排队时延与 buffer 占率成正比】

基于这个新模型，我准备用欧拉法求数值解，先简单介绍欧拉法。给定一个常微分方程：

$\frac{dy}{dt}=f(t,y),y(t_0)=y_0$

欧拉法求解步骤如下：

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

我一开始硬编码了一个 python matplotlib 实现，非常不雅，如果涉及共享瓶颈的流太多，欧拉法的写法将非常令人烦躁，所以不得不重构代码。

但我真不会编程，稍微会一点，就这个代码吭哧折腾快一下午，但还是没有做到可扩展：

#!/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 * (...)
...

假努力，工整而恶心。

浙江温州皮鞋湿，下雨进水不会胖。