Pyomo optimization RuntimeError

Question

I'm trying to implement an optimization problem on gurobi unsuccessfully. I have this script below in which I cut out the upper part because it would be huge.

\##CONCRETE MODEL
ONR = ConcreteModel(name="Radialità+Connettività")

\#COSTANTI

\#VARIABILI CONTINUE
\#Ci dovrei inserrire Fij e Fsorgente ma alla fine dovrei inserirle come variabili binarie
ONR.f_szk = Var(NomiBus, within=Reals) #per vedere se funziona nella cmd eseguire "ONR.f_s.pprint()"  
ONR.fijk = Var(NomiBus, nomi_switch, within=Reals)
ONR.lamb_k = Var(NomiBus, within=Reals) #NonNegativeReals(?)
ONR.ENS = Var(within=Reals) #e che set di eq devo dargli?
ONR.U_k = Var(NomiBus, within=Reals)

\#VARIABILI BINARIE
ONR.y_ij = Var(nomi_switch, within=Binary)
ONR.Z_act_i = Var(zona_bus, NomiBus, within=Binary)
ONR.Z_act_j = Var(zona_bus, NomiBus, within=Binary)
ONR.w_zk_i = Var(zona_bus, NomiBus, within=Binary)
ONR.w_zk_j = Var(zona_bus, NomiBus, within=Binary)
ONR.Z_act = Var(zona_bus, NomiBus, within=Binary)

\#Definisco le funzioni costo in questo caso della lunghezza della linea solo per verificare il corretto funzionamento del codice
Ca = 10
c_ENS = 10
c_SAIDI = 10
c_SAIFI = 10

\#VINCOLI

\#26
\#27
def FlussiArtificiali_rule(ONR,bus,zona):

    membroA = 0
    membroB = 0
    
    for sw in nomi_switch:
        
        i = dict_bus_from_sw[sw]
        j = dict_bus_to_sw[sw]
        
        # se il bus k è nella zona z, allora fdzk = 1:
        if bus in dict_bus_zone[zona]:
            membroB = 1
        
        # se il busfrom dello switch sta nella zona z prendi il flusso col segno meno:
        if i in dict_bus_zone[zona]:
            membroA = membroA - ONR.fijk[bus,sw]
        # se il busto dello switch sta nella zona z prendi il flusso col segno più:
        elif j in dict_bus_zone[zona]:
            membroA = membroA + ONR.fijk[bus,sw] 
           
    # Solo per la zona contenente la slack, aggiungi fskz:           
    if zona == 'Z0':
        membroA = membroA + ONR.f_szk[bus]
    
    # Equazione di bilancio dei flussi artificiali:        
    return membroA == membroB

ONR.FlussiArtificiali = Constraint(NomiBus,zona_bus,rule=FlussiArtificiali_rule)

\#Rubata da Roberto per garantire la radialità
def LimitY_rule(ONR,y_ij):
return sum(ONR.y_ij\[sw\] for sw in nomi_switch) \<= 7
ONR.LimitY = Constraint(nomi_switch,rule=LimitY_rule)

\#28
def stato_flusso_ij_switch_UB_rule(ONR, switch, bus):
return ONR.fijk\[bus,switch\] \<= ONR.y_ij\[switch\]
ONR.stato_flusso_ij_switch_UB = Constraint(nomi_switch, NomiBus, rule=stato_flusso_ij_switch_UB_rule)

def stato_flusso_ij_switch_LB_rule(ONR, switch, bus):
return ONR.fijk\[bus,switch\] \>= - ONR.y_ij\[switch\]
ONR.stato_flusso_ij_switch_LB = Constraint(nomi_switch, NomiBus, rule=stato_flusso_ij_switch_LB_rule)

\#29
def Zact_rule(ONR,sw,bus):
i = dict_bus_from_swfrom\[sw\]
j = dict_bus_from_sw\[sw\]
if i in dict_bus_zone\[zona\] or j in dict_bus_zone\[zona\]:
return (ONR.Z_act\[zona,bus\] \>= ONR.fijk\[bus, sw\]) or (- ONR.Z_act\[z,k\] \<= ONR.fijk\[bus, sw\])
ONR.Z_act = Constraint(nomi_switch,NomiBus,rule=Zact_rule)  \

So I'm trying to implement Z_act_rule in which I should know if my zone is active for the shortest path. I look for the fij in both i and j nodes and then I want to know if one of them is active. But I do get the following error

I even tried to do it like this

def zona_attiva_j_UB_rule(ONR, switch, bus):
    for zona in dict_bus_zone.keys(): 
        return ONR.Z_act_j[zona, bus] <= - ONR.fijk[bus, switch] 
ONR.zona_attiva_j_UB = Constraint(nomi_switch, NomiBus, rule=zona_attiva_j_UB_rule)

But it does not loop for every zone but only for one. Instead I would like to run it for every zone and to know which one is active and part of the shortest path