problem with warnings
AnsweredHello everyone,
I haven't much experience with gurobi. I am trying to build this model but i take this warnings. I searched and i found somethings but i don't understand what i must do. I have also problem with voltage limit. I want to set for as a constraint with low and upper bound but from my research i saw that is possible only if i put 2 constraints with different name. I did it but doesn't work. For now i just put only one to run the code and i take this warnigs.
Can any one help me please?
Warning: zero or small (< 1e-13) coefficients in quadratic constraints, ignored
Initializing OPF model...
Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3096 rows, 8428 columns and 41280 nonzeros
Model fingerprint: 0x2434bfde
Model has 7568 quadratic constraints
Coefficient statistics:
Matrix range [4e-02, 4e+07]
QMatrix range [1e+00, 4e+15]
QLMatrix range [1e+00, 2e+03]
Objective range [1e+00, 1e+03]
Bounds range [3e-01, 1e+00]
RHS range [3e+03, 6e+03]
QRHS range [1e+02, 8e+09]
Warning: Quadratic constraints contain large coefficient range
Warning: Model contains large rhs on quadratic constraints
Consider reformulating model or setting NumericFocus parameter
to avoid numerical issues.
Presolve time: 0.01s
Barrier solved model in 0 iterations and 0.01 seconds (0.00 work units)
Model is infeasible
Solving OPF model...
Extracting data
Code:
import gurobipy as grb
def OPF_model_creator(pd,pyo,math,Time_sim,V_init_pu,Load_data,PV_set,PV_Gen_data,Bus_Vnom,V_statutory_lim,min_Cosphi,Inverter_S_oversized,Cable_data,Transformer_rating,security_margin_current,security_margin_Transformer_S):
model = grb.Model() # Defines an abstract model which data can be imported
## Core sets definition ##
# Loaded from csv in OPF_model\Pyomo_OPF_dat
b = pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Buses_Pyomo.csv").values
ph = pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Phases_Pyomo.csv").values
ln = pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Lines_Pyomo.csv").values
lds = pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Loads_Pyomo.csv").values
gn = pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Generator_Pyomo.csv").values
ph_gen = pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Generator_phases.csv").values
Bus_g_b = pd.read_csv("C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Y_bus_Pyomo.csv").values
conn = pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Connectivity_Pyomo.csv").values
Lines_data= pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Lines_data_Pyomo.csv").values
Loads_data= pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Loads_data_Pyomo.csv").values
Gen_data= pd.read_csv(
"C:/Users/titos/PycharmProjects/pythonProject11/Open-DSOPF-master/OPF_model/Pyomo_OPF_data/Generator_data_Pyomo.csv").values
def read_data_conn(g):
connection = {}
for i in range(len(g[:, 0])):
connection[str(g[i, 0]), str(g[i, 1])] = int(g[i, 2])
return connection
def loads_dt(g):
Load_bus_conn = {}
Load_phase_conn = {}
Load_Vnom={}
for i in range(len(g[:, 0])):
Load_bus_conn[str(g[i, 0])] = str(g[i, 1])
Load_phase_conn[str(g[i, 0])] = int(g[i, 2])
Load_Vnom[str(g[i, 0])] = int(g[i, 3])
return [Load_bus_conn,Load_phase_conn,Load_Vnom]
def Generator_dt(g):
Generator_bus_conn = {}
Generator_phase_conn = {}
Generator_Vnom={}
for i in range(len(g[:, 0])):
Generator_bus_conn[str(g[i, 0])] = str(g[i, 1])
Generator_phase_conn[str(g[i, 0])] = int(g[i, 2])
Generator_Vnom[str(g[i, 0])] = int(g[i, 3])
return [Generator_bus_conn,Generator_phase_conn,Generator_Vnom]
def lines_dt(g):
Lines_k = {}
Lines_i = {}
Lines_cable={}
for i in range(len(g[:, 0])):
Lines_k[str(g[i, 0])] = str(g[i, 1])
Lines_i[str(g[i, 0])] = str(g[i, 2])
Lines_cable[str(g[i, 0])] = str(g[i, 3])
return [Lines_k,Lines_i,Lines_cable]
def read_data_2(g):
bus_G = {}
bus_B = {}
for i in range(len(g[:, 0])):
bus_G[str(g[i, 0]),str(g[i, 1]), int(g[i, 2]), int(g[i, 3])] = float(g[i, 4])
bus_B[str(g[i, 0]),str(g[i, 1]), int(g[i, 2]), int(g[i, 3])] = float(g[i, 5])
return [bus_G,bus_B]
def read(dt,type):
return {type(i[0]) for i in dt}
Buses = (read(b,str))
Phases_abc = (read(ph,int))
Lines = (read(ln,str))
Loads = (read(lds,str))
PVs = (read(gn,str))
Phases = (read(ph_gen,int))
# Defined h
time = {int(i) for i in Time_sim}
for i in Phases_abc:
print(type(i))
# PVs = {i for i in PV_set}
## Parameters ##
if 1 == 1:
## Loaded from csv in OPF_model\Pyomo_OPF_data\
# Network parameters
[Bus_G,Bus_B] =read_data_2(Bus_g_b)
Connectivity = read_data_conn(conn)
[Lines_k,Lines_i,Lines_cable]=lines_dt(Lines_data)
# Load characteristics
[Load_bus_conn,Load_phase_conn,Load_Vnom]=loads_dt(Loads_data)
# pvgen characteristics
[Generator_bus_conn, Generator_phase_conn, Generator_Vnom]=Generator_dt(Gen_data)
# Variables ##
if 1 == 1:
# PV Control variables
if 1 == 1:
P_control = model.addVars(PVs,time,vtype='C',lb=0.0,ub= 1.0,name='P_control')
Thanphi_control = model.addVars(PVs,time,vtype='C',lb=0.0,ub=math.tan(math.acos(min_Cosphi)))
# Voltages - Real and Imaginary part - in V
if 1 == 1:
# Initialization rules
def init_Vre_rule(k, s, t):
bus_check = V_init_pu['Bus_k'] == k
phase_check = V_init_pu['Phase_k'] == str(s)
return float(V_init_pu.loc[bus_check & phase_check, 'Vinitre_pu'].values[0]) * float(
V_init_pu.loc[bus_check & phase_check, 'Vinit'].values[0])
def init_Vim_rule(k, s, t):
bus_check = V_init_pu['Bus_k'] == k
phase_check = V_init_pu['Phase_k'] == str(s)
return float(V_init_pu.loc[bus_check & phase_check, 'Vinitim_pu'].values[0]) * float(
V_init_pu.loc[bus_check & phase_check, 'Vinit'].values[0])
# Variables definition
Vre = model.addVars(Buses, Phases_abc,time,vtype='C') # Real and imaginary phase voltages
Vim = model.addVars(Buses, Phases_abc,time,vtype='C')
Isp_re = model.addVars(Buses, Phases_abc, time, vtype='C')
Isp_im = model.addVars(Buses, Phases_abc, time, vtype='C')
P_flow_sending= model.addVars(Lines, Phases_abc, time, vtype='C')
Q_flow_sending= model.addVars(Lines, Phases_abc, time, vtype='C')
P_flow_receiving = model.addVars(Lines, Phases_abc, time, vtype='C')
Q_flow_receiving= model.addVars(Lines, Phases_abc, time, vtype='C')
## Expressions ##
for k in Buses:
for s in Phases_abc:
for t in time:
#initialize
Vre[k,s,t].Start=init_Vre_rule(k, s, t)
Vim[k,s,t].Start=init_Vim_rule(k, s, t)
Isp_re[k,s,t].Start=0.0
Isp_im[k,s,t].Start=0.0
#initialize p_control && tan
for pv in PVs:
for t in time:
P_control[pv,t].Start=0.001
Thanphi_control[pv,t].Start=0.001
## Expressions ##
if 1 == 1:
# Voltage magnitude
def V_mag_2(k, s, t):
return Vre[k, s, t]**2+Vim[k, s, t]**2
# PV P and Q
if 1 == 1:
def PV_P(pv, t):
return PV_Gen_data[pv, 'Profile_P'].loc[t] * P_control[pv, t]
def PV_Q(pv, t):
return PV_P(pv,t) * -Thanphi_control[pv, t]
# Power injection per node
if 1 == 1:
# Active Power demand load
def PDem(k, s, t):
return sum(Load_data[h, 'P_Profile'].loc[t] for h in Loads if(Load_bus_conn[h] == k and Load_phase_conn[h] == s))
# Reactive Power demand load
def QDem(k, s, t):
return sum(Load_data[h, 'Q_Profile'].loc[t] for h in Loads if(Load_bus_conn[h] == k and Load_phase_conn[h] == s))
# Active power generation
def PGen(k, s, t):
return sum(PV_P(pv, t) for pv in PVs if (Generator_bus_conn[pv] == k and Generator_phase_conn[pv] == s))
def QGen(k, s, t):
return sum(PV_Q(pv, t) for pv in PVs if(Generator_bus_conn[pv] == k and Generator_phase_conn[pv] == s))
# Power injections per node (generation - demand)
def Psp(k, s, t):
return PGen(k, s, t) - PDem(k, s, t)
def Qsp(k, s, t):
return QGen(k, s, t) - QDem(k, s, t)
# Current injections
if 1 == 1:
# Specified current injection (from specified PQ injections)
if 1 == 1:
for k in Buses:
for s in Phases_abc:
for t in time:
model.addConstr(
Psp(k, s, t) == Vre[k, s, t] * Isp_re[k, s, t] + Vim[k, s, t] * Isp_im[k, s, t])
model.addConstr(
Qsp(k, s, t) == Vim[k, s, t] * Isp_re[k, s, t] - Vre[k, s, t] * Isp_im[k, s, t])
# Calculated current injection (from voltages and admittances)
if 1 == 1:
# Expressions for currents
def Icalc_re(k, s, t): # Real calculated current injection
return sum(sum(
Bus_G[k, i, s, a] * Vre[i, a, t] - Bus_B[k, i, s, a] * Vim[
i, a, t] for a in Phases_abc) for i in Buses if
Connectivity[k, i] == 1)
def Icalc_im(k, s, t): # Imaginary calculated current injection
return sum(sum(
Bus_G[k, i, s, a] * Vim[i, a, t] + Bus_B[k, i, s, a] * Vre[
i, a, t] for a in Phases_abc) for i in Buses if
Connectivity[k, i] == 1)
# Current line flows
if 1 == 1:
def Iflow_re(l, s, t): # Real line flows
return sum(Bus_G[Lines_k[l], Lines_i[l], s, a] * (
Vre[Lines_i[l], a, t] - Vre[Lines_k[l], a, t]) -
Bus_B[Lines_k[l], Lines_i[l], s, a] * (
Vim[Lines_i[l], a, t] - Vim[Lines_k[l], a, t])
for a in Phases_abc)
def Iflow_im(l, s, t): # Imaginary line flows
return sum(Bus_B[Lines_k[l], Lines_i[l], s, a] * (
Vre[Lines_i[l], a, t] - Vre[Lines_k[l], a, t]) +
Bus_G[Lines_k[l], Lines_i[l], s, a] * (
Vim[Lines_i[l], a, t] - Vim[Lines_k[l], a, t])
for a in Phases_abc)
# PQ line flows
if 1 == 1:
for l in Lines:
for s in Phases_abc:
for t in time:
model.addConstr(P_flow_sending[l, s, t] == Vre[Lines_k[l], s, t] * Iflow_re(l, s, t) + Vim[Lines_k[l], s, t] * Iflow_im(l, s, t))
model.addConstr(P_flow_receiving[l, s, t] == -1.0 * (Vre[Lines_i[l], s, t] * Iflow_re(l, s, t) + Vim[
Lines_i[l], s, t] * Iflow_im(l, s, t)))
model.addConstr(Q_flow_sending[l, s, t] == Vim[Lines_k[l], s, t] * Iflow_re(l, s, t) - Vre[
Lines_k[l], s, t] * Iflow_im(l, s, t))
model.addConstr(Q_flow_receiving[l, s, t] ==-1.0 * (Vim[Lines_i[l], s, t] * Iflow_re(l, s, t) - Vre[
Lines_i[l], s, t] * Iflow_im(l, s, t)))
# Losses
# Losses fotime step per line and phase
def P_losses(l, s, t): # Active power losses Watts!
return P_flow_sending[l, s, t] + P_flow_receiving[l, s, t]
# Total losses fotime step
def P_losses_total(s, t): # Total active power losses Watts!
return sum(P_losses(l, s, t) for l in Lines)
## Constraints ##
# Power flow constraints
if 1 == 1:
for k in Buses:
for s in Phases_abc:
for t in time:
model.addConstr(Isp_re[k, s, t] - Icalc_re(k, s, t) == 0.0)
model.addConstr(Isp_im[k, s, t] - Icalc_im(k, s, t) == 0.0)
if k == 'slack':
bus_check = V_init_pu['Bus_k'] == k
phase_check = V_init_pu['Phase_k'] == str(s)
model.addConstr(Vim[k, s, t] == float(
V_init_pu.loc[bus_check & phase_check, 'Vinitim_pu'].values[0]) * float(
V_init_pu.loc[bus_check & phase_check, 'Vinit'].values[0]))
model.addConstr(Vre[k, s, t] == float(
V_init_pu.loc[bus_check & phase_check, 'Vinitre_pu'].values[0]) * float(
V_init_pu.loc[bus_check & phase_check, 'Vinit'].values[0]))
# Network operational constraints
# Votlage limits
if 1 == 1:
for k in Buses:
for s in Phases_abc:
for t in time:
if k == 'slack':
continue
else:
V_nom = Bus_Vnom.loc[k, 'Vnom_pn']
# model.addConstr((V_nom * V_statutory_lim[0])** 2 <= V_mag_2(k, s, t),name="left_hand_side")
model.addConstr(V_mag_2(k, s, t)<=(V_nom * V_statutory_lim[1])**2,name="right_hand_side")
if 1 == 1:
for l in Lines:
for s in Phases_abc:
for t in time:
model.addConstr(Iflow_re(l, s, t) ** 2 + Iflow_im(l, s, t) ** 2 <= (
security_margin_current * Cable_data.loc[
Lines_cable[l], 'Current rating [A]']) ** 2)
# Transformer power limit
if 1 == 1:
for t in time:
lhs=(sum(sum((P_flow_sending[l, s, t]**2 + Q_flow_sending[l, s, t]**2) for s in Phases_abc)for l in Lines if Lines_k[l] == 'secondary'))
rhs=(security_margin_Transformer_S *1000* Transformer_rating ) ** 2
model.addConstr(lhs <= rhs)
# PV operational limits
if 1 == 1:
for pv in PVs:
for t in time:
max_P = PV_Gen_data[pv, 'Profile_P'].loc[t]
max_S = (max_P ** 2 + (max_P * math.tan(math.acos(min_Cosphi))) ** 2) ** 0.5
if max_S > PV_Gen_data[pv, 'Rating'] * 1000 * Inverter_S_oversized:
lhs = PV_P(pv, t) ** 2 + PV_Q(pv, t) ** 2
rhs = (PV_Gen_data[pv, 'Rating'] * 1000 * Inverter_S_oversized) ** 2
model.addConstr(lhs <= rhs)
""" Objective function """
if 1 == 1:
def Objective_func(model):
return sum(
sum((-1e3 *P_control[pv, t] + Thanphi_control[pv, t]) for pv in PVs) for t in time)
model.setObjective(sense=grb.GRB.MINIMIZE, expr=Objective_func(model))
return model # the problem is solved at the main script
-
Please note that it is not possible to execute your code, because it uses multiple excel sheets. Moreover, the code is very long and it is not clear where the source of your issue lies. It is always best to generate a minimal reproducible example to get the best help possible.
The warning
Warning: zero or small (< 1e-13) coefficients in quadratic constraints, ignored
means that you have coefficient with value lesser than \(10^{-13}\) in your model and that Gurobi is going to treat all such coefficient as \(0\). In general, it is strongly recommended that all meaningful coefficients are at least one order of magnitude greater than the FeasibilityTol value.
The warnings
Coefficient statistics:
Matrix range [4e-02, 4e+07]
QMatrix range [1e+00, 4e+15]
QLMatrix range [1e+00, 2e+03]
Objective range [1e+00, 1e+03]
Bounds range [3e-01, 1e+00]
RHS range [3e+03, 6e+03]
QRHS range [1e+02, 8e+09]
Warning: Quadratic constraints contain large coefficient range
Warning: Model contains large rhs on quadratic constraints
Consider reformulating model or setting NumericFocus parameter
to avoid numerical issues.mean that the range of coefficients in your linear and quadratic constraints is very large. Indeed, the ranges have differences of \(9\) and \(15\) orders of magnitude, which very likely leads to numerical issues within the optimization algorithm. You can try experimenting with the NumericFocus parameter to overcome these problems but the best way to deal with bad scaling and numerical issues is to re-scale the coefficients of your model. Please have a look at our Guidelines for Numerical Issues for more information.
I want to set for as a constraint with low and upper bound but from my research i saw that is possible only if i put 2 constraints with different name. I did it but doesn't work.
I don't fully understand your issue here. You can add 2 constraints with different names or just don't provide any names at all, e.g.,
model.addConstr(x <= 1, name="x_ub")
model.addConstr(x >= 0, name="x_lb")Alternatively, you can use the addRange method. Note that addRange accepts linear terms only. For quadratic terms, you have to add two constraints manually.
Best regards,
Jaromił0
Please sign in to leave a comment.
Comments
1 comment