if temperature 17 increase until 23, swich on,Q=15*10^6 and if temperature 23 decrease untill17 switch off, Q = 0 how can I plot the temperature with time when T(0)=20degree

dT/dt =(Q/15172.5)+29.66+2.966T

[Sam Chak]

Well, @annop butmarn. Your intention is just to plot the temperature in your description, and you want to know the command in MATLAB.

Please clarify that this is not a temperature control design problem and it has nothing to do with the differential equation with the initial condition C, such that

Are you asking how to model the hysteresis?

[annop butmarn]

first of all, I'm sorry for my english skill. And I want my graph to look like red line (for example) .

can you help me what command I have to use in matlab?

[annop butmarn]

% I try this but it doesn't look like what I want.
clear; clc;
syms s t y(t) Y Q Tmax Tmin;
Tmax = 23
Tmin = 17
eqn1 = diff(y)==(Q/15172.5)+29.66-2.966*y ;
lteqn1 = laplace(eqn1,t,s) ;
neweqn1 = subs(lteqn1,[y(0),laplace(y(t),t,s)],[20 Y]) ;
ytrans1 = simplify(solve(neweqn1,Y)) ;
y = ilaplace(ytrans1,s,t) ;

for y =[Tmin,Tmax]
    if y<=Tmin
      e1 = subs(y,Q,0)
fplot(y)
    elseif y>=Tmax
       e2= subs(y,Q,(1.5*10^6))
    end
fplot(y,[0 20])
axis([0,100,5,25])
title 'room temperature corresponding to the heater status'
xlabel 'time'
ylabel 'temperature'
grid on
end

[Sam Chak]

Hi @annop butmarn

There is nothing wrong with your English. I can understand what you wrote. It was just I don't truly get what you want at the beginning, because plotting the results requires solving the temperature control differential equation, and modeling the hysteresis behavior in the thermostat (Q). Yet, you did not mention about any issue with the ODE. That's why I requested clarification from you.

Luckily you cleared things up. However, I have discovered two things:

1. Since the thermostat exhibit nonlinearity behavior (from the hysteresis), the Laplace and Inverse Laplace Transforms cannot be used in the nonlinear system.
2. Your modeling of the thermostat (Q) appears to be incorrect, and no action is defined between . Please confirm if this is intended.

Answer 1

[Davide Masiello]

You could try this

clear,clc

tf   = 3;               % Final time of integration
T0   = 20;              % Initial temperature

tend = 0;               % Initial time of integration
idx  = 0;               % Iteration index

while tend < tf
    
    if T0 < 23          % Sets positive heat input when T < 23
        Tx = 23;
        Q  = 1.5e+6;
    else                % Sets negative heat input when T >= 23
        Tx = 17;
        Q  = -1.5e+6;
    end
    
    idx      = idx+1;
    tspan    = [tend,tf];
    opts     = odeset('Event',@(t,T)heatSwitch(t,T,Tx));
    sol(idx) = ode15s(@(t,T) tempControl (t,T,Q), tspan, T0, opts);
    tend     = sol(idx).x(end);
    T0       = sol(idx).y(end);
    
end

plot(   [sol(:).x],[sol(:).y],'k',...
        [sol(:).x],ones(1,length([sol(:).x]))*23,'r--',...
        [sol(:).x],ones(1,length([sol(:).x]))*17,'r--')
xlabel('Time, s')
ylabel('Temperature, C')
axis([-inf +inf 16 24])

function dTdt = tempControl(t,T,Q)
dTdt =(Q/15172.5)+29.66+2.966*T;
end

function [position,isterminal,direction] = heatSwitch(t,T,Tx)
  position = T-Tx; % The value that we want to be zero
  isterminal = 1;  % Halt integration 
  direction = 0;   % The zero can be approached from either direction
end

However, please notice that when T = 23, I switch Q to -Q otherwise the temperature wouldn't go down back to 17. In fact, if Q = 0, dT/dt in your espression would still be positive.

[Sam Chak]

Hi @Davide Masiello

Thanks for your initiative to solve the Temperature Control Problem.

I think the thermostat (Q) is modeled incorrectly in @annop butmarn's original code, because it has only two switching logics, as shown in the figure posted in my comment above.

At the initial condition, is inside the range , what is the initial ? Is the temperature rising or dropping from ?

The hysteresis behavior of the thermostat (Q) is detailed in the "Title" above. Noticed you have added the HeatSwtich event. Agreed that Q cannot be positive at 23­°C in the original switch logic. Else, the sign in the ODE shall be negative, . By the way, can you plot Q vs. T? I'd like to see how to model the HeatSwtich event in mathematics.

[annop butmarn]

At the initial condition,T(0) = 20 degreeC. Temperature is dropping from 20 degreeC.

[annop butmarn]

@Davide Masiello and @Sam Chak thank you very much for helping me. I learn from you all answers. finally, I can solve it.

[Sam Chak]

Hi @annop butmarn

Kudos to you and thanks for updating us. @Davide Masiello did helped you, while me raising some query related to the hysteresis behavior of the thermostat.

If you don't mind, would you post your solution in a new answer as a reference for other people? Then you can "Accept" your own answer as a closure to this question.

[Davide Masiello]

Hi @Sam Chak - you have a good point there when you wonder what Q(0) we should consider, I wrongly assumed temperature would be incresing first, then hit 23 and go down, while @annop butmarn pointed out that it goes the other way round. This can be easily fixed in my code above, as I am sure @annop butmarn has already done by now.

Glad I could contribute to the discussion.

@Sam Chak Regarding plotting Q, did you mean Q vs time?

[Sam Chak]

Hi @Davide Masiello

Thanks for your reply and your contribution. Yes, I was referring the plot of Q vs, time. By the way, could you briefly advise me on how the heatSwitch(t,T,Tx) works in the ODE?

[Davide Masiello]

Sure. For a general overview of ODE Event Location functions, I recommend checking the following link

https://ch.mathworks.com/help/matlab/math/ode-event-location.html

In summary, the ode event locator checks if a condition is met at each time step of integration of a ODE.

You can then choose whether to stop the integration or not. In any case, you can store the values of t and y when the event happens.

In the particular case above, heatSwitch checks if the temperature reaches a target temperature Tx, which I have to be either 17 or 23 depending on what stage the process is at (quenching or heating).

[Sam Chak]

Hey @Davide Masiello

Thanks for info and the explanation. I'll look up the info.

Grazie mille.