I think I have implemented the changes you suggested and below I attach the updated code.
From the first results I’ve seen that the prediction has improved dramatically, but I need some clarification regarding some parameters and functions in use in the neural network.
Note: I have a dataset of 16400 samples (.csv file consisting of a matrix 16400x6)
from numpy import loadtxt
import numpy as np
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers import Dense
from numpy.core.einsumfunc import einsum
import tensorflow as tf
from keras.layers import LayerNormalization
from numpy import array


# Model definition
def get_model(I_N_L_1, I_N_L_2, I_N_L_3, I_N_L_4, I_N_L_5, I_N_L_6, n_inputs, n_outputs):

	model = Sequential()
	model.add(Dense(I_N_L_1, input_dim=n_inputs, kernel_initializer='he_uniform', activation='relu'))
	model.add(LayerNormalization(axis=-1 , center=True , scale=True))
	model.add(Dense(I_N_L_2, activation='relu'))
	model.add(Dense(I_N_L_3, activation='relu'))
	model.add(Dense(I_N_L_4, activation='relu'))
	model.add(Dense(I_N_L_5, activation='relu'))
	model.add(Dense(I_N_L_6, activation='relu'))
	model.add(Dense(n_outputs))    
	model.compile(loss='mae', optimizer='adam', metrics=["mae"])
	return model


print('start reading CSV')
# Load Training set csv
dataset_EF = loadtxt('EF_from_Welding2.csv', delimiter=',')
x_train = dataset_EF[0:12000,0:6]
print('shape: ',x_train.shape)

# Load Label set csv
dataset_joints = loadtxt('Welding2.csv', delimiter=',')
y_train = dataset_joints[0:12000,0:6]
print('shape: ',y_train.shape)

# Validation set definition
x_val = dataset_EF[12001:14000,0:6]
print('shape: ',x_val.shape)

# Label of the validation set definition
y_val = dataset_joints[12001:14000,0:6]
print('shape: ',y_val.shape)


# Test set definition
x_test = dataset_EF[14001:,0:6]
print('shape: ',x_test.shape)
	
# Label of the test set definition
y_test = dataset_joints[14001:,0:6]
print('shape: ',y_test.shape)


print('end reading CSV')

# Number of nodes in the hidden layers
I_N_L_1 = 700
I_N_L_2 = 450
I_N_L_3 = 300
I_N_L_4 = 150
I_N_L_5 = 75
I_N_L_6 = 15

# Number of nodes in the input and output layers
n_inputs = 6
n_outputs = 6

print('start model and training')
# calling the "get_model" function
model = get_model(I_N_L_1, I_N_L_2, I_N_L_3, I_N_L_4, I_N_L_5, I_N_L_6, n_inputs, n_outputs)

# fit model
es = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=20, mode='min', restore_best_weights=True, verbose=1)
model.fit(x_train, y_train, validation_data=(x_val, y_val), verbose=1, epochs=100, callbacks=[es]) 
print('end model and training')

print('start saving model')
# saving the model
model.save("ModelloDiProva1.h5")
print('end saving model')


print('start predictions')
# Testing procedure
pred = []

# Computing the Prediction on the Test Set
for i in range(len(x_test)-1):
	b = [x_test[i][0], x_test[i][1], x_test[i][2], x_test[i][3], x_test[i][4], x_test[i][5]]
	ToBePredicted = array([b])
	Prediction = model.predict(ToBePredicted)
	a = [Prediction[0][0], Prediction[0][1], Prediction[0][2], Prediction[0][3], Prediction[0][4], Prediction[0][5]]
	pred.append(a)

print('end predictions')

print('start validation')
# Computing the mean vector of the error for each predicted joint trajectory
average_vector = []
sum = 0
average = 0
for j in range(6): # column
	for i in range(len(y_test)-1): #row
	
		sum = sum + (pred[i][j] - y_test[i][j])
	
	average = sum/len(y_test)
	average_vector.append(average)

	average = 0
	sum = 0

print('average_vector: ', average_vector)

# Computing the standard deviation vector of the error for each predicted joint trajectory
sum = 0
std_vector = []
for j in range(6): # colonne
	for i in range(len(y_test)-1): #righe
		sum = sum + ((pred[i][j] - y_test[i][j]) - average_vector[j])**2
	
	std = (sum/len(y_test))**(0.5)
	std_vector.append(std)

	std = 0
	sum = 0

print('std_vector: ', std_vector)
print('end validation')


as per your advice I have changed the “shape” of the neural network, starting with an initial layer of 700 nodes to arrive at the penultimate layer with 15 nodes. Do you have any advice or rule of thumb to use on the number of layers and the number of nodes per layer, in relation to the position of the layer within the network?
I have divided my dataset in the following way: 12000 values for the training set, 2000 values for the validation set and the remaining values for the test set. Is this a sensible choice or is it better to have all 3 sets with the same size?
Until now, each layer has the activation function “relu”. I have tried to use other activation functions, getting less precise results; in several examples on the internet I see that, for example, in the output layer a different activation function is used than in the previous layers. Is there a way to choose the best activation function based on the problem you are using? Why is a different activation function used in the output layer than in the other layers?
always searching between the examples in internet I have seen that the function “DropOut()” is used. I have understood that it is used to avoid the overfitting of the network by acting randomly on the weights stored in a particular layer of the network. Could it be useful to insert it also in my network? If yes, is it necessary to insert it between two specific layers or is it necessary to “go by attempts”?

Relatively to the normalization of the input, I have used the function “LayerNormalization”:

is it necessary to insert it only once, or in multiple layers of the network?
There is also a normalization function called “BatchNormalization”, but I could not understand the difference between the first and the second.

Thanks in advance for your attention