Multilayer Ballistic Missile Defence
A good deal of information about BMD program one can find here: http://www.mda.mil/system/system.html.
library(expm)## Loading required package: Matrix##
## Attaching package: 'expm'## Следующий объект скрыт от 'package:Matrix':
##
## expmlibrary(diagram)## Loading required package: shapelibrary(ggplot2)
#From Rasmus Baath
#Bayesian First Aid Alternative to the Binomial Test
#http://www.sumsar.net/blog/2014/06/bayesian-first-aid-prop-test/
library(rjags)## Loading required package: coda## Linked to JAGS 4.1.0## Loaded modules: basemod,bugslibrary(mcmc)
library(stringr)## Warning: package 'stringr' was built under R version 3.2.5source(file = "binom_test.R")
source(file = "r_jags.R")
source(file = "generic.R")In previous topic about BMD (http://rpubs.com/alex-lev/42458) we used simplified case for point estimation of overall probability for ICBM to pass through multilayer defense system: boost, midterm, terminal. In order to make situation close to real we will consider the same Markov chain applying Bayesian statistics for evaluating posterior probability densities through binomial test with real data.
To make our case more real we use the following data for ICBM and BMD flight tests: https://en.wikipedia.org/wiki/Intercontinental_ballistic_missile,
https://en.wikipedia.org/wiki/Pukkuksong-1,https://en.wikipedia.org/wiki/Terminal_High_Altitude_Area_Defense, https://en.wikipedia.org/wiki/Aegis_Ballistic_Missile_Defense_System.
icbm_F<-bayes.binom.test(6,29,0.95)# ICBM failed to start successfully
icbm_S<-bayes.binom.test(23,29,0.95)# ICBM started successfully
thaad_F<-bayes.binom.test(5,11,0.8)# THAAD failed to hit ICBM RV
thaad_H<-bayes.binom.test(6,11,0.8)# THAAD hit ICBM RV
aegis_F<-bayes.binom.test(6,37,0.95)# AEGIS failed to hit ICBM
aegis_H<-bayes.binom.test(31,37,0.95)# AEGIS hit ICBM
summary(icbm_F)## Data
## number of successes = 6, number of trials = 29
##
## Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
##
## Measures
## mean sd HDIlo HDIup %<comp %>comp
## theta 0.225 0.073 0.09 0.371 1.000 0.000
## x_pred 6.506 3.046 1.00 12.000 0.006 0.994
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
##
## Quantiles
## q2.5% q25% median q75% q97.5%
## theta 0.099 0.172 0.219 0.271 0.383
## x_pred 1.000 4.000 6.000 8.000 13.000summary(icbm_S)## Data
## number of successes = 23, number of trials = 29
##
## Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
##
## Measures
## mean sd HDIlo HDIup %<comp %>comp
## theta 0.775 0.073 0.629 0.91 1 0
## x_pred 22.494 3.046 16.000 27.00 0 1
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
##
## Quantiles
## q2.5% q25% median q75% q97.5%
## theta 0.617 0.729 0.781 0.828 0.901
## x_pred 16.000 21.000 23.000 25.000 28.000summary(thaad_F)## Data
## number of successes = 5, number of trials = 11
##
## Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
##
## Measures
## mean sd HDIlo HDIup %<comp %>comp
## theta 0.462 0.133 0.212 0.724 0.997 0.003
## x_pred 5.090 2.153 1.000 9.000 0.010 0.990
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.8.
##
## Quantiles
## q2.5% q25% median q75% q97.5%
## theta 0.208 0.367 0.46 0.554 0.721
## x_pred 1.000 4.000 5.00 7.000 9.000summary(thaad_H)## Data
## number of successes = 6, number of trials = 11
##
## Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
##
## Measures
## mean sd HDIlo HDIup %<comp %>comp
## theta 0.538 0.133 0.276 0.788 0.979 0.021
## x_pred 5.910 2.153 1.000 9.000 0.004 0.996
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.8.
##
## Quantiles
## q2.5% q25% median q75% q97.5%
## theta 0.279 0.446 0.54 0.633 0.792
## x_pred 2.000 4.000 6.00 7.000 10.000summary(aegis_F)## Data
## number of successes = 6, number of trials = 37
##
## Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
##
## Measures
## mean sd HDIlo HDIup %<comp %>comp
## theta 0.178 0.060 0.07 0.298 1.000 0.000
## x_pred 6.588 3.158 1.00 12.000 0.006 0.994
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
##
## Quantiles
## q2.5% q25% median q75% q97.5%
## theta 0.078 0.135 0.174 0.216 0.31
## x_pred 1.000 4.000 6.000 9.000 13.00summary(aegis_H)## Data
## number of successes = 31, number of trials = 37
##
## Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
##
## Measures
## mean sd HDIlo HDIup %<comp %>comp
## theta 0.822 0.060 0.702 0.93 0.998 0.002
## x_pred 30.412 3.158 25.000 36.00 0.000 1.000
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
##
## Quantiles
## q2.5% q25% median q75% q97.5%
## theta 0.69 0.784 0.826 0.865 0.922
## x_pred 24.00 28.000 31.000 33.000 36.000icbm_F_post<-as.data.frame(icbm_F$mcmc_samples[1:5000,])
icbm_S_post<-as.data.frame(icbm_S$mcmc_samples[1:5000,])
thaad_F_post<-as.data.frame(thaad_F$mcmc_samples[1:5000,])
thaad_H_post<-as.data.frame(thaad_H$mcmc_samples[1:5000,])
aegis_F_post<-as.data.frame(aegis_F$mcmc_samples[1:5000,])
aegis_H_post<-as.data.frame(aegis_H$mcmc_samples[1:5000,])
pS1S5<-icbm_F_post$theta # ICBM failed to start successfully
pS1S2<-icbm_S_post$theta # ICBM started successfully and went boosting
pS2S3<-aegis_F_post$theta # AEGIS failed to hit ICBM on boost and ICBM went to midterm phase
pS2S5<-aegis_H_post$theta # AEGIS hit ICBM on boost
pS3S4<-pS2S3 # AEGIS failed to hit ICBM on midterm and ICBM went reentry phase
pS3S5<-pS2S5 # AEGIS hit ICBM on midterm
pS4S6<-thaad_F_post$theta # THAAD failed to hit ICBM RV on impact - ICBM passed BMD
pS4S5<-thaad_H_post$theta # THAAD hit ICBM RV on impact - ICBM failed passing BMD
hist(pS1S2,col="lightblue", main="ICBM started successfully")
hist(pS1S5,col="red", main="ICBM failed to start successfully")
hist(pS2S5,col="lightblue", main="AEGIS hit ICBM")
hist(pS2S3,col="red", main="AEGIS failed to hit ICBM")
hist(pS4S5,col="lightblue", main="THAAD hit ICBM RV")
hist(pS4S6,col="red", main="THAAD failed to hit ICBM RV")
Our goal is to calculate the probability distribution of hitting the target by missile carrying one RV or warhead \(P(S1,S6)\). We suppose that the flight of missile is composed by consequent events: \(S1\) - start, \(S2\) - boosting (boost stage), \(S3\) - orbiting (midterm stage), \(S4\) - reentry (terminal), \(S5\) - fail, \(S6\) - impact (hit the target). The events listed and described above (\(S1\), \(S2\), \(S3\), \(S4\), \(S5\), \(S6\)) can be considered as a Markov chain - a mathematical system that undergoes transitions from one state to another on a state space. It is a random process usually characterized as memoryless: the next state depends only on the current state and not on the sequence of events that preceded it. (see http://en.wikipedia.org/wiki/Markov_chain).
\(MC\) is a transition matrix with 6 states (\(S1, S2, S3, S4, S5, S6\)). For instance, for \(S1\) we get two paths to leave: missile would start (be launched) but go down due to technical problems (\(P(S1,S5)\)), missile would start and take a successful boost (\(P(S1,S2)\)). An absorbing state is a state from which it is impossible to leave: \(S5\) (missile technically failed or hit by BMD), \(S6\) (missile hit target). Here we consider RV has 100% reliability.
The most energy-efficient standard trajectory for flying a ballistic missile over a given range carries it high above the atmosphere; the maximum range for a given missile is achieved by traveling on such a standard trajectory passing through all three layers of BMD.
summary(pS1S2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4803 0.7295 0.7790 0.7749 0.8277 0.9743summary(pS1S5)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02569 0.17230 0.22100 0.22510 0.27050 0.51970summary(pS2S3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02055 0.13530 0.17470 0.17850 0.21500 0.45180summary(pS2S5)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5482 0.7850 0.8253 0.8215 0.8647 0.9795summary(pS4S6)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03834 0.36740 0.46250 0.46180 0.55410 0.89400summary(pS4S5)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1060 0.4459 0.5375 0.5382 0.6326 0.9617res1<-c()
res2<-c()
for (i in 1:5000) {
MC<-matrix(c(0.00,pS1S2[i],0,0,pS1S5[i],0,0,0,pS2S3[i],0,pS2S5[i],0,0,0,0,pS3S4[i],pS3S5[i],0,0,0,0,0,pS4S5[i],pS4S6[i],0,0,0,0,1,0,0,0,0,0,0,1),nrow=6,ncol=6,byrow=T)
# MC
MCT<-t(MC)
# MCT
res<-solve(diag(nrow=4,ncol=4)-MC[1:4,1:4])%*%MC[1:4,5:6]
# print(res[1,2])
res2[i]<-res[1,2]
res1[i]<-res[1,1]
}
i## [1] 5000#res2
summary(res2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0001827 0.0059540 0.0102800 0.0126100 0.0164500 0.0984200summary(res1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.9016 0.9836 0.9897 0.9874 0.9940 0.9998hist(res2, col="lightblue", main="Posterior probability density for successfull impact of ICBM RV", xlab = "Probability", freq = F)
hist(res1, col="red", main="Posterior probability density for ICBM failed to pass through BMD", xlab = "Probability", freq = F)
stateNames <- c("S1","S2","S3", "S4", "S5", "S6")
row.names(MCT) <- stateNames; colnames(MC) <- stateNames
plotmat(MCT,pos = NULL,
lwd = 1, box.lwd = 2,
cex.txt = 0.8,
box.size = 0.1,
box.type = "circle",
box.prop = 0.5,
box.col = "light yellow",
arr.length=.1,
arr.width=.1,
self.cex = .4,
self.shifty = -.01,
self.shiftx = .13,
main = "Markov chain for BMD: ICBM case")
So we see very low chances for ICBM to start and pass through all layers of BMD even with 50% of THAAD kill probability. But don’t forget that current AEGIS capabilities are very modest in terms of velocity of potential medium range ICBM target. In case of long range ICBM with special capabilities to pass through BMD AEGIS would hardly hit any target at all.
SLBM
In case of SLBM there is possibility for short flight time with depressed trajectory (low-apogee). Missiles flown on a depressed trajectory (DT) can have significantly shorter flight paths, and therefore significantly shorter flight times, than those flown on a standard trajectory of the same range. For DT there would be effective only two layers of BMD: boost and terminal while midterm phase with 8 km per second velocity of flying SLBM would be unreachable for AEGIS. In that case we have the same Markov chain as above with \(P(S3,S4)=1\) and \(P(S3,S5)=0\). Let’s do it!
res1<-c()
res2<-c()
for (i in 1:5000) {
MC<-matrix(c(0.00,pS1S2[i],0,0,pS1S5[i],0,0,0,pS2S3[i],0,pS2S5[i],0,0,0,0,1,0,0,0,0,0,0,pS4S5[i],pS4S6[i],0,0,0,0,1,0,0,0,0,0,0,1),nrow=6,ncol=6,byrow=T)
# MC
MCT<-t(MC)
# MCT
res<-solve(diag(nrow=4,ncol=4)-MC[1:4,1:4])%*%MC[1:4,5:6]
# print(res[1,2])
res2[i]<-res[1,2]
res1[i]<-res[1,1]
}
i## [1] 5000#res2
summary(res2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.005906 0.042680 0.059370 0.063600 0.079970 0.224900summary(res1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7751 0.9200 0.9406 0.9364 0.9573 0.9941hist(res2, col="lightblue", main="Posterior probability density for successfull impact of SLBM RV", xlab = "Probability", freq = F)
hist(res1, col="red", main="Posterior probability density for SLBM failed to pass through BMD", xlab = "Probability", freq = F)
stateNames <- c("S1","S2","S3", "S4", "S5", "S6")
row.names(MCT) <- stateNames; colnames(MC) <- stateNames
plotmat(MCT,pos = NULL,
lwd = 1, box.lwd = 2,
cex.txt = 0.8,
box.size = 0.1,
box.type = "circle",
box.prop = 0.5,
box.col = "light yellow",
arr.length=.1,
arr.width=.1,
self.cex = .4,
self.shifty = -.01,
self.shiftx = .13,
main = "Markov chain for BMD: SLBM DT case")
For SLBM flying on DT we got just a little better probability to pass through two layers of BMD given midterm is unreachable for AEGIS. All in all there are few chances for successful impact of RV.
Hypersonic weapon system
For hypersonic weapon system (HWS) there would be effective only one layer of BMD if any - terminal while boost and midterm phase would be unreachable for AEGIS. In that case we have the same Markov chain as above with \(P(S2,S3)=1\) and \(P(S2,S5)=0\). Let’s do it!
res1<-c()
res2<-c()
for (i in 1:5000) {
MC<-matrix(c(0.00,pS1S2[i],0,0,pS1S5[i],0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,pS4S5[i],pS4S6[i],0,0,0,0,1,0,0,0,0,0,0,1),nrow=6,ncol=6,byrow=T)
# MC
MCT<-t(MC)
# MCT
res<-solve(diag(nrow=4,ncol=4)-MC[1:4,1:4])%*%MC[1:4,5:6]
# print(res[1,2])
res2[i]<-res[1,2]
res1[i]<-res[1,1]
}
i## [1] 5000#res2
summary(res2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03326 0.28030 0.35440 0.35760 0.43260 0.74190summary(res1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2581 0.5674 0.6456 0.6424 0.7197 0.9667hist(res2, col="lightblue", main="Posterior probability density for successfull impact of HSW", xlab = "Probability", freq = F)
hist(res1, col="red", main="Posterior probability density for HSW failed to pass through BMD", xlab = "Probability", freq = F)
stateNames <- c("S1","S2","S3", "S4", "S5", "S6")
row.names(MCT) <- stateNames; colnames(MC) <- stateNames
plotmat(MCT,pos = NULL,
lwd = 1, box.lwd = 2,
cex.txt = 0.8,
box.size = 0.1,
box.type = "circle",
box.prop = 0.5,
box.col = "light yellow",
arr.length=.1,
arr.width=.1,
self.cex = .4,
self.shifty = -.01,
self.shiftx = .13,
main = "Markov chain for BMD: HSW case")
As we see HWS makes great difference in capability of passing through current BMD (35% mean probability) compared to ICBM (1% mean probability) and SLBM (6% mean probability) cases. The point is that HWS has been tested only few times so far and has far more problems in implementation compared to ICBM and SLBM real systems. That is!
Although this example of Markov Chain is for demonstration purpose only it can be used for real estimates of BMD effectiveness if one should expand state (events) space and would be able to get relevant probabilities proved by test flights (see http://rpubs.com/alex-lev/40511).