Bianchi模型、python计算及ns3验证

由于项目与学习需要，最近学习了bianchi模型，并在python中进行了公式->代码的转化，仿真结果与ns3结果对比。

本文更多的是理解模型各个部分的含义、把各个简单的推导过程转化为python、ns3对比：

1 理论吞吐与传输概率、传输成功概率、包长、速率、排队时间、成功传输时间、碰撞耗时的关系

2 传输概率、传输成功概率、包长、速率、排队时间、成功传输时间、碰撞耗时的关系这几个参数的简单推导

3 各个公式都有最初的已知量、推导出的中间量、最终的目标吞吐，一级一级的整理就可获得最初已知量和目标吞吐的关系，并将其画图

4 ns3在无实测环境的情况下，是一个比较好的验证环境，对比第三步的图以及仿真结果图，可以验证模型是否符合实测情况

而事实上，建模的标准流程就是以上步骤：根据初始变量、各种场景，列出基本的推导->根据中间量导出吞吐或者延时->算出各个情况下吞吐、延时结果，描点画出折线->ns3设置同样的条件，仿真对应的延时吞吐结果，描点并证明自己的模型可靠

分为3部分

本文为第一部分

第二部分是Bianchi模型、python计算及ns3验证_关于E[P*]的补充-CSDN博客

第三部分是验证2~10 STA的情况下，bianchi和ns3对比结果Bianchi模型、python计算及ns3验证_关于2~10 STA验证的补充-CSDN博客

bianchi模型简要说明

来源：

https://zhuanlan.zhihu.com/p/24302361https://zhuanlan.zhihu.com/p/24302361

结合2步骤中的python代码看公式：

预设

self.bitrate = 56E6#传输速率，如连续传输56Mbit，需要1s
self.n = 30#STA数目self.b00 = 0#马尔科夫链回到 0 0状态的概率self.ACK = 14 * 8#ack包长度 bit
self.SIFS = 10E-6#sifs时间 10us
self.slot = 9E-6#9us
self.DIFS = 3 * 9E-6#difs时间 3*slot
self.H = 2 * 8#头长度 bitself.E_P = 80#平均包长
self.E_P_star = 100#碰撞包最长平均包长self.W = 12#初始化窗口大小
self.m = 7#最大退避轮次self.prop_delay = 0#传播时延

 在开始计算前再初始化一次：

bitrate = 1E6
ACK = 112 + 128
SIFS = 28E-6
slot = 50E-6
DIFS = 128E-6
E_P = 8184
E_P_star = E_P
WW = [32, 128]
mm = [3, 5]
H = 272 + 128
prop_delay = 0nn = range(5, 50)

self.calculate_p_t()

计算冲突概率p和某个STA的发送概率τ：

def calculate_p_t(self):def equations(x):p, t = xmy_sum = 0.0for i in range(0, self.m):my_sum += (2.0 * p) ** i
##p - 1.0 + (1.0 - t) ** (self.n - 1.0),
##2.0 / (1.0 + self.W + p * self.W * my_sum) - t return ( p - 1.0 + (1.0 - t) ** (self.n - 1.0), 2.0 / (1.0 + self.W + p * self.W * my_sum) - t )self.p, self.t = fsolve(equations, (0.1, 0.1))print("System solved, error (p, tau): " + str(equations((self.p, self.t))))

其中起始窗口Wmin、最大重试次数m、STA数目n都是已知的。

my_sum 是Σ求和结果，calculate_p_t联立以上公式并求出p与t的结果——上式只有p和τ两个未知数，可以解出结果。

self.calculate_Ptr()

def calculate_Ptr(self):self.P_tr = 1.0 - (1.0 - self.t) ** self.n

 

计算n个STA，除了大家都没尝试传输以外的情况——有一个或多个STA尝试传输——其中包含传输成功和碰撞两种情况。

中间值已知

self.calculate_Ps()

def calculate_Ps(self):self.Ps = self.n * self.t * (1.0 - self.t) ** (self.n - 1) / self.P_tr

 计算在有STA传输的情况下，只有一个STA尝试传输而其他STA都没尝试传输的概率，共有n个STA，所以是n倍。

self.calculate_Ts()

    def calculate_Ts(self):self.T_s = self.H / self.bitrate + self.E_P / self.bitrate + self.SIFS + self.prop_delay + self.ACK / self.bitrate + self.DIFS + self.prop_delay

self.calculate_Tc()

def calculate_Tc(self):self.T_c = self.H / self.bitrate + self.E_P_star / self.bitrate + self.DIFS + self.prop_delay

 

Ts是成功传输耗时，Tc是碰撞传输耗时，其中H是头用时，EP是负载也就是数据部分用时，ack是回ACK用时，SIFS是负载和ack之间的等待延时，DIFS是传输成功或传输收不到ack之后的等待延时 、δ是传播时延——不过事实上SIFS和DIFS本来就是等传播时延的时间，传播时延含了，所以self.prop_delay = 0

其中Tc的计算Ep*是碰撞包的最长包的平均值，如果假设所有的包长都相等，则Ep=Ep*，我们这里取值（如不假设包长相等，参看Bianchi模型、python计算及ns3验证_关于E[P*]的补充-CSDN博客文章浏览阅读21次。bianchi的原文，在包长都选一样的情况下，P=E[P]=E[P*]，也就是说正常传输、碰撞的payload耗时都是一样的——不过bianchi分析了如果引入不同的长度会造成何种影响：k个STA包碰撞的情况下，决定E[P*]的是k个STA中时间占用最长的那个,事实上此公式很让人不解，左侧需要求值是平均最大包长，右边分母是概率，分子左半边也是概率，分母右半边如果按照PDF看的话也是概率，如果按照CDF看就是概率的积分，再乘以1作为取值，就合理了。https://blog.csdn.net/Mr_liu_666/article/details/142748859）

        self.E_P = 80self.E_P_star = 100

Ep  Ep* H ACK的长度都已知，除以已知的bitrate就是对应时间，δ DIFS SIFS都已知

self.calculate_S()

def calculate_S(self):self.S = self.Ps * self.P_tr * (self.E_P / self.bitrate) / ((1.0 - self.P_tr) * self.slot + self.P_tr * self.Ps * self.T_s + self.P_tr * (1.0 - self.Ps) * self.T_c)

 

中间值都已知

self.calculate_b00()

    def calculate_b00(self):self.b00 = 2.0 * (1.0 - 2.0 * self.p) * (1.0 - self.p) / ((1.0 - 2.0 * self.p) * (self.W + 1) + self.p * self.W * (1.0 - (2.0 * self.p)**self.m))

 

初始窗口W、最大退避重试次数m已知，中间量已知

python公式化

来源：

https://github.com/segmentation-fault/BianchiPy/blob/master/calculate_bianchi_tr.pyhttps://github.com/segmentation-fault/BianchiPy/blob/master/calculate_bianchi_tr.py

源码内容：

__author__ = 'antonio franco''''
Copyright (C) 2018  Antonio Franco (antonio_franco@live.it)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
'''from scipy.optimize import fsolve
import matplotlib.pyplot as pltclass Bianchi:# Calculates the throughput of a saturated IEEE 802.11 WLAN basic scheme according to:# G.Bianchi, "Performance analysis of the IEEE 802.11 distributed coordination function," in IEEE# Journal on Selected Areas in Communications, vol. 18, no. 3, pp. 535 - 547, March 2000.# doi: 10.1109 / 49.840210def __init__(self, bitrate, n, ACK, SIFS, slot, DIFS, E_P, E_P_star, W, m, H, prop_delay):# INPUT:# bitrate: raw bitrate in bps# n: number of STAs# ACK: ACK length in bits# SIFS: SIFS duration in seconds# slot: slot duration in seconds# DIFS: DIFS duration in seconds# E_P: average packet payload size in bits# E_P_star:  average  length  of  the  longest packet payload involved in a collision in bit (for an example, see eq. 16 in the paper)# W: Minimum contention window size in slots# m: Retry limit# H: Header size in bits# prop_delay: Propagation delay in seconds## OUTPUT:# S: normalized system throughput, defined as the fraction of time the channel is used to successfully transmit payload bits.# p: the probability of a collision seen by a packet being transmitted on the channel# t: the probability that a station transmits in a randomly chosen slot time# Ps: the probability that a transmission occurring on the channel is successful is given by the probability that exactly one station transmits on the channel, conditioned on the fact that at least one station transmits# P_tr: the probability that there is at least one transmission in the considered slot time# T_s: the average time the channel is sensed busy, in seconds# T_c: the average time the channel is sensed busy by each station during a collision in seconds.# independent varzself.bitrate = 56E6self.n = 30self.b00 = 0self.ACK = 14 * 8self.SIFS = 10E-6self.slot = 9E-6self.DIFS = 3 * 9E-6self.H = 2 * 8self.E_P = 80self.E_P_star = 100self.W = 12self.m = 7self.prop_delay = 0# dependent varzself.p = 0self.t = 0self.Ps = 0self.P_tr = 0self.T_s = 0self.T_c = 0self.S = 0self.bitrate = bitrateself.n = nself.ACK = ACKself.SIFS = SIFSself.slot = slotself.DIFS = DIFSself.E_P = E_Pself.E_P_star = E_P_starself.W = Wself.m = mself.H = Hself.prop_delay = prop_delayself.calculate_p_t()self.calculate_Ptr()self.calculate_Ps()self.calculate_Ts()self.calculate_Tc()self.calculate_S()self.calculate_b00()def calculate_b00(self):self.b00 = 2.0 * (1.0 - 2.0 * self.p) * (1.0 - self.p) / ((1.0 - 2.0 * self.p) * (self.W + 1) + self.p * self.W * (1.0 - (2.0 * self.p)**self.m))def calculate_b(self, i, k):if i < self.m:W_i = 2.0 ** i * self.Wb_i0 = self.p ** i * self.b00else:W_i = 2.0 ** self.m * self.Wb_i0 = self.p ** self.m / (1 - self.p) * self.b00return (W_i - k) / W_i * b_i0def check_p_t(self):c1 = self.p - 1.0 + (1.0 - self.t) ** (self.n - 1.0) <= 1.49012e-08my_sum = 0.0for i in range(0, self.m):my_sum += (2.0 * self.p) ** ic2 = 2.0 / (1.0 + self.W + self.p * self.W * my_sum) - self.t <= 1.49012e-08return c1 and c2def calculate_p_t(self):def equations(x):p, t