Имитационная модель оценки и прогнозирования эффективности поиска подводной лодки

Курсовая работа

Имитационная модель оценки и прогнозирования эффективности поиска подводной лодки

1. Постановка задачи на решение

2. Q — схема процесса функционирования системы

3. Описание переменных и констант

4. Математическая модель

5. Накопление результатов моделирования

6. Обработка результатов моделирования

7. Блок-схема обобщенного алгоритма модели

8. Блок-схемы процедур алгоритма

9. Листинг программы

10. Результаты (количество обнаруженных пл из 100 прогонов)

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

12. Дисперсионный анализ результатов моделирования

13. Регрессионный анализ

14. Покоординатный спуск

15. Наискорейший спуск

Заключение

программа имитация модель поиск подводная лодка

Перед началом поиска подводной лодки корабельно-поисковой ударной группой (КПУГ) по вызову важно спрогнозировать эффективность поиска. Эффективность поиска зависит от расстояния между КПУГ и точкой потери контакта (ТПК).

Кроме того, эффективность поиска пл зависит от:

— интенсивности обнаружения кораблей КПУГ подводной лодкой;

— интенсивности обнаружения пл средствами КПУГ без уклонения пл;

интенсивности обнаружения пл средствами КПУГ при уклонении пл;

— предполагаемой скорости пл;

— скорости выхода кораблей КПУГ в позицию начала поиска.

В работе приведена достаточно упрощенная имитационная модель процессов:

— выхода кораблей КПУГ в позицию начала поиска пл;

— движения пл до момента выхода кораблей КПУГ в позицию начала поиска пл;

— обнаружения подводной лодки кораблями КПУГ после начала поиска без уклонения пл;

— обнаружения подводной лодкой кораблей КПУГ и ее уклонения;

— обнаружения подводной лодки кораблями КПУГ после начала поиска при уклонении пл.

Время поиска пл ограничено запасом топлива на кораблях, количеством и времени жизни радиогидроакустических буев (РГАБ), уменьшением вероятности обнаружения пл при окончании обследования КПУГ области возможного местоположения (ОВМ) пл. Поэтому в имитационной модели введено время продолжительности поиска.

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

Если пл обнаружит КПУГ раньше, то она будет уклоняться. В этом случае интенсивностью поиска уменьшится. И момент времени обнаружения пл будет вычисляться заново, как значение случайной величины времени обнаружения пл с начала поиска, распределенной по экспоненциальному закону с новой интенсивностью поиска.

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

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

Ограничениями и допущениями модели являются:

— время устаревания данных от ТПК равно «0»;

— предполагаемая скорость пл меньше скорости выхода кораблей КПУГ в позицию начала поиска;

— эти скорости имитационной модели не будут меняться;

— поисковая скорость кораблей КПУГ не учитывается;

— предполагается, что область возможного местоположения подводной лодки ограждена барьерами РГАБ, выставленных корабельными вертолетами КПУГ;

— интенсивности обнаружения кораблей КПУГ подводной лодкой, обнаружения пл средствами КПУГ с уклонением и без уклонения известны.

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

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

— расстояние между КПУГ и ТПК;

— скорость выхода кораблей КПУГ в позицию начала поиска;

— интенсивность обнаружения пл средствами КПУГ без уклонения пл;

— интенсивность обнаружения пл средствами КПУГ при уклонении пл.

Имитационная модель создана с использованием математической схемы (Q-схемы), применимой для исследования систем массового обслуживания. Программа модели написана на языке Borland Pascal 7.0. Произведена отладка и верификация программы, определено необходимое количество повторений (прогонов) основной процедуры программы для обеспечения заданной точности и достоверности результатов моделирования (вероятности обнаружения пл кораблями КПУГ), доказана возможность принятия гипотезы о распределении случайной величины отклика (количества повторений, в которых пл обнаружена) по нормальному закону.

В результате проведения экспериментов с программой имитационной модели построено уравнение регрессии, показывающее зависимость отклика (вероятности обнаружения пл кораблями КПУГ) от факторов: расстояние между КПУГ и ТПК, интенсивностей обнаружения кораблей КПУГ подводной лодкой, обнаружения пл средствами КПУГ без уклонения пл и при уклонении. Проверена адекватность уравнения регрессии и программной модели. С использованием методологии поверхности отклика и уравнения регрессии выявлены значения расстояния между КПУГ и ТПК, интенсивностей обнаружения кораблей КПУГ подводной лодкой, обнаружения пл средствами КПУГ без уклонения пл и при уклонении, при которых вероятность обнаружения пл кораблями КПУГ является максимальной.

Таким образом, предложен способ оценки и прогнозирования эффективности поиска пл при известных интенсивностях поиска и расстоянии между КПУГ и ТПК при помощи имитационного моделирования.

1. Постановка задачи на решение

Корабельная поисковая ударная группа осуществляет поиск подводной лодки противника по вызову.

С получением доклада о координатах точки потери контакта с пл корабли КПУГ начинают движение в позицию начала поиска, которая находится на дистанции 15 км от ТПК. ТПК в момент получения КПУГ доклада находится на дистанции 25±10 км от КПУГ.

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

Если к моменту выхода кораблей КПУГ в позицию начала поиска радиуса ОВМ пл оказывается не больше 15 км, то корабли КПУГ начинают поиск пл.

Если во время поиска пл обнаруживает корабли КПУГ раньше, чем будет обнаружена сама, то она предпринимает маневр уклонения.

Интенсивность обнаружения кораблей КПУГ пл 0,03 1/мин.

Интенсивность обнаружения пл средствами КПУГ без уклонения пл 0,02 1/мин.

Интенсивность обнаружения пл средствами КПУГ при уклонения пл 0,01 1/мин.

Определить вероятность обнаружения пл за 2 часа с момента получения доклада о координатах ТПК, если скорость пл 10 м/с, а скорость выхода кораблей КПУГ в позицию начала поиска 20 м/с.

Считать, что время устаревания координат ТПК на момент доклада равно «0», а закон распределения моментов времени обнаружения с начала поиска — экспоненциальный.

2. Q — схема процесса функционирования системы

Описание Q-схемы

Источники:

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

Источники U₁, U₂ и U₃— предназначены для генерации величин, моментов времен обнаружения на основании интенсивностей обнаружения пл в первом, втором и обнаружения КПУГ в третьем канале соответственно. Эти величины распределены по экспоненциальному закону :

— 0,03 1/мин для источника U₃,

— 0,02 1/мин для источника U₁,

— 0,01 1/мин для источника U₂,

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

Каналы:

Канал К₁ — предназначен для обслуживания заявки, поступающей в систему из источника U_z. Момент времени выхода КПУГ в ОВМ пл.

Канал К₂ — предназначен для обслуживания заявки, поступающей в систему из источника U_z. Момент времени выхода пл из радиуса ОВМ.

Канал К₃ — предназначен для обслуживания заявки, поступающей в него после окончания обслуживания в канале К₂. Момент времени обнаружения пл средствами КПУГ без уклонения.

Канал К₄ — предназначен для обслуживания заявки, поступающей в него после окончания обслуживания в канале К₂. Момент времени обнаружения пл средствами КПУГ с уклонением.

Канал К₅ — предназначен для обслуживания заявки, поступающей в него после окончания обслуживания в канале К₂. Момент времени обнаружения кораблей КПУГ пл.

Канал К₆ — предназначен для обслуживания заявки, поступающей в систему из источника U_z. Момент времени, после которого обнаружения пл средствами КПУГ не производится.

Клапана:

Клапан 1 — предназначен для пропускания заявки из канала К₁. В исходном состоянии клапан открыт. Клапан закрывается в момент времени выхода пл из радиуса ОВМ.

Клапан 2 — предназначен для пропускания заявки из канала К₁, в момент времени выхода пл из радиуса ОВМ. В исходном состоянии клапан закрыт.

Клапан 3 — предназначен для пропускания заявки из канала К₁, в момент времени обнаружения пл средствами КПУГ без уклонения. В исходном состоянии клапан открыт. Клапан закрывается в момент времени обнаружения пл средствами КПУГ с уклонением.

Клапан 4 — предназначен для пропускания заявки из канала К₁, в момент времени обнаружения пл средствами КПУГ с уклонением. В исходном состоянии клапан открыт. Клапан закрывается в момент времени обнаружения пл средствами КПУГ без уклонения.

Клапан 5 — предназначен для пропускания заявки из источника U₁, в момент времени обнаружения пл средствами КПУГ без уклонения.

Клапан 6 — предназначен для пропускания заявки из источника U₂, в момент времени обнаружения пл средствами КПУГ с уклонением.

Клапан 7 — предназначен для пропускания заявки из источника U₃, в момент времени обнаружения кораблей КПУГ пл.

Клапан 8 — предназначен для пропускания заявки из источников U₁, U₂ или U₃, в момент времени, после которого обнаружения пл средствами КПУГ не производится, клапан открывается. В исходном состоянии клапан закрыт.

Клапан 9 — предназначен для пропускания заявки из источников U₁, U₂ или U₃, в момент времени до истечении заданного. В исходном состоянии клапан открыт. Закрывается если время станет больше 2 часов.

3. Описание переменных и констант

Константы:

1. R=15 000 — радиус ОВМ пл в метрах;

2. Ttreb=7200 — время одного прогона в секундах;

3. Vpl=15 — скорость пл в м/с;

4. Vnk=15 — скорость КПУГ в м/с;

5. Y1=0,03/60 — интенсивность обнаружения КПУГ пл в 1/сек;

6. Y2=0,02/60 — интенсивность обнаружения пл средствами КПУГ без уклонения в 1/сек;

7. Y3=0,01/60 — интенсивность обнаружения ПЛ средствами КПУГ с уклонением в 1/сек;

8. D0min=15 000 -минимальное расстояние от КПУГ до ТПК в метрах;

9. D0max=35 000 — максимальное расстояние от КПУГдо ТПК в метрах;

10 .M=100 — количество больших прогонов;

11. N=100 -количество малых прогонов;

12. Tb=1,96 -квантиль;

13. е0=0,05 — заданная относительная точность;

Переменные:

1. D0 — расстояние от КПУГ до ТПК;

2. t — текущее время;

3. Nobn — количество обнаруженных пл;

4. Nnеobn — количество необнаруженных пл;

5. tk1, tk2, tk3, tk4, tk5 — моменты выхода заявки из каналов 1, 2, 3, 4 и 5;

6. Pobnвероятность обнаружения пл;

7. SMO — сумма математического ожидания;

8. SM2 — квадрат суммы математического ожидания;

9. МО — математическое ожидание;

10. D — дисперсия;

11. E — относительная точность;

12. Nneobh — необходимое количество прогонов.

4. Математическая модель

1. Pobn:=Nobn/N — основная формула. Вероятность обнаружения пл;

2. Nobn:=Nobn+1, если (t=tk3) или (t=tk4) — накопление обнаруженных пл;

3. tk3:=t-ln (Random)/Y2, если (t=tk1) и (tk2>tk1) — расчет момента вре-мени обнаружения пл средствами КПУГ без уклонения;

4. tk4:=t-ln (Random)/Y3, если (t=tk5) и (tk5

5. tk5:=t-ln (Random)/Y1, если (t=tk1) и (tk2>tk1) — расчет момента вре-мени обнаружения КПУГ пл;

6. tk1:=t+(D0-R)/Vnk, если t=0 — расчет момента времени достижения КПУГ до начала поиска пл;

7. tk2:=t+R/Vpl, если t=0 — расчет времени выхода пл из радиуса ОВМ;

t:=tk1, если (tk2>tk1)

t:=tk5, если (tk5<=Ttreb)

t:=tk3, если (tk3<=tk5) и (tk3<=Ttreb)

t:=tk4, если (tk4<=Ttreb)

5. Накопление результатов моделирования

1. Pobn:=Nobn/N — вероятность обнаружения пл;

2. SMO:=SMO+Pobn — сумма вероятностей обнаружения пл;

3. SM2:=SM2+Pobn*Pobn — сумма квадратов вероятностей обнаружения пл.

6. Обработка результатов моделирования

1. MO:=SMO/N; 2. D:=(SM2-(SMO*SMO)/N)/(N-1); 3. E:=Tb*sqrt (D/N); 4. e0r:=E/MO или e0r:=E/(1-MO); 5. Nneobh:=Tb*Tb*D/(MO*MO*e0*e0);

7. Блок-схема обобщенного алгоритма модели

8. Блок-схемы процедур алгоритма

Процедура определения факта обнаружения пл и накопление количества обнаружения пл

Процедура задания и вычисления начальных условий одного малого прогона

Процедура задания и вычисления начальных условий одного большого прогон

Процедура обработки результатов моделирования

Накопление результатов моделирования

Ввод исходных данных

Вывод результатов моделирования данных

9. Листинг программы

Program kursovik;

Uses Crt;

Const

R=15 000;{Padius OVM pl v metrax}

Ttreb=7200;{vrem9 odnogo progona}

Vpl=15;{skorost' pl}

Vnk=15;{skorost KPYG}

Y1=0.03/60;{intensivnost' obnary8eni9 KPUG Plom}

Y2=0.02/60;{intensivnost' obnary8eni9 PL sredstvami KPYG bez ykloneni9}

Y3=0.01/60;{intensivnost' obnary8eniz PL sredstvami KPYG s ykloneniem}

D0min=15 000;{min rasstoenie ot KPYG do TPK}

D0max=35 000;{max rasstoanie ot KPYG do TPK}

Nzad=100;

m=100;

Tb=1.96;

e0=0.05;

Var

Nobn, Nneobn, t, tk1,tk2,tk3,tk4,tk5,i, j, Nneobx:integer;

SMO, SM2, Pobn, MO, D, SKO, E, E1,E0r, D0: real;

Procedure NU_BP;

Begin

Nobn:=0;

Nneobn:=0;

End;

Procedure NU_MP;

Begin

t:=0;

D0:=D0min+Random (D0max-D0min);

tk1:=t+Round ((D0-R)/Vnk);

tk2:=t+Round (R/Vpl);

End;

Procedure Opr_obn;

Begin

If tk2>tk1 then

begin

t:=tk1;

tk3:=t-Round (ln (Random)/Y2);

tk5:=t-Round (ln (Random)/Y1);

If tk5

begin

If tk5

begin

t:=tk5;

tk4:=t-Round (ln (Random)/Y3);

If tk4<=Ttreb then

begin

t:=tk4;

Nobn:=Nobn+1

end

else

Nneobn:=Nneobn+1

end

else

Nneobn:=Nneobn+1

end

else

If tk3<=Ttreb then

begin

t:=tk3;

Nobn:=Nobn+1

end

else

Nneobn:=Nneobn+1;

end

else Nneobn:=Nneobn+1

END;

Procedure NAK_Rez;

Begin

Pobn:=Nobn/Nzad;

SMO:=SMO+Pobn;

SM2:=SM2+Pobn*Pobn;

End;

Procedure Obrabotka;

Begin

MO:=SMO/m;

D:=(SM2-(SMO*SMO)/m)/(m-1);

SKO:=sqrt (D);

E:=Tb*Sqrt (D/m);

If MO<=0.5 then

begin

e0r:=E/MO;

Nneobx:=Round ((Tb*Tb*d)/(MO*MO*e0*e0));

end

else

begin

e0r:=E/(1-MO);

Nneobx:=Round ((tb*Tb*D)/((1-MO)*(1-MO)*e0*e0));

end;

End;

Procedure VIVOD;

Begin

Writeln;

Write (' Mo=', mo:1:4,' SKO=', sko:1:4);

Write (' D=', D:1:6,' E=', e:1:5,' E0r=', e0r:1:6,' Nneobx=', Nneobx);

End;

begin

clrscr;

Randomize;

Writeln;

SMO:=0;

SM2:=0;

For i:=1 to m do

begin

NU_BP;

For j:=1 to Nzad do

begin

NU_MP;

Opr_obn;

end;

NAK_REZ;

Write (' ', Nobn);

end;

Obrabotka;

VIVOD;

End.

10. Результаты (количество обнаруженных пл из 100 прогонов)

MO=0,8172, SKO=0,0397, D=0,1 574, E=0,778, E_0r=0,42 537, N_neobh=72.

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

Определение эмпирического закона распределения случайной величины отклика и поверить гипотезу о его согласованности с нормальным законом распределения Пирсона и Колмогорова-Смирнова, если в результате 100 прогонов программной модели получена следующая статистика значений отклика:

y_min=47, y_max=67;

MO=56,91, SKO=4,38, D=19,21, E=0,86, E_0r=0,019;

оценка математического ожидания г=56,91;

Оценка дисперсии D=19,21;

у=4,38 — среднеквадратическое отклонение.

Разобьем интервал значений отклика на L равных интервалов. Если y — дискретная переменная, то каждый интервал должен соответствовать одному из возможных значений y_i. Если же переменная y — непрерывная, то весь диапазон ее значений разбивается на 5 — 20 равных интервалов. Один из возможных способов определения количества и размера интервалов:

L .

Определим количество прогонов N_эk (частоту), в которых значения отклика попадают в каждый k-й интервал.

Определим наблюдаемую относительную частоту попадания значения отклика в каждый k-й интервал.

Определим среднюю точку k-го интервала M_k, как среднее арифметическое значений отклика, попавших в данный интервал.

Построим гистограмму, откладывая по оси абсцисс значения отклика, а по оси ординат частоту N_эk или относительную частоту f_эk (в последнем случае гистограмма оказывается нормированной — по оси ординат максимальное значение равно единице).

После определения значений выборочного закона распределения F_э(y) (или функции плотности f_э(y)) и выдвигают нулевую гипотезу Н₀ о том, что полученное эмпирическое распределение согласуется с каким-либо теоретическим распределением.

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

Значения f_еk(y) для k-го интервала [y_k; y_k+1] непрерывной случайной величины определял по формуле где y_k — левая граница интервала, а y_k+1 — правая граница интервала, изначения функции распределения вероятностей для значений отклика, соответствующих левой и правой границе интервала, соответственно.

Значение функции распределения вероятностей можно вычислить с использованием пакета Excel путем выбора функции НОРМРАСП категории «Статистические» и заданием следующих параметров — X: значение y (например, 51), Среднее: = 56,91, Стандартное отклонение: = 4,38, Интегральная: 1 (интегральная функция распределения).

Тогда значения f_еk(y) для нормального закона распределения примут вид, в приведенной выше таблице. Значения в таблице округлили.

Значения f_еk(y), например, для нормального закона можно вычислить, подставляя значения y_k в формулу плотности вероятностей

.

Значения кумулятивных (накопленных) вероятностей F_эk(y) и F_еk(y) находятся как суммы f_эk(y) или f_еk(y) в данном интервале и во всех интервалах, предшествующих данному.

Значения N_еk(y) находятся путем умножения значений f_еk(y) на N =100.

Проверка гипотезы о незначимом отличии эмпирического и теоретического распределения значений отклика.

При проверке гипотезы по критерию Пирсона (хи-квадрат) необходимо пользоваться значениями N_эk(y) и N_еk(y),

В каждом интервале значение N_эk(y) или N_еk(y) должно быть не меньше 5, в противном случае интервал объединяется с соседним.

Проверить гипотезу по критерию согласия Пирсона. Для этого вычислим значение Определим количество степеней свободы = L — r — l=5−2-1=2, где L — количество интервалов (после объединения), r — количество параметров теоретического закона распределения (для нормального распределения r = 2, так, как параметров нормального закона распределения два: математическое ожидание и дисперсия ²);

Из таблицы выберем значение ²_крит по доверительной вероятности =0,95 и количеству степеней свободы =2. Табличное значение будет равно 5,99

Таким образом, табличное значение при в=0.95 и v=2 ч²_крит=7,81, а полученный ч²_крит=1,182. Так как ч²_крит? ч²_, то гипотеза согласованности эмпирического и теоретического распределения случайной величины поверкой по критерию Пирсона не опровергается. Проверим гипотезу по критерию согласия Колмогорова-Смирнова. Для этого вычислить значения модулей разностей для каждого из интервалов:K=|F_эk(y) — F_еk(y)|;

1. | F_эk1(y) — F_еk1(y)|=0,029;

2. | F_эk2(y) — F_еk2(y)|=0,031;

3. | F_эk3(y) — F_еk3(y)|=0,033;

4. | F_эk4(y) — F_еk4(y)|=0,008;

5. | F_эk5(y) — F_еk5(y)|=0,01.

Определим значение К=max =0,033.

Из таблицы «Критические числа Колмогорова-Смирнова», по уровню значимости = 0,05 для проверки единичной выборки и количеству степеней свободы =100 определить критическое значение K_крит по формуле K_крит=.Сравним К=0,033 и К_крит.=0,136. Так как K K_крит, то гипотеза Н₀ о согласованности эмпирического и теоретического распределения случайной величины проверкой по критерию согласия Колмогорова-Смирнова не опровергается.

Представлен график зависимостей плотности функции распределения и самой функции распределения на интервалах:

Корреляционный анализ

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

Корреляционный анализ результатов моделирования сводится к оценке разброса значений Y относительно среднего значения (отклика) в зависимости от разброса значений Х относительно среднего значения (фактора), т. е. к оценке силы корреляционной связи между случайной величиной Y отклика и случайной величиной фактора X.

— если =0, то это свидетельствует о взаимной независимости случайных переменных? и ?, исследуемых при моделировании;

— если | |=1, то имеет место функциональная (т. е. не стохастическая) линейная зависимость вида у=b0+b1х, причем если r??>0, то говорят о положительной корреляции, т. е. большие значения одной случайной величины соответствуют большим значениям другой;

— если 0<| |<1, то этот случай соответствует, либо наличию ли-нейной корреляции с рассеянием, либо наличию нелинейной корреляции результатов моделирования.

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

причем случайная величина функции w() приближенно подчиняется нормальному распределению.

Для этого нам необходимо в программу ввести дополнительные переменные:

Sxy = ?(x_k*y_k); =MO; =MOx;

SMOx, SM2, SM2x, Sxy, MOx, SD0, Xk, r1, w:real;

тогда наша формула в Pascal будет выглядеть следующим образом:

1.Процедура накопления результатов:

Procedure NAK_Rez;

Begin

Pobn:=Nobn/Nzad;

SMO:=SMO+Pobn;

SM2:=SM2+Pobn*Pobn;

Xk:=SD0/Nzad;

SMOx:=SMOx+Xk;

SM2x:=SM2x+Xk*Xk;

Sxy:=Sxy+Xk*Pobn

End;

2.Процедура обработки результатов моделирования:

Procedure Obrabotka;

Begin

MO:=SMO/m;

MOx:=SMOx/m;

D:=(SM2-(SMO*SMO)/m)/(m-1);

SKO:=sqrt (D);

E:=Tb*Sqrt (D/m);

If MO<=0.5 then

begin

e0r:=E/MO;

Nneobx:=Round ((Tb*Tb*d)/(MO*MO*e0*e0));

end

else

begin

e0r:=E/(1-MO);

Nneobx:=Round ((tb*Tb*D)/((1-MO)*(1-MO)*e0*e0));

end;

r1:=(Sxy-m*MOx*MO)/sqrt ((SM2x-m*MOx*MOx)*(SM2-m*MO*MO));

w:=(sqrt (m-3)*ln ((1+r1)/(1-r1)))/2

End;

2.Процедура вывод результатов моделирования:

Procedure VIVOD;

Begin

Writeln;

Write (' Mo=', mo:1:4,' SKO=', sko:1:4);

Write (' D=', D:1:6,' E=', e:1:5,' E0r=', e0r:1:6,' Nneobx=', Nneobx);

Writeln;

Write (' r1=', r1:1:4,' w=', w:1:4);

End;

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

Procedure NU_BP;

Begin

Nobn:=0;

Nneobn:=0;

SD0:=0

End;

Procedure NU_MP;

Begin

t:=0;

D0:=D0min+Random (D0max-D0min);

tk1:=t+Round ((D0-R)/Vnk);

tk2:=t+Round (R/Vpl);

SD0:=SD0+D0;

End;

4.В теле программы необходимо обнулять накопление:

begin

clrscr;

Randomize;

Writeln;

SMO:=0;

SM2:=0;

SMOx:=0;

SM2x:=0;

Sxy:=0;

For i:=1 to m do

begin

NU_BP;

For j:=1 to Nzad do

begin

NU_MP;

Opr_obn;

end;

NAK_REZ;

Write (' ', Nobn);

end;

Obrabotka;

VIVOD;

End.

В итоге мы получили результат:

R₁=-0,5118, W=-9,7407 — коэффициент корреляции.

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

12. Дисперсионный анализ результатов моделирования

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

Для того, что бы выявить влияние фактора на отклик, необходимо определить: МО и D на верхнем и нижнем уровнях, вычислить D_му, МО_му и вычислить D0 по формуле D0=(D_v+D_n)/2, и не забываем вводить в исходную программу L (200 на нижнем уровне, 0 на верхнем уровне и 100 между уровнями.

Листинг измененной программы для фактора D0:

Const

R=15 000;{Padius OVM pl v metrax}

Ttreb=7200;{vrem9 odnogo progona}

Vpl=15;{skorost' pl}

Vnk=15;{skorost' nk}

Y1=0.03/60;{intensivnost' obnary8eni9 KPUG Plom}

Y2=0.02/60;{intensivnost' obnary8eni9 PL sredstvami KPYG bez ykloneni9}

Y3=0.01/60;{intensivnost' obnary8eniz PL sredstvami KPYG s ykloneniem}

D0minN=5000; минимальное значение D0 на нижнем уровне

D0maxN=25 000; максимальное значение D0 на нижнем уровне

D0minV=25 000; минимальное значение D0 на верхнем уровне

D0maxV=45 000; максимальное значение D0 на верхнем уровне

Nzad=100;

m=200;

l=100;

Tb=1.96;

e0=0.05;

{====================================}

Var

Nobn, Nneobn, t, tk1,tk2,tk3,tk4,tk5,i, j, Nneobx, D0min, D0max: longint;

SMO, SM2, Pobn, MO, D, SKO, E, E0r, D0: real;

далее изменения в самом теле программы

begin

clrscr;

Randomize;

Writeln;

SMO:=0;

SM2:=0;

For i:=1 to m do

begin

If i<=l then

Begin

D0min:=D0minN;

D0max:=D0maxN

End else

Begin

D0min:=D0minV;

D0max:=D0maxV

end;

NU_BP;

For j:=1 to Nzad do

begin

NU_MP;

Opr_obn;

end;

NAK_REZ;

Write (' ', Nobn);

end;

Obrabotka;

VIVOD;

End.

Далее находим F= D_му/ D0, из таблицы находим F_таб=1,23.

Фaктор будет влиять на отклик если F>F_кр.

D0 — расстояние от КПУГ до ТПК;

MO_н=0,7725 D_н=0,1 665

MO_в=0,1908 D_в=0,1 559

MO_му=0,4736 D_му=0,86 982

D0=0,1 612

F=5,3959 F> F_кр

F_кр=1,23. Фактор D0 влияет на отклик.

V_nk — скорость КПУГ

MO_н=0,4881 D_н=0,2 074

MO_в=0,6371 D_в=0,2 479

MO_му=0,1831 D_му=0,6 115

D0=0,002

F=2,686 F> F_кр

Fкр=1,23. Фактор V_nk влияет на отклик.

Y₂ — интенсивность обнаружения пл средствами КПУГ без уклонения

MO_н=0,5212 D_н=0,1

MO_в=0,7201 D_в=0,29

MO_му=0,2070 D_му=0,2 861

D0=0,0025

F=1,312 F> F_кр

F_кр=1,23. Фактор Y₂ влияет на отклик.

Y₃ — интенсивность обнаружения пл средствами КПУГ с уклонения

MO_н=0,5305 D_н=0,2 560

MO_в=0,7453 D_в=0,1 740

MO_му=0,2186 D_му=0,2 807

D0=0,215

F=70,3 F> F_кр

F_кр=1,23. Фактор Y₃ влияет на отклик.

13. Регрессионный анализ

Алгоритм проведения регрессионного анализа

— В результате проведения корреляционного и дисперсионного анализа определить факторы x_i, которые значимо влияют на отклик.

Для того, чтобы уравнение регрессии получилось адекватным (соответствующим) имитационной модели в нем должны быть учтены все факторы. Однако, в этом случае размерность плана полного факторного эксперимента может оказаться очень большой, так как если даже принять число уровней всех факторов одинаковым, то количество экспериментов будет равно N_c=q^k, где k — количество факторов, q — количество уровней каждого фактора.

Поэтому для уменьшения размерности эксперимента следует:

- либо уменьшить количество факторов, пренебрегая теми факторами, которые не значимо влияют на отклик;

- либо уменьшить количество уровней каждого фактора;

- либо проводить не полный, а дробный факторный эксперимент.

— Определить количество и значения каждого уровня каждого фактора x_i. Если каждый фактор будет изменяться только на двух уровнях, то определить значения верхнего x_iВ и нижнего уровня x_iН каждого фактора.

— Построить матрицу планирования и план полного факторного эксперимента (ПФЭ) или дробного факторного эксперимента (ДФЭ).

План проведения полного факторного эксперимента.

Обозначим наши факторы V_nk (скорость выхода КПУГ), Y₂ (интенсивность обнаружения пл средствами КПУГ без уклонения), Y₃ (интенсивность обнаружения пл средствами КПУГ с уклонением) и D₀ (расстоянием между КПУГ и ТПК) на Х₁, Х₂, Х₃ и Х₄ соответственно.

Матрица планирования Значения всех факторов на пяти уровнях сведены в таблицу

Вычислить коэффициенты b₀, b₁, b₂, b₃, …, b₁₂, b₁₃, b₂₃, …, b₁₂₃ и т. д. в размерности отклика.

В результате уравнение регрессии примет вид:

где, , N — количество экспериментов (т. е. строк в ПФЭ или ДФЭ),

— кодированное значение каждого фактора (элемент матрицы планирования), при двух уровнях может принимать значение «-1», если фактор находится на нижнем уровне, или «+1», если фактор находится на верхнем уровне. Коэффициент b₀ всегда будет положительным, а остальные коэффициенты b_i могут получаться как положительными, так и отрицательными.

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

?х₁=5; ?х₂=0,18/600; ?х₃=0,019/600; ?х₄=1000;

b0= 0,5 678 125 b1= 0,18 725

b2= 0,38 275 b3= 0,51 375

b4= -0,672 875 b1,2= -0,256 375

b1,3= 0,171 125 b1,4= 0,14 725

b2,3= -0,49 125 b2,4= -0,0061

b3,4= -0,0074 b1,2,3=-0,1 975

b1,2,4=0,388 125 b1,3,4=0,27 625

b1,2,3,4= -0,0013 b2,3,4= 0,4 201 875

Вычислить коэффициенты k₀, k₁, k₂, k₃, …, k₁₂, k₁₃, k₂₃, …, k₁₂₃ и т. д. в размерности факторов.

, , ,

где, x_iВ — фактическое значение верхнего уровня фактора, x_iН — фактическое значение нижнего уровня фактора. Если фактор не детерминированный, а стохастический, то за x_iВ и x_iН принимаются математические ожидания верхнего и нижнего уровня соответственно.

В результате уравнение регрессии примет вид:

где dx_i = x_iф — x_i0 — отклонение фактического значения фактора x_i от его среднего значения, x_iф — фактическое значение фактора x_i в размерности этого фактора, x_i0 = x_iВ — x_i = x_iН + x_i — среднее значение фактора x_i.,, ,. В «Excel» вычисляем y. В исходную программу вводим новые переменные DR и SR. Вычислить SR по формуле SR=SR+(Pobn-y)* (Pobn-y), а в обработке вычисляем DR=SR/(M-1). В телепрограммы необходимо обнулить SR.

k₀=0,5 678 125 k₁=0,3 745

k₂=127,5 833 333 k₃=342,5

k₄=-0,3 364 k_1,2=-17,9 166 667

k_1,3=22,8166 k_1,4=0,1 472

k_2,3=-109 166,6667 k_2,4=-0,10 166 667

k_3,4=0,937 k_1,2,3=3,19444E-05

k_1,2,4=8,36806E-05 k_1,3,4=0,298 177

k_1,2,3,4=9,29 977E-05 k_2,3,4=0,194 531

Подробно опишем как в «Excel» найти все коэффициенты и y.

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

Для проверки адекватности уравнения регрессии имитационной модели следует выполнить следующую последовательность действий:

провести пробный эксперимент (на N прогонов) с любыми значениями факторов из области их определения (аналогично новой строке в матрице планирования ПФЭ);

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

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

определить дисперсию разности откликов из эксперимента и из уравнения регрессии по формуле ;

по критерию согласия Фишера проверить гипотезу H₀ о незначимом отличии дисперсии воспроизводимости (соответствующей ошибке эксперимента) от дисперсии разности (свидетельствующей об адекватности уравнения регрессии).

Если гипотеза H₀ принимается, то уравнение регрессии адекватно модели. В противном случае следует корректировать уравнение: например, путем проведения ПФЭ вместо ДФЭ, увеличения уровней факторов, учета степенных коэффициентов и т. д.

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

Условно обозначим верхний уровень за 1, нижний — -1, между нижним и нулевым — -½, между верхним и нулевым — ½.

Сравнив F c F_крит видим, что уравнение регрессии неадекватно во всех точках, поэтому необходим МПО.

14. Покоординатный спуск

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

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

Х₁=20±1, Х₂=0,38/600, Х₃=0,19/600, Х₄=23 000±100 000

15. Наискорейший спуск

Для начала необходимо определить нижний и верхний уровни, не далеко отстоящие от нулевого уровня, а затем провести полный факторный эксперимент. На нулевом уровне значение отклика P_obn=0,5539.

Первая итерация:

Наибольшее значение отклика: 0,8545

Комбинация значений факторов: Х₁=19±1, Х₂=0,03/60, Х₃=0,015/60, Х₄=24 000±10 000 это будет новый нулевой уровень для следующей итерации.

Вторая итерация:

Наибольшее значение отклика: 0,9215

Комбинация значений факторов:

Х₁=21±1, Х₂=0,034/60, Х₃=0,019/60, Х₄=23 500±10 000 это будет новый нулевой уровень для следующей итерации.

Третья итерация:

Наибольшее значение отклика: 0,9423

Комбинация значений факторов: Х₁=21±1, Х₂=0,036/60, Х₃=0,021/60, Х₄=230 000±10 000

Таким образом, комбинация этих факторов даёт наилучшее значение отклика. В последующем, если повышать фактор Y₂ и Y₃, а фактор D₀ уменьшать, то вероятность будет расти до 1. Фактор Vnk повышать больше не рекомендуется, т.к. скорость кораблей на данный момент времени является предельным.

Заключение

Для обеспечения эффективности поиска подводной лодки кораблями КПУГ при возможности уклонения пл, разработана имитационная модель процесса обнаружения пл противника.

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

Проведены корреляционный и дисперсионный анализы результатов моделирования, показавшие наличие значимой и неслучайной зависимости вероятности обнаружения пл кораблями КПУГ от расстояние между КПУГ и ТПК, скорости выхода кораблей КПУГ в позицию начала поиска, интенсивности обнаружения пл средствами КПУГ без уклонения пл, интенсивности обнаружения пл средствами КПУГ при уклонении пл.

В результате корреляционного и дисперсионного анализов выявлены управляемые факторы для проведения регрессионного анализа.

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

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

Таким образом, комбинация 4 факторов даёт наилучшее значение отклика при повышении интенсивностей обнаружения пл средствами КПУГ при уклонении и без уклонения пл. Так же, с уменьшением дистанции между КПУГ и ТПК вероятность обнаружения пл будет расти. Скорость кораблей КПУГ увеличивать не рекомендуется, т.к. скорость кораблей в имитационной модели принята максимально возможной. Значения всех факторов будут таковыми: