Сообщество Engee
Engee Community
/
Примеры и проекты
/
Зимняя школа Julia - День 5

Зимняя школа Julia - День 5

Автор
avatar-daryadarya
Соавторы
avatar-alexevsalexevs
Notebook

Работа с ControlSystems.jl

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

Для работы с системами управления на Julia мы будем пользоваться пакетом ControlSystems.jl. Если у вас в вашем пространстве Engee не был установлен этот пакет, то установить его можно с помощью пакетного менеджера Pkgи функции Pkg.add. Затем подключите пакет используя using.

In [ ]: 
import Pkg; Pkg.add("ControlSystems")
In [ ]: 
using ControlSystems

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

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

In [ ]: 
tf([1], [1, 2])   #   коэффициенты числителя и знаменателя
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Int64}}
  1
-----
s + 2

Continuous-time transfer function model

Задать передаточную функцию можно используя переменную Лапласа.

In [ ]: 
s = tf("s")     #   переменная Лапласа
P = (s + 0.5) / (s^2 + 2s + 1)
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
    1.0s + 0.5
-------------------
1.0s^2 + 2.0s + 1.0

Continuous-time transfer function model

Создадим передаточную функцию с помощью еще одной стандартной функции пакета pid().

In [ ]: 
C = pid(1,2)
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0s + 0.5
----------
   1.0s

Continuous-time transfer function model

Если у нас уже есть система, которая описана некоторой передаточной функцией и нам нужно узнать числитель или знаменатель. Мы можем использовать функцию numvec() для числителя и denvec() для знаменателя. Эти функции возвращают матрицу векторов, поскольку системы зачастую могут иметь несколько входов и выходов. Чтобы получить доступ к вектору коэффициентов нашей передаточной функции используем такую запись: numvec(P)[1].

In [ ]: 
numvec(P)       #   коэффициенты числителя ПФ
Out[0]:
1×1 Matrix{Vector{Float64}}:
 [1.0, 0.5]
In [ ]: 
denvec(P)       #   коэффициенты знаменателя ПФ
Out[0]:
1×1 Matrix{Vector{Float64}}:
 [1.0, 2.0, 1.0]

Для получения фактического представления многочлена числителя передаточной функции используем функцию numpoly(P). Она возвращает матрицу многочленов. Однако члены здесь представлены в обратном порядке. Сначала идет член с меньшим порядком, затем с большим.

In [ ]: 
numpoly(P)
Out[0]:
1×1 Matrix{Polynomials.Polynomial{Float64, :x}}:
 Polynomial(0.5 + 1.0*x)

Операции и структурные предобразования

Мы также можем умножать передаточную функцию на скаляр.

In [ ]: 
2P              #   умножение на скаляр
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
    2.0s + 1.0
-------------------
1.0s^2 + 2.0s + 1.0

Continuous-time transfer function model

И получать обратное значение, используя специальную функцию inv или запись 1/P.

In [ ]: 
inv(P)          #   обратная функция
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0s^2 + 2.0s + 1.0
-------------------
    1.0s + 0.5

Continuous-time transfer function model

Для структурных преобразований можно использовать стандартные функции пакета или использовать операторы “+” и “*”. Выражение P + C реализует параллельное подключение моделей. А P*C - последовательное подключение передаточных функций.

gsconnectingmodels_02.png

In [ ]: 
P + C           #   подключение моделей в параллельной форме
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0s^3 + 3.5s^2 + 2.5s + 0.5
----------------------------
   1.0s^3 + 2.0s^2 + 1.0s

Continuous-time transfer function model

gsconnectingmodels_01.png

In [ ]: 
P * C           #   подключение моделей в последовательной форме
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
 1.0s^2 + 1.0s + 0.25
----------------------
1.0s^3 + 2.0s^2 + 1.0s

Continuous-time transfer function model

Для получения системы с обратной связью в пакете есть функция feedback(). Однако можно использовать формулу для замкнутой системы. Но в результате мы не получим минимальной реализации и вычисляться такое выражение будет дольше.

In [ ]: 
G1 = P*C / (1 + P*C)    #   создание соединиения обратной связи с использованием "/"
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
    1.0s^5 + 3.0s^4 + 3.25s^3 + 1.5s^2 + 0.25s
---------------------------------------------------
1.0s^6 + 5.0s^5 + 9.0s^4 + 7.25s^3 + 2.5s^2 + 0.25s

Continuous-time transfer function model

Получить минимальную реализацию из G1 можно используя функции: minreal и sminreal.

In [ ]: 
minreal(G1)     #   Получим структурно минимальную реализацию функции G1
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
            1.0s^4 + 2.9999999999999996s^3 + 3.250000000000001s^2 + 1.500000000000002s + 0.2500000000000005
-----------------------------------------------------------------------------------------------------------------------
1.0s^5 + 5.000000000000001s^4 + 9.000000000000002s^3 + 7.250000000000003s^2 + 2.4999999999999987s + 0.24999999999999986

Continuous-time transfer function model
In [ ]: 
G2 = feedback(P*C)      #   также реализует модель с обратной связью; рекомендуемый способ образование замкнутого контура
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
    1.0s^2 + 1.0s + 0.25
-----------------------------
1.0s^3 + 3.0s^2 + 2.0s + 0.25

Continuous-time transfer function model

Если объединить две системы, используя вертикальную конкатенацию [P; C], будет создана MIMO система с двумя выходами. Если объединить системы горизонтально, то получится MIMO система с двумя входами. И если объединим системы по вертикали и по горизонтали, то получим систему с двумя входами и двумя выходами.

In [ ]: 
[P; C]          #   конкатенация моделей; создание MIMO системы
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
Input 1 to output 1
    1.0s + 0.5
-------------------
1.0s^2 + 2.0s + 1.0

Input 1 to output 2
1.0s + 0.5
----------
   1.0s

Continuous-time transfer function model

Результирующая модель имеет два выхода и один вход и соответствует следующей блок-схеме:

gsconnectingmodels_06.png

In [ ]: 
[P C]
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
Input 1 to output 1
    1.0s + 0.5
-------------------
1.0s^2 + 2.0s + 1.0

Input 2 to output 1
1.0s + 0.5
----------
   1.0s

Continuous-time transfer function model

Результирующая модель имеет два входа и один выход, что соответствует следующей схеме:

gsconnectingmodels_05.png

In [ ]: 
[P C; P C]
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
Input 1 to output 1
    1.0s + 0.5
-------------------
1.0s^2 + 2.0s + 1.0

Input 1 to output 2
    1.0s + 0.5
-------------------
1.0s^2 + 2.0s + 1.0

Input 2 to output 1
1.0s + 0.5
----------
   1.0s

Input 2 to output 2
1.0s + 0.5
----------
   1.0s

Continuous-time transfer function model

Результирующая модель имеет два входа и два выхода и соответствует структурной схеме:
gsconnectingmodels_07.png

Пространство состояний

Перейдем к следующему способу создания систем - описание в пространстве состояний.

Мы можем конвертировать нашу передаточную функцию в объект пространства состояний. Для этого передадим передаточную функцию P в функцию ss(). Таким образом мы получаем объект StateSpace.

In [ ]: 
P = ss(P)       #   конвертация в вид пространства состояний
Out[0]:
StateSpace{Continuous, Float64}
A = 
  0.0   1.0
 -1.0  -2.0
B = 
 0.0
 1.0
C = 
 0.5  1.0
D = 
 0.0

Continuous-time state-space model

Создать новую систему можно используя матрицы A, B, C, D, описывающие состояния нашей системы.

In [ ]: 
A = [0 1; -1 -2]
B = [0; 1]
C = [1 0]
D = 0

P = ss(A,B,C,D)         #   описание пространства состояний с помощью матриц
Out[0]:
StateSpace{Continuous, Int64}
A = 
  0   1
 -1  -2
B = 
 0
 1
C = 
 1  0
D = 
 0

Continuous-time state-space model

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

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

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

In [ ]: 
import Pkg; Pkg.add("LinearAlgebra")
In [ ]: 
using LinearAlgebra
In [ ]: 
P = ss(A,B,I,D)
Out[0]:
StateSpace{Continuous, Int64}
A = 
  0   1
 -1  -2
B = 
 0
 1
C = 
 1  0
 0  1
D = 
 0
 0

Continuous-time state-space model

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

In [ ]: 
P.A                     #   доступ к матрице А
Out[0]:
2×2 Matrix{Int64}:
  0   1
 -1  -2
In [ ]: 
(; A, B, C, D) = P      #   деструктуризация системы по матрицам
Out[0]:
StateSpace{Continuous, Int64}
A = 
  0   1
 -1  -2
B = 
 0
 1
C = 
 1  0
 0  1
D = 
 0
 0

Continuous-time state-space model

Аналогично мы можем узнать число состояний, входов и выходов системы. Например, обращаясь P.nx мы получим число входов.

In [ ]: 
P.nx            #   число состояний
P.nu            #   число входов
P.ny            #   число выходов
Out[0]:
2

Дискретные системы

In [ ]: 
Ts = 0.1;
Pd1 = tf([1], [1, 2], Ts)
Out[0]:
TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Int64}}
  1
-----
z + 2

Sample Time: 0.1 (seconds)
Discrete-time transfer function model
In [ ]: 
Pd2 = ss(A, B, C, D, Ts)    # получение дискретной системы из матриц с периодом дискретизации Ts
Out[0]:
StateSpace{Discrete{Float64}, Int64}
A = 
  0   1
 -1  -2
B = 
 0
 1
C = 
 1  0
 0  1
D = 
 0
 0

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

Функция c2d() позволяет создать из непрерывной системы дискретную. В качестве аргументов передаем непрерывную систему и период дискретизации.

In [ ]: 
C = pid(1, 2)
Cd = c2d(C, 0.1)             #   конвертация непрерывной ПФ в дискретную
Out[0]:
TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Float64}}
1.0z - 0.95
-----------
1.0z - 1.0

Sample Time: 0.1 (seconds)
Discrete-time transfer function model
In [ ]: 
Cd.Ts                 #   доступ к периоду дискретизации
Out[0]:
0.1

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

  • дискретна ли эта система?
  • устойчива ли она?
In [ ]: 
isdiscrete(Cd)              #   проверка: дискретная ли система
Out[0]:
true
In [ ]: 
isstable(Cd) #   проверка: система устойчива
Out[0]:
false

Анализ систем

При проектировании систем управления нам важно анализировать систему. Начнем с временных характеристик.

Во временной области

In [ ]: 
sys = tf([8, 18, 32],[1, 6, 14, 24])
Out[0]:
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Int64}}
   8s^2 + 18s + 32
---------------------
s^3 + 6s^2 + 14s + 24

Continuous-time transfer function model

В библиотеке ControlSystems.jl есть функции для построения откликов на стандартные входные сигналы. Например, такие как step() - ступенчатое воздействие и impulse() - импульсное воздействие. Эти функции возвращают экземпляр SimResult. Далее его можно использовать, чтобы построить переходную и импульсную характеристики.

In [ ]: 
step(sys, 10)
Out[0]:
ControlSystemsBase.SimResult{Matrix{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Matrix{Float64}, Matrix{Float64}, StateSpace{Continuous, Float64}}([0.0 0.16153638760006078 … 1.3333102517356628 1.3333091364345453], 0.0:0.021:9.996, [0.0 4.785984063879122e-5 … 1.3333866425850898 1.3333877482064391; 0.0 0.0033826376850362256 … 2.858422567558489e-5 2.4084302134424855e-5; 0.0 0.07884153018187812 … -5.427405165606136e-5 -5.2853305897493374e-5], [1.0 1.0 … 1.0 1.0], StateSpace{Continuous, Float64}
A = 
  0.0   2.0   0.0
  0.0   0.0   4.0
 -3.0  -3.5  -6.0
B = 
 0.0
 0.0
 4.0
C = 
 1.0  1.125  2.0
D = 
 0.0

Continuous-time state-space model)
In [ ]: 
impulse(sys, 10)
Out[0]:
ControlSystemsBase.SimResult{Matrix{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Matrix{Float64}, Matrix{Float64}, StateSpace{Continuous, Float64}}([8.0 7.391488176280441 … -5.57212515463793e-5 -5.0510955124734026e-5], 0.0:0.021:9.996, [0.0 0.0067652754108830345 … 5.716845141929535e-5 4.816860433415205e-5; 0.0 0.31536612065294245 … -0.00021709620669415274 -0.00021141322366153976; 4.0 3.5149680075674987 … 6.56717647826236e-5 6.958015858017308e-5], [0.0 0.0 … 0.0 0.0], StateSpace{Continuous, Float64}
A = 
  0.0   2.0   0.0
  0.0   0.0   4.0
 -3.0  -3.5  -6.0
B = 
 0.0
 0.0
 4.0
C = 
 1.0  1.125  2.0
D = 
 0.0

Continuous-time state-space model)
In [ ]: 
using Plots
In [ ]: 
plot(
    plot(step(sys, 10), title="Реакция на ступенчатое воздействие"),
    plot(impulse(sys, 10), title="Реакция на импульс"),
    layout=(2,1)
)
Out[0]:

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

In [ ]: 
A = [-0.8 3.6 -2.1;-3 -1.2 4.8;3 -4.3 -1.1];
B = [0; -1.1; -0.2];
C = [1.2 0 0.6];
D = -0.6;
G = ss(A,B,C,D)

t1 = 0:0.001:18;
u1 = zeros(1,length(t1))
x_0 = [-1.0; 0.0; 2.0]  
plot(lsim(G, u1, t1, x0 = x_0))
Out[0]:

Получить характеристики переходного процесса можно с помощью stepinfo(). На вход функция принимает объект SimResult, а на выход возвращает список параметров, которые помогают оценить перерегулирование системы, время установившегося процесса и тд. Эти данные можно отобразить и на графике, используя функцию plot().

In [ ]: 
h = stepinfo(step(sys, 10))
Out[0]:
StepInfo:
Initial value:     0.000
Final value:       1.333
Step size:         1.333
Peak:              1.687
Peak time:         0.609 s
Overshoot:         26.55 %
Undershoot:         0.00 %
Settling time:     3.507 s
Rise time:         0.210 s
In [ ]: 
plot(h)
Out[0]:

В частотной области

Анализ в частотной области является ключом к пониманию устойчивости и эксплуатационных характеристик систем управления.

Получить амплитудную и фазовую характеристики системы можно используя функцию bode(). А функция margin() возвращает запасы по амплитуде и по фазе, а также значения частот для них. Полученные данные можно отразить на графике используя встроенные функции bodeplot(), marginplot().

In [ ]: 
bode(sys)
Out[0]:
([1.3331771710590545;;; 1.3331661744267964;;; 1.3331544215579518;;; … ;;; 0.08569751877696492;;; 0.08277727446912024;;; 0.079956316147604], [-0.11688671654325213;;; -0.12083092230236496;;; -0.12489421682435445;;; … ;;; -87.69815518360896;;; -87.77660830569647;;; -87.8523925997141], [0.1, 0.10353218432956621, 0.10718913192051278, 0.11097524964120718, 0.11489510001873092, 0.11895340673703195, 0.12315506032928256, 0.12750512407130132, 0.1320088400831418, 0.13667163564620066  …  73.16807143427197, 75.75250258771912, 78.42822061337682, 81.19844993184013, 84.06652885618325, 87.03591361485161, 90.11018251665018, 93.29304026284686, 96.58832241158703, 100.0])
In [ ]: 
margin(sys)
Out[0]:
(wgm = [NaN;;], gm = [Inf;;], wpm = [7.375411262505239;;], pm = [117.19782190490633;;])
In [ ]: 
bodeplot(sys, label="sys",title="Диаграмма Боде")
Out[0]:
In [ ]: 
marginplot(sys)
Out[0]:

Функции nicholsplot() и nyquistplot() также помогут определить устойчивость системы. Отображение графиков можно настраивать. Так, например, на диаграмму Найквиста можно добавить единичную окружность, критическую точку (-1).

In [ ]: 
nyquistplot(
    sys, 
    unit_circle = true,
    Mt_circles = false,
    hz = true,
    label="sys",
    title="Диаграмма Найквиста"
)
Out[0]:
In [ ]: 
nicholsplot(sys, sat = 0, val = 1)
xlims!(-180, 180)
Out[0]:

Чтобы оценить устойчивость системы, полезно знать ее расположение полюсов и нулей. Значение полюсов и нулей можно получить с помощью функции zpkdata().

In [ ]: 
G = -(2s+1)/(s^2+3s+2);
k = 0.7;
T = feedback(G*k,1);

zpkdata(T)
Out[0]:
(Vector{ComplexF64}[[-0.5 + 0.0im];;], Vector{ComplexF64}[[-0.8 + 0.8124038404635958im, -0.8 - 0.8124038404635958im];;], [-1.4;;])

Отразить графически их нахождение на плоскости можно используя pzmap(). Полюса отмечены крестиками, а нули кружочками.

In [ ]: 
G = -(2s+1)/(s^2+3s+2);
k = 0.7;
T = feedback(G*k,1);

pzmap(T; hz = true)
Out[0]:

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

In [ ]: 
plot(rlocus(G,3))
Out[0]:
Группа видимости
ПубличноПросматривать контент могут все пользователи сообщества

Тип

Проекты

Краткое описание

Обсудим прикладные пакеты языка Julia. Использование сторонних языков в скриптах Julia.

Связанные материалы

Зимняя школа Julia - День 4

Теги

Категории

МероприятияЗимняя школа Julia 2025

Языки

    C/C++
    Julia
    MATLAB
    Python
    Другой

Форматы

    engee
    ngscript

Уровень

    Средний

Статус

Опубликовано

Обновлено

Источник

Сообщество