Constraining fit via function value instead of coefficient values

I am modelling some electric parameters of a battery as a function of the state of charge of the battery. This state of charge is always a value between 0 and 1. The electric parameters such as the internal resistance can only be positive (negative resistance is not physically possible).

I'd like to fit a function through the measured datapoints and make sure that the functionvalues of this fit for x in [0 1] are all positive.

Is there a way to constrain the fit function to do this?

The following code gives an example of a dataset that is fitted correctly and one that results in an interpolating function that becomes negative for x>0.93 which is not physically possible.

x1 = [0.9573 0.8974 0.8286 0.7748 0.7190 0.6586 0.5952 0.5359 0.4761 0.4161 0.3553 0.2969 0.2345 0.1812 0.1160 0.0606 0.0000];
y1 = [1.0015 0.5321 0.7086 0.9221 0.3032 0.2380 0.3100 0.3532 0.3183 0.3900 1.5636 0.2052 0.1584 0.1523 0.1993 0.1433 0.1107]*10^3;
x2 = [0.9133 0.7934 0.5838 0.4764 0.3748 0.1682 0.0630 0.0000];
y2 = [0.0197 0.0541 0.0232 0.0116 0.0965 0.0454 0.0355 0.0182]*10^4;

figure
subplot(211)
hold on
scatter(x1, y1)
plot(fitting_function(x1, y1))

subplot(212)
hold on
scatter(x2, y2)
plot(fitting_function(x2, y2))


function [coefs] = fitting_function(x, y)
    [coefs, ~] = fit(x', y', 'poly5', 'Robust', 'Bisquare');
end

This code uses a simple polynomial fit but the actual code also uses other models.

Thanks

[Torsten]

Are the data values you try to fit also in the range [0 1] ?

[Simon]

Yes the x-values are all within [0 1]

[the cyclist]

If possible, it would probably be helpful for you to post

• the code that does the fit
• an example of a dataset the fits well by default
• an example of a dataset the fits poorly by default

[Torsten]

@Simon

I thought you wanted the y-values to be in [0 1] ?

That's why I asked: Are your original y-data (those that are to be fitted) already in [0 1] ?

[Simon]

@Torsten@the cyclist I rewrote my question to make it a bit clearer and realised that what I wanted has not really to do with the initial guess for the coefficients but rather constraining the function values in general within a certain range. Hope this makes it clearer.

[the cyclist]

It is probably important here to understand your goal for this fit. What are you going to do with the result?

High-order polynomials are almost never the best way to perform a fit, unless you know that the underlying physical process leads to that specific functional form. They tend to lead to wildly inaccurate fits for some points, especially (but not exclusively) when extrapolating outside the data range. Is there are reason you chose a fifth-degree polynomial?

[Simon]

@the cyclist I'm trying to show that some models which are used in literature cannot be used in generally for a certain testing procedure and only in selected cases. One of these models closely resembles this fifth-degree polynomial. For my purpose it's no problem that the fit is bad (that's exactly what I want to show) but it should at least by physically meaningfull.

Answer 1

[John D'Errico]

No. You cannot use fit to perform such a fit, where you place a constraint on the function values.

And, yes, a polynomial is a bad thing to use for such a fit, but you don't seem to care.

Regardless, you cannot put a constraint that the MAXIMUM value of the polynomial (or minimum) be any specific value. The problem is, the maximum is a rather nonlinear function. And you don't know where the maximum lies either. It may be at some interior point on the curve, or it may be at one of the endpoints of the interval. The thing is, the maximum is quite difficult to pin down. It will lie at one of the roots of the first derivative of the polynomial, OR at an endpoint. And that means you will not be able to use it for a constraint, even for fmincon, since that maximum will not even be a differentiable function of the parameters.

So the best solver here is problemy a tool like GA, which seriously does not care about crap like differentiability. Throw stiuff at it, yeah, right, it can handle anything.

Having said all of that, there are some tricks one can employ. The simplest one is to recognize the difference between a necessary and sufficient condition on the maximum. For example, a NECESSARY condition for the maximum is that the maximum of a function will be at least as large as the maximum of the function, sampled at many locations. As such, we can constrain the polyno,mial to NEVER exceed 1 at the points

xtest = linspace(0,1,100);

But that is not truly sufficient to define the maximum, since the max may lie between one of those points. That is, the max may be slighttly larger then the value at any of those 100 points. I'll show an example.

x = rand(100,1)*.6 + .2;
y = sin(x*pi) + randn(size(x))/4; % a noisy portion of a sine wave

Now, we know that the true maximum of the function can never be greater than 1 here, since this is a sine wave, right? But the computer has no such knowledge. And if we desire to fit the curve with a polynomial, then we cannot tell fit to obey something it does not understand.

poly5 = fit(x,y,'poly5')

poly5 = 
     Linear model Poly5:
     poly5(x) = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6
     Coefficients (with 95% confidence bounds):
       p1 =      -396.6  (-980.3, 187.1)
       p2 =        1003  (-462.2, 2468)
       p3 =      -982.1  (-2403, 438.6)
       p4 =       456.5  (-206.6, 1120)
       p5 =      -97.71  (-246.2, 50.8)
       p6 =       8.325  (-4.42, 21.07)

plot(x,y,'o')
hold on
plot(poly5)
hold off

Honestly, this is a bit of a mathematical obscenity to use a 5th degree polynomial on this complete crap for data. So sue me.

What is the maximum of that function on the interval? I can use fminbnd to give me that information, or I could have differentiated the polynomial, find the roots, then eveluate the polynomial at each root location, plus the interval endpoints.

format long g
[xmax,fmax] = fminbnd(@(x) -poly5(x),min(x),max(x));
xmax

xmax = 
         0.484298335981994

fmax = -fmax

fmax = 
          1.12124540162064

As you can see, the fit does extend above 1. But again, I cannot use fit to constrain the polynomial.

I could use lsqlin to do the fit, to constrain the polynomial at a specific set of points. Again, this is only a necessary thing. It is NOT sufficient to insure the polynomial ALWAYS lies below 1. So I might do this:

xcon = linspace(min(x),max(x),100)';
V = xcon.^(5:-1:0);

A = x.^(5:-1:0);
poly5con = lsqlin(A,y,V,ones(size(xcon)))

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

poly5con = 6×1

-309.118860734151 717.050526979835 -638.026054657531 266.190330572932 -49.6524178368595 3.87581914363508

So that polynomial was constrained to NEVER exceed 1, but only at the 100 test points. It is NOT a sufficient condition that the polynomial does not exceed 1. In fact, the maximum of that polynomial lies at (I'll use roots this time instead of fminbnd)

poly5conder = polyder(poly5con);
xpotential = [min(x),roots(poly5conder)',max(x)];
xpotential(imag(xpotential) ~= 0) = [];
polymax = max(polyval(poly5con,xpotential))

polymax = 
          1.00002258830354

So even though we constrained the polynomial to NOT exceed 1, it did anyway!

The point is, by setting such a constraint at 100 specific points, we missed the exact point we needed. That was only a necessary condition on the maximum. It was not sufficient to insure the polynomila never exceeds 1. This unfortunately, is about the best you can do, if you are using a tool like lsqlin. Again, fit has no chance to constrain anything like that. Ar the two polynomials different? They both were fit to the same set of data.

[coeffvalues(poly5)',poly5con]

ans = 6×2

1.0e+00 * -396.604966647345 -309.118860734151 1002.87384163545 717.050526979835 -982.055587920609 -638.026054657531 456.474479569381 266.190330572932 -97.7068727744401 -49.6524178368595 8.32507269015535 3.87581914363508

So you should see the coefficients of the two polynomials are actually seriously different, even though the maximum value was not that far above 1 even for the first fit.

fplot(@(x) poly5(x),[min(x),max(x)],'r')

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on
plot(x,y,'go')
fplot(@(x) polyval(poly5con,x),[min(x),max(x)],'b')
hold off

You can see the difference here. The red curve had no constraint. The blue curve had only the necessary constraints, at 100 points. Have I used more points, say 1000, then the

xcon = linspace(min(x),max(x),1000)';
V = xcon.^(5:-1:0);
poly5con1000 = lsqlin(A,y,V,ones(size(xcon)))

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

poly5con1000 = 6×1

-313.620756447606 728.200012493809 -648.577300793215 270.931498896521 -50.6598802609233 3.95681479356138

[coeffvalues(poly5)',poly5con,poly5con1000]

ans = 6×3

1.0e+00 * -396.604966647345 -309.118860734151 -313.620756447606 1002.87384163545 717.050526979835 728.200012493809 -982.055587920609 -638.026054657531 -648.577300793215 456.474479569381 266.190330572932 270.931498896521 -97.7068727744401 -49.6524178368595 -50.6598802609233 8.32507269015535 3.87581914363508 3.95681479356138

The last curve is refined yet more, with the coefficients being not terribly different.

[xmax,fmax] = fminbnd(@(x) -polyval(poly5con1000,x),min(x),max(x));
xmax

xmax = 
         0.484790757503588

fmax = -fmax

fmax = 
          1.00000023590079

So, with 100 points for the constraints, the max value was 1.000022..., with 1000 points, I got two more zeros in there, even though the max STILL exceeded 1. Is the tiny ecess over 1 a problem? Only you know that. Such is life. I hope that at least somebody reads this far. :)

[Simon]

Thank you for this very elaborated answer! This does help me a lot and I hope others as well.

And no, the tiny excess over 1 is not a problem for me ;)