Demystifying Minuteman III
Alexander Levakov, Senior Research Fellow, Ph.D
November 21, 2018
Minuteman III in silo
Minuteman III
The LGM-30 Minuteman is a U.S. land-based intercontinental ballistic missile (ICBM), in service with the Air Force Global Strike Command. As of 2018, the LGM-30G Minuteman III version[a] is the only land-based ICBM in service in the United States. The current U.S. force consists of 399 Minuteman-III missiles as of September 2017.
For more see:https://en.wikipedia.org/wiki/LGM-30_Minuteman
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/bayesplotoptions(mc.cores = parallel::detectCores())
load("slbm.dat")
fit.slbm.bs<-stan_glm(sqrt(R)~S+D+L+W+log(M)+log(P),
data=slbm,chains=4,iter=10000,seed=12345)
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) -40.6 106.7
## S 9.4 5.1
## D 37.4 19.8
## L 0.9 1.7
## W 3.3 5.6
## log(M) 6.5 15.0
## log(P) -7.3 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: 20000 (posterior sample size)
## observations: 26
## predictors: 7
##
## Estimates:
## mean sd 2.5% 25% 50% 75% 97.5%
## (Intercept) -40.0 109.1 -254.9 -112.3 -40.6 31.7 175.7
## S 9.4 5.4 -1.2 5.9 9.4 12.9 20.2
## D 37.4 20.3 -2.6 24.1 37.4 50.9 77.4
## L 0.9 1.8 -2.7 -0.3 0.9 2.0 4.4
## W 3.2 5.8 -8.5 -0.5 3.3 7.1 14.8
## log(M) 6.4 15.6 -24.3 -3.6 6.5 16.6 37.0
## log(P) -7.3 6.2 -19.5 -11.3 -7.3 -3.3 5.0
## sigma 10.4 1.8 7.6 9.2 10.2 11.4 14.6
## mean_PPD 69.0 2.9 63.2 67.1 69.0 70.9 74.7
## log-posterior -110.5 2.4 -116.3 -111.8 -110.1 -108.8 -107.1
##
## Diagnostics:
## mcse Rhat n_eff
## (Intercept) 1.2 1.0 8809
## S 0.0 1.0 12972
## D 0.2 1.0 9063
## L 0.0 1.0 12118
## W 0.1 1.0 11991
## log(M) 0.2 1.0 8316
## log(P) 0.1 1.0 13099
## sigma 0.0 1.0 10055
## mean_PPD 0.0 1.0 19069
## log-posterior 0.0 1.0 6234
##
## 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).Minuteman III Operational Range
Minuteman III MIRV path
slbm.pp <- posterior_predict(fit.slbm.bs,
newdata = data.frame(L=18.2,D=1.85,S=3,W=1,
M=34467,P=1000),seed=12345)
#Note. Payload is an expert estimate for the model
Range.km <- slbm.pp^2
Range.km <- Range.km[1:10000,]
quantile(Range.km,probs = c(0.05,0.5,0.95))## 5% 50% 95%
## 4327.626 8674.922 14653.407library(grid)
Range.km <- slbm.pp^2
Range.km <- Range.km[1:10000,]
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.05),
xmax=quantile(Range.km,0.95)),
data=as.data.frame(Range.km), col="#0094EA", size=3) +
ggtitle(label="Probability density of Minuteman III range",
subtitle = "P(4327<Range.km<14653)=0.9")
Minuteman survivability
Here we apply formula from Bruice G.Blair “Strategic Command and Control. Redefining the Nuclear Threat”, page 305.
tkprob <- function(OAR,Y,H,CEP){
#OAR - Overall reliability of attacking ICBM
#Y - RV Yield, MT
#CEP - Circular error probable of attacking RV, km
#H - Hardness of target (silo), psi (pound per sq.inch)
OAR*(1-0.5^(8.41*Y^(2/3)/((H^0.7)*((CEP*1.86)^2))))
}
H1 = seq(100,by=100,length.out = 30)
test_data <-
data.frame(
H1,
TK1 = tkprob(OAR=0.95,Y=0.15,H1,CEP=0.09),
TK2 = tkprob(OAR=0.95,Y=0.3,H1,CEP=0.09),
TK3 = tkprob(OAR=0.95,Y=0.45,H1,CEP=0.09),
TK4 = tkprob(OAR=0.95,Y=0.5,H1,CEP=0.09),
TK5 = tkprob(OAR=0.95,Y=1,H1,CEP=0.09)
)
ggplot(test_data, aes(x=H1)) +
geom_line(aes(y = TK1, colour = "0.15"))+
geom_line(aes(y = TK2, colour = "0.3"))+
geom_line(aes(y = TK3, colour = "0.45"))+
geom_line(aes(y = TK4, colour = "0.5"))+
geom_line(aes(y = TK5, colour = "1"))+
ggtitle(label="Kill probability for CEP=0.09 km")+
labs(x="Hardness,psi",y="Probability", colour="Yield,MT",
caption = "based on “Strategic Command and Control.
Redefining the Nuclear Threat” by Bruice G.Blair")
test_data_2 <-
data.frame(
H1,
TK1 = tkprob(OAR=0.95,Y=1,H1,CEP=0.09),
TK2 = tkprob(OAR=0.95,Y=1,H1,CEP=0.15),
TK3 = tkprob(OAR=0.95,Y=1,H1,CEP=0.2),
TK4 = tkprob(OAR=0.95,Y=1,H1,CEP=0.3),
TK5 = tkprob(OAR=0.95,Y=1,H1,CEP=0.5)
)
ggplot(test_data_2, aes(x=H1)) +
geom_line(aes(y = TK1, colour = "0.09"))+
geom_line(aes(y = TK2, colour = "0.15"))+
geom_line(aes(y = TK3, colour = "0.2"))+
geom_line(aes(y = TK4, colour = "0.3"))+
geom_line(aes(y = TK5, colour = "0.5"))+
ggtitle(label="Kill probability for Yield=1MT")+
labs(x="Hardness,psi",y="Probability", colour="CEP,km",
caption = "based on “Strategic Command and Control.
Redefining the Nuclear Threat” by Bruice G.Blair")
#Terminal kill probability of target (silo) by one RV block: Y=1MT,CEP=0.2km,H=1000psi,OAR=0.95
TKP <-tkprob(0.95,1,1000,0.2)
TKP2 <- 1-(1-TKP)^2 #Terminal kill probability of target (silo) by two RV blocks
TKP2## [1] 0.4673365seed=12345
S1 <- round(TKP2*100,0)
S2 <- 100-S1
PK<-rbeta(10000, shape1 = S1, shape2 = S2)
quantile(PK,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.3727763 0.4704085 0.5663720ggplot(data=as.data.frame(PK), aes(PK)) +
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(PK), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(PK,0.025), xmax=quantile(PK,0.975)),
data=as.data.frame(PK), col="#0094EA", size=3) +
ggtitle(label="Probability density of Minituman III silo kill by two RV")
Conclusions
1.Applying Bayesian linear regression model for SLBM data (M,P,D,L,S,W) we can produce posterior distribution sample given Minuteman III data (M=34467 kg,P=1000 kg,D=1.85 m,L=18.2 m,S=3,W=1) with 90% credible interval \(P(4327\le Range\le14653)=0.9\) and the mean \(Range=8674\)
2.Minuteman III still has operational Range corresponding to ICBM which could survive bolt-out-of-the blue attack and hit any possible target world wide.