Updates Simulink code for easier reuse in GP identification/predictions

Simulink/PolydomeCARNOT.prj

@@ -1,2 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<MATLABProject xmlns="http://www.mathworks.com/MATLABProjectFile" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0"/>

Simulink/main.m

@@ -1,93 +0,0 @@
-clear all
-close all
-clc
-%%%%%%%%%%%%%%%%%%%%%%%
-
-%% Load the experimental data
-Exp1 = load("../Data/Luca_experimental_data/Exp1.mat");
-py_Exp1 = load("../Data/Exp1_WDB.mat");
-
-tin = py_Exp1.Exp1_WDB(:,1);
-
-% The power trick: when the setpoint is larger than the actual temperature
-% the HVAC system is heating the room, otherwise it is cooling the room
-Setpoint = Exp1.Exp1.Setpoint.values;
-InsideTemp = mean([Exp1.Exp1.InsideTemp.values, Exp1.Exp1.LakeTemp.values], 2);
-OutsideTemp = Exp1.Exp1.OutsideTemp.values;
-HVAC_COP = 4.5;
-Heating_coeff = sign(Setpoint - InsideTemp);
-Heating_coeff(Heating_coeff == -1) = -1 * HVAC_COP;
-
-
-%% Set the model parameters
-
-% Large side windows
-window_size = [2 25];
-window_roof_size = [5 5];
-surface_part = 0.1;
-U = 1.8; % heat transfer coefficient [W/m2K]
-g = 0.7; % total solar energy transmittance
-v_g = 0.65; % transmittance in visible range of the sunlight
-
-% Roof
-wall_size = [25 25];
-roof_position = [0 0 0];
-% The roof is supposed to be made of [5cm wood, 10cm insulation, 5cm wood]
-node_thickness = [0.05 0.10 0.05]; % Data from 03.03 email with Manuel
-% Data from https://simulationresearch.lbl.gov/modelica/releases/latest/help
-% /Buildings_HeatTransfer_Data_Solids.html#Buildings.HeatTransfer.Data.Solids.Plywood
-node_conductivity = [0.12 0.03 0.12]; 
-node_capacity = [1210 1200 1210];
-node_density = [540 40 540];
-
-% Floor
-ceiling_size = [25 25];
-% The floor is supposed to be made of []
-layer_thickness = [0.05 0.10 0.20];
-layer_conductivity = [0.12 0.03 1.4];
-layer_capacity = [1210 1200 840];
-layer_density = [540 40 2240];
-
-
-
-%% Set the run parameters
-
-air_exchange_rate = tin;
-air_exchange_rate(:,2) = 2.0;
-t0 = 24;
-
-power = [tin Heating_coeff .* (Exp1.Exp1.Power.values - 1.67 * 1000)];
-
-Te = 60*60*24*365;
-
-%% Run the simulation
-% Note: The simlulink model loads the data separately, includes the
-% calculated solar position and radiations from pvlib
-simout = sim("polydome_model_1");
-
-%% Compare the simulation results with the measured values
-SimulatedTemp = simout.SimulatedTemp.Data;
-
-figure; hold on; grid minor;
-plot(tin, InsideTemp);
-plot(tin, OutsideTemp);
-plot(simout.tout, SimulatedTemp, 'LineWidth', 2);
-plot(tin, Setpoint);
-legend('InsideTemp', 'OutsideTemp', 'SimulatedTemp', 'Setpoint');
-hold off;
-
-% calculation notes for furniture wall parameters
-
-% surface:
-% 1/4 * 1.8 [m2/m2 of floor space] * 625 m2 surface = 140 m2
-% 140 m2 = [7 20] m [height width]
-
-% mass:
-% 1/4 * 40 [kg/m2 of floor space] * 625 m2 surface = 6250 kg
-
-% volume:
-% 6250[kg]/600[kg/m3] = 10.41 [m3]
-
-% thickness:
-%10.41[m3]/140[m2] = 0.075m = 7.5cm
-

Simulink/model_identification.m

@@ -0,0 +1,78 @@
+clear all
+close all
+clc
+%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Load the experimental data
+exp_id = "Exp1";
+exp_path = strcat("../Data/Luca_experimental_data/", exp_id,".mat");
+wdb_path = strcat("../Data/Experimental_data_WDB/", exp_id, "_WDB.mat");
+
+Exp_data = load(exp_path);
+load(wdb_path);
+
+% Save the current WDB to the Simulink model import (since Carnot's input file is hardcoded)
+save('../Data/input_WDB.mat', 'Exp_WDB');
+
+tin = Exp_WDB(:,1);
+
+% The power trick: when the setpoint is larger than the actual temperature
+% the HVAC system is heating the room, otherwise it is cooling the room
+Setpoint = Exp_data.(exp_id).Setpoint.values;
+InsideTemp = mean([Exp_data.(exp_id).InsideTemp.values, Exp_data.(exp_id).LakeTemp.values], 2);
+OutsideTemp = Exp_data.(exp_id).OutsideTemp.values;
+
+HVAC_COP = 3;
+Heating_coeff = sign(Setpoint - InsideTemp);
+Heating_coeff(Heating_coeff == -1) = -1 * HVAC_COP;
+
+
+%% Set the model parameters
+
+surface_part = 0.1;
+
+%% Set the run parameters
+
+air_exchange_rate = tin;
+air_exchange_rate(:,2) = 2.0;
+
+% Set the initial temperature to be the measured initial temperature
+t0 = Exp_data.(exp_id).InsideTemp.values(1);
+
+% 
+power = [tin Heating_coeff .* (Exp_data.(exp_id).Power.values - 1.67 * 1000)];
+
+% Turn down the air exchange rate when the HVAC is not running
+night_air_exchange_rate = 0;
+air_exchange_rate(abs(power(:, 2)) < 100, 2) = night_air_exchange_rate;
+
+
+%% Run the simulation
+% Note: The simlulink model loads the data separately, includes the
+% calculated solar position and radiations from pvlib
+load_system("polydome");
+set_param('polydome', 'StopTime', int2str(tin(end)));
+simout = sim("polydome");
+
+%% Compare the simulation results with the measured values
+SimulatedTemp = simout.SimulatedTemp.Data;
+
+figure; hold on; grid minor;
+plot(tin, InsideTemp);
+plot(tin, OutsideTemp);
+plot(simout.tout, SimulatedTemp, 'LineWidth', 2);
+
+
+legend('InsideTemp', 'OutsideTemp', 'SimulatedTemp');
+
+x0=500;
+y0=300;
+width=1500;
+height=500;
+set(gcf,'position',[x0,y0,width,height]);
+title(exp_id);
+%title(sprintf('Night Air exchange rate %f', night_air_exchange_rate));
+
+hold off;
+
+saveas(gcf, strcat(exp_id, '_simulation'), 'svg')