I have the following python code to perform a straight line movement of a robot end-effector in operational space. How could I set the end-effector velocity?
"""
This file demonstrates the use of the Cartesian planner to generate a trajectory
"""
import sys
import os
if os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "src") not in sys.path:
    sys.path.append(os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "src"))
import numpy as np
import time
## pin
import pinocchio as pin
from pinocchio.visualize import MeshcatVisualizer
from pinocchio.robot_wrapper import RobotWrapper 
##  POLYTHESIS
from polythesis.core.utils import generate_mod_file
from polythesis.robot.robot_loader import loadABB_IRB7600325, loadABB_IRB7600325v2, loadABB_IRB4600255, loadABB_IRB4600255v2
# from polythesis.ik.differential_ik import DifferentialIk, DifferentialIkOptions
from polythesis.ik.mydifferential_ik import DLS_Options, DLS_Ik
# from polythesis.planning.cartesian_planner import (CartesianPlanner,CartesianPlannerOptions)
from polythesis.planning.my_cartesian_planner import (CartesianPlanner,CartesianPlannerOptions)
from polythesis.visualization.meshcat_utils import visualize_frames
from polythesis.plot.plot_results import PlotResults
import polythesis.EGM.Core as egm
import matplotlib.pyplot as plt

##   Add access if it is not in the system path.


# CONSTANTS
VERBOSE=True                      # Print out debug information
VISUAL=True                       # Visualize the robot in MeshCat
PLOT=True                         # Plot results
# Constants for the simulation
dt=0.04                          # [s] Time step for the simulation
SIMU_REAL_TIME=True               # Set to True if you want to run the simulation in real time
ndt=10                            # [s] Number of integration steps for each control loop
# Constants for the REAL robot controller
CONTROL_REAL_ROBOT=False           # Set to True if you want to control the real robotABB
IP_CONTROLLER='172.18.206.202'    # IP address of the ABB controller
PORT_CONTROLLER=6511              # Port of the ABB controller
# Constants for the robot max values
MAX_LINEAR_VELOCITY=0.4           # [m/s] Maximum TCP linear velocity
MAX_LINEAR_ACCELERATION=0.4       # [m/s^2] Maximum TCP linear acceleration 
MAX_ANGULAR_VELOCITY=0.5          # [rad/s] Maximum TCP angular velocity
MAX_ANGULAR_ACCELERATION=0.4      # [rad/s^2] Maximum TCP angular acceleration
# Other constants
LINE_WIDTH = 80

# Load your robot's URDF model
robot = loadABB_IRB7600325v2() # Available models: loadABB_IRB4600255,loadABB_IRB4600255v2, loadABB_IRB7600325,loadABB_IRB7600325v2
model = robot.model
data = robot.data

# Retrieve joint limits from the robot model
lower_limits = model.lowerPositionLimit
upper_limits = model.upperPositionLimit

# Specify the end-effector frame name as defined in your URDF
end_effector_frame_name = 'tool0'  # Replace with your end-effector frame name
frame_id = model.getFrameId(end_effector_frame_name)

# Initial joint configuration q0 (robot's home position)
# Ensure the initial configuration is within joint limits
q0 = pin.neutral(model)
q0 = np.minimum(np.maximum(q0, lower_limits), upper_limits)

# Time parameters
dt = 0.01  # Time step (seconds)
T = 7.0    # Total time for the motion (seconds)
num_steps = int(T / dt)

## FOR THE ABB IRB 7600-325
start_TCP_pose = pin.SE3(pin.rpy.rpyToMatrix(np.array([0,np.pi/2,0])), np.array([2.21, -1.0, 2.02]))
end_TCP_pose = pin.SE3(pin.rpy.rpyToMatrix(np.array([0,np.pi/2,0])), np.array([2.21, +1.0, 2.02]))

# Function to interpolate between two poses
def interpolate_pose(start_pose, end_pose, alpha):
    # Linear interpolation for position
    position = (1 - alpha) * start_pose.translation + alpha * end_pose.translation
    # Spherical linear interpolation (slerp) for orientation
    quat_start = pin.Quaternion(start_pose.rotation)
    quat_end = pin.Quaternion(end_pose.rotation)
    quat_interp = quat_start.slerp(alpha, quat_end)
    orientation = quat_interp.matrix()
    return pin.SE3(orientation, position)

# Function to compute error between current and desired poses
def compute_pose_error(current_pose, desired_pose):
    # Position error
    e_translation = desired_pose.translation - current_pose.translation
    # Orientation error (using log of rotation matrix)
    R_current = current_pose.rotation
    R_desired = desired_pose.rotation
    R_error = R_desired @ R_current.T
    e_rotation = pin.log3(R_error)
    # Concatenate translation and rotation errors
    e = np.concatenate((e_translation, e_rotation))
    return e, e_translation, e_rotation

# Inverse Kinematics to find initial joint configuration q0
max_ik_iterations = 1000
tolerance = 1e-6
q0 = pin.neutral(model)  # Initial guess for IK

for iteration in range(max_ik_iterations):
    # Compute current TCP pose using forward kinematics
    pin.forwardKinematics(model, data, q0)
    pin.updateFramePlacement(model, data, frame_id)
    oMf = data.oMf[frame_id]  # Current end-effector pose
    x_current = pin.SE3(oMf.rotation, oMf.translation)

    # Compute error between current and desired pose
    e,_,_ = compute_pose_error(x_current, start_TCP_pose)
    error_norm = np.linalg.norm(e)

    if error_norm < tolerance:
        print(f"IK Converged in {iteration} iterations.")
        break

    # Compute Jacobian at q0
    J = pin.computeFrameJacobian(model, data, q0, frame_id, pin.ReferenceFrame.LOCAL_WORLD_ALIGNED)
    J = J[:6, :]

    # Damped least squares
    U, S, Vh = np.linalg.svd(J, full_matrices=False)
    lambda_damping = 0.01
    S_damped_inv = np.array([s / (s**2 + lambda_damping**2) for s in S])
    J_damped_pinv = Vh.T @ np.diag(S_damped_inv) @ U.T

    # Compute delta_q
    delta_q = J_damped_pinv @ e

    # Update q0
    q0 = pin.integrate(model, q0, delta_q)

    # Apply joint limits
    q0 = np.minimum(np.maximum(q0, lower_limits), upper_limits)
else:
    print(f"IK did not converge in {max_ik_iterations} iterations.")

# Initialize joint configuration
q = q0.copy()
q_list = [q.copy()]        # List to record joint positions
time_list = [0.0]          # List to record time steps
# Initialize lists to record errors, damping factors, and singular values
e_translation_norm_list = []
e_rotation_norm_list = []
lambda_damping_list = []
S_min_list = []
x_list = []

# Main control loop
for i in range(num_steps):
    # Compute time and interpolation factor
    t = i * dt
    alpha = t / T  # Interpolation ratio from 0 to 1

    # Interpolate to get desired TCP pose at current time
    x_desired = interpolate_pose(start_TCP_pose, end_TCP_pose, alpha)

    # Compute current TCP pose using forward kinematics
    pin.forwardKinematics(model, data, q)
    pin.updateFramePlacement(model, data, frame_id)
    oMf = data.oMf[frame_id]  # Current transform of the end-effector frame
    x_current = pin.SE3(oMf.rotation, oMf.translation)
    x_list.append(x_current)

    # Compute error between current and desired poses
    e, e_translation, e_rotation = compute_pose_error(x_current, x_desired)
    # Record the norms of translation and rotation errors
    e_translation_norm_list.append(np.linalg.norm(e_translation))
    e_rotation_norm_list.append(np.linalg.norm(e_rotation))

    # Compute the Jacobian of the end-effector
    J = pin.computeFrameJacobian(model, data, q, frame_id, pin.ReferenceFrame.LOCAL_WORLD_ALIGNED)
    nn = J.shape[0] # Number of rows 
    nm = J.shape[1] # Number of columns

    # Compute SVD of the Jacobian
    U, S, Vh = np.linalg.svd(J, full_matrices=False)

    # Determine the damping factor lambda based on singular values
    singular_threshold = 0.05
    # Minimum singular value
    S_min = S[-1]
    S_min_list.append(S_min)  # Record S_min

    # If near singularity, increase damping
    if S_min < singular_threshold:
        lambda_damping = (1 - (S_min / singular_threshold)) * 0.01
    else:
        lambda_damping = 0.0  # No damping needed
    # Small regularization to avoid division by zero
    #lambda_damping = max(lambda_damping, 1e-6)
    # Record lambda_damping
    lambda_damping_list.append(lambda_damping)  
    
    # Normalizing matrices
    Nx = np.identity(nm)
    Nq = np.identity(nn)
    # Weight matrices
    W1 = Nx.T @ Nx
    W2 = lambda_damping**2 * (Nq.T @ Nq)

    JJT = J @(W1@ J.T) + W2
    delta_q = J.T@ (W1 @ np.linalg.solve(JJT, e))

    # Optional: Limit the maximum joint update step to avoid too large movements
    max_delta_q = 0.05  # Maximum joint update (radians)
    norm_delta_q = np.linalg.norm(delta_q)
    if norm_delta_q > max_delta_q:
        delta_q *= max_delta_q / norm_delta_q

    # Update joint positions using integration
    q_new = pin.integrate(model, q, delta_q)

    # Apply joint limits
    q_new = np.minimum(np.maximum(q_new, lower_limits), upper_limits)

    # Check if the joint update was feasible
    if not np.allclose(q_new, q + delta_q, atol=1e-6):
        # Joint limits have modified the joint update, which may affect the end-effector pose
        # Consider recomputing delta_q with joint limits projection (advanced topic)
        pass  # For this simple implementation, we proceed with the limited q_new

    # Update joint positions
    q = q_new

    # Record the joint positions and time
    q_list.append(q.copy())
    time_list.append(t + dt)

    # (Optional) Print progress
    if i % 100 == 0:
        print(f"Step {i}/{num_steps}, Time {t:.2f}s, Damping {lambda_damping:.4f}")

# The q_list now contains the joint configurations over time, considering joint limits

## PLOT THE RESULTS
# Convert lists to numpy arrays for easier indexing
q_array = np.array(q_list)  # Shape: (num_steps + 1, num_joints)
e_translation_norm_array = np.array(e_translation_norm_list)
e_rotation_norm_array = np.array(e_rotation_norm_list)
lambda_damping_array = np.array(lambda_damping_list)
S_min_array = np.array(S_min_list)

# Compute joint velocities
q_dot = np.diff(q_array, axis=0) / dt  # Shape: (num_steps, num_joints)
time_list_qdot = time_list[1:]  # Corresponding time steps

# Compute joint accelerations
q_ddot = np.diff(q_dot, axis=0) / dt  # Shape: (num_steps - 1, num_joints)
time_list_qddot = time_list[2:]  # Corresponding time steps

# Plot TCP positions over time
plt.figure(figsize=(10, 6))
plt.plot(time_list[1:], [x.translation[0] for x in x_list], label="X")
plt.plot(time_list[1:], [x.translation[1] for x in x_list], label="Y")
plt.plot(time_list[1:], [x.translation[2] for x in x_list], label="Z")
plt.xlabel("Time [s]")
plt.ylabel("Position [m]")
plt.title("End-Effector Position Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Plot joint positions over time
plt.figure(figsize=(10, 6))
num_joints = q_array.shape[1]
for joint_index in range(num_joints):
    plt.plot(time_list, q_array[:, joint_index], label=f"Joint {joint_index + 1}")
plt.xlabel("Time [s]")
plt.ylabel("Joint Angles [rad]")
plt.title("Joint Angles Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Plot joint velocities over time
plt.figure(figsize=(10, 6))
for joint_index in range(num_joints):
    plt.plot(time_list_qdot[1:], q_dot[1:, joint_index], label=f"Joint {joint_index + 1}")
plt.xlabel("Time [s]")
plt.ylabel("Joint Velocities [rad/s]")
plt.title("Joint Velocities Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Plot joint accelerations over time
plt.figure(figsize=(10, 6))
for joint_index in range(num_joints):
    plt.plot(time_list_qddot[1:], q_ddot[1:, joint_index], label=f"Joint {joint_index + 1}")
plt.xlabel("Time [s]")
plt.ylabel("Joint Accelerations [rad/s²]")
plt.title("Joint Accelerations Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Plot position error norm over time
plt.figure(figsize=(10, 6))
plt.plot(time_list[1:-1], e_translation_norm_array[1:], label="Position Error Norm")
max_pos_error_norm=np.max(e_translation_norm_array[1:])
plt.axhline(max_pos_error_norm, color='r', linestyle='--', label=f"Max Value")
plt.text(time_list[-1], max_pos_error_norm, f'{max_pos_error_norm:.5f}', color='r', ha='right')
plt.xlabel("Time [s]")
plt.ylabel("Position Error Norm [m]")
plt.title("Position Error Norm Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Plot orientation error norm over time
plt.figure(figsize=(10, 6))
plt.plot(time_list[:-1], e_rotation_norm_array, label="Orientation Error Norm")
max_rot_error_norm=np.max(e_rotation_norm_array)
plt.axhline(max_rot_error_norm, color='r', linestyle='--', label=f"Max Value")
plt.text(time_list[-1], max_rot_error_norm, f'{max_rot_error_norm:.5f}', color='r', ha='right')
plt.xlabel("Time [s]")
plt.ylabel("Orientation Error Norm [rad]")
plt.title("Orientation Error Norm Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Plot lambda damping values over time
plt.figure(figsize=(10, 6))
plt.plot(time_list[:-1], lambda_damping_array, label="Lambda Damping")
max_lambda=np.max(lambda_damping_array)
plt.axhline(max_lambda, color='r', linestyle='--', label=f"Max Value")
plt.text(time_list[-1], max_lambda, f'{max_lambda:.5f}', color='r', ha='right')
plt.xlabel("Time [s]")
plt.ylabel("Lambda Damping")
plt.title("Lambda Damping Values Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Plot smallest singular value over time
plt.figure(figsize=(10, 6))
plt.plot(time_list[:-1], S_min_array, label="Smallest Singular Value")
plt.xlabel("Time [s]")
plt.ylabel("Singular Value")
plt.title("Smallest Singular Value Over Time")
plt.legend()
plt.grid(True)
plt.show()

# Optional: Visualize the robot's motion
for idx in range(q_array.shape[0]):
    robot.viz.display(q_array[idx, :])
    time.sleep(dt)

input("Press 'q' to exit the program...")