How to Obtain State Space Model from Bode Plot Data

Hi Everyone,

I have Experiemental Data (Bode Plot), using MATLAB I found the Transfer Function of my plant. I have attached my E7i_CSV file.

T2 = readtable('E7i_CSV.csv');   
                                                                       % Fill With Actual Sampling Frequency
FHz = T2.F;
Ts = 1/(2*(max(FHz)+10000))
for k = 1:size(T2,1)-1
    if FHz(k+1) == FHz(k)
        FHz(k+1) = FHz(k+1)+0.5;                                                % 'Brute Force' Interpolation 
    end
end
Mag = T2.G;
PhDeg = T2.P;
Response = Mag.*exp(1j*deg2rad(PhDeg));                                         % Complex Vector
sysfr = idfrd(Response, FHz, Ts, 'FrequencyUnit','Hz')
%bode(sysfr)
%tfsys = tfest(sysfr,4,3) This is much better
tfsys = tfest(sysfr,4,3) %82.69%
%tfsys = tfest(sysfr,6,3)
figure
compare(sysfr, tfsys)
is_stable = isstable(tfsys);
disp(is_stable);

I want to convert this transfer function to statespace equations, which will be used for Model Predictive Control Design. Simple tf2ss command give me TF but it doesn't look very accrurate.

 %Assuming you have a tfsys transfer function object named 'tf'
num = tfsys.Numerator;  % Extract numerator coefficients
den = tfsys.Denominator;  % Extract denominator coefficients
[A, B, C, D] = tf2ss(num,den)

So, I was reading a paper which i will follow (also attached), they have mentioned that they had FRF (Frequency Response Function) data and they find state space model/equations using Prediction Error Method (PEM). But I am unable to understand how I can used Bode Plot data and convert it state space using PEM. I will appreciate your help.

Secondy I verified the paper StateSpace Model by first converting it to Transfer Function and then using tf2ss command. which shows inaccurate results similiar to my results. I am looking for your help. Following is the verification code for paper i am following

%Rana Paper State Space Equations Coversion to
                            %Transfer Function for Validation Reason
% Define state-space matrices in the Paper 
Ax = [-49.3277, -980.7648, 618.2198;
      567.3557, -120.2273, 651.7342;
      -45.8885, -415.4194, -115.7522];

Bx = [-1.8122; 3.901370; 4.3011];

Cx = [14.9368, 16.9208, -23.6827];

Dx = 0;

% Create a state-space model
sys = ss(Ax, Bx, Cx, Dx);

% Display the state-space model
disp('State-Space Model:');
disp(sys);
% Convert the state-space model to a transfer function
tf_sys = tf(sys);

% Display the transfer function
disp('Transfer Function:');
disp(tf_sys);
                        %TF to StateSpace Rana Model
num_rana = [ -62.92 3.625*10^4 -1.065*10^8];
den_rana = [ 1 285.3 8.811*10^5 1.982*10^8];
tff_rana = tf(num_rana,den_rana);
tffss_rana = tf2ss(num_rana,den_rana);
[Arana,Brana,Crana,Drana] = tf2ss(num_rana,den_rana);

Answer 1

[Paul]

Hi Muhammad,

Here is the first part of the code:

T2 = readtable('E7i_CSV.csv');   
                                                                       % Fill With Actual Sampling Frequency
FHz = T2.F;
Ts = 1/(2*(max(FHz)+10000))

Ts = 1.6667e-05

for k = 1:size(T2,1)-1
    if FHz(k+1) == FHz(k)
        FHz(k+1) = FHz(k+1)+0.5;                                                % 'Brute Force' Interpolation 
    end
end
Mag = T2.G;
PhDeg = T2.P;
Response = Mag.*exp(1j*deg2rad(PhDeg));                                         % Complex Vector
sysfr = idfrd(Response, FHz, Ts, 'FrequencyUnit','Hz')

sysfr =

IDFRD model.
Contains Frequency Response Data for 1 output(s) and 1 input(s).
Response data is available at 1001 frequency points, ranging from 100 Hz to 2e+04 Hz.
 
Sample time: 1.6667e-05 seconds
Status:                                                         
Created by direct construction or transformation. Not estimated.
 

%bode(sysfr)
%tfsys = tfest(sysfr,4,3) This is much better
tfsys = tfest(sysfr,4,3) %82.69%

tfsys =
     -6.928 s^3 + 5.189e05 s^2 - 3.755e10 s + 3.039e15
  -------------------------------------------------------
  s^4 + 3.085e04 s^3 + 5.188e09 s^2 + 1.2e14 s + 4.462e18
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 4   Number of zeros: 3
   Number of free coefficients: 8
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                  
Estimated using TFEST on frequency response data "sysfr".
Fit to estimation data: 79.61%                           
FPE: 1.574e-08, MSE: 1.549e-08                           
 

%tfsys = tfest(sysfr,6,3)
figure
compare(sysfr, tfsys)

is_stable = isstable(tfsys);
disp(is_stable);

   1

num = tfsys.Numerator;  % Extract numerator coefficients
den = tfsys.Denominator;  % Extract denominator coefficients
[A, B, C, D] = tf2ss(num,den);

Here, make a new plot adding the frequency response of the state space model. As expected, it exactly (to withing numerical error) overlays the response of the transfer function model.

figure
compare(sysfr, tfsys, ss(A,B,C,D))

Before going any further, can you explain why "it doesn't look very accrurate."?

[Muhammad]

Dear Paul,

Thanks for answering to my Question, let me explain you, If you see my Matric A, only A13 and A14 have numerical values, all other values are 0 in matric A. It made me to worry about it because I want to design kalman Estimator and Model Predictive Control based on that experimental data and state space model.

Then I wanted to Evaluate the Paper Matrices as I explained in my previous post, >>>>> Secondly I verified the paper StateSpace Model by first converting it to Transfer Function and then using tf2ss command. which shows inaccurate results similiar to my results. I am looking for your help. Following is the verification code for paper i am following.

The above screenshot has paper A,B,C,D matrices, and when I use following code to convert it back to the state-space form to evaluate my Arana (please picture attached below code, and read last line code) matric doesn't match with Ax.

%Rana Paper State Space Equations Coversion to
                            %Transfer Function for Validation Reason
% Define state-space matrices in the Paper 
Ax = [-49.3277, -980.7648, 618.2198;
      567.3557, -120.2273, 651.7342;
      -45.8885, -415.4194, -115.7522];
Bx = [-1.8122; 3.901370; 4.3011];
Cx = [14.9368, 16.9208, -23.6827];
Dx = 0;
% Create a state-space model
sys = ss(Ax, Bx, Cx, Dx);
% Display the state-space model
disp('State-Space Model:');
disp(sys);
% Convert the state-space model to a transfer function
tf_sys = tf(sys);
% Display the transfer function
disp('Transfer Function:');
disp(tf_sys);
                        %TF to StateSpace Rana Model
num_rana = [ -62.92 3.625*10^4 -1.065*10^8];
den_rana = [ 1 285.3 8.811*10^5 1.982*10^8];
tff_rana = tf(num_rana,den_rana);
tffss_rana = tf2ss(num_rana,den_rana);
[Arana,Brana,Crana,Drana] = tf2ss(num_rana,den_rana)

Please Help me I am very stressed.

[Paul]

A linear system represented by a transfer function can also be represented by an infinite number of state space models that have the exact same input-output response. However, tf2ss can only return one of those infinite possibilities. So tf2ss chooses to always return the result in a canonical form. Notice that, in both cases shown above, the A and B matrices have a specific form. The first element of the B-matrix is exactly one, and the others are exactly zero. The first row of the A-matrix are the negatives of the coefficients of the transfer function denominator (assuming the leading coefficient is 1), and the first lower diagonal of A are all 1's. Also, compare the C-matrix to the numerator coefficients.

Ax = [-49.3277, -980.7648, 618.2198;
      567.3557, -120.2273, 651.7342;
      -45.8885, -415.4194, -115.7522];
Bx = [-1.8122; 3.901370; 4.3011];
Cx = [14.9368, 16.9208, -23.6827];
Dx = 0;
% Create a state-space model
sys1 = ss(Ax, Bx, Cx, Dx);

sys2 = tf(sys1)   % convert to transfer function

sys2 =
 
    -62.92 s^2 + 3.625e04 s - 1.065e08
  ---------------------------------------
  s^3 + 285.3 s^2 + 8.811e05 s + 1.982e08
 
Continuous-time transfer function.

format short e
num_rana = sys2.Numerator{1}

num_rana = 1×4

1.0e+00 * 0 -6.2916e+01 3.6251e+04 -1.0646e+08

den_rana = sys2.Denominator{1}

den_rana = 1×4

1.0e+00 * 1.0000e+00 2.8531e+02 8.8111e+05 1.9824e+08

Note the pattern described above.

format short e
[A,B,C,D] = tf2ss(num_rana,den_rana)

A = 3×3

1.0e+00 * -2.8531e+02 -8.8111e+05 -1.9824e+08 1.0000e+00 0 0 0 1.0000e+00 0

B = 3×1

1 0 0

C = 1×3

1.0e+00 * -6.2916e+01 3.6251e+04 -1.0646e+08

D = 
     0

sys3 = ss(A,B,C,D);

Identical responses

bode(sys1,sys2,sys3)

If you are starting with a state space model for which you want to design a Kalman Filter, then there shouldn't be a need to convert to transfer function and then back to state space. I'm not familiar with MPC, so can't really comment on that.

[Muhammad]

Thank You Paul for helping me. I got it that even the values came out to be different but when we compare them by making bode plot each state space model has similiar frequency response.

Stay Blessed

[Muhammad]

One thing more @Paul can you further convert above state-space model of sys3 = ss(A,B,C,D); to discrete state space please.

[Paul]

See c2d for appoximating a continuous time system with a discrete time model. That function offers a few options to apply for the approximation.