%% Monte Carlo Optimization Study
clear all
syms x1 x2 n n2
tic
m_tlife = 33; 
std_tlife = 11; 
dist_tlife = makedist('Normal', m_tlife, std_tlife);
a_drate = 2; 
b_drate = 0.006;
dist_drate = makedist('Gamma', a_drate, b_drate);
a_OM = 8;
b_OM = 10;
c_OM = 20;
dist_OM = makedist('Triangular', a_OM, b_OM, c_OM);
m_Ee = 4.56;
std_Ee = 0.27;
dist_Ee = makedist('Normal', m_Ee, std_Ee);
a_pi20 = 1.50*10^3;
b_pi20 = 2.25*10^3; 
c_pi20 = 3*10^3;
dist_pi20 = makedist('Triangular', a_pi20, b_pi20, c_pi20);
m = 34;
ITC = 0.3;
pi08 = 8.50*10^3;
Pmod = 0.205; 
Asp = 40.1; 
em = 0.162;
Np = 1000;
for i = 1:Np
    tlife = round(random(dist_tlife));
    sigma = random(dist_drate);
    OM = random(dist_OM);
    Ee = random(dist_Ee);
    pi20 = random(dist_pi20);
    dpi = (pi20 - pi08)/(12*pi08);
    pi_t = @(x)dpi*x(1)*pi08 + pi08;
    Ci0 = pi08*Pmod*m;
    Ci_t = @(x)pi_t(x)*Pmod*x(2);
    Ep = Asp*em*Ee*365/m;
    Elife = matlabFunction(symsum(Ep*m*(1 - sigma)^n, n, 1, round(x1)) + symsum(Ep*x2*(1 - sigma)^n, n, 1, round(tlife - x1)) + symsum(Ep*(m - x2)*(1 - sigma)^n, n, round(x1), round(tlife)));
    Com = OM*Pmod*m*tlife;
    fun = @(x)((Ci0 + Ci_t(x))*(1 - ITC) + Com)/(Elife(x(1),x(2)));
    x0 = [0,0];
    lb = [0,0];
    up = [tlife m];
    A = [];
    b = [];
    Aeq = [];
    beq = [];
    x(i,:) = fmincon(fun, x0, A, b, Aeq, beq, lb, up);
    LCOE(i) = fun(x)*10^2;
end 

Mt = mean(x(:,1))
SD_t = std(x(:,1))
SE_Mt = SD_t/sqrt(length(x(:,1)))

MRsp = mean(x(:,2))
SD_Rsp = std(x(:,2))
SE_Rsp = SD_Rsp/sqrt(length(x(:,2)))

MLCOE = mean(LCOE)
SD_LCOE = std(LCOE)
SE_MLCOE = SD_LCOE/sqrt(length(LCOE))

% Checking upper integer and lower integer for Rsp
xl = [Mt floor(MRsp)];
LCOE_l = fun(xl)*10^2
xu = [Mt ceil(MRsp)];
LCOE_u = fun(xu)*10^2

Tsim = toc