diff --git a/Matlab_scripts/MPCforSonja/MPCcasadi_v1_0.m b/Matlab_scripts/MPCforSonja/MPCcasadi_v1_0.m new file mode 100644 index 0000000..930ff00 --- /dev/null +++ b/Matlab_scripts/MPCforSonja/MPCcasadi_v1_0.m @@ -0,0 +1,225 @@ +classdef MPCcasadi_v1_0 < matlab.System + + % Public, tunable properties + properties(Nontunable) + TimeStep = 0; % Time step MPC + N = 0; % Planning and control horizon N + R = 1; % Weights for control cost R + T = 1; % Weights for slack variable for output constraints T + nState = 0; % Number of states X + nOut = 0; % Number of outputs Y + nIn = 0; % Number of controlled inputs U + nDst = 0; % Number of disturbance inputs + A = 0; % A + Bd = 0; % Bd (disturbance) + Bu = 0; % Bu (control) + C = 0; % C + D = 0; % D + uMin = 0; % Lower control constraints uMin + uMax = 0; % Upper control constraints uMax + yMin = 0; % Lower output constraints yMin + yMax = 0; % Upper output constraints yMax + end + + properties(DiscreteState) + + end + + % Pre-computed constants + properties(Access = private) + casadi_solver + lbg + ubg + end + + methods(Access = protected) + function sts = getSampleTimeImpl(obj) + sts = createSampleTime(obj, 'Type', 'Controllable', 'TickTime', obj.TimeStep); % Time step + end + + function num = getNumInputsImpl(~) % Number of inputs + num = 4; + end + function num = getNumOutputsImpl(~) % Number of outputs + num = 5; + end + function [dt1, dt2, dt3, dt4, dt5] = getOutputDataTypeImpl(~) % Output data type + dt1 = 'double'; + dt2 = 'double'; + dt3 = 'double'; + dt4 = 'double'; + dt5 = 'double'; + end + function dt1 = getInputDataTypeImpl(~) % Input data type + dt1 = 'double'; + end + function [sz1, sz2, sz3, sz4, sz5] = getOutputSizeImpl(obj) % OUtput dimensions + sz1 = [1, obj.nIn]; % mv + sz2 = [obj.N+1, obj.nState]; % xStar + sz3 = [obj.N, obj.nOut]; % sStar + sz4 = [obj.N, obj.nIn]; % uStar + sz5 = [1, obj.nOut]; % yStarOut + end + function [sz1, sz2, sz3, sz4] = getInputSizeImpl(obj) % Input dimensions + sz1 = [obj.nState, 1]; % xHat + sz2 = [obj.N, obj.nDst]; % disturbances + sz3 = [obj.N, 1]; % elec price + sz4 = [1, 1]; % on + end + function cp1 = isInputComplexImpl(~) % Inputs are complex numbers? + cp1 = false; + end + function [cp1, cp2, cp3, cp4, cp5] = isOutputComplexImpl(~) % Outputs are complex numbers? + cp1 = false; + cp2 = false; + cp3 = false; + cp4 = false; + cp5 = false; + end + function fz1 = isInputFixedSizeImpl(~) % Input fixed size? + fz1 = true; + end + function [fz1, fz2, fz3, fz4, fz5] = isOutputFixedSizeImpl(~) % Output fixed size? + fz1 = true; + fz2 = true; + fz3 = true; + fz4 = true; + fz5 = true; + end + + function setupImpl(obj) + % Perform one-time calculations, such as computing constants + import casadi.* + + %% Parameters + nState = obj.nState; + nIn = obj.nIn; + nOut = obj.nOut; + nDst = obj.nDst; + N = obj.N; + R = obj.R; + T = obj.T; + A = obj.A; + Bd = obj.Bd; + Bu = obj.Bu; + C = obj.C; + D = obj.D; + + %% Prepare variables + U = MX.sym('U', nIn, N); + P = MX.sym('P', nState + N + nDst*N); % Initial values, costElec, disturbances + X = MX.sym('X', nState, (N+1)); + S = MX.sym('S', nOut, N); % First state free + + J = 0; % Objective function + g = []; % constraints vector + + %% P indices + iX0 = [1:nState]; + iCoEl = [nState+1:nState+N]; + iDist = [nState+N+1:nState+N+nDst*N]; + + %% Disassemble P + pX0 = P(iX0); + pCoEl = P(iCoEl); + pDist = reshape(P(iDist), [nDst N]); % Prone to shaping error + + %% Define variables + states = MX.sym('states', nState); + controls = MX.sym('controls', nIn); + disturbances = MX.sym('disturbances', nDst); + + %% Dynamics + f = Function('f',{P, states, controls, disturbances},{A*states + Bu*controls + Bd*disturbances}); + + %% Compile all constraints + g = [g; X(:,1) - pX0]; + + for i = 1:N + g = [g; C*X(:,i+1) - S(:,i)]; % State/output constraints, first state free + g = [g; U(:,i)]; % Control constraints + g = [g; X(:,i+1) - f(P, X(:,i), U(:,i), pDist(:,i))]; % System dynamics + + % Cost function, first state given -> not punished + J = J + R * U(:,i) * pCoEl(i) + S(:,i)'*T*S(:,i); + end + + %% Reshape variables + OPT_variables = veccat(X, S, U); + + %% Optimization + nlp_mhe = struct('f', J, ... + 'x', OPT_variables, ... + 'g', g, ... + 'p', P); + + opts = struct; + opts.ipopt.print_level = 0; %5; + solver = nlpsol('solver', 'ipopt', nlp_mhe, opts); + + %% Pack opj + obj.casadi_solver = solver; + + end + + function [mv, xStar, sStar, uStar, yStarOut] = stepImpl(obj, xHat, dist, cE, on) + % Implement algorithm. Calculate y as a function of input u and + % discrete states. + if on > 0.5 + %% Parameters + nState = obj.nState; + N = obj.N; + nOut = obj.nOut; + nDst = obj.nDst; + nIn = obj.nIn; + yMin = obj.yMin; + yMax = obj.yMax; + uMin = obj.uMin; + uMax = obj.uMax; + C = obj.C; + solver = obj.casadi_solver; + Pdata = [xHat; cE; reshape(dist', [nDst*N, 1])]; % Prone to shaping error!!! + + %% Constraints + lbg = zeros(nState,1); % x0 constraints + ubg = zeros(nState,1); + + % Output, control and dynamics constraints + for i = 1:N + lbg = [lbg; yMin]; + lbg = [lbg; uMin]; + lbg = [lbg; zeros(nState,1)]; + + ubg = [ubg; yMax]; + ubg = [ubg; uMax]; + ubg = [ubg; zeros(nState,1)]; + end + + %% Solver + sol = solver('x0', 0, ... % x0 = x* from before, shift one time step, double last time step + 'lbg', lbg, ... + 'ubg', ubg, ... + 'p', Pdata); + + %% Outputs + xStar = reshape(full(sol.x(1 :nState*(N+1))), [nState, (N+1)])'; + sStar = reshape(full(sol.x(nState*(N+1)+1 :nState*(N+1)+nOut*N)), [nOut, N])'; + uStar = reshape(full(sol.x(nState*(N+1)+nOut*N+1:end)), [nIn, N])'; + + mv = full(sol.x(nState*(N+1)+nOut*N+1:nState*(N+1)+nOut*N+nIn))'; + yStarOut = C*xStar(2,:)'; % Second value is the target + + else % Zero output if MPC is disabled + mv = zeros(1, obj.nIn); + xStar = zeros(obj.N+1, obj.nState); + uStar = zeros(obj.N, obj.nIn); + sStar = zeros(obj.N, obj.nOut); + yStarOut = zeros(1, obj.nOut); + end % \if on + end % \stepImpl + + function resetImpl(obj) + % Initialize / reset discrete-state properties + end + end +end diff --git a/Matlab_scripts/MPCforSonja/MPCsimulink.slx b/Matlab_scripts/MPCforSonja/MPCsimulink.slx new file mode 100644 index 0000000..056e3fc Binary files /dev/null and b/Matlab_scripts/MPCforSonja/MPCsimulink.slx differ diff --git a/Matlab_scripts/MPCforSonja/setupMPCcasadi_v1_0.m b/Matlab_scripts/MPCforSonja/setupMPCcasadi_v1_0.m new file mode 100644 index 0000000..79b597e --- /dev/null +++ b/Matlab_scripts/MPCforSonja/setupMPCcasadi_v1_0.m @@ -0,0 +1,30 @@ + + +%% Settings +TimeStep = 900; % Step time +nHor = 4*24; % Length of ontrol and planning horizon +%tSmp = 0:TimeStep:nHor*TimeStep-1; + +nStt = 1; % Number of states +chY = 1; % Number of observed variables +nDst = 1; % Number of disturbance variables +nMV = 1; % Number of controlled variables + +%% System matrices +A = 1; +B = [-1, 1]/(3000*4182/TimeStep); +Bd = B(:, 1:nDst); +Bu = B(:, nDst+1:end); +C = 1; +D = 0; + +%% Constraints and normalization +uMin = 0; +uMax = 7500; +yMin = 40; +yMax = 50; + +%% Weights +R = 1/uMax/0.1; +T = 1e5*eye(chY); + diff --git a/Matlab_scripts/casadi_gp_mpc.m b/Matlab_scripts/casadi_gp_mpc.m new file mode 100644 index 0000000..14a196d --- /dev/null +++ b/Matlab_scripts/casadi_gp_mpc.m @@ -0,0 +1,86 @@ +clear all +close all +clc + +%% Import CasADi +addpath('/home/radu/Media/MATLAB/casadi-linux-matlabR2014b-v3.5.5') +import casadi.* + +load("gpr_model.mat") +%% Initialize casadi callback +cs_model = gpCallback('model'); + +T_set = 20; +N_horizon = 5; +n_states = 7; + + +COP = 5; %cooling +EER = 5; %heating +Pel = 6300; % Electric Power Consumption of the HVAC + +u_min = - COP * Pel; +u_max = EER * Pel; + +J = 0; % optimization objective +g = []; % constraints vector + +% Set up the symbolic variables +U = MX.sym('U', N_horizon, 1); +W = MX.sym('W', N_horizon, 2); +x0 = MX.sym('x0', 1, n_states - 3); + +% setup the first state +wk = W(1, :); +uk = U(1); % scaled input +xk = [wk, Pel*uk, x0]; +yk = cs_model(xk); +J = J + (yk - T_set).^2; + +% Setup the rest of the states +for idx = 2:N_horizon + wk = W(idx, :); + uk_1 = uk; uk = U(idx); + xk = [wk, Pel*uk, Pel*uk_1, yk, xk(5:6)]; + yk = cs_model(xk); + J = J + (yk - T_set).^2; +end + +p = [vec(W); vec(x0)]; + +nlp_prob = struct('f', J, 'x', vec(U), 'g', g, 'p', p); + +opts = struct; +%opts.ipopt.max_iter = 5000; +opts.ipopt.max_cpu_time = 15*60; +opts.ipopt.hessian_approximation = 'limited-memory'; +%opts.ipopt.print_level =1;%0,3 +opts.print_time = 0; +opts.ipopt.acceptable_tol =1e-8; +opts.ipopt.acceptable_obj_change_tol = 1e-6; + +solver = nlpsol('solver', 'ipopt', nlp_prob,opts); + +real_x0 = [0, 23, 23, 23]; +real_W = [[57.9261000000000;54.9020333333334;73.8607000000000;76.0425333333333;64.9819666666667], [22; 22; 22; 22; 22]]; + +real_p = vertcat(vec(DM(real_W)), vec(DM(real_x0))); + + +res = solver('p', real_p, 'ubx', EER, 'lbx', -COP); +%% Interpret the optimization result +x = Pel * full(res.x); +X = [real_W, x, [real_x0; zeros(N_horizon -1, size(real_x0, 2))]]; + +X(2:end, 4) = X(1:end-1, 3); + +for idx=2:N_horizon + X(idx, 5) = full(cs_model(X(idx - 1, :))); + X(idx, 6:7) = X(idx - 1, 5:6); +end +T_horizon = cs_model(X'); + + +figure; hold on; +plot(1:N_horizon, full(T_horizon)); +plot(1:N_horizon, T_set*ones(1, N_horizon)); \ No newline at end of file diff --git a/Matlab_scripts/gpCallback.m b/Matlab_scripts/gpCallback.m new file mode 100644 index 0000000..180b34e --- /dev/null +++ b/Matlab_scripts/gpCallback.m @@ -0,0 +1,39 @@ +classdef gpCallback < casadi.Callback + properties + model + end + methods + function self = gpCallback(name) + self@casadi.Callback(); + construct(self, name, struct('enable_fd', true)); + end + + % Number of inputs and outputs + function v=get_n_in(self) + v=1; + end + function v=get_n_out(self) + v=1; + end + % Function sparsity + function v=get_sparsity_in(self, i) + v=casadi.Sparsity.dense(7, 1); + end + + % Initialize the object + function init(self) + disp('initializing gpCallback') + gpr = load('gpr_model.mat', 'model'); + self.model = gpr.model; + end + + % Evaluate numerically + function arg = eval(self, arg) + x = full(arg{1}); + % Transpose x since gp predictor takes row by row, and casadi gives + % colum by column + [mean, ~] = predict(self.model, x'); + arg = {mean}; + end + end +end \ No newline at end of file diff --git a/Matlab_scripts/gp_casadi_test.m b/Matlab_scripts/gp_casadi_test.m new file mode 100644 index 0000000..3f122d8 --- /dev/null +++ b/Matlab_scripts/gp_casadi_test.m @@ -0,0 +1,83 @@ +clear all +close all +clc +%%%%%%%%%%%% +addpath('/home/radu/Media/MATLAB/casadi-linux-matlabR2014b-v3.5.5') +import casadi.* + +%% Generate GP data +size = 500; n_samples = 15; +X = linspace(-2, 2, size); +Y = 3 * X .^2; + +% Add noise to the output +mean = 0; std = 0.5; +noise = mean + std.*randn(1, n_samples); + +idx_samples = randperm(size, n_samples); + +X_sampled = X(idx_samples); +Y_sampled = Y(idx_samples); + +Y_sampled = Y_sampled + noise; + +figure; hold on; +plot(X, Y); +scatter(X_sampled, Y_sampled); + +tbl_gpr_in = array2table([X_sampled', Y_sampled']); +tbl_gpr_in.Properties.VariableNames = {'X', 'Y'}; + +tic; +model = fitrgp(tbl_gpr_in, 'Y', 'KernelFunction', 'ardsquaredexponential', ... + 'FitMethod', 'sr', 'PredictMethod', 'fic', 'Standardize', 1); +toc; + +%% Predict stuff +[yhat_test, sigma_test] = predict(model, X'); +std_test = sqrt(sigma_test); + +% prepare it for the fill function +x_ax = X'; +X_plot = [x_ax; flip(x_ax)]; +Y_plot = [yhat_test-1.96.*std_test; flip(yhat_test+1.96.*std_test)]; + +% plot a line + confidence bands +figure(); hold on; +title("GP performance on test data"); +plot(x_ax, Y, 'red', 'LineWidth', 1.2); +plot(x_ax, yhat_test, 'blue', 'LineWidth', 1.2) +fill(X_plot, Y_plot , 1,.... + 'facecolor','blue', ... + 'edgecolor','none', ... + 'facealpha', 0.3); +legend({'data','prediction_mean', '95% confidence'},'Location','Best'); +hold off + +%% Save the model +save('test_gpr_model.mat', 'model') + +%% CasADi optimization problem +cs_model = test_gpCallback('model'); +cs_x = MX.sym('x'); +cs_y = 2 * cs_model(cs_x) + 5; +f = Function('f', {cs_x}, {cs_y}); + + +nlp_prob = struct('f', f(cs_x), 'x', cs_x); + + +opts = struct; +opts.ipopt.max_iter = 2000; +opts.ipopt.hessian_approximation = 'limited-memory'; +%opts.ipopt.print_level =1;%0,3 +opts.print_time = 0; +opts.ipopt.acceptable_tol =1e-8; +opts.ipopt.acceptable_obj_change_tol = 1e-6; + +solver = nlpsol('solver', 'ipopt', nlp_prob,opts); + +res = solver('lbx', -2, 'ubx', 2); + +res + diff --git a/Matlab_scripts/gp_identification.m b/Matlab_scripts/gp_identification.m new file mode 100644 index 0000000..68bb67e --- /dev/null +++ b/Matlab_scripts/gp_identification.m @@ -0,0 +1,71 @@ +clear all +close all +clc +%%%%%%%%%%% + +load("gpr_carnot.mat"); + +%% Format the train/test data arrays +tbl_gpr_train = array2table(gpr_train); +tbl_gpr_train.Properties.VariableNames = cellstr(table_cols); +tbl_gpr_train = removevars(tbl_gpr_train,{'u'}); +tbl_gpr_train_x = removevars(tbl_gpr_train, {'y'}); + +tbl_gpr_test = array2table(gpr_test); +tbl_gpr_test.Properties.VariableNames = cellstr(table_cols); +tbl_gpr_test = removevars(tbl_gpr_test,{'u'}); +tbl_gpr_test_x = removevars(tbl_gpr_test, {'y'}); + + + +%% Train the GP model +OutputName = 'y'; + +tic; +model = fitrgp(tbl_gpr_train, OutputName, 'KernelFunction', 'ardsquaredexponential', ... + 'FitMethod', 'sr', 'PredictMethod', 'fic', 'Standardize', 1); +toc; +%% Validate the model using training data +[yhat_train, sigma_train] = predict(model, tbl_gpr_train_x); +std_train = sqrt(sigma_train); + +% prepare it for the fill function +x_ax = (1:size(tbl_gpr_train, 1))'; +X_plot = [x_ax; flip(x_ax)]; +Y_plot = [yhat_train-1.96.*std_train; flip(yhat_train+1.96.*std_train)]; + +% plot a line + confidence bands +figure(); hold on; +title("GP performance on training data"); +plot(x_ax, tbl_gpr_train.y, 'red', 'LineWidth', 1.2); +plot(x_ax, yhat_train, 'blue', 'LineWidth', 1.2) +fill(X_plot, Y_plot , 1,.... + 'facecolor','blue', ... + 'edgecolor','none', ... + 'facealpha', 0.3); +legend({'data','prediction_mean', '95% confidence'},'Location','Best'); +hold off + +%% Validate the model using test data +[yhat_test, sigma_test] = predict(model, tbl_gpr_test_x); +std_test = sqrt(sigma_test); + +% prepare it for the fill function +x_ax = (1:size(tbl_gpr_test, 1))'; +X_plot = [x_ax; flip(x_ax)]; +Y_plot = [yhat_test-1.96.*std_test; flip(yhat_test+1.96.*std_test)]; + +% plot a line + confidence bands +figure(); hold on; +title("GP performance on test data"); +plot(x_ax, tbl_gpr_test.y, 'red', 'LineWidth', 1.2); +plot(x_ax, yhat_test, 'blue', 'LineWidth', 1.2) +fill(X_plot, Y_plot , 1,.... + 'facecolor','blue', ... + 'edgecolor','none', ... + 'facealpha', 0.3); +legend({'data','prediction_mean', '95% confidence'},'Location','Best'); +hold off + +%% Export the final GP model +save('gpr_model.mat', 'model') \ No newline at end of file diff --git a/Matlab_scripts/gpr_carnot.mat b/Matlab_scripts/gpr_carnot.mat new file mode 100644 index 0000000..ff67d5c Binary files /dev/null and b/Matlab_scripts/gpr_carnot.mat differ diff --git a/Matlab_scripts/gpr_model.mat b/Matlab_scripts/gpr_model.mat new file mode 100644 index 0000000..3f9319e Binary files /dev/null and b/Matlab_scripts/gpr_model.mat differ diff --git a/Matlab_scripts/test_gpCallback.m b/Matlab_scripts/test_gpCallback.m new file mode 100644 index 0000000..443786e --- /dev/null +++ b/Matlab_scripts/test_gpCallback.m @@ -0,0 +1,39 @@ +classdef test_gpCallback < casadi.Callback + properties + model + end + methods + function self = test_gpCallback(name) + self@casadi.Callback(); + construct(self, name, struct('enable_fd', true)); + end + + % Number of inputs and outputs + function v=get_n_in(self) + v=1; + end + function v=get_n_out(self) + v=1; + end + % Function sparsity + function v=get_sparsity_in(self, i) + v=casadi.Sparsity.dense(1, 1); + end + + % Initialize the object + function init(self) + disp('initializing gpCallback') + gpr = load('test_gpr_model.mat', 'model'); + self.model = gpr.model; + end + + % Evaluate numerically + function arg = eval(self, arg) + x = full(arg{1}); + % Transpose x since gp predictor takes row inputs, and casadi gives + % colum by column + [mean, ~] = predict(self.model, x'); + arg = {mean}; + end + end +end \ No newline at end of file diff --git a/Matlab_scripts/test_gpr_model.mat b/Matlab_scripts/test_gpr_model.mat new file mode 100644 index 0000000..af245e5 Binary files /dev/null and b/Matlab_scripts/test_gpr_model.mat differ diff --git a/Simulink/Exp_CARNOT.mat b/Simulink/Exp_CARNOT.mat new file mode 100644 index 0000000..530d6ba Binary files /dev/null and b/Simulink/Exp_CARNOT.mat differ diff --git a/Simulink/Exp_datasets_simid.sid b/Simulink/Exp_datasets_simid.sid new file mode 100644 index 0000000..e4600e0 Binary files /dev/null and b/Simulink/Exp_datasets_simid.sid differ diff --git a/Simulink/Images/Exp1_simulation (1).svg b/Simulink/Images/Exp1_simulation (1).svg new file mode 100644 index 0000000..4d604d5 --- /dev/null +++ b/Simulink/Images/Exp1_simulation (1).svg @@ -0,0 +1,242 @@ + + + diff --git a/Simulink/gp_mpc_system.m b/Simulink/gp_mpc_system.m index 56a7a7b..0a0148a 100644 --- a/Simulink/gp_mpc_system.m +++ b/Simulink/gp_mpc_system.m @@ -133,14 +133,17 @@ classdef gp_mpc_system < matlab.System & matlab.system.mixin.Propagates real_p = vertcat(vec(DM(w)), vec(DM([obj.u_lags obj.y_lags]))); disp("Starting optimization") tic - %res = obj.casadi_solver('p', real_p, 'ubx', obj.ubx, 'lbx', obj.lbx); + res = obj.casadi_solver('p', real_p, 'ubx', obj.ubx, 'lbx', obj.lbx); t = toc; disp(t) u = obj.Pel * full(res.x(1)); + u = 15000 * (20 - x); % Update the u lags obj.u_lags = [u, obj.u_lags(2:end-1)]; + + end function resetImpl(obj) diff --git a/Simulink/model_identification.m b/Simulink/model_identification.m index ff04de0..0be1205 100644 --- a/Simulink/model_identification.m +++ b/Simulink/model_identification.m @@ -3,74 +3,58 @@ 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 run parameters -air_exchange_rate = tin; -air_exchange_rate(:,2) = 1.0; - % Set the initial temperature to be the measured initial temperature -t0 = Exp_data.(exp_id).InsideTemp.values(1); +t0 = 23; -power = Exp_data.(exp_id).Power.values - 1.67 * 1000; +runtime1 = 161400; +runtime2 = 136200; +runtime3 = 208200; +runtime4 = 208200; +runtime5 = 208200; +runtime6 = 208200; +runtime7 = 553800; -power = [tin Heating_coeff .* power]; +runtime = 24 * 3600; +set_param('polydome', 'StopTime', int2str(runtime)) +Tsample = 900; +steps = runtime/Tsample; +tin = Tsample *(0:steps)'; + +prbs_sig = 2*prbs(8, steps+1)' - 1; +COP = 5.0; +Pel = 6300; -% Turn down the air exchange rate when the HVAC is not running -night_air_exchange_rate = 0.5; -air_exchange_rate(abs(power(:, 2)) < 100, 2) = night_air_exchange_rate; +power = [tin COP*Pel*prbs_sig(1:steps+1)]; -%% 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"); +%% Simulate the model +out = sim('polydome'); -SimulatedTemp = simout.SimulatedTemp; -%% Compare the simulation results with the measured values -figure; hold on; grid minor; -plot(tin, InsideTemp); -plot(tin, OutsideTemp); -plot(SimulatedTemp, 'LineWidth', 2); -legend('InsideTemp', 'OutsideTemp', 'SimulatedTemp'); +%% For manual simulation running +WeatherMeasurement = struct; +WeatherMeasurement.data = squeeze(out.WeatherMeasurement.data)'; +WeatherMeasurement.time = out.WeatherMeasurement.time; +input = [power(:, 2:end) WeatherMeasurement.data]; -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)); +Exp7_data = iddata(out.SimulatedTemp.data, input); -hold off; +Exp7_table = array2table([input out.SimulatedTemp.data], 'VariableNames', {'Power', 'SolRad', 'OutsideTemp', 'SimulatedTemp'}); -saveas(gcf, strcat(exp_id, '_simulation'), 'svg') +writetable(Exp7_table, 'Exp7_table.csv') -%% Export simulated temperature to a .mat file for further use -carnot_output_dir = strcat("../Data/CARNOT_output/",exp_id,"_carnot_temp.mat"); -save(carnot_output_dir, 'SimulatedTemp'); +%% +save('Exp_CARNOT.mat', ... + 'Exp1_data', 'Exp1_table', ... + 'Exp2_data', 'Exp2_table', ... + 'Exp3_data', 'Exp3_table', ... + 'Exp4_data', 'Exp4_table', ... + 'Exp5_data', 'Exp5_table', ... + 'Exp6_data', 'Exp6_table', ... + 'Exp7_data', 'Exp7_table' ... +) + +data_train = merge(Exp1_data, Exp3_data, Exp5_data); +data_test = merge(Exp2_data, Exp4_data, Exp6_data, Exp7_data); diff --git a/Simulink/mpc_simulink/casadi_block.m b/Simulink/mpc_simulink/casadi_block.m deleted file mode 100644 index 9bebd96..0000000 --- a/Simulink/mpc_simulink/casadi_block.m +++ /dev/null @@ -1,177 +0,0 @@ -classdef casadi_block < matlab.System & matlab.system.mixin.Propagates - % untitled Add summary here - % - % This template includes the minimum set of functions required - % to define a System object with discrete state. - - properties - % Public, tunable properties. - - end - - properties (DiscreteState) - end - - properties (Access = private) - % Pre-computed constants. - casadi_solver - x0 - lbx - ubx - lbg - ubg - end - - methods (Access = protected) - function num = getNumInputsImpl(~) - num = 2; - end - function num = getNumOutputsImpl(~) - num = 1; - end - function dt1 = getOutputDataTypeImpl(~) - dt1 = 'double'; - end - function dt1 = getInputDataTypeImpl(~) - dt1 = 'double'; - end - function sz1 = getOutputSizeImpl(~) - sz1 = [1,1]; - end - function sz1 = getInputSizeImpl(~) - sz1 = [1,1]; - end - function cp1 = isInputComplexImpl(~) - cp1 = false; - end - function cp1 = isOutputComplexImpl(~) - cp1 = false; - end - function fz1 = isInputFixedSizeImpl(~) - fz1 = true; - end - function fz1 = isOutputFixedSizeImpl(~) - fz1 = true; - end - function setupImpl(obj,~,~) - % Implement tasks that need to be performed only once, - % such as pre-computed constants. - - import casadi.* - - T = 10; % Time horizon - N = 20; % number of control intervals - - % Declare model variables - x1 = SX.sym('x1'); - x2 = SX.sym('x2'); - x = [x1; x2]; - u = SX.sym('u'); - - % Model equations - xdot = [(1-x2^2)*x1 - x2 + u; x1]; - - % Objective term - L = x1^2 + x2^2 + u^2; - - % Continuous time dynamics - f = casadi.Function('f', {x, u}, {xdot, L}); - - % Formulate discrete time dynamics - % Fixed step Runge-Kutta 4 integrator - M = 4; % RK4 steps per interval - DT = T/N/M; - f = Function('f', {x, u}, {xdot, L}); - X0 = MX.sym('X0', 2); - U = MX.sym('U'); - X = X0; - Q = 0; - for j=1:M - [k1, k1_q] = f(X, U); - [k2, k2_q] = f(X + DT/2 * k1, U); - [k3, k3_q] = f(X + DT/2 * k2, U); - [k4, k4_q] = f(X + DT * k3, U); - X=X+DT/6*(k1 +2*k2 +2*k3 +k4); - Q = Q + DT/6*(k1_q + 2*k2_q + 2*k3_q + k4_q); - end - F = Function('F', {X0, U}, {X, Q}, {'x0','p'}, {'xf', 'qf'}); - - % Start with an empty NLP - w={}; - w0 = []; - lbw = []; - ubw = []; - J = 0; - g={}; - lbg = []; - ubg = []; - - % "Lift" initial conditions - X0 = MX.sym('X0', 2); - w = {w{:}, X0}; - lbw = [lbw; 0; 1]; - ubw = [ubw; 0; 1]; - w0 = [w0; 0; 1]; - - % Formulate the NLP - Xk = X0; - for k=0:N-1 - % New NLP variable for the control - Uk = MX.sym(['U_' num2str(k)]); - w = {w{:}, Uk}; - lbw = [lbw; -1]; - ubw = [ubw; 1]; - w0 = [w0; 0]; - - % Integrate till the end of the interval - Fk = F('x0', Xk, 'p', Uk); - Xk_end = Fk.xf; - J=J+Fk.qf; - - % New NLP variable for state at end of interval - Xk = MX.sym(['X_' num2str(k+1)], 2); - w = {w{:}, Xk}; - lbw = [lbw; -0.25; -inf]; - ubw = [ubw; inf; inf]; - w0 = [w0; 0; 0]; - - % Add equality constraint - g = {g{:}, Xk_end-Xk}; - lbg = [lbg; 0; 0]; - ubg = [ubg; 0; 0]; - end - - % Create an NLP solver - prob = struct('f', J, 'x', vertcat(w{:}), 'g', vertcat(g{:})); - options = struct('ipopt',struct('print_level',0),'print_time',false); - solver = nlpsol('solver', 'ipopt', prob, options); - - obj.casadi_solver = solver; - obj.x0 = w0; - obj.lbx = lbw; - obj.ubx = ubw; - obj.lbg = lbg; - obj.ubg = ubg; - end - - function u = stepImpl(obj,x,t) - disp(t) - tic - w0 = obj.x0; - lbw = obj.lbx; - ubw = obj.ubx; - solver = obj.casadi_solver; - lbw(1:2) = x; - ubw(1:2) = x; - sol = solver('x0', w0, 'lbx', lbw, 'ubx', ubw,... - 'lbg', obj.lbg, 'ubg', obj.ubg); - - u = full(sol.x(3)); - toc - end - - function resetImpl(obj) - % Initialize discrete-state properties. - end - end -end diff --git a/Simulink/mpc_simulink/mpc_demo.slx b/Simulink/mpc_simulink/mpc_demo.slx deleted file mode 100644 index 6d54e37..0000000 Binary files a/Simulink/mpc_simulink/mpc_demo.slx and /dev/null differ diff --git a/Simulink/mpc_simulink/mpc_demo.slx.r2014b b/Simulink/mpc_simulink/mpc_demo.slx.r2014b deleted file mode 100644 index 3827f32..0000000 Binary files a/Simulink/mpc_simulink/mpc_demo.slx.r2014b and /dev/null differ diff --git a/Simulink/polydome.slx b/Simulink/polydome.slx index 75ae951..69713eb 100644 Binary files a/Simulink/polydome.slx and b/Simulink/polydome.slx differ diff --git a/Simulink/polydome.slx.autosave b/Simulink/polydome.slx.autosave deleted file mode 100644 index e70cc7b..0000000 Binary files a/Simulink/polydome.slx.autosave and /dev/null differ diff --git a/Simulink/polydome_params.mat b/Simulink/polydome_params.mat index 242bdf3..ab9b739 100644 Binary files a/Simulink/polydome_params.mat and b/Simulink/polydome_params.mat differ diff --git a/Simulink/polydome_python.slx b/Simulink/polydome_python.slx new file mode 100644 index 0000000..39d6597 Binary files /dev/null and b/Simulink/polydome_python.slx differ diff --git a/Simulink/test_casadi_mpc.m b/Simulink/test_casadi_mpc.m deleted file mode 100644 index 2e0a7da..0000000 --- a/Simulink/test_casadi_mpc.m +++ /dev/null @@ -1,10 +0,0 @@ -clear all -close all -clc -%%%%%%%%%%%%%%%% - - -%% Load the existing GP -addpath("../../Gaussiandome/Identification/Computation results/") -load("Identification_Validation.mat") -load("Gaussian_Process_models.mat") \ No newline at end of file diff --git a/Simulink/weather_predictor.m b/Simulink/weather_predictor.m index 822934d..34e424c 100644 --- a/Simulink/weather_predictor.m +++ b/Simulink/weather_predictor.m @@ -29,17 +29,20 @@ classdef weather_predictor < matlab.System num = 2; end function num = getNumOutputsImpl(~) - num = 1; + num = 2; end - function dt1 = getOutputDataTypeImpl(~) + function [dt1, dt2] = getOutputDataTypeImpl(~) dt1 = 'double'; + dt2 = 'double'; + end function [dt1, dt2] = getInputDataTypeImpl(~) dt1 = 'double'; dt2 = 'double'; end - function sz1 = getOutputSizeImpl(obj) - sz1 = [obj.N 2]; + function [sz1, sz2] = getOutputSizeImpl(obj) + sz1 = [1 2]; + sz2 = [obj.N 2]; end function sz1 = getInputSizeImpl(~) sz1 = 1; @@ -47,14 +50,16 @@ classdef weather_predictor < matlab.System function cp1 = isInputComplexImpl(~) cp1 = false; end - function cp1 = isOutputComplexImpl(~) + function [cp1, cp2] = isOutputComplexImpl(~) cp1 = false; + cp2 = false; end function fz1 = isInputFixedSizeImpl(~) - fz1 = true; + fz1 = true; end - function fz1 = isOutputFixedSizeImpl(~) + function [fz1, fz2] = isOutputFixedSizeImpl(~) fz1 = true; + fz2 = true; end @@ -63,13 +68,29 @@ classdef weather_predictor < matlab.System % Perform one-time calculations, such as computing constants end - function w = stepImpl(obj,wdb_mat,timestamp) + function [w, w_hat] = stepImpl(obj,wdb_mat,timestamp) disp(timestamp) - % Implement algorithm. Calculate y as a function of input u and - % discrete states. - curr_idx = find(wdb_mat(:, 1) == timestamp); - N_idx = (1:obj.N) + curr_idx; - w = [wdb_mat(N_idx, 18) + wdb_mat(N_idx, 19), wdb_mat(N_idx, 7)]; + + + forecast_start = timestamp + obj.TimeStep; + forecast_end = timestamp + obj.N * obj.TimeStep; + + xq = forecast_start:obj.TimeStep:forecast_end; + + weather_start_idx = find(wdb_mat(:, 1) <= timestamp, 1); + weather_end_idx = find(wdb_mat(:, 1) >= forecast_end, 1); + weather_idx = weather_start_idx:weather_end_idx; + + solar_direct = interp1(wdb_mat(weather_idx, 1), wdb_mat(weather_idx, 18), timestamp); + solar_diffuse = interp1(wdb_mat(weather_idx, 1), wdb_mat(weather_idx, 19), timestamp); + outside_temp = interp1(wdb_mat(weather_idx, 1), wdb_mat(weather_idx, 7), timestamp); + w = [solar_direct + solar_diffuse, outside_temp]; + + solar_direct = interp1(wdb_mat(weather_idx, 1), wdb_mat(weather_idx, 18), xq)'; + solar_diffuse = interp1(wdb_mat(weather_idx, 1), wdb_mat(weather_idx, 19), xq)'; + outside_temp = interp1(wdb_mat(weather_idx, 1), wdb_mat(weather_idx, 7), xq)'; + + w_hat = [solar_direct + solar_diffuse, outside_temp]; end function resetImpl(obj) diff --git a/Simulink/weather_predictor2.m b/Simulink/weather_predictor2.m deleted file mode 100644 index f890417..0000000 --- a/Simulink/weather_predictor2.m +++ /dev/null @@ -1,8 +0,0 @@ -function w = weather_predictor2(wdb_mat,timestamp, N) -%WEATHER_PREDICTOR2 Summary of this function goes here -% Detailed explanation goes here -curr_idx = find(wdb_mat(:, 1) == timestamp); -N_idx = (1:N) + curr_idx; -w = [wdb_mat(N_idx, 18) + wdb_mat(N_idx, 19), wdb_mat(N_idx, 7)]; - -end \ No newline at end of file diff --git a/carnot_vars_nomenclature b/carnot_vars_nomenclature new file mode 100644 index 0000000..bfefe8e --- /dev/null +++ b/carnot_vars_nomenclature @@ -0,0 +1,42 @@ +Model Parameters +================ + +Side Windows +------------ +windows_size - [height width] of the large side-windows +window_roof_size - [height width] of the roof windows (approximated as one large window as opposed to + the four smaller windows in + reality) +surface_part - ??? +U - heat transfer coefficient [W/m2K] (of the window?) +g - total solar energy transmittance +v_g = 0.65 - transmittance in visible range of the sunlight + + +Roof +---- +wall_size - [length width] of the roof (considered flat) +roof_position - position angles of the roof + + +Roof structure +~~~~~~~~~~~~~~ +The roof is supposed to be made of + - 5cm wood + - 10cm insulation + - 5cm wood + +Data from https://simulationresearch.lbl.gov/modelica/releases/latest/help/Buildings_HeatTransfer_Data_Solids.html#Buildings.HeatTransfer.Data.Solids.Plywood + +node_thickness +node_conductivity +node_capacity +node_density + +Floor +----- +ceiling_size +layer_thickness +layer_conductivity +layer_capacity +layer_density% Large side windows diff --git a/research-notes b/research-notes index 52263c1..c033119 100644 --- a/research-notes +++ b/research-notes @@ -6,3 +6,23 @@ Info on windows (transmittance, U-factor, etc.) Computing DNI/DHI from GHI and solar position: https://plantpredict.com/algorithm/irradiance-radiation/ + + +Gaussian Process Kernel Cookbook: + https://www.cs.toronto.edu/~duvenaud/cookbook/ + +GPflow Manipulating Kernels: + https://gpflow.readthedocs.io/en/develop/notebooks/advanced/kernels.html + +Neural Network Gaussian Process: + https://en.wikipedia.org/wiki/Neural_network_Gaussian_process + +Sparse Gaussian Process Tutorial (github): + http://krasserm.github.io/2018/03/19/gaussian-processes/ + http://krasserm.github.io/2020/12/12/gaussian-processes-sparse/ + +Sparse Gaussian Processes (youtube): + https://www.youtube.com/watch?v=sQmsQq_Jfi8 + +Deep Gaussian Processes (youtube): + https://www.youtube.com/watch?v=750fRY9-uq8