Add files via upload

bASOCP.m

@@ -0,0 +1,102 @@
+%ECE686 Project: implementation of 
+%"Simulation-Based stochastic optimal control design and its application to
+%building control problems" Lee et al. (2018)
+%by Nurken Tuktibayev, University of Waterloo, 2019
+
+clear all;
+clc;
+tic
+%picking random day between 01 Jan 2018 and 30 Dec 2018 
+dayn=randi([1 364],1,1);
+dayn=80;%{or you can insert any day manually between 1 and 364}
+
+[inp]=inputData(dayn); %initialize state-space model for stochastic system and load weather data for dayn
+if inp.noData==1 
+    return
+else
+end;
+%%%%%%Stochastic Algorithm for beta-ASOCP
+%Algorithm settings and initial data
+N_of_iterations=10000;%NOTE: THIS CAN BE LOWERED TO SAVE TIME, BUT 10000 ITERATIONS PROVIDES "GOOD ENOUGH" RESULT
+beta = 0.001;
+Nj=20;
+Nn=1;
+delta=0.1;
+gamma=0.0001;
+Theta1=zeros(50,6);
+Theta2=zeros(1,51);
+psy1=Theta1*0;
+psy2=Theta2*0;
+
+for ii=1:N_of_iterations
+%%%%computing stochastic gradient estimate%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+g1=0;
+g2=0;
+    for j=1:Nn
+        sumplus=0;
+        summinus=0;
+        nu1=randn(50,6); %generate i-th sample of nu
+        nu2=randn(1,51); %generate i-th sample of nu
+        %computing the sample average for (Theta+beta*nu)
+        for i=1:Nj
+            w1=randi(inp.as,1,1);
+            w2=randi(inp.ds,1,1);
+            s=simulation(inp,w1,w2,Theta1+beta*nu1,Theta2+beta*nu2);
+            sumplus=sumplus+s.cost;
+        end
+        sumplus=sumplus/Nj;
+        %computing the sample average for (Theta-beta*nu)
+        for i=1:Nj
+            w1=randi(inp.as,1,1);
+            w2=randi(inp.ds,1,1);
+            s=simulation(inp,w1,w2,Theta1-beta*nu1,Theta2-beta*nu2);
+            summinus=summinus+s.cost;
+        end
+        summinus=summinus/Nj;
+        g1=g1+(1/(2*beta))*nu1*(sumplus-summinus); %sum of differences
+        g2=g2+(1/(2*beta))*nu2*(sumplus-summinus); %sum of differences
+    end
+    g1=g1/Nn; %taking average of the sum
+    g2=g2/Nn; %taking average of the sum
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%heuristic modification of Algorithm 1 
+    psy1=psy1+delta*(g1-psy1);
+    psy2=psy2+delta*(g2-psy2);
+    Theta1=Theta1-gamma*psy1;
+    Theta2=Theta2-gamma*psy2;
+    IterationsLeftInCycle1=N_of_iterations-ii %countdown
+    toc
+end
+disp('New control policy computed based on data for day');
+num2str(dayn) %training day number
+
+%Test:
+%comparing beta-ASOCP control policy with LQG 
+%for the (dayn+1) day of 2018
+[inp]=inputData(dayn);
+if inp.noData==1 
+    return
+else
+end;
+tN=1000;%number of samples
+histN=0;
+histL=0;
+for i=1:tN
+w1=randi(inp.as,1,1);
+w2=randi(inp.ds,1,1);
+%Simulation with new control policy (LQG+calculated parameter theta)
+sN=simulation(inp,w1,w2,Theta1,Theta2);
+histN(i)=sN.cost;
+%Simulation with LQG only
+sL=simulation(inp,w1,w2,Theta1*0,Theta2*0);
+histL(i)=sL.cost;
+IterationsLeftInCycle2=tN-i
+end
+%output plots for the last simulation and histogram for tN samples
+outputData(inp,sN.x,sN.u,sN.z,sN.M,'New',sN.cost,histN,max(max(histN),max(histL))+20);
+outputData(inp,sL.x,sL.u,sL.z,sL.M,'LQG',sL.cost,histL,max(max(histN),max(histL))+20);
+
+
+
+
+

inputData.m

@@ -0,0 +1,61 @@
+function [inp] = inputData(dayn)
+
+inp.N=96; %discrete time horizon
+dt=15*60; %sampling time in seconds
+
+inp.as=[32 36];  %arrival time distribution set (8-9AM)
+inp.ds=[64 76]; %departure time distribution set (4-7PM)
+
+%importing thermal probability density function
+if (exist('pdf_data.mat')==2)
+    load('pdf_data.mat','pdf_x','pdf_y');
+    inp.pdf_x=pdf_x;
+    inp.pdf_y=pdf_y;
+    inp.noData=0;
+else disp('Please load thermal pdf to pdf_data.mat and restart the program');
+    inp.noData=1;
+    return
+end
+
+%importing weather data weather2018.cvs (Source: http://weather.uwaterloo.ca/)
+if (exist('weather2018.mat')==2)
+    load('weather2018.mat','To','qsolar');
+    inp.noData=0;
+else disp('Please load weather data to weather2018.mat and restart the program');
+    inp.noData=1;
+    return
+end
+
+%loading data for selected day
+inp.dayn=dayn;
+for i=1:inp.N
+inp.To(i)=To((inp.dayn-1)*96+i);
+inp.qsolar(i)=qsolar((inp.dayn-1)*96+i);
+end
+
+%Data from Table I for 3x3 sqm room with 2.5 sqm window (from EnergyBuild)
+R1=0.0084197;
+R2=0.044014;
+R3=4.38;
+C1=9861100;
+C2=128560;
+a=0.55;
+
+%cost function data
+inp.R=0.00001;%(by experiments)
+inp.Q=[1 0 -1 0]'*[1 0 -1 0];
+
+%STATE-SPACE MODEL
+inp.x0=[21 21 21 To(1)]';
+inp.A=[1-dt/(C2*R2)-dt/(C2*R1) dt/(C2*R1) 0 dt/(C2*R2); dt/(C1*R1) 1-dt/(C1*R1)-dt/(C1*R3) 0 dt/(C1*R3); 0 0 1 0; 0 0 0 1];
+inp.B=[dt/C2; 0; 0; 0];
+inp.D=[(dt*(1-a))/C2 dt/C2 0 0;(dt*a)/C1 0 0 0;0 0 1 0;0 0 0 1];
+
+%FINITE-HORIZON LQR
+G(:,:,inp.N)=inp.Q;
+for i=inp.N-1:-1:1
+inp.F(i,:)=inv(inp.R+inp.B'*G(:,:,i+1)*inp.B)*inp.B'*G(:,:,i+1)*inp.A;
+G(:,:,i)=inp.Q+inp.A'*G(:,:,i+1)*inp.A-inp.F(i,:)'*(inp.R+inp.B'*G(:,:,i+1)*inp.B)*inp.F(i,:);
+end
+end
+

outputData.m

@@ -0,0 +1,32 @@
+function outputData(inp,x,u,z,M,method,cost,h,h_xlim)
+figure();
+subplot(5,1,1);
+hold on
+plot(x(1,:));
+plot(x(3,:),'--');
+plot(x(4,:),'-.');
+axis([1, inp.N, -40, 50]);
+xlabel('k'); ylabel('Ta Tref(-) To(-.)');
+if (method == 'New') 
+title(['beta-ASOCP control policy for day=' num2str(inp.dayn) '. Cost=' num2str(round(cost))]);
+else 
+title(['LQR control policy for day=' num2str(inp.dayn) '. Cost=' num2str(round(cost))]);
+end
+subplot(5,1,2);
+stairs(z(:),'--');
+axis([1, inp.N, 0.5, 3.5]);
+xlabel('k'); ylabel('User feel');
+subplot(5,1,3);
+stairs(M(:,1));
+axis([1, inp.N, -0.5, 1.5]);
+xlabel('k'); ylabel('Occupancy');
+subplot(5,1,4);
+plot(u(:));
+axis([1, inp.N, -1050, 1050]);
+xlabel('k'); ylabel('Control');
+subplot(5,1,5);
+histogram(h);
+xlim([0, h_xlim]);
+hold off
+end
+

readme.txt

@@ -0,0 +1,19 @@
+MATLAB Files:
+1. inputData.m:
+        Loads pdf for user actions from pdf_data.mat file
+        Loads weather data for dayn from weather2018.mat file
+        Initializes all parameters for dynamic state-space model
+        Computes state-feedback for LQR
+2. outputData.m - generates plots and histogram
+3. simulation.m - 24 hour simulation for a given weights and day
+4. bASOP.m - main script:
+        Picks random day within [1 364] range as dayn
+	Calls inputData.m for dayn
+	Sets parameters for �-ASOCP algorithm and runs it N_of_iteration times for initila theta=0
+	Simulates (dayn+1) day using new control policy with computed theta and calculates cost
+      	Simulates (dayn+1) day using state-feedback for LQR and calculates cost
+	Compares both methods in histogram for 1000 samples
+WARNING: number of iterations set to 10000, which takes ~15 min. Can be lowered for test purposes.
+
+5. weather2018.mat - real weather data from UW weather station for 2018 with step=15 minutes
+6. pdf_data.mat - probability density function of occupant actions based on indoor temperature

simulation.m

@@ -0,0 +1,72 @@
+function s = simulation(inp,w1,w2,weight1,weight2)
+%initializing data
+s.x(:,1)=inp.x0;
+s.cost=0;
+s.u=0;
+s.z=0;
+s.M=[0;0];
+for k=1:inp.N-1 %24 hours cycle (96*15 minutes)
+Pi=weight2*[1./(1+exp(-(weight1*[s.x(:,k);k;1])));1]; %new parameter
+s.u(k)=1000*satlins((-inp.F(k,:)*s.x(:,k)+Pi)/1000); %control
+%calculating day number%%%%%%%%%%%%%%
+if (k/96)<(round(k/96)) 
+dayn=round(k/96)-1;
+else
+dayn=round(k/96);
+end
+%cost function selection
+if ((k<w1+dayn*96) || (k>w2+dayn*96))
+s.cost=s.cost+s.u(k)'*inp.R*s.u(k);
+else
+s.cost=s.cost+s.x(:,k)'*inp.Q*s.x(:,k)+s.u(k)'*inp.R*s.u(k);
+end 
+%occupancy map: 1=occupant in the office, 2=occupant not in the office
+if ((k<w1+dayn*96) || (k>w2+dayn*96))
+s.M(k,1)=0;
+else
+s.M(k,1)=1;
+end
+%computing probablity of occupant action based indoor temperature using pdf
+p_index=round(s.x(1,k))-inp.pdf_x(1)+1; %scaling current temperature into [10 30] scale
+if p_index<1
+    p_index=1;
+end
+if p_index>21
+    p_index=21;
+end
+prob(1)=inp.pdf_y(p_index,1); %probability of feeling cold
+prob(2)=inp.pdf_y(p_index,2); %probability of comfort feeling
+prob(3)=inp.pdf_y(p_index,3); %probability of feeling hot
+prob = prob /sum(prob); %scaling all 3 into [0 1] scale
+random=rand;
+lowest=0;
+for i=1:3 %finding correct interval for the random value
+if((random>=0)&&(random<lowest+prob(i)))
+s.z(k) = i;
+break;
+end    
+lowest=lowest+prob(i);
+end
+%selecting ocupants action based on occupancy and temperature feelings
+w3=randi([0 1],1);
+if s.z(k)==1
+    s.M(k,2)=w3*s.M(k,1);
+elseif s.z(k)==2
+   s.M(k,2)=0;
+elseif s.z(k)==3
+    s.M(k,2)=-w3*s.M(k,1);
+end
+%state-space model
+w=[inp.qsolar(k);75+70*s.M(k,1);s.M(k,2);inp.To(k+1)-inp.To(k);];%noise
+s.x(:,k+1)=inp.A*s.x(:,k)+inp.B*s.u(k)+inp.D*w; %state
+
+if s.x(3,k+1)>30 %setpoint limit
+s.x(3,k+1)=30;
+end
+if s.x(3,k+1)<15
+s.x(3,k+1)=15;
+end
+
+end
+end
+