如何使用Pyomo结合Gurobi建模Frobenius范数矩阵最小二乘问题

问题对应公式参考：

Pyomo原生不支持直接传入矩阵形式的代数表达式，不需要依赖第三方矩阵运算扩展，通过索引逐元素展开表达式的方式即可完成建模，配合Gurobi的非凸求解参数即可完成求解，具体实现流程如下：

实现步骤 #

• 配置Gurobi非凸求解环境
  Gurobi默认关闭非凸二次问题求解能力，必须先手动开启对应参数，否则会直接报非凸模型不可解的错误：

  from pyomo.environ import *
  # 初始化具体模型
  model = ConcreteModel()
  # 调用Gurobi求解器
  solver = SolverFactory('gurobi')
  # 设置非凸求解模式：参数值2代表全局求解非凸二次约束/目标问题
  solver.options['NonConvex'] = 2
  # 可按需设置求解时间上限、收敛gap等参数
  solver.options['TimeLimit'] = 300
  solver.options['MIPGap'] = 1e-4
  

• 用索引集替代矩阵维度定义
  不要直接把numpy矩阵/向量整体传入Pyomo表达式，先把矩阵对应的行、列维度定义为Pyomo原生集合，再按索引绑定参数值：

  # 替换为你自己的问题实际维度
  n = 10  # 决策变量维度
  m = 5   # 线性约束数量
  # 定义索引集合
  model.var_idx = RangeSet(1, n)  # 决策变量索引，对应向量长度
  model.con_idx = RangeSet(1, m)  # 约束索引，对应矩阵行数
  # 绑定系数参数：把原矩阵、向量拆为按索引取值的形式
  # 如果你提前用numpy生成了系数矩阵A（m*n）、向量b（m维），可以直接按索引映射
  # 注意Pyomo RangeSet默认从1开始计数，numpy索引从0开始，需要做偏移
  model.A = Param(model.con_idx, model.var_idx, 
                  initialize=lambda model,j,i: A_np[j-1, i-1])
  model.b = Param(model.con_idx, 
                  initialize=lambda model,j: b_np[j-1])
  

• 定义决策变量
  按索引逐元素定义变量即可，不需要定义向量/矩阵类型的变量：

  # 替换为你自己的变量类型、上下界要求，比如Binary是0-1变量，NonNegativeReals是非负变量
  model.x = Var(model.var_idx, within=Reals, bounds=(-10, 10))
  

• 逐索引展开目标与约束，替代矩阵代数运算
  所有矩阵乘法、二次项、双线性项都用Pyomo原生的sum()函数按索引累加实现，完全规避矩阵运算的依赖：

  以常见的二次目标+线性约束形式举例：若目标为最小化 $x^TQx + c^Tx$，约束为 $Ax \leq b$，展开写法如下

  # 绑定二次项、线性项系数，Q为n*n二次系数矩阵，c为n维线性系数向量
  model.Q = Param(model.var_idx, model.var_idx,
                  initialize=lambda model,i,k: Q_np[i-1, k-1])
  model.c = Param(model.var_idx,
                  initialize=lambda model,i: c_np[i-1])
  
  # 构造目标函数：二次项双重求和+线性项单重求和
  def obj_rule(model):
      quad_part = sum(model.Q[i,k] * model.x[i] * model.x[k] for i in model.var_idx for k in model.var_idx)
      linear_part = sum(model.c[i] * model.x[i] for i in model.var_idx)
      return quad_part + linear_part
  model.obj = Objective(rule=obj_rule, sense=minimize)
  
  # 构造约束：逐行展开矩阵乘法
  def linear_con_rule(model, j):
      return sum(model.A[j,i] * model.x[i] for i in model.var_idx) <= model.b[j]
  model.linear_con = Constraint(model.con_idx, rule=linear_con_rule)
  

• 求解与结果提取
  直接调用求解器即可，Pyomo会自动把展开后的二次模型传递给Gurobi，不需要手动做格式转换：

  # tee=True表示把求解日志打印到控制台
  res = solver.solve(model, tee=True)
  # 提取决策变量的求解结果
  x_solution = [value(model.x[i]) for i in model.var_idx]
  

注意事项 #

• 上述索引展开的写法是Pyomo原生支持的语法，不需要安装任何额外扩展包，生成的模型可以直接被Gurobi识别为二次规划模型，不存在接口兼容问题
• 如果你的问题包含两个不同变量相乘的双线性项，只要保持NonConvex=2的参数设置，Gurobi就可以正常执行全局求解
• 如果模型维度很高，双重求和的写法不会影响求解效率，Pyomo在编译模型时会自动把累加表达式转换成高效的求解器输入格式，和直接传入矩阵的性能没有差异

内容的提问来源于stack exchange，提问作者lyh458