Fitting several curves simultaneously in Python with correlated parameters

Question

I am trying to fit several curves simultaneously in Python. The parameters of these curves are related to each other. To explain what I meant, here is a minimal working example:

import numpy as np
from scipy import optimize

def linear(x, slope, y0):
    return slope * x + y0

slope, y0 = [1, 3, 6], [6, 4, 1]
x = np.linspace(0, 1, 101)

for i in range(len(y0)):
    y = linear(x, slope[i], y0[i]) + np.random.random(len(x)) * 0.1
    popt = optimize.curve_fit(linear, x, y, sigma=x*0.+0.1)[0]
    print(popt)

These are simply straight lines. The slopes and the y-intersections nicely follow an ascending / descending pattern. For simple case such as this one, fitting the curves separately works at a reasonably accurate level. However, I wonder if there is any way to fit these curves simultaneously, knowing that the parameters of one curve are loosely correlated to the adjacent curves, in order to improve the quality of the fit.

I looked up a little bit but failed to find the resources that I need. Most of the examples of simultaneous fitting require the parameters to be the same / to be fixed, which is not what I am looking for here. A simple way is to set some boundaries for the parameter based on the results of adjacent curves. In this case I will have to run the fit several times, maybe starting from the curve with the highest S/N. This seems very inefficient to me and I feel like there should be a neat way to tackle with problem like this.

Thank you very much!