Моделирование управления AC двигателя — Field oriented control of PMSM с помощью opensource решений / forpes.ru

Главная
Моделирование управления AC двигателя — Field oriented control of PMSM с помощью opensource решений

Моделирование управления AC двигателя — Field oriented control of PMSM с помощью opensource решений +25

14.02.2025 14:47

vladipirogov 4 1100 Источник

В этой статье я хочу поделиться результатом своих исследований в области моделирования систем управления двигателями переменного тока. В качестве объекта управления был выбран синхронный двигатель с постоянными магнитами PMSM (Permanent Magnet Synchronous Machine) как наиболее распространенная машина в современных транспортных средствах. Основное внимание будет уделено построению математической модели системы, объекта управления, и алгоритмов для симуляции. Для реализации модели я выбрал open source решения: Python control, Scilab. Мне было интересно, возможно ли использование свободных средств моделирования для построения более-менее сложных и реальных систем. Далее я поделюсь своими впечатлениями. В первой части статьи приводится теоретический материал, где описываются основные уравнения двигателя и элементы теории управления. Для теоретической части необходимы базовые понимания электротехники, ниже приложу ссылки, где можно обновить знания. Я постарался проработать разные источники литературы, чтобы взять необходимый минимум, с которым самому пришлось столкнуться для понимания сути процессов управления двигателем. Читатель вправе пропустить матчасть и перейти сразу к описанию реализации, и при необходимости вернуться к некоторым теоретическим аспектам в этом материале, или других источниках. Реализация алгоритмов управления построена по классическому принципу с помощью диаграммы потоков.

Содержание

Векторное управление: концепции и принципы
Источники для быстрого старта
Обобщенная электрическая машина
Пространственно-векторное представление
Преобразование Кларка - ABC to αβ (Clarke Transformation)
Преобразование Парка - αβ to dq (Park Transformation)
Синхронная машина с постоянными магнитами - PMSM
Уравнения PMSM в координатах dq
Основы теории управления
Диаграмма алгоритма управления FOC
Scilab
Python Control
Реализация модели системы управления PMSM
Параметры модели
Симуляция системы управления
- Результаты симуляции Scilab / Xcos
- Симуляция Python control
Заключение
Используемые источники

Векторное управление: концепции и принципы

Электрические машинами переменного тока: асинхронные двигатели (IM), синхронные двигатели с постоянными магнитами (PMSM) являются основными в современных электрических приводах.
Существует два основных подхода к управлению машинами переменного тока: скалярное и векторное. Скалярное управление проще в реализации и обеспечивает стабильные рабочие характеристики, но его динамика ограничена, поскольку оно не управляет переходными процессами. Для достижения высокой точности и быстрого отклика были разработаны векторные методы управления, использующие замкнутые контуры обратной связи.

Использование регулируемых приводов позволяет значительно снизить энергопотребление и обеспечить быстрый и точный отклик в промышленных приложениях. В 1970-х годах были предложены принципы управления моментом и потоком, получившие название Field-Oriented Control, FOC или векторное управление. Изначально они применялись для асинхронных машин с короткозамкнутым ротором, а позже — для синхронных машин. Основная идея векторного управления заключается в управлении векторами тока статора аналогично тому, как это делается в двигателях постоянного тока, но с более сложными алгоритмами.

Метод векторного управления был впервые предложен Blaschke в 1970-х годах как способ решения этой проблемы. Основной принцип управления по ориентации поля (FOC) позволил добиться независимого управления моментом и потоком.

Ключевая идея метода заключается в представлении переменных в произвольной системе координат. Если система координат вращается вместе с вектором магнитного потока, то метод управления называется поток-ориентированным (Field-Oriented Control), ротор-ориентированным и т. д. Это позволяет применять линейные методы управления, аналогичные подходам для двигателей постоянного тока, но с учетом особенностей машин переменного тока.

Источники для быстрого старта

Действия с матрицами - Данное методическое пособие поможет Вам научиться выполнять действия с матрицами.
Лекции и уроки по высшей математике
PMSM motor vector control за две минуты
Электротехник Равилов - Канал посвящен задачам и теориям по электротехнике
Field Oriented Control of Permanent Magnet Motors - объяснения FOC от Dave Wilson Texas Instruments
Motor Control - онлайн-курс от STMicroelectronics
Репозиторий с примерами статьи

Обобщенная электрическая машина

В конце 1950-х годов американские исследователи Дэвид Уайт и Херберт Вудсон предложили метод представления любых электромеханических преобразователей в виде универсальной электромагнитной машины. В основе этой модели лежит концепция с двухфазным явнополюсным статором и неявнополюсным ротором. Они доказали, что любую n-фазную электрическую машину можно эквивалентно преобразовать в двухфазную систему с помощью фазовых преобразований.

На момент создания данной модели приводы постоянного тока преобладали в промышленности, поэтому конструкция предполагала наличие явнополюсного статора, а ротор, напротив, был гладким. В случае явнополюсных синхронных машин, требующих четко выраженных полюсов на роторе, ученые предложили использовать обратную модель машины.

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

На рисунке представлена схема машины. Оси α и β жестко связаны со статором, на котором размещены две ортогональные обмотки. Эти обмотки показаны в виде концентрированных катушек. Входящий в плоскость рисунка ток обозначен крестиком в кружке, а выходящий — точкой. Катушка фазы α расположена так, что формирует магнитное поле вдоль оси α. На рисунке система изображена в виде пространственных электрических цепей. Индуктивность обмотки статора по оси α также генерирует магнитное поле вдоль этой оси. Для наглядности активные сопротивления обмоток в схеме не показаны, однако их влияние учитывается в уравнениях:

$u_{sα} = R_s i_{sα} + \frac{dΨ_{sα}}{dt}, \quad u_{sβ} = R_s i_{sβ} + \frac{dΨ_{sβ}}{dt}$

где $u_{sα}$ , $u_{sβ}$ — напряжения фаз статора, $i_{sα}$ , $i_{sβ}$ — токи фаз статора, R_s — сопротивление обмоток (считается одинаковым для обеих фаз), $Ψ_{sα}$ , $Ψ_{sβ}$ — потокосцепления фаз статора.

Ротор может свободно вращаться, а его обмотки ориентированы вдоль осей d и q, которые привязаны к положению ротора. Уравнения для ротора принимают следующий вид:

$u_{rd} = R_r i_{rd} + \frac{dΨ_{rd}}{dt}, \ {rq} = R_r i_{rq} + \frac{dΨ_{rq}}{dt}$

где $u_{rd}$ , $u_{rq}$ — напряжения ротора, $i_{rd}$ , $i_{rq}$ — токи обмоток ротора, R_r — сопротивление (одинаковое для обеих фаз), $Ψ_{rd}$ , $Ψ_{rq}$ — потокосцепления фаз ротора.

Таким образом, модель Уайта и Вудсона позволяет универсально описывать поведение различных электромагнитных машин, сводя их к двухфазному представлению, что упрощает анализ и расчеты.

Пространственно-векторное представление

Трехфазная машина переменного тока может быть описана с использованием метода пространственных векторов. Переменные двигателя KA(t), KB(t), KC(t) для симметричной машины удовлетворяют условию:

Суммирование этих переменных дает пространственный вектор:

где $a = e^{j2/3\pi}$ и $a^2 = e^{j4/3\pi}$
На рисунке ниже показан комплексный пространственный вектор тока статора.

Пространственный вектор K может представлять переменные двигателя (например, ток, напряжение и поток). Принцип векторного управления машин переменного тока использует преобразования переменных из физической трехфазной системы ABC в неподвижную систему отсчета αβ или вращающуюся систему отсчета dq. Первое преобразование может быть выполнено в неподвижную систему отсчета αβ - преобразование Кларка (Clark transformation), а затем во вращающуюся систему отсчета dq - преобразование Парка (Park transformation).

Преобразование Кларка - ABC to αβ (Clarke Transformation)

Синусоидальные трехфазные переменные можно представить в виде пространственного вектора, выраженного на двух ортогональных осях (αβ), как показано на рисунке.

Вектор тока статора описывается следующем комплексным выражением:

$\bar{i_s} = i_{s\alpha} + ji_{s\beta}$

где:

$\begin{aligned}i_{s\alpha} = Re\{2/3[i_{sA} + \alpha i_{sB} + \alpha^2i_{sC}]\} \\i_{s\beta} = Im\{2/3[i_{sA} + \alpha i_{sB} + \alpha^2i_{sC}]\}\end{aligned}$

Тогда

$\begin{aligned}i_{s\alpha} &= i_{sA} \\i_{s\beta} &= \frac{1}{\sqrt{3}} (i_{sA} + 2i_{sB})\end{aligned}$

Или эквивалентно в матричной форме:

$\begin{bmatrix}i_{s\alpha} \\i_{s\beta}\end{bmatrix} =\begin{bmatrix}1 & 0 & 0 \\\frac{1}{\sqrt{3}} & \frac{2}{\sqrt{3}} & 0 \end{bmatrix}\begin{bmatrix}i_{sA} \\ i_{sB} \\ i_{sC}\end{bmatrix}$

Реализация на Python для токов:

    def clarc_transform(i_a, i_b):
        i_y = (1/math.sqrt(3))*(i_a + 2*i_b)
        i_x = i_a
        return i_x, i_y

Реализация на Python для напряжений

    def clark_transform(u_a, u_b):
        u_y = 1/math.sqrt(3)*u_a + 2/math.sqrt(3)*u_b
        u_x = u_a
        return u_x, u_y

Обратное преобразование Кларка из αβ в ABC:

$\begin{bmatrix}i_{sA} \\i_{sB} \\i_{sC}\end{bmatrix} =\begin{bmatrix}1 & 0 \\-\frac{1}{2} & -\frac{\sqrt{3}}{2}\\-\frac{1}{2} & \frac{\sqrt{3}}{2}\end{bmatrix}\begin{bmatrix}i_{s\alpha} \\ i_{s\beta}\end{bmatrix}$ Реализация на Python:```python def inverse_clark(i_x, i_y): i_b = -1/2*i_x + math.sqrt(3)/2*i_y i_c = -1/2*i_x - math.sqrt(3)/2*i_y return i_x, i_b, i_c```

Реализация на Python:```python def inverse_clark(i_x, i_y): i_b = -1/2*i_x + math.sqrt(3)/2*i_y i_c = -1/2*i_x - math.sqrt(3)/2*i_y return i_x, i_b, i_c```

Преобразование Парка - αβ to dq (Park Transformation)

Здесь проецируются фазовые компоненты αβ в неподвижной системе отсчета на вращающуюся систему отсчета dq с угловой скоростью ωk. Если мы рассмотрим ось d, совмещенную с потоком ротора по оси d, то мы говорим о системе, ориентированной на поток ротора. В качестве примера, пространственный вектор тока статора и его компонент в (α,β) и в (d,q) показаны на рисунке.

Вектор тока в dq-координатах:

$\bar{i_s} = i_{sd} + ji_{sq}$

Что эквивалентно:

$\bar{i_s} = (i_{s\alpha}cos(\gamma_s) + i_{s\beta}sin(\gamma_s)) + j(i_{s\beta}cos(\gamma_s) - i_{s\alpha} sin(\gamma_s))$

В матричной форме:

$\bar{i_s} =\begin{bmatrix}i_{sd} \\i_{sq}\end{bmatrix} =\begin{bmatrix}cos(\gamma_s) & sin(\gamma_s) \\-sin(\gamma_s) & cos(\gamma_s)\end{bmatrix}\begin{bmatrix}i_{s\alpha} \\ i_{s\beta}\end{bmatrix}$

Реализация на Python для токов:

    def park_transform(i_sa, i_sb, gamma):
        i_d = i_sa*math.cos(gamma) + i_sb*math.sin(gamma)
        i_q = -i_sa*math.sin(gamma) + i_sb*math.cos(gamma)
        return i_d, i_q

Обратное преобразование Парка:

$\bar{i_s} =\begin{bmatrix}i_{s\alpha*} \\i_{s\beta*}\end{bmatrix} =\begin{bmatrix}cos(\gamma_s) & -sin(\gamma_s) \\sin(\gamma_s) & cos(\gamma_s)\end{bmatrix}\begin{bmatrix}i_{sd*} \\ i_{sq*}\end{bmatrix}$

Реализация на Python:

    def inverse_park(i_d, i_q, gamma):
        i_sa = i_d*math.cos(gamma) - i_q*math.sin(gamma)
        i_sb = i_d*math.sin(gamma) + i_q*math.cos(gamma)
        return i_sa, i_sb

Преобразование может быть применено ко всем пространственным векторам (например, напряжению статора, потоку ротора и т. д.).

Синхронная машина с постоянными магнитами - PMSM

В англоязычной литературе используется термин PMSM - Permanent Magnet Synchronous Machine. В дальнейшем будет использоваться эта аббревиатура.
На рисунке угол $\theta$ то же самое, что угол $\phi$ в уравнениях.

Электрические уравнения для для трехфазной синхронной машины:

$\begin{aligned}u_A = Ri_A + \frac{d\Psi_A}{dt}, \\u_B = Ri_B + \frac{d\Psi_B}{dt}, \\u_C = Ri_C + \frac{d\Psi_C}{dt},\end{aligned}\tag{1}$ $\begin{aligned}Ψ_A &= L_A \cdot i_A + L_{AB} \cdot i_B + L_{AC} \cdot i_C + Ψ_{0A}(\phi),\\Ψ_B &= L_B \cdot i_B + L_{AB} \cdot i_A + L_{BC} \cdot i_C + Ψ_{0B}(\phi),\\Ψ_C &= L_C \cdot i_C + L_{AC} \cdot i_A + L_{BC} \cdot i_B + Ψ_{0C}(\phi)\end{aligned}\tag{2}$

где L_A , L_B , L_C – собственная индуктивность фаз обмоток статора; $L_{АВ}$ , $L_{AC}$ … – взаимная индуктивность между обмотками статора; $Ψ_{0A}(φ)$ – поток сцепления от постоянного магнита, пронизывающий катушку фазы А, зависящий от угла поворота ротора . Взаимные индуктивности статорных обмоток равны половине от максимальной взаимной индуктивности L_m между катушками фаз А, В, С статора, если бы эти катушки можно было совместить (они расположены под углом 120°), то $L_{AB} = L_mcos(120) = -L_m/2$
Собственная индуктивность фаз обмоток статора складывается из индуктивности L_m и индуктивности рассеяния L_p . Если учесть, что i_a + i_b + i_C = 0 , то

$\begin{aligned}Ψ_A &= (L_p + L_m)i_A - \frac{L_m}{2}i_B - \frac{L_m}{2}i_C + Ψ_{0A}(\phi) =\\&= (L_p - L_m)i_A + \frac{L_m}{2}i_A + Ψ_{0A}(\phi) =\\&= (L_p + \frac{3L_m}{2})i_A + Ψ_{0A}(\phi) =\\&= Li_A + Ψ_{0A}(\phi)\end{aligned}\tag{3}$

Аналогичные вычисления для других фаз; индуктивность фазы $L = L_p + \frac{3L_m}{2}$
Изменение Ψ_0(φ) с углом имеет вид

$\begin{aligned}Ψ_{0A}(\phi) &= Ψ_mcos(\phi), \\Ψ_{0B}(\phi) &= Ψ_mcos(\phi - \frac{2\pi}{3}), \\Ψ_{0C}(\phi) &= Ψ_mcos(\phi + \frac{2\pi}{3})\end{aligned}\tag{4}$

После подстановки (4) в (3) и в (1)

$\begin{aligned}u_A &= Ri_A + L\frac{di_A}{dt} - ω \cdot Ψ_m sin(\phi), \\u_B &= Ri_B + L\frac{di_B}{dt} - ω \cdot Ψ_m sin(\phi - \frac{2\pi}{3}), \\u_C &= Ri_C + L\frac{di_C}{dt} - ω \cdot Ψ_m sin(\phi + \frac{2\pi}{3}),\end{aligned}\tag{5}$

где $\frac{d\phi}{dt} = \omega$ .
Момент на валу определяется как производная от электромагнитной энергии W по углу поворота φ при постоянном магнитном потоке.
Электромагнитная энергия:

$\begin{aligned}W &= \frac{L_A \cdot i^2_A}{2} + \frac{L_B \cdot i^2_B}{2} + \frac{L_C \cdot i^2_C}{2} + L_{AB} \cdot i_A \cdot i_B + L_{AC} \cdot i_C \cdot i_B + L_{BC} \cdot i_B \cdot i_C + \\&+ Ψ_m cos(\phi) \cdot i_A + Ψ_m cos(\phi - \frac{2\pi}{3}) \cdot i_B + Ψ_m cos(\phi + \frac{2\pi}{3}) \cdot i_C,\end{aligned}$

момент на валу:

$M_r = - \frac{∂W}{∂\phi} = Ψ_m \cdot [sin(\phi) \cdot i_A + sin(\phi - \frac{2\pi}{3}) \cdot i_B + sin(\phi + \frac{2\pi}{3}) \cdot i_C] \tag{6}$

Уравнения механической части двигателя:

$J\frac{d\omega}{dt} = M - M_C \tag{7}$ $\frac{d\phi}{dt} = \omega \tag{8}$

Угол $\phi$ - важная часть всех уравнений, он отсчитывается от оси катушки А. Поля совпадут, если токи изменяются по закону косинуса с нулевой начальной фазой $\phi_0$ и частотой ω_0

$\begin{aligned}i_A(t) &= I_m cos(ω_0t + \phi_0), \\i_B(t) &= I_m cos(ω_0t - \frac{2\pi}{3} + \phi_0), \\i_C(t) &= I_m cos(ω_0t + \frac{2\pi}{3} + \phi_0). \\\end{aligned}\tag{9}$

Уравнения PMSM в координатах dq

Чтобы получить уравнения в координатах dq, необходимо произвести преобразования Clarke-Park в системе уравнений (5) из описания PMSM.
После всех алгебраических приведений с учетом числа пар полюсов p уравнения принимают следующий вид:

$\begin{aligned}u_d &= R \cdot i_d + L \frac{di_d}{dt} - p\omega L \cdot i_q \\u_q &= R \cdot i_q + L \frac{di_q}{dt} + p\omega L \cdot i_d + p\omega\Psi_m \\M &= \frac{3}{2} \Psi_m \cdot i_q\end{aligned}$

Схема PMSM для оси dq представлена на рисунке. В PMSM основное магнитное поле создается постоянными магнитами. Эти магниты размещаются на роторе. Результирующий поток постоянен во времени, при условии, что ток статора не оказывает никакого влияния. В действительности ток статора создает собственное магнитное поле, влияющее на исходное, которое называется реакцией якоря.

Потокосцепление вдоль осей d и q

$\begin{aligned}ψ_d &= ψ_{ad} + ψ_f = L_di_d + ψ_f \\ψ_q &= ψ_{aq} = L_qi_q\end{aligned}$

Где $ψ_{ad}$ и $ψ_{aq}$ — это поток рассеяния статора вдоль осей d и q соответственно;
ψ_f - тоже что и ψ_m - потокосцепление постоянного магнита, в различных источниках указываются разные индексы

Искажение магнитного поля зависит от положения и величины тока статора. В режиме холостого хода реакция якоря считается равной нулю из-за незначительного тока статора.

ЭДС, индуцированная в обмотках двигателя, согласно уравнений и приведенной схеме представлена в dq-координатах следующим образом:

$\begin{aligned}e_d &= -L_q\omega_ri_q \\e_q &= L_d\omega_ri_d + \omega_r \Psi_f\end{aligned}$

Электромагнитный момент машины может быть получен путем анализа баланса мощности. Необходимо учитывать электрическую мощность, поставляемую машине от трехфазного источника напряжения, механическую мощность на валу ротора и мощность потерь. Исходя из выражения для электрической мощности, и используя баланс напряжения получаем следующее соотношение:

$\begin{aligned}p_e &= \frac{3}{2}(u_di_d - u_qi_q) \\&= \frac{3}{2}R_s(i^2_d + i^2_q) + \frac{3}{2}(i_d\frac{d\Psi_d}{dt} + i_q\frac{\Psi_q}{dt}) + \frac{3}{2}\omega_m(\Psi_di_d - \Psi_qi_d) \\&= p_{Cu1} + \frac{dW_m}{dt} + p_{em}\end{aligned}$

$p_{Cu1}$ - мощность потерь в меди статора. Второе слагаемое $\frac{dW_m}{dt}$ определяет скорость изменения энергии, накопленной в магнитном цепи.
$p_{em}$ - электромагнитная мощность, она получается путем вычитания потерь из входной электрической мощности. Следовательно $p_{em}$ представляет собой меру преобразования электрической энергии в механическую работу. Электромагнитная мощность равна произведению внутреннего электромагнитного момента $T_{em}$ и скорости ротора $\Omega = \omega_m/p$ . Этот момент является механическим взаимодействием между статором и ротором, вызванным электромагнитными силами, создаваемыми полем связи. Электромагнитный момент рассчитывается по выражению:

$T_{em} = \frac{p_{em}}{\Omega_m} = p \frac{p_{em}}{\omega_m} = \frac{3p}{2}(\Psi_di_d - \Psi_qi_d)$

Вводя соотношения: $\Psi_d = L_si_d + L_mi_r = L_si_d + \Psi_{Rm}$ и $\Psi_q = L_si_q$ в приведенное выше выражение крутящего момента, получаем выражение для момента изотропной машины, в которой магнитная цепь имеет одинаковое магнитное сопротивление во всех направлениях:

$T_{em} = \frac{3p}{2}(\Psi_di_d - \Psi_qi_d) = \frac{3p}{2}L_mi_Ri_q = \frac{3p}{2}(L_mi_R)i_q = \frac{3p}{2}\Psi_{Rm}i_q$

Компонент потока $\Psi_{Rm} = L_mi_R$ представляет собой часть потока ротора, вызванного постоянными магнитами, которая охватывает обмотки статора. Небольшая часть потока постоянных магнитов не достигает сердечника статора и не вносит вклад в процесс электромеханического преобразования.
Магнитная цепь ротора может иметь нецилиндрическую форму, что вносит разницу в магнитные сопротивления вдоль осей d и q. Это приводит к разным собственным индуктивностям L_d и L_q обмоток виртуальных фаз статора. Поток виртуальной фазной обмотки статора по оси d равен $\Psi_d = L_di_d + L_m$ , в то время как поток в виртуальной фазной обмотке статора q равен $\Psi_q = L_qi_q$ .
Выражение крутящего момента содержит дополнительный компонент, пропорциональный разнице самоиндукций по осям d и q. Этот компонент крутящего момента называется реактивным моментом, поскольку он появляется как следствие различий в магнитных сопротивлениях.

На рисунке показана роторная магнитная цепь с постоянными магнитами, встроенными во внутреннюю часть ротора. Магниты вложены вдоль оси d. Дифференциальная проницаемость постоянных магнитов близка к $\mu_0$ , и их присутствие на пути потока по оси d увеличивает магнитное сопротивление $R_{md}$ . При $R_{md} > R_{mq}$ индуктивности фазных обмоток L_d < L_q . Потокосцепления машины равны:

$\begin{aligned}\Psi_d &= L_di_d + L_mi_r \\\Psi_q &= L_qi_q \\\Psi_R &= L_Ri_R + L_mi_d\end{aligned}$

Различия между самоиндукциями L_d и L_q виртуальных фаз статора d и q влияют на выражение для электромагнитного момента:

$T_{em} = \frac{3p}{2}(\Psi_di_d - \Psi_qi_d) = \frac{3p}{2}\Psi_{Rm}i_q + \frac{3p}{2}(L_d - L_q)i_di_q$

Для построения автоматического управления PMSM модель двигателя удобно представить в виде дифференциальных уравнений состояний:

$\begin{aligned}\frac{di_d}{dt} &= -\frac{R_s}{L_d}i_d + \frac{L_q}{L_d}\omega_ri_q + \frac{1}{L_d}u_d \\\frac{di_q}{dt} &= -\frac{R_s}{L_q}i_q - \frac{L_d}{L_q}\omega_ri_d -\frac{1}{L_q}\omega_r\Psi_f + \frac{1}{L_q}u_q \\\frac{d\omega_r}{dt} &= \frac{1}{J}[\Psi_fi_q + (L_d - L_q)i_di_q] \\\frac{d\theta_r}{dt} &= \omega_r\end{aligned}$

где R_s — сопротивление статора; L_d и L_q — индуктивности ротора относительно осей d и q; ω_r — угловая скорость ротора; — момент инерции

Модель PMSM может быть представлена в виде следующей диаграммы потоков

Основы теории управления

Здесь я приведу только некоторые понятия теории управления, которые далее используются в алгоритмах системы. Я упускаю описание анализа системы на устойчивость и частотные характеристики. Для полного погружения в теорию систем управления можно использовать рекомендованную литературу или посмотреть курс лекций по ТАУ, который можно найти на просторах youtube.

Представление систем в пространстве состояний

Для удобства анализа моделей объектов управления их целесообразно приводить к стандартному виду, позволяющему применять существующие методы решения. В теории управления таким стандартом является система дифференциальных уравнений первого порядка, известная как нормальная форма Коши.

Рассмотрим модель электродвигателя при отсутствии внешней нагрузки, то есть M_H(t) = 0 . Учитывая, что угловая скорость ω(t) связана с углом поворота вала θ(t) соотношением:

$\dot{\theta}(t) = \omega(t)$

можно записать уравнение движения двигателя в виде системы:

$\dot{\omega}(t) = -\frac{k_1k_2}{J R} \omega(t) + \frac{k_1}{J R} u(t)$

Запишем эту систему в матричной форме:

$\begin{bmatrix} \dot{\theta}(t) \\ \dot{\omega}(t) \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & -\frac{k_1k_2}{J R} \end{bmatrix} \begin{bmatrix} \theta(t) \\ \omega(t) \end{bmatrix} + \begin{bmatrix} 0 \\ \frac{k_1}{J R} \end{bmatrix} u(t)$

Здесь переменные θ(t) и ω(t) определяют текущее состояние системы. Если их значения известны в момент времени t_0 , а также задан входной сигнал u(t) для всех $t \geq t_0$ , можно однозначно определить поведение системы на будущее, вне зависимости от предшествующих значений. В связи с этим, θ(t) и ω(t) называются переменными состояния, а вектор:

$\mathbf{x}(t) = \begin{bmatrix} \theta(t) \\ \omega(t) \end{bmatrix}$

– вектором состояния.

Приняв обозначения $\mathbf{x}(t)$ для вектора состояния и u(t) для входного сигнала, модель можно записать в стандартной форме:

$\dot{\mathbf{x}}(t) = A \mathbf{x}(t) + B u(t)$

где матрицы:

$A = \begin{bmatrix} 0 & 1 \\ 0 & -\frac{k_1k_2}{J R} \end{bmatrix}, B = \begin{bmatrix} 0 \\ \frac{k_1}{J R} \end{bmatrix}$

Эта модель описывает связь между входом системы u(t) и вектором состояния $\mathbf{x}(t)$ , поэтому она называется моделью вход-состояние.

Для полной характеристики системы вводится уравнение выхода, определяющее зависимость выходного сигнала y(t) от состояния и входа:

$y(t) = C \mathbf{x}(t) + D u(t)$

Совместно с уравнением состояния эта модель записывается в виде:

$\begin{cases} \dot{\mathbf{x}}(t) = A \mathbf{x}(t) + B u(t), \\ y(t) = C \mathbf{x}(t) + D u(t) \end{cases}$

где матрицы C и D определяют, какие переменные состояния и входные сигналы формируют выход системы.

Например, если выходной переменной является угол поворота вала, то:

$C = \begin{bmatrix} 1 & 0 \end{bmatrix}, D = 0$

Если же выходом является угловая скорость двигателя, то:

$C = \begin{bmatrix} 0 & 1 \end{bmatrix}, D = 0$

Благодаря гибкости модели (возможности варьировать матрицы C и D), можно анализировать различные параметры системы, выбирая в качестве выхода любую линейную комбинацию переменных состояния и входного сигнала.

Если параметры системы (например, момент инерции J, сопротивление якоря R, коэффициенты k_1 и k_2 ) остаются неизменными во времени, модель является стационарной. В противном случае она называется нестационарной.

Переходная функция

Один из распространённых способов моделирования динамических систем методом «вход-выход» — анализ реакции системы на стандартные тестовые сигналы. Одним из таких сигналов является так называемый единичный скачок (или ступенчатый сигнал), который представляет собой мгновенное изменение входного сигнала с 0 до 1 в момент времени t = 0:

$x(t) = \begin{cases} 0, & t < 0, \\ 1, & t \geq 0. \end{cases}$

Реакция системы на этот сигнал называется переходной функцией h(t). При этом предполагается, что начальное состояние системы — нулевое, то есть все её внутренние переменные равны нулю, а накопленная энергия отсутствует.

Если начальные условия отличны от нуля, для расчёта выходного сигнала при любом входе требуется учитывать дифференциальные уравнения системы или её модель в пространстве состояний. Таким образом, переходная характеристика содержит меньше информации, чем полное математическое описание системы.

Переходная характеристика системы первого порядка

Рассмотрим систему, описываемую линейным дифференциальным уравнением первого порядка:

$T \frac{dy(t)}{dt} + y(t) = kx(t),$

где:

k — коэффициент усиления (безразмерный),
T — постоянная времени, характеризующая инерционность системы (измеряется в секундах).

Рассчитаем переходную характеристику для данной системы, предполагая, что входной сигнал x(t) представляет собой единичный скачок x(t) = 1 при t > 0. Решение уравнения при начальном условии y(0) = 0 даёт:

$h(t) = k \left(1 - e^{-t/T} \right).$

График переходной функции показывает, что при увеличении T выход y(t) достигает установившегося значения медленнее. Таким образом, постоянная времени T характеризует инерционность системы: чем больше её значение, тем медленнее система реагирует на изменение входного сигнала.

Практически ступенчатый сигнал можно легко воспроизвести в экспериментальных условиях, что делает переходную характеристику удобной для исследования динамических свойств системы.

Импульсная характеристика (весовая функция)

В качестве тестового входного воздействия можно использовать не только ступенчатый сигнал, но и импульсный. Рассмотрим, что произойдёт, если использовать прямоугольный импульс малой ширины, но фиксированной площади.

По мере уменьшения ширины импульса его амплитуда будет расти, и в пределе, когда ширина стремится к нулю, амплитуда станет бесконечной. Это приводит к понятию единичного импульса или дельта-функции Дирака δ(t), которая определяется следующим образом:

$\delta(t) = \begin{cases} 0, & t \neq 0, \\ \infty, & t = 0. \end{cases}$

при этом её интеграл равен единице:

$\int_{-\infty}^{\infty} \delta(t) dt = 1.$

Так как дельта-функция является формально бесконечным импульсом, на графиках её принято изображать в виде стрелки.

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

Реакция системы на единичный импульс δ(t) называется импульсной характеристикой и обозначается w(t). Она связана с переходной функцией следующим соотношением:

$w(t) = \frac{dh(t)}{dt}.$

Наоборот, переходную функцию можно выразить через импульсную характеристику:

$h(t) = \int_{0}^{t} w(\tau) d\tau.$

Для системы первого порядка с переходной функцией $h(t) = k (1 - e^{-t/T})$ , импульсная характеристика равна:

$w(t) = \frac{k}{T} e^{-t/T}.$

Импульсная характеристика также называется весовой функцией, так как выходной сигнал y(t) при любом входном воздействии x(t) можно выразить через свёртку:

$y(t) = \int_{0}^{\infty} w(\tau) x(t - \tau) d\tau.$

Этот интеграл показывает, что весовая функция w(t) определяет вклад каждого значения входного сигнала в формирование выходного отклика системы.

В отличие от ступенчатого сигнала, мгновенный импульс бесконечной величины невозможно реализовать физически, поэтому импульсную характеристику нельзя получить напрямую в эксперименте, а её вычисляют косвенными методами.

Преобразование Лапласа

Одной из задач, поставленных а теории управления - вычисление выхода системы при известном входе. Для этого необходимо решать дифференциальные уравнения. Чтобы упростить процедуру, математики придумали преобразование, которое позволило заменить решение дифференциальных уравнений алгебраическими вычислениями, то есть, операциями с полиномами и рациональными функциями.
Для функции f(t) вводится преобразование Лапласа, которое обозначается как $\mathcal{L} \{ f(t) \}:$

$F(s) = \mathcal{L} \{ f(t) \} = \int^{\infty}_0 f(t) e^{-st} dt$

Функция F(s) - называется изображением для функции оригинала f(t) . Здесь - это
комплексная переменная, которая выбирается так, чтобы интеграл сходился.
Обратное преобразование Лапласа $\mathcal{L}^{-1} \{ f(t) \}$ позволяет вычислить оригинал f(t) по известному изображению F(s) :

$f(t) = \mathcal{L}^{-1} \{ f(t) \} = \frac{1}{2\pi j} \int^{\sigma + j \infty} _{\sigma - j \infty} F(s) e^{st}ds,$

где $j = \sqrt{-1}$ , а постоянная $\sigma$ выбирается так, чтобы интеграл сходился. На практике вместо интеграла чаще всего используют готовые таблицы, по которым можно сразу определить изображение по оригиналу и наоборот. Для оператора Лапласа свойственны линейные преобразования, такие как суперпозиция, т.е. оператор суммы равняется сумме операторов. При нулевых начальных условиях изображение производной равно изображению самой функции, умноженному на .

Передаточная функция

Выходной сигнал системы можно представить как результат действия некоторого оператора на ее вход. Для линейных моделей такой оператор можно записать следующим образом.
Пусть модель объекта задана линейным дифференциальным уравнением второго порядка, связывающим вход x(t) и выход y(t) :

$b_2\frac{d^2y(t)}{dt^2} + b1\frac{dy(t)}{dt} + b_0y(t) = a_1\frac{dx(t)}{dt} + a_0x(t),$

где a_i(i=0, 1) и b_i(i=0, 1, 2) - постоянные.
Применим к левой и правой частям преобразование Лапласа, считая, что все начальные условия нулевые. Получается уравнение в изображениях, связывающее преобразования Лапласа входа X(s) и выхода Y(s) :

$b_2 \cdot s^2Y(s) + b_1 \cdot sY(s) + b_0 \cdot Y(s) = a_1 \cdot sX(s) + a_0 \cdot X(s)$

Можно вынести за скобки Y(s) в левой части и X(s) в правой части:

$(b_2s^2 + b_1s + b_0) \cdot Y(s) = (a_1s + a_0) \cdot X(s).$

Разделив обе части этого равенства на b_2s^2 + b_1s + b_0 , получаем

$Y(s) = \frac{a_1s + a_0}{b_2s^2 + b_1s + b_0} \cdot X(s) = W(s) \cdot X(s),$

где

$W(s) = \frac{a_1s + a_0}{b_2s^2 + b_1s + b_0}$

W(s) - это это передаточная функция объекта, записанная в виде функции от комплексной переменной .
Таким образом, при нулевых начальных условиях изображение выхода линейного объекта вычисляется как произведение его передаточной функции на изображение входного сигнала. Из этого следует и другой важный вывод: передаточная функция равна отношению изображений по Лапласу выхода и входа при нулевых начальных условиях.
Нулями передаточной функции называются корни ее числителя, а полюсами – корни
знаменателя. Например, функция $W(λ) = \frac{λ-1}{λ^2 + 3λ + 2}$ имеет нуль в точке λ = 1 и два полюса в точках λ = -1 и λ = -2.

Основные типы динамических звеньев

Современные системы управления состоят из отдельных функциональных блоков, каждый из которых описывается математическими уравнениями низкого порядка (чаще всего первого или второго). Чтобы эффективно анализировать поведение системы в целом, необходимо понимать, как ведут себя ее компоненты. При анализе характеристик сложной системы передаточную функцию разбивают на простейшие составляющие

Усилительное звено

Компоненты, обладающие ненулевым коэффициентом усиления постоянного сигнала, называются позиционными. Это означает, что их передаточная функция имеет ненулевые свободные члены в числителе и знаменателе:

Здесь коэффициент усиления определяет, во сколько раз входной сигнал усиливается на выходе. В идеальном случае выход мгновенно повторяет вход, но с усилением. Для ступенчатого входного воздействия выход также является ступенчатым, но увеличенным в раз. На синусоидальный входной сигнал усилитель действует без изменения фазы, а амплитудная частотная характеристика остается постоянной

Апериодическое звено

Апериодическое звено – одно из наиболее часто встречающихся в системах управления. Оно описывается дифференциальным уравнением:

$T \frac{dy(t)}{dt} + y(t) = k \cdot x(t)$

Его передаточная функция имеет вид:

$W(s) = \frac{k}{Ts+1},$

где – коэффициент усиления, а – постоянная времени, характеризующая инерционность элемента. Чем больше значение , тем медленнее выход реагирует на изменения входного сигнала.

Интегрирующее звено

Простым примером интегрирующего звена является резервуар с водой, где входной сигнал – расход воды, а выходной – уровень жидкости. Математически оно описывается уравнением:

$\frac{d y(t)}{d t} = k x(t)$

Передаточная функция имеет вид:

$W(s) = \frac{k}{s}$

При ступенчатом входе выходное значение возрастает линейно, что подтверждает накапливающий характер интегрирующего звена.

$h(t) = k \cdot t$

Интеграл от дельта-функции на любом интервале, включающем t=0, равен 1. Поэтому w(t) = k (при $t \geq 0$ )

Дифференцирующее звено

Дифференцирующее звено формирует на выходе производную входного сигнала:

$\frac{d y(t)}{d t} = k \frac{d x(t)}{d t}$

Передаточная функция:

производная единичного ступенчатого сигнала 1(t) в точке t=0 – это
дельта-функция δ(t) . Поэтому переходная и весовая функции дифференцирующего звена

$h(t) = k\delta(t), w(t) = k \frac{d \delta(t)}{dt}$

На практике идеальные дифференцирующие звенья недостижимы из-за бесконечных значений производной, однако приближенные реализации широко используются для предсказания быстрых изменений в системах управления.

Запаздывающее звено

Запаздывание в системе проявляется в виде временного смещения сигнала. Примером может служить нагрев воздуха в трубе. Очевидно, что при изменении температуры воздуха датчик обнаружит это не сразу, а через время t = L/v , где L - длина трубы, а v - скорость потока воздуха. В этом случае говорят, что в системе есть транспортное запаздывание на величину t.
Передаточная функция запаздывающего звена имеет экспоненциальную форму:

$W(s) = e^{-\tau s}$

Запаздывание не изменяет форму сигнала, а лишь смещает его во времени.

Структурные схемы

Для наглядного представления структуры системы управления применяются блок-схемы. В таких схемах каждый элемент изображается в виде прямоугольника с его передаточной функцией, а сигналы обозначаются входными и выходными стрелками. В сложных системах элементы соединяются через сумматоры, развилки и обратные связи.

Для суммирующих элементов используют специальное обозначение – круг, разбитый на сектора. Если сектор залит черным цветом, или в нем стоит знак "минус", поступающий в него сигнал вычитается, а не складывается с другими. Разветвление сигнала обозначается точкой.

На следующем рисунке показана типичная схема системы управления кораблем по курсу. Здесь вход x – заданная величина, выход y – фактическая величина. Сигналы e , u и δ обозначают соответственно ошибку регулирования, сигнал управления и управляющее воздействие привода на объект. Сигнал g – это возмущение (влияние внешних факторов]), а m – шум измерений.

Основными способами соединения звеньев являются.
Последовательное, в котором выходной сигнал предыдущего звена является входным сигналом последующего звена

Передаточная функция последовательного соединения звеньев равна
произведению передаточных функций звеньев, входящих в это соединение

$W(s) = \frac{Xout(s)}{Xin(s)} = W1(s) \cdot W2(s) \cdot ... \cdot Wn(s)$

Параллельное, в котором на вход всех звеньев подается одновременно
один и тот же входной сигнал, а выходные сигналы этих звеньев алгебраически
суммируются, образуя общую выходную величину

Передаточная функция параллельного соединения звеньев равна алгебраической сумме передаточных функций звеньев, входящих в это соединение

$W(s) = \frac{Xout(s)}{Xin(s)} = W1(s) + W2(s) + ... + Wn(s)$

Системы с обратной связью

Динамическая система — это система, чье поведение изменяется во времени, зачастую в ответ на внешние воздействия. Обратная связь — это механизм, при котором две (или более) динамические системы связаны таким образом, что они оказывают влияние друг на друга, создавая тесную взаимосвязь в их динамике.

Открытые и замкнутые системы

На практике различают замкнутые (closed-loop) и разомкнутые (open-loop) системы. Если компоненты системы соединены в кольцо, то она считается замкнутой. Если же одно из звеньев разорвано, то система становится разомкнутой. В замкнутом контуре разделение на систему 1 и систему 2 носит условный характер — точка отсчета зависит от того, с какого элемента мы начинаем описание.

Положительная и отрицательная обратная связь

Отрицательная обратная связь - система стремится компенсировать отклонения и стабилизировать свое состояние. Однако существует и положительная обратная связь, при которой усиление сигнала вызывает его дальнейшее увеличение. Такое поведение может приводить к неустойчивости, но также используется в системах, требующих быстрого реагирования.

Примеры положительной обратной связи:

Усиление аплодисментов в аудитории — чем громче хлопают одни, тем больше людей присоединяются.
Лавинообразный рост в электронике, например, в триггерах или генераторах.
В биологии — процесс свертывания крови, где один активированный белок запускает цепную реакцию активации других.

Преимущества и недостатки обратной связи

Обратная связь дает мощные инструменты для управления системами:

Устойчивость к внешним воздействиям (например, поддержание температуры тела).
Линеаризация нелинейных процессов (широко применяется в электронике).
Автоматическое регулирование параметров (как в круиз-контроле автомобиля).

Однако у обратной связи есть и недостатки:

Риск неустойчивости (системы могут самопроизвольно колебаться или даже выйти из-под контроля).
Чувствительность к шумам (требуются специальные фильтры и компенсационные механизмы). Для контура с отрицательной обратной связью имеем следующее представление передаточной функций

$W(s) = \frac{W_1(s)}{1 + W_1(s)W_2(s)}$

Для доказательства заметим, что Y(s) = W_1(s) E(s) , а а изображение ошибки равно

Поэтому

Перенося X(s) в левую часть, получаем

$\begin{aligned}Y(s)[1 + W_1(s)W_2(s)] = W_1(s)X(s) \\Y(s) = \frac{W_1(s)}{1 + W_1(s)W_2(s)} X(s)\end{aligned}$

Если обратная связь – положительная (сигналы x и f складываются), в знаменателе будет стоять знак «минус»:

$W(s) = \frac{W_1(s)}{1 - W_1(s)W_2(s)}$

Дискретные системы автоматического управления

Физический мир — это непрерывная аналоговая система, описываемая дифференциальными уравнениями. Компьютер же работает в дискретном времени. Дискретность по уровню можно игнорировать при достаточной разрядности АЦП, но дискретность по времени критична для быстродействующих систем управления.

Дискретизация сигнала

Аналоговый сигнал преобразуется в дискретный с помощью АЦП. При синхронизации с ШИМ и цифровым регулятором выборки определяются значениями сигнала в моменты kT:

При оверсемплинге (увеличении частоты дискретизации) формула усреднения приобретает вид:

$x[k] = \frac{1}{N} \sum_{i=0}^{N-1} x\left(kT + \frac{i}{N}T \right)$

где N — кратность оверсемплинга. Например, при N=8:

$x_8[k] = \frac{1}{8} \sum_{i=0}^{7} x\left(kT + \frac{i}{8}T \right)$

Оверсемплинг улучшает точность, усредняя данные за период ШИМ.

Для восстановления непрерывного сигнала используется цифро-аналоговое преобразование (ЦАП), чаще всего реализуемое как фиксатор нулевого порядка (ZOH). Он удерживает значение y[k] до следующего такта (k+1)T. В электроприводах роль ЦАП выполняет силовой инвертор с ШИМ, создающий сглаженное напряжение на обмотках двигателя.

После получения значения x[k] АЦП, цифровому регулятору требуется время на вычисления. Поэтому управляющее воздействие y[k+1] применяется с задержкой. Чтобы компенсировать этот сдвиг, возможно использование прогнозирования: регулятор вычисляет y[k], используя предсказанное значение x[k+1].

Дискретное описание системы

Дискретные процессы описываются разностными уравнениями:

y[k-1] + ... + b_N y[k-N] = a_0 x[k] + a_1 x[k-1] + ... + a_N x[k-N]

Используя оператор сдвига z (zx[k] = x[k+1]), это уравнение преобразуется:

$(1 + b_1 z^{-1} + ... + b_N z^{-N}) y[k] = (a_0 + a_1 z^{-1} + ... + a_N z^{-N}) x[k]$

Передаточная функция в области z-преобразования:

$W(z) = \frac{a_0 + a_1 z^{-1} + ... + a_N z^{-N}}{1 + b_1 z^{-1} + ... + b_N z^{-N}}$

Умножая на z^N , избавляемся от отрицательных степеней:

$W(z) = \frac{a_0 z^N + a_1 z^{N-1} + ... + a_N}{z^N + b_1 z^{N-1} + ... + b_N}$

Это позволяет формировать дискретную модель непрерывного объекта.

Процесс дискретизации включает два этапа:

Квантование – деление сигнала по уровню или по времени:
- При квантовании по уровню фиксируются значения, кратные шагу квантования .
- При квантовании по времени значения принимаются в моменты, кратные (периоду дискретизации ).
- Возможно комбинированное квантование.
Модуляция – преобразование квантованного сигнала в дискретную последовательность (АИМ - амплитудно-импульсная модуляция, ШИМ - широтно-импульсная модуляция и т.п.)

На этапе преобразования непрерывного сигнала в дискретный возможно использовать общую математическую форму в виде решетчатой функции, представляемая через $\delta$ -функции:

$f[t] = \sum{f(iT) \delta(t - iT)},$

где где T — период дискретизации. $\delta$ -функция - идеальный сигнал который равен нулю во всех точках, кроме t=0, где он уходит к бесконечность, причем его площадь (интеграл по всей оси времени) равен единице (см ранее).

Для анализа дискретных систем используются:

Z-преобразование – аналог дискретного преобразования Лапласа.
Дискретное преобразование Лапласа

Выбор интервала дискретизации опирается на теорему Котельникова-Шеннона, суть которой заключается в следующем. Пусть спектр сигнала x(t) ограничен максимальной частотой $ω_{max}$ , тогда точное восстановление функции x(t) воз-
можно при частоте квантования , более чем в 2 раза превышающей максимальную частоту $ω_{max}$ , т. е.

$\omega \geq 2\omega_{max}, T < \frac{\pi}{\omega_{max}}$

Определением дискретного преобразования Лапласа можно считать обычное преобразование Лапласа от решетчатой функции, то есть:

$\begin{aligned}F(s) &= \int^{\infty}_0 \sum^{\infty}_{i = 0}f(iT)\delta(t-iT)e^{-st}dt \\&= \sum^{\infty}_{i=0}[f(iT) \int^{\infty}_0 \delta(t - iT)e^{-st}dt] = \\&= \sum^{\infty}_{i=0}f(iT)e^{-st}, \\z &= e^{sT}, \\F(z) &= \sum^{\infty}_{i=0}f(iT)z^{-1}\end{aligned}$

Также существуют таблицы переходов между функциями-оригиналами, их преобразованиями Лапласа и Z-преобразованиями
С помощью программных пакетов как Matlab, Scilab, Python control и т.п. можно выполнить z-преобразования. Вот пример скрипта в Scilab

s = poly(0, 's');  // комплексная величина
W_s = syslin('c', 1/s);  // Интегральное звено в непрерывном времени

T = 0.001;  // Период дискретизации
//Поскольку dscr принимает данные в формате пространства состояний, 
передаточные функции должны быть сначала преобразованы 
в описания пространства состояний.
W_ss = tf2ss(W_s);
W_zss = dscr(W_ss, T);  // Дискретизация интегрального звена
//и преобразовать обратно в передаточные функции
W_z = ss2tf(W_zss)
disp(W_z, "Передаточная функция в Z-пространстве:");

0.001  
-----  
-1 +z  
"Передаточная функция в Z-пространстве:"

Диаграмма алгоритма управления FOC

PMSM — управляемый объект (plant) - синхронный двигатель
F - программная фильтрация сигналов фазных токов
P - датчик положения (инкрементальный, синусно-косинусный и др.)
PI(w) - регулятор угловой скорости
PI(Ψp) - регулятор потока
PI(Id), PI(Iq) - регуляторы продольной и квадратурной составляющих тока
e^+jα, e^-jα - преобразование Clarke
3/2, 2/3 - преобразование Park
SVPWM - (Space Vector Pulse Width Modulation) пространственно-векторной ШИМ
Gate driver - модуль управления затворами силовых ключей инвертора
PWM Inverter - трехфазный силовой инвертер

Рассмотрим основные принципы алгоритма управления FOC на примере диаграммы управления потока.

В нижней части диаграммы находится двигатель PMSM, который выступает в качестве объекта управления. Управление двигателем осуществляется трехфазным инвертором на основе ШИМ (PWM). Программное обеспечение, будет регулировать скорость или положение ротора в соответствии с заданными параметрами.

Измерение токов и преобразования координат

Первый этап алгоритма FOC — измерение фазных токов двигателя с помощью шунтовых резисторов. Для подавления помех используются фильтры нижних частот, обозначенные на диаграмме буквой F. Можно измерять либо все три фазных тока, либо только два и вычислять третий.

Затем выполняется преобразование Clarke (Clarke transformation), переводящее измеренные трехфазные токи в двухфазную систему координат (XY). Это позволяет представить переменные токи в удобной для анализа форме.

Следующий шаг — преобразование Park (Park transformation), где, используя электрический угол потока ротора $\alpha$ , мы определяем продольную (Id) и квадратурную (Iq) составляющие тока. Эти компоненты являются абстрактным представлением трехфазной системы, облегчающим работу регуляторов, так как переменные тока в DQ-координатах становятся аналогами DC-компонент (постоянного тока).

Определение скорости и регулирование тока

Скорость вращения ротора вычисляется как производная его углового положения, измеряемого с помощью энкодера. В диаграмме этот процесс представлен блоком d/dt. Полученный сигнал скорости используется в контуре регулирования скорости для определения требуемого момента (в ньютонах-метрах). На практике для измерения угловой скорости вращения используют ABI-инкрементальный энкодер или синусно-косинусный абсолютный энкодер. В нашей модели мы будем программно вычислять скорость в блоке PMSM.

Зная параметры двигателя, можно рассчитать требуемую квадратурную составляющую тока Iq, поскольку именно она отвечает за генерацию момента. Для рассматриваемой модели Id принимается равным нулю, однако возможные эффекты ненулевой продольной составляющей будут рассмотрены позже.

PI-регуляторы для Id и Iq вычисляют соответствующие напряжения, которые затем используются для управления инвертором. В системе электроприводов напряжение является основным управляющим воздействием, и задача FOC состоит в том, чтобы корректировать его так, чтобы достичь заданных значений момента, скорости или положения.

Инверсия преобразований и ШИМ-модуляция

После расчета напряжений в системе DQ выполняется обратное преобразование Park-Clarke, чтобы получить опорные фазные напряжения Va, Vb, Vc для двигателя.
Оптимальное распределение напряжений в трехфазной системе обеспечивается методом пространственно-векторной ШИМ - SVPWM (Space Vector Pulse Width Modulation), который эффективно использует доступное напряжение шины постоянного тока.

Аппаратное обеспечение

В современных системах управления электродвигателями используются три основных типа силовых транзисторов:

MosFET (Metal-Oxide-Semiconductor Field-Effect Transistor) — отличается высокой скоростью переключения (до сотен килогерц) и низкими потерями при работе на малых и средних мощностях. Идеален для применения в высокочастотных инверторах.
IGBT (Insulated Gate Bipolar Transistor) — биполярный транзистор с изолированным затвором, который сочетает преимущества высоковольтных биполярных транзисторов и управляемости полевых транзисторов. Работает в диапазоне частот от единиц до десятков килогерц, оптимален для мощных промышленных приводов.
GTO (Gate Turn-Off Thyristor) — запираемый тиристор, предназначенный для управления большими токами и напряжениями. Ограничен по частоте коммутаций (~500 Гц), но применяется в мощных преобразовательных системах.

Для управления электродвигателями транзисторы чаще всего используются в составе модулей силовой электроники, которые включают в себя:

IPM (Intelligent Power Modules) — интеллектуальные силовые модули с интегрированными драйверами и защитой.
ASIPM (Application Specific Intelligent Power Modules) — специализированные модули, дополнительно оснащенные датчиками тока, термоконтролем и выпрямительными блоками.

Транзисторные стойки в силовых схемах

Ключевым элементом силовых схем преобразователей частоты является транзисторная стойка — конфигурация, состоящая из двух последовательно соединенных транзисторов и параллельно включенных защитных диодов. Такие стойки формируют основу для построения инверторов, выпрямителей и других силовых цепей управления электродвигателями.

Рисунок - варианты исполнения стойки

Существуют различные конфигурации транзисторных модулей:

Одиночный транзистор с защитным диодом — используется в мощных драйверах.
Стойка из двух транзисторов — применяется в однофазных инверторах.
Модуль из 6-7 транзисторов — формирует трехфазный инвертор для управления асинхронными и синхронными двигателями.

Рисунок - стойка в составе преобразователя с питанием от сети однофазного переменного тока

Трехфазный инвертор — это силовой преобразователь, предназначенный для формирования переменного напряжения из постоянного. Инверторы работают по принципу широтно-импульсной модуляции ШИМ (PWM), при которой выходное напряжение переключается между напряжением шины DC-link и «землей» с высокой частотой. Среднее значение выходного напряжения регулируется за счет изменения скважности ШИМ-сигнала, то есть соотношения времени включенного и выключенного состояний силовых ключей.

Современные инверторы поддерживают управление в четырех квадрантах, что позволяет электродвигателю работать в обоих направлениях как в моторном, так и в генераторном режимах. Это достигается за счет возможности формирования как положительных, так и отрицательных фазных напряжений.

Модель PI-регулятора

Математически PI-регулятор описывается следующим уравнением:

$u(t) = K_p e(t) + K_i\int^t_0 e(t) dt$

где:

- выходной сигнал регулятора (опорный момент или ток в квадратурной оси);
e(t) - ошибка сигнала;
- коэффициент пропорционального усиления;
коэффициент интегрального усиления.

Выбор физических единиц входных и выходных сигналов критически важен, так как коэффициенты усиления K_p и K_i зависят от используемой системы измерений.

Управление скоростью вращения двигателя PMSM

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

Сам регулятор скорости реализован в виде пропорционально-интегрального (PI) контроллера. Это означает, что ошибка скорости умножается на коэффициент пропорционального усиления и одновременно интегрируется. В результате формируется управляющий момент, позволяющий корректировать фактическую скорость ротора в соответствии с заданным значением.

Ниже представлена блок-схема регулятора скорости.

Его выходом является опорное значение тока в квадратурном направлении.
Составляющие регулятора-скорости

Входная физическая величина регулятора: ??[???] || ??[deg/s] || ??[???/?]
Выходная физическая величина регулятора: ???[?] || ???[??]
Значение обратной связи регулятора: ?[???] || ??[deg/s] || ??[???]
PI-регулятор: $u(t) = K_p e(t) + K_i\int^t_0 e(t) dt$

Логика формирования опорного тока

В большинстве случаев входом системы управления двигателем является опорное значение момента, рассчитанное на основе логики управления механической системой. Исходя из этого момента, необходимо вычислить соответствующее значение квадратурного тока, которое сможет обеспечить требуемый электромеханический момент.

При управлении двигателем PMSM продольный и квадратурный индуктивные параметры равны, что упрощает уравнение момента, позволяя учитывать только синхронный момент при расчетах опорного тока.

Регулятор скорости определяет требуемый момент, который затем преобразуется в опорное значение квадратурного тока. Регулятор квадратурного тока поддерживает необходимый уровень этого тока, чтобы обеспечить заданное значение момента. В то же время регулятор продольного тока корректирует его значение для поддержания требуемого магнитного потока в двигателе.
Принцип регулирования:

В случае опорного крутящего момента me[Nm] на выходе нужно указать соответствующие значения ???[?]
Опорный квадратурный ток ???[?] зависит от фактической конструкции PMSM
В случае SPMSMs:
$m_e = \frac{3}{2}p\Psi_pi_q$
В случае IPMSMs:

$m_e = \frac{3}{2}p(\Psi_pi_q + (L_d - L_q)i_di_q) \rightarrow i_q = \frac{2m_e}{3p(\Psi_p + (L-_d - l_q)i_d)}$

Регуляторы тока

В алгоритме FOC используются два регулятора тока:

Регулятор продольного тока (D-ось) — регулирует значение продольного тока, отвечающего за управление магнитным потоком.
Регулятор квадратурного тока (Q-ось) — управляет величиной крутящего момента.

Оба регулятора генерируют выходные управляющие напряжения, которые затем преобразуются в трехфазную систему координат и подаются на инвертор для создания необходимых средних значений напряжений с помощью ШИМ.

Квадратурный ток отвечает за генерацию эффективного момента, а продольный ток регулирует уровень магнитного потока.

Принцип регулирования:

Регулятор продольного тока ()
- Входная физическая величина: $i_{da} \xrightarrow{yields} \Psi_p[Wb]$
- Выходная физическая величина: $u_{vd}[V]$
- Обратная связь c сигналом:
Регулятор квадратурного тока ():
- Входная физическая величина: $i_{qa} \xrightarrow{yields} m_e[Nb]$
- Выходная физическая величина: $u_{vq}[V]$
- Обратная связь c сигналом:

Ограничения и насыщение сигналов

Для корректной работы регуляторов необходимо учитывать ограничения выходных сигналов. Например, управляющие напряжения ограничены доступным напряжением источника питания (аккумулятора или сети). Аналогичным образом, регулятор скорости не должен запрашивать значения выше максимально допустимой скорости вращения двигателя PMSM.

Векторная широтно-импульсная модуляция (SVPWM) в системе управления FOC

Основные принципы работы PWM

Синусоидальная широтно-импульсная модуляция (ШИМ) широко применяется в приводах малой мощности.
Однако в современных промышленных преобразователях частоты используются более совершенные алгоритмы управления шестиплечевым инвертором, среди которых особое место занимает пространственно-векторная ШИМ (SVPWM).

Основное преимущество SVPWM перед классической синусоидальной ШИМ заключается в увеличении амплитуды выходного линейного напряжения, что позволяет более эффективно использовать питающее напряжение. Вторая важная особенность – снижение количества коммутаций транзисторов, что уменьшает потери на переключение и повышает КПД системы.

При использовании синусоидальной ШИМ фазное напряжение представляет собой синусоидальный сигнал с постоянной составляющей 0.5E, где E – напряжение питания инвертора. В линейном напряжении эта постоянная составляющая исчезает, а максимальная амплитуда составляет 0.866E. Чтобы добиться увеличения амплитуды без искажения синусоидальной формы, в управляющий сигнал вводят третью гармонику с коэффициентом порядка 0.15.
Линейное импульсное и усредненное напряжение АB в относительных единицах представлена на рисунке ниже

На входе алгоритма SVPWM подаются управляющие фазные напряжения, а на выходе — оптимизированные фазные напряжения, учитывающие доступное напряжение шины постоянного тока (DC-link). Данный этап преобразует вычисленные управляющие напряжения в сигналы для непосредственного управления силовыми ключами (IGBT или MOSFET) в трехфазном инверторе напряжения.

Описание алгоритма SVPWM

Принцип работы SVPWM основан на представлении вращающегося магнитного поля, создаваемого тремя симметрично расположенными катушками. Если через эти катушки пропускать синусоидальные токи со сдвигом по фазе на 120°, формируется вращающееся магнитное поле. Вектор напряжения в данном случае удобно анализировать с помощью шести секторов по 60°.

В каждый момент времени в шестиплечевом инверторе включены три транзистора. Коммутация транзисторов происходит так, чтобы результирующее линейное напряжение максимально приближалось к синусоидальному. Верхние ключи обозначаются как 1 (включен) или 0 (выключен), а нижние работают в противофазе. Например, в секторе 0–60° включены ключи 1, 4, 6, в секторе 60–120° – 1, 3, 6 и так далее.

Алгоритм SVPWM предусматривает вычисление длительности включения транзисторов на основе проекций вектора напряжения на базисные векторы. Если обозначить длительности включенного состояния как γ1 и γ2, а длительность выключенного состояния как γ0, то они удовлетворяют выражению: γ0 = 1 - γ1 - γ2.

Это время делится между состояниями 111 и 000 пополам. Так как транзистор
Т1 должен обеспечить работу и в состоянии 100, и в состоянии 110, то он будет включен по времени

$γ_1 + γ_2 + \frac{γ_0}{2} = γ_1 + γ_2 + \frac{1-γ_1-γ_2}{2} = \frac{1+γ_1+γ_2}{2},$

транзистор Т3 поддерживает состояние 110, он будет включен по времени

$γ_2 + \frac{γ_0}{2} = γ_2 + \frac{1-γ_1-γ_2}{2} = \frac{1-γ_1+γ_2}{2}$

транзистор Т5 поддерживает состояние 111, он будет включен

$\frac{γ_0}{2} = \frac{1-γ_1-γ_2}{2}$

Схемы замыкания ключей трехфазного инвертора и относительные напряжения

Через 180° формулы повторяются, т. е. в секторах 1 и 4, 2 и 5, 3 и 6 одинаковые вычисления, поэтому можно в два раза сократить модель, определяя не 6 секторов, а только 3

На рисунке приведены: трехфазный сигнал, по которому вычисляется угол и задается
амплитуда управляющего сигнала; по углу определяются секторы,
сформированные сигналы управления
Таким образом, общий алгоритм SVPWM выглядит следующим образом:

Определение текущего сектора угла θ.
Вычисление угла внутри сектора.
Расчет временных интервалов для каждой фазы.
Формирование управляющих сигналов для транзисторов и их распределение по ключам инвертора.

Новый алгоритм SVPWM

Традиционный алгоритм SVPWM требует преобразований координат, тригонометрической преобразований, векторного разложения в сектор и т. д., что требует больших затрат ресурсов CPU при управлении в реальном времени. В статье A Novel SVPWM Modulation Scheme предлагается новый алгоритм SVPWM, основанный на традиционном анализе SVPWM. В течение периода управления время пребывания базового вектора напряжения может быть рассчитано без преобразований координат и векторного разложения, что упрощает традиционный алгоритм SVPWM. Предложенный метод эффективен по результатам моделирования и эксперимента. Он показывает простоту программы кода и сокращение времени выполнения ЦП в системе управления в реальном времени.

В классическом методе трехфазные координаты преобразуются в систему координат α-β:

$\begin{bmatrix} V_{\alpha ref} \\ V_{\beta ref} \end{bmatrix} = \frac{2}{3}\begin{bmatrix} 1 & -\frac{1}{2} & -\frac{1}{2} \\ 0 & \frac{\sqrt{3}}{2} & -\frac{\sqrt{3}}{2} \end{bmatrix} \begin{bmatrix} V_a \\ V_b \\ V_c \end{bmatrix} \tag{1}$

V_a, V_b, V_c — трехфазное напряжение в трех координатах, $V_{αref}, V_{βref}$ - это $V_{ref}$ в координатах α-β . $V_{ref}$ можно получить через смежные векторы V_k и $V_{k+1}$ :

$V_{ref} = \frac{T_k}{T_s} V_k + \frac{T_{k+1}}{T_s}V_{k+1} \tag{2}$

где T_k и $T_{k+1}$ - время пребывания V_k и $V_{k+1}$ в течение каждого контрольного периода T_s , - номер сектора, угол вектора можно получить из двух оставшихся координатных систем. Время пребывания оценивается с помощью следующих уравнений в секторе I:

$\begin{aligned}T_1 &= \frac{T_s}{V_{dc}} |V_{ref}|(cosθ - \frac{sinθ}{\sqrt{3}}) \\T_2 &= \frac{2}{\sqrt{3}} \frac{T_s}{V_{dc}} |V_{ref}| sinθ\end{aligned} \tag{3}$

Из вышеприведенного следует, что большое количество вычислений потребует большого количества процессорного времени и снизит точность управления.
Амплитуда $V_{ref} (V_r)$ составляет 2/3 от фазных напряжений, которые выражаются следующим образом:

$\begin{aligned}V_a &= V_r \cos(\theta) \\V_b &= V_r \cos(\theta - \frac{2\pi}{3}) \\V_c &= V_r \cos(\theta + \frac{2\pi}{3})\end{aligned} \tag{4}$

Формулы из тригонометрии - функции суммы и разности углов

$\begin{aligned}cos(α – β) &= cos α · cos β + sin α · sin β \\sin(\frac{2}{3}\pi) &= \frac{\sqrt{3}}{3} \\cos(\frac{2}{3}\pi) &= -\frac{1}{2}\end{aligned} \tag{5}$

Пример для сектора I - Va > Vb > Vc, разность фазных напряжений:

$\begin{aligned}V_a - V_b &= \frac{2}{3}V_r(cosθ - cos(θ - \frac{2}{3}\pi)) \\&= \frac{2}{3}V_r(cosθ + \frac{1}{2} cosθ - \frac{\sqrt{3}}{3} sinθ) \\&= V_r(cosθ - \frac{1}{\sqrt{3}}sinθ) \\cosθ - \frac{1}{\sqrt{3}}sinθ &= \frac{1}{V_r}(V_a - V_b)\end{aligned} \tag{6}$

Таким же образом могут быть получены и другие выражения:

$\begin{aligned}V_b - V_c &= \frac{2}{3} - (cos(\theta - \frac{2}{3}\pi) - cos(\theta + \frac{2}{3}\pi)) \\&= \frac{2}{\sqrt{3}}V_r sin\theta \\sin\theta &=\frac{\sqrt{3}}{2V_r}(V_b- V_c)\end{aligned} \tag{7}$

Таким образом время пребывания может быть вычислено:

$\begin{aligned}T_1 = \frac{T_s }{V_{dc}}(V_a - V_b) \\T_2 = \frac{T_s }{V_{dc}}(V_b - V_c)\end{aligned} \tag{8}$

Улучшенным метод нахождения секторов использует простой критерий сравнения фазных напряжений, который можно свести в таблицу

Sector	θ	Va, Vb, Vc
I	[0, π/3)	Va > Vb > V
II	[π/3, 2π/3)	Vb > Va > V
III	[2π/3, π)	Vb > Vc > Va
IV	[π, 4π/3)	Vc > Vb > V
V	[4π/3, 5π/3)	Vc > Va > Vb
VI	[5π/3, 2π]	Va > Vc > Vb

В классическом методе полагая, что время открытого состояния вектора верхнего плеча taon, tbon, tcon соответствует фазам U, V, W, для цифрового преобразования они могут быть выражены как:

$\begin{aligned}taon &= tcon + T_k + T_{k+1}\\tbon &= tcon + T_{k+1} \\rcon &= (T_s - T_k - T_{k+1}) /2\end{aligned} \tag{9}$

Подставляя (8) в (9), получим:

$\begin{aligned}taon &= \frac{T_s}{2V_{dc}}[V_{dc} + (V_a - V_c)] \\tbon &= \frac{T_s}{2V_{dc}}[V_{dc} -(V_a - V_b) + (V_b - V_c)] \\rcon &= \frac{T_s}{2V_{dc}}[V_{dc} + (V_a - V_c)]\end{aligned}$

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

Sector	Phase duties
I \|\| IV	$t_{aon} = \frac{T_s}{2U_{dc}}(U_{dc} + (U_{aref} - U_{cref}))$
	$t_{bon} = \frac{T_s}{2U_{dc}}(U_{dc} + (2U_{bref} - U_{aref} - U_{cref}))$
	$t_{con} = \frac{T_s}{2U_{dc}}(U_{dc} + (U_{cref} - U_{aref}))$
II \|\| V	$t_{aon} = \frac{T_s}{2U_{dc}}(U_{dc} + (2U_{aref} - U_{bref} - U_{cref}))$
	$t_{bon} = \frac{T_s}{2U_{dc}}(U_{dc} + (U_{bref} - U_{cref}))$
	$t_{con} = \frac{T_s}{2U_{dc}}(U_{dc} + (U_{cref} - U_{bref}))$
III \|\| VI	$t_{aon} = \frac{T_s}{2U_{dc}}(U_{dc} + (U_{aref} - U_{bref}))$
	$t_{bon} = \frac{T_s}{2U_{dc}}(U_{dc} + (U_{bref} - U_{aref}))$
	$t_{con} = \frac{T_s}{2U_{dc}}(U_{dc} + (2U_{cref} - U_{aref} - U_{bref}))$

Частота коммутации инвертора должна быть значительно выше времени реакции двигателя, иначе скачкообразные изменения напряжения приведут к нестабильной работе. Как правило, частота ШИМ должна быть не менее чем в пять раз выше постоянной времени двигателя, чтобы избежать резких механических колебаний и пульсаций.

Пример реализации алгоритма SVPWM на Python

Для наглядности кода и уменьшения числа if-операторов в Python-скрипте используется функционал numpy.select (возвращает массив, составленный из элементов списка выбора, в зависимости от условий.) и numpy.where (возвращает элементы, выбранные из x или y в зависимости от условия.). Примеры из документации. Создается массив целых чисел от 0 до 5 (включительно), элементы меньше 3 отбрасываются, элементы больше 3 возводятся в квадрат, а элементы, не соответствующие ни одному из этих условий (ровно 3), заменяются значением по умолчанию 42.

import numpy as np

x = np.arange(6)
condlist = [x < 3, x > 3]
choicelist = [x, x**2]
np.select(condlist, choicelist, 42)

на выходе получается

array([ 0, 1, 2, 42, 16, 25])

Для элементов со значение менее 5 выводится значение самого элемента, а со значением более 5 выводится 10 умноженное на элемент

import numpy as np
a = np.arange(10)
np.where(a < 5, a, 10*a)

на выходе:

array([ 0,  1,  2,  3,  4, 50, 60, 70, 80, 90])

Реализация алгоритма SVPWM выполнена в классе SVPWM. Про NonlinearIOSystem и параметры t, x, u, params в duty_ref я расскажу позже, когда речь будет идти про Python control.
Я любезно попросил chat GPT перевести Python-скрипт в Scilab-скрипт, синтаксис которого во многом напоминает Matlab. В дальнейшем я использую Scilab-скрипт для функции пользовательского блока в графическом редакторе Xcos. В Matlab/ SImulink существует удобный визуальный редактор условий переходов и состояний Stateflow. В Scilab/Xcos для подобных целей с множеством условий переходов удобно использовать блок scifunc_block_m с которым сопоставляется скрипт. В арсенале есть блок AUTOMAT для реализации логики конечных автоматов, но я с ним не работал.

Scilab

Scilab — это бесплатное программное обеспечение с открытым исходным кодом для численных вычислений, предоставляющее вычислительную среду для инженерных и научных приложений.

Согласно описанию с оффсайта Scilab предоставляет следующую функциональность.

Maths & Simulation. Для инженерных и научных приложений, включая математические операции и анализ данных.
2-D & 3-D Visualization. Графические функции для визуализации, аннотирования и экспорта данных, а также множество способов создания и настройки различных типов графиков и диаграмм.
Optimization. Алгоритмы решения задач непрерывной и дискретной оптимизации с ограничениями и без ограничений.
Statistics. Инструменты для анализа данных и моделирования
Control Systems. Стандартные алгоритмы и инструменты для исследования систем управления
Signal Processing. Визуализация, анализ и фильтрация сигналов во временной и частотной областях.
Application Development. Расширение функционала Scilab и управление обменом данными с помощью внешних инструментов.
Xcos - Dynamic systems modeling. Моделирование механических систем, гидравлических цепей, систем управления.

Я выбрал Scilab/Xcos как ближайший свободный аналог Matlab/Simulink. В первую очередь меня интересовала возможность визуального проектирования систем управления, а также удобство использования данного инструмента. В целом Scilab справляется с поставленными задачами, и позволяет строить системы управления различной сложности.

Приведу субъективную оценку преимуществ и недостатков работы со Scilab.

Что понравилось:

свободное ПО;
работает на разных системах: Win, Lin, ARM;
размер на диске с некоторыми тулбоксами у меня занял 765 MB;
скорость запуска приложения и симуляции;
базовый функционал из-под коробки позволяет строить комплексные системы;
базовая концепция интерфейса и скриптовый язык похожи на Matlab;
Быстрое вхождение, особенно если есть опыт с подобными инструментами;
встроенная справка с примерами, по клику меню на блоке или операторе в тексте можно тут же вызвать документацию;
есть визуальный редактор диаграмм Xcos.

Что не понравилось:

на первом месте - это баги программы, с которыми я столкнулся
- вылет (закрытие) приложения во время симуляции;
- после запуска симуляции не выводятся графики результатов моделирования - помогает рестарт приложения;
- иногда после редактирования в Xcos возле всех блоков начинает мигать иконка с восклицательным знаком, иногда это не имеет значения, иногда симуляция не запускается - помогает рестарт;
- был случай когда в диаграмме Xcos вложенный суперблок просто не открывался на редактирование, но симуляция работала, пришлось заново перерисовать;
- на диаграмме связь между блоками может внезапно выглядеть в виде какой-то ломанной кривой, ни на что не влияет, но портит вид диаграммы - связь приходится перерисовывать
- баг с текстовым описанием диалога Format -> Edit -> Block legend, туда лучше ничего не писать, иначе текст улетает куда-то на диаграмме, и потом даже обычные текстовые метки куда-то пропадают
- иногда после редактирования блоков может появляться ошибка в консоли, связанное с индексами, симуляция не запускается - помогает рестарт;
на втором месте - устаревший пользовательский интерфейс, даже не сколько внешний вид, а именно удобство построения диаграмм по сравнению с подобными приложениями, можно долго рассказывать и придираться, но лучше попробовать и самому сделать вывод;
В суперблоках (аналог subsystem Matlab) отсутствует возможность текстовой маркировки портов (возможно я что-то не понял), поэтому я отдельно создавал текстовые метки, а затем перемещал их на соответствующие порты;
есть возможность установки тулбоксов, но большинство из них давно устарело и даже не устанавливается по причине несовместимости API;
если вы перешли с Matlab или других подобных коммерческих систем, вам будет не хватать тех или иных инструментов, но это больше как придирка к свободному ПО
отдельная установка и настройка C++ компилятора со своими нюансами для возможности моделирования электрических схем

Внешний вид основного окна приложения. Фон серого цвета я настроил отдельно.

Отдельное слово хочется сказать по источникам информации. Если поискать в сети книги, то основной период изданий - это 2003 - 2014 годы. Хотя недавно вышла книга Алексеев Е. Р "Scilab: Решение инженерных и математических задач", год издания 2024. По Xcos мне понравилась книга Arvind Kumar "Introduction to Xcos A Scilab Tool for Modeling Dynamical Systems" 2020. Судя по всему, Scilab пользуется популярностью в учебных заведениях, поэтому можно найти кучу материала в виде методических пособий, лабораторных работ и статей. Например, Учебное пособие по имитационному моделированию Scilab. Блог по системам управления Scilab Control Engineering Basics На youtube знаниями щедро делятся коллеги из Индии.

Python Control

Пакет Python Control предназначен для анализа и проектирования динамических систем в целом и систем управления с обратной связью в частности. Пакет напоминает Control System Toolbox в Matlab.
Пакет разработан в Калифорнийском технологическом институте (Caltech), США, проф. Ричардом М. Мюрреем и коллегами.
Пакет требует Numpy, Scipy и Matplotlib
Домашняя страница пакета Python Control https://pypi.org/project/control/
Примеры использования библиотеки https://github.com/python-control/python-control/tree/main/examples

Основные возможности:

Линейные системы ввода/вывода в пространстве состояний и частотной области
Моделирование, имитация и анализ нелинейных систем ввода/вывода
Алгебра блок-схем: последовательные, параллельные и обратные связи
Временная характеристика
Частотная характеристика: диаграммы Боде и Найквиста
Анализ управления: устойчивость

Установить пакет можно с помощью команды

pip install control

Создание передаточных функций

Способ 1: Использование переменной Лапласа

Один из наиболее интуитивных способов задания передаточной функции — использование переменной s, представляющей оператор Лапласа. Рассмотрим передаточную функцию:

$H(s) = \frac{2}{5s + 1}$

В Python код для создания такой функции выглядит так:

import control
s = control.tf('s')
H = 2 / (5 * s + 1)
print("H(s) =", H)

Вывод в консоли:

H(s) =
     2
-------------
5 s + 1

Этот способ удобен, так как он близок к математической записи и облегчает понимание.

Способ 2: Использование массивов коэффициентов

Передаточные функции также можно задавать массивами коэффициентов числителя и знаменателя. Например, тот же самый объект можно создать следующим образом:

import numpy as np
import control
num = np.array([2])
den = np.array([5, 1])
H = control.tf(num, den)
print("H(s) =", H)

Этот метод полезен, когда коэффициенты вычисляются динамически или загружаются из данных.

Комбинирование передаточных функций

Последовательное соединение

При последовательном соединении передаточных функций результирующая передаточная функция получается их произведением:

H1 = 2 / s
H2 = 3 / (4 * s + 1)
H = H1 * H2
print("H(s) =", H)

Вывод:

H(s) =
     6
---------------
4 s^2 + s

Также можно использовать функцию control.series():

H = control.series(H1, H2)

Параллельное соединение

При параллельном соединении передаточных функций их передаточные функции суммируются:

H = H1 + H2

Или с помощью control.parallel():

H = control.parallel(H1, H2)

Обратная связь

Для обратной связи с отрицательной обратной связью используется выражение:

H = H1 / (1 + H1 * H2)

Можно также воспользоваться встроенной функцией:

H = control.feedback(H1, H2)

По умолчанию control.feedback() предполагает отрицательную обратную связь.

Симуляция с передаточными функциями

Для моделирования систем с передаточными функциями в control используется функция control.forced_response(), а также специальные функции для типовых воздействий:

control.step_response() — отклик на ступенчатое воздействие
control.impulse_response() — импульсная характеристика

Пример моделирования системы с передаточной функцией $H(s) = \frac{2}{5s + 1}$ при входном ступенчатом сигнале амплитуды 4:

import numpy as np
import control
import matplotlib.pyplot as plt

s = control.tf('s')
H = 2 / (5 * s + 1)

# Setting the simulation time
t = np.linspace(0, 20, 2000)
u = 4 * np.ones_like(t)

# Modelling
_, y = control.forced_response(H, t, u)

# Plotting a graph
plt.subplot(2, 1, 1)
plt.plot(t, y, 'b', label='y')
plt.grid()
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(t, u, 'g', label='u')
plt.xlabel('t [s]')
plt.grid()
plt.legend()

plt.show()

График показывает, как система реагирует на входной сигнал во времени.

Дискретные передаточные функции

Передаточные функции во временной области z создаются аналогично, но с указанием шага дискретизации Ts:

Ts = 0.1
num = np.array([0.1, 0])
den = np.array([1, -1])
H = control.tf(num, den, Ts)
print("H(z) =", H)

Вывод:

0.1 z
-----
z - 1

dt = 0.1

Дискретизация передаточной функции

Для преобразования передаточной функции из непрерывного времени в дискретное используется control.sample_system():

Ts = 0.001
s = control.tf('s')
H = 1 / s
H_disc = control.sample_system(H, Ts, method='zoh')
print("H_disc(z) =", H_disc)

Вывод:

0.001
-----
z - 1

dt = 0.001

Метод zoh (zero-order hold) является стандартным способом дискретизации систем управления.

Динамическая система, реализующая сложные взаимосвязи

Для построения сложных взаимосвязанных систем удобно использовать функцию control.interconnect или класс control.InterconnectedSystem

Варианты использования этой функции можно посмотреть здесь
Например, сигнал автоматически направляется в каждую систему, имеющую вход с таким же именем.

     u        y1
u  +--> sys1 --->
---|          
   +--> sys2 --->
     u        y2

# arbitrary example systems
sys1 = ct.tf(1, [10, 1], inputs='u', outputs='y1')
sys2 = ct.tf(1, [1, 0], inputs='u', outputs='y2')

# create interconnected system
interconnected = ct.interconnect([sys1, sys2], inputs='u', outputs=['y1', 'y2'])
display(interconnected) # 1-input, 2-output system

Для этой системы входом служит один сигнал , в то время, как выход представлен в виде вектора с двумя элементами y1, y2 .

Второй пример, где складываются системы с одноименными выходными сигналами.

  u1        y
---> sys1 ---+ 
             |
           + V    y
             O----->
           + ^
             |
---> sys2 ---+
  u2        y

sys1 = ct.tf(1, [10, 1], inputs='u1', outputs='y')
sys2 = ct.tf(1, [1, 0], inputs='u2', outputs='y')

# create interconnected system
interconnected = ct.interconnect([sys1, sys2], inplist=['u1', 'u2'], outlist='y')
display(interconnected) #  2-input, 1-output system

Для суммирования входных сигналов или изменения знака сигнала можно использовать функцию control.summing_junction

 u      w  
---> O --->    
     ^
     | -v
     |

summer = ct.summing_junction(['u', '-v'], 'w') # w = u - v
display(summer)

Пример построения комплексной системы

# constants
K = 10
zc = 0.001 # controller zero location
pc = 2 # controller pole location
tau = 1
J = 100
b = 1

# systems
C = control.tf([K, K*zc],[1, pc], inputs='e', outputs='u')
lopass = control.tf(1, [tau, 1], inputs='u', outputs='v')
P = control.tf(1, [J, b], inputs='w', outputs='thetadot')
integrator = control.tf(1, [1, 0], inputs='thetadot', outputs='theta')
error = control.summing_junction(['thetaref', '-theta'], 'e') # e = thetaref-theta
disturbance = control.summing_junction(['d', 'v'], 'w') # w = d+v

# interconnect everything based on signal names
sys = control.interconnect([C, lopass, P, integrator, error, disturbance], 
                       inputs=['thetaref', 'd'], outputs='theta')
display(sys)

plt.figure(figsize=(7,3))
sys_thetaref_to_theta = sys[0, 0] 
sys_d_to_theta = sys[0, 1]
t, y = control.step_response(sys_thetaref_to_theta) # step response
plt.plot(t,y)
plt.figure(figsize=(7,3))
t, yd = control.step_response(sys_d_to_theta) # disturbance response
plt.plot(t,yd);
plt.figure(figsize=(7,5))
control.bode_plot(sys_thetaref_to_theta, plot=True);

Библиотека Python control позволяет строить нелинейные динамические системы с помощью класса classcontrol.NonlinearIOSystem(updfcn, outfcn=None, params=None, **kwargs). С анализом таких систем можно ознакомится на примере круиз-контроля авто здесь.
Практически любую Python-функцию возможно преобразовать в NonlinearIOSystem.
для этого сигнатура функции должна иметь следующий вид

def foo(t, x, u, params):
	pass

где

t (float) - время измерения
x (array_like) - текущее состояние
u (array_like) - вход
params (dict_,_ optional) - параметры системы

Рассмотрим, как используется NonlinearIOSystem на примере модели PMSM двигателя. В модели производиться динамический расчет продольного, квадратурного тока и угловой скорости. Входными сигналами являются задание скважности напряжений, обратная связь в виде текущих значения токов, угловой скорости и угла ротора duty_a, duty_b, i_d, i_q, w_m, alpha_e. В дискретной системе нам необходимо вычислить приращения токов и угловой скорости для дальнейшего интегрирования. Для этого из уравнений необходимо вычесть текущее состояния этих величин - x[i], и вернуть полученный результат, т.е.

$\Delta y = f(u) - x[i], \ где \ i = 0, 1, 2...$

Порядок индексов состояний соответствует тому, в каком порядке возвращаются величины, например return di_d, di_q, dw_m, p*dw_m - значит, что состояние x[0] = di_d в предыдущей итерации, x[1] = di_q и т.п.

def update(self, t, x, u, params):
        L_d = params.get('L_d', 0.0003)
        L_q = params.get('L_q', 0.0003)
        Psi_p = params.get('Psi_p', 0.095)
        p = params.get('p', 3)
        F = params.get("F", 0.0048)
        J = params.get("J", 0.01)
        R = params.get('R', 0.01)
        u_bat = params.get('u_bat', 400)

        duty_a, duty_b, i_d, i_q, w_m, alpha_e = u
        u_a = self.inverter(duty_a, u_bat)
        u_b = self.inverter(duty_b, u_bat)
        u_x, u_y = self.clark_transform(u_a, u_b)
        u_d, u_q = self.park_transform(u_x, u_y, alpha_e)

        w_e = p*w_m
        load_torque = 0 #TODO
        di_d = 1/L_d*(u_d - R*i_d + w_e*L_q*i_q) - x[0]
        di_q= 1/L_q*(u_q - R*i_q - w_e*L_d*i_d - Psi_p*w_e) - x[1]
        dw_m = (3/2*p*(Psi_p*i_q + (L_d - L_q)*i_d*i_q) - abs(w_m)*F - load_torque)/J - x[2]

        return di_d, di_q, dw_m, p*dw_m

Объект системы строиться следующим образом

self.diff = control.NonlinearIOSystem(
            self.update, self.state_output,
            inputs=["u_a", "u_b", "i_d", "i_q", "w_m", "alpha_e"], 
            outputs=('di_d', 'di_q', 'dw_m', 'w_e'), 
            states=('di_d', 'di_q', 'dw_m', 'w_e'), 
            name='diff',
            dt=self.Ts
        )

Если в объекте необходимо использовать состояния, то вычисляемая функция (в нашем примере - update) должна быть в качестве аргумента updfcn (по умолчанию первая позиция), аргумент outfcn (по умолчанию - вторая позиция) можно оставить None или задать некую функцию. Если состояния не предполагаются, то вычисляемую функцию можно использовать как аргумент outfcn, а оставить None.
В inputs перечисляются названия входных сигналов, порядок индексов входов вычисляемой функции (в нашем случае update) соответствует порядку следования названий в inputs.
Созданный объект NonlinearIOSystem затем используется для построения комплексной системы с обратными связями с помощью функции control.interconnect, или класса InterconnectedSystem. Связи между системами строятся с помощью аргумента connections, в котором в виде списков передаются названия входных сигналов одной системы и выходных сигналов другой системы.

self.sys = control.InterconnectedSystem(
            (
                self.sys_d,
                self.sys_q,
                self.ref_trans_sys,
                self.voltage_trans_sys
            ),
            connections = (
                ['sys_d.i_ref', 'ref_trans.i_d_ref'],
                ['sys_d.current', 'ref_trans.i_d'],
                ['sys_d.u_decoupl', 'ref_trans.u_d_decoup'],

                ['sys_q.i_ref', 'ref_trans.i_q_ref'],
                ['sys_q.current', 'ref_trans.i_q'],
                ['sys_q.u_decoupl', 'ref_trans.u_q_decoup'],

                ['voltage_trans.u_d', 'sys_d.u'],
                ['voltage_trans.u_q', 'sys_q.u'],
                ['voltage_trans.alpha_e', 'ref_trans.alpha_e'],
            ),
            inplist = ["ref_trans.t_mot_ref", "ref_trans.alpha_e", "ref_trans.w_e", "ref_trans.i_a", "ref_trans.i_b"],
            inputs = ["t_mot_ref", "alpha_e", "w_e", "i_a", "i_b"],
            outlist = ["voltage_trans.u_ref_a", "voltage_trans.u_ref_b", "voltage_trans.u_ref_c", "ref_trans.i_d", "ref_trans.i_q"],
            outputs = ["u_ref_a", "u_ref_b", "u_ref_c", "i_d", "i_q"],
            name = "current_control"
        )

Выводы использования Python control

Библиотека control предоставляет удобные инструменты для работы с передаточными функциями, их комбинирования и моделирования. Использование этой библиотеки облегчает анализ и проектирование динамических систем управления. Возможно сочетание с различными инструментами в Python. В процессе работы с багами не сталкивался. Не хватает только средств визуального проектирования диаграммы потоков.

Реализация модели системы управления PMSM в Scilab/Xcos и Python control

На рисунке ниже показана общая схема модели в Scilab/Xcos в виде связанных блоков. Диаграмма Xcos сохраняется в файл с расширением .zcos, который представляет собой обычный архив.

В качестве управляющего сигнала с заданной дискретизацией измерений является источник CCLOCK. Для использования единого источника в разных местах используются событийные порты (events ports) - ромб с красной каёмкой.

Период дискретизации источника CCLOCK задан в скрипте параметров. Для данной системы его значение должно быть не менее 0.001 с.

В блоках системы в качестве интегральных звеньев используются дискретные передаточные функции. В Xcos этот блок называется DLR - Discrete transfer function.

DLR содержит управляющий порт, на который необходимо подавать сигнал CLOCK. В числителе необходимо умножение на параметр дискретизации Ts.

Для отображения результатов симуляции в виде графиков используются блоки CSCOPE. Я расположил осциллоскопы в верхней части диаграммы. Cигналы подаются через референсные порты для уменьшения числа линий связи. Затем сигналы объединяются на мультиплексоре MUX и подключаются с осциллоскопу.

Рисунок - Общая схема системы управления

Параметры системы вынесены в отдельный скрипт data.sci для удобства централизованного изменения. Скрипт необходимо выполнить перед симуляцией, чтобы эти параметры загрузились как переменные окружения среды для дальнейшего использования в диаграмме.

Рисунок - скрипт с параметрами

В отдельный скрипт svpwm.sci вынесена логика SVPWM, который также необходимо загрузить перед запуском симуляции.

Генерация опорного тока - Current reference generation

Генерация опорного тока - это операция, в ходе которой заданный крутящий момент двигателя преобразуется в составляющие продольного i_d и квадратурного i_q тока. Этот процесс основан на аналитическом выражении для момента двигателя PMSM

$\begin{aligned}m_e &= \frac{3}{2}p(\Psi_pi_q + (L_d - L_q)i_di_q) \rightarrow i_q = \frac{2m_e}{3p(\Psi_p + (L_d - L_q)i_d)}{} \\m_e &= \frac{3}{2}p\Psi_pi_q \rightarrow i_q = \frac{2m_e}{3p\Psi_p}\end{aligned}$

Выходными параметрами блока являются опорные значения токов в прямом и квадратурном направлении, обозначаемые как i_d_ref и i_q_ref. В случае конструкции SPMSM, влияние продольного тока на выходной момент отсутствует, поэтому его опорное значение принимается равным нулю (i_d_ref = 0).
Алгоритм вычисления i_q_ref сводится к следующему:

делим заданный момент 2 x T_mot_ref на 3 × P × Ψ
Используем блок saturation для предотвращения превышения допустимого значения с пределами +/- I_s_max

Контроллер тока - Current controller

Рисунок - Схема блока контроллера тока

Контроллер тока требует следующие входные параметры:

Опорные значения тока в продольном и квадратурном направлениях.
Измеренные значения тока в продольном и квадратурном направлениях (они являются управляемыми физическими величинами).
Опция включения контроллера тока (переменная логического типа).
Измеренная электрическая скорость вращения двигателя (необходима для реализации функции развязки напряжения).

Выходными параметрами контроллера являются:

Управляемые значения продольного и квадратурного напряжений.
Эти значения в дальнейшем преобразуются обратно в трехфазную систему координат.

PI-регулятор

Для управления токами в данной системе применяется пропорционально-интегральный (PI) регулятор. Входным значением является разница между заданным и текущем значением продольного или квадратурного тока (на схеме присутствует два отдельных блока) $e = i_{ref} - i$ . Интегратор представлен в виде дискретной передаточной функции Ts/(z-1). В операторной форме

U(s) = Kp_c * e + Ki_c * 1/s * e

Функция anti-windup служит для предотвращения перегрузки интегратора, если контролируемый сигнал не может достичь требуемого значения, а выходная величина с интегратора превышает предел насыщения. Сигнал с выхода регулятора через дискретное звено задержки Ts/z поступает на вход anti-windup.

Пропорциональная часть представляет собой усилитель разницы токов с коэффициентом усиления Kp. Интегральная часть позволяет достичь плавного и устойчивого состояния контролируемых значений токов. На выходе PI-регулятора расположен блок saturation, в котором устанавливаются лимиты +/- dq_max_voltage контролируемого сигнала.

Реализация функции развязки напряжения

Для корректной работы PI-регулятора необходимо реализовать логику развязки напряжений. Данный процесс включает:

Получение измеренной электрической скорости двигателя.
Вычисление индуцированных напряжений в продольном и квадратурном направлениях.
Добавление рассчитанных значений к пропорциональной части PI-регулятора.

Входные данные для функции развязки:

Электрическая скорость двигателя (ωe).
Измеренные значения продольного (Id) и квадратурного (Iq) токов.
Значения индуктивностей в продольном (Ld) и квадратурном (Lq) направлениях.
Параметр поток статора (ψp).

Выходные данные:

Развязанное напряжение в продольном направлении (Ud_decoupled).
Развязанное напряжение в квадратурном направлении (Uq_decoupled).

Рисунок - Схема блока преобразования тока

Контроллер скорости - Frequency Control

Основной целью блока является стабилизация заданной скорости и вычисление опорного момента T_mot_ref для блока генерации опорного тока.

Рисунок - схема блока контроллера скорости

Входные параметры:

Заданная угловая скорость в RPM - w_m_ref
текущая угловая скорость скорость (в нашем случае вычисленная программно) w_m[rad/s]

PI-регулятор скорости

Сигналом ошибки является разница между заданной скоростью w_m_ref [RPM] и текущей угловой скоростью w_m[rad/s], умноженной на коэффициент rad_s_to_rpm = 60/2pi = 9.5493
Интегратор представлен в виде дискретной передаточной функции 1/s = T/2(z+1)/(z-1)*, которая производит дискретизацию трапецеидальным способом, теоретические более точнее, но можно применять и другие типы функций. В операторной форме результат с необходимыми коэффициентами выглядит

T_ref(s) = Kp_w * e + Ki_w * 1/s * e

Как и в регуляторе тока присутствует дискретная anti-windup функция, контролирующее переполнение и выходная функция насыщения.

Преобразователь напряжения - Voltage transformation

Основная цель преобразования напряжения — преобразование регулируемых значений продольного и квадратурного напряжения в симметричные трехфазные напряжения A, B и C. После успешного преобразования полученные опорные напряжения подаются на обмотки двигателя, замыкая контур управления.

Рисунок - блок преобразования напряжения

Входными параметрами для этой подсистемы будут:

Продольное и квадратурное напряжения, поступающие от контроллера тока.
Электрический угол ротора, необходимый для обратного преобразования Парка. Электрический угол представлен в радианах и обозначается как alpha_e.

Выходными значениями алгоритма преобразования напряжения будут:

Опорные напряжения для статорных фаз A, B и C. Для получения этих значений необходимы последовательные этапы преобразований: обратное преобразование Парка и обратное преобразование Кларка.

Обратное преобразование Парка - inverse_park_transform

Входами являются:

Продольное и квадратурное напряжения.

Выходными данными этой подсистемы станут напряжения в стационарной двухфазной системе координат (X, Y).

Для реализации алгоритма вычисляется синус и косинус электрического угла. Затем по матричным формулам выполняем преобразование:

Продольное напряжение умножаем на косинус угла.
Квадратурное напряжение умножаем на синус угла с учетом знаков согласно матричной формуле.
Складываем и вычитаем полученные значения для получения компонентов X и Y.

Обратное преобразование Кларка - inverse_clarke_transform

Входами являются напряжения в системе координат X и Y. Выходами — опорные напряжения фаз A, B и C.

Реализация преобразования включает следующие этапы:
2. Определение напряжения в фазе A: оно равно напряжению в направлении X.
3. Вычисление напряжения в фазе B:
- Умножаем X на коэффициент 1/2.
- Умножаем Y на коэффициент √(3)/2.
- Суммируем полученные значения.
4. Аналогично определяем напряжение в фазе C, используя разности полученных значений.

Векторная широтно-импульсная модуляция (SVPWM)

Модуль SVPWM представлен в виде блока scifunc_block_m, логика которого реализована в виде Scilab листинга scpwm.sci. Пример работы с данным блоком можно посмотреть в этом видео

Логику работы SVPWM я описывал ранее, здесь можно посмотреть Scilab листинг

Модель PMSM

Входные параметры модели

Средние фазные напряжения, поступающие от модели инвертора.
Нагрузочный момент двигателя

Выходные параметры модели

Измеренные (или вычисленные) фазные токи двигателя.
Механическую и электрическую скорость вращения.
Электрический угол двигателя, который рассчитывается в рамках модели и необходим для алгоритма управления FOC (Field Oriented Control).

Первый этап обработки данных — преобразование средних фазных напряжений в продольные (d) и квадратурные (q) составляющие. Для этого используется преобразование Кларка и Парка, переводящее координаты из трехфазной системы (A, B, C) в систему вращающихся координат (d, q). Эти преобразования являются неотъемлемой частью реализации алгоритма FOC. Зная эти компоненты, можно вычислить электромеханический момент, развиваемый двигателем.

Основные уравнения модели

Для описания PMSM используются четыре фундаментальных уравнения, выведенных на основе электрической модели двигателя. В этой модели примем эти уравнения как данность и будем оперировать ими с помощью математических преобразований для упрощения их реализации.

$\begin{aligned}u_d &= R_s i_d + s L_d i_d - w_e L_q i_q \\u_q &= R_si_q + s L_qi_q + w_e L_d i_d + w_e \Psi_p \\m_e &= \frac{3}{2}p(\Psi_pi_q + (L_d - L_q)) \\\sum{m} &= sJw_m\end{aligned}$

Первыми двумя уравнениями являются уравнения напряжения в продольной (d) и квадратурной (q) осях

Первые слагаемые выражают омическое сопротивление.
Вторые представляют собой наведённое напряжение, возникающее при изменении токов.
В установившемся режиме, когда двигатель вращается с постоянной скоростью, производные токов равны нулю, а значит, и соответствующие составляющие напряжения обнуляются.
В режиме разгона или торможения токи изменяются, что приводит к возникновению наведённого напряжения, которое необходимо учитывать.

Кроме того, в уравнениях присутствуют члены, зависящие от электрической угловой скорости двигателя , которая определяется как произведение механической скорости на число пар полюсов. Эти компоненты описывают наведённое напряжение в статоре вследствие вращения ротора.

Третье уравнение представляет собой выражение для электромагнитного момента двигателя, а четвёртое — уравнение движения. В последнем:

момент инерции двигателя.
механическая угловая скорость.
суммарный момент, действующий на двигатель (включает электромагнитный момент, момент нагрузки).

Для упрощения вычислений применяем дискретное преобразование Лапласа, которое позволяет заменить дифференциальные операторы на алгебраические выражения. Интегрирование представляется как дискретная передаточная функция $\frac{KT_s}{z-1}$
После преобразования мы можем выразить токи в осях d и q через известные параметры:

$\begin{aligned}i_d &= \frac{1}{s} \frac{u_d - R_si_d + w_eL_qi_q}{L_d} = \frac{KT_s}{z-1} \frac{u_d - R_si_d + w_eL_qi_q}{L_d} \\i_q &= \frac{1}{s} \frac{u_q - R_si_q - w_eL_di_d - \Psi_pw_e}{L_q} = \frac{KT_s}{z-1} \frac{u_q - R_si_q - w_eL_di_d - \Psi_pw_e}{L_q} \\m_e &= \frac{3}{2}p(\Psi_pi_q + (L_d - L_q)i_di_q) \\w_m &= \frac{1}{3} \frac{\sum{m}}{J} = \frac{KT_s}{z-1} \frac{\sum{m}}{J}\end{aligned}$

Электромеханический момент двигателя включает:

Компоненту реактивного момента, зависящую от разности индуктивностей в продольном и квадратурном направлениях.
Синхронную компоненту, определяемую произведением магнитного потока и квадратурного тока.

Для двигателей PMSM индуктивности в осях d и q могут быть равны, что приводит к нулевому значению реактивного момента. Однако модель позволяет изменять эти параметры для анализа их влияния на работу двигателя.

После вычисления электромеханического момента с помощью общего уравнения движения определяется механическая скорость вращения двигателя. Учитывая количество пар полюсов двигателя, рассчитывается электрическая скорость вращения.

Электрический угол определяется с помощью дискретного интегратора электрической скорости. Этот угол используется в обратном преобразовании Кларка и Парка для получения трехфазных значений токов (A, B, C).

Исходный код реализации модели PMSM на Python control можно посмотреть здесь.

Параметры модели

В качестве примера выберем PMSM с параметрами:

Максимальная мощность: 190 кВт.
Индуктивности = 0,3 мГн.
Сопротивление статора: 10 мОм.
Число пар полюсов: 3.
Номинальное напряжение: 400 В.
Номинальная скорость: 11 000 об/мин.
Максимальный крутящий момент: ≈170 Н·м.
Постоянная обратной ЭДС: 36,9 мВ/об/мин.
Максимальный ток: ≈390 А.

Эти параметры закладываются в модель для дальнейшего анализа и тестирования алгоритмов управления.

В результате расчета получаем значения:

Фазных токов A, B и C.
Механической и электрической скорости вращения.
Электрического угла.

Симуляция системы управления

В Scilab / Xcos запуск симуляции можно запустить, нажав кнопку Start на панели инструментов.

Или из консоли, выполнив команды:

загрузка диаграммы

importXcosDiagram("FieldOrientedControl.zcos")

запуск симуляции
```
xcos_simulate(scs_m, 4);
```

Настройка параметров симуляции производится в диалоговом окне вызовом меню Simulation -> setup

В настройках необходимо задать период симуляции - Final integration time, в данном примере - 2 с. Остальные параметры чаще всего - без изменений.

Результаты симуляции Scilab / Xcos

Время симуляции - 2с
Сигнал ступенчатый:
t = [0, 1, 2]
in = [6000, 3000, 0]
Сигнал формируется с помощью блока Signal builder

Результаты моделирования с заданными параметрами.

Рисунок - Сигнал Duty SPWM за период 0 - 0.75 с

Угловые скорости за период 0 - 2 с

Время симуляции - 2с
Сигнал - линейно-нарастающий:
t = [0...2]
in = [0...6000]
параметры регуляторов те же

Рисунок - Угловые скорости за период 0 - 2 с

Время симуляции с данными параметрами занимает около 30 с.
Меняя параметры PI-регуляторов, структуру звеньев регуляторов и обратной связи можно достигнуть оптимального режима работы двигателя и преобразователя, уменьшить перерегулирование, пусковые токи.

Симуляция Python control

Для запуска симуляции необходимо задать временной интервал в виде списка numpy array. Создать комплексную систему с обратной связью с выхода PMSM на вход контроллера тока. Запустить симуляцию, выполнив функцию control.input_output_response. В результате получим массив с временным интервалом и набор массивов выходных сигналов. Построение графиков производиться стандартными методами Matplotlib. Время симуляции у меня заняло около 2.5 минуты.

Заключение

В данной статье рассмотрен процесс моделирования системы векторного управления синхронным двигателем с постоянными магнитами (PMSM) с использованием открытых программных инструментов, таких как Python control и Scilab. Были изучены основные принципы векторного управления (Field Oriented Control, FOC), включая преобразования координат, регулирование тока и скорости, а также алгоритмы пространственно-векторной ШИМ - SVPWM.
Судя по научным публикациям исследования в данном направлении продолжаются, что позволяет оптимизировать алгоритмы управления системой.
Параметры регуляторов значительно влияют на качество управления. Важно учитывать характеристики двигателя и настраивать коэффициенты регуляторов в зависимости от динамических требований системы.
Использование Python и Scilab позволяет построить полноценную модель системы управления, однако необходимо учитывать удобство инструментов, их развитие и поддержку.
Модель системы полезна для прототипирования и отладки алгоритмов управления перед их внедрением в микроконтроллеры. Кроме того, такой подход позволяет быстро тестировать различные стратегии управления без необходимости сборки физического стенда.
Дальнейшие исследования могут быть направлены на имитационное моделирование в среде OpenModelica. Где в отличие от традиционного подхода диаграммы потоков используется объектно-ориентированное моделирование, что упрощает построение сложных многокомпонентных систем.

Используемые источники

Комментарии (4)

almaz1c
14.02.2025 19:40
#27923174
Плюсанул статью и карму. Пишите ещё!

П.с. как ни пытался понять ТАУ в универе, так и не смог. Незакрытый гештальд.
1. vladipirogov Автор
  14.02.2025 19:40
  #27924316
  Спасибо). Планирую привести пример в Open Modelica. Там другой подход к моделированию. От знания физики процессов не избавляет, но вместо диаграммы потоков используется имитационное моделирование. В нем модель системы строиться по форме близко к реальной системе. Максимум, придется строить уравнения состояний, а в ряде случаев просто использовать готовые библиотеки моделей

Andy_U
14.02.2025 19:40
#27924672
А вы сравнивали результаты вашего моделирования с каким-нибудь реальным электродвигателем?
1. vladipirogov Автор
  14.02.2025 19:40
  #27924726
  Да, я смотрел графики реальных систем управления. Многое зависит от структуры системы, выбора коэффициентов. Мои результаты симуляции не достаточно оптимальные. Возможно стоит рассмотреть вопрос оптимизации систем, исследовать, какие сейчас разрабатываются методики в этом направлении. В предыдущей статье "Анализ кривой падения добычи нефтяных и газовых скважин" я постарался показать, как с помощью инструмента SciPy optimization curve_fit производиться подгонка коэффициентов уравнения для прогноза. Классическими методами настройки ПИД-регулятора считаются правила Зиглера-Никольса или Коэна-Куна. Вопрос, действительно, интересный и важный