在学术圈里有时候最磨人的不是“想不出模型”而是“模型摆在面前却跑不出结果”。供应链领域尤其如此一个供应商带两个零售商的博弈模型看似只是把经典的单渠道问题往双渠道推了一步真动手用Mathematica复现的时候才会发现里面全是细节逆需求函数怎么写、批发价和零售价谁先动、符号解和数值解在什么条件下互相打架、利润函数在什么情况下不再是凹函数……这些问题不踩一遍光靠看论文很难有体感。这篇文章我想复盘一次完整的论文复现过程。目标很明确把“一个供应商 两个零售商”的经典二级供应链模型在Mathematica里从建模到求解全部跑通包括Stackelberg博弈的逆向归纳、一阶条件的符号推导、均衡解的数值验证以及那些论文里通常不会写、但实际复现时一定会遇到的坑。适合正在做供应链方向研究、准备用Mathematica做实验、或者想快速验证某篇论文结果的同学参考。1. 模型结构与博弈逻辑先把“谁决策”想清楚1.1 一个供应商两个零售商的典型设定这类模型的供应链层级不复杂上游是唯一的供应商制造商用批发价 (w) 向两个零售商供货下游是两个零售商分别面对自己的终端市场决定零售价或订货量。经典文献里通常采用线性逆需求函数[ p_i a - q_i - d q_j, \quad i,j \in {1,2}, i \neq j ]其中 (a) 是市场规模参数(q_i) 是零售商 (i) 的销量(d) 表示两个零售商之间的产品替代程度(0 \leq d 1)越接近1说明产品同质化越高。这个函数的意思很直观零售商 (i) 的价格不仅取决于自己的销量还会被对手的销量压制因为两个零售商卖的是有一定替代性的产品消费者会互相分流。我和很多刚开始复现这类模型的人一样最容易在第一眼就忽略一个关键点这个模型里的决策顺序不是同时的而是分层的。供应商先定批发价两个零售商看到批发价之后再各自决定自己的零售价或订货量。这就是典型的Stackelberg博弈而且是一个领导者供应商面对两个跟随者零售商的博弈结构。决策时序拆开看是这样的供应商基于对零售商反应行为的预判选择批发价 (w)目标是最大化自己的利润 (\pi_s w(q_1 q_2))。两个零售商观察到 (w) 后分别选择 (q_i) 最大化自己的利润 (\pi_{r_i} (p_i - w)q_i)。因为两个零售商在同一层、同时决策所以零售商之间存在一个Nash博弈而供应商在更上层需要用“零售商会对我的批发价做出什么反应”来倒推自己的最优决策。这种“上层领导者 下层Nash博弈”的结构在博弈论里叫双层博弈在供应链管理里就是最常见的分散决策场景。复现论文的第一步不是写代码而是把这个决策顺序用数学语言固化下来否则后面求导、代入解析解的时候顺序一乱结果全乱。1.2 逆向归纳法为什么从零售商开始推Stackelberg博弈的标准解法是逆向归纳Backward Induction。思路听起来简单——先从最后一步往回推但实际推导时很多人容易犯一个错误把零售商利润对 (q_i) 求导之后直接代进供应商的利润函数却忘了这一步求出来的是“反应函数”不是一个确定数值。具体来说零售商 (i) 的利润是[ \pi_{r_i} (a - q_i - d q_j - w) q_i ]对 (q_i) 求一阶条件[ \frac{\partial \pi_{r_i}}{\partial q_i} a - 2q_i - d q_j - w 0 ]于是得到零售商 (i) 对 (q_j) 的反应函数[ q_i \frac{a - w - d q_j}{2} ]两个零售商联立这个方程组解出来的才是均衡销量[ q_1^* q_2^* \frac{a - w}{2 d} ]注意这里得到的是关于 (w) 的表达式不是具体数字。它们才是供应商做决策时看到的“零售商的反应”。把 (q_1^, q_2^) 代进供应商利润[ \pi_s w(q_1^* q_2^*) \frac{2w(a - w)}{2 d} ]对 (w) 求一阶条件[ \frac{d\pi_s}{dw} \frac{2(a - w) - 2w}{2 d} 0 ]得到[ w^* \frac{a}{2} ]这个结果非常干净最优批发价只取决于市场规模 (a)与替代程度 (d) 完全无关。我第一次跑出来的时候也有点意外仔细想想又合理——供应商面对的是两个对称零售商的总需求而总需求 (q_1 q_2 2(a - w)/(2d)) 对 (w) 的敏感度并不受 (d) 影响所以供应商不用关心下游竞争有多激烈只管定一半的价格水平。把 (w^* a/2) 带回零售商反应函数可以得到均衡零售价和利润这一步我会在后面的Mathematica实现里完整展开。2. Mathematica 求解思路与代码实现从符号推导到数值验证2.1 为什么选Mathematica而不是手动推导这类供应链模型的代数推导其实不算难但有一个非常现实的问题当你给模型加入非对称参数比如两个零售商的成本不同、或者把需求函数改成非线性形式时手推的出错率会大幅上升而且一旦推到一半发现前面某一步符号错了回头看非常痛苦。Mathematica的核心优势在于符号计算Symbolic Computation。你只需要把利润函数写对求导、联立方程、代入消元这些机械但易错的操作全部交给内核。我复现这篇论文时几乎全部流程都是用Solve、D、Reduce这几个基础函数完成的不需要任何高级包代码量也压得很低。这里多说一句Mathematica的符号计算和MATLAB的符号工具箱不是一个量级的存在。MATLAB更偏数值符号运算只是附加功能Mathematica的内核设计就是围绕符号表达式展开的处理这种“先求反应函数、再联立求均衡、再代入求利润”的嵌套推导Mathematica的表达方式几乎和数学表达式一一对应调试心智负担小很多。2.2 核心代码完整的求解流程下面这段代码是我整理后的完整复现脚本包含了从模型设定到均衡验证的全部流程。我用的是Wolfram Language自然语法在Mathematica 12及以上版本中可以直接运行。(* 清空环境避免变量残留干扰 *) ClearAll[Global*] (* 1. 模型设定 *) (* 逆需求函数p_i a - q_i - d*q_j *) p[i_, j_] : a - q[i] - d*q[j] (* 零售商利润销售收入 - 进货成本 *) piR[i_, j_, w_] : (p[i, j] - w)*q[i] (* 供应商利润批发价格 * 总销量 *) piS[w_] : w*(q[1] q[2]) (* 2. 第二步零售商博弈Nash *) (* 对每个零售商的利润求一阶条件 *) foc1 D[piR[1, 2, w], q[1]]; foc2 D[piR[2, 1, w], q[2]]; (* 联立求解两个零售商的最优反应 *) retailerSolution Solve[{foc1 0, foc2 0}, {q[1], q[2]}] // Simplify; (* 提取均衡销量表达式 *) q1Star[w_] q[1] /. retailerSolution[[1]]; q2Star[w_] q[2] /. retailerSolution[[1]]; (* 3. 第一步供应商决策Stackelberg 领导者 *) (* 供应商利润 w * (q1* q2*) *) supplierProfit[w_] piS[w] /. {q[1] - q1Star[w], q[2] - q2Star[w]} // Simplify; (* 对 w 求一阶条件 *) focS D[supplierProfit[w], w]; (* 求最优批发价 *) wStar Solve[focS 0, w][[1, 1]] // Simplify; (* 4. 回代得到全部均衡值 *) q1Eq q1Star[w] /. wStar // Simplify; q2Eq q2Star[w] /. wStar // Simplify; p1Eq p[1, 2] /. {q[1] - q1Eq, q[2] - q2Eq} // Simplify; p2Eq p[2, 1] /. {q[1] - q1Eq, q[2] - q2Eq} // Simplify; piREq piR[1, 2, w] /. {q[1] - q1Eq, q[2] - q2Eq, w - (w /. wStar)} // Simplify; piSEq supplierProfit[w] /. wStar // Simplify; (* 输出关键结果 *) Print[最优批发价 w* , w /. wStar]; Print[零售商1均衡销量 q1* , q1Eq]; Print[零售商2均衡销量 q2* , q2Eq]; Print[零售商1均衡价格 p1* , p1Eq]; Print[零售商2均衡价格 p2* , p2Eq]; Print[零售商1均衡利润 , piREq]; Print[零售商2均衡利润 , piREq]; Print[供应商均衡利润 , piSEq]; (* 5. 数值验证 *) (* 设定参数并验证一阶条件与二阶条件 *) params {a - 100, d - 0.6}; numericW w /. wStar /. params; (* 验证零售商一阶条件在该批发价下是否为0 *) q1Num q1Eq /. params; q2Num q2Eq /. params; foc1Num (D[(a - q[1] - d*q[2] - w)*q[1], q[1]] /. {q[1] - q1Num, q[2] - q2Num, w - numericW}) /. params; Print[零售商1一阶条件数值 , foc1Num]; (* 验证供应商利润的二阶条件应小于0 *) socS D[supplierProfit[w], {w, 2}] /. wStar /. params; Print[供应商二阶条件数值 , socS];这段代码跑完之后你会得到一组非常漂亮的符号解[ w^* \frac{a}{2}, \quad q_1^* q_2^* \frac{a}{2(2d)}, \quad p_1^* p_2^* \frac{a(3d)}{2(2d)} ]再把均衡值代回利润函数[ \pi_{r_1}^* \pi_{r_2}^* \frac{a^2}{4(2d)^2}, \quad \pi_s^* \frac{a^2}{2(2d)} ]我建议你跑完代码后先别急着改参数做实验而是对着论文里的结果逐项核对这几个表达式的分母和系数。我复现的时候发现这类论文最常出现排版错误的地方就是利润表达式——分母多了个平方、或者系数从1/4写成1/2符号计算的好处就是你可以完全不依赖论文的最终公式自己推一遍就有绝对把握。2.3 代码里容易被忽略的三个细节这段代码虽然短但我实际调试时踩过不少坑分享三个最典型的细节一Solve返回的解列表顺序可能不固定。我在代码里用retailerSolution[[1]]提取第一个解但如果你换一个版本的Mathematica或者改了方程形式解的排列顺序有可能变化。稳妥的做法是加一个条件筛选比如选择解中 (q[1] 0) 的那个。处理办法可以这样写positiveSolution Select[retailerSolution, (q[1] 0 q[2] 0) /. # ];细节二变量代入顺序错位导致结果恒等于0。在计算零售商利润时如果你先把 (q_1^, q_2^) 代入利润函数再代入 (w^*)顺序反了会出现“利润表达式里还有 (w)但 (w) 已经没地方可代”的尴尬。我的建议是统一用/ .规则一次性替换不要让中间步骤里的变量分叉。细节三Simplify的时机决定结果的可读性。如果每次求导后不Simplify直接代入下一步表达式会膨胀到完全没法看。我在每个关键步骤后面都加了对齐的Simplify代价是多花一点点计算时间但结果非常干净。对于这种规模的问题计算开销完全可以忽略。3. 数值实验与结果对比验证模型行为的直觉3.1 固定参数的基准算例模型推导完毕接下来用数值实验验证均衡结果的合理性。设定市场规模 (a 100)替代系数 (d 0.6)代入公式得到变量符号表达式数值结果a100, d0.6最优批发价 (w^*)(\frac{a}{2})50零售商均衡销量 (q^*)(\frac{a}{2(2d)})19.23零售商均衡价格 (p^*)(\frac{a(3d)}{2(2d)})69.23零售商单方利润 (\pi_r^*)(\frac{a^2}{4(2d)^2})369.82供应商利润 (\pi_s^*)(\frac{a^2}{2(2d)})1923.08供应链总利润(\frac{a^2(3d)}{4(2d)^2} \frac{a^2}{2(2d)})2662.72这些数字本身并不神秘但有一个点很值得注意在这个基准算例里供应商的利润是单个零售商的5倍以上。这个差距不是偶然的。供应商作为“先动者”可以把批发价定在成本之上很高的一档而两个零售商因为互相竞争价格和销量都被压住了利润自然薄。这正是Stackelberg领导者的“先手优势”。3.2 替代系数对均衡结果的影响比基准算例更有意思的是参数敏感性分析。固定 (a 100)让 (d) 从0变化到0.9步长0.1观察各均衡值的变化趋势替代系数 d批发价 w*零售价 p*零售商销量 q*零售商利润 πr*供应商利润 πs*0.050.0075.0025.00625.002500.000.250.0071.4322.73516.532272.730.450.0068.7520.83434.032083.330.650.0066.6719.23369.821923.080.850.0065.0017.86318.881785.710.950.0064.2917.24297.261724.14这张表信息量很大。首先确认了之前的结论随着 (d) 增大批发价始终稳定在50完全不受下游竞争强度影响。零售价格持续下降从完全差异化时的75降到高度同质化时的64.29。零售商的利润几乎被腰斩从625降到297。更有意思的是供应链总利润的变化总利润从3125逐步降到2318。这意味着在分散决策下下游竞争加剧不仅伤害零售商也伤害整个供应链。原因在于双重边际效应——供应商定价时只看自己的利润不考虑这个批发价会如何扭曲零售商的定价决策。当 (d) 增大时零售商对批发价的反应变得更激烈这种扭曲被放大总利润损失就更大。这个结果给管理实践提供了很直接的启示如果上游供应商能影响零售商的竞争强度比如通过品牌区隔、产品线设计来降低替代性即使批发价不变整个供应链的利润空间也会被撑大。3.3 用Mathematica做敏感性分析的三种方式做这类参数扫描我比较推荐三种方式按适用场景排序方式一直接生成表格数据。适合快速、一次性看全貌。Table[{d, w /. wStar, (q1Eq /. {a - 100}), (p1Eq /. {a - 100}), (piREq /. {a - 100}), (piSEq /. {a - 100})} /. {d - dval}, {dval, 0, 0.9, 0.1}] // TableForm方式二用Manipulate做动态演示。适合在汇报或答辩时展示参数变化的连续趋势。我一般用Manipulate[ Plot[{p1Eq, p2Eq}, {d, 0, 1}, PlotRange - {50, 80}, AxesLabel - {d, 价格}, PlotLegends - {零售商1价格, 零售商2价格}], {a, 50, 200, 10}, {w, 20, 80, 5}]方式三求解析导数。在写论文时你需要给出“随着 (d) 增大(p^*) 单调递减”之类的严格结论那就直接求导数D[p1Eq, d] // Simplify跑出来是 (-a/(2(2d)^2))显然小于0说明零售价格确实随 (d) 递增而递减。4. 复现论文时的高频坑与排查思路4.1 符号解推出来了但数值验证对不上这是我复现这类论文时遇到最多的一个问题。符号求解非常顺利公式也漂亮但一旦代入具体参数去算和论文给出的数值结果有个别位数对不上。排查思路其实很清晰先检查是不是零售商的二阶条件不满足。这听起来很基础但在非对称模型里非常容易翻车。比如某个零售商的利润函数对 (q_i) 求二阶导[ \frac{\partial^2 \pi_{r_i}}{\partial q_i^2} -2 0 ]这个负号恒成立所以Nash均衡确实存在且唯一。但如果是非线性需求函数或者零售商之间存在成本不对称二阶条件就可能在某些参数区域失效。我的经验是每次跑完均衡解先顺手把一阶和二阶条件的数值验算写在代码里不要等到结果对不上了再回头查。另一个隐蔽的原因是参数范围越界。有些论文做了归一化处理比如把市场规模 (a) 设为1但你可能在实验里用了更大规模的参数比如 (a100)虽然模型是线性的很多比例关系不会变但如果你把两个不同量纲的参数混在一起比较就会出现看似矛盾的结论。建议复现时先把论文里的算例原封不动地跑一遍确认数值完全一致后再扩展。4.2 逆向归纳的顺序最容易错这个问题值得再说一遍Stackelberg博弈的“先定”和“后定”不代表求解顺序也是先供应商后零售商。很多初学者会先对供应商利润求导把 (w) 解出来再代到零售商的利润里——这个顺序完全反了。正确的做法是永远从产业链最下游开始推先把零售商的反应函数 (q_i(w)) 解出来再把它们代入供应商的利润函数。Mathematica里面这个顺序是通过“先算零售商问题后把结果作为规则代入供应商问题”实现的代码层面并不复杂但如果你在推导思路里没有把顺序理清就会在某个角落出现“(w) 被消掉了供应商利润变成常数”这类诡异结果。我自己的习惯是在代码开头用注释把决策时序写清楚像这样(* 博弈时序 1. 供应商选择 w 2. 零售商同时选择 q[i] 3. 利润实现 求解方式逆向归纳先解步骤2再解步骤1 *)这句话看起来多余但当你改模型加参数、或者在Mathematica里反复调试的时候它能帮你在思路混乱时迅速找回坐标。4.3 读论文公式和代码输出相差一个负号符号问题也不少见。我遇到过一次论文公式里写的均衡价格比代码输出的结果多了个负号差的怪事查了半天发现是论文里的逆需求函数定义和正文描述不一致——公式里写的是 (p a q_i d q_j)正文却说“价格随销量上升而下降”显然是排版或笔误。这种问题只有一种可靠解法不轻信论文的公式把模型设定里的函数形式和符号逻辑核对一遍。具体到这个模型里(p a - q_i - d q_j) 才能保证价格随销量下降如果是加号需求定律直接违背后面所有结果都失去经济意义。建议复现前先花5分钟从头验一遍模型的基本定性性质比出问题时瞎猜高效得多。4.4 Fast Solve 遇到多解怎么办当模型稍微扩展比如两个零售商的成本参数不同联立方程就会出现两组甚至多组解。论文里通常会直接跳过这一细节用Word排版写“解得”但你在复现时会发现Solve返回了两个解手工选错组就会得到负价格、负利润之类的荒谬结果。我的筛选策略是用物理意义约束销量、价格、利润都必须是正的。可以这样写validSolution Select[solutions, (q1 0 q2 0 p1 0 p2 0) /. # ]对于经济模型来说负均衡价格或负产量在大多数情况下意味着该解落在定义域之外可能是因为替代系数超过了某个临界值。只要把正约束加上通常能干净地去掉伪解。5. 扩展思考从复现到创新复现一篇论文的终点不应该只是让代码跑通、结果能对上。真正有价值的复现是你在跑通之后开始琢磨这个模型还能怎么改以“一个供应商 两个零售商”为起点几条自然的扩展方向值得试一试扩展一零售商之间由数量竞争改为价格竞争Bertrand竞争。把零售商的决策变量从 (q_i) 换成 (p_i)需求函数改为 (q_i a - p_i d p_j)整个推导过程的结构类似但均衡结果会不同。价格竞争通常比数量竞争更激烈零售商的利润会被进一步压缩供应商的批发价决策会不会因此受影响这个问题本身就够写一篇小论文。扩展二引入两个零售商的不对称性。比如一个零售商有线下渠道优势、另一个有线上流量优势体现在模型里就是两个零售商的需求函数截距不同(a_1 \neq a_2)或者在竞争项前的系数不同。这时供应商面对的不再是对称反应函数(w) 的最优解会发生什么样的偏移我试过把 (a_1100, a_280) 代入模型结果批发价会往强势零售商的方向倾斜但倾斜幅度并非直觉上的线性关系。扩展三加入供应商的公平关切或者风险规避。这意味着目标函数不再是简单的利润最大化而是效用最大化。比如用均值-方差效用 (U_s \pi_s - \lambda \mathrm{Var}(\pi_s)) 替换原来的利润函数这种改动在Mathematica里只需要把目标函数改一行但推导复杂度会明显上升——刚好可以检验Mathematica符号计算在更大规模问题上是否依然顺手。我在复现完基准模型后尝试了上述第一种扩展Bertrand竞争额外花了两天时间收获非常大。你会发现价格竞争下的零售商最优反应函数变成[ p_i \frac{a w d p_j}{2} ]注意这里的 (w) 出现在分子的加号位置与数量竞争下 (q_i (a - w - d q_j)/2) 的减号位置形成对比。表面上只是符号差异但代入后得到的批发价不再是 (a/2)而是和 (d) 相关的更复杂表达式。为什么价格竞争和数量竞争会产生如此本质的区别这个“为什么”比结论本身更值得写进论文里。如果你准备在这个方向深耕我建议拿到论文的第一时间先不看结论而是自己在Mathematica里把模型推一遍再对照论文结果找出分歧点。这样得到的理解深度远远超过“跑通代码、截个图”式的复现。根据我个人经验Mathematica处理供应链模型最舒服的地方在于你不需要在“推导公式”和“数值实验”之间切换两套工具。符号推完直接赋值就能跑数值画图、敏感性分析、动态演示也都在同一个环境里完成。把论文复现当成一次“思维演练”而不是“交作业”收获会大得多。最后再分享一个我常用的技巧在跑敏感分析前先用Simplify把均衡表达式的结构看清比如判断 (q_i^*) 表达式中 (a) 和 (d) 是否可分离这通常会帮你预判图表趋势省去不少可视化调参的时间。
用Mathematica复现供应链Stackelberg博弈：从建模到数值验证
在学术圈里有时候最磨人的不是“想不出模型”而是“模型摆在面前却跑不出结果”。供应链领域尤其如此一个供应商带两个零售商的博弈模型看似只是把经典的单渠道问题往双渠道推了一步真动手用Mathematica复现的时候才会发现里面全是细节逆需求函数怎么写、批发价和零售价谁先动、符号解和数值解在什么条件下互相打架、利润函数在什么情况下不再是凹函数……这些问题不踩一遍光靠看论文很难有体感。这篇文章我想复盘一次完整的论文复现过程。目标很明确把“一个供应商 两个零售商”的经典二级供应链模型在Mathematica里从建模到求解全部跑通包括Stackelberg博弈的逆向归纳、一阶条件的符号推导、均衡解的数值验证以及那些论文里通常不会写、但实际复现时一定会遇到的坑。适合正在做供应链方向研究、准备用Mathematica做实验、或者想快速验证某篇论文结果的同学参考。1. 模型结构与博弈逻辑先把“谁决策”想清楚1.1 一个供应商两个零售商的典型设定这类模型的供应链层级不复杂上游是唯一的供应商制造商用批发价 (w) 向两个零售商供货下游是两个零售商分别面对自己的终端市场决定零售价或订货量。经典文献里通常采用线性逆需求函数[ p_i a - q_i - d q_j, \quad i,j \in {1,2}, i \neq j ]其中 (a) 是市场规模参数(q_i) 是零售商 (i) 的销量(d) 表示两个零售商之间的产品替代程度(0 \leq d 1)越接近1说明产品同质化越高。这个函数的意思很直观零售商 (i) 的价格不仅取决于自己的销量还会被对手的销量压制因为两个零售商卖的是有一定替代性的产品消费者会互相分流。我和很多刚开始复现这类模型的人一样最容易在第一眼就忽略一个关键点这个模型里的决策顺序不是同时的而是分层的。供应商先定批发价两个零售商看到批发价之后再各自决定自己的零售价或订货量。这就是典型的Stackelberg博弈而且是一个领导者供应商面对两个跟随者零售商的博弈结构。决策时序拆开看是这样的供应商基于对零售商反应行为的预判选择批发价 (w)目标是最大化自己的利润 (\pi_s w(q_1 q_2))。两个零售商观察到 (w) 后分别选择 (q_i) 最大化自己的利润 (\pi_{r_i} (p_i - w)q_i)。因为两个零售商在同一层、同时决策所以零售商之间存在一个Nash博弈而供应商在更上层需要用“零售商会对我的批发价做出什么反应”来倒推自己的最优决策。这种“上层领导者 下层Nash博弈”的结构在博弈论里叫双层博弈在供应链管理里就是最常见的分散决策场景。复现论文的第一步不是写代码而是把这个决策顺序用数学语言固化下来否则后面求导、代入解析解的时候顺序一乱结果全乱。1.2 逆向归纳法为什么从零售商开始推Stackelberg博弈的标准解法是逆向归纳Backward Induction。思路听起来简单——先从最后一步往回推但实际推导时很多人容易犯一个错误把零售商利润对 (q_i) 求导之后直接代进供应商的利润函数却忘了这一步求出来的是“反应函数”不是一个确定数值。具体来说零售商 (i) 的利润是[ \pi_{r_i} (a - q_i - d q_j - w) q_i ]对 (q_i) 求一阶条件[ \frac{\partial \pi_{r_i}}{\partial q_i} a - 2q_i - d q_j - w 0 ]于是得到零售商 (i) 对 (q_j) 的反应函数[ q_i \frac{a - w - d q_j}{2} ]两个零售商联立这个方程组解出来的才是均衡销量[ q_1^* q_2^* \frac{a - w}{2 d} ]注意这里得到的是关于 (w) 的表达式不是具体数字。它们才是供应商做决策时看到的“零售商的反应”。把 (q_1^, q_2^) 代进供应商利润[ \pi_s w(q_1^* q_2^*) \frac{2w(a - w)}{2 d} ]对 (w) 求一阶条件[ \frac{d\pi_s}{dw} \frac{2(a - w) - 2w}{2 d} 0 ]得到[ w^* \frac{a}{2} ]这个结果非常干净最优批发价只取决于市场规模 (a)与替代程度 (d) 完全无关。我第一次跑出来的时候也有点意外仔细想想又合理——供应商面对的是两个对称零售商的总需求而总需求 (q_1 q_2 2(a - w)/(2d)) 对 (w) 的敏感度并不受 (d) 影响所以供应商不用关心下游竞争有多激烈只管定一半的价格水平。把 (w^* a/2) 带回零售商反应函数可以得到均衡零售价和利润这一步我会在后面的Mathematica实现里完整展开。2. Mathematica 求解思路与代码实现从符号推导到数值验证2.1 为什么选Mathematica而不是手动推导这类供应链模型的代数推导其实不算难但有一个非常现实的问题当你给模型加入非对称参数比如两个零售商的成本不同、或者把需求函数改成非线性形式时手推的出错率会大幅上升而且一旦推到一半发现前面某一步符号错了回头看非常痛苦。Mathematica的核心优势在于符号计算Symbolic Computation。你只需要把利润函数写对求导、联立方程、代入消元这些机械但易错的操作全部交给内核。我复现这篇论文时几乎全部流程都是用Solve、D、Reduce这几个基础函数完成的不需要任何高级包代码量也压得很低。这里多说一句Mathematica的符号计算和MATLAB的符号工具箱不是一个量级的存在。MATLAB更偏数值符号运算只是附加功能Mathematica的内核设计就是围绕符号表达式展开的处理这种“先求反应函数、再联立求均衡、再代入求利润”的嵌套推导Mathematica的表达方式几乎和数学表达式一一对应调试心智负担小很多。2.2 核心代码完整的求解流程下面这段代码是我整理后的完整复现脚本包含了从模型设定到均衡验证的全部流程。我用的是Wolfram Language自然语法在Mathematica 12及以上版本中可以直接运行。(* 清空环境避免变量残留干扰 *) ClearAll[Global*] (* 1. 模型设定 *) (* 逆需求函数p_i a - q_i - d*q_j *) p[i_, j_] : a - q[i] - d*q[j] (* 零售商利润销售收入 - 进货成本 *) piR[i_, j_, w_] : (p[i, j] - w)*q[i] (* 供应商利润批发价格 * 总销量 *) piS[w_] : w*(q[1] q[2]) (* 2. 第二步零售商博弈Nash *) (* 对每个零售商的利润求一阶条件 *) foc1 D[piR[1, 2, w], q[1]]; foc2 D[piR[2, 1, w], q[2]]; (* 联立求解两个零售商的最优反应 *) retailerSolution Solve[{foc1 0, foc2 0}, {q[1], q[2]}] // Simplify; (* 提取均衡销量表达式 *) q1Star[w_] q[1] /. retailerSolution[[1]]; q2Star[w_] q[2] /. retailerSolution[[1]]; (* 3. 第一步供应商决策Stackelberg 领导者 *) (* 供应商利润 w * (q1* q2*) *) supplierProfit[w_] piS[w] /. {q[1] - q1Star[w], q[2] - q2Star[w]} // Simplify; (* 对 w 求一阶条件 *) focS D[supplierProfit[w], w]; (* 求最优批发价 *) wStar Solve[focS 0, w][[1, 1]] // Simplify; (* 4. 回代得到全部均衡值 *) q1Eq q1Star[w] /. wStar // Simplify; q2Eq q2Star[w] /. wStar // Simplify; p1Eq p[1, 2] /. {q[1] - q1Eq, q[2] - q2Eq} // Simplify; p2Eq p[2, 1] /. {q[1] - q1Eq, q[2] - q2Eq} // Simplify; piREq piR[1, 2, w] /. {q[1] - q1Eq, q[2] - q2Eq, w - (w /. wStar)} // Simplify; piSEq supplierProfit[w] /. wStar // Simplify; (* 输出关键结果 *) Print[最优批发价 w* , w /. wStar]; Print[零售商1均衡销量 q1* , q1Eq]; Print[零售商2均衡销量 q2* , q2Eq]; Print[零售商1均衡价格 p1* , p1Eq]; Print[零售商2均衡价格 p2* , p2Eq]; Print[零售商1均衡利润 , piREq]; Print[零售商2均衡利润 , piREq]; Print[供应商均衡利润 , piSEq]; (* 5. 数值验证 *) (* 设定参数并验证一阶条件与二阶条件 *) params {a - 100, d - 0.6}; numericW w /. wStar /. params; (* 验证零售商一阶条件在该批发价下是否为0 *) q1Num q1Eq /. params; q2Num q2Eq /. params; foc1Num (D[(a - q[1] - d*q[2] - w)*q[1], q[1]] /. {q[1] - q1Num, q[2] - q2Num, w - numericW}) /. params; Print[零售商1一阶条件数值 , foc1Num]; (* 验证供应商利润的二阶条件应小于0 *) socS D[supplierProfit[w], {w, 2}] /. wStar /. params; Print[供应商二阶条件数值 , socS];这段代码跑完之后你会得到一组非常漂亮的符号解[ w^* \frac{a}{2}, \quad q_1^* q_2^* \frac{a}{2(2d)}, \quad p_1^* p_2^* \frac{a(3d)}{2(2d)} ]再把均衡值代回利润函数[ \pi_{r_1}^* \pi_{r_2}^* \frac{a^2}{4(2d)^2}, \quad \pi_s^* \frac{a^2}{2(2d)} ]我建议你跑完代码后先别急着改参数做实验而是对着论文里的结果逐项核对这几个表达式的分母和系数。我复现的时候发现这类论文最常出现排版错误的地方就是利润表达式——分母多了个平方、或者系数从1/4写成1/2符号计算的好处就是你可以完全不依赖论文的最终公式自己推一遍就有绝对把握。2.3 代码里容易被忽略的三个细节这段代码虽然短但我实际调试时踩过不少坑分享三个最典型的细节一Solve返回的解列表顺序可能不固定。我在代码里用retailerSolution[[1]]提取第一个解但如果你换一个版本的Mathematica或者改了方程形式解的排列顺序有可能变化。稳妥的做法是加一个条件筛选比如选择解中 (q[1] 0) 的那个。处理办法可以这样写positiveSolution Select[retailerSolution, (q[1] 0 q[2] 0) /. # ];细节二变量代入顺序错位导致结果恒等于0。在计算零售商利润时如果你先把 (q_1^, q_2^) 代入利润函数再代入 (w^*)顺序反了会出现“利润表达式里还有 (w)但 (w) 已经没地方可代”的尴尬。我的建议是统一用/ .规则一次性替换不要让中间步骤里的变量分叉。细节三Simplify的时机决定结果的可读性。如果每次求导后不Simplify直接代入下一步表达式会膨胀到完全没法看。我在每个关键步骤后面都加了对齐的Simplify代价是多花一点点计算时间但结果非常干净。对于这种规模的问题计算开销完全可以忽略。3. 数值实验与结果对比验证模型行为的直觉3.1 固定参数的基准算例模型推导完毕接下来用数值实验验证均衡结果的合理性。设定市场规模 (a 100)替代系数 (d 0.6)代入公式得到变量符号表达式数值结果a100, d0.6最优批发价 (w^*)(\frac{a}{2})50零售商均衡销量 (q^*)(\frac{a}{2(2d)})19.23零售商均衡价格 (p^*)(\frac{a(3d)}{2(2d)})69.23零售商单方利润 (\pi_r^*)(\frac{a^2}{4(2d)^2})369.82供应商利润 (\pi_s^*)(\frac{a^2}{2(2d)})1923.08供应链总利润(\frac{a^2(3d)}{4(2d)^2} \frac{a^2}{2(2d)})2662.72这些数字本身并不神秘但有一个点很值得注意在这个基准算例里供应商的利润是单个零售商的5倍以上。这个差距不是偶然的。供应商作为“先动者”可以把批发价定在成本之上很高的一档而两个零售商因为互相竞争价格和销量都被压住了利润自然薄。这正是Stackelberg领导者的“先手优势”。3.2 替代系数对均衡结果的影响比基准算例更有意思的是参数敏感性分析。固定 (a 100)让 (d) 从0变化到0.9步长0.1观察各均衡值的变化趋势替代系数 d批发价 w*零售价 p*零售商销量 q*零售商利润 πr*供应商利润 πs*0.050.0075.0025.00625.002500.000.250.0071.4322.73516.532272.730.450.0068.7520.83434.032083.330.650.0066.6719.23369.821923.080.850.0065.0017.86318.881785.710.950.0064.2917.24297.261724.14这张表信息量很大。首先确认了之前的结论随着 (d) 增大批发价始终稳定在50完全不受下游竞争强度影响。零售价格持续下降从完全差异化时的75降到高度同质化时的64.29。零售商的利润几乎被腰斩从625降到297。更有意思的是供应链总利润的变化总利润从3125逐步降到2318。这意味着在分散决策下下游竞争加剧不仅伤害零售商也伤害整个供应链。原因在于双重边际效应——供应商定价时只看自己的利润不考虑这个批发价会如何扭曲零售商的定价决策。当 (d) 增大时零售商对批发价的反应变得更激烈这种扭曲被放大总利润损失就更大。这个结果给管理实践提供了很直接的启示如果上游供应商能影响零售商的竞争强度比如通过品牌区隔、产品线设计来降低替代性即使批发价不变整个供应链的利润空间也会被撑大。3.3 用Mathematica做敏感性分析的三种方式做这类参数扫描我比较推荐三种方式按适用场景排序方式一直接生成表格数据。适合快速、一次性看全貌。Table[{d, w /. wStar, (q1Eq /. {a - 100}), (p1Eq /. {a - 100}), (piREq /. {a - 100}), (piSEq /. {a - 100})} /. {d - dval}, {dval, 0, 0.9, 0.1}] // TableForm方式二用Manipulate做动态演示。适合在汇报或答辩时展示参数变化的连续趋势。我一般用Manipulate[ Plot[{p1Eq, p2Eq}, {d, 0, 1}, PlotRange - {50, 80}, AxesLabel - {d, 价格}, PlotLegends - {零售商1价格, 零售商2价格}], {a, 50, 200, 10}, {w, 20, 80, 5}]方式三求解析导数。在写论文时你需要给出“随着 (d) 增大(p^*) 单调递减”之类的严格结论那就直接求导数D[p1Eq, d] // Simplify跑出来是 (-a/(2(2d)^2))显然小于0说明零售价格确实随 (d) 递增而递减。4. 复现论文时的高频坑与排查思路4.1 符号解推出来了但数值验证对不上这是我复现这类论文时遇到最多的一个问题。符号求解非常顺利公式也漂亮但一旦代入具体参数去算和论文给出的数值结果有个别位数对不上。排查思路其实很清晰先检查是不是零售商的二阶条件不满足。这听起来很基础但在非对称模型里非常容易翻车。比如某个零售商的利润函数对 (q_i) 求二阶导[ \frac{\partial^2 \pi_{r_i}}{\partial q_i^2} -2 0 ]这个负号恒成立所以Nash均衡确实存在且唯一。但如果是非线性需求函数或者零售商之间存在成本不对称二阶条件就可能在某些参数区域失效。我的经验是每次跑完均衡解先顺手把一阶和二阶条件的数值验算写在代码里不要等到结果对不上了再回头查。另一个隐蔽的原因是参数范围越界。有些论文做了归一化处理比如把市场规模 (a) 设为1但你可能在实验里用了更大规模的参数比如 (a100)虽然模型是线性的很多比例关系不会变但如果你把两个不同量纲的参数混在一起比较就会出现看似矛盾的结论。建议复现时先把论文里的算例原封不动地跑一遍确认数值完全一致后再扩展。4.2 逆向归纳的顺序最容易错这个问题值得再说一遍Stackelberg博弈的“先定”和“后定”不代表求解顺序也是先供应商后零售商。很多初学者会先对供应商利润求导把 (w) 解出来再代到零售商的利润里——这个顺序完全反了。正确的做法是永远从产业链最下游开始推先把零售商的反应函数 (q_i(w)) 解出来再把它们代入供应商的利润函数。Mathematica里面这个顺序是通过“先算零售商问题后把结果作为规则代入供应商问题”实现的代码层面并不复杂但如果你在推导思路里没有把顺序理清就会在某个角落出现“(w) 被消掉了供应商利润变成常数”这类诡异结果。我自己的习惯是在代码开头用注释把决策时序写清楚像这样(* 博弈时序 1. 供应商选择 w 2. 零售商同时选择 q[i] 3. 利润实现 求解方式逆向归纳先解步骤2再解步骤1 *)这句话看起来多余但当你改模型加参数、或者在Mathematica里反复调试的时候它能帮你在思路混乱时迅速找回坐标。4.3 读论文公式和代码输出相差一个负号符号问题也不少见。我遇到过一次论文公式里写的均衡价格比代码输出的结果多了个负号差的怪事查了半天发现是论文里的逆需求函数定义和正文描述不一致——公式里写的是 (p a q_i d q_j)正文却说“价格随销量上升而下降”显然是排版或笔误。这种问题只有一种可靠解法不轻信论文的公式把模型设定里的函数形式和符号逻辑核对一遍。具体到这个模型里(p a - q_i - d q_j) 才能保证价格随销量下降如果是加号需求定律直接违背后面所有结果都失去经济意义。建议复现前先花5分钟从头验一遍模型的基本定性性质比出问题时瞎猜高效得多。4.4 Fast Solve 遇到多解怎么办当模型稍微扩展比如两个零售商的成本参数不同联立方程就会出现两组甚至多组解。论文里通常会直接跳过这一细节用Word排版写“解得”但你在复现时会发现Solve返回了两个解手工选错组就会得到负价格、负利润之类的荒谬结果。我的筛选策略是用物理意义约束销量、价格、利润都必须是正的。可以这样写validSolution Select[solutions, (q1 0 q2 0 p1 0 p2 0) /. # ]对于经济模型来说负均衡价格或负产量在大多数情况下意味着该解落在定义域之外可能是因为替代系数超过了某个临界值。只要把正约束加上通常能干净地去掉伪解。5. 扩展思考从复现到创新复现一篇论文的终点不应该只是让代码跑通、结果能对上。真正有价值的复现是你在跑通之后开始琢磨这个模型还能怎么改以“一个供应商 两个零售商”为起点几条自然的扩展方向值得试一试扩展一零售商之间由数量竞争改为价格竞争Bertrand竞争。把零售商的决策变量从 (q_i) 换成 (p_i)需求函数改为 (q_i a - p_i d p_j)整个推导过程的结构类似但均衡结果会不同。价格竞争通常比数量竞争更激烈零售商的利润会被进一步压缩供应商的批发价决策会不会因此受影响这个问题本身就够写一篇小论文。扩展二引入两个零售商的不对称性。比如一个零售商有线下渠道优势、另一个有线上流量优势体现在模型里就是两个零售商的需求函数截距不同(a_1 \neq a_2)或者在竞争项前的系数不同。这时供应商面对的不再是对称反应函数(w) 的最优解会发生什么样的偏移我试过把 (a_1100, a_280) 代入模型结果批发价会往强势零售商的方向倾斜但倾斜幅度并非直觉上的线性关系。扩展三加入供应商的公平关切或者风险规避。这意味着目标函数不再是简单的利润最大化而是效用最大化。比如用均值-方差效用 (U_s \pi_s - \lambda \mathrm{Var}(\pi_s)) 替换原来的利润函数这种改动在Mathematica里只需要把目标函数改一行但推导复杂度会明显上升——刚好可以检验Mathematica符号计算在更大规模问题上是否依然顺手。我在复现完基准模型后尝试了上述第一种扩展Bertrand竞争额外花了两天时间收获非常大。你会发现价格竞争下的零售商最优反应函数变成[ p_i \frac{a w d p_j}{2} ]注意这里的 (w) 出现在分子的加号位置与数量竞争下 (q_i (a - w - d q_j)/2) 的减号位置形成对比。表面上只是符号差异但代入后得到的批发价不再是 (a/2)而是和 (d) 相关的更复杂表达式。为什么价格竞争和数量竞争会产生如此本质的区别这个“为什么”比结论本身更值得写进论文里。如果你准备在这个方向深耕我建议拿到论文的第一时间先不看结论而是自己在Mathematica里把模型推一遍再对照论文结果找出分歧点。这样得到的理解深度远远超过“跑通代码、截个图”式的复现。根据我个人经验Mathematica处理供应链模型最舒服的地方在于你不需要在“推导公式”和“数值实验”之间切换两套工具。符号推完直接赋值就能跑数值画图、敏感性分析、动态演示也都在同一个环境里完成。把论文复现当成一次“思维演练”而不是“交作业”收获会大得多。最后再分享一个我常用的技巧在跑敏感分析前先用Simplify把均衡表达式的结构看清比如判断 (q_i^*) 表达式中 (a) 和 (d) 是否可分离这通常会帮你预判图表趋势省去不少可视化调参的时间。