Элементы математического моделирования в программных средах MATLAB 5 и Scilab (Андриевский Фрадков)

объектов (манипуляционных роботов, летательных аппаратов и т.д.). Обозначим входное напряжение через u(t)y выходное - через 2/(f), токи в ветвях - через г ' г г г ( 0 » г з ( 0 - этой цепи составим в символической форме уравнения Кирхгофа, используя оператор дифференцирования р = -jjj. Получим

у = Л1Ч + Lipii + щл г, < - ^ - L p i 3 + y = u!

1ргъ + Riz +	= О,
г2 = i\ +г3.

Начнем с исследования системы (3.8) в общем виде. Используя очевидные обозначения для переменных, составим следующую MATLAB-программу для получения из (3.8) операторного уравнения относительно выходного сигнала y(t):

S = s o l v e ( , y = R l * i l + L l * p * i l + l / C l / p * i l ' ,:

'-R*i3-L*p*i3+y=u',...

'L*p*i3+R*i3+l/C/p*i2=0',:

'i 2 = i l + i 3 1 , ' у ' , ' i l ' , ' i 2 ' , ^ З 1 ) y=simple(S.y);

y=collect(y,'p')

latex(y)

WyssimpleCy/'u');

[n,d]=numden(Wy);

n=collect(n);

n=simple(n)

d=collect(d);

Wy=n/d

Оператор collect служит для упрощения выражения приведением подобных членов относительно указанного вторым параметром символьного аргумента. Остальные операторы рассмотрены в п. 2.2. В результате выполнения программы, в формате M g X получается выражение 1

ujl + Ljp'Cj + Rt	С,р){Ьр7С+ ДСр + 1)	3.9
	Л(р)

1 Знаменатель этой дроби выделен в многочлен А(р) "вручную", чтобы выражение уместилось по ширине страницы.

61

где А(р) = L, р4С, LC+	(Л, С, 1С + Lt	С, RC) pz + (Я, Ct RC+
+ Ct L + L, Ct + LC)p2	+ (С, Л + ЯС +	C i ) p + 1. Последую-

щие операторы служат для получения выражения для передаточной функции W(p) от входного сигнала и к выходу у. В программе эта передаточная функция обозначается идентификатором Wy. Фактически требуемое выражение получается после выполнения оператора Wy=simple(y/'и'). Последующие операторы оказываются полезными для упрощения выражения, приведения его к стандартному виду рациональной функции с упорядоченными по степеням р многочленами в числителе и знаменателе. Лля этого из найденной дроби с помощью оператора numden выделяются числитель п и знаменатель d,B которых затем производится приведение подобных членов оператором collect. Кроме того, для данной задачи удается получить аналитически разложение числителя

на множители с помощью оператора simple. 1

В результате

для числителя получаем выражение

n= (l+Ll*p~2+Rl*Cl*p)* (L*p~2*C+R*Op+l).

Окончательно

оператором Wy=n/d находится

выражение

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

Таким

образом, фильтр (3.8) описывается

передаточной

функцией

а0р* + ахр3 + а2р2 + а3р + 1

v

;

где

постоянные

времени Т}ТХ определяются

выражениями

Т

=

у/LC,

Тх

=

\/LxCx\ относительные коэффициенты

демп-

фирования

- выражениями f

=

=

а коэффициенты знаменателя - выражениями

а0 = LLXCCX,

ах

=

(LRXC1

+ LXRCX)C> а2 = RRXCCX

+ LCX

+ LXCX + LC,

а3

= RCX + RC + RXCX. Найденные выражения используем для

расчета временных и частотных характеристик

фильтра.

Расчет статического режима для данной системы выпол-

няется просто:

как известно, для

определения

статическо-

го

коэффициента

передачи

достаточно

в операторной

фор-

ме

записи

положить р =

0. Поэтому

используем

оператор

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

62

yO=simple(subs(y, 'p' ,0)). Как и следовало ожидать, получаем у0 =и. Следовательно, статический коэффициент.передачи фильтра равен единице.

Перейдем к численному исследованию динамики фильтра.

Примем следующие значения				параметров [10]:		Я = 1 2 кОм,
Rx = 22 кОм, L = 103		Гн, ЬХ = Ю3		Гн, С = 0.172 мФ, Сх =0.078 мФ.
Получим частотные		характеристики			фильтра,	которые вы-
ведем в виде диаграмм Воде			(логарифмическая			амплитудно-
частотная и фазочастотная			характеристики, ЛАХ и ЛФЧХ)
[3, 6,	75].

Продолжим программу. Подставим в полученное выражение для передаточной функции (3.10) числовые значения па-

раметров. Затем выполним переход к частотной	передаточ-
ной функции, произведя замену аргумента р на ju.	Послед-

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

63

Рис. 3.4. S-модель фильтра (3.10).

Программа расчета частотных характеристик фильтра

R=12e3; Rl=22e3;

L=le3; Ll=le3; C=0.172e-6; Cl=0.078e-6; W=subs(Wy);

Wj =subs(W,'p',1 (j*omega)'); omegaslogspace(0,3,500); Wi=subs(Wj,'omega1,omega); A=abs(Wi);

La=20*logl0(A);

phi=angle(Wi);

subplot(211),semilogx(omega,La),grid

set(gca,'FontSize',12)

ylabel('L(\omega),	дБ')
xlabelC1\omega,	1 / c ' )
subplot(212),semilogx(omega,phi*57.3),grid
set(gca,'FontSize',12)
ylabeK'XphKXomega) , град') ,xlabel('\omega, 1 / c O

В первых трех строчках вводятся значения параметров фильтра. Оператор W=subs(Wy) служит для подстановки численных значений в выражение для передаточной функции. Затем подстановкой р = ju получается частотная передаточная функция Wj. Оператором logspace вводится вектор из 500 значений частоты элементы которого образуют геометрическую прогрессию и лежат в пределах от 101 до 103. Следующий оператор подставляет значения ш в выражение для

частотной передаточной функции, в результате	чего находит-
ся массив значений частотной характеристики	W(ju).	Най-
дя модуль и аргумент этой характеристики функциями		abs

64

и angle, получим амплитудно- и фазочастотную характеристики (АЧХ и ФЧХ). Перейдя к логарифмическому масштабу для АЧХ, получим значение JI АХ (в дБ). Вывод графика выполняем в полулогарифмическом масштабе по х-оси оператором semilogx. Остальные операторы аналогичны операторам, рассмотренным в п. 2.2, и служат для повышения наглядности графика. Результат работы программы представлен на рис. 3.3. Как видно из графика, фильтр обеспечивает подавление сигнала в диапазоне от 70 до 100 [рад/с] не ниже 10 [дБ].

Получим теперь реакцию фильтра на входное воздействие. Пусть входной процесс u(t) представляет собой сумму двух гармоник с одинаковыми (единичными) амплитудами и разными частотами: u(t) = sinu>\t + sinu>2f. Первое слагаемое будем трактовать как "полезный" сигнал, второе - как помеху измерений, вызванную, например, упругостью конструкции объекта. Положим u>i = 10 с- 1 , и>7 = 80 с"1. Составим SIMLINK - программу моделирования. Ее вид представлен на рис. 3.4. В программу входит два блока Sine Wave генераторов полезного сигнала и помехи, блок сумматора Sum и блок решения линейного ОЛУ, заданного передаточной функцией Transfer Fen. Параметры этого блока - Numerator (чи-

слитель) и Denominator (знаменатель) - в программе		заданы
в общем виде: массивами num и den соответственно.		Соглас-
но общим правилам системы MATLAB, многочлены предста-
вляются вектор-строками своих коэффициентов,	расположен-
ными в порядке убывания степени аргумента. В	рассматри-

ваемом примере это должны быть массивы коэффициентов передаточной функции. Как и выше, используем систему SIMULINK для описания модели и собственно процедуры интегрирования, а подготовку исходных данных и отображение результатов оформим в виде MATLAB-программы, представленной ниже. 1

Программа расчета переходных характеристик фильтра

num=conv( [L1*C1 ,R1*C1,1] , [L*C,R*C,1])

den=[Ll*Cl*L*C, R1*C1*L*C+L1*C1*R*C,...

1 Данная программа является продолжением программы расчета частотных характеристик и использует заданные в ней номиналы электротехнических элементов.

з Б. Р. Андриевский и др.

65

R1*C1*R*C+C1*L+L1*C1+L*C,...

C1*R+R*C+R1*C1, 1]; sim('s2 _ l')

f igure plot(t,y),grid

set(gca,'FontSize',12) x l a b e l O t , с ' )

l e g e n d ( ' y ~ * ' , ' y _ c ' , ' y . c + y . n ' , 0 )

text(t(125) ,y(125,l),	' < - - y ~ * ' )
text (t (200) ,y(200,2) ,9 <		y . c ' )
text(t(150) ,y(150,3),	' < — y . c + y . n 1 )

В первых строках программы формируются массивы num и den. Используется оператор свертки векторов conv. В применении к многочленам это соответствует получению вектора коэффициентов их произведения. Заметим, что в программе выражения для многочленов получены несложной корректировкой программы расчета частотных характеристик. Оператор sim запускает на выполнение 5-модель s2_l, после чего оператором figure создается новое графическое окно, в которое выводятся результаты моделирования. Поскольку в монохромном изображении не видны цвета линий графика, для него могут быть использованы операторы text. В результате появляются указатели около соответствующих кривых графика (рис. 3.5). Символом у* обозначен выход фильтра, ус - полезный сигнал, а уп - сигнал помехи (на рисунке пока-

зан входной сигнал ус + уп).		В 5-модели, в меню			Simulation,
опция Parameters,	окно Workspace		I/O, поле	Save	to workspace
установлены параметры:		Time	- t, Output	- у.	Под этими
именами значения	результатов		моделирования		сохраняют-

ся в рабочей области (Workspace) и используются MATLAB-

программой.	На вкладке Solver в окне Solver options		пара-
метр Max step size (максимальная		величина шага) установлен
равным 0.005.	Это сделано не для повышения точности вы-
числений (которая и так высока),		а для получения плавных
кривых при выводе результатов.		Входной и выходной	про-

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

66

2

1

О

- 1
-2.О				t, с
	0.2	0.4	0.6	0.8

Рис. 3.5. Переходные характеристики фильтра (3.10). Выделение сигнала на фоне помехи.

3.2.2.Гармонический анализ процессов

При исследовании колебательных процессов часто применяются их энергетические характеристики, в первую очередь - мощность и энергия [25, 74].

Мгновенная мощность p(t) сигнала y(t) определяется как квадрат его мгновенного значения: p(t) = у(t)2. Энергия Р сигнала на интервале [Ji,^] находится как интеграл от мгновен-

	д <2				р	1	*2
ной мощности Р = J y(t)2dt.			Отношение ^ _ ^			= ^ _ ^ f		y(t)2dt
		h		2	1	2	1 «1
выражает среднюю (на интервале [fi,^])					мощность		сигнала.
Обозначим ее через y(t)2.			Получить представление об					этих
характеристиках		процесса	можно на основе			преобразования
Фурье [25]. Рассмотрим этот метод более					подробно.
Для периодических процессов y(t) с периодом Т							можно за-
писать ряд	Фурье в виде
		оо
y(t)	= ty	+ ^(akcosk^t		+ bk3\nk%l!-ty				(3.11)
где коэффициенты		разложения		находятся из формул
	т			т

67

	рТ
Ьк = т ]	2/(Osin	(k^rt^dt	{к = 1,2,...).
Совокупность величин s0 = |flo|/2,			=	yja\ + b\	(к = 1,2,...)
называется амплитудным		частотным		спектром	периодиче-
ской функции y(t).	Значения s* представляют собой ампли-
				еу
туды гармоник с частотой		= Ш	(Q = уг) в разложении
процесса в ряд (3.11). Они зависят			от	номера	гармоники к

и обычно графически представляются в виде отрезков высотой Sjt, проведенных в точках и>к оси частот (поэтому спектр периодической функции называют линейчатым, или дискретным). Он несет в себе информацию о частотных свойствах сигнала: если сигнал имеет выраженные колебания на неко-

торых частотах, то его спектр				на этих частотах содержит
пики.
Обобщением ряда Фурье на непериодические процессы
является интеграл		Фурье, при котором используется						предста-
вление		гоо
У(0 =	i /	(V(u;)sinu;f		+	U(u)	cos u>t) dt,		(3.12)
где	Jo
	оо					гоо
/	•00y{t)cos(ut)dt,			V(u)=		J-oo	y(t)sin(ut)dt.
Аналогично вводится			частотный		спектр		процесса	y(t) как
S(w) = y/UW	+ V(u)2.		Это функция от непрерывного аргу-
мента и.

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

(ДПФ) Лля этого исследуемый процесс y(t) (длительностью Т) заменяется выборочной дискретной функцией (т.е. после-

довательностью) у[к] =	где tk = кТ0 (к = 1,2,... , N), N -
заданное число точек, Т0 =	j j - шаг дискретности	(кванто-
вания). Далее вычисляется	функция
N
У(*) = £ > [ п ] е х	l<k<N	(3.13)
n=l

68

("изображение по Фурье"), имеющая комплексные значения. В
пакете MATLAB дискретное преобразование Фурье выполня-
ется	процедурой fft. Для ускорения вычислений рекоменду-
ется	брать число точек N = 2", где и - некоторое натуральное

число. В этом случае программой реализуется	так называ-
емое "быстрое преобразование Фурье" (БПФ). 1	Обратный

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

N
2/[п] = - ^ У ( * ) е х р ( ; 2 7 г ( * - 1 ) ^ ) ,	(3-14)
П= 1

Для вычисления спектральной плотности с помощью процедуры fft исходная реализация процесса разбивается на N точек, соответствующих равноотстоящим моментам времени с интервалом Т0. Если исходная последовательность получена в результате моделирования с автоматическим выбором шага (как обычно бывает при использовании стандартных процедур численного интегрирования), то ее следует предварительно преобразовать. Для этой цели можно использовать входящие в пакет MATLAB процедуры интерполяции. К ним относятся, например, функции i n t e r p l или interplq . 2

При выборе параметров вычисления спектра (длина реализации Т, число точек N и связанный с ними интервал квантования Т0 = T/N) следует учитывать, что диапазон существенных частот исследуемого процесса не должен выходить

за частоту Найквиста	= Д-. Несоблюдение этого усло-
	1 о

вия может привести к значительным ошибкам при опреде-

лении	характеристик процесса.	Данное требование вытека-
ет из	известной теоремы отсчетов		Котельникова-Шеннона
[25].	Смысл этого требования	можно пояснить следующим

образом.	Рассмотрим	два гармонических сигнала с частота-
ми и и и	+ u>N: y(t) =	sin(u;J), yx(t) = sin ((u> + и ^ ) 0 . Образу-

ем из них дискретные последовательности с интервалом Т0 :

!/[*] = У{кТ0) = sin(u>fc7o) и			yx[k]	= yi(kT0)	= sin ((си + cuN)kTQ)
(к = 0,1,2,...).		Учитывая,	что	= тг, получим		\у\[к]\ =
1	По-английски - "Fast Fourier		Transform" у сокращенно FFT.
2	П р о ц е д у р а	i n t e r p l q выполняется		быстрее	процедуры	i n t e r p l при

неравномерно заданных значениях аргумента, что характерно для рассматриваемой задачи.

69

= |sin(cjfcTo + 7гА:)| = |y[fc]| (для всех к = 0,1,2,...). Таким образом, полученные в результате дискретизации процессы у[к] и 2/i [к] оказались (с точностью до знака, или начальной фазы) неразличимыми, в то время как исходные непрерывные процессы 2/(f) и j/i(*)i конечно, различны.

Поскольку рассматриваемый процесс у(к) в общем случае не является периодическим с частотой П = 2ж/Т, вычисление его спектра с помощью рассматриваемой процедуры являет-

ся приближенным.		Как видно из		формулы	(3.13), соседние
точки по	частоте отстоят на		величину 6и		=	Учиты-
				о _		^
вая, что	N = Т/То,	получим	6ш =	njr. Таким		образом, дли-

тельность исследуемой реализации должна быть достаточно большой для получения спектра с заданным шагом дискретности по частоте (Т 1/<5и>.)

Задача оценивания спектра по выборке имеет еще одну важную особенность: оценка спектральной плотности, получаемая по ЛПФ, как правило, является смещенной. Это значит, что даже неограниченное уменьшение шага дискретности не приводит к неограниченному убыванию погрешности оценивания. На самом деле средним значением оценки является не точное значение спектральной плотности на заданной частоте, а некоторое сглаженное значение, получаемое усреднением функции - сверткой с некоторым ядром ("окном") [37, 119]. Удачный выбор функции-"окна" позволяет уменьшить погрешность оценивания при заданном конечном

шаге дискретности.
П р и м е р	3.2.1. Выполним гармонический			анализ процес-
сов на входе и выходе рассмотренного в 3.2.1				частотно-изби-
рательного фильтра.		Исследуемые процессы получим с по-
мощью программы, приведенной на с.			65.	Для	обработки
процессов воспользуемся следующей MATLAB-программой.
Программа	расчета	спектров входного	и выходного		сигналов
	фильтра (фрагмент)		1

T=max(t) ;

Nt=2~12

T0=T/Nt;

1 Э т а п р о г р а м м а является продолжением программы, приведенной на с. 65.

70