現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
How to approximate a multivariable arctan function?
10 ビュー (過去 30 日間)
古いコメントを表示
Hi,
let the following function be given:
Is there a way in Matlab to approximate this function as a multivariable (x1,x2,x3,u1) rational polynomial?
lv and bv are positive constants.
x2 should be defined within -pi/4 and pi/4 [rad].
x3 should be defined within (-pi/4) and (pi/4) [rad/sec].
x1 > 0 [meter/sec].
u1 should be defined within (-pi/12) and (pi/12) [rad].
Thanks a lot in advance!
19 件のコメント
Jannis Erz
2022 年 9 月 18 日
I'd like to represent the function as a sum-of-squares to perform lyapunov stability analysis.
I need the function to be represented at least as a rational polynomial including four variables to be able to do that.
Jannis Erz
2022 年 9 月 18 日
This was just an example of a rational model! I'm looking for a rational model incl. 4 independent variables of course!
Torsten
2022 年 9 月 18 日
編集済み: Torsten
2022 年 9 月 18 日
Ok.
Then evaluate your function for a number of points over the relevant intervals and make a least-square fit for the coefficient of the rational function you want to use.
A MATLAB program you could use for this task is lsqcurvefit.
But atan has a very restricted codomain. I wonder how rational functions in 4 variables should be able to approximate it.
Sam Chak
2022 年 9 月 18 日
編集済み: Sam Chak
2022 年 9 月 18 日
Hi @Jannis Erz
Do you think it is feasible to work out something like this
for the approximation of
where
and then using the Symbolic Math Toolbox to transform it into a your desired Rational Polynomial function?
Jannis Erz
2022 年 9 月 19 日
Dear @Sam Chak
thanks a lot for the effort calculating the series expansion! However, I'm afraid that your proposed solution might be unfeasible for my application.
Let's discuss that down below!
Jannis Erz
2022 年 9 月 19 日
I tried to implement a first version regarding lsqcurvefit to a slightly different atan related function. However, I always get an error message as soon as I add rational parts to the approximate function. Could you be so kind and have a look to my code (see below)?
The rational function to approximate F_q has to look somehow like this one:
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
F_z = 3000; % Actual vertical wheel load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = [-20:0.5:20];
alpha = deg2rad(alpha_deg);
for idx1 = 1:length(alpha)
Dy0 = muy0 * (F_z/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
df_z = ((F_z/1000)-F_z0)/F_z0; % Normalized change in vertical load % Normalized change in vertical load
F_q(idx1) = (Dy0 * sin(Cy*atan(By0*(alpha(idx1)) - Ey*(By0*(alpha(idx1))-atan(By0*(alpha(idx1)))))))*1000; % Pure Lateral Force [N]
end
fun = @(p,alpha)((p(1).*alpha.^3+p(2).*alpha.^2+p(3).*alpha)./((alpha.^4)+p(4).*alpha.^2)+p(5));
% Rational Function to approximate Lateral Tyre Force
p0 = [0.1,0.1,0.1,0.1,0.1];
p = lsqcurvefit(fun,p0,alpha,F_q);
plot(rad2deg(alpha),F_q,'o')
hold on
plot(rad2deg(alpha),fun(p,alpha))
xlabel('Tyre Side Slip Angle [deg]')
ylabel('Pure Lateral Tyre Force [N]')
legend('Pacejka Data','Fitted exponential')
Torsten
2022 年 9 月 19 日
fun = @(p,alpha)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5));
instead of
fun = @(p,alpha)((p(1).*alpha.^3+p(2).*alpha.^2+p(3).*alpha)./((alpha.^4)+p(4).*alpha.^2)+p(5));
Torsten
2022 年 9 月 19 日
編集済み: Torsten
2022 年 9 月 19 日
The poles will make trouble.
Maybe you can formulate a constraint on p(4) and p(5) such that the roots are outside the interval of approximation. Then use fmincon instead of lsqcurvefit and define this constraint in the nonlcon function.
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
F_z = 3000; % Actual vertical wheel load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = [-20:0.5:20];
alpha = deg2rad(alpha_deg);
for idx1 = 1:length(alpha)
Dy0 = muy0 * (F_z/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
df_z = ((F_z/1000)-F_z0)/F_z0; % Normalized change in vertical load % Normalized change in vertical load
F_q(idx1) = (Dy0 * sin(Cy*atan(By0*(alpha(idx1)) - Ey*(By0*(alpha(idx1))-atan(By0*(alpha(idx1)))))))*1000; % Pure Lateral Force [N]
end
fun = @(p)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)-F_q.*(alpha.^4+p(4)*alpha.^2+p(5));
fun1 = @(p)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5));
%fun = @(p,alpha)((p(1).*alpha.^3+p(2).*alpha.^2+p(3).*alpha)./((alpha.^4)+p(4).*alpha.^2)+p(5));
% Rational Function to approximate Lateral Tyre Force
p0 = [0.1,0.1,0.1,0.1,0.1];
format long
p = lsqnonlin(fun,p0)
Local minimum possible.
lsqnonlin stopped because the size of the current step is less than
the value of the step size tolerance.
p = 1×5
1.0e+03 *
1.573707413267325 0.000000000311668 -0.034491567121686 0.000043504411553 -0.000001348685846
rad2deg(roots([1 0 p(4) 0 p(5)]))
ans =
-0.000000000000003 +14.544020290876990i
-0.000000000000003 -14.544020290876990i
-8.289268225905060 + 0.000000000000000i
8.289268225905060 + 0.000000000000000i
plot(rad2deg(alpha),F_q,'o')
hold on
plot(rad2deg(alpha),fun1(p))
xlabel('Tyre Side Slip Angle [deg]')
ylabel('Pure Lateral Tyre Force [N]')
legend('Pacejka Data','Fitted exponential')
Torsten
2022 年 9 月 19 日
編集済み: Torsten
2022 年 9 月 19 日
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
F_z = 3000; % Actual vertical wheel load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = [-20:0.1:20];
alpha = deg2rad(alpha_deg);
for idx1 = 1:length(alpha)
Dy0 = muy0 * (F_z/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
df_z = ((F_z/1000)-F_z0)/F_z0; % Normalized change in vertical load % Normalized change in vertical load
F_q(idx1) = (Dy0 * sin(Cy*atan(By0*(alpha(idx1)) - Ey*(By0*(alpha(idx1))-atan(By0*(alpha(idx1)))))))*1000; % Pure Lateral Force [N]
end
fun = @(p)sum(((p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5))-F_q).^2);
fun1 = @(p)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5));
p0 = [0.1,0.1,0.1,0.1,0.1];
options = optimset('MaxFunEvals',100000,'MaxIter',100000);
format long
p = fmincon(fun,p0,[],[],[],[],[],[],@nonlcon,options)
Local 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.
p = 1×5
1.0e+03 *
1.320317528765901 0.000000004546610 0.059807026912929 0.000090641897411 0.000002053988392
rad2deg(roots([1 0 p(4) 0 p(5)]))
ans =
-0.000001716290581 +12.197536562689598i
-0.000001716290581 -12.197536562689598i
0.000001716290582 +12.197536562700607i
0.000001716290582 -12.197536562700607i
plot(rad2deg(alpha),F_q,'o')
hold on
plot(rad2deg(alpha),fun1(p))
xlabel('Tyre Side Slip Angle [deg]')
ylabel('Pure Lateral Tyre Force [N]')
legend('Pacejka Data','Fitted exponential')
function [c,ceq] = nonlcon(p)
c = p(4)^2/4-p(5);
ceq = [];
end
Jannis Erz
2022 年 9 月 20 日
Thank you so much for your effort handing me over such a good solution!
There is one thing still unclear to me regarding the nonlinear equality constraint of the denominator (x^4 + p4*x^2 + p5). I use x instead of alpha here for the sake of convenience. Why do we want the inner sqrt to become negative?
I do understand that the roots of the first solution are causing the serrated process output when the denominator gets close to zero, i.e., close to the root values (+- 8.3°) on the real axis.
Torsten
2022 年 9 月 20 日
One way to avoid real poles in the interval of interest is to make all roots of the denominator complex. That's what I did. Of course, this is quite arbitrary and might exclude a better solution. Feel free to experiment.
Jannis Erz
2022 年 9 月 23 日
@Torsten Again thanks a lot for the hints! I played a little with the roots and the best solutions seems to make all roots of the denominator complex. I also tried out the Curve Fitter App using the Rational Option (only available at R2022b). With one varying variable it delivers perfect fitting using rationals with the degree 3 and 4 respectively. Unfortunately it is limited to one variable only!
Now I'm trying to extend the problem to two varying variables (see code below). However, I'm still not satisfied with the results (absolute error around 500N). I'd like to force the absolute error below 200 N. @John D'Errico and @Torsten do you have any ideas how to systematically vary with the rational expressions and/or starting points to get better results?
%% Version 2 of Rational Approximation of tyre lateral force "F_q" with two varying variables side slip angle "a" and wheel load "F_z"
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = -20:0.1:20;
a_vec = deg2rad(alpha_deg);
F_z_vec = linspace(500,7500,length(a_vec));
[a, F_z] = meshgrid(a_vec,F_z_vec);
for idx1 = 1:length(a)
for idx2 = 1:length(F_z)
Dy0 = muy0 * (F_z(idx1,idx2)/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z(idx1,idx2)/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
F_q(idx1,idx2) = (Dy0 * sin(Cy*atan(By0*(a(idx1,idx2)) - Ey*(By0*(a(idx1,idx2))-atan(By0*(a(idx1,idx2)))))))*1000; % Pure Lateral Force [N]
end
end
fun = @(p)sum((calcrational(p,a,F_z)-F_q).^2,'all'); % Problem to be minimized
fun1 = @(p)(calcrational(p,a,F_z)); % Rational polynomial approximation of F_q
p0(1,1:10) = 0.1;
options = optimset('MaxFunEvals',1000000,'MaxIter',1000000);
format long
p = fmincon(fun,p0,[],[],[],[],[],[],[],options)
F_q_approx = calcrational(p,a,F_z);
abs_err = abs(F_q_approx)-abs(F_q); % Absolute error
%mesh(a,F_z,F_q)
%hold on
%mesh(a,F_z,(calcrational(p,a,F_z)))
mesh(a,F_z,abs_err) % Display absolute error
xlabel('Side Slip Angle [rad]')
ylabel('Wheel Load [N]')
zlabel('Absolute Error [N]')
% Rational approximation function
function rational = calcrational(p,x,y)
num = p(1) + p(2)*x + p(3)*y + p(4)*x.^2 + p(5)*y.^2 + p(6)*x.*y;
denom = p(7) + p(8)*x + p(9)*y + p(10)*x.^2;
rational = num./denom;
end
Torsten
2022 年 9 月 23 日
Your problem is overfitted.
As you can see in MATLAB's rational models section, one of the parameters to be fitted must be normalized to 1.
Jannis Erz
2022 年 9 月 23 日
@Torsten what do you mean by normalized? Sorry, but I'm a little bit clueless here :)
How would I implement that in the code? Enforce a nonlinear equality constraint on, e.g., p(10) to be equal to 1?
Torsten
2022 年 9 月 23 日
編集済み: Torsten
2022 年 9 月 23 日
No, in the same way as the rational polynomials in MATLAB are defined - by setting p(10) = 1 and fitting only 9 parameters.
Note that if p(1),...,p(10) is a solution for the Least-squares problem, then also p(1)*c,p(2)*c,...,p(10)*c is a solution for every constant c not equal to 0. This is prohibited by normalizing one of the parameters to 1.
回答 (1 件)
John D'Errico
2022 年 9 月 18 日
編集済み: John D'Errico
2022 年 9 月 18 日
@Sam Chak has suggested an interesting idea in a comment. However, while it is not a bad idea at first glance, I suspect it will be problematic. The issue is, that classic atan series is not very strongly convergent. Expect to need many terms before you get anything good. And that means the polynomial approximation will be poor at best.
For example, how many terms do you need for convergence as a function of x to k significant digits, in the simple atan series?
For example, when x == 1, how many terms are required? We can look at it easily, since that is just the Leibniz formula for pi/4. Thus...
N = 0:1000;
S = mod(N+1,2)*2 - 1;
approx = cumsum(S./(2*N+1));
truth = pi/4;
semilogy(N,abs(approx - truth))
grid on
So we need thousands of terms for convergence near x==1.
However, if you can use range reduction methods to force the argument x to the atan to be in a small interval near zero, you get much faster convergence. The problem is, that in itself may take some work, and it is not at all trivial to get good convergence.
Instead, you may gain from using a direct rational polynomial approximation, perhaps something like those found in the classic, by JF Hart, et al, "Computer Approximations". (Start reading around page 120 in my edition from 1978. The tables there give some pretty good approxmations, though they still employ range reduction.)
Note: Even though Hart is an old text, it is still a book I love dearly, but that might apply only to a real gearhead numerical analyst like me. I think I recall it is available as a Dover reprint.
7 件のコメント
Sam Chak
2022 年 9 月 19 日
@John D'Errico, thanks for the heads-up. @Jannis Erz has mentioned about using the approximated Rational Polynomial function for some kind of Lyapunov Stability analysis.
Since the range for and are given, I guess that the stability analysis may be applied on the vicinity of some equilibrium points like and , and therefore
.
Does @Jannis Erz need to approximate until for the purpose of showing the stability proof around the equilibrium points?
Jannis Erz
2022 年 9 月 19 日
Dear @John D'Errico thanks a lot for the recommendations! I already got myself the 1978 version of your mentioned book! Let's see how I can get along with it :)
@Sam Chak You're right that I'm examining the equilibrium points for stability proof. In addition to that I'd like to analyse the domain of attraction. Therefore I need the approximation to be within -pi/4 and + pi/4.
John D'Errico
2022 年 9 月 19 日
編集済み: John D'Errico
2022 年 9 月 19 日
@Sam Chak - the point is, there are arguably multiple ways to generate a rational polynomial approximation to a function. However, if you START with a poorly convergent Taylor series, you will likely not be happy with the result. This is my fear, and why I would suggest using a rational polynomial approximation itself for the atan which can have better properties over the desired domain.
@Jannis Erz - You can learn a lot from the Hart text. It remains one of the books I leave on my shelf next to my desk. (Abramowitz & Stegun is another, where my copy has tabs attached to the sections I have used regularly over the years.)
Bruno Luong
2022 年 9 月 19 日
編集済み: Bruno Luong
2022 年 9 月 19 日
@John D'Errico: not necessary. It is well known that Padé fractional series converges better much than Taylor series, but the Padé series is derived directly from the Taylor series.
That being say, the best approximation on a given interval is just from numerical calculation.
Sam Chak
2022 年 9 月 19 日
OP has changed the function to:
F_z0 = 3;
F_z = 3000;
muy0 = 1;
Dy0 = muy0*(F_z/1000);
c1 = 8;
c2 = 1.33;
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0)));
Cy = 1.3;
By0 = CFa0/(Cy*Dy0);
Ey = -1;
alpha_deg = linspace(-20, 20, 401);
alpha = deg2rad(alpha_deg);
Fq = 1000*Dy0*sin(Cy*atan(By0*alpha - Ey*(By0*alpha - atan(By0*alpha))));
plot(alpha_deg, Fq, 'linewidth', 1.5), ylim([-4e3 4e3])
grid on, xlabel('\alpha, deg'), ylabel('F_{q}')
Bruno Luong
2022 年 9 月 19 日
Yes, Padé fractional series is my loosly wording, rather use Padé approximant please.
Sam Chak
2022 年 9 月 19 日
help \sym\pade
PADE approximation of a symbolic expression
PADE(f <, x> <,options>) computes the third order Pade approximant of f at x = 0.
If x is not given, it is determined using symvar.
PADE(f, x, a <, options>) computes the third order Pade approximant of f at x = a.
A different order can be specified using the option 'Order' (see below).
The following options can be given as name-value pairs:
Parameter Value
'ExpansionPoint' a
Compute the Pade approximation about the point a. It is also
possible to specify the expansion point as third argument
without explicitly using a parameter value pair.
'Order' [m, n]
Compute the Pade approximation with numerator order m and
denominator order n.
'Order' m
Compute the Pade approximation with both numerator and
denominator order equal to m. By default, 3 is used.
'OrderMode' 'Absolute' or 'Relative'.
If the order mode is 'Relative' and f has a zero or pole at the
expansion point, add the multiplicity of that zero to the order.
If f has no zero or pole at the expansion point, this option has
no effect. Default is 'Absolute'.
Examples:
syms x
pade(sin(x))
returns (60*x - 7*x^3)/(3*(x^2 + 20))
pade(cos(x), x, 'Order', [1, 2])
returns 2/(x^2 + 2)
See also TAYLOR.
Documentation for pade
doc pade
参考
カテゴリ
Help Center および File Exchange で Polynomials についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)