The Negligible Lab

Astray in the forest of electrical and electronic circuits. Adrift in the gap between time and frequency domains. 独立独歩を是とする。

続・LTspiceによる電流制御のシミュレーション

はじめに

先日の記事で,LTspiceによる電流制御のシミュレーションについて述べました。 negligible.hatenablog.com

この記事の中で私個人としてひとつの発見と言えるのは,

  • 操作量に対してむだ時間T,フィードバック量に対して窓Tの移動平均フィルタを施した連続時間制御システム
  • 操作量に対してむだ時間TとPWM,フィードバック量に対してサンプリング周期Tのサンプル & ホールドを設けた離散時間制御システム

がほぼ同じ挙動を示すということです。 図1に上記前者のブロック図を示します。

f:id:s-inoue2010:20200714180250p:plain:w500
図1: 操作量にむだ時間,フィードバック量に移動平均フィルタを施した連続時間システム

言い換えれば,離散時間制御 + PWMの電流制御系を,移動平均フィルタとむだ時間のある連続制御系で近似できると言えるのではないでしょうか(PWMの過変調が起きない範囲で)。 例えば図2ですね。

f:id:s-inoue2010:20200714180654p:plain
図2: 連続時間制御と離散時間制御 + PWM

ほぼ挙動は一致しています。 そこで,連続系である図1の安定性について検討してみます。 また,最近はLTspiceばかり触っているので,JupyterLabでボード線図を描いてみました。

f:id:s-inoue2010:20200702023510p:plain:w260
図3: RL回路

回路定数は先日の記事と同じくR = 20 mΩ,L = 5 mHとします。 また,同じくサンプリング時間 = 移動平均フィルタの窓 = むだ時間T = 100 μsとします。

何が問題か

先日の記事では,むだ時間も移動平均もないブロック図(図4)について, 開ループ伝達関数(一巡伝達関数)Go(s)が次式で表されること,

G_{o}(s) = \displaystyle \left (K_{P} + \frac{K_{I}}{s} \right ) \frac{1}{s L + R} = \frac{K_{P}}{s L} \frac{s + \frac{K_{I}}{K_{P}}}{s + \frac{R}{L}} \tag{1}

および,

\displaystyle \frac{K_{I}}{K_{P}} = \frac{R}{L} \tag{2}

を満たすように比例ゲインKP積分ゲインKIを設定できれば,極零相殺されて

G_{o}(s) =\displaystyle \frac{K_{P}}{s L} \tag{3}

となることを述べました。

f:id:s-inoue2010:20200702030157p:plain:w520
図4: むだ時間も移動平均もない理想的な連続時間制御システム

(3)式は当然ですが,角周波数ωに関係なく,その位相は−90°となります。 (3)式のゲインが1 (= 0 dB)になるのはω = KP / Lの場合であり,位相余裕は90°となります。 したがって安定です。比例ゲインKPをいくつに設定しても安定になります。 しかも極零相殺のおかげで綺麗な1次遅れ系です。

しかし,離散時間(ディジタル)制御システムでサンプル & ホールドしたり,制御演算に1サンプリング時間を要するような場合は, 比例ゲインKPを高く設定し過ぎると振動したり,発散したりするおそれがあります。 先日の記事のように,比例ゲインKP ≤ L / (4 T)に設定しておけば,ステップ応答はオーバーシュートしないことが知られています。 KP > L / (4 T)となった場合はステップ応答にオーバーシュートが生じますが,直ちに不安定にはなりません。 それでは,KPはどこまで増やせるのでしょうか?

ボード線図を描く

という訳で,開ループ伝達関数のボード線図を描いてみました。 自作制御ライブラリContrailleを使えばLTspiceでもボード線図を簡単に描けるのですが,今回はJupyterLabを使ってみました。 図4と図1の開ループ伝達関数Gi(s),Gs(s)は次式となります。なお,(4)式は(1)式の再掲です。

G_{i}(s) =  \displaystyle \frac{K_{P}}{s L} \frac{s + \frac{K_{I}}{K_{P}}}{s + \frac{R}{L}} \tag{4}
G_{s}(s) = \displaystyle \frac{K_{P}}{s L} \frac{s + \frac{K_{I}}{K_{P}}}{s + \frac{R}{L}} \frac{e^{-s T} - e^{- 2 s  T}}{s T} \tag{5}

ここで,

K_{P} = K \displaystyle \frac{L}{4 T}, \quad K_{I} = K_{P} \frac{R}{L} \tag{6}

として,パラメータKを変更しながらGi(s),Gs(s)のボード線図を描いてみましょう。 なお,図5 ~ 7において"Ideal"はGi(s),"Sampled"はGs(s)を示します。

f:id:s-inoue2010:20200714183658p:plain
図5: パラメータK = 1の場合のボード線図

f:id:s-inoue2010:20200714183805p:plain
図6: パラメータK = 3の場合のボード線図

f:id:s-inoue2010:20200714183908p:plain
図7: パラメータK = 4.5の場合のボード線図

図5 ~ 7は例としてそれぞれパラメータK = 1,K = 3,K = 4.5の場合を描いています。 それぞれGs(s)の位相余裕PMが68.57°,26.98°,−2.12°となりました。 そう! K = 4.5の場合には不安定になります! 安定不安定の境界はK = 4.2付近にあるようです。

Pythonソースコードは後述します。

LTspiceでのシミュレーション

さて,それではK = 1, 3, 4.5の場合をシミュレーションしてみます。 モデル(schematic)は同じく図7となります。

  • 連続時間制御 + PWMなし
  • 離散時間制御 + PWMあり

の2つを比較します。 PWMする場合,直流電圧250 Vの単相フルブリッジインバータを想定します。

f:id:s-inoue2010:20200714201436p:plain
図7: 電流制御のschematic

f:id:s-inoue2010:20200714201815p:plain
図8: K = 1の場合のシミュレーション波形

K = 1の場合はまさにKP = L / (4 T)の場合であり,理論上はステップ応答にオーバーシュートが発生しない波形となります。 位相余裕は60°以上あり,申し分ない(?)安定した波形となります。

f:id:s-inoue2010:20200715001653p:plain
図9: K = 3の場合のシミュレーション波形

K = 3の場合はオーバーシュートが発生しますが,位相余裕が約27°あり,まだ安定です。 ただし,文献*1の6.3節の最後,p. 223,表6.1にて,きなみ・みゆう姉妹の会話の中で示されている通り,サーボ系としては位相余裕40° ~ 60°は欲しいところです。

f:id:s-inoue2010:20200715002633p:plain
図10: K = 4.5の場合のシミュレーション波形

K = 4.5となると,位相余裕が負となってしまい,不安定になります。 シミュレーションからは,連続系が発散して行く様子が分かりますね。 ただし,PWMでは操作量 = 電圧振幅に制限があるため,発散というよりは持続振動になります。

以上より,連続時間での電流制御系として近似した開ループ伝達関数のボード線図から得られる位相余裕にで,離散時間での電流制御系の安定・不安定を判別できることが分かりました。

(注)位相余裕が0°に近い,安定か不安定かの境界付近では, 連続時間系での近似は安定でも,離散時間での電流制御系が不安定になる場合もあるようです。

ボード線図を描くPythonソースコード

最後に,位相余裕を求め,ボード線図を描いたPythonソースコードを示します。

# -*- coding: utf-8 -*-
"""
current-control-bode.py
Copyright (c) 2020 @RR_Inyo
Released under the MIT license.
https://opensource.org/licenses/mit-license.php
"""

# 電流制御系の安定性を確認する
# 準備
import numpy as np
from scipy.optimize import newton
import matplotlib.pyplot as plt

# パラメータの定義
## 主回路定数
R = 20e-3 # 抵抗器, [Ω]
L = 5e-3 # インダクタ, [H]

## 制御定数
T = 100e-6 # サンプリング周期, [s]
K = 3.8 # 次行のKP = L / (4 * T)との比, K = 1でオーバーシュートなし
KP = K * L / (4 * T) # 比例ゲイン, [V/A]
KI = KP * (R / L) # 積分ゲイン, [V/As]

print(f'比例ゲインKP = {KP} [V/A]')
print(f'積分ゲインKI = {KI} [V/As]')

# 開ループ周波数伝達関数の定義
## 理想
def Gi(s):
    G1 = KP + KI / s # PI制御
    G3 = 1 / (R + s * L)
    return G1 * G3

## サンプリング考慮
def Gs(s):
    G1 = KP + KI / s # PI制御
    G2 = np.exp(-s * T) # 1サンプル遅れ
    G3 = 1 / (R + s * L)
    G4 = (1 - np.exp(-s * T)) / (s * T) # サンプル & ホールド
    return G1 * G2 * G3 * G4

# 周波数応答計算
omega = np.logspace(2, 5, 1000)
Gi_FR = Gi(1j * omega)
Gs_FR = Gs(1j * omega)

# 理想的な開ループ伝達関数Gi(s)のゲインが1 = 0 dBとなる角周波数
omega_0_i = KP / L
print(f'Gi(jω)のゲインが1 = 0 dBとなる角周波数ω0 = {omega_0_i} [rad/s]')
print(f'検算|Gi(jω0)| = {20 * np.log10(np.abs(Gi(1j * omega_0_i)))} [dB]')

# サンプリングを考慮したGs(s)のゲインが1 = 0 dBとなる角周波数

omega_0_s = newton(lambda omega: np.abs(Gs(1j * omega)) - 1, KP / L)
print(f'Gs(jω)のゲインが1 = 0 dBとなる角周波数ω0 = {omega_0_s} [rad/s]')
print(f'検算|Gs(jω0)| = {20 * np.log10(np.abs(Gs(1j * omega_0_s)))} [dB]')

# Gs(s)のゲインが1 = 0 dBとなった場合の位相
phase_0_s = np.angle(Gs(1j * omega_0_s)) / np.pi * 180 # [deg]
if phase_0_s > 0:
    phase_0_s -= 360
print(f'そのときの位相 = {phase_0_s} [rad/s]')

# 位相余裕PM
line_omega_0 = [omega_0_s, omega_0_s]
line_PM = [-180, phase_0_s]
PM = phase_0_s + 180
print(f'位相余裕PM = {PM} [deg]')
if PM > 0:
    color_PM = 'deepskyblue'
    pos = 'right'
else:
    color_PM = 'tomato'
    pos = 'left'

# ボード線図
fig, ax = plt.subplots(2, figsize = (12, 6))
fig.patch.set_facecolor('lavender')

## ゲイン線図
ax[0].set_title(f'Comparison of the delay-less and sampled controllers, K = {K}')
ax[0].semilogx(omega, 20 * np.log10(np.abs(Gi_FR)), label = 'Ideal')
ax[0].semilogx(omega, 20 * np.log10(np.abs(Gs_FR)), label = 'Sampled')
ax[0].legend()
ax[0].set_xlim(1e2, 1e5)
ax[0].set_ylabel('Gain [dB]')
ax[0].set_ylim(-80, 40)
ax[0].set_ylabel('Gain')
ax[0].grid()

## 位相線図
ax[1].semilogx(omega, 180 * np.unwrap(np.angle(Gi_FR)) / np.pi, label = 'Ideal')
ax[1].semilogx(omega, 180 * np.unwrap(np.angle(Gs_FR)) / np.pi, label = 'Sampled')
ax[1].semilogx(line_omega_0, line_PM, color = color_PM, alpha = 0.8)
ax[1].text(omega_0_s, (phase_0_s - 180) / 2, f'PM = {PM:.2f} [rad] ', ha = pos, alpha = 0.8)
ax[1].legend()
ax[1].set_xlim(1e2, 1e5)
ax[1].set_xlabel('Angular frequency [rad/s]')
ax[1].set_ylim(-270, 30)
ax[1].set_yticks(range(-270, 60, 30))
ax[1].set_ylabel('Phase [deg]')
ax[1].grid()

ラムダ式,初めて使いました。 こういう時にあると便利です。

*1:南裕樹:「Pythonのよる制御工学入門」,オーム社,2019年