Test launch of a Peacekeeper ICBM

Test launch of a Peacekeeper ICBM

About ICBM

The LGM-118 Peacekeeper, also known as the MX missile (for Missile-eXperimental), was a land-based ICBM deployed by the United States starting in 1986. The Peacekeeper was a MIRV missile that could carry up to 10 re-entry vehicles, each armed with a 300-kiloton W87 warhead in a Mk.21 reentry vehicle (RV). A total of 50 missiles were deployed starting in 1986, after a long and contentious development program that traced its roots into the 1960s. Since 2005 has been phased out of service as obsolete. Spare parts of the second stage are used for Antares project.

The Peacekeeper program began in 1971 as the missile experimental (MX) system as a way to increase the U.S. counterstrike capabilities against the Soviet Union, which at the time was focusing on constructing hardened shelters and missile defense systems. The Peacekeeper employed an advanced guidance system, a Multiple Independently Targetable Reentry Vehicle (MIRV) system of a dozen warheads, and a cold launch system to allow for silo reuse. The LGM-118 Peacekeeper missile started development in 1971. The full-scale development of the Peacekeeper began in 1974 and the first flight test occurred in 1983. While there were plans for 100 of the silo-based missiles, that number was cut to 50 in 1984. The remaining missiles were intended for deployment on railcar launchers. The Peacekeeper was first deployed in 1986 at Warren Air Force Base in Wyoming. A total of 114 missiles were produced by the end of 1988.

For more see:

1.https://missilethreat.csis.org/missile/lgm-118-peacekeeper-mx/

2.https://en.m.wikipedia.org/wiki/LGM-118_Peacekeeper

3.http://nuclearweaponarchive.org/Usa/Weapons/Mx.html

SLBM Data

load("slbm.dat")
library(DT)
datatable(slbm)
library(ggplot2)
library(GGally)
#Here we use function from 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

Bayesian Linear Regression Model

library(rstanarm)
## Loading required package: Rcpp
## rstanarm (Version 2.17.3, packaged: 2018-02-17 05:11:16 UTC)
## - Do not expect the default priors to remain the same in future rstanarm versions.
## Thus, R scripts should specify priors explicitly, even if they are just the defaults.
## - For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores())
## - Plotting theme set to bayesplot::theme_default().
library(bayesplot)
## This is bayesplot version 1.4.0
## - Plotting theme set to bayesplot::theme_default()
## - Online documentation at 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).

ICBM Operational Range Estimation

Here we use open source data (see above) for the Bayesian linear regression model.

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")

ICBM Reliability

As far as we have got no relevant data about MX flight test statistics (success, fail), we can use MTBF metric i.e mean time before failure. There is one open source that gives us a clue to solve this problem - “ICBM MODERNIZATION. Availability Problems and Flight Test Delays in Peacekeeper Program” that is March 1989 GAO Report. Guidance system availability in terms of MTBF metric is crucial to the total ICBM reliability during flight time. As we know \(P_{GS}=e^{-T/MTBF}\) where \(T\) - is the time of ICBM flight from launch pad (silo) to the target.

For more about MTBF see: https://en.wikipedia.org/wiki/Mean_time_between_failures

According to the GAO Report (see above) we can assume: \(MTBF=2839\) and \(T=0.5\).

“As of January 3,1989, the program office reported that the MGCS cumulative mean time between recycles was 2,444 hours, or about 19 percent less than the planned level of 3,000 hours. The mean time between recycle rate has improved, however, with a recycle rate of 2,839 hours achieved during the 10-month period of March through December 1988.” - quote from the GAO Report, page 22.

ICBM Warhead Kill Probability

Estimating SSKP by Lethal Radius

The Lethal Radius (LR) is defined as the distance from the point of the nuclear explosion that the warhead will be able to destroy its target. The formula for LR (in meters) as a function of Yield (Y–in Megatons) and silo hardness (H –in overpressure pounds per square inch or psi) is given by:

Lethal Radius im meters

Lethal Radius im meters

Nuclear blast overpressure

Nuclear blast overpressure

The Effects of Nuclear Weapons Samuel Glasstone and Philip J. Dolan, Third edition, 1977

The height of burst and energy yield of the nuclear explosion are important factors in determining the extent of damage at the surface.

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 reliability during the flight time = 0.9998239
PGS50 <- PGS^50
cat("All 50 MX ICBMs are reliabale during the flight time =", PGS50)
## All 50 MX ICBMs are reliabale during the flight time = 0.9912327
PK1 <- PGS*SSKP(0.350,10000,90)
cat ("Probability of hardend silo (10000 psi) kill by one warhead = ", PK1)
## Probability of hardend silo (10000 psi) kill by one warhead =  0.8546648
PK2 <- 1-(1-PK1)^2
cat ("Probability of hardend silo (10000 psi) kill by two warheads = ", PK2)
## Probability of hardend silo (10000 psi) kill by two warheads =  0.9788777

Conclusions

1.Applying Bayesian linear regression model for SLBM data (M,P,D,L,S,W) we can produce posterior distribution sample given MX data (L=21.1,D=2.34,S=3,W=2,M=87750,P=3950) with 80% credible interval \(P(8883\le Range_{MX}\le17403)=0.8\) and the mean \(M[Range_{MX}]=12754 km\).
2.MX had operational Range and Reliability (0.9998239) corresponding to ICBM which could hit any hardened target (\(H\le10000psi\)) worldwide \(P_{kill,1,0.35}=0.855\) and \(P_{kill,2,0.35}=0.979\).