[Python] 限界感度法で調整されたPIDコントローラの目標値追従性能が低いことを示すスクリプト

目的

ジーグラ・ニコルス(Ziegler-Nichols)の限界感度法は有名である。一方、限界感度法で調整されたPIDコントローラの制御結果はオーバーシュートが大きく、収束までに時間がかかるため目標値追従性能は低い。『Pythonによる制御工学入門』を参考にして限界感度法で調整されたPIDコントローラの目標値追従性能が低いことを示すスクリプトを作成したので、制御結果を確認していただき、モデルマッチング等の他の手法を試すきっかけになれば幸いである。

尚、限界感度法は電子計算機(ENIAC)誕生よりも早い1943年に考案されており時代背景から、目標値が変化する自動制御ではなく目標値が常に一定であるレギュレータ用として考案されている可能性がある。

シミュレーション条件

・プラント：２次遅れむだ時間系
・ゲイン調整方法：PID制御(※1) と No Overshoot (※2)

　(※1)はZN(Original)、(※2)はZN(No overshort)と記載する。

P制御で持続発振するPゲインと発振周期の確認

プラントのパラメータやゲインは『Pythonによる制御工学入門』に記載のものを使用したが念のため持続発振するPゲインと発振周期に誤りがないことを確認した。
持続発振するPゲイン(Kp0) = 2.9, 発振周期(T0) = 0.3 [s]

制御結果 ( Step応答 )

以下に示すように、ZN(Original)はオーバーシュートも大きく収束までに時間がかかる。また、ZN(No overshort)はオーバーシュートがないものの収束までに時間がかかる。

参考情報

Pythonによる制御工学入門
限界感度法の利点・欠点・使い所。実はアテにしちゃダメ！？

スクリプト

a05_Ziegler_Nichols.ipynb

d12_Pytorch_CNN_Time_Sequence.ipynb

from control.matlab import*
import numpy as np
import matplotlib.pyplot as plt
import time
import torch
import torch.nn as nn
import torch.optim as optim
from torchinfo import summary
from torch.utils.data import DataLoader, TensorDataset
torch.manual_seed(0);

#--------プラントの作成----------------
#https://nbviewer.org/github/373yuki/pyctrl2/blob/main/chap5.ipynbより
g  = 9.81                # 重力加速度[m/s^2]
l  = 0.2                 # アームの長さ[m]
M  = 0.5                 # アームの質量[kg]
mu = 1.5e-2              # 粘性摩擦係数[kg*m^2/s]
J  = 1.0e-2              # 慣性モーメント[kg*m^2]

P_tmp = tf( [0,1], [J, mu, M*g*l] ) #2次遅れ系
num_delay, den_delay = pade( 0.005, 1) # 1次のパデ近似
P = P_tmp * tf(num_delay, den_delay) #２次遅れむだ時間系

#---------時間刻み-------------------
# 時間軸の作成
duration = 10  # 信号の継続時間 (秒)
sampling_rate = 1000  # サンプリングレート (Hz)
time = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)

#-------------目標値信号の生成---------
def make_signal():
    step_time = duration/3 #ステップするタイミングをdurationの1/3毎とする
    u = np.zeros_like(time)
    for i in range(len(time)):
        if (i // int(step_time * sampling_rate)) % 2 == 1:
            u[i] = 1 #目標角度 [deg]
    return u

r_step = make_signal(); #目標値信号 

#----------------応答の計算する関数---------------
def calc_y(P, u, t):
    x0 = 0 #初期値
    y, time_tmp, x_tmp = lsim(P, U=u, T=t, X0=x0)
    return y

#--------限界感度法で定常コントローラの作成-------------
Kp0 = 2.9 #持続振動する比例ゲイン

C = Kp0
G = feedback(P*C, 1) #閉ループ伝達関数の作成 ( 1はフィードバックゲイン )#--------コントローラの作成-------------
y_act = calc_y(G, r_step, time) 

#---グラフ表示---------
plt.figure(figsize=(10, 5))

#全体表示
plt.subplot(2, 1, 1)
plt.axhline(0, color="gray", linestyle="--", linewidth=1)
plt.plot(time, r_step, label="Target", color='b', linewidth=1)
plt.plot(time, y_act, label="Actual", color='r', linewidth=1)
plt.ylabel("Signal [-]")
plt.title("Step responce")
plt.legend()
plt.grid(True)

#拡大表示
plt.subplot(2, 1, 2)
plt.axhline(0, color="gray", linestyle="--", linewidth=1)
plt.plot(time, r_step, label="Target", color='b', linewidth=1)
plt.plot(time, y_act, label="Actual", color='r', linewidth=1)
plt.ylabel("Signal [-]")
plt.legend()
plt.grid(True)
plt.xlim(3.5,4.1)

#---------Ziegler_Nicholsの限界感度法で決めたPID-----------
Kp0 = 2.9 #持続振動する比例ゲイン
T0 = 0.3 #持続振動の周期

#ZN オリジナル
Kp = 0.6 * Kp0
Ki = Kp / (0.5 * T0)
Kd = Kp * (0.125 * T0)
C = tf([Kd, Kp, Ki], [1, 0])
G = feedback(P*C, 1) #閉ループ伝達関数の作成 ( 1はフィードバックゲイン )
y_act_org = calc_y(G, r_step, time)

#ZN改良版(オーバーシュートなし)
Kp = 0.2 * Kp0
Ki = Kp / (0.5 * T0)
Kd = Kp * (0.33 * T0)
C = tf([Kd, Kp, Ki], [1, 0])
G = feedback(P*C, 1) #閉ループ伝達関数の作成 ( 1はフィードバックゲイン )
y_act_no_overshort = calc_y(G, r_step, time)

#---グラフ表示---------
plt.figure(figsize=(10, 5))

#全体表示
plt.subplot(2, 1, 1)
plt.axhline(0, color="gray", linestyle="--", linewidth=1)
plt.plot(time, r_step, label="Target", color='b', linewidth=1)
plt.plot(time, y_act_org, label="ZN(Original)", color='r', linewidth=1)
plt.plot(time, y_act_no_overshort, label="ZN(No overshort)", color='g', linewidth=1)
plt.ylabel("Signal [-]")
plt.title("Step responce")
plt.legend()
plt.grid(True)

#拡大表示
plt.subplot(2, 1, 2)
plt.axhline(0, color="gray", linestyle="--", linewidth=1)
plt.plot(time, r_step, label="Target", color='b', linewidth=1)
plt.plot(time, y_act_org, label="ZN(Original)", color='r', linewidth=1)
plt.plot(time, y_act_no_overshort, label="ZN(No overshort)", color='g', linewidth=1)
plt.ylabel("Signal [-]")
plt.legend()
plt.grid(True)
plt.xlim(3,6)