There are mutiple approaches possible, depending on the type of known relations and expression to evaluate. If they are linear over the values of $U()$ you can employ a linear algebra over the vectors of used $U()$ values. For a given example, we can work with vectors of coefficients of $U(10),U(16),U(20)$, and substitution of $U(16)$ into $EU_1$ can be performed by matrix multiplication:
eu_1 = vector([0.1,0.4,0.5])
eu_2 = Matrix([ [1,0.2,0], [0, 0, 0], [0, 0.8, 1] ]) * eu_1
Here we get eu_2 equal (0.180000000000000, 0.000000000000000, 0.820000000000000) with the coefficient of $U(16)$ (i.e., the second component) nullified.
If relations and expressions are (not necessarily linear) polynomials, a possible approach is to represent each used $U()$ value as a polynomial variable, and each known relation as a polynomial over those variables. Then defined a polynomial ring with the lexicographical term order (with variables to be substituted coming first), define an ideal generated by the relations, and perform reduction of a polynomial of interest with respect to this ideal:
K.<u16,u10,u20> = PolynomialRing(RR, order='lex')
J = K.ideal( [u16 - .2*u10 - .8*u20] )
eu_1 = 0.1*u10 + 0.4*u16 + 0.5*u20
eu_2 = eu_1.reduce(J)
Here we get eu_2 equal 0.180000000000000*u10 + 0.820000000000000*u20, that is, u16 got substituted.
