Функциональность нелинейности

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

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

      ┌─────────┐
 y◄───┤         │◄────u
      │    P    │
Δy┌───┤         │◄───┐Δu
  │   └─────────┘    │
  │                  │
  │      ┌───┐       │
  └─────►│ f ├───────┘
         └───┘

то есть как дробно-линейное преобразование (ДЛП) между линейной системой и диагональной матрицей с помощью скалярных нелинейных функций . Это представление идентично применяемому для систем с запаздыванием и доступно пользователям в аналогичной форме. Основной точкой входа является функция nonlinearity, которая принимает нелинейную функцию : nonlinearity(f). В результате создается примитивная система, содержащая только нелинейность, но ведущая себя как обычная система LTISystem во время алгебраических операций. Мы проиллюстрируем ее использование на ряде примеров.

Примеры

Насыщение управляющего сигнала

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

using ControlSystems, Plots
using ControlSystemsBase: nonlinearity # Эта функциональность не экспортируется из-за статуса бета-версии

C    = pid(1, 0.1, form=:parallel)                  # Стандартный ПИ-регулятор
nl   = nonlinearity(x->clamp(x, -0.7, 0.7)) # насыщающаяся нелинейность
satC = nl*C # Подключаем насыщение на выходе C

ControlSystemsBase.HammersteinWienerSystem{Float64}

P: StateSpace{Continuous, Float64}
A =
 0.0
B =
 0.25  0.0
C =
 0.0
 0.4
D =
 0.0  1.0
 1.0  0.0

Continuous-time state-space model

Nonlinearities: Function[Main.var"#1#2"()]

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

P   = tf(1, [1, 1])    # объект
G   = feedback(P*C)    # замкнутый контур без нелинейности
Gnl = feedback(P*satC) # замкнутый контур с насыщением

Gu   = feedback(C, P)    # замкнутый контур от опорного узла к управляющему сигналу без нелинейности
Gunl = feedback(satC, P) # замкнутый контур от опорного узла к управляющему сигналу с насыщением

plot(step([G; Gu], 5), lab = ["Linear y" "Linear u"])
plot!(step([Gnl; Gunl], 5), lab = ["Nonlinear y" "Nonlinear u"])

Так как насыщающаяся нелинейность часто встречается на практике, мы предоставляем конструктор ControlSystemsBase.saturation, который автоматически создает эквивалент nonlinearity(x->clamp(x, -0.7, 0.7)), в то же время обеспечивая легко идентифицируемое имя функции при выводе системы:

using ControlSystemsBase: saturation
saturation(0.7)

ControlSystemsBase.HammersteinWienerSystem{Float64}

P: StateSpace{Continuous, Float64}
D =
 0.0  1.0
 1.0  0.0

Continuous-time state-space model

Nonlinearities: Function[saturation(0.7)]

См. также описание функции ControlSystemsBase.ratelimit, которая насыщает производную сигнала.

Ненулевая рабочая точка

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

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

Ниже мы демонстрируем совместное использование конструкторов offset и saturation. Рассматриваемая модель — это линеаризованная модель процесса четырех емкостей.

Система линеаризуется вокруг рабочей точки:

xr = [10, 10, 4.9, 4.9] # базовое состояние
ur = [0.263, 0.263]     # управляющий вход в точке управления

и задается следующим образом:

using LinearAlgebra
kc, k1, k2, g = 0.5, 1.6, 1.6, 9.81
A1 = A3 = A2 = A4 = 4.9
a1, a3, a2, a4 = 0.03, 0.03, 0.03, 0.03
h01, h02, h03, h04 = xr
T1, T2 = (A1/a1)sqrt(2*h01/g), (A2/a2)sqrt(2*h02/g)
T3, T4 = (A3/a3)sqrt(2*h03/g), (A4/a4)sqrt(2*h04/g)
c1, c2 = (T1*k1*kc/A1), (T2*k2*kc/A2)
γ1, γ2 = 0.3, 0.3

# Определяем динамику процесса
A = [-1/T1     0 A3/(A1*T3)          0
        0     -1/T2          0 A4/(A2*T4)
        0         0      -1/T3          0
        0         0          0      -1/T4]
B = [γ1*k1/A1     0
        0                γ2*k2/A2
        0                (1-γ2)k2/A3
        (1-γ1)k1/A4 0              ]

C = kc*[I(2) 0*I(2)] # Измеряем уровень в первых двух емкостях
D = 0
G = ss(A,B,C,D)

ПИД-регулятор с фильтром задается следующим образом:

F = tf(1, [0.63, 1.12, 1])
Cpid = pid(0.26, 0.001, 15.9, form=:parallel)*F |> ss

Чтобы сделать регулятор многомерным, мы добавляем статический предкомпенсатор, который разъединяет систему при нулевой частоте.

iG0 = dcgain(G)
iG0 ./= maximum(abs, iG0)
C = (Cpid .* I(2)) * iG0

Насосы (два), обслуживающие емкости, могут только подавать жидкость в емкости, но не высасывать ее. Таким образом, насос насыщается снизу при 0 и сверху при максимальной мощности насоса 0,4.

using ControlSystemsBase: offset
umin = [0.0, 0.0]
umax = [0.4, 0.4]

yr    = G.C*xr  # Опорный выходной сигнал
Gop   = offset(yr) * G * offset(-ur) # Заставляем объект работать в Δ-координатах,
C_sat = saturation(umin, umax) * C   # в то время как регулятор и насыщение работают в исходных координатах

ControlSystemsBase.HammersteinWienerSystem{Float64}

P: StateSpace{Continuous, Float64}
A =
 0.0   0.015625             0.0                0.0   0.0                  0.0
 0.0   0.0                  1.0                0.0   0.0                  0.0
 0.0  -1.5873015873015872  -1.777777777777778  0.0   0.0                  0.0
 0.0   0.0                  0.0                0.0   0.015625             0.0
 0.0   0.0                  0.0                0.0   0.0                  1.0
 0.0   0.0                  0.0                0.0  -1.5873015873015872  -1.777777777777778
B =
 0.0                 0.0                 0.0  0.0
 0.0                 0.0                 0.0  0.0
 1.7142857142857144  4.0                 0.0  0.0
 0.0                 0.0                 0.0  0.0
 0.0                 0.0                 0.0  0.0
 4.0                 1.7142857142857144  0.0  0.0
C =
 0.0                   0.0                  0.0                0.0                   0.0                  0.0
 0.0                   0.0                  0.0                0.0                   0.0                  0.0
 0.025396825396825397  0.10317460317460318  6.309523809523809  0.0                   0.0                  0.0
 0.0                   0.0                  0.0                0.025396825396825397  0.10317460317460318  6.309523809523809
D =
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

Continuous-time state-space model

Nonlinearities: Function[saturation(0.0, 0.4), saturation(0.0, 0.4)]

Теперь мы моделируем замкнутую систему. Начальное состояние объекта корректируется по рабочей точке x0-xr, так как объект работает в Δ-координатах.

x0 = [2, 1, 8, 3] # Начальные уровни в емкостях
plot(
    plot(lsim(feedback(Gop*C_sat), yr, 0:1:3000, x0=[x0-xr; zeros(C.nx)]), layout=1, sp=1, title="Outputs", ylabel=""),
    plot(lsim(feedback(C_sat, Gop), yr, 0:1:3000, x0=[zeros(C.nx); x0-xr]), layout=1, sp=1, title="Control signals", ylabel="")
)
hline!([yr[1]], label="Reference", l=:dash, sp=1, c=1)

Осциллятор Дуффинга

В этом примере мы моделируем нелинейную систему и управляем ею.

Для этого сначала нарисуем блок-схему.

10u    ┌───┐
──────►│+  │   ┌───┐   ┌───┐
 ┌────►│-  │ ẍ │ 1 │ ẋ │ 1 │ x
 │ ┌──►│-  ├──►│ - ├┬─►│ - ├─┬──►
 │ │ ┌►│-  │   │ s ││  │ s │ │
 │ │ │ └───┘   └───┘│  └───┘ │
 │ │ │              │        │
 │ │ │   ┌───┐      │        │
 │ │ └───┤ c │◄─────┘        │
 │ │     └───┘               │
 │ │                         │
 │ │     ┌───┐               │
 │ └─────┤ k │◄──────────────┤
 │       └───┘               │
 │                           │
 │       ┌───┐   ┌───┐       │
 └───────┤ k³│◄──┤ x³│◄──────┘
         └───┘   └───┘

Как мы видим, входной сигнал проходит через внутренний контур скорости, прежде чем попасть на выход . Мы можем построить передаточную функцию этого внутреннего замкнутого контура с помощью feedback(1/s, c), то есть замкнуть контур вокруг интегратора по . Затем этот внутренний контур последовательно соединяется с другим интегратором, и контур обратной связи замыкается с помощью pos_loop_feedback по пути обратной связи. Обратите внимание: для получения правильного усиления на входе итоговая система умножается на 10 справа, так как для нелинейных систем 10*sys и sys*10 — это не всегда одно и то же!

using ControlSystems, Plots
using ControlSystemsBase: nonlinearity
k  = 10
k3 = 2
c  = 1

s = tf("s")

cube = nonlinearity(x->x^3)
vel_loop = feedback(1/s, c)
pos_loop_feedback = (k3*cube + k)
duffing = feedback(vel_loop/s, pos_loop_feedback)*10

plot(step(duffing, 20), title="Duffing oscillator open-loop step response")

Теперь покажем, как с помощью кругового критерия можно подтвердить устойчивость замкнутого контура. Приведенная ниже функция circle_criterion строит кривую Найквиста передаточной функции и определяет исключаемую окружность, находя границы сектора для статической нелинейности . Затем мы выбираем регулятор и проверяем, остается ли он за пределами окружности. Для нахождения границ сектора мы выбираем область, в которой будет вычисляться нелинейность. Функция стремится к бесконечности быстрее, чем любая линейная функция, поэтому верхняя граница сектора равна ∞. Однако если мы ограничим нелинейность более узкой областью, то получим конечную границу сектора:

function circle_criterion(L::ControlSystemsBase.HammersteinWienerSystem, domain::Tuple; N=10000)
    fun = x->L.f[](x)/x
    x = range(domain[1], stop=domain[2], length=N)
    0 ∈ x && (x = filter(!=(0), x)) # Делить на ноль нельзя
    k1, k2 = extrema(fun, x)

    f1 = plot(L.f[], domain[1], domain[2], title="Nonlinearity", lab="f(x)", xlab="x")
    plot!(x, [k1.*x k2.*x], lab=["k1 = $(round(k1, sigdigits=2))" "k2 = $(round(k2, sigdigits=2))"], l=(:dash), legend=:bottomright)

    p1 = -1/k2 # Близко к исходному
    p2 = -1/k1 # Далеко от исходного

    c = (p1 + p2)/2
    r = (p2 - p1)/2

    Lnominal = sminreal(ss(L.A, L.B1, L.C1, L.D11, L.P.timeevol))
    f2 = nyquistplot(Lnominal)
    if p2 < -1000 # Из-за ошибки на графиках
        vspan!([-1000, p1], fillalpha=0.7, c=:red, primary=false)
    else
        th = 0:0.01:2pi
        Cs,Ss = cos.(th), sin.(th)
        plot!(r.*Cs .+ c, r.*Ss, fill=true, fillalpha=0.7, c=:red, primary=false)
    end

    plot(f1,f2)
end


C = pid(2, 0, 1, form=:parallel)*tf(1, [0.01,1])
f1 = circle_criterion(duffing*C, (-1, 1))
plot!(sp=2, ylims=(-10, 3), xlims=(-5, 11))
f2 = plot(step(feedback(duffing, C), 8), plotx=true, plot_title="Controlled oscillator disturbance step response", layout=4)
plot(f1,f2, size=(1300,800))

Так как мы вычислили нелинейность в небольшой области, выход за ее пределы чреват рисками.

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

для получения фактической окружности на плоскости Найквиста.

wiggly = nonlinearity(x->x+sin(x)) # Эта функция менее линейна.
vel_loop = feedback(1/s, c)
pos_loop_feedback = (k3*wiggly + k)
duffing = feedback(vel_loop/s, pos_loop_feedback)*10

C = pid(2, 5, 1, form=:parallel)*tf(1,[0.1, 1])
f1 = circle_criterion(duffing*C, (-2pi, 2pi))
plot!(sp=2, ylims=(-5, 2), xlims=(-2.1, 0.1))
f2 = plot(step(feedback(duffing, C), 8), plotx=true, plot_title="Controlled wiggly oscillator disturbance step response", layout=5)
plot(f1,f2, size=(1300,800))

Ограничения

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

Возможные пути дальнейшего развития

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

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

Для упрощения более сложных задач нелинейного моделирования существуют пакеты ModelingToolkit.jl (MTK) и ModelingToolkitStandardLibrary.jl. В руководствах

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

Строки docstring

# ControlSystemsBase.nonlinearity — Function

nonlinearity(f)
nonlinearity(T, f)

Создает чистую нелинейность. f считается статической (без памяти) нелинейной функцией из .

Типом T по умолчанию является Float64.

ПРИМЕЧАНИЕ. Нелинейная функциональность в ControlSystemsBase.jl в настоящее время является экспериментальной и может быть изменена без соблюдения семантического управления версиями Используйте на свой страх и риск.

Пример:

Создает LTI-систему со статической входной нелинейностью, насыщающей вход до [-1,1].

tf(1, [1, 1])*nonlinearity(x->clamp(x, -1, 1))

См. также описание предварительно определенных нелинейностей saturation, offset.

Примечание. При создании линейных систем с нелинейностями часто важно правильно обращаться с рабочими точками. Сведения об обработке рабочих точек см. в описании ControlSystemsBase.offset.

# ControlSystemsBase.offset — Function

offset(val)

Создает нелинейность с постоянным смещением x -> x + val.

Пример:

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

offset_sys = offset(y₀) * sys * offset(-u₀)

обратите внимание на знак смещения u₀. Это гарантирует, что sys работает в координатах Δu = u-u₀, Δy = y-y₀, а входные и выходные данные для системы смещения находятся в своей системе координат без смещения. Если система линеаризована вокруг x₀, y₀ задается с помощью C*x₀. Дополнительную информацию и пример можно найти здесь: https://juliacontrol.github.io/ControlSystemsBase.jl/latest/lib/nonlinear/#Non-zero-operating-point.

# ControlSystemsBase.saturation — Function

saturation(val)
saturation(lower, upper)

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

Csat = saturation(val) * C

           y▲   ────── upper
            │  /
            │ /
            │/
  ──────────┼────────► u
           /│
		  / │
         /  │
lower────

# ControlSystemsBase.ratelimit — Function

ratelimit(val; Tf)
ratelimit(lower, upper; Tf)

Создает нелинейность, ограничивающую скорость изменения сигнала, примерно равна . Tf управляет постоянной времени фильтра для производной, используемой для вычисления скорости. ПРИМЕЧАНИЕ. Нелинейная функциональность в ControlSystemsBase.jl в настоящее время является экспериментальной и может быть изменена без соблюдения семантического управления версиями Используйте на свой страх и риск.

# ControlSystemsBase.deadzone — Function

deadzone(val)
deadzone(lower, upper)

Создает нелинейность мертвой зоны.

       y▲
        │     /
        │    /
  lower │   /
─────|──┼──|───────► u
    /   │   upper
   /    │
  /     │

# ControlSystemsBase.linearize — Function

linearize(sys::HammersteinWienerSystem, Δy)

Линеаризовывает нелинейную систему sys вокруг рабочей точки, заданной Δy.

      ┌─────────┐
 y◄───┤         │◄────u
      │    P    │
Δy┌───┤         │◄───┐Δu
  │   └─────────┘    │
  │                  │
  │      ┌───┐       │
  │      │   │       │
  └─────►│ f ├───────┘
         │   │
         └───┘

A, B = linearize(f, x, u, args...)

Линеаризовывает динамику вокруг рабочей точки с помощью ForwardDiff. args может быть пустой или содержать, например, параметры и время (p, t), как в интерфейсе SciML. Эту функцию также можно использовать для линеаризации выходного уравнения C, D = linearize(h, x, u, args...).