주파수 영역에서의 PID 제어기 설계 (SISO)

대충 써서 한번 정리해야할 것 같다.

 

들어가며

제어 공학에서 어떤 시스템 모델이 선형 시불변 시스템(LTI, Linear Time-Invariant system)이라면 푸리에 변환을 통해 주파수 영역에서의 해석 기법들을 사용할 수 있다.

전통적으로 단일 입력, 단일 출력(SISO, Single-Input Single-Output) 인 전달함수에 대해서 모델 특성을 보면 주로 이득 여유와 위상 여유를 본다.

Gain/Phase margin on Bode plot

사전적인 의미로

이득 여유는 위상 곡선이 -180도를 지나는 지점에서의 크기가 0 dB으로부터의 간격(여유)

위상 여유는 이득 곡선이 0 dB를 지나는 지점에서의 위상 -180도까지의 간격이다.

체감 의미를 대충 요약하면 모델 오차가 저 정도의 여유가 나도 괜찮다는 그런 의미이다.

 

그렇다면 이를 이용해서 여유를 가지는 제어기를 설계할 수 있다.

 

PID 제어기의 특성

PID 제어기의 각 특성을 요약해보면 다음과 같다.

종류	특성
P	시스템 전체 이득을 조절한다.
I	크기는 우하강이며, 위상을 -90도 떨어뜨린다.
D	크기는 우상승이며, 위상을 90도 올린다.
PI	시스템 전체 이득을 조절하면서 (P) 이득 비율 지점을 기준으로  저주파 영역에서 크기를 좌상승 (I) 저주파 영역에서  위상을 -90도 떨어뜨린다. (I) $$G_c (s) = K_p + K_i \frac{1}{s} = K_p \left( 1+ \frac{1}{\tau_I s} \right)$$
PD	시스템 전체 이득을 조절하면서 (P) 이득 비율 지점을 기준으로 고주파 영역에서 크기를 우상승 (D) 고주파 영역에서 위상을 90도 올린다. (D) $$G_c (s) = K_p + K_d s = K_p \left( 1+ \tau_d s \right)$$
PID	PI, PD의 특성을 다 가진다.

 

문제 상황

내용 무

 

제어기 설계

내용 무

 

코드 구현

메인 파일 : control_test.py

import matplotlib.pyplot as plt
# from control.matlab import *
import control as ctrl
import numpy as np
from control_lib import *

# 1st-order model
mag = 2
tau = 0.1
# model = ctrl.tf([mag], [tau, 1])

# 2nd-order model
# mag = 2
# wn = 10
# zeta = 1
# model = ctrl.tf([mag*wn*wn], [1, 2*zeta*wn, wn*wn])

# 3rd-order model
mag = 2
tau = 1
wn = 10
zeta = 1
model = ctrl.tf([mag*wn*wn], [1, 2*zeta*wn, wn*wn]) * ctrl.tf([mag], [tau, 1])

mag, phase, omega = ctrl.bode_plot(model, plot=False, margins=True)
gm, pm, wcg, wcp = ctrl.margin(model)
print("Gain cross  : {:6.2f} dB  at {:6.2f} rad/s".format(20*np.log10(gm), wcg))
print("Phase cross : {:6.2f} deg at {:6.2f} rad/s".format(pm, wcp))

wco_target=wcg/6
Gc, kp, ki = design_pid(model=model, pid_type='pi', wco_target=wco_target, pm_target=60)
mag, phase, omega = ctrl.bode_plot(model*Gc, plot=False, margins=True)

show_bode_plot(model, fignum=1, margin=True)
show_bode_plot(model*Gc, fignum=2, margin=True)

#####################################################################
# Step response
end_time = 5
t1, y1 = ctrl.step_response(model, end_time)
t2, y2 = ctrl.step_response(ctrl.feedback(model*Gc), end_time)
fig = plt.figure(10)
ax = fig.add_subplot(111)
ax.plot(t1, y1)
ax.plot(t2, y2)
ax.set_xlabel('Time [sec]'); ax.set_ylabel('Amplitude'); ax.grid()
ax.legend(['Model','Feedback'])

plt.tight_layout()
plt.show()

 

라이브러리 파일 : control_lib.py

import matplotlib.pyplot as plt
import control as ctrl
import numpy as np

def show_bode_plot(model, fignum = 1, margin:bool = False):
    mag, phase, omega = ctrl.bode_plot(model, plot=False)

    fig = plt.figure(fignum)
    fig.clf()
    ax1 = fig.add_subplot(211)
    ax1.semilogx(omega, 20*np.log10(mag))
    ax1.set_xlim([omega[0], omega[-1]])
    ax1.set_ylabel('Magnitude [dB]'); ax1.grid()

    ax2 = fig.add_subplot(212)
    ax2.semilogx(omega, np.rad2deg(phase))
    ax2.set_xlim([omega[0], omega[-1]])
    ax2.set_xlabel('Omega [rad/s]'); ax2.set_ylabel('Phase [deg]'); ax2.grid()

    ax1.semilogx([omega[0], omega[-1]], [0, 0],'k:') # Zero-Crossing Line

    if(margin == True):
        gm, pm, wcg, wcp = ctrl.margin(model)
        if(not np.isinf(pm) or not np.isinf(gm)):
            phase_limit = np.ceil(np.min(phase)/(2*np.pi))*360 - 180
            ax2.semilogx([omega[0], omega[-1]],[phase_limit, phase_limit], 'k:') # Imaginary line
        if(not np.isinf(gm)):
            ax1.semilogx([wcg], [-20*np.log10(gm)], 'or', markersize=2)
            ax1.semilogx([wcg, wcg], [0, -20*np.log10(gm)], 'k:') # perpendicular line @ phase margin
            ax2.semilogx([wcg], [phase_limit], 'or', markersize=2) # perpendicular line @ phase margin
            ax1.set_title("Gm {:.2f} dB @ {:.1f} rad/s".format(20*np.log10(gm), wcg))
        if(not np.isinf(pm)):
            ax1.semilogx([wcp], [0], 'or', markersize=2) # Gain cross-over point
            ax2.semilogx([wcp], [phase_limit+pm], 'or', markersize=2) # Phase margin point
            ax2.semilogx([wcp, wcp], [phase_limit, phase_limit+pm], 'k:') # perpendicular line @ phase margin
            ax2.set_title("Pm {:.2f} deg @ {:.1f} rad/s".format(pm, wcp))
    plt.tight_layout()


def design_pid(model, pid_type, wco_target, pm_target=None):
    '''
        model : Linear system model
        pid_type ::
            p
            pi
            pd
            pid
        wco_target : Gain cross-over frequency
        pm_target : Phase margin for pi, pd
    '''
    gm, pm, wcg, wcp = ctrl.margin(model) # _, deg, rad/s, rad/s
    mag, phase, omega = ctrl.freqresp(model, wco_target)
    _, phase_end, _ = ctrl.freqresp(model, wco_target*100)
    print("phase {:.2f} deg @ wco {:.2f} rad/s".format(np.rad2deg(phase[0]), wco_target))

    if(pid_type.lower() == 'p'):
        kp  = 1 / mag[0]
        Gc  = kp
        print("[Designer] P    = kp {:.3f}".format(kp))
        print(Gc)
        return Gc, kp
    
    elif(pid_type.lower() == 'pi'):
        phase_limit = np.ceil(np.min(phase_end)/(2*np.pi))*360 - 180
        taui= 1 / (wco_target * np.tan(np.deg2rad((np.rad2deg(phase[0]) + phase_limit) - pm_target)))
        kp  = 1 / (mag[0] * np.linalg.norm([1, 1/(taui * wco_target)]))
        ki  = kp / taui
        Gc  = kp * ctrl.tf([taui, 1],[taui, 0])
        print("[Designer] PI   = kp {:.3f} ki {:.3f} taui {:.3f}".format(kp, ki, taui))
        print(Gc)
        if(taui < 0):
            _, now_pm, _ = ctrl.freqresp(model/mag[0], wco_target)
            assert False, "[Designer] Failed to design PI : PM should be under {:.2f} deg".format(np.rad2deg(now_pm[0])+phase_limit)
        return Gc, kp, ki

    elif(pid_type.lower() == 'pd'):
        phase_limit = np.ceil(np.min(phase_end)/(2*np.pi))*360 - 180
        taud= np.tan(np.deg2rad((-np.rad2deg(phase[0]) + phase_limit) + pm_target))/wco_target
        kp  = 1 / (mag[0] * np.linalg.norm([1, (taud * wco_target)]))
        kd  = kp / taud
        Gc  = kp * ctrl.tf([taud, 1],[1])
        print("[Designer] PD   = kp {:.3f} kd {:.3f} taud {:.3f}".format(kp, kd, taud))
        print(Gc)
        if(taud < 0):
            _, now_pm, _ = ctrl.freqresp(model/mag[0], wco_target)
            assert False, "[Designer] Failed to design PD : PM should be under {:.2f} deg".format(np.rad2deg(now_pm[0])+phase_limit)
        return Gc, kp, ki

 

 

 

Reference

 

** EOF **