Matlab求解Stackelberg博弈：双层优化与嵌套fmincon实现指南

简介：这套代码基于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 = 0
q2*(q1) = (a - c2 - b·q1) / (2b)

领导者把q2*(q1)代入π1，再对q1求导：

q1_Stack = (a + c2 - 2c1) / (2b)

代入a=100, b=2, c1=5, c2=8：Nash产量是q1=16.33、q2=14.83，而Stackelberg产量是q1=24.50、q2=10.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不需要额外依赖，对中低维连续问题足够可靠，风险是外层目标函数非凸时结果会很依赖初值。启发式搜索（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*(q1+q2) 的截距 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; % 价格非负，隐含 q1+q2 <= 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，这就是前面讲的嵌套结构。运行结果应当是q1=24.5、q2=10.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*(q1+q2)) - 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)) / ((N+1)*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*(q1+q2)) - 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*(Q1+Q2)) - 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代码真正用起来。

本文还有配套的精品资源，点击获取