Mismatch between the frequency response of a transfer function and bode plot

I generate a state space model as follows (The details are not important. This part is used to get matrix A_xi and B_xi and the state space model is xi(k+1) = A_xi\*xi(k)+B1_xi\*u0(k); y(k) = x(k)):

% vehicle dynimcs

h = 0.8; % time headway

N = 4; % number of vehicles 0 leading

tau = 0.2; % actuator lag

Ts = 0.1;% sampling interval

% control gain

kp = 0.45;

kv = 0.84;

ka = 0.5;

% vehicle dynamics

% x_i = [e_i;dv_i;a_i], ei=d_i-h*v_i, dvi = v_i-1-v_i

A1 = [0,1,-h;0 0 -1; 0 0 -1/tau]; % x_i to x_i

A2 = [0 0 0; 0 0 1;0 0 0]; % x_i-1 to x_i

B = [0; 0; 1/tau]; % u_i to x_i

% x0 = [q0;v0;a0]

A0 = [0 1 0;0 0 1; 0 0 -1/tau]; % x0 to x0

B0 = B; % u0 to x0

Ac = zeros(3*N,3*N);

Ac(1:3,1:3) = A0;

Ac(5,3) = 1;

for n = 1:N-1

Ac(3*n+1:3*n+3,3*n+1:3*n+3) = A1;

if n < N-1

Ac(3*n+4:3*n+6,3*n+1:3*n+3) = A2;

end

end

Bc = zeros(3*N,N);

Bc(1:3,1) = B0;

for n = 1:N-1

Bc(3*n+1:3*n+3,n+1) = B;

end

veh_dyn_c = ss(Ac,Bc,eye(3*N),zeros(3*N,N));

veh_dyn = c2d(veh_dyn_c,Ts);

% kalman filter

% x = [d,dv,da]

Ak = [0 1 0; 0 0 1; 0 0 -1/tau];

Ck = [0; 0; -1/tau]; % ui to x known input

Bk = [0; 0; 1/tau]; % ui-1 to x disturbance

H = [1 0 0; 0 1 0];

sysk_c = ss(Ak, [Ck, Bk], H, zeros(2,2));

sysk = c2d(sysk_c,Ts);

delta_d = 0.5; % resolution

delta_v = 0.6;%0.6;%

SNR_db = 20;

SNR = 10^(SNR_db/10);

sigma_d = delta_d*sqrt(3/(2*pi^2*SNR));% standard deviation

sigma_v = delta_v*sqrt(3/(2*pi^2*SNR));% standard deviation (don't forget to square)

var_jerk = 1;%0.15499;% process noise variance (jerk)

[KEST,~,~,L0,~,~] = kalman(sysk,var_jerk, diag([sigma_d^2,sigma_v^2])); % L: innovation gain

L = L0;

% whole dynamics

% xi = [x0;x_1;d_1;u_1^-;\hat_{x_1^-};...;x_N-1;d_N-1;u_N-1^-;\hat_{x_N-1^-}]

A_all = zeros(3+8*(N-1),3+8*(N-1)); % xi to xi^+

A_f_all = zeros(3+8*(N-1),3+8*(N-1)); % xi^+ to xi^+

B_all = zeros(3+8*(N-1),1); % u0 to to xi^+

C_all = zeros(3+8*(N-1),2*(N-1)); % epsilon = [epsilon_d1;epsilon_v1;...;epsilon_dN-1;epsilon_vN-1] to xi^+

% x0

A_all(1:3,1:3) = veh_dyn.A(1:3,1:3); % x0 to x0+

B_all(1:3) = veh_dyn.B(1:3,1); % u0 to x0+

% x_i

A_all(4:6,1:3) = veh_dyn.A(4:6,1:3); % x_0 to x_1^+

B_all(4:6) = veh_dyn.B(4:6,1); % u0 to x_1^+

for id1 = 1:N-1

for id2 = 1:N-1

A_all(8*id1-4:8*id1-2,8*id2-4:8*id2-2) = veh_dyn.A(3*id1+1:3*id1+3,3*id2+1:3*id2+3); % x_id2 to x_id1^+

A_f_all(8*id1-4:8*id1-2,8*id2) = veh_dyn.B(3*id1+1:3*id1+3,id2+1); % u_id2 to x_id1^+

end

end

% u_i

for id = 1:N-1

A_all(8*id,8*id-4) = kp; % e_i to u_i

A_all(8*id,8*id-3) = kv; % dv_i to u_i

A_all(8*id,8*id-2) = ka; % a_i to u_i

A_f_all(8*id,8*id+3) = ka; % \hat_da_i to u_i

C_all(8*id,2*id-1) = kp; % epsilon_d_i to u_i

C_all(8*id,2*id) = kv; % epsilon_v_i to u_i

end

% d_i

for id = 1:N-1

A_f_all(8*id-1,8*id-4) = 1; % e_i to d_i

A_f_all(8*id-1,2) = h; % v0 to d_i

for id2 = 1:id

A_f_all(8*id-1,8*id2-3) = -h;

end

end

% \hat{x}_i

for id = 1:N-1

A_all(8*id+1:8*id+3, 8*id-1) = L(:,1); % d_i to \hat{x}_i

A_all(8*id+1:8*id+3, 8*id-3) = L(:,2); % dv_i to \hat{x}_i

C_all(8*id+1:8*id+3, 2*id-1) = L(:,1); % epsilon_d_i to \hat{x}_i

C_all(8*id+1:8*id+3, 2*id) = L(:,2); % epsilon_v_i to \hat{x}_i

A_all(8*id+1:8*id+3, 8*id+1:8*id+3) = (eye(3)-L*H)*sysk.A; % \hat{x}_i^- to \hat{x}_i

A_all(8*id+1:8*id+3, 8*id) = (eye(3)-L*H)*sysk.B(:,1); % u_i^- to \hat{x}_i

end

% feedback

A_xi = (eye(length(A_f_all))-A_f_all)\A_all;

B1_xi = (eye(length(A_f_all))-A_f_all)\B_all; % u0

B2_xi = (eye(length(A_f_all))-A_f_all)\C_all; % epsilon

Then I want to obtain the tranfer function between the 13th and 5th output as follows:

[NUM_u,DEN_u] = ss2tf(A_xi,B1_xi,eye(length(A_xi)),zeros(length(A_xi),size(B1_xi,2)),1);

Gamma_v21 = tf(NUM_u(13,:),NUM_u(5,:),Ts); % v2/v1

figure

bode(Gamma_v21)

The result is as follows:

However, when I use a sinunoid signal as u0 and plot these two outputs, the result does not match the Bode plot. For example, when the frequency of u0 is 0.1 rad/s, the magnitude of Bode diagram is -37.4dB. The simulated result is generated with the following code:

L_t = 5000;

t = (1:L_t)*Ts;

u0 = sin(0.1*t);

xi = zeros(3+8*(N-1),L_t);

xi(1:3,1)=[0;20;0];

for id = 1:N-1

xi(8*id-1,1) = h*xi(2,1);% d_i

%[id size(xi)]

xi(8*id+1,1) = xi(8*id-1,1); % hat d_i

%[id size(xi)]

end

for l_t = 2:L_t

xi(:,l_t) = A_xi*xi(:,l_t-1)+B1_xi*u0(l_t-1);

end

figure

hold on

plot(t,xi(5,:))

plot(t,xi(13,:))

legend('5-th output','13-th output')

The result is as below. The ratio of the amplitudes is clearly not -37.4dB. I also tried the frequency 1.11 rad/s (corresponding to the peak of bode amplitude), there is also a mismatch. I wonder what causes this mismatch?

Answer 1

[Paul]

Couldn't run the code becase sigma_d (and possibly other variables) were not defined.

Can you clarify what this code snippet is trying to do:

%Then I want to obtain the tranfer function between the 13th and 5th output as follows:
    [NUM_u,DEN_u] = ss2tf(A_xi,B1_xi,eye(length(A_xi)),zeros(length(A_xi),size(B1_xi,2)),1);

First of all, what is meant by: "the transfer function between the 13th and 5th output"? Was this supposed to mean "between the 13th input and the 5th output"? But that can't be the case, because (I think) B1_xi has only one column. So what are you actually trying to compute?

Presumably the outputs of the system are the state variables themselves?

The call to ss2tf doesn't seem correct. It has 6 inputs, when ss2tf takes at most 5 inputs.

[Hang Ruan]

Hi, Paul. Thank you for your generous help.

I have edited the code so that the variables are defined.

I am trying to get the transfer function between the 13th and 5th output according to

which has been done as equation (34) in this paper.

The number of inputs of ss2tf is 5, namely, A_xi, B1_xi, an identity matrix, a zero matrix and the index of input 1.

[Paul]

My mistake on the ss2tf. I miscounted the arguments.

I'm going to assume that the inputs to the ss2tf() call define the plant for a discrete time system with sample time Ts.

Maybe somehting like this is what you're looking for?

% vehicle dynimcs
    h = 0.8; % time headway
    N = 4; % number of vehicles 0 leading
    tau = 0.2; % actuator lag
    Ts = 0.1;% sampling interval
    % control gain
    kp = 0.45;
    kv = 0.84;
    ka = 0.5;
    % vehicle dynamics
    % x_i = [e_i;dv_i;a_i], ei=d_i-h*v_i, dvi = v_i-1-v_i
    A1 = [0,1,-h;0 0 -1; 0 0 -1/tau]; % x_i to x_i
    A2 = [0 0 0; 0 0 1;0 0 0]; % x_i-1 to x_i
    B = [0; 0; 1/tau]; % u_i to x_i
    % x0 = [q0;v0;a0]
    A0 = [0 1 0;0 0 1; 0 0 -1/tau]; % x0 to x0
    B0 = B; % u0 to x0
    Ac = zeros(3*N,3*N);
    Ac(1:3,1:3) = A0;
    Ac(5,3) = 1;
    for n = 1:N-1
        Ac(3*n+1:3*n+3,3*n+1:3*n+3) = A1;
        if n < N-1
            Ac(3*n+4:3*n+6,3*n+1:3*n+3) = A2;
        end
    end
    
    Bc = zeros(3*N,N);
    Bc(1:3,1) = B0;
    for n = 1:N-1
        Bc(3*n+1:3*n+3,n+1) = B;
    end
    
    veh_dyn_c = ss(Ac,Bc,eye(3*N),zeros(3*N,N));
    veh_dyn = c2d(veh_dyn_c,Ts);
    
    % kalman filter
    % x = [d,dv,da]
    Ak = [0 1 0; 0 0 1; 0 0 -1/tau];
    Ck = [0; 0; -1/tau]; % ui to x known input
    Bk = [0; 0; 1/tau]; % ui-1 to x disturbance
    H = [1 0 0; 0 1 0];
    sysk_c = ss(Ak, [Ck, Bk], H, zeros(2,2));
    sysk = c2d(sysk_c,Ts);
    delta_d = 0.5; % resolution
    delta_v = 0.6;%0.6;%
    SNR_db = 20;
    SNR = 10^(SNR_db/10);
    sigma_d = delta_d*sqrt(3/(2*pi^2*SNR));% standard deviation
    sigma_v = delta_v*sqrt(3/(2*pi^2*SNR));% standard deviation (don't forget to square)
    var_jerk = 1;%0.15499;% process noise variance (jerk)
    [KEST,~,~,L0,~,~] = kalman(sysk,var_jerk, diag([sigma_d^2,sigma_v^2])); % L: innovation gain
    L = L0;
    
    % whole dynamics
    % xi = [x0;x_1;d_1;u_1^-;\hat_{x_1^-};...;x_N-1;d_N-1;u_N-1^-;\hat_{x_N-1^-}]
    A_all = zeros(3+8*(N-1),3+8*(N-1)); % xi to xi^+
    A_f_all = zeros(3+8*(N-1),3+8*(N-1)); % xi^+ to xi^+
    B_all = zeros(3+8*(N-1),1); % u0 to to xi^+
    C_all = zeros(3+8*(N-1),2*(N-1)); % epsilon = [epsilon_d1;epsilon_v1;...;epsilon_dN-1;epsilon_vN-1] to xi^+
    % x0
    A_all(1:3,1:3) = veh_dyn.A(1:3,1:3); % x0 to x0+
    B_all(1:3) = veh_dyn.B(1:3,1); % u0 to x0+
    % x_i
    A_all(4:6,1:3) = veh_dyn.A(4:6,1:3); % x_0 to x_1^+
    B_all(4:6) = veh_dyn.B(4:6,1); % u0 to x_1^+
    for id1 = 1:N-1
        for id2 = 1:N-1
        A_all(8*id1-4:8*id1-2,8*id2-4:8*id2-2) = veh_dyn.A(3*id1+1:3*id1+3,3*id2+1:3*id2+3); % x_id2 to x_id1^+
        A_f_all(8*id1-4:8*id1-2,8*id2) = veh_dyn.B(3*id1+1:3*id1+3,id2+1); % u_id2 to x_id1^+
        end
    end
    % u_i
    for id = 1:N-1
        A_all(8*id,8*id-4) = kp; % e_i to u_i
        A_all(8*id,8*id-3) = kv; % dv_i to u_i
        A_all(8*id,8*id-2) = ka; % a_i to u_i
        A_f_all(8*id,8*id+3) = ka; % \hat_da_i to u_i
        C_all(8*id,2*id-1) = kp; % epsilon_d_i to u_i
        C_all(8*id,2*id) = kv; % epsilon_v_i to u_i
    end
    % d_i
    for id = 1:N-1
        A_f_all(8*id-1,8*id-4) = 1; % e_i to d_i
        A_f_all(8*id-1,2) = h; % v0 to d_i
        for id2 = 1:id
            A_f_all(8*id-1,8*id2-3) = -h;
        end
    end
    % \hat{x}_i
    for id = 1:N-1
        A_all(8*id+1:8*id+3, 8*id-1) = L(:,1); % d_i to \hat{x}_i
        A_all(8*id+1:8*id+3, 8*id-3) = L(:,2); % dv_i to \hat{x}_i
        C_all(8*id+1:8*id+3, 2*id-1) = L(:,1); % epsilon_d_i to \hat{x}_i
        C_all(8*id+1:8*id+3, 2*id) = L(:,2); % epsilon_v_i to \hat{x}_i
        A_all(8*id+1:8*id+3, 8*id+1:8*id+3) = (eye(3)-L*H)*sysk.A; % \hat{x}_i^- to \hat{x}_i
        A_all(8*id+1:8*id+3, 8*id) = (eye(3)-L*H)*sysk.B(:,1); % u_i^- to \hat{x}_i
    end
    
    % feedback
    A_xi = (eye(length(A_f_all))-A_f_all)\A_all;
    B1_xi = (eye(length(A_f_all))-A_f_all)\B_all; % u0
    B2_xi = (eye(length(A_f_all))-A_f_all)\C_all; % epsilon
% %Then I want to obtain the tranfer function between the 13th and 5th output as follows:
%     [NUM_u,DEN_u] = ss2tf(A_xi,B1_xi,eye(length(A_xi)),zeros(length(A_xi),size(B1_xi,2)),1);
%     Gamma_v21 = tf(NUM_u(13,:),NUM_u(5,:),Ts); % v2/v1
%     figure
%     bode(Gamma_v21)
    
% Obtain the ratio Y13(s)/Y5(s) when the systme is responding to an input
% U(s). The output Y is the state vector. Also, size() is better to use
% than length() for this use.

% get the ss model
    sys = ss(A_xi,B1_xi,eye(size(A_xi)),zeros(size(A_xi,1),size(B1_xi,2)),Ts);
    % zpk form of the transfer function to Y13 and h5
    h13 = zpk(sys(13,1));
    h5  = zpk(sys(5,1));
    % take the ratio
    R = h13/h5

R =
 
                                                                                                                                                                                                               
  0.08138 z^6 (z+0.8491) (z+0.8467) (z-0.8869) (z-0.8884) (z-0.9018) (z-1)^2 (z-0.7467) (z-0.6065) (z-0.03409) (z^2 - 1.757z + 0.7719) (z^2 - 1.834z + 0.8409) (z^2 - 1.807z + 0.8167) (z^2 - 0.2432z + 0.101) 
                                                                                                                                                                                                               
                                                                                                                                                                          (z^2 - 0.3533z + 0.1294)             
                                                                                                                                                                                                               
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                                                              
                   z^6 (z-0.6065) (z-0.7467)^2 (z-0.8869)^2 (z-0.9018) (z-1)^2 (z+0.8467) (z^2 - 1.757z + 0.7719) (z^2 - 1.834z + 0.8409)^2 (z^2 - 0.2432z + 0.101) (z^2 - 0.3533z + 0.1294)^2
                                                                                                                                                                                              
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

    % clean up the answer by cancelling common pole/zero pairs
    R1 = minreal(R)

R1 =
 
  0.08138 (z+0.8491) (z-0.8884) (z-0.7467) (z-0.03409) (z^2 - 1.834z + 0.8409) (z^2 - 1.807z + 0.8167)
  ----------------------------------------------------------------------------------------------------
               (z-0.7467)^2 (z-0.8869) (z^2 - 1.834z + 0.8409)^2 (z^2 - 0.3533z + 0.1294)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

It appears that additional reduction can be achieved by using a larger tolerance if believed to be appropriate, but this would need to carefully considered.

Verify R1 at 0.1 rad/sec

 L_t = 5000;
 t = (0:L_t)*Ts;
 y = lsim(sys,sin(0.1*t),t);
 max(y(:,13))/max(y(:,5))

ans = 0.9957

 m = bode(R1,0.1)

m = 0.9957

[Hang Ruan]

Update: The below mentioned problem has solved. It is caused by the error of simulated magnitude-frequency response due to transient state. Once I use the following code to get the simulated magnitude-frequency response, the simulated response fits the Bode plot well.

for id = 1:length(omega)

y = lsim(sys,sin(omega(id)*t),t);

m_sim(id) = max(y(200:L_t,13))/max(y(200:L_t,5));

end

------------------------------------------------------------------------------------------------------------------------------------

Thank you very much for your generous help! It helps a lot.

May I bother you with one more question. I try your method and test more frequencies. However, the magnitudes of the simulated response and Bode plot are different in higher frequencies. This still happens when I set other values of TOL for function minreal. I wonder if you have some idea on how to deal with it.

Test with different frequencies:

% get the ss model
sys = ss(A_xi,B1_xi,eye(size(A_xi)),zeros(size(A_xi,1),size(B1_xi,2)),Ts);
% zpk form of the transfer function to Y13 and h5
h13 = zpk(sys(13,1));
h5  = zpk(sys(5,1));
% take the ratio
R = h13/h5;
% clean up the answer by cancelling common pole/zero pairs
R1 = minreal(R);
TOL = 1e-3;
R2 = minreal(R,TOL);
L_t = 5000;
t = (0:L_t)*Ts;
omega = 10.^(-2:0.1:1);
m_sim = zeros(length(omega),1);
for id = 1:length(omega)
    y = lsim(sys,sin(omega(id)*t),t);
    m_sim(id) = max(y(:,13))/max(y(:,5));
end
m_bode = squeeze(bode(R1,omega));
m_bode2 = squeeze(bode(R2,omega));
figure
hold on
plot(omega,20*log10(m_sim));
plot(omega,20*log10(m_bode));
plot(omega,20*log10(m_bode2));
xlabel('$\omega$ (rad/s)','interpreter','latex')
set(gca,'xscale','log')
ylabel('Amplitude (dB)')
legend('Simulated','Bode',['Bode (Tol = ',num2str(TOL), ')'])

The result is as below

[Paul]

Yes, I was a bit loose there. The comparison needs to be done with the steady state output of the system after the transient response has (sufficiently) decayed.

Also, if you look carefully at the plot, you'll see that the low frequency gain is not quite the same between R1 and R2. That might be because the tol input into minreal is allowing some pole/zero cancellations that really shouldn't be. That's the type of thing I was referring to in my statement about "carefully considered."

One thing I did learn here is that the zpk() computation already does some type of reduction. After all, sys.a has 27 states, but nether h13 nor h5 have 27th order denominators, and, IIRC, they both have different order denominators. Maybe the zpk() calls are doing an sminreal() as part of that process?