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

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.