Code covered by the BSD License

# Numerical Computing with Simulink, Vol. 1

08 Sep 2007 (Updated 19 May 2009)

This sequel to Numerical Computing with MATLAB explores the mathematics of simulation.

TwoRoomHouse24.m
```% This m-file generates the state space model for the 2-Room house heating
% system model in Chapter 6 of NCS.

%  The Thermal Cpacitance Data:

CR1     = 62.67;
CR2     = 83.03;
CHE1    = 62.5;
CHE2    = 62.5;

% The Thermal Resistance Data:

RHE1    = 0.01*3600;
RHE2    = 0.01*3600;
RHE3    = 0.01*3600;
Rwall1  = 0.0264*3600;
Rflceil1= 0.1*3600;
Rwind1  = 0.0417*3600;
Rwall2  = 0.0242*3600;
Rflceil2= 0.075*3600;
Rwind2  = 0.0167*3600;
RL1     = 0.1*3600;
RL2     = 0.5*3600;
RL3     = 0.5*3600;

% The Values in the Model:

Req2    = 1/(1/Rwall1+1/Rflceil1+1/Rwind1);
Req1    = 1/(1/Req2+1/RHE1);
Req4    = 1/(1/Rwall2+1/Rflceil2+1/Rwind2);
Req3    = 1/(1/Req4+1/RHE2);
Req5    = 1/(1/RL1+1/RL2+1/RL3);
alpha1  = 1-Req5/RL2;
alpha2  = 1-Req5/RL3;

% The C Matrix:
C   = diag([CR1 CR2 CHE1 CHE2]);
G   = [-1/Req1    0             1/RHE1               0
0    -1/Req3            0                1/RHE2
1/RHE1    0     -(1/RHE1+alpha1/RL2)         0
0     1/RHE2            0         -(1/RHE2+alpha2/RL3)];
A   = C\G;

%  Heat Lost to the Outside Ambient Temperature:
B1  = C\[1/Req2 1/Req4       0         0   ]';
%  Heat Flow into Zone for Room 1:
B2  = C\[  0      0     Req5/RL2       0   ]';
%  Heat Flow into Zone for Room 2:
B3  = C\[  0      0          0    Req5/RL3 ]';

%  Sources in the model:

Eff         = 0.8;                  % Boiler Efficiency
HperGal     = 100000;               % BTU per gallon for Oil
NozzleFlow  = 0.7;                   % Nozzle FLow Rate (gal./hr.)
FlowRate    = NozzleFlow/3600;      % Flow rate (Gal/Sec.)
Qf          = Eff*HperGal*FlowRate; % Maximum Heat flow from the furnace

% Linear Control System PID Gains and Sample Time:

delt         = 5;            % Sample time of the controller
Control_Gain = 1;            % PID Proportional Control Gain
Ki           = 6.25e-6;         % PID Integral Gain
Kd           = 1000;            % PID derivative Gain

% Simulation Data

Tend = 24;                   % Simulationn stop time in Hours
Tend = 3600*Tend;            %  --- Converted to Seconds.

%  Create 24 hour outside temperature profile
Weather_data;

Tend1 = 24*3600;            % Create a stop time for the 24 hour data set. ```