Problem with constraint tolerance when using MPCC formulation in GEKKO

Question

I am using GEKKO to solve a nonconvex optimization problem. However, I have encountered an issue with one of the model constraints. This constraint involves a sign function that we have modeled using MPCC. Although the solver (IPOPT) provides an optimal solution, when we evaluate this answer in the constraints with the sign function, we find an error that is greater than the tolerance by a factor of approximately 1E-6, despite the tolerance being set to 1E-7.

The restriction in question is the so-called Weymouth equation, an expression relating the pressures at the ends of a pipeline to the flow through it:

Could this problem be related to the use of MPCC, or could there be another cause?

The constraint is added and evaluated using the weymouth_MPCC and weymouth_eval methods, respectively. Essentially, the numerical values of the solution (pressures and flows) are used to evaluate them in the Weymouth equation, where the sign is no longer modeled with MPCC, but instead, the numpy function is used.

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

To replicate a result I leave the following data. In that case the output obtained has values with orders between 1E-5 and 1E-6 which are larger than the tolerance set in the solver.

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)

The solver tolerance was established as follows:

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',]

score 0 · Answer 1 · answered Apr 02 '23 at 04:46

The Gekko tolerances are for equations, not for the accuracy of an individual variable. There is also scaling applied to the variables based on the initial guess value that may affect the convergence tolerance of each equation.

If you need a more accurate solution based on the signum function, try using the sign3() function instead with the m.options.SOLVER=1 (APOPT solver).

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.sign3(f[i]) * f[i]**2 == p*K)

It may take longer to solve, but the solution is often more accurate because it relies on binary variables instead of MPCCs.

Problem with constraint tolerance when using MPCC formulation in GEKKO

1 Answers1