I am trying to find a solution for my equation with a command better than fminsearch.

I have a code as below. I am trying to modelling the relationship between solar radiation and temperature extremes and humidity extremes. I was using fminsearch, but the mean percentage error of my model is resulted as %9. I want to decrease this error. What can I do? Which alternatively commands will be suitable for me? I am very new in Matlab.

%My matlab code
  v = xlsread('I:\Yeni klasör\Tez 30 temmuz\Tez kilitlenme\Türkiye\Kahramanmaraş\matlab model.xlsx')
z=v(:,1)
y=v(:,2)
x=v(:,7)
u=v(:,5)
syms a, syms b, syms c
p(1) = a,  p(2) = b,  p(3) = c
z = @(p,v) v(:,2).*(p(1).*(1-exp(v(:,7).^p(2)))+p(3).*v(:,5))
P0 = [-0,3564; 0,0596; -0,4732];                                                                     % Choose Appropriate Initial Parameter Estimates
SSECF = @(p) sum((v(:,1) - z(p,v)).^2);                                             % Sum-Squared-Error Cost Function
options = optimset('MaxFunEvals',100000000);
options = optimset('MaxIter',100000000);
[abc, SSE] = fminsearch(SSECF, P0, options);                                                 % Estimate Parameters
a = abc(1)
b = abc(2)
c = abc(3)

[Walter Roberson]

If you post your xlsx file then I could try optimizing with the hybrid fitting algorithm that I wrote a few months ago.

[Can Ekici]

I added my xlsx file in attachments. "," is used for decimals in the file.

[Can Ekici]

Dear Walter, Could you look at my xlsx file? Thank you.

[Walter Roberson]

Could you confirm that your formula is

 z - y*(a*(1-exp(x^b))+c*u)

[Can Ekici]

Yes my formula is z = y*(a*(1-exp(x^b))+c*u)

Answer 1

[Star Strider]

If you have the Global Optimization Toolbox, use the patternsearch function to see if you can estimate a better set of parameters. If not, you can write your own genetic algorithm to do the same sort of search.

[Can Ekici]

I am very new in Matlab. And I don't know how to use optimization toolbox. Would you help me to create a algorithm on my code?

[Star Strider]

In its simplest form (without constraints), patternsearch works just like fminsearch:

[abc,fval] = patternsearch(SSECF, P0);

Also — Delete this line:

syms a, syms b, syms c

You do not need the Symbolic Math Toolbox, since you are not doing a symbolic computation.

[Can Ekici]

I was using the pattersearch command as you write. It gives the same results. An it gives a warning "Maximum number of iterations exceeded: increase options.MaxIter."

%My code
 v = xlsread('I:\Yeni klasör\Tez 30 temmuz\Tez kilitlenme\Türkiye\Kahramanmaraş\matlab model.xlsx')
z=v(:,1)
y=v(:,2)
x=v(:,7)
u=v(:,5)
syms a, syms b, syms c
p(1) = a,  p(2) = b,  p(3) = c
z = @(p,v) v(:,2).*(p(1).*(1-exp(v(:,7).^p(2)))+p(3).*v(:,5))
P0 = [-0,3564; 0,0596; -0,4732];                                                                     % Choose Appropriate Initial Parameter Estimates
SSECF = @(p) sum((v(:,1) - z(p,v)).^2);                                             % Sum-Squared-Error Cost Function
options = psoptimset('MaxIter',1000000);
[abc,fval] = patternsearch(SSECF, P0);
options = optimset('MaxFunEvals',1000000);
options = optimset('MaxIter',1000000);
[abc, SSE] = fminsearch(SSECF, P0, options);                                                 % Estimate Parameters
a = abc(1)
b = abc(2)
c = abc(3)

[Star Strider]

In your code, you are overwriting the ‘abc’ parameter vector patternsearch returns with the vector fminsearch returns.

Use this code instead:

v = xlsread('I:\Yeni klasör\Tez 30 temmuz\Tez kilitlenme\Türkiye\Kahramanmaraş\matlab model.xlsx')
z=v(:,1)
y=v(:,2)
x=v(:,7)
u=v(:,5)
% % % MAPPING: p(1) = a,  p(2) = b,  p(3) = c
z = @(p,v) v(:,2).*(p(1).*(1-exp(v(:,7).^p(2)))+p(3).*v(:,5))
P0 = [-0,3564; 0,0596; -0,4732];                                                                     % Choose Appropriate Initial Parameter Estimates
SSECF = @(p) sum((v(:,1) - z(p,v)).^2);                                             % Sum-Squared-Error Cost Function
options = psoptimset('MaxIter',1000000);
[abc,fval] = patternsearch(SSECF, P0);
% options = optimset('MaxFunEvals',1000000);
% options = optimset('MaxIter',1000000);
% [abc, SSE] = fminsearch(SSECF, P0, options);                                                 % Estimate Parameters
a = abc(1)
b = abc(2)
c = abc(3)

If patternsearch cannot find a better solution, the solution you have is the best you can get. We all have to live with the best fit to noisy data not being as good as we would like it to be.

[Can Ekici]

It gives the same warning again "Maximum number of iterations exceeded: increase options.MaxIter." Results are changed, but the mean percentage error is increased.

[Star Strider]

Continue to experiment with the patternsearch function. Increase 'MaxIter' and the appropriate other options.

It may take the function a while to discover the best parameter set, so give it that opportunity. Remember, it is doing an iterative search. It is probably the most efficient search algorithm available in the Toolbox.

[Can Ekici]

I am trying to increase maxiter. The warning is still same. How can I use this code in Toolbox?

SSECF = @(p) sum((v(:,1) - z(p,v)).^2);                                             % Sum-Squared-Error Cost Function
options = psoptimset('MaxIter',99999999999999);
[abc,fval] = patternsearch(SSECF, P0);
% options = optimset('MaxFunEvals',1000000);
% options = optimset('MaxIter',99999999999999);
% [abc, SSE] = fminsearch(SSECF, P0, options);                                                 % Estimate Parameters

[Alan Weiss]

You need to pass the options to patternsearch.

 [abc,fval] = patternsearch(SSECF, P0,[],[],[],[],[],[],[],options);

Alan Weiss

MATLAB mathematical toolbox documentation

[Can Ekici]

How can I get a better code to find my coefficients? Would you give me an example code?

[Star Strider]

We’ve already given you all the applicable MATLAB built-in functions and user-written File Exchange functions available. You know your data and the model you want to fit to it. You have to write the code to implement your model, then choose the solver best suited to estimate the parameters. The fminsearch function is a good place to begin, providing you don’t have too many parameters. The patternsearch function is much more robust in that it searches a greater parameter space.

You will likely have to settle for a less-than-perfect fit of your model function to your data. This is the usual outcome.

I would definitely not use a symbolic approach. The Symbolic Math Toolbox is not designed for iterative problems.

[Walter Roberson]

My code, the one I used to find the coefficients I gave below, is not publishable quality. I have continued to make significant changes to it since I used it to examine your coefficients. For example I improved bounds checking a fair bit.

I am in the middle of improving its memory behavior for handling more variables by breaking the search into smaller pieces and patching the results together, which gets especially important for curve-fitting (non-linear least-squares fitting) situations like you have.

I am outlining practices for replicating variables for bulk fitting for least-squared, but I am not nearly ready to automate that; currently this requires hand-manipulating the input functions.

The user interface to call my program is not very friendly, having a large number of parameters whose meaning is sometimes difficult to describe. I started cleaning that up a bit -- the program now accepts an "options" style of structure... whose meaning is not especially obvious to anyone not intimately familiar with the code.

It is, in short, a piece of research code. I would not put it out to beta test; I would hesitate to put it out to alpha test. It is, at best, at the "collaborating programmer in the same office" stage (if not for the fact that I work alone.)

It progresses if I am not busy answering Questions, and I am inspired, and I am feeling healthy enough, and I can't sleep, and I am not working on my other long term project. [That one is a program that primarily validates .fig files but is growing from there.]

That is just how things happen sometimes. Research ~= Development.

Answer 2

[Matt J]

Fminsearch is not guaranteed convergent for 3 variables, but I'm doubtful that the inaccuracy you see would be due to failure of the optimization. Still, you might try lsqcurvefit instead. Another good option might be fminspleas ( Download ), since your model function has a non-linear dependence on only one of your variables, b.

funlist={@(b,v) (1-exp(v(:,7).^b,  @(b,v) v(:,5)};

[b,ac]=fminspleas(funlist,P0(2),v,v(:,1) );

a=ac(1);
c=ac(2);

[Can Ekici]

I am trying to adapt my code, but it gives an error: "Undefined function or variable 'fminspleas'."

 %My code
   v = xlsread('I:\Yeni klasör\Tez 30 temmuz\Tez kilitlenme\Türkiye\Kahramanmaraş\matlab model.xlsx')
z=v(:,1)
y=v(:,2)
x=v(:,7)
u=v(:,5)
syms a, syms b, syms c
p(1) = a,  p(2) = b,  p(3) = c
z = @(p,v) v(:,2).*(p(1).*(1-exp(v(:,7).^p(2)))+p(3).*v(:,5))
P0 = [-0,3564; 0,0596; -0,4732];                                                                     % Choose Appropriate Initial Parameter Estimates
funlist={@(b,v) 1-exp(v(:,7).^b),  @(b,v) v(:,5)};
[b,ac]=fminspleas(funlist,P0(2),v,v(:,1) );
a=ac(1);
c=ac(2);

[Walter Roberson]

Matt J included a link to a location you can download fminspleas . Or you could use "Get Add-Ons" to install it.

Answer 3

[Walter Roberson]

The best location I could find was

[-0.24036462711882056186851741585996933281421661376953125, 0.13543162823623988710863841333775781095027923583984375, -0.00345339601881113934800371367828120128251612186431884765625]

where the sum-of-squared residue was 17101.81404290700083947740495204925537109375

[Can Ekici]

Thank you. I was changing the parameters as given. But the MPE results are seen same. Would you give the code to me?

[Walter Roberson]

Sorry, my search routine is unreleased code that is still under improvement. Generally speaking it combines grid searches and minimizations. It is for those times that brute force helps.

[Can Ekici]

I have too many datasets like this xlsx. I want to find a, b, c constants for each dataset with a more accurate code that is better than my code. Would you share an example code better than mine except of your unreleased code?

[Walter Roberson]

To optimize well, you need really good information about what the minimum and maximum ranges are for variables, and what the more "plausible" ranges are. For example you used

P0 = [-0,3564; 0,0596; -0,4732];

What that intentionally an array with 6 values, or did you mean

P0 = [-0.3564; 0.0596; -0.4732];

? If you did, then can we rely upon the pattern that the first value, a will always be negative, that b will always be positive, that c will always be negative? If so then that boundary information can be crucial: if you run searches over negative and positive ranges and you do not happen to have your granularity set fine enough, then you are very likely to get caught in a local minima near [4.42524074094453, -1.41837290431822, 0.0109500127989109] which has a residue of about 83603.5121110216 -- a region that you could eliminate from consideration entirely if you knew for sure that (for example) a must be negative.

Likewise, your formula includes (1-exp(v(:,7).^p(2))) . Do we know for sure that v(:.7) will be positive? If there is any possibility that it might be negative, then we can eliminate b from being negative, because b is continuous not integral, and negative to a non-integral power would be complex, and trying to find a mix of components that exactly balances complex terms out to be real-valued is more or less a waste of time (in such a case you are better off rewriting the formula into real and imaginary terms and force a solution where the overall imaginary term is 0.)

What else do we know about the range of v(:,7) ? If any v(:,7) can be greater than 1, then we need not check any p(2) that are larger than the point that exp(v(:,7).^p(2)) exceeds realmax and so we would not need to allow any p(2) larger than log(realmax)/log(max possible v(:,7)) . This is not just an academic discussion: in the searches I was doing, most of the residue results were coming out as infinite.

[Walter Roberson]

One technique that can be useful is to use just the initial few values of your vectors, just enough to have simultaneous equations instead of trying to do a fitting, and then to use the "exact" (to within numeric precision) results to be approximations of the eventual values. Provided that your initial values are somewhat representative of your other values, this can provide a fair bit of guidance as to what the expected values are to be.

When I do that with the first few values from your file, and using the formula z(k) - y(k)*(a*(1-exp(x(k)^b))+c*u(k)) that I extract from your variable names and code, I find that the model literally cannot fit those first few values -- or at least cannot if one restricts to real coefficients. Even extending to complex coefficients does not help. I can get values by taking the three (initial implied) equations two at a time, solving for a and c, substituting those in to the third equation, and solve()'ing for b, but back substitution shows the result is pure nonsense, producing overflows and undefined results.

When the initial values in the file cannot be fit by the model, it is time to ask whether the model is correct and whether the measurements are correct (including, for example, the possibility that the wrong columns are being extracted from the file.)

[Can Ekici]

I think, the measurements are correct. They are taken from Turkish Meteorology Agency. Parameters were measured with calibrated devices that are traceable to national standards. I am trying to find a new model which explains the relationship between these parameters. In the literature there is a model like "z = y*(a*(1-exp(x^b))". It can be said; my model is the modified version of that model. Would you give any recommendation to find different new models that explain the parameters well?

[Walter Roberson]

Correction to what I wrote above: if v(:,7) can be negative, that does not mean that p(2) cannot be negative, it means that p(2) can only be integers.

[Can Ekici]

v(:,7) cannot be negative. Because this value is the difference between daily maximum temperature and daily minimum temperature. It must be always positive. How can I enter these limitations to my matlab code?

[Can Ekici]

I have used Buckingham theorem in my dataset. And make the variables as dimensionless. Will it be better to solve this problem?

[Walter Roberson]

I am not familiar with the Buckingham theorem; I glanced at it just now, but I did not get enough "feel" for it to say whether it would help or not.

[Walter Roberson]

I noticed that you have some duplicate x and y that lead to different z.

The curvefitting toolbox cftool is able to quickly find parameters for the original model z = y*(a*(1-exp(x^b)), probably using fmincon with Trust-Region.

It might be worth exploring predicting with the classic model, subtracting that off from z, and then seeing how that residue varies with the variables, perhaps by plotting two at a time versus the residue. Sometimes visualization is much more useful than just trying random formulae.

[Can Ekici]

Thank you Mr. Robenson. I am trying to find the coefficients without code. You say and I think the best option is Curvefit toolbox to solve this equation.

[Star Strider]

‘I am trying to find the coefficients without code.’

Ouija Board?

[Walter Roberson]

Curvefit is not able to go beyond two independent variables and one dependent variable, but within that limitation it can be fast.