Замечания по производительности
Числовая точность
Передаточные функции, как и многочлены в целом, известны плохими числовыми характеристиками. По этой причине не рекомендуется использовать передаточные функции высокого порядка. Высоким может считаться порядок уже начиная с шестого. Когда передаточная функция преобразуется в представление пространства состояний с помощью ss(G), автоматически производится уравновешивание с целью улучшить числовые характеристики модели.
Эта проблема иллюстрируется ниже: сначала мы создаем систему пространства состояний и преобразуем ее в передаточную функцию . Затем мы вносим возмущение в один элемент матрицы динамики , добавляя машинный эпсилон для Float64 (eps() = 2.22044e-16), и преобразуем эту возмущенную систему пространства состояний в передаточную функцию . Разница между этими двумя передаточными функциями огромная: векторы коэффициентов в знаменателях различаются на порядок .
sys = ssrand(1,1,100);
G1 = tf(sys);
sys.A[1,1] += eps();
G2 = tf(sys);
norm(denvec(G1)[] - denvec(G2)[])
6.270683106765845e96Если построить графики полюсов этих двух систем, они также будут существенно различаться.
scatter(poles(G1)); scatter!(poles(G2))
Если мы вместо этого вычислим полюса модели пространства состояний до и после внесения возмущения, они будут почти неотличимы друг от друга.
Уравновешивание пространства состояний
С помощью функции balance_statespace можно вычислить балансирующее преобразование , которое пытается масштабировать систему так, чтобы нормы строк и столбцов
были примерно равны. Как правило, это повышает численную производительность ряда алгоритмов, включая вычисление частотной характеристики и моделирование в непрерывном времени. При построении графиков частотных характеристик с помощью любой из встроенных функций, например bodeplot или nyquistplot, уравновешивание производится автоматически. Однако при вызове bode и nyquist напрямую за уравновешивание отвечает пользователь. Уравновешивание требует сравнительно немного затрат, однако:
-
оно изменяет представление состояний системы (но не сопоставление входов и выходов); если уравновешивание производится перед моделированием, выход будет соответствовать выходу исходной системы, но траектория состояний не будет;
-
выделяется некоторый объем памяти.
Уравновешивание также выполняется автоматически, когда передаточная функция преобразуется в систему пространства состояний с помощью ss(G). Чтобы выполнить преобразование без уравновешивания, вызовите convert(StateSpace, G, balance=false).
На интуитивном уровне (и немного упрощая) уравновешивание может быть полезно, когда величина элементов матрицы значительно отличается от величины элементов матрицы или когда матрица содержит очень большие коэффициенты. Вот пример, демонстрирующий все эти особенности.
using ControlSystemsBase, LinearAlgebra
A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]
B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]
C = [0.0 0.0 0.0 1.0 0.0 0.0]
D = [0.0]
linsys = ss(A,B,C,D)
norm(linsys.A, Inf), norm(linsys.B, Inf), norm(linsys.C, Inf)(2.824060629317756e8, 297809.51426114445, 1.0)После преобразования получается следующее.
bsys, T = balance_statespace(linsys)
norm(bsys.A, Inf), norm(bsys.B, Inf), norm(bsys.C, Inf)(6.537773175952662, 0.032156093066689026, 0.0625)Если мы построим график частотной характеристики этих двух систем с помощью bodeplot, то увидим, что они существенно различаются (правильным будет график уравновешенной системы).
Вычисление частотной характеристики
Для небольших систем (с небольшим количеством входов, выходов и состояний) частотная характеристика передаточной функции вычисляется достаточно точно и очень быстро.
G = tf(1, [1, 1])
w = exp10.(LinRange(-2, 2, 200));
@btime freqresp($G, $w);
# 4,351 мкс (2 выделения: 3,31 КиБ)Для вычисления частотной характеристики эквивалентной системы пространства состояний потребуется выделение дополнительной памяти из-за разложения матрицы Хессенберга.
sys = ss(G);
@btime freqresp($sys, $w);
# 20,820 мкс (16 выделений: 37,20 КиБ)В случае с системами большего размера вычисление пространства состояний будет существенно более точным при условии, что реализация хорошо сбалансирована.
Для оптимальной производительности можно предварительно выделить место в памяти под возвращаемый массив.
ny,nu = size(G)
R = zeros(ComplexF64, ny, nu, length(w));
@btime freqresp!($R, $G, $w);
# 4,214 мкс (1 выделение: 64 байта)Другие функции, принимающие предварительно выделенные рабочие пространства:
Пример использования bodemag!:
using ControlSystemsBase
G = tf(ssrand(2,2,5))
w = exp10.(LinRange(-2, 2, 20000))
@btime bode($G, $w);
# 55,120 мс (517 957 выделений: 24,42 МиБ)
@btime bode($G, $w, unwrap=false); # восстановление фазы выполняется медленно
# 3,624 мс (7 выделений: 2,44 МиБ)
ws = ControlSystemsBase.BodemagWorkspace(G, w)
@btime bodemag!($ws, $G, $w);
# 2,991 мс (1 выделение: 64 байта)Моделирование временной области
Временная шкала
При моделировании динамической системы в непрерывном времени применяется интегратор дифференциального уравнения. Такие интеграторы чувствительны к масштабу уравнений и могут показывать низкую производительность в случае с жесткими задачами или задачами с неоптимальной временной шкалой. Например, в электронике часто моделируются системы, в которых постоянные времени главного динамического процесса имеют величину порядка микросекунд. Для обеспечения точности часто бывает полезно моделировать такие системы в микросекундах, а не секундах. С помощью функции time_scale можно автоматически изменять временной масштаб системы.
Передаточные функции
Передаточные функции автоматически преобразуются в форму пространства состояний перед моделированием временной области. Если вы хотите в точности контролировать внутреннее представление, можно с самого начала моделировать систему как систему пространства состояний.
Моделирование в дискретном времени
Линейные системы с входами с фиксацией нулевого порядка можно точно моделировать в дискретном времени. Вы можете указать дискретизацию ZoH в вызове lsim с помощью аргумента method=:zoh или выполнить дискретизацию вручную с помощью c2d. Моделирование в дискретном времени часто выполняется гораздо быстрее, чем интегрирование в непрерывном времени.
Для систем с дискретным временем функция lsim! принимает предварительно выделенные объекты рабочих пространств, которые позволяют избежать выделения памяти при повторном моделировании.
Численное уравновешивание
Если вас интересуют только смоделированные выходы, но не траектории состояний, вы можете применить уравновешивание к модели пространства состояний с помощью balance_statespace перед моделированием. См. раздел State-space balancing выше. Если траектории состояний также представляют интерес, перед моделированием можно также произвести уравновешивание, а после него применить к траекториям состояний обратное преобразование.