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

numpy.piecewise can do this.

piecewise(x, condlist, funclist, *args, **kw)

Evaluate a piecewise-defined function.

Given a set of conditions and corresponding functions, evaluate each function on the input data wherever its condition is true.

An example is given on SO here. For completeness, here is an example:

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

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, 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))

You could directly copy the segments_fit implementation

from scipy import optimize

def segments_fit(X, Y, count):
    xmin = X.min()
    xmax = X.max()

    seg = np.full(count - 1, (xmax - xmin) / count)

    px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
    py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init])

    def func(p):
        seg = p[:count - 1]
        py = p[count - 1:]
        px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        return px, py

    def err(p):
        px, py = func(p)
        Y2 = np.interp(X, px, py)
        return np.mean((Y - Y2)**2)

    r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')
    return func(r.x)

Then you apply it as follows

import numpy as np;

# mimic your data
x = np.linspace(0, 50)
y = 50 - np.clip(x, 10, 40)

# apply the segment fit
fx, fy = segments_fit(x, y, 3)

This will give you (fx,fy) the corners your piecewise fit, let's plot it

import matplotlib.pyplot as plt

# show the results
plt.figure(figsize=(8, 3))
plt.plot(fx, fy, 'o-')
plt.plot(x, y, '.')
plt.legend(['fitted line', 'given points'])

EDIT: Introducing constant segments

As mentioned in the comments the above example doesn't guarantee that the output will be constant in the end segments.

Based on this implementation the easier way I can think is to restrict func(p) to do that, a simple way to ensure a segment is constant, is to set y[i+1]==y[i]. Thus I added xanchor and yanchor. If you give an array with repeated numbers you can bind multiple points to the same value.

from scipy import optimize

def segments_fit(X, Y, count, xanchors=slice(None), yanchors=slice(None)):
    xmin = X.min()
    xmax = X.max()
    seg = np.full(count - 1, (xmax - xmin) / count)

    px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
    py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init])

    def func(p):
        seg = p[:count - 1]
        py = p[count - 1:]
        px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        py = py[yanchors]
        px = px[xanchors]
        return px, py

    def err(p):
        px, py = func(p)
        Y2 = np.interp(X, px, py)
        return np.mean((Y - Y2)**2)

    r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')
    return func(r.x)

I modified a little the data generation to make it more clear the effect of the change

import matplotlib.pyplot as plt
import numpy as np;

# mimic your data
x = np.linspace(0, 50)
y = 50 - np.clip(x, 10, 40) + np.random.randn(len(x)) + 0.25 * x
# apply the segment fit
fx, fy = segments_fit(x, y, 3)
plt.plot(fx, fy, 'o-')
plt.plot(x, y, '.k')
# apply the segment fit with some consecutive points having the 
# same anchor
fx, fy = segments_fit(x, y, 3, yanchors=[1,1,2,2])
plt.plot(fx, fy, 'o--r')
plt.legend(['fitted line', 'given points', 'with const segments'])

It is a problem with types, you have to change the following line, so that the x is given as floats:

x = np.array([-3, -2, -1, 0, 1, 2, 3]).astype(np.float)

otherwise the piecewise_linear will might end up casting the types.

Just to be on the safe side you could also make the initial points float here:

popt_piecewise, pcov = curve_fit(piecewise_linear, x, y, p0=[0.1, 0.1, -5., 5.])

To finish this up here, I'll share my own final solution to the problem. In order to stay close to my original question, you just have to define the vectorized function yourself and not use np.vectorize.

import scipy.optimize as so
import numpy as np

def fitfunc(x,p):
   if x>p:
      return x-p
   else:
      return -(x-p)

fitfunc_vec = np.vectorize(fitfunc) #vectorize so you can use func with array

def fitfunc_vec_self(x,p):
  y = np.zeros(x.shape)
  for i in range(len(y)):
    y[i]=fitfunc(x[i],p)
  return y


x=np.arange(1,10)
y=fitfunc_vec_self(x,6)+0.1*np.random.randn(len(x))

popt, pcov = so.curve_fit(fitfunc_vec_self, x, y) #fitting routine that gives error
print popt
print pcov

Output:

[ 6.03608994]
[[ 0.00124934]]

Within your optimize.curve_fit() you need to specify some initial guess for the fitting using p0 = [] where you input your initial guesses into p0, the documentation of which can be found here.

In terms of the example you have given, you already have the values of x0, y0, k1, k2 as you used them to calculate y, therefore just input these into your curve_fit:

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

x=np.arange(0.0,100.0,1.0)
y=piecewise_linear( x, 45.0, 2025.0, 1.0, 2.0)

#insert the initial guesses into curve_fit below using p0 = [...]
popt , pcov = optimize.curve_fit(piecewise_linear, x, y, p0=[45, 2000, 1, 2])
tau = np.linspace(x[0], x[-1], 200)
perr = np.sqrt(np.diag(pcov))

print (popt)
print (perr)
print (pcov)

plt.plot(x, y, 'b+')
plt.plot(tau, piecewise_linear(tau, *popt),'r')
plt.loglog()
plt.show()

This give the following graph:

You could wrap the curve fit problem into an outer optimization problem where you optimize x0 and x1.

Try:

def piecewise_func(p,x):
    a0,a1,c0,c1,x0,x1 = p

    b3 = ((c1-(a0+a1*x0-c0-c1*x1)/(x0-x1) + (x0+x1)/(2*(x0-x1)) * (a1-c1) -x1/(x0-x1)*(a1-c1))/(3.*x1**2))/ (1.-(-1.5*(x0+x1)**2 + (x0**3-x1**3)/(x0-x1) + 3.*x1*(x0+x1))/(3.*x1**2))
    b2 = (a1 - c1 - 3.*b3*(x0**2 - x1**2))/(2.*(x0-x1))
    b1 = (a0+a1*x0-c0-c1*x1-b2*(x0**2-x1**2)-b3*(x0**3-x1**3))/(x0-x1)
    b0 = a0 + a1*x0 - b1*x0 - b2*x0**2 - b3*x0**3

    condlist = [x < x0, (x >= x0) & (x < x1), x >= x1]
    funclist = [lambda x: a0 + a1*x,
                lambda x: b0 + b1*x + b2*x**2 + b3*x**3,
                lambda x: c0 + c1*x]
    return np.piecewise(x, condlist, funclist)


def optFunc(p,extra):
    x,y = extra
    return piecewise_func(p,x)-y


fit,ier = optimize.leastsq(optFunc,[0.,1.,0.,1.,0.01,0.1],[x,np.log(y)])
print fit

plt.plot(x,np.log(y),'k+')
plt.plot(x,piecewise_func(fit,x))

This should yield the following fit: click here