ploton=1;
simdays = 1; %controls the amount of days to simulate for.
yolo = 1;
%Default arguments: 
%Reads the function arguments. If not specified the defualt args are used
mode = 1; V_max = 200; ploton=0; backoff=0;
warning off
%clc
%clear all
close all

%% Initializing
global par
par.c_p  = 4.19;                      % heat capacity of water [kJ/(kgC)]
par.rho  = 1;                                   % density of water [kg/l]
par.Q_loss = 0;                              % loss of heat from the tank

%disturbances and set point limitations
T_hws  = 50;      % Temperature wanted on the hot water reaching the user
T_cw   = 5;                               % Temperature of the cold water
T_smax = 90;                             % Maximum tepmperature set point
T_smin = T_hws+backoff;           % adjusts minimum temperature set point  
V_min  = 50;                                  % minimum water in the tank
%V_max is defined on line 24

totfailures = [];
Ekfailures = [];
Eendy= [];

%simulation parameters (initial conditions)
T_s = T_smax;        %For base case
par.V = V_max;       %volume starting point
par.T = T_smax;      %Temperature starting point

%Controller tunings
Kc1  = 6;    % gain for controller of volume
Kc2  = 0.21; % gain for controller of temperature
tau1 = 1/80; % Really 80 but simulink uses tuning paramer I=1/tau_i in its
tau2 = 1/80; % PI controllers.
%input limitations
q_in_min = 0; % dm^3/s
q_in_max = 10;% dm^3/s
Q_min    = 0; % kJ/s
Q_max    = 5; % kJ/s

%Generating the demand/price and the setpoints:
%
% In general we prefer to load a pregenerated profiles (longprofile)
% to ensure that we get comparable results for all experiments.
% For the hot water setpoint T_hws and cold water temperature T_cw a 
% constant value is wanted. These are created using the "ext" variable. 
% The electricity price are loaded from a price.csv file.
%
test=load('langprofil.mat');
a = test.testlong;   %contains the pregrnerated demand and price profiles
profile.time=(a.time'+0.167)*3600; 
ext = ones(length(a.time),1);    % used to extend the constant set points
sttime= 2;
ext2= ones(length(a.time)-sttime,1);
ext3= zeros(sttime,1);
demand = [ext3; 2.5*ext2]
%Generating variable containing all disturbances for the given period:
profile.signals.values =  [demand, T_hws.*ext, T_cw.*ext, a.price];

simlength = 24*simdays;                 % number of hours to simulate for

% Load pregenerated average energy demand and average price profiles, 
% this is just used in case II:
av.demand = csvread('average.csv',0,0,[0,0,23,0]);       % q_hw, dm^3/min
predictedprice = getPriceHour(4,1,4,31);    % 31 because we need the full
                      % prediction horizon at the end of the 30th day too
switch mode  
    case 1                      % CONSTANT TEMPERATURE CASE, T_s = T_smax
        T_s =T_smax;                             % constant max set point
        
        % Define the matrices used to store state variables
        ptime_base=[]; %changing time from seconds to hours.
        sp_base=[]; % extracting setpoints
        y_base=[]; % T, V and modified cost fx
        d_base=[]; % disturbances
        u_base=[]; % inputs
        
        mdl = 'hwtank_plot.mdl';
        load_system(mdl);
        set_param('hwtank_plot/Temp_s', 'Value', num2str(T_s));
        out1 = sim(mdl,'SrcWorkspace','current','StopTime', ...
            (num2str(3600*10.7)), 'SaveFinalState', 'on', ...
            'LoadInitialState', 'off', 'SaveCompleteFinalSimState', ...
            'on','FinalStateName', 'xFinal');
        
        %update the state variables
        ptime_base=[out1.get('time')/3600]; %time from seconds to hours.
        sp_base=[out1.get('sp')]; %extracting setpoints
        y_base=[out1.get('y')]; % T, V and modified cost fx
        d_base=[out1.get('d')]; % disturbances
        u_base=[out1.get('u')]; % inputs
        
            figure(2)
            subplot(611)
            plot(ptime_base, sp_base(:,1),'green', ptime_base, y_base(:,1))
            set(gca,'FontSize',14)
            title('Model implementation ')
            ylabel('V, dm^3')
            
            subplot(612)
            plot(ptime_base,u_base(:,2),'red')
            set(gca,'FontSize',14)
            ylabel('q_i_n, dm^3/min')
            hold on
            
            subplot(613)
            plot(ptime_base, sp_base(:,2),'green', ptime_base,y_base(:,2))
            set(gca,'FontSize',14)
            ylabel('T, ^oC')
            hold on
            
            subplot(614)
            plot(ptime_base,u_base(:,1),'red')
            set(gca,'FontSize',14)
            ylabel('Q, kW')
            hold on
            
            subplot(615)
            plot(ptime_base,d_base(:,4))
            set(gca,'FontSize',14)
            ylabel('P, EUR/MWh')
            hold on
            
            subplot(616)
            plot(ptime_base,y_base(:,3))
            set(gca,'FontSize',14)
            ylabel('J, EUR')
            % hold on      
            figure(3)
            subplot(311)
            plot(ptime_base,d_base(:,1))
            set(gca,'FontSize',14)
            title('Hot water demand')
            ylabel('q_o_u_t, dm^3/min')
            hold on
            
            subplot(312)
            plot(ptime_base,u_base(:,2),'red')
            set(gca,'FontSize',14)
            ylabel('q_i_n, dm^3/min')
            hold on
            
            subplot(313)
            plot(ptime_base,y_base(:,1))
            set(gca,'FontSize',14)
            ylabel('V, dm^3')
        
        
        % If a 30 day simulation, save the data points to .mat file
        if yolo == 1
            if exist(['switchplot.mat'],'file')
                save(['switchplot.mat'],'ptime_base',...
                'sp_base','y_base','d_base','u_base','-append');
            else 
                 save(['switchplot.mat'],'ptime_base',...
                'sp_base','y_base','d_base','u_base');
            end
        end
                
        % Variables for calculating T_violation and J_end
        listlength= length(y_base(:,2));
        yvar = y_base;
end
