Standard form of Transfer function.

Hi there,

is there a way in MATLAB (Symbolic math toolbox) to go from this form:

to this form:

knowing that 'A' is big coefficient itself:

Answer 1

[Sam Chak]

Hi @Hamza Harbi

This display of the transfer function is not referred to as the 'Standard form', or rather, it is known as the 'Natural Frequency form.' However, I am not aware of a single function in the Symbolic Math Toolbox that can achieve this. In fact, it is significantly easier to obtain this display using the zpk() function from the Control System Toolbox.

s   = tf('s');

%% Transfer function in raw form
G   = (5*s^2 + 25*s + 30)/(6*s^3 + 42*s^2 + 72*s + 60)

G =
 
      5 s^2 + 25 s + 30
  --------------------------
  6 s^3 + 42 s^2 + 72 s + 60
 
Continuous-time transfer function.

%% Convert transfer function to Frequency form
G   = zpk(G);
G.DisplayFormat = 'frequency'

G =
 
             0.5 (1+s/3) (1+s/2)
  ------------------------------------------
  (1+s/5) (1 + 1.414(s/1.414) + (s/1.414)^2)
 
Continuous-time zero/pole/gain model.

[Sam Chak]

Hi @Hamza Harbi

Due to the unique structure of the 4th-order transfer function , we can perform pole-zero cancellation to obtain a minimized 2nd-order transfer function. Since both the numerator and denominator are 2nd-degree polynomials, we can easily perform algebraic calculations for the coefficients of the zero-pole-gain model and express the transfer function in natural frequency form.

%% Symbolic approach
syms s  A  R0 V0 C1 D0 D1 L1 L2

%% terms
n1  = V0/D0^2;
n2  = D1/(D0*L1);
n3  = C1/D0;
n4  = V0/(D0*D1*R0);
d1  = (D1^2)/(L1*D0^2);
d2  = 1/L2;
d3  = C1/(D0^2);
d4  = 1/L2;

%% numerator & denominator
num = n1*(n2/s + n3*s) - n4;
den = A*(d1/s + d2/s + d3*s) - d4;

%% original transfer function for the RLC circuit
Gvd = num/den

Gvd = 

Gvd = simplify(Gvd, 'Steps', 100)

Gvd = 

[N, D]  = numden(Gvd);
expN    = expand(N);
expD    = expand(D);
[Ncoeffs, Nterms] = coeffs(expN, s)

Ncoeffs = 

Nterms = 

[Dcoeffs, Dterms] = coeffs(expD, s)

Dcoeffs = 

Dterms = 

%% parameters (supplied by the user)
A   = sym(2);
R0  = sym(3);    % Ro
V0  = sym(5);    % Vo
C1  = sym(7);
D0  = sym(11);   % D
D1  = sym(13);   % D' (Dprime)
L1  = sym(17);
L2  = sym(19);

sympref("FloatingPointOutput", true);

%% coefficients of system zeros (by formulas in natural frequency form)
an  =   C1*D1*L1*L2*R0*V0;                  % a*s^2
bn  = - D0^2*L1*L2*V0;                      % b*s
cn  =   D1^2*L2*R0*V0;                      % c
wn  = sqrt(cn/an);                          % omega (natural frequency)
disp(vpa(wn))

0.3305

zn  = wn*bn/cn;                             % coefficient of s/w
disp(vpa(zn))

%% coefficients of system poles (by formulas in natural frequency form)
ad  =   A*C1*D0*D1*L1*L2*R0;                % a*s^2
bd  = - D0^3*D1*L1*R0;                      % b*s
cd  =   A*L1*R0*D0^3*D1 + A*L2*R0*D0*D1^3;  % c
wd  = sqrt(cd/ad);                          % omega (natural frequency)
disp(vpa(wd))

1.5264

zd  = wd*bd/cd;                             % coefficient of s/w
disp(vpa(zd))

%% system gain
k   = cn/cd;
disp(vpa(k))

0.0107

%% Manually express the transfer function in natural frequency form
Gnf = @(s) k*(1 + zn*(s/wn) + (s/wn)^2)/(1 + zd*(s/wd) + (s/wd)^2)

Gnf = function_handle with value:
    @(s)k*(1+zn*(s/wn)+(s/wn)^2)/(1+zd*(s/wd)+(s/wd)^2)

f   = symfun(Gnf(s), s)

f(s) = 

formula(Gnf(s))

ans = 

At this point, I am unsure which functions in the MATLAB or Symbolic Math Toolbox can be utilized to prevent expansion while preserving the exact raw form of the natural frequency form that allows further symbolic computation.

However, we can verify the computated coefficients of the symbolic zero-pole-gain model using the zpk() function from the Control System Toolbox to obtain this display of the transfer function in natural frequency form.

%% Declare Laplace variable, s
s   = tf('s');

%% parameters
A   = 2;
R0  = 3;    % Ro
V0  = 5;    % Vo
C1  = 7;
D0  = 11;   % D
D1  = 13;   % D'
L1  = 17;
L2  = 19;

%% terms
n1  = V0/D0^2;
n2  = D1/(D0*L1);
n3  = C1/D0;
n4  = V0/(D0*D1*R0);
d1  = (D1^2)/(L1*D0^2);
d2  = 1/L2;
d3  = C1/(D0^2);
d4  = 1/L2;

%% numerator & denominator
num = n1*(n2/s + n3*s) - n4;
den = A*(d1/s + d2/s + d3*s) - d4;

%% original transfer function
Gvd = num/den

Gvd =
 
  0.0263 s^4 - 0.01166 s^3 + 0.002873 s^2
  ---------------------------------------
   0.1157 s^4 - 0.05263 s^3 + 0.2696 s^2
 
Continuous-time transfer function.

%% after pole-zero cancellation at the origin
Gvd = minreal(Gvd)

Gvd =
 
  0.2273 s^2 - 0.1007 s + 0.02483
  -------------------------------
       s^2 - 0.4549 s + 2.33
 
Continuous-time transfer function.

%% Display transfer function in Frequency form
Gvd = zpk(Gvd);
Gvd.DisplayFormat = 'frequency'

Gvd =
 
  0.010656 (1 - 1.341(s/0.3305) + (s/0.3305)^2)
  ---------------------------------------------
       (1 - 0.298(s/1.526) + (s/1.526)^2)
 
Continuous-time zero/pole/gain model.

[Hamza Harbi]

Thank you @Sam Chak for your effort, that actually helped me alot.

Answer 2

[Walter Roberson]

An intermediate step in converting to natural frequency form is finding the roots of the denominator.

It happens that it is possible to solve the denoninator for G_vd symbolically. However, the solution involves finding the four roots of

(C_1*C_o*L_1*L_2*R_c + C_1*C_o*L_1*L_2*R_o)*z^4 + (C_1*L_1*L_2 + C_1*C_o*L_1*R_c*R_o)*z^3 + (C_1*L_1*R_o + C_o*D^2*L_1*R_c + C_o*D^2*L_1*R_o + C_o*D_prime^2*L_2*R_c + C_o*D_prime^2*L_2*R_o)*z^2 + (D^2*L_1 + D_prime^2*L_2 + C_o*D_prime^2*R_c*R_o)*z + D_prime^2*R_o

in z. That is possible explicitly, but each of the roots is over 37000 characters.

So although it is theoretically possible to create the symbolic natural frequency form of that equation, in practice doing so is crap, a meaningless rather long expression. Very very few humans can make sense of the resulting expression. If I were to change the sign at a random location in the expression, it is very doubtful that you would be able to pick out the location by reading the expression.