Разработка математической модели и синтез системы управления барабанной сушилки

Сглаживание четвертыми разностями производится путем аппроксимации  пяти соседних значений последовательности x(n) параболой с помощью метода наименьших квадратов. В качестве элемента сглаженной последовательности y(n) принимается точка параболы, наилучшим образом аппроксимирующей значения сглаживаемой последовательности x(n) в пяти точках.

e(n)

,

где - центральная четвёртая разность.

Величину сглаженной функции  можно найти по формуле:

y(1), y(2), y(n-1), y(n) рассчитываются по формулам:

 

- центральная третья разность, рассчитываемая по формуле:

Для упрощения обработки  данных необходимо произвести сглаживание, что позволит отфильтровать высокочастотную  составляющую, которая ошибочно может  быть принята, за изменения характера  трафика.

Усреднение переходных функций  методами четвертых разностей дает лучшие результаты, потому что этот метод основан на определении, что на интервале усреднения переходная функция имеет вид параболы. Также достоинством этого метода является наличие не смещенных оценок. Недостаток – он не очень эффективный.

Рассмотрим задачу сглаживания  зашумленной переходной функции  объекта с передаточной функцией вида:

 

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

clc

clear

W=tf([0.1],[ ])

T_end=4000;% интервал измерений

dt=10;% шаг дискретизации

t=0:dt:T_end;% массив дискретного  времени

N=length(t);% размер выборки

u=ones(N,1);% моделирование единичного входного воздействие

x=lsim(W,u,t);%полезный сигнал

v=normrnd(0,1,N,1); % моделирование помехи

y=lsim(W,u,t)+v/250;% моделирование выходного воздействия с учетом аддитивной выходной помехи

plot(t,y),grid ; xlabel('t-время'); ylabel('y');

m=17; % задание числа точек  для усреднения (глубина сглаживания)

h(1)=y(1);

for i=2:m % ycpeднение начального участка

del=i-1;

h(i)=sum(y(1:i+del))/( 2*del+1);

end;

for i=m+1:N-m % основной алгоритм ycpeднения «скользящим средним»

h(i)=sum(y(i-m:i+m))/(2*m+1);

end;

for i=N-m+1:N % ycpeднение конечного участка

del=N-i;

h(i)=sum(y(i-del:N))/( 2*del+1);

end;

figure(2)

plot(t,y,t,h,'r',t,x,'g');

xlabel('t-время')

ylabel('y(t),h(t),x(t)')

grid;

 

Рисунок 3.1 – Результат сглаживания переходной

Характеристики

1 – зашумленный сигнал; 2 – идеальный сигнал; 3 – сглаженный  сигнал;

Графическое представление  результатов сглаживания переходной характеристики в условиях действия аддитивной помехи показывает (рисунок 3.1) удовлетворительное качество рассмотренного алгоритма усреднения.

3.2 Нахождение коэффициентов  переходной функции методом площадей

В основе метода площадей лежит  предположение, что объект может  быть описан линейным дифференциальным уравнением с постоянными  коэффициентами, а его нормированная (приведенная  к 1) переходная характеристика может быть аппроксимирована передаточной функцией вида:

                              (24)

На  первом  этапе  осуществляют  нормирование  переходной  характеристики и входного воздействия:

                                       (25)

Искомые коэффициенты W(p) определяются из системы уравнений:

                                                (26)

Входящие в это уравнение  коэффициенты S₁, S₂, S₃ связаны с кривой  разгона интегральными соотношениями (3.7) – относительное время.

Для расчета S₁, S₂ S₃ используют численные методы (метод прямоугольников, метод трапеций и др.):

                         (27)

                              

Зная коэффициенты S₁, S₂, S₃, находим a₁, a₂, а₃ и _, получим передаточную функцию

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

clc

clear

W=tf([0.1],[ ])

T_end=4000;% интервал измерений

dt=10;% шаг дискретизации

t=0:dt:T_end;% массив дискретного  времени

N=length(t);% размер выборки

u=ones(N,1);% моделирование единичного входного воздействие

x=lsim(W,u,t);%полезный сигнал

v=normrnd(0,1,N,1); % моделирование помехи

y=lsim(W,u,t)+v/450;% моделирование выходного воздействия с учетом аддитивной выходной помехи

m=17; % задание числа точек  для усреднения (глубина сглаживания)

h(1)=y(1);

for i=2:m % ycpeднение начального участка

del=i-1;

h(i)=sum(y(1:i+del))/( 2*del+1);

end;

for i=m+1:N-m % основной алгоритм ycpeднения «скользящим средним»

h(i)=sum(y(i-m:i+m))/(2*m+1);

end;

for i=N-m+1:N % ycpeднение конечного участка

del=N-i;

h(i)=sum(y(i-del:N))/( 2*del+1);

end;

T_end=12000;% интервал измерений

dt=10;% шаг дискретизации

t=0:dt:T_end;% массив дискретного  времени

N=length(t);% размер выборки

u=ones(N,1);% моделирование единичного входного воздействие

figure(1)

%h=lsim(W,u,t);%полезный сигнал

plot(t,h),grid

k=h(N)%найдём коэффициент  усиления(передачи)

h0=h/k;%найдём нормированное  значение переходной характеристики

f=1-h0;%подынтегральная функция

%trapz-Инетгрирование методом трапеций

A1 = trapz(t,f)%z должно быть приблизительно равно коэффициенту при s в передаточной функции

x=t/A1;%относительное время

x=x;

f2=(1-h0).*(1-x);%подынтегральная  функция

A2=A1^2*trapz(x,f2);

f3=(1-h0).*(1-2.*x+(x.^2)./2);

A3=(A1^3)*trapz(x,f3);

A=[1 0 -1;

   0 1 -A1;

   0 0 -A2];

B=[A1;A2;A3];

z=A^(-1)*B;

a1=z(1);

a2=z(2);

b1=z(3);

Wr=tf([b1*k k],[a2 a1 1])

h1=lsim(Wr,u,t);

plot(t,h,t,h1),grid

Рисунок 3.2 -Экспериментальная и расчётная переходная характеристика

2 – получнная характеристика методом площадей; 1 – реальная характеристика

 

k =  0.1000

интервал измерений: T_end=4000

шаг дискретизации: dt=10

размер выборки: N=401

Вывод: По результатам  получили следующую передаточную функцию

 

 

Изм.

Лист

№ докум.

Подпись

Дата

Лист

1

КП 04.04.06. ПЗ

Разраб.

Валенчиц А.А

Провер.

Лялько А.А.

Реценз.

 

Н. Контр.

 

 Утверд.

Лялько А.А.

Синтез оптимальной  системы управления

Лит.

Листов

5

БГТУ 2013

4 Синтез оптимальной системы управления

4.1 Постановка задачи и  описание метода синтеза

Задана система:

     (28)

Известен вид задающего  воздействия . Решение аналогично решению задачи стабилизации, при этом функция Беллмана имеет вид:

  (29)

где q(t) - вектор-функция, связанная с ;1(t) - слагаемое, гарантирующее положительную определённость функции.

Введём функцию:

      (30)

Вектор сопряженных переменных:

        (31)

Функционал качества имеет  вид:

 (32)

Получаем два дифференциальных уравнения (уравнение Риккати и  уравнение отностельно q(t)):

  (33)

Начальные условия:

(34)

 

Из условия трансверсальности  получаем:

(35)

Приравнивая два последних  уравнения, получаем:

(36)

Оптимальное управление имеет  вид:

(37)

4.2 Нахождение параметров  наблюдателя

Предположим, что задана непрерывная система с одним  выходом и одним входом, описываемая  уравнениями:

      (38)

 

Необходимо получить оценку вектора состояния системы x(t) которую мы будем обозначать x(t). В процессе оценки мы будем использовать всю доступную информацию, т.е. входной сигнал u(t), измененное значение выхода y(t) и матрица А,В,С. Устройство оценки состояния, называемое также наблюдателем состояния, имеет ту же динамику, что и сама система.

Рисунок 4.1 Структурная схема

Тогда наблюдатель можно  записать в виде:

    (39)

Отсюда характеристическое уравнение наблюдателя:

     (40)

Один из методов наблюдателя  состоит в том, чтобы сделать  его в 2-4 раза более быстродействующим , чем замкнутая система, следовательно мы можем выбрать такое характеристическое уравнение наблюдателя, которое содержит информацию о желаемом быстродействии:

  (41)

Тогда матрица G должна удовлетворять уравнению:

    (42)

Отсюда видно сходство с синтезом путем размещения полюсов, поэтому для синтеза наблюдателя  может быть использована формула  Аккермана и получим:

                                             (43)

Это уравнение позволяет  вычислить матрицу G по заданному характеристическому полиному наблюдателя и матрицам A и С.

Уравнения Риккати:

F=A's+sA-sB1(R^(-1))B1's+Q;

F=(S1B1R^(-1)B1'-A')z+QXz;

Матрицы наблюдателя состояния  определяются в виде:

F=A – GC      (44)

H=B         (45)

Уравнения Риккати:

F=A'*s+s*A-s*B1*(R^(-1))*B1'*s+Q;

F=(S1*B1*R^(-1)*B1'-A')*z+Q*Xz;

 

clc

clear all

global A B1 R Q S1  Xz

W=tf([0.1],[ ])

 [A B C D]=ssdata(W);

   Sy=ss(A,B, C, D);

   step(Sy)

     grid

   [m n]=size(A);

   R=0.5*eye(m)

   Q=0.5*eye(n)

   B1=B*ones(1,n)

   [K,S,e]=lqr(A, B1, Q, R)

   Sn=zeros(5,n);

   Z2=zeros(5,1);

S1=fsolve('rikkati_val',Sn)

K1=(R^-1)*(B1'*S1)

t(1)=0;    T=2.5;     X(1,1)=0;    X(2,1)=0;    X(3,1)=0; X(4,1)=0; X(5,1)=0;   dt=0.01;     N=T/dt ;

for i=1:N

    U(:,i)=1;

    X(:,i+1)=dt*(A*X(:,i)+B*U(:,i))+X(:,i);

    Y(i)=C*X(:,i)+D*U(:,i);

    t(i+1)=t(i)+dt;

end

  figure(2);

    plot(t(1:N),Y(1:N))

    grid

         Xz=X(:,N)/Y(N)

         eden=C*Xz

         global A B1 R Q S1  Xz

    Z1=fsolve('rikkati_val_two',Z2)

    K1=(R^(-1))*B1'*S1

    K2=(R^-1)*B1'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    p=2*e;

t(1)=0;  x(1,1)=0; x(2,1)=0; x(3,1)=0; x(4,1)=0; x(5,1)=0; dt=0.01; N=T/dt; xe(1,1)=0; xe(2,1)=0; xe(5,1)=0; xe(4,1)=0; xe(3,1)=0; % начальные условия

L=acker(A',C',p)

L=L';

for i=1:N

    u(:,i)=-K1(1,:)*xe(:,i)-K2(1,:)*Z1;

     x(:,i+1)=dt*(A*x(:,i)+B*u(:,i))+x(:,i);

    y(i)=C*x(:,i);

     xe(:,i+1)=dt*(A*xe(:,i)+B*u(:,i))+L*(y(i)-C*xe(:,i))+xe(:,i);

    t(i+1)=t(i)+dt;

    ye(i)=C*xe(:,i);

end

Ke=1/ye(N)

Xz=Xz*Ke;

         global A B1 R Q S1  Xz

    Z1=fsolve('rikkati_val_two',Z2)

    K1=(R^(-1))*B1'*S1

    K2=(R^-1)*B1'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    p=2*e;

t(1)=0;  x(1,1)=0; x(2,1)=0; x(3,1)=0; x(4,1)=0; x(5,1)=0; dt=0.01; N=T/dt; xe(1,1)=0; xe(2,1)=0; xe(5,1)=0; xe(4,1)=0; xe(3,1)=0; % начальные условия

L=acker(A',C',p)

L=L';

for i=1:N

    u(:,i)=-K1(1,:)*xe(:,i)-K2(1,:)*Z1;

     x(:,i+1)=dt*(A*x(:,i)+B*u(:,i))+x(:,i);

    y(i)=C*x(:,i);

     xe(:,i+1)=dt*(A*xe(:,i)+B*u(:,i))+L*(y(i)-C*xe(:,i))+xe(:,i);

    t(i+1)=t(i)+dt;

    ye(i)=C*xe(:,i);

end

   figure(2);

plot(t(1:N),xe(1,1:N))

xlabel('t-время')

ylabel('xe')

grid

figure(7);

plot(t(1:N),xe(4,1:N))

xlabel('t-время')

ylabel('xe')

grid

figure(8);

plot(t(1:N),xe(5,1:N))

xlabel('t-время')

ylabel('xe')

grid

figure(3);

plot(t(1:N),xe(2,1:N))

xlabel('t-время')

ylabel('xe')

grid

figure(4);

plot(t(1:N),xe(3,1:N))

xlabel('t-время')

ylabel('xe')

grid

figure(5);

plot(t(1:N),ye(1:N))

xlabel('t-время')

ylabel('ye')

grid

figure(6);

plot(t(1:N),u(1:N))

xlabel('t-время')

ylabel('u')

grid

 

В результате получим:

 

Xz =

   -0.0234

    0.5991

   -5.2755

   10.6704

   60.3497

 

 

Ke =

    1.4665

 

Z1 =

   1.0e+04 *

   -0.1263

   -0.3166

   -0.5962

   -1.9410

   -3.2716

 

K1 =

    0.5502    0.7505    0.6466    0.6634    0.4153

    0             0             0             0             0

    0             0             0             0             0

    0             0             0             0             0

    0             0             0             0             0

 

K2 =

    0.0313         0         0         0         0

    0.0313         0         0         0         0

    0.0313         0         0         0         0

    0.0313         0         0         0         0

    0.0313         0         0         0         0

 

L =

  319.4667 -455.6276   15.6255  138.8374    7.3207

Рисунок 4.2 Переходная характеристики х

Рисунок 4.3 Переходная характеристика y

Рисунок 4.4 Переходная характеристика управляющего воздействия

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

 

Заключение

В результате синтеза данной системы,  мы определили передаточную функцию объекта по каналу «изменения количества влажного материала –  изменение влагосодержания материала», в нашем случае древесных стружек, а также представили её в пространстве состояний. Получили дискретную модель системы. Провели сглаживание данных методом четвертых разностей. Определили передаточную функцию объекта по экспериментальной переходной характеристике методом площадей. Провели синтез оптимальной системы на основе квадратичного критерия качества, сформировали наблюдатель,  рассчитали коэффициенты оптимального регулятора, построили структурную схему системы.

 

Список использованных источников

1.Системы управления с  обратной связью. Ч.Филлипс, Р. Харбор, 2001.

2.Теория оптимизации систем  автоматического управления К.А.  Пупков., Н.Д. Егупов, 2004