MATLAB,Python,Scilab,Julia比較 その78【PID制御⑦】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia-backnumber/

はじめに
登場人物
【再掲】PID制御器のブロック線図と全体構成
Pythonコード
シミュレーション結果
まとめ

はじめに

前回は、MATLABによるDCモータ状態空間モデルをPID制御のシミュレーションを実施。
問題無く動作確認はできた。

今回は、これのPython(Numpy)版

登場人物

博識フクロウのフクさん

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

【再掲】PID制御器のブロック線図と全体構成

太郎くん

今回からPython(Numpy)版のシミュレーションだね。

フクさん

前回と同じく、PID制御器のブロック線図と全体構成を再掲しておこう。

PID制御器(変形式)ブロック線図、目標値、実値、1/zは前回値保持の意、Ki、Kp、Kd、T1、ΔT、e(t)、u(t)

PID制御器を加味した構成図、指令器0~1の範囲で出力、指令器、制御器、PID制御器、DCモータ、状態空間モデル、target(t)、actual(t)、E(t)、u(t)、ω(t)

太郎くん

コードだけだと接続で迷子になっちゃうもんね。

Pythonコード

フクさん

そしてPythonコードは以下となる。

import numpy as np
import matplotlib.pyplot as plt


def statespacemodel(A,B,C,D,u,dt,x):
  # 状態方程式
  x = x + (A@x + B@u) * dt
  
  # 出力方程式
  y = C@x + D@u
  return x,y

class PID_state:
  def __init__(self):
    self.pzi = 0  # P項前回値
    self.izi = 0  # I項前回値
    self.dzi = 0  # D項前回値


def PIDController(state, target, actual, Kp, Ki, Kd, dt, t1):
  e_i = target - actual
  e_p = (t1/dt)*(e_i - state.pzi)
  e_d = (t1/dt)*(e_p - state.dzi)
  u = (dt/t1)*(e_p*Kp + e_i*Ki + e_d*Kd) + state.izi
  
  state.pzi = e_i
  state.dzi = e_p
  state.izi = u
  
  return state, u


def statespacemodel_pid():
  K=0.016
  J=0.000000919
  R=1.34
  L=0.00012
  
  A=np.array([[0,1,0],[0,0,K/J],[0,-K/L,-R/L]])
  B=np.array([[0],[0],[1/L]])
  C=np.array([[1,0,0],[0,1,0],[0,0,1]])
  D=np.array([[0],[0],[0]])
  
  dt = 0.0001
  t = np.linspace(0, 1, 10000) # 時間(横)軸
  u = np.zeros((1,10000))	# 入力信号生成
  u[0][5000:10000]=1	# 0.5秒後に0から1へ
  y = np.zeros((3,len(t)))
  x = np.zeros((3,1))
  
  state = PID_state()
  ratio = 1/60
  Kp = 0.80
  Ki = 0.45
  Kd = 0.0
  t1 = 0.005
  
  omega = 0
  uPID = np.zeros((1,10000))
  
  for i in range(0, len(t)):
    state,uPID[0,[i]] = PIDController( state, u[:,[i]], omega*ratio, Kp, Ki, Kd, dt, t1 )
    x,y[:,[i]] = statespacemodel(A,B,C,D,uPID[:,[i]],dt,x)
    omega = y[[1],[i]]
  
  fig = plt.figure()
  ax1 = fig.add_subplot(3, 1, 1)
  ax2 = fig.add_subplot(3, 1, 2)
  ax3 = fig.add_subplot(3, 1, 3)
  
  ax1.plot(t,uPID.T,'-r', t,u.T, '--b')
  ax1.set_xlim(0,1)
  ax2.plot(t,y[0:2,:].T)
  ax2.set_xlim(0,1)
  #ax2.set_ylim(-1,80)
  ax3.plot(t,y[2,:].T)
  ax3.set_xlim(0,1)
  
  plt.show()

if __name__ == "__main__":
  statespacemodel_pid()

太郎くん

MATLABの時は構造体使ってたけど、
Pythonの場合はclassを使うのね。