Lagrange Multipliers /w sympy

Question

I am currently trying to find the maximum radius of a circle I can manifest between existing circles around it.

i.e. I'm trying to find not only the maximum radius, but the center point most suited for it over a specific given straight line.

In order to find said maxima I'm trying to implement a generalized Lagrange multipliers solution using sympy.

If "n" represents the amount of constraints I have, then I was able to:

• Create n symbols generator.
• Perform the necessary nth-gradient over the Lagrange function
• Manifest the required inequalities (from constraints) to achieve the list of equalities and inequalities needed to be solved.

The code:

from sympy import S
from sympy import *
import sympy as smp


#Lagrange Multipliers

def sympy_distfun(cx,cy,radius):
       
       x,y=smp.symbols('x y',real=True)
       return sqrt((x-cx)**2+(y-cy)**2)-radius

def sympy_circlefun(cx,cy,radius):
       
       x,y=smp.symbols('x y',real=True)
       return (x-cx)**2+(y-cy)**2-radius**2

def sympy_linefun(slope,b):
       
       x,y=smp.symbols('x y',real=True)
       return slope*x+b-y

def lagrange_multiplier(objective,constraints):
     
     x,y=smp.symbols('x y',real=True)
     a=list(smp.symbols('a0:%d'%len(constraints),real=True))
     
     cons=[constraints[i]*a[i] for i in range(len(a))]
     L=objective+(-1)*sum(cons)
     gradL=[smp.diff(L,var) for var in [x,y]+a]
     constraints=[(con)>= 0 for con in constraints]
     eqs=gradL+constraints
     
     
     vars=a+[x,y]
     solution=smp.solve(eqs[0],vars)
     #solution=smp.solveset(eqs,vars)
     print(solution)
       
       
       
line=sympy_linefun(0.66666,-4.3333)
dist=sympy_distfun(11,3,4)
circlefunc1=sympy_circlefun(11,3,4)
circlefunc2=sympy_circlefun(0,0,3)
lagrange_multiplier(dist,[line,circlefunc1,circlefunc2])

But, when using smp.solveset(eqs,vars) I encounter the error message:

ValueError: [-0.66666*a0 - a1*(2*x - 22) - 2*a2*x + (x - 11)/sqrt((x - 11)**2 + (y - 3)**2), a0 - a1*(2*y - 6) - 2*a2*y + (y - 3)/sqrt((x - 11)**2 + (y - 3)**2), -0.66666*x + y + 4.3333, -(x - 11)**2 - (y - 3)**2 + 16, -x**2 - y**2 + 9, 0.66666*x - y - 4.3333 >= 0, (x - 11)**2 + (y - 3)**2 - 16 >= 0, x**2 + y**2 - 9 >= 0] is not a valid SymPy expression

When using: solution=smp.solve(eqs[0],vars) to try and solve one equation, it sends sympy into a CPU crushing frenzy and it obviously fails to complete the calculation. I made sure to declare all variables as real so i fail to see why it takes so long to solve.

Would like to understand what I'm missing when it comes to handling multiple inequalities with sympy, and if there is a more optimized faster way to solve Lagrange multiplication I'd love to give it a try