Я пытаюсь преобразовать модель по этой ссылке в скрипт MATLAB.Я сделал систему 3D-моделью для dTroom / dt, dQlosses / dt и dQheate / dt, как предлагает модель.Вот код для ODE:
function dX_dt = odes_thermal2(t ,X, Tout_t)
% TinIC = 26;
r2d = 180/pi;
Troom=X(1);
Qlosses=X(2);
Qheater=X(3);
% stat=X(4);
Tout=Tout_t(t);
stat(1)=0;
% -------------------------------
% Define the house geometry
% -------------------------------
% House length = 30 m
lenHouse = 30;
% House width = 10 m
widHouse = 10;
% House height = 4 m
htHouse = 4;
% Roof pitch = 40 deg
pitRoof = 40/r2d;
% Number of windows = 6
numWindows = 6;
% Height of windows = 1 m
htWindows = 1;
% Width of windows = 1 m
widWindows = 1;
windowArea = numWindows*htWindows*widWindows;
wallArea = 2*lenHouse*htHouse + 2*widHouse*htHouse + ...
2*(1/cos(pitRoof/2))*widHouse*lenHouse + ...
tan(pitRoof)*widHouse - windowArea;
% -------------------------------
% Define the type of insulation used
% -------------------------------
% Glass wool in the walls, 0.2 m thick
% k is in units of J/sec/m/C - convert to J/hr/m/C multiplying by 3600
kWall = 0.038*3600; % hour is the time unit
LWall = .2;
RWall = LWall/(kWall*wallArea);
% Glass windows, 0.01 m thick
kWindow = 0.78*3600; % hour is the time unit
LWindow = .01;
RWindow = LWindow/(kWindow*windowArea);
% -------------------------------
% Determine the equivalent thermal resistance for the whole building
% -------------------------------
Req = RWall*RWindow/(RWall + RWindow);
% c = cp of air (273 K) = 1005.4 J/kg-K
c = 1005.4;
% -------------------------------
% Enter the temperature of the heated air
% -------------------------------
% The air exiting the heater has a constant temperature which is a heater
% property. THeater =20 deg C
THeater = 20;
% Air flow rate Mdot = 1 kg/sec = 3600 kg/hr
Mdot = 3600; % hour is the time unit
% -------------------------------
% Determine total internal air mass = M
% -------------------------------
% Density of air at sea level = 1.2250 kg/m^3
densAir = 1.2250;
M = (lenHouse*widHouse*htHouse+tan(pitRoof)*widHouse*lenHouse)*densAir;
%%
t1=22;t2=26;
% if(Troom(t-1)-Troom(t-2)>0&&Troom(t-1)>t1 && Troom(t-1)<t2)
% status(t)=1;
% end
% if(Troom(t-1)-Troom(t-2)<=0&&Troom(t-1)>t1 && Troom(t-1)<t2)
% status(t)=0;
% end
% if(Troom(t-1)>t2)
% status(t)=0;
% end
% if(Troom(t-1)<t1)
% status(t)=1;
% end
% status=0;
%%
dQheater_dt= (THeater-Troom)*Mdot*c*stat;
dQlosses_dt=(Troom-Tout)/Req;
dTroom_dt=(1/(M*c))*( dQheater_dt - dQlosses_dt );
if dTroom_dt >=0
stat=1;
end
if dTroom_dt<=0
stat=0;
end
if Troom >t2
stat=0;
end
if Troom<t1
stat =1;
end
dX_dt=[dTroom_dt;dQlosses_dt;dQheater_dt];
end
, а вот код для вызова системы:
t=1:1:48;
T_const = 40; k = 2; Tout = @(t) T_const + k * sin(t);
Troom_0=25;Qlosses_0=0;Qheater_0=0;
X0=[Troom_0;Qlosses_0;Qheater_0];
tspan=[0,48];
[t,X] = ode45(@(t,y) odes_thermal3(t,y,Tout),tspan,X0);
Troom=X(:,1);
Qlosses=X(:,2);
Qheater=X(:,3);
% stat=X(:,4);
figure(1);
plot(t,Troom);
ylabel('Troom');
xlabel('t');
Вывод, который я получаю, показан на рисунке ниже:
Теперь я знаю, что это не то, как должен выглядеть выходной сигнал, это может быть из-за синусоидальной волны, но я сделал это в точности какчто они сделали в модели.Кроме того, добавление параметра status было немного сложным для меня, и я не знаю, правильно ли я сделал то, что сделал.Я не смог получить статус и нарисовать его, чтобы проверить, правильно ли работает система (логика говорит, что если Troom возрастал и в пределах лимита, то статус должен быть включен, если Troom превышает лимит, он определенно должен статусбыть 0). Как получить вектор состояния из функции ODE?