Test launch of a Peacekeeper 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/
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
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).
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")
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.
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
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
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\).