Constraining fit via function value instead of coefficient values

2 次查看（过去 30 天）
Simon
Simon 2023-4-13
评论： Simon 2023-4-13
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
  7 个评论
显示 5更早的评论
the cyclist
the cyclist 2023-4-13
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
Simon 2023-4-13
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2023-4-13
编辑：John D'Errico 2023-4-13
Ran in:
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. :)
  1 个评论
Simon
Simon 2023-4-13
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 ;)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by