Regression for function approximation I have a program for heat exchanger calculations which uses correlations that are complex and highly non-linear. I need to come up with an approximation of the function using regression. 
The function takes two inputs : inlet air temperature $T_{a,i}$ and valve control signal $u$. The output is $\Delta T_a$. I have plotted the output values from the program for a range of input values. How can the function be approximated with regression? 
EDIT : The data is generated from the following python program which implements the heat exchanger model in [2] which is also used in the Building Simulation Library 3. I need to find aproximations for the efficiency and outlet-inlet temperature difference for the air and water streams. The efficiency one is simple and can be approximated by an exponential function but I am struggling with the others.
# Air-water finned heat exchanger model
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pylab import meshgrid,show,figure
exp = np.exp
rand = np.random.rand

def fcumodel(T_a_in, T_w_in, fs, val):
    # Nominal parameters
    T_a_in_0 = 20.0  # Nominal air inlet temperature
    T_w_in_0 = 40.0  # Nominal water inlet temperature
    T_a_out_0 = 27.0 # Nominal air outlet temperature
    T_w_out_0 = 32.0 # Nominal water outlet temperature

    c_a = 1005.0     # specific heat capacity of air
    c_w = 4180.0     # specific heat capacity of water

    rh = 0.5         # ratio of heat transfers (at nominal conditions)

    mdota_0 = 5.0    # Nominal air mass flow rate
    mdotw_0 = 1.0    # Nominal water mass flow rate

    #NTU0 = calculate_nominal(T_a_in_0, T_w_in_0, T_a_out_0, T_w_out_0, mdota_0, mdotw_0)
    #print NTU0

    NTU0 = 0.747     # precomputed externally for above parameters

    UA_0 = NTU0 * min(c_a * mdota_0, c_w * mdotw_0)   # Nominal conductance

    hw_0 = UA_0 * (rh+1)/rh   # water side heat transfer

    ha_0 = rh * hw_0          # air side heat transfer

    # operating values
    mdota = fs * mdota_0      # air mass flow rate  
    mdotw = val * mdotw_0     # water mass flow rate

    Cdota = c_a * mdota       # air capacitance rate
    Cdotw = c_w * mdotw       # water capacitance rate

    x_a = 1 + 4.769 * 1e-3 * (T_a_in - T_a_in_0)   # air side heat transfer temperature variation

    ha = ha_0 * x_a * (mdota/mdota_0)**0.7         # air side heat transfer

    s_w = 0.014/(1 + 0.014 * T_w_in)               # water side heat transfer temperature sensitivity

    x_w = 1 + s_w * (T_w_in - T_w_in_0)            # water side heat transfer temperature dependence

    hw = hw_0 * x_w * (mdotw/mdotw_0)**0.85        # water side heat transfer

    UA = 1/(1/hw + 1/ha)         # overall heat exchanger conductance             

    Cdotmin = min(Cdota, Cdotw) 
    Cdotmax = max(Cdota, Cdotw)
    Z = Cdotmin/Cdotmax

    if Cdota < Cdotw:
        small = 'air'
    else :
        small = 'water'

    NTU = UA/Cdotmin

    eff = 1 - exp((exp(-Z * (NTU ** 0.78) ) - 1) * (NTU ** 0.22) / Z)    # Cross flow heat exchanger eff

    Qdot = eff * Cdotmin * (T_w_in - T_a_in)

    T_a_out = T_a_in + Qdot/Cdota
    T_w_out = T_w_in - Qdot/Cdotw

    print "eff", eff 
    #print "det:", eff*(Z+1)

    return Qdot, T_a_out, T_w_out, 1/UA, NTU, eff, small


R_a = np.arange(15,21,0.5)
R_u = np.arange(0.1,1,0.05)

N_a = R_a.shape[0]
N_v = R_u.shape[0]

EFF = np.zeros((N_a,N_v),dtype = 'float64')
Ta = np.zeros((N_a,N_v), dtype = 'float64')
Tw = np.zeros((N_a,N_v), dtype = 'float64')


for i in range(N_a):
    for j in range(N_v):
        T_a_in = R_a[i]
        T_w_in = 35

        fs = R_u[j]
        val = R_u[j]

        Qdot, T_a_out, T_w_out, R, NTU, eff, small = fcumodel(T_a_in, T_w_in, fs, val)

        EFF[i][j] = eff
        Ta[i][j] = T_a_out - T_a_in
        Tw[i][j] = T_w_out - T_w_in


[x, y] = meshgrid(R_u, R_a)

fig1 = figure(1)
ax = Axes3D(fig1)
ax.plot_surface(x,y,EFF)


fig2 = figure(2)
ax = Axes3D(fig2)
ax.plot_surface(x,y,Ta)

fig3 = figure(3)
ax = Axes3D(fig3)
ax.plot_surface(x,y,Tw)

show()


[2]: Wetter, M. (1999). Simulation model finned water-air-coil withoutcondensation (No. LBNL--42355). Ernest Orlando Lawrence Berkeley National Laboratory, Berkeley, CA (US).
 A: Your problem is much more straightforward than a typical regression problem because 1) You know the actual ground truth, and can generate as many data points as you'd like, 2) There's no noise, and 3) The function is very well behaved (smooth and low dimensional).
Simple interpolation should work well in this setting (also see here). For example, the surface plot you showed was almost certainly generated using interpolation. In this approach, you'd evaluate the function at a fixed set of sample points. The function is then approximated as a set of local pieces. Each piece is a simple function defined over the space between neighboring sample points, and runs through the true function values at these points.
The type of local function used determines the type of interpolation. For example, linear and cubic spline interpolation are popular choices. Cubic splines require more computation, but produce smoother approximations. Since you're using Python, check out scipy.interpolate.
Interpolating using a greater number of sample points gives a more accurate approximation at the expense of more computation. In some cases it's possible to balance between accuracy and computational expense by allocating sample points adaptively, such that they're more densely distributed in regions where the function is more complicated. In the simplest case, sample points can be chosen on a regular grid. Extrapolating outside the range of the sample points may not work well, so they should be chosen to span the largest range you anticipate needing.
There are many fancier regression methods that would also give good results (e.g. Gaussian process regression, kernel methods, tree/forest-based methods, various basis expansions, etc.). But, because of the simplicity of your function and the lack of noise, I don't think the added complexity would buy you anything over simple interpolation.
Another alternative would be to come up with a simple, parametric form for the approximation, then fit the parameters. This would have the advantages of being highly computationally efficient and easily interpretable. But, it would require that such a form exists (with high enough accuracy), and that you can find it.
