GEKKO中采用MPCC建模时的约束容差异常问题咨询

我使用GEKKO求解非凸优化问题，其中约束涉及符号函数，采用MPCC（数学规划互补约束）建模。尽管求解器IPOPT返回了最优解，但用numpy的符号函数验证Weymouth方程时，发现误差约为1e-6，比设定的1e-7容差高出一个数量级。

Weymouth方程 #

该方程描述管道两端压力与流量的关系，形式如下：
$$|q| \cdot q = K \cdot (p_i^2 - p_j^2)$$
其中$q$为管道流量，$p_i$、$p_j$分别为管道起点和终点压力，$K$为管道特性系数。

问题疑问 #

此误差问题是否与MPCC的使用相关，或是存在其他诱因？

约束建模与验证代码 #

约束的添加与验证分别通过weymouth_MPCC和weymouth_eval方法实现：验证时使用求解得到的压力和流量数值，不再采用MPCC建模符号函数，改用numpy函数。

def weymouth_MPCC(self):
    f = (self.net_flow(self.flow).reshape(-1,))
    press2 = np.array(self.press) ** 2
    for i in range(self.P):
        i_index = (self.pipe['fnode'] - 1)[i]
        j_index = (self.pipe['tnode'] - 1)[i]
        p = press2[i_index] - press2[j_index]
        K = self.Kij[i]
        self.m.Equation(self.m.sign2(f[i]) * f[i]**2 == p*K)

def weymouth_eval(self):
    self.x_sol = []
    for val in self.flow:
        self.x_sol.append(val.value[0])
    f_plus = np.array(self.x_sol[self.W: self.W+self.P])
    f_minus = np.array(self.x_sol[self.W+self.P : self.W+self.P+self.P])
    f_net = f_plus + f_minus
    p2 = np.array([ (p.value[0])**2 for p in self.press])
    self.pipe_f_node = self.pipe['fnode'] - 1
    self.pipe_t_node = self.pipe['tnode'] - 1
    K = self.pipe['Kij'].values
    p2_sub = p2[self.pipe_f_node] - p2[self.pipe_t_node]
    self.wey_eval = (K * (p2_sub)) - np.sign(f_net) * f_net**2
    return self.wey_eval

误差复现代码 #

运行以下代码可复现误差，输出误差值处于1e-5至1e-6量级，超出求解器设定的容差：

press = np.array([1200.5497888, 1196.1393509, 1196.7560686, 1197.5013672, 
                  1196.1935696, 1196.5249472, 1195.6770031, 1192.2243072])
f = np.array([38.63350002, 19.49458569, 19.13891432, 
              -6.92031431, 12.2186    , 26.41490001])
K = np.array([0.1412, 0.1214, 0.1567, 0.0604, 0.0736, 0.0736])
index_i = np.array([0, 3, 3, 4, 5, 4])
index_j = np.array([1, 4, 5, 5, 6, 7])

p2_sub = (press[index_i])**2 - (press[index_j])**2
error = K * p2_sub - np.sign(f) * f**2
print(error)

求解器容差设置 #

self.m.options.OTOL = 1e-7
self.m.options.MAX_ITER = 1e9
self.m.solver_options = ['mu_strategy adaptive',
                         'constr_viol_tol 1e-7',
                         'acceptable_tol 1e-7',]

可能的原因分析 #

1. MPCC近似特性：Gekko的sign2函数通过MPCC实现符号函数，本质是用互补约束做近似，并非严格等价于numpy的np.sign。尤其是当流量$f$接近0时，MPCC的近似解与严格符号函数的计算结果会存在偏差，导致验证误差。
2. 求解器实际收敛精度：虽然设置了1e-7的容差，但IPOPT的约束违反容差（constr_viol_tol）是约束残差的最大值，实际求解中可能因非线性问题的收敛特性，部分约束的残差接近但略低于该值；同时，数值计算中的舍入误差累积也会让验证时的误差被放大。
3. 数值计算精度差异：Gekko求解时采用自动微分，而numpy是数值计算。压力值本身较大（约1200），平方后数值量级达$10^6$，微小的压力差会被放大，两种计算方式的精度累积差异最终导致误差超出预期。

解决建议 #

1. 替换符号函数实现：用Gekko内置的abs函数替代sign2，将约束改写为f[i] * abs(f[i]) == p*K。Gekko的abs同样基于MPCC实现，但对这类非线性约束的处理可能更稳定。
2. 调整求解器参数：适当增大迭代次数（如设置MAX_ITER=10000，1e9的设置无实际意义，IPOPT默认上限远低于此），或设置更严格的求解器容差（如tol 1e-8），尝试获取更高精度的解。
3. 统一计算方式：验证时直接在Gekko模型中添加约束残差计算变量，查看求解器内部的残差值，避免将变量值提取到numpy中计算带来的精度损失。

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