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, or you apply a wrong cost function when fitting the model.
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_{i,j+1}-f_{ij})) \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. Maybe you would prefer quantile estimates of Random Forest (like quantregForest in R), which are always consistent.
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()
