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

Примеры

Проектирование ЛКР

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

где  — вектор состояния, а  — входной вектор.

В приведенном ниже примере выполняется простое проектирование ЛКР для двойного интегратора в дискретном времени с использованием функции lqr. В данном примере мы будем использовать метод lsim, принимающий в качестве входа функцию . Это позволяет легко моделировать систему как с входным сигналом управления, так и с входным сигналом возмущения. Сведения о более сложном проектировании ЛКР и ЛКГ см. в разделе с описанием типа LQGProblem в RobustAndOptimalControl.

using ControlSystemsBase
using LinearAlgebra # Для матрицы тождественности I
using Plots

# Создание системы
Ts      = 0.1
A       = [1 Ts; 0 1]
B       = [0; 1]
C       = [1 0]
sys     = ss(A,B,C,0,Ts)

# Проектирование регулятора
Q       = I # Весовая матрица для состояния
R       = I # Весовая матрица для входа
L       = lqr(Discrete,A,B,Q,R) # Также можно использовать lqr(sys,Q,R)

# Моделирование
u(x,t)  = -L*x .+ 1.5(t>=2.5) # Формирование закона управления (u — функция для t и x), постоянное входное возмущение воздействует на систему с t≧2.5
t       = 0:Ts:5              # Вектор времени
x0      = [1,0]               # Начальное условие
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]")

save_docs_plot("lqrplot.svg"); # скрыть
lqrplot

Для проектирования ЛКГ-регулятора (ЛКГ с фильтром Калмана) см. описание функций

См. также следующий обучающий видеоролик по проектированию ЛКР и ЛКГ.

Функции проектирования ПИД-регулятора

Базовый ПИД-регулятор можно построить с помощью конструктора pid. В ControlSystems.jl мы часто ссылаемся на три различные формулировки ПИД-регулятора, которые определяются следующим образом.

  • Стандартная форма:

  • Последовательная форма:

  • Параллельная форма:

Большинство функций для построения ПИД-регуляторов позволяют выбирать нужную для использования форму.

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

В следующих примерах показаны основные рабочие процессы для проектирования ПИ-/ПИД-регуляторов.

Пример формирования контура ПИ

При построении графика Большой четверки при единичной обратной связи для процесса

using ControlSystemsBase, Plots
P = tf(1, [1,1])^4
gangoffourplot(P, tf(1))

можно заметить, что в районе частот ω = 0,8 рад/с функция чувствительности несколько завышена. Поскольку мы хотим управлять процессом с помощью простого ПИ-регулятора, используем функцию loopshapingPI и указываем, что требуется запас устойчивости по фазе в 60 градусов на данной частоте. График результирующей Большой четверки формируется как для построенного регулятора, так и для единичной обратной связи.

using ControlSystemsBase, Plots
P = tf(1, [1,1])^4
ωp = 0.8
C,kp,ki,fig = loopshapingPI(P,ωp,phasemargin=60,form=:parallel, doplot=true)
fig

Можно также рассмотреть ситуацию, когда нужно создать замкнутую систему с полосой пропускания ω = 2 рад/с. В этом случае мы напишем что-то вроде следующего:

ωp = 2
C60,kp,ki,fig = loopshapingPI(P,ωp,rl=1,phasemargin=60,form=:standard,doplot=true)
fig

Здесь мы указываем, что нужно, чтобы кривая Найквиста L(iω) = P(iω)C(iω) проходила через точку |L(iω)| = rl = 1, arg(L(iω)) = -180 + phasemargin = -180 + 60. По данным Большой четверки видно, что с помощью этого метода проектирования действительно можно получить очень надежный и быстрый регулятор, но удвоение полосы пропускания всех четырех полюсов будет стоить значительного управляющего воздействия.

Пример формирования контура ПИД

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

using ControlSystemsBase, Plots
P  = tf(1, [1,0,0]) # Двойной интегратор
Mt = 1.3            # Максимальная амплитуда комплементарной чувствительности
ϕt = 75             # Угол в точке касания
ω  = 1              # Частота, при которой сохраняется спецификация
C, kp, ki, kd, fig = loopshapingPID(P, ω; Mt, ϕt, doplot=true)
fig

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

Поскольку мы проектируем ПИД-регулятор, мы ожидаем большого коэффициента усиления регулятора для высоких частот. Как правило, это нежелательно как по соображениям надежности, так и по причинам, связанным с шумами. Обычно этот вопрос решается последовательным добавлением фильтра нижних частот к регулятору. В приведенном ниже примере передача именованного аргумента Tf=1/20ω указывает на то, что требуется добавить фильтр нижних частот второго порядка с частотой среза, в 20 раз превышающей расчетную частоту.

Tf = 1/20ω
C, kp, ki, kd, fig, CF = loopshapingPID(P, ω; Mt, ϕt, doplot=true, Tf)
fig

Как видно, добавление фильтра увеличивает спад на высокой частоте как в , так и в , что, как правило, нам и необходимо.

Чтобы лучше контролировать фильтр, его можно предварительно спроектировать и передать в метод loopshapingPID с именованным аргументом F:

F = tf(1, [Tf^2, 2*Tf/sqrt(2), 1])
C, kp, ki, kd, fig, CF = loopshapingPID(P, ω; Mt, ϕt, doplot=true, F)

Размещение полюсов и нулей повышенной сложности

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

В следующем примере показано, как можно выполнить размещение полюсов и нулей повышенной сложности с помощью функции rstc (rstd в дискретном времени). Задача состоит в том, чтобы немного ускорить процесс и погасить плохо затухающие полюса.

Определим процесс

ζ = 0.2
ω = 1

B = [1]
A = [1, 2ζ*ω, ω^2]
P = tf(B,A)

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

using ControlSystems
import DSP: conv
# Схема выполнения
ζ0 = 0.7
ω0 = 2
Am = [1, 2ζ0*ω0, ω0^2]
Ao = conv(2Am, [1/2, 1]) # Многочлен-наблюдатель, добавление дополнительного полюса из-за интегратора
AR = [1,0] # Назначение регулятору интегратора ( 1/(s+0) )

B⁺  = [1] # Многочлен-числитель процесса можно представить в виде = B⁺B⁻, где B⁻ содержит нули, которые не нужно отменять (неминимальная фаза и плохо затухающие нули)
B⁻  = [1]
Bm  = conv(B⁺, B⁻) # В этом случае следует сохранить весь многочлен-числитель процесса

R,S,T = rstc(B⁺,B⁻,A,Bm,Am,Ao,AR) # Вычисление многочленов-регуляторов с двумя степенями свободы

Gcl = tf(conv(B,T),zpconv(A,R,B,S)) # Формирование многочлена замкнутого контура от эталонного сигнала к выходному, характеристический многочлен замкнутого контура равен AR + BS, функция zpconv выполняет перемножение многочленов и следит за тем, чтобы векторы коэффициентов были одинаковой длины

plot(step(P, 20))
plot!(step(Gcl, 20)) # Визуализация откликов открытого и замкнутого контуров.
save_docs_plot("ppstepplot.svg") # hide
gangoffourplot(P, tf(-S,R)) # Построение графика Большой четверки, чтобы убедиться, что все передаточные функции работают
save_docs_plot("ppgofplot.svg"); # hide

ppstepplot ppgofplot

Границы устойчивости для ПИД-регуляторов

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

P1 = s -> exp(-sqrt(s))
doplot = true
form = :parallel
kp, ki, f1 = stabregionPID(P1,exp10.(range(-5, stop=1, length=1000)); doplot, form); f1
P2 = s -> 100*(s+6).^2. /(s.*(s+1).^2. *(s+50).^2)
kp, ki, f2 = stabregionPID(P2,exp10.(range(-5, stop=2, length=1000)); doplot, form); f2
P3 = tf(1,[1,1])^4
kp, ki, f3 = stabregionPID(P3,exp10.(range(-5, stop=0, length=1000)); doplot, form); f3

save_docs_plot(f1, "stab1.svg") # hide
save_docs_plot(f2, "stab2.svg") # hide
save_docs_plot(f3, "stab3.svg"); # hide

stab1 stab2 stab3

Графики ПИД

В данном примере используется функция pidplots, которая принимает векторы ПИД-параметров и строит соответствующие графики. Задача состоит в том, чтобы взять систему с полосой пропускания 1 рад/с и создать замкнутую систему с полосой пропускания 0,1 рад/с. Если не проявить осторожность и приступить к размещению полюсов, можно легко получить систему с очень низкой надежностью.

using ControlSystemsBase
P = tf([1.], [1., 1])

ζ = 0.5 # Желаемое демпфирование
ws = exp10.(range(-1, stop=2, length=8)) # Вектор полос пропускания замкнутых контуров
kp = 2*ζ*ws .- 1 # Простое размещение полюсов с ПИ при заданной полосе пропускания замкнутого контура, полюса размещаются по схеме Баттерворта
ki = ws.^2

ω = exp10.(range(-3, stop = 2, length = 500))
pidplots(
    P,
    :nyquist;
    params_p = kp,
    params_i = ki,
    ω = ω,
    ylims = (-2, 2),
    xlims = (-3, 3),
    form = :parallel,
)
save_docs_plot("pidplotsnyquist1.svg") # hide
pidplots(P, :gof; params_p = kp, params_i = ki, ω = ω, legend = false, form=:parallel, legendfontsize=6, size=(1000, 1000))
# Вы также можете запросить графики Найквиста и Большой четверки (доступны дополнительные графики, см. ?pidplots):
# pidplots(P,:nyquist,:gof;kps=kp,kis=ki,ω=ω);
save_docs_plot("pidplotsgof1.svg"); # hide

pidplotsnyquist1 pidplotsgof1

Теперь попробуем другую стратегию, в которой мы задали частоту единичного усиления в 0,1 рад/с

kp = range(-1, stop=1, length=8) #
ki = sqrt.(1 .- kp.^2)/10

pidplots(P,:nyquist,;params_p=kp,params_i=ki,ylims=(-1,1),xlims=(-1.5,1.5), form=:parallel)
save_docs_plot("pidplotsnyquist2.svg") # hide
pidplots(P,:gof,;params_p=kp,params_i=ki,legend=false,ylims=(0.08,8),xlims=(0.003,20), form=:parallel, legendfontsize=6, size=(1000, 1000))
save_docs_plot("pidplotsgof2.svg"); # hide

pidplotsnyquist2 pidplotsgof2

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