I agree with Mikael's answer that you should try to make your expression as concise as possible. I have an example of your real data shape, that combines linear functions and a constant function.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define step function
def step_func(x, a, c, f1, f2):
    result = []
    for xi in x:
        if xi >= 0 and xi < f1:
            result.append(a * xi)
        elif xi >= f1 and xi < f2:
            result.append(a * f1)
        elif xi >= f2:
            result.append(c * xi + a * f1 - f2 * c)
        else:
            result.append(0)
    return np.asarray(result)

# Generating some demo data
x_data = np.linspace(0, 10, 100)
y_data = step_func(x_data, 2, 1, 3, 7) + np.random.normal(0, 0.5, x_data.size)  # 添加噪声

# Fitting function
def fitting_func(x, a, c, f1, f2):
    return step_func(x, a, c, f1, f2)

# Fitting for data
popt, pcov = curve_fit(fitting_func, x_data, y_data, p0=[1, -1, 4, 6])

# Obtain fitting parameters
a_fit, c_fit, f1_fit, f2_fit = popt
print(f'Fitting param: a={a_fit}, c={c_fit}, f1={f1_fit}, f2={f2_fit}')

# Generat fitting results
y_fit = fitting_func(x_data, *popt)

# Plot
plt.scatter(x_data, y_data, label='Data', s=10)
plt.plot(x_data, y_fit, color='red', label='Fitted')
plt.vlines([f1_fit, f2_fit], ymin = min(x_data), ymax = max(x_data) )
# plt.xlim(0,16)
# plt.ylim(0,1)
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Step func fitting')
plt.show()

The above returns:

Fitting param: a=1.9532567200579396, c=1.085588494015425, f1=3.084102835100397, f2=7.1628504872999095