How can simulate The system with Transfer function in ODE45?

Hi guys, I have the trouble with solving the ODE with a trasfer function.

I would like to simulate using ode solver in matlab as following block diagram:

where and the plant is described by:

.

Kis the LQ gain.

Here is the my code.

clear all
close all

% Set Parameters for Simulation
cal_step = 0.001;
end_time = 200;
tspan = [0:cal_step:end_time];
x0 = [10;0];

%% Set Parameters of the a(s)%%
alpha = 100;
beta = 0.1;
up_num = [alpha 0];
under_num = [beta 1];
[A,B,C,D] = tf2ss(up_num, under_num);
b_s = tf(up_num, under_num);

%% Set Parameters of the Plant %%
m = 2; k = 0.5; c = 0.05;
A_plant = [0 1; -k/m -c/m];
B_plant = [0; 1/m];
C_plant = eye(2);
D_plant = [0;0];

%% Calculate the LQ Gain %%
Q = eye(2);
R = 1;
P_riccati = icare(A_plant, B_plant, eye(2), 1);
K_ipd = -inv(R)*transpose(B_plant)*P_riccati;

%% Make the Structure for the ode function %%
para.A = A;
para.B = B;
para.C = C;
para.D = D;
para.A_plant = A_plant;
para.B_plant = B_plant;
para.K_ipd = K_ipd;

%% Simulation ode45 and simulink %%
[t x_plant] = ode45(@(t,x_plant) odeplant(t, x_plant, para), tspan, x0);
sim('simulink_model')

simresult = ans;
x_simulink = simresult.x.Data;
t_simulink = simresult.x.Time;

%% Plot the Result %%
figure(1)
subplot(3,1,1)
plot(t_simulink, x_simulink)
title("Simulink")

subplot(3,1,2)
plot(tspan,x_plant)
title("M file Program")

subplot(3,1,3)
plot(tspan,x_plant-x_simulink)
title("Error between M file and Simulink")


function dxdt = odeplant(t,x,para)
  A_plant = para.A_plant;
  B_plant = para.B_plant;
  A = para.A;
  B = para.B;
  C = para.C;
  D = para.D;

  K_ipd = para.K_ipd;

  u = K_ipd*x;
  u_dash = C*integral(@(tau) exp(A*(t-tau))*B.*u, 0, t) + D.*u;  % Calculating the transfer function by convolution

  dxdt = A_plant*x+B_plant*u_dash;
end

This is the result figure:

In my code, I calculate the trasfer function as a ordinary differential equation.

Thank you!

Answer 1

[Sam Chak]

Hi @Kaito Sato-kun

Most likely the culprits come from this part

up_num = [alpha 0];

and this part

u_dash = C*integral(@(tau) exp(A*(t-tau))*B.*u, 0, t) + D.*u;

Anyhow, this is what I would probably do to analyze the system in continuous-time domain.

Let's begin with your function

which can be simplified to

alpha = 100;
beta  = 0.1;
num   = [alpha 1];
den   = [beta 1];
astf  = minreal(tf(num, den))

astf =
 
  1000 s + 10
  -----------
    s + 10
 
Continuous-time transfer function.

Note that the 1st-order transfer function

when converted to state-space form, it becomes

or

For some unknown reasons, since your P(s) transfer function is described by state-space

then we can combine both of them by first defining the state variables , , and we have

.

Thus, your original Q has to be modified to have the same size as the state matrix A, in order to determine the gain K from your icare() function

m = 2; 
k = 0.5; 
c = 0.05;
A = [0 1 0; -k/m -c/m 1/m; 0 0 -10];
B = [0; (1/m)*1000; -9990];
Q = [1 0 0; 0 1 0; 0 0 0];
R = 1;
K = lqr(A, B, Q, R)         % I prefer LQR though :)

K = 1×3

1.0005 0.9975 0.0008

% Solving the system with ode45
fv1 = @(t, x, y, z) y;
fv2 = @(t, x, y, z) - (k/m)*x - (c/m)*y + (1/m)*z + (1/m)*1000*(- K(1)*x - K(2)*y - K(3)*z);
fv3 = @(t, x, y, z) - 10*z - 9990*(- K(1)*x - K(2)*y - K(3)*z);
opt = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);
tspan = [0 20];
x0  = [10 0 0];
[t, v] = ode45(@(t, x) ([fv1(t, x(1), x(2), x(3)); fv2(t, x(1), x(2), x(3)); fv3(t, x(1), x(2), x(3))]), tspan, x0, opt);

% Plotting the result
plot(t, v(:,1:2), 'linewidth', 1.5)
grid on
xlabel({'$t$'}, 'Interpreter', 'latex')
ylabel({'$\mathbf{x}(t)$'}, 'Interpreter', 'latex')
legend({'$x(t)$', '$\dot{x}(t)$'}, 'Interpreter', 'latex', 'FontSize', 14)

[Paul]

Hi @Sam Chak

Once a, P, and K are known, why not use interconnection functions to get the ss realization of the closed-loop system and then implement odefun as xdot = A_closedloop * x.

[Sam Chak]

Yes @Paul, you're right! 👍 Somehow, I didn't think deeper about that. Thanks for highting that to me and OP.

[Kaito Sato]

@Sam Chak Thank you so much for your help!!! We can solve the transfer function like using not a convolution but ODE function, corresct? Therefore you create one state-space equation which includes the plant and .

Thanks a lot! I got it!

[Sam Chak]

You are welcome, @Kaito Sato-san. If you find the mini tutorial about simulating systems with mixed transfer function and ODE is helpful, please consider accepting ✔ and voting 👍 the Answer. Thanks!