Problem with reciprocal constraints in Gurobi

How to implement P0 formula from the below link in Gurobi Python?

https://people.revoledu.com/kardi/tutorial/Queuing/MMs-Queuing-System.html

import gurobipy as gp
from gurobipy import GRB
import math

# Create a new model
m = gp.Model()

# Create variables
γ = m.addVar(vtype=GRB.CONTINUOUS, name="Arrival_rate")
μ = m.addVar(vtype=GRB.CONTINUOUS, name="Service_rate")
W = m.addVar(vtype=GRB.CONTINUOUS, name="Avg_time_a_truck_waits_in_the_system")

# Define an auxiliary variables to handle the constraints
ρ = m.addVar(vtype=GRB.CONTINUOUS, name="Ratio_of_traffic_intensity")
U = m.addVar(lb=0.0, ub=1.0, vtype=GRB.CONTINUOUS, name="Utilization_factor")
S_t = m.addVar(vtype=GRB.CONTINUOUS, name="Service_time")
L = m.addVar(vtype=GRB.INTEGER, name="Avg_no_of_trucks_in_the_terminal")
L_q = m.addVar(vtype=GRB.INTEGER, name="Avg_trucks_in_line_for_service")
W_q = m.addVar(vtype=GRB.CONTINUOUS, name="Avg_trucks_waiting_in_line_for_service")
P_0 = m.addVar(lb=0.0, ub=1.0, vtype=GRB.CONTINUOUS, name="Probability_of_trucks_in_the_terminal")
P_n = m.addVar(lb=0.0, ub=1.0, vtype=GRB.CONTINUOUS, name="Probability_of_n_trucks_in_the_terminal")
a = m.addVar(vtype=GRB.CONTINUOUS, name='a')

# Define possible values for δ and ρ
delta_values = [1, 2, 3, 4]
rho_values = {
    2: [1.5, 1.6, 1.7, 1.8, 1.9],
    3: [1.5, 1.6, 1.7, 1.8, 1.9, 2.25, 2.4, 2.55, 2.7, 2.85],
    4: [1.5, 1.6, 1.7, 1.8, 1.9, 2.25, 2.4, 2.55, 2.7, 2.85, 3, 3.2, 3.4, 3.6, 3.8]
}

# Define binary variables for δ
delta_bin_vars = m.addVars(delta_values, vtype=GRB.BINARY, name="delta_bin")

# Define binary variables for ρ for each δ
rho_bin_vars = {}
for d in delta_values:
    if d in rho_values:
        rho_bin_vars[d] = m.addVars(len(rho_values[d]), vtype=GRB.BINARY, name=f"rho_bin_{d}")

# Constraint to ensure exactly one δ is selected
m.addConstr(gp.quicksum(delta_bin_vars[d] for d in delta_values) == 1, "one_delta")

# Constraints for ρ based on selected δ
for d in delta_values:
    if d in rho_values:
        # Constraint to ensure exactly one ρ is selected for the given δ
        m.addConstr(gp.quicksum(rho_bin_vars[d][i] for i in range(len(rho_values[d]))) == delta_bin_vars[d], f"one_rho_{d}")
        # Define ρ as a weighted sum of the selected values
        m.addConstr(ρ == gp.quicksum(rho_values[d][i] * rho_bin_vars[d][i] for i in range(len(rho_values[d]))), f"rho_value_{d}")

# Add constraints
m.addConstr(ρ * μ == γ, "Traffic_intensity")
m.addConstr(U * gp.quicksum(d * delta_bin_vars[d] for d in delta_values) == ρ, "Utilization_factor")
m.addConstr(S_t * μ == 1, "Service_time")
m.addConstr(L == γ * W, "Trucks_in_terminal")
m.addConstr(W_q == L_q * γ, "Trucks_in_waiting_line")

# Constraints for P_0
# Define the expression inside the summation
sum_expression = gp.LinExpr()

for i in range(1, max(delta_values)):
    # Calculate ρ^i
    term_expression = gp.LinExpr()
    term_expression.add(ρ, i)  # Equivalent to ρ^i

    # Calculate (ρ^i) / i!
    coefficient = 1 / math.factorial(i)
    term_expression *= coefficient

    # Add term to sum_expression
    sum_expression.add(term_expression)

for d in delta_values:
    # Define the expression involving δ, μ, and γ
    m.addGenConstrPow(ρ, a, d)

    # Create an auxiliary variable for the fraction part (d * μ) / (d * μ - γ)
    aux_var_frac_part = m.addVar(vtype=GRB.CONTINUOUS, name=f"aux_var_frac_part_{d}")
    
    # Define the coefficient for the constraint
    coefficient = a / math.factorial(d)
    
    # Calculate (d * μ)
    numerator = d * μ
    
    # Calculate (d * μ - γ)
    denominator = numerator - γ
    
    # Introduce an auxiliary variable to represent 1 / denominator
    inv_denominator = m.addVar(vtype=GRB.CONTINUOUS, name=f"inv_denominator_{d}")
    
    # Add a constraint: inv_denominator * denominator = 1
    m.addConstr(inv_denominator * denominator == 1, f"inv_denominator_def_{d}")

    # Now, aux_var_frac_part = numerator * inv_denominator
    m.addConstr(aux_var_frac_part == numerator * inv_denominator, f"aux_var_frac_part_def_{d}")
    
    # Define the term expression
    term_expression = coefficient * aux_var_frac_part
    
    # Create an auxiliary variable for the final term
    aux_var_frac = m.addVar(vtype=GRB.CONTINUOUS, name=f"aux_var_frac_{d}")
    
    # Add the constraint defining the auxiliary variable for the final term
    m.addConstr(aux_var_frac == term_expression, f"aux_var_frac_def_{d}")
    
    # Add to frac_expression
    frac_expression += delta_bin_vars[d] * aux_var_frac

# Define the entire right-hand side expression raised to the power of -1
rhs_expression = sum_expression + frac_expression

rhs_numerator = m.addVar(vtype=GRB.CONTINUOUS, name="rhs_numerator")
# Set numerator to 1
m.addConstr(rhs_numerator == 1, "numerator_def")

# Create an auxiliary variable to represent the right-hand side expression
aux_rhs_expression = m.addVar(vtype=GRB.CONTINUOUS, name="aux_rhs_expression")

# Add a constraint to define the auxiliary variable for the right-hand side expression
#m.addConstr(aux_rhs_expression, GRB.EQUAL, rhs_expression, "aux_rhs_expression_def")

# Introduce an auxiliary variable to represent 1 / rhs_expression
inv_rhs_expression = m.addVar(vtype=GRB.CONTINUOUS, name="inv_rhs_expression")

# Add the reciprocal constraint
#m.addGenConstrInv(aux_rhs_expression, inv_rhs_expression, "inv_rhs_constraint")
m.addConstr(inv_rhs_expression * rhs_expression == rhs_numerator, "inv_rhs_expression_def")

m.addConstr(P_0 == inv_rhs_expression)

This is the code that I've written.
Error line:

m.addConstr(inv_rhs_expression * rhs_expression == rhs_numerator, "inv_rhs_expression_def")

GurobiError: Invalid argument to QuadExpr multiplication