Comparing model to data to find the best fitting variable with scipy.optimize.minimize() [python]

Question

I need to compare a model to data to find the best fitting input. My data is surface reflectance of a body (= scattered light + thermal emission) at different wavelengths, and I simulate the reflectance by modeling the thermal emission part with Planck's law. I need to find the temperature at which the modelled emission (and thus reflectance) fits the observations.

To do it, I define Planck's law as:

h = 6.626070e-34 # the Planck constant
c = 2.997924e+8  # the speed of light in vacuum
k = 1.380649e-23 # the Boltzmann constant

def planck(T):  # Planck's law for black body radiation
    intensity = (2*h*c**2) / ( ((wvleng*1e-9)**5) * (np.exp(h*c/((wvleng*1e-9)*k*T)) - 1.0) )
    thermal = 0.95 * (intensity*np.pi*D**2)/(SSI*1e9)
    return thermal

which returns the thermal emission of the body at different wavelengths. Then, I defined the function that will qualify the fit (sum of squared difference between model and data at each wavelength):

def residuals(T): # 
  S = scat_light + 0.95*planck(T) # simulated reflectance
  return np.sum((reflectance - S)**2)

where T is the black body temperature. Finally, I use optimize.minimize to find which T leads to the best fit between S (simulated reflectance) and data reflectance:

x0 = 345.0  # initial guess, temperature is the only variable
res = scipy.optimize.minimize(residuals, x0, method='SLSQP') # optimization
fitted_temperature = res.x

Here's the problem: it gives me a temperature way below (~220K) than what it should be (around 340K), not fitting the data at all. Does anyone have an idea about why it doesn't converge properly? Thanks in advance.

Answer 1

According to the scipy.optimize.minimize documentation x0 should be, "Initial guess. Array of real elements of size (n,), where ‘n’ is the number of independent variables." So you're guessing that a temperature of T=345 will minimize the sum of squared differences. This is then returning a single best T for all T's in your data (assuming you have more than one sample) to minimize this function. Since apparently you know the answer should be around 340K I would imagine that there is an error in your formulas leading into or out of the planck() function somewhere.

It also wouldn't hurt to look at the other attributes of the res output to see if there can be additional clues garnered there.

Answer 2

Ok I had several problems in my code:

1) In my planck function, my thermal emission relative to solar spectral irradiance is thermal = 0.95*(intensity*4*np.pi)/(SSI*1e9) instead of thermal = 0.95 * (intensity*np.pi*D**2)/(SSI*1e9).

2) The biggest mistake, which gave me ~200K instead of ~300, was in my residuals function. The idea is to minimize the difference between the observed 'overall' reflectance spectra (= back-scattered light + thermal emission) and the simulated 'overall' spectra (correctly calculated as S = scat_light + 0.95*planck(T) in residuals. The error comes from the fact that this function returns the sum of squared differences between the BACK-SCATTERED light and the simulated OVERALL reflectance spectrum (=baskscaterred + thermal) : as the 'amount' of back-scattered light is much 'smaller' that the thermal emission at the studied wavelengths, it tried to fit a 'smaller' data (and thus explaining the low temperature fit).

I changed that and I now obtain a temperature of ~ 313 K that matches my observed overall reflectance curve very well.

Thank you all for your responses!