%% IDEAS: Inverse Dynamics Energy Assessment and Simulation
%  IDEAS_V5_1_Uof1.m
%  IDEAS Version 5.1 - IDEAS ESRU Baseline

%%% External Files Required
% * *For Disturbances:          To.mat, Qfree.mat 
% * *For Temperature Setpoint:  SAPsetpoint.mat

% Standard Test Case Benchmark for a dwelling with Structure 
% U-Value of 1.048 W/(m^2*K)
% Purpose : Run time now in the region of 11 seconds

%% Begin code
% clears command window, all variables, and starts simulation timer
clc;clear all;

%% Building parameters
% Define all building parameters here

% Mass of Internal Air in Zone { kg }
% SOURCE: Volume of Air * Pa (Density of Air STANDARD VALUE)
ma = 222*1.205; %267.51 kg

% Heat Capacity of Air @ -50 to 40-degC { J/(kg*K) }
% SOURCE: STANDARD VALUE
% <http://www.engineeringtoolbox.com/air-properties-d_156.html>
Ca = 1005;

% Furniture assumed to made of wood w/ a volume fraction defined by ffrac
ffrac= 0.05; Cwood = 2.25e+3; rwood = 600;
mf = rwood*ffrac*ma/1.205;  %fraction of zone that is wood.

% System Model Constants
mzone = ((1-ffrac)*ma) + mf;
Czone =  (1-ffrac)*Ca + ffrac*Cwood;

% Mass Flow Rate of Air through Internal Zone { kg/s }
% SOURCE: = based upon SAP methods
mdotv = (1-ffrac)*0.039662493237701;

% U-value of Windows { W/(m^2*K) } 
% SOURCE: FROM BREDEM CELL AB53
Uw = 1.8520;

% Surface Area of Windows { m^2 }
Aw = 16.9;

% U-value of Floor { W/(m^2*K) }
Uf = 0.7;

% Surface Area of Floor { m^2 }
Af = 44.4;

% U-value of Roof { W/(m^2*K) }
Ur = 2.3;

% Surface Area of Roof { m^2 }
Ar = 44.4;

% Wall Construction
% C { J/(kg*K) }  |  t { m }   |  r { kg/m^3 }
% ThermalConductivityCalc_MasonryCavityWall.docx
Cplaster = 1000; tplaster = 12.5e-3; rplaster = 700; 
Cconcrete = 1000; tconcrete = 100e-3; rconcrete = 2000; 
Ccavity = 1450; tcavity = 20e-3; rcavity = 30; 
Cbrick = 800; tbrick = 105e-3; rbrick = 1700;
%N.B. rbrick is density of brick, r = rho

% U-value of Building Structure { W/(m^2*K) }
Us = 1/0.954; % 1.048 W/(m^2*K)

% U-value of Internal Mass { W/(m^2*K) }
% SOURCE: based on 2 sided single brick  
% (1/Uim) = (1/7.7)+(0.12/0.8)+(1/7.7)
Uim = 2.5;

% Door U value W/m2K
Udoor = 2.1;

% Area of doors m2.
Adoor = 3.8;

% Heat Capacity of Internal Mass { J/(kg*K) }
% SOURCE: = based upon SAP methods
Cim = 1000;

% Heat Transfer Coefficient of Still Internal Air { W/(m^2*K) }
% SOURCE: BS EN ISO 6946:2007
hi = 7.6923;

% Heat Transfer Coefficient of External Air { W/(m^2*K) }
% SOURCE: BS EN ISO 6946:2007
he = 25;

% Surface Area of Structure { m^2 }
% SOURCE:
As = 81.8;

% Comfort Temperature Ratios
% CIBSE Guide A: 
% enivronmental temp = (1/3)*Tair + (2/3)*Tmean.surface
% operative temp = (1/2)*Tair + (1/2)*Tmean.surface
% SAP calibration: fzone = 1; fms = 0;
fzone  = (1);
fms = (0);

%% Calculated parameters
% All calculated building parameters

% Thickness of Structure { m }
tw = tplaster+tbrick+tconcrete+tcavity;

% Heat Capacity of Structure { J/(kg*K) }
Cs = tplaster/tw*Cplaster + tbrick/tw*Cbrick + tconcrete/tw*Cconcrete +...
    tcavity/tw*Ccavity;

% Heat Capacity of Structure {J/m^2*K }
Kt = tplaster*rplaster*Cplaster+tbrick*rbrick*Cbrick+...
    tconcrete*rconcrete*Cconcrete+tcavity*rcavity*Ccavity;

% Mass of Internal Mass in Zone { kg }
% SOURCE: 200kg/m2 of floor area
Mim = 2000*sqrt(Af);

% Surface Area of Internal Mass { m^2 }
% Approximately 2*floor area 
Aim = 20*sqrt(Af);

% Mass of Structure { kg }
ms = Kt*As/Cs;

% Mass of Internal Structure { kg }
% SOURCE: - ms
msi = ms/2;

% Mass of External Structure { kg }
% SOURCE: - ms
mse = ms/2;

% Equivalent Thermal Conductivity of Structure { W/(m*K) }
% Old model calculation still valid for new structure 
kw = he*hi*tw*Us/(he*hi-(he+hi)*Us);

%% Heater parameters
% Heating system parameters

% Maximum Heat Output of Heater { W }
Qmax = 20e+3;

% Heater Ratios
ra = 1;
rsi = 0;
rim = 0;
if ra+rsi+rim ~= 1;error('IDEAS:rChk','Incorrect Heater Ratios');end

%% State-space parameters
% Define all state-space parameters
%
% * State Differential Equation: $\dot{x}(t)=Ax(t)+Bu(t)+Ed(t)$
% * Output Equation: $y(t)=Cx(t)$

%global A B E
% System model constants
J1 = 1/(mzone*Czone);
J2 = 1/(msi*Cs);
J3 = 1/(mse*Cs);
J4 = 1/(Mim*Cim);
J5 = Uw*Aw + Uf*Af + Ur*Ar + Udoor*Adoor+ mdotv*Ca; 
J6 = Uim*Aim;
J7 = ((4*kw*hi)/((4*kw)+(hi*tw)))*As;
J8 = (2*kw/tw)*As;
J9 = ((4*kw*he)/((4*kw)+(he*tw)))*As;

% System Matrix, A (4x4)
a11 = -J1*(J5 + J6 + J7); a12 = J1*J7; a13 = 0; a14 = J1*J6;
a21 = J2*J7; a22 = -J2*(J7 + J8); a23 = J2*J8; a24 = 0;
a31 = 0; a32 = J3*J8; a33 = -J3*(J8 + J9); a34 = 0;
a41 = J4*J6; a42 = 0; a43 = 0; a44 = -J4*J6;
A = [a11 a12 a13 a14;a21 a22 a23 a24;a31 a32 a33 a34;a41 a42 a43 a44];

% Input Matrix, B (4x1)
b11 = J1*ra; b21 = J2*rsi; b31 = 0; b41 = J4*rim;
% Create column vector for Matrix F
B = [b11;b21;b31;b41];

% Output Matrix, C (1x4)
c11 = fzone + fms*(As/(Aim+As)*tw*hi/(4*kw+tw*hi)); 
c12 = fms*(As/(Aim+As)*4*kw/(4*kw+tw*hi)); 
c13 = 0; 
c14 = fms*(Aim/(Aim+As));
% Create row vector for Matrix C
C = [c11 c12 c13 c14] %#ok<NOPTS>
if  round(sum(C)*100000)/100000 ~= 1;
    error('IDEAS:cChk','Check C Matrix Ratios');end
% Feedforward Matrix, D (1x1)
d11 = 0;
D = d11;
% System Disturbance Matrix, E (4x2)
e11 = J1*J5; e12 = J1;
e21 = 0; e22 = 0;
e31 = J3*J9; e32 = 0;
e41 = 0; e42 = 0;
%Create Matrix E
E = [e11 e12;e21 e22;e31 e32;e41 e42];

% Output Disturbance Matrix, F (1x2)
f11 = 0; f12 = 0;
% Create row vector for Matrix F
F = [f11 f12];

%% RIDE controller parameters
% Define all controller parameters

% Controller Matrix, K = g*(CB)^-1
g = 1/(300);
CB = C*B; CBinv = 1/CB;
K = g*CBinv;
% ueq = -CBinv*(CA*x + CE*d)
CA = C*A;
CE = C*E;

%% Small perturbation
% Define parameters required for small perturbation

load('To.mat'); load('Qfree.mat');
Toi = To(2,1); Qfreei = Qfree(2,1);
d0 = [Toi; Qfreei];
% method: u0 = 0 (no initial heat input, calculates initial temperatures
% based upon initial disturbances)
x0 = -A\E*d0
u0 = 0;

%initial comfort
y0 = C*x0

%% Simulation parameters
% Define simulation parameters

% start time, length of simulation, time step
tstart = 0; tsim = 365*24*3600; ts = 60;
trange = tstart:ts:(tsim-ts);

% Temperature setpoint input and interpolation
load('SAPsetpoint');
Tset = SAPsetpoint;

% linear interpolation for disturbances
d = [interp1(To(1,:),To(2,:),trange,[],'extrap');...
    interp1(Qfree(1,:),Qfree(2,:),trange,[],'extrap')];

%% Simulation
% simulation solver

% Preallowcate variables
dx = zeros(4,length(trange));
ueq = zeros(1,length(trange));
uc = zeros(1,length(trange)); 
u = zeros(1,length(trange));
dy = zeros(1,length(trange));

% Small perturbation
dTset = Tset - y0;
dTset(1) = 0;

dd(1,:) = d(1,:) - d0(1);
dd(2,:) = d(2,:) - d0(2);
% simulation loop
tic;

neq = length(dx(:,1));

Y = zeros(neq,2);
F = zeros(neq,4);

%simulation loop
tic;    

for n = 1:length(trange)-1
    tspan = tstart+(n-1)*ts:ts:tstart+n*ts;
    
    % RIDE inputs
    ueq(n) = -CBinv*(CA*dx(:,n)+CE*dd(:,n));
    uc(n) = g*CBinv*(dTset(n)-dy(n));
    u(n) = min(Qmax,max(uc(n) + ueq(n),0));
    
    % Classic Runge-Kutta 4th order method
    h = diff(tspan);
    
    ti = tspan(1);
    hi = h(1);
    yi = dx(:,n);  
    
    F(:,1) = A*yi+B*u(n)+E*dd(:,n);
    F(:,2) = A*(yi+0.5*hi*F(:,1))+B*u(n)+E*dd(:,n);
    F(:,3) = A*(yi+0.5*hi*F(:,2))+B*u(n)+E*dd(:,n);
    F(:,4) = A*(yi+hi*F(:,3))+B*u(n)+E*dd(:,n);
   
    dx(:,n+1) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));

    dy(n+1) = C*dx(:,n);
    
end
toc;

%% Model check and parameter calcs
% check small perturbation has been correctly implemented

if (uc(1) == 0 && ueq(1) == 0 && u(1) == 0 && dy(1) == 0)
    disp('Small perturbation has been satisfied')
else
    error('IDEAS:sptChk','Small perturbation has failed, please check model')
end

% small perturbation + initial conditions
y = dy+y0;
x(1,:) = x0(1) + dx(1,:);
x(2,:) = x0(2) + dx(2,:);
x(3,:) = x0(3) + dx(3,:);
x(4,:) = x0(4) + dx(4,:);