Терморегуляция
Стандартная модель системы с терморегуляцией имеет следующий вид:
где — температура, — управляющий сигнал, а — постоянное смещение, например связанное с температурой среды, окружающей регулируемую систему. Временная постоянная отражает связь между потенциальной энергией и теплостойкостью и определяет то, насколько быстро меняется температура. Эту систему можно записать в форме передаточной функции (опустив ) следующим образом:
Это простая передаточная функция первого порядка, которую можно оценить, например, с помощью функции arx или plr. С целью иллюстрации мы создадим такую систему и смоделируем данные на ее основе.
using ControlSystemsBase, ControlSystemIdentification, Plots
w = 2pi .* exp10.(LinRange(-3, log10(0.5), 500))
G0 = tf(1, [10, 1]) # Истинная система, 10ẋ = -x + u
G = c2d(G0, 1) # Производим дискретизацию с интервалом выборки 1 с
println("True system")
display(G0)
u = sign.(sin.((0:0.01:20) .^ 2))' # Выполняем выборку входного управляющего сигнала для идентификации
y, t, x = lsim(ss(G), u) # Моделируем истинную систему, чтобы получить тестовые данные
yn = y .+ 0.2 .* randn.() # Добавляем шум измерения
data = iddata(yn, u, t[2] - t[1]) # Создаем объект данных
plot(data)
Мы видим, что данные, которые будут использоваться для идентификации, представляют собой вход с линейной частотной модуляцией. Сигналы с линейной частотной модуляцией отлично подходят для идентификации, так как они имеют для этого четко определенный и легко контролируемый интервал частот. Сначала изучим график когерентности, чтобы убедиться в том, что данные пригодны для идентификации линейной системы.
coherenceplot(data, hz=true)
Когерентность высока для всех частот, охватываемых сигналом с линейной частотной модуляцией, после чего существенно снижается. Из этого следует, что полагаться на точность идентифицированной модели можно только до максимальной частоты, имеющейся на входе с линейной частотной модуляцией.
Далее зададим параметры для оценки: в числителе и знаменателе имеется по одному параметру, поэтому мы задаем и оцениваем обе модели.
na, nb = 1, 1 # Количество параметров в числителе и знаменателе
Gh = arx(data, na, nb, estimator = wtls_estimator(data.y, na, nb)) # Оцениваем модель ARX
Gh2, noise_model = plr(data, na, nb, 1) # Пробуем другой метод идентификации
Gh, Gh2(TransferFunction{ControlSystemsBase.Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
0.09393634941260658
-------------------------
1.0z - 0.9057691033924613
Sample Time: 1 (seconds)
Discrete-time transfer function model, TransferFunction{ControlSystemsBase.Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
0.09650004560797866
-------------------------
1.0z - 0.8999349368352745
Sample Time: 1 (seconds)
Discrete-time transfer function model)Оценка моделей ARX методом наименьших квадратов на основе данных с высоким уровнем шума измерения, как известно, обеспечивает низкий уровень соответствия моделей на низких частотах. Поэтому мы использовали вызов wtls_estimator(data.y, na, nb), выполняющий оценку методом наименьших полных квадратов.
Построить график результатов можно несколькими разными способами:
julia> # Строим график результатов
println("Estimated system in continuous time")
Estimated system in continuous time
julia> display(d2c(Gh)) # Преобразуем дискретную систему в непрерывнуюbp = bodeplot(G, w, lab = "G (true)", hz = true, l = 5)
bodeplot!(Gh, w, lab = "arx", hz = true)
bodeplot!(Gh2, w, lab = "plr", hz = true, ticks = :default)
sp = plot(step(G, 150), lab="G (true)")
plot!(step(Gh, 150), lab = "arx")
plot!(step(Gh2, 150), lab = "plr", ticks = :default)
hline!([1], primary = false, l = (:black, :dash))
lp = plot(lsim(ss(G), u), lab="G (true)")
plot!(lsim(ss(Gh), u), lab = "arx")
plot!(lsim(ss(Gh2), u), lab = "plr", ticks = :default)
plot!(data.t, yn[:], lab = "Estimation data", alpha=0.3)
plot(bp, sp, lp, layout = @layout([[a b]; c]))