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:
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?
