Linearization of quadratic constraint

I encountered a quadratic constraint in my pyomo model. It worked more or less solving it with gurobi but it often gave me memory issues. So I linearized this quadratic constraint. However, now I run into another problem: Link to Stackoverflow. So I am wondering if I made a mistake in the linearization.

The problem is a heat pump that can deliver cooling or heating but not both at the same time. Quadratic version:

h(t) = p(t)*bh(t)*COPh  #quadratic const
c(t) = p(t)*bc(t)*COPc  #quadratic const
h(t) <= cap_hp
c(t) <= cap_hp*(COPc/(COPc+1)
bh(t) + bc(t) <= 1

Linearised version:

p(t) = h(t)/COPh + c(t)/COPc
h(t) <= cap_hp
c(t) <= cap_hp*(COPc/(COPc+1)
h(t) <= M * b(t) + 0 * (1-b(t))
c(t) <= 0 * b(t) + M * (1-b(t))

Decision variables: p(t): electricity consumption of heating or cooling; h(t): heating output; c(t): cooling output; cap_hp: maximal heating capacity (size of heat pump); b(t): binary variable that is 1 during heating or 0 during cooling. Or in quadratic version bh(t) or bc(t) binary variables for heating and cooling respectively. Input parameters: COPh/c: "efficiency of heatpump in heating or cooling mode; M: large number

Below pyomo code (ignore technicalities of a heat pump eg COPs, those are changed by now). Quadratic:

m.b_hph = Var(year_i, ts_i, within=Binary)
m.b_hpc = Var(year_i, ts_i, within=Binary)
   
def ashpHeat_rule(m,y,ts):
    return sum(m.heat["heat_pump_air", hCons, y, ts] for hCons in hIn) == \
        sum(m.power[elSup, "heat_pump_air",y, ts] for elSup in elOut)\
        * m.b_hph[y,ts] * 3
m.const_ashpHeat = Constraint(year_i, ts_i, rule = ashpHeat_rule)

def ashpCool_rule(m,y,ts):
    return sum(m.cool["heat_pump_air", cCons, y, ts] for cCons in cIn) == \
        sum(m.power[elSup, "heat_pump_air",y, ts] for elSup in elOut)\
        * m.b_hpc[y,ts] * 2
m.const_ashpCool = Constraint(year_i, ts_i, rule = ashpCool_rule)

def ashpCapah_rule(m,y, ts):
    return sum(m.heat["heat_pump_air", hCons, y, ts] for hCons in hIn) <=\
        m.c_c["heat_pump_air"]
m.const_ashphCapa = Constraint(year_i, ts_i, rule = ashpCapah_rule)

def ashpCapac_rule(m,y, ts):
    return sum(m.cool["heat_pump_air", cCons, y, ts] for cCons in cIn) <=\
        m.c_c["heat_pump_air"] * (2/3)
m.const_ashpcCapa = Constraint(year_i, ts_i, rule = ashpCapac_rule)

def hpbin_rule(m,y,ts):
    return m.b_hph[y,ts] + m.b_hpc[y,ts] <= 1
m.const_hpbin = Constraint(year_i, ts_i, rule = hpbin_rule)

Linear:

m.b_hph = Var(year_i, ts_i, within=Binary)

def ashpcons_rule(m,y,ts):
    return sum(m.power[elSup, "heat_pump_air",y, ts] for elSup in elOut) == \
        sum(m.heat["heat_pump_air", hCons, y, ts] for hCons in hIn)/4.5 +\
        sum(m.cool["heat_pump_air", cCons, y, ts] for cCons in cIn)/3.5
m.const_ashpcons = Constraint(year_i,ts_i, rule = ashpcons_rule)

def ashpheatdecision_rule(m,y,ts):
    return sum(m.heat["heat_pump_air", hCons, y, ts] for hCons in hIn) <=\
        99999 * m.b_hph[y,ts] + 0 * (1-m.b_hph[y,ts])
m.const_ashpheatdecision = Constraint(year_i, ts_i, rule = ashpheatdecision_rule)

def ashpcooldecision_rule(m,y,ts):
    return sum(m.cool["heat_pump_air", cCons, y, ts] for cCons in cIn) <=\
        0 * m.b_hph[y,ts] + 99999 * (1-m.b_hph[y,ts])
m.const_ashpcooldecision = Constraint(year_i, ts_i, rule = ashpcooldecision_rule)

def ashpCapah_rule(m,y, ts):
    return sum(m.heat["heat_pump_air", hCons, y, ts] for hCons in hIn) <=\
        m.c_c["heat_pump_air"]
m.const_ashphCapa = Constraint(year_i, ts_i, rule = ashpCapah_rule)

def ashpCapac_rule(m,y, ts):
    return sum(m.cool["heat_pump_air", cCons, y, ts] for cCons in cIn) <=\
        m.c_c["heat_pump_air"] * (3.5/4.5)
m.const_ashpcCapa = Constraint(year_i, ts_i, rule = ashpCapac_rule)

Thanks a lot for your help

Axel