안정도 판별¶

1. Routh Array¶

In [1]:
!pip install tbcontrol

Requirement already satisfied: tbcontrol in /usr/local/lib/python3.9/dist-packages (0.2.1)
Requirement already satisfied: numpy in /usr/local/lib/python3.9/dist-packages (from tbcontrol) (1.22.3)
Requirement already satisfied: scipy in /usr/local/lib/python3.9/dist-packages (from tbcontrol) (1.8.1)
Requirement already satisfied: tqdm in /usr/local/lib/python3.9/dist-packages (from tbcontrol) (4.66.1)
Requirement already satisfied: packaging in /usr/local/lib/python3.9/dist-packages (from tbcontrol) (21.3)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /usr/local/lib/python3.9/dist-packages (from packaging->tbcontrol) (3.0.9)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
In [2]:
import sympy
from sympy import solve
from scipy.optimize import fsolve
import control as ct
import numpy as np
import matplotlib.pyplot as plt
from tbcontrol.symbolic import routh
In [3]:
s = sympy.Symbol('s')
K = sympy.Symbol('K')
G = K/(s*(s+1)*(2*s+1))
CLTF = G /(1+G)
CLTF = sympy.cancel(CLTF)
CLTF = sympy.expand(CLTF)
In [4]:
CLTF
Out[4]:
$\displaystyle \frac{K}{K + 2 s^{3} + 3 s^{2} + s}$
In [5]:
p = K + 1*s**1 + 3*s**2 + 2*s**3
p = sympy.Poly(p, s)
p
Out[5]:
$\displaystyle \operatorname{Poly}{\left( 2 s^{3} + 3 s^{2} + s + K, s, domain=\mathbb{Z}\left[K\right] \right)}$
In [6]:
routh(p)
Out[6]:
$\displaystyle \left[\begin{matrix}2 & 1\\3 & K\\1 - \frac{2 K}{3} & 0\\K & 0\end{matrix}\right]$

2.Root Locus¶

In [7]:
K = 1
G = K/(s*(s+1)*(2*s+1))
CLTF = G /(1+G)
CLTF = sympy.cancel(CLTF)
CLTF = sympy.expand(CLTF)
OLTF = K*G
In [8]:
num, den = [[float(i) for i in sympy.Poly(i, s).all_coeffs()] for i in CLTF.as_numer_denom()]
CLTF_system = ct.TransferFunction(num, den)
num_ol, den_ol = [[float(i) for i in sympy.Poly(i, s).all_coeffs()] for i in OLTF.as_numer_denom()]
OLTF_system = ct.TransferFunction(num_ol, den_ol)
In [9]:
CLTF_system
Out[9]:
$$\frac{1}{2 s^3 + 3 s^2 + s + 1}$$
In [10]:
ct.root_locus(CLTF_system, plot = True)
# Customize the plot (optional)
plt.title('Root Locus Plot')
plt.xlabel('Real Axis')
plt.ylabel('Imaginary Axis')
plt.grid(True)
plt.show()

root locus를 그릴 때 CLTF를 사용하면 안됨
open loop transfer function을 사용

In [11]:
OLTF_system
Out[11]:
$$\frac{1}{2 s^3 + 3 s^2 + s}$$
In [12]:
ct.root_locus(OLTF_system, plot = True)
# Customize the plot (optional)
plt.title('Root Locus Plot')
plt.xlabel('Real Axis')
plt.ylabel('Imaginary Axis')
plt.grid(True)
plt.show()
In [13]:
# Find the gain K where the root locus crosses the imaginary axis
# Look for poles with zero real part
# Generate root locus data
poles, gain_values = ct.root_locus(OLTF_system, plot=True)

# Initialize variables
crossing_gain = None
crossing_poles = None
check = np.isclose(poles.real, 0, atol=0.02)[:,1:3]
index = np.where(np.all(check, axis=1))[0]



crossing_gain = gain_values[index]
crossing_poles = poles[index]


# Check if a crossing point was found and print the results
if crossing_gain is not None and crossing_poles is not None:
    print("Gain K at which the root locus crosses the imaginary axis:", crossing_gain)
    print("Poles at this gain:\n", crossing_poles)
else:
    print("No crossing of the imaginary axis was found within the calculated range.")
real_num = np.real(poles[index][0])
img_num = np.imag(poles[index][0])

plt.plot(real_num, img_num, marker='o', markersize=10, color='red', linestyle='')  # 'o' is a circle marker
plt.show()

Gain K at which the root locus crosses the imaginary axis: [1.34879279 1.68599098]
Poles at this gain:
 [[-1.47163858+0.j         -0.01418071+0.67680215j -0.01418071-0.67680215j]
 [-1.53264156+0.j          0.01632078+0.741459j    0.01632078-0.741459j  ]]

3.Nyquist plot¶

In [14]:
# Create a Nyquist plot
ct.nyquist_plot(OLTF_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.xlim([-1.2, 0])      # X축의 범위: [xmin, xmax]
plt.ylim([-0.1, 0.1]) 
plt.show()
In [15]:
# Calculate margins and associated crossover frequencies
gm, pm, wg, wp = ct.margin(OLTF_system)

# Display the gain margin (GM) and phase margin (PM)
print("Gain Margin (in dB):", 20*np.log10(gm))
print("Gain Margin:", gm)

print("Phase Margin (in degrees):", pm)

Gain Margin (in dB): 3.521825181113623
Gain Margin: 1.4999999999999998
Phase Margin (in degrees): 11.424981844921405