def METHOD_NAME(num_locations=50, num_customers=50):
    N = list(range(num_locations))  # warehouse locations
    M = list(range(num_customers))  # customers
    d = dict()  # distances from warehouse locations to customers
    for n in N:
        for m in M:
            d[n, m] = np.random.randint(low=1, high=100)
    max_num_warehouses = 2
    model = pyo.ConcreteModel(name="(WL)")
    model.P = pyo.Param(initialize=max_num_warehouses, mutable=True)
    model.x = pyo.Var(N, M, bounds=(0, 1))
    model.y = pyo.Var(N, bounds=(0, 1))
    def obj_rule(mdl):
        return sum(d[n, m] * mdl.x[n, m] for n in N for m in M)
    model.obj = pyo.Objective(rule=obj_rule)
    def demand_rule(mdl, m):
        return sum(mdl.x[n, m] for n in N) == 1
    model.demand = pyo.Constraint(M, rule=demand_rule)
    def warehouse_active_rule(mdl, n, m):
        expr = LinearExpression(
            constant=0, linear_coefs=[1, -1], linear_vars=[mdl.x[n, m], mdl.y[n]]
        )
        return expr <= 0
    model.warehouse_active = pyo.Constraint(N, M, rule=warehouse_active_rule)
    def num_warehouses_rule(mdl):
        return sum(mdl.y[n] for n in N) <= model.P
    model.num_warehouses = pyo.Constraint(rule=num_warehouses_rule)
    return model