install.packages("copula")
install.packages("plot3D")
library(copula)
library(plot3D)



#====== a1:=Alpha_1, a1:=Alpha_2,  b1:=Beta1_1,   b2:=Beta1_2    
#====== S1:=omega_1, S2:=omega_2, p:=lambda   ,   t0:=t
#====== m1:=M_1    , m2:=M_2  
#====== x1 and x2 are the initial states of components in each element. 
 
  x1=0
  x2=0
  t0=4
  p=1
  
  a1=2
  b1=1
  a2=3
  b2=1
  S1=4
  S2=7
  c1=5
  c2=5
  cr=0.05
  cp=0.01
  cu=10
  
  
  Fy1=function(y1,t){
    pgamma(y1,(a1*t),1/b1)
  }
  
  fy1=function(y1,t){
    dgamma(y1,(a1*t),1/b1)
  }
  
  
  Fy2=function(y2,t){
    pgamma(y2,(a2*t),1/b2)
  }
  
  fy2=function(y2,t){
    dgamma(y2,(a2*t),1/b2)
  }
  
  
  Fc=function(y1,y2,t){
    pCopula(c(Fy1(y1,t),Fy2(y2,t)),fgmCopula(p))
  }
  
  
  
  fc=function(y1,y2,t){
    if (Fy1(y1,t)<1*10^(-4) | Fy2(y2,t)<1*10^(-4)) 
      0  else
        dCopula(c(Fy1(y1,t),Fy2(y2,t)),fgmCopula(p))*fy1(y1,t)*fy2(y2,t)
  }
  
  
  Fbc=function(y1,y2,t){
    Fc(y1,y2,t)-Fy2(y2,t)+1-Fy1(y1,t)
  }
  
  #========================  Cost function finite time   =========================
  
  I0=function(y1,t){
    if (y1==0) 1
    else 0
  }
  
  Iy1=function(y1,t,m1){
    if (y1<=m1 & x1<=y1) 1
    else 0
  }
  
  Q11=function(y1,t,m1){
    (1-Fy1(m1-x1,t))*I0(y1,t)+fy1(y1-x1,t)*Iy1(y1,t,m1) 
  }
  
  
  
  
  
  IIy2=function(y2,t,m2){
    if (y2<=m2 & x2<=y2) 1
    else 0
  }
  
  
  Q22=function(y2,t,m2){
    (1-Fy2(m2-x2,t))*I0(y2,t)+fy2(y2-x2,t)*IIy2(y2,t,m2) 
  }
  
  
  
  #===   Density function  ===========================
  
  
  J1=function(y1,y2,t,m1,m2){
    if (y1<=m1 & x1<=y1 & y2<=m2 & x2<=y2) 1
    else 0
  }
  
  
  Q1=function(y1,y2,t,m1,m2){
    J1(y1,y2,t,m1,m2)*fc(y1-x1,y2-x2,t)
  }
  
  
  J2=function(y1,y2,t,m2){
    if (y2<=m2 & x2<=y2 & y1==0) 1
    else 0
  }
  
  
  D=function(y2,t,m1){
    dd=function(u) {fc(u-x1,y2-x2,t) }
    dd=Vectorize(dd)
    integrate(dd,m1,30)$value
  }
  
  
  Q2=function(y1,y2,t,m1,m2){
    J2(y1,y2,t,m2)*D(y2,t,m1)
  }
  
  J3=function(y1,y2,t,m1){
    if (y1<=m1 & x1<=y1 & y2==0) 1
    else 0
  }
  
  
  Dq=function(y1,t,m2){
    dd=function(u) {fc(y1-x1,u-x2,t) }
    dd=Vectorize(dd)
    integrate(dd,m2,30)$value
  }
  
  
  Q3=function(y1,y2,t,m1,m2){
    J3(y1,y2,t,m1)*Dq(y1,t,m2)
  }
  
  
  
  J4=function(y1,y2,t){
    if (y1==0 & y2==0) 1
    else 0
  }
  
  Q4=function(y1,y2,t,m1,m2){
    Fbc(m1-x1,m2-x2,t)*J4(y1,y2,t)  
  }
  
  Q=function(y1,y2,t,m1,m2){
    Q1(y1,y2,t,m1,m2)+Q2(y1,y2,t,m1,m2)+Q3(y1,y2,t,m1,m2)+Q4(y1,y2,t,m1,m2)
  }
 
  
 
  #===============================   n>=2   ======================================
  
  
  f1=function(T,m1) {
    a=c()
    a[1]=1-Fy1(m1-x1,T) 
    if (floor(t0/T)>=2){
      for(i in 2:floor(t0/T))
        a[i]=integrate(Vectorize(function(y1) a[i-1]*Q11(y1,T,m1)),0.0001,m1)$v}
    return(a)
  }
  
  
  f2=function(T,m2) {
    b=c()
    b[1]=1-Fy2(m2-x2,T) 
    if (floor(t0/T)>=2){
      for(i in 2:floor(t0/T))
        b[i]=integrate(Vectorize(function(y2) b[i-1]*Q22(y2,T,m2)),0.0001,m2)$v}
    return(b)
  }
  
  
  
  f12=function(T,m1,m2) {
    a=c()
    a[1]=Fbc(m1-x1,m2-x2,T)
    if (floor(t0/T)>=2){
      for(i in 2:floor(t0/T))
        a[i]=integrate(Vectorize(function(y2) 
          integrate(Vectorize(function(y1) a[i-1]*Q(y1,y2,T,m1,m2)),0.005,m1)$v),0.005,m2)$v}
    return(a)
  }
  
  
  
  
  ff1=function(T,m1){
    n=floor(t0/T)
    return(sum(f1(T,m1)[1:n]))
  }
  
  
  ff2=function(T,m2){
    n=floor(t0/T)
    return(sum(f2(T,m2)[1:n]))
  }
  
  
  
  
  ff12=function(T,m1,m2){
    n=floor(t0/T)
    return(sum(f12(T,m1,m2)[1:n]))
  }
  
  #=======================     Cost function     ==============================
  
  
  UT1=function(T,x11,x12){
    LL=function(u) {Fbc(S1-x11,S2-x12,u)}
    LL=Vectorize(LL)
    integrate(LL,0.0001,T)$value
  }
  
  UT2=function(T,x21,x22){
    LL=function(u) {Fbc(S1-x21,S2-x22,u)}
    LL=Vectorize(LL)
    integrate(LL,0.0001,T)$value
  }
  
  
  
  
  
  UU0=function(y11,T,m1,m2){
    LL=function(y12) {UT1(t0-T,y11,y12)*Q(y11,y12,T,m1,m2)}
    LL=Vectorize(LL)
    (integrate(LL,0.0001,m2)$value)
  }
  
  
  UU1=function(T,m1,m2){
    LL=function(y11) {UU0(y11,T,m1,m2)}
    LL=Vectorize(LL)
    integrate(LL,0.0001,m1)$value
  }
  

  
  UU00=function(y21,T,m1,m2){
    LL=function(y22) {UT2(t0-T,y21,y22)*Q(y21,y22,T,m1,m2)}
    LL=Vectorize(LL)
    (integrate(LL,0.0001,m2)$value)
  }
  
  
  UU2=function(T,m1,m2){
    LL=function(y21) {UU00(y21,T,m1,m2)}
    LL=Vectorize(LL)
    integrate(LL,0.0001,m1)$value
  }
  
  
  U=function(T,m1,m2){
    UT1(T,0,0)*UT2(T,0,0)+(UU1(T,m1,m2)*UU2(T,m1,m2))
  }
  

  Cost=function(T,m1,m2){
    (cu*U(T,m1,m2)+2*((c1+cr)*ff1(T,m1)+(c2+cr)*ff2(T,m2)-cr*ff12(T,m1,m2))+
       cp*floor(t0/T))/t0
  }
 
  
Cost=Vectorize(Cost)



x = seq(0.1,4,0.25)
length(x)
y = seq(0.1,7,0.25)

z = matrix(data=NA, nrow=length(x), ncol=length(y))
for(i in 1:length(x))
{
  for(j in 1:length(y))
  {
    
    z[i,j] = Cost(3,x[i],y[j])
  }
}




Xx=rep(x, each=length(y)) 
Yy=rep(y,times=length(x))
xxx=matrix(Xx,length(x),length(y),byrow = TRUE)
yyy=matrix(Yy,length(x),length(y),byrow = TRUE)
zzz=matrix(z,length(x),length(y),byrow = TRUE)


surf3D(x=xxx, y=yyy, z,  plot = TRUE, bty="b2", phi = 10, theta = 110,
       xlab = "", ylab="", zlab="Cost function for FGM copula",
       ticktype = "detailed",colkey = FALSE,
       cex.axis=0.8, cex.lab=1)
text3D(72, 67, 20, labels = expression(M[2]), add = TRUE, adj = 1,cex=1.1)


text3D(75, 1, 15, labels = expression(M[1]), add = TRUE, adj = 1,cex=1.1)


contour(x,y,z,col = "darkblue",method = "flattest",xlab=expression(omega[1]),
        ylab=expression(omega[2]),cex.lab=1.2)
grid( ny = 6, col = "lightgray", lty = 1,
      lwd = 0.01, equilogs = TRUE)<pre>An Error occurred while handling another error:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/ErrorHandler.php(136): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/ErrorHandler.php(135): yii\web\ErrorHandler-&gt;renderException()
#3 [internal function]: yii\base\ErrorHandler-&gt;handleException()
#4 {main}
Previous exception:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(656): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/models/Files.php(696): yii\base\Application-&gt;end()
#3 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/controllers/FilesController.php(484): faravaghi\filemanager\models\Files-&gt;getFile()
#4 [internal function]: faravaghi\filemanager\controllers\FilesController-&gt;actionGetFile()
#5 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/InlineAction.php(57): call_user_func_array()
#6 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Controller.php(180): yii\base\InlineAction-&gt;runWithParams()
#7 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Module.php(528): yii\base\Controller-&gt;runAction()
#8 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Application.php(103): yii\base\Module-&gt;runAction()
#9 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(386): yii\web\Application-&gt;handleRequest()
#10 /var/www/html/prof-homepages/frontend/web/index.php(18): yii\base\Application-&gt;run()
#11 {main}</pre>