How can I tune an existing model to fit two correlations with multiple (>2) independent variables?

function(T_out1, P_out1, T_out2,P_out2) = thermal(T_in1, P_in1, T_in2, P_in2)

k1 = fn(T_in1, P_in1, T_in2,P_in2) and of the form a*x1^b*x2^c k2 = fn(T_in1, P_in1, T_in2,P_in2) and of the form d*x1^e*x2^f*x3^g

T_out1 = f(T_in1, P_in1, T_in2,P_in2, k1.k2) T_out2 = f(T_in1, P_in1, T_in2,P_in2, k1,k2) end

I have 20 data points (T_outs, P_outs for corresponding T_ins, and P_ins) and I want to use these to calculate the constants(a,b,c,d,e,f,g) in the equations for k1, and k2. I tried using curve fitting for k1 but k2 has more than two independent variables. Another concern here is that both of these equations need to be satisfied simultaneously for the given outputs.