DF41 survivability under Trident II D5 attack
Alexander Levakov, Senior Research Fellow, Ph.D
February, 2020
DF41 ICBM
DF41
China - Heilongjiang
The American website Washington Free Beacon published an article on 05 June 2018 saying that China recently completed the 10th flight test of the DF-41 intercontinental ballistic missile (ICBM), marking a substantial step forward for the actual deployment of China’s latest and most powerful strategic weapon. An official from the U.S. Department of Defense said that the DF-41 missile was launched from Taiyuan City of northern China’s Shanxi Province on May 27, 2018, and travelled thousands of miles before it hit a simulated target in western China’s Gobi desert.
For more see: https://www.globalsecurity.org/wmd/world/china/df-41.htm
Some Hong Kong and Taiwan media reported that pictures of China’s Dongfeng-41 ballistic missile were exposed on Chinese mainland websites. It was revealed that the pictures were taken in Heilongjiang Province. Military analysts believe that this is perhaps the second Dongfeng-41 strategic missile brigade and it should be deployed in northeastern China.
According to reports, the Dongfeng-41 is a nuclear solid-fuel road-mobile intercontinental ballistic missile. With a range of 14,000 kilometers and a payload of 10-12 nuclear warheads, it can target anywhere in the world and is widely considered one of the most advanced intercontinental ballistic missiles.
For more see: http://www.globaltimes.cn/content/1030353.shtml
SLBM data
load("slbm.dat")
library(DT)
datatable(slbm)library(ggplot2)
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2#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(rstanarm)## Loading required package: Rcpp## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## rstanarm (Version 2.18.2, packaged: 2018-11-08 22:19:38 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.7.0## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme settingoptions(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
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 10.2 1.7
##
## Sample avg. posterior predictive distribution of y:
## Median MAD_SD
## mean_PPD 69.0 2.8
##
## ------
## * 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)##
## 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 8838
## S 0.0 1.0 12972
## D 0.2 1.0 9062
## L 0.0 1.0 12215
## W 0.1 1.0 11906
## log(M) 0.2 1.0 8349
## log(P) 0.1 1.0 13073
## sigma 0.0 1.0 10097
## mean_PPD 0.0 1.0 19001
## log-posterior 0.0 1.0 6237
##
## 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).Trident II D5 Operational Range
Trident II D5 SLBM
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/
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%
## 6238.101 9636.941 13706.396library(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 Trident II D5 range",
subtitle = "P(6238<Range.km<13706)=0.9")
Current Reliability of Trident II D5
Now we want to estimate Trident II D5 SLBM reliability, assuming 177 successful and 6 failed historical DASO trials.
seed=12345
PD5<-rbeta(10000, shape1 = 177, shape2 = 6)
bnt <- binom.test(171,177,0.95)
print(bnt)##
## Exact binomial test
##
## data: 171 and 177
## number of successes = 171, number of trials = 177, p-value = 0.3913
## alternative hypothesis: true probability of success is not equal to 0.95
## 95 percent confidence interval:
## 0.9276809 0.9874602
## sample estimates:
## probability of success
## 0.9661017bnt.low <- bnt$conf.int[1]
bnt.up <- bnt$conf.int[2]
quantile(PD5,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.9364940 0.9686164 0.9878033ggplot(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, 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 DASO 30")
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.006,H1,CEP=0.09),
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 = TK0, colour = "0.006"))+
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")
Hardened silo based DF41
DF41 test silo
#Terminal kill probability of target (silo) by one RV block:
#M[Y}=0.4MT & SD[Y]=0.01MT,M[CEP]=0.09km & SD[CEP]=0.01km,H=1000psi,OAR=PD5
seed=12345
TKP <-tkprob(PD5[1:5000],rnorm(5000,0.45,0.01),1000,rnorm(5000,0.09,0.01))
TKP2 <- 1-(1-TKP)^2 #Terminal kill probability of target (silo) by two RV blocks
quantile(TKP2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.7902819 0.8972163 0.9726899ggplot(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 DF41 silo (1000 psi) kill by two D5 RV 0.45MT")
Mobile fixed DF41 without tactical warning
https://www.princeton.edu/sgs/faculty-staff/bruce-blair/One-Hundred-Nuclear-Wars-B-Blair.pdf
“One Hundred Nuclear Wars: Stable Deterrence between the United States and Russia at Reduced Nuclear Force Levels Off Alert in the Presence of Limited Missile Defenses, Victor Esin, Matthew Mckinzie, Valery Yarynich, Pavel Zolotarev,”Science & Global Security 19, no. 3 (2011): 167-194.
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
The motion of an object moving near the surface of the earth can be described using the equations: \(x = x_o + vx_ot\) (1)
\(y = y_o + vy_ot - 0.5gt^2\) (2)
We should solve these two simultaneous equations to obtain a description of the ballistic trajectory. The horizontal \(x\) and vertical \(y\) components of initial velocity are determined from:
\(v_{xo} = v_ocos(\theta)\)
\(v_{yo} = v_osin(\theta)\)
Then we compute the time to reach the maximum height \(h\) by finding the time \(t_{rise}\) at which vertical velocity \(v_{y}\) becomes zero:
\(v_y = v_{yo} - gt\)
\(t_{rise} = v_{yo}/g\)
Maximum height is obtained by substitution of this time into equation (2).
\(h_{max} = y_o + v_{yo}t - 0.5gt^2\)
Next, the time to fall from the maximum height is computed by solving equation (2) for an object dropped from the maximum height with zero initial velocity.
\(h - 0.5gt^2=0\)
\(t_{fall} = \sqrt{2h/g}\)
The total flight time of the projectile is then:
\(t_{flight} = t_{rise} + t_{fall}\)
From this, equation (1) gives the maximum range:
\(R = v_{xo}t_{flight}\)
The projectile speed at impact \(v_f\) is determined by applying the Pythagorean Theorem:
\(v_f = \sqrt{(v_{xf}^2 + v_{yf}^2)}\)
In which:
\(v_{xf} = v_{xo}\)
\(v_{yf} = -gt_{fall}\)
The question remains what is the surface distance L between the launch and impact point in this situation. We know from the conservation of angular momentum that \(v_{\theta}=RV_0cos(\beta)/r\) and for \(H/R<<1\) that \(r\) is well approximated by \(R\). Thus we have the range \(L=tV_0 cos(\beta)\) and \(t=L/(V_0cos(\beta))\). For \(V_0=5.5 km/sec\), \(\beta=30\) and \(L=2500 km\) we got \(t=2500/(5.5*\sqrt3/2)=525 sec = 525/60=8.75 min=9min\)
#Terminal kill probability of fixed land soft target by one RV block:
#M[Y}=0.1MT & SD[Y]=0.01MT,M[CEP]=0.09km & SD[CEP]=0.01km,H=30psi,OAR=PD5
#L=21m,R=2.25m S=2*PI*R*L=6.28*2.25*21=297sq.m=2970000sq.sm F=H*0.43/(0.25*0.25)*S/2=30*0.07*148.5*10000=3118500kg
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))
}
#### 0.1 MT
seed=12345
TKP <-tkprob(PD5[1:5000],rnorm(5000,0.1,0.01),30,rnorm(5000,0.09,0.01))
TKP2 <- 1-(1-TKP)^2 #Terminal kill probability of target (silo) by two RV blocks
quantile(TKP2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.9920810 0.9981175 0.9996784ggplot(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 fixed mobile DF41 (30 psi) kill by two D5 RV 0.1MT")+labs(x="Pkill") 
#### 6 KT
seed=12345
TKP <-tkprob(PD5[1:5000],rnorm(5000,0.006,0.0005),30,rnorm(5000,0.09,0.01))
TKP2 <- 1-(1-TKP)^2 #Terminal kill probability of target (silo) by two RV blocks
quantile(TKP2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.5985115 0.7405592 0.8895959ggplot(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 fixed mobile DF41 (30 psi) kill by two D5 RV 6 KT")+labs(x="Pkill") 
Moving mobile DF41 with tactical warning
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 im- portant factors in determining the extent of damage at the surface.
#Terminal kill probability of fixed land soft target by one RV block:
#M[Y}=0.1MT & SD[Y]=0.01MT,M[CEP]=0.09km & SD[CEP]=0.01km,H=30psi,OAR=PD5
#VDF41=60km/h, TW=2 min, LD5=2500 km, VD5=5 km/s, TD5=LD5/VD5=500 sec (8.3 min),
#LDF41=VDF41*(TD5-TW)=60*(9-2)/60=7 km, SDF41=3.14*LDF41*LDF41=113 sq.km
# D=D1*Y^1/3 D=600[1kt,30psi]*0.3*100^1/3=6000m=6km[100kt,30psi]
#(D/LDF41)^2=(6/7)^2=0.735
seed=12345
PKSQ <- function(Y,H,V,TT){
#Lethal area versus moving target area
(LRM(Y,H)*0.001)^2/(V*TT/60)^2
}
#With tactical warning DF41 starts moving with V=30 km/h after 2 minutes having got 9-2=7 minutes before D5 RV impact.
TKPA <-PD5[1:3000]*PKSQ(1,15,rnorm(3000,30,1),rnorm(3000,7,0.1))
TKPA2 <- 1-(1-TKPA)^2 #Terminal kill probability of target (silo) by two RV blocks
quantile(TKPA2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.6485256 0.7170054 0.7912849ggplot(data=as.data.frame(TKPA2), aes(TKPA2)) +
geom_histogram(bins = 20,col="black",fill="green") +
geom_vline(xintercept = mean(TKPA2), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(TKPA2,0.025), xmax=quantile(TKPA2,0.975)),
data=as.data.frame(TKPA2), col="#0094EA", size=3) +
ggtitle(label="Probability density of moving mobile DF41 kill by two D5 RV 1MT")+labs(x="Pkill", caption = "Hardness 15 psi, V=30 km/h, Distance=2550 km") 
#With tactical warning DF41 starts moving with V=40 km/h after 2 minutes having got 9-2=7 minutes before D5 RV impact.
TKPA <-PD5[1:3000]*PKSQ(1,15,rnorm(3000,40,1),rnorm(3000,7,0.1))
TKPA2 <- 1-(1-TKPA)^2 #Terminal kill probability of target (silo) by two RV blocks
quantile(TKPA2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.4150260 0.4568069 0.5039178ggplot(data=as.data.frame(TKPA2), aes(TKPA2)) +
geom_histogram(bins = 20,col="black",fill="green") +
geom_vline(xintercept = mean(TKPA2), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(TKPA2,0.025), xmax=quantile(TKPA2,0.975)),
data=as.data.frame(TKPA2), col="#0094EA", size=3) +
ggtitle(label="Probability density of moving mobile DF41 kill by two D5 RV 1MT")+labs(x="Pkill", caption = "Hardness - 15 psi, V=40 km/h, Distance=2500 km") 
#With tactical warning DF41 starts moving with V=50 km/min after 2 minutes having got 9-2=7 minutes before D5 RV impact.
TKPA <-PD5[1:3000]*PKSQ(1,15,rnorm(3000,50,1),rnorm(3000,7,0.1))
TKPA2 <- 1-(1-TKPA)^2 #Terminal kill probability of target (silo) by two RV blocks
quantile(TKPA2,probs = c(0.025,0.5,0.975))## 2.5% 50% 97.5%
## 0.2811364 0.3087687 0.3388762ggplot(data=as.data.frame(TKPA2), aes(TKPA2)) +
geom_histogram(bins = 20,col="black",fill="green") +
geom_vline(xintercept = mean(TKPA2), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(TKPA2,0.025), xmax=quantile(TKPA2,0.975)),
data=as.data.frame(TKPA2), col="#0094EA", size=3) +
ggtitle(label="Probability density of moving mobile DF41 kill by two D5 RV 1MT")+labs(x="Pkill", caption = "Hardness - 15 psi, V=50 km/h, Distance=2500 km") 
Conclusions
1.DF41 ICBM based in fixed silo (H=1000 psi) would be destroyed by two Trident D5 RV 0.45 MT each with the 95% credible interval \(P(0.790\leq P_{kill} \leq0.973)=0.95\) with the mean value \(M(P_{kill})=0.897\).
2.Ground mobile DF41 in fixed position (H=30 psi) without tactical warning would be destroyed by two Trident D5 RV 0.1 MT each with the 95% credible interval \(P(0.992\leq P_{kill} \leq0.999)=0.95\) with the mean value \(M(P_{kill})=0.998\).
3.Mobile DF41 (H=15 psi) with tactical warning (2 minutes) would be destroyed by two Trident D5 RV 1 MT each from the distance 2500 km after 7 minutes run (V=30 km/h) with the 95% credible interval \(P(0.651\leq P_{kill} \leq0.793)=0.95\) with the mean value \(M(P_{kill})=0.718\).
4.Mobile DF41 (H=15 psi) with tactical warning (2 minutes) would be destroyed by two Trident D5 RV 1 MT each from the distance 2500 km after 7 minutes run (V=40 km/h) with the 95% credible interval \(P(0.415\leq P_{kill} \leq0.506)=0.95\) with the mean value \(M(P_{kill})=0.457\).
5.Mobile DF41 (H=15 psi) with tactical warning (2 minutes) would be destroyed by two Trident D5 RV 1 MT each from the distance 2500 km after 7 minutes run (V=50 km/h) with the 95% credible interval \(P(0.282\leq P_{kill} \leq0.338)=0.95\) with the mean value \(M(P_{kill})=0.308\).
6.US could not be able to destroy DF41 ICBM both fixed and mobile without collateral damage to the neighboring country (Russia).