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

    0.0053

    0.0081

R =

  -20.0000    0.0011    0.0015   -0.0123

  -19.0000    0.0015   -0.0021   -0.0221

  -18.0000    0.0017    0.0007   -0.0370

  -17.0000    0.0017    0.0069   -0.0287

  -16.0000    0.0013    0.0123    0.0080

  -15.0000    0.0005    0.0074    0.0289

  -14.0000   -0.0003    0.0051    0.0470

  -13.0000   -0.0010    0.0092    0.0236

  -12.0000   -0.0018   -0.0070    0.0419

  -11.0000   -0.0019    0.0064    0.0221

  -10.0000   -0.0010   -0.0008    0.0000

   -9.0000   -0.0005    0.0004   -0.0054

   -8.0000    0.0001    0.0005    0.0018

   -7.0000    0.0011   -0.0003   -0.0124

   -6.0000    0.0031    0.0001   -0.0299

   -5.0000    0.0065    0.0005   -0.0161

   -4.0000    0.0110    0.0001   -0.0167

   -3.0000    0.0163   -0.0001    0.0021

   -2.0000    0.0261   -0.0007    0.0152

   -1.0000    0.0393    0.0001    0.0259

         0    0.0479    0.2477    0.0304

    1.0000    0.0393    0.0001    0.3341

    2.0000    0.0261   -0.0007    0.5129

    3.0000    0.0163   -0.0001    0.4239

    4.0000    0.0110    0.0001    0.2174

    5.0000    0.0065    0.0005    0.1442

    6.0000    0.0031    0.0001    0.1040

    7.0000    0.0011   -0.0003    0.0382

    8.0000    0.0001    0.0005    0.0150

    9.0000   -0.0005    0.0004    0.0121

   10.0000   -0.0010   -0.0008    0.0105

   11.0000   -0.0019    0.0064    0.0066

   12.0000   -0.0018   -0.0070    0.0154

   13.0000   -0.0010    0.0092   -0.0155

   14.0000   -0.0003    0.0051   -0.0225

   15.0000    0.0005    0.0074   -0.0224

   16.0000    0.0013    0.0123   -0.0038

   17.0000    0.0017    0.0069    0.0131

   18.0000    0.0017    0.0007    0.0341

   19.0000    0.0015   -0.0021    0.0119

   20.0000    0.0011    0.0015    0.0185 

cl =

     0.0343

     

  На рис. 2. 3 приведены результаты расчета автокорреляционной функции выходного сигнала (Covf for filtered y); автокорреляционной функции входного сигнала (Covf for prewhitened u); взаимная корреляционная функция (Correlation from u to y); импульсная характеристика (Impulse response estimate).

 

 
 

       Можно получить более подробный график импульсной характеристики, если выполнить функцию cra с одним аргументом zdan (Рис. 2. 4):

                        >> cra(zdan)

ans =

   0.0134

    0.1469

    0.2256

    0.1864

    0.0956

    0.0634

    0.0457

    0.0168

    0.0066

    0.0053

    0.0046

    0.0029

    0.0068

   -0.0068

   -0.0099

   -0.0099

   -0.0017

    0.0058

    0.0150

    0.0053

    0.0081              

 

      Необходимо отметить, что на графиках по оси абсцисс откладываются промежутки времени τ = t_i – t_i-1, а по оси ординат значения корреляционных функций для входного u2 и выходного у2 сигналов; значения взаимокорреляционой функции и импульсной характеристики. Из полученных характеристик следует, что с увеличением τ наблюдается резкий спад корреляционной зависимости входного сигнала, что свидетельствует о слабой взаимосвязи между сечениями процесса, соответствующими произвольным моментам времени (процесс более близок к белому шуму, а автокорреляционная функция к дельта-функции). Выходная величина наоборот более плавно изменяет свои состояния от одного момента времени к другому и, следовательно, взаимосвязь между предыдущим и последующим значениями выходного сигнала более тесная, чем у входного.

        Для получения частотных характеристик  экспериментальных данных воспользуемся  функциями оценивания частотных  характеристик. Функция spa возвращает частотные характеристики одномерного объекта и оценки спектральной плотности его сигналов для обобщенной линейной модели объекта:

[g, phiv] = spa(z)

[g, phiv, z_spe] = spa(z,M,w,maxsize,T)

Аргументы:

• z – матрица исходных данных;
• М – ширина временного окна (по умолчанию М = min(30, length(z) /10), где length(z) – число строк матрицы z);
• w – вектор частот для расчета частотных характеристик (по умолчанию              [1: 128]/128*pi/T);
• Т – интервал дискретизации;
• maxsize – параметр, определяющий максимальный размер матриц, создаваемых в процессе вычислений.

Возвращаемые  величины:

• g – оценка W(e ^jωT) в частотном формате;
• phiv – оценка спектральной плотности шума v(t);
• z_spe – матрица спектральных плотностей входного и выходного сигналов.

      Построим диаграмму Боде (АЧХ, ФЧХ), используя функции spa и bodeplot и данные, полученные при изучении теплового объекта, и содержащиеся в файле dryer2

>> load project24;

             >> z=[y2 u2];

            >> g=spa(z);

>> bodeplot(g)

      Результаты моделирования (АЧХ построена в логарифмическом масштабе) без доверительного коридора представлены на рис. 2. 5.

             

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

      Для получения частотных характеристик  вместе с доверительным коридором  шириной в три среднеквадратических  отклонения в пакете System Identification Toolbox MATLAB имеется следующие возможности:

• установка  границ изменения частот с помощью команды

>> w=logspace(w₁,w₂,N),

где w₁ – нижняя граница диапазона частот (10^w1), w₂ – верхняя граница диапазона частот (10^w2) и N – количество точек графика.

• построение АФХ, ФЧХ и S(ω) – функции спектральной плотности шума e(t)
• вычисление g – оценки АФХ и ФЧХ в частотном формате и phiv – оценки спектральной плотности шума с помощью команды

>> [g,phiv]=spa(z,[],w);

      Графики АФХ, ФЧХ и S(ω) строятся с доверительным коридором в три среднеквадратических отклонения с помощью команды

>> bodeplot([g p],'sd',3,'fill'),

где 'sd' – указывает на сплошную линию доверительного коридора (по умолчанию эта линия штриховая); 3 – величина доверительного коридора в три среднеквадратических отклонения; 'fill' – способ заливки доверительного коридора (желтым цветом).

       Построим АЧХ, ФЧХ, используя  функции spa, bodeplot, logspace и данные, полученные в файле Project24 с соответствующим доверительным коридором:

        >> w = logspace(-2,pi,128);

   >> [g,phiv]=spa(z,[],w);

   >> bodeplot([g,phiv],3,'fill')

           Результаты моделирования представлены на рис. 2. 6.

       Для построения графика оценки спектральной плотности шума с доверительным коридором выполним следующую команду:

>> bodeplot([phiv],'sd',3,'fill')

       Результаты моделирования представлены  на рис. 2. 7.

       Полученный график оценки спектральной плотности шума с доверительным коридором показывает наличие равномерного распределения мощности сигнала по частотному спектру с последующим спадом мощности на частоте выше 1, 1 рад/с.

        Далее необходимо выполнить параметрическое  оценивание ТОУ.

2. 8. Параметрическое оценивание данных

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

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

              >> zdanv=zdan(1:500);

>> zdane=zdan(501:1000).

        Первая часть массива данных  будет использоваться для параметрического  оценивания и построения модели  ТОУ. Вторая часть необходима  будет для верификации (проверки  качества) модели, определения адекватности  полученной модели и определения  погрешностей идентификации. Необходимо  отметить, что параметрическая идентификация  в пакете System Identification Toolbox MATLAB выполняется в дискретном виде и полученные модели, являются дискретными.

        В пакете System Identification Toolbox рассмотрены различные виды моделей ТОУ,  которые с различной степенью достоверности описывают объект автоматизации. Для выбора наиболее приемлемой структуры и вида моделей при параметрическом оценивании экспериментальных данных в пакете System Identification Toolbox MATLAB имеются специальные функции

• параметрического оценивания,
• задания структуры модели,
• изменения и уточнения структуры модели и выбора структуры модели.

      Функция оценивания ar оценивает параметры модели авторегресии:

A(z) y(t) = e(t),

где: A(z) = 1 + a₁z ^{– 1} + a₂z ^{– 2} +...+ a _naz ^{– na} , т.е. коэффициенты полинома A(z), при моделировании скалярных временных последовательностей. Функция имеет синтаксис:

th = ar(y,n)

       Или другое написание, позволяющее  изменять параметры моделирования: 

[th,refl]=ar(y,n,approach,win,maxsize,T)

где аргументы:

        y – вектор-столбец данных, содержащий N элементов;

        n – порядок модели (число оцениваемых коэффициентов);

        approach – аргумент (строковая переменная) определяет метод оценивания:

    • 'ls' – метод наименьших квадратов;

    • 'yw' – метод Юла-Уокера;