KOCW design lag compensator¶

예제 9.1 풀이

In [1]:
import sympy
from sympy import solve
from scipy.optimize import fsolve
import control as ct
import numpy as np
import matplotlib.pyplot as plt

root locus가 damping ration와 만나는 점 확인¶

In [2]:
zeta = 0.174
theta = np.arcsin(0.174)
radian_ = np.pi/2 - theta
np.rad2deg(radian_)
Out[2]:
79.97953045136717
In [3]:
tan_factor = np.tan(np.pi/2 - theta)
tan_factor
Out[3]:
5.659457772645189

wd/σ = 5.66
wd = 5.66 * σ

In [4]:
import scipy.optimize as opt
import numpy as np

def equation(x):
    return -np.arctan(tan_factor*x/(10+x)) - np.arctan(tan_factor*x/(2+x))\
- np.arctan(tan_factor*x/(1+x)) - np.radians(180)

# Initial guess for x
initial_guess = 0.01

# Solve the equation numerically
solution = opt.fsolve(equation, initial_guess)

print("real number:", solution[0])
print("imaginary number:{}i".format(solution[0]*tan_factor))

real number: -0.6936132551253392
imaginary number:-3.925474927928832i
In [5]:
from sympy import I 
s0 = solution[0]+I*solution[0]*tan_factor
print("s0 is {}".format(s0))
K_notcompansated = abs(s0+10)*abs(s0+2)*abs(s0+1)
K_notcompansated

s0 is -0.693613255125339 - 3.92547492792883*I
Out[5]:
$\displaystyle 164.532228522218$
In [6]:
from sympy import I, atan2, im, re, Abs
# Function to calculate the angle of a complex number in radians
def complex_angle(expr):
    return atan2(im(expr), re(expr))

# Calculate the angles
angle1 = complex_angle(s0 + 10)
angle2 = complex_angle(s0 + 2)
angle3 = complex_angle(s0 + 1)

# Sum of angles
total_angle = angle1 + angle2 + angle3

print("Total angle in radians:", total_angle)

# If you want to convert the angle to degrees
total_angle_degrees = total_angle * 180 / sympy.pi
print("Total angle in degrees:", total_angle_degrees)

Total angle in radians: -3.14159265358979
Total angle in degrees: -565.486677646163/pi
In [7]:
K = K_notcompansated
s = sympy.Symbol('s')
G = 1/((s+10)*(s+1)*(s+2))
CLTF = K * G /(1+K*G)
CLTF = sympy.cancel(CLTF)
OLTF = K*G
In [8]:
#유리식을 약분
#cancle이 없으면 약분하지 않음.
#3차식이 6차식으로 바뀜.
#sympy.cancel(CLTF)
In [9]:
num, den = [[float(i) for i in sympy.Poly(i, s).all_coeffs()] for i in CLTF.as_numer_denom()]
num_ol, den_ol = [[float(i) for i in sympy.Poly(i, s).all_coeffs()] for i in OLTF.as_numer_denom()]
CLTF_system = ct.TransferFunction(num, den)
OLTF_system = ct.TransferFunction(num_ol, den_ol)
In [10]:
CLTF_system
Out[10]:
$$\frac{164.5 s^3 + 2139 s^2 + 5265 s + 3291}{s^6 + 26 s^5 + 233 s^4 + 1037 s^3 + 3683 s^2 + 6545 s + 3691}$$
In [11]:
OLTF_system
Out[11]:
$$\frac{164.5}{s^3 + 13 s^2 + 32 s + 20}$$
In [12]:
wn, zetas, poles = ct.damp(CLTF_system)

    Eigenvalue (pole)       Damping     Frequency
               -11.61             1         11.61
                  -10             1            10
   -0.6936    +3.925j         0.174         3.986
   -0.6936    -3.925j         0.174         3.986
                   -2             1             2
                   -1             1             1
In [13]:
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()
In [14]:
# Time vector
time, response = ct.step_response(CLTF_system)
# Plot the step response
plt.plot(time, response)
plt.title('Step Response')
plt.xlabel('Time (seconds)')
plt.ylabel('Amplitude')
plt.grid(True)
# Set x-axis limits
plt.xlim(0, 10)  # x-axis will range from 0 to 10
plt.show()
In [15]:
ct.step_info(CLTF_system)
Out[15]:
{'RiseTime': 0.32087752924572477,
 'SettlingTime': 5.1459248208666235,
 'SettlingMin': 0.6149955074584611,
 'SettlingMax': 1.3734915686374334,
 'Overshoot': 54.04487150851232,
 'Undershoot': 0,
 'Peak': 1.3734915686374334,
 'PeakTime': 0.8913264701270133,
 'SteadyStateValue': 0.8916178482199818}

ideal compensator를 추가했을 경우 새로운 Gain¶

compensator = K(s+0.1)/s
s0는 compensator를 추가하기 전에 구했음.

uncompansated system에서 s0를 고정하여, compansated에 삽입¶

In [16]:
s0
Out[16]:
$\displaystyle -0.693613255125339 - 3.92547492792883 i$
In [17]:
## compansated system의 gain
K_compansated = abs(s0+10)*abs(s0+2)*abs(s0+1)*abs(s0+0.1) /abs(s0)
print("compansated system의 새로운 gain은 {}".format(K_compansated))

compansated system의 새로운 gain은 163.864466237033
In [18]:
from sympy import I, atan2, im, re, Abs
import sympy
# Function to calculate the angle of a complex number in radians
def complex_angle(expr):
    return atan2(im(expr), re(expr))

# Calculate the angles
angle1 = complex_angle(s0 + 10)
angle2 = complex_angle(s0 + 2)
angle3 = complex_angle(s0 + 1)
angle4 = complex_angle(s0)
angle5 = complex_angle(s0+0.1)
# Sum of angles
total_angle = angle1 + angle2 + angle3 + angle4 -angle5

print("Total angle in radians:", total_angle)

# If you want to convert the angle to degrees
total_angle_degrees = total_angle * 180 / sympy.pi
print("Total angle in degrees:", total_angle_degrees)

Total angle in radians: -3.16639922089135
Total angle in degrees: -569.951859760442/pi

ideal compensator 추가¶

In [19]:
ideal_compensator = (s+0.1)/(s)
ideal_compensator
Out[19]:
$\displaystyle \frac{s + 0.1}{s}$
In [20]:
K2 = K_compansated
CLTF_compensated = K2 * G * ideal_compensator /(1 + K2 * G * ideal_compensator)
CLTF_compensated = sympy.cancel(CLTF_compensated)
num_com, den_com = [[float(i) for i in sympy.Poly(i, s).all_coeffs()] for i in CLTF_compensated.as_numer_denom()]
CLTF_compensated_system = ct.TransferFunction(num_com, den_com)
CLTF_compensated
Out[20]:
$\displaystyle \frac{163.864466237033 s^{4} + 2146.62450770514 s^{3} + 5456.68672569321 s^{2} + 3801.65561669917 s + 327.728932474067}{1.0 s^{7} + 26.0 s^{6} + 233.0 s^{5} + 1035.86446623703 s^{4} + 3690.62450770514 s^{3} + 6736.68672569321 s^{2} + 4201.65561669917 s + 327.728932474067}$
In [21]:
wn_com, zeta_com, poles_com = ct.damp(CLTF_compensated_system)

    Eigenvalue (pole)       Damping     Frequency
                -11.6             1          11.6
                  -10             1            10
   -0.6561    +3.897j         0.166         3.951
   -0.6561    -3.897j         0.166         3.951
                   -2             1             2
                   -1             1             1
              -0.0905             1        0.0905

damping_ratio = 0.166

In [24]:
#compansated system
# Time vector
time, response = ct.step_response(CLTF_compensated_system)
# Plot the step response
plt.plot(time, response)
plt.title('Step Response')
plt.xlabel('Time (seconds)')
plt.ylabel('Amplitude')
plt.grid(True)
# Set x-axis limits
plt.xlim(0, 10)  # x-axis will range from 0 to 10

#uncompansated system
# Time vector
time, response = ct.step_response(CLTF_system)
# Plot the step response
plt.plot(time, response)
plt.title('Step Response')
plt.xlabel('Time (seconds)')
plt.ylabel('Amplitude')
plt.grid(True)
# Set x-axis limits
plt.xlim(0, 10)  # x-axis will range from 0 to 10

plt.show()
In [23]:
ct.step_info(CLTF_compensated_system)
Out[23]:
{'RiseTime': 0.3359298355175059,
 'SettlingTime': 17.39200375701996,
 'SettlingMin': 0.6214103019543759,
 'SettlingMax': 1.4131621197095117,
 'Overshoot': 41.316211970951166,
 'Undershoot': 0,
 'Peak': 1.4131621197095117,
 'PeakTime': 0.9009027407060384,
 'SteadyStateValue': 1.0}