Как стать автором
Обновить

Оптимизация настроек ПИД регулятора по интегральному критерию качества регулирования

Время на прочтение 12 мин
Количество просмотров 12K


Постановка задачи оценки качества регулирования


Интегральные оценки качества регулирования характеризует суммарное отклонение реального переходного процесса в системе от идеализированного переходного процесса.

В качестве идеализированного процесса обычно принимается ступенчатый (скачкообразный) переходный процесс или экспоненциальный процесс с заданными параметрами экспоненты.

До настоящего времени остаётся не выясненным, какой из интегральных критериев линейный или квадратичный выбрать для оптимизации настроек регуляторов.

В данной публикации определение оптимальных настроек регулятора осуществляется на основе критерия оптимальности в виде суммы линейного и квадратичного интегральных показателей качества регулирования.

Чтобы осуществить такой выбор, вначале необходимо провести серию предварительных расчетов для различных значений отношения постоянных времени дифференцирования и интегрирования α=Td/Ti и выбрать несколько значений α, которые обеспечивает наибольшее значение отношения Kr/Ki.

После этого, имея несколько вариантов настроек регулятора, построить для них переходные процессы в замкнутой системе и по предложенному критерию отобрать оптимальные настройки.

Программная реализация расчёта оптимальных настроек регулятора


Расчёт можно выполнить как в математическом пакете, так и на любом языке программирования. Однако, по причине наличия мощных библиотек для численного интегрирования и простой реализации поиска локальных экстремумов в списках, я выбрал высокоуровневый язык программирования Python.

Для упрощения поисковых процедур в начале программы целесообразно установить исходные данные и привести готовые списки частот отдельно для каждой решаемой задачи:

Подготовка исходных данных
# -*- coding: utf8 -*-    
import numpy as np # модуль для работы со списками
from scipy.integrate import quad # модуль для численного  интегрирования
import matplotlib.pyplot as plt # модуль для построения графиков
T1=14; T2=18; T3=28; K=0.9; tau=6.4# Постоянные времени, коэффициент усиления, запаздывание
m=0.366; m1=0.221# Запас устойчивости
n= np.arange(0.001,0.15,0.0002)# Список частот для плоскости Kr-Ki
n1=np.arange(0.00001,0.12,0.0001) # Список частот для графика Ki=f(w)
n2=np.arange(0.0002,0.4,0.0001) #Массив частот для построения АЧХ


Вводим передаточную функцию для водяного теплообменника с учётом корневого показателя колебательности m и передаточную функцию ПИД регулятора с учётом Kr, Ti, Td, остальные функции вспомогательные:

Ввод основных и вспомогательных функций с учётом комплексной арифметики
def WO(m,w):#Передаточная функция объекта
         j=(-1)**0.5
         return K*np.exp(-tau*(-m+j)*w)/((T1*(-m+j)*w+1)*(T2*(-m+j)*w+1)*(T3*(-m+j)*w+1))
def WR(w,Kr,Ti,Td):#Передаточная функция регулятора
         j=(-1)**0.5
         return Kr*(1+1/(j*w*Ti)+j*w*Td)
def ReW(m,w):#Действительная часть передаточной функции
          j=(-1)**0.5
          return WO(m,w).real
def ImW(m,w):#Мнимая часть передаточной функции
          j=(-1)**0.5
          return WO(m,w).imag
def A0(m,w):#Вспомогательная функция
         return -(ImW(m,w)*m/(w+w*m**2)+ReW(m,w)/(w+w*m**2))
def Ti(alfa,m,w):#Коэффициент регулятора
         return  (-ImW(m,w)-(ImW(m,w)**2-4*((ReW(m,w)*alfa*w-ImW(m,w)*alfa*w*m)*A0(m,w)))**0.5)/(2*(ReW(m,w)*alfa*w-ImW(m,w)*alfa*w*m))
def Ki(alfa,m,w):#Коэффициент регулятора
         return 1/(w*Ti(alfa,m,w)**2*alfa*(m*ReW(m,w)+ImW(m,w))-Ti(alfa,m,w)*ReW(m,w)+(m*ReW(m,w)-ImW(m,w))/(w+w*m**2))
def Kr(alfa,m,w):#Коэффициент регулятора
         if Ki(alfa,m,w)*Ti(alfa,m,w)<0:
                  z=0
         else:
                  z=Ki(alfa,m,w)*Ti(alfa,m,w)
         return z         
def Kd(alfa,m,w):#Коэффициент регулятора
         return alfa*Kr(alfa,m,w)*Ti(alfa,m,w)


Построим три плоскости настройки ПИД регулятора для отношений alfa=Td/Ti=0.2,alfa=Td/Ti=0.7 alfa=Td/Ti=1.2.

Поиск оптимальных значений alfa
alfa=0.2
Ki_1=[Ki(alfa,m1,w) for w in n]
Kr_1=[Kr(alfa,m1,w) for w in n]
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
Ki_3=[Ki(alfa,0,w) for w in n]
Kr_3=[Kr(alfa,0,w) for w in n]    
plt.figure()
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для alfa=%s"%alfa)
plt.axis([0.0, round(max(Kr_3),4), 0.0, round(max(Ki_3),4)])
plt.plot(Kr_1, Ki_1, label='Линия запаса устойчивости m=%s'%m1)
plt.plot(Kr_2, Ki_2, label='Линия запаса устойчивости m=%s'%m)
plt.plot(Kr_3, Ki_3, label='Линия  границы устойчивости m=0')
plt.xlabel("Коэффициенты - Ki")
plt.ylabel("Коэффициенты - Kr")
plt.legend(loc='best')
plt.grid(True)
alfa=0.7
Ki_1=[Ki(alfa,0.221,w) for w in n]
Kr_1=[Kr(alfa,0.221,w) for w in n]
Ki_2=[Ki(alfa,0.366,w) for w in n]
Kr_2=[Kr(alfa,0.366,w) for w in n]
Ki_3=[Ki(alfa,0,w) for w in n]
Kr_3=[Kr(alfa,0,w) for w in n]
plt.figure()
plt.axis([0.0, round(max(Kr_3),3), 0.0, round(max(Ki_3),3)])
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для alfa=%s"%alfa)
plt.plot(Kr_1, Ki_1, label='Линия запаса устойчивости m=%s'%m1)
plt.plot(Kr_2, Ki_2, label='Линия запаса устойчивости m=%s'%m)
plt.plot(Kr_3, Ki_3, label='Линия  границы устойчивости m=0')
plt.xlabel("Коэффициенты - Ki")
plt.ylabel("Коэффициенты - Kr")
plt.legend(loc='best')
plt.grid(True)
alfa=1.2
Ki_1=[Ki(alfa,0.221,w) for w in n]
Kr_1=[Kr(alfa,0.221,w) for w in n]
Ki_2=[Ki(alfa,0.366,w) for w in n]
Kr_2=[Kr(alfa,0.366,w) for w in n]
Ki_3=[Ki(alfa,0,w) for w in n]
Kr_3=[Kr(alfa,0,w) for w in n]
plt.figure()
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для alfa=%s"%alfa)
plt.axis([0.0, round(max(Kr_3),3), 0.0, round(max(Ki_3),3)])
plt.plot(Kr_1, Ki_1, label='Линия запаса устойчивости m=%s'%m1)
plt.plot(Kr_2, Ki_2, label='Линия запаса устойчивости m=%s'%m)
plt.plot(Kr_3, Ki_3, label='Линия  границы устойчивости m=0')
plt.xlabel("Коэффициенты - Ki")
plt.ylabel("Коэффициенты - Kr")
plt.legend(loc='best')
plt.grid(True)
plt.figure(
)

Для запаса устойчивости m=0.386 при значениях alfa -0.2,0,4,0.7 построим D –разбиения для определения критических значения alfa.

Общая плоскость всех настроечных параметров
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для запаса устойчивости m=%s"%m)
alfa=0.2
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
alfa=0.4
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
alfa=0.7
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
alfa=1.2
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
plt.axis([0.0, round(max(Kr_2),3), 0.0, round(max(Ki_2),3)])
plt.legend(loc='best')
plt.grid(True)


Строим график зависимости Ki от частоты w для двух выявленных из предыдущего графика значений alfa: 0.2; 0.7. Настройки Ki определяются по резонансным частотам.

plt.figure()
plt.title("График Ki=f(w)")
Ki_1=[Ki(0.2,m,w) for w in n1]
Ki_2=[Ki(0.7,m,w) for w in n1]
Ky=max([round(max(Ki_1),4),round(max(Ki_2),4)])
plt.axis([0.0,round(max(n1),4),0.0, Ky])
plt.plot(n1, Ki_1,label='allfa=Td/Ti =0.2, запас устойчивости m=0.366')
plt.plot(n1, Ki_2,label='allfa=Td/Ti =0.7, запас устойчивости m=0.366')
plt.legend(loc='best')
plt.grid(True)

Определяем три варианта настроек ПИД регулятора:

Получение численных значений настроек для выбора оптимальных по предложенному критерию
maxKi=max( [Ki(0.7,m,w) for w in n1])
wa=round([w for w in n1 if Ki(0.7,m,w)==maxKi][0],3)
Ki1=round(Ki(0.7,m,wa),3)
Kr1=round(Kr(0.7,m,wa),3)
Ti1=round(Kr1/Ki1,3)
Td1=round(0.7*Ti1,3)
d=[]
d[0]= "Настройки №1 ПИД регулятора (wa =%s,m=0.366,alfa=0.7): Kr=%s; Ti=%s; Ki=%s; Td=%s "%(wa,Kr1,Ti1,Ki1,Td1)
print(d[0])
maxKi=max( [Ki(0.2,m,w) for w in n1])
wa=round([w for w in n1 if Ki(0.2,m,w)==maxKi][0],3)
Ki2=round(Ki(0.2,m,wa),3)
Kr2=round(Kr(0.2,m,wa),3)
Ti2=round(Kr2/Ki2,3)
Td2=round(0.2*Ti2,3)
d[1]= "Настройки №2 ПИД регулятора(wa =%s,m=0.366,alfa=0.2): Kr=%s; Ti=%s; Ki=%s; Td=%s "%(wa,Kr2,Ti2,Ki2,Td2)
print(d[1])
wa=fsolve(lambda w:Ki(0.7,m,w)-0.14,0.07)[0]
wa=round(wa,3)
Ki3=round(Ki(0.7,m,wa),3)
Kr3=round(Kr(0.7,m,wa),3)
Ti3=round(Kr3/Ki3,3)
Td3=round(0.7*Ti3,3)
d[2]= ("Настройки №3 ПИД регулятора(wa =%s,m=0.366,alfa=0.7): Kr=%s; Ti=%s; Ki=%s; Td=%s "%(wa,Kr3,Ti3,Ki3,Td3)
print(d[2])
def Wsys(w,Kr,Ti,Td):
         j=(-1)**0.5
         return (WO(0,w)*WR(w,Kr,Ti,Td)/(1+WO(0,w)*WR(w,Kr,Ti,Td)))
Wsys_1=[abs(Wsys(w,Kr1,Ti1,Td1)) for w in n2]
Wsys_2=[abs(Wsys(w,Kr2,Ti2,Td2)) for w in n2]
Wsys_3=[abs(Wsys(w,Kr3,Ti3,Td3)) for w in n2]
plt.figure()


Определяем динамику АСР по значению частотного показателя колебательности:

plt.title("Амплитудно-частотные характеристики замкнутой АСР \n с ПИД регулятором")
plt.plot(n2, Wsys_1,label=' Для настройки №1 ПИД регулятора')
plt.plot(n2, Wsys_2,label=' Для настройки №2 ПИД регулятора')
plt.plot(n2, Wsys_3,label=' Для настройки №3 ПИД регулятора')
plt.legend(loc='best')
plt.grid(True)

Определяем оптимальные настройки регулятора по предложенному интегральному критерию качества регулирования:



Численный и графический анализ оптимальности полученных настроек
def ReWsys(w,t,Kr,Ti,Td):
         return(2/np.pi)* (WO(0,w)*WR(w,Kr,Ti,Td)/(1+WO(0,w)*WR(w,Kr,Ti,Td))).real*(np.sin(w*t)/w)
def h(t,Kr,Ti,Td):
         return quad(lambda w: ReWsys(w,t,Kr,Ti,Td),0,0.6)[0] 
tt=np.arange(0,300,3)
h1=[h(t,Kr1,Ti1,Td1) for t in tt]
h2=[h(t,Kr2,Ti2,Td2) for t in tt]
h3=[h(t,Kr3,Ti3,Td3) for t in tt]
I1=round(quad(lambda t: h(t,Kr1,Ti1,Td1), 0,200)[0],3)
I11=round(quad(lambda t: h(t,Kr1,Ti1,Td1)**2,0, 200)[0],3)
I2=round(quad(lambda t: h(t,Kr2,Ti2,Td2), 0,200)[0],3)
I21=round(quad(lambda t: h(t,Kr2,Ti2,Td2)**2,0, 200)[0],3)
I3=round(quad(lambda t: h(t,Kr3,Ti3,Td3), 0,200)[0],3)
I31=round(quad(lambda t: h(t,Kr3,Ti3,Td3)**2,0, 200)[0],3)
print("Линейный интегральный  критерий качества I1 =%s (настройки №1)"%I1)
print("Квадратичный интегральный  критерий качества I2 =%s (настройки №1"%I11)
print("Линейный интегральный критерий качества I1 =%s (настройки №2 )"%I2)
print("Квадратичный интегральный критерий качества I2 =%s (настройки №2)"%I21)
print("Линейный интегральный критерий качества I1 =%s (настройки №3 )"%I3)
print("Квадратичный интегральный критерий качества I2 =%s (настройки №3)"%I31)
Rez=[I1+I11,I2+I21,I3+I31]
In=Rez.index(min(Rez))
print("Оптимальные параметры по интегральным критериям :\n %s"%d[In])
plt.figure()
plt.title("Переходные характеристики замкнутой АСР \n с ПИД регулятором")
plt.plot(tt,h1,'r',linewidth=1,label=' Для настройки №1 ПИД регулятора')
plt.plot(tt,h2,'b',linewidth=1,label=' Для настройки №2 ПИД регулятора')
plt.plot(tt,h3,'g',linewidth=1,label=' Для настройки №3 ПИД регулятора')
plt.legend(loc='best')
plt.grid(True)
plt.show()


Результат работы программы –текстовый вывод:

Настройки №1 ПИД регулятора (wa =0.066,m=0.366,alfa=0.7): Kr=4.77; Ti=21.682; Ki=0.22; Td=15.177
Настройки №2 ПИД регулятора(wa =0.056,m=0.366,alfa=0.2): Kr=2.747; Ti=50.87; Ki=0.054; Td=10.174
Настройки №3 ПИД регулятора(wa =0.085,m=0.366,alfa=0.7): Kr=3.747; Ti=26.387; Ki=0.142; Td=18.471
Линейный интегральный критерий качества I1 =194.65 (настройки №1)
Квадратичный интегральный критерий качества I2 =222.428 (настройки №1
Линейный интегральный критерий качества I1 =179.647 (настройки №2 )
Квадратичный интегральный критерий качества I2 =183.35 (настройки №2)
Линейный интегральный критерий качества I1 =191.911 (настройки №3 )
Квадратичный интегральный критерий качества I2 =204.766 (настройки №3)
Оптимальные параметры по интегральным критериям:
Настройки №2 ПИД регулятора(wa =0.056,m=0.366,alfa=0.2): Kr=2.747; Ti=50.87; Ki=0.054; Td=10.174

Результат работы программы, графический вывод:















Полный текст программы
# -*- coding: utf8 -*-    
import numpy as np
from scipy.integrate import quad
from scipy.optimize import *
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.family'] = 'fantasy'
mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial'
T1=14;T2=18;T3=28;K=0.9;tau=6.4# Постоянные времени, коэффициент, запаздывание
m=0.366;m1=0.221# Запас устойчивости
n= np.arange(0.001,0.15,0.0002)#Массив частот для плоскости Kr-Ki
n1=np.arange(0.00001,0.12,0.0001)#Массив частот для графика Ki=f(w)
n2=np.arange(0.0002,0.4,0.0001)#Массив частот для построения АЧХ
def WO(m,w):#Передаточная функция объекта
         j=(-1)**0.5
         return K*np.exp(-tau*(-m+j)*w)/((T1*(-m+j)*w+1)*(T2*(-m+j)*w+1)*(T3*(-m+j)*w+1))
def WR(w,Kr,Ti,Td):#Передаточная функция регулятора
         j=(-1)**0.5
         return Kr*(1+1/(j*w*Ti)+j*w*Td)
def ReW(m,w):#Действительная часть передаточной функции
          j=(-1)**0.5
          return WO(m,w).real
def ImW(m,w):#Мнимая часть передаточной функции
          j=(-1)**0.5
          return WO(m,w).imag
def A0(m,w):#Вспомогательная функция
         return -(ImW(m,w)*m/(w+w*m**2)+ReW(m,w)/(w+w*m**2))
def Ti(alfa,m,w):#Коэффициент регулятора
         return  (-ImW(m,w)-(ImW(m,w)**2-4*((ReW(m,w)*alfa*w-ImW(m,w)*alfa*w*m)*A0(m,w)))**0.5)/(2*(ReW(m,w)*alfa*w-ImW(m,w)*alfa*w*m))
def Ki(alfa,m,w):#Коэффициент регулятора
         return 1/(w*Ti(alfa,m,w)**2*alfa*(m*ReW(m,w)+ImW(m,w))-Ti(alfa,m,w)*ReW(m,w)+(m*ReW(m,w)-ImW(m,w))/(w+w*m**2))
def Kr(alfa,m,w):#Коэффициент регулятора
         if Ki(alfa,m,w)*Ti(alfa,m,w)<0:
                  z=0
         else:
                  z=Ki(alfa,m,w)*Ti(alfa,m,w)
         return z         
def Kd(alfa,m,w):#Коэффициент регулятора
         return alfa*Kr(alfa,m,w)*Ti(alfa,m,w)
alfa=0.2
Ki_1=[Ki(alfa,m1,w) for w in n]
Kr_1=[Kr(alfa,m1,w) for w in n]
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
Ki_3=[Ki(alfa,0,w) for w in n]
Kr_3=[Kr(alfa,0,w) for w in n]    
plt.figure()
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для alfa=%s"%alfa)
plt.axis([0.0, round(max(Kr_3),4), 0.0, round(max(Ki_3),4)])
plt.plot(Kr_1, Ki_1, label='Линия запаса устойчивости m=%s'%m1)
plt.plot(Kr_2, Ki_2, label='Линия запаса устойчивости m=%s'%m)
plt.plot(Kr_3, Ki_3, label='Линия  границы устойчивости m=0')
plt.xlabel("Коэффициенты - Ki")
plt.ylabel("Коэффициенты - Kr")
plt.legend(loc='best')
plt.grid(True)
alfa=0.7
Ki_1=[Ki(alfa,0.221,w) for w in n]
Kr_1=[Kr(alfa,0.221,w) for w in n]
Ki_2=[Ki(alfa,0.366,w) for w in n]
Kr_2=[Kr(alfa,0.366,w) for w in n]
Ki_3=[Ki(alfa,0,w) for w in n]
Kr_3=[Kr(alfa,0,w) for w in n]
plt.figure()
plt.axis([0.0, round(max(Kr_3),3), 0.0, round(max(Ki_3),3)])
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для alfa=%s"%alfa)
plt.plot(Kr_1, Ki_1, label='Линия запаса устойчивости m=%s'%m1)
plt.plot(Kr_2, Ki_2, label='Линия запаса устойчивости m=%s'%m)
plt.plot(Kr_3, Ki_3, label='Линия  границы устойчивости m=0')
plt.xlabel("Коэффициенты - Ki")
plt.ylabel("Коэффициенты - Kr")
plt.legend(loc='best')
plt.grid(True)
alfa=1.2
Ki_1=[Ki(alfa,0.221,w) for w in n]
Kr_1=[Kr(alfa,0.221,w) for w in n]
Ki_2=[Ki(alfa,0.366,w) for w in n]
Kr_2=[Kr(alfa,0.366,w) for w in n]
Ki_3=[Ki(alfa,0,w) for w in n]
Kr_3=[Kr(alfa,0,w) for w in n]
plt.figure()
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для alfa=%s"%alfa)
plt.axis([0.0, round(max(Kr_3),3), 0.0, round(max(Ki_3),3)])
plt.plot(Kr_1, Ki_1, label='Линия запаса устойчивости m=%s'%m1)
plt.plot(Kr_2, Ki_2, label='Линия запаса устойчивости m=%s'%m)
plt.plot(Kr_3, Ki_3, label='Линия  границы устойчивости m=0')
plt.xlabel("Коэффициенты - Ki")
plt.ylabel("Коэффициенты - Kr")
plt.legend(loc='best')
plt.grid(True)
plt.figure()
plt.title("Плоскость настроечных параметров ПИД регулятора\n  для запаса устойчивости m=%s"%m)
alfa=0.2
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
alfa=0.4
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
alfa=0.7
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
alfa=1.2
Ki_2=[Ki(alfa,m,w) for w in n]
Kr_2=[Kr(alfa,m,w) for w in n]
plt.plot(Kr_2, Ki_2,label=' Линия для allfa=Td/Ti =%s'%alfa)
plt.axis([0.0, round(max(Kr_2),3), 0.0, round(max(Ki_2),3)])
plt.legend(loc='best')
plt.grid(True)
plt.figure()
plt.title("График Ki=f(w)")
Ki_1=[Ki(0.2,m,w) for w in n1]
Ki_2=[Ki(0.7,m,w) for w in n1]
Ky=max([round(max(Ki_1),4),round(max(Ki_2),4)])
plt.axis([0.0,round(max(n1),4),0.0, Ky])
plt.plot(n1, Ki_1,label='allfa=Td/Ti =0.2, запас устойчивости m=0.366')
plt.plot(n1, Ki_2,label='allfa=Td/Ti =0.7, запас устойчивости m=0.366')
plt.legend(loc='best')
plt.grid(True)
maxKi=max( [Ki(0.7,m,w) for w in n1])
wa=round([w for w in n1 if Ki(0.7,m,w)==maxKi][0],3)
Ki1=round(Ki(0.7,m,wa),3)
Kr1=round(Kr(0.7,m,wa),3)
Ti1=round(Kr1/Ki1,3)
Td1=round(0.7*Ti1,3)
d={}
d[0]= "Настройки №1 ПИД регулятора (wa =%s,m=0.366,alfa=0.7): Kr=%s; Ti=%s; Ki=%s; Td=%s "%(wa,Kr1,Ti1,Ki1,Td1)
print(d[0])
maxKi=max( [Ki(0.2,m,w) for w in n1])
wa=round([w for w in n1 if Ki(0.2,m,w)==maxKi][0],3)
Ki2=round(Ki(0.2,m,wa),3)
Kr2=round(Kr(0.2,m,wa),3)
Ti2=round(Kr2/Ki2,3)
Td2=round(0.2*Ti2,3)
d[1]= "Настройки №2 ПИД регулятора(wa =%s,m=0.366,alfa=0.2): Kr=%s; Ti=%s; Ki=%s; Td=%s "%(wa,Kr2,Ti2,Ki2,Td2)
print(d[1])
wa=fsolve(lambda w:Ki(0.7,m,w)-0.14,0.07)[0]
wa=round(wa,3)
Ki3=round(Ki(0.7,m,wa),3)
Kr3=round(Kr(0.7,m,wa),3)
Ti3=round(Kr3/Ki3,3)
Td3=round(0.7*Ti3,3)
d[2]= "Настройки №3 ПИД регулятора(wa =%s,m=0.366,alfa=0.7): Kr=%s; Ti=%s; Ki=%s; Td=%s "%(wa,Kr3,Ti3,Ki3,Td3)
print(d[2])
def Wsys(w,Kr,Ti,Td):
         j=(-1)**0.5
         return (WO(0,w)*WR(w,Kr,Ti,Td)/(1+WO(0,w)*WR(w,Kr,Ti,Td)))
Wsys_1=[abs(Wsys(w,Kr1,Ti1,Td1)) for w in n2]
Wsys_2=[abs(Wsys(w,Kr2,Ti2,Td2)) for w in n2]
Wsys_3=[abs(Wsys(w,Kr3,Ti3,Td3)) for w in n2]
plt.figure()
plt.title("Амплитудно-частотные характеристики замкнутой АСР \n с ПИД регулятором")
plt.plot(n2, Wsys_1,label=' Для настройки №1 ПИД регулятора')
plt.plot(n2, Wsys_2,label=' Для настройки №2 ПИД регулятора')
plt.plot(n2, Wsys_3,label=' Для настройки №3 ПИД регулятора')
plt.legend(loc='best')
plt.grid(True)
def ReWsys(w,t,Kr,Ti,Td):
         return(2/np.pi)* (WO(0,w)*WR(w,Kr,Ti,Td)/(1+WO(0,w)*WR(w,Kr,Ti,Td))).real*(np.sin(w*t)/w)
def h(t,Kr,Ti,Td):
         return quad(lambda w: ReWsys(w,t,Kr,Ti,Td),0,0.6)[0] 
tt=np.arange(0,300,3)
h1=[h(t,Kr1,Ti1,Td1) for t in tt]
h2=[h(t,Kr2,Ti2,Td2) for t in tt]
h3=[h(t,Kr3,Ti3,Td3) for t in tt]
I1=round(quad(lambda t: h(t,Kr1,Ti1,Td1), 0,200)[0],3)
I11=round(quad(lambda t: h(t,Kr1,Ti1,Td1)**2,0, 200)[0],3)
I2=round(quad(lambda t: h(t,Kr2,Ti2,Td2), 0,200)[0],3)
I21=round(quad(lambda t: h(t,Kr2,Ti2,Td2)**2,0, 200)[0],3)
I3=round(quad(lambda t: h(t,Kr3,Ti3,Td3), 0,200)[0],3)
I31=round(quad(lambda t: h(t,Kr3,Ti3,Td3)**2,0, 200)[0],3)
print("Линейный интегральный  критерий качества I1 =%s (настройки №1)"%I1)
print("Квадратичный интегральный  критерий качества I2 =%s (настройки №1"%I11)
print("Линейный интегральный критерий качества I1 =%s (настройки №2 )"%I2)
print("Квадратичный интегральный критерий качества I2 =%s (настройки №2)"%I21)
print("Линейный интегральный критерий качества I1 =%s (настройки №3 )"%I3)
print("Квадратичный интегральный критерий качества I2 =%s (настройки №3)его"%I31)
Rez=[I1+I11,I2+I21,I3+I31]
In=Rez.index(min(Rez))
print("Оптимальные параметры по интегральным критериям :\n %s"%d[In])
plt.figure()
plt.title("Переходные характеристики замкнутой АСР \n с ПИД регулятором")
plt.plot(tt,h1,'r',linewidth=1,label=' Для настройки №1 ПИД регулятора')
plt.plot(tt,h2,'b',linewidth=1,label=' Для настройки №2 ПИД регулятора')
plt.plot(tt,h3,'g',linewidth=1,label=' Для настройки №3 ПИД регулятора')
plt.legend(loc='best')
plt.grid(True)
plt.show()


Выводы


Разработано программное обеспечение, для реализации в полном объёме поискового метода определения оптимальных настроек ПИД регуляторов с использованием интегрального критерия, объединяющего преимущества его линейного и квадратичного представления.
Теги:
Хабы:
+10
Комментарии 2
Комментарии Комментарии 2

Публикации

Истории

Работа

Ближайшие события

Московский туристический хакатон
Дата 23 марта – 7 апреля
Место
Москва Онлайн
Геймтон «DatsEdenSpace» от DatsTeam
Дата 5 – 6 апреля
Время 17:00 – 20:00
Место
Онлайн