It might help to “stack” your problem into one big tall y vector, instead of casting it as y_ij.
I think then your data y would have a shape of (n x 1) and look like
y = (
(dut_1, T_1),
(dut_1, T_2),
...
(dut_1, T_n),
(dut_2, T_1),
(dut_2, T_2),
...
(dut_2, T_n),
...
)
Your X variable would be set up similarly, (n x 2). Then you have a 2d GP with inputs time and condition. Depending on how many data points you have this could be slow. If it’s possible for you to use a covariance k(t_i, S_i) = k1(t_i, t_i*) * k2(S_i,S_i*) (multiply instead of add), then it’s separable and you can use the Kronecker implementation, which should greatly speed things up for you. Please let me know if I didn’t quite understand your question!