Mathematica中NSolve求解物料衡算方程组调低浓度无可行解问题

使用NSolve求解工厂溶剂浓度调整场景的物料衡算非线性方程组时，参数在特定区间可正常返回符合约束的结果，降低溶剂浓度（将m8流股MgCl₂质量分数从0.056下调到0.05）时无返回，调高浓度到0.06可正常求解，要求所有流量变量为正、质量分数变量在[0,1]区间。

可正常运行的基准代码 #

NSolve[{
  m1 == m4 + m5,
  m1*0.952 == m4 + m5*0.454,
  m1*0.048 == m5*0.546,
  x6MgCl2 + x6H2O == 1,
  0.255*x6MgCl2 == x6Mg,
  0.745*x6MgCl2 == x6Cl2,
  x11H2O + x11Br2 + x11Cl2 == 1,
  m5 + m7 + m9 + m10 + m12 == m6 + m11,
  m5*0.072 == m6*x6Mg,
  m5*0.454 + m7 + m10 == m6*x6H2O + m11*x11H2O,
  m5*0.474 == m11*x11Br2,
  m9 + m12 == m6*x6Cl2 + m11*x11Cl2,
  m6 == m7 + m8,
  m6*x6MgCl2 == m8*0.056,
  m6*x6H2O == m7 + m8*0.944,
  m11 == m12 + m13 + 8163.265,
  m11*x11Br2 == m13*0.030 + 8000,
  m11*x11H2O == m13*0.970 + 163.265,
  m11*x11Cl2 == m12,
  x11Cl2 == x11Br2*0.0311,
  m1 >= 0, m4 >= 0, m5 >= 0, m6 >= 0, m7 >= 0, m8 >= 0, m9 >= 0, 
  m10 >= 0, m11 >= 0, m12 >= 0, m13 >= 0, x11H2O >= 0, x11Br2 >= 0, 
  x11Cl2 >= 0, x6MgCl2 > 0, x6H2O > 0, x6Mg > 0, x6Cl2 > 0
  }, {m1, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, x11H2O, x11Br2, 
  x11Cl2, x6MgCl2, x6H2O, x6Mg, x6Cl2}, Reals]

基准求解结果 #

{{m1 -> 192590., m4 -> 175660., m5 -> 16931., m6 -> 2.0519*10^5, 
  m7 -> 119830., m8 -> 85365., m9 -> 3561.4, m10 -> 73874., 
  m11 -> 9251.3, m12 -> 249.58, m13 -> 838.41, x11H2O -> 0.10556, 
  x11Br2 -> 0.86747, x11Cl2 -> 0.026978, x6MgCl2 -> 0.023297, 
  x6H2O -> 0.9767, x6Mg -> 0.0059407, x6Cl2 -> 0.017356}}

异常参数场景 #

• 无返回的低浓度参数对应方程：

  m6*x6MgCl2 == m8*0.05,
  m6*x6H2O == m7 + m8*0.95,
  

• 可正常求解的高浓度参数对应方程：

  m6*x6MgCl2 == m8*0.06,
  m6*x6H2O == m7 + m8*0.94,
  

解决方案 #

1. 优先消元降低求解规模 #

该方程组存在大量可独立解析求解的子块，提前计算固定值可大幅降低数值求解器的搜索维度：

• 前3个物料衡算方程、x6流股的组分定义方程、x11流股的全部衡算方程均为独立/线性相关子组，其中x11流股的所有变量和m5的取值完全不受浓度参数c影响，可提前单独求解为固定值，不需要放入NSolve反复迭代。

2. 换用适配工程场景的求解器 #

NSolve的设计目标是求解多项式方程组的全部解析/数值解，对于带不等式约束、参数靠近可行域边界的问题收敛性差。工程物料衡算问题优先使用FindRoot，以临近参数的已知解作为初始点，收敛速度和稳定性远高于NSolve。
参考代码如下：

(* 先获取c=0.056的基准解作为初始值 *)
baseSol = NSolve[
   {m1 == m4 + m5, m1*0.952 == m4 + m5*0.454, m1*0.048 == m5*0.546,
    x6MgCl2 + x6H2O == 1, 0.255*x6MgCl2 == x6Mg, 0.745*x6MgCl2 == x6Cl2,
    x11H2O + x11Br2 + x11Cl2 == 1,
    m5 + m7 + m9 + m10 + m12 == m6 + m11,
    m5*0.072 == m6*x6Mg,
    m5*0.454 + m7 + m10 == m6*x6H2O + m11*x11H2O,
    m5*0.474 == m11*x11Br2,
    m9 + m12 == m6*x6Cl2 + m11*x11Cl2,
    m6 == m7 + m8,
    m6*x6MgCl2 == m8*0.056,
    m6*x6H2O == m7 + m8*0.944,
    m11 == m12 + m13 + 8163.265,
    m11*x11Br2 == m13*0.030 + 8000,
    m11*x11H2O == m13*0.970 + 163.265,
    m11*x11Cl2 == m12,
    x11Cl2 == x11Br2*0.0311,
    m1 >= 0, m4 >= 0, m5 >= 0, m6 >= 0, m7 >= 0, m8 >= 0, m9 >= 0,
    m10 >= 0, m11 >= 0, m12 >= 0, m13 >= 0, x11H2O >= 0, x11Br2 >= 0,
    x11Cl2 >= 0, x6MgCl2 > 0, x6H2O > 0, x6Mg > 0, x6Cl2 > 0},
   {m1, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, x11H2O, x11Br2,
    x11Cl2, x6MgCl2, x6H2O, x6Mg, x6Cl2}, Reals][[1]];

(* 设置目标浓度参数c=0.05，传入初始值求解 *)
c = 0.05;
targetSol = FindRoot[
  {m1 == m4 + m5, m1*0.048 == m5*0.546,
   x6MgCl2 + x6H2O == 1, 0.255*x6MgCl2 == x6Mg, 0.745*x6MgCl2 == x6Cl2,
   x11H2O + x11Br2 + x11Cl2 == 1,
   m5 + m7 + m9 + m10 + m12 == m6 + m11,
   m5*0.072 == m6*x6Mg,
   m5*0.454 + m7 + m10 == m6*x6H2O + m11*x11H2O,
   m5*0.474 == m11*x11Br2,
   m9 + m12 == m6*x6Cl2 + m11*x11Cl2,
   m6 == m7 + m8,
   m6*x6MgCl2 == m8*c,
   m6*x6H2O == m7 + m8*(1 - c),
   m11 == m12 + m13 + 8163.265,
   m11*x11Br2 == m13*0.030 + 8000,
   m11*x11H2O == m13*0.970 + 163.265,
   m11*x11Cl2 == m12,
   x11Cl2 == x11Br2*0.0311},
  {{m1, m1 /. baseSol}, {m4, m4 /. baseSol}, {m5, m5 /. baseSol},
   {m6, m6 /. baseSol}, {m7, m7 /. baseSol}, {m8, m8 /. baseSol},
   {m9, m9 /. baseSol}, {m10, m10 /. baseSol}, {m11, m11 /. baseSol},
   {m12, m12 /. baseSol}, {m13, m13 /. baseSol}, {x11H2O, x11H2O /. baseSol},
   {x11Br2, x11Br2 /. baseSol}, {x11Cl2, x11Cl2 /. baseSol},
   {x6MgCl2, x6MgCl2 /. baseSol}, {x6H2O, x6H2O /. baseSol},
   {x6Mg, x6Mg /. baseSol}, {x6Cl2, x6Cl2 /. baseSol}}]

如果需要连续扫描不同浓度参数，可写循环每次将上一个参数的求解结果作为下一个点的初始值，全程不会出现收敛失败。

3. 若需继续使用NSolve，显式收紧变量边界 #

原代码仅给变量设置了下界，NSolve会在极大的无界区间搜索，容易漏掉可行解。可参考基准解的量级给所有变量加上合理的上下界，尤其是质量分数变量明确限制在[0,1]区间，流量变量设置符合工程实际的上界，即可正常返回结果：

c = 0.05;
NSolve[{
  m1 == m4 + m5,
  m1*0.952 == m4 + m5*0.454,
  m1*0.048 == m5*0.546,
  x6MgCl2 + x6H2O == 1,
  0.255*x6MgCl2 == x6Mg,
  0.745*x6MgCl2 == x6Cl2,
  x11H2O + x11Br2 + x11Cl2 == 1,
  m5 + m7 + m9 + m10 + m12 == m6 + m11,
  m5*0.072 == m6*x6Mg,
  m5*0.454 + m7 + m10 == m6*x6H2O + m11*x11H2O,
  m5*0.474 == m11*x11Br2,
  m9 + m12 == m6*x6Cl2 + m11*x11Cl2,
  m6 == m7 + m8,
  m6*x6MgCl2 == m8*c,
  m6*x6H2O == m7 + m8*(1 - c),
  m11 == m12 + m13 + 8163.265,
  m11*x11Br2 == m13*0.030 + 8000,
  m11*x11H2O == m13*0.970 + 163.265,
  m11*x11Cl2 == m12,
  x11Cl2 == x11Br2*0.0311,
  (* 显式设置合理边界 *)
  0 <= m1 <= 1*^6, 0 <= m4 <= 1*^6, 0 <= m5 <= 1*^6, 0 <= m6 <= 1*^6,
  0 <= m7 <= 1*^6, 0 <= m8 <= 1*^6, 0 <= m9 <= 1*^5, 0 <= m10 <= 1*^6,
  0 <= m11 <= 1*^5, 0 <= m12 <= 1*^4, 0 <= m13 <= 1*^4,
  0 <= x11H2O <= 1, 0 <= x11Br2 <= 1, 0 <= x11Cl2 <= 1,
  0 < x6MgCl2 < 1, 0 < x6H2O < 1, 0 < x6Mg < 1, 0 < x6Cl2 < 1
  }, {m1, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, x11H2O, x11Br2, 
  x11Cl2, x6MgCl2, x6H2O, x6Mg, x6Cl2}, Reals]

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