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

Скрипт Open Live Script

Создайте пользовательские основанные на уравнении компоненты для Библиотеки 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.

(J ∂∂t w(t)=Ki I(t)−Dr w(t)L ∂∂t I(t)+R I(t)=V(t)−Kb w(t))

Параметризация компонента 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 = 

(DrJKbKiLR0.0000030.0000010.0000720.0000720.0000123.91 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."1 kg"kilogram - a physical unit of mass." m"meter - a physical unit of length."21 V"volt - a physical unit of electric potential." s"second - a physical unit of time."rad"radian - a physical unit of plane angle."1 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.")[Доктор, J, Кбит, Ки, L, R; vpa ('0.000003'), vpa ('0.000001'), vpa ('0.000072'), vpa ('0.000072'), vpa ('0.000012'), vpa ('3.9'); 1* ((symunit ('N')*symunit ('m') *symunit ('))/symunit ('рад')), 1*symunit ('кг') *symunit ('m') ^2, 1* ((symunit ('V') *symunit ('))/symunit ('рад')), 1* ((symunit ('N')*symunit ('m'))/symunit), symunit ('H'), symunit ('Ом')]

Добавьте названия параметра к рабочему пространству 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 = 

(I(t)V(t)Vb(t)τ(t)torque(t)w(t)000000A"ampere - a physical unit of electric current."V"volt - a physical unit of electric potential."V"volt - a physical unit of electric potential."1 N"newton - a physical unit of force." m"meter - a physical unit of length."1 N"newton - a physical unit of force." m"meter - a physical unit of length."1 rad"radian - a physical unit of plane angle."s"second - a physical unit of time.")

Добавьте имена переменных в рабочее пространство 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\it представляет ток и w(t) скорость вращения.

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

J ∂∂t w(t)+torque(t)=τ(t)-Dr w(t)J*diff (w (t), t) + крутящий момент (t) == tau (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 ∂∂t I(t)+R I(t)=V(t)-Vb(t)
eq4 = Vb(t) == Kb*w(t)
eq4 = Vb(t)=Kb w(t)Vb(t) == Kb*w(t)

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

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

(J ∂∂t w(t)+torque(t)=τ(t)-Dr w(t)τ(t)=Ki I(t)L ∂∂t I(t)+R I(t)=V(t)-Vb(t)Vb(t)=Kb w(t))

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

operands = children(eqs);
operList = [ operands{:} ];
lhs = operList(1:2:end)
lhs=1×4 cell array
    {1x1 sym}    {1x1 sym}    {1x1 sym}    {1x1 sym}


rhs = operList(2:2:end)
rhs=1×4 cell array
    {1x1 sym}    {1x1 sym}    {1x1 sym}    {1x1 sym}

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

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

J ∂∂t w(t)+torque(t)=Ki I(t)-Dr w(t)
equations(2) = subs(eqs(3), lhs(4), rhs(4))
equations = 

(J ∂∂t w(t)+torque(t)=Ki I(t)-Dr w(t)L ∂∂t I(t)+R I(t)=V(t)-Kb w(t))
equations.'
ans = 

(J ∂∂t w(t)+torque(t)=Ki I(t)-Dr w(t)L ∂∂t I(t)+R I(t)=V(t)-Kb w(t))

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

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

(0.000001 ∂∂t w(t)=0.000072 I(t)-0.000003 w(t)0.000012 ∂∂t I(t)+3.9 I(t)=1.0-0.000072 w(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.151120734vpa ('6.151120734')
steadyStateI = vpa(result.ii,10)
steadyStateI = 0.2562966973vpa ('0.2562966973')

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

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

Создайте допустимый вход для ode45 от символьных уравнений. Используйте odeToVectorField создать процедуру MATLAB, которая представляет динамическую систему dYdt=f(t,Y) с начальным условием Y(t0)=Y0.

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

(1475739525896764129281770887431076117-6 Y2-2877692075498690052096 Y1885443715538058583010348331692984375 Y11152921504606846976-27670116110564328125 Y29223372036854775808)[sym ('147573952589676412928/1770887431076117') - 6*Y (2) - (sym ('2877692075498690052096') *Y (1))/sym ('8854437155380585'); (sym ('83010348331692984375') *Y (1))/sym ('1152921504606846976') - (sym ('27670116110564328125') *Y (2))/sym ('9223372036854775808')]
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]. Первым значением является текущий 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.151120734vpa ('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')

Figure contains 2 axes. Axes 1 with title DC Motor - angular velocity contains 2 objects of type line. These objects represent w(t), w(t): steady state. Axes 2 with title DC Motor - current contains 2 objects of type line. These objects represent I(t), I(t): steady state.

Создайте компонент 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/BR2021ad_1657350_232851/mlx_to_docbook18/tpffc34045/symbolic-ex98670381' ...

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

Поддержка

Для просмотра документации необходимо авторизоваться на сайте