PID function equivalent from MATLAB to Python

I manage to correctly simulate a closed loop with Simulink using pure python "control" library. However, I'm now trying to do the same with a new closed loop that has uses a PID controller:

import numpy as np
    import control as ctl
    
    def nuclear_params(array=[0.0001889, 0, -1.289, 0.2891]):
    num = array[:2]
    den = [1.0] + array[2:]
    ref_start = 0
    ref_value = 6.6
    steps = 700
    timestep = 0.5
    Kp, Ki, Kd = 348.52, 17.25, 10.79
    return simulation_nuclear(
    Kp, Ki, Kd, num, den, ref_start, ref_value, timestep, steps
    )
    
    def simulation_nuclear(
    Kp, Ki, Kd, num_plant, den_plant, ref_start, ref_value, timestep=0.5, steps=700
    
    ):
    Ts = timestep
    
    Gc = pid_discrete(Kp, Ki, Kd, Ts, N=1000)
    Gp = ctl.TransferFunction(num_plant, den_plant, Ts)
    
    fdback = ctl.feedback(Gc * Gp, 1)
    
    # Time vector and input reference step
    t = np.linspace(0, steps, int(round(steps / Ts)) + 1)
    ref = np.ones_like(t) * ref_value
    
    # System response
    t_out, y_out = ctl.forced_response(fdback, T=t, U=ref)
    
    # Retrieve control output
    e = ref - y_out
    _, u_m = ctl.forced_response(Gc, T=t, U=e)
    
    return t_out, y_out, u_m

With u_m giving:

4610.50075588, 4022.47644969, 3334.11898869, 2705.64516918, 2175.79419241, 1743.04805239

The main problem now is that the control output is more or less DOUBLE of what Matlab gives stumbler number:

2327.93749436396, 2267.04211542563, 2134.18318770832, 1995.12935675866, 1854.49534374506,

This result is most likely the fault of PID approximation function pid_discrete which converts the continuous into discrete and yield the best result for now.

I have tried the following implementations without success:

def pid_continuous(Kp, Ki, Kd, N=1000):
Ti = Kp / Ki
Td = Kd / Kp
s = ctl.TransferFunction.s
Gc = Kp * (1 + 1 / (Ti * s) + (Td * N * s) / (1 + Td * N * s))
return Gc


def pid_discrete(Kp, Ki, Kd, Ts, N=1000, method="tustin"):
Gc_s = pid_continuous(Kp, Ki, Kd, N=N)
return ctl.sample_system(Gc_s, Ts, method=method)


def pid_forward_euler(Kp, Ki, Kd, Ts, N=1000):
z = ctl.TransferFunction.z
if z is None:
    raise ValueError("Cannot create z-domain transfer function")
    # Forward Euler integrator: IF(z) = Ts*z/(z-1)
    integrator = Ts * (z / (z - 1))
    # Forward Euler derivative with filter:
    # Derivative term: Kd * N * (z-1) / (z + N*Ts - 1)
    derivative = (Kd * N * (z - 1)) / (z + N * Ts - 1)
    # Parallel form: Kp + Ki*integrator + derivative
    pid_controller = Kp + Ki * integrator + derivative
    
    return pid_controller
    
    
    def pid_backward_euler(Kp, Ki, Kd, Ts, N=1000):
    z = ctl.TransferFunction.z
    if z is None:
        raise ValueError("Cannot create z-domain transfer function")
        # Backward Euler integrator: IF(z) = Ts/(z-1)
        integrator = Ts / (z - 1)
        # Backward Euler derivative with filter
        derivative = Kd * N * (z - 1) / (z - 1 + N * Ts)
        pid_controller = Kp + Ki * integrator + derivative
        return pid_controller
        
        
        def pid_trapezoidal(Kp, Ki, Kd, Ts, N=1000):
        z = ctl.TransferFunction.z
        if z is None:
            raise ValueError("Cannot create z-domain transfer function")
            
            # Trapezoidal integrator: IF(z) = Ts*(z+1)/(2*(z-1))
            integrator = Ts * (z + 1) / (2 * (z - 1))
            # Trapezoidal derivative with filter
            # This is more complex for trapezoidal - simplified version:
            derivative = Kd * 2 * N / (2 + N * Ts) * (z - 1) / z
            pid_controller = Kp + Ki * integrator + derivative
            return pid_controller

Except for pid_discrete, all of the other functions either explode or gives some error. My question is: there any working implementation of a PID function with python + control library?

[Sam Chak]

Hi @Fahad

Is the plant modeled in continuous time or discrete time? I have the same question regarding the PID controller. If possible, could you demonstrate how the control system behaves in MATLAB?

Ts  = 0.5;
num = [0.0001889, 0];
den = [1, -1.289, 0.2891];

%% continuous-time transfer function
Gp  = tf(num, den)

Gp =
 
       0.0001889 s
  ----------------------
  s^2 - 1.289 s + 0.2891
 
Continuous-time transfer function.

%% discrete-time transfer function
Gpd = tf(num, den, Ts)

Gpd =
 
       0.0001889 z
  ----------------------
  z^2 - 1.289 z + 0.2891
 
Sample time: 0.5 seconds
Discrete-time transfer function.