경직성(Stiff) 트랜지스터 미분대수 방정식 풀기

이 예제에서는 전기 회로를 나타내는 경직성 미분대수 방정식(DAE)을 푸는 방법을 보여줍니다[1]. 예제 파일 amp1dae.m에 작성된 단일 트랜지스터 증폭기 문제는 반명시적(semi-explicit) 형식으로 재작성할 수 있지만, 이 예제에서는 원래 형식 $$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 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.

질량 행렬 코딩하기

키르히호프의 법칙(Kirchoff's law)에 따라 각 노드(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;

연립방정식 풀기

ODEFcn, MassMatrix, InitialValue 속성의 값을 설정하여 문제를 설명하는 ode 객체를 만듭니다. 문제를 푸는 데 사용할 솔버로 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: [1x1 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;

결과 플로팅하기

노드 5인 $$U_e(t)$$까지 초기 전압 $$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 풀기