%% ECE 602 PROJECT: Operation and Configuration of Storage Portfolio via Convex Optimization


%% Problem formulation - Configuration Problem
% The objective of the configuration problem is to find the optimal
% combination of storage units to minimize the average long-term utility cost. For
% this purpose, the paper suggests to evaluate the long-term cost of all possible portfolio congigurations  
% consisting of 3 storage units, being the available choices: 0,1,2 or 3 units of each type (Small S, Medium M, and Large L), Hence
% $$4^3$ different possible combinations.
% The average operational cost will be optimizing for one year which
% requires solving the optimal operational problem (described in the
% Operation.m) file every half hour for one year. In fact, from the equation of
% the average operational cost: 
%%
% $$\mathrm{minimize\hspace{0.5cm}} \frac{1}{T} \sum_{\tau=t}^{t+T-1}\hat{p}_{\tau|t}\hat{s}_t+\alpha(\hat{r}_{\tau|t}-\hat{d}_t)_{+}$
%%
% it is implied that, the longer the period of simulation, the more exact
% the average operational cost is, as $$T\to\infty$. The anual-levelized captial
% cost of the units are given, so the total anaul cost of the portfolio can be found
% as the multiplication of the number of units of each type by their respective 
% levelized-cost. Finally, a plot of the the operational cost vs the
% capital cost for each possible portfolio configuraiton is presented, and the
% pareto optimal boundary is highlighted which can be interpretaed as the optimal solutions. 
%%
% Since most of the code used to solve this problem is the same as for the
% operational prblem, most of the comments will be the same as in the operation.m file; however, we would 
% like to point out that some additional code is include in this file, so additional comments have been 
% added as well.

%% Input data generation
% In this section, the demand and prices for a complete year of operation
% of the storage portfolio will be generated, using the stochastic models
% provided in the papaer. Likewise, some general parameters defined in the model are declared here.
% All these values presented have been taken directly from the original paper.
% 
%%
% The variable Ts contains the total simulation time used to comptue the
% average operational cost. In this case, we will use 1 month rather than 1
% year to find this value. The reason for this is that the simulation time
% required to solve this problem is extremely high, as highlighted by the
% authors of the original paper. Indeed, they use a special software to
% solve this problem since it can take more than 15 hours to solve it using
% basic solvers a CVX. Nevertheless, as we will demonstrate in the results, the differences obtained 
% based on this asumption does not drastically affect the solution when
% compared to the one provided in the original paper. The reason for this
% is that since the cyclic behaviour of the demand remains consistent in time, i.e., 
% does not change very much during the whole year, the value of the operational
% cost tends to stabilize. 
%%
Ts=2*24*30;        %Simulation time: 1 month simulations
%%
% The variable T stablishes the period for future predictions of demand and price rquired for the RHC
% control algorithm, i.e., 24-hours beyond "t" or 48 periods of half-hour.
%%
T=48;                   %Period for future prediction in the RHC control
%%
% The variable Tt defines the total number of periods for which we want to
% forecast the price and demand. We can still forecast the demand for one
% year, since the computational burden for this is relatively low.
Tt=2*12*30*24;          %2 subperiods of 30 min for each hour
%%
% T_in contains the period that we will later consider as the starting point for the operatinoal problem 
T_in=100;      
%%
% The vectors "r" and "p" will store the actual demand and price. At every
% period "t" the RCH algorithm will assume that the values stored in "r" and "p"
% ar the historical data from which the forecast will be performed. However,
% "r" and "p" will be artificially generated using the procedure, also
% suggested in the paper.
r=zeros(Tt,1);          %Actual Request (demand)
p=zeros(Tt,1);          %Actual Price 
%%
% The vector "v" will contain the optimal values of the primal varibales for the
% whole period of simulation
v=zeros(Tt,11);         %Vector of answers [dt st up_s1 up_s2 up_s3 um_s1 um_s2 um_s q(t+1)_s1 q(t+1)_s2 q(t+1)_s3]
%%
% The variables $$J_{op}$ and $$J_{cap}$ will store the average operational
% cost and captial cost of each of the portfolio alternatives.
Jop=zeros(64,1);
Jcap=zeros(64,1);

%%
% The stochastic processes of the price and demand are modeled as:
% 
% $$\log(r_t)=a_r+b_r\cos(2\pi t/48 - 5\pi/4) + u_t+x_t$$
%%
% $$\log(p_t)=a_p+b_p\cos(2\pi t/48 - 3\pi/2) + u_t+y_t$$
% 
%%
% $$u_t=hu_{t-1} + z_t$$
% 
%%
% Where $$x_t$, $$y_t$, and $$z_t$ are random variables normally
% distributed with $$\mathcal{N} (0,\sigma^2)$.
% The code to set the parameters of these processes is presented below. These parameters have been taken from the original papaer.
ar=ones(Tt,1)*0.2;
br=ones(Tt,1)*0.4;
ap=ones(Tt,1)*0.15;
bp=ones(Tt,1)*0.4;
h=0.9;

%%
% The term $$u_t$ is a first order autoregressive proces. The code to
% generate the values of $$u_t$ is presented below
x=randn(Tt,1)*(0.01)^(1/2);
y=randn(Tt,1)*(0.01)^(1/2);
z=randn(Tt,1)*(0.01)^(1/2);
u=zeros(Tt,1);
u(1)=z(1); 
    for i=2:Tt
        u(i)=h*u(i-1)+z(i);
    end
  
%% 
% The historical values of "p" and "r" are generated using the equations presented before. The variable "t" contins the index of the total simulation time, and will be
% used for ploting.
t=[1:Tt]';
r=exp(ar+br.*cos(2*pi*t/T-5*pi/4)+u+x);
p=exp(ap+bp.*cos(2*pi*t/T-3*pi/2)+u+y);

%%
% The code to calculate all possible portfolio configurations is presented below 
%Configuration of all possible portafolios:
opt=0:1:3;
New_D=kron(ones(1,3),opt);
Units=unique(nchoosek(New_D,3),'rows');
Units=Units+0.005;

%%
% Once all the possible portafolio configurations are found, the levelized-average cost of each
% is found using the code shown next.
Jcap=round(Units*[5;3;2]);

%% Proposed solution
% In this section, the code to solve the problem will be presented. The
% mathematical expressions that describe the optimization problem will be
% commented.
%% 
% The following loop is used to calculate average operational cost 
% for each one of the 64 possible portfolio configurations.
 
    for k=1:64
%%
% The follwing code set the initial state of charge for each of the storage
% units in the portfolio
    v(T_in-1,9:11)=[Units(k,1)*2 Units(k,2)*2 Units(k,3)*1]; 
%%
% The loop over the variable Tau solves the optimization problem for one
% month (1440 periods) for the operation problem.
%

        for Tau=T_in:T_in+Ts   
    
%%
% *Demand and price forecast at time t*
%%
% This portion of the code generates the demand and and price forcast at
% time=t. The forecast uses the information of the latest 24 hours (48 values) befor "t".
% The forecast estimate the demand and price for the next 24 hours after "t". The results are
% stored in two new variables: $$\hat{r_{\tau \mid t}}$ and $$\hat{p_{\tau
% \mid t}}$ respectively. 
            
        r_p=r; %r_p contains the projected demand
        p_p=p; %r_p contains the projected price
%%
 % The expressions used to forecast the value of the demand and price are
 % presented below:
%%
 % $$\hat{r}_{\tau|t}=\mathrm{exp}(\mu_{r_{\tau|t}}+\sigma_{r_{\tau|t}}^2/2)$$
%%
% $$\hat{p}_{\tau|t}=\mathrm{exp}(\mu_{p_{\tau|t}}+\sigma_{p_{\tau|t}}^2/2)$$
%% 
% The code use to reproduce these two equations is shwon below:
        
            for i=1:T-1 
                mean_ur=mean(log(r_p(Tau+i-T:Tau+i)));
                mean_up=mean(log(p_p(Tau+i-T:Tau+i)));
                var_ur=var(log(r_p(Tau+i-T:Tau+i)));
                var_up=var(log(p_p(Tau+i-T:Tau+i)));
                r_p(Tau+i)= exp(mean_ur+var_ur/2);
                p_p(Tau+i)= exp(mean_up+var_up/2);
            end

            %%
% *Optimization Problem*
% The formulation of the convex optimization problem to solve the operational cost is presented below. 
%%
% 
% $$\mathrm{minimize\hspace{0.5cm}} \frac{1}{T} \sum_{\tau=t}^{t+T-1}\hat{p}_{\tau|t}\hat{s}_t+\alpha(\hat{r}_{\tau|t}-\hat{d}_t)_{+}$$
%
%%
% $$\mathrm{subject\hspace{0.1cm}to\hspace{0.1cm}}
% \hat{q}_{\tau+1}=\eta^l\circ\hat{q}_{\tau}+\eta^c\circ\hat{u}_{\tau}^{+}-(1/\eta^d)\circ\hat{u}_{\tau}^{-}+\hat{w}_{\tau|t}$$
% 
%%
% $$\hat{d}_{\tau}-\hat{s}_{\tau}+\textbf{1}^T
% u_{\tau}^{+}-\textbf{1}^T u_{\tau}^{-}=0$$
%%
% $$0\leq\hat{q}_{\tau}\leq Q,\hspace{0.3cm}0\leq u_{\tau}^{+}\leq C,$$
%
%%
% $$0\leq u_{\tau}\leq D,\hspace{0.3cm}\tau = t,...,t+T-1$$
%
%%
% $$\hat{q}_t=q_t, \hspace{0.3cm}\hat{q}_{t+T}=q_{final}$$
%
%%
% $$0\leq s_t \leq s_{max}$$
%
%%
% Where:
%%
% 
% * $$T$ is the number of steps in the future that the algorithm projects the price and the demand
% * $$\hat{p}_{\tau|t}$ is the forecasted price at time $$\tau$ based on the information available at time $$t$
% * $$\hat{r}_{\tau|t}$ is the forecasted demand at time $$\tau$ based on the information available at time $$t$
% * $$\hat{s}_{t}$ is the estimated energy purchased from the system at time $$t$
% * $$\hat{d}_{t}$ is the estimated energy injected into the grid at time $$t$
% * $$\alpha$ is the cost of the unsuplied energy
% * $$\hat{q}_{\tau}$ is the estimated charging level of each storage unit at time $$t$
% * $$q$ is the actual storage level of each storage unit at time $$t$, 
% * $$\eta^l,\eta^c,1(/\eta^d)$ are the leakage, charging and discharging efficiencies of the storage unit respectively
% * $$\hat{u}_{\tau}^{+}$ and $$\hat{u}_{\tau}^{-}$ are the estimated charging and
% discharging rates (inflows and outflows) of the storage units at time $$t$
% * $$q_{final}$ is the state of charge of the portfolio at the last step of the simulation
% * $$Q,C,D$ are the unit storage capacity, maximum charging, and maximum discharging rates of the storage
% units respectively.
% * $$s_{max}$ is the maximum energy that can be purchased from the system
% * $$\hat{w}_{\tau|t}$ is the an additional imput of energy to the storage units. For the simulations, the authors neglect this term.

%%
% In this problem, the optimization variables are:  $$\hat{d}$, $$\hat{s}$, $$\hat{q}$, $$\hat{q}$, $$\hat{u}_{\tau}^{+}$, and $$\hat{u}_{\tau}^{-}$ 

%%
% As it can be seen in the problem formulation, the objective function
% includes two terms: the $$\hat{p}_{\tau|t}\hat{s}_t$ term represents the
% cost of purchasing energy at price $$s_t$, while the term
% $$\alpha(\hat{r}_{\tau|t}-\hat{d}_t)_{+}$ represent the cost of not
% supplying the demand. If the whole demand is supplied, the term
% $$\hat{r}_{\tau|t}-\hat{d}_t$ becomes zero; otherwise, it contains the
% unsupplied demand. Since the unsupplied deamnd can only be positive, the
% this later term has been replaced in the cost function with the variable
% "g" which has been declared as nonegative; hence, an additional
% constraint has be included in the model to defined the value of "g" as
% $$g=\hat{r}_{\tau|t}-\hat{d}_t$. Furthermore, the operational cost is the devided by the total simulation time T 
% so that an average hourly cost is found. From the convex optimization theroy, it is 
% clear that the factor $$(1/T)$ does not affect the optimal solution.The code used to solve the optimization problem
% is presented below. The solver CVX has been used for this purpose.

%%    
            
            Smax=1.5*ones(T,1);
            alpha=20;
            ef_1=[0.98*ones(T,1) 0.99*ones(T,1) 0.995*ones(T,1)];
            ef_c=[0.8*ones(T,1) 0.9*ones(T,1) 1*ones(T,1)];
            ef_d=[(0.8)*ones(T,1) 0.9*ones(T,1) 1*ones(T,1)]; %This vector is 1/ef_d
            Q=[Units(k,1)*5*ones(T,1) Units(k,2)*2*ones(T,1) Units(k,3)*1*ones(T,1)];
            C=[Units(k,1)*0.75*ones(T,1) Units(k,2)*0.5*ones(T,1) Units(k,3)*0.5*ones(T,1)];
            D=[Units(k,1)*0.75*ones(T,1) Units(k,2)*0.5*ones(T,1) Units(k,3)*0.5*ones(T,1)];
            qfinal=[Units(k,1)*1 Units(k,2)*1 Units(k,3)*1];

            cvx_solver SDPT3
            cvx_begin quiet

            variables d(T,1) s(T,1) q(T+1,3) up(T,3) um(T,3)
            variable g(T,1) nonnegative

            minimize ((1/T)*(p_p(Tau:Tau+T-1)'*s + alpha*sum(g)))
            r_p(Tau:Tau+T-1)-d==g;
            q(2:T+1,:)==q(1:T,:).*ef_1+ef_c.*up-ef_d.*um;
            d-s+up*ones(3,1)-um*ones(3,1)==0;
            zeros(T,3)<=q(1:T,:)<=Q;
            zeros(T,3)<=up<=C;
            zeros(T,3)<=um<=D;
            q(1,:)==v(Tau-1,9:11);
            q(T+1,:)==qfinal;
            0<=s<=Smax;

            cvx_end
%%
% The values of the parameters have been taken from the original papaer.
%% 
% At every iteration, as the RCH algorithm requires, the optimal value of
% the state variables at time $$t$ is recorder in the vector
% $$v_t^*=(d_t,s_t,u_t^{+},u_t^{-})$, while the other optimal values for
% the rest of the period $$t=t+1,t+2,...t+T-1$ are discared. This vector
% will record the optimal value of the state variables for every iteration,
% i.e., for every step in the whole simulation period. The code introduced
% to record the results is presented below.
%%
            v(Tau,:)=[d(1,1) s(1,1) up(1,:) um(1,:) q(2,:)];
        end
        
 %%
 % The following code stores the optimal operational cost of the portfolio
 % configuration $$k$ in the vector $$Jop^k$ by plugging in the optimal variables
 % in the objective function:
        Jop(k)=(1/Ts)*(p(T_in:T_in+Ts)'*v(T_in:T_in+Ts,2) + alpha*sum(max(r(T_in:T_in+Ts)-v(T_in:T_in+Ts,1),0)));
 %%
 % End of the portfolio configuration loop
    end

%% Results
% Once the operational cost has been found for all the 64 configurations, the Pareto optimal is found
% using the code shown below:
%%
Opt_sol=sortrows([Jcap Jop]);

Boundary=ones(length(Opt_sol(:,1)),1);
Boundary(1)=1;
pivot=1;
c=1;

    for i=2:length(Opt_sol(:,1))
        if (Opt_sol(i,2)<=Opt_sol(pivot,2))
        pivot=i;
        c=c+1;
        Boundary(c)=i;
        end
    end
 %%
 % The pareto frontier stablishes configurations that cannot be improved without affecting the optimal 
 % cost (it cannot be reduced), so every point in this frontier is optimial in a multivariable
 % optimization problem.
    
    All_Points_Index=[1:length(Opt_sol(:,1))]';
    Other_Points_Index = setdiff(All_Points_Index,Boundary);
    figure('Color',[1 1 1]);
    scatter(Opt_sol(Other_Points_Index,1),Opt_sol(Other_Points_Index,2));
    hold on
    scatter(Opt_sol(Boundary,1),Opt_sol(Boundary,2),'filled');
    scatter(Opt_sol(22,1),Opt_sol(22,2),'filled','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',[0 0 0]);
    hold off
    xlabel('Jcap');
    ylabel('Jop');
    title('Capital cost versus average operating cost');
    axis([-inf,+inf,-inf,+inf]);

 %%
 % The points in red show the portafolio configurations that are part of the pareto frontier, while the other points represent
 % the rest of configurations. The point in black is the optimal solution for the configuration presented for the the
 % operation problem (1 small, 1 medium, and 1 large unit). 
 %%
 % As it can be seen, the costs asociated to each option in the pareto frontier is very close to the original study as well. 
 % However, the configurations that are part of this frontier are slightly
 % different from the original study. The reason for this is the assumption
 % made of llenght of T (number of problmes solved to find the average oeprational cost).
 %% Conclusions
%
% * As demonstrated in the simulations, we were able to succesfully reproduce the results presented in the
% papaer for the configuration problem.
% * The use of the CVX toolboxk allowed us to solve the problem; however, the number of optimization problems 
% required to calculate the average operational cost were reduced from 1 year (17520) to 1 month (1440)  due to the computational burden it represented (as indicated in the paper). 
% Nevertheless, the results we obtained were very close to the shown in the
% original papaer.
% * A possible improvement in the algorithm could be the introduction of a
% variable that controls the convergence of the operational cost in such a
% way that if it detects that there is not an important variation of the
% operation cost, it stops the loop and records the optimal cost. In this
% way, important amount of computational time could be saved.

