import gurobipy as gp
from gurobipy import GRB


def solve(points):
    pts_min_x = min(z[0] for z in points)
    pts_min_y = min(z[1] for z in points)
    offset_points = tuple((x - pts_min_x, y - pts_min_y) for x, y in points)

    model = gp.Model(f"geometric_median_{len(offset_points)}")
    x, y = model.addVar(name="x"), model.addVar(name="y")
    model.addConstr(0 <= x), model.addConstr(0 <= y)
    D = tuple(model.addVar(name=f"d_{i}") for i in range(len(offset_points)))

    for i, (point, d) in enumerate(zip(offset_points, D)):
        model.addConstr((x - point[0]) ** 2 + (y - point[1]) ** 2 <= d ** 2, f"dist_const_{i}")
        model.addConstr(0 <= d)

    model.setObjective(gp.quicksum(D), GRB.MINIMIZE)
    model.optimize()

    x_opt, y_opt = model.getVarByName("x").X, model.getVarByName("y").X
    return x_opt + pts_min_x, y_opt + pts_min_y


points = (
    (0, 0),
    (1, 4),
    (3, 2)
)
print(solve(points))