Настройте и расширьте библиотеки Simscape для пользовательского двигателя постоянного тока

Создайте пользовательские основанные на уравнении компоненты для Библиотеки Simscape с помощью Symbolic Math Toolbox.

Введение

Symbolic Math Toolbox обеспечивает гибкий способ разработать модели из первых технических принципов в любой пространственной размерности. Можно разработать математические модели для устойчивого состояния или переходной физики.

Можно разработать и решить уравнения, требуемые представлять физику, необходимую для компонента; и выполните свою собственную модель уменьшаемого порядка, сопоставляющую между входными параметрами x и количеством интереса f (x).

Здесь f является собственным компонентом, который может представлять управляющие уравнения в форме:

Математические формулы
Числовые симуляции ОДУ и УЧП

Шаги в этом примеры

Параметризация компонента Simscape использование symReadSSCVariables
Определите пользовательские уравнения для своего компонента Simscape использование diff
Решите установившиеся уравнения аналитически с помощью solve и subs
Решите зависящие от времени уравнения численно в MATLAB при помощи matlabFunction и ode45
Создайте компонент Simscape использование symWriteSSC

Чтобы запустить этот пример, у вас должны быть лицензии на Simscape и Symbolic Math Toolbox.

Модель двигателя постоянного тока

Двигатель постоянного тока является устройством для преобразования электроэнергии в механическую энергию и наоборот. Схематический из двигателя постоянного тока показывают ниже (оставленный фигуру). Блоки, которые моделируют двигатели постоянного тока, обеспечиваются в Simscape Electrical™ (правильная фигура), который является дополнительным продуктом для Simscape.

В этом примере мы выведем представление модели уменьшаемого порядка двигателя постоянного тока с помощью управляющих Обыкновенных дифференциальных уравнений (ОДУ). Для двигателя постоянного тока электрическое напряжение и текущий выведено из законов Кирхгоффа, и формула механического крутящего момента выведена из законов Ньютона. Используя эти уравнения мы можем реализовать пользовательский и параметрический компонент Simscape.

$(\begin{matrix} J \frac{\partial}{\partial t} w (t) = Ки I (t) - Доктор w (t) \\ L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Кбит w (t) \end{matrix})$

Параметризация компонента Simscape

Импортируйте параметры и переменные компонента шаблона

Предположим, что у вас есть компонент Simscape MyMotorTemplate.ssc в вашей текущей папке или на пути MATLAB по умолчанию. Этот компонент еще не имеет никаких уравнений. Шаблон записывает параметры и переменные, которые будут использоваться, чтобы разработать наш двигатель. Можно использовать 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{matrix} Доктор & J & Кбит & Ки & L & R \\ 0.000003 & 0.000001 & 0.000072 & 0.000072 & 0.000012 & 3.9 \\ 1 \frac{N "ньютон - физическая единица измерения силы". m "метр - физическая единица измерения длины". s "второй - физическая единица измерения времени".}{рад "радиан - физическая единица измерения плоского угла".} & 1 kg "килограмм - физическая единица измерения массы". {m "метр - физическая единица измерения длины".}^{2} & 1 \frac{V "вольт - физическая единица измерения электрического потенциала". s "второй - физическая единица измерения времени".}{рад "радиан - физическая единица измерения плоского угла".} & 1 \frac{N "ньютон - физическая единица измерения силы". m "метр - физическая единица измерения длины".}{A "ампер - физическая единица измерения электрического тока".} & H "henry - физическая единица измерения индуктивности". & Ω "Ом - физическая единица измерения сопротивления". \end{matrix})$

Добавьте названия параметра к рабочему пространству MATLAB при помощи функции syms. Параметры появляются как символьные переменные в рабочей области. Можно использовать 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{matrix} I (t) & V (t) & Vb (t) & τ (t) & крутящий момент (t) & w (t) \\ 0 & 0 & 0 & 0 & 0 & 0 \\ A "ампер - физическая единица измерения электрического тока". & V "вольт - физическая единица измерения электрического потенциала". & V "вольт - физическая единица измерения электрического потенциала". & 1 N "ньютон - физическая единица измерения силы". m "метр - физическая единица измерения длины". & 1 N "ньютон - физическая единица измерения силы". m "метр - физическая единица измерения длины". & 1 \frac{рад "радиан - физическая единица измерения плоского угла".}{s "второй - физическая единица измерения времени".} \end{matrix})$

Добавьте имена переменных в рабочее пространство MATLAB при помощи функции syms. Переменные появляются как символьные функции в рабочей области. Проверьте, что вы объявили все необходимые символьные переменные и функции с помощью syms.

syms(varFuns)
syms

Your symbolic variables are:

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

Определите пользовательские уравнения для своего Компонента Simscape

Задайте систему уравнений для моделирования двигателя постоянного тока

Дифференциальные уравнения для механического крутящего момента определены как eq1 и eq2. I(t) представляет текущий и w(t) угловой velocity.

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

eq1(t) = 
 $J \frac{\partial}{\partial t} w (t) + крутящий момент (t) = τ (t) - Доктор w (t)$

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

eq2 =  $τ (t) = Ки 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) = Кбит w (t)$

Мы можем перечислить их вместе. Здесь, крутящий момент двигателя принят, чтобы быть пропорциональным току.

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

eqs = 
 $(\begin{matrix} J \frac{\partial}{\partial t} w (t) + крутящий момент (t) = τ (t) - Доктор w (t) \\ τ (t) = Ки I (t) \\ L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Vb (t) \\ Vb (t) = Кбит w (t) \end{matrix})$

Извлеките левые и правые стороны уравнений.

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

lhs = 
 $(\begin{matrix} J \frac{\partial}{\partial t} w (t) + крутящий момент (t) & τ (t) & L \frac{\partial}{\partial t} I (t) + R I (t) & Vb (t) \end{matrix})$

rhs = operList(2:2:end)

rhs =  $(\begin{matrix} τ (t) - Доктор w (t) & Ки I (t) & V (t) - Vb (t) & Кбит w (t) \end{matrix})$

Вторые и четвертые уравнения задают значения tau(t) и Vb(t). Чтобы упростить систему четырех уравнений к системе двух уравнений, замените этими значениями в первое и третье уравнение.

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

equations = 
 $J \frac{\partial}{\partial t} w (t) + крутящий момент (t) = Ки I (t) - Доктор w (t)$

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

equations = 
 $(\begin{matrix} J \frac{\partial}{\partial t} w (t) + крутящий момент (t) = Ки I (t) - Доктор w (t) & L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Кбит w (t) \end{matrix})$

equations.'

ans = 
 $(\begin{matrix} J \frac{\partial}{\partial t} w (t) + крутящий момент (t) = Ки I (t) - Доктор w (t) \\ L \frac{\partial}{\partial t} I (t) + R I (t) = V (t) - Кбит w (t) \end{matrix})$

Прежде, чем решить уравнения, замените параметрами с их числовыми значениями. Кроме того, используйте V(t) = 1.

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

ans = 
 $(\begin{matrix} 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{matrix})$

Решите установившиеся уравнения аналитически

Решите уравнения для устойчивого состояния.

Для этого удалите зависимость времени для функций 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$

Решите зависящие от времени уравнения численно

Численно моделируйте символьное выражение в MATLAB при помощи `matlabFunction` и `ode45`.

Создайте допустимый вход для ode45 от символьных уравнений. Используйте odeToVectorField, чтобы создать процедуру MATLAB, которая представляет динамическую систему $\frac{dY}{dt} = f (t, Y)$ с начальным условием $Y (t_{0}) = Y_{0}$ .

[vfEquations, tVals] = odeToVectorField(equations)

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

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

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]. Первым значением является текущий I(t), и второе значение является угловой скоростью 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.

Сгенерируйте код Simscape с помощью исходных уравнений eqs.

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/BR2019ad_1062519_57051/mlx_to_docbook10/tp0bc82050/symbolic-ex98670381' ...

Документация