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

Создание систем

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

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

tf — рациональное представление

Синтаксис для создания передаточной функции с помощью tf имеет следующий вид.

tf(num, den)     # Система с непрерывным временем
tf(num, den, Ts) # Система с дискретным временем

Здесь num и den — это полиномиальные коэффициенты числителя и знаменателя многочлена, а Ts (если указано) — интервал дискретизации для системы с дискретным временем.

Пример:

tf([1.0],[1,2,1])

# output

TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
        1.0
-------------------
1.0s^2 + 2.0s + 1.0

Continuous-time transfer function model

Создаваемые таким способом передаточные функции имеют тип TransferFunction{SisoRational}.

zpk — представление полюсов, нулей и усиления

Иногда передаточную функцию удобнее представлять в виде ее полюсов, нулей и коэффициента усиления. Для этого служит функция zpk.

zpk(zeros, poles, gain)     # Система с непрерывным временем
zpk(zeros, poles, gain, Ts) # Система с дискретным временем

Здесь zeros и poles — это векторы (Vectors) нулей и полюсов системы, а gain — коэффициент усиления.

Пример

zpk([-1.0,1], [-5, -10], 2)

# output

TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, Float64}}
   (1.0s + 1.0)(1.0s - 1.0)
2.0-------------------------
   (1.0s + 5.0)(1.0s + 10.0)

Continuous-time transfer function model

Создаваемые таким способом передаточные функции имеют тип TransferFunction{SisoZpk}.

Системы пространства состояний

Система пространства состояний создается следующим образом:

ss(A,B,C,D)    # Система с непрерывным временем
ss(A,B,C,D,Ts) # Система с дискретным временем

Работать с такой системой можно так же, как с передаточными функциями.

Конструктор ss позволяет делать следующее:

  • передавать 0 вместо матрицы ; нулевая матрица соответствующего размера создается автоматически;

  • передавать I вместо матрицы ; матрица тождественности соответствующего размера создается автоматически. Оператор I UniformScaling предоставляется стандартной библиотекой LinearAlgebra, которую нужно предварительно загрузить.

Также доступны системы пространства состояний с неоднородными типами матриц, которые можно использовать для создания систем со статическими или размерными матрицами, например:

using ControlSystemsBase, StaticArrays
sys = ss([-5 0; 0 -5],[2; 2],[3 3],[0])
HeteroStateSpace(sys, to_sized)
HeteroStateSpace(sys, to_static)
StaticStateSpace{Continuous, 1, 1, 2, StaticArraysCore.SMatrix{2, 2, Int64, 4}, StaticArraysCore.SMatrix{2, 1, Int64, 2}, StaticArraysCore.SMatrix{1, 2, Int64, 2}, StaticArraysCore.SMatrix{1, 1, Int64, 1}}
A =
 -5   0
  0  -5
B =
 2
 2
C =
 3  3
D =
 0

Continuous-time state-space model

Обратите внимание на использование матриц различных типов.

Сведения о связывании имен с состояниями, входами и выходами см. в описании named_ss из пакета RobustAndOptimalControl.jl.

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

Иногда бывает полезно выполнить преобразование из одного представления в другое. Это можно сделать с помощью конструкторов tf, zpk, ss, например:

tf(zpk([-1], [1], 2, 0.1))

# output

TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Int64}}
2z + 2
------
z - 1

Sample Time: 0.1 (seconds)
Discrete-time transfer function model

Системы с запаздыванием

Конструктор delay создает чистое запаздывание, которое можно ввести в систему путем умножения:

delay(1.2)               # Чистое запаздывание или 1,2 с
tf(1, [1, 1])*delay(1.2) # Запаздывание входа
delay(1.2)*tf(1, [1, 1]) # Запаздывание выхода

Системы с запаздыванием также можно создавать следующим образом:

s = tf("s")
L = 1.2 # Время запаздывания
tf(1, [1, 1]) * exp(-L*s)

Аппроксимацию Паде для запаздывания можно создать с помощью функции pade.

Руководство по системам с запаздыванием доступно здесь:

Нелинейные системы

См. раздел Nonlinear functionality.

Упрощение систем

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

Примеры:

julia> using ControlSystemsBase


julia> G = tf([1, 1], [1, 1])
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Int64}}
s + 1
-----
s + 1

Continuous-time transfer function model

julia> minreal(G) # Производит устранение нулевых значений сигнала
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0
---
1.0

Continuous-time transfer function model

julia> P = tf(1, [1, 1]) |> ss
StateSpace{Continuous, Float64}
A =
 -1.0
B =
 1.0
C =
 1.0
D =
 0.0

Continuous-time state-space model

julia> G = P / (1 + P) # здесь создается неминимальная реализация, поэтому используем вместо этого feedback(P)
StateSpace{Continuous, Float64}
A =
 -1.0  -1.0
  0.0  -2.0
B =
 1.0
 1.0
C =
 1.0  0.0
D =
 0.0

Continuous-time state-space model

julia> feedback(P) # Создает минимальную реализацию напрямую
StateSpace{Continuous, Float64}
A =
 -2.0
B =
 1.0
C =
 1.0
D =
 0.0

Continuous-time state-space model

julia> Gmin = minreal(G) # таким образом реализация упрощается до минимальной
StateSpace{Continuous, Float64}
A =
 -1.9999999999999993
B =
 -1.4142135623730951
C =
 -0.7071067811865472
D =
 0.0

Continuous-time state-space model

julia> norm(Gmin - feedback(P), Inf) # Нет отличий
1.110224134848181e-16

julia> bodeplot([G, Gmin, feedback(P)]) # Все одинаково
Plot{Plots.GRBackend() n=6}

Умножение систем

Две системы можно соединить последовательно путем умножения.

using ControlSystemsBase
P1 = ss(-1,1,1,0)
P2 = ss(-2,1,1,0)
P2*P1
StateSpace{Continuous, Int64}
A =
 -2   1
  0  -1
B =
 0
 1
C =
 1  0
D =
 0

Continuous-time state-space model

Если входное измерение P2 не соответствует выходному измерению P1, выдается ошибка. Если одна из систем одномерная, а другая — двухмерная, транслированное умножение расширяет одномерную систему так, чтобы она соответствовала входному или выходному измерению многомерной, например:

Pmimo = ssrand(2,2,1)
Psiso = ss(-2,1,1,0)
# Psiso * Pmimo # ошибка
Psiso .* Pmimo ≈ [Psiso 0; 0 Psiso] * Pmimo # Транслированное умножение расширяет одномерную систему до диагональной
true

Транслированное умножение системы и массива допустимо только для диагональных массивов.

using LinearAlgebra
Psiso .* I(2)
StateSpace{Continuous, Int64}
A =
 -2   0
  0  -2
B =
 1  0
 0  1
C =
 1  0
 0  1
D =
 0  0
 0  0

Continuous-time state-space model

Многомерные системы и массивы систем

В результате конкатенации систем образуется многомерная система, которая отличается от массива систем. Например:

using ControlSystemsBase
P = ss(-1,1,1,0)
P_MIMO = [P 2P]
StateSpace{Continuous, Int64}
A =
 -1   0
  0  -1
B =
 1  0
 0  1
C =
 1  2
D =
 0  0

Continuous-time state-space model

это многомерная система 1×2, а не массив 1×2.

Из одномерной системы в многомерную

Одномерные системы не умножаются на многомерные напрямую:

using Test
siso = ss(-1,1,1,0)
mimo = ssrand(2,2,2)
@test_throws DimensionMismatch siso * mimo
Test Passed
      Thrown: DimensionMismatch

Чтобы в приведенном выше примере умножить siso на каждый выходной канал mimo, используйте трансляцию:

siso .* mimo
StateSpace{Continuous, Float64}
A =
 -1.0   0.0   0.6480001475444149  -0.5837277961193962
  0.0  -1.0  -1.7039316164712754   1.5401479520865413
  0.0   0.0  -2.115494351059544   -1.8761537345884445
  0.0   0.0  -1.3413422781604374  -1.4565866982414075
B =
 -0.9508832623720646   -1.5870294046946989
  0.37005783730851793   1.069217230146147
 -0.04176136936679013   0.1320280925972221
  0.15655430568685796  -1.2664162819294336
C =
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
D =
 0.0  0.0
 0.0  0.0

Continuous-time state-space model

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

using LinearAlgebra
(siso .* I(2)) * mimo
StateSpace{Continuous, Float64}
A =
 -1.0   0.0   0.6480001475444149  -0.5837277961193962
  0.0  -1.0  -1.7039316164712754   1.5401479520865413
  0.0   0.0  -2.115494351059544   -1.8761537345884445
  0.0   0.0  -1.3413422781604374  -1.4565866982414075
B =
 -0.9508832623720646   -1.5870294046946989
  0.37005783730851793   1.069217230146147
 -0.04176136936679013   0.1320280925972221
  0.15655430568685796  -1.2664162819294336
C =
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
D =
 0.0  0.0
 0.0  0.0

Continuous-time state-space model

Преобразование массива систем в многомерную систему

Диагональную многомерную систему можно создать из вектора систем с помощью функции append.

P1 = ssrand(1,1,1)
P2 = ssrand(1,1,1)
append(P1, P2)
StateSpace{Continuous, Float64}
A =
 -0.031445085752152435   0.0
  0.0                   -0.11431097869771967
B =
 -0.6661593733640553   0.0
  0.0                 -0.2871155475506572
C =
 1.2299213185896871   0.0
 0.0                 -0.5861950130060937
D =
 -1.7142289926543897  0.0
  0.0                 2.0816093031690217

Continuous-time state-space model

Более общие массивы систем можно преобразовывать в многомерные системы с помощью метода array2mimo.

sys_array = fill(P, 2, 2) # Создает массив систем
mimo_sys = array2mimo(sys_array)
StateSpace{Continuous, Int64}
A =
 -1   0   0   0
  0  -1   0   0
  0   0  -1   0
  0   0   0  -1
B =
 1  0
 0  1
 1  0
 0  1
C =
 1  1  0  0
 0  0  1  1
D =
 0  0
 0  0

Continuous-time state-space model

Преобразование многомерной системы в массив систем

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

P = ssrand(2,3,1) # Случайная многомерная система 2×3
sys_array = getindex.(Ref(P), 1:P.ny, (1:P.nu)')
2×3 Matrix{StateSpace{Continuous, Float64}}:
 StateSpace{Continuous, Float64}
A =
 -0.9561216469448534
B =
 -0.060187565566125686
C =
 -0.8266381980246545
D =
 1.1416543199147604

Continuous-time state-space model   …  StateSpace{Continuous, Float64}
A =
 -0.9561216469448534
B =
 -0.8979659719080596
C =
 -0.8266381980246545
D =
 -0.41416330898904435

Continuous-time state-space model
 StateSpace{Continuous, Float64}
A =
 -0.9561216469448534
B =
 -0.060187565566125686
C =
 -0.9174184153688885
D =
 -1.5772540039273821

Continuous-time state-space model     StateSpace{Continuous, Float64}
A =
 -0.9561216469448534
B =
 -0.8979659719080596
C =
 -0.9174184153688885
D =
 -0.3751748796599928

Continuous-time state-space model

Создание массивов с различными типами систем

При вызове hcat/vcat Julia пытается автоматически продвинуть типы до наименьшего общего супертипа, то есть создать массив с одной непрерывной и одной дискретной системами не удастся.

P_cont = ssrand(2,3,1)
P_disc = ssrand(2,3,1, Ts=1)
@test_throws ErrorException [P_cont, P_disc] # ERROR: Sampling time mismatch
Test Passed
      Thrown: ErrorException

Вы можете явным образом сообщить Julia нужный супертип, например:

StateSpace[P_cont, P_disc]
2-element Vector{StateSpace}:
 StateSpace{Continuous, Float64}
A =
 -1.3716716788894652
B =
 -0.3624875727629689  0.9644915820182179  0.6987358102174509
C =
  0.8623674722193923
 -0.2178013242366765
D =
 0.5138344957212424   -0.054813542456689734  -1.5355805182590119
 0.09198259556814194   1.2387723604810115    -0.09329636744953426

Continuous-time state-space model
 StateSpace{Discrete{Int64}, Float64}
A =
 0.9
B =
 0.310562326219272  -1.5299520100783581  0.6793697708186918
C =
 -0.3827715076233991
  0.3127849910381858
D =
 -0.9475730778435572   0.025282874604404964  -0.11929165323769904
 -0.34601203395549096  0.18186138919749933    0.7078840093774816

Sample Time: 1 (seconds)
Discrete-time state-space model

Тип StateSpace является абстрактным, так как его параметры не заданы.

Демонстрационные системы

Модуль ControlSystemsBase.DemoSystems содержит ряд систем, демонстрирующих различные виды динамических процессов.

От блок-схем к коду

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

Функцию feedback(G1, G2) можно представить так: первый аргумент G1 — это система, которая находится непосредственно между входом и выходом (прямой путь), а второй аргумент G2 (по умолчанию 1, если он не указан) содержит все остальные системы, из которых состоит замкнутый контур (путь обратной связи). Предполагается, что обратная связь отрицательная, если не передан аргумент pos_feedback = true (исключением является функция lft, которая согласно соглашению по умолчанию обладает положительной обратной связью). Это означает, что feedback(G, 1) даст отрицательную единичную обратную связь, а feedback(G, -1) или feedback(G, 1, pos_feedback = true) — положительную.

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

r   ┌─────┐     ┌─────┐
───►│     │  u  │     │ y
    │  C  ├────►│  P  ├─┬─►
 -┌►│     │     │     │ │
  │ └─────┘     └─────┘ │
  │                     │
  └─────────────────────┘

Код: feedback(P*C) или, что то же самое, comp_sensitivity(P, C). В данном случае система находится непосредственно между входом и выходом , а контур обратной связи представляет собой отрицательное тождество. Поэтому мы вызываем feedback(P*C) = feedback(P*C, 1).

d     ┌───┐   y
───+─►│ P ├─┬───►
  -▲  └───┘ │
   │        │
   │  ┌───┐ │
   └──┤ C │◄┘
      └───┘

Код: feedback(P, C) или, что то же самое, G_PS(P, C). В данном случае только находится непосредственно между а , в то время как находится на первом месте в контуре обратной связи. Поэтому мы вызываем feedback(P, C).

Функция чувствительности на входе объекта

d    e┌───┐
───+─►│ P ├─┬───►
  -▲  └───┘ │
   │        │
   │  ┌───┐ │
   └──┤ C │◄┘
      └───┘

Код: feedback(1, C*P) или, что то же самое, input_sensitivity(P, C). В данном случае непосредственно между входом и выходом нет систем, поэтому мы вызываем feedback(1, C*P). Обратите внимание на порядок в выражении C*P: он важен для многомерных систем. Таким образом вычисляется функция чувствительности на входе объекта. Функция чувствительности чаще анализируется на выходе объекта, как показано ниже (для одномерных систем ситуация аналогичная).

Функция чувствительности на выходе объекта

      ┌───┐
───+─►│ P ├─+◄── e
  -▲  └───┘ │
   │        │y
   │  ┌───┐ │
   └──┤ C │◄┘
      └───┘

Код: feedback(1, P*C) или, что то же самое, output_sensitivity(P, C). Обратите внимание на то, что порядок в выражении обратный в сравнении с приведенной выше функцией чувствительности на входе.

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

              d
     ┌─────┐  │  ┌─────┐
r    │     │u ▼  │     │ y
──+─►│  C  ├──+─►│  P  ├─┬─►
  ▲  │     │     │     │ │
 -│  └─────┘     └─────┘ │
  │                      │
  └──────────────────────┘

Код: feedback(C, P, W2=:, Z2=:, Zperm=[(1:P.ny).+P.nu; 1:P.nu]) # y,u from r,d. Для получения управляющего сигнала с правильным знаком в данном случае порядок следования P и C был изменен на обратный. Мы также использовали именованные аргументы W2 и Z2, чтобы указать, что входы и выходы P нужно включить как внешние входы и выходы, а также аргумент Zperm, чтобы указать порядок выходов ( перед ).

Дробно-линейное преобразование

     ┌─────────┐
z◄───┤         │◄────w
     │    P    │
y┌───┤         │◄───┐u
 │   └─────────┘    │
 │                  │
 │      ┌───┐       │
 │      │   │       │
 └─────►│ K ├───────┘
        │   │
        └───┘

Код: lft(P, K)

      z1          z2
      ▲  ┌─────┐  ▲      ┌─────┐
      │  │     │  │      │     │
w1──+─┴─►│  C  ├──┴───+─►│  P  ├─┐
    │    │     │      │  │     │ │
    │    └─────┘      │  └─────┘ │
    │                 w2         │
    └────────────────────────────┘

Передаточная функция от к содержит все передаточные функции, которые обычно называются «большой четверкой» (см. также описание метода gangoffour).

Код: Для этой функции требуется пакет RobustAndOptimalControl.jl.

RobustAndOptimalControl.extended_gangoffour(P, C, pos=true)
# Для одномерного P
S  = G[1, 1]
PS = G[1, 2]
CS = G[2, 1]
T  = G[2, 2]

# Для многомерного P
S  = G[1:P.ny,     1:P.nu]
PS = G[1:P.ny,     P.nu+1:end]
CS = G[P.ny+1:end, 1:P.nu]
T  = G[P.ny+1:end, P.nu+1:end]

Дополнительные материалы

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

                 yF
              ┌────────────────────────────────┐
              │                                │
    ┌───────┐ │  ┌───────┐ yR   ┌─────────┐    │    ┌───────┐
uF  │       │ │  │       ├──────►         │ yC │  uP│       │    yP
────►   F   ├─┴──►   R   │      │    C    ├────+────►   P   ├────┬────►
    │       │    │       │   ┌──►         │         │       │    │
    └───────┘    └───────┘   │- └─────────┘         └───────┘    │
                             │                                   │
                             └───────────────────────────────────┘

См. пример кода complicated_feedback.jl.

Синтезирование фильтров

Фильтры можно синтезировать с помощью DSP.jl. В результате получаются объекты фильтров, типы которых определены в пакете DSP. Их можно преобразовывать в передаточные функции с помощью функции tf из ControlSystemsBase.

using DSP, ControlSystemsBase, Plots

fs = 100
df = digitalfilter(Bandpass(5, 10; fs), Butterworth(2))
G = tf(df, 1/fs) # Для получения правильной шакалы частот на графике Боде необходимо указать интервал дискретизации.
bodeplot(G, xscale=:identity, yscale=:identity, hz=true)

Дополнительные материалы