CODEWITHSUNDEEP
pythonconstraintspyomo
Christos Antonak
Christos Antonak

Reputation: 19

Error evaluating constraint in Pyomo Modelling Language

Hello to the community,

So I have the following error:

Error evaluating constraint 5965: can't evaluate sqrt'(0). ERROR:pyomo.opt:Solver (asl) returned non-zero return code (1) ERROR:pyomo.opt:See the solver log above for diagnostic information.

Solver (asl) did not exit normally

Actually the problem is in the following part of code: 

def thermal_lim_ineq_con(model, i, j, k):
    if line_loading_limit == True:
        P_branch_ij = (model.V_magn[i, k]**2) * np.float(real(Y_matrix[i, j])) - model.V_magn[i, k] * model.V_magn[j, k] * (np.float(real(Y_matrix[i, j])) * cos(model.V_angle[i, k]-model.V_angle[j, k]) + np.float(imag(Y_matrix[i,j])) * sin(model.V_angle[i, k]-model.V_angle[j, k]))
        Q_branch_ij = (model.V_magn[i, k]**2) * np.float(imag(Y_matrix[i, j])) + model.V_magn[i, k] * model.V_magn[j, k] * (np.float((real(Y_matrix[i, j]))) * sin(model.V_angle[i, k]-model.V_angle[j, k]) - np.float(imag(Y_matrix[i,j])) * cos(model.V_angle[i, k]-model.V_angle[j, k]))
        return (sqrt((P_branch_ij**2) + (Q_branch_ij**2)) <= limits_flows[i, j]) ########### ATTENTION
    else: 
        return Constraint.Skip 

model.therlim_ineq_con = Constraint(branch_from_to, scenario_set, rule = thermal_lim_ineq_con) # run this constraint for all branches

Any ideas?? Removing the sqrt and squaring the variable limits_flows on the right, solves this problem but then the limits are very low ( order -6, -7) and the problem gets infeasible.

Thank you.

Edit:

Finally the problem could be solved, by giving random initial conditions to 

model.V_magn[i, k]

Upvotes: 0

Views: 2188

Answers (1)

jsiirola
jsiirola

Reputation: 2589

This is technically not a problem with Pyomo – but rather a problem with your formulation.

The error is being tossed by the ASL (part of the solver executable) when the solver drives the values of P_branch_ij and Q_branch_ij to 0 (this may be happening because of an inactive line, or bad initial values). At zero, the derivative of sqrt(0) is undefined. This numerical error causes the solver to exit without returning a solution.

There are several possible workarounds:

  1. As you pointed out in your edit, you can try different initial variable values. When the solver starts from a different starting point, it may not take a path that drives P_branch_ij and Q_branch_ij to zero. The downside of this approach is that you may have to try many different initial values before getting one that does not have the solver pass through sqrt(0), and you are not guaranteed to find one that works (especially if the solution to the problem really is at 0!).

  2. Bound the variables so that the total branch power (P_branch_ij**2 + Q_branch_ij**2) never hits zero. This is a bit trickier, as you cannot do it with a simple constraint: most solvers are allowed to violate constraints during the solution process and the fatal error is if the solver ever attempts to evaluate the derivative of sqrt(0). Fortunately, most solvers will always respect variable bounds. Therefore, to prevent the solver from ever evaluating sqrt'(0) you would need to define intermediate Pyomo variables (not just Python expressions) for PQ2_branch_ij, and then set the variable lower bound to something greater than 0. The downside of this approach is you are now solving a different problem, and may force yourself to a suboptimal solution if the solution to the original problem really was at 0.

model.PQ2_branch_ij = Var(branch_from_to, scenario_set, bounds=(1e-8, None))
def compute_PQ2_branch_ij(m, i, j, k):
    P_branch_ij = (m.V_magn[i, k]**2) * np.float(real(Y_matrix[i, j])) - m.V_magn[i, k] * m.V_magn[j, k] * (np.float(real(Y_matrix[i, j])) * cos(m.V_angle[i, k]-m.V_angle[j, k]) + np.float(imag(Y_matrix[i,j])) * sin(m.V_angle[i, k]-m.V_angle[j, k]))
    Q_branch_ij = (m.V_magn[i, k]**2) * np.float(imag(Y_matrix[i, j])) + m.V_magn[i, k] * m.V_magn[j, k] * (np.float((real(Y_matrix[i, j]))) * sin(m.V_angle[i, k]-m.V_angle[j, k]) - np.float(imag(Y_matrix[i,j])) * cos(m.V_angle[i, k]-m.V_angle[j, k]))
    return m.PQ2_branch_ij[i,j,k] == P_branch_ij**2 + Q_branch_ij**2
def thermal_lim_ineq_con(m, i, j, k):
    return sqrt(m.PQ2_branch_ij) <= limits_flows[i, j]
if line_loading_limit:
    model.compute_PQ2_branch_ij = Constraint(branch_from_to, scenario_set, rule=compute_PQ2_branch_ij)
    model.thermal_lim_ineq_con = Constraint(branch_from_to, scenario_set, rule=thermal_lim_ineq_con)`
  1. Remove the sqrt() from the formulation by squaring limits_flows[i,j]. This is usually the best option: your problem is mathematically equivalent, but no longer contains the sqrt() function and will not have the problem if the branch flow goes to 0. The infeasibility that you observed when you tried this approach may be a separate issue entirely. For example, with nonconvex problems (like the ACOPF) you can initialize the solver in such a way that it cannot move to a feasible position and will "converge to a point of local infeasibility" (for more explanation, see your specific solver's documentation; e.g. IPOPT, Knitro).
def thermal_lim_ineq_con(model, i, j, k):
    if line_loading_limit == True:
        P_branch_ij = (model.V_magn[i, k]**2) * np.float(real(Y_matrix[i, j])) - model.V_magn[i, k] * model.V_magn[j, k] * (np.float(real(Y_matrix[i, j])) * cos(model.V_angle[i, k]-model.V_angle[j, k]) + np.float(imag(Y_matrix[i,j])) * sin(model.V_angle[i, k]-model.V_angle[j, k]))
        Q_branch_ij = (model.V_magn[i, k]**2) * np.float(imag(Y_matrix[i, j])) + model.V_magn[i, k] * model.V_magn[j, k] * (np.float((real(Y_matrix[i, j]))) * sin(model.V_angle[i, k]-model.V_angle[j, k]) - np.float(imag(Y_matrix[i,j])) * cos(model.V_angle[i, k]-model.V_angle[j, k]))
        return P_branch_ij**2 + Q_branch_ij**2 <= limits_flows[i, j]**2
    else: 
        return Constraint.Skip

Upvotes: 3

Related Questions