コンバータ制御設計に役立つ状態平均化モデルをsimulinkで書いてみた1

このシリーズの記事一覧

記事その２https://qiita.com/griffin921/items/53dc2f3218bcd8d8c188
記事その３https://qiita.com/griffin921/items/4078d4bb6e1302d6237c
記事その４https://qiita.com/griffin921/items/87dcda9ba1879cca1181

コンバータってどうやって制御すんの？

制御の基本はまず制御対象のモデリングから始まると言っても過言ではない（犬調べ）。
対象の挙動がよくわかっていれば、それに合わせた制御構成・パラメータの設定が可能になるし、
制御対象がもつばらつきなどの変動要素に対して構築しようとしているシステムがどれだけ堅牢なのかを試作する前に机上で検討することができてとても嬉しい。

しかし。コンバータはいわゆる非線形なプラントだと言われる。
これはパワー半導体のスイッチ切り替えによって回路の経路が常に高速に切り替わっているためである。
(もちろん非線形な要素はこれだけではなくて、ダイオードやインダクタ等の特性、負荷特性によっても伝達関数が変わるのですべて列挙するとキリがないが　ここで言う非線形性は切り替わっていることによる非線形性のことを言ってます）
なのでみんなFRA (Frequency Response Analyzer) で回路の特性を計測したりする。

最近はもっと安い機材も出てきているが、
きっとエヌエフ回路ブロックのFRAにはお世話になった人もいるはずである。
https://www.nfcorp.co.jp/pro/category/mi/fra/

でも、基本的にはこういう計測はあくまで事後評価であるので、実機を作る前に伝達関数を取得することができない。
それは困るわけである。我々は試作する前に大体の設計可能範囲を見積もりたい。

状態平均化法

そこで登場するのが状態平均化法という理論である。
ググるといっぱい出てくるので調べてみていただきたい。

後述するReferenceに記載がある[1],[2]の書籍で詳細は確認していただきたいのですが、
ざっくりいうと

• パワー半導体がスイッチONのときとOFFのときの回路微分方程式を求める
• それぞれの微分方程式をDuty比率で重み付けして足し合わせると平均化された式が取得できる
• そこで求めた式を制御対象として制御組むと　あら不思議！うまくいく

・・・はい。ざっくりしすぎて専門家に見られるとぶん殴られると思いますが。
できることはこういうことです。

で、ここではなぜこの方法でうまくいくのかが非常に大事ですが、
そこを気になる方は引用してる本をぜひ買ってください。買ったほうがいいですほんと。

基本的には本来電源が制御したいのは電圧や電流のリプルを無視した平均値であるので
そこに着目したモデルを立てることが大切なのかなとざっくり理解してますが、
一度理論を抑えてから設計したほうが、この例題とは違う回路を検討する際に、うまくいかなかったときの考察がやりやすいのでおすすめです。

ここはあくまで方法論をお伝えするだけという風にご理解いただければ幸いです。

状態平均化モデルの算出

MATALBのSymbolic Mathを使ってバックコンバータの状態平均化モデルを算出してみました。
色々理論を書き綴っていますが参考にしてくだしあ！

この記事は３章構成になりそうで、最後の章はAdvent Calenderに乗っけようと思います。よろしくお願いいたします。

1．対象回路の定義

図１にBuckコンバータの回路図を示す。入力電圧は$V_{\textrm{in}}$とし、出力電圧は負荷抵抗(Ro)の両端電圧を$v_{\textrm{out}}$とする。

図１：Buckコンバータ回路図

2．スイッチがON状態での状態空間モデルを求める

まず、スイッチがON状態の微分方程式を立てる。

スイッチQ1がONの際に回路に流れる電流は下記図２のとおりになる。

図２：スイッチON時におけるBuckコンバータ回路と電流経路

この回路を元に微分方程式を立てていく。

まず入力電圧$V_{\textrm{in}}$はインダクタ両端に発生する電圧$v_L$と出力電圧$v_{\textrm{out}}$との和であるため、

V_{\textrm{in}\;} =v_L +v_{\textrm{out}} \cdots \left(2\ldotp 1\right)

$v_L$はインダクタ電流$i_L$を微分してインダクタンス$L$をかけたものと等しい。そのため下記式に書き換えることができる。

V_{\textrm{in}} =L_1 \frac{{\textrm{di}}_L }{\textrm{dt}}+v_{\textrm{out}} \cdots \left(2\ldotp 2\right)

次にコンデンサに流入する電流$i_C$を求める。コンデンサの両端電圧と流れる電流の関係式から

i_C =C_1 \frac{{\textrm{dv}}_C }{\textrm{dt}}\cdots \left(2\ldotp 3\right)

図２から、コンデンサには流入する電流$i_L$と流出する電流$I_{\textrm{out}}$が存在する。

キルヒホッフの電流則に従うと、コンデンサに流入する電流$i_C$は流出電流と流入電流の差で求められる。よって

i_C =i_L -I_{\textrm{out}} \cdots \left(2\ldotp 4\right)

$I_{\textrm{out}}$は

I_{\textrm{out}} =\frac{v_{\textrm{out}} }{R_o }\cdots \left(2\ldotp 5\right)

また、 出力電圧はコンデンサ電圧と等しいため

v_{\textrm{out}\;} =v_C \cdots \left(2\ldotp 6\right)

上記の式(2.2),(2.3),(2.4),(2.5),(2.6)を元に次式を得る。

\frac{{\textrm{di}}_L }{\textrm{dt}}=\frac{1}{L_1 }V_{\textrm{in}} -\frac{1}{L_1 }v_C \cdots \left(2\ldotp 7\right)

\frac{{d\;v}_{C\;} }{\textrm{dt}}=\frac{1}{C_1 }i_L -\frac{1}{{C_1 R}_O }v_C \cdots \left(2\ldotp 8\right)

インダクタ電流$i_L$とコンデンサ電圧$v_C$を状態変数として状態空間モデルとして整理し直す。

出力となる状態方程式はコンバータの出力電圧である$v_{\textrm{out}}$であるが式(2.6)よりすでに導出済みである。

微分方程式を元に各状態空間行列を得る。

\mathit{\mathbf{x}}\left(t\right)=\left\lbrack \begin{array}{c}
i_L \\
v_C 
\end{array}\right\rbrack　\cdots \left(2\ldotp 9\right)

\mathit{\mathbf{u}}\left(t\right)=V_{\textrm{in}}　\cdots \left(2\ldotp 10\right)

\begin{array}{l}
{\mathit{\mathbf{A}}}_{\textrm{on}} =\left\lbrack \begin{array}{cc}
0 & -\frac{1}{L_1 }\\
\frac{1}{C_1 } & -\frac{1}{{C_1 R}_O }
\end{array}\right\rbrack \cdots \left(2\ldotp 11\right)\\
{\mathit{\mathbf{B}}}_{\textrm{on}} =\left\lbrack \begin{array}{cc}
\frac{1}{L_1 } & 0\\
0 & 0
\end{array}\right\rbrack \cdots \left(2\ldotp 12\right)\\
\begin{array}{l}
{\mathit{\mathbf{C}}}_{\textrm{on}} =\left\lbrack \begin{array}{c}
0\\
1
\end{array}\right\rbrack \cdots \left(2\ldotp 13\right)\\
{\mathit{\mathbf{D}}}_{\textrm{on}} =\mathrm{０}\cdots \left(2\ldotp 14\right)
\end{array}
\end{array}

よって状態空間モデルは下記式となる

\begin{array}{l}
\overset{\cdot }{\mathit{\mathbf{x}}} \left(t\right)={\mathit{\mathbf{A}}}_{\textrm{on}} \mathit{\mathbf{x}}\left(t\right)+{\mathit{\mathbf{B}}}_{\textrm{on}} \mathit{\mathbf{u}}\left(t\right)\cdots \left(2\ldotp 15\right)\\
v_{\textrm{out}} \left(t\right)={\mathit{\mathbf{C}}}_{\textrm{on}} \mathit{\mathbf{x}}\left(t\right)+{\mathit{\mathbf{D}}}_{\textrm{on}} \mathit{\mathbf{u}}\left(t\right)\cdots \left(2\ldotp 16\right)
\end{array}

式(2.15),(2.16)に状態空間行列を代入すると

\frac{\mathrm{d}}{\mathrm{d}t}\left\lbrack \begin{array}{c}
i_L \\
v_C 
\end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
0 & -\frac{1}{L_1 }\\
\frac{1}{C_1 } & -\frac{1}{{C_1 R}_O }
\end{array}\right\rbrack \left\lbrack \begin{array}{c}
i_L \\
v_C 
\end{array}\right\rbrack +\left\lbrack \begin{array}{cc}
\frac{1}{L_1 } & 0\\
0 & 0
\end{array}\right\rbrack \textrm{Vin}\cdots \left(2\ldotp 17\right)

v_{\textrm{out}} \left(t\right)=v_c \left(t\right)\cdots \left(2\ldotp 18\right)

となる。

MATLABではSynbolic MathToolboxを用いると代数式を記述することが可能になる。

状態空間行列を記述すると下記のようになる。

今回はBuck Converterという比較的シンプルな回路に適応したため、検算するほどでもないかもしれないが、

より複雑な回路を同様のスキームで解析する際に、理論式の妥当性を確認するために計算機の計算結果を参考にすることは有効であると考える。

Code

clear variables
%代数シンボルを定義
syms L1 C1 Ro;

%状態空間行列を定義
Aon = [0 -1/L1;1/C1 -1/C1/Ro]

Aon =

 \left(\begin{array}{cc}
0 & -\frac{1}{L_1 }\\
\frac{1}{C_1 } & -\frac{1}{C_1 \,\textrm{Ro}}
\end{array}\right)

Code

Bon = [1/L1; 0]

Bon =

 \left(\begin{array}{c}
\frac{1}{L_1 }\\
0
\end{array}\right)

Code

Con = [0 1]

Output

Con = 1x2    
     0     1

Code

Don = 0

Output

Don = 0

3．スイッチがOFF状態での状態空間モデルを求める

次に、スイッチがOFF状態の微分方程式を立てる。

スイッチQ1がOFFの際に回路に流れる電流は下記図３のとおりになる。

入力電圧がスイッチQ1によって切断され、代わりにダイオードD1が導通するためインダクタ電流経路が変化する。

図3：スイッチOFF時におけるBuckコンバータ回路と電流経路

この回路を元に微分方程式を立てていく。

まず入力電圧$V_{\textrm{in}}$はスイッチによって切断されているため、

0=v_L +v_{\textrm{out}} \cdots \left(3\ldotp 1\right)

次にコンデンサに流入する電流$i_C$を求める。

コンデンサの両端電圧と流れる電流の関係式から(2.3)式と同様の式を得る。

i_C =C_1 \frac{{\textrm{dv}}_C }{\textrm{dt}}\cdots \left(3\ldotp 2\right)

図３から、コンデンサには流入する電流$i_L$と流出する電流$I_{\textrm{out}}$の関係は図２と変わらないため

以降立式すべき式は(2.4),(2.5),(2.6)とおなじになる。

これらの式をもとに微分方程式を立てると

\frac{{\textrm{di}}_L }{\textrm{dt}}=-\frac{v_C }{L_1 }\cdots \left(3\ldotp 3\right)

\frac{{d\;v}_{C\;} }{\textrm{dt}}=\frac{1}{C_1 }i_L -\frac{1}{C_1 R_o }v_C \cdots \left(3\ldotp 4\right)

インダクタ電流$i_L$とコンデンサ電圧$v_C$を状態変数として

一階微分方程式を導出する。

\frac{\mathrm{d}}{\mathrm{d}t}\left\lbrack \begin{array}{c}
i_L \\
v_C 
\end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
0 & -\frac{1}{L_1 }\\
\frac{1}{C_1 } & -\frac{1}{{C_1 R}_O }
\end{array}\right\rbrack \left\lbrack \begin{array}{c}
i_L \\
v_C 
\end{array}\right\rbrack \cdots \left(3\ldotp 5\right)

出力となる状態方程式はコンバータの出力電圧である$v_{\textrm{out}}$であるが式(2.6)よりすでに導出済みである。

式(3.5)を元に各状態空間行列を得る。

\mathit{\mathbf{x}}\left(t\right)=\left\lbrack \begin{array}{c}
i_L \\
v_C 
\end{array}\right\rbrack 
\cdots \left(3\ldotp 6\right)

\mathit{\mathbf{u}}\left(t\right)=V_{\textrm{in}}
\cdots \left(3\ldotp 7\right)

Symbolic Mathで代数式を記述すると下記のようになる。

Code

%状態空間行列を定義
Aoff=[0 -1/L1;1/C1 -1/C1/Ro]

Aoff =

 \left(\begin{array}{cc}
0 & -\frac{1}{L_1 }\\
\frac{1}{C_1 } & -\frac{1}{C_1 \,\textrm{Ro}}
\end{array}\right)

Code

Boff = 0

Output

Boff = 0

Code

Coff = [0 1]

Output

Coff = 1x2    
     0     1

Code

Doff = 0

Output

Doff = 0

4．Dutyで重み付けした状態平均モデルの算出

Symbolic MathでOn時、Off時の状態空間モデルが取得できているので

ON時間およびOFF時間の時比率$\textrm{Duty}$で平均化した各状態空間行列を計算する。

これにより、平均化した状態平均化方程式、即ち状態平均モデル式(4.1)と出力方程式(4.2)を下記の通り取得することができる。

\begin{array}{l}
\frac{d}{\mathrm{d}t}\mathit{\mathbf{x}}\left(t\right)=\textrm{Ax}\left(t\right)+\textrm{Bu}\left(t\right)\cdots \left(4\ldotp 1\right)\\
v_{\textrm{out}} \left(t\right)=\textrm{Cx}\left(t\right)+\textrm{Du}\left(t\right)\cdots \left(4\ldotp 2\right)
\end{array}

Code

%シンボルを追加定義
syms Duty dxdt x vout u iL vc Vin

x = [iL; vc] 

x =

 \left(\begin{array}{c}
\textrm{iL}\\
\textrm{vc}
\end{array}\right)

Code

u = Vin

u =

 \textrm{Vin}

Code

A = simplify(Aon*Duty + Aoff*(1-Duty))

A =

 \left(\begin{array}{cc}
0 & -\frac{1}{L_1 }\\
\frac{1}{C_1 } & -\frac{1}{C_1 \,\textrm{Ro}}
\end{array}\right)

Code

B = simplify(Bon*Duty + Boff*(1-Duty))

B =

 \left(\begin{array}{c}
\frac{\textrm{Duty}}{L_1 }\\
0
\end{array}\right)

Code

C = simplify(Con*Duty + Coff*(1-Duty))

C =

 \left(\begin{array}{cc}
0 & 1
\end{array}\right)

Code

D = simplify(Don*Duty + Doff*(1-Duty))

D =

 0

Code

dxdt = simplify(A*x + B*u)

dxdt =

 \left(\begin{array}{c}
-\frac{\textrm{vc}-\textrm{Duty}\,\textrm{Vin}}{L_1 }\\
-\frac{\textrm{vc}-\textrm{Ro}\,\textrm{iL}}{C_1 \,\textrm{Ro}}
\end{array}\right)

Code

vout = C*x + D*u

vout =

 \textrm{vc}

続きます！

次はこちら。
https://qiita.com/griffin921/items/53dc2f3218bcd8d8c188
作成した状態平均化モデルは後で使うので保存しておきましょう

Code

%変数を保存
save('A','B','C','D','dxdt');

Reference

[1] 安部　征哉, 財津　俊行, 上松 武.　デジタル電源の基礎と設計法 -スイッチング電源のデジタル制御-. 科学情報出版株式会社, 2020,246p, ISBN978-4-904774-84-7.

[2] Frede Blaabjerg, ed. Yanfelng Shen, Zian Qi and Huai Wang. Control of Power Electronic Converters And Systems. Volume 1, Chapter 3 "Modeling and Control of DC-DC Converters" p69-p92, Academic Press, 2018, ISBN 978-0-12-804254-7.

[3]MathWorks."Symbolic Math Toolbox 入門". MathWorks ヘルプセンター　ヘルプドキュメント,2022/10/19. https://jp.mathworks.com/help/symbolic/getting-started-with-symbolic-math-toolbox.html