Solving a Squared-root constraint in Pyomo with Gurobi

I'm using Gurobi 9.1.2 and Pyomo 6.1.2

I got a MILP model in pyomo created with pyomo.environ modeling layer.

In order to advance in my research, I need to implement a new constraint to my model. This new constraint got a square-root and I want to solve it with Gurobi.

import pyomo.environ as pyo
...
model = pyo.AbstractModel()
#Set
model.J = pyo.Set(doc='Generation set')
model.t = pyo.Set(doc='datetime Set')
model.S = pyo.Set(doc='Storage Systems Set')
#Params
model.Delt = pyo.Param(model.t)
model.Us = pyo.Param()
model.Ud = pyo.Param()
model.Ss = pyo.Param()
model.Sd = pyo.Param()
model.erf = pyo.Param()
#Variables
model.pnomj = pyo.Var(model.J, domain=pyo.NonNegativeReals)
model.pjtResUp = pyo.Var(model.J, model.T, domain=pyo.NonNegativeReals)
model.pbesstResUp = pyo.Var(model.S, model.T, domain=pyo.NonNegativeReals)

The constraint I want to implement is this one Quadratic Constraint

where Δp^{SR+}_{j,t} is model.pjtResUp and Δp^{SR+}__{B,t} is model.pbesstResUp

def upReserves(model, t):
    return model.Us*model.pnomj['PV'] - model.Ud*model.Delt[t] + \
           model.erf*pyo.sqrt((model.Ss*model.pnomj['PV'])**2 + (model.Sd*model.Delt[t])**2) <= \
           sum(model.pjtResUp[j,t] for j in model.J) + model.pbesstResUp['BESS',t]

This gives me the next ValueError:

ValueError: Cannot write legal LP file.  Constraint 'upReserves[2020-01-01 00:00:00]' has a body with nonlinear terms.

In This Answer they say there's a Pyomo problem when creating the LP file with CPLEX, since Gurobi can solve even NonConvex problem (Passing the argument options={'NonConvex':2} to SolverFactory)

Is there some recommendation which I can check? Thanks in advances

I get to solved this one.

In order to avoid non-convexity in the conic-constraint-style constraint, this is what I did:

1. To avoid Non-Convexity in the SOC constraint, I used a new var (let's call it g) such that:

    model.g = pyo.Var(model.T, doc='non-squared part of conic constraint')
   

2. Add a new constraint of the non-squared-root part of the conic constraint:

    def RHS_constraint(model, t):
     return model.g[t] == sum(model.pjtResUp[j,t] for j in model.J) + model.pbesstResUp['BESS',t] + (model.Ud*model.Delt[t] - model.Us*model.pnomj['PV'])
   

   model.RHS_constraint = pyo.Constraint(model.T, rule=RHS_constraint)

3. To avoid the ValueError in non-linear terms by Pyomo LP file creation, I squared both sides:

    def upReserves(model, t):
        return model.g[t]**2 >= (model.erf**2)*((model.Ss*model.pnomj['PV'])**2 + (model.Sd*model.Delt[t])**2)
    model.upReserves= pyo.Constraint(model.T, rule=upReserves)
   

Thanks for those who at least saw this post. If you think this approach is wrong, I'm also seeking for new approach