news.shamcode.ru | Модель ПИД регулятора на Python

Модель ПИД регулятора на Python

понедельник, 15 мая 2017 г. в 03:13:50

Python

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

Моделированию работы ПИД регулятора посвящено большое количество публикаций в сети. Лидирует проектирование моделей ПИД регулятора с применением 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].

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

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

Единичное входное воздействие 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 сведём в таблицу.

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

Код программы

#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def mnkLIN(x,y,q): # аппроксимация характеристик статических объектов МНК                            
              a=round((len(x)*sum([x[i]*y[i] for i in range(0,len(x))])-sum(x)*sum(y))/(len(x)*sum([x[i]**2 for i in range(0,len(x))])-sum(x)**2),3)
              b=round((sum(y)-a*sum(x))/len(x) ,3)
              y1=[round(a*w+b ,3) for w in x]         
              s=[round((y1[i]-y[i])**2,3) for i in range(0,len(x))]                          
              sko=round((sum(s)/(len(x)-1))**0.5,3)
              if q==1:
                       plt.subplot(221)
                       plt.plot(x, y, color='r', linewidth=2, marker='o', label='Коэффициент a= %s'%str(a))
                       plt.plot(x, y1, color='g', linestyle=' ', marker='o', label='СКО =%s'%str(sko))
                       plt.legend(loc='best')
                       plt.grid(True)
              else:
                      plt.subplot(222)
                      plt.plot(x, y, color='r', linewidth=2, marker='o', label='Коэффициент a= %s'%str(a))
                      plt.plot(x, y1, color='g', linestyle=' ', marker='o', label='СКО =%s'%str(sko))
                      plt.legend(loc='best')
                      plt.grid(True)                      
              return a             
T1=7.0;T2=5.0;s=0.4;k1=5.5;k2=5.5# параметры  объекта регулирования.
x=[0,1,2,3,4,5]
y=[0,-11,-22,-33,-44,-55]
a2=mnkLIN(x,y,1)
x=[0,1,2,3,4,5]
y=[0,17,34,51,68,85]
a3=mnkLIN(x,y,2)
#m1=0.5;m2=2.0;m3=0.3# параметры  ПИД регулятора.
m1=0.0;m2=0.0;m3=0.0
a=a2+a3
A=round((s*T1*T2+T1**2)/(T1*T2**2),3)
B=round((T2+s*T1+k1*k2*m2*a)/(T1*T2**2),3)
C=round((1+k1*k2*m1*a)/(T1*T2**2),3)
D=round(k1*k2*m3*a/(T1*T2**2),3)
E=round(k1*k2*a/(T1*T2**2),3)
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] # вектор значений производной
plt.subplot(223)
plt.plot(t,y1,linewidth=1, label=' p- %s,d- %s,i- %s,'%(str(m1),str(m2),str(m3)))
plt.ylabel("z")
plt.xlabel("t")
plt.legend(loc='best')
plt.grid(True)
plt.show()

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

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

Изменим табличные данные для третьего звена, увеличим а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 может эффективно использоваться для численного моделирования контуров автоматического регулирования содержащих динамические и статические звенья.