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

Открыть Live Script

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

Введение

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

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

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

  • Математические формулы

  • Численные симуляции ОДУ и PDE

Шаги в этих примерах:

  • Параметризация компонента Simscape с помощью symReadSSCVariables

  • Задайте пользовательские уравнения для компонента Simscape с помощью diff

  • Решите статические уравнения аналитически, используя solve и subs

  • Решите зависящие от времени уравнения численно в MATLAB при помощи matlabFunction и ode45

  • Создайте компонент Simscape с помощью symWriteSSC

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

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

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

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

(J ∂∂t w(t)=Ki Я(t)−Dr w(t)L ∂∂t Я(t)+R Я(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.")[Dr, J, Kb, Ki, 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 ('s ') )/symunit (' rad ')), 1 * symunit (' kg ') * symunit (' m ') ^ 2, 1 * (symunit (' V ') * symunit (' s ') )/symunit (' rad '

Добавьте имена параметров в Рабочее пространство 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. I(t) представляет ток и 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 и я (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

Поддержка

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