以下的代码绘制了 bbr 的收敛相图：

#!/opt/homebrew/bin/python3

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

def model(vars, t, C, g):
    x, y = vars
    dxdt = C * (g * x) / (g * x + y) - x
    dydt = C * (g * y) / (g * y + x) - y
    # 下面是 inflight 守恒算法的模型代码
    #dxdt = C * (x + g) / (x + y + g) - x
    #dydt = C * (y + g) / (y + x + g) - y
    return [dxdt, dydt]

def curve_length(trajectory, cycles):
    length = 0
    for i in range(1, cycles):
        x1, y1 = trajectory[i - 1]
        x2, y2 = trajectory[i]
        dx = x2 - x1
        dy = y2 - y1
        length += np.sqrt(dx**2 + dy**2)
    return length

def convergence_time(trajectory):
    length = 0
    for i in range(1, len(trajectory)):
        x1, y1 = trajectory[i - 1]
        x2, y2 = trajectory[i]
        if np.abs(x1 - x2) < 0.001 and np.abs(y1 - y2) < 0.001:
            return i
    return len(trajectory)

C = 10.0
g_values = [0.6, 1.25, 2.25, 5]
# 下面是 inflight 守恒算法的 I 参数，为了不改 model 代码，仍用 g_values 名称
#g_values = [0.8, 2, 5, 10, 20]
initial_values = [(2, 8), (3, 7), (6, 4)]

t = np.linspace(0, 100, 1000)
plt.figure(figsize=(8, 8))
for g in g_values:
    for x0, y0 in initial_values:
        solution = odeint(model, [x0, y0], t, args=(C, g))
        trajectory = [(solution[i, 0], solution[i, 1]) for i in range(len(solution))]
        cycles = convergence_time(trajectory)
        length = curve_length(trajectory, cycles)
        print(cycles, length, g)
        if g < 1:
            plt.plot(solution[:, 0], solution[:, 1], label=f'g={g}', linestyle = 'dashed')
        else:
            plt.plot(solution[:, 0], solution[:, 1], label=f'g={g}')


x = np.linspace(0, 100, 1000)
plt.plot(x, x, label='x = y')
plt.plot(x, C - x, label='x + y = C')
plt.axis('equal')

for x0, y0 in initial_values:
    plt.annotate(f'({x0}, {y0})', (x0, y0), textcoords="offset points", xytext=(0,10), ha='center')

plt.xlabel('x')
plt.ylabel('y')
x_min, x_max = plt.xlim()
y_min, y_max = plt.ylim()
plt.xlim(0, 10 + 2)
plt.ylim(0, 10 + 2)
plt.title(f'bbr convergence, C = {C}')
plt.legend()
plt.grid()
plt.show()

绘图如下：

虚线表示 g < 1 的场景，而 g = 1 时就是 x +y = C 本身，可见 g <=1 时，系统都不会收敛，当 g > 1 时，收敛速度与 g 正相关。

代码中的 convergence_time 和 curve_length 函数可以定量计算收敛速度：

• convergence_time，当 x，y 均不再变化时，步骤数越少即用时越短，收敛越快；
• curve_time，相同或更短的时间内，相图轨迹 “跑” 得越远，收敛越快。

我们需要从定量计算中总结出定性规律，企图可以通过相图的几何特征中，即相图轨迹的走势与公平线夹角越小，收敛越快，相图轨迹越平直(而不弯曲迂回)，收敛越快。

但收敛快的代价是更大的 buffer 占用，时延增加，这又是前面谈到的效率和公平不可兼得的问题。

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