How to solve an overdetermined non linear set of equations numerically in python?

I am trying to solve a system of 4 exponential equations with two variables. However If I use fsolve python will only allow me two use as many equations as I have variables. But as I have infinitely many pairs of solutions (if only two equations are used) and I need to find the pair of variables that fits not only two but all four equations, fsolve does not seem to work.

So right know my code look something like this:

from scipy.optimize import fsolve
c_1 = w+(y/(x*R))
c_2 = offset - (c_1/(x*R)) #offset(1.05), w(10) and R(0.35) are constants 
def equations(p):
    x,y = p
    eq1 = (c_1*sp.exp(-y*R*1.017))/(y*R)+c_2-(x*1.017)/(y*R)-(5.1138*2*np.pi)
    eq2 = (c_1*sp.exp(-y*R*2.35))/(y*R)+c_2-(x*2.35)/(y*R)-(2.02457*4*np.pi)
    eq3 = (c_1*sp.exp(-y*R*2.683))/(y*R)+c_2-(x*2.683)/(y*R)-(6.0842178*4*np.pi)
    return (eq1,eq2,eq3)

x, y =  fsolve(equations, (1,1))

print (equations((x, y)))

However this will give me wildly different result depending on the initial guesses I give. I now want to so edit this code so that I can put additional equations to guarantee that I get the correct pair of x,y as a solution. Simply adding a third equation here in the into the function will return a TypeError of course so I was wondering how this can be done.

To add an arbitrary number of equations, I don't think you can use fsolve at all. Just run a least-squares minimisation, and make sure that you vectorise properly instead of rote repetition. The following does run and generate a result, though I have not assessed its accuracy.

import numpy as np
from scipy.optimize import minimize

offset = 1.05
w = 10
R = 0.35
a = np.array((-1.017, -2.350, -2.683))
b = np.array((
    -5.1138000*2,
    -2.0245700*4,
    -6.0842178*4,
))*np.pi


def equations(x: float, y: float) -> np.ndarray:
    c_1 = w + y/x/R
    c_2 = offset - c_1/x/R

    return (c_1*np.exp(a*y*R) + a*x)/y/R + c_2 + b


def error(xy: np.ndarray) -> np.ndarray:
    eqs = equations(*xy)
    return eqs.dot(eqs)


result = minimize(
    fun=error,
    x0=(1, 1),
    method='nelder-mead',
)
print(result.message)
print('x, y =', result.x)
print('equations:', equations(*result.x))

Optimization terminated successfully.
x, y = [0.19082078 0.13493941]
equations: [ 27.4082168   13.91087443 -42.00388728]