【入門】状態空間モデルをPID制御(Python)【数値計算】

MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/

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

はじめに

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

を書き直したもの。

DCモータの状態空間モデルは作成したが、
そこにPID制御器を追加。
これをPython(Numpy)で実現する話。

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

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を使用。
classではあるが、構造体的な使い方しかしてない。
Pythonの場合は、事前にclassの構造を明確に定義してあげる必要がある。
(一般的なプログラムの仕様)