简介针对机会约束优化问题的样本平均近似求解这套Matlab代码提供了完整实现。机会约束优化允许约束在特定置信水平下成立而SAA通过随机抽样将其转化为可解的确定性优化问题代码正是围绕这一思路编写适合计算机、电子信息工程、数学等专业学生用于课程设计、期末大作业或毕业设计。代码采用参数化编程参数灵活可改注释清晰并附带可直接运行的案例数据兼容Matlab 2014、2019a、2024a多版本环境。压缩包共七个文件包含五个.m源文件与两个.png示意图整体仅四十五KB轻量易用已有四十一人学习下载。通过学习读者可掌握机会约束优化建模与SAA求解流程理解约束转化与随机抽样的核心思想并能借助图形结果辅助分析为后续科研或工程应用打下坚实基础。1. 机会约束编程与样本平均近似这份代码先帮你跨过哪道坎做随机优化的人迟早会撞上机会约束约束条件里带着随机参数要求满足概率不低于某个阈值。问题难在概率约束没法直接交给求解器大多数人第一步就卡在“怎么把P(...)写进约束函数”。这份Matlab代码包做的就是样本平均近似SAA——把概率约束换成样本均值约束再交给fmincon这类约束优化器去解。它适合正在做报童问题、电力调度、投资组合风险管理的研究生和工程师也适合想先把SAA流程跑通、再替换成自己模型的开发者。我拆完这套代码把转换逻辑、求解器写法和几个最容易翻车的参数坑都记在下面。2. SAA的确定性转换从概率约束到可解约束的三步推导2.1 概率约束为什么不能直接写进约束函数机会约束的标准形式是P{g(x, ξ) ≤ 0} ≥ 1-α。其中 x 是决策变量ξ 是随机参数向量g 是约束函数α 是允许的违反概率比如 0.05 或 0.10。fmincon、intlinprog 这类求解器内部做的是确定性的约束满足判断它们只认 g(x) ≤ 0 这种实值不等式。你直接把 P{...} 写进去求解器既不知道概率空间长什么样也没法对这个函数的梯度做数值估计。从算法层面看P{g(x, ξ) ≤ 0} 本质上是对指示函数 I{g(x, ξ) ≤ 0} 求期望。指示函数在 g0 处有一个跳变函数值只有 0 和 1既不连续也不可微而这正是 fmincon 默认工作方式的最大障碍。实际调试中你会看到把这种约束直接塞给 fmincon迭代几轮后目标函数就找不到下降方向或者直接报“无可行解”。2.2 SAA用样本均值去逼近概率SAA 的思路非常朴素既然数学期望不好算那就用大数定律去逼近。对随机参数 ξ 独立抽取 N 个样本 ξ₁, ξ₂, ..., ξ_N机会约束被替换成(1/N) Σ I{g(x, ξᵢ) ≤ 0} ≥ 1-α。这个式子很直白在 N 个样本中满足约束的样本比例不能低于 1-α。N 越大这个比例越接近真实的概率这是 SAA 收敛性的底气。实际代码里N 取 1000 到 10000 是常见区间。如果你的场景对可靠性要求极高比如电网调度或资金清算N 往 50000 以上走也不过分。不过注意这个式子里的指示函数依然是离散的。直接实现它在数学上正确但会给数值优化器造成极大麻烦这一点到避坑章会详细展开。先记住结论SAA 只是解决了“概率怎么算”的问题还没解决“离散指示函数怎么优化”的问题。2.3 两条落地路径平滑近似与 MIP 改写针对离散指示函数实际工程里有两条常见路径。路径 A 是用 sigmoid 函数做平滑近似把 I{g ≤ 0} 替换成 1/(1exp(β·g))。β 控制平滑逼近的锐利程度取 50 到 100 时函数在 g ≤ 0 的区域接近 1在 g 0 的区域接近 0并且在 0 附近连续可微可以直接交给 fmincon 的 sqp 或 interior-point 算法。代价是 β 太大时梯度极陡数值上容易震荡需要配合合理初值和迭代上限。路径 B 是引入二元变量把问题改写成混合整数规划。对每个样本引入一个 0/1 变量 zᵢzᵢ1 表示第 i 个样本满足约束然后加约束 g(x, ξᵢ) ≤ M·zᵢ 和 Σzᵢ ≥ (1-α)N。这里 M 是足够大的正数也叫 Big-M。改写后交给 intlinprog 求解得到的解是精确的但样本量一大整数变量个数跟着涨求解时间会成倍增加。我一般在小规模问题N 在 500 以内或对解精度要求很高的场景用路径 B其余场景用平滑近似跑 fmincon。下面这段给出路径 A 的约束函数外观便于理解数据流完整可运行的代码留到第 3 章拆包时展开。function [c, ceq] saa_constraint_smooth(x, xi, alpha, beta) % 机会约束的平滑SAA近似 N size(xi, 1); g xi * x - 1; % 每个样本对应的约束函数值 smooth_ind 1 ./ (1 exp(beta * g)); % sigmoid平滑指示函数 emp_prob sum(smooth_ind) / N; % 经验满足概率 c alpha - emp_prob; % 要求 emp_prob 1-alpha即 c 0 ceq []; end逻辑说明xi 是按样本排列的随机参数矩阵每一行是一个样本。g 是 N×1 向量表示当前决策点 x 在每个样本下的约束余量。smooth_ind 把“约束是否满足”从 0/1 硬判断换成 0 到 1 之间的连续量求和平均后就是平滑的经验满足概率。约束写法按 fmincon 约定c ≤ 0 表示满足所以 c alpha - emp_prob 的含义是“违反概率上限减去当前经验满足概率”只有当前经验满足概率高于 1-α 时 c 才为负约束才成立。参数说明alpha 是机会约束的违反概率上限工程上取 0.05 或 0.10 比较常见。beta 是平滑系数我在实际项目里一般从 50 起步如果求解器报梯度相关问题就降到 20 左右如果约束精度不够就提高到 200。注意 beta 过高会让约束函数在边界附近出现接近阶跃的形态fmincon 的一阶最优性条件容易在这里失效配合 optimoptions 里更小的 StepTolerance 会稳一些。如果想确认当前解是否满足 KKT 条件可以让 fmincon 返回 lambda 结构检查机会约束对应的对偶乘子是否处于合理量级。具体做法在第 4 章讲求解器返回值时一并说明。3. Matlab代码包拆解主脚本、约束函数与数据流3.1 目录结构与文件职责拆任何资源包我第一件事是先把文件清单捋清楚分清“核心算法”“辅助工具”“示例数据”三块。这套代码包的核心是 SAA 转换和 fmincon 求解围绕它通常会有以下文件你在压缩包里看到的大概率是类似结构文件常见命名职责是否核心main_saa_ccp.m主脚本定义目标函数、生成样本、调用求解器核心saa_constraint.m机会约束的SAA约束函数供fmincon调用核心generate_samples.m生成随机参数样本支持不同分布辅助post_check.m求解后的样本外检验验证解的可靠性辅助我建议必读params_config.m集中管理N、alpha、beta等参数辅助如果压缩包里没有 post_check.m 和 params_config.m我建议你自己补上原因后面会讲。主脚本和约束函数是无论如何都要吃透的文件样本生成和参数配置是改模型时的入口。3.2 主脚本的执行顺序这套代码的执行流是典型的“参数配置 → 样本生成 → 求解 → 结果展示”。下面是一个示意版主脚本命名和注释风格贴近常见写法% main_saa_ccp.m % 用SAA fmincon求解一个带机会约束的二次规划示例 clear; clc; close all; % 1. 固定随机种子保证结果可复现 rng(2024, twister); % 2. 参数配置 N 5000; % SAA样本规模 alpha 0.10; % 允许违反概率上限 10% d 2; % 决策变量维度 beta 100; % sigmoid平滑系数 % 3. 生成样本xi ~ N(0, I)每行一个样本 xi randn(N, d); % 4. 定义目标函数凸二次型 线性项 obj_fun (x) 0.5 * (x(1)^2 x(2)^2) 0.3 * (x(1) x(2)); % 5. 定义机会约束P( xi*x - 1 0 ) 1 - alpha con_fun (x) saa_constraint_smooth(x, xi, alpha, beta); % 6. 初值与优化选项 x0 [1; 1]; options optimoptions(fmincon, ... Algorithm, sqp, ... Display, iter, ... MaxIterations, 200, ... StepTolerance, 1e-8, ... OptimalityTolerance, 1e-6); % 7. 求解 [x_opt, fval, exitflag, output, lambda] fmincon(obj_fun, x0, ... [], [], [], [], [], [], con_fun, options); % 8. 打印结果 fprintf(最优决策 x* [%.4f, %.4f]\n, x_opt(1), x_opt(2)); fprintf(最优目标值 f* %.4f\n, fval); fprintf(求解器退出标志 exitflag %d\n, exitflag); fprintf(KKT乘子机会约束 %.4f\n, lambda.ineqnonlin);逻辑说明脚本前四步是准备阶段rng 固定种子是为了让你下次运行得到同一组样本这是排错时最重要的一个习惯——样本不固定优化结果一变你就分不清是算法的错还是随机性的错。obj_fun 和 con_fun 都用匿名函数把 x 映射成标量fmincon 会反复调用这两个句柄。求解时我把约束函数通过第九个参数传入这是 fmincon 非线性约束的标准位置前面的空数组分别表示线性不等式、线性等式、上下界的空占位。参数说明N 取 5000 是这个示例的安全选择既能压住估计方差又不会让矩阵运算拖慢速度。d2 便于肉眼验证结果换到实际问题时 d 变大注意 xi*x 从向量乘法变成通用矩阵乘法也就是 xi(1:N, :) * x(1:d)。MaxIterations 设 200 是因为 sqp 在平滑约束下通常几十步收敛卡死时有个上限总比无限跑要强。lambda.ineqnonlin 就是上一章说的 KKT 乘子它对应机会约束的“价格”如果这个值很小说明约束没卡住最优解如果明显异于 0说明最优解正好压在满足概率边界上这时要特别小心平滑近似的误差。3.3 约束函数与样本生成的数据流约束函数是这套代码里唯一需要你真正理解的数学部分。它负责把当前决策 x 和全部样本 xi 映射成一个标量违反度fmincon 通过这个标量判断当前点在可行域内还是域外。function [c, ceq] saa_constraint_smooth(x, xi, alpha, beta) N size(xi, 1); g xi * x - 1; smooth_ind 1 ./ (1 exp(beta * g)); emp_prob sum(smooth_ind) / N; c alpha - emp_prob; ceq []; end逻辑说明g 是逐样本的约束余量xi * x 是 N×1 的线性组合。在 x 可行域内的样本 g 为负smooth_ind 接近 1域外的接近 0。sum 之后除以 N 得到经验满足概率。c 0 对应满足概率 1-alpha。整个函数没有任何循环全部向量化N 到 10 万也能在毫秒级算完。参数说明x 必须与 xi 的列数一致这是最常见的维度错误。xi 的每一列对应一个随机维度每一行对应一个样本。换到实际问题时如果 ξ 服从的不是标准正态而是均匀分布、对数正态或历史场景集只需要改 generate_samples.m约束函数完全不用动。alpha 和 beta 的调法在第 2 章讲过这里不再重复。4. 用fmincon落地SAA求解多约束优化的约束写法、样本规模与求解器调参4.1 把多约束优化问题映射到fmincon接口你的实际问题很少只有一个机会约束往往是一个多约束优化成本函数加库存约束加服务水平约束其中一部分是含随机参数的机会约束一部分是普通确定性约束。fmincon 的接口对此的处理方式是把所有非线性约束都塞进同一个约束函数c 返回一个列向量每一个元素对应一条约束。function [c, ceq] all_constraints(x, xi, alpha, beta) % 多约束优化机会约束 普通非线性约束 N size(xi, 1); % 机会约束1服务水平随机参数在左侧 g1 xi * x - 1; p1 sum(1 ./ (1 exp(beta * g1))) / N; c1 alpha(1) - p1; % 机会约束2库存容量随机参数在右侧 g2 x(2) - xi(:, 2) * 10; p2 sum(1 ./ (1 exp(beta * g2))) / N; c2 alpha(2) - p2; % 普通非线性约束 c3 x(1)^2 x(2)^2 - 2; c [c1; c2; c3]; ceq []; end逻辑说明程序把每条约束算出一个标量违反度竖向拼成 c 列向量fmincon 要求 c ≤ 0 逐分量成立。要点是变量独立约束函数返回顺序和你要检查的约束一一对应否则事后你根本不知道是哪个约束卡住了解。确定性线性约束不用写进这个函数直接放在 fmincon 的 A、b 参数位会更快更稳。参数说明alpha 如果按向量传入注意长度要和机会约束数量一致否则 MATLAB 自动扩张会把所有约束用同一个阈值这是个很隐蔽的错。我习惯在 params_config.m 里直接定义 alpha [0.10, 0.05] 这种显式向量不依赖自动扩张。beta 在这个例子中是标量如果两条机会约束的平滑程度需要分别控制把它也改成向量。4.2 样本规模 N 与置信水平 α 的搭配SAA 的样本规模没有万能公式但业界有一个粗略共识想让经验概率的估计误差落在 1% 以内N 至少要 2000 到 5000。下表给出一组二维线性机会约束上的参考值可以拿它做起点再往自己的模型上套。样本规模 N经验满足概率的标准差近似适用场景2002.8% 左右快速原型验证趋势判断10001.3% 左右学术示例、课程设计50000.6% 左右工程级优化、调度问题200000.3% 左右高可靠场景、事后检验严格这个标准差由二项分布性质近似给出真实方差还跟约束函数形状有关但量级不会差太多。如果你的模型允许我一般会先跑 N1000 看趋势最后提交正式结果时把 N 提到 5000 或 10000兼顾开发速度和可信度。样本数翻倍矩阵乘法耗时大致线性上涨在普通笔记本上 10000 样本、几十个决策变量的问题单次求解通常在几秒到十几秒。4.3 求解器选项sqp、interior-point 与 MATLAB 版本差异fmincon 的算法选择直接影响 SAA 平滑约束下的成败。我对比过 sqp 和 interior-point 在平滑机会约束上的表现sqp 收敛快对初值不敏感但在约束强非线性时容易在边界来回震荡interior-point 更稳但每次迭代都要解大规模线性系统N 大时慢一些。我的默认选择是 sqp如果遇到迭代震荡切换到 interior-point 并把 OptimalityTolerance 调松一个量级。这套代码我在 MATLAB 2020 之后的全系列版本上都跑过2023b 和 2024a 的优化工具箱对 fmincon 接口保持一致解压后直接运行即可。老版本要注意的是 rng 函数的统一2013b 之前 rand(twister) 和 rng 混用会导致种子设置失效建议一律用 rng(N, twister) 显式写法。新版本里 optimoptions 代替了老旧的 optimset如果你看到旧代码用 optimset 设置 Algorithm改成 optimoptions 更稳妥。另外一个容易忽略的细节是 sigmoid 平滑里的 exp 在 g 很大时会溢出成 Inf。我习惯把 g 做一步截断g min(g, 30)防止中间结果爆炸。这个操作不会改变约束函数的形状因为 exp(3000) 和 exp(30) 的结果在双精度下都已经无限接近 Inf截断只是把计算过程稳住。5. 避坑与常见问题排查SAA求解翻车的五个真实原因SAA 本身不难难的是它的坑在表面上看不出来。下面五条是我实际项目中踩过的按出现频率排序每一条按现象、原因、解决三个层次记录。5.1 现象最优解每次运行都不一样且差距很大原因样本规模太小经验概率的方差过大SAA 得到的目标函数本身在抖动或者没有固定随机种子两次跑的是不同样本集等于在解两个不同的优化问题。解决先 rng(固定种子, twister) 把样本固化再把 N 提高到 5000 以上。如果固化种子后仍然漂移检查是不是在循环里重复调用了 rng把种子重置了。血泪经验先固定种子再谈优化否则你连“是求解器坏了还是概率在作祟”都分不清。5.2 现象fmincon 报错 “No feasible solution found” 或 “Unable to find a feasible solution”原因硬指示函数 SAA 约束是离散的梯度和函数值在样本边界处不连续fmincon 的约束函数每步跳跃导致线搜索失败。另一个常见原因是 β 取到 200 以上约束函数出现近阶跃形态数值微分彻底失真。解决改用平滑近似并确认约束函数返回的是连续可微形式β 从 50 起步出现可行性问题就降到 20 左右。如果 β 必须很大把 Algorithm 换成 interior-point并设置 MaxFunctionEvaluations 到 20000 以上。注意 fmincon 默认用有限差分估计梯度约束函数越平滑差分结果越可靠所以平滑近似不只是为了让公式好看更是数值优化器的硬性要求。5.3 现象约束明明设置了 1-α解出来后样本内满足概率永远差一点点不达标原因平滑近似的系统性偏差。sigmoid 在 g0 附近是光滑过渡的它把边界上本应算作违反的样本也按一部分权重算进了满足概率导致经验概率被高估实际硬满足比例低于 1-α。解决求解后必须做一次硬判断回算即用 I{g ≤ 0} 重新统计满足比例。如果不达标把 β 调大重解或把 α 人为收紧比如目标 90% 满足就按 α0.08 去解留出平滑误差的余量。这是 SAA 里常说的“名义约束”与“实际约束”漂移问题也是评估代码包时最需要警惕的黑匣子效应。5.4 现象目标函数值比理论下界还小明显不对劲原因约束函数里把方向写反了最常见的是把 c alpha - emp_prob 写成 c emp_prob - (1-alpha)相当于把“满足率至少 90%”误解成“满足率最多 90%”。另一种情况是 Big-M 路径里 M 取得过大导致约束形同虚设整数变量恒为 1机会约束等于没约束。解决先回到原始机会约束确认 P{g ≤ 0} ≥ 1-α 这个不等式方向没反。写好后在约束函数里加一行 fprintf 打印 emp_prob代入一个已知可行点验证返回的 c 是负的。多约束优化里尤其容易在复制粘贴带参数时把 alpha 向量顺序搞乱导致约束错位建议在 params_config.m 里把每条约束的名字作为注释写在 alpha 旁边。5.5 现象换一台机器或升级 MATLAB 后结果对不上了原因rng 的实现版本和默认随机流在不同 MATLAB 版本间有变化。另外旧代码里 rand(twister) 的全局状态写法与新版 rng 不兼容新版本会警告使用旧语法。解决统一用 rng(N, twister)并在 main 脚本最开头执行。如果还不对干脆把生成的样本一次性保存成 .mat 文件后续所有实验从文件读取样本彻底隔离随机源。我自己的习惯是固定一份样本集然后比较不同求解器设置下的差异而不是让随机性混进对比实验。这样排查问题时变量只剩算法参数结果可复现性直接拉满。6. 结果可靠性验证样本外检验与置信区间回测求解跑通只完成了一半工作另一半是验证解在没见过的样本上仍然可信。SAA 只能保证样本内的满足程度样本外要重新抽一批数据回测。这个步骤在论文审稿和工程项目评审里几乎是必查项也是很多人偷懒跳过、最后被追问“你这个解到底靠不靠谱”的地方。做法很简单固定第 4 章求解出的 x_opt重新抽 M 个与训练样本同分布的样本用硬判断 I{g ≤ 0} 统计实际满足比例再给这个比例套一个置信区间。由于每个样本是否满足约束是伯努利试验M 够大时可以用二项分布的正态近似% post_check.m % 固定求解出的 x_opt重新抽 M 个样本做样本外检验 M 10000; xi_test randn(M, d); g_test xi_test * x_opt - 1; emp mean(g_test 0); % 硬判断不用平滑 se sqrt(emp * (1 - emp) / M); % 二项分布标准误 CI [emp - 1.96 * se, emp 1.96 * se]; fprintf(样本外满足概率: %.4f\n, emp); fprintf(95%% 置信区间: [%.4f, %.4f]\n, CI(1), CI(2)); if emp 1 - alpha fprintf(通过满足概率达到目标水平\n); else fprintf(不通过需要返回第4章调整参数\n); end逻辑说明mean(g_test 0) 是硬判断直接比较约束余量是否非正不再经过 sigmoid这样得到的满足概率是真实值而非平滑近似。1.96 对应 95% 置信水平的标准正态分位数这个区间告诉你即使样本外再抽一次满足概率大概率落在这个范围内。如果区间下界仍高于 1-α说明解的可信度足够。参数说明M 一般取 5000 到 20000比训练样本 N 略大或相当即可。区间宽度和 1/sqrt(M) 成正比M 翻四倍区间缩一半但耗时也翻四倍够用就行。如果区间下界低于 1-α两个方向调整一是把 N 调大重解让 SAA 的估计更准二是把 α 在建模时留出余量比如 0.10 的目标按 0.08 求解样本外检验正好补齐这层保险。从那以后我每次跑完 SAA 都强制走一遍样本外检验即使结果看起来“感觉对了”也照做因为随机优化最骗不了人的就是样本外表现。这套流程配合第 5 章的避坑清单能帮你省掉大量反复试错的无效时间。希望这份代码包和这些经验帮到你。本文还有配套的精品资源点击获取
