用 python scipy 库模拟拥塞控制模型

接着昨天的继续说，参见 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=\dfrac{numflows\cdot I}{C}+微小扰动$

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

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

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as pltclass E_best_model:def __init__(self, R, C, I):self.R = Rself.C = Cself.I = Iself.last_z = Nonedef tcp_e_best(self, state, t, C, R, I, num_pairs):dstate_dt = []tsum = 0if self.last_z is None:z = state[num_pairs*2]else:z = self.last_zfor i in range(num_pairs):xi = state[2*i]yi = state[2*i + 1]dxi_dt = yi / (z + R) - xidyi_dt = C * (yi * z + I) / (sum(yj * z + I for yj in state[1::2])) - yitsum += yidstate_dt.extend([dxi_dt, dyi_dt])dz_dt = 0self.last_z = num_pairs*I / tsumdstate_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()