カスタム DC モーターの Simscape ライブラリのカスタマイズと拡張

ライブスクリプトを開く

Symbolic Math Toolbox™ を使用して、Simscape™ ライブラリ用のカスタム式ベースのコンポーネントを作成します。

はじめに

Symbolic Math Toolbox は、任意の空間次元において、工学第一原理からモデル開発を行うための柔軟な方法を提供しています。定常状態や過渡現象に関する数学モデルを開発できます。

コンポーネントに不可欠な物理特性を表す際に必要となる方程式を作成し、解を求めることができます。また、独自の低次元化されたモデルの、入力 x と対象となる量 f(x) との間のマッピングが可能です。

ここで、f は、次の形式で支配方程式を表すことができるカスタムコンポーネントです。

数式
ODE と PDE の数値シミュレーション

この例の手順は、以下です。

symReadSSCVariables を使用した Simscape™ コンポーネントのパラメーター化
diff を使用した Simscape コンポーネントのカスタム式の定義
solve と subs を使用した定常状態方程式の解析的な求解
matlabFunction と ode45 を使用した MATLAB® における時間依存方程式の数値的な求解
symWriteSSC を使用した Simscape コンポーネントの作成

この例を実行するには、Simscape および Symbolic Math Toolbox のライセンスが必要です。

DC モーターモデル

DC モーターは電気エネルギーを機械エネルギーに (またその逆方向にも) 変換するデバイスです。DC モーターの概略図を以下に示します (左図)。DC モーターのシミュレーションブロックはSimscape Electrical™ (右図) で与えられており、これは Simscape のオプション製品です。

この例では、支配方程式の常微分方程式 (ODE) を使用して DC モーターの低次元化モデル表現を導出します。DC モーターの電圧と電流はキルヒホッフの法則から、機械トルクの式はニュートンの法則から導かれます。これらの方程式を使用することで、カスタムのパラメトリックな Simscape コンポーネントを実装できます。

$(\begin{matrix} J \frac{\partial}{\partial t} w (t) = Ki I (t) - Dr w (t) \\ L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Kb w (t) \end{matrix})$

テンプレートコンポーネントのパラメーターと変数のインポート

現在のフォルダーまたは既定の MATLAB パスに Simscape コンポーネント MyMotorTemplate.ssc があるとします。このコンポーネントの方程式はまだありません。テンプレートにはモーターの開発に使用するパラメーターと変数が記録されています。type を使用すると、そのテンプレートのプレビューを表示できます。

type MyMotorTemplate.ssc

component MyMotorTemplate
% Custom DC Motor
% This block implements a custom DC motor
  nodes
    p = foundation.electrical.electrical; % +:left
    n = foundation.electrical.electrical; % -:left
    r = foundation.mechanical.rotational.rotational; % R:right
    c = foundation.mechanical.rotational.rotational; % C:right
  end
  parameters
      R = {3.9, 'Ohm'};                 %Armature resistance
      L = {0.000012, 'H'};              %Armature inductance
      J = {0.000001, 'kg*m^2'};         %Inertia
      Dr = {0.000003, '(N*m*s)/rad'};   %Rotor damping
      Ki = {0.000072, '(N*m)/A'};       %Torque constant 
      Kb = {0.000072, '(V*s)/rad'};     %Back-emf constant
  end
  variables
      torque = {0, 'N*m'};  %Total Torque
      tau = {0, 'N*m'};     %Electric Torque
      w = {0, 'rad/s'};     %Angular Velocity 
      I = {0, 'A'};         %Current
      V = {0, 'V'};         %Applied voltage
      Vb = {0, 'V'};        %Counter electromotive force
  end
  function setup
    if(R<=0)
        error('Winding resistance must be greater than 0.');
    end
  end
  branches
    torque : r.t -> c.t; % Through variable tau from r to c
    I   : p.i -> n.i; % Through variable i from p to n
  end
  equations
    w == r.w -c.w; % Across variable w from r to c
    V == p.v -n.v; % Across variable v from p to n
  end
end

テンプレートコンポーネントからパラメーターの名前、値および単位を読み取ります。

[parNames, parValues, parUnits] = symReadSSCParameters('MyMotorTemplate');

パラメーター、その値および対応する単位をベクトルで表示します。

vpa([parNames; parValues; parUnits],10)

ans = 
 $(\begin{array}{c} Dr & J & Kb & Ki & L & R \\ 0.000003 & 0.000001 & 0.000072 & 0.000072 & 0.000012 & 3.9 \\ \frac{N "newton - a physical unit of force." m "meter - a physical unit of length." s "second - a physical unit of time."}{rad "radian - a physical unit of plane angle."} & kg "kilogram - a physical unit of mass." {m "meter - a physical unit of length."}^{2} & \frac{V "volt - a physical unit of electric potential." s "second - a physical unit of time."}{rad "radian - a physical unit of plane angle."} & \frac{N "newton - a physical unit of force." m "meter - a physical unit of length."}{A "ampere - a physical unit of electric current."} & H "henry - a physical unit of inductance." & Ω "ohm - a physical unit of resistance." \end{array})$

関数 syms を使用してパラメーター名を MATLAB ワークスペースに追加します。パラメーターはワークスペースでシンボリック変数として表示されます。who を使用すると、ワークスペースの変数のリストを表示できます。

syms(parNames)
syms

Your symbolic variables are:

Dr   J    Kb   Ki   L    R    ans

コンポーネント変数の名前を読み取って表示します。ReturnFunction を使用して、これらの変数を変数 t の関数に同時に変換します。

[varFuns, varValues, varUnits] = symReadSSCVariables('MyMotorTemplate', 'ReturnFunction', true);
vpa([varFuns; varValues; varUnits],10)

ans = 
 $(\begin{array}{c} I (t) & V (t) & Vb (t) & τ (t) & torque (t) & w (t) \\ 0 & 0 & 0 & 0 & 0 & 0 \\ A "ampere - a physical unit of electric current." & V "volt - a physical unit of electric potential." & V "volt - a physical unit of electric potential." & N "newton - a physical unit of force." m "meter - a physical unit of length." & N "newton - a physical unit of force." m "meter - a physical unit of length." & \frac{rad "radian - a physical unit of plane angle."}{s "second - a physical unit of time."} \end{array})$

関数 syms を使用して変数名を MATLAB ワークスペースに追加します。変数はワークスペースでシンボリック関数として表示されます。必要なすべてのシンボリック変数とシンボリック関数を syms で宣言したことを確認します。

syms(varFuns)
syms

Your symbolic variables are:

Dr      I       J       Kb      Ki      L       R       V       Vb      ans     t       tau     torque  w

DC モーターをモデル化するための方程式系の定義

機械トルクの微分方程式は、eq1 および eq2 のように定義されます。I(t) は電流、w(t) は角速度を表します。

eq1 = torque + J*diff(w(t)) == -Dr*w(t) + tau(t)

eq1(t) = 
 $J \frac{\partial}{\partial t} w (t) + torque (t) = τ (t) - Dr w (t)$

eq2 = tau(t) == Ki*I(t)

eq2 =  $τ (t) = Ki I (t)$

電圧と電流の方程式は、eq3 および eq4 です。V(t) および Vb(t) はそれぞれ印加電圧と逆起電力を表します。

eq3 = L*diff(I(t)) + R*I(t) == V(t) - Vb(t)

eq3 = 
 $L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Vb (t)$

eq4 = Vb(t) == Kb*w(t)

eq4 =  $Vb (t) = Kb w (t)$

これらをまとめてリストにできます。ここで、モーターのトルクは、電流に比例すると仮定します。

eqs = formula([eq1; eq2; eq3; eq4])

eqs = 
 $(\begin{array}{c} J \frac{\partial}{\partial t} w (t) + torque (t) = τ (t) - Dr w (t) \\ τ (t) = Ki I (t) \\ L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Vb (t) \\ Vb (t) = Kb w (t) \end{array})$

方程式の左辺と右辺を抽出します。

operands = children(eqs);
operList = [ operands{:} ];
lhs = operList(1:2:end)

lhs=1×4 cell array
    {[J*diff(w(t), t) + torque(t)]}    {[tau(t)]}    {[L*diff(I(t), t) + R*I(t)]}    {[Vb(t)]}

rhs = operList(2:2:end)

rhs=1×4 cell array
    {[tau(t) - Dr*w(t)]}    {[Ki*I(t)]}    {[V(t) - Vb(t)]}    {[Kb*w(t)]}

2 番目と 4 番目の方程式は、tau(t) と Vb(t) の値を定義します。4 つの方程式の系を 2 つの方程式の系に単純化するには、これらの値を 1 番目と 3 番目の方程式に代入します。

equations(1) = subs(eqs(1), lhs(2), rhs(2))

equations = 
 $J \frac{\partial}{\partial t} w (t) + torque (t) = Ki I (t) - Dr w (t)$

equations(2) = subs(eqs(3), lhs(4), rhs(4))

equations = 
 $(\begin{array}{c} J \frac{\partial}{\partial t} w (t) + torque (t) = Ki I (t) - Dr w (t) & L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Kb w (t) \end{array})$

equations.'

ans = 
 $(\begin{array}{c} J \frac{\partial}{\partial t} w (t) + torque (t) = Ki I (t) - Dr w (t) \\ L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Kb w (t) \end{array})$

方程式を解く前にパラメーターに数値を代入します。また、V(t) = 1 を使用します。

equations = subs(equations, [parNames,V(t)], [parValues,1]);
equations = subs(equations, torque, 0);
vpa(equations.',10)

ans = 
 $(\begin{array}{c} 0.000001 \frac{\partial}{\partial t} w (t) = 0.000072 I (t) - 0.000003 w (t) \\ 0.000012 \frac{\partial}{\partial t} I (t) + 3.9 I (t) = 1.0 - 0.000072 w (t) \end{array})$

定常状態の方程式の解法

ここでは、関数 w(t) と関数 I(t) の時間依存性を取り除きます。たとえば、それらにシンボリック変数 ww および ii を代入します。

syms ww ii
equations_steady = subs(equations, [w(t),I(t)], [ww,ii]);

result = solve(equations_steady,ww,ii);
steadyStateW = vpa(result.ww,10)

steadyStateW =  $6.151120734$

steadyStateI = vpa(result.ii,10)

steadyStateI =  $0.2562966973$

`matlabFunction` と `ode45` を使用した MATLAB でのシンボリック式の数値的なシミュレーション

ode45 に有効な入力をシンボリック方程式から作成します。odeToVectorField を使用して、力学系 $\frac{dY}{dt} = f (t, Y)$ (初期条件 $Y (t_{0}) = Y_{0}$ ) を表す MATLAB プロシージャを作成します。

[vfEquations, tVals] = odeToVectorField(equations)

vfEquations = 
 $(\begin{array}{c} \frac{147573952589676412928}{1770887431076117} - 6 Y_{2} - \frac{2877692075498690052096 Y_{1}}{8854437155380585} \\ \frac{83010348331692984375 Y_{1}}{1152921504606846976} - \frac{27670116110564328125 Y_{2}}{9223372036854775808} \end{array})$

tVals = 
 $(\begin{array}{c} I \\ w \end{array})$

M = matlabFunction(vfEquations,'Vars', {'t','Y'})

M = function_handle with value:
    @(t,Y)[Y(1).*(-3.25e+5)-Y(2).*6.0+8.333333333333333e+4;Y(1).*7.2e+1-Y(2).*3.0]

初期条件 w(0) = 0 と I(0) = 0 を使用して微分方程式を解きます。

solution = ode45(M,[0 3],[0 0])

solution = struct with fields:
     solver: 'ode45'
    extdata: [1x1 struct]
          x: [0 2.4114e-09 1.4468e-08 7.4754e-08 3.7618e-07 1.8833e-06 4.0507e-06 6.6351e-06 9.6952e-06 1.3333e-05 1.7740e-05 2.3237e-05 3.0371e-05 4.0045e-05 5.2842e-05 6.2514e-05 7.2187e-05 8.3013e-05 9.4295e-05 1.0467e-04 ... ] (1x293775 double)
          y: [2x293775 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

t=[0.5,0.75,1] の時点で解を評価します。1 番目の値は電流 I(t)、2 番目の値は角速度 w(t) です。角速度の解が定常状態 steadyStateW に近づき始めていることがわかります。

deval(solution,0.5), deval(solution,.75), deval(solution,1)

ans = 2×1

    0.2563
    4.7795

ans = 2×1

    0.2563
    5.5034

ans = 2×1

    0.2563
    5.8453

steadyStateW

steadyStateW =  $6.151120734$

解をプロットします。

time  = linspace(0,2.5);
iValues = deval(solution, time, 1);
wValues = deval(solution, time, 2);
steadyStateValuesI = vpa(steadyStateI*ones(1,100),10);
steadyStateValuesW = vpa(steadyStateW*ones(1,100),10);

figure;
plot1 = subplot(2,1,1);
plot2 = subplot(2,1,2);

plot(plot1, time, wValues, 'blue',  time, steadyStateValuesW, '--red', 'LineWidth', 1)
plot(plot2, time, iValues, 'green', time, steadyStateValuesI, '--red', 'LineWidth', 1)

title(plot1, 'DC Motor - angular velocity')
title(plot2, 'DC Motor - current')
ylabel(plot1, 'Angular velocity [rad/s]')
ylabel(plot2, 'Current [A]')
xlabel(plot1, 'time [s]')
xlabel(plot2, 'time [s]')
legend(plot1, 'w(t)', 'w(t): steady state', 'Location', 'northeastoutside')
legend(plot2, 'I(t)', 'I(t): steady state', 'Location', 'northeastoutside')

Simscape で使用するための数学モデルの保存

作成した方程式 eqs を使用して Simscape コードを生成します。

symWriteSSC('MyMotor.ssc', 'MyMotorTemplate.ssc', eqs, ...
    'H1Header', '% Custom DC Motor', ...
    'HelpText', {'% This block implements a custom DC motor'})

生成されたコンポーネントを type コマンドで表示します。

type MyMotor.ssc

component MyMotor
% Custom DC Motor
% This block implements a custom DC motor
  nodes
    p = foundation.electrical.electrical; % +:left
    n = foundation.electrical.electrical; % -:left
    r = foundation.mechanical.rotational.rotational; % R:right
    c = foundation.mechanical.rotational.rotational; % C:right
  end
  parameters
      R = {3.9, 'Ohm'};                 %Armature resistance
      L = {0.000012, 'H'};              %Armature inductance
      J = {0.000001, 'kg*m^2'};         %Inertia
      Dr = {0.000003, '(N*m*s)/rad'};   %Rotor damping
      Ki = {0.000072, '(N*m)/A'};       %Torque constant 
      Kb = {0.000072, '(V*s)/rad'};     %Back-emf constant
  end
  variables
      torque = {0, 'N*m'};  %Total Torque
      tau = {0, 'N*m'};     %Electric Torque
      w = {0, 'rad/s'};     %Angular Velocity 
      I = {0, 'A'};         %Current
      V = {0, 'V'};         %Applied voltage
      Vb = {0, 'V'};        %Counter electromotive force
  end
  function setup
    if(R<=0)
        error('Winding resistance must be greater than 0.');
    end
  end
  branches
    torque : r.t -> c.t; % Through variable tau from r to c
    I   : p.i -> n.i; % Through variable i from p to n
  end
  equations
    w == r.w -c.w; % Across variable w from r to c
    V == p.v -n.v; % Across variable v from p to n
    torque+J*w.der == tau-Dr*w;
    tau == Ki*I;
    L*I.der+R*I == V-Vb;
    Vb == Kb*w;
  end
end

生成されたコンポーネントから Simscape ライブラリを作成します。

if ~isdir('+MyLib')
    mkdir +MyLib;
end
copyfile MyMotor.ssc +MyLib;
ssc_build MyLib;

Generating Simulink library 'MyLib_lib' in the current directory '/tmp/Bdoc24a_2550221_459077/tp3293a62c/symbolic-ex98670381' ...

参考

関数

diff | subs | solve | matlabFunction | ode45 | symReadSSCParameters | symReadSSCVariables | symWriteSSC