Software & Systems

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

Для вычисления степени затухания переходных процессов в simulink-моделях удобно использовать стандартные функции MATLAB. Запуск моделей при этом осуществляется из М-файлов. Алгоритм примерно следующий:

–      запуск модели и передача в рабочую среду MATLAB массива выходного сигнала модели;

–      обработка данного массива с использованием стандартных функций MATLAB и вычисление необходимых показателей;

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

Получение массива в рабочую среду MATLAB и построение переходного процесса

Рассмотрим простейшую модель замкнутой системы (рис. 1).

Передаточная функция объекта (блок Transfer Fcn) имеет вид: .

Для регулятора PID введем следующие настроечные параметры: Proportional=2,1, Integral=2, Derivative=1.

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

Для этого щелкнем по блоку и на вкладке Data history в поле Variable name введем какую-либо переменную, например G. В раскрывающемся списке Format выберем значение Array.

Запустим модель. Теперь, если в командном окне MATLAB ввести переменную G, выведется двумерный массив (первый столбец – время, второй – значение на выходе). Фрагмент массива приведен в таблице.

Фрагмент выходного массива G

Чтобы построить график переходного процесса, необходимо из полученного массива выделить первый и второй столбцы. Для этого напишем небольшой М-файл:

% устанавливаем путь к модели

cd 'D:\2012\ASUTP\SZ'

% открываем модель

open sz;

% запускаем модель

sim sz;

% получаем первый столбец массива

t=G(:,1);

% получаем второй столбец массива

y=G(:,2);

% строим график

plot(t,y);

% наносим координатную сетку

grid on;

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

Вычисление степени затухания (способ 1)

Чтобы вычислить степень затухания, необходимо определить первый и второй максимумы переходного процесса [1]. Готовой функции, определяющей два экстремума, в MATLAB нет, поэтому выполним следующее.

Добавим к нашему М-файлу фрагмент:

% определяем значение первого максимума

ma1=max(y);

% определяем значение минимума в переменной mi1 и номер точки, в которой он находится в переменной i

[mi1,i]=min(y);

% формируем массив из номеров точек, начиная с точки i и до конца массива

x=[i:length(y)];

% формируем массив y1 из значений ординат в полученных точках, в результате получаем хвост кривой от первого минимума и до конца

y1=y(x);

% определяем максимум в получившемся массиве

ma2=max(y1);

% вычисляем степень затухания

SZ1=(ma1-ma2)/ma1;

Выполнив М-файл и набрав в командном окне SZ1, получим значение степени затухания 0,8152.

Визуализация вычисленного значения степени затухания

Добавим в представленную модель блоки Constant и Display, как показано на рисунке 3.

В настройках блока Constant введем переменную SZ1.

Изменим М-файл следующим образом:

cd 'D:\2012\ASUTP\SZ'

open sz;

SZ1=1;

sim sz;

t=G(:,1);

y=G(:,2);

plot(t,y);

grid on;

ma1=max(y);

[mi1,i]=min(y);

x=[i:length(y)];

y1=y(x);

ma2=max(y1);

SZ1=(ma1-ma2)/ma1;

sim sz;

Как видно из М-файла, до первого запуска модели присваиваем переменной SZ1 значение, равное единице. После вычисления степени затухания присваиваем переменной SZ1 вычисленное значение и повторно запускаем модель (sim sz).

Вычисление степени затухания (способ 2)

При вычислении степени затухания рассмотренным выше способом следует иметь в виду, что переходный процесс сам по себе должен быть затухающим. Если переходный процесс расходящийся, вычисленное значение будет неверным. Изменим настроечный параметр PID-регулятора Integral=0,5.

Запустим М-файл. В результате получим расходящийся переходный процесс, приведенный на рисунке 4.

Значение степени затухания составляет 1,69, что невозможно, так как второй максимум больше первого и значение должно быть отрицательным. Это происходит потому, что максимальное значение находится в конце переходного процесса (1,1257), как и минимум (–1,4173). В оставшемся после минимума фрагменте максимальное значение составляет –0,7872. В результате получим .

Второй способ заключается в формировании массива разностей с использованием функции diff(y). Экстремальные значения будут находиться в точках, где полученные разности меняют знак на противоположный. Нам достаточно первых трех экстремумов, которые можно определить в цикле.

Добавим к М-файлу следующий фрагмент:

% получаем массив разностей

y1=diff(y);

c=1; j=0;

% определяем длину массива

n=length(y1)

for i=1 :1:n

  if c>0

    if y1(i)<0

      % добавляем в массив d элемент при смене знака

      c=c*(-1); j=j+1;d(j)=i;

    end

  end

  if c<0

    if y1(i)>0

      c=c*(-1); j=j+1;d(j)=i;

    end

  end

  if j³3

    % завершаем цикл, если количество элементов равно трем

    break

  end

end

SZ1=(y(d(1))-y(d(3)))/y(d(1));

sim sz1md;

В результате получим массив d из номеров экстремумов: d=[112  235  358].

Сами экстремумы будут иметь значения: Y(112)=0,2602, Y(235)= –0,3485, Y(358)=0,4429.

В итоге определится степень затухания, равная –0,7022, что соответствует действительности.

Получение массива экстремумов средствами simulink (способ 3)

Блок Triggered To Workspace формирует массив экстремумов и передает их в рабочую область MATCAD. Один из входов блока подключается к выходному сигналу, а управляющий вход подключается через блок Detect Decrease.

Соберем модель, приведенную на рисунке 5.

Для настройки блока Trigger To Workspace кликнем по блоку.

В раскрывающемся списке Triggered Type выберем тип экстремумов – в нашем случае Rising edge (максимальные значения). В поле Variable name введем имя переменной, значения которой будут передаваться в рабочую область MATLAB.

Запустим модель и введем в командной строке G2.

В результате получим массив экстремумов G2= и приблизительно такие же, как в первом примере, значения максимумов: Ma1=0,2901, Ma2=0,0564.

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

cd 'D:\2012\ASUTP\SZ'

open sz;

SZ1=1;

sim sz;

SZ1=(G2(1)-G2(2))/G2(1);

sim sz;

Полученная таким образом степень затухания составляет 0,8057, а в первом примере она составляла 0,8058.

Заметим, что в командном окне MATLAB выводится предупреждение, что блок Detect Decrease требует на входе дискретный сигнал.

Чтобы устранить данное предупреждение, подадим на вход блока Detect Decrease сигнал через блок Zero Order Hold, который осуществляет дискретизацию аналогового сигнала. Создадим модель, приведенную на рисунке 6.

Визуализация экстремумов

Для визуализации экстремумов можно использовать блок Maximum.

Создадим модель, приведенную на рисунке 7.

Оба входа данного блока подключаются через блок Zero Order Hold. Для настройки блока Maximum кликнем по блоку. В поле Mode введем значение Running, а в поле Reset ports – значение Rising elge.

После запуска модели на осциллографе Scope1 получим значения максимумов, как показано на рисунке 8.

Объединим все рассмотренные способы в одной модели, представленной на рисунке 9.

Напишем М-файл для управления данной моделью:

cd 'D:\2012\ASUTP\SZ'

open sz;

SZ1=1; SZ2=1; SZ3=1; Kp=2; Ti=0.5;

sim sz;

% построение графика

t=G(:,1); y=G(:,2); plot(t,y); grid on;

% способ 1

[ma1]=max(y); [mi1,i]=min(y); x=[i:length(y)];

y1=y(x); ma2=max(y1); SZ1=(ma1-ma2)/ma1;

% способ 2

y1=diff(y); c=1; j=0; n=length(y1)

for i=1 :1:n

  if c>0

    if y1(i)<0 c=c*(-1); j=j+1;d(j)=i; end

  end

  if c<0

    if y1(i)>0 c=c*(-1); j=j+1;d(j)=i; end

  end

  if j³3

    break

  end

end

SZ2=(y(d(1))-y(d(3)))/y(d(1));

% способ 3

SZ3=(G2(1)-G2(2))/G2(1);

sim sz;

Как видно из рисунка 9, при настроечном параметре регулятора Integral=0,5 получим одинаковую степень затухания переходного процесса, вычисленную тремя различными способами.

Если изменить настроечный параметр Inte­gral=1, получится расходящийся переходный процесс. В этом случае одинаковые степени затухания будут только при расчете вторым и третьим способами, как показано на рисунке 10.

В результате можно сделать следующие выводы.

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

Литература

1.     Мерзлякова О.В., Бильфельд Н.В. Получение массива в рабочую среду MATLAB и построение переходного процесса // Молодежная наука в развитии регионов. 2011. С. 14–20.

2.     Ануфриев И. Самоучитель MatLab 5.3/6.x. СПб: БХВ-Петербург, 2002. 512 с.

3.     Черных И.В. Моделирование электротехнических устройств в MATLAB, SymPowerSystem и Simulink. СПб: Питер, 2008. 288 с.

4.     Затонский А.В. Программная оптимизация настроек систем автоматического регулирования // Новый университет, 2011. № 1. С. 28–33.