%% Figure 2
%1
v_b=7 
mu_b=1 
alpha=-0.5
beta=1 
C=0.5
x = []
syms  theta
for lambda  = 0.0: 0.1: 3.0
    syms theta
  s=max(solve(v_b+alpha*mu_b-alpha*lambda/theta-2*sqrt(C*alpha/theta)-2*beta,theta))
    x = [x,s]
end 
lambda  = 0: 0.1: 3;
figure;
subplot(221)
plot(lambda,x)
legend('\theta_{b}^{SR}')
xlabel('\lambda')
ylabel('\theta')
title('i')
%2
v_b=7 
mu_b=1 
alpha=0.5
lambda=1 
C=0.5
x = []
syms  theta
for beta = 0.0: 0.1: 3
    syms theta
    s=max(solve(v_b+alpha*mu_b-alpha*lambda/theta-2*sqrt(C*alpha/theta)-2*beta,theta))
    x = [x,s]
end 
beta  = 0: 0.1: 3;
subplot(222)
plot(beta,x)
legend('\theta_{b}^{SR}')
xlabel('\beta')
ylabel('\theta')
title('ii')
%3
v_b=7 
mu_b=1 
beta=1 
lambda=1 
C=0.5

x = []
syms  theta
for alpha = 0.0: 0.1: 3
    syms theta
   s=max(solve(v_b+alpha*mu_b-alpha*lambda/theta-2*sqrt(C*alpha/theta)-2*beta,theta))
    x = [x,s]
end 
alpha  = 0: 0.1: 3;
alpha(1)=[]
subplot(223)
plot(alpha,x)
legend('\theta_{b}^{SR}')
xlabel('\alpha')
ylabel('\theta')
title('iii')
%4 
v_b=2
mu_b=7
beta=1 
lambda=1 
C=0.5
x = []
syms  theta
for alpha = 0: 0.1: 3
    syms theta
    s=max(solve(v_b+alpha*mu_b-alpha*lambda/theta-2*sqrt(C*alpha/theta)-2*beta,theta))
    x = [x,s]
end 
alpha  = 0: 0.1: 3;
subplot(224)
plot(alpha,x)
legend('\theta_{b}^{SR}')
xlabel('\alpha')
ylabel('\theta')
title('iv')
%% Figure 3
clc;clear;
alpha=2.5
v_b=2.5
mu_b=2
C=0.5
lambda=0.4
theta=0.9
%
x_q=[]
x_p=[]
x_mu=[]
x_SR=[]
x_SW=[]
x_u=[]
for beta  = 1.4: 0.05: 2.4
syms q
q=solve(v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta)-beta*log(q/(1-q))-beta/(1-q),q)
p=v_b+alpha*mu_b-alpha*lambda*q/theta-2*sqrt(C*alpha/theta)-beta*log(q./(1-q))
mu=lambda*q/theta+sqrt(C/(alpha*theta))
SR=lambda*q.*p
SW=SR+(lambda*q.*beta).*log(q./(1-q))
u=beta*log(q./(1-q))
 x_q=[x_q,q]
x_p=[x_p,p]
x_mu=[x_mu,mu]
x_SR=[x_SR,SR]
x_SW=[x_SW,SW]
x_u=[x_u,u]
end 
%
(v_b+alpha*mu_b-alpha*lambda-2*sqrt(C*alpha))/2
%
q=1
C*theta/(lambda^2*q^2)
%
beta  = 1.4: 0.05: 2.4
beta(1)=[]
figure;
subplot(121)
plot(beta,x_p,'ks-')%4-5
hold on
plot(beta,x_u,'ko-')%-0.4-0.5
hold on
xlabel('\beta')
legend('P','U')
subplot(122)
plot(beta,x_q,'rs-')%0.46-0.6
hold on
plot(beta,x_mu,'bs-')%0.67-0.73
hold on
plot(beta,x_SW,'ro-')%0.85-1.05
hold on
plot(beta,x_SR,'bo-')%0.915-0.94
xlabel('\beta')
legend('q','\mu','SW','SR')
%% Figure 4
clc;clear;
alpha=2.5
v_b=2.5
mu_b=1
C=1
lambda=0.4
theta=0.9
%
x_q=[]
x_p=[]
x_mu=[]
x_SR=[]
x_SW=[]
x_u=[]
for beta  = 0: 0.1: 3
syms q
q=solve(v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta),q)
p=v_b+alpha*mu_b-alpha*lambda*q/theta-2*sqrt(C*alpha/theta)-beta*log(q./(1-q))
mu=lambda*q/theta+sqrt(C/(alpha*theta))
SR=lambda*q.*p
SW=SR+(lambda*q.*beta).*log(q./(1-q))
u=beta*log(q./(1-q))
 x_q=[x_q,q]
x_p=[x_p,p]
x_mu=[x_mu,mu]
x_SR=[x_SR,SR]
x_SW=[x_SW,SW]
x_u=[x_u,u]
end 
% >0
q=0
v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta)
% <0
q=1
v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta)
%
beta  = 0: 0.1: 3
plot(beta,x_p,'ks-')%4-5
hold on
plot(beta,x_u,'ko-')%-0.4-0.5
hold on

plot(beta,x_q,'rs-')%0.46-0.6
hold on
plot(beta,x_mu,'bs-')%0.67-0.73
hold on
plot(beta,x_SW,'ro-')%0.85-1.05
hold on
plot(beta,x_SR,'bo-')%0.915-0.94
xlabel('\beta')
legend('P','U','q','\mu','SW','SR')
%% Figure 5
clc;clear;
alpha=2.5
v_b=2.5
mu_b=1
C=1
lambda=0.4
for theta=0.89:0.02:0.93
%
x_SWR=[]
for beta  = 0: 0.2: 3
syms q
q=solve(v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta)-beta*log(q/(1-q))-beta/(1-q),q)
SW=lambda*q*(v_b+alpha*mu_b-alpha*lambda*q/theta-2*sqrt(C*alpha/theta))
x_SWR=[x_SWR,SW]
end
%
x_SWW=[]
for beta  = 0: 0.2: 3
syms q
q=solve(v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta),q)
p=v_b+alpha*mu_b-alpha*lambda*q/theta-2*sqrt(C*alpha/theta)-beta*log(q./(1-q))
mu=lambda*q/theta+sqrt(C/(alpha*theta))
SW=lambda*q.*p+(lambda*q.*beta).*log(q./(1-q))
x_SWW=[x_SWW,SW]
end 
%
poa=x_SWR./x_SWW
beta  = 0: 0.2: 3
xlabel('\beta')
ylabel('PoA')
plot(beta,poa,'s-')
hold on
end
legend('\theta=0.89','\theta=0.91','\theta=0.93')
%% Figure 6
clc;clear;
alpha=2.5
v_b=2.5
mu_b=1
C=1
lambda=0.4
theta=1
beta=2
syms q
q=solve(v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta),q)
p=v_b+alpha*mu_b-alpha*lambda*q/theta-2*sqrt(C*alpha/theta)-beta*log(q./(1-q))
u=beta*log(q./(1-q))
mu=lambda*q/theta+sqrt(C/(alpha*theta))
U=v_b+alpha*mu_b-alpha*mu-p-C/(mu*theta-lambda*q)
%theta=1 q=0.9189
n = 10000; %
lambda  = 0.4*0.9189; muA = 1;
tarr = cumsum(exprnd(1/lambda,1,n));
tsrv = exprnd(muA,1,n);
tsta = zeros(1,n);
tlea = zeros(1,n);
twat = zeros(1,n);
tsta(1) = tarr(1);
 tlea(1) = tsta(1) + tsrv(1);
    twat(1)=tlea(1)-tarr(1);
for i = 2:n
    tsta(i) = max(tarr(i),tlea(i-1));
    tlea(i) = tsta(i) + tsrv(i);
    twat(i)=tlea(i)-tarr(i);
end
mean(twat)
%
u_0=v_b+alpha*mu_b-alpha*mu-p-C*twat
mean(u_0)
plot(u_0)
hold on
A= U*ones(1,n);
plot(A)
legend('U_a','U')
axis([0 1000,-5 10])
%% Table 1
clc;clear;
alpha=2.5
v_b=2.5
mu_b=1
C=1
lambda=0.4
theta=0.99
beta=2
syms q
q=solve(v_b+alpha*mu_b-2*alpha*lambda*q/theta-2*sqrt(C*alpha/theta),q)
p=v_b+alpha*mu_b-alpha*lambda*q/theta-2*sqrt(C*alpha/theta)-beta*log(q./(1-q))
u=beta*log(q./(1-q))
mu=lambda*q/theta+sqrt(C/(alpha*theta))
U=v_b+alpha*mu_b-alpha*mu-p-C/(mu*theta-lambda*q)
%
n = 10000;
lambda  = 0.4*0.9018; muA = 1;
tarr = cumsum(exprnd(1/lambda,1,12000));
tsrv = exprnd(muA,1,12000);
tsta = zeros(1,12000);
tlea = zeros(1,12000);
twat = zeros(1,12000);
theta=unifrnd(0,1,1,12000)
tsta(1) = tarr(1);
 tlea(1) = tsta(1) + tsrv(1);
    twat(1)=tlea(1)-tarr(1);
count =1
for i = 2:length(tarr)
    if theta(i)<0.99
    if count == 10000
        break;
    end
    tsta(i) = max(tarr(i),tlea(i-1));
    tlea(i) = tsta(i) + tsrv(i);
    twat(i)=tlea(i)-tarr(i);
    count=count+1
    else
    tsta(i) = max(tarr(i),tlea(i-1));
    tlea(i) = tsta(i) + tsrv(i);
    twat(i)=tlea(i)-tarr(i);
     tarr=[tarr,tlea(i)]  ;
     tarr=sort(tarr);
    end
end
ave=sum(twat)/10000
u_0=v_b+alpha*mu_b-alpha*mu-p-C*ave