Calculating inverse laplace transform using python or matlab

Question

I'm trying to simulate simple closed-loop system with PID controller in Python or MATLAB. In both cases I have problems with calculating time domain response of the system using inverse Laplace transform.

To ilustrate the problem better, I'm following the names on the picture below:

Depending on the transfer function of the system the result either calculates only in MATLAB or doesn't calculate at all.

Here is what I have:

syms s t k_p T_i T_d D_d v real;
% controller transfer function
C_s = k_p * (1 + 1 / (T_i * s) + (T_d * s) / (D_d * s + 1));

% plant transfer function
G_s = 1 / ((s + 1) * (s + 2));

% closed loop transfer function
W_s = C_s  * G_s / (1 + C_s  * G_s);

% input signal
R_t = 4 * heaviside(t - 2) * (1 - heaviside(t - 5)) + ...
      6 * heaviside(t - 5) * (1 - heaviside(t - 6)) + ...
      2 * heaviside(t - 6);

R_s = laplace(R_t, t, s, noconds=True)

% PID-controller system response
C_res_s = C_s * R_s;
w_t = ilaplace(C_res_s , s, t)

% Closed loop system response
Y_s = W_s * R_s
y_t_closed_loop = ilaplace(Y_s, s, t)

% Parameters
k_p_val = 15;
T_i_val = 5;
T_d_val = 2;
D_d_val = 0.1;
time = linspace(0, 10, 1000);

% Convert symbolic expressions to MATLAB functions
w_t_param = matlabFunction(w_t, 'Vars', {t, k_p, T_i, T_d, D_d});

matlabFunction(y_t_closed_loop, 'File', 'y_t_closed_loop_func', 'Vars', {t, k_p, T_i, T_d, D_d});
w_t_param_closed_loop = str2func('y_t_closed_loop_func');

R_t_lam = matlabFunction(R_t, 'Vars', {t});

% Calculate system responses
solution = zeros(size(time));
solution_closed_loop = zeros(size(time));
for i = 1:length(time)
    solution(i) = w_t_param(time(i), k_p_val, T_i_val, T_d_val, D_d_val);
    solution_closed_loop(i) = w_t_param_closed_loop(time(i), k_p_val, T_i_val, T_d_val, D_d_val);
end

For G_s = 1 / ((s + 1) * (s + 2)) I have (seems fine):

For G_s = 1 / ((s + 1)) I have (doesn't seem fine):

For G_s = 1 / ((s + 1) * (s + 2) * (s + 3)) I have nothing since the inverse laplace transform never calculates.

Here is the same code in Python:

s, t, k_p, T_i, T_d, D_d, v = symbols("s t k_p T_i T_d D_d v", real=True)

#controller transfer function
C_s = k_p * (1 + 1 / (T_i * s) + (T_d * s) / (D_d * s + 1))

#plant transfer function
G_s = 1 / ((s + 1) * (s + 2))

#closed loop transfer function
W_s = C_s  * G_s / (1 + C_s  * G_s)

#input signal
R_t = 4 * Heaviside(t - 2) * (1 - Heaviside(t - 5)) + 6 * Heaviside(t - 5) * (1 - Heaviside(t - 6)) + 2 * Heaviside(t - 6)
R_s = laplace_transform(R_t, t, s, noconds=True)

# PID-controller system response
C_res_s = C_s * R_s
w_t = inverse_laplace_transform(C_res_s , s, t)

# Closed loop system response
Y_s = W_s * R_s
y_t_closed_loop = inverse_laplace_transform(Y_s, s, t)

Python code falls into endless loop and never finds a solution to y_t_closed_loop.

Additionaly, for some reason the MATLAB code doesn't run on my machine with Windows. It only works on linux.

Am I missing something? Is there any way I can be sure that inverse laplace transform is calculated? The transfer funstions seem pretty normal to me and yet I have a feeling that something is wrong.

Any help very appreciated!!

EDIT: I also tried using sympy and control library in python, hoping it will change something:

import sympy
import numpy
import matplotlib.pyplot as plt
from tbcontrol.loops import feedback

s = sympy.Symbol('s')
t = sympy.Symbol('t', positive=True)
tau = sympy.Symbol('tau', positive=True)


K_p = sympy.Symbol('K_p')
T_i = sympy.Symbol('T_i')
T_d = sympy.Symbol('T_d')
D_d = sympy.Symbol('D_d')

G_p = 1/(s+1)
G_c = K_p * (1 + 1 / (T_i * s) + (T_d * s) / (D_d * s + 1))
G_OL = G_p*G_c
G_CL = feedback(G_OL, 1).cancel()
general_timeresponse = sympy.inverse_laplace_transform(sympy.simplify(G_CL/s), s, t)

As before, general_timeresponse never executes.

Answer 1

NOTE: I can't share my code, but I can give you the steps to reproduce it (which requires a little bit of work). The advantage of this procedure, which uses sympy's inverse_laplace_transform, is that you don't have to deal with Pade approximation of time delays. The main disadvantage is that you need to spend some time coding and testing it.

1. Once you have your output, Y_s, insert the numerical values:

   sd = { k_p: 15, T_i: 5, T_d: 2, D_d: 0.1 }
   tmp = Y_s.subs(sd).nsimplify().simplify().expand().together()
   tmp
   

   

   Note the last time delay, exp(-6*s), it is rendered on the numerator, but it is actually located on the denominator of the expression.

2. Extract the numerator and denominator. Time delays should be all on the numerator:

   n, d = fraction(tmp)
   n = (n * exp(-6*s)).expand()
   d = d / exp(6*s)
   display(n, d)
   

   

3. Notice the numerator is an addition. We can apply linearity and works on the single addends:

   args = [a / d for a in n.args]
   display(args)
   

   

4. This is where you have to put in your work. We need to create a function, func(expr), that does the following steps:

   i. Given an addend, func(a), look if a contains a time delay. If it does, remove it. (You can use delay = a.find(exp).pop()....

   ii. Extract the numerator and denominator: n, d = fraction(a).

   iii. Extract the coefficients from the numerator and denominator. For example, for the numerator: cn = Poly(n, s).all_coeffs().

   iv. Use the residue function from scipy.signal, which is going to create a numerical partial fraction expansion.

   v. Build a new symbolic expression from the residues, new_expr. This is likely the hardest step.

   vi. Compute the inverse laplace transform of new_expr. Be sure to use new_expr.nsimplify(). If you don't, the computation might takes much longer.

   vii. If a time delay was found on point i., apply the time delay with out.subs(t, t+time_delay), where time_delay=delay.args[0]/s.

   viii. Return the result.

5. Loops over the addends and run the function created in step 3. outs = [func(a) for a in args].

6. Create the output signal, and plot it (I'm using SymPy Plotting Backed:

   from spb import plot
   out = sum(outs)
   plot((out, "closed loop res"), (R_t, "R_t"), (t, -1, 10))