0.1. 第一阶段：环境搭建与认知重构（第1-3天）
0.1.1. 第一天：理解”优化”本质与工具链定位
0.1.1.1. 上午（2小时）：概念认知重塑
1.1 什么是优化问题？（从生活实例讲起）
想象你在超市购物，手里有200元预算，要买10种商品。每种商品有：
- 价格（成本）
- 营养价值（收益）
- 保质期（约束条件）
目标：在预算内买到营养最高的组合，且每种商品至少买1件。
这就是优化问题的三要素：
- 决策变量：每种商品买多少（x₁, x₂, …, x₁₀）
- 目标函数：最大化营养价值（或最小化成本）
- 约束条件：总预算≤200元，每种商品≥1件
数学表达：
maximize: Σ(营养值ᵢ × xᵢ)
subject to: Σ(价格ᵢ × xᵢ) ≤ 200
xᵢ ≥ 1, ∀i ∈ {1,2,...,10}优化求解器（Gurobi）的作用就是：自动帮你找到满足所有条件的最优解，而不是让你穷举所有可能组合。
1.2 工具链角色分工（保姆级比喻）
┌─────────────────────────────────────────────────────────┐
│ 你的办公桌 (MATLAB) │
│ │
│ ┌──────────────┐ ┌──────────────┐ │
│ │ 秘书小姐 │◄─────►│ 数学苦力 │ │
│ │ (YALMIP) │ 翻译 │ (Gurobi) │ │
│ └──────────────┘ └──────────────┘ │
└─────────────────────────────────────────────────────────┘- MATLAB：你的”操作系统”。所有操作都在这里完成，负责数据预处理、结果可视化、脚本管理
- YALMIP：你的”智能翻译秘书”。她懂数学语言（≤, ≥, Σ），会把你的数学公式翻译成Gurobi能听懂的机器语言
- Gurobi：你的”超级计算苦力”。他只做一件事：快速解方程，但需要你告诉他方程长什么样
关键认知：千万不要跳过YALMIP直接学Gurobi！ 就像你不会为了学英语直接去读莎士比亚原著，而是应该先学语法。YALMIP就是那个”语法老师”。
0.1.1.2. 下午（3小时）：软件安装实战
1.3 MATLAB安装与基础配置
【安装步骤】
1. 获取安装包：
- 学校版：通过学校软件中心，使用校园邮箱注册MathWorks账号
- 试用版：官网下载R2024a/R2024b（30天试用）
- 注意：建议使用R2021b及以上版本，对YALMIP兼容性最好
2. 安装注意事项：
- 勾选"MATLAB Compiler"和"Optimization Toolbox"（虽然我们用Gurobi，但这两个工具箱有时辅助用）
- 安装路径不要有中文和空格，例如：D:\MATLAB\R2024b
- 安装后首次启动，选择"激活"，输入许可证或登录MathWorks账号
3. 基础配置（关键！）：
>> cd('D:\MATLAB\R2024b\bin') % 进入安装目录
>> matlab -register % 注册环境变量（Windows）
% 设置默认工作路径（避免每次打开都在C盘）
>> userpath('D:\MyMATLABProjects')
>> savepath1.4 Gurobi学术许可证申请（全程图文）
这是最容易卡壳的环节，务必按步骤操作：
【申请流程（学生版）】
1. 访问 Gurobi 官网：https://www.gurobi.com/
2. 点击 "Free Trial" → "Academic Licenses"
3. 选择 "Gurobi Academic License"（不是Gurobi Cloud！）
4. 填写信息：
- First Name: 你的名（拼音，如 Xiaoming）
- Last Name: 你的姓（拼音，如 Zhang）
- Email: 教育网邮箱（关键！必须是 .edu.cn 或学校邮箱）
- Institution: 你的大学全称
5. 提交后，邮箱会收到验证邮件，点击链接
6. 登录Gurobi账户，进入 "Downloads" 页面
7. 下载对应版本（建议 11.0.3 及以上）
8. 安装时勾选 "Install for all users"（避免权限问题）
【激活许可证】
1. 安装完成后，打开 "Gurobi License Center"
2. 复制你的许可证密钥（格式：xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx）
3. 打开MATLAB，输入：
>> grbgetkey xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx
4. 等待激活成功提示
【验证安装】
在MATLAB命令行输入：
>> g = gurobi()
如果显示 "Gurobi Optimizer version 11.0.3..."，说明安装成功1.5 YALMIP安装（最简单但最重要）
【安装步骤】
1. 访问 YALMIP 官网：http://yalmip.github.io/
2. 点击 "Download" → "Download ZIP"
3. 解压到：D:\MATLAB\R2024b\toolbox\yalmip（推荐路径）
4. 在MATLAB中添加路径：
>> addpath('D:\MATLAB\R2024b\toolbox\yalmip')
>> savepath % 永久保存路径
【验证安装】
在MATLAB命令行输入：
>> yalmiptest预期输出应包含：
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| Test| CPLEX| GUROBI| MOSEK| SDPT3| SEDUMI| FMINCON|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| Linear | OK | OK | OK | OK | OK | OK |
| SOCP | OK | OK | OK | OK | OK | OK |
| InfeasLP| OK | OK | OK | OK | OK | OK |
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++关键：Gurobi栏必须显示”OK”，否则需检查许可证或路径。
0.1.1.3. 晚上（1小时）：环境测试与第一个程序
1.6 运行你的第一个优化问题
% 文件名：test_optimization.m
% 问题描述：minimize x^2 + y^2，满足 x + y = 1
clear all; clc; % 清除变量和命令行
% 步骤1：定义变量
x = sdpvar(1,1); % 实数变量x
y = sdpvar(1,1); % 实数变量y
% 步骤2：设置约束
Constraints = [x + y == 1]; % 注意：== 表示等式约束
% 步骤3：设置目标函数
Objective = x^2 + y^2; % 最小化平方和
% 步骤4：求解设置
options = sdpsettings('solver','gurobi','verbose',1); % 显示求解过程
% 步骤5：求解
diary('solver_log.txt'); % 记录日志
result = optimize(Constraints, Objective, options);
diary off; % 关闭日志
% 步骤6：提取结果
if result.problem == 0 % 求解成功
x_opt = value(x);
y_opt = value(y);
fprintf('最优解: x = %.4f, y = %.4f\n', x_opt, y_opt);
fprintf('最优目标值: %.4f\n', value(Objective));
else
fprintf('求解失败，错误代码: %d\n', result.problem);
end常见问题排查：
- 错误”Undefined function ‘sdpvar’”：YALMIP路径未添加成功，重新addpath
- 错误”Gurobi license error”：许可证未激活或过期，重新运行grbgetkey
- 错误”Infeasible model”：约束矛盾，检查==和<=的使用
0.1.2. 第二天：数学基础补全与YALMIP核心语法
0.1.2.1. 上午（2.5小时）：优化问题的数学分类
2.1 问题的分类体系（必须理解）
优化问题分类树：
优化问题
│
┌────────┴────────┐
线性规划(LP) 非线性规划
│
┌──────────┴──────────┐
二次规划(QP) 一般非线性(NLP)
│
┌─────┴─────┐
凸QP 非凸QP详细解释（逐项拆解）：
① 线性规划 (LP - Linear Programming)
- 数学形式：所有表达式都是线性的（变量的1次方）
minimize: c₁x₁ + c₂x₂ + ... + cₙxₙ subject to: a₁₁x₁ + a₁₂x₂ + ... + a₁ₙxₙ ≤ b₁ a₂₁x₁ + a₂₂x₂ + ... + a₂ₙxₙ ≤ b₂ xᵢ ≥ 0 - YALMIP特点：使用
sdpvar定义变量，目标函数和约束都是线性表达式 - 求解难度：★☆☆☆☆（Gurobi可在毫秒级求解百万变量问题）
② 混合整数线性规划 (MILP - Mixed Integer LP)
- 数学形式：LP + 部分变量必须是整数（通常是0或1）
xᵢ ∈ {0,1} 表示"是否选择i"，或 xᵢ ∈ ℤ 表示"整数个" - YALMIP特点：使用
binvar或intvar定义变量 - 求解难度：★★★★☆（NP-hard问题，Gurobi的最强项）
③ 二次规划 (QP - Quadratic Programming)
- 数学形式：目标函数是二次的，约束是线性的
minimize: xᵀQx + cᵀx (Q是二次项系数矩阵) subject to: Ax ≤ b - YALMIP特点：目标函数可写为
x'*Q*x + c'*x - 求解难度：★★★☆☆（凸QP易解，非凸QP极难）
④ 二阶锥规划 (SOCP - Second Order Cone Programming)
- 数学形式：包含范数约束
||Ax + b||₂ ≤ cᵀx + d - YALMIP特点：可写为
norm(A*x + b) <= c'*x + d - 求解难度：★★★☆☆
2.2 YALMIP核心语法深度解析
语法1：定义变量（sdpvar/binvar/intvar）
% 1. 实数变量（sdpvar）
x = sdpvar(1); % 标量实数变量
x_vec = sdpvar(3,1); % 3×1列向量 [x1;x2;x3]
x_mat = sdpvar(2,3); % 2×3矩阵（每个元素都是变量）
% 2. 二元变量（binvar）- 只能取0或1
y = binvar(1); % 开关变量
y_vec = binvar(5,1); % 5个开关的组合
% 3. 整数变量（intvar）
z = intvar(1); % 整数变量（如购买数量）
z_mat = intvar(3,3); % 3×3整数矩阵
% 4. 特殊形式：对称矩阵变量（用于SDP问题）
X = sdpvar(3,3,'symmetric'); % 3×3对称矩阵，X = X'
% 5. 多维变量：决策变量是三维张量（高阶问题）
W = sdpvar([2,3,4]); % 2×3×4三维变量数组
% 重要技巧：变量命名规范
P_gen = sdpvar(10,1); % 发电机功率（10台机组）
u_on = binvar(10,1); % 机组启停状态（1=运行）语法2：构建约束（Constraints）
% 方式1：逐个添加（适合小型问题）
C1 = [x >= 0]; % 单个约束
C2 = [x + y <= 10]; % 线性约束
C = [C1, C2]; % 约束集合
% 方式2：批量添加（推荐，效率高）
n = 100;
x = sdpvar(n,1);
A = randn(50,n); b = rand(50,1);
C = [A*x <= b, x >= 0]; % 一次性添加50个约束
% 方式3：逻辑约束（需要最新版YALMIP）
y = binvar(1);
C = [implies(y==1, x>=5)]; % 如果y=1，则x必须≥5
% 方式4：绝对值约束（自动转化为线性约束）
C = [abs(x) <= 5]; % 等价于 -5 <= x <= 5
% 方式5：分段约束（分段线性）
C = [-5 <= x <= 5, -3 <= y <= 3];
% 常见错误示例：
% 错误1：使用等号赋值约束
C = [x = 5]; % 错误！必须用 ==
% 错误2：未初始化的约束数组
C = []; % 正确做法：先初始化为空数组
for i = 1:10
C = [C, x(i) <= i]; % 逐个添加
end语法3：设置目标函数（Objective）
% 类型1：最小化（最常用）
obj = c'*x; % 线性目标
obj = x'*Q*x + c'*x; % 二次目标
obj = sum(x.^2); % 平方和
% 类型2：最大化
obj = -c'*x; % 最大化c'*x 等价于 最小化 -c'*x
% 或在求解时指定：
ops = sdpsettings('solver','gurobi','solvefor','max');
% 类型3：多目标（加权和）
obj1 = c1'*x; obj2 = c2'*x;
weights = [0.7, 0.3];
obj = weights(1)*obj1 + weights(2)*obj2;
% 重要：目标函数必须是标量
obj = sum(x); % 正确，标量
obj = x; % 错误，当x是向量时，这是向量目标语法4：求解与结果提取（optimize & value）
% 基础求解
ops = sdpsettings('solver','gurobi');
optimize(C, obj, ops); % 最小化obj
optimize(C, -obj, ops); % 最大化obj
% 高级设置
ops = sdpsettings(...
'solver','gurobi',... % 指定求解器
'verbose',1,... % 显示求解过程
'debug',1,... % 开启调试模式
'saveinput','model.mps'); % 保存模型到文件
% 结果诊断
result = optimize(C, obj, ops);
disp(result); % 显示求解信息
% result.problem 取值含义：
% 0: 求解成功
% 1: 问题无可行解（约束矛盾）
% 2: 问题无界（目标函数可无限优化）
% 3: 达到迭代限制
% 4: 求解器错误
% 提取结果
if result.problem == 0
x_opt = value(x); % 提取最优解
obj_opt = value(obj); % 提取最优目标值
dual_opt = dual(C{1}); % 提取约束的对偶变量（灵敏度分析）
% 格式化输出
fprintf('最优目标值: %.4f\n', obj_opt);
disp('最优解:');
disp(x_opt);
else
% 错误处理
switch result.problem
case 1
error('模型无可行解，请检查约束条件！');
case 2
error('模型无界，请检查目标函数和约束！');
otherwise
error('求解失败，错误代码: %d', result.problem);
end
end0.1.2.2. 下午（2小时）：YALMIP调试技巧与常见错误
2.3 调试工具箱
你的问题非常好！你提供的代码片段反映了你在使用 YALMIP 进行建模和调试，但其中确实 混入了一些非 YALMIP 的函数或用法错误。下面我将：
- 逐条分析你例子中的问题；
- 系统介绍 YALMIP 正确的调试方法；
- 提供修正后的可运行示例。
0.1.3. 一、你例子中的问题分析
0.1.3.1. ❌ analyze(C)
- 问题：YALMIP 没有名为
analyze的函数。 - 正确做法：可用
sdisplay(C)查看约束表达式，或用depends(C)查看涉及的变量，但 没有一键显示变量数、约束数、问题类型的函数。 - 若想查看模型结构，推荐使用：
% 查看约束内容 sdisplay(C) % 查看所有变量 all_vars = allvariables(C, obj); % obj 是目标函数 disp(['Number of variables: ', num2str(length(all_vars))]);
0.1.3.2. ❌ assign(x,5); check(C)
- 问题 1：你定义的是
n = binvar(10...，但后面用assign(x,5)，变量名不一致（x未定义）。 - 问题 2：
check不是 YALMIP 函数。 - 正确做法：
assign(n, some_binary_vector); % 必须是 0/1 向量，长度1000 value(C) % 返回每个约束的数值（如 sum(n) <= 500 的左边值） % 检查是否满足：所有 value(C) <= 0（因为 YALMIP 把约束表示为 <= 0 形式）
0.1.3.3. ❌ feasibility = optimize(C,[],ops);
- 语法问题：YALMIP 的
optimize不能接受空目标（[]）。 - 正确做法：使用
0作为目标（常数）：feasibility = optimize(C, 0, ops); % 仅求可行解
0.1.3.4. ❌ disp(value(obj)); 但 obj 未定义
- 在工具1和工具4中你都没有定义
obj，但工具5却用了它。 - 这会导致运行错误。
0.1.3.5. ❌ saveinput 参数拼写错误
- 正确参数是
'savesolverinput'或'saveduals'，但 YALMIP 没有'saveinput'。 - 若想导出
.lp文件，应使用：
或更推荐用：ops = sdpsettings('solver','gurobi', 'savesolverinput', 'debug.lp');export(C, obj, 'debug.lp'); % 直接导出模型，不求解
0.1.4. 二、YALMIP 正确的调试方法（系统总结）
0.1.4.1. ✅ 1. 检查模型结构
n = binvar(1000,1);
C = [sum(n) <= 500];
obj = sum(n); % 示例目标
% 查看约束表达式
sdisplay(C)
% 查看所有变量
vars = allvariables(C, obj);
fprintf('变量数量: %d (应为1000)\n', length(vars));
% 查看问题是否为混合整数
is_mip = any([vars(:).binvar]) || any([vars(:).intvar]);0.1.4.2. ✅ 2. 导出模型为 LP/MPS 文件
export(C, obj, 'model.lp'); % 不求解，仅导出
% 然后用文本编辑器打开 model.lp0.1.4.3. ✅ 3. 手动赋值并验证约束
% 构造一个可行解：前500个为1，其余为0
test_val = [ones(500,1); zeros(500,1)];
assign(n, test_val);
con_vals = value(C); % 例如返回 -1（表示 sum(n)-500 = -1 <= 0）
fprintf('约束值: %.2f (应 <= 0)\n', con_vals);⚠️ 注意：
assign只用于调试，不会影响 optimize 的求解。求解时 YALMIP 会忽略手动赋值。
0.1.4.4. ✅ 4. 仅求可行解（Feasibility Check）
ops = sdpsettings('verbose', 1); % 显示求解器输出
feas = optimize(C, 0, ops); % 目标为 0
if feas.problem == 0
disp('找到可行解');
else
disp('问题不可行');
end0.1.4.5. ✅ 5. 检查目标函数值（在某赋值下）
assign(n, test_val);
obj_val = value(obj);
disp(['目标值: ', num2str(obj_val)]);0.1.4.6. ✅ 6. 查看求解器返回信息
sol = optimize(C, obj, ops);
if sol.problem ~= 0
disp(['求解失败，错误代码: ', num2str(sol.problem)]);
% 常见：1=不可行，2=无界，3=求解失败
else
optimal_n = value(n);
optimal_obj = value(obj);
end0.1.4.7. ✅ 7. 使用 spy 查看稀疏结构（对大型问题）
% 查看约束矩阵结构（需先构造内部模型）
diagnostic = optimize(C, obj, sdpsettings('solver','gurobi','debug',1));
% 或使用 internal model
[F, f] = model(C, obj); % 不推荐直接用，但可用于分析0.1.5. 三、修正后的完整示例
% 定义变量
n = binvar(1000, 1);
% 约束
C = [sum(n) <= 500];
% 目标（假设最小化 sum(n)）
obj = sum(n);
% 1. 检查模型
fprintf('变量数: %d\n', length(allvariables(C, obj)));
sdisplay(C)
% 2. 导出模型
export(C, obj, 'debug.lp');
% 3. 可行性检查
ops = sdpsettings('solver','gurobi','verbose',0);
feas = optimize(C, 0, ops);
if feas.problem == 0
disp('模型可行');
else
disp('模型不可行！');
end
% 4. 手动赋值验证
assign(n, [ones(300,1); zeros(700,1)]);
con_val = value(C);
obj_val = value(obj);
fprintf('手动赋值下：约束值=%.1f (<=0?), 目标值=%.1f\n', con_val, obj_val);
% 5. 正常求解
sol = optimize(C, obj, ops);
if sol.problem == 0
fprintf('最优目标: %.1f\n', value(obj));
else
error('求解失败');
end0.1.6. 总结
analyze,check不是 YALMIP 函数；- 空目标
[]应替换为0； - 变量名要一致（你用了
n但赋值x）； - 使用
assign+value是调试赋值的有效方法； - 用
export导出.lp文件，比saveinput更可靠； - 始终检查
sol.problem判断求解状态。
如果你有具体的模型不可行、求解慢、结果异常等问题，也可以贴出来，我可以帮你进一步诊断。
2.4 十个新手必犯错误与解决方案
| 错误现象 | 原因分析 | 解决方案 |
|---|---|---|
Undefined function 'sdpvar' |
YALMIP路径未添加 | addpath('.../yalmip'); savepath |
Gurobi license error 10009 |
学术许可证未激活 | 重新运行grbgetkey，检查系统时间 |
Model is infeasible |
约束矛盾 | 用relax(C)或逐个注释约束排查 |
Matrix dimensions must agree |
矩阵维度不匹配 | 用size()检查变量和系数矩阵维度 |
Objective must be scalar |
目标函数是向量 | 使用sum()或square_pos()转换 |
Cannot convert sdpvar to double |
在optimize前使用value() | 只在求解后使用value()提取结果 |
Solver not found |
Gurobi未正确安装 | 运行gurobi_setup，检查MATLAB路径 |
Out of memory |
问题规模太大 | 使用sparse()创建稀疏矩阵 |
Undefined operator '*' |
非线性项错误 | 检查二次项是否写成x'*Q*x形式 |
Warning: Solver not applicable |
问题类型与求解器不匹配 | 检查是否是MILP却用了不支持整数求解器 |
0.1.7. 第三天：问题类型实战与参数调优
0.1.7.1. 上午（2.5小时）：参数调优详解
3.1 Gurobi参数体系（性能调优关键）
% 参数设置结构
ops = sdpsettings(...
'solver','gurobi',...
'gurobi.MIPGap',0.01,... % MIP间隙容忍度（1%）
'gurobi.TimeLimit',3600,... % 时间限制（秒）
'gurobi.Threads',8,... % 使用CPU核数
'gurobi.OutputFlag',1,... % 显示迭代过程
'gurobi.Heuristics',0.5,... % 启发式算法强度
'gurobi.CutPasses',2,... % 割平面法迭代次数
'gurobi.Presolve',2,... % 预求解级别（0-2）
'gurobi.Method',-1,... % 算法选择（-1=自动）
'gurobi.LogFile','gurobi.log');% 日志文件
% 常用参数组合模板：
% 模板1：快速可行解（适合初次运行）
ops_fast = sdpsettings(...
'solver','gurobi',...
'gurobi.MIPFocus',1,... % 重点找可行解
'gurobi.Heuristics',0.8,... % 加强启发式
'gurobi.TimeLimit',300); % 5分钟上限
% 模板2：精确最优解（适合最终求解）
ops_accurate = sdpsettings(...
'solver','gurobi',...
'gurobi.MIPGap',0.001,... % 0.1%精度
'gurobi.MIPFocus',2,... % 重点优化界
'gurobi.Cuts',3,... % 激进割平面
'gurobi.NodeMethod',2); % 节点用对偶单纯形
% 模板3：大规模问题（内存受限）
ops_large = sdpsettings(...
'solver','gurobi',...
'gurobi.Threads',4,... % 限制线程数
'gurobi.NodefileStart',0.5,... % 内存占用50%时转硬盘
'gurobi.Method',2,... % 使用内点法
'gurobi.BarConvTol',1e-4); % 降低精度换速度3.2 参数调优案例：背包问题性能对比
% 问题：100个物品中选择，价值最大化
n = 100;
value = rand(n,1)*100;
weight = rand(n,1)*50;
capacity = 1000;
% 定义变量
x = binvar(n,1);
% 约束
C = [sum(weight.*x) <= capacity];
% 目标
obj = -value'*x; % 最大化价值
% 测试不同参数
param_sets = {
{'MIPGap',0.01,'Threads',4};
{'MIPGap',0.05,'Heuristics',0.8};
{'MIPGap',0.001,'Cuts',3,'MIPFocus',2};
};
for i = 1:length(param_sets)
fprintf('\n===== 参数组合 %d =====\n', i);
ops = sdpsettings('solver','gurobi',...
param_sets{i}{:},'verbose',1,'savesolveroutput',1);
tic;
result = optimize(C, obj, ops);
time = toc;
if result.problem == 0
fprintf('求解时间: %.2f 秒\n', time);
fprintf('目标值: %.2f\n', -value(value(x)));
fprintf('节点数: %d\n', result.solveroutput.IterCount);
end
end0.1.7.2. 下午（2.5小时）：第一阶段总结与实战检验
3.3 第一阶段综合测试题
% 任务：生产排程问题
% 工厂生产A、B两种产品，需要经过机器M1、M2
% 数据：
% - 产品A利润：30元/件，M1耗时2h，M2耗时3h
% - 产品B利润：50元/件，M1耗时4h，M2耗时2h
% - M1可用时间：120h/周，M2可用时间：100h/周
% - 产品A最大销量：25件/周
% 求：利润最大的生产方案
% 请完成以下代码（填空）：
% 1. 定义变量
x_A = ______(1,1); % 产品A产量
x_B = ______(1,1); % 产品B产量
% 2. 约束条件
C = [2*x_A + 4*x_B <= ______, % M1时间约束
3*x_A + 2*x_B <= ______, % M2时间约束
______ <= 25, % 销量约束
x_A >= 0, x_B >= 0];
% 3. 目标函数
Profit = ______*x_A + ______*x_B;
% 4. 求解
ops = sdpsettings('solver','gurobi');
result = optimize(C, -Profit, ops); % 最大化利润
% 5. 结果输出
if result.problem == ____
fprintf('最优生产：A=%.0f件, B=%.0f件\n', ______, ______);
fprintf('最大利润：%.0f元\n', ______);
end答案与解析：
% 答案：
x_A = sdpvar(1,1);
x_B = sdpvar(1,1);
C = [2*x_A + 4*x_B <= 120,
3*x_A + 2*x_B <= 100,
x_A <= 25,
x_A >= 0, x_B >= 0];
Profit = 30*x_A + 50*x_B;
result = optimize(C, -Profit, ops);
if result.problem == 0
fprintf('最优生产：A=%.0f件, B=%.0f件\n', value(x_A), value(x_B));
fprintf('最大利润：%.0f元\n', value(Profit));
end
% 最优解：A=20件, B=20件, 利润=1600元3.4 第一阶段学习检查清单
- [ ] 能独立安装MATLAB、Gurobi、YALMIP
- [ ] 运行yalmiptest显示Gurobi为OK
- [ ] 理解sdpvar/binvar/intvar的区别
- [ ] 能写出至少3种约束条件
- [ ] 会设置目标函数（最小化和最大化）
- [ ] 会调用optimize并提取value
- [ ] 能读懂Gurobi求解日志
- [ ] 遇到”Infeasible”知道如何排查
0.2. 第二阶段：经典模型实战（第2周）
0.2.1. 第4天：线性规划（LP）深度实战
0.2.1.1. 上午（3小时）：线性规划理论体系
4.1 线性规划标准形式与转换
标准形式：
minimize cᵀx
subject to Ax = b
x ≥ 0
实际问题的转换技巧：
1. 不等式转等式：添加松弛变量
aᵢᵀx ≤ bᵢ → aᵢᵀx + sᵢ = bᵢ, sᵢ ≥ 0
2. 自由变量处理：x无约束 → x = x⁺ - x⁻, x⁺,x⁻ ≥ 0
3. 最大化转最小化：
max cᵀx → min -cᵀx4.2 生产排程问题完整案例
% 问题描述：
% 某工厂生产4种产品(P1-P4)，使用3种原材料(R1-R3)
% 数据：
% 利润向量（万元/吨）：[5, 8, 6, 7]
% 原材料消耗矩阵（吨原材料/吨产品）：
% P1 P2 P3 P4
% R1 2 3 1 2
% R2 1 2 3 1
% R3 3 1 2 3
% 原材料可用量（吨）：R1=150, R2=180, R3=200
% 市场需求上限（吨）：P1=30, P2=40, P3=35, P4=25
% MATLAB实现（工业级代码）：
%% 第一步：数据准备（通常从Excel读取）
profit = [5; 8; 6; 7]; % 利润向量
A = [2 3 1 2; % 原材料消耗矩阵
1 2 3 1;
3 1 2 3];
b = [150; 180; 200]; % 原材料上限
demand = [30; 40; 35; 25]; % 市场需求
%% 第二步：定义决策变量
x = sdpvar(4,1); % x(i)表示产品i的产量
%% 第三步：构建约束条件（矩阵化风格）
Constraints = [
A*x <= b, % 原材料约束（3个约束，一行代码）
x >= 0, % 非负约束
x <= demand % 市场需求约束
];
%% 第四步：设置目标函数（最大化利润）
Objective = -profit'*x; % 负号因为YALMIP默认最小化
%% 第五步：求解设置（工业级参数）
ops = sdpsettings(...
'solver','gurobi',...
'verbose',1,...
'gurobi.Method',0,... % 自动选择算法
'gurobi.FeasibilityTol',1e-6,... % 可行性容差
'gurobi.OptimalityTol',1e-6); % 最优性容差
%% 第六步：求解与结果输出
result = optimize(Constraints, Objective, ops);
%% 第七步：结果分析与可视化
if result.problem == 0
x_opt = value(x);
profit_opt = -value(Objective);
% 创建结果表格
fprintf('\n========== 生产排程最优方案 ==========\n');
fprintf('%-10s %-15s %-15s\n', '产品', '产量(吨)', '利润(万元)');
fprintf('----------------------------------------\n');
for i = 1:4
fprintf('%-10s %-15.2f %-15.2f\n', ['P',num2str(i)], x_opt(i), profit(i)*x_opt(i));
end
fprintf('----------------------------------------\n');
fprintf('总利润: %.2f 万元\n', profit_opt);
% 可视化
figure;
subplot(2,1,1);
bar(x_opt);
set(gca,'xticklabel',{'P1','P2','P3','P4'});
title('最优生产方案');
ylabel('产量(吨)');
subplot(2,1,2);
pie(x_opt, {'P1','P2','P3','P4'});
title('产量分布');
else
error('求解失败，错误代码: %d', result.problem);
end4.3 灵敏度分析（Shadow Price）
% 在求解后提取对偶变量
if result.problem == 0
% 原材料约束的对偶变量（影子价格）
shadow_price = dual(Constraints(1:3));
fprintf('\n原材料影子价格（万元/吨）:\n');
fprintf('R1: %.4f\n', shadow_price(1));
fprintf('R2: %.4f\n', shadow_price(2));
fprintf('R3: %.4f\n', shadow_price(3));
% 影子价格>0表示该原材料是瓶颈资源
end0.2.1.2. 下午（3小时）：运输问题与指派问题
4.4 运输问题（Transportation Problem）
% 问题：3个工厂 → 4个仓库，最小化运输成本
% 数据
supply = [100; 150; 200]; % 工厂供应量
demand = [80; 90; 120; 160]; % 仓库需求量
cost = [4 5 6 8; % 单位运输成本矩阵（工厂×仓库）
6 4 3 5;
7 8 5 6];
% 决策变量：x(i,j)表示从工厂i到仓库j的运量
x = sdpvar(3,4,'full'); % 'full'表示非对称矩阵
% 约束
Constraints = [
sum(x,2) == supply, % 行和=供应（每个工厂运出量）
sum(x,1)' == demand, % 列和=需求（每个仓库接收量）
x >= 0 % 非负
];
% 目标
Objective = sum(sum(cost.*x)); % 总成本
% 求解
optimize(Constraints, Objective, sdpsettings('solver','gurobi'));
% 结果展示
x_opt = value(x);
fprintf('\n最优运输方案:\n');
for i = 1:3
for j = 1:4
if x_opt(i,j) > 0.1
fprintf('工厂%d → 仓库%d: %.1f 吨\n', i, j, x_opt(i,j));
end
end
end4.5 指派问题（Assignment Problem）——纯0-1问题
% 问题：5个工人 → 5个任务，最小化总工时
% 数据：工时矩阵
time_cost = [9 2 7 8 2;
6 4 3 7 6;
5 8 1 4 7;
7 6 9 8 5;
2 3 4 5 6];
n = 5;
% 决策变量：x(i,j)=1表示工人i做任务j
x = binvar(n,n); % 5×5二元矩阵
% 约束：每个工人只做1个任务，每个任务只由1个工人做
Constraints = [
sum(x,2) == 1, % 行和=1（每个工人1个任务）
sum(x,1)' == 1 % 列和=1（每个任务1个工人）
];
% 目标
Objective = sum(sum(time_cost.*x));
% 求解（Gurobi在此类问题上性能卓越）
optimize(Constraints, Objective, sdpsettings('solver','gurobi'));
% 结果
x_opt = value(x);
fprintf('\n任务分配方案:\n');
for i = 1:n
for j = 1:n
if x_opt(i,j) > 0.5
fprintf('工人%d ←→ 任务%d (工时:%d)\n', i, j, time_cost(i,j));
end
end
end0.2.2. 第5天：混合整数线性规划（MILP）——Gurobi最擅长的领域
0.2.2.1. 上午（3小时）：MILP核心理论与建模技巧
5.1 为什么MILP这么难？（NP-hard问题）
MILP的求解复杂度：
- 若有n个二元变量，理论上需要检查2ⁿ种组合
- 当n=100时，2¹⁰⁰ ≈ 1.27×10³⁰，比宇宙中原子总数还多
- Gurobi的Branch-and-Cut算法：通过剪枝和割平面，实际只需检查极小一部分组合
Gurobi在MILP上的优势：
1. 先进的预求解技术（Presolve）：可消去30-50%的变量和约束
2. 启发式算法（Heuristics）：快速找到高质量可行解
3. 割平面（Cutting Planes）：强化松弛问题，逼近整数解
4. 并行计算：利用多核CPU加速分支定界5.2 背包问题（Knapsack Problem）——MILP入门
% 问题：容量15kg的背包，选哪些物品使总价值最大
% 数据：10个物品
weights = [2,3,4,5,6,1,2,3,4,5]; % 重量(kg)
values = [3,5,6,8,9,2,3,4,6,7]; % 价值
capacity = 15;
n = length(weights);
% 决策变量：x(i)=1表示选第i个物品
x = binvar(n,1);
% 约束
Constraints = [weights*x <= capacity];
% 目标
Objective = -values'*x; % 最大化价值
% 求解
ops = sdpsettings('solver','gurobi','gurobi.MIPGap',0.01);
result = optimize(Constraints, Objective, ops);
% 结果
x_opt = value(x);
fprintf('\n选中物品: ');
for i = 1:n
if x_opt(i) > 0.5
fprintf('%d ', i);
end
end
fprintf('\n总价值: %.0f\n', -value(Objective));5.3 设施选址问题（Facility Location）——经典MILP
% 问题：在5个候选地建仓库，服务10个客户，最小化总成本（建设+运输）
% 数据
build_cost = [100,150,120,90,110]; % 建设成本（万元）
transport_cost = randn(5,10); % 运输成本矩阵（5仓库×10客户）
demand = ones(10,1)*50; % 每个客户需求50吨
% 决策变量1：y(i)=1表示在i地建仓库
y = binvar(5,1);
% 决策变量2：x(i,j)表示仓库i供应客户j的量
x = sdpvar(5,10);
% 约束
Constraints = [
sum(x,1)' == demand, % 满足所有客户需求
sum(x,2) <= 1000*y, % 只有在建的仓库才能供应（大M法）
x >= 0
];
% 目标
Objective = build_cost'*y + sum(sum(transport_cost.*x));
% 求解（Gurobi会自动处理大M法）
optimize(Constraints, Objective, sdpsettings('solver','gurobi'));0.2.2.2. 下午（3小时）：MILP高级建模技巧
5.4 逻辑约束建模（Indicator Constraints）
% 场景：如果工厂开启(y=1)，则产量x≥100；否则x=0
y = binvar(1); % 开关变量
x = sdpvar(1); % 产量
% 传统大M法（可能数值不稳定）
M = 1000;
Constraints = [
x <= M*y,
x >= 100*y
];
% YALMIP高级语法（自动转化为Indicator）
Constraints = [implies(y==1, x>=100), implies(y==0, x==0)];
% YALMIP会自动转换为Gurobi的indicator约束，数值更稳定5.5 分段线性成本函数建模
% 场景：用电成本分段计价
% 0-100度：0.5元/度
% 100-200度：0.8元/度
% 200度以上：1.2元/度
x = sdpvar(1); % 用电量
cost = sdpvar(1); % 成本
% 方法1：引入辅助变量和约束
x1 = sdpvar(1); x2 = sdpvar(1); x3 = sdpvar(1);
Constraints = [
x == x1 + x2 + x3,
0 <= x1 <= 100,
0 <= x2 <= 100,
0 <= x3 <= inf,
cost == 0.5*x1 + 0.8*x2 + 1.2*x3
];
% 方法2：使用YALMIP的pwlf函数（推荐）
cost = pwlf(x, [0 100 200], [0 50 130]); % 断点[0,100,200]，累计成本[0,50,130]0.2.3. 第6天：二次规划（QP）与投资组合优化
0.2.3.1. 上午（3小时）：QP理论与凸性判断
6.1 二次规划的数学本质
% 标准形式：
% minimize 0.5*x'*Q*x + c'*x
% subject to Ax <= b
% 关键概念：矩阵Q必须是半正定（PSD）才是凸QP
% Gurobi只能高效求解凸QP
% 判断凸性方法：
Q = [2 -1; -1 2]; % 正定矩阵
eig_Q = eig(Q);
if all(eig_Q >= -1e-6) % 容忍小负特征值
disp('是凸QP，Gurobi可解');
else
error('是非凸QP，Gurobi可能无法找到全局最优');
end6.2 投资组合优化（Markowitz模型）——经典QP
% 问题：在10支股票中选择投资组合，期望收益≥12%，风险最小
% 数据准备（通常从CSV/Excel读取）
returns = [0.1;0.12;0.09;0.11;0.13;0.08;0.1;0.09;0.12;0.14]; % 期望收益
cov_matrix = cov(randn(100,10)); % 10×10协方差矩阵
% 决策变量：x(i)表示投资比例
x = sdpvar(10,1);
% 约束
Constraints = [
sum(x) == 1, % 总投资100%
x >= 0, % 不允许做空
returns'*x >= 0.12 % 期望收益≥12%
];
% 目标：最小化风险（方差）
Objective = x'*cov_matrix*x; % 二次项
% 求解（Gurobi对凸QP有专门算法）
ops = sdpsettings('solver','gurobi','gurobi.BarHomogeneous',1);
optimize(Constraints, Objective, ops);
% 结果分析
x_opt = value(x);
risk_opt = sqrt(value(Objective)); % 标准差
return_opt = returns'*x_opt;
fprintf('\n最优投资组合:\n');
for i = 1:10
if x_opt(i) > 0.01 % 只显示占比>1%的股票
fprintf('股票%d: %.2f%%\n', i, x_opt(i)*100);
end
end
fprintf('预期收益: %.2f%%\n', return_opt*100);
fprintf('组合风险: %.4f\n', risk_opt);0.2.3.2. 下午（3小时）：二次约束规划（QCP）与鲁棒优化
6.3 带二次约束的问题
% 场景：投资组合中，要求90%的置信水平下风险不超过15%
% 数据同上
% 约束中引入二次项
VaR_constraint = [x'*cov_matrix*x <= 0.15^2]; % 风险约束
Constraints = [
sum(x) == 1,
x >= 0,
returns'*x >= 0.12,
VaR_constraint
];
% Gurobi会自动识别为QCP问题
optimize(Constraints, -returns'*x, ops); % 最大化收益6.4 鲁棒优化初步（仿射可调鲁棒）
% 场景：收益不确定，在±0.02范围内波动，求最坏情况下的最优解
% 定义不确定集
delta = sdpvar(10,1); % 收益波动
Gamma = 2; % 预算不确定参数（控制保守度）
% 鲁棒约束
Robust_Constraints = [
sum(x) == 1,
x >= 0,
(returns+delta)'*x >= 0.12, % 收益约束需满足所有delta
sum(abs(delta)) <= Gamma, % 不确定预算约束
-0.02 <= delta <= 0.02 % 波动范围
];
% YALMIP会自动将鲁棒约束转化为确定性等价形式
optimize(Robust_Constraints, -min((returns+delta)'*x), ops);0.2.4. 第7天：第二阶段综合实战与复盘
0.2.4.1. 全天（6小时）：完整项目——微电网调度
%% 微电网经济调度问题（MILP+QP混合）
% 系统组成：
% - 柴油发电机：成本函数 0.1P² + 5P + 100（二次成本）
% - 光伏：发电成本0，但不确定
% - 电池：充放电效率90%，容量限制
% - 负荷：必须满足
% 预测数据（24小时）
load_profile = [50 45 40 38 35 40 55 70 85 90 88 85 82 80 78 75 80 88 92 85 78 70 62 55];
solar_forecast = [0 0 0 0 0 5 15 30 45 60 70 75 80 75 65 50 35 20 8 0 0 0 0 0];
%% 参数设置
n = 24; % 24小时调度
P_diesel = sdpvar(n,1); % 柴油发电机功率
P_battery = sdpvar(n,1); % 电池充放电（正=放电，负=充电）
P_grid = sdpvar(n,1); % 主网购售电（正=购电，负=售电）
SoC = sdpvar(n+1,1); % 电池荷电状态（0-100%）
u_diesel = binvar(n,1); % 柴油机启停状态
% 技术参数
P_diesel_max = 100; % 柴油机最大功率
P_battery_max = 50; % 电池功率限制
SoC_min = 0.2; SoC_max = 0.9; % SoC限制
SoC_initial = 0.5; % 初始SoC
battery_eff = 0.9; % 充放电效率
diesel_cost_coeff = [0.1, 5, 100]; % 成本系数[a,b,c]对应aP²+bP+c
grid_buy_price = 0.8; % 购电价格
grid_sell_price = 0.5; % 售电价格
%% 约束条件
Constraints = [];
% 1. 功率平衡约束（最关键的等式约束）
for t = 1:n
Constraints = [Constraints,
P_diesel(t) + P_battery(t) + solar_forecast(t) + P_grid(t) == load_profile(t)];
end
% 2. 柴油机约束
Constraints = [Constraints,
0 <= P_diesel <= P_diesel_max*u_diesel, % 不运行时功率为0
u_diesel <= 1, u_diesel >= 0
];
% 3. 电池约束（动力学模型）
Constraints = [Constraints,
-P_battery_max <= P_battery <= P_battery_max,
SoC(1) == SoC_initial % 初始荷电状态
];
for t = 1:n
Constraints = [Constraints,
SoC(t+1) == SoC(t) - P_battery(t)/P_battery_max*battery_eff,
SoC(t+1) >= SoC_min,
SoC(t+1) <= SoC_max
];
end
% 4. 主网约束（设双向功率限制）
P_grid_max = 200;
Constraints = [Constraints, -P_grid_max <= P_grid <= P_grid_max];
%% 目标函数（分段成本）
% 柴油机成本 = 0.1*P² + 5*P + 100*u
diesel_cost = diesel_cost_coeff(1)*sum(P_diesel.^2) + ...
diesel_cost_coeff(2)*sum(P_diesel) + ...
diesel_cost_coeff(3)*sum(u_diesel);
% 电网交互成本（购电正成本，售电负成本）
grid_cost = grid_buy_price*sum(max(P_grid,0)) - grid_sell_price*sum(max(-P_grid,0));
% 总成本
Objective = diesel_cost + grid_cost;
%% 求解
ops = sdpsettings('solver','gurobi', 'gurobi.TimeLimit',600, 'verbose',1);
result = optimize(Constraints, Objective, ops);
%% 结果分析与绘图
if result.problem == 0
% 提取结果
P_diesel_opt = value(P_diesel);
P_battery_opt = value(P_battery);
P_grid_opt = value(P_grid);
SoC_opt = value(SoC(1:n)); % 取前24小时的
% 绘制调度曲线
figure('Position',[100 100 1200 800]);
subplot(3,1,1);
plot(1:n, load_profile, 'k-', 'LineWidth',2); hold on;
plot(1:n, P_diesel_opt, 'r--');
plot(1:n, solar_forecast, 'y:');
plot(1:n, P_battery_opt, 'b-.');
plot(1:n, P_grid_opt, 'g--');
legend('负荷','柴油机','光伏','电池','主网');
title('功率调度曲线');
xlabel('时间(h)'); ylabel('功率(kW)');
grid on;
subplot(3,1,2);
stairs(1:n, value(u_diesel), 'r-', 'LineWidth',2);
title('柴油机启停状态');
ylabel('状态'); ylim([-0.1 1.1]); grid on;
subplot(3,1,3);
plot(0:n, value(SoC), 'b-o', 'LineWidth',2);
title('电池荷电状态');
xlabel('时间(h)'); ylabel('SoC(%)');
grid on;
% 成本分析
total_cost = value(Objective);
fprintf('\n总运行成本: %.2f 元\n', total_cost);
fprintf('柴油机成本: %.2f 元\n', value(diesel_cost));
fprintf('电网交互成本: %.2f 元\n', value(grid_cost));
end0.3. 第三阶段：进阶技巧与参数调优（第3周）
0.3.1. 第8天：大规模问题处理技巧
0.3.1.1. 上午（3小时）：稀疏矩阵与向量化编程
8.1 为什么不能用for循环？（性能对比）
% 场景：创建1000个发电机的出力约束
n = 1000;
P = sdpvar(n,1);
P_max = rand(n,1);
%% 方法1：for循环（不推荐）
tic;
Constraints = [];
for i = 1:n
Constraints = [Constraints, P(i) <= P_max(i)];
end
toc; % 耗时约 0.5-2秒，内存占用高
%% 方法2：向量化（推荐）
tic;
Constraints = [P <= P_max]; % 一行代码
toc; % 耗时约 0.01秒，内存占用低
% 性能对比：
% n=1000时，方法2快50-200倍
% n=10000时，方法1可能内存溢出，方法2依然流畅8.2 稀疏矩阵在电力系统中的应用
%% 场景：IEEE 118节点系统（大规模潮流约束）
% 导入系统数据（导纳矩阵）
load('case118.mat'); % 包含bus, branch, gen等结构体
nb = 118; % 节点数
ng = 54; % 发电机数
% 决策变量
P_gen = sdpvar(ng,1);
V = sdpvar(nb,1); % 电压幅值
% 构建节点功率平衡（稀疏方式）
% 传统稠密矩阵：A = zeros(nb,nb); 占用内存 118*118*8 ≈ 105KB
% 稀疏矩阵：A = sparse(nb,nb); 占用内存 < 10KB
% 导纳矩阵G和B通常是稀疏的
G = sparse(gen_G_matrix); % 电导矩阵
B = sparse(gen_B_matrix); % 电纳矩阵
% 功率注入表达式（向量化计算）
P_inj = V.*(G*V) + V.*(B*V); % 节点注入功率
% 负荷向量
P_load = bus(:,3); % 第3列是负荷
% 功率平衡（稀疏约束）
Constraints = [P_inj == P_load + A_gen*P_gen]; % A_gen是稀疏的关联矩阵8.3 参数化的约束生成
%% 场景：不同场景下的鲁棒优化
% 定义场景数
n_scenarios = 100;
P_load = sdpvar(n_scenarios,1); % 每个场景下的负荷
% 基础约束
Constraints = [];
for s = 1:n_scenarios
% 从数据文件读取场景s的负荷
load_data = load(['scenario_',num2str(s),'.mat']);
% 参数化约束：根据不同场景调整约束右端项
Constraints = [Constraints,
P_gen(s) + P_grid(s) == load_data.P_load];
end
% 等价的高效写法（如果需要）：
% 使用cell数组预先生成所有约束
constraint_cell = cell(n_scenarios,1);
for s = 1:n_scenarios
constraint_cell{s} = (P_gen(s) + P_grid(s) == load_data.P_load);
end
Constraints = [constraint_cell{:}];0.3.1.2. 下午（3小时）：求解器日志解读与性能调优
8.4 Gurobi日志深度解读
典型日志输出分析：
Optimize a model with 100 rows, 200 columns and 500 nonzeros
Model fingerprint: 0x12345678
Variable types: 150 continuous, 50 integer (50 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+03]
Objective range [1e-01, 2e+02]
Bounds range [1e+00, 1e+00]
RHS range [1e+02, 1e+03]
Presolve removed 25 rows and 50 columns (presolve time = 0.05s)...
【解读】
- 100 rows = 100个约束
- 200 columns = 200个变量
- 500 nonzeros = 500个非零系数（稀疏度=500/(100*200)=2.5%）
- Presolve removed：预求解消去了25个冗余约束和50个变量
- 50 binary：50个二元变量 → MILP问题
Presolve time: 0.12s
Presolved: 75 rows, 150 columns, 375 nonzeros
Explored 0 nodes (0 simplex iterations) in 0.05 seconds
Thread count was 8 (of 8 available processors)
Solution count 2: 1234.5 1567.8
Optimal solution found (tolerance 1.00e-04)
Best objective 1.234500000000e+03, best bound 1.234500000000e+03, gap 0.0000%
【解读】
- Explored 0 nodes：根节点就找到了最优解（幸运！）
- Solution count 2：找到2个可行解
- gap 0.0000%：最优性间隙为0，证明全局最优
- Best objective = Best bound：上下界重合，最优性得证8.5 性能调优实战案例
%% 场景：混合整数规划求解太慢
% 问题：机组组合（Unit Commitment），1440个二元变量
% 目标：24小时调度，最小化成本
% 初始设置（慢，需600秒）
ops_slow = sdpsettings('solver','gurobi');
% optimize(...) % 600秒
% 调优后设置（快，仅需45秒）
ops_fast = sdpsettings(...
'solver','gurobi',...
'gurobi.MIPGap',0.005,... % 接受0.5%间隙
'gurobi.Heuristics',0.8,... % 加强启发式
'gurobi.VarBranch',3,... % 最大冲突分支策略
'gurobi.CutPasses',5,... % 增加割平面迭代
'gurobi.Presolve',2,... % 最高级预求解
'gurobi.NumFocus',1,... % 侧重找到可行解
'gurobi.Threads',12); % 使用更多线程
% 性能提升原理：
% 1. MIPGap放宽：从默认0.01%到0.5%，减少枚举
% 2. Heuristics加强：更快找到优质初始解
% 3. VarBranch策略：选择对目标影响大的变量分支
% 4. Cuts增加：收紧松弛问题，更接近整数解
% 5. 并行加速：12线程 vs 4线程，理论加速3倍（实际约2.5倍）0.3.2. 第9天：不可行性诊断与模型修复
0.3.2.1. 上午（3小时）：冲突约束检测
9.1 使用conflict函数定位矛盾
%% 场景：模型无可行解，但不知道哪个约束冲突
% 构造一个故意矛盾的模型
x = sdpvar(1);
C1 = [x >= 5];
C2 = [x <= 3];
C3 = [x == 4]; % 与前两个矛盾
Constraints = [C1, C2, C3];
% 求解
result = optimize(Constraints, x, sdpsettings('solver','gurobi'));
% 检测冲突
if result.problem == 1 % 无可行解
[conflict_constraints, conflict_flags] = conflict(Constraints);
% conflict_flags显示哪些约束参与冲突
disp('冲突约束索引:');
find(conflict_flags)
% 详细分析
analyze(conflict_constraints)
end
% 输出示例：
% Conflict set size 2
% Constraint 1: x >= 5
% Constraint 2: x <= 3
% 说明这两个约束直接矛盾，C3是冗余的9.2 软约束（Soft Constraints）与罚函数
%% 场景：负荷平衡约束可能无法满足，允许少量切负荷
P_gen = sdpvar(10,1);
P_load = ones(10,1)*50; % 预测负荷
% 硬约束（可能导致无解）
% Constraints = [sum(P_gen) == sum(P_load)]; % 可能不可行
% 软约束（允许偏差，但惩罚）
violation = sdpvar(1); % 偏差变量
penalty_weight = 1000; % 惩罚系数
% 修改约束
Constraints = [
sum(P_gen) + violation == sum(P_load), % 允许偏差
violation >= 0
];
% 目标函数增加惩罚项
Base_Cost = sum(P_gen);
Penalty_Cost = penalty_weight * violation;
Objective = Base_Cost + Penalty_Cost;
% 求解后会优先满足硬约束，软约束尽量满足
optimize(Constraints, Objective, sdpsettings('solver','gurobi'));0.3.2.2. 下午（3小时）：数值稳定性处理
9.3 数据缩放（Scaling）
%% 场景：模型系数范围太大导致数值问题
% 原始数据（极差巨大）
coeff = [1e-6, 1e3, 1e-3, 1e6]; % 矩阵范围[1e-6,1e6]
% Gurobi会警告：Matrix coefficient range exceeds 1e+06
% 解决方法：缩放
scale_factor = 1e3;
coeff_scaled = coeff / scale_factor; % 新范围[1e-9,1e3]
% 求解后再反缩放结果
% 或在YALMIP中自动处理：
ops = sdpsettings('solver','gurobi','gurobi.ScaleFlag',1);
% ScaleFlag=1启用自动缩放9.4 大M法的问题与改进
%% 场景：逻辑约束中的大M值选择
y = binvar(1);
x = sdpvar(1);
% 错误做法：M过大导致数值不稳定
M_bad = 1e6;
Constraints_bad = [x <= M_bad*y];
% 正确做法：M尽量紧
% 先估算x的最大可能值
x_max = 100; % 基于业务逻辑估算
Constraints_good = [x <= x_max*y];
% 最佳做法：使用YALMIP的implies（自动避免大M）
Constraints_best = [implies(y==0, x==0), implies(y==1, x>=0)];
% YALMIP会生成Gurobi的indicator约束，无需手动大M0.3.3. 第10天：结果后处理与可视化
0.3.3.1. 上午（2小时）：提取对偶变量与灵敏度分析
10.1 对偶变量（Shadow Price）的经济学意义
%% 场景：电力系统节点电价计算
% 节点功率平衡约束
P_inj = sdpvar(n,1); % 节点注入
P_load = bus(:,3); % 节点负荷
Gen_incidence = ...; % 发电机-节点关联矩阵
% 约束
Constraints = [P_inj == P_load + Gen_incidence*P_gen];
% 求解
optimize(Constraints, Cost, sdpsettings('solver','gurobi'));
% 提取对偶变量（节点电价）
if result.problem == 0
nodal_price = dual(Constraints); % 每个节点的影子价格
% 对偶变量表示该节点每增加1MW负荷，总成本增加量（$/MW）
% 绘制电价地图
figure;
scatter(bus(:,2), bus(:,3), 100, nodal_price, 'filled');
colorbar;
title('节点电价分布');
xlabel('经度'); ylabel('纬度');
end10.2 参数灵敏度分析
%% 场景：天然气价格变化对机组调度影响
gas_price_range = 2:0.5:5; % 天然气价格范围（$/MMBtu）
total_cost_curve = zeros(size(gas_price_range));
for i = 1:length(gas_price_range)
% 重新计算目标函数（系数含gas_price）
gas_price = gas_price_range(i);
Fuel_Cost = gas_price * Heat_Rate * P_gen;
% 求解
optimize(Constraints, Fuel_Cost, ops);
if result.problem == 0
total_cost_curve(i) = value(Fuel_Cost);
end
end
% 绘制成本曲线
plot(gas_price_range, total_cost_curve, '-o');
xlabel('天然气价格 ($/MMBtu)');
ylabel('总运行成本 ($)');
title('价格灵敏度分析');
grid on;0.3.3.2. 下午（2小时）：结果可视化高级技巧
10.3 Gantt图（调度问题必备）
%% 场景：5台机组24小时启停计划
% 求解后得到u_on矩阵（5台×24小时）
u_on_opt = value(u_on); % 5×24的二元矩阵
% 绘制Gantt图
figure('Position',[100 100 1200 400]);
y_tics = 1:5;
for i = 1:5
on_periods = find(u_on_opt(i,:)==1);
if ~isempty(on_periods)
for t = on_periods
x_start = t-1;
y_bottom = i-0.3;
rectangle('Position',[x_start,y_bottom,1,0.6],...
'FaceColor',[0.2 0.6 0.9],...
'EdgeColor','none');
end
end
end
xlim([0 24]);
ylim([0.5 5.5]);
set(gca,'YTick',y_tics,'YTickLabel',{'G1','G2','G3','G4','G5'});
xlabel('时间 (h)');
ylabel('机组');
title('机组启停计划甘特图');
grid on;10.4 3D曲面图（展示多参数影响）
%% 场景：电池容量 vs 光伏容量 对系统成本影响
battery_cap_range = 0:10:100; % 电池容量范围
solar_cap_range = 0:5:50; % 光伏容量范围
cost_matrix = zeros(length(battery_cap_range), length(solar_cap_range));
for i = 1:length(battery_cap_range)
for j = 1:length(solar_cap_range)
% 重新设置容量参数
battery_capacity = battery_cap_range(i);
solar_capacity = solar_cap_range(j);
% 更新约束（电池容量约束）
Constraints = [SoC <= battery_capacity/1000, ...];
% 求解
optimize(Constraints, Objective, ops);
if result.problem == 0
cost_matrix(i,j) = value(Objective);
end
end
end
% 绘制3D曲面
[B,S] = meshgrid(solar_cap_range, battery_cap_range);
surf(B,S,cost_matrix/1e6);
xlabel('光伏容量 (MW)');
ylabel('电池容量 (MWh)');
zlabel('年化成本 (百万元)');
title('容量配置优化分析');
colorbar;
shading interp; % 平滑着色0.4. 第四阶段：综合项目实战（第4周）
0.4.1. 第11天：电力系统——机组组合问题（Unit Commitment）
0.4.1.1. 全天（6小时）：IEEE 30节点系统UC问题
%% IEEE 30节点系统机组组合问题（完整工业级代码）
clear; clc; close all;
%% 1. 系统数据加载
% 数据包括：机组参数、负荷曲线、网络参数
load('case30data.mat'); % 假设已准备好
n_gen = 6; % 6台发电机
n_hour = 24; % 24小时调度
% 机组参数（通常从Excel读取）
gen_param = [
10 20 100 50 0.1 0.2 5 100; % [Pmin Pmax a b c start_cost min_down min_up]
20 50 150 80 0.08 0.15 8 120;
15 35 120 60 0.09 0.18 6 110;
12 40 130 70 0.11 0.22 7 130;
18 45 140 75 0.12 0.24 9 140;
25 60 160 90 0.13 0.26 10 150
];
P_min = gen_param(:,1);
P_max = gen_param(:,2);
cost_a = gen_param(:,3);
cost_b = gen_param(:,4);
cost_c = gen_param(:,5);
start_cost = gen_param(:,6);
min_down = gen_param(:,7);
min_up = gen_param(:,8);
% 24小时负荷曲线（MW）
load_curve = 100*[2.4 2.3 2.2 2.1 2.0 2.1 2.3 2.8 3.2 3.5 3.6 3.4 3.2 3.0 3.1 3.3 3.5 3.8 3.6 3.4 3.2 2.9 2.6 2.5];
%% 2. 决策变量定义
P = sdpvar(n_gen, n_hour); % 发电功率
u = binvar(n_gen, n_hour); % 启停状态
v_start = binvar(n_gen, n_hour); % 启动变量
v_shutdown = binvar(n_gen, n_hour); % 停机变量
%% 3. 约束条件
Constraints = [];
% 3.1 功率平衡（每个小时）
for t = 1:n_hour
Constraints = [Constraints,
sum(P(:,t)) == load_curve(t)];
end
% 3.2 机组出力限制
for t = 1:n_hour
for g = 1:n_gen
Constraints = [Constraints,
P(g,t) >= P_min(g)*u(g,t),
P(g,t) <= P_max(g)*u(g,t)
];
end
end
% 3.3 最小启停时间约束（逻辑复杂，需仔细）
% 假设初始状态：已运行3小时
initial_status = ones(n_gen,1);
initial_runtime = 3*ones(n_gen,1);
% 启动/停机变量定义
for g = 1:n_gen
for t = 1:n_hour
if t == 1
Constraints = [Constraints,
v_start(g,t) >= u(g,t) - initial_status(g),
v_shutdown(g,t) >= initial_status(g) - u(g,t)
];
else
Constraints = [Constraints,
v_start(g,t) >= u(g,t) - u(g,t-1),
v_shutdown(g,t) >= u(g,t-1) - u(g,t)
];
end
Constraints = [Constraints,
v_start(g,t) <= u(g,t),
v_shutdown(g,t) <= 1 - u(g,t)
];
end
end
% 最小运行/停机时间
for g = 1:n_gen
for t = 1:n_hour
% 最小运行时间
for tau = max(1,t-min_up(g)+1):t
Constraints = [Constraints,
u(g,t) <= sum(v_start(g,tau:t)) + (initial_runtime(g) >= min_up(g))];
end
% 最小停机时间
for tau = max(1,t-min_down(g)+1):t
Constraints = [Constraints,
1-u(g,t) <= sum(v_shutdown(g,tau:t)) + (initial_status(g)==0)];
end
end
end
%% 4. 目标函数
Generation_Cost = 0;
Start_Cost = 0;
for g = 1:n_gen
for t = 1:n_hour
Generation_Cost = Generation_Cost + ...
cost_c(g)*P(g,t)^2 + cost_b(g)*P(g,t) + cost_a(g)*u(g,t);
Start_Cost = Start_Cost + start_cost(g)*v_start(g,t);
end
end
Objective = Generation_Cost + Start_Cost;
%% 5. 求解设置与运行
ops = sdpsettings(...
'solver','gurobi',...
'verbose',1,...
'gurobi.MIPGap',0.01,... % 1%间隙
'gurobi.TimeLimit',1800,... % 30分钟
'gurobi.MIPFocus',1,... % 重点找可行解
'gurobi.Heuristics',0.8); % 强化启发式
% 添加初始解（暖启动）
% 假设全部机组以最小出力运行
P_init = repmat(P_min, 1, n_hour);
u_init = ones(n_gen, n_hour);
ops = sdpsettings(ops,'gurobi.Start',struct('P',P_init,'u',u_init));
% 求解
tic;
result = optimize(Constraints, Objective, ops);
solve_time = toc;
%% 6. 结果分析与验证
if result.problem == 0
P_opt = value(P);
u_opt = value(u);
v_start_opt = value(v_start);
total_cost = value(Objective);
fprintf('\n========== 机组组合求解成功 ==========\n');
fprintf('求解时间: %.2f 秒\n', solve_time);
fprintf('总成本: $%.2f\n', total_cost);
% 绘制调度曲线
figure('Position',[100 100 1400 800]);
% 发电堆叠图
subplot(2,1,1);
area(1:n_hour, P_opt', 'LineStyle','none');
hold on;
plot(1:n_hour, load_curve, 'k-', 'LineWidth',2);
legend('G1','G2','G3','G4','G5','G6','负荷');
title('机组出力曲线');
xlabel('时间(h)'); ylabel('功率(MW)');
grid on;
% 启停状态图
subplot(2,1,2);
imagesc(1:n_hour, 1:n_gen, u_opt);
set(gca,'YTick',1:n_gen,'YTickLabel',{'G1','G2','G3','G4','G5','G6'});
xlabel('时间(h)'); ylabel('机组');
title('启停状态（蓝色=运行）');
colorbar;
% 成本构成分析
gen_cost_val = value(Generation_Cost);
start_cost_val = value(Start_Cost);
figure;
pie([gen_cost_val, start_cost_val], {'发电成本','启停成本'});
title('成本构成');
else
error('机组组合求解失败，错误代码: %d', result.problem);
end0.4.2. 第12天：控制工程——模型预测控制（MPC）
0.4.2.1. 全天（6小时）：MPC完整实现
%% 线性MPC控制器设计（以双积分系统为例）
% 系统模型：x(k+1) = A*x(k) + B*u(k)
% 状态：x = [位置;速度]
% 控制输入：u = 加速度
%% 1. 系统参数
Ts = 0.1; % 采样时间(s)
A = [1 Ts; 0 1]; % 离散状态矩阵
B = [0.5*Ts^2; Ts]; % 离散输入矩阵
nx = 2; % 状态数
nu = 1; % 控制输入数
% 约束
u_min = -5; u_max = 5; % 加速度限制
x_max = [10; 5]; % 位置和速度限制
% MPC参数
N = 20; % 预测时域
Q = eye(nx); % 状态权重
R = 0.1; % 输入权重
%% 2. 决策变量定义
X = sdpvar(nx, N+1); % 预测时域内所有状态
U = sdpvar(nu, N); % 预测时域内所有输入
x0 = sdpvar(nx,1); % 当前状态（已知参数）
%% 3. MPC模型约束
Constraints = [];
% 初始条件
Constraints = [X(:,1) == x0];
% 系统动力学（对每个预测步）
for k = 1:N
Constraints = [Constraints,
X(:,k+1) == A*X(:,k) + B*U(:,k)];
end
% 状态约束
for k = 1:N+1
Constraints = [Constraints,
-x_max <= X(:,k) <= x_max];
end
% 输入约束
for k = 1:N
Constraints = [Constraints,
u_min <= U(:,k) <= u_max];
end
%% 4. 目标函数
Objective = 0;
for k = 1:N
Objective = Objective + X(:,k)'*Q*X(:,k) + R*U(:,k)^2;
end
% 终端代价
Objective = Objective + X(:,N+1)'*Q*X(:,N+1);
%% 5. 创建MPC求解器（离线编译）
% 将x0作为参数
Controller = optimizer(Constraints, Objective, sdpsettings('solver','gurobi'), x0, U(:,1));
%% 6. 仿真闭环控制
sim_time = 50; % 仿真步数
x_current = [5; 0]; % 初始状态[位置;速度]
x_history = zeros(nx, sim_time);
u_history = zeros(nu, sim_time);
for t = 1:sim_time
% 求解MPC（在线计算）
u_opt = Controller{x_current};
% 应用第一个控制输入
u_current = value(u_opt);
% 记录
x_history(:,t) = x_current;
u_history(:,t) = u_current;
% 模拟系统演化（真实系统）
x_current = A*x_current + B*u_current;
% 可视化
clf;
subplot(2,1,1);
plot(1:t, x_history(1,1:t), 'b-', 'LineWidth',2);
hold on;
plot(t, x_current(1), 'ro');
title('位置轨迹'); xlabel('步数'); ylabel('位置');
grid on;
subplot(2,1,2);
stairs(1:t, u_history(:,1:t), 'r-', 'LineWidth',2);
title('控制输入'); xlabel('步数'); ylabel('加速度');
grid on;
ylim([u_min u_max]);
drawnow;
pause(0.1);
end
%% 7. 性能分析
figure;
subplot(2,1,1);
plot(x_history(1,:), x_history(2,:), 'b-', 'LineWidth',2);
title('状态相轨迹');
xlabel('位置'); ylabel('速度');
grid on;
subplot(2,1,2);
plot(u_history, 'r-', 'LineWidth',2);
title('控制输入时序');
xlabel('步数'); ylabel('加速度');
grid on;0.4.3. 第13天：路径规划与旅行商问题（TSP）
0.4.3.1. 全天（6小时）：TSP建模与求解
%% 旅行商问题（TSP）—— 经典组合优化
% 问题：访问10个城市，每个城市一次，返回起点，路程最短
% 城市坐标
cities = [
0 0; 10 0; 5 8.66; 15 8.66; 20 0;
25 8.66; 30 0; 35 8.66; 40 0; 45 8.66
];
n = size(cities,1);
% 计算距离矩阵D(i,j)
D = zeros(n,n);
for i = 1:n
for j = 1:n
D(i,j) = norm(cities(i,:) - cities(j,:));
end
end
%% 1. Miller-Tucker-Zemlin (MTZ) 建模法
% 决策变量：x(i,j)=1表示路径i→j
x = binvar(n,n,'full'); % 'full'表示非对称
% 辅助变量：u(i)用于消除子回路
u = sdpvar(n,1);
% 约束
Constraints = [];
% 每个城市恰好一次入边和出边
Constraints = [Constraints, sum(x,1)' == 1]; % 出度=1
Constraints = [Constraints, sum(x,2) == 1]; % 入度=1
% 消除子回路（MTZ约束）
M = n; % 大M
for i = 2:n
for j = 2:n
if i ~= j
Constraints = [Constraints, u(i) - u(j) + M*x(i,j) <= M-1];
end
end
end
% 固定起点（城市1）
Constraints = [Constraints, x(1,1) == 0]; % 不自我连接
% 目标函数
Objective = sum(sum(D.*x));
%% 2. 求解（TSP是MILP，Gurobi强项）
ops = sdpsettings(...
'solver','gurobi',...
'verbose',1,...
'gurobi.MIPGap',0.001,... % 0.1%精度
'gurobi.Heuristics',0.8,... % 强化启发式
'gurobi.CutPasses',5); % 加强割平面
result = optimize(Constraints, Objective, ops);
%% 3. 结果可视化
if result.problem == 0
x_opt = value(x);
% 提取路径
tour = zeros(1,n+1);
tour(1) = 1; % 起点
for i = 2:n+1
next_city = find(x_opt(tour(i-1),:) > 0.5);
tour(i) = next_city;
end
% 绘制路径图
figure('Position',[100 100 800 600]);
plot(cities(:,1), cities(:,2), 'ko', 'MarkerSize',10, 'LineWidth',2);
hold on;
% 绘制路线
plot(cities(tour,1), cities(tour,2), 'b-', 'LineWidth',2);
plot([cities(tour(end),1), cities(tour(1),1)],...
[cities(tour(end),2), cities(tour(1),2)], 'b-', 'LineWidth',2);
% 标注城市
for i = 1:n
text(cities(i,1), cities(i,2)+1, sprintf('C%d',i),...
'HorizontalAlignment','center');
end
title(sprintf('TSP最优路径 (距离=%.2f)', value(Objective)));
xlabel('X坐标');
ylabel('Y坐标');
grid on;
axis equal;
% 绘制时间序列（如果是动态TSP）
figure;
plot(1:n+1, tour, 'b-o', 'LineWidth',2);
yticks(1:n);
title('访问顺序');
xlabel('步骤');
ylabel('城市编号');
grid on;
else
error('TSP求解失败');
end
%% 4. 对比：使用Gurobi的预热启动（Warm Start）
% 提供初始可行解（如最近邻启发式）
initial_tour = nearest_neighbor_tour(cities); % 自定义函数
x_initial = zeros(n,n);
for k = 1:n-1
x_initial(initial_tour(k), initial_tour(k+1)) = 1;
end
x_initial(initial_tour(n), initial_tour(1)) = 1;
% 设置初始解
ops_warm = sdpsettings(ops, 'gurobi.Start', struct('x',x_initial));
% 再次求解
tic;
result_warm = optimize(Constraints, Objective, ops_warm);
time_warm = toc;
fprintf('\n预热启动效果:\n');
fprintf('原求解时间: %.2f秒\n', result.solverinfo.solvetime);
fprintf('预热后时间: %.2f秒\n', time_warm);
fprintf('加速比: %.2f%%\n', (1-time_warm/result.solverinfo.solvetime)*100);0.4.4. 第14天：考试与项目答辩
0.4.4.1. 全天（6小时）：模拟项目答辩
14.1 答辩题目（随机抽取）
题目1：微电网调度优化
- 需求：设计一个含5台柴油发电机、2个储能、光伏和负荷的微电网日前调度模型
- 约束：功率平衡、机组爬坡、储能动力学、旋转备用
- 目标：最小化运行成本
- 要求：提交代码、结果曲线、经济性分析报告
题目2：无人机路径规划
- 需求：5架无人机从基地出发，访问20个目标点，返回基地
- 约束：每架无人机航程<50km，避障（禁飞区），时间窗
- 目标：最小化总航程
- 要求：3D可视化路径，冲突检测分析
题目3：智能制造排程
- 需求：10台机器，20个订单，每个订单多道工序
- 约束：机器能力、工序先后、交付时间窗
- 目标：最小化延迟交付惩罚
- 要求：Gantt图，机器利用率分析14.2 评分标准
| 维度 | 权重 | 评分细则 |
|---|---|---|
| 模型正确性 | 40% | 约束完整、逻辑正确、变量定义合理 |
| 代码质量 | 30% | 结构化、注释清晰、矩阵化（无for循环） |
| 结果分析 | 20% | 可视化专业、灵敏度分析、经济性评估 |
| 创新性 | 10% | 额外功能（如鲁棒性、多目标） |
14.3 项目报告模板
%% 项目报告生成脚本（自动生成Word/PDF）
% 1. 设置报告模板
import mlreportgen.report.*;
import mlreportgen.dom.*;
rpt = Report('optimization_project','pdf');
% 2. 添加标题页
tp = TitlePage();
tp.Title = '微电网经济调度优化项目';
tp.Author = '学习者';
tp.Date = date;
add(rpt,tp);
% 3. 添加问题描述
chapter1 = Chapter('问题描述');
add(chapter1, Paragraph('本研究考虑一个含...'));
add(rpt, chapter1);
% 4. 添加数学模型
chapter2 = Chapter('数学模型');
eq1 = FormalEq('$min \sum_{t=1}^{24} \sum_{g=1}^{6} (a_g P_{g,t}^2 + b_g P_{g,t} + c_g u_{g,t})$');
add(chapter2, eq1);
add(rpt, chapter2);
% 5. 添加代码
chapter3 = Chapter('MATLAB实现');
code = ProgrammaticCode(which('microgrid_uc.m'));
add(chapter3, code);
add(rpt, chapter3);
% 6. 添加结果图
chapter4 = Chapter('结果分析');
fig = Figure();
fig.Snapshot = 'results_figure1.png';
add(chapter4, fig);
add(rpt, chapter4);
% 7. 生成报告
close(rpt);
rptview(rpt);0.5. 第五阶段：避坑指南与进阶学习（第4周补充）
0.5.1. 第15天：新手必踩的100个坑（精选TOP 20）
0.5.1.1. 上午（3小时）：代码级错误
TOP 1-5：变量类型混淆
% ❌ 错误1：数值变量和sdpvar混用
x = sdpvar(1);
if x > 0 % 错误！x是符号变量，不能用于逻辑判断
y = 1;
end
% ✅ 正确：用YALMIP的逻辑约束
y = binvar(1);
Constraints = [implies(x>=0, y==1), implies(x<0, y==0)];
% ❌ 错误2：赋值与约束混淆
x = sdpvar(1);
x = 5; % 这是赋值，不是约束！x不再是sdpvar
% ✅ 正确：
Constraints = [x == 5];
% ❌ 错误3：使用非线性函数
x = sdpvar(1);
obj = sin(x); % 错误！YALMIP不支持sin
% ✅ 正确：用分段线性近似
obj = pwlf(x, [-pi -pi/2 0 pi/2 pi], [-1 0 1 0 -1]); % 近似sin(x)
% ❌ 错误4：矩阵乘法维度错误
x = sdpvar(3,1);
A = rand(2,3);
obj = x'*A*x; % 错误！x'是1×3，A是2×3
% ✅ 正确：
obj = x'*(A'*A)*x; % A'*A是3×3
% 或
obj = sum((A*x).^2);
% ❌ 错误5：未初始化的约束数组
Constraints = []; % 正确初始化
for i = 1:10
Constraints = [Constraints, x(i) <= i]; % 逐个添加
end
% 如果忘记初始化：Constraints未定义，报错TOP 6-10：求解器设置错误
% ❌ 错误6：求解器名称拼写错误
ops = sdpsettings('solver','gurobi'); % 正确
ops = sdpsettings('solver','Gurobi'); % 错误！大小写敏感
% ❌ 错误7：参数名称错误
ops = sdpsettings('solver','gurobi','MIPGap',0.01); % 错误！应加前缀
% ✅ 正确：
ops = sdpsettings('solver','gurobi','gurobi.MIPGap',0.01);
% ❌ 错误8：设置了求解器不支持的参数
ops = sdpsettings('solver','gurobi','barrier_iterations',100);
% Gurobi没有barrier_iterations参数，正确的是'gurobi.BarIterLimit'
% ❌ 错误9：时间单位混淆
ops = sdpsettings('solver','gurobi','gurobi.TimeLimit',100);
% 单位是秒！设置为100以为100分钟，实际是100秒
% ❌ 错误10：verbose和OutputFlag混淆
sdpsettings('verbose',0) % 抑制YALMIP输出
sdpsettings('gurobi.OutputFlag',0) % 抑制Gurobi输出
% 两者都要设才能完全静默0.5.1.2. 下午（3小时）：建模级错误
TOP 11-15：约束建模错误
% ❌ 错误11：严格不等式
Constraints = [x > 0]; % 错误！MILP不支持严格不等式
% ✅ 正确：
epsilon = 1e-6;
Constraints = [x >= epsilon]; % 或x >= 0（如果允许零）
% ❌ 错误12：双向约束写错
% 意图：10 <= x <= 20
Constraints = [10 <= x <= 20]; % MATLAB解析错误
% ✅ 正确：
Constraints = [10 <= x, x <= 20];
% ❌ 错误13：非凸二次约束
x = sdpvar(2,1);
Constraints = [x'*[1 0;0 -1]*x >= 0]; % 非凸！Gurobi报错
% 解决方案：用Cplex或Bonmin求解器，或改写约束
% ❌ 错误14：循环约束索引错误
x = sdpvar(10,1);
Constraints = [];
for i = 1:10
Constraints = [Constraints, x(i) <= x(i+1)]; % i+1=11越界！
end
% ✅ 正确：
for i = 1:9
Constraints = [Constraints, x(i) <= x(i+1)];
end
% ❌ 错误15：使用绝对值导致非线性
x = sdpvar(1);
obj = abs(x); % 这不是线性！
% ✅ 正确：引入辅助变量
y = sdpvar(1);
Constraints = [y >= x, y >= -x];
obj = y; % 最小化y等价于最小化|x|TOP 16-20：算法与理论错误
% ❌ 错误16：QP问题Q矩阵非正定
Q = [1 2; 2 1]; % 特征值=-1,3，非正定
x = sdpvar(2,1);
obj = x'*Q*x; % Gurobi警告"Q matrix is not positive semi-definite"
% 解决方案：Q = 0.5*(Q+Q') + epsilon*eye(n); 强制对称正定
% ❌ 错误17：未检查求解状态直接使用结果
result = optimize(C,obj,ops);
x_opt = value(x); % 如果求解失败，value返回NaN
% ✅ 正确：
if result.problem == 0
x_opt = value(x);
else
error('求解未成功，结果无效');
end
% ❌ 错误18：MILP问题设置了连续算法
ops = sdpsettings('solver','gurobi','gurobi.Method',0); % 对MILP，Method只在根节点有效
% 对MILP，主要控制分支定界的是'MIPFocus'、'Heuristics'等
% ❌ 错误19：混合精度问题
% 数据：load = 100.000000000001（CSV导入产生的微小误差）
% 约束：sum(gen) == load
% 结果：可行解可能因数值误差被判为不可行
% ✅ 正确：
load = round(load, 6); % 保留6位小数
% 或设置容差：
ops = sdpsettings('gurobi.FeasibilityTol',1e-5);
% ❌ 错误20：未使用暖启动重复求解
% 场景：滚动优化中，每次从头求解
for k = 1:100
optimize(C,obj,ops); % 每次都冷启动，慢！
end
% ✅ 正确：利用前次解
x_previous = [];
for k = 1:100
if ~isempty(x_previous)
ops.Start = x_previous; % 暖启动
end
result = optimize(C,obj,ops);
x_previous = value(x); % 保存当前解
end0.5.2. 第16天：进阶主题与前沿技术
0.5.2.1. 上午（3小时）：鲁棒优化与分布鲁棒
16.1 经典鲁棒优化（Box Uncertainty）
%% 场景：需求不确定，在±10%范围内波动
P_load_nominal = 100; % 标称负荷
uncertainty = 0.1; % 10%波动
% 传统确定性优化
P_gen = sdpvar(1);
Constraints_det = [P_gen == P_load_nominal];
optimize(Constraints_det, P_gen^2);
% 鲁棒优化（最坏情况）
delta = sdpvar(1); % 不确定性变量
Constraints_robust = [
P_gen >= P_load_nominal + delta, % 必须覆盖最坏情况
-uncertainty*P_load_nominal <= delta <= uncertainty*P_load_nominal
];
optimize(Constraints_robust, P_gen^2); % 最大化最坏情况下的性能16.2 分布鲁棒优化（Wasserstein球）
%% 场景：基于历史数据的分布鲁棒
% 假设有100个历史负荷场景
load_scenarios = load_data_hist; % 100×1向量
% Wasserstein球半径
epsilon = 0.1;
% 定义概率分布变量
p = sdpvar(100,1);
Constraints = [
p >= 0,
sum(p) == 1
];
% 构造分布鲁棒约束（较复杂，需对偶转换）
% YALMIP的distribute函数可简化
% 此处为概念性代码0.5.2.2. 下午（3小时）：多目标优化与帕累托前沿
16.3 加权和方法
%% 场景：成本 vs 排放 双目标
Cost = ...; % 经济性目标
Emission = ...; % 环境目标
% 加权和
weights = linspace(0,1,11); % 11个权重组合
pareto_front = zeros(2,length(weights));
for i = 1:length(weights)
Objective = weights(i)*Cost + (1-weights(i))*Emission;
optimize(Constraints, Objective, ops);
pareto_front(1,i) = value(Cost);
pareto_front(2,i) = value(Emission);
end
% 绘制帕累托前沿
plot(pareto_front(1,:), pareto_front(2,:), '-o');
xlabel('成本');
ylabel('排放');
title('帕累托前沿');
grid on;16.4 epsilon约束法
%% 更精确的方法：epsilon约束
cost_target = linspace(min_cost, max_cost, 20);
emission_min = zeros(size(cost_target));
for i = 1:length(cost_target)
Constraints_epsilon = [Constraints, Cost <= cost_target(i)];
optimize(Constraints_epsilon, Emission, ops);
emission_min(i) = value(Emission);
end
plot(cost_target, emission_min, '-x');
xlabel('成本上限');
ylabel('最小排放');
title('epsilon约束法帕累托前沿');0.5.3. 第17天：工业界最佳实践
0.5.3.1. 上午（3小时）：代码工程化
17.1 函数封装与模块化
%% 规范的优化项目结构
% 项目文件夹结构：
% microgrid_project/
% ├── main.m % 主脚本
% ├── data/
% │ ├── case30.mat % 系统数据
% │ └── load_forecast.csv % 负荷预测
% ├── src/
% │ ├── build_model.m % 模型构建
% │ ├── solve_model.m % 求解器调用
% │ ├── post_process.m % 后处理
% │ └── visual_results.m % 可视化
% └── results/
% ├── figures/
% └── logs/
%%% build_model.m 示例
function [Constraints, Objective, vars] = build_model(system_data, params)
% 构建微电网调度模型
% 输入：
% system_data - 系统数据结构体
% params - 参数结构体
% 输出：
% Constraints - 约束集合
% Objective - 目标函数
% vars - 变量结构体
% 解包数据
n_gen = system_data.n_gen;
n_hour = params.n_hour;
% 定义变量
vars.P_gen = sdpvar(n_gen, n_hour);
vars.u_on = binvar(n_gen, n_hour);
% 构建约束（此处简化，实际更复杂）
Constraints = [];
Constraints = [Constraints, vars.P_gen >= 0];
Constraints = [Constraints, vars.P_gen <= system_data.P_max*vars.u_on];
% 目标函数
cost_coeff = system_data.cost_coeff;
Objective = sum(sum(cost_coeff(1)*vars.P_gen.^2 + ...
cost_coeff(2)*vars.P_gen + ...
cost_coeff(3)*vars.u_on));
end17.2 错误处理与日志系统
%%% solve_model.m 示例（工业级错误处理）
function [result, vars_opt] = solve_model(Constraints, Objective, vars, ops)
% 求解模型并处理错误
% 创建日志文件
log_file = sprintf('solver_log_%s.txt', datestr(now,'yyyymmdd_HHMMSS'));
diary(log_file);
try
% 求解
result = optimize(Constraints, Objective, ops);
% 检查求解状态
switch result.problem
case 0
fprintf('求解成功！\n');
% 提取结果
vars_opt = extract_results(vars);
case 1
warning('模型无可行解');
vars_opt = [];
% 调用冲突分析
conflict_analysis(Constraints);
case 2
warning('模型无界');
vars_opt = [];
otherwise
error('求解失败，错误代码: %d', result.problem);
end
catch ME
% 捕获异常
fprintf('求解过程中发生异常: %s\n', ME.message);
vars_opt = [];
result.problem = -1;
end
diary off;
end
%%% 辅助函数
function vars_opt = extract_results(vars)
% 递归提取所有变量的值
fields = fieldnames(vars);
for i = 1:length(fields)
field = fields{i};
if isa(vars.(field), 'sdpvar')
vars_opt.(field) = value(vars.(field));
elseif isstruct(vars.(field))
vars_opt.(field) = extract_results(vars.(field));
end
end
end0.5.3.2. 下午（3小时）：性能优化与并行计算
17.3 YALMIP + Parallel Computing Toolbox
%% 场景：多场景并行优化（场景间独立）
% 创建并行池
parpool(8); % 启动8个worker
% 定义场景
scenarios = 1:50;
% 并行求解
results = parallel_fcn(@solve_single_scenario, scenarios);
% 关闭并行池
delete(gcp);
%%% 辅助函数
function result = solve_single_scenario(scenario_id)
% 每个worker独立求解一个场景
% 加载场景数据
load(['scenario_',num2str(scenario_id),'.mat']);
% 构建并求解模型（与主程序隔离）
Constraints = ...;
Objective = ...;
ops = sdpsettings('solver','gurobi','gurobi.OutputFlag',0); % 抑制输出
result = optimize(Constraints, Objective, ops);
end17.4 Gurobi分布式优化（Tuning Tool）
%% 使用Gurobi Tuning Tool自动调参
% 保存模型
savemodel(Constraints, Objective, 'model.mps');
% 调用gurobi调参工具
% 命令行：gurobi_cl TuneTimeLimit=3600 TuneCriterion=2 model.mps
% 在MATLAB中：
system('gurobi_cl TuneTimeLimit=3600 TuneCriterion=2 model.mps');
% 调参后生成.prm文件，加载使用
ops = sdpsettings('solver','gurobi','gurobi.Params','tuned.prm');0.6. 附录：学习资源与工具链
0.6.1. A-1 必备资源清单
【软件安装】
- MATLAB R2024b（官网：mathworks.com）
- Gurobi 11.0（官网：gurobi.com，学术许可证免费）
- YALMIP（GitHub：yalmip/yalmip）
【文档教程】
- YALMIP Wiki：yalmip.github.io（最好的教材，所有函数说明）
- Gurobi Documentation：www.gurobi.com/documentation（参数详解）
- MATLAB Optimization Toolbox：mathworks.com/help/optim/
【书籍推荐】
- 《Convex Optimization》（Boyd & Vandenberghe）- 理论基石
- 《Optimization in Engineering》（Kirby）- 工程应用
- 《Model Predictive Control》（Rawlings & Mayne）- MPC圣经
【在线课程】
- Coursera: Discrete Optimization（含MILP）
- edX: Optimization Methods for Business Analytics
- B站：搜索"YALMIP教程"（有中文视频）
【论文复现】
- 电力系统：《Unit commitment in electric power systems》（Wood & Wollenberg）
- 路径规划：《A survey on the vehicle routing problem》
- MPC：《Model predictive control: theory and practice》
【社区支持】
- Stack Overflow: 标签[yalmip], [gurobi]
- YALMIP Google Groups: 提问和查看历史问题
- GitHub Issues: 报告bug和请求功能0.6.2. A-2 自测题与答案
%% 自测题1：判断模型类型
% 根据以下目标函数和约束，判断问题类型（LP/MILP/QP/QCP）
% (1)
% minimize x1 + 2*x2
% s.t. x1 + x2 >= 5
% x1, x2 >= 0
% 答案：LP（线性目标，线性约束，连续变量）
% (2)
% minimize x1^2 + x2^2
% s.t. x1 + x2 >= 5
% x1, x2 ∈ {0,1}
% 答案：非凸QP（二元变量+二次目标）→ Gurobi可处理但可能慢
% (3)
% minimize x1 + y1
% s.t. x1 + 2*y1 >= 3
% x1 >= 0, y1 ∈ {0,1}
% 答案：MILP（线性混合整数）
% (4)
% minimize x1'*Q*x1 + c'*x1 + y1'*d
% s.t. ||x1||_2 <= 10
% y1 ∈ {0,1}
% 答案：MIQCP（混合整数二次约束）
%% 自测题2：代码纠错
% 以下代码有3处错误，请找出并修正
x = sdpvar(5,1);
y = binvar(5,1);
Constraints = [];
for i = 1:5
Constraints = [x(i) + y(i) <= 10]; % 错误1：覆盖而非追加
end
obj = sum(x.^2 + y);
ops = sdpsettings('solver','gurobi','MIPGap',0.01); % 错误2：参数名
optimize(Constraints, obj, ops);
disp('最优值:', value(obj)); % 错误3：语法
% 修正后：
x = sdpvar(5,1);
y = binvar(5,1);
Constraints = [];
for i = 1:5
Constraints = [Constraints, x(i) + y(i) <= 10]; % 追加
end
obj = sum(x.^2 + y);
ops = sdpsettings('solver','gurobi','gurobi.MIPGap',0.01); % 加前缀
result = optimize(Constraints, obj, ops);
if result.problem == 0
fprintf('最优值: %.4f\n', value(obj)); % 正确语法
end
%% 自测题3：建模练习
% 建模问题：仓库选址，最小化建设+运输成本
% 数据：
% - 5个候选仓库位置，建设成本[100,150,120,90,110]万元
% - 10个客户，需求量[30,20,15,25,18,22,30,12,28,35]吨
% - 运输成本：距离(km)×0.5元/吨km
% 要求：写出完整MATLAB代码
% 答案框架：
clear;
% 数据
build_cost = [100,150,120,90,110];
demand = [30,20,15,25,18,22,30,12,28,35];
distance = rand(5,10)*50; % 随机生成距离
transport_rate = 0.5;
% 变量
y = binvar(5,1); % 选址
x = sdpvar(5,10); % 运输量
% 约束
C = [];
C = [C, sum(x,2) <= 1000*y]; % 只有建的仓库才能供应
C = [C, sum(x,1)' == demand]; % 满足需求
C = [C, x >= 0];
% 目标
Build_Cost = build_cost'*y;
Transport_Cost = sum(sum(distance*transport_rate.*x));
Objective = Build_Cost + Transport_Cost;
% 求解
optimize(C, Objective, sdpsettings('solver','gurobi'));
disp(value(Objective));0.7. 结语：从入门到精通
恭喜你完成了这个万字学习路线图！回顾四周的学习：
- 第1周：环境搭建 + LP/MILP/QP基础 → 掌握YALMIP语法
- 第2周：经典模型实战 → 理解建模思维
- 第3周：参数调优 + 调试技巧 → 解决实际问题
- 第4周：综合项目 + 避坑指南 → 具备独立项目能力
下一步建议：
- 立即行动：明天就开始第1天的环境安装，不要拖延
- 项目驱动：结合你的研究方向，1个月内完成一个真实项目
- 代码复用：保存所有写过的代码，建立个人函数库
- 持续学习：关注YALMIP GitHub更新，每年Gurobi新版本都有性能提升
- 社区贡献：在Stack Overflow回答YALMIP问题，教学相长
记住：优化建模是艺术 + 科学的结合。多练习、多调试、多总结，你很快就能成为团队中的”优化专家”！
祝学习顺利！
文档信息：
- 版本：v2.5（2025年12月更新）
- 作者：基于YALMIP官方文档与工业实践整理
- 许可证：CC BY-NC-SA 4.0（非商业用途可自由分享）
