From POLARIS to TRIDENT
Alexander Levakov, Senior Research Fellow, Ph.D
August 2021
Introduction
For more then sixty years US nuclear triad has been based on bombers, ICBM and SLBM. The latter - Trident II D5 is of utmost importance for national security as a real prompt global strike asset. Here we estimate evolution of SLBM reliability in terms of test flight program results and operational range.
SLBMs specifications by radar chart
We apply common sense approach for comparison two SLBMs at once with some transformation of data i.e. kg - tons, m - km etc.
Specifications see here:
https://missilethreat.csis.org/missile/trident/#easy-footnote-bottom-5-668
https://web.archive.org/web/20170123163450/http://www.designation-systems.net/dusrm/m-27.html
library(fmsb)
A3.spec <- data.frame(L=9.86,D=1.37,S=2,WH=3,
M=16.200,P=0.760,R=4.500,CEP=0.6)
D5.spec <- data.frame(L=13.6,D=2.11,S=3,WH=12,
M=59.000,P=2.800,R=12.000,CEP=0.090)
slbm.spec <- rbind(A3.spec,D5.spec)
rownames(slbm.spec) <- c("Polaris A3","Trident II D5")
#We have to add 2 lines to the dataframe: the max and min of each variable to show on the plot!
slbm.spec.2 <- rbind(c(14,3,3,14,60,3.0,13.0,0.6),c(9,1.0,2,1,16.0,0.7,4.0,0.09),slbm.spec)
# Color vector
colors_border=c( rgb(0.2,0.5,0.5,0.9), rgb(0.8,0.2,0.5,0.9) , rgb(0.7,0.5,0.1,0.9) )
colors_in=c( rgb(0.2,0.5,0.5,0.4), rgb(0.8,0.2,0.5,0.4) , rgb(0.7,0.5,0.1,0.4) )
# plot with default options:
radarchart( slbm.spec.2 , axistype=1 ,
#custom polygon
pcol=colors_border , pfcol=colors_in , plwd=4 , plty=1,
#custom the grid
cglcol="grey", cglty=1, axislabcol="grey", caxislabels=seq(0,20,5), cglwd=0.8,
#custom labels
vlcex=0.8
)
# Add a legend
legend(x=0.7, y=1, legend = rownames(slbm.spec.2[-c(1,2),]), bty = "n", pch=20 , col=colors_in , text.col = "grey", cex=1.2, pt.cex=3)
Reliability
Polaris A3
Polaris A3 SLBM
See https://rpubs.com/alex-lev/774763
library(ggplot2)
seed=12345
PD3<-rbeta(10000, shape1 = 20, shape2 = 5)
bnt3 <- binom.test(15,20,0.9)
print(bnt3)##
## Exact binomial test
##
## data: 15 and 20
## number of successes = 15, number of trials = 20, p-value = 0.04317
## alternative hypothesis: true probability of success is not equal to 0.9
## 95 percent confidence interval:
## 0.5089541 0.9134285
## sample estimates:
## probability of success
## 0.75quantile(PD3,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.6244974 0.8071248 0.9283659ggplot(data=as.data.frame(PD3), aes(PD3)) +
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(PD3), color = "red") +
geom_errorbarh(aes(y=-5, x=mean(PD3), xmin=quantile(PD3,0.025),
xmax=quantile(PD3,0.975)),
data=as.data.frame(PD3), col="#0094EA", size=3) +
ggtitle(label="Probability density of Polaris A3 successful trial after test flight program") +
labs(x="Probability",y="Density")
Trident II D5
Trident II D5 SLBM
See https://rpubs.com/alex-lev/377492
seed=12345
PD5<-rbeta(10000, shape1 = 23, shape2 = 4)
bnt5 <- binom.test(23,27,0.9)
print(bnt5)##
## Exact binomial test
##
## data: 23 and 27
## number of successes = 23, number of trials = 27, p-value = 0.3403
## alternative hypothesis: true probability of success is not equal to 0.9
## 95 percent confidence interval:
## 0.6626891 0.9581126
## sample estimates:
## probability of success
## 0.8518519quantile(PD5,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.7019555 0.8603876 0.9553582ggplot(data=as.data.frame(PD5), aes(PD5)) +
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(PD5), color = "red") +
geom_errorbarh(aes(y=-5, x=mean(PD5), xmin=quantile(PD5,0.025), xmax=quantile(PD5,0.975)),
data=as.data.frame(PD5), col="#0094EA", size=3) +
ggtitle(label="Probability density of Trident II D5 successful trial after test flight program")+
labs(x="Probability",y="Density")
Comparison of initial reliability
quantile(PD5,probs = c(0.025,0.5,0.975)) - quantile(PD3,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.07745811 0.05326275 0.02699235A3 <- data.frame(prob=PD3,name='Polaris A3')
D5 <- data.frame(prob=PD5,name='Trident II D5')
slbm <- rbind(A3,D5)
ggplot(slbm, aes(prob, fill = name)) + geom_density(alpha = 0.2)+
ggtitle(label="Evolution of US SLBM Initial Reliability")+
labs(x="Probability",y="Density",subtitle = "Test Flight Program Results") + scale_fill_discrete(name="SLBM")
Operational Range
SLBM data
load("slbm.dat")
library(DT)
#library(knitr)
datatable(slbm)#kable(slbm,format="latex", caption = "SLBM data")Correlations
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
heatmap(cor(slbm))
Bayesian Linear Regression Model
library(rstan)
library(rstanarm)
library(bayesplot)
options(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) -41.1 107.2
## S 9.4 5.2
## D 37.3 19.6
## L 0.8 1.7
## W 3.3 5.6
## log(M) 6.7 15.0
## log(P) -7.3 6.0
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 10.2 1.7
##
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanregplot(fit.slbm.bs)
posterior_vs_prior(fit.slbm.bs)
fit.slbm.bs2<-as.array(fit.slbm.bs)
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
## sample: 20000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 26
## predictors: 7
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) -41.4 110.0 -178.1 -41.1 96.4
## S 9.4 5.3 2.7 9.4 16.1
## D 37.2 20.5 11.7 37.3 62.6
## L 0.8 1.8 -1.4 0.8 3.1
## W 3.2 5.8 -4.1 3.3 10.5
## log(M) 6.7 15.8 -13.0 6.7 26.2
## log(P) -7.3 6.2 -15.2 -7.3 0.5
## sigma 10.4 1.8 8.3 10.2 12.8
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 69.0 2.9 65.3 69.0 72.7
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 1.1 1.0 9162
## S 0.0 1.0 13466
## D 0.2 1.0 9822
## L 0.0 1.0 12332
## W 0.1 1.0 13064
## log(M) 0.2 1.0 8766
## log(P) 0.1 1.0 15290
## sigma 0.0 1.0 11985
## mean_PPD 0.0 1.0 19056
## log-posterior 0.0 1.0 6342
##
## 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).Polaris A3
A3.slbm.pp <- posterior_predict(fit.slbm.bs,
newdata = data.frame(L=9.86,D=1.37,S=2,W=2,
M=16200,P=760),seed=12345)
#Note. Input parameters for the model
#http://www.astronautix.com/p/polarisa3.html
A3.Range.km <- A3.slbm.pp^2
A3.Range.km <- A3.Range.km[1:10000,]
quantile(A3.Range.km,probs = c(0.05,0.5,0.95))## 5% 50% 95%
## 1553.025 3478.665 6211.902library(grid)
ggplot(data=as.data.frame(A3.Range.km), aes(A3.Range.km)) +
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(A3.Range.km), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(A3.Range.km,0.05),
xmax=quantile(A3.Range.km,0.95)),
data=as.data.frame(A3.Range.km), col="#0094EA", size=3) +
ggtitle(label="Probability density of Polaris A3 range")+labs(y="Density", x="Range km")
Trident II D5
D5.slbm.pp <- posterior_predict(fit.slbm.bs,
newdata = data.frame(L=13.6,D=2.11,S=3,W=2,
M=59000,P=2800),seed=12345)
#Note. Payload is an expert estimate for the model
#https://missilethreat.csis.org/missile/trident/
D5.Range.km <- D5.slbm.pp^2
D5.Range.km <- D5.Range.km[1:10000,]
quantile(D5.Range.km,probs = c(0.05,0.5,0.95))## 5% 50% 95%
## 6289.607 9632.449 13627.746ggplot(data=as.data.frame(D5.Range.km), aes(D5.Range.km)) +
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(D5.Range.km), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(D5.Range.km,0.05),
xmax=quantile(D5.Range.km,0.95)),
data=as.data.frame(D5.Range.km), col="#0094EA", size=3) +
ggtitle(label="Probability density of Trident II D5 range")+labs(y="Density",x="Range km")
Comparison of operational range
A3.R <- data.frame(Range=A3.Range.km,name='Polaris A3')
D5.R <- data.frame(Range=D5.Range.km,name='Trident II D5')
slbm.R <- rbind(A3.R,D5.R)
ggplot(slbm.R, aes(Range, fill = name)) + geom_density(alpha = 0.2)+
ggtitle(label="Evolution of US SLBM Range")+
labs(x="Range km",y="Density",) + scale_fill_discrete(name="SLBM")
Terminal kill probability
Here we apply formula from Bruice G.Blair “Strategic Command and Control. Redefining the Nuclear Threat”, page 305.
\(P_{kill}=OAR(1-0.5^{8.41Y^{2/3}/H^{2/3}CEP^2})\)
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)
return(OAR*(1-0.5^(8.41*Y^(2/3)/((H^(2/3))*((CEP*1.86)^2)))))
}
H1 = seq(100,by=100,length.out = 30)
test_data <-
data.frame(
H1,
TK0 = tkprob(OAR=0.95,Y=0.2,H1,CEP=0.6),
TK1 = tkprob(OAR=0.95,Y=0.15,H1,CEP=0.09)
)
ggplot(test_data, aes(x=H1)) +
geom_line(aes(y = TK0, colour = "Polaris A3"))+
geom_line(aes(y = TK1, colour = "Trident II D5"))+
ggtitle(label="Kill probability by single warhead")+
labs(x="Hardness of target, psi",y="Probability", colour="SLBM",
caption = "based on “Strategic Command and Control.
Redefining the Nuclear Threat” by Bruice G.Blair")
ggplot(test_data, aes(x=H1)) +
geom_line(aes(y = 1-(1-TK0)^3, colour = "Polaris A3"))+
geom_line(aes(y = 1-(1-TK1)^3, colour = "Trident II D5"))+
ggtitle(label="Kill probability by three warheads")+
labs(x="Hardness of target, psi",y="Probability", colour="SLBM",
caption = "based on “Strategic Command and Control.
Redefining the Nuclear Threat” by Bruice G.Blair")