Demystifying TITAN-II
Alexander Levakov, Senior Research Fellow, Ph.D
March, 2019
TITAN-II
About ICBM
The Titan II was the largest and heaviest missile ever built by the United States. The missile was 31.3 m long and 3.05 m wide. It weighed 149,700 kilograms when fully fueled and had a range of 15,000 km. Its inertial guidance system gave an accuracy of 900 meters CEP and was capable of making in-flight corrections without ground control input. The missile was upgraded with an improved guidance system in 1979. The missile used a single Mk 6 Reentry Vehicle (RV) which carried a W-53 9.0 MT nuclear warhead. The missile had a diameter of 3.05 m, a length of 31.30 m and a launch weight of 149,700 kg. The missiles had a two-stage liquid propellant design and reached a speed of 25 times the speed of sound by the time the engines cut off.
For more see:
1.https://missilethreat.csis.org/missile/titan-ii/
2.https://en.wikipedia.org/wiki/LGM-25C_Titan_II
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)
library(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=31,D=3,S=2,W=1,
M=149700,P=3700),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%
## 9504.752 18711.820 30817.618mean(Range.km>15000)## [1] 0.689ggplot(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 TITAN-II Range")
ICBM Flight Test Reliability
Here we use open source data on TITAN-II ICBM flight test events (77 success and 5 failure) and Bayesian Binomial model.
library(rstan)
x <- 77
n <- 77+5
data_list <- list(n=n,x=x,a=2,b=1)
set.seed(12345)
fit.bin <- stan(file = "binom.stan",data = data_list,chains = 2,
iter = 5000,warmup = 1000)stancode <- rstan::get_stancode(fit.bin)
cat(stancode)## data{
## int n; // Number of trials
## int x; // Number of success
## int a; // Alpha
## int b; // Beta
## }
##
## parameters{
## real<lower=0,upper=1> p; // unknown parameter
## }
##
## model{
## //p~uniform(0,1);
## p~beta(a+1,b+1);
## x~binomial(n,p);
## //x ~ beta_binomial(n, a+1, b+1);
##
##
## }
##
## generated quantities {
## real pp;
## pp = beta_binomial_rng(n, a+1, b+1);
## }fit.bin## Inference for Stan model: binom.
## 2 chains, each with iter=5000; warmup=1000; thin=1;
## post-warmup draws per chain=4000, total post-warmup draws=8000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## p 0.92 0.00 0.03 0.86 0.90 0.92 0.94 0.97 2487 1
## pp 49.17 0.19 16.76 15.00 37.00 50.00 62.00 77.00 8000 1
## lp__ -24.86 0.01 0.71 -26.86 -25.01 -24.59 -24.40 -24.35 3075 1
##
## Samples were drawn using NUTS(diag_e) at Sun Mar 17 16:00:26 2019.
## For each parameter, 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).stan_hist(fit.bin,pars = "p")
stan_trace(fit.bin)
ICBM Warhead Kill Probability
Here we apply formula from Bruice G.Blair “Strategic Command and Control. Redefining the Nuclear Threat”, page 305. We want to estimate kill probability for TITAN-II warhead given CEP=900 m and Y=9 Mt.
\(P_{kill}=OAR(1-0.5^{8.41Y^{2/3}/H^{2/3}CEP^2})\)
Single Shot
library(ggplot2)
#Single shot kill probability
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(10,by=10,length.out = 30)
Y1 <- 9
OAR1 <- 0.97
test_data <-
data.frame(
H1,
TK1 = tkprob(OAR1,Y1,H1,CEP=0.9)
)
ggplot(test_data, aes(x=H1)) +
geom_line(aes(y = TK1, colour = "9"))+
ggtitle(label="Single Shot Kill Probability for CEP=0.9 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")
Two Shots
seed=12345
PT2<-rbeta(10000, shape1 = 77, shape2 = 5)
bnt <- binom.test(77,82,0.95)
print(bnt)##
## Exact binomial test
##
## data: 77 and 82
## number of successes = 77, number of trials = 82, p-value = 0.6075
## alternative hypothesis: true probability of success is not equal to 0.95
## 95 percent confidence interval:
## 0.8634132 0.9799070
## sample estimates:
## probability of success
## 0.9390244bnt.low <- bnt$conf.int[1]
bnt.up <- bnt$conf.int[2]
quantile(PT2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.8803911 0.9428891 0.9792855ggplot(data=as.data.frame(PT2), aes(PT2)) +
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(PT2), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(PT2,0.025), xmax=quantile(PT2,0.975)),
data=as.data.frame(PT2), col="#0094EA", size=3) +
ggtitle(label="Probability density of TITAN-II successful trial after 82 test flights")
seed=12345
TKP <-tkprob(PT2[1:5000],rnorm(5000,9,0.05),30,rnorm(5000,0.9,0.05))
TKP2 <- 1-(1-TKP)^2 #Terminal kill probability of target by two warheads
quantile(TKP2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.7427380 0.8149041 0.8813279ggplot(data=as.data.frame(TKP2), aes(TKP2)) +
geom_histogram(bins = 20,col="black",fill="green") +
geom_vline(xintercept = mean(TKP2), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(TKP2,0.025), xmax=quantile(TKP2,0.975)),
data=as.data.frame(TKP2), col="#0094EA", size=3) +
ggtitle(label="Probability density of target (30 psi) kill by two warheads")+labs(x="Pkill")
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
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))
}
1-(1-SSKP(9,1000,900))^2## [1] 0.8055555Conclusions
1.Applying Bayesian linear regression model for SLBM data (M,P,D,L,S,W) we can produce posterior distribution sample given TITAN-II data (M=149700 kg,P=3700 kg,D=3 m,L=30 m,S=2,W=1) with 80% credible interval \(P(9505\le Range\le30818)=0.8\) and the mean \(Range=18712\).
2.Applying Bayesian binomial model for ICBM historical flight test events record (x=77, n=82) we can get 95% credible interval for TITAN-II reliability as \(P(0.86\le Reliability\le0.97)=0.95\) and the mean \(M[Reliability]=0.92\).
3.TITAN-II had operational Range and Reliability corresponding to ICBM which could hit any medium hardened target (\(H\le30psi\)) worldwide \(P(0.743\le P_{kill}\le0.881)=0.95\).