Loss of precision in calculations involving fractions and floating-point numbers in Python

import sympy as sp
from fractions import Fraction

new_matrix_aux = [['X',   'B', 'X1', 'X2', 'X3', 'X4', 'X5', 'U1', 'U2'], 
                  ['X2', 10/3,    0,    1,  2/3,  1/3, -2/3, -1/3,  2/3],  
                  ['X1',  4/3,    1,    0,  2/3, -2/3,  1/3,  2/3, -1/3]]

z_function = "Z = 1 * X1 + 2 * X2 + 3 * X3 + 0 * X4 + 0 * X5 + M * U1 + M * U2"


lista_de_operaciones_zj = []

for j in range(1, len(new_matrix_aux[0])):
    z = z_function
    for i in range(1, len(new_matrix_aux)):
        z = z.replace(new_matrix_aux[i][0], str(new_matrix_aux[i][j]))
    for k in range(1, len(new_matrix_aux[0])):
        z = z.replace(new_matrix_aux[0][k], '0')
    z = z.replace(z[0], z[0] + '_' + new_matrix_aux[0][j])
    lista_de_operaciones_zj.append(z)
    print(f'"{z}"')

variables, resultados = [], []

for operacion in lista_de_operaciones_zj:
    variable = operacion.split("Z_")[1].split(" = ")[0]
    variables.append(variable)
    expression = operacion.split("= ")[1]
    resultado = sp.sympify(expression)
    resultados.append(resultado)

output = [variables, resultados]
print(output)

This is the incorrect output that I am getting, where for some reason throughout the mathematical processes of the code the precision of the calculated values is lost (especially in those operations that involved fractions as can be seen below):

[['B', 'X1', 'X2', 'X3', 'X4', 'X5', 'U1', 'U2'], [8.0000000000000003, 1, 2, 2.000000000000000, 0, -0.9999999999999999, 0, 0.9999999999999999]]

This should be the correct output matrix where there was no loss of precision during the calculations.

[['B', 'X1', 'X2', 'X3', 'X4', 'X5', 'U1', 'U2'], [8, 1, 2, 2, 0, -1, 0, 1]]

I assume this problem is due to the nature of floating-point calculations in Python. Floating point numbers are internally represented in binary and can introduce small rounding errors.

Solution

So as others have said, you're not using Fraction at all – also, if I understood correctly what your new_matrix_aux stuff should do, here's a version that

uses proper Fractions
substitutes the desired variable values into a series of expressions and sympy-evaluates them:

from fractions import Fraction

import sympy as sp


def generate_funcs(function_template, variables):
    variable_names, variable_values = zip(*variables.items())
    for this_values in zip(*variable_values, strict=True):
        replacements = dict(zip(variable_names, this_values))
        func = function_template
        for var, value in replacements.items():
            func = func.replace(var, str(value) if isinstance(value, int) else f"({value})")
        yield func


out_names = ["B", "X1", "X2", "X3", "X4", "X5", "U1", "U2"]
vals = {
    "U1": [0] * 8,
    "U2": [0] * 8,
    "X1": [Fraction(4, 3), 1, 0, Fraction(2, 3), Fraction(-2, 3), Fraction(1, 3), Fraction(2, 3), Fraction(-1, 3)],
    "X2": [Fraction(10, 3), 0, 1, Fraction(2, 3), Fraction(1, 3), Fraction(-2, 3), Fraction(-1, 3), Fraction(2, 3)],
    "X3": [0] * 8,
    "X4": [0] * 8,
    "X5": [0] * 8,
}

z_function = "1 * X1 + 2 * X2 + 3 * X3 + 0 * X4 + 0 * X5 + M * U1 + M * U2"

results = {}

for out_name, operation in zip(out_names, generate_funcs(z_function, vals)):
    print(out_name, "=", operation)
    results[out_name] = sp.sympify(operation)

print(results)

The output is e.g.

B = 1 * (4/3) + 2 * (10/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X1 = 1 * 1 + 2 * 0 + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X2 = 1 * 0 + 2 * 1 + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X3 = 1 * (2/3) + 2 * (2/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X4 = 1 * (-2/3) + 2 * (1/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X5 = 1 * (1/3) + 2 * (-2/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
U1 = 1 * (2/3) + 2 * (-1/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
U2 = 1 * (-1/3) + 2 * (2/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
{'B': 8, 'X1': 1, 'X2': 2, 'X3': 2, 'X4': 0, 'X5': -1, 'U1': 0, 'U2': 1}