# Problem when utilizing the bode function to plot the multiplication of several transfer functions

3 views (last 30 days)
Tsai Han Hao on 18 Nov 2022
Edited: Tsai Han Hao on 18 Nov 2022
Hi all,
I'm right now trying to plot the bode plot of a system that contains multiple transfer fucntions, for example: G1*G2*G3.
(all the tranfer functions are declared by using tf() )
However, it just occurred to me that the bode plot of G1*G2*G3 is different from G2*G1*G3. And I felt quite confused about the results since mathematically they should be the same?
==> This is the results of the code below, as you can see, the bode plots are different.
The example code is just similar to:
%% Initialization
omega_D_Val = 50;
M_Val = 4;
Kf_Val = 40;
Mn_Val = 4;
Kfn_Val = 40;
Ts_Val = 0.0005;
wR_Larger = 3000;
Cf_Val = 1;
Pure_gain = (M_Val/70) * wR_Larger * (1/M_Val);
D_env_Val = 1.5;
K_env_Val = 2500;
Integrator = tf(Ts_Val,[1 -1],Ts_Val);
%% Plotting Part
figure;
set(gcf,'Units','centimeters','position',[0 0 25 30]);
[~, Tenv_None, Nenv_None, ~, OSenv_None]...
= OuterLoop_TransferFcn_withEnv(omega_D_Val, wR_Larger, Ts_Val, (M_Val/70), Mn_Val, M_Val, Kfn_head_Val, Kfn_Val, Kf_Val, D_env_Val, K_env_Val, Cf_Val);
[Gain, ~, Frequency] = bode(Pure_gain * Integrator * Tenv_None * OSenv_None, {1e-4,6283.19});
semilogx(reshape(Frequency,[],1), 20*log10(reshape(Gain,[],1)), "-", 'Color', 'k', 'LineWidth',1.5);
hold on
[Gain, ~, Frequency] = bode(Integrator * Tenv_None * OSenv_None * Pure_gain, {1e-4,6283.19});
semilogx(reshape(Frequency,[],1), 20*log10(reshape(Gain,[],1)), "--", 'Color', 'k', 'LineWidth',1.5);
%% Function for tranfer functions build up
function [G_env, T_env, N_env, Compensator_env, OuterLoop_Sensitivity_Function_env]...
= OuterLoop_TransferFcn_withEnv(omega_D_Val, omega_R_Val, Ts_Val, Mn_head_Val, Mn_Val, M_Val, Kfn_head_Val, Kfn_Val, Kf_Val, D_env_Val, K_env_Val, Cf_Val)
syms omega_R omega_D Ts Mn_head Mn M Kfn_head Kfn Kf D_env K_env
syms z
Cf = Cf_Val;
% Alpha and Beta Define
alpha = (Mn*Kf)/(M*Kfn);
Assigned_Parameters = [omega_D, omega_R, Ts, Mn_head, Mn, M, Kfn_head, Kfn, Kf, D_env, K_env];
Assigned_Values = [omega_D_Val, omega_R_Val, Ts_Val, Mn_head_Val, Mn_Val, M_Val, Kfn_head_Val, Kfn_Val, Kf_Val, D_env_Val, K_env_Val];
Denominator = (2*(z - 1)*((2 + omega_R*Ts)*z + (-2 + omega_R*Ts))*(M*z^2 + (M*alpha*omega_D*Ts + D_env*Ts - 2*M)*z + (M + K_env*Ts^2 - D_env*Ts - M*alpha*omega_D*Ts)))...
+ (Mn_head*Cf*omega_R*Ts*((2 + omega_D*Ts)*z + (-2 + omega_D*Ts))*(M*beta*z^3 + (D_env*Ts*beta - 2*M*alpha - M*beta)*z^2 + (K_env*beta*Ts^2 + 4*M*alpha -M*beta)*z + (K_env*beta*Ts^2 - D_env*beta*Ts - 2*M*alpha + M*beta)));
Denominator = subs(Denominator, Assigned_Parameters, Assigned_Values);
G_tf_Numerator = (Mn_head*omega_R*Ts*((2 + omega_D*Ts)*z + (-2 + omega_D*Ts))*(M*beta*z^3 + (D_env*Ts*beta - 2*M*alpha - M*beta)*z^2 + (K_env*beta*Ts^2 + 4*M*alpha -M*beta)*z + (K_env*beta*Ts^2 - D_env*beta*Ts - 2*M*alpha + M*beta)));
G_tf_Denominator = (2*(z - 1)*((2 + omega_R*Ts)*z + (-2 + omega_R*Ts))*(M*z^2 + (M*alpha*omega_D*Ts + D_env*Ts - 2*M)*z + (M + K_env*Ts^2 - D_env*Ts - M*alpha*omega_D*Ts)));
G_tf_Numerator = subs(G_tf_Numerator, Assigned_Parameters, Assigned_Values);
G_tf_Denominator = subs(G_tf_Denominator, Assigned_Parameters, Assigned_Values);
[G_tf_Numerator, G_tf_Denominator] = numden(G_tf_Numerator/G_tf_Denominator);
[bG,~] = coeffs(G_tf_Numerator, z, 'All');
[bG_D,~] = coeffs(G_tf_Denominator, z, 'All');
T_env_Numerator = M*z^2 + (alpha*omega_D*Ts*M - 2*M)*z + (M - alpha*omega_D*Ts*M);
T_env_Denominator = M*z^2 + (D_env*Ts - 2*M + alpha*omega_D*Ts*M)*z + (K_env*Ts^2 - D_env*Ts + M - alpha*omega_D*Ts*M);
T_env_Numerator = subs(T_env_Numerator, Assigned_Parameters, Assigned_Values);
T_env_Denominator = subs(T_env_Denominator, Assigned_Parameters, Assigned_Values);
[T_env_Numerator, T_env_Denominator] = numden(T_env_Numerator/T_env_Denominator);
[bT,~] = coeffs(T_env_Numerator, z, 'All');
[bT_D,~] = coeffs(T_env_Denominator, z, 'All');
N_env_Numerator = M*z^2 + (D_env*Ts - 2*M)*z + (M + K_env*Ts^2 - D_env*Ts);
N_env_Denominator = M*z^2 + (D_env*Ts - 2*M + M*alpha*omega_D*Ts)*z + (M + K_env*Ts^2 - D_env*Ts - M*alpha*omega_D*Ts);
N_env_Numerator = subs(N_env_Numerator, Assigned_Parameters, Assigned_Values);
N_env_Denominator = subs(N_env_Denominator, Assigned_Parameters, Assigned_Values);
[N_env_Numerator, N_env_Denominator] = numden(N_env_Numerator/N_env_Denominator);
[bN,~] = coeffs(N_env_Numerator, z, 'All');
[bN_D,~] = coeffs(N_env_Denominator, z, 'All');
C_tf_Numerator = (Mn_head*Cf*omega_R*Ts*((2 + omega_D*Ts)*z + (-2 + omega_D*Ts))*(M*beta*z^3 + (D_env*Ts*beta - 2*M*alpha - M*beta)*z^2 + (K_env*beta*Ts^2 + 4*M*alpha -M*beta)*z + (K_env*beta*Ts^2 - D_env*beta*Ts - 2*M*alpha + M*beta)));
C_tf_Numerator = subs(C_tf_Numerator, Assigned_Parameters, Assigned_Values);
[C_tf_Numerator, C_tf_Denominator] = numden(C_tf_Numerator/Denominator);
[bC,~] = coeffs(C_tf_Numerator, z, 'All');
[bC_D,~] = coeffs(C_tf_Denominator, z, 'All');
OS_tf_Numerator = (2*(z - 1)*((2 + omega_R*Ts)*z + (-2 + omega_R*Ts))*(M*z^2 + (M*alpha*omega_D*Ts + D_env*Ts - 2*M)*z + (M + K_env*Ts^2 - D_env*Ts - M*alpha*omega_D*Ts)));
OS_tf_Numerator = subs(OS_tf_Numerator, Assigned_Parameters, Assigned_Values);
[OS_tf_Numerator, OS_tf_Denominator] = numden(OS_tf_Numerator/Denominator);
[bOS,~] = coeffs(OS_tf_Numerator, z, 'All');
[bOS_D,~] = coeffs(OS_tf_Denominator, z, 'All');
% Transfer Function Declarations
G_env = tf(double(bG)/double(bG_D(1)), double(bG_D)/double(bG_D(1)), Ts_Val);
T_env = tf(double(bT)/double(bT_D(1)), double(bT_D)/double(bT_D(1)), Ts_Val);
N_env = tf(double(bN)/double(bN_D(1)), double(bN_D)/double(bN_D(1)), Ts_Val);
Compensator_env = tf(double(bC)/double(bC_D(1)), double(bC_D)/double(bC_D(1)), Ts_Val); % (G*Cf)/(1 + G*Cf)
OuterLoop_Sensitivity_Function_env = tf(double(bOS)/double(bOS_D(1)), double(bOS_D)/double(bOS_D(1)), Ts_Val); % 1/(1 + G*Cf)
end

Tsai Han Hao on 18 Nov 2022
Edited: Tsai Han Hao on 18 Nov 2022
Hi all,
I've found out that the tf data type is the main problem that causes the difference between the Bode plot of bode(G1*G2*G3) and bode(G2*G3*G1).
If you look through the Algorithms on how Matlab do the bode calculations, you can find out that zpk is the data type that Matlad do the Bode plot.
https://www.mathworks.com/help/ident/ref/lti.bode.html
Therefore, by utilizing bode(zpk(G1)*zpk(G2)*zpk(G3)), the results of the Bode plot will be correct no matter how you arrange the tranfer functions.