yalmip使用

1. yalmip

1.1. 常用流程

1. 创建变量

   x = sdpvar(1, 2, 3, 'full') % 这里全部使用'

   ? 为什么maltab 中’’ 与 “” 不同呢？

2. 优化

   d = optimize(Constraints, Objective, options)

3. 提取求解值

   xopt = value(x);
   yopt = value(y);

4. 创建options

   ops = sdpsettings('solver', 'gurobi');
   ops = sdpsettings(ops, 'verbose', 0) ;% 追加设置

1.2. 创建变量

1. sdpvar
2. binvar
3. intvar

1.3. 优化

diagnostics = optimize(Constraints, Objective, options)

1.3.1. 优化器

提供参数，用于调参

P = optimizer(Con, Obj, Options, Parameters, WantedVariables)

1.4. Constraints

1. 单个约束

2. 多个约束组合,使用[]包括

3. 矩阵约束，使用矩阵的方式

4. 逻辑约束

   1. 指示约束

      `c = implies(x >0 , y==1)`, 如果x >0 ， y 必须等于1

   2. 互斥约束

   c_mutual = xor(x_A > 0 , x_B > 0);

   1. 整数约束

      x_int = intvar(1, 1);
      c_mod = mod(x_int , 2) == 0; % x 必须是整数

1.5. 求解options

options = sdpsettings('solver', 'gurobi'); % 指定求解器
options = sdpsettings(options, 'verbose', 1)% 0 静默，1 详细， 2 更详细
options = sdpsettings(options, 'timeout', 60) % 用s 计算

optiosn = sdpsettings(options, 'FeasibilityTol', 1e-6); % 可行性容忍度
options = sdpsettings('gurobi.OptimalityTol',1e-6);   % 最优性容忍度
options = sdpsettings('gurobi.MIPGap',0.01);          % MIP间隙（1%）

1.6. 求解状态

sol = optimize(Constraits, Objective, options);

xOpt = value(x); % 提取x 
optimal_value = value(objective); % 目标函数

sol.problem % 状态
sol.info % 错误信息

求解状态：

1. 0 ：最优
2. 1： 不可行
3. 2 ： 问题无解
4. 3： 求解超时，返回最优解
5. 4： 内存不足
6. 5： 用户中断
7. 10-100：求解器错误

2. examples

2.1. 工厂排班

问题描述：
某工厂生产两种产品A和B。生产每单位A需要2小时机器时间，3小时人工时间，利润为40元；生产每单位B需要4小时机器时间，2小时人工时间，利润为50元。工厂每天最多有100小时机器时间和80小时人工时间。问如何安排生产计划使利润最大？

1. 变量： A, B 生成量
2. 目标函数：最大化利润
3. 约束条件：
   1. 机器用时
   2. 人工用时

$$max \space 40x_A + 50x_b\\ s.b \space \begin{array} 2x_A + 4x_A \leq 100 ;机器时间\\ 3x_A + 2 x_B \leq 80 \end{array}$$

factory_produce

2.2. 背包问题

你有一个容量为15kg的背包，有5件物品可选，每件物品有重量和价值。问如何选择物品使总价值最大，且总重量不超过15kg？

物品信息：

• 物品1：重量2kg，价值3元，有1个
• 物品2：重量3kg，价值4元 ，有2个
• 物品3：重量4kg，价值5元， 3个
• 物品4：重量5kg，价值8元， 1个
• 物品5：重量9kg，价值10元， 1个

1. 决策变量： $y_I$ 表示选择个数

2. 目标函数：最大价值v

   $max \space v*y$

3. 约束条件：

   1. 重量约束
   2. 个数约束

代码路径：

weight_bag

2.3. 二次规划(QP)

QP matlab

3. gurobi

3.1. 求解参数策略

该章节记住有这些方法即可，使用时再进行查阅，详见Qcode

1. MIPFocus， 求解策略， 2= 可行解， 3= 最优性
2. NodeMethod， 节点松弛方法
3. Heuristics， 启发式搜索强度， 0-1
4. Threads 使用多个线程求解
5. ConcurrentMIP, 并发求解MIP 问题

% MIP策略
options = sdpsettings('gurobi.MIPFocus',2);  % 1=侧重找到可行解，2=侧重证明最优性
options = sdpsettings('gurobi.NodeMethod',2); % 2=barrier方法，用于节点松弛
options = sdpsettings('gurobi.Heuristics',0.1); % 启发式搜索强度(0-1)

% 并行计算
options = sdpsettings('gurobi.Threads',8); % 使用8个线程
options = sdpsettings('gurobi.ConcurrentMIP',2); % 并行MIP求解

2. 收敛性参数

% 精度控制
options = sdpsettings('gurobi.FeasibilityTol',1e-6); % 可行性容忍度
options = sdpsettings('gurobi.OptimalityTol',1e-6);  % 最优性容忍度
options = sdpsettings('gurobi.IntegralityTol',1e-5); % 整数容忍度

% MIP间隙控制
options = sdpsettings('gurobi.MIPGap',0.01);  % 1%的相对间隙
options = sdpsettings('gurobi.MIPGapAbs',0.001); % 绝对间隙

3. 内存与时间管理

% 内存控制， 与分支定界有关
options = sdpsettings('gurobi.NodefileStart',0.5); % 0.5GB后开始写节点文件
options = sdpsettings('gurobi.NodefileDir','/tmp'); % 节点文件目录

% 时间控制
options = sdpsettings('gurobi.TimeLimit',300); % 300秒超时
options = sdpsettings('timeout',300); % YALMIP级别的超时		

4. 预处理与切割平面

% 预处理
options = sdpsettings('gurobi.Presolve',2); % 2=激进预处理

% 切割平面
options = sdpsettings('gurobi.Cuts',2); % 2=激进切割
options = sdpsettings('gurobi.GomoryPasses',2); % Gomory切割强度

3.2. 根据参数规模调优

1. 小规模问题： 1. 精确解

   MIPGap = 1e-4;
   TimeLimit = 60

2. 中等规模问题： 1. 解适度精确；2. 时间放宽

   MIPGap = 1e-2;
   TimeLimite = 300;

3. 大规模问题： 1. 优先找到可行解; 2. 使用启发算法； 3. 增大解经读 ； 4. 增大求解时间

   MIPFocus = 1;
   Heuristics = 0.5;
   MIPGap = 0.05;
   TimeLimit = 1200;

4. 约束密度调整

   约束数量 超过 变量 一定程度

   density = n_constraits/(n_vars^2);
   if density > 0.1
   	Presolve = 2 ; % 使用激进预处理
   end
   

5. MIP 问题：1. 增加线程数量 2. 使用ConcurrentMIP 优化

   threads = min(16, max(1, floor(nCores * 0.8)));
   ConcurrentMIP = 2; % 2 中求解方案并行

4. debug

4.1. infeasible

sol == 1

4.1.1. deBug 步骤

1. 无目标函数，求解一次，查看是否依旧是不可行的
2. 使用方法检测不可行的约束
3. 适当减少部分约束

4.1.2. debug 工具

4.2. unbounded

sol == 2

1. 查看目标函数
2. 人为 增加一个约束，添加合理的上界和下届

5. 编程技巧

5.1. 矩阵

使用循环表示矩阵相乘，替换为矩阵，对于0-1变量时，非常有效

eg : 变量两两互斥，非0 即 1

% 低效：使用嵌套循环
n = 1000;
x = sdpvar(n,1);
Constraints = [];
for i = 1:n
    Constraints = [Constraints, x(i) >= 0];
    for j = 1:n
        if i ~= j
            Constraints = [Constraints, x(i) + x(j) <= 1];
        end
    end
end

$o(n^2)$个约束

使用矩阵

sum(x) <= 1;

$O(1)$ 约束， 简单有效

1. 使用Yalmip 的check 工具

   [infeasible_constraints, ~] = check(Constraints)
   for i = 1: length(infeasible_constraints):
   	if infearible_constraints(i):
   		 fprintf("约束 %d 不可行 \n")
   		 fprintf("输入约束")
   	end
   end

2. 使用Gurobi中IIS

   % 重新求解并要求IIS
   options_iis = sdpsettings('solver','gurobi','gurobi.IISMethod',1);
   sol_iis = optimize(Constraints, Objective, options_iis);
   
   % 获取IIS信息
   if isfield(sol_iis.info, 'IIS')
       fprintf('\nGurobi IIS分析结果:\n');
       iis_constraints = sol_iis.info.IIS;
       for i = 1:length(iis_constraints)
           if iis_constraints(i)
               fprintf('  约束 %d 属于IIS\n', i);
           end
       end
   end

3. 引入松弛变量

   ? 方法原理不懂

   slack = sdpvar(length(Constraints),1);
   Constraints_relaxed = Constraints + diag(sign(slack))*slack;
   Objective_relaxed = sum(abs(slack));
   options_relaxed = sdpsettings('solver','gurobi');
   sol_relaxed = optimize(Constraints_relaxed, Objective_relaxed, options_relaxed);
   
   if sol_relaxed.problem == 0
       slack_values = value(slack);
       fprintf('约束松弛量:\n');
       for i = 1:length(slack_values)
           if abs(slack_values(i)) > 1e-6
               fprintf('  约束 %d 需要松弛 %.4f\n', i, slack_values(i));
           end
       end
   end