Sympy无法将代入值后的克罗内克δ不等式函数求值为数值

我已尝试代入所有对应值评估函数，但SymPy无法将其进一步求值为数值。该函数是矩阵中的一个元素，所有类似元素均无法完成求值。我试过simplify()、doit()等多种组合，都没有效果。

以下是复现该问题的代码：

from sympy import summation, Min, IndexedBase, symbols, pretty_print, Matrix, diff, ImmutableDenseMatrix

K1, K2 = symbols("K1 K2")
N_T, N_f, i, j = symbols("N_T N_f i j", integer=True, positive=True)

s_x = IndexedBase("s_x")
s_y = IndexedBase("s_y")
s_z = IndexedBase("s_z")
s_star_x = IndexedBase("s_star_x")
s_star_y = IndexedBase("s_star_y")
s_star_z = IndexedBase("s_star_z")
w = IndexedBase("w")
h = IndexedBase("h")


def norm2(vec):
    if isinstance(vec, (list, Matrix)):
        return sum(comp**2 for comp in vec)
    else:
        return vec**2


s_diff_sum = summation(
    summation(
        K1
        * norm2(
            [
                s_x[i] - s_x[j] - (s_star_x[i] - s_star_x[j]),
                s_y[i] - s_y[j] - (s_star_y[i] - s_star_y[j]),
                s_z[i] - s_z[j] - (s_star_z[i] - s_star_z[j]),
            ]
        ),
        (j, i + 1, Min(N_f, i + 3)),
    ),
    (i, 1, N_f),
)


matrix_elements = []
N_f = 8
Nt_value = 4

for t in range(N_f):
    row1 = [s_x[t]]
    row2 = [s_y[t]]
    row3 = [s_z[t]]
    matrix_elements.append(row1)
    matrix_elements.append(row2)
    matrix_elements.append(row3)

S = Matrix(matrix_elements)


matrix_elements = []


for t in range(Nt_value):
    row1 = [h[t]]
    row2 = [w[t, 1]]
    row3 = [w[t, 2]]
    row4 = [w[t, 3]]

    matrix_elements.append(row1)
    matrix_elements.append(row2)
    matrix_elements.append(row3)
    matrix_elements.append(row4)


matrix_elements.append(S)

U = Matrix(matrix_elements)

pretty_print(U)


cost_diff_U_2 = diff(s_diff_sum, U, 2)

N_T, N_f, i, j = symbols("N_T N_f i j", integer=True, positive=True)

cost_diff_U_2 = cost_diff_U_2.subs({N_f: 8, K1: 0.1})

cost_diff_U_2_reshaped = cost_diff_U_2.reshape(40 * 1 * 40 * 1)

cost_diff_U_2_reshaped = cost_diff_U_2_reshaped.reshape(40, 40)

cost_diff_U_2_reshaped = ImmutableDenseMatrix(
    cost_diff_U_2_reshaped.tolist()
)

pretty_print(cost_diff_U_2_reshaped.doit(), num_columns=200)

内容的提问来源于stack exchange，提问作者Rhyn155