Plotting script of Digital Control System in NTU. Mainly with SciPy and Matplotlib.

Project 1.py

from scipy import signal
import matplotlib.pyplot as plt
import math
import numpy as np
plt.rc('text', usetex=True)
W_RANGE = np.arange(-4, 6, 0.001).tolist()
W_RANGE = [10**i for i in W_RANGE]

def plot_bode(sys, legend, tittle = ''):
    fig, axs = plt.subplots(2, 1, constrained_layout=True)
    fig.suptitle(tittle)
    n = len(sys)
    phase_plot = [0] * n
    mag_plot = [0] * n
    for i in range(n):
        w, mag, phase = signal.bode(sys[i], W_RANGE)
        mag = [10.0**(data / 10.0) for data in mag]
        w = [data / (2.0 * math.pi) for data in w]

        mag_plot[i] = axs[0].semilogx(w, mag, label=legend[i])
        axs[0].axvline(x=50, color='k', linestyle='--')
        axs[0].set_title('Gain')
        axs[0].set_xlabel('Hz')

        phase_plot[i] = axs[1].semilogx(w, phase, label=legend[i])
        axs[1].axvline(x=50, color='k', linestyle='--')
        axs[1].set_title('Phase')
        axs[1].set_xlabel('Hz')
        axs[1].set_ylabel('degree')
    axs[1].legend(loc=3)

# ProblemA
J = 0.01
B = 3e-2
#BW = 50.0 * 2 * math.pi
BW = 50.0 * 2.0 * math.pi
W_max = 300
f = 100.0

## Problem A1
PI_ctrl = signal.lti([BW], [1.0, BW])
plot_bode([PI_ctrl], ['PI'], 'Fig A1. Bode Plot of PI-Controlled System')

## Problem A2
t = np.linspace(0, 0.1, num=100000)
u = (W_max / 2.0) * (-1.0 * np.cos((f * 2 * math.pi) * t) + 1.0)

tout, y, x = signal.lsim(system=PI_ctrl, U=u, T=t)
plt.figure()
plt.xlabel("time(s)")
plt.ylabel("magnitude")
plt.title("Fig A2. Time Response for PI-Controlled System")
plt.plot(t, u, linestyle='-.')
plt.plot(t, y)

## Problem A3
def CFF(J_hat, B_hat):
    num = [J_hat, BW * J + B_hat, BW * B]
    denum = [J, B + BW * J, BW * B]
    return signal.lti(num, denum)
CFF_ctrl = CFF(J, B)
_, y_cff, _ = signal.lsim(system=CFF_ctrl, U=u, T=t)
plt.figure()
plt.xlabel("time(s)")
plt.ylabel("magnitude")
plt.title("Fig A3. Time Response for PI-Controlled System and CFF-Controlled System")
plt.plot(t, u, linestyle='-.')
plt.plot(t, y)
plt.plot(t, y_cff)

## Problem A4
plot_bode(\
[PI_ctrl, CFF(J, B), CFF(2.0*J, 1.0*B), CFF(0.5*J, 1.0*B)],\
[r'$PI without CFF$', '$CFF, \hat{J} = J$', r'$CFF, \hat{J} = 2J$', r'$CFF, \hat{J} = 0.5J$'], \
"Fig A4. Bode Plots w.r.t. Different J")

## Problem A5
plot_bode (\
[PI_ctrl, CFF(J, B), CFF(J, 2.0*B), CFF(J, 0.5*B), CFF(0.5*J, B), CFF(2*J, B)], \
[r'$PI without CFF$', \
r"$CFF, \hat{B} = B$", \
r"$CFF, \hat{B} = 2B$", \
r"$CFF, \hat{B} = 0.5B$", \
r"$CFF, \hat{J} = 2J$",\
r"$CFF, \hat{J} = 0.5J$",\
], \
"Fig A5. Bode Plots w.r.t. Different Bs and Js")

# Probrlem A6
def A6(J_hat, B_hat):
    return signal.lti([J_hat, J * BW + B_hat], [J, J * BW + B]);
plot_bode (\
[PI_ctrl, A6(J, B), A6(J, 2.0*B), A6(J, 0.5*B), A6(0.5*J, B), A6(2*J, B)], \
[r'$PIi\ without\ CFF$', \
r"$CFF, \hat{B} = B$", \
r"$CFF, \hat{B} = 2B$", \
r"$CFF, \hat{B} = 0.5B$", \
r"$CFF, \hat{J} = 2J$",\
r"$CFF, \hat{J} = 0.5J$",\
], \
"Fig A6. Bode Plots w.r.t. Different Bs and Js")

# Problem A Bonus
sys1 = CFF(0.5*J, B)
sys2 = A6(0.5*J, B)
t = np.linspace(0, 0.1, num=10000)
plt.figure()
plt.xlabel("time(s)")
plt.ylabel("magnitude")
plt.title("Fig Bonus A. ")
t1, y1 = sys1.step(T=t)
plt.plot(t1, y1)
t2, y2 = sys1.step(T=t)
plt.plot(t2, y2)

# Problem B
K_p = BW * J
K_i = BW * B
t = np.linspace(0, 0.1, num=10000)
SIGFREQHZ = 60.0
u = -0.2 * signal.square((SIGFREQHZ * 2.0 * math.pi) * t)

## Problem B1
P_disturb = signal.lti([1], [J, B + K_p])
_, y_p, _ = signal.lsim(system=P_disturb, U=u, T=t)
plt.figure()
plt.xlabel("time(s)")
plt.ylabel("magnitude")
plt.title("Fig B1. Time Response for P-Controlled System with AC Torque Disturbance")
plt.plot(t, u)
plt.plot(t, y_p)

## Problem B2
PI_disturb = signal.lti([1, 0], [J, B + K_p, K_i])
_, y_pi, _ = signal.lsim(system=PI_disturb, U=u, T=t)
plt.figure()
plt.xlabel("time(s)")
plt.ylabel("magnitude")
plt.title("Fig B2. Time Response for PI-Controlled System with AC Torque Disturbance")
plt.plot(t, u)
plt.plot(t, y_pi)

# Problem B3
plot_bode([P_disturb, PI_disturb], \
['P', 'PI'],\
"Fig B3. Disturbance Rejection of P and PI Controllers")

# Problem B4
with_did = signal.lti([0.0, 0.0], [10000000])
without_did = signal.lti([1.0], [J, B])
_, y_no_did, _ = signal.lsim(system=without_did, U=u, T=t)
_, y_did, _ = signal.lsim(system=with_did, U=u, T=t)
plt.figure()
plt.xlabel("time(s)")
plt.ylabel("magnitude")
plt.title("Fig B4 before and after did with zero-command.")
plt.plot(t, u)
plt.plot(t, y_no_did)
plt.plot(t, y_did)


# bonus
t = np.linspace(0, 0.1, num=10000)
SIGFREQHZ = 60.0
u = -0.2 * signal.square((SIGFREQHZ * 2.0 * math.pi) * t)
ignore_regulator = signal.lti([3.0e-4, 0], [3.0e-6, 0.01 + 3e-6, 3e-2])
_, y_ignore, _ = signal.lsim(system=ignore_regulator, U=u, T=t)
plt.figure()
plt.xlabel("time(s)")
plt.ylabel("magnitude")
plt.title("Fig B6. time response ignoring regulator, with disturbance.")
plt.plot(t, u)
plt.plot(t, y_ignore)


plt.show()