why didn't the event function detect the events?

Hi all! I'm solving the thermostat model, which presents the characteristic of hybrid dynamical systems, when the ith (here i=1,2) room's temperature reduces to T_ref-0.5, the thermal is off, when the ith (here i=1,2) room's temperature raises to T_ref+0.5, the thermal is on, I use the ode suits with events, but it failed, plus the program also runs a long time, about 40 seconds(performance of computer is not bad), can anyone help me to fix the problem?

This is ode45 codes, but failed.

tic
 %options=odeset('Events',@Events1,'Events',@Events2,'AbsTol',1e-8,'RelTol',1e-8);
 options=odeset('AbsTol',1e-8,'RelTol',1e-8);
 y0 = [21;25;17];
 tspan=0:1000;
 [tout,yout]=ode45(@Tq_Tj,tspan, y0,options);
 figure(1)
plot(tout,yout(:,1),'k');hold on
plot(tout,yout(:,2),'-.');
plot(tout,yout(:,3),'--');
legend('room1','room2','wall')
toc

function f=Tq_Tj(~,y)
C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;Ta=32;p1=14;p2=14;T_ref=20;

    if(y(1)>T_ref+0.5)
        p1=14;
    elseif(y(1)<T_ref-0.5)
        p1=0;
    end
    if(y(2)>T_ref+0.5)
        p2=14;
    elseif(y(2)<T_ref-0.5)
        p2=0;
    end
      f=[((Ta-y(1))/R1+(y(3)-y(1))/R1w-p1)/C1; ...
          ((Ta-y(2))/R2+(y(3)-y(2))/R2w-p2)/C2;...
          ((y(1)-y(3))/R1w+(y(2)-y(3))/R2w)/Cw];
end
 
 function [g,isterminal,direction]=Events1(~,y)
 T_ref=20;
    g=[y(1)-(T_ref-0.5);y(1)-(T_ref+0.5)];
    isterminal=[0;0];
    direction=[0;0];
 end
  function [g,isterminal,direction]=Events2(~,y)
  T_ref=20;
    g=[y(2)-(T_ref-0.5);y(2)-(T_ref+0.5)];
    isterminal=[0;0];
    direction=[0;0];
  end

this is another way to solve the problem, you can run it to get the right results.

  clc
clear
tic
%为各变量赋初值
delta_q=1e-4;
q1=21;q2=25;q3=17;
x1=q1;x2=q2;x3=q3;

t=0;delta_t=0;
A=zeros(10,6);
n=0;
C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;Ta=32;p1=14;p2=14;m1=1;m2=1;T_ref=20;
%C1=0.56;C2=0.28;Cw=0.2;R1=5;R2=6;R1w=3;R2w=3;Ta=20;p1=4;p2=4;m1=1;m2=1;T_ref=20;
%开始while循环
while (t<1000)
    if(x1>T_ref+0.5)
        m1=1;
    elseif(x1<T_ref-0.5)
        m1=0;
    end
    if(x2>T_ref+0.5)
        m2=1;
    elseif(x2<T_ref-0.5)
        m2=0;
    end
    Dx1=((Ta-q1)/R1+(q3-q1)/R1w-m1*p1)/C1;  
    Dx2=((Ta-q2)/R2+(q3-q2)/R2w-m2*p2)/C2; 
    Dx3=((q1-q3)/R1w+(q2-q3)/R2w)/Cw;
    DDx1=(-Dx1/R1+(Dx3-Dx1)/R1w)/C1;
    DDx2=(-Dx2/R2+(Dx3-Dx2)/R2w)/C2;
    DDx3=((Dx1-Dx3)/R1w+(Dx2-Dx3)/R2w);
    DDDx1=(-DDx1/R1+(DDx3-DDx1)/R1w)/C1;
    DDDx2=(-DDx2/R2+(DDx3-DDx2)/R2w)/C2;
    DDDx3=((DDx1-DDx3)/R1w+(DDx2-DDx3)/R2w);
    
   %求Δt1和Δt2的值
       delta_t1=sqrt(2*delta_q/abs(DDx1));
       delta_t2=sqrt(2*delta_q/abs(DDx2));
       delta_t3=sqrt(2*delta_q/abs(DDx3));
     delta_tmin=min([delta_t1 delta_t2 delta_t3]);


   %比较Δt1和Δt2的大小，进而确定q1和q2谁先跃迁
   if (delta_t1==delta_tmin)
       delta_t=delta_t1;
       t=t+delta_t;
       caribe_1=-0.125*DDDx1*(exp(-2*delta_t)-1)+0.5*(DDx1+0.5*DDDx1)*delta_t^2+(Dx1-0.25*DDDx1)*delta_t;
       caribe_2=-0.125*DDDx2*(exp(-2*delta_t)-1)+0.5*(DDx2+0.5*DDDx2)*delta_t^2+(Dx2-0.25*DDDx2)*delta_t;
       caribe_3=-0.125*DDDx3*(exp(-2*delta_t)-1)+0.5*(DDx3+0.5*DDDx3)*delta_t^2+(Dx3-0.25*DDDx3)*delta_t;
       x1=x1+caribe_1;
       x2=x2+caribe_2;
       x3=x3+caribe_3;
       q1=x1;
       q2=q2+caribe_2;
       q3=q3+caribe_3;
   elseif (delta_t2==delta_tmin)
       delta_t=delta_t2;
       t=t+delta_t;
       caribe_1=-0.125*DDDx1*(exp(-2*delta_t)-1)+0.5*(DDx1+0.5*DDDx1)*delta_t^2+(Dx1-0.25*DDDx1)*delta_t;
       caribe_2=-0.125*DDDx2*(exp(-2*delta_t)-1)+0.5*(DDx2+0.5*DDDx2)*delta_t^2+(Dx2-0.25*DDDx2)*delta_t;
       caribe_3=-0.125*DDDx3*(exp(-2*delta_t)-1)+0.5*(DDx3+0.5*DDDx3)*delta_t^2+(Dx3-0.25*DDDx3)*delta_t;
       x1=x1+caribe_1;
       x2=x2+caribe_2;
       x3=x3+caribe_3;
       q2=x2;
       q1=q1+caribe_1;
       q3=q3+caribe_3;
   elseif (delta_t3==delta_tmin)
       delta_t=delta_t3;
       t=t+delta_t;
       caribe_1=-0.125*DDDx1*(exp(-2*delta_t)-1)+0.5*(DDx1+0.5*DDDx1)*delta_t^2+(Dx1-0.25*DDDx1)*delta_t;
       caribe_2=-0.125*DDDx2*(exp(-2*delta_t)-1)+0.5*(DDx2+0.5*DDDx2)*delta_t^2+(Dx2-0.25*DDDx2)*delta_t;
       caribe_3=-0.125*DDDx3*(exp(-2*delta_t)-1)+0.5*(DDx3+0.5*DDDx3)*delta_t^2+(Dx3-0.25*DDDx3)*delta_t;
       x1=x1+caribe_1;
       x2=x2+caribe_2;
       x3=x3+caribe_3;
       q3=x3;
       q1=q1+caribe_1;
       q2=q2+caribe_2;
   end
   n=n+1;
   A(n,1)=t;
   A(n,2)=x1;
   A(n,3)=x2;
   A(n,4)=x3;
   A(n,5)=m1;
   A(n,6)=m2;
end
% A;
 figure(2);
  plot(A(:,1),A(:,2),'--');hold on
  plot(A(:,1),A(:,3),'-.');
  plot(A(:,1),A(:,4),'-');
  grid on

%   plot(A(:,1),A(:,5)*22);hold on
%   plot(A(:,1),A(:,6)*22);
  legend('room1','room2','wall')
toc
       

Answer 1

[Mischa Kim]

In your ode45 call you are not returning any event information. Replace with

[tout,yout,te,ye,ie] = ode45(@Tq_Tj,tspan, y0,options);

and you should receive events information such as event times te, solution value at each of the event times, etc.

Regarding execution time, this also depends on the tolerances you set in the ode options. I do not know what your requirements are, however, try lowering the tolerances (or removing them entirely, just to see what happens) to see if the results are still satisfactory.

[汉武 沈]

Hello, I just did what you said, but the results is not satisfactory, the result is below.

The model is in fact that two air conditionings placed in two rooms, the wall is in the medium, I write it down in a paper, hope you can understand it.

And another thing What shoud I do with the 'te' and 'ie' information? Is the event function right, how to correct it?

Thank you very much!

Answer 2

[Jan]

Just my standard interjection: Matlab's ODE integrators are designed to integrate smooth functions. Your function Tq_Tj() conatins discontinuities. This can reduce the step size until the processing time gets huge and the result is dominated by rounding errors.

In

options=odeset('Events',@Events1,'Events',@Events2, ...    

you overwrite the Event function. Only the 2nd one is considered. Instead of trying to provide 2 functions, implement both events in 1 function. You have to stop the integration, when an event occurs to avoid the discontinuity. Restart the integration with the final values of the former interval.

[汉武 沈]

Hi，though I have understood your meaning, but sadly I don't know how to code, forgive my poor ability :), if you are free, can you help me get it done?

Thank you, anyway!