Я работаю над симуляцией переменного тока так же, как показано в этой ссылке .Мне удалось преобразовать его в переменный ток, перевернув некоторые уравнения и повторно инициализировав некоторые значения.Теперь мне нужно сделать это в Matlab, чтобы моя система работала плавно и быстро, без необходимости запускать simulink для каждой итерации.Я прочитал ODE45 и решил использовать его, во-первых, я преобразовал модель на веб-странице в функцию следующим образом:
function dTroom_dt = odes_thermal(Troom,Tout)
TinIC = 26;
r2d = 180/pi;
% -------------------------------
% 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;
%%
dQlosses_dt= (THeater-Troom)*Mdot*c;
dQlosses_dt=(Troom-Tout)/Req;
dTroom_dt=(1/(M*c))*(dQlosses_dt);
end
Затем я просто вызвал ее, используя следующие строки:
clear;
t=1:1:1440;
Troom_0=25;
[t,Troom] = ode45('odes_thermal',[t(1) t(3)],Troom_0);
figure(1);
plot(t,Troom);
когда я запускаю его, вот что я получаю: 
Мой вопрос заключается в том, как включить в него состояние или отключить функции.В Simulink они использовали реле, которое будет выдавать сигнал состояния, который будет 1 или 0 в зависимости от температуры, вот как они это сделали:

Я знаю, что это может быть представлено как "dQlosses_dt = (THeater-Troom) Mdot c * status;"но я не знаю, как обновить состояние для каждого временного образца, так как нет возможности проверить Troom в каждом временном интервале. Можете ли вы мне помочь в этом или хотя бы с чего начать?
Редактировать:Следуя комментарию Pacta, я попытался превратить его в многовариантную систему с двумя выходами Troom и Qlosses.Вот новые функции, но я получаю сообщение об ошибке (слишком много входных аргументов в odes_thermal2)
clear;
t=1:1:1440;
Troom_0=25;Qlosses_0=10;
X0=[Troom_0;Qlosses_0];
[t,X] = ode45('odes_thermal2',[t(1) t(end)],X0);
Troom=X(1);
Qlosses=X(2);
figure(1);
plot(t,Troom);
ylabel('Troom');
xlabel('t');
И функция ODE имеет следующие изменения:
function dX_dt = odes_thermal2(X)
TinIC = 26;
r2d = 180/pi;
и
X(2)= (THeater-Troom)*Mdot*c;
dQlosses_dt=(Troom-Tout)/Req;
X(1)=(1/(M*c))*(dQlosses_dt);
Спасибо