Использование обратного преобразования Лапласа для анализа динамических звеньев систем управления
До настоящего времени в арсенале средств высокоуровневого языка программирования Python отсутствовали модули для численного преобразования передаточных функций элементов САУ из частотной области во временную.
Поскольку функции обратного преобразования Лапласа широко используются при анализе динамических систем контроля измерения и управления, использование Python для указанных целей было весьма затруднительно, поскольку приходилось использовать менее точное обратное Фурье преобразование [1].
Указанную проблему решает модуль mpmath библиотеки Python свободного распространения (под лицензией BSD), предназначенный для решения задач вещественной и комплексной арифметики с плавающей точкой и заданной точностью.
Работу над модулем ещё в 2007 году начал Fredrik Johansson [2], и, благодаря помощи многих участников проекта, в настоящее время mpmath приобрёл возможности серьёзного математического пакета. Однако нас будет интересовать только заявленная в статье тема, реализуемая при помощи одношагового алгоритма invertlaplace. В качестве примера рассмотрим сравнение результатов обратного преобразования Лапласа передаточной функции w(p)=1/(1+p)** 2 при помощи invertlaplace с известной переходной функцией h(t)=e**-t от указанной передаточной функции:
Значение тестовой функции: 0.000999000499833375 Полученное значение функции: 0.000999000499833375 Значение h(t): 0.000999000499833375. Абсолютная погрешность :8.57923043561212e-20. Значение тестовой функции: 0.00990049833749168 Полученное значение функции: 0.00990049833749168 Значение h(t): 0.00990049833749168. Абсолютная погрешность :3.27007646698047e-19. Значение тестовой функции: 0.090483741803596 Полученное значение функции: 0.090483741803596 Значение h(t): 0.090483741803596. Абсолютная погрешность: -1.75215800052168e-18. Значение тестовой функции: 0.367879441171442 Полученное значение функции: 0.367879441171442 Значение h(t): 0.367879441171442. Абсолютная погрешность :1.2428864009344e-17. Значение тестовой функции: 0.000453999297624849 Полученное значение функции: 0.000453999297624849 Значение h(t): 0.000453999297624849. Абсолютная погрешность :4.04513489306658e-20. Время, затраченное на обратное преобразование Лапласа: 0.18808794021606445
Тестовый пример ограничен по объёму, но и на нём видно, что одношаговый алгоритм invertlaplace имеет высокую точность и не критичное время выполнения для ограниченного числа значений времени.
В рассмотренном примере был использован метод talbot, об особенностях других методов можно прочесть в документации [3].
Однако следует учитывать, что все численные методы обратного преобразования Лапласа требуют, чтобы их абсцисса сдвигалась ближе к началу координат для больших времен. Если абсцисса сдвинется влево от самой правой координаты в области Лапласа, ответ будет совершенно неверным.
Поэтому необходимо провести исследование применимости численного обратного преобразования Лапласа для определённых передаточных функций и оценить погрешность в сравнении с другим методом, в данной публикации это метод обратного преобразования Фурье.
1. Построение переходной характеристики объекта управления по его передаточной функции с использованием invertlaplaceДопустим, у нас есть водно-водяной теплообменник с передаточной функцией по каналу, температура нагреваемой воды – расход греющей воды.Преобразование Лапласа выходного сигнала — передаточная функция объекта управления умноженная на 1/p (единичного возмущения по расходу греющей воды) с учётом запаздывания τ имеет следующий вид:
где: Ti – постоянные времени звеньев; K – статический коэффициент передачи; p ¬ оператор Лапласа.
Получим переходную характеристику объекта с запаздыванием и самовыравниванием:
2. Построение переходной характеристики ПИД- регулятора по его передаточной функции с использованием invertlaplaceИзображение по Лапласу передаточной функции ПИД регулятора имеет вид:
где: Td, Ti – постоянные времени, дифференцирующего и интегрирующего звеньев; Kp, Kd – статические коэффициенты передачи пропорционального и дифференцирующего звеньев; p – оператор Лапласа.
По приведенному листингу можно получить не только переходную характеристику ПИД- регулятора, но и ПИ-регулятора, приняв Kd=0:
3. Оценка точности численного метода обратного преобразования Лапласа invertlaplace применительно к типовым объектам САУС учётом предостережений относительно применимости invertlaplace, приведенных в начале статьи, в дальнейшем изложении была доказана применимость метода к объектам с запаздыванием и самовыравниванием, а так же к регуляторам.
Остаётся до конца не выясненным вопрос точности численного решения. Для внесения ясности воспользуемся передаточной функцией (2) и следующим точным её преобразованием в переходную функцию, приведённую в [1]:
Получим следующий график:
Из графика видно, что погрешность численного метода для приведенного класса передаточных функций пренебрежимо мала (3*10^-15). Кроме того, погрешностью численного метода можно управлять, устанавливая значения mp.dps и шаг tt[i+1]-tt[i] (см. листинг).
4. Сравнительная оценка точности численного метода обратного преобразования Лапласа invertlaplace и численного метода обратного преобразования Фурье применительно к типовым объектам САУПостроение переходной характеристики можно произвести на основе формулы обратного преобразования Фурье [1]:
где X(j∙ω)— Фурье-изображение оригинала x(t)
где Re(W(j∙ω))— вещественная частотная характеристика объекта регулирования.
В качестве верхнего предела интегрирования θв в расчете берется значение частоты ωс, при котором модуль Re(W(j∙ω)) уменьшается до некоторого малого значения (например, 0,05*K) и не превосходит это значение при дальнейшем возрастании ω.
Вначале определим верхний предел интегрирования в соотношении (5). Для этого воспользуемся передаточной функцией (1), предварительно избавившись от единичного воздействия, умножив правую часть на оператор p.
Из графика видно, что верхний предел интегрирования в соотношении (5) может быть принят равным 0.6, поскольку при дальнейшем увеличении частоты вещественная часть передаточной функции объекта управления сохраняет нулевое значение.
Используя (5), получим переходную характеристику объекта управления методом обратного преобразования Фурье:
Для сравнительной оценки точности обратного преобразования Лапласа и Фурье, воспользуемся точным преобразованием передаточной функции (1) объекта управления приведенном в [1]:
Из точного преобразования (6) вычтем по очереди обратное преобразование (1) Лапласа и обратное преобразование Фурье. Для построения графика, характеризующего сравнение точности обеих методов, составим следующий листинг программы:
Из графика следует, что численное обратное преобразование Лапласа является более стабильным, чем Фурье и, кроме этого, оно совпадает с точным решением.
Выводы1. Проведен сравнительный анализ численных методов обратного преобразования Лапласа и Фурье, показана большая точность и стабильность обратного преобразования Лапласа. 2. Показаны возможности применения библиотеки mpmath Python для обратного преобразования Лапласа основных передаточных функций объектов и элементов САУ.