Анализ временных и частотных характеристик

Частотная характеристика

Частотные характеристики рассчитываются с помощью методов freqresp, bode, sigma и nyquist. Для получения графиков частотных характеристик служат функции bodeplot, sigmaplot, nyquistplot, marginplot и nicholsplot.

Значение любой функции TransferFunction в определенной точке можно вычислить с помощью вызовов F(s), F(omega, true) и F(z, false).

  • F(s) вычисляет передаточную функцию с непрерывным временем F в s.

  • F(omega,true) вычисляет передаточную функцию с непрерывным временем F в exp(i*Ts*omega).

  • F(z,false) вычисляет передаточную функцию с непрерывным временем F в z.

Ниже приводится видео, в котором демонстрируется анализ частотных характеристик в ControlSystems.jl.

Временная характеристика (моделирование)

За моделирование с произвольными входами в первую очередь отвечает функция lsim, а step и impulse служат вспомогательными функциями для моделирования реакции на определенные входные сигналы.

Функция lsim может принимать входной вектор u, содержащий дискретизированную входную траекторию, или входную функцию, принимающую состояние и время в качестве аргументов: u(x,t). С помощью этой функции можно, например, легко моделировать линейно нарастающие характеристики, управление с насыщенной обратной связью по состоянию и т. д. Дополнительные сведения см. в docstring метода lsim.

Сведения о расширенных возможностях нелинейного моделирования см. в примечаниях к ModelingToolkit и DifferentialEquations в разделе The wider Julia ecosystem for control.

Пример ступенчатой характеристики:

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

using ControlSystemsBase, Plots
G = tf(1, [1, 1, 1])
res = step(G, 20) # Моделируем ступенчатую характеристику с интервалом 20 секунд
plot(res)

С помощью функции stepinfo можно вычислить свойства ступенчатой характеристики:

si = stepinfo(res)
StepInfo:
Initial value:     0.000
Final value:       1.000
Step size:         1.000
Peak:              1.163
Peak time:         3.652 s
Overshoot:         16.30 %
Undershoot:         0.00 %
Settling time:     8.134 s
Rise time:         1.660 s

Можно также построить график объекта StepInfo.

plot(si)

Пример lsim:

Функция lsim может принимать управляющий входной сигнал в одной из следующих форм:

  1. Массив равномерно дискретизированных значений; в этом случае предполагается, что аргумент u имеет форму nu × n_time;

  2. Функция от состояния и времени u(x,t). Данная форма позволяет производить моделирование обратной связи по состоянию со ступенчатой характеристикой в момент времени : u(x, t) = amplitude * (t > t0), с линейно нарастающей характеристикой: u(x, t) = t и т. д.

В примере ниже моделируется обратная связь по состоянию со ступенчатым возмущением в момент времени путем передачи функции u(x,t) = -L*x .+ (t > 4) в lsim:

using ControlSystems
using LinearAlgebra: I
using Plots

A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)

u(x,t) = -L*x .+ (t > 4) # Обратная связь по состоянию со ступенчатым возмущением
t  = 0:0.1:12
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]"); vline!([4], lab="Step disturbance", l=(:black, :dash, 0.5))

Ниже приводится видео, в котором демонстрируется моделирование во временной области в ControlSystems.jl.

Строки docstring

# ControlSystemsBase.lsimMethod

res = lsim(sys::DelayLtiSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), abstol=1e-6, reltol=1e-6, force_dtmin=true, kwargs...)

Моделирует систему sys в течение времени t с использованием входного сигнала u, с начальным состоянием x0 и с помощью метода alg.

Аргументы

t: должен быть AbstractVector с равноудаленными временными выборками (постоянная t[i] - t[i-1]) u: функция определения управляющего сигнала ut в момент времени t в любой из следующих форм.

  • u: функция определения управляющего сигнала uₜ в момент времени t в любой из следующих форм.

    • Константа Number или Vector, интерпретируемая как входная константа.

    • Функция u(x, t), принимающая внутреннее состояние и время. Заметим, что представление состояния для систем с запаздыванием не такое, как для рациональных систем.

    • Функция на месте u(uₜ, x, t). (Чуть более эффективна.)

alg, abstol, reltol и kwargs...: передаются в DelayDiffEq.solve.

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

Возвращает экземпляр SimResult, который можно построить напрямую или деструктурировать в y, t, x, u = res.

# ControlSystemsBase.lsimMethod

lsim(sys::HammersteinWienerSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=Tsit5(), abstol=1e-6, reltol=1e-6, kwargs...)

Моделирует систему sys в течение времени t с использованием входного сигнала u, с начальным состоянием x0 и с помощью метода alg.

Аргументы

  • t: должен быть AbstractVector с равноудаленными временными выборками (постоянная t[i] - t[i-1])

  • u: функция определения управляющего сигнала uₜ в момент времени t в любой из следующих форм.

    • Константа Number или Vector, интерпретируемая как uₜ .= u, или

    функция uₜ .= u(x, t) или функция на месте u(uₜ, x, t). (Чуть более эффективна.)

  • alg, abstol, reltol и kwargs...: передаются в OrdinaryDiffEq.solve.

Возвращает экземпляр SimResult.

# ControlSystems.SimulatorType

Simulator

Поля

P::StateSpace
f = (x,p,t) -> x
y = (x,t)   -> y

# ControlSystems.SimulatorMethod

Simulator(P::StateSpace, u = (x,t) -> 0)

Используется для моделирования систем с непрерывным временем. Для получения дополнительных сведений см. описание функции ?solve.

Использование:

using OrdinaryDiffEq, Plots
dt             = 0.1
tfinal         = 20
t              = 0:dt:tfinal
P              = ss(tf(1,[2,1])^2)
K              = 5
reference(x,t) = [1.]
s              = Simulator(P, reference)
x0             = [0.,0]
tspan          = (0.0,tfinal)
sol            = solve(s, x0, tspan, Tsit5())
plot(t, s.y(sol, t)[:], lab="Open loop step response")

# Base.stepMethod

y, t, x = step(sys[, tfinal])
y, t, x = step(sys[, t])

Вычисляет отклик системы sys на единичный шаг в момент времени t = 0. Если конечное время tfinal или вектор времени t не указаны, вычисление выполняется на основе расположения полюсов системы.

Возвращаемое значение — структура типа SimResult. График SimResul можно построить с помощью plot(result) или деструктурировать как y, t, x = result.

y имеет размер (ny, length(t), nu), x имеет размер (nx, length(t), nu).

См. также описание stepinfo и lsim.

# ControlSystemsBase.impulseMethod

y, t, x = impulse(sys[, tfinal])
y, t, x = impulse(sys[, t])

Вычисляет отклик системы sys на импульс в момент времени t = 0. Для систем с непрерывным временем импульс является единичным импульсом Дирака. Для систем с дискретным временем импульс длится одну выборку и имеет амплитуду 1/Ts. Если конечное время tfinal или вектор времени t не указаны, вычисление выполняется на основе расположения полюсов системы.

Возвращаемое значение — структура типа SimResult. График SimResul можно построить с помощью plot(result) или деструктурировать как y, t, x = result.

y имеет размер (ny, length(t), nu), x имеет размер (nx, length(t), nu).

См. также описание lsim.

# ControlSystemsBase.lsim!Method

res = lsim!(ws::LsimWorkspace, sys::AbstractStateSpace{<:Discrete}, u, [t]; x0)

Версия lsim на месте, которая принимает объект рабочей области, созданный при вызове LsimWorkspace. Обратите внимание, что если u является функцией, res.u === ws.u. Если u является массивом, res.u === u.

# ControlSystemsBase.lsimMethod

result = lsim(sys, u[, t]; x0, method])
result = lsim(sys, u::Function, t; x0, method)

Вычисляет временной отклик системы sys на входной сигнал u. Если x0 пропущен, используется нулевой вектор.

Результирующая структура результата содержит поля y, t, x, u и может быть деструктурирована автоматически с помощью итерации, например

y, t, x, u = result

график result::SimResult также можно построить напрямую:

plot(result, plotu=true, plotx=false)

y, x, u имеют время во втором измерении. По умолчанию начальное состояние x0 равно нулю.

Моделирование систем с непрерывным временем осуществляется с помощью решателя ODE, если u является функцией (требуется использование ControlSystems). Если u представляет собой массив, перед моделированием система дискретизируется (по умолчанию с помощью method=:zoh). Сведения об интерфейсе никого уровня см. в описании ?Simulator и ?solve. Для систем с непрерывным временем именованные аргументы передаются решателю ODE. По умолчанию используется параметр dtmax = t[2]-t[1], не позволяющий решателю перешагивать через разрывы в u(x, t). Так решатель не может делать слишком большие шаги, но может замедлить моделирование при плавности u. Для отключения этого поведения следует задать dtmax = Inf.

u может быть функцией или матрицей предварительно вычисленных управляющих сигналов и должна иметь измерения (nu, length(t)). Если u — функция, то u(x,i) (для дискретных систем) или u(x,t) (для непрерывных) вызывается для вычисления управляющего сигнала при каждой итерации (момент времени, используемый решателем). Это можно использовать для установки закона управления, например обратной связи по состоянию u(x,t) = -L*x, вычисляемой с помощью lqr. Для моделирования единичного шага при t=t₀ используйте (x,t)-> t ≥ t₀, для линейного изменения — (x,t)-> t, для шага при t=5 — (x,t)-> (t >= 5) и т. д.

Примечание. Функция u будет вызвана один раз перед моделированием, чтобы проверить, что она возвращает массив правильных измерений. Это может вызвать проблемы, если u отслеживает состояние. Чтобы отключить эту проверку, передайте check_u = false.

Для достижения максимальной производительности см. описание функции lsim!, доступной только для систем с дискретным временем.

Пример использования:

using ControlSystems
using LinearAlgebra: I
using Plots

A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)

u(x,t) = -L*x # Строим закон управления
t  = 0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]")

# Альтернативный способ построения графика
res = lsim(sys,u,t,x0=x0)
plot(res)

# ControlSystemsBase.stepinfoMethod

stepinfo(res::SimResult; y0 = nothing, yf = nothing, settling_th = 0.02, risetime_th = (0.1, 0.9))

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

  • y0: начальное значение отклика.

  • yf: конечное значение отклика.

  • stepsize: размер шага.

  • peak: пиковое значение отклика.

  • peaktime: время возникновения пика.

  • overshoot: процентное превышение отклика.

  • undershoot: процентное недостижение отклика. Если ступенчатый отклик никогда не опускается ниже начального значения, недостижение равно нулю.

  • settlingtime: время, в которое отклик установился в пределах settling_th от конечного значения.

  • settlingtimeind: индекс, в котором отклик установился в пределах settling_th от конечного значения.

  • risetime: время, в которое отклик вырос с risetime_th[1] до risetime_th[2] от конечного значения.

Аргументы

  • res: результат моделирования с использованием step (или lsim).

  • y0: начальное значение. Если оно не указано, используется первое значение из отклика.

  • yf: конечное значение. Если оно не указано, используется первое значение из отклика. Чтобы автоматически вычисленное значение имело смысл, моделирование должно достичь установившегося состояния. В противном случае конечное значение можно указать вручную.

  • settling_th: пороговое значение для вычисления времени установления. Время установления — это время, в которое отклик установился в пределах settling_th от конечного значения.

  • risetime_th: нижний и верхний порог для вычисления времени нарастания. Время нарастания — это время, в которое отклик вырос с risetime_th[1] до risetime_th[2] от конечного значения.

Пример:

G = tf([1], [1, 1, 1])
res = step(G, 15)
si = stepinfo(res)
plot(si)

# ControlSystemsBase.LsimWorkspaceMethod

LsimWorkspace(sys::AbstractStateSpace, N::Int)
LsimWorkspace(sys::AbstractStateSpace, u::AbstractMatrix)
LsimWorkspace{T}(ny, nu, nx, N)

Генерирует объект рабочей области для использования с функцией на месте lsim!. sys — система с дискретным временем, подлежащая моделированию, N — количество временных шагов. В качестве альтернативы вместо N можно указать входной u. Примечание. Для многопоточных приложений следует создавать по одному объекту рабочей области на поток.

# ControlSystemsBase.StepInfoType

StepInfo

Вычисляется с помощью stepinfo.

Поля

  • y0: начальное значение ступенчатого отклика.

  • yf: конечное значение ступенчатого отклика.

  • stepsize: размер шага.

  • peak: пиковое значение ступенчатого отклика.

  • peaktime: время возникновения пика.

  • overshoot: превышение ступенчатого отклика.

  • settlingtime: время, в которое ступенчатый отклик установился в пределах settling_th от конечного значения.

  • settlingtimeind::Int: индекс, при котором ступенчатый отклик установился в пределах settling_th от конечного значения.

  • risetime: время, в которое отклик вырос с risetime_th[1] до risetime_th[2] от конечного значения.

  • i10::Int: индекс, при котором отклик достигает risetime_th[1].

  • i90::Int: индекс, при котором отклик достигает risetime_th[2].

  • res::SimResult{SR}: результат моделирования, используемый для расчета характеристик ступенчатого отклика.

  • settling_th: пороговое значение, используемое для вычисления settlingtime и settlingtimeind.

  • risetime_th: пороговые значения, используемые для вычисления risetime, i10 и i90.

# ControlSystemsBase.bodeMethod

mag, phase, w = bode(sys[, w]; unwrap=true)

Вычисляет амплитудную и фазовую части частотной характеристики системы sys на частотах w. См. также описание bodeplot.

mag и phase имеют размер (ny, nu, length(w)). Если unwrap задано значение true (по умолчанию), то к фазовым углам будет применена функция unwrap!. Эта процедура является дорогостоящей, и от нее можно отказаться, если развертывание не требуется.

Для повышения производительности см. описание функции bodemag!, которая вычисляет только амплитуду.

# ControlSystemsBase.bodemag!Method

mag = bodemag!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector)

Вычисляет амплитуду Боде, действующую на месте, для экземпляра BodemagWorkspace. Обратите внимание, что возвращаемому массиву амплитуд присваивается псевдоним ws.mag. Выходным массивом mag является ∈ xD835__xDC11(ny, nu, nω).

# ControlSystemsBase.evalfrMethod

evalfr(sys, x)

Вычисляет передаточную функцию LTI-системы sys при комплексном числе s=x (непрерывное время) или z=x (дискретное время).

Для многих значений x следует использовать freqresp.

# ControlSystemsBase.freqresp!Method

freqresp!(R::Array{T, 3}, sys::LTISystem, w_vec::AbstractVector{<:Real})

Версия freqresp на месте, которая принимает предварительно выделенный массив R размером (ny, nu, nw)`.

# ControlSystemsBase.freqrespMethod

sys_fr = freqresp(sys, w)

Вычисляет частотную характеристику линейной системы

w -> C*((iw*im*I - A)^-1)*B + D

системы sys по вектору частоты w.

# ControlSystemsBase.nyquistMethod

re, im, w = nyquist(sys[, w])

Вычисляет действительную и мнимую части частотной характеристики системы sys на частотах w. См. также описание nyquistplot.

re и im имеют размер (ny, nu, length(w)).

# ControlSystemsBase.sigmaMethod

sv, w = sigma(sys[, w])

Вычисляет сингулярные значения sv частотной характеристики для систем sys на частотах w. См. также описание sigmaplot.

sv имеет размер (min(ny, nu), length(w)).

# ControlSystemsBase.BodemagWorkspaceMethod

BodemagWorkspace(sys::LTISystem, N::Int)
BodemagWorkspace(sys::LTISystem, ω::AbstractVector)
BodemagWorkspace(R::Array{Complex{T}, 3}, mag::Array{T, 3})
BodemagWorkspace{T}(ny, nu, N)

Генерирует объект рабочей области для использования с функцией на месте bodemag!. N — это количество частотных точек, в качестве альтернативы вместо N можно указать входной ω. Примечание. Для многопоточных приложений следует создавать по одному объекту рабочей области на поток.

Аргументы

  • mag: выходной массив ∈ xD835__xDC11(ny, nu, nω).

  • R: массив частотных характеристик ∈ xD835__xDC02(ny, nu, nω)

# ControlSystemsBase.TransferFunctionMethod

F(s), F(omega, true), F(z, false)

Нотация для вычисления частотных характеристик.

  • F(s) вычисляет передаточную функцию с непрерывным временем F в s.

  • F(omega,true) вычисляет передаточную функцию с непрерывным временем F в exp(imTsomega)

  • F(z,false) вычисляет передаточную функцию с непрерывным временем F в z