# Revision history [back]

It's ok, I found the solution using MILP:

def center_inscribed_sphere_polytope(polytope, solver=None):
r"""
Return the center and radius of maximal inscribed sphere

INPUT:

- polytope -- polytope

OUTPUT:

EXAMPLES::

sage: P = polytopes.associahedron(['A',3])
sage: center_inscribed_sphere_polytope(P)
([0.03553390593273766, 0.5355339059327378, 0.03553390593273766],
1.4644660940672622)
sage: center_inscribed_sphere_polytope(P + vector((10,10,10)))
([10.035533905932738, 10.535533905932738, 10.035533905932738],
1.4644660940672622)

::

sage: P = Polyhedron([(0,0), (1,0), (0,10)])
sage: center_inscribed_sphere_polytope(P)
([0.47506218943955486, 0.47506218943955486], 0.47506218943955486)

"""
from math import sqrt
from sage.numerical.mip import MixedIntegerLinearProgram
p = MixedIntegerLinearProgram(solver=solver)
x = p.new_variable()
d = p.new_variable()
for ineq in polytope.inequalities_list():
constant = ineq
coeffs = ineq[1:]
norm = sqrt(sum(a**2 for a in coeffs))
S = p.sum(a/norm*x[i] for (i,a) in enumerate(coeffs))