ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

Symbolic Math Toolbox を使用して、Simscape ライブラリにカスタム方程式に基づいたコンポーネントを作成します。

はじめに

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

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

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

  • 数式

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

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

  • symReadSSCVariables を使用した Simscape コンポーネントのパラメーター化

  • diff を使用した Simscape コンポーネントのカスタム式の定義

  • solvesubs を使用した定常状態方程式の解析的な求解

  • matlabFunctionode45 を使用した MATLAB における時間依存方程式の数値的な求解

  • symWriteSSC を使用した Simscape コンポーネントの作成

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

DC モーター モデル

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

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

(Jt w(t)=KiI(t)Drw(t)Lt I(t)+RI(t)=V(t)Kbw(t))

Simscape コンポーネントのパラメーター化

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

現在のフォルダーまたは既定の 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 = 

(DrJKbKiLR0.0000030.0000010.0000720.0000720.0000123.91N"ニュートン - 力の物理単位です。"m"メートル - 長さの物理単位です。"s"秒 - 時間の物理単位です。"rad"ラジアン - 平面角の物理単位です。"1kg"キログラム - 質量の物理単位です。"m"メートル - 長さの物理単位です。"21V"ボルト - 電位の物理単位です。"s"秒 - 時間の物理単位です。"rad"ラジアン - 平面角の物理単位です。"1N"ニュートン - 力の物理単位です。"m"メートル - 長さの物理単位です。"A"アンペア - 電流量の物理単位です。"H"ヘンリー - インダクタンスの物理単位です。"Ω"オーム - 抵抗の物理単位です。")

関数 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 = 

(I(t)V(t)Vb(t)τ(t)torque(t)w(t)000000A"アンペア - 電流量の物理単位です。"V"ボルト - 電位の物理単位です。"V"ボルト - 電位の物理単位です。"1N"ニュートン - 力の物理単位です。"m"メートル - 長さの物理単位です。"1N"ニュートン - 力の物理単位です。"m"メートル - 長さの物理単位です。"1rad"ラジアン - 平面角の物理単位です。"s"秒 - 時間の物理単位です。")

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

syms(varFuns)
syms
Your symbolic variables are:

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

Simscape コンポーネントのカスタム式の定義

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

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

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

Jt w(t)+torque(t)=τ(t)-Drw(t)

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

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

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

Lt I(t)+RI(t)=V(t)-Vb(t)

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

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

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

(Jt w(t)+torque(t)=τ(t)-Drw(t)τ(t)=KiI(t)Lt I(t)+RI(t)=V(t)-Vb(t)Vb(t)=Kbw(t))

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

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

(Jt w(t)+torque(t)τ(t)Lt I(t)+RI(t)Vb(t))

rhs = operList(2:2:end)
rhs = (τ(t)-Drw(t)KiI(t)V(t)-Vb(t)Kbw(t))

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

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

Jt w(t)+torque(t)=KiI(t)-Drw(t)

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

(Jt w(t)+torque(t)=KiI(t)-Drw(t)Lt I(t)+RI(t)=V(t)-Kbw(t))

equations.'
ans = 

(Jt w(t)+torque(t)=KiI(t)-Drw(t)Lt I(t)+RI(t)=V(t)-Kbw(t))

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

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

(0.000001t w(t)=0.000072I(t)-0.000003w(t)0.000012t I(t)+3.9I(t)=1.0-0.000072w(t))

定常状態方程式の解析的な求解

定常状態の方程式を解く

ここでは、関数 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

時間依存方程式の数値的な求解

matlabFunctionode45 を使用して MATLAB でシンボリック式を数値的にシミュレートします。

ode45 に有効な入力をシンボリック方程式から作成します。odeToVectorField を使用して、力学系 dYdt=f(t,Y) (初期条件 Y(t0)=Y0) を表す MATLAB プロシージャを作成します。

[vfEquations, tVals] = odeToVectorField(equations)
vfEquations = 

(1475739525896764129281770887431076117-6Y2-2877692075498690052096Y1885443715538058583010348331692984375Y11152921504606846976-27670116110564328125Y29223372036854775808)

tVals = 

(Iw)

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: [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 コンポーネントの作成

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/Bdoc19a_1054962_126006/tpc4e40eb3/symbolic-ex98670381' ...