Python code

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy
numpy.random.seed(3)
nMess = 18
x0 = numpy.random.rand(nMess)
x0 = x0-min(x0); x0 *= 10 / max(x0)
x1 = numpy.random.rand(nMess)*10
x1 = x1-min(x1); x1 *= 10 / max(x1)
y = numpy.random.randn(nMess)*10
x = numpy.array(list(zip(x0,x1)))
def base(x0, x1, n):
 return numpy.array([x0**i0.ravel() * x1**i1.ravel() for i0, i1 in numpy.mgrid[:n,:n].reshape(2,n*n).T]).T
nBasisProAchse = 3
A = base(x0, x1, nBasisProAchse)
c = numpy.linalg.inv(A.T @ A) @ A.T @ y
yModell = A @ c
x0neu, x1neu = numpy.meshgrid(numpy.linspace(0, 10, num=20), numpy.linspace(0, 10, num=20))
yneu = base(x0neu, x1neu, nBasisProAchse) @ c
fig = plt.figure(figsize=(4.5,4.5))
ax = fig.add_subplot(111, projection='3d')
p2 = ax.plot_wireframe(x1neu, x0neu, yneu, color='r', rstride=1, cstride=1, label='Modellfunktion', linewidth=1, zorder=0)
p1 = ax.scatter(x0, x1, y, label='Messpunkte', zorder=2)
for i in range(len(y)):
 ax.plot([x0[i]]*2, [x1[i]]*2, [y[i], yModell[i]], color='#60c060', linewidth=1.5, label='Residuum' if i==0 else None, zorder=1)
lx0 = ax.set_xlabel('x1')
lx1 = ax.set_ylabel('x2')
ly = ax.set_zlabel('y')
leg = plt.legend(frameon=True)
plt.tight_layout()
plt.show()
plt.savefig('MDKQ_3D.svg')