Загребаев Методы обработки статистической информации в задачах контроля 2008

ri+1 = mod(k ri + b, M ) ,	(4.1.3)
где M – модуль ( M > 0 ); k – множитель (0 < k < M ) ;	b – прира-

щение (0 < b < M ) ; r0 – начальное значение (0 < r0 < M ) .

Рис. 4.2. Схема метода серединных произведений

Последовательность случайных чисел, полученных с помощью данной формулы, называется линейной конгруэнтной последовательностью. Многие авторы называют линейную конгруэнтную последовательность при b = 0 мультипликативным конгруэнтным методом, а при b ≠ 0 смешанным конгруэнтным методом.

Для качественного генератора требуется подобрать подходящие коэффициенты. Необходимо, чтобы число М было довольно большим, так как период не может иметь больше, чем М элементов.

Встроенный генератор случайных чисел. На практике для моделирования случайных величин часто пользуются встроенными в среду разработки генераторами случайных чисел. Например, встроенный в Borland C++ Builder 5.0. датчик случайных чисел генерирует случайные числа с равномерным законом распределения. Ниже приведены результаты работы этого датчика при формирова-

291

нии случайной последовательности из 10000 чисел. На рис. 4.3 приведена гистограмма распределения, а в табл. 4.1 – ряд распределения.

в интервал	0,10
	0,08
попадания	0,06
	0,04
Вероятность
	0,02
	0,00
	0,0	0,2	0,4	0,6	0,8	1,0
			Значения случайной величины

Рис. 4.3. Гистограмма распределения случайных чисел, полученных с помощью датчика случайных чисел

Таблица 4.1

Сравнительная таблица попаданий случайной величины в заданный интервал и значения mi

Интервал Ii		0 – 0,1					0,1 – 0,2	0,2 – 0,3		0,3 – 0,4	0,4 – 0,5
mi		1018					1018	1023		993	974
pin		1000					1000	1000		1000	1000
										Окончание табл. 4.1
Интервал Ii		0,5 – 0,6					0,6 – 0,7	0,7 – 0,8		0,8 – 0,9	0,9 – 1
mi		977					1005	1029		1021	942
pin		1000					1000	1000		1000	1000
В таблице		p* =	mi		,		m – число значений,			попавших в i -й ин-
		i		n			i

тервал, n – общее число испытаний, pi – теоретическая вероятность попадания в i -й интервал.

292

Определим значение меры расхождения:

χ 2	k	( p* − p )2		k	(m − np )2
	= n∑	i	i	= ∑	i	i	.	(4.1.1)
			p
	i=1			i=1		np
			i			i

Исходя из данных табл. 4.1, получим χ2 = 7,102 . Число степе-

ней свободы 7, тогда из таблицы значений χ2 получаем при уровне

значимости 5 % интервал (1,56 – 16,62). Рассчитанное значение критерия попадает в указанный интервал, таким образом, гипотеза о виде закона распределения не противоречит опытным данным.

4.2. Моделирование случайных величин с нормальным законом распределения

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

Алгоритм вычисления можно представить следующим образом. Берется датчик случайной величины, которая имеет равномерную плотность распределения вероятности в интервале от 0 до 1. При этом математическое ожидание случайной величины равно 0.5, а ее дисперсия равна 1/12. Если просуммировать порядка 8 – 10 независимых реализаций СВ, то результатом суммирования будет случайная величина с плотностью распределения вероятности, очень близкой к нормальному закону.

N
x = ∑ yi ,	(4.2.1)

i=1

293

где yi	– случайные величины с равномерной плотностью распре-
деления вероятности в интервале от 0 до 1. Математическое ожи-
дание и дисперсия величины x соответственно равны:
			M[Y ] =	N ;	D[Y ] =	N .
				2		12
Ниже приведены результаты, полученные при использовании
формулы (4.2.1)		при	N =10 . Гистограмма распределения,						полу-
чившаяся при моделировании последовательности из 10000 слу-
чайных чисел, показана на рис. 4.4.
	0,20
интервал	0,18
	0,16
	0,14
в
попадания	0,12
	0,10
	0,08
Вероятность	0,06
	0,04
	0,02
	0,00
	-4	-3	-2	-1	0	1	2	3	4
			Значения случайной величины

Рис.4.4. Гистограмма распределения при генерации случайных чисел по нормальному закону

Рассчитанное значение критерия χ2 =14,996 . Число интервалов

равно 16, при уровне значимости 5 % получим интервал (4,76 – 25,5), т.е. гипотеза не противоречит опытным данным.

294

4.3.Моделирование случайных величин

сзаданным видом корреляционной функции (метод формирующего фильтра)

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

a		d 2 x(t)	+ a	dx(t)	+ a x(t) = b		dξ(t)	+ b	ξ(t),	(4.3.1)
	2 dt2
			1 dt		0	1 dt		0

где x(t) – случайная величина с заданной корреляционной функцией; ξ(t) – белый шум.

Для такого вида уравнений формирующего фильтра аналитически известен вид корреляционной функции. Если корни характеристического уравнения действительные, то

λ

b2

a0

+ b a2a

exp(λ τ) −

K (τ) =

2

a2

+ b2

2

0

1

a

1 0 1

1

2

λ

b2

a0

exp(λ

t

−

a2

+ b2

+ b a2a

τ)

;

2

1 2

0

1

1 0 1

2a

a a2

(λ

− λ

)

a2

2

1

0 1 2

τ = m

t.

(4.3.2)

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

G G	N	(4.3.3)
J (a,b) = ∑(K (t j ) − K (t j ))dt,

j=1

где K (t) – корреляционная функция формирующего фильтра, оп-

~

ределяемая исходя из формулы (4.3.2), K (t) – заданная корреляци-

онная функция.

Эмпирическая корреляционная функция рассчитывается по формуле:

		1	N −m−1
	K (τ) = K (m t) =		∑	x(i)x(i + m) .	(4.3.4)
		N − m
			i=0

295

Например, для значений коэффициентов: a2 =1; a1 = 3; a0 = 2; b1 = 0; b0 =1 результаты моделирования представлены на рис. 4.5,

на котором приведены графики теоретической корреляционной функции, рассчитанной по формуле (4.3.2), и эмпирической корреляционной функции, рассчитанной по формуле (4.3.4), исходя из значений случайной величины, полученных при решении уравнения (4.3.1) методом Эйлера.

Ошибка, рассчитанная по формуле

N

(K(t

) − K

(t

2

∑

j

j

))

δ =

j=1

100 % ,

(4.3.5)

N

K(t

)

2

j

∑

j=1

составляет ≈ 4 %.

Рис. 4.5. Эмпирическая и теоретическая корреляционные функции

Экспоненциальная корреляционная функция. Если корреля-

ционная функция задана в виде:

296

K (τ) = D exp(−ατ); τ = m t ,

(4.3.6)

то случайную величину можно рассчитывать по следующей формуле:

x(tn ) = x(n t) = x(n) = k1e(n) + k2 x(n −1) ;				(4.3.7)
k =	D(1− k2 ); k	2	= exp(−αΔt),	(4.3.8)
1	2

где e(n) – дискретный белый гауссов шум с нулевым математиче-

ским ожиданием и единичной дисперсией.

Например, при D =1 , α = 0,05 , =1 графики эмпирической и теоретической корреляционных функций представлены на рис. 4.6.

Рис. 4.6. Корреляционные функции

Погрешность, рассчитанная исходя из формулы (4.3.5) составляет ≈ 3 %.

Если корреляционная функция имеет вид:

K (τ) = D exp(−a2τ2 ); τ = m t ,

(4.3.9)

то последовательность этапов моделирования следующая.

297

1.Необходимо получить реализацию дискретного белого шума длительностью N (где N достаточно большое – порядка 1000 и более отсчетов) с нулевым математическим ожиданием и единичной дисперсией. Для получения данной реализации необходимо N раз обратиться к датчику, выдающему независимые случайные числа, распределенные по нормальному закону. Эту реализацию в дальнейшем будем обозначать e(n).

2.Далее выполняются следующие преобразования:

P		2Dα
x(n) = ∑ C(k)e(n − k);	C(k) =			exp(−2α2k2 ).	(4.3.10)
		4
k=−P			π

Здесь неопределенным остается предел суммирования P , и для его определения может служить рекомендация: P – целая часть от деления 2 на α ( α <1 , иначе моделирование не имело бы смысла). После этого в получившейся реализации необходимо отбросить первые и последние P отсчетов и оставить только N − 2P отсчетов. Дело в том, что стационарным фрагментом моделируемого случайного процесса (с постоянной дисперсией) является именно центральная часть. Таким образом, длительность реализации равна

N = N − 2P .

На рис. 4.7 показан пример расчета при D =1 , α = 0,05 , =1.

Рис. 4.7. Корреляционная функция

298

Погрешность, рассчитанная по формуле (4.3.5), составляет

1,603 %.

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

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

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

Вкачестве основы для создания математической модели с флюктуирующими параметрами был принят подход Шефа и Альбрехта, использующий для получения статистических характеристик плотности потока нейтронов стохастические уравнения. Ниже приводятся стохастические уравнения пространственной динамики, включая обратные связи по температуре топлива, теплоносителя, замедлителя, паросодержания (для кипящего реактора), ксенонового отравления и системы регулирования [2 – 4]. Источником случайных возмущений при этом считается флюктуация в размножающих свойствах среды.

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

299

щениями являлись флюктуации температуры топлива и теплоносителя и плотности теплоносителя.

4.4.1. Математическая модель и статистические исследования параметров ячейки реактора

Для исследования статистических характеристик ячейки реактора использовалась классическая модель нейтронно-физического расчета ячейки тепловых реакторов, основанная на следующих предположениях. В резонансной области энергии – это приближение плоского потока и эффективных резонансных интегралов. В области тепловых нейтронов энергетическая зависимость описывалась распределением Максвелла с рассчитываемой температурой нейтронного «газа». Пространственное распределение тепловых нейтронов рассчитывалось по методу вероятностей первых столкновений.

Схема проведения статистических исследований выглядит следующим образом: задаются характерные для данного типа реактора геометрические параметры ячейки и изотопный состав, номинальные значения плотности теплоносителя, температуры топлива и температуры теплоносителя. Затем с помощью генератора случайных чисел разыгрываются отклонения от номинальных значений по заданному закону распределения, и для этих значений параметров производится расчет нейтронно-физических параметров ячейки. Результаты расчетов подвергаются статистической обработке на предмет получения закона распределения, его числовых характеристик, максимальных и минимальных отклонений параметров. На рис. 4.8 схематично представлена логика исследования.

плотность				k, Σa ,ν f Σ f , M 2	,τ
температура

Рис. 4.8. Схема статистических исследований ячейки

300