Fitting data in a non-linear function with constraints

Hello,

I need to fit some data (y_data = f(x1_data,x2_data)) in a certain function. I used 'objfcn', but I have to apply some boundary conditions (i.e. the curve must intersect a certain point, the second derivative must be positive, etc.).

Does anybody know if there is a way to do it with Matlab? (Well, I could artificially add the desired points to my data, but I guess this is not cientifically correct?)

I also checked lsqcurvefit, but, if I understood it well, you can only specify the ranges of the coefficients, but not the data.

In the following you find the case I need to solve, just in case my question was not clear.

Thank you in advance!

===================================================================================

I would like to fit some data in the following equation:

y_data = f(x1_data, x2_data)

Where f(x1, x2) = e^((a*x2^3 + b*x2^2 + c*x2 + d )*(-1/x1)) - z

and a,b,c,d and z are the coefficients to be determined.

There are many constraints:

1) 0 =< x2 =< 1

2) 0 =< y =< 1

(My data ranges from x2 = [0.3, 0.95] and y = [0.02,0.95])

3) x2 = 0 --> f = 0

(This is way I added the coefficient z. It is a kind of correction, since e^0 = 1, and there is otherwise no way to make it 0).

4) x2 = 1 --> f = 1

5) x2 = 1 --> f' = 1

6) x2 = 0 --> f'' > 0

To accomplish constraints 3, 4 and 5 I resolved the equations, so that dependent coefficients are obtained:

z = -e^(-d/x1)

c = -x1/(1+e^(-d/x1)) - 3*a -2*b

I tried to solve it as follows:

objfcn = @(b,x1,x2) exp(-1./x1.*( b(1).*x2.^3 + b(2).*x2.^2 + (-x1./(1+exp(b(3).*(-1./(x1))))-3*b(1)-2*b(2)).*x2 + b(3))) - exp(b(3).*(-1./(x1)));

[B,resnorm] = fminsearch(@(b) norm(y - objfcn(b,x1,x2)), [1;1;1]);

However, I don't know how to implement 1, 2 and 6.

*Here my questions:

- Is it possible to introduce constraints 3,4 and 5 with a Matlab function, instead of solving equations?

- Do you know if it is possible to apply constraints 1,2 and 6?*

Thanks!!!

[Matt J]

I don't understand (1) and (2). What is "y"? And since your only unknowns are a,b,c,d, & z, how can there be constraints on anything else?

[Vogel]

Hi Matt J,

thanks a lot for answering.

'y' is a variable, which is a function of x1 and x2. y = f(x1, x2). I have measurements of y,x1 and x2 that I need to correlate.

(1) and (2) mean that 'x2' and 'y' can just have values between 0 and 1.

You are right with the number of coefficients. Actually the function is much longer. I just tried to write a shorter example for readers.

The complete function is:

y = exp((-a./x1).*(b.*x2.^3 + c.*x2.^2 + d.*x2 + f)) - exp((p./x1).*(q.*x2.^3 + r.*x2.^2 + s.*x2 + t));

Where a,b,c,d,f,p,q,r,s,t are the coefficients to be determined.

[Matt J]

So you want y to be bounded between 0 and 1 for 0<=x1<=1 and 0<=x2<=1?

[Vogel]

Only for variables 'y' and 'x2'. But not for 'x1'.

[Vogel]

It's a 3D fitting

[Matt J]

But y is an exponential function depending on both x1 and x2. How can you bound an exponential over all x1? Exponentials are bounded only over finite intervals.

[Vogel]

Mmmm my mathematical knowledge is not so good. I will try to explain shortly the physical relationship:

x1 = temperature

x2 = mass fraction

y = thermodynamic property, which is a function of temperature and mass fraction

x2 and y can only get values between 0 and 1. It is impossible that they have other values. This should occur theoretically for any temperature, so for any x1. If mathematically I need to constrain x1, then I could give my temperature work range: 273.15 - 400 K.

Attached you see the resulting curves for two temperatures and the whole mass fraction range evaluated in the fit function that I managed to generate by using the procedure described in my first post (and adding artificially as measurements points [x2=0, y=0] and [x2=1,y=1] for temperatures 273:10:400).

It is quiet good, the problem is that for low values of x2 I get negative values of y (very small though). Moreover, the fit underestimates slightly the value of 'y'.

This is why I am looking for another way to fit the data.

Thanks again

[Matt J]

Based on what you've shown y is bounded as a function of both x1 and x2. Is there any combination of values x1,x2 such that y(x1,x2)>1 or y(x1,x2)<0?

[Vogel]

No, there is no possible combination of values x1,x2 that gives y(x1,x2)>1 or y(x1,x2)<0.

[Matt J]

The function

 y = exp((-a./x1).*(b.*x2.^3 + c.*x2.^2 + d.*x2 + f)) - exp((p./x1).*(q.*x2.^3 + r.*x2.^2 + s.*x2 + t));

doesn't look like a very physical model. If you don't already have a physical model, why is this difference of exponentials a compelling choice for a hypothesized model? Why not instead use some kind of sigmoid (e.g., arctan()) which is inherently bounded?

If it is a known physical model, why doesn't the existing physics theory that it came from tell you for what range of parameters the function is inherently bounded?

[Vogel]

It is a mixture of a physical and mathematical model.

The physics say:

1) ln(y)=(-1/x1).*a.*B(x2), where 'a' is a constant and B(x) introduces the influence of mass fraction x2. This is approximately true for some temperature ranges, which can differ for every substance.

2) y = f(x1, x2=0)=0 and y = f(x1, x2=1)=1

3) 0<=y(x1,x2)<=1

4) y'= f'(x1, x2=1) = 1

5) At x2=0, the slope of the curve must be a stright line with a positive slope.

I need to find B(x2). However, there are not physical models to determine it, so here I just looked for mathematical functions that permits implementing all the boundary conditions. I tried many: a*x2^b, polynoms of second, third and fourth order, ....etc. The best one I found was this: B(x2) = b.*x2.^3 + c.*x2.^2 + d.*x2 + f

So that ln(y)=(-a/x1).*(b.*x2.^3 + c.*x2.^2 + d.*x2 + f).

Therefore: y = exp((-a/x1).*(b.*x2.^3 + c.*x2.^2 + d.*x2 + f))

It is observed that, due to the exponential, 'y' can never be eaxct zero. This is why I had to add a correction factor z, which needes to be a function of temperature and mass fraction. I just applied the same form as the principal function.

z = exp((-p/x1).*(q.*x2.^3 + r.*x2.^2 + s.*x2 + t))

Do you know any other function that could pass better? And, in any case, do you know how to apply constraints 2 to 5 to a fit function?

Thanks

[Matt J]

Well (1) and (4) look impossible to satisfy simultaneously. I assume the f'() is the derivative with respect to x2, not x1. If so, then using (1) we can compute the form of f'() explicitly

y'(x1,x2=1) = (a/x1)^2*B(1)*B'(1)

There is no way the right hand side can equal 1 independently of x1, as demanded by constraint (4).

[Vogel]

Yes, f' is the derivative respect x2 and I calculated it my self and applied the constraint. However, for complicated functions, it is an arduous task. I was wondering if there was a way for Matlab to do it during the fitting process.

Answer 1

[Matt J]

I need to find B(x2). However, there are not physical models to determine it, so here I just looked for mathematical functions that permits implementing all the boundary conditions.

If you aren't committed to a particular model, it might be useful for you to look at SLM. It allows you to do spline fitting subject to various boundary conditions and montononicity conditions both on the function and on slope.

[Vogel]

This is very interesting. Thank you.