Best optimization protocol for minimization

Question

I need to do a simultaneous fit of the following spectral data. I need to fit skewed gaussians to a spectrum A, take the first and second derivative of this fitted skewed_gaussian spectrum and fit that two spectrums B and C each with their own coefficients. I would like to determine the residuals between the actual and fitted data and use this as a Cost function to minimize. I am trying to find the optimal way to do this fit. I thought of using minimize from scipy but I know it can be a bit slow. What is the best way to approach this?

Let me clarify a couple of things. Spectrum A is an absorption spectrum. Spectrum B and C are Stark Absorption Spectra (Spectra obtained once an electric field is passed through it). The lineshape of the Stark Spectra are dependent on the normal Absorption Spectra.

At the moment I have a Find_Residue function which takes an initial guess and builds the Absorption Spectrum for skewed gaussians, computes the first and second order derivatives and uses initial guesses of the coefficients to fit the Stark Spectra as follows:

SA = c1A + c2(derv A) + c3*(derv2 A) SA2 = d1A + d2(derv A) + d3*(derv2 A)

The find_residue function then determines the difference between these fitted functions and actual results and computes the residue which I attempt to minimize using scipy.optimize.minimize.

Since the fitted Skewed Gaussians yield physically relevant properties, I would like to avoid inputting bounds or guesses if at all possible to determine a more accurate fitting of the data. At the moment the fitting does not yield good results and I believe it's due to my bounds which is based on visual subjectivity