トランジスタのスティッフな微分代数方程式の求解

この例では、電気回路を記述するスティッフな微分代数方程式 (DAE) を解く方法を説明します [1]。例のファイル amp1dae.m にコード記述されている 1 石トランジスタ アンプの問題は、半陽的な形式に書き換えることもできますが、この例では元の形式 $$M{u}^{\prime }=\varphi (u)$$ で解きます。この問題には、定数の特異質量行列 $\mathit{M}$ が含まれます。

このトランジスタ アンプ回路には 6 つの抵抗、3 つのコンデンサ、1 つのトランジスタが含まれます。

• 初期の電圧信号は $$U_e(t)=0.4\sin \left(200\pi t\right)$$ です。

• 動作電圧は ${\mathit{U}}_{\mathit{b}}=6$ です。

• ノードの電圧は ${\mathit{U}}_{\mathit{i}}\left(\mathit{t}\right)\left(\mathit{i}=1,2,3,4,5\right)$ によって与えられます。

• 抵抗の値 ${\mathit{R}}_{\mathit{i}}\left(\mathit{i}=1,2,3,4,5,6\right)$ は定数であり、各抵抗を通る電流は $\mathit{I}=\mathit{U}/\mathit{R}$ を満たします。

• コンデンサの値 ${\mathit{C}}_{\mathit{i}}\left(\mathit{i}=1,2,3\right)$ は定数であり、各コンデンサを通る電流は $\mathit{I}=\mathit{C}\cdot \mathrm{dU}/\mathrm{dt}$ を満たします。

目的はノード 5 にかかる出力電圧 ${\mathit{U}}_5\left(\mathit{t}\right)$ を求めることです。

MATLAB® でこの方程式を解くには、ode オブジェクトを使用し、オブジェクトのプロパティを設定して方程式、質量行列、初期状態を定義します。その後、solve メソッドを使用して、連立方程式を経過的にシミュレートします。(ここで行ったように) 必要な関数をスクリプト内に関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パスのディレクトリに保存することもできます。

質量行列のコード化

キルヒホッフの法則を使用して各ノード (1 ～ 5) を通る電流を等しくすると、回路を記述する 5 つの方程式からなる連立方程式が得られます。

$$\begin{array}{l}\mathrm{node}1:\frac{{\mathit{U}}_{\mathit{e}}\left(\mathit{t}\right)}{{\mathit{R}}_0}-\frac{{\mathit{U}}_1}{{\mathit{R}}_0}+{\mathit{C}}_1\left({{\mathit{U}}_2}^{\prime }-{{\mathit{U}}_1}^{\prime }\right)=0,\\ \mathrm{node}2:\frac{{\mathit{U}}_{\mathit{b}}}{{\mathit{R}}_2}-{\mathit{U}}_2\left(\frac{1}{{\mathit{R}}_1}+\frac{1}{{\mathit{R}}_2}\right)+{\mathit{C}}_1\left({{\mathit{U}}_1}^{\prime }-{{\mathit{U}}_2}^{\prime }\right)-0.01\mathit{f}\left({\mathit{U}}_2-{\mathit{U}}_3\right)=0,\\ \mathrm{node}3:\mathit{f}\left({\mathit{U}}_2-{\mathit{U}}_3\right)-\frac{{\mathit{U}}_3}{{\mathit{R}}_3}-{\mathit{C}}_2{{\mathit{U}}_3}^{\prime }=0,\\ \mathrm{node}4:\frac{{\mathit{U}}_{\mathit{b}}}{{\mathit{R}}_4}-\frac{{\mathit{U}}_4}{{\mathit{R}}_4}+{\mathit{C}}_3\left({{\mathit{U}}_5}^{\prime }-{{\mathit{U}}_4}^{\prime }\right)-0.99\mathit{f}\left({\mathit{U}}_2-{\mathit{U}}_3\right)=0,\\ \mathrm{node}5:-\frac{{\mathit{U}}_5}{{\mathit{R}}_5}+{\mathit{C}}_3\left({{\mathit{U}}_4}^{\prime }-{{\mathit{U}}_5}^{\prime }\right)=0.\end{array}$$

この連立方程式の質量行列は、方程式の左辺の微分項を集めることによって求められ、次の形式になります。

$$M=\left(\begin{array}{ccccc}-c_1& c_1& 0& 0& 0\\ c_1& -c_1& 0& 0& 0\\ 0& 0& -c_2& 0& 0\\ 0& 0& 0& -c_3& c_3\\ 0& 0& 0& c_3& -c_3\end{array}\right),$$

ここで、$$k=1,2,3$$ に対して $$c_k=k\times 10^{-6}$$ です。

適切な定数 ${\mathit{c}}_{\mathit{k}}$ を使って質量行列を作成し、その質量行列を格納するための odeMassMatrix オブジェクトを作成します。質量行列が特異行列であることは明らかですが、ソルバーによる DAE 問題の自動検出をテストするために、Singular プロパティを既定値 "maybe" のままにします。

c = 1e-6 * (1:3);
M = zeros(5,5);
M(1,1) = -c(1);
M(1,2) =  c(1);
M(2,1) =  c(1);
M(2,2) = -c(1);
M(3,3) = -c(2);
M(4,4) = -c(3);
M(4,5) =  c(3);
M(5,4) =  c(3);
M(5,5) = -c(3);
Mass = odeMassMatrix(MassMatrix=M);

方程式のコード作成

関数 transampdae にはこの例の連立方程式が含まれています。この関数は、すべての電圧および定数パラメーターの値を定義します。方程式の左辺に集められた導関数は質量行列にコード化され、transampdae は方程式の右辺をコード化します。

function dudt = transampdae(t,u)
% Define voltages and parameters
Ue = @(t) 0.4*sinpi(200*t);
Ub = 6;
R0 = 1000;
R15 = 9000;
alpha = 0.99;
beta = 1e-6;
Uf = 0.026;

% Define system of equations
f23 = beta*(exp((u(2) - u(3))/Uf) - 1);
dudt = [ -(Ue(t) - u(1))/R0
    -(Ub/R15 - u(2)*2/R15 - (1-alpha)*f23)
    -(f23 - u(3)/R15)
    -((Ub - u(4))/R15 - alpha*f23)
    (u(5)/R15) ];
end

初期条件のコード化

初期条件を設定します。この例では、[1] で計算された各ノードを通る電流に対し、整合性のある初期条件を使用します。

Ub = 6;
u0(1) = 0;
u0(2) = Ub/2;
u0(3) = Ub/2;
u0(4) = Ub;
u0(5) = 0;

連立方程式を解く

ode オブジェクトを作成して問題を記述し、ODEFcn、MassMatrix、および InitialValue の各プロパティの値を設定します。さらに、問題を解くためのソルバーとして ode23t を選択します。

F = ode(ODEFcn=@transampdae,MassMatrix=Mass,InitialValue=u0,Solver="ode23t")

F = 
  ode with properties:

   Problem definition
               ODEFcn: @transampdae
          InitialTime: 0
         InitialValue: [0 3 3 6 0]
             Jacobian: []
           MassMatrix: [1×1 odeMassMatrix]
         EquationType: standard

   Solver properties
    AbsoluteTolerance: 1.0000e-06
    RelativeTolerance: 1.0000e-03
               Solver: ode23t

  Show all properties

solve メソッドを使用して、時間区間 [0 0.05] での DAE 系の解を求めます。

S = solve(F,0,0.05);
t = S.Time;
u = S.Solution;

結果のプロット

初期電圧 $$U_e(t)$$ とノード 5 にかかる出力電圧 $$U_5(t)$$ をプロットします。

Ue = @(t) 0.4*sinpi(200*t);
plot(t,Ue(t),"o",t,u(5,:),".")
axis([0 0.05 -3 2]);
legend("Input Voltage U_e(t)","Output Voltage U_5(t)",Location="NorthWest");
title("One Transistor Amplifier DAE Problem Solved by ODE23T");
xlabel("t");

参照

[1] Hairer, E., and Gerhard Wanner. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer Berlin Heidelberg, 1991, p. 377.

参考

ode23t | ode15s

トピック

• ODE ソルバーの選択
• ODE オプションの概要
• 微分代数方程式 (DAE) の求解
• 状態依存の強い質量行列による ODE の求解