I implemented an ACOPF problem to maximise the PV exports of a LV network, to validate the answers and also to compare the time it takes to solve each I implemented the same problem in pyomo, while the powermodels.jl solves the problem rather quickly, pyomo throws an infeasible error. This is the pyomo model:
model = ConcreteModel()
model.N = Set(initialize = AllNodes)
model.L = Set(within = model.N*model.N, initialize=modified_AllLines)
model.Slackbus = Set(initialize = [AllNodes[0],AllNodes[1],AllNodes[2]])
model.T = Set(initialize=Time_sim)
model.PV_Nodes = Set(initialize=Load_nodes)
slack_bus_dict = {node: (1 if node in model.Slackbus else 0) for node in model.N}
model.IsSlackBus = Param(model.N, initialize=slack_bus_dict, within=Binary)
model.SlackBusVoltageMagnitude = Param(initialize = 22000**2)
model.R = Param(model.L,initialize = filtered_R)
model.X = Param(model.L,initialize = filtered_X)
model.S_max = Param(model.L,initialize = S_max)
model.p_d = Param(model.N,model.T,initialize=P_d)
model.q_d = Param(model.N,model.T,initialize=Q_d)
Vmin=216**2
Vmax=253**2
model.p = Var(model.L,model.T,within=Reals) # Active power flow (kW)
model.q = Var(model.L,model.T, within=Reals) # Reactive power flow (kVAr)
model.I = Var(model.L,model.T, within=Reals) # Current squared (A^2)
model.v = Var(model.N,model.T, within=Reals, bounds =(Vmin,Vmax)) # Voltage squared (V^2)
model.p_exp = Var(model.PV_Nodes,model.T,within=NonNegativeReals,bounds=(0,10000))
R and X is modelled for each node, for example (nodei,nodej):rij, p_d = pgen-pexport(actual pv generation). Objective is to find the maximum PV that can be exported(eventhough it is actually generated or not)