В поисках простой модели ПИД регулятора с объектом

Моделированию работы ПИД регулятора посвящено большое количество публикаций в сети. Лидирует проектирование моделей ПИД регулятора с применением Matlab Simulink [1,2] (134 миллиона ссылок в yandex). Сам процесс создания модели какой-то однообразный. В модель переносят всё новые и новые блоки. Одно движение ручного манипулятора и нате вам ПИД контролер, ещё одно и вот передаточная функция объекта. Соединяешь блоки, настраиваешь параметры, готовишь вычислитель. Да, возможностей много, но как-то слишком всё искусственно. И уже становится совсем непонятным к чему тут дифференциальные уравнения, методы их решения и то операционное исчисление, которым долго морочили голову. Ищу реализацию ПИД в Mathcad, тут ссылок в том же yandex, поменьше, всего то 81 миллион, а математики и формул побольше. Рассматриваю пример ПИД, поставляемый вместе с пакетом Mathcad 14.


В качестве объекта колебательное звено. Много умных объяснений, но в итоге два оператора laplace и invlaplace. Общая передаточная функция имеет в числителе вторую степень оператора, а в знаменателе четвёртую. Чтобы операторы laplace и invlaplace сработали, когда подключены все три составляющих ПИД, находят ещё и корни знаменателя передаточной функции, эти корни комплексно сопряжённые.



Теперь ищу реализацию ПИД на Python. Тихо радовался 97 миллионам результатов, но не долго. О Python 2.7 только применительно к прошивке Arduino на примере ESP32. Но и это переполняет сердце гордостью за Python.

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

Структурная схема объекта управления

Постановку задачи взял из [3], в ней меня привлекло разнообразие структуры объекта регулирования. По ходу исправил несколько ошибок. Так, в конечном дифференциальном уравнении был потерян квадрат при T2, данные таблицы для a2 не соответствовали результатам аппроксимации. Кроме того, ввёл установку управления в виде единичного бесконечного импульса, чтобы оживить картинку переходного процесса. Уж как-то не естественно для ПИД регулятора она выглядит в [3].


Будем исследовать систему регулирования со следующей структурной схемой.

1. Колебательное звено [4] описывается дифференциальным уравнением второго порядка.
2. Линейное статическое звено, определяющие зависимость выхода от входа которого определяется по табличным данным. Данные могут быть получены экспериментальным путём. Из-за погрешностей снятия данных, следует применять линейную аппроксимацию .
3. Второе линейное статическое звено — зависимость выхода от входа определяется по табличным данным, которые могут быть получены экспериментальным путём.. Из-за погрешностей снятия данных, следует применять линейную аппроксимацию .
4. Линейное динамическое звено, на вход которого поступает сума сигналов от блоков 2,3
   Звено имеет передаточную функцию, которая описывается линейным дифференциальным уравнением первого порядка.
   
5. ПИД регулятор с известным законом регулирования, при котором для настройки на объект регулирования используются коэффициент пропорциональности при линейной функции, коэффициент дифференцирования при первой производной и коэффициент интегрирования при неопределённом интеграле.
   

Единичное входное воздействие 1(t) преобразуется в операторную форму как 1 / p. Переведём в операторную форму уравнения 1,4 в результате получим.



Продифференцируем обе части (5), для этого просто умножим на оператор p, получим.



Переведём в операторную форму уравнения (2), (3) и выразим y(y1).



Подставим (7) в (6) получим:



Перейдём (8) во временную область, для этого операторную форму заменим дифференциальным уравнением.

Решение задачи определения настроек ПИД регулятора на Python

Зададим следующие параметры объекта регулирования:
T1=7.0; T2=5.0; s=0.4; k1=5.5; k2=5.5. Данные для аппроксимации характеристик звеньев 2,3 сведём в таблицу.



Отставим объект без регулятора, для этого примем .



Результат роботы программы.



Без ПИД регулятора объект пошёл в «разнос». После предварительного подбора параметров ПИД получим.



Изменим табличные данные для третьего звена, увеличим а3 в два раза.



Регулятор стабилен.

Выбор функции для решения уравнения (9) на Python и сравнение результатов с решением на Matcad

Я использовал функцию odeint из библиотеки scipy. Integrate.

def f(y, t):#  решение дифференциального уравнения системы регулирования.
                y1,y2,y3,y4 = y 
                return [y2,y3,y4,-A*y4-B*y3-C*y2-D*y1+E]
t = np. linspace(0,100,10000)
y0 = [0,1,0,0]#начальные условия
w = odeint(f, y0, t)
y1=w[:,0] # вектор значений решения
y2=w[:,1] # вектор значений производной

Если уравнение (9) перестроить так, чтобы в левой части была производная высшего порядка, а в правой всё остальное, то получим.



Так и строится функция def f (y, t). Функция odeint имеет три обязательных аргумента и выглядит вот так odeint (func, y0, t[,args=(), ...]).

Для сравнения выведем используемые в приведенной функции def f (y, t) значения переменных A, B, C, D, E, получим 0.36 7.996 1.994 1.193 3.976. Составим на Mathcad 14 программу для решения дифференциального уравнения (9) и построим переходную характеристику.


Графики и числовые значения полностью идентичны. Таким образом, функции Odesolve Mathcad и odeint Python дают одинаковые численные решения приведённой в статье задачи.

Вывод

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

Ссылки

1. Автоматизированная настройка ПИД-регулятора для объекта управления следящей системы с использованием программного пакета MATLAB Simulink Филиппов А. В.1, Косолапов М. А.2, Маслов И. А.3, Тарасова Г. И.
2. Методы настройки ПИД-регулятора в matlab simulink.
3. Подбор ПИД регулятора.
4. Модель колебательного звена в режиме резонансных колебаний на Python.