Bilinear Interpolation Help

# Interpolation variables
u = model.addVar(lb=min(water_head_values), ub=max(water_head_values), name="water_head")
v = model.addVar(lb=min(flow_values), ub=max(flow_values), name="flow")
power_output = model.addVar(name="power_output")

# Spill flow variable
spill = model.addVar(lb=0, name="spill")

# Weight variables
lambda_u = model.addVars(len(water_head_values), lb=0, ub=1, name="lambda_u")
lambda_v = model.addVars(len(flow_values), lb=0, ub=1, name="lambda_v")

# Constraint: sum of weights equals 1
model.addConstr(lambda_u.sum() == 1, name="lambda_u_sum")
model.addConstr(lambda_v.sum() == 1, name="lambda_v_sum")

# Flow limit constraint
for i, u_val in enumerate(water_head_values):
    if u_val in flow_limits:
        max_flow = flow_limits[u_val]
        # Calculate effective power generation flow when water head is u
        model.addConstr((v - spill) <= max_flow, f"flow_limit_{u_val}")

# Interpolation formula: calculate power output using effective flow
interpolation_terms = []
for i, u_val in enumerate(water_head_values):
    for j, v_val in enumerate(flow_values):
        power_val = df.iloc[j, i]
        interpolation_terms.append(lambda_u[i] * lambda_v[j] * power_val)

# Output constraint
model.addConstr(power_output == quicksum(interpolation_terms), name="power_interpolation")

# Set non-convex optimization mode
model.setParam('NonConvex', 2)

# Objective function: maximize output and minimize spill
model.setObjective(power_output - 10 * spill, GRB.MAXIMIZE)

# Add SOS2 constraints
model.addSOS(GRB.SOS_TYPE2, [lambda_u[i] for i in range(len(water_head_values))])
model.addSOS(GRB.SOS_TYPE2, [lambda_v[j] for j in range(len(flow_values))])

# Start optimization
model.optimize()

I am currently facing an issue with a third variable in bivariate interpolation. I am interpolating power output based on a relationship table of water head, flow, and power output. However, there are varying flow limits at different water head levels, which results in empty values in the table. Additionally, when performing the interpolation for power output, it consistently indicates non-convexity. How can I linearize this relationship table for interpolation?