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

фильтра, границы полосы пропускания и вид фильтрации: значение ' s t o p ' задает заграждающий фильтр (ЗФ); значения 'high1 , 'low', 'bandpass' задают фильтр верхних частот

(ВЧ), нижних частот (НЧ) и полосовой фильтр (ПФ), соответственно. Чтобы получить аналоговый фильтр, следует указать четвертый параметр ' s \ 1 Заметим, что порядок по- лученных ПФ и ЗФ равен удвоенному значению первого параметра. В рассматриваемом примере рассчитываются матрицы разностных уравнений состояния заграждающего фильтра четвертого порядка с полосой задерживания от 70 до 100 рад/с. Результат присваивается массивам а, Ь, с, d. Следующим оператором создается JIC-объект с именем fb. Для этого используется функция ss тулбокса CONTROL SYSTEMS (см. Приложение и [61]). Выполнив оператор zpk(fb), получим выражение для передаточной функции W(z) фильтра в виде произведения сомножителей:

Wlz) =	0.9(z2 - l-83z +	l)2
v ;	(z2 - 1.78г + 0.91 ){z2 -	1.68* + 0.89)"

Для расчета частотных и переходных характеристик фильтра используем рассмотренные выше функции bode step. Результаты приведены в следующем пункте на рис. 3.14 а.

3.4.3.Континуализация моделей

Рассмотрим теперь задачу континуализации дискретных моделей. В качестве примера построим сначала непрерывный фильтр, используя полученный выше цифровой фильтр Баттерворта в качестве дискретного фильтра-прототипа.2 Для получения непрерывной модели фильтра используем функцию преобразования d2c тулбокса CONTROL SYSTEMS:

cfb=d2c(fb)

1 Третий параметр может принимать и значения 'в* или , z ) (установлено по умолчанию). Этот параметр показывает тогда, какой фильтр - аналоговый или цифровой - должен быть получен. При этом, если второй параметр - скаляр, то рассчитывается НЧ-фильтр, а если вектор из Двух элементов - то ПФ.

2 Заметим, что пример носит иллюстративный характер и данный подход не характерен для проектирования частотно-избирательных фильтров.

91

Получим JIC-объект cfb, который соответствует модели дискретной системы fb. При переходе к этой модели учтено требование совпадения переходных характеристик. Выполнив затем оператор zpk(cfb), получим передаточную функ-

цию		непрерывной системы в виде
	w (	v=	0.9(д2 + 8.36s + 6522)(s2			- 10.15 + 8.4 • 103)
		[ S )	(52	+ 18.75 + 5433) (в2		+ 23.75 + 9054) '
На	рис.	3.14,	а показаны		амплитудные характеристики ци-
фрового фильтра (кривая					1) и его непрерывной модели (кри-
вая	2). Ниже		(рис.	3.14,	в) приведены графики переходных
характеристик.

Как известно, дискретные ЛС обладают рядом свойств, которые не могут быть адекватно отражены непрерывными моделями. Например, переходный процесс у стационарных непрерывных Л С всегда носит асимптотический характер, а у дискретных систем может иметь конечную длительность. Кроме того, дискретные системы первого порядка в отличие от непрерывных могут иметь колебательную переходную характеристику. Рассмотрим на примерах, какая модель полу-

чается в этих			случаях.
Пример		3.4.1.		Создадим дискретную систему операто-
ром	s = t f ( l , [ l		0	0] ДО). Это система с передаточной функ-
цией	W(z)	=	\/z2	и интервалом дискретности Г0 = 0.005 с.

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

работы	(т.е.	при t	= 0.01 с). Оператор c=d2c(s)	приводит к
сообщению об ошибке ("D2C with ZOH: cannot handle					systems
with poles at		: "Функция d2c с экстраполятором			нулево-
го порядка		не может использоваться для систем, имеющих
нулевые	полюсы").		Можно использовать другой метод пе-
рехода,	например метод Тастина [6, 24, 19, 74], и написать
c=d2c(s, ' t u s t i n ' ) .			Получим ответ: W„(«) =	~ f 5 +	f . Пе-
			5 +	4 5 +	4

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

Пример 3.4.2. Возьмем теперь дискретную систему, заданную передаточной функцией W(z) = z q g. Несложно

92

Рис. 3.14. Континуализация моделей.

убедиться, что такая система имеет колебательную переходную характеристику. Получим непрерывную модель этой системы. Выполним операторы:

s l = t f ( 1 . 5 , [1 0.5],Т0) csl=d2c(sl)

В результате появляется сообщение "System order was increased to handle real negative poles": "Порядок системы был увеличен, чтобы обработать отрицательные вещественные полюсы". Непрерывная система csl имеет второй порядок и описывается передаточной функцией WJs) = s*+ 2ns + 4.14^ • ,10 • Переходные характеристики исходной дискретной системы и ее непрерывной модели в моменты tk = кТ0)к = 0,1,2,... совпадают (рис. 3.14, г).

З а м е ч а н и е . Стоит обратить внимание на то, что в некоторых версиях тулбокса CONTROL SYSTEMS функция d2c для систем, заданных в НПК-форме (в виде совокупности нулей, полюсов и коэффициента передачи), работает не-

верно.	Например, выполнив оператор s2=zpk(sl), получим
ту же дискретную систему, что и исходная, W(z)		= z	q ^,
однако для нее оператор cs2=d2c(s2) дает ошибочный			ответ
WH(s)	= s ^ q ^. В версии 4.2 тулбокса эта ошибка	устранена.

93

3.5.Детерминированные и стохастические модели

Модели систем, о которых мы говорили до сих пор, были детерминированными (определенными), т.е. задание входного воздействия определяло выход системы однозначно. Однако на практике так бывает редко: описанию реальных систем обычно присуща неопределенность. Например, для статической модели неопределенность можно учесть, записывая вместо (3.1) соотношение

y(t) = F(u(t)) + <p(t),

(3.44)

где ip(t) - погрешность, приведенная к выходу системы. Причины неопределенности разнообразны:

-погрешности и помехи измерений входов и выходов системы (естественные погрешности);

-неточность самой модели системы, учитываемая путем искусственного введения в модель погрешности;

-неполнота информации о параметрах системы и т.д. Среди различных способов уточнения и формализации не-

определенности наибольшее распространение получил стохастический (вероятностный) подход, при котором неопределенные величины считаются случайными. Развитый понятийный и вычислительный аппараты теории вероятностей и математической статистики позволяют дать конкретные рекомендации по выбору структуры системы и оценке ее параметров. Классификация стохастических моделей систем и методов их исследования представлена в табл. 3.3.

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

ном	больших	чисел	и центральной	предельной	теоремой. За-
кон	больших	чисел	гласит, что	если	~ случайные

величины с математическим ожиданием (среднее значение)

= а и дисперсией М(£,	-	а)2	= <72, то
+	-	+	0	(3.45)
при достаточно больших N. Это говорит о принципиальной
возможности сколь угодно точной оценки				по измерениям.

94

Т а б л и ц а З . З . Стохастические модели систем

	Статические		Дискретные по Т		Динамические
					Непрерывные по Т
	Дискретные	Непрерывные	Дискретные	Непрерывные	Дискретные	Непрерывные
Математи-	по С/, У	по l/,Y	по U,Y	по UyY	по U, Y	по U,Y
	Схема незави-	Регрессион-	Марковские	Стохастичес-	Системы массо-	Стохастичес-
ческий ап-	симых испыта-	ные модели	цепи; стохас-	кие разност-	вого обслужи-	кие диффе-
парат опи-	ний [101, 102]	[101, 102]	тические ав-	ные уравне-	вания; полумар-	ренциальные
сания			томаты [88]	ния [75]	ковские про-	уравнения
					цессы [15, 94]	[53]
Методы	Статистические	Регрессион-	Оценка пере-	Статистическое	Теория массо-	Теория устой-
оценки па-	оценки вероят-	ный анализ	ходных веро-	оценивание со-	вого обслу-	чивости [75]
раметров	ности; дисперси-	[1, 86, 78, 110, 119]	ятностей; ста-	стояний и пара-	живания; ими-
и анализа	онный анализ		тистическое	метров; анализ	тационное мо-
	[46, 110]		моделирова-	стох астической	делирование
			ние [15, 27]	устойчивости [75, 86] [103]
Методы	Стохастическое	Планирование	Динамическое	Динамическое	Перебор; мето-	Оптимальное
синтеза	программирова-	эксперимента;	программиро-	программиро-	ды оптималь-	и адаптивное
	ние [80]	стохастическое	вание [88]	вание [75, 95]	ного управле-	управление
		программиро-			ния [103]	[75, 95]
		вание [78, 110]
Области	Задачи выбора	Обработка	Компьютеры	Импульсные	Системы обслужи-	САУ; механи-
примене-	из конечного	результатов		и цифровые	вания (вычисл.	ческие, электронные,
ния	числа вариантов	измерений		САУ	системы, про-	тепловые и другие
	(испытания,	и испытаний			изводственные	процессы
	управление)				системы и т.д.)

Центральная предельная теорема, уточняя (3.45), утверждает, что

	^ K i	+ ... + Ы - в	» ^ ,		(3.46)
где	стандартная	= О, М£2 =	1) нормально		распреде-
ленная случайная величина. (Если			величина	а2 неизвестна,
то следует заменить в (3.46) о на оценку s= yjjjzГ					]Cili(&"~02i
где £ =	jj'YliLi^i- При этом величина f будет			распределена

уже не нормально, а по закону Стьюдента с числом степней свбоды N — 1, который, к счастью, при 7V > 20 практически неотличим от нормального).

	Функция распределения		нормальной случайной	величины
£	хорошо	известна и детально затабулирована.		Имеются
в	книгах	по статистике и	таблицы функций распределения

других часто встречающихся случайных величин, например стьюдентовой. Однако нет необходимости добывать статистические таблицы, если под рукой есть компьютер с установленной системой MATLAB. В составе тулбокса STATISTICS есть модули, вычисляющие функции распределения и другие характеристики более двадцати распространенных типов случайных величин. Узнать перечень всех функций тулбокса можно с помощью команды help tooIbox\stats или встроенной системы помощи через меню Help, окно Help Window. Проиллюстрируем применение MATLAB на простой статистической задаче.

Пусть требуется найти, при каком числе измерений N по-

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

меньше чем 0.02, с вероятностью не меньше чем 0.99. Пусть среднеквадратическая погрешность каждого измерения из-

вестна: а = 0.5, а в качестве оценки берется среднее		арифме-
тическое величин	• • •
Допустим, что N > 30, так что можно пользоваться цен-
тральной предельной	теоремой. Из (3.46) следует, что иско-
мое N должно удовлетворять соотношению Р		< 0.02^ >
>0.99, или эквивалентному соотношению Р j|£| <		j >
> 0.99, где £ - нормальная случайная величина со		средним

96

О и дисперсией 1. Вероятность, стоящую в левой части последнего соотношения, легко определить с помощью функции norminv, вычисляющей по величине р значение х) удовлетворяющее равенству Р {£ < х) = р. В силу симметрии нормального распределения Р{|f| < х} = 2р — 1. Поэтому значение р, удовлетворяющее заданному условию 2р— 1 = 0.99, равно 0.995. Выполнение оператора norminv(0.995) дает результат

ans = 2.578, откуда y/N	>	или N> ( 0 5 ц (ji>578)2 «
« 4154. Поскольку 4154 >	30, применение центральной пре-

дельной теоремы законно (в противном случае вместо функции norminv следовало бы воспользоваться функцией tinv, обращающей функцию распределения Стьюдента). В заключение приведем полезное равенство Р{|£| < 1.96} = 0.05, позволяющее в прикидочных расчетах обходиться и без таблиц, и без компьютера. Интересно, что для решения задачи нам не понадобилась точная формулировка центральной предельной теоремы. Разумеется, формулировкам (3.45), (3.46) можно придать более строгий вид и это легко достижимо с помощью понятий вероятностной сходимости. Однако при попытке проверить условия этих строгих утверждений могут возникнуть трудности. В частности, в законе больших чисел и центральной предельной теореме требуется независимость отдельных измерений (реализаций) случайной вели-

чины и конечность ее дисперсии. Если		эти условия		наруша-
ются, то могут нарушаться и выводы.			Например,	если все
измерения совпадают: ^ = ... =	то,		хотя все остальные

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

(с плотностью распределения р(х) = 1 /тг(1 + я2)), не обладающему конечным математическим ожиданием и дисперсией. А ведь такой закон встречается в жизни! Например, по Коши распределена интегральная освещенность точек прямолинейного берега равномерно вращающимся прожектором, находящимся в море (на корабле) и включающимся в случайные моменты времени.

Но еще большие трудности вызывает проверка обоснованности самого употребления термина "случайный". Что такое случайная величина, случайное событие и т.д.? Часто говорят, что событие А случайно, если в результате экспе-

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

римента оно может наступить (с вероятностью р) или не наступить (с вероятностью 1— р). Все, однако, не так просто. Сама вероятность события может быть связана с результатами экспериментов лишь через частоту наступления события в некотором ряде (серии) экспериментов: vN = NA/N, где NA - число экспериментов, в которых событие наступило; N - общее число экспериментов. Если числа и^ при достаточно большом N приближаются к некоторому постоянному числу

РА-
^ « Р л ,	(3.47)

то событие А можно назвать случайным, а число р - его вероятностью. При этом частоты, наблюдавшиеся в различных сериях экспериментов, должны быть близки между собой

(это свойство называется статистической	устойчивостью,
или однородностью). Сказанное относится	и к понятию слу-

чайной величины, поскольку величина f является случайной, если случайными являются события {а < £ < 6} для любых чисел а, 6. Частоты наступления таких событий в длинных сериях экспериментов должны группироваться около некоторых постоянных значений.

Итак, для применимости стохастического подхода должны выполняться следующие требования:

1)массовость проводимых экспериментов, т.е. достаточно большое их число;

2)повторяемость условий экспериментов, оправдывающая сравнение результатов различных экспериментов;

3)статистическая устойчивость.

Стохастический подход заведомо нельзя применять к единичным экспериментам: бессмысленны выражения типа "вероятность того, что завтра будет дождь", "с вероятностью 0.8 « З е н и т » выиграет кубок" и т.п. Но даже если массовость и повторяемость экспериментов имеются, статистической устойчивости может и не быть, а проверить это - непростое дело. Известные оценки допустимого отклонения частоты от вероятности основаны на центральной предельной теореме1 или неравенстве Чебышева и требуют дополнительных гипотез о независимости или слабой зависимости измерений.

	1 Легко	заметить, что	(3.47) есть частный	случай (3.45), когда берет-
ся	= 1,	если событие	А наступило в i-м	эксперименте,	и £, = 0 - в
противном случае. При этом Л/& = рд,				- рА)2 = рл(1 -	р^) .

98

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

Более подробно методология и практические рецепты применения теории вероятностей изложены в поучительной книге В.Н. Тутубалина [101], представление о которой дают приводимые ниже цитаты.

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

"По современным представлениям область применения те- оретико-вероятностных методов ограничена явлениями, которым присуща статистическая устойчивость. Однако проверка статистической устойчивости трудна и всегда неполна; к тому же часто она дает отрицательный вывод. В результате в целых областях знания, например в геологии, нормой стал такой подход, при котором статистическая устойчивость вообще не проверяется, что неизбежно приводит к серьезным ошибкам. К тому же пропаганда кибернетики,1 предпринятая нашими ведущими учеными, дала (в некоторых случаях) несколько неожиданный результат: теперь считается, что только машина (а не человек) способна получать объективные научные результаты.

В таких обстоятельствах долг каждого преподавателя - вновь и вновь пропагандировать ту старую истину, которую еще Петр I пытался (безуспешно) внушить русским купцам, что торговать надо честно, без обмана, так как в конечном счете это для самих же себя выгоднее."

Как же построить модель системы, если неопределенность

взадаче есть, но стохастический подход неприменим? Ниже,

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

Но сначала рассмотрим практический пример.

1 А сегодня можно добавить "... и информатики" (прим. авторов).

99

3.6» Пример статистической обработки данных программами тулбокса STATISTICS

Рассмотрим в качестве примера курс акций компании "Coca- ColaR" за 1995 г. График исходного процесса представлен на рис. 3.15, а. 1 Пусть данные о курсе акций содержатся в массиве у, находящемся в файле cocdat.mat. Выполним следующие операторы:

load cocdat plot(tdat,у) xlabel('Дата') ylabel('Цена, дол . ')

title('Kypc акций компании Coca-Cola~R') timeser(datenum('01-Jan-1995')); dateaxis('x',2)

Процедура timeser входит в тулбокс FINANCIAL. Она преобразует разметку оси X графика в соответствующие календарные даты. При обращении к ней указывается начальная дата, которая с помощью процедуры datenum предварительно преобразуется из строкового формата в числовой. Функция dateaxis выполняет обратное преобразование, в результате которого разметка оси X приобретает требуемый символьный вид (формат записи дат задается вторым параметром процедуры). Выполним обработку процесса. Сначала вычислим и покажем на этом же графике изменение скользящего среднего рассматриваемого процесса. Чтобы иметь более полное представление о процессе, получим упреждающее и запаздывающее скользящие средние. С этой целью используем функцию movavg тулбокса FINANCIAL. Ее параметром является число дней, за которое производится усреднение. Вызов этой функции имеет в общем случае вид

[short,long] = movavg(asset,lead,lag,alpha)

Здесь short и long - соответственно упреждающее и запаздывающее скользящие средние; asset - подлежащие обработке исходные данные; lead и lag - количества точек, по которым вычисляются упреждающее и запаздывающее скользящие средние; alpha - управляющий параметр, который

1 Л л я простоты отображения исходные данные в 245 точках интерполированы на 365 точек, т.е. каждая точка соответствует суткам.

100