【回路×Python】Lcapyを使って回路の伝達関数を求める方法

はじめに

　Pythonには線形回路を解析できるLcapyというパッケージがあります。
Lcapyには様々な機能がありますが、この記事では回路の伝達関数を記号的（数式的に）に求める方法について解説します。

なお、Lcapyを使えるようにする方法については以下の記事で解説しています。
【回路×Python】線形回路解析パッケージLcapyを使えるようにする方法

環境

Python: 3.7.4、SymPy: 1.6.2、Lcapy: 0.67.0

Lcapyはまだ発展途上のパッケージなので、上と同じverにしてもこの記事のコードが動作しない可能性はあります。
私の環境でも、公式ドキュメントのコードが動かないとかありますからね・・・。
あなたのコードに悪い箇所が無くても動かない可能性があるので、注意しましょう！

例１）RC回路

この記事ではCircuitメソッドを使用して伝達関数を求めます。
第１の例として下図のRC回路の伝達関数を解いてみましょう。

伝達関数を求めるコード：

例１－１）

from lcapy import *
cct = Circuit("""
Vi 1 0 
R 1 2 RF
C 2 0 CF
""")

H = (cct.C.V(s) / cct.Vi.V(s)).simplify()
H

実行結果：

$$\frac{1}{C_{F} R_{F} s + 1}$$

解説：
Circuit(""" """)の内部にネットリスト入力しています。
cct.C.V(s)でCの印加電圧を求め、これを入力電圧cct.Vi.V(s)で割ることで伝達関数を求めています。
伝達関数はsymplifyで単純化しています。
このように伝達関数が記号的に（数式的に）求められたことが分かります。

---

素子の値（例１－１におけるRFとCF）を省略して実行することも可能です。
この場合、素子値＝素子名（RとC）になります。

例１－２）

from lcapy import *
cct = Circuit("""
Vi 1 0 
R 1 2
C 2 0
""")

H = (cct.C.V(s) / cct.Vi.V(s)).simplify()
H

実行結果：

$$\frac{1}{C R s + 1}$$

---

当然、素子値を記号ではなく数値で入力することも可能です。

コード：

例１－３）

from lcapy import *
cct = Circuit("""
Vi 1 0 
R 1 2 1e3
C 2 0 1e-6
""")

H = (cct.C.V(s) / cct.Vi.V(s)).simplify()
H

実行結果：
$$\frac{1000}{s + 1000}$$

伝達関数を記号的に求め、その後に値を代入することも可能です。

例１－４）

from lcapy import *
cct = Circuit("""
Vi 1 0 
R 1 2 RF
C 2 0 CF
""")

H = (cct.C.V(s) / cct.Vi.V(s)).simplify()
H.subs('CF',1e-6).subs('RF',1e3)

実行結果：
$$\frac{1}{\frac{s}{1000} + 1}$$

例２）より複雑な受動回路

例１）のような単純な回路であればネットリストを自力で入力しても時間はかからないでしょう。
しかし、下図のような少し複雑な回路になると、自力入力は面倒です。

回路設計ツールで回路図を作成し、そこでネットリストを生成し、これをLcapyのコードに渡すことでネットリストを自力で書かずに済みます。

上図はLTSPICEで描きましたが、LTSPICEの場合View => SPICE NETLISTよりネットリストを取得できます。

ネットリストは以下のようにLcapyコードにペーストします。

例１－３）RC回路

from lcapy import *
cct = Circuit("""
V1 IN 0 s 1
R1 OUT N001 RP
C1 OUT 0 CM
L1 IN N001 LP
R2 OUT 0 RM
C2 OUT N001 CP
""")

H = cct["OUT"].V(s).simplify()
H

伝達関数は例１）とは異なる方法で求めています。まず、コピペしたネットリストにおいてV1を以下のように修正しています。
　　V1 IN 0 s 1
s領域において振幅=1の信号という意味です。
伝達関数Hはノード"OUT"の電圧を測定することで求めています。

実行結果：
$$\frac{R_{M} \left(C_{P} R_{P} s + 1\right)}{C_{M} C_{P} L_{P} R_{M} R_{P} s^{3} + L_{P} s^{2} \left(C_{M} R_{M} + C_{P} R_{P}\right) + R_{M} + R_{P} + s \left(C_{M} R_{M} R_{P} + C_{P} R_{M} R_{P} + L_{P}\right)}$$

例３）　有限ゲインのオペアンプ

有限ゲイン、シングル出力、周波数依存なし

下図のような非反転増幅器もLcapyで解析できます。

オペアンプとしてはVCVS（電圧制御電圧信号源）を使用します。
LTSPICEで描くと下図のようになります。

これも例２）と同様にネットリストをコピペして、V1を修正すれば伝達関数が求められます。

例３－１）

from lcapy import *
cct = Circuit("""
R1 N001 N002 RG
R2 OUT N001 RF
V1 N002 0 s 1
C1 N001 OUT CF
E1 OUT 0 0 N001 AOL
C2 N001 0 CG
""")

H = cct["OUT"].V(s).simplify()
H

実行結果：
$$- \frac{A_{\mathrm{OL}} R_{F}}{A_{\mathrm{OL}} R_{G} + R_{F} R_{G} s \left(A_{\mathrm{OL}} C_{F} + C_{F} + C_{G}\right) + R_{F} + R_{G}}$$

有限ゲイン、シングル出力、シングルポール

極を1つもつオペアンプで解析したい場合は下図のような回路図を作成すれば良いです。

オペアンプの特性は以下の式で表されるものとします。
$$A(f) = \frac{A_{OL}}{1+s/p}$$
オペアンプはVCCS（電圧制御電流源）、抵抗、容量、そしてVCVSを組み合わせて作成しています。
極は容量で決めており、$C={1}/p$です。
容量の値は数式で入力することが可能であり、数式は{}の内部に入力します。
このためLTSPICE上で容量値を{1/p}としています。

例３－２）

from lcapy import *
cct = Circuit("""
E1 0 OUT N002 0 1
C1 N002 0 {1/p}
R1 N002 0 1
R2 OUT INN RF
R3 INN N001 RG
V1 N001 0 s 1
G1 0 N002 0 INN AOL
""")

H = cct["OUT"].V(s).simplify()
H

実行結果：
$$- \frac{A_{\mathrm{OL}} R_{F} p}{A_{\mathrm{OL}} R_{G} p + R_{F} p + R_{G} p + s \left(R_{F} + R_{G}\right)}$$

有限ゲイン、完全差動、周波数依存なし

完全差動アンプは下図のような回路を使えばLcapyで計算することが可能です。

終わりに

回路方程式を解くことなしに簡単に伝達関数が求められるLcapyパッケージは非常に便利です。
伝達関数が得られた後はそれを展開したり（極を求めるなど）、数値解析をしたくなりますが、これについては別の記事で解説します。

【回路×Python】Lcapyを使って伝達関数を展開・計算する方法

※ついでに、伝達関数を近似によって単純化したくなることもあるかと思いますが、Lcapyには近似に関する機能は無いようです。このため私は自力で近似しています。