★ FEATURED ARTICLE 简介这套代码基于Stackelberg博弈模型使用Matlab实现了对应的求解流程适用于计算机、电子信息工程、数学等专业学生的课程设计、期末大作业与毕业设计。代码采用参数化编程参数可灵活修改注释详细便于理解与二次开发压缩包内附案例数据解压后即可在Matlab 2014/2019a/2024a等版本中直接运行。资源包共包含24个文件主体为23个.m脚本另有1个.md说明文档整体大小仅8KB结构紧凑。从文件构成来看包含了多个功能函数如func_cA、func_u_2、P_ode等和主程序入口script.m覆盖了Stackelberg博弈中常见的状态方程、目标函数求解与结果输出等环节能够帮助用户快速掌握模型搭建与数值仿真方法。目前已有71人学习下载适合需要完成相关课题设计或博弈论仿真作业的初学者参考。通过阅读md说明文档及对照脚本中的注释可以了解每个函数的作用并根据自身需求调整参数快速获得仿真结果。1. Stackelberg-Game是什么一个主从决策问题为什么值得用Matlab写Stackelberg-Game斯塔克伯格博弈不是电子游戏而是博弈论里最经典的主从决策模型领导者先出牌追随者看清楚之后再做出自己的最优应对。定价、电网需求响应、供应链契约、无线资源分配现实中大量工程问题都带这种先后顺序。这类问题落到代码上是一个双层优化而Matlab的fmincon与优化工具箱正好能把整套逻辑快速搭起来。下面沿着一个产量竞争案例从逆向归纳法讲到嵌套fmincon的写法、参数调整和排查经验让你拿到这类Matlab代码包后能照着最小工程骨架改出自己的场景而不是把时间耗在查错上。2. 把Stackelberg博弈翻译成可求解的数学结构领导者先手、追随者响应的均衡逻辑2.1 先想清楚Stackelberg均衡和Nash均衡的差别谁先动谁说了算拿到代码先别急着跑第一步是搞清楚模型里的“决策顺序”是怎么被体现的。Nash均衡的世界里两个玩家同时出招每个人都只能基于对对方的预期做最优反应而Stackelberg博弈里领导者先动并且这个决策是公开、可观察、不可反悔的。追随者不需要猜领导者的策略它把领导者的产量当已知数解一个普通的单变量最大化问题。领导者则多算一步它知道追随者会怎么反应于是把追随者的反应函数嵌进自己的目标函数里再最大化。用产量竞争模型看最直观。设逆需求函数P a - b(q1 q2)领导者利润π1 q1·P - c1·q1追随者利润π2 q2·P - c2·q2。先看同时决策的Nash均衡。π1对q1求导、π2对q2求导令导数为0得到两个反应函数的交点q1_Nash (a - 2c1 c2) / (3b)q2_Nash (a - 2c2 c1) / (3b)再看Stackelberg。追随者先做最优反应dπ2/dq2 a - b·q1 - 2b·q2 - c2 0q2*(q1) (a - c2 - b·q1) / (2b)领导者把q2*(q1)代入π1再对q1求导q1_Stack (a c2 - 2c1) / (2b)代入a100, b2, c15, c28Nash产量是q116.33、q214.83而Stackelberg产量是q124.50、q210.75总产量从31.17升到35.25市场价格随之下降。先动不是简单地抢先一步而是通过公开承诺改变了对手的决策环境。领导者比Nash时多产了约50%这就是Stackelberg先动优势。判断一份代码实现的到底是哪一套博弈最简单的方法就看追随者的决策是不是在领导者的每次候选解之下实时求解。如果代码是两个优化器交替迭代到收敛那得到的是Nash而不是Stackelberg。这一条是看代码和写代码时的第一道分水岭。2.2 逆向归纳法把双层优化拆成下层最优反应上层套壳逆向归纳法Backward Induction从博弈树的最后一个决策节点开始倒着求解。Stackelberg里最后一个决策者是追随者所以先解追随者的最大化问题得到最优反应函数再把这个函数代入领导者的目标函数。这个逻辑映射到代码上就是嵌套结构内层求解器计算q2*(q1)外层求解器在这个函数之上再优化q1。这里有个容易被忽略的工程代价外层fmincon在数值求梯度时会对q1做微小扰动每扰动一次就触发一次完整的内层优化所以总计算量是外层迭代次数乘以内层迭代次数的量级。低维连续问题完全能接受但上层决策变量一旦多起来这个嵌套结构会慢到没法用。因此我一般会把能做解析的内层全部解析化。内层只有一个决策变量、目标函数可导时手推一阶条件只需要几分钟却能换来外层目标函数的光滑性和收敛速度。实际项目里我习惯先把解析反应函数写在注释里再用数值fmincon做交叉验证两边结果一致才说明代码没实现错。另外注意fmincon内部实际是在解KKT条件当追随者被边界约束卡住时内层返回的拉格朗日乘子会告诉你约束有多紧。如果你的内层结果经常压在下界0上说明追随者在当前参数下根本没有生产意愿这种边界情形要单独讨论不能当普通内点解处理。2.3 三条求解路线选型解析求导、fmincon嵌套、还是启发式在Matlab里做Stackelberg博弈求解常见的有三套做法选哪套取决于模型结构和探索阶段。求解路线适用场景实现成本精度主要风险解析求导函数可导、约束简单、结构稳定低高模型一变就要重新推导嵌套fmincon中低维连续变量双层优化中较高非凸时易陷局部最优启发式搜索套内层外层非凸、不光滑或离散变量高中计算量大、结果有随机性解析求导是最理想的验证基准它要求逆需求函数、成本函数都能解析处理一旦加约束或改非线性成本就得重新推所以只适合模型稳定后做理论分析不适合边探索边改的场景。嵌套fmincon是默认首选优化工具箱自带的fmincon不需要额外依赖对中低维连续问题足够可靠风险是外层目标函数非凸时结果会很依赖初值。启发式搜索ga或particleswarm套内层适合外层有多个峰甚至离散变量的情况但可复现性差必须固定随机种子并多次运行。我的选型标准很固定先数外层变量个数。小于5个、连续、目标函数大体光滑的直接用嵌套fmincon发现多峰就上多起点再不行才换ga。有人问我为什么不用Python做这套Python在数据预处理上确实方便但Matlab的优化工具箱在中小规模双层优化上的调试体验更顺手尤其fmincon的输出信息和可视化检查非常直接。3. 在Matlab里把Stackelberg博弈跑通嵌套fmincon的主程序与函数写法3.1 最小可运行骨架上层主程序下层最优反应函数下面这个脚本是完整的Stackelberg产量竞争实现内层用fmincon求追随者最优产量外层用fmincon求领导者最优产量。代码可以直接复制运行。% Stackelberg_main.m % 产量竞争场景领导者先定 q1追随者看到 q1 后再定 q2 clear; clc; % ------- 模型参数 ------- a 100; % 逆需求函数 P a - b*(q1q2) 的截距 b 2; % 价格对总产量的敏感系数 c1 5; % 领导者单位成本 c2 8; % 追随者单位成本 % ------- 优化器公共设置 ------- opts optimoptions(fmincon, ... Display, off, ... % 调试时改成 iter OptimalityTolerance, 1e-8, ... % 内层解要够准否则外层会抖动 ConstraintTolerance, 1e-8, ... MaxIterations, 1000); % ------- 内层追随者的最优反应函数 ------- % 输入 q1返回追随者的最优产量 q2 follower_obj (q2, q1) -( q2 * (a - b*(q1 q2)) - c2*q2 ); lb_q2 0; % 产量非负 ub_q2 a/b - q1; % 价格非负隐含 q1q2 a/b q2_opt (q1) fmincon((q2) follower_obj(q2, q1), ... 0.5*(a - c2 - b*q1)/b, ... % 初值取解析反应函数近似 [], [], [], [], lb_q2, ub_q2, [], opts); % ------- 外层领导者利润内层被包进来------- leader_obj (q1) -( q1 * (a - b*(q1 q2_opt(q1))) - c1*q1 ); lb_q1 0; ub_q1 a/b; q1_star fmincon(leader_obj, a/(2*b), ... [], [], [], [], lb_q1, ub_q1, [], opts); % ------- 结果输出 ------- q2_star q2_opt(q1_star); fprintf(领导者产量 q1 %.6f\n, q1_star); fprintf(追随者产量 q2 %.6f\n, q2_star); fprintf(市场价格 P %.4f\n, a - b*(q1_star q2_star));逻辑上要看清两个负号fmincon只能做最小化所以追随者利润和领导者利润都取负号再传入。外层每次评估leader_obj时都会完整运行一次内层的fmincon这就是前面讲的嵌套结构。运行结果应当是q124.5、q210.75与解析解一致。这里的参数说明值得细看。ub_q2写成a/b - q1把“价格必须非负”这个经济语义直接编成约束边界而不是靠目标函数里的负数去惩罚这样内层不会在负价格区间浪费迭代。内层初值我特意用了0.5*(a - c2 - b*q1)/b这是无约束解析反应函数的一半保证初值落在可行域内部且远离边界fmincon从接近最优的位置启动会快得多。OptimalityTolerance设到1e-8是为了让内层解足够干净否则外层看到的是一个带噪声的函数。3.2 参数怎么设初值、容差、迭代上限初值是嵌套优化里最影响成败的一项。外层初值我习惯放a/(2b)这个位置基本落在常见参数的吸引域内内层初值最好不要用固定常数而是随当前q1变化的解析近似这样q1在可行域内移动时内层fmincon始终从“接近最优”的地方起步。更激进的做法是用上一轮外层迭代的q2结果做初值但要注意记录顺序别把内外层变量搞混。容差设置要分内外层区别对待。外层只做决策输出OptimalityTolerance用1e-6就够内层是给外层提供函数值的必须更严格我一般设1e-8。如果把外层也设到1e-10嵌套计算量会迅速爆炸换来的精度提升在图形上都看不出来。还要注意StepTolerance它控制相邻两步决策变量的最小变化量嵌套收敛后期如果一直在这个阈值附近来回程序会判定无法继续下降而提前退出。迭代上限是最容易踩的坑。fmincon默认的MaxIterations是100乘变量数但在嵌套结构里外层一次迭代可能触发几十次内层优化默认值远远不够。我通常直接把外层MaxIterations设1000、MaxFunctionEvaluations设5000内层同理会因为反复被调用而需要更大的评估预算。不同Matlab版本的默认容差和迭代上限并不一致一个在旧版本上正常收敛的模型换到新版本可能因为默认MaxIterations变化而中途退出所以这些参数永远要显式写出来不要依赖默认值。3.3 内层函数怎么组织匿名函数与独立function的边界嵌套结构的代码组织方式会影响调试效率。小规模场景里匿名函数够用但一旦进入多追随者或多次改模型的阶段匿名函数捕获外部变量的行为会带来隐蔽Bug。我一般会把内层优化封装成独立的function文件参数通过结构体透传。function q2 follower_fmincon(q1, p) % p 是结构体包含 a b c2 opts 等字段 obj (q2) -( q2 * (p.a - p.b*(q1q2)) - p.c2*q2 ); q2 fmincon(obj, 0.5*(p.a - p.c2 - p.b*q1)/p.b, ... [], [], [], [], 0, p.a/p.b - q1, [], p.opts); end这样主脚本只需要写成q2_star follower_fmincon(q1_star, p)调试时可以直接在这个函数文件里加断点观察每次外层调用时内层走了多少步、收敛到哪个点。用结构体p传参比传一堆散装参数干净得多也比全局变量安全。匿名函数因为会捕获工作区变量在循环里频繁重建时很难排查独立function没有这个问题每次调用只创建一次临时匿名函数性能也更稳定。这种结构还有一个额外好处当你需要把内层换成解析解或fsolve时只需要改这一个文件的内部实现外层代码完全不用动。Stackelberg代码的迭代核心在内外层接口是否干净而不是在某个具体求解器上。4. 从报错到结果不对Stackelberg代码的5个排查记录这一章是实打实的踩坑记录。以下5个问题我都在不同模型里遇到过按“现象→原因→解决”来写每一条都能直接对照自查。4.1 坑1内层目标不是凹的最优反应函数取到了局部极值现象代码逻辑和公式一致但算出来的均衡和解析解差一截换个内层初值结果又变了。原因fmincon只能保证找到局部最优。产量竞争里线性逆需求函数让利润函数对q2是凹的一般不会出事但改成非线性价格P a - b·Q²或者成本函数带凸项之后凹性就未必保得住。内层一旦多峰fmincon返回的只是一个局部最优反应领导者基于错误的反应函数做决策整个均衡都是错的。解决先画内层目标函数曲线把q2当作横轴、利润当作纵轴肉眼确认有没有多峰。有峰就内层也做多起点或者换成fmincon配合全局搜索。我还会记录每个起点的内层最优值如果多个起点给出的数值不同说明内层本身就不可信这时候的数值均衡必须废弃先回头改模型。4.2 坑2两层嵌套不收敛外层目标函数在抖动现象外层迭代中目标函数值忽大忽小收敛后打印q2_opt(q1_star)和上一次迭代传出来的q2对不上。原因内层容差太松返回的q2精度不够外层看到的不是一个确定性函数而是一个带噪声的函数。此时外层用有限差分求梯度微分步长小于内层解的误差算出来的梯度全是噪声迭代自然震荡。解决把内层OptimalityTolerance和ConstraintTolerance都收到1e-8外层有限差分步长用默认值。如果还抖就直接手推内层一阶条件用解析反应函数替换内层fmincon。内层解析化是嵌套优化里最有效的后悔药它把整个问题变成一个光滑的单层优化振荡问题会直接消失。4.3 坑3最大化写成最小化解直接飞去了边界现象解出来的产量不是0就是上界打印出来的“利润”是一个巨大的负数或者fmincon提示目标函数返回了异常值。原因Matlab的fmincon、fminunc全都是最小化器要求最大化利润必须取负号。初学者经常在层与层之间漏掉一个负号或者领导者取负了、追随者没取负内外层符号不一致。解决给每个目标函数都单独检查一遍利润收入-成本传入fmincon的目标负利润。我养成了一个调试习惯——在主脚本里手动给定一个q1单独调用q2_opt打印追随者利润确认它是正数再调用leader_obj(q1)确认领导者利润也是正数。数值量级和符号都对再交给优化器。4.4 坑4目标机器报“msvcp140.dll无法继续执行代码”现象代码在自己电脑上跑得好好的用Matlab Compiler打包成exe发给同事对方双击就弹窗提示“由于找不到msvcp140.dll无法继续执行代码”。原因编译后的exe依赖Matlab Runtime而Runtime依赖Microsoft Visual C 2015-2022 Redistributable。目标机器没装过Visual Studio或对应运行库就会报这个错。这跟博弈代码本身没关系纯粹是部署环境缺组件。解决在目标机器上装一次VC Redistributable如果使用人数多直接把运行库安装包和exe放同一目录写个安装脚本一并处理。注意Runtime版本要和打包用的Matlab版本匹配混装不同年份的Runtime容易引发奇怪行为。4.5 坑5嵌套循环里反复重建句柄慢到怀疑人生现象外层迭代次数一多程序越跑越慢用profiler看大量时间耗在匿名函数构造和options创建上。原因如果把q2_opt、leader_obj定义在循环内部每轮循环都重新创建匿名函数和options对象。Matlab是解释型平台重复构造的开销很可观而且匿名函数在循环里捕获循环变量还容易引入隐蔽Bug。解决把options和函数句柄提到循环外只建一次内层目标函数写成独立function文件见3.3参数用结构体传。性能优化放最后做先保证逻辑正确但“句柄只建一次”这个习惯建议从一开始就保持。5. 从单追随者到多追随者内层方程组与全局求解器升级5.1 多追随者内层从单点最优变成方程组求解单追随者模型的内部结构是一个最大化问题但换成多个追随者时内层变成一个小型Nash博弈——追随者之间是同时行动的它们看到q1之后各自定产量最终收敛到彼此反应函数的交点。这个交点不是一个极值点而是一个方程组的根用fsolve求解是最直接的做法。function q2_vec follower_response_multi(q1, params) % 多追随者各自最优产量的联立方程组求解 a params.a; b params.b; c2 params.c2; N length(c2); function F foc(q2) Q q1 sum(q2); F zeros(N, 1); for i 1:N % 追随者 i 的一阶条件a - b*Q - b*q2(i) - c2(i) 0 F(i) a - b*Q - b*q2(i) - c2(i); end end % 初值从对称平均产量附近起解避免全 0 导致的雅可比奇异 q0 max((a - max(c2)) / ((N1)*b), 0.1) * ones(N, 1); q2_vec fsolve(foc, q0, params.fsolve_opts); end这里的数学结构要拎清楚多追随者之间不再是Stackelberg关系而是Cournot关系它们同时行动各自的一阶条件联立成方程组。外层仍然由领导者主导但每个追随者的决策同时依赖于q1和其他追随者的产量。fsolve默认用信赖域算法N在50以内一般很快初值别给全0否则雅可比矩阵在边界处可能奇异。需要说明的是内层追随者之间存在纳什均衡的唯一性问题理论上要验证反应函数是压缩映射才保证收敛实际操作里可以多换几组初值确认解的稳定性。5.2 外层求解器升级多起点与ga的取舍外层目标函数如果只有单一峰直接fmincon就够了一旦出现多峰就需要多起点搜索。多起点的思路很简单随机生成一批初始q1各自用fmincon收敛最后取负利润最小即真实利润最大的那个解。% 多起点外层搜索对领导者目标函数做 20 次随机起始点 N_starts 20; q1_best NaN; neg_pi1_best inf; rng(0); % 固定随机种子保证结果可复现 for k 1:N_starts q1_0 rand * a / b; % 随机初始产量 [q1_try, fval] fmincon(leader_obj, q1_0, ... [], [], [], [], lb_q1, ub_q1, [], opts); if fval neg_pi1_best neg_pi1_best fval; q1_best q1_try; end end如果想全局性更强可以直接换ga。ga不需要梯度对外层不光滑的函数更稳健但每次运行有随机性必须固定rng且加大种群规模来摊平方差。实际使用时要看清工具箱版本ga属于Global Optimization Toolbox不是所有Matlab安装都带。我的习惯是先多起点跑完把最优解作为初值再喂一次fmincon做局部精修这样既保住全局性又有梯度法的精度。不管用哪种方法外层目标函数对领导者来说仍然是一个黑匣子你只能通过输入输出观察它的行为所以每次改完参数后固定随机种子是必须养成的纪律。5.3 敏感性分析扫追随者成本看均衡往哪儿移模型调通之后下一步通常是做敏感性分析看均衡如何随参数漂移。比如扫追随者成本系数c2观察q1_star和q2_star的变化路径。这类扫描能快速验证模型是否有经济直觉。% 扫描追随者成本系数 c2观察均衡产量移动 c2_range linspace(4, 12, 30); q1_hist zeros(size(c2_range)); q2_hist zeros(size(c2_range)); for k 1:length(c2_range) c2 c2_range(k); q2_opt_c (q1) fmincon((q2) -(q2*(a - b*(q1q2)) - c2*q2), ... (a - c2 - b*q1)/(2*b), [], [], [], [], 0, a/b - q1, [], opts); leader_obj_c (q1) -(q1*(a - b*(q1 q2_opt_c(q1))) - c1*q1); q1_hist(k) fmincon(leader_obj_c, a/(2*b), ... [], [], [], [], 0, a/b, [], opts); q2_hist(k) q2_opt_c(q1_hist(k)); end figure; yyaxis left; plot(c2_range, q1_hist, -o); ylabel(q1); yyaxis right; plot(c2_range, q2_hist, --s); ylabel(q2); xlabel(追随者成本 c2);扫描结果通常是q1随c2上升而上升、q2下降——追随者成本越高越不想生产领导者趁机扩产抢占市场。如果你的模型扫出反直觉变化先别急着接受回到第2章的博弈逻辑里找原因是内层出现了多解还是外层陷入了局部最优。6. 验证解是否站得住等高线地形图与解析解交叉检验6.1 把领导者目标函数地形和最优反应曲线叠在一张图上数值解跑完后第一个验证动作是画地形图。把q1和q2扫成网格计算领导者利润并画等高线再把追随者的最优反应曲线叠上去均衡点应当落在等高线与反应曲线的切点附近。q1g linspace(0.1, 60, 300); q2g linspace(0.1, 30, 300); [Q1, Q2] meshgrid(q1g, q2g); PI1 Q1 .* (a - b*(Q1Q2)) - c1*Q1; contour(Q1, Q2, PI1, 40); hold on; fplot((q1) max(0, (a - c2 - b*q1)/(2*b)), [0 60], r, LineWidth, 1.5); plot(q1_star, q2_star, ko, MarkerFaceColor, k); xlabel(q1); ylabel(q2);如果黑点明显偏离切点要么内层没收敛要么外层掉进了局部最优这张图能一眼看出问题。6.2 解析解对比误差超过1e-3就回第4章排查线性逆需求下Stackelberg均衡有闭式解直接比对数值解和解析解是最快的检验手段。q1_ana (a c2 - 2*c1) / (2*b); q2_ana (a - c2 - b*q1_ana) / (2*b); assert(abs(q1_star - q1_ana) 1e-3, 数值解与解析解偏差过大);我把这套验证脚本存成固定模板每次改完模型先跑一遍解析对比再做敏感性扫描。这几年下来多数冤枉路都发生在省掉这一步的时候。算例里跑出来的误差超过阈值不要急着调容差先用第4章的清单逐条排除问题多半藏在符号、初值或内层收敛性里而不是精度设置上。希望这些经验能帮你少走一段弯路也希望你带着这套验证习惯把Stackelberg代码真正用起来。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
