Nyquist Plot TIme Delayed¶

기출문제 풀이¶

20년 3교시 6번

In [1]:
import control as ctrl
import matplotlib.pyplot as plt
import numpy as np
import sympy
import math
In [2]:
s = sympy.Symbol('s')
ce = (s+1)*(s+1)*(s+1)
ce_polynomial = sympy.expand(ce)
ce_polynomial
Out[2]:
$\displaystyle s^{3} + 3 s^{2} + 3 s + 1$
In [3]:
num = [2*math.sqrt(2)]
den = [1, 3, 3, 1]
system = ctrl.TransferFunction(num, den)
system
Out[3]:
$$\frac{2.828}{s^3 + 3 s^2 + 3 s + 1}$$
In [4]:
# Create a Nyquist plot
ctrl.nyquist_plot(system, omega=np.logspace(-5, 5, 1000), plot=True)

# Customize the plot (optional)
plt.title('Nyquist Plot')
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.show()

time dealy 는 Pad´e Approximation으로
1-s tau/2
------------
1+s tau/2
로 근사

In [5]:
tau = 0.785
#numerator
numerator_delay = (1-s*tau/2)*2*math.sqrt(2)
numerator_delay
Out[5]:
$\displaystyle 2.82842712474619 - 1.11015764646288 s$
In [6]:
coef_num = []
coef_num.append(float(numerator_delay.coeff(s,0)))
coef_num.append(float(numerator_delay.coeff(s,1)))
In [7]:
#denominator
ce_delay = (s+1)*(s+1)*(s+1)*(1+s*tau/2)
ce_polynomial_delay = sympy.expand(ce_delay)
ce_polynomial_delay
Out[7]:
$\displaystyle 0.3925 s^{4} + 2.1775 s^{3} + 4.1775 s^{2} + 3.3925 s + 1$
In [8]:
coef = []
coef.append(float(ce_polynomial_delay.coeff(s,0)))
coef.append(float(ce_polynomial_delay.coeff(s,1)))
coef.append(float(ce_polynomial_delay.coeff(s,2)))
coef.append(float(ce_polynomial_delay.coeff(s,3)))
coef.append(float(ce_polynomial_delay.coeff(s,4)))
In [9]:
num_delay = [coef_num[1], coef_num[0]]
den_delay = [coef[4], coef[3], coef[2], coef[1], coef[0]]
system_delay = ctrl.TransferFunction(num_delay, den_delay)
system_delay
Out[9]:
$$\frac{-1.11 s + 2.828}{0.3925 s^4 + 2.178 s^3 + 4.178 s^2 + 3.393 s + 1}$$
In [10]:
# Create a Nyquist plot
ctrl.nyquist_plot(system_delay, omega=np.logspace(-5, 5, 1000), plot=True)

# Customize the plot (optional)
plt.title('Nyquist Plot Delayed')
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.show()
In [16]:
#-1 Nyquist Plot을 확대
print("tau: {}".format(tau))
ctrl.nyquist_plot(system_delay, omega=np.logspace(-5, 5, 1000), plot=True)
plt.xlim([-1.2, -0.8])      # X축의 범위: [xmin, xmax]
plt.ylim([-0.1, 0.1])     # Y축의 범위: [ymin, ymax]
plt.show()

tau: 1.0

원본 Nyquist와 시간 지연된 Nyquist를 한번에 표시¶

In [12]:
system
Out[12]:
$$\frac{2.828}{s^3 + 3 s^2 + 3 s + 1}$$
In [13]:
system_delay
Out[13]:
$$\frac{-1.11 s + 2.828}{0.3925 s^4 + 2.178 s^3 + 4.178 s^2 + 3.393 s + 1}$$
In [14]:
ctrl.nyquist_plot(system, omega=np.logspace(-5, 5, 1000), plot=True, color="Blue")
ctrl.nyquist_plot(system_delay, omega=np.logspace(-5, 5, 1000), plot=True, color ="Red")
# Customize the plot (optional)
plt.title('Nyquist Plot not Delayed and Delayed(Blue = original)')
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
print("time delayed: {}".format(tau))
plt.show()

time delayed: 0.785

tau를 0.1초 단위로 증가 시키면서 Nuquist plot

In [22]:
import matplotlib.patches as patches

for i in range (0,6):
    tau = i*0.2
    i = i + 1
    #numerator
    numerator_delay = (1-s*tau/2)*2*math.sqrt(2)
    coef_num = []
    coef_num.append(float(numerator_delay.coeff(s,0)))
    coef_num.append(float(numerator_delay.coeff(s,1)))
    #denominator
    ce_delay = (s+1)*(s+1)*(s+1)*(1+s*tau/2)
    ce_polynomial_delay = sympy.expand(ce_delay)
    coef = []
    coef.append(float(ce_polynomial_delay.coeff(s,0)))
    coef.append(float(ce_polynomial_delay.coeff(s,1)))
    coef.append(float(ce_polynomial_delay.coeff(s,2)))
    coef.append(float(ce_polynomial_delay.coeff(s,3)))
    coef.append(float(ce_polynomial_delay.coeff(s,4)))
    num_delay = [coef_num[1], coef_num[0]]
    den_delay = [coef[3], coef[2], coef[1], coef[0]]
    system_delay = ctrl.TransferFunction(num_delay, den_delay)
    ctrl.nyquist_plot(system_delay, omega=np.logspace(-5, 5, 1000), plot=True)
#plt.xlim([-1.5, 0.5])      # X축의 범위: [xmin, xmax]
#plt.ylim([-3, 0.5])     # Y축의 범위: [ymin, ymax]
#circle
#shp=patches.Circle((0,0), radius=1, fill = 0)
#plt.gca().add_patch(shp)