My function creates a high oscillating unusual plot

7 次查看（过去 30 天）
Wathsala Karunarathne
回答： Christopher Creutzig 2021-8-6
I need to plot and find the optimum point of the function below. But it creates a weird plot which I have inserted. First, I thought this is because the function has coefficients less than matlab eps. But setting them to zero did not change the plot. Also, the optimization alogorithms (fmincon,ga,particle swarm,PSO) gives wrong answers. Does anyone have an explaination for this behaviour?
Thank you very much.
f =17796509464316297.363895748*exp(-11.0*t1) - 17796509464316296.574074074074074*exp(-10.0*t1) - 0.68982167392592592592592592592593*exp(-1.0*t1) + 1.2222222222222222222222222222222*t1*exp(-1.0*t1) + 6174450735868760.1111111111111111*t1*exp(-10.0*t1) + 11622058728447497.43895748*t1*exp(-11.0*t1) - 982200909334137.5*t1^2*exp(-10.0*t1) + 93586137679683.333333333333333333*t1^3*exp(-10.0*t1) - 5859613763979.1666666666666666667*t1^4*exp(-10.0*t1) + 252901537750.0*t1^5*exp(-10.0*t1) - 8145447708.3333333333333333333333*t1^6*exp(-10.0*t1) + 219516865.07936507936507936507937*t1^7*exp(-10.0*t1) - 2683903.7698412698412698412698413*t1^8*exp(-10.0*t1) - 222938.71252204585537918871252205*t1^9*exp(-10.0*t1) + 12125.220458553791887125220458554*t1^10*exp(-10.0*t1) + 3706004905623298.2447874*t1^2*exp(-11.0*t1) + 767474314438508.23262466666666667*t1^3*exp(-11.0*t1) + 115820035783760.539895*t1^4*exp(-11.0*t1) + 13549259810802.40979*t1^5*exp(-11.0*t1) + 1276017179970.1820111111111111111*t1^6*exp(-11.0*t1) + 99090055157.583646825396825396825*t1^7*exp(-11.0*t1) + 6431882521.8889831349206349206349*t1^8*exp(-11.0*t1) + 349879383.20745425485008818342152*t1^9*exp(-11.0*t1) + 15772035.50003169091710758377425*t1^10*exp(-11.0*t1) + 573815.32070005611672278338945006*t1^11*exp(-11.0*t1) + 16540.373890295982309871198760088*t1^12*exp(-11.0*t1) + 427.06496366609213831436053658276*t1^13*exp(-11.0*t1) + 14.861012365197633054775911918769*t1^14*exp(-11.0*t1) + 0.44769060580072484834389596294358*t1^15*exp(-11.0*t1) + 0.0045334735715336740469015601290734*t1^16*exp(-11.0*t1) + 0.0044089949942699028444437719258759*t1^17*exp(-11.0*t1) + 0.00072758477778833722136969361785214*t1^18*exp(-11.0*t1) + 0.000077871208689571203517652792619439*t1^19*exp(-11.0*t1) + 0.000005672202354891583523257115371024*t1^20*exp(-11.0*t1) + 0.0000002583191156593700955608594990967*t1^21*exp(-11.0*t1) + 0.0000000065841638862926674918780324157471*t1^22*exp(-11.0*t1) + 0.00000000007601031748692706747794486283681*t1^23*exp(-11.0*t1) + 0.00000000000029107980533995897383819762284709*t1^24*exp(-11.0*t1) - 0.00000000000000031590056393483919641354780774103*t1^25*exp(-11.0*t1) - 2.0

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

采纳的回答

Walter Roberson
Walter Roberson 2021-8-3
Ran in:
The summary of the below investigation is that what you are observing appears to be a bug or limitation in MATLAB.
It might be this bug: https://www.mathworks.com/support/bugreports/2422167?ref=ts_R2020b_Update_6 but that is supposed to be fixed by R2021a Update 1.
I will refer this over to Mathworks as it looks like a bug I reported earlier might not be fixed after all.
syms t1
f =17796509464316297.363895748*exp(-11.0*t1) - 17796509464316296.574074074074074*exp(-10.0*t1) - 0.68982167392592592592592592592593*exp(-1.0*t1) + 1.2222222222222222222222222222222*t1*exp(-1.0*t1) + 6174450735868760.1111111111111111*t1*exp(-10.0*t1) + 11622058728447497.43895748*t1*exp(-11.0*t1) - 982200909334137.5*t1^2*exp(-10.0*t1) + 93586137679683.333333333333333333*t1^3*exp(-10.0*t1) - 5859613763979.1666666666666666667*t1^4*exp(-10.0*t1) + 252901537750.0*t1^5*exp(-10.0*t1) - 8145447708.3333333333333333333333*t1^6*exp(-10.0*t1) + 219516865.07936507936507936507937*t1^7*exp(-10.0*t1) - 2683903.7698412698412698412698413*t1^8*exp(-10.0*t1) - 222938.71252204585537918871252205*t1^9*exp(-10.0*t1) + 12125.220458553791887125220458554*t1^10*exp(-10.0*t1) + 3706004905623298.2447874*t1^2*exp(-11.0*t1) + 767474314438508.23262466666666667*t1^3*exp(-11.0*t1) + 115820035783760.539895*t1^4*exp(-11.0*t1) + 13549259810802.40979*t1^5*exp(-11.0*t1) + 1276017179970.1820111111111111111*t1^6*exp(-11.0*t1) + 99090055157.583646825396825396825*t1^7*exp(-11.0*t1) + 6431882521.8889831349206349206349*t1^8*exp(-11.0*t1) + 349879383.20745425485008818342152*t1^9*exp(-11.0*t1) + 15772035.50003169091710758377425*t1^10*exp(-11.0*t1) + 573815.32070005611672278338945006*t1^11*exp(-11.0*t1) + 16540.373890295982309871198760088*t1^12*exp(-11.0*t1) + 427.06496366609213831436053658276*t1^13*exp(-11.0*t1) + 14.861012365197633054775911918769*t1^14*exp(-11.0*t1) + 0.44769060580072484834389596294358*t1^15*exp(-11.0*t1) + 0.0045334735715336740469015601290734*t1^16*exp(-11.0*t1) + 0.0044089949942699028444437719258759*t1^17*exp(-11.0*t1) + 0.00072758477778833722136969361785214*t1^18*exp(-11.0*t1) + 0.000077871208689571203517652792619439*t1^19*exp(-11.0*t1) + 0.000005672202354891583523257115371024*t1^20*exp(-11.0*t1) + 0.0000002583191156593700955608594990967*t1^21*exp(-11.0*t1) + 0.0000000065841638862926674918780324157471*t1^22*exp(-11.0*t1) + 0.00000000007601031748692706747794486283681*t1^23*exp(-11.0*t1) + 0.00000000000029107980533995897383819762284709*t1^24*exp(-11.0*t1) - 0.00000000000000031590056393483919641354780774103*t1^25*exp(-11.0*t1) - 2.0
f = 
f2 = collect(f, t1)
f2 = 
fplot(f2, [0 1e-15])
df1 = diff(f2, t1)
df1 = 
fplot(df1, [0 1e-15])
fplot(df1, [0 3e-17]); ylim auto
T1 = linspace(0, 10, 1000).';
f2 = subs(f, t1, T1);
vpa(f2(1:10))
ans = 
plot(T1, double(f2))
  2 个评论
Walter Roberson
Walter Roberson 2021-8-3
I updated my case with Mathworks.
The summary of the workaround I used here is to choose specific values and subs() them in and then double() the results.
The alternative workaround would be to wrap the function like
F = @(T1) double(subs(f, t1, T1))

请先登录，再进行评论。

更多回答（1 个）

Christopher Creutzig
Ran in:
What you are observing is a limitation of computing with floating point numbers. When computing a value of your function, say, for t1=1e-15, the computer needs to subtract several intermediate results that are close to one another, and that will introduce artefacts.
As Walter has shown, you can alleviate that by doing these numerical calculations in the Symbolic Math Toolbox, where you can use higher-precision numbers than the doubles MATLAB is working in. But you should understand that just pushes the problem further out, and at some point you may need to increase the precision used, by calling digits(30), digits(50), or even higher values. Also, using variable precision arithmetic is much slower than using the double-precision numbers your computer has special hardware for.
syms t1
f = 17796509464316297.363895748*exp(-11.0*t1) - 17796509464316296.574074074074074*exp(-10.0*t1) - 0.68982167392592592592592592592593*exp(-1.0*t1) + 1.2222222222222222222222222222222*t1*exp(-1.0*t1) + 6174450735868760.1111111111111111*t1*exp(-10.0*t1) + 11622058728447497.43895748*t1*exp(-11.0*t1) - 982200909334137.5*t1^2*exp(-10.0*t1) + 93586137679683.333333333333333333*t1^3*exp(-10.0*t1) - 5859613763979.1666666666666666667*t1^4*exp(-10.0*t1) + 252901537750.0*t1^5*exp(-10.0*t1) - 8145447708.3333333333333333333333*t1^6*exp(-10.0*t1) + 219516865.07936507936507936507937*t1^7*exp(-10.0*t1) - 2683903.7698412698412698412698413*t1^8*exp(-10.0*t1) - 222938.71252204585537918871252205*t1^9*exp(-10.0*t1) + 12125.220458553791887125220458554*t1^10*exp(-10.0*t1) + 3706004905623298.2447874*t1^2*exp(-11.0*t1) + 767474314438508.23262466666666667*t1^3*exp(-11.0*t1) + 115820035783760.539895*t1^4*exp(-11.0*t1) + 13549259810802.40979*t1^5*exp(-11.0*t1) + 1276017179970.1820111111111111111*t1^6*exp(-11.0*t1) + 99090055157.583646825396825396825*t1^7*exp(-11.0*t1) + 6431882521.8889831349206349206349*t1^8*exp(-11.0*t1) + 349879383.20745425485008818342152*t1^9*exp(-11.0*t1) + 15772035.50003169091710758377425*t1^10*exp(-11.0*t1) + 573815.32070005611672278338945006*t1^11*exp(-11.0*t1) + 16540.373890295982309871198760088*t1^12*exp(-11.0*t1) + 427.06496366609213831436053658276*t1^13*exp(-11.0*t1) + 14.861012365197633054775911918769*t1^14*exp(-11.0*t1) + 0.44769060580072484834389596294358*t1^15*exp(-11.0*t1) + 0.0045334735715336740469015601290734*t1^16*exp(-11.0*t1) + 0.0044089949942699028444437719258759*t1^17*exp(-11.0*t1) + 0.00072758477778833722136969361785214*t1^18*exp(-11.0*t1) + 0.000077871208689571203517652792619439*t1^19*exp(-11.0*t1) + 0.000005672202354891583523257115371024*t1^20*exp(-11.0*t1) + 0.0000002583191156593700955608594990967*t1^21*exp(-11.0*t1) + 0.0000000065841638862926674918780324157471*t1^22*exp(-11.0*t1) + 0.00000000007601031748692706747794486283681*t1^23*exp(-11.0*t1) + 0.00000000000029107980533995897383819762284709*t1^24*exp(-11.0*t1) - 0.00000000000000031590056393483919641354780774103*t1^25*exp(-11.0*t1) - 2.0;
f1 = @(x) subs(f,t1,x); % force MATLAB to use vpa for fplot
fplot(f1,[0 1])
You may or may not be able to use the same workaround when calling optimization functions - but it will cost you a lot of time, and it is probably better to find a numerically more stable form of your input. You may even want to think about approximating it. Here is one way you can do that, using a plot to assess graphically how good the approximation is. I am assuming you already decided to focus on the range from 0 to 1, but the concept should easily generalize:
f2 = series(f,t1,1/2,'Order',5);
f3 = series(f,t1,1/2,'Order',10);
fplot(f1,[0,1],'-.','LineWidth',1.4)
hold on
fplot(f2,[0,1])
fplot(f3,[0,1],'k')
legend(["Original","Order 5 approx","Order 10 approx"])
If, looking at this, you decide the order 10 approximation is a close enough fit in the range you are interested in, using that for the optimization should give you results very quickly.

类别

在 Help Center 和 File Exchange 中查找有关 Solver-Based Optimization Problem Setup 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by