рекомендации

вторник, 18 августа 2020 г.

Раскрываем тайны МБР LGM-118 Peacekeeper (MX) с помощью R



Немного об МБР

LGM-118 Peacekeeper, также известная как MX (Missile-eXperimental), была наземной МБР, развернутой Соединенными Штатами начиная с 1986 года. Peacekeeper была ракетой с разделяющимися боеголовками, которая могла нести до 10 боеголовок W87 мощностью 300 кт на платформе разведения Mk.21 (RV). В общей сложности были развернуты 50 ракет, после долгой и спорной программы разработки, которая уходила корнями в 1960-е годы. С 2005 года она была выведена из эксплуатации как устаревшая. Некоторые части второй ступени используются для проекта Antares.

Программа «Peacekeeper» началась в 1971 году как экспериментальная система ракет (MX), способная увеличить возможности США по противодействию Советскому Союзу, который в то время занимался строительством укрепленных шахт и систем противоракетной обороны. Peacekeeper использовал передовую систему наведения, систему индивидуального нацеливания дюжины боеголовок и систему холодного пуска для повторного использования шахты. Проектирование ракеты LGM-118 Peacekeeper было начато в 1971 году. Полномасштабная разработка Peacekeeper началась в 1974 году, а первые летные испытания состоялись в 1983 году. Несмотря на то, что было запланировано развертывания 100 ракет в шахтных пусковых установках, это число было сокращено до 50 в 1984 году. Остальные ракеты предназначались для размещения на железнодорожных пусковых установках. Peacekeeper был впервые развернут в 1986 году на базе ВВС Warren в Вайоминге. Всего к концу 1988 года было произведено 114 ракет.

Более подробная информация:


SLBM Data

load("slbm.dat")
library(DT)
datatable(slbm)

Show entries
Search:
  S M L D R P W
1 1 5400 10.4 0.58 150 975 1
2 1 13700 11.8 1.3 650 1597 1
3 1 19650 14.2 1.3 1420 1179 1
4 1 14200 8.89 1.5 2400 650 1
5 1 14200 8.89 1.5 3000 650 1
6 1 14200 9.65 1.5 3000 650 2
7 2 33300 13 1.8 7800 1100 1
8 2 33300 13 1.8 9100 1100 1
9 2 26900 10.6 1.54 3900 450 1
10 2 35300 14.1 1.8 6500 1600 2
Showing 1 to 10 of 26 entries

library(ggplot2)
library(GGally)
#Здесь мы используем функцию из https://www.r-bloggers.com/multiple-regression-lines-in-ggpairs/
my_fn = function(data, mapping, ...){
  p = ggplot(data = data, mapping = mapping) + 
    geom_point() + 
    geom_smooth(method=loess, fill="red", color="red", ...) +
    geom_smooth(method=lm, fill="blue", color="blue", ...)
  p
}

g = ggpairs(slbm,columns = 2:6, lower = list(continuous = my_fn))
g



Байесовская модель линейной регрессии

library(rstanarm)
## Загрузка необходимого пакета: Rcpp
## rstanarm (Version 2.17.3, packaged: 2018-02-17 05:11:16 UTC)
#Не ожидайте, что приоритеты по умолчанию останутся такими же в будущих версиях rstanarm
##Поэтому скрипты R должны явно указывать приоритеты, даже если они являются просто значениями по умолчанию.
##Для выполнения на многоядерном процессоре с избытком оперативной памяти:

## options(mc.cores = parallel::detectCores())
## - Тема графика bayesplot::theme_default().
library(bayesplot)
## Это bayesplot версии 1.4.0
## - Тема графика bayesplot::theme_default()
## - Онлайн-документация на mc-stan.org/bayesplot
options(mc.cores = parallel::detectCores())
set.seed(12345)
fit.slbm.bs=stan_glm(sqrt(R)~S+D+L+W+log(M)+log(P),
                        data=slbm,chains=2,iter=10000)
fit.slbm.bs
## stan_glm
##  family:       gaussian [identity]
##  formula:      sqrt(R) ~ S + D + L + W + log(M) + log(P)
##  observations: 26
##  predictors:   7
## ------
##             Median MAD_SD
## (Intercept) -43.1  105.2 
## S             9.5    5.2 
## D            36.7   19.5 
## L             0.8    1.7 
## W             3.3    5.6 
## log(M)        6.9   15.0 
## log(P)       -7.5    6.0 
## sigma        10.2    1.7 
## 
## Sample avg. posterior predictive distribution of y:
##          Median MAD_SD
## mean_PPD 69.0    2.8  
## 
## ------
## For info on the priors used see help('prior_summary.stanreg').
plot(fit.slbm.bs)











posterior_vs_prior(fit.slbm.bs)
## 
## Drawing from prior...

fit.slbm.bs2=as.array(fit.slbm.bs)
mcmc_hist(fit.slbm.bs2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.



mcmc_trace(fit.slbm.bs2)







summary(fit.slbm.bs)
## 
## Model Info:
## 
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      sqrt(R) ~ S + D + L + W + log(M) + log(P)
##  algorithm:    sampling
##  priors:       see help('prior_summary')
##  sample:       10000 (posterior sample size)
##  observations: 26
##  predictors:   7
## 
## Estimates:
##                 mean   sd     2.5%   25%    50%    75%    97.5%
## (Intercept)    -41.8  109.6 -259.0 -112.4  -43.1   29.2  176.0 
## S                9.5    5.3   -1.0    6.0    9.5   13.0   20.1 
## D               36.8   20.1   -3.6   23.5   36.7   49.8   75.9 
## L                0.8    1.8   -2.7   -0.3    0.8    2.0    4.4 
## W                3.2    5.8   -8.2   -0.5    3.3    6.9   14.6 
## log(M)           6.8   15.6  -24.1   -3.2    6.9   17.0   37.6 
## log(P)          -7.4    6.3  -19.9  -11.5   -7.5   -3.4    5.0 
## sigma           10.5    1.8    7.6    9.2   10.2   11.5   14.6 
## mean_PPD        69.0    2.9   63.3   67.1   69.0   70.9   74.9 
## log-posterior -110.5    2.4 -116.2 -111.8 -110.1 -108.8 -107.1 
## 
## Diagnostics:
##               mcse Rhat n_eff
## (Intercept)   1.7  1.0   4072
## S             0.1  1.0   5834
## D             0.3  1.0   4403
## L             0.0  1.0   5603
## W             0.1  1.0   6586
## log(M)        0.2  1.0   3928
## log(P)        0.1  1.0   6953
## sigma         0.0  1.0   5258
## mean_PPD      0.0  1.0  10000
## log-posterior 0.0  1.0   3008
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Оценка дальности действия МБР

Здесь мы используем данные с открытым исходным кодом (см. выше) для модели байесовской линейной регрессии.

set.seed(12345)
slbm.pp = rstanarm::posterior_predict(fit.slbm.bs,
          newdata = data.frame(L=21.1,D=2.34,S=3,W=2,
                               M=87750,P=3950),seed=12345)

Range.km = slbm.pp^2
Range.km = Range.km[1:10000,]
quantile(Range.km,probs = c(0.1,0.5,0.9))
##       10%       50%       90% 
##  8883.001 12754.707 17403.278
mean(Range.km>15000)
## [1] 0.2525
ggplot(data=as.data.frame(Range.km), aes(Range.km)) + 
  geom_histogram(bins = 30,col="black",fill="green") + 
  geom_vline(xintercept = mean(Range.km), color = "red") + 
  geom_errorbarh(aes(y=-5, xmin=quantile(Range.km,0.1),
                     xmax=quantile(Range.km,0.9)), 
                 data=as.data.frame(Range.km), col="#0094EA", size=3) +
  ggtitle(label="Probability density of MX Range")



Надежность МБР

Поскольку у нас нет соответствующих данных о статистике летных испытаний MX (успех, неудача), мы можем использовать метрику MTBF, т.е. среднее время до отказа. Существует один открытый источник, который дает нам ключ к решению этой проблемы - «ICBM MODERNIZATION. Availability Problems and Flight Test Delays in Peacekeeper Program» - доклад GAO за март 1989 года. Надежность системы наведения в показателях MTBF имеет решающее значение для общей надежности МБР во время полета. Как известно, PGS = e − T / MTBF, где T - время полета МБР от стартовой площадки (шахты) к цели.

Для получения дополнительной информации о MTBF: 

Согласно отчету GAO (см. выше) мы можем принять: MTBF = 2839 и T = 0,5.

Вероятность поражения цели боеголовкой МБР

Оценка вероятности поражения цели одной боеголовкой (SSKP) по радиусу поражения

Радиус поражения (LR) определяется как расстояние от эпицентра ядерного взрыва до цели, при котором боеголовка обеспечивает уничтожение цели. Формула для LR (в метрах) как функция мощности (Y - в мегатоннах) и прочности шахты (H - в фунтах избыточного давления на квадратный дюйм или фунт / кв. дюйм) определяется как:

Радиус поражения в метрах



Избыточное давление при ядерном взрыве
  
The Effects of Nuclear Weapons Samuel Glasstone and Philip J. Dolan, Third edition, 1977

Высота взрыва и выход энергии ядерного взрыва являются важными факторами при определении степени разрушений на поверхности.

LRM = function(Y,H){
  LRM1 = 4540*(Y^(1/3))/(H^(1/3))
  LRM2 = (sqrt(1+2.79/H)+1.67/sqrt(H))^(2/3)
  return(LRM1*LRM2)
}

SSKP = function(Y,H,CEP) {
  LRM1 = 4540*(Y^(1/3))/(H^(1/3))
  LRM2 = (sqrt(1+2.79/H)+1.67/sqrt(H))^(2/3)
  return(1-0.5^((LRM1*LRM2)^2/CEP^2))
}

PGS = exp(-0.5/2839)
cat("MX reliability during the flight time =", PGS)
## Надежность MX во время полета = 0.9998239
PGS50 = PGS^50
cat("All 50 MX ICBMs are reliabale during the flight time =", PGS50)
## Вероятность надежной работы всех 50 МБР MX во время полета = 0.9912327
PK1 = PGS*SSKP(0.350,10000,90)
cat ("Probability of hardend silo (10000 psi) kill by one warhead = ", PK1)
## Вероятность поражения защищенной шахты (10000 psi) одной боеголовкой =  0.8546648
PK2 = 1-(1-PK1)^2
cat ("Probability of hardend silo (10000 psi) kill by two warheads = ", PK2)
## Вероятность поражения защищенной шахты (10000 psi) двумя боеголовками =  0.9788777

Выводы

1. Используя байесовскую модель линейной регрессии для данных БРПЛ (M, P, D, L, S, W), мы можем получить выборку апостериорного распределения с данными MX (L = 21,1, D = 2,34, S = 3, W = 2, M = 87750, P = 3950) с вероятностью 80% P (8883≤RangeMX≤17403) = 0,8 и средним значением M [RangeMX] = 12754 км.
2. МБР MX имела дальность и надежность (0,99998239), соответствующие МБР, которая могла поразить любую защищенную цель (H≤10000psi) по всему миру: Pkill, 1,0,35 = 0,855 и Pkill, 2,0,35 = 0,979.

Комментариев нет:

Отправить комментарий