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

Автоматическое дифференцирование

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

Существует два общих метода автоматического дифференцирования: прямой и обратный режимы. Прямой режим алгоритмически более предпочтителен для функций с малым количеством входов, но большим количеством выходов, а обратный режим более эффективен для функций с большим количеством параметров, но малым количеством выходов (как в глубоком обучении). В Julia автоматическое дифференцирование в прямом режиме обеспечивается пакетом ForwardDiff.jl, а автоматическое дифференцирование в обратном режиме поддерживается несколькими различными пакетами, такими как Zygote.jl и ReverseDiff.jl. Автоматическое дифференцирование в прямом режиме обычно связано с меньшими затратами, чем автоматическое дифференцирование в обратном режиме, поэтому для функций с небольшим числом параметров, скажем, менее 10 или 100, прямой режим обычно является наиболее эффективным. ForwardDiff.jl также поддерживает дифференцирование большей части языка Julia, поэтому вероятность успеха здесь выше, чем при использовании других пакетов, и обычно рекомендуется сначала попробовать ForwardDiff.jl.

Линеаризация нелинейной динамики

Нелинейная динамика вида

легко приводится к линейному виду в точке с помощью ForwardDiff.jl:

using ControlSystemsBase, ForwardDiff

"An example of nonlinear dynamics"
function f(x, u)
    x1, x2 = x
    u1, u2 = u
    [x2; u1*x1 + u2*x2]
end

x0 = [1.0, 0.0] # Рабочая точка, вокруг которой выполняется линеаризация
u0 = [0.0, 1.0]

A = ForwardDiff.jacobian(x -> f(x, u0), x0)
B = ForwardDiff.jacobian(u -> f(x0, u), u0)

"An example of a nonlinear output (measurement) function"
function g(x, u)
    y = [x[1] + 0.1x[1]*u[2]; x[2]]
end

C = ForwardDiff.jacobian(x -> g(x, u0), x0)
D = ForwardDiff.jacobian(u -> g(x0, u), u0)

linear_sys = ss(A, B, C, D)
StateSpace{Continuous, Float64}
A =
 0.0  1.0
 0.0  1.0
B =
 0.0  0.0
 1.0  0.0
C =
 1.1  0.0
 0.0  1.0
D =
 0.0  0.1
 0.0  0.0

Continuous-time state-space model

В приведенном примере f линеаризуется в точке для получения линейных матриц пространства состояний и , а также линеаризуется g для получения линейных выходных матриц и . Вместо того чтобы вручную вызывать ForwardDiff.jl для линеаризации динамики, пользователь может вызвать функцию ControlSystemsBase.linearize, которая включает в себя необходимые вызовы ForwardDiff.jl.

Оптимизационная настройка — ПИД-регулятор

В этом примере демонстрируется простое использование автоматического дифференцирования с помощью ForwardDiff.jl для автоматической настройки ПИД-регулятора на основе оптимизации.

Система, которой мы будем управлять, представляет собой систему с удвоенной массой, в которой две массы (или инерции) соединены гибкой передачей.

Начнем с определения модели системы и начального предположения для параметров ПИД-регулятора.

using ControlSystemsBase, ForwardDiff, Plots

P = DemoSystems.double_mass_model()

bodeplot(P, title="Bode plot of Double-mass system \$P(s)\$")

Ω = exp10.(-2:0.03:3)
kp,ki,kd,Tf =  1, 0.1, 0.1, 0.01 # параметры регулятора

C  = pid(kp, ki, kd; Tf, form=:parallel, state_space=true) # Построение ПИД-регулятора с фильтром
G  = feedback(P*C) # Замкнутая система
S  = 1/(1 + P*C)   # Функция чувствительности
Gd = c2d(G, 0.1)   # Дискретизация системы
res = step(Gd,15)  # Ступенчатый отклик

mag = bodev(S, Ω)[1]
plot(res, title="Time response", layout = (1,3), legend=:bottomright)
plot!(Ω, mag, title="Sensitivity function", xscale=:log10, yscale=:log10, subplot=2, legend=:bottomright, ylims=(3e-2, Inf))
Ms, _ = hinfnorm(S)
hline!([Ms], l=(:black, :dash), subplot=2, lab="\$M_S = \$ $(round(Ms, digits=3))", sp=2)
nyquistplot!(P*C, Ω, sp=3, ylims=(-2.1,1.1), xlims=(-2.1,1.2), size=(1200,400))

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

Начнем с определения вспомогательной функции plot_optimized, которая будет оценивать производительность настроенного регулятора. Далее определим функцию systems, которая строит группу из четырех передаточных функций (extended_gangoffour) и выполняет моделирование во временной области передаточных функций и , т. е. передаточных функций от эталона к ошибке управления и передаточной функции от входного возмущения по нагрузке к ошибке управления . Оптимизируя эти ступенчатые отклики с учетом параметров ПИД-регулятора, мы получим регулятор с хорошей производительностью. Для повышения надежности замкнутого контура, а также ограничения усиления шумов измерений в управляющем сигнале мы снижаем пик функции чувствительности и (приближенно) частотно-взвешенную норму передаточной функции .

Ограничительная функция constraints делает так, что пик функции чувствительности будет находиться ниже Msc. Наконец, мы применяем Optimization.jl для оптимизации функции стоимости и указываем ей использовать ForwardDiff.jl для вычисления градиента функции стоимости. В этом примере используется оптимизатор Ipopt.

using Optimization, Statistics, LinearAlgebra
using Ipopt, OptimizationMOI; MOI = OptimizationMOI.MOI

function plot_optimized(P, params, res, systems)
    fig = plot(layout=(1,3), size=(1200,400), bottommargin=2Plots.mm)
    for (i,params) = enumerate((params, res))
        ls = (i == 1 ? :dash : :solid)
        lab = (i==1 ? "Initial" : "Optimized")
        C, G = systems(params, P)
		r1, r2 = sim(G)
        mag = reshape(bode(G, Ω)[1], 4, :)'[:, [1, 2, 4]]
        plot!([r1, r2]; title="Time response", subplot=1,
            lab = lab .* [" \$r → e\$" " \$d → e\$"], legend=:bottomright, ls,
            fillalpha=0.05, linealpha=0.8, seriestype=:path, c=[1 3])
        plot!(Ω, mag; title="Sensitivity functions \$S(s), CS(s), T(s)\$",
            xscale=:log10, yscale=:log10, subplot=2, lab, ls,
            legend=:bottomright, fillalpha=0.05, linealpha=0.8, c=[1 2 3], linewidth=i)
        nyquistplot!(P*C, Ω; Ms_circles=Msc, sp=3, ylims=(-2.1,1.1), xlims=(-2.1,1.2), lab, seriescolor=i, ls)
    end
    hline!([Msc], l=:dashdot, c=1, subplot=2, lab="Constraint", ylims=(9e-2, Inf))
    fig
end

"A helper function that creates a PID controller and closed-loop transfer functions"
function systemspid(params, P)
    kp,ki,kd,Tf = params # Оптимизируем параметры в
    C    = pid(kp, ki, kd; form=:parallel, Tf, state_space=true)
    G    = extended_gangoffour(P, C) # [S PS; CS T]
    C, G
end

"A helper function that simulates the closed-loop system"
function sim(G)
    Gd = c2d(G, 0.1, :tustin)   # Дискретизация системы
    res1 = step(Gd[1, 1], 0:0.1:15) # Моделирование S
    res2 = step(Gd[1, 2], 0:0.1:15) # Моделирование PS
    res1, res2
end

"The cost function to optimize"
function cost(params::AbstractVector{T}, (P, systems)) where T
    CSweight = 0.001 # Снижение усиления шума
    C, G = systems(params, P)
    res1, res2 = sim(G)
    R, _ = bodev(G[2, 1], Ω; unwrap=false)
    CS = sum(R .*= Ω)               # частотно-взвешенная чувствительность к шуму
    perf = mean(abs2, res1.y .*= res1.t') + mean(abs2, res2.y .*= res2.t')
    return perf + CSweight * CS # Объединим все цели
end

"The sensitivity constraint to enforce robustness"
function constraints(res, params::AbstractVector{T}, (P, systems)) where T
    C, G = systems(params, P)
    S, _ = bodev(G[1, 1], Ω; unwrap=false)
    res .= maximum(S) # максимальная чувствительность
    nothing
end

Msc = 1.3        # Ограничение для Ms

params  = [kp, ki, kd, 0.01] # Начальное предположение для параметров

solver = Ipopt.Optimizer()
MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), 0)
MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), 200)
MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_tol"), 1e-1)
MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_constr_viol_tol"), 1e-2)
MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_iter"), 5)
MOI.set(solver, MOI.RawOptimizerAttribute("hessian_approximation"), "limited-memory")

fopt = OptimizationFunction(cost, Optimization.AutoForwardDiff(); cons=constraints)

prob = OptimizationProblem(fopt, params, (P, systemspid);
    lb    = fill(-10.0, length(params)),
    ub    = fill(10.0, length(params)),
    ucons = fill(Msc, 1),
    lcons = fill(-Inf, 1),
)

res = solve(prob, solver)
plot_optimized(P, params, res.u, systemspid)

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

Оптимизационная настройка — ЛКГ-регулятор

Можно попытаться провести аналогичную автоматическую настройку ЛКГ-регулятора. На этот раз мы выбираем для оптимизации весовые матрицы задачи ЛКР и ковариационную матрицу состояния шума. Синтез ЛКР-регулятора предполагает решение уравнения Рикатти, которое, в свою очередь, предполагает выполнение разложения Шура. Эти этапы трудно продифференцировать обычным способом, но можно воспользоваться неявным дифференцированием, используя теорему о неявной функции. Для этого загружаем пакет ImplicitDifferentiation и определяем условия, которые выполняются при решении уравнения Рикатти:

При загрузке ImplicitDifferentiation автоматически загружаются дифференцируемые версии lqr и kalman, использующие неявную функцию.

using ImplicitDifferentiation, ComponentArrays # Оба эти пакета необходимы для загрузки правил неявного дифференцирования.

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

function triangular(x)
    m = length(x)
    n = round(Int, sqrt(2m-1))
    T = zeros(eltype(x), n, n)
    k = 1
    for i = 1:n, j = i:n
        T[i,j] = x[k]
        k += 1
    end
    T
end
invtriangular(T) = [T[i,j] for i = 1:size(T,1) for j = i:size(T,1)]

function systemslqr(params::AbstractVector{T}, P) where T
    n2 = length(params) ÷ 2
    Qchol = triangular(params[1:n2])
    Rchol = triangular(params[n2+1:2n2])
    Q = Qchol'Qchol
    R = Rchol'Rchol
    L = lqr(P, Q, one(T)*I(1)) # Важно, чтобы последняя матрица имела правильный тип.
    K = kalman(P, R, one(T)*I(1))
    C = observer_controller(P, L, K)
    G = extended_gangoffour(P, C) # [S PS; CS T]
    C, G
end

Q0 = diagm([1.0, 1, 1, 1]) # Начальное предположение штрафа по состоянию ЛКР
R0 = diagm([1.0, 1, 1, 1]) # Начальное предположение ковариации состояния Калмана
params2 = [invtriangular(cholesky(Q0).U); invtriangular(cholesky(R0).U)]

prob2 = OptimizationProblem(fopt, params2, (P, systemslqr);
    lb    = fill(-10.0, length(params2)),
    ub    = fill(10.0, length(params2)),
    ucons = fill(Msc, 1),
    lcons = fill(-Inf, 1),
)

res2 = solve(prob2, solver)
plot_optimized(P, params2, res2.u, systemslqr)

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

Анализ надежности

Для проверки надежности спроектированного ЛКГ-регулятора в отношении параметрической неопределенности в объекте загрузим пакет MonteCarloMeasurements и воссоздадим модель объекта с 20-процентной неопределенностью в коэффициенте упругости.

using MonteCarloMeasurements
Pu = DemoSystems.double_mass_model(k = Particles(32, Uniform(80, 120))) # Создание модели с неопределенностью в жесткости упругого элемента k ~ U(80, 120)
unsafe_comparisons(true) # Чтобы график Боде работал

C,_ = systemslqr(res2.u, P)             # Получение регулятора с предположением, что P без неопределенности
Gu = extended_gangoffour(Pu, C)     # Из Большой четверки с неопределенностью
w = exp10.(LinRange(-1.5, 2, 500))
bodeplot(Gu, w, plotphase=false, ri=false, N=32, ylims=(1e-1, 30), layout=1, sp=1, c=[1 2 4 3], lab=["S" "CS" "PS" "T"])
hline!([Msc], l=:dashdot, c=1, lab="Constraint", ylims=(9e-2, Inf))