机场用户选择套餐的数学建模

前段时间,有个机场主找到我,希望我提供数学上的帮助,帮助他们计算用户会选择什么套餐。我回复,需要必要的数据来建立数学模型。但他们拒绝给我提供数据,可能原因是没有收集或者脱敏麻烦。

于是,经过协商,我为他们建立不需要数据、且可以根据数据调整的模型,并且顺带给出一些降低成本的方案;他们给我更高的报酬。

尽管我强烈怀疑他们能够通过这个模型获得到的利润需要多久才能填补聘用我的费用,但我还是给他们做了。并且,他们允许我把这个模型部分开源出来。

他们的机场用的是 PHP,但我为了方便计算,用的是 pytorch,我想他们似乎还需要找个会 Python 的来跑起来这个模型。

模型参数

假设用户的预算服从的分布的概率密度函数为b(x1)b(x_1),流量需求服从的分布的概率密度函数为t(x2)t(x_2),节点质量需求服从的概率密度函数为q(x3)q(x_3)

假设用户对以上 3 种因素侧重比例的概率密度函数为a(y1,y2)a(y_1,y_2),定义为:

  • y1(0,1)y_1\in(0,1),y2(0,1)y_2\in(0,1)y1+y2<1y_1+y_2<1,否则a(y1,y2)=0a(y_1,y_2)=0
  • y1y_1为预算的权重,y2y_2为流量的权重,1y1y21-y_1-y_2为质量的权重

每种套餐都有价格、流量、质量 3 个属性。假设当用户思考买哪一种套餐时,会按照层次分析法决策。上述 bb, tt, qq 的参数需要参数估计。

P=(pr1,pr2,pr3,,prn)P=(p_{r1},p_{r2},p_{r3},\cdots,p_{rn})', T=(t1,t2,t3,,tn)T=(t_1,t_2,t_3,\cdots,t_n)', Q=(q1,q2,q3,,qn)Q=(q_1,q_2,q_3,\cdots,q_n)' 为产品数组向量,分别为产品的价格、流量、质量。产品数组向量由产品性质决定,不需要参数估计。

import torch

# 定义套餐的价格、流量、质量
packageNumber = 3
packagePrice = [1, 2, 3]
packageTrafficLimit = [1, 2, 3]
packageQuality = [1, 2, 3]

层次分析法量化

对于预算的满意度,可以定义为f(pr,x1)=10prx1f(p_r,x_1)=10^{\displaystyle-\frac{p_r}{x_1}},其中prp_r为套餐价格,x1x_1为预算

对于流量的满意度,可以定义为g(t,x2)=4tx234g(t,x_2)=\sqrt{\displaystyle\frac{4t}{x_2}-\frac{3}{4}}(当根号下为 0 以下时定义为 0),其中tt为套餐流量,x2x_2为套餐需求

对于节点质量的满意度,可以定义为h(q,x3)=qx33h(q,x_3)=\sqrt[3]{\displaystyle\frac{q}{x_3}},其中qq为节点质量量化指标,x3x_3为节点质量需求。

用户会选择满意度

y110prx1+y24tx234+y3qx33y_1\cdot10^{\displaystyle-\frac{p_r}{x_1}}+y_2\sqrt{\displaystyle\frac{4t}{x_2}-\frac{3}{4}} + y_3\sqrt[3]{\displaystyle\frac{q}{x_3}}

最大的套餐

并且,如果对套餐的满意度都小于 0.1,那么他将不会选择这个套餐。

写成向量形式,则有:

SATISFY=[f(P,x1),g(T,x2),h(Q,x3)]Y\text{SATISFY} = [f(P,x_1),g(T,x_2),h(Q,x_3)]'Y SELECT={maxSATISFYPACKAGE, SATISFY>0.1, SATISFY0.1SELECT=\begin{cases} \max_{\text{SATISFY}} \text{PACKAGE},\ \text{SATISFY}>0.1\\ \varnothing,\ \text{SATISFY}\leq 0.1 \end{cases}

对应的 Python 代码为

# 定义用户参数与套餐到满意度的映射
def user_satisfaction(vector_x, vector_y, pr, t, q):
    # 确保 X 和 Y 是 torch 张量
    vector_x = torch.tensor(vector_x, dtype=torch.float32)
    vector_y = torch.tensor(vector_y, dtype=torch.float32)
    
    term1 = vector_y[:, 0] * 10 ** (-pr / vector_x[:, 0])
    term2 = vector_y[:, 1] * torch.sqrt(4 * t / vector_x[:, 1] - 0.75)
    term3 = vector_y[:, 2] * (q / vector_x[:, 2]) ** (1 / 3)
    
    return term1 + term2 + term3

用户建模

用户预算分布

预算可以认为大致是服从帕累托分布的,假设初始参数为x1Par(1,0.6)x_1\sim Par(1,0.6),即

b(x1)=0.6x11.6(x1>1)b(x_1) = \frac{0.6}{x_1^{1.6}} (x_1>1)

其中 0.6 需要参数估计以获得准确值。

帕累托分布通常用来描述不平等的分布,例如财富或收入,以及对应的预算。在现实世界中,财富分布往往是高度不平等的,少数人拥有大部分财富。因此,为了反映用户预算,使用帕累托分布应该是合适的。

此外,帕累托分布与所谓的“80/20 规则”密切相关,即 80%的效果(如消费)来自 20%的原因(如消费者)。在许多经济模型中,这种现象是普遍存在的,比如少数消费者贡献了大部分的消费支出。

帕累托分布的两个参数中,第一个参数为截止值,当x1<1x_1<1时,b(x1)=0b(x_1)=0。截止值无法参数估计,但假设为 1 是合适的。

流量需求分布

流量需求可以认为近似服从正态分布,参数为x2N(130,302)x_2\sim N(130,30^2),即

t(x2)=1402πe12(x212040)2t(x_2)=\frac{1}{40\sqrt{2\pi}}e^{-\displaystyle\frac{1}{2}\left(\frac{x_2-120}{40}\right)^2}

这两个参数都需要参数估计。

质量需求分布

如果量化质量为1101\sim10,其中 1 为完全灵车机场,10 为全 IPLC 多入口智能解析机场(性能接近游戏加速器),那么质量需求也可以认为是服从帕累托分布的,参数为x3Par(1,2)x_3\sim Par(1,2),即

q(x3)=2x32q(x_3)=\frac{2}{x_3^2}

这个帕累托分布同样只有第二个参数需要参数估计。

权重分布

注意到y1y_1, y2y_2, y3y_3的分布实际上被均匀地放在了一个等边三角形上。这个等边三角形是一个三棱锥的底面。这个三棱锥的 3 个侧面都是直角边边上为 1 的直角三角形。

等边三角形的 3 个顶点分别表示其中一个权重为 1 而另外两个权重为 0 的情况。

在使用(y1,y2)(y_1,y_2)表示时,3 种情况分别对应(0,0)(0,0), (1,0)(1,0), (0,1)(0,1)。而为了让y1,y2,y3y_1,y_2,y_3变成更加均匀的等边三角形,需要将这个区域线性变换为(3+1,0)(-\sqrt{3}+1,0), (31,1)(\sqrt{3}-1,1), (31,1)(\sqrt{3}-1,-1)。这个线性变换很容易求得。

首先把(3+1,0)(-\sqrt{3}+1,0)(0,0)(0,0)对齐,那么其他两个点是(3,1)(\sqrt{3},1), (3,1)(\sqrt{3},-1)。于是,变换矩阵为

[3311][1001]1=[3311]\begin{bmatrix} \sqrt{3} & \sqrt{3} \\ 1 & -1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}^{-1}= \begin{bmatrix} \sqrt{3} & \sqrt{3} \\ 1 & -1 \end{bmatrix}

平移向量为(1,0)(-1,0),即

[y1y2]=[3311][y1y2]+[10]\begin{bmatrix} y_1' \\ y_2' \end{bmatrix}= \begin{bmatrix} \sqrt{3} & \sqrt{3} \\ 1 & -1 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}+ \begin{bmatrix} -1 \\ 0 \end{bmatrix}

不过,对于 3 种因素的侧重,可以认为一般人根本不在乎节点质量,而只考虑节点的价格和流量。所以,可以选择截半边的二元正态分布。截取之前的参数为

(y1,y2)N([1/21/2],[1/9001/9])(y_1,y_2)\sim \mathcal{N}\left( \begin{bmatrix} 1/2 \\ 1/2 \end{bmatrix}, \begin{bmatrix} 1/9 && 0 \\ 0 && 1/9 \end{bmatrix} \right)

即概率密度函数为

a(y1,y2)=9πe92[(x11/2)2+(x21/2)2]a(y_1,y_2)= \frac{9}{\pi}e^{\displaystyle -\frac{9}{2}\left[(x_1-1/2)^2+(x_2-1/2)^2\right]}

这两个参数应该也是不需要参数估计的。

from torch.distributions import Pareto, Normal
from torch.distributions.multivariate_normal import MultivariateNormal

# 定义用户预算分布
user_budget_distribution = Pareto(1, 0.6)
# 定义用户流量需求分布
user_traffic_demand_distribution = Normal(130, 30)
# 定义用户质量需求分布
user_quality_demand_distribution = Pareto(1, 2)
# 定义用户侧重分布
user_weight_mean = torch.tensor([0.5, 0.5])
user_weight_covariance = torch.tensor([[1 / 9, 0], [0, 1 / 9]])
user_weight_distribution = MultivariateNormal(user_weight_mean, user_weight_covariance)

# 生成样本
sample_size = 10000     # 样本数量
user_budget_sample = user_budget_distribution.sample((sample_size,))
user_traffic_demand_sample = user_traffic_demand_distribution.sample((sample_size,))
user_quality_demand_sample = user_quality_demand_distribution.sample((sample_size,))
user_weight_sample = user_weight_distribution.sample((sample_size,))

# 模拟用户选择套餐
y_1 = user_weight_sample[:, 0]
y_2 = user_weight_sample[:, 1]
y_3 = 1 - y_1 - y_2
user_sample = torch.stack([user_budget_sample, user_traffic_demand_sample, user_quality_demand_sample], dim=1)
user_weight_sample = torch.stack([y_1, y_2, y_3], dim=1)
user_satisfaction_of_packages = []
for i in range(packageNumber):
    pr = packagePrice[i]
    t = packageTrafficLimit[i]
    q = packageQuality[i]
    satisfactions = user_satisfaction(user_sample, user_weight_sample, pr, t, q)
    user_satisfaction_of_packages.append(satisfactions)
    
user_satisfaction_of_packages = torch.stack(user_satisfaction_of_packages, dim=1)
max_values, max_indices = torch.max(user_satisfaction_of_packages, dim=1)
indices = torch.where(max_values < 0.1, torch.tensor(-1), max_indices)

参数估计

综上,需要估计的参数有:

  • 用户预算的帕累托分布kk参数:kbk_b
  • 流量需求分布的均值和方差:μt\mu_t,σt2\sigma^2_t
  • 质量需求的帕累托分布kk参数:kqk_q

似乎是可以通过神经网络估计出这几个参数,不过我没去研究,初始参数又不是不能用。

最小化成本

参数定义

假设用户的节点偏好为

R=[r1r2rn]R = \begin{bmatrix} r_1 \\ r_2 \\ \vdots \\ r_n \end{bmatrix}

且规定i=1nri=1\sum_{i=1}^n r_i = 1

用户在各节点实际使用的流量为

U=[u1u2un]=uRU= \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix} = uR

其中,uu为用户实际使用的总流量。

这里假设用户不会因为差点用超流量而改用不喜欢的节点。

不考虑每个节点能够承载的用户数量限制,假设用户会先用完节点的流量而不是因为请求导致节点超载。于是,假设各个节点的边际成本为

C=[c1c2cn]C=\begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{bmatrix}

则成本为CU=uCRC'U=uC'R

这里假设了边际成本为常数,但一般机场是规模经济(边际成本随用量降低)的。

若定价为prp_r,则边际利润为pruCRp_r-uC'R

于是,得到利润数据,平均每个套餐的利润为

rˉ=pruˉCR\bar{r} = p_r - \bar{u}C'R

超售问题

为了确定节点的流量边际成本,需要计算超售比,因为cn=ctn/Rosc_n=c_{tn}/R_{os},其中ctnc_{tn}为实际流量边际成本,RosR_{os}为超售比。

为了确保不会因为超售太多而导致可用性下降,定义RamaxRa_{max}为超售可用性,意义为:有RamaxRa_{max}的概率不会因为超售而不可用。

这里定义超参数Ramax=1104=99.99%Ra_{max}=1-10^{-4}=99.99\%

对于套餐g1,g2,,gng_1,g_2,\cdots,g_n,设用户实际使用总流量的概率密度函数为fn(u)f_{n}(u),则对于某一个节点,使用的流量的概率密度函数为rnfn(urn)r_nf_n(ur_n)

设节点实际的流量为trt_r,则为满足可用性,需要

Pr{u<tr}RamaxPr\{u<t_r\} \geq Ra_{max}

n0trrnfn(urn)durn=n0tr/rnfn(u)duRamax\sum_n\int_0^{t_r} r_nf_n(ur_n) dur_n=\sum_n\int_0^{t_r/r_n}f_n(u)du\geq Ra_{max}

例如,假设只有1个套餐,这个套餐的限额为150G,对于日本节点,倍率为1,用户的偏好程度为0.5,使用的流量服从正态分布uN(80,100)u\sim N(80,100),那么,trt_r应该为58.59558.595,超售比为2.54232.5423

用户使用的流量在与之前参数估计得到的用户需求流量大致相同,但有一定误差。这是因为,用户往往会高估自己需要使用的流量。不过,假设用户实际使用的流量等于用户流量需求可能是可行的,因为留一点流量上的冗余是问题不大的。

因为

Ros=trnrntR_{os}=\frac{t_r}{\sum_n r_nt}

trt_r可以求解方程得到,

所以RosR_{os}可以由数据分析得到。

最大期望利润问题

最大期望利润问题,也就是使得

nfnrnˉ=fn(pr,nunˉCR)\sum_n f_n\bar{r_n} =f_n\left( p_{r,n} - \bar{u_n}C'R\right)

最大的问题。(其中fnf_n为选择率)

由于CCRR都是基本不变的量,而且在参数估计出用户的需求分布后,fnf_n只取决于套餐的属性,unˉ\bar{u_n}也是稳定的,所以最大利润问题就是确定套餐属性的问题。

如果假设套餐的质量属性是常数,那么变量就只有各个套餐的价格和流量了。

所以,计算最大期望利润的方法为:

  • 统计得到CRC'R
  • 确定各个套餐的质量
  • 确定各个套餐的价格或者流量,把另一个设置为变量。
  • 对变量进行梯度下降,使上面那个和式取最大值
    • fnf_n可以在有用户的数学模型之后,蒙特卡罗方法得到
    • prp_r为价格,是输入变量
    • unu_n与用户的建模相关