How can I add a time controlled variable in an ODE system?

Question

Isra Haroun 2018-10-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/422658-how-can-i-add-a-time-controlled-variable-in-an-ode-system

评论： Isra Haroun 2018-10-10

Hello

I am trying to convert the model in this link into a matlab script. I made the system a 3D model for dTroom/dt, dQlosses/dt, and dQheate/dt as the model suggests. Here is the code for the 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

and here is the code to call the system:

        clear;
        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');

The output I get is as shown in the figure bellow:

Now,I know that this is not how the output is supposed to look like, it could be due to the sine wave but I did it exactly like what they did in the model. Also, adding the status parameter was a bit tricky for me, and I don't know if what I did is correct or not. Can someone help me in this, I couldn't get status out and plot it to verify if the system works correctly,(the logic says, if Troom was increasing and within the limit then status must be on, if Troom exceeds the limit them definately status must be 0).

Can someone PLEASE help me in this matter?