Для преобразования математической модели в дискретное время использовалась функция программного пакета Matlab c2d. При этом шаг дискретности нужно выбирать с учетом того что процессы в замкнутой системе будут проходить в 10 раз быстрее чем в объекте.
dt=0.01/max(abs(eig(A)))
t=0:dt:999;
[Ad,Bd]=c2d(A,B,dt);
dt=0.4500
Проверить найденную модель в дискретном времени следует с помощью расчета разгонных характеристик. Для этого следует использовать функцию dstep. Для вывода графиков следует использовать функции: subplot, plot, grid.
Ad =
Columns 1 through 8
0.9872 0 0 0 0 0 0 0
0 0.9965 0 0 0 0 0 0
0 0.0371 0.9814 0 0 0 0 0
0 0 0 0.9882 0 0 0 0
0 0 0 0.1892 0.9048 0 0 0
0 0 0 0 0 0.9959 0 0
0 0 0 0 0 0.0134 0.9933 0
0 0 0 0 0 0 0 0.9672
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
Columns 9 through 14
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0.9954 0 0 0 0 0
0.0135 0.9933 0 0 0 0
0 0 0.9910 0 0 0
0 0 0.1439 0.9277 0 0
0 0 0 0 0.9976 0
0 0 0 0 0.0119 0.9940
Bd =
-0.5429 0
-0.2525 0
-0.0047 0
20.2483 0
1.9628 0
2.9803 0
0.0200 0
0 0.0326
0 0.0021
0 0.0000
0 -0.0618
0 -0.0045
0 -0.0123
0 -0.0001
Построим разгонные характеристики с помощью функций dstep, subplot, plot, grid.
Рисунок 4.Кривые разгона.
В результате анализа кривых разгона можно сделать вывод, что значения полученные на выходе каналов регулирования описанных инерционными звеньями 1-го порядка совпадают со значением коэффициента К инерционного звена, а на выходе каналов регулирования представленных интегрирующим звеном, кривые разгона направлены в отрицательную сторону, если имеют знак «-» в передаточной функции звена и наоборот. Если сравнить матрицу передаточных функций и полученные разгонные характеристики, видно, что Кр совпадают, можно сделать вывод: построение модели и преобразование выполнены верно.
4 Синтез многомерного ПИ-регулятора
Для синтеза ПИ-регулятора полученные матрицы должны быть расширены в матрицы A1, B1, C1:
A1=[Ad zeros(n,l); C eye(l)];
B1=[Bd;zeros(m)];
C1=[C eye(l)];
Матрицы параметров регулятора должны быть расчитаны с помощью функции dlqr.
K=dlqr(A1,B1,Q,R)
L=dlqr(A1',C1',Q1,R1)'
Весовые матрицы Q1,R1,Q,R выбраны как единичные (для простоты матрицы генерирует функция eye).
Матрицы имеют вид:
A1 =
Columns 1 through 8
0.9872 0 0 0 0 0 0 0
0 0.9965 0 0 0 0 0 0
0 0.0371 0.9814 0 0 0 0 0
0 0 0 0.9882 0 0 0 0
0 0 0 0.1892 0.9048 0 0 0
0 0 0 0 0 0.9959 0 0
0 0 0 0 0 0.0134 0.9933 0
0 0 0 0 0 0 0 0.9672
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
1.0000 -1.0000 1.0000 0 0 0 0 1.0000
0 0 0 -1.0000 1.0000 -1.0000 1.0000 0
Columns 9 through 16
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0.9954 0 0 0 0 0 0 0
0.0135 0.9933 0 0 0 0 0 0
0 0 0.9910 0 0 0 0 0
0 0 0.1439 0.9277 0 0 0 0
0 0 0 0 0.9976 0 0 0
0 0 0 0 0.0119 0.9940 0 0
-1.0000 1.0000 0 0 0 0 1.0000 0
0 0 -1.0000 1.0000 -1.0000 1.0000 0 1.0000
B1 =
-0.5429 0
-0.2525 0
-0.0047 0
20.2483 0
1.9628 0
2.9803 0
0.0200 0
0 0.0326
0 0.0021
0 0.0000
0 -0.0618
0 -0.0045
0 -0.0123
0 -0.0001
0 0
0 0
C1 =
Columns 1 through 13
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
Columns 14 through 16
0 1 0
0 0 1
K =
Columns 1 through 8
0.0367 -0.0578 0.0407 0.0634 0.0239 -0.0813 0.1013 0.0485
21.0412 -24.2138 21.6345 4.5472 11.1495 -21.5375 25.0390 22.7474
Columns 9 through 16
-0.0419 0.0319 0.0446 0.0349 -0.0865 0.1026 0.0000 0.0001
-21.8436 20.3135 0.9756 13.2017 -22.4572 25.2658 0.0195 0.0270
L =
0.0925 -0.0000
0.1180 0.0000
0.3752 0.0000
0.0000 0.0568
0.0000 0.1807
0.0000 0.0279
-0.0000 0.2379
0.0694 -0.0000
0.0275 0.0000
0.2971 0.0000
-0.0000 0.0629
-0.0000 0.1964
0.0000 0.0673
0.0000 0.3243
1.6139 0.0000
0.0000 1.6702
5 Моделирование замкнутой системы и оценка качества переходных процессов
Рисунок 5 – Структурная схема системы в виде переменных состояния с учетом запаздывания.
Для получения переходных процессов следует сформировать матрицы замкнутой системы и получить переходные процессы с помощью программы dstep .
AA=[A1 -B1*K;L*C1 A1-B1*K-L*C1];
BB=[B1;zeros(n+l,m)];
CC=[C zeros(l) zeros(l,n+l)];
При оценке качества переходных процессов необходимо чтоб управляющее воздействие не превышало 100% открытия.
Максимальное возмущение следует принять на уровне 10% номинального значения соответствующих параметров. Допустимое значение урегулированных переменных нужно принять равными 20% номинального значения. Если качество не соответствует нужно сменить весовые матрицы и повторить расчет.
Рисунок 6.Переходные процессы замкнутой системы.
Проанализировав переходные процессы можно сделать вывод, что значения регулируемых параметров не превышают допустимых. Для определения времени регулирования нужно количество итераций цикла умножить на шаг:
.6 Преобразование модели регулятора в форму, отвечающую ее реализации в програмном обеспечении
ПИ закон регулирования вычисляется по формулам:
K1=K(:,1:14);
K2=K(:,15:16);
L1=L(1:14,:);
L2=L(15:16,:);
Ar=[Ad-Bd*K1 -Bd*K2-L1 L1; C eye(2)-L2 L2; zeros(2,14) zeros(2) eye(2)];
Br=[zeros(14,2); zeros(2); eye(2)];
Cr=[-K zeros(2)];
Az=[Ad Bd*Cr; Br*C Ar];
Bf=[Bd; zeros(18,2)];
Bz=[zeros(14,2); Br];
Cz=[C zeros(2,18)];
и записать замкнутую систему в вид в котором она будет реализована в программном обеспечении:
x=zeros(14,1);xr=zeros(18,1); u=zeros(2,1);
yy=[]; uu=[];f=[.0010;.0010];z=[0;0];
for i=1:2000,
y=C*x; e=-z+y;
u=Cr*xr; xr=Ar*xr+Br*e;
y=C*x; x=Ad*x+Bd*(u+f);
yy=[yy; y']; uu=[uu; u'];
end
x1=x; xr1=xr; u1=u;
В результате выполнения программного кода будут получены переходные процессы изменения возмущения, которое поступает на каждый из каналов регулирования и переходные процессы на выходе системы.
7 Выбор технических средств реализации системы управления
Технические средства реализации системы правления включают датчики ругулированых параметров, исполнительные механизмы и регулирующие органы, преобразователи,рабочая станция
Общая структурная схема рабочей станции изображена на рисунке
Рабочая станция имеет вид:
Рисунок 7 – Схема рабочей станции
Таблица 3
Спецификация технических средств:
№ | Тип | К-во | Предназначение |
1б,2б | ADAM 4012 | 2 | Модуль аналогового ввода с датчиков давления, спиртометра, тип входного сигнала: mV, V или mA, диапазон: ±150мВ, ±20мА, ±5В, ±10В |
ADAM 4561 | 1 | Преобразователь интерфейса USB в RS-232/422/485 | |
1а | Сапфир-22М-ДА2050 | 1 | Датчик давления в магистрали, верхние пределы измерений: 1.6 МПа, Предел допускаемой основной погрешности, 0,5 0,25; 0,5%, вых. 0-5; 4-20; 5-0; 20-4мА, |
2а | Alcolyzer Plus Spirits | 1 | Спиртомер для крепких спиртных напитков; Диапазон измерения: Спирт: от 35 до 65 об. % (значения отображаются от 0 до 90 об. % спирта, однако, при содержании спирта менее 35 об. % и более 65 об. % точность измерения уменьшается); Значение рН (опционально): от 0 до 14; Цвет (опционально): от 0 до 120 EBC; Плотность (опционально): от 0 до 3 г/см3; |
2д,3д | ADAM 4069 | 2 | Модуль c релейными выходами, 8 реле с нормально разомкнутым контактом, нагрузочная способность контактов: 250 В/ 5 A для перем. тока, 30 В/ 5 A для пост. тока, время включения 5 мс , время выключения 5,6 мс |
2г,2в, 3в,3г | МЭО 40/25-0,25 | 4 | Механизм исполнительный одно-оборотный , номинальный крутя-щий момент 40кгс/м, номинальный ход выходного органа 0,25 оборота за 25с, Напряжение питания 220В. Частота 50Гц |
2г,2в | КРП-100 | 2 | Клапан регулирующий, Ду=50мм |
3г,3в | КРП-100 | 2 | Клапан регулирующий, Ду=200мм |
Вывод
В курсовом проекте было выполнено математическое обеспечение АСУТП колонны перегонки коньячного спирта. Были построены характеристики объекта: кривые разгона и при 10-и процентном возмущении. Результаты показали, что качество отвечает требуемому. В результате выбора технического обеспечения: разработана функциональная схема автоматизации, подобрано оборудование для технической реализации данной системы. Разработано программное обеспечение: программа, которая моделирует поведение системы. При тестировании данной программы было показано регулятор работает адекватно.
Литература
1.Стопакевич А.А.Теория систем и системный анализ.Учебник для вузов.-Киев: ВИПОЛ,1996.-200с.
2.Демченко В.А.Автоматизация и моделирование технологичных процессов АЭС и ТЭС.-Одесса:Астроприт,2001.-308с.
3.Стопакевич А.А.Матлаб. Методические указания к лабораторным работам,курсового и дипломного проектирования.-Одесса,2000.-18с.
Приложение 1
Текст программы
A=[-1/35 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 -1/129 0 0 0 0 0 0 0 0 0 0 0 0;
0 4/48 -2/48 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 -1/38 0 0 0 0 0 0 0 0 0 0;
0 0 0 4/9 -2/9 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 -1/110 0 0 0 0 0 0 0 0;
0 0 0 0 0 4/134 -2/134 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 -1/13.5 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 -1/98 0 0 0 0 0;
0 0 0 0 0 0 0 0 4/133 -2/133 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 -1/50 0 0 0;
0 0 0 0 0 0 0 0 0 0 4/12 -2/12 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 -1/186 0;
0 0 0 0 0 0 0 0 0 0 0 0 4/150 -2/150];
B=[-42.5/35 0;
-72.5/129 0;
0 0;
1720/38 0;
0 0;
730/110 0;
0 0;
0 0.994/13.5;
0 0.459/98;
0 0;
0 -6.9/50;
0 0;
0 -5.1/186;
0 0];
C=[1 -1 1 0 0 0 0 1 -1 1 0 0 0 0;
0 0 0 -1 1 -1 1 0 0 0 -1 1 -1 1];
D=[0 0;0 0];
dt=0.1/max(abs(eig(A)));
t=0:dt:999;
G=length(t);
[Ad Bd]=c2d(A,B,dt);
y=dstep(Ad,Bd,C,D,1,G);
figure(1)
subplot(2,2,1);plot(y(:,1));grid;ylabel('u1,МПа');title('Razgon u1,1%');
subplot(2,2,3);plot(y(:,2));grid;ylabel('u2,');
%u2
y=dstep(Ad,Bd,C,D,2,G);
subplot(2,2,2);plot(y(:,1));grid;ylabel('y1,М3/с');title('Razgon u2,1%');
subplot(2,2,4);plot(y(:,2));grid;ylabel('y2,кг/с');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A1=[Ad zeros(14,2);C eye(2)];
B1=[Bd; zeros(2)];
C1=[zeros(2,14) eye(2)];
Q2=1e6*[1 0;0 1];
V=C'*Q2*C;
Q=[V zeros(14,2);zeros(2,14) eye(2)];
R=eye(2);
Q1=eye(16);
R1=eye(2);
K=dlqr(A1,B1,Q,R);
L=dlqr(A1',C1',Q1,R1)';
K1=K(:,1:14);
K2=K(:,15:16);
L1=L(1:14,:);
L2=L(15:16,:);
Ar=[Ad-Bd*K1 -Bd*K2-L1 L1; C eye(2)-L2 L2; zeros(2,14) zeros(2) eye(2)];
Br=[zeros(14,2); zeros(2); eye(2)];
Cr=[-K zeros(2)];
Az=[Ad Bd*Cr; Br*C Ar];
Bf=[Bd; zeros(18,2)];
Bz=[zeros(14,2); Br];
Cz=[C zeros(2,18)];
x=zeros(14,1);xr=zeros(18,1); u=zeros(2,1);
yy=[]; uu=[];f=[.0010;.0010];z=[0;0];
for i=1:2000,
y=C*x; e=-z+y;
u=Cr*xr; xr=Ar*xr+Br*e;
y=C*x; x=Ad*x+Bd*(u+f);
yy=[yy; y']; uu=[uu; u'];
end
x1=x; xr1=xr; u1=u;
figure(2)
subplot(2,2,1); plot(yy(:,1));grid;ylabel('y1,MPa');title('Perehod proces Braga ');
subplot(2,2,3); plot(yy(:,2));grid;ylabel('y2,');
subplot(2,2,2); plot(uu(:,1));grid;ylabel('u1,M3/c');title('Perehod proces Par ');
subplot(2,2,4); plot(uu(:,2));grid;ylabel('u2,kg/c');