How to approximate a multivariable arctan function?

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!

[Torsten]

Is there a way in Matlab to approximate this function as a multivariable (x1,x2,x3,u1) rational polynomial?

What's the reason you want to do that ?

What's a rational polynomial ? A polynomial with rational coefficients ?

[Jannis Erz]

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.

[Torsten]

Strange that the exponent of x in the denominator of y is m-1 and not m-i ...

I don't see how you want to cope with the multidimensionality of your function. The rational model only has 1 independent variable while your function has 4.

[Jannis Erz]

This was just an example of a rational model! I'm looking for a rational model incl. 4 independent variables of course!

[Torsten]

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]

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]

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]

Dear @Torsten,

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]

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]

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]

% 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

[Sam Chak]

@Torsten 👍. Should move that solution for the Rational Polynomial function to the Answer section.

[Jannis Erz]

@Torsten

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]

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]

@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]

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]

@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]

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.

[Jannis Erz]

Thanks for the explanation! Really appreciate it!

Answer 1

[John D'Errico]

@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.

[Sam Chak]

@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]

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]

@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]

@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.

So I wouldn't reject @Sam Chak proposition even if he sarts with poor Taylor series.

That being say, the best approximation on a given interval is just from numerical calculation.

[Sam Chak]

@Bruno Luong, thanks. Learned new things. Is Padé fractional series also known as Padé approximant?

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]

Yes, Padé fractional series is my loosly wording, rather use Padé approximant please.

[Sam Chak]

Found the pade() command. Maybe @Jannis Erz can work out the desired Rational Polynomial function.

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