It looks like the method used to do the fits is stochastic. From the pwlf source code, fit calls the function scipy.optimize.differential_evolution. The Scipy docs for this function describe it as a stochastic method to determine a global minimum of a function.

It appears that with the default settings used by pwlf, its not quite converged or at least its approaching convergence from different directions each time. You maybe able to edit the keywords so that the optimization will do a more stringent search for the true minimum. For example, atol sets the absolute tolerance for convergence. You also set seed to a consistent value to at least get the same result every time, even if it may not have found the "true" minimum. Just be aware that if you change one of the keywords, you need to explicitly write out all the others, as the **kwargs passed into fit appear to overwrite the defaults that would be passed to differential_evolution.

The pwlf library has an example specifying bounds for the breakpoint search, which should work for your example.

Unknown breakpoints is a tough optimization problem, I could foresee using soft constraints to modify the loss function to penalize small N breaks, hard constraints I think will be difficult. So I think restricting search space is the best bet.

Some changepoint formulations do this explicitly by marginalizing out the changepoint (so only search over a discrete set of change point locations), here is an example in pystan I have written.

You can use numpy.piecewise() to create the piecewise function and then use curve_fit(), Here is the code

from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15], dtype=float)
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03])

def piecewise_linear(x, x0, y0, k1, k2):
    return np.piecewise(x, [x < x0], [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0])

p , e = optimize.curve_fit(piecewise_linear, x, y)
xd = np.linspace(0, 15, 100)
plt.plot(x, y, "o")
plt.plot(xd, piecewise_linear(xd, *p))

the output:

For an N parts fitting, please reference segments_fit.ipynb