Задачи исследования динамических систем

Листинг 1

input u,v;   //входы  блока (граничные условия)

init T=N#20; //начальные условия

const alfa=0.02, K=alfa*(N-1)^2,

x=linspace(0,1,N);

//x - массив узлов  равномерной сетки

T[1]=u; T[N]=v;

for (i=2,N-1) T'[i]=K*(T[i-1]-2*T[i]+T[i+1]);

output x[N],T[N]; //векторные  выходы блока

(комментарии следуют за двойным  слешем — "//"). Число узлов N задается в окне глобальных параметров. Для отображения распределения температуры используется типовой блок График Y от X, на входы которого подаются векторные сигналы x и T.

Рассмотрим еще один пример –  фильтрацию сигналов с использованием быстрого преобразования Фурье (БПФ) (файл \Demo\Язык_программирования\Фильтр_БПФ.mrj). Пусть исходный сигнал – сумма двух гармоник, а на него накладывается нормальный белый шум, сопоставимый по величине с исходным сигналом. Тогда алгоритм фильтра БПФ, записанный с помощью языка программирования, можно представить в виде

Листинг 2

t=linspace(0,1023,1024)/4;

x=sin(5*pi*t/128)+0.9*cos(pi*t/8);

xn=1024#0;

for (k=1,1024) xn[k]=x[k]+2*randg(0,1);

// t - массив значений  времени

// x - исходный  сигнал

// xn - зашумленный сигнал

y=fft(xn);  //прямое  БПФ

py=abs(y).*abs(y);

Porog=4e4;

// py - спектральная  плотность сигнала

// Porog - порог  для спектральной плотности

pz=1024#(0,0);

for (k=1,1024)

if py[k]<Porog then pz[k]=0+0i else pz[k]=y[k];

z=real(ifft(pz)); //обратное БПФ

// z - фильтрованный  сигнал

f=linspace(1,511,511)/256;

// массив значений  частоты

for (k=1,511) E[k]=py[k+1]*1e-4;

E_Porog=511#(Porog*1e-4);

// нормированные  значения спектр. плотности и  порога

output t[1024],xn[1024],x[1024],f[511],

       E[511],E_Porog[511],z[1024];

На выходе блока получаем массивы  значений сигналов, а также массив значений спектральной плотности. Результаты, отображенные с помощью блоков График Y от X, показаны на рис. 2.

Рис. 2. Фильтрация сигнала с помощью  БПФ: а - зашумленный сигнал (зеленая  линия) и исходный сигнал (красная  линия); б –спектральная плотность (красная линия) и порог для  спектральной плотности (синяя линяя); в – исходный сигнал (красная линия) и фильтрованный сигнал(синяя линия).

 

Как правило, сложная техническая  система является гибридной системой, поведение которой описывается  не только дифференциально-алгебраическими  уравнениями, но и логическими условиями, задающими переходы из одного дискретного состояния в другое. Для моделирования таких систем применяют специальные программные средства, такие как Simulink+Stateflow или MVS [6]. Язык программирования и большой набор дискретных, переключательных и логических элементов позволяют реализовывать гибридные модели также и в ПК «МВТУ».

Рассмотрим пример, включенный в  состав демонстрационных примеров пакета Stateflow под названием Stick-Slip Friction Demonstration. Моделируется движение бруска, прикрепленного к пружине, под действием внешней силы и с учетом сухого трения (такая модель подробно рассмотрена в [7]). Эта же модель, реализованная в ПК «МВТУ» (файл \Demo\Язык_программирования\ Сухое_трение_анимация.mrj), показана на рис. 3. В этом примере моделирование выполняется с использованием типового блока Анимация, демонстрирующего движение бруска в реальном масштабе времени.

Рис. 3. Модель движения бруска с сухим  трением.

 

Параметры модели: M – масса бруска; K – упругость пружины; Fstatic –  сила, которую нужно приложить, чтобы  сдвинуть брусок с места; Fsliding – сила трения скольжения (Fstatic ³Fsliding). Изменение состояния бруска описывается следующими логическими условиями: если брусок остановился или был в состоянии покоя и при этом суммарная сила Fsum не превышает по модулю Fstatic, то брусок остается в состоянии покоя; в противном случае он находится в состоянии движения. В соответствии с этими условиями алгоритм вычисления скорости бруска, реализованный в блоке Язык программирования, имеет вид

Листинг 3

input Fsum;

init V=0;

var Vold=V;

if (V*Vold<=0) and (abs(Fsum)<=Fstatic)

then begin V'=0; V=0 end           //покой

else V'=(Fsum-Fsliding*sign(V))/M; //движение

if goodstep then Vold=V;    

//Vold - скорость  на предыдущем успешном шаге

output V;

Здесь goodstep – системная логическая переменная, которая принимает значение true в случае успешного шага интегрирования и false в противном случае. На каждом шаге интегрирования производится оценка ошибки. Шаг считается успешным (goodstep=true), если эта оценка меньше допустимой. В противном случае (goodstep=false) происходит возврат к предыдущему шагу, после чего выполняется шаг меньшего размера. Изменение состояния (покой или движение) фиксируется только при успешном шаге. Благодаря использованию языка программирования модель оказалась более простой и наглядной, чем аналогичная модель в среде Simulink+Stateflow. Отметим, что и время счета такой модели в ПК «МВТУ» в несколько десятков раз меньше, чем в Simulink (для определения времени счета следует использовать пример без анимации – Сухое_трение.mrj).

Язык программирования ПК «МВТУ» позволяет каждому пользователю создавать собственные функции на основе стандартных конструкций языка. Такие функции могут быть объединены в библиотеки, вызов которых осуществляется оператором вида

Листинг 4

include "filename.txt";

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

Листинг 5

include "procedures.txt";

input B[M],A[N];

//массивы коэффициентов  передаточной функции

w=logspace(-1,3,200);

//массив частот  в логарифмическом масштабе

Nyquist(B,A,w,X,Y);

output X[200],Y[200]; //массивы  точек годографа

где процедура расчета точек  годографа содержится в файле procedures.txt и имеет вид

Листинг 6

procedure Nyquist(B[1],A[1],w[1],out Re[1],out Im[1])

  var z:complex,i:integer;

  for (i=1,cols(w))

  begin

    z=1j*w[i]; z=polyval(B,z)/polyval(A,z);

    Re[i]=real(z); Im[i]=imag(z);

  end;

end;

Размерности массивов B и A (M и N) задаются в окне глобальных параметров модели. Для отображения годографа следует использовать типовой блок График Y от X, на вход которого подаются векторные сигналы X и Y. Расчет производится на каждом шаге интегрирования, что позволяет наблюдать изменение годографа в процессе моделирования нестационарных и нелинейных систем.

1.3. Реализация дифференциально-алгебраических  уравнений

 

При использовании ПК «МВТУ», а  также аналогичных программных  комплексов (Simulink, VisSim и др.) наиболее просто реализуются модели, элементы которых описываются соотношениями  «вход-выход» либо «вход-состояние-выход», что соответствует явным алгебраическим зависимостям и дифференциальным уравнениям в форме Коши. Однако часто модели технических устройств описываются системами дифференциально-алгебраических уравнений (ДАУ), приведение которых к форме Коши затруднено или невозможно. Поэтому в состав библиотеки ПК «МВТУ» включен блок F(Y)=0 (решение алгебраических уравнений). Блок векторизован (вход и выход могут быть векторами), а его выход рассчитывается исходя из условия равенства нулю входного сигнала.

Система ДАУ может быть представлена в полуявной форме

(1.1)

x' = f(x, y, t), 0 = g(x, y, t)

либо в неявной форме

(1.2)

F(x', x, y, t) = 0

Схемы реализации уравнений вида (1.1) и (1.2) в ПК «МВТУ» показаны на рис. 4, где предполагается, что все линии связи – векторные. Здесь для формирования функций f и g используется блок Язык программирования, а для формирования функции F – блок Новый. При реализации уравнения (1.2) используется также блок Демультиплексор, который «расщепляет» вектор неизвестных размерности n + m на вектор x' размерности n и вектор y размерности m.

Рис. 4. Схемы реализации систем ДАУ: а – в полуявной форме (1.1); б  – в неявной форме (1.2).

 

Таким образом, в ПК «МВТУ» можно  реализовать практически любую  модель, описываемую системой ДАУ. Например, можно построить модель электрической  схемы, придав ей вид, показанный на рис. 5 (файл \Demo\Электротехника\Выпрямитель.mrj). Все элементы этой схемы являются замаскированными макроблоками, соединенными между собой векторными линиями связи. Векторные переменные имеют две компоненты: напряжение и ток. Модели некоторых элементов схемы показаны на рис. 6.

Рис. 5. Пример построения модели электрической  схемы: а – схема выпрямителя; б – эквивалентная схема трансформатора (макроблок).

Рис. 6. Модели элементов электрической  схемы (синяя линия – напряжение, красная – ток).

 

Аналогичным образом строятся модели линейного и нелинейного сопротивления, а также независимых и управляемых  источников напряжения и тока. К  числу элементов относятся узлы, которые могут быть двух видов: узлы-замыкания и узлы-ответвления. Узлы-замыкания имеют входов на 1 больше, чем выходов, а узлы-ответвления имеют выходов на 1 больше, чем входов. К узлам относятся также и заземления. В каждом узле-замыкании формируется алгебраическое уравнение из условия равенства напряжений на его входах (для заземления – из условия равенства напряжения нулю). Из этих уравнений находятся контурные токи, которые с помощью блоков В память передаются в узлы-ответвления.

Приведенный пример показывает, что  с помощью макроблоков и блоков В память, Из памяти и F(Y)=0 в ПК «МВТУ» достаточно просто реализуется контурный метод расчета электрических схем. Реализация некоторых других физических моделей затрудняется тем, что линии связи между блоками могут быть только направленными. В последующих версиях ПК «МВТУ» предполагается ввести также и ненаправленные (бинаправленные) линии связи, что, в частности, удобно для реализации узлового метода расчета электрических схем.

2. Моделирование

Когда модель построена, можно произвести ее моделирование. Предварительно задаются параметры расчета: метод и интервал интегрирования, минимальное и максимальное значения шага интегрирования, шаг  вывода результатов, точность. Дополнительно  можно установить режим масштабирования времени с заданием «множителя ускорения» модельного времени. Если задать этот множитель равным 1, то скорости протекания модельного и реального времени будут совпадать. Такой режим позволяет обеспечить синхронизацию обмена данных с внешними устройствами при моделировании в реальном времени.

Для решения дифференциальных уравнений  в ПК «МВТУ» реализованы 10 явных  и 6 неявных методов, среди которых  есть новые оригинальные методы [8, 9]. Известно, что классические явные методы (Рунге-Кутты, Адамса и др.) неэффективны при интегрировании жестких систем, поэтому для решения таких задач обычно используют неявные методы. В ПК «МВТУ» реализованы новые явные адаптивные методы [8], параметры которых автоматически настраиваются на решаемую задачу, благодаря чему они позволяют эффективно решать многие жесткие системы. Кроме этого, один из адаптивных методов (Адаптивный 5) обеспечивает точное воспроизведение быстро осциллирующих решений при наличии собственных значений якобиана вблизи мнимой оси. Результаты решения множества тестовых и прикладных задач показали, что реализованные в ПК «МВТУ» методы позволяют быстро и качественно правильно решать задачи разных типов (жесткие, осциллирующие, локально неустойчивые, разрывные). Некоторые результаты тестовых испытаний приведены в [10].

ПК «МВТУ» позволяет также эффективно решать системы ДАУ. При использовании явных методов алгебраическая подсистема решается независимо от дифференциальной, при этом можно применять один из трех методов (простых итераций, Ньютона-Рафсона, Бройдена). При использовании неявных методов алгебраическая и дифференциальная подсистемы решаются совместно, что позволяет решать системы ДАУ высших индексов. Система ДАУ (1.1) имеет индекс 1, если матрица частных производных ∂g(x, y, t) / ∂y обратима в любой точке на траектории решения. В этом случае можно аналитически либо численно исключить из уравнений вектор y, приведя таким образом систему к форме Коши. Если матрица ∂g / ∂y вырождена, то приведение к форме Коши невозможно, такие системы имеют индекс 2 и выше и называются системами ДАУ высших индексов (подробнее об индексе ДАУ см. в [11]). Простейшая схема, описываемая системой ДАУ индекса 2, показана на рис. 7. Эта схема осуществляет дифференцирование входного сигнала, обеспечивая более высокую точность, чем типовой блок Дифференцирование. Подобные схемы позволяют строить обратные модели (данная схема является обратной моделью для интегратора).

Рис. 7. Модель, описываемая системой ДАУ индекса 2 (дифференцирование  входного сигнала).

 

Системы ДАУ высших индексов часто  возникают при решении задач  механики, теории управления, электротехники и т.д. [11]. Например, электрическая схема, представленная на рис. 5, описывается системой ДАУ индекса 2. Однако Simulink, Vissim и многие другие известные системы моделирования позволяют решать ДАУ только индекса 1. ПК «МВТУ» имеет очевидное преимущество перед этими ПК, позволяя решать системы ДАУ высших индексов.