続・LTspiceによる電流制御のシミュレーション
はじめに
先日の記事で,LTspiceによる電流制御のシミュレーションについて述べました。 negligible.hatenablog.com
この記事の中で私個人としてひとつの発見と言えるのは,
- 操作量に対してむだ時間T,フィードバック量に対して窓Tの移動平均フィルタを施した連続時間制御システム
- 操作量に対してむだ時間TとPWM,フィードバック量に対してサンプリング周期Tのサンプル & ホールドを設けた離散時間制御システム
がほぼ同じ挙動を示すということです。 図1に上記前者のブロック図を示します。
言い換えれば,離散時間制御 + PWMの電流制御系を,移動平均フィルタとむだ時間のある連続制御系で近似できると言えるのではないでしょうか(PWMの過変調が起きない範囲で)。 例えば図2ですね。
ほぼ挙動は一致しています。 そこで,連続系である図1の安定性について検討してみます。 また,最近はLTspiceばかり触っているので,JupyterLabでボード線図を描いてみました。
回路定数は先日の記事と同じくR = 20 mΩ,L = 5 mHとします。 また,同じくサンプリング時間 = 移動平均フィルタの窓 = むだ時間T = 100 μsとします。
何が問題か
先日の記事では,むだ時間も移動平均もないブロック図(図4)について, 開ループ伝達関数(一巡伝達関数)Go(s)が次式で表されること,
および,
を満たすように比例ゲインKP,積分ゲインKIを設定できれば,極零相殺されて
となることを述べました。
(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)式の再掲です。
ここで,
として,パラメータKを変更しながらGi(s),Gs(s)のボード線図を描いてみましょう。 なお,図5 ~ 7において"Ideal"はGi(s),"Sampled"はGs(s)を示します。
図5 ~ 7は例としてそれぞれパラメータK = 1,K = 3,K = 4.5の場合を描いています。 それぞれGs(s)の位相余裕PMが68.57°,26.98°,−2.12°となりました。 そう! K = 4.5の場合には不安定になります! 安定不安定の境界はK = 4.2付近にあるようです。
LTspiceでのシミュレーション
さて,それではK = 1, 3, 4.5の場合をシミュレーションしてみます。 モデル(schematic)は同じく図7となります。
- 連続時間制御 + PWMなし
- 離散時間制御 + PWMあり
の2つを比較します。 PWMする場合,直流電圧250 Vの単相フルブリッジインバータを想定します。
K = 1の場合はまさにKP = L / (4 T)の場合であり,理論上はステップ応答にオーバーシュートが発生しない波形となります。 位相余裕は60°以上あり,申し分ない(?)安定した波形となります。
K = 3の場合はオーバーシュートが発生しますが,位相余裕が約27°あり,まだ安定です。 ただし,文献*1の6.3節の最後,p. 223,表6.1にて,きなみ・みゆう姉妹の会話の中で示されている通り,サーボ系としては位相余裕40° ~ 60°は欲しいところです。
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()
ラムダ式,初めて使いました。 こういう時にあると便利です。