接着昨天的继续说，参见 inflight 守恒建模。

欧拉数值解看起来不够优雅，所以我打算找个别的方式试一下，顺便学一下 python，我不会编程，但也不是一点也不会，我稍微会一点，所以想进一步学习一点。

但我连 python 环境都都没，于是开始装环境，python3 装好了之后 pip 不能用，折腾 pip 好久，发现 pip 很慢且频繁 timeout，于是改了 source 为:

# 进入 ~/.pip 目录，编辑 pip.conf 文件，没有就新建 
[global] 
timeout = 6000 
index-url = https://pypi.tuna.tsinghua.edu.cn/simple
trusted-host = pypi.tuna.tsinghua.edu.cn

快了很多。

开始装 matplotlib，numpy … 我之前一直使用在线版本 https://matplotlib.online/，也挺好用。
最后发现还要装 scipy：

pip install scipy

然后开始学习 scipy，求助 chatgpt，很快就上道了，昨天的 E_best 模型的 z 表达式更新如下：

$z=\frac{numflows\cdot I}{C}+微小扰动$

更合理的应该是：

$z=\frac{numflows\cdot I}{C\cdot \frac{C}{tsum}}$

其中 tsum 为所有 flow 的发送速率之和，即所有 y 的求和。

接下来所有问题就归结为一个问题，用负反馈压住 I 即可，确保它时收缩的，这就是最佳的拥塞控制。

花了好几个小时编了下面的代码，连垃圾也没顾得上扔：

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

class E_best_model:
    def __init__(self, R, C, I):
        self.R = R
        self.C = C
        self.I = I
        self.last_z = None

    def tcp_e_best(self, state, t, C, R, I, num_pairs):
        dstate_dt = []

        tsum = 0
        if self.last_z is None:
            z = state[num_pairs*2]
        else:
            z = self.last_z

        for i in range(num_pairs):
            xi = state[2*i]
            yi = state[2*i + 1]

            dxi_dt = yi / (z + R) - xi
            dyi_dt = C * (yi * z + I) / (sum(yj * z + I for yj in state[1::2])) - yi

            tsum += yi
            dstate_dt.extend([dxi_dt, dyi_dt])

        dz_dt = 0
        self.last_z = num_pairs*I / tsum
        dstate_dt.append(dz_dt)

        return dstate_dt

# 参数
R = 2.0  # 传播时延
C = 18.0  # 瓶颈带宽
I = 1.0  # 余量，用于公平收敛，没有余量就不会收敛

# 初始条件
num_pairs = 3
initial_state = [0.2, 4, 0.9, 10, 0.1, 5, 1]  # x1, y1, x2, y2, x3, y3, z 的初始值

t = np.linspace(0, 200, 2000)

model_instance = E_best_model(R, C, I)

solution = odeint(model_instance.tcp_e_best, initial_state, t, args=(C, R, I, num_pairs))

x_solution = solution[:, :num_pairs]
y_solution = solution[:, num_pairs:-1]
z_solution = solution[:, -1]

plt.figure(figsize=(12, 6))
for i in range(num_pairs):
    plt.subplot(2, 1, 1)
    plt.plot(t, x_solution[:, i], label=f'x{i+1}(t)')
    plt.plot(t, y_solution[:, i], label=f'y{i+1}(t)')
    plt.legend()

plt.subplot(2, 1, 2)
plt.plot(t, z_solution, label='z(t)')
plt.legend()

plt.tight_layout()
plt.show()

这个代码有几个要点：

• 涉及 n 条流共存时，我不想一下子写很多微分方程，所以要用 num_pairs + for loop 来迭代；
• 关于 z 的赋值，每次的计算都要用上一次的 z 值，所以要保存一个全局变量，我用 class 封装；
• R，C，I 参数，本应该封装在 class 里，但我好像写得比较乱，封装在 class 里的没有用，用了传参；
• 依旧没有负反馈，因为代码写不好…

挺好，挺好，效果如下：

我也终于明白了 “人生苦短，我用 python” 的意义，真的是太好用了，只可惜我无能，编程编得不好，有太多东西要学习了。

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