In general case, the problem of non-crossing linear quantile regresion *has no solutions*. 

Indeed, either your quantile lines are parallel (and trivial), or they do cross *somewhere*. 

By placing restrictions, you could ensure that the lines do not cross in the training sample, but it will in no way guarantee that crossing will not occur at the next observation you see after training the model. And if the quantile lines intersect *within the training sample*, it most probably means that your model is specified incorrectly: either mean or standard deviation change nonlinearly.

If you still want to estimate such a model, I would propose fitting it as a neural network. 

Your model should take input $X$, multiply it by matrix(!) $\beta$ to get a matrix of forecasts $f=x\beta$ with size $n \times k$, where $n$ is number of observations and $k$ is number of estimated quantiles. I assume that quantile percentages $q$ are in increasing order. You should minimize the function 

$$
L = \sum_{i=1}^n \left(\sum_{j=1}^k  \max(q_j(y_i-f_{ij}), (q_j-1)(y_i-f_{ij})) +  \sum_{j=1}^{k-1} \alpha \max(0, \delta - (f_{j+1}-f_{j})) \right)
$$

The first term in the inner sum is just the ordinary quantile regression loss. The second term is the penalty that is applied if two consecutive quantile predictions differ by less then $\delta$.

Minimizing this function by gradient descent will give you your non-crossing quantile lines (if $\alpha$ is large enough). But I still warn you that if "natural" quantile lines intersect, then ther might be problems with the functional form of your model.

There is a Python example. As for me, both restricted and unrestricted versions look ugly.

    # import everything
    import keras
    from keras import backend as K
    from keras import Sequential
    from keras.layers import Dense
    import numpy as np
    import matplotlib.pyplot as plt

    # create the dataset
    n = 200
    np.random.seed(1)
    X = np.sort(np.random.normal(size=n))[:, np.newaxis] + 4
    y = 5 + X.ravel()*(1 + np.random.normal(size=n)*0.2)
    quantiles = np.array([0.1, 0.25, 0.5, 0.75, 0.9])

    # define loss function
    def quantile_ensemble_loss(q, y, f, alpha=100, margin=0.1):
        error = (y - f)
        quantile_loss = K.mean(K.maximum(q*error, (q-1)*error))
        diff = f[:, 1:] - f[:, :-1]
        penalty = K.mean(K.maximum(0.0, margin - diff)) * alpha
        return quantile_loss + penalty

    # fit two models
    for i, alpha, name in [(1, 0, 'w/o penalty'), (2, 10, 'with_penalty')]:
        model = Sequential()
        model.add(Dense(len(quantiles), input_dim=X.shape[1]))
        model.compile(loss=lambda y,f: quantile_ensemble_loss(quantiles,y,f,alpha), optimizer=keras.optimizers.RMSprop(lr=0.003))
        model.fit(X, y, epochs=3000, verbose=0, batch_size=100);
        plt.subplot(1,2,i)
        plt.scatter(X.ravel(), y, s=1)
        plt.plot(X.ravel(), model.predict(X))
        plt.title(name)

    plt.show()

[![enter image description here][1]][1]


  [1]: https://i.sstatic.net/SaByH.png